###############################
### fitten op IC-incidentie ###
###############################

##################################################
### reporting delay to correct incidence curve ###
##################################################

nicefiles <- grep("2020", grep("NICE_", list.files("data"), value = T), value = T)
prevdata <- read_csv(paste0("data/", nicefiles[2])) %>%
  mutate(ic_opname = if(exists("ic_opname", where = .)) ic_opname else adm_date_icu,
         adm_date_icu = as.Date(ic_opname, "%d/%m/%Y")) %>%
  group_by(adm_date_icu) %>% 
  summarise(aantal = n())
repdelays <- rep(0, 20)
for(fi in sort(nicefiles)[3:11]) {
  newdata <- read_csv(paste0("data/", fi)) %>%
    mutate(ic_opname = if(exists("ic_opname", where = .)) ic_opname else adm_date_icu,
           covid19status = if(exists("covid19status", where = .)) covid19status else "lab",
           adm_date_icu = as.Date(ic_opname, "%d/%m/%Y")) %>%
    filter(covid19status %in% c("lab", "ct")) %>%
    group_by(adm_date_icu) %>% 
    summarise(aantal = n())
  repdelays <- repdelays + 
    full_join(prevdata, newdata, by = "adm_date_icu") %>% 
    arrange(adm_date_icu) %>% 
    mutate(
      aantal.x = if_else(is.na(aantal.x), 0L, aantal.x),
      newentries = aantal.y - aantal.x) %>%
    pull(newentries) %>%
    tail(20)
  prevdata <- newdata
}

repdelayincdist <- c(cumsum(rev(repdelays)/sum(repdelays)), rep(1,80))

###################################################
### reporting delay to correct prevalence curve ###
###################################################

nicefiles <- grep("2020", grep("NICE_", list.files("data"), value = T), value = T)
prevdata <- read_csv(paste0("data/", nicefiles[2])) %>%
  mutate(ic_opname = if(exists("ic_opname", where = .)) ic_opname else adm_date_icu,
         ic_ontslag = if(exists("ic_ontslag", where = .)) ic_ontslag else dis_date_icu,
         adm_date_icu = as.Date(ic_opname, "%d/%m/%Y"),
         dis_date_icu = as.Date(ic_ontslag, "%d/%m/%Y")) 
prevprevs <- sapply(seq(
  as.Date("2020-02-12"), as.Date(substr(nicefiles[2], 6, 13), "%d%m%Y"), by = 1), function(day) {
    sum(prevdata$adm_date_icu <= day & (is.na(prevdata$dis_date_icu) | prevdata$dis_date_icu > day))
  })
prevdprevs <- tail(prevprevs, -1) - head(prevprevs, -1)

repdelays <- rep(0, 20)
for(fi in sort(nicefiles)[3:11]) {
  newdata <- read_csv(paste0("data/", fi)) %>%
    mutate(ic_opname = if(exists("ic_opname", where = .)) ic_opname else adm_date_icu,
           ic_ontslag = if(exists("ic_ontslag", where = .)) ic_ontslag else dis_date_icu,
           covid19status = if(exists("covid19status", where = .)) covid19status else "lab",
           adm_date_icu = as.Date(ic_opname, "%d/%m/%Y"),
           dis_date_icu = as.Date(ic_ontslag, "%d/%m/%Y")) %>%
    filter(covid19status %in% c("lab", "ct"))
  newprevs <- sapply(seq(
    as.Date("2020-02-12"), as.Date(substr(fi, 6, 13), "%d%m%Y"), by = 1), function(day) {
      sum(newdata$adm_date_icu <= day & (is.na(newdata$dis_date_icu) | newdata$dis_date_icu > day))
    })
  newdprevs <- tail(newprevs, -1) - head(newprevs, -1)
  
  repdelays <- repdelays + rev(newdprevs)[1:20] - c(0, rev(prevdprevs)[1:19])
  prevdprevs <- newdprevs
}

repdelayprevdist <- c(cumsum(repdelays/sum(repdelays)), rep(1,80))


