`sitmo`

within C++ CodeWithin 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.

`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.

`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`

.

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. |

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

```
// Generate engine called eng_org.
sitmo::prng eng_org;
// Generate engine called eng_org.
sitmo::threefry eng_tf;
// Generate engine called eng_vc.
sitmo::vandercorput eng_vc;
```

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. |

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 |

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.

Below is the most rudimentary example using `sitmo`

. It describes the process of creating a `sitmo`

engine and obtaining a draw.

```
#include <Rcpp.h>
#include <sitmo.h> // SITMO PPRNG
// [[Rcpp::depends(sitmo)]]
//' Example RNG Draws with sitmo
//'
//' Shows a basic setup and use case for sitmo.
//'
//' @param n A \code{unsigned int} is a .
//' @return A \code{vec} with random sequences.
//' @examples
//' n = 10
//' sitmo_draws(n)
// [[Rcpp::export]]
Rcpp::NumericVector sitmo_draws(unsigned int n) {
Rcpp::NumericVector o(n);
// Create a prng engine
sitmo::prng eng;
// Draw from base engine
for (unsigned int i=0; i< n ; ++i){
o(i) = eng();
}
return o;
}
```

Here the ability for a seed to be set is used.

```
#include <Rcpp.h>
#include <sitmo.h> // SITMO PPRNG
// [[Rcpp::depends(sitmo)]]
//' Example Seed Set and RNG Draws with sitmo
//'
//' Shows how to set a seed in sitmo.
//'
//' @param n An \code{unsigned int} that dictates how many realizations occur.
//' @param seed An \code{unsigned int} that controls the rng seed.
//' @return A \code{vector} with random sequences.
//' @examples
//' n = 10
//' a = sitmo_engine_seed(n, 1337)
//' b = sitmo_engine_seed(n, 1337)
//' c = sitmo_engine_seed(n, 1338)
//'
//' isTRUE(all.equal(a,b))
//' isTRUE(all.equal(a,c))
// [[Rcpp::export]]
Rcpp::NumericVector sitmo_engine_seed(unsigned int n, unsigned int seed) {
// Create Rcpp Matrix
Rcpp::NumericVector o(n);
// Create a prng engine with a specific seed
sitmo::prng eng(static_cast<uint32_t>(seed));
// Draw from base engine
for (unsigned int i=0; i < n; ++i){
o(i) = eng();
}
return o;
}
```

The code used here can be found to work when the initial state of the engine needs to be reverted.

```
#include <Rcpp.h>
#include <sitmo.h> // SITMO PPRNG
// [[Rcpp::depends(sitmo)]]
//' Example Seed Set and RNG Draws with sitmo
//'
//' Shows how to set a seed in sitmo.
//'
//' @param n An \code{unsigned int} that dictates how many realizations occur.
//' @param seed An \code{unsigned int} that controls the rng seed.
//' @return A \code{matrix} with random sequences.
//' @examples
//' n = 10
//' a = sitmo_engine_seed(n, 1337)
//'
//' isTRUE(all.equal(a[,1],a[,2]))
// [[Rcpp::export]]
Rcpp::NumericMatrix sitmo_engine_reset(unsigned int n, unsigned int seed) {
// Create Rcpp Vector
Rcpp::NumericMatrix o(n,2);
// Create a prng engine with a specific seed
sitmo::prng eng(static_cast<uint32_t>(seed));
// Draw from base engine
for (unsigned int i=0; i < n ; ++i){
o(i,0) = eng();
}
// Reset seed
eng.seed();
// Draw from base engine
for (unsigned int i=0; i< n ; ++i){
o(i,1) = eng();
}
return o;
}
```

This example displays the ability of `sitmo`

to handle parallel streams of rng when a predefined number of streams is known.

```
#include <Rcpp.h>
#include <sitmo.h> // SITMO PPRNG
// [[Rcpp::depends(sitmo)]]
//' Two RNG engines running side-by-side
//'
//' Shows how to create two separate RNGs and increase them together.
//'
//' @param n An \code{unsigned int} that dictates how many realizations occur.
//' @param seeds A \code{vec} containing two integers greater than 0.
//' @return A \code{matrix} with random sequences.
//' @examples
//' n = 10
//' a = sitmo_two_seeds(n, c(1337,1338))
//'
//' b = sitmo_two_seeds(n, c(1337,1337))
//'
//' isTRUE(all.equal(a[,1],a[,2]))
//'
//' isTRUE(all.equal(b[,1],b[,2]))
//'
//' isTRUE(all.equal(a[,1],b[,1]))
// [[Rcpp::export]]
Rcpp::NumericMatrix sitmo_two_seeds(unsigned int n, Rcpp::NumericVector seeds) {
if(seeds.size() != 2) Rcpp::stop("Need exactly two seeds for this example.");
// Create Rcpp Matrix
Rcpp::NumericMatrix o(n,2);
// Create a prng engine with a specific seed
sitmo::prng eng1;
eng1.seed(seeds(0));
sitmo::prng eng2;
eng2.seed(seeds(1));
// Draw from base engine
for (unsigned int i=0; i< n ; ++i){
o(i,0) = eng1();
o(i,1) = eng2();
}
return o;
}
```

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”).

```
#include <Rcpp.h>
#include <sitmo.h> // SITMO PPRNG
// [[Rcpp::depends(sitmo)]]
//' Random Uniform Number Generator with sitmo
//'
//' The function provides an implementation of sampling from a random uniform distribution
//'
//' @param n An \code{unsigned integer} denoting the number of realizations to generate.
//' @param min A \code{double} indicating the minimum \eqn{a} value
//' in the uniform's interval \eqn{\left[a,b\right]}
//' @param max A \code{double} indicating the maximum \eqn{b} value
//' in the uniform's interval \eqn{\left[a,b\right]}
//' @param seed A special \code{unsigned integer} containing a single seed.
//' @return A \code{vec} containing the realizations.
//' @export
//' @examples
//' a = runif_sitmo(10)
// [[Rcpp::export]]
Rcpp::NumericVector runif_sitmo(unsigned int n, double min = 0.0, double max = 1.0, uint32_t seed = 1) {
Rcpp::NumericVector o(n);
// Create a prng engine
sitmo::prng eng(seed);
// Obtain the range between max and min
double dis = max - min;
for(int i = 0; i < n; ++i) {
// Sample from the RNG and divide it by the maximum value possible (can also use SITMO_RAND_MAX, which is 4294967295)
// Apply appropriate scale (MAX-MIN)
o[i] = min + ((double) eng() / (sitmo::prng::max())) * (dis);
}
return o;
}
```

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:

```
#ifdef _OPENMP
// multithreaded OpenMP version of code
#else
// single-threaded version of code
#endif
```

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;
}
```