###########################################################
### Analyse voor ___UITZONDERINGSGROND_2___, data tot 5 juni         ###
###   SEEIIR-model, SI = 4 dagen                        ###
###   Leeftijdsmodel via kans op symptomen (Symp)       ###
###   Opmerkingen      :                                ###
###     * Contactmatrices obv EpiPose                   ###
###           - matrix 'lockdown' vanaf 16 maart        ###
###                 => EpiPose data, bijgeplust tot     ###
###                       reductie van covid-pienter    ###
###           - matrices 'batch1' en 'batch2' opnieuw   ###
###               ingeschat door Jantien en Don         ###
###     * Periode 2 tot dag 44 ipv 40                   ###
###           - vr 28 maart ipv ma 24 maart             ###
###           - geen reden meer voor 24 maart, omdat    ###
###               er geen daling in R te zien was       ###
###           - iets betere fit dan dag 40              ###
###########################################################

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

################################
### eerst alle data-analyses ###
################################
analysisdate <- as.Date("2020-06-05")
source("R/code4ode/Reportingdelays4ode.R")
source("R/code4ode/NICEanalyses4ode_v3.R")
source("R/code4ode/OSIRISanalyses4ode_v1.R")



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

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

### contacts, relative infectivities, natural history and control
source("R/code4ode/ContactsInfectivities4ode_v2.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")
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)
              )
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), 4)),
  list(
    `00` = 1,
    `99` = 0.946,
    IsoMild = 0.98,  
    HQMild = 0.86
  ))
rm(devs)


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

### model code
source("R/code4ode/modelcode4ode_v1.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", "lockdown", "batch1"),
                               infectivitycontrol = c("99", "IsoMild", "IsoMild", "IsoMild", "IsoMild"),
                               endtimes = c(32, 44, 
                                            as.numeric(filedate - 28 - as.Date("2020-02-12")),
                                            as.numeric(as.Date("2020-05-10") - 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)])) %>%
    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))
}
logLikijk_1 <- function(parms) {
  siminc <- covidcontrolsimple(contactcontrol = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                               infectivitycontrol = c("99", "IsoMild", "IsoMild"),
                               endtimes = c(32, 44, 
                                            as.numeric(filedate - as.Date("2020-02-12"))),
                               seeding = exp(5.823), SI = 4,
                               fittedrelinfs = exp(parms[c(1,2,3)])) %>%
    filter(time > 0 & time <= as.numeric(filedate - as.Date("2020-02-12"))) %>%
    pull(totincICU)
  siminc <- siminc * rev(repdelayICinc[1:length(siminc)])
  - sum(dpois(observedICinc, siminc, log = T))
}
optim(log(c(330,1,.5,1)), logLikijk) # fri = 23444, minloglik = 270.68, ests = 6.142; relinf1=.0337; relinf2 = -0.592; relinf3=-0.005
optim(log(c(330,1,.5,1,1)), logLikijk) # fri = 23455, minloglik = 268.73, ests = 6.189; relinf1=.0273; relinf2 = -0.559; relinf3=-0.028, relinf4=0.182
# 1000 samples met beste model
ijkres <- optim(c(5,  0, -1, 0), 
                logLikijk, hessian = T) # IsoMild, IsoMild, IsoMild
library(mvtnorm)
set.seed(filedate)
rndpars <- rmvnorm(200, ijkres$par, solve(ijkres$hessian))



### 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")
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/ContactmatricesD3praktijk_midpoint_24mrt2020.rds"),
              list(lockdown = Reduce(`+`, ctrlmatslockdown)/200,
                   batch1 = Reduce(`+`, ctrlmatsbatch1)/200)
)
for(i in 1:200) {
  for(j in 0:2) {
    ctrlmats[[i + 200 * j]] <- ctrlmats[[i + 200 * j]] * fittedrelinfmatrix(exp(rndpars[i, 4]), exp(rndpars[i, 2]), rep(1, 9))
  }
}

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


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"),
    kleur == 2 ~ c("99", "lockdown", "lockdown", "batch1", "RND", "RND")
  )
}

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

eindtijden <- function() c(32, 44, 
                           as.numeric(as.Date("2020-05-10") - as.Date("2020-02-12")),
                           as.numeric(filedate - as.Date("2020-02-12")) - 14, 
                           as.numeric(as.Date("2020-06-01") - 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 = list(exp(rndpars[repnr, c(2,2,2,2,2,2)]),
                                                      c(exp(rndpars[repnr, c(2,3,4,4)]), 1, 1))[[k]], 
                                 nr = c(repnr, repnr, repnr, repnr, repnr + 200, repnr + 400))
        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____05062020")
allfigdata <- readRDS("results/allfigdata___UITZONDERINGSGROND_2____05062020")

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

# IC prevented
figuredata %>%
  filter(tijd <= as.Date("2020-06-05")) %>%
  group_by(repl, kleur) %>%
  summarise(totICU = sum(incICU)) %>%
  group_by(repl) %>%
  mutate(
    zonder = totICU[kleur == "geen bestrijding"],
    met = totICU[kleur == "huidige pakket"],
    verschil = zonder - met
  ) %>%
  ungroup() %>%
  select(repl,zonder,met,verschil) %>%
  distinct() %>%
  select(-repl) %>%
  as.matrix() %>%
  apply(2, quantile, probs = c(.025,.5,.975))

# two weeks symptomatic incidence
figuredata %>%
  filter(tijd >= as.Date("2020-05-23") & tijd <= as.Date("2020-06-05") & kleur == "huidige pakket") %>%
  group_by(repl, tijd) %>%
  summarise(totsymp = sum(incsymp)) %>%
  group_by(tijd) %>%
  summarise(
    lower = quantile(totsymp, .025),
    median = median(totsymp),
    upper = quantile(totsymp, .975)
  )
figuredata %>%
  filter(tijd <= as.Date("2020-06-05") & kleur == "huidige pakket") %>%
  mutate(weeknr = 7+floor((time+2)/7)) %>% 
  group_by(repl, weeknr) %>%
  summarise(totsymp = sum(incsymp)) %>%
  group_by(weeknr) %>%
  summarise(
    lower = quantile(totsymp, .025),
    median = median(totsymp),
    upper = quantile(totsymp, .975)
  ) 
figuredata %>%
  filter(tijd >= as.Date("2020-05-23") & tijd <= as.Date("2020-06-05") & kleur == "huidige pakket") %>%
  group_by(repl) %>%
  summarise(totsymp = sum(incsymp)) %>%
  mutate(totinf = totsymp)%>%# * 2^runif(200,1,3)) %>%
  summarise(
    lower = quantile(totinf, .025),
    median = median(totinf),
    upper = quantile(totinf, .975)
  )



###############
### 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 < 155) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "met gevolgde\n maatregelen",
                         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 = 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 < 155) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "met gevolgde\n maatregelen",
                         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-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_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 < 155) %>%
  mutate(kleur = if_else(kleur == "huidige pakket",
                         "met gevolgde\n maatregelen",
                         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 = 200, height = 120, dpi = 300)


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

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(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 = 200, 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) %>%
  summarise(totsimulated = sum(incmort)/200) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, mort) %>%
              distinct(), by = c("tijd")) %>% 
  mutate(datapunten = "sterfte (OSIRIS)") %>%
  ggplot(aes(x = tijd, y = totsimulated)) +
  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)