################
### ICU-data ###
################

# read data
ICUdata <- read_csv("data/NICE_29032020.csv") %>%
  mutate(age = if_else(age == "NULL", NA_real_, as.numeric(age)),
         adm_date_icu = as.Date(ic_opname, "%d/%m/%Y"),
         dis_date_icu = as.Date(ic_ontslag, "%d/%m/%Y"),
         discharged_to = case_when(
           discharged_to == "NULL" ~ "inICU",
           discharged_to %in% as.character(c(1, 3, 4, 6, 8)) ~ "recovery",
           discharged_to == "7" ~ "death",
           discharged_to == "9" ~ "unknown",
           discharged_to %in% as.character(c(2, 5)) ~ "otherICU",
           TRUE ~ "inICU"
         ))

ICinc <- function(day, lastday = today() - 1, corr = F) {
    selection <- ICUdata %>%
      filter(covid19status %in% c("lab", "ct")) %>%
      select(adm_date_icu, dis_date_icu, discharged_to)
    if(corr) {
      (sum(selection$adm_date_icu == day) - sum(selection$dis_date_icu == day & selection$discharged_to == "otherICU")) / 
        repdelayincdist[as.numeric(lastday - day + 1)]
    } else {
      sum(selection$adm_date_icu == day) - sum(selection$dis_date_icu == day & selection$discharged_to == "otherICU")
    }
  }

ICinccurve <- function(lastday = today() - 1, corr = F) {
  round(sapply(seq(as.Date("2020-02-13"), lastday, 1), ICinc, corr = corr, lastday = lastday))
}

ICprev <- function(day) {
  selection <- NICEdata %>%
    filter(covid19status %in% c("lab", "ct")) %>%
    select(adm_date_icu, dis_date_icu)
  sum(selection$adm_date_icu <= day & (is.na(selection$dis_date_icu) | selection$dis_date_icu > day))
}
ICprevcurve <- function(firstday = as.Date("2020-02-13"), lastday = today() - 1, corr = F) {
  uncorrected <- sapply(seq(firstday, lastday, 1), ICprev)
  if(corr) {
    dprevuncorrected <- tail(uncorrected, -1) - head(uncorrected, -1)
    dprevcorrected <- dprevuncorrected / rev(repdelayprevdist[1:length(dprevuncorrected)])
    return(round(cumsum(c(uncorrected[1], dprevcorrected))))
  } else {
    return(uncorrected)
  }
}

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(...) {
  17180000
}

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

### relative infectiousness by age. "99" is < 13 March 
ctrlinfs <- list(
  `00` = 1,
  `99` = 0.946,#0.77,#0.967,
  IsoMild = 0.98,  #0.968*0.83,#0.96,
  HQMild = 0.86,#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", fittedrelinf = 1, ...) {
  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]] * fittedrelinf)
  } else {
    return(ctrlinfs[[scenario]] * fittedrelinf)
  }
}

# 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))
  }
}

### IN NL VEEL MINDER IC VOOR 80+ (NICE)
# MRC: (0.243*0.432)/(0.273*0.709) = 0.542 keer 'grotere' kans voor 70-80 naar IC dan voor 80+
# NICE: 497 in 70-80, 61 in 80+. Omgerekend naar demografie: (149/.0808)/(29/.0447) = 4.5 grotere kans. 
# verhouding tussen MRC en NICE is 4.5/.542 = 8.3. Gebruik nu: 6 (infectie-inc ook lager in die groep)
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 / 6))
  } 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
