graph_pipeline.R

suppressPackageStartupMessages({
  library(dplyr)
  library(forcats)
  library(fs)
  library(ggplot2)
  library(ggrepel)
  library(here)
  library(lubridate)
  library(patchwork)
  library(readr)
  library(scales)
  library(slider)
  library(stringr)
  library(tibble)
  library(tidyr)
})

model_palette <- c(
  annual_seasonal_naive = wmata_colors[["light_grey"]],
  lag_7 = wmata_colors[["dark_grey"]],
  linear = wmata_colors[["turquoise"]],
  glmnet = wmata_colors[["true_blue"]],
  xgboost = wmata_colors[["dodger_blue"]],
  fallback_hierarchy = wmata_colors[["dark_blue"]],
  station_aggregated = wmata_colors[["green"]]
)

mode_palette <- c(
  Bus = wmata_colors[["dodger_blue"]],
  Rail = wmata_colors[["green"]]
)

label_model_name <- function(model_name) {
  dplyr::recode(
    model_name,
    annual_seasonal_naive = "Annual seasonal naive",
    lag_7 = "Lag 7 baseline",
    linear = "Linear regression",
    glmnet = "GLMNET",
    xgboost = "XGBoost",
    fallback_hierarchy = "Fallback hierarchy",
    station_aggregated = "Station-aggregated",
    .default = model_name
  )
}

filter_visual_models <- function(model_summary, metric_col = "mae", multiplier = 20) {
  metric_values <- model_summary[[metric_col]]
  finite_values <- metric_values[is.finite(metric_values) & metric_values > 0]

  if (length(finite_values) == 0) {
    return(model_summary %>% mutate(visual_keep = TRUE))
  }

  threshold <- min(finite_values, na.rm = TRUE) * multiplier

  model_summary %>%
    mutate(
      visual_keep = is.finite(.data[[metric_col]]) & .data[[metric_col]] <= threshold
    )
}

feature_group <- function(feature_name) {
  case_when(
    str_detect(feature_name, "^lag_|^roll_mean_|same_weekday") ~ "Recent ridership",
    str_detect(feature_name, "day_of_week|is_weekend|weekday_saturday_sunday|service_type|holiday") ~ "Calendar structure",
    str_detect(feature_name, "month|week_of_year|week_number") ~ "Seasonality",
    str_detect(feature_name, "time_index|year|station_age_days|new_station_flag") ~ "Trend and age",
    str_detect(feature_name, "^station_name") ~ "Station identity",
    TRUE ~ "Other"
  )
}

is_winter_date <- function(date) {
  lubridate::month(date) %in% c(12L, 1L, 2L)
}

winter_season_label <- function(date) {
  if_else(
    lubridate::month(date) == 12L,
    paste0(lubridate::year(date), "-", lubridate::year(date) + 1L),
    paste0(lubridate::year(date) - 1L, "-", lubridate::year(date))
  )
}

format_day_name <- function(date) {
  factor(
    lubridate::wday(date, label = TRUE, abbr = TRUE, week_start = 1),
    levels = lubridate::wday(seq(as.Date("2026-01-05"), by = "day", length.out = 7), label = TRUE, abbr = TRUE, week_start = 1)
  )
}

build_rail_system_history <- function(silver_data) {
  silver_data$ridership_rail_station_daily %>%
    group_by(date) %>%
    summarise(
      ridership = sum(ridership, na.rm = TRUE),
      holiday = dplyr::first(holiday),
      day_of_week = dplyr::first(day_of_week),
      weekday_saturday_sunday = dplyr::first(weekday_saturday_sunday),
      service_type = dplyr::first(service_type),
      .groups = "drop"
    ) %>%
    arrange(date)
}

build_rail_future_context <- function(rail_results) {
  bind_rows(
    rail_results$main_results$final_future_forecasts,
    rail_results$fallback_results$future_forecasts
  ) %>%
    group_by(date) %>%
    summarise(
      day_of_week = dplyr::first(day_of_week),
      weekday_saturday_sunday = dplyr::first(weekday_saturday_sunday),
      service_type = dplyr::first(service_type),
      .groups = "drop"
    ) %>%
    arrange(date)
}

selected_backtest_predictions <- function(bus_results, rail_results, bus_data, rail_system_history) {
  bind_rows(
    bus_results$backtest_predictions %>%
      filter(model == bus_results$selected_model_name) %>%
      left_join(bus_data %>% select(date, holiday, day_of_week, weekday_saturday_sunday, service_type), by = "date") %>%
      mutate(mode = "Bus"),
    rail_results$system_backtest_predictions %>%
      filter(model == rail_results$main_results$selected_model_name) %>%
      left_join(rail_system_history %>% select(date, holiday, day_of_week, weekday_saturday_sunday, service_type), by = "date") %>%
      mutate(mode = "Rail")
  )
}

selected_holdout_predictions <- function(bus_results, rail_results, bus_data, rail_system_history) {
  bind_rows(
    bus_results$holdout_predictions %>%
      left_join(bus_data %>% select(date, holiday, day_of_week, weekday_saturday_sunday, service_type), by = "date") %>%
      mutate(mode = "Bus"),
    rail_results$system_holdout_predictions %>%
      left_join(rail_system_history %>% select(date, holiday, day_of_week, weekday_saturday_sunday, service_type), by = "date") %>%
      mutate(mode = "Rail")
  )
}

build_future_forecasts_by_mode <- function(bus_results, rail_results) {
  rail_future_context <- build_rail_future_context(rail_results)

  bind_rows(
    add_prediction_intervals(
      bus_results$final_future_forecasts %>% select(date, prediction),
      bus_results$backtest_predictions %>% filter(model == bus_results$selected_model_name)
    ) %>%
      left_join(
        bus_results$final_future_forecasts %>% select(date, day_of_week, weekday_saturday_sunday, service_type),
        by = "date"
      ) %>%
      mutate(mode = "Bus"),
    add_prediction_intervals(
      rail_results$system_future_forecasts,
      rail_results$system_backtest_predictions %>% filter(model == rail_results$main_results$selected_model_name)
    ) %>%
      left_join(rail_future_context, by = "date") %>%
      mutate(mode = "Rail")
  )
}

compute_coverage_summary <- function(predictions) {
  bind_rows(
    predictions %>%
      group_by(mode) %>%
      summarise(
        interval = "80%",
        nominal_coverage = 0.80,
        observed_coverage = mean(actual >= lower_80 & actual <= upper_80, na.rm = TRUE),
        average_width = mean(upper_80 - lower_80, na.rm = TRUE),
        .groups = "drop"
      ),
    predictions %>%
      group_by(mode) %>%
      summarise(
        interval = "95%",
        nominal_coverage = 0.95,
        observed_coverage = mean(actual >= lower_95 & actual <= upper_95, na.rm = TRUE),
        average_width = mean(upper_95 - lower_95, na.rm = TRUE),
        .groups = "drop"
      )
  )
}

compute_success_rate <- function(actual, prediction, tolerance = 0.10) {
  mean(abs(prediction - actual) / pmax(actual, 1) <= tolerance, na.rm = TRUE)
}

build_model_leaderboard <- function(backtest_predictions) {
  backtest_predictions %>%
    filter(is.finite(actual), is.finite(prediction)) %>%
    group_by(model) %>%
    summarise(
      mae = mean(abs(prediction - actual), na.rm = TRUE),
      rmse = sqrt(mean((prediction - actual) ^ 2, na.rm = TRUE)),
      mape = 100 * mean(abs(prediction - actual) / pmax(actual, 1), na.rm = TRUE),
      success_rate_10 = 100 * compute_success_rate(actual, prediction, tolerance = 0.10),
      .groups = "drop"
    ) %>%
    filter_visual_models(metric_col = "mae") %>%
    filter(visual_keep) %>%
    mutate(
      model_label = label_model_name(model)
    )
}

build_feature_importance_groups <- function(feature_importance) {
  if (nrow(feature_importance) == 0) {
    return(tibble(feature_group = character(), importance = numeric()))
  }

  feature_column <- if ("term" %in% names(feature_importance)) "term" else if ("Variable" %in% names(feature_importance)) "Variable" else names(feature_importance)[[1]]
  importance_column <- if ("importance" %in% names(feature_importance)) "importance" else if ("Importance" %in% names(feature_importance)) "Importance" else names(feature_importance)[[2]]

  feature_importance %>%
    transmute(
      feature = .data[[feature_column]],
      importance = abs(.data[[importance_column]])
    ) %>%
    filter(!is.na(feature), !is.na(importance)) %>%
    mutate(feature_group = feature_group(feature)) %>%
    group_by(feature_group) %>%
    summarise(importance = sum(importance, na.rm = TRUE), .groups = "drop") %>%
    arrange(desc(importance))
}

