Armadillo conflicts with Boost Odeint: Odeint resizes the state vector to zero during integration

I just began to use Boost Odeint to integrate an ODE system. For convenience, I would like to use it with Armadillo, since both are modern C++ libraries with a convenient API. However, if I specify arma::vec as state type, it turns out immediately that in the first step of the integration, integrate_adaptive() resizes the state vector to 0x1. I post a trivial example here:

#include <iostream>
#include <armadillo>

#include <boost/numeric/odeint.hpp>

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

typedef vec state_type;

class harm_osc
{
    private:
    mat my_A;

    public:
        harm_osc(double gam)
        {
            my_A = { {0.0, 1.0}, {-gam*gam, 0.0} };
        }

        harm_osc()
        {
            my_A = { {0.0, 1.0}, {-1.0, 0.0} };
        }

        void operator() (const vec& x, vec& dxdt, const double t)
        {
            cout << "size of x: " << size(x) << endl;
            cout << "size of dxdt: " << size(dxdt) << endl;
            dxdt = my_A*x;
        }
};

class observer
{
    private:
        mat& my_states;
        vec& my_times;                        

    public:
        observer(mat& states, vec& times):
            my_states(states),
            my_times(times)
        {}    

        void operator() (const vec& x, double t)
        {
            my_states.insert_rows(my_states.n_rows, x);
            my_times.insert_rows(my_times.n_rows, t);
            cout << t << '\t';
            for(auto elem : x)
                cout << elem << '\t';
            cout << endl;
        }
};

typedef runge_kutta_cash_karp54<state_type> error_stepper_type;
typedef controlled_runge_kutta<error_stepper_type> controlled_stepper_type;

int main()
{
    state_type x = {0.0, 1.0};

    vec t;
    mat x_full;

    integrate_adaptive(make_controlled<error_stepper_type>(1e-5, 1e-5), harm_osc(1.0), x, 0.0, 200.0, 0.01, observer(x_full, t));
}

If I specify arma::vec::fixed<2> instead of arma::vec as state_type, this easy demo runs correctly. My problem is that, in the current project I'm working on, I don't know the size of my state vector at compile time, therefore I cannot fix it with a template parameter as mentioned just before.

Is there any solution to use Armadillo with Boost Odeint without fixing the size of the state vector at compile time? In the other parts of my project, I can use Armadillo properly. I would like to use it in my whole project. Right now, when integrating the ODE, I use Boost Ublas to define the ODE system.

Solution

Odeint stores internally several intermediate states. In the your case I think it does not know how to correctly resize the intermediate states. This can easily be fixed by introducing the correct resize adapters:

namespace boost { namespace numeric { namespace odeint {

template <>
struct is_resizeable<arma::vec>
{
    typedef boost::true_type type;
    const static bool value = type::value;
};

template <>
struct same_size_impl<arma::vec, arma::vec>
{
    static bool same_size(const arma::vec& x, const arma::vec& y)
    {
        return x.size() == y.size();   // or use .n_elem attributes
    }
};

template<>
struct resize_impl<arma::vec, arma::vec>
{
    static void resize(arma::vec &v1, const arma::vec& v2)
    {
        v1.resize(v2.size());     // not sure if this is correct for arma
    }
};

} } } // namespace boost::numeric::odeint

Have a look here: http://www.boost.org/doc/libs/1_63_0/libs/numeric/odeint/doc/html/boost_numeric_odeint/odeint_in_detail/state_types__algebras_and_operations.html