diff --git a/404.html b/404.html index ae63d0f..481d3b6 100644 --- a/404.html +++ b/404.html @@ -6,7 +6,7 @@
vignettes/cpp-api.Rmd
cpp-api.Rmd
Rcpp::IntegerVector dqrng::dqsample_int(int n, int size, bool replace = false,
- Rcpp::Nullable<Rcpp::NumericVector> probs = R_NilValue,
- int offset = 0)
-Rcpp::NumericVector dqrng::dqsample_num(double n, double size, bool replace = false,
- Rcpp::Nullable<Rcpp::NumericVector> probs = R_NilValue,
- int offset = 0)
Rcpp::IntegerVector dqrng::dqsample_int(int n, int size, bool replace = false,
+ Rcpp::Nullable<Rcpp::NumericVector> probs = R_NilValue,
+ int offset = 0)
+Rcpp::NumericVector dqrng::dqsample_num(double n, double size, bool replace = false,
+ Rcpp::Nullable<Rcpp::NumericVector> probs = R_NilValue,
+ int offset = 0)
n
The two functions are used for “normal” and “long-vector” support in R.
+std::vector<std::string> dqrng::dqrng_get_state()
+void dqrng::dqrng_set_state(std::vector<std::string> state)
state
std::vector<std::string>
as produced by
+dqrng_get_state()
+Direct usage of this method is discouraged. The
+prefered way of accesing the global RNG is to instantiate
+dqrng::random_64bit_accessor
within your function. Note
+that you MUST NOT delete the global RNG. Using
+dqrng::random_64bit_accessor
makes this impossible. In
+addition, you SHOULD NOT store a reference to the RNG permanently,
+because it can be invalidated by calls to dqRNGkind
.
+Therefore, instances of dqrng::random_64bit_accessor
SHOULD
+be stored as (non-static) variables in functions.
Note that dqrng::random_64bit_accessor
supports UniformRandomBitGenerator
+and can therefore be used together with any C++11 distribtion function.
+In addition, the follwing functions are supported, since they are
+inherited from the abstract parent class
+random_64bit_generator
:
// clone RNG and seelct a different stream
+std::unique_ptr<random_64bit_generator> clone(uint64_t stream)
+// uniform doubles in [0,1) and double-int-pairs
+double uniform01()
+std::pair<double, int> generate_double_8bit_pair()
+// uniform integers in a range
+uint32_t operator() (uint32_t range)
+uint64_t operator() (uint64_t range)
+ResultType variate<ResultType, DistTmpl>(param1, ... paramN)
+generate<DistTmpl>(container, param1, ... paramN)
+generate<DistTmpl>(start, end, param1, ... paramN)
+Dist::result_type variate<Dist>(param1, ... paramN)
+generate<Dist>(container, param1, ... paramN)
+generate<Dist>(start, end, param1, ... paramN)
stream
range
[0, range]
+ResultType
DistTmpl
+DistTmpl
std::uniform_distribution
.
+DistTmpl<ResultType>
defines the full distribution.
+Dist
std::uniform_distribution<double>
or
+dqrng::normal_distriubtion
.
+param1, ... paramN
container
std::begin
and
+std::end
.
+start, end
Site built with pkgdown 2.0.7.
+Site built with pkgdown 2.0.8.
diff --git a/articles/dqrng.html b/articles/dqrng.html index 69a374c..1c4976c 100644 --- a/articles/dqrng.html +++ b/articles/dqrng.html @@ -6,7 +6,7 @@vignettes/dqrng.Rmd
dqrng.Rmd
Using dqrng is about three times faster:
library(dqrng)
@@ -174,46 +174,46 @@ Usage from Rsystem.time(cat("pi ~= ", piR(N, rng = dqrng::dqrunif), "\n"))
#> pi ~= 3.141457
#> user system elapsed
-#> 0.222 0.041 0.262
Since the calculations add a constant off-set, the speed-up for the RNGs alone has to be even greater:
system.time(stats::runif(N))
#> user system elapsed
-#> 0.078 0.013 0.090
+#> 0.087 0.003 0.090
system.time(dqrng::dqrunif(N))
#> user system elapsed
-#> 0.049 0.008 0.056
Similar for the exponential distribution:
system.time(stats::rexp(N))
#> user system elapsed
-#> 0.338 0.004 0.342
+#> 0.342 0.005 0.346
system.time(dqrng::dqrexp(N))
#> user system elapsed
-#> 0.063 0.008 0.071
And for the normal distribution:
system.time(stats::rnorm(N))
#> user system elapsed
-#> 0.367 0.000 0.367
+#> 0.366 0.000 0.366
system.time(dqrng::dqrnorm(N))
#> user system elapsed
-#> 0.086 0.008 0.093
As well as for sampling with and without replacement:
system.time(for (i in 1:100) sample.int(N, N/100, replace = TRUE))
#> user system elapsed
-#> 0.573 0.004 0.577
+#> 0.576 0.000 0.576
system.time(for (i in 1:100) dqsample.int(N, N/100, replace = TRUE))
#> user system elapsed
-#> 0.031 0.016 0.047
+#> 0.031 0.017 0.047
system.time(for (i in 1:100) sample.int(N, N/100))
#> user system elapsed
-#> 1.533 0.304 1.837
+#> 1.430 0.293 1.722
system.time(for (i in 1:100) dqsample.int(N, N/100))
#> user system elapsed
-#> 0.061 0.008 0.069
It is also possible to register the supplied generators as
user-supplied RNGs. This way set.seed()
and
dqset.seed()
influence both (dq)runif
and
@@ -373,15 +373,14 @@
Since SplitMix is a very fast RNG, the speed of this function is
comparable to dqrnorm
. Generally speaking you can use any
C++11 compliant RNG with 64 bit output size. For example, here the 64
@@ -395,16 +394,15 @@
Note that for the (recommended) Threefry engine with 20 rounds some
additional integration is provided in the dqrng_threefry.h
header file.
Typically this is not as fast as dqrnorm
, but the
technique is useful to support distributions not (yet) included in
dqrng. Note however, that the algorithms used for the distributions from
@@ -560,7 +557,7 @@
Site built with pkgdown 2.0.7.
+Site built with pkgdown 2.0.8.
diff --git a/articles/index.html b/articles/index.html index b8d0e76..26dfae2 100644 --- a/articles/index.html +++ b/articles/index.html @@ -1,5 +1,5 @@ -vignettes/parallel.Rmd
parallel.Rmd
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:
#include <Rcpp.h>
+// [[Rcpp::depends(dqrng, BH)]]
+#include <dqrng.h>
+#include <dqrng_distribution.h>
+// [[Rcpp::plugins(cpp11)]]
+// [[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 that the thread local RNG uses std::unique_ptr
+instead of Rcpp::Xptr
since it is used in parallel
+code.
dqrng::generator
. The resulting correlation
matrix should be close to the identity matrix if the different streams
are independent:
-#include <Rcpp.h>
-// [[Rcpp::depends(dqrng, BH)]]
-#include <pcg_random.hpp>
-#include <dqrng_sample.h>
-// [[Rcpp::plugins(cpp11)]]
-// [[Rcpp::depends(RcppParallel)]]
-#include <RcppParallel.h>
-
-struct RandomFill : public RcppParallel::Worker {
- RcppParallel::RMatrix<int> output;
- uint64_t seed;
-
- RandomFill(Rcpp::IntegerMatrix output, const uint64_t seed) : output(output), seed(seed) {};
-
- void operator()(std::size_t begin, std::size_t end) {
- auto rng = dqrng::generator<pcg64>(seed, end);
- for (std::size_t col = begin; col < end; ++col) {
- auto sampled = dqrng::sample::sample<std::vector<int>, uint32_t>(*rng, 100000, output.nrow(), true);
- RcppParallel::RMatrix<int>::Column column = output.column(col);
- std::copy(sampled.begin(), sampled.end(), column.begin());
- }
- }
-};
-
-// [[Rcpp::export]]
-Rcpp::IntegerMatrix parallel_random_matrix(const int n, const int m, const int ncores) {
- Rcpp::IntegerMatrix 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)
-*/
#include <Rcpp.h>
+// [[Rcpp::depends(dqrng, BH)]]
+#include <dqrng_generator.h>
+#include <dqrng_sample.h>
+// [[Rcpp::plugins(cpp11)]]
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallel.h>
+
+struct RandomFill : public RcppParallel::Worker {
+ RcppParallel::RMatrix<int> output;
+ uint64_t seed;
+
+ RandomFill(Rcpp::IntegerMatrix output, const uint64_t seed) : output(output), seed(seed) {};
+
+ void operator()(std::size_t begin, std::size_t end) {
+ auto rng = dqrng::generator<pcg64>(seed, end);
+ for (std::size_t col = begin; col < end; ++col) {
+ auto sampled = dqrng::sample::sample<std::vector<int>, uint32_t>(*rng, 100000, output.nrow(), true);
+ RcppParallel::RMatrix<int>::Column column = output.column(col);
+ std::copy(sampled.begin(), sampled.end(), column.begin());
+ }
+ }
+};
+
+// [[Rcpp::export]]
+Rcpp::IntegerMatrix parallel_random_matrix(const int n, const int m, const int ncores) {
+ Rcpp::IntegerMatrix 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,] 67984 16279 69262 7126 21441 37720 51107 51045
-[2,] 69310 21713 82885 81157 54051 5261 91165 17833
-[3,] 76742 31232 78953 4626 94939 29416 85652 78296
-[4,] 76349 47427 1770 37957 33888 59134 94591 65793
-[5,] 85008 89224 43493 7925 60866 2464 14080 10763
-[6,] 38017 88509 51195 73086 1883 68193 75259 62216
+[1,] 16144 89512 16144 53103 16144 41760 16144 99788
+[2,] 21200 86552 21199 38327 21198 90293 21198 45713
+[3,] 124 60045 97553 79574 99977 84423 87126 35715
+[4,] 6590 93044 18536 98945 66969 53338 38426 53560
+[5,] 80618 42051 15815 79380 71590 98006 87220 63036
+[6,] 49954 90168 12614 32488 54145 68029 34733 48989
Correlation matrix:
[1,] 1
[2,] 1
-[3,] ? 1
-[4,] ? 1
-[5,] 1
-[6,] ? ? ? 1
-[7,] ? 1
-[8,] ? 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.
+As in the previous section, we can again switch from seeding our own +RNG to starting from the global RNG and selecting separate substreams +from there. RNG selection and seeding can then be handled from R:
+#include <Rcpp.h>
+// [[Rcpp::depends(dqrng, BH)]]
+#include <dqrng.h>
+#include <dqrng_sample.h>
+// [[Rcpp::plugins(cpp11)]]
+// [[Rcpp::depends(RcppParallel)]]
+#include <RcppParallel.h>
+
+struct RandomFill : public RcppParallel::Worker {
+ RcppParallel::RMatrix<int> output;
+ dqrng::random_64bit_accessor rng{};
+
+ RandomFill(Rcpp::IntegerMatrix output) : output(output) {};
+
+ void operator()(std::size_t begin, std::size_t end) {
+ auto lrng = rng.clone(end);
+ for (std::size_t col = begin; col < end; ++col) {
+ auto sampled = dqrng::sample::sample<std::vector<int>, uint32_t>(*lrng, 100000, output.nrow(), true);
+ RcppParallel::RMatrix<int>::Column column = output.column(col);
+ std::copy(sampled.begin(), sampled.end(), column.begin());
+ }
+ }
+};
+
+// [[Rcpp::export]]
+Rcpp::IntegerMatrix parallel_random_matrix(const int n, const int m, const int ncores) {
+ Rcpp::IntegerMatrix res(n, m);
+ RandomFill randomFill(res);
+ RcppParallel::parallelFor(0, m, randomFill, m/ncores + 1);
+ return res;
+}
+
+/*** R
+dqRNGkind("pcg64")
+dqset.seed(42)
+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)
+*/
vignettes/sample.Rmd
sample.Rmd
Site built with pkgdown 2.0.7.
+Site built with pkgdown 2.0.8.
diff --git a/authors.html b/authors.html index 5b680d0..41e3ef5 100644 --- a/authors.html +++ b/authors.html @@ -1,5 +1,5 @@ -sigma <- matrix(c(4,2,2,3), ncol=2)
x <- dqrmvnorm(n=500, mean=c(1,2), sigma=sigma)
colMeans(x)
-#> [1] 0.9358251 2.0251444
+#> [1] 0.8614361 2.0263967
var(x)
#> [,1] [,2]
-#> [1,] 3.909685 1.884089
-#> [2,] 1.884089 2.976024
+#> [1,] 3.930309 1.917535
+#> [2,] 1.917535 2.881480
plot(x)