build_horizon_decay <- function(predictions) {
  predictions %>%
    filter(horizon_day >= 1, horizon_day <= 30, is.finite(actual), is.finite(prediction)) %>%
    group_by(mode, horizon_day) %>%
    summarise(
      mae = mean(abs(prediction - actual), na.rm = TRUE),
      .groups = "drop"
    )
}

build_error_by_month <- function(predictions) {
  predictions %>%
    mutate(month_label = factor(month(date, label = TRUE, abbr = TRUE), levels = month.abb)) %>%
    group_by(mode, month_label) %>%
    summarise(
      mae = mean(abs(prediction - actual), na.rm = TRUE),
      .groups = "drop"
    )
}

build_error_by_ridership_level <- function(predictions) {
  predictions %>%
    group_by(mode) %>%
    mutate(
      ridership_band = ntile(actual, 4),
      ridership_band = factor(ridership_band, labels = c("Lowest quartile", "Lower-middle", "Upper-middle", "Highest quartile"))
    ) %>%
    group_by(mode, ridership_band) %>%
    summarise(
      mae = mean(abs(prediction - actual), na.rm = TRUE),
      .groups = "drop"
    )
}

build_calibration_curve <- function(predictions) {
  predictions %>%
    group_by(mode) %>%
    mutate(calibration_bin = ntile(prediction, 10)) %>%
    group_by(mode, calibration_bin) %>%
    summarise(
      avg_prediction = mean(prediction, na.rm = TRUE),
      avg_actual = mean(actual, na.rm = TRUE),
      .groups = "drop"
    )
}

build_directional_accuracy <- function(predictions) {
  predictions %>%
    arrange(mode, date) %>%
    group_by(mode) %>%
    mutate(
      previous_actual = lag(actual),
      actual_direction = sign(actual - previous_actual),
      predicted_direction = sign(prediction - previous_actual)
    ) %>%
    filter(!is.na(previous_actual)) %>%
    summarise(
      directional_accuracy = mean(actual_direction == predicted_direction, na.rm = TRUE),
      .groups = "drop"
    )
}

build_accuracy_traffic_light <- function(predictions) {
  predictions %>%
    mutate(horizon_bucket = horizon_bucket(horizon_day)) %>%
    group_by(mode, horizon_bucket) %>%
    summarise(
      mape = 100 * mean(abs(prediction - actual) / pmax(actual, 1), na.rm = TRUE),
      .groups = "drop"
    ) %>%
    mutate(
      status = case_when(
        mape <= 10 ~ "Green",
        mape <= 15 ~ "Yellow",
        TRUE ~ "Red"
      )
    )
}

build_worst_prediction_days <- function(predictions, top_n = 5) {
  predictions %>%
    mutate(
      abs_error = abs(prediction - actual)
    ) %>%
    group_by(mode) %>%
    slice_max(abs_error, n = top_n, with_ties = FALSE) %>%
    mutate(
      error = prediction - actual,
      error_direction = if_else(error >= 0, "Overpredicted", "Underpredicted"),
      likely_explanation = case_when(
        !is.na(holiday) & holiday != "No" ~ "Holiday or service-calendar effect",
        !is.na(weekday_saturday_sunday) & weekday_saturday_sunday != "Weekday" ~ "Weekend demand shift",
        abs_error / pmax(actual, 1) >= 0.20 ~ "Possible disruption or special event",
        TRUE ~ "Needs manual review"
      )
    ) %>%
    ungroup()
}

build_success_rate_summary <- function(predictions) {
  prediction_splits <- predictions %>%
    group_by(mode) %>%
    summarise(data = list(pick(everything())), .groups = "drop")

  crossing(
    mode = unique(predictions$mode),
    threshold = c(0.05, 0.10, 0.15)
  ) %>%
    left_join(prediction_splits, by = "mode") %>%
    rowwise() %>%
    mutate(
      success_rate = compute_success_rate(data$actual, data$prediction, tolerance = threshold)
    ) %>%
    ungroup() %>%
    transmute(
      mode,
      threshold_label = paste0("Within ±", round(threshold * 100), "%"),
      success_rate
    )
}

build_smoothed_horizon_decay <- function(predictions) {
  build_horizon_decay(predictions) %>%
    arrange(mode, horizon_day) %>%
    group_by(mode) %>%
    mutate(
      smoothed_mae = slider::slide_dbl(mae, mean, .before = 3, .after = 3, .complete = FALSE, na.rm = TRUE)
    ) %>%
    ungroup()
}

build_decision_matrix <- function(future_forecasts) {
  future_forecasts %>%
    mutate(
      interval_width = upper_80 - lower_80,
      day_type = if_else(weekday_saturday_sunday == "Weekday", "Weekday", "Weekend")
    ) %>%
    group_by(mode) %>%
    mutate(
      demand_midpoint = median(prediction, na.rm = TRUE),
      uncertainty_midpoint = median(interval_width, na.rm = TRUE),
      action_bucket = case_when(
        prediction >= demand_midpoint & interval_width < uncertainty_midpoint ~ "Plan capacity",
        prediction >= demand_midpoint & interval_width >= uncertainty_midpoint ~ "Prepare contingency",
        prediction < demand_midpoint & interval_width < uncertainty_midpoint ~ "Optimize resources",
        TRUE ~ "Manual review"
      )
    ) %>%
    ungroup()
}

build_forecast_extremes <- function(future_forecasts) {
  bind_rows(
    future_forecasts %>%
      group_by(mode) %>%
      slice_max(prediction, n = 1, with_ties = FALSE) %>%
      mutate(extreme_type = "Highest forecast day"),
    future_forecasts %>%
      group_by(mode) %>%
      slice_min(prediction, n = 1, with_ties = FALSE) %>%
      mutate(extreme_type = "Lowest forecast day")
  ) %>%
    ungroup()
}

build_recent_growth_table <- function(rail_results, silver_data) {
  recent_actual <- silver_data$ridership_rail_station_daily %>%
    filter(date >= max(date) - 27) %>%
    group_by(station_name) %>%
    summarise(recent_avg_ridership = mean(ridership, na.rm = TRUE), .groups = "drop")

  future_station_forecasts <- bind_rows(
    rail_results$main_results$final_future_forecasts %>% mutate(forecast_group = "main"),
    rail_results$fallback_results$future_forecasts %>% mutate(forecast_group = "fallback")
  ) %>%
    group_by(station_name) %>%
    summarise(
      forecast_avg_ridership = mean(prediction, na.rm = TRUE),
      .groups = "drop"
    )

  future_station_forecasts %>%
    left_join(recent_actual, by = "station_name") %>%
    mutate(
      recent_avg_ridership = coalesce(recent_avg_ridership, 0),
      absolute_growth = forecast_avg_ridership - recent_avg_ridership,
      pct_growth = if_else(recent_avg_ridership > 0, absolute_growth / recent_avg_ridership, NA_real_)
    ) %>%
    arrange(desc(absolute_growth))
}

build_rail_station_risk <- function(rail_results, silver_data) {
  station_growth <- build_recent_growth_table(rail_results, silver_data)

  station_accuracy <- rail_results$main_results$holdout_predictions %>%
    group_by(station_name) %>%
    summarise(
      holdout_mae = mean(abs(prediction - actual), na.rm = TRUE),
      .groups = "drop"
    )

  station_growth %>%
    left_join(station_accuracy, by = "station_name") %>%
    filter(is.finite(holdout_mae))
}

build_data_freshness <- function(bus_data, rail_system_history) {
  tibble(
    mode = c("Bus", "Rail"),
    latest_actual_date = c(max(bus_data$date, na.rm = TRUE), max(rail_system_history$date, na.rm = TRUE))
  ) %>%
    mutate(
      report_date = Sys.Date(),
      days_old = as.integer(report_date - latest_actual_date)
    )
}

