From 050326f024aeb46d8dfa3c58c0da7285a6536968 Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Tue, 13 Feb 2024 16:42:15 +0000 Subject: [PATCH 1/2] Add a `restore` flag to the Event constructor. By default, when restoring the simulation state all previous schedule on events are cleared and restored from the saved state. This goes against use cases that wish to resume a simulation with different intervention schedules to compare effects. In those use cases, a different initialization sequence is used when creating the simulation, and we do not want that to be cleared and overwritten. The new `restore` flag, when set to false, overrides this default behaviour and the state of an Event is (mostly) unaffected by a restore. Thanks to this, a new event schedule, that is unrelated to the schedule of the original run, can be configured. --- R/RcppExports.R | 4 +- R/event.R | 14 +++-- inst/include/Event.h | 16 +++++- src/RcppExports.cpp | 9 ++-- src/event.cpp | 4 +- vignettes/Checkpoint.Rmd | 107 +++++++++++++++++++++------------------ 6 files changed, 90 insertions(+), 64 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 18b51bd..f67a6e6 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -165,8 +165,8 @@ double_variable_queue_shrink_bitset <- function(variable, index) { invisible(.Call(`_individual_double_variable_queue_shrink_bitset`, variable, index)) } -create_event <- function() { - .Call(`_individual_create_event`) +create_event <- function(restoreable) { + .Call(`_individual_create_event`, restoreable) } create_targeted_event <- function(size) { diff --git a/R/event.R b/R/event.R index c7d3443..dc60cd5 100644 --- a/R/event.R +++ b/R/event.R @@ -38,18 +38,24 @@ EventBase <- R6Class( #' @export Event <- R6Class( 'Event', - inherit=EventBase, + inherit = EventBase, public = list( #' @description Initialise an Event. - initialize = function() { - self$.event <- create_event() + #' @param restore if true, the schedule of this event is restored when restoring from a saved + #' simulation. + initialize = function(restore = TRUE) { + self$.event <- create_event(restore) }, #' @description Schedule this event to occur in the future. #' @param delay the number of time steps to wait before triggering the event, #' can be a scalar or a vector of values for events that should be triggered #' multiple times. - schedule = function(delay) event_schedule(self$.event, delay), + schedule = function(delay) { + if (!is.null(delay)) { + event_schedule(self$.event, delay) + } + }, #' @description Stop a future event from triggering. clear_schedule = function() event_clear_schedule(self$.event), diff --git a/inst/include/Event.h b/inst/include/Event.h index 9d32d9b..73596e7 100644 --- a/inst/include/Event.h +++ b/inst/include/Event.h @@ -70,8 +70,10 @@ inline size_t EventBase::get_time() const { class Event : public EventBase { std::set simple_schedule; + bool restoreable; public: + Event(bool restoreable); virtual ~Event() = default; virtual void process(Rcpp::XPtr listener); @@ -90,6 +92,9 @@ inline void Event::process(Rcpp::XPtr listener) { (*listener)(get_time()); } +inline Event::Event(bool restoreable) : restoreable(restoreable) { +} + //' @title should first event fire on this timestep? inline bool Event::should_trigger() { if (simple_schedule.empty()) { @@ -124,8 +129,15 @@ inline std::vector Event::checkpoint() { //' @title restore this event's state from a previous checkpoint inline void Event::restore(size_t time, std::vector schedule) { t = time; - simple_schedule.clear(); - simple_schedule.insert(schedule.begin(), schedule.end()); + if (restoreable) { + simple_schedule.clear(); + simple_schedule.insert(schedule.begin(), schedule.end()); + } else { + // We don't restore the event, but it is possible that the resume time + // is beyond some already scheduled timesteps. These need to be removed. + auto it = simple_schedule.lower_bound(time); + simple_schedule.erase(simple_schedule.begin(), it); + } } //' @title a targeted event in the simulation diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 7cb6374..0106ad3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -512,12 +512,13 @@ BEGIN_RCPP END_RCPP } // create_event -Rcpp::XPtr create_event(); -RcppExport SEXP _individual_create_event() { +Rcpp::XPtr create_event(bool restoreable); +RcppExport SEXP _individual_create_event(SEXP restoreableSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; - rcpp_result_gen = Rcpp::wrap(create_event()); + Rcpp::traits::input_parameter< bool >::type restoreable(restoreableSEXP); + rcpp_result_gen = Rcpp::wrap(create_event(restoreable)); return rcpp_result_gen; END_RCPP } @@ -1487,7 +1488,7 @@ static const R_CallMethodDef CallEntries[] = { {"_individual_double_variable_queue_extend", (DL_FUNC) &_individual_double_variable_queue_extend, 2}, {"_individual_double_variable_queue_shrink", (DL_FUNC) &_individual_double_variable_queue_shrink, 2}, {"_individual_double_variable_queue_shrink_bitset", (DL_FUNC) &_individual_double_variable_queue_shrink_bitset, 2}, - {"_individual_create_event", (DL_FUNC) &_individual_create_event, 0}, + {"_individual_create_event", (DL_FUNC) &_individual_create_event, 1}, {"_individual_create_targeted_event", (DL_FUNC) &_individual_create_targeted_event, 1}, {"_individual_event_base_tick", (DL_FUNC) &_individual_event_base_tick, 1}, {"_individual_event_base_get_timestep", (DL_FUNC) &_individual_event_base_get_timestep, 1}, diff --git a/src/event.cpp b/src/event.cpp index 522291f..ef8fd7a 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -9,8 +9,8 @@ #include "utils.h" //[[Rcpp::export]] -Rcpp::XPtr create_event() { - return Rcpp::XPtr(new Event(), true); +Rcpp::XPtr create_event(bool restoreable) { + return Rcpp::XPtr(new Event(restoreable), true); } //[[Rcpp::export]] diff --git a/vignettes/Checkpoint.Rmd b/vignettes/Checkpoint.Rmd index 73ce393..5ba9e14 100644 --- a/vignettes/Checkpoint.Rmd +++ b/vignettes/Checkpoint.Rmd @@ -26,7 +26,8 @@ knitr::asis_output( 2. [Usage](#usage) 3. [Example](#example) 4. [Caveats](#caveats) - 1. [Restoring random number generator state](#rng)" + 1. [Restoring random number generator state](#rng) + 2. [Saving and restoring events](#events)" ) ``` @@ -41,14 +42,14 @@ Individual allows the user to run a simulation for a number of time steps, save The typical way to use this feature is to define a simulation function which creates all the relevant simulation data and then calls `simulation_loop`. The function we define takes in an optional `state` parameter that is passed through to `simulation_loop`. ```{r} -run_simulation <- function(timesteps, state=NULL) { +run_simulation <- function(timesteps, state = NULL) { health <- CategoricalVariable$new(c("S", "I"), rep("S", 10)) process <- bernoulli_process(health, "S", "I", 0.01) simulation_loop( - variables=list(health), - processes=list(process), - timesteps=timesteps, - state=state) + variables = list(health), + processes = list(process), + timesteps = timesteps, + state = state) } ``` @@ -59,7 +60,7 @@ state <- run_simulation(timesteps = 50) Finally, the simulation is resumed with a larger number of time steps, passing in the state object as an argument. The `timesteps` argument refers to the total number of time steps, including both the original simulation run and the new one. In this case, `run_simulation` will only simulate 50 extra steps. Before running the actual simulation, `simulation_loop` will reload the simulation state from its argument, overwriting any values we had set when initializing the variables. ```{r, results='hide'} -run_simulation(timesteps = 100, state=state) +run_simulation(timesteps = 100, state = state) ``` ## Practical example {#example} @@ -76,7 +77,7 @@ make_variables <- function(N, I0) { vaccinated <- CategoricalVariable$new(categories = c("Y", "N"), initial_values = rep("N", N)) - list(health=health, vaccinated=vaccinated) + list(health = health, vaccinated = vaccinated) } ``` @@ -100,19 +101,18 @@ make_infection_process <- function(health, vaccinated, N, beta, vaccine_efficacy } ``` -For a while at the start of the simulation no vaccination takes place. Only after a number of time steps, determined by the `vaccination_start` parameter, does the intervention begin. Periodically, an event will fire and every individual has a fixed probability of becoming vaccinated. If `vaccination_start` is `NULL`, the intervention never begins. +Vaccination happens at fixed mass-vaccination event, according to a pre-determined schedule defined by the `vaccination_times` parameter. Each vaccination event affects a percentage of the population, as determined by the `vaccination_coverage` parameter. The coverage may vary from one instance to the next. ```{r} -make_vaccination_event <- function(vaccinated, vaccination_start, vaccination_interval, vaccination_rate) { - e <- Event$new() +make_vaccination_event <- function(vaccinated, vaccination_times, vaccination_coverage) { + e <- Event$new(restore=FALSE) + e$schedule(vaccination_times - 1) e$add_listener(function(t) { - vaccinated$queue_update(value = "Y", - vaccinated$get_index_of("N")$sample(vaccination_rate)) - e$schedule(vaccination_interval) + index <- which(vaccination_times == t) + coverage <- vaccination_coverage[[index]] + targets <- vaccinated$get_index_of("N")$sample(coverage) + vaccinated$queue_update(value = "Y", targets) }) - if (!is.null(vaccination_start)) { - e$schedule(vaccination_start - 1) - } e } ``` @@ -123,13 +123,12 @@ We will define our simulation as a function, taking the simulation parameters as run_simulation <- function( steps, N = 1e3, - I0 = 5, - beta = 0.1, # S -> I + I0 = 10, + beta = 0.08, # S -> I gamma = 0.05, # I -> R - xi = 0.03, # R -> S - vaccination_start = NULL, - vaccination_interval = 10, - vaccination_rate = 0.05, # N -> Y + xi = 0.02, # R -> S + vaccination_times = numeric(0), + vaccination_coverage = rep(0.2, length(vaccination_times)), # N -> Y vaccine_efficacy = 1, ...) { @@ -141,8 +140,7 @@ run_simulation <- function( recovery_process <- bernoulli_process(variables$health, "I", "R", gamma) return_process <- bernoulli_process(variables$health, "R", "S", xi) vaccination_event <- make_vaccination_event( - variables$vaccinated, vaccination_start, vaccination_interval, vaccination_rate) - + variables$vaccinated, vaccination_times, vaccination_coverage) renderer <- Render$new(timesteps = steps) health_render_process <- categorical_count_renderer_process( renderer = renderer, @@ -163,86 +161,96 @@ run_simulation <- function( timesteps = steps, ...) - list(result=renderer$to_dataframe(), state=final_state) + list(result = renderer$to_dataframe(), state = final_state) } ``` We will start by running and plotting our baseline simulation, with the intervention disabled. ```{r} -data <- run_simulation(steps=1500)$result +data <- run_simulation(steps = 1500)$result colours <- c("royalblue3","firebrick3","darkorchid3") matplot( - x=data["timestep"], - y=data[c("S_count","I_count", "R_count")], - xlab="Time", ylab="Count", - type="l", lwd=2, lty = 1, col = colours + x = data["timestep"], + y = data[c("S_count","I_count", "R_count")], + xlab = "Time", ylab = "Population count", + type = "l", lwd = 2, lty = 1, col = colours ) legend( x = "topright", pch = rep(16,3), col = colours, - legend = c("S", "I", "R"), cex = 1.5, + legend = c("S", "I", "R"), cex = 1, bg='white' ) ``` -We see that the simulation takes some time to settle from its initial parameters to its steady-state conditions. -We will now enable the vaccine intervention, but only starting at a point after the simulation has settled, for example at `t=500`. +We see that the simulation takes some time to settle from its initial parameters to its steady-state conditions. We will now enable the vaccine intervention, but only starting at a point after the simulation has settled, for example at `t=500`, and every 250 steps after that. ```{r} -data <- run_simulation(steps=1500, vaccination_start = 500, vaccine_efficacy = 1)$result +vaccination_times <- seq(500,1500,250) +data <- run_simulation(steps = 1500, vaccination_times = vaccination_times)$result colours <- c("royalblue3","firebrick3","darkorchid3") matplot( - x=data["timestep"], - y=data[c("S_count","I_count", "R_count")], - xlab="Time", ylab="Count", - type="l", lwd=2, lty = 1, col = colours + x = data["timestep"], + y = data[c("S_count","I_count", "R_count")], + xlab = "Time", ylab = "Count", + type = "l", lwd = 2, lty = 1, col = colours ) +abline(v = vaccination_times, col = "snow4") legend( x = "topright", pch = rep(16,3), col = colours, - legend = c("S", "I", "R"), cex = 1.5, + legend = c("S", "I", "R"), + cex = 1, bg='white' ) ``` The simulation above clearly shows the effect of the vaccination campaign, starting at `t=500`. However, it made the optimistic assumption of a 100% vaccine efficacy. We wish to run the simulation again but with varying levels of efficacy, in order the compare its impact. -While we could run the code above many times over, each simulation would repeat the first 499 timesteps, despite the result being identical each time. Instead we start by running only these timesteps, and saving the result. We do need to specify the start of the intervention, as it is necessary to schedule the first vaccination event. However the details of the intervention (ie. `vaccine_efficacy`) are irrelevant and can be omitted. +While we could run the code above many times over, each simulation would repeat the first 499 timesteps, despite the result being identical each time. Instead we start by running only these timesteps, and saving the result. The details of the intervention (ie. `vaccination_times` and `vaccine_efficacy`) are irrelevant and can be omitted. ```{r} -initial <- run_simulation(steps=499, vaccination_start = 500) +initial <- run_simulation(steps = 499) ``` From this initial result, we can resume the simulation, but using different values of vaccine efficacy each time. We also include a control simulation, in which no vaccination takes place. Each of these simulation will skip the first 499 steps and only run the next 1001 time steps. ```{r} -control <- run_simulation(steps=1500, vaccination_start = 500, vaccine_efficacy=0.0, state=initial$state) -vaccine30 <- run_simulation(steps=1500, vaccination_start = 500, vaccine_efficacy=0.3, state=initial$state) -vaccine50 <- run_simulation(steps=1500, vaccination_start = 500, vaccine_efficacy=0.5, state=initial$state) -vaccine100 <- run_simulation(steps=1500, vaccination_start = 500, vaccine_efficacy=1.0, state=initial$state) +control <- run_simulation(steps = 1500, state = initial$state) +vaccine50 <- run_simulation( + steps = 1500, vaccination_times = vaccination_times, vaccine_efficacy = 0.5, + state = initial$state) +vaccine80 <- run_simulation( + steps = 1500, vaccination_times = vaccination_times, vaccine_efficacy = 0.8, + state = initial$state) +vaccine100 <- run_simulation( + steps = 1500, vaccination_times = vaccination_times, vaccine_efficacy = 1.0, + state = initial$state) ``` Finally we aggregate and plot the results from all these simulations. We also need to include the data from our initial run, which we will plot the same colour as our control simulation. ```{r} -colours <- c("royalblue3","firebrick3","darkorchid3", "seagreen3") +colours <- c("royalblue3", "firebrick3","darkorchid3", "seagreen3") # Pad initial out to ensure it has the same shape as other series. initial$result[500:1500,] <- NA matplot( data.frame( initial$result[,"I_count"], - vaccine30$result[,"I_count"], vaccine50$result[,"I_count"], + vaccine80$result[,"I_count"], vaccine100$result[,"I_count"], control$result[,"I_count"]), xlab = "Time", ylab = "Infected count", type = "l", lwd = 1.5, lty = 1, col = colours, ) +abline(v = vaccination_times, col = "snow4") legend( x = "topright", pch = rep(16,3), col = colours, - legend = c("Control", "30%", "50%", "100%"), cex = 1.5, + legend = c("Control", "50%", "80%", "100%"), + cex = 1, bg='white' ) ``` @@ -253,7 +261,6 @@ Saving and restoring the simulation state comes with a number of caveats. - All simulation state must be represented as objects managed by individual. Any state maintained externally will not be saved nor restored. - The state object's structure is not stable and is expected to change. One should not expect to serialize the state to disk and have it work with future versions of the individual package. - The simulation must be re-created in an identical way. Variables and events may not be added or removed, variable sizes must remain constant, the list of categories in a `CategoricalVariable` cannot be modified, etc. The order of variables and events passed to the `run_simulation` function must remain stable. -- If an event is scheduled before the checkpoint, the time at which it will execute cannot be changed when resuming, even if that time is in the future. For example in the SIRS model above, we would not be able to resume the simulation with different values for `vaccination_start`; changing that parameter would have no effect. While parameters of the simulation can be changed between the initial run and the subsequent runs (as demonstrated with the `vaccine_efficacy` parameter above), in general you should not modify parameters that would have been already had an impact on the first part of the simulation. Doing so would produce results that can only be produced through checkpoint and resume, and not as a single simulation. From 9e79ea3222abf387053d687a8b5b814294432ceb Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Wed, 14 Feb 2024 16:39:18 +0000 Subject: [PATCH 2/2] Update tests and vignette. --- tests/testthat/test-events.R | 52 +++++++++++++++++++++++++++++++++++- vignettes/Checkpoint.Rmd | 51 ++++++++++++++++++++++++++++++----- 2 files changed, 96 insertions(+), 7 deletions(-) diff --git a/tests/testthat/test-events.R b/tests/testthat/test-events.R index af371e6..2045b2e 100644 --- a/tests/testthat/test-events.R +++ b/tests/testthat/test-events.R @@ -127,6 +127,8 @@ test_that("events can be saved and restored", { old_event$.timestep(), old_event$.checkpoint()) + expect_equal(new_event$.timestep(), 3) + #time = 3 new_event$.process() mockery::expect_called(listener, 1) @@ -164,7 +166,7 @@ test_that("events are cleared when restored", { mockery::expect_called(listener, 0) new_event$.tick() - #time=2 + #time=3 new_event$.process() mockery::expect_called(listener, 0) new_event$.tick() @@ -176,6 +178,54 @@ test_that("events are cleared when restored", { new_event$.tick() }) +test_that("event is not cleared when restore is disabled", { + old_event <- Event$new() + old_event$schedule(3) # t=4 + + # time = 1 + old_event$.process() + old_event$.tick() + + # time = 2 + old_event$.process() + old_event$.tick() + + listener <- mockery::mock() + new_event <- Event$new(restore=FALSE) + new_event$add_listener(listener) + new_event$schedule(1) # t=2, this is before the restore timestep and will not fire. + new_event$schedule(5) # t=6, this will trigger + + # This should restore the timestep, but not any of the old event's schedule. + new_event$.restore( + old_event$.timestep(), + old_event$.checkpoint()) + + expect_equal(new_event$.timestep(), 3) + + # time = 3 + new_event$.process() + mockery::expect_called(listener, 0) + new_event$.tick() + + # time = 4 + # The saved event from the `old_event` is not restored. + new_event$.process() + mockery::expect_called(listener, 0) + new_event$.tick() + + # time = 5 + new_event$.process() + mockery::expect_called(listener, 0) + new_event$.tick() + + # time = 6 + new_event$.process() + mockery::expect_called(listener, 1) + mockery::expect_args(listener, 1, t = 6) + new_event$.tick() +}) + test_that("empty event never triggers", { event <- Event$new() listener <- mockery::mock() diff --git a/vignettes/Checkpoint.Rmd b/vignettes/Checkpoint.Rmd index 5ba9e14..4d0fae7 100644 --- a/vignettes/Checkpoint.Rmd +++ b/vignettes/Checkpoint.Rmd @@ -103,9 +103,11 @@ make_infection_process <- function(health, vaccinated, N, beta, vaccine_efficacy Vaccination happens at fixed mass-vaccination event, according to a pre-determined schedule defined by the `vaccination_times` parameter. Each vaccination event affects a percentage of the population, as determined by the `vaccination_coverage` parameter. The coverage may vary from one instance to the next. +Because we intend on saving and restoring the simulation with various vaccination schedules, we pass `restore = FALSE` to the `Event` constructor. See the [Saving and restoring events](#events) section for details. + ```{r} make_vaccination_event <- function(vaccinated, vaccination_times, vaccination_coverage) { - e <- Event$new(restore=FALSE) + e <- Event$new(restore = FALSE) e$schedule(vaccination_times - 1) e$add_listener(function(t) { index <- which(vaccination_times == t) @@ -127,7 +129,7 @@ run_simulation <- function( beta = 0.08, # S -> I gamma = 0.05, # I -> R xi = 0.02, # R -> S - vaccination_times = numeric(0), + vaccination_times = NULL, vaccination_coverage = rep(0.2, length(vaccination_times)), # N -> Y vaccine_efficacy = 1, ...) @@ -306,15 +308,52 @@ piecewise_run <- rbind(piecewise_run_initial$result, piecewise_run_final$result[ all.equal(uninterrupted_run, piecewise_run) ``` -We can try the same again, but this time set `restore_random_state = TRUE` to enable restoring the simulation state. This time we've successfully managed to reproduce the data from our uninterrupted run. +We can try to resume the simulation again, but this time set `restore_random_state = TRUE` to enable restoring the simulation state. This time we've successfully managed to reproduce the data from our uninterrupted run. ```{r} -set.seed(123) -piecewise_run_initial <- run_simulation(steps = 499) -piecewise_run_final <- run_simulation(steps = 1500, state = piecewise_run_initial$state, restore_random_state = TRUE) +piecewise_run_final <- run_simulation( + steps = 1500, + state = piecewise_run_initial$state, + restore_random_state = TRUE) + piecewise_run <- rbind(piecewise_run_initial$result, piecewise_run_final$result[500:1500,]) all.equal(uninterrupted_run, piecewise_run) ``` Using `restore_random_state = TRUE` resets the global random number generator's state, which could have surprising and undesirable side effects. It is generally useful in tests, but should be used carefully elsewhere. + +### Saving and restoring events {#events} + +By default, the `Event` and `TargetedEvent` classes save their schedule in the checkpoint and restore them when loading from a previous simulation state. +When restoring its schedule, the event will clear its existing schedule and overwrite any pending time steps with the schedule found in the saved state. + +For example in the code below, the simulation is first initialized with an event scheduled at `t=10` and is run for only 5 steps, hence the event does not yet trigger. On the second run, the event is seemingly scheduled for `t=15`. However, by resuming the simulation, the event's schedule is overwritten with the saved state, clearing the newly scheduled time. The event is triggered at `t=10` only. + +```{r} +e <- Event$new() +e$schedule(9) +state <- simulation_loop(timesteps = 5, events = list(e)) + +e <- Event$new() +e$schedule(14) +e$add_listener(function(t) cat("Triggered at timestep", t)) +simulation_loop(timesteps = 30, events = list(e), state = state) +``` + +Such a behaviour may prevent some modeling scenarios from being used with the save and restore effectively. For example, in our prior SIRS model, we wanted to be able to simulate historical data, where no intervention took place, and then resume it with different intervention parameters. If the vaccination event were to always restore its state from the original run, we would not be able to set a new schedule on the intervention runs. + +This behaviour can be modified on a per-event basis by passing `restore = FALSE` as an argument to the `Event` constructor. In that case, the schedule from the original run is ignored, and instead the schedule of the new run is only determined from its new initialization. + +In the modified example below, the event triggers at `t=15` as intended. When restoring the simulation, we had also schedules the event to trigger at `t=3`, but this never happens as it is before the point where the simulation is restored. + +```{r} +e <- Event$new() +e$schedule(9) +state <- simulation_loop(timesteps = 5, events = list(e)) + +e <- Event$new(restore=FALSE) +e$schedule(c(2, 14)) +e$add_listener(function(t, index) cat("Triggered at timestep", t)) +simulation_loop(timesteps = 30, events = list(e), state = state) +```