########################################################
### Ad hoc rapport over opbouw immuniteit            ###
###   terug naar D3-praktijkmatrices                 ###
###     omdat nog niet uitgekristalliseerd           ###
########################################################

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

################################
### eerst alle data-analyses ###
################################
analysisdate <- as.Date("2020-04-15")
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")

### possibility to choose  matrices and infectivities
# ctrlmats <- c(readRDS("data/ContactmatrixP3_13apr2020.rds"))
# ctrlmats <- c(readRDS("data/ContactmatricesP3like_midpoint_15apr2020.rds"))
ctrlmats <- c(readRDS("data/ContactmatricesD3praktijk_leisure70_24mrt2020.rds"),
              readRDS("data/ContactmatricesD3praktijk_midpoint_24mrt2020.rds")
)
ctrlinfs <- list(
  `00` = 1,
  `99` = 0.946,
  IsoMild = 0.98,  
  HQMild = 0.86
)
### contacts, relative infectivities, natural history and control
source("R/code4ode/ContactsInfectivities4ode_v1.R")

set.seed(filedate)
devs <- runif(1000, -1, 1)
relinfsamples <- list(
  `00` = 1 + rep(0, 1000),
  `99` = 1 + 0.01 * devs,
  IsoMild = 1 + 0.01 * devs,
  HQMild = 1 + 0.02 * devs
)
rm(devs)

Rstarts <- c(2.2, runif(1000, 2, 2.4))
doublings <- c(2.9, runif(1000, 2.8, 3))


### 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()
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
# endtimes: 32 = Sun 15 March, 40 = Mon 23 March, 81 = Mon 3 May
logLikijk <- function(parms, normalisedby = "Symp", pattern = c(1, 1, 1)) {
  siminc <- covidcontrolsimple(contactcontrol = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                               infectivitycontrol = c("99", "IsoMild", "HQMild"),
                               endtimes = c(32, 40, as.numeric(filedate - as.Date("2020-02-12"))),
                               seeding = exp(parms[1]), 
                               fittedrelinfs = exp(parms[1 + pattern]),
                               normalisedby = normalisedby) %>%
    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))
}
agecaseGOF <- function(parms, normalisedby = "Symp", pattern = c(1, 1, 1)) {
  siminc <- covidcontrolsim(contactcontrol = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                               infectivitycontrol = c("99", "IsoMild", "HQMild"),
                               endtimes = c(32, 40, as.numeric(filedate - as.Date("2020-02-12"))),
                               seeding = exp(parms[1]), 
                               fittedrelinfs = exp(parms[1 + pattern]),
                               normalisedby = normalisedby)
  simsummaries <- siminc %>%
    group_by(ageclass) %>%
    summarise(
      hospital = sum(inchosp),
      icu = sum(incICU),
      dead = sum(incmort)
    )
  print(simsummaries)
  obssummaries <- tibble(
    hospitalNICE = sapply(inc_nice_hosp(ages = TRUE), sum),
    hospitalOSIRIS = sapply(inc_osiris_hosp(ages = TRUE), sum),
    icu = sapply(inc_nice_ic(ages = TRUE), sum),
    dead = sapply(inc_osiris_mort(ages = TRUE), sum)
  ) 
  c(sum(dpois(obssummaries$hospitalNICE, sum(obssummaries$hospitalNICE) * simsummaries$hospital / sum(simsummaries$hospital), log = T)),
    sum(dpois(obssummaries$hospitalOSIRIS, sum(obssummaries$hospitalOSIRIS) * simsummaries$hospital / sum(simsummaries$hospital), log = T)),
    sum(dpois(obssummaries$icu, simsummaries$icu, log = T)),
    sum(dpois(obssummaries$dead, sum(obssummaries$dead) * simsummaries$dead / sum(simsummaries$dead), log = T)))
}

