detailed PAK scenarios

setup/PAK_detail/Makefile

@@ -0,0 +1,41 @@
+
+-include ../../local.makefile
+
+COVIDMPATH ?= ../../../covidm
+DATADIR ?= ~/Dropbox/covidm_reports
+#DATADIR ?= ~/Dropbox/covidm_test
+SAMPN ?= 500
+INDIR := ${DATADIR}/generation
+OUTDIR := ${DATADIR}/hpc_inputs
+DETDIR := ${DATADIR}/hpc_detailed
+
+R = Rscript $^ $@
+
+# TODO other directory structure dependencies
+${DETDIR}:
+	mkdir -p $@
+
+params: create_params_set.R ${DETDIR}/worldpop5yr.lfs.csv
+	Rscript $^ ${COVIDMPATH} ${DETDIR}
+
+# currently, no explicit province contact matrices, so use national one
+${DETDIR}/contact_matrices.rds: ${OUTDIR}/PAK/contact_matrices.rds
+	cp $< $@
+
+# TODO update dependencies
+# ${DETDIR}/%/params_set.rds: create_params_set.R ${INDIR}/data_contacts_missing.csv
+#	mkdir -p $(@D)
+#	Rscript $^ ${COVIDMPATH} $* $@
+
+${DETDIR}/timing.rds: ${OUTDIR}/PAK/timing.rds
+	cp $< $@
+
+${DETDIR}/run_options.rds: ${OUTDIR}/run_options.rds
+	cp $< $@
+
+${DETDIR}/alt_scenarios.rds: generate_scenario_parameters.R
+	Rscript $^ $@
+
+${DETDIR}/report_ref.rds: ${DETDIR}/alt_scenarios.rds
+
+dosetup: params ${DETDIR}/alt_scenarios.rds

setup/PAK_detail/create_params_set.R