build_ridership_history <- function(bus_data, rail_system_history) {
  bind_rows(
    bus_data %>%
      select(date, ridership) %>%
      mutate(mode = "Bus"),
    rail_system_history %>%
      select(date, ridership) %>%
      mutate(mode = "Rail")
  ) %>%
    arrange(mode, date) %>%
    group_by(mode) %>%
    mutate(
      rolling_28 = slider::slide_dbl(ridership, mean, .before = 27, .complete = FALSE, na.rm = TRUE)
    ) %>%
    ungroup()
}

build_recovery_index <- function(history_data) {
  baseline <- history_data %>%
    filter(date >= as.Date("2019-01-01"), date <= as.Date("2019-12-31")) %>%
    group_by(mode) %>%
    summarise(baseline_ridership = mean(ridership, na.rm = TRUE), .groups = "drop")

  history_data %>%
    mutate(month_date = floor_date(date, unit = "month")) %>%
    group_by(mode, month_date) %>%
    summarise(monthly_ridership = mean(ridership, na.rm = TRUE), .groups = "drop") %>%
    left_join(baseline, by = "mode") %>%
    mutate(index_to_2019 = monthly_ridership / baseline_ridership)
}

plot_actual_vs_predicted <- function(predictions, title_text, subtitle_text) {
  interval_layer <- if (all(c("lower_80", "upper_80") %in% names(predictions))) {
    geom_ribbon(
      data = predictions,
      aes(x = date, ymin = lower_80, ymax = upper_80),
      inherit.aes = FALSE,
      fill = wmata_colors[["sky_blue"]],
      alpha = 0.35
    )
  } else {
    NULL
  }

  ggplot(predictions, aes(x = date)) +
    interval_layer +
    geom_line(aes(y = actual, color = "Actual"), linewidth = 0.95) +
    geom_line(aes(y = prediction, color = "Predicted"), linewidth = 0.95, linetype = "dashed") +
    scale_color_manual(values = c(Actual = wmata_colors[["midnight_blue"]], Predicted = wmata_colors[["dodger_blue"]])) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = subtitle_text,
      x = NULL,
      y = "Daily Ridership",
      color = NULL,
      caption = "What this shows: held-out actual ridership versus model predictions and the 80% planning band when available | Why it matters: it is the clearest out-of-sample reliability view | What WMATA should do: use it as the main credibility check before operational use"
    ) +
    theme_wmata()
}

plot_decision_workflow <- function(title_text) {
  workflow_steps <- tibble(
    step = factor(
      c("Ridership history", "Forecast model", "30-day outlook", "WMATA actions"),
      levels = c("Ridership history", "Forecast model", "30-day outlook", "WMATA actions")
    ),
    x = c(1, 2, 3, 4),
    y = c(1, 1, 1, 1),
    detail = c(
      "Daily Bus and Rail demand history",
      "Forecast-safe features\nand shared production pipeline",
      "Expected demand,\npeaks, troughs, uncertainty",
      "Staffing, service assumptions,\nmanual review, communications"
    )
  )

  arrows_tbl <- tibble(
    x = c(1.2, 2.2, 3.2),
    xend = c(1.8, 2.8, 3.8),
    y = 1,
    yend = 1
  )

  ggplot() +
    geom_segment(data = arrows_tbl, aes(x = x, xend = xend, y = y, yend = yend), arrow = arrow(length = grid::unit(0.15, "inches")), linewidth = 0.9, color = wmata_colors[["dark_grey"]]) +
    geom_label(
      data = workflow_steps,
      aes(x = x, y = y, label = step),
      fill = wmata_colors[["midnight_blue"]],
      color = "white",
      fontface = "bold",
      size = 4.2,
      linewidth = 0
    ) +
    geom_text(data = workflow_steps, aes(x = x, y = y - 0.2, label = detail), color = wmata_colors[["dark_blue"]], size = 3.5, lineheight = 1.1, vjust = 1) +
    coord_cartesian(xlim = c(0.5, 4.5), ylim = c(0.3, 1.25), clip = "off") +
    labs(
      title = title_text,
      subtitle = "The forecast is meant to help planners prioritize attention, not automate service changes by itself"
    ) +
    theme_wmata() +
    theme(
      axis.text = element_blank(),
      axis.title = element_blank(),
      axis.ticks = element_blank(),
      panel.grid = element_blank()
    )
}

plot_pattern_capture <- function(predictions, title_text) {
  window_n <- min(42L, nrow(predictions))

  window_data <- predictions %>%
    arrange(date) %>%
    slice_tail(n = window_n) %>%
    mutate(
      abs_error = abs(prediction - actual),
      label_type = case_when(
        date == date[which.max(if_else(weekday_saturday_sunday == "Weekday", actual, -Inf))] ~ "Weekday peak",
        date == date[which.min(if_else(weekday_saturday_sunday != "Weekday", actual, Inf))] ~ "Weekend trough",
        date == date[which.max(abs_error)] ~ "Largest miss",
        TRUE ~ NA_character_
      )
    )

  label_points <- window_data %>% filter(!is.na(label_type))

  ggplot(window_data, aes(x = date)) +
    geom_ribbon(aes(ymin = lower_80, ymax = upper_80), fill = wmata_colors[["sky_blue"]], alpha = 0.28) +
    geom_line(aes(y = actual, color = "Actual"), linewidth = 0.95) +
    geom_line(aes(y = prediction, color = "Predicted"), linewidth = 0.95, linetype = "dashed") +
    geom_point(data = label_points, aes(y = actual), color = wmata_colors[["red"]], size = 2.4) +
    geom_text_repel(
      data = label_points,
      aes(y = actual, label = label_type),
      size = 3.3,
      color = wmata_colors[["midnight_blue"]],
      nudge_y = 0.04 * max(window_data$actual, na.rm = TRUE)
    ) +
    scale_color_manual(values = c(Actual = wmata_colors[["midnight_blue"]], Predicted = wmata_colors[["dodger_blue"]])) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "A shorter 6-week window makes the weekly operating rhythm and the most visible miss easier to interpret",
      x = NULL,
      y = "Daily Ridership",
      color = NULL
    ) +
    theme_wmata()
}

plot_forecast_trend <- function(history_data, forecast_data, title_text) {
  history_tail <- history_data %>% filter(date >= max(date) - 60)

  ggplot() +
    geom_line(data = history_tail, aes(x = date, y = ridership), color = wmata_colors[["dark_grey"]], linewidth = 0.8) +
    geom_ribbon(
      data = forecast_data,
      aes(x = date, ymin = lower_80, ymax = upper_80),
      fill = wmata_colors[["sky_blue"]],
      alpha = 0.35
    ) +
    geom_line(data = forecast_data, aes(x = date, y = prediction), color = wmata_colors[["green"]], linewidth = 1) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Recent actuals are shown in gray and the 30-day production forecast is shown in green",
      x = NULL,
      y = "Daily Ridership",
      caption = "What this shows: the near-term planning outlook with empirical uncertainty | Why it matters: helps planners distinguish expected demand from upside risk | What WMATA should do: use it for staffing, service planning, and rider communications"
    ) +
    theme_wmata()
}

plot_residual_distribution <- function(predictions, title_text) {
  residual_data <- predictions %>%
    mutate(
      residual = prediction - actual,
      ape = abs(prediction - actual) / pmax(actual, 1)
    )

  stats_label <- glue::glue(
    "Mean: {scales::comma(round(mean(residual_data$residual, na.rm = TRUE)))}\n",
    "Median: {scales::comma(round(median(residual_data$residual, na.rm = TRUE)))}\n",
    "Within 10%: {scales::percent(mean(residual_data$ape <= 0.10, na.rm = TRUE), accuracy = 1)}\n",
    "Within 15%: {scales::percent(mean(residual_data$ape <= 0.15, na.rm = TRUE), accuracy = 1)}"
  )

  ggplot(residual_data, aes(x = residual)) +
    geom_histogram(fill = wmata_colors[["dodger_blue"]], color = "white", bins = 30) +
    geom_vline(xintercept = 0, color = wmata_colors[["red"]], linewidth = 0.8, linetype = "dashed") +
    annotate("label", x = Inf, y = Inf, label = stats_label, hjust = 1.02, vjust = 1.1, size = 3.2, fill = "white") +
    scale_x_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "The summary callout converts the residual shape into an operational performance statement",
      x = "Prediction Error",
      y = "Count",
      caption = "What this shows: the spread and bias of remaining error | Why it matters: systematic drift creates planning risk | What WMATA should do: watch for widening residuals before reusing the model operationally"
    ) +
    theme_wmata()
}

