```{r, results='hide',warning=FALSE, message=FALSE}
library(data.table)
library(dplyr)
library(lubridate)
library(ggplot2)
library(viridis)
library(reshape2)
library(knitr)
library(kableExtra)
```

**New version: Uses DTs from Jeremy**


#### Load conditional severity rates at county-level for CA counties.
Get expected proportion of hospital census in ICU by calculating LOS-weighted expected proportion of hospitalizations in ICU based on county demographics and hospitalization severity by age from Imperial report (calculated earlier in separate code).

```{r}
# Constants from state request https://tinyurl.com/v4u5vj6
ICU_LOS = 9
ACUTE_LOS = 7

df <- read.csv('../data/county_age_severity_rates_v6.csv', stringsAsFactors = FALSE)
df$County <- gsub('city', 'City', df$County)
df = df %>% filter(State == "California") %>% group_by(FIPS, County) %>%
    summarise(
      total_population = sum(population_in_age_group),
      wtd_critical_case_rate =  weighted.mean(critical_case_rate, population_in_age_group),
      wtd_acute_case_rate =  weighted.mean(severe_cases_rate, population_in_age_group)) %>%
  mutate(
      wtd_prop_hosp_in_icu = ICU_LOS*wtd_critical_case_rate/(ICU_LOS*wtd_critical_case_rate+ACUTE_LOS*wtd_acute_case_rate))

df = df %>% ungroup() %>% select(FIPS, County, wtd_prop_hosp_in_icu)
df = df %>% mutate(County = gsub(' County', '', County))
df[1:10,] %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```

#### Join in bed (Acute, ICU) data
```{r}
acute_beds_dt = fread('../data/acute_byFIPS.csv')
icu_beds_dt = fread('../data/icu_byFIPS.csv')
bed_dt = merge(acute_beds_dt, icu_beds_dt, by = "FIPS")
df_w_beds = merge(df, bed_dt, by="FIPS", all.x = TRUE)
df_w_beds[1:10,] %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```


#### Merge in current hospitalizations and DTs from Jeremy (4/4)

* Assume that hospitalizations numbers provided are current census. 
* **NOTE: I do not use the reported ICU Census by County because it is not clear whether the reported hospital census includes ICU census. Their website says it does, but some counties (e.g. Shasta) have a higher ICU census than overall census.**
* I break overall hospital census into ICU and Acute given the expected breakdown based on the county's demographics and the conditional ICU rates by age estimated by Imperial.
* I then calculate the current estimated cumulative Floor and ICU hospitalizations given a past doubling time (using the short and long DTs estimated by Jeremy) and LOS values provided by the state.

```{r}
# CA Hospitalizations on 3/30
hosp = fread('data/hosp-est-dts-2020-04-04.csv')[,.(County, 
                                                    current_census = Hospitalized_4_4,
                                                    dt_short = bound_double_short,
                                                    dt_long = bound_double_long)]
df_w_beds_hosp = as.data.table(merge(df_w_beds, hosp, by = 'County'))

# break census into icu and floor
df_w_beds_hosp[,`:=`(
  est_current_icu_census = wtd_prop_hosp_in_icu * current_census,
  est_current_floor_census = (1-wtd_prop_hosp_in_icu) * current_census
)]

# estimate cumulative floor/icu hospitalizations with los and dt

df_w_beds_hosp[,`:=`(
  est_current_cum_icu_hosp_short = 
    est_current_icu_census/(1 - 2^(-1*ICU_LOS/dt_short)),
  est_current_cum_icu_hosp_long = 
    est_current_icu_census/(1 - 2^(-1*ICU_LOS/dt_long)),
  
  est_current_cum_floor_hosp_short = 
    est_current_floor_census/(1 - 2^(-1*ACUTE_LOS/dt_short)),
  est_current_cum_floor_hosp_long = 
    est_current_floor_census/(1 - 2^(-1*ACUTE_LOS/dt_long))
)]

df_w_beds_hosp %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```

