###################################################################
# Functions to use with estimationRt_OSIRIS script
###################################################################

library(tidyverse)
library(lubridate)

get_reportingdelay_Osiris <- function(var_name = NULL, type = NULL, start_date = as.Date("2020-03-13"), end_date = today()) {
  
  if(length(var_name) == 0) stop("provide Osiris name to analyse")
  if(length(type) == 0) stop("provide type of reports (\"reported\", \"hosp\" or \"IC\")")
  
  # get filenames of all Osiris downloads
  files <- list.files("/___UITZONDERINGSGROND_6___Previous", full.names = TRUE)
  files <- c(files, list.files("/___UITZONDERINGSGROND_6___Geschoond", full.names = TRUE))
  files <- files[grepl(files, pattern = "rds")]
  files <- tibble(file_name = files, 
                  date = files %>% str_extract_all(pattern = "(Data_[0-9]+)") %>% unlist %>% gsub(pattern = "Data_", replacement = "") %>% as.Date(format = "%Y%m%d"),
                  time = files %>% file.info %>% pull(ctime))

  report_dates <- seq(start_date, end_date, by = "day") %>% as.character

  rep_delay <- rep(0,40)
  for(report_date in report_dates) {
    file_name <- files %>% 
      filter(date == report_date) %>%
      filter(time == max(time)) %>% # get most recent file (if multiple files exist) 
      pull(file_name)
    
    # if today's data is not yet available
    if(length(file_name)==0) {
      print(paste("only files up to ", report_date, "used"))
      return(rev(rep_delay))
    }
    
    newdata <- read_rds(file_name) %>% 
    {if(type == "IC") filter(., NCOVopnameICU == "J") else if(type == "hosp") filter(., NCOVpatZhs == "J") else .} %>% 
      mutate(var_date = pull(., var_name)) %>%
      mutate(var_date = as.Date(var_date, tryFormats = c("%d-%m-%Y", "%Y-%m-%d"))) %>% 
      filter(!is.na(var_date) & var_date >= as.Date("2020-02-01") & var_date <= end_date) %>%
      group_by(var_date) %>%
      summarise(incidence = n()) %>% 
      full_join(tibble(var_date = seq(min(min(.$var_date), as.Date(report_date)-length(rep_delay)+1), as.Date(report_date), by = "day"))) %>%
      replace_na(list(incidence = 0)) %>%
      arrange(var_date)
    
    if(as.Date(report_date) > start_date) {
      newdata <- full_join(prevdata, newdata, by = "var_date") %>% 
        replace_na(list(incidence.x = 0, incidence.y = 0)) %>%
        arrange(var_date) %>% 
        mutate(newentries = incidence.y - incidence.x) %>%
        mutate(newentries = ifelse(newentries < 0, 0, newentries)) # later corrections can cause decrease in cumulative numbers, ignored here
      rep_delay <- rep_delay + (newdata %>% pull(newentries) %>% tail(length(rep_delay)))
      newdata <- newdata %>% rename(incidence = incidence.y)  %>% select(var_date, incidence)
    } else {
      startdata <- newdata
    }
    
    prevdata <- newdata
  }
  
  tmp <- full_join(startdata %>% rename(start_incidence = incidence), newdata) %>%
    replace_na(list(start_incidence = 0)) %>%
    mutate(diff = incidence - start_incidence) %>%
    mutate(diff = ifelse(diff < 0, 0, diff)) %>% pull(diff) %>% rev
  
  exp_cases <- tmp[1:length(rep_delay)]/(cumsum(rev(rep_delay))/sum(rep_delay)) - tmp[1:length(rep_delay)]
  
  return(list(repdelay = rev(rep_delay), exp_cases = exp_cases))
}

