Eigen Vectors as ODEINT integrate parameters

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
    }
  }
}

Solution

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.