Issue with random numbers and R calls from C++

[Enchufa2]

It may have something to do with #465, which was unresolved back in 2016:

Rcpp::cppFunction("bool comp(Function fun) { return runif(1)[0] > as<double>(fun()); }")
set.seed(1)
sum(replicate(1e3, comp(function() runif(1, 0, 1))))
# 0
sum(replicate(1e3, comp(function() runif(1, 0.0001, 0.9999))))
# 481

This is very odd. I'm not sure what's happening here. I've spent hours debugging a package of mine due to a very strange behaviour that has to do with this, but I finally managed to get this simple example that only involves Rcpp. I don't have much energy left today to investigate further. Any help would be appreciated.

[Enchufa2]

It may also have something to do with #823 and associated issues.

And before I get the usual calling R from C++ is not a good idea, the package in question is this one, and calling R from C++ is the whole point, because it's a simulator and the user needs to be able to insert custom functions into the model.

[Enchufa2]

Equivalent reprex, but arguably more informative:

Rcpp::cppFunction("bool comp(Function fun) { return runif(1)[0] == as<double>(fun()); }")
set.seed(1)
sum(replicate(1e3, comp(function() runif(1, 0, 1))))
# 1000
sum(replicate(1e3, comp(function() runif(1, 0.0001, 0.9999))))
# 0

[eddelbuettel]

I am not sure what else to add besides all that was said in #823:

• for normal use cases, Rcpp does the right thing of setting/resetting RNGs so for standard calls from R to C++ to R you get the same streams in either R or C++; this is important
• when you do something else, for whichever reason, you have to be careful -- that was the Shiny use case in RNG behaves differently in nested calls to Rcpp functions #823, it may be yours here.
• the PR unconditionally get/put RNG state #825 made changes
• we also always had the ability to say Rcpp::export(rng=false) if you want things to be left alone
• similarly, you could use Rcpp and not use our interface generator if they don't suit you

Hope this helps. I understand your simulator framework may have very special needs. But there are also 1830 other packages using Rcpp and they do just fine.

[Enchufa2]

I understand, and I understand that this is not a normal case, statistically speaking; but it is definitely a valid case.

The thing is,

• I'm not sure whether I fully understand why this happens. Rcpp continues the RNG stream and then restores it, and that's why I get the same number?
• Is there any simple workaround using the functions to get/set the RNG state?
• I suppose, then, that if I call the R function stats::runif from Rcpp, instead of Rcpp::runif, I will get two different numbers, because Rcpp will restore/set, restore/set the RNG state for both R calls.
• I'm not sure whether I fully understand the implications of setting Rcpp::export(rng=false). Does this mean that the random stuff called from C++ will continue the RNG stream when coming from R and the other way around?

[Enchufa2]

BTW, the same example with Rcpp::export(rng=false) returns exactly the same, so I'm even more confused than before. :(

[Enchufa2]

And:

Rcpp::cppFunction(
  "bool comp(Function fun) {
    double x = as<double>(fun());
    return runif(1)[0] == x;
  }")
set.seed(1)
sum(replicate(1e3, comp(function() runif(1, 0, 1))))
# 0
sum(replicate(1e3, comp(function() runif(1, 0.0001, 0.9999))))
# 0

seems to work just ok (well, just as I expected at the very least).

[eddelbuettel]

Maybe your use case requires not using our as<>() and wrap(). They are offered as a courtesy, if they don't work for you (in the very special cases of repeatedly crossing the R/C++ border) then maybe you either need to help us debug why, or not use them.

[eddelbuettel]

It would also help if you could restrict yourself to maybe one post at time. These can be edited. I am not sure how I am going this thread up with the parallel questions from you.

The base use of RNGs is described in our intro paper shipped as a vignette. We ensure we have the same streams in normal use cases. You construct a non-normal use case. Things break. Sorry. But that did not disprove the normal use case.

[Enchufa2]

Ok, sorry, one question at a time. Here's an important one: how calling as<>() inline makes any difference? E.g.,

Rcpp::cppFunction(
  "bool comp(Function fun) {
    double x = as<double>(fun());
    return runif(1)[0] == x;
  }")
set.seed(1)
sum(replicate(1e3, comp(function() runif(1, 0, 1))))
# 0
sum(replicate(1e3, comp(function() runif(1, 0.0001, 0.9999))))
# 0

but

Rcpp::cppFunction(
  "bool comp(Function fun) {
    return runif(1)[0] == as<double>(fun());
  }")
set.seed(1)
sum(replicate(1e3, comp(function() runif(1, 0, 1))))
# 1000
sum(replicate(1e3, comp(function() runif(1, 0.0001, 0.9999))))
# 0

[eddelbuettel]

Well it is TMP so other code is involved. My first instinct would be to try a counter example. Do all you are doing with the C API of R. Then add Rcpp one feature at a time and see where it breaks. Or do the inverse and remove Rcpp use until it no longer breaks.

Also, I am not sure how often I have to repeat this: Your function comp(), being passed back and forth 1e3 times via replicate(), is not really what our setup was built for. I understand you are very frustrated, but if you step back a just a bit maybe you can see that you are constructing a special case. Which as it turns out is not covered. Which we knew already.

[Enchufa2]

Sorry, TMP?

replicate has nothing to do with this. In a simple call, without replicate the issue is the same.

[eddelbuettel]

Sorry, TMP is template meta-programming.

[romainfrancois]

@Enchufa2 the difference is the order in which the various runif() are called, i.e. if you save the result of fun() in an object before calling as<> you'll see that as<> does not affect the result:

library(Rcpp)

Rcpp::cppFunction(
  "bool comp1(Function fun) {
    Rcpp::Shield<SEXP> res(fun());
    double x = as<double>(res);
    return runif(1)[0] == x;
  }")
set.seed(1)
sum(replicate(1e3, comp1(function() runif(1, 0, 1))))
#> [1] 0
sum(replicate(1e3, comp1(function() runif(1, 0.0001, 0.9999))))
#> [1] 0

Rcpp::cppFunction(
  "bool comp2(Function fun) {
    Rcpp::Shield<SEXP> res(fun());
    return runif(1)[0] == as<double>(res);
  }")
set.seed(1)
sum(replicate(1e3, comp2(function() runif(1, 0, 1))))
#> [1] 0
sum(replicate(1e3, comp2(function() runif(1, 0.0001, 0.9999))))
#> [1] 0

^{Created on 2019-11-08 by the reprex package (v0.3.0.9000)}

[eddelbuettel]

Nice work. I wanted to convert his example to just doing (scalar) SEXP which is close in spirit, but I didn't think about the various "concurrent" transformations happening when the return object is computed in (what was in his case a compound expression). Also interesting that comp2 works that way via the protection.