Correct community infectious and remove between-sector infections

DESCRIPTION

@@ -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 (>= 4.3)
+    R (>= 2.10)
 Imports: 
     checkmate,
     cli,

NEWS.md

@@ -1,3 +1,25 @@
+# daedalus 0.0.20
+
+This patch updates the state variable in `daedalus()` to substantially reduce the number of compartments; all empty compartments have been removed. This improves the speed of model runs as there are fewer derivatives to calculate.
+
+1. State variable is a 3D array with the dimensions: `c(49, 9, 3)`. The working-age economic sectors are added on as rows after the age-groups. Epidemiological compartments are unchanged. The third dimension holds data on vaccination status, including new vaccinations.
+
+2. Helper functions are updated to support the new state dimensions.
+
+3. $R_t$ calculation in `r_eff()` has been simplified to the proportional susceptible times $R_0$, pending a more accurate calculation.
+
+4. The internal helper function `prepare_output()` has been refactored to avoid use of `base::array2DF()` to remove the dependency on newer versions of R.
+
+5. Removes unused between-sector contacts in force-of-infection calculations.
+
+6. Moves some contact scaling to the parameter preparation stage to reduce operations during the model run.
+
+7. Correctly sums the total number of infectious individuals in the community; previously, the implementation resulted in segregation by economic sector.
+
+# daedalus 0.0.19
+
+This patch adds functionality to run `daedalus()` with the `country` argument passed as one of `country_codes_iso2c` or `country_codes_iso3c` for 2 and 3 letter ISO country codes.
+
 # daedalus 0.0.18
 
 This patch adds logging of daily new vaccinations and provides an output helper function, `get_new_vaccinations()`, to get daily new vaccinations.

R/class_country.R

@@ -430,12 +430,13 @@ prepare_parameters.daedalus_country <- function(x, ...) {
   singv <- max(Re(svd(cmcw)[["d"]]))
   cmcw <- cmcw / singv
 
+  workers <- get_data(x, "workers")
+
   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
+    contacts_workplace = cmw / workers,
+    contacts_consumer_worker = cmcw
   )
 }

R/class_infection.R

@@ -394,5 +394,15 @@ prepare_parameters.daedalus_infection <- function(x, ...) {
 
   validate_daedalus_infection(x)
   x <- unclass(x)
+
+  age_varying_params <- c("eta", "omega", "gamma_H")
+
+  x[age_varying_params] <- lapply(
+    x[age_varying_params], function(p) {
+      p <- c(p, rep(p[i_WORKING_AGE], N_ECON_SECTORS))
+      p
+    }
+  )
+
   x[names(x) != "name"]
 }

R/class_vaccination.R

