I'm trying to write a general-purpose integrator in C++ leveraging Odeint and Eigen, but im stuck in trying to make Odeint happy about my custom State vector class that inherits from Eigen. To provide a minimal reproducible code, first im going to define my custom vector class, thats actually copy-pasted from my project:
#include <boost/numeric/odeint.hpp>
#include <eigen3/Eigen/Dense>
#include <functional>
#include <iostream>
#include <vector>
using namespace Eigen;
namespace odeint = boost::numeric::odeint;
template <size_t N>
using VReal = Eigen::Vector<double, N>;
using VReal3 = VReal<3>;
template <size_t N>
class State : public Eigen::Vector<double, 2 * N> {
private:
double t = 0;
public:
/////////////////////// Constructors ///////////////////////
/// ///
// default constructor
State() : Eigen::Vector<double, 2 * N>() { this->setZero(); }
// constructor with constant value
State(const double val) : Eigen::Vector<double, 2 * N>() {
this->setConstant(val);
}
// Constructor with initializer list
State(std::initializer_list<double> init_list)
: Eigen::Vector<double, 2 * N>() {
assert(init_list.size() ==
2 * N); // Ensure the initializer list is the correct size
int i = 0;
for (double val : init_list) {
(*this)[i++] = val;
}
}
// Constructor with VReal of 2N elements
State(const Eigen::Vector<double, 2 * N>& vec)
: Eigen::Vector<double, 2 * N>(vec) {}
// Constructor with 2 VReal of N rows
State(VReal<N> head, VReal<N> tail) {
this->template head<N>() = head;
this->template tail<N>() = tail;
}
// move constructor
State(State&& other) noexcept
: Eigen::Vector<double, 2 * N>(std::move(other)) {
t = other.t;
}
// copy constructor
State(const State& other)
: Eigen::Vector<double, 2 * N>(other), t(other.t) {}
/// ///
////////////////////////////////////////////////////////////
/////////////////////////// Operators ///////////////////////////
/// ///
// operator=
State& operator=(const State& other) {
if (this != &other) {
this->Eigen::Vector<double, 2 * N>::operator=(other);
t = other.t;
}
return *this;
}
// move operator=
State& operator=(State&& other) noexcept {
if (this != &other) {
this->Eigen::Vector<double, 2 * N>::operator=(std::move(other));
t = other.t;
}
return *this;
}
// Copy with Eigen::Matrix
State& operator=(
const Eigen::Matrix<double, 2 * N, 1, 0, 2 * N, 1>& other) {
Eigen::Vector<double, 2 * N>::operator=(other);
return *this;
}
//// Operators (+=, -=, *=, /=) for State (scalar) and State (vector)
////
// Vectorial
State operator+(const State& other) {
return State(this->Eigen::Vector<double, 2 * N>::operator+(other));
}
State& operator+=(const State& other) {
this->Eigen::Vector<double, 2 * N>::operator+=(other);
return *this;
}
State& operator-(const State& other) {
return State(this->Eigen::Vector<double, 2 * N>::operator-(other));
}
State& operator-=(const State& other) {
this->Eigen::Vector<double, 2 * N>::operator-=(other);
return *this;
}
// Scalar
State& operator*(const double scalar) {
return State(this->Eigen::Vector<double, 2 * N>::operator*(scalar));
}
State& operator*=(const double scalar) {
this->Eigen::Vector<double, 2 * N>::operator*=(scalar);
return *this;
}
State& operator/(const double scalar) {
return State(this->Eigen::Vector<double, 2 * N>::operator/(scalar));
}
State& operator/=(const double scalar) {
this->Eigen::Vector<double, 2 * N>::operator/=(scalar);
return *this;
}
// Vector product
State& operator*(const State& other) {
return State(this->Eigen::Vector<double, 2 * N>::operator*(other));
}
State& operator*=(const State& other) {
this->Eigen::Vector<double, 2 * N>::operator*=(other);
return *this;
}
/// ///
//////////////////////////////////////////////////////////////////
/////////////////////////// Methods ///////////////////////////
VReal<N> X() const { return this->template head<N>(); }
VReal<N> X_dot() const { return this->template tail<N>(); }
};
And some more classes to do the integration steps, these follow the same philosophy of my project but they are minimal, and some odd choices you see here make sense in the whole code:
template <size_t N>
class Model {
// Ballistic model without drag force
private:
public:
Model() {}
void operator()(const State<N>& x, State<N>& dxdt, double t) {
dxdt.X() = x.X_dot();
dxdt.X_dot() = Eigen::Vector<double, N>::Zero();
}
};
template <size_t N>
class RK4Stepper {
public:
RK4Stepper() {}
~RK4Stepper() {}
void do_step(Model<N> h, State<N>& in, double t, State<N>& out, double dt) {
stepper.do_step(h, in, t, out, dt);
}
private:
odeint::runge_kutta4<State<N>&, double, State<N>&, double,
odeint::vector_space_algebra>
stepper;
};
template <size_t N>
class Sim {
private:
RK4Stepper<N> stepper;
public:
Sim() {}
void run(Model<N> h, State<N> S0) {
std::cout << "Sim::run() called" << std::endl;
// std::reference_wrapper<Model<N>> h_ref = *h;
for (int i = 0; i < 10; i++) {
std::cout << "i = " << i << std::endl;
State<N> S1;
stepper.do_step(h, S0, 0.0, S1, 0.1);
S0 = S1;
}
}
};
And here is the main function:
int main() {
Model<3> model;
Sim<3> sim;
sim.run(model, State<3>(Eigen::Vector3d(1.0, 2.0, 3.0),
Eigen::Vector3d(4.0, 5.0, 6.0)));
return 0;
}
Now if I try to compile these, I get the following errors:
In file included from /usr/include/boost/numeric/odeint/stepper/runge_kutta4.hpp:27,
from /usr/include/boost/numeric/odeint.hpp:29,
from main.cpp:2:
/usr/include/boost/numeric/odeint/stepper/explicit_generic_rk.hpp: In instantiation of ‘boost::numeric::odeint::explicit_generic_rk<StageCount, Order, State, Value, Deriv, Time, Algebra, Operations, Resizer>::explicit_generic_rk(const coef_a_type&, const coef_b_type&, const coef_c_type&, const algebra_type&) [with long unsigned int StageCount = 4; long unsigned int Order = 4; State = State<3>&; Value = double; Deriv = State<3>&; Time = double; Algebra = boost::numeric::odeint::vector_space_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; coef_a_type = boost::fusion::vector<const boost::array<double, 1>, const boost::array<double, 2>, const boost::array<double, 3> >; coef_b_type = boost::array<double, 4>; coef_c_type = boost::array<double, 4>; algebra_type = boost::numeric::odeint::vector_space_algebra]’:
/usr/include/boost/numeric/odeint/stepper/runge_kutta4.hpp:141:81: required from ‘boost::numeric::odeint::runge_kutta4<State, Value, Deriv, Time, Algebra, Operations, Resizer>::runge_kutta4(const algebra_type&) [with State = State<3>&; Value = double; Deriv = State<3>&; Time = double; Algebra = boost::numeric::odeint::vector_space_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; algebra_type = boost::numeric::odeint::vector_space_algebra]’
main.cpp:154:18: required from ‘RK4Stepper<N>::RK4Stepper() [with long unsigned int N = 3]’
main.cpp:173:11: required from ‘Sim<N>::Sim() [with long unsigned int N = 3]’
main.cpp:192:12: required from here
/usr/include/boost/numeric/odeint/stepper/explicit_generic_rk.hpp:141:64: error: use of deleted function ‘boost::numeric::odeint::state_wrapper<State<3>&, void>::state_wrapper()’
141 | : stepper_base_type( algebra ) , m_rk_algorithm( a , b , c )
| ^
In file included from /usr/include/boost/numeric/odeint/util/ublas_wrapper.hpp:33,
from /usr/include/boost/numeric/odeint.hpp:25:
/usr/include/boost/numeric/odeint/util/state_wrapper.hpp:36:8: note: ‘boost::numeric::odeint::state_wrapper<State<3>&, void>::state_wrapper()’ is implicitly deleted because the default definition would be ill-formed:
36 | struct state_wrapper
| ^~~~~~~~~~~~~
/usr/include/boost/numeric/odeint/util/state_wrapper.hpp:36:8: error: uninitialized reference member in ‘struct boost::numeric::odeint::state_wrapper<State<3>&, void>’
/usr/include/boost/numeric/odeint/util/state_wrapper.hpp:40:7: note: ‘State<3>& boost::numeric::odeint::state_wrapper<State<3>&, void>::m_v’ should be initialized
40 | V m_v;
| ^~~
/usr/include/boost/numeric/odeint/stepper/explicit_generic_rk.hpp:141:64: error: use of deleted function ‘boost::numeric::odeint::state_wrapper<State<3>&, void>::state_wrapper()’
141 | : stepper_base_type( algebra ) , m_rk_algorithm( a , b , c )
| ^
/usr/include/boost/numeric/odeint/stepper/explicit_generic_rk.hpp:141:64: error: use of deleted function ‘boost::numeric::odeint::state_wrapper<State<3>&, void>::state_wrapper()’
In file included from /usr/include/boost/numeric/odeint/stepper/euler.hpp:23,
from /usr/include/boost/numeric/odeint.hpp:27:
/usr/include/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp: In instantiation of ‘boost::numeric::odeint::explicit_stepper_base<Stepper, Order, State, Value, Deriv, Time, Algebra, Operations, Resizer>::explicit_stepper_base(const algebra_type&) [with Stepper = boost::numeric::odeint::explicit_generic_rk<4, 4, State<3>&, double, State<3>&, double, boost::numeric::odeint::vector_space_algebra, boost::numeric::odeint::default_operations, boost::numeric::odeint::initially_resizer>; short unsigned int Order = 4; State = State<3>&; Value = double; Deriv = State<3>&; Time = double; Algebra = boost::numeric::odeint::vector_space_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; algebra_type = boost::numeric::odeint::vector_space_algebra]’:
/usr/include/boost/numeric/odeint/stepper/explicit_generic_rk.hpp:141:64: required from ‘boost::numeric::odeint::explicit_generic_rk<StageCount, Order, State, Value, Deriv, Time, Algebra, Operations, Resizer>::explicit_generic_rk(const coef_a_type&, const coef_b_type&, const coef_c_type&, const algebra_type&) [with long unsigned int StageCount = 4; long unsigned int Order = 4; State = State<3>&; Value = double; Deriv = State<3>&; Time = double; Algebra = boost::numeric::odeint::vector_space_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; coef_a_type = boost::fusion::vector<const boost::array<double, 1>, const boost::array<double, 2>, const boost::array<double, 3> >; coef_b_type = boost::array<double, 4>; coef_c_type = boost::array<double, 4>; algebra_type = boost::numeric::odeint::vector_space_algebra]’
/usr/include/boost/numeric/odeint/stepper/runge_kutta4.hpp:141:81: required from ‘boost::numeric::odeint::runge_kutta4<State, Value, Deriv, Time, Algebra, Operations, Resizer>::runge_kutta4(const algebra_type&) [with State = State<3>&; Value = double; Deriv = State<3>&; Time = double; Algebra = boost::numeric::odeint::vector_space_algebra; Operations = boost::numeric::odeint::default_operations; Resizer = boost::numeric::odeint::initially_resizer; algebra_type = boost::numeric::odeint::vector_space_algebra]’
main.cpp:154:18: required from ‘RK4Stepper<N>::RK4Stepper() [with long unsigned int N = 3]’
main.cpp:173:11: required from ‘Sim<N>::Sim() [with long unsigned int N = 3]’
main.cpp:192:12: required from here
/usr/include/boost/numeric/odeint/stepper/base/explicit_stepper_base.hpp:94:42: error: use of deleted function ‘boost::numeric::odeint::state_wrapper<State<3>&, void>::state_wrapper()’
94 | : algebra_stepper_base_type( algebra )
| ^
Now, I get what the compiler is trying to tell me. Odeint can't figure out how to digest my custom State class and doesn't know how to perform algebraic operations on it. But as in this answer, I thought that all I needed to do is to specify vector_space_algebra in the stepper template definition, but that obviously doesn't seem to work. So what should I do to make Odeint happy about my custom class?
As per the comment from chtz, changing
odeint::runge_kutta4<State<N>&, double, State<N>&, double, odeint::vector_space_algebra> stepper;
to
odeint::runge_kutta4<State<N>, double, State<N>, double, odeint::vector_space_algebra> stepper;
fixes the error.