@@ -0,0 +1,186 @@
+#' default simulation parameters
+suppressPackageStartupMessages({
+  require(data.table)
+})
+
+.args <- if (interactive()) c(
+  "~/Dropbox/covidm_reports/hpc_detailed/worldpop5yr.lfs.csv",
+  "../covidm",
+  "~/Dropbox/covidm_reports/hpc_detailed" # going to be writing admin#/param_set.rds
+) else commandArgs(trailingOnly = TRUE)
+#' @examples 
+#' .args <- gsub("ZWE","guineabissau",.args)
+#' .args <- gsub("ZWE","palestine",.args)
+
+reference = fread(
+  .args[1], strip.white = FALSE)[
+    level == 1 & country == "PAK"
+]
+country <- matref <- "Pakistan"
+
+cm_path = .args[2]
+outfilefmt = file.path(tail(.args, 1),"%s","params_set.rds")
+
+cm_force_rebuild = F;
+cm_build_verbose = F;
+cm_force_shared = T;
+cm_version = 2;
+
+suppressPackageStartupMessages({
+  source(file.path(cm_path, "R", "covidm.R"))
+})
+
+# country <- cm_populations[
+#   country_code == popcode,
+#   unique(as.character(name))
+# ]
+
+#set up node-runs
+probs = fread(
+  "Age_low,Age_high,Prop_symptomatic,IFR,Prop_inf_hosp,Prop_inf_critical,Prop_critical_fatal,Prop_noncritical_fatal,Prop_symp_hospitalised,Prop_hospitalised_critical
+0,9,0.66,8.59E-05,0.002361009,6.44E-05,0.5,0,0,0.3
+10,19,0.66,0.000122561,0.003370421,9.19E-05,0.5,9.47E-04,0.007615301,0.3
+20,29,0.66,0.000382331,0.010514103,0.000286748,0.5,0.001005803,0.008086654,0.3
+30,39,0.66,0.000851765,0.023423527,0.000638823,0.5,0.001231579,0.009901895,0.3
+40,49,0.66,0.001489873,0.0394717,0.001117404,0.5,0.002305449,0.018535807,0.3
+50,59,0.66,0.006933589,0.098113786,0.005200192,0.5,0.006754596,0.054306954,0.3
+60,69,0.66,0.022120421,0.224965092,0.016590316,0.5,0.018720727,0.150514645,0.3
+70,79,0.66,0.059223786,0.362002579,0.04441784,0.5,0.041408882,0.332927412,0.3
+80,100,0.66,0.087585558,0.437927788,0.065689168,0.5,0.076818182,0.617618182,0.3"
+)
+
+#increase CFR
+cfr_RR <- 1.5
+probs[, Prop_critical_fatal := Prop_critical_fatal * cfr_RR]
+probs[, Prop_noncritical_fatal := Prop_noncritical_fatal * cfr_RR]
+
+#min(1, x) does not work for noncritical_fatal, for some reason
+probs[Prop_critical_fatal > 1, Prop_critical_fatal := 1]
+probs[Prop_noncritical_fatal > 1, Prop_noncritical_fatal := 1]
+
+reformat = function(P, lmic_adjust=TRUE) {
+  # no info to re-weight these, so assume 70-74 is like 70-79, and 75+ is like 80+
+  if(lmic_adjust){
+    P <- P[2:length(P)]
+    return (rep(P[1:8], each = 2))
+  } else {
+    return (c(rep(P[1:7], each = 2),P[8:9])) 
+  }
+}
+
+P.icu_symp     = reformat(probs[, Prop_symp_hospitalised * Prop_hospitalised_critical], lmic_adjust = TRUE);
+P.nonicu_symp  = reformat(probs[, Prop_symp_hospitalised * (1 - Prop_hospitalised_critical)], lmic_adjust = TRUE);
+P.death    = reformat(probs[, Prop_noncritical_fatal], lmic_adjust = TRUE);
+
+max_time <- 60
+tres <- 0.25
+
+ponset2hosp <- cm_delay_gamma(7, 7, max_time, tres)$p
+pignore <- cm_delay_skip(max_time, tres)$p
+icustay <- cm_delay_gamma(10, 10, max_time, tres)$p
+nonicustay <- cm_delay_gamma(8, 8, max_time, tres)$p
+ponset2death <- cm_delay_gamma(22, 22, max_time, tres)$p
+
+cm_multinom_process <- function(
+  src, outcomes, delays,
+  report = ""
+) {
+  if ("null" %in% names(outcomes)) {
+    if (length(report) != length(outcomes)) report <- rep(report, length(outcomes))
+    report[which(names(outcomes)=="null")] <- ""
+    if (!("null" %in% names(delays))) {
+      delays$null <- c(1, rep(0, length(delays[[1]])-1))
+    }
+  } else if (!all(rowSums(outcomes)==1)) {
+    report <- c(rep(report, length(outcomes)), "")
+    outcomes$null <- 1-rowSums(outcomes)
+    delays$null <- c(1, rep(0, length(delays[[1]])-1))
+  }
+  nrow <- length(outcomes)
+  list(
+    source = src, type="multinomial", names=names(outcomes), report = report,
+    prob = t(as.matrix(outcomes)), delays = t(as.matrix(delays))
+  )
+}
+
+cm_track_process <- function(src, name, delays, agecats = 16, report = "p") {
+  list(
+    source = src, type="multinomial", names = name, report = report,
+    prob = matrix(1, nrow = 1, ncol = agecats),
+    delays = t(delays)
+  )
+}
+
+burden_processes = list(
+  # process of sending symptomatic cases to hospital icu or ward
+  cm_multinom_process(
+    "Ip",
+    outcomes = data.frame(to_icu = P.icu_symp, to_nonicu = P.nonicu_symp),
+    delays = data.frame(to_icu = ponset2hosp, to_nonicu = ponset2hosp)
+  ),
+  # track icu prevalance
+  cm_track_process("to_icu", "icu", icustay),
+  # track ward prevalence
+  cm_track_process("to_nonicu", "nonicu", nonicustay),
+  # track infections - get from delta R prevalence
+#  cm_track_process("S", "infection", pignore, report="i"),
+  # send some cases to death, tracking outcidence
+  cm_multinom_process(
+    "Ip",
+    outcomes = data.table(death=P.death),
+    delays = data.table(death=ponset2death),
+    report = "o"
+  )
+)
+
+popnorm <- function(x, seed_cases = 50){
+  #no contact matrix - will be overwritten by empirical ones
+  x$matrices$home <- x$matrices$work <- x$matrices$school <- x$matrices$other <- x$matrices$other*0
+
+  #age-specific probability of being symptomatic
+  #x$y <- c(rep(0.056, 3), rep(0.49, 8), rep(0.74, 8))
+  #new values proposed by nick
+  x$y <- c(
+    rep(0.2973718, 2), rep(0.2230287, 2), rep(0.4191036, 2),
+    rep(0.4445867, 2), rep(0.5635720, 2), rep(0.8169443, 6)
+  )
+
+  #no cases in empty compartments
+  x$dist_seed_ages <- as.numeric(!(x$size == 0))
+
+  #seed cases
+  x$seed_times <- rep(0, seed_cases)
+
+  return(x)
+}
+
+for (admin_code in reference[, key]) {
+  params_set <- list()
+
+  params1 <- cm_parameters_SEI3R(
+    country, matref,
+    deterministic=FALSE, date_end = 365*2
+  )
+
+  extractpop <- melt(reference[
+    key == admin_code, .SD, .SDcols = grep("(f|m)_", colnames(reference), value = TRUE)
+  ], measure.vars = grep("(f|m)_", colnames(reference), value = TRUE))
+  extractpop[, agelb := as.integer(gsub(".+_(\\d+)","\\1", variable)) ]
+  redage <- extractpop[,.(value = sum(value)), keyby=agelb ]
+  redage[16]$value <- sum(redage[16:.N]$value)
+
+  params1$pop[[1]]$size <- redage[1:16, value]
+
+  params1$processes = burden_processes
+
+  params1$pop <- lapply(params1$pop, popnorm)
+  #params1$time1 <- as.Date(params1$time1)
+
+  params_set[[1]] <- params1
+
+  dir.create(dirname(sprintf(outfilefmt, admin_code)))
+  saveRDS(params_set, sprintf(outfilefmt, admin_code))
+
+}
+
+#general population parameters

