Cpp model for IfG workshop

.Rbuildignore

@@ -29,3 +29,4 @@
 ^data-raw$
 ^doc$
 ^Meta$
+^.*_cpp.sh$

.gitignore

@@ -18,3 +18,4 @@ pkgdown
 scratch*
 /doc/
 /Meta/
+*_cpp.sh

.lintr

@@ -32,5 +32,6 @@ exclusions: list(
       missing_package_linter = Inf,
       namespace_linter = Inf
     ),
-    "R/cpp11.R"
+    "R/cpp11.R",
+    "R/RcppExports.R"
     )

DESCRIPTION

@@ -34,6 +34,8 @@ Imports:
     cli,
     data.table,
     deSolve,
+    Rcpp,
+    RcppEigen,
     rlang
 Suggests:
     countrycode,
@@ -46,6 +48,10 @@ Suggests:
     spelling,
     testthat (>= 3.0.0),
     tidyr
+LinkingTo: 
+    BH,
+    Rcpp,
+    RcppEigen
 VignetteBuilder: 
     knitr
 Config/Needs/website: jameel-institute/jameelinst.rpkg.theme

NAMESPACE

@@ -14,16 +14,20 @@ S3method(set_data,daedalus_vaccination)
 export(daedalus)
 export(daedalus_country)
 export(daedalus_infection)
+export(daedalus_rtm)
 export(daedalus_vaccination)
 export(get_costs)
 export(get_data)
 export(get_epidemic_summary)
 export(get_incidence)
+export(get_life_years_lost)
 export(get_new_vaccinations)
 export(is_daedalus_country)
 export(is_daedalus_infection)
 export(is_daedalus_vaccination)
 export(set_data)
+import(RcppEigen)
+importFrom(Rcpp,sourceCpp)
 importFrom(data.table,":=")
 importFrom(data.table,.BY)
 importFrom(data.table,.EACHI)
@@ -33,3 +37,4 @@ importFrom(data.table,.N)
 importFrom(data.table,.NGRP)
 importFrom(data.table,.SD)
 importFrom(data.table,data.table)
+useDynLib(daedalus, .registration = TRUE)

R/RcppExports.R

@@ -0,0 +1,7 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+.model_daedalus_cpp <- function(initial_state, params, contact_matrix, contacts_work, contacts_consumers, openness, hospital_capacity, t_start, t_end, auto_social_distancing = FALSE, social_distancing_mandate = 1.0, time_end = 100.0, increment = 1.0) {
+    .Call(`_daedalus_model_daedalus_internal`, initial_state, params, contact_matrix, contacts_work, contacts_consumers, openness, hospital_capacity, t_start, t_end, auto_social_distancing, social_distancing_mandate, time_end, increment)
+}
+

R/class_country.R

@@ -162,6 +162,7 @@ daedalus_country <- function(country,
   params$workers <- params$workers + 1
 
   # add value of statistical life (VSL)
+  life_expectancy <- daedalus::life_expectancy[[name]]
   vsl <- daedalus::life_value[[name]]
   gni <- daedalus::country_gni[[name]]
 
@@ -180,7 +181,8 @@ daedalus_country <- function(country,
       contacts_between_sectors =
         daedalus::economic_contacts[["contacts_between_sectors"]],
       vsl = vsl,
-      gni = gni
+      gni = gni,
+      life_expectancy = life_expectancy
     )
   )
   parameters <- Filter(parameters, f = function(x) !is.null(x))
