Random number generation with Rcpp and OpenMP
Page content
The following code shows how to write some simple code to draw random numbers from a normal and a binomial distribution. Notice that instead of declaring A as a numeric matri
Serial
Double loop
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix my_matrix(int I) {
NumericMatrix A(I,2);
for(int i = 0; i < I; i++){
A(i,0) = R::rnorm(2,1) ;
A(i,1) = R::rbinom(1,0.5) ;
}
colnames(A) = CharacterVector::create("Normal", "Bernoulli");
return A;
}
set.seed(42)
my_matrix(5)
## Normal Bernoulli
## [1,] 3.370958 0
## [2,] 2.955936 1
## [3,] 2.632863 1
## [4,] 2.539024 1
## [5,] 3.511522 0
Vectorized
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix my_matrix2(int I) {
NumericMatrix A(I,2);
A(_,0) = Rcpp::rnorm(I,2,1) ;
A(_,1) = Rcpp::rbinom(I,1,0.5) ;
colnames(A) = CharacterVector::create("Normal", "Bernoulli");
return A;
}
set.seed(42)
my_matrix2(5)
## Normal Bernoulli
## [1,] 3.370958 0
## [2,] 1.435302 1
## [3,] 2.363128 1
## [4,] 2.632863 0
## [5,] 2.404268 0
Parallel
R’s RNG is another “don’t do that” when doing parallel computing. Instead, I’m using dqrng for RNGs which supports parallel code.
// [[Rcpp::depends(dqrng, BH, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <boost/random/binomial_distribution.hpp>
#include <xoshiro.h>
#include <dqrng_distribution.h>
#include <omp.h>
// [[Rcpp::plugins(openmp)]]
// [[Rcpp::plugins(cpp11)]]
using namespace Rcpp;
// aliases for the used ditributions
using binomial = boost::random::binomial_distribution<int>;
using normal = dqrng::normal_distribution;
// [[Rcpp::export]]
NumericMatrix my_parallel_matrix(int I, int nthreads, double p=0.5) {
dqrng::xoshiro256plus rng(42); // SET SEED
arma::mat A(I,2);
#pragma omp parallel num_threads(nthreads)
{
int i;
dqrng::xoshiro256plus lrng(rng); // make thread local copy of rng
lrng.jump(omp_get_thread_num() + 1); // advance rng by 1 ... nthreads jumps
// distributions with default parameters
binomial bernoulli;
normal normal1;
#pragma omp for
for(i = 0; i < I; i++){
// p could be a function of i. Same thing for mu and sigma.
A(i,0) = normal1(rng, normal::param_type(p, 1.0));
A(i,1) = bernoulli(rng, binomial::param_type(1, p));
}
}
NumericMatrix out(I,2);
out = wrap(A);
colnames(out) = CharacterVector::create("Normal", "Bernoulli");
return out;
}
my_parallel_matrix(5,4)
## Normal Bernoulli
## [1,] 0.9632770 0
## [2,] 0.5000870 0
## [3,] 0.6958893 0
## [4,] 0.9632770 0
## [5,] 0.3481794 0
Benchmark
library(microbenchmark)
library(ggplot2)
mb <- microbenchmark(my_matrix(100000),
my_matrix2(100000),
my_parallel_matrix(100000,5))
autoplot(mb)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
