Scenario hub round3

src/main/R/Covid19ScenarioHubSubmissionRound3.R

@@ -0,0 +1,114 @@
+
+library(lubridate)
+library(tidyverse)
+library(readr)
+
+
+
+rm(list=ls())
+source("/Users/jakob/git/matsim-episim/src/main/R/masterJR-utils.R", encoding = 'utf-8')
+
+# Global variables
+directory <- "/Users/jakob/git/public-svn/matsim/scenarios/countries/de/episim/battery/jakob/2022-09-15/1-eu-b-noAgg/"
+
+origin_date <- ymd("2022-07-24") # first day (sunday) of epiweek 30, 2022
+end_date <- ymd("2023-07-29") # last day (saturday) of epiweek 30, 2023 (latest possible value)
+
+
+
+## read & prep infections
+infections_raw <- read_combine_episim_output_zipped(directory,"infections.txt.csv")
+
+infections_incidence <- convert_infections_into_incidence(directory,infections_raw,FALSE) %>%
+  select(-c(infections_week,nShowingSymptomsCumulative, district, incidence))
+
+
+population_cologne <- infections_raw[1, "nSusceptible"]
+
+population_germany <- 83695430 #https://www.destatis.de/EN/Themes/Society-Environment/Population/Current-Population/Tables/liste-current-population.html
+scale_factor <- population_germany / population_cologne # pop germany / pop koelln
+
+infections_ready <- infections_incidence %>%
+  filter(date >= origin_date & date <= end_date) %>%
+  filter(campDuration == 91) %>%
+  filter(vacPool == "vaccinated") %>%
+  filter(vacFreq!="withStrain") %>%
+  mutate(optimistic = (mutEsc==3.7)) %>%
+  # mutate(weekday = lubridate::wday(date, label = TRUE)) %>%
+  mutate(year = epiyear(date)) %>%
+  mutate(epiweek = epiweek(date)) %>%
+  group_by(seed,vacCamp,vacType,year,epiweek) %>%
+  summarise(value = sum(infections) * scale_factor, target_end_date = last(date) ) %>%
+  mutate(target_variable = "inc infection") %>%
+  select(year, epiweek, target_end_date, target_variable, value, seed, vacCamp, vacType)
+
+
+## read & prep hospitalizations
+hosp_raw <- read_combine_episim_output_zipped(directory,"post.hospital.tsv")
+
+hosp_ready <- hosp_raw %>%
+  filter(date >= origin_date & date <= end_date) %>%
+  mutate(year = year(date)) %>%
+  mutate(wkday = lubridate::wday(date, label = TRUE)) %>%
+  filter(wkday == "Sat") %>%
+  filter(measurement == "intakesHosp") %>%
+  filter(severity == "Omicron") %>%
+  mutate(epiweek = epiweek(date)) %>%
+  rename("target_end_date" = date, value = n) %>%
+  mutate(value = value * population_cologne / 100000 * scale_factor) %>%
+  mutate(target_variable = "inc hosp") %>%
+  select(year, epiweek, target_end_date, target_variable, value, seed, vacCamp, vacType)
+
+
+
+
+# combine two dataframes and modify columns to match specs
+combined <- rbind(infections_ready, hosp_ready)
+# combined <- infections_ready #todo revert
+
+
+
+
+seed <- unique(combined$seed)
+sample <- seq(length(seed))
+map <- data.frame(seed,sample)
+
+final <- combined %>% filter(vacCamp!="off") %>%
+  mutate(scenario_id = paste0(vacCamp,"_",vacType)) %>%
+  mutate(scenario_id = case_when(scenario_id == "60plus_omicronUpdate"~"A-2022-07-24",
+                              scenario_id == "18plus_omicronUpdate"~"B-2022-07-24",
+                              scenario_id == "60plus_mRNA"~"C-2022-07-24",
+                              scenario_id == "18plus_mRNA"~"D-2022-07-24")) %>%
+  merge(map,by ="seed") %>%
+  mutate(horizon =  case_when(year == 2022 ~ epiweek - 29, year == 2023~ (52-29) + epiweek)) %>%
+  mutate("origin_date" = "2022-07-24") %>%
+  mutate("location" = "DE") %>%
+  mutate(value = round(value)) %>%
+  select(origin_date,scenario_id, horizon, target_end_date, location, sample,target_variable, value) %>%
+  arrange(scenario_id,sample,horizon) %>%
+  mutate(horizon = paste0(horizon," wk"))
+
+write.csv(final,"/Users/jakob/git/covid19-scenario-hub-europe/data-processed/MODUS_Covid-Episim/2022-07-24-MODUS_Covid-Episim.csv", row.names = FALSE)
+
+
+vac <- read_combine_episim_output_zipped(directory,"vaccinationsPerAgeGroup.tsv")
+
+vac_filtered <- vac %>% filter(seed == 4711) %>%
+  filter(vacFreq=="annual") %>%
+  filter(vacCamp == "18plus") %>%
+  filter(vacPool=="vaccinated") %>%
+  filter(campDuration==91.) %>% 
+  filter(mutEsc==3.7)
+
+
+# xxx <- read.delim("/Users/jakob/antibodies_2022-07-23.tsv",sep = "\t")
+#
+# yyy <- xxx %>% filter(nVaccinations == 0 & nInfections == 0)
+#
+# nrow(yyy)/nrow(xxx)
+
+
+
+
+
+

