CODEWITHSUNDEEP
rrcpp
Daniel
Daniel

Reputation: 41

Inconsistent results in R using .Call inside lapply

I've created a package using the Rcpp package's Rcpp.package.skeleton function containing one .cpp file with a C++ function returning 0:

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
RcppExport SEXP just_zero() {
BEGIN_RCPP
    Rcpp::RNGScope __rngScope;
    return wrap(0.0);
END_RCPP
}

When I install and load the package from R I can call the function with .Call via lapply. As expected, it (appears to) always return 0:

> x <- lapply(seq(10000), function(i) { x <- .Call('just_zero'); stopifnot(x == 0); x } )
#*no errors!*

However, apparently the value which is returned by lapply contains non-zeros:

> range(simplify2array(x))
[1] 0 3

Unfortunately, using set.seed doesn't make these returned values reproducible, and sometimes I do get [1] 0 0, sometimes other values, e.g. [1] "0" "TRUE". An additional clue is that removing the line Rcpp::RNGScope __rngScope; fixes the problem.

Why are there non-zero elements in the object returned by lapply (particularly as we've checked the values returned by .Call), and how come using RNGScope is causing it?

I've reproduced this behaviour on Linux and OS X. Session info from OS X pasted below:

> sessionInfo()
R Under development (unstable) (2016-08-03 r71023)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X Mavericks 10.9.5

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] bug_1.0         devtools_1.12.0

loaded via a namespace (and not attached):
[1] tools_3.4.0   withr_1.0.2   memoise_1.0.0 Rcpp_0.12.6   git2r_0.15.0
[6] digest_0.6.10

Upvotes: 3

Views: 233

Answers (2)

Kevin Ushey
Kevin Ushey

Reputation: 21315

Here's an example that reproduces the issue consistently for me.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
void dummy() {}

extern "C" SEXP just_zero() {
  Rcpp::RNGScope rngScope;
  return wrap(0.0);
}

/*** R
n <- 1E5
x <- unlist(lapply(seq(n), function(i) {
  .Call('just_zero')
}))

unique(x)
*/

Calling Rcpp::sourceCpp() on this gives me e.g.

> unique(x)
 [1]     0  8371 16021 20573 25109 43103 47563 56438 60852 78413 82773 95765

What's causing this issue? To understand this, we need to understand the definition of RNGScope: we can see it here:

https://github.com/RcppCore/Rcpp/blob/9fbdf78fe225607dc2c8af7267a6840af0aba10e/inst/include/Rcpp/stats/random/random.h#L27-L31

and the methods it uses are defined here:

https://github.com/RcppCore/Rcpp/blob/9fbdf78fe225607dc2c8af7267a6840af0aba10e/src/api.cpp#L63-L75

with the most important function, PutRNGState, being defined here:

https://github.com/wch/r-source/blob/1c88a057594a0348f2bf75514a8015caeedbff93/src/main/RNG.c#L424-L446

Now, this is effectively what happens when just_zero is called:

  1. RNGScope object is created, with its associated destructor primed to run.
  2. The result of wrap(0.0), being a REALSXP with length 1 and value 0, is created. Note that this object is unprotected and so is eligible for garbage collection.
  3. The function returns, and the RNGScope destructor is called.
  4. This calls the R routine PutRNGstate(), which itself will call allocVector, thereby triggering the garbage collector. This implies the SEXP object you wanted to return can be collected, and hence will be garbage!

So, in sum -- use Rcpp attributes, as it will take care of doing this all safely for you.

To understand why Rcpp attributes makes this 'safe', look at the generated code:

// just_zero
SEXP just_zero();
RcppExport SEXP sourceCpp_0_just_zero() {
BEGIN_RCPP
    Rcpp::RObject __result;
    Rcpp::RNGScope __rngScope;
    __result = Rcpp::wrap(just_zero());
    return __result;
END_RCPP
}

Note that the output result is assigned to Rcpp::RObject, which protects the object, and we ensure that this object is constructed before the RNGScope object, which ensures it will remain protected while the RNGScope destructor is run.

Upvotes: 5

coatless
coatless

Reputation: 20746

This is not fully an answer, but I'll expand it as it goes. However, I'm not able to reproduce this behavior on macOS.

I believe the issue you are encountering is due to your use of the RcppExport alongside Rcpp Attributes, e.g. :

// [[Rcpp::export]]
RcppExport SEXP just_zero()

Should be done via:

library("Rcpp")
cppFunction("double just_zero() {
    return 0.0;
}")

x <- lapply(seq(10000), function(i) { x <- just_zero(); stopifnot(x == 0); x } )

all(range(simplify2array(x)) == 0)

or place the following into file.cpp:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
double just_zero() {
    return 0.0;
}

/*** R
x <- lapply(seq(10000), function(i) { x <- just_zero(); stopifnot(x == 0); x } )
all(range(simplify2array(x)) == 0)
*/

If the intent is to use sourceCpp().

Upvotes: 3

Related Questions