###########################################################
# 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_nowcast.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-05-15")
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()

nowcast_adm <- get_nowcast_Osiris(var_name = "NCOVdat1ezkhopn", type = "hosp", end_date = report_date)

epicurve <- extract_totcurve_Osiris(data = dataOsiris, IC = FALSE, report_date = report_date, now_cast = nowcast_adm)

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

epicurve <- calculate_Ru_Osiris(epicurve, SOtoAdm = SOtoAdm)

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

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

ggsave(filename = paste0("/___UITZONDERINGSGROND_6___epicurveReff_Osiris_ZKH_", as.Date(report_date), ".png"), width = 7, height = 6, dpi = 300)

# Per provincie

provincies <- unique(dataOsiris$Provincie)

epicurves <- tibble()
for(provincie in provincies[4:12]) {
  print(provincie)
  AdmNowcast <- get_nowcast_Osiris_province(province = provincie, var_name = "NCOVdat1ezkhopn", type = "hosp", start_date = as.Date("2020-03-13"), end_date = report_date)
  epicurve <- extract_totcurve_Osiris(data = dataOsiris %>% filter(Provincie == provincie), IC = FALSE, report_date = report_date, now_cast = AdmNowcast)
  SOtoAdm <- get_symptomonset2admission_Osiris(IC = FALSE, report_date = report_date)
  epicurve <- calculate_Ru_Osiris(epicurve, SOtoAdm = SOtoAdm)
  epicurve <- epicurve %>% mutate(provincie = provincie)
  epicurves <- bind_rows(epicurves, epicurve)
}


plot_totcurve_Osiris(epicurve = epicurves) + facet_wrap(facets = vars(provincie)) + guides(fill = FALSE, color = FALSE)

ggsave(filename = paste0("/___UITZONDERINGSGROND_6___epicurve_Osiris_ZKH_provinces", as.Date(report_date), ".png"), width = 15, height = 9, dpi = 300)


plot_Reff_Osiris(epicurve = epicurves  %>% mutate(caseR = ifelse(dates > report_date - 14, NA, caseR)), IC = FALSE) + facet_wrap(facets = vars(provincie))

ggsave(filename = paste0("/___UITZONDERINGSGROND_6___Reff_Osiris_ZKH_provinces", as.Date(report_date), ".png"), width = 15, height = 9, dpi = 300)