---
title: "Parallel RNG usage"
author: "Ralf Stubner"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Parallel RNG usage}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

When you want to use random number generators (RNG) for parallel computations, you need to make sure that the sequences of random numbers used by the different processes do not overlap.
There are two main approaches to this problem:^[See for example http://www.pcg-random.org/posts/critiquing-pcg-streams.html.]

* Partition the complete sequence of random numbers produced for one seed into non-overlapping sequences and assign each process one sub-sequence.
* Re-parametrize the generator to produce independent sequences for the same seed.

The RNGs included in `dqrng` offer at least one of these methods for parallel RNG usage.
When using the R or C++ interface independent streams can be accessed using the two argument `dqset.seed(seed, stream)` or `dqset_seed(seed, stream)` functions.

# Threefry: usage from R

The Threefry engine uses internally a counter with $2^{256}$ possible states, which can be split into different substreams.
When used from R or C++ with the two argument `dqset.seed` or `dqset_seed` this counter space is split into $2^{64}$ streams with $2^{192}$ possible states each.
This is equivalent to $2^{64}$ streams with a period of $2^{194}$ each.

In the following example a matrix with random numbers is generated in parallel using the parallel package.
The resulting correlation matrix should be close to the identity matrix if the different streams are independent:

```{r}
library(parallel)
cl <- parallel::makeCluster(2)
res <- clusterApply(cl, 1:8, function(stream, seed, N) {
  library(dqrng)
  dqRNGkind("Threefry")
  dqset.seed(seed, stream)
  dqrnorm(N)
}, 42, 1e6)
stopCluster(cl)

res <- matrix(unlist(res), ncol = 8)
symnum(x = cor(res), cutpoints = c(0.001, 0.003, 0.999),
       symbols = c(" ", "?", "!", "1"),
       abbr.colnames = FALSE, corr = TRUE)
```

As expected the correlation matrix for the different columns is almost equal to the identity matrix.


# Xo(ro)shiro: jump ahead with OpenMP

The Xoshiro256+ generator has a period of $2^{256} -1$ and offers $2^{128}$ sub-sequences that are $2^{128}$ random draws apart as well as $2^{64}$ streams that are $2^{192}$ random draws appart.
The Xoroshiro128+ generator has a period of $2^{128} -1$ and offers $2^{64}$ sub-sequences that are $2^{64}$ random draws apart as well as $2^{32}$ streams that are $2^{98}$ random draws appart. 
You can go from one sub-sequence to the next using the `jump()` or `long_jump()` method and the convenience wrapper `jump(int n)` or `long_jump(int n)`, which  advances to the `n`th sub-sequence.
When used from R or C++ with the two argument `dqset.seed` and `dqset_seed` you get $2^{64}$ streams that are $2^{192}$ and $2^{64}$ random draws appart for Xoshiro256+ and Xoroshiro128+, respectively.

As an example using C++ we draw and sum a large number of uniformly distributed numbers.
This is done several times using OpenMP for parallelisation.
Care has been taken to keep the global RNG `rng` usable outside of the parallel block.

```{r, eval=FALSE, engine='Rcpp'}
#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH, sitmo)]]
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

// [[Rcpp::export]]
std::vector<double> parallel_random_sum(int n, int m, int ncores) {
  dqrng::uniform_distribution dist(0.0, 1.0); // Uniform distribution [0,1)
  dqrng::xoshiro256plus rng(42);              // properly seeded rng
  std::vector<double> res(m);
  // ok to use rng here
  
  #pragma omp parallel num_threads(ncores)
  {
    dqrng::xoshiro256plus lrng(rng);      // make thread local copy of rng 
    lrng.long_jump(omp_get_thread_num() + 1);  // advance rng by 1 ... ncores jumps

    #pragma omp for
    for (int i = 0; i < m; ++i) {
      double lres(0);
      for (int j = 0; j < n; ++j) {
        lres += dist(lrng);
      }
      res[i] = lres / n;
    }
  }
  // ok to use rng here
  return res;
}

/*** R
parallel_random_sum(1e7, 8, 4)
*/
```


Result:

    [1] 0.5001591 0.5000428 0.4999855 0.4999706 0.5000061 0.4999447 0.4999188 0.5001192


# PCG: multiple streams with RcppParallel

From the PCG family we will look at pcg64, a 64-bit generator with a period of $2^{128}$.
It offers the function [`advance(int n)`](http://www.pcg-random.org/using-pcg-cpp.html#void-advance-state-type-delta), which is equivalent to `n` random draws but scales as $O(ln(n))$ instead of $O(n)$. 
In addition, it offers $2^{127}$ separate streams that can be enabled via the function [`set_stream(int n)`](http://www.pcg-random.org/using-pcg-cpp.html#void-pcg32-set-stream-state-type-stream) or the [two argument constructor](http://www.pcg-random.org/using-pcg-cpp.html#pcg32-pcg32-constructor) with `seed` and `stream`.
When used from R or C++ with the two argument `dqset.seed` and `dqset_seed` you get $2^{64}$ streams out of the possible $2^{127}$ separate streams.

In the following example a matrix with random numbers is generated in parallel using RcppParallel.
The resulting correlation matrix should be close to the identity matrix if the different streams are independent:


```{r, eval=FALSE, engine='Rcpp'}
#include <Rcpp.h>
// [[Rcpp::depends(dqrng, BH, sitmo)]]
#include <pcg_random.hpp>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel)]]
#include <RcppParallel.h>

struct RandomFill : public RcppParallel::Worker {
  RcppParallel::RMatrix<double> output;
  uint64_t seed;
  dqrng::normal_distribution dist{0.0, 1.0};

  RandomFill(Rcpp::NumericMatrix output, const uint64_t seed) : output(output), seed(seed) {};

  void operator()(std::size_t begin, std::size_t end) {
    pcg64 rng(seed, end);
    auto gen = std::bind(dist, std::ref(rng));
    for (std::size_t col = begin; col < end; ++col) {
      RcppParallel::RMatrix<double>::Column column = output.column(col);
      std::generate(column.begin(), column.end(), std::ref(gen));
    }
  }
};

// [[Rcpp::export]]
Rcpp::NumericMatrix parallel_random_matrix(const int n, const int m, const int ncores) {
  Rcpp::NumericMatrix res(n, m);
  RandomFill randomFill(res, 42);
  RcppParallel::parallelFor(0, m, randomFill, m/ncores + 1);
  return res;
}

/*** R
res <- parallel_random_matrix(1e6, 8, 4)
head(res)
symnum(x = cor(res), cutpoints = c(0.001, 0.003, 0.999),
       symbols = c(" ", "?", "!", "1"),
       abbr.colnames = FALSE, corr = TRUE)
*/
```

Head of the random matrix:

               [,1]       [,2]        [,3]       [,4]       [,5]       [,6]       [,7]       [,8]
    [1,]  0.7114429  1.2642926 -0.47149983  0.4277034 -0.3709571  0.4336255  0.8185977  0.3505209
    [2,]  0.8721661 -1.1908067  1.10411136  0.4160170 -1.3276661 -0.4182534 -1.2437521  1.0084814
    [3,] -1.4959624 -0.1726986 -0.54343828 -0.5635330 -1.1779352  0.7539401 -0.4341252 -1.2256560
    [4,]  0.5087201  0.1116120  0.19007581 -0.8220532  0.9672267 -1.1246750 -0.5283977 -1.3198088
    [5,] -0.8191448 -1.1856318 -0.03458304  0.8027613  1.0594094 -0.4456480 -0.5456669  0.2593546
    [6,]  1.2289518 -0.2576456  0.40222707  0.9757078 -0.5796549 -1.0148114  2.8961294  0.6391352

Correlation matrix:

    [1,] 1              
    [2,]   1            
    [3,] ?   1          
    [4,] ? ? ? 1        
    [5,]         1      
    [6,]       ?   1    
    [7,]         ? ? 1  
    [8,]   ?     ?     1
    attr(,"legend")
    [1] 0 ‘ ’ 0.001 ‘?’ 0.003 ‘!’ 0.999 ‘1’ 1

So as expected the correlation matrix is almost equal to the identity matrix.
