Matteo,
This may reveal a subtle bug somewhere in the handling of temporary
expressions, with possible interactions with the garbage collection by R
itself.
That said, a trivial rewrite (below) which via a single sourceCpp() compiles
and runs, or re-runs it, has not lead to a single segfault here, and I even
increased N from 10^6 to 10^7.
The only other change I made was to move your params() matrix (you could have
used a named vector too; also why the log() with subsequent exp() ?) outside
the loop. It just works here...
Dirk
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix genRickerCpp(int days, int nSimul, int nBurn, NumericMatrix params) {
int nParams = params.ncol();
int totDays = nBurn + days;
bool multiParams = false;
if(nParams != 4) stop("Wrong number of parameters");
if(params.nrow() > 1) { multiParams = true; }
if(multiParams == true && params.nrow() != nSimul)
stop("Number of parameters vectors is different from the number of simulations");
double r = exp(params(0, 0));
double theta = exp(params(0, 1));
double sigma = exp(params(0, 2));
double phi = exp(params(0, 3));
NumericVector procNoise( rnorm( totDays * nSimul ) );
NumericVector initState( runif( nSimul ) );
NumericMatrix output( nSimul, days );
NumericVector::iterator noiseIter = procNoise.begin();
NumericVector::iterator initIter = initState.begin();
double currState;
for(int iRow = 0; iRow < nSimul; iRow++, initIter++) {
if( multiParams == true ) {
r = exp(params(iRow, 0));
theta = exp(params(iRow, 1));
sigma = exp(params(iRow, 2));
phi = exp(params(iRow, 3));
}
currState = *initIter;
for(int iCol = 1; iCol <= nBurn; iCol++, noiseIter++){
currState = r * currState * exp( - pow( currState, theta ) + *noiseIter * sigma );
}
output(iRow, 0) = rpois(1, phi * currState)[0];
for(int iCol = 1; iCol < days; iCol++, noiseIter++){
currState = r * currState * exp( - pow( currState, theta ) + *noiseIter * sigma );
output(iRow, iCol) = rpois(1, phi * currState)[0];
}
}
return output;
}
// #the function seems to work well, I tried to compare the output the an
// #equivalent R function
// #and I get the same results. The problem is that if I run it a lot of times:
// #library(Rcpp)
// #sourceCpp("~/Desktop/genRickerCpp.cpp")
/*** R
params <- matrix(log(c(r = exp(3.8), theta = 1, sigma = 0.3, phi = 10)), 1, 4)
for (ii in 1:10^7) {
data <- genRickerCpp(days = 1, nSimul = 1, nBurn = 1, params)
}
cat("Done\n")
*/