# not yet adapted to new rep_delay
get_reportingdelay_Osiris_province <- function(province = NULL, var_name = NULL, type = NULL, start_date = as.Date("2020-03-13"), end_date = today()) {
  
  if(length(province) == 0) stop("provide province name to analyse")
  if(length(var_name) == 0) stop("provide Osiris name to analyse")
  if(length(type) == 0) stop("provide type of reports (\"reported\", \"hosp\" or \"IC\")")
  
  # get filenames of all Osiris downloads
  files <- list.files("/___UITZONDERINGSGROND_6___Previous", full.names = TRUE)
  files <- c(files, list.files("/___UITZONDERINGSGROND_6___Geschoond", full.names = TRUE))
  files <- files[grepl(files, pattern = "rds")]
  files <- tibble(file_name = files, 
                  date = files %>% str_extract_all(pattern = "(Data_[0-9]+)") %>% unlist %>% gsub(pattern = "Data_", replacement = "") %>% as.Date(format = "%Y%m%d"),
                  time = files %>% file.info %>% pull(ctime))
  
  report_dates <- seq(start_date, end_date, by = "day") %>% as.character
  
  rep_delay <- rep(0,40)
  for(report_date in report_dates) {
    file_name <- files %>% 
      filter(date == report_date) %>%
      filter(time == max(time)) %>% # get most recent file (if multiple files exist) 
      pull(file_name)
    
    # if today's data is not yet available
    if(length(file_name)==0) {
      print(paste("only files up to ", report_date, "used"))
      return(rev(rep_delay))
    }
    
    newdata <- read_rds(file_name) %>% 
      filter(Provincie == province) %>%
      {if(type == "IC") filter(., NCOVopnameICU == "J") else if(type == "hosp") filter(., NCOVpatZhs == "J") else .} %>% 
      mutate(var_date = pull(., var_name)) %>%
      mutate(var_date = as.Date(var_date, tryFormats = c("%d-%m-%Y", "%Y-%m-%d"))) %>% 
      filter(!is.na(var_date) & var_date >= as.Date("2020-02-01") & var_date <= end_date) %>%
      group_by(var_date) %>%
      summarise(incidence = n()) %>% 
      full_join(tibble(var_date = seq(min(min(.$var_date), as.Date(report_date)-length(rep_delay)+1), as.Date(report_date), by = "day"))) %>%
      replace_na(list(incidence = 0)) %>%
      arrange(var_date)
    
    if(as.Date(report_date) > start_date) {
      newdata <- full_join(prevdata, newdata, by = "var_date") %>% 
        replace_na(list(incidence.x = 0, incidence.y = 0)) %>%
        arrange(var_date) %>% 
        mutate(newentries = incidence.y - incidence.x) %>%
        mutate(newentries = ifelse(newentries < 0, 0, newentries)) # later corrections can cause decrease in cumulative numbers, ignored here
      rep_delay <- rep_delay + (newdata %>% pull(newentries) %>% tail(length(rep_delay)))
      newdata <- newdata %>% rename(incidence = incidence.y)  %>% select(var_date, incidence)
    }
    
    prevdata <- newdata
  }
  return(rev(rep_delay))
}

get_symptomonset2admission_Osiris <- function(IC = TRUE, report_date = today()) {
  files <- list.files("/___UITZONDERINGSGROND_6___Geschoond", full.names = TRUE)
  files <- c(files, list.files("/___UITZONDERINGSGROND_6___Previous", full.names = TRUE))
  files <- grep(files, pattern = "rds", value = TRUE)
  file_name <- grep(files, pattern = report_date %>% gsub(pattern = "-", replacement = ""), value = TRUE)
  file_name <- file_name[file_name %>% file.info %>% pull(ctime) %>% which.max]
  data <- read_rds(file_name)
  
  data <- data %>% 
    mutate(NCOVdat1ezkhopn = NCOVdat1ezkhopn %>% as.Date(format = "%Y-%m-%d")) %>% 
    filter((if(IC) NCOVopnameICU else NCOVpatZhs) == "J") %>%
    {if(IC) mutate(., adm_date = NCOVopnamedatumICU) else mutate(., adm_date = NCOVdat1ezkhopn)}
  
  interval <- data %>% filter(ZIE1eZiekteDt >= as.Date("2020-02-01") & ZIE1eZiekteDt <= report_date &
                                adm_date >= as.Date("2020-02-27") & adm_date <= report_date) %>%
    mutate(interval = as.integer(adm_date - ZIE1eZiekteDt)) %>%
    pull(interval)
  
  interval_cdf <- ecdf(interval[interval >= 0]) # don't trust symptom onset after day of admission
  interval_cdf <- interval_cdf(0:100)
  interval_pdf <- diff(c(0, interval_cdf))
  interval_pdf <- interval_pdf[1:max(which(interval_pdf > 0))]
  
  return(interval_pdf)
}