setup/PAK_detail/generate_scenario_parameters.R

@@ -0,0 +1,103 @@
+suppressPackageStartupMessages({
+  require(data.table)
+})
+
+.args <- if (interactive()) c(
+  "~/Dropbox/covidm_reports/hpc_detailed/alt_scenarios.rds"
+) else commandArgs(trailingOnly = TRUE)
+
+outfile <- tail(.args, 1)
+reffile <- gsub("/(\\w+)\\.rds", "/report_ref.rds", outfile)
+
+# say that 11 March was the 50 active infections date
+mar23 <- as.integer(as.Date("2020-03-23") - as.Date("2020-03-11"))
+
+dg <- function(...) data.table(expand.grid(..., stringsAsFactors = FALSE))
+
+ref <- function(
+  age_split = NA,
+  self_iso = 0,
+  population = -1, coverage = 1, # coverage generally only applies to multi-pop models
+  school = 0, home = 0, work = 0, other = 0, # recall == these are reductions
+  travel = 0,
+  start_day = NA_integer_, end_day = Inf,
+  trigger_type = NA_character_,
+  trigger_value = NA_real_,
+  ...
+) do.call(dg, c(as.list(environment()), list(...)))
+
+#' TODO: ignore? going to just iterate over rows where scen_id = 1, and if it's empty, nothing to do
+unmitigated <- ref(scen_type = "unmitigated")[, scen_id := 1 ]
+
+reportref <- data.table(scen_id=integer(), label=character())
+append_scenarios <- function(dt, new_scen_ids, new_labels, scen_type) {
+  newdt <- data.table(scen_id = new_scen_ids, label = new_labels)
+  rbind(dt, newdt)
+}
+
+scen_counter <- 1
+tagscenario <- function(dt, sc) { 
+  dt[, scen_id := (1:.N) + sc ]
+  return(dt[.N, scen_id])
+}
+
+# cat(sprintf("unmitigated: %i", 1))
+reportref <- append_scenarios(reportref, 1L, "Unmitigated")
+
+#' all remaining scenarios consider some level of self-isolation with symptoms
+#' and potentially different school effects
+self_iso_lvls <- c(0.25)
+schoolslowhi <- c(low=0.8, hi=1) # on = 0 reduction, off = 100% reduction
+
+#'  - generic social distancing in work / other, levels 20-80%, by 10%
+#' social distancing levels; use minimum for combination interventions
+
+work_dist <- c(0.2, 0.4)
+home_dist <- c(0.2, 0.4)
+other_dist <- c(0, 0.2)
+
+low_soc_all <- ref(
+  self_iso = self_iso_lvls[1],
+  home = home_dist,
+  school = schoolslowhi,
+  other = other_dist,
+  work = work_dist,
+  trigger_type = "day",
+  trigger_value = mar23
+)[, start_day := 0 ][, end_day := 90 ]
+scen_counter <- tagscenario(low_soc_all, scen_counter)
+
+# cat(sprintf("40%% elder sheilding at home + low level all other measures: %i", elder_shielding[home == 0.4, scen_id]))
+reportref <- append_scenarios(
+  reportref,
+  low_soc_all[, scen_id],
+  sprintf("generic options")
+)
+
+#'  - schools shutdown + continuous low level social distance (20%) + weekly alternating extra social distancing (40-50-60-70) => that + 20%
+alternating <- ref(
+  self_iso = self_iso_lvls[1],
+  school = schoolslowhi["hi"],
+  home = max(home_dist),
+  work = max(work_dist),
+  other = max(other_dist),
+  start_day = mar23, end_day = mar23+30,
+  trigger_type = "stride",
+  trigger_value = 90
+)
+scen_counter <- tagscenario(alternating, scen_counter)
+
+# cat("no school + intermittent 40%% lockdown: %i", alternating[work == 0.5, scen_id])
+reportref <- append_scenarios(
+  reportref,
+  alternating[, scen_id],
+  sprintf("alternating")
+)
+
+all_scenarios <- setkey(rbind(
+  low_soc_all,
+  alternating
+), scen_id)
+
+saveRDS(reportref, reffile)
+saveRDS(all_scenarios, outfile)

