The held-out test, and the justification

import textwrap

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from lensml import (
    FEATURES,
    GOLD_CMAP,
    RNG_SEED,
    RUNG0,
    cite,
    credit_own,
    decile_edges,
    ladder,
    rebuild_register,
    setup,
    show_fig,
    wet_by_bin,
)
from sklearn.metrics import (
    ConfusionMatrixDisplay,
    accuracy_score,
    average_precision_score,
    classification_report,
    confusion_matrix,
    precision_score,
    recall_score,
)

pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 30)

S = setup()  # the cache, the test members, the transforms (chapters 1-2)
train, test = S.train, S.test
X_tr, y_tr = train[FEATURES].values, train["wet"].values
g_tr = train["member_id"].values
print(
    f"train: {S.n_train_members} members / {len(train)} winters   "
    f"test: {test.member_id.nunique()} members / {len(test)} winters   "
    f"N_FOLDS = {S.n_folds}"
)

YOUR_TRIALS = []  # the same list as chapter 5's Your-trial cell
reg = rebuild_register(S, YOUR_TRIALS)
print(f"incumbent after {len(reg.rows)} trials: {reg.incumbent_name}")
reg.bonferroni_line()
cols = [
    "trial",
    "candidate",
    "incumbent",
    "recall_gap",
    "called_wet_cand",
    "called_wet_inc",
    "recall_p",
    "precision_gap",
    "precision_p",
    "verdict",
]
reg.table()[cols]

FINAL = True  # the shipped register is final; set False while adding trials

if not FINAL:  # nothing reads the test labels
    print(
        f"test not run: set FINAL = True after your last trial (register so far: "
        f"{len(reg.rows)} tested rows, incumbent {reg.incumbent_name})"
    )
