---
title: "Random number generation"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Random number generation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  error = FALSE,
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5)
lang_output <- function(x, lang) {
  cat(c(sprintf("```%s", lang), x, "```"), sep = "\n")
}
cc_output <- function(x) lang_output(x, "cc")
r_output <- function(x) lang_output(x, "r")
plain_output <- function(x) lang_output(x, "plain")
```

The core of `dust` is a very fast parallel random number generation system.  This vignette describes the details of the system and why it is needed.

We provide an interface to the ["Xoshiro" / Blackman-Vigna](https://prng.di.unimi.it/) family of generators (Xoshiro is derived from XOR/shift/rotate). These are designed to allow use in parallel by "jumping ahead" in the sequence and we use this below to interleave generators. The stream is unrelated to and unaffected by R's random number generation. `set.seed` has no effect, for example. **The random numbers are not cryptographically secure**; for that see the excellent [sodium](https://cran.r-project.org/package=sodium) package.

Ordinarily this is used from C++; the model as discussed in `vignette("dust")` uses `rng_state_type` and a function from `dust::random` to interact with the generator. However, an R interface is provided for debugging and testing purposes.

```{r}
rng <- dust::dust_rng$new(seed = 1)
rng
```

Currently only a few distributions are supported, with an interface that somewhat mimics R's interface (we avoid abbreviations so that `rnorm` becomes `normal` and there are no default arguments).  For example, with the generator above we can generate 100 standard normal samples:

```{r}
rng$normal(100, 0, 1)
```

One feature that we use here is to allow multiple streams of random numbers within a rng object. If running in parallel (either via a dust model or by passing `n_threads`) these different random number streams can be given to different threads.

```{r}
rng1 <- dust::dust_rng$new(seed = 1, n_streams = 1)
rng2 <- dust::dust_rng$new(seed = 1, n_streams = 2)
rng1$random_real(5)
rng2$random_real(5)
```

Notice here how in the output from `rng2`, the first column matches the series of numbers out of `rng1`.

This is achieved by "jumping" the random number streams forward.  Here are the second column of numbers from the output of `rng2`:

```{r}
rng3 <- dust::dust_rng$new(seed = 1)$jump()
rng3$random_real(5)
```

With the default generator (`xoshoro256+`, see below), a jump is equivalent to 2^128 draws from the random number generator (about 10^38). There are 2^128 of these non-overlapping subsequences in the generator, which is quite a lot. If this feels too close together, then the `$long_jump()` method jumps even further (2^192 draws, or about 10^57). There are 2^64 (10^20) of these sequences.

## Supported distributions

We do not yet support the full set of distributions provided by R, and the distributions that we do support represent our immediate needs.  Contributions of further distributions is welcome. Note that in the R version the first argument represents the number of draws requested, but in the C++ interface we always return a single number.

**Uniform distribution** between `min` and `max`:

```{r}
rng$uniform(10, 3, 6)
```

**Normal distribution** with parameters `mean` and `sd`

```{r}
rng$normal(10, 3, 6)
```

This function supports using the ziggurat algorithm, which will generally be faster:

```{r}
rng$normal(10, 3, 6, algorithm = "ziggurat")
```

**Poisson distribution** with mean `lambda`

```{r}
rng$poisson(10, 4.5)
```

**Binomial distribution** with mean `size` and `prob`

```{r}
rng$binomial(10L, 10, 0.3)
```

**Negative binomial distribution** with mean `size` and `prob`

```{r}
rng$nbinomial(10L, 10, 0.3)
```

**Hypergeometric distribution** with parameters `n1` (number of white balls), `n2` (number of black balls) and `k` (number of samples)

```{r}
rng$hypergeometric(10L, 3, 10, 4)
```

**Exponential distribution** with rate `rate`


```{r}
rng$exponential(10L, 4)
```


**Multinomial distribution** with parameters `size` and `prob`

```{r}
rng$multinomial(10L, 20, c(0.3, 0.5, 0.2))
```

There are also two special case distributions, the **Standard uniform distribution** (between 0 and 1) - faster than using `$uniform(n, 0, 1)`:

```{r}
rng$random_real(10)
```

```cc
real_type u = dust::random::random_real<real_type>(state);
```

the **Standard normal distribution** (between mean 0 and sd 1) - faster than using `$normal(n, 0, 1)`:

```{r}
rng$random_normal(10)
```

```cc
real_type u = dust::random::random_normal<real_type>(state);
```

**Gamma distribution** parameterized by `shape` and `scale`
```{r}
rng$gamma(10L, 2, 5)
```

All of these distributions are available in C++, and this is documented [in the C++ reference documentation](https://mrc-ide.github.io/dust/cpp/random-distributions.html)

## Performance

Performance should be on par or better with R's random number generator, though here the timings are likely to be mostly due to allocations and copies of memory:

```{r}
bench::mark(
  rng1$random_real(1000),
  rng1$uniform(1000, 0, 1),
  runif(1000),
  time_unit = "us",
  check = FALSE)
