Feature/parallel generate (#82)

* Add header file for parallel variate generation
* Make clone(0) work also for PCG
* Move seeding from R into dqrng::generator<RNG>()
* Rework Parallel RNG vignette
    * use dqrng::generator to create the (main) RNG and rng->clone() for the local versions
    * remove the first two examples that make use of the global RNG
    * add a section on global RNG usage which discusses dqrng_extra/parallel_generate.h

Author: rstub

DESCRIPTION

@@ -1,16 +1,17 @@
 Package: dqrng
 Type: Package
 Title: Fast Pseudo Random Number Generators
-Version: 0.3.2.5
+Version: 0.3.2.6
 Authors@R: c(
-    person("Ralf", "Stubner", email = "ralf.stubner@gmail.com", role = c("aut", "cre")),
+    person("Ralf", "Stubner", email = "ralf.stubner@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0009-0009-1908-106X")),
     person("daqana GmbH", role = "cph"),
     person("David Blackman", role = "cph", comment = "Xoroshiro / Xoshiro family"),
     person("Melissa O'Neill", email = "oneill@pcg-random.org", role = "cph", comment = "PCG family"),
     person("Sebastiano Vigna", email = "vigna@acm.org", role = "cph", comment = "Xoroshiro / Xoshiro family"),
     person("Aaron", "Lun", role="ctb"), 
     person("Kyle", "Butts", role = "ctb", email = "kyle.butts@colorado.edu"),
-    person("Henrik", "Sloot", role = "ctb")
+    person("Henrik", "Sloot", role = "ctb"),
+    person("Philippe", "Grosjean", role = c("ctb"), comment = c(ORCID = "0000-0002-2694-9471"))
   )
 Description: Several fast random number generators are provided as C++
   header only libraries: The PCG family by O'Neill (2014

NEWS.md