# old
extract_epicurve_Osiris <- function(data, IC = FALSE, report_date = today(), SOtoRep = NULL) {
  data <- data %>%
  {if(IC) filter(., NCOVopnameICU == "J") else filter(., NCOVpatZhs == "J")} %>%
    mutate(bron = if("Bron" %in% names(data)) Bron else Bron_land)

  epicurve <- data %>% 
    filter(!is.na(ZIE1eZiekteDt) & ZIE1eZiekteDt <= today()) %>% # omit symptom onset dates in future
    group_by(ZIE1eZiekteDt) %>%
    summarise(incidenceABROAD = sum(bron == "Buitenland", na.rm = TRUE),
              incidenceTOT = n()) %>%
    mutate(incidenceNL = incidenceTOT - incidenceABROAD) %>% 
    full_join(tibble(ZIE1eZiekteDt = seq(min(.$ZIE1eZiekteDt), as.Date(report_date), by = "day"))) %>%
    replace_na(list(incidenceNL = 0, incidenceTOT = 0, incidenceABROAD = 0)) %>%
    arrange(ZIE1eZiekteDt) %>%
    mutate(dates = ZIE1eZiekteDt) %>%
    select(dates, incidenceNL, incidenceABROAD, incidenceTOT)

  if(length(SOtoRep) == 0) {
    # default SO to reporting delay over last 7 days
    reporting_delay <- get_symptomonset2reporting(start_date = as.Date(report_date) - 7, end_date = as.Date(report_date), IC = IC)
  } else {
    reporting_delay = SOtoRep
  }
  
  reporting_delay <- cumsum(reporting_delay)/sum(reporting_delay)
  
  if(nrow(epicurve) > length(reporting_delay)) {
    reporting_delay <- c(rep(1, nrow(epicurve)-length(reporting_delay)), rev(reporting_delay))
  } else {
    reporting_delay <- rev(reporting_delay[1:nrow(epicurve)])
  }
  
  epicurve <- epicurve %>% mutate(p = reporting_delay)
  
  # if most recent days are zero, give them one case and a p to reflect total number of last 3 non-zero days
  max_nonzero_inc_day <- max(which(epicurve$incidenceTOT != 0))
  if(max_nonzero_inc_day < nrow(epicurve)) {
    epicurve$incidenceTOT[(max_nonzero_inc_day + 1):nrow(epicurve)] <- 1
    epicurve$incidenceNL[(max_nonzero_inc_day + 1):nrow(epicurve)] <- 1
    #epicurve$p[(max_nonzero_inc_day + 1):nrow(epicurve)] <- epicurve$p[max_nonzero_inc_day]/epicurve$incidenceTOT[max_nonzero_inc_day]
    epicurve$p[(max_nonzero_inc_day + 1):nrow(epicurve)] <- 1/mean(epicurve$incidenceTOT[max_nonzero_inc_day - 0:2]/epicurve$p[max_nonzero_inc_day - 0:2])
  }
  return(epicurve)
}

