Gradient-Boosting-Modelle über Trading-Fallstudien vergleichen

# 12_case_study_insights.py


```py
# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: tags,-all
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.19.3
#   kernelspec:
#     display_name: Python 3 (ipykernel)
#     language: python
#     name: python3
# ---

# %% [markdown]
# # Case Study Insights: Gradient Boosting
#
# **Docker image**: `ml4t`
#
# **Purpose**: assemble the cross-case-study view of LightGBM and contrast it
# with the linear baseline (Ch11) and the tabular deep-learning variant (TabM).
# Per-case-study deep dives live in `case_studies/{cs}/13_model_analysis.py`;
# this notebook is the comparative view across all nine case studies.
#
# **Learning objectives**
#
# - For each case study, read the highest-IC GBM configuration's average daily
#   Spearman IC with HAC interval on the primary label
# - Compare GBM design choices - loss function, tree depth, optimal iteration
#   count - and locate the operating regime that achieves the highest IC
# - Inspect per-fold IC distributions and the validation to holdout decay where
#   the holdout retrain has been run
# - Compare full-coverage GBM and linear daily-IC point estimates without
#   treating fold summaries as an uncertainty estimator, faceted across labels
# - Extend the comparison across labels (horizon view) and across the
#   classification ↔ regression metric symmetry
# - Inspect feature-importance rank shift versus the linear baseline, per-fold rank
#   stability, and the TabM-vs-GBM-vs-linear three-way picture
#
# **Book reference**: Section 12.6 - Case study insights.
#
# **Prerequisites**: each case study's `07_gbm.py` pipeline has populated
# `run_log/registry.db` for the GBM family. Where present, `tabular_dl.py`
# adds TabM rows. Teaching notebooks NB01-NB11 cover the underlying techniques.

# %%
"""Case Study Insights: Gradient Boosting cross-case-study registry aggregation.

Every comparison below ranks only configurations that covered the same folds and the
same number of days.
"""

import sqlite3
import warnings

import joblib
import matplotlib.pyplot as plt
import numpy as np
import polars as pl

# ml4t.diagnostic dlopens cudart; load torch first so its bundled CUDA
# runtime wins. Same precedence pattern as case_studies/utils/model_analysis.py.
import torch  # noqa: F401
from IPython.display import Markdown, display
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.lines import Line2D
from scipy.stats import rankdata
from sklearn.exceptions import InconsistentVersionWarning
from sklearn.metrics import roc_auc_score

from case_studies.utils.analytics import (
    CASE_STUDY_IDS,
    PRIMARY_LABELS,
    SHORT_NAMES,
)
from case_studies.utils.insight_chapter import (
    SYMMETRY_TABLE_SCHEMA,
    collect_fold_ic_per_cs,
    collect_gbm_checkpoint_trajectories,
    collect_grid_per_cs,
    collect_multi_label_per_cs,
    collect_rank1_per_cs,
    discover_symmetry_pairs,
    load_gbm_feature_importance,
    parse_gbm_config,
    plot_cross_cs_forest,
)
from case_studies.utils.model_analysis import (
    load_metrics_from_registry,
    load_predictions,
)
from utils.paths import get_case_study_dir
from utils.reproducibility import set_global_seeds
from utils.style import COLORS, ml4t_diverging, ml4t_palette, show_with_alt

# LightGBM records synthetic feature names when fitted on an array with an eval_set,
# and sklearn then warns at every predict on an array that has none to compare. One
# message, not the category: the fit and the predictions are unaffected.
warnings.filterwarnings(
    "ignore",
    message="X does not have valid feature names",
    category=UserWarning,
    module="sklearn.utils.validation",
)

# %% tags=["parameters"]
FAMILY = "gbm"
BASELINE_FAMILY = "linear"
SEED = 42


# %%
set_global_seeds(SEED)

# %% [markdown]
# ## 1. Scope and Coverage
#
# The GBM grid spans four tree-depth profiles
# (7 / 15 / 31 / 63 leaves) × 3 regression loss functions (MSE / MAE /
# Huber), evaluated at 10 boosting checkpoints per configuration. Direction
# labels add a binary-logistic variant. The headline metric is average daily
# cross-sectional Spearman IC with a HAC confidence interval on the primary label
# (`prediction_metrics.ic_mean_daily`, `ic_ci_lo`, `ic_ci_hi`, `ic_t_hac`).
# The linear family (Ch11) and TabM (`tabular_dl`) are loaded as baselines.

# %%
coverage_rows = []
for cs in CASE_STUDY_IDS:
    primary = PRIMARY_LABELS[cs]
    gbm = load_metrics_from_registry(cs, families=[FAMILY])
    tabm = load_metrics_from_registry(cs, families=["tabular_dl"])
    lin = load_metrics_from_registry(cs, families=[BASELINE_FAMILY])
    coverage_rows.append(
        {
            "case_study": SHORT_NAMES[cs],
            "primary_label": primary,
            "gbm_labels": gbm["label"].n_unique() if not gbm.is_empty() else 0,
            "gbm_configs_primary": (
                gbm.filter(pl.col("label") == primary)["config_name"].n_unique()
                if not gbm.is_empty()
                else 0
            ),
            "tabm_configs": tabm["config_name"].n_unique() if not tabm.is_empty() else 0,
            "linear_present": not lin.is_empty(),
        }
    )

coverage_df = pl.DataFrame(coverage_rows)
print("GBM / TabM / Linear coverage per case study (primary label and labels trained):")
coverage_df

# %% [markdown]
# ## 2. Cross-CS Forest of Highest-IC GBM Configurations
#
# For each case study, the GBM configuration with the highest average daily IC
# on the primary label is plotted with its HAC interval. Filled markers
# indicate $|t_{HAC}| > 2$ (CI excludes zero); open markers indicate the CI
# overlaps zero.

# %%
gbm_rank1 = collect_rank1_per_cs(
    CASE_STUDY_IDS,
    family=FAMILY,
)
print(
    "Highest-IC GBM configuration per case study (primary label, average daily IC ± HAC 95 % CI):"
)
gbm_rank1.select(
    "short_name",
    "label",
    "config_name",
    pl.col("checkpoint_value").cast(pl.Int64).alias("trees"),
    pl.col("ic_mean_daily").round(4).alias("ic"),
    pl.col("ic_ci_lo").round(4).alias("ci_lo"),
    pl.col("ic_ci_hi").round(4).alias("ci_hi"),
    pl.col("ic_t_hac").round(2).alias("t_hac"),
    pl.col("ic_n_days").cast(pl.Int64).alias("n_days"),
)

# %%
fig, forest_ax = plot_cross_cs_forest(
    gbm_rank1,
    family=FAMILY,
    title="Highest-IC GBM per case study, with HAC intervals",
)
forest_ax.set_xlabel("Average daily IC (HAC 95 % CI)")
show_with_alt(
    fig,
    "Forest plot of average daily IC with a HAC interval, one row per case study on its "
    "primary label, against a vertical line at zero. A filled marker is an interval that "
    "excludes zero and an open one an interval that does not.",
)

# %% [markdown]
# The next cell derives the cross-case conclusion from the selected rows, naming
# no case study or value that it has not just computed.

# %%
clear_zero = gbm_rank1.filter((pl.col("ic_ci_lo") > 0) | (pl.col("ic_ci_hi") < 0))[
    "short_name"
].to_list()
overlap_zero = gbm_rank1.filter((pl.col("ic_ci_lo") <= 0) & (pl.col("ic_ci_hi") >= 0))[
    "short_name"
].to_list()
display(
    Markdown(
        "**Computed reading.** "
        f"The GBM HAC interval excludes zero for {len(clear_zero)} of "
        f"{gbm_rank1.height} case studies ({', '.join(clear_zero) or 'none'}). "
        f"It overlaps zero for {', '.join(overlap_zero) or 'none'}."
    )
)

# %% [markdown]
# ## 3. Within-Family Comparison
#
# How much does each design knob change the IC at the primary label? Three
# subsections trace the structured grid: loss function (3a), tree depth (3b),
# and checkpoint dynamics (3c). Each is presented with HAC intervals so the
# resolution between operating points stays explicit.


# %%
grid_primary = collect_grid_per_cs(
    CASE_STUDY_IDS,
    FAMILY,
    config_parser=parse_gbm_config,
)
grid_regression = grid_primary.filter(pl.col("objective_kind") == "regression")
print(
    f"Per-(CS, config) GBM grid: {grid_primary.height} rows total, "
    f"{grid_regression.height} regression rows."
)

# %% [markdown]
# ### 3a. Loss function (MSE / MAE / Huber)
#
# For each case study, the highest-IC configuration is selected within each
# regression loss family, then the three (case study, loss) IC values are plotted
# side by side with HAC interval bars.

# %%
loss_best = (
    grid_regression.sort("ic_mean_daily", descending=True, nulls_last=True)
    .unique(subset=["case_study", "loss"], keep="first")
    .filter(pl.col("loss").is_in(["mse", "mae", "huber"]))
)

# %%
losses_present = ["mse", "mae", "huber"]
loss_colors = {
    "mse": COLORS["neutral"],
    "mae": COLORS["blue"],
    "huber": COLORS["copper"],
}
cs_order = sorted(loss_best["short_name"].unique().to_list())

# %% [markdown]
# The helper adds bars and intervals without rendering an unfinished figure
# between notebook cells.


# %%
def add_loss_bars(ax: plt.Axes, x: np.ndarray, width: float) -> None:
    for i, loss in enumerate(losses_present):
        sub = loss_best.filter(pl.col("loss") == loss)
        ic, err_lo, err_hi = [], [], []
        for cs in cs_order:
            row = sub.filter(pl.col("short_name") == cs)
            if row.height == 0:
                ic.append(np.nan)
                err_lo.append(0.0)
                err_hi.append(0.0)
            else:
                r = row.row(0, named=True)
                ic.append(r["ic_mean_daily"])
                err_lo.append(r["ic_mean_daily"] - r["ic_ci_lo"])
                err_hi.append(r["ic_ci_hi"] - r["ic_mean_daily"])
        ax.bar(
            x + (i - 1) * width,
            np.array(ic, dtype=float),
            width=width,
            yerr=np.vstack([err_lo, err_hi]),
            capsize=2,
            color=loss_colors[loss],
            alpha=0.9,
            label=loss.upper(),
        )


# %% [markdown]
# Shared labels and a zero line make the grouped confidence intervals comparable.

# %%
fig, ax = plt.subplots(figsize=(11, 5))
x = np.arange(len(cs_order))
width = 0.26
add_loss_bars(ax, x, width)
ax.set_xticks(x)
ax.set_xticklabels(cs_order, rotation=35, ha="right")
ax.axhline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
ax.set_ylabel("Average daily IC (HAC 95 % CI)")
ax.set_title("Highest-IC GBM by loss function and case study")
ax.legend(frameon=False, fontsize=9, loc="best")
show_with_alt(
    fig,
    "Grouped bars of average daily IC, one group per case study and one bar per loss "
    "function, against a line at zero.",
)

# %%
# Aggregate: which loss achieves the highest of the three on each CS?
loss_top_per_cs = (
    loss_best.sort("ic_mean_daily", descending=True, nulls_last=True)
    .group_by("short_name", maintain_order=True)
    .first()
    .group_by("loss", maintain_order=True)
    .len()
    .rename({"len": "n_cs_with_highest_ic"})
    .sort(["n_cs_with_highest_ic", "loss"], descending=[True, False])
)
print(
    "Loss function achieving the highest IC per case study (count across regression-primary CSs):"
)
loss_top_per_cs

# %%
loss_count_text = ", ".join(
    f"{row['loss'].upper()}: {row['n_cs_with_highest_ic']}"
    for row in loss_top_per_cs.iter_rows(named=True)
)
display(
    Markdown(
        f"**Computed loss comparison.** Highest-IC loss counts are {loss_count_text}. "
        "The plotted HAC intervals show whether within-panel loss differences are resolved."
    )
)

# %% [markdown]
# ### 3b. Tree depth heatmap
#
# Within each case study, the highest IC achieved by each leaf profile -
# 7 / 15 / 31 / 63 leaves - is shown as a heatmap cell. Panels with a
# clean diagonal favor a specific depth; panels that are nearly flat across
# leaves indicate the depth knob has no resolution at this signal-to-noise
# ratio.

# %%
depth_pivot = (
    grid_regression.filter(pl.col("leaves").is_not_null())
    .group_by(["short_name", "leaves"], maintain_order=True)
    .agg(pl.col("ic_mean_daily").max().alias("ic"))
    .sort(["short_name", "leaves"])
)
ic_matrix = depth_pivot.pivot(index="short_name", on="leaves", values="ic").sort("short_name")
leaf_cols = sorted(
    [c for c in ic_matrix.columns if c != "short_name"],
    key=lambda s: int(s),
)
matrix_values = ic_matrix.select(leaf_cols).to_numpy()
cs_labels = ic_matrix["short_name"].to_list()

fig, ax = plt.subplots(figsize=(7.5, 5))
vmax = float(np.nanmax(np.abs(matrix_values))) if np.isfinite(matrix_values).any() else 0.05
diverging_cmap = LinearSegmentedColormap.from_list("ml4t_diverging", ml4t_diverging())
im = ax.imshow(matrix_values, cmap=diverging_cmap, vmin=-vmax, vmax=vmax, aspect="auto")
ax.set_xticks(np.arange(len(leaf_cols)))
ax.set_xticklabels([f"{c} leaves" for c in leaf_cols])
ax.set_yticks(np.arange(len(cs_labels)))
ax.set_yticklabels(cs_labels)
for i in range(len(cs_labels)):
    for j in range(len(leaf_cols)):
        v = matrix_values[i, j]
        if np.isfinite(v):
            ax.text(
                j,
                i,
                f"{v:+.3f}",
                ha="center",
                va="center",
                fontsize=8,
                color=COLORS["silver"] if abs(v) > 0.6 * vmax else COLORS["neutral"],
            )
ax.set_title("Highest-IC GBM by leaf profile and case study")
fig.colorbar(im, ax=ax, fraction=0.045, pad=0.04, label="Average daily IC")
show_with_alt(
    fig,
    "Heatmap of average daily IC with one row per case study and one column per leaf "
    "count, each cell shaded by its value against the colour bar.",
)

# %%
depth_spread = (
    depth_pivot.group_by("short_name", maintain_order=True)
    .agg(
        min_ic=pl.col("ic").min(),
        max_ic=pl.col("ic").max(),
        n_leaf_profiles=pl.col("leaves").n_unique(),
    )
    .with_columns(spread=pl.col("max_ic") - pl.col("min_ic"))
    .sort("spread", descending=True)
)
widest_depth = depth_spread.row(0, named=True)
display(
    Markdown(
        f"**Computed depth comparison.** {widest_depth['short_name']} has the widest "
        f"leaf-profile range in this snapshot ({widest_depth['spread']:+.4f} IC across "
        f"{widest_depth['n_leaf_profiles']} profiles). Read each range against that panel's "
        "HAC interval rather than treating the largest point estimate as a resolved difference."
    )
)

# %% [markdown]
# ### 3c. Checkpoint dynamics
#
# Where in the boosting trajectory does the highest-IC checkpoint sit? For
# each case study we plot the highest-validation-IC GBM configuration's mean
# cross-sectional IC against the boosting checkpoint (number of trees,
# 50–500). The trajectory is read from `learning_curves.parquet` written by
# the boosting runner during training; the final early-stopped checkpoint is
# the IC that propagates to the selected row in `prediction_metrics`.

# %%
ckpt_df = collect_gbm_checkpoint_trajectories(gbm_rank1)
if ckpt_df.is_empty() or "short_name" not in ckpt_df.columns:
    msg = (
        "collect_gbm_checkpoint_trajectories returned no rows - "
        "learning_curves.parquet is missing for every case study. "
        "Re-run the GBM boosting sweep before this notebook."
    )
    raise RuntimeError(msg)
print(
    f"Per-checkpoint IC for {ckpt_df['short_name'].n_unique()} case studies "
    f"({ckpt_df.height} (CS, iteration) rows)."
)

# %%
# Ordered by the iteration of each case study's peak IC, ties by peak magnitude.
peak_table = (
    ckpt_df.group_by("short_name", maintain_order=True)
    .agg(
        pl.col("iteration")
        .filter(pl.col("ic_mean") == pl.col("ic_mean").max())
        .first()
        .alias("argmax_iter"),
        pl.col("ic_mean").max().alias("max_ic"),
    )
    .sort(["argmax_iter", "max_ic"], descending=[False, True])
)
cs_order = peak_table["short_name"].to_list()

# %% [markdown]
# Small multiples keep each trajectory on its own axis while sharing the tree budget.

# %%
fig, axes = plt.subplots(3, 3, figsize=(11, 7.5), sharex=True)
for i, cs in enumerate(cs_order):
    ax = axes.flat[i]
    sub = ckpt_df.filter(pl.col("short_name") == cs).sort("iteration")
    cfg = sub["config_name"].first()
    x = sub["iteration"].to_numpy()
    ic = sub["ic_mean"].to_numpy()
    color = COLORS["blue"]
    argmax = int(x[ic.argmax()])
    ax.plot(x, ic, "o-", color=color, linewidth=1.6, markersize=4)
    ax.axvline(argmax, color=color, linewidth=0.6, linestyle=":", alpha=0.7)
    ax.axhline(0, color=COLORS["neutral"], linewidth=0.6, linestyle="--")
    ax.set_title(f"{cs}  ({cfg.replace('_', ' ')})", fontsize=9)
    if i % 3 == 0:
        ax.set_ylabel("IC (mean across folds)")
    if i // 3 == 2:
        ax.set_xlabel("Boosting iteration")
min_peak = int(peak_table["argmax_iter"].min())
max_peak = int(peak_table["argmax_iter"].max())
fig.suptitle("Validation IC against boosting iteration, per case study", fontsize=10)
show_with_alt(
    fig,
    "One small panel per case study, each showing validation IC against boosting "
    "iteration for that case study's selected configuration, with the peak marked.",
)

# %%
peak_table

# %%
early_peaks = peak_table.filter(pl.col("argmax_iter") <= 150)["short_name"].to_list()
budget_peaks = peak_table.filter(pl.col("argmax_iter") == max_peak)["short_name"].to_list()
display(
    Markdown(
        f"**Computed checkpoint reading.** {len(early_peaks)} of {peak_table.height} panels "
        f"peak by 150 trees ({', '.join(early_peaks) or 'none'}). "
        f"{len(budget_peaks)} peak at the observed {max_peak}-tree boundary "
        f"({', '.join(budget_peaks) or 'none'})."
    )
)

# %% [markdown]
# ## 4. Stability and Uncertainty
#
# Average daily IC with HAC CI is the headline metric. Per-fold IC is the
# stability diagnostic, and the validation→holdout decay is the
# generalization diagnostic.

# %% [markdown]
# ### 4a. Per-fold IC distribution
#
# For each case study's highest-IC GBM configuration, the per-fold IC
# distribution is shown as a box-plus-scatter. The reference comparison is
# Ch11 §4 - the linear panels frame the GBM panels' fold-stability picture.

# %%
gbm_fold = collect_fold_ic_per_cs(gbm_rank1)
gbm_fold_summary = (
    gbm_fold.group_by(["case_study", "short_name"], maintain_order=True)
    .agg(
        n_folds=pl.col("ic").count(),
        median=pl.col("ic").median(),
        std=pl.col("ic").std(),
        pct_positive=(pl.col("ic") > 0).mean(),
    )
    .sort("median", descending=True)
)
print("Per-fold IC summary for the highest-IC GBM configuration (primary label):")
gbm_fold_summary

# %%
order = gbm_rank1.sort("ic_mean_daily", descending=True)["short_name"].to_list()
present = [c for c in order if c in gbm_fold["short_name"].unique().to_list()]

fig, ax = plt.subplots(figsize=(11, 4.5))
data = [gbm_fold.filter(pl.col("short_name") == cs)["ic"].to_numpy() for cs in present]
positions = np.arange(len(present))
ax.boxplot(data, positions=positions, widths=0.55, showfliers=True)
for i, arr in enumerate(data):
    if len(arr):
        ax.scatter(np.full(len(arr), i), arr, alpha=0.5, s=14, color=COLORS["blue"])
ax.axhline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
ax.set_xticks(positions)
ax.set_xticklabels(present, rotation=30, ha="right")
ax.set_ylabel("Per-fold Spearman IC")
ax.set_title("Per-fold IC of the highest-IC GBM configuration")
show_with_alt(
    fig,
    "Box plots of per-fold Spearman IC, one box per case study with the individual folds "
    "overlaid as points, against a dashed line at zero.",
)

# %%
gbm_positive_majority = gbm_fold_summary.filter(pl.col("pct_positive") > 0.5)[
    "short_name"
].to_list()
gbm_min_folds = int(gbm_fold_summary["n_folds"].min())
display(
    Markdown(
        f"**Computed fold diagnostic.** {len(gbm_positive_majority)} of "
        f"{gbm_fold_summary.height} selected GBM rows have a positive-fold majority "
        f"({', '.join(gbm_positive_majority) or 'none'}). The smallest exact fold panel "
        f"contains {gbm_min_folds} folds; inference remains attached to the daily HAC series."
    )
)

# %% [markdown]
# ### 4b. Validation → holdout decay dumbbell
#
# Where the holdout retrain has been run, the validation-fold IC and the
# nested-holdout IC are linked by a dumbbell. The case studies without GBM
# holdout rows in the registry are excluded - the gap is explicit, not
# silent.


# %%
HOLDOUT_QUERY = """
    SELECT p.prediction_hash, t.training_hash, t.config_name,
           pm.ic_mean_daily, pm.ic_ci_lo, pm.ic_ci_hi, pm.ic_n_days
    FROM prediction_metrics pm
    JOIN prediction_sets p ON pm.prediction_hash = p.prediction_hash
    JOIN training_runs t ON p.training_hash = t.training_hash
    WHERE t.family = ? AND t.label = ? AND t.config_name = ? AND p.split = 'holdout'
"""


# %%
def load_selected_holdout(case_study: str, family: str, config_name: str) -> dict | None:
    """Load one unique holdout row for the validation-selected configuration."""
    db_path = get_case_study_dir(case_study) / "run_log" / "registry.db"
    db = sqlite3.connect(f"file:{db_path}?mode=ro", uri=True)
    primary = PRIMARY_LABELS[case_study]
    rows = db.execute(
        HOLDOUT_QUERY,
        (family, primary, config_name),
    ).fetchall()
    db.close()
    valid = [
        row
        for row in rows
        if row[6] is not None
        and row[6] > 0
        and all(value is not None and np.isfinite(value) for value in row[3:7])
    ]
    if not valid:
        return None
    if len(valid) != 1:
        raise RuntimeError(f"Ambiguous holdout rows for {case_study}/{family}/{config_name}")
    row = valid[0]
    return {
        "case_study": case_study,
        "short_name": SHORT_NAMES[case_study],
        "label": primary,
        "prediction_hash": row[0],
        "training_hash": row[1],
        "config_name": row[2],
        "holdout_ic": row[3],
        "holdout_ci_lo": row[4],
        "holdout_ci_hi": row[5],
        "holdout_n_days": row[6],
    }


# %%
holdout_rows = []
for selected in gbm_rank1.iter_rows(named=True):
    holdout = load_selected_holdout(
        selected["case_study"],
        FAMILY,
        selected["config_name"],
    )
    if holdout is not None:
        holdout_rows.append(holdout)
holdout_df = pl.DataFrame(holdout_rows) if holdout_rows else pl.DataFrame()
print(
    f"GBM holdout for the validation-selected configuration is available for "
    f"{holdout_df.height} of {gbm_rank1.height} case studies."
)

# %% [markdown]
# Available rows are aligned to the validation-selected configurations before plotting.

# %% [markdown]
# Segments and endpoint intervals show the direction and uncertainty of each change.

# %%
if not holdout_df.is_empty():
    val_lookup = {r["case_study"]: r for r in gbm_rank1.iter_rows(named=True)}
    decay_rows = []
    for r in holdout_rows:
        v = val_lookup.get(r["case_study"])
        if v is None:
            continue
        decay_rows.append(
            {
                "short_name": r["short_name"],
                "val_ic": v["ic_mean_daily"],
                "val_lo": v["ic_ci_lo"],
                "val_hi": v["ic_ci_hi"],
                "holdout_ic": r["holdout_ic"],
                "holdout_lo": r["holdout_ci_lo"],
                "holdout_hi": r["holdout_ci_hi"],
            }
        )
    decay_df = pl.DataFrame(decay_rows).sort("val_ic", descending=True)


# %% [markdown]
# Legend elements distinguish endpoints from the direction of each segment.


# %%
def holdout_legend() -> list[Line2D]:
    return [
        Line2D([0], [0], marker="o", color=COLORS["blue"], label="Validation"),
        Line2D([0], [0], marker="D", color=COLORS["amber"], label="Holdout"),
        Line2D([0], [0], color=COLORS["positive"], linewidth=2, label="Holdout ≥ Val"),
        Line2D([0], [0], color=COLORS["negative"], linewidth=2, label="Holdout < Val"),
    ]


# %% [markdown]
# The complete dumbbell is assembled in one rendering cell so Jupyter never
# captures a partial plot.


# %%
def plot_holdout_decay(decay_df: pl.DataFrame) -> plt.Figure:
    fig, ax = plt.subplots(figsize=(10, max(2.5, 0.5 * decay_df.height + 1)))
    y = np.arange(decay_df.height)
    val, ho = decay_df["val_ic"].to_numpy(), decay_df["holdout_ic"].to_numpy()
    val_lo, val_hi = decay_df["val_lo"].to_numpy(), decay_df["val_hi"].to_numpy()
    ho_lo, ho_hi = decay_df["holdout_lo"].to_numpy(), decay_df["holdout_hi"].to_numpy()
    for i in range(decay_df.height):
        color = COLORS["positive"] if ho[i] >= val[i] else COLORS["negative"]
        ax.plot([val[i], ho[i]], [i, i], color=color, linewidth=2.2, zorder=1)
    ax.errorbar(val, y, xerr=[val - val_lo, val_hi - val], fmt="o", color=COLORS["blue"], capsize=3)
    ax.errorbar(ho, y, xerr=[ho - ho_lo, ho_hi - ho], fmt="D", color=COLORS["amber"], capsize=3)
    ax.set_yticks(y)
    ax.set_yticklabels(decay_df["short_name"].to_list())
    ax.invert_yaxis()
    ax.axvline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
    ax.set_xlabel("Average daily IC (HAC 95 % CI)")
    ax.set_title("Validation and holdout IC for the highest-IC GBM configuration")
    ax.legend(handles=holdout_legend(), loc="best", frameon=False, fontsize=9)
    fig.tight_layout()
    return fig


# %% [markdown]
# The completed chart is rendered only when an exact selected holdout row exists.

# %%
if holdout_df.is_empty():
    display(
        Markdown(
            "**No holdout chart**: no case study in this run has a holdout result for its "
            "selected GBM configuration, so there is nothing to join a validation IC to."
        )
    )
else:
    fig = plot_holdout_decay(decay_df)
    show_with_alt(
        fig,
        "One horizontal line per case study joining its validation IC to its holdout IC, "
        "with a marker at each end, against a dashed line at zero.",
    )

# %%
holdout_names = holdout_df["short_name"].to_list() if not holdout_df.is_empty() else []
display(
    Markdown(
        f"**Computed holdout coverage.** The validation-selected GBM configuration has an "
        f"exact primary-label holdout row for {len(holdout_names)} of {gbm_rank1.height} "
        f"case studies ({', '.join(holdout_names) or 'none'}). No aggregate decay claim is "
        "made when this coverage is sparse."
    )
)

# %% [markdown]
# ## 5. GBM versus Linear
#
# Per case study, how much does GBM raise the average daily IC over the
# strongest full-coverage linear configuration? Two descriptive views compare
# the primary-label point estimates (5a) and the same delta across labels (5b).
# We do not treat a pair of fold summaries as an uncertainty estimator; each
# family's HAC interval comes from its chronological daily IC series.

# %%
linear_rank1 = collect_rank1_per_cs(
    CASE_STUDY_IDS,
    family=BASELINE_FAMILY,
)
linear_fold = collect_fold_ic_per_cs(linear_rank1)


# %% [markdown]
# `select_rank1` makes a family's own candidates comparable - same folds, same number of
# days - but it does that within a family. Two families can each be internally
# comparable and still be scored over different windows, which is what the day counts
# printed below reveal. Subtracting across such a pair charges a window difference to
# the model family, so those pairs are dropped.

# %%
delta_all = (
    gbm_rank1.select(
        "case_study",
        "short_name",
        pl.col("ic_mean_daily").alias("gbm_ic"),
        pl.col("ic_n_days").alias("gbm_days"),
        pl.col("prediction_hash").alias("gbm_prediction_hash"),
    )
    .join(
        linear_rank1.select(
            "case_study",
            pl.col("ic_mean_daily").alias("linear_ic"),
            pl.col("ic_n_days").alias("linear_days"),
            pl.col("prediction_hash").alias("linear_prediction_hash"),
        ),
        on="case_study",
        how="inner",
    )
    .with_columns(delta=pl.col("gbm_ic") - pl.col("linear_ic"))
    .sort("delta", descending=True)
)
delta_excluded = delta_all.filter(pl.col("gbm_days") != pl.col("linear_days"))
delta_primary = delta_all.filter(pl.col("gbm_days") == pl.col("linear_days"))

if not delta_excluded.is_empty():
    display(
        Markdown(
            f"**Excluded for unequal coverage:** {delta_excluded.height} of "
            f"{delta_all.height} case studies have a selected GBM and a selected Linear configuration scored "
            "over a different number of days. Their difference would mix a window "
            "effect into a family effect, so they are left out of the chart below."
        )
    )
    display(delta_excluded.select("short_name", "gbm_days", "linear_days"))

if delta_primary.is_empty():
    display(
        Markdown(
            "**No comparison is possible.** No case study has a selected GBM and a "
            "selected Linear configuration scored over the same number of days, so there "
            "is no matched primary-label delta to report."
        )
    )
else:
    print("Descriptive GBM minus Linear daily-IC delta at matched full coverage:")
    display(
        delta_primary.select(
            "short_name",
            pl.col("gbm_ic").round(4),
            pl.col("linear_ic").round(4),
            pl.col("delta").round(4),
            "gbm_days",
            "linear_days",
        )
    )

# %%
if not delta_primary.is_empty():
    fig, ax = plt.subplots(figsize=(9, 4.5))
    y = np.arange(delta_primary.height)
    delta = delta_primary["delta"].to_numpy()
    colors = [COLORS["blue"] if value >= 0 else COLORS["amber"] for value in delta]
    ax.barh(y, delta, color=colors, alpha=0.9)
    ax.axvline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
    ax.set_yticks(y)
    ax.set_yticklabels(delta_primary["short_name"].to_list())
    ax.invert_yaxis()
    ax.set_xlabel("Average daily IC point-estimate delta (GBM - Linear)")
    ax.set_title("Matched-coverage GBM minus Linear at the primary label")
    show_with_alt(
        fig,
        "Horizontal bars of the difference in average daily IC between the GBM and the "
        "linear model on matched coverage, one bar per case study, against a line at zero.",
    )

# %%
n_positive = delta_primary.filter(pl.col("delta") > 0).height
display(
    Markdown(
        f"**Computed comparison.** GBM has the higher daily-IC point estimate in "
        f"{n_positive} of the {delta_primary.height} case studies whose two selections were "
        f"scored over the same number of days, out of {delta_all.height} with a selection in "
        "both families. This chart is descriptive: the two model families have separate "
        "daily-series HAC intervals, so no paired-fold confidence claim is attached to "
        "their difference."
    )
)

# %% [markdown]
# ### 5b. GBM minus Linear, faceted by label
#
# The primary-label delta is a single horizon per case study. Faceting
# across all labels where both families trained shows whether the GBM
# advantage (or absence thereof) is consistent across horizons within a
# panel.


# %%
HORIZON_EXCLUSIONS = {
    (
        "nasdaq100_microstructure",
        "gbm",
        "fwd_ret_5m",
    ): "two complete GBM candidates tie for rank one in the provisional registry",
    (
        "sp500_equity_option_analytics",
        "linear",
        "fwd_ret_risk_adj_5d",
    ): "five complete Linear candidates tie for rank one",
}


def regression_labels(cs: str, family: str) -> list[str]:
    df = load_metrics_from_registry(cs, families=[family])
    if df.is_empty():
        return []
    return [
        lbl
        for lbl in df["label"].unique().to_list()
        if lbl is not None
        and lbl.startswith("fwd_ret_")
        and "spot" not in lbl
        and (cs, family, lbl) not in HORIZON_EXCLUSIONS
    ]


# %% [markdown]
# Both families pass through the same complete-coverage selector, which makes
# each family's selection comparable against its own alternatives. It does not make
# the two selections comparable to each other, so the join below keeps only the
# cells where both were scored over the same number of days.

# %%
gbm_horizon = collect_multi_label_per_cs(
    CASE_STUDY_IDS,
    family=FAMILY,
    labels=lambda cs: regression_labels(cs, FAMILY),
)
lin_horizon = collect_multi_label_per_cs(
    CASE_STUDY_IDS,
    family=BASELINE_FAMILY,
    labels=lambda cs: regression_labels(cs, BASELINE_FAMILY),
)

# %% [markdown]
# Two optional family-label cells are ambiguous: NASDAQ-100 GBM
# `fwd_ret_5m` has two complete rank-one candidates, while the S&P
# equity-option Linear `fwd_ret_risk_adj_5d` cell has five. Those family-label
# cells are excluded rather than resolved arbitrarily. Primary labels and all
# unambiguous family-label cells remain in the notebook.

# %% [markdown]
# Matched case-study and label rows form the descriptive family-difference panel.

# %%
facet_all = (
    gbm_horizon.select(
        "case_study",
        "short_name",
        "label",
        pl.col("ic_mean_daily").alias("gbm_ic"),
        pl.col("ic_n_days").alias("gbm_days"),
    )
    .join(
        lin_horizon.select(
            "case_study",
            "label",
            pl.col("ic_mean_daily").alias("lin_ic"),
            pl.col("ic_n_days").alias("lin_days"),
        ),
        on=["case_study", "label"],
        how="inner",
    )
    .with_columns(delta=pl.col("gbm_ic") - pl.col("lin_ic"))
    .sort(["short_name", "label"])
)
# Same cross-family coverage rule as 5a: each family's selection is comparable
# within its own family, which does not make the two comparable to each other.
facet_excluded = facet_all.filter(pl.col("gbm_days") != pl.col("lin_days"))
facet_df = facet_all.filter(pl.col("gbm_days") == pl.col("lin_days"))

if not facet_excluded.is_empty():
    display(
        Markdown(
            f"**Excluded for unequal coverage:** {facet_excluded.height} of "
            f"{facet_all.height} (case study, label) cells pair selections scored over a "
            "different number of days."
        )
    )
    display(facet_excluded.select("short_name", "label", "gbm_days", "lin_days"))

if facet_df.is_empty():
    display(
        Markdown(
            "**No comparison is possible.** No (case study, label) cell pairs a selected "
            "GBM and a selected Linear configuration scored over the same number of days."
        )
    )
else:
    print(f"GBM-minus-Linear deltas across {facet_df.height} matched-coverage (CS, label) cells:")
    display(
        facet_df.select(
            "short_name",
            "label",
            pl.col("gbm_ic").round(4).alias("gbm"),
            pl.col("lin_ic").round(4).alias("lin"),
            pl.col("delta").round(4).alias("delta"),
            "gbm_days",
            "lin_days",
        )
    )

# %%
if not facet_df.is_empty():
    cs_sorted = sorted(facet_df["short_name"].unique().to_list())
    fig, ax = plt.subplots(figsize=(11, max(3.5, 0.45 * facet_df.height + 1.5)))
    y_off = 0
    yticks, ylabels = [], []
    for cs in cs_sorted:
        sub = facet_df.filter(pl.col("short_name") == cs)
        for r in sub.iter_rows(named=True):
            color = COLORS["blue"] if r["delta"] >= 0 else COLORS["amber"]
            ax.barh(y_off, r["delta"], color=color, height=0.65, alpha=0.9)
            yticks.append(y_off)
            ylabels.append(f"{cs} · {r['label']}")
            y_off += 1
        y_off += 0.5  # gap between case studies
    ax.set_yticks(yticks)
    ax.set_yticklabels(ylabels, fontsize=8)
    ax.invert_yaxis()
    ax.axvline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
    ax.set_xlabel("Average daily IC delta (GBM − Linear)")
    ax.set_title("GBM minus Linear, by label within each case study")
    show_with_alt(
        fig,
        "Horizontal bars of the IC difference between the GBM and the linear model, one "
        "bar per label grouped by case study, against a line at zero.",
    )

# %%
if facet_df.is_empty():
    display(
        Markdown(
            "**No horizon comparison.** Every (case study, label) cell paired selections "
            "scored over a different number of days, so there is no matched range to "
            "report."
        )
    )
else:
    facet_positive = facet_df.filter(pl.col("delta") > 0).height
    largest_facet = facet_df.sort("delta", descending=True).row(0, named=True)
    smallest_facet = facet_df.sort("delta").row(0, named=True)
    display(
        Markdown(
            f"**Computed horizon comparison.** GBM has the higher point estimate in "
            f"{facet_positive} of {facet_df.height} matched case-study/label cells. The range "
            f"runs from {smallest_facet['short_name']} {smallest_facet['label']} "
            f"({smallest_facet['delta']:+.4f}) to {largest_facet['short_name']} "
            f"{largest_facet['label']} ({largest_facet['delta']:+.4f})."
        )
    )

# %% [markdown]
# ## 6. Multi-Label Horizon and Metric Symmetry

# %% [markdown]
# ### 6a. GBM IC across regression labels per case study

# %%
HORIZON_DAYS = {
    "fwd_ret_5m": 5 / (6.5 * 60),
    "fwd_ret_15m": 15 / (6.5 * 60),
    "fwd_ret_60m": 60 / (6.5 * 60),
    "fwd_ret_8h": 1.0 / 3,
    "fwd_ret_24h": 1.0,
    "fwd_ret_1d": 1.0,
    "fwd_ret_5d": 5.0,
    "fwd_ret_10d": 10.0,
    "fwd_ret_21d": 21.0,
    "fwd_ret_1m": 21.0,
    "fwd_ret_3m": 63.0,
    "fwd_ret_1m_win": 21.0,
    "fwd_ret_risk_adj_5d": 5.0,
}

# %% [markdown]
# Only panels with at least two mapped horizons enter the log-scale comparison.
#
# The sibling figure in Chapter 13 also drops cells whose scored folds cover only part of
# the case study's modelling grid, because joining a partial-grid point to a full-grid one
# on a shared IC axis reads as one quantity moving with horizon. `gbm_horizon` carries the
# same `covers_fold_grid` column and this cell does not filter on it, because no gbm or
# linear run can be partial: `us_equities_panel/12_dl_weekly` is the only notebook in the
# corpus with a non-zero `MAX_FOLDS`, and it is a deep-learning notebook. Give any gbm or
# linear notebook a fold reduction and this cell needs Chapter 13's retain-and-report
# split before it is read again.

# %%
plot_horizon = gbm_horizon.with_columns(
    horizon_days=pl.col("label").replace_strict(HORIZON_DAYS, default=None).cast(pl.Float64),
).filter(pl.col("horizon_days").is_not_null())
multi_cs = (
    plot_horizon.group_by("short_name", maintain_order=True)
    .len()
    .filter(pl.col("len") >= 2)["short_name"]
    .to_list()
)
plot_horizon = plot_horizon.filter(pl.col("short_name").is_in(multi_cs))

# %% [markdown]
# Each line carries its own HAC band; marker and line style supplement the compact palette.

# %%
if plot_horizon.height > 0:
    fig, ax = plt.subplots(figsize=(10, 5))
    palette = ml4t_palette(5, categorical=True)
    cs_sorted = sorted(plot_horizon["short_name"].unique().to_list())
    markers = ["o", "s", "D", "^", "v", "P", "X", "*"]
    linestyles = ["-", "--", "-.", ":", "-", "--", "-.", ":"]
    for idx, cs in enumerate(cs_sorted):
        # Sorted on (horizon_days, label): two targets can share a horizon, and their tie
        # order decides which point the line reaches first, so sorting on the horizon
        # alone made the figure differ run to run.
        sub = plot_horizon.filter(pl.col("short_name") == cs).sort(["horizon_days", "label"])
        if sub.height < 2:
            continue
        x = sub["horizon_days"].to_numpy()
        ic = sub["ic_mean_daily"].to_numpy()
        lo = sub["ic_ci_lo"].to_numpy()
        hi = sub["ic_ci_hi"].to_numpy()
        color = palette[idx % len(palette)]
        ax.fill_between(x, lo, hi, color=color, alpha=0.12)
        ax.plot(
            x,
            ic,
            marker=markers[idx % len(markers)],
            linestyle=linestyles[idx % len(linestyles)],
            color=color,
            label=cs,
            linewidth=1.6,
            markersize=6,
            alpha=0.9,
        )
    ax.set_xscale("log")
    ax.set_xlabel("Horizon (trading days, log scale)")
    ax.set_ylabel("Average daily IC (HAC 95 % CI band)")
    ax.axhline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
    ax.set_title("Highest-IC GBM configuration across regression horizons")
    ax.legend(loc="best", frameon=False, fontsize=8, ncol=2)
    show_with_alt(
        fig,
        "Average daily IC against label horizon on a log axis, one line per case study, "
        "against a line at zero.",
    )

# %%
horizon_ranges = (
    plot_horizon.group_by("short_name", maintain_order=True)
    .agg(
        n_horizons=pl.len(),
        min_ic=pl.col("ic_mean_daily").min(),
        max_ic=pl.col("ic_mean_daily").max(),
    )
    .with_columns(ic_range=pl.col("max_ic") - pl.col("min_ic"))
    .sort("ic_range", descending=True)
)
display(
    Markdown(
        f"**Computed horizon coverage.** {horizon_ranges.height} case studies have at least "
        "two comparable regression horizons. Their exact IC ranges are shown above; no "
        "cross-horizon trend is asserted for single-horizon panels."
    )
)

# %% [markdown]
# ### 6b. Classification ↔ regression metric symmetry
#
# A GBM trained on a continuous return label has a regression score; the
# same model class can be trained on the binary direction label of the
# same horizon. Two symmetric questions arise on the case studies that
# carry binary direction labels paired to the regression labels:
#
# - **Direction A** - the GBM classification model's score, evaluated as
#   IC against the *continuous* return, asks whether the directional
#   classifier is also a useful *cross-sectional ranker*. Read from
#   `prediction_metrics.ic_mean_daily` for `task_type='classification'` rows.
# - **Direction B** - the GBM regression model's score, evaluated as AUC
#   against the *binary* direction, asks whether the continuous regression
#   score is also a useful *binary classifier*. Computed on the fly here
#   from raw OOF predictions.
#
# The pairs are discovered from the registry rather than listed here. A case study
# qualifies when this family has registered validation predictions for a
# `fwd_ret_<horizon>` label and for a `fwd_dir_<horizon>` or `fwd_class_<horizon>`
# label at the same horizon, and the direction label's own surface is binary. A
# ternary label needs a multiclass score and is out of scope, so it is skipped by
# measuring its domain and naming it, not by being absent from a list: the previous
# hand-written literal had lost `us_firm_characteristics` and reported a full count
# of itself, and the comment that stood here called `fwd_class_1m` ternary when its
# domain is {0, 1}. This mirrors Ch11 §6b for the linear family.

# %%
SYMMETRY_PAIRS, SYMMETRY_SKIPS = discover_symmetry_pairs(CASE_STUDY_IDS, FAMILY)
for line in SYMMETRY_SKIPS:
    print(f"  skipped {line}")
print(
    f"{sum(len(pairs) for pairs in SYMMETRY_PAIRS.values())} matched label pairs across "
    f"{len(SYMMETRY_PAIRS)} case studies"
)


# %% [markdown]
# Binary labels align on canonical timestamp and symbol before the regression
# score is evaluated as a direction classifier.

# %%


def load_binary_label(cs: str, dir_label: str) -> pl.DataFrame:
    p = get_case_study_dir(cs) / "labels" / f"{dir_label}.parquet"
    if not p.exists():
        return pl.DataFrame()
    df = pl.read_parquet(p)
    return df.rename({dir_label: "y_dir"}).select("timestamp", "symbol", "y_dir")


# %%
def gbm_direction_b_auc(selected: dict, dir_label: str) -> dict | None:
    cs = selected["case_study"]
    reg_label = selected["label"]
    preds = load_predictions(
        cs,
        family=FAMILY,
        label=reg_label,
        config_name=selected["config_name"],
        checkpoint_value=selected["checkpoint_value"],
        split="validation",
    ).filter(pl.col("prediction_hash") == selected["prediction_hash"])
    if preds.height == 0:
        return None
    dir_df = load_binary_label(cs, dir_label)
    if dir_df.is_empty():
        return None
    canonical_types = {"timestamp": pl.Datetime("ms"), "symbol": pl.Utf8}
    for column, dtype in canonical_types.items():
        if preds[column].dtype != dir_df[column].dtype:
            preds = preds.with_columns(pl.col(column).cast(dtype))
            dir_df = dir_df.with_columns(pl.col(column).cast(dtype))
    merged = preds.join(dir_df, on=["timestamp", "symbol"], how="inner")
    domain = set(merged["y_dir"].unique().drop_nulls().to_list())
    if not domain.issubset({0, 1}):
        print(f"  SKIP {cs}/{dir_label}: y_dir domain {domain} is not binary {{0,1}}")
        return None
    merged = merged.filter(pl.col("y_dir").is_in([0, 1]) & pl.col("y_score").is_not_null())
    if merged.height == 0 or merged["y_dir"].n_unique() < 2:
        return None
    auc = float(roc_auc_score(merged["y_dir"].to_numpy(), merged["y_score"].to_numpy()))
    return {
        "case_study": cs,
        "short_name": SHORT_NAMES[cs],
        "reg_label": reg_label,
        "dir_label": dir_label,
        "reg_config": selected["config_name"],
        "reg_prediction_hash": selected["prediction_hash"],
        "reg_score_auc": auc,
        "n": merged.height,
    }


# %%
direction_labels = {
    cs: [direction_label for _, direction_label in pairs] for cs, pairs in SYMMETRY_PAIRS.items()
}
direction_rank1 = collect_multi_label_per_cs(
    SYMMETRY_PAIRS,
    family=FAMILY,
    labels=lambda cs: direction_labels[cs],
)

# %% [markdown]
# Direction A comes from the selected classification row; Direction B uses the
# exact selected regression prediction hash against the paired binary label.

# %%
sym_rows = []
for cs, pairs in SYMMETRY_PAIRS.items():
    for reg_lbl, dir_lbl in pairs:
        cls_ic = cls_lo = cls_hi = cls_t = cls_cfg = None
        cls_auc = cls_auc_lo = cls_auc_hi = None
        cls_selected = direction_rank1.filter(
            (pl.col("case_study") == cs) & (pl.col("label") == dir_lbl)
        )
        if not cls_selected.is_empty():
            top = cls_selected.row(0, named=True)
            cls_ic = top["ic_mean_daily"]
            cls_lo = top.get("ic_ci_lo")
            cls_hi = top.get("ic_ci_hi")
            cls_t = top.get("ic_t_hac")
            cls_cfg = top["config_name"]
            # The classifier's AUC against its own direction label. It is the third
            # column of the book's Table 12.4 and sits in the same prediction_metrics
            # row as the IC above; until 2026-09-18 the selector did not read it, so
            # this notebook published a table it could not reproduce from its own output.
            cls_auc = top.get("auc_mean_daily")
            cls_auc_lo = top.get("auc_ci_lo")
            cls_auc_hi = top.get("auc_ci_hi")
        reg_selected = gbm_horizon.filter(
            (pl.col("case_study") == cs) & (pl.col("label") == reg_lbl)
        )
        b = (
            gbm_direction_b_auc(reg_selected.row(0, named=True), dir_lbl)
            if not reg_selected.is_empty()
            else None
        )
        sym_rows.append(
            {
                "short_name": SHORT_NAMES[cs],
                "reg_label": reg_lbl,
                "dir_label": dir_lbl,
                "cls_config": cls_cfg,
                "cls_score_ic": cls_ic,
                "cls_score_ic_lo": cls_lo,
                "cls_score_ic_hi": cls_hi,
                "cls_score_ic_t": cls_t,
                "cls_score_auc": cls_auc,
                "cls_score_auc_lo": cls_auc_lo,
                "cls_score_auc_hi": cls_auc_hi,
                "reg_config": (b or {}).get("reg_config"),
                "reg_score_auc": (b or {}).get("reg_score_auc"),
                "n_b": (b or {}).get("n"),
            }
        )

# %% [markdown]
# The combined table keeps both metric directions and their selected identities visible.

# %%

# A full schema, not schema_overrides. Discovery can legitimately return nothing - every
# declared pair skipped, or a registry with no classification runs for this family - and a
# frame built from an empty list with partial overrides has no string columns at all, so
# the selection below raises ColumnNotFoundError and the skip reasons this section exists
# to print never reach the reader. An empty frame with the right columns renders as an
# empty table, which is the correct answer.
sym_df = pl.DataFrame(sym_rows, schema=SYMMETRY_TABLE_SCHEMA)
print(
    "Direction A (GBM classification score → IC), the classifier's own AUC, "
    "and Direction B (GBM regression score → AUC):"
)
sym_df.select(
    "short_name",
    "reg_label",
    "dir_label",
    pl.col("cls_score_ic").round(4).alias("A_ic"),
    pl.col("cls_score_ic_lo").round(4).alias("A_lo"),
    pl.col("cls_score_ic_hi").round(4).alias("A_hi"),
    pl.col("cls_score_ic_t").round(2).alias("A_t"),
    pl.col("cls_score_auc").round(4).alias("native_auc"),
    pl.col("cls_score_auc_lo").round(4).alias("native_lo"),
    pl.col("cls_score_auc_hi").round(4).alias("native_hi"),
    pl.col("reg_score_auc").round(4).alias("B_auc"),
)

# %%
if sym_df.is_empty():
    raise RuntimeError(
        "no declared regression/direction pair qualified, so there is nothing to compare. "
        f"Skipped: {'; '.join(SYMMETRY_SKIPS) or 'nothing'}"
    )

fig, axes = plt.subplots(1, 2, figsize=(13, 4.0))
labels_y = [f"{r['short_name']} · {r['dir_label']}" for r in sym_df.iter_rows(named=True)]
y = np.arange(sym_df.height)

ax = axes[0]
ic = sym_df["cls_score_ic"].to_numpy()
lo = sym_df["cls_score_ic_lo"].to_numpy()
hi = sym_df["cls_score_ic_hi"].to_numpy()
ax.errorbar(
    ic,
    y,
    xerr=[ic - lo, hi - ic],
    fmt="o",
    color=COLORS["blue"],
    capsize=3,
    lw=1,
)
ax.axvline(0, color=COLORS["neutral"], lw=0.7, linestyle="--")
ax.set_yticks(y)
ax.set_yticklabels(labels_y)
ax.invert_yaxis()
ax.set_xlabel("Average daily IC (HAC 95 % CI)")
ax.set_title("Classification score against the continuous return")

ax = axes[1]
auc = sym_df["reg_score_auc"].to_numpy()
ax.scatter(auc, y, color=COLORS["copper"], s=60, zorder=3)
ax.axvline(0.5, color=COLORS["neutral"], lw=0.7, linestyle="--")
ax.set_yticks(y)
ax.set_yticklabels([])
ax.invert_yaxis()
ax.set_xlabel("AUC (regression score, pooled OOF)")
ax.set_title("Regression score against the binary direction")

show_with_alt(
    fig,
    "Two panels. Left: the classification score's IC against the continuous return, one "
    "row per pair with its HAC interval, against a line at zero. Right: the regression "
    "score's pooled out-of-fold AUC against the binary direction, against a line at one "
    "half.",
)

# %%
direction_a_positive = sym_df.filter(pl.col("cls_score_ic") > 0).height
direction_b_valid = sym_df.filter(pl.col("reg_score_auc").is_not_null())
native_valid = sym_df.filter(pl.col("cls_score_auc").is_not_null())
native_clears = native_valid.filter(pl.col("cls_score_auc_lo") > 0.5).height
# A count over the pairs that happen to be present is not a count over the corpus, so the
# skipped candidates are named in the same sentence rather than left to the log above.
skip_note = (
    f", with {len(SYMMETRY_SKIPS)} candidate pair(s) skipped ({'; '.join(SYMMETRY_SKIPS)})"
    if SYMMETRY_SKIPS
    else ", with no candidate pair skipped"
)
max_auc_distance = (
    float((direction_b_valid["reg_score_auc"] - 0.5).abs().max())
    if not direction_b_valid.is_empty()
    else float("nan")
)
display(
    Markdown(
        f"**Computed metric symmetry.** Direction A is positive in {direction_a_positive} of "
        f"{sym_df.height} matched cells, discovered from the registry rather than declared"
        f"{skip_note}. The classifier's own AUC clears one half on "
        f"{native_clears} of {native_valid.height} cells. Direction B is available in "
        f"{direction_b_valid.height} cells, with maximum absolute distance from chance of "
        f"{max_auc_distance:.4f}."
    )
)

# %% [markdown]
# ## 7. Interpretability
#
# Three diagnostics read the fold models each selected run saved: feature-importance
# rank shift against the linear baseline (7a), per-fold rank stability of GBM
# importances (7b), and the TabM-vs-GBM-vs-Linear three-way picture (7c).


# %%
NON_FEATURE_COLS = {"timestamp", "symbol", "stock_id", "product", "position", "instrument_id"}
IMPORTANCE_CASES = ["etfs", "sp500_options", "us_firm_characteristics", "us_equities_panel"]

# A penalized fit drives coefficients to exactly zero, and a feature the model does not
# use has no rank. Ranking those alongside the rest would fill the comparison with ties
# broken by feature name, so they come out and the surviving count is reported.
ZERO_TOL = 1e-12


def _load_linear_importance(cs: str, training_hash: str) -> dict[str, float]:
    """Mean absolute coefficient per feature, over the folds of one linear training run.

    The pipeline stores each fold's fitted estimator under the run's `models`
    directory; the coefficients live on the estimator rather than in a table of
    their own. Features the model zeroed in every fold are dropped.
    """
    models_dir = get_case_study_dir(cs) / "run_log" / "training" / training_hash / "models"
    fold_paths = sorted(models_dir.glob("fold_*.joblib")) if models_dir.exists() else []
    rows = []
    for fold_path in fold_paths:
        with warnings.catch_warnings():
            # Read for `coef_` only, never to predict, so an estimator pickled by a
            # different scikit-learn is safe to unpickle here.
            warnings.filterwarnings(
                "ignore", category=InconsistentVersionWarning, module="sklearn.base"
            )
            payload = joblib.load(fold_path)
        names = list(payload["feature_names"])
        coef = np.asarray(payload["model"].coef_, dtype=float).ravel()
        if coef.size != len(names):
            raise RuntimeError(
                f"{cs} {training_hash}: {fold_path.name} carries {coef.size} coefficients "
                f"against {len(names)} feature names, so the two cannot be paired."
            )
        rows.extend(
            {"feature": name, "abs_coef": abs(float(value))}
            for name, value in zip(names, coef, strict=True)
        )
    if not rows:
        return {}
    return dict(
        pl.DataFrame(rows)
        .group_by("feature", maintain_order=True)
        .agg(pl.col("abs_coef").mean())
        .filter(pl.col("abs_coef") > ZERO_TOL)
        .sort("abs_coef", descending=True)
        .iter_rows()
    )


# %% [markdown]
# Rank dictionaries use the same feature intersection for both families.


# %%
def _feature_ranks(
    gbm_imp_df: pl.DataFrame,
    linear_imp: dict[str, float],
) -> tuple[dict[str, int], dict[str, int]]:
    gbm_imp = dict(
        gbm_imp_df.group_by("feature", maintain_order=True)
        .agg(pl.col("importance").mean().alias("imp"))
        .iter_rows()
    )
    # Ties break by feature name rather than by set iteration order: gain importances tie
    # readily, and with per-process string hashing the same registries produced different
    # ranks from one run to the next.
    common = sorted(set(gbm_imp) & set(linear_imp))
    gbm_order = sorted(common, key=lambda feature: (-gbm_imp[feature], feature))
    linear_order = sorted(common, key=lambda feature: (-linear_imp[feature], feature))
    return (
        {feature: rank for rank, feature in enumerate(gbm_order, 1)},
        {feature: rank for rank, feature in enumerate(linear_order, 1)},
    )


# %% [markdown]
# ### 7a. Feature importance - GBM versus the linear baseline
#
# For each case study that saved both boosters and linear fold models, features are
# ranked by GBM gain importance and by the linear model's mean |coefficient|, both
# averaged across folds. The rank shift `linear_rank − gbm_rank` is positive for
# features the GBM promotes over the linear model (typically interaction or regime
# features) and negative for features the linear model promotes (typically monotonic
# predictors). The baseline is whichever linear configuration validation selected, so
# a sparse fit such as a lasso contributes only the features it kept; `n_common_features`
# says how many that leaves.


# %%
def feature_rank_shift(cs: str) -> dict | None:
    gbm_selected = gbm_rank1.filter(pl.col("case_study") == cs)
    linear_selected = linear_rank1.filter(pl.col("case_study") == cs)
    if gbm_selected.is_empty() or linear_selected.is_empty():
        return None
    gbm_row = gbm_selected.row(0, named=True)
    linear_row = linear_selected.row(0, named=True)
    gbm_imp_df = load_gbm_feature_importance(
        cs,
        gbm_row["training_hash"],
        gbm_row["config_name"],
        top_n=50,
        num_iteration=gbm_row["checkpoint_value"],
    )
    if gbm_imp_df.is_empty():
        return None
    linear_imp = _load_linear_importance(cs, linear_row["training_hash"])
    gbm_ranks, linear_ranks = _feature_ranks(gbm_imp_df, linear_imp)
    if not gbm_ranks:
        return None
    shifts = pl.DataFrame(
        [
            {
                "feature": f,
                "gbm_rank": gbm_ranks[f],
                "linear_rank": linear_ranks[f],
                "rank_shift": linear_ranks[f] - gbm_ranks[f],
            }
            for f in gbm_ranks
        ]
    )
    min_shift = shifts["rank_shift"].min()
    return {
        "short_name": SHORT_NAMES[cs],
        "n_common_features": shifts.height,
        "median_abs_shift": float(shifts["rank_shift"].abs().median()),
        "max_gbm_promotion": int(shifts["rank_shift"].max() or 0),
        "max_linear_promotion": int(-min_shift if min_shift else 0),
        "_shifts": shifts,
    }


# %% [markdown]
# Apply the comparison only where both exact selected artifact sets exist.

# %%
rank_shift_summary = [
    entry for cs in IMPORTANCE_CASES if (entry := feature_rank_shift(cs)) is not None
]

print(f"Computed GBM-vs-linear rank shifts for {len(rank_shift_summary)} case studies.")
shift_summary_df = (
    pl.DataFrame(
        [{k: v for k, v in r.items() if not k.startswith("_")} for r in rank_shift_summary]
    )
    if rank_shift_summary
    else pl.DataFrame()
)
shift_summary_df

# %%
if not rank_shift_summary:
    display(
        Markdown(
            "**No promotion chart**: no case study has both a GBM importance artifact and "
            "linear fold models, so there is no pair of rankings to difference."
        )
    )
else:
    # Each panel ranks its own case study's features, so the y axes carry different
    # categories and cannot be shared: one shared axis would draw the last panel's
    # feature names beside every panel's bars.
    fig, axes = plt.subplots(
        1,
        len(rank_shift_summary),
        figsize=(4.5 * len(rank_shift_summary), 4.5),
    )
    if len(rank_shift_summary) == 1:
        axes = [axes]
    for ax, entry in zip(axes, rank_shift_summary, strict=False):
        # `rank_shift` is an integer and features tie on it readily, so head/tail would
        # otherwise pick an arbitrary member of a tied block: which feature names the
        # panel draws would change between executions on identical data.
        s = entry["_shifts"].sort(["rank_shift", "feature"], descending=[True, False])
        n_show = min(15, s.height)
        top_promotions = s.head(n_show // 2)
        bot_promotions = s.tail(n_show - top_promotions.height)
        plot_set = pl.concat([top_promotions, bot_promotions]).sort(
            ["rank_shift", "feature"], descending=[False, True]
        )
        y = np.arange(plot_set.height)
        colors = [
            COLORS["blue"] if v > 0 else COLORS["amber"] for v in plot_set["rank_shift"].to_list()
        ]
        ax.barh(y, plot_set["rank_shift"].to_numpy(), color=colors, height=0.6, alpha=0.9)
        ax.set_yticks(y)
        ax.set_yticklabels(plot_set["feature"].to_list(), fontsize=7)
        ax.axvline(0, color=COLORS["neutral"], linewidth=0.7, linestyle="--")
        ax.set_xlabel("Linear rank − GBM rank (positive = GBM promotion)")
        ax.set_title(entry["short_name"])
    fig.suptitle("Feature rank difference between the linear baseline and the GBM")
    show_with_alt(
        fig,
        "One panel per case study of horizontal bars, each p