dev/main comparison

.Rbuildignore

@@ -1,2 +1,4 @@
+^renv$
+^renv\.lock$
 ^.*\.Rproj$
 ^\.Rproj\.user$

DESCRIPTION

@@ -8,15 +8,13 @@ Description: Provide real-time revision forecasts.
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Suggests: 
     testthat (>= 3.0.0)
 Config/testthat/edition: 3
-Imports: 
+Imports:
     arrow,
-    covidcast,
     dplyr,
-    evalcast,
     english,
     jsonlite,
     lubridate,

NAMESPACE

@@ -21,6 +21,7 @@ export(WEEK_ISSUES)
 export(Y7DAV)
 export(YITL)
 export(add_7davs)
+export(aux_feature_names)
 export(add_dayofweek)
 export(add_lagged_terms)
 export(add_log_transformed)
@@ -33,6 +34,7 @@ export(add_weights_related)
 export(assert)
 export(create_dir_not_exist)
 export(create_params_list)
+export(create_target_lookup)
 export(cv_revision_forecast)
 export(data_filteration)
 export(data_preprocessing)
@@ -68,7 +70,6 @@ importFrom(dplyr,slice_max)
 importFrom(dplyr,starts_with)
 importFrom(dplyr,ungroup)
 importFrom(english,english)
-importFrom(evalcast,weighted_interval_score)
 importFrom(jsonlite,read_json)
 importFrom(lubridate,days_in_month)
 importFrom(lubridate,make_date)

R/forecast.R

