Memory leaks in a simple Rcpp function

Question

I am developing a package in R that I would like to convert to Rcpp for better performance. I'm new to Rcpp (and C++ in general.) My problem is that the Rcpp function I've written works fine if I run it many times with one set of arguments, but if I try to loop it over many combinations of arguments, it springs memory leaks and causes the R session to abort.

Here is the code in R, which holds up well to any test I throw at it:

raw_noise <- function(timesteps, mu, sigma, phi) {
   delta <- mu * (1 - phi)
   variance <- sigma^2 * (1 - phi^2)
   noise <- vector(mode = "double", length = timesteps)
   noise[1] <- c(rnorm(1, mu, sigma))
   for (i in (1:(timesteps - 1))) {
     noise[i + 1] <- delta + phi * noise[i] + rnorm(1, 0, sqrt(variance))
  }
   return(noise)
}

Here is the code in Rcpp, using three Rcpp sugar functions (pow, sqrt, rnorm):

NumericVector raw_noise(int timesteps, double mu, double sigma, double phi) {
  double delta = mu * (1 - phi);
  double variance = pow(sigma, 2.0) * (1 - pow(phi, 2.0));
  NumericVector noise(timesteps);
  noise[0] = R::rnorm(mu, sigma);
  for(int i = 0; i < timesteps; ++i) {
    noise[i+1] = delta + phi*noise[i] + R::rnorm(0, sqrt(variance));
  }
  return noise;
}

What really confuses me is that this code runs without problems:

library(purrr)
rerun(10000, raw_noise(timesteps = 30, mu = 0.5, sigma = 0.2, phi = 0.3))

But when I run this code:

test_loop <- function(timesteps, mu, sigma, phi, replicates) {
  params <- cross_df(list(timesteps = timesteps, phi = phi, mu = mu, sigma = 
sigma))
  for (i in 1:nrow(params)) {
    print(params[i,])
    pmap(params[i,], raw_noise)
  }
}
library(purrr)
test_loop(timesteps=c(5, 6, 7, 8, 9, 10), mu=c(0.2, 0.5), sigma=c(0.2, 0.5),
      phi=c(0, 0.1))

More often than not, the R session aborts and RStudio crashes altogether. But sometimes I manage to catch this error message before the R session aborts:

Error in match(x, table, nomatch = 0L) : GC encountered a node (0x10db7af50) with an unknown SEXP type: NEWSXP at memory.c:1692

As I understand it, NEWSXP is an exotic object type in R that doesn't come up very often. What's happening looks to me like a memory leak, but I'm not at all sure how to fix it. Like I said, I'm new to Rcpp and C++ generally so I'd appreciate any nudges in the right direction.

Answer 1

You have an out of bounds error:

for(int i = 0; i < timesteps; ++i)

causes

noise[i+1]

to exceed the defined range since C++ indices start at 0 and not 1.

---

For example, 0 to timesteps - 1 has a length of timesteps and, thus, is okay.

but

0 to timesteps would have a length of timesteps + 1

This can be seen if you change noise[i+1] to noise(i+1), which performs a bounds check on the requested index.

Error in raw_noise(100, 2, 3, 0.2) : 
  Index out of bounds: [index=100; extent=100].

---

To address this, make the following change:

NumericVector raw_noise(int timesteps, double mu, double sigma, double phi) {
  double delta = mu * (1 - phi);
  double variance = pow(sigma, 2.0) * (1 - pow(phi, 2.0));
  NumericVector noise(timesteps);
  noise[0] = R::rnorm(mu, sigma);

  // change here
  for(int i = 0; i < timesteps - 1; ++i) { // 1 less time step

    noise[i+1] = delta + phi*noise[i] + R::rnorm(0, sqrt(variance));
  }
  return noise;
}