@@ -4,7 +4,7 @@
 
 * The default RNG has changed from Xoroshiro128+ to Xoroshiro128++. The older generators Xoroshiro128+ and Xoshiro256+ are still available but should only be used for backward compatibility or for generating floating point numbers, i.e. not sampling etc.  ([#57](https://github.com/daqana/dqrng/pull/57) fixing [#56](https://github.com/daqana/dqrng/issues/56))
 * The `dqrng::rng64_t` type has been changed to use `Rcpp::XPtr` instead of `std::shared_ptr` and the functions from `dqrng_sample.h` now expect a reference to `dqrng::random_64bit_generator` instead of `dqrng::rng64_t` ([#70](https://github.com/daqana/dqrng/pull/70) fixing [#63](https://github.com/daqana/dqrng/issues/63))
-* The two argument constructor and `seed`  function from PCG has [surprising properties](https://github.com/imneme/pcg-cpp/issues/91): it is not identical to the one argument version followed by `set_stream(stream)`. For consistency with the new `clone(stream)` method, the two argument versions are no longer used. This influences code the uses multiple stream with PCG together with the tooling from this package, e.g. the example code in the vignette on parallel RNG usage. 
+* The two argument constructor and `seed`  function from PCG has [surprising properties](https://github.com/imneme/pcg-cpp/issues/91): it is not identical to the one argument version followed by `set_stream(stream)`. For consistency with the new `clone(stream)` method, the two argument versions are no longer used. This influences code that uses multiple streams with PCG together with the tooling from this package, e.g. the example code in the vignette on parallel RNG usage. In addition, setting the stream on PCG64 via `dqset.seed(seed, stream)` or at the C++ level using the interface provided by dqrng will be relative to the current stream, i.e. setting `stream=0` will not change the RNG. This is for consistency with the other provided RNGs. You still get the standard behaviour if you are using the C++ classes for PCG directly.
 
 ## Other changes
 
@@ -18,6 +18,7 @@
     * A `clone(stream)` method to allow using the global RNG state for parallel computation
     * New methods `variate<dist>(param)`, `generate<dist>(container, param)` etc. using and inspired by [`randutils`](https://www.pcg-random.org/posts/ease-of-use-without-loss-of-power.html).
 * The scalar functions `dqrng::runif`, `dqrng::rnorm` and `dqrng::rexp` available from `dqrng.h` have been deprecated and will be removed in a future release. Please use the more flexible and faster `dqrng::random_64bit_accessor` together with `variate<Dist>()` instead. The same applies to `dqrng::uniform01` from `dqrng_distribution.h`, which can be replaced by the member function `dqrng::random_64bit_generator::uniform01`. 
+* New template function `dqrng::extra::parallel_generate` in `dqrng_extra/parallel_generate.h` as an example for using the global RNG in a parallel context (fixing [#77](https://github.com/daqana/dqrng/issues/77) in [#82](https://github.com/daqana/dqrng/issues/82) together with Philippe Grosjean)
 
 
 # dgrng 0.3.2

README.Rmd

@@ -38,10 +38,9 @@ Intermediate releases can also be obtained via
 [r-universe](https://rstub.r-universe.dev/dqrng):
 
 ```r
-options(repos = c(
+install.packages('dqrng', repos = c(
   rstub = 'https://rstub.r-universe.dev',
   CRAN = 'https://cloud.r-project.org'))
-install.packages('dqrng')
 ```
 
 ## Example

README.md

@@ -32,10 +32,9 @@ Intermediate releases can also be obtained via
 [r-universe](https://rstub.r-universe.dev/dqrng):
 
 ``` r
-options(repos = c(
+install.packages('dqrng', repos = c(
   rstub = 'https://rstub.r-universe.dev',
   CRAN = 'https://cloud.r-project.org'))
-install.packages('dqrng')
 ```
 
 ## Example
@@ -61,8 +60,8 @@ bm[, 1:4]
 #> # A tibble: 2 × 4
 #>   expression      min   median `itr/sec`
 #>   <bch:expr> <bch:tm> <bch:tm>     <dbl>
-#> 1 rnorm(N)    608.9µs  692.8µs     1342.
-#> 2 dqrnorm(N)   82.7µs   86.2µs    10227.
+#> 1 rnorm(N)    606.2µs  654.7µs     1456.
+#> 2 dqrnorm(N)   82.8µs   85.9µs    10984.
 ```
 
 This is also true for the provided sampling functions with replacement:
@@ -79,10 +78,10 @@ bm[, 1:4]
 #> # A tibble: 4 × 4
 #>   expression                                     min   median `itr/sec`
 #>   <bch:expr>                                <bch:tm> <bch:tm>     <dbl>
-#> 1 sample.int(m, n, replace = TRUE)            7.03ms   7.46ms      133.
-#> 2 sample.int(1000 * m, n, replace = TRUE)     8.87ms    9.4ms      104.
-#> 3 dqsample.int(m, n, replace = TRUE)        292.69µs 316.75µs     2821.
-#> 4 dqsample.int(1000 * m, n, replace = TRUE) 408.94µs 446.19µs     1987.
+#> 1 sample.int(m, n, replace = TRUE)            6.88ms   7.07ms      140.
+#> 2 sample.int(1000 * m, n, replace = TRUE)     8.57ms   8.81ms      112.
+#> 3 dqsample.int(m, n, replace = TRUE)        289.69µs 296.86µs     2834.
+#> 4 dqsample.int(1000 * m, n, replace = TRUE) 407.45µs 489.33µs     1645.
 ```
 
 And without replacement:
@@ -100,11 +99,11 @@ bm[, 1:4]
 #> # A tibble: 5 × 4
 #>   expression                            min   median `itr/sec`
 #>   <bch:expr>                       <bch:tm> <bch:tm>     <dbl>
-#> 1 sample.int(m, n)                  25.43ms   27.8ms      31.8
-#> 2 sample.int(1000 * m, n)           12.99ms   17.4ms      57.0
-#> 3 sample.int(m, n, useHash = TRUE)  10.44ms   12.8ms      77.2
-#> 4 dqsample.int(m, n)               617.83µs  762.9µs     882. 
-#> 5 dqsample.int(1000 * m, n)          1.51ms    2.5ms     334.
+#> 1 sample.int(m, n)                  40.98ms   42.8ms      23.1
+#> 2 sample.int(1000 * m, n)           12.01ms   13.3ms      66.9
+#> 3 sample.int(m, n, useHash = TRUE)   9.35ms   10.4ms      92.4
+#> 4 dqsample.int(m, n)               616.34µs  679.1µs    1265. 
+#> 5 dqsample.int(1000 * m, n)          1.42ms    1.7ms     501.
 ```
 
 Note that sampling from `10^10` elements triggers “long-vector support”

inst/include/dqrng_extra/parallel_generate.h

@@ -0,0 +1,81 @@
+// Copyright 2024 Ralf Stubner
+// Copyright 2024 Philippe Grosjean
+//
+// This file is part of dqrng.
+//
+// dqrng is free software: you can redistribute it and/or modify it
+// under the terms of the GNU Affero General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// dqrng is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Affero General Public License for more details.
+//
+// You should have received a copy of the GNU Affero General Public License
+// along with dqrng.  If not, see <http://www.gnu.org/licenses/>.
+
+#include <sstream>
+#include <dqrng.h>
+#include <RcppParallel/RVector.h>
+#include <omp.h>
+
+namespace dqrng {
+namespace extra {
+template<typename Dist, typename... Params>
+Rcpp::NumericVector parallel_generate(std::size_t n,
+                                      std::size_t threads,
+                                      std::size_t streams,
+                                      Params&&... params) {
+  if (n < streams)
+    streams = n;
+  std::size_t stream_size = n / streams;
+  std::size_t remainder = n % streams;
+
+  // use RcppParallel::RVector as thread safe accessor
+  Rcpp::NumericVector res(Rcpp::no_init(n));
+  RcppParallel::RVector<double> work(res);
+
+  // use global RNG from dqrng
+  dqrng::random_64bit_accessor rng{};
+  std::stringstream buffer;
+
+#ifdef _OPENMP
+  std::size_t maxthreads = omp_get_num_procs();
+  if (threads > maxthreads)
+    threads = maxthreads;
+  // No need for more threads than there are streams
+  if (threads > streams)
+    threads = streams;
+#endif
+
+#pragma omp parallel num_threads(threads)
+{
+  std::size_t start,end;
+
+#pragma omp for schedule(static,1)
+  for (std::size_t i = 0; i < streams; ++i) {
+    if (i < remainder) {
+      start = i * stream_size + i;
+      end = start + stream_size + 1;
+    } else {
+      start = i * stream_size + remainder;
+      end = start + stream_size;
+    }
+    // private RNG in each stream; RNG with i == 0 is identical to global RNG
+    auto prng = rng.clone(i);
+    prng->generate<Dist>(std::begin(work) + start, std::begin(work) + end,
+                         std::forward<Params>(params)...);
+    if (i == 0) {// Save the state of the global RNG's clone
+      buffer << *prng;
+    }
+  }
+}
+  // Make sure that the global RNG advances as well by applying the state
+  // of the global RNG's clone to the global RNG
+  buffer >> rng;
+  return res;
+}
+} // namespace extra
+} // namespace dqrng

inst/include/dqrng_generator.h

@@ -24,7 +24,10 @@
 #include <dqrng_types.h>
 #include <xoshiro.h>
 #include <pcg_random.hpp>
+#include <boost/multiprecision/cpp_int.hpp>
 #include <Rcpp.h>
+#include <convert_seed.h>
+#include <R_randgen.h>
 
 #if defined(__cpp_lib_make_unique) && (__cpp_lib_make_unique >= 201304)
 using std::make_unique;
@@ -98,19 +101,35 @@ inline void random_64bit_wrapper<::dqrng::xoshiro256starstar>::set_stream(result
 
 template<>
 inline void random_64bit_wrapper<pcg64>::set_stream(result_type stream) {
-  gen.set_stream(stream);
+  // set the stream relative to the current stream, i.e. stream = 0 does not change the RNG
+  boost::multiprecision::uint128_t number;
+  std::vector<boost::multiprecision::uint128_t> state;
+  std::stringstream iss;
+  iss << gen;
+  while (iss >> number)
+    state.push_back(number);
+  // state[1] is the current stream
+  // PCG will do 2*stream + 1 to make sure stream is odd; need to revert that here
+  gen.set_stream(pcg_extras::pcg128_t(state[1]/2) + stream);
 }
 
+uint64_t get_seed_from_r() {
+  Rcpp::RNGScope rngScope;
+  Rcpp::IntegerVector seed(2, dqrng::R_random_int);
+  return dqrng::convert_seed<uint64_t>(seed);
+}
+
+
 template<typename RNG = default_64bit_generator>
 typename std::enable_if<!std::is_base_of<random_64bit_generator, RNG>::value, rng64_t>::type
 generator () {
-  return rng64_t(new random_64bit_wrapper<RNG>());
+  return rng64_t(new random_64bit_wrapper<RNG>(get_seed_from_r()));
 }
 
 template<typename RNG = default_64bit_generator>
 typename std::enable_if<std::is_base_of<random_64bit_generator, RNG>::value, rng64_t>::type
 generator () {
-  return rng64_t(new RNG());
+  return rng64_t(new RNG(get_seed_from_r()));
 }
 
 template<typename RNG = default_64bit_generator>

man/dqrng-package.Rd

@@ -24,18 +24,10 @@
 #include <xoshiro.h>
 #include <pcg_random.hpp>
 #include <dqrng_threefry.h>
-#include <convert_seed.h>
-#include <R_randgen.h>
 
 namespace {
-dqrng::rng64_t rng = dqrng::generator();
+dqrng::rng64_t rng = dqrng::generator(56478348);
 std::string rng_kind = "default";
-
-void init() {
-  Rcpp::RNGScope rngScope;
-  Rcpp::IntegerVector seed(2, dqrng::R_random_int);
-  rng->seed(dqrng::convert_seed<uint64_t>(seed));
-}
 }
 
 // [[Rcpp::interfaces(r, cpp)]]
@@ -44,7 +36,7 @@ void init() {
 void dqset_seed(Rcpp::Nullable<Rcpp::IntegerVector> seed,
                 Rcpp::Nullable<Rcpp::IntegerVector> stream = R_NilValue) {
   if (seed.isNull()) {
-    init();
+    rng = dqrng::generator();
   } else {
     uint64_t _seed = dqrng::convert_seed<uint64_t>(seed.as());
     if (stream.isNotNull()) {

src/dqrng.cpp

@@ -72,8 +72,7 @@ test_that("cloned external Xoshiro256++ gives different result only when a diffe
 test_that("cloned external PCG64 gives different result", {
   dqrng::dqRNGkind("PCG64")
   dqset.seed(use_seed)
-  # PCG64 gives a different result as long as you are not selecting the default stream
-  expect_false(cloned_calls(stream = 0))
+  expect_true(cloned_calls(stream = 0))
   expect_false(cloned_calls(stream = 1))
   dqrng::dqRNGkind("default")
 })

tests/testthat/test-external-generator.R

@@ -203,8 +203,6 @@ Here a minimal [SplitMix](https://xoroshiro.di.unimi.it/splitmix64.c) generator
 #include <Rcpp.h>
 // [[Rcpp::depends(dqrng, BH)]]
 #include <dqrng_distribution.h>
-// [[Rcpp::plugins(cpp11)]]
-#include <cstdint>
 
 class SplitMix {
 public:
@@ -254,7 +252,6 @@ For example, here the 64 bit Threefry engine with 13 rounds from package sitmo i
 // [[Rcpp::depends(dqrng, BH, sitmo)]]
 #include <dqrng_distribution.h>
 #include <threefry.h>
-// [[Rcpp::plugins(cpp11)]]
 
 // [[Rcpp::export]]
 Rcpp::NumericVector threefry_rnorm(const int n, const double mean = 0.0, const double sd = 1.0) {
@@ -281,11 +278,10 @@ For example, this function generates random numbers according to the normal dist
 // [[Rcpp::depends(dqrng)]]
 #include <dqrng_generator.h>
 #include <xoshiro.h>
-// [[Rcpp::plugins(cpp11)]]
 
 // [[Rcpp::export]]
 Rcpp::NumericVector std_rnorm(const int n, const double mean = 0.0, const double sd = 1.0) {
-  auto rng = dqrng::generator<dqrng::xoroshiro128plus>(42);
+  auto rng = dqrng::generator<dqrng::xoroshiro128plusplus>(42);
   Rcpp::NumericVector result(n);
   rng->generate<std::normal_distribution>(result, mean, sd);
   return result;
@@ -334,25 +330,17 @@ The parameters of the distributions are adjusted for every draw using `<distribu
 // [[Rcpp::depends(dqrng, BH)]]
 #include <boost/random/binomial_distribution.hpp>
 #include <dqrng_distribution.h>
-#include <xoshiro.h>
-// [[Rcpp::plugins(cpp11)]]   
 
-// aliases for the used distributions
-using binomial = boost::random::binomial_distribution<int>;
-using normal = dqrng::normal_distribution;
 // [[Rcpp::export]]
 Rcpp::NumericMatrix multiple_distributions(int n) {
   auto rng = dqrng::generator<dqrng::xoshiro256plusplus>(42);
-  // distributions with default parameters
-  binomial bernoulli;
-  normal normal;
   Rcpp::NumericMatrix out(n, 3);
   double p = 0.0;
   for (int i = 0; i < n; ++i) {
     p = double(i) / double(n);
-    out(i,0) = bernoulli(*rng, binomial::param_type(1, p));
-    out(i,1) = normal(*rng, normal::param_type(p, 1.0));
-    out(i,2) = normal(*rng, normal::param_type(4.0, 3.0 - p));
+    out(i,0) = rng->variate<boost::random::binomial_distribution<int>>(1, p);
+    out(i,1) = rng->variate<dqrng::normal_distribution>(p, 1.0);
+    out(i,2) = rng->variate<dqrng::normal_distribution>(4.0, 3.0 - p);
   }
   Rcpp::colnames(out) = Rcpp::CharacterVector::create("Bernoulli", "Normal1", "Normal2");
   return out;
@@ -375,24 +363,17 @@ RNG engine of `dqrng`. Please note that the included RNG will be invalidated if
 #include <boost/random/binomial_distribution.hpp>
 #include <dqrng.h>
 #include <dqrng_distribution.h>
-// [[Rcpp::plugins(cpp11)]]   
 
-// aliases for the used distributions
-using binomial = boost::random::binomial_distribution<int>;
-using normal = dqrng::normal_distribution;
 // [[Rcpp::export]]
 Rcpp::NumericMatrix multiple_distributions(int n) {
   auto rng = dqrng::random_64bit_accessor{};
-  // distributions with default parameters
-  binomial bernoulli;
-  normal normal;
   Rcpp::NumericMatrix out(n, 3);
   double p = 0.0;
   for (int i = 0; i < n; ++i) {
     p = double(i) / double(n);
-    out(i,0) = bernoulli(rng, binomial::param_type(1, p));
-    out(i,1) = normal(rng, normal::param_type(p, 1.0));
-    out(i,2) = normal(rng, normal::param_type(4.0, 3.0 - p));
+    out(i,0) = rng.variate<boost::random::binomial_distribution<int>>(1, p);
+    out(i,1) = rng.variate<dqrng::normal_distribution>(p, 1.0);
+    out(i,2) = rng.variate<dqrng::normal_distribution>(4.0, 3.0 - p);
   }
   Rcpp::colnames(out) = Rcpp::CharacterVector::create("Bernoulli", "Normal1", "Normal2");
   return out;