Re: [eigen] FFT update

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

Hi,

Here's an example: see attached file.
Enjoy!!
Benoit

2010/1/22 Benoit Jacob <jacob.benoit.1@xxxxxxxxx>:
> ok let me write an example...
>
> 2010/1/22 Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>:
>> On Fri, 22 Jan 2010, Benoit Jacob wrote:
>>
>>> Our "policy" is now:
>>> - try to return by value, optimizing with the ReturnByValue class,
>>> wherever possible;
>>> - otherwise you can just use references.
>>
>> Can somebody please explain how to use the ReturnByValue class?
>> Or should I just extrapolate from the source code?
>> I could not find any documentation.
>>
>> Jitse
>>
>>
>>
>

/*** Bon appetit bien sur! With chef Benoit.
 ***
 *** Today we'll learn to cook a closure-like function with optimized returning-by-value
 *** thanks to the ReturnByValue class.
 ***/

#include <Eigen/Core>

using namespace Eigen;
using namespace std;

// here is a function naively returning a matrix or array by value.

template<typename Derived>
typename Derived::PlainMatrixType // get the plain matrix/array type corresponding to this expression type.
                                  // yes, this typedef is from the days when we didn't have Array.
rotate_naive(const DenseBase<Derived>& src)
{
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
  typename Derived::PlainMatrixType result(src.size());
  const int n = src.size()-1;
  result.tail(n) = src.head(n);
  result.coeffRef(0) = src.coeff(n);
  return result; 
}

// now the problem is that when we do:
//   y = rotate_naive(x)
// the result is first computed in a temporary variable that is then copied into y.
// In order to optimize that, we use the ReturnByValue class. It is a return-value-proxy, or a
// "closure" if you want. The idea is that the rotate() function should not return the whole result
// array/matrix, instead it should return a lightweight ReturnByValue object that only describes
// the operation, and the operation is only done by the subsequent "=". So this is lazy evaluation.
// The reason why this is nontrivial to implement is that the user might instead be doing something like
//   rotate(x) + something_else
// So the ReturnByValue class is a real expression in the sense of Eigen's expression templates, and
// upon nesting into a larger expression, it immediately evaluates its result.

// The ReturnByValue class is not used directly, instead you derive it to specialize it to your problem.
// So let's define a derived class rotate_retval.

// First we declare it.
template<typename Derived> struct rotate_retval;

// Then we specialize ei_traits for it. After all, it's an expression class like any other.
// Don't forget that ei_traits belongs to namespace Eigen!
namespace Eigen {
template<typename Derived>
struct ei_traits<rotate_retval<Derived> >
{
  typedef typename Derived::PlainMatrixType ReturnMatrixType;
    // get the plain matrix/array type corresponding to this expression type.
    // yes, this typedef is from the days when we didn't have Array. Need to rename it.
};
}

// Now we can define the rotate_retval class. It inherits ReturnByValue<rotate_retval>:
// yes, this is again the famous "curiously recursive template pattern".
// The ReturnByValue base class is an expression class in the sense of expression templates:
// it inherits DenseBase. As such it brings us all the expression template facilities that we need,
// so we only have minimal stuff to implement here.
template<typename Derived> struct rotate_retval
 : public ReturnByValue<rotate_retval<Derived> >
{
  // The constructor is all what's called by the rotate() function that we implement below.
  // All it does is store a reference to the source matrix/array.
  rotate_retval(const Derived& src) : m_src(src) {}

  // almost all we have to do is explain how the "rotate" expression will be evaluated. This is where
  // the implementation goes! This time, the return value is a reference, Dest& dst, so we don't
  // actually return the result by value.
  template<typename Dest> inline void evalTo(Dest& dst) const
  {
    dst.resize(m_src.size());
    const int n = m_src.size()-1;
    dst.tail(n) = m_src.head(n);
    dst.coeffRef(0) = m_src.coeff(n);
  }

  // oh right and we must also say how many rows/columns our expression has. This is needed to be
  // able to create the destination matrix/array when we need to evaluate into a temporary.
  int rows() const { return m_src.rows(); }
  int cols() const { return m_src.cols(); }

  protected:
    const Derived& m_src;
};

// Finally we can implement our rotate() function.
template<typename Derived>
rotate_retval<Derived> rotate(const DenseBase<Derived>& src)
{
  // this is the right place for assertions, so errors are caught are the place where they exist.
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived)
  return rotate_retval<Derived>(src.derived());
     // do you know how I lost 30 minutes debugging this shit? I forgot the derived() in the above line.
     // So my src object of type DenseBase<Derived> was being deep-copied into a temporary of type Derived.
     // The temporary then died, and I was reading undefined values. So, don't forget the derived().
     // Thanks for listening.
}

/*** Bake in GCC for 40 minutes... ***/

// Now let's taste it:
int main()
{
  RowVector3f v(1,2,3);
  cout << "Here's v: " << v << endl;

  // Test rotate_naive()
  cout << "Here's rotate_naive(v): " << rotate_naive(v) << endl;
  
  // Test rotate()
  cout << "Here's rotate(v): " << rotate(v) << endl;
  // Wait! You ask. Since rotate(v) returns an expression of type rotate_retval that can only be evaluated
  // and on which we may not call coeff(), how was it possible to print it? Didn't it try to call coeff()
  // on it to print the coeffs? No! Because the first thing that cout<< does, is to _nest_ the expression
  // that is passed to it. ReturnByValue says "when nesting me, you must first evaluate me". So the expression
  // gets evaluated into a temporary before being printed.

  // But this would crash (infinite recursion):
  //   rotate(v).coeff(0)
  // or
  //   rotate(v)(0).

  // Just to show off, let's nest rotate() into a bigger expression:
  cout << "Here's 2*rotate(v): " << 2*rotate(v) << endl;
  // see? no problemo.
}

// Bon appetit bien sur!!!