diff --git a/Makefile b/Makefile index 75386db13..82bf4558f 100644 --- a/Makefile +++ b/Makefile @@ -4,11 +4,11 @@ WD := ../shared-svn/projects/episim/matsim-files # All available Scenarios -ALL := BerlinWeek MunichWeek +ALL := BerlinBrandenburg JAR := matsim-episim-*.jar # Shortcut to the scenario creation tool -sc = java -Xmx20G -cp $(JAR) org.matsim.run.ScenarioCreation +sc = java -Xmx30G -cp $(JAR) org.matsim.run.ScenarioCreation .PHONY: all clean $(ALL) @@ -24,7 +24,7 @@ clean: # Includes all the scenarios with local variables # https://stackoverflow.com/questions/32904790/can-i-have-local-variables-in-included-makefiles -SUBDIRS := scenarios/Cologne.mk +SUBDIRS := scenarios/BerlinBrandenburg.mk define INCLUDE_FILE path = $S include $S diff --git a/README.md b/README.md index bb5cac238..399c1f142 100644 --- a/README.md +++ b/README.md @@ -56,3 +56,4 @@ The **MATSim input files, output files, analysis data and visualizations** are l For more information about the methodology and preliminary results, see VSP working paper https://dx.doi.org/10.14279/depositonce-9835 . For more information about MATSim, see here: https://www.matsim.org/. + diff --git a/scenarios/BerlinBrandenburg.mk b/scenarios/BerlinBrandenburg.mk new file mode 100644 index 000000000..ee49af71d --- /dev/null +++ b/scenarios/BerlinBrandenburg.mk @@ -0,0 +1,160 @@ +# To run this, you need a newer version of make. 4.4.1 works (in macOS under the command of gmake) +# To run, execute "gmake BerlinBrandenburg" in terminal +in = $(WD)/snz/BerlinBrandenburg/original-data +out = $(WD)/snz/BerlinBrandenburg/episim-input +tmp = $(WD)/snz/BerlinBrandenburg/processed-data + +BerlinBrandenburg: $(JAR) $(out)/bb_2020-week_snz_episim_events_wt_25pt_split.xml.gz $(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_25pt_split.xml.gz \ +$(out)/bb_2020-week_snz_episim_events_wt_100pt_split.xml.gz $(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt_split.xml.gz + echo "Building Berlin Brandenburg Week scenario" + +BerlinBrandenburgSamples: $(JAR) $(out)/samples/bb_2020-week_snz_episim_events_wt_10pt_split.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_5pt_split.xml.gz \ +$(out)/samples/bb_2020-week_snz_episim_events_wt_1pt_split.xml.gz + echo "Building Berlin Brandenburg Week samples" + +$(tmp)/personIds.diluted.txt.gz: + $(sc) filterPersons $(in)/de2020gsmwt_events_reduced.xml.gz\ + --facilities $(in)/facilities_assigned_simplified.xml.gz\ + --attributes $(in)/populationAttributes.xml.gz\ + --shape-file $(out)/shape-File/BerlinBrandenburg.shp\ + --shape-crs EPSG:25833\ + --output $@ + + + +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_100pt.xml.gz: $(tmp)/personIds.diluted.txt.gz + $(sc) convertPersonAttributes $(in)/populationAttributes.xml.gz\ + --ids $<\ + --requireAttribute "microm:modeled:age"\ + --output $@ + +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz: $(out)/bb_2020-week_snz_entirePopulation_emptyPlans_100pt.xml.gz + $(sc) districtLookup $<\ + --output $@\ + --shp ../public-svn/matsim/scenarios/countries/de/episim/original-data/landkreise-in-germany/landkreise-in-germany.shp + +######## +# 25pct +######## + +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_25pt.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_25pt.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz + $(sc) downSample 0.25\ + --population $<\ + --events $(in)/de2020gsmwt_events_reduced.xml.gz\ + --events $(in)/de2020gsmsa_events_reduced.xml.gz\ + --events $(in)/de2020gsmso_events_reduced.xml.gz\ + --output $(tmp) + + mv $(tmp)/population0.25.xml.gz $(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_25pt.xml.gz + mv $(tmp)/de2020gsmwt_events_reduced-0.25.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_25pt.xml.gz + mv $(tmp)/de2020gsmsa_events_reduced-0.25.xml.gz $(out)/bb_2020-week_snz_episim_events_sa_25pt.xml.gz + mv $(tmp)/de2020gsmso_events_reduced-0.25.xml.gz $(out)/bb_2020-week_snz_episim_events_so_25pt.xml.gz + +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_25pt_split.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_25pt_split.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_25pt.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_25pt.xml.gz + $(sc) splitHomeFacilities $<\ + --events $(out)/bb_2020-week_snz_episim_events_wt_25pt.xml.gz\ + --events $(out)/bb_2020-week_snz_episim_events_sa_25pt.xml.gz\ + --events $(out)/bb_2020-week_snz_episim_events_so_25pt.xml.gz\ + --shape-file $(out)/../shape-File/BerlinBrandenburg.shp\ + --output $(out) + +$(out)/wLeisure/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_25pt_split.xml.gz $(out)/wLeisure/bb_2020-week_snz_episim_events_wt_25pt_split.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_25pt.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_25pt.xml.gz + $(sc) splitHomeFacilities $<\ + --events $(out)/bb_2020-week_snz_episim_events_wt_25pt.xml.gz\ + --events $(out)/bb_2020-week_snz_episim_events_sa_25pt.xml.gz\ + --events $(out)/bb_2020-week_snz_episim_events_so_25pt.xml.gz\ + --shape-file $(out)/../shape-File/dilutionArea.shp\ + --remap visit --remap leisure\ + --output $(out)/wLeisure + +########### +# 100 pct +########### + +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt_filtered.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_100pt.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz + $(sc) downSample 1.0\ + --population $<\ + --events $(in)/de2020gsmwt_events_reduced.xml.gz\ + --events $(in)/de2020gsmsa_events_reduced.xml.gz\ + --events $(in)/de2020gsmso_events_reduced.xml.gz\ + --output $(tmp) + + mv $(tmp)/population1.0.xml.gz $(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt_filtered.xml.gz + mv $(tmp)/de2020gsmwt_events_reduced-1.0.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_100pt.xml.gz + mv $(tmp)/de2020gsmsa_events_reduced-1.0.xml.gz $(out)/bb_2020-week_snz_episim_events_sa_100pt.xml.gz + mv $(tmp)/de2020gsmso_events_reduced-1.0.xml.gz $(out)/bb_2020-week_snz_episim_events_so_100pt.xml.gz + +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt_split.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_100pt_split.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt_filtered.xml.gz $(out)/bb_2020-week_snz_episim_events_wt_100pt.xml.gz + $(sc) splitHomeFacilities $<\ + --events $(out)/bb_2020-week_snz_episim_events_wt_100pt.xml.gz\ + --events $(out)/bb_2020-week_snz_episim_events_sa_100pt.xml.gz\ + --events $(out)/bb_2020-week_snz_episim_events_so_100pt.xml.gz\ + --shape-file $(out)/../shape-File/dilutionArea.shp\ + --output $(out) + +################# +# 10pct 5pct 1pct +################# + + +$(out)/samples/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_10pt.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_10pt.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz + $(sc) downSample 0.1\ + --population $<\ + --events $(in)/de2020gsmwt_events_reduced.xml.gz\ + --events $(in)/de2020gsmsa_events_reduced.xml.gz\ + --events $(in)/de2020gsmso_events_reduced.xml.gz\ + --output $(tmp) + + mkdir $(out)/samples + + mv $(tmp)/population0.1.xml.gz $(out)/samples//bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_10pt.xml.gz + mv $(tmp)/de2020gsmwt_events_reduced-0.1.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_10pt.xml.gz + mv $(tmp)/de2020gsmsa_events_reduced-0.1.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_sa_10pt.xml.gz + mv $(tmp)/de2020gsmso_events_reduced-0.1.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_so_10pt.xml.gz + + +$(out)/samples/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_5pt_split.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_5pt_split.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz + $(sc) downSample 0.05\ + --population $<\ + --events $(in)/de2020gsmwt_events_reduced.xml.gz\ + --events $(in)/de2020gsmsa_events_reduced.xml.gz\ + --events $(in)/de2020gsmso_events_reduced.xml.gz\ + --output $(tmp) + + mkdir $(out)/samples + mv $(tmp)/population0.05.xml.gz $(out)/samples//bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_5pt_split.xml.gz + mv $(tmp)/de2020gsmwt_events_reduced-0.05.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_5pt_split.xml.gz + mv $(tmp)/de2020gsmsa_events_reduced-0.05.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_sa_5pt_split.xml.gz + mv $(tmp)/de2020gsmso_events_reduced-0.05.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_so_5pt_split.xml.gz + +$(out)/samples/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_1pt_split.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_1pt_split.xml.gz &: \ +$(out)/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz + $(sc) downSample 0.01\ + --population $<\ + --events $(in)/de2020gsmwt_events_reduced.xml.gz\ + --events $(in)/de2020gsmsa_events_reduced.xml.gz\ + --events $(in)/de2020gsmso_events_reduced.xml.gz\ + --output $(tmp) + + mkdir $(out)/samples + mv $(tmp)/population0.01.xml.gz $(out)/samples//bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_1pt_split.xml.gz + mv $(tmp)/de2020gsmwt_events_reduced-0.01.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_1pt_split.xml.gz + mv $(tmp)/de2020gsmsa_events_reduced-0.01.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_sa_1pt_split.xml.gz + mv $(tmp)/de2020gsmso_events_reduced-0.01.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_so_1pt_split.xml.gz + + +$(out)/samples/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_10pt_split.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_10pt_split.xml.gz &: \ +$(out)/samples/bb_2020-week_snz_entirePopulation_emptyPlans_withDistricts_10pt.xml.gz $(out)/samples/bb_2020-week_snz_episim_events_wt_10pt.xml.gz + $(sc) splitHomeFacilities $<\ + --events $(out)/samples/bb_2020-week_snz_episim_events_wt_10pt.xml.gz\ + --events $(out)/samples/bb_2020-week_snz_episim_events_sa_10pt.xml.gz\ + --events $(out)/samples/bb_2020-week_snz_episim_events_so_10pt.xml.gz\ + --shape-file $(out)/../shape-File/dilutionArea.shp\ + --output $(out)/samples \ No newline at end of file diff --git a/src/main/R/analyseCases.R b/src/main/R/analyseCases.R index e552fce58..7fa46b31c 100644 --- a/src/main/R/analyseCases.R +++ b/src/main/R/analyseCases.R @@ -39,13 +39,13 @@ imm_inf <- imm_inf_raw %>% # mutate(vax = generic + mRNA + vector + ba1Update + ba5Update + natural) ggplot() + #nShowingSymptoms # SARS_CoV_2 + geom_line(snap_inf, mapping = aes(date, nShowingSymptoms , group = seed, col = "base")) + geom_line(imm_inf, mapping = aes(date, nShowingSymptoms , group = seed, col = "imm-hist")) + - geom_line(snap_inf, mapping = aes(date, nShowingSymptoms , group = seed, col = "snapshot")) + scale_color_manual(name='Regression Model', - breaks=c('snapshot', 'imm-hist'), - values=c('snapshot'='red', 'imm-hist'='blue'))+ - labs(alt = "hello world") + - ggtitle("Infections") + breaks=c('base', 'imm-hist'), + values=c('base'='red', 'imm-hist'='blue'))+ + ggtitle("Cologne Cases") + + labs(x="Date", y="New Cases") diff --git a/src/main/R/gpsAnalyses.R b/src/main/R/gpsAnalyses.R new file mode 100644 index 000000000..a67431f61 --- /dev/null +++ b/src/main/R/gpsAnalyses.R @@ -0,0 +1,230 @@ +library("osmextract") +library("plyr") +library("tidyverse") +library("lubridate") +library("stringr") +library("sf") + + +setwd("/Users/sydney/root/svn/shared-svn/projects/episim/data/netcheck/data2") + +# read osm points +osm_points <- oe_read("/Users/sydney/Downloads/koeln-regbez-latest.osm.pbf", layer = "points", extra_tags = c("amenity", "shop", "leisure")) %>% + select(osm_id, shop, leisure, amenity) %>% + filter(!is.na(osm_id)) + +# read osm polygons +osm_multipolygons <- oe_read("/Users/sydney/Downloads/koeln-regbez-latest.osm.pbf", layer = "multipolygons") %>% + select(osm_way_id, shop, leisure, amenity, building, landuse) %>% + filter(!is.na(osm_way_id)) %>% + mutate(osm_way_id = as.numeric(osm_way_id)) + +# filter shop, amenity, leisure points +osm_points_shop <- osm_points %>% + filter(!is.na(shop)) %>% + select(c(osm_id, shop)) + +osm_points_amenity <- osm_points %>% + filter(!is.na(amenity)) %>% + select(c(osm_id, amenity)) + +osm_points_leisure <- osm_points %>% + filter(!is.na(leisure)) %>% + select(c(osm_id, leisure)) + +# filter building polygons +osm_multipolygons_buildings <- osm_multipolygons %>% + filter(!is.na(building)) %>% + select(osm_way_id) + +# join shop, amenity, leisure ponts with building polygons +sf_use_s2(FALSE) + +osm_multipolygons_with_shop_points <- osm_multipolygons_buildings %>% + st_join(osm_points_shop) %>% + filter(!is.na(osm_way_id)) %>% + filter(!is.na(shop)) %>% + select(c(osm_way_id, shop)) %>% + st_drop_geometry() + +osm_multipolygons_with_amenity_points <- osm_multipolygons_buildings %>% + st_join(osm_points_amenity) %>% + filter(!is.na(osm_way_id)) %>% + filter(!is.na(amenity)) %>% + select(c(osm_way_id, amenity)) %>% + st_drop_geometry() + +osm_multipolygons_with_leisure_points <- osm_multipolygons_buildings %>% + st_join(osm_points_leisure) %>% + filter(!is.na(osm_way_id)) %>% + filter(!is.na(leisure)) %>% + select(c(osm_way_id, leisure)) %>% + st_drop_geometry() + +osm_multipolygons <- osm_multipolygons %>% + st_drop_geometry() + +# merge with all polygons +osm_multipolygons_with_points <- osm_multipolygons %>% left_join(osm_multipolygons_with_shop_points, by = c("osm_way_id" = "osm_way_id")) +osm_multipolygons_with_points <- osm_multipolygons_with_points %>% left_join(osm_multipolygons_with_amenity_points, by = c("osm_way_id" = "osm_way_id")) +osm_multipolygons_with_points <- osm_multipolygons_with_points %>% left_join(osm_multipolygons_with_leisure_points, by = c("osm_way_id" = "osm_way_id")) + + +# remove some stuff from memory +rm(osm_multipolygons_buildings) +rm(osm_multipolygons) +rm(osm_points) +rm(osm_points_shop) +rm(osm_points_amenity) +rm(osm_points_leisure) +rm(osm_multipolygons_with_shop_points) +rm(osm_multipolygons_with_amenity_points) +rm(osm_multipolygons_with_leisure_points) + + +# osm_lines <- oe_read("/Users/sebastianmuller/git/koeln-regbez-latest.osm.pbf", layer = "lines") + +# osm_multilinestrings <- oe_read("/Users/sebastianmuller/git/koeln-regbez-latest.osm.pbf", layer = "multilinestrings") + +# osm_other_relations <- oe_read("/Users/sebastianmuller/git/koeln-regbez-latest.osm.pbf", layer = "other_relations") + +#read bank holidays +bankHolidays = read.csv("/Users/sydney/git/matsim-episim-libs/src/main/resources/bankHolidays.csv", header = TRUE) %>% + mutate(date = as.Date(bankHoliday)) + +#read senozon activity reductions +snz <- read.delim("/Users/sydney/root/svn/shared-svn/projects/episim/matsim-files/snz/Cologne/episim-input/CologneSnzData_daily_until20221205.csv", header = TRUE, sep = "\t") %>% + select(-c(notAtHomeExceptLeisureAndEdu, notAtHomeExceptEdu, notAtHome_22, accomp, traveling, undefined, visit, home)) %>% + pivot_longer(!date, names_to = "act", values_to = "reduction") %>% + mutate(newDate = as.Date(strptime(date, "%Y%m%d"))) %>% + left_join(bankHolidays, by = c("newDate" = "date")) %>% + mutate(weekday = wday(newDate, week_start = 1)) %>% + filter(weekday < "6") %>% + mutate(holiday = "no") %>% + mutate(holiday = ifelse(str_count(Bundesland, "Germany") > 0, "yes", holiday )) %>% + mutate(holiday = ifelse(str_count(Bundesland, "Nordrhein") > 0, "yes", holiday )) %>% + mutate(holiday = ifelse(is.na(holiday), "no", holiday )) %>% + filter(holiday == "no") %>% + mutate( week = paste0(isoweek(newDate), "-", isoyear(newDate))) %>% + group_by( week, act ) %>% + summarize( actReduction=mean(reduction) * 0.01, newDate=mean(newDate)) %>% + ungroup() + + + +#read all netcheck csv files +netcheck <- ldply( .data = list.files(path ="/Users/sydney/root/svn/shared-svn/projects/episim/data/netcheck/data2", pattern="*_v2.csv"), + .fun = read.csv, + header = TRUE) + +# add an ID +netcheck <- tibble::rowid_to_column(netcheck, "ID") + +# merge netcheck with bankholidays +netcheck <- netcheck %>% left_join(bankHolidays, by = c("day" = "bankHoliday")) + +netcheck <- netcheck %>% + mutate(date = as.Date(day)) %>% + mutate(weekday = "Mon-Fri") %>% + mutate(weekday = ifelse(wday(date) == 1, "Sun", weekday )) %>% + mutate(weekday = ifelse(wday(date) == 7, "Sat", weekday )) %>% + mutate(holiday = "no") %>% + mutate(holiday = ifelse(str_count(Bundesland, "Germany") > 0, "yes", holiday )) %>% + mutate(holiday = ifelse(str_count(Bundesland, "Nordrhein") > 0, "yes", holiday )) %>% + mutate(holiday = ifelse(is.na(holiday), "no", holiday)) %>% + mutate(nTimeStamps = 1 + str_count(timestamps, ",")) + +# calcualte duration +netcheck <- netcheck %>% + mutate(maxTime = strptime(substr(timestamps, nchar(timestamps) - 7, nchar(timestamps)), format="%H:%M:%S")) %>% + mutate(minTime = strptime(substr(timestamps, 1, 8), format="%H:%M:%S")) %>% + mutate(duration = maxTime - minTime) + +sumDurations <- as.numeric(sum(netcheck$duration)) +sumTimePeriods <- sum(netcheck$nTimeStamps) - nrow(netcheck) +timePerTimePeriod <- sumDurations / sumTimePeriods + +netcheck <- netcheck %>% + mutate(duration = timePerTimePeriod + as.numeric(maxTime - minTime)) + + +merged_multipolygons <- netcheck %>% inner_join(osm_multipolygons_with_points, by = c("osm_id" = "osm_way_id")) + +merged_multipolygons <- merged_multipolygons %>% + group_by(ID) %>% + add_count(name = "id_occurrence") %>% + mutate(weightedTimeStamps = nTimeStamps / id_occurrence) %>% + mutate(weightedDurations = duration / id_occurrence) %>% + ungroup() + +# remove some stuff from memory +#rm(netcheck) +#rm(osm_multipolygons_with_points) + +# number of pings per day, needed for scaling +all <- merged_multipolygons %>% + group_by(date, weekday, holiday) %>% + summarise(sumTimeStamps = sum(weightedTimeStamps), sumDurations = sum(weightedDurations)) %>% + ungroup() + +# not at home +notAtHome <- merged_multipolygons %>% + # filter(distance_to_home > 0) + filter(landuse != "residential" | is.na(landuse)) + +colnames(notAtHome)[10] <- "date" +colnames(all)[1] <- "date" + +notAtHome <- notAtHome %>% + group_by(date) %>% + summarise(sumTimeStamps = sum(weightedTimeStamps), sumDurations = sum(weightedDurations)) %>% + ungroup() + +notAtHome <- notAtHome %>% full_join(all, by = c("date" = "date")) %>% + mutate(sumTimeStamps.x = ifelse(is.na(sumTimeStamps.x), 0, sumTimeStamps.x)) %>% + mutate(sumDurations.x = ifelse(is.na(sumDurations.x), 0, sumDurations.x)) + +notAtHome$normTimeStamps <- notAtHome$sumTimeStamps.x / notAtHome$sumTimeStamps.y +notAtHome$normDurations <- notAtHome$sumDurations.x / notAtHome$sumDurations.y + +baseline <- notAtHome %>% + filter(date >= as.Date("2020-09-01") & date <= as.Date("2020-09-30")) %>% + # filter(date >= as.Date("2020-03-01")) %>% + filter(holiday == "no") %>% + group_by(weekday) %>% + summarise(baseTimeStamps = median(normTimeStamps), baseDurations = median(normDurations)) %>% + # summarise(baseTimeStamps = quantile(normTimeStamps, 0.9), baseDurations = quantile(normDurations, 0.9)) %>% + ungroup() + +scaleMonFriTimeStamps <- baseline$baseTimeStamps[1] +scaleSatTimeStamps <- baseline$baseTimeStamps[2] +scaleSunTimeStamps <- baseline$baseTimeStamps[3] +scaleMonFriDurations <- baseline$baseDurations[1] +scaleSatDurations <- baseline$baseDurations[2] +scaleSunDurations <- baseline$baseDurations[3] + +notAtHome <- notAtHome %>% + mutate(normTimeStamps2 = -1) %>% + mutate(normTimeStamps2 = ifelse(weekday == "Mon-Fri", (normTimeStamps - scaleMonFriTimeStamps) / scaleMonFriTimeStamps, normTimeStamps2)) %>% + mutate(normTimeStamps2 = ifelse(weekday == "Sat", (normTimeStamps - scaleSatTimeStamps) / scaleSatTimeStamps, normTimeStamps2)) %>% + mutate(normTimeStamps2 = ifelse(weekday == "Sun", (normTimeStamps - scaleSunTimeStamps) / scaleSunTimeStamps, normTimeStamps2)) %>% + mutate(normDurations2 = -1) %>% + mutate(normDurations2 = ifelse(weekday == "Mon-Fri", (normDurations - scaleMonFriDurations) / scaleMonFriDurations, normDurations2)) %>% + mutate(normDurations2 = ifelse(weekday == "Sat", (normDurations - scaleSatDurations) / scaleSatDurations, normDurations2)) %>% + mutate(normDurations2 = ifelse(weekday == "Sun", (normDurations - scaleSunDurations) / scaleSunDurations, normDurations2)) + +notAtHomeWeekly <- notAtHome %>% + filter(weekday == "Mon-Fri") %>% + filter(holiday == "no") %>% + mutate( week = paste0(isoweek(date), "-", isoyear(date))) %>% + group_by( week ) %>% + summarize( sumTimeStamps.x=mean(sumTimeStamps.x), sumDurations.x=mean(sumDurations.x), date=mean(date), normTimeStamps=mean(normTimeStamps), normDurations=mean(normDurations), normTimeStamps2=mean(normTimeStamps2), normDurations2=mean(normDurations2)) + + + + +sort(table(merged_multipolygons$name, useNA="always"), decreasing = TRUE) + +table(merged_multipolygons$ID, useNA="always") + + diff --git a/src/main/R/gpsPlots.R b/src/main/R/gpsPlots.R new file mode 100644 index 000000000..34a2ebdb6 --- /dev/null +++ b/src/main/R/gpsPlots.R @@ -0,0 +1,564 @@ +library(gridExtra) +library(ggiraphExtra) +library(xlsx) +library(tidyverse) + +source("/Users/sydney/git/matsim-episim-libs/src/main/R/gpsAnalyses.R") + +options(scipen = 999) #Use this, to NOT display axis labeling of plots in scientific notation + +#Renaming of column and +colnames(netcheck)[6] <- "Cologne Resident" +netcheck$`Cologne Resident`[netcheck$`Cologne Resident` == ""] <- "Unknown" +netcheck$holiday[netcheck$holiday == "no"] <- "False" +netcheck$holiday[netcheck$holiday == "yes"] <- "True" + +# First exploratory plots +# Counting no of timestamps over time, differentiatd by whether data point is provided by a resident of the city of Cologne +netcheck %>% + group_by(date, `Cologne Resident`) %>% + summarise(count = sum(nTimeStamps)) %>% + ungroup() %>% + ggplot(aes(x = date, y = count, color = `Cologne Resident`)) + + geom_point(size = 2) + + theme_minimal() + + theme(legend.position = "bottom") + + theme(text = element_text(size = 13)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + scale_x_date(date_breaks = "3 month", date_labels = "%d/%b/%y") + + xlab("Date") + + ylab("Counts") + + scale_color_brewer(palette = "Set2") + +#Similar plot as above, here differentiation by weekday and holiday (TRUE/FALSE) +netcheck %>% + group_by(date, weekday, holiday) %>% + summarise(count = sum(nTimeStamps)) %>% + ungroup() %>% + ggplot(aes(x = date, y = count, color = weekday)) + + geom_point(size = 1.5) + + theme_minimal() + + theme(legend.position = "bottom", legend.title = element_blank()) + + theme(text = element_text(size = 17)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + theme(axis.text.x = element_text(angle = 45, vjust = 1)) + + scale_x_date(date_breaks = "2 month", date_labels = "%d/%b/%y", limit=c(as.Date("2019-10-01"),as.Date("2021-03-01"))) + + scale_y_continuous(labels = scales::comma) + + xlab("Date") + + ylab("Counts") + + scale_color_brewer(palette = "Dark2") + +ggsave("timepoints.pdf", dpi = 500, w = 9, h = 4.5) + +#Plots for results section +filteredForPlot <- merged_multipolygons %>% + #restaurant + # filter(amenity.x == "restaurant" | amenity.x == "bar" | amenity.x == "biergarten" | amenity.y == "restaurant" | amenity.y == "bar" | amenity.y == "biergarten") + + # other leisure + # filter(amenity == "theatre" | amenity == "cinema") + + # private leisure + # filter(distance_to_home > 0 & landuse == "residential") + # filter(distance_to_home > 0 & (building == "apartments" | building == "house" | building == "semidetached_house")) + + #school + filter(amenity.x == "school" | amenity.y == "school") +# filter(landuse == "education") + +#university +# filter(amenity.x == "university" | amenity.y == "university") + +#kindergarten +# filter(amenity == "kindergarten") + +# work +# filter(distance_to_work < 1) +# filter(landuse == "commercial" | landuse == "industrial") +# filter(landuse == "commercial") +# filter(building == "office" | building == "commercial" | building == "industrial") + +# home +# filter(distance_to_home != 0) +# filter(landuse == "residential") + +# not home +# filter(distance_to_home > 0) + +# shop +#filter(!is.na(shop.x) | !is.na(shop.y)) %>% +# filter(shop.x != "supermarket" | shop.y != "supermarket") +# filter(shop.x == "mall" | shop.y == "mall") +# filter(landuse == "retail") +# filter(shop.x == "hairdresser" | shop.y == "hairdresser") + +# park +#filter(leisure.x == "park") + +#errands + +# S: I believe the following two lines are unnecessary and can be removed +#shareTimeStamps = 100 * sum(filteredForPlot$weightedTimeStamps) / sum(merged_multipolygons$weightedTimeStamps) +#shareDurations = 100 * sum(filteredForPlot$weightedDurations) / sum(merged_multipolygons$weightedDurations) + +#For each day the no. of time stamps are counted +filteredForPlot <- filteredForPlot %>% + group_by(date) %>% + summarise(sumTimeStamps = sum(weightedTimeStamps), sumDurations = sum(weightedDurations)) %>% + ungroup() + +#If for a certain date no data exists, then the corresponding entry of the data frame is set equal to 0 +filteredForPlot <- filteredForPlot %>% full_join(all, by = c("date" = "date")) %>% + mutate(sumTimeStamps.x = ifelse(is.na(sumTimeStamps.x), 0, sumTimeStamps.x)) %>% + mutate(sumDurations.x = ifelse(is.na(sumDurations.x), 0, sumDurations.x)) + +#For a given category, its share (of time stamps and durations) of all activity types is computed +filteredForPlot$normTimeStamps <- filteredForPlot$sumTimeStamps.x / filteredForPlot$sumTimeStamps.y +filteredForPlot$normDurations <- filteredForPlot$sumDurations.x / filteredForPlot$sumDurations.y + +baseline <- filteredForPlot %>% + filter(date >= as.Date("2020-09-01") & date <= as.Date("2020-09-30")) %>% + # filter(date >= as.Date("2020-03-01")) %>% + filter(holiday == "no") %>% + group_by(weekday) %>% + summarise(baseTimeStamps = median(normTimeStamps), baseDurations = median(normDurations)) %>% + #summarise(baseTimeStamps = quantile(normTimeStamps, 0.9), baseDurations = quantile(normDurations, 0.9)) %>% + ungroup() + +scaleMonFriTimeStamps <- baseline$baseTimeStamps[1] +scaleSatTimeStamps <- baseline$baseTimeStamps[2] +scaleSunTimeStamps <- baseline$baseTimeStamps[3] +scaleMonFriDurations <- baseline$baseDurations[1] +scaleSatDurations <- baseline$baseDurations[2] +scaleSunDurations <- baseline$baseDurations[3] + +filteredForPlot <- filteredForPlot %>% + mutate(normTimeStamps2 = -1) %>% + mutate(normTimeStamps2 = case_when(weekday == "Mon-Fri" ~ (normTimeStamps - scaleMonFriTimeStamps) / scaleMonFriTimeStamps, + weekday == "Sat" ~ (normTimeStamps - scaleSatTimeStamps) / scaleSatTimeStamps, + weekday == "Sun" ~ (normTimeStamps - scaleSunTimeStamps) / scaleSunTimeStamps, + .default = normTimeStamps2)) %>% + mutate(normDurations2 = case_when(weekday == "Mon-Fri" ~ (normDurations - scaleMonFriDurations) / scaleMonFriDurations, + weekday == "Sat" ~ (normDurations - scaleSatDurations) / scaleSatDurations, + weekday == "Sun" ~ (normDurations - scaleSunDurations) / scaleSunDurations, + .default = normTimeStamps2)) + +#Aggregation of daily data on a weekly level +filteredForPlotWeekly <- filteredForPlot %>% + filter(weekday == "Mon-Fri") %>% + filter(holiday == "no") %>% + mutate(week = paste0(isoweek(date), "-", isoyear(date))) %>% + group_by(week) %>% + summarize(sumTimeStamps.x=mean(sumTimeStamps.x), sumDurations.x=mean(sumDurations.x), date=mean(date), normTimeStamps=mean(normTimeStamps), normDurations=mean(normDurations), normTimeStamps2=mean(normTimeStamps2), normDurations2=mean(normDurations2)) %>% + ungroup() + + +# NETCHECK PLOTS ----------------------------------------------------------- + +ggplot() + + geom_point(data=filteredForPlot, aes(x=date, y=sumTimeStamps.x, color = weekday, shape = holiday)) + + theme_minimal() + + theme(legend.position = "bottom") + + theme(text = element_text(size = 13)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y", limit=c(as.Date("2020-03-01"),as.Date("2021-03-24"))) + + xlab("Date") + + ylab("Counts") + + scale_color_brewer(palette = "Set2") + +ggplot() + + geom_point(data=filteredForPlotWeekly, aes(x=date, y=normTimeStamps2), color = "#1b9e77") + + geom_point(data=notAtHomeWeekly, aes(x=date, y=normTimeStamps2), color = "#d95f02") + + #geom_point(data=snz, aes(x=newDate, y=actReduction), color="darkgrey") + + labs( + # title="Time Stamps", + # caption="red: act (nc), green: out of home (nc), blue: out of home (snz)", + caption = "Time stamps: activity (green), out of home (orange)", + # x="date", y="Red. vs. last week in Feb. 2020") + + x = "Date", y = "Red. vs. 90th percentile") + + scale_y_continuous(labels = scales::percent, limit=c(-1.0, 0.3)) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y", limit=c(as.Date("2020-03-01"),as.Date("2021-03-24"))) + + theme_minimal() + + theme(legend.position = "bottom") + + theme(text = element_text(size = 13)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + +ggsave("ts22.png", width = 4, height = 3.2) + +# NETCHECK AND SENOZON PLOTS ----------------------------------------------------------- +colnames(snz)[2] <- "Activity" + +#Plot of +#1) Activity reduction of chosen activity type (netcheck and senozon) +#2) Reduction of out of home duration (senozon) +colors <- c("Not At Home (cell-based)" = "dimgrey", "School (cell-based)" = "darkgrey", "School (GPS-based)" = "#1b9e77") + +#School Vacation +date_breaks <- data.frame(start = c(as.Date("2020-04-06"), as.Date("2020-06-29"), as.Date("2020-10-12"), as.Date("2020-12-21")), + end = c(as.Date("2020-04-18"), as.Date("2020-08-11"), as.Date("2020-10-24"), as.Date("2021-01-06")), + colors = c("Vacation", "Vacation", "Vacation", "Vacation")) + + +school <- ggplot() + + geom_rect(data = date_breaks, + aes(xmin = start, + xmax = end, + ymin = - Inf, + ymax = Inf, + fill = colors), + alpha = 0.3) + + scale_fill_manual(values = c("#B9BBB6", "#B9BBB6", "#B9BBB6", "#B9BBB6")) + + geom_point(data=filteredForPlotWeekly, aes(x=date, y=normTimeStamps2, color ="School (GPS-based)"), size = 2) + + geom_line(data=snz %>% filter(Activity == "notAtHome"), aes(x=newDate, y=actReduction, color = "Not At Home (cell-based)"), size = 2) + + geom_point(data=snz %>% filter(Activity == "education"), aes(x=newDate, y=actReduction, color = "School (cell-based)"), size = 2) + + labs( + caption = "", + x = "Date", +# y = "") + + y = "Activity Reduction") + + geom_vline(xintercept = as.Date("2020-03-13"), linetype = "dotted") + #School closure + geom_vline(xintercept = as.Date("2020-12-14"), linetype = "dotted") + #School closure + scale_y_continuous(labels = scales::percent, limit = c(-1.0, 0.3)) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y", limit=c(as.Date("2020-03-01"),as.Date("2021-03-01"))) + + theme_minimal() + + theme(legend.position = "bottom", legend.title = element_blank()) + + theme(text = element_text(size = 20)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + theme(axis.text.x = element_text(angle = 45, vjust = 1)) + + scale_color_manual(values = colors) + + guides(color=guide_legend(nrow=2, byrow = TRUE)) + +p <- arrangeGrob(school,work, nrow=1) + +ggsave("SchoolWorkActivityRed.pdf", p, dpi = 500, width = 15.75, height = 5) + +#Plot of +#1) "Not at home" reduction via reduction of # of time stamps (netcheck) +#2) "Not at home" reduction via reduction of out of home duration (netcheck) +#3) "Not at home" reduction via reduction of out of home duration (senozon) +ggplot() + + geom_point(data = snz, aes(x = newDate, y = actReduction), color = "#666666") + + geom_point(data=notAtHomeWeekly, aes(x=date, y=normDurations2), color="#66C2A5") + + geom_point(data=notAtHomeWeekly, aes(x=date, y=normTimeStamps2), color="#FC8D62") + + labs( + caption="Out of home gps time stamps (orange) and durations (green), out of home mobile phone (grey)", + x = "Date", + y = "Red. vs. last week in Feb. 2020") + + scale_y_continuous(labels = scales::percent, limit=c(-0.7, 0.3)) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y", limit=c(as.Date("2020-03-01"),as.Date("2021-03-24"))) + + theme_minimal() + + theme(legend.position = "bottom") + + theme(text = element_text(size = 13)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + +ggsave("outOfHome.png", width = 8, height = 4) + +# SENOZON PLOTS ----------------------------------------------------------- + +# Data prep +# Include column that contains information whether or not a given date is a school holiday +snz <- snz %>% + mutate(schoolVac = "no") %>% + mutate(schoolVac = ifelse(newDate > as.Date("2020-04-05") & newDate < as.Date("2020-04-19"), "yes", schoolVac)) %>% + mutate(schoolVac = ifelse(newDate > as.Date("2020-06-28") & newDate < as.Date("2020-08-12"), "yes", schoolVac)) %>% + mutate(schoolVac = ifelse(newDate > as.Date("2020-10-11") & newDate < as.Date("2020-10-25"), "yes", schoolVac)) %>% + mutate(schoolVac = ifelse(newDate > as.Date("2020-12-22") & newDate < as.Date("2021-01-07"), "yes", schoolVac)) + +# Plot depicting reduction of out of home duration (line) and reduction of school activity (dots) +colnames(snz)[5] <- "School Vacation" +snz <- snz %>% mutate(`School Vacation` = case_when(`School Vacation` == "yes" ~ "education, during school vacation", + `School Vacation` == "no" ~ "education, no school vacation")) +ggplot() + + geom_line(data = snz %>% filter(Activity == "notAtHome") , aes(x=newDate, y=actReduction, linetype = Activity), color="#666666", size = 1.5) + + geom_point(data = snz %>% filter(Activity == "education") , aes(x=newDate, y=actReduction, color=`School Vacation`), size = 1.5) + + geom_vline(xintercept = as.Date("2020-03-13"), linetype = "dotted") + + geom_vline(xintercept = as.Date("2020-12-14"), linetype = "dotted") + + labs(x="Date", y="Red. vs. before pandemic") + + scale_y_continuous(labels = scales::percent, limit = c(-0.6, 0.1)) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y", limit=c(as.Date("2020-03-01"),as.Date("2021-03-01"))) + + theme_minimal() + + theme(legend.position = "bottom", legend.title = element_blank()) + + theme(text = element_text(size = 17)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + theme(legend.position = "bottom") + + scale_color_brewer(palette = "Dark2") + + theme(axis.text.x = element_text(angle = 45, vjust = 1)) + + guides(color = guide_legend(nrow = 2)) + +ggsave("snzSchool.pdf", width = 8, height = 5, dpi = 500) + +# Plot depicting reduction of out of home duration (line) and reduction of all other activities (dots) +snz$Activity <- factor(snz$Activity, levels = c("business", "education", "errands", "leisure", "shop_daily", "shop_other", "work", "notAtHome")) +ggplot(snz, aes(x = newDate, y = actReduction, colour = Activity, + linetype = Activity, shape = Activity), size = 1.5) + +geom_point(size = 1.5) + +geom_line(size = 1.5) + + scale_linetype_manual(values=c(NA, NA, NA, NA, NA, NA, NA, "solid")) + + scale_shape_manual(values=c(16, 16, 16, 16, 16, 16, 16, NA)) + + #geom_point(data = snz %>% filter(Activity != "notAtHome") , aes(x = newDate, y = actReduction, color = Activity)) + + # geom_line(data = snz %>% filter(Activity == "notAtHome"), aes(x = newDate, y = actReduction), color = "#666666", size = 1.2) + + labs(x = "Date", y = "Reduction vs. before pandemic") + + scale_y_continuous(labels = scales::percent, limit=c(-0.6, 0.1)) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y", limit=c(as.Date("2020-03-01"),as.Date("2021-03-01"))) + + theme_minimal() + + theme(legend.position = "bottom", legend.title = element_blank()) + + theme(text = element_text(size = 17)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + scale_color_brewer(palette = "Dark2") + + theme(axis.text.x = element_text(angle = 45, vjust = 1)) + +ggsave("snzAll.pdf", width = 8, height = 5.5, dpi = 500) + + +# GOOGLE MOBILITY REPORTS PLOTS ----------------------------------------------------------- + +# Google Mobility Report +googleMobilityReportRaw <- read_csv("/Users/sebastianmuller/Downloads/Global_Mobility_Report.csv") + +googleMobilityReport <- googleMobilityReportRaw %>% + filter(sub_region_1 == "North Rhine-Westphalia") %>% + select(-c(country_region_code, country_region, sub_region_1, sub_region_2, metro_area, iso_3166_2_code, census_fips_code, place_id)) %>% + mutate(notAtHome = -2 * residential_percent_change_from_baseline) %>% + pivot_longer(!date, names_to = "type", values_to = "restriction") %>% + # filter(type == "retail_and_recreation_percent_change_from_baseline" | type == "workplaces_percent_change_from_baseline" | type == "grocery_and_pharmacy_percent_change_from_baseline" | type == "parks_percent_change_from_baseline" ) %>% + filter(type == "workplaces_percent_change_from_baseline") %>% + #filter(type == "retail_and_recreation_percent_change_from_baseline" | type == "workplaces_percent_change_from_baseline" | type == "notAtHome") %>% + #filter(type == "retail_and_recreation_percent_change_from_baseline" | type == "parks_percent_change_from_baseline") %>% + mutate(weekday = wday(date, week_start = 1)) %>% + filter(weekday < "6") %>% + mutate(week = paste0(isoweek(date), "-", isoyear(date))) %>% + group_by(week, type) %>% + summarize(restriction = 0.01 * mean(restriction), date = mean(date)) + +# ggplot(data = googleMobilityReport, mapping=aes(x = date)) + +# labs( +# title="Google mobility report (NRW)", +# caption="Source: Google", +# x="date", y="Reduction in %") + +# geom_point(mapping=aes(y = restriction, colour = type)) + +# ylim(-75, 50) + +# xlim(c(as.Date("2020-03-01"), as.Date("2021-02-28"))) + +# theme(legend.position = "bottom") + +ggplot () + + geom_point(data=filteredForPlotWeekly, aes(x=date, y=normDurations2), color="red") + + # geom_point(data=notAtHomeWeekly, aes(x=date, y=normDurations2), color="blue") + + geom_point(data=googleMobilityReport, aes(x=date, y=restriction), color="darkgreen") + + # geom_point(data=snz, aes(x=newDate, y=actReduction), color="darkgrey") + + labs( + # title="Duration", + # caption="red: act (nc), green: out of home (nc), blue: out of home (snz)", + caption = "durations: activity (red), work google (green)", + x = "date", y = "Red. vs. last week in Feb. 2020") + + scale_x_date(date_labels = "%m-%Y", limit=c(as.Date("2020-01-01"),as.Date("2021-03-24")), date_breaks = "3 months") + + scale_y_continuous(labels = scales::percent, limit = c(-0.8, 0.4)) + +ggsave("dur.png", width = 4, height = 3.2) + +# WEATHER PLOTS ----------------------------------------------------------- +weatherData <- read_delim("https://bulk.meteostat.net/daily/10513.csv.gz", delim = ",", col_names = FALSE, col_types = cols( + X1 = col_date(format = ""), + X2 = col_double(), + X3 = col_double(), + X4 = col_double(), + X5 = col_double(), + X6 = col_double(), + X7 = col_double(), + X8 = col_double(), + X9 = col_double(), + X10 = col_double(), + X11 = col_double() +)) + +colnames(weatherData) <- c("date", "tavg", "tmin", "tmax", "prcp", "snow", "wdir", "wspd", "wpgt", "pres", "tsun") + +weatherDataByWeek <- weatherData %>% + mutate(week = paste0(isoweek(date), "-", isoyear(date))) %>% + group_by(week) %>% + summarize( date=mean(date), tavg=mean(tavg), tmin=mean(tmin), tmax=mean(tmax), prcp=mean(prcp), snow=mean(snow), wdir=mean(wdir), wspd=mean(wspd), wpgt=mean(wpgt), pres=mean(pres), tsun=mean(tsun)) + + +#Plot of Daily maximum temperature over time and reduction of chosen activity type (GPS based data) +ggplot() + + geom_point(data=filteredForPlotWeekly, aes(x=date, y=normTimeStamps2), color="#FC8D62") + + # geom_point(data=notAtHomeWeekly, aes(x=date, y=normDurations2), color="blue") + + geom_point(data=weatherDataByWeek, aes(x=date, y= -0.7 + tmax * 0.03), color="darkgreen") + + # geom_point(data=snz, aes(x=newDate, y=actReduction), color="darkgrey") + + labs( + caption = "time stamps: activity (red), tmax (green, scaled)", + x = "Date", + y = "Red. vs. 90th percentile") + + scale_y_continuous(labels = scales::percent, limit = c(-1.0, 0.3)) + + scale_x_date(date_breaks = "3 month", date_labels = "%d/%b/%y", limit=c(as.Date("2020-03-01"),as.Date("2021-03-24"))) + + theme_minimal() + + theme(legend.position = "bottom") + + theme(text = element_text(size = 13)) + + theme(axis.ticks.x = element_line(), + axis.ticks.y = element_line(), + axis.ticks.length = unit(5, "pt")) + + theme(legend.position = "bottom") + + scale_color_brewer(palette = "Set2") + +ggsave("weather.png", width = 4, height = 3.2) + + + +#Figure 5 --> Time stamps at OSM facilities +amenity <- read.xlsx("/Users/sydney/Downloads/timestampsOSM.xlsx", sheetIndex = 1) +amenity <- amenity[1:11,] +amenity <- amenity[,-1] +amenity <- amenity[,-6] +colnames(amenity) <- c("amenity", "timepointsDatasetA", "timepointsDatasetB", "timepointsAll", "share") +amenity$amenity <- factor(amenity$amenity, levels = c("other", "social_facility", "pharmacy", "doctors", "cafe", "school", "fast_food", "atm", "restaurant", "hospital", "parking")) +amenityplot <- ggplot(amenity) + + geom_col(aes(timepointsAll/1000, amenity), fill = "#1b9e77", width = 0.6) + + theme_minimal() + + theme(text = element_text(size = 19)) + + scale_x_continuous() + + ylab("Amenity") + + scale_x_continuous(labels = scales::comma, breaks = c(0,1000,2000,3000)) + + xlab("Thousand Timestamps") +ggsave("amenity.pdf", amenityplot, dpi = 500, w = 6, h = 3) + +shop <- read.xlsx("/Users/sydney/Downloads/timestampsOSM.xlsx", sheetIndex = 2) +shop <- shop[1:11,] +shop <- shop[,-6] +colnames(shop) <- c("shop", "timepointsDatasetA", "timepointsDatasetB", "timepointsAll", "share") +shop$shop <- factor(shop$shop, levels = c("other", "electronics", "department__store", "hairdresser", "do_it_yourself", "kiosk", "bakery", "furniture", "clothes", "mall", "supermarket")) +shopplot <- ggplot(shop) + + geom_col(aes(timepointsAll/1000, shop), fill = "#1b9e77", width = 0.6) + + theme_minimal() + + theme(text = element_text(size = 19)) + + scale_x_continuous() + + ylab("Shop") + + scale_x_continuous(labels = scales::comma, breaks = c(0,500, 1000,1500,2000)) + + xlab("Thousand Timestamps") +ggsave("shop.pdf", shopplot, dpi = 500, w = 6, h = 3) + +leisure <- read.xlsx("/Users/sydney/Downloads/timestampsOSM.xlsx", sheetIndex = 3) +leisure <- leisure[1:11,] +leisure <- leisure[,-6] +colnames(leisure) <- c("leisure facility", "timepointsDatasetA", "timepointsDatasetB", "timepointsAll", "share") +leisure$`leisure facility` <- factor(leisure$`leisure facility`, levels = c("other", "common", "pitch", "dance", "golf_course", "garden", "stadium", "playground", "fitness_centre", "sports_centre", "park")) +leisureplot <- ggplot(leisure) + + geom_col(aes(timepointsAll/1000, `leisure facility`), fill = "#1b9e77", width = 0.6) + + theme_minimal() + + theme(text = element_text(size = 19)) + + scale_x_continuous() + + ylab("Leisure") + + scale_x_continuous(labels = scales::comma, breaks = c(0,250, 500, 750, 1000,1250, 1500, 1750,2000)) + + xlab("Thousand Timestamps") +ggsave("leisure.pdf", leisureplot, dpi = 500, w = 6, h = 3) + + +landuse <- read.xlsx("/Users/sydney/Downloads/timestampsOSM.xlsx", sheetIndex = 4) +landuse <- landuse[1:11,] +landuse <- landuse[,-6] +colnames(landuse) <- c("landuse", "timepointsDatasetA", "timepointsDatasetB", "timepointsAll", "share") +landuse$landuse <- factor(landuse$landuse, levels = c("other", "farmyard", "farmland", "allotments", "forest", "grass", "railway", "retail", "industrial", "commercial", "residential")) +ggplot(landuse) + + geom_col(aes(timepointsAll/1000, landuse), fill = "#1b9e77", width = 0.6) + + theme_minimal() + + theme(text = element_text(size = 19)) + + scale_x_continuous() + + ylab("Landuse Activity") + + scale_x_continuous(labels = scales::comma, breaks = c(0,25000, 50000, 75000, 100000)) + + xlab("Thousand Timestamps") + +ggsave("landuse.pdf", dpi = 500, w = 6, h = 3) + +building <- read.xlsx("/Users/sydney/Downloads/timestampsOSM.xlsx", sheetIndex = 5) +building <- building[1:11,] +building <- building[,-6] +colnames(building) <- c("building", "timepointsDatasetA", "timepointsDatasetB", "timepointsAll", "share") +building$building <- factor(building$building, levels = c("other", "commercial", "industrial", "hospital", "semidetachedhouse", "office", "house", "retail", "terrace", "apartments", "yes")) +buildingplot <- ggplot(building) + + geom_col(aes(timepointsAll/1000, building), fill = "#1b9e77", width = 0.6) + + theme_minimal() + + theme(text = element_text(size = 19)) + + scale_x_continuous() + + ylab("Building") + + scale_x_continuous(labels = scales::comma, breaks = c(0,3000,6000,9000)) + + xlab("Thousand Timestamps") +ggsave("building.pdf", buildingplot, dpi = 500, w = 6, h = 3) + + +#Figure 6 --> Demonstration of the data standardization process +TimeStamps <- read.xlsx("/Users/sydney/Downloads/Standardizationtimestmaps.xlsx", sheetIndex = 1) +TimeStamps <- TimeStamps[1:13, 1:5] +TimeStamps <- pivot_longer(TimeStamps, cols = c("Act.A", "Act.B", "Other", "All")) +TimeStamps <- TimeStamps %>% mutate(name = case_when(name == "Act.A" ~ "Activity A", + name == "Act.B" ~ "Activity B", + name == "Other" ~ "Other", + name == "All" ~ "All")) +#TimeStamps$Date <- as.Date(TimeStamps$Date) +#TimeStamps$Date <- TimeStamps$Date + 25 +ggplot(TimeStamps %>% filter(name != "All"), aes(fill = name, y = value, x = Date)) + + geom_bar(position="stack", stat="identity") + + theme_minimal() + + theme(text = element_text(size = 21)) + + theme(legend.position = "bottom", legend.title = element_blank()) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y") + + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + + ylab("Timestamps") + + xlab("Date") + + scale_fill_brewer(palette = "Dark2") + +ggsave("TimestampsLeft.pdf", dpi = 500, w = 5.25, h = 5.25) + +TimeStampsRelative <- read.xlsx("/Users/sydney/Downloads/Standardizationtimestmaps.xlsx", sheetIndex = 2) +TimeStampsRelative <- pivot_longer(TimeStampsRelative, cols = c("Activity.A", "Activity.B", "Other")) +TimeStampsRelative <- TimeStampsRelative %>% mutate(name = case_when(name == "Activity.A" ~ "Activity A", + name == "Activity.B" ~ "Activity B", + name == "Other" ~ "Other")) + +ggplot(TimeStampsRelative, aes(fill = name, y = value, x = Date)) + + geom_bar(position="fill", stat="identity") + + theme_minimal() + + theme(text = element_text(size = 21)) + + theme(legend.position = "bottom", legend.title = element_blank()) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y") + + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + + scale_y_continuous(labels=scales::percent) + + ylab("") + + xlab("Date") + + scale_fill_brewer(palette = "Dark2") + +ggsave("TimestampsMiddle.pdf", dpi = 500, w = 5.25, h = 5.25) + +TimeStampsFinal <- read.xlsx("/Users/sydney/Downloads/Standardizationtimestmaps.xlsx", sheetIndex = 3) +TimeStampsFinal <- pivot_longer(TimeStampsFinal, cols = c("Activity.A", "Activity.B", "Other")) +TimeStampsFinal <- TimeStampsFinal %>% mutate(name = case_when(name == "Activity.A" ~ "Activity A", + name == "Activity.B" ~ "Activity B", + name == "Other" ~ "Other")) +ggplot(TimeStampsFinal, aes(color = name, y = value, x = Date)) + + geom_hline(yintercept=0) + + geom_point(size = 4) + + theme_minimal() + + theme(text = element_text(size = 21)) + + theme(legend.position = "bottom", legend.title = element_blank()) + + scale_x_date(date_breaks = "1 month", date_labels = "%d/%m/%y") + + theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) + + scale_y_continuous(labels=scales::percent) + + ylab("") + + xlab("Date") + + scale_color_brewer(palette = "Dark2") + +ggsave("TimestampsRight.pdf", dpi = 500, w = 5.25, h = 5.25) + diff --git a/src/main/java/org/matsim/episim/EpisimRunner.java b/src/main/java/org/matsim/episim/EpisimRunner.java index e091c3744..ccd7801fb 100644 --- a/src/main/java/org/matsim/episim/EpisimRunner.java +++ b/src/main/java/org/matsim/episim/EpisimRunner.java @@ -24,9 +24,6 @@ import com.google.inject.Provider; import org.apache.commons.compress.archivers.*; import org.apache.commons.compress.archivers.zip.ZipArchiveEntry; -import org.apache.commons.compress.compressors.gzip.GzipCompressorOutputStream; -import org.apache.commons.csv.CSVFormat; -import org.apache.commons.csv.CSVPrinter; import org.apache.logging.log4j.LogManager; import org.apache.logging.log4j.Logger; import org.matsim.api.core.v01.events.Event; diff --git a/src/main/java/org/matsim/episim/EpisimUtils.java b/src/main/java/org/matsim/episim/EpisimUtils.java index 78cd0cb8f..8e6bda5f9 100644 --- a/src/main/java/org/matsim/episim/EpisimUtils.java +++ b/src/main/java/org/matsim/episim/EpisimUtils.java @@ -547,8 +547,9 @@ public static Map getOutdoorFractions2(File weatherCSV, File return outdoorFractions; } - public static Map getOutDoorFractionFromDateAndTemp2(File weatherCSV, File avgWeatherCSV, double rainThreshold, Double TmidSpring2020, Double TmidFall2020, Double TmidSpring, Double TmidFall, Double Trange, Double alpha, double maxOutdoorFraction) throws IOException { + public static Map getOutDoorFractionFromDateAndTemp2(File weatherCSV, File avgWeatherCSV, double rainThreshold, Double TmidSpring2020, Double TmidFall2020, Double TmidSpring, Double TmidFall, Double Trange, Double alpha, double maxOutdoorFraction) throws IOException { + // // 18.5 25 18.5 18.5 -> move to 15? Reader in = new FileReader(weatherCSV); Iterable records = CSVFormat.DEFAULT.withFirstRecordAsHeader().withCommentMarker('#').parse(in); DateTimeFormatter fmt = DateTimeFormatter.ofPattern("yyyy-MM-dd"); @@ -604,6 +605,68 @@ public static Map getOutDoorFractionFromDateAndTemp2(File wea } } +// System.exit(-1); + return outdoorFractions; + } + + public static Map getOutDoorFractionFromDateAndTemp2Fall2022Override(File weatherCSV, File avgWeatherCSV, double rainThreshold, Double TmidSpring2020, Double TmidFall2020, Double TmidSpring, Double TmidFall, Double TmidFall2022, Double Trange, Double alpha, double maxOutdoorFraction) throws IOException { + + Reader in = new FileReader(weatherCSV); + Iterable records = CSVFormat.DEFAULT.withFirstRecordAsHeader().withCommentMarker('#').parse(in); + DateTimeFormatter fmt = DateTimeFormatter.ofPattern("yyyy-MM-dd"); + + LocalDate lastDate = null; + final Map outdoorFractions = new TreeMap<>(); + for (CSVRecord record : records) { +// System.out.println( record ); + LocalDate date = LocalDate.parse(record.get("date"), fmt); + if (record.get("tmax").isEmpty() || record.get("prcp").isEmpty()) { +// System.out.println("Skipping day because tmax or prcp data is not available. Date: " + date.toString()); + continue; + } + + double tMax = Double.parseDouble(record.get("tmax")); + double prcp = Double.parseDouble(record.get("prcp")); + + if (date.isBefore(LocalDate.parse("2021-01-01"))) { + outdoorFractions.put(date, maxOutdoorFraction * getOutDoorFractionFromDateAndTemp(date, TmidSpring2020, TmidFall2020, Trange, tMax, prcp, rainThreshold, alpha)); + } else if (date.isAfter(LocalDate.parse("2022-01-01")) && date.isBefore(LocalDate.parse("2023-01-01"))) { + outdoorFractions.put(date, maxOutdoorFraction * getOutDoorFractionFromDateAndTemp(date, TmidSpring, TmidFall2022, Trange, tMax, prcp, rainThreshold, alpha)); + } else { + outdoorFractions.put(date, maxOutdoorFraction * getOutDoorFractionFromDateAndTemp(date, TmidSpring, TmidFall, Trange, tMax, prcp, rainThreshold, alpha)); + } + + lastDate = date; + } + + in = new FileReader(avgWeatherCSV); + records = CSVFormat.DEFAULT.withFirstRecordAsHeader().withCommentMarker('#').parse(in); + HashMap tmaxPerDay = new HashMap(); + HashMap prcpPerDay = new HashMap(); + + for (CSVRecord record : records) { + String monthDay = record.get("monthDay"); + double tMax = Double.parseDouble(record.get("tmax")); + double prcp = Double.parseDouble(record.get("prcp")); + tmaxPerDay.put(monthDay, tMax); + prcpPerDay.put(monthDay, prcp); + } + + for (int i = 1; i < 365*3; i++) { + LocalDate date = lastDate.plusDays(i); + int month = date.getMonth().getValue(); + int day = date.getDayOfMonth(); + String monthDay = month + "-" + day; + double tMax = tmaxPerDay.get(monthDay); + double prcp = prcpPerDay.get(monthDay); + if (date.isBefore(LocalDate.parse("2021-01-01"))) { + outdoorFractions.put(date, getOutDoorFractionFromDateAndTemp(date, TmidSpring2020, TmidFall2020, Trange, tMax, prcp, rainThreshold, alpha)); + } + else { + outdoorFractions.put(date, getOutDoorFractionFromDateAndTemp(date, TmidSpring, TmidFall, Trange, tMax, prcp, rainThreshold, alpha)); + } + } + // System.exit(-1); return outdoorFractions; } diff --git a/src/main/java/org/matsim/episim/InfectionEventHandler.java b/src/main/java/org/matsim/episim/InfectionEventHandler.java index 81f3641bd..66ff74195 100644 --- a/src/main/java/org/matsim/episim/InfectionEventHandler.java +++ b/src/main/java/org/matsim/episim/InfectionEventHandler.java @@ -1013,7 +1013,7 @@ void initImmunization(Path history) { if (handler.isContinueProcessingEvents()) { - throw new RuntimeException("Immunisation history is not long enough (only contains " + days.size() + "days)"); + throw new RuntimeException("Immunisation history is not long enough (only contains " + days.size() + " days)"); } } diff --git a/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEvents.java b/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEvents.java index 6dd17f5f1..c294ad59c 100644 --- a/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEvents.java +++ b/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEvents.java @@ -48,6 +48,7 @@ import java.nio.file.Path; import java.time.LocalDate; import java.util.*; + import java.util.stream.Collectors; /** @@ -60,7 +61,9 @@ public class HospitalNumbersFromEvents implements OutputAnalysis { - @CommandLine.Option(names = "--output", defaultValue = "/Users/jakob/git/matsim-episim/A_originalImmHist") + @CommandLine.Option(names = "--output", defaultValue = "/Users/jakob/git/matsim-episim/2023-10-27/events_hosp") +// @CommandLine.Option(names = "--output", defaultValue = "/Users/jakob/git/matsim-episim/2023-10-06/1/output/") +// @CommandLine.Option(names = "--output", defaultValue = "/Users/jakob/git/matsim-episim/A_originalImmHist") // @CommandLine.Option(names = "--output", defaultValue = "/Users/jakob/git/matsim-episim/B_startedFromImmHist") // @CommandLine.Option(names = "--output", defaultValue = "/Users/jakob/git/public-svn/matsim/scenarios/countries/de/episim/battery/jakob/2022-10-18/3-meas/analysis/") private Path output; @@ -93,97 +96,110 @@ public class HospitalNumbersFromEvents implements OutputAnalysis { // TODO: check age or strain based lags in literature // source: incidence wave vs. hospitalization wave in cologne/nrw (see https://docs.google.com/spreadsheets/d/1jmaerl27LKidD1uk3azdIL1LmvHuxazNQlhVo9xO1z8/edit?usp=sharing) - private static final Object2IntMap lagBetweenInfectionAndHospitalisation = new Object2IntAVLTreeMap<>( - Map.of(VirusStrain.SARS_CoV_2, 14, - VirusStrain.ALPHA, 14, - VirusStrain.DELTA, 14, - VirusStrain.OMICRON_BA1, 14, - VirusStrain.OMICRON_BA2, 14, - VirusStrain.OMICRON_BA5, 14, - VirusStrain.STRAIN_A, 14, - VirusStrain.STRAIN_B, 14 - )); - - // source: hospitalization wave vs. ICU wave in cologne/nrw (see https://docs.google.com/spreadsheets/d/1jmaerl27LKidD1uk3azdIL1LmvHuxazNQlhVo9xO1z8/edit?usp=sharing) - private static final Object2IntMap lagBetweenHospitalizationAndICU = new Object2IntAVLTreeMap<>( - Map.of(VirusStrain.SARS_CoV_2, 6, - VirusStrain.ALPHA, 6, - VirusStrain.DELTA, 6, - VirusStrain.OMICRON_BA1, 6, - VirusStrain.OMICRON_BA2, 6, - VirusStrain.OMICRON_BA5, 6, - VirusStrain.STRAIN_A, 6, - VirusStrain.STRAIN_B, 6 - )); + private static final Object2IntMap lagBetweenInfectionAndHospitalisation = setLagBetweenInfectionAndHospitalisation(); - // Austria study in https://docs.google.com/spreadsheets/d/1jmaerl27LKidD1uk3azdIL1LmvHuxazNQlhVo9xO1z8/edit#gid=0 - private static final Object2IntMap daysInHospitalGivenNoICU = new Object2IntAVLTreeMap<>( - Map.of(VirusStrain.SARS_CoV_2, 12, - VirusStrain.ALPHA, 12, - VirusStrain.DELTA, 12, - VirusStrain.OMICRON_BA1, 7, - VirusStrain.OMICRON_BA2, 7, - VirusStrain.OMICRON_BA5,7, - VirusStrain.STRAIN_A, 7, - VirusStrain.STRAIN_B, 7 - )); - - private static final Object2IntMap daysInICU = new Object2IntAVLTreeMap<>( - Map.of(VirusStrain.SARS_CoV_2, 15, // Debeka & Ireland studies - VirusStrain.ALPHA, 15, // Debeka & Ireland studies - VirusStrain.DELTA, 15, // this and following values come from nrw analysis on Tabellenblatt 5 - VirusStrain.OMICRON_BA1, 10, // TODO: Where does this number come from? - VirusStrain.OMICRON_BA2, 10, - VirusStrain.OMICRON_BA5,10, - VirusStrain.STRAIN_A, 10, - VirusStrain.STRAIN_B, 10 - )); + private static Object2IntAVLTreeMap setLagBetweenInfectionAndHospitalisation() { + Object2IntAVLTreeMap lagBetweenInfectionAndHospitalisation = new Object2IntAVLTreeMap<>(); - // ?? - private static final Object2IntMap daysInHospitalGivenICU = new Object2IntAVLTreeMap<>( - Map.of(VirusStrain.SARS_CoV_2, 60, // TODO: Where does this number come from? - VirusStrain.ALPHA, 60, - VirusStrain.DELTA, 60, - VirusStrain.OMICRON_BA1, 60, - VirusStrain.OMICRON_BA2, 60, - VirusStrain.OMICRON_BA5,60, - VirusStrain.STRAIN_A, 60, - VirusStrain.STRAIN_B, 60 - )); + for (VirusStrain strain : VirusStrain.values()) { + lagBetweenInfectionAndHospitalisation.put(strain, 14); + } + return lagBetweenInfectionAndHospitalisation; + } + // source: hospitalization wave vs. ICU wave in cologne/nrw (see https://docs.google.com/spreadsheets/d/1jmaerl27LKidD1uk3azdIL1LmvHuxazNQlhVo9xO1z8/edit?usp=sharing) + private static final Object2IntMap lagBetweenHospitalizationAndICU = setLagBetweenHospitalizationAndICU(); - private static final double beta = 1.2; + private static Object2IntAVLTreeMap setLagBetweenHospitalizationAndICU() { + Object2IntAVLTreeMap lagBetweenHospitalizationAndICU = new Object2IntAVLTreeMap<>(); - private static final double hospitalFactor = 0.3; // Based on "guess & check", accounts for unreported cases TODO: Potential follow-up + for (VirusStrain strain : VirusStrain.values()) { + lagBetweenHospitalizationAndICU.put(strain, 6); + } - private static final double factorWild = 1.0; + return lagBetweenHospitalizationAndICU; - private static final double factorAlpha = 1.0 * factorWild; + } - // delta: 2.3x more severe than alpha - Hospital admission and emergency care attendance risk for SARS-CoV-2 delta (B.1.617.2) compared with alpha (B.1.1.7) variants of concern: a cohort study - private static final double factorDelta = 1.2 * factorWild; //1.6 * factorWild; + // Austria study in https://docs.google.com/spreadsheets/d/1jmaerl27LKidD1uk3azdIL1LmvHuxazNQlhVo9xO1z8/edit#gid=0 + private static final Object2IntMap daysInHospitalGivenNoICU = setDaysInHospitalGivenNoICU(); - // omicron: approx 0.3x (intrinsic) severity of delta - Comparative analysis of the risks of hospitalisation and death associated with SARS-CoV-2 omicron (B.1.1.529) and delta (B.1.617.2) variants in England: a cohort study - private static final double factorOmicron = 0.45 * factorDelta; // reportedShareOmicron / reportedShareDelta -// private static final double factorOmicron = 0.6 * factorDelta;// reportedShareOmicron / reportedShareDelta + private static Object2IntAVLTreeMap setDaysInHospitalGivenNoICU() { + Object2IntAVLTreeMap daysInHospitalGivenNoICU = new Object2IntAVLTreeMap<>(); + for (VirusStrain strain : VirusStrain.values()) { + daysInHospitalGivenNoICU.put(strain, 7); + } - private static final double factorBA5 = 1.0 * factorOmicron; // old: 1.5 + daysInHospitalGivenNoICU.put(VirusStrain.SARS_CoV_2, 12); + daysInHospitalGivenNoICU.put(VirusStrain.ALPHA, 12); + daysInHospitalGivenNoICU.put(VirusStrain.DELTA, 12); - private static final double factorScen2 = factorBA5; // reportedShareOmicron / reportedShareDelta - private static final double factorScen3 = factorBA5 * 3;// reportedShareOmicron / reportedShareDelta + return daysInHospitalGivenNoICU; + } + private static final Object2IntMap daysInICU = setDaysInICU(); -// private static final double factorBA5 = 1.0 * factorOmicron; + private static Object2IntAVLTreeMap setDaysInICU() { + Object2IntAVLTreeMap daysInICU = new Object2IntAVLTreeMap<>(); + for (VirusStrain strain : VirusStrain.values()) { + daysInICU.put(strain, 10); + } + + daysInICU.put(VirusStrain.SARS_CoV_2, 15); + daysInICU.put(VirusStrain.ALPHA, 15); + daysInICU.put(VirusStrain.DELTA, 15); + return daysInICU; + } // ?? - private static final double factorWildAndAlphaICU = 1.; // TODO : Check literature for reasonable values - private static final double factorDeltaICU = 1.; - private static final double factorOmicronICU = 1.; - private static final double factorBA5ICU = 1.; + private static final Object2IntMap daysInHospitalGivenICU = setDaysInHospitalGivenICU(); + + private static Object2IntAVLTreeMap setDaysInHospitalGivenICU() { + Object2IntAVLTreeMap daysInHospitalGivenICU = new Object2IntAVLTreeMap<>(); + + for (VirusStrain strain : VirusStrain.values()) { + daysInHospitalGivenICU.put(strain, 60); + } + + return daysInHospitalGivenICU; + } + + + private static final double beta = 1.2; + + private static final double hospitalFactor = 0.3; // Based on "guess & check", accounts for unreported cases TODO: Potential follow-up + + private static final Map seriouslySickFactorModifier_BASE = Map.of( + VirusStrain.DELTA, 1.2, + VirusStrain.OMICRON_BA1, 0.45 + ); + + private static final Map seriouslySickFactorModifier_MILD = Map.of( + VirusStrain.DELTA, 1.2, + VirusStrain.OMICRON_BA1, 0.45, + VirusStrain.OMICRON_BA5, 1.2 + ); +// private static final Map seriouslySickFactorModifier_MILD = Map.of( +// VirusStrain.DELTA, 1.2, +// VirusStrain.OMICRON_BA1, 0.45, +// VirusStrain.OMICRON_BA5, 1.2 +// ); + + + private static final Map seriouslySickFactorModifier_SEVERE = Map.of( + VirusStrain.DELTA, 1.2, + VirusStrain.OMICRON_BA1, 0.45, + VirusStrain.OMICRON_BA5, 1.2, + VirusStrain.A_1, 1.5 + ); + + + // ICU: so far, we assume no difference between strains + private static final double factorICU = 1.; // TODO : Check literature for reasonable values public static void main(String[] args) { System.exit(new CommandLine(new HospitalNumbersFromEvents()).execute(args)); } @@ -225,8 +241,10 @@ public Integer call() throws Exception { // HospitalNumbersFromEventsPlotter.aggregateAndProducePlots(output, pathList); //TODO: move to other class -// HospitalNumbersFromEventsPlotter.aggregateAndProducePlots(output, pathList, "_Omicron", startDate, "Omicron"); -// HospitalNumbersFromEventsPlotter.aggregateAndProducePlots(output, pathList, "_Delta", startDate, "Delta"); +// HospitalNumbersFromEventsPlotter.aggregateAndProducePlots(output, pathList, "_Base", startDate, "Base"); +// HospitalNumbersFromEventsPlotter.aggregateAndProducePlots(output, pathList, "_Waning", startDate, "Waning"); +// HospitalNumbersFromEventsPlotter.aggregateAndProducePlots(output, pathList, "_Mild", startDate, "Mild"); +// HospitalNumbersFromEventsPlotter.aggregateAndProducePlots(output, pathList, "_Severe", startDate, "Severe"); // } @@ -266,12 +284,15 @@ private void calculateHospitalizationsAndWriteOutput(Path pathToScenario, Path t BufferedWriter bw = Files.newBufferedWriter(tsvPath); bw.write(AnalysisCommand.TSV.join(DAY, DATE,"measurement", "severity", "n")); // + "\thospNoImmunity\thospBaseImmunity\thospBoosted\tincNoImmunity\tincBaseImmunity\tincBoosted")); - ConfigHolder holderOmicron = configure(factorScen2, 1.0); // scenario 2 - ConfigHolder holderDelta = configure(factorScen3, 1.0); //scenario 3 + + ConfigHolder holderBase = configure(seriouslySickFactorModifier_BASE); +// ConfigHolder holderMild = configure(seriouslySickFactorModifier_MILD); + ConfigHolder holderSevere = configure(seriouslySickFactorModifier_SEVERE); List handlers = List.of( - new Handler("Omicron", population, holderOmicron), - new Handler("Delta", population, holderDelta) + new Handler("Base", population, holderBase), +// new Handler("Mild", population, holderMild), + new Handler("Severe", population, holderSevere) ); // feed the output events file to the handler, so that the hospitalizations may be calculated @@ -284,8 +305,21 @@ private void calculateHospitalizationsAndWriteOutput(Path pathToScenario, Path t // calculates the number of agents in the scenario's population (25% sample) who live in Cologne // this is used to normalize the hospitalization values double popSize = (int) population.getPersons().values().stream() - .filter(x -> x.getAttributes().getAttribute("district").equals(district)).count(); + .filter(x -> x.getAttributes().getAttribute("district").equals(district)).count(); + + + // calcualtes population in each age bin. + Int2LongAVLTreeMap popSizeByAge = new Int2LongAVLTreeMap(Collections.reverseOrder()); + long popAboveUpperBound = 0; + for (int lowerBound : handler.postProcessHospitalAdmissionsByAge.keySet()) { + long popAboveLowerBound = (long) population.getPersons().values().stream() + .filter(x -> x.getAttributes().getAttribute("district").equals(district)) + .filter(x -> (int) x.getAttributes().getAttribute("microm:modeled:age") >= lowerBound).count(); + long popInBin = popAboveLowerBound - popAboveUpperBound; + popAboveUpperBound = popAboveLowerBound; + popSizeByAge.put(lowerBound, popInBin); + } for (int day = 0; day < eventFiles.size(); day++) { LocalDate date = startDate.plusDays(day); @@ -300,7 +334,19 @@ private void calculateHospitalizationsAndWriteOutput(Path pathToScenario, Path t double occupancyIcu = handler.postProcessHospitalFilledBedsICU.getOrDefault(day, 0) * 100_000. / popSize; bw.newLine(); - bw.write(AnalysisCommand.TSV.join(day, date, HospitalNumbersFromEventsPlotter.INTAKES_HOSP, handler.name , intakesHosp)); + bw.write(AnalysisCommand.TSV.join(day, date, HospitalNumbersFromEventsPlotter.INTAKES_HOSP, handler.name, intakesHosp)); + + List ages = handler.postProcessHospitalAdmissionsByAge.keySet().stream().sorted().collect(Collectors.toList()); + + for (int i = 0; i < ages.size(); i++) { + int lowerBound = ages.get(i); + String lab = String.valueOf(lowerBound) + (i < ages.size() -1 ? "to" + (ages.get(i + 1) - 1) : "+"); + double incidenceForAgeBin = getWeeklyHospitalizations(handler.postProcessHospitalAdmissionsByAge.get(lowerBound), day) * 100_000. / popSizeByAge.get(lowerBound); + + bw.newLine(); + bw.write(AnalysisCommand.TSV.join(day, date, HospitalNumbersFromEventsPlotter.INTAKES_HOSP + "_" + lab, handler.name, incidenceForAgeBin)); + } + bw.newLine(); bw.write(AnalysisCommand.TSV.join(day, date, HospitalNumbersFromEventsPlotter.INTAKES_ICU, handler.name, intakesIcu)); bw.newLine(); @@ -341,6 +387,8 @@ public static final class Handler implements EpisimInfectionEventHandler, Episim final Int2IntSortedMap postProcessHospitalFilledBedsICU; private final AgeDependentDiseaseStatusTransitionModel transitionModel; + private final Int2ObjectAVLTreeMap postProcessHospitalAdmissionsByAge; + Handler(String name, Population population, ConfigHolder holder) { @@ -353,6 +401,15 @@ public static final class Handler implements EpisimInfectionEventHandler, Episim // key : iteration, value : admissions/filled beds this.postProcessHospitalAdmissions = new Int2IntAVLTreeMap(); + this.postProcessHospitalAdmissionsByAge = new Int2ObjectAVLTreeMap<>(Collections.reverseOrder()); + + Integer[] ageBins = {0, 18, 60, 80}; + + for (int ageLowerBound : ageBins) { + this.postProcessHospitalAdmissionsByAge.put(ageLowerBound, new Int2IntAVLTreeMap()); + } + + this.postProcessICUAdmissions = new Int2IntAVLTreeMap(); this.postProcessHospitalFilledBeds = new Int2IntAVLTreeMap(); this.postProcessHospitalFilledBedsICU = new Int2IntAVLTreeMap(); @@ -471,6 +528,14 @@ private void updateHospitalizationsPost(ImmunizablePerson person, VirusStrain st int inHospital = infectionIteration + lagBetweenInfectionAndHospitalisation.getInt(strain); postProcessHospitalAdmissions.mergeInt(inHospital, 1, Integer::sum); + + for (int lowerBound : postProcessHospitalAdmissionsByAge.keySet()) { + if (person.age >= lowerBound) { + postProcessHospitalAdmissionsByAge.get(lowerBound).merge(inHospital, 1, Integer::sum); + break; + } + } + if (goToICU(person, inHospital)) { // newly admitted to ICU @@ -685,7 +750,7 @@ public int getAge() { * necessary for post processing. * @param */ - private static ConfigHolder configure(double facA, double facAICU) { + private static ConfigHolder configure(Map seriouslySickFactorModifier) { Config config = ConfigUtils.createConfig(new EpisimConfigGroup()); @@ -695,26 +760,16 @@ private static ConfigHolder configure(double facA, double facAICU) { // configure strainConfig: add factorSeriouslySick for each strain VirusStrainConfigGroup strainConfig = ConfigUtils.addOrGetModule(config, VirusStrainConfigGroup.class); - strainConfig.getOrAddParams(VirusStrain.SARS_CoV_2).setFactorSeriouslySick(factorWild); - strainConfig.getOrAddParams(VirusStrain.SARS_CoV_2).setFactorCritical(factorWildAndAlphaICU); - strainConfig.getOrAddParams(VirusStrain.ALPHA).setFactorSeriouslySick(factorAlpha); - strainConfig.getOrAddParams(VirusStrain.ALPHA).setFactorCritical(factorAlpha); - - strainConfig.getOrAddParams(VirusStrain.DELTA).setFactorSeriouslySick(factorDelta); - strainConfig.getOrAddParams(VirusStrain.DELTA).setFactorCritical(factorDeltaICU); - - strainConfig.getOrAddParams(VirusStrain.OMICRON_BA1).setFactorSeriouslySick(factorOmicron); - strainConfig.getOrAddParams(VirusStrain.OMICRON_BA1).setFactorCritical(factorOmicronICU); - strainConfig.getOrAddParams(VirusStrain.OMICRON_BA2).setFactorSeriouslySick(factorOmicron); - strainConfig.getOrAddParams(VirusStrain.OMICRON_BA2).setFactorCritical(factorOmicronICU); - strainConfig.getOrAddParams(VirusStrain.OMICRON_BA5).setFactorSeriouslySick(factorBA5); - strainConfig.getOrAddParams(VirusStrain.OMICRON_BA5).setFactorCritical(factorBA5ICU); - - strainConfig.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySick(facA); - strainConfig.getOrAddParams(VirusStrain.STRAIN_A).setFactorCritical(facAICU); - - strainConfig.getOrAddParams(VirusStrain.STRAIN_B).setFactorSeriouslySick(facA); - strainConfig.getOrAddParams(VirusStrain.STRAIN_B).setFactorCritical(facAICU); + + for (VirusStrain strain : VirusStrain.values()) { + double seriouslySickFactorParent = 1.0; + if (strain.parent != null) { + seriouslySickFactorParent = strainConfig.getParams(strain.parent).getFactorSeriouslySick(); + } + strainConfig.getOrAddParams(strain).setFactorSeriouslySick(seriouslySickFactorParent * seriouslySickFactorModifier.getOrDefault(strain, 1.0)); + + strainConfig.getOrAddParams(strain).setFactorCritical(factorICU); + } // configure vaccinationConfig: set beta factor VaccinationConfigGroup vaccinationConfig = ConfigUtils.addOrGetModule(config, VaccinationConfigGroup.class); diff --git a/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEventsPlotter.java b/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEventsPlotter.java index 383c405a1..1e284b179 100644 --- a/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEventsPlotter.java +++ b/src/main/java/org/matsim/episim/analysis/HospitalNumbersFromEventsPlotter.java @@ -16,13 +16,18 @@ import tech.tablesaw.table.TableSliceGroup; import java.io.*; +import java.net.URL; import java.nio.charset.StandardCharsets; import java.nio.file.Files; import java.nio.file.Path; +import java.time.DayOfWeek; import java.time.LocalDate; import java.time.format.DateTimeFormatter; +import java.time.format.DateTimeFormatterBuilder; +import java.time.temporal.ChronoField; import java.time.temporal.ChronoUnit; import java.util.List; +import java.util.Locale; public class HospitalNumbersFromEventsPlotter { private static final String DATE = "date"; @@ -38,7 +43,7 @@ public class HospitalNumbersFromEventsPlotter { - static void aggregateAndProducePlots(Path output, List pathList, String outputAppendix, LocalDate startDate, String strainToPlot) throws IOException { + static void aggregateAndProducePlots(Path output, List pathList, String outputAppendix, LocalDate startDate, String scenarioToPlot) throws IOException { // read hospitalization tsv for all seeds and aggregate them! @@ -69,7 +74,7 @@ static void aggregateAndProducePlots(Path output, List pathList, String ou for (CSVRecord record : parser) { - if (!record.get("severity").equals(strainToPlot)) { + if (!record.get("severity").equals(scenarioToPlot)) { continue; } @@ -109,9 +114,45 @@ static void aggregateAndProducePlots(Path output, List pathList, String ou } } + // read rki COVID-SARI data (new) and add to tsv + + Int2DoubleMap covidSariHospIncidence = new Int2DoubleAVLTreeMap(); + URL url = new URL("https://raw.githubusercontent.com/robert-koch-institut/COVID-SARI-Hospitalisierungsinzidenz/main/COVID-SARI-Hospitalisierungsinzidenz.tsv"); + try (CSVParser parser = new CSVParser(new BufferedReader(new InputStreamReader(url.openStream())), + CSVFormat.DEFAULT.withDelimiter('\t').withFirstRecordAsHeader())) { + + for (CSVRecord record : parser) { + if (!record.get("agegroup").equals("00+")) { + continue; + } + + DateTimeFormatter formatter = new DateTimeFormatterBuilder() + .appendPattern("YYYY-'W'ww") + .parseDefaulting(ChronoField.DAY_OF_WEEK, DayOfWeek.THURSDAY.getValue()) + .toFormatter(Locale.GERMAN); + LocalDate date = LocalDate.parse(record.get("date"), formatter); + System.out.println(date); + + + int day = (int) startDate.until(date, ChronoUnit.DAYS); + System.out.println(day); + + double incidence; + try { + incidence = Double.parseDouble(record.get("sari_covid19_incidence")); + } catch (NumberFormatException e) { + incidence = Double.NaN; + } + + covidSariHospIncidence.put(day, incidence); + } + } + + // read rki data and add to tsv Int2DoubleMap rkiHospIncidence = new Int2DoubleAVLTreeMap(); Int2DoubleMap rkiHospIncidenceAdj = new Int2DoubleAVLTreeMap(); + // try (CSVParser parser = new CSVParser(Files.newBufferedReader(Path.of("../covid-sim/COVID-Hospitalization/Aktuell_Deutschland_adjustierte-COVID-19-Hospitalisierungen.csv")), try (CSVParser parser = new CSVParser(Files.newBufferedReader(Path.of("../covid-sim/src/assets/rki-deutschland-hospitalization.csv")), CSVFormat.DEFAULT.withDelimiter(',').withFirstRecordAsHeader())) { @@ -140,24 +181,6 @@ static void aggregateAndProducePlots(Path output, List pathList, String ou incidenceAdj = incidence / 3; } - -// if (date.isBefore(LocalDate.of(2020, 12, 10))) { -// incidenceAdj = incidence; -// } else if (date.isBefore(LocalDate.of(2021, 1, 11))) { -// incidenceAdj = 23. / 16. * incidence; -// } else if (date.isBefore(LocalDate.of(2021, 3, 22))) { -// incidenceAdj = 8. / 6. * incidence; -// } else if (date.isBefore(LocalDate.of(2021, 5, 3))) { -// incidenceAdj = 15./11. * incidence; -// } else if (date.isBefore(LocalDate.of(2021, 11, 8))) { -// incidenceAdj = incidence; -// } else if (date.isBefore(LocalDate.of(2021, 12, 6))) { -// incidenceAdj = 16. / 13. * incidence; -// } else if (date.isBefore(LocalDate.of(2022, 1, 24))) { -// incidenceAdj = incidence; -// } else { -// incidenceAdj = 11./14 * incidence; -// } rkiHospIncidenceAdj.put(day, incidenceAdj); } } @@ -360,6 +383,7 @@ static void aggregateAndProducePlots(Path output, List pathList, String ou date, intakeHosp.get(day), intakeIcu.get(day), + covidSariHospIncidence.get(day), occupancyHosp.get(day), occupancyIcu.get(day), rkiHospIncidence.get(day), @@ -391,67 +415,75 @@ static void aggregateAndProducePlots(Path output, List pathList, String ou groupings.append("model: intakeHosp"); } - // model: intakeIcu - for (Int2DoubleMap.Entry entry : intakeIcu.int2DoubleEntrySet()) { + for (Int2DoubleMap.Entry entry : covidSariHospIncidence.int2DoubleEntrySet()) { int day = entry.getIntKey(); records.append(day); recordsDate.append(startDate.plusDays(day)); - values.append(intakeIcu.getOrDefault(day, Double.NaN)); - groupings.append("model: intakeICU"); + values.append(covidSariHospIncidence.getOrDefault(day, Double.NaN)); + groupings.append("obs: covid-sari"); } + // model: intakeIcu +// for (Int2DoubleMap.Entry entry : intakeIcu.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// values.append(intakeIcu.getOrDefault(day, Double.NaN)); +// groupings.append("model: intakeICU"); +// } - for (Int2DoubleMap.Entry entry : rkiHospIncidence.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - final double value = rkiHospIncidence.getOrDefault(day, Double.NaN); - if (Double.isNaN(value)) { - values.appendMissing(); - } else { - values.append(value); - } - groupings.append("reported: intakeHosp (rki, nrw adjusted) WITH Covid"); - } - - for (Int2DoubleMap.Entry entry : rkiHospIncidenceAdj.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - final double value = rkiHospIncidenceAdj.getOrDefault(day, Double.NaN); - if (Double.isNaN(value)) { - values.appendMissing(); - } else { - values.append(value); - } - groupings.append("reported: intakeHosp (rki, nrw adjusted) FROM Covid"); - } - - for (Int2DoubleMap.Entry entry : hospIncidenceKoeln.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - final double value = hospIncidenceKoeln.getOrDefault(day, Double.NaN); - if (Double.isNaN(value)) { - values.appendMissing(); - } else { - values.append(value); - } - groupings.append("reported: intakeHosp (köln)"); - } - - for (Int2DoubleMap.Entry entry : reportedIcuIncidence.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - final double value = reportedIcuIncidence.getOrDefault(day, Double.NaN); - if (Double.isNaN(value)) { - values.appendMissing(); - } else { - values.append(value); - } - groupings.append("reported: intakeIcu (divi, nrw)"); - } +// +// for (Int2DoubleMap.Entry entry : rkiHospIncidence.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// final double value = rkiHospIncidence.getOrDefault(day, Double.NaN); +// if (Double.isNaN(value)) { +// values.appendMissing(); +// } else { +// values.append(value); +// } +// groupings.append("reported: intakeHosp (rki, nrw adjusted) WITH Covid"); +// } +// +// for (Int2DoubleMap.Entry entry : rkiHospIncidenceAdj.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// final double value = rkiHospIncidenceAdj.getOrDefault(day, Double.NaN); +// if (Double.isNaN(value)) { +// values.appendMissing(); +// } else { +// values.append(value); +// } +// groupings.append("reported: intakeHosp (rki, nrw adjusted) FROM Covid"); +// } +// +// for (Int2DoubleMap.Entry entry : hospIncidenceKoeln.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// final double value = hospIncidenceKoeln.getOrDefault(day, Double.NaN); +// if (Double.isNaN(value)) { +// values.appendMissing(); +// } else { +// values.append(value); +// } +// groupings.append("reported: intakeHosp (köln)"); +// } +// +// for (Int2DoubleMap.Entry entry : reportedIcuIncidence.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// final double value = reportedIcuIncidence.getOrDefault(day, Double.NaN); +// if (Double.isNaN(value)) { +// values.appendMissing(); +// } else { +// values.append(value); +// } +// groupings.append("reported: intakeIcu (divi, nrw)"); +// } producePlot(recordsDate, values, groupings, "", "7-Tage Hospitalisierungsinzidenz", "HospIncidence" + outputAppendix + ".html", output); @@ -459,109 +491,109 @@ static void aggregateAndProducePlots(Path output, List pathList, String ou // PLOT 2: People taking up beds in hospital (regular and ICU) - { - IntColumn records = IntColumn.create("day"); - DateColumn recordsDate = DateColumn.create("date"); - DoubleColumn values = DoubleColumn.create("hospitalizations"); - StringColumn groupings = StringColumn.create("scenario"); - - - for (Int2DoubleMap.Entry entry : occupancyHosp.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - - values.append(occupancyHosp.get(day)); - groupings.append("generalBeds"); - } - - for (Int2DoubleMap.Entry entry : occupancyIcu.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - - values.append(occupancyIcu.get(day)); - groupings.append("ICUBeds"); - } - - - for (Int2DoubleMap.Entry entry : reportedBeds.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - - values.append(entry.getDoubleValue()); - groupings.append("Reported: General Beds"); - - } - - for (Int2DoubleMap.Entry entry : reportedBedsAdj.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - - values.append(entry.getDoubleValue()); - groupings.append("Reported: General Beds (SARI)"); - } - - - for (Int2DoubleMap.Entry entry : reportedBedsICU.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - - values.append(entry.getDoubleValue()); - groupings.append("Reported: ICU Beds"); - - } - - -// for (Int2DoubleMap.Entry entry : reportedBedsNrw.int2DoubleEntrySet()) { +// { +// IntColumn records = IntColumn.create("day"); +// DateColumn recordsDate = DateColumn.create("date"); +// DoubleColumn values = DoubleColumn.create("hospitalizations"); +// StringColumn groupings = StringColumn.create("scenario"); +// +// +// for (Int2DoubleMap.Entry entry : occupancyHosp.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// +// values.append(occupancyHosp.get(day)); +// groupings.append("generalBeds"); +// } +// +// for (Int2DoubleMap.Entry entry : occupancyIcu.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// +// values.append(occupancyIcu.get(day)); +// groupings.append("ICUBeds"); +// } +// +// +// for (Int2DoubleMap.Entry entry : reportedBeds.int2DoubleEntrySet()) { // int day = entry.getIntKey(); // records.append(day); // recordsDate.append(startDate.plusDays(day)); // // values.append(entry.getDoubleValue()); -// groupings.append("Reported: General Beds (NRW)"); +// groupings.append("Reported: General Beds"); // // } - - -// for (Int2DoubleMap.Entry entry : reportedBedsIcuNrw.int2DoubleEntrySet()) { +// +// for (Int2DoubleMap.Entry entry : reportedBedsAdj.int2DoubleEntrySet()) { // int day = entry.getIntKey(); // records.append(day); // recordsDate.append(startDate.plusDays(day)); // // values.append(entry.getDoubleValue()); -// groupings.append("Reported: ICU Beds (NRW)"); +// groupings.append("Reported: General Beds (SARI)"); +// } +// +// +// for (Int2DoubleMap.Entry entry : reportedBedsICU.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// +// values.append(entry.getDoubleValue()); +// groupings.append("Reported: ICU Beds"); // // } - - for (Int2DoubleMap.Entry entry : reportedBedsNrw2.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - - values.append(entry.getDoubleValue()); - groupings.append("Reported: General Beds (NRW2)"); - - } - - for (Int2DoubleMap.Entry entry : reportedBedsIcuNrw2.int2DoubleEntrySet()) { - int day = entry.getIntKey(); - records.append(day); - recordsDate.append(startDate.plusDays(day)); - - values.append(entry.getDoubleValue()); - groupings.append("Reported: ICU Beds (NRW2)"); - - } - - - - // Make plot - producePlot(recordsDate, values, groupings, "Filled Beds", "Beds Filled / 100k Population", "FilledBeds" + outputAppendix + ".html", output); - } +// +// +//// for (Int2DoubleMap.Entry entry : reportedBedsNrw.int2DoubleEntrySet()) { +//// int day = entry.getIntKey(); +//// records.append(day); +//// recordsDate.append(startDate.plusDays(day)); +//// +//// values.append(entry.getDoubleValue()); +//// groupings.append("Reported: General Beds (NRW)"); +//// +//// } +// +// +//// for (Int2DoubleMap.Entry entry : reportedBedsIcuNrw.int2DoubleEntrySet()) { +//// int day = entry.getIntKey(); +//// records.append(day); +//// recordsDate.append(startDate.plusDays(day)); +//// +//// values.append(entry.getDoubleValue()); +//// groupings.append("Reported: ICU Beds (NRW)"); +//// +//// } +// +// for (Int2DoubleMap.Entry entry : reportedBedsNrw2.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// +// values.append(entry.getDoubleValue()); +// groupings.append("Reported: General Beds (NRW2)"); +// +// } +// +// for (Int2DoubleMap.Entry entry : reportedBedsIcuNrw2.int2DoubleEntrySet()) { +// int day = entry.getIntKey(); +// records.append(day); +// recordsDate.append(startDate.plusDays(day)); +// +// values.append(entry.getDoubleValue()); +// groupings.append("Reported: ICU Beds (NRW2)"); +// +// } +// +// +// +// // Make plot +// producePlot(recordsDate, values, groupings, "Filled Beds", "Beds Filled / 100k Population", "FilledBeds" + outputAppendix + ".html", output); +// } } diff --git a/src/main/java/org/matsim/episim/model/AntibodyModel.java b/src/main/java/org/matsim/episim/model/AntibodyModel.java index bbda67d6f..b84f52bf4 100644 --- a/src/main/java/org/matsim/episim/model/AntibodyModel.java +++ b/src/main/java/org/matsim/episim/model/AntibodyModel.java @@ -46,12 +46,14 @@ static AntibodyModel.Config newConfig(Map> initialAntibodies; final Map> antibodyRefreshFactors; private double immuneReponseSigma = 0.; + double hlMultiForInfected = 1.0; + public Config() { //initial antibodies @@ -203,6 +205,13 @@ public Config(Map> initialAntibodies, Ma } + public Config(Map> initialAntibodies, Map> antibodyRefreshFactors, double hlMultiForInfected) { + this.initialAntibodies = initialAntibodies; + this.antibodyRefreshFactors = antibodyRefreshFactors; + this.hlMultiForInfected = hlMultiForInfected; + + } + public double getImmuneReponseSigma() { return immuneReponseSigma; } diff --git a/src/main/java/org/matsim/episim/model/DefaultAntibodyModel.java b/src/main/java/org/matsim/episim/model/DefaultAntibodyModel.java index b040e796f..24880ad6a 100644 --- a/src/main/java/org/matsim/episim/model/DefaultAntibodyModel.java +++ b/src/main/java/org/matsim/episim/model/DefaultAntibodyModel.java @@ -13,6 +13,9 @@ public class DefaultAntibodyModel implements AntibodyModel { public static final double HALF_LIFE_DAYS = 60; // todo: would 40 work better? + // optimistic: 46 + // + private final AntibodyModel.Config antibodyConfig; private final SplittableRandom localRnd; @@ -94,10 +97,19 @@ public void updateAntibodies(EpisimPerson person, int day) { } } + double halflifeDays = HALF_LIFE_DAYS; + + if ((person.getNumInfections() > 0 && person.getNumVaccinations() > 0) + || person.getNumInfections() > 4 + || person.getNumVaccinations() > 3) { + + halflifeDays *= antibodyConfig.hlMultiForInfected; // 1, 2, 5 + } + // if no immunity event: exponential decay, day by day: for (VirusStrain strain : VirusStrain.values()) { double oldAntibodyLevel = person.getAntibodies(strain); - person.setAntibodies(strain, oldAntibodyLevel * Math.pow(0.5, 1 / HALF_LIFE_DAYS)); + person.setAntibodies(strain, oldAntibodyLevel * Math.pow(0.5, 1 / halflifeDays)); } } diff --git a/src/main/java/org/matsim/episim/model/InfectionModelWithAntibodies.java b/src/main/java/org/matsim/episim/model/InfectionModelWithAntibodies.java index d70712dfb..3b2ffbccc 100644 --- a/src/main/java/org/matsim/episim/model/InfectionModelWithAntibodies.java +++ b/src/main/java/org/matsim/episim/model/InfectionModelWithAntibodies.java @@ -108,14 +108,18 @@ public double calcInfectionProbability(EpisimPerson target, EpisimPerson infecto // igaFactor = Math.exp( - target.daysSinceInfection(lastInfectionWithStrain, iteration) / 120.0); igaFactor = 1.0 / (1.0 + Math.exp(-2.0 * (1.0 - target.daysSinceInfection(lastInfectionWithStrain, iteration) / igaTimePeriod))); - } - + } + ArrayList strainsLineA = new ArrayList(); strainsLineA.add(VirusStrain.OMICRON_BA1); strainsLineA.add(VirusStrain.OMICRON_BA2); strainsLineA.add(VirusStrain.OMICRON_BA5); strainsLineA.add(VirusStrain.STRAIN_A); strainsLineA.add(VirusStrain.STRAIN_B); + strainsLineA.add(VirusStrain.BQ); + strainsLineA.add(VirusStrain.XBB_15); + strainsLineA.add(VirusStrain.XBB_19); + strainsLineA.add(VirusStrain.EG); ArrayList strainsLineB = new ArrayList(); strainsLineB.add(VirusStrain.OMICRON_BA1); @@ -123,9 +127,13 @@ public double calcInfectionProbability(EpisimPerson target, EpisimPerson infecto strainsLineB.add(VirusStrain.OMICRON_BA5); strainsLineB.add(VirusStrain.STRAIN_A); strainsLineB.add(VirusStrain.STRAIN_B); - + strainsLineB.add(VirusStrain.BQ); + strainsLineB.add(VirusStrain.XBB_15); + strainsLineB.add(VirusStrain.XBB_19); + strainsLineB.add(VirusStrain.EG); + if (vaccinationConfig.getUseIgA()) { - + for (VirusStrain str : VirusStrain.values()) { if (str.toString().startsWith("A_")) strainsLineA.add(str); @@ -135,7 +143,7 @@ public double calcInfectionProbability(EpisimPerson target, EpisimPerson infecto strainsLineB.add(str); } } - + // if (vaccinationConfig.getUseIgA()) { if(strainsLineA.contains(infector.getVirusStrain())){ @@ -154,7 +162,7 @@ public double calcInfectionProbability(EpisimPerson target, EpisimPerson infecto igaFactor = Math.max(fac, igaFactor); } } - + if(strainsLineB.contains(infector.getVirusStrain())){ int lastInfectionWithStrain = 0; boolean targetHadStrain = false; @@ -171,7 +179,7 @@ public double calcInfectionProbability(EpisimPerson target, EpisimPerson infecto igaFactor = Math.max(fac, igaFactor); } } - + // if(strainsLineA.contains(infector.getVirusStrain()) || strainsLineB.contains(infector.getVirusStrain())){ // int lastInfectionWithStrain = 0; // boolean targetHadStrain = false; @@ -188,9 +196,9 @@ public double calcInfectionProbability(EpisimPerson target, EpisimPerson infecto // igaFactor = Math.max(fac, igaFactor); // } // } - + susceptibility = susceptibility * (1.0 - igaFactor); -// } +// } } diff --git a/src/main/java/org/matsim/episim/model/VaccinationType.java b/src/main/java/org/matsim/episim/model/VaccinationType.java index 57f9be945..8022ac0cc 100644 --- a/src/main/java/org/matsim/episim/model/VaccinationType.java +++ b/src/main/java/org/matsim/episim/model/VaccinationType.java @@ -10,6 +10,7 @@ public enum VaccinationType implements ImmunityEvent { vector, ba1Update, ba5Update, + xbbUpdate, /** * Not a real vaccination, but used to describe the profile for persons that have been infected and gained a natural immunity. diff --git a/src/main/java/org/matsim/episim/model/VirusStrain.java b/src/main/java/org/matsim/episim/model/VirusStrain.java index cc278b5c8..515230772 100644 --- a/src/main/java/org/matsim/episim/model/VirusStrain.java +++ b/src/main/java/org/matsim/episim/model/VirusStrain.java @@ -8,117 +8,133 @@ public enum VirusStrain implements ImmunityEvent { /** * This describes the base virus strain. */ - SARS_CoV_2, + SARS_CoV_2(null), /** * More "infectious" variant B.1.1.7 that has been prevalent in the UK, starting during end of 2020. * Also known as VOC-202012/01. */ // B117, - ALPHA, + ALPHA(SARS_CoV_2), /** * South-african variant also known as auch 501Y.V2. */ - B1351, + B1351(SARS_CoV_2), //todo? /** * unknown mutation */ // MUTB, - DELTA, + DELTA(ALPHA), /** * VoC B.1.1.529, first reported to WHO from South Africa on 24 November 2021 */ // OMICRON, - OMICRON_BA1, - - OMICRON_BA2, - - OMICRON_BA5, - - STRAIN_A, - - STRAIN_B, - - A_1, - - A_2, - - A_3, - - A_4, - - A_5, - - A_6, - - A_7, - - A_8, - - A_9, - - A_10, - - A_11, - - A_12, - - A_13, - - A_14, - - A_15, - - A_16, - - A_17, - - A_18, - - A_19, - - A_20, - - B_1, - - B_2, - - B_3, - - B_4, - - B_5, - - B_6, - - B_7, - - B_8, - - B_9, - - B_10, - - B_11, - - B_12, - - B_13, - - B_14, - - B_15, - - B_16, - - B_17, - - B_18, - - B_19, - - B_20 + OMICRON_BA1(DELTA), + + OMICRON_BA2(OMICRON_BA1), + + OMICRON_BA5(OMICRON_BA2), + + XBB_15(OMICRON_BA2), + + XBB_19(OMICRON_BA2), + + BQ(OMICRON_BA5), + + EG(XBB_19), + + STRAIN_A(OMICRON_BA5), + + STRAIN_B(OMICRON_BA5), + + A_1(EG), + + A_2(A_1), + + A_3(A_2), + + A_4(A_3), + + A_5(A_4), + + A_6(A_5), + + A_7(A_6), + + A_8(A_7), + + A_9(A_8), + + A_10(A_9), + + A_11(A_10), + + A_12(A_11), + + A_13(A_12), + + A_14(A_13), + + A_15(A_14), + + A_16(A_15), + + A_17(A_16), + + A_18(A_17), + + A_19(A_18), + + A_20(A_19), + + B_1(null), + + B_2(null), + + B_3(null), + + B_4(null), + + B_5(null), + + B_6(null), + + B_7(null), + + B_8(null), + + B_9(null), + + B_10(null), + + B_11(null), + + B_12(null), + + B_13(null), + + B_14(null), + + B_15(null), + + B_16(null), + + B_17(null), + + B_18(null), + + B_19(null), + + B_20(null); + + public final VirusStrain parent; + + VirusStrain(VirusStrain parent) { + this.parent = parent; + } + + } diff --git a/src/main/java/org/matsim/episim/model/progression/AgeDependentDiseaseStatusTransitionModel.java b/src/main/java/org/matsim/episim/model/progression/AgeDependentDiseaseStatusTransitionModel.java index 6f72b99c0..8ce626d16 100644 --- a/src/main/java/org/matsim/episim/model/progression/AgeDependentDiseaseStatusTransitionModel.java +++ b/src/main/java/org/matsim/episim/model/progression/AgeDependentDiseaseStatusTransitionModel.java @@ -74,12 +74,14 @@ protected double getProbaOfTransitioningToShowingSymptoms (EpisimPerson person) @Override public double getProbaOfTransitioningToSeriouslySick(Immunizable person) { + // https://docs.google.com/spreadsheets/d/1jmaerl27LKidD1uk3azdIL1LmvHuxazNQlhVo9xO1z8/edit#gid=802030488 + double proba; int age = person.getAge(); - // source 3 + // source 3: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8353925/ if (age < 5) { proba = 4.0 / 100; } else if (age < 15) { diff --git a/src/main/java/org/matsim/episim/model/progression/AntibodyDependentTransitionModel.java b/src/main/java/org/matsim/episim/model/progression/AntibodyDependentTransitionModel.java index 63dc0dc36..ead6bc7fc 100644 --- a/src/main/java/org/matsim/episim/model/progression/AntibodyDependentTransitionModel.java +++ b/src/main/java/org/matsim/episim/model/progression/AntibodyDependentTransitionModel.java @@ -128,8 +128,17 @@ public double getSeriouslySickFactor(Immunizable person, VaccinationConfigGroup abNoWaning *= 4; } // b) if strain is omicron, an additional factor of 3.7 is applied - if (strain.equals(VirusStrain.OMICRON_BA1) || strain.equals(VirusStrain.OMICRON_BA2) || strain.equals(VirusStrain.OMICRON_BA5) || strain.equals(VirusStrain.STRAIN_A) || strain.equals(VirusStrain.STRAIN_B) - || strain.toString().startsWith("A_") || strain.toString().startsWith("B_")) { + if (strain.equals(VirusStrain.OMICRON_BA1) || + strain.equals(VirusStrain.OMICRON_BA2) || + strain.equals(VirusStrain.OMICRON_BA5) || + strain.equals(VirusStrain.BQ) || + strain.equals(VirusStrain.EG) || + strain.equals(VirusStrain.XBB_15) || + strain.equals(VirusStrain.XBB_19) || + strain.equals(VirusStrain.STRAIN_A) || + strain.equals(VirusStrain.STRAIN_B) || + strain.toString().startsWith("A_") || + strain.toString().startsWith("B_")) { abNoWaning *= 3.7; } diff --git a/src/main/java/org/matsim/episim/model/testing/DefaultTestingModel.java b/src/main/java/org/matsim/episim/model/testing/DefaultTestingModel.java index 2036656dc..8bbee11cc 100644 --- a/src/main/java/org/matsim/episim/model/testing/DefaultTestingModel.java +++ b/src/main/java/org/matsim/episim/model/testing/DefaultTestingModel.java @@ -203,13 +203,17 @@ protected boolean testAndQuarantine(EpisimPerson person, int day, TestingConfigG double rate = params.getFalseNegativeRate(); // TODO: configurable - if (params.getType().equals(TestType.RAPID_TEST) && (person.getVirusStrain() == VirusStrain.OMICRON_BA1 || + if (params.getType().equals(TestType.RAPID_TEST) && ( + person.getVirusStrain() == VirusStrain.OMICRON_BA1 || person.getVirusStrain() == VirusStrain.OMICRON_BA2 || person.getVirusStrain() == VirusStrain.OMICRON_BA5 || + person.getVirusStrain() == VirusStrain.BQ || + person.getVirusStrain() == VirusStrain.XBB_15 || person.getVirusStrain() == VirusStrain.STRAIN_A || person.getVirusStrain() == VirusStrain.STRAIN_B || person.getVirusStrain().toString().startsWith("A_") || - person.getVirusStrain().toString().startsWith("B_"))) { + person.getVirusStrain().toString().startsWith("B_")) + ) { rate = 0.5; } diff --git a/src/main/java/org/matsim/episim/model/vaccination/VaccinationStrategyBMBF202310XX.java b/src/main/java/org/matsim/episim/model/vaccination/VaccinationStrategyBMBF202310XX.java new file mode 100644 index 000000000..d57904dc1 --- /dev/null +++ b/src/main/java/org/matsim/episim/model/vaccination/VaccinationStrategyBMBF202310XX.java @@ -0,0 +1,151 @@ +package org.matsim.episim.model.vaccination; + +import com.google.inject.Inject; +import it.unimi.dsi.fastutil.ints.Int2DoubleMap; +import org.matsim.api.core.v01.Id; +import org.matsim.api.core.v01.population.Person; +import org.matsim.episim.EpisimPerson; +import org.matsim.episim.EpisimUtils; +import org.matsim.episim.model.VaccinationType; + +import java.time.LocalDate; +import java.util.*; +import java.util.stream.Collectors; + +/** + * Update vaccination campaign. + */ +public class VaccinationStrategyBMBF202310XX implements VaccinationModel { + + private final SplittableRandom rnd; + private final Config config; + + @Inject + public VaccinationStrategyBMBF202310XX(SplittableRandom rnd, Config config) { + this.rnd = rnd; + this.config = config; + } + + @Override + public void handleVaccination(Map, EpisimPerson> persons, LocalDate date, int iteration, double now) { + + if (date.isAfter(config.start) && date.isBefore(config.start.plusDays(config.campaignDuration))) { + + if (config.complianceByAge.values().doubleStream() + .reduce(0, Double::sum) > 0) { + + List filteredPersons = new ArrayList<>(persons.values()); + + List candidates = filteredPersons.stream() + .filter(EpisimPerson::isVaccinable) // todo: what determines who is vaccinable? + .filter(p -> p.getDiseaseStatus() == EpisimPerson.DiseaseStatus.susceptible) + .filter(p -> p.getNumVaccinations() >= config.vaccinationPool.vaxCnt) + .filter(p -> p.getNumVaccinations() == 0 || p.daysSinceVaccination(p.getNumVaccinations() - 1, iteration) > config.minDaysAfterVaccination) // only people who've had their last vaccination more than 90 days ago + .filter(p -> p.getNumInfections() == 0 || p.daysSinceInfection(p.getNumInfections() - 1, iteration) > config.minDaysAfterInfection) // only people who've had their last vaccination more than 90 days ago + .collect(Collectors.toList()); + + int num0 = (int) filteredPersons.stream().filter(p -> (p.getAge() >= 0 && p.getAge() < 12)).count(); + int num12 = (int) filteredPersons.stream().filter(p -> (p.getAge() >= 12 && p.getAge() < 18)).count(); + int num18 = (int) filteredPersons.stream().filter(p -> (p.getAge() >= 18 && p.getAge() < 60)).count(); + int num60 = (int) filteredPersons.stream().filter(p -> p.getAge() >= 60).count(); + + double numVaccines = num0 * config.complianceByAge.get(0) + + num12 * config.complianceByAge.get(12) + + num18 * config.complianceByAge.get(18) + + num60 * config.complianceByAge.get(60); + + int vaccinationsLeft = (int) (numVaccines / config.campaignDuration); + + + List people0 = candidates.stream().filter(p -> (p.getAge() >= 0 && p.getAge() < 12)).collect(Collectors.toList()); + List people12 = candidates.stream().filter(p -> (p.getAge() >= 12 && p.getAge() < 18)).collect(Collectors.toList()); + List people18 = candidates.stream().filter(p -> (p.getAge() >= 18 && p.getAge() < 60)).collect(Collectors.toList()); + List people60 = candidates.stream().filter(p -> p.getAge() >= 60).collect(Collectors.toList()); + + final List[] perAge = new List[4]; + perAge[0] = people0; + perAge[1] = people12; + perAge[2] = people18; + perAge[3] = people60; + + int ageIndex = 3; + while (ageIndex >= 0 && vaccinationsLeft > 0) { + + List candidatesForAge = perAge[ageIndex]; + + // list is shuffled to avoid eventual bias + if (candidatesForAge.size() > vaccinationsLeft) + Collections.shuffle(perAge[ageIndex], new Random(EpisimUtils.getSeed(rnd))); + + + int vaccinesForDayAndAgeGroup = Math.min(candidatesForAge.size(), vaccinationsLeft); + for (int i = 0; i < vaccinesForDayAndAgeGroup; i++) { + EpisimPerson person = candidatesForAge.get(i); + vaccinate(person, iteration, config.vaccinationType); + vaccinationsLeft--; + } + + ageIndex--; + + } + } + + } + + } + + + public static class Config { + + /** + * Start of vaccination campaign. + */ + private final LocalDate start; + /** + * Duration of vaccination campaign. + */ + private final int campaignDuration; + + private final VaccinationType vaccinationType; + + private final Int2DoubleMap complianceByAge; + + private final VaccinationPool vaccinationPool; + + private final int minDaysAfterVaccination; + + private final int minDaysAfterInfection; + + public enum VaccinationPool { + + unvaccinated(0), + + vaccinated(1), + + boostered(2); + + private final int vaxCnt; + + VaccinationPool(int vaxCnt) { + this.vaxCnt = vaxCnt; + } + + } + + + public Config(LocalDate start, int campaignDuration, VaccinationType vaccinationType, + Int2DoubleMap complianceByAge, VaccinationPool vaccinationPool, + int minDaysAfterInfection, int minDaysAfterVaccination) { + + this.start = start; + this.campaignDuration = campaignDuration; + this.vaccinationType = vaccinationType; + this.complianceByAge = complianceByAge; + this.vaccinationPool = vaccinationPool; + this.minDaysAfterInfection = minDaysAfterInfection; + this.minDaysAfterVaccination = minDaysAfterVaccination; + + } + } + +} diff --git a/src/main/java/org/matsim/run/CreateBatteryForCluster.java b/src/main/java/org/matsim/run/CreateBatteryForCluster.java index 0f59fc630..0a9adc16c 100644 --- a/src/main/java/org/matsim/run/CreateBatteryForCluster.java +++ b/src/main/java/org/matsim/run/CreateBatteryForCluster.java @@ -84,10 +84,10 @@ public class CreateBatteryForCluster implements Callable { @CommandLine.Option(names = "--jvm-opts", description = "Additional options for JVM", defaultValue = "-Xms82G -Xmx82G -XX:+UseParallelGC") private String jvmOpts; - @CommandLine.Option(names = "--setup", defaultValue = "org.matsim.run.batch.CologneBMBF202212XX_bq1") + @CommandLine.Option(names = "--setup", defaultValue = "org.matsim.run.batch.StarterBatchCologne") private Class> setup; - @CommandLine.Option(names = "--params", defaultValue = "org.matsim.run.batch.CologneBMBF202212XX_bq1$Params") + @CommandLine.Option(names = "--params", defaultValue = "org.matsim.run.batch.StarterBatchCologne$Params") private Class params; @SuppressWarnings("rawtypes") diff --git a/src/main/java/org/matsim/run/batch/CologneBMBF202212XX_soup.java b/src/main/java/org/matsim/run/batch/CologneBMBF20221202_soup.java similarity index 99% rename from src/main/java/org/matsim/run/batch/CologneBMBF202212XX_soup.java rename to src/main/java/org/matsim/run/batch/CologneBMBF20221202_soup.java index 0ce1f29ec..3cc1f048b 100644 --- a/src/main/java/org/matsim/run/batch/CologneBMBF202212XX_soup.java +++ b/src/main/java/org/matsim/run/batch/CologneBMBF20221202_soup.java @@ -30,7 +30,7 @@ /** * Batch for Bmbf runs */ -public class CologneBMBF202212XX_soup implements BatchRun { +public class CologneBMBF20221202_soup implements BatchRun { boolean DEBUG_MODE = true; int runCount = 0; @@ -652,7 +652,7 @@ public static final class Params { public static void main(String[] args) { String[] args2 = { - RunParallel.OPTION_SETUP, CologneBMBF202212XX_soup.class.getName(), + RunParallel.OPTION_SETUP, CologneBMBF20221202_soup.class.getName(), RunParallel.OPTION_PARAMS, Params.class.getName(), RunParallel.OPTION_TASKS, Integer.toString(1), RunParallel.OPTION_ITERATIONS, Integer.toString(1000), diff --git a/src/main/java/org/matsim/run/batch/CologneBMBF202212XX_bq1.java b/src/main/java/org/matsim/run/batch/CologneBMBF2023.java similarity index 54% rename from src/main/java/org/matsim/run/batch/CologneBMBF202212XX_bq1.java rename to src/main/java/org/matsim/run/batch/CologneBMBF2023.java index a29b25077..a8263d8cd 100644 --- a/src/main/java/org/matsim/run/batch/CologneBMBF202212XX_bq1.java +++ b/src/main/java/org/matsim/run/batch/CologneBMBF2023.java @@ -26,7 +26,7 @@ /** * Batch for Bmbf runs */ -public class CologneBMBF202212XX_bq1 implements BatchRun { +public class CologneBMBF2023 implements BatchRun { boolean DEBUG_MODE = false; @@ -60,36 +60,36 @@ protected void configure() { Map startDateToVaccination = new HashMap<>(); startDateToVaccination.put(start, vaccinationType); - if (params != null) { - if (params.vacCamp.equals("base")) { // + - - } else if(params.vacCamp.equals("ph1_90")){ - minDaysAfterInfection = 90; - minDaysAfterVaccination = 90; - - emergencyDate = restrictionDatePhase1; - } else if(params.vacCamp.equals("ph1_90vax180")){ - minDaysAfterInfection = 90; - emergencyDate = restrictionDatePhase1; - } else if(params.vacCamp.equals("ph1_180")){ // + - emergencyDate = restrictionDatePhase1; - } else if (params.vacCamp.equals("ph1_180_ph2_inf90vax180")) { - emergencyDate = restrictionDatePhase1; - dateToTurnDownMinDaysAfterInfection = restrictionDatePhase2; - // same as ifsg180 but after phase 2 date, minDaysAfterInfection = 90; - }else if(params.vacCamp.equals("ph2_90")){ - minDaysAfterInfection = 90; - minDaysAfterVaccination = 90; - emergencyDate = restrictionDatePhase2; - } else if(params.vacCamp.equals("ph2_inf90vax180")){ - minDaysAfterInfection = 90; - emergencyDate = restrictionDatePhase2; - } else if(params.vacCamp.equals("ph2_180")) { - emergencyDate = restrictionDatePhase2; - }else { - throw new RuntimeException(); - } - } +// if (params != null) { +// if (params.vacCamp.equals("base")) { // + +// +// } else if(params.vacCamp.equals("ph1_90")){ +// minDaysAfterInfection = 90; +// minDaysAfterVaccination = 90; +// +// emergencyDate = restrictionDatePhase1; +// } else if(params.vacCamp.equals("ph1_90vax180")){ +// minDaysAfterInfection = 90; +// emergencyDate = restrictionDatePhase1; +// } else if(params.vacCamp.equals("ph1_180")){ // + +// emergencyDate = restrictionDatePhase1; +// } else if (params.vacCamp.equals("ph1_180_ph2_inf90vax180")) { +// emergencyDate = restrictionDatePhase1; +// dateToTurnDownMinDaysAfterInfection = restrictionDatePhase2; +// // same as ifsg180 but after phase 2 date, minDaysAfterInfection = 90; +// }else if(params.vacCamp.equals("ph2_90")){ +// minDaysAfterInfection = 90; +// minDaysAfterVaccination = 90; +// emergencyDate = restrictionDatePhase2; +// } else if(params.vacCamp.equals("ph2_inf90vax180")){ +// minDaysAfterInfection = 90; +// emergencyDate = restrictionDatePhase2; +// } else if(params.vacCamp.equals("ph2_180")) { +// emergencyDate = restrictionDatePhase2; +// }else { +// throw new RuntimeException(); +// } +// } bind(VaccinationStrategyReoccurringCampaigns.Config.class).toInstance(new VaccinationStrategyReoccurringCampaigns.Config(startDateToVaccination, campaignDuration, vaccinationPool, minDaysAfterInfection, minDaysAfterVaccination, emergencyDate, dateToTurnDownMinDaysAfterInfection)); @@ -103,22 +103,28 @@ protected void configure() { double mutEscStrainA = 0.; double mutEscStrainB = 0.; + double mutEscBq = 0.; + double mutEscXbb = 0.; + if (params != null) { - mutEscStrainA = Double.parseDouble(params.StrainA); + mutEscBa5 = params.escBa5; +// mutEscStrainA = Double.parseDouble(params.escStrA); + mutEscBq = params.escBq; + mutEscXbb = params.escXbb; } //initial antibodies Map> initialAntibodies = new HashMap<>(); Map> antibodyRefreshFactors = new HashMap<>(); - configureAntibodies(initialAntibodies, antibodyRefreshFactors, mutEscDelta, mutEscBa1, mutEscBa5, mutEscStrainA, mutEscStrainB); + configureAntibodies(initialAntibodies, antibodyRefreshFactors, mutEscDelta, mutEscBa1, mutEscBa5, mutEscBq, mutEscXbb, mutEscStrainA, mutEscStrainB); AntibodyModel.Config antibodyConfig = new AntibodyModel.Config(initialAntibodies, antibodyRefreshFactors); - double immuneSigma = 3.0; if (params != null) { + double immuneSigma = params.immuneSigma;//3.0; antibodyConfig.setImmuneReponseSigma(immuneSigma); } @@ -130,7 +136,7 @@ protected void configure() { // designates a 35% of households as super safe; the susceptibility of that subpopulation is reduced to 1% wrt to general population. bind(HouseholdSusceptibility.Config.class).toInstance( HouseholdSusceptibility.newConfig() - .withSusceptibleHouseholds(0.35, 0.01) + .withSusceptibleHouseholds(params.pHouseholds, 0.01) //0.35 // .withNonVaccinableHouseholds(params.nonVaccinableHh) // .withShape(SnzCologneProductionScenario.INPUT.resolve("CologneDistricts.zip")) // .withFeature("STT_NAME", vingst, altstadtNord, bickendorf, weiden) @@ -140,7 +146,9 @@ protected void configure() { private void configureAntibodies(Map> initialAntibodies, Map> antibodyRefreshFactors, - double mutEscDelta, double mutEscBa1, double mutEscBa5, double mutEscStrainA, double mutEscStrainB) { + double mutEscDelta, double mutEscBa1, double mutEscBa5, + double mutEscBq, double mutEscXbb, + double mutEscStrainA, double mutEscStrainB) { for (VaccinationType immunityType : VaccinationType.values()) { initialAntibodies.put(immunityType, new EnumMap<>( VirusStrain.class ) ); for (VirusStrain virusStrain : VirusStrain.values()) { @@ -177,8 +185,10 @@ else if (immunityType == VaccinationType.vector) { initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.SARS_CoV_2, 0.01); initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.SARS_CoV_2, 0.01); initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.SARS_CoV_2, 0.01); - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.SARS_CoV_2, 0.01); - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.SARS_CoV_2, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.SARS_CoV_2, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.SARS_CoV_2, 0.01); initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.SARS_CoV_2, mRNAAlpha / mutEscDelta / mutEscBa1); initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.SARS_CoV_2, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscBa5); @@ -191,8 +201,10 @@ else if (immunityType == VaccinationType.vector) { initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.ALPHA, 0.01); initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.ALPHA, 0.01); initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.ALPHA, 0.01); - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.ALPHA, 0.01); - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.ALPHA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.ALPHA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.ALPHA, 0.01); initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.ALPHA, mRNAAlpha / mutEscDelta / mutEscBa1); initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.ALPHA, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscBa5); @@ -206,8 +218,10 @@ else if (immunityType == VaccinationType.vector) { initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.DELTA, 0.01); initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.DELTA, 0.01); initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.DELTA, 0.01); - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.DELTA, 0.01); - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.DELTA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.DELTA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.DELTA, 0.01); initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.DELTA, mRNADelta / mutEscBa1); initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.DELTA, mRNADelta / mutEscBa1 / mutEscBa5); @@ -221,8 +235,10 @@ else if (immunityType == VaccinationType.vector) { initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA1, 64.0 / 300.); initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / 1.4); initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5); //todo: is 1.4 - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5 / mutEscStrainA); - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscStrainB); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscXbb); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscStrainB); initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA1, mRNAAlpha); initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA1, mRNAAlpha / mutEscBa5); @@ -236,8 +252,10 @@ else if (immunityType == VaccinationType.vector) { initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / 1.4); initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA2, 64.0 / 300.); initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5); - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5 / mutEscStrainA); - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscStrainB); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscXbb); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscStrainB); initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA2, mRNAAlpha); initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA2, mRNAAlpha / mutEscBa5); @@ -252,41 +270,81 @@ else if (immunityType == VaccinationType.vector) { initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBa5);// todo: do we need 1.4? initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscBa5); initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA5, 64.0 / 300.); - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscStrainA); //todo ??? - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBa5 / mutEscStrainB); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBq); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscXbb / mutEscBa5 / mutEscBq); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscStrainA); //todo ??? +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBa5 / mutEscStrainB); initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA5, mRNAAlpha / mutEscBa5); initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA5, mRNAAlpha); + + //BQ1.1 + double mRNABq = mRNABa5 / mutEscBq; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.BQ, mRNABq); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.BQ, mRNABq * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.BQ, mRNABq * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.BQ, mRNABq * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.BQ, mRNABq * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.BQ, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.BQ, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.BQ, 64.0 / 300. / mutEscBq); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.BQ, 64.0 / 300.); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.BQ, 64.0 / 300. / mutEscXbb / mutEscBa5 / mutEscBq); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.BQ, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.BQ, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.BQ, mRNAAlpha / mutEscBa5 / mutEscBq); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.BQ, mRNAAlpha / mutEscBa5); + + + //XBB1.5 + + double mRNAXbb = mRNABA2 / mutEscXbb; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.XBB_15, mRNAXbb); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.XBB_15, mRNAXbb * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.XBB_15, mRNAXbb * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.XBB_15, mRNAXbb * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.XBB_15, mRNAXbb * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXbb); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXbb); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscBa5 / mutEscXbb); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscBq / mutEscBa5 / mutEscXbb); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.XBB_15, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.XBB_15, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.XBB_15, mRNAAlpha / mutEscXbb); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.XBB_15, mRNAAlpha / mutEscBa5 / mutEscXbb); + + //StrainA - double mRNAStrainA = mRNABa5 / mutEscStrainA; - initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.STRAIN_A, mRNAStrainA); - initialAntibodies.get(VaccinationType.vector).put(VirusStrain.STRAIN_A, mRNAStrainA * 4./20.); - initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); - initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); - initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.STRAIN_A, mRNAStrainA * 8./20.); - initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscBa5 /mutEscStrainA); - initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.STRAIN_A, 64.0 / 300./ mutEscBa5 /mutEscStrainA); - initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA); - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.STRAIN_A, 64.0 / 300.); - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); - initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscBa5 / mutEscStrainA); - initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscStrainA); +// double mRNAStrainA = mRNABa5 / mutEscStrainA; +// initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.STRAIN_A, mRNAStrainA); +// initialAntibodies.get(VaccinationType.vector).put(VirusStrain.STRAIN_A, mRNAStrainA * 4./20.); +// initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); +// initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); +// initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.STRAIN_A, mRNAStrainA * 8./20.); +// initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscBa5 /mutEscStrainA); +// initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.STRAIN_A, 64.0 / 300./ mutEscBa5 /mutEscStrainA); +// initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.STRAIN_A, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscStrainA); //StrainB - double mRNAStrainB = mRNABA2 / mutEscStrainB; - initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.STRAIN_B, mRNAStrainB); - initialAntibodies.get(VaccinationType.vector).put(VirusStrain.STRAIN_B, mRNAStrainB * 4./20.); - initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.STRAIN_B, mRNAStrainB * 6./20.); - initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.STRAIN_B, mRNAStrainB * 6./20.); - initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.STRAIN_B, mRNAStrainB * 8./20.); - initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.STRAIN_B, 64.0 / 300. / mutEscStrainB); - initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.STRAIN_B, 64.0 / 300./ mutEscStrainB); - initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.STRAIN_B, 64.0 / 300. / mutEscStrainB / mutEscBa5); - initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.STRAIN_B, 64.0 / 300./ mutEscStrainA / mutEscStrainB / mutEscBa5); - initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.STRAIN_B, 64.0 / 300.); - initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.STRAIN_B, mRNAAlpha / mutEscStrainB); - initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.STRAIN_B, mRNAAlpha / mutEscStrainB / mutEscBa5); - +// double mRNAStrainB = mRNABA2 / mutEscStrainB; +// initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.STRAIN_B, mRNAStrainB); +// initialAntibodies.get(VaccinationType.vector).put(VirusStrain.STRAIN_B, mRNAStrainB * 4./20.); +// initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.STRAIN_B, mRNAStrainB * 6./20.); +// initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.STRAIN_B, mRNAStrainB * 6./20.); +// initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.STRAIN_B, mRNAStrainB * 8./20.); +// initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.STRAIN_B, 64.0 / 300. / mutEscStrainB); +// initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.STRAIN_B, 64.0 / 300./ mutEscStrainB); +// initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.STRAIN_B, 64.0 / 300. / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.STRAIN_B, 64.0 / 300./ mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.STRAIN_B, 64.0 / 300.); +// initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.STRAIN_B, mRNAAlpha / mutEscStrainB); +// initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.STRAIN_B, mRNAAlpha / mutEscStrainB / mutEscBa5); +// for (VaccinationType immunityType : VaccinationType.values()) { antibodyRefreshFactors.put(immunityType, new EnumMap<>( VirusStrain.class ) ); @@ -375,7 +433,7 @@ public Config prepareConfig(int id, Params params) { EpisimConfigGroup episimConfig = ConfigUtils.addOrGetModule(config, EpisimConfigGroup.class); - episimConfig.setCalibrationParameter(episimConfig.getCalibrationParameter() * 1.2 * 1.7); + episimConfig.setCalibrationParameter(episimConfig.getCalibrationParameter() * 1.2 * 1.7 * params.theta); // start from immunization history @@ -395,14 +453,29 @@ public Config prepareConfig(int id, Params params) { VirusStrainConfigGroup virusStrainConfigGroup = ConfigUtils.addOrGetModule(config, VirusStrainConfigGroup.class); + // BQ.1.1 double ba5Inf = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA5).getInfectiousness(); double ba5Hos = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA5).getFactorSeriouslySick(); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setInfectiousness(ba5Inf); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorSeriouslySick(ba5Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorSeriouslySickVaccinated(ba5Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorCritical(ba5Hos); + - virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setInfectiousness(ba5Inf); - virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySick(ba5Hos); - virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySickVaccinated(ba5Hos); - virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorCritical(ba5Hos); + // XBB.1.5 + double ba2Inf = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA2).getInfectiousness(); + double ba2Hos = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA2).getFactorSeriouslySick(); + + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setInfectiousness(ba2Inf); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorSeriouslySick(ba2Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorSeriouslySickVaccinated(ba2Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorCritical(ba2Hos); + +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setInfectiousness(ba5Inf); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySick(ba5Hos); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySickVaccinated(ba5Hos); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorCritical(ba5Hos); //--------------------------------------- // I M P O R T @@ -434,6 +507,7 @@ public Config prepareConfig(int id, Params params) { + //--------------------------------------- // R E S T R I C T I O N S //--------------------------------------- @@ -443,94 +517,94 @@ public Config prepareConfig(int id, Params params) { //ifsg - if ("base".equals(params.ifsg)) { - - } else if ("45".equals(params.ifsg) || "90".equals(params.ifsg)) { - double compliance = Double.parseDouble(params.ifsg) / 100.; - builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of( - FaceMask.CLOTH, 0.0, - FaceMask.SURGICAL, compliance)), - "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); - builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, compliance)), "leisPublic"); - builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, compliance)), "shop_daily", "shop_other", "errands"); - builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, 0.9)), "pt"); // pt has 90 compliance either way - } else { - throw new RuntimeException(); - } +// if ("base".equals(params.ifsg)) { +// +// } else if ("45".equals(params.ifsg) || "90".equals(params.ifsg)) { +// double compliance = Double.parseDouble(params.ifsg) / 100.; +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, compliance)), +// "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, compliance)), "leisPublic"); +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, compliance)), "shop_daily", "shop_other", "errands"); +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, 0.9)), "pt"); // pt has 90 compliance either way +// } else { +// throw new RuntimeException(); +// } // EMERGENCY RESTRICTIONS //work - builder.restrict(LocalDate.parse("2022-10-15"), 0.88, "work", "business"); - double homeOfficeFactor = 0.5; - switch (params.work) { - case "base": - break; - case "half": - builder.restrict(restrictionDatePhase2, 0.88 * homeOfficeFactor, "work"); // dont include business bc harder to do from home office - builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> rf * homeOfficeFactor, "work"); - break; - case "half&mask": - builder.restrict(restrictionDatePhase2, 0.88 * homeOfficeFactor, "work"); // dont include business bc harder to do from home office - builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> rf * homeOfficeFactor, "work"); - - builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of(FaceMask.SURGICAL, 0.0, FaceMask.N95, 0.9)), "work", "business"); - break; - default: - throw new RuntimeException("invalid parameter"); - } +// builder.restrict(LocalDate.parse("2022-10-15"), 0.88, "work", "business"); +// double homeOfficeFactor = 0.5; +// switch (params.work) { +// case "base": +// break; +// case "half": +// builder.restrict(restrictionDatePhase2, 0.88 * homeOfficeFactor, "work"); // dont include business bc harder to do from home office +// builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> rf * homeOfficeFactor, "work"); +// break; +// case "half&mask": +// builder.restrict(restrictionDatePhase2, 0.88 * homeOfficeFactor, "work"); // dont include business bc harder to do from home office +// builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> rf * homeOfficeFactor, "work"); +// +// builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of(FaceMask.SURGICAL, 0.0, FaceMask.N95, 0.9)), "work", "business"); +// break; +// default: +// throw new RuntimeException("invalid parameter"); +// } // leisure public + private - switch (params.leis) { - case "base": - break; - case "pub50": - builder.restrict(restrictionDatePhase2, 0.88 * 0.5, "leisPublic"); - break; - case "pubPriv50": - builder.restrict(restrictionDatePhase2, 0.88 * 0.5, "leisPublic", "leisPrivate"); - break; - default: - throw new RuntimeException("invalid parameter"); - } +// switch (params.leis) { +// case "base": +// break; +// case "pub50": +// builder.restrict(restrictionDatePhase2, 0.88 * 0.5, "leisPublic"); +// break; +// case "pubPriv50": +// builder.restrict(restrictionDatePhase2, 0.88 * 0.5, "leisPublic", "leisPrivate"); +// break; +// default: +// throw new RuntimeException("invalid parameter"); +// } //school - switch (params.edu) { - case "base": - break; - case "mask": - builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of( - FaceMask.CLOTH, 0.0, - FaceMask.SURGICAL, 0.0, - FaceMask.N95, 0.90)), - "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); - break; - case "half&mask": - builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of( - FaceMask.CLOTH, 0.0, - FaceMask.SURGICAL, 0.0, - FaceMask.N95, 0.90)), - "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); - - builder.restrict(restrictionDatePhase2, 0.5, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); - builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> Math.min(0.5, rf), "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); - break; - default: - throw new RuntimeException("invalid parameter"); - } - - if ("base".equals(params.maskPt)) { - } else { - builder.restrict(LocalDate.parse(params.maskPt), - Restriction.ofMask(Map.of( - FaceMask.CLOTH, 0.0, - FaceMask.SURGICAL, 0.0, - FaceMask.N95, 0.0)), - "pt"); - } +// switch (params.edu) { +// case "base": +// break; +// case "mask": +// builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, 0.0, +// FaceMask.N95, 0.90)), +// "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// break; +// case "half&mask": +// builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, 0.0, +// FaceMask.N95, 0.90)), +// "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// +// builder.restrict(restrictionDatePhase2, 0.5, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> Math.min(0.5, rf), "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// break; +// default: +// throw new RuntimeException("invalid parameter"); +// } + +// if ("base".equals(params.maskPt)) { +// } else { +// builder.restrict(LocalDate.parse(params.maskPt), +// Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, 0.0, +// FaceMask.N95, 0.0)), +// "pt"); +// } // vary amount of "school" activity that takes place during vacation - builder.restrict(LocalDate.parse("2022-06-27"), 0.8, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); +// builder.restrict(LocalDate.parse("2022-06-27"), params.eduRfVacation, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); episimConfig.setPolicy(builder.build()); @@ -550,13 +624,45 @@ public Config prepareConfig(int id, Params params) { private void configureFutureDiseaseImport(Params params, EpisimConfigGroup episimConfig) { Map infPerDayBa1 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA1, new TreeMap<>())); Map infPerDayBa2 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA2, new TreeMap<>())); - Map infPerDayBa5 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA5, new TreeMap<>())); +// Map infPerDayBa5 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA5, new TreeMap<>())); + Map infPerDayBq = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.BQ, new TreeMap<>())); + Map infPerDayXbb = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.XBB_15, new TreeMap<>())); Map infPerDayStrA = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.STRAIN_A, new TreeMap<>())); + // Override BA5 + + Map infPerDayBa5 = new HashMap<>(); + infPerDayBa5.put(LocalDate.parse("2020-01-01"), 0); + LocalDate ba5Date = LocalDate.parse(params.ba5Date); + + for (int i = 0; i < 7; i++) { + infPerDayBa5.put(ba5Date.plusDays(i), 4); + } + infPerDayBa5.put(ba5Date.plusDays(7), 1); + + + LocalDate startDate = LocalDate.parse(params.bqXbbStartDate); // 2022-09-19 + + + // BQ.1.1 + infPerDayBq.put(LocalDate.parse("2020-01-01"), 0); + for (int i = 0; i < 7; i++) { + infPerDayBq.put(startDate.plusDays(i), 4); + } + infPerDayBq.put(startDate.plusDays(7), 1); + + // XBB.1.5 + infPerDayXbb.put(LocalDate.parse("2020-01-01"), 0); + for (int i = 0; i < 7; i++) { + infPerDayXbb.put(startDate.plusDays(i), 4); + } + infPerDayXbb.put(startDate.plusDays(7), 1); + + //StrainA - if (!params.StrainA.equals("off")) { + if (!params.escStrA.equals("off")) { infPerDayStrA.put(LocalDate.parse("2020-01-01"), 0); - LocalDate strADate = LocalDate.parse(params.strainADate); + LocalDate strADate = LocalDate.parse(params.strainADate); // 2022-09-19 for (int i = 0; i < 7; i++) { infPerDayStrA.put(strADate.plusDays(i), 4); @@ -568,66 +674,98 @@ private void configureFutureDiseaseImport(Params params, EpisimConfigGroup episi episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA1, infPerDayBa1); episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA2, infPerDayBa2); episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA5, infPerDayBa5); + episimConfig.setInfections_pers_per_day(VirusStrain.BQ, infPerDayBq); + episimConfig.setInfections_pers_per_day(VirusStrain.XBB_15, infPerDayXbb); - if (!params.StrainA.equals("off")) { + if (!params.escStrA.equals("off")) { episimConfig.setInfections_pers_per_day(VirusStrain.STRAIN_A, infPerDayStrA); } } public static final class Params { // general - @GenerateSeeds(20) + @GenerateSeeds(5) // 5 public long seed; - @StringParameter({"sepSeeds", "4711"}) - public String startFromImm; + @Parameter({1.8, 3.0, 4.0}) // 3 + public double escBq; - // BQ 1 -// @StringParameter({"off", "2.0", "2.25", "2.5", "2.75", "3.0"}) -// @StringParameter({"1.0", "1.2", "1.4", "1.6", "1.8", "2.0"}) - @StringParameter({"2.0"}) - public String StrainA; + @Parameter({1.8 * 5, 3.5 * 5, 4.5 * 5}) // 3 + public double escXbb; + + @StringParameter({"2022-09-11", "2022-09-15", "2022-09-19"}) // 3 + public String bqXbbStartDate; -// @StringParameter({"base", "2022-11-15", "2022-12-01", "2022-12-15", "2023-01-01"}) - @StringParameter({"base"}) - public String maskPt; +// @Parameter({0.0,0.2,0.4,0.6,0.8,1.0}) +// public double eduRfVacation; +// @Parameter({0.0, 0.3, 0.6}) + @Parameter({0.3}) + public double immuneSigma; +// @Parameter({4.2, 4.6, 5.0}) + @Parameter({4.6}) + public double escBa5; + +// @StringParameter({"2022-04-10", "2022-04-13", "2022-04-16", "2022-04-19"}) + @StringParameter({"2022-04-13"}) + public String ba5Date; + +// @Parameter({0.35, 0.5, 0.65}) + @Parameter({0.5}) + public double pHouseholds; + +// @Parameter({1.0, 1.1, 1.2}) + @Parameter({1.2}) + public double theta; + + + @StringParameter({"sepSeeds"}) + public String startFromImm; + +// @Parameter({1.0}) +// public double infBa5; + + // BQ 1 +// @StringParameter({"off", "2.0", "2.25", "2.5", "2.75", "3.0"}) + @StringParameter({"off"}) //1.8 + public String escStrA; -// @StringParameter({"2022-08-24", "2022-08-29", "2022-09-04", "2022-09-09", "2022-09-14", "2022-09-19"}) @StringParameter({"2022-09-19"}) public String strainADate; + // @StringParameter({"2022-08-24", "2022-08-29", "2022-09-04", "2022-09-09", "2022-09-14", "2022-09-19"}) + // @StringParameter({"base", "2022-11-15", "2022-12-01", "2022-12-15", "2023-01-01"}) +// @StringParameter({"base"}) +// public String maskPt; //IFSG - @StringParameter({"base"}) - public String ifsg; +// @StringParameter({"base"}) +// public String ifsg; // Vaccination Campaign - @StringParameter({"base"}) - String vacCamp; +// @StringParameter({"base"}) +// String vacCamp; // NEW RESTRICTIONS - @StringParameter({"base"}) - public String work; +// @StringParameter({"base"}) +// public String work; // leisure Public - @StringParameter({"base"}) - public String leis; +// @StringParameter({"base"}) +// public String leis; //edu - @StringParameter({"base"}) - public String edu; - - +// @StringParameter({"base"}) +// public String edu; } public static void main(String[] args) { String[] args2 = { - RunParallel.OPTION_SETUP, CologneBMBF202212XX_bq1.class.getName(), + RunParallel.OPTION_SETUP, CologneBMBF2023.class.getName(), RunParallel.OPTION_PARAMS, Params.class.getName(), RunParallel.OPTION_TASKS, Integer.toString(1), RunParallel.OPTION_ITERATIONS, Integer.toString(1000), diff --git a/src/main/java/org/matsim/run/batch/CologneBMBF202302XX.java b/src/main/java/org/matsim/run/batch/CologneBMBF202302XX.java new file mode 100644 index 000000000..ffd592acd --- /dev/null +++ b/src/main/java/org/matsim/run/batch/CologneBMBF202302XX.java @@ -0,0 +1,809 @@ +package org.matsim.run.batch; + +import com.google.inject.AbstractModule; +import com.google.inject.Module; +import com.google.inject.Singleton; +import com.google.inject.multibindings.Multibinder; +import com.google.inject.util.Modules; +import org.matsim.core.config.Config; +import org.matsim.core.config.ConfigUtils; +import org.matsim.episim.BatchRun; +import org.matsim.episim.EpisimConfigGroup; +import org.matsim.episim.VirusStrainConfigGroup; +import org.matsim.episim.analysis.*; +import org.matsim.episim.model.*; +import org.matsim.episim.model.listener.HouseholdSusceptibility; +import org.matsim.episim.model.vaccination.VaccinationModel; +import org.matsim.episim.model.vaccination.VaccinationStrategyReoccurringCampaigns; +import org.matsim.episim.policy.FixedPolicy; +import org.matsim.run.RunParallel; +import org.matsim.run.modules.SnzCologneProductionScenario; +import org.matsim.run.modules.SnzProductionScenario; + +import javax.annotation.Nullable; +import java.time.LocalDate; +import java.util.*; + + +/** + * Batch for Bmbf runs + */ +public class CologneBMBF202302XX implements BatchRun { + + boolean DEBUG_MODE = false; + + String START_DATE = "2022-04-01"; + int runCount = 0; + + LocalDate restrictionDatePhase1 = LocalDate.parse("2022-12-01"); + LocalDate restrictionDatePhase2 = restrictionDatePhase1.plusDays(10); + + @Nullable + @Override + public Module getBindings(int id, @Nullable Params params) { + return Modules.override(getBindings(0.0, params)).with(new AbstractModule() { + @Override + protected void configure() { + + // VACCINATION MODEL + Multibinder set = Multibinder.newSetBinder(binder(), VaccinationModel.class); + set.addBinding().to(VaccinationStrategyReoccurringCampaigns.class).in(Singleton.class); + // fixed values + LocalDate start = LocalDate.parse("2022-10-15"); + VaccinationType vaccinationType = VaccinationType.ba5Update; + int campaignDuration = 300000; + + // default values, to be changed if params != null + int minDaysAfterInfection = 180; + int minDaysAfterVaccination = 180; + VaccinationStrategyReoccurringCampaigns.Config.VaccinationPool vaccinationPool = VaccinationStrategyReoccurringCampaigns.Config.VaccinationPool.vaccinated; + LocalDate emergencyDate = LocalDate.MAX; + LocalDate dateToTurnDownMinDaysAfterInfection = LocalDate.MAX; + Map startDateToVaccination = new HashMap<>(); + startDateToVaccination.put(start, vaccinationType); + +// if (params != null) { +// if (params.vacCamp.equals("base")) { // + +// +// } else if(params.vacCamp.equals("ph1_90")){ +// minDaysAfterInfection = 90; +// minDaysAfterVaccination = 90; +// +// emergencyDate = restrictionDatePhase1; +// } else if(params.vacCamp.equals("ph1_90vax180")){ +// minDaysAfterInfection = 90; +// emergencyDate = restrictionDatePhase1; +// } else if(params.vacCamp.equals("ph1_180")){ // + +// emergencyDate = restrictionDatePhase1; +// } else if (params.vacCamp.equals("ph1_180_ph2_inf90vax180")) { +// emergencyDate = restrictionDatePhase1; +// dateToTurnDownMinDaysAfterInfection = restrictionDatePhase2; +// // same as ifsg180 but after phase 2 date, minDaysAfterInfection = 90; +// }else if(params.vacCamp.equals("ph2_90")){ +// minDaysAfterInfection = 90; +// minDaysAfterVaccination = 90; +// emergencyDate = restrictionDatePhase2; +// } else if(params.vacCamp.equals("ph2_inf90vax180")){ +// minDaysAfterInfection = 90; +// emergencyDate = restrictionDatePhase2; +// } else if(params.vacCamp.equals("ph2_180")) { +// emergencyDate = restrictionDatePhase2; +// }else { +// throw new RuntimeException(); +// } +// } + + bind(VaccinationStrategyReoccurringCampaigns.Config.class).toInstance(new VaccinationStrategyReoccurringCampaigns.Config(startDateToVaccination, campaignDuration, vaccinationPool, minDaysAfterInfection, minDaysAfterVaccination, emergencyDate, dateToTurnDownMinDaysAfterInfection)); + + + // ANTIBODY MODEL + // default values + double mutEscDelta = 29.2 / 10.9; + double mutEscBa1 = 10.9 / 1.9; + double mutEscBa5 = 5.0; + + double mutEscStrainA = 0.; + double mutEscStrainB = 0.; + + double mutEscBq = 0.; + double mutEscXbb = 0.; + + + if (params != null) { + + mutEscBa5 = params.escBa5; +// mutEscStrainA = Double.parseDouble(params.escStrA); + mutEscBq = params.escBq; + mutEscXbb = params.escXbb; + + } + + //initial antibodies + Map> initialAntibodies = new HashMap<>(); + Map> antibodyRefreshFactors = new HashMap<>(); + configureAntibodies(initialAntibodies, antibodyRefreshFactors, mutEscDelta, mutEscBa1, mutEscBa5, mutEscBq, mutEscXbb, mutEscStrainA, mutEscStrainB); + + AntibodyModel.Config antibodyConfig = new AntibodyModel.Config(initialAntibodies, antibodyRefreshFactors); + + if (params != null) { + double immuneSigma = params.immuneSigma;//3.0; + antibodyConfig.setImmuneReponseSigma(immuneSigma); + } + + bind(AntibodyModel.Config.class).toInstance(antibodyConfig); + + if (params == null) return; + + // HOUSEHOLD SUSCEPTIBILITY + // designates a 35% of households as super safe; the susceptibility of that subpopulation is reduced to 1% wrt to general population. + bind(HouseholdSusceptibility.Config.class).toInstance( + HouseholdSusceptibility.newConfig() + .withSusceptibleHouseholds(params.pHouseholds, 0.01) //0.35 +// .withNonVaccinableHouseholds(params.nonVaccinableHh) +// .withShape(SnzCologneProductionScenario.INPUT.resolve("CologneDistricts.zip")) +// .withFeature("STT_NAME", vingst, altstadtNord, bickendorf, weiden) + ); + + } + + private void configureAntibodies(Map> initialAntibodies, + Map> antibodyRefreshFactors, + double mutEscDelta, double mutEscBa1, double mutEscBa5, + double mutEscBq, double mutEscXbb, + double mutEscStrainA, double mutEscStrainB) { + for (VaccinationType immunityType : VaccinationType.values()) { + initialAntibodies.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + + if (immunityType == VaccinationType.mRNA) { + initialAntibodies.get(immunityType).put(virusStrain, 29.2); //10.0 + } + else if (immunityType == VaccinationType.vector) { + initialAntibodies.get(immunityType).put(virusStrain, 6.8); //2.5 + } + else { + initialAntibodies.get(immunityType).put(virusStrain, 5.0); + } + } + } + + for (VirusStrain immunityType : VirusStrain.values()) { + initialAntibodies.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + initialAntibodies.get(immunityType).put(virusStrain, 5.0); + } + } + + //mRNAAlpha, mRNADelta, mRNABA1 comes from Sydney's calibration. + //The other values come from Rössler et al. + + //Wildtype + double mRNAAlpha = 29.2; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.SARS_CoV_2, mRNAAlpha); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.SARS_CoV_2, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.SARS_CoV_2, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.SARS_CoV_2, mRNAAlpha / mutEscDelta / mutEscBa1); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.SARS_CoV_2, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscBa5); + + //Alpha + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.ALPHA, mRNAAlpha); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.ALPHA, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.ALPHA, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.ALPHA, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.ALPHA, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.ALPHA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.ALPHA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.ALPHA, mRNAAlpha / mutEscDelta / mutEscBa1); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.ALPHA, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscBa5); + + //DELTA + double mRNADelta = mRNAAlpha / mutEscDelta; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.DELTA, mRNADelta); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.DELTA, mRNADelta * 150./300.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.DELTA, mRNADelta * 64./300.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.DELTA, mRNADelta * 64./300.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.DELTA, mRNADelta * 450./300.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.DELTA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.DELTA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.DELTA, mRNADelta / mutEscBa1); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.DELTA, mRNADelta / mutEscBa1 / mutEscBa5); + + //BA.1 + double mRNABA1 = mRNADelta / mutEscBa1; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.OMICRON_BA1, mRNABA1); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.OMICRON_BA1, mRNABA1 * 4./20.); //??? + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.OMICRON_BA1, mRNABA1 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.OMICRON_BA1, mRNABA1 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.OMICRON_BA1, mRNABA1 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA1, 64.0 / 300.); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / 1.4); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5); //todo: is 1.4 + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscXbb); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscStrainB); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA1, mRNAAlpha); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA1, mRNAAlpha / mutEscBa5); + + //BA.2 + double mRNABA2 = mRNABA1; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.OMICRON_BA2, mRNABA2); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.OMICRON_BA2, mRNABA2 * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.OMICRON_BA2, mRNABA2 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.OMICRON_BA2, mRNABA2 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.OMICRON_BA2, mRNABA2 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / 1.4); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA2, 64.0 / 300.); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscXbb); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscStrainB); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA2, mRNAAlpha); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA2, mRNAAlpha / mutEscBa5); + + + //BA.5 + double mRNABa5 = mRNABA2 / mutEscBa5; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.OMICRON_BA5, mRNABa5); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.OMICRON_BA5, mRNABa5 * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.OMICRON_BA5, mRNABa5 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.OMICRON_BA5, mRNABa5 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.OMICRON_BA5, mRNABa5 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBa5);// todo: do we need 1.4? + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscBa5); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA5, 64.0 / 300.); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBq); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscXbb / mutEscBa5 / mutEscBq); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscStrainA); //todo ??? +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBa5 / mutEscStrainB); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA5, mRNAAlpha / mutEscBa5); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA5, mRNAAlpha); + + + //BQ1.1 + double mRNABq = mRNABa5 / mutEscBq; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.BQ, mRNABq); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.BQ, mRNABq * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.BQ, mRNABq * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.BQ, mRNABq * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.BQ, mRNABq * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.BQ, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.BQ, 64.0 / 300. / mutEscBa5 / mutEscBq); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.BQ, 64.0 / 300. / mutEscBq); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.BQ, 64.0 / 300.); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.BQ, 64.0 / 300. / mutEscXbb / mutEscBa5 / mutEscBq); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.BQ, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.BQ, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.BQ, mRNAAlpha / mutEscBa5 / mutEscBq); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.BQ, mRNAAlpha / mutEscBa5); + + + //XBB1.5 + + double mRNAXbb = mRNABA2 / mutEscXbb; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.XBB_15, mRNAXbb); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.XBB_15, mRNAXbb * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.XBB_15, mRNAXbb * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.XBB_15, mRNAXbb * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.XBB_15, mRNAXbb * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXbb); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXbb); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscBa5 / mutEscXbb); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscBq / mutEscBa5 / mutEscXbb); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.XBB_15, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.XBB_15, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.XBB_15, mRNAAlpha / mutEscXbb); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.XBB_15, mRNAAlpha / mutEscBa5 / mutEscXbb); + + + //StrainA +// double mRNAStrainA = mRNABa5 / mutEscStrainA; +// initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.STRAIN_A, mRNAStrainA); +// initialAntibodies.get(VaccinationType.vector).put(VirusStrain.STRAIN_A, mRNAStrainA * 4./20.); +// initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); +// initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); +// initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.STRAIN_A, mRNAStrainA * 8./20.); +// initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscBa5 /mutEscStrainA); +// initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.STRAIN_A, 64.0 / 300./ mutEscBa5 /mutEscStrainA); +// initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.STRAIN_A, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscStrainA); + + //StrainB +// double mRNAStrainB = mRNABA2 / mutEscStrainB; +// initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.STRAIN_B, mRNAStrainB); +// initialAntibodies.get(VaccinationType.vector).put(VirusStrain.STRAIN_B, mRNAStrainB * 4./20.); +// initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.STRAIN_B, mRNAStrainB * 6./20.); +// initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.STRAIN_B, mRNAStrainB * 6./20.); +// initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.STRAIN_B, mRNAStrainB * 8./20.); +// initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.STRAIN_B, 64.0 / 300. / mutEscStrainB); +// initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.STRAIN_B, 64.0 / 300./ mutEscStrainB); +// initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.STRAIN_B, 64.0 / 300. / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.STRAIN_B, 64.0 / 300./ mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.STRAIN_B, 64.0 / 300.); +// initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.STRAIN_B, mRNAAlpha / mutEscStrainB); +// initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.STRAIN_B, mRNAAlpha / mutEscStrainB / mutEscBa5); +// + + for (VaccinationType immunityType : VaccinationType.values()) { + antibodyRefreshFactors.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + + if (immunityType == VaccinationType.mRNA) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } + else if (immunityType == VaccinationType.vector) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 5.0); + } + else if (immunityType == VaccinationType.ba1Update) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } + else if (immunityType == VaccinationType.ba5Update) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } + else { + antibodyRefreshFactors.get(immunityType).put(virusStrain, Double.NaN); + } + + } + } + + for (VirusStrain immunityType : VirusStrain.values()) { + antibodyRefreshFactors.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } + } + + +// UtilsJR.printInitialAntibodiesToConsole(initialAntibodies); + + } + }); + + } + + private SnzCologneProductionScenario getBindings(double pHousehold, Params params) { + return new SnzCologneProductionScenario.Builder() + .setCarnivalModel(SnzCologneProductionScenario.CarnivalModel.yes) + .setSebastianUpdate(true) + .setLeisureCorrection(1.3) //params == null ? 0.0 : params.actCorrection) + .setScaleForActivityLevels(1.3) + .setSuscHouseholds_pct(pHousehold) + .setWeatherModel(params != null && params.outdoorLeisure.equals("false") ? SnzProductionScenario.WeatherModel.midpoints_500_500 : SnzProductionScenario.WeatherModel.midpoints_185_250 ) + .setActivityHandling(EpisimConfigGroup.ActivityHandling.startOfDay) +// .setTestingModel(params != null ? FlexibleTestingModel.class : DefaultTestingModel.class) + .setInfectionModel(InfectionModelWithAntibodies.class) + .build(); + } + + @Override + public Metadata getMetadata() { + return Metadata.of("cologne", "calibration"); + } + + @Override + public Collection postProcessing() { + return List.of( + new VaccinationEffectiveness().withArgs("--start-date", START_DATE), + new RValuesFromEvents().withArgs("--start-date", START_DATE), +// new VaccinationEffectivenessFromPotentialInfections().withArgs("--remove-infected"), + new FilterEvents().withArgs("--output", "./output/"), + new HospitalNumbersFromEvents().withArgs("--output", "./output/", "--input", "/scratch/projects/bzz0020/episim-input", "--start-date", START_DATE) +// new SecondaryAttackRateFromEvents().withArgs() + ); + } + + @Override + public Config prepareConfig(int id, Params params) { + + if (DEBUG_MODE) { + if (runCount == 0) { + runCount++; + } else { + return null; + } + } + + SnzCologneProductionScenario module = getBindings(0.0, params); + + Config config = module.config(); + + config.global().setRandomSeed(params.seed); + + EpisimConfigGroup episimConfig = ConfigUtils.addOrGetModule(config, EpisimConfigGroup.class); + + episimConfig.setCalibrationParameter(episimConfig.getCalibrationParameter() * 1.2 * 1.7 * params.theta); + + + // start from immunization history +// episimConfig.setStartDate(LocalDate.parse(START_DATE)); + +// if (params.startFromImm.equals("sepSeeds")) { +// episimConfig.setStartFromImmunization("/scratch/projects/bzz0020/episim-input/imm-cologne-2022-11-24/" + params.seed + "-events_reduced.tar"); +// } else if (params.startFromImm.equals("4711")) { +// episimConfig.setStartFromImmunization("/scratch/projects/bzz0020/episim-input/imm-cologne-2022-11-24/4711-events_reduced.tar"); +// } else { +// throw new RuntimeException("invalid param"); +// } + + //--------------------------------------- + // S T R A I N S + //--------------------------------------- + + VirusStrainConfigGroup virusStrainConfigGroup = ConfigUtils.addOrGetModule(config, VirusStrainConfigGroup.class); + + // BQ.1.1 + double ba5Inf = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA5).getInfectiousness(); + double ba5Hos = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA5).getFactorSeriouslySick(); + + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setInfectiousness(ba5Inf); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorSeriouslySick(ba5Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorSeriouslySickVaccinated(ba5Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorCritical(ba5Hos); + + + // XBB.1.5 + double ba2Inf = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA2).getInfectiousness(); + double ba2Hos = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA2).getFactorSeriouslySick(); + + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setInfectiousness(ba2Inf); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorSeriouslySick(ba2Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorSeriouslySickVaccinated(ba2Hos); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorCritical(ba2Hos); + +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setInfectiousness(ba5Inf); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySick(ba5Hos); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySickVaccinated(ba5Hos); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorCritical(ba5Hos); + + //--------------------------------------- + // I M P O R T + //--------------------------------------- + + configureFutureDiseaseImport(params, episimConfig); + + // modify import: + LocalDate impModDate = LocalDate.parse("2022-01-31"); + double impRedBa1 = 0.0; + double impRedBa2 = 0.0; + if (impRedBa1 != 1.0) { + NavigableMap impBa1 = episimConfig.getInfections_pers_per_day().get(VirusStrain.OMICRON_BA1); + for (Map.Entry entry : impBa1.entrySet()) { + if (entry.getKey().isAfter(impModDate)) { + impBa1.put(entry.getKey(), (int) (entry.getValue() * impRedBa1)); + } + } + } + + if (impRedBa2 != 1.0) { + NavigableMap impBa2 = episimConfig.getInfections_pers_per_day().get(VirusStrain.OMICRON_BA2); + for (Map.Entry entry : impBa2.entrySet()) { + if (entry.getKey().isAfter(impModDate)) { + impBa2.put(entry.getKey(), (int) (entry.getValue() * impRedBa2)); + } + } + } + + + + + //--------------------------------------- + // R E S T R I C T I O N S + //--------------------------------------- + + FixedPolicy.ConfigBuilder builder = FixedPolicy.parse(episimConfig.getPolicy()); + + + + //ifsg +// if ("base".equals(params.ifsg)) { +// +// } else if ("45".equals(params.ifsg) || "90".equals(params.ifsg)) { +// double compliance = Double.parseDouble(params.ifsg) / 100.; +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, compliance)), +// "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, compliance)), "leisPublic"); +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, compliance)), "shop_daily", "shop_other", "errands"); +// builder.restrict(restrictionDatePhase1, Restriction.ofMask(Map.of(FaceMask.N95, 0.9)), "pt"); // pt has 90 compliance either way +// } else { +// throw new RuntimeException(); +// } + + // EMERGENCY RESTRICTIONS + //work +// builder.restrict(LocalDate.parse("2022-10-15"), 0.88, "work", "business"); +// double homeOfficeFactor = 0.5; +// switch (params.work) { +// case "base": +// break; +// case "half": +// builder.restrict(restrictionDatePhase2, 0.88 * homeOfficeFactor, "work"); // dont include business bc harder to do from home office +// builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> rf * homeOfficeFactor, "work"); +// break; +// case "half&mask": +// builder.restrict(restrictionDatePhase2, 0.88 * homeOfficeFactor, "work"); // dont include business bc harder to do from home office +// builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> rf * homeOfficeFactor, "work"); +// +// builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of(FaceMask.SURGICAL, 0.0, FaceMask.N95, 0.9)), "work", "business"); +// break; +// default: +// throw new RuntimeException("invalid parameter"); +// } + + // leisure public + private +// switch (params.leis) { +// case "base": +// break; +// case "pub50": +// builder.restrict(restrictionDatePhase2, 0.88 * 0.5, "leisPublic"); +// break; +// case "pubPriv50": +// builder.restrict(restrictionDatePhase2, 0.88 * 0.5, "leisPublic", "leisPrivate"); +// break; +// default: +// throw new RuntimeException("invalid parameter"); +// } + + //school +// switch (params.edu) { +// case "base": +// break; +// case "mask": +// builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, 0.0, +// FaceMask.N95, 0.90)), +// "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// break; +// case "half&mask": +// builder.restrict(restrictionDatePhase2, Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, 0.0, +// FaceMask.N95, 0.90)), +// "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// +// builder.restrict(restrictionDatePhase2, 0.5, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// builder.applyToRf(restrictionDatePhase2.plusDays(1).toString(), restrictionDatePhase2.plusDays(1000).toString(), (d, rf) -> Math.min(0.5, rf), "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other", "educ_higher"); +// break; +// default: +// throw new RuntimeException("invalid parameter"); +// } + +// if ("base".equals(params.maskPt)) { +// } else { +// builder.restrict(LocalDate.parse(params.maskPt), +// Restriction.ofMask(Map.of( +// FaceMask.CLOTH, 0.0, +// FaceMask.SURGICAL, 0.0, +// FaceMask.N95, 0.0)), +// "pt"); +// } + + + // vary amount of "school" activity that takes place during vacation +// builder.restrict(LocalDate.parse("2022-06-27"), params.eduRfVacation, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); + + switch (params.edu2020) { + case "baseLateRes": + break; + case "earlyRes": + builder.restrict(LocalDate.parse("2020-12-16"), 0.2, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); + break; + case "noRes": + builder.restrict(LocalDate.parse("2021-01-06"), 1.0, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); + break; + case "noResNoXmas": + builder.restrict(LocalDate.parse("2020-12-23"), 1.0, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); + builder.restrict(LocalDate.parse("2021-01-06"), 1.0, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); + break; + default: + throw new RuntimeException(); + + } + + + episimConfig.setPolicy(builder.build()); + + + //--------------------------------------- + // M I S C + //--------------------------------------- + + + if (DEBUG_MODE) { + UtilsJR.produceDiseaseImportPlot(episimConfig.getInfections_pers_per_day()); + } + + return config; + } + + private void configureFutureDiseaseImport(Params params, EpisimConfigGroup episimConfig) { + Map infPerDayBa1 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA1, new TreeMap<>())); + Map infPerDayBa2 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA2, new TreeMap<>())); +// Map infPerDayBa5 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA5, new TreeMap<>())); + Map infPerDayBq = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.BQ, new TreeMap<>())); + Map infPerDayXbb = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.XBB_15, new TreeMap<>())); + Map infPerDayStrA = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.STRAIN_A, new TreeMap<>())); + + // Override BA5 + + Map infPerDayBa5 = new HashMap<>(); + infPerDayBa5.put(LocalDate.parse("2020-01-01"), 0); + LocalDate ba5Date = LocalDate.parse(params.ba5Date); + + for (int i = 0; i < 7; i++) { + infPerDayBa5.put(ba5Date.plusDays(i), 4); + } + infPerDayBa5.put(ba5Date.plusDays(7), 1); + + + LocalDate startDate = LocalDate.parse(params.bqXbbStartDate); // 2022-09-19 + + + // BQ.1.1 + infPerDayBq.put(LocalDate.parse("2020-01-01"), 0); + for (int i = 0; i < 7; i++) { + infPerDayBq.put(startDate.plusDays(i), 4); + } + infPerDayBq.put(startDate.plusDays(7), 1); + + // XBB.1.5 + infPerDayXbb.put(LocalDate.parse("2020-01-01"), 0); + for (int i = 0; i < 7; i++) { + infPerDayXbb.put(startDate.plusDays(i), 4); + } + infPerDayXbb.put(startDate.plusDays(7), 1); + + + //StrainA + if (!params.escStrA.equals("off")) { + infPerDayStrA.put(LocalDate.parse("2020-01-01"), 0); + LocalDate strADate = LocalDate.parse(params.strainADate); // 2022-09-19 + + for (int i = 0; i < 7; i++) { + infPerDayStrA.put(strADate.plusDays(i), 4); + } + infPerDayStrA.put(strADate.plusDays(7), 1); + } + + // save disease import + episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA1, infPerDayBa1); + episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA2, infPerDayBa2); + episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA5, infPerDayBa5); + episimConfig.setInfections_pers_per_day(VirusStrain.BQ, infPerDayBq); + episimConfig.setInfections_pers_per_day(VirusStrain.XBB_15, infPerDayXbb); + + if (!params.escStrA.equals("off")) { + episimConfig.setInfections_pers_per_day(VirusStrain.STRAIN_A, infPerDayStrA); + } + } + + public static final class Params { + // general + @GenerateSeeds(5) // 5 + public long seed; + +// @StringParameter({"true", "false"}) + @StringParameter({"false"}) + public String outdoorLeisure; + + @StringParameter({"baseLateRes", "noRes","noResNoXmas", "earlyRes"}) + public String edu2020; + + @Parameter({3.0}) // 3 + public double escBq; + + @Parameter({3.5 * 5}) // 3 + public double escXbb; + + @StringParameter({"2022-09-15"}) // 3 + public String bqXbbStartDate; + +// @Parameter({0.0,0.2,0.4,0.6,0.8,1.0}) +// public double eduRfVacation; + + // @Parameter({0.0, 0.3, 0.6}) + @Parameter({0.3}) + public double immuneSigma; + + // @Parameter({4.2, 4.6, 5.0}) + @Parameter({4.6}) + public double escBa5; + + // @StringParameter({"2022-04-10", "2022-04-13", "2022-04-16", "2022-04-19"}) + @StringParameter({"2022-04-13"}) + public String ba5Date; + + // @Parameter({0.35, 0.5, 0.65}) + @Parameter({0.5}) + public double pHouseholds; + + // @Parameter({1.0, 1.1, 1.2}) + @Parameter({1.2}) + public double theta; + + + @StringParameter({"sepSeeds"}) + public String startFromImm; + +// @Parameter({1.0}) +// public double infBa5; + + // BQ 1 +// @StringParameter({"off", "2.0", "2.25", "2.5", "2.75", "3.0"}) + @StringParameter({"off"}) //1.8 + public String escStrA; + + + @StringParameter({"2022-09-19"}) + public String strainADate; + // @StringParameter({"2022-08-24", "2022-08-29", "2022-09-04", "2022-09-09", "2022-09-14", "2022-09-19"}) + + // @StringParameter({"base", "2022-11-15", "2022-12-01", "2022-12-15", "2023-01-01"}) +// @StringParameter({"base"}) +// public String maskPt; + + //IFSG +// @StringParameter({"base"}) +// public String ifsg; + + // Vaccination Campaign +// @StringParameter({"base"}) +// String vacCamp; + + // NEW RESTRICTIONS +// @StringParameter({"base"}) +// public String work; + + // leisure Public +// @StringParameter({"base"}) +// public String leis; + + //edu +// @StringParameter({"base"}) +// public String edu; + + } + + + public static void main(String[] args) { + String[] args2 = { + RunParallel.OPTION_SETUP, CologneBMBF202302XX.class.getName(), + RunParallel.OPTION_PARAMS, Params.class.getName(), + RunParallel.OPTION_TASKS, Integer.toString(1), + RunParallel.OPTION_ITERATIONS, Integer.toString(150), + RunParallel.OPTION_METADATA + }; + + RunParallel.main(args2); + } + + +} + diff --git a/src/main/java/org/matsim/run/batch/CologneBMBF202310XX_soup.java b/src/main/java/org/matsim/run/batch/CologneBMBF202310XX_soup.java new file mode 100644 index 000000000..1c2eced77 --- /dev/null +++ b/src/main/java/org/matsim/run/batch/CologneBMBF202310XX_soup.java @@ -0,0 +1,1016 @@ +package org.matsim.run.batch; + +import com.google.inject.AbstractModule; +import com.google.inject.Module; +import com.google.inject.Singleton; +import com.google.inject.multibindings.Multibinder; +import com.google.inject.util.Modules; +import it.unimi.dsi.fastutil.ints.Int2DoubleAVLTreeMap; +import it.unimi.dsi.fastutil.ints.Int2DoubleMap; +import org.matsim.core.config.Config; +import org.matsim.core.config.ConfigUtils; +import org.matsim.episim.*; +import org.matsim.episim.analysis.*; +import org.matsim.episim.model.*; +import org.matsim.episim.model.listener.HouseholdSusceptibility; +import org.matsim.episim.model.vaccination.VaccinationModel; +import org.matsim.episim.model.vaccination.VaccinationStrategyBMBF202310XX; +import org.matsim.episim.policy.FixedPolicy; +import org.matsim.episim.policy.Restriction; +import org.matsim.run.RunParallel; +import org.matsim.run.modules.SnzCologneProductionScenario; +import org.matsim.run.modules.SnzProductionScenario; + +import javax.annotation.Nullable; +import java.io.IOException; +import java.io.UncheckedIOException; +import java.nio.file.Files; +import java.nio.file.Path; +import java.time.LocalDate; +import java.time.temporal.ChronoUnit; +import java.util.*; +import java.util.Map.Entry; +import java.util.stream.Collectors; +import java.util.stream.Stream; + + +/** + * Batch for Bmbf runs + */ +public class CologneBMBF202310XX_soup implements BatchRun { + + boolean DEBUG_MODE = false; + int runCount = 0; + + + @Nullable + @Override + public Module getBindings(int id, @Nullable Params params) { + return Modules.override(getBindings(0.0, params)).with(new AbstractModule() { + @Override + protected void configure() { + + // VACCINATION MODEL + + Multibinder set = Multibinder.newSetBinder(binder(), VaccinationModel.class); + + set.addBinding().to(VaccinationStrategyBMBF202310XX.class).in(Singleton.class); + + LocalDate start = LocalDate.parse("2023-10-01"); + + VaccinationType vaccinationType = VaccinationType.xbbUpdate; + + VaccinationStrategyBMBF202310XX.Config.VaccinationPool vaccinationPool = VaccinationStrategyBMBF202310XX.Config.VaccinationPool.vaccinated; + Int2DoubleMap compliance = new Int2DoubleAVLTreeMap(); + compliance.put(60, 0.0); // 60+ + compliance.put(18, 0.0); // 18 - 59 + compliance.put(12, 0.0); // 12 - 18 + compliance.put(0, 0.0); // 0 - 12 + + int minDaysAfterInfectionOrVaccination = 365; + + if (params != null) { + switch (params.vacCamp) { + case "base": + break; + case "pessimistic": + compliance.put(60, 0.19); + compliance.put(18, 0.); + compliance.put(12, 0.); + compliance.put(0, 0.); + break; + case "optimistic": + compliance.put(60, 0.54); + compliance.put(18, 0.); + compliance.put(12, 0.); + compliance.put(0, 0.); + break; + case "upperBound": + compliance.put(60, 1.0); + compliance.put(18, 0.); + compliance.put(12, 0.); + compliance.put(0, 0.); + vaccinationPool = VaccinationStrategyBMBF202310XX.Config.VaccinationPool.unvaccinated; + minDaysAfterInfectionOrVaccination = 182; + break; + default: + throw new RuntimeException("Not a valid option for vaccinationCampaignType"); + } + } + + + + bind(VaccinationStrategyBMBF202310XX.Config.class).toInstance(new VaccinationStrategyBMBF202310XX.Config(start, 90, vaccinationType, compliance, vaccinationPool, minDaysAfterInfectionOrVaccination, minDaysAfterInfectionOrVaccination)); + + + // ANTIBODY MODEL + // default values + double mutEscDelta = 29.2 / 10.9; + double mutEscBa1 = 10.9 / 1.9; + double mutEscBa5 = 5.0; + +// double mutEscStrainA = 1.0; +// double mutEscStrainB = 1.0; + double mutEscBqq = 5.0; // w/ respect to BA5 + double mutEscXBB_19 = 5.0; // wrt BA2 + double mutEscEG = 4.0; // wrt XBB_19 + double mutEscXBB_15 = 8.0; // wrt BA2 + + + double escape = 12.; + int days = 30; + String strainSeed = "no"; + LocalDate soupStartDate = LocalDate.parse("2020-01-01"); + boolean lineB = false; + double escapeBetweenLines = 1.0; + double hlMultiForInfected = 1.0; + + if (params != null) { +// mutEscBa1 = params.ba1Esc; +// mutEscBa5 = params.ba5Esc; + +// String StrainA = "6.0"; +// String StrainB = "off"; + + +// if (!params.StrainA.equals("off")) { +// mutEscStrainA = Double.parseDouble(params.StrainA); +// } +// if (!StrainB.equals("off")) { +// mutEscStrainB = Double.parseDouble(StrainB); +// } + escape = params.esc; + days = params.days; + strainSeed = params.strainRnd; + soupStartDate = LocalDate.parse(params.soupStartDate); + hlMultiForInfected = params.hlMultiForInfected; + mutEscBqq = params.escBqq; + mutEscEG = params.escEg; + + } + + //initial antibodies + Map> initialAntibodies = new HashMap<>(); + Map> antibodyRefreshFactors = new HashMap<>(); + configureAntibodies(initialAntibodies, antibodyRefreshFactors, mutEscDelta, mutEscBa1, mutEscBa5, mutEscBqq, mutEscXBB_15, mutEscXBB_19, mutEscEG, escape, days, strainSeed, soupStartDate, lineB, escapeBetweenLines); + + AntibodyModel.Config antibodyConfig = new AntibodyModel.Config(initialAntibodies, antibodyRefreshFactors, hlMultiForInfected); + + double immuneSigma = 3.0; + if (params != null) { + antibodyConfig.setImmuneReponseSigma(immuneSigma); + } + + bind(AntibodyModel.Config.class).toInstance(antibodyConfig); + + +// UtilsJR.printInitialAntibodiesToConsole(initialAntibodies, true); + + + if (params == null) return; + + // HOUSEHOLD SUSCEPTIBILITY + // designates a 35% of households as super safe; the susceptibility of that subpopulation is reduced to 1% wrt to general population. + bind(HouseholdSusceptibility.Config.class).toInstance( + HouseholdSusceptibility.newConfig() + .withSusceptibleHouseholds(0.35, 0.01) +// .withNonVaccinableHouseholds(params.nonVaccinableHh) +// .withShape(SnzCologneProductionScenario.INPUT.resolve("CologneDistricts.zip")) +// .withFeature("STT_NAME", vingst, altstadtNord, bickendorf, weiden) + ); + + } + + private void configureAntibodies(Map> initialAntibodies, + Map> antibodyRefreshFactors, + double mutEscDelta, + double mutEscBa1, + double mutEscBa5, + double mutEscBQ, + double mutEscXBB_15, + double mutEscXBB_19, + double mutEscEG, + double escapePerYear, + int days, + String strainSeed, + LocalDate soupStartDate, + boolean lineB, + double escapeBetweenLines) { + for (VaccinationType immunityType : VaccinationType.values()) { + initialAntibodies.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + + if (immunityType == VaccinationType.mRNA) { + initialAntibodies.get(immunityType).put(virusStrain, 29.2); //10.0 + } + else if (immunityType == VaccinationType.vector) { + initialAntibodies.get(immunityType).put(virusStrain, 6.8); //2.5 + } + else { + initialAntibodies.get(immunityType).put(virusStrain, 5.0); + } + } + } + + for (VirusStrain immunityType : VirusStrain.values()) { + initialAntibodies.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + initialAntibodies.get(immunityType).put(virusStrain, 5.0); + } + } + + //mRNAAlpha, mRNADelta, mRNABA1 comes from Sydney's calibration. + //The other values come from Rössler et al. + + //Wildtype + double mRNAAlpha = 29.2; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.SARS_CoV_2, mRNAAlpha); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.SARS_CoV_2, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.SARS_CoV_2, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.SARS_CoV_2, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.SARS_CoV_2, mRNAAlpha / mutEscDelta / mutEscBa1); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.SARS_CoV_2, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.SARS_CoV_2, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscXBB_15); + + + //Alpha + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.ALPHA, mRNAAlpha); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.ALPHA, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.ALPHA, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.ALPHA, mRNAAlpha * 300. / 700.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.ALPHA, mRNAAlpha * 210. / 700.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.ALPHA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.ALPHA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.ALPHA, mRNAAlpha / mutEscDelta / mutEscBa1); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.ALPHA, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.ALPHA, mRNAAlpha / mutEscDelta / mutEscBa1 / mutEscXBB_15); + + //DELTA + double mRNADelta = mRNAAlpha / mutEscDelta; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.DELTA, mRNADelta); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.DELTA, mRNADelta * 150./300.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.DELTA, mRNADelta * 64./300.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.DELTA, mRNADelta * 64./300.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.DELTA, mRNADelta * 450./300.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.DELTA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.DELTA, 0.01); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.DELTA, 0.01); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.DELTA, mRNAAlpha / mutEscBa1); // TODO: shouldn't this be mRNAALpha? + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.DELTA, mRNAAlpha / mutEscBa1 / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.DELTA, mRNAAlpha / mutEscBa1 / mutEscXBB_15); + + //BA.1 + double mRNABA1 = mRNADelta / mutEscBa1; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.OMICRON_BA1, mRNABA1); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.OMICRON_BA1, mRNABA1 * 4./20.); //??? + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.OMICRON_BA1, mRNABA1 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.OMICRON_BA1, mRNABA1 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.OMICRON_BA1, mRNABA1 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA1, 64.0 / 300.); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / 1.4); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5); //todo: is 1.4 +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscStrainB); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscBa5 / mutEscBQ); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscXBB_15); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscXBB_19); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscXBB_19 / mutEscEG); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA1, mRNAAlpha); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA1, mRNAAlpha / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.OMICRON_BA1, mRNAAlpha / mutEscXBB_15); + + //BA.2 + double mRNABA2 = mRNABA1; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.OMICRON_BA2, mRNABA2); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.OMICRON_BA2, mRNABA2 * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.OMICRON_BA2, mRNABA2 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.OMICRON_BA2, mRNABA2 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.OMICRON_BA2, mRNABA2 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / 1.4); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA2, 64.0 / 300.); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscStrainB); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscBa5 / mutEscBQ); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscXBB_15); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscXBB_19); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscXBB_19 / mutEscEG); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA2, mRNAAlpha); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA2, mRNAAlpha / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.OMICRON_BA2, mRNAAlpha / mutEscXBB_15); + + //BA.5 + double mRNABa5 = mRNABA2 / mutEscBa5; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.OMICRON_BA5, mRNABa5); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.OMICRON_BA5, mRNABa5 * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.OMICRON_BA5, mRNABa5 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.OMICRON_BA5, mRNABa5 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.OMICRON_BA5, mRNABa5 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBa5);// todo: do we need 1.4? + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscBa5); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.OMICRON_BA5, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscStrainA); //todo ??? +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBa5 / mutEscStrainB); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscBQ); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscXBB_15 / mutEscBa5); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscXBB_19 / mutEscBa5); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.OMICRON_BA5, 64.0 / 300. / mutEscEG / mutEscXBB_19 / mutEscBa5); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.OMICRON_BA5, mRNAAlpha / mutEscBa5); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.OMICRON_BA5, mRNAAlpha); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.OMICRON_BA5, mRNAAlpha / mutEscXBB_15 / mutEscBa5); +// +// //StrainA +// double mRNAStrainA = mRNABa5 / mutEscStrainA; +// initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.STRAIN_A, mRNAStrainA); +// initialAntibodies.get(VaccinationType.vector).put(VirusStrain.STRAIN_A, mRNAStrainA * 4./20.); +// initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); +// initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.STRAIN_A, mRNAStrainA * 6./20.); +// initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.STRAIN_A, mRNAStrainA * 8./20.); +// initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscBa5 /mutEscStrainA); +// initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.STRAIN_A, 64.0 / 300./ mutEscBa5 /mutEscStrainA); +// initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.STRAIN_A, 64.0 / 300.); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.STRAIN_A, 64.0 / 300. / mutEscStrainA / mutEscStrainB / mutEscBa5); +// +// initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscBa5 / mutEscStrainA); +// initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.STRAIN_A, mRNAAlpha / mutEscStrainA); +// initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.STRAIN_A, ... ) + + //XBB_15 + double mRNAXBB_15 = mRNABA2 / mutEscXBB_15; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.XBB_15, mRNAXBB_15); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.XBB_15, mRNAXBB_15 * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.XBB_15, mRNAXBB_15 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.XBB_15, mRNAXBB_15 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.XBB_15, mRNAXBB_15 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXBB_15); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.XBB_15, 64.0 / 300./ mutEscXBB_15); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXBB_15 / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.XBB_15, 64.0 / 300./ mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.XBB_15, 64.0 / 300.); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXBB_15 / mutEscBa5 / mutEscBQ); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.XBB_15, 64.0 / 300. ); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscXBB_19 / mutEscXBB_15); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscEG / mutEscXBB_19 / mutEscXBB_15); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.XBB_15, mRNAAlpha / mutEscXBB_15); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.XBB_15, mRNAAlpha / mutEscXBB_15 / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.XBB_15, mRNAAlpha); + + //XBB_19 + double mRNAXBB_19 = mRNABA2 / mutEscXBB_19; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.XBB_19, mRNAXBB_19); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.XBB_19, mRNAXBB_19 * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.XBB_19, mRNAXBB_19 * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.XBB_19, mRNAXBB_19 * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.XBB_19, mRNAXBB_19 * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.XBB_19, 64.0 / 300. / mutEscXBB_19); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.XBB_19, 64.0 / 300./ mutEscXBB_19); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.XBB_19, 64.0 / 300. / mutEscXBB_19 / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.XBB_19, 64.0 / 300./ mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.XBB_19, 64.0 / 300.); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.XBB_19, 64.0 / 300. / mutEscBQ / mutEscBa5 / mutEscXBB_19); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.XBB_19, 64.0 / 300. / mutEscXBB_15 / mutEscXBB_19); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.XBB_19, 64.0 / 300. ); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.XBB_19, 64.0 / 300. / mutEscEG); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.XBB_19, mRNAAlpha / mutEscXBB_19); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.XBB_19, mRNAAlpha / mutEscXBB_19 / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.XBB_19, mRNAAlpha / mutEscXBB_15 / mutEscXBB_19); + + //BQ + double mRNABQ = mRNABa5 / mutEscBQ; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.BQ, mRNABQ); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.BQ, mRNABQ * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.BQ, mRNABQ * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.BQ, mRNABQ * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.BQ, mRNABQ * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.BQ, 64.0 / 300. / mutEscBa5 / mutEscBQ); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.BQ, 64.0 / 300./ mutEscBa5 / mutEscBQ); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.BQ, 64.0 / 300. / mutEscBQ); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.BQ, 64.0 / 300./ mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.BQ, 64.0 / 300.); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.BQ, 64.0 / 300. ); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.BQ, 64.0 / 300. / mutEscXBB_15 / mutEscBa5 / mutEscBQ); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.BQ, 64.0 / 300. / mutEscXBB_19 /mutEscBa5 / mutEscBQ); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.BQ, 64.0 / 300. / mutEscEG / mutEscXBB_19 / mutEscBa5 / mutEscBQ); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.BQ, mRNAAlpha / mutEscBa5 / mutEscBQ); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.BQ, mRNAAlpha / mutEscBQ); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.BQ, mRNAAlpha / mutEscXBB_15 / mutEscBa5 / mutEscBQ); + + //EG + double mRNAEG = mRNAXBB_19 / mutEscEG; + initialAntibodies.get(VaccinationType.mRNA).put(VirusStrain.EG, mRNAEG); + initialAntibodies.get(VaccinationType.vector).put(VirusStrain.EG, mRNAEG * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(VirusStrain.EG, mRNAEG * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(VirusStrain.EG, mRNAEG * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(VirusStrain.EG, mRNAEG * 8./20.); + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(VirusStrain.EG, 64.0 / 300. / mutEscXBB_19 / mutEscEG); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(VirusStrain.EG, 64.0 / 300./ mutEscXBB_19 / mutEscEG); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(VirusStrain.EG, 64.0 / 300. / mutEscXBB_19 / mutEscEG / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_A).put(VirusStrain.EG, 64.0 / 300./ mutEscStrainA / mutEscStrainB / mutEscBa5); +// initialAntibodies.get(VirusStrain.STRAIN_B).put(VirusStrain.EG, 64.0 / 300.); + initialAntibodies.get(VirusStrain.BQ).put(VirusStrain.EG, 64.0 / 300. / mutEscBQ / mutEscBa5 / mutEscXBB_19 / mutEscEG); + initialAntibodies.get(VirusStrain.XBB_15).put(VirusStrain.EG, 64.0 / 300. / mutEscXBB_15 / mutEscXBB_19 / mutEscEG); + initialAntibodies.get(VirusStrain.XBB_19).put(VirusStrain.EG, 64.0 / 300. / mutEscEG); + initialAntibodies.get(VirusStrain.EG).put(VirusStrain.EG, 64.0 / 300. ); + initialAntibodies.get(VaccinationType.ba1Update).put(VirusStrain.EG, mRNAAlpha / mutEscXBB_19 / mutEscEG); + initialAntibodies.get(VaccinationType.ba5Update).put(VirusStrain.EG, mRNAAlpha / mutEscXBB_19 / mutEscEG / mutEscBa5); + initialAntibodies.get(VaccinationType.xbbUpdate).put(VirusStrain.EG, mRNAAlpha / mutEscXBB_15 / mutEscXBB_19 / mutEscEG); + + // strains A_1, A_2, ... + { + + ArrayList strains = getNewStrains(Boolean.valueOf(lineB)); + + ArrayList dates = getDatesNewStrains(strains, days, strainSeed, soupStartDate); + + for (int i = 0; i < strains.size(); i++) { + long daysSince = ChronoUnit.DAYS.between(soupStartDate, dates.get(i)); + double escape = 1. + (escapePerYear - 1.0) * daysSince / 365.; //factor 6, if variant appears 6 months later + VirusStrain strain = strains.get(i); + + initialAntibodies.get(strain).put(VirusStrain.SARS_CoV_2, 0.01); + initialAntibodies.get(strain).put(VirusStrain.ALPHA, 0.01); + initialAntibodies.get(strain).put(VirusStrain.DELTA, 0.01); + + initialAntibodies.get(strain).put(VirusStrain.OMICRON_BA1, 64.0 / 300. / mutEscXBB_19 / mutEscEG / escape); + initialAntibodies.get(strain).put(VirusStrain.OMICRON_BA2, 64.0 / 300. / mutEscXBB_19 / mutEscEG / escape); + initialAntibodies.get(strain).put(VirusStrain.OMICRON_BA5, 64.0 / 300./ mutEscBa5 / mutEscXBB_19 / mutEscEG / escape); + + initialAntibodies.get(strain).put(VirusStrain.XBB_15, 64.0 / 300. / mutEscEG / mutEscXBB_19 / mutEscXBB_15 / escape); + initialAntibodies.get(strain).put(VirusStrain.XBB_19, 64.0 / 300. / mutEscEG / escape); + initialAntibodies.get(strain).put(VirusStrain.BQ, 64.0 / 300. / mutEscEG / mutEscXBB_19 / mutEscBa5 / mutEscBQ / escape); + initialAntibodies.get(strain).put(VirusStrain.EG, 64.0 / 300. / escape); + +// initialAntibodies.get(strain).put(VirusStrain.STRAIN_A, 64.0 / 300. / escape); + + double mRNAStrain = mRNAEG / escape; + initialAntibodies.get(VaccinationType.mRNA).put(strain, mRNAStrain); + initialAntibodies.get(VaccinationType.vector).put(strain, mRNAStrain * 4./20.); + initialAntibodies.get(VirusStrain.SARS_CoV_2).put(strain, mRNAStrain * 6./20.); + initialAntibodies.get(VirusStrain.ALPHA).put(strain, mRNAStrain * 6./20.); + initialAntibodies.get(VirusStrain.DELTA).put(strain, mRNAStrain * 8./20.); + + initialAntibodies.get(VirusStrain.OMICRON_BA1).put(strain, 64.0 / 300. / mutEscXBB_19 / mutEscEG / escape); + initialAntibodies.get(VirusStrain.OMICRON_BA2).put(strain, 64.0 / 300./ mutEscXBB_19 / mutEscEG / escape); + initialAntibodies.get(VirusStrain.OMICRON_BA5).put(strain, 64.0 / 300. / mutEscXBB_19 / mutEscEG / mutEscBa5 / escape); + + initialAntibodies.get(VirusStrain.BQ).put(strain, 64.0 / 300. / mutEscBQ / mutEscBa5 / mutEscXBB_19 / mutEscEG / escape); + initialAntibodies.get(VirusStrain.XBB_15).put(strain, 64.0 / 300. / mutEscXBB_15 / mutEscXBB_19 / mutEscEG / escape); + initialAntibodies.get(VirusStrain.XBB_19).put(strain, 64.0 / 300. / mutEscEG / escape); + initialAntibodies.get(VirusStrain.EG).put(strain, 64.0 / 300. / escape); + initialAntibodies.get(VaccinationType.ba1Update).put(strain, mRNAAlpha / mutEscXBB_19 / mutEscEG / escape); + initialAntibodies.get(VaccinationType.ba5Update).put(strain, mRNAAlpha / mutEscXBB_19 / mutEscEG / mutEscBa5 / escape); + initialAntibodies.get(VaccinationType.xbbUpdate).put(strain, mRNAAlpha / mutEscXBB_15 / mutEscXBB_19 / mutEscEG / escape); + + for (int j = 0; j < strains.size(); j++) { + LocalDate date1 = dates.get(i); + LocalDate date2 = dates.get(j); + long daysBetweenStrains = Math.abs(ChronoUnit.DAYS.between(date1, date2)); + double escapeBetweenStrains = 1. + (escapePerYear - 1.0) * daysBetweenStrains / 365.; //factor 6, if variant appears 6 months later + VirusStrain strain2 = strains.get(j); + + if (strain.toString().charAt(0) != strain2.toString().charAt(0)) + escapeBetweenStrains = escapeBetweenStrains * escapeBetweenLines; + + initialAntibodies.get(strain).put(strain2, 64.0 / 300. / escapeBetweenStrains); + } + + } + + } + + + + for (VaccinationType immunityType : VaccinationType.values()) { + antibodyRefreshFactors.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + + if (immunityType == VaccinationType.mRNA) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } else if (immunityType == VaccinationType.vector) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 5.0); + } else if (immunityType == VaccinationType.ba1Update) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } else if (immunityType == VaccinationType.ba5Update) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } else if (immunityType == VaccinationType.xbbUpdate) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } else { + antibodyRefreshFactors.get(immunityType).put(virusStrain, Double.NaN); + } + + } + } + + for (VirusStrain immunityType : VirusStrain.values()) { + antibodyRefreshFactors.put(immunityType, new EnumMap<>( VirusStrain.class ) ); + for (VirusStrain virusStrain : VirusStrain.values()) { + antibodyRefreshFactors.get(immunityType).put(virusStrain, 15.0); + } + } + + +// UtilsJR.printInitialAntibodiesToConsole(initialAntibodies, false); + + } + }); + + } + + private SnzCologneProductionScenario getBindings(double pHousehold, Params params) { + return new SnzCologneProductionScenario.Builder() + .setCarnivalModel(SnzCologneProductionScenario.CarnivalModel.yes) + .setFutureVacations(SnzCologneProductionScenario.FutureVacations.yes)//params != null ? params.futureVacations : SnzCologneProductionScenario.FutureVacations.no) + .setSebastianUpdate(true) + .setLeisureCorrection(1.3) //params == null ? 0.0 : params.actCorrection) + .setScaleForActivityLevels(1.3) + .setSuscHouseholds_pct(pHousehold) + .setActivityHandling(EpisimConfigGroup.ActivityHandling.startOfDay) +// .setTestingModel(params != null ? FlexibleTestingModel.class : DefaultTestingModel.class) + .setInfectionModel(InfectionModelWithAntibodies.class) + .build(); + } + + @Override + public Metadata getMetadata() { + return Metadata.of("cologne", "calibration"); + } + + @Override + public Collection postProcessing() { + return List.of( +// new VaccinationEffectiveness().withArgs(), +// new RValuesFromEvents().withArgs(), +// new VaccinationEffectivenessFromPotentialInfections().withArgs("--remove-infected"), +// new FilterEvents().withArgs("--output","./output/"), + new HospitalNumbersFromEvents().withArgs("--output","./output/","--input","/scratch/projects/bzz0020/episim-input") +// new SecondaryAttackRateFromEvents().withArgs() + ); + } + + @Override + public Config prepareConfig(int id, Params params) { + + if (DEBUG_MODE) { + if (runCount == 0) { //&& params.strAEsc != 0.0 && params.ba5Inf == 0. && params.eduTest.equals("true")) { + runCount++; + } else { + return null; + } + } + + SnzCologneProductionScenario module = getBindings(0.0, params); + + Config config = module.config(); + + config.global().setRandomSeed(params.seed); + + EpisimConfigGroup episimConfig = ConfigUtils.addOrGetModule(config, EpisimConfigGroup.class); + + episimConfig.setCalibrationParameter(episimConfig.getCalibrationParameter() * 1.2 * 1.7); + + //--------------------------------------- + // S N A P S H O T + //--------------------------------------- + + // start from immunization history + String directoryForImmHist = "/scratch/projects/bzz0020/runs/jakob/2023-10-31/1-bmbf/output/seed_" + params.seed + "-TmidFall2022_" + params.TmidFall2022 + "-vacCamp_base-soupStartDate_2023-09-01-esc_12.0-days_30-strainRnd_no-lineB_false-iga_true-seasonal_true-rf2023_base-hlMultiForInfected_" + params.hlMultiForInfected + "-escBqq_2.0/"; + episimConfig.setStartFromImmunization(directoryForImmHist); + episimConfig.setStartDate(LocalDate.parse("2023-04-01")); + + //--------------------------------------- + // S T R A I N S + //--------------------------------------- + + VirusStrainConfigGroup virusStrainConfigGroup = ConfigUtils.addOrGetModule(config, VirusStrainConfigGroup.class); + + // BQ (from BA5) + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setInfectiousness(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA5).getInfectiousness()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorSeriouslySick(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA5).getFactorSeriouslySick()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorSeriouslySickVaccinated(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA5).getFactorSeriouslySickVaccinated()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.BQ).setFactorCritical(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA5).getFactorCritical()); + + // XBB 1.5 (from BA2) + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setInfectiousness(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getInfectiousness()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorSeriouslySick(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getFactorSeriouslySick()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorSeriouslySickVaccinated(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getFactorSeriouslySickVaccinated()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_15).setFactorCritical(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getFactorCritical()); + + // XBB 1.9 (from BA2) + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_19).setInfectiousness(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getInfectiousness()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_19).setFactorSeriouslySick(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getFactorSeriouslySick()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_19).setFactorSeriouslySickVaccinated(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getFactorSeriouslySickVaccinated()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.XBB_19).setFactorCritical(virusStrainConfigGroup.getParams(VirusStrain.OMICRON_BA2).getFactorCritical()); + + + // EG + virusStrainConfigGroup.getOrAddParams(VirusStrain.EG).setInfectiousness(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getInfectiousness()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.EG).setFactorSeriouslySick(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getFactorSeriouslySick()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.EG).setFactorSeriouslySickVaccinated(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getFactorSeriouslySickVaccinated()); + virusStrainConfigGroup.getOrAddParams(VirusStrain.EG).setFactorCritical(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getFactorCritical()); + + +// double ba5Inf = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA5).getInfectiousness(); +// double ba5Hos = virusStrainConfigGroup.getOrAddParams(VirusStrain.OMICRON_BA5).getFactorSeriouslySick(); + +// STRAIN_A +// if (!params.StrainA.equals("off")) { +// +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setInfectiousness(ba5Inf); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySick(ba5Hos); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorSeriouslySickVaccinated(ba5Hos); +// virusStrainConfigGroup.getOrAddParams(VirusStrain.STRAIN_A).setFactorCritical(ba5Hos); +// } + + for (VirusStrain strain : getNewStrains(Boolean.valueOf(params.lineB))) { + virusStrainConfigGroup.getOrAddParams(strain).setInfectiousness(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getInfectiousness()); + virusStrainConfigGroup.getOrAddParams(strain).setFactorSeriouslySick(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getFactorSeriouslySick()); + virusStrainConfigGroup.getOrAddParams(strain).setFactorSeriouslySickVaccinated(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getFactorSeriouslySickVaccinated()); + virusStrainConfigGroup.getOrAddParams(strain).setFactorCritical(virusStrainConfigGroup.getParams(VirusStrain.XBB_19).getFactorCritical()); + } + + //--------------------------------------- + // I M P O R T + //--------------------------------------- + + configureFutureDiseaseImport(params, episimConfig); + + // modify import: + LocalDate impModDate = LocalDate.parse("2022-01-31"); + double impRedBa1 = 0.0; + double impRedBa2 = 0.0; + if (impRedBa1 != 1.0) { + NavigableMap impBa1 = episimConfig.getInfections_pers_per_day().get(VirusStrain.OMICRON_BA1); + for (Entry entry : impBa1.entrySet()) { + if (entry.getKey().isAfter(impModDate)) { + impBa1.put(entry.getKey(), (int) (entry.getValue() * impRedBa1)); + } + } + } + + if (impRedBa2 != 1.0) { + NavigableMap impBa2 = episimConfig.getInfections_pers_per_day().get(VirusStrain.OMICRON_BA2); + for (Entry entry : impBa2.entrySet()) { + if (entry.getKey().isAfter(impModDate)) { + impBa2.put(entry.getKey(), (int) (entry.getValue() * impRedBa2)); + } + } + } + + + + //--------------------------------------- + // R E S T R I C T I O N S + //--------------------------------------- + + FixedPolicy.ConfigBuilder builder = FixedPolicy.parse(episimConfig.getPolicy()); + + + // vary amount of "school" activity that takes place during summer vacation 2022 + builder.restrict(LocalDate.parse("2022-06-27"), 0.8, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); + + // vary amount of "school" activity that takes place during summer vacation 2023 + builder.restrict(LocalDate.parse("2023-06-22"), params == null ? 0.2 : params.schoolVac, "educ_primary", "educ_kiga", "educ_secondary", "educ_tertiary", "educ_other"); + + // increase work, leisure Rf after SNZ Data runs out + if (params.rf2023.equals("base")) { + builder.restrict(LocalDate.parse("2023-01-01"), 0.73, "work", "leisure", "leisPublic", "leisPrivate"); + } else if (params.rf2023.equals("leis")) { + builder.restrict(LocalDate.parse("2023-01-01"), 0.73, "work"); + builder.interpolate(LocalDate.parse("2023-01-01"), LocalDate.parse("2023-03-31"), Restriction.of(0.73), Restriction.of(1.0), "leisure", "leisPublic", "leisPrivate"); + } else if (params.rf2023.equals("work_leis")) { + builder.interpolate(LocalDate.parse("2023-01-01"), LocalDate.parse("2023-03-31"), Restriction.of(0.73), Restriction.of(1.0), "work", "leisure", "leisPublic", "leisPrivate"); + } else { + throw new RuntimeException(); + } + + + episimConfig.setPolicy(builder.build()); + + + //--------------------------------------- + // M I S C + //--------------------------------------- + + + // weather + double maxOutdoorFraction = 0.8; + double midpoint1 = 18.5; + double midpoint2 = 25.0; + double TmidFall2022 = 18.5; + + if (params != null) { + TmidFall2022 = params.TmidFall2022; + } + + try { + Map outdoorFractions = EpisimUtils.getOutDoorFractionFromDateAndTemp2Fall2022Override( + SnzCologneProductionScenario.INPUT.resolve("cologneWeather.csv").toFile(), + SnzCologneProductionScenario.INPUT.resolve("weatherDataAvgCologne2000-2020.csv").toFile(), + 0.5, + midpoint1, + midpoint2, + midpoint1, + midpoint1, + TmidFall2022, + 5., + 1.0, + maxOutdoorFraction); + + episimConfig.setLeisureOutdoorFraction(outdoorFractions); + } catch (IOException e) { + throw new UncheckedIOException(e); + } + + // add vaccination + + VaccinationConfigGroup vaccinationConfig = ConfigUtils.addOrGetModule(config, VaccinationConfigGroup.class); + vaccinationConfig.setUseIgA(Boolean.valueOf(params.iga)); + + vaccinationConfig.getOrAddParams(VaccinationType.xbbUpdate); + + if (!Boolean.valueOf(params.seasonal)) { + + Map fractionsOld = episimConfig.getLeisureOutdoorFraction(); + Map fractionsNew = new HashMap(); + + for (Entry e : fractionsOld.entrySet()) { + if (e.getKey().isBefore(LocalDate.parse("2022-12-01"))) + fractionsNew.put(e.getKey(), e.getValue()); + + } + episimConfig.setLeisureOutdoorFraction(fractionsNew); + } + + + if (DEBUG_MODE) { + UtilsJR.produceDiseaseImportPlot(episimConfig.getInfections_pers_per_day()); + } + + return config; + } + + private void configureFutureDiseaseImport(Params params, EpisimConfigGroup episimConfig) { + Map infPerDayBa1 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA1, new TreeMap<>())); + Map infPerDayBa2 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA2, new TreeMap<>())); + Map infPerDayBa5 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.OMICRON_BA5, new TreeMap<>())); + + + Map infPerDayXbb15 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.XBB_15, new TreeMap<>())); + Map infPerDayXbb19 = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.XBB_19, new TreeMap<>())); + Map infPerDayEg = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.EG, new TreeMap<>())); + Map infPerDayBq = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.BQ, new TreeMap<>())); + +// Map infPerDayStrA = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.STRAIN_A, new TreeMap<>())); +// Map infPerDayStrB = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(VirusStrain.STRAIN_B, new TreeMap<>())); + + //StrainA +// if (!params.StrainA.equals("off")) { +// infPerDayStrA.put(LocalDate.parse("2020-01-01"), 0); +//// LocalDate strADate = LocalDate.parse("2022-11-01"); +// LocalDate strADate = LocalDate.parse(params.strainADate); +// +// for (int i = 0; i < 7; i++) { +// infPerDayStrA.put(strADate.plusDays(i), 4); +// } +// infPerDayStrA.put(strADate.plusDays(7), 1); +// } + + // BQQ + LocalDate bqDate = LocalDate.parse("2022-10-01"); + infPerDayBq.put(LocalDate.parse("2020-01-01"), 0); + for (int j = 0; j < 7; j++) { + infPerDayBq.put(bqDate.plusDays(j), 4); + } + infPerDayBq.put(bqDate.plusDays(7), 1); + episimConfig.setInfections_pers_per_day(VirusStrain.BQ, infPerDayBq); + + // XBB15 + LocalDate xbb15Date = LocalDate.parse("2023-01-01"); + infPerDayXbb15.put(LocalDate.parse("2020-01-01"), 0); + for (int j = 0; j < 7; j++) { + infPerDayXbb15.put(xbb15Date.plusDays(j), 4); + } + infPerDayXbb15.put(xbb15Date.plusDays(7), 1); + episimConfig.setInfections_pers_per_day(VirusStrain.XBB_15, infPerDayXbb15); + + + LocalDate xbb19Date = LocalDate.parse("2023-03-01"); + infPerDayXbb19.put(LocalDate.parse("2020-01-01"), 0); + for (int j = 0; j < 7; j++) { + infPerDayXbb19.put(xbb19Date.plusDays(j), 4); + } + infPerDayXbb19.put(xbb19Date.plusDays(7), 1); + episimConfig.setInfections_pers_per_day(VirusStrain.XBB_19, infPerDayXbb19); + + + LocalDate egDate = LocalDate.parse("2023-06-01"); + infPerDayEg.put(LocalDate.parse("2020-01-01"), 0); + + if (params != null) { + egDate = LocalDate.parse(params.startDateEg); + } + + for (int j = 0; j < 7; j++) { + infPerDayEg.put(egDate.plusDays(j), 4); + } + infPerDayEg.put(egDate.plusDays(7), 1); + episimConfig.setInfections_pers_per_day(VirusStrain.EG, infPerDayEg); + + + for (int i = 0; i < getNewStrains(Boolean.valueOf(params.lineB)).size(); i++) { + LocalDate date = getDatesNewStrains(getNewStrains(Boolean.valueOf(params.lineB)), params.days, params.strainRnd, LocalDate.parse(params.soupStartDate)).get(i); + VirusStrain strain = getNewStrains(Boolean.valueOf(params.lineB)).get(i); + + Map infPerDayStrainX = new HashMap<>(episimConfig.getInfections_pers_per_day().getOrDefault(strain, new TreeMap<>())); + infPerDayStrainX.put(LocalDate.parse("2020-01-01"), 0); + for (int j = 0; j < 7; j++) { + infPerDayStrainX.put(date.plusDays(j), 4); + } + infPerDayStrainX.put(date.plusDays(7), 1); + episimConfig.setInfections_pers_per_day(strain, infPerDayStrainX); + } + + // save disease import + episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA1, infPerDayBa1); + episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA2, infPerDayBa2); + episimConfig.setInfections_pers_per_day(VirusStrain.OMICRON_BA5, infPerDayBa5); + +// if (!params.StrainA.equals("off")) { +// episimConfig.setInfections_pers_per_day(VirusStrain.STRAIN_A, infPerDayStrA); +// } +// i f (!params.StrainB.equals("off")) { +// episimConfig.setInfections_pers_per_day(VirusStrain.STRAIN_B, infPerDayStrB); +// } + } + + public static final class Params { + // general + @GenerateSeeds(5) + public long seed; + + @Parameter({18.5, 22.0}) + public double TmidFall2022; + + // future vacations +// @EnumParameter(SnzCologneProductionScenario.FutureVacations.class) +// public SnzCologneProductionScenario.FutureVacations futureVacations; + + // Vaccination Campaign + @StringParameter({"base", "upperBound"}) + String vacCamp; + + @StringParameter({"2023-09-01"}) + public String soupStartDate; + + @Parameter({12., 24.}) + public double esc; + + @IntParameter({30}) + public int days; + + @StringParameter({"no"}) + public String strainRnd; + + @StringParameter({"false"}) + public String lineB; + + @StringParameter({"true"}) + public String iga; + + @StringParameter({"true"}) + public String seasonal; + + @StringParameter({"base"}) + public String rf2023; + + @Parameter({1., 1.5, 2.}) + public double hlMultiForInfected; + + @Parameter({2.}) + public double escBqq; + + @Parameter({0.2, 0.8}) + public double schoolVac; + + @Parameter({4.0, 6.0, 8.0}) // 3 + public double escEg; + + @StringParameter({"2023-06-01", "2023-05-15", "2023-05-01"}) // 3 + public String startDateEg; + } + + + public static void main(String[] args) { + String[] args2 = { + RunParallel.OPTION_SETUP, CologneBMBF202310XX_soup.class.getName(), + RunParallel.OPTION_PARAMS, Params.class.getName(), + RunParallel.OPTION_TASKS, Integer.toString(1), + RunParallel.OPTION_ITERATIONS, Integer.toString(1000), + RunParallel.OPTION_METADATA + }; + + RunParallel.main(args2); + } + + private static ArrayList getNewStrains(boolean lineB) { + ArrayList strains = new ArrayList(); + strains.add(VirusStrain.A_1); + strains.add(VirusStrain.A_2); + strains.add(VirusStrain.A_3); + strains.add(VirusStrain.A_4); + strains.add(VirusStrain.A_5); + strains.add(VirusStrain.A_6); + strains.add(VirusStrain.A_7); + strains.add(VirusStrain.A_8); + strains.add(VirusStrain.A_9); + strains.add(VirusStrain.A_10); + strains.add(VirusStrain.A_11); + strains.add(VirusStrain.A_12); + strains.add(VirusStrain.A_13); + strains.add(VirusStrain.A_14); + strains.add(VirusStrain.A_15); + strains.add(VirusStrain.A_16); + strains.add(VirusStrain.A_17); + strains.add(VirusStrain.A_18); + strains.add(VirusStrain.A_19); + strains.add(VirusStrain.A_20); + + + if (lineB) { + strains.add(1, VirusStrain.B_1); + strains.add(3, VirusStrain.B_2); + strains.add(5, VirusStrain.B_3); + strains.add(7, VirusStrain.B_4); + strains.add(9, VirusStrain.B_5); + strains.add(11, VirusStrain.B_6); + strains.add(13, VirusStrain.B_7); + strains.add(15, VirusStrain.B_8); + strains.add(17, VirusStrain.B_9); + strains.add(19, VirusStrain.B_10); + strains.add(21, VirusStrain.B_11); + strains.add(23, VirusStrain.B_12); + strains.add(25, VirusStrain.B_13); + strains.add(27, VirusStrain.B_14); + strains.add(29, VirusStrain.B_15); + strains.add(31, VirusStrain.B_16); + strains.add(33, VirusStrain.B_17); + strains.add(35, VirusStrain.B_18); + strains.add(37, VirusStrain.B_19); + strains.add(39, VirusStrain.B_20); + } + + return strains; + } + + private static ArrayList getDatesNewStrains(ArrayList strains, int days, String seed, LocalDate start) { + ArrayList dates = new ArrayList(); + + if (seed.equals("no")) { + for (LocalDate date = start; ; date = date.plusDays(1)) { + long daysBetween = ChronoUnit.DAYS.between(start, date); + if (daysBetween % days == 0) + dates.add(date); + if (dates.size() == strains.size()) + break; + } + return dates; + } + + else { + Random rand = new Random(Integer.parseInt(seed)); + for (LocalDate date = LocalDate.parse("2022-11-15"); ; date = date.plusDays(1)) { + if (rand.nextDouble() < 1. / days) + dates.add(date); + if (dates.size() == strains.size()) + break; + } + return dates; + } + + } +} + diff --git a/src/main/java/org/matsim/run/batch/StarterBatchBerlin.java b/src/main/java/org/matsim/run/batch/StarterBatchBerlin.java new file mode 100644 index 000000000..0dbcc9778 --- /dev/null +++ b/src/main/java/org/matsim/run/batch/StarterBatchBerlin.java @@ -0,0 +1,120 @@ +package org.matsim.run.batch; + +import com.google.inject.Module; +import org.matsim.core.config.Config; +import org.matsim.core.config.ConfigUtils; +import org.matsim.episim.BatchRun; +import org.matsim.episim.EpisimConfigGroup; +import org.matsim.episim.VirusStrainConfigGroup; +import org.matsim.episim.analysis.OutputAnalysis; +import org.matsim.episim.model.InfectionModelWithAntibodies; +import org.matsim.run.RunParallel; +import org.matsim.run.modules.SnzBerlinProductionScenario; +import org.matsim.run.modules.SnzCologneProductionScenario; + +import javax.annotation.Nullable; +import java.util.Collection; +import java.util.List; + + +/** + * boilerplate batch for berlin + */ +public class StarterBatchBerlin implements BatchRun { + + /* + * here you can swap out vaccination model, antibody model, etc. + * See CologneBMBF202310XX_soup.java for an example + */ + @Nullable + @Override + public Module getBindings(int id, @Nullable Params params) { + return getBindings(params); + } + + + /* + * here you select & modify models specified in the SnzCologneProductionScenario & SnzProductionScenario. + */ + private SnzBerlinProductionScenario getBindings(Params params) { + return new SnzBerlinProductionScenario.Builder() + .setActivityHandling(EpisimConfigGroup.ActivityHandling.startOfDay) + .setInfectionModel(InfectionModelWithAntibodies.class) + .setEasterModel(SnzBerlinProductionScenario.EasterModel.no) + .setChristmasModel(SnzBerlinProductionScenario.ChristmasModel.no) + .setSample(1) + .build(); + } + + /* + * Metadata is needed for covid-sim. + */ + @Override + public Metadata getMetadata() { + return Metadata.of("berlin", "calibration"); + } + + + /* + * Here you can add post-processing classes, that will be executed after the simulation. + */ + @Override + public Collection postProcessing() { + return List.of(); + } + + /* + * Here you can specify configuration options + */ + @Override + public Config prepareConfig(int id, Params params) { + + // Level 1: General (matsim) config. Here you can specify number of iterations and the seed. + Config config = getBindings(params).config(); + + config.global().setRandomSeed(params.seed); + + // Level 2: Episim specific configs: + // 2a: general episim config + EpisimConfigGroup episimConfig = ConfigUtils.addOrGetModule(config, EpisimConfigGroup.class); + + episimConfig.setCalibrationParameter(episimConfig.getCalibrationParameter() * 1.2 * 1.7 * params.thetaFactor); + + // 2b: specific config groups, e.g. virusStrainConfigGroup + VirusStrainConfigGroup virusStrainConfigGroup = ConfigUtils.addOrGetModule(config, VirusStrainConfigGroup.class); + + return config; + } + + + /* + * Specify parameter combinations that will be run. + */ + public static final class Params { + // general + @GenerateSeeds(1) + public long seed; + + @Parameter({1.0, 2.0}) + public double thetaFactor; + } + + + + /* + * top-level parameters for a run on your local machine. + */ + public static void main(String[] args) { + String[] args2 = { + RunParallel.OPTION_SETUP, StarterBatchBerlin.class.getName(), + RunParallel.OPTION_PARAMS, Params.class.getName(), + RunParallel.OPTION_TASKS, Integer.toString(8), + RunParallel.OPTION_ITERATIONS, Integer.toString(50), + RunParallel.OPTION_METADATA + }; + + RunParallel.main(args2); + } + +} + diff --git a/src/main/java/org/matsim/run/batch/StarterBatchCologne.java b/src/main/java/org/matsim/run/batch/StarterBatchCologne.java new file mode 100644 index 000000000..173c6d563 --- /dev/null +++ b/src/main/java/org/matsim/run/batch/StarterBatchCologne.java @@ -0,0 +1,136 @@ +package org.matsim.run.batch; + +import com.google.inject.AbstractModule; +import com.google.inject.Module; +import com.google.inject.Singleton; +import com.google.inject.multibindings.Multibinder; +import com.google.inject.util.Modules; +import it.unimi.dsi.fastutil.ints.Int2DoubleAVLTreeMap; +import it.unimi.dsi.fastutil.ints.Int2DoubleMap; +import org.matsim.core.config.Config; +import org.matsim.core.config.ConfigUtils; +import org.matsim.episim.*; +import org.matsim.episim.analysis.HospitalNumbersFromEvents; +import org.matsim.episim.analysis.OutputAnalysis; +import org.matsim.episim.model.*; +import org.matsim.episim.model.listener.HouseholdSusceptibility; +import org.matsim.episim.model.vaccination.VaccinationModel; +import org.matsim.episim.model.vaccination.VaccinationStrategyBMBF202310XX; +import org.matsim.episim.policy.FixedPolicy; +import org.matsim.episim.policy.Restriction; +import org.matsim.run.RunParallel; +import org.matsim.run.modules.SnzCologneProductionScenario; + +import javax.annotation.Nullable; +import java.io.IOException; +import java.io.UncheckedIOException; +import java.time.LocalDate; +import java.time.temporal.ChronoUnit; +import java.util.*; +import java.util.Map.Entry; + + +/** + * boilerplate batch for cologne + */ +public class StarterBatchCologne implements BatchRun { + + /* + * here you can swap out vaccination model, antibody model, etc. + * See CologneBMBF202310XX_soup.java for an example + */ + @Nullable + @Override + public Module getBindings(int id, @Nullable Params params) { + return getBindings(params); + } + + + /* + * here you select & modify models specified in the SnzCologneProductionScenario & SnzProductionScenario. + */ + private SnzCologneProductionScenario getBindings(Params params) { + return new SnzCologneProductionScenario.Builder() + .setCarnivalModel(SnzCologneProductionScenario.CarnivalModel.yes) + .setFutureVacations(SnzCologneProductionScenario.FutureVacations.yes) + .setSebastianUpdate(true) + .setLeisureCorrection(1.3) //params == null ? 0.0 : params.actCorrection) + .setScaleForActivityLevels(1.3) + .setSuscHouseholds_pct(0.35) + .setActivityHandling(EpisimConfigGroup.ActivityHandling.startOfDay) + .setInfectionModel(InfectionModelWithAntibodies.class) + .build(); + } + + /* + * Metadata is needed for covid-sim. + */ + @Override + public Metadata getMetadata() { + return Metadata.of("cologne", "calibration"); + } + + + /* + * Here you can add post-processing classes, that will be executed after the simulation. + */ + @Override + public Collection postProcessing() { + return List.of(); + } + + /* + * Here you can specify configuration options + */ + @Override + public Config prepareConfig(int id, Params params) { + + // Level 1: General (matsim) config. Here you can specify number of iterations and the seed. + Config config = getBindings(params).config(); + + config.global().setRandomSeed(params.seed); + + // Level 2: Episim specific configs: + // 2a: general episim config + EpisimConfigGroup episimConfig = ConfigUtils.addOrGetModule(config, EpisimConfigGroup.class); + + episimConfig.setCalibrationParameter(episimConfig.getCalibrationParameter() * 1.2 * 1.7 * params.thetaFactor); + + // 2b: specific config groups, e.g. virusStrainConfigGroup + VirusStrainConfigGroup virusStrainConfigGroup = ConfigUtils.addOrGetModule(config, VirusStrainConfigGroup.class); + + return config; + } + + + /* + * Specify parameter combinations that will be run. + */ + public static final class Params { + // general + @GenerateSeeds(1) + public long seed; + + @Parameter({1.0, 2.0}) + public double thetaFactor; + } + + + + /* + * top-level parameters for a run on your local machine. + */ + public static void main(String[] args) { + String[] args2 = { + RunParallel.OPTION_SETUP, StarterBatchCologne.class.getName(), + RunParallel.OPTION_PARAMS, Params.class.getName(), + RunParallel.OPTION_TASKS, Integer.toString(8), + RunParallel.OPTION_ITERATIONS, Integer.toString(10), + RunParallel.OPTION_METADATA + }; + + RunParallel.main(args2); + } + +} + diff --git a/src/main/java/org/matsim/run/batch/UtilsJR.java b/src/main/java/org/matsim/run/batch/UtilsJR.java index 9883e471d..d40b2b810 100644 --- a/src/main/java/org/matsim/run/batch/UtilsJR.java +++ b/src/main/java/org/matsim/run/batch/UtilsJR.java @@ -29,7 +29,7 @@ public class UtilsJR { static void produceDiseaseImportPlot(Map> infections_pers_per_day) { LocalDate startDate = LocalDate.of(2020, 2, 1); - LocalDate endDate = LocalDate.of(2023, 2, 25); + LocalDate endDate = LocalDate.of(2024, 6, 25); DateColumn recordsDate = DateColumn.create("date"); DoubleColumn values = DoubleColumn.create("import"); diff --git a/src/main/java/org/matsim/run/modules/CologneStrainScenario.java b/src/main/java/org/matsim/run/modules/CologneStrainScenario.java deleted file mode 100644 index ff2d9018d..000000000 --- a/src/main/java/org/matsim/run/modules/CologneStrainScenario.java +++ /dev/null @@ -1,27 +0,0 @@ -package org.matsim.run.modules; - -import org.matsim.episim.EpisimConfigGroup; -import org.matsim.episim.model.InfectionModelWithAntibodies; -import org.matsim.episim.model.vaccination.NoVaccination; -import org.matsim.scenarioCreation.RunTrial; - -/** - * Cologne scenario for calibration of different strains. - */ -public class CologneStrainScenario extends SnzCologneProductionScenario { - - - public CologneStrainScenario() { - super((Builder) new Builder() - .setScaleForActivityLevels(1.3) - .setSuscHouseholds_pct(0.0) - .setLeisureCorrection(RunTrial.parseParam("leisureCorrection", 1.9)) - .setVaccinations(Vaccinations.no) - .setActivityHandling(EpisimConfigGroup.ActivityHandling.startOfDay) - .setVaccinationModel(NoVaccination.class) - .setInfectionModel(InfectionModelWithAntibodies.class) - ); - - } - -} diff --git a/src/main/java/org/matsim/run/modules/SnzCologneProductionScenario.java b/src/main/java/org/matsim/run/modules/SnzCologneProductionScenario.java index df5a084eb..29dadac59 100644 --- a/src/main/java/org/matsim/run/modules/SnzCologneProductionScenario.java +++ b/src/main/java/org/matsim/run/modules/SnzCologneProductionScenario.java @@ -54,7 +54,7 @@ /** * Scenario for Cologne using Senozon events for different weekdays. */ - public class SnzCologneProductionScenario extends SnzProductionScenario { + public final class SnzCologneProductionScenario extends SnzProductionScenario { public static class Builder extends SnzProductionScenario.Builder { @@ -67,6 +67,8 @@ public static class Builder extends SnzProductionScenario.Builder inputDays, */ public static void configureWeather(EpisimConfigGroup episimConfig, WeatherModel weatherModel, File weather, File avgWeather, double maxOutdoorFraction) { if (weatherModel != WeatherModel.no) { - double midpoint1 = 0.1 * Double.parseDouble(weatherModel.toString().split("_")[1]); - double midpoint2 = 0.1 * Double.parseDouble(weatherModel.toString().split("_")[2]); + double midpoint1 = 0.1 * Double.parseDouble(weatherModel.toString().split("_")[1]); // 18.5 + double midpoint2 = 0.1 * Double.parseDouble(weatherModel.toString().split("_")[2]); // 25.0 try { Map outdoorFractions = EpisimUtils.getOutDoorFractionFromDateAndTemp2(weather, avgWeather, 0.5, midpoint1, midpoint2, midpoint1, midpoint1, 5., 1.0, maxOutdoorFraction); episimConfig.setLeisureOutdoorFraction(outdoorFractions); @@ -521,7 +521,7 @@ public enum Vaccinations {yes, no} public enum ChristmasModel {no, restrictive, permissive} - public enum WeatherModel {no, midpoints_175_175, midpoints_185_250, midpoints_175_250, midpoints_200_250, midpoints_175_200, midpoints_200_200} + public enum WeatherModel {no, midpoints_175_175, midpoints_185_250, midpoints_175_250, midpoints_200_250, midpoints_175_200, midpoints_200_200, midpoints_500_500} public enum AdjustRestrictions {yes, no} diff --git a/src/main/java/org/matsim/scenarioCreation/ConvertPersonAttributes.java b/src/main/java/org/matsim/scenarioCreation/ConvertPersonAttributes.java index 0330a86af..71bba8ad8 100644 --- a/src/main/java/org/matsim/scenarioCreation/ConvertPersonAttributes.java +++ b/src/main/java/org/matsim/scenarioCreation/ConvertPersonAttributes.java @@ -116,7 +116,7 @@ public Integer call() throws Exception { String attributesFileForConversion = input.toString(); if (personIds != null) { String outputPath = output.toAbsolutePath().toString(); - attributesFileForConversion = outputPath.substring(0, outputPath.lastIndexOf('\\')) + "/filtered_" + input.getFileName(); + attributesFileForConversion = outputPath.substring(0, outputPath.lastIndexOf('/')) + "/filtered_" + input.getFileName(); new ObjectAttributesXmlWriter(attributes).writeFile(attributesFileForConversion); } diff --git a/src/main/resources/postProcess.sh b/src/main/resources/postProcess.sh index 0451a5fc3..303cda7c5 100644 --- a/src/main/resources/postProcess.sh +++ b/src/main/resources/postProcess.sh @@ -27,7 +27,7 @@ else fi -classpath="matsim-episim-*-SNAPSHOT.jar" +classpath="matsim-episim-*.jar" # main main="analysis extractInfectionsByAge --population=$EPISIM_INPUT/ --district Berlin" diff --git a/src/main/resources/run.sh b/src/main/resources/run.sh index 1202f6a08..72a1a4366 100644 --- a/src/main/resources/run.sh +++ b/src/main/resources/run.sh @@ -11,7 +11,7 @@ date hostname -classpath="matsim-episim-*-SNAPSHOT.jar" +classpath="matsim-episim-*.jar" echo "***" echo "classpath: $classpath" diff --git a/src/main/resources/runParallel.sh b/src/main/resources/runParallel.sh index 1e3c9eb6f..b078a694c 100644 --- a/src/main/resources/runParallel.sh +++ b/src/main/resources/runParallel.sh @@ -8,7 +8,7 @@ hostname cd $SLURM_SUBMIT_DIR -classpath="matsim-episim-*-SNAPSHOT.jar" +classpath="matsim-episim-*.jar" echo "***" echo "classpath: $classpath" diff --git a/src/main/resources/runSlurm.sh b/src/main/resources/runSlurm.sh index c5ad98905..c68e3ea1d 100644 --- a/src/main/resources/runSlurm.sh +++ b/src/main/resources/runSlurm.sh @@ -16,7 +16,7 @@ hostname cd $SLURM_SUBMIT_DIR -classpath="matsim-episim-*-SNAPSHOT.jar" +classpath="matsim-episim-*.jar" echo "***" echo "classpath: $classpath" diff --git a/src/test/java/org/matsim/episim/model/DefaultAntibodyModelTest.java b/src/test/java/org/matsim/episim/model/DefaultAntibodyModelTest.java index 186883c51..45f0d8384 100644 --- a/src/test/java/org/matsim/episim/model/DefaultAntibodyModelTest.java +++ b/src/test/java/org/matsim/episim/model/DefaultAntibodyModelTest.java @@ -10,6 +10,7 @@ import org.apache.log4j.Logger; import org.assertj.core.data.Offset; import org.junit.*; +import org.matsim.core.utils.io.IOUtils; import org.matsim.episim.EpisimConfigGroup; import org.matsim.episim.EpisimPerson; import org.matsim.episim.EpisimTestUtils; @@ -29,6 +30,7 @@ import java.io.*; import java.nio.charset.StandardCharsets; +import java.nio.file.Path; import java.util.*; import static com.google.common.math.Quantiles.*; @@ -179,13 +181,16 @@ public void testNoImmunityEvents() { @Test public void testMixOfVaccinesAndInfections() { - List immunityEvents = List.of(VirusStrain.DELTA, VaccinationType.mRNA, VirusStrain.OMICRON_BA1, VirusStrain.OMICRON_BA2, VirusStrain.OMICRON_BA5,VirusStrain.OMICRON_BA5,VirusStrain.OMICRON_BA5); - IntList immunityEventDays = IntList.of(538, 644,720,736,845,958,979); + + List immunityEvents = List.of(VirusStrain.SARS_CoV_2, VaccinationType.mRNA, VirusStrain.DELTA); + IntList immunityEventDays = IntList.of(50, 200, 600); +// List immunityEvents = List.of(VirusStrain.DELTA, VaccinationType.mRNA, VirusStrain.OMICRON_BA1, VirusStrain.OMICRON_BA2, VirusStrain.OMICRON_BA5,VirusStrain.OMICRON_BA5,VirusStrain.OMICRON_BA5); +// IntList immunityEventDays = IntList.of(538, 644,720,736,845,958,979); // List immunityEvents = List.of(VaccinationType.mRNA); // IntList immunityEventDays = IntList.of(1); EpisimPerson person = EpisimTestUtils.createPerson(); - person.setImmuneResponseMultiplier(0.1); +// person.setImmuneResponseMultiplier(0.1); Int2ObjectMap> antibodyLevels = simulateAntibodyLevels(immunityEvents, immunityEventDays, 1000, person); @@ -210,6 +215,18 @@ public void testMixOfVaccinesAndInfections() { } producePlot(records, values, groupings, "nAb", "nAb: " + immunityEvents, "nAb.html"); + BufferedWriter writer = IOUtils.getBufferedWriter("nAb-table.csv"); + try { + writer.write("day,antibodies,scenario"); + for (int i = 0; i < records.size(); i++) { + writer.newLine(); + writer.write(records.get(i) + "," + values.get(i) + "," + groupings.get(i)); + } + } catch (IOException e) { + throw new RuntimeException(e); + } + + } // Plot 2: ve diff --git a/src/test/java/org/matsim/episim/model/vaccination/VaccinationFromDataTest.java b/src/test/java/org/matsim/episim/model/vaccination/VaccinationFromDataTest.java index 2dd2b986c..328369ce4 100644 --- a/src/test/java/org/matsim/episim/model/vaccination/VaccinationFromDataTest.java +++ b/src/test/java/org/matsim/episim/model/vaccination/VaccinationFromDataTest.java @@ -25,7 +25,7 @@ public class VaccinationFromDataTest { - private static final String URL = "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/master/Aktuell_Deutschland_Landkreise_COVID-19-Impfungen.csv"; + private static final String URL = "https://raw.githubusercontent.com/robert-koch-institut/COVID-19-Impfungen_in_Deutschland/main/Deutschland_Landkreise_COVID-19-Impfungen.csv"; private static File input;