library(deSolve)
library(tidyverse)

################################
### INPUT voor de simulaties ###
################################
### 2017 population in 10 year age bands (till 80+)
agedist <- function() {
  c(0.10546312, 0.11800504, 0.12623769, 0.12003092, 0.13795451, 0.14478945, 0.12200229, 0.08077414, 0.04474283)
} 
popsize <- function(N = 17180000) return(N)

### contactmatrix of contact hours (Pienter3), and all control matrices. "99" is < 13 March (no contact reduction)
ctrlmats <- readRDS("data/Contactmatrices_18mrt2020.rds")
ctrlmats <- readRDS("data/Contactmatrices_19mrt2020.rds")
contactmatrix <- function(scenario = "99") {
  if(scenario %in% names(ctrlmats)) {
    return(ctrlmats[[scenario]])
  } else {
    return(ctrlmats[["99"]])
  }
} 

### relative infectiousness by age. "99" is < 13 March 
ctrlinfs <- list(
  `00` = 1,
  `99` = 0.95,#0.77,#0.967,
  IsoMild = 0.96,  #0.968*0.83,#0.96,
  HQMild = 0.84,#0.83*0.84,
  Elderly = c(1, 1, 0.98, 0.98, 0.98, 0.98, 0.98, 0.885, 0.895),
  IsoMildElderly = c(0.96, 0.96, 0.94, 0.94, 0.94, 0.94, 0.94, 0.85, 0.86),
  UK_CI = 0.825,
  UK_CI_HQ = 0.74,
  UK_SDO = 1,
  UK_SD = 1,
  UK_PC = 1
)
relinfectivity <- function(unequal = FALSE, scenario = "00", ...) {
  if(unequal) {
    # inschatting obs China CDC Attack Rate, aangenomen dat AR = kans op symptomen = relinf
    return(c(0.05, 0.1, 0.25, 0.55, 0.75, 0.9, 0.95, 0.95, 0.95) * ctrlinfs[[scenario]])
  } else {
    return(ctrlinfs[[scenario]])
  }
}
    
# default is what we think MRC assumed, based on peak IC rate
probsymptomatic <- function(MRC = TRUE, probsymp = 0.5, ...) {
  if(MRC) {
    return(rep(probsymp, 9))
  } else {
    return(relinfectivity(TRUE))
  }
}

# de default in deze en volgende is % in OSIRIS vroeg in epidemie
probhospitalisation <- function(MRC = TRUE, probhosp = 0.15, ...) {
  if(MRC) {
    return(c(0.001, 0.003, 0.012, 0.032, 0.049, 0.102, 0.166, 0.243, 0.273))
  } else {
    return(rep(probhosp, 9))
  }
}

probICU <- function(MRC = TRUE, probICU = 0.15, ...) {
  if(MRC) {
    return(c(0.05, 0.05, 0.05, 0.05, 0.063, 0.122, 0.274, 0.432, 0.709))
  } else {
    return(rep(probICU, 9))
  }
}

# natural history parameters
parmsNatHist <- function(Rstart = 2.2, SI = 5, ...) {
  return(list(betaA = Rstart / relinfectivity(scenario = "99") / SI,
              gammaA = 1 / SI))
}

### delays
# modelled as negative binomial distributions, shapes give SDs of 4-5 with default means.
shapeI2R <- function(shI2R = 8, ...) return(shI2R) 
shapeI2H <- function(shI2H = 8, ...) return(shI2H) 
shapeI2ICU <- function(shI2ICU = 8, ...) return(shI2ICU) 
# symptom to reporting from OSIRIS
meanI2R <- function(meanS2R = 4, incubation = 5, ...) return(meanS2R + incubation)
# symptom to hospitalisation from OSIRIS
meanI2H <- function(meanS2H = 5, incubation = 5, ...) return(meanS2H + incubation)
# symptom to ICU = hosp + 1 dag
meanI2ICU <- function(meanS2ICU = 6, incubation = 5, ...) return(meanS2ICU + incubation)
# infection to immune = 21 dagen, assumption, not relevant until serosurveys will be carried out
meanI2Immune <- function(meanI2Imm = 21, ...) return(meanI2Imm)

