###########################################################
# Epicurve symptom onset of hospitalised cases in Osiris data
# Rt development over time
# corrected for SO to first report of hospital admission
# from historic Osiris downloads
# 16-04-2020 ___UITZONDERINGSGROND_2___
###########################################################

library(tidyverse)
library(lubridate)
library(cowplot)

source("/___UITZONDERINGSGROND_6___functions_Osiris.r")

# get all case data files from Osiris
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_date <- as.Date("2020-04-16")
report_date <- today()

dataOsiris <- files %>% 
  filter(date == report_date) %>%
  filter(time == max(time)) %>% # get most recent file (if multiple files exist) 
  pull(file_name) %>%
  read_rds()


symptomonset2reporting <- get_reportingdelay_Osiris(var_name = "ZIE1eZiekteDt", type = "hosp", start_date = report_date - 7, end_date = report_date)

epicurve <- extract_epicurve_Osiris(data = dataOsiris, IC = FALSE, report_date = report_date, SOtoRep = symptomonset2reporting)

epicurve <- calculate_Ru_Osiris_SO(epicurve)

plot_epicurve_Osiris(epicurve, IC = FALSE)

# lelijke NA op 21 feb weggehaald door te middelen
epicurve[28, "caseRlower"] <- epicurve[c(27,29), ] %>% pull("caseRlower") %>% mean
epicurve[28, "caseRupper"] <- epicurve[c(27,29), ] %>% pull("caseRupper") %>% mean

plot_Reff_Osiris(epicurve, caseReff = TRUE, IC = FALSE)

plot_grid(plot_epicurve_Osiris(epicurve, IC = FALSE),
          plot_Reff_Osiris(epicurve, caseR = TRUE, IC = FALSE) + labs(subtitle = NULL),
          ncol = 1,
          align = "hv")

ggsave(filename = "/___UITZONDERINGSGROND_6___OLDepicurveReff_Osiris_ZKH_16apr2020.tiff", width = 6, height = 6, dpi = 150)


AdmToRep <- get_reportingdelay_Osiris(var_name = "NCOVdat1ezkhopn", type = "hosp", start_date = report_date - 7, end_date = report_date)

epicurve <- extract_totcurve_Osiris(data = dataOsiris, IC = FALSE, report_date = report_date, rep_delay = AdmToRep)

SOtoAdm <- get_symptomonset2admission_Osiris(IC = FALSE, report_date = report_date)

epicurve <- calculate_Ru_Osiris(epicurve, SOtoAdm = SOtoAdm)

plot_Reff_Osiris(epicurve %>% mutate(caseR = ifelse(dates > report_date - 14, NA, caseR),
                                    caseRlower = ifelse(dates > report_date - 14, NA, caseRlower),
                                    caseRupper = ifelse(dates > report_date - 14, NA, caseRupper)), caseReff = TRUE, IC = FALSE) 


which(cumsum(AdmToRep)/sum(AdmToRep) > 0.8)
qgamma(p = 0.8, scale = 1.0, shape = 4)

plot_grid(plot_totcurve_Osiris(epicurve = epicurve, IC = FALSE),
          plot_Reff_Osiris(epicurve #%>% mutate(caseRlower = ifelse(dates > report_date - 14, NA, caseRlower),
                                     #          caseRupper = ifelse(dates > report_date - 14, NA, caseRupper),
                                      #         caseR = ifelse(dates > report_date - 14, NA, caseR))
                           , caseR = TRUE, IC = FALSE) + labs(subtitle = NULL),
          ncol = 1,
          align = "hv")


ggsave(filename = "/___UITZONDERINGSGROND_6___epicurveReff_Osiris_ZKH_16apr2020_vol.tiff", width = 6, height = 6, dpi = 150)

#ggsave(filename = "/___UITZONDERINGSGROND_6___epicurve_Osiris_ZKH_9apr2020.tiff", dpi = 150)

#ggsave(filename = "/___UITZONDERINGSGROND_6___Reff_Osiris_ZKH_9apr2020.tiff", dpi = 150)



# Per provincie (ook nog delay per provincie)

provincies <- unique(dataOsiris$Provincie)

epicurves <- tibble()
for(provincie in provincies) {
  print(provincie)
  symptomonset2reporting <- get_reportingdelay_Osiris_province(province = provincie, var_name = "ZIE1eZiekteDt", type = "hosp", start_date = report_date - 7, end_date = report_date)
  epicurve <- extract_epicurve(data = dataOsiris %>% filter(Provincie == provincie), IC = FALSE, report_date = report_date, SOtoRep = symptomonset2reporting)
  epicurve <- calculate_Ru_Osiris(epicurve)
  epicurve <- epicurve %>% mutate(provincie = provincie)
  epicurves <- bind_rows(epicurves, epicurve)
}

plot_epicurve_Osiris(epicurve = epicurves) + facet_wrap(facets = vars(provincie))

plot_Reff(epicurve = epicurves, source = "OSIRIS", IC = FALSE) + facet_wrap(facets = vars(provincie))