Eigen Vectors as ODEINT integrate parameters

Question

I am trying to port the Harmonic Oscillator tutorial from ODEINT to Eigen, so that I could use VectorXd for state vectors. I cannot, however, make it compile.

I've been following some questions, for instance this is the closest except that I don't use any stepper here.

This is the code:

#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <boost/numeric/odeint.hpp>

typedef Eigen::VectorXd state_type;
// was vector<double>

const double gam = 0.15;

/* The rhs of x' = f(x) defined as a class */
class harm_osc
{

public:

    void operator() ( const state_type &x , state_type &dxdt , const double /* t */ )
    {
        dxdt(0) =  x(1);
        dxdt(1) = -x(0) - gam*x(1);
//        dxdt[0] = x[1];
//        dxdt[1] = -x[0] - gam*x[1];
    }
};
/* The rhs of x' = f(x) */
void harmonic_oscillator(const state_type &x, state_type &dxdt, const double /* t */ )
{
    dxdt(0) =  x(1);
    dxdt(1) = -x(0) - gam*x(1);
//    dxdt[0] = x[1];
//    dxdt[1] = -x[0] - gam*x[1];
}

void printer(const state_type &x , const double t)
{
//    std::cout << t << "," << x[0] << "," << x[1] << std::endl;
    std::cout << t << "," << x(0) << "," << x(1) << std::endl;
};

int main(int argc, const char * argv[])
{
    state_type x(2);
    x(0) = 1.0;
    x(1) = 0.0;
//    x[0] = 1.0;
//    x[1] = 0.0;

    std::cout << ">>>> FUNCTION" << std::endl;
//    boost::numeric::odeint::integrate(harmonic_oscillator, x, 0.0, 10.0, 0.1, printer);
//    boost::numeric::odeint::runge_kutta4<state_type, double, state_type, double, boost::numeric::odeint::vector_space_algebra> stepper();
    boost::numeric::odeint::integrate<double, decltype(harmonic_oscillator), state_type, double, decltype(printer)>(harmonic_oscillator, x, 0.0, 10.0, 0.1, printer);

    std::cout << ">>>> CLASS" << std::endl;
    x(0) = 1.0;
    x(1) = 0.0;
//    x[0] = 1.0;
//    x[1] = 0.0;
    harm_osc ho;
    boost::numeric::odeint::integrate<double, decltype(harmonic_oscillator), state_type, double, decltype(printer)>(ho, x, 0.0, 10.0, 0.1, printer);

    return 0;
}

The compiler complains about No matching function for call to 'begin' in range_algebra.hpp from ODEINT in both class and function versions of integrate. As a matter of fact, Eigen matrices have no begin/end.

I've tried to play with the template parameters (as you see) with no avail.

Any hint?

Assertion failed using Eigen from repo

Using the latest Eigen from the repo (not the latest stable version), I can, as suggested, compile it and run. However, it fails an assertion in the integrate call tree:

Assertion failed: (index >= 0 && index < size()), function operator(), file eigen/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h, line 427.

The call that fails is dxdt(0) = x(1); and subsequently in DenseCoeffsBase.h:

    EIGEN_DEVICE_FUNC
    EIGEN_STRONG_INLINE Scalar&
    operator()(Index index)
    {
      eigen_assert(index >= 0 && index < size()); // <---- INDEX IS 0, SIZE IS 0
      return coeffRef(index);
    }

Is it possible that ODEINT is trying to default-construct a VectorXd? I followed the path to my ODE call and dxdt is actually NULL:

(lldb) e dxdt
(state_type) $1 = {
  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
    m_storage = {
      m_data = 0x0000000000000000
      m_rows = 0
    }
  }
}

What is worse is that when using resizeLike to allow resizing dxdt, in the second step (so the first real computation of integrate) I have a x with NULL values:

