library(tidyverse)
library(readxl)
library(cowplot)
par(mar = c(4,4,1,1), oma = c(0,0,0,0), las = 1, xaxs = "i", yaxs = "i")

setwd("./___UITZONDERINGSGROND_6___ContactAnalysis")

source("./___UITZONDERINGSGROND_6___functions_new.r")
source("./___UITZONDERINGSGROND_6___modelContacts.r")

#list.files(path = "functions_new", full.names = TRUE) %>% walk(source)


# contact3 <- getContactData(dir = "data/Pienter3")
# population3 <- getPopulationData(dir = "data/", year = 2017)

# saveRDS(contact3, "../___UITZONDERINGSGROND_6___contact3.rds")
# saveRDS(population3, "../___UITZONDERINGSGROND_6___population3.rds")

contact3 <- readRDS("../___UITZONDERINGSGROND_6___contact3.rds")
population3 <- readRDS("../___UITZONDERINGSGROND_6___population3.rds")

contact3$cnt_location %>% table

contact3 %>% filter(part_age >= 18 & cnt_location == "Other place") %>% nrow
contact3 %>% filter(part_age >= 18 & !duplicated(part_id)) %>% nrow

1726/910

3863/910
2379/910

demo <- population3 %>% 
  mutate(age = cut(age, breaks = c(seq(0,80,10), Inf), right = FALSE, include.lowest = TRUE)) %>%
  group_by(age) %>%
  summarize(ntot = sum(population)) %>%
  ungroup %>%
  mutate(frac = ntot/sum(ntot)) %>%
  pull(frac)


sum(pull(filter(population3, age >= 10 & age < 18), population))/sum(pull(filter(population3, age >= 10 & age < 20), population))



locations <- levels(contact3$cnt_location)

contactpattern <- ContactPatternLocation(contact.data = contact3, pop.data = population3, countWhat = countContactsLoc)
#contactpattern2 <- ContactPatternLocation(contact.data = contact3, pop.data = population3, countWhat = countExposureContactsLoc)


plotContactMatrix(contactpattern, var.name = "m_smt", min = 0.1, max = 20)

# contactpattern$m_smt_inf <- as.vector(t(relinf*contactpattern2matrix(contactpattern.data = contactpattern, var.name = "m_smt")))
# plotContactMatrix(contactpattern.data = contactpattern, var.name = "m_smt_inf", min = 0.01, max = 10)


# relative infectiousness -> nog even netter!

# infectiousness ~ probability severe
# relinf <- c(0, 0.0133, 0.0141, 0.0173, 0.0321, 0.0916, 0.2361, 0.4563, 0.6909)
# 
# 10*c(0, 0.0133, 0.0141, 0.0173, 0.0321, 0.0916, 0.2361, 0.4563, 0.6909) + (1-c(0, 0.0133, 0.0141, 0.0173, 0.0321, 0.0916, 0.2361, 0.4563, 0.6909))
# 
# mild <- c(416, 540, 3553, 7431, 8215, 8818, 5915, 1498, 31)
# severe <- c(0, 7, 51, 131, 275, 916, 2026, 1788, 973)
# critical <- c(0, 2, 15, 38, 80, 274, 642, 682, 405)
# 
# ratio <- 10
# relinf <- (mild/(ratio) + (severe+critical))/(mild + severe + critical)

#relinf <- c(0.05, 0.05, 0.2, 0.5, 0.75, 0.9, 0.95, 0.95, 0.95) # guess symptomatic, comparable to symptomatic AR in China
relinf <- c(0.05, 0.1, 0.25, 0.55, 0.75, 0.9, 0.95, 0.95, 0.95) # guess symptomatic, comparable to symptomatic AR in China
obsChina <- c(416, 549, 3619, 7600, 8571, 10008, 8583, 3918, 1408)
popChina <- c(68313+63314, 60727+59251, 73185+100701, 88959+82553, 87713+105476, 96760+59823, 68044+51552, 32590+22553, 14708+6606+1964+455)
plot(obsChina/popChina, ylim = c(0, 0.08))
lines(0.075*relinf)



