###########################################################
### Analyse voor ___UITZONDERINGSGROND_2___, data tot 24 september   ###
###   SEEIIR-model, SI = 4 dagen                        ###
### Met nieuwe NICE-analyses (veranderende kansen)      ###
###   en PIENTER-serologie,                             ###
###   Leeftijdsmodel gescheiden in deel tot infecties   ###
###     via vatbaarheid/infectiviteit (SuscInf) en deel ###
###     tot symptomen/ziekenhuisopname. Scheiding       ###
###     mogelijk door serologie                         ###
### Toegevoegd: regionale ROAZ-prognoses                ###
### EXTRA op 28 sept: R = 0.9 vanaf wo 30 sept, m
###########################################################

library(tidyverse)
library(lubridate)
library(deSolve)

################################
### eerst alle data-analyses ###
################################
analysisdate <- as.Date("2020-09-24")
source("R/code4ode/Reportingdelays4ode.R")
source("R/code4ode/NICEanalyses4ode_v6.R")
source("R/code4ode/OSIRISanalyses4ode_v2.R")



################################
### INPUT voor de simulaties ###
################################

### population
source("R/code4ode/populationdata4ode.R")

### infection to symptoms, contacts, relative infectivities, natural history and control
source("R/code4ode/ContactsInfectivities4ode_v3.R")

# possibility to add alternative matrices and infectivities
ctrlmatslockdown <- readRDS("data/Contactmatrix_D3asEpiPose1_residualincreased_27mei2020.rds")
ctrlmatsbatch1 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch1_29mei2020.rds")
ctrlmatsbatch2 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch2_29mei2020.rds")
ctrlmatsbatch3 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3_25juni2020.rds")
ctrlmatssep2020 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3vanaf1sep_6aug2020.rds")
ctrlmats <- c(readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch1_29mei2020.rds"),             
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch2_29mei2020.rds"),
              readRDS("data/ContactmatricesD3praktijk_midpoint_24mrt2020.rds"),
              list(lockdown = Reduce(`+`, ctrlmatslockdown)/200,
                   batch1 = Reduce(`+`, ctrlmatsbatch1)/200,
                   batch2 = Reduce(`+`, ctrlmatsbatch2)/200,
              batch3 = Reduce(`+`, ctrlmatsbatch3)/200,
              batchsep = Reduce(`+`, ctrlmatssep2020)/200)
)
ctrlinfs <- list(
  `00` = 1,
  `99` = 0.946,
  IsoMild = 0.98,  
  HQMild = 0.86
)
set.seed(filedate)
devs <- runif(200, -1, 1)
ctrlinfs <- c(as.list(
  rep(0.98 * (1 + 0.02 * devs), 6)),
  list(
    `00` = 1,
    `99` = 0.946,
    IsoMild = 0.98,  
    HQMild = 0.86
  ))
rm(devs)


### downstream numbers: symptomatic, hospitalised, ICU, mortality
source("R/code4ode/DelaysProbabilities4ode_v2.R")

### model code
source("R/code4ode/modelcode4ode_v4.R")

### reporting delays
repdelayICinc <- get_reportingdelay_NICE(end_date = analysisdate)
disdelayICinc <- get_reportingdelay_NICE(discharge = T, end_date = analysisdate)
repdelayhospinc <- get_reportingdelay_NICE(IC = F, end_date = analysisdate)
disdelayhospinc <- get_reportingdelay_NICE(IC = F, discharge = T, end_date = analysisdate)
observedICinc <- inc_nice_ic()

### Calibration with control up to "today"
# parameters to change: seeding, fittedrelinfs, probdeath
# endtimes: 32 = Sun 15 March, 44 = Fri 27 March (based on better fit than 23 March)
logLikijk <- function(parms) {
  siminc <- covidcontrolsimple(contactcontrol = c("99", "lockdown", "lockdown",  "batch1", "batch2", "batch3", "batchsep"),
                               infectivitycontrol = c("99", "IsoMild", "IsoMild",  "IsoMild", "IsoMild", "IsoMild", "IsoMild"),
                               endtimes = c(32, 50, 
                                            as.numeric(as.Date("2020-05-10") - as.Date("2020-02-12")), 
                                            as.numeric(as.Date("2020-06-01") - as.Date("2020-02-12")),
                                            as.numeric(as.Date("2020-07-05") - as.Date("2020-02-12")),
                                            # as.numeric(filedate - 31 - as.Date("2020-02-12")),
                                            as.numeric(as.Date("2020-08-30") - as.Date("2020-02-12")),
                                            as.numeric(filedate - as.Date("2020-02-12"))),
                               seeding = exp(parms[1]), SI = 4,
                               fittedrelinfs = exp(parms[c(2,3,4,4,4,4,5)])) %>%
    filter(time > 0 & time <= as.numeric(filedate - as.Date("2020-02-12"))) %>%
    pull(totincICU) 
  # return(siminc)
  siminc[seq(to = length(siminc), length.out = 100)] <- 
    siminc[seq(to = length(siminc), length.out = 100)] * rev(repdelayICinc)
  - sum(dpois(observedICinc, siminc, log = T))
}
logLikijk_1 <- function(parms) {
  siminc <- covidcontrolsimple(contactcontrol = c("99", "lockdown", "lockdown",  "batch1", "batch2", "batch3", "batchsep"),
                               infectivitycontrol = c("99", "IsoMild", "IsoMild",  "IsoMild", "IsoMild", "IsoMild", "IsoMild"),
                               endtimes = c(32, 50, 
                                            as.numeric(as.Date("2020-05-10") - as.Date("2020-02-12")), 
                                            as.numeric(as.Date("2020-06-01") - as.Date("2020-02-12")),
                                            as.numeric(as.Date("2020-07-05") - as.Date("2020-02-12")),
                                            # as.numeric(filedate - 31 - as.Date("2020-02-12")),
                                            as.numeric(as.Date("2020-08-30") - as.Date("2020-02-12")),
                                            as.numeric(filedate - as.Date("2020-02-12"))),
                               seeding = exp(5.04), SI = 4,
                               fittedrelinfs = exp(parms[c(1,2,3,3,3,3,3)])) %>%
    filter(time > 0 & time <= as.numeric(filedate - as.Date("2020-02-12"))) %>%
    pull(totincICU)
  siminc[seq(to = length(siminc), length.out = 100)] <- 
    siminc[seq(to = length(siminc), length.out = 100)] * rev(repdelayICinc)
  - sum(dpois(observedICinc, siminc, log = T))
}
optim(c(4.79,.083,-.342,-.110), logLikijk) # fri = 2344444 (32,50), minloglik = 477.36, ests = 4.88; relinf1=.070; relinf2 = -.296; relinf3=-.110
optim(c(4.88,.069,-.293,-.113,-.113), logLikijk) # fri = 2344445 (32,50), minloglik = 469.91, ests = 4.92; relinf1=.065; relinf2 = -.285; relinf3=-.117; relinf4=0.086
optim(c(4.88,.069,-.293,-.113,-.113), logLikijk) # fri = 2344455 (32,50), minloglik = 477.22, ests = 4.93; relinf1=.065; relinf2 = -.291; relinf3=-.113; relinf4=-.105
optim(c(5.00,.057,-.286,-.116,-.116,-.116), logLikijk) # fri = 2344566 (32,50), minloglik = 447.14, ests = 4.98; relinf1=.061; relinf2 = -.307; relinf3=-.095; relinf4=-.128; relinf5=-.131


