discovery_layer.R

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

build_discovery_mode_history <- function(silver_data) {
  rail_system_daily <- 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"
    ) %>%
    mutate(mode = "Rail")

  bind_rows(
    silver_data$ridership_bus_daily %>%
      select(date, ridership, holiday, day_of_week, weekday_saturday_sunday, service_type) %>%
      mutate(mode = "Bus"),
    rail_system_daily
  ) %>%
    arrange(mode, date) %>%
    mutate(
      month_num = lubridate::month(date),
      month_label = factor(month.abb[month_num], levels = month.abb),
      year = lubridate::year(date),
      day_of_week = factor(day_of_week, levels = c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")),
      weekday_saturday_sunday = factor(weekday_saturday_sunday, levels = c("Weekday", "Saturday", "Sunday")),
      winter_flag = factor(if_else(month_num %in% c(12L, 1L, 2L), "Winter", "Non-winter"), levels = c("Non-winter", "Winter")),
      holiday = if_else(holiday == "Yes", "Yes", "No"),
      covid_period = case_when(
        date < as.Date("2020-03-13") ~ "Pre COVID",
        date < as.Date("2023-01-01") ~ "COVID/Recovery",
        TRUE ~ "Post COVID"
      )
    )
}

add_discovery_baseline <- function(mode_history) {
  baseline_reference <- mode_history %>%
    filter(covid_period == "Post COVID") %>%
    group_by(mode, month_num, day_of_week, holiday) %>%
    summarise(
      expected_ridership = median(ridership, na.rm = TRUE),
      baseline_n = n(),
      .groups = "drop"
    )

  mode_history %>%
    left_join(baseline_reference, by = c("mode", "month_num", "day_of_week", "holiday")) %>%
    mutate(
      ridership_lift = ridership / pmax(expected_ridership, 1) - 1
    )
}

build_discovery_selected_predictions <- function(bus_results, rail_results, silver_data) {
  rail_history <- silver_data$ridership_rail_station_daily %>%
    group_by(date) %>%
    summarise(
      actual = 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"
    )

  bus_history <- silver_data$ridership_bus_daily %>%
    transmute(
      date,
      actual = ridership,
      holiday,
      day_of_week,
      weekday_saturday_sunday,
      service_type
    )

  bind_rows(
    bus_results$backtest_predictions %>%
      filter(model == bus_results$selected_model_name) %>%
      left_join(bus_history, by = c("date", "actual")) %>%
      mutate(mode = "Bus", evaluation_split = "backtest"),
    bus_results$holdout_predictions %>%
      left_join(bus_history, by = c("date", "actual")) %>%
      mutate(mode = "Bus", evaluation_split = "holdout"),
    rail_results$system_backtest_predictions %>%
      filter(model == rail_results$main_results$selected_model_name) %>%
      left_join(rail_history, by = c("date", "actual")) %>%
      mutate(mode = "Rail", evaluation_split = "backtest"),
    rail_results$system_holdout_predictions %>%
      left_join(rail_history, by = c("date", "actual")) %>%
      mutate(mode = "Rail", evaluation_split = "holdout")
  ) %>%
    mutate(
      month_num = lubridate::month(date),
      month_label = factor(month.abb[month_num], levels = month.abb),
      year = lubridate::year(date),
      winter_flag = factor(if_else(month_num %in% c(12L, 1L, 2L), "Winter", "Non-winter"), levels = c("Non-winter", "Winter")),
      residual = actual - prediction,
      abs_error = abs(residual),
      interval_width_80 = upper_80 - lower_80,
      outside_80 = actual < lower_80 | actual > upper_80
    )
}

