Re: [eigen] Eigen Types as Parameters for Functions

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi all,

I've implemented a solution fully described there:
http://eigen.tuxfamily.org/bz/show_bug.cgi?id=481

feedback more than welcome!

For the lazy ones, here is a self-explanatory example of the proposed solution:

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;

// read-write access to 'a'
void foo1(Ref<VectorXf> a) { }

// read-only access to 'a'
void foo2(Ref<const VectorXf> a) { }

int main()
{
  VectorXf a(10);
  const VectorXf& ac(a);
  VectorBlock<VectorXf> ab(a,0,3);
  MatrixXf A(10,10);
  const VectorBlock<VectorXf> abc(a,0,3);

  foo1(a);
  //foo1(ac);       // does not compile because ac is const
  foo1(ab);
  foo1(a.head(4));
  foo1(abc);
  foo1(A.col(3));
  foo1(A.row(3));   // copied into a temp because innerstride!=1

  foo2(A*A.col(1)); // evaluated into a temp
  foo2(ac.head(5));
  foo2(ac);
  foo2(a);
  foo2(ab);
  foo2(a.head(4));
  foo2(a+a);        // evaluated into a temp

  return 0;
}

cheers,
Gaël

On Mon, Feb 27, 2012 at 7:44 PM, Gael Guennebaud
<gael.guennebaud@xxxxxxxxx> wrote:
> btw,
>
> this more or less what Christian suggested but more restrictions because:
> 1 - the storage order really has to be known at compile time otherwise
> we cannot construct an expression
> 2 - same for the alignment, especially for small fixed sizes.
> Otherwise the alignment is already determined at runtime by checking
> the pointers, so no need to specify it.
> 3 - a runtime inner-stride can really kill the performance so it is
> important to let the user decide whether he can tolerate that or not.
> 4 - distinguishing at compile-time between vectors and matrices is
> fundamental too
>
> So if we go for a specific class (not Map) it could had 5 template parameters:
> - typename scalar
> - int rows
> - int cols
> - int storage_order
> - int inner_stride = 1
>
> cheers
> gael
>
> On Mon, Feb 27, 2012 at 7:34 PM, Gael Guennebaud
> <gael.guennebaud@xxxxxxxxx> wrote:
>> Alright, I thought I had already replied to this thread but I haven't
>> so here I go.
>>
>> My general idea to offer the ability to write non template function
>> while still offering some flexibility on the input is really simple:
>> simply add constructors to Map<> that can take Eigen's objects such as
>> Matrix, Block, Map, etc. Then it is up to user to specify the level of
>> flexibility he wants for his arguments: dynamic/fixed sizes,
>> outer-stride? inner-stride? matrix/vector?
>>
>> For instance:
>>
>> void foo(const Map<MatrixXd,0,OuterStride<> >& arg);
>>
>> should accept any MatrixXd, Matrix4d, Block, Map, etc.
>>
>> There are some limitations though:
>> - No template is possible at all, in particular it is not possible to
>> templatize the scalar type, but this is fine to me since if the user
>> can tolerate some template parameter then why not allowing everything.
>> - The storage order has to be fixed too. I understand a runtime
>> storage order might be convenient, but this is really not the Eigen
>> philosophy, so this limitation will have to stay.
>>
>> However, I have question regarding the handling of non compatible
>> expressions. For instance what should we do with:
>>
>> foo(A+B);  // this is a pure expression without direct access
>>
>> or
>>
>> MatrixXcf C;
>> foo(C.real());  // C.real() has an inner-stride, but foo does not accept that
>>
>> If we stick with Map<> I think that there is only one possible
>> solution: compilation error. Then it is up to the user to explicitly
>> evaluate it.
>>
>> Another approach would be to add a variant of Map tailored for this
>> use case that would allow to silently evaluate the argument into a
>> temporary. This is what happen if the user declare his argument as a
>> const MatrixXd&.
>>
>> Any opinion on which approach we should follow?
>>
>> cheers,
>> gael
>>
>> On Thu, Jan 26, 2012 at 3:12 PM, Christoph Hertzberg
>> <chtz@xxxxxxxxxxxxxxxxxxxxxxxx> wrote:
>>> On 26.01.2012 14:50, Benoit Jacob wrote:
>>>>
>>>> 2012/1/26 Christoph Hertzberg<chtz@xxxxxxxxxxxxxxxxxxxxxxxx>:
>>>>>
>>>>> On 26.01.2012 14:33, Benoit Jacob wrote:
>>>>>>
>>>>>>
>>>>>> 2012/1/26 Christian Seiler<christian@xxxxxxxx>:
>>>>>>
>>>>>>>  - subclass of MatrixBase, i.e. can be used inside the function like a
>>>>>>>   Matrix object
>>>>>>>  - has private members containing the above-mentioned bits of
>>>>>>>   information (pointer to first element, dimension, ...)
>>>>>>>  - can be implicitly converetd to from any compatible MatrixBase (i.e.
>>>>>>>   foo(mat) and foo(mat.block(...)) can be used directly)
>>>>>>
>>>>>>
>>>>>>
>>>>>> Just note that this is only possible for those matrix expressions
>>>>>> whose coefficients are stored in memory with a simple enough layout
>>>>>> that can be described in terms of strides. That is what we call a
>>>>>> "Direct Access" matrix expression; the test is ExpressionType::Flags&
>>>>>> DirectAccessBit.
>>>>>
>>>>>
>>>>>
>>>>> Yes, and I would say this would be the main application, if you have
>>>>> complicated-enough functions (i.e. where you don't really bother about
>>>>> evaluating expressions to a temporary if required).
>>>>> And more importantly, when having parameters that are actually output
>>>>> parameters, you are more or less stuck to direct access types anyway --
>>>>> and
>>>>> this is a very important point, where currently const-correctness is
>>>>> ignored
>>>>> (almost) entirely.
>>>>>
>>>>>
>>>>>> With that said, yes, we have been considering doing this, that's bug 58:
>>>>>> http://eigen.tuxfamily.org/bz/show_bug.cgi?id=58
>>>>>>
>>>>>> If someone wants to do this and needs mentoring, I should find time to
>>>>>> do at least that mentoring.
>>>>>
>>>>>
>>>>>
>>>>> I'd be interested at least to join the "specification committee" ;)
>>>>
>>>>
>>>> Really, I haven't thought about this in a long time so your thoughts
>>>> on this subject are probably a lot fresher than mine. From the top of
>>>> my head, the starting point is to look at MapBase (which is inherited
>>>> both by Map and by Block<direct access type>) and ask: how much of the
>>>> templated stuff there can we replace by runtime variables (data
>>>> members), and try to make all direct-access types inherit the
>>>> resulting base class (in particular, Matrix should inherit it).
>>>
>>>
>>> Yes, my thoughts were about the same. Actually, the only difference between
>>> a Map/Block and a Matrix should be that the Matrix "owns" its memory, i..e.
>>> has to deallocate it on destruction (with a special case of fixed size
>>> matrices having static memory).
>>>
>>> Furthermore, if one allows strides in the Matrix-class, some specialized
>>> types, such as AlignedVector3 would become just a typedef. Currently, if
>>> someone needs an aligned vector7, he'd have to to the same thing again --
>>> also it might be interesting to store an aligned 3x3 matrix as the top 3
>>> rows of a 4x3 matrix allowing a much more efficient matrix-vector product.
>>>
>>>
>>> Christoph
>>>
>>>
>>>
>>> --
>>> ----------------------------------------------
>>> Dipl.-Inf. Christoph Hertzberg
>>> Cartesium 0.051
>>> Universität Bremen
>>> Enrique-Schmidt-Straße 5
>>> 28359 Bremen
>>>
>>> Tel: (+49) 421-218-64252
>>> ----------------------------------------------
>>>
>>>