Run an ensemble simulation of a stock-and-flow model, varying initial conditions and/or parameters in the range specified in range. The ensemble can be run in parallel using multiple threads by first setting use_threads(). The results are returned as a data.frame with summary statistics and optionally individual simulations.
Usage
ensemble(
sfm,
n = 10,
return_sims = FALSE,
range = NULL,
cross = TRUE,
quantiles = c(0.025, 0.975),
only_stocks = TRUE,
keep_nonnegative_flow = TRUE,
keep_nonnegative_stock = FALSE,
keep_unit = TRUE,
verbose = TRUE
)Arguments
- sfm
Stock-and-flow model, object of class
sdbuildR_xmile.- n
Number of simulations to run in the ensemble. When range 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 10.
- return_sims
If TRUE, return the individual simulations in the ensemble. Set to FALSE to save memory. Defaults to FALSE.
- range
A named list specifying parameter ranges for ensemble conditions. Names must correspond to existing stock or constant variable names 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 range 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).
- only_stocks
If TRUE, only return stocks in output, discarding flows and auxiliaries. If FALSE, flows and auxiliaries are saved, which slows down the simulation. Defaults to FALSE.
- keep_nonnegative_flow
If TRUE, keeps original non-negativity setting of flows. Defaults to TRUE.
- keep_nonnegative_stock
If TRUE, keeps original non-negativity setting of stocks Defaults to FALSE.
- keep_unit
If TRUE, keeps units of variables. Defaults to TRUE.
- verbose
If TRUE, print details and duration of simulation. Defaults to TRUE.
Value
Object of class sdbuildR_ensemble, 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 return_sims is TRUE. The iteration number is indicated by column "i". If range was specified, the condition is indicated by column "j".
- summary
data.frame with summary statistics of the ensemble, including quantiles specified in quantiles. If range was specified, summary statistics are calculated for each condition (j) in the ensemble.
- n
Number of simulations run in the ensemble (per condition j if range is specified).
- n_total
Total number of simulations run in the ensemble (across all conditions if range is specified).
- n_conditions
Total number of conditions.
- conditions
data.frame with the conditions used in the ensemble, if range is specified.
- init
List with df (if return_sims = TRUE) and summary, containing data.frame with the initial values of the stocks used in the ensemble.
- constants
List with df (if return_sims = TRUE) and summary, containing data.frame with the constant parameters used in the ensemble.
- script
Julia script used for the ensemble simulation.
- duration
Duration of the simulation in seconds.
- ...
Other parameters passed to ensemble
Details
To run large simulations, it is recommended to limit the output size by saving fewer values. To create a reproducible ensemble simulation, set a seed using sim_specs().
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.
Examples
# Load example and set simulation language to Julia
sfm <- xmile("predator_prey") |> sim_specs(language = "Julia")
# Set random initial conditions
sfm <- build(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 every 1 time units and discard
# the first 100 time units, use:
sfm <- sim_specs(sfm, save_at = 1, save_from = 100)
# Run ensemble simulation with 100 simulations
sims <- ensemble(sfm, n = 100)
#> Running a total of 100 simulations
#> Starting Julia ...
#> Connecting to Julia TCP server at localhost:11980 ...
#> Setting up Julia environment for sdbuildR...
#> Simulation took 9.0726 seconds
plot(sims)
# Plot individual trajectories
sims <- ensemble(sfm, n = 10, return_sims = TRUE)
#> Running a total of 10 simulations
#> Simulation took 0.5612 seconds
plot(sims, type = "sims")
# Specify which trajectories to plot
plot(sims, type = "sims", i = 1)
# Plot the median with lighter individual trajectories
plot(sims, central_tendency = "median", type = "sims", alpha = 0.1)
# 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 hunderd samples per
# condition. By default, the parameters are crossed, meaning that all
# combinations of the parameters are run.
sims <- ensemble(sfm,
n = 50,
range = list("predator" = c(10, 50), "delta" = c(.025, .05))
)
#> Running a total of 200 simulations for 4 conditions (50 simulations per condition)
#> Simulation took 2.3831 seconds
plot(sims)
# By default, a maximum of nine conditions is plotted.
# Plot specific conditions:
plot(sims, j = c(1, 3), nrows = 1)
# Generate a non-crossed design, where the length of each range needs to be
# equal:
sims <- ensemble(sfm,
n = 10, cross = FALSE,
range = list(
"predator" = c(10, 20, 30),
"delta" = c(.020, .025, .03)
)
)
#> Running a total of 30 simulations for 3 conditions (10 simulations per condition)
#> Simulation took 0.5252 seconds
plot(sims, nrows = 3)
# Run simulation in parallel
use_threads(4)
#> Warning: n is set to 4, which is higher than the number of available cores minus 1. Setting it to 3.
sims <- ensemble(sfm, n = 10)
#> Running a total of 10 simulations
#> Simulation took 1.3652 seconds
# Stop using threads
use_threads(stop = TRUE)
# Close Julia
use_julia(stop = TRUE)
#> Julia session closed.