Hi Mark, many thanks for all the info, I will certainly have a detailed look at what you are doing.
I think it would be nice to have a package that uses C++ level parallelism (OpenMP or not) to speed up random number generation at R level. For example: > library("microbenchmark") > library("mvtnorm") > library("mvnfast") > library("MASS") > # Generating 1e4 50-dimensional multivariate normal r.v. > N <- 10000 > d <- 50 > mu <- 1:d > tmp <- matrix(rnorm(d^2), d, d) > mcov <- tcrossprod(tmp, tmp) > myChol <- chol(mcov) > > microbenchmark(rmvn(N, mu, mcov, ncores = 4), + rmvn(N, mu, mcov), + rmvnorm(N, mu, mcov, method = "chol"), + mvrnorm(N, mu, mcov)) Unit: milliseconds expr min lq median uq max neval rmvn(N, mu, mcov, ncores = 4) 8.339534 12.85869 13.23915 13.82640 38.04688 100 rmvn(N, mu, mcov) 31.645092 32.32220 33.70345 34.03561 57.78102 100 rmvnorm(N, mu, mcov, method = "chol") 57.956183 58.99730 59.84668 82.11610 90.71133 100 mvrnorm(N, mu, mcov) 60.621391 62.06400 82.86829 85.42350 90.50662 100 Here I followed Matt's suggestion and used C++11 mt19937_64 rng, with seeds generated by R::runif. Obviously what you have done by using RngStreams with Boost's distributions is much more rigorous, so probably that the approach rights to provide faster version of rnorm(), rpois() etc. Thanks, Matteo On Wed, Apr 30, 2014 at 9:10 AM, Mark Clements <mark.cleme...@ki.se> wrote: > My two cents worth: > > For the microsimulation package, we needed uniform random number streams > and sub-streams at the C++ level, while supporting R's non-uniform > random number distributions[*]. For this, we used the C++ RngStreams > library and provided "double *user_unif_rand ()" for user-defined RNGs. > Essentially, this reproduces the "L'Ecuyer-CMRG" RNG provided by R, with > seed manipulation at the C++ level (cf. using .Random.seed and R's C > machinery). We also allow for the C++ seeds to be set to and from R. > > See: > > https://github.com/mclements/microsimulation/blob/master/src/microsimulation.cc > > https://github.com/mclements/microsimulation/blob/master/src/microsimulation.h > (circa lines 290-355) > > https://github.com/mclements/microsimulation/blob/master/R/rcpp_hello_world.R > (circa lines 6-36) > > This is not implemented using OpenMP. For our purpose, where we are > doing 10^7 small simulations, we chunk the simulations at the R level > and use the parallel package to call the C++ code on each chunk. This > approaches scales well with more processors. > > We have also looked at using the C++ RngStream library with > Boost.Random's and C++11's non-uniform random number distributions. > Again, this has not been implemented using OpenMP (todo?). For a simple > wrapper for Boost, see: > > https://github.com/mclements/microsimulation/blob/master/src/rngstream-boost.hpp > > with an example: > > https://github.com/mclements/microsimulation/blob/master/src/rngstream-example.cpp > > Two brief notes: first, we have put the RngStream library in the ssim > namespace for use with the microsimulation package. Other uses of the > Boost wrapper would probably want to change the namespace. > > Second, the implementation for C++11 random number generators and > distributions required a small change to the RngStream library - > specifically, to output the random uniform as an unsigned long. > Interestingly, the resulting C++11 random numbers drop every second > value compared with R and Boost.Random. > > For the C+++11 wrapper, see the code RngStream* and rngstream* in: > https://github.com/mclements/microsimulation/tree/master/test > > Sincerely, Mark. > > [*] The BH package was not available, C++11 compiler requirements were > not accepted on CRAN, reimplementing non-uniform random number > distributions using UNURAN would have taken too long, and the R > non-uniform random number distributions are extensive and well tested. >
_______________________________________________ Rcpp-devel mailing list Rcpp-devel@lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel