Based on @Piotr Skotnicki answer here I have defined the following data type:
#ifndef TENSOR_HPP
#define TENSOR_HPP
#include <vector>
#include <array>
#include <boost/numeric/odeint.hpp>
namespace nex {
template <size_t...> struct seq {};
template <size_t N, size_t... Is> struct gen_seq : gen_seq<N-1, N-1, Is...> {};
template <size_t... Is> struct gen_seq<0, Is...> { using type = seq<Is...>; };
template<typename T, typename S>
class tensor;
template<typename T, size_t... Is>
class tensor<T, seq<Is...>>
{
typedef std::vector<T> vector;
public:
// constructor
tensor ();
tensor (decltype(Is)...);
// iterators
using iterator = typename vector::iterator;
using const_iterator = typename vector::const_iterator;
iterator begin()
{ return m_v.begin(); }
const_iterator begin() const
{ return m_v.begin(); }
iterator end()
{ return m_v.end(); }
const_iterator end() const
{ return m_v.end(); }
// resizing
size_t size() const;
void resize( const size_t );
// operators
const T& operator() (decltype(Is)...) const;
T& operator() (decltype(Is)...);
private:
// helper functions
template <size_t I>
auto index(const std::array<size_t, sizeof...(Is)>& a) const
-> typename std::enable_if<I == 0, size_t>::type
{
return a[I];
}
template <size_t I>
auto index(const std::array<size_t, sizeof...(Is)>& a) const
-> typename std::enable_if<I != 0, size_t>::type
{
return index<I-1>(a) * s[I] + a[I];
}
size_t mult(size_t N)
{
return N;
}
template <typename... S>
size_t mult(size_t N, S... Ns)
{
return N * mult(Ns...);
}
vector m_v;
const std::array<size_t, sizeof...(Is)> s;
}; // class tensor
template<typename T, size_t... Is
tensor<T, seq<Is...>> :: tensor ()
: m_v(), s{ { 0 } }
{}
template<typename T, size_t... Is>
tensor<T, seq<Is...>> :: tensor (decltype(Is)... size)
: m_v( mult( size... ) ), s{ { size... } }
{}
template<typename T, size_t... Is>
size_t tensor<T, seq<Is...>> :: size () const
{ return m_v.size(); }
template<typename T, size_t... Is>
void tensor<T, seq<Is...>> :: resize (const size_t n)
{ m_v.resize( n ); }
template<typename T, size_t... Is>
const T& tensor<T, seq<Is...>> :: operator() (decltype(Is)... n) const
{ return m_v.at(index<sizeof...(Is)-1>( { { n... } } ) ); }
template<typename T, size_t... Is>
T& tensor<T, seq<Is...>> :: operator() (decltype(Is)... n)
{ return m_v.at(index<sizeof...(Is)-1>( { { n... } } ) ); }
} // namespace nex
namespace boost { namespace numeric { namespace odeint {
template<typename T, size_t... Is>
struct is_resizeable< nex::tensor<T, nex::seq<Is...>> >
{
typedef boost::true_type type;
static const bool value = type::value;
};
} // namespace odeint
} // namespace numeric
} // namespace boost
template <typename T, std::size_t N>
using tensor = nex::tensor<T, typename nex::gen_seq<N>::type>;
#endif // TENSOR_HPP
Then I wanted to use this data type in boost::odeint to integrate very basic problem
#include "tensor.hpp"
typedef tensor<double, 1> state_type;
void lorenz( const state_type &x , state_type &dxdt , const double t )
{
dxdt(0) = x(0);
dxdt(1) = x(1);
dxdt(2) = x(2);
}
void write_lorenz( const state_type &x , const double t )
{
std::cout << t << '\t' << x(0) << '\t' << x(1) << '\t' << x(2) << std::endl;
}
using namespace boost::numeric::odeint;
int main()
{
state_type x(3);
x(0) = 1.0 ; x(1) = 10.0 ; x(2) = 10.0;
integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 1.0 , 0.1, write_lorenz );
}
and it worked fine, even for different set of equations (like a lorenz attractor). Unfortunatelly everyting went wrong when I wanted to use the same capacity of the tensor but defined differently
#include "tensor.hpp"
typedef tensor<double, 2> state_type;
void lorenz( const state_type &x , state_type &dxdt , const double t )
{
dxdt(0,0) = x(0,0);
dxdt(1,0) = x(1,0);
dxdt(2,0) = x(2,0);
}
void write_lorenz( const state_type &x , const double t )
{
std::cout << t << '\t' << x(0,0) << '\t' << x(1,0) << '\t' << x(2,0) << std::endl;
}
using namespace boost::numeric::odeint;
int main()
{
state_type x(3,1);
x(0,0) = 1.0 ; x(1,0) = 10.0 ; x(2,0) = 10.0;
integrate_const( runge_kutta4< state_type >() , lorenz , x , 0.0 , 1.0 , 0.1, write_lorenz );
}
You can run the code for yourself and check that x(1,0) x(2,0) are not changing at all. I though that appropriate definition of operator() is enough for odeint to work properly. I have rather poor knowledge of boost::odeint so maybe this is just a basic bug. I am counting on your help.
UPDATE
When I define tensor as a fixed size array using std::array it works fine. But then I don't need to add boost::is_resizeable. I think resizing needs to be defined differently (I don't really know exactly).
