Additional support needed for "integrate" wrappers with Eigen matrices

[jefftrull]

I'm trying to convert a some existing code from using a std::vector wrapped with an Eigen Map to use an Eigen Matrix directly as a state variable with odeint. I found that the added support in 1.56 for Eigen was helpful but still incomplete. In particular:

1. The integrate wrapper functions expect the state type to have a member value_type, which Eigen::Matrix lacks, so I had to unroll them, creating a stepper typedef with the right value_type and calling integrate_adaptive.
2. Without an algebra_dispatcher specialization for Eigen::Matrix, range_algebra gets used in some places instead of vector_space_algebra
3. We need an overload for abs that ADL can pick up so default_operations::rel_error works properly.

A version of the harmonic oscillator example with the necessary changes made is here: https://gist.github.com/jefftrull/5625b77c0f86c439f29f

[headmyshoulder]

Hi Jeff,

point 2 and 3 should be easy to implement. Point 1 is a bit more complicated. But it should be doable with a meta function like

template< typename T > struct get_value_type { typedef typename T::value_type type; };

which then can be specialized for any Eigen type. I think I will come up with a patch very soon.

On 15.08.2014 08:03, Jeff Trull wrote:

I'm trying to convert a some existing code from using a std::vector wrapped with an Eigen Map to use an Eigen Matrix directly as a state variable with odeint. I found that the added support in 1.56 for Eigen was helpful but still incomplete. In particular:

1. The |integrate| wrapper functions expect the state type to have a member |value_type|, which Eigen::Matrix lacks, so I had to unroll them, creating a stepper typedef with the right value_type and calling |integrate_adaptive|.
2. Without an |algebra_dispatcher| specialization for Eigen::Matrix, |range_algebra| gets used in some places instead of |vector_space_algebra|
3. We need an overload for |abs| that ADL can pick up so |default_operations::rel_error| works properly.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/issues/139.

[mariomulansky]

Adding this metafunction just to get the integrate function running seems a bit like overkill to me. Maybe we can provide another version of integrate that takes the Time type as value_type, together with the standard way of deducing State::value_type if available this should cover 99.99% of all cases - and for whom this doesnt work should use integrate_adaptive directly...

[headmyshoulder]

@mariomulansky Yes, you are right. This metafunction looks easy at the first sight but turns out to be complicated, for examples when extracting the value type from std::vector< std::complex< double > >. I found a nice and easy solution: I added a second function template integrate has the value type as an explicit template parameter. An example is

integrate< double >( sys , x , t0 , t1 , dt );

I also fixed the other two points: I introduced the specialzation of the algebra dispatcher for Eigen matrices and an overload of the abs function.

And there is an include

boost/numeric/odeint/external/eigen/eigen.hpp

now which includes all headers needed for eigen.

@jefftrull can you check if everyhing works as expected?

Thanks you for reporting this,

Karsten

[jefftrull]

The integrate<value_type> overload is working for me now, with no more hacks. Thanks for the quick response!