# obsSing <- c(3, 3, 17, 36, 26, 30, 28, 7, 0)
# popSing <- c(186+199, 207+227, 256+292, 281+304, 304 + 308, 309+304, 272+212, 136+93, 57+50)
# sum(obsSing)
# lines(obsSing/popSing, type = "l", col = 2)

all_measures <- read_excel("./___UITZONDERINGSGROND_6___MeasureInputs_21apr2020.xlsx")


# Pakketten voor ___UITZONDERINGSGROND_2___
measures____UITZONDERINGSGROND_2___ <- all_measures %>% filter(measureID %in% c(99, "D1", "D2", "D3", "D3praktijk", "D4", "E", "D1 + E", "D2 + E", "UK_SDO", "UK_SD", "UK_PC"))
measures____UITZONDERINGSGROND_2___ <- all_measures %>% filter(measureID %in% c(99, "D3praktijk_leisure80", "D3praktijk_leisure70"))
measures____UITZONDERINGSGROND_2___ <- all_measures %>% filter(measureID %in% c(99, "D3praktijk_leisure80", "D3praktijk_leisure70", "D3praktijk_leisure80_plusprimschool", "D3praktijk_leisure70_plusprimschool"))
measures____UITZONDERINGSGROND_2___ <- all_measures %>% filter(measureID %in% c(99, "P3like", "P3like_plusprimschool", "D3praktijk_leisure70"))
measures____UITZONDERINGSGROND_2___ <- all_measures %>% filter(measureID %in% c(99, "D3praktijk_leisure70", "D3praktijk_withprimschool", "D3praktijk_withprimschool_upto12", "D3praktijk_withprimschool_pessimistic", "D3praktijk_withprimschool_pessimistic_withyouths"))

eigenmaxinf <- max(eigen(relinf*contactpattern2matrix(contactpattern, var.name = "m_smt"))$value)
# from IC report 
symphos <- c(0.1,0.3,1.2,3.2,4.9,10.2,16.6,24.3,27.3)/100
hospcrit <- c(5,5,5,5,6.3,12.2,27.4,43.2,70.9)/100


res <- tibble(measureID = measures____UITZONDERINGSGROND_2___$measureID, effect = 0, check = 0, fsIC = 0)
for(i in 1:nrow(measures____UITZONDERINGSGROND_2___)) {
  tmp <-EffectOnTransmission(contactpattern, measures = measures____UITZONDERINGSGROND_2___[i,], relinf = rep(1, 9), locations = locations)
  res$effect[i] <- tmp
  tmp <- EffectOfMeasureOnContacts(contactpattern = contactpattern, par = measures____UITZONDERINGSGROND_2___[i,])
  tmp$n_cnt <- rowSums(tmp %>% select(paste0("n_cnt_", locations)))
  tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = FALSE)
  infvector <- as.numeric(measures____UITZONDERINGSGROND_2___[i, paste0("I_",1:9)])*relinf
  tmp <- infvector*contactpattern2matrix(tmp, var = "m_smt")
  res$check[i] <- 100*(max(eigen(tmp)$value)/eigenmaxinf - 1)
  fs <- final.size(R0 = 2.2, contact.mat = tmp, eigenmax = eigenmaxinf)
  res$fsIC[i] <- 17120000*sum(demo*fs*symphos*hospcrit)
}
res <- res %>% mutate(effectIC = 100*(fsIC/res$fsIC[1] - 1))




eigenmax <- max(eigen(contactpattern2matrix(contactpattern, var.name = "m_smt"))$value)
contactmatrix___UITZONDERINGSGROND_2___ <- list()
for(i in 1:nrow(measures____UITZONDERINGSGROND_2___)) {
  # adjust number of contacts per location per age group
  tmp <- EffectOfMeasureOnContacts(contactpattern = contactpattern, par = measures____UITZONDERINGSGROND_2___[i,])
  # count all contacts on different locations
  tmp$n_cnt <- rowSums(tmp %>% select(paste0("n_cnt_", locations)))
  # model new 'reported' contacts
  tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = FALSE)
  # get symmetric contact matrix
  tmp <- contactpattern2matrix(tmp, var.name = "c_smt")
  # make more symmetric :-)
  tmp <- (tmp + t(tmp))/2
  # scale with eigenmax of full contact matrix (note: no relinf, just the contacts)
  tmp <- tmp/eigenmax
  print(c(measures____UITZONDERINGSGROND_2___[i, "measureID"], max(eigen(demo*tmp)$value)))
  contactmatrix___UITZONDERINGSGROND_2___[[i]] <- tmp
}


