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_26mei2020.r")
source("./___UITZONDERINGSGROND_6___modelContacts.r")

contact3 <- readRDS("../../nCoV/ContactAnalysis/data/contact3.rds") %>%
  mutate(cnt_location = as.character(cnt_location),
         cnt_location = ifelse(cnt_location == "Home", ifelse(hh_member == "Yes", "HomeHH", cnt_location), cnt_location),
         cnt_location = factor(cnt_location))
  
population3 <- readRDS("../../nCoV/ContactAnalysis/data/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)


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

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

# cnt_location EpiPose_wave1 Pienter_2017    frac
# <chr>                <dbl>        <dbl>   <dbl>
# 1 Home               0.566          0.911 0.621  
# 2 HomeHH             1.18           1.65  0.714  # this stays at 1
# 3 Leisure            0.0138         2.78  0.00497
# 4 Other place        0.721          1.92  0.376  
# 5 School             0.00966        0.440 0.0219 
# 6 Transport          0.0229         0.407 0.0563 
# 7 Work               0.464          4.20  0.111  

# contact matrices with uncertainty

fct <- 0 # to retain EpiPose reductions (0.26 residual mev)
fct <- 0.16 # to mimic reduction found in retrospecive Pienters (0.38 residual mev) by reducing reduction

fracH <- 0.621*(1-fct) + fct
fracW <- 0.111*(1-fct) + fct
fracS <- 0.0219*(1-fct) + fct
fracT <- 0.0563*(1-fct) + fct
fracL <- 0.00497*(1-fct) + fct
fracO <- 0.376*(1-fct) + fct


fct <- 1 # to retain EpiPose reductions (0.26 residual mev)
fct <- 2.1 # to mimic reduction found in retrospecive Pienters (0.38 residual mev) by increasing residual

fracH <- min(0.621*fct, 1)
fracW <- 0.111*fct
fracS <- 0.0219*fct
fracT <- 0.0563*fct
fracL <- 0.00497*fct
fracO <- 0.376*fct


# D3_leisure70
fracH <- 1.225^2
fracW <- 0.55^2
fracS <- 0.2^2
fracT <- 0.15^2
fracL <- 0.3^2
fracO <- 0.5^2

sqrt(c(fracH, fracW, fracS, fracT, fracL, fracO))


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

# D3 as EpiPose1
r <- runif(nsens)
r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
#for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := sqrt(fracH) - 0.1 + 0.2*r) # only for D3_leisure70!
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := sqrt(fracH) - 0.1 + 0.2*r)
r <- runif(nsens)
r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := sqrt(fracW) - 0.1 + 0.2*r)
r <- runif(nsens)
r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := sqrt(fracS) - 0.1 + 0.2*r)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := sqrt(fracT))
r <- runif(nsens)
r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := sqrt(fracL) - 0.1 + 0.2*r)
r <- runif(nsens)
r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := sqrt(fracO) - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus + batch1
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 - 0.1 + 0.2*r)
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.53 - 0.1 + 0.2*rwork)
r <- runif(nsens)
#r <- 0.5
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.61 - 0.1 + 0.2*r)
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.32 - 0.1 + 0.2*r)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.3 - 0.1 + 0.2*rleis)
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.16 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.102 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.53 - 0.1 + 0.2*rwork) + 0.35*((0.3 + 0.16 + 7*0.102)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 0.9 - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus + batch2
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 - 0.1 + 0.2*r)
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.63 - 0.1 + 0.2*rwork)
r <- runif(nsens)
#r <- 0.5
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) := 0.64 - 0.1 + 0.2*r)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.45 - 0.1 + 0.2*rleis)
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.42 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.30 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.63 - 0.1 + 0.2*rwork) + 0.35*((0.45 + 0.42 + 7*0.3)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 0.95 - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus + batch3
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 - 0.1 + 0.2*r)
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.63 - 0.1 + 0.2*rwork)
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.80 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.60 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.63 - 0.1 + 0.2*rwork) + 0.35*((2*0.8 + 7*0.6)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus from 1 sep (only schools open, as first approximation)
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 - 0.1 + 0.2*r)
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.63 - 0.1 + 0.2*rwork)
r <- runif(nsens)
#r <- 0.5
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) := 0.84 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*0.28 + (1-0.214)*0.52
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.80 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.60 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.63 - 0.1 + 0.2*rwork) + 0.35*((2*0.8 + 7*0.6)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus from 1 sep (schools open, more work, more leisure)
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 1 - 0.1 + 0.2*r)
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.77 - 0.1 + 0.2*rwork)
r <- runif(nsens)
#r <- 0.5
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) := 0.84 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*0.28 + (1-0.214)*0.52
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:4) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.80 - 0.1 + 0.2*rleis)
for(i in 5:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.60 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.77 - 0.1 + 0.2*rwork) + 0.35*((4*0.8 + 5*0.6)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 1 - 0.1 + 0.2*r)


