私は C++ にかなり慣れていないので、Rcpp 経由で使用して R コードを高速化しようとしています。
以下のコードは、t0 から t1 までを統合します。これは「lorenz」関数で行われます。Test4 は、「lorenz」「counts」回数を使用して統合します。ただし、時間「t1」で、システムが再実行される前にシステムの状態が「write_lorenz」で変更され、これが問題の場所です。Rからtest4を呼び出して同じプログラムを何度も実行すると、画面への出力は常に同じ結果になりますが、返された行列「u」はそうではなく、最終的に「t1」が何であるかに収束するようです問題。
入力値が変わらないので、メモリリークがあるのか、それとも何か他のことが起こっているのか、どうすれば修正できるのか疑問に思っています。また、「u」の初期化が正しくないので、「new」コマンドを使用する必要があるのではないかと考えています。
私が試したのは NumericMatrix* u = NULL; です。*u = 新しい NumericMatrix; 次に、たとえば *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
}