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

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

   If this is indeed the case, then I think replacing this with the following 
will allow greater extensibility:     setCoeff(row,col,value)

   in short, we can have expressions like:
 
  MatrixXcd A(N,N);
  A.conjugate().setCoeff(i, j, value);
  A.selfadjointView<Lower>()(j, j+1, value); // note: write to upper triangle

  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