estI0list <- list()
estrelinfslist <- list()
logLikslist <- list()
for(normby in c("Symp", "Susc", "SympSusc",
                "SympInf", "SuscInf", "SympSuscInf")) {
  for(i in c(2,4)) {
    print(paste0(normby, i))
    pat <- list(c(1,1,1), c(1,2,2), c(1,1,2), c(1,2,3))[[i]]
    ijkres <- optim(log(c(100, 1, 1, 1)), logLikijk, normalisedby = normby,
                    pattern = pat)
    estI0list[[paste0(normby, i)]] <- ijkres$par[1]
    estrelinfslist[[paste0(normby, i)]] <- ijkres$par[2:4][pat]
    logLikslist[[paste0(normby, i)]] <- ijkres$value
  }
}
# for(normby in c("Symp")) {
#   for(i in 1:4) {
#     print(paste0(normby, i))
#     pat <- list(c(1,1,1), c(1,2,2), c(1,1,2), c(1,2,3))[[i]]
#     ijkres <- optim(log(c(100, 1, 1, 1)), logLikijk, normalisedby = normby, 
#                     pattern = pat)
#     estI0list[[paste0(normby, i)]] <- ijkres$par[1]
#     estrelinfslist[[paste0(normby, i)]] <- ijkres$par[2:4][pat]
#     logLikslist[[paste0(normby, i)]] <- ijkres$value
#     print(ijkres)
#   }
# }
# for(normby in c("Susc", 
#                 "SympInf")) {
#   for(i in 4:4) {
#     print(paste0(normby, i))
#     pat <- list(c(1,1,1), c(1,2,2), c(1,1,2), c(1,2,3))[[i]]
#     ijkres <- optim(log(c(100, 1, 1, 1)), logLikijk, normalisedby = normby, 
#                     pattern = pat)
#     estI0list[[paste0(normby, i)]] <- ijkres$par[1]
#     estrelinfslist[[paste0(normby, i)]] <- ijkres$par[2:4][pat]
#     logLikslist[[paste0(normby, i)]] <- ijkres$value
#   }
# }
agecaseGOF(c(estI0list$Symp2, estrelinfslist$Symp2), pattern = c(1,2,2))
agecaseGOF(c(estI0list$Symp4, estrelinfslist$Symp4), pattern = c(1,2,3))
agecaseGOF(c(estI0list$Susc2, estrelinfslist$Susc2), normalisedby = "Susc", pattern = c(1,2,2))
agecaseGOF(c(estI0list$Susc4, estrelinfslist$Susc4), normalisedby = "Susc", pattern = c(1,2,3))
agecaseGOF(c(estI0list$SympSusc2, estrelinfslist$Symp2), normalisedby = "SympSusc",  pattern = c(1,2,2))
agecaseGOF(c(estI0list$SympSusc4, estrelinfslist$Symp4), normalisedby = "SympSusc", pattern = c(1,2,3))
agecaseGOF(c(estI0list$SympInf2, estrelinfslist$SympInf2), normalisedby = "SympInf", pattern = c(1,2,2))
agecaseGOF(c(estI0list$SympInf4, estrelinfslist$SympInf4), normalisedby = "SympInf", pattern = c(1,2,3))
agecaseGOF(c(estI0list$SuscInf2, estrelinfslist$SuscInf2), normalisedby = "SuscInf", pattern = c(1,2,2))
agecaseGOF(c(estI0list$SuscInf4, estrelinfslist$SuscInf4), normalisedby = "SuscInf", pattern = c(1,2,3))
agecaseGOF(c(estI0list$SympSuscInf2, estrelinfslist$SympSuscInf2), normalisedby = "SympSuscInf", pattern = c(1,2,2))
agecaseGOF(c(estI0list$SympSuscInf4, estrelinfslist$SympSuscInf4), normalisedby = "SympSuscInf", pattern = c(1,2,3))

save.image(paste0("environments/datalogliks",today(),".RData"))
load(paste0("environments/datalogliks",today()-1,".RData"))

R_effective <- function(contactcontrol = "99", infectivitycontrol = "99",
                        normalisedby = "Symp", fittedlogrelinf = 0) {
  relbeta <- 1/eigen(t(contactmatrix("99") * relinfectivity(normalisedby = normalisedby, scenario = "00")) *
                       agedist() * relsusceptibility(normalisedby = normalisedby), only.values = T)[[1]][1]
  
  relbeta * eigen(t(ctrlmats[[contactcontrol]] * 
                      relinfectivity(normalisedby = normalisedby, scenario = infectivitycontrol)) * 
                    agedist() * relsusceptibility(normalisedby = normalisedby) * 
                    exp(fittedlogrelinf) *
                    (parmsNatHist()$betaC/parmsNatHist()$gammaC +
                       parmsNatHist()$betaD/parmsNatHist()$gammaD), only.values = T)[[1]][1]
}
R_effective("99", "99", "Symp", estrelinfslist[["Symp2"]][1])

