ArmadilloはBoost Odeintと競合します：Odeintは統合中に状態ベクトルをゼロにサイズ変更します

Question

ODEをBoost Odeintで使用するようになりました。便宜上、便利なAPIを備えた最新のC++ライブラリであるため、Armadilloで使用したいと思います。しかし、arma::vecを状態タイプとして指定すると、統合の最初のステップでintegrate_adaptive()が状態ベクトルのサイズを0x1に変更することがすぐに分かります。私はここで簡単な例を投稿：

#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)); 
} 

私はstate_typeとしてarma::vec::fixed<2>の代わりarma::vecを指定した場合は、この簡単なデモが正常に実行されます。私の問題は、私が作業している現在のプロジェクトでは、コンパイル時の状態ベクトルのサイズがわからないため、直前に述べたテンプレートパラメータで修正できないということです。

ArmadilloをBoost Odeintで使用する方法はありますか？コンパイル時に状態ベクトルのサイズを固定しないでください。私のプロジェクトの他の部分では、私は適切にアルマジロを使うことができます。私は私のプロジェクト全体でそれを使用したいと思います。今、ODEを統合するとき、私はBoost Ublasを使用してODEシステムを定義します。

Answer 1

Odeintは内部的にいくつかの中間状態を格納します。あなたのケースでは、私はそれが中間状態を正しくサイズ変更する方法を知らないと思う。これは、簡単に正しいサイズ変更アダプターを導入することで固定することができます。

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::ver& x, const arma::vec& y) 
    { 
     return x.size() == y.size(); // not sure if this is correct for arma 
    } 
}; 

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 

はこちらをご覧：標準例（https://github.com/headmyshoulder/odeint-v2/blob/master/examples/my_vector.cpp）の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

Answer 2

少し微調整はあなたの作業コードを与える：

#include <boost/numeric/odeint.hpp> 
#include <armadillo> 
#include <iostream> 

namespace boost { 
namespace numeric { 
namespace odeint { 

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

} } } 

typedef arma::vec state_type; 

void lorenz(const state_type &x , state_type &dxdt , const double t) 
{ 
const double sigma(10.0); 
const double R(28.0); 
const double b(8.0/3.0); 

dxdt[0] = sigma * (x[1] - x[0]); 
dxdt[1] = R * x[0] - x[1] - x[0] * x[2]; 
dxdt[2] = -b * x[2] + x[0] * x[1]; 
} 

using namespace boost::numeric::odeint; 

int main() 
{ 
state_type x(3); 
x[0] = 5.0 ; x[1] = 10.0 ; x[2] = 10.0; 

// make sure resizing is ON 
BOOST_STATIC_ASSERT(is_resizeable<state_type>::value == true); 
// no further work is required 
std::cout << "initial x : " << x << std::endl; 
integrate_const(runge_kutta4<state_type>() , lorenz , x , 
0.0 , 10.0 , 0.1); 
std::cout << "final x : " << x << std::endl; 
return 0; 
}