```

For normally distributed random numbers, the different algorithms will perform differently:

```{r}
bench::mark(
  rng1$random_normal(1000),
  rng1$random_normal(1000, algorithm = "ziggurat"),
  rng1$normal(1000, 0, 1),
  rnorm(1000),
  time_unit = "us",
  check = FALSE)
```

On reasonable hardware (10-core i9 at 2.8 GHz) we see throughput of ~1 billion U(0, 1) draws per second (1ns/draw) scaling linearly to 10 cores if the numbers are immediately discarded.  Doing anything with the numbers or trying to store them comes with a non-trivial overhead.

The difference between `random_real` and `uniform` here is the cost of recycling the parameters, not the actual generation!

Binomial distribution, small `n * p`, which uses an inversion algorithm

```{r}
rng1 <- dust::dust_rng$new(seed = 1)
n <- as.numeric(rep(9:10, length.out = 1000))
p <- rep(c(0.1, 0.11), length.out = 1000)
bench::mark(
  rng1$binomial(1000, n, p),
  rbinom(1000, n, p),
  time_unit = "us",
  check = FALSE)
```

(note we vary `n` and `p` here as we've optimised things for random parameter access).

Large `n * p` uses a rejection sampling algorithm

```{r}
n <- as.numeric(rep(9999:10000, length.out = 1000))
p <- rep(c(0.3, 0.31), length.out = 1000)
bench::mark(
  rng1$binomial(1000, n, p),
  rbinom(1000, n, p),
  time_unit = "us",
  check = FALSE)
```

Practically, performance will vary based on the parameters of the distributions and the underlying algorithms, and the principle performance gain we hope to get here is driven by the ability to parallelise safely rather than the speed of the draws themselves.

## Underlying random number engine

Under the hood, the random number generators work by first creating a random integer, which we then convert to a floating point number, then for the distributions above we apply algorithms that convert one or more uniformly distributed (i.e., U(0, 1)) floating point numbers to a given distribution through techniques such as inversion (for a random exponential) or rejection sampling (for a random binomial).

```
[random integer] -> [random real] -> [random draw from a distribution]
```

We include 12 different algorithms for creating the underlying random integer, from the [same family of generators](https://prng.di.unimi.it/).  These provide different underlying storage types (either 32 bit or 64 bit integers) and different state types.

Normally you do not need to worry about these details and declaring

```cc
using rng_state_type = dust::random::generator<real_type>;
```

in your model will select a reasonable generator.

If for some reason you want to try a different generator you can directly specify one of the types, for example

```cc
using rng_state_type = dust::random::xoshiro256starstar;
```

which means that whatever real type you use, you want to use the Xoshiro256** generator.

The supported types are:

* `xoshiro256starstar`, `xoshiro256plusplus`, `xoshiro256plus` (4 x 64 bit state)
* `xoroshiro128starstar`, `xoroshiro128plusplus`, `xoroshiro128plus` (2 x 64 bit state)
* `xoshiro128starstar`, `xoshiro128plusplus`, `xoshiro128plus` (4 x 32 bit state)
* `xoshiro512starstar`, `xoshiro512plusplus`, `xoshiro512plus` (8 x 64 bit state; this is far more state space than typically needed)

The "starstar", "plusplus" or "plus" refers to the final scrambling operation (two multiplications, two additions or one addition, respectively); the speeds of these might vary depending on your platform.  The "plus" version will be the fastest but produces slightly less randomness in the the lower bits of the underlying integers, which theoretically is not a problem when converting to a real number.

If you are generating single precision `float` numbers for a GPU, then you may want to use the `xoshiro128` family as these will be faster than the 64 bit generators on that platform.  On a CPU you will likely not see a difference

```{r}
rng_f <- dust::dust_rng$new(seed = 1, real_type = "float")
rng_d <- dust::dust_rng$new(seed = 1, real_type = "double")
bench::mark(
  rng_f$random_real(1000),
  rng_d$random_real(1000),
  time_unit = "us",
  check = FALSE)