R_effective("P3like", "IsoMild", "Symp", estrelinfslist[["Symp2"]][2])
R_effective("P3like", "HQMild", "Symp", estrelinfslist[["Symp2"]][3])
R_effective("P3like_plusprimschool", "HQMild", "Symp", estrelinfslist[["Symp2"]][3])
R_effective("P3like_plusprimschool", "IsoMild", "Symp", estrelinfslist[["Symp2"]][3])



panels_rijen <- c(NA)
panels_kolommen <- c("SympSuscInf", "SympSusc", "SympInf", "SuscInf", "Susc", "Symp")
panels_kleuren <- c("Geen bestrijding", "Huidige programma")#, "lockdown", "thuisquarantaine")

contactkeus <- function(rij = 1, kolom = 1, kleur = 1) {
  case_when(
    kleur == 1 ~ c("99", "99", "99"),
    kleur == 2 ~ c("99", "D3praktijk_leisure70", "D3praktijk_leisure70")
  )
}

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 0:200) {
  allsimulresults <- list()
  for(i in 1:1) {
    for(j in 1:6) {
      for(k in 1:2) {
        simul <- covidcontrolsim(contactcontrol = contactkeus(kleur = k), 
                                 infectivitycontrol = infectiviteitkeus(kleur = k),
                                 endtimes = eindtijden(),  
                                 seeding = exp(estI0list[[paste0(panels_kolommen[j],2)]]),  
                                 Rstart = Rstarts[repnr + 1], 
                                 SI = 3.5*(Rstarts[repnr + 1] - 1) * doublings[repnr + 1]/log(2)/5, 
                                 fittedrelinfs = list(exp(rep(estrelinfslist[[paste0(panels_kolommen[j],2)]][1], 3)),
                                                      exp(estrelinfslist[[paste0(panels_kolommen[j],2)]][c(1,2,3)]))[[k]], 
                                 relinfsamples = relinfsamples,
                                 nr = repnr,
                                 normalisedby = panels_kolommen[j])
        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_16042020")
allfigdata <- readRDS("results/allfigdata_16042020")

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


figuredata %>% 
  filter(tijd <= as.Date("2020-04-17"), kolom == "Symp") %>%
  group_by(repl) %>%
  summarise(difICU = sum(incICU[kleur == "Geen bestrijding"]) - sum(incICU[kleur == "Huidige programma"])) %>%
  ungroup() %>%
  select(repl, difICU) %>%
  distinct() %>%
  summarise(
    lowdif = quantile(difICU, .025),
    meddif = median(difICU),
    uppdif = quantile(difICU, .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(kleur != "Geen bestrijding") %>%
  filter(time >= 25 & time < 150) %>%
  # mutate(kleur = if_else(kleur == "huidige pakket",
  #                        "met gevolgde\n maatregelen",
  #                        kleur)) %>%
  select(-kleur) %>% rename (kleur = kolom) %>%
  group_by(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(incICU)
    ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = quantile(totsimulated, .025),
    upperbound = 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 = 1) +
  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,160)) 
ggsave(paste0("results/sims_ICincall_ribbon_", today(), ".jpg"), units = "mm", width = 200, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 75 & time < 150) %>%
  # 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 = quantile(totsimulated, .025),
    upperbound = 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 = 1) +
  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,10)) 
ggsave(paste0("results/sims_ICincallzoom_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 = quantile(totsimulated, .025, na.rm = T),
    upperbound = 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 < 150) %>%
  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 = quantile(totsimulated, .025),
    upperbound = quantile(totsimulated, .975),
    midbound = median(totsimulated)
  ) %>%
  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,3000))  +
  geom_hline(aes(yintercept = 1208), color = "grey30", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-07-01"), 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_ICprevall_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(tijd, kleur, repl) %>%
  summarise(
    totsimulated = sum(prevICU)
  ) %>%
  group_by(tijd, kleur) %>%
  summarise(
    lowerbound = quantile(totsimulated, .025),
    upperbound = quantile(totsimulated, .975),
    midbound = median(totsimulated)
  ) %>%
  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,600))  +
  geom_hline(aes(yintercept = 1208), color = "grey30", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-07-01"), 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_ICprevallzoom_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 < 150) %>%
  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 = quantile(totsimulated, .025),
    upperbound = quantile(totsimulated, .975),
    midbound = median(totsimulated)
  ) %>%
  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,500)) 
