Allow vectors as distribution parameters for the RNGs

[ChristK]

Thank you for all your hard work in this excellent package.

Please consider an enhancement to the RNG functions to allow vectors for all arguments in the RNGs. This would make dq RNGs a direct replacement to R RNGs.

For example runif(c(1, 2), c(-1, 0), c(0, 1)) works but dqrunif(c(1, 2), c(-1, 0), c(0, 1)) is not

[rstub]

That sounds reasonable. Implementation note: The meaning of using a vector in the different places is quite different in base R:

• If the n argument is a scalar, its value is treated as the desired output length. If it is a vector, its length is treated as the desired output length.
• The distribution parameter vectors are then recycled to fit the desired output length.

[ChristK]

The function rscalar_direct2 from #40 (comment) is useful for the suggested enhancement here. For example below with a minor modification to accept a vector for the rate of the exponential distribution:

// [[Rcpp::export]]
Rcpp::NumericVector rscalar_direct2(int n, Rcpp::NumericVector rate) {
  // TODO handle case where n is a vector
  if (rate.length() < n) rate = rep_len(rate, n); // recycling
  Rcpp::NumericVector result(Rcpp::no_init(n));
  for (int i = 0; i < n; ++i) {
    using parm_t = decltype(exponential)::param_type;
    exponential.param(parm_t(rate[i]));
    result[i] = rexp_impl();
  }
  return result;
}

[rstub]

@ChristK Interesting. Renaming the function to dqrexp2 gives the following comparison

> bench::mark(dqrexp(1e5, 1),
+             dqrexp2(1e5, 1),
+             check = FALSE)
# A tibble: 2 x 13
  expression             min median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time
  <bch:expr>        <bch:tm> <bch:>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm>
1 dqrexp(1e+05, 1)  790.75µs 1.05ms      797.   781.3KB     2.07   384     1      482ms
2 dqrexp2(1e+05, 1)   1.66ms 1.82ms      497.    1.53MB     6.27   238     3      479ms
# … with 4 more variables: result <list>, memory <list>, time <list>, gc <list>

So a special case for rate.length() == 1 would be needed from my point of view.

[ChristK]

This might work better?

// [[Rcpp::export]]
Rcpp::NumericVector dqrexp3(Rcpp::IntegerVector n, Rcpp::NumericVector rate) {
 int n_ = n[0];
 if (n.length() > 1) n_ = n.length();
 Rcpp::NumericVector result(Rcpp::no_init(n_));

 if (rate.length() == 1)
   {
   using parm_t = decltype(exponential)::param_type;
   exponential.param(parm_t(rate[0]));
   for (int i = 0; i < n_; ++i)
     result[i] = rexp_impl();
 }
 else
   {
   if (rate.length() < n_) rate = rep_len(rate, n_); // recycling
   for (int i = 0; i < n_; ++i) {
     using parm_t = decltype(exponential)::param_type;
     exponential.param(parm_t(rate[i]));
     result[i] = rexp_impl();
   }
 };
 return result;
}

for large vectors (n> 1e6) in my machine that gives better timings than currently implemented dqrexp

[rstub]

The performance difference is interesting. And since the function pointer rexp_impl was introduced to use a simple Rcpp::NumericVector(n, rexp_impl), one could get rid of it and use exponential(*rng) to gain even more speed. However, limiting the first argument to be an Rcpp::IntegerVector is to restrictive. We will need Rcpp::NumericVector for long-vector support. And in principle other data types should be possible as well, e.g. runif(complex(real = 1:5, imaginary = 1:5)) works. @ChristK, I am not asking you to implement this. I am just recording the situation.

[rstub]

I am considering to close this with appropriate documentation since the current development version allows the RNG to be registered as user defined RNG. This way, one can use the functions from stats that do support vectors as distribution parameters but still get the faster RNGs from dqrng, c.f.:

library(dqrng)
bench::mark(stats::rnorm(1e6), dqrng::dqrnorm(1e6), check = FALSE)[ , 1:4]
#> # A tibble: 2 × 4
#>   expression                 min   median `itr/sec`
#>   <bch:expr>            <bch:tm> <bch:tm>     <dbl>
#> 1 stats::rnorm(1e+06)    21.11ms  21.49ms      46.7
#> 2 dqrng::dqrnorm(1e+06)   4.96ms   5.29ms     188.
register_methods()
bench::mark(stats::rnorm(1e6), dqrng::dqrnorm(1e6), check = FALSE)[ , 1:4]
#> # A tibble: 2 × 4
#>   expression                 min   median `itr/sec`
#>   <bch:expr>            <bch:tm> <bch:tm>     <dbl>
#> 1 stats::rnorm(1e+06)     7.64ms   8.21ms      120.
#> 2 dqrng::dqrnorm(1e+06)   3.52ms    3.7ms      268.
set.seed(42)
dqrunif(6)
#> [1] 0.90837042 0.33061240 0.95093543 0.20447122 0.34282983 0.09195553
set.seed(42)
runif(6, c(0,1), c(1,2))
#> [1] 0.9083704 1.3306124 0.9509354 1.2044712 0.3428298 1.0919555
restore_methods()

^{Created on 2024-04-23 with reprex v2.1.0}

@ChristK: What do you think about that possibility?

[ChristK]

Yes, of course, feel free to close it. It was a mere suggestion and I totally understand if that's not a priority. Many thanks for all your hard work with this package. I use it daily for my simulations.

[ChristK]

Closing as out of scope (see discussion above)