final.size <- function(R0, contact.mat) {
  G <- contact.mat
  G <- G/max(eigen(G)$values) # normalised to max eigenvalue
  s0 <- rep(1, nrow(contact.mat))
  
  fnr <- function(zzs) s0*(1-exp(-R0*apply(G*zzs, MARGIN = 2, sum))) - zzs
  tG <- t(G)
  zs <- s0*(1-exp(-R0*apply(G*s0, MARGIN = 2, sum)))
  fz <- fnr(zs)
  
  ctr <- 1
  while(max(abs(fz)) > 1e-09 & ctr < 200) {
    tmp <- apply(G*zs, MARGIN = 2, sum)
    Jac <- R0*s0*exp(-R0*tmp) 
    Jac <- tG*Jac - diag(nrow(tG))
    
    invJ <- solve(Jac)
    zs <- zs - as.vector(invJ%*%fz)
    fz <- fnr(zs)
    ctr <- ctr + 1
  }
  
  return(attack = zs)
}

contactpattern2matrix <- function(contactpattern.data, var.name = "c_smt") {
  sex <- contactpattern.data %>% names %>% str_detect("gender") %>% any
  if(sex) {
    contactpattern.data <- contactpattern.data %>%   
      select(part_gender, cnt_gender, part_age, cnt_age, var.name) %>%   
      mutate(part = paste(part_gender, part_age),
             cnt = paste(cnt_gender, cnt_age)) %>%
      select(-part_age, - part_gender, -cnt_age, -cnt_gender) 
  } else {
    contactpattern.data <- contactpattern.data %>% 
      select(part_age, cnt_age, var.name) %>%
      mutate(part = part_age,
             cnt = cnt_age) %>%
      select(-part_age, -cnt_age) 
  }
  part_order <- contactpattern.data$part %>% unique
  contact.mat <- contactpattern.data %>% spread(key = part, value = var.name)
  rownames <- contact.mat$cnt
  contact.mat <- contact.mat %>% select(-cnt) %>% as.matrix 
  rownames(contact.mat) <- rownames
  contact.mat <- contact.mat[part_order, part_order]
  return(t(contact.mat))
}

