【问题标题】:Rcpp NumericalMatrix data type, global declarationRcpp NumericalMatrix 数据类型,全局声明
【发布时间】:2015-05-17 06:55:39
【问题描述】:

我是 C++ 的新手,我正在尝试通过 Rcpp 使用它来加速我的 R 代码。

下面的代码从 t0 集成到 t1——这是在“lorenz”函数中完成的。 Test4 使用“lorenz”“counts”次数进行积分。然而,在时间“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 



}

【问题讨论】:

  • What I tried was NumericMatrix* u = NULL; *u = new NumericMatrix; 不要。您的 C++ 代码中没有内存泄漏,除非 NumericMatrix 是一个编码不佳的类。我不知道你的 C++ 代码上面是什么东西。唯一需要注意的是返回值——您正在按值返回 NumericMatrix,这在 C++ 中是可以的,但不知道您的其他代码期望什么。
  • 这里是包的共同作者。它不是一个编码不良的类——这些是底层 R SEXP 类型的代理对象——带有内部指针。 Rcpp Gallery 有 90 多个工作示例,Rcpp package 中有 8 个 pdf 小插图,my website under paper, talks, blog, ... 有更多示例,然后是 Rcpp book。最后,按值复制在这里没有问题,因为一切都与SEXP 指针对象相关。

标签: c++ r boost rcpp odeint


【解决方案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> 

【讨论】:

    猜你喜欢
    • 2019-04-02
    • 2020-07-25
    • 2019-01-04
    • 2012-04-30
    • 2020-05-15
    • 2021-11-22
    • 1970-01-01
    • 2015-07-25
    • 2020-01-15
    相关资源
    最近更新 更多