simulate/run_scenarios_detail.R

@@ -4,17 +4,17 @@ suppressPackageStartupMessages({
 })
 
 .args <- if (interactive()) c(
-  "../covidm", "ZWE", "001", 
-  sprintf("~/Dropbox/covidm_reports/hpc_inputs"),
-  "simulate/ZWE/001.qs"
+  "../covidm", "2428", "001", 
+  sprintf("~/Dropbox/covidm_reports/hpc_detailed"),
+  "simulate/detail/001.qs"
 ) else commandArgs(trailingOnly = TRUE)
 
 # load covidm
 cm_path = .args[1]
 cm_force_rebuild = F;
 cm_build_verbose = F;
 cm_force_shared = T
-cm_version = 1
+cm_version = 2
 
 suppressPackageStartupMessages({
   source(paste0(cm_path, "/R/covidm.R"))
@@ -102,17 +102,6 @@ uref <- unname(as.matrix(yu[, .SD, .SDcols = grep("u_",colnames(yu))]))
 
 has_age_split <- iv_data[, any(!is.na(age_split))]
 
-reduce_ages <- function (dt) {
-  fctr <- function(i, lvls = c("<14", "15-29", "30-44","45-59", "60+")) factor(
-    lvls[i], levels = lvls, ordered = T
-  )
-  dt[between(as.integer(group), 1, 3), age := fctr(1) ]
-  dt[between(as.integer(group), 4, 6), age := fctr(2) ]
-  dt[between(as.integer(group), 7, 9), age := fctr(3) ]
-  dt[between(as.integer(group), 10, 12), age := fctr(4) ]
-  dt[as.integer(group) >= 13, age := fctr(5) ]
-}
-
 cm_calc_R0_extended <- function(
   params
 ){
@@ -354,13 +343,11 @@ for(i in 1:nrow(run_options)){
     model_seed = run_options[i, model_seed]
   )$dynamics
 
-  result <- reduce_ages(
-    sim[
+  result <- sim[
       compartment %in% c("cases","death_o","icu_p","nonicu_p","R")
-    ]
-  )[,
+    ][,
     .(value = sum(value)),
-    keyby = .(run, t, age, compartment)
+    keyby = .(run, t, group, compartment)
   ][value != 0][, run := i ]
 
   rm(params)