# durations
# modelled as negative binomial distributions, shapes give highly variable stays
shapestayhosp <- function(shDhosp = 2, ...) return(shDhosp)
shapestayICU <- function(shDICU = 2, ...) return(shDICU)
meanstayhosp <- function(meanDhosp = 16, ...) return(meanDhosp)
meanstayICU <- function(meanDICU = 10, ...) return(meanDICU)


#################
### modelcode ###
#################
# compartments: S, I, and cumulative incidence
# I has letter A, to allow more I-compartments later, each with own letter
ODEcompartments <- function() return(c("S", "IA", "cumI"))

### rate of new infections into different age classes
covidlambda <- function(y, contacts, infectiousness, parms) {
  res <- rep(0, 9)
  for(i in ODEcompartments()) {
    if(substr(i, 1, 1) == "I") {
      res <- res + parms[[paste0("beta", substr(i, 2, 2))]] * 
        colSums(y[paste0(i, 1:9)] * infectiousness * contacts) * y[paste0("S", 1:9)]
    }
  }
  names(res) <- paste0("la", 1:9)
  return(res)
}

### model equations of SIR model: S, I, and cumulative incidence
covidODEmodel <- function(t, y, parms) {
  la <- covidlambda(y, parms$contacts, parms$infectiousness, parms)
  
  res <- rep(0, length(y))
  names(res) <- names(y)
  
  # dS/dt
  for(i in 1:9) {
    res[paste0("S", i)] <- -la[paste0("la", i)]
  }
  
  # dIA/dt
  for(i in 1:9) {
    res[paste0("IA", i)] <- la[paste0("la", i)] - parms[["gammaA"]] * y[paste0("IA", i)]
  }

  # dcumI/dt
  for(i in 1:9) {
    res[paste0("cumI", i)] <- la[paste0("la", i)]
  }
  
  # finished
  return(list(res))
}

### function for a single simulation, returning time series of all compartments
covidsim <- function(contacts = c("99", "99"),
                     infectivities = c("99", "99"),
                     endtimes = c(30, 250),
                     seeding = 10, ...) {
  ### scale transmission rates
  relbeta <- 1/eigen(t(contactmatrix(contacts[1]) * agedist()) * relinfectivity(...), only.values = T)[[1]][1]
  nathistparms <- parmsNatHist(...)
  for(i in names(nathistparms)) {
    if(substr(i, 1, 4) == "beta") nathistparms[[i]] <- nathistparms[[i]] * relbeta
  }
  
  ### define all compartments
  compartments <- t(outer(ODEcompartments(), 1:9, paste0))
  
  ### initialize compartments
  initialstate <- rep(0, length(compartments))
  names(initialstate) <- compartments
  imports <- rep(seeding/9, 9) / popsize()
  initialstate[paste0("S", 1:9)] <- agedist() - imports
  initialstate[paste0("IA", 1:9)] <- imports
  
  ### simulate no control
  ressimul <- ode(
    func = covidODEmodel,
    y = initialstate,
    times = seq(0, endtimes[1], 1),
    parms = c(list(contacts = contactmatrix(contacts[1]),
                   infectiousness = relinfectivity(scenario = infectivities[1], ...)),
              nathistparms)
  )
  
  ### simulate controls
  for(i in 2:length(contacts)) {
    ressimul <- rbind(
      head(ressimul, -1),
      ode(
        func = covidODEmodel,
        y = tail(ressimul, 1)[, compartments],
        times = seq(endtimes[i - 1], endtimes[i], 1),
        parms = c(list(contacts = contactmatrix(contacts[i]),
                       infectiousness = relinfectivity(scenario = infectivities[i], ...)),
                  nathistparms)
      )
    )
  }
  
  return(ressimul)
}  
  
