<div dir="ltr">Hi All, <div><br></div><div style> I am trying to simulate from the following model:</div><div style><br></div><div style>Y_t = Pois( phi * N_t ); </div><div style>N_t = r * N_{t-1} * exp( -N_{t-1}^theta + e_t )</div>
<div style>e_t ~ Norm(0, sigma)</div><div style><br></div><div style>so I have written a Rcpp function with prototype:</div><div style><br></div><div style>genRickerCpp(int days, int nSimul, int nBurn, NumericMatrix params)<br>
</div><div style><br></div><div style>where:</div><div style><br></div><div style>days = length of each simulated path</div><div style>nSimul = number of simulated paths</div><div style>nBurn = number of simulations I'm going to discard before storing </div>
<div style> each path</div><div style>params = matrix of input parameters that can be either 1 by 4 (in which</div><div style> case all the paths are simulated using the same parameters)</div><div style>
or nSimul by 4 (in which case each path uses a different </div><div style> vector of parameters)</div><div style><br></div><div style>This is my code:</div><div style><br></div><div style><div>#include <Rcpp.h></div>
<div><br></div><div>using namespace Rcpp;</div><div><br></div><div>// [[Rcpp::export]]</div><div>NumericMatrix genRickerCpp(int days, int nSimul, int nBurn, NumericMatrix params)</div><div>{</div><div> int nParams = params.ncol(); </div>
<div> int totDays = nBurn + days;</div><div> bool multiParams = false;</div><div> </div><div> if(nParams != 4) stop("Wrong number of parameters");</div><div> if(params.nrow() > 1) { multiParams = true; }</div>
<div> if(multiParams == true && params.nrow() != nSimul) </div><div> stop("Number of parameters vectors is different from the number of simulations");</div><div> </div><div> double r = exp(params(0, 0));</div>
<div> double theta = exp(params(0, 1));</div><div> double sigma = exp(params(0, 2));</div><div> double phi = exp(params(0, 3));</div><div> </div><div> NumericVector procNoise( rnorm( totDays * nSimul ) );</div><div>
NumericVector initState( runif( nSimul ) ); </div><div> NumericMatrix output( nSimul, days );</div><div> </div><div> NumericVector::iterator noiseIter = procNoise.begin();</div><div> NumericVector::iterator initIter = initState.begin();</div>
<div> </div><div> double currState;</div><div> </div><div> for(int iRow = 0; iRow < nSimul; iRow++, initIter++)</div><div> {</div><div> </div><div> if( multiParams == true )</div><div> {</div><div> r = exp(params(iRow, 0));</div>
<div> theta = exp(params(iRow, 1));</div><div> sigma = exp(params(iRow, 2));</div><div> phi = exp(params(iRow, 3));</div><div> }</div><div> </div><div> currState = *initIter;</div><div> </div><div> for(int iCol = 1; iCol <= nBurn; iCol++, noiseIter++){</div>
<div> currState = r * currState * exp( - pow( currState, theta ) + *noiseIter * sigma );</div><div> }</div><div> </div><div> output(iRow, 0) = rpois(1, phi * currState)[0];</div><div> </div><div> for(int iCol = 1; iCol < days; iCol++, noiseIter++){</div>
<div> currState = r * currState * exp( - pow( currState, theta ) + *noiseIter * sigma );</div><div> output(iRow, iCol) = rpois(1, phi * currState)[0];</div><div> }</div><div> </div><div> }</div><div> </div>
<div> return output;</div><div> </div><div>}</div><div><br></div><div><br></div><div style>the function seems to work well, I tried to compare the output the an equivalent R function</div><div style>and I get the same results. The problem is that if I run it a lot of times:</div>
<div style><br></div><div style><div>library(Rcpp)</div><div>sourceCpp("~/Desktop/genRickerCpp.cpp")</div><div><br></div><div>for(ii in 1:10^6){</div><div>data <- genRickerCpp(days = 1, nSimul = 1, nBurn = 1,</div>
<div> params = matrix(log(c(r = exp(3.8), theta = 1, sigma = 0.3, phi = 10)), 1, 4))</div><div>data <- as.numeric(data)</div><div>}</div><div><br></div><div style>occasionally R crashes with error:</div>
<div style><br></div><div style><div> *** caught segfault ***</div><div>address 0x28, cause 'memory not mapped'</div><div><br></div><div style>the strange thing is that in most cases I can call it 10^6 times without any error.</div>
<div style>I tried to go through the code in gdb, but I didn't see anything wrong. I also ran the previous</div><div style>R code in valgrind, and there I get the following errors while the code is running:</div><div style>
<br></div><div style><div>==3031== Invalid read of size 1</div><div>==3031== at 0x4EF7BF1: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)</div><div>==3031== Address 0xebd20c8 is 1,688 bytes inside a block of size 1,968 free'd</div><div>==3031== at 0x4C2A82E: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)</div>
<div>==3031== by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0xDD2017D: sourceCpp_74073_genRickerCpp (random.h:71)</div><div>==3031== by 0x4EEBEC5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25BCC: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>
==3031== by 0x4F28E0C: Rf_applyClosure (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2588F: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== </div><div>==3031== Invalid read of size 8</div><div>==3031== at 0x4EF7F89: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)</div><div>==3031== Address 0xebd20e8 is 1,720 bytes inside a block of size 1,968 free'd</div><div>==3031== at 0x4C2A82E: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)</div>
<div>==3031== by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0xDD2017D: sourceCpp_74073_genRickerCpp (random.h:71)</div><div>==3031== by 0x4EEBEC5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25BCC: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>
==3031== by 0x4F28E0C: Rf_applyClosure (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2588F: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== </div><div>==3031== Invalid read of size 8</div><div>==3031== at 0x4EF7FA8: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)</div><div>==3031== Address 0xebd20f8 is 1,736 bytes inside a block of size 1,968 free'd</div><div>==3031== at 0x4C2A82E: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)</div>
<div>==3031== by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0xDD2017D: sourceCpp_74073_genRickerCpp (random.h:71)</div><div>==3031== by 0x4EEBEC5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25BCC: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>
==3031== by 0x4F28E0C: Rf_applyClosure (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2588F: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== </div><div>==3031== Invalid read of size 8</div><div>==3031== at 0x4EF7FC6: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)</div><div>==3031== Address 0xebd20d0 is 1,696 bytes inside a block of size 1,968 free'd</div><div>==3031== at 0x4C2A82E: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)</div>
<div>==3031== by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0xDD2017D: sourceCpp_74073_genRickerCpp (random.h:71)</div><div>==3031== by 0x4EEBEC5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25BCC: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>
==3031== by 0x4F28E0C: Rf_applyClosure (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2588F: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== </div><div>==3031== Invalid read of size 8</div><div>==3031== at 0x4EF7F7C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)</div><div>==3031== Address 0xebd20f0 is 1,728 bytes inside a block of size 1,968 free'd</div><div>==3031== at 0x4C2A82E: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)</div>
<div>==3031== by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0xDD2017D: sourceCpp_74073_genRickerCpp (random.h:71)</div><div>==3031== by 0x4EEBEC5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25BCC: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>
==3031== by 0x4F28E0C: Rf_applyClosure (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2588F: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== </div><div>==3031== Invalid read of size 1</div><div>==3031== at 0x4EF7E34: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)</div><div>==3031== Address 0xebd20c8 is 1,688 bytes inside a block of size 1,968 free'd</div><div>==3031== at 0x4C2A82E: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)</div>
<div>==3031== by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0xDD2017D: sourceCpp_74073_genRickerCpp (random.h:71)</div><div>==3031== by 0x4EEBEC5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25BCC: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>
==3031== by 0x4F28E0C: Rf_applyClosure (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2588F: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== </div><div>==3031== Invalid read of size 1</div><div>==3031== at 0x4EF7C27: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)</div><div>==3031== Address 0xebd20c9 is 1,689 bytes inside a block of size 1,968 free'd</div><div>==3031== at 0x4C2A82E: free (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)</div>
<div>==3031== by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== by 0xDD2017D: sourceCpp_74073_genRickerCpp (random.h:71)</div><div>==3031== by 0x4EEBEC5: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F25BCC: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>
==3031== by 0x4F28E0C: Rf_applyClosure (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2588F: Rf_eval (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)</div><div>==3031== by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)</div>
<div>==3031== </div><div><br></div><div><br></div><div style>The only thing that I understand is that probably there is an invalid read somewhere,</div><div style>but I went through the code several times and I even re-wrote everything from </div>
<div style>scratch and I can't find anything wrong. Given that I'm a C++/Rcpp beginner </div><div style>I guess that probably I'm doing a systematic error in my code, maybe you can point it</div><div style>out to me.</div>
<div style>Finally I'm running my code on Ubuntu 12.04 (precise) 64-bit and R version 2.15.3.<br>Thanks a lot.</div><div style><br></div><div style>Matteo</div></div></div></div></div></div>