From 3e811ab948050df5010336bc481a2674efd34311 Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratikgupte16@gmail.com>
Date: Mon, 7 Oct 2024 11:18:26 +0100
Subject: [PATCH 1/7] Prevent closure activation if epidemic is not growing

---
 R/closure.R | 29 +++++++++++++++++++++++------
 1 file changed, 23 insertions(+), 6 deletions(-)

diff --git a/R/closure.R b/R/closure.R
index 81fd06f..1d2ee10 100644
--- a/R/closure.R
+++ b/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)
   }

From 69778eaf196bec6ba1aabfffad923341484fcdc2 Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratikgupte16@gmail.com>
Date: Mon, 7 Oct 2024 11:18:47 +0100
Subject: [PATCH 2/7] Manually turn off closures if epidemic is not growing

---
 R/daedalus.R | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)

diff --git a/R/daedalus.R b/R/daedalus.R
index 0abfc44..faeeb66 100644
--- a/R/daedalus.R
+++ b/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
 

From 71ceda18e3dc7496c5ed44fd8bd11f42047b0eae Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratikgupte16@gmail.com>
Date: Mon, 7 Oct 2024 11:19:24 +0100
Subject: [PATCH 3/7] Test that closures are not activated in 'herd immunity'
 scenarios

---
 tests/testthat/test-closures.R | 39 ++++++++++++++++++++++++++++++++++
 1 file changed, 39 insertions(+)

diff --git a/tests/testthat/test-closures.R b/tests/testthat/test-closures.R
index ee55f6a..e35cfb6 100644
--- a/tests/testthat/test-closures.R
+++ b/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
+        )
+      }
+    )
+  )
 })

From 6d8ed9cd68c3efbffbd85ac3970831450cbfcf85 Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratikgupte16@gmail.com>
Date: Mon, 7 Oct 2024 11:19:36 +0100
Subject: [PATCH 4/7] Bump to v0.0.15

---
 DESCRIPTION |  4 ++--
 NEWS.md     | 14 ++++++++++++++
 2 files changed, 16 insertions(+), 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 5cb5f2c..101610f 100644
--- a/DESCRIPTION
+++ b/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,
diff --git a/NEWS.md b/NEWS.md
index a1c58c5..a70f83c 100644
--- a/NEWS.md
+++ b/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".

From 80a785a18de0874ac63448eeafa0cccff8331798 Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratikgupte16@gmail.com>
Date: Tue, 8 Oct 2024 10:48:39 +0100
Subject: [PATCH 5/7] Add unscaled contacts to parameter list