plot_error_by_horizon <- function(horizon_summary, title_text) {
  horizon_visual <- horizon_summary %>%
    filter_visual_models(metric_col = "mae") %>%
    filter(visual_keep) %>%
    mutate(model_label = label_model_name(model))

  ggplot(horizon_visual, aes(x = horizon_bucket, y = mae, color = model, group = model)) +
    geom_line(linewidth = 0.9) +
    geom_point(size = 2.4) +
    scale_color_manual(values = model_palette, labels = label_model_name) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Only stable candidate models are shown to avoid broken scale distortion",
      x = NULL,
      y = "MAE",
      color = NULL
    ) +
    theme_wmata()
}

plot_error_by_service_type <- function(predictions, title_text) {
  service_summary <- predictions %>%
    mutate(service_group = coalesce(as.character(service_type), as.character(weekday_saturday_sunday), "Unknown")) %>%
    group_by(service_group) %>%
    summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), .groups = "drop")

  ggplot(service_summary, aes(x = reorder(service_group, mae), y = mae)) +
    geom_col(fill = wmata_colors[["true_blue"]]) +
    coord_flip() +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Segments with higher MAE may need separate monitoring",
      x = NULL,
      y = "MAE"
    ) +
    theme_wmata()
}

plot_station_error_rank <- function(predictions, title_text) {
  station_summary <- predictions %>%
    group_by(station_name) %>%
    summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), .groups = "drop") %>%
    arrange(desc(mae)) %>%
    slice_head(n = 20)

  ggplot(station_summary, aes(x = mae, y = reorder(station_name, mae))) +
    geom_col(fill = wmata_colors[["dark_blue"]]) +
    geom_text_repel(aes(label = scales::comma(round(mae))), direction = "y", size = 3.1, color = wmata_colors[["midnight_blue"]]) +
    scale_x_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Highest-error stations in the holdout window",
      x = "MAE",
      y = NULL
    ) +
    theme_wmata()
}

plot_model_leaderboard <- function(leaderboard_data, title_text) {
  metric_labels <- c(
    mae = "MAE",
    rmse = "RMSE",
    mape = "MAPE",
    success_rate_10 = "Share Within 10%"
  )

  leaderboard_long <- leaderboard_data %>%
    select(model, model_label, mae, rmse, mape, success_rate_10) %>%
    pivot_longer(cols = c(mae, rmse, mape, success_rate_10), names_to = "metric", values_to = "value") %>%
    mutate(
      metric = factor(metric, levels = names(metric_labels), labels = unname(metric_labels))
    )

  ggplot(leaderboard_long, aes(x = value, y = fct_reorder(model_label, value, .fun = min), color = model)) +
    geom_point(size = 3) +
    facet_wrap(~ metric, scales = "free_x", ncol = 2) +
    scale_color_manual(values = model_palette, guide = "none") +
    scale_x_continuous(labels = scales::label_number(accuracy = 0.1, big.mark = ",")) +
    labs(
      title = title_text,
      subtitle = "The leaderboard removes visually broken candidates so the meaningful comparisons stay readable",
      x = NULL,
      y = NULL
    ) +
    theme_wmata()
}

plot_feature_importance <- function(feature_importance, title_text) {
  feature_groups <- build_feature_importance_groups(feature_importance)

  if (nrow(feature_groups) == 0) {
    return(
      ggplot() +
        annotate("text", x = 0, y = 0, label = "Feature importance unavailable for this run", color = wmata_colors[["dark_grey"]], size = 5) +
        labs(title = title_text, subtitle = "No feature importance output was available") +
        theme_wmata() +
        theme(axis.text = element_blank(), axis.title = element_blank(), panel.grid = element_blank())
    )
  }

  feature_groups <- feature_groups %>%
    mutate(share = importance / sum(importance, na.rm = TRUE))

  ggplot(feature_groups, aes(x = share, y = fct_reorder(feature_group, share))) +
    geom_col(fill = wmata_colors[["green"]]) +
    geom_text(aes(label = scales::percent(share, accuracy = 1)), hjust = -0.1, size = 3.3) +
    scale_x_continuous(labels = scales::percent, expand = expansion(mult = c(0, 0.12))) +
    labs(
      title = title_text,
      subtitle = "The grouped view makes the model story legible: recent ridership and calendar timing drive most of the forecast movement",
      x = "Share of grouped importance",
      y = NULL,
      caption = "What this shows: the model is driven mostly by recent ridership and calendar structure, not opaque black-box effects | Why it matters: easier auditability and easier stakeholder explanation | What WMATA should do: use this to explain why the forecast moves on specific days"
    ) +
    theme_wmata()
}

plot_winter_actual_vs_predicted <- function(predictions, title_text) {
  winter_predictions <- predictions %>%
    filter(is_winter_date(date)) %>%
    mutate(
      winter_season = winter_season_label(date)
    )

  ggplot(winter_predictions, aes(x = date)) +
    geom_line(aes(y = actual, color = "Actual"), linewidth = 0.9) +
    geom_line(aes(y = prediction, color = "Predicted"), linewidth = 0.9, linetype = "dashed") +
    facet_grid(mode ~ winter_season, scales = "free_x") +
    scale_color_manual(values = c(Actual = wmata_colors[["midnight_blue"]], Predicted = wmata_colors[["dodger_blue"]])) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Winter months are isolated explicitly so stakeholders can judge cold-season performance rather than annual averages alone",
      x = NULL,
      y = "Daily Ridership",
      color = NULL
    ) +
    theme_wmata(base_size = 13)
}

plot_winter_mae_comparison <- function(predictions, title_text) {
  winter_summary <- predictions %>%
    mutate(season_group = if_else(is_winter_date(date), "Winter months", "Non-winter months")) %>%
    group_by(mode, season_group) %>%
    summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), .groups = "drop")

  ggplot(winter_summary, aes(x = season_group, y = mae, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    geom_text(aes(label = scales::comma(round(mae))), position = position_dodge(width = 0.7), vjust = -0.3, size = 3.2) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.12))) +
    labs(
      title = title_text,
      subtitle = "This directly answers whether winter is materially harder than the rest of the year",
      x = NULL,
      y = "MAE",
      fill = NULL
    ) +
    theme_wmata()
}