contactmatrix___UITZONDERINGSGROND_2___[["D3praktijk_leisure70"]]

contactmatrix___UITZONDERINGSGROND_2___[["D3praktijk_withprimschool"]]
contactmatrix___UITZONDERINGSGROND_2___[["D3praktijk_withprimschool_upto12"]]
names(contactmatrix___UITZONDERINGSGROND_2___) <- measures_Don$measureID
# check
max(eigen(demo*contactmatrix___UITZONDERINGSGROND_2___[["D3praktijk_leisure70"]])$value)
max(eigen(demo*contactmatrix___UITZONDERINGSGROND_2___[["D3praktijk_leisure80"]])$value)
max(eigen(demo*contactmatrix___UITZONDERINGSGROND_2___[["D1 + E"]])$value)
max(eigen(demo*contactmatrix___UITZONDERINGSGROND_2___[["P3like"]])$value)
max(eigen(demo*contactmatrix___UITZONDERINGSGROND_2___[["P3like_plusprimschool"]])$value)


measures____UITZONDERINGSGROND_2___ %>% filter(measureID %in% c("D1", "D1 + E")) %>% View

saveRDS(contactmatrix___UITZONDERINGSGROND_2___, "ContactMatrixSchoolOpening_midpoint_21apr2020.rds")
1-0.4427217
mean(mev)
quantile(mev, probs = c(0.025,0.975))


# sensitivity matrices D3praktijk (22 mrt 2020)

nsens <- 200
rwork <- runif(nsens, min = 0.3, max = 0.6)
rschool <- runif(nsens, min = 0.7, max = 0.9)
#rleisure <- runif(nsens, min = 0.7, max = 0.9) # used for ContactmatricesD3praktijk_22mrt2020.rds (wrong) and ContactmatricesD3praktijk_leisure80_24mrt2020.rds
rleisure <- runif(nsens, min = 0.6, max = 0.8)  # used for ContactmatricesD3praktijk_23mrt2020.rds (wrong) and ContactmatricesD3praktijk_leisure70_24mrt2020.rds
rother <- runif(nsens, min = 0.3, max = 0.7)
rinf <- runif(nsens, min = 0.02, max = 0.08)

measures_sens <- tibble(measureID = 1:200)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*rwork)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - rwork)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - rschool)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.15)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - rleisure)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - rother)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("I_", i) := 1 - rinf)

# school opening
measures_sens$W_1 <- 1
measures_sens$S_1 <- 1

# measure sens per age category
nsens <- 1000
measures_sens <- tibble(measureID = 1:nsens)

# D3praktijk_withprimschool
r <- runif(nsens)
# home and work related
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.3 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.3 + 0.3*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.15)
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.6 + 0.2*r))
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.3 + 0.4*r))

# D3praktijk_withprimschool_upto12
r <- runif(nsens)
# home and work related
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.22 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.3 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.22 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.3 + 0.3*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.48 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.15)
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.6 + 0.2*r))
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.3 + 0.4*r))

# D3praktijk_withprimschool_pessimistic
r <- runif(nsens)
# home and work related
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.22 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.3 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.22 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.3 + 0.3*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.48 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.15)
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.4 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.6 + 0.2*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.26 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.3 + 0.4*r))

# D3praktijk_withprimschool_pessimistic_withyouths
r <- runif(nsens)
# home and work related
for(i in 1:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.22 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.3 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.22 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.3 + 0.3*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.48 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.15)
r <- runif(nsens)
for(i in 1:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (-0.1 + 0.2*r))
for(i in 3:3) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.4 + 0.2*r))
for(i in 4:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.6 + 0.2*r))
r <- runif(nsens)
for(i in 1:2) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (-0.1 + 0.2*r))
for(i in 3:3) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.26 + 0.2*r))
for(i in 4:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.3 + 0.4*r))

