使用openMP和Rcpp设置RNG状态 [英] set RNG state with openMP and Rcpp

问题描述

我有一个澄清的问题. 据我了解，sourceCpp自动传递RNG状态，因此set.seed(123)在调用Rcpp代码时为我提供了可重现的随机数.编译程序包时，我必须添加一个设置的RNG语句. 现在，这一切如何在sourceCpp或程序包中与openMP一起使用?

I have a clarification question. It is my understanding, that sourceCpp automatically passes on the RNG state, so that set.seed(123) gives me reproducible random numbers when calling Rcpp code. When compiling a package, I have to add a set RNG statement. Now how does this all work with openMP either in sourceCpp or within a package?

考虑以下Rcpp代码

#include <Rcpp.h>
#include <omp.h>

// [[Rcpp::depends("RcppArmadillo")]]

// [[Rcpp::export]]
Rcpp::NumericVector rnormrcpp1(int n, double mu, double sigma  ){
  Rcpp::NumericVector out(n);
        for (int i=0; i < n; i++) {
          out(i) =R::rnorm(mu,sigma);
        }

  return(out);
}


// [[Rcpp::export]]
Rcpp::NumericVector rnormrcpp2(int n, double mu, double sigma, int cores=1  ){
  omp_set_num_threads(cores);
  Rcpp::NumericVector out(n);
    #pragma omp parallel for schedule(dynamic)   
        for (int i=0; i < n; i++) {
          out(i) =R::rnorm(mu,sigma);
        }

  return(out);
}

然后运行

   set.seed(123)
    a1=rnormrcpp1(100,2,3,2)
    set.seed(123)
    a2=rnormrcpp1(100,2,3,2)
    set.seed(123)
    a3=rnormrcpp2(100,2,3,2)
    set.seed(123)
    a4=rnormrcpp2(100,2,3,2)
    all.equal(a1,a2)
    all.equal(a3,a4)

虽然a1和a2相同，但a3和a4不同.如何使用openMP循环调整RNG状态?我可以吗?

While a1 and a2 are identical, a3 and a4 are not. How can I adjust the RNG state with the openMP loop? Can I?

推荐答案

是的，sourceCpp()等以及RNGScope的实例化，因此RNG处于适当的状态.

Yes, sourceCpp() etc and an instantiation of RNGScope so the RNGs are left in a proper state.

是的，可以执行OpenMP.但是在OpenMP段中，您无法控制线程执行的顺序-因此您需要延长相同的序列.我在开发一个软件包时遇到了同样的问题，我想在其中使用可重现的抽签，但仍使用OpenMP.但是看来你做不到.

And yes one can do OpenMP. But inside of OpenMP segment you cannot control in which order the threads are executed -- so you longer the same sequence. I have the same problem with a package under development where I would like to have reproducible draws yet use OpenMP. But it seems you can't.

这篇关于使用openMP和Rcpp设置RNG状态的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！