countContactsLoc <- function(contact.data, 
                          pop.data,
                          loc = NULL,
                          age.breaks = c(seq(from = 0, to = 99, by = 1), Inf),
                          sex = TRUE) {
  
  # check if cnt_age_cor available
  if(!("cnt_age_cor" %in% names(contact.data))) stop("Variable ", "\"", "cnt_age_cor", "\"", " not found: run ", "\"", "correctContactAgeDigitPreference", "\""," first")
  
  # Count participants, contacts and population ----
  
  # Participants
  part_count.data <- contact.data %>%
    filter(!is.na(cnt_location)) %>% # ignore 108 contacts without location
    # Filter out duplicated participants (i.e. the contacts), because we are counting participants here
    filter(!duplicated(part_id)) %>% 
    # Cut part_age in provided age categories, open on the right [left, right)
    mutate(part_age = part_age %>% cut(breaks = age.breaks,right = FALSE, include.lowest = TRUE)) %>% 
    # Group by age categories, drop none if they appear to be empty
    group_by(part_age, .drop = FALSE) %>% 
    # Group by part_gender if sex = TRUE
    {if (sex) group_by(., part_gender, add = TRUE, .drop = FALSE) else .} %>% 
    # Count number of participants
    count(name = "n_part")

  # Contacts
  cont_count.data <- contact.data %>%
    filter(!is.na(cnt_location)) %>% # ignore 108 contacts without location
    {if(length(loc) != 0) filter(., cnt_location == loc) else .} %>%
    # Filter out missing cnt_age_cor or cnt_gender. These contacts are useless
    # Then drop empty level "Missing" from cnt_gender
    filter(!is.na(cnt_age_cor) & cnt_gender != "(Missing)") %>% 
    droplevels %>% 
    # Cut part_age and cont_age_cor in provided age categories, open on the right [left, right)
    mutate(
      part_age = part_age    %>% cut(breaks = age.breaks, right = FALSE, include.lowest = TRUE),
      cnt_age  = cnt_age_cor %>% cut(breaks = age.breaks, right = FALSE, include.lowest = TRUE)) %>% 
    # Group by age categories, drop none if they appear to be empty
    group_by(part_age, cnt_age, .drop = FALSE) %>% 
    # Group by part_gender and cnt_gender if sex = TRUE
    {if (sex) group_by(., part_gender, cnt_gender, add = TRUE, .drop = FALSE) else .} %>% 
    # Count number of contacts
    count(name = "n_cnt")
  
  # Population
  pop_count.data <- pop.data %>% 
    # Cut age in provided age categories, open on the right [left, right)
    mutate(age = age %>% cut(breaks = age.breaks, right = FALSE, include.lowest = TRUE)) %>% 
    # Group by age categories, drop none if they appear to be empty
    group_by(age, .drop = FALSE) %>% 
    # Group by gender if sex = TRUE, and make sure gender factor levels from contact.data and pop.data are aligned
    {if (sex) group_by(., gender = gender %>% fct_relevel(levels(part_count.data$part_gender)), add = TRUE, .drop = FALSE)
      else .} %>% 
    # Count the population number
    summarize(n_pop = sum(population)) %>%
    ungroup %>%
    # Calculate population fraction
    mutate(f_pop = n_pop/sum(n_pop))
  
  # Join the tibbles above into one tibble ----
  # Calculate observed contact intensities and rates
  
  # Create contactpattern.data by joining part_cont.data and cont_count.data
  contactpattern.data <- 
    full_join(part_count.data, cont_count.data) %>%
    # The next three steps depend on sex. Therefore we use an if statemant
    # 1. Join with pop_count.data
    #    Note: because the number of contacts is proportional with population,
    #          explicitly join age and gender in pop_count.data with that in cont_count.data
    # 2. For modelling the contact rates, be sure the records are arranged as given
    # 3. Reorder columns as given
    {if (sex) {
      full_join(., pop_count.data, by = c("cnt_age" = "age", "cnt_gender" = "gender")) %>% 
        arrange(part_gender, cnt_gender, part_age, cnt_age) %>% 
        select(part_gender, cnt_gender, part_age, cnt_age, n_part, n_pop, f_pop, n_cnt)
    } else {
      full_join(., pop_count.data, by = c("cnt_age" = "age")) %>% 
        arrange(part_age, cnt_age) %>% 
        select(part_age, cnt_age, n_part, n_pop, f_pop, n_cnt)
    }} %>%
    mutate(
      # Calculate observed contact intensties (m_obs) and rates (c_obs)
      # Note: based on population fractions and not population numbers!
      m_obs = n_cnt/n_part,
      c_obs = m_obs/f_pop) %>% 
    # Ungroup data
    ungroup
  
  # for age categories of 1 year, use numbers for participant and contact age
  if(all(diff(age.breaks) == Inf | diff(age.breaks) == 1)) {
    levels(contactpattern.data$part_age) <- age.breaks[-length(age.breaks)]
    levels(contactpattern.data$cnt_age) <- age.breaks[-length(age.breaks)]
  }
  
  return(contactpattern.data)  
}