extract_totcurve_Osiris <- function(data, IC = FALSE, report_date = today(), rep_delay = NULL) {
  data <- data %>%
    {if(IC) filter(., NCOVopnameICU == "J") else filter(., NCOVpatZhs == "J")} %>%
    mutate(bron = if("Bron" %in% names(data)) Bron else Bron_land,
           adm_date = if(IC) NCOVopnamedatumICU else NCOVdat1ezkhopn)

  admcurve <- data %>% 
    filter(!is.na(adm_date) & adm_date >= as.Date("2020-02-01") & adm_date <= report_date) %>% # omit symptom onset dates in future
    group_by(adm_date) %>%
    summarise(adm_incidenceABROAD = sum(bron == "Buitenland", na.rm = TRUE),
              adm_incidenceTOT = n()) %>%
    mutate(adm_incidenceNL = adm_incidenceTOT - adm_incidenceABROAD) %>% 
    full_join(tibble(adm_date = seq(min(.$adm_date), as.Date(report_date), by = "day"))) %>%
    replace_na(list(adm_incidenceNL = 0, adm_incidenceTOT = 0, adm_incidenceABROAD = 0)) %>%
    arrange(adm_date) %>%
    select(adm_date, adm_incidenceNL, adm_incidenceABROAD, adm_incidenceTOT)
 
  # choose to ignore whether case is infected abroad or not (11 in total)
  admcurve_SO <- data %>% 
    filter(!is.na(adm_date) & adm_date >= as.Date("2020-02-01") & adm_date <= report_date) %>% # omit symptom onset dates in future
    filter(!is.na(ZIE1eZiekteDt)) %>%
    group_by(adm_date) %>%
    summarise(adm_SO_incidenceTOT = n()) %>%
    full_join(tibble(adm_date = seq(min(.$adm_date), as.Date(report_date), by = "day"))) %>%
    replace_na(list(adm_SO_incidenceTOT = 0)) %>%
    arrange(adm_date) %>%
    select(adm_date, adm_SO_incidenceTOT)
  
  admcurve <- full_join(admcurve, admcurve_SO)
  
  if(length(rep_delay) == 0) {
    # default SO to reporting delay over last 7 days
    reporting_delay <- get_reportingdelay_Osiris(var_name = if(IC) "NCOVopnamedatumICU" else "NCOVdat1ezkhopn", type = if(IC) "IC" else "hosp", start_date = as.Date(report_date) - 7, end_date = as.Date(report_date))
  } else {
    reporting_delay = rep_delay
  }

  exp_cases <- reporting_delay$exp_cases
  
  exp_cases <- c(rep(0, max(0,nrow(admcurve)-length(exp_cases))), rev(exp_cases[1:min(length(exp_cases), nrow(admcurve))]))

  admcurve <- admcurve %>% mutate(adm_incidenceEXP = exp_cases)
  
  # if most recent days are zero, give them one case and a p to reflect total number of last 3 non-zero days, also to max_nonzero_inc_day itself
  max_nonzero_inc_day <- max(which(admcurve$adm_incidenceTOT != 0))
  if(max_nonzero_inc_day < nrow(admcurve)) {
    admcurve$adm_incidenceTOT[(max_nonzero_inc_day + 1):nrow(admcurve)] <- 1
    admcurve$adm_incidenceNL[(max_nonzero_inc_day + 1):nrow(admcurve)] <- 1
    admcurve$adm_incidenceEXP[(max_nonzero_inc_day):nrow(admcurve)] <- mean(admcurve$adm_incidenceTOT[max_nonzero_inc_day - 1:3] + admcurve$adm_incidenceEXP[max_nonzero_inc_day - 1:3]) - admcurve$adm_incidenceTOT[(max_nonzero_inc_day):nrow(admcurve)]
  }

    #reporting_delay <- cumsum(reporting_delay)/sum(reporting_delay)
  
  # if(nrow(admcurve) > length(reporting_delay)) {
  #   reporting_delay <- c(rep(1, nrow(admcurve)-length(reporting_delay)), rev(reporting_delay))
  # } else {
  #   reporting_delay <- rev(reporting_delay[1:nrow(admcurve)])
  # }
  # 
  # admcurve <- admcurve %>% mutate(p = reporting_delay)
  # 
  # # if most recent days are zero, give them one case and a p to reflect total number of last 3 non-zero days
  # max_nonzero_inc_day <- max(which(admcurve$adm_incidenceTOT != 0))
  # if(max_nonzero_inc_day < nrow(admcurve)) {
  #   admcurve$adm_incidenceTOT[(max_nonzero_inc_day + 1):nrow(admcurve)] <- 1
  #   admcurve$adm_incidenceNL[(max_nonzero_inc_day + 1):nrow(admcurve)] <- 1
  #   admcurve$p[(max_nonzero_inc_day + 1):nrow(admcurve)] <- 1/mean(admcurve$adm_incidenceTOT[max_nonzero_inc_day - 0:2]/admcurve$p[max_nonzero_inc_day - 0:2])
  # }
  
  epicurve <- data %>% 
    filter(!is.na(ZIE1eZiekteDt) & ZIE1eZiekteDt <= today()) %>% # omit symptom onset dates in future
    group_by(ZIE1eZiekteDt) %>%
    summarise(incidenceABROAD = sum(bron == "Buitenland", na.rm = TRUE),
              incidenceTOT = n()) %>%
    mutate(incidenceNL = incidenceTOT - incidenceABROAD) %>% 
    full_join(tibble(ZIE1eZiekteDt = seq(min(.$ZIE1eZiekteDt), as.Date(report_date), by = "day"))) %>%
    replace_na(list(incidenceNL = 0, incidenceTOT = 0, incidenceABROAD = 0)) %>%
    arrange(ZIE1eZiekteDt) %>%
    mutate(dates = ZIE1eZiekteDt) %>%
    select(dates, incidenceNL, incidenceABROAD, incidenceTOT)
  
  totcurve <- full_join(admcurve, epicurve, by = c("adm_date" = "dates")) %>%
    rename(dates = adm_date) %>%
    arrange(dates) %>%
    replace_na(list(adm_incidenceNL = 0, adm_incidenceTOT = 0, adm_incidenceABROAD = 0, adm_SO_incidenceTOT = 0, adm_incidenceEXP = 0, incidenceNL = 0, incidenceTOT = 0, incidenceABROAD = 0)) %>%
    mutate(adm_incidenceEXP = round(adm_incidenceEXP))
  
  return(totcurve)
}