nsens <- 200
measures_sens <- tibble(measureID = 1:nsens)

# Batch0
r <- runif(nsens)
# home and work related
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.19 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.25 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.19 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.25 + 0.3*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.48 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.25)
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.3 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.45 + 0.2*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.26 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.3 + 0.4*r))

# Batch1
r <- runif(nsens)
# home and work related
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.15 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.20 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.15 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.20 + 0.3*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.48 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.25)
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.26 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.40 + 0.2*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.22 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.25 + 0.4*r))

# Batch2
r <- runif(nsens)
# home and work related
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.05 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.10 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.05 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.10 + 0.3*r))
r <- runif(nsens)
# for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
# for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.06 + 0.2*r))
# summer holiday
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.4 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.46 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.4)
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.10 + 0.2*r))
for(i in 3:7) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.30 + 0.2*r))
for(i in 8:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.40 + 0.2*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.10 + 0.2*r))
for(i in 3:7) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.20 + 0.4*r))
for(i in 8:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.25 + 0.4*r))

# Batch2+
r <- runif(nsens)
# home and work related
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(-0.05 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*(0.10 + 0.3*r))
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (-0.05 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - (0.10 + 0.3*r))
r <- runif(nsens)
# for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (-0.1 + 0.2*r))
# for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.06 + 0.2*r))
# summer holiday
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.4 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.46 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - (0.7 + 0.2*r))
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.4)
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.10 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - (0.30 + 0.2*r))
r <- runif(nsens)
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (-0.1 + 0.2*r))
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.0 + 0.2*r))
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - (0.0 + 0.4*r))

nsens <- 200
measures_sens <- tibble(measureID = 1:nsens)

# D3 as EpiPose1
r <- runif(nsens)
#for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 - 0.1 + 0.2*r) # fullhome
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 0.88 - 0.1 + 0.2*r) # fullhh
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.108 - 0.1 + 0.2*r)
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.012 - 0.1 + 0.2*r)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.0482)
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.00435 - 0.1 + 0.2*r)
r <- runif(nsens)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 0.386 - 0.1 + 0.2*r)



eigenmax <- max(eigen(contactpattern2matrix(contactpattern, var.name = "m_smt"))$value)
contactmatrix <- list()
mev <- rep(0, nrow(measures_sens))
for(i in 1:nrow(measures_sens)) {
  # adjust number of contacts per location per age group
  tmp <- EffectOfMeasureOnContacts(contactpattern = contactpattern, par = measures_sens[i,])
  # count all contacts on different locations
  tmp$n_cnt <- rowSums(tmp %>% select(paste0("n_cnt_", locations)))
  # model new 'reported' contacts
  tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = FALSE)
  # get symmetric contact matrix
  tmp <- contactpattern2matrix(tmp, var.name = "c_smt")
  # make more symmetric :-)
  tmp <- (tmp + t(tmp))/2
  # scale with eigenmax of full contact matrix (note: no relinf, just the contacts)
  tmp <- tmp/eigenmax
  mev[i] <- max(eigen(demo*tmp)$value)
  contactmatrix[[i]] <- tmp
}

mean(mev)
mean(mev_batch1)

mean(mev_batch2_sh)

mev_batch0 <- mev
mev_batch1 <- mev
mev_batch2 <- mev
mev_batch2_sh <- mev
mev_batch2plus <- mev
mev_batch2plus_sh <- mev
mevupto12 <-  mev
mev_EpiPose1_fullhome <- mev
mev_EpiPose1_fullhh <- mev



mean(mev1)
mean(mev2)
mean(mev3)

hist(mev3)
hist(mev1, add = TRUE, col = 2)
median(2.2*mev)

saveRDS(contactmatrix, "Contactmatrix_D3_asEpiPose1_fullhh_26mei2020.rds")

tmp <- c()
for(i in 1:1){
  tmp <- tmp + contactmatrix[[i]]
}

tmp/1000

contactmatrix[[1]] + contactmatrix[[2]]


#sens_original <- sens
sens <- measures_sens %>% 
  select(measureID, H_1, W_1, S_1, T_1, L_1, O_1) %>%
  bind_cols(mev = mev)

