Random number generation with Rcpp and OpenMP

Random number generation with Rcpp and OpenMP

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.

Ignacio Martinez avatar
About Ignacio Martinez
research economist, tech enthusiast
comments powered by Disqus