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.
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).
Indeed, the problem lies in the resizing. Currently, resizing your tensor only changes the size of the underlying vector m_v
. Although this ensures correct memory allocation, it is not sufficient. You also need to set the dimensions correctly, i.e. the member s
. You can also think about it that way: it is impossible to resize a multi-dimensional tensor given just a single integer, you need to set the size of each dimension. The only exception is a 1D tensor, for which things worked fine for you, for any higher dimension this had to fail.
The following tensor.hpp addresses this problem and seems to work well in my quick test:
#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 tensor<T, seq<Is...>>& );
// 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;
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 tensor<T, seq<Is...>> &x)
{
m_v.resize( x.size() );
s = x.s;
}
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;
};
template<typename T, size_t... Is>
struct resize_impl< nex::tensor<T, nex::seq<Is...>>, nex::tensor<T, nex::seq<Is...>> >
{
typedef nex::tensor<T, nex::seq<Is...>> state_type;
static void resize( state_type &x1 , const state_type &x2 )
{
x1.resize(x2);
}
};
} // 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
Collected from the Internet
Please contact [email protected] to delete if infringement.
Comments