[Rcpp-devel] Segfault error during simulation in Rcpp

[Dirk Eddelbuettel]

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")

*/



-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com