(lldb) e dxdt
(state_type) $0 = {
  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
    m_storage = {
      m_data = 0x0000000000000000
      m_rows = 0
    }
  }
}
(lldb) e x
(state_type) $1 = {
  Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > = {
    m_storage = {
      m_data = 0x0000000000000000
      m_rows = 0
    }
  }
}

Try the master branch of Eigen, `begin()`/`end()` have been introduced a while ago. — chtz, Apr 14 '20 at 17:52
And N.B. `dxdt[0] = x[1];` etc, works for Eigen vector types (not for matrices, though), so no need to rewrite these lines. — chtz, Apr 14 '20 at 17:55
The requirements for `boost::odeint` StateType are described in detail [here](https://www.boost.org/doc/libs/1_72_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html); TL;DR: for non-fixed size types, apart from providing `begin()` and `end()` your type must also support `boost::size( x )` and `x.resize( boost::size( y ) )` — pptaszni, Apr 14 '20 at 20:31
Mmmh... so my choices are using a non-stable branch, or if I can, make a wrapper? @chtz do you have any info on why were begin/end excluded in the release? — senseiwa, Apr 15 '20 at 09:51
It just has been a while since the last release -- that's mostly a time issue of the maintainers. You can of course write a wrapper, which just returns `this->data()` for `.begin()` and `this->data()+this->size()` for `end()` (plus const/non-const variants and some typedefs) — chtz, Apr 15 '20 at 10:42
@chtz I've updated the question. Now it compiles but Eigen crashes with an ominous zero size of the vector. — senseiwa, Apr 16 '20 at 11:32
If this fails in your code, you could just `dxdt.resize(2)` (not sure why it would work with `std::vector` but not with `Eigen::VectorXd` in that case. But actually, if the size is known at compiletime to be 2, just use `Eigen::Vector2d`. — chtz, Apr 16 '20 at 16:09
@chtz unfortunately, I know just at runtime what vector space I'll be using. — senseiwa, Apr 20 '20 at 13:45
Did you try `dxdt.resize(2)` or `dxdt.resizeLike(x)` before setting its members? — chtz, Apr 20 '20 at 13:48
@chtz interestingly, in the second step I have both `x` and `dxdt` invalid, and the assertion fires again. I've updated the question. That's weird. On other questions seemed that it should work out of the box... — senseiwa, Apr 21 '20 at 10:25
Yea, no idea what happens inside odeint. And I don't want to debug that. — chtz, Apr 23 '20 at 12:28

score 1 · Accepted Answer · answered Jun 20 '20 at 05:48

I found that ODEINT actually works fine with Eigen... only it is not documented, as far as I can see.

Digging around ODEINT's code, I have found a promising eigen.hpp header deep in the external directory.

And lo and behold, it works flawlessly:

#include <fstream>
#include <iostream>
#include <vector>

#include <boost/numeric/odeint/external/eigen/eigen.hpp>
#include <Eigen/Eigenvalues>

#define FMT_HEADER_ONLY
#include "fmt/core.h"
#include "fmt/format.h"
#include "fmt/format-inl.h"
#include "fmt/printf.h"

using namespace std;

int main(int argc, char *argv[])
{
    Eigen::VectorXd v;

    v.resize(2);

    typedef Eigen::VectorXd state_type;

    const double gam = 0.15;

    v(0) = 1.0;
    v(1) = 1.1;

    auto harmonic_oscillator = [&](const state_type &x, state_type &dxdt, const double t)
    {
        dxdt[0] = x[1];
        dxdt[1] = -x[0] - gam*x[1];
    };

    auto printer = [&](const state_type &x, const double t)
    {
        fmt::print(out, "time: {} state: {}\n", t, x);
    };

    odeint::integrate(harmonic_oscillator, v, 0.0 , 10.0 , 0.01, printer);

    return 0;
}

Hope it helps others.

Eigen Vectors as ODEINT integrate parameters

Assertion failed using Eigen from repo

1 Answers1