I'm trying to learn how to draw random numbers from a normal and a binomial distribution from within a Rcpp OpenMP loop.
I wrote the following code using R::rnorm and R::rbinom which I understand to be a be a don't do that.
#include <RcppArmadillo.h>
#include <omp.h>
#include <dqrng.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(dqrng)]]
// [[Rcpp::plugins(openmp)]]
using namespace Rcpp;
// [[Rcpp::export]]
arma::mat my_matrix3(const arma::mat& A,
const arma::mat& B) {
dqrng::dqRNGkind("Xoroshiro128+");
dqrng::dqset_seed(42);
const int nObservations = A.n_cols;
const int nDraws = B.n_rows;
const int nRows = nObservations * nDraws;
// Show initialization information
Rcpp::Rcout << "nObservations: " << nObservations << std::endl
<< "nDraws: " << nDraws << std::endl
<< "nRows: " << nRows << std::endl;
arma::mat out(nRows, 5);
// Show trace of matrix construction
Rcpp::Rcout << "out - rows: " << out.n_rows << std::endl
<< "out - columns: " << out.n_cols << std::endl;
int i,n,iter,row;
double dot_product, random_number, p;
omp_set_num_threads(2);
#pragma omp parallel for private(i, n, iter, row)
for(i = 0; i < nDraws; ++i){
for(n = 0; n < nObservations; ++n) {
row = i * nObservations + n;
dot_product = arma::as_scalar(A.col(n).t() * B.row(i).t());
// Show trace statement of index being accessed
out(row, 0) = i + 1;
out(row, 1) = n + 1;
out(row, 2) = R::rnorm(2.0, 1.0) ;//dqrng::dqrnorm(1, 2.0, 1.0)[1];
out(row, 3) = 1 / (1 + std::exp(-out(row, 3) - std::exp(dot_product)));
out(row, 4) = R::rbinom(1,out(row, 3));
}
}
return out;
}
/*** R
set.seed(9782)
A <- matrix(rnorm(10), nrow = 5)
B <- matrix(rnorm(10), ncol = 5)
test <- my_matrix3(A = A, B = B)
test
mean(test[,5])
*/
As suggested by @ralf-stubner I'm trying to replace that code with dqrng. If I understood the documentation correctly, I can replace R::rnorm(2.0, 1.0) with dqrng::dqrnorm(1, 2.0, 1.0)[1]. Is that right? What about replacing R::rbinom(1,out(row, 3))? I was not able to find any reference to drawing from a binomial in the documentation
I was able to write the following code that draws from a binomial distribution in parallel:
#include <RcppArmadillo.h>
// [[Rcpp::depends(dqrng, BH, RcppArmadillo)]]
#include <pcg_random.hpp>
#include <boost/random/binomial_distribution.hpp>
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
arma::mat parallel_random_matrix(int n, int m, int ncores, double p=0.5) {
//dqrng::uniform_distribution dist(0.0, 1.0);
boost::random::binomial_distribution<int> dist(1,p);
dqrng::xoshiro256plus rng(42);
arma::mat out(n,m);
// ok to use rng here
#pragma omp parallel num_threads(ncores)
{
dqrng::xoshiro256plus lrng(rng); // make thread local copy of rng
lrng.jump(omp_get_thread_num() + 1); // advance rng by 1 ... ncores jumps
auto gen = std::bind(dist, std::ref(lrng));
#pragma omp for
for (int i = 0; i < m; ++i) {
double lres(0);
for (int j = 0; j < n; ++j) {
out(j,i) = gen(); /// CAN I MAKE P BE A FUNCTION OF i and j? how???
}
}
}
// ok to use rng here
return out;
}
/*** R
parallel_random_matrix(5, 5, 4, 0.75)
*/
> parallel_random_matrix(5, 5, 4, 0.75)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 1 1 1 1
[2,] 0 1 1 1 1
[3,] 0 1 0 1 0
[4,] 1 1 1 1 1
[5,] 1 1 1 1 1
Is there a way to call out(j,i) = gen(); and make the probability be a function of i and j??
Is there anything wrong with the code I wrote?
A simple solution would be to create a new distribution object within the loop:
// [[Rcpp::depends(dqrng, BH, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <boost/random/binomial_distribution.hpp>
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::export]]
arma::mat parallel_random_matrix(int n, int m, int ncores, double p=0.5) {
dqrng::xoshiro256plus rng(42);
arma::mat out(n,m);
// ok to use rng here
#pragma omp parallel num_threads(ncores)
{
dqrng::xoshiro256plus lrng(rng); // make thread local copy of rng
lrng.jump(omp_get_thread_num() + 1); // advance rng by 1 ... ncores jumps
#pragma omp for
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
// p can be a function of i and j
boost::random::binomial_distribution<int> dist(1,p);
auto gen = std::bind(dist, std::ref(lrng));
out(j,i) = gen();
}
}
}
// ok to use rng here
return out;
}
/*** R
parallel_random_matrix(5, 5, 4, 0.75)
*/
This way you can make p depend on i and j. It might be more efficient to keep one global dist object and reconfigure it within the loop, see here for a similar question.