sens <- sens %>% pivot_longer(cols = c(H_1, W_1, S_1, T_1, L_1, O_1), names_to = "contact")

p <- ggplot(data = sens,
       mapping = aes(x = value, y = mev)) + 
  geom_point(pch = 21, bg = "white", col = "black") + 
  facet_wrap(facets = vars(contact),
             scales = "free_x") +
  labs(x = "multiplication of contact",
       y = "maximum eigenvalue") +
  theme_minimal()

ggsave(plot = p, "sens_percontacttype_D3praktijk_leisure70_24mrt2020.tiff", dpi = 150)


# which percentile gives mean(mev) with uncertainties?

rwork <- seq(0.3, 0.6, length.out = 11)
rschool <- seq(0.7, 0.9, length.out = 11)
rleisure <- seq(0.6, 0.8, length.out = 11)  # used for ContactmatricesD3praktijk_23mrt2020.rds
rother <- seq(0.3, 0.7, length.out = 11)
rinf <- seq(0.02, 0.08, length.out = 11) # not used

measures_sens <- tibble(measureID = 1:11)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 + 0.5*rwork)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 1 - rwork)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 1 - rschool)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.15)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 1 - rleisure)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - rother)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("I_", i) := 1 - rinf)


eigenmax <- max(eigen(contactpattern2matrix(contactpattern, var.name = "m_smt"))$value)
contactmatrix_D3praktijk <- list()
mev <- rep(0, nrow(measures_sens))
for(i in 1:nrow(measures_sens)) {
  # adjust number of contacts per location per age group
  tmp <- EffectOfMeasureOnContacts(contactpattern = contactpattern, par = measures_sens[i,])
  # count all contacts on different locations
  tmp$n_cnt <- rowSums(tmp %>% select(paste0("n_cnt_", locations)))
  # model new 'reported' contacts
  tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = FALSE)
  # get symmetric contact matrix
  tmp <- contactpattern2matrix(tmp, var.name = "c_smt")
  # make more symmetric :-)
  tmp <- (tmp + t(tmp))/2
  # scale with eigenmax of full contact matrix (note: no relinf, just the contacts)
  tmp <- tmp/eigenmax
  mev[i] <- max(eigen(demo*tmp)$value)
  contactmatrix_D3praktijk[[i]] <- tmp
}

cbind(rwork, mev)


bind_rows(measures____UITZONDERINGSGROND_2___ %>% filter(measureID == "D3praktijk_leisure70") %>% select(-measure, -measureID), measures_sens[6,]) %>% View


# get D3 praktijk matrix to compare to P3 (13 apr 2020)

measure <- all_measures %>% filter(measureID %in% c("P3like"))

# adjust number of contacts per location per age group
tmp <- EffectOfMeasureOnContacts(contactpattern = contactpattern, par = measure)
# count all contacts on different locations
tmp$n_cnt <- rowSums(tmp %>% select(paste0("n_cnt_", locations)))
# model new 'reported' contacts
tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = FALSE)

saveRDS(tmp, "/___UITZONDERINGSGROND_6___ContactmatrixP3like.rds")

plotContactMatrix(tmp, var.name = "m_smt")




tmp <- sapply(which(all_measures$measureID %in% c(28)), function(i) EffectOnTransmission(contactpattern, measures = all_measures[i,], relinf = relinf, locations = locations))
options(tibble.print_max = 50, tibble.print_min = 50)
tibble(measure = all_measures$measure[which(all_measures$measureID %in% c(28))], effect = tmp)

tmp <- sapply(1:nrow(all_measures), function(i) EffectOnTransmission(contactpattern, measures = all_measures[i,], relinf = relinf, locations = locations))
options(tibble.print_max = 50, tibble.print_min = 50)
tibble(measure = all_measures$measure, effect = tmp)[1:36,]


# tmp <- sapply(1:nrow(all_measures), function(i) EffectOnTransmission(contactpattern, measures = all_measures[i,], relinf = rep(1,9), locations = locations))
# options(tibble.print_max = 50, tibble.print_min = 50)
# tibble(measure = all_measures$measure, effect = tmp)[1:34,]

