如何在 R 中使用带有 Rcpp 的 C++ ODE 求解器? [英] How to use C++ ODE solver with Rcpp in R?
问题描述
为了评估 R 和 C++ 之间求解 ODE 的速度差异，我在 R 中构建了以下 ODE 系统:
To assess the difference in speed to solve ODEs between R and C++, I build the following ODE system in 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 包.
To solve it, I used the deSolve package.
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))
这有效.我试图通过在 R 中使用 Rcpp 库，使用 C++ 求解器解决相同的系统.这是我添加的内容:
This works. I tried to solve the same system with C++ solver, by using Rcpp library in R. Here is what I add:
#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_control > >::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 );
In function 'bool my_fun10_solver()': ex3.cpp:114:64: error: no matching function for call to '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 );
我的代码有什么问题?
以下是我采用并适应我的问题的脚本.此脚本运行良好.
Below is the script I took and adapted to my problem. This scripts works well.
#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;
}
推荐答案
您可以尝试以下操作吗?
Could you try the following please?
integrate_adaptive(make_controlled( 1E-12 , 1E-12 , stepper_type () ) ,
eqsir , x , 1.0 , 2.0 , 1.0/52.0 , write_cout_2 );
不过有一些优化建议.您正在求解 ODE，这表明您是物理学家或工程师，这太棒了，我也是物理学家，而物理学家只是摇滚.
A couple of optimisation suggestions though. You are solving an ODE, which suggests that you are a physicist or engineer, which is awesome, I am a physicist too, and physicists just rock.
但计算机科学家也是如此，他们已尽最大努力尽可能快地完成工作.他们构建了编译器，这消除了很多相关的思考.尽你所能帮助他们.
But so do computer scientists, and they have tried their best to do thing as fast as possible. They build compilers, which take away a lot of the thinking involved. Try to help them where ever you can.
我为什么要告诉你这个?回到 ODE 的事情.boost 的自适应积分器可能会调用 eqsir() 1e9 次.尝试使该功能尽可能高效.考虑重写 my_fun2 以便 x 被覆盖而不是创建一个新的并返回它.x 是最后一个状态.除非您想绘制动态图，否则您不会太在意它.
Why am I telling you this? Back to the ODE thing. boost's adaptive intergrator is maybe calling eqsir() 1e9 times. Try to make that function as performant as possible. Consider rewriting my_fun2 so that x gets overwritten rather than create a new one and return it. x is the last state. You wouldn't care less about it unless you wanted to plot the dynamics.
void my_fun2(Rcpp::NumericVector &x, const double t);
还有
Rcpp::NumericVector nvec=boost_array_to_nvec(x);
Rcpp::NumericVector nvec2(3);
必须在每次调用时分配新内存.最后，考虑使用 nvec 和 state_type 的 2 个转换器以及我作为第二个选项编写的参考资料.您的新 eqsir 看起来像这样并且运行速度可能会快很多.
have to allocate new memory on every call. Finally, consider using the 2 converters for nvec and state_type with the references I wrote as second option. Your new eqsir would look like this and run possibly a lot faster.
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);
}
这篇关于如何在 R 中使用带有 Rcpp 的 C++ ODE 求解器?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
