##################################################
### contacts, infectivities,                   ###
###   susceptibilities, symptomatic fractions, ###
###   natural history and control              ###
###     v1: SEEIIR-model instead of SIR        ###
##################################################

### contact matrices
# contactmatrix of contact hours (Pienter3), and all control matrices. "99" is < 13 March (no contact reduction)
ctrlmats <- if(!exists("ctrlmats")) readRDS("data/ContactmatricesD3praktijk_midpoint_24mrt2020.rds") else ctrlmats
contactmatrix <- function(scenario = "99", nr = 1, altscenario = "99", ...) {
  if(scenario %in% names(ctrlmats)) {
    return(ctrlmats[[scenario]])
  } else if(nr == 0) {
    return(ctrlmats[["D3praktijk_leisure70"]])
  } else {
    return(ctrlmats[[nr]])
  }
} 

### from contact matrix to case age distribution
# normaliseren via (1) probI2S, (2) relsusceptibility, (3) probI2S + relsusceptibility,
# (4) relinfectivity + probI2S, (5) relinfectivity + relsusceptibility, (6) alle drie voorgaande
# in geval (2) en (3) moet een probI2S worden gekozen. Default 50%.
aantallensymptomatisch <- dataOsiris %>%
  filter(ZIE1eZiekteDt < as.Date("2020-03-20")) %>%
  select(Leeftijdsgroep) %>%
  separate(Leeftijdsgroep, c("onder", "boven")) %>% 
  filter(onder != "Niet") %>%
  mutate(onder = as.numeric(onder),
         age_class = floor(onder / 10),
         age_class = if_else(age_class > 8, 8, age_class)) %>%
  group_by(age_class) %>%
  summarise(aantal = n()) %>% pull()
verhoudingensymptomatisch <- aantallensymptomatisch/sum(aantallensymptomatisch)


#######
# (1) #
#######
infectieverhoudingen <- abs(eigen(ctrlmats$`99` * agedist())$vectors[, 1])
vectorfit <- verhoudingensymptomatisch/infectieverhoudingen
vectorfit <- vectorfit/max(vectorfit)
relincidence <- abs(eigen(ctrlmats$`99` * agedist())$vectors[, 1])
# relsymp <- relincidence * vectorfit
# relsymp <- relsymp/sum(relsymp)
# verhoudingensymptomatisch/relsymp

AgeNormList <- list(
  Symp = vectorfit
)
MaxSymptFractionList <- list(
  Symp = sum(relincidence * vectorfit)/sum(relincidence)
)

#######
# (2) #
#######
agenormfit <- function(agenormvector) {
  ev <- abs(eigen(ctrlmats$`99` * agedist() * c(exp(agenormvector), 1))$vectors[, 1])
  ev <- ev/(sum(ev))
  sum(log(ev / verhoudingensymptomatisch)^2)
}
vectorfit <- c(exp(optim(rep(0,8), agenormfit, method = "CG")$par), 1)
vectorfit <- vectorfit/(max(vectorfit))
relincidence <- abs(eigen(ctrlmats$`99` * agedist() * vectorfit)$vectors[, 1])
# relsymp <- relincidence * 1
# relsymp <- relsymp/sum(relsymp)
# verhoudingensymptomatisch/relsymp

AgeNormList <- c(AgeNormList,
                 list(
                   Susc = vectorfit
                 ))
MaxSymptFractionList <- c(MaxSymptFractionList,
                          list(
                            Susc = 1
                          ))

#######
# (3) #
#######
agenormfit <- function(agenormvector) {
  ev <- abs(eigen(ctrlmats$`99` * agedist() * c(exp(agenormvector), 1))$vectors[, 1])
  ev <- ev/(sum(ev))
  sv <- verhoudingensymptomatisch/c(exp(agenormvector), 1)
  sv <- sv/(sum(sv))
  sum(log(ev / sv)^2)
}
vectorfit <- c(exp(optim(rep(0,8), agenormfit, method = "CG")$par), 1)
vectorfit <- vectorfit/(max(vectorfit))
relincidence <- abs(eigen(ctrlmats$`99` * agedist() * vectorfit)$vectors[, 1])
# relsymp <- relincidence * vectorfit
# relsymp <- relsymp/sum(relsymp)
# verhoudingensymptomatisch/relsymp

AgeNormList <- c(AgeNormList,
                 list(
                   SympSusc = vectorfit
                 ))
MaxSymptFractionList <- c(MaxSymptFractionList,
                          list(
                            SympSusc = sum(relincidence * vectorfit)/sum(relincidence)
                          ))

#######
# (4) #
#######
agenormfit <- function(agenormvector) {
  ev <- abs(eigen(t(ctrlmats$`99` * c(exp(agenormvector), 1)) * agedist())$vectors[, 1])
  ev <- ev/(sum(ev))
  sv <- verhoudingensymptomatisch/c(exp(agenormvector), 1)
  sv <- sv/(sum(sv))
  sum(log(ev / sv)^2)
}
vectorfit <- c(exp(optim(rep(0,8), agenormfit, method = "CG")$par), 1)
vectorfit <- vectorfit/(max(vectorfit))
relincidence <- abs(eigen(t(ctrlmats$`99` * vectorfit) * agedist())$vectors[, 1])
# relsymp <- relincidence * vectorfit
# relsymp <- relsymp/sum(relsymp)
# verhoudingensymptomatisch/relsymp