@@ -41,6 +41,7 @@ revision_forecast <- function(train_data, test_data, taus,
                               smoothed_target=TRUE,
                               lagged_term_list=NULL,
                               params_list=NULL,
+                              extra_params=NULL,
                               temporal_resol="daily",
                               lambda = 0.1, gamma = 0.1,
                               lp_solver=LP_SOLVER, test_lag_group="",
@@ -49,9 +50,11 @@ revision_forecast <- function(train_data, test_data, taus,
                               indicator="testdata", signal="",
                               geo_level="state", signal_suffix="",
                               training_end_date="",
-                              training_days =365,
+                              training_days=365,
                               train_models = TRUE,
-                              make_predictions=TRUE) {
+                              make_predictions=TRUE,
+                              onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun")),
+                              time_limit = NULL) {
 
 
 
@@ -77,7 +80,7 @@ revision_forecast <- function(train_data, test_data, taus,
   }
 
   if (is.null(params_list)) {
-    params_list <- create_params_list(train_data, lagged_term_list, temporal_resol)
+    params_list <- create_params_list(train_data, lagged_term_list, temporal_resol, onehot_weekdays, extra_params)
   }
 
   if (smoothed_target) {
@@ -94,17 +97,42 @@ revision_forecast <- function(train_data, test_data, taus,
 
   test_data_list <- list()
 
-  if (train_models) {
-    sqrt_max_raw <- sqrt(max(train_data$value_7dav, na.rm=TRUE))
-    train_result <- add_sqrtscale(train_data, sqrt_max_raw)
-    train_data <- train_result$data
-    kept_bins <- train_result$kept_bins
-    #for (col in kept_bins) {
-    #  proportion <- sum(train_data[[col]]) / nrow(train_data)
-    #  cat(sprintf("Sum of %s: %.4f\n", col, proportion))
-    #}
-    train_data <- train_data[, c(basic_cols, params_list, extra_cols, kept_bins, response)] %>% drop_na()
+  sqrt_max_raw <- sqrt(max(train_data$value_7dav, na.rm=TRUE))
+  train_result <- add_sqrtscale(train_data, sqrt_max_raw)
+  train_data <- train_result$data
+  kept_bins <- train_result$kept_bins
+  train_data <- train_data[, c(basic_cols, params_list, extra_cols, kept_bins, response)] %>% drop_na()
+
+  .response_sd <- stats::sd(train_data[[response]], na.rm = TRUE)
+  if (is.na(.response_sd)) {
+    # drop_na() removed every row — no usable training data at all.
+    warning(sprintf(
+      "No training rows after preprocessing [geo=%s lag_group=%s]; skipping",
+      geo, test_lag_group
+    ))
+    return(data.frame())
   }
+  if (.response_sd < 1e-8) {
+    # Constant response (sd ~ 0): fill_missing_updates synthesised the entire lag
+    # group from forward-filled zeros. GLPK cycles indefinitely on degenerate LPs,
+    # so skip the solver and return the constant as the prediction directly.
+    warning(sprintf(
+      "Constant training response [geo=%s lag_group=%s]; predicting constant — likely all synthetic data",
+      geo, test_lag_group
+    ))
+    if (!make_predictions) return(data.frame())
+    .constant_val <- mean(train_data[[response]], na.rm = TRUE)
+    test_out <- test_data[, intersect(c(basic_cols, response), colnames(test_data)), drop = FALSE]
+    test_out <- tidyr::drop_na(test_out, dplyr::all_of(basic_cols))
+    test_out[paste0("predicted_tau", taus)] <- .constant_val
+    if (response %in% colnames(test_out)) {
+      test_out <- evaluate(test_out, taus, response = response)
+    }
+    test_out$gamma <- gamma[1]
+    test_out$lambda <- lambda[1]
+    return(as.data.frame(test_out))
+  }
+  rm(.response_sd)
 
   # pre-process the test data with max_raw
   if (make_predictions) {
@@ -120,7 +148,7 @@ revision_forecast <- function(train_data, test_data, taus,
     # Get the trained_model
     obj <- get_model(model_path, train_data, params_list, response, taus,
                      sqrt_max_raw, kept_bins,
-                     lambda[1], gamma[1], lp_solver, train_models)
+                     lambda[1], gamma[1], lp_solver, train_models, time_limit)
 
     sqrt_max_raw <- attr(obj, "sqrt_max_raw")
     kept_bins <- attr(obj, "kept_bins")
@@ -150,7 +178,7 @@ revision_forecast <- function(train_data, test_data, taus,
 
       # Get the trained_model
       obj <- get_model(model_path, train_data, params_list, response, taus, sqrt_max_raw,
-                       l, g, lp_solver, train_models)
+                       l, g, lp_solver, train_models, time_limit)
 
       if (make_predictions) {
         test_data <- get_prediction(test_data, taus, params_list, response, obj,
@@ -317,6 +345,7 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
                      smoothed_target=TRUE,
                      lagged_term_list=NULL,
                      params_list=NULL,
+                     extra_params=NULL,
                      lambda=LAMBDA, gamma=GAMMA, lag_pad=LAG_PAD,
                      temporal_resol="daily",
                      lp_solver=LP_SOLVER,
@@ -327,7 +356,9 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
                      training_end_date="",
                      training_days=365,
                      train_models = TRUE,
-                     make_predictions = TRUE) {
+                     make_predictions = TRUE,
+                     onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun")),
+                     time_limit = NULL) {
 
   testing_start_date <- as.Date(testing_start_date)
 
@@ -340,8 +371,8 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
 
       # Detect weekly spacing
       if (length(lag_diffs) == 1 && lag_diffs == 7) {
+        if (temporal_resol != "weekly") message("Auto-detected weekly temporal resolution from lag spacing.")
         temporal_resol <- "weekly"
-        message("Auto-detected weekly temporal resolution from lag spacing.")
       }
     }
   }
@@ -388,13 +419,14 @@ DelphiRF <- function(df, testing_start_date, taus=TAUS,
 
     results <- revision_forecast(train_data, test_data, taus,
                                  smoothed_target, lagged_term_list,
-                                 params_list, temporal_resol,
+                                 params_list, extra_params, temporal_resol,
                                  l, g, lp_solver, test_lag_group,
                                  geo, value_type, model_save_dir,
                                  indicator, signal, geo_level,
                                  signal_suffix, as.character(testing_start_date),
                                  training_days, train_models,
-                                 make_predictions)
+                                 make_predictions, onehot_weekdays,
+                                 time_limit = time_limit)
 
     test_data_list <- append(test_data_list, list(results))
   }

R/model.R

@@ -142,17 +142,26 @@ get_prediction <- function(test_data, taus, covariates, response, obj,
   return (as.data.frame(test_data))
 }
 
+#' Weighted interval score for a single observation
+#'
+#' Inlined from the evalcast package.
+#'
+#' @param taus Numeric vector of quantile levels.
+#' @param residuals Numeric vector of (quantile_prediction - actual) values.
+#' @param point_pred Unused; kept for interface compatibility.
+#' @keywords internal
+weighted_interval_score <- function(taus, residuals, point_pred) {
+  alpha <- 2 * pmin(taus, 1 - taus)
+  mean(alpha * (abs(residuals) + (residuals) * (2 * (taus >= 0.5) - 1)))
+}
+
 #' Evaluation of the test results based on WIS score
-#' The WIS score calculation is based on the weighted_interval_score function
-#' from the `evalcast` package from Delphi
 #'
 #' @param test_data dataframe with a column containing the prediction results of
 #'    each requested quantile. Each row represents an update with certain
 #'    (reference_date, report_date, location) combination.
 #' @template taus-template
 #'
-#' @importFrom evalcast weighted_interval_score
-#'
 #' @export
 evaluate <- function(test_data, taus, response) {
   n_row <- nrow(test_data)
@@ -209,11 +218,10 @@ exponentiate_preds <- function(test_data, taus) {
 #' @importFrom quantgen quantile_lasso
 get_model <- function(model_path, train_data, covariates, response, tau,
                       sqrt_max_raw, kept_bins,
-                      lambda, gamma, lp_solver, train_models) {
+                      lambda, gamma, lp_solver, train_models, time_limit = NULL) {
   if (train_models || !file.exists(model_path)) {
     if (!train_models && !file.exists(model_path)) {
-      warning(str_interp("user requested use of cached model but file {model_path}"),
-        " does not exist; training new model")
+      warning(str_interp("user requested use of cached model but file ${model_path} does not exist; training new model"))
     }
     # Quantile regression
     vec_7dav <- train_data[["value_7dav_diff"]]
@@ -225,11 +233,15 @@ get_model <- function(model_path, train_data, covariates, response, tau,
     } else {
       weights <- NULL
     }
-    obj <- quantile_lasso(as.matrix(train_data[covariates]),
-                         train_data[[response]], # - train_data[["log_value_7dav"]],
-                         tau = tau,
-                         lambda = lambda, standardize = TRUE, lp_solver = lp_solver, intercept=TRUE,
-                         weights = weights)
+    lasso_args <- list(
+      as.matrix(train_data[covariates]),
+      train_data[[response]],
+      tau = tau,
+      lambda = lambda, standardize = TRUE, lp_solver = lp_solver, intercept = TRUE,
+      weights = weights
+    )
+    if (!is.null(time_limit)) lasso_args$time_limit <- time_limit
+    obj <- do.call(quantile_lasso, lasso_args)
 
     # Save model to cache.
     create_dir_not_exist(dirname(model_path))
@@ -243,7 +255,7 @@ get_model <- function(model_path, train_data, covariates, response, tau,
   } else {
     # Load model from cache invisibly. Object has the same name as the original
     # model object, `obj`.
-    print(str_interp("Loading from ${model_path}"))
+    message(str_interp("Loading from ${model_path}"))
     obj <- readRDS(model_path)
   }
 
@@ -332,27 +344,26 @@ generate_filename <- function(indicator, signal,
 #'
 #' @importFrom dplyr mutate select
 #'
-create_params_list <- function(train_data, lagged_term_list, temporal_resol) {
+create_params_list <- function(train_data, lagged_term_list, temporal_resol,
+                               onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun")),
+                               extra_params = NULL) {
   params_list <- c(
     WEEK_ISSUES[1],
     Y7DAV,
     paste0("log_value_7dav_lag", lagged_term_list),
     paste0("log_delta_value_7dav_lag", lagged_term_list)
   )
-  # Include log lag adjustments if multiple lags exist
-  if (length(unique(train_data$lag)) > 1){
+  if (length(unique(train_data$lag)) > 1) {
     params_list <- c(params_list, LOG_LAG)
   }
 
-  dayofweek <- c("Mon", "Weekends")
-  extra_params_for_daily <- c(
-    paste0(dayofweek, "_ref"),
-    paste0(dayofweek, "_issue")
-  )
-
-  if (temporal_resol == "daily"){
-    return (c(params_list, extra_params_for_daily))
+  group_names <- if (!is.null(names(onehot_weekdays))) {
+    names(onehot_weekdays)
   } else {
-    return(params_list)
+    vapply(onehot_weekdays, function(grp) paste0(grp, collapse = ""), character(1))
   }
+  extra_params_for_daily <- c(paste0(group_names, "_ref"), paste0(group_names, "_issue"))
+
+  base_params <- if (temporal_resol == "daily") c(params_list, extra_params_for_daily) else params_list
+  if (!is.null(extra_params)) c(base_params, extra_params) else base_params
 }