library(tidyverse)
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")

list.files(path = "functions", 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")

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)

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(contactpattern2, var.name = "m_smt_tot", 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)

# contactpattern <- InfectiousContactPattern(contactpattern, locations = locations, relinf = relinf)
# contactpattern2 <- InfectiousContactPattern(contactpattern2, locations = locations, relinf = relinf)


all_measures <- read_tsv("./___UITZONDERINGSGROND_6___MeasureInputs_12mrt2020.txt")


EffectOnTransmission(contactpattern, measures = tibble(effect_Home = 0), relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Work = 0), relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_School = 0), relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Transport = 0), relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Leisure = 0), relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Other = 0), relinf = relinf, locations = locations)

EffectOnTransmission(contactpattern, measures = tibble(effect_Home = 0), relinf = rep(1, 9), locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Work = 0), relinf = rep(1, 9), locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_School = 0), relinf = rep(1, 9), locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Transport = 0), relinf = rep(1, 9), locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Leisure = 0), relinf = rep(1, 9), locations = locations)
EffectOnTransmission(contactpattern, measures = tibble(effect_Other = 0), relinf = rep(1, 9), locations = locations)


EffectOnTransmission(contactpattern2, measures = all_measures[c(6, 23),])

all_measures$effect_Home <- 1
all_measures$effect_Leisure <- 1
all_measures$effect_Work <- 1

# Pakketten

all_measures$measure[which(all_measures$measureID %in% c(33))]
EffectOnTransmission(contactpattern, measures = all_measures[which(all_measures$measureID %in% c(33)),], relinf = relinf, locations = locations)

all_measures$measure[which(all_measures$measureID %in% c(33,9))]
EffectOnTransmission(contactpattern, measures = all_measures[which(all_measures$measureID %in% c(33,9)),], relinf = relinf, locations = locations)

all_measures$measure[which(all_measures$measureID %in% c(33,9,13))]
EffectOnTransmission(contactpattern, measures = all_measures[which(all_measures$measureID %in% c(33,9,13)),], relinf = relinf, locations = locations)

all_measures$measure[which(all_measures$measureID %in% c(33,9,13,23))]
EffectOnTransmission(contactpattern, measures = all_measures[which(all_measures$measureID %in% c(33,9,13,23)),], relinf = relinf, locations = locations)

all_measures$measure[which(all_measures$measureID %in% c(33,9,13,23,6))]
EffectOnTransmission(contactpattern, measures = all_measures[which(all_measures$measureID %in% c(33,9,13,23,6)),], relinf = relinf, locations = locations)


# Pakketten voor ___UITZONDERINGSGROND_2___
measures____UITZONDERINGSGROND_2___ <- all_measures %>% filter(measureID %in% c(99, "D1", "D2", "D3"))
EffectOnTransmission(contactpattern, measures = measures____UITZONDERINGSGROND_2___[1,], relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = measures____UITZONDERINGSGROND_2___[2,], relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = measures____UITZONDERINGSGROND_2___[3,], relinf = relinf, locations = locations)
EffectOnTransmission(contactpattern, measures = measures____UITZONDERINGSGROND_2___[4,], relinf = relinf, locations = locations)

eigenmax <- max(eigen(demo*tmp)$value)
contactmatrix___UITZONDERINGSGROND_2___ <- list()
for(i in 1:4) {
  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)
  tmp <- contactpattern2matrix(tmp, var.name = "c_smt")
  tmp <- (tmp + t(tmp))/2
  tmp <- tmp/eigenmax
  print(max(eigen(demo*tmp)$value))
  contactmatrix___UITZONDERINGSGROND_2___[[i]] <- tmp
}

names(contactmatrix___UITZONDERINGSGROND_2___) <- c("none","Work from home + limit social (-20%)", "Work from home + limit social (-40%)", "Work from home + limit social (-40%) + school closure")

names(contactmatrix___UITZONDERINGSGROND_2___) <- all_measures %>% filter(measureID %in% c(99, "D1", "D2", "D3")) %>% pull(measure)

contactmatrix___UITZONDERINGSGROND_2___[["none"]]

saveRDS(contactmatrix___UITZONDERINGSGROND_2___, "Scontactmatrix_4interventions_13mrt2020.rds")


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 <- relsus*t(relsus*(t(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)