Merge branch 'main' into implement-rerun

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: monty
 Title: Monte Carlo Models
-Version: 0.2.19
+Version: 0.2.22
 Authors@R: c(person("Rich", "FitzJohn", role = c("aut", "cre"),
                     email = "[email protected]"),
              person("Wes", "Hinsley", role = "aut"),
@@ -31,7 +31,6 @@ Suggests:
     coda,
     decor,
     knitr,
-    dust,
     mockery,
     mvtnorm,
     numDeriv,
@@ -43,5 +42,3 @@ Suggests:
 Config/testthat/edition: 3
 Language: en-GB
 VignetteBuilder: knitr
-Remotes:
-    mrc-ide/dust

R/cpp11.R

@@ -40,8 +40,12 @@ monty_rng_binomial <- function(ptr, n, r_size, r_prob, n_threads, is_float) {
   .Call(`_monty_monty_rng_binomial`, ptr, n, r_size, r_prob, n_threads, is_float)
 }
 
-monty_rng_nbinomial <- function(ptr, n, r_size, r_prob, n_threads, is_float) {
-  .Call(`_monty_monty_rng_nbinomial`, ptr, n, r_size, r_prob, n_threads, is_float)
+monty_rng_negative_binomial_prob <- function(ptr, n, r_size, r_prob, n_threads, is_float) {
+  .Call(`_monty_monty_rng_negative_binomial_prob`, ptr, n, r_size, r_prob, n_threads, is_float)
+}
+
+monty_rng_negative_binomial_mu <- function(ptr, n, r_size, r_mu, n_threads, is_float) {
+  .Call(`_monty_monty_rng_negative_binomial_mu`, ptr, n, r_size, r_mu, n_threads, is_float)
 }
 
 monty_rng_hypergeometric <- function(ptr, n, r_n1, r_n2, r_k, n_threads, is_float) {
@@ -76,6 +80,18 @@ monty_rng_state <- function(ptr, is_float) {
   .Call(`_monty_monty_rng_state`, ptr, is_float)
 }
 
+cpp_monty_random_real <- function(ptr) {
+  .Call(`_monty_cpp_monty_random_real`, ptr)
+}
+
+cpp_monty_random_binomial <- function(r_size, r_prob, ptr) {
+  .Call(`_monty_cpp_monty_random_binomial`, r_size, r_prob, ptr)
+}
+
+cpp_monty_random_exponential_rate <- function(r_rate, ptr) {
+  .Call(`_monty_cpp_monty_random_exponential_rate`, r_rate, ptr)
+}
+
 monty_rng_pointer_init <- function(n_streams, seed, long_jump, algorithm) {
   .Call(`_monty_monty_rng_pointer_init`, n_streams, seed, long_jump, algorithm)
 }

R/domain.R

@@ -54,7 +54,7 @@ monty_domain_expand <- function(domain, packer) {
       arg = "domain")
   }
 
-  nms_full <- packer$parameters
+  nms_full <- packer$names()
   nms_map <- packer$unpack(nms_full)
   nms_logical <- names(nms_map)
 
@@ -76,5 +76,5 @@ monty_domain_expand <- function(domain, packer) {
                     domain[!i, , drop = FALSE])
   }
 
-  domain[order(match(rownames(domain), packer$parameters)), , drop = FALSE]
+  domain[order(match(rownames(domain), packer$names())), , drop = FALSE]
 }

R/dsl-distributions.R

