Constructing Eigen expression templates with Boost.Proto

Question

I'd like to use Boost.Proto to transform an embedded domain-specific language into a series of matrix operations implemented with the Eigen library. Since efficiency is important, I want proto to generate Eigen expression templates and avoid premature evaluation.

I've implemented a simple grammar that can generate matrix multiplication expressions. The code below compiles without warnings (on g++ 4.8.0 and Intel C++ 2013.3, with Boost 1.54.0 and Eigen 3.1.3) and works as long as my expression only has a single multiplication operation. As soon as I add more multiplications to the chain, it crashes. Valgrind tells me that this is because one of the Eigen::GeneralProduct expression template temporaries gets destroyed before the evaluation is completed.

I don't understand why this happens, or what I can do to prevent it. All help is appreciated!

#include <iostream>

#include <boost/fusion/container.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/void.hpp>
#include <boost/proto/proto.hpp>
#include <boost/ref.hpp>
#include <boost/type_traits/remove_const.hpp>
#include <boost/type_traits/remove_reference.hpp>
#include <boost/utility.hpp>

#include <Eigen/Dense>

namespace fusion = boost::fusion;
namespace mpl = boost::mpl;
namespace proto = boost::proto;

typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> matrix;

// Placeholders

const proto::terminal<mpl::int_<0> >::type I1 = {{}};
const proto::terminal<mpl::int_<1> >::type I2 = {{}};
const proto::terminal<mpl::int_<2> >::type I3 = {{}};

// Grammar

template<class Rule, class Callable = proto::callable>
struct External :
    proto::when<Rule, proto::external_transform> {};

struct matmul_transform : proto::callable {
    template<class Sig> struct result;

    template<class This, class MatrixExpr1, class MatrixExpr2>
    struct result<This(MatrixExpr1, MatrixExpr2)> {
            typedef typename Eigen::ProductReturnType<
                    typename boost::remove_const<typename boost::remove_reference<MatrixExpr1>::type>::type,
                    typename boost::remove_const<typename boost::remove_reference<MatrixExpr2>::type>::type>::Type
                    type;
    };

    template<class MatrixExpr1, class MatrixExpr2>
    typename result<matmul_transform(MatrixExpr1, MatrixExpr2)>::type
    operator()(const MatrixExpr1 &a, const MatrixExpr2 &b) const {
            return a * b;
    }
};


struct MatmulGrammar;

struct InputPlaceholder : proto::terminal<proto::_> {};

struct MatrixMultiplication :
    proto::multiplies<MatmulGrammar, MatmulGrammar> {};

struct MatmulGrammar : proto::or_<
    External<InputPlaceholder>,
    External<MatrixMultiplication> > {};

struct matmul_transforms : proto::external_transforms<
    proto::when<MatrixMultiplication, matmul_transform(MatmulGrammar(proto::_left), MatmulGrammar(proto::_right))>,
    proto::when<InputPlaceholder, proto::functional::at(proto::_data, proto::_value)> > {};

int main() {
    matrix mat1(2,2), mat2(2,2), mat3(2,2), result(2,2);

    mat1 << 1, 2, 3, 4;
    mat2 << 5, 6, 7, 8;
    mat3 << 1, 3, 6, 9;

    MatmulGrammar mmg;

    // THIS WORKS:
    result = mmg(I1 * I2,
            mpl::void_(),
            (proto::data = fusion::make_vector(boost::cref(mat1), boost::cref(mat2), boost::cref(mat3)),
             proto::transforms = matmul_transforms()));

    std::cout << result << std::endl;

    // THIS CRASHES:
    result = mmg(I1 * I2 * I3,
            mpl::void_(),
            (proto::data = fusion::make_vector(boost::cref(mat1), boost::cref(mat2), boost::cref(mat3)),
             proto::transforms = matmul_transforms()));

    std::cout << result << std::endl;

    return 0;
}

Remember

There was a talk in this years C++Now about using eigen with proto. If I remember correctly they talked explicitly about this problem and explained how they solved it. You can find the slides here, the video of the talk here and the source code here if you are interested. Unfortunately the audio is missing in the last minutes, but everything works fine until then. – Birdbath

Answer 1

