<div dir="ltr"><div>Hi Mark,</div><div><br></div><div>many thanks for all the info, I will certainly have a detailed look at what you are doing.</div><div><br></div><div>I think it would be nice to have a package that uses C++ level parallelism (OpenMP or not)</div>
<div>to speed up random number generation at R level. For example:</div><div><br></div><div><div>> library("microbenchmark")</div><div>> library("mvtnorm")</div><div>> library("mvnfast")</div>
<div>> library("MASS")</div></div><div><br></div><div>> # Generating 1e4 50-dimensional multivariate normal r.v.</div><div>> N <- 10000</div><div>> d <- 50</div><div>> mu <- 1:d</div><div>
> tmp <- matrix(rnorm(d^2), d, d)</div><div>> mcov <- tcrossprod(tmp, tmp)</div><div>> myChol <- chol(mcov)</div><div>> </div><div>> microbenchmark(rmvn(N, mu, mcov, ncores = 4),</div><div>+ rmvn(N, mu, mcov),</div>
<div>+ rmvnorm(N, mu, mcov, method = "chol"),</div><div>+ mvrnorm(N, mu, mcov))</div><div><br></div><div><span class="" style="border-collapse:separate;border-spacing:0px;background-color:rgb(225,226,229)"><pre tabindex="0" class="" style="color:rgb(0,0,0);font-family:'Ubuntu Mono';line-height:1.2;outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px;font-size:10.4pt!important;white-space:pre-wrap!important">
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</pre><pre tabindex="0" class="" style="color:rgb(0,0,0);font-family:'Ubuntu Mono';line-height:1.2;outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px;font-size:10.4pt!important;white-space:pre-wrap!important">
<br></pre><pre tabindex="0" class="" style="outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px"><span style="color:rgb(34,34,34);font-family:arial;line-height:normal;background-color:rgb(255,255,255);font-size:10.4pt!important;white-space:pre-wrap!important"> Here I followed Matt's suggestion and used C++11 </span><span style="background-color:rgb(255,255,255);white-space:normal"><font face="arial">mt19937_64 rng, with seeds generated by R::runif.</font></span></pre>
<pre tabindex="0" class="" style="outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px"><span style="background-color:rgb(255,255,255);white-space:normal"><font face="arial">Obviously what you have done by using </font></span><span style="background-color:rgb(255,255,255);font-family:arial,sans-serif;font-size:13px">RngStreams with Boost's distributions is much more rigorous, so probably</span></pre>
<pre tabindex="0" class="" style="outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px"><span style="background-color:rgb(255,255,255);font-family:arial,sans-serif;font-size:13px">that the approach rights to provide faster version of rnorm(), rpois() etc.</span></pre>
<pre tabindex="0" class="" style="outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px"><span style="background-color:rgb(255,255,255);font-family:arial,sans-serif;font-size:13px"><br></span></pre>
<pre tabindex="0" class="" style="outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px"><span style="background-color:rgb(255,255,255);font-family:arial,sans-serif;font-size:13px">Thanks,</span></pre>
<pre tabindex="0" class="" style="outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px"><span style="background-color:rgb(255,255,255);font-family:arial,sans-serif;font-size:13px"><br></span></pre>
<pre tabindex="0" class="" style="outline:none;border:none;word-break:break-all;margin-top:0px;margin-bottom:0px"><span style="background-color:rgb(255,255,255);font-family:arial,sans-serif;font-size:13px">Matteo </span></pre>
</span></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Apr 30, 2014 at 9:10 AM, Mark Clements <span dir="ltr"><<a href="mailto:mark.clements@ki.se" target="_blank">mark.clements@ki.se</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">My two cents worth:<br>
<br>
For the microsimulation package, we needed uniform random number streams<br>
and sub-streams at the C++ level, while supporting R's non-uniform<br>
random number distributions[*]. For this, we used the C++ RngStreams<br>
library and provided "double *user_unif_rand ()" for user-defined RNGs.<br>
Essentially, this reproduces the "L'Ecuyer-CMRG" RNG provided by R, with<br>
seed manipulation at the C++ level (cf. using .Random.seed and R's C<br>
machinery). We also allow for the C++ seeds to be set to and from R.<br>
<br>
See:<br>
<a href="https://github.com/mclements/microsimulation/blob/master/src/microsimulation.cc" target="_blank">https://github.com/mclements/microsimulation/blob/master/src/microsimulation.cc</a><br>
<a href="https://github.com/mclements/microsimulation/blob/master/src/microsimulation.h" target="_blank">https://github.com/mclements/microsimulation/blob/master/src/microsimulation.h</a><br>
(circa lines 290-355)<br>
<a href="https://github.com/mclements/microsimulation/blob/master/R/rcpp_hello_world.R" target="_blank">https://github.com/mclements/microsimulation/blob/master/R/rcpp_hello_world.R</a><br>
(circa lines 6-36)<br>
<br>
This is not implemented using OpenMP. For our purpose, where we are<br>
doing 10^7 small simulations, we chunk the simulations at the R level<br>
and use the parallel package to call the C++ code on each chunk. This<br>
approaches scales well with more processors.<br>
<br>
We have also looked at using the C++ RngStream library with<br>
Boost.Random's and C++11's non-uniform random number distributions.<br>
Again, this has not been implemented using OpenMP (todo?). For a simple<br>
wrapper for Boost, see:<br>
<a href="https://github.com/mclements/microsimulation/blob/master/src/rngstream-boost.hpp" target="_blank">https://github.com/mclements/microsimulation/blob/master/src/rngstream-boost.hpp</a><br>
<br>
with an example:<br>
<a href="https://github.com/mclements/microsimulation/blob/master/src/rngstream-example.cpp" target="_blank">https://github.com/mclements/microsimulation/blob/master/src/rngstream-example.cpp</a><br>
<br>
Two brief notes: first, we have put the RngStream library in the ssim<br>
namespace for use with the microsimulation package. Other uses of the<br>
Boost wrapper would probably want to change the namespace.<br>
<br>
Second, the implementation for C++11 random number generators and<br>
distributions required a small change to the RngStream library -<br>
specifically, to output the random uniform as an unsigned long.<br>
Interestingly, the resulting C++11 random numbers drop every second<br>
value compared with R and Boost.Random.<br>
<br>
For the C+++11 wrapper, see the code RngStream* and rngstream* in:<br>
<a href="https://github.com/mclements/microsimulation/tree/master/test" target="_blank">https://github.com/mclements/microsimulation/tree/master/test</a><br>
<br>
Sincerely, Mark.<br>
<br>
[*] The BH package was not available, C++11 compiler requirements were<br>
not accepted on CRAN, reimplementing non-uniform random number<br>
distributions using UNURAN would have taken too long, and the R<br>
non-uniform random number distributions are extensive and well tested.<br>
</blockquote></div><br></div>