else:
    TEST_RUNS = (
        globals().get("TEST_RUNS", 0) + 1
    )  # a second run in one session is a second test
    X_te, y_te = test[FEATURES].values, test["wet"].values
    final = reg.incumbent_make(RNG_SEED).fit(X_tr, y_tr)
    yhat_te = final.predict(X_te)
    cm_te = confusion_matrix(y_te, yhat_te, labels=[0, 1])
    score_te = final.predict_proba(X_te)[:, 1]
    tn, fp, fn, tp = (int(v) for v in cm_te.ravel())

    print(
        f"test run #{TEST_RUNS}"
        + ("" if TEST_RUNS == 1 else "   — NOT a test: report run #1")
    )
    print(
        f"final model: {reg.incumbent_name}   test members: {test.member_id.nunique()}   "
        f"test winters: {len(y_te)}   base rate: {y_te.mean():.3f}"
    )
    acc_final = accuracy_score(y_te, yhat_te)
    print(
        f"precision (wet) = {precision_score(y_te, yhat_te, zero_division=0):.3f}   "
        f"recall (wet) = {recall_score(y_te, yhat_te, zero_division=0):.3f}   "
        f"wet-call share = {yhat_te.mean():.3f}   accuracy = {acc_final:.3f}   "
        f"AP = {average_precision_score(y_te, score_te):.3f}"
    )
    print(
        classification_report(
            y_te, yhat_te, target_names=["dry", "wet"], zero_division=0
        )
    )
    print(
        "test confusion matrix (rows = true, columns = predicted); the figure below "
        "draws the same matrix:"
    )
    print(
        pd.DataFrame(
            cm_te,
            index=["true dry", "true wet"],
            columns=["pred dry", "pred wet"],
        )
    )
    print(
        f"false alarms (true dry, called wet) {fp} · hits (true wet, called wet) {tp} · "
        f"false alarms per hit {fp / tp:.2f} · misses {fn} · correct rejections {tn}"
    )
    print(
        f"break-even: saying 'wetter' on these calls pays only if a false alarm "
        f"costs less than hits / false alarms = {tp / fp:.2f} of a hit's value"
    )
    base = y_te.mean()
    print(
        f"always-wet's break-even on the same winters, for comparison: "
        f"base rate / (1 − base rate) = {base / (1 - base):.2f}"
    )

    acc = {}
    # the two bars and the physical prior, on the same test members
    for name in RUNG0:
        yb = ladder(RNG_SEED)[name].fit(X_tr, y_tr).predict(X_te)
        acc[name] = accuracy_score(y_te, yb)
        print(
            f"{name:11s} precision {precision_score(y_te, yb, zero_division=0):.3f}  "
            f"recall {recall_score(y_te, yb, zero_division=0):.3f}  "
            f"wet-call share {yb.mean():.3f}  accuracy {acc[name]:.3f}"
        )
    print(
        f"accuracy, final minus always-dry: {acc_final - acc['always-dry']:+.3f} · "
        f"final minus nino-sign: {acc_final - acc['nino-sign']:+.3f}"
    )

    N_BOOT = 1000
    boot_rng = np.random.RandomState(RNG_SEED)
    member_rows = {
        m: np.flatnonzero(test["member_id"].values == m)
        for m in test["member_id"].unique()
    }
    names = list(member_rows)
    point = {
        "precision − base rate": precision_score(
            y_te, yhat_te, zero_division=0
        )
        - y_te.mean(),
        "AP − base rate": average_precision_score(y_te, score_te)
        - y_te.mean(),
    }
    draws = {k: [] for k in point}
    for _ in range(N_BOOT):
        picked = boot_rng.choice(names, size=len(names), replace=True)
        pick = np.concatenate([member_rows[m] for m in picked])
        yb, hb, sb = y_te[pick], yhat_te[pick], score_te[pick]
        draws["precision − base rate"].append(
            precision_score(yb, hb, zero_division=0) - yb.mean()
        )
        draws["AP − base rate"].append(
            average_precision_score(yb, sb) - yb.mean()
        )
    for k, v in draws.items():
        lo, hi = (
            round(float(b), 3) for b in np.percentile(v, [2.5, 97.5])
        )  # as printed
        reading = (
            "above no skill"
            if lo > 0
            else (
                "below no skill"
                if hi < 0
                else "NOT distinguishable from no skill"
            )
        )
        print(
            f"{reg.incumbent_name}: {k:22s} {point[k]:+.3f}   95% interval over test "
            f"members [{lo:+.3f}, {hi:+.3f}]  →  {reading}"
        )

    fig, ax = plt.subplots(figsize=(3.6, 3.4))
    ConfusionMatrixDisplay(cm_te, display_labels=["dry", "wet"]).plot(
        ax=ax, colorbar=False, cmap=GOLD_CMAP
    )
    ax.set_xlim(-0.85, 1.85)
    ax.set_ylim(1.85, -0.85)
    ax.set_title(
        f"Test: {reg.incumbent_name} on {test.member_id.nunique()} held-out members",
        fontsize=9,
    )
    cite(fig, S.access_date)
    plt.tight_layout()
    show_fig(
        fig,
        (
            "The final model's confusion matrix on the held-out test members:"
            " true dry and wet rows against predicted dry and wet columns, "
            "counts in each cell."
        ),
    )

if FINAL:  # reads the test labels, so it sits behind the same switch
    edges = decile_edges(train)  # the TRAINING members' edges, as in chapter 2
    test_bins, chi2_te, p_te, dof_te = wet_by_bin(test, edges)
    print(
        "wet rate by TRAINING decile of the autumn Niño 3.4 anomaly, test "
        f"members (base rate {test['wet'].mean():.3f}):"
    )
    print(test_bins.round(3).to_string())
    print(
        f"chi-square test of wet against decile, test members: "
        f"chi2 = {chi2_te:.1f}, dof = {dof_te}, p = {p_te:.2g}"
    )
    # Post hoc, chosen after the decile table above was seen: the two extreme
    # bins pooled against the eight middle bins, on the test members.
    x = test["nino34_anom"]
    extreme = (x <= edges[1]) | (x > edges[-2])
    parts = {"two extreme bins": extreme, "eight middle bins": ~extreme}
    rates = {}
    for label, mask in parts.items():
        w = test.loc[mask, "wet"]
        rates[label] = (w.mean(), np.sqrt(w.mean() * (1 - w.mean()) / len(w)))
        print(
            f"test members, {label}: wet rate {rates[label][0]:.3f} "
            f"± {rates[label][1]:.3f} (binomial SE, n = {len(w)})"
        )
    (r_e, se_e), (r_m, se_m) = rates.values()
    print(
        f"extremes minus middle: {r_e - r_m:+.3f} ± {np.hypot(se_e, se_m):.3f} "
        "(binomial SEs treat winters as independent: untested for clustering "
        "within members)"
    )