build_discovery_future_signals <- function(bus_results, rail_results, mode_history) {
  bus_future <- add_prediction_intervals(
    bus_results$final_future_forecasts %>% select(date, prediction),
    bus_results$backtest_predictions %>% filter(model == bus_results$selected_model_name)
  ) %>%
    mutate(mode = "Bus")

  rail_future <- 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")

  history_thresholds <- mode_history %>%
    filter(covid_period == "Post COVID") %>%
    transmute(
      mode,
      weekday_saturday_sunday,
      ridership
    ) %>%
    group_by(mode, weekday_saturday_sunday) %>%
    summarise(
      demand_p10 = quantile(ridership, probs = 0.10, na.rm = TRUE),
      demand_p90 = quantile(ridership, probs = 0.90, na.rm = TRUE),
      .groups = "drop"
    )

  future_signals <- bind_rows(bus_future, rail_future) %>%
    mutate(
      day_of_week = factor(wday(date, label = TRUE, abbr = TRUE, week_start = 1), levels = c("Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun")),
      weekday_saturday_sunday = factor(
        case_when(
          day_of_week %in% c("Mon", "Tue", "Wed", "Thu", "Fri") ~ "Weekday",
          day_of_week == "Sat" ~ "Saturday",
          TRUE ~ "Sunday"
        ),
        levels = c("Weekday", "Saturday", "Sunday")
      ),
      interval_width_80 = upper_80 - lower_80
    ) %>%
    left_join(history_thresholds, by = c("mode", "weekday_saturday_sunday"))

  uncertainty_thresholds <- future_signals %>%
    group_by(mode) %>%
    summarise(
      uncertainty_p90 = quantile(interval_width_80, probs = 0.90, na.rm = TRUE),
      .groups = "drop"
    )

  future_signals %>%
    left_join(uncertainty_thresholds, by = "mode") %>%
    mutate(
      high_demand_flag = prediction >= demand_p90,
      low_demand_flag = prediction <= demand_p10,
      high_uncertainty_flag = interval_width_80 >= uncertainty_p90
    )
}

analyse_weekday_structure <- function(mode_history) {
  weekday_index <- mode_history %>%
    filter(covid_period %in% c("Pre COVID", "Post COVID")) %>%
    group_by(mode, covid_period, day_of_week) %>%
    summarise(avg_ridership = mean(ridership, na.rm = TRUE), .groups = "drop_last") %>%
    mutate(weekday_index = avg_ridership / mean(avg_ridership, na.rm = TRUE)) %>%
    ungroup()

  weekday_tests <- mode_history %>%
    filter(covid_period == "Post COVID", day_of_week %in% c("Mon", "Tue", "Wed", "Thu", "Fri")) %>%
    mutate(tue_thu_flag = day_of_week %in% c("Tue", "Wed", "Thu")) %>%
    group_split(mode) %>%
    map_dfr(function(df) {
      weekday_model <- lm(log1p(ridership) ~ tue_thu_flag + factor(month_num) + holiday, data = df)
      tue_thu_row <- broom::tidy(weekday_model) %>% filter(term == "tue_thu_flagTRUE")
      weekday_ttest <- t.test(ridership ~ tue_thu_flag, data = df)

      tibble(
        mode = unique(df$mode),
        mon_fri_avg_ridership = mean(df$ridership[!df$tue_thu_flag], na.rm = TRUE),
        tue_thu_avg_ridership = mean(df$ridership[df$tue_thu_flag], na.rm = TRUE),
        raw_pct_gap = tue_thu_avg_ridership / pmax(mon_fri_avg_ridership, 1) - 1,
        adjusted_pct_gap = exp(tue_thu_row$estimate) - 1,
        model_p_value = tue_thu_row$p.value,
        t_test_p_value = weekday_ttest$p.value
      )
    })

  list(
    weekday_index = weekday_index,
    weekday_tests = weekday_tests
  )
}