```

If you do not need to run many parallel realisations, then the `xoroshiro256` (note the extra `ro`)  generators may be preferable as they carry half the state and may be slightly faster.  The `xoshiro512` generators are included for completeness and offer an excessive state space).

The table below summarises the properties of the underlying generators (each of these exists with three different scrambling options).

Name          | Size    | State    | Period | Jump   | Long jump
--------------+---------+----------+--------+--------+-----------
`xoshiro128`  | 32 bits | 2 uint32 | 2^128  | 2^64   | 2^96
`xoroshiro128`| 64 bits | 2 uint64 | 2^128  | 2^64   | 2^96
`xoshiro256`  | 64 bits | 4 uint64 | 2^256  | 2^128  | 2^192
`xoshiro512`  | 64 bits | 4 uint64 | 2^512  | 2^512  | 2^384

Size is the size of the returned random integer (32 or 64 bits, though with the `+` scrambler not all bits will be of high quality). The State is the number and size of the internal state of the generator.  Period refers to the number of states that the generator will move through, Jump and Long jump are the number of steps (equivalent) that a jump operation will take through this sequence.

Note that the period is the two size of number of bits in the model state (e.g., xoshiro256 uses 4 x 64 bit integers, for a total of 256 bits of state, 2^256 possible combinations, each of which will be visited once in the cycle).  The Jump coefficients have been computed to have size `sqrt(Period)`.

## Reusing the random random number generator in other projects

Our random number library can be reused in other projects without using anything else from dust; either in an R package or in a standalone project.

### In a package

A minimal package looks like

```{r pkg_tree, echo = FALSE}
path_pkg <- dust:::dust_file("random/package")
withr::with_dir(path_pkg, fs::dir_tree())
```

with the core of the C++ file containing a small program that uses the dust random number generator to draw a series of normally distributed random number with a single mean and standard deviation.

```{r echo = FALSE, results = "asis"}
code <- readLines(file.path(path_pkg, "src/rnguse.cpp"))
cc_output(code[!grepl("^//", code)])
```

To complete the package, the `DESCRIPTION` includes `dust` and `cpp11` in the `LinkingTo` section:

```{r echo = FALSE, results = "asis"}
plain_output(readLines(file.path(path_pkg, "DESCRIPTION")))
```

You must remember to update your `NAMESPACE` to include `useDynLib` (either directly or via roxygen)

```{r echo = FALSE, results = "asis"}
plain_output(readLines(file.path(path_pkg, "NAMESPACE")))
```

Finally, run `cpp11::cpp_register()` before compiling your package so that the relevant interfaces are created (`R/cpp11.R` and `cpp11/cpp11.cpp`).  A similar process would likely work with Rcpp without any dependency on cpp11.

The issue with this example is that every call to `random_normal` must provide its own `seed` argument, so the random number streams are not continued with each call to the function. This is is not very useful in practice and we describe more fully how to do this properly in `vignette("rng_package.Rmd")`.

### Standalone, parallel with OpenMP

*This is somewhat more experimental, so let us know if you have success using the library this way.*

It is possible to include the dust random library in a standalone C++ program (or one embedded from another language) without using any R support.  The library is a header only library and `<dust/random/random.hpp>` is the main file to include.

The simplest way to get started is to copy the contents of `inst/include/dust/random` into your project.  You can do this by downloading and unpacking [the latest release of dust-random](https://github.com/mrc-ide/dust/releases/latest/download/dust-random.tar.gz).  Then after including `<dust/random/random.hpp>` you can use the contents of the random number library.

For example, below is a small program that computes the sum of a series of random numbers, in parallel using OpenMP to parallelise across a set of generators.

```{r echo = FALSE, results = "asis"}
cc_output(readLines(dust:::dust_file("random/openmp/rnguse.cpp")))
```

This program can be compiled with

```
g++ -I$(PATH_DUST_INCLUDE) -O2 -std=c++11 -fopenmp -o rnguse rnguse.cpp
```

where `$PATH_DUST_INCLUDE` is the path to the header only library.

### Standalone, parallel on a GPU

This is considerably more complicated and will depend on your aims with your GPU-accelerated program.

We include [examples in a repository showing `dust::random` vs `curand` benchmarks](https://github.com/mrc-ide/dust-random-bench/blob/main/dustrand.cu). This covers setting up the RNG state so that it can be set from the host and retrieved, interleaving as appropriate, plus samples from the uniform, normal, exponential Poisson and binomial distributions. We overlap with `curand` only for uniform, normal, and Poisson.  See [the repository](https://github.com/mrc-ide/dust-random-bench) for more details.

## Other packages with similar functionality

There are many packages offering similar functionality to `dust`:

* [`sitmo`](https://cran.r-project.org/package=sitmo) offers several modern RNGs (sitmo, threefry, vandercorput) for generating standard uniform numbers (U(0, 1)), designed to be usable from an Rcpp-using package.
* [`dqrng`](https://cran.r-project.org/package=dqrng) offers an overlapping set and generators for normal and exponential random numbers, also designed to be used with Rcpp and BH.
* [`rTRNG`](https://cran.r-project.org/package=rTRNG) exposes "Tina’s Random Number Generator", among others, and support for several distributions. It can be used from other packages but requires linking (i.e., not header only) which does not work on macOS. Generator jumps typically require a known number of draws per thread.
* [`miranda`](https://coolbutuseless.github.io/2020/07/09/introducing-miranda-a-package-of-fast-modern-uniform-pseudo-random-number-generators/) includes support for several underlying generators (including `xoshiro256+`) but with a single global state and no support for use from other packages' compiled code.

Why does this package exist?  We had a set of needs that meant using the a

* generation of single precision `float` numbers with the same interface
* absolutely no global state
* nice interface for reproducible parallel random number generation where the generated number does not depend on the number of threads used
* possible to use on GPUs; this exposes some exotic issues around what types of structures can be used
* minimal dependencies and impact on the C++ code (e.g., no boost/BH) and header-only in order to satisfy the above
* support for additional distributions, especially binomial, in a thread-safe way, optimised for rapidly changing parameters

On this last point; we noticed that many distribution packages, including Boost.Random (and to a degree R's random functions) assume that you want to sample many random numbers from a distribution with fixed parameters.  This is probably reasonable in many cases, but in the use case we had (stochastic simulation models) we expected that all parameters were likely to change at every iteration.  These different approaches have important tradeoffs - if you will take many samples from a single distribution you might compute tables of coefficients once and use them many times which will make the sampling faster, but this will be wasteful if only a single sample is taken before the parameters change.

The situation is slightly more complex on a GPU where we need to make sure that different threads within a block do not get needlessly onto different branches of conditional logic.  Things like early exits need to be avoided and we have carefully profiled to make sure that threads within a warp end up re-synchronised at the optimal spot (see for example [this blog post on `__syncwarp`](https://developer.nvidia.com/blog/using-cuda-warp-level-primitives/)).