### function for a single simulation, starting 12 Feb (day 30 = 13 March, start of control), 
# returning observed numbers of cases in age classes:
#   * symptomatic incidence (could be reported)
#   * hospitalisation incidence
#   * ICU incidence
#   * immunity prevalence
#   * hospitalisation prevalence
#   * ICU prevalence
covidcontrolsim <- function(contacts = c("99", "99"),
                            infectivities = c("99", "99"),
                            endtimes = c(30, 250),
                            seeding = 10, ...) {
  result <- covidsim(contacts = contacts,
                     infectivities = infectivities,
                     endtimes = endtimes,
                     seeding = seeding, ...)
  
  # get daily incidence (new infections)
  cumincidencecurves <- result[, paste0("cumI", 1:9)] * popsize() + seeding/9
  incidencecurves <- rbind(rep(seeding/9/popsize(), 9) ,
                             tail(cumincidencecurves, -1) - head(cumincidencecurves, -1))
  
  # calculate numbers that will be observed (after delay)
  tobesymptomaticcurves <- t(incidencecurves) * probsymptomatic(...)
  tobehospitalisedcurves <- tobesymptomaticcurves * probhospitalisation(...)
  tobeICUcurves <- tobehospitalisedcurves * probICU(...)
  
  # calculate observation incidences (applying delays)
  incsymptomaticcurves <- sapply(1:(1 + tail(endtimes, 1)), 
                                 function(x) colSums(t(tobesymptomaticcurves[, 1:x, drop = FALSE]) * 
                                                       dnbinom((x-1):0, shapeI2R(...), mu = meanI2R(...))))
  inchospitalisedcurves <- sapply(1:(1 + tail(endtimes, 1)), 
                                  function(x) colSums(t(tobehospitalisedcurves[, 1:x, drop = FALSE]) * 
                                                        dnbinom((x-1):0, shapeI2H(...), mu = meanI2H(...))))
  incICUcurves <- sapply(1:(1 + tail(endtimes, 1)), 
                           function(x) colSums(t(tobeICUcurves[, 1:x, drop = FALSE]) * 
                                                 dnbinom((x-1):0, shapeI2ICU(...), mu = meanI2ICU(...))))
  immunecurves <- sapply(1:(1 + tail(endtimes, 1)), 
                         function(x) colSums(cumincidencecurves[1:x, , drop = FALSE] * 
                                               dpois((x-1):0, meanI2Immune()))) 
  
  # calculate observation prevalences (applying durations)
  prevhospitalisedcurves <- sapply(1:(1 + tail(endtimes, 1)), 
                                   function(x) colSums(t(inchospitalisedcurves[, max(1, x - meanstayhosp(...) + 1):x, drop = FALSE])))
  prevICUcurves <- sapply(1:(1 + tail(endtimes, 1)), 
                          function(x) colSums(t(incICUcurves[, max(1, x - meanstayICU(...) + 1):x, drop = FALSE])))
  

  # combine and add names
  toreturn <- as_tibble(t(incsymptomaticcurves)) %>%
    rename_all(~rownames(contactmatrix())) %>%
    mutate(time = 0:tail(endtimes, 1)) %>%
    pivot_longer(1:9, names_to = "ageclass", values_to = "incsymp") %>%
    full_join(
      as_tibble(t(inchospitalisedcurves)) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "inchosp"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(t(incICUcurves)) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "incICU"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(t(immunecurves)) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "immune"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(t(prevhospitalisedcurves)) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "prevhosp"),
      by = c("time", "ageclass")
    ) %>%
    full_join(
      as_tibble(t(prevICUcurves)) %>%
        rename_all(~rownames(contactmatrix())) %>%
        mutate(time = 0:tail(endtimes, 1)) %>%
        pivot_longer(1:9, names_to = "ageclass", values_to = "prevICU"),
      by = c("time", "ageclass")
    ) 
    
  return(toreturn)
}



### Calibration with control up to "today"
# change 'seeding' until cumulative ICU patients
# are similar to observations in NICE
# On 18 March: approx 180 known.
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3", "D3"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 42, unequal = F, MRC = T, Rstart = 2.2)


ijksimulatie %>% 
  group_by(time) %>%
  summarise(totprevICU = sum(prevICU)) %>%
  filter(time == 35)
