##########################################################################
# Script for Rt estimation from Osiris hospital data
# - from case data (symptom onset, admission date, admission reported, infected abroad)
# - extract epicurve (admissions adjusted with nowcast for expected admissions)
# - calculate instantaneous R and case R (with 95% confidence interval)
# - save epicurve
# - save plot admission curve, epicurve and case R for Jaap
# - export json file for dashboard
##########################################################################

library(jsonlite)
library(tidyverse)
library(lubridate)
library(cowplot)
library(odbc)
library(runjags)
library(coda)


# Source functions
list.files(path = "./___UITZONDERINGSGROND_6___functions", full.names = TRUE) %>% walk(.f = source)


analysis_date <- today()
analysis_date <- as.Date("2020-03-30")

# Generate / retrieve case data
case_data <- generate_case_data(date = analysis_date)

# Extract epicurve: expected admissions by nowcast, distinction between those infected abroad or not, and between known and unknown symptom onset
epicurve <- extract_epicurve(data = case_data %>% filter(admission_date_reported > min(case_data$admission_date_reported) & admission_date >= as.Date("2020-02-27")), date = analysis_date, start_date = analysis_date - 100)

# distribution symptom onset to admission
SO2adm <- determine_SO2admission(data  = case_data, date = analysis_date, start_date = analysis_date - 100)

# calculate reproduction numbers (instantaneous R and case R)
tmp <- calculate_R(epicurve = epicurve, SO2adm = SO2adm)
epicurve <- tmp$epicurve
incidence <- tmp$incidence
saveRDS(object = epicurve, file = paste0("/___UITZONDERINGSGROND_6___epicurve_", as.Date(analysis_date), ".rds"))
#saveRDS(object = incidence, file = paste0("/___UITZONDERINGSGROND_6___incidence_", as.Date(analysis_date), ".rds"))


# plot epicurve and reproduction number
plot_grid(plot_totcurve_Osiris(epicurve = epicurve),
          plot_Reff_Osiris(epicurve %>% mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR)), caseR = TRUE) + labs(subtitle = NULL),
          ncol = 1,
          align = "hv")
ggsave(filename = paste0("/___UITZONDERINGSGROND_6___epicurveReff_Osiris_ZKH_", as.Date(analysis_date), ".png"), width = 7, height = 6, dpi = 300)


epicurve %>% filter(dates == analysis_date - 14) %>% select(dates, caseR, caseRlower, caseRupper)

# png aan ___UITZONDERINGSGROND_2___ mailen; hij maakt presentatie voor OMT en/of Jaap

# json file for dashboard
# reproduction.json <- epicurve %>% 
#   select(dates, caseRlower, caseR, caseRupper) %>% 
#   filter(dates >= as.Date("2020-02-17")) %>%
#   mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR),
#          caseR = round(caseR, digits = 2),
#          caseRlower = round(caseRlower, digits = 2),
#          caseRupper = round(caseRupper, digits = 2)) %>%
#   rename(Date = dates,
#          Rt_low = caseRlower,
#          Rt_avg = caseR,
#          Rt_up = caseRupper) %>%
#   toJSON 

# write(x = reproduction.json, file = paste0("/___UITZONDERINGSGROND_6___COVID-19_reproductiegetal.json"))
# write(x = reproduction.json, file = paste0("/___UITZONDERINGSGROND_6___COVID-19_reproductiegetal_", as.Date(analysis_date), ".json"))

# eerste file (COVID-19_Reproductiegetal.json) op maandagmiddag (na 15:00) of dinsdagochtend (voor 13:00) op R:\Projecten\COVID-19\Surveillance\Open_data zetten


epicurve_ZKH <- readRDS(file = paste0("/___UITZONDERINGSGROND_6___epicurve_", as.Date(analysis_date), ".rds"))
epicurve_reports <- readRDS(file = paste0("/___UITZONDERINGSGROND_6___epicurve_reports_", as.Date(analysis_date), ".rds"))

ggplot(data = epicurve_ZKH %>% filter(dates >= as.Date("2020-07-01")) %>% mutate(caseR = ifelse(dates >= analysis_date - 14, NA, caseR)),
       mapping = aes(x = dates, 
                     y = caseR , 
                     ymin = caseRlower , 
                     ymax = caseRupper)) + 
  geom_ribbon(fill = adjustcolor("#007bc7", alpha = 0.3), col = NA) +
  geom_line(col = "#007bc7") +
  geom_ribbon(data = epicurve_reports %>% filter(dates >= as.Date("2020-07-01")) %>% mutate(caseR = ifelse(dates >= analysis_date - 14, NA, caseR)),
              mapping = aes(x = dates, 
                            y = caseR , 
                            ymin = caseRlower , 
                            ymax = caseRupper),
              fill = adjustcolor("#42145f", alpha = 0.3), col = NA) +
  geom_line(data = epicurve_reports %>% filter(dates >= as.Date("2020-07-01")) %>% mutate(caseR = ifelse(dates >= analysis_date - 14, NA, caseR)),
            mapping = aes(x = dates, 
                          y = caseR , 
                          ymin = caseRlower , 
                          ymax = caseRupper),
            col = "#42145f") +
  geom_line(y = 1, lty = 2) +
  coord_cartesian(
    ylim = c(0, 3)) +
  scale_x_date(
    limits = c(as.Date("2020-07-01"), last(epicurve$dates)+1),
    expand = expansion(add = -0.5),
    breaks = "4 days",
    minor_breaks = "1 day",
    date_labels = "%e %b") +
  labs(x = NULL,
       y = "effectieve R",
       subtitle = paste("gebaseerd op ziekenhuisopnames (blauw) en meldingen (paars) uit OSIRIS data", last(epicurve$dates))) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
        panel.grid.minor = element_blank())


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