# 1000 samples met beste model: fri = 2344445 (32,50).
ijkres <- optim(c(4.88,.069,-.293,-.113,-.113), 
                logLikijk, hessian = T) # IsoMild, IsoMild, IsoMild
# ijkres <- optim(c(.083,-.342,-.110), 
#                 logLikijk_1, hessian = T) # IsoMild, IsoMild, IsoMild
library(mvtnorm)
set.seed(filedate)
rndpars <- rmvnorm(200, ijkres$par, solve(ijkres$hessian))
# rndpars <- cbind(rep(5.04, 200), rndpars)


### relinfectivity in de contactmatrices integreren
fittedrelinfmatrix <- function(fri_ctrl, fri_base, ageprops, method = 1) {
  frivector <- (fri_ctrl / fri_base) * ageprops + 1 - ageprops
  frimatrix <- matrix(rep(frivector, 9), nrow = 9)
  if(method == 1) {
    ### both age classes contribute
    frimatrix <- sqrt(frimatrix) * sqrt(t(frimatrix))
  } else if(method == 2) {
    ### active age class dominates
    frimatrix <- pmax(frimatrix, t(frimatrix))
  } else {
    ### inactive age class dominates
    frimatrix <- pmin(frimatrix, t(frimatrix))
  }
  return(frimatrix * fri_base)
}



## nieuwe contactmatrices maken,
ctrlmatslockdown <- readRDS("data/Contactmatrix_D3asEpiPose1_residualincreased_27mei2020.rds")
ctrlmatsbatch1 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch1_29mei2020.rds")
ctrlmatsbatch2 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch2_29mei2020.rds")
ctrlmatsbatch3 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3_25juni2020.rds")
ctrlmatssep2020 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3vanaf1sep_6aug2020.rds")
ctrlmats <- c(readRDS("data/Contactmatrix_D3asEpiPose1_residualincreased_27mei2020.rds"),
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch1_29mei2020.rds"),             
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch2_29mei2020.rds"),
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3_25juni2020.rds"),
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3vanaf1sep_6aug2020.rds"),
              readRDS("data/ContactmatricesD3praktijk_midpoint_24mrt2020.rds"),
              list(lockdown = Reduce(`+`, ctrlmatslockdown)/200,
                   batch1 = Reduce(`+`, ctrlmatsbatch1)/200,
                   batch2 = Reduce(`+`, ctrlmatsbatch2)/200,
                   batch3 = Reduce(`+`, ctrlmatsbatch3)/200,
                   batchsep2020 = Reduce(`+`, ctrlmatssep2020)/200)
)
for(i in 1:200) {
  for(j in 0) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 4]), exp(rndpars[i, 2]), rep(1, 9))
  }
  for(j in 1) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 4]), exp(rndpars[i, 2]), rep(1, 9))
  }
  for(j in 2:3) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 4]), exp(rndpars[i, 2]), rep(1, 9))
  }
  for(j in 4) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 5]), exp(rndpars[i, 2]), rep(1, 9))
  }
}

relbeta <- 1/eigen(t(contactmatrix() * 
                       agedist() * relsusceptibility()) * 
                     relinfectivity(), only.values = T)[[1]][1]