ijksimulatie %>% 
  filter(time <= 30) %>%
  summarise(totinchosp = sum(inchosp))
ijksimulatie %>% 
  filter(time <= 30) %>%
  summarise(totincsymp = sum(incsymp))



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"),
    kleur == 3 ~ c("99", "D3", "D4"),
    kleur == 4 ~ 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", "IsoMild"),
    kleur == 3 ~ c("99", "IsoMild", "IsoMild"),
    kleur == 4 ~ c("99", "IsoMild", "HQMild")
  )
}

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

allsimulresults <- list()
for(i in 1:1) {
  for(j in 1:1) {
    for(k in 1:4) {
      simul <- covidcontrolsim(contacts = contactkeus(kleur = k), 
                               infectivities = infectiviteitkeus(kleur = k),
                               endtimes = eindtijden(), unequal = F, MRC = T, seeding = 42)
      allsimulresults <- c(allsimulresults,
                           list(simul %>%
                             mutate(rij = panels_rijen[i],
                                    kolom = panels_kolommen[j],
                                    kleur = panels_kleuren[k]))
      )
    }
  }
}


figuredata <- do.call(rbind, allsimulresults)

figuredata %>% 
  filter(time >= 25 & time < 90) %>%
  group_by(time, rij, kolom, kleur) %>%
  summarise(ICUbezetting = sum(prevICU)) %>%
  # filter(scenario %in% c("zonder bestrijding", "doorgaan met huidige pakket")) %>%
  mutate(tijd = as.Date("2020-02-12") + time) %>%
  # mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
  ggplot(aes(x = tijd, y = ICUbezetting, color = kleur)) +
  geom_line(size = 2) +
  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)) +
  # 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_packhqlock_short.tiff", units = "mm", width = 180, height = 120, dpi = 150)

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(time, rij, kolom, kleur) %>%
  summarise(ICUbezetting = sum(prevICU)) %>%
  # filter(scenario %in% c("zonder bestrijding", "doorgaan met huidige pakket")) %>%
  mutate(tijd = as.Date("2020-02-12") + time) %>%
  # mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
  ggplot(aes(x = tijd, y = ICUbezetting, color = kleur)) +
  geom_line(size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
  # scale_x_date(date_labels = "%d-%m") + 
  scale_x_date(date_breaks = "1 month", date_labels = "%d-%m") +
  scale_y_continuous(limits = c(0,5000)) +
  # 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)
gggsave("results/sims_packhqlock_long.tiff", units = "mm", width = 180, height = 120, dpi = 150)

figuredata %>% 
  filter(tijd >= 25 & tijd < 400) %>%
  mutate(tijd = as.Date("2020-02-11") + tijd) %>%
  mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
  ggplot(aes(x = tijd, y = intensivecare, color = scenario)) +
  geom_line(size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
  # scale_x_date(date_labels = "%d-%m") +
  scale_x_date(date_breaks = "1 month", date_labels = "%d-%m") +
  scale_y_continuous(limits = c(0,5000)) +
  facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 1500, 
           label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)
ggsave("results/sims_agedep_year.tiff", units = "mm", width = 180, height = 240, dpi = 150)

figuredata %>% 
  filter(tijd >= 25 & tijd < 400) %>%
  mutate(tijd = as.Date("2020-02-11") + tijd) %>%
  mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
  ggplot(aes(x = tijd, y = hersteld, color = scenario)) +
  geom_line(size = 2) +
  theme_light() + 
  labs(x = "Dag", y = "Aantal herstelde (immune) personen") +
  # scale_x_date(date_labels = "%d-%m") +
  scale_x_date(date_breaks = "1 month", date_labels = "%d-%m") +
  scale_y_continuous(breaks = seq(0, 1.2e7, 2e6), labels = c(0,paste0(c(2,4,6,8,10,12), "mln"))) +
  facet_wrap(~ maatregelen, ncol = 1) +
  geom_hline(aes(yintercept = 17180000 * (1 - 1/1.9)), color = "black", linetype = "dashed") +
  annotate(geom = "text", x = as.Date("2020-03-07"), y = 8800000, 
           label = "Groepsimmuniteitsgrens", hjust = 0, size = 3)
ggsave("results/sims_agedep_recovered.tiff", units = "mm", width = 180, height = 240, dpi = 150)



# maatregelenpakketten <- c(
#   "1: Huidige pakket",
#   "1 + 2: Contactbeperking ouderen (Restaurants en scholen open)",
#   "1 + 2 + 3: Contactbeperking ouderen (Scholen open)")
# huisgenoten <- c("geen quarantaine huisgenoten", "met quarantaine huisgenoten")
# allsimulresults <- list()
# for(i in 1:3) {
#   for(j in 1:2) {
#     simulwithout <- covidcontrolsim(controlmatrix = contactmatrix(), reducedinfectiousness = 1)
#     simulwithtracing <- covidcontrolsim(controlmatrix = controlmatrices[[i]], 
#                                       reducedinfectiousness = controlinfectivities[[1]])
#     simulwithouttracing <- covidcontrolsim(controlmatrix = controlmatrices[[i]], 
#                                          reducedinfectiousness = controlinfectivities[[2]])
#     allsimulresults <- c(allsimulresults,
#                          list(tibble(
#                            tijd = rep(0:500, 3),
#                            incidentie = c(simulwithout[,1], simulwithtracing[,1], simulwithouttracing[,1]),
#                            ziekenhuis = c(simulwithout[,2], simulwithtracing[,2], simulwithouttracing[,2]),
#                            intensivecare = c(simulwithout[,4], simulwithtracing[,4], simulwithouttracing[,4]),
#                            hersteld = c(simulwithout[,5], simulwithtracing[,5], simulwithouttracing[,5]),
#                            scenario = c(rep("geen ingrijpen,\n met contacttracering", 501), rep("wel ingrijpen,\n met contracttracering", 501), 
#                                         rep("wel ingrijpen,\n zonder contacttracering", 501)),
#                            maatregelen = maatregelenpakketten[i],
#                            quarantaine = huisgenoten[j]
#                          ))
#     )
#   }
#   
# }
# figuredata <- bind_rows(allsimulresults[[1]], allsimulresults[[2]]) %>%
#   bind_rows(allsimulresults[[3]]) 
# %>%
#   bind_rows(allsimulresults[[4]]) %>%
#   bind_rows(allsimulresults[[5]]) %>%
#   bind_rows(allsimulresults[[6]])
# 
# figuredata %>% 
#   filter(tijd >= 25 & tijd < 90) %>%
#   mutate(tijd = as.Date("2020-02-12") + tijd) %>%
#   mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
#   ggplot(aes(x = tijd, y = intensivecare, color = scenario)) +
#   geom_line(size = 2) +
#   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,2000)) +
#   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_agedep_short.tiff", units = "mm", width = 180, height = 240, dpi = 150)
# 
# figuredata %>% 
#   filter(tijd >= 25 & tijd < 250) %>%
#   mutate(tijd = as.Date("2020-02-11") + tijd) %>%
#   mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
#   ggplot(aes(x = tijd, y = intensivecare, color = scenario)) +
#   geom_line(size = 2) +
#   theme_light() + 
#   labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
#   # scale_x_date(date_labels = "%d-%m") +
#   scale_x_date(date_breaks = "1 month", date_labels = "%d-%m") +
#   scale_y_continuous(limits = c(0,5000)) +
#   facet_wrap(~ maatregelen, ncol = 1) +
#   geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
#   annotate(geom = "text", x = as.Date("2020-03-07"), y = 1500, 
#            label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)
# ggsave("results/sims_agedep_long.tiff", units = "mm", width = 180, height = 240, dpi = 150)
# 
# figuredata %>% 
#   filter(tijd >= 25 & tijd < 400) %>%
#   mutate(tijd = as.Date("2020-02-11") + tijd) %>%
#   mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
#   ggplot(aes(x = tijd, y = intensivecare, color = scenario)) +
#   geom_line(size = 2) +
#   theme_light() + 
#   labs(x = "Dag", y = "Aantal bezette IC-plaatsen") +
#   # scale_x_date(date_labels = "%d-%m") +
#   scale_x_date(date_breaks = "1 month", date_labels = "%d-%m") +
#   scale_y_continuous(limits = c(0,5000)) +
#   facet_wrap(~ maatregelen, ncol = 1) +
#   geom_hline(aes(yintercept = 1208), color = "black", linetype = "dashed") +
#   annotate(geom = "text", x = as.Date("2020-03-07"), y = 1500, 
#            label = "Nederlandse IC-capaciteit", hjust = 0, size = 3)
# ggsave("results/sims_agedep_year.tiff", units = "mm", width = 180, height = 240, dpi = 150)
# 
# figuredata %>% 
#   filter(tijd >= 25 & tijd < 400) %>%
#   mutate(tijd = as.Date("2020-02-11") + tijd) %>%
#   mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
#   ggplot(aes(x = tijd, y = hersteld, color = scenario)) +
#   geom_line(size = 2) +
#   theme_light() + 
#   labs(x = "Dag", y = "Aantal herstelde (immune) personen") +
#   # scale_x_date(date_labels = "%d-%m") +
#   scale_x_date(date_breaks = "1 month", date_labels = "%d-%m") +
#   scale_y_continuous(breaks = seq(0, 1.2e7, 2e6), labels = c(0,paste0(c(2,4,6,8,10,12), "mln"))) +
#   facet_wrap(~ maatregelen, ncol = 1) +
#   geom_hline(aes(yintercept = 17180000 * (1 - 1/1.9)), color = "black", linetype = "dashed") +
#   annotate(geom = "text", x = as.Date("2020-03-07"), y = 8800000, 
#            label = "Groepsimmuniteitsgrens", hjust = 0, size = 3)
# ggsave("results/sims_agedep_recovered.tiff", units = "mm", width = 180, height = 240, dpi = 150)



figuredata %>% 
  filter(tijd > 15 & tijd <= 45) %>%
  mutate(tijd = as.Date("2020-02-11") + tijd) %>%
  mutate(maatregelen = factor(maatregelen, levels = maatregelenpakketten)) %>%
  ggplot(aes(x = tijd, y = incidentie, color = scenario)) +
  geom_line() +
  theme_light() + 
  labs(x = "tijd in dagen", "meldingen per dag") +
  annotate(geom = "text", x = 30, y = 120, label = "vandaag: 78 meldingen/dag", angle = 90, hjust = 0) +
  annotate(geom = "text", x = 37, y = 240, label = "na 1 week: 207 meldingen/dag", angle = 90, hjust = 0) 

  

simul <- covidcontrolsim(controlmatrix = contactmatrix * 0.6, reducedinfectiousness = 1)
plot(simul)
plot(sapply(1:151, function(x) sum(simul[1:x] * dnbinom((x-1):0, 8, mu = 8))))

xstart = c(sus = (1-1E-09)*pop, inf = 1E-09*pop)

parms <- list(beta=.2,gamma=.1, C = contactmatrix)
times <- seq(from=0,to=250,by=.1)

closed.sir.age.model <- function (t, y, parms) {
  S <- y[grepl("sus", names(y))]
  I <- y[grepl("inf", names(y))]
  beta <- parms$beta
  gamma <- parms$gamma
  pop <- parms$pop
  contactmatrix <- parms$C
  #beta <- beta/(eigen(pop*contactmatrix)$values[1])
  #beta <- beta/14.24827
  beta <- beta/14.34231 # eigen(G)$values[1] for population2 and contactmatrix2w
  #print(beta)
  ## now code the model equations
  lambda <- beta*rowSums(S*t(I*contactmatrix))
  dSdt <- -lambda
  dIdt <- lambda - gamma*I
  ## combine results into a single vector
  dxdt <- c(dSdt,dIdt)
  ## return result as a list!
  list(dxdt)
}


res <-
  ode(
    func=closed.sir.age.model,
    y=xstart,
    times=times,
    parms=parms
  ) %>%
  as.data.frame()