Re: [eigen] how to create a "virtual" array or matrix?

[Manuel Yguel <manuel.yguel@xxxxxxxxx>]

Hello everyone,

I would like to continue this thread, starting from the remark of Benoît.

Case study:
    a friend of mine want to build the following matrix (using Matlab
syntax) with Eigen
                    M([1,3,2,5,8],1:3:9)
    given M is a 9x9 matrix
note the tricky fact that rows 3 and 2 are inverted.

I used the pointer of Benoît (that I repeat for consistency):
http://eigen.tuxfamily.org/dox-devel/classEigen_1_1DenseBase.html#a84c2897928216cfb21ab81917e03b3a0
Stop reading further if you think this is totally not the good approach.

I hit the fact that isRepeatable is not defined in the default
implementation of ei_functor_traits in Core/util/XprHelper.h

template<typename T> struct ei_functor_traits
{
  enum
  {
    Cost = 10,
    PacketAccess = false
  };
};

I think that something like ei_functor_has_linear_access<NullaryOp>
but for isRepeatable is needed, or am I missing something?
So I just specialized ei_functor_traits for my use case.
It works fine, at the end of this mail is a code compiling and executing fine.

I understand that this approach works fine only if one wants to
read-only the matrix, am I correct?
If one wants to modify the matrix, one wants something close to a map,
is there a constructor of a map taking a functor ?
I guess not, since map assume linear or stride base access.

I am willing to add this functionality (for any kind of functor, then
especially for the matlab syntax) to push up my Eigen skills.
I am also willing to document the approach proposed by Benoît since it
needs a bit of Eigen skills (not too much I agree) and because I think
it might be interesting for a lot of people.

Note that I am totally aware that this kind of approach is very likely
to kill performances and that a correct ordering of the matrix or a
copy might be the best solution.
This problem came to me as a friend request who needed to avoid copy.
However I think that for prototyping, this feature might be very
usefull, especially to matlab users.

Thanks,

Manuel




#include <Eigen/Core>
#include <iostream>

using namespace Eigen;
using namespace std;


template< class SCALAR, class ORIGINAL_MATRIX >
struct CustomNicolasAccessor
{
  static const int map[];

  inline
  const SCALAR operator()( int& row, int& col ) const
  {
    int coeff = col*5+col;
    return mat[ map[coeff] ];
  }

  inline
  SCALAR& operator()( int& row, int& col )
  {
    int coeff = col*5+row;
    return mat[ map[coeff] ];
  }

  inline
  const SCALAR operator()( int& idx ) const
  {
    return mat[ map[idx] ];
  }

  inline
  SCALAR& operator()( int& idx )
  {
    return mat[ map[idx] ];
  }



  ORIGINAL_MATRIX& mat;

  CustomNicolasAccessor( ORIGINAL_MATRIX& m )
  :mat( m )
  {}
};

template< class SCALAR, class ORIGINAL_MATRIX >
const int CustomNicolasAccessor<SCALAR,ORIGINAL_MATRIX>::map[] = { 0,
2, 1, 4, 7, 18, 20, 19, 22, 25, 45, 47, 46, 49, 52, 72, 74, 73, 76, 79
};

namespace Eigen {
  template<typename SCALAR, class ORIGINAL_MATRIX>
    struct ei_functor_traits< CustomNicolasAccessor<SCALAR,ORIGINAL_MATRIX> >
    {
      enum {
        Cost = 10,
        PacketAccess = false,
        IsRepeatable = false
      };
    };
}


int main()
{

  MatrixXi orig(9,9);
  int k=0;
  for( int i=0; i<9; ++i )
  {
    for( int j=0; j<9; ++j ){
      orig(i,j) = k++; }
  }
  cout << orig << endl;
  cout << endl;

  CustomNicolasAccessor<int, MatrixXi> acc( orig );
  cout << MatrixXi::NullaryExpr( 5, 4, acc );
  cout << endl;
}








P.S. I came back from holidays yesterday (3 weeks ;D) and I am just
amazed by how much work you did for the beta during that time. Cheers
folks !