diff --git a/.Rbuildignore b/.Rbuildignore
index 91114bf..d821302 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,2 +1,4 @@
+^renv$
+^renv\.lock$
 ^.*\.Rproj$
 ^\.Rproj\.user$
diff --git a/DESCRIPTION b/DESCRIPTION
index db43f7e..f929e62 100644
--- a/DESCRIPTION
+++ b/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,
diff --git a/NAMESPACE b/NAMESPACE
index 0c0d5a1..33a664f 100644
--- a/NAMESPACE
+++ b/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)
diff --git a/R/feature_engineering.R b/R/feature_engineering.R
index c253065..afde37d 100644
--- a/R/feature_engineering.R
+++ b/R/feature_engineering.R
@@ -32,6 +32,37 @@ add_dayofweek <- function(df, time_col, suffix, wd = WEEKDAYS_ABBR) {
 }
 
 
+#' Add grouped day-of-week one-hot columns
+#'
+#' Creates one binary column per group in `onehot_weekdays`. Each group is a
+#' character vector of day abbreviations (from `WEEKDAYS_ABBR`); a row is 1 if
+#' the date falls on any day in the group.  Column names are derived from the
+#' list names when present, otherwise by concatenating the day abbreviations.
+#'
+#' @param df A data frame containing the date column.
+#' @param time_col Name of the date column.
+#' @param suffix Column name suffix (e.g. `"_ref"` or `"_issue"`).
+#' @param onehot_weekdays Named or unnamed list of character vectors, each
+#'   specifying a group of days (e.g. `list(Mon=c("Mon"), Weekends=c("Sat","Sun"))`).
+#'
+#' @return `df` with one additional integer column per group.
+#' @export
+add_grouped_dayofweek <- function(df, time_col, suffix, onehot_weekdays) {
+  df <- df %>% mutate({{ time_col }} := as.Date(.data[[time_col]]))
+  dayofweek <- as.numeric(format(df[[time_col]], format = "%u"))
+  group_names <- if (!is.null(names(onehot_weekdays))) {
+    names(onehot_weekdays)
+  } else {
+    vapply(onehot_weekdays, function(grp) paste0(grp, collapse = ""), character(1))
+  }
+  for (ii in seq_along(onehot_weekdays)) {
+    day_indices <- match(onehot_weekdays[[ii]], WEEKDAYS_ABBR)
+    df[[paste0(group_names[ii], suffix)]] <- as.integer(dayofweek %in% day_indices)
+  }
+  df
+}
+
+
 #' Add one-hot encoding for week of the month based on issue date
 #'
 #' This function calculates the week of the month for each date in the specified