plot_winter_weekday_pattern <- function(predictions, title_text) {
  weekday_pattern <- predictions %>%
    filter(is_winter_date(date)) %>%
    mutate(day_name = format_day_name(date)) %>%
    group_by(mode, day_name) %>%
    summarise(
      actual = mean(actual, na.rm = TRUE),
      prediction = mean(prediction, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    pivot_longer(cols = c(actual, prediction), names_to = "series", values_to = "ridership") %>%
    mutate(series = recode(series, actual = "Actual", prediction = "Predicted"))

  ggplot(weekday_pattern, aes(x = day_name, y = ridership, color = series, group = series)) +
    geom_line(linewidth = 0.9) +
    geom_point(size = 2.1) +
    facet_wrap(~ mode, scales = "free_y") +
    scale_color_manual(values = c(Actual = wmata_colors[["midnight_blue"]], Predicted = wmata_colors[["dodger_blue"]])) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Average winter weekday and weekend structure is a useful sanity check for service planning relevance",
      x = NULL,
      y = "Average Daily Ridership",
      color = NULL
    ) +
    theme_wmata()
}

plot_winter_coverage <- function(predictions, title_text) {
  coverage_tbl <- predictions %>%
    filter(is_winter_date(date)) %>%
    compute_coverage_summary()

  ggplot(coverage_tbl, aes(x = interval, y = observed_coverage, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    geom_point(aes(y = nominal_coverage), position = position_dodge(width = 0.7), color = wmata_colors[["red"]], size = 2.5) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
    labs(
      title = title_text,
      subtitle = "Bars show observed winter coverage and red points show the nominal interval target",
      x = NULL,
      y = "Coverage",
      fill = NULL
    ) +
    theme_wmata()
}

plot_backtest_mae_over_time <- function(selected_backtests, title_text) {
  fold_summary <- selected_backtests %>%
    group_by(mode, origin_date) %>%
    summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), .groups = "drop") %>%
    group_by(mode) %>%
    mutate(
      avg_mae = mean(mae, na.rm = TRUE),
      season_group = if_else(is_winter_date(origin_date), "Winter fold", "Non-winter fold")
    ) %>%
    ungroup()

  label_points <- bind_rows(
    fold_summary %>% group_by(mode) %>% slice_max(mae, n = 1, with_ties = FALSE) %>% mutate(label = "Worst month"),
    fold_summary %>% group_by(mode) %>% slice_min(mae, n = 1, with_ties = FALSE) %>% mutate(label = "Best month")
  )

  ggplot(fold_summary, aes(x = origin_date, y = mae)) +
    geom_hline(aes(yintercept = avg_mae), linetype = "dashed", color = wmata_colors[["dark_grey"]]) +
    geom_line(color = wmata_colors[["midnight_blue"]], linewidth = 0.9) +
    geom_point(aes(fill = season_group), color = wmata_colors[["midnight_blue"]], size = 3, shape = 21) +
    geom_text_repel(data = label_points, aes(label = label), size = 3.1, color = wmata_colors[["midnight_blue"]]) +
    scale_fill_manual(values = c("Winter fold" = wmata_colors[["dodger_blue"]], "Non-winter fold" = wmata_colors[["light_grey"]])) +
    scale_y_continuous(labels = scales::comma) +
    scale_x_date(date_labels = "%b", date_breaks = "1 month") +
    facet_wrap(~ mode, scales = "free_y", ncol = 1) +
    labs(
      title = title_text,
      subtitle = "Each point is a simulated future forecast from that month, with winter folds highlighted and the dashed line marking average MAE",
      x = "Validation origin",
      y = "MAE",
      fill = NULL
    ) +
    theme_wmata()
}

plot_horizon_decay_curve <- function(horizon_decay, title_text) {
  zone_labels <- tibble(
    xmin = c(1, 8),
    xmax = c(7, 30),
    ymin = -Inf,
    ymax = Inf,
    zone = c("Days 1-7:\nOperational planning", "Days 8-30:\nService and staffing outlook")
  )

  ggplot(horizon_decay, aes(x = horizon_day, y = smoothed_mae, color = mode)) +
    geom_rect(data = zone_labels, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), inherit.aes = FALSE, fill = wmata_colors[["light_grey"]], alpha = 0.25) +
    geom_line(linewidth = 1) +
    scale_color_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::comma) +
    scale_x_continuous(breaks = c(1, 7, 14, 21, 30)) +
    annotate("text", x = 4, y = Inf, label = "Operational planning", vjust = 1.4, color = wmata_colors[["midnight_blue"]], size = 3.2) +
    annotate("text", x = 19, y = Inf, label = "Directional staffing/service view", vjust = 1.4, color = wmata_colors[["midnight_blue"]], size = 3.2) +
    labs(
      title = title_text,
      subtitle = "The line is a 7-day rolling MAE by horizon, which makes the planning value of near-term versus farther-out forecasts easier to read",
      x = "Forecast horizon day",
      y = "MAE",
      color = NULL
    ) +
    theme_wmata()
}

plot_residual_distribution_by_mode <- function(predictions, title_text) {
  residuals <- predictions %>%
    mutate(residual = prediction - actual)

  ggplot(residuals, aes(x = residual, fill = mode)) +
    geom_histogram(bins = 30, alpha = 0.85, color = "white") +
    geom_vline(xintercept = 0, linetype = "dashed", color = wmata_colors[["midnight_blue"]]) +
    facet_wrap(~ mode, scales = "free_y") +
    scale_fill_manual(values = mode_palette) +
    scale_x_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "The two residual distributions can be compared directly for spread and bias",
      x = "Residual",
      y = "Count",
      fill = NULL
    ) +
    theme_wmata()
}

plot_predicted_vs_actual_scatter <- function(predictions, title_text) {
  ggplot(predictions, aes(x = actual, y = prediction, color = mode)) +
    geom_abline(intercept = 0, slope = 1, color = wmata_colors[["dark_grey"]], linetype = "dashed") +
    geom_point(alpha = 0.6, size = 2) +
    facet_wrap(~ mode, scales = "free") +
    scale_color_manual(values = mode_palette) +
    scale_x_continuous(labels = scales::comma) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Points near the diagonal indicate better calibration",
      x = "Actual ridership",
      y = "Predicted ridership",
      color = NULL
    ) +
    theme_wmata()
}

plot_coverage_chart <- function(coverage_summary, title_text) {
  ggplot(coverage_summary, aes(x = interval, y = observed_coverage, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    geom_point(aes(y = nominal_coverage), position = position_dodge(width = 0.7), size = 2.8, color = wmata_colors[["red"]]) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
    labs(
      title = title_text,
      subtitle = "Bars are observed holdout coverage and red points are the nominal interval targets",
      x = NULL,
      y = "Coverage",
      fill = NULL
    ) +
    theme_wmata()
}

plot_error_by_month <- function(error_by_month, title_text) {
  ggplot(error_by_month, aes(x = month_label, y = mae, group = mode, color = mode)) +
    geom_line(linewidth = 0.95) +
    geom_point(size = 2.2) +
    scale_color_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "This is a simple seasonality robustness check on the selected models",
      x = NULL,
      y = "MAE",
      color = NULL
    ) +
    theme_wmata()
}

plot_error_by_ridership_level <- function(error_by_ridership_level, title_text) {
  ggplot(error_by_ridership_level, aes(x = ridership_band, y = mae, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "This checks whether the forecast struggles disproportionately on very low or very high demand days",
      x = NULL,
      y = "MAE",
      fill = NULL
    ) +
    theme_wmata()
}

plot_calibration_curve <- function(calibration_curve, title_text) {
  ggplot(calibration_curve, aes(x = avg_prediction, y = avg_actual, color = mode)) +
    geom_abline(intercept = 0, slope = 1, color = wmata_colors[["dark_grey"]], linetype = "dashed") +
    geom_line(linewidth = 0.85) +
    geom_point(size = 2.4) +
    scale_color_manual(values = mode_palette) +
    scale_x_continuous(labels = scales::comma) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Prediction bins near the diagonal indicate better calibration",
      x = "Average predicted ridership",
      y = "Average actual ridership",
      color = NULL
    ) +
    theme_wmata()
}

plot_directional_accuracy <- function(directional_accuracy, title_text) {
  ggplot(directional_accuracy, aes(x = mode, y = directional_accuracy, fill = mode)) +
    geom_col(width = 0.6) +
    geom_text(aes(label = scales::percent(directional_accuracy, accuracy = 1)), vjust = -0.3, size = 3.6) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
    labs(
      title = title_text,
      subtitle = "Managers often care whether the model gets the direction right even before the exact magnitude",
      x = NULL,
      y = "Directional accuracy",
      fill = NULL
    ) +
    theme_wmata()
}

plot_combined_mode_forecast <- function(future_forecasts, title_text) {
  forecast_extremes <- build_forecast_extremes(future_forecasts)

  ggplot(future_forecasts, aes(x = date, y = prediction, color = mode, fill = mode)) +
    geom_ribbon(aes(ymin = lower_80, ymax = upper_80), alpha = 0.18, color = NA) +
    geom_line(linewidth = 1) +
    geom_point(data = forecast_extremes, aes(x = date, y = prediction), inherit.aes = FALSE, color = wmata_colors[["red"]], size = 2.3) +
    geom_text_repel(
      data = forecast_extremes,
      aes(x = date, y = prediction, label = paste(extreme_type, format(date, "%b %d"))),
      inherit.aes = FALSE,
      size = 3,
      color = wmata_colors[["midnight_blue"]]
    ) +
    facet_wrap(~ mode, scales = "free_y", ncol = 1) +
    scale_color_manual(values = mode_palette) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "This is the main 30-day planning view for staffing, service, and communications",
      x = NULL,
      y = "Predicted Daily Ridership",
      color = NULL,
      fill = NULL
    ) +
    theme_wmata()
}