eigen(t(ctrlmats[["99"]] * agedist() * relsusceptibility()) * relinfectivity(scenario =  "99") * exp(ijkres$par[2]) *
        2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta
eigen(t(ctrlmats[["lockdown"]] * agedist() * relsusceptibility()) * relinfectivity(scenario = "IsoMild") * exp(ijkres$par[3]) *
        2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta
eigen(t(ctrlmats[["lockdown"]] * agedist() * relsusceptibility()) * relinfectivity(scenario = "IsoMild") * exp(ijkres$par[4]) *
        2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta
mean(sapply(1:200, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                         2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(201:400, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                         2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(401:600, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                         2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(601:800, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                         2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(801:1000, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                          2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))



panels_rijen <- c(NA)
panels_kolommen <- c(NA)
panels_kleuren <- c("geen bestrijding", "huidige pakket")#, "lockdown", "thuisquarantaine")

contactkeus <- function(rij = 1, kolom = 1, kleur = 1) {
  case_when(
    kleur == 1 ~ c("99", "99", "99", "99", "99", "99", "99", "99"),
    kleur == 2 ~ c("99", "lockdown", "lockdown", "batch1", "batch2", "batch3", "batchsep2020", "RND")
  )
}

infectiviteitkeus <- function(rij = 1, kolom = 1, kleur = 1) {
  case_when(
    kleur == 1 ~ c("99", "99", "99", "99", "99", "99", "99", "99"),
    kleur == 2 ~ c("99", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "RND")
  )
}

eindtijden <- function() c(32, 50, 
                           as.numeric(as.Date("2020-05-10") - as.Date("2020-02-12")),
                           as.numeric(as.Date("2020-06-01") - as.Date("2020-02-12")),
                           as.numeric(as.Date("2020-07-05") - as.Date("2020-02-12")),
                           as.numeric(as.Date("2020-08-30") - as.Date("2020-02-12")),
                           as.numeric(filedate - as.Date("2020-02-12")) - 14, 
                           400)

allfigdata <- list()
for(repnr in 1:200) {
  allsimulresults <- list()
  for(i in 1:1) {
    for(j in 1:1) {
      for(k in 1:2) {
        simul <- covidcontrolsim(contactcontrol = contactkeus(kleur = k), 
                                 infectivitycontrol = infectiviteitkeus(kleur = k),
                                 endtimes = eindtijden(),  
                                 seeding = exp(rndpars[repnr, 1]),  
                                 Rstart = 2.2, 
                                 SI = 4, 
                                 fittedrelinfs = list(exp(rndpars[repnr, c(2,2,2,2,2,2,2,2)]),
                                                      c(exp(rndpars[repnr, c(2,3,4,4,4,4,4)]), 1))[[k]], 
                                 nr = c(repnr, repnr, repnr, repnr, repnr, repnr, repnr, repnr + 800))
        allsimulresults <- c(allsimulresults,
                             list(simul %>%
                                    mutate(rij = panels_rijen[i],
                                           kolom = panels_kolommen[j],
                                           kleur = panels_kleuren[k]))
        )
      }
    }
  }
  allfigdata <- c(allfigdata,
                  list(do.call(rbind, allsimulresults) %>%
                          mutate(tijd = as.Date("2020-02-12") + time) 
                  ))
  cat(repnr)
}

saveRDS(allfigdata, "results/allfigdata___UITZONDERINGSGROND_2____25092020")
allfigdata <- readRDS("results/allfigdata___UITZONDERINGSGROND_2____25092020")

allfigdata <- lapply(1:200,
                     function(x) allfigdata[[x]] %>% 
                       mutate(repl = x))
figuredata <- do.call(rbind, allfigdata) 

################
### 28 SEPT: VINDEN VAN RELATIEVE INFECTIVITEIT VOOR RT = 0.9
###################
immunitybyage <- figuredata %>%
  filter(tijd == as.Date("2020-04-29") & kleur == "huidige pakket") %>%
  select(ageclass, repl, immune) %>%
  pivot_wider(names_from = repl, values_from = immune) %>%
  select(-ageclass) %>%
  as.matrix()
currentRt <- mean(sapply(801:1000, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility() * (1 - immunitybyage[,x-800])) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                          2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
requiredreduction <- 0.9/currentRt


##########
### NEW SIMULATIONS (28 SEPT)
###########
## nieuwe contactmatrices maken,
ctrlmatslockdown <- readRDS("data/Contactmatrix_D3asEpiPose1_residualincreased_27mei2020.rds")
ctrlmatsbatch1 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch1_29mei2020.rds")
ctrlmatsbatch2 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch2_29mei2020.rds")
ctrlmatsbatch3 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3_25juni2020.rds")
ctrlmatssep2020 <- readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3vanaf1sep_6aug2020.rds")
ctrlmats <- c(readRDS("data/Contactmatrix_D3asEpiPose1_residualincreased_27mei2020.rds"),
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch1_29mei2020.rds"),             
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch2_29mei2020.rds"),
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3_25juni2020.rds"),
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3vanaf1sep_6aug2020.rds"),
              readRDS("data/Contactmatrix_D3asEpiPose1_residualplus_batch3vanaf1sep_6aug2020.rds"),
              readRDS("data/ContactmatricesD3praktijk_midpoint_24mrt2020.rds"),
              list(lockdown = Reduce(`+`, ctrlmatslockdown)/200,
                   batch1 = Reduce(`+`, ctrlmatsbatch1)/200,
                   batch2 = Reduce(`+`, ctrlmatsbatch2)/200,
                   batch3 = Reduce(`+`, ctrlmatsbatch3)/200,
                   batchsep2020 = Reduce(`+`, ctrlmatssep2020)/200)
)
for(i in 1:200) {
  for(j in 0) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 4]), exp(rndpars[i, 2]), rep(1, 9))
  }
  for(j in 1) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 4]), exp(rndpars[i, 2]), rep(1, 9))
  }
  for(j in 2:3) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 4]), exp(rndpars[i, 2]), rep(1, 9))
  }
  for(j in 4) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 5]), exp(rndpars[i, 2]), rep(1, 9))
  }
  for(j in 5) {
    ctrlmats[[i + 200 * j]] <- requiredreduction *
      ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 5]), exp(rndpars[i, 2]), rep(1, 9))
  }
}

relbeta <- 1/eigen(t(contactmatrix() * 
                       agedist() * relsusceptibility()) * 
                     relinfectivity(), only.values = T)[[1]][1]