calculate_Ru_Osiris <- function(epicurve, SOtoAdm, nsample = 1000) {
  
  #pdfSI <- dnorm(x = 1:7, mean = 6.8/2, sd = 1.27)
  pdfSI <- dgamma(x = 1:10, scale = 1.0, shape = 4)
  pdfSI <- pdfSI/sum(pdfSI)
  g <- pdfSI

  # sample non-observed cases due to reporting delay from neg binomial
  incidence_adm <- matrix(0, nrow = nrow(epicurve), ncol = nsample)
  for(j in 1:nrow(epicurve)) {
    if(epicurve$adm_incidenceTOT[j] != 0) {
      incidence_adm[j, ] <- rnbinom(n = nsample, size = epicurve$adm_incidenceTOT[j], prob = epicurve$adm_incidenceTOT[j]/(epicurve$adm_incidenceTOT[j]+epicurve$adm_incidenceEXP[j])) + epicurve$adm_incidenceTOT[j]
    }
  }
  
  # known SO make up fixed part of incidence matrix
  incidence <- matrix(rep(epicurve$incidenceTOT, nsample), ncol = nsample)

  # sampling symtpom onsets
  for(j in 1:nrow(incidence)) {
    for(i in 1:nsample) {
      tmp <- table(j - sample(0:(length(SOtoAdm)-1), size = incidence_adm[j,i] - epicurve$adm_SO_incidenceTOT[j] , prob = SOtoAdm, replace = TRUE))
      tmp <- tmp[names(tmp) > 0]
      if(length(tmp) > 0) incidence[as.integer(names(tmp)), i] <- incidence[as.integer(names(tmp)), i] + tmp
    }
  }
  # en dan nog eind uppen met sum(SOtoAdm2)
  for(j in (nrow(incidence)-length(SOtoAdm)+2):nrow(incidence)) {
    SOtoAdm2 <- SOtoAdm[1:min(length(SOtoAdm), nrow(incidence)-j+1)]
    incidence[j,] <- incidence[j,]/sum(SOtoAdm2)
  }
  
  # construct instantaneous Rt from incidence
  Rt <- matrix(0, nrow = nrow(incidence), ncol = nsample)
  for(t in 1:nrow(incidence)) {
    g2 <- g[1:min(length(g), t-1)]
    Rt[t,] <- as.numeric((incidence[t, ] - epicurve$incidenceABROAD[t])/colSums(incidence[t - (1:length(g2)), , drop = FALSE]*g2))
  }
  
  # construct case Ru from instantaneous Rt 
  Ru <- matrix(0, nrow = nrow(incidence), ncol = nsample)
  for(u in 1:(nrow(Rt)-1)) {
    g2 <- g[1:min(length(g), nrow(Rt)-u)]
    Ru[u,] <- colSums(g2*Rt[u+1:length(g2), ,drop = FALSE])/sum(g2) # divide by sum(g2) to account for not observed Rt's in future
  }
  
  # construct confidence bounds by sampling from poisson
  Rtupper <- rep(0, nrow(incidence))
  Rtlower <- rep(0, nrow(incidence))
  for(j in 1:nrow(incidence)) {
    tmp <- sapply(1:nsample, function(i) if(incidence[j,i] != 0 & Rt[j,i] !=Inf) rpois(100, incidence[j,i]*Rt[j,i])/incidence[j,i])
    Rtlower[j] <- quantile(unlist(tmp), probs = 0.025)
    Rtupper[j] <- quantile(unlist(tmp), probs = 0.975)
  }
  
  # construct confidence bounds by sampling from poisson
  Ruupper <- rep(0, nrow(incidence))
  Rulower <- rep(0, nrow(incidence))
  for(j in 1:nrow(incidence)) {
    tmp <- sapply(1:nsample, function(i) if(incidence[j,i] != 0 & Ru[j,i] !=Inf) rpois(100, incidence[j,i]*Ru[j,i])/incidence[j,i])
    Rulower[j] <- quantile(unlist(tmp), probs = 0.025)
    Ruupper[j] <- quantile(unlist(tmp), probs = 0.975)
  }
  
  # last Ru's are 0, change to NA to prevent plotting them
  Ru[nrow(incidence),] <- NA
  Rulower[nrow(incidence)] <- NA
  Ruupper[nrow(incidence)] <- NA
  
  epicurve <- epicurve %>%
    mutate(incidence_mean = apply(incidence, MARGIN = 1, mean),
           incidence_lower = apply(incidence, MARGIN = 1, quantile, probs = 0.025),
           incidence_upper = apply(incidence, MARGIN = 1, quantile, probs = 0.975),
           instR = apply(Rt, MARGIN = 1, mean),
           instRlower = Rtlower,
           instRupper = Rtupper,
           caseR = apply(Ru, MARGIN = 1, mean),
           caseRlower = Rulower,
           caseRupper = Ruupper)
  
  return(epicurve)
}

