diff --git a/DESCRIPTION b/DESCRIPTION index 2569d935..ceaac2c3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: monty Title: Monte Carlo Models -Version: 0.3.17 +Version: 0.3.18 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"), email = "rich.fitzjohn@gmail.com"), person("Wes", "Hinsley", role = "aut"), diff --git a/R/cpp11.R b/R/cpp11.R index 4e1d6618..818b87bd 100644 --- a/R/cpp11.R +++ b/R/cpp11.R @@ -184,6 +184,14 @@ cpp_monty_random_n_truncated_normal <- function(n_samples, mean, sd, min, max, p .Call(`_monty_cpp_monty_random_n_truncated_normal`, n_samples, mean, sd, min, max, ptr) } +cpp_monty_random_multinomial <- function(r_size, r_prob, ptr) { + .Call(`_monty_cpp_monty_random_multinomial`, r_size, r_prob, ptr) +} + +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) } diff --git a/R/random.R b/R/random.R index 2e4f6888..6d0e1069 100644 --- a/R/random.R +++ b/R/random.R @@ -661,3 +661,12 @@ monty_random_log_normal <- function(meanlog, sdlog, state) { monty_random_n_log_normal <- function(n_samples, meanlog, sdlog, state) { cpp_monty_random_n_log_normal(n_samples, meanlog, sdlog, state) } + + +monty_random_multinomial <- function(size, prob, state) { + cpp_monty_random_multinomial(size, prob, state) +} + +monty_random_n_multinomial <- function(n_samples, size, prob, state) { + cpp_monty_random_n_multinomial(n_samples, size, prob, state) +} diff --git a/inst/include/monty/random/multinomial.hpp b/inst/include/monty/random/multinomial.hpp index 77e60171..f64b8f80 100644 --- a/inst/include/monty/random/multinomial.hpp +++ b/inst/include/monty/random/multinomial.hpp @@ -53,8 +53,8 @@ namespace random { template __host__ __device__ -void multinomial(rng_state_type& rng_state, int size, const T& prob, - int prob_len, U& ret) { +void multinomial(rng_state_type& rng_state, real_type size, const T& prob, + int prob_len, U ret) { real_type p_tot = 0; for (int i = 0; i < prob_len; ++i) { if (prob[i] < 0) { diff --git a/src/cpp11.cpp b/src/cpp11.cpp index 638d9075..e1f8814a 100644 --- a/src/cpp11.cpp +++ b/src/cpp11.cpp @@ -330,6 +330,20 @@ extern "C" SEXP _monty_cpp_monty_random_n_truncated_normal(SEXP n_samples, SEXP return cpp11::as_sexp(cpp_monty_random_n_truncated_normal(cpp11::as_cpp>(n_samples), cpp11::as_cpp>(mean), cpp11::as_cpp>(sd), cpp11::as_cpp>(min), cpp11::as_cpp>(max), cpp11::as_cpp>(ptr))); END_CPP11 } +// random.cpp +cpp11::doubles cpp_monty_random_multinomial(cpp11::doubles r_size, cpp11::doubles r_prob, cpp11::sexp ptr); +extern "C" SEXP _monty_cpp_monty_random_multinomial(SEXP r_size, SEXP r_prob, SEXP ptr) { + BEGIN_CPP11 + return cpp11::as_sexp(cpp_monty_random_multinomial(cpp11::as_cpp>(r_size), cpp11::as_cpp>(r_prob), cpp11::as_cpp>(ptr))); + END_CPP11 +} +// random.cpp +cpp11::doubles cpp_monty_random_n_multinomial(size_t n_samples, cpp11::doubles r_size, cpp11::doubles r_prob, cpp11::sexp ptr); +extern "C" SEXP _monty_cpp_monty_random_n_multinomial(SEXP n_samples, SEXP r_size, SEXP r_prob, SEXP ptr) { + BEGIN_CPP11 + 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) { @@ -544,6 +558,7 @@ static const R_CallMethodDef CallEntries[] = { {"_monty_cpp_monty_random_gamma_scale", (DL_FUNC) &_monty_cpp_monty_random_gamma_scale, 3}, {"_monty_cpp_monty_random_hypergeometric", (DL_FUNC) &_monty_cpp_monty_random_hypergeometric, 4}, {"_monty_cpp_monty_random_log_normal", (DL_FUNC) &_monty_cpp_monty_random_log_normal, 3}, + {"_monty_cpp_monty_random_multinomial", (DL_FUNC) &_monty_cpp_monty_random_multinomial, 3}, {"_monty_cpp_monty_random_n_beta", (DL_FUNC) &_monty_cpp_monty_random_n_beta, 4}, {"_monty_cpp_monty_random_n_beta_binomial_ab", (DL_FUNC) &_monty_cpp_monty_random_n_beta_binomial_ab, 5}, {"_monty_cpp_monty_random_n_beta_binomial_prob", (DL_FUNC) &_monty_cpp_monty_random_n_beta_binomial_prob, 5}, @@ -555,6 +570,7 @@ static const R_CallMethodDef CallEntries[] = { {"_monty_cpp_monty_random_n_gamma_scale", (DL_FUNC) &_monty_cpp_monty_random_n_gamma_scale, 4}, {"_monty_cpp_monty_random_n_hypergeometric", (DL_FUNC) &_monty_cpp_monty_random_n_hypergeometric, 5}, {"_monty_cpp_monty_random_n_log_normal", (DL_FUNC) &_monty_cpp_monty_random_n_log_normal, 4}, + {"_monty_cpp_monty_random_n_multinomial", (DL_FUNC) &_monty_cpp_monty_random_n_multinomial, 4}, {"_monty_cpp_monty_random_n_negative_binomial_mu", (DL_FUNC) &_monty_cpp_monty_random_n_negative_binomial_mu, 4}, {"_monty_cpp_monty_random_n_negative_binomial_prob", (DL_FUNC) &_monty_cpp_monty_random_n_negative_binomial_prob, 4}, {"_monty_cpp_monty_random_n_normal_box_muller", (DL_FUNC) &_monty_cpp_monty_random_n_normal_box_muller, 4}, diff --git a/src/random.cpp b/src/random.cpp index bfa3723d..c7a81e3d 100644 --- a/src/random.cpp +++ b/src/random.cpp @@ -41,6 +41,46 @@ struct input { } }; +struct input_array { + const double * const data; + size_t len; + bool fixed; + + input_array(cpp11::doubles r_data, size_t expected, const char *name) : + data(REAL(r_data)) { + const auto dim = r_data.attr("dim"); + if (dim == R_NilValue) { + fixed = true; + size_ = r_data.size(); + } else { + if (Rf_length(dim) != 2) { + cpp11::stop("Expected '%s' to be a matrix", name); + } + size_ = INTEGER(dim)[0]; + const auto ncol = INTEGER(dim)[1]; + fixed = ncol == 1; + if (!fixed && static_cast(ncol) != expected) { + if (expected == 1) { + cpp11::stop("Expected '%s' to have 1 column, not %d", name, ncol); + } else { + cpp11::stop("Expected '%s' to have %d or 1 columns, not %d", + name, static_cast(expected), ncol); + } + } + } + } + + size_t size() const { + return size_; + } + + auto operator[](size_t i) const { + return data + (fixed ? 0 : i * size_); + } +private: + size_t size_; +}; + bool preserve_stream_dimension(size_t n, cpp11::sexp ptr) { return n > 1 || INTEGER(ptr.attr("preserve_stream_dimension"))[0] == 1; } @@ -158,7 +198,14 @@ void set_dimensions(size_t n_samples, size_t n_streams, cpp11::sexp ptr, cpp11:: if (preserve_stream_dimension(n_streams, ptr)) { ret.attr("dim") = cpp11::writable::integers{static_cast(n_samples), static_cast(n_streams)}; } +} +void set_dimensions(size_t n_state, size_t n_samples, size_t n_streams, cpp11::sexp ptr, cpp11::sexp ret) { + if (preserve_stream_dimension(n_streams, ptr)) { + ret.attr("dim") = cpp11::writable::integers{static_cast(n_state), static_cast(n_samples), static_cast(n_streams)}; + } else { + ret.attr("dim") = cpp11::writable::integers{static_cast(n_state), static_cast(n_samples)}; + } } template @@ -219,8 +266,8 @@ cpp11::doubles monty_random_sample_n_2(Fn fn, size_t n_samples, double * y = REAL(r_y); for (size_t i = 0; i < n_streams; ++i) { + auto& state_i = rng->state(i); for (size_t j = 0; j < n_samples; ++j) { - auto& state_i = rng->state(i); auto y_i = y + n_samples * i; y_i[j] = fn(state_i, a[i], b[i]); } @@ -756,3 +803,54 @@ cpp11::doubles cpp_monty_random_n_truncated_normal(size_t n_samples, // Other // multinomial +[[cpp11::register]] +cpp11::doubles cpp_monty_random_multinomial(cpp11::doubles r_size, + cpp11::doubles r_prob, + cpp11::sexp ptr) { + auto * rng = safely_read_externalptr(ptr, "multinomial"); + const size_t n_streams = rng->size(); + auto size = input(r_size, n_streams, "size"); + auto prob = input_array(r_prob, n_streams, "prob"); + const auto len = prob.size(); + + cpp11::writable::doubles r_y = cpp11::writable::doubles(n_streams * len); + double * y = REAL(r_y); + + for (size_t i = 0; i < n_streams; ++i) { + auto &state = rng->state(i); + monty::random::multinomial(state, size[i], prob[i], len, y + i * len); + } + + set_dimensions(len, n_streams, ptr, r_y); + + return r_y; +} + + +[[cpp11::register]] +cpp11::doubles cpp_monty_random_n_multinomial(size_t n_samples, + cpp11::doubles r_size, + cpp11::doubles r_prob, + cpp11::sexp ptr) { + auto * rng = safely_read_externalptr(ptr, "multinomial"); + const size_t n_streams = rng->size(); + auto size = input(r_size, n_streams, "size"); + auto prob = input_array(r_prob, n_streams, "prob"); + const auto len = prob.size(); + + cpp11::writable::doubles r_y = + cpp11::writable::doubles(n_samples * n_streams * len); + double * y = REAL(r_y); + + for (size_t i = 0; i < n_streams; ++i) { + auto &state = rng->state(i); + for (size_t j = 0; j < n_samples; ++j) { + auto y_ij = y + len * (n_samples * i + j); + monty::random::multinomial(state, size[i], prob[i], len, y_ij); + } + } + + set_dimensions(len, n_samples, n_streams, ptr, r_y); + + return r_y; +} diff --git a/tests/testthat/test-rng.R b/tests/testthat/test-rng.R index 3f034b9b..d2365cb7 100644 --- a/tests/testthat/test-rng.R +++ b/tests/testthat/test-rng.R @@ -502,23 +502,27 @@ test_that("multinomial algorithm is correct", { size <- 20 n <- 5 - res <- monty_rng$new(1, seed = 1L)$multinomial(n, size, prob) - ## Separate implementation of the core algorithm: cmp_multinomial <- function(rng, size, prob) { p <- prob / (1 - cumsum(c(0, prob[-length(prob)]))) ret <- numeric(length(prob)) for (i in seq_len(length(prob) - 1L)) { - ret[i] <- rng$binomial(1, size, p[i]) + ret[i] <- monty_random_binomial(size, p[i], rng) size <- size - ret[i] } ret[length(ret)] <- size ret } - rng2 <- monty_rng$new(1, seed = 1L) - cmp <- replicate(n, cmp_multinomial(rng2, size, prob)) - expect_equal(res, cmp) + r1 <- monty_rng_create(seed = 1) + r2 <- monty_rng_create(seed = 1) + res1 <- monty_random_multinomial(size, prob, r1) + res2 <- cmp_multinomial(r2, size, prob) + expect_equal(res1, res2) + + res1 <- monty_random_n_multinomial(n, size, prob, r1) + res2 <- replicate(n, cmp_multinomial(r2, size, prob)) + expect_equal(res1, res2) }) @@ -526,7 +530,8 @@ test_that("multinomial expectation is correct", { p <- runif(10) p <- p / sum(p) n <- 10000 - res <- monty_rng$new(1, seed = 1L)$multinomial(n, 100, p) + rng <- monty_rng_create(seed = 1) + res <- monty_random_n_multinomial(n, 100, p, rng) expect_equal(dim(res), c(10, n)) expect_equal(colSums(res), rep(100, n)) expect_equal(rowMeans(res), p * 100, tolerance = 1e-2) @@ -539,7 +544,8 @@ test_that("multinomial allows zero probs", { p <- p / sum(p) n <- 500 size <- 100 - res <- monty_rng$new(1, seed = 1L)$multinomial(n, size, p) + rng <- monty_rng_create(seed = 1) + res <- monty_random_n_multinomial(n, size, p, rng) expect_equal(res[4, ], rep(0, n)) expect_equal(colSums(res), rep(size, n)) @@ -549,104 +555,89 @@ test_that("multinomial allows zero probs", { test_that("multinomial allows non-normalised prob", { p <- runif(10, 0, 10) n <- 50 - res1 <- monty_rng$new(1, seed = 1L)$multinomial(n, 100, p) - res2 <- monty_rng$new(1, seed = 1L)$multinomial(n, 100, p / sum(p)) + res1 <- monty_random_n_multinomial(n, 100, p, monty_rng_create(seed = 1)) + res2 <- monty_random_n_multinomial(n, 100, p / sum(p), + monty_rng_create(seed = 1)) expect_equal(res1, res2) }) test_that("Invalid prob throws an error", { - r <- monty_rng$new(1, seed = 1L) + r <- monty_rng_create(seed = 1) expect_error( - r$multinomial(1, 10, c(0, 0, 0)), + monty_random_multinomial(10, c(0, 0, 0), r), "No positive prob in call to multinomial") expect_error( - r$multinomial(1, 10, c(-0.1, 0.6, 0.5)), + monty_random_multinomial(10, c(-0.1, 0.6, 0.5), r), "Negative prob passed to multinomial") }) -test_that("Can vary parameters for multinomial, single generator", { - np <- 7L - ng <- 1L - size <- 13 - n <- 17L - prob <- matrix(runif(np * n), np, n) - prob <- prob / rep(colSums(prob), each = np) - - rng <- monty_rng$new(1, seed = 1L) - cmp <- vapply(seq_len(n), function(i) rng$multinomial(1, size, prob[, i]), - numeric(np)) - res <- monty_rng$new(1, seed = 1L)$multinomial(n, size, prob) - expect_equal(res, cmp) - - expect_error( - monty_rng$new(1, seed = 1L)$multinomial(n, size, prob[, -5]), - "If 'prob' is a matrix, it must have 17 columns") - expect_error( - monty_rng$new(1, seed = 1L)$multinomial(n, size, prob[0, ]), - "Input parameters imply length of 'prob' of only 0 (< 2)", - fixed = TRUE) - expect_error( - monty_rng$new(1, seed = 1L)$multinomial(n, size, prob[1, , drop = FALSE]), - "Input parameters imply length of 'prob' of only 1 (< 2)", - fixed = TRUE) -}) - - -test_that("Can vary parameters by generator for multinomial", { +test_that("Can vary scalar parameters by generator for multinomial", { np <- 7L ng <- 3L - size <- 13 n <- 17L - prob <- array(runif(np * ng), c(np, 1, ng)) - prob <- prob / rep(colSums(prob), each = np) + # prob <- matrix(runif(np * ng), np, ng) + prob <- runif(np) + size <- 10 * seq_len(ng) + + r <- monty_rng_create(seed = 1, n_streams = ng) + state <- matrix(monty_rng_state(r), ncol = ng) - state <- matrix(monty_rng$new(ng, seed = 1L)$state(), ncol = ng) cmp <- vapply(seq_len(ng), function(i) { - monty_rng$new(1, seed = state[, i])$multinomial(n, size, prob[, , i]) + ri <- monty_rng_create(seed = state[, i]) + monty_random_n_multinomial(n, size[i], prob, ri) }, matrix(numeric(), np, n)) - res <- monty_rng$new(ng, seed = 1L)$multinomial(n, size, prob) + res <- monty_random_n_multinomial(n, size, prob, r) expect_equal(res, cmp) }) -test_that("Can vary parameters for multinomial, multiple generators", { +test_that("Can vary vector parameters by generator for multinomial", { np <- 7L ng <- 3L size <- 13 n <- 17L - prob <- array(runif(np * n * ng), c(np, n, ng)) - prob <- prob / rep(colSums(prob), each = np) - ## Setting up the expectation here is not easy, we need a set of - ## generators. This test exploits the fact that we alredy worked out - ## we could vary a parameter over draws with a single generator. - state <- matrix(monty_rng$new(ng, seed = 1L)$state(), ncol = ng) + prob <- matrix(runif(np * ng), np, ng) + + r <- monty_rng_create(seed = 1, n_streams = ng) + state <- matrix(monty_rng_state(r), ncol = ng) + cmp <- vapply(seq_len(ng), function(i) { - monty_rng$new(1, seed = state[, i])$multinomial(n, size, prob[, , i]) + ri <- monty_rng_create(seed = state[, i]) + monty_random_n_multinomial(n, size, prob[, i], ri) }, matrix(numeric(), np, n)) - res <- monty_rng$new(ng, seed = 1L)$multinomial(n, size, prob) + res <- monty_random_n_multinomial(n, size, prob, r) expect_equal(res, cmp) +}) + +test_that("can validate multinomial inputs", { + r <- monty_rng_create(seed = 1, n_streams = 3) expect_error( - monty_rng$new(ng, seed = 1L)$multinomial(n, size, prob[, -5, ]), - "If 'prob' is a 3d array, it must have 1 or 17 columns") - expect_error( - monty_rng$new(ng, seed = 1L)$multinomial(n, size, prob[, , -1]), - "If 'prob' is a 3d array, it must have 3 layers") + monty_random_multinomial(5, matrix(1, 7, 2), r), + "Expected 'prob' to have 3 or 1 columns, not 2") expect_error( - monty_rng$new(ng, seed = 1L)$multinomial(n, size, prob[0, , ]), - "Input parameters imply length of 'prob' of only 0 (< 2)", - fixed = TRUE) - ## Final bad inputs: - p4 <- array(prob, c(dim(prob), 1)) + monty_random_multinomial(5, array(1, c(7, 2, 3)), r), + "Expected 'prob' to be a matrix") + + r <- monty_rng_create(seed = 1, n_streams = 1) expect_error( - monty_rng$new(ng, seed = 1L)$multinomial(n, size, p4), - "'prob' must be a vector, matrix or 3d array") + monty_random_multinomial(5, matrix(1, 7, 2), r), + "Expected 'prob' to have 1 column, not 2") +}) + + +test_that("deterministic multinomial returns expectation", { + r <- monty_rng_create(seed = 1, deterministic = TRUE) + p <- runif(10) + p <- p / sum(p) + res <- monty_random_multinomial(100, p, r) + expect_equal(res, 100 * p) }) @@ -1366,11 +1357,11 @@ test_that("can draw weibull random numbers", { shape <- 5 scale <- 3 n <- 10000000 - + ans1 <- monty_random_n_weibull(n, shape, scale, monty_rng_create(seed = 1)) ans2 <- monty_random_n_weibull(n, shape, scale, monty_rng_create(seed = 1)) expect_identical(ans1, ans2) - + expect_equal(mean(ans1), scale * gamma(1 + 1 / shape), tolerance = 1e-3) true_var <- scale^2 * (gamma(1 + 2 / shape) - gamma(1 + 1 / shape)^2) expect_equal(var(ans1), true_var, tolerance = 1e-3) @@ -1381,10 +1372,10 @@ test_that("deterministic weibull returns mean", { n_reps <- 10 shape <- as.numeric(sample(10, n_reps, replace = TRUE)) scale <- as.numeric(sample(10, n_reps, replace = TRUE)) - + rng <- monty_rng_create(seed = 1, deterministic = TRUE) state <- monty_rng_state(rng) - + expect_equal( mapply(monty_random_weibull, shape, scale, MoreArgs = list(rng)), scale * gamma(1 + 1 / shape)) @@ -1396,7 +1387,7 @@ test_that("weibull random numbers prevent bad inputs", { r <- monty_rng_create(seed = 1) expect_equal(monty_random_weibull(0, 0, r), 0) expect_equal(monty_random_weibull(Inf, Inf, r), Inf) - + expect_error( monty_random_weibull(-1.1, 5.1, r), "Invalid call to Weibull with shape = -1.1, scale = 5.1") @@ -1410,13 +1401,13 @@ test_that("can draw log-normal random numbers", { meanlog <- 1.5 sdlog <- 0.5 n <- 10000000 - - ans1 <- + + ans1 <- monty_random_n_log_normal(n, meanlog, sdlog, monty_rng_create(seed = 1)) - ans2 <- + ans2 <- monty_random_n_log_normal(n, meanlog, sdlog, monty_rng_create(seed = 1)) expect_identical(ans1, ans2) - + expect_equal(mean(ans1), exp(meanlog + sdlog^2 / 2), tolerance = 1e-3) true_var <- (exp(sdlog^2) - 1) * exp(2 * meanlog + sdlog^2) expect_equal(var(ans1), true_var, tolerance = 1e-2) @@ -1427,10 +1418,10 @@ test_that("deterministic log-normal returns mean", { n_reps <- 10 meanlog <- as.numeric(sample(seq(-10, 10), n_reps, replace = TRUE)) sdlog <- as.numeric(sample(10, n_reps, replace = TRUE)) - + rng <- monty_rng_create(seed = 1, deterministic = TRUE) state <- monty_rng_state(rng) - + expect_equal( mapply(monty_random_log_normal, meanlog, sdlog, MoreArgs = list(rng)), exp(meanlog + sdlog^2 / 2)) @@ -1440,7 +1431,7 @@ test_that("deterministic log-normal returns mean", { test_that("log-normal random numbers prevent bad inputs", { r <- monty_rng_create(seed = 1) - + expect_error( monty_random_log_normal(1.1, -5.1, r), "Invalid call to log_normal with meanlog = 1.1, sdlog = -5.1")