eddelbuettel / rcppziggurat

Rcpp bindings for different Ziggurat RNG implementations
Other
12 stars 4 forks source link

sourceCpp seems to reset random state #13

Closed luc-vermeylen closed 4 years ago

luc-vermeylen commented 4 years ago

Hi,

I found that when using sourceCpp, the random state is reset. For example, let's take this function:

#include <Rcpp.h>
#include <Ziggurat.h>

using namespace Rcpp;
static Ziggurat::Ziggurat::Ziggurat zigg;

// [[Rcpp::export]]
double rand_number() {
  zigg.norm();
  return zigg.norm();
}

If i then source via sourceCpp and run the function i get the same random number each time. Only when i re-run the function without re-sourcing, i get different numbers.

How can i change this behavior so that also upon resourcing i get a new random number rather than the same? Is this intentional? If so, whats the reasoning and how can i change it?

Kind regards, Luc

eddelbuettel commented 4 years ago

Good catch. That is normal (desired) behaviour for the "usual" R RNGs to ensure proper handshake. But because there can be corner cases where one does not want it, or if one wants to shave a microsecond, or ... we accomodated this added an option.

But I am once again having a hard time finding where exactly we documented that :-/ Hold one for a sec, and lemme poke @kevinushey ...

eddelbuettel commented 4 years ago

No, wait, that is something else. I misunderstood the issue.

Let me think about this for a second, but I have the suspicion that your use of the static instance may enforce just that.

eddelbuettel commented 4 years ago

Well the vignette show the usage as you have it plus the added few lines of

// [[Rcpp::export]]
void zsetseed(unsigned long int s) {
    zigg.setSeed(s);
    return;
}

You could simply accomodate your non-package use via sourceCpp() by adding a block such as this

/*** R
zsetseed(as.numeric(Sys.time())
*/

Demo

R> Rcpp::sourceCpp("/tmp/zigg.cpp")

R> zsetseed(as.numeric(Sys.time()))

R> rand_number()
[1] 0.145456
R> Rcpp::sourceCpp("/tmp/zigg.cpp")

R> zsetseed(as.numeric(Sys.time()))

R> rand_number()
[1] 0.357142
R> Rcpp::sourceCpp("/tmp/zigg.cpp")

R> zsetseed(as.numeric(Sys.time()))

R> rand_number()
[1] 0.115685
R> 

Good enough?

luc-vermeylen commented 4 years ago

Works for me!

Thank you very much for the quick response!

Kind regards, Luc

eddelbuettel commented 4 years ago

And the "real" real reason is, of course, that Ziggurat has a fixed seed in its code.

Hence the need for a setter. Glad you're set!