#######################################################
### contacts, infectivities,                        ###
###   susceptibilities, symptomatic fractions,      ###
###   natural history and control                   ###
###     v2: v1 + relinfectivity as contactmatrix    ###
###     v3: serology_pienter directly               ###
###                 to hospitalisations NICE        ###
###           => OSIRIS only for symptomatic        ###
#######################################################

### 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[["99"]])
  } else {
    return(ctrlmats[[nr]])
  }
} 

### from contact matrix to case age distribution in hospital
# normaliseren in twee stappen: contact matrix naar serologie, via (1) relsusceptibility of 
# (2) relinfectivity + relsusceptibility (default); vervolgens serologie naar ziekenhuis in NICE
# of symptomatisch in OSIRIS.
# Aanname 1: PIENTERsera op één (gemiddelde) dag afgenomen, reflecteren OSIRIS van 14 en NICE van 7 dagen eerder
#             (5 dagen incubatie, 7 dagen symp-to-hosp, 19 dagen infectie-seropositief)
# Aanname 2: leeftijdsverdelinge p(sero->symptoom) = constante * p(symptoom->OSIRIS), dus evenredig aan sqrt(OSIRIS/sero)
#             (hiermee krijg je ongeveer 50% symptomatic rate (ietsje lager door aftoppen in 80+))
########################################################################################
serodata_pico1 <- read_csv("data/pico1____UITZONDERINGSGROND_2___.csv")

meansampleday <- serodata_pico1 %>%
  filter(region < 6) %>%
  mutate(day = as.Date(algdatuminvul, "%d%b%Y")) %>% 
  pull(day) %>% mean(na.rm = T) %>% round()

aantallensymptomatisch <- 
  dataOsiris %>%
  filter(ZIE1eZiekteDt <= meansampleday - 14) %>%
  select(Leeftijdsgroep) %>%
  separate(Leeftijdsgroep, c("onder", "boven")) %>% 
  mutate(onder = if_else(onder == "Niet", "-1", onder)) %>%
  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(totaalS2Hp = n()) %>% 
  pull(totaalS2Hp)

aantallenhospitalisaties <-
  pddata4fit %>%
  filter(ti <= as.numeric(meansampleday - 7 - as.Date("2020-02-12")) & ti > 0) %>%
  mutate(weight = mapply(function(x,y) 1/NICEprobabilities$probSe2A[x, y], x = ti, y = age + 1)) %>%
  group_by(age) %>%
  summarise(aantal.weighted = sum(weight),
            aantal.unweighted = n()) %>% pull(aantal.weighted)

seroverhoudingen <- serodata_pico1 %>%
  filter(region < 6) %>%
  mutate(lftgroep = floor(leeftijd/10),
         lftgroep = pmin(lftgroep, 8)) %>%
  group_by(lftgroep) %>%
  summarise(seropos = sum(positive) / n(),
            stderror = sqrt(seropos * (1 - seropos) / n()),
            lower = seropos - 1.96 * stderror,
            upper = seropos + 1.96 * stderror) %>%
  ungroup() %>%
  mutate(popnrs = agedist() * popsize(),
         infnrs = seropos * popnrs,
         # verhouding in geinfecteerden: eigenvector van transmissiematrix
         infprop = infnrs/sum(infnrs),
         sympnrs = aantallensymptomatisch,
         probsymp = sympnrs/infnrs,
         hospnrs = aantallenhospitalisaties,
         probhosp = hospnrs/infnrs)


AgeNormList <-
  list(
    I2S = sqrt(seroverhoudingen$probsymp),
    I2Se = seroverhoudingen$probhosp
  )
MeanSymptFractionList <- list(
  I2S = sum(seroverhoudingen$infprop * AgeNormList$I2S)/sum(seroverhoudingen$infprop),
  I2Se = sum(seroverhoudingen$infprop * AgeNormList$I2Se)/sum(seroverhoudingen$infprop)
)

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

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

#######
# (2) #
#######
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 <- seroverhoudingen$infprop
  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))

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




#############################
### Infection to symptoms ###
#############################
probI2S <- function(meanprI2S = 0.5, ...) {
  toreturn <- pmin(1, AgeNormList$I2S * meanprI2S / MeanSymptFractionList$I2S)
  return(toreturn)
}

############################################
### Infection to severe (hospitalisable) ###
############################################
probI2Se <- function() {
  return(AgeNormList$I2Se)
}

###############################
### Relative Susceptibility ###
###############################
relsusceptibility <- function(normalisedby = "SuscInf", ...) {
  toreturn <-  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 = "SuscInf", scenario = "00", nr = 0, ...) {
  toreturn <-  if(normalisedby == "Susc") 
    1 else AgeNormList[[normalisedby]] 
  if(scenario %in% names(ctrlinfs)) {
    return(toreturn * ctrlinfs[[scenario]])
  } else if(nr == 0) {
    return(toreturn * ctrlinfs[[`99`]])
  } else {
    return(toreturn * ctrlinfs[[nr]])
  }
}


# 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, 
   aantallenhospitalisaties, 
    agenormfit,
   vectorfit)