For now, only project hospitalizations for counties with nonzero hospitalizations.
Projecting hospitalizations for other counties requires stronger assumptions about testing and/or fatality rates in order to estimate hospitalizations.

#### Project 14-day time series of cumulative hospitalizations (Floor and ICU) with lower and upper doubling time estimates from Jeremy.

```{r}
start_date = '2020-04-04'
num_days_to_project = 14

dates_dt = data.table(date = as.Date(start_date) + seq(0, num_days_to_project))

long_dt = merge(
  dates_dt[, i:=1],
  df_w_beds_hosp[,i:=1],
  by = 'i',
  allow.cartesian = TRUE
  )[,i:=NULL]

long_dt[, days_since_start_date := as.numeric(date - as.Date(start_date)), by=.I]

long_dt[, `:=`(
  est_cum_floor_hosp_short = est_current_cum_floor_hosp_short * 2^(days_since_start_date/dt_short),
  est_cum_floor_hosp_long = est_current_cum_floor_hosp_long * 2^(days_since_start_date/dt_long),
  
  est_cum_icu_hosp_short = est_current_cum_icu_hosp_short * 2^(days_since_start_date/dt_short),
  est_cum_icu_hosp_long = est_current_cum_icu_hosp_long * 2^(days_since_start_date/dt_long)
)]
long_dt[County == "Santa Clara", .(County, date, est_cum_floor_hosp_short, est_cum_floor_hosp_long, est_cum_icu_hosp_short, est_cum_icu_hosp_long)] %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```


#### Estimate hospital census using LOS values from State
Main idea: `Census[t] <- CumulativeHospitalizations[t] - CumulativeHospitalizations[t-LOS]`

```{r}

# function from Shiny app
get_hospitalizations = function(cumulative_cases, los, doubling_time) {
  
  days_to_hospitalization = 0 # not modeling infections so no need to account for this 
  
  # project back los + days to hospitalization days
  back_vec = c(rep(NA, los + days_to_hospitalization), cumulative_cases)
  for (i in (los + days_to_hospitalization):1) {
    back_vec[i] = back_vec[i + 1]/2^(1/doubling_time)
  }
  
  # get indices of original vectors
  original_start = los + days_to_hospitalization + 1
  original_end = los + days_to_hospitalization + length(cumulative_cases)
  stopifnot(all.equal(back_vec[original_start:original_end], cumulative_cases))
  stopifnot(length(back_vec) == original_end)
  
  # get indices of vectors shifted by days to hospitalization
  shifted_start = original_start - days_to_hospitalization
  shifted_end = original_end - days_to_hospitalization
  
  # subtract off for length of stay
  return_vec = back_vec[shifted_start:shifted_end] - back_vec[(shifted_start - los):(shifted_end - los)]
  
  return(list(result = return_vec, back_vec = back_vec[1:los]))
}


# use on our data
long_dt[, est_icu_census_short := get_hospitalizations(est_cum_icu_hosp_short, ICU_LOS, doubling_time=min(dt_short)), by = c('County')]
long_dt[, est_icu_census_long := get_hospitalizations(est_cum_icu_hosp_long, ICU_LOS, doubling_time=min(dt_long)), by = c('County')]

long_dt[, est_floor_census_short := get_hospitalizations(est_cum_floor_hosp_short, ACUTE_LOS, doubling_time=min(dt_short)), by = c('County')]
long_dt[, est_floor_census_long := get_hospitalizations(est_cum_floor_hosp_long, ACUTE_LOS, doubling_time=min(dt_long)), by = c('County')]

long_dt[County == "Santa Clara", .(County, date, est_icu_census_long, est_icu_census_short, est_floor_census_long, est_floor_census_short)] %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```


Example plot of ICU demand
```{r}
ggplot(melt(long_dt[County %in% c('Alameda', 'Santa Clara', 'San Mateo', 'San Francisco'), 
                    .(County, date, est_icu_census_long, est_icu_census_short)],
            id.vars = c('County', 'date')),
       aes(x=date, y=value, color=factor(variable))) + geom_line() +
       facet_wrap(~County, ncol=2, scales="free")
```


