How to set random seed using RcppParallel or OpenMP to repeatable
xiangpin opened this issue · 1 comments
xiangpin commented
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
coatless commented
Closing as this is a cross-post and Ralf was kind to provide a worked example.