analyse_winter_penalty <- function(mode_history, selected_predictions) {
  winter_ridership <- mode_history %>%
    filter(covid_period == "Post COVID") %>%
    group_split(mode) %>%
    map_dfr(function(df) {
      winter_model <- lm(log1p(ridership) ~ winter_flag + day_of_week + factor(year) + holiday, data = df)
      winter_row <- broom::tidy(winter_model) %>% filter(term == "winter_flagWinter")

      winter_summary <- df %>%
        group_by(winter_flag) %>%
        summarise(avg_ridership = mean(ridership, na.rm = TRUE), .groups = "drop") %>%
        tidyr::pivot_wider(names_from = winter_flag, values_from = avg_ridership)

      tibble(
        mode = unique(df$mode),
        non_winter_avg_ridership = winter_summary$`Non-winter`,
        winter_avg_ridership = winter_summary$Winter,
        adjusted_winter_effect = exp(winter_row$estimate) - 1,
        ridership_p_value = winter_row$p.value
      )
    })

  winter_error <- selected_predictions %>%
    group_split(mode) %>%
    map_dfr(function(df) {
      error_ttest <- t.test(abs_error ~ winter_flag, data = df)

      error_summary <- df %>%
        group_by(winter_flag) %>%
        summarise(
          mae = mean(abs_error, na.rm = TRUE),
          .groups = "drop"
        ) %>%
        tidyr::pivot_wider(names_from = winter_flag, values_from = mae)

      tibble(
        mode = unique(df$mode),
        non_winter_mae = error_summary$`Non-winter`,
        winter_mae = error_summary$Winter,
        winter_error_premium = winter_mae / pmax(non_winter_mae, 1) - 1,
        error_p_value = error_ttest$p.value
      )
    })

  list(
    ridership = winter_ridership,
    error = winter_error
  )
}

analyse_mode_divergence <- function(mode_history) {
  lifted_history <- add_discovery_baseline(mode_history) %>%
    filter(covid_period == "Post COVID")

  bus_lift <- lifted_history %>%
    filter(mode == "Bus") %>%
    transmute(
      date,
      day_of_week,
      weekday_saturday_sunday,
      holiday,
      month_label,
      bus_lift = ridership_lift
    )

  rail_lift <- lifted_history %>%
    filter(mode == "Rail") %>%
    transmute(date, rail_lift = ridership_lift)

  divergence_days <- bus_lift %>%
    inner_join(rail_lift, by = "date") %>%
    mutate(
      divergence_gap = bus_lift - rail_lift,
      divergence_type = case_when(
        bus_lift >= 0.15 & rail_lift <= -0.15 ~ "Bus above / Rail below",
        rail_lift >= 0.15 & bus_lift <= -0.15 ~ "Rail above / Bus below",
        TRUE ~ NA_character_
      )
    ) %>%
    arrange(desc(abs(divergence_gap)))

  divergence_summary <- divergence_days %>%
    filter(!is.na(divergence_type)) %>%
    group_by(divergence_type, weekday_saturday_sunday) %>%
    summarise(
      days_flagged = n(),
      avg_gap = mean(abs(divergence_gap), na.rm = TRUE),
      .groups = "drop"
    )

  list(
    divergence_days = divergence_days,
    divergence_summary = divergence_summary
  )
}

