fix+refactor: cleanup pipelines (#169)

* fix covid and flu explore to functioning state
* move shared utilities to new directory
* switch to new convention of marking targets global variables with g_
* share convention with prod
* move prod to using params_grid idiom
* lots of small tweaks to utilities and bug fixes
* use_crew = yes default
* simplify targets options
* remove flu_hosp_tiny
* update notebooks

Author: dshemetov

Makefile

@@ -13,12 +13,12 @@ run:
 	Rscript scripts/run.R
 
 prod-covid:
-	export TAR_RUN_PROJECT=covid_hosp_prod; \
-	Rscript scripts/run_prod.R
+	export TAR_PROJECT=covid_hosp_prod; \
+	Rscript scripts/run.R
 
 prod-flu:
-	export TAR_RUN_PROJECT=flu_hosp_prod; \
-	Rscript scripts/run_prod.R
+	export TAR_PROJECT=flu_hosp_prod; \
+	Rscript scripts/run.R
 
 prod: prod-covid prod-flu update_site netlify
 

R/aux_data_utils.R

@@ -176,7 +176,17 @@ gen_pop_and_density_data <-
       )
   }
 
-daily_to_weekly <- function(epi_df, agg_method = c("sum", "mean"), day_of_week = 4L, day_of_week_end = 7L, keys = "geo_value", values = c("value")) {
+#' Aggregate a daily archive to a weekly archive.
+#'
+#' By default, aggregates from Sunday to Saturday and labels with the Wednesday
+#' of that week.
+#'
+#' @param epi_df the archive to aggregate.
+#' @param agg_method the method to use to aggregate the data, one of "sum" or "mean".
+#' @param keys the columns to group by.
+#' @param values the columns to aggregate.
+daily_to_weekly <- function(epi_df, agg_method = c("sum", "mean"), keys = "geo_value", values = c("value")) {
+  agg_method <- arg_match(agg_method)
   epi_df %>%
     mutate(epiweek = epiweek(time_value), year = epiyear(time_value)) %>%
     group_by(across(any_of(c(keys, "epiweek", "year")))) %>%
@@ -299,6 +309,13 @@ drop_non_seasons <- function(epi_data, min_window = 12) {
     )
 }
 
+get_nwss_coarse_data <- function(disease = c("covid", "flu")) {
+  disease <- arg_match(disease)
+  aws.s3::get_bucket_df(prefix = glue::glue("exploration/aux_data/nwss_{disease}_data"), bucket = "forecasting-team-data") %>%
+    slice_max(LastModified) %>%
+    pull(Key) %>%
+    aws.s3::s3read_using(FUN = readr::read_csv, object = ., bucket = "forecasting-team-data")
+}
 
 #' add a column summing the values in the hhs region
 #' @param hhs_region_table the region table

R/forecasters/data_transforms.R

@@ -234,8 +234,7 @@ data_whitening <- function(epi_data, colname, learned_params, nonlin_method = c(
     join_cols <- key_colnames(epi_data, exclude = "time_value")
   }
   nonlin_method <- arg_match(nonlin_method)
-  res <- epi_data %>%
-    left_join(learned_params, by = join_cols)
+  res <- epi_data %>% left_join(learned_params, by = join_cols)
   if (nonlin_method == "quart_root") {
     res %<>% mutate(across(all_of(colname), ~ (.x + 0.01)^(1 / 4)))
   }

R/forecasters/data_validation.R

@@ -69,7 +69,7 @@ confirm_sufficient_data <- function(epi_data, ahead, args_input, outcome, extra_
   # TODO: Buffer should probably be 2 * n(lags) * n(predictors). But honestly,
   # this needs to be fixed in epipredict itself, see
   # https://github.com/cmu-delphi/epipredict/issues/106.
-  if (extra_sources == c("")) {
+  if (identical(extra_sources, "")) {
     extra_sources <- character(0L)
   }
   has_no_last_nas <- epi_data %>%

R/forecasters/ensemble_linear_climate.R

@@ -14,7 +14,7 @@
 #' @param other_weights if non null, it should be a tibble giving a list of weights by forecaster and geo_value
 #' @importFrom rlang %||%
 #' @export
-ensemble_linear_climate <- function(forecasts,
+ensemble_climate_linear <- function(forecasts,
                                     aheads,
                                     other_weights = NULL,
                                     probs = covidhub_probs(),

R/forecasters/forecaster_climatological.R

@@ -59,7 +59,7 @@ climate_linear_ensembled <- function(epi_data,
   pred_geo_climate <- climatological_model(epi_data, ahead, geo_agg = FALSE) %>% mutate(forecaster = "climate_geo")
   pred_linear <- forecaster_baseline_linear(epi_data, ahead, residual_tail = residual_tail, residual_center = residual_center) %>% mutate(forecaster = "linear")
   pred <- bind_rows(pred_climate, pred_linear, pred_geo_climate) %>%
-    ensemble_linear_climate((args_list$aheads[[1]]) / 7) %>%
+    ensemble_climate_linear((args_list$aheads[[1]]) / 7) %>%
     ungroup()
   # undo whitening
   pred_final <- pred %>%

R/forecasters/forecaster_no_recent_outcome.R

@@ -62,13 +62,13 @@ no_recent_outcome <- function(epi_data,
   args_input[["quantile_levels"]] <- quantile_levels
   args_list <- do.call(default_args_list, args_input)
   # if you want to hardcode particular predictors in a particular forecaster
-  predictors <- c(outcome, extra_sources[[1]])
-  c(args_list, tmp_pred, trainer) %<-% sanitize_args_predictors_trainer(epi_data, outcome, predictors, trainer, args_list)
-  if (extra_sources[[1]] == "") {
+  if (identical(extra_sources[[1]], "")) {
     predictors <- character()
   } else {
     predictors <- extra_sources[[1]]
   }
+  c(args_list, tmp_pred, trainer) %<-% sanitize_args_predictors_trainer(epi_data, outcome, predictors, trainer, args_list)
+
   # end of the copypasta
   # finally, any other pre-processing (e.g. smoothing) that isn't performed by
   # epipredict

R/forecasters/forecaster_scaled_pop_seasonal.R

@@ -59,6 +59,7 @@ scaled_pop_seasonal <- function(epi_data,
 
   epi_data <- validate_epi_data(epi_data)
 
+  # TODO: handle this when creating param grid?
   if (typeof(seasonal_method) == "list") {
     seasonal_method <- seasonal_method[[1]]
   }

R/looping.R

@@ -84,24 +84,27 @@ slide_forecaster <- function(epi_archive,
   )
 }
 
-epix_slide_simple <- function(epi_archive, forecaster, ref_time_values, before, cache_key = NULL) {
+
+epix_slide_simple <- function(epi_archive, forecaster, ref_time_values, before = Inf, cache_key = NULL) {
   # this is so that changing the object without changing the name doesn't result in pulling the wrong cache
   cache_hash <- rlang::hash(epi_archive)
   dir.create(".exploration_cache/slide_cache", showWarnings = FALSE, recursive = TRUE)
-  purrr::map(ref_time_values, function(tv) {
+  out <- purrr::map(ref_time_values, function(tv) {
     if (is.null(cache_key)) {
       epi_df <- epi_archive %>%
-        epix_as_of(tv, min_time_value = tv - before)
+        epix_as_of(min(tv, .$versions_end), min_time_value = tv - before)
     } else {
       file_path <- glue::glue(".exploration_cache/slide_cache/{cache_key}_{cache_hash}_{before}_{tv}.parquet")
       if (file.exists(file_path)) {
         epi_df <- qs::qread(file_path)
       } else {
         epi_df <- epi_archive %>%
-          epix_as_of(tv, min_time_value = tv - before)
+          epix_as_of(min(tv, .$versions_end), min_time_value = tv - before)
         qs::qsave(epi_df, file_path)
       }
     }
     epi_df %>% forecaster()
   }) %>% bind_rows()
+  gc()
+  return(out)
 }

R/scoring.R

@@ -22,9 +22,24 @@ evaluate_predictions <- function(forecasts, truth_data) {
 
   scores <- forecast_obj %>%
     scoringutils::score(metrics = get_metrics(.)) %>%
-    as_tibble() %>%
+    as_tibble()
+  missing_metrics <- setdiff(
+    c("model", "geo_value", "forecast_date", "target_end_date", "wis", "ae_median", "interval_coverage_50", "interval_coverage_90"),
+    names(scores)
+  )
+  if (length(missing_metrics) > 0) {
+    cli::cli_abort(c(
+      "scoring error",
+      "i" = "missing metrics: {missing_metrics}",
+      "i" = "if wis is missing, then likely quantile monotonicity was violated"
+    ))
+  }
+  scores %>%
     select(
-      model, geo_value, forecast_date, target_end_date,
+      model,
+      geo_value,
+      forecast_date,
+      target_end_date,
       wis,
       ae = ae_median,
       coverage_50 = interval_coverage_50,