使用Armadillo和boost :: numeric :: odeint进行模板实例化

Posted

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了使用Armadillo和boost :: numeric :: odeint进行模板实例化相关的知识，希望对你有一定的参考价值。

我试图将boost :: numeric :: odeint与我自己的System类的实现结合起来（参见System.hpp）。

（模板）System对象在BatchFilter类方法中使用，如下所示：

# BatchFilter.cpp
# template <typename state_type> BatchFilter class {...}

System<state_type> dynamics(this -> args,
            N_true,
            this -> true_dynamics_fun );

        typedef boost::numeric::odeint::runge_kutta_cash_karp54< state_type > error_stepper_type;
        auto stepper = boost::numeric::odeint::make_controlled<error_stepper_type>( 1.0e-13 , 1.0e-16 );

        auto tbegin = T_obs.begin();
        auto tend = T_obs.end();

        boost::numeric::odeint::integrate_times(stepper, dynamics, X0_true_copy, tbegin, tend,0.1,
            Observer::push_back_state(this -> true_state_history));

BatchFilter是一个模板派生类，我在BatchFilter.cpp的底部显式实例化。这两个显式实例是

template class BatchFilter< arma::vec >; 
template class BatchFilter< arma::vec::fixed<2> >;

Base类也是显式实例化的。我必须使用arma::vec::fixed<2>，因为在arma::vec中使用odeint会导致运行时崩溃，因为状态将没有适当的大小。

什么不起作用？

使用arma::vec::fixed<2>的实例化失败，而arma::vec的实例成功。编译器抱怨非法绑定：

cannot bind non-const lvalue reference of type'arma::Col<double>::fixed<2>&' to an rvalue of type'arma::Col<double>::fixed<2>' sys( x , m_dxdt.m_v ,t );

令我困惑的是，当BatchFilter< arma::vec >的显式实例化成功而BatchFilter< arma::vec::fixed<2> >失败时，一切正常。

对于发生了什么有任何见解？

