ensemble: Run ensemble simulations

Description

Run large-scale (i.e., ensemble) simulations of stock-and-flow models, varying initial conditions and/or constants specified in conditions.

Usage

ensemble(
  object,
  n = 10,
  conditions = NULL,
  cross = TRUE,
  quantiles = c(0.025, 0.975),
  verbose = TRUE,
  ...
)

Value

Object of class ensemble_stockflow, which is a list containing:

success: If TRUE, simulation was successful. If FALSE, simulation failed.
error_message: If success is FALSE, contains the error message.
df: data.frame with simulation results in long format, if save_sims is TRUE. The iteration number is indicated by column "sim". If conditions was specified, the condition is indicated by column "condition".
summary: data.frame with summary statistics of the ensemble, including quantiles specified in quantiles. If conditions was specified, summary statistics are calculated for each condition in the ensemble.
n: Number of simulations run in the ensemble (per condition if conditions is specified).
n_total: Total number of simulations run in the ensemble (across all conditions if conditions is specified).
n_conditions: Total number of conditions.
conditions: data.frame with the conditions used in the ensemble, if conditions was specified.
init: List with df (if save_sims = TRUE) and summary, containing data.frame with the initial values of the stocks used in the ensemble.
constants: List with df (if save_sims = TRUE) and summary, containing data.frame with the constant parameters used in the ensemble.
script: Script used for the ensemble simulation.
duration: Duration of the simulation in seconds.
...: Other parameters passed to ensemble

Arguments

object

Stock-and-flow model, object of class stockflow.

n

Number of simulations to run in the ensemble. When conditions is specified, n defines the number of simulations to run per condition. If each condition only needs to be run once, set n = 1. Defaults to n = 10.

conditions

A named list specifying fixed values for ensemble conditions. Names must correspond to stocks or constants in the model. Each list element should be a numeric vector of values to test.

If cross = TRUE (default), all combinations of values are generated. For example, list(param1 = c(1, 2), param2 = c(10, 20)) creates 4 conditions: (1,10), (1,20), (2,10), (2,20).

If cross = FALSE, values are paired element-wise, requiring all vectors to have equal length. For example, list(param1 = c(1, 2, 3), param2 = c(10, 20, 30)) creates 3 conditions: (1,10), (2,20), (3,30). Defaults to NULL (no parameter variation).

cross

If TRUE, cross the parameters in the conditions list to generate all possible combinations of parameters. Defaults to TRUE.

quantiles

Quantiles to calculate in the summary, e.g., c(0.025, 0.975).

verbose

If TRUE (default), print details and duration of simulation.

...

Optional arguments passed to sim_settings(); these can be used to override the simulation specifications set in the model object.

Details

It is strongly recommended to reduce the size of the simulation output by saving fewer values with save_at or save_n in sim_settings().

By default, only summary statistics across simulations are returned. To return individual simulations, set save_sims = TRUE in sim_settings() or pass save_sims = TRUE via ... in ensemble(). Note that returning individual simulations can consume a lot of memory for large ensembles.

For simulations in Julia, the ensemble can be run in parallel using multiple threads by setting nthreads in use_julia(). For simulations in R, use future::plan() to control parallel execution.

To create a reproducible ensemble simulation, set a seed using sim_settings().

If you do not see any variation within a condition of the ensemble (i.e., the confidence bands are virtually non-existent), there are likely no random elements in your model. Without these, there can be no variability in the model. Try specifying a random initial condition or adding randomness to other model elements (see examples).

Examples

Run this code

# Ensemble simulation in R (no parallelization)
# Load example
sfm <- stockflow("predator_prey")

# Set random initial conditions
sfm <- update(sfm, c(predator, prey),
  eqn = runif(1, min = 20, max = 80)
)

# For ensemble simulations, it is highly recommended to reduce the
# returned output. For example, to save only 20 values per simulation:
sfm <- sim_settings(sfm, save_n = 20)

# Run ensemble simulation with a small number of simulations
sims <- ensemble(sfm, n = 3)
if (interactive()) plot(sims)

if (FALSE) { # Sys.getenv("NOT_CRAN") == "true"
# To plot individual trajectories, rerun the ensemble with save_sims = TRUE.
# Note that this can consume a lot of memory for large simulations.
sims <- ensemble(sfm, n = 10, save_sims = TRUE)
plot(sims, which = "sims")

# Specify which trajectories to plot
plot(sims, which = "sims", sim = 1)

# Plot the median with lighter individual trajectories
plot(sims, central_tendency = "median", which = "sims", alpha = 0.1)

# For larger ensembles, we can use parallelization with future
if (requireNamespace("future", quietly = TRUE) &&
  requireNamespace("future.apply", quietly = TRUE)) {
  future::plan(future::multisession, workers = 4)
}

# Ensembles can also be run with exact values for the initial conditions
# and parameters. Below, we vary the initial values of the predator and the
# birth rate of the predators (delta). We generate a hundred samples per
# condition. By default, the parameters are crossed, meaning that all
# combinations of the parameters are run.
sims <- ensemble(sfm,
  n = 50,
  conditions = list(predator = c(10, 50), delta = c(.025, .05))
)

plot(sims)

# By default, a maximum of nine conditions is plotted.
# Plot specific conditions:
plot(sims, condition = c(1, 3), nrows = 1)

# Generate a non-crossed design, where the length of each conditions vector
# needs to be equal:
sims <- ensemble(sfm,
  n = 10, cross = FALSE,
  conditions = list(
    predator = c(10, 20, 30),
    delta = c(.020, .025, .03)
  )
)
plot(sims, nrows = 3)

# Stop parallelization after use
if (requireNamespace("future", quietly = TRUE) &&
  requireNamespace("future.apply", quietly = TRUE)) {
  future::plan(future::sequential)
}
}