analyse_error_concentration <- function(selected_predictions, rail_results) {
  holdout_predictions <- selected_predictions %>%
    filter(evaluation_split == "holdout") %>%
    arrange(mode, desc(abs_error))

  system_curve <- holdout_predictions %>%
    group_by(mode) %>%
    mutate(
      day_rank = row_number(),
      day_share = day_rank / n(),
      cumulative_error_share = cumsum(abs_error) / sum(abs_error, na.rm = TRUE),
      top_decile_flag = day_rank <= ceiling(n() * 0.10)
    ) %>%
    ungroup()

  system_summary <- system_curve %>%
    group_by(mode) %>%
    summarise(
      top_decile_error_share = sum(abs_error[top_decile_flag], na.rm = TRUE) / sum(abs_error, na.rm = TRUE),
      all_day_mae = mean(abs_error, na.rm = TRUE),
      normal_day_mae = mean(abs_error[!top_decile_flag], na.rm = TRUE),
      .groups = "drop"
    )

  top_error_days <- holdout_predictions %>%
    mutate(
      likely_explanation = case_when(
        holiday == "Yes" ~ "Holiday or special calendar effect",
        weekday_saturday_sunday != "Weekday" ~ "Weekend-like demand shift",
        TRUE ~ "Possible disruption or unmodeled event"
      )
    ) %>%
    group_by(mode) %>%
    slice_max(abs_error, n = 10, with_ties = FALSE) %>%
    ungroup()

  rail_station_holdout <- bind_rows(
    rail_results$main_results$holdout_predictions,
    rail_results$fallback_results$holdout_predictions
  ) %>%
    mutate(abs_error = abs(actual - prediction))

  station_pareto <- rail_station_holdout %>%
    group_by(station_name) %>%
    summarise(
      total_abs_error = sum(abs_error, na.rm = TRUE),
      average_abs_error = mean(abs_error, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    arrange(desc(total_abs_error)) %>%
    mutate(
      station_rank = row_number(),
      error_share = total_abs_error / sum(total_abs_error, na.rm = TRUE),
      cumulative_error_share = cumsum(error_share)
    )

  list(
    system_curve = system_curve,
    system_summary = system_summary,
    top_error_days = top_error_days,
    station_pareto = station_pareto
  )
}

analyse_station_sensitivity <- function(rail_station_model_frame) {
  station_effects <- rail_station_model_frame %>%
    filter(date >= as.Date("2023-01-01")) %>%
    mutate(
      weekend_flag = weekday_saturday_sunday != "Weekday",
      tue_thu_flag = day_of_week %in% c("Tue", "Wed", "Thu"),
      holiday_yes = holiday == "Yes"
    ) %>%
    group_split(station_name) %>%
    map_dfr(function(df) {
      if (nrow(df) < 365) {
        return(tibble())
      }

      station_model <- lm(log1p(ridership) ~ weekend_flag + tue_thu_flag + holiday_yes + factor(month), data = df)

      broom::tidy(station_model) %>%
        filter(term %in% c("weekend_flagTRUE", "tue_thu_flagTRUE", "holiday_yesTRUE")) %>%
        transmute(
          station_name = unique(df$station_name),
          effect_type = recode(
            term,
            weekend_flagTRUE = "Weekend sensitivity",
            tue_thu_flagTRUE = "Tue-Thu commuter premium",
            holiday_yesTRUE = "Holiday sensitivity"
          ),
          effect_pct = 100 * (exp(estimate) - 1),
          p_value = p.value
        )
    })

  top_sensitivity <- station_effects %>%
    group_by(effect_type) %>%
    slice_max(abs(effect_pct), n = 15, with_ties = FALSE) %>%
    ungroup()

  list(
    station_effects = station_effects,
    top_sensitivity = top_sensitivity
  )
}

analyse_forecast_triggers <- function(future_signals, selected_predictions) {
  trigger_days <- future_signals %>%
    mutate(
      trigger_count = rowSums(cbind(high_demand_flag, low_demand_flag, high_uncertainty_flag)),
      trigger_summary = str_squish(paste(
        if_else(high_demand_flag, "High demand", ""),
        if_else(low_demand_flag, "Low demand", ""),
        if_else(high_uncertainty_flag, "High uncertainty", ""),
        sep = " | "
      )),
      trigger_summary = str_replace_all(trigger_summary, "\\|\\s+\\|", "|"),
      trigger_summary = str_replace_all(trigger_summary, "^\\s*\\|\\s*|\\s*\\|\\s*$", ""),
      recommended_action = case_when(
        high_demand_flag & high_uncertainty_flag ~ "Prepare capacity and contingency",
        high_demand_flag ~ "Staff and capacity planning",
        low_demand_flag ~ "Maintenance or resource optimization",
        high_uncertainty_flag ~ "Manual review before acting",
        TRUE ~ "Routine monitoring"
      )
    ) %>%
    filter(trigger_count > 0) %>%
    arrange(mode, date)

  trigger_counts <- trigger_days %>%
    transmute(
      mode,
      date,
      high_demand_flag,
      low_demand_flag,
      high_uncertainty_flag
    ) %>%
    pivot_longer(cols = ends_with("_flag"), names_to = "trigger_type", values_to = "flagged") %>%
    filter(flagged) %>%
    mutate(
      trigger_type = recode(
        trigger_type,
        high_demand_flag = "High-demand day",
        low_demand_flag = "Low-demand day",
        high_uncertainty_flag = "High-uncertainty day"
      )
    ) %>%
    count(mode, trigger_type, name = "days_flagged")

  monitoring_examples <- selected_predictions %>%
    filter(evaluation_split == "holdout") %>%
    arrange(mode, date) %>%
    group_by(mode) %>%
    mutate(
      outside_change = outside_80 != dplyr::lag(outside_80, default = FALSE),
      run_id = cumsum(outside_change),
      rolling_bias_7 = slider::slide_dbl(residual, mean, .before = 6, .complete = TRUE, na.rm = TRUE)
    ) %>%
    group_by(mode, run_id) %>%
    summarise(
      trigger_name = if_else(all(outside_80) & n() >= 3, "Pattern break", "No trigger"),
      start_date = min(date),
      end_date = max(date),
      days_in_run = n(),
      .groups = "drop"
    ) %>%
    filter(trigger_name == "Pattern break")

  list(
    trigger_days = trigger_days,
    trigger_counts = trigger_counts,
    monitoring_examples = monitoring_examples
  )
}

build_discovery_insight_cards <- function(discovery_results) {
  weekday_card <- discovery_results$weekday_structure$weekday_tests %>%
    arrange(desc(adjusted_pct_gap)) %>%
    slice(1)

  winter_card <- discovery_results$winter_penalty$error %>%
    arrange(desc(winter_error_premium)) %>%
    slice(1)

  error_card <- discovery_results$error_concentration$system_summary %>%
    arrange(desc(top_decile_error_share)) %>%
    slice(1)

  station_card <- discovery_results$error_concentration$station_pareto %>%
    slice_head(n = 10) %>%
    summarise(top_station_share = sum(error_share, na.rm = TRUE)) %>%
    mutate(mode = "Rail")

  divergence_card <- discovery_results$mode_divergence$divergence_summary %>%
    arrange(desc(days_flagged), desc(avg_gap)) %>%
    slice(1)

  trigger_card <- discovery_results$forecast_triggers$trigger_counts %>%
    arrange(desc(days_flagged)) %>%
    slice(1)

  cards <- list(
    list(
      finding = paste0(weekday_card$mode, " ridership is materially stronger Tuesday to Thursday than Monday and Friday."),
      evidence = paste0(
        "Post-COVID Tue-Thu ridership is ",
        percent(weekday_card$adjusted_pct_gap, accuracy = 0.1),
        " above Mon/Fri after controlling for month and holidays (p = ",
        signif(weekday_card$model_p_value, 3),
        ")."
      ),
      operational_meaning = "Weekday service assumptions should not treat all weekdays the same.",
      action = "Use separate weekday demand profiles for Mon/Fri versus Tue-Thu.",
      graph = "outputs/figures/discovery_weekday_index.png"
    ),
    list(
      finding = paste0(winter_card$mode, " forecast error is materially worse in winter than outside winter."),
      evidence = paste0(
        "Winter MAE is ",
        percent(winter_card$winter_error_premium, accuracy = 0.1),
        " higher than non-winter MAE (p = ",
        signif(winter_card$error_p_value, 3),
        ")."
      ),
      operational_meaning = "Cold-season planning should carry more caution and review.",
      action = "Add a winter review gate for forecast interpretation and staffing assumptions.",
      graph = "outputs/figures/discovery_winter_penalty.png"
    ),
    list(
      finding = paste0("A small number of ", error_card$mode, " days explain a disproportionate share of total forecast error."),
      evidence = paste0(
        "The top 10% highest-error days account for ",
        percent(error_card$top_decile_error_share, accuracy = 0.1),
        " of holdout absolute error."
      ),
      operational_meaning = "Most days are forecastable; the biggest misses are concentrated in abnormal periods.",
      action = "Prioritize better event, disruption, and service-context inputs before changing the core model.",
      graph = "outputs/figures/discovery_error_concentration.png"
    ),
    list(
      finding = "MetroRail forecast error is concentrated in a relatively small set of stations.",
      evidence = paste0(
        "The top 10 stations account for ",
        percent(station_card$top_station_share, accuracy = 0.1),
        " of station-level absolute error."
      ),
      operational_meaning = "Forecast improvement is more valuable at a limited number of high-impact stations.",
      action = "Target monitoring and local diagnostics at the highest-contribution stations first.",
      graph = "outputs/figures/discovery_station_error_pareto.png"
    ),
    list(
      finding = paste0(divergence_card$divergence_type, " is the most common Bus/Rail divergence pattern."),
      evidence = paste0(
        divergence_card$days_flagged,
        " post-COVID days were flagged, with an average lift gap of ",
        percent(divergence_card$avg_gap, accuracy = 0.1),
        "."
      ),
      operational_meaning = "Mode divergence likely reflects different trip purposes rather than broad regional demand shifts.",
      action = "Use divergence days as prompts to check events, tourism, airport travel, or local access factors.",
      graph = "outputs/figures/discovery_mode_divergence.png"
    ),
    list(
      finding = paste0(trigger_card$mode, " has the largest upcoming concentration of ", tolower(trigger_card$trigger_type), "."),
      evidence = paste0(
        trigger_card$days_flagged,
        " forecast days in the next 30 days are flagged for that trigger."
      ),
      operational_meaning = "The discovery layer can be turned into an operational review queue instead of just a model artifact.",
      action = "Surface trigger days in a recurring planning dashboard with recommended actions.",
      graph = "outputs/figures/discovery_forecast_triggers.png"
    )
  )

  cards
}

build_discovery_results <- function(bus_results, rail_results, gold_data, silver_data) {
  mode_history <- build_discovery_mode_history(silver_data)
  selected_predictions <- build_discovery_selected_predictions(bus_results, rail_results, silver_data)
  future_signals <- build_discovery_future_signals(bus_results, rail_results, mode_history)

  weekday_structure <- analyse_weekday_structure(mode_history)
  winter_penalty <- analyse_winter_penalty(mode_history, selected_predictions)
  mode_divergence <- analyse_mode_divergence(mode_history)
  error_concentration <- analyse_error_concentration(selected_predictions, rail_results)
  station_sensitivity <- analyse_station_sensitivity(gold_data$rail_station_model_frame)
  forecast_triggers <- analyse_forecast_triggers(future_signals, selected_predictions)

  discovery_results <- list(
    mode_history = mode_history,
    selected_predictions = selected_predictions,
    future_signals = future_signals,
    weekday_structure = weekday_structure,
    winter_penalty = winter_penalty,
    mode_divergence = mode_divergence,
    error_concentration = error_concentration,
    station_sensitivity = station_sensitivity,
    forecast_triggers = forecast_triggers
  )

  discovery_results$insight_cards <- build_discovery_insight_cards(discovery_results)
  discovery_results
}

plot_discovery_weekday_index <- function(weekday_index) {
  ggplot(weekday_index, aes(x = day_of_week, y = weekday_index, color = covid_period, group = covid_period)) +
    geom_hline(yintercept = 1, color = wmata_colors[["dark_grey"]], linetype = "dashed") +
    geom_line(linewidth = 1) +
    geom_point(size = 2.5) +
    facet_wrap(~ mode, scales = "free_y") +
    scale_color_manual(values = c("Pre COVID" = wmata_colors[["dark_blue"]], "Post COVID" = wmata_colors[["dodger_blue"]])) +
    scale_y_continuous(labels = percent_format(accuracy = 1)) +
    labs(
      title = "Weekday demand structure has shifted since the pre-COVID commute pattern",
      subtitle = "Index values above 1 indicate days stronger than each period's own weekly average",
      x = NULL,
      y = "Weekday ridership index",
      color = NULL
    ) +
    theme_wmata()
}

plot_discovery_winter_penalty <- function(winter_error) {
  winter_long <- winter_error %>%
    select(mode, non_winter_mae, winter_mae) %>%
    pivot_longer(cols = c(non_winter_mae, winter_mae), names_to = "season_group", values_to = "mae") %>%
    mutate(
      season_group = recode(season_group, non_winter_mae = "Non-winter", winter_mae = "Winter")
    )

  ggplot(winter_long, aes(x = season_group, y = mae, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    geom_text(aes(label = comma(round(mae))), position = position_dodge(width = 0.7), vjust = -0.25, size = 3.3) +
    scale_fill_manual(values = c(Bus = wmata_colors[["dodger_blue"]], Rail = wmata_colors[["green"]])) +
    scale_y_continuous(labels = comma, expand = expansion(mult = c(0, 0.12))) +
    labs(
      title = "Winter raises forecast error for both modes, but the premium is not equal",
      subtitle = "This is a planning-confidence view, not just a model metric",
      x = NULL,
      y = "MAE",
      fill = NULL
    ) +
    theme_wmata()
}

plot_discovery_mode_divergence <- function(divergence_days) {
  plot_data <- divergence_days %>%
    mutate(
      divergence_flag = !is.na(divergence_type)
    )

  labels_tbl <- plot_data %>%
    filter(divergence_flag) %>%
    slice_max(abs(divergence_gap), n = 10, with_ties = FALSE)

  ggplot(plot_data, aes(x = bus_lift, y = rail_lift)) +
    geom_hline(yintercept = c(-0.15, 0, 0.15), color = wmata_colors[["light_grey"]], linetype = c("dashed", "solid", "dashed")) +
    geom_vline(xintercept = c(-0.15, 0, 0.15), color = wmata_colors[["light_grey"]], linetype = c("dashed", "solid", "dashed")) +
    geom_point(aes(color = divergence_flag), alpha = 0.7, size = 2.2) +
    geom_text_repel(data = labels_tbl, aes(label = format(date, "%Y-%m-%d")), size = 3, color = wmata_colors[["midnight_blue"]]) +
    scale_color_manual(values = c(`TRUE` = wmata_colors[["red"]], `FALSE` = wmata_colors[["dodger_blue"]]), labels = c(`TRUE` = "Flagged divergence", `FALSE` = "Typical day")) +
    scale_x_continuous(labels = percent_format(accuracy = 1)) +
    scale_y_continuous(labels = percent_format(accuracy = 1)) +
    labs(
      title = "Bus and Rail usually move together, but a subset of days shows clear mode divergence",
      subtitle = "Points in opposite quadrants indicate days where one mode outperformed its baseline while the other underperformed",
      x = "Bus lift versus expected same day type",
      y = "Rail lift versus expected same day type",
      color = NULL
    ) +
    theme_wmata()
}

plot_discovery_error_concentration <- function(system_curve) {
  ggplot(system_curve, aes(x = day_share, y = cumulative_error_share, color = mode)) +
    geom_line(linewidth = 1) +
    geom_abline(intercept = 0, slope = 1, color = wmata_colors[["dark_grey"]], linetype = "dashed") +
    scale_color_manual(values = c(Bus = wmata_colors[["dodger_blue"]], Rail = wmata_colors[["green"]])) +
    scale_x_continuous(labels = percent_format(accuracy = 1)) +
    scale_y_continuous(labels = percent_format(accuracy = 1)) +
    labs(
      title = "A small share of days explains a disproportionate share of total forecast error",
      subtitle = "The farther the line bends above the diagonal, the more concentrated the error is in abnormal days",
      x = "Share of holdout days",
      y = "Cumulative share of absolute error",
      color = NULL
    ) +
    theme_wmata()
}

plot_discovery_station_error_pareto <- function(station_pareto) {
  top_stations <- station_pareto %>% slice_head(n = 20)

  ggplot(top_stations, aes(x = reorder(station_name, total_abs_error), y = total_abs_error)) +
    geom_col(fill = wmata_colors[["dark_blue"]]) +
    coord_flip() +
    scale_y_continuous(labels = comma) +
    labs(
      title = "Rail forecast error is concentrated in a limited number of stations",
      subtitle = "These are the highest-value stations for targeted monitoring and future feature improvement",
      x = NULL,
      y = "Total holdout absolute error"
    ) +
    theme_wmata()
}

plot_discovery_station_sensitivity <- function(top_sensitivity) {
  ggplot(top_sensitivity, aes(x = effect_pct, y = reorder(station_name, effect_pct), fill = effect_type)) +
    geom_col(show.legend = FALSE) +
    facet_wrap(~ effect_type, scales = "free_y") +
    scale_fill_manual(values = c(
      "Weekend sensitivity" = wmata_colors[["green"]],
      "Tue-Thu commuter premium" = wmata_colors[["dodger_blue"]],
      "Holiday sensitivity" = wmata_colors[["red"]]
    )) +
    scale_x_continuous(labels = function(x) paste0(round(x), "%")) +
    labs(
      title = "Stations do not respond equally to calendar conditions",
      subtitle = "The largest sensitivity effects identify where generic system assumptions are least reliable",
      x = "Estimated ridership effect",
      y = NULL
    ) +
    theme_wmata(base_size = 12)
}

plot_discovery_forecast_triggers <- function(trigger_counts) {
  ggplot(trigger_counts, aes(x = trigger_type, y = days_flagged, fill = mode)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.65) +
    geom_text(aes(label = days_flagged), position = position_dodge(width = 0.7), vjust = -0.25, size = 3.3) +
    scale_fill_manual(values = c(Bus = wmata_colors[["dodger_blue"]], Rail = wmata_colors[["green"]])) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
    labs(
      title = "The next 30 days can be turned into a trigger queue rather than a passive forecast",
      subtitle = "Triggered days should drive manual review, staffing checks, or resource optimization decisions",
      x = NULL,
      y = "Forecast days flagged",
      fill = NULL
    ) +
    theme_wmata() +
    theme(axis.text.x = element_text(angle = 12, hjust = 1))
}

write_discovery_tables <- function(discovery_results) {
  create_project_directories()

  output_map <- list(
    discovery_weekday_index = discovery_results$weekday_structure$weekday_index,
    discovery_weekday_tests = discovery_results$weekday_structure$weekday_tests,
    discovery_winter_ridership = discovery_results$winter_penalty$ridership,
    discovery_winter_error = discovery_results$winter_penalty$error,
    discovery_mode_divergence_days = discovery_results$mode_divergence$divergence_days,
    discovery_mode_divergence_summary = discovery_results$mode_divergence$divergence_summary,
    discovery_error_concentration = discovery_results$error_concentration$system_summary,
    discovery_top_error_days = discovery_results$error_concentration$top_error_days,
    discovery_station_error_pareto = discovery_results$error_concentration$station_pareto,
    discovery_station_sensitivity = discovery_results$station_sensitivity$station_effects,
    discovery_station_sensitivity_top = discovery_results$station_sensitivity$top_sensitivity,
    discovery_trigger_days = discovery_results$forecast_triggers$trigger_days,
    discovery_trigger_counts = discovery_results$forecast_triggers$trigger_counts,
    discovery_monitoring_examples = discovery_results$forecast_triggers$monitoring_examples
  )

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

  insight_lines <- c("# WMATA Discovery Insight Cards", "")

  for (i in seq_along(discovery_results$insight_cards)) {
    card <- discovery_results$insight_cards[[i]]
    insight_lines <- c(
      insight_lines,
      paste0("## Insight ", i),
      "",
      paste0("Finding: ", card$finding),
      paste0("Evidence: ", card$evidence),
      paste0("Operational meaning: ", card$operational_meaning),
      paste0("Action: ", card$action),
      paste0("Graph: ", card$graph),
      ""
    )
  }

  insight_path <- here::here("outputs/tables/discovery_insight_cards.md")
  writeLines(insight_lines, con = insight_path, useBytes = TRUE)

  c(csv_files, insight_path)
}

write_discovery_figures <- function(discovery_results) {
  create_project_directories()

  figure_map <- list(
    discovery_weekday_index = plot_discovery_weekday_index(discovery_results$weekday_structure$weekday_index),
    discovery_winter_penalty = plot_discovery_winter_penalty(discovery_results$winter_penalty$error),
    discovery_mode_divergence = plot_discovery_mode_divergence(discovery_results$mode_divergence$divergence_days),
    discovery_error_concentration = plot_discovery_error_concentration(discovery_results$error_concentration$system_curve),
    discovery_station_error_pareto = plot_discovery_station_error_pareto(discovery_results$error_concentration$station_pareto),
    discovery_station_sensitivity = plot_discovery_station_sensitivity(discovery_results$station_sensitivity$top_sensitivity),
    discovery_forecast_triggers = plot_discovery_forecast_triggers(discovery_results$forecast_triggers$trigger_counts)
  )

  purrr::imap_chr(figure_map, function(plot_obj, name) {
    save_plot_wmata(plot_obj, here::here(glue::glue("outputs/figures/{name}.png")))
  })
}

write_discovery_artifact <- function(discovery_results) {
  create_project_directories()
  safe_qsave(discovery_results, here::here("outputs/diagnostics/discovery_artifacts.qs"))
}