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

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Hi,

2010/7/7 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.

OK, so, I see where you are coming from. You can't use unaryExpr since
this is not just a coeff-wise operation on M (and the result even has
a different size).

So I see, NullaryExpr is all you can do if you don't want to write
your own expression class. It's a good way of prototyping this. But I
think that you should really implement your own new expression class
for that, rather than using NullaryExpr.

So a good starting point would be the Minor class from Eigen2, since
it's a very simple expression class that does a special case of what
you are doing here. In Eigen3, you can find it in the Eigen2Support
module (it has been removed from the standard Eigen3 API and requires
one defining EIGEN2_SUPPORT):

Eigen/src/Eigen2Support/Minor.h

What it does is provide an expression of a matrix with one row and one
column taken out.

>
> I think that something like ei_functor_has_linear_access<NullaryOp>
> but for isRepeatable is needed, or am I missing something?

See? By not using NullaryExpr, you avoid having to deal with all that ;-)

> I understand that this approach works fine only if one wants to
> read-only the matrix, am I correct?

Indeed, NullaryExpr doesn't allow usage as a lvalue, since it was
never meant to have any underlying matrix storing coefficients! The
fact that you pass the matrix M here as argument is a hack ;-)

Again the fix for this is to write your own expression class, like Minor.

> If one wants to modify the matrix, one wants something close to a map,
> is there a constructor of a map taking a functor ?

Map is just mapping a plain matrix or array. The closest thing to what
you are asking is CwiseUnaryView. It is the lvalue variant of
CwiseUnaryOp. However, you can't use it in your case because it is
only for coefficient-wise operations... whence the 'cwise' in its name
:)

> 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.

The expression class in question, allowing to take a sub-matrix of a
matrix by specifying a range or a set of rows/cols, would be very
welcome in Eigen 3.1. If you want to try doing it, you are very
welcome. Yes, it would be slow because non-vectorizable. The only
vectorizable case would be the case of contiguous rows/cols but then
we have Block already for that. So it's OK, in my opinion, to be
non-vectorizable here.

The API should look like

   VectorXi rowIndices(3);
   rowIndices << 4, 8, 10;
   VectorXi colIndices(4);
   rowIndices << 0, 2, 5;

   MatrixXd matrix(100, 100)

   matrix.goodFunctionName(rowIndices, colIndices)  // returns a 3x4
matrix expression that can be used as rvalue or as lvalue

Benoit

>
> 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 !
>
>
>