@@ -188,8 +219,13 @@ add_lagged_terms <- function(df, value_col, refd_col, lag_col, lagged_term_list=
 #' @importFrom dplyr rename
 #' @export
 add_targets <- function(df, value_col, refd_col, lag_col, ref_lag, temporal_resol) {
-  # Add target
-  target_df <- df[df[[lag_col]]==ref_lag, c(refd_col, "report_date", value_col, "value_7dav")]
+  available_lags <- sort(unique(df[[lag_col]]))
+  effective_ref_lag <- min(available_lags[available_lags >= ref_lag])
+  if (is.infinite(effective_ref_lag))
+    stop(sprintf("ref_lag %d exceeds all available lags (max %d)", ref_lag, max(available_lags)))
+  if (effective_ref_lag != ref_lag)
+    message(sprintf("ref_lag %d not available; using next lag %d", ref_lag, effective_ref_lag))
+  target_df <- df[df[[lag_col]] == effective_ref_lag, c(refd_col, "report_date", value_col, "value_7dav")]
   # Rename columns for clarity
   target_df <- target_df %>%
     dplyr::rename(
@@ -202,6 +238,90 @@ add_targets <- function(df, value_col, refd_col, lag_col, ref_lag, temporal_reso
   return (as.data.frame(backfill_df))
 }
 
+#' Construct one raw-aware target per reference date
+#'
+#' Selects the latest genuine value-changing revision in
+#' `[ref_lag - lower_tolerance, ref_lag + upper_tolerance]`. If none exists,
+#' it falls back to the latest raw value at or before the lower boundary.
+#' The returned table is sparse (one row per reference date) and is intended
+#' to be joined after feature-grid filling.
+#'
+#' @param df Raw reporting-triangle data frame.
+#' @param value_col Name of the value column used to identify revisions.
+#' @param refd_col Name of the reference-date column.
+#' @param lag_col Name of the day-based reporting-lag column.
+#' @param ref_lag Central target lag, in days.
+#' @param lower_tolerance Non-negative days before `ref_lag` included in the
+#'   genuine-revision search window.
+#' @param upper_tolerance Non-negative days after `ref_lag` included in the
+#'   genuine-revision search window.
+#' @param temporal_resol Either `"daily"` or `"weekly"`.
+#' @return A compact data frame with at most one target per reference date.
+#' @export
+create_target_lookup <- function(df, value_col, refd_col, lag_col, ref_lag,
+                                 lower_tolerance = 0, upper_tolerance = 0,
+                                 temporal_resol = "daily") {
+  if (lower_tolerance < 0 || upper_tolerance < 0) {
+    stop("Target lag tolerances must be non-negative.")
+  }
+  if (nrow(df) == 0) {
+    return(data.frame(reference_date = as.Date(character()),
+                      target_date = as.Date(character()),
+                      target_lag = numeric(), target_type = character()))
+  }
+  raw <- df[, unique(c(refd_col, lag_col, value_col)), drop = FALSE]
+  raw[[refd_col]] <- as.Date(raw[[refd_col]])
+  raw$report_date <- raw[[refd_col]] + raw[[lag_col]]
+  if (temporal_resol == "weekly") {
+    raw <- normalize_weekly_observations(raw, refd_col, lag_col)
+  } else if (temporal_resol != "daily") {
+    stop("Invalid temporal_resol. Choose either 'daily' or 'weekly'.")
+  }
+  raw <- raw[raw[[lag_col]] >= 0, , drop = FALSE]
+  raw <- raw[order(raw[[refd_col]], raw[[lag_col]], raw$report_date), , drop = FALSE]
+  lower <- ref_lag - lower_tolerance
+  upper <- ref_lag + upper_tolerance
+
+  chosen <- lapply(split(raw, raw[[refd_col]]), function(g) {
+    values <- g[[value_col]]
+    changed <- c(TRUE, (is.na(values[-1]) != is.na(values[-length(values)])) |
+      (!is.na(values[-1]) & !is.na(values[-length(values)]) &
+         values[-1] != values[-length(values)]))
+    candidates <- which(changed & g[[lag_col]] >= lower & g[[lag_col]] <= upper)
+    if (length(candidates) > 0) {
+      idx <- tail(candidates, 1)
+      type <- "revision"
+    } else {
+      fallback <- which(g[[lag_col]] <= lower)
+      if (length(fallback) == 0) return(NULL)
+      idx <- tail(fallback, 1)
+      type <- "fallback"
+    }
+    data.frame(
+      reference_date = as.Date(g[[refd_col]][idx]),
+      target_date = as.Date(g$report_date[idx]),
+      target_lag = as.numeric(g[[lag_col]][idx]),
+      target_type = type
+    )
+  })
+  result <- dplyr::bind_rows(chosen)
+  if (nrow(result) == 0) {
+    return(data.frame(reference_date = as.Date(character()),
+                      target_date = as.Date(character()),
+                      target_lag = numeric(), target_type = character()))
+  }
+  result
+}
+
+attach_target_lookup <- function(df, target_lookup) {
+  target_values <- df %>%
+    select(reference_date, target_date = report_date,
+           value_target = value_raw, value_target_7dav = value_7dav)
+  df %>%
+    left_join(target_lookup, by = "reference_date") %>%
+    left_join(target_values, by = c("reference_date", "target_date"))
+}
+
 
 #' Add log-transformed columns for specified numerical variables
 #'
@@ -243,31 +363,98 @@ add_log_transformed <- function(df, lagged_term_list) {
 #' @param lag_col Column name representing the lag between the reference and issue date.
 #' @param temporal_resol A string indicating the temporal resolution ("daily" or "weekly").
 #'                       Defaults to "daily".
+#' @param onehot_weekdays Named or unnamed list of character vectors defining day groups
+#'   for one-hot encoding.  Each group becomes one binary column per date column.
+#'   Names, when present, are used as column prefixes; otherwise day abbreviations are
+#'   concatenated.  Defaults to `list(Mon=c("Mon"), Weekends=c("Sat","Sun"))`.
 #'
 #' @details
-#' - If `temporal_resol` is "daily", one-hot encoded day-of-week columns are added
+#' - If `temporal_resol` is "daily", one-hot encoded day-of-week group columns are added
 #'   for both `refd_col` (reference date) and `"report_date"`.
 #' - One-hot encoded week-of-month columns are added for `"report_date"` in all cases.
 #'
 #' @return A modified data frame with additional date-related feature columns.
 #'
 #' @export
-add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol="daily") {
+add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol = "daily",
+                                  onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun"))) {
   df$report_date <- df[[refd_col]] + df[[lag_col]]
-  if (temporal_resol=="daily"){
-    # Add columns for day-of-week effect
-    df <- add_dayofweek(df, refd_col, "_ref", WEEKDAYS_ABBR)
-    df <- add_dayofweek(df, "report_date", "_issue", WEEKDAYS_ABBR)
-    # Add columns for weekends
-    df$Weekends_issue <- as.integer(df$Sat_issue == 1 | df$Sun_issue == 1)
-    df$Weekends_ref <- as.integer(df$Sat_ref == 1 | df$Sun_ref == 1)
+  if (temporal_resol == "daily") {
+    df <- add_grouped_dayofweek(df, refd_col, "_ref", onehot_weekdays)
+    df <- add_grouped_dayofweek(df, "report_date", "_issue", onehot_weekdays)
   }
-  # Add columns for week-of-month effect
   df <- add_weekofmonth(df, "report_date", WEEK_ISSUES)
+  return(as.data.frame(df))
+}
 
-  return (as.data.frame(df))
+#' Process an auxiliary reporting triangle into prefixed feature columns
+#'
+#' Runs a single auxiliary data frame (one geo, one signal) through the
+#' fill → 7-day-average → lagged-terms → log-transform pipeline and
+#' renames all value-derived columns with a `{name}_` prefix so they can
+#' be safely joined to the primary preprocessed data frame.
+#'
+#' @param df Data frame with columns `reference_date`, `report_date`, `lag`,
+#'   and `value` (the signal values).
+#' @param name Character scalar used as the column prefix (e.g. `"nssp"`).
+#' @param lagged_term_list Numeric vector of lag values (same as used for the
+#'   primary signal).
+#' @param temporal_resol `"daily"` or `"weekly"`.
+#' @param smoothed Logical; if `FALSE` and `temporal_resol == "daily"`, a 7-day
+#'   moving average is computed. Otherwise `value_7dav` is set equal to
+#'   `value_raw`.
+#'
+#' @return Data frame with columns `reference_date`, `report_date`, `lag`, and
+#'   all value/log columns prefixed with `{name}_`.  Useful predictors are
+#'   `{name}_log_value_7dav_lag{N}` and `{name}_log_delta_value_7dav_lag{N}`;
+#'   see [aux_feature_names()].
+process_aux_triangle <- function(df, name, lagged_term_list, temporal_resol, smoothed,
+                                max_report_override = NULL) {
+  filled_df <- fill_missing_updates(df, "value", "reference_date", "lag", temporal_resol,
+                                    max_report_override = max_report_override)
+  if (nrow(filled_df) == 0) {
+    expected_cols <- c("reference_date", "report_date", "lag",
+                       aux_feature_names(name, lagged_term_list))
+    return(setNames(data.frame(matrix(ncol = length(expected_cols), nrow = 0)),
+                    expected_cols))
+  }
+  if (!smoothed && temporal_resol == "daily") {
+    filled_df <- add_7davs(filled_df, "value_raw", "reference_date", "lag")
+  } else {
+    filled_df$value_7dav <- filled_df$value_raw
+  }
+  filled_df <- add_lagged_terms(
+    filled_df, "value_7dav", "reference_date", "lag", lagged_term_list, temporal_resol
+  )
+  filled_df <- add_log_transformed(filled_df, lagged_term_list)
+
+  value_cols <- grep("^(value_|log_)", colnames(filled_df), value = TRUE)
+  colnames(filled_df)[colnames(filled_df) %in% value_cols] <- paste0(name, "_", value_cols)
+
+  filled_df[, c("reference_date", "report_date", "lag", paste0(name, "_", value_cols))]
 }
 
+
+#' Return the feature column names produced by an auxiliary triangle
+#'
+#' Gives the column names that [process_aux_triangle()] adds for a given
+#' auxiliary signal, matching the log-value and log-delta features used by
+#' the primary signal in [create_params_list()].
+#'
+#' @param name Character scalar; the aux triangle name (must match what was
+#'   passed to [data_preprocessing()]).
+#' @param lagged_term_list Numeric vector of lag values.
+#'
+#' @return Character vector of feature column names.
+#' @export
+aux_feature_names <- function(name, lagged_term_list) {
+  c(
+    paste0(name, "_log_value_7dav_lag", lagged_term_list),
+    paste0(name, "_log_delta_value_7dav_lag", lagged_term_list)
+  )
+}
+
+
 #' Data Preprocessing Function
 #'
 #' This function processes input data by handling missing values, computing lagged terms,
@@ -287,14 +474,29 @@ add_params_for_dates <- function(df, refd_col, lag_col, temporal_resol="daily")
 #' @param value_type Character indicating the type of values ('count' or 'fraction').
 #' @param temporal_resol Character specifying temporal resolution ('daily' or 'weekly').
 #' @param smoothed Logical indicating whether smoothing should be applied.
-#'
-#' @importFrom dplyr full_join distinct
+#' @param target_lag_lower_tolerance Non-negative days before `ref_lag` to
+#'   search for genuine target revisions.
+#' @param target_lag_upper_tolerance Non-negative days after `ref_lag` to
+#'   search for genuine target revisions.
+#' @param aux_triangles Named list of auxiliary reporting-triangle data frames.
+#'   Each element must have columns `reference_date`, `report_date`, `lag`, and
+#'   `value` (already filtered to the same single geo as `df`).  Each is run
+#'   through the same fill → lag → log pipeline as the primary signal and
+#'   joined to the result; columns are prefixed with the list element name.
+#'   Use [aux_feature_names()] to obtain the resulting predictor column names
+#'   for [create_params_list()].
+#'
+#' @importFrom dplyr full_join left_join distinct
 #' @importFrom english english
 #'
 #' @export
 data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag,
                                suffixes=c(""), lagged_term_list = NULL, value_type="count",
-                               temporal_resol="daily", smoothed=FALSE) {
+                               temporal_resol="daily", smoothed=FALSE,
+                               target_lag_lower_tolerance = 0,
+                               target_lag_upper_tolerance = 0,
+                               aux_triangles = NULL,
+                               onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun"))) {
   if (value_type == "count") {
     if (length(value_col) > 1) warning("Multiple value column names provided; only the first one will be used.")
     if (length(unique(suffixes)) > 1) warning("Multiple suffixes provided; only the first one will be used.")
@@ -332,6 +534,11 @@ data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag,
     lagged_term_list <- c(lagged_term_list, 7)
   } # Make sure we always have the lagged term from last week
 
+  target_lookup <- create_target_lookup(
+    df, value_col[1], refd_col, lag_col, ref_lag,
+    target_lag_lower_tolerance, target_lag_upper_tolerance, temporal_resol
+  )
+
   dfList <- lapply(value_col, function(value_col) {
     filled_df <- fill_missing_updates(df, value_col, refd_col, lag_col, temporal_resol)
     if (!(smoothed) & (temporal_resol == "daily")){
@@ -340,14 +547,15 @@ data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag,
       filled_df$value_7dav = filled_df$value_raw
     }
     filled_df <- add_lagged_terms(filled_df, "value_7dav", "reference_date", "lag", lagged_term_list, temporal_resol)
-    filled_df <- add_targets(filled_df, "value_raw", "reference_date", "lag", ref_lag, temporal_resol)
+    filled_df <- attach_target_lookup(filled_df, target_lookup)
     add_log_transformed(filled_df, lagged_term_list)
   })
 
   merged_df <- Reduce(
     function(x, y) full_join(
       x, y,
-      by = c("reference_date", "report_date", "lag", "target_date"),
+      by = c("reference_date", "report_date", "lag", "target_date",
+             "target_lag", "target_type"),
       suffix = suffixes
     ),
     dfList
@@ -367,7 +575,26 @@ data_preprocessing <- function(df, value_col, refd_col, lag_col, ref_lag,
   }
 
   merged_df$inv_log_lag <- 1/(merged_df$lag + 1)
-  merged_df <- add_params_for_dates(merged_df, "reference_date", "lag", temporal_resol)
+  merged_df <- add_params_for_dates(merged_df, "reference_date", "lag", temporal_resol, onehot_weekdays)
+
+  if (!is.null(aux_triangles)) {
+    primary_max_report <- max(merged_df$report_date)
+    for (nm in names(aux_triangles)) {
+      aux_processed <- process_aux_triangle(
+        aux_triangles[[nm]], nm, lagged_term_list, temporal_resol, smoothed,
+        max_report_override = primary_max_report
+      )
+      if (nrow(aux_processed) == 0L) {
+        merged_df <- merged_df[0L, ]
+        break
+      }
+      merged_df <- dplyr::left_join(
+        merged_df, aux_processed,
+        by = c("reference_date", "report_date", "lag")
+      )
+      merged_df <- merged_df[!is.na(merged_df[[paste0(nm, "_value_raw")]]), ]
+    }
+  }
 
   merged_df <- merged_df %>%
     filter(.data$lag < ref_lag)
diff --git a/R/forecast.R b/R/forecast.R
index 5ab1b21..cbbecaa 100644
--- a/R/forecast.R
+++ b/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,10 @@ 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"))) {
 
 
 
@@ -77,7 +79,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 +96,11 @@ 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()
 
   # pre-process the test data with max_raw
   if (make_predictions) {
@@ -317,6 +313,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 +324,8 @@ 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"))) {
 
   testing_start_date <- as.Date(testing_start_date)
 
@@ -340,8 +338,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 +386,13 @@ 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)
 
     test_data_list <- append(test_data_list, list(results))
   }
diff --git a/R/model.R b/R/model.R
index 1af536f..d814b0b 100644
--- a/R/model.R
+++ b/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)
@@ -210,11 +219,17 @@ exponentiate_preds <- function(test_data, taus) {
 get_model <- function(model_path, train_data, covariates, response, tau,
                       sqrt_max_raw, kept_bins,
                       lambda, gamma, lp_solver, train_models) {
-  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")
+  if (!train_models && !file.exists(model_path)) {
+    candidates <- list.files(dirname(model_path), pattern = "\\.rds$", full.names = TRUE)
+    if (length(candidates) > 0L) {
+      alt_path <- sort(candidates, decreasing = TRUE)[1L]
+      message(sprintf("Exact model not found; loading most recent cached model: %s", basename(alt_path)))
+      return(readRDS(alt_path))
     }
+    warning(str_interp("user requested use of cached model but file ${model_path} does not exist and no alternatives found; training new model"))
+    train_models <- TRUE
+  }
+  if (train_models || !file.exists(model_path)) {
     # Quantile regression
     vec_7dav <- train_data[["value_7dav_diff"]]
     vec_slope <- train_data[["value_slope_diff"]]
@@ -332,27 +347,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
 }
diff --git a/R/preprocessing.R b/R/preprocessing.R
index 4371cec..85d1d3a 100644
--- a/R/preprocessing.R
+++ b/R/preprocessing.R
@@ -36,6 +36,21 @@ fill_rows <- function(df, refd_col, lag_col, min_refd, max_refd, ref_lag) {
   return (df_new)
 }
 
+normalize_weekly_observations <- function(df, refd_col, lag_col) {
+  if (nrow(df) == 0) return(df)
+  original_reference_date <- as.Date(df[[refd_col]])
+  original_report_date <- original_reference_date + df[[lag_col]]
+  epiweek_end <- function(x) x + ((6L - as.POSIXlt(x)$wday) %% 7L)
+  df[[refd_col]] <- epiweek_end(original_reference_date)
+  df$report_date <- epiweek_end(original_report_date)
+  df[[lag_col]] <- as.numeric(df$report_date - df[[refd_col]])
+  ordering <- order(df[[refd_col]], df$report_date,
+                    original_report_date, original_reference_date)
+  df <- df[ordering, , drop = FALSE]
+  key <- paste(df[[refd_col]], df$report_date, sep = "\r")
+  df[!duplicated(key, fromLast = TRUE), , drop = FALSE]
+}
+
 #' Fill missing updates in a time series dataset with lagged values
 #' Get pivot table, filling NANs. If there is no update on issue date D but
 #' previous reports exist for issue date D_p < D, all the dates between
@@ -56,7 +71,8 @@ fill_rows <- function(df, refd_col, lag_col, min_refd, max_refd, ref_lag) {
 #' @importFrom tidyr fill pivot_wider pivot_longer replace_na expand_grid
 #' @importFrom dplyr %>% select left_join mutate arrange distinct filter everything
 #' @export
-fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_resol="daily") {
+fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_resol = "daily",
+                              max_report_override = NULL) {
   df <- df %>% distinct()  # Remove duplicates if any
 
   if (nrow(df) == 0) {
@@ -64,6 +80,10 @@ fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_reso
   }
   df$report_date <- df[[refd_col]] + df[[lag_col]]
 
+  if (temporal_resol == "weekly") {
+    df <- normalize_weekly_observations(df, refd_col, lag_col)
+  }
+
   # Generate a sequence of all possible dates
   if (temporal_resol == "daily"){
     all_reference_dates <- seq(min(df[[refd_col]]), max(df[[refd_col]]), by = "day")
@@ -73,21 +93,19 @@ fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_reso
     all_reference_dates <- seq(min(df[[refd_col]]), max(df[[refd_col]]), by = "7 days")
     all_report_dates <- seq(min(df[["report_date"]]), max(df[["report_date"]]), by = "7 days")
 
-    # Check if all reference dates in df are within the generated sequence
-    if (!all(df[[refd_col]] %in% all_reference_dates)) {
-      stop("The reference dates do not regularly have a gap of 7 days. Some reference dates will be ignored. Please check your input data.")
-    }
-
-    # Check if all reference dates in df are within the generated sequence
-    if (!all(df$report_date %in% all_report_dates)) {
-      stop("The report dates do not regularly have a gap of 7 days. Some report dates will be ignored. Please check your input data.")
-    }
-
     gap <- 7
   } else {
     stop("Invalid temporal_resol. Choose either 'daily' or 'weekly'.")
   }
 
+  if (!is.null(max_report_override)) {
+    max_report_override <- as.Date(max_report_override)
+    if (max_report_override > max(all_report_dates)) {
+      extra <- seq(max(all_report_dates) + gap, max_report_override, by = gap)
+      all_report_dates <- c(all_report_dates, extra)
+    }
+  }
+
   # Create a complete grid of all combinations of reference_date and report_date
   complete_grid <- tidyr::expand_grid(
     !!refd_col := all_reference_dates,
@@ -118,10 +136,10 @@ fill_missing_updates <- function(df, value_col, refd_col, lag_col, temporal_reso
   backfill_df <- pivot_df %>%
     tidyr::pivot_longer(-report_date, values_to = "value_raw", names_to = refd_col) %>%
     mutate(
-      reference_date := as.Date(.data[[refd_col]]),
+      reference_date = as.Date(.data[[refd_col]]),
       lag = as.numeric(report_date - reference_date)
     ) %>%
-    filter(lag>=0)
+    filter(lag >= 0)
   return (as.data.frame(backfill_df))
 }
 
@@ -177,6 +195,3 @@ get_weekofmonth <- function(date) {
 
   return(week_number)
 }
-
-
-
diff --git a/R/utils.R b/R/utils.R
index 57f88f8..7103110 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -200,19 +200,13 @@ training_days_check <- function(report_date, training_days) {
 
 #' Subset list of counties to those included in the 200 most populous in the US
 #'
+#' Requires the covidcast package, which is not installed by default.
+#'
 #' @importFrom dplyr select %>% arrange desc pull
 #' @importFrom rlang .data
 #' @importFrom utils head
 get_populous_counties <- function() {
-  return(
-    covidcast::county_census %>%
-      dplyr::select(pop = .data$POPESTIMATE2019, fips = .data$FIPS) %>%
-      # Drop megacounties (states)
-      filter(!endsWith(.data$fips, "000")) %>%
-      arrange(desc(.data$pop)) %>%
-      pull(.data$fips) %>%
-      head(n=200)
-  )
+  stop("get_populous_counties() requires the covidcast package. Install it with renv::install(\"cmu-delphi/covidcast/R-packages/covidcast\").")
 }
 
 #' Write a message to the console with the current time
diff --git a/man/create_target_lookup.Rd b/man/create_target_lookup.Rd
new file mode 100644
index 0000000..65ab270
--- /dev/null
+++ b/man/create_target_lookup.Rd
@@ -0,0 +1,44 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/feature_engineering.R
+\name{create_target_lookup}
+\alias{create_target_lookup}
+\title{Construct one raw-aware target per reference date}
+\usage{
+create_target_lookup(
+  df,
+  value_col,
+  refd_col,
+  lag_col,
+  ref_lag,
+  lower_tolerance = 0,
+  upper_tolerance = 0,
+  temporal_resol = "daily"
+)
+}
+\arguments{
+\item{df}{Raw reporting-triangle data frame.}
+
+\item{value_col}{Name of the value column used to identify revisions.}
+
+\item{refd_col}{Name of the reference-date column.}
+
+\item{lag_col}{Name of the day-based reporting-lag column.}
+
+\item{ref_lag}{Central target lag, in days.}
+
+\item{lower_tolerance}{Non-negative days before \code{ref_lag} included in the genuine-revision search window.}
+
+\item{upper_tolerance}{Non-negative days after \code{ref_lag} included in the genuine-revision search window.}
+
+\item{temporal_resol}{Either \code{"daily"} or \code{"weekly"}.}
+}
+\value{
+A compact data frame with at most one target per reference date.
+}
+\description{
+Selects the latest genuine value-changing revision in
+\verb{[ref_lag - lower_tolerance, ref_lag + upper_tolerance]}. If none exists,
+it falls back to the latest raw value at or before the lower boundary.
+The returned table is sparse (one row per reference date) and is intended
+to be joined after feature-grid filling.
+}
diff --git a/man/data_preprocessing.Rd b/man/data_preprocessing.Rd
index 16c9ef9..3faf482 100644
--- a/man/data_preprocessing.Rd
+++ b/man/data_preprocessing.Rd
@@ -14,7 +14,11 @@ data_preprocessing(
   lagged_term_list = NULL,
   value_type = "count",
   temporal_resol = "daily",
-  smoothed = FALSE
+  smoothed = FALSE,
+  target_lag_lower_tolerance = 0,
+  target_lag_upper_tolerance = 0,
+  aux_triangles = NULL,
+  onehot_weekdays = list(Mon = c("Mon"), Weekends = c("Sat", "Sun"))
 )
 }
 \arguments{
@@ -37,6 +41,14 @@ data_preprocessing(
 \item{temporal_resol}{Character specifying temporal resolution ('daily' or 'weekly').}
 
 \item{smoothed}{Logical indicating whether smoothing should be applied.}
+
+\item{target_lag_lower_tolerance}{Non-negative days before \code{ref_lag} included in the genuine-revision target window.}
+
+\item{target_lag_upper_tolerance}{Non-negative days after \code{ref_lag} included in the genuine-revision target window.}
+
+\item{aux_triangles}{Optional named list of auxiliary reporting triangles.}
+
+\item{onehot_weekdays}{Optional named weekday groups used for daily features.}
 }
 \description{
 This function processes input data by handling missing values, computing lagged terms,
diff --git a/tests/testthat/test-feature_engineering.R b/tests/testthat/test-feature_engineering.R
index c4a577f..2a11ad6 100644
--- a/tests/testthat/test-feature_engineering.R
+++ b/tests/testthat/test-feature_engineering.R
@@ -158,6 +158,89 @@ test_that("add_targets correctly adds target columns", {
   expect_equal(unique(df_new[df_new$ref_date == as.Date("2022-01-10"), "value_target"]), NA_real_)
 })
 
+test_that("create_target_lookup selects latest genuine revision in window", {
+  raw <- data.frame(
+    ref_date = as.Date("2024-01-01"),
+    lag = c(50, 59, 61, 62, 70),
+    value = c(10, 10, 12, 15, 20)
+  )
+  target <- create_target_lookup(raw, "value", "ref_date", "lag", 60, 1, 2)
+  expect_equal(target$target_lag, 62)
+  expect_equal(target$target_date, as.Date("2024-03-03"))
+  expect_equal(target$target_type, "revision")
+})
+
+test_that("create_target_lookup falls back to latest raw value at lower bound", {
+  raw <- data.frame(
+    ref_date = as.Date("2024-01-01"),
+    lag = c(40, 55, 59, 70),
+    value = c(8, 10, 10, 20)
+  )
+  target <- create_target_lookup(raw, "value", "ref_date", "lag", 60, 1, 2)
+  expect_equal(target$target_lag, 59)
+  expect_equal(target$target_type, "fallback")
+})
+
+test_that("create_target_lookup never selects negative reporting lags", {
+  raw <- data.frame(
+    ref_date = as.Date("2024-01-10"),
+    lag = c(-7, 5),
+    value = c(5, 8)
+  )
+  target <- create_target_lookup(raw, "value", "ref_date", "lag", 7)
+  expect_equal(target$target_lag, 5)
+  expect_equal(target$target_type, "fallback")
+})
+
+test_that("weekly raw targets retain Friday and keep day-based lags", {
+  raw <- data.frame(
+    ref_date = as.Date("2024-01-06"),
+    lag = c(60, 62), # Wednesday and Friday in the same epiweek
+    value = c(10, 14)
+  )
+  target <- create_target_lookup(
+    raw, "value", "ref_date", "lag", 63, 0, 0, "weekly"
+  )
+  expect_equal(target$target_lag, 63)
+  expect_equal(target$target_date, as.Date("2024-03-09"))
+  expect_equal(target$target_type, "revision")
+})
+
+test_that("create_target_lookup returns a typed empty frame when all lags exceed ref_lag", {
+  # Simulate bulk-loaded historical data: every reference_date has only one
+  # observation with lag >> ref_lag (e.g. loaded months after the fact).
+  # Before the fix, dplyr::bind_rows(NULL, NULL, ...) returned a 0x0 tibble,
+  # causing attach_target_lookup to fail with "reference_date not in y".
+  raw <- data.frame(
+    ref_date = as.Date(c("2020-01-01", "2020-01-02", "2020-01-03")),
+    lag = c(250, 249, 248),  # all >> ref_lag=60; no fallback <= 60 exists
+    value = c(10, 12, 11)
+  )
+  result <- create_target_lookup(raw, "value", "ref_date", "lag", 60)
+  expect_s3_class(result, "data.frame")
+  expect_equal(nrow(result), 0L)
+  expect_true("reference_date" %in% colnames(result))
+  expect_true("target_date"    %in% colnames(result))
+  expect_true("target_lag"     %in% colnames(result))
+  expect_true("target_type"    %in% colnames(result))
+})
+
+test_that("data_preprocessing attaches raw-aware target after filling", {
+  raw <- data.frame(
+    ref_date = rep(as.Date("2024-01-01"), 4),
+    lag = c(50, 59, 61, 62),
+    value = c(10, 10, 12, 15)
+  )
+  result <- data_preprocessing(
+    raw, "value", "ref_date", "lag", 60,
+    lagged_term_list = c(1, 7), temporal_resol = "daily", smoothed = TRUE,
+    target_lag_lower_tolerance = 1, target_lag_upper_tolerance = 2
+  )
+  expect_true(all(result$target_lag == 62))
+  expect_true(all(result$target_type == "revision"))
+  expect_true(all(result$value_target == 15))
+})
+
 
 test_that("add_log_transformed correctly applies log transformation", {
   # Create a test dataframe
@@ -189,57 +272,73 @@ test_that("add_log_transformed correctly applies log transformation", {
 })
 
 
-test_that("add_params_for_dates correctly adds date-related features", {
-  # Create a sample data frame
+test_that("add_grouped_dayofweek creates one column per group with correct values", {
+  # Mon=1, Tue=2, ..., Sun=7 via %u format
+  df <- data.frame(
+    date = as.Date(c("2024-02-26", "2024-02-27", "2024-02-28", "2024-03-01", "2024-03-02", "2024-03-03"))
+    # Mon,        Tue,           Wed,           Fri,           Sat,           Sun
+  )
+  groups <- list(Mon = c("Mon"), Weekends = c("Sat", "Sun"), Other = c("Tue", "Wed", "Thurs", "Fri"))
+  result <- add_grouped_dayofweek(df, "date", "_ref", groups)
+
+  expect_true(all(c("Mon_ref", "Weekends_ref", "Other_ref") %in% colnames(result)))
+  expect_equal(result$Mon_ref,      c(1L, 0L, 0L, 0L, 0L, 0L))
+  expect_equal(result$Weekends_ref, c(0L, 0L, 0L, 0L, 1L, 1L))
+  expect_equal(result$Other_ref,    c(0L, 1L, 1L, 1L, 0L, 0L))
+  # every row sums to exactly 1 (exhaustive, non-overlapping groups)
+  expect_true(all(rowSums(result[, c("Mon_ref", "Weekends_ref", "Other_ref")]) == 1L))
+})
+
+test_that("add_grouped_dayofweek derives column names from day abbreviations when unnamed", {
+  df <- data.frame(date = as.Date(c("2024-02-26", "2024-03-02")))  # Mon, Sat
+  result <- add_grouped_dayofweek(df, "date", "_ref", list(c("Mon"), c("Sat", "Sun")))
+  expect_true("Mon_ref" %in% colnames(result))
+  expect_true("SatSun_ref" %in% colnames(result))
+})
+
+test_that("add_params_for_dates uses default Mon/Weekends groups in daily mode", {
   test_df <- data.frame(
     ref_date = as.Date(c("2022-01-01", "2022-01-05", "2022-01-10", "2022-02-01", "2022-02-15")),
     lag = c(0, 2, 5, 7, 10)
   )
-
-  # Run the function with daily resolution
   df_with_params <- add_params_for_dates(test_df, "ref_date", "lag", "daily")
 
-  # Check that report_date column is correctly created
   expect_true("report_date" %in% colnames(df_with_params))
-
-  # Verify day-of-week encoding is added for both reference and issue date
-  expect_true(all(paste0(WEEKDAYS_ABBR, "_ref") %in% colnames(df_with_params)))
-  expect_true(all(paste0(WEEKDAYS_ABBR, "_issue") %in% colnames(df_with_params)))
-
-  expect_true(all(c("Weekends_issue", "Weekends_ref") %in% colnames(df_with_params)))
-
-  # Verify that exactly one column per row is 1 for each set of one-hot encoded days
-  expect_true(all(rowSums(df_with_params[, paste0(WEEKDAYS_ABBR, "_ref")]) == 1))
-  expect_true(all(rowSums(df_with_params[, paste0(WEEKDAYS_ABBR, "_issue")]) == 1))
-
-  # Verify week-of-month encoding is added for report_date
+  expect_true(all(c("Mon_ref", "Weekends_ref", "Mon_issue", "Weekends_issue") %in% colnames(df_with_params)))
+  # old per-day columns should not be present
+  expect_false(any(c("Tue_ref", "Wed_ref", "Thurs_ref", "Fri_ref", "Sat_ref", "Sun_ref") %in% colnames(df_with_params)))
+  # each row is either Mon (1,0), Weekends (0,1), or Other (0,0) — never (1,1)
+  expect_true(all(df_with_params$Mon_ref + df_with_params$Weekends_ref <= 1L))
+  expect_true(all(df_with_params$Mon_issue + df_with_params$Weekends_issue <= 1L))
   expect_true(all(WEEK_ISSUES %in% colnames(df_with_params)))
-
-  # Check that only one week column per row has a value of 1
   expect_true(all(rowSums(df_with_params[, WEEK_ISSUES]) <= 1))
 })
 
+test_that("add_params_for_dates respects custom onehot_weekdays", {
+  # Dates: 2024-02-26=Mon, 2024-02-28=Wed, 2024-03-01=Fri, 2024-03-02=Sat
+  test_df <- data.frame(
+    ref_date = as.Date(c("2024-02-26", "2024-02-28", "2024-03-01", "2024-03-02")),
+    lag = c(0, 0, 0, 0)
+  )
+  groups <- list(WedFri = c("Wed", "Fri"), Other = c("Mon", "Tue", "Thurs", "Sat", "Sun"))
+  result <- add_params_for_dates(test_df, "ref_date", "lag", "daily", onehot_weekdays = groups)
+
+  expect_true(all(c("WedFri_ref", "Other_ref") %in% colnames(result)))
+  expect_equal(result$WedFri_ref, c(0L, 1L, 1L, 0L))
+  expect_equal(result$Other_ref,  c(1L, 0L, 0L, 1L))
+})
+
 test_that("add_params_for_dates correctly handles weekly resolution", {
-  # Create a sample data frame
   test_df <- data.frame(
     ref_date = as.Date(c("2022-03-01", "2022-03-08", "2022-03-15")),
     lag = c(0, 7, 14)
   )
-
-  # Run the function with weekly resolution
   df_with_params <- add_params_for_dates(test_df, "ref_date", "lag", "weekly")
 
-  # Check that report_date column is correctly created
   expect_true("report_date" %in% colnames(df_with_params))
-
-  # Verify that day-of-week encoding is NOT added in weekly mode
-  expect_false(any(paste0(WEEKDAYS_ABBR, "_ref") %in% colnames(df_with_params)))
-  expect_false(any(paste0(WEEKDAYS_ABBR, "_issue") %in% colnames(df_with_params)))
-
-  # Verify week-of-month encoding is still added for report_date
+  # no day-of-week columns in weekly mode
+  expect_false(any(c("Mon_ref", "Weekends_ref", "Mon_issue", "Weekends_issue") %in% colnames(df_with_params)))
   expect_true(all(WEEK_ISSUES %in% colnames(df_with_params)))
-
-  # Check that only one week column per row has a value of 1
   expect_true(all(rowSums(df_with_params[, WEEK_ISSUES]) <= 1))
 })
 
@@ -305,10 +404,13 @@ test_that("data_preprocessing handles multiple value columns correctly", {
   expect_true("log_delta_value_7dav_lag7" %in% colnames(result_df))
   expect_true("log_delta_value_7dav_lag7" %in% colnames(result_df))
 
-  expect_error(data_preprocessing(df, value_col = c("cases", "deaths"), suffixes=c("_num", "_denom"),
-                                  refd_col = "ref_date", lag_col = "lag", ref_lag = 7, value_type = "fraction",
-                                  temporal_resol = "weekly"),
-               "The reference dates do not regularly have a gap of 7 days. Some reference dates will be ignored. Please check your input data.")
+  weekly_result <- data_preprocessing(
+    df, value_col = c("cases", "deaths"), suffixes = c("_num", "_denom"),
+    refd_col = "ref_date", lag_col = "lag", ref_lag = 7,
+    value_type = "fraction", temporal_resol = "weekly"
+  )
+  expect_true(all(weekdays(weekly_result$reference_date) == "Saturday"))
+  expect_true(all(weekly_result$lag %% 7 == 0))
 
   expect_true(max(result_df$lag) < 7)
   expect_true("reference_date" %in% colnames(result_df))
@@ -341,3 +443,90 @@ test_that("Testing add weighted related features", {
   expect_true(all(expected_columns %in% colnames(result)))
 
 })
+
+
+make_daily_tri <- function(ref_dates, lags) {
+  do.call(rbind, lapply(ref_dates, function(rd) {
+    data.frame(
+      reference_date = rd,
+      lag = lags,
+      value = as.numeric(as.Date(rd) - as.Date("2022-12-31"))
+    )
+  }))
+}
+
+test_that("data_preprocessing drops rows where aux has no coverage (earlier min)", {
+  ref_all  <- seq(as.Date("2023-01-01"), as.Date("2023-01-15"), by = "day")
+  ref_late <- seq(as.Date("2023-01-08"), as.Date("2023-01-15"), by = "day")
+
+  primary <- make_daily_tri(ref_all, 1:3)
+  aux     <- make_daily_tri(ref_late, 1:3)
+
+  result <- data_preprocessing(
+    primary,
+    value_col = "value", refd_col = "reference_date", lag_col = "lag",
+    ref_lag = 5L, temporal_resol = "daily",
+    aux_triangles = list(beds = aux)
+  )
+
+  expect_true(all(result$reference_date >= as.Date("2023-01-08")))
+  expect_true("beds_value_raw" %in% colnames(result))
+  expect_false(anyNA(result[["beds_value_raw"]]))
+})
+
+test_that("data_preprocessing forward-fills aux to primary max report_date (later max)", {
+  ref_dates    <- seq(as.Date("2023-01-01"), as.Date("2023-01-10"), by = "day")
+  primary <- make_daily_tri(ref_dates, 1:5)
+  aux     <- make_daily_tri(ref_dates, 1:3)
+  # primary max report_date = 2023-01-10 + 5 = 2023-01-15
+  # aux    max report_date  = 2023-01-10 + 3 = 2023-01-13
+
+  result <- data_preprocessing(
+    primary,
+    value_col = "value", refd_col = "reference_date", lag_col = "lag",
+    ref_lag = 6L, temporal_resol = "daily",
+    aux_triangles = list(beds = aux)
+  )
+
+  # rows with lag 4 and 5 (report_dates past aux max) should be present
+  expect_true(any(result$lag == 4L))
+  expect_true(any(result$lag == 5L))
+  expect_true("beds_value_raw" %in% colnames(result))
+  expect_false(anyNA(result[["beds_value_raw"]]))
+})
+
+test_that("data_preprocessing returns 0 rows when aux has no data for the geo", {
+  ref_dates <- seq(as.Date("2023-01-01"), as.Date("2023-01-15"), by = "day")
+  primary   <- make_daily_tri(ref_dates, 1:4)
+  empty_aux <- data.frame(
+    reference_date = as.Date(character()),
+    report_date    = as.Date(character()),
+    lag            = integer(),
+    value          = numeric()
+  )
+
+  result <- data_preprocessing(
+    primary,
+    value_col = "value", refd_col = "reference_date", lag_col = "lag",
+    ref_lag = 5L, temporal_resol = "daily",
+    aux_triangles = list(beds = empty_aux)
+  )
+
+  expect_equal(nrow(result), 0L)
+})
+
+test_that("process_aux_triangle returns zero-row df with correct columns when aux is empty", {
+  empty_aux <- data.frame(
+    reference_date = as.Date(character()),
+    report_date    = as.Date(character()),
+    lag            = integer(),
+    value          = numeric()
+  )
+  lagged_term_list <- c(1L, 7L)
+
+  result <- process_aux_triangle(empty_aux, "beds", lagged_term_list, "daily", TRUE)
+
+  expect_equal(nrow(result), 0L)
+  expect_true(all(aux_feature_names("beds", lagged_term_list) %in% colnames(result)))
+  expect_true(all(c("reference_date", "report_date", "lag") %in% colnames(result)))
+})
diff --git a/tests/testthat/test-model.R b/tests/testthat/test-model.R
index 94727f4..381d005 100644
--- a/tests/testthat/test-model.R
+++ b/tests/testthat/test-model.R
@@ -307,3 +307,60 @@ test_that("testing data_filteration", {
   expect_equal(result$lag, expected_lags)
 
 })
+
+test_that("revision_forecast with train_models=FALSE loads cached model and produces identical predictions", {
+  tmpdir <- tempfile()
+  dir.create(tmpdir)
+  on.exit(unlink(tmpdir, recursive = TRUE))
+
+  set.seed(42)
+  nn_train <- 200
+  nn_test  <- 30
+
+  make_rf_data <- function(nn, start_date) {
+    dates <- seq(as.Date(start_date), by = "day", length.out = nn)
+    data.frame(
+      reference_date   = dates,
+      report_date      = dates + 3L,
+      lag              = 3L,
+      value_7dav       = runif(nn, 0, 1),
+      log_value_7dav   = rnorm(nn),
+      log_value_target = rnorm(nn),
+      value_7dav_diff  = rnorm(nn),
+      value_slope_diff = rnorm(nn),
+      Mon_ref          = sample(c(0L, 1L), nn, replace = TRUE)
+    )
+  }
+
+  train_data <- make_rf_data(nn_train, "2021-01-01")
+  test_data  <- make_rf_data(nn_test,  "2021-07-20")
+
+  rf_args <- list(
+    taus             = 0.5,
+    smoothed_target  = FALSE,
+    params_list      = c("log_value_7dav", "Mon_ref"),
+    temporal_resol   = "daily",
+    lambda           = 0.1,
+    gamma            = 0.1,
+    model_save_dir   = tmpdir,
+    indicator        = "test",
+    signal           = "sig",
+    geo_level        = "state",
+    geo              = "pa",
+    training_days    = nn_train,
+    make_predictions = TRUE
+  )
+
+  result_fit <- do.call(
+    revision_forecast,
+    c(list(train_data = train_data, test_data = test_data, train_models = TRUE), rf_args)
+  )
+  expect_gt(nrow(result_fit), 0)
+
+  result_cached <- do.call(
+    revision_forecast,
+    c(list(train_data = train_data, test_data = test_data, train_models = FALSE), rf_args)
+  )
+  expect_gt(nrow(result_cached), 0)
+  expect_equal(result_fit$predicted_tau0.5, result_cached$predicted_tau0.5)
+})
diff --git a/tests/testthat/test-preprocessing.R b/tests/testthat/test-preprocessing.R
index 00b380f..91d3ec6 100644
--- a/tests/testthat/test-preprocessing.R
+++ b/tests/testthat/test-preprocessing.R
@@ -63,28 +63,32 @@ test_that("fill_missing_updates correctly processes weekly data", {
   expect_equal(nrow(filled_df), 2 + 3)  # Two rows expected since it's weekly data
 })
 
-test_that("fill_missing_updates raises an error for irregular gaps", {
+test_that("weekly data are normalized to epiweek-ending Saturdays", {
   df <- data.frame(
     ref_date = as.Date(c("2023-01-01", "2023-01-05")),
     lag = c(0, 4),
     value = c(10, 20)
   )
 
-  expect_error(
-    fill_missing_updates(df, "value", "ref_date", "lag", "weekly"),
-    "The reference dates do not regularly have a gap of 7 days. Some reference dates will be ignored. Please check your input data."
-  )  # Function should stop due to irregular gap
+  filled_df <- fill_missing_updates(df, "value", "ref_date", "lag", "weekly")
+  expect_true(all(weekdays(filled_df$reference_date) == "Saturday"))
+  expect_true(all(weekdays(filled_df$report_date) == "Saturday"))
+  expect_true(all(filled_df$lag %% 7 == 0))
+})
 
+test_that("weekly data keep the latest revision in each epiweek", {
   df <- data.frame(
-    ref_date = as.Date(c("2023-01-01", "2023-01-08")),
-    lag = c(0, 4),
+    ref_date = as.Date(c("2023-01-07", "2023-01-07")),
+    lag = c(4, 6), # Wednesday and Friday of the following epiweek
     value = c(10, 20)
   )
 
-  expect_error(
-    fill_missing_updates(df, "value", "ref_date", "lag", "weekly"),
-    "The report dates do not regularly have a gap of 7 days. Some report dates will be ignored. Please check your input data."
-  )  # Function should stop due to irregular gap
+  filled_df <- fill_missing_updates(df, "value", "ref_date", "lag", "weekly")
+  expect_equal(nrow(filled_df), 1)
+  expect_equal(filled_df$reference_date, as.Date("2023-01-07"))
+  expect_equal(filled_df$report_date, as.Date("2023-01-14"))
+  expect_equal(filled_df$lag, 7)
+  expect_equal(filled_df$value_raw, 20)
 })
 
 
@@ -173,5 +177,3 @@ test_that("testing the calculation of week of a month", {
 
 })
 
-
-