---
 NEWS.md                        | 4 +++-
 R/class_country.R              | 4 +++-
 R/closure.R                    | 9 +++------
 R/daedalus.R                   | 2 +-
 R/pkg_generics.R               | 3 +++
 man/prepare_parameters.Rd      | 2 ++
 tests/testthat/test-daedalus.R | 1 +
 7 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index a70f83c..a61b4b3 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -10,7 +10,9 @@ 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()`.
+- The package now requires a minimum R version >= v4.3 due to the use of the new function `array2DF()`;
+
+- `prepare_parameters.daedalus_country()` now provides a contact matrix scaled by the leading eigenvalue but not the demography for use in `r_eff()`.
 
 # daedalus 0.0.14
 
diff --git a/R/class_country.R b/R/class_country.R
index 2787e4c..5b53fb1 100644
--- a/R/class_country.R
+++ b/R/class_country.R
@@ -418,7 +418,8 @@ prepare_parameters.daedalus_country <- function(x, ...) {
   # scale contacts by largest real eigenvalue
   cm <- get_data(x, "contact_matrix")
   eigv <- max(Re(eigen(cm)$values))
-  cm <- (cm / eigv) %*% diag(1 / demography)
+  cm_unscaled <- cm / eigv
+  cm <- (cm_unscaled) %*% diag(1 / demography)
 
   cmw <- get_data(x, "contacts_workplace")
   cmw <- cmw / max(cmw) # max(cmw) is leading eigenvalue of diag matrix cmw
@@ -431,6 +432,7 @@ prepare_parameters.daedalus_country <- function(x, ...) {
   list(
     demography = demography,
     contact_matrix = cm,
+    cm_unscaled = cm_unscaled, # for use in Rt calculations
     contacts_workplace = cmw,
     contacts_consumer_worker = cmcw,
     contacts_between_sectors = get_data(x, "contacts_between_sectors") # 0s
diff --git a/R/closure.R b/R/closure.R
index 1d2ee10..1a0f4ac 100644
--- a/R/closure.R
+++ b/R/closure.R
@@ -23,8 +23,8 @@ make_response_threshold_event <- function(response_threshold) {
       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)
+
+    rt <- r_eff(parameters[["r0"]], state, parameters[["cm_unscaled"]])
 
     # NOTE: to ensure only first hosp threshold crossing is logged
     is_hosp_switch_on <- rlang::env_get(parameters[["mutables"]], "hosp_switch")
@@ -63,11 +63,8 @@ make_rt_end_event <- function() {
       c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_STRATA)
     )
 
-    # because contacts are divided by demography during parameter prep
-    cm <- parameters[["contact_matrix"]] %*% diag(parameters[["demography"]])
-
     # arbitrary precision, may not be hit!
-    r_eff(parameters[["r0"]], state, cm) - 0.99
+    rt <- r_eff(parameters[["r0"]], state, parameters[["cm_unscaled"]]) - 0.99
   }
 
   event_function <- function(time, state, parameters) {
diff --git a/R/daedalus.R b/R/daedalus.R
index faeeb66..b27e701 100644
--- a/R/daedalus.R
+++ b/R/daedalus.R
@@ -286,7 +286,7 @@ daedalus <- function(country,
   is_epidemic_growing <- r_eff(
     parameters[["r0"]],
     state_temp,
-    parameters[["contact_matrix"]] %*% diag(parameters[["demography"]])
+    parameters[["cm_unscaled"]]
   ) >= 1.0
 
   if (!is_epidemic_growing) {
diff --git a/R/pkg_generics.R b/R/pkg_generics.R
index 228e13b..f8ff92a 100644
--- a/R/pkg_generics.R
+++ b/R/pkg_generics.R
@@ -102,6 +102,9 @@ set_data <- function(x, ...) {
 #' - `contact_matrix`: the scaled contact matrix suitable for multiplication
 #' with \eqn{R_0} in force of infection calculations;
 #'
+#' - `cm_unscaled`: the contact matrix multiplied column-wise by the demography
+#' vector, suitable for calculations of \eqn{R_{\text{eff}}}; see [r_eff()].
+#'
 #' - `contacts_workplace`: the contacts in workplaces scaled by their largest
 #' value (which is the leading eigenvalue of the diagonal matrix of contacts);
 #'
diff --git a/man/prepare_parameters.Rd b/man/prepare_parameters.Rd
index 596b81d..c8e4892 100644
--- a/man/prepare_parameters.Rd
+++ b/man/prepare_parameters.Rd
@@ -45,6 +45,8 @@ The returned parameter list consists of:
 \item \code{demography}: the demography vector;
 \item \code{contact_matrix}: the scaled contact matrix suitable for multiplication
 with \eqn{R_0} in force of infection calculations;
+\item \code{cm_unscaled}: the contact matrix multiplied column-wise by the demography
+vector, suitable for calculations of \eqn{R_{\text{eff}}}; see \code{\link[=r_eff]{r_eff()}}.
 \item \code{contacts_workplace}: the contacts in workplaces scaled by their largest
 value (which is the leading eigenvalue of the diagonal matrix of contacts);
 \item \code{contacts_consumer_worker}: contacts in workplaces distributed in
diff --git a/tests/testthat/test-daedalus.R b/tests/testthat/test-daedalus.R
index a1d4409..3cae26c 100644
--- a/tests/testthat/test-daedalus.R
+++ b/tests/testthat/test-daedalus.R
@@ -65,6 +65,7 @@ test_that("daedalus: basic expectations", {
 })
 
 # test that daedalus runs for all epidemic infection parameter sets
+skip()
 test_that("daedalus: Runs for all country x infection x response", {
   country_infection_combos <- data.table::CJ(
     country = daedalus::country_names,

From b860eb099d792c10478bc474079bab5413cbf2ab Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratikgupte16@gmail.com>
Date: Tue, 8 Oct 2024 10:52:48 +0100
Subject: [PATCH 6/7] Skip no tests

---
 tests/testthat/test-daedalus.R | 1 -
 1 file changed, 1 deletion(-)

diff --git a/tests/testthat/test-daedalus.R b/tests/testthat/test-daedalus.R
index 3cae26c..a1d4409 100644
--- a/tests/testthat/test-daedalus.R
+++ b/tests/testthat/test-daedalus.R
@@ -65,7 +65,6 @@ test_that("daedalus: basic expectations", {
 })
 
 # test that daedalus runs for all epidemic infection parameter sets
-skip()
 test_that("daedalus: Runs for all country x infection x response", {
   country_infection_combos <- data.table::CJ(
     country = daedalus::country_names,

From c653b1c8e0402ec0ba5f7c5331208b25d003c6b0 Mon Sep 17 00:00:00 2001
From: Pratik Gupte <pratikgupte16@gmail.com>
Date: Tue, 8 Oct 2024 11:03:34 +0100
Subject: [PATCH 7/7] Fix lintr errors

---
 R/closure.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/closure.R b/R/closure.R
index 1a0f4ac..5736fb7 100644
--- a/R/closure.R
+++ b/R/closure.R
@@ -64,7 +64,7 @@ make_rt_end_event <- function() {
     )
 
     # arbitrary precision, may not be hit!
-    rt <- r_eff(parameters[["r0"]], state, parameters[["cm_unscaled"]]) - 0.99
+    r_eff(parameters[["r0"]], state, parameters[["cm_unscaled"]]) - 0.99
   }
 
   event_function <- function(time, state, parameters) {