diff --git a/DESCRIPTION b/DESCRIPTION index ceaac2c3..def22610 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.3.18 +Version: 0.3.19 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), diff --git a/NAMESPACE b/NAMESPACE index 5378f687..0db1ef65 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -71,13 +71,9 @@ export(monty_random_real) export(monty_random_truncated_normal) export(monty_random_uniform) export(monty_random_weibull) -export(monty_rng) export(monty_rng_create) -export(monty_rng_distributed_pointer) -export(monty_rng_distributed_state) export(monty_rng_jump) export(monty_rng_long_jump) -export(monty_rng_pointer) export(monty_rng_set_state) export(monty_rng_state) export(monty_runner_callr) diff --git a/R/cpp11.R b/R/cpp11.R index 818b87bd..018cbd6a 100644 --- a/R/cpp11.R +++ b/R/cpp11.R @@ -1,5 +1,9 @@ # Generated by cpp11: do not edit by hand +monty_rng_alloc <- function(r_seed, n_streams, deterministic) { + .Call(`_monty_monty_rng_alloc`, r_seed, n_streams, deterministic) +} + cpp_monty_rng_state <- function(ptr) { .Call(`_monty_cpp_monty_rng_state`, ptr) } @@ -192,114 +196,6 @@ cpp_monty_random_n_multinomial <- function(n_samples, r_size, r_prob, ptr) { .Call(`_monty_cpp_monty_random_n_multinomial`, n_samples, r_size, r_prob, ptr) } -monty_rng_alloc <- function(r_seed, n_streams, deterministic) { - .Call(`_monty_monty_rng_alloc`, r_seed, n_streams, deterministic) -} - -monty_legacy_rng_jump <- function(ptr) { - invisible(.Call(`_monty_monty_legacy_rng_jump`, ptr)) -} - -monty_legacy_rng_long_jump <- function(ptr) { - invisible(.Call(`_monty_monty_legacy_rng_long_jump`, ptr)) -} - -monty_rng_random_real <- function(ptr, n, n_threads) { - .Call(`_monty_monty_rng_random_real`, ptr, n, n_threads) -} - -monty_rng_random_normal <- function(ptr, n, n_threads, algorithm) { - .Call(`_monty_monty_rng_random_normal`, ptr, n, n_threads, algorithm) -} - -monty_rng_uniform <- function(ptr, n, r_min, r_max, n_threads) { - .Call(`_monty_monty_rng_uniform`, ptr, n, r_min, r_max, n_threads) -} - -monty_rng_exponential_rate <- function(ptr, n, r_rate, n_threads) { - .Call(`_monty_monty_rng_exponential_rate`, ptr, n, r_rate, n_threads) -} - -monty_rng_exponential_mean <- function(ptr, n, r_mean, n_threads) { - .Call(`_monty_monty_rng_exponential_mean`, ptr, n, r_mean, n_threads) -} - -monty_rng_normal <- function(ptr, n, r_mean, r_sd, n_threads, algorithm) { - .Call(`_monty_monty_rng_normal`, ptr, n, r_mean, r_sd, n_threads, algorithm) -} - -monty_rng_binomial <- function(ptr, n, r_size, r_prob, n_threads) { - .Call(`_monty_monty_rng_binomial`, ptr, n, r_size, r_prob, n_threads) -} - -monty_rng_beta_binomial_ab <- function(ptr, n, r_size, r_a, r_b, n_threads) { - .Call(`_monty_monty_rng_beta_binomial_ab`, ptr, n, r_size, r_a, r_b, n_threads) -} - -monty_rng_beta_binomial_prob <- function(ptr, n, r_size, r_prob, r_rho, n_threads) { - .Call(`_monty_monty_rng_beta_binomial_prob`, ptr, n, r_size, r_prob, r_rho, n_threads) -} - -monty_rng_negative_binomial_prob <- function(ptr, n, r_size, r_prob, n_threads) { - .Call(`_monty_monty_rng_negative_binomial_prob`, ptr, n, r_size, r_prob, n_threads) -} - -monty_rng_negative_binomial_mu <- function(ptr, n, r_size, r_mu, n_threads) { - .Call(`_monty_monty_rng_negative_binomial_mu`, ptr, n, r_size, r_mu, n_threads) -} - -monty_rng_hypergeometric <- function(ptr, n, r_n1, r_n2, r_k, n_threads) { - .Call(`_monty_monty_rng_hypergeometric`, ptr, n, r_n1, r_n2, r_k, n_threads) -} - -monty_rng_gamma_scale <- function(ptr, n, r_shape, r_scale, n_threads) { - .Call(`_monty_monty_rng_gamma_scale`, ptr, n, r_shape, r_scale, n_threads) -} - -monty_rng_gamma_rate <- function(ptr, n, r_shape, r_rate, n_threads) { - .Call(`_monty_monty_rng_gamma_rate`, ptr, n, r_shape, r_rate, n_threads) -} - -monty_rng_poisson <- function(ptr, n, r_lambda, n_threads) { - .Call(`_monty_monty_rng_poisson`, ptr, n, r_lambda, n_threads) -} - -monty_rng_cauchy <- function(ptr, n, r_location, r_scale, n_threads) { - .Call(`_monty_monty_rng_cauchy`, ptr, n, r_location, r_scale, n_threads) -} - -monty_rng_beta <- function(ptr, n, r_a, r_b, n_threads) { - .Call(`_monty_monty_rng_beta`, ptr, n, r_a, r_b, n_threads) -} - -monty_rng_multinomial <- function(ptr, n, r_size, r_prob, n_threads) { - .Call(`_monty_monty_rng_multinomial`, ptr, n, r_size, r_prob, n_threads) -} - -monty_rng_truncated_normal <- function(ptr, n, r_mean, r_sd, r_min, r_max, n_threads) { - .Call(`_monty_monty_rng_truncated_normal`, ptr, n, r_mean, r_sd, r_min, r_max, n_threads) -} - -monty_legacy_rng_state <- function(ptr) { - .Call(`_monty_monty_legacy_rng_state`, ptr) -} - -monty_legacy_rng_set_state <- function(ptr, r_state) { - invisible(.Call(`_monty_monty_legacy_rng_set_state`, ptr, r_state)) -} - -monty_rng_pointer_init <- function(n_streams, seed, long_jump, algorithm) { - .Call(`_monty_monty_rng_pointer_init`, n_streams, seed, long_jump, algorithm) -} - -monty_rng_pointer_sync <- function(obj, algorithm) { - invisible(.Call(`_monty_monty_rng_pointer_sync`, obj, algorithm)) -} - -test_rng_pointer_get <- function(obj, n_streams) { - .Call(`_monty_test_rng_pointer_get`, obj, n_streams) -} - -test_xoshiro_run <- function(obj) { - .Call(`_monty_test_xoshiro_run`, obj) +test_xoshiro_run <- function(ptr) { + .Call(`_monty_test_xoshiro_run`, ptr) } diff --git a/R/dsl.R b/R/dsl.R index f0e67648..1373cbad 100644 --- a/R/dsl.R +++ b/R/dsl.R @@ -117,7 +117,7 @@ monty_dsl_parse <- function(x, type = NULL, gradient = NULL, fixed = NULL) { ##' provided! ##' ##' * `sample`: A function to sample from the distribution, given (as -##' a first argument) a rng object (see [monty_rng]) +##' a first argument) a `monty_rng` object (see [monty_rng_create]) ##' ##' @export ##' @examples diff --git a/R/model.R b/R/model.R index 80123b7e..e9fdf238 100644 --- a/R/model.R +++ b/R/model.R @@ -118,13 +118,13 @@ monty_model_properties <- function(has_gradient = NULL, ##' `(-Inf, Inf)`. ##' ##' * `direct_sample`: A function to sample directly from the -##' parameter space, given a [monty_rng] object to sample from. -##' In the case where a model returns a posterior (e.g., in Bayesian -##' inference), this is assumed to be sampling from the prior. -##' We'll use this for generating initial conditions for MCMC where -##' those are not given, and possibly other uses. If not given then -##' when using [monty_sample()] the user will have to provide a -##' vector of initial states. +##' parameter space, given a `monty_rng` object to sample from (see +##' [monty_rng_create]). In the case where a model returns a +##' posterior (e.g., in Bayesian inference), this is assumed to be +##' sampling from the prior. We'll use this for generating initial +##' conditions for MCMC where those are not given, and possibly +##' other uses. If not given then when using [monty_sample()] the +##' user will have to provide a vector of initial states. ##' ##' * `gradient`: A function to compute the gradient of `density` with ##' respect to the parameter vector; takes a parameter vector and @@ -141,11 +141,11 @@ monty_model_properties <- function(has_gradient = NULL, ##' look after their own stream, and that they may need many ##' streams). Models that provide this method are assumed to be ##' stochastic; however, you can use the `is_stochastic` property -##' (via [monty_model_properties()]) to override this (e.g., to -##' run a stochastic model with its deterministic expectation). -##' This function takes a raw vector of random number state from -##' [monty_rng] and uses it to set the random number state for -##' your model; this is derived from the random number stream for a +##' (via [monty_model_properties()]) to override this (e.g., to run +##' a stochastic model with its deterministic expectation). This +##' function takes a raw vector of random number state from a +##' `monty_rng` and uses it to set the random number state for your +##' model; this is derived from the random number stream for a ##' particular chain, jumped ahead. ##' ##' * `get_rng_state`: A function to get the RNG state; must be @@ -304,8 +304,8 @@ monty_model_gradient <- function(model, parameters, named = FALSE) { ##' ##' @inheritParams monty_model_gradient ##' -##' @param rng Random number state, created by [monty_rng]. Use of -##' an RNG with more than one stream may or may not work as +##' @param rng Random number state, created by [monty_rng_create]. +##' Use of an RNG with more than one stream may or may not work as ##' expected; this is something we need to tidy up (`mrc-5292`) ##' ##' @return A vector or matrix of sampled parameters diff --git a/R/random.R b/R/random.R index 6d0e1069..a9efbe4b 100644 --- a/R/random.R +++ b/R/random.R @@ -156,6 +156,20 @@ monty_rng_long_jump <- function(state, n = 1) { } +## This was used previously, we have a variant of this in dust2, as +## well. +monty_rng_distributed_state <- function(seed = NULL, n_nodes = 1L) { + rng <- monty_rng_create(seed = seed) + ret <- vector("list", n_nodes) + for (i in seq_len(n_nodes)) { + ret[[i]] <- monty_rng_state(rng) + if (i < n_nodes) { + monty_rng_long_jump(rng) + } + } + ret +} + ##' @export print.monty_rng_state <- function(x, ...) { diff --git a/R/rng.R b/R/rng.R deleted file mode 100644 index 1aac6dfb..00000000 --- a/R/rng.R +++ /dev/null @@ -1,478 +0,0 @@ -##' @title Monty Random Number Generator -##' -##' @description Create an object that can be used to generate random -##' numbers with the same RNG as monty uses internally. This is -##' primarily meant for debugging and testing the underlying C++ -##' rather than a source of random numbers from R. -##' -##' # Warning -##' -##' This interface is subject to change in the near future, we do not -##' recommend its use in user code. -##' -##' # Running multiple streams, perhaps in parallel -##' -##' The underlying random number generators are designed to work in -##' parallel, and with random access to parameters (see -##' `vignette("rng")` for more details). However, this is usually -##' done within the context of running a model where each particle -##' sees its own stream of numbers. We provide some support for -##' running random number generators in parallel, but any speed -##' gains from parallelisation are likely to be somewhat eroded by -##' the overhead of copying around a large number of random numbers. -##' -##' All the random distribution functions support an argument -##' `n_threads` which controls the number of threads used. This -##' argument will *silently* have no effect if your installation -##' does not support OpenMP. -##' -##' Parallelisation will be performed at the level of the stream, -##' where we draw `n` numbers from *each* stream for a total of `n * -##' n_streams` random numbers using `n_threads` threads to do this. -##' Setting `n_threads` to be higher than `n_streams` will therefore -##' have no effect. If running on somebody else's system (e.g., an -##' HPC, CRAN) you must respect the various environment variables -##' that control the maximum allowable number of threads. -##' -##' With the exception of `random_real`, each random number -##' distribution accepts parameters; the interpretations of these -##' will depend on `n`, `n_streams` and their rank. -##' -##' * If a scalar then we will use the same parameter value for every draw -##' from every stream -##' -##' * If a vector with length `n` then we will draw `n` random -##' numbers per stream, and every stream will use the same parameter -##' value for every stream for each draw (but a different, -##' shared, parameter value for subsequent draws). -##' -##' * If a matrix is provided with one row and `n_streams` -##' columns then we use different parameters for each stream, but -##' the same parameter for each draw. -##' -##' * If a matrix is provided with `n` rows and `n_streams` -##' columns then we use a parameter value `[i, j]` for the `i`th -##' draw on the `j`th stream. -##' -##' The rules are slightly different for the `prob` argument to -##' `multinomial` as for that `prob` is a vector of values. As such -##' we shift all dimensions by one: -##' -##' * If a vector we use same `prob` every draw from every stream -##' and there are `length(prob)` possible outcomes. -##' -##' * If a matrix with `n` columns then vary over each draw (the -##' `i`th draw using vector `prob[, i]` but shared across all -##' streams. There are `nrow(prob)` possible outcomes. -##' -##' * If a 3d array is provided with 1 column and `n_streams` -##' "layers" (the third dimension) then we use then we use different -##' parameters for each stream, but the same parameter for each -##' draw. -##' -##' * If a 3d array is provided with `n` columns and `n_streams` -##' "layers" then we vary over both draws and streams so that with -##' use vector `prob[, i, j]` for the `i`th draw on the `j`th -##' stream. -##' -##' The output will not differ based on the number of threads used, -##' only on the number of streams. -##' -##' @return A `monty_rng` object, which can be used to drawn random -##' numbers from monty's distributions. -##' -##' @export -##' @examples -##' rng <- monty::monty_rng$new(42) -##' -##' # Shorthand for Uniform(0, 1) -##' rng$random_real(5) -##' -##' # Shorthand for Normal(0, 1) -##' rng$random_normal(5) -##' -##' # Uniform random numbers between min and max -##' rng$uniform(5, -2, 6) -##' -##' # Normally distributed random numbers with mean and sd -##' rng$normal(5, 4, 2) -##' -##' # Binomially distributed random numbers with size and prob -##' rng$binomial(5, 10, 0.3) -##' -##' # Negative binomially distributed random numbers with size and prob -##' rng$negative_binomial_prob(5, 10, 0.3) -##' -##' # Negative binomially distributed random numbers with size and mean mu -##' rng$negative_binomial_mu(5, 10, 25) -##' -##' # Hypergeometric distributed random numbers with parameters n1, n2 and k -##' rng$hypergeometric(5, 6, 10, 4) -##' -##' # Gamma distributed random numbers with parameters shape and scale -##' rng$gamma_scale(5, 0.5, 2) -##' -##' # Gamma distributed random numbers with parameters shape and rate -##' rng$gamma_rate(5, 0.5, 2) -##' -##' # Poisson distributed random numbers with mean lambda -##' rng$poisson(5, 2) -##' -##' # Exponentially distributed random numbers with rate -##' rng$exponential_rate(5, 2) -##' -##' # Exponentially distributed random numbers with mean -##' rng$exponential_mean(5, 0.5) -##' -##' # Multinomial distributed random numbers with size and vector of -##' # probabiltiies prob -##' rng$multinomial(5, 10, c(0.1, 0.3, 0.5, 0.1)) -monty_rng <- R6::R6Class( - "monty_rng", - cloneable = FALSE, - - private = list( - ptr = NULL, - n_streams = NULL - ), - - public = list( - ##' @field info Information about the generator (read-only) - info = NULL, - - ##' @description Create a `monty_rng` object - ##' - ##' @param seed The seed, as an integer, a raw vector or `NULL`. - ##' If an integer we will create a suitable seed via the "splitmix64" - ##' algorithm, if a raw vector it must the correct length (a multiple - ##' of 32). If `NULL` then we create a seed using R's random - ##' number generator. - ##' - ##' @param n_streams The number of streams to use (see Details) - ##' - ##' @param deterministic Logical, indicating if we should use - ##' "deterministic" mode where distributions return their - ##' expectations and the state is never changed. - initialize = function(seed = NULL, n_streams = 1L, - deterministic = FALSE) { - private$ptr <- monty_rng_alloc(seed, n_streams, deterministic) - private$n_streams <- n_streams - - size_int_bits <- 64L - name <- "xoshiro256plus" - size_int_bits <- 64L - self$info <- list( - int_type = sprintf("uint%s_t", size_int_bits), - name = name, - deterministic = deterministic, - ## Size, in bits, of the underlying integer - size_int_bits = size_int_bits, - ## Number of integers used for state - size_state_ints = 4L, - ## Total size in bytes of the state - size_state_bytes = 4L * size_int_bits / 8L) - lockBinding("info", self) - }, - - ##' @description Number of streams available - size = function() { - private$n_streams - }, - - ##' @description The jump function updates the random number state for - ##' each stream by advancing it to a state equivalent to - ##' 2^128 numbers drawn from each stream. - jump = function() { - monty_legacy_rng_jump(private$ptr) - invisible(self) - }, - - ##' @description Longer than `$jump`, the `$long_jump` method is - ##' equivalent to 2^192 numbers drawn from each stream. - long_jump = function() { - monty_legacy_rng_long_jump(private$ptr) - invisible(self) - }, - - ##' @description Generate `n` numbers from a standard uniform distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param n_threads Number of threads to use; see Details - random_real = function(n, n_threads = 1L) { - monty_rng_random_real(private$ptr, n, n_threads) - }, - - ##' @description Generate `n` numbers from a standard normal distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param n_threads Number of threads to use; see Details - ##' - ##' @param algorithm Name of the algorithm to use; currently `box_muller` - ##' and `ziggurat` are supported, with the latter being considerably - ##' faster. - random_normal = function(n, n_threads = 1L, algorithm = "box_muller") { - monty_rng_random_normal(private$ptr, n, n_threads, algorithm) - }, - - ##' @description Generate `n` numbers from a uniform distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param min The minimum of the distribution (length 1 or n) - ##' - ##' @param max The maximum of the distribution (length 1 or n) - ##' - ##' @param n_threads Number of threads to use; see Details - uniform = function(n, min, max, n_threads = 1L) { - monty_rng_uniform(private$ptr, n, min, max, n_threads) - }, - - ##' @description Generate `n` numbers from a normal distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param mean The mean of the distribution (length 1 or n) - ##' - ##' @param sd The standard deviation of the distribution (length 1 or n) - ##' - ##' @param n_threads Number of threads to use; see Details - ##' - ##' @param algorithm Name of the algorithm to use; currently `box_muller` - ##' and `ziggurat` are supported, with the latter being considerably - ##' faster. - normal = function(n, mean, sd, n_threads = 1L, algorithm = "box_muller") { - monty_rng_normal(private$ptr, n, mean, sd, n_threads, algorithm) - }, - - ##' @description Generate `n` numbers from a binomial distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param size The number of trials (zero or more, length 1 or n) - ##' - ##' @param prob The probability of success on each trial - ##' (between 0 and 1, length 1 or n) - ##' - ##' @param n_threads Number of threads to use; see Details - binomial = function(n, size, prob, n_threads = 1L) { - monty_rng_binomial(private$ptr, n, size, prob, n_threads) - }, - - ##' @description Generate `n` numbers from a beta-binomial distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param size The number of trials (zero or more, length 1 or n) - ##' - ##' @param a The first shape parameter (zero or more, length 1 or n) - ##' - ##' @param b The second shape parameter (zero or more, length 1 or n) - ##' - ##' @param n_threads Number of threads to use; see Details - beta_binomial_ab = function(n, size, a, b, n_threads = 1L) { - monty_rng_beta_binomial_ab(private$ptr, n, size, a, b, n_threads) - }, - - ##' @description Generate `n` numbers from a beta-binomial distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param size The number of trials (zero or more, length 1 or n) - ##' - ##' @param prob The mean probability of success on each trial - ##' (between 0 and 1, length 1 or n) - ##' - ##' @param rho The dispersion parameter (between 0 and 1, length 1 or n) - ##' - ##' @param n_threads Number of threads to use; see Details - beta_binomial_prob = function(n, size, prob, rho, n_threads = 1L) { - monty_rng_beta_binomial_prob(private$ptr, n, size, prob, rho, n_threads) - }, - - ##' @description Generate `n` numbers from a negative binomial distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param size The target number of successful trials - ##' (zero or more, length 1 or n) - ##' - ##' @param prob The probability of success on each trial - ##' (between 0 and 1, length 1 or n) - ##' - ##' @param n_threads Number of threads to use; see Details - negative_binomial_prob = function(n, size, prob, n_threads = 1L) { - monty_rng_negative_binomial_prob(private$ptr, n, size, prob, n_threads) - }, - - ##' @description Generate `n` numbers from a negative binomial distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param size The target number of successful trials - ##' (zero or more, length 1 or n) - ##' - ##' @param mu The mean - ##' (zero or more, length 1 or n) - ##' - ##' @param n_threads Number of threads to use; see Details - negative_binomial_mu = function(n, size, mu, n_threads = 1L) { - monty_rng_negative_binomial_mu(private$ptr, n, size, mu, n_threads) - }, - - ##' @description Generate `n` numbers from a hypergeometric distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param n1 The number of white balls in the urn (called n in - ##' R's [rhyper]) - ##' - ##' @param n2 The number of black balls in the urn (called m in - ##' R's [rhyper]) - ##' - ##' @param k The number of balls to draw - ##' - ##' @param n_threads Number of threads to use; see Details - hypergeometric = function(n, n1, n2, k, n_threads = 1L) { - monty_rng_hypergeometric(private$ptr, n, n1, n2, k, n_threads) - }, - - ##' @description Generate `n` numbers from a gamma distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param shape Shape - ##' - ##' @param scale Scale - ##'' - ##' @param n_threads Number of threads to use; see Details - gamma_scale = function(n, shape, scale, n_threads = 1L) { - monty_rng_gamma_scale(private$ptr, n, shape, scale, n_threads) - }, - - ##' @description Generate `n` numbers from a gamma distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param shape Shape - ##' - ##' @param rate Rate - ##'' - ##' @param n_threads Number of threads to use; see Details - gamma_rate = function(n, shape, rate, n_threads = 1L) { - monty_rng_gamma_rate(private$ptr, n, shape, rate, n_threads) - }, - - ##' @description Generate `n` numbers from a Poisson distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param lambda The mean (zero or more, length 1 or n). Only valid for - ##' lambda <= 10^7 - ##' - ##' @param n_threads Number of threads to use; see Details - poisson = function(n, lambda, n_threads = 1L) { - monty_rng_poisson(private$ptr, n, lambda, n_threads) - }, - - ##' @description Generate `n` numbers from a exponential distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param rate The rate of the exponential - ##' - ##' @param n_threads Number of threads to use; see Details - exponential_rate = function(n, rate, n_threads = 1L) { - monty_rng_exponential_rate(private$ptr, n, rate, n_threads) - }, - - ##' @description Generate `n` numbers from a exponential distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param mean The mean of the exponential - ##' - ##' @param n_threads Number of threads to use; see Details - exponential_mean = function(n, mean, n_threads = 1L) { - monty_rng_exponential_mean(private$ptr, n, mean, n_threads) - }, - - ##' @description Generate `n` draws from a Cauchy distribution. - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param location The location of the peak of the distribution - ##' (also its median) - ##' - ##' @param scale A scale parameter, which specifies the distribution's - ##' "half-width at half-maximum" - ##' - ##' @param n_threads Number of threads to use; see Details - cauchy = function(n, location, scale, n_threads = 1L) { - monty_rng_cauchy(private$ptr, n, location, scale, n_threads) - }, - - ##' @description Generate `n` draws from a multinomial distribution. - ##' In contrast with most of the distributions here, each draw is a - ##' *vector* with the same length as `prob`. - ##' - ##' @param n The number of samples to draw (per stream) - ##' - ##' @param size The number of trials (zero or more, length 1 or n) - ##' - ##' @param prob A vector of probabilities for the success of each - ##' trial. This does not need to sum to 1 (though all elements - ##' must be non-negative), in which case we interpret `prob` as - ##' weights and normalise so that they equal 1 before sampling. - ##' - ##' @param n_threads Number of threads to use; see Details - multinomial = function(n, size, prob, n_threads = 1L) { - monty_rng_multinomial(private$ptr, n, size, prob, n_threads) - }, - - ##' @description Generate `n` numbers from a beta distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param a The first shape parameter - ##' - ##' @param b The second shape parameter - ##' - ##' @param n_threads Number of threads to use; see Details - beta = function(n, a, b, n_threads = 1L) { - monty_rng_beta(private$ptr, n, a, b, n_threads) - }, - - ##' @description Generate `n` numbers from a truncated normal distribution - ##' - ##' @param n Number of samples to draw (per stream) - ##' - ##' @param mean The mean of the parent (untruncated) normal distribution - ##' - ##' @param sd The standard deviation of the parent (untruncated) - ##' normal distribution. - ##' - ##' @param min The lower bound - ##' - ##' @param max The upper bound - ##' - ##' @param n_threads Number of threads to use; see Details - truncated_normal = function(n, mean, sd, min, max, n_threads = 1L) { - monty_rng_truncated_normal(private$ptr, n, mean, sd, min, max, n_threads) - }, - - ##' @description - ##' Returns the state of the random number stream. This returns a - ##' raw vector of length 32 * `n_streams`. - state = function() { - monty_legacy_rng_state(private$ptr) - }, - - ##' @description - ##' Sets the state of the random number stream. - ##' @param state Raw vector of state, with length 32 * `n_streams`. - set_state = function(state) { - monty_legacy_rng_set_state(private$ptr, state) - } - )) diff --git a/R/rng_distributed.R b/R/rng_distributed.R deleted file mode 100644 index 04c61ecb..00000000 --- a/R/rng_distributed.R +++ /dev/null @@ -1,58 +0,0 @@ -##' Create a set of initial random number seeds suitable for using -##' within a distributed context (over multiple processes or nodes) at -##' a level higher than a single group of synchronised threads. -##' -##' See `vignette("rng_distributed")` for a proper introduction to -##' these functions. -##' -##' @title Create a set of distributed seeds -##' -##' @param seed Initial seed to use. As for [monty::monty_rng], this can -##' be `NULL` (create a seed using R's generators), an integer or a -##' raw vector of appropriate length. -##' -##' @param n_streams The number of streams to create per node. -##' -##' @param n_nodes The number of separate seeds to create. Each will -##' be separated by a "long jump" for your generator. -##' -##' @param algorithm The name of an algorithm to use. -##' -##' @return A list of either raw vectors (for -##' `monty_rng_distributed_state`) or of [monty::monty_rng_pointer] -##' objects (for `monty_rng_distributed_pointer`) -##' -##' @export -##' @rdname monty_rng_distributed -##' @examples -##' monty::monty_rng_distributed_state(n_nodes = 2) -##' monty::monty_rng_distributed_pointer(n_nodes = 2) -monty_rng_distributed_state <- function(seed = NULL, - n_streams = 1L, - n_nodes = 1L, - algorithm = "xoshiro256plus") { - p <- monty_rng_pointer$new(seed, n_streams, algorithm = algorithm) - - ret <- vector("list", n_nodes) - for (i in seq_len(n_nodes)) { - s <- p$state() - ret[[i]] <- s - if (i < n_nodes) { - p <- monty_rng_pointer$new(s, n_streams, 1L, algorithm = algorithm) - } - } - - ret -} - - -##' @export -##' @rdname monty_rng_distributed -monty_rng_distributed_pointer <- function(seed = NULL, - n_streams = 1L, - n_nodes = 1L, - algorithm = "xoshiro256plus") { - state <- monty_rng_distributed_state(seed, n_streams, n_nodes, algorithm) - lapply(state, monty_rng_pointer$new, - n_streams = n_streams, algorithm = algorithm) -} diff --git a/R/rng_pointer.R b/R/rng_pointer.R deleted file mode 100644 index a9a17bc9..00000000 --- a/R/rng_pointer.R +++ /dev/null @@ -1,81 +0,0 @@ -##' @title Create pointer to random number generator stream -##' -##' @description This function exists to support use from other -##' packages that wish to use monty's random number support, and -##' creates an opaque pointer to a set of random number streams. -##' -##' @export -##' @examples -##' monty::monty_rng_pointer$new() -monty_rng_pointer <- R6::R6Class( - "monty_rng_pointer", - cloneable = FALSE, - - private = list( - ptr_ = NULL, - state_ = NULL, - is_current_ = NULL - ), - - public = list( - ##' @field algorithm The name of the generator algorithm used (read-only) - algorithm = NULL, - - ##' @field n_streams The number of streams of random numbers provided - ##' (read-only) - n_streams = NULL, - - ##' @description Create a new `monty_rng_pointer` object - ##' - ##' @param seed The random number seed to use (see [monty::monty_rng] - ##' for details) - ##' - ##' @param n_streams The number of independent random number streams to - ##' create - ##' - ##' @param long_jump Optionally an integer indicating how many - ##' "long jumps" should be carried out immediately on creation. - ##' This can be used to create a distributed parallel random number - ##' generator (see [monty::monty_rng_distributed_state]) - ##' - ##' @param algorithm The random number algorithm to use. The default is - ##' `xoshiro256plus` which is a good general choice - initialize = function(seed = NULL, n_streams = 1L, long_jump = 0L, - algorithm = "xoshiro256plus") { - dat <- monty_rng_pointer_init(n_streams, seed, long_jump, algorithm) - private$ptr_ <- dat[[1L]] - private$state_ <- dat[[2L]] - private$is_current_ <- TRUE - - self$algorithm <- algorithm - self$n_streams <- n_streams - lockBinding("algorithm", self) - lockBinding("n_streams", self) - }, - - ##' @description Synchronise the R copy of the random number state. - ##' Typically this is only needed before serialisation if you have - ##' ever used the object. - sync = function() { - monty_rng_pointer_sync(private, self$algorithm) - invisible(self) - }, - - ##' @description Return a raw vector of state. This can be used to - ##' create other generators with the same state. - state = function() { - if (!private$is_current_) { - self$sync() - } - private$state_ - }, - - ##' @description Return a logical, indicating if the random number - ##' state that would be returned by `state()` is "current" (i.e., the - ##' same as the copy held in the pointer) or not. This is `TRUE` on - ##' creation or immediately after calling `$sync()` or `$state()` - ##' and `FALSE` after any use of the pointer. - is_current = function() { - private$is_current_ - } - )) diff --git a/_pkgdown.yml b/_pkgdown.yml index 8f10a95a..81731d32 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -119,14 +119,6 @@ reference: - monty_random_n_log_normal - monty_random_weibull - monty_random_n_weibull - - subsection: Legacy interface - desc: >- - Functions imported from dust; these will be retired soon - contents: - - monty_rng - - monty_rng_distributed_pointer - - monty_rng_distributed_state - - monty_rng_pointer articles: - title: Basics diff --git a/inst/include/monty/r/random.hpp b/inst/include/monty/r/random.hpp index c374f6f8..f8f641eb 100644 --- a/inst/include/monty/r/random.hpp +++ b/inst/include/monty/r/random.hpp @@ -109,84 +109,6 @@ cpp11::raws rng_state_vector(prng* rng) { } } - -template -SEXP rng_pointer_init(int n_streams, cpp11::sexp r_seed, int long_jump) { - auto seed = as_rng_seed(r_seed); - auto *rng = new prng(n_streams, seed); - for (int i = 0; i < long_jump; ++i) { - rng->long_jump(); - } - auto r_ptr = cpp11::external_pointer>(rng); - auto r_state = rng_state_vector(rng); - return cpp11::writable::list({r_ptr, r_state}); -} - -/// Recieve and check the pointer to rng state. This checks that the -/// object is valid, is of the correct state type, has sufficient -/// streams and has not been invalidated by serialisation. -/// -/// @tparam rng_state_type The random number state type to use -/// -/// @param obj A `monty_rng_pointer` object, created in R with -/// ``monty::monty_rng_pointer`` -/// -/// @param n_streams The number of required streams. Set this to 0 to -/// disable the check. If you are going to use 100 streams pass 100 -/// here, and a runtime error will be thrown if the object does not -/// contain enough streams, which is nicer than a crash when -/// `prng::state` fails. -template -prng* rng_pointer_get(cpp11::environment obj, - int n_streams = 0) { - // We could probably do this more efficiently if we store an enum - // in the object but this is probably ok. - const auto algorithm_given = cpp11::as_cpp(obj["algorithm"]); - const auto algorithm_expected = algorithm_name(); - if (algorithm_given != algorithm_expected) { - cpp11::stop("Incorrect rng type: given %s, expected %s", - algorithm_given.c_str(), algorithm_expected.c_str()); - } - - cpp11::environment env_enclos = - cpp11::as_cpp(obj[".__enclos_env__"]); - cpp11::environment env = - cpp11::as_cpp(env_enclos["private"]); - - using ptr_type = cpp11::external_pointer>; - auto ptr = cpp11::as_cpp(env["ptr_"]); - - auto * rng = ptr.get(); - if (rng == nullptr) { - if (!cpp11::as_cpp(env["is_current_"])) { - cpp11::stop("Can't unserialise an rng pointer that was not synced"); - } - cpp11::raws seed_data = cpp11::as_cpp(env["state_"]); - auto seed = raw_seed(seed_data); - const auto n_streams_orig = seed.size() / rng_state_type::size(); - rng = new prng(n_streams_orig, seed); - env["ptr_"] = cpp11::external_pointer>(rng); - } - - if (n_streams > 0 && static_cast(rng->size()) < n_streams) { - cpp11::stop("Requested a rng with %d streams but only have %d", - n_streams, rng->size()); - } - env["is_current_"] = cpp11::as_sexp(false); - - return rng; -} - -template -void rng_pointer_sync(cpp11::environment obj) { - using ptr_type = cpp11::external_pointer>; - if (!cpp11::as_cpp(obj["is_current_"])) { - auto ptr = cpp11::as_cpp(obj["ptr_"]); - obj["state_"] = rng_state_vector(ptr.get()); - obj["is_current_"] = cpp11::as_sexp(true); - } -} - } } } diff --git a/man/monty_dsl_parse_distribution.Rd b/man/monty_dsl_parse_distribution.Rd index ebd80470..9ffc5c2b 100644 --- a/man/monty_dsl_parse_distribution.Rd +++ b/man/monty_dsl_parse_distribution.Rd @@ -34,7 +34,7 @@ likely change once we support creation of differentiable models because we will want to do something with the arguments provided! \item \code{sample}: A function to sample from the distribution, given (as -a first argument) a rng object (see \link{monty_rng}) +a first argument) a \code{monty_rng} object (see \link{monty_rng_create}) } } \description{ diff --git a/man/monty_model.Rd b/man/monty_model.Rd index a553a0cd..9e759468 100644 --- a/man/monty_model.Rd +++ b/man/monty_model.Rd @@ -73,13 +73,13 @@ parameters. If named, then you can provide a subset, with parameters that are not included assumed to have a domain of \verb{(-Inf, Inf)}. \item \code{direct_sample}: A function to sample directly from the -parameter space, given a \link{monty_rng} object to sample from. -In the case where a model returns a posterior (e.g., in Bayesian -inference), this is assumed to be sampling from the prior. -We'll use this for generating initial conditions for MCMC where -those are not given, and possibly other uses. If not given then -when using \code{\link[=monty_sample]{monty_sample()}} the user will have to provide a -vector of initial states. +parameter space, given a \code{monty_rng} object to sample from (see +\link{monty_rng_create}). In the case where a model returns a +posterior (e.g., in Bayesian inference), this is assumed to be +sampling from the prior. We'll use this for generating initial +conditions for MCMC where those are not given, and possibly +other uses. If not given then when using \code{\link[=monty_sample]{monty_sample()}} the +user will have to provide a vector of initial states. \item \code{gradient}: A function to compute the gradient of \code{density} with respect to the parameter vector; takes a parameter vector and returns a vector the same length. For efficiency, the model may @@ -94,11 +94,11 @@ as that is the \emph{sampler's} rng stream, but we assume models will look after their own stream, and that they may need many streams). Models that provide this method are assumed to be stochastic; however, you can use the \code{is_stochastic} property -(via \code{\link[=monty_model_properties]{monty_model_properties()}}) to override this (e.g., to -run a stochastic model with its deterministic expectation). -This function takes a raw vector of random number state from -\link{monty_rng} and uses it to set the random number state for -your model; this is derived from the random number stream for a +(via \code{\link[=monty_model_properties]{monty_model_properties()}}) to override this (e.g., to run +a stochastic model with its deterministic expectation). This +function takes a raw vector of random number state from a +\code{monty_rng} and uses it to set the random number state for your +model; this is derived from the random number stream for a particular chain, jumped ahead. \item \code{get_rng_state}: A function to get the RNG state; must be provided if \code{set_rng_state} is present. Must return the random diff --git a/man/monty_model_direct_sample.Rd b/man/monty_model_direct_sample.Rd index 954e313e..27784c9c 100644 --- a/man/monty_model_direct_sample.Rd +++ b/man/monty_model_direct_sample.Rd @@ -9,8 +9,8 @@ monty_model_direct_sample(model, rng, named = FALSE) \arguments{ \item{model}{A \link{monty_model} object} -\item{rng}{Random number state, created by \link{monty_rng}. Use of -an RNG with more than one stream may or may not work as +\item{rng}{Random number state, created by \link{monty_rng_create}. +Use of an RNG with more than one stream may or may not work as expected; this is something we need to tidy up (\code{mrc-5292})} \item{named}{Logical, indicating if the output should be named diff --git a/man/monty_rng.Rd b/man/monty_rng.Rd deleted file mode 100644 index f00ef473..00000000 --- a/man/monty_rng.Rd +++ /dev/null @@ -1,715 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/rng.R -\name{monty_rng} -\alias{monty_rng} -\title{Monty Random Number Generator} -\value{ -A \code{monty_rng} object, which can be used to drawn random -numbers from monty's distributions. -} -\description{ -Create an object that can be used to generate random -numbers with the same RNG as monty uses internally. This is -primarily meant for debugging and testing the underlying C++ -rather than a source of random numbers from R. -} -\section{Warning}{ -This interface is subject to change in the near future, we do not -recommend its use in user code. -} - -\section{Running multiple streams, perhaps in parallel}{ -The underlying random number generators are designed to work in -parallel, and with random access to parameters (see -\code{vignette("rng")} for more details). However, this is usually -done within the context of running a model where each particle -sees its own stream of numbers. We provide some support for -running random number generators in parallel, but any speed -gains from parallelisation are likely to be somewhat eroded by -the overhead of copying around a large number of random numbers. - -All the random distribution functions support an argument -\code{n_threads} which controls the number of threads used. This -argument will \emph{silently} have no effect if your installation -does not support OpenMP. - -Parallelisation will be performed at the level of the stream, -where we draw \code{n} numbers from \emph{each} stream for a total of \code{n * n_streams} random numbers using \code{n_threads} threads to do this. -Setting \code{n_threads} to be higher than \code{n_streams} will therefore -have no effect. If running on somebody else's system (e.g., an -HPC, CRAN) you must respect the various environment variables -that control the maximum allowable number of threads. - -With the exception of \code{random_real}, each random number -distribution accepts parameters; the interpretations of these -will depend on \code{n}, \code{n_streams} and their rank. -\itemize{ -\item If a scalar then we will use the same parameter value for every draw -from every stream -\item If a vector with length \code{n} then we will draw \code{n} random -numbers per stream, and every stream will use the same parameter -value for every stream for each draw (but a different, -shared, parameter value for subsequent draws). -\item If a matrix is provided with one row and \code{n_streams} -columns then we use different parameters for each stream, but -the same parameter for each draw. -\item If a matrix is provided with \code{n} rows and \code{n_streams} -columns then we use a parameter value \verb{[i, j]} for the \code{i}th -draw on the \code{j}th stream. -} - -The rules are slightly different for the \code{prob} argument to -\code{multinomial} as for that \code{prob} is a vector of values. As such -we shift all dimensions by one: -\itemize{ -\item If a vector we use same \code{prob} every draw from every stream -and there are \code{length(prob)} possible outcomes. -\item If a matrix with \code{n} columns then vary over each draw (the -\code{i}th draw using vector \code{prob[, i]} but shared across all -streams. There are \code{nrow(prob)} possible outcomes. -\item If a 3d array is provided with 1 column and \code{n_streams} -"layers" (the third dimension) then we use then we use different -parameters for each stream, but the same parameter for each -draw. -\item If a 3d array is provided with \code{n} columns and \code{n_streams} -"layers" then we vary over both draws and streams so that with -use vector \code{prob[, i, j]} for the \code{i}th draw on the \code{j}th -stream. -} - -The output will not differ based on the number of threads used, -only on the number of streams. -} - -\examples{ -rng <- monty::monty_rng$new(42) - -# Shorthand for Uniform(0, 1) -rng$random_real(5) - -# Shorthand for Normal(0, 1) -rng$random_normal(5) - -# Uniform random numbers between min and max -rng$uniform(5, -2, 6) - -# Normally distributed random numbers with mean and sd -rng$normal(5, 4, 2) - -# Binomially distributed random numbers with size and prob -rng$binomial(5, 10, 0.3) - -# Negative binomially distributed random numbers with size and prob -rng$negative_binomial_prob(5, 10, 0.3) - -# Negative binomially distributed random numbers with size and mean mu -rng$negative_binomial_mu(5, 10, 25) - -# Hypergeometric distributed random numbers with parameters n1, n2 and k -rng$hypergeometric(5, 6, 10, 4) - -# Gamma distributed random numbers with parameters shape and scale -rng$gamma_scale(5, 0.5, 2) - -# Gamma distributed random numbers with parameters shape and rate -rng$gamma_rate(5, 0.5, 2) - -# Poisson distributed random numbers with mean lambda -rng$poisson(5, 2) - -# Exponentially distributed random numbers with rate -rng$exponential_rate(5, 2) - -# Exponentially distributed random numbers with mean -rng$exponential_mean(5, 0.5) - -# Multinomial distributed random numbers with size and vector of -# probabiltiies prob -rng$multinomial(5, 10, c(0.1, 0.3, 0.5, 0.1)) -} -\section{Public fields}{ -\if{html}{\out{
}} -\describe{ -\item{\code{info}}{Information about the generator (read-only)} -} -\if{html}{\out{
}} -} -\section{Methods}{ -\subsection{Public methods}{ -\itemize{ -\item \href{#method-monty_rng-new}{\code{monty_rng$new()}} -\item \href{#method-monty_rng-size}{\code{monty_rng$size()}} -\item \href{#method-monty_rng-jump}{\code{monty_rng$jump()}} -\item \href{#method-monty_rng-long_jump}{\code{monty_rng$long_jump()}} -\item \href{#method-monty_rng-random_real}{\code{monty_rng$random_real()}} -\item \href{#method-monty_rng-random_normal}{\code{monty_rng$random_normal()}} -\item \href{#method-monty_rng-uniform}{\code{monty_rng$uniform()}} -\item \href{#method-monty_rng-normal}{\code{monty_rng$normal()}} -\item \href{#method-monty_rng-binomial}{\code{monty_rng$binomial()}} -\item \href{#method-monty_rng-beta_binomial_ab}{\code{monty_rng$beta_binomial_ab()}} -\item \href{#method-monty_rng-beta_binomial_prob}{\code{monty_rng$beta_binomial_prob()}} -\item \href{#method-monty_rng-negative_binomial_prob}{\code{monty_rng$negative_binomial_prob()}} -\item \href{#method-monty_rng-negative_binomial_mu}{\code{monty_rng$negative_binomial_mu()}} -\item \href{#method-monty_rng-hypergeometric}{\code{monty_rng$hypergeometric()}} -\item \href{#method-monty_rng-gamma_scale}{\code{monty_rng$gamma_scale()}} -\item \href{#method-monty_rng-gamma_rate}{\code{monty_rng$gamma_rate()}} -\item \href{#method-monty_rng-poisson}{\code{monty_rng$poisson()}} -\item \href{#method-monty_rng-exponential_rate}{\code{monty_rng$exponential_rate()}} -\item \href{#method-monty_rng-exponential_mean}{\code{monty_rng$exponential_mean()}} -\item \href{#method-monty_rng-cauchy}{\code{monty_rng$cauchy()}} -\item \href{#method-monty_rng-multinomial}{\code{monty_rng$multinomial()}} -\item \href{#method-monty_rng-beta}{\code{monty_rng$beta()}} -\item \href{#method-monty_rng-truncated_normal}{\code{monty_rng$truncated_normal()}} -\item \href{#method-monty_rng-state}{\code{monty_rng$state()}} -\item \href{#method-monty_rng-set_state}{\code{monty_rng$set_state()}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-new}{}}} -\subsection{Method \code{new()}}{ -Create a \code{monty_rng} object -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$new(seed = NULL, n_streams = 1L, deterministic = FALSE)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{seed}}{The seed, as an integer, a raw vector or \code{NULL}. -If an integer we will create a suitable seed via the "splitmix64" -algorithm, if a raw vector it must the correct length (a multiple -of 32). If \code{NULL} then we create a seed using R's random -number generator.} - -\item{\code{n_streams}}{The number of streams to use (see Details)} - -\item{\code{deterministic}}{Logical, indicating if we should use -"deterministic" mode where distributions return their -expectations and the state is never changed.} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-size}{}}} -\subsection{Method \code{size()}}{ -Number of streams available -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$size()}\if{html}{\out{
}} -} - -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-jump}{}}} -\subsection{Method \code{jump()}}{ -The jump function updates the random number state for -each stream by advancing it to a state equivalent to -2^128 numbers drawn from each stream. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$jump()}\if{html}{\out{
}} -} - -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-long_jump}{}}} -\subsection{Method \code{long_jump()}}{ -Longer than \verb{$jump}, the \verb{$long_jump} method is -equivalent to 2^192 numbers drawn from each stream. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$long_jump()}\if{html}{\out{
}} -} - -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-random_real}{}}} -\subsection{Method \code{random_real()}}{ -Generate \code{n} numbers from a standard uniform distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$random_real(n, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-random_normal}{}}} -\subsection{Method \code{random_normal()}}{ -Generate \code{n} numbers from a standard normal distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$random_normal(n, n_threads = 1L, algorithm = "box_muller")}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{n_threads}}{Number of threads to use; see Details} - -\item{\code{algorithm}}{Name of the algorithm to use; currently \code{box_muller} -and \code{ziggurat} are supported, with the latter being considerably -faster.} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-uniform}{}}} -\subsection{Method \code{uniform()}}{ -Generate \code{n} numbers from a uniform distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$uniform(n, min, max, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{min}}{The minimum of the distribution (length 1 or n)} - -\item{\code{max}}{The maximum of the distribution (length 1 or n)} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-normal}{}}} -\subsection{Method \code{normal()}}{ -Generate \code{n} numbers from a normal distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$normal(n, mean, sd, n_threads = 1L, algorithm = "box_muller")}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{mean}}{The mean of the distribution (length 1 or n)} - -\item{\code{sd}}{The standard deviation of the distribution (length 1 or n)} - -\item{\code{n_threads}}{Number of threads to use; see Details} - -\item{\code{algorithm}}{Name of the algorithm to use; currently \code{box_muller} -and \code{ziggurat} are supported, with the latter being considerably -faster.} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-binomial}{}}} -\subsection{Method \code{binomial()}}{ -Generate \code{n} numbers from a binomial distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$binomial(n, size, prob, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{size}}{The number of trials (zero or more, length 1 or n)} - -\item{\code{prob}}{The probability of success on each trial -(between 0 and 1, length 1 or n)} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-beta_binomial_ab}{}}} -\subsection{Method \code{beta_binomial_ab()}}{ -Generate \code{n} numbers from a beta-binomial distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$beta_binomial_ab(n, size, a, b, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{size}}{The number of trials (zero or more, length 1 or n)} - -\item{\code{a}}{The first shape parameter (zero or more, length 1 or n)} - -\item{\code{b}}{The second shape parameter (zero or more, length 1 or n)} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-beta_binomial_prob}{}}} -\subsection{Method \code{beta_binomial_prob()}}{ -Generate \code{n} numbers from a beta-binomial distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$beta_binomial_prob(n, size, prob, rho, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{size}}{The number of trials (zero or more, length 1 or n)} - -\item{\code{prob}}{The mean probability of success on each trial -(between 0 and 1, length 1 or n)} - -\item{\code{rho}}{The dispersion parameter (between 0 and 1, length 1 or n)} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-negative_binomial_prob}{}}} -\subsection{Method \code{negative_binomial_prob()}}{ -Generate \code{n} numbers from a negative binomial distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$negative_binomial_prob(n, size, prob, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{size}}{The target number of successful trials -(zero or more, length 1 or n)} - -\item{\code{prob}}{The probability of success on each trial -(between 0 and 1, length 1 or n)} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-negative_binomial_mu}{}}} -\subsection{Method \code{negative_binomial_mu()}}{ -Generate \code{n} numbers from a negative binomial distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$negative_binomial_mu(n, size, mu, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{size}}{The target number of successful trials -(zero or more, length 1 or n)} - -\item{\code{mu}}{The mean -(zero or more, length 1 or n)} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-hypergeometric}{}}} -\subsection{Method \code{hypergeometric()}}{ -Generate \code{n} numbers from a hypergeometric distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$hypergeometric(n, n1, n2, k, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{n1}}{The number of white balls in the urn (called n in -R's \link{rhyper})} - -\item{\code{n2}}{The number of black balls in the urn (called m in -R's \link{rhyper})} - -\item{\code{k}}{The number of balls to draw} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-gamma_scale}{}}} -\subsection{Method \code{gamma_scale()}}{ -Generate \code{n} numbers from a gamma distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$gamma_scale(n, shape, scale, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{shape}}{Shape} - -\item{\code{scale}}{Scale -'} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-gamma_rate}{}}} -\subsection{Method \code{gamma_rate()}}{ -Generate \code{n} numbers from a gamma distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$gamma_rate(n, shape, rate, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{shape}}{Shape} - -\item{\code{rate}}{Rate -'} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-poisson}{}}} -\subsection{Method \code{poisson()}}{ -Generate \code{n} numbers from a Poisson distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$poisson(n, lambda, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{lambda}}{The mean (zero or more, length 1 or n). Only valid for -lambda <= 10^7} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-exponential_rate}{}}} -\subsection{Method \code{exponential_rate()}}{ -Generate \code{n} numbers from a exponential distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$exponential_rate(n, rate, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{rate}}{The rate of the exponential} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-exponential_mean}{}}} -\subsection{Method \code{exponential_mean()}}{ -Generate \code{n} numbers from a exponential distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$exponential_mean(n, mean, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{mean}}{The mean of the exponential} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-cauchy}{}}} -\subsection{Method \code{cauchy()}}{ -Generate \code{n} draws from a Cauchy distribution. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$cauchy(n, location, scale, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{location}}{The location of the peak of the distribution -(also its median)} - -\item{\code{scale}}{A scale parameter, which specifies the distribution's -"half-width at half-maximum"} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-multinomial}{}}} -\subsection{Method \code{multinomial()}}{ -Generate \code{n} draws from a multinomial distribution. -In contrast with most of the distributions here, each draw is a -\emph{vector} with the same length as \code{prob}. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$multinomial(n, size, prob, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{The number of samples to draw (per stream)} - -\item{\code{size}}{The number of trials (zero or more, length 1 or n)} - -\item{\code{prob}}{A vector of probabilities for the success of each -trial. This does not need to sum to 1 (though all elements -must be non-negative), in which case we interpret \code{prob} as -weights and normalise so that they equal 1 before sampling.} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-beta}{}}} -\subsection{Method \code{beta()}}{ -Generate \code{n} numbers from a beta distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$beta(n, a, b, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{a}}{The first shape parameter} - -\item{\code{b}}{The second shape parameter} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-truncated_normal}{}}} -\subsection{Method \code{truncated_normal()}}{ -Generate \code{n} numbers from a truncated normal distribution -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$truncated_normal(n, mean, sd, min, max, n_threads = 1L)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{n}}{Number of samples to draw (per stream)} - -\item{\code{mean}}{The mean of the parent (untruncated) normal distribution} - -\item{\code{sd}}{The standard deviation of the parent (untruncated) -normal distribution.} - -\item{\code{min}}{The lower bound} - -\item{\code{max}}{The upper bound} - -\item{\code{n_threads}}{Number of threads to use; see Details} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-state}{}}} -\subsection{Method \code{state()}}{ -Returns the state of the random number stream. This returns a -raw vector of length 32 * \code{n_streams}. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$state()}\if{html}{\out{
}} -} - -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng-set_state}{}}} -\subsection{Method \code{set_state()}}{ -Sets the state of the random number stream. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng$set_state(state)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{state}}{Raw vector of state, with length 32 * \code{n_streams}.} -} -\if{html}{\out{
}} -} -} -} diff --git a/man/monty_rng_distributed.Rd b/man/monty_rng_distributed.Rd deleted file mode 100644 index fda963f7..00000000 --- a/man/monty_rng_distributed.Rd +++ /dev/null @@ -1,51 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/rng_distributed.R -\name{monty_rng_distributed_state} -\alias{monty_rng_distributed_state} -\alias{monty_rng_distributed_pointer} -\title{Create a set of distributed seeds} -\usage{ -monty_rng_distributed_state( - seed = NULL, - n_streams = 1L, - n_nodes = 1L, - algorithm = "xoshiro256plus" -) - -monty_rng_distributed_pointer( - seed = NULL, - n_streams = 1L, - n_nodes = 1L, - algorithm = "xoshiro256plus" -) -} -\arguments{ -\item{seed}{Initial seed to use. As for \link{monty_rng}, this can -be \code{NULL} (create a seed using R's generators), an integer or a -raw vector of appropriate length.} - -\item{n_streams}{The number of streams to create per node.} - -\item{n_nodes}{The number of separate seeds to create. Each will -be separated by a "long jump" for your generator.} - -\item{algorithm}{The name of an algorithm to use.} -} -\value{ -A list of either raw vectors (for -\code{monty_rng_distributed_state}) or of \link{monty_rng_pointer} -objects (for \code{monty_rng_distributed_pointer}) -} -\description{ -Create a set of initial random number seeds suitable for using -within a distributed context (over multiple processes or nodes) at -a level higher than a single group of synchronised threads. -} -\details{ -See \code{vignette("rng_distributed")} for a proper introduction to -these functions. -} -\examples{ -monty::monty_rng_distributed_state(n_nodes = 2) -monty::monty_rng_distributed_pointer(n_nodes = 2) -} diff --git a/man/monty_rng_pointer.Rd b/man/monty_rng_pointer.Rd deleted file mode 100644 index f7fef23a..00000000 --- a/man/monty_rng_pointer.Rd +++ /dev/null @@ -1,104 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/rng_pointer.R -\name{monty_rng_pointer} -\alias{monty_rng_pointer} -\title{Create pointer to random number generator stream} -\description{ -This function exists to support use from other -packages that wish to use monty's random number support, and -creates an opaque pointer to a set of random number streams. -} -\examples{ -monty::monty_rng_pointer$new() -} -\section{Public fields}{ -\if{html}{\out{
}} -\describe{ -\item{\code{algorithm}}{The name of the generator algorithm used (read-only)} - -\item{\code{n_streams}}{The number of streams of random numbers provided -(read-only)} -} -\if{html}{\out{
}} -} -\section{Methods}{ -\subsection{Public methods}{ -\itemize{ -\item \href{#method-monty_rng_pointer-new}{\code{monty_rng_pointer$new()}} -\item \href{#method-monty_rng_pointer-sync}{\code{monty_rng_pointer$sync()}} -\item \href{#method-monty_rng_pointer-state}{\code{monty_rng_pointer$state()}} -\item \href{#method-monty_rng_pointer-is_current}{\code{monty_rng_pointer$is_current()}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng_pointer-new}{}}} -\subsection{Method \code{new()}}{ -Create a new \code{monty_rng_pointer} object -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng_pointer$new( - seed = NULL, - n_streams = 1L, - long_jump = 0L, - algorithm = "xoshiro256plus" -)}\if{html}{\out{
}} -} - -\subsection{Arguments}{ -\if{html}{\out{
}} -\describe{ -\item{\code{seed}}{The random number seed to use (see \link{monty_rng} -for details)} - -\item{\code{n_streams}}{The number of independent random number streams to -create} - -\item{\code{long_jump}}{Optionally an integer indicating how many -"long jumps" should be carried out immediately on creation. -This can be used to create a distributed parallel random number -generator (see \link{monty_rng_distributed_state})} - -\item{\code{algorithm}}{The random number algorithm to use. The default is -\code{xoshiro256plus} which is a good general choice} -} -\if{html}{\out{
}} -} -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng_pointer-sync}{}}} -\subsection{Method \code{sync()}}{ -Synchronise the R copy of the random number state. -Typically this is only needed before serialisation if you have -ever used the object. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng_pointer$sync()}\if{html}{\out{
}} -} - -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng_pointer-state}{}}} -\subsection{Method \code{state()}}{ -Return a raw vector of state. This can be used to -create other generators with the same state. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng_pointer$state()}\if{html}{\out{
}} -} - -} -\if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-monty_rng_pointer-is_current}{}}} -\subsection{Method \code{is_current()}}{ -Return a logical, indicating if the random number -state that would be returned by \code{state()} is "current" (i.e., the -same as the copy held in the pointer) or not. This is \code{TRUE} on -creation or immediately after calling \verb{$sync()} or \verb{$state()} -and \code{FALSE} after any use of the pointer. -\subsection{Usage}{ -\if{html}{\out{
}}\preformatted{monty_rng_pointer$is_current()}\if{html}{\out{
}} -} - -} -} diff --git a/src/cpp11.cpp b/src/cpp11.cpp index e1f8814a..6348c5cb 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -5,6 +5,13 @@ #include "cpp11/declarations.hpp" #include +// random.cpp +SEXP monty_rng_alloc(cpp11::sexp r_seed, int n_streams, bool deterministic); +extern "C" SEXP _monty_monty_rng_alloc(SEXP r_seed, SEXP n_streams, SEXP deterministic) { + BEGIN_CPP11 + return cpp11::as_sexp(monty_rng_alloc(cpp11::as_cpp>(r_seed), cpp11::as_cpp>(n_streams), cpp11::as_cpp>(deterministic))); + END_CPP11 +} // random.cpp cpp11::sexp cpp_monty_rng_state(cpp11::sexp ptr); extern "C" SEXP _monty_cpp_monty_rng_state(SEXP ptr) { @@ -344,204 +351,11 @@ extern "C" SEXP _monty_cpp_monty_random_n_multinomial(SEXP n_samples, SEXP r_siz return cpp11::as_sexp(cpp_monty_random_n_multinomial(cpp11::as_cpp>(n_samples), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(ptr))); END_CPP11 } -// random_legacy.cpp -SEXP monty_rng_alloc(cpp11::sexp r_seed, int n_streams, bool deterministic); -extern "C" SEXP _monty_monty_rng_alloc(SEXP r_seed, SEXP n_streams, SEXP deterministic) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_alloc(cpp11::as_cpp>(r_seed), cpp11::as_cpp>(n_streams), cpp11::as_cpp>(deterministic))); - END_CPP11 -} -// random_legacy.cpp -void monty_legacy_rng_jump(SEXP ptr); -extern "C" SEXP _monty_monty_legacy_rng_jump(SEXP ptr) { - BEGIN_CPP11 - monty_legacy_rng_jump(cpp11::as_cpp>(ptr)); - return R_NilValue; - END_CPP11 -} -// random_legacy.cpp -void monty_legacy_rng_long_jump(SEXP ptr); -extern "C" SEXP _monty_monty_legacy_rng_long_jump(SEXP ptr) { - BEGIN_CPP11 - monty_legacy_rng_long_jump(cpp11::as_cpp>(ptr)); - return R_NilValue; - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_random_real(SEXP ptr, int n, int n_threads); -extern "C" SEXP _monty_monty_rng_random_real(SEXP ptr, SEXP n, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_random_real(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_random_normal(SEXP ptr, int n, int n_threads, std::string algorithm); -extern "C" SEXP _monty_monty_rng_random_normal(SEXP ptr, SEXP n, SEXP n_threads, SEXP algorithm) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_random_normal(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(algorithm))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_uniform(SEXP ptr, int n, cpp11::doubles r_min, cpp11::doubles r_max, int n_threads); -extern "C" SEXP _monty_monty_rng_uniform(SEXP ptr, SEXP n, SEXP r_min, SEXP r_max, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_uniform(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_min), cpp11::as_cpp>(r_max), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_exponential_rate(SEXP ptr, int n, cpp11::doubles r_rate, int n_threads); -extern "C" SEXP _monty_monty_rng_exponential_rate(SEXP ptr, SEXP n, SEXP r_rate, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_exponential_rate(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_rate), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_exponential_mean(SEXP ptr, int n, cpp11::doubles r_mean, int n_threads); -extern "C" SEXP _monty_monty_rng_exponential_mean(SEXP ptr, SEXP n, SEXP r_mean, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_exponential_mean(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_mean), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_normal(SEXP ptr, int n, cpp11::doubles r_mean, cpp11::doubles r_sd, int n_threads, std::string algorithm); -extern "C" SEXP _monty_monty_rng_normal(SEXP ptr, SEXP n, SEXP r_mean, SEXP r_sd, SEXP n_threads, SEXP algorithm) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_normal(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_mean), cpp11::as_cpp>(r_sd), cpp11::as_cpp>(n_threads), cpp11::as_cpp>(algorithm))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_binomial(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads); -extern "C" SEXP _monty_monty_rng_binomial(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_binomial(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_beta_binomial_ab(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_a, cpp11::doubles r_b, int n_threads); -extern "C" SEXP _monty_monty_rng_beta_binomial_ab(SEXP ptr, SEXP n, SEXP r_size, SEXP r_a, SEXP r_b, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_beta_binomial_ab(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_a), cpp11::as_cpp>(r_b), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_beta_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, cpp11::doubles r_rho, int n_threads); -extern "C" SEXP _monty_monty_rng_beta_binomial_prob(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP r_rho, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_beta_binomial_prob(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(r_rho), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads); -extern "C" SEXP _monty_monty_rng_negative_binomial_prob(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_negative_binomial_prob(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_mu, int n_threads); -extern "C" SEXP _monty_monty_rng_negative_binomial_mu(SEXP ptr, SEXP n, SEXP r_size, SEXP r_mu, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_negative_binomial_mu(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_mu), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_hypergeometric(SEXP ptr, int n, cpp11::doubles r_n1, cpp11::doubles r_n2, cpp11::doubles r_k, int n_threads); -extern "C" SEXP _monty_monty_rng_hypergeometric(SEXP ptr, SEXP n, SEXP r_n1, SEXP r_n2, SEXP r_k, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_hypergeometric(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_n1), cpp11::as_cpp>(r_n2), cpp11::as_cpp>(r_k), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_gamma_scale(SEXP ptr, int n, cpp11::doubles r_shape, cpp11::doubles r_scale, int n_threads); -extern "C" SEXP _monty_monty_rng_gamma_scale(SEXP ptr, SEXP n, SEXP r_shape, SEXP r_scale, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_gamma_scale(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_shape), cpp11::as_cpp>(r_scale), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_gamma_rate(SEXP ptr, int n, cpp11::doubles r_shape, cpp11::doubles r_rate, int n_threads); -extern "C" SEXP _monty_monty_rng_gamma_rate(SEXP ptr, SEXP n, SEXP r_shape, SEXP r_rate, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_gamma_rate(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_shape), cpp11::as_cpp>(r_rate), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_poisson(SEXP ptr, int n, cpp11::doubles r_lambda, int n_threads); -extern "C" SEXP _monty_monty_rng_poisson(SEXP ptr, SEXP n, SEXP r_lambda, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_poisson(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_lambda), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_cauchy(SEXP ptr, int n, cpp11::doubles r_location, cpp11::doubles r_scale, int n_threads); -extern "C" SEXP _monty_monty_rng_cauchy(SEXP ptr, SEXP n, SEXP r_location, SEXP r_scale, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_cauchy(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_location), cpp11::as_cpp>(r_scale), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_beta(SEXP ptr, int n, cpp11::doubles r_a, cpp11::doubles r_b, int n_threads); -extern "C" SEXP _monty_monty_rng_beta(SEXP ptr, SEXP n, SEXP r_a, SEXP r_b, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_beta(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_a), cpp11::as_cpp>(r_b), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_multinomial(SEXP ptr, int n, cpp11::doubles r_size, cpp11::doubles r_prob, int n_threads); -extern "C" SEXP _monty_monty_rng_multinomial(SEXP ptr, SEXP n, SEXP r_size, SEXP r_prob, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_multinomial(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_rng_truncated_normal(SEXP ptr, int n, cpp11::doubles r_mean, cpp11::doubles r_sd, cpp11::doubles r_min, cpp11::doubles r_max, int n_threads); -extern "C" SEXP _monty_monty_rng_truncated_normal(SEXP ptr, SEXP n, SEXP r_mean, SEXP r_sd, SEXP r_min, SEXP r_max, SEXP n_threads) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_truncated_normal(cpp11::as_cpp>(ptr), cpp11::as_cpp>(n), cpp11::as_cpp>(r_mean), cpp11::as_cpp>(r_sd), cpp11::as_cpp>(r_min), cpp11::as_cpp>(r_max), cpp11::as_cpp>(n_threads))); - END_CPP11 -} -// random_legacy.cpp -cpp11::sexp monty_legacy_rng_state(SEXP ptr); -extern "C" SEXP _monty_monty_legacy_rng_state(SEXP ptr) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_legacy_rng_state(cpp11::as_cpp>(ptr))); - END_CPP11 -} -// random_legacy.cpp -void monty_legacy_rng_set_state(SEXP ptr, cpp11::raws r_state); -extern "C" SEXP _monty_monty_legacy_rng_set_state(SEXP ptr, SEXP r_state) { - BEGIN_CPP11 - monty_legacy_rng_set_state(cpp11::as_cpp>(ptr), cpp11::as_cpp>(r_state)); - return R_NilValue; - END_CPP11 -} -// rng_pointer.cpp -cpp11::sexp monty_rng_pointer_init(int n_streams, cpp11::sexp seed, int long_jump, std::string algorithm); -extern "C" SEXP _monty_monty_rng_pointer_init(SEXP n_streams, SEXP seed, SEXP long_jump, SEXP algorithm) { - BEGIN_CPP11 - return cpp11::as_sexp(monty_rng_pointer_init(cpp11::as_cpp>(n_streams), cpp11::as_cpp>(seed), cpp11::as_cpp>(long_jump), cpp11::as_cpp>(algorithm))); - END_CPP11 -} -// rng_pointer.cpp -void monty_rng_pointer_sync(cpp11::environment obj, std::string algorithm); -extern "C" SEXP _monty_monty_rng_pointer_sync(SEXP obj, SEXP algorithm) { - BEGIN_CPP11 - monty_rng_pointer_sync(cpp11::as_cpp>(obj), cpp11::as_cpp>(algorithm)); - return R_NilValue; - END_CPP11 -} -// rng_pointer.cpp -double test_rng_pointer_get(cpp11::environment obj, int n_streams); -extern "C" SEXP _monty_test_rng_pointer_get(SEXP obj, SEXP n_streams) { - BEGIN_CPP11 - return cpp11::as_sexp(test_rng_pointer_get(cpp11::as_cpp>(obj), cpp11::as_cpp>(n_streams))); - END_CPP11 -} // test_rng.cpp -std::vector test_xoshiro_run(cpp11::environment obj); -extern "C" SEXP _monty_test_xoshiro_run(SEXP obj) { +std::vector test_xoshiro_run(cpp11::sexp ptr); +extern "C" SEXP _monty_test_xoshiro_run(SEXP ptr) { BEGIN_CPP11 - return cpp11::as_sexp(test_xoshiro_run(cpp11::as_cpp>(obj))); + return cpp11::as_sexp(test_xoshiro_run(cpp11::as_cpp>(ptr))); END_CPP11 } @@ -595,33 +409,7 @@ static const R_CallMethodDef CallEntries[] = { {"_monty_cpp_monty_rng_long_jump", (DL_FUNC) &_monty_cpp_monty_rng_long_jump, 2}, {"_monty_cpp_monty_rng_set_state", (DL_FUNC) &_monty_cpp_monty_rng_set_state, 2}, {"_monty_cpp_monty_rng_state", (DL_FUNC) &_monty_cpp_monty_rng_state, 1}, - {"_monty_monty_legacy_rng_jump", (DL_FUNC) &_monty_monty_legacy_rng_jump, 1}, - {"_monty_monty_legacy_rng_long_jump", (DL_FUNC) &_monty_monty_legacy_rng_long_jump, 1}, - {"_monty_monty_legacy_rng_set_state", (DL_FUNC) &_monty_monty_legacy_rng_set_state, 2}, - {"_monty_monty_legacy_rng_state", (DL_FUNC) &_monty_monty_legacy_rng_state, 1}, {"_monty_monty_rng_alloc", (DL_FUNC) &_monty_monty_rng_alloc, 3}, - {"_monty_monty_rng_beta", (DL_FUNC) &_monty_monty_rng_beta, 5}, - {"_monty_monty_rng_beta_binomial_ab", (DL_FUNC) &_monty_monty_rng_beta_binomial_ab, 6}, - {"_monty_monty_rng_beta_binomial_prob", (DL_FUNC) &_monty_monty_rng_beta_binomial_prob, 6}, - {"_monty_monty_rng_binomial", (DL_FUNC) &_monty_monty_rng_binomial, 5}, - {"_monty_monty_rng_cauchy", (DL_FUNC) &_monty_monty_rng_cauchy, 5}, - {"_monty_monty_rng_exponential_mean", (DL_FUNC) &_monty_monty_rng_exponential_mean, 4}, - {"_monty_monty_rng_exponential_rate", (DL_FUNC) &_monty_monty_rng_exponential_rate, 4}, - {"_monty_monty_rng_gamma_rate", (DL_FUNC) &_monty_monty_rng_gamma_rate, 5}, - {"_monty_monty_rng_gamma_scale", (DL_FUNC) &_monty_monty_rng_gamma_scale, 5}, - {"_monty_monty_rng_hypergeometric", (DL_FUNC) &_monty_monty_rng_hypergeometric, 6}, - {"_monty_monty_rng_multinomial", (DL_FUNC) &_monty_monty_rng_multinomial, 5}, - {"_monty_monty_rng_negative_binomial_mu", (DL_FUNC) &_monty_monty_rng_negative_binomial_mu, 5}, - {"_monty_monty_rng_negative_binomial_prob", (DL_FUNC) &_monty_monty_rng_negative_binomial_prob, 5}, - {"_monty_monty_rng_normal", (DL_FUNC) &_monty_monty_rng_normal, 6}, - {"_monty_monty_rng_pointer_init", (DL_FUNC) &_monty_monty_rng_pointer_init, 4}, - {"_monty_monty_rng_pointer_sync", (DL_FUNC) &_monty_monty_rng_pointer_sync, 2}, - {"_monty_monty_rng_poisson", (DL_FUNC) &_monty_monty_rng_poisson, 4}, - {"_monty_monty_rng_random_normal", (DL_FUNC) &_monty_monty_rng_random_normal, 4}, - {"_monty_monty_rng_random_real", (DL_FUNC) &_monty_monty_rng_random_real, 3}, - {"_monty_monty_rng_truncated_normal", (DL_FUNC) &_monty_monty_rng_truncated_normal, 7}, - {"_monty_monty_rng_uniform", (DL_FUNC) &_monty_monty_rng_uniform, 5}, - {"_monty_test_rng_pointer_get", (DL_FUNC) &_monty_test_rng_pointer_get, 2}, {"_monty_test_xoshiro_run", (DL_FUNC) &_monty_test_xoshiro_run, 1}, {NULL, NULL, 0} }; diff --git a/src/random.cpp b/src/random.cpp index c7a81e3d..a8611ed1 100644 --- a/src/random.cpp +++ b/src/random.cpp @@ -337,7 +337,22 @@ cpp11::doubles monty_random_sample_n_4(Fn fn, size_t n_samples, return r_y; } + +template +SEXP monty_rng_alloc(cpp11::sexp r_seed, int n_streams, bool deterministic) { + auto seed = monty::random::r::as_rng_seed(r_seed); + T *rng = new T(n_streams, seed, deterministic); + return cpp11::external_pointer(rng); +} + + // Real functions that we export +[[cpp11::register]] +SEXP monty_rng_alloc(cpp11::sexp r_seed, int n_streams, bool deterministic) { + return monty_rng_alloc(r_seed, n_streams, deterministic); +} + + [[cpp11::register]] cpp11::sexp cpp_monty_rng_state(cpp11::sexp ptr) { auto * rng = safely_read_externalptr(ptr, "state"); diff --git a/src/random_legacy.cpp b/src/random_legacy.cpp deleted file mode 100644 index 53127b87..00000000 --- a/src/random_legacy.cpp +++ /dev/null @@ -1,1091 +0,0 @@ -#include - -#ifdef _OPENMP -#include -#endif - -#include -#include -#include -#include - -#include -#include -#include - -using default_rng = monty::random::prng>; - -template -T* safely_read_externalptr(SEXP ptr, const char * context) { - if (!R_ExternalPtrAddr(ptr)) { - cpp11::stop("Pointer has been serialised, cannot continue safely (%s)", - context); - } - return cpp11::as_cpp>(ptr).get(); -} - -template -SEXP monty_rng_alloc(cpp11::sexp r_seed, int n_streams, bool deterministic) { - auto seed = monty::random::r::as_rng_seed(r_seed); - T *rng = new T(n_streams, seed, deterministic); - return cpp11::external_pointer(rng); -} - -template -void monty_legacy_rng_jump(SEXP ptr) { - T *rng = safely_read_externalptr(ptr, "jump"); - rng->jump(); -} - -template -void monty_legacy_rng_long_jump(SEXP ptr) { - T *rng = safely_read_externalptr(ptr, "long_jump"); - rng->long_jump(); -} - -// Little helper for returning x as a vector (m == 1) or matrix (n * -// m) by setting the dimension attribute. -cpp11::sexp sexp_matrix(cpp11::sexp x, int n, int m) { - if (m > 1) { - x.attr("dim") = cpp11::writable::integers{n, m}; - } - return x; -} - -template -cpp11::sexp monty_rng_random_real(SEXP ptr, int n, int n_threads) { - T *rng = safely_read_externalptr(ptr, "random_real"); - const int n_streams = rng->size(); - - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - auto &state = rng->state(i); - auto y_i = y + n * i; - for (size_t j = 0; j < (size_t)n; ++j) { - y_i[j] = monty::random::random_real(state); - } - } - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_random_normal(SEXP ptr, int n, int n_threads) { - T *rng = safely_read_externalptr(ptr, "random_normal"); - const int n_streams = rng->size(); - - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - auto &state = rng->state(i); - auto y_i = y + n * i; - for (size_t j = 0; j < (size_t)n; ++j) { - y_i[j] = monty::random::random_normal(state); - } - } - - return sexp_matrix(ret, n, n_streams); -} - -struct input_vary { - size_t len; - size_t offset; - bool draw; - bool generator; -}; - -// See notes in R/rng.R or ?rng -input_vary check_input_type(cpp11::doubles x, int n, int m, const char *name) { - input_vary ret {1, 1, false, false}; - if (Rf_isMatrix(x)) { - if (Rf_ncols(x) != m) { - cpp11::stop("If '%s' is a matrix, it must have %d columns", name, m); - } - ret.generator = true; - if (Rf_nrows(x) == n) { - ret.draw = true; - } else if (Rf_nrows(x) != 1) { - cpp11::stop("If '%s' is a matrix, it must have 1 or %d rows", name, n); - } - } else { - if (x.size() == n) { - ret.draw = true; - } else if (x.size() != 1) { - cpp11::stop("If '%s' is a vector, it must have 1 or %d elements", - name, n); - } - } - - if (ret.draw) { - ret.offset = n; - } - - return ret; -} - -// See notes in R/rng.R or ?rng -input_vary check_input_type2(cpp11::doubles x, int n, int m, const char *name) { - input_vary ret {1, 1, false, false}; - cpp11::sexp r_dim = x.attr("dim"); - if (r_dim == R_NilValue) { - ret.len = x.size(); - } else if (LENGTH(r_dim) == 2) { // matrix - auto dim = cpp11::as_cpp(r_dim); - ret.len = dim[0]; - if (dim[1] == n) { - ret.draw = true; - } else { - // TODO: must be n, not 1 surely? - cpp11::stop("If '%s' is a matrix, it must have %d columns", name, n); - } - } else if (LENGTH(r_dim) == 3) { - auto dim = cpp11::as_cpp(r_dim); - ret.len = dim[0]; - if (dim[1] == n) { - ret.draw = true; - } else if (dim[1] != 1) { - cpp11::stop("If '%s' is a 3d array, it must have 1 or %d columns", - name, n); - } - if (dim[2] != m) { - cpp11::stop("If '%s' is a 3d array, it must have %d layers", name, m); - } - ret.generator = true; - } else { - cpp11::stop("'%s' must be a vector, matrix or 3d array", name); - } - - if (ret.len < 2) { - cpp11::stop("Input parameters imply length of '%s' of only %d (< 2)", - name, ret.len); - } - - if (ret.draw) { - ret.offset = n * ret.len; - } else { - ret.offset = ret.len; - } - - return ret; -} - -// Below here is very repetitive, and could probably be deduplicated -// with some clever template magic. Most of the faff is because we -// want to support 4 modes of taking 1 or 2 parameters (each varying -// or not over draws and generators) -template -cpp11::sexp monty_rng_uniform(SEXP ptr, int n, - cpp11::doubles r_min, - cpp11::doubles r_max, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "uniform"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * min = REAL(r_min); - const double * max = REAL(r_max); - auto min_vary = check_input_type(r_min, n, n_streams, "min"); - auto max_vary = check_input_type(r_max, n, n_streams, "max"); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto min_i = min_vary.generator ? min + min_vary.offset * i : min; - auto max_i = max_vary.generator ? max + max_vary.offset * i : max; - for (size_t j = 0; j < (size_t)n; ++j) { - auto min_ij = min_vary.draw ? min_i[j] : min_i[0]; - auto max_ij = max_vary.draw ? max_i[j] : max_i[0]; - y_i[j] = monty::random::uniform(state, min_ij, max_ij); - } - } - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_exponential_rate(SEXP ptr, int n, cpp11::doubles r_rate, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "exponential_rate"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * rate = REAL(r_rate); - auto rate_vary = check_input_type(r_rate, n, n_streams, "rate"); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto rate_i = rate_vary.generator ? rate + rate_vary.offset * i : rate; - for (size_t j = 0; j < (size_t)n; ++j) { - auto rate_ij = rate_vary.draw ? rate_i[j] : rate_i[0]; - y_i[j] = monty::random::exponential_rate(state, rate_ij); - } - } - - return sexp_matrix(ret, n, n_streams); -} - - -template -cpp11::sexp monty_rng_exponential_mean(SEXP ptr, int n, cpp11::doubles r_mean, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "exponential_mean"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * mean = REAL(r_mean); - auto mean_vary = check_input_type(r_mean, n, n_streams, "mean"); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto mean_i = mean_vary.generator ? mean + mean_vary.offset * i : mean; - for (size_t j = 0; j < (size_t)n; ++j) { - auto mean_ij = mean_vary.draw ? mean_i[j] : mean_i[0]; - y_i[j] = monty::random::exponential_mean(state, mean_ij); - } - } - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_normal(SEXP ptr, int n, - cpp11::doubles r_mean, cpp11::doubles r_sd, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "normal"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * mean = REAL(r_mean); - const double * sd = REAL(r_sd); - auto mean_vary = check_input_type(r_mean, n, n_streams, "mean"); - auto sd_vary = check_input_type(r_sd, n, n_streams, "sd"); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto mean_i = mean_vary.generator ? mean + mean_vary.offset * i : mean; - auto sd_i = sd_vary.generator ? sd + sd_vary.offset * i : sd; - for (size_t j = 0; j < (size_t)n; ++j) { - auto mean_ij = mean_vary.draw ? mean_i[j] : mean_i[0]; - auto sd_ij = sd_vary.draw ? sd_i[j] : sd_i[0]; - y_i[j] = monty::random::normal(state, mean_ij, sd_ij); - } - } - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_binomial(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_prob, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "binomial"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * size = REAL(r_size); - const double * prob = REAL(r_prob); - auto size_vary = check_input_type(r_size, n, n_streams, "size"); - auto prob_vary = check_input_type(r_prob, n, n_streams, "prob"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto size_i = size_vary.generator ? size + size_vary.offset * i : size; - auto prob_i = prob_vary.generator ? prob + prob_vary.offset * i : prob; - for (size_t j = 0; j < (size_t)n; ++j) { - auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; - auto prob_ij = prob_vary.draw ? prob_i[j] : prob_i[0]; - y_i[j] = monty::random::binomial(state, size_ij, prob_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - - -template -cpp11::sexp monty_rng_beta_binomial_prob(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_prob, - cpp11::doubles r_rho, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "beta_binomial_prob"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * size = REAL(r_size); - const double * prob = REAL(r_prob); - const double * rho = REAL(r_rho); - auto size_vary = check_input_type(r_size, n, n_streams, "size"); - auto prob_vary = check_input_type(r_prob, n, n_streams, "prob"); - auto rho_vary = check_input_type(r_rho, n, n_streams, "rho"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto size_i = size_vary.generator ? size + size_vary.offset * i : size; - auto prob_i = prob_vary.generator ? prob + prob_vary.offset * i : prob; - auto rho_i = rho_vary.generator ? rho + rho_vary.offset * i : rho; - for (size_t j = 0; j < (size_t)n; ++j) { - auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; - auto prob_ij = prob_vary.draw ? prob_i[j] : prob_i[0]; - auto rho_ij = rho_vary.draw ? rho_i[j] : rho_i[0]; - y_i[j] = monty::random::beta_binomial_prob(state, size_ij, prob_ij, rho_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - - -template -cpp11::sexp monty_rng_beta_binomial_ab(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_a, - cpp11::doubles r_b, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "beta_binomial_ab"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * size = REAL(r_size); - const double * a = REAL(r_a); - const double * b = REAL(r_b); - auto size_vary = check_input_type(r_size, n, n_streams, "size"); - auto a_vary = check_input_type(r_a, n, n_streams, "a"); - auto b_vary = check_input_type(r_b, n, n_streams, "b"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto size_i = size_vary.generator ? size + size_vary.offset * i : size; - auto a_i = a_vary.generator ? a + a_vary.offset * i : a; - auto b_i = b_vary.generator ? b + b_vary.offset * i : b; - for (size_t j = 0; j < (size_t)n; ++j) { - auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; - auto a_ij = a_vary.draw ? a_i[j] : a_i[0]; - auto b_ij = b_vary.draw ? b_i[j] : b_i[0]; - y_i[j] = monty::random::beta_binomial_ab(state, size_ij, a_ij, b_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - - -template -cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_prob, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "negative_binomial_prob"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * size = REAL(r_size); - const double * prob = REAL(r_prob); - auto size_vary = check_input_type(r_size, n, n_streams, "size"); - auto prob_vary = check_input_type(r_prob, n, n_streams, "prob"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto size_i = size_vary.generator ? size + size_vary.offset * i : size; - auto prob_i = prob_vary.generator ? prob + prob_vary.offset * i : prob; - for (size_t j = 0; j < (size_t)n; ++j) { - auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; - auto prob_ij = prob_vary.draw ? prob_i[j] : prob_i[0]; - y_i[j] = monty::random::negative_binomial_prob(state, size_ij, prob_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - - -template -cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_mu, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "negative_binomial_mu"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * size = REAL(r_size); - const double * mu = REAL(r_mu); - auto size_vary = check_input_type(r_size, n, n_streams, "size"); - auto mu_vary = check_input_type(r_mu, n, n_streams, "mu"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto size_i = size_vary.generator ? size + size_vary.offset * i : size; - auto mu_i = mu_vary.generator ? mu + mu_vary.offset * i : mu; - for (size_t j = 0; j < (size_t)n; ++j) { - auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; - auto mu_ij = mu_vary.draw ? mu_i[j] : mu_i[0]; - y_i[j] = monty::random::negative_binomial_mu(state, size_ij, mu_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - - -template -cpp11::sexp monty_rng_poisson(SEXP ptr, int n, cpp11::doubles r_lambda, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "poisson"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * lambda = REAL(r_lambda); - auto lambda_vary = check_input_type(r_lambda, n, n_streams, "lambda"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto lambda_i = lambda_vary.generator ? lambda + lambda_vary.offset * i : - lambda; - for (size_t j = 0; j < (size_t)n; ++j) { - auto lambda_ij = lambda_vary.draw ? lambda_i[j] : lambda_i[0]; - y_i[j] = monty::random::poisson(state, lambda_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_multinomial(SEXP ptr, int n, - cpp11::doubles r_size, - cpp11::doubles r_prob, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "multinomial"); - const int n_streams = rng->size(); - - const double * size = REAL(r_size); - const double * prob = REAL(r_prob); - auto size_vary = check_input_type(r_size, n, n_streams, "size"); - auto prob_vary = check_input_type2(r_prob, n, n_streams, "prob"); - const int len = prob_vary.len; - - // Normally we return a block of doubles with the first 'n' entries - // being the results for the first generator, the second 'n' for the - // second, and so on. Here the first n * len are the first generator - // (with the first 'len' being the first sample. - cpp11::writable::doubles ret = - cpp11::writable::doubles(len * n * n_streams); - double * y = REAL(ret); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + len * n * i; - auto size_i = size_vary.generator ? size + size_vary.offset * i : size; - auto prob_i = prob_vary.generator ? prob + prob_vary.offset * i : prob; - for (size_t j = 0; j < (size_t)n; ++j) { - auto size_ij = size_vary.draw ? size_i[j] : size_i[0]; - auto prob_ij = prob_vary.draw ? prob_i + j * len : prob_i; - auto y_ij = y_i + j * len; - monty::random::multinomial(state, size_ij, prob_ij, len, - y_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - errors.report("generators", 4, true); - - if (n_streams == 1) { - ret.attr("dim") = cpp11::writable::integers{len, n}; - } else { - ret.attr("dim") = cpp11::writable::integers{len, n, n_streams}; - } - return ret; -} - -template -cpp11::sexp monty_rng_hypergeometric(SEXP ptr, int n, - cpp11::doubles r_n1, cpp11::doubles r_n2, - cpp11::doubles r_k, int n_threads) { - T *rng = safely_read_externalptr(ptr, "hypergeometric"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * n1 = REAL(r_n1); - const double * n2 = REAL(r_n2); - const double * k = REAL(r_k); - auto n1_vary = check_input_type(r_n1, n, n_streams, "n1"); - auto n2_vary = check_input_type(r_n2, n, n_streams, "n1"); - auto k_vary = check_input_type(r_k, n, n_streams, "k"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto n1_i = n1_vary.generator ? n1 + n1_vary.offset * i : n1; - auto n2_i = n2_vary.generator ? n2 + n2_vary.offset * i : n2; - auto k_i = k_vary.generator ? k + k_vary.offset * i : k; - for (size_t j = 0; j < (size_t)n; ++j) { - auto n1_ij = n1_vary.draw ? n1_i[j] : n1_i[0]; - auto n2_ij = n2_vary.draw ? n2_i[j] : n2_i[0]; - auto k_ij = k_vary.draw ? k_i[j] : k_i[0]; - y_i[j] = monty::random::hypergeometric(state, n1_ij, n2_ij, k_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - - -template -cpp11::sexp monty_rng_gamma_scale(SEXP ptr, int n, - cpp11::doubles r_shape, - cpp11::doubles r_scale, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "gamma_scale"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * shape = REAL(r_shape); - const double * scale = REAL(r_scale); - auto shape_vary = check_input_type(r_shape, n, n_streams, "shape"); - auto scale_vary = check_input_type(r_scale, n, n_streams, "scale"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto shape_i = shape_vary.generator ? shape + shape_vary.offset * i : shape; - auto scale_i = scale_vary.generator ? scale + scale_vary.offset * i : scale; - for (size_t j = 0; j < (size_t)n; ++j) { - auto shape_ij = shape_vary.draw ? shape_i[j] : shape_i[0]; - auto scale_ij = scale_vary.draw ? scale_i[j] : scale_i[0]; - y_i[j] = monty::random::gamma_scale(state, shape_ij, scale_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_gamma_rate(SEXP ptr, int n, - cpp11::doubles r_shape, - cpp11::doubles r_rate, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "gamma_rate"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * shape = REAL(r_shape); - const double * rate = REAL(r_rate); - auto shape_vary = check_input_type(r_shape, n, n_streams, "shape"); - auto rate_vary = check_input_type(r_rate, n, n_streams, "rate"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto shape_i = shape_vary.generator ? shape + shape_vary.offset * i : shape; - auto rate_i = rate_vary.generator ? rate + rate_vary.offset * i : rate; - for (size_t j = 0; j < (size_t)n; ++j) { - auto shape_ij = shape_vary.draw ? shape_i[j] : shape_i[0]; - auto rate_ij = rate_vary.draw ? rate_i[j] : rate_i[0]; - y_i[j] = monty::random::gamma_rate(state, shape_ij, rate_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - -// Below here is very repetitive, and could probably be deduplicated -// with some clever template magic. Most of the faff is because we -// want to support 4 modes of taking 1 or 2 parameters (each varying -// or not over draws and generators) -template -cpp11::sexp monty_rng_cauchy(SEXP ptr, int n, - cpp11::doubles r_location, - cpp11::doubles r_scale, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "cauchy"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * location = REAL(r_location); - const double * scale = REAL(r_scale); - auto location_vary = check_input_type(r_location, n, n_streams, "location"); - auto scale_vary = check_input_type(r_scale, n, n_streams, "scale"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto location_i = location_vary.generator ? location + location_vary.offset * i : location; - auto scale_i = scale_vary.generator ? scale + scale_vary.offset * i : scale; - for (size_t j = 0; j < (size_t)n; ++j) { - auto location_ij = location_vary.draw ? location_i[j] : location_i[0]; - auto scale_ij = scale_vary.draw ? scale_i[j] : scale_i[0]; - y_i[j] = monty::random::cauchy(state, location_ij, scale_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_beta(SEXP ptr, int n, - cpp11::doubles r_a, - cpp11::doubles r_b, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "beta"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * a = REAL(r_a); - const double * b = REAL(r_b); - auto a_vary = check_input_type(r_a, n, n_streams, "a"); - auto b_vary = check_input_type(r_b, n, n_streams, "b"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto a_i = a_vary.generator ? a + a_vary.offset * i : a; - auto b_i = b_vary.generator ? b + b_vary.offset * i : b; - for (size_t j = 0; j < (size_t)n; ++j) { - auto a_ij = a_vary.draw ? a_i[j] : a_i[0]; - auto b_ij = b_vary.draw ? b_i[j] : b_i[0]; - y_i[j] = monty::random::beta(state, a_ij, b_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_rng_truncated_normal(SEXP ptr, int n, - cpp11::doubles r_mean, - cpp11::doubles r_sd, - cpp11::doubles r_min, - cpp11::doubles r_max, - int n_threads) { - T *rng = safely_read_externalptr(ptr, "truncated_normal"); - const int n_streams = rng->size(); - cpp11::writable::doubles ret = cpp11::writable::doubles(n * n_streams); - double * y = REAL(ret); - - const double * mean = REAL(r_mean); - const double * sd = REAL(r_sd); - const double * min = REAL(r_min); - const double * max = REAL(r_max); - - auto mean_vary = check_input_type(r_mean, n, n_streams, "mean"); - auto sd_vary = check_input_type(r_sd, n, n_streams, "sd"); - auto min_vary = check_input_type(r_min, n, n_streams, "min"); - auto max_vary = check_input_type(r_max, n, n_streams, "max"); - - monty::utils::openmp_errors errors(n_streams); - -#ifdef _OPENMP -#pragma omp parallel for schedule(static) num_threads(n_threads) -#endif - for (int i = 0; i < n_streams; ++i) { - try { - auto &state = rng->state(i); - auto y_i = y + n * i; - auto mean_i = mean_vary.generator ? mean + mean_vary.offset * i : mean; - auto sd_i = sd_vary.generator ? sd + sd_vary.offset * i : sd; - auto min_i = min_vary.generator ? min + min_vary.offset * i : min; - auto max_i = max_vary.generator ? max + max_vary.offset * i : max; - for (size_t j = 0; j < (size_t)n; ++j) { - auto mean_ij = mean_vary.draw ? mean_i[j] : mean_i[0]; - auto sd_ij = sd_vary.draw ? sd_i[j] : sd_i[0]; - auto min_ij = min_vary.draw ? min_i[j] : min_i[0]; - auto max_ij = max_vary.draw ? max_i[j] : max_i[0]; - y_i[j] = monty::random::truncated_normal(state, mean_ij, sd_ij, min_ij, max_ij); - } - } catch (std::exception const& e) { - errors.capture(e, i); - } - } - - errors.report("generators", 4, true); - - return sexp_matrix(ret, n, n_streams); -} - -template -cpp11::sexp monty_legacy_rng_state(SEXP ptr) { - T *rng = safely_read_externalptr(ptr, "rng_state"); - auto state = rng->export_state(); - size_t len = sizeof(typename T::int_type) * state.size(); - cpp11::writable::raws ret(len); - std::memcpy(RAW(ret), state.data(), len); - return ret; -} - - -template -void monty_legacy_rng_set_state(SEXP ptr, cpp11::raws r_state) { - T *rng = safely_read_externalptr(ptr, "rng_set_state"); - - using int_type = typename T::int_type; - const auto len = rng->state_size() * sizeof(int_type); - if ((size_t)r_state.size() != len) { - cpp11::stop("'state' must be a raw vector of length %d (but was %d)", - len, r_state.size()); - } - std::vector state(len); - std::memcpy(state.data(), RAW(r_state), len); - rng->import_state(state); -} - - -[[cpp11::register]] -SEXP monty_rng_alloc(cpp11::sexp r_seed, int n_streams, bool deterministic) { - return monty_rng_alloc(r_seed, n_streams, deterministic); -} - -[[cpp11::register]] -void monty_legacy_rng_jump(SEXP ptr) { - monty_legacy_rng_jump(ptr); -} - -[[cpp11::register]] -void monty_legacy_rng_long_jump(SEXP ptr) { - monty_legacy_rng_long_jump(ptr); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_random_real(SEXP ptr, int n, int n_threads) { - return monty_rng_random_real(ptr, n, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_random_normal(SEXP ptr, int n, int n_threads, - std::string algorithm) { - cpp11::sexp ret; - if (algorithm == "box_muller") { - constexpr auto a = monty::random::algorithm::normal::box_muller; - ret = monty_rng_random_normal(ptr, n, n_threads); - } else if (algorithm == "polar") { - constexpr auto a = monty::random::algorithm::normal::polar; - ret = monty_rng_random_normal(ptr, n, n_threads); - } else if (algorithm == "ziggurat") { - constexpr auto a = monty::random::algorithm::normal::ziggurat; - ret = monty_rng_random_normal(ptr, n, n_threads); - } else { - cpp11::stop("Unknown normal algorithm '%s'", algorithm.c_str()); - } - return ret; -} - -[[cpp11::register]] -cpp11::sexp monty_rng_uniform(SEXP ptr, int n, - cpp11::doubles r_min, - cpp11::doubles r_max, - int n_threads) { - return monty_rng_uniform(ptr, n, r_min, r_max, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_exponential_rate(SEXP ptr, int n, cpp11::doubles r_rate, - int n_threads) { - return monty_rng_exponential_rate(ptr, n, r_rate, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_exponential_mean(SEXP ptr, int n, cpp11::doubles r_mean, - int n_threads) { - return monty_rng_exponential_mean(ptr, n, r_mean, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_normal(SEXP ptr, int n, cpp11::doubles r_mean, - cpp11::doubles r_sd, int n_threads, - std::string algorithm) { - cpp11::sexp ret; - if (algorithm == "box_muller") { - constexpr auto a = monty::random::algorithm::normal::box_muller; - ret = monty_rng_normal(ptr, n, r_mean, r_sd, n_threads); - } else if (algorithm == "polar") { - constexpr auto a = monty::random::algorithm::normal::polar; - ret = monty_rng_normal(ptr, n, r_mean, r_sd, n_threads); - } else if (algorithm == "ziggurat") { - constexpr auto a = monty::random::algorithm::normal::ziggurat; - ret = monty_rng_normal(ptr, n, r_mean, r_sd, n_threads); - } else { - cpp11::stop("Unknown normal algorithm '%s'", algorithm.c_str()); - } - return ret; -} - -[[cpp11::register]] -cpp11::sexp monty_rng_binomial(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_prob, - int n_threads) { - return monty_rng_binomial(ptr, n, r_size, r_prob, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_beta_binomial_ab(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_a, - cpp11::doubles r_b, - int n_threads) { - return monty_rng_beta_binomial_ab(ptr, n, r_size, r_a, r_b, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_beta_binomial_prob(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_prob, - cpp11::doubles r_rho, - int n_threads) { - return monty_rng_beta_binomial_prob(ptr, n, r_size, r_prob, r_rho, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_negative_binomial_prob(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_prob, - int n_threads) { - return monty_rng_negative_binomial_prob(ptr, n, r_size, r_prob, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_negative_binomial_mu(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_mu, - int n_threads) { - return monty_rng_negative_binomial_mu(ptr, n, r_size, r_mu, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_hypergeometric(SEXP ptr, int n, - cpp11::doubles r_n1, - cpp11::doubles r_n2, - cpp11::doubles r_k, - int n_threads) { - return monty_rng_hypergeometric(ptr, n, r_n1, r_n2, r_k, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_gamma_scale(SEXP ptr, int n, - cpp11::doubles r_shape, - cpp11::doubles r_scale, - int n_threads) { - return monty_rng_gamma_scale(ptr, n, r_shape, r_scale, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_gamma_rate(SEXP ptr, int n, - cpp11::doubles r_shape, - cpp11::doubles r_rate, - int n_threads) { - return monty_rng_gamma_rate(ptr, n, r_shape, r_rate, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_poisson(SEXP ptr, int n, - cpp11::doubles r_lambda, - int n_threads) { - return monty_rng_poisson(ptr, n, r_lambda, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_cauchy(SEXP ptr, int n, - cpp11::doubles r_location, - cpp11::doubles r_scale, - int n_threads) { - return monty_rng_cauchy(ptr, n, r_location, r_scale, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_beta(SEXP ptr, int n, - cpp11::doubles r_a, - cpp11::doubles r_b, - int n_threads) { - return monty_rng_beta(ptr, n, r_a, r_b, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_multinomial(SEXP ptr, int n, - cpp11::doubles r_size, cpp11::doubles r_prob, - int n_threads) { - return monty_rng_multinomial(ptr, n, r_size, r_prob, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_rng_truncated_normal(SEXP ptr, int n, - cpp11::doubles r_mean, - cpp11::doubles r_sd, - cpp11::doubles r_min, - cpp11::doubles r_max, - int n_threads) { - return monty_rng_truncated_normal(ptr, n, r_mean, r_sd, r_min, r_max, n_threads); -} - -[[cpp11::register]] -cpp11::sexp monty_legacy_rng_state(SEXP ptr) { - return monty_legacy_rng_state(ptr); -} - -[[cpp11::register]] -void monty_legacy_rng_set_state(SEXP ptr, cpp11::raws r_state) { - return monty_legacy_rng_set_state(ptr, r_state); -} diff --git a/src/rng_pointer.cpp b/src/rng_pointer.cpp deleted file mode 100644 index 74746754..00000000 --- a/src/rng_pointer.cpp +++ /dev/null @@ -1,84 +0,0 @@ -#include - -#ifdef _OPENMP -#include -#endif - -#include -#include - -[[cpp11::register]] -cpp11::sexp monty_rng_pointer_init(int n_streams, cpp11::sexp seed, - int long_jump, std::string algorithm) { - cpp11::sexp ret; - - using namespace monty::random; - if (algorithm == "xoshiro256starstar") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro256plusplus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro256plus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro128starstar") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro128plusplus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro128plus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoroshiro128starstar") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoroshiro128plusplus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoroshiro128plus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro512starstar") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro512plusplus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else if (algorithm == "xoshiro512plus") { - ret = r::rng_pointer_init(n_streams, seed, long_jump); - } else { - cpp11::stop("Unknown algorithm '%s'", algorithm.c_str()); - } - - return ret; -} - -[[cpp11::register]] -void monty_rng_pointer_sync(cpp11::environment obj, std::string algorithm) { - using namespace monty::random; - if (algorithm == "xoshiro256starstar") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro256plusplus") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro256plus") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro128starstar") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro128plusplus") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro128plus") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoroshiro128starstar") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoroshiro128plusplus") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoroshiro128plus") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro512starstar") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro512plusplus") { - r::rng_pointer_sync(obj); - } else if (algorithm == "xoshiro512plus") { - r::rng_pointer_sync(obj); - } -} - -// This exists to check some error paths in rng_pointer_get; it is not -// for use by users. -[[cpp11::register]] -double test_rng_pointer_get(cpp11::environment obj, int n_streams) { - using namespace monty::random; - auto rng = r::rng_pointer_get(obj, n_streams); - return random_real(rng->state(0)); -} diff --git a/src/test_rng.cpp b/src/test_rng.cpp index 5eb40e22..96b2d573 100644 --- a/src/test_rng.cpp +++ b/src/test_rng.cpp @@ -4,8 +4,10 @@ #include + #include -#include +#include + template std::string to_string(const T& t) { std::ostringstream ss; @@ -13,9 +15,10 @@ std::string to_string(const T& t) { return ss.str(); } -template -std::vector test_xoshiro_run1(cpp11::environment ptr) { - auto rng = monty::random::r::rng_pointer_get(ptr, 1); +[[cpp11::register]] +std::vector test_xoshiro_run(cpp11::sexp ptr) { + using default_rng64 = monty::random::prng; + auto rng = cpp11::as_cpp>(ptr).get(); auto& state = rng->state(0); constexpr int n = 10; @@ -33,36 +36,3 @@ std::vector test_xoshiro_run1(cpp11::environment ptr) { return ret; } - -[[cpp11::register]] -std::vector test_xoshiro_run(cpp11::environment obj) { - const auto algorithm = cpp11::as_cpp(obj["algorithm"]); - std::vector ret; - if (algorithm == "xoshiro256starstar") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro256plusplus") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro256plus") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro128starstar") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro128plusplus") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro128plus") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoroshiro128starstar") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoroshiro128plusplus") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoroshiro128plus") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro512starstar") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro512plusplus") { - ret = test_xoshiro_run1(obj); - } else if (algorithm == "xoshiro512plus") { - ret = test_xoshiro_run1(obj); - } - - return ret; -} diff --git a/tests/testthat/test-random.R b/tests/testthat/test-random.R index 2dbd2c9c..43a6dee8 100644 --- a/tests/testthat/test-random.R +++ b/tests/testthat/test-random.R @@ -1,12 +1,3 @@ -## This file will eventually replace test-rng.R, but until it does we -## use the legacy interface as the source of truth. -test_that("can generate a random number", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_real(s), cmp$random_real(1)) -}) - - test_that("can print rng state", { s <- monty_rng_create() res <- evaluate_promise(withVisible(print(s))) @@ -24,59 +15,28 @@ test_that("can get the length of rng state as the number of streams", { }) -test_that("can generate a random number from multiple streams at once", { - s <- monty_rng_create(seed = 42, n_streams = 10) - cmp <- monty_rng$new(seed = 42, n_streams = 10) - u <- monty_random_real(s) - expect_equal(u, drop(cmp$random_real(1))) - expect_null(dim(u)) -}) - - -test_that("can sample multiple random numbers from a single stream", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - u <- monty_random_n_real(10, s) - expect_null(dim(u)) - expect_equal(drop(u), cmp$random_real(10)) -}) - - -test_that("can sample multiple random numbers from multiple streams", { - s <- monty_rng_create(n_streams = 3, seed = 42) - cmp <- monty_rng$new(n_streams = 3, seed = 42) - u <- monty_random_n_real(10, s) - expect_equal(dim(u), c(10, 3)) - expect_equal(u, cmp$random_real(10)) -}) - - test_that("Preserve stream dimension on generation", { s <- monty_rng_create(n_streams = 1, seed = 42, preserve_stream_dimension = TRUE) - cmp <- monty_rng$new(n_streams = 1, seed = 42) + cmp <- monty_rng_create(n_streams = 1, seed = 42) u <- monty_random_n_real(10, s) expect_equal(dim(u), c(10, 1)) - expect_equal(u, cbind(cmp$random_real(10))) + expect_equal(u, cbind(monty_random_n_real(10, cmp))) u <- monty_random_n_real(1, s) expect_equal(dim(u), c(1, 1)) - expect_equal(u, cbind(cmp$random_real(1))) + expect_equal(u, cbind(monty_random_n_real(1, cmp))) }) test_that("can get and set rng state", { s <- monty_rng_create(seed = 42, n_streams = 10) - cmp <- monty_rng$new(seed = 42, n_streams = 10) r <- monty_rng_state(s) expect_type(r, "raw") expect_length(r, 320) - expect_equal(r, cmp$state()) monty_rng_set_state(rev(r), s) - cmp$set_state(rev(r)) expect_equal(monty_rng_state(s), rev(r)) - expect_equal(monty_random_real(s), drop(cmp$random_real(1))) expect_error( monty_rng_set_state(r[-1], s), @@ -132,16 +92,6 @@ test_that("Can vary parameters across streams", { }) -test_that("can sample from binomial distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_binomial(10, 0.3, s), - cmp$binomial(1, 10, 0.3)) - expect_equal(monty_random_n_binomial(10, 7, 0.3, s), - cmp$binomial(10, 7, 0.3)) -}) - - test_that("protect against pointer serialisation", { s <- monty_rng_create(seed = 42) s2 <- unserialize(serialize(s, NULL)) @@ -168,153 +118,3 @@ test_that("validate input size with multiple streams", { expect_error(monty_random_binomial(10, c(0.1, 0.2), s), "Expected 'prob' to have length 3 or 1, not 2") }) - - -test_that("can sample from exponential distribution (rate)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_exponential_rate(0.5, s), - cmp$exponential_rate(1, 0.5)) - expect_equal(monty_random_n_exponential_rate(10, 0.5, s), - cmp$exponential_rate(10, 0.5)) -}) - - -test_that("can sample from exponential distribution (mean)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_exponential_mean(0.5, s), - cmp$exponential_mean(1, 0.5)) - expect_equal(monty_random_n_exponential_mean(10, 0.5, s), - cmp$exponential_mean(10, 0.5)) -}) - - -test_that("can sample from poisson distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_poisson(0.5, s), - cmp$poisson(1, 0.5)) - expect_equal(monty_random_n_poisson(10, 0.5, s), - cmp$poisson(10, 0.5)) -}) - - -test_that("can sample from beta distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_beta(2, 3.1, s), - cmp$beta(1, 2, 3.1)) - expect_equal(monty_random_n_beta(10, 2, 3.1, s), - cmp$beta(10, 2, 3.1)) -}) - - -test_that("can sample from cauchy distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_cauchy(-1.5, 5, s), - cmp$cauchy(1, -1.5, 5)) - expect_equal(monty_random_n_cauchy(10, -1.5, 5, s), - cmp$cauchy(10, -1.5, 5)) -}) - - -test_that("can sample from gamma distribution (scale)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_gamma_scale(2, 5, s), - cmp$gamma_scale(1, 2, 5)) - expect_equal(monty_random_n_gamma_scale(10, 2, 5, s), - cmp$gamma_scale(10, 2, 5)) -}) - - -test_that("can sample from gamma distribution (rate)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_gamma_rate(2, 5, s), - cmp$gamma_rate(1, 2, 5)) - expect_equal(monty_random_n_gamma_rate(10, 2, 5, s), - cmp$gamma_rate(10, 2, 5)) -}) - - -test_that("can sample from negative binomial distribution (prob)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_negative_binomial_prob(200, 0.05, s), - cmp$negative_binomial_prob(1, 200, 0.05)) - expect_equal(monty_random_n_negative_binomial_prob(10, 200, 0.05, s), - cmp$negative_binomial_prob(10, 200, 0.05)) -}) - - -test_that("can sample from negative binomial distribution (mu)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_negative_binomial_mu(200, 65.1, s), - cmp$negative_binomial_mu(1, 200, 65.1)) - expect_equal(monty_random_n_negative_binomial_mu(10, 200, 65.1, s), - cmp$negative_binomial_mu(10, 200, 65.1)) -}) - - -test_that("can sample from normal distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_normal(-1.5, 5, s), - cmp$normal(1, -1.5, 5)) - expect_equal(monty_random_n_normal(10, -1.5, 5, s), - cmp$normal(10, -1.5, 5)) -}) - - -test_that("can sample from uniform distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_uniform(-1.5, 5, s), - cmp$uniform(1, -1.5, 5)) - expect_equal(monty_random_n_uniform(10, -1.5, 5, s), - cmp$uniform(10, -1.5, 5)) -}) - - -test_that("can sample from beta binomial distribution (prob)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_beta_binomial_prob(100, 0.3, 0.125, s), - cmp$beta_binomial_prob(1, 100, 0.3, 0.125)) - expect_equal(monty_random_n_beta_binomial_prob(10, 100, 0.3, 0.125, s), - cmp$beta_binomial_prob(10, 100, 0.3, 0.125)) -}) - - -test_that("can sample from beta binomial distribution (ab)", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_beta_binomial_ab(100, 1.5, 8.5, s), - cmp$beta_binomial_ab(1, 100, 1.5, 8.5)) - expect_equal(monty_random_n_beta_binomial_ab(10, 100, 1.5, 8.5, s), - cmp$beta_binomial_ab(10, 100, 1.5, 8.5)) -}) - - -test_that("can sample from hypergeometric distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_hypergeometric(7, 10, 8, s), - cmp$hypergeometric(1, 7, 10, 8)) - expect_equal(monty_random_n_hypergeometric(10, 7, 10, 8, s), - cmp$hypergeometric(10, 7, 10, 8)) -}) - - -test_that("can sample from truncated normal distribution", { - s <- monty_rng_create(seed = 42) - cmp <- monty_rng$new(seed = 42) - expect_equal(monty_random_truncated_normal(1, 4, -2, 3, s), - cmp$truncated_normal(1, 1, 4, -2, 3)) - expect_equal(monty_random_n_truncated_normal(10, 1, 4, -2, 3, s), - cmp$truncated_normal(10, 1, 4, -2, 3)) -}) diff --git a/tests/testthat/test-rng-distributed.R b/tests/testthat/test-rng-distributed.R deleted file mode 100644 index 16832bac..00000000 --- a/tests/testthat/test-rng-distributed.R +++ /dev/null @@ -1,35 +0,0 @@ -test_that("Can create a set of distributed random number states", { - s <- monty_rng_distributed_state(1L, n_nodes = 2) - expect_type(s, "list") - expect_length(s, 2) - cmp <- monty_rng_create(seed = 1) - expect_equal(s[[1]], monty_rng_state(cmp)) - monty_rng_long_jump(cmp) - expect_equal(s[[2]], monty_rng_state(cmp)) -}) - - -test_that("Can create a set of distributed random number states", { - s <- monty_rng_distributed_state(1L, n_streams = 3, n_nodes = 2) - expect_type(s, "list") - expect_length(s, 2) - cmp <- monty_rng_create(seed = 1, n_streams = 3) - expect_equal(s[[1]], monty_rng_state(cmp)) - monty_rng_long_jump(cmp) - expect_equal(s[[2]], monty_rng_state(cmp)) -}) - - -test_that("can create distributed rng pointers", { - algorithm <- "xoroshiro128starstar" - p <- monty_rng_distributed_pointer(1L, n_streams = 3, n_nodes = 2, - algorithm = algorithm) - expect_type(p, "list") - expect_length(p, 2) - expect_equal( - p[[1]]$state(), - monty_rng_pointer$new(1, 3, 0, algorithm = algorithm)$state()) - expect_equal( - p[[2]]$state(), - monty_rng_pointer$new(1, 3, 1, algorithm = algorithm)$state()) -}) diff --git a/tests/testthat/test-rng-interface.R b/tests/testthat/test-rng-interface.R deleted file mode 100644 index 444890b8..00000000 --- a/tests/testthat/test-rng-interface.R +++ /dev/null @@ -1,99 +0,0 @@ -test_that("Can create pointer object", { - obj <- monty_rng_pointer$new() - expect_true(obj$is_current()) - expect_type(obj$state(), "raw") - expect_length(obj$state(), 32) - expect_equal(obj$algorithm, "xoshiro256plus") - r <- obj$state() - obj$sync() - expect_identical(r, obj$state()) -}) - - -test_that("Using object requires sync", { - obj <- monty_rng_pointer$new() - expect_true(obj$is_current()) - r <- obj$state() - test_xoshiro_run(obj) - expect_false(obj$is_current()) - obj$sync() - expect_true(obj$is_current()) - expect_false(identical(obj$state(), r)) -}) - - -test_that("Fetching state syncs", { - obj <- monty_rng_pointer$new() - expect_true(obj$is_current()) - r <- obj$state() - test_xoshiro_run(obj) - expect_false(obj$is_current()) - expect_false(identical(obj$state(), r)) - expect_true(obj$is_current()) -}) - - -test_that("Invalidated pointers can be rebuilt", { - obj1 <- monty_rng_pointer$new() - obj2 <- corrupt_pointer(obj1) - expect_equal( - test_xoshiro_run(obj2), - test_xoshiro_run(obj1)) - - ## This saves that the pointer is invalid: - obj3 <- corrupt_pointer(obj2) - expect_error( - test_xoshiro_run(obj3), - "Can't unserialise an rng pointer that was not synced") - - ## But if we sync things it's ok: - obj2$sync() - obj4 <- corrupt_pointer(obj2) - expect_equal( - test_xoshiro_run(obj4), - test_xoshiro_run(obj1)) -}) - - -test_that("can't create invalid pointer types", { - expect_error( - monty_rng_pointer$new(algorithm = "mt19937"), - "Unknown algorithm 'mt19937'") -}) - - -test_that("Validate pointers on fetch", { - obj <- monty_rng_pointer$new(algorithm = "xoshiro256starstar") - expect_error( - test_rng_pointer_get(obj, 1), - "Incorrect rng type: given xoshiro256starstar, expected xoshiro256plus") - obj <- monty_rng_pointer$new(algorithm = "xoshiro256plus", n_streams = 4) - expect_error( - test_rng_pointer_get(obj, 20), - "Requested a rng with 20 streams but only have 4") - expect_silent( - test_rng_pointer_get(obj, 0)) - expect_silent( - test_rng_pointer_get(obj, 1)) -}) - - -test_that("Create pointer with a long jump", { - s0 <- monty_rng_pointer$new(1, 4, 0)$state() - s1 <- monty_rng_pointer$new(1, 4, 1)$state() - s2 <- monty_rng_pointer$new(1, 4, 2)$state() - - cmp <- monty_rng_create(seed = 1, n_streams = 4) - expect_equal(s0, monty_rng_state(cmp)) - expect_equal(s1, monty_rng_long_jump(s0)) - expect_equal(s2, monty_rng_long_jump(s0, 2)) -}) - - -test_that("can summarise errors", { - r <- monty_rng$new(n_streams = 10) - err <- expect_error( - r$binomial(1, 1, -1), - "10 generators reported errors") - expect_match(err$message, "and 6 more") -}) diff --git a/tests/testthat/test-rng.R b/tests/testthat/test-rng.R index d2365cb7..15fee6f6 100644 --- a/tests/testthat/test-rng.R +++ b/tests/testthat/test-rng.R @@ -333,6 +333,18 @@ test_that("rexp agrees with stats::rexp", { }) + +test_that("can sample from exponential distribution (mean)", { + r1 <- monty_rng_create(seed = 42) + r2 <- monty_rng_create(seed = 42) + expect_equal(monty_random_exponential_mean(0.5, r1), + monty_random_exponential_rate(2, r2)) + expect_equal(monty_random_n_exponential_mean(10, 0.5, r1), + monty_random_n_exponential_rate(10, 2, r2)) +}) + + + test_that("continue stream", { rng1 <- monty_rng_create(seed = 1) rng2 <- monty_rng_create(seed = 1) @@ -1070,6 +1082,9 @@ test_that("negative_binomial_mu follows from negative_binomial_prob", { prob <- 0.3 size <- 20 mu <- size * (1 - prob) / prob + expect_identical( + monty_random_negative_binomial_mu(size, mu, rng1), + monty_random_negative_binomial_prob(size, prob, rng2)) expect_identical( monty_random_n_negative_binomial_mu(100, size, mu, rng1), monty_random_n_negative_binomial_prob(100, size, prob, rng2)) @@ -1125,6 +1140,9 @@ test_that("beta_binomial_prob follows from beta_binomial_ab", { b <- 5 prob <- a / (a + b) rho <- 1 / (a + b + 1) + expect_identical( + monty_random_beta_binomial_prob(size, prob, rho, rng1), + monty_random_beta_binomial_ab(size, a, b, rng2)) expect_identical( monty_random_n_beta_binomial_prob(100, size, prob, rho, rng1), monty_random_n_beta_binomial_ab(100, size, a, b, rng2)) diff --git a/tests/testthat/test-xoshiro.R b/tests/testthat/test-xoshiro.R index 059b4f23..a9cb0175 100644 --- a/tests/testthat/test-xoshiro.R +++ b/tests/testthat/test-xoshiro.R @@ -1,14 +1,9 @@ test_that("xoshiro output agrees with reference data", { - path <- "xoshiro-ref" - status <- list() - for (name in dir(path)) { - obj <- monty_rng_pointer$new(seed = 42, algorithm = name) - res <- test_xoshiro_run(obj) - obj$sync() - len <- if (grepl("^xoshiro128", name)) 4 else 8 - s <- matrix(obj$state(), len)[rev(seq_len(len)), ] - s_str <- apply(s, 2, paste, collapse = "") - cmp <- readLines(sprintf("%s/%s", path, name)) - expect_equal(c(res, s_str), cmp) - } + obj <- monty_rng_create(seed = 42) + res <- test_xoshiro_run(obj) + len <- 8 + s <- matrix(monty_rng_state(obj), len)[rev(seq_len(len)), ] + s_str <- apply(s, 2, paste, collapse = "") + cmp <- readLines("xoshiro-ref/xoshiro256plus") + expect_equal(c(res, s_str), cmp) })