@@ -103,7 +103,7 @@ distr_gamma_rate <- distribution(
 
 distr_gamma_scale <- distribution(
   name = "Gamma",
-  variant = "rate",
+  variant = "scale",
   density = function(x, shape, scale) {
     dgamma(x, shape, scale = scale, log = TRUE)
   },
@@ -125,6 +125,36 @@ distr_hypergeometric <- distribution(
     mean = quote(k * n1 / (n1 + n2))),
   cpp = list(density = "hypergeometric", sample = "hypergeometric"))
 
+distr_negative_binomial_prob <- distribution(
+  name = "NegativeBinomial",
+  variant = "prob",
+  density = function(x, size, prob) {
+    dnbinom(x, size, prob = prob, log = TRUE)
+  },
+  domain = c(0, Inf),
+  sample = function(rng, size, prob) rng$negative_binomial_prob(1, size, prob),
+  expr = list(
+    density = quote(lgamma(x + size) - lgamma(size) - lgamma(x + 1) +
+                      x * log(1 - prob) + size * log(prob)),
+    mean = quote(size * (1 - prob) / prob)),
+  cpp = list(density = "negative_binomial_prob",
+             sample = "negative_binomial_prob"))
+
+distr_negative_binomial_mu <- distribution(
+  name = "NegativeBinomial",
+  variant = "mu",
+  density = function(x, size, mu) {
+    dnbinom(x, size, mu = mu, log = TRUE)
+  },
+  domain = c(0, Inf),
+  sample = function(rng, size, mu) rng$negative_binomial_mu(1, size, mu),
+  expr = list(
+    density = quote(lgamma(x + size) - lgamma(size) - lgamma(x + 1) +
+                      size * log(size) + x * log(mu) -
+                      (size + x) * log(size + mu)),
+    mean = quote(mu)),
+  cpp = list(density = "negative_binomial_mu", sample = "negative_binomial_mu"))
+
 distr_normal <- distribution(
   name = "Normal",
   density = function(x, mean, sd) dnorm(x, mean, sd, log = TRUE),
@@ -166,6 +196,8 @@ dsl_distributions <- local({
     distr_gamma_rate,
     distr_gamma_scale,
     distr_hypergeometric,
+    distr_negative_binomial_prob,
+    distr_negative_binomial_mu,
     distr_normal,
     distr_poisson,
     distr_uniform)

R/dsl-generate.R

@@ -56,7 +56,7 @@ dsl_generate_direct_sample <- function(dat, env, meta) {
   exprs <- lapply(dat$exprs, dsl_generate_sample_expr, meta)
   body <- c(call("<-", meta[["pars"]], quote(list())),
             exprs,
-            bquote(unlist(.(meta[["pars"]])[packer$parameters], FALSE, FALSE)))
+            bquote(unlist(.(meta[["pars"]])[packer$names()], FALSE, FALSE)))
   as_function(alist(rng = ), body, env)
 }
 

R/model-function.R

@@ -63,7 +63,7 @@ monty_model_function <- function(density, packer = NULL, fixed = NULL) {
   }
 
   monty_model(
-    list(parameters = packer$parameters,
+    list(parameters = packer$names(),
          density = function(x) {
            rlang::inject(density(!!!packer$unpack(x)))
          }))

R/packer.R

@@ -1,9 +1,11 @@
-##' Build a parameter packer, which can be used in models to translate
-##' between an unstructured vector of numbers (the vector being
-##' updated by an MCMC for example) to a structured list of named
-##' values, which is easier to program against.  We refer to the
-##' process of taking a named list of scalars, vectors and arrays and
-##' converting into a single vector "packing" and the inverse
+##' Build a packer, which can be used to translate between an
+##' unstructured vector of numbers (the vector being updated by an
+##' MCMC for example) and a structured list of named values, which is
+##' easier to program against.  This is useful for the bridge between
+##' model parameters and a models's implementation, but it is also
+##' useful for the state vector in a state-space model.  We refer to
+##' the process of taking a named list of scalars, vectors and arrays
+##' and converting into a single vector "packing" and the inverse
 ##' "unpacking".
 ##'
 ##' There are several places where it is most convenient to work in an
@@ -145,57 +147,64 @@
 ##' vector vs packing an array where all dimensions are 1.  See the
 ##' examples, and please let us know if the behaviour needs changing.
 ##'
-##' @title Build a parameter packer
+##' @title Build a packer
 ##'
-##' @param scalar Names of scalar parameters.  This is similar for
-##'   listing elements in `array` with values of 1, though elements in
+##' @param scalar Names of scalars.  This is similar for listing
+##'   elements in `array` with values of 1, though elements in
 ##'   `scalar` will be placed ahead of those listed in `array` within
 ##'   the final parameter vector, and elements in `array` will have
 ##'   generated names that include square brackets.
 ##'
-##' @param array A list, where names correspond to the names of array
-##'   parameters and values correspond to the lengths of parameters.
-##'   Multiple dimensions are allowed (so if you provide an element
-##'   with two entries these represent dimensions of a matrix).
-##'   Zero-length integer vectors or `NULL` values are counted as
-##'   scalars, which allows you to put scalars at positions other than
-##'   the front of the packing vector. In future, you may be able to
-##'   use *strings* as values for the lengths, in which case these
-##'   will be looked for within `fixed`.
-##'
-##' @param fixed A named list of fixed parameters; these will be added
-##'   into the final list directly.  These typically represent
-##'   additional pieces of data that your model needs to run, but
-##'   which you are not performing inference on.
+##' @param array A list, where names correspond to the names of arrays
+##'   and values correspond to their lengths.  Multiple dimensions are
+##'   allowed (so if you provide an element with two entries these
+##'   represent dimensions of a matrix).  Zero-length integer vectors
+##'   or `NULL` values are counted as scalars, which allows you to put
+##'   scalars at positions other than the front of the packing
+##'   vector. In future, you may be able to use *strings* as values
+##'   for the lengths, in which case these will be looked for within
+##'   `fixed`.
+##'
+##' @param fixed A named list of fixed data to be inserted into the
+##'   final unpacked list; these will be added into the final list
+##'   directly.  In the parameter packer context, these typically
+##'   represent additional pieces of data that your model needs to
+##'   run, but which you are not performing inference on.
 ##'
 ##' @param process An arbitrary R function that will be passed the
-##'   final assembled parameter list; it may create any *additional*
-##'   entries, which will be concatenated onto the original list.  If
-##'   you use this you should take care not to return any values with
-##'   the same names as entries listed in `scalar`, `array` or
-##'   `fixed`, as this is an error (this is so that `pack()` is
-##'   not broken).  We will likely play around with this process in
-##'   future in order to get automatic differentiation to work.
-##'
-##' @return An object of class `monty_packer`, which has three
-##'   elements:
-##'
-##' * `parameters`: a character vector of computed parameter names;
-##'   these are the names that your statistical model will use.
+##'   final assembled list; it may create any *additional* entries,
+##'   which will be concatenated onto the original list.  If you use
+##'   this you should take care not to return any values with the same
+##'   names as entries listed in `scalar`, `array` or `fixed`, as this
+##'   is an error (this is so that `pack()` is not broken).  We will
+##'   likely play around with this process in future in order to get
+##'   automatic differentiation to work.
+##'
+##' @return An object of class `monty_packer`, which has elements:
+##'
+##' * `names`: a function that returns a character vector of computed
+##'   names; in the parameter packer context these are the names that
+##'   your statistical model will use.
+##'
 ##' * `unpack`: a function that can unpack an unstructured vector
 ##'   (say, from your statistical model parameters) into a structured
 ##'   list (say, for your generative model)
-##' * `pack`: a function that can pack your structured list of
-##'   parameters back into a numeric vector suitable for the
+##'
+##' * `pack`: a function that can pack your structured list of data
+##'   back into a numeric vector, for example suitable for a
 ##'   statistical model.  This ignores values created by a
-##'   `preprocess` function.
+##'   `preprocess` function and present in `fixed`.
+##'
 ##' * `index`: a function which produces a named list where each
-##'   element has the name of a value in `parameters` and each value
-##'   has the indices within an unstructured vector where these values
-##'   can be found.
-##' * `subset`: an experimental interface which can be used to subset a
-##'   packer to a packer for a subset of contents. Documentation will
-##'   be provided once the interface settles.
+##'   element has the name of a value in `scalar` or `array` and each
+##'   value has the indices within an unstructured vector where these
+##'   values can be found, in the shape of the data that would be
+##'   unpacked.  This is of limited most use to most people.
+##'
+##' * `subset`: an experimental interface which can be used to subset
+##'   a packer to a packer for a subset of contents. Documentation
+##'   will be provided once the interface settles, but this is for
+##'   advanced use only!
 ##'
 ##' @export
 ##'
@@ -260,7 +269,7 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL,
                          process = NULL) {
   call <- environment()
 
-  parameters <- character(0)
+  nms <- character(0)
   idx <- list()
   shape <- list()
 
@@ -279,15 +288,15 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL,
     }
     shape[scalar] <- rep(list(integer()), length(scalar))
     idx[scalar] <- as.list(seq_along(scalar))
-    parameters <- c(parameters, scalar)
+    nms <- c(nms, scalar)
   }
 
   len <- length(scalar) # start arrays after scalars
   if (!is.null(array)) {
     assert_named(array, unique = TRUE, call = call)
     for (nm in names(array)) {
       tmp <- prepare_pack_array(nm, array[[nm]], call)
-      parameters <- c(parameters, tmp$names)
+      nms <- c(nms, tmp$names)
       shape[[nm]] <- tmp$shape
       idx[[nm]] <- seq_len(tmp$n) + len
       len <- len + tmp$n
@@ -325,9 +334,9 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL,
 
   unpack <- function(x) {
     if (is.null(dim(x))) {
-      unpack_vector(x, parameters, len, idx, shape, fixed, process)
+      unpack_vector(x, nms, len, idx, shape, fixed, process)
     } else {
-      unpack_array(x, parameters, len, idx, shape, fixed, process)
+      unpack_array(x, nms, len, idx, shape, fixed, process)
     }
   }
 
@@ -379,7 +388,7 @@ monty_packer <- function(scalar = NULL, array = NULL, fixed = NULL,
     list(index = index, packer = packer)
   }
 
-  ret <- list(parameters = parameters,
+  ret <- list(names = function() nms,
               unpack = unpack,
               pack = pack,
               index = function() idx,
@@ -438,7 +447,7 @@ array_indices <- function(shape) {
 print.monty_packer <- function(x, ...) {
   cli::cli_h1("<monty_packer>")
   cli::cli_alert_info(
-    "Packing {length(x$parameters)} parameter{?s}: {squote(x$parameters)}")
+    "Packing {length(x$names())} parameter{?s}: {squote(x$names())}")
   cli::cli_alert_info(
     "Use '$pack()' to convert from a list to a vector")
   cli::cli_alert_info(
@@ -448,10 +457,10 @@ print.monty_packer <- function(x, ...) {
 }
 
 
-unpack_vector <- function(x, parameters, len, idx, shape, fixed, process) {
+unpack_vector <- function(x, nms, len, idx, shape, fixed, process) {
   call <- parent.frame()
   if (!is.null(names(x))) {
-    if (!identical(names(x), parameters)) {
+    if (!identical(names(x), nms)) {
       ## Here, we could do better I think with this message; we
       ## might pass thropuigh empty names, and produce some summary
       ## of different names.  Something for later though.
@@ -476,7 +485,7 @@ unpack_vector <- function(x, parameters, len, idx, shape, fixed, process) {
     err <- intersect(names(extra), names(res))
     if (length(err) > 0) {
       cli::cli_abort(
-        c("'process()' is trying to overwrite entries in parameters",
+        c("'process()' is trying to overwrite entries in your list",
           i = paste("The 'process()' function should only create elements",
                     "that are not already present in 'scalar', 'array'",
                     "or 'fixed', as this lets us reverse the transformation",
@@ -494,10 +503,10 @@ unpack_vector <- function(x, parameters, len, idx, shape, fixed, process) {
 }
 
 
-unpack_array <- function(x, parameters, len, idx, shape, fixed, process) {
+unpack_array <- function(x, nms, len, idx, shape, fixed, process) {
   call <- parent.frame()
   dn <- dimnames(x)
-  if (!is.null(dn) && !is.null(dn[[1]]) && !identical(dn[[1]], parameters)) {
+  if (!is.null(dn) && !is.null(dn[[1]]) && !identical(dn[[1]], nms)) {
     ## See comment above about reporting on this better
     cli::cli_abort("Incorrect rownames in input")
   }

R/rng.R

@@ -96,7 +96,10 @@
 ##' rng$binomial(5, 10, 0.3)
 ##'
 ##' # Negative binomially distributed random numbers with size and prob
-##' rng$nbinomial(5, 10, 0.3)
+##' 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)
@@ -283,9 +286,25 @@ monty_rng <- R6::R6Class(
     ##'   (between 0 and 1, length 1 or n)
     ##'
     ##' @param n_threads Number of threads to use; see Details
-    nbinomial = function(n, size, prob, n_threads = 1L) {
-      monty_rng_nbinomial(private$ptr, n, size, prob, n_threads,
-                          private$float)
+    negative_binomial_prob = function(n, size, prob, n_threads = 1L) {
+      monty_rng_negative_binomial_prob(private$ptr, n, size, prob, n_threads,
+                                       private$float)
+    },
+
+    ##' @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,
+                                     private$float)
     },
 
     ##' @description Generate `n` numbers from a hypergeometric distribution