@@ -282,12 +282,9 @@ scale_nu <- function(state, nu, uptake_limit) {
   # NOTE: state must be a 4D array with only vaccinated and unvaccinated layers
   # in dim 4
   total <- sum(
-    state[
-      , i_EPI_COMPARTMENTS, ,
-      c(i_VACCINATED_STRATUM, i_UNVACCINATED_STRATUM)
-    ]
+    state[, i_EPI_COMPARTMENTS, c(i_VACCINATED_STRATUM, i_UNVACCINATED_STRATUM)]
   )
-  total_vaccinated <- sum(state[, i_EPI_COMPARTMENTS, , i_VACCINATED_STRATUM])
+  total_vaccinated <- sum(state[, i_EPI_COMPARTMENTS, i_VACCINATED_STRATUM])
   prop_vaccinated <- total_vaccinated / total
 
   # NOTE: simplified scaling works only for uniform rates and start times

R/closure.R

@@ -8,10 +8,7 @@ make_response_threshold_event <- function(response_threshold) {
   # NOTE: input checking at top level
   ## event triggered when thresholds are crossed
   root_function <- function(time, state, parameters) {
-    state <- array(
-      state,
-      c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_STRATA)
-    )
+    state <- values_to_state(state)
 
     get_hospitalisations(state) - response_threshold
   }

R/constants.R

@@ -43,6 +43,10 @@
 #' @keywords model_constant
 N_AGE_GROUPS <- 4L
 
+#' @name model_constants
+#' @keywords model_constant
+i_AGE_GROUPS <- seq_len(N_AGE_GROUPS)
+
 #' @name model_constants
 #' @keywords model_constant
 N_VACCINE_STRATA <- 2L
@@ -63,6 +67,10 @@ i_WORKING_AGE <- 3L
 #' @keywords model_constant
 N_ECON_SECTORS <- 45L
 
+#' @name model_constants
+#' @keywords model_constant
+i_ECON_SECTORS <- seq(N_AGE_GROUPS + 1, N_AGE_GROUPS + N_ECON_SECTORS)
+
 #' @name model_constants
 #' @keywords model_constant
 i_EDUCATION_SECTOR <- 41L

R/costs.R

@@ -63,7 +63,7 @@ get_costs <- function(x, summarise_as = c("none", "total", "domain")) {
   )
 
   # absences due to infection, hospitalisation, death
-  model_data <- x$model_data
+  model_data <- get_data(x)
   worker_absences <- model_data[
     model_data$compartment %in%
       c("infect_symp", "infect_asymp", "hospitalised", "dead") &

R/daedalus.R

@@ -276,11 +276,12 @@ daedalus <- function(country,
   # 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(
+  rt <- r_eff(
     parameters[["r0"]],
     state_temp,
     parameters[["cm_unscaled"]]
-  ) >= 1.0
+  )
+  is_epidemic_growing <- rt >= 1.0
 
   if (!is_epidemic_growing) {
     rlang::env_bind(

R/model_helpers.R

@@ -1,13 +1,21 @@
 #' @title Generate a default initial state for DAEDALUS
-#' @description Function to prepare the model initial state.
+#' @description Function to prepare the model initial state. Assumes that
+#' 1 in every million individuals is initially infected, and that 60% are
+#' asymptomatic infections. This does not affect the actual probability of
+#' asymptomatic infections in the simulation, which is a property of a
+#' `<daedalus_infection>`.
 #'
 #' @inheritParams daedalus
 #'
-#' @return An array with as many dimensions as `N_ECON_STRATA` (currently, 46)
-#' with each layer giving the proportion of individuals of each group in each
-#' epidemiological compartment.
+#' @return An array with as many dimensions as `N_VACCINE_DATA_GROUPS`
+#' (currently, 3); rows specify the age and economic groups, columns specify the
+#' epidemiological compartments (including new infections and hospitalisations),
+#' and array layers hold information on vaccination status (including new
+#' vaccinations).
 #' @keywords internal
-make_initial_state <- function(country, initial_state_manual) {
+make_initial_state <- function(
+    country,
+    initial_state_manual) {
   # NOTE: country checked in daedalus()
   initial_infect_state <- list(
     p_infectious = 1e-6,
@@ -49,38 +57,32 @@ make_initial_state <- function(country, initial_state_manual) {
     dE = 0.0, dH = 0.0
   )
 
-  # build for all age groups
-  initial_state <- array(
-    rep(initial_state, each = N_AGE_GROUPS),
-    c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_DATA_GROUPS)
+  # build for all age groups and economic sectors (working age only)
+  initial_state <- matrix(
+    initial_state,
+    N_AGE_GROUPS + N_ECON_SECTORS,
+    N_MODEL_COMPARTMENTS,
+    byrow = TRUE
   )
-  # set vaccinated to zero
-  initial_state[, , , i_VACCINATED_STRATUM] <- 0.0
-  initial_state[, , , i_NEW_VAX_STRATUM] <- 0.0
 
-  # get demography and sector workforce
+  # get demography and sector workforce, including non-working
   demography <- get_data(country, "demography")
   sector_workforce <- get_data(country, "workers")
-
-  # calculate non-working working-age, and prepare initial state accounting
-  # for distribution of working age into economic sectors
   inactive_workers <- demography[i_WORKING_AGE] - sum(sector_workforce)
+  demography[i_WORKING_AGE] <- inactive_workers
 
-  initial_state[-i_WORKING_AGE, , i_NOT_WORKING, i_UNVACCINATED_STRATUM] <-
-    initial_state[-i_WORKING_AGE, , i_NOT_WORKING, i_UNVACCINATED_STRATUM] *
-      demography[-i_WORKING_AGE]
-
-  initial_state[i_WORKING_AGE, , i_NOT_WORKING, i_UNVACCINATED_STRATUM] <-
-    initial_state[i_WORKING_AGE, , i_NOT_WORKING, i_UNVACCINATED_STRATUM] *
-      inactive_workers
+  # multiply by demography and sector workforce, row-wise
+  initial_state[seq_len(N_AGE_GROUPS), ] <-
+    initial_state[seq_len(N_AGE_GROUPS), ] * demography
 
-  # explicit col-wise multiplication as R tries to guess interpretation of `*`
-  initial_state[i_WORKING_AGE, , -i_NOT_WORKING, i_UNVACCINATED_STRATUM] <-
-    initial_state[i_WORKING_AGE, , -i_NOT_WORKING, i_UNVACCINATED_STRATUM] %*%
-    diag(sector_workforce)
+  initial_state[(N_AGE_GROUPS + 1):nrow(initial_state), ] <-
+    initial_state[(N_AGE_GROUPS + 1):nrow(initial_state), ] * sector_workforce
 
-  # set all economic sector non-working age values to 0
-  initial_state[-i_WORKING_AGE, , -i_NOT_WORKING, ] <- 0.0
+  # add strata for vaccination groups and set to zero
+  initial_state <- array(
+    initial_state, c(dim(initial_state), N_VACCINE_DATA_GROUPS)
+  )
+  initial_state[, , -i_UNVACCINATED_STRATUM] <- 0
 
   initial_state
 }
@@ -145,8 +147,9 @@ get_closure_info <- function(mutables) {
 #' @return An array of dimensions (4, 9, 46, 3).
 #' @keywords internal
 values_to_state <- function(x) {
-  array(
-    x,
-    c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_DATA_GROUPS)
+  dim(x) <- c(
+    N_AGE_GROUPS + N_ECON_SECTORS, N_MODEL_COMPARTMENTS, N_VACCINE_DATA_GROUPS
   )
+
+  x
 }