plot_success_rate_chart <- function(success_rate_summary, title_text) {
  ggplot(success_rate_summary, aes(x = threshold_label, y = success_rate, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    geom_text(aes(label = scales::percent(success_rate, accuracy = 1)), position = position_dodge(width = 0.7), vjust = -0.25, size = 3.3) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::percent, limits = c(0, 1)) +
    labs(
      title = title_text,
      subtitle = "This is easier to interpret operationally than MAE because it answers how often the forecast is close enough for planning",
      x = NULL,
      y = "Share of holdout days",
      fill = NULL
    ) +
    theme_wmata()
}

plot_demand_uncertainty_matrix <- function(decision_matrix, title_text) {
  threshold_tbl <- decision_matrix %>%
    group_by(mode) %>%
    summarise(
      demand_midpoint = first(demand_midpoint),
      uncertainty_midpoint = first(uncertainty_midpoint),
      .groups = "drop"
    )

  label_points <- decision_matrix %>%
    slice_max(interval_width, n = 4, with_ties = FALSE)

  ggplot(decision_matrix, aes(x = prediction, y = interval_width, color = day_type)) +
    geom_vline(data = threshold_tbl, aes(xintercept = demand_midpoint), linetype = "dashed", color = wmata_colors[["dark_grey"]]) +
    geom_hline(data = threshold_tbl, aes(yintercept = uncertainty_midpoint), linetype = "dashed", color = wmata_colors[["dark_grey"]]) +
    geom_point(alpha = 0.75, size = 2.6) +
    geom_text_repel(data = label_points, aes(label = format(date, "%b %d")), size = 3.1, color = wmata_colors[["midnight_blue"]]) +
    facet_wrap(~ mode, scales = "free") +
    scale_color_manual(values = c(Weekday = wmata_colors[["midnight_blue"]], Weekend = wmata_colors[["green"]])) +
    scale_x_continuous(labels = scales::comma) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "High-demand or high-uncertainty days deserve the most operational attention; dashed lines split each mode at its own median demand and uncertainty",
      x = "Predicted ridership",
      y = "80% interval width",
      color = NULL
    ) +
    theme_wmata()
}

plot_weekday_weekend_forecast <- function(future_forecasts, title_text) {
  summary_tbl <- future_forecasts %>%
    group_by(mode, weekday_saturday_sunday) %>%
    summarise(prediction = mean(prediction, na.rm = TRUE), .groups = "drop")

  ggplot(summary_tbl, aes(x = weekday_saturday_sunday, y = prediction, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "This translates the forecast into a simpler service-frequency planning view",
      x = NULL,
      y = "Average predicted ridership",
      fill = NULL
    ) +
    theme_wmata()
}

plot_top_station_growth <- function(station_growth, title_text) {
  top_stations <- station_growth %>%
    slice_max(absolute_growth, n = 15, with_ties = FALSE)

  ggplot(top_stations, aes(x = absolute_growth, y = reorder(station_name, absolute_growth))) +
    geom_col(fill = wmata_colors[["green"]]) +
    geom_text(aes(label = scales::comma(round(absolute_growth))), hjust = -0.05, size = 3.2) +
    scale_x_continuous(labels = scales::comma, expand = expansion(mult = c(0, 0.2))) +
    labs(
      title = title_text,
      subtitle = "Growth is measured against the trailing 28-day station baseline",
      x = "Forecasted average daily growth",
      y = NULL
    ) +
    theme_wmata()
}

plot_rail_station_residual_heatmap <- function(predictions, title_text) {
  top_stations <- predictions %>%
    group_by(station_name) %>%
    summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), .groups = "drop") %>%
    slice_max(mae, n = 20, with_ties = FALSE) %>%
    pull(station_name)

  heatmap_data <- predictions %>%
    filter(station_name %in% top_stations) %>%
    mutate(
      residual = prediction - actual,
      station_name = forcats::fct_reorder(station_name, abs(residual), .fun = mean)
    )

  ggplot(heatmap_data, aes(x = date, y = station_name, fill = residual)) +
    geom_tile() +
    scale_fill_gradient2(low = wmata_colors[["red"]], mid = "white", high = wmata_colors[["dodger_blue"]], midpoint = 0, labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "Red cells are underpredictions and blue cells are overpredictions for the highest-error stations",
      x = NULL,
      y = NULL,
      fill = "Residual"
    ) +
    theme_wmata(base_size = 12)
}

plot_rail_station_risk_matrix <- function(station_risk, title_text) {
  top_labels <- station_risk %>%
    arrange(desc(forecast_avg_ridership * holdout_mae)) %>%
    slice_head(n = 12)

  ggplot(station_risk, aes(x = forecast_avg_ridership, y = holdout_mae)) +
    geom_point(color = wmata_colors[["dark_blue"]], alpha = 0.7, size = 2.4) +
    geom_text_repel(data = top_labels, aes(label = station_name), size = 3, color = wmata_colors[["midnight_blue"]]) +
    scale_x_continuous(labels = scales::comma) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "This combines forecasted volume with historical holdout error to identify stations that deserve closer monitoring",
      x = "Forecasted average daily ridership",
      y = "Historical holdout MAE"
    ) +
    theme_wmata()
}

plot_worst_prediction_days <- function(worst_days, title_text) {
  ggplot(worst_days, aes(x = abs_error, y = reorder(format(date, "%Y-%m-%d"), abs_error), color = mode)) +
    geom_segment(aes(x = 0, xend = abs_error, yend = reorder(format(date, "%Y-%m-%d"), abs_error)), linewidth = 1.1, alpha = 0.6) +
    geom_point(size = 3) +
    facet_wrap(~ mode, scales = "free_y") +
    scale_color_manual(values = mode_palette) +
    scale_x_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "The worst misses are useful for exception review and operational post-mortems",
      x = "Absolute error",
      y = NULL,
      color = NULL
    ) +
    theme_wmata()
}

plot_history_with_covid_shading <- function(history_data, title_text) {
  covid_rect <- tibble(
    xmin = as.Date("2020-03-13"),
    xmax = as.Date("2022-12-31"),
    ymin = -Inf,
    ymax = Inf
  )

  ggplot(history_data, aes(x = date, y = ridership)) +
    geom_rect(data = covid_rect, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), inherit.aes = FALSE, fill = wmata_colors[["light_grey"]], alpha = 0.6) +
    geom_line(color = wmata_colors[["sky_blue"]], linewidth = 0.35, alpha = 0.7) +
    geom_line(aes(y = rolling_28), color = wmata_colors[["midnight_blue"]], linewidth = 0.95) +
    facet_wrap(~ mode, scales = "free_y", ncol = 1) +
    scale_y_continuous(labels = scales::comma) +
    labs(
      title = title_text,
      subtitle = "The shaded period marks the pandemic and early recovery era, which is structurally different from normal operations",
      x = NULL,
      y = "Daily Ridership"
    ) +
    theme_wmata()
}

plot_recovery_normalization <- function(recovery_index, title_text) {
  ggplot(recovery_index, aes(x = month_date, y = index_to_2019, color = mode)) +
    geom_hline(yintercept = 1, linetype = "dashed", color = wmata_colors[["dark_grey"]]) +
    geom_line(linewidth = 1) +
    scale_color_manual(values = mode_palette) +
    scale_y_continuous(labels = scales::number_format(accuracy = 0.1)) +
    labs(
      title = title_text,
      subtitle = "Indexing to the 2019 baseline makes the post-COVID demand regime change visible",
      x = NULL,
      y = "Ridership index vs 2019",
      color = NULL
    ) +
    theme_wmata()
}

plot_data_freshness <- function(data_freshness, title_text) {
  ggplot(data_freshness, aes(x = mode, y = days_old, fill = mode)) +
    geom_col(width = 0.6) +
    geom_text(aes(label = paste(days_old, "days")), vjust = -0.3, size = 3.5) +
    scale_fill_manual(values = mode_palette) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
    labs(
      title = title_text,
      subtitle = "A stale input feed weakens confidence in any forecast, even when the model is strong",
      x = NULL,
      y = "Days since latest actual ridership input",
      fill = NULL
    ) +
    theme_wmata()
}

plot_accuracy_traffic_light <- function(traffic_light_data, title_text) {
  status_palette <- c(Green = wmata_colors[["green"]], Yellow = "#f3c744", Red = wmata_colors[["red"]])

  ggplot(traffic_light_data, aes(x = horizon_bucket, y = mode, fill = status)) +
    geom_tile(color = "white", linewidth = 0.8) +
    geom_text(aes(label = paste0(round(mape, 1), "%")), color = "white", fontface = "bold", size = 3.5) +
    scale_fill_manual(values = status_palette) +
    labs(
      title = title_text,
      subtitle = "The colors provide a management-friendly accuracy readout by holdout horizon bucket",
      x = NULL,
      y = NULL,
      fill = NULL
    ) +
    theme_wmata()
}