eigen(t(ctrlmats[["99"]] * agedist() * relsusceptibility()) * relinfectivity(scenario =  "99") * exp(ijkres$par[2]) *
        2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta
eigen(t(ctrlmats[["lockdown"]] * agedist() * relsusceptibility()) * relinfectivity(scenario = "IsoMild") * exp(ijkres$par[3]) *
        2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta
eigen(t(ctrlmats[["lockdown"]] * agedist() * relsusceptibility()) * relinfectivity(scenario = "IsoMild") * exp(ijkres$par[4]) *
        2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta
mean(sapply(1:200, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                       2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(201:400, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                         2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(401:600, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                         2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(601:800, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                         2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(801:1000, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                          2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))
mean(sapply(1001:1200, function(x) eigen(t(ctrlmats[[x]] * agedist() * relsusceptibility()) * ctrlinfs[[x]] * relinfectivity(scenario = "IsoMild") * 
                                          2*parmsNatHist(SI = 4)$betaC/parmsNatHist(SI = 4)$gammaC, only.values = T)[[1]][1] * relbeta))


panels_rijen <- c(NA)
panels_kolommen <- c(NA)
panels_kleuren <- c("huidige pakket", "Rt = 0.9 vanaf 30 sept")#, "lockdown", "thuisquarantaine")

contactkeus <- function(rij = 1, kolom = 1, kleur = 1) {
  case_when(
    kleur == 1 ~ c("99", "lockdown", "lockdown", "batch1", "batch2", "batch3", "batchsep2020", "RND", "RND"),
    kleur == 2 ~ c("99", "lockdown", "lockdown", "batch1", "batch2", "batch3", "batchsep2020", "RND", "RND")
  )
}

infectiviteitkeus <- function(rij = 1, kolom = 1, kleur = 1) {
  case_when(
    kleur == 1 ~ c("99", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "RND", "RND"),
    kleur == 2 ~ c("99", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "IsoMild", "RND", "RND")
  )
}

eindtijden <- function() c(32, 50, 
                           as.numeric(as.Date("2020-05-10") - as.Date("2020-02-12")),
                           as.numeric(as.Date("2020-06-01") - as.Date("2020-02-12")),
                           as.numeric(as.Date("2020-07-05") - as.Date("2020-02-12")),
                           as.numeric(as.Date("2020-08-30") - as.Date("2020-02-12")),
                           as.numeric(filedate - as.Date("2020-02-12")) - 14, 
                           as.numeric(as.Date("2020-09-29") - as.Date("2020-02-12")),
                           400)

allfigdata <- list()
for(repnr in 1:200) {
  allsimulresults <- list()
  for(i in 1:1) {
    for(j in 1:1) {
      for(k in 1:2) {
        simul <- covidcontrolsim(contactcontrol = contactkeus(kleur = k), 
                                 infectivitycontrol = infectiviteitkeus(kleur = k),
                                 endtimes = eindtijden(),  
                                 seeding = exp(rndpars[repnr, 1]),  
                                 Rstart = 2.2, 
                                 SI = 4, 
                                 fittedrelinfs = c(exp(rndpars[repnr, c(2,3,4,4,4,4,4)]), 1, 1), 
                                 nr = c(repnr, repnr, repnr, repnr, repnr, repnr, repnr, repnr + 800, repnr + 600 + 200*k))
        allsimulresults <- c(allsimulresults,
                             list(simul %>%
                                    mutate(rij = panels_rijen[i],
                                           kolom = panels_kolommen[j],
                                           kleur = panels_kleuren[k]))
        )
      }
    }
  }
  allfigdata <- c(allfigdata,
                  list(do.call(rbind, allsimulresults) %>%
                         mutate(tijd = as.Date("2020-02-12") + time) 
                  ))
  cat(repnr)
}

saveRDS(allfigdata, "results/allfigdata___UITZONDERINGSGROND_2____28092020")
allfigdata <- readRDS("results/allfigdata___UITZONDERINGSGROND_2____25092020")

allfigdata <- lapply(1:200,
                     function(x) allfigdata[[x]] %>% 
                       mutate(repl = x))
figuredata <- do.call(rbind, allfigdata) 


###########################
### projection per ROAZ ###
###########################
ROAZprior <- sapply(1:11, function(x) tail(inc_nice_hosp(ROAZ = x), 15)) %>% colSums() %>% `+`(1)
set.seed(filedate)
ROAZdists <- MCMCpack::rdirichlet(200, ROAZprior)
colnames(ROAZdists) <- paste0("ROAZ_", substr(seq(.011,.111,.01), 3, 4))

figuredataROAZ <- figuredata %>%
  filter(kleur == "huidige pakket" & tijd > analysisdate - 15 & tijd < analysisdate + 15) %>%
  select(time, ageclass, inchosp, incICU, kleur, tijd, repl) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totICU = sum(incICU),
    tothosp = sum(inchosp)
  ) %>%
  ungroup() %>%
  full_join(as_tibble(ROAZdists) %>%
              mutate(repl = 1:200), by = "repl") %>%
  pivot_longer(ROAZ_01:ROAZ_11, names_to = "ROAZ", values_to = "weight") %>% 
  mutate(simhosp = tothosp * weight,
         simICU =  totICU * weight) 
figuredataROAZ <- figuredataROAZ %>%
  group_by(tijd, ROAZ, kleur) %>%
  summarise(
    med = round(median(simICU)),
    low = qpois(.025, quantile(simICU, .025)),
    upp = qpois(.975, quantile(simICU, .975)),
    location = "ICU"
  ) %>%
  ungroup() %>%
  bind_rows(figuredataROAZ %>%
              group_by(tijd, ROAZ, kleur) %>%
              summarise(
                med = round(median(simhosp)),
                low = qpois(.025, quantile(simhosp, .025)),
                upp = qpois(.975, quantile(simhosp, .975)),
                location = "hosp"
              ) %>%
              ungroup()) %>%
  mutate(
    ROAZ = sort(unique(ROAZdata$Regio_ROAZ))[-9][as.numeric(substr(ROAZ, 6, 7))]
  )
  
obsdata <- bind_cols(
  tibble(tijd = seq(as.Date("2020-02-13"), filedate, 1)),
  as_tibble(matrix(sapply(1:11, function(x) inc_nice_hosp(ROAZ = x)), ncol = 11, 
                   dimnames = list(NULL, paste0("ROAZ_", substr(seq(.011,.111,.01), 3, 4)))))) %>%
  pivot_longer(ROAZ_01:ROAZ_11, names_to = "ROAZ", values_to = "obs") %>%
  mutate(location = "hosp") %>%
  bind_rows(bind_cols(
    tibble(tijd = seq(as.Date("2020-02-13"), filedate, 1)),
    as_tibble(matrix(sapply(1:11, function(x) inc_nice_ic(ROAZ = x)), ncol = 11, 
                     dimnames = list(NULL, paste0("ROAZ_", substr(seq(.011,.111,.01), 3, 4)))))) %>%
      pivot_longer(ROAZ_01:ROAZ_11, names_to = "ROAZ", values_to = "obs") %>%
      mutate(location = "ICU")) %>%
  mutate(
    ROAZ = sort(unique(ROAZdata$Regio_ROAZ))[-9][as.numeric(substr(ROAZ, 6, 7))]
  ) %>%
  filter(tijd > analysisdate - 15 & tijd < analysisdate + 15)


# figuredataROAZ %>%
#   filter(location == "ICU") %>%
#   ggplot(aes(x = tijd, y = med)) +
#   geom_line(size = 2, color = "red") +
#   geom_ribbon(aes(ymin = low, ymax = upp, color = NULL), alpha = 0.3, fill = "red") +
#   geom_point(aes(y = obs), data = obsdata %>% filter(location == "ICU")) +
#   geom_text(aes(y = ytxt, label = dag), data = tibble(tijd = (-14:14) + analysisdate, 
#                                                       ytxt = 0,
#                                                       dag = c("z","m","d","w","d","v","z")[wday(analysisdate + (-14:14))]),
#             vjust = 1.3, size = 2.5) +
#   facet_wrap(~ ROAZ, ncol = 3) + 
#   scale_y_continuous(expand = c(.12,.1)) +
#   theme_light() +
#   labs(title = "Prognose IC-opnames per ROAZ-regio",
#        x = "Datum", y = "Aantal IC-opnames")
# ggsave("results/figROAZIC23092020.pdf", units = "mm", width = 200, height = 180, dpi = 600)

figuredataROAZ %>%
  filter(location == "ICU") %>%
  mutate(colorvar = "Modelprognose\n(mediaan en 95% interval)") %>%
  ggplot(aes(x = tijd, y = med)) +
  geom_crossbar(aes(ymin = low, ymax = upp, color = colorvar), fatten = 5, size = 0.2) +
  geom_point(aes(y = obs, shape = shapevar), data = obsdata %>% filter(location == "ICU") %>% mutate(shapevar = "NICE-data")) +
  geom_text(aes(y = ytxt, label = dag), data = tibble(tijd = (-14:14) + analysisdate, 
                                                      ytxt = 0,
                                                      dag = c("z","m","d","w","d","v","z")[wday(analysisdate + (-14:14))]),
            vjust = 1.3, size = 2.5) +
  facet_wrap(~ ROAZ, ncol = 3, scales = "free_y") + 
  scale_y_continuous(expand = c(.12,.1), breaks = function(x) seq(max(0, ceiling(x[1])),ceiling(x[2]),
                                                                  c(1,1,1,2,5,10,20,50,100,200,500)[1+floor(log(x[2] - x[1], 10^(1/3)))])) +
  theme_minimal() +
  labs(title = paste0("Prognose IC-opnames per ROAZ-regio (data t/m ", filedate,")"),
       x = "Datum", y = "Aantal IC-opnames") +
  scale_shape_discrete(name = NULL) +
  scale_color_manual(name = NULL, values = "black") +
  theme(panel.grid = element_blank(), axis.line.y = element_line(colour = "black"), 
        axis.ticks.y = element_line(colour = "black"), panel.grid.major.y = element_line(colour = "lightgrey"),
        legend.position = c(1, 0),
        legend.justification = c(1.2, 0)) 
ggsave(paste0("results/figROAZIC_", filedate, ".jpg"), units = "mm", width = 200, height = 180, dpi = 600)


figuredataROAZ %>%
  filter(location == "hosp") %>%
  mutate(colorvar = "Modelprognose\n(mediaan en 95% interval)") %>%
  ggplot(aes(x = tijd, y = med)) +
  geom_crossbar(aes(ymin = low, ymax = upp, color = colorvar), fatten = 5, size = 0.2) +
  geom_point(aes(y = obs, shape = shapevar), data = obsdata %>% filter(location == "hosp") %>% mutate(shapevar = "NICE-data")) +
  geom_text(aes(y = ytxt, label = dag), data = tibble(tijd = (-14:14) + analysisdate, 
                                                      ytxt = 0,
                                                      dag = c("z","m","d","w","d","v","z")[wday(analysisdate + (-14:14))]),
            vjust = 1.3, size = 2.5) +
  facet_wrap(~ ROAZ, ncol = 3, scales = "free_y") + 
  scale_y_continuous(expand = c(.12,.1), breaks = function(x) seq(max(0, ceiling(x[1])),ceiling(x[2]),
                                                                  c(1,1,1,2,5,10,20,50,100,200,500)[1+floor(log(x[2] - x[1], 10^(1/3)))])) +
  theme_minimal() +
  theme(panel.grid = element_blank(), axis.line.y = element_line(colour = "black"), 
        axis.ticks.y = element_line(colour = "black"), panel.grid.major.y = element_line(colour = "lightgrey")) +
  labs(title = paste0("Prognose ziekenhuisopnames per ROAZ-regio (data t/m ", filedate,")"),
       x = "Datum", y = "Aantal ziekenhuisopnames") +
  scale_shape_discrete(name = NULL) +
  scale_color_manual(name = NULL, values = "black") +
  theme(panel.grid = element_blank(), axis.line.y = element_line(colour = "black"), 
        axis.ticks.y = element_line(colour = "black"), panel.grid.major.y = element_line(colour = "lightgrey"),
        legend.position = c(1, 0),
        legend.justification = c(1.2, 0)) 
ggsave(paste0("results/figROAZhosp_", filedate, ".jpg"), units = "mm", width = 200, height = 180, dpi = 600)

figuredataROAZ %>%
  filter(location == "ICU") %>%
  mutate(predint = paste0(med, " (", low, " ; ", upp, ")")) %>%
  select(tijd, ROAZ, predint) %>%
  pivot_wider(names_from = ROAZ, values_from = predint) %>%
  rename(opnamedatum = tijd) %>%
  write_csv2(paste0("results/ICtabel_", filedate, ".csv"))
figuredataROAZ %>%
  filter(location == "hosp") %>%
  mutate(predint = paste0(med, " (", low, " ; ", upp, ")")) %>%
  select(tijd, ROAZ, predint) %>%
  pivot_wider(names_from = ROAZ, values_from = predint) %>%
  rename(opnamedatum = tijd) %>%
  write_csv2(paste0("results/hosptabel_", filedate, ".csv"))


%>%
  knitr::kable(format = "html", caption = "Prognose IC-opnames (mediaan en 95% interval)") %>%
  kableExtra::save_kable("results/ICtabel.pdf")
  

###############
### figuren ### 
############### 

### IC-incidence
obsdata <- tibble(
  tijd = seq(as.Date("2020-02-13"), filedate, 1),
  ICinc = inc_nice_ic()
) %>% bind_cols(
as_tibble(inc_nice_ic(ages=T))) %>%
            pivot_longer(3:11, names_to = "ageclass", values_to = "observed")  %>%
  group_by(ageclass) %>%
  mutate(ICcuminc = cumsum(ICinc),
         cumobserved = cumsum(observed)) %>%
  ungroup()

figuredata %>% 
  filter(time >= 25 & time < 350) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "geen verdere\nmaatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(incICU)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025)),
    upperbound = qpois(0.975, quantile(totsimulated, .975)),
    midbound = median(totsimulated)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICinc) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-opnames") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = ICinc, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijks aantal IC-opnames") +
  labs(title = "Aantal IC-opnames per dag") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,150)) 
ggsave(paste0("results/sims_ICincallshort_ribbon_", today(), ".jpg"), units = "mm", width = 240, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 110 & time < 350) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "geen verdere\nmaatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(incICU)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025)),
    upperbound = qpois(0.975, quantile(totsimulated, .975)),
    midbound = median(totsimulated)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICinc) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-opnames") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = ICinc, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijks aantal IC-opnames") +
  labs(title = "Aantal IC-opnames per dag") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,75)) 