This is my attempt at merging your approach with the solution linked in the comment. I have copied stored_result_expression, do_wrap_expression and wrap_expression from here. The changes I've made to either your code or the one from the talk are marked with //CHANGED.

#include <iostream>

#include <boost/fusion/container.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/void.hpp>
#include <boost/proto/proto.hpp>
#include <boost/ref.hpp>
#include <boost/type_traits/remove_const.hpp>
#include <boost/type_traits/remove_reference.hpp>
#include <boost/utility.hpp>

#include <Eigen/Dense>

namespace fusion = boost::fusion;
namespace mpl = boost::mpl;
namespace proto = boost::proto;

typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> matrix;

// Placeholders

const proto::terminal<mpl::int_<0> >::type I1 = {{}};
const proto::terminal<mpl::int_<1> >::type I2 = {{}};
const proto::terminal<mpl::int_<2> >::type I3 = {{}};

// Grammar

template<class Rule, class Callable = proto::callable>
struct External :
    proto::when<Rule, proto::external_transform> {};

struct matmul_transform : proto::callable {
    template<class Sig> struct result;

    template<class This, class Expr, class MatrixExpr1, class MatrixExpr2>
    struct result<This(Expr, MatrixExpr1, MatrixExpr2)> {
            typedef typename Eigen::MatrixBase<
                        typename Eigen::ProductReturnType<
                            typename boost::remove_const<typename boost::remove_reference<MatrixExpr1>::type>::type,
                            typename boost::remove_const<typename boost::remove_reference<MatrixExpr2>::type>::type
                        >::Type 
                    >::PlainObject&
                    type; //CHANGED - THIS IS THE TYPE THAT IS USED IN THE CODE OF THE TALK
    };

    template<class Expr, class MatrixExpr1, class MatrixExpr2>
    typename result<matmul_transform(Expr, MatrixExpr1, MatrixExpr2)>::type
    operator()(Expr& expr, const MatrixExpr1 &a, const MatrixExpr2 &b) const { //CHANGED - ADDED THE expr PARAMETER
            expr.value = a*b; 
            return expr.value; 
    }
};


struct MatmulGrammar;

struct InputPlaceholder : proto::terminal<proto::_> {};

struct MatrixMultiplication :
    proto::multiplies<MatmulGrammar, MatmulGrammar> {};

struct MatmulGrammar : proto::or_<
    External<InputPlaceholder>,
    External<MatrixMultiplication> > {};

struct matmul_transforms : proto::external_transforms<
    proto::when<MatrixMultiplication, matmul_transform(proto::_, MatmulGrammar(proto::_left), MatmulGrammar(proto::_right))>, //CHANGED - ADAPTED TO THE NEW SIGNATURE OF matmul_transform
    proto::when<InputPlaceholder, proto::functional::at(proto::_data, proto::_value)> > {};

// THE FOLLOWING CODE BLOCK IS COPIED FROM https://github.com/barche/eigen-proto/blob/master/eigen_calculator_solution.cpp
//----------------------------------------------------------------------------------------------
/// Wraps a given expression, so the value that it represents can be stored inside the expression itself
template<typename ExprT, typename ValueT>
struct stored_result_expression :
  proto::extends< ExprT, stored_result_expression<ExprT, ValueT> >
{
  EIGEN_MAKE_ALIGNED_OPERATOR_NEW

  typedef proto::extends< ExprT, stored_result_expression<ExprT, ValueT> > base_type;

  explicit stored_result_expression(ExprT const &expr = ExprT())
    : base_type(expr)
  {
  }

  /// Temporary storage for the result of the expression
  mutable ValueT value;
};

struct do_wrap_expression : proto::transform< do_wrap_expression >
{
  template<typename ExprT, typename StateT, typename DataT>
  struct impl : proto::transform_impl<ExprT, StateT, DataT>
  {
    typedef typename boost::result_of<MatmulGrammar(ExprT, StateT, DataT)>::type result_ref_type; //CHANGED - TO USE YOUR GRAMMAR
    typedef typename boost::remove_reference<result_ref_type>::type value_type;
    typedef typename boost::remove_const<typename boost::remove_reference<ExprT>::type>::type expr_val_type;
    typedef stored_result_expression<expr_val_type, value_type> result_type;

