Как использовать C ++ ODE Solver с Rcpp в R?

Question

Чтобы оценить разницу в скорости для решения ODE между R и C ++, я создаю следующую систему ODE в R:

modelsir_cpp =function(t,x){
S = x[1]
I1 = x[2]
I2 = x[3]
N=S+I1+I2
with(as.list(parm), {
   dS=B*I1-mu*S-beta*(S*(I1+I2)/N)
   dI1=beta*(S*(I1+I2)/N)-B*I1-lambda12*I1
   dI2=lambda12*I1
   res=c(dS,dI1,dI2)
   return(res)
 })
}

Чтобы решить эту проблему, я использовал пакет deSolve.

times = seq(0, 10, by = 1/52)
parm=c(B=0.01,mu=0.008,beta=10,lambda12=1)
xstart=c(S=900,I1=100,I2=0)
out = as.data.frame(lsoda(xstart, times, modelsir, parm))

Это работает. Я попытался решить ту же систему с помощью решателя C ++, используя библиотеку Rcpp в R. Вот что я добавляю:

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;

typedef boost::array< double ,3 > state_type;

// [[Rcpp::export]]
 Rcpp::NumericVector my_fun2(const Rcpp::NumericVector &x, const double t){
 Function f("modelsir_cpp");
    return f(_["t"]=t,_["x"]=x);
}
    void eqsir(const state_type &x, state_type &dxdt, const double t){
      Rcpp::NumericVector nvec=boost_array_to_nvec(x);
      Rcpp::NumericVector nvec2(3);
      nvec2=my_fun2(nvec,t);
      dxdt=nvec_to_boost_array(nvec2);
        }


void write_cout_2( const state_type &x , const double t ) {
  // use Rcpp's stream
  Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
}
typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
bool my_fun10_solver() {
  state_type x = { 900 , 100, 0 }; // initial conditions
  integrate_adaptive(make_controlled( 1E-12 , 1E-12 , stepper_type () ) ,
                     eqsir , x , 1.0 , 2 , 1/52 , write_cout_2 );
  return true;
}

Но появилось сообщение об ошибке:

В функции 'bool my_fun10_solver ()': ex3.cpp: 114: 64: ошибка: нет соответствующей функции для вызова 'integrate_adaptive (boost :: numeric :: odeint :: result_of :: make_controlled>> :: type, void (&) (const state_type &, state_type &, double) , state_type &, double, int, int, void (&) (const state_type &, double)) ' eqsir, x, 1,0, 2, 1/52, write_cout_2);

Что не так с моим кодом?

Ниже приведен сценарий, который я взял и адаптировал к своей проблеме. Этот скрипт хорошо работает.

---

#include <Rcpp.h>
#include <boost/array.hpp>

// include Boost's odeint
#include <boost/numeric/odeint.hpp>

// [[Rcpp::depends(BH)]]

using namespace Rcpp;
using namespace std;
using namespace boost::numeric::odeint;


typedef boost::array< double ,3 > state_type;

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

  dxdt[0] = 3.0/(2.0*t*t) + x[0]/(2.0*t);
  dxdt[1] = 3.0/(2.0*t*t) + x[1]/(2.0*t);
  dxdt[2] = 3.0/(2.0*t*t) + x[1]/(2.0*t);
}

void write_cout( const state_type &x , const double t ) {
  // use Rcpp's stream
  Rcpp::Rcout << t << '\t' << x[0] << '\t' << x[1] << '\t' << x[2] <<  endl;
}

 typedef runge_kutta_dopri5< state_type > stepper_type;

// [[Rcpp::export]]
bool boostExample() {
  state_type x = { 1.0 , 1.0, 1.0 }; // initial conditions
  integrate_adaptive(make_controlled( 1E-12 , 1E-12 , stepper_type () ) ,
                     rhs , x , 1.0 , 10.0 , 0.1 , write_cout );
  return true;
}

Answer 1

Не могли бы вы попробовать следующее, пожалуйста?

integrate_adaptive(make_controlled( 1E-12 , 1E-12 , stepper_type () ) ,
                 eqsir , x , 1.0 , 2.0 , 1.0/52.0 , write_cout_2 );

Пара предложений по оптимизации. Вы решаете ODE, что говорит о том, что вы физик или инженер, и это здорово, я тоже физик, а физики просто потрясающие.

Но так же делают и ученые-компьютерщики, и они изо всех сил старались сделать все как можно быстрее. Они создают компиляторы, которые отнимают уйму мышления. Постарайся помочь им, где только сможешь.

Почему я говорю вам это? Вернемся к ODE. Адаптивный интегратор Boost может вызывать eqsir() 1e9 раз. Постарайтесь сделать эту функцию максимально эффективной. Попробуйте переписать my_fun2, чтобы x перезаписывался, а не создавал новый и возвращал его. x является последним состоянием. Тебя это не волнует, если ты не хочешь построить динамику.

void my_fun2(Rcpp::NumericVector &x, const double t);

Также

Rcpp::NumericVector nvec=boost_array_to_nvec(x);
Rcpp::NumericVector nvec2(3);

должен выделять новую память при каждом вызове. Наконец, рассмотрите возможность использования 2 конвертеров для nvec и state_type со ссылками, которые я написал как второй вариант. Ваш новый eqsir будет выглядеть так и работать, возможно, намного быстрее.

Rcpp::NumericVector nvec(3); // declared outside
void eqsir(const state_type &x, state_type &dxdt, const double t){
  boost_array_to_nvec(x, nvec);
  my_fun2(nvec,t);
  nvec_to_boost_array(nvec, dxdt);
}