ggsave(paste0("results/sims_ICincallnopeak_ribbon_", today(), ".jpg"), units = "mm", width = 200, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 25 & time < 155) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "met gevolgde\n maatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(incICU_C)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025)),
    upperbound = qpois(0.975, quantile(totsimulated, .975)),
    midbound = median(totsimulated)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICinc) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-opnames") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = ICinc, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijks aantal IC-opnames") +
  labs(title = "Aantal IC-opnames per dag (hypothetisch: bij gelijkgebleven opnamebeleid)") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,200)) 
ggsave(paste0("results/sims_ICincallshort_counterfact_ribbon_", today(), ".jpg"), units = "mm", width = 200, height = 120, dpi = 300)


figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "met gevolgde\n maatregelen",
                         kleur)) %>%
  group_by(kleur, ageclass, repl) %>%
  mutate(
    totsimulated = cumsum(incICU),
    totsimulated = if_else(kleur == "geen bestrijding", NA_real_, totsimulated)
  ) %>%
  group_by(tijd, kleur, ageclass) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "NICE IC-opnames") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 1) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = cumobserved, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Cumulatief aantal IC-opnames") +
  facet_wrap(~ ageclass, scales = "free_y") +
  labs(title = "Cumulatief aantal IC-opnames per leeftijdsgroep")
