A different way of bridging Eigen and odeint

[cshelton]

This is a single header file that adds some ops to the Eigen namespace and a vector_space_algebra to the boost::numeric::odeint namespace so that Eigen matrices can be used more or less "out of the box" with odeint. It works for me with dynamic allocation and static allocation Eigen matrices and vectors.

[headmyshoulder]

@cshelton Short question: It is an addition to eigen_resize in order to work with the controlled stepper?

[cshelton]

I'm not sure I understand your question fully. My code doesn't reference resize. Eigen matrices automatically resize if necessary when assigned to. Is there a line in the code in particular that illustrates your question?

[headmyshoulder]

On 03/25/2013 07:05 PM, Christian Shelton wrote:

I'm not sure I understand your question fully. My code doesn't reference resize. Eigen matrices automatically resize if necessary when assigned to. Is there a line in the code in particular that illustrates your question?

Ok, I understand. You are correct. Resizing happens if the assign operator is called.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/pull/74#issuecomment-15410258.

[headmyshoulder]

Hi Christian, I tried to include your header and build a unit test, see https://github.com/headmyshoulder/odeint-v2/blob/eigen/libs/numeric/odeint/test_external/eigen/runge_kutta_dopri5.cpp .

Unfortunately it does not compile, because the Eigen traits are not found. Any ideas? I think there is only one header mssing, but I am not so familiar with Eigen. Or I need a newer version of Eigen.

[cshelton]

