Merge pull request #80 from matsim-org/develop

Develop

Makefile

@@ -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/Heinsberg.mk
+SUBDIRS := scenarios/Germany.mk
 define INCLUDE_FILE
 path = $S
 include $S

pom.xml

@@ -4,7 +4,7 @@
     <modelVersion>4.0.0</modelVersion>
     <groupId>com.github.matsim-org</groupId>
     <artifactId>matsim-episim</artifactId>
-    <version>21.4</version>
+    <version>21.5-SNAPSHOT</version>
 
     <name>MATSim Episim</name>
     <description>Epidemic simulation for MATSim</description>

scenarios/Dresden.mk

@@ -0,0 +1,54 @@
+
+in = $(WD)/snz/Dresden/original-data
+out = $(WD)/snz/Dresden/episim-input
+tmp = $(WD)/snz/Dresden/processed-data
+
+Dresden: $(out)/dresden_snz_episim_events_wt_100pt_split.xml.gz $(out)/dresden_snz_entirePopulation_emptyPlans_withDistricts_100pt_split.xml.gz
+	echo "Building Dresden scenario"
+
+$(tmp)/personIds.diluted.txt.gz:
+	$(sc) filterPersons $(in)/de2020gsmwt_events_reduced.xml.gz\
+	 --facilities $(in)/facilities_assigned_simplified.xml.gz\
+	 --shape-file $(out)/../shape-file/case-study_Dresden_PLZ.shp\
+	 --output $@
+
+$(out)/dresden_snz_entirePopulation_emptyPlans_100pt.xml.gz: $(tmp)/personIds.diluted.txt.gz
+	$(sc) convertPersonAttributes $(in)/populationAttributes.xml.gz\
+   --ids $<\
+   --requireAttribute "microm:modeled:age"\
+   --output $@
+
+$(out)/dresden_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz: $(out)/dresden_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
+
+###########
+# 100 pct
+###########
+# https://www.destatis.de/DE/Themen/Gesellschaft-Umwelt/Bevoelkerung/Haushalte-Familien/Tabellen/1-2-privathaushalte-bundeslaender.html
+
+$(out)/dresden_snz_entirePopulation_emptyPlans_withDistricts_100pt_filtered.xml.gz $(out)/dresden_snz_episim_events_wt_100pt.xml.gz &: \
+$(out)/dresden_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)/dresden_snz_entirePopulation_emptyPlans_withDistricts_100pt_filtered.xml.gz
+	mv $(tmp)/de2020gsmwt_events_reduced-1.0.xml.gz $(out)/dresden_snz_episim_events_wt_100pt.xml.gz
+	mv $(tmp)/de2020gsmsa_events_reduced-1.0.xml.gz $(out)/dresden_snz_episim_events_sa_100pt.xml.gz
+	mv $(tmp)/de2020gsmso_events_reduced-1.0.xml.gz $(out)/dresden_snz_episim_events_so_100pt.xml.gz
+
+$(out)/dresden_snz_entirePopulation_emptyPlans_withDistricts_100pt_split.xml.gz $(out)/dresden_snz_episim_events_wt_100pt_split.xml.gz &: \
+$(out)/dresden_snz_entirePopulation_emptyPlans_withDistricts_100pt_filtered.xml.gz $(out)/dresden_snz_episim_events_wt_100pt.xml.gz
+	$(sc) splitHomeFacilities $<\
+	 --events $(out)/dresden_snz_episim_events_wt_100pt.xml.gz\
+	 --events $(out)/dresden_snz_episim_events_sa_100pt.xml.gz\
+	 --events $(out)/dresden_snz_episim_events_so_100pt.xml.gz\
+	 --target="44.9 35.2 10.4 7.4 2.1"\
+	 --shape-file $(out)/../shape-file-utm32/case-study_Dresden_PLZ.shp\
+	 --output $(out)
+

scenarios/Germany.mk

@@ -0,0 +1,27 @@
+
+in = $(WD)/snz/Dresden/original-data
+out = $(WD)/snz/Dresden/episim-input
+tmp = $(WD)/snz/Dresden/processed-data
+
+Germany: $(out)/germany_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz
+	echo "Building Germany scenario"
+
+$(out)/germany_snz_entirePopulation_emptyPlans_100pt.xml.gz:
+	$(sc) convertPersonAttributes $(in)/populationAttributes.xml.gz\
+   --requireAttribute "microm:modeled:age"\
+   --output $@
+
+$(out)/germany_snz_entirePopulation_emptyPlans_withDistricts_100pt.xml.gz: $(out)/germany_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
+
+$(out)/germany_snz_entirePopulation_emptyPlans_withDistricts_25pt.xml.gz: $(out)/germany_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)/germany_snz_entirePopulation_emptyPlans_withDistricts_25pt.xml.gz

src/main/python/analysis/dispersion.R → src/main/R/dispersion.R

