########################################################
### verse schattingen van alle intervallen en kansen ###
########################################################

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

################################
### eerst alle data-analyses ###
################################
analysisdate <- as.Date("2020-04-09")
source("R/code4ode/Reportingdelays4ode.R")
source("R/code4ode/NICEanalyses4ode.R")
source("R/code4ode/OSIRISanalyses4ode.R")




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

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

### contacts, relative infectivities, natural history and control
source("R/code4ode/ContactsInfectivities4ode.R")

# possibility to add alternative matrices and infectivities
ctrlmats <- readRDS("data/ContactmatricesD3praktijk_midpoint_24mrt2020.rds")
ctrlinfs <- list(
  `00` = 1,
  `99` = 0.946,
  IsoMild = 0.98,  
  HQMild = 0.86
)

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

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

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

### Calibration with control up to "today"
# parameters to change: seeding, fittedrelinfs, probdeath
# endtimes: 32 = Sun 15 March, 40 = Mon 23 March
logLikijk <- function(parms) {
  siminc <- covidcontrolsim(contactcontrol = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                            infectivitycontrol = c("99", "IsoMild", "HQMild"),
                            endtimes = c(32, 40, 400),
                            seeding = exp(parms[1]), 
                            fittedrelinfs = exp(parms[c(2,3,3)])) %>%
    filter(time > 0 & time <= as.numeric(filedate - as.Date("2020-02-12"))) %>%
    group_by(time) %>%
    summarise(totincIC = sum(incICU)) %>% pull(totincIC)
  siminc <- siminc * rev(repdelayICinc[1:length(siminc)])
  - sum(dpois(observedICinc, siminc, log = T))
}
optim(log(c(71,1)), logLikijk) # minloglik = 148.90, ests = 6.46; relinf123 = -.110
optim(log(c(71,1,1)), logLikijk) # minloglik = 133.25, ests = 5.407; relinf1=-.0150; relinf23=-0.303
optim(log(c(71,1,1)), logLikijk) # minloglik = 136.99, ests = 5.934; relinf12=-0.0641; relinf3=-0.707
optim(log(c(71,1,1,1)), logLikijk) # minloglik = 133.25, ests = 5.403; relinf1=-0.0145; relinf2 = -0.306; relinf3=-0.297

eigen(t(ctrlmats[["99"]] * agedist()) * relinfectivity(F, "99") * exp(-.0150) *
        parmsNatHist()$betaA/parmsNatHist()$gammaA, only.values = T)[[1]][1]
eigen(t(ctrlmats[["D3"]] * agedist()) * relinfectivity(F, "IsoMild") * exp(-.303) *
        parmsNatHist()$betaA/parmsNatHist()$gammaA, only.values = T)[[1]][1]
eigen(t(ctrlmats[["D3praktijk_leisure70"]] * agedist()) * relinfectivity(F, "HQMild") * exp(-.303) *
        parmsNatHist()$betaA/parmsNatHist()$gammaA, 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"),
    kleur == 2 ~ c("99", "D3", "D3")
  )
}

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

eindtijden <- function() c(32, 40, 400)

allsimulresults <- list()
allfigdata <- list()
for(repnr in 1:1) {
  allsimulresults <- list()
  for(i in 1:1) {
    for(j in 1:1) {
      for(k in 1:2) {
        simul <- covidcontrolsim(contacts = contactkeus(kleur = k), 
                                 infectivities = infectiviteitkeus(kleur = k),
                                 endtimes = eindtijden(),  
                                 seeding = exp(5.407),  
                                 fittedrelinfs = list(exp(c(-.0191,-.0191,-.0191)),
                                                      exp(c(-.01910,-.458,-.458)))[[k]], 
                                 fittedrelSIs = c(1,1,1))
        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) 
                  ))
}




figuredata <- do.call(rbind, allfigdata) 

###############
### 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 < 150) %>%
  filter(kleur == "huidige pakket") %>%
  group_by(tijd) %>%
  summarise(totsimulated = sum(incICU)) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICinc) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-opnames") %>%
  ggplot(aes(x = tijd, y = totsimulated)) +
  geom_line(size = 2) +
  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")
ggsave(paste0("results/sims_ICincall_", 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(incICU)) %>%
  ungroup() %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "NICE IC-opnames") %>%
  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 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 < 150) %>%
  filter(kleur == "huidige pakket") %>%
  group_by(tijd) %>%
  summarise(totsimulated = sum(prevICU)) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, ICprev) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE IC-bezetting") %>%
  ggplot(aes(x = tijd, y = totsimulated)) +
  geom_line(size = 1) +
  geom_point(mapping = aes(y = ICprev, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "IC-bezetting") +
  labs(title = "Aantal bezette IC-bedden")
ggsave(paste0("results/sims_ICprevall_", today(), ".jpg"), units = "mm", width = 180, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  filter(kleur == "huidige pakket") %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "NICE IC-bezetting") %>%
  ggplot(aes(x = tijd, y = prevICU)) +
  geom_line(size = 1) +
  geom_point(mapping = aes(y = observed, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal IC-bedden") +
  facet_wrap(~ ageclass, scales = "free_y") +
  labs(title = "Bezette IC-bedden per leeftijdsgroep")
ggsave(paste0("results/sims_ICprevlft_", today(), ".jpg"), units = "mm", width = 180, 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 < 150) %>%
  filter(kleur == "huidige pakket") %>%
  group_by(tijd) %>%
  summarise(totsimulated = sum(inchosp)) %>%
  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 = totsimulated)) +
  geom_line(size = 2) +
  geom_point(mapping = aes(y = hospinc, shape = datapunten, color = databron), size = 2) +
  scale_color_discrete(labels = c("NICE", "OSIRIS"), na.translate = F) +
  theme_light() + 
  labs(x = "Dag", y = "Dagelijks aantal ziekenhuisopnames") +
  labs(title = "Aantal ziekenhuisopnames per dag")
ggsave(paste0("results/sims_hospincall_", 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(inchosp)) %>%
  ungroup() %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "NICE ziekenhuisopnames") %>%
  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 ziekenhuisopnames") +
  facet_wrap(~ ageclass, scales = "free_y") +
  labs(title = "Cumulatief aantal ziekenhuisopnames per leeftijdsgroep")
ggsave(paste0("results/sims_hospinclft_", today(), ".jpg"), units = "mm", width = 180, 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 < 150) %>%
  filter(kleur == "huidige pakket") %>%
  group_by(tijd) %>%
  summarise(totsimulated = sum(prevhosp + prevICU)) %>%
  ungroup() %>%
  left_join(obsdata %>%
              select(tijd, hospprev) %>%
              distinct(), by = c("tijd")) %>%
  mutate(datapunten = "NICE ziekenhuisbezetting") %>%
  ggplot(aes(x = tijd, y = totsimulated)) +
  geom_line(size = 1) +
  geom_point(mapping = aes(y = hospprev, shape = datapunten), size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "ziekenhuisbezetting") +
  labs(title = "Aantal bezette ziekenhuisbedden")
ggsave(paste0("results/sims_hospprevall_", today(), ".jpg"), units = "mm", width = 180, 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)) %>%
  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 = 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)