# old
calculate_Ru_Osiris_SO <- function(epicurve, nsample = 1000) {
  
  #pdfSI <- dnorm(x = 1:7, mean = 6.8/2, sd = 1.27)
  pdfSI <- dgamma(x = 1:10, scale = 1.0, shape = 4)
  pdfSI <- pdfSI/sum(pdfSI)
  g <- pdfSI
  
  # sample non-observed cases due to reporting delay from neg binomial
  incidence <- matrix(0, nrow = nrow(epicurve), ncol = nsample)
  for(j in 1:nrow(epicurve)) {
    if(epicurve$incidenceTOT[j] != 0) {
      incidence[j, ] <- rnbinom(n = nsample, size = epicurve$incidenceTOT[j], prob = epicurve$p[j]) + epicurve$incidenceTOT[j]
    }
  }

  # construct instantaneous Rt from incidence
  Rt <- matrix(0, nrow = nrow(epicurve), ncol = nsample)
  for(t in 2:nrow(epicurve)) {
    g2 <- g[1:min(length(g), t-1)]
    Rt[t,] <- as.numeric((incidence[t, ] - epicurve$incidenceABROAD[t])/colSums(incidence[t - (1:length(g2)), , drop = FALSE]*g2))
  }
  
  # construct case Ru from instantaneous Rt 
  Ru <- matrix(0, nrow = nrow(epicurve), ncol = nsample)
  for(u in 1:(nrow(Rt)-1)) {
    g2 <- g[1:min(length(g), nrow(Rt)-u)]
    Ru[u,] <- colSums(g2*Rt[u+1:length(g2), ,drop = FALSE])/sum(g2) # divide by sum(g2) to account for not observed Rt's in future
  }
  
  # construct confidence bounds by sampling from poisson
  Rtupper <- rep(0, nrow(incidence))
  Rtlower <- rep(0, nrow(incidence))
  for(j in 1:nrow(incidence)) {
    tmp <- sapply(1:nsample, function(i) if(incidence[j,i] != 0 & Rt[j,i] !=Inf) rpois(100, incidence[j,i]*Rt[j,i])/incidence[j,i])
    Rtlower[j] <- quantile(unlist(tmp), probs = 0.025)
    Rtupper[j] <- quantile(unlist(tmp), probs = 0.975)
  }

  # construct confidence bounds by sampling from poisson
  Ruupper <- rep(0, nrow(incidence))
  Rulower <- rep(0, nrow(incidence))
  for(j in 1:nrow(incidence)) {
    tmp <- sapply(1:nsample, function(i) if(incidence[j,i] != 0 & Ru[j,i] !=Inf) rpois(100, incidence[j,i]*Ru[j,i])/incidence[j,i])
    Rulower[j] <- quantile(unlist(tmp), probs = 0.025)
    Ruupper[j] <- quantile(unlist(tmp), probs = 0.975)
  }

  # last Ru's are 0, change to NA to prevent plotting them
  Ru[nrow(incidence)] <- NA
  Rulower[nrow(incidence)] <- NA
  Ruupper[nrow(incidence)] <- NA
  
  epicurve <- epicurve %>% mutate(mean_incidence = apply(incidence, MARGIN = 1, mean),
                                  instR = apply(Rt, MARGIN = 1, mean),
                                  caseR = apply(Ru, MARGIN = 1, mean),
                                  instRlower = Rtlower,
                                  instRupper = Rtupper,
                                  caseRlower = Rulower,
                                  caseRupper = Ruupper)
  return(epicurve)
}