SPINE_WALKED = [
    (
        "framing",
        "chapter 1",
        "a target that does not answer the decision anyone will make",
    ),
    (
        "data",
        "chapters 1 · 2",
        "inputs whose origin and quality were never recorded",
    ),
    (
        "split",
        "chapter 3",
        "training and scoring sets that share ensemble members, so the "
        '"unseen" score is inflated',
    ),
    (
        "train",
        "chapter 4",
        "a fit that cannot be rerun because its settings and seed were "
        "not recorded",
    ),
    (
        "validate",
        "chapter 4",
        "tuning with no held-out score to show whether the model is "
        "memorising",
    ),
    ("refine", "chapters 4 · 5", "stopping at the first configuration tried"),
    (
        "test",
        "chapter 6",
        "a final number that was tuned toward through repeated use of the "
        "same held-out cases",
    ),
    (
        "justify",
        "chapter 6",
        "a model chosen for its sophistication rather than its measured "
        "performance and cost",
    ),
]
fig, ax = plt.subplots(figsize=(12, 3.2))
ax.set(xlim=(0, len(SPINE_WALKED)), ylim=(0.22, 1))
ax.axis("off")
for k, (step, where, failure) in enumerate(SPINE_WALKED):
    ax.add_patch(
        plt.Rectangle(
            (k + 0.06, 0.6),
            0.88,
            0.32,
            facecolor="#E7CD82",
            edgecolor="#8F7022",
            lw=1.2,
        )
    )
    ax.text(
        k + 0.5,
        0.82,
        step,
        ha="center",
        va="center",
        fontsize=10.5,
        fontweight="bold",
    )
    ax.add_patch(
        plt.Rectangle(
            (k + 0.12, 0.62),
            0.76,
            0.12,
            facecolor="#F5F3ED",
            edgecolor="#E4E0D6",
            lw=0.8,
        )
    )
    ax.text(k + 0.5, 0.68, where, ha="center", va="center", fontsize=8)
    ax.text(
        k + 0.5,
        0.54,
        textwrap.fill("guards against " + failure, 21),
        ha="center",
        va="top",
        fontsize=6.6,
        linespacing=1.25,
    )
    if k < len(SPINE_WALKED) - 1:
        ax.annotate(
            "",
            xy=(k + 1.06, 0.76),
            xytext=(k + 0.94, 0.76),
            arrowprops=dict(arrowstyle="->", color="#8F7022", lw=1.3),
        )
ax.set_title(
    "The supervised-learning spine as this cookbook walked it — each step, the "
    "chapter that took it, and the failure it guards against",
    fontsize=9,
)
credit_own(
    fig,
    "a diagram of the analysis spine (shape: NSF Unidata, Apache-2.0); "
    "no CESM2 LENS data drawn",
)
plt.tight_layout()
show_fig(
    fig,
    (
        "A row of eight gold boxes joined by arrows, reading framing, "
        "data, split, train, validate, refine, test, justify; under each "
        "step's name is the chapter that took it (chapter 1; chapters 1 "
        "and 2; chapter 3; chapter 4 for train and validate; chapters 4 "
        "and 5; chapter 6 for test and justify) and below the box the "
        "failure the step guards against."
    ),
)