src/main/R/extendingDiseaseImport.R

@@ -0,0 +1,57 @@
+library(lubridate)
+library(ggplot2)
+raw <- read.csv("~/git/shared-svn/projects/episim/matsim-files/snz/Cologne/episim-input/cologneDiseaseImport_Projected.csv")
+
+x <- raw
+
+# first date after previous projections. 
+date <-dmy(tail(raw$date,1)) +1
+
+# we iterate day for day until the end of 2032. Each day, we check whether to 
+# add a import wave from 2022 or 2023 into the current year. 
+while(date<ymd("2032-12-31")){
+
+  date <- date+1
+
+  # default value for cases
+  cases <- 1
+
+  # if the date is before Jan 25th, we copy 2023 import into the current year 
+  if(month(date) < 2 && mday(date)<=24){
+    cases <- x[dmy(x$date)==ymd(paste0("2023","-",month(date),"-",mday(date))),]$cases
+  }
+  # if the date is after July 4th, we copy 2022's import into the current year 
+  if((month(date) == 7 && mday(date) >4) || month(date)>7){
+    cases <- x[dmy(x$date)==ymd(paste0("2022","-",month(date),"-",mday(date))),]$cases
+  }
+
+
+  x[nrow(x) + 1,] = c(format(date,"%d.%m.%y"), cases)
+}
+
+
+
+ggplot(x) + geom_line(aes(dmy(date),as.numeric(cases)))  
+
+
+# Here we add the hypothetical easter disease import
+# This always begins on april 1 of the year x and continues the same number of 
+# days as the autumn import (absolute values are also same, copied to array below)
+easter_import <-c(3,11,20,27,34,40,45,50,54,57,59,61,62,63,62,61,59,57,54,50,45,40,34,27,20,11,3)
+
+for(year in 2023:2032){
+
+  date <- ymd(paste0(year,"-04-01"))
+  for(import_for_day in easter_import){
+    x[which(dmy(x$date)==date),"cases"] <- import_for_day
+    date <- date + 1
+  }
+
+}
+
+
+ggplot(x) + geom_line(aes(dmy(date),as.numeric(cases)))  
+
+
+write.csv(x,file='~/git/shared-svn/projects/episim/matsim-files/snz/Cologne/episim-input/cologneDiseaseImport_Projected_2032.csv', row.names=FALSE, quote = FALSE)
+

src/main/java/org/matsim/episim/EpisimPerson.java

