Add clon(stream) method for the RNGs

inst/include/dqrng_generator.h

@@ -53,18 +53,25 @@ class random_64bit_wrapper : public random_64bit_generator {
     has_cache = true;
     return random >> 32;
   }
+  void set_stream(result_type stream) {throw std::runtime_error("Stream handling not supported for this RNG!");}
 
 protected:
   virtual void output(std::ostream& ost) const override {ost << gen;}
   virtual void input(std::istream& ist) override {ist >> gen;}
 
 public:
   random_64bit_wrapper() : gen() {};
+  random_64bit_wrapper(RNG _gen) : gen(_gen) {};
   random_64bit_wrapper(result_type seed) : gen(seed) {};
   random_64bit_wrapper(result_type seed, result_type stream) : gen(seed, stream) {};
   virtual result_type operator() () override {return this->bit64();}
   virtual void seed(result_type seed) override {cache = false; gen.seed(seed);}
   virtual void seed(result_type seed, result_type stream) override {throw std::runtime_error("Stream handling not supported for this RNG!");}
+  virtual std::unique_ptr<random_64bit_generator> clone(result_type stream) override {
+    auto rng = std::make_unique<random_64bit_wrapper<RNG>>(gen);
+    rng->set_stream(stream);
+    return rng;
+  }
 
   /*
    * https://raw.githubusercontent.com/imneme/bounded-rands/3d71f53c975b1e5b29f2f3b05a74e26dab9c3d84/bounded32.cpp
@@ -143,50 +150,84 @@ template<>
 inline void random_64bit_wrapper<::dqrng::xoroshiro128plus>::seed(result_type seed, result_type stream) {
     gen.seed(seed);
     gen.jump(stream);
-    cache = false;
+    has_cache = false;
+}
+
+template<>
+inline void random_64bit_wrapper<::dqrng::xoroshiro128plus>::set_stream(result_type stream) {
+  gen.jump(stream);
 }
 
 template<>
 inline void random_64bit_wrapper<::dqrng::xoroshiro128plusplus>::seed(result_type seed, result_type stream) {
   gen.seed(seed);
   gen.jump(stream);
-  cache = false;
+  has_cache = false;
+}
+
+template<>
+inline void random_64bit_wrapper<::dqrng::xoroshiro128plusplus>::set_stream(result_type stream) {
+  gen.jump(stream);
 }
 
 template<>
 inline void random_64bit_wrapper<::dqrng::xoroshiro128starstar>::seed(result_type seed, result_type stream) {
   gen.seed(seed);
   gen.jump(stream);
-  cache = false;
+  has_cache = false;
+}
+
+template<>
+inline void random_64bit_wrapper<::dqrng::xoroshiro128starstar>::set_stream(result_type stream) {
+  gen.jump(stream);
 }
 
 template<>
 inline void random_64bit_wrapper<::dqrng::xoshiro256plus>::seed(result_type seed, result_type stream) {
   gen.seed(seed);
   gen.long_jump(stream);
-  cache = false;
+  has_cache = false;
+}
+
+template<>
+inline void random_64bit_wrapper<::dqrng::xoshiro256plus>::set_stream(result_type stream) {
+  gen.long_jump(stream);
 }
 
 template<>
 inline void random_64bit_wrapper<::dqrng::xoshiro256plusplus>::seed(result_type seed, result_type stream) {
   gen.seed(seed);
   gen.long_jump(stream);
-  cache = false;
+  has_cache = false;
+}
+
+template<>
+inline void random_64bit_wrapper<::dqrng::xoshiro256plusplus>::set_stream(result_type stream) {
+  gen.long_jump(stream);
 }
 
 template<>
 inline void random_64bit_wrapper<::dqrng::xoshiro256starstar>::seed(result_type seed, result_type stream) {
     gen.seed(seed);
     gen.long_jump(stream);
-    cache = false;
+    has_cache = false;
+}
+
+template<>
+inline void random_64bit_wrapper<::dqrng::xoshiro256starstar>::set_stream(result_type stream) {
+  gen.long_jump(stream);
 }
 
 template<>
 inline void random_64bit_wrapper<pcg64>::seed(result_type seed, result_type stream) {
     gen.seed(seed, stream);
-    cache = false;
+    has_cache = false;
 }
 
+template<>
+inline void random_64bit_wrapper<pcg64>::set_stream(result_type stream) {
+  gen.set_stream(stream);
+}
 
 template<typename RNG = default_64bit_generator>
 typename std::enable_if<!std::is_base_of<random_64bit_generator, RNG>::value, rng64_t>::type

inst/include/dqrng_types.h

@@ -22,6 +22,7 @@
 
 #include <mystdint.h>
 #include <stdexcept>
+#include <memory>
 #include <Rcpp/XPtr.h>
 
 namespace dqrng {
@@ -37,6 +38,7 @@ class random_64bit_generator {
   virtual result_type operator() () = 0;
   virtual void seed(result_type seed) = 0;
   virtual void seed(result_type seed, result_type stream) = 0;
+  virtual std::unique_ptr<random_64bit_generator> clone(result_type stream) = 0;
   static constexpr result_type min() {return 0;};
   static constexpr result_type max() {return UINT64_MAX;};
   virtual uint32_t operator() (uint32_t range) = 0;
@@ -68,6 +70,8 @@ class random_64bit_accessor : public random_64bit_generator {
 
 public:
   explicit random_64bit_accessor();
+  // defined in dqrng.h:
+  // random_64bit_accessor() : gen(dqrng::get_rng()) {}
 
   virtual result_type operator() () override {
     return (*gen)();
@@ -81,6 +85,10 @@ class random_64bit_accessor : public random_64bit_generator {
     throw std::runtime_error("Seed handling not supported for this class!");
   };
 
+  virtual std::unique_ptr<random_64bit_generator> clone(result_type stream) override {
+    return gen->clone(stream);
+  };
+
   virtual uint32_t operator() (uint32_t range) override {
     return (*gen)(range);
   }

vignettes/parallel.Rmd

@@ -105,7 +105,7 @@ std::vector<double> random_sum(int n, int m) {
 // [[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::xoshiro256plusplus rng(42);              // properly seeded rng
+  dqrng::xoshiro256plusplus rng(42);          // properly seeded rng
   std::vector<double> res(m);
   // ok to use rng here
   
@@ -145,6 +145,71 @@ Result:
     1 parallel_random_sum(1e+07, 8, 4)   98.3ms     99ms     10.1 
     2 random_sum(1e+07, 8)              270.2ms    271ms      3.68
 
+In the previous example it was necessary to create our own RNG object. It would be nice if it were possible to use dqrng's global RNG instead. To showcase this, we implement the same procedure as before, but this time make use of `dqrng::random_64bit_accessor` to access the global RNG and its `clone(stream)` method that creates a cloned version with a different stream:
+
+```{r, eval=FALSE, engine='Rcpp'}
+#include <Rcpp.h>
+// [[Rcpp::depends(dqrng, BH)]]
+#include <dqrng.h>
+#include <dqrng_distribution.h>
+// [[Rcpp::plugins(cpp14)]]
+// [[Rcpp::plugins(openmp)]]
+#include <omp.h>
+
+// [[Rcpp::export]]
+std::vector<double> random_sum(int n, int m) {
+  dqrng::uniform_distribution dist(0.0, 1.0); // Uniform distribution [0,1)
+  auto rng = dqrng::random_64bit_accessor{};  // properly seeded rng
+  std::vector<double> res(m);
+  for (int i = 0; i < m; ++i) {
+    double lres(0);
+    for (int j = 0; j < n; ++j) {
+      lres += dist(rng);
+    }
+    res[i] = lres / n;
+  }
+  return res;
+}
+
+// [[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)
+  auto rng = dqrng::random_64bit_accessor{};  // properly seeded rng
+  std::vector<double> res(m);
+  // ok to use rng here
+
+#pragma omp parallel num_threads(ncores)
+{
+  // make thread local copy of rng and advance it by 1 ... ncores jumps
+  auto lrng = rng.clone(omp_get_thread_num() + 1);      
+
+#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
+bm <- bench::mark(
+  parallel_random_sum(1e7, 8, 4),
+  random_sum(1e7, 8),
+  check = FALSE
+)
+bm[,1:4]
+*/
+```
+
+Notes:
+
+* The thread local RNG does not make use of `Rcpp::Xptr` since it is used in parallel code. Instead a `std::unique_ptr` is used that is created using `std::make_unique`. Therefore, using this feature requires at least C++14. 
+* At the moment this kind of stream handling is not supported for the Threefry engine.
 
 # PCG: multiple streams with RcppParallel