write_pipeline_tables <- function(bus_results, rail_results, silver_data) {
  create_project_directories()

  rolling_summary <- bind_rows(
    bus_results$backtest_predictions %>%
      group_by(model, origin_date) %>%
      summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), rmse = sqrt(mean((prediction - actual) ^ 2, na.rm = TRUE)), bias = mean(prediction - actual, na.rm = TRUE), .groups = "drop") %>%
      mutate(mode = "Bus", level = "system"),
    rail_results$main_results$backtest_predictions %>%
      group_by(model, origin_date) %>%
      summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), rmse = sqrt(mean((prediction - actual) ^ 2, na.rm = TRUE)), bias = mean(prediction - actual, na.rm = TRUE), .groups = "drop") %>%
      mutate(mode = "Rail", level = "station_main_cohort"),
    rail_results$system_backtest_predictions %>%
      group_by(model, origin_date) %>%
      summarise(mae = mean(abs(prediction - actual), na.rm = TRUE), rmse = sqrt(mean((prediction - actual) ^ 2, na.rm = TRUE)), bias = mean(prediction - actual, na.rm = TRUE), .groups = "drop") %>%
      mutate(mode = "Rail", level = "system_aggregated")
  )

  overall_model_comparison <- bind_rows(
    bus_results$backtest_summary %>% mutate(mode = "Bus", level = "system"),
    rail_results$main_results$backtest_summary %>% mutate(mode = "Rail", level = "station_main_cohort"),
    rail_results$system_backtest_summary %>% mutate(mode = "Rail", level = "system_aggregated")
  )

  holdout_performance <- bind_rows(
    bus_results$holdout_summary %>% mutate(mode = "Bus", level = "system"),
    rail_results$main_results$holdout_summary %>% mutate(mode = "Rail", level = "station_main_cohort"),
    rail_results$fallback_results$holdout_summary %>% mutate(mode = "Rail", level = "station_fallback"),
    rail_results$system_holdout_summary %>% mutate(mode = "Rail", level = "system_aggregated")
  )

  horizon_bucket_performance <- bind_rows(
    bus_results$horizon_summary %>% mutate(mode = "Bus", level = "system"),
    rail_results$main_results$horizon_summary %>% mutate(mode = "Rail", level = "station_main_cohort"),
    rail_results$system_horizon_summary %>% mutate(mode = "Rail", level = "system_aggregated")
  )

  bus_forecasts <- add_prediction_intervals(
    bus_results$final_future_forecasts %>% transmute(date, prediction),
    bus_results$backtest_predictions %>% filter(model == bus_results$selected_model_name)
  ) %>%
    mutate(mode = "Bus", model = bus_results$selected_model_name)

  rail_station_forecasts <- bind_rows(
    add_prediction_intervals(
      rail_results$main_results$final_future_forecasts,
      rail_results$main_results$backtest_predictions %>% filter(model == rail_results$main_results$selected_model_name)
    ) %>% mutate(forecast_group = "main_cohort", model = rail_results$main_results$selected_model_name),
    add_prediction_intervals(
      rail_results$fallback_results$future_forecasts,
      rail_results$fallback_results$backtest_predictions
    ) %>% mutate(forecast_group = "fallback", model = "fallback_hierarchy")
  )

  rail_system_forecasts <- add_prediction_intervals(
    rail_results$system_future_forecasts,
    rail_results$system_backtest_predictions %>% filter(model == rail_results$main_results$selected_model_name)
  ) %>%
    mutate(mode = "Rail", model = "station_aggregated")

  station_completeness <- rail_results$station_completeness

  unassigned_qa <- silver_data$ridership_rail_unassigned_daily %>%
    arrange(desc(ridership))

  output_map <- list(
    overall_model_comparison = overall_model_comparison,
    rolling_validation_summary = rolling_summary,
    holdout_performance = holdout_performance,
    horizon_bucket_performance = horizon_bucket_performance,
    bus_forecast_output = bus_forecasts,
    rail_station_forecast_output = rail_station_forecasts,
    rail_system_forecast_output = rail_system_forecasts,
    station_completeness_fallback = station_completeness,
    rail_unassigned_qa = unassigned_qa
  )

  purrr::imap_chr(output_map, function(data, name) {
    write_csv_portable(data, here::here(glue::glue("outputs/tables/{name}.csv")))
  })
}