ggsave(paste0("results/sims_hospincall_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(inchosp),
    totsimulated = if_else(kleur == "geen bestrijding", NA_real_, totsimulated)
  ) %>%
  group_by(tijd, kleur, ageclass) %>%
  summarise(
    lowerbound = quantile(totsimulated, .025, na.rm = T),
    upperbound = quantile(totsimulated, .975, na.rm = T),
    midbound = median(totsimulated, na.rm = T)
  ) %>%
  ungroup() %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "NICE ziekenhuisopnames") %>%
  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 ziekenhuisopnames") +
  facet_wrap(~ ageclass, scales = "free_y") +
  labs(title = "Cumulatief aantal ziekenhuisopnames per leeftijdsgroep")

# 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) %>%
  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 = quantile(totsimulated, .025),
    upperbound = quantile(totsimulated, .975),
    midbound = median(totsimulated)
  ) %>%
  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,4000))  
ggsave(paste0("results/sims_hospprevall_ribbon_", today(), ".jpg"), units = "mm", width = 200, height = 120, dpi = 300)

figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  filter(kleur != "geen bestrijding") %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "NICE ziekenhuisbezetting") %>%
  ggplot(aes(x = tijd, y = prevhosp + prevICU)) +
  geom_line(aes(color = kleur), 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 != "geen bestrijding") %>%
  group_by(tijd, kleur) %>%
  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(aes(color = kleur), 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 < 200) %>%
  filter(kleur != "geen bestrijding") %>%
  group_by(ageclass, kleur) %>%
  mutate(cumsimulated = cumsum(incmort)) %>%
  ungroup() %>%
  left_join(obsdata, by = c("tijd", "ageclass")) %>%
  mutate(datapunten = "sterfte (OSIRIS)") %>%
  ggplot(aes(x = tijd, y = cumsimulated)) +
  geom_line(aes(color = kleur), size = 2) +
  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(kolom != "Symp") %>%
  filter(kleur != "Geen bestrijding") %>%
  select(-kleur) %>% rename (kleur = kolom) %>%
  ggplot(aes(x = tijd, y = 100 * immune)) +
  geom_line(aes(color = kleur), 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(time == 400 & kleur != "Geen bestrijding") %>%
  select(ageclass, kolom, immune) %>%
  pivot_wider(names_from = kolom, values_from = immune)

figuredata %>%
  filter(time == 400 & kleur != "Geen bestrijding") %>%
  group_by(ageclass, kolom) %>%
  summarise(lower = quantile(immune, .025),
            middle = median(immune),
            upper = quantile(immune, .975)) %>%
  mutate(resultaat = paste0(round(100 * middle, 1), "% (", 
                            round(100 * lower, 1), " - ",
                            round(100 * upper, 1), ")")) %>%
  select(ageclass, kolom, resultaat) %>%
  pivot_wider(names_from = kolom, values_from = resultaat)


figuredata %>%
  filter(time == 400 & kleur != "Geen bestrijding" & repl == 0) %>%
  select(ageclass, kolom, immune) %>%
  group_by(kolom) %>%
  mutate(immuneaantal = immune * agedist() * popsize()) %>%
  summarise(pctimmune = sum(immuneaantal)/popsize())

figuredata %>%
  filter(time == 400 & kleur != "Geen bestrijding") %>%
  select(ageclass, kolom, immune, repl) %>%
  group_by(kolom, repl) %>%
  mutate(immuneaantal = immune * agedist() * popsize()) %>%
  summarise(pctimmune = sum(immuneaantal)/popsize()) %>%
  group_by(kolom) %>%
  summarise(lower = quantile(pctimmune, .025),
            middle = median(pctimmune),
            upper = quantile(pctimmune, .975))

pivot_wider(names_from = kolom, values_from = immune)


figuredata %>%
  filter(time == 400 & kleur != "Geen bestrijding") %>%
  select(ageclass, kolom, immune) %>%
  ggplot(aes(x = ageclass, y = immune)) +
  geom_col() +
  facet_wrap(~kolom)


figuredata %>% 
  filter(time >= 75 & time < 250) %>%
  filter(kleur != "geen bestrijding") %>%
  ggplot(aes(x = tijd, y = 100 * immune)) +
  geom_line(aes(color = kleur), 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)