Ability to save and load state and simulation history

[pwinskill]

Adding the ability to save (and re-load) a simulation state and history may be the most low-hanging fruit to improve model run times. For example, for a country-site file, once calibrated, we could run and save the historical part of the simulation run. Then anyone using that site file would just need to load it and stitch on the future part of any runs. This may cut effective run time in half for lots of the runs people are doing.

[pwinskill]

I'm just wondering @giovannic - should this be a feature request for individual? It seems like something that could possibly be generically implemented and I'm sure would be useful outside of just this model

[giovannic]

Yes, it would be nice for individual to provide a universal way to (de-)serialize Variables and Events. Although for malariasimulation, we would also have to serialize History (the storage class for the delay differential equations) and the solver state.

Then there'd have to be some logic for deciding what things to deserialize and what to make new. Malariasimulation has different solver states, variables and events depending on which interventions are enabled and if we're using individual-based mosquitoes. You could imagine a scenario where people want to load a calibrated run and then add events for rtss (which were never initialised in the original run). I think that's the most complicated case... but would have to think about it.

[cwhittaker1000]

Just wanted to +1 and re-up the above :smile: It'd be great to have the ability to output the current state at the end of a simulation, and then use that as the input to running the model further forwards in time.

I think it'd be really useful both from a time-saving perspective, but also beginning to open up the possibility of inference (via particle-filter based methods), which typically require a lot of stop-starting.

[pwinskill]

Adding some additional notes and a simple example:

To save the state, we would need to record (at least):

• Variables
• Scheduled events
• Solver states (for ODE mosquito model)
• Parameters?

Challenges include:

1. If future scenarios includes new events or require new variables that aren’t initialised in the original run.
2. Coping with difference between ODE and individual mosquito model options.
3. Working with delay-differential equations for mosquito ODE model.
4. Ensuring that parameters defined for future scenarios do not conflict with those in the original run.

I include a simple motivating example below. In this example we have a 2 year simulation and can consider year 1 as being “historical” and therefore only needed to be run once. Two future scenarios are also defined. Both historical and future scenarios have an intervention implemented (MDA = mass drug administration), chosen as it is one that includes a scheduled event.

# Example of motivating case to be able to save and re-load model runs

library(malariasimulation)

# Run the historical scenario (t = 1:365)
## This is the part we would like to run once and save

p_historical <- get_parameters(
  overrides = list(
    human_population = 1000
  )
) |>
  set_drugs(
    drugs = list(SP_AQ_params)
  ) |>
  set_mda(
    drug = 1,
    timesteps = 100,
    coverages =  0.5,
    min_ages = 0,
    max_ages = 80 * 365
  ) |>
  set_equilibrium(
    init_EIR = 5
  )

set.seed(123)
s_historical <- run_simulation(
  timesteps = 365,
  parameters = p_historical
)

s_historical$prevalence <- s_historical$n_detect_730_3650 / s_historical$n_730_3650
plot(
  s_historical$prevalence ~ s_historical$timestep,
  t = "l",
  ylab = "Prevalence", xlab = "Time",
  ylim = c(0, 1), xlim = c(0, 365 * 2)
)
abline(v = 365, lty = 2)

# Run two future scenarios (t = 366:730).
## We would like to initialise these at the end of the historical scenario,
## but currently we have to re-run everything again from t = 0.

p_future1 <- get_parameters(
  overrides = list(
    human_population = 1000
  )
) |>
  set_drugs(
    drugs = list(SP_AQ_params)
  ) |>
  set_mda(
    drug = 1,
    timesteps = c(100, 365 + 100),
    coverages =  c(0.5, 0.1),
    min_ages = c(0, 0),
    max_ages = c(80 * 365, 80 * 365)
  ) |>
  set_equilibrium(
    init_EIR = 5
  )

set.seed(123)
s_future1 <- run_simulation(
  timesteps = 365 * 2,
  parameters = p_future1
)
s_future1$prevalence <- s_future1$n_detect_730_3650 / s_future1$n_730_3650
lines(s_future1$prevalence ~ s_future1$timestep, col = adjustcolor("deeppink", alpha.f = 0.5))

p_future2 <- get_parameters(
  overrides = list(
    human_population = 1000
  )
) |>
  set_drugs(
    drugs = list(SP_AQ_params)
  ) |>
  set_mda(
    drug = 1,
    timesteps = c(100, 365 + 100),
    coverages =  c(0.5, 0.8),
    min_ages = c(0, 0),
    max_ages = c(80 * 365, 80 * 365)
  ) |>
  set_equilibrium(
    init_EIR = 5
  )

set.seed(123)
s_future2 <- run_simulation(
  timesteps = 365 * 2,
  parameters = p_future2
)
s_future2$prevalence <- s_future2$n_detect_730_3650 / s_future2$n_730_3650
lines(s_future2$prevalence ~ s_future2$timestep, col = adjustcolor("dodgerblue", alpha.f = 0.5))