Hi!

I am new to C++ and have trouble specializing typename template parameter for Eigen::Matrix types. I have to do this because of the mandatory use of .cwise() for component-wise operations.

Consider the following silly sample code:

[code]const int SIZE = 5;

template
class Array {
public:
// ... Ctor ...

const T& operator[](int ind) const { return data[ind]; }

void someFunc(const Array& array);

private:
T data[SIZE];
};

template
void Array::someFunc(const Array& array) {
for (int i = 0; i -s or any Eigen::Matrix type (such as Vector3cd). What I want is, that *= in someFunc() is ordinary multiplication for T = double, complex, etc. and .cwise() component-wise multiplication for T = Eigen types.

Do I have to specialize the template class Array to do this, and add .cwise() to the definition? If yes, please help, how do I do that?
Or do you have any advanced idea how to solve this problem?

Thank you very much!

You'd just declare the function with an empty template parameter and replace it with the type in the function itself. E.g.

template
void Array::someFunc(const MatrixBase& array) {
return data.cwise() * array;
}

But let me tell you that I find this to be quite strange. Do you really need such a wrapper? It likely makes the resulting code quite hard to read for other people.

It is just an example of the problem, of course, the real code doesn't look like this (it is a physical optics simultation of laser beams).

What I would like to do, is to represent an ordered set of mathematical objects (real/complex scalar field values, real/complex vector field components).
Fundamentally, they are stored in a std::vector container.
For e.g. complex vector field components
std::vector
with the Eigen/StdVector header included, as advised by the documentation.

But std::vector is pretty dumb in the sense that it is only for storage, doesn't provide any mathematical operations...

That's why I created a
template class Array
class, which stores the component values in a std::vector container, and also includes an operator *= for the element-wise multiplication with another Array (and also other related member functions to export and import such an array to/from a file, etc.).
*= is desired to do the multiplication on the corresponding elements component-wise. And that needs different syntax when e.g. T = double ("plain *=") vs. T = Eigen::Vector3cd (".cwise() *=").

Can you give a better solution to the whole problem?

Btw, I haven't tried the code yet you have suggested but doesn't MatrixBase need any template parameters? Or am I wrong?

Andre wrote:You'd just declare the function with an empty template parameter and replace it with the type in the function itself. E.g.

template
void Array::someFunc(const MatrixBase& array) {
return data.cwise() * array;
}

But let me tell you that I find this to be quite strange. Do you really need such a wrapper? It likely makes the resulting code quite hard to read for other people.

If you are not into template programming and want to write "readable" code, what about:

Code: Select all

   template
   class MyArray
   {
   public:
      typedef T* TArrayType;
      typedef unsigned int TIndexType;
   protected:
      TArrayType  _data;
      TIndexType _size;
   public:
      MyArray();
      MyArray( MyArray & other);
      MyArray & operator=(const MyArray &other)
      {
         _size = other.size();
         for (TIndexType dim(0);dim!=other.size();++dim) _data[dim]=other._data[dim];
         return *this;
      }
      const T & operator[](const TIndexType & index) const
      {
         return _data[index];
      }
      T & operator[](const TIndexType &  index)
      {
         return _data[index];
      }
      MyArray & operator+=(const MyArray & other)
      {
         return *this = *this + other;
      }
      const MyArray operator+(const MyArray & other) const
      {
         here: create a temporary MyArray object, set its data and return it
         return tmp;
      }
      void swap(MyArray& other); //  may be useful as well?
      const TArrayType & data() const {return _data;}
      TArrayType & data() {return _data;}
   };


First of all, thanks a lot to Andre and Seb, it's great to see people helping each other, that's the whole point of this forum. Here are a couple more ideas (didn't reply earlier as I'm very busy so it's really awesome that others do!)

If all your problem is that Eigen objects require the cwise() prefix for some operators, then the easiest solution is probably that you write a MatrixBase plugin allowing normal operators to work without cwise(). This is explained here:
http://eigen.tuxfamily.org/api/Customiz ... MatrixBase

This has also been discussed in some threads on this forum. We don't do that by default in Eigen but in your case since you want this semantic anyway in your Array class, why not extend MatrixBase.

Another solution is that you do for your Array class the same as we did in Eigen with cwise(). So you have a Cwise wrapper, and implement the cwise operations there. Of course you may not want to do that if you find this heavy.

Another solution is to write template wrappers around these operators and then specialize them for Eigen types. What you probably don't want to do, is to specialize your whole Array class, as you would have to rewrite everything for each specialization! You can't specialize just 1 member function, and you can't partially specialize functions, so all advanced specialization has to go through helper structs.

Here's an example.

Code: Select all

// the generic case
template
struct operator_plus_wrapper
{
  typedef T ReturnType;
  ReturnType run(const T& t, const U& u)
  {
    return t+u;
  }
};

// Partial specialization : adding a number to a vector/matrix expression.
// Since vector/matrix expressions are exactly the subclasses of MatrixBase
// for some type T, we use that to "filter" and do the partial specialization on.
// Notice that in MatrixBase, T is the "derived type" so we can use it directly --
// the only reason why we mention MatrixBase at all is in order to do the partial
// specialization.
template
struct operator_plus_wrapper, U>
{
  typedef T ReturnType;
  ReturnType run(const T& t, const U& u)
  {
    return t.cwise()+u;
  }
};

// now we can use our wrapper....

template
class Array
{
  public:
  template
  void operator+(const U& u)
  {
    for(int i = 0; i ::run(m_array[i], u);
  }

  protected:
  int size;
  T *m_array;
};


bjacob wrote:If all your problem is that Eigen objects require the cwise() prefix for some operators, then the easiest solution is probably that you write a MatrixBase plugin allowing normal operators to work without cwise(). This is explained here:
http://eigen.tuxfamily.org/api/Customiz ... MatrixBase

unfortunately this won't work for operator* or operator*= since they already correspond to the matrix-matrix product.

Another solution is to write template wrappers around these operators and then specialize them for Eigen types. What you probably don't want to do, is to specialize your whole Array class, as you would have to rewrite everything for each specialization! You can't specialize just 1 member function, and you can't partially specialize functions, so all advanced specialization has to go through helper structs.

Here's an example.

Code: Select all

// the generic case
template
struct operator_plus_wrapper
{
  typedef T ReturnType;
  ReturnType run(const T& t, const U& u)
  {
    return t+u;
  }
};

// Partial specialization : adding a number to a vector/matrix expression.
// Since vector/matrix expressions are exactly the subclasses of MatrixBase
// for some type T, we use that to "filter" and do the partial specialization on.
// Notice that in MatrixBase, T is the "derived type" so we can use it directly --
// the only reason why we mention MatrixBase at all is in order to do the partial
// specialization.
template
struct operator_plus_wrapper, U>
{
  typedef T ReturnType;
  ReturnType run(const T& t, const U& u)
  {
    return t.cwise()+u;
  }
};

// now we can use our wrapper....

template
class Array
{
  public:
  template
  void operator+(const U& u)
  {
    for(int i = 0; i ::run(m_array[i], u);
  }

  protected:
  int size;
  T *m_array;
};


again this is not going to work because your objects are not of the type MatrixBase but of the type Derived and so your specialization of op_plus_wrapper will never been called. However, if you know that T cannot be a generic eigen expression but always a Matrix type (that seems to be the case), then you can do it with Matrix instead of MatrixBase.

Another solution: if you can restrict yourself to scalar or vector types (no matrix), then you can directly use a Matrix instead of your Array class. So if you store your data in column major, instead of writing:

array.at(i), array.size(), etc...

you would write:

array.col(i), array.cols(), etc...

Then, you might also add some convenient functions to Matrix, for instance "at(int)" as an alias for "col(int)", push_back, etc...


Actually, this remind me that a while ago (before the .cwise() solution) I hacked Eigen so that you can have an "array mode" for which operator* would mean cwise product and allowing you to use a Matrix as the scalar type of a Matrix: Matrix, ... >, and this with complete lazy evaluation. However this was really really complicated, heavy for the compiler, and a bit dangerous from the user point of view.

I forgot to say that the solution I proposed (using Matrix instead of Array) give you full vectorization with lazy evaluation.

Also might write your own Array class inheriting MatrixBase >, just like our Matrix but with a nicer API matching your specific use case. Actually, I'm thinking such a class would make a very useful module.

Thank you very much for your ideas. I have constructed a similar solution to the one of bjacob's, using a wrapper class and specializing only for Eigen matrices to use the cwise-formalism in that case.

Question: does it mean any overhead to use these kind of wrappers, when all function calls are inlined in the wrapper definition?

ggael:
As I also need "matrix fields" (to do matrix-vector multiplications point-wise --> rotations), for now it's simpler to represent these fields as an 1d array (std::vector at the moment) of the corresponding math type (scalar, vector, matrix). How efficient it is, I don't know. I compile the code with g++ 4.3.x, with vectorization enabled, hoping that it will vectorize the loops running over all field points, too, if beneficial.

But your idea of representing the field as a Matrix of Matrices is valuable, and perhaps worth experimenting with (unfortunately I don't have much time for this now, 'cause I have to finish my MSc thesis in time)... Blitz++ also proposes an Array ... > construct as a vector-field. In that case it is simpler because all arithmetic operators mean coeff-wise operations by default (I have picked Eigen over Blitz++, because it is lightweight and actively maintained).

If you have more thoughts to share, please post them