Your code compiled for me (I had to remove the eigen_resize.h #include, as I don't use that header, but I don't think that's a problem). I am using Eigen3 and that might be the difference.

I notice that you are also testing with Eigen arrays. I don't think my code handles that case. I've found the distinction between matrices and arrays annoying in Eigen and have largely tried to stick with matrices. If the code doesn't work with arrays, I can try to modify it so that it does.

[headmyshoulder]

That is strange. I obtain errors like this:

../../../../../boost/numeric/odeint/external/eigen/eigen_algebra.hpp: In function ‘const Eigen::CwiseUnaryOp<Eigen::internal::scalar_add_op<typename Eigen::internal::traits::Scalar>, const Derived> Eigen::operator+(const Eigen::MatrixBase&, const typename Eigen::internal::traits::Scalar&)’: ../../../../../boost/numeric/odeint/external/eigen/eigen_algebra.hpp:43:23: error: ‘traits’ was not declared in this scope ../../../../../boost/numeric/odeint/external/eigen/eigen_algebra.hpp:43:23: note: suggested alternatives: /home/karsten/boost/boost_1_53_0/boost/fusion/support/tag_of_fwd.hpp:14:5: note: ‘boost::fusion::traits’ /home/karsten/boost/eigen/Eigen/src/Core/util/ForwardDeclarations.h:31:29: note: ‘Eigen::internal::traits’ ../../../../../boost/numeric/odeint/external/eigen/eigen_algebra.hpp:43:31: error: template argument 1 is invalid ../../../../../boost/numeric/odeint/external/eigen/eigen_algebra.hpp:42:18: error: expected primary-expression before ‘(’ token ../../../../../boost/numeric/odeint/external/eigen/eigen_algebra.hpp:43:41: error: expected primary-expression before ‘,’ token ../../../../../boost/numeric/odeint/external/eigen/eigen_algebra.hpp:44:19: error: expected primary-expression before ‘const’ /home/karsten/boost/boost_1_53_0/boost/test/unit_test_log.hpp: At global scope:

Do you have an idea what is happening here and how I can fix this. I use Eigen version 3.0.6

[cshelton]

Your line numbers don't match up with mine, so I expect you've adjusted the file at least slightly (maybe just added comments?). I cannot reproduce your errors on my system (with Eigen 3.0.6) using the file you attached and the current version of odeint and the eigen_algebra.hpp that I included. I'm using gcc 4.7.2 if that matters.

I notice that your errors do mention Eigen::internal::traits as a possibility, so I'm not sure how to parse the error. Can you supply the whole source tree that you have?

As a side note: this is one (of many) areas where c++11 helps: I could just use auto/decltype for the return type and avoid this whole mess (and be more robust to changes in Eigen).

[headmyshoulder]

You can checkout the eigen branch:

git branch -b eigen origin/eigen

or

git clone -b eigen git@github.com:headmyshoulder/odeint-v2.git odeint-v2-eigen

The executable which fails is in libs/numeric/odeint/test_external/runge_kutta_dopri5.cpp

I also use g++-4.7 and I also tried with std=c++11, without success. Thank you for your help.

[cshelton]

It looks like you edited the file to remove the second operator+. In doing that (or perhaps something else), two lines were inserted at lines 43 and 44 Remove those two lines and you'll have the correct code.

[headmyshoulder]

On 03/27/2013 10:54 PM, Christian Shelton wrote:

It looks like you edited the file to remove the second operator+. In doing that (or perhaps something else), two lines were inserted at lines 43 and 44 Remove those two lines and you'll have the correct code.

Ok, I got it. Thank you again. I will try to finish the tests tomorrow and then merge everything.

— Reply to this email directly or view it on GitHub https://github.com/headmyshoulder/odeint-v2/pull/74#issuecomment-15555610.

[headmyshoulder]

I used another branch from odeint, but all your changes are now in the trunk. I also needed the resizing stuff, since it allows you to write equations like

dxdt = 1.0

inside the system function.

Thank you for your fixes.

[luminoctum]

Dear cshelton,

Could you give more instructions on how to use the head file? I've copied your head file but it fails to compile and outputs a bunch of wield error message.

[cshelton]

@luminoctum Could you give more details? What version of Eigen are you using? What are the compilation errors you are getting? The usage is reasonably straight-forward, I think: you need to include "eigen_algebra.hpp" (order of includes shouldn't matter -- but let me know if that's not true). Then you can use Eigen vectors and matrices as the types for odeint's states and derivatives.

[luminoctum]

Thanks for the quick response.

I'm using boost version 1.53 and Eigen version 3.1.3

My test code is a simple test for odeint:

include

include <boost/numeric/odeint.hpp>

include <boost/numeric/odeint/external/eigen/eigen_algebra.hpp>

using namespace std; using namespace boost::numeric::odeint; using namespace Eigen;

const double sigma = 10.0; const double R = 28.0; const double b = 8.0 / 3.0;

typedef Vector3f state_type;

void lorenz( const state_type &x, state_type &dxdt, double){ dxdt[0] = sigma * (x[1] - x[0]); dxdt[1] = sigma * x[0] - x[1] - x[0]x[2]; dxdt[2] = -b * x[2] + x[0] \ x[1]; }

void write_lorenz( const state_type &x, const double t){ cout << t << "\t" << x[0] << "\t" << x[1] << "\t" << x[2] << endl; }

int main(){ state_type x(10.0, 1.0, 1.0); integrate(lorenz, x, 0.0, 25.0, 0.1, write_lorenz); cout << x << endl; }

I ended up with long error messages. First few lines are:

/home/cli/lib/boost_1_53_0/boost/mpl/eval_if.hpp: In instantiation of ‘boost::mpl::eval_if_c<true, boost::range_const_iterator<Eigen::Matrix<float, 3, 1, 0, 3, 1> >, boost::range_mutable_iterator<const Eigen::Matrix<float, 3, 1, 0, 3, 1> > >’: /home/cli/lib/boost_1_53_0/boost/range/iterator.hpp:63: instantiated from ‘boost::range_iterator<const Eigen::Matrix<float, 3, 1, 0, 3, 1> >’ /home/cli/lib/boost_1_53_0/boost/numeric/odeint/algebra/range_algebra.hpp:52: instantiated from ‘static void boost::numeric::odeint::range_algebra::for_each3(S1&, S2&, S3&, Op) [with S1 = Eigen::Matrix<float, 3, 1, 0, 3, 1>, S2 = const Eigen::Matrix<float, 3, 1, 0, 3, 1>, S3 = const Eigen::Matrix<float, 3, 1, 0, 3, 1>, Op = boost::numeric::odeint::default_operations::rel_error]’ /home/cli/lib/boost_1_53_0/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:80: instantiated from ‘Value boost::numeric::odeint::default_error_checker<Value, Algebra, Operations>::error(Algebra&, const State&, const Deriv&, Err&, Time) const [with State = Eigen::Matrix<float, 3, 1, 0, 3, 1>, Deriv = Eigen::Matrix<float, 3, 1, 0, 3, 1>, Err = Eigen::Matrix<float, 3, 1, 0, 3, 1>, Time = double, Value = double, Algebra = boost::numeric::odeint::range_algebra, Operations = boost::numeric::odeint::default_operations]’ ...

I appreciate your help

[luminoctum]

Sorry I'm new to github. The "# include" seems not working.

They are actually:
"# include iostream" "# include boost/numeric/odeint.hpp" "# include boost/numeric/odeint/external/eigen/eigen_algebra.hpp"

[cshelton]

Eigen vectors are not considered containers (to be iterated over by odeint) in my version. Rather, I tell odeint to use Eigen's versions of addition and such. This makes things a little faster (at least in theory). However, this also means that you cannot use the default "integrate" function. You have to tell the integration what stepper you want and the stepper has to be associated with vector_space_algebra.

Try replacing the line

integrate(lorenz, x, 0.0, 25.0, 0.1, write_lorenz);

with

typedef runge_kutta_dopri5<state_type, double, state_type, double, vector_space_algebra> stepper;
integrate_adaptive(stepper(), lorenz, x, 0.0, 25.0, 0.1, write_lorenz);

The "integrate functions" section of the documentation gives a few more details about the type of integrators available. The "vector space algebra" part of the docs shows an example for ublas and a user-defined Point class. This is similar.

[luminoctum]

Ha, It works！ Thank you very much for the help!

[headmyshoulder]

@cshelton Thanks for your help and explanations.