countExposureContactsLoc <- function(contact.data, 
                                  pop.data,
                                  loc = NULL,
                                  age.breaks = c(seq(from = 0, to = 80, by = 10), Inf),
                                  sex = FALSE) {
  
  # check if cnt_age_cor available
  if(!("cnt_age_cor" %in% names(contact.data))) stop("Variable ", "\"", "cnt_age_cor", "\"", " not found: run ", "\"", "correctContactAgeDigitPreference", "\""," first")
  
  # Count participants, contacts and population ----
  
  # Participants
  part_count.data <- contact.data %>%
    filter(!is.na(cnt_location)) %>% # ignore 108 contacts without location
    # Filter out duplicated participants (i.e. the contacts), because we are counting participants here
    filter(!duplicated(part_id)) %>% 
    # Cut part_age in provided age categories, open on the right [left, right)
    mutate(part_age = part_age %>% cut(breaks = age.breaks,right = FALSE, include.lowest = TRUE)) %>% 
    # Group by age categories, drop none if they appear to be empty
    group_by(part_age, .drop = FALSE) %>% 
    # Group by part_gender if sex = TRUE
    {if (sex) group_by(., part_gender, add = TRUE, .drop = FALSE) else .} %>% 
    # Count number of participants
    count(name = "n_part")
  
  # Contacts
  cont_count.data <- contact.data %>%
    filter(!is.na(cnt_location)) %>% # ignore 108 contacts without location
    {if(length(loc) != 0) filter(., cnt_location == loc) else .} %>%
    # Filter out missing cnt_age_cor or cnt_gender. These contacts are useless
    # Then drop empty level "Missing" from cnt_gender
    filter(!is.na(cnt_age_cor) & cnt_gender != "(Missing)") %>% 
    droplevels %>% 
    # Cut part_age and cont_age_cor in provided age categories, open on the right [left, right)
    mutate(
      duration = duration_multi %>% fct_recode("5" = "less than 5 minutes", 
                                               "10" = "5 - 15 mins",
                                               "30" = "15 mins - 1 hour",
                                               "120" =  "1 hour - 4 hours",
                                               "240" = "4 hours or more",
                                               "60" = "(Missing)") %>% as.character %>% as.numeric,
      part_age = part_age    %>% cut(breaks = age.breaks, right = FALSE, include.lowest = TRUE),
      cnt_age  = cnt_age_cor %>% cut(breaks = age.breaks, right = FALSE, include.lowest = TRUE)) %>% 
    # Group by age categories, drop none if they appear to be empty
    group_by(part_age, cnt_age, .drop = FALSE) %>% 
    # Group by part_gender and cnt_gender if sex = TRUE
    {if (sex) group_by(., part_gender, cnt_gender, add = TRUE, .drop = FALSE) else .} %>% 
    # Count number of contacts
    # count(name = "n_cnt")
    summarise(n_cnt = sum(duration)/60)
  
  # Population
  pop_count.data <- pop.data %>% 
    # Cut age in provided age categories, open on the right [left, right)
    mutate(age = age %>% cut(breaks = age.breaks, right = FALSE, include.lowest = TRUE)) %>% 
    # Group by age categories, drop none if they appear to be empty
    group_by(age, .drop = FALSE) %>% 
    # Group by gender if sex = TRUE, and make sure gender factor levels from contact.data and pop.data are aligned
    {if (sex) group_by(., gender = gender %>% fct_relevel(levels(part_count.data$part_gender)), add = TRUE, .drop = FALSE)
      else .} %>% 
    # Count the population number
    summarize(n_pop = sum(population)) %>%
    ungroup %>%
    # Calculate population fraction
    mutate(f_pop = n_pop/sum(n_pop))
  
  # Join the tibbles above into one tibble ----
  # Calculate observed contact intensities and rates
  
  # Create contactpattern.data by joining part_cont.data and cont_count.data
  contactpattern.data <- 
    full_join(part_count.data, cont_count.data) %>%
    # The next three steps depend on sex. Therefore we use an if statemant
    # 1. Join with pop_count.data
    #    Note: because the number of contacts is proportional with population,
    #          explicitly join age and gender in pop_count.data with that in cont_count.data
    # 2. For modelling the contact rates, be sure the records are arranged as given
    # 3. Reorder columns as given
    {if (sex) {
      full_join(., pop_count.data, by = c("cnt_age" = "age", "cnt_gender" = "gender")) %>% 
        arrange(part_gender, cnt_gender, part_age, cnt_age) %>% 
        select(part_gender, cnt_gender, part_age, cnt_age, n_part, n_pop, f_pop, n_cnt)
    } else {
      full_join(., pop_count.data, by = c("cnt_age" = "age")) %>% 
        arrange(part_age, cnt_age) %>% 
        select(part_age, cnt_age, n_part, n_pop, f_pop, n_cnt)
    }} %>%
    mutate(
      # Calculate observed contact intensties (m_obs) and rates (c_obs)
      # Note: based on population fractions and not population numbers!
      m_obs = n_cnt/n_part,
      c_obs = m_obs/f_pop) %>% 
    # Ungroup data
    ungroup
  
  # for age categories of 1 year, use numbers for participant and contact age
  if(all(diff(age.breaks) == Inf | diff(age.breaks) == 1)) {
    levels(contactpattern.data$part_age) <- age.breaks[-length(age.breaks)]
    levels(contactpattern.data$cnt_age) <- age.breaks[-length(age.breaks)]
  }
  
  return(contactpattern.data)  
}