# ICUstay from data
shapestayhosp <- function(shDhosp = 2, ...) return(shDhosp)
shapestayICU <- function(shDICU = 1, ...) return(shDICU)
meanstayhosp <- function(meanDhosp = 16, ...) return(meanDhosp)
meanstayICU <- function(meanDICU = 24, ...) 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, fittedrelinf = 1, fittedrelinf2 = 1, ...) {
  ### 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 * fittedrelinf
  }

  ### 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
  nathistparms$betaA <- nathistparms$betaA * fittedrelinf2
  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[, 1:x, drop = FALSE]) * 
                                                         pnbinom((x-1):0, shapestayhosp(...), 
                                                                 mu = meanstayhosp(...), lower.tail = F)))
  prevICUcurves <- sapply(1:(1 + tail(endtimes, 1)), 
                          function(x) colSums(t(incICUcurves[, 1:x, drop = FALSE]) * 
                                                pnbinom((x-1):0, shapestayICU(...), 
                                                        mu = meanstayICU(...), lower.tail = F)))
  # 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 current ICU patients
# are similar to observations in NICE
# NL: On 26 March: approx 826 known.
logLikijk <- function(parms) {
  # siminc <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
  #                           infectivities = c("99", "IsoMild", "HQMild"),
  #                           endtimes = c(30, 41, 400),
  #                           seeding = exp(parms[1]), unequal = F, MRC = T, Rstart = 2.2, meanDICU = 24,
  #                           fittedrelinf = exp(parms[2])) %>%
  #   filter(time > 0 & time <= 46) %>%
  #   group_by(time) %>%
  #   summarise(totincIC = sum(incICU)) %>% pull(totincIC)
  # - sum(dpois(ICinccurve(corr = T), siminc, log = T))
  simprev <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                            infectivities = c("99", "IsoMild", "HQMild"),
                            endtimes = c(30, 41, 400),
                            seeding = exp(parms[1]), unequal = F, MRC = T, Rstart = 2.2, meanDICU = 24,
                            fittedrelinf = exp(parms[2]), fittedrelinf2 = exp(parms[3])) %>%
    filter(time > 0 & time <= 46) %>%
    group_by(time) %>%
    summarise(totprevICU = sum(prevICU)) %>% pull(totprevICU)
  simdprev <- c(simprev[1], tail(simprev, -1) - head(simprev, -1))
  dataset <- ICprevcurve(corr = T)
  dataset <- c(dataset[1], tail(dataset, -1) - head(dataset, -1))
  - sum(dpois(dataset, simdprev, log = T))
}
optim(log(c(71,1,1)), logLikijk)

eigen(t(ctrlmats[["99"]] * agedist()) * relinfectivity(F, "99", 1.124) * 
        parmsNatHist()$betaA/parmsNatHist()$gammaA, only.values = T)[[1]][1]
eigen(t(ctrlmats[["D3praktijk_leisure70"]] * agedist()) * relinfectivity(F, "IsoMild", 1.124) * 
        parmsNatHist()$betaA/parmsNatHist()$gammaA, only.values = T)[[1]][1]
eigen(t(ctrlmats[["D3praktijk_leisure70"]] * agedist()) * relinfectivity(F, "HQMild", 1.124) * 
        parmsNatHist()$betaA/parmsNatHist()$gammaA, only.values = T)[[1]][1]

ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "HQMild"),
                                endtimes = c(30, 41, 400),
                                seeding = 16.4, unequal = F, MRC = T, Rstart = 2.2, shDICU = 1, meanDICU = 24,
                                fittedrelinf = 1.125, fittedrelinf2 = 1.2)
ijksimulatie %>% 
  group_by(time) %>%
  summarise(totprevICU = sum(prevICU)) %>%
  filter(time == 46)
ijksimulatie %>% 
  filter(time <= 80) %>% 
  group_by(time) %>%
  summarise(totincIC = sum(incICU)) %>% pull(totincIC) %>% round()
ijksimulatie %>% 
  filter(time > 0 & time <= 70) %>% 
  group_by(time) %>%
  summarise(totprevICU = sum(prevICU)) %>% pull(totprevICU) %>% round()
ijksimulatie %>% 
  filter(time == 46) %>% 
  summarise(totimmune = sum(immune)) %>% pull(totimmune) %>% round()
ijksimulatie %>% 
  filter(time <= 46) %>% 
  group_by(time) %>%
  summarise(totinc = sum(incsymp)) %>% pull(totinc) %>% round() %>% sum()
