[Rcpp-devel] Mersenne Twister in RcppArmadillo?

[Dirk Eddelbuettel]

On 14 February 2013 at 22:24, Yan Zhou wrote:
| It runs in parallel does not mean it works correctly.

Well put :)

| Simply increase seed certainly not avoid overlapping. The best algorithm that
| use MT19937 in parallel environment without overlapping has a complexity O(log
| n) AFAIK.
| 
| Second, R's RNG has a single instance within the whole program, which means
| your program does not work as you thought. There are obvious race conditions in
| your approach.

Precisely -- First rule of OpenMP with R: never ever call back into R from a
thread. Which implies to not access Rcpp::Function etc.
 
| I suggest you google a little bit of how RNG really works before purse how to
| generate them in parallel. For example section 6 of the R's parallel package's
| document.

I meant to make the same suggestion about the same carefully written
document. Prof Ripley also has an (older) bug on this, he _really_ knows this
stuff. 

Regards, Dirk
 
| Best,
| 
| Yan Zhou
| 
| On Feb 14, 2013, at 09:40 PM, Simon Zehnder <szehnder at uni-bonn.de> wrote:
| 
| 
|     Well, in OpenMP it seems to work, if I use inside the #pragma parallel for:
| 
|     Rcpp::Environment base("package:base");
|     Rcpp::Function SetSeed = base["set.seed"];
|     SetSeed(1 + i * nobs_intern * 3);
| 
|     I used OMP_NUMTHREADS = 4 and 8 iterations.
| 
|     Best Simon
| 
|     On Feb 14, 2013, at 10:10 PM, Simon Zehnder <szehnder at uni-bonn.de> wrote:
| 
|     > I made all my simple tests now by using Rcpp::Environment and then
|     Rcpp::Function for "set.seed" when calling R::rnorm. The next step would be
|     to parallelize the iterations via OpenMP. The suggestion of Yan makes quite
|     sense to produce a RNG via std::mt19937. This is also the way I did it when
|     I used Scythe Statistical Library.
|     > As I use now the R RNG it should be able to use it in parallel, if
|     Rcpp::Function is a kind of wrapper which contains the code of "set.seed".
|     But I would rather guess, that it is a call object, that just calls the R
|     function "set.seed" in R, which then would not be threadsafe (as also the
|     RNG would then just be called when using R::rnorm).....
|     >
|     > So I guess falling back to std::mt19937 is not a bad idea. For no
|     overlapping I just increase the seed in every iteration dependent on the
|     iteration and the number of random numbers to generate.
|     >
|     > Best Simon
|     >
|     >
|     >
|     >
|     > On Feb 11, 2013, at 2:16 PM, Chris Jefferson <chris at bubblescope.net>
|     wrote:
|     >
|     >> On 11/02/13 10:23, c s wrote:
|     >>> On Sun, Feb 10, 2013 at 12:32 AM, Yan Zhou <zhouyan at me.com> wrote:
|     >>>> To have Armadillo randn use MT19937 is not easy.
|     >>>> Since it use srand() for seed, I guess it also use C rand(),
|     >>>> whose quality is quite questionable.
|     >>> The quality of the rand() function from C depends on the
|     >>> implementation in libc, which varies from system to system.
|     >>>
|     >>> While I'd like to incorporate a Mersenne-Twister random number
|     >>> generator directly in Armadillo, it would either add a dependency on
|     >>> Boost, or on C++11. Boost might not be available on a user's system,
|     >>> and the degree of support for C++11 features varies from compiler to
|     >>> compiler. We also have to bear in mind that R folks currently
|     >>> disallow CRAN packages that use C++11.
|     >>
|     >> If you just want a mersenne-twister random number generator, I will
|     extract it from boost (which is under a fairly free licence).
|     >>
|     >> Chris
|     >> _______________________________________________
|     >> Rcpp-devel mailing list
|     >> Rcpp-devel at lists.r-forge.r-project.org
|     >> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
|     >
|     > _______________________________________________
|     > Rcpp-devel mailing list
|     > Rcpp-devel at lists.r-forge.r-project.org
|     > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
| 
|     _______________________________________________
|     Rcpp-devel mailing list
|     Rcpp-devel at lists.r-forge.r-project.org
|     https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
| 
| 
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com