ContactPatternLocation <- function(contact.data, pop.data, age.breaks = c(seq(from = 0, to = 80, by = 10), Inf), sex = FALSE, countWhat = countContactsLoc) {
  contactpattern <- countWhat(contact.data = contact.data,
                              pop.data = pop.data,
                              age.breaks = age.breaks,
                              sex = sex)
  contactpattern <- modelContacts(contactpattern.data = contactpattern, posterior.sample = 0, verbose = FALSE)
  for(loc in levels(contact.data$cnt_location)) {
    # tmp <- contact.data %>% filter(cnt_location == loc)
    tmp <- countWhat(contact.data = contact.data,
                     pop.data = pop.data,
                     loc = loc,
                     age.breaks = age.breaks,
                     sex = sex)
    # tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = TRUE)
    contactpattern[[paste0("n_cnt_",loc)]] <- tmp$n_cnt
  }
  # contactpattern[["m_smt_tot"]] <- rowSums(contactpattern %>% select(paste0("m_smt_", levels(contact3$cnt_location))))
  return(contactpattern)
}

ContactPatternLocation_old <- function(contact.data, pop.data, age.breaks = c(seq(from = 0, to = 80, by = 10), Inf), sex = FALSE, countWhat = countContactsLoc) {
  contactpattern <- countWhat(contact.data = contact.data,
                              pop.data = pop.data,
                              age.breaks = age.breaks,
                              sex = sex)
  contactpattern <- modelContacts(contactpattern.data = contactpattern, posterior.sample = 0, verbose = FALSE)
  for(loc in levels(contact.data$cnt_location)) {
    # tmp <- contact.data %>% filter(cnt_location == loc)
    tmp <- countWhat(contact.data = contact.data,
                     pop.data = pop.data,
                     loc = loc,
                     age.breaks = age.breaks,
                     sex = sex)
    tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = TRUE)
    contactpattern[[paste0("m_smt_",loc)]] <- tmp$m_smt
  }
  contactpattern[["m_smt_tot"]] <- rowSums(contactpattern %>% select(paste0("m_smt_", levels(contact3$cnt_location))))
  return(contactpattern)
}

