c++ - Rcpp NumericalMatrix 数据类型，全局声明

Question

我是 C++ 的新手，我正尝试通过 Rcpp 使用它来 boost 我的 R 代码。

下面的代码从 t0 整合到 t1 - 这是在“lorenz”函数中完成的。 Test4 使用“lorenz”“计数”次数进行积分。然而，在时间“t1”，系统状态在系统重新运行之前在“write_lorenz”中被修改，这就是问题所在。如果我通过从 R 调用 test4 一遍又一遍地运行相同的程序，打印到屏幕总是会产生相同的结果，但是，我返回的矩阵“u”不会，并且似乎最终会收敛到任何“t1”问题。

我的输入值没有改变，所以我想知道是否存在内存泄漏，或者是否发生了其他问题，如何修复。另外我想知道我对“u”的初始化是否不正确，我应该使用“new”命令。

我试过的是 数值矩阵* 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 



}

Answer 1

这是您链接到的纯 C++ 文件的简单改编。我们只需定义一个包含三个 vector 的 struct，我们将每次迭代的元素推送到其中——而不是将它们打印到标准输出。

对于增长的数据结构，我们更喜欢 C++ 标准库类型(因为我们的类型必须像 R 类型，它们的内部结构不能很好地匹配一个一个地增加，这对它们来说是昂贵的)。所以最后，我们只是复制到一个 R data.frame。

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

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

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

const double sigma = 10.0, r = 28.0, b = 8.0 / 3.0;
typedef boost::array< double , 3 > state_type;

void lorenz( const state_type &x , state_type &dxdt , double t ) {
    dxdt[0] = sigma * ( x[1] - x[0] );
    dxdt[1] = r * x[0] - x[1] - x[0] * x[2];
    dxdt[2] = -b * x[2] + x[0] * x[1];
}

struct foo { std::vector<double> a, b, c; };
struct foo f;

void append_lorenz(const state_type &x , const double t ) {
    f.a.push_back(x[0]);
    f.b.push_back(x[1]);
    f.c.push_back(x[2]);
}

using namespace Rcpp;

// [[Rcpp::export]]
DataFrame callMain() {
    state_type x = { 10.0 , 1.0 , 1.0 }; // initial conditions
    integrate( lorenz , x , 0.0 , 25.0 , 0.1 , append_lorenz );
    return DataFrame::create(Named("a") = f.a,
                             Named("b") = f.b,
                             Named("c") = f.c);
}

/*** R
res <- callMain()
print(head(res))
*/

这里是 R 代码的输出(有意限制为几行):

R> sourceCpp("/tmp/lorenz.cpp")

R> res <- callMain()

R> print(head(res))
         a        b       c
1 10.00000  1.00000 1.00000
2  9.40816  2.99719 1.12779
3  8.92164  5.35684 1.46991
4  8.68193  7.82671 2.05762
5  8.73730 10.42718 2.94783
6  9.11080 13.10452 4.18849
R>