contactpattern2matrix(contactpattern.data = contactpattern2, var.name = "m_smt_inf_tot")


EffectOnTransmission(contactpattern2, measures = all_measures[c(23,9,10,11,12,14),])
EffectOnTransmission(contactpattern2, measures = all_measures[8:12,])


EffectOnTransmission(contactpattern = contactpattern2, effect_Home = 0)
EffectOnTransmission(contactpattern = contactpattern2, effect_Work = 0)
EffectOnTransmission(contactpattern = contactpattern2, effect_School = 0)
EffectOnTransmission(contactpattern = contactpattern2, effect_Transport = 0)
EffectOnTransmission(contactpattern = contactpattern2, effect_Leisure = 0)
EffectOnTransmission(contactpattern = contactpattern2, effect_Other = 0)


EffectOnTransmission(contactpattern = contactpattern2, 
                     effect_Home = 1.2, min_Home = 1, max_Home = 5, 
                     effect_School = 0.2, min_School = 2, max_School = 2,
                     effect_Leisure = 1.2, min_Leisure = 1, max_Leisure = 5)

EffectOnTransmission(contactpattern = contactpattern2, 
                     effect_Home = 1.2, min_Home = 3, max_Home = 6, 
                     effect_Work = 0.8, min_Work = 3, max_Work = 6,
                     effect_Leisure = 1.2, min_Leisure = 3, max_Leisure = 6)




EffectOnTransmission <- function(contactpattern, 
                                 effect_Home = 1, min_Home = 1, max_Home = 9,
                                 effect_Work = 1, min_Work = 1, max_Work = 9,
                                 effect_School = 1, min_School = 1, max_School = 9,
                                 effect_Transport = 1, min_Transport = 1, max_Transport = 9, 
                                 effect_Leisure = 1, min_Leisure = 1, max_Leisure = 9,
                                 effect_Other = 1, min_Other = 1, max_Other = 9) {
  
  eigenmax <- max(eigen(contactpattern2matrix(contactpattern.data = contactpattern, var.name = "m_smt_inf_tot"))$value)
  
  tmp <- contactpattern %>% mutate(m_smt_inf_Home = ifelse(as.integer(part_age) >= min_Home & as.integer(part_age) <= max_Home, effect_Home*m_smt_inf_Home, m_smt_inf_Home),
                                   m_smt_inf_Home = ifelse(as.integer(cnt_age) >= min_Home & as.integer(cnt_age) <= max_Home, effect_Home*m_smt_inf_Home, m_smt_inf_Home),
                                   m_smt_inf_Work = ifelse(as.integer(part_age) >= min_Work & as.integer(part_age) <= max_Work, effect_Work*m_smt_inf_Work, m_smt_inf_Work),
                                   m_smt_inf_Work = ifelse(as.integer(cnt_age) >= min_Work & as.integer(cnt_age) <= max_Work, effect_Work*m_smt_inf_Work, m_smt_inf_Work),
                                   m_smt_inf_School = ifelse(as.integer(part_age) >= min_School & as.integer(part_age) <= max_School, effect_School*m_smt_inf_School, m_smt_inf_School),
                                   m_smt_inf_School = ifelse(as.integer(cnt_age) >= min_School & as.integer(cnt_age) <= max_School, effect_School*m_smt_inf_School, m_smt_inf_School),
                                   m_smt_inf_Transport = ifelse(as.integer(part_age) >= min_Transport & as.integer(part_age) <= max_Transport, effect_Transport*m_smt_inf_Transport, m_smt_inf_Transport),
                                   m_smt_inf_Transport = ifelse(as.integer(cnt_age) >= min_Transport & as.integer(cnt_age) <= max_Transport, effect_Transport*m_smt_inf_Transport, m_smt_inf_Transport),
                                   m_smt_inf_Leisure = ifelse(as.integer(part_age) >= min_Leisure & as.integer(part_age) <= max_Leisure, effect_Leisure*m_smt_inf_Leisure, m_smt_inf_Leisure),
                                   m_smt_inf_Leisure = ifelse(as.integer(cnt_age) >= min_Leisure & as.integer(cnt_age) <= max_Leisure, effect_Leisure*m_smt_inf_Leisure, m_smt_inf_Leisure),
                                   `m_smt_inf_Other place` = ifelse(as.integer(part_age) >= min_Other & as.integer(part_age) <= max_Other, effect_Other*`m_smt_inf_Other place`, `m_smt_inf_Other place`),
                                   `m_smt_inf_Other place` = ifelse(as.integer(cnt_age) >= min_Other & as.integer(cnt_age) <= max_Other, effect_Other*`m_smt_inf_Other place`, `m_smt_inf_Other place`))
  
  tmp$m_smt_inf_tot <- rowSums(tmp %>% select(paste0("m_smt_inf_", levels(contact3$cnt_location))))
  neweigenmax <- max(eigen(contactpattern2matrix(tmp, var.name = "m_smt_inf_tot"))$value)
  
  return(100*(neweigenmax - eigenmax)/eigenmax)
}