EffectOfMeasureOnContacts <- function(contactpattern, par = NULL) {
  
  effect_Home <- if("effect_Home" %in% names(par)) par$effect_Home else 1
  min_Home <- if("min_Home" %in% names(par)) par$min_Home else 1
  max_Home <- if("max_Home" %in% names(par)) par$max_Home else 9
  effect_Work <- if("effect_Work" %in% names(par)) par$effect_Work else 1
  min_Work <- if("min_Work" %in% names(par)) par$min_Work else 1
  max_Work <- if("max_Work" %in% names(par)) par$max_Work else 9
  effect_School <- if("effect_School" %in% names(par)) par$effect_School else 1
  min_School <- if("min_School" %in% names(par)) par$min_School else 1
  max_School <- if("max_School" %in% names(par)) par$max_School else 9
  effect_Transport <- if("effect_Transport" %in% names(par)) par$effect_Transport else 1
  min_Transport <- if("min_Transport" %in% names(par)) par$min_Transport else 1
  max_Transport <- if("max_Transport" %in% names(par)) par$max_Transport else 9
  effect_Leisure <- if("effect_Leisure" %in% names(par)) par$effect_Leisure else 1
  min_Leisure <- if("min_Leisure" %in% names(par)) par$min_Leisure else 1
  max_Leisure <- if("max_Leisure" %in% names(par)) par$max_Leisure else 9
  effect_Other <- if("effect_Other" %in% names(par)) par$effect_Other else 1
  min_Other <- if("min_Other" %in% names(par)) par$min_Other else 1
  max_Other <- if("max_Other" %in% names(par)) par$max_Other else 9
  
  tmp <- contactpattern %>% 
    # effect on contact rate
    mutate(n_cnt_Home = ifelse(as.integer(part_age) >= min_Home & as.integer(part_age) <= max_Home, effect_Home*n_cnt_Home, n_cnt_Home),
           n_cnt_Home = ifelse(as.integer(cnt_age) >= min_Home & as.integer(cnt_age) <= max_Home, effect_Home*n_cnt_Home, n_cnt_Home),
           n_cnt_Work = ifelse(as.integer(part_age) >= min_Work & as.integer(part_age) <= max_Work, effect_Work*n_cnt_Work, n_cnt_Work),
           n_cnt_Work = ifelse(as.integer(cnt_age) >= min_Work & as.integer(cnt_age) <= max_Work, effect_Work*n_cnt_Work, n_cnt_Work),
           n_cnt_School = ifelse(as.integer(part_age) >= min_School & as.integer(part_age) <= max_School, effect_School*n_cnt_School, n_cnt_School),
           n_cnt_School = ifelse(as.integer(cnt_age) >= min_School & as.integer(cnt_age) <= max_School, effect_School*n_cnt_School, n_cnt_School),
           n_cnt_Transport = ifelse(as.integer(part_age) >= min_Transport & as.integer(part_age) <= max_Transport, effect_Transport*n_cnt_Transport, n_cnt_Transport),
           n_cnt_Transport = ifelse(as.integer(cnt_age) >= min_Transport & as.integer(cnt_age) <= max_Transport, effect_Transport*n_cnt_Transport, n_cnt_Transport),
           n_cnt_Leisure = ifelse(as.integer(part_age) >= min_Leisure & as.integer(part_age) <= max_Leisure, effect_Leisure*n_cnt_Leisure, n_cnt_Leisure),
           n_cnt_Leisure = ifelse(as.integer(cnt_age) >= min_Leisure & as.integer(cnt_age) <= max_Leisure, effect_Leisure*n_cnt_Leisure, n_cnt_Leisure),
           `n_cnt_Other place` = ifelse(as.integer(part_age) >= min_Other & as.integer(part_age) <= max_Other, effect_Other*`n_cnt_Other place`, `n_cnt_Other place`),
           `n_cnt_Other place` = ifelse(as.integer(cnt_age) >= min_Other & as.integer(cnt_age) <= max_Other, effect_Other*`n_cnt_Other place`, `n_cnt_Other place`))
  #%>%
  #tmp$m_obs_tot <- rowSums(tmp %>% select(paste0("m_obs_", levels(contact3$cnt_location))))
  
  return(tmp)
}

EffectOfMeasureOnInfSus <- function(mat, par = NULL) {
 
  Ieffect <- if("Ieffect" %in% names(par)) par$Ieffect else 1
  Imin <- if("Imin" %in% names(par)) par$Imin else 1
  Imax <- if("Imax" %in% names(par)) par$Imax else 9
  
  Ieff <- rep(1, nrow(mat))
  Ieff[Imin:Imax] <- Ieff[Imin:Imax]*Ieffect
  
  Seffect <- if("Seffect" %in% names(par)) par$Seffect else 1
  Smin <- if("Smin" %in% names(par)) par$Smin else 1
  Smax <- if("Smax" %in% names(par)) par$Smax else 9
  
  Seff <- rep(1, nrow(mat))
  Seff[Smin:Smax] <- Seff[Smin:Smax]*Seffect
  
  return(t(Seff*t(Ieff*mat)))
}

EffectOnTransmission <- function(contactpattern, measures, relinf, locations) {
  
  eigenmax <- max(eigen(relinf*contactpattern2matrix(contactpattern.data = contactpattern, var.name = "m_smt"))$value)
  
  tmp <- contactpattern
  for(i in 1:nrow(measures)) {
    input <- measures[i,] #%>% select(-measure, -measureID)
    tmp <- EffectOfMeasureOnContacts(tmp, par = input)
  }
  tmp$n_cnt <- rowSums(tmp %>% select(paste0("n_cnt_", locations)))
  #tmp$c_obs <- tmp$m_obs/tmp$f_pop
  tmp <- modelContacts(contactpattern.data = tmp, posterior.sample = 0, verbose = TRUE)
  tmp <- relinf*contactpattern2matrix(tmp, var.name = "m_smt")
  
  for(i in 1:nrow(measures)) {
    input <- measures[i,] #%>% select(-measure, -measureID)
    tmp <- EffectOfMeasureOnInfSus(tmp, par = input)
  }
  
  neweigenmax <- max(eigen(tmp)$value)
  
  return(100*(neweigenmax - eigenmax)/eigenmax)
}