    result_type operator()(typename impl::expr_param expr, typename impl::state_param state, typename impl::data_param data)
    {
      return result_type(expr);
    }
  };
};

/// Wrap multiplies expressions so they can store a temporary result
struct wrap_expression :
  proto::or_
  <
    proto::terminal<proto::_>,
    proto::when
    <
      proto::multiplies<proto::_, proto::_>,
      do_wrap_expression(
        proto::functional::make_multiplies
        (
            wrap_expression(proto::_left), wrap_expression(proto::_right)
        ),
        proto::_state, //CHANGED - THESE EXTRA PARAMETERS ARE NEEDED TO CALCULATE result_ref_type IN do_wrap_expression
        proto::_env
      )
    >,
    proto::nary_expr< proto::_, proto::vararg<wrap_expression> >
  >
{
};
//--------------------------------------------------------------------------------------------------

int main() {
    matrix mat1(2,2), mat2(2,2), mat3(2,2), result(2,2);

    mat1 << 1, 1, 0, 1;
    mat2 << 1, 1, 0, 1;
    mat3 << 1, 1, 0, 1;

    MatmulGrammar mmg;
    wrap_expression wrap;

    //THIS WORKS:
     result = mmg( //THIS IS REALLY HORRIBLE, BUT IT WORKS. IT SHOULD PROBABLY BE HIDDEN BEHIND A FUNCTION
                wrap(
                    I1 * I2,
                    mpl::void_(),
                    ( proto::data = fusion::make_vector(boost::cref(mat1), boost::cref(mat2), boost::cref(mat3)),
                      proto::transforms = matmul_transforms() )
                ),
                mpl::void_(),
                ( proto::data = fusion::make_vector(boost::cref(mat1), boost::cref(mat2), boost::cref(mat3)),
                  proto::transforms = matmul_transforms() )
            );

    std::cout << result << std::endl;

    // THIS DOESN'T CRASH ANYMORE:
    result = mmg(
                wrap(
                    I1 * I2 * I3 * I1 * I2 * I3,
                    mpl::void_(),
                    ( proto::data = fusion::make_vector(boost::cref(mat1), boost::cref(mat2), boost::cref(mat3)),
                      proto::transforms = matmul_transforms() )
                ),
                mpl::void_(),
                ( proto::data = fusion::make_vector(boost::cref(mat1), boost::cref(mat2), boost::cref(mat3)),
                  proto::transforms = matmul_transforms() )
            );

    std::cout << result << std::endl;

    return 0;
}

Unsociable

Thank you very much for the useful links and for taking the time to edit my code! Your suggestions brought me quite a bit on the way! What I don't quite like about the solution from the C++Now talk is that it forces evaluation and creates a temporary at every matrix product operation. I think I managed to avoid this by storing shared pointers to the expression template objects instead of fully precomputed matrices in stored_result_expressions. Difficult to say if it's actually more efficient without benchmarking, of course. – Remember

Answer 2

Here's another solution that seems to work. Instead of interfering with proto's expression objects, it saves shared pointers to the intermediate Eigen objects as part of the state. Compared to the solution inspired by the C++Now talk, it has the following advantages:

• It doesn't force early evaluation of Eigen's expression templates.
• It requires fewer changes to the grammar, so it's less intrusive in the syntax of the domain-specific language.
• Responsibility for keeping the intermediate objects alive is assumed by the state, which is arguably where it belongs. In particular, I believe this makes the grammar thread-safe (if proto is).
• It returns an expression template, not a matrix. You should be safe even if you store this template in a variable and evaluate it later at your leisure, since all parts are included.

Drawbacks:

• Instead of returning a neat matrix, you get an unwieldy structure from which you have to extract the part that you're actually interested in.
• Temporary objects are allocated on the heap instead of the stack.
• You have to provide shared pointers to your matrices if you like it or not.

#include <iostream>

#include <boost/fusion/include/container.hpp>
#include <boost/fusion/include/join.hpp>
#include <boost/make_shared.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/void.hpp>
#include <boost/proto/proto.hpp>
#include <boost/ref.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/type_traits/remove_const.hpp>
#include <boost/type_traits/remove_reference.hpp>

#include <Eigen/Dense>