Create wide views with time series by county

First, ICU time series by county
```{r}
PROP_OF_AHA_BEDS_AVAILABLE = 0.3

icu_time_series_long = long_dt %>% group_by(County, date) %>%
  summarise(
    est_icu_beds_available = max(round(num_icu_beds * PROP_OF_AHA_BEDS_AVAILABLE)),
    icu_census_lb = round(est_icu_census_long),
    icu_census_m = round((est_icu_census_long+est_icu_census_short)/2),
    icu_census_ub = round(est_icu_census_short),
  ) %>% ungroup()

icu_time_series_long = icu_time_series_long %>% mutate(icu_census_str = paste0(icu_census_m, ' [', icu_census_lb, ', ', icu_census_ub, ']'))

# icu_time_series_long %>% filter(County == "Santa Clara")

wide_icu_str = as.data.table(dcast(icu_time_series_long, County + est_icu_beds_available ~ date, value.var = c('icu_census_str')))
wide_icu_str %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```

Now create a version that just shows excess census over 30% of AHA ICU beds
```{r}
icu_time_series_long_f = icu_time_series_long %>%
  mutate(excess_demand = paste0(pmax(0, icu_census_m - est_icu_beds_available), 
                                ' [', pmax(0, icu_census_lb - est_icu_beds_available), ', ', 
                                icu_census_ub - est_icu_beds_available, ']')) %>%
  mutate(excess_demand = ifelse(icu_census_ub>est_icu_beds_available, excess_demand, NA))

wide_icu_excess_demand = as.data.table(dcast(icu_time_series_long_f, County + est_icu_beds_available ~ date, value.var = c('excess_demand')))
wide_icu_excess_demand %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```

Write the ICU tables to CSV for now
```{r}
fwrite(wide_icu_str, 'data/wide_icu_str_v5.csv')
fwrite(wide_icu_excess_demand, 'data/wide_icu_excess_demand_v5.csv')
```


And let's do the same for Acute/Floor beds

```{r}
acute_time_series_long = long_dt %>% group_by(County, date) %>%
  summarise(
    est_floor_beds_available = max(round(num_acute_beds * PROP_OF_AHA_BEDS_AVAILABLE)),
    acute_census_lb = round(est_floor_census_long),
    acute_census_m = round((est_floor_census_long + est_floor_census_short)/2),
    acute_census_ub = round(est_floor_census_short),
  ) %>% ungroup()

acute_time_series_long = acute_time_series_long %>% mutate(acute_census_str = paste0(acute_census_m, ' [', acute_census_lb, ', ', acute_census_ub, ']'))

# acute_time_series_long %>% filter(County == "Santa Clara")

wide_acute_str = as.data.table(dcast(acute_time_series_long, County + est_floor_beds_available ~ date, value.var = c('acute_census_str')))
wide_acute_str %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```

```{r}
acute_time_series_long_f = acute_time_series_long %>%
  mutate(excess_demand = paste0(pmax(0, acute_census_m - est_floor_beds_available), 
                                ' [', pmax(0, acute_census_lb - est_floor_beds_available), ', ', 
                                acute_census_ub - est_floor_beds_available, ']')) %>%
  mutate(excess_demand = ifelse(acute_census_ub > est_floor_beds_available, excess_demand, NA))

wide_acute_excess_demand = as.data.table(dcast(acute_time_series_long_f, County + est_floor_beds_available ~ date, value.var = c('excess_demand')))
wide_acute_excess_demand %>% kable() %>% kable_styling(bootstrap_options = c("striped", "hover"))
```


```{r}
fwrite(wide_acute_str, 'data/wide_acute_str_v5.csv')
fwrite(wide_acute_excess_demand, 'data/wide_acute_excess_demand_v5.csv')
```