# ggsave(paste0("results/sims_ICinclft_", today(), ".jpg"), units = "mm", width = 180, height = 120, dpi = 300)

# IC-prevalence
obsdata <- tibble(
  tijd = seq(as.Date("2020-02-13"), filedate, 1),
  ICprev = prev_nice_ic()
) %>% bind_cols(
  as_tibble(prev_nice_ic(ages=T))) %>%
  pivot_longer(3:11, names_to = "ageclass", values_to = "observed")  

figuredata %>% 
  filter(time >= 25 & time < 350) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "geen verdere\nmaatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(prevICU)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICprev) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-bezetting") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = ICprev, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "IC-bezetting") +
  labs(title = "Aantal bezette IC-bedden") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,1500))  +
  geom_hline(aes(yintercept = 1208), color = "grey30", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-09-15"), y = 1280, 
           label = "Gewoonlijke Nederlandse IC-capaciteit", hjust = 1, size = 3, color = "grey30") 
# +
#   geom_hline(aes(yintercept = 2400), color = "black", linetype = "dashed") +
#   annotate(geom = "text", x = as.Date("2020-07-01"), y = 2470, 
#            label = "Opgeschaalde Nederlandse IC-capaciteit", hjust = 1, size = 3)
ggsave(paste0("results/sims_ICprevallshort_ribbon_", today(), ".jpg"), units = "mm", width = 240, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 110 & time < 350) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "geen verdere\nmaatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(prevICU)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICprev) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-bezetting") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = ICprev, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "IC-bezetting") +
  labs(title = "Aantal bezette IC-bedden") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,800))  +
  geom_hline(aes(yintercept = 1208), color = "grey30", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-06-15"), y = 1280, 
           label = "Gewoonlijke Nederlandse IC-capaciteit", hjust = 1, size = 3, color = "grey30") 
# +
#   geom_hline(aes(yintercept = 2400), color = "black", linetype = "dashed") +
#   annotate(geom = "text", x = as.Date("2020-07-01"), y = 2470, 
#            label = "Opgeschaalde Nederlandse IC-capaciteit", hjust = 1, size = 3)
ggsave(paste0("results/sims_ICprevallnopeak_ribbon_", today(), ".jpg"), units = "mm", width = 200, height = 120, dpi = 300)


figuredata %>% 
  filter(time >= 25 & time < 155) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "met gevolgde\n maatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(prevICU_C)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICprev) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-bezetting") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = ICprev, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "IC-bezetting") +
  labs(title = "Aantal bezette IC-bedden (hypothetisch: bij gelijkgebleven opnamebeleid)") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,2000))  +
  geom_hline(aes(yintercept = 1208), color = "grey30", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-06-15"), y = 1280, 
           label = "Gewoonlijke Nederlandse IC-capaciteit", hjust = 1, size = 3, color = "grey30") 
# +
#   geom_hline(aes(yintercept = 2400), color = "black", linetype = "dashed") +
#   annotate(geom = "text", x = as.Date("2020-07-01"), y = 2470, 
#            label = "Opgeschaalde Nederlandse IC-capaciteit", hjust = 1, size = 3)
ggsave(paste0("results/sims_ICprevallshort_counterfact_ribbon_", today(), ".jpg"), units = "mm", width = 200, height = 120, dpi = 300)


# hospital incidence
obsdata <- tibble(
  tijd = seq(as.Date("2020-02-13"), filedate, 1),
  hospincNICE = inc_nice_hosp(),
  hospincOSIRIS = inc_osiris_hosp()
) %>% bind_cols(
  as_tibble(inc_nice_hosp(ages=T))) %>%
  pivot_longer(4:12, names_to = "ageclass", values_to = "observed")  %>%
  group_by(ageclass) %>%
  mutate(hospcumincNICE = cumsum(hospincNICE),
         hospcumincOSIRIS = cumsum(hospincOSIRIS),
         cumobserved = cumsum(observed)) %>%
  ungroup()

figuredata %>% 
  filter(time >= 25 & time < 350) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "geen verdere\nmaatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(inchosp)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, hospincNICE, hospincOSIRIS) %>%
              distinct() %>%
              pivot_longer(2:3,
                           names_to = "databron",
                           values_to = "hospinc"), by = c("tijd")) %>% 
  mutate(datapunten = "ziekenhuisopnames") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = hospinc, shape = databron), size = 2) +
  scale_shape_discrete(labels = c("NICE", "OSIRIS"), na.translate = F) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijks aantal ziekenhuisopnames") +
  labs(title = "Aantal ziekenhuisopnames per dag") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,600)) 
ggsave(paste0("results/sims_hospincallshort_ribbon_", today(), ".jpg"), units = "mm", width = 240, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 110 & time < 350) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "geen verdere\nmaatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(inchosp)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, hospincNICE, hospincOSIRIS) %>%
              distinct() %>%
              pivot_longer(2:3,
                           names_to = "databron",
                           values_to = "hospinc"), by = c("tijd")) %>% 
  mutate(datapunten = "ziekenhuisopnames") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = hospinc, shape = databron), size = 2) +
  scale_shape_discrete(labels = c("NICE", "OSIRIS"), na.translate = F) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijks aantal ziekenhuisopnames") +
  labs(title = "Aantal ziekenhuisopnames per dag") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,350)) 
ggsave(paste0("results/sims_hospincallnopeak_ribbon_", today(), ".jpg"), units = "mm", width = 200, height = 120, dpi = 300)


