#################################################################################
# Extract epicurve of reported admissions of hospitalised patients from case data
# - with nowcast to account for expected admissions
# - distinction between infected abroad or not
# - distinction between patients with known and unknown symptom onset
#################################################################################

extract_epicurve <- function(data, date, start_date = NULL, nowcst = TRUE) {
  
  data <- data %>% filter(admission_date_reported <= date)
  
  nowcast_data <- if(nowcst) {
    nowcast(
      case.data = data %>%
        filter(admission_date <= admission_date_reported) %>%
        transmute(admission.date = admission_date, report.date = admission_date_reported) %>% 
        drop_na,
      start.date = start_date,
      nowcast.date = date)
  } else {
    tibble(admission.date = data %>% pull(admission_date) %>% unique, nowcast_p50 = NA, trend_med = NA)
  }
  
  # note: symptom onset is used in definition for infected_abroad (Bron_buitenland mogelijk of waarschijnlijk) -> all infected abroad have symptoms and symptom onset
  admission_curve <- data %>% 
    filter(admission_date <= date & !is.na(admission_date)) %>%
    group_by(admission_date) %>% 
    summarise(admission_infectedabroad = sum(infected_abroad == "Yes"),
              admission_knownSO = sum(!is.na(symptom_onset)),
              admission_nosymptoms = sum(no_symptoms == 1),
              admission_knownSO_infectedabroad = sum(!is.na(symptom_onset) & infected_abroad == "Yes"),
              admission_nosymptoms_infectedabroad = sum(no_symptoms == 1 & infected_abroad == "Yes"),
              admission = n()) %>% 
    full_join(nowcast_data %>% rename(admission_expected = nowcast_p50), by = c("admission_date" = "admission.date")) %>%
    full_join(tibble(admission_date = seq(min(.$admission_date), date, by = "day"))) %>%
    arrange(admission_date) %>%
    replace_na(replace = list(admission = 0, admission_infectedabroad = 0, admission_knownSO = 0, admission_nosymptoms = 0, admission_knownSO_infectedabroad = 0, admission_nosymptoms_infectedabroad = 0, admission_expected = 0)) %>%
    select(admission_date, admission, admission_knownSO, admission_nosymptoms, admission_infectedabroad, admission_knownSO_infectedabroad, admission_nosymptoms_infectedabroad, admission_expected, trend_med) %>%
    mutate(admission_expected = ifelse(admission > admission_expected, admission, admission_expected)) # should not be necessary
  
  # check whether nowcast overshoots (when twice higher than expected from reporting delay)
  # interval <- data %>% mutate(interval = as.integer(admission_date_reported - admission_date)) %>% filter(interval >= 0) %>% pull(interval)
  # interval_cdf <- ecdf(interval)(0:(nrow(admission_curve)-1)) 
  # admission_curve <- admission_curve %>% mutate(admission_check = round(admission/rev(interval_cdf)))
  # 
  # if(any(admission_curve$admission_expected > 2*admission_curve$admission_check)) {
  #   print("correct nowcast overshooting")
  #   admission_curve <- admission_curve %>% mutate(admission_expected_nc = admission_expected,
  #                                                 admission_expected = ifelse(admission_expected > 2*admission_check, admission_check, admission_expected))
  # }
  
  epicurve <- data %>% 
    filter(symptom_onset <= date & !is.na(symptom_onset)) %>%
    group_by(symptom_onset) %>%
    summarise(incidence_knownSO_infectedabroad = sum(infected_abroad == "Yes"),
              incidence_knownSO = n()) %>%
    select(symptom_onset, incidence_knownSO, incidence_knownSO_infectedabroad)
  
  epicurve <- full_join(admission_curve, epicurve, by = c("admission_date" = "symptom_onset")) %>%
    rename(dates = admission_date) %>%
    arrange(dates) %>%
    replace_na(list(admission = 0, admission_infectedabroad = 0, admission_knownSO = 0, admission_nosymptoms = 0, admission_knownSO_infectedabroad = 0, admission_expected = 0, incidence_knownSO = 0, incidence_knownSO_infectedabroad = 0, admission_nosymptoms_infectedabroad = 0))
  
  return(epicurve)
  
}