[Rcpp-devel] Mersenne Twister in RcppArmadillo?

Dirk Eddelbuettel edd at debian.org
Fri Feb 15 13:31:03 CET 2013 

Hi Simon,

On 15 February 2013 at 08:56, Simon Zehnder wrote:
| that is good to know. At this point I would like to know: Functions like R::qbinom or R::qnorm from R's C library can still be used right? They are not called, but the code is read and worked through in compilation?

That should work. You just can't get to the single-threaded "R engine" from
multiple-threads. Anything stateful will go funky, or even just blow up --
and the RNG is stateful.

Linking through should be fine. 
 
| Best
| 
| Simon
| 
| P.S.: How is it with code using RcppArmadillo: Is there a general interest to share reusable code (e.g. I programmed for my GMM estimation a Newey West covariance estimation procedure, reproducing results from R's package sandwich with prewhitening, bandwith bwNeweyWest and weights weightsAndrews. Also sandwich estimation is possible. It becomes much faster than the R equivalent when working with datasets bigger than 50.000 observations)?

Sounds good. 

You may want to think about packaging that. Already 30 packages on CRAN use
RcppArmadillo and provide examples of using RcppArmadillo in a package.

Dirk
 
| 
| On Feb 15, 2013, at 12:08 AM, Dirk Eddelbuettel <edd at debian.org> wrote:
| 
| > 
| > 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  
| 

-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com

More information about the Rcpp-devel mailing list