Как использовать решатель ODE C++ с Rcpp в R? ⇐ C++

[Anonymous]

Чтобы оценить разницу в скорости решения ОДУ между R и C++, я создаю следующую систему ОДУ на 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
#include

// include Boost's odeint
#include

// [[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

Подробнее здесь: https://stackoverflow.com/questions/500 ... -rcpp-in-r