From e534db7e9903b5b64d07d6a77d956def875c509d Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Tue, 13 Feb 2024 17:26:26 +0000 Subject: [PATCH] Allow simulation to be restored with new interventions. (#286) This uses a few of individual's new features that make restoring more flexible. It also fixes a bug when restoring the mosquitto solvers, by correctly restoring the timestep, which is needed to model seasonality. The intervention events use the new `restore = FALSE` flag to make sure their schedule can be modified when resuming. Instead of having the events re-schedule themselves everytime they fire, we setup the entire schedule upfront when initialising the simulation. It adds end-to-end testing of this feature, across a range of scenarios. For each scenario, the outcomes of the simulation with and without restoring are compared and we make sure they are equivalent. --- DESCRIPTION | 2 +- R/RcppExports.R | 4 +- R/compartmental.R | 8 +- R/correlation.R | 9 +- R/events.R | 18 ++-- R/lag.R | 2 +- R/mda_processes.R | 7 -- R/model.R | 17 ++-- R/pev.R | 5 - R/stateful.R | 47 ---------- R/tbv.R | 5 - R/utils.R | 24 +++++ man/CorrelationParameters.Rd | 7 +- src/RcppExports.cpp | 9 +- src/solver.cpp | 3 +- tests/testthat/test-mda.R | 15 --- tests/testthat/test-pev.R | 7 -- tests/testthat/test-resume.R | 174 +++++++++++++++++++++++++++++------ 18 files changed, 217 insertions(+), 146 deletions(-) delete mode 100644 R/stateful.R diff --git a/DESCRIPTION b/DESCRIPTION index 76eb8165..74445497 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -72,7 +72,7 @@ Remotes: Additional_repositories: https://mrc-ide.r-universe.dev Imports: - individual (>= 0.1.15), + individual (>= 0.1.16), malariaEquilibrium (>= 1.0.1), Rcpp, statmod, diff --git a/R/RcppExports.R b/R/RcppExports.R index 678a50bf..ad658f2c 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -53,8 +53,8 @@ solver_get_states <- function(solver) { .Call(`_malariasimulation_solver_get_states`, solver) } -solver_set_states <- function(solver, state) { - invisible(.Call(`_malariasimulation_solver_set_states`, solver, state)) +solver_set_states <- function(solver, t, state) { + invisible(.Call(`_malariasimulation_solver_set_states`, solver, t, state)) } solver_step <- function(solver) { diff --git a/R/compartmental.R b/R/compartmental.R index 7df83a6e..4e16c7c1 100644 --- a/R/compartmental.R +++ b/R/compartmental.R @@ -155,8 +155,8 @@ Solver <- R6::R6Class( save_state = function() { solver_get_states(private$.solver) }, - restore_state = function(state) { - solver_set_states(private$.solver, state) + restore_state = function(t, state) { + solver_set_states(private$.solver, t, state) } ) ) @@ -173,7 +173,7 @@ AquaticMosquitoModel <- R6::R6Class( # state of the ODE is stored separately). We still provide these methods to # conform to the expected interface. save_state = function() { NULL }, - restore_state = function(state) { } + restore_state = function(t, state) { } ) ) @@ -187,7 +187,7 @@ AdultMosquitoModel <- R6::R6Class( save_state = function() { adult_mosquito_model_save_state(self$.model) }, - restore_state = function(state) { + restore_state = function(t, state) { adult_mosquito_model_restore_state(self$.model, state) } ) diff --git a/R/correlation.R b/R/correlation.R index 458a3015..2a59fdbc 100644 --- a/R/correlation.R +++ b/R/correlation.R @@ -125,11 +125,16 @@ CorrelationParameters <- R6::R6Class( }, #' @description Restore the correlation state. + #' #' Only the randomly drawn weights are restored. The object needs to be #' initialized with the same rhos. + #' + #' @param timestep the timestep at which simulation is resumed. This + #' parameter's value is ignored, it only exists to conform to a uniform + #' interface. #' @param state a previously saved correlation state, as returned by the - #' save_state method. - restore_state = function(state) { + #' save_state method. + restore_state = function(timestep, state) { private$.mvnorm <- state$mvnorm } ) diff --git a/R/events.R b/R/events.R index d7ad3416..de0c7c1f 100644 --- a/R/events.R +++ b/R/events.R @@ -1,11 +1,11 @@ create_events <- function(parameters) { events <- list( # MDA events - mda_administer = individual::Event$new(), - smc_administer = individual::Event$new(), + mda_administer = individual::Event$new(restore=FALSE), + smc_administer = individual::Event$new(restore=FALSE), # TBV event - tbv_vaccination = individual::Event$new(), + tbv_vaccination = individual::Event$new(restore=FALSE), # Bednet events throw_away_net = individual::TargetedEvent$new(parameters$human_population) @@ -21,7 +21,7 @@ create_events <- function(parameters) { seq_along(parameters$mass_pev_booster_spacing), function(.) individual::TargetedEvent$new(parameters$human_population) ) - events$mass_pev <- individual::Event$new() + events$mass_pev <- individual::Event$new(restore=FALSE) events$mass_pev_doses <- mass_pev_doses events$mass_pev_boosters <- mass_pev_boosters } @@ -63,16 +63,16 @@ initialise_events <- function(events, variables, parameters) { # Initialise scheduled interventions if (!is.null(parameters$mass_pev_timesteps)) { - events$mass_pev$schedule(parameters$mass_pev_timesteps[[1]] - 1) + events$mass_pev$schedule(parameters$mass_pev_timesteps - 1) } if (parameters$mda) { - events$mda_administer$schedule(parameters$mda_timesteps[[1]] - 1) + events$mda_administer$schedule(parameters$mda_timesteps - 1) } if (parameters$smc) { - events$smc_administer$schedule(parameters$smc_timesteps[[1]] - 1) + events$smc_administer$schedule(parameters$smc_timesteps - 1) } if (parameters$tbv) { - events$tbv_vaccination$schedule(parameters$tbv_timesteps[[1]] - 1) + events$tbv_vaccination$schedule(parameters$tbv_timesteps - 1) } } @@ -158,7 +158,6 @@ attach_event_listeners <- function( if (parameters$mda == 1) { events$mda_administer$add_listener(create_mda_listeners( variables, - events$mda_administer, parameters$mda_drug, parameters$mda_timesteps, parameters$mda_coverages, @@ -174,7 +173,6 @@ attach_event_listeners <- function( if (parameters$smc == 1) { events$smc_administer$add_listener(create_mda_listeners( variables, - events$smc_administer, parameters$smc_drug, parameters$smc_timesteps, parameters$smc_coverages, diff --git a/R/lag.R b/R/lag.R index ca858ac7..de8d1653 100644 --- a/R/lag.R +++ b/R/lag.R @@ -20,7 +20,7 @@ LaggedValue <- R6::R6Class( timeseries_save_state(private$history) }, - restore_state = function(state) { + restore_state = function(t, state) { timeseries_restore_state(private$history, state) } ) diff --git a/R/mda_processes.R b/R/mda_processes.R index 42960be4..72cbe779 100644 --- a/R/mda_processes.R +++ b/R/mda_processes.R @@ -1,6 +1,5 @@ #' @title Create listeners for MDA events #' @param variables the variables available in the model -#' @param administer_event the event schedule for drug administration #' @param drug the drug to administer #' @param timesteps timesteps for each round #' @param coverages the coverage for each round @@ -14,7 +13,6 @@ #' @noRd create_mda_listeners <- function( variables, - administer_event, drug, timesteps, coverages, @@ -78,11 +76,6 @@ create_mda_listeners <- function( variables$drug$queue_update(drug, to_move) variables$drug_time$queue_update(timestep, to_move) } - - # Schedule next round - if (time_index < length(timesteps)) { - administer_event$schedule(timesteps[[time_index + 1]] - timestep) - } } } diff --git a/R/model.R b/R/model.R index def6db23..67f79758 100644 --- a/R/model.R +++ b/R/model.R @@ -135,16 +135,19 @@ run_resumable_simulation <- function( lagged_eir <- create_lagged_eir(variables, solvers, parameters) lagged_infectivity <- create_lagged_infectivity(variables, parameters) - stateful_objects <- unlist(list( + stateful_objects <- list( RandomState$new(restore_random_state), correlations, vector_models, solvers, lagged_eir, - lagged_infectivity)) + lagged_infectivity) if (!is.null(initial_state)) { - restore_state(initial_state$malariasimulation, stateful_objects) + individual::restore_object_state( + initial_state$timesteps, + stateful_objects, + initial_state$malariasimulation) } individual_state <- individual::simulation_loop( @@ -161,16 +164,16 @@ run_resumable_simulation <- function( timesteps ), variables = variables, - events = unlist(events), + events = events, timesteps = timesteps, state = initial_state$individual, restore_random_state = restore_random_state ) final_state <- list( - timesteps=timesteps, - individual=individual_state, - malariasimulation=save_state(stateful_objects) + timesteps = timesteps, + individual = individual_state, + malariasimulation = individual::save_object_state(stateful_objects) ) data <- renderer$to_dataframe() diff --git a/R/pev.R b/R/pev.R index bc27fc9f..e460e5c3 100644 --- a/R/pev.R +++ b/R/pev.R @@ -130,11 +130,6 @@ create_mass_pev_listener <- function( parameters, events$mass_pev_doses ) - if (time_index < length(parameters$mass_pev_timesteps)) { - events$mass_pev$schedule( - parameters$mass_pev_timesteps[[time_index + 1]] - timestep - ) - } } } diff --git a/R/stateful.R b/R/stateful.R deleted file mode 100644 index 05c614a0..00000000 --- a/R/stateful.R +++ /dev/null @@ -1,47 +0,0 @@ -#' @title Save the state of a list of \emph{stateful objects} -#' @description The state of each element is saved and stored into a single -#' object, representing them in a way that can be exported and re-used later. -#' @param objects A list of stateful objects to be saved. Stateful objects are -#' instances of R6 classes with a pair of \code{save_state} and -#' \code{restore_state} methods. -#' @noRd -save_state <- function(objects) { - lapply(objects, function(o) o$save_state()) -} - -#' @title Restore the state of a collection of stateful objects -#' @description This is the counterpart of \code{save_state}. Calling it -#' restores the collection of objects back into their original state. -#' @param state A state object, as returned by the \code{save_state} function. -#' @param objects A collection of stateful objects to be restored. -#' @noRd -restore_state <- function(state, objects) { - stopifnot(length(state) == length(objects)) - for (i in seq_along(state)) { - objects[[i]]$restore_state(state[[i]]) - } -} - -#' @title a placeholder class to save the random number generator class. -#' @description the class integrates with the simulation loop to save and -#' restore the random number generator class when appropriate. -#' @noRd -RandomState <- R6::R6Class( - 'RandomState', - private = list( - restore_random_state = NULL - ), - public = list( - initialize = function(restore_random_state) { - private$restore_random_state <- restore_random_state - }, - save_state = function() { - random_save_state() - }, - restore_state = function(state) { - if (private$restore_random_state) { - random_restore_state(state) - } - } - ) -) diff --git a/R/tbv.R b/R/tbv.R index 27c29dfb..d0f75de1 100644 --- a/R/tbv.R +++ b/R/tbv.R @@ -75,11 +75,6 @@ create_tbv_listener <- function(variables, events, parameters, correlations, ren to_vaccinate ) } - if (time_index < length(parameters$tbv_timesteps)) { - events$tbv_vaccination$schedule( - parameters$tbv_timesteps[[time_index + 1]] - timestep - ) - } } } diff --git a/R/utils.R b/R/utils.R index f25244bd..2262c939 100644 --- a/R/utils.R +++ b/R/utils.R @@ -64,3 +64,27 @@ rtexp <- function(n, m, t) { itexp(runif(n), m, t) } match_timestep <- function(ts, t) { min(sum(ts <= t), length(ts)) } + +#' @title a placeholder class to save the random number generator class. +#' @description the class integrates with the simulation loop to save and +#' restore the random number generator class when appropriate. +#' @noRd +RandomState <- R6::R6Class( + 'RandomState', + private = list( + restore_random_state = NULL + ), + public = list( + initialize = function(restore_random_state) { + private$restore_random_state <- restore_random_state + }, + save_state = function() { + random_save_state() + }, + restore_state = function(t, state) { + if (private$restore_random_state) { + random_restore_state(state) + } + } + ) +) diff --git a/man/CorrelationParameters.Rd b/man/CorrelationParameters.Rd index 2898aee5..8bb43d70 100644 --- a/man/CorrelationParameters.Rd +++ b/man/CorrelationParameters.Rd @@ -121,15 +121,20 @@ Save the correlation state. \if{latex}{\out{\hypertarget{method-CorrelationParameters-restore_state}{}}} \subsection{Method \code{restore_state()}}{ Restore the correlation state. + Only the randomly drawn weights are restored. The object needs to be initialized with the same rhos. \subsection{Usage}{ -\if{html}{\out{
}}\preformatted{CorrelationParameters$restore_state(state)}\if{html}{\out{
}} +\if{html}{\out{
}}\preformatted{CorrelationParameters$restore_state(timestep, state)}\if{html}{\out{
}} } \subsection{Arguments}{ \if{html}{\out{
}} \describe{ +\item{\code{timestep}}{the timestep at which simulation is resumed. This +parameter's value is ignored, it only exists to conform to a uniform +interface.} + \item{\code{state}}{a previously saved correlation state, as returned by the save_state method.} } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 439cb12b..e2a6b9d6 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -203,13 +203,14 @@ BEGIN_RCPP END_RCPP } // solver_set_states -void solver_set_states(Rcpp::XPtr solver, std::vector state); -RcppExport SEXP _malariasimulation_solver_set_states(SEXP solverSEXP, SEXP stateSEXP) { +void solver_set_states(Rcpp::XPtr solver, double t, std::vector state); +RcppExport SEXP _malariasimulation_solver_set_states(SEXP solverSEXP, SEXP tSEXP, SEXP stateSEXP) { BEGIN_RCPP Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Rcpp::XPtr >::type solver(solverSEXP); + Rcpp::traits::input_parameter< double >::type t(tSEXP); Rcpp::traits::input_parameter< std::vector >::type state(stateSEXP); - solver_set_states(solver, state); + solver_set_states(solver, t, state); return R_NilValue; END_RCPP } @@ -363,7 +364,7 @@ static const R_CallMethodDef CallEntries[] = { {"_malariasimulation_rainfall", (DL_FUNC) &_malariasimulation_rainfall, 5}, {"_malariasimulation_exponential_process_cpp", (DL_FUNC) &_malariasimulation_exponential_process_cpp, 2}, {"_malariasimulation_solver_get_states", (DL_FUNC) &_malariasimulation_solver_get_states, 1}, - {"_malariasimulation_solver_set_states", (DL_FUNC) &_malariasimulation_solver_set_states, 2}, + {"_malariasimulation_solver_set_states", (DL_FUNC) &_malariasimulation_solver_set_states, 3}, {"_malariasimulation_solver_step", (DL_FUNC) &_malariasimulation_solver_step, 1}, {"_malariasimulation_create_timeseries", (DL_FUNC) &_malariasimulation_create_timeseries, 2}, {"_malariasimulation_timeseries_at", (DL_FUNC) &_malariasimulation_timeseries_at, 3}, diff --git a/src/solver.cpp b/src/solver.cpp index 0a5f7401..4c3f182f 100644 --- a/src/solver.cpp +++ b/src/solver.cpp @@ -48,7 +48,8 @@ std::vector solver_get_states(Rcpp::XPtr solver) { } //[[Rcpp::export]] -void solver_set_states(Rcpp::XPtr solver, std::vector state) { +void solver_set_states(Rcpp::XPtr solver, double t, std::vector state) { + solver->t = t; solver->state = state; } diff --git a/tests/testthat/test-mda.R b/tests/testthat/test-mda.R index bf92ae3a..1ee440f6 100644 --- a/tests/testthat/test-mda.R +++ b/tests/testthat/test-mda.R @@ -81,11 +81,8 @@ test_that('MDA moves the diseased and non-diseased population correctly', { drug = mock_double(c(1, 2, 1, 2)) ) - events$mda_administer <- mock_event(events$mda_administer) - listener <- create_mda_listeners( variables, - events$mda_administer, parameters$mda_drug, parameters$mda_timesteps, parameters$mda_coverages, @@ -141,12 +138,6 @@ test_that('MDA moves the diseased and non-diseased population correctly', { 50, c(3, 4) ) - - mockery::expect_args( - events$mda_administer$schedule, - 1, - 100 - ) }) test_that('MDA moves the diseased and non-diseased population correctly - second round, varying age range', { @@ -176,11 +167,8 @@ test_that('MDA moves the diseased and non-diseased population correctly - second drug = mock_double(c(1, 2, 1, 2)) ) - events$mda_administer <- mock_event(events$mda_administer) - listener <- create_mda_listeners( variables, - events$mda_administer, parameters$mda_drug, parameters$mda_timesteps, parameters$mda_coverages, @@ -265,11 +253,8 @@ test_that('MDA ignores non-detectable asymptomatics', { drug = mock_double(c(1, 2, 1, 2)) ) - events$mda_administer <- mock_event(events$mda_administer) - listener <- create_mda_listeners( variables, - events$mda_administer, parameters$mda_drug, parameters$mda_timesteps, parameters$mda_coverages, diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R index 88b95b35..09c4cf7d 100644 --- a/tests/testthat/test-pev.R +++ b/tests/testthat/test-pev.R @@ -203,7 +203,6 @@ test_that('Mass vaccinations update vaccination time', { c(-1, -1, -1, 50, 50) ) - events$mass_pev <- mock_event(events$mass_pev) events$mass_pev_doses <- lapply(events$mass_pev_doses, mock_event) listener <- create_mass_pev_listener( @@ -241,12 +240,6 @@ test_that('Mass vaccinations update vaccination time', { c(1, 3), parameters$pev_doses[[3]] ) - - mockery::expect_args( - events$mass_pev$schedule, - 1, - 365 - ) }) test_that('Mass vaccinations ignore EPI individuals', { diff --git a/tests/testthat/test-resume.R b/tests/testthat/test-resume.R index 83bada35..4d4a572d 100644 --- a/tests/testthat/test-resume.R +++ b/tests/testthat/test-resume.R @@ -1,41 +1,161 @@ -test_that('Simulation can be resumed', { - initial_timesteps <- 50 - total_timesteps <- 100 - - parameters <- get_parameters() +#' Test simulation saving and restoring for a given parameter set. +#' +#' This function runs the simulation twice. A first, continuous and uninterrupted +#' run of the simulation is used as a reference. The second run is split into +#' two phases. Between the two phases, the simulation state is saved and +#' restored. Optionally, the initial warmup phase can use a different set of +#' parameters, by specifying a value for warmup_parameters. +test_resume <- function( + parameters, + timesteps = 200, + warmup_parameters = parameters, + warmup_timesteps = 50 + ) { + set.seed(123) + uninterrupted_run <- run_simulation(timesteps, parameters=parameters) - set.seed(1) - first_phase <- run_resumable_simulation(initial_timesteps, parameters) + set.seed(123) + first_phase <- run_resumable_simulation(warmup_timesteps, warmup_parameters) second_phase <- run_resumable_simulation( - total_timesteps, + timesteps, parameters, initial_state=first_phase$state, restore_random_state=TRUE) - set.seed(1) - uninterrupted_run <- run_simulation(total_timesteps, parameters=parameters) + expect_equal(nrow(first_phase$data), warmup_timesteps) + expect_equal(nrow(second_phase$data), timesteps - warmup_timesteps) + + expect_mapequal( + second_phase$data, + uninterrupted_run[-(1:warmup_timesteps),]) + + invisible(second_phase$data) +} + +test_that('Simulation can be resumed', { + test_resume(get_parameters(overrides=list( + individual_mosquitoes = FALSE, + model_seasonality = TRUE + ))) + test_resume(get_parameters(overrides=list( + individual_mosquitoes = TRUE, + model_seasonality = TRUE + ))) + test_resume(get_parameters(overrides=list( + individual_mosquitoes = FALSE, + model_seasonality = TRUE + ))) + test_resume(get_parameters(overrides=list( + individual_mosquitoes = TRUE, + model_seasonality = TRUE + ))) +}) + +test_that('PEV intervention can be added when resuming', { + set_default_mass_pev <- function(parameters, timesteps) { + n <- length(timesteps) + set_mass_pev( + parameters, + profile = rtss_profile, + timesteps = timesteps, + coverages = rep(0.5, n), + min_wait = 5, + min_ages = 365*10, + max_ages = 365*60, + booster_spacing = NULL, + booster_coverage = NULL, + booster_profile = NULL) + } + base <- get_parameters(overrides=list(pev_doses=c(0, 45, 90))) + + data <- test_resume( + warmup_parameters = base, + parameters = base %>% set_default_mass_pev(100)) + expect_equal(data[data$n_pev_mass_dose_1 > 0, "timestep"], 100) + expect_equal(data[data$n_pev_mass_dose_2 > 0, "timestep"], 145) + expect_equal(data[data$n_pev_mass_dose_3 > 0, "timestep"], 190) + + # Add a second mass PEV intervention when resuming the simulation. + data <- test_resume( + warmup_parameters = base %>% set_default_mass_pev(25), + parameters = base %>% set_default_mass_pev(c(25, 100))) - expect_equal(nrow(first_phase$data), initial_timesteps) - expect_equal(nrow(second_phase$data), total_timesteps - initial_timesteps) - expect_equal(rbind(first_phase$data, second_phase$data), uninterrupted_run) + # The first dose, at time step 25, happens during the warmup and is not + # returned by test_resume, hence why we don't see it in the data. Follow-up + # doses do show up, even though they we scheduled during warmup. + expect_equal(data[data$n_pev_mass_dose_1 > 0, "timestep"], c(100)) + expect_equal(data[data$n_pev_mass_dose_2 > 0, "timestep"], c(70, 145)) + expect_equal(data[data$n_pev_mass_dose_3 > 0, "timestep"], c(115, 190)) }) -test_that('Intervention parameters can be changed when resuming', { - initial_timesteps <- 50 - total_timesteps <- 100 - tbv_timesteps <- 70 +test_that("TBV intervention can be added when resuming", { + set_default_tbv <- function(parameters, timesteps) { + set_tbv( + parameters, + timesteps=timesteps, + coverage=rep(1, length(timesteps)), + ages=5:60) + } + + base <- get_parameters() + + data <- test_resume( + warmup_parameters = base, + parameters = base %>% set_default_tbv(100)) + expect_equal(data[!is.na(data$n_vaccinated_tbv), "timestep"], 100) + + data <- test_resume( + warmup_parameters = base %>% set_default_tbv(25), + parameters = base %>% set_default_tbv(c(25, 100))) + expect_equal(data[!is.na(data$n_vaccinated_tbv), "timestep"], 100) +}) + +test_that("MDA intervention can be added when resuming", { + set_default_mda <- function(parameters, timesteps) { + parameters %>% set_drugs(list(SP_AQ_params)) %>% set_mda( + drug = 1, + timesteps = timesteps, + coverages = rep(0.8, length(timesteps)), + min_ages = rep(0, length(timesteps)), + max_ages = rep(60*365, length(timesteps))) + } + + base <- get_parameters() + + data <- test_resume( + warmup_parameters = base, + parameters = base %>% set_default_mda(100)) + expect_equal(data[data$n_mda_treated > 0, "timestep"], 100) + + data <- test_resume( + warmup_parameters = base %>% set_default_mda(25), + parameters = base %>% set_default_mda(c(25, 100))) + expect_equal(data[data$n_mda_treated > 0, "timestep"], 100) +}) - # Because of how event scheduling works, we must enable TBV even in the inital phase. - # We set a coverage of 0 to act as-if it was disabled. - initial_parameters <- get_parameters() |> set_tbv(timesteps=tbv_timesteps, coverage=0, ages=5:60) +test_that("Bednets intervention can be added when resuming", { + set_default_bednets <- function(parameters, timesteps) { + n <- length(timesteps) + set_bednets( + parameters, + timesteps = timesteps, + coverages = rep(0.5, n), + retention = 25, + dn0 = matrix(rep(0.533, n), ncol=1), + rn = matrix(rep(0.56, n), ncol=1), + rnm = matrix(rep(0.24, n), ncol=1), + gamman = rep(2.64 * 365, n)) + } - tbv_parameters <- initial_parameters |> - set_tbv(timesteps=tbv_timesteps, coverage=1, ages=5:60) + base <- get_parameters() - initial_run <- run_resumable_simulation(initial_timesteps, initial_parameters) - control_run <- run_resumable_simulation(total_timesteps, initial_parameters, initial_state = initial_run$state) - tbv_run <- run_resumable_simulation(total_timesteps, tbv_parameters, initial_state = initial_run$state) + data <- test_resume( + warmup_parameters = base, + parameters = base %>% set_default_bednets(100)) + expect_equal(data[diff(data$n_use_net) > 0, "timestep"], 100) - expect_equal(control_run$data$n_vaccinated_tbv[tbv_timesteps - initial_timesteps], 0) - expect_gt(tbv_run$data$n_vaccinated_tbv[tbv_timesteps - initial_timesteps], 0) + data <- test_resume( + warmup_parameters = base %>% set_default_bednets(25), + parameters = base %>% set_default_bednets(c(25, 100))) + expect_equal(data[diff(data$n_use_net) > 0, "timestep"], 100) })