Vectorise ebola model

[pratikunterwegs]

This PR is WIP on a project to release {epidemics} v0.3.0 focusing on functionality related to the Ebola model.

This PR:

1. Fixes #169 by allowing the Ebola model to accept vectors of infection parameters following Tidyverse recycling rules;
2. Fixes #170 by allowing the Ebola model to accept lists of rate intervention sets (time-dependence is the only other composable element and is not allowed to be vectorised);
3. Fixes #172 by using {withr} to reset the random seed across parameter-set and scenario combinations, while pre-generating and using local seeds for each replicate within each parameter-scenario combination; this ensures that each $i$-th replicate of each scenario uses the same seed;
4. Fixes #173 by lifting the dynamic parameter calculation and compartmental transitions code into an internal, 'unsafe' function without input checks called .model_ebola_internal(); the user facing model_ebola() holds code for input checking, parameter-scenario combination, and model output handling;
5. Fixes #174 by introducing the replicates argument to model_ebola() with a default of 100 replicates;
6. Fixes #163 by showing how to use {withr} for seed management in the Ebola vignette;
7. Fixes #164 by showing multiple replicates for each model run in the Ebola vignette.

Note that the continuous benchmarking check is expected to fail (or show reduced performance) as the default for model_ebola() is now 100 replicates.

A small example showing seed preservation:

library(epidemics)
library(ggplot2)

# Prepare population and parameters
demography_vector <- 67000 # small population
contact_matrix <- matrix(1)

# manual case counts divided by pop size rather than proportions as small sizes
# introduce errors when converting to counts in the model code; extra
# individuals may appear
infectious <- 1
exposed <- 10
initial_conditions <- matrix(
  c(demography_vector - infectious - exposed, exposed, infectious, 0, 0, 0) /
    demography_vector,
  nrow = 1
)
rownames(contact_matrix) <- "full_pop"
pop <- population(
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  initial_conditions = initial_conditions
)

compartments <- c(
  "susceptible", "exposed", "infectious", "hospitalised", "funeral", "removed"
)

# prepare integer values for expectations on data length
time_end <- 100L
replicates <- 10L

npi_list <- list(
  scenario_baseline = NULL,
  scenario_00 = list(
    transmission_rate = intervention(
      "transmission", "rate", 50, 100, 0.1
    )
  ),
  scenario_01 = list(
    transmission_rate = intervention(
      "transmission", "rate", 50, 100, 1.0
    )
  )
)
output <- withr::with_seed(
  1,
  model_ebola(
    population = pop,
    intervention = npi_list,
    replicates = 5L, # replicates
    time_end = time_end
  )
)
data <- output[, unlist(data, recursive = FALSE), by = "scenario"]

ggplot(data[compartment == "exposed"], aes(time, value)) +
  geom_line(
    aes(
      col = as.factor(scenario),
      group = interaction(scenario, replicate)
    )
  ) +
  facet_grid(~replicate)

^{Created on 2024-04-08 with reprex v2.0.2}

[github-actions[bot]]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if e99365afa2cc4fe1a21952917ae8fa32df97b6b1 is merged into main:

• :heavy_check_mark:default_ode: 15ms -> 15.1ms [-0.51%, +1.92%]
• :heavy_check_mark:default_ode_interventions: 73.4ms -> 73ms [-1.05%, +0.1%]
• :heavy_check_mark:default_ode_param_vec: 838ms -> 841ms [-0.1%, +0.89%]
• :heavy_check_mark:default_ode_paramvec_intervs: 6.25s -> 6.27s [-0.22%, +0.72%]
• :exclamation::snail:ebola: 43ms -> 615ms [+1325.2%, +1337.77%] Further explanation regarding interpretation and methodology can be found in the documentation.

[github-actions[bot]]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if e99365afa2cc4fe1a21952917ae8fa32df97b6b1 is merged into main:

• :heavy_check_mark:default_ode: 14.8ms -> 14.9ms [-0.45%, +1.88%]
• :rocket:default_ode_interventions: 73.8ms -> 73ms [-1.94%, -0.28%]
• :heavy_check_mark:default_ode_param_vec: 843ms -> 845ms [-1.37%, +1.72%]
• :rocket:default_ode_paramvec_intervs: 6.29s -> 6.24s [-1.23%, -0.31%]
• :exclamation::snail:ebola: 37.2ms -> 607ms [+1522.66%, +1535.58%] Further explanation regarding interpretation and methodology can be found in the documentation.

[github-actions[bot]]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 6258c0c61d6d0b976b19bd860947731789403956 is merged into main:

• :heavy_check_mark:default_ode: 15.6ms -> 15.5ms [-2.1%, +0.56%]
• :rocket:default_ode_interventions: 74.2ms -> 73.8ms [-1.11%, -0.09%]
• :heavy_check_mark:default_ode_param_vec: 851ms -> 849ms [-1.05%, +0.55%]
• :heavy_check_mark:default_ode_paramvec_intervs: 6.34s -> 6.33s [-0.69%, +0.36%]
• :exclamation::snail:ebola: 46.6ms -> 628ms [+1238.71%, +1257.64%] Further explanation regarding interpretation and methodology can be found in the documentation.

[pratikunterwegs]

Thanks @joshwlambert for looking through this, you've already managed to catch a few issues that I'll fix (func docs, notes, captions). Will go through the comments and answer there.

[github-actions[bot]]

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 6258c0c61d6d0b976b19bd860947731789403956 is merged into main:

• :heavy_check_mark:default_ode: 15.2ms -> 15.1ms [-2.41%, +1.07%]
• :heavy_check_mark:default_ode_interventions: 73.7ms -> 73.2ms [-1.6%, +0.23%]
• :heavy_check_mark:default_ode_param_vec: 837ms -> 836ms [-0.88%, +0.76%]
• :heavy_check_mark:default_ode_paramvec_intervs: 6.38s -> 6.34s [-1.38%, +0.1%]
• :exclamation::snail:ebola: 47.8ms -> 617ms [+1183.9%, +1199.66%] Further explanation regarding interpretation and methodology can be found in the documentation.

[pratikunterwegs]

Thanks @joshwlambert for reviewing. I'm merging this PR now so the next one can be up for review; happy to take more feedback on the Ebola model as issues.