namespace fusion = boost::fusion;
namespace mpl = boost::mpl;
namespace proto = boost::proto;

typedef Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic> matrix;

// Placeholders

const proto::terminal<mpl::int_<0> >::type I1 = {{}};
const proto::terminal<mpl::int_<1> >::type I2 = {{}};
const proto::terminal<mpl::int_<2> >::type I3 = {{}};

// Grammar

template<class Rule, class Callable = proto::callable>
struct External :
    proto::when<Rule, proto::external_transform> {};

struct matmul_transform : proto::callable {
    template<class Sig> struct result;

    template<class This, class ExprList1, class ExprList2>
    struct result<This(ExprList1, ExprList2)> {
            typedef typename boost::remove_reference<
                    typename fusion::result_of::front<ExprList1>::type>::type::element_type M1;
            typedef typename boost::remove_reference<
                    typename fusion::result_of::front<ExprList2>::type>::type::element_type M2;
            typedef typename Eigen::ProductReturnType<
                    typename boost::remove_const<typename boost::remove_reference<M1>::type>::type,
                    typename boost::remove_const<typename boost::remove_reference<M2>::type>::type>::Type
                    product_return_type;
            typedef typename fusion::result_of::push_front<
                            const typename fusion::result_of::join<const ExprList1, const ExprList2>::type,
                            boost::shared_ptr<product_return_type> >::type
                    type;

    };

    template<class ExprList1, class ExprList2>
    typename result<matmul_transform(ExprList1, ExprList2)>::type
    operator()(const ExprList1 &a, const ExprList2 &b) const {
            typedef typename result<matmul_transform(ExprList1, ExprList2)>::product_return_type product_return_type;
            return push_front(join(a, b), boost::make_shared<product_return_type>(*front(a) * *front(b)));
    }
};

struct placeholder_transform : proto::callable {
    template<class Sig> struct result;

    template<class This, class Data, class Value>
    struct result<This(Data, Value)> {
            typedef typename boost::remove_const<typename boost::remove_reference<
                    typename fusion::result_of::at<Data, typename boost::remove_reference<Value>::type>::type>
                            ::type>::type ptr_type;
            typedef typename fusion::list<ptr_type> type;
    };

    template<class Data, class Value>
    typename result<placeholder_transform(Data, Value)>::type
    operator()(Data &data, Value value) const {
            return fusion::make_list(fusion::at<Value>(data));
    }
};

struct MatmulGrammar;

struct InputPlaceholder : proto::terminal<proto::_> {};

struct MatrixMultiplication :
    proto::multiplies<MatmulGrammar, MatmulGrammar> {};

struct MatmulGrammar : proto::or_<
    External<InputPlaceholder>,
    External<MatrixMultiplication> > {};

struct matmul_transforms : proto::external_transforms<
    proto::when<MatrixMultiplication, matmul_transform(MatmulGrammar(proto::_left), MatmulGrammar(proto::_right))>,
    proto::when<InputPlaceholder, placeholder_transform(proto::_data, proto::_value)> > {};

int main() {
    boost::shared_ptr<matrix> mat1 = boost::make_shared<matrix>(2,2);
    boost::shared_ptr<matrix> mat2 = boost::make_shared<matrix>(2,2);
    boost::shared_ptr<matrix> mat3 = boost::make_shared<matrix>(2,2);
    boost::shared_ptr<matrix> result = boost::make_shared<matrix>(2,2);

    *mat1 << 1, 1, 0, 1;
    *mat2 << 1, 1, 0, 1;
    *mat3 << 1, 1, 0, 1;

    MatmulGrammar mmg;

    // THIS WORKS:
    *result = *front(
            mmg(I1 * I2, mpl::void_(),
            (proto::data = fusion::make_vector(mat1, mat2, mat3),
             proto::transforms = matmul_transforms())));

    std::cout << *result << std::endl;

    // THIS WORKS, TOO:
    *result = *front(
            mmg(I1 * I2 * I3 * I3 * I2 * I1, mpl::void_(),
            (proto::data = fusion::make_vector(mat1, mat2, mat3),
             proto::transforms = matmul_transforms())));

    std::cout << *result << std::endl;

    return 0;
}

Remember