coatless-rpkg/sitmo

How to set random seed using RcppParallel or OpenMP to repeatable

xiangpin opened this issue · 1 comments

Good Day

Thanks for developing the package, I using shuffle of Armadillo with RcppParallel or OpenMP , but the result is not repeatable. I found this package might solve the problem. To be reproducible, is there a function like set.seed of R base to do this? Thank you.

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <omp.h>

using namespace arma;


arma::rowvec calkld(arma::rowvec z, arma::rowvec bg){
    arma::rowvec x = z % log(z / bg);
    return x;

}


arma::rowvec calmeankld(arma::rowvec x, arma::rowvec bg, int permutation = 100){

    arma::vec bootkld(permutation);

    for (int j = 0; j < permutation; j++){
        arma::rowvec z = arma::shuffle(x);
        arma::rowvec k = calkld(z, bg);
        bootkld[j] = log(sum(k));
    }

    double bmean = arma::mean(bootkld);
    double bsd = arma::stddev(bootkld);
    arma::rowvec res = {bmean, bsd};

    return res;
}

//[[Rcpp::export]]
arma::mat testparallel(arma::mat m, int permutation = 100){
  m = m + 1e-300;
  int n = m.n_rows;
  arma::mat out(m.n_rows, 2);

  arma::rowvec bg = mean(m, 0);
  bg = bg + 1e-300;

  // Parallelize the Loop
  #pragma omp parallel for schedule(static)
  for (unsigned int i = 0; i < n; i++){
    out.row(i) = calmeankld(m.row(i), bg, permutation);
  }

  return out;
}

> da <- matrix(abs(rnorm(200)), ncol=10)
> withr::with_seed(123, testparallel(da, permutation=100)) -> t1
> withr::with_seed(123, testparallel(da, permutation=100)) -> t2
> identical(t1, t2)

[1] FALSE

> RhpcBLASctl::omp_set_num_threads(1)
> withr::with_seed(123, testparallel(da, permutation=100)) -> t3
> withr::with_seed(123, testparallel(da, permutation=100)) -> t4
> identical(t3, t4)

[1] TRUE

Closing as this is a cross-post and Ralf was kind to provide a worked example.

daqana/dqrng#83