# old
plot_epicurve_Osiris <- function(epicurve, IC = FALSE) {
p <- 
  ggplot(data = epicurve %>% mutate(incidenceNL = incidenceTOT - incidenceABROAD,
                                    incidenceExpected = incidenceTOT/p - incidenceTOT) %>%
           pivot_longer(cols = c(incidenceABROAD, incidenceNL, incidenceExpected), names_to = "incidenceType", values_to = "incidence") %>%
           mutate(incidenceType = factor(incidenceType) %>% fct_relevel("incidenceExpected", "incidenceNL", "incidenceABROAD")),
         mapping = aes(x = dates, y = incidence, fill = incidenceType)) + 
  geom_bar(stat = "identity") +
  scale_fill_manual(
    values = c("grey", "#01689b", "#ca005d"),
    labels = c("verwacht op basis van\nrapportagevertraging", "geïnfecteerd in Nederland", "geïnfecteerd buiten Nederland")) +
  scale_x_date(
    limits = c(as.Date("2020-02-27"), NA),
    breaks = as.Date("2020-02-17") + seq(0, 100, 7),
    labels = format(as.Date("2020-02-17") + seq(0, 100, 7), "%d %b")) +
  labs(x = "eerste ziektedag",
       y = "aantal per dag",
       subtitle = paste("gebaseerd op", if(IC) "IC" else "ziekenhuis", "opnames uit OSIRIS data", last(epicurve$dates))) +
  theme_minimal() +
  theme(legend.position = c(0.2, 0.85),
        legend.title = element_blank())

  return(p)
}

