我对 C++ 很陌生，我正在尝试通过 Rcpp 使用它来加速我的 R 代码。
下面的代码从 t0 集成到 t1 - 这是在“lorenz”函数中完成的。Test4 使用“lorenz”“counts”次数进行积分。然而，在时间“t1”，系统状态在系统重新运行之前在“write_lorenz”中被修改，这就是问题所在。如果我通过从 R 调用 test4 一遍又一遍地运行相同的程序，打印到屏幕上总是会产生相同的结果，但是，我返回的矩阵“u”不会，并且似乎最终会收敛到任何“t1”是问题。
我的输入值没有改变，所以我想知道是否存在内存泄漏，或者是否发生了其他事情，如何解决它。另外我想知道我对“u”的初始化是否不正确，我应该使用“new”命令。
我尝试的是 NumericMatrix* u = NULL; *u = 新数值矩阵；然后我尝试以 *u(1,2) 的形式访问矩阵的元素，但是以这种方式访问元素会导致错误说 u 不是函数。
任何帮助将不胜感激
我从以下站点修改了此代码
http://headmyshoulder.github.io/odeint-v2/examples.html
所以我可以将它与 Rcpp 一起使用
//###################################R Code ###############################
library(Rcpp)
sourceCpp("test4.cpp")
sigma <- 10.0
R <-28.0
b <- 8.0 / 3.0
a2 <- c(10.0 , 1.0 , 1.0) #initial conditions X0,I0,Y0
L2 <- c(0.0 , 2.0 , 0.1) #initial time, kick time, error
counts <- 2
kick <-1.0; # kick size
pars <-c(sigma,R,b,kick)
test4(a=a,L2=L2,counts=counts,pars= pars)
// C ++ code
//[[Rcpp::depends(BH)]]
//[[Rcpp::depends(RcppEigen)]]
//[[Rcpp::plugins("cpp11")]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <math.h>
#include <boost/array.hpp>
#include <boost/numeric/odeint.hpp>
#include <boost/algorithm/string.hpp>
using namespace boost::numeric::odeint;
using namespace std;
using namespace Rcpp;
using namespace Eigen;
double sigma =0;
double e =0;
double b =0 ;
double t0 =0;
double t1 = 0;
double kick =0;
double X0 = 0;
double I0 =0;
double Y0 =0;
double N = 0;
int counter =0;
typedef boost::array< double , 3 > state_type;
NumericMatrix u(4,5);
void lorenz( const state_type &x , state_type &dxdt , double t )
{
dxdt[0] = sigma * ( x[1] - x[0] );
dxdt[1] = e * x[0] - x[1] - x[0] * x[2];
dxdt[2] = -b * x[2] + x[0] * x[1];
}
void write_lorenz( const state_type &x , const double t )
{
if(t==t1){
X0 = x[0];
I0 = x[1];
Y0 = x[2];
N = X0+I0+Y0;
X0 = X0 + exp(kick*N);
counter++;
//for some reason cout and u don't match for multiple runs of the
//program
cout << t << '\t' << X0 << '\t' << I0 << '\t' << Y0 << endl;
u(counter,0) = t;
u(counter,1) = X0;
u(counter,2) = I0;
u(counter,3) = Y0;
}
}
// [[Rcpp::export]]
NumericMatrix test4(NumericVector a, NumericVector L2,int counts,NumericVector pars)
{
double error; // control integration error
// initialize model parameters
//maybe these should be parenthesis?
sigma = pars[0]; //# average per capita birh rate 1-10(the mean is 4.7)
e = pars[1]; // # per capita average growth rate
b= pars[2];
kick = pars[3]; // kick size
//cout << sigma << R << b << kick << endl;
//myfile.open (ST);
X0 = a[0]; I0 =a[1]; Y0 = a[2];
int i = 0;
t0 = L2[0];
t1 = L2[1];
error = L2[2];
u(0,0) = t0;
u(0,1) = X0;
u(0,2) = I0;
u(0,3) = Y0;
// initial conditions
for(i=0;i<counts;i++)
{
state_type x = { X0 , I0 , Y0 };
integrate( lorenz , x , t0 , t1 , error , write_lorenz );
}
return u; // the variable I hope will be global
}