tmp$m_smt_inf_School 
contactpattern$m_smt_inf_School

contactpattern$part_age %>% as.integer()

# TODO
# 2 plots van 6 locaties, white label, infector age, [variable legend] number of (infectious) contact(hours) per infector
# percentage transmissiereductie per locatie (in labels?)
# percentage transmissiereductie per age class?
# age distribution met vergelijking China obs

contactpattern2$m

tmp <- sapply(1:6, function(i) max(eigen(contactpattern2matrix(contactpattern2[[i]], var.name = "m_smt_inf"))$value))
tmp/sum(tmp)

varname <- "m_smt_inf"
maxeigen <- max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname))$value)
tmp[1] <- maxeigen - max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname)-contactpattern2matrix(contactpattern2[[1]], var.name = varname))$value)
tmp[2] <- maxeigen - max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname)-contactpattern2matrix(contactpattern2[[2]], var.name = varname))$value)
tmp[3] <- maxeigen - max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname)-contactpattern2matrix(contactpattern2[[3]], var.name = varname))$value)
tmp[4] <- maxeigen - max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname)-contactpattern2matrix(contactpattern2[[4]], var.name = varname))$value)
tmp[5] <- maxeigen - max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname)-contactpattern2matrix(contactpattern2[[5]], var.name = varname))$value)
tmp[6] <- maxeigen - max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname)-contactpattern2matrix(contactpattern2[[6]], var.name = varname))$value)
tmp/maxeigen

maxeigen <- max(eigen(contactpattern2matrix(contactpattern2[[7]], var.name = varname))$value)
partR0 <- rep(0, 9)
for(i in 1:9) {
  tmp <- contactpattern2matrix(contactpattern2[[7]], var.name = varname) 
  tmp[i,] <- 0
  partR0[i] <- maxeigen-max(eigen(tmp)$values)
}
partR0/maxeigen
sum(partR0)/maxeigen
partR0/sum(partR0)

contactpattern2matrix(contactpattern2[[5]], var.name = "m_smt")

library(cowplot)


for(i in 1:6) {
  loc <- locations[i]
  tmp <- contactpattern
  tmp$n_cnt <- tmp[[paste0("n_cnt_", loc)]]
  tmp <- modelContacts(contactpattern.data = tmp)
  contactpattern[[paste0("m_smt_", loc)]] <- tmp$m_smt
  contactpattern[[paste0("m_smt_inf_", loc)]] <- as.vector(t(relinf*contactpattern2matrix(contactpattern.data = tmp, var.name = paste0("m_smt_",loc))))
  
}


max(contactpattern %>% select(paste0("m_smt_", locations)))
min(contactpattern %>% select(paste0("m_smt_", locations)))

max(contactpattern %>% select(paste0("m_smt_inf_", locations)))
min(contactpattern %>% select(paste0("m_smt_inf_", locations)))




plots <- list()
for(i in 1:6) {
  loc <- locations[i]
  plots[[i]] <- plotContactMatrix(contactpattern, var.name = paste0("m_smt_inf_", loc), min = 0.0001, max = 2) + 
    geom_label(data = data.frame(x = 5, y = 8, label = loc),
               mapping = aes(x = x, y = y, label = label),
               inherit.aes = FALSE)
}