plot_totcurve_Osiris <- function(epicurve, IC = FALSE) {
  
  perc_reportedSO <- round(100*sum(epicurve$incidenceTOT)/sum(epicurve$incidence_mean))
  perc_estimatedSO <- 100 - perc_reportedSO
  
  p <- ggplot(data = epicurve,
              mapping = aes(x = dates)) + 
    geom_bar(data = . %>% #mutate(adm_incidenceExpected = adm_incidenceTOT/p - adm_incidenceTOT) %>%
               pivot_longer(cols = c(adm_incidenceABROAD, adm_incidenceNL, adm_incidenceEXP), names_to = "incidenceType", values_to = "incidence") %>%
               mutate(incidenceType = factor(incidenceType) %>% fct_relevel("adm_incidenceEXP", "adm_incidenceNL", "adm_incidenceABROAD")),
             mapping = aes(y = incidence, fill = incidenceType),
             stat = "identity") +
    geom_ribbon(mapping = aes(ymin = incidence_lower, ymax = incidence_upper),
                fill = adjustcolor("#d52b1e", alpha = 0.5)) +
    geom_line(mapping = aes(y = incidence_mean, col = cut(incidence_lower, breaks = c(-1, 1000, 2000)))) +
    scale_fill_manual(
      values = c("grey", "#01689b", "#ffb612"),
      labels = c("verwacht op basis van rapportagevertraging", paste0("opnamedatum ", if(IC) "IC" else "ziekenhuis", ", geïnfecteerd\nin Nederland"), paste0("opnamedatum ", if(IC) "IC" else "ziekenhuis", ", geïnfecteerd\nbuiten Nederland"))) +
    scale_color_manual(
      values = c("#d52b1e", 1),
      labels = c(paste0("eerste ziektedag\n (", perc_reportedSO, "% gerapporteerd, ", perc_estimatedSO, "% geschat)"), "")) +
    scale_x_date(
      limits = c(as.Date("2020-02-14"), last(epicurve$dates)+1),
      breaks = as.Date("2020-02-17") + seq(0, 100, 7),
      labels = format(as.Date("2020-02-17") + seq(0, 100, 7), "%d %b")) +
    labs(x = NULL,
         y = "aantal per dag",
         subtitle = paste("gebaseerd op", if(IC) "IC" else "ziekenhuis", "opnames uit OSIRIS data", last(epicurve$dates))) +
    theme_minimal() +
    theme(legend.position = c(0, 1),
          legend.justification = c(0, 1),
          legend.title = element_blank()) +
    guides(fill = guide_legend(order = 1),
           color = guide_legend(order = 2))
  
  return(p)
}


plot_Reff_Osiris <- function(epicurve, caseReff = TRUE, IC = TRUE) {
  p <- 
    ggplot(data = epicurve,
           mapping = aes(x = dates, 
                         y = if(caseReff) caseR else instR, 
                         ymin = if(caseReff) caseRlower else instRlower, 
                         ymax = if(caseReff) caseRupper else instRupper)) + 
    geom_ribbon(fill = "#c6b8cf", col = NA) +
    geom_line(col = "#42145f") +
    geom_line(y = 1, lty = 2) +
    coord_cartesian(
      ylim = c(0,4)) +
    scale_x_date(
      limits = c(as.Date("2020-02-14"), last(epicurve$dates)+1),
      breaks = as.Date("2020-02-17") + seq(0, 100, 7),
      labels = format(as.Date("2020-02-17") + seq(0, 100, 7), "%d %b")) +
    labs(x = NULL,
         y = if(caseReff) "effectieve R" else "instantane R",
         subtitle = paste("gebaseerd op", if(IC) "IC" else "ziekenhuis", "opnames uit OSIRIS data", last(epicurve$dates))) +
    theme_minimal()
  
  return(p)
}