##########################################################################
# Script for Rt estimation from Osiris hospital data and reports
# - combine epicurve ZKH Osiris (2020-06-26) with most recent epicurve reports Osiris
# - 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 <- as.Date("2020-07-03")
analysis_date <- today()

# import most recent epicurve reports and epicurve_ZKH of June 26 (last reported Reff_ZKH on June 12)
epicurve_reports <- readRDS(file = paste0("/___UITZONDERINGSGROND_6___epicurve_reports_", analysis_date, ".rds"))
epicurve_ZKH <- readRDS(file = paste0("/___UITZONDERINGSGROND_6___epicurve_", as.Date("2020-06-26"), ".rds"))
#epicurve_ZKH <- readRDS(file = paste0("/___UITZONDERINGSGROND_6___epicurve_", analysis_date, ".rds"))

# transition from ZKH to reports on June 12
epicurve <- bind_rows("ziekenhuis" = epicurve_ZKH %>% filter(dates <= as.Date("2020-06-12")),
                      "meldingen" = epicurve_reports %>% filter(dates > as.Date("2020-06-12")),
                      .id = "source")


# plot effective R with cut on June 12
plot_Reff_Osiris(epicurve %>% mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR)), caseR = TRUE) +
  labs(subtitle = paste("OSIRIS data", last(epicurve$dates))) +
  geom_vline(xintercept = as.Date("2020-06-12") + 0.5, lty = 2) +
  annotate("text", x = as.Date("2020-06-10"), y = 3.5, hjust = 1, label = "gebaseerd op\nziekenhuisopnames") +
  annotate("text", x = as.Date("2020-06-14"), y = 3.5, hjust = 0, label = "gebaseerd op\nmeldingen")


plot_totcurve_Osiris2(epicurve) +
  geom_vline(xintercept = as.Date("2020-06-12") + 0.5, lty = 2)


# plot epicurve and reproduction number
plot_grid(plot_totcurve_Osiris2(epicurve) +
            geom_vline(xintercept = as.Date("2020-06-12") + 0.5, lty = 2),
          plot_Reff_Osiris(epicurve %>% mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR)), caseR = TRUE) +
            labs(subtitle = paste("OSIRIS data", last(epicurve$dates))) +
            geom_vline(xintercept = as.Date("2020-06-12") + 0.5, lty = 2) +
            annotate("text", x = as.Date("2020-06-10"), y = 3.5, size = 3, hjust = 1, label = "gebaseerd op\nziekenhuisopnames") +
            annotate("text", x = as.Date("2020-06-14"), y = 3.5, size = 3, hjust = 0, label = "gebaseerd op\nmeldingen") + 
            labs(subtitle = NULL),
          ncol = 1,
          align = "hv") 

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

# plot epicurve and reproduction number
plot_grid(plot_totcurve_Osiris3(epicurve),
          plot_Reff_Osiris(epicurve %>% mutate(caseR = ifelse(dates > analysis_date - 14, NA, caseR)), caseR = TRUE) +
            labs(subtitle = paste("OSIRIS data", last(epicurve$dates))) +
            coord_cartesian(
              ylim = c(0,2)) +
            scale_x_date(
              limits = c(as.Date("2020-06-12"), last(epicurve$dates)+1),
              expand = expansion(add = -0.5),
              breaks = "4 days",
              minor_breaks = "1 day",
              date_labels = "%e %b") +
            labs(subtitle = NULL),
          ncol = 1,
          align = "hv") 

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

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


# json file for dashboard
reproduction.json <- epicurve %>% 
  select(dates, caseRlower, caseR, caseRupper, source) %>% 
  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),
         source = ifelse(source == "ziekenhuis", "hosp", "testpos")) %>%
  rename(Date = dates,
         Rt_low = caseRlower,
         Rt_avg = caseR,
         Rt_up = caseRupper,
         population = source) %>%
  toJSON 

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