Re: [eigen] Feature idea: coeffRef() ---> setCoeff()

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

2010/7/29 Manoj Rajagopalan <rmanoj@xxxxxxxxx>:
> Hi eigen developers,
>
>   Grepping for coeffRef() in src/Core/*.h seems to show that all occurrences
> seem to be returning Scalar& for one of two purposes:
> (1) for use as lvalue in assignments,
> (2) in return expressions within member functions so that users may do the
> same at a higher level.

and (3) when one really wants the memory address where a given
coefficient is stored, which is useful when doing low-level
coefficient access or interfacing with other libraries.

So coeffRef() is certainly not going away; we could still consider
adding setCoeff() alongside if there was a good reason to.

>  A.selfadjointView<Lower>()(j, j+1, value); // note: write to upper triangle

oh... i see. I guess that you omitted setCoeff in the above line, by
the way, so you meant:

A.selfadjointView<Lower>().setCoeff(j, j+1, value); // note: write to
upper triangle

so this would write 2 coeffs at once in A ? Why not. Indeed that cant
be done with coeffRef. But I don't adhere to your rather complex plan
of generally allowing to write to Conjugate expressions, this is not a
prerequisite for that, we could just implement setCoeff in
selfadjointView, where it is useful, without implementing it in stuff
like conjugate() which IMO just don't want to be writable expressions.

This would be much simpler --- the patch could be just a few kB's.

But again, coeffRef() is definitely not going away.

Benoit

>
>  The nice '=' syntax is not present at the coefficient-wise level but,
> *  it is still present at matrix-expression level,
> *  but we gain a few more advantages (see below)
>
>   In some detail:
>
> 1. The Transpose<> class has been created to allow writes to the nested
> MatrixType even though the structure has been changed. This is because the
> Transpose class *logically* represents the structure transformation in a
> *reversible* way so as to allow coherent definitions of coeff() and
> coeffRef().
>
> 2. This idea can be extended to conjugates. A Conjugate<> wrapper, like
> Transpose<>, will reversibly encapsulate the coefficient-wise conjugate
> operation. IMPORTANT: it will enable writes because conjugation is a
> bijection on real/complex fields but this will depend on the existence of
> either of two alternatives:
>   (a)  a setCoeff() method as opposed to a "coeffRef() =" kind of assignment
> which reverse-transforms values before writing to memory, OR
>
>   (b) a wrapper class for Scalar& that is returned in coeffRef (instead of
> Scalar&), whose operator=() is overloaded to perform the
> reverse-value-transformation.
>
>
> 3. The benefit here is that we may define a HermitianMatrix or
> SelfAdjointMatrix class for which writes to either triangle will now be
> coherent from the user's view: Eigen can internally map the write down to the
> stored/viewed triangular portion (I am working on CompactSelfAdjointMatrix
> and related classes that do this - will release soon).
>
>   (a) We may even define SkewHermitianMatrix where the elements must be
> conjugated and negated between the two triangles. Not possible easily with
> existing architecture.
>
>   (b) Such specialized types capture important physical information about the
> problem at reduced memory requirements and improve bug-safety.
>
>
> 4. Generally, in any view of an underlying dense matrix type, coherent writes
> to the coefficients can be performed if the view transforms each coefficient
> in a 1-1 lossless manner. Structural rearrangements like transpose and
> permutations allow direct access of Scalar& because only a computation on the
> indices is involved. But value-transformations like conjugate() will require
> a setCoeff() method because the value is transformed (reversibly).
>
>
> 5. By allowing both structural and value reversible transformations
> simultaneously, we can extend Eigen by composing these features. For example,
> the return type of adjoint() can now be Conjugate<Transpose<MatrixType> >
> which can now *allow writes* through the setCoeff() call-chain. I imagine
> this will be useful in some inplace-solve functions. Currently,
> AdjointReturnType is built similarly, but because conjugate() is a "const"
> method, no writes are possible.
>
>
> 6. FAR-FETCHED OPTIONAL IDEA: The above idea can be extended to other
> reversible value-transformations if needed (by the user). For example, the
> following contrived example becomes possible for free:
>
>   VectorXd v(N);
>   (5.0*v).setConstant(1.0); // v[i]=0.2 for all i
>
>  This form of expression might be natural in some contexts where scales are
> important (eg. performing energy computations in eV instead of J).
>
>   (a) why the user would do the above is upto the user but the functionality
> is available for free because ei_scalar_product_op is 1-1 for non-zero scales
>
>   (b) For many ops like scalar_sum_op and scalar_difference_op, the "inverse"
> op is easy to derive and in fact can be templated in the traits class. Eg:
>
>    template<typename Scalar>
>       struct ei_traits<ei_scalar_add_op> {
>           //...
>           typedef ei_scalar_subtract_op inverse_op; // user has to
> instantiate with negative of add-op's parameter though
>      };
>
>   If for each op, an inverse op can be so defined then the architecture will
> naturally allow synthesis of write-able cwise-transformed views of dense
> matrix storage types.
>
>   (c) it would be up to the user to avoid non-sensical parameters like 0 for
> the product which would render the inverse undefined
>
>   (d) (5.0*v) is a lossy transformation: bits are lost in finite-precision
> arithmetic - the losslessness is only in the algebraic definition of the *
> operator but as long as the user is OK with this, he/she can use these
> features.
>
>
> 7.  I've included the above far-fetched idea as an example of how users can
> extend features beyond what the developers intended. Template Metaprogramming
> itself is a result of such serendipity. The above example might be useless
> but who knows someone might figure out something really useful :-)
>
> Thanks,
> Manoj
>
>
>