[  3%] Building CXX object CMakeFiles/ASPEN.dir/source/BatchFilter.cpp.o
In file included from /usr/local/include/boost/numeric/odeint.hpp:35:0,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/Filter.hpp:5,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/BatchFilter.hpp:5,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/source/BatchFilter.cpp:1:
/usr/local/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp: In instantiation of 'boost::numeric::odeint::controlled_step_result boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::try_step_v1(System, StateInOut&, boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::time_type&, boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::time_type&) [with System = System<arma::Col<double>::fixed<2>, arma::Mat<double> >; StateInOut = arma::Col<double>; ErrorStepper = boost::numeric::odeint::runge_kutta_cash_karp54<arma::Col<double> >; ErrorChecker = boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>; StepAdjuster = boost::numeric::odeint::default_step_adjuster<double, double>; Resizer = boost::numeric::odeint::initially_resizer; boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::time_type = double]':
/usr/local/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:283:27:   required from 'boost::numeric::odeint::controlled_step_result boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::try_step(System, StateInOut&, boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::time_type&, boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::time_type&) [with System = System<arma::Col<double>::fixed<2>, arma::Mat<double> >; StateInOut = arma::Col<double>; ErrorStepper = boost::numeric::odeint::runge_kutta_cash_karp54<arma::Col<double> >; ErrorChecker = boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>; StepAdjuster = boost::numeric::odeint::default_step_adjuster<double, double>; Resizer = boost::numeric::odeint::initially_resizer; boost::numeric::odeint::controlled_runge_kutta<ErrorStepper, ErrorChecker, StepAdjuster, Resizer, boost::numeric::odeint::explicit_error_stepper_tag>::time_type = double]'
/usr/local/include/boost/numeric/odeint/integrate/detail/integrate_times.hpp:101:81:   required from 'size_t boost::numeric::odeint::detail::integrate_times(Stepper, System, State&, TimeIterator, TimeIterator, Time, Observer, boost::numeric::odeint::controlled_stepper_tag) [with Stepper = boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<arma::Col<double> >, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_tag>; System = System<arma::Col<double>::fixed<2>, arma::Mat<double> >; State = arma::Col<double>; TimeIterator = __gnu_cxx::__normal_iterator<const double*, std::vector<double> >; Time = double; Observer = Observer::push_back_state<arma::Col<double> >; size_t = long unsigned int]'
/usr/local/include/boost/numeric/odeint/integrate/integrate_times.hpp:129:35:   required from 'size_t boost::numeric::odeint::integrate_times(Stepper, System, State&, TimeIterator, TimeIterator, Time, Observer) [with Stepper = boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<arma::Col<double> >, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_tag>; System = System<arma::Col<double>::fixed<2>, arma::Mat<double> >; State = arma::Col<double>; TimeIterator = __gnu_cxx::__normal_iterator<const double*, std::vector<double> >; Time = double; Observer = Observer::push_back_state<arma::Col<double> >; size_t = long unsigned int]'
/Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/source/BatchFilter.cpp:231:42:   required from void BatchFilter<state_type>::compute_reference_state_history(const std::vector<double>&, std::vector<_RealType>&, std::vector<arma::Mat<double> >&) [with state_type = arma::Col<double>::fixed<2>]'
/Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/source/BatchFilter.cpp:327:16:   required from here
/usr/local/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:481:12: error: no match for call to '(boost::numeric::odeint::unwrap_reference<System<arma::Col<double>::fixed<2>, arma::Mat<double> > >::type {aka System<arma::Col<double>::fixed<2>, arma::Mat<double> >}) (arma::Col<double>&, arma::Col<double>&, boost::numeric::odeint::controlled_runge_kutta<boost::numeric::odeint::runge_kutta_cash_karp54<arma::Col<double> >, boost::numeric::odeint::default_error_checker<double, boost::numeric::odeint::range_algebra, boost::numeric::odeint::default_operations>, boost::numeric::odeint::default_step_adjuster<double, double>, boost::numeric::odeint::initially_resizer, boost::numeric::odeint::explicit_error_stepper_tag>::time_type&)'
         sys( x , m_dxdt.m_v ,t );
         ~~~^~~~~~~~~~~~~~~~~~~~~
