Reputation: 11
I have an RcppArmadillo code that runs sequentially and I am trying to adapt it for parallel compilation with RcppParellel. The code is a copy paste of sequential code and actually runs faster. Yet the results are not correct (probably due to the misuse of the random seeds in parallel cores computing). Can anyone help me?
// Worker
struct UpdateOmegaWorker : public Worker {
const arma::mat X;
const arma::mat Omega_old;
const arma::vec delta;
const double W_slab;
const double W_spike;
const double seed_for_start;
arma::mat& Omega_ret;
// Costruttore
UpdateOmegaWorker(const arma::mat X, const arma::mat Omega_old, const arma::vec delta,
const double W_slab, const double W_spike, arma::mat& Omega_ret, const double seed_for_start): X(X), Omega_old(Omega_old), delta(delta), W_slab(W_slab), W_spike(W_spike), Omega_ret(Omega_ret), seed_for_start(seed_for_start) {}
// Operator
void operator()(std::size_t begin, std::size_t end) {
int n = X.n_rows;
int p = X.n_cols;
arma::mat Xs_it = X;
arma::colvec x_it(n);
arma::colvec phi(n);
arma::colvec pg_aux(n);
arma::colvec kappa_it(n);
arma::colvec zeta_it(n);
arma::colvec xs_tilde(n);
arma::mat node_V(p, p);
arma::colvec node_m(p);
arma::mat Winv = arma::diagmat(1 / (W_slab * delta + W_spike * (1 - delta)));
for (std::size_t s = begin; s < end; s++) {
// HERE IS PROBABILY THE ERROR ?!
arma::arma_rng::set_seed( seed_for_start ) ;
//cores Code (copy Paste of a working code)
// Parametrizzazione
Xs_it.col(s).fill(1.0);
xs_tilde = X.col(s);
// Pòlya Gamma Augmentation
zeta_it = Xs_it * Omega_old.col(s);
pg_aux = cpp_rpg_z_vec(n, 1, zeta_it, 200);
// Ising-Regression
Xs_it.col(s) = xs_tilde;
kappa_it = xs_tilde - 0.5;
node_V = arma::inv(Xs_it.t() * (X.each_col() % pg_aux) + Winv);
node_m = node_V * Xs_it.t() * kappa_it;
Omega_ret.col(s) = cpp_mvrnormArma1(node_m, node_V);
}
}
};
// \[\[Rcpp::export\]\]
arma::mat update_Omega_parallel(const arma::mat X, const arma::mat Omega_old,
const arma::vec delta, const double W_slab, const double W_spike, const double seed_for_start) {
int p = X.n_cols;
arma::mat Omega_ret(p, p, arma::fill::zeros);
UpdateOmegaWorker worker(X, Omega_old, delta, W_slab, W_spike, Omega_ret, seed_for_start);
parallelFor(0, p, worker);
return Omega_ret;
}
I want to fix the seed within cores, using arma functions. E.g., before the operator
std::atomic\<std::size_t\> counter(seed_for_start);
and inside the operator:
arma_rng::set_seed(seed_for_start + counter++);
Upvotes: 1
Views: 74