Fix logging of closure limits

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: daedalus
 Title: Model health, social, and economic costs of a pandemic using
     _DAEDALUS_
-Version: 0.0.14
+Version: 0.0.15
 Authors@R: c(
     person("Pratik", "Gupte", , "p.gupte24@imperial.ac.uk", role = c("aut", "cre"),
            comment = c(ORCID = "0000-0001-5294-7819")),
@@ -16,7 +16,7 @@ License: MIT + file LICENSE
 URL: https://github.com/j-idea/daedalus
 BugReports: https://github.com/j-idea/daedalus/issues
 Depends: 
-    R (>= 2.10)
+    R (>= 4.3)
 Imports: 
     checkmate,
     cli,

NEWS.md

@@ -1,3 +1,17 @@
+# daedalus 0.0.15
+
+This patch fixes an issue where closures were not ended in model runs where the closure began after the epidemic had stopped growing. This mostly affected edge cases of countries with large spare hospital capacity, and relatively late `response_time`s. In such cases the closure end time was assigned as the simulation end time, inflating costs related to closures. The fix:
+
+- Prevents closures from being activated by the `hospital_capacity` trigger if the epidemic is not growing, even if the hospital capacity threshold is crossed;
+
+- Prevents closures from being activated by the `response_time` trigger if the epidemic is not growing between stage 01 and stage 02. Closures are manually turned off if the epidemic is not growing ($R_t < 1.0$).
+
+Tests for different response times check that the model behaves as expected.
+
+## Miscellaneous changes
+
+- The package now requires a minimum R version >= v4.3 due to the use of the new function `array2DF()`.
+
 # daedalus 0.0.14
 
 This patch adds vaccine investment scenarios to _daedalus_. All model runs must now include an assumption about the level of advance vaccine investment, defaulting to "none".

R/closure.R

@@ -17,15 +17,32 @@ make_response_threshold_event <- function(response_threshold) {
   }
 
   event_function <- function(time, state, parameters) {
-    # prevent flipping switch when checkEventFunc runs
-    # turn closure and excess hospitalisation switch on
+    # NOTE: prevent flipping switch when checkEventFunc runs
+    # NOTE: prevent response activation if epdedemic is not growing
+    state <- array(
+      state,
+      c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_STRATA)
+    )
+    cm <- parameters[["contact_matrix"]] %*% diag(parameters[["demography"]])
+    rt <- r_eff(parameters[["r0"]], state, cm)
+
+    # NOTE: to ensure only first hosp threshold crossing is logged
     is_hosp_switch_on <- rlang::env_get(parameters[["mutables"]], "hosp_switch")
+
     if (time != parameters[["min_time"]] && !is_hosp_switch_on) {
-      rlang::env_bind(
-        parameters[["mutables"]],
-        switch = TRUE, hosp_switch = TRUE,
-        closure_time_start = time
+      rlang::env_poke(
+        parameters[["mutables"]], "hosp_switch", TRUE
       )
+
+      # NOTE: trigger response and log closure start time only if
+      # epidemic is growing
+      if (rt >= 1.0) {
+        rlang::env_bind(
+          parameters[["mutables"]],
+          switch = TRUE,
+          closure_time_start = time
+        )
+      }
     }
     as.numeric(state)
   }

R/daedalus.R

@@ -275,6 +275,28 @@ daedalus <- function(country,
   ]
   initial_state <- as.numeric(initial_state)
 
+  ### cancel closures if epidemic is not growing ###
+  state_temp <- array(
+    initial_state,
+    c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_STRATA)
+  )
+  # close intervention IFF epidemic is not growing
+  # NOTE: this step placed here as a conditional on r_eff < 1.0 is not
+  # practical in a root-finding function (Error: 'root too near initial point')
+  is_epidemic_growing <- r_eff(
+    parameters[["r0"]],
+    state_temp,
+    parameters[["contact_matrix"]] %*% diag(parameters[["demography"]])
+  ) >= 1.0
+
+  if (!is_epidemic_growing) {
+    rlang::env_bind(
+      parameters[["mutables"]],
+      switch = FALSE,
+      closure_time_end = response_time
+    )
+  }
+
   # reset min time
   parameters[["min_time"]] <- response_time
 

tests/testthat/test-closures.R

@@ -201,4 +201,43 @@ test_that("Closures: correct logging of time limits", {
     response_data[["closure_info"]][["closure_time_end"]],
     time_end
   )
+
+  # Check that closures are terminated at the start time when
+  # the epidemic is not growing
+  dvx <- daedalus_vaccination("low", vax_start_time = 90)
+  tend <- 100
+  response_time <- seq(10, 80, 10)
+
+  # run scenario with very few susceptibles
+  invisible(
+    lapply(
+      response_time, function(x) {
+        data <- daedalus(
+          "United States", "influenza_1957",
+          initial_state_manual = list(
+            p_infectious = 0.9
+          ),
+          response_strategy = "elimination",
+          response_time = x,
+          vaccine_investment = dvx,
+          time_end = tend
+        )
+
+        closure_info <- data$response_data$closure_info
+
+        expect_identical(
+          closure_info$closure_time_end,
+          closure_info$closure_time_start
+        )
+        expect_identical(
+          closure_info$closure_time_end,
+          x
+        )
+        expect_identical(
+          closure_info$closure_duration,
+          0.0
+        )
+      }
+    )
+  )
 })