Skip to content

[Rcpp-devel] Mersenne Twister in RcppArmadillo?

7 messages · Yan Zhou, Dirk Eddelbuettel, Darren Cook +1 more

Original
Yan Zhou
# 
It runs in parallel does not mean it works correctly.

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.

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.

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:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        
_______________________________________________
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130214/b8a53809/attachment.html>

  
  


  


    

    

      
      

        


  

        

      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

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Darren Cook
    

    

      
      #
    

  


  

          

      
      
      
    

        
BTW, the O'Reilly Parallel R book mentions random number generation in
quite a few places. Reading it, you really feel nostalgic for the simple
life of good 'ole single-threaded random number generators :-)

Darren

  
  


  


      

    

    

      
      

        


  

        

      Simon Zehnder
    

    

      
      #
    

  


  

      
Hi Yan,

thank you very much for your comment. I read now the section 6 of the parallel package documentation. 

It's better to give programming a little more time than to get at the end undefined behavior. At this point only a short question: I use the nlopt library, constructing in every thread an nlopt::opt object. That should be thread-safe right as I have an object for every thread (the void data pointer is constructed in every thread and it is not shared)?

Best

Simon

        
On Feb 14, 2013, at 11:24 PM, Yan Zhou <zhouyan at me.com> wrote:

        

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20130215/41250d42/attachment.html>

  
  


  


      

    

    

      
      

        


  

        

      Simon Zehnder
    

    

      
      #
    

  


  

      
Hi Dirk,

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?

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)?

        
On Feb 15, 2013, at 12:08 AM, Dirk Eddelbuettel <edd at debian.org> wrote:

        

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Simon Zehnder
    

    

      
      #
    

  


  

      
Hi Darren,

I haven't known about this book. I own the R cookbook from O'Reilly and this is quite well written. I will give it a closer view soon! Thanks for this info!

Best

Simon

        
On Feb 15, 2013, at 1:58 AM, Darren Cook <darren at dcook.org> wrote:

        

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Dirk Eddelbuettel
    

    

      
      #
    

  


  

      
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  
|