write_pipeline_figures <- function(bus_results, rail_results, bus_data, silver_data) {
  create_project_directories()

  rail_system_history <- build_rail_system_history(silver_data)
  selected_backtests <- selected_backtest_predictions(bus_results, rail_results, bus_data, rail_system_history)
  selected_holdout <- selected_holdout_predictions(bus_results, rail_results, bus_data, rail_system_history)
  evaluation_predictions <- bind_rows(selected_backtests, selected_holdout)
  future_by_mode <- build_future_forecasts_by_mode(bus_results, rail_results)
  history_data <- build_ridership_history(bus_data, rail_system_history)
  recovery_index <- build_recovery_index(history_data)
  coverage_summary <- compute_coverage_summary(selected_holdout)
  horizon_decay <- build_smoothed_horizon_decay(selected_backtests)
  error_by_month <- build_error_by_month(evaluation_predictions)
  error_by_ridership_level <- build_error_by_ridership_level(evaluation_predictions)
  calibration_curve <- build_calibration_curve(evaluation_predictions)
  directional_accuracy <- build_directional_accuracy(selected_holdout)
  traffic_light_data <- build_accuracy_traffic_light(selected_holdout)
  worst_days <- build_worst_prediction_days(selected_holdout, top_n = 10)
  success_rate_summary <- build_success_rate_summary(selected_holdout)
  station_growth <- build_recent_growth_table(rail_results, silver_data)
  station_risk <- build_rail_station_risk(rail_results, silver_data)
  data_freshness <- build_data_freshness(bus_data, rail_system_history)
  decision_matrix <- build_decision_matrix(future_by_mode)

  workflow_plot <- plot_decision_workflow(
    "The forecast is designed to support short-term operating decisions, not replace planners"
  )

  bus_pattern_plot <- plot_pattern_capture(
    selected_holdout %>% filter(mode == "Bus"),
    "The model captures the weekly ridership rhythm that drives most planning decisions"
  )

  bus_holdout_plot <- plot_actual_vs_predicted(
    bus_results$holdout_predictions,
    "MetroBus holdout predictions track the systemwide daily pattern with a usable planning band",
    "January 1, 2026 to March 31, 2026"
  )

  bus_forecast_plot <- plot_forecast_trend(
    history_data = bus_data %>% select(date, ridership),
    forecast_data = future_by_mode %>% filter(mode == "Bus"),
    title_text = "MetroBus 30-day forecast provides a planning view for the next month"
  )

  bus_residual_plot <- plot_residual_distribution(
    bus_results$holdout_predictions,
    "MetroBus residuals now show both shape and practical error tolerance"
  )

  bus_horizon_plot <- plot_error_by_horizon(
    bus_results$horizon_summary,
    "MetroBus accuracy decays with horizon, but the stable candidates remain readable"
  )

  bus_service_plot <- plot_error_by_service_type(
    selected_holdout %>% filter(mode == "Bus"),
    "MetroBus errors vary by service period and should be monitored by operating context"
  )

  rail_system_holdout_plot <- plot_actual_vs_predicted(
    rail_results$system_holdout_predictions,
    "MetroRail station forecasts aggregate into a stable systemwide holdout forecast with an uncertainty band",
    "January 1, 2026 to March 31, 2026"
  )

  rail_station_rank <- plot_station_error_rank(
    rail_results$main_results$holdout_predictions,
    "A small group of stations still drives most of the remaining MetroRail holdout error"
  )

  bus_leaderboard <- plot_model_leaderboard(
    build_model_leaderboard(bus_results$backtest_predictions),
    "MetroBus model leaderboard shows the credible accuracy comparison, not the broken scale artifacts"
  )

  rail_leaderboard <- plot_model_leaderboard(
    build_model_leaderboard(rail_results$main_results$backtest_predictions),
    "MetroRail model leaderboard focuses on stable candidates in the pooled station framework"
  )

  comparison_plot <- bus_leaderboard + rail_leaderboard + patchwork::plot_layout(ncol = 2)

  importance_plot <- plot_feature_importance(
    if (nrow(rail_results$main_results$feature_importance) > 0) rail_results$main_results$feature_importance else bus_results$feature_importance,
    "The selected production model is driven mostly by recent ridership and calendar structure"
  )

  winter_actual_vs_predicted_plot <- plot_winter_actual_vs_predicted(
    evaluation_predictions,
    "Winter validation shows whether the forecast still tracks ridership in the hardest seasonal regime"
  )

  winter_mae_plot <- plot_winter_mae_comparison(
    evaluation_predictions,
    "Winter error versus non-winter error is now explicit for both Bus and Rail"
  )

  winter_weekday_pattern_plot <- plot_winter_weekday_pattern(
    evaluation_predictions,
    "Winter weekday patterns show whether the model preserves the operating rhythm of the system"
  )

  winter_coverage_plot <- plot_winter_coverage(
    selected_holdout,
    "Winter interval coverage shows whether the planning ranges stay credible during the cold season"
  )

  backtest_mae_plot <- plot_backtest_mae_over_time(
    selected_backtests,
    "Backtesting shows the model was tested across many months, not one favorable split"
  )

  horizon_decay_plot <- plot_horizon_decay_curve(
    horizon_decay,
    "Forecast value is strongest in the near term and should be treated as directional farther out"
  )

  residual_by_mode_plot <- plot_residual_distribution_by_mode(
    selected_holdout,
    "Bus and Rail residual distributions can now be compared directly for spread and bias"
  )

  scatter_plot <- plot_predicted_vs_actual_scatter(
    selected_holdout,
    "Predicted versus actual scatterplots show the holdout calibration more clearly than one summary metric"
  )

  coverage_plot <- plot_coverage_chart(
    coverage_summary,
    "Holdout interval coverage shows whether the planning ranges are too narrow or too wide"
  )

  error_by_month_plot <- plot_error_by_month(
    error_by_month,
    "Error by month shows whether the selected models stay stable across the calendar"
  )

  error_by_ridership_level_plot <- plot_error_by_ridership_level(
    error_by_ridership_level,
    "Error by ridership level tests whether high-demand and low-demand days behave differently"
  )

  calibration_plot <- plot_calibration_curve(
    calibration_curve,
    "Calibration curves show whether average predicted demand matches average actual demand"
  )

  directional_accuracy_plot <- plot_directional_accuracy(
    directional_accuracy,
    "Directional accuracy shows whether the model gets the daily up-or-down move right"
  )

  success_rate_plot <- plot_success_rate_chart(
    success_rate_summary,
    "Most forecasts are close enough for planning, but not precise enough for automatic decisions"
  )

  combined_forecast_plot <- plot_combined_mode_forecast(
    future_by_mode,
    "The 30-day forecast gives planners a mode-by-mode demand outlook"
  )

  weekday_weekend_plot <- plot_weekday_weekend_forecast(
    future_by_mode,
    "Forecasted weekday versus weekend demand helps translate ridership into service planning decisions"
  )

  top_station_growth_plot <- plot_top_station_growth(
    station_growth,
    "A small set of stations is projected to drive the next 30 days of MetroRail growth"
  )

  station_residual_heatmap_plot <- plot_rail_station_residual_heatmap(
    rail_results$main_results$holdout_predictions,
    "Station residual heatmaps quickly surface which MetroRail stations miss together and when"
  )

  station_risk_plot <- plot_rail_station_risk_matrix(
    station_risk,
    "High-ridership and high-error stations deserve closer operational monitoring"
  )

  worst_days_plot <- plot_worst_prediction_days(
    worst_days,
    "The biggest misses reveal when human review is needed"
  )

  decision_matrix_plot <- plot_demand_uncertainty_matrix(
    decision_matrix,
    "The forecast becomes actionable when WMATA prioritizes high-demand and high-uncertainty days"
  )

  covid_history_plot <- plot_history_with_covid_shading(
    history_data,
    "Historical ridership makes the COVID and recovery disruption period visually undeniable"
  )

  recovery_plot <- plot_recovery_normalization(
    recovery_index,
    "Ridership indexed to 2019 shows that current demand lives in a different post-COVID regime"
  )

  freshness_plot <- plot_data_freshness(
    data_freshness,
    "Data freshness is part of forecast trust, not just a pipeline housekeeping detail"
  )

  traffic_light_plot <- plot_accuracy_traffic_light(
    traffic_light_data,
    "Accuracy traffic lights give managers a simple Bus and Rail readiness view by holdout horizon"
  )

  figure_map <- list(
    decision_support_workflow = list(plot = workflow_plot, width = 13, height = 6),
    bus_pattern_capture = list(plot = bus_pattern_plot, width = 12, height = 7),
    top_forecast_misses = list(plot = worst_days_plot, width = 13, height = 8),
    success_rate_thresholds = list(plot = success_rate_plot, width = 11, height = 7),
    horizon_planning_value = list(plot = horizon_decay_plot, width = 12, height = 7),
    backtest_operational_stability = list(plot = backtest_mae_plot, width = 12, height = 8),
    demand_uncertainty_matrix = list(plot = decision_matrix_plot, width = 12, height = 8),
    mode_forecast_outlook = list(plot = combined_forecast_plot, width = 12, height = 8),
    bus_actual_vs_predicted = list(plot = bus_holdout_plot, width = 11, height = 7),
    bus_forecast_trend = list(plot = bus_forecast_plot, width = 11, height = 7),
    bus_residual_distribution = list(plot = bus_residual_plot, width = 11, height = 7),
    bus_error_by_horizon = list(plot = bus_horizon_plot, width = 11, height = 7),
    bus_error_by_service_type = list(plot = bus_service_plot, width = 11, height = 7),
    rail_actual_vs_predicted = list(plot = rail_system_holdout_plot, width = 11, height = 7),
    rail_station_error_ranking = list(plot = rail_station_rank, width = 11, height = 8),
    model_comparison = list(plot = comparison_plot, width = 16, height = 8),
    feature_importance = list(plot = importance_plot, width = 11, height = 7),
    winter_actual_vs_predicted = list(plot = winter_actual_vs_predicted_plot, width = 15, height = 8),
    winter_mae_vs_nonwinter = list(plot = winter_mae_plot, width = 11, height = 7),
    winter_weekday_pattern = list(plot = winter_weekday_pattern_plot, width = 12, height = 7),
    winter_interval_coverage = list(plot = winter_coverage_plot, width = 11, height = 7),
    backtest_mae_over_time = list(plot = backtest_mae_plot, width = 12, height = 7),
    horizon_decay_curve = list(plot = horizon_decay_plot, width = 12, height = 7),
    residual_distribution_by_mode = list(plot = residual_by_mode_plot, width = 12, height = 7),
    predicted_vs_actual_scatter = list(plot = scatter_plot, width = 12, height = 7),
    coverage_chart = list(plot = coverage_plot, width = 11, height = 7),
    error_by_month = list(plot = error_by_month_plot, width = 12, height = 7),
    error_by_ridership_level = list(plot = error_by_ridership_level_plot, width = 12, height = 7),
    calibration_curve = list(plot = calibration_plot, width = 12, height = 7),
    directional_accuracy = list(plot = directional_accuracy_plot, width = 11, height = 7),
    demand_forecast_by_mode = list(plot = combined_forecast_plot, width = 12, height = 8),
    forecast_weekday_weekend_mix = list(plot = weekday_weekend_plot, width = 11, height = 7),
    rail_top_stations_by_growth = list(plot = top_station_growth_plot, width = 12, height = 8),
    rail_station_residual_heatmap = list(plot = station_residual_heatmap_plot, width = 14, height = 8),
    rail_station_risk_matrix = list(plot = station_risk_plot, width = 12, height = 8),
    worst_forecast_days = list(plot = worst_days_plot, width = 12, height = 7),
    ridership_history_with_covid_shading = list(plot = covid_history_plot, width = 12, height = 8),
    recovery_normalization = list(plot = recovery_plot, width = 12, height = 7),
    data_freshness_monitor = list(plot = freshness_plot, width = 10, height = 6),
    accuracy_traffic_light = list(plot = traffic_light_plot, width = 11, height = 6)
  )

  purrr::imap_chr(figure_map, function(spec, name) {
    save_plot_wmata(
      spec$plot,
      here::here(glue::glue("outputs/figures/{name}.png")),
      width = spec$width,
      height = spec$height
    )
  })
}

write_diagnostic_artifacts <- function(bus_results, rail_results) {
  create_project_directories()

  artifacts <- list(bus_results = bus_results, rail_results = rail_results)
  safe_qsave(artifacts, here::here("outputs/diagnostics/wmata_model_artifacts.qs"))
}