Make vignettes buildable

.github/workflows/R-CMD-check.yaml

@@ -70,6 +70,11 @@ jobs:
           remotes::install_cran("rcmdcheck")
         shell: Rscript {0}
 
+      - name: Move real vignettes
+        run: |
+          cp vignettes_src/rng_package.Rmd vignettes
+          cp vignettes/rng_pi*.cpp vignettes
+
       - name: Check
         env:
           _R_CHECK_CRAN_INCOMING_REMOTE_: false

.gitignore

@@ -16,7 +16,7 @@ inst/doc
 pkgdown
 inst/include/cub
 *.gcov
-vignettes_src/gpu.md
+vignettes_src/*.md
 
 .vscode/
 

Makefile

@@ -49,7 +49,10 @@ clean:
 		src/*.gcov src/*.gcda src/*.gcno
 
 vignettes/gpu.Rmd: vignettes_src/gpu.Rmd
-	./scripts/build_gpu_vignette
+	./scripts/build_vignette gpu
+
+vignettes/rng_package.Rmd: vignettes_src/rng_package.Rmd
+	./scripts/build_vignette rng_package
 
 vignettes: vignettes/dust.Rmd vignettes/rng.Rmd
 	${RSCRIPT} -e 'tools::buildVignettes(dir = ".")'

scripts/build_vignette

@@ -0,0 +1,10 @@
+#!/usr/bin/env bash
+set -e
+SRC="$1.Rmd"
+DEST="$1.md"
+
+(cd vignettes_src && Rscript -e "knitr::knit(\"${SRC}\")")
+header="DO NOT EDIT THIS FILE - see vignettes_src and make changes there"
+sed -s 's/[[:space:]]*$//' vignettes_src/$DEST |
+    sed 's/\r//g' |
+    sed "s/HEADER/$header/" > vignettes/$SRC

vignettes/rng_package.Rmd

@@ -0,0 +1,226 @@
+---
+title: "Using RNGs from packages"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Using RNGs from packages}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+<!-- DO NOT EDIT THIS FILE - see vignettes_src and make changes there -->
+
+
+
+The dust random number generators are suitable for use from other packages and we provide a few helpers in both R and C++ to make this easier.
+
+For illustration purposes we will assume that you want to estimate pi using rejection sampling.  The basic idea of this algorithm is simple; sample two U(0, 1) points `x` and `y` and test if they lie in the unit circle (i.e. `sqrt(x^2 + x^2) < 1`) giving the ratio of the area of a unit circle to a square. Multiplying the fraction of the points that we accept by 4 yields our estimate of pi.
+
+## Background using R's random number generator
+
+First, an example that uses R's API for example (note that while the R API is C, we're using the cpp11 package interface here so that the following examples are similar):
+
+```cc
+#include <cpp11.hpp>
+#include <R_ext/Random.h>
+
+[[cpp11::register]]
+double pi_r(int n) {
+  int tot = 0;
+  GetRNGstate();
+  for (int i = 0; i < n; ++i) {
+    const double u1 = unif_rand();
+    const double u2 = unif_rand();
+    if (u1 * u1 + u2 * u2 < 1) {
+      tot++;
+    }
+  }
+  PutRNGstate();
+  return tot / static_cast<double>(n) * 4.0;
+}
+```
+
+With cpp11 we can load this with `cpp11::cpp_source`
+
+
+```r
+cpp11::cpp_source("rng_pi_r.cpp")
+```
+
+and then run it wih
+
+
+```r
+pi_r(1e6)
+#> [1] 3.143772
+```
+
+The key bits within the code above are that we:
+
+* included random number support from R via the `R_ext/Random.h` header
+* initialised the state of the global random number stream within our C++ context with `GetRNGstate` before drawing any random numbers
+* drew some possible large number of numbers using `unif_rand`
+* restored the random number state back to R.
+
+Failure to run the `GetRNGstate` / `PutRNGstate` will result in the stream not behaving properly.  This is explained in detail in the "Writing R Extensions" manual.
+
+## Basic implementation using dust
+
+One of the design ideas in dust is that there is no single global source of random numbers, so we need to create a source
+
+
+```r
+rng <- dust:::dust_rng_pointer$new(seed = 42)
+rng
+#> <dust_rng_pointer>
+#>   Public:
+#>     algorithm: xoshiro256plus
+#>     initialize: function (seed = NULL, n_streams = 1L, algorithm = "xoshiro256plus")
+#>     is_current: function ()
+#>     n_streams: 1
+#>     state: function ()
+#>     sync: function ()
+#>   Private:
+#>     is_current_: TRUE
+#>     ptr_: externalptr
+#>     state_: 95 6e eb 2f 26 32 d7 bd 04 72 10 65 ba fa e1 57 46 7f 20 ...
+```
+
+This object acts as a "handle" to some random number state that can be passed safely to C++ programs; the state will be updated when the program runs as a side effect.  Unlike the `dust::dust_rng` object there are no real useful methods on this object and from the R side we'll treat it as a black box. Importantly the `rng` object knows which algorithm it has been created to use
+
+
+```r
+rng$algorithm
+#> [1] "xoshiro256plus"
+```
+
+The default will be suitable for most purposes.
+
+We can rewrite the pi approximation function as:
+
+```cc
+#include <cpp11.hpp>
+#include <dust/r/random.hpp>
+
+[[cpp11::linking_to(dust)]]
+[[cpp11::register]]
+double pi_dust(int n, cpp11::sexp ptr) {
+  auto rng =
+    dust::random::r::rng_pointer_get<dust::random::xoshiro256plus_state>(ptr);
+  auto& state = rng->state(0);
+  int tot = 0;
+  for (int i = 0; i < n; ++i) {
+    const double u1 = dust::random::random_real<double>(state);
+    const double u2 = dust::random::random_real<double>(state);
+    if (u1 * u1 + u2 * u2 < 1) {
+      tot++;
+    }
+  }
+  return tot / static_cast<double>(n) * 4.0;
+}
+```
+
+This snippet looks much the same as above:
+
+* We've added `[[cpp::linking_to(dust)]]` and included the dust random interface (`dust/r/random.hpp`)
+* The first line of the function safely creates a pointer to the random state data. The template argument here (`<dust::random::xoshiro256plus_state>`) refers to the rng algorithm and matches `rng$algorithm`
+* The second line extracts a reference to the first (C++ indexing starting at 0) random number stream - this pair of lines is roughly equivalent to `GetRNGstate()` except that that the random numbers do not come from some global source
+* After that the listing proceeds as before proceeds as before, except there is no equivalent to `PutRNGstate()`
+
+
+
+```r
+cpp11::cpp_source("rng_pi_dust.cpp")
+```
+
+and then run it wih
+
+
+```r
+pi_dust(1e6, rng)
+#> [1] 3.14098
+```
+
+## Parallel implementation with dust and OpenMP
+
+Part of the point of dust's random number generators is that they create independent streams of random numbers that can be safely used in parallel.
+
+```cc
+#include <cpp11.hpp>
+#include <dust/r/random.hpp>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+[[cpp11::linking_to(dust)]]
+[[cpp11::register]]
+double pi_dust_parallel(int n, cpp11::sexp ptr, int n_threads) {
+  auto rng =
+    dust::random::r::rng_pointer_get<dust::random::xoshiro256plus_state>(ptr);
+  const auto n_streams = rng->size();
+  int tot = 0;
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static)  num_threads(n_threads) \
+  reduction(+:tot)
+#endif
+  for (size_t i = 0; i < n_streams; ++i) {
+    auto& state = rng->state(0);
+    int tot_i = 0;
+    for (int i = 0; i < n; ++i) {
+      const double u1 = dust::random::random_real<double>(state);
+      const double u2 = dust::random::random_real<double>(state);
+      if (u1 * u1 + u2 * u2 < 1) {
+        tot_i++;
+      }
+    }
+    tot += tot_i;
+  }
+  return tot / static_cast<double>(n * n_streams) * 4.0;
+}
+```
+
+
+```r
+cpp11::cpp_source("rng_pi_parallel.cpp")
+```
+
+Here we've made a number of decisions about how to split the problem up subject to a few constraints about using OpenMP together with R:
+
+* Generally speaking we want the answer to be independent of the number of threads used to run it, as this will vary in different sessions. As such avoid the temptation to do a loop over the threads; here we instead iterate over *streams* with the idea that there will be one or more streams used per threads. If we ran a single thread we'd get the same answer as if we ran one thread per stream.
+* Each thread is going to do it's own loop of length `n` so we need to divide by `n * n_streams` at the end as that's many attempts we have made.
+* We use OpenMP's `reduction` clause to safely accumulate the different subtotals (the `tot_i` values) into one `tot` value.
+* In order to compile gracefully on machines that do not have OpenMP support both the `#include <omp.h>` line and the `#pragma omp` line are wrapped in guards that test for `_OPENMP` (see "Writing R Extensions").
+* We let the generator tell us how many streams it has (`n_streams = rng->size()`) but we could as easily specify an ideal number of streams as an argument here and then test that the generator has *at least that many* by adding an argument to the call to `rng_pointer_get` (e.g., if we wanted `m` streams the call would be `rng_pointer_get<type>(ptr, m)`)
+
+
+```r
+rng <- dust:::dust_rng_pointer$new(seed = 42, n_streams = 20)
+pi_dust_parallel(1e6, rng, 4)
+#> [1] 3.141703
+```
+
+Unfortunately [`cpp11::cpp_source` does not support using OpenMP](https://github.com/r-lib/cpp11/issues/243) so in the example above the code will run in serial and we can't see if parallelisation will help.
+
+In order to compile with support, we need to build a little package and set up an appropriate `Makevars` file
+
+
+
+Once we have a parallel version we can see a speed-up as we add threads:
+
+
+```r
+rng <- dust:::dust_rng_pointer$new(n_streams = 20)
+bench::mark(
+  pi_dust_parallel(1e6, rng, 1),
+  pi_dust_parallel(1e6, rng, 2),
+  pi_dust_parallel(1e6, rng, 3),
+  pi_dust_parallel(1e6, rng, 4),
+  check = FALSE)
+#> # A tibble: 4 x 6
+#>   expression                           min   median `itr/sec` mem_alloc `gc/sec`
+#>   <bch:expr>                      <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
+#> 1 pi_dust_parallel(1e+06, rng, 1)   44.3ms   44.7ms      22.4        0B        0
+#> 2 pi_dust_parallel(1e+06, rng, 2)   22.5ms     23ms      43.4        0B        0
+#> 3 pi_dust_parallel(1e+06, rng, 3)   15.8ms   15.9ms      62.7        0B        0
+#> 4 pi_dust_parallel(1e+06, rng, 4)   11.5ms   11.6ms      85.5        0B        0
+```

vignettes_src/rng_use.Rmd → vignettes_src/rng_package.Rmd

@@ -61,7 +61,7 @@ Failure to run the `GetRNGstate` / `PutRNGstate` will result in the stream not b
 One of the design ideas in dust is that there is no single global source of random numbers, so we need to create a source
 
 ```{r}
-rng <- dust:::dust_rng_pointer$new()
+rng <- dust:::dust_rng_pointer$new(seed = 42)
 rng
 ```
 
@@ -118,7 +118,7 @@ Here we've made a number of decisions about how to split the problem up subject
 * We let the generator tell us how many streams it has (`n_streams = rng->size()`) but we could as easily specify an ideal number of streams as an argument here and then test that the generator has *at least that many* by adding an argument to the call to `rng_pointer_get` (e.g., if we wanted `m` streams the call would be `rng_pointer_get<type>(ptr, m)`)
 
 ```{r}
-rng <- dust:::dust_rng_pointer$new(n_streams = 20)
+rng <- dust:::dust_rng_pointer$new(seed = 42, n_streams = 20)
 pi_dust_parallel(1e6, rng, 4)
 ```