Rcpp NumericalMatrix数据类型,全局声明

Question

我对C++很新,我试图通过Rcpp来加速我的R代码.

下面的代码从t0到t1集成 - 这是在"lorenz"函数中完成的.Test4使用"lorenz""计数"次数进行整合.但是,在时间"t1",系统的状态在系统重新运行之前在"write_lorenz"中被修改,这就是问题所在.如果我通过从R调用test4反复运行相同的程序,打印到屏幕总是产生相同的结果,但是,我返回的矩阵"u"没有,并且似乎最终收敛到任何"t1"是什么问题.

我的输入值没有改变所以我想知道是否有内存泄漏,或者是否有其他事情发生,如何解决它.另外我想知道我的"u"初始化是否不正确我应该使用"new"命令.

我试过的是NumericMatrix*u = NULL;*u =新的NumericMatrix; 然后我尝试访问矩阵的元素为*u(1,2),但是以这种方式访问​​元素会导致错误,说你不是函数.

任何帮助将不胜感激

我从以下网站修改了此代码

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++文件的简单修改.我们简单地定义了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>