# regio 1: On 23 March: approx 20 known.
ICprevcurve(1, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 2.1, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20,
                                ROAZ = 1)
# regio 2: On 23 March: approx 32 known.
ICprevcurve(2, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 3.4, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20,
                                ROAZ = 2)
# regio 3: On 23 March: approx 49 known.
ICprevcurve(3, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 5.2, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20,
                                ROAZ = 3)
# regio 4: On 23 March: approx 83 known.
ICprevcurve(4, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 9.4, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20, ROAZ = 4)
# regio 1: On 23 March: approx 20 known.
ICprevcurve(1, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 60, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20)
# regio 1: On 23 March: approx 20 known.
ICprevcurve(1, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 60, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20)
# regio 1: On 23 March: approx 20 known.
ICprevcurve(1, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 60, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20)
# regio 1: On 23 March: approx 20 known.
ICprevcurve(1, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 60, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20)
# regio 1: On 23 March: approx 20 known.
ICprevcurve(1, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 60, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20)
# regio 1: On 23 March: approx 20 known.
ICprevcurve(1, finalday = as.Date("2020-03-23"))
ijksimulatie <- covidcontrolsim(contacts = c("99", "D3praktijk_leisure70", "D3praktijk_leisure70"),
                                infectivities = c("99", "IsoMild", "IsoMild"),
                                endtimes = c(30, 40, 400),
                                seeding = 60, unequal = F, MRC = T, Rstart = 2.2, nr=1, meanDICU = 20)
ijksimulatie %>% 
  group_by(time) %>%
  summarise(totprevICU = sum(prevICU)) %>%
  filter(time == 40)


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

plot(ijksimulatie$immune)

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", "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(30, 41, 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(), unequal = F, MRC = T, 
                                 seeding = 16.4, fittedrelinf = 1.125)
        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) %>%
                         group_by(time, rij, kolom, kleur) %>%
                         summarise(ICUbezetting = sum(prevICU),
                                   ICUopnames = sum(incICU)) %>%
                         ungroup() %>%
                         mutate(tijd = as.Date("2020-02-12") + time) %>%
                         select(time, tijd, ICUbezetting, ICUopnames, kleur) %>%
                         mutate(replicate = repnr)
                  ))
}
# saveRDS(allfigdata, "results/allfigdata_leisure80_IC20")
# saveRDS(allfigdata, "results/allfigdata_leisure70_IC20")
# saveRDS(allfigdata, "results/allfigdata_leisure70")
# saveRDS(allfigdata, "results/allfigdata_leisure80")
allfigdata <- readRDS("results/allfigdata_leisure70_IC20")






figuredata <- do.call(rbind, allfigdata) %>%
  rename(repl = `replicate`)

obsdata <- tibble(
  tijd = seq(as.Date("2020-02-13"), today() - 1 , 1),
  ICobserved = ICprevcurve(lastday = today() - 1),#ICprevcurve("all", finalday = today() - 1),
  ICinc = ICinccurve(lastday = today() - 1)#c(rep(0,15),1,0,1,0,0,3,1,2,3,5,4,8,6,21,14,17,26,33,34,51,58,86,69,69,102,108,108,92,92,33)
)

figuredata %>% 
  filter(time >= 25 & time < 150) %>%
  left_join(obsdata, by = "tijd") %>%
  mutate(datapunten = "NICE IC-bezetting NL") %>%
  # filter(kleur == "geen bestrijding") %>%
  # group_by(tijd, kleur) %>%
  ggplot(aes(x = tijd, y = ICUbezetting)) +
  geom_line(mapping = aes(color = kleur), size = 2) +
  geom_point(mapping = aes(y = ICobserved, shape = datapunten), 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") +
  coord_cartesian(ylim = c(0,2500)) +
  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_fitted_", today(), ".tiff"), units = "mm", width = 180, height = 120, dpi = 150)



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_fittedinc_", 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)