# hospital prevalence
obsdata <- tibble(
  tijd = seq(as.Date("2020-02-13"), filedate, 1),
  hospprev = prev_nice_hosp(maxligduur = 100)
) %>% bind_cols(
  as_tibble(prev_nice_hosp(ages=T))) %>%
  pivot_longer(3:11, names_to = "ageclass", values_to = "observed")  

figuredata %>% 
  filter(time >= 25 & time < 350) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "geen verdere\nmaatregelen",
                         kleur)) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(prevhosp + prevICU)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, hospprev) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE ziekenhuisbezetting") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(aes(color = kleur), size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur, color = NULL), alpha = 0.3) +
  geom_point(mapping = aes(y = hospprev, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "ziekenhuisbezetting") +
  labs(title = "Aantal bezette ziekenhuisbedden (inclusief IC)") +
  scale_color_discrete(name = "scenario") +
  scale_fill_discrete(name = "scenario") +
  scale_x_date(date_breaks = "1 month", date_labels = "1 %b") +
  coord_cartesian(ylim = c(0,4500))  
ggsave(paste0("results/sims_hospprevallshort_ribbon_", today(), ".jpg"), units = "mm", width = 240, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  filter(kleur == "huidige pakket") %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "NICE ziekenhuisbezetting") %>%
  ggplot(aes(x = tijd, y = prevhosp + prevICU)) +
  geom_line(size = 1) +
  geom_point(mapping = aes(y = observed, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal ziekenhuisbedden") +
  facet_wrap(~ ageclass, scales = "free_y") +
  labs(title = "Bezette ziekenhuisbedden per leeftijdsgroep")
ggsave(paste0("results/sims_hospprevlft_", today(), ".jpg"), units = "mm", width = 180, height = 120, dpi = 300)

# mortaliteit
obsdata <- tibble(
  tijd = seq(as.Date("2020-02-13"), filedate  , 1),
  mort = inc_osiris_mort()
) %>% bind_cols(
  as_tibble(inc_osiris_mort(ages=T))) %>%
  pivot_longer(3:11, names_to = "ageclass", values_to = "observed") %>%
  group_by(ageclass) %>%
  mutate(cummort = cumsum(mort),
         cumobserved = cumsum(observed)) %>%
  ungroup()

figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  filter(kleur == "huidige pakket") %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(incmort)# + sum(incmorthosp)# + sum(incmortsevere)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = qpois(0.025, quantile(totsimulated, .025, na.rm = T)),
    upperbound = qpois(0.975, quantile(totsimulated, .975, na.rm = T)),
    midbound = median(totsimulated, na.rm = T)
  ) %>% pull(midbound) %>% sum()
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, mort) %>%
              distinct(), by = c("tijd")) %>% 
  mutate(datapunten = "sterfte (OSIRIS)") %>%
  ggplot(aes(x = tijd, y = midbound)) +
  geom_line(size = 2) +
  geom_point(mapping = aes(y = mort, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal sterfgevallen") +
  labs(title = "Sterfte per dag")
ggsave(paste0("results/sims_mortall_", today(), ".jpg"), units = "mm", width = 180, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  filter(kleur == "huidige pakket") %>%
  group_by(ageclass) %>%
  mutate(cumsimulated = cumsum(incmort)) %>%
  ungroup() %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "sterfte (OSIRIS)") %>%
  ggplot(aes(x = tijd, y = cumsimulated)) +
  geom_line(size = 1) +
  geom_point(mapping = aes(y = cumobserved, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Cumulatief aantal sterfgevallen") +
  facet_wrap(~ ageclass, scales = "free_y") +
  labs(title = "Cumulatieve sterfte per leeftijdsgroep")
ggsave(paste0("results/sims_mortlft_", today(), ".jpg"), units = "mm", width = 180, height = 120, dpi = 300)

# immuniteit
figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  filter(kleur == "huidige pakket") %>%
  ggplot(aes(x = tijd, y = 100 * immune)) +
  geom_line(size = 1) +
  theme_light() + 
  labs(x = "Dag", y = "Percentage met doorgemaakte infectie") +
  facet_wrap(~ ageclass, scales = "free_y") +
  labs(title = "Percentage doorgemaakte infecties per leeftijdsgroep")
ggsave(paste0("results/sims_immlft_", today(), ".jpg"), units = "mm", width = 180, height = 120, dpi = 300)



figuredata %>% 
  filter(kleur == "huidige pakket") %>%
  filter(time >= 25 & time < 150) %>%
  group_by(tijd) %>%
  summarise(totmort = sum(incmort)) %>%
  left_join(obsdata, by = "tijd") %>%
  mutate(datapunten = "OSIRIS mortaliteit") %>%
  ggplot(aes(x = tijd, y = totmort)) +
  geom_line(size = 2) +
  geom_point(mapping = aes(y = mort, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijkse mortaliteit") +
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  coord_cartesian(ylim = c(0,400)) 

figuredata %>% 
  filter(kleur == "huidige pakket") %>%
  filter(time >= 25 & time < 150) %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "OSIRIS mortaliteit") %>%
  ggplot(aes(x = tijd, y = incmort)) +
  geom_line(size = 1) +
  geom_point(mapping = aes(y = mortality, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijkse mortaliteit") +
  facet_wrap(~ ageclass, scales = "free_y")
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  coord_cartesian(ylim = c(0,200)) +
  scale_color_discrete(name = "scenario") +
  guides(color = guide_legend(override.aes = list(size = 1))) 



+ 
  # facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208 * popsize()/popsize()), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300 * popsize()/popsize(),
           label = "Nederlandse IC-capaciteit", hjust = 0, size = 4)


ggsave(paste0("results/sims_fitted2_", today(), ".tiff"), units = "mm", width = 180, height = 120, dpi = 150)
obsdata <- tibble(
  tijd = seq(as.Date("2020-02-13"), filedate  , 1),
  ICobserved = prev_nice_ic(),
  ICinc = inc_nice_ic(),
  hosp = inc_osiris_hosp(),
  mort = inc_osiris_mort(),
  report = inc_osiris_report()
) %>% bind_cols(
  as_tibble(inc_osiris_mort(ages=T))) %>%
  pivot_longer(7:15, names_to = "ageclass", values_to = "mortality") 



figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  left_join(obsdata, by = "tijd") %>%
  mutate(datapunten = "NICE IC-opnames NL") %>%
  # filter(kleur == "geen bestrijding") %>%
  # group_by(tijd, kleur) %>%
  ggplot(aes(x = tijd, y = ICUopnames)) +
  geom_line(mapping = aes(color = kleur), size = 2) +
  geom_point(mapping = aes(y = ICinc, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal nieuwe IC-opnames") +
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  coord_cartesian(ylim = c(0,250)) +
  scale_color_discrete(name = "scenario") +
  guides(color = guide_legend(override.aes = list(size = 1))) 
+ 
  # facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208 * popsize()/popsize()), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300 * popsize()/popsize(),
           label = "gem NL IC-capaciteit", hjust = 0, size = 4)
ggsave(paste0("results/sims_fittedinc2_", today(), ".tiff"), units = "mm", width = 180, height = 120, dpi = 150)





eigenvalues <- sapply(1:1000, 
                      function(x) eigen(t(ctrlmats[[x]] * agedist(...)) * 0.9525 * relinfectivity(F)/.946, only.values = T)[[1]][1])

figuredata %>% 
  filter(time >= 25 & time < 250 & kleur == "huidige pakket") %>%
  mutate(
    R0 = Rstarts[repl],
    control = eigenvalues[repl],
    lambda = doublings[repl]) %>%
  # filter(kleur == "geen bestrijding") %>%
  # group_by(tijd, kleur) %>%
  ggplot(aes(x = tijd, y = ICUbezetting, group = repl, color = lambda)) +
  geom_line(size = 0.1) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  scale_y_continuous(limits = c(0,5000)) +
  scale_color_gradient2(name = "doubling\ntime", midpoint = 2.9, low = "red", high = "blue") +
  # guides(color = guide_legend(override.aes = list(size = 1))) + 
  # facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, 
           label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)
ggsave("results/sims_spaghettileisure70_ICU20_doubling.tiff", units = "mm", width = 180, height = 120, dpi = 150)
figuredata %>% 
  filter(time >= 25 & time < 250 & kleur == "huidige pakket") %>%
  mutate(
    R0 = Rstarts[repl],
    control = eigenvalues[repl],
    lambda = doublings[repl]) %>%
  # filter(kleur == "geen bestrijding") %>%
  # group_by(tijd, kleur) %>%
  ggplot(aes(x = tijd, y = ICUbezetting, group = repl, color = R0)) +
  geom_line(size = 0.1) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  scale_y_continuous(limits = c(0,5000)) +
  scale_color_gradient2(name = "Reproduction\nnumber", midpoint = 2.2, high = "red", low = "blue") +
  # guides(color = guide_legend(override.aes = list(size = 1))) + 
  # facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, 
           label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)
ggsave("results/sims_spaghettileisure70_ICU20_reprratio.tiff", units = "mm", width = 180, height = 120, dpi = 150)
figuredata %>% 
  filter(time >= 25 & time < 250 & kleur == "huidige pakket") %>%
  mutate(
    R0 = Rstarts[repl],
    control = eigenvalues[repl],
    lambda = doublings[repl]) %>%
  # filter(kleur == "geen bestrijding") %>%
  # group_by(tijd, kleur) %>%
  ggplot(aes(x = tijd, y = ICUbezetting, group = repl, color = control)) +
  geom_line(size = 0.1) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  scale_y_continuous(limits = c(0,5000)) +
  scale_color_gradient2(name = "Contact\nreduction", midpoint = 0.469, low = "blue", high = "red") +
  # guides(color = guide_legend(override.aes = list(size = 1))) + 
  # facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, 
           label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)
ggsave("results/sims_spaghettileisure70_ICU20_control.tiff", units = "mm", width = 180, height = 120, dpi = 150)


median(eigenvalues)

distsdata <- allsimulsummaries[[1]]
distsdata %>% 
  # filter(scenario %in% c("zonder bestrijding", "doorgaan met huidige pakket")) %>%
  ggplot(aes(x = agecat, y = immuniteit, fill = scenario)) +
  geom_col(position = position_dodge())
ggsave("results/sims_packhqlock_immunity.tiff", units = "mm", width = 180, height = 120, dpi = 150)
distsdata %>% 
  # filter(scenario %in% c("zonder bestrijding", "doorgaan met huidige pakket")) %>%
  ggplot(aes(x = agecat, y = ICopnames, fill = scenario)) +
  geom_col(position = position_dodge())
ggsave("results/sims_packhqlock_ICU.tiff", units = "mm", width = 180, height = 120, dpi = 150)

figuredata %>%
  filter(time >= 25 & time < 250) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = quantile(ICUbezetting, .025),
    upperbound = quantile(ICUbezetting, .975),
    midbound = median(ICUbezetting)
  ) %>%
  ungroup() %>%
  ggplot(aes(x = tijd, y = midbound, color = kleur)) +
  geom_line(size = 2) +
  geom_ribbon(aes(ymin = lowerbound, ymax = upperbound, fill = kleur), alpha = 0.5) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  coord_cartesian(ylim = c(0,5000)) +
  # guides(color = guide_legend(override.aes = list(size = 1))) + 
  # facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, 
           label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)

figuredata %>%
  filter(time >= 25 & time < 250) %>%
  group_by(kleur, repl) %>%
  mutate(maxIC = max(ICUbezetting)) %>%
  group_by(tijd, kleur) %>%
  mutate(ICpct = dense_rank(maxIC)) %>% 
  filter(ICpct %in% c(seq(50,950,50))) %>%
  ungroup() %>%
  mutate(repl = repl + 1000 * (kleur == "huidige pakket")) %>%
  ggplot(aes(x = tijd, y = ICUbezetting, color = kleur, group = repl)) +
  geom_line(size = 1) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
  # scale_x_date(date_labels = "%d-%m") + 
  # scale_x_date(date_breaks = "1 week", date_labels = "%d-%m") +
  coord_cartesian(ylim = c(0,5000)) +
  # guides(color = guide_legend(override.aes = list(size = 1))) + 
  # facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1300, 
           label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)
figuredata %>%
  filter(time >= 25 & time < 250) %>%
  group_by(kleur, repl) %>%
  mutate(maxIC = max(ICUbezetting)) %>%
  group_by(tijd, kleur) %>%
  mutate(ICpct = dense_rank(maxIC)) %>% 
  filter(ICpct %in% c(seq(50,950,50))) %>%
  filter(kleur == "huidige pakket") %>%
  ungroup() %>%
  select(ICpct, maxIC) %>% distinct() %>% arrange(ICpct)