In file included from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/Filter.hpp:6:0,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/BatchFilter.hpp:5,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/source/BatchFilter.cpp:1:
/Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/System.hpp:32:7: note: candidate: void System<state_type, jacobian_type>::operator()(const state_type&, state_type&, double) [with state_type = arma::Col<double>::fixed<2>; jacobian_type = arma::Mat<double>] <near match>
  void operator() (const state_type & x , state_type & dxdt , const double t ){
       ^~~~~~~~
/Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/System.hpp:32:7: note:   conversion of argument 2 would be ill-formed:
In file included from /usr/local/include/boost/numeric/odeint.hpp:35:0,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/Filter.hpp:5,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/BatchFilter.hpp:5,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/source/BatchFilter.cpp:1:
/usr/local/include/boost/numeric/odeint/stepper/controlled_runge_kutta.hpp:481:25: error: cannot bind non-const lvalue reference of type 'arma::Col<double>::fixed<2>&' to an rvalue of type 'arma::Col<double>::fixed<2>'
         sys( x , m_dxdt.m_v ,t );
                  ~~~~~~~^~~
In file included from /usr/local/include/armadillo:568:0,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/include/BatchFilter.hpp:4,
                 from /Users/bbercovici/GDrive/CUBoulder/Research/code/ASPEN_gui_less/lib/source/BatchFilter.cpp:1:
/usr/local/include/armadillo_bits/Col_meat.hpp:1145:1: note:   after user-defined conversion: arma::Col<eT>::fixed<fixed_n_elem>::fixed(const arma::Base<eT, T1>&) [with T1 = arma::Mat<double>; long long unsigned int fixed_n_elem = 2; eT = double]
 Col<eT>::fixed<fixed_n_elem>::fixed(const Base<eT,T1>& A)
 ^~~~~~~
make[2]: *** [CMakeFiles/ASPEN.dir/source/BatchFilter.cpp.o] Error 1
make[1]: *** [CMakeFiles/ASPEN.dir/all] Error 2
make: *** [all] Error 2

作为参考，System.hpp在这里：

# System.hpp
template <typename state_type,typename jacobian_type = arma::mat> class System {
public:

    System(const Args & args,
        unsigned int N_est,
        state_type (*estimate_dynamics_fun)(double, const state_type & , const Args & args) ,
        jacobian_type (*jacobian_estimate_dynamics_fun)(double, const state_type & , const Args & args),
        unsigned int N_true = 0,
        state_type (*true_dynamics_fun)(double, const state_type & , const Args & args) = nullptr) 
    : N_est(N_est), N_true(N_true){

        this -> estimate_dynamics_fun = estimate_dynamics_fun;
        this -> true_dynamics_fun = estimate_dynamics_fun;
        this -> jacobian_estimate_dynamics_fun = jacobian_estimate_dynamics_fun;
        this -> args = args;

    }

    System(const Args & args,
        unsigned int N_true,
        state_type (*true_dynamics_fun)(double, const state_type & , const Args & args)) 
    : N_est(0), N_true(N_true){
        this -> true_dynamics_fun = true_dynamics_fun;
        this -> args = args;
    }

    void operator() (const state_type & x , state_type & dxdt , const double t ){


        if (this -> true_dynamics_fun != nullptr){

            dxdt.rows(this -> N_est + this -> N_est * this -> N_est,
                this -> N_est + this -> N_est * this -> N_est + this -> N_true - 1) = this -> true_dynamics_fun(t,

                x.rows(this -> N_est + this -> N_est * this -> N_est,
                    this -> N_est + this -> N_est * this -> N_est + this -> N_true - 1),args);


            }   

            if (this -> estimate_dynamics_fun != nullptr){

                arma::mat Phi = arma::reshape(x.rows(this -> N_est,
                    this -> N_est + this -> N_est * this -> N_est - 1), this -> N_est, this -> N_est );

                arma::mat A = this -> jacobian_estimate_dynamics_fun(t,x.rows(0,this -> N_est - 1),this -> args);

                dxdt.rows(0,this -> N_est - 1) = this -> estimate_dynamics_fun(t,x.rows(0,this -> N_est - 1),this -> args);
                dxdt.rows(this -> N_est,
                    this -> N_est + this -> N_est * this -> N_est - 1) = arma::vectorise(A * Phi);

            }



        }

    protected:
        const unsigned int N_est;
        const unsigned int N_true;
        state_type (*estimate_dynamics_fun)(double, const state_type &  , const Args & args) = nullptr;
        state_type (*true_dynamics_fun)(double, const state_type &  , const Args & args) = nullptr;
        jacobian_type (*jacobian_estimate_dynamics_fun)(double, const state_type &  , const Args & args) = nullptr;
        Args args;
    };
答案

我在这里找到了解决调整大小问题的方法：

https://reformatcode.com/code/c/armadillo-conflicts-with-boost-odeint-odeint-resizes-the-state-vector-to-zero-during-integration

解决方案包括使用以下适配器扩展boost：numeric :: odeint，这些适配器处理arma::vec的自动调整大小以便安全地使用odeint

#include <armadillo>
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.n_rows == y.n_rows;
    }
};

template<>
struct resize_impl<arma::vec, arma::vec>
{
    static void resize(arma::vec &v1, const arma::vec & v2)
    {
        v1.resize(v2.n_rows,1); 
    }
};

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

以上是关于使用Armadillo和boost :: numeric :: odeint进行模板实例化的主要内容，如果未能解决你的问题，请参考以下文章

Armadillo 与 Boost Odeint 冲突：Odeint 在积分期间将状态向量的大小调整为零

比较 blitz++、犰狳、boost::MultiArray

Armadillo + OpenBLAS 多线程

使用 Armadillo 和 OpenBLAS 进行多线程处理时性能不一致

C++ Armadillo 和 OpenMp：外积求和的并行化 - 定义 Armadillo 矩阵的约简

RcppArmadillo 和 Armadillo 之间的性能差异