InfectiousContactPattern <- function(contactpattern, locations, relinf) {
  for(loc in locations) {
    contactpattern[[paste0("m_smt_inf_", loc)]] <- as.vector(t(relinf*contactpattern2matrix(contactpattern, var.name = paste0("m_smt_", loc))))
  }
  contactpattern[["m_smt_inf"]] <- as.vector(t(relinf*contactpattern2matrix(contactpattern, var.name = "m_smt")))
  contactpattern[["m_smt_inf_tot"]] <- as.vector(t(relinf*contactpattern2matrix(contactpattern, var.name = "m_smt_tot")))
  return(contactpattern)
}

EffectOfMeasure_old <- function(contactpattern, par = NULL) {
  
  effect_Home <- if("effect_Home" %in% names(par)) par$effect_Home else 1
  min_Home <- if("min_Home" %in% names(par)) par$min_Home else 1
  max_Home <- if("max_Home" %in% names(par)) par$max_Home else 9
  effect_Work <- if("effect_Work" %in% names(par)) par$effect_Work else 1
  min_Work <- if("min_Work" %in% names(par)) par$min_Work else 1
  max_Work <- if("max_Work" %in% names(par)) par$max_Work else 9
  effect_School <- if("effect_School" %in% names(par)) par$effect_School else 1
  min_School <- if("min_School" %in% names(par)) par$min_School else 1
  max_School <- if("max_School" %in% names(par)) par$max_School else 9
  effect_Transport <- if("effect_Transport" %in% names(par)) par$effect_Transport else 1
  min_Transport <- if("min_Transport" %in% names(par)) par$min_Transport else 1
  max_Transport <- if("max_Transport" %in% names(par)) par$max_Transport else 9
  effect_Leisure <- if("effect_Leisure" %in% names(par)) par$effect_Leisure else 1
  min_Leisure <- if("min_Leisure" %in% names(par)) par$min_Leisure else 1
  max_Leisure <- if("max_Leisure" %in% names(par)) par$max_Leisure else 9
  effect_Other <- if("effect_Other" %in% names(par)) par$effect_Other else 1
  min_Other <- if("min_Other" %in% names(par)) par$min_Other else 1
  max_Other <- if("max_Other" %in% names(par)) par$max_Other else 9
  
  Ieffect_Home <- if("Ieffect_Home" %in% names(par)) par$Ieffect_Home else 1
  Imin_Home <- if("Imin_Home" %in% names(par)) par$Imin_Home else 1
  Imax_Home <- if("Imax_Home" %in% names(par)) par$Imax_Home else 9
  Ieffect_Work <- if("Ieffect_Work" %in% names(par)) par$Ieffect_Work else 1
  Imin_Work <- if("Imin_Work" %in% names(par)) par$Imin_Work else 1
  Imax_Work <- if("Imax_Work" %in% names(par)) par$Imax_Work else 9
  Ieffect_School <- if("Ieffect_School" %in% names(par)) par$Ieffect_School else 1
  Imin_School <- if("Imin_School" %in% names(par)) par$Imin_School else 1
  Imax_School <- if("Imax_School" %in% names(par)) par$Imax_School else 9
  Ieffect_Transport <- if("Ieffect_Transport" %in% names(par)) par$Ieffect_Transport else 1
  Imin_Transport <- if("Imin_Transport" %in% names(par)) par$Imin_Transport else 1
  Imax_Transport <- if("Imax_Transport" %in% names(par)) par$Imax_Transport else 9
  Ieffect_Leisure <- if("Ieffect_Leisure" %in% names(par)) par$Ieffect_Leisure else 1
  Imin_Leisure <- if("Imin_Leisure" %in% names(par)) par$Imin_Leisure else 1
  Imax_Leisure <- if("Imax_Leisure" %in% names(par)) par$Imax_Leisure else 9
  Ieffect_Other <- if("Ieffect_Other" %in% names(par)) par$Ieffect_Other else 1
  Imin_Other <- if("Imin_Other" %in% names(par)) par$Imin_Other else 1
  Imax_Other <- if("Imax_Other" %in% names(par)) par$Imax_Other else 9
  
  Seffect_Home <- if("Seffect_Home" %in% names(par)) par$Seffect_Home else 1
  Smin_Home <- if("Smin_Home" %in% names(par)) par$Smin_Home else 1
  Smax_Home <- if("Smax_Home" %in% names(par)) par$Smax_Home else 9
  Seffect_Work <- if("Seffect_Work" %in% names(par)) par$Seffect_Work else 1
  Smin_Work <- if("Smin_Work" %in% names(par)) par$Smin_Work else 1
  Smax_Work <- if("Smax_Work" %in% names(par)) par$Smax_Work else 9
  Seffect_School <- if("Seffect_School" %in% names(par)) par$Seffect_School else 1
  Smin_School <- if("Smin_School" %in% names(par)) par$Smin_School else 1
  Smax_School <- if("Smax_School" %in% names(par)) par$Smax_School else 9
  Seffect_Transport <- if("Seffect_Transport" %in% names(par)) par$Seffect_Transport else 1
  Smin_Transport <- if("Smin_Transport" %in% names(par)) par$Smin_Transport else 1
  Smax_Transport <- if("Smax_Transport" %in% names(par)) par$Smax_Transport else 9
  Seffect_Leisure <- if("Seffect_Leisure" %in% names(par)) par$Seffect_Leisure else 1
  Smin_Leisure <- if("Smin_Leisure" %in% names(par)) par$Smin_Leisure else 1
  Smax_Leisure <- if("Smax_Leisure" %in% names(par)) par$Smax_Leisure else 9
  Seffect_Other <- if("Seffect_Other" %in% names(par)) par$Seffect_Other else 1
  Smin_Other <- if("Smin_Other" %in% names(par)) par$Smin_Other else 1
  Smax_Other <- if("Smax_Other" %in% names(par)) par$Smax_Other else 9
  
  tmp <- contactpattern %>% 
    # effect on contact rate
    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`)) %>%
    # effect on infectiousness infectors (participants)
    mutate(m_smt_inf_Home = ifelse(as.integer(part_age) >= Imin_Home & as.integer(part_age) <= Imax_Home, Ieffect_Home*m_smt_inf_Home, m_smt_inf_Home),
           m_smt_inf_Work = ifelse(as.integer(part_age) >= Imin_Work & as.integer(part_age) <= Imax_Work, Ieffect_Work*m_smt_inf_Work, m_smt_inf_Work),
           m_smt_inf_School = ifelse(as.integer(part_age) >= Imin_School & as.integer(part_age) <= Imax_School, Ieffect_School*m_smt_inf_School, m_smt_inf_School),
           m_smt_inf_Transport = ifelse(as.integer(part_age) >= Imin_Transport & as.integer(part_age) <= Imax_Transport, Ieffect_Transport*m_smt_inf_Transport, m_smt_inf_Transport),
           m_smt_inf_Leisure = ifelse(as.integer(part_age) >= Imin_Leisure & as.integer(part_age) <= Imax_Leisure, Ieffect_Leisure*m_smt_inf_Leisure, m_smt_inf_Leisure),
           `m_smt_inf_Other place` = ifelse(as.integer(part_age) >= Imin_Other & as.integer(part_age) <= Imax_Other, Ieffect_Other*`m_smt_inf_Other place`, `m_smt_inf_Other place`)) %>%
    # effect on susceptibility infectees (contacts)
    mutate(m_smt_inf_Home = ifelse(as.integer(cnt_age) >= Smin_Home & as.integer(cnt_age) <= Smax_Home, Seffect_Home*m_smt_inf_Home, m_smt_inf_Home),
           m_smt_inf_Work = ifelse(as.integer(cnt_age) >= Smin_Work & as.integer(cnt_age) <= Smax_Work, Seffect_Work*m_smt_inf_Work, m_smt_inf_Work),
           m_smt_inf_School = ifelse(as.integer(cnt_age) >= Smin_School & as.integer(cnt_age) <= Smax_School, Seffect_School*m_smt_inf_School, m_smt_inf_School),
           m_smt_inf_Transport = ifelse(as.integer(cnt_age) >= Smin_Transport & as.integer(cnt_age) <= Smax_Transport, Seffect_Transport*m_smt_inf_Transport, m_smt_inf_Transport),
           m_smt_inf_Leisure = ifelse(as.integer(cnt_age) >= Smin_Leisure & as.integer(cnt_age) <= Smax_Leisure, Seffect_Leisure*m_smt_inf_Leisure, m_smt_inf_Leisure),
           `m_smt_inf_Other place` = ifelse(as.integer(cnt_age) >= Smin_Other & as.integer(cnt_age) <= Smax_Other, Seffect_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))))
  
  return(tmp)
}

EffectOnTransmission_old <- function(contactpattern, measures) {
  
  eigenmax <- max(eigen(contactpattern2matrix(contactpattern.data = contactpattern, var.name = "m_smt_inf_tot"))$value)
  
  tmp <- contactpattern
  for(i in 1:nrow(measures)) {
    input <- measures[i,] #%>% select(-measure, -measureID)
    tmp <- EffectOfMeasure(tmp, par = input)
  }
  neweigenmax <- max(eigen(contactpattern2matrix(tmp, var.name = "m_smt_inf_tot"))$value)
  
  return(100*(neweigenmax - eigenmax)/eigenmax)
}

plotContactMatrix <- function(contactpattern.data, var.name = "c_smt", n.ticks = 20, plot.contour = FALSE, min = NA, max = NA) {
  
  sex <- contactpattern.data %>% names %>% str_detect("gender") %>% any
  
  n.levels <- length(levels(contactpattern.data$part_age))
  which2plot <- rep(FALSE, n.levels)
  which2plot[seq(1, n.levels, by = ceiling(n.levels/n.ticks))] <- TRUE
  
  rate <- str_detect(var.name, pattern = "c_")
  smoothed <- str_detect(var.name, pattern = "smt")
  
  data = contactpattern.data %>% 
    mutate(variable = contactpattern.data[[var.name]])
  
  log10max <- if(is.na(max)) ceiling(log10(max(data$variable, na.rm = TRUE))) else log10(max)
  log10min <- if(is.na(min)) floor(log10(min(data$variable[data$variable > 0], na.rm = TRUE))) else log10(min)
  
  p <- ggplot(
    data = data,
    mapping = aes(x = part_age, y = cnt_age, fill = variable, z = variable)) +
    geom_tile() +
    # scale_fill_distiller(
    #   breaks = c(0, outer(c(1, 2, 5), 10^(log10min:log10max))),
    #   palette = "YlOrRd",
    #   na.value = NA,
    #   trans = if(smoothed) "log" else "sqrt",
    #   direction = 1) +
    scale_fill_viridis_c(
      limits = 10^c(log10min, log10max),
      breaks = c(0, outer(c(1, 2, 5), 10^(floor(log10min):ceiling(log10max)))),
      na.value = "#450256",
      trans = if(smoothed) "log" else "sqrt",
      direction = 1) +
    coord_equal() +
    labs(
      x = "Participant age",
      y = "Contact age",
      #fill = if(rate) expression("Contact\nrate") else expression("Contact\nintensity")) +
      fill = if(rate) expression("Contact\nrate") else expression("Number of\ncontacts per\nparticipant")) +
    scale_x_discrete(
      breaks = levels(contactpattern.data$part_age)[which2plot]) +
    scale_y_discrete(
      breaks = levels(contactpattern.data$cnt_age)[which2plot]) +
      {if(sex) facet_grid(
        rows = vars(cnt_gender),
        cols = vars(part_gender %>% fct_rev))} +
    theme_bw() +
    theme(
      legend.position = "none",
      panel.grid = element_blank(),
      legend.key.height = unit(0.09, "npc"),
      axis.text.x = element_text(
        angle = 90,
        hjust = 1,
        vjust = 0.5))
  
  if(plot.contour) {
    p <- p + 
      geom_contour(
        aes(x = as.integer(part_age), 
            y = as.integer(cnt_age)),
        size = 0.1,
        colour = "black",
        breaks = c(outer(c(1, 2, 5), 10^(log10min:log10max))))
  }
  return(p)
}