@@ -1,97 +1,99 @@
-
-library(fitdistrplus)
-library(dplyr)
-library(purrr)
-library(ggplot2)
-library(reticulate)
-
-use_miniconda(required = TRUE)
-
-np <- import("numpy")
-
-f <- "C:/home/Development/matsim-org/matsim-episim/output-berlin-25pct-superSpreader-calibrParam-1.65E-5/infectionEvents.txt"
-f <- "C:/home/Development/matsim-org/matsim-episim/output-berlin-25pct-superSpreader-calibrParam-2.7E-5/infectionEvents.txt"
-f <- "C:/home/Development/matsim-org/matsim-episim/output-base/infectionEvents.txt"
-
-counts <- function(f) {
-
-  data <- read.csv(f, sep = "\t")
-
-  df <- table(data[[2]])  
-
-  v <- as.numeric(df)
-
-  # Persons who did not infect other persons
-  no_inf <- setdiff(data$infected, data$infector)
-
-  res <- sort(c(rep(0, length(no_inf)), v))
-}
-
-est_disp <- function(f) {
-
-  matrix <- np$load(f)
-
-  est <- list()
-  inf80 <- list()
-
-  for(row in 1:nrow(matrix)) {
-
-      sprintf("Processing row %d\n", row)
-      v <- matrix[row,]
-      v <- v[!is.na(v)]
-
-      if (length(v) == 0) {
-          next
-      }
-
-      fit <- try(fitdist(v, fix.arg=list(mu=2.5), "nbinom"), silent = T)
-      if(inherits(fit, "try-error")) {
-        fit <- list(estimate=NaN)
-      }
-
-      est <- c(est, as.numeric(fit$estimate))
-
-      total = sum(v)
-      s80 <- v[0:-length(v) * 0.8]
-
-      inf80 <- c(inf80, sum(s80) * 100 / total)
-
-  }
-
-  df <- data.frame(as.numeric(est), as.numeric(inf80))
-  colnames(df) <- c("est", "top20")
-
-  return(df)
-}
-
-
-est_disp_zip <- function(f) {
-
-  tmpdir <- tempdir()
-  unzip(f, exdir = tmpdir )
-
-  res <- data.frame()
-
-  for (f in list.files(tmpdir, pattern = "*.npy", full.names = T)) {
-    df <- est_disp(f)
-    means <- colMeans(df, na.rm = T)
-    df <- data.frame(t(means))
-
-    row.names(df) <- c(basename(f)) 
-
-    res <- rbind(res, df)
-  }
-
-  return(res) 
-}
-
-setwd("C:/home/Development/matsim-org/matsim-episim/src/main/python/analysis")
-
-
-df <- est_disp_zip("data/infections.zip")
-
-# Testing the distribution
-df <- sort(rnbinom(15000, size=1, mu=2.5))
-plot(hist(df, breaks = 30))
-
-ggplot() + aes(df) + geom_histogram(binwidth=1, colour="black", fill="white", alpha=0.2)     
+
+library(fitdistrplus)
+library(dplyr)
+library(purrr)
+library(ggplot2)
+library(reticulate)
+
+use_miniconda(required = TRUE)
+
+np <- import("numpy")
+
+f <- "C:/home/Development/matsim-org/matsim-episim/output-berlin-25pct-superSpreader-calibrParam-1.65E-5/infectionEvents.txt"
+f <- "C:/home/Development/matsim-org/matsim-episim/output-berlin-25pct-superSpreader-calibrParam-2.7E-5/infectionEvents.txt"
+f <- "C:/home/Development/matsim-org/matsim-episim/output-base/infectionEvents.txt"
+
+counts <- function(f) {
+
+  data <- read.csv(f, sep = "\t")
+
+  df <- table(data[[2]])  
+
+  v <- as.numeric(df)
+
+  # Persons who did not infect other persons
+  no_inf <- setdiff(data$infected, data$infector)
+
+  res <- sort(c(rep(0, length(no_inf)), v))
+}
+
+est_disp <- function(f, mu=1) {
+
+  matrix <- np$load(f)
+
+  est <- list()
+  inf80 <- list()
+
+  # Loaded vector is one row
+  v <- as.numeric(matrix)
+  #for(row in 1:nrow(matrix)) {
+
+      #sprintf("Processing row %d\n", row)
+      #v <- matrix[row,]
+      v <- v[!is.na(v)]
+
+      if (length(v) == 0) {
+          next
+      }
+
+      fit <- try(fitdist(v, fix.arg=list(mu=mu), "nbinom"), silent = T)
+      if(inherits(fit, "try-error")) {
+        fit <- list(estimate=NaN)
+      }
+
+      est <- c(est, as.numeric(fit$estimate))
+
+      total = sum(v)
+      s80 <- v[0:-length(v) * 0.8]
+
+      inf80 <- c(inf80, sum(s80) * 100 / total)
+
+  #}
+
+  df <- data.frame(as.numeric(est), as.numeric(inf80))
+  colnames(df) <- c("est", "top20")
+
+  return(df)
+}
+
+
+est_disp_zip <- function(f) {
+
+  tmpdir <- tempdir()
+  unzip(f, exdir = tmpdir )
+
+  res <- data.frame()
+
+  for (f in list.files(tmpdir, pattern = "*.npy", full.names = T)) {
+    df <- est_disp(f)
+    means <- colMeans(df, na.rm = T)
+    df <- data.frame(t(means))
+
+    row.names(df) <- c(basename(f)) 
+
+    res <- rbind(res, df)
+  }
+
+  return(res) 
+}
+
+setwd("C:/Users/chris/Development/matsim-org/matsim-episim/biggest-jan")
+
+est_disp("0.46.npy")
+est_disp("0.56.npy")
+
+# Testing the distribution
+df <- sort(rnbinom(15000, size=1, mu=2.5))
+plot(hist(df, breaks = 30))
+
+ggplot() + aes(df) + geom_histogram(binwidth=1, colour="black", fill="white", alpha=0.2)