@@ -549,7 +549,7 @@ public void setVaccinationStatus(VaccinationStatus vaccinationStatus, Vaccinatio
 		vaccinations.add(type);
 		vaccinationDates.add(iteration);
 
-		reporting.reportVaccination(personId, iteration, type, vaccinations.size());
+		reporting.reportVaccination(personId, iteration, type, vaccinations.size(), age);
 	}
 
 	public TestStatus getTestStatus() {

src/main/java/org/matsim/episim/EpisimReporting.java

@@ -23,10 +23,9 @@
 import com.google.common.collect.Lists;
 import com.google.inject.Inject;
 import com.typesafe.config.ConfigRenderOptions;
-import it.unimi.dsi.fastutil.objects.Object2DoubleMap;
-import it.unimi.dsi.fastutil.objects.Object2IntMap;
-import it.unimi.dsi.fastutil.objects.Object2IntOpenHashMap;
-import it.unimi.dsi.fastutil.objects.ObjectIntPair;
+import it.unimi.dsi.fastutil.ints.Int2ObjectAVLTreeMap;
+import it.unimi.dsi.fastutil.ints.Int2ObjectMap;
+import it.unimi.dsi.fastutil.objects.*;
 import org.apache.commons.compress.archivers.tar.TarArchiveEntry;
 import org.apache.commons.compress.archivers.tar.TarArchiveOutputStream;
 import org.apache.logging.log4j.LogManager;
@@ -64,6 +63,23 @@
  */
 public final class EpisimReporting implements BasicEventHandler, Closeable, Externalizable {
 
+	/**
+	 * Age groups used for various outputs. AgeGroup -> minimum age of age group.
+	 * Important: age groups must be in descending order
+	 */
+	public enum AgeGroup {
+		age_60_plus(60),
+		age_18_59(18),
+		age_12_17(12),
+		age_0_11(0);
+
+		public final int lowerBoundAge;
+
+		AgeGroup(int lowerBoundAge) {
+			this.lowerBoundAge = lowerBoundAge;
+		}
+	}
+
 	private static final Logger log = LogManager.getLogger(EpisimReporting.class);
 	private static final AtomicInteger specificInfectionsCnt = new AtomicInteger(300);
 
@@ -103,6 +119,10 @@ public final class EpisimReporting implements BasicEventHandler, Closeable, Exte
 	 * Map of (VaccinationType, nth Vaccination) -> Number per day
 	 */
 	public final Object2IntMap<ObjectIntPair<VaccinationType>> vaccinationStats = new Object2IntOpenHashMap<>();
+	/**
+	 * Map of Age group -> Number of Vaccinations administered per day
+	 */
+	public final Object2IntMap<AgeGroup> vaccinationPerAgeGroupStats = new Object2IntOpenHashMap<>();
 
 	/**
 	 * Number format for logging output. Not static because not thread-safe.
@@ -146,6 +166,8 @@ public final class EpisimReporting implements BasicEventHandler, Closeable, Exte
 	private BufferedWriter vaccinationsPerType;
 	private BufferedWriter vaccinationsPerTypeAndNumber;
 
+	private BufferedWriter vaccinationsPerAgeGroup;
+
 	private final Map<String, BufferedWriter> externalWriters = new HashMap<>();
 
 	private String memorizedDate = null;
@@ -209,6 +231,7 @@ public final class EpisimReporting implements BasicEventHandler, Closeable, Exte
 		antibodiesPerPerson = EpisimWriter.prepare(base + "antibodies.tsv", "day", "date", (Object[]) VirusStrain.values());
 		vaccinationsPerType = EpisimWriter.prepare(base + "vaccinations.tsv", "day", "date", (Object[]) VaccinationType.values());
 		vaccinationsPerTypeAndNumber = EpisimWriter.prepare(base + "vaccinationsDetailed.tsv", "day", "date", "type", "number", "amount");
+		vaccinationsPerAgeGroup = EpisimWriter.prepare(base + "vaccinationsPerAgeGroup.tsv", "day", "date", "age_group","amount");
 
 		sampleSize = episimConfig.getSampleSize();
 		writeEvents = episimConfig.getWriteEvents();
@@ -281,6 +304,7 @@ void append(String date) throws IOException {
 		antibodiesPerPerson = EpisimWriter.prepare(base + "antibodies.tsv");
 		vaccinationsPerType = EpisimWriter.prepare(base + "vaccinations.tsv");
 		vaccinationsPerTypeAndNumber = EpisimWriter.prepare(base + "vaccinationsDetailed.tsv");
+		vaccinationsPerAgeGroup = EpisimWriter.prepare(base + "vaccinationsPerAgeGroup.tsv");
 		// cpu time is overwritten
 		cpuTime = EpisimWriter.prepare(base + "cputime.tsv", "iteration", "where", "what", "when", "thread");
 		memorizedDate = date;
@@ -585,9 +609,22 @@ void reporting(Map<String, InfectionReport> reports, int iteration, String date)
 
 			writer.append(vaccinationsPerTypeAndNumber, vacOut);
 		}
-
 		vaccinationStats.clear();
 
+		// vaccinations per age group
+		String[] vacPerAgeGroupOut = new String[4];
+		vacPerAgeGroupOut[0] = String.valueOf(iteration);
+		vacPerAgeGroupOut[1] = date;
+		for (AgeGroup ageGroup : AgeGroup.values()) {
+			vacPerAgeGroupOut[2] = ageGroup.name();
+			vacPerAgeGroupOut[3] = Integer.toString(vaccinationPerAgeGroupStats.getInt(ageGroup));
+
+			writer.append(vaccinationsPerAgeGroup, vacPerAgeGroupOut);
+
+		}
+
+		vaccinationPerAgeGroupStats.clear();
+
 		// Write all reports for each district
 		for (InfectionReport r : reports.values()) {
 			if (r.name.equals("total")) continue;
@@ -771,11 +808,26 @@ void reportPersonStatus(EpisimPerson person, EpisimPersonStatusEvent event) {
 	/**
 	 * Report the vaccination of a person.
 	 */
-	void reportVaccination(Id<Person> personId, int iteration, VaccinationType type, int n) {
+	void reportVaccination(Id<Person> personId, int iteration, VaccinationType type, int n, int age) {
 
 		vaccinations.merge(type, 1, Integer::sum);
 		vaccinationStats.merge(ObjectIntPair.of(type, n), 1, Integer::sum);
 
+		int previousMinAge = 1000;
+
+		for (AgeGroup ageGroup : AgeGroup.values()) {
+
+			int minAge = ageGroup.lowerBoundAge;
+			if (minAge >= previousMinAge) {
+				throw new RuntimeException("AGE_GROUP Map not sorted in descending order");
+			} else if (age >= minAge) {
+				vaccinationPerAgeGroupStats.merge(ageGroup, 1, Integer::sum);
+				break;
+			}
+
+			previousMinAge = minAge;
+		}
+
 		manager.processEvent(new EpisimVaccinationEvent(EpisimUtils.getCorrectedTime(episimConfig.getStartOffset(), 0, iteration), personId, type, n));
 	}
 
@@ -879,6 +931,7 @@ public void close() {
 		writer.close(cpuTime);
 		writer.close(vaccinationsPerType);
 		writer.close(vaccinationsPerTypeAndNumber);
+		writer.close(vaccinationsPerAgeGroup);
 
 		for (BufferedWriter v : externalWriters.values()) {
 			writer.close(v);

src/main/java/org/matsim/episim/EpisimUtils.java

@@ -565,14 +565,14 @@ public static Map<LocalDate, Double> getOutDoorFractionFromDateAndTemp2(File wea
 
 			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 {
 				outdoorFractions.put(date, maxOutdoorFraction * getOutDoorFractionFromDateAndTemp(date, TmidSpring, TmidFall, Trange, tMax, prcp, rainThreshold, alpha));
 			}
-			
+
 			lastDate = date;
 		}
 
@@ -589,7 +589,8 @@ public static Map<LocalDate, Double> getOutDoorFractionFromDateAndTemp2(File wea
 			prcpPerDay.put(monthDay, prcp);
 		}
 
-		for (int i = 1; i < 365*3; i++) {
+		// project weather 11 years into the future
+		for (int i = 1; i < 365*11; i++) {
 			LocalDate date = lastDate.plusDays(i);
 			int month = date.getMonth().getValue();
 			int day = date.getDayOfMonth();

src/main/java/org/matsim/episim/analysis/ExtractInfectionsByAge.java

@@ -53,7 +53,7 @@ public class ExtractInfectionsByAge implements Callable<Integer> {
 	private boolean detailed = false;
 
 	@CommandLine.Option(names = "--age-groups", description = "Age groups as list of bin edges")
-	private List<Integer> ageGroups = List.of(0, 5, 10, 15, 20, 25, 30, 40, 50, 60, 70, 80, 90);
+	private List<Integer> ageGroups = List.of(0, 5, 10, 15, 18, 20, 25, 30, 40, 50, 60, 70, 80, 90);
 
 	@CommandLine.Option(names = "--attr", defaultValue = "microm:modeled:age", description = "Name of the age attribute")
 	private String ageAttr;