pdf(file = "locations_inf.pdf", height = 8, width = 12, pointsize = 20)
plot_grid(plotlist = plots, ncol = 3)
dev.off()

pdf(file = "Mmatrix_perlocation_inf2.pdf", height = 18, width = 12, pointsize = 12)
plot_grid(plotlist = plots_inf, ncol = 2)
dev.off()

max(contactpattern2[[3]]$m_smt)

A <- relinf*contactpattern2matrix(contactpattern, var.name = "m_smt")
#A[,5] <- 0
v <- eigen(A)$vectors[,1]
w <- (eigen(A)$vectors %>% solve)[1,]
plot(abs(v), type = "l", ylim = c(0,1))
lines(abs(w), col = 2)
lines(relinf*abs(w), col = 4)
obs <- c(416, 549, 3619, 7600, 8571, 10008, 8583, 3918, 1408)
points(2*obs/sum(obs))

ggplot() + 
  geom_bar(data = tibble(x = 1:9, y = abs(w)),
           mapping = aes(x=x, y = y),
           stat = "identity",
           fill = "Grey") +
  scale_x_continuous(
    breaks = 1:9,
    labels = levels(contactpattern$part_age)) +
  scale_y_continuous(
    expand = expansion(c(0, 0.05))
  ) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.title = element_blank(),
        axis.text.y = element_blank(),
        axis.line.y = element_blank(),
        axis.ticks.y = element_blank())
# + 
#   geom_line(data = tibble(x = 1:9, y = 0.075*relinf),
#             mapping = aes(x=x, y = y),
#             lwd = 2,
#             col = "Darkblue")

plotContactMatrix(contactpattern, var.name = "m_smt_inf")


(416+549)/sum(obs)
(3619 + 7600 + 8571)/sum(obs)
(10008 + 8583)/sum(obs)
(3918 + 1408)/sum(obs)


plot(obs/sum(obs), ylim = c(0,0.25))
lines(1.8*relinf*demo)

popChina <- c(68313+63314, 60727+59251, 73185+100701, 88959+82553, 87713+105476, 96760+59823, 68044+51552, 32590+22553, 14708+6606+1964+455)
plot(obs/popChina, ylim = c(0, 0.08))
lines(0.075*relinf)

v/sum(v)
relinf/sum(relinf)




contact.mat <- contactpattern2matrix(contactpattern[[7]], var.name = "m_smt_inf")
R0 <- 1.2

plot(final.size(R0 = 1.8, contact.mat = contactpattern2matrix(contactpattern, var.name = "m_smt")), type = "l", ylim = c(0,1))
lines(final.size(R0 = 1.8, contact.mat = contactpattern2matrix(contactpattern, var.name = "m_smt")), col = 2)

plot(demo*final.size(R0 = 1.8, contact.mat = contactpattern2matrix(contactpattern, var.name = "m_smt")), col = 1, type = "l", ylim = c(0, 0.13))
lines(relinf*demo*final.size(R0 = 1.8, contact.mat = contactpattern2matrix(contactpattern, var.name = "m_smt")), col = 4)
lines(relinf*demo*final.size(R0 = 1.8, contact.mat = contactpattern2matrix(contactpattern2[[7]], var.name = "m_smt_inf")), col = 4)


symphos <- c(0.1,0.3,1.2,3.2,4.9,10.2,16.6,24.3,27.3)/100
hospcrit <- c(5,5,5,5,6.3,12.2,27.4,43.2,70.9)/100
IFR <- c(0.002,0.006,0.03,0.08,0.15,0.6,2.2,5.1,9.3)/100
sum(demo*IFR)
0.044/sum(demo*symphos)
sum(demo*symphos*hospcrit)
2/3

plot(symphos*hospcrit)


IFR/(0.5*symphos*hospcrit)

sum(demo*relinf)

DPpop <- c(16,23,347,429,333,398,924,1015,215+11)
DPpos <- c(1,2,20,23,25,49,129,228,52+2)
plot(DPpos/DPpop)
sum(DPpos)/sum(DPpop)
relsus <- DPpos/DPpop

# relsus Hiroshi
relsus <- c(0.05, 0.08, 0.1, 0.1, 0.12, 0.22, 0.19, 0.13, 0.08)