Rcpp Тип данных NumericalMatrix, глобальная декларация

Question

Я довольно новичок в C++, и я пытаюсь использовать его через Rcpp для ускорения моего R-кода.

Подключенный код от t0 до t1- это делается в функции «lorenz». Test4 объединяет, используя «lorenz» «counts» количество раз. Однако во время «t1» состояние системы изменяется в «write_lorenz» до повторной перезагрузки системы, и это проблема. Если я запускаю одну и ту же программу снова и снова, вызывая test4 из R, печать на экране всегда дает тот же результат, однако моя возвращенная матрица «u» не делает и, похоже, в конечном итоге сходится к тому, что «t1» есть проблема.

Мои значения ввода не меняются, поэтому мне интересно, есть ли утечка памяти или что-то еще происходит, как ее исправить. Также мне интересно, неверна ли моя инициализация «u», и я должен использовать «новую» команду.

Что я пытался: NumericMatrix * u = NULL; * u = new Числовая матрица; , а затем я попытался получить доступ к элементам матрицы как * 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++, с которым вы связались. Мы просто определяем 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 (intentially ограничено кулак несколько строк):

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>