@@ -216,7 +218,7 @@ validate_daedalus_country <- function(x) {
     "name", "demography", "contact_matrix",
     "contacts_workplace", "contacts_consumer_worker",
     "contacts_between_sectors", "workers", "gva",
-    "vsl", "hospital_capacity", "gni"
+    "vsl", "hospital_capacity", "gni", "life_expectancy"
   )
   has_invariants <- checkmate::test_names(
     attributes(x)$names,
@@ -308,6 +310,12 @@ validate_daedalus_country <- function(x) {
       checkmate::test_count(
         x$gni,
         positive = TRUE
+      ),
+    "Country `life_expectancy` must be a 4-element vector of positive values" =
+      checkmate::test_numeric(
+        x$life_expectancy,
+        len = N_AGE_GROUPS, any.missing = FALSE,
+        lower = 0, upper = 100 # reasonable upper limit
       )
   )
 

R/closure_data.R

@@ -20,3 +20,16 @@
 #' @source Multiple sources; to be updated shortly. See processing details in
 #' `data-raw/closure_data.R
 "closure_data"
+
+#' @title Economic sector names
+#'
+#' @description Names or descriptions of openness of economic sectors modelled
+#' in _daedalus_, to which closures apply.
+#'
+#' @format ## `econ_sector_names`
+#' A character vector of length `N_ECON_SECTORS` (45).
+#' @source Multiple sources; to be updated shortly. See processing details in
+#' `data-raw/closure_data.R
+#' @examples
+#' econ_sector_names
+"econ_sector_names"

R/costs.R

@@ -29,13 +29,20 @@
 #' (`education_cost_closures`), and pandemic-related absences due to illness and
 #' death (`education_cost_absences`).
 #'
-#' ## Life-years lost
+#' ## Life-value lost
 #'
 #' A four-element vector (for the number of age groups) giving the value of
 #' life-years lost per age group. This is calculated as the life-expectancy of
 #' each age group times the value of a statistical life, with all years assumed
 #' to have the same value.
 #'
+#' ## Life-years lost
+#'
+#' A four-element vector (for the number of age groups) giving the value of
+#' life-years lost per age group. This is calculated as the life-expectancy of
+#' each age group times the number of deaths in that age group. No quality
+#' adjustment is applied.
+#'
 #' @examples
 #' output <- daedalus("Canada", "influenza_1918")
 #'
@@ -61,10 +68,9 @@ get_costs <- function(x, summarise_as = c("none", "total", "domain")) {
   vsd <- vsd / 1e6 # for uniformity with daily GVA
   n_students <- x$country_parameters$demography[i_SCHOOL_AGE]
 
-  # cost of closures
-  economic_cost_closures <- sum(
-    gva * (1 - openness) * closure_duration
-  )
+  # cost of closures, per sector and total
+  sector_cost_closures <- gva * (1 - openness) * closure_duration
+  economic_cost_closures <- sum(sector_cost_closures)
 
   education_cost_closures <- sum(
     vsd * n_students * (1 - openness[i_EDUCATION_SECTOR]) *
@@ -85,8 +91,9 @@ get_costs <- function(x, summarise_as = c("none", "total", "domain")) {
   )
   workforce <- x$country_parameters$workers
 
-  # calculate daily GVA loss; scale GVA loss by openness when closures are
-  # active
+  # calculate daily GVA loss due to illness-related absencees;
+  # scale GVA loss by openness when closures are active
+  # assuming no working from home
   gva_loss <- worker_absences %*% diag(gva / workforce)
 
   gva_loss[seq(closure_start, closure_end), ] <-
@@ -113,22 +120,29 @@ get_costs <- function(x, summarise_as = c("none", "total", "domain")) {
     levels = unique(total_deaths$age_group)
   )
   total_deaths <- tapply(total_deaths$value, total_deaths$age_group, sum)
+  life_years_lost <- x$country_parameters$life_expectancy * total_deaths
 
   # NOTE: in million $s
-  life_years_lost <- x$country_parameters$vsl * total_deaths / 1e6
+  life_value_lost <- x$country_parameters$vsl * total_deaths / 1e6
 
   cost_list <- list(
     total_cost = NA,
     economic_costs = list(
       economic_cost_total = economic_cost_closures + economic_cost_absences,
       economic_cost_closures = economic_cost_closures,
-      economic_cost_absences = economic_cost_absences
+      economic_cost_absences = economic_cost_absences,
+      sector_cost_closures = sector_cost_closures,
+      sector_cost_absences = gva_loss
     ),
     education_costs = list(
       education_cost_total = education_cost_closures + education_cost_absences,
       education_cost_closures = education_cost_closures,
       education_cost_absences = education_cost_absences
     ),
+    life_value_lost = list(
+      life_value_lost_total = sum(life_value_lost),
+      life_value_lost_age = life_value_lost
+    ),
     life_years_lost = list(
       life_years_lost_total = sum(life_years_lost),
       life_years_lost_age = life_years_lost
@@ -138,7 +152,7 @@ get_costs <- function(x, summarise_as = c("none", "total", "domain")) {
   # probably a neater way of doing this
   cost_list$total_cost <- economic_cost_closures + economic_cost_absences +
     education_cost_closures + education_cost_absences +
-    sum(life_years_lost)
+    sum(life_value_lost)
 
   # return summary if requested, defaults to no summary
   summarise_as <- rlang::arg_match(summarise_as)
@@ -149,7 +163,7 @@ get_costs <- function(x, summarise_as = c("none", "total", "domain")) {
     domain = {
       cost_list[["total_cost"]] <- NULL
       vec_costs <- vapply(cost_list, `[[`, 1L, FUN.VALUE = numeric(1))
-      names(vec_costs) <- c("economic", "education", "life_years")
+      names(vec_costs) <- c("economic", "education", "life_value", "life_years")
 
       vec_costs
     }

R/daedalus-package.R

@@ -2,6 +2,7 @@
 "_PACKAGE"
 
 ## usethis namespace: start
+#' @import RcppEigen
 #' @importFrom data.table :=
 #' @importFrom data.table .BY
 #' @importFrom data.table .EACHI
@@ -11,5 +12,7 @@
 #' @importFrom data.table .NGRP
 #' @importFrom data.table .SD
 #' @importFrom data.table data.table
+#' @importFrom Rcpp sourceCpp
+#' @useDynLib daedalus, .registration = TRUE
 ## usethis namespace: end
 NULL