AgeNormList <- c(AgeNormList,
                 list(
                   SympInf = vectorfit
                 ))
MaxSymptFractionList <- c(MaxSymptFractionList,
                          list(
                            SympInf = sum(relincidence * vectorfit)/sum(relincidence)
                          ))


#######
# (5) #
#######
agenormfit <- function(agenormvector) {
  ev <- abs(eigen(t(ctrlmats$`99` * c(exp(agenormvector), 1)) * agedist() * c(exp(agenormvector), 1))$vectors[, 1])
  ev <- ev/(sum(ev))
  sv <- verhoudingensymptomatisch
  sv <- sv/(sum(sv))
  sum(log(ev / sv)^2)
}
vectorfit <- c(exp(optim(rep(0,8), agenormfit, method = "CG")$par), 1)
vectorfit <- vectorfit/(max(vectorfit))
relincidence <- abs(eigen(t(ctrlmats$`99` * vectorfit) * agedist() * vectorfit)$vectors[, 1])
# relsymp <- relincidence
# relsymp <- relsymp/sum(relsymp)
# verhoudingensymptomatisch/relsymp

AgeNormList <- c(AgeNormList,
                 list(
                   SuscInf = vectorfit
                 ))
MaxSymptFractionList <- c(MaxSymptFractionList,
                          list(
                            SuscInf = 1
                          ))


#######
# (6) #
#######
agenormfit <- function(agenormvector) {
  ev <- abs(eigen(t(ctrlmats$`99` * c(exp(agenormvector), 1)) * agedist() * c(exp(agenormvector), 1))$vectors[, 1])
  ev <- ev/(sum(ev))
  sv <- verhoudingensymptomatisch/c(exp(agenormvector), 1)
  sv <- sv/(sum(sv))
  sum(log(ev / sv)^2)
}
vectorfit <- c(exp(optim(rep(0,8), agenormfit, method = "CG")$par), 1)
vectorfit <- vectorfit/(max(vectorfit))
relincidence <- abs(eigen(t(ctrlmats$`99` * vectorfit) * agedist() * vectorfit)$vectors[, 1])
# relsymp <- relincidence * vectorfit
# relsymp <- relsymp/sum(relsymp)
# verhoudingensymptomatisch/relsymp

AgeNormList <- c(AgeNormList,
                 list(
                   SympSuscInf = vectorfit
                 ))
MaxSymptFractionList <- c(MaxSymptFractionList,
                          list(
                            SympSuscInf = sum(relincidence * vectorfit)/sum(relincidence)
                          ))



#############################
### Infection to symptoms ###
#############################
probI2S <- function(normalisedby = "Symp", meanprI2S = 0.5, ...) {
  toreturn <- if(normalisedby %in% c("Susc", "SuscInf")) 
    meanprI2S else AgeNormList[[normalisedby]] * meanprI2S / MaxSymptFractionList[[normalisedby]]
  if(max(toreturn) > 1) {
    toreturn <- toreturn/max(toreturn)
  }
  return(toreturn)
}

###############################
### Relative Susceptibility ###
###############################
relsusceptibility <- function(normalisedby = "Symp", ...) {
  toreturn <-  if(normalisedby %in% c("Symp", "SympInf"))
    1 else AgeNormList[[normalisedby]] 
  return(toreturn)
}


###############################
### Relative Infectiousness ###
###############################
# reduction due to control; "99" is < 13 March
ctrlinfs <- list(
  `00` = 1,
  `99` = 0.946,
  IsoMild = 0.98,  
  HQMild = 0.86
)
# age-dependent infectivity
relinfectivity <- function(normalisedby = "Symp", scenario = "00", nr = 0, relinfsamples = c(), ...) {
  toreturn <-  if(normalisedby %in% c("Symp", "Susc", "SympSusc")) 
    1 else AgeNormList[[normalisedby]] 
  toreturn <- toreturn * ctrlinfs[[scenario]]
  if(nr == 0) {
    rndcorrection <- 1
  } else {
    rndcorrection <- relinfsamples[[scenario]][nr]
  } 
  return(toreturn * rndcorrection)
}


# natural history parameters
parmsNatHist <- function(Rstart = 2.2, SI = 3.5, relinfstart = "99", ...) {
  return(list(betaA = 0,
              betaB = 0,
              betaC = (Rstart / ctrlinfs[[relinfstart]]) * 3.5 / SI / 2,
              betaD = (Rstart / ctrlinfs[[relinfstart]]) * 3.5 / SI / 2,
              gammaA = 3.5 / SI,
              gammaB = 3.5 / SI,
              gammaC = 3.5 / SI,
              gammaD = 3.5 / SI))
}

rm(aantallensymptomatisch, verhoudingensymptomatisch,
   infectieverhoudingen, agenormfit,
   vectorfit, relincidence)