Deployment of sitmo within C++ Code

2018-06-28

Within this vignette, details on how to use sitmo’s header will be detailed. First, the background on sitmo will be provided. Secondly, function calls will be shown alongside of a description. Thirdly, examples will be provided of how one can use the sitmo header.

What is sitmo and can I eat it?

sitmo is the consultancy agency founded by Thijs van den Berg. They first released a Parallel Psuedo Random Number Generator (PPRNG) under the same name using work in Salmon, K., et al.’s “Parallel Random Numbers: As Easy as 1, 2, 3” in the conference proceedings of the 2011 International Conference for High Performance Computing, Networking, Storage and Analysis. Support for sitmo exists for both C++ standards: C++98 and C++11. Furthermore, there are many different PPRNGs that are available: trng, SPRNG, RngStreams, OMPRNG. However, none are as appealing in my eyes than sitmo, which provides a straight forward interface to generating psuedo-random numbers (RNG), the least restrictive license (MIT), and speed.

Over the span of the last few years, the sitmo agency has released two other engines of interest for C++11: threefry and vandercorput. The threefry engine is a rewritten PPRNG version of sitmo for C++11 whereas the vandercorput engine provides one dimensional low-discrepancy sequencing. The latter engines are also available under the MIT license.

Accessing and using engines in sitmo

The header files for sitmo, threefry, and vandercorput engines are contained within this package. To use one of these engine header files within your own package, you can link to the sitmo package within your description file. e.g.

LinkingTo: Rcpp, sitmo
Imports:
    Rcpp (>= 0.12.11)

To use C++11’s statistical distributions, you may want to add the following to your src/Makevars and src/Makevars.win file:

CXX_STD = CXX11

Within a C++ file in src/, then add:

#include <Rcpp.h>
#include <sitmo.h>      // SITMO for C++98 & C++11 PPRNG
#include <threefry.h>   // THREEFRY C++11-only PPRNG
#include <vandercorput> // VANDERCORPUT C++11-only Low-discrepancy sequence

You do not need to add each header file. Pick and choose the appropriate engine for your needs.

Or you can do a direct embed in your application. I would advise for the prior though and, hence, the reason for this package.

Below is a breakdown of functions that are available for the engines. Please note, that the engine predominantly highlight is the original: sitmo::prng.

Construct an engine

Expression Description
prng() Creates an engine with a default initial state.
prng(prng& x) Creates an engine with the same initial state as the engine x.
prng(uint32_t s) Creates an engine with initial state determined by s. Engines created with different initial states have the guarantee to generate independent non-overlapping random sequences of length \(2^128\).
prng(SeedSeq q) Creates an engine with an initial state that depends on a sequence produced by one call to q.generate.

Seed modifiers

To use the seed modifiers, one must first construct an engine using a method detailed in the previous table.

From there, the engine state can be modified using:

Expression Description
e.seed() Returns the random engine to the default state. The same prng().
e.seed(uint32_t s) Set the engine to a state determined by s. Same as prng(uint32_t s)
e.seed(SeedSeq q) Set the engine to a state that depends on a sequence produced by one call to q.generate. Same as prng(SeedSeq q)
e() Advances the internal state and returns a 32 bit random number.
e.discard(uint64_t n) Advances the internal state with n steps in constant time.

Misc Seed

Using the same engine created above, one can access additional state information using the following:

Expression Description
e1 == e2 Test for equivalence of two prng’s. Two engines are the same if they generate the exact same random sequence.
e1 != e2 Test for non-equivalence of two prng’s. Two engines are different if they generate different random sequences.
e.version() The current version of the engine, returns the value 2

Examples

The examples displayed in the vignette are taken directly from the project’s src directory that is found here https://github.com/coatless/sitmo. Additional commentary is added.

Uniform Random Number Generator

Under C++98, one does not have access to the C++11 implementation of the Uniform distribution. This is particularly problematic as a lot of the distribution RNG rely upon being able to sample from \(\left[0,1\right]\) ala the Probability Integral Transformation Theorem. Additional details are discussed in a separate vignette (“Making a Uniform PRNG with sitmo”).

OpenMP Example

One of the primary reasons why sitmo is desirable is because it can be used under parallelization via OpenMP and MPI. Below is an example where it is used in a parallel setting to generate numbers. Note, to ensure that code works cross-platform, please protect against OpenMP includes as the package will otherwise fail on OS X.

To protect against a lack of OpenMP headers use:

When writing sections of parallelized code, also protect that code using:

Furthermore, add the following to your Makevars and Makevars.win:

With this being said, let’s take a look at an example parallelization using sitmo:

#include <Rcpp.h>
#include <sitmo.h> // SITMO PPRNG

// [[Rcpp::depends(sitmo)]]

#ifdef _OPENMP
#include <omp.h>
#endif

// [[Rcpp::plugins(openmp)]]

//' Test Generation using sitmo and C++11
//' 
//' The function provides an implementation of creating realizations from the default engine.
//' 
//' @param n An \code{unsigned integer} denoting the number of realizations to generate.
//' @param seeds A \code{vec} containing a list of seeds. Each seed is run on its own core.
//' @return A \code{vec} containing the realizations.
//' @details
//' The following function's true power is only accessible on platforms that support OpenMP (e.g. Windows and Linux).
//' However, it does provide a very good example as to how to make ones code applicable across multiple platforms.
//' 
//' With this being said, how we determine how many cores to split the generation to is governed by the number of seeds supplied.
//' In the event that one is using OS X, only the first seed supplied is used. 
//' 
//' @export
//' @examples
//' a = sitmo_parallel(10, 5.0, c(1))
//' 
//' b = sitmo_parallel(10, 5.0, c(1))
//' 
//' c = sitmo_parallel(10, 5.0, c(2))
//' 
//' isTRUE(all.equal(a,b))
//' 
//' isTRUE(all.equal(a,c))
// [[Rcpp::export]]
Rcpp::NumericVector sitmo_parallel(unsigned int n, Rcpp::NumericVector& seeds){ 
  
  unsigned int ncores = seeds.size();

  Rcpp::NumericVector q(n);
  
  #ifdef _OPENMP
  #pragma omp parallel num_threads(ncores) if(ncores > 1)
  {
  #endif
  
    // Engine requires uint32_t inplace of unsigned int
    uint32_t active_seed;
      
    // Write the active seed per core or just write one of the seeds.
    #ifdef _OPENMP
      active_seed = static_cast<uint32_t>(seeds[omp_get_thread_num()]);
    #else
      active_seed = static_cast<uint32_t>(seeds[0]);
    #endif
    
    sitmo::prng eng( active_seed );
  
    // Parallelize the Loop
    #ifdef _OPENMP
    #pragma omp for schedule(static)
    #endif
    for (unsigned int i = 0; i < n; i++){
      q[i] = eng(). 
    }

  #ifdef _OPENMP
  }
  #endif
  
  return q;
}