# D3asEpiPose1_residualplus from 1 okt (schools open, less work, less leisure (bars early closing, max group sizes))
# SENSITIVITY home = 1 D3asEpiPose1_residualplus from 1 okt (schools open, less work, less leisure (bars early closing, max group sizes))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 0.90 - 0.1 + 0.2*r) # max 3 people at home outside HH
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.69 - 0.1 + 0.2*rwork) 
r <- runif(nsens)
#r <- 0.5
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) := 0.84 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*0.28 + (1-0.214)*0.52
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.60 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.30 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.69 - 0.1 + 0.2*rwork) + 0.35*((2*0.6 + 7*0.3)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 0.95 - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus from 1 okt with holiday (schools half open, less work for teachers (compared to from 1 okt without holiday))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 0.9 - 0.1 + 0.2*r) # max 3 people at home outside HH
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.66 - 0.1 + 0.2*rwork) 
r <- runif(nsens)
#r <- 0.5
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.61 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*0.5
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.53 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*(0.28 + 0.52)*0.5
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.60 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.30 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.66 - 0.1 + 0.2*rwork) + 0.35*((2*0.6 + 7*0.3)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 0.95 - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus from 1 okt EXTREME (schools open, even less work, even less leisure (bars closed, smaller max group sizes))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 0.90 - 0.1 + 0.2*r) # max 3 people at home outside HH
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.60 - 0.1 + 0.2*rwork) 
r <- runif(nsens)
#r <- 0.5
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) := 0.84 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*0.28 + (1-0.214)*0.52
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.60 - 0.1 + 0.2*rleis)
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.40 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.20 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.60 - 0.1 + 0.2*rwork) + 0.35*((0.6 + 0.4 + 7*0.2)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 0.95 - 0.1 + 0.2*r)

# D3asEpiPose1_residualplus from 1 okt EXTREME with holiday (schools half open, even less work, even less leisure (bars closed, smaller max group sizes))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("HH_", i) := 1) # no changes to contacts with household members at home
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("H_", i) := 0.90 - 0.1 + 0.2*r) # max 3 people at home outside HH
rwork <- runif(nsens)
#rwork <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("W_", i) := 0.57 - 0.1 + 0.2*rwork) 
r <- runif(nsens)
#r <- 0.5
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.61 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*0.5
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.53 - 0.1 + 0.2*r) # 0.214 + (1-0.214)*(0.28 + 0.52)*0.5
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("S_", i) := 0.214 - 0.1 + 0.2*r)
rleis <- runif(nsens)
#rleis <- 0.5
for(i in 1:1) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.60 - 0.1 + 0.2*rleis)
for(i in 2:2) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.40 - 0.1 + 0.2*rleis)
for(i in 3:9) measures_sens <- measures_sens %>% mutate(!!paste0("L_", i) := 0.20 - 0.1 + 0.2*rleis)
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("T_", i) := 0.65*(0.57 - 0.1 + 0.2*rwork) + 0.35*((0.6 + 0.4 + 7*0.2)/9 - 0.1 + 0.2*rleis))
r <- runif(nsens)
#r <- 0.5
for(i in 1:9) measures_sens <- measures_sens %>% mutate(!!paste0("O_", i) := 0.95 - 0.1 + 0.2*r)



# i <- 1
# tmp <- EffectOfMeasureOnContacts(contactpattern = contactpattern, par = measures_sens[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
# max(eigen(demo*tmp)$value)

# i <- 1
# tmp <- EffectOfMeasureOnContacts(contactpattern = contactpattern, par = measures_sens[i,])
# tmp$n_cnt <- rowSums(tmp %>% select(paste0("n_cnt_", locations)))
# tmp <- tmp %>%
#  select(-n_cnt, -m_obs, -c_obs, -m_smt, -c_smt) %>%
#  pivot_longer(cols = grep(names(.), pattern = "n_cnt_", value = TRUE), names_to = "location", values_to = "n_cnt") %>%
#  mutate(location = gsub(location, pattern = "n_cnt_", replacement = ""),
#         m_obs = n_cnt/n_part,
#         c_obs = m_obs/f_pop)
# saveRDS(tmp, "contactpattern_D3_leisure70.rds")

contactpattern2matrix(contactpattern, var.name = "c_smt")

locations <- levels(contact3$cnt_location)
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)

saveRDS(contactmatrix, "Contactmatrix_D3asEpiPose1_residualplus_extremeautumn_holiday_8okt2020.rds")


# contactmatrix for ___UITZONDERINGSGROND_2___
tmp <- contactpattern2matrix(contactpattern, var.name = "c_smt")
tmp <- (tmp + t(tmp))/2
# scale with eigenmax of full contact matrix (note: no relinf, just the contacts)
tmp <- tmp/max(eigen(demo*tmp)$value)
saveRDS(tmp, file = "contactmatrix2017_9x9_voor___UITZONDERINGSGROND_2___.rds")


0.483^2
sqrt(0.5)
0.607^2

0.214+(1-0.214)*0.5
0.214+(1-0.214)*0.28*0.5


(0.2*0.214 + 0.8)
sqrt(0.2*0.214^2 + 0.8)
0.8^2

sqrt(0.8 + 0.2*0.214^2)
sqrt(0.5 + 0.5*0.214^2)

sqrt(0.5*0.28 + (1-0.5*0.28)*0.214^2)

5/40
sqrt(0.3 + 0.7*0.102^2)
0.65*0.53 + 0.35*(0.3 + 0.16 + 7*0.102)/9

0.102 + 0.2*(1-0.102)
0.214 + (1-0.214)*0.28 + (1-0.214)*0.52*0.5
0.64*0.64

0.102 + 0.30*(1-0.102)
0.72*0.102 + 0.28*0.3
0.8*0.45 + 0.2*0.3

0.65*0.63 + 0.35*(0.45 + 0.42 + 7*0.3)/9