An out-of-period test, under a protocol fixed before the run

import contextlib
import hashlib
import io
import math

import matplotlib.pyplot as plt
import numpy as np
from lensml import (
    FEATURES,
    FUTURE_ACCESS_DATE,
    FUTURE_YEARS,
    RNG_SEED,
    TRIALS,
    apply_transforms,
    at_base_rate,
    cite,
    decile_edges,
    fit_transforms,
    ladder,
    load_future_cache,
    score_frozen,
    setup,
    show,
    show_fig,
    wet_by_bin,
)

show(score_frozen)

PROTOCOL = 'data/ssp370-test-protocol.md'
# The hash of the file as committed with the test (851d900, 2026-10-29),
# recorded here on 2026-10-30: a later edit shows here as 'differs'.
PROTOCOL_SHA256 = (
    '11d291795a00bb6d10c818c27a40f9a2a45dff1d8641dfacf931a044489d9628'
)
with open(PROTOCOL, 'rb') as f:
    sha = hashlib.sha256(f.read()).hexdigest()
print(f'{PROTOCOL} sha256 {sha} (first 16: {sha[:16]})')
print(
    'protocol file: '
    + ('matches' if sha == PROTOCOL_SHA256 else 'DIFFERS from')
    + ' the hash committed with the test'
)

S = (
    setup()
)  # the historical table: training members, test members, chapter 6's model inputs
fut, fut_header = load_future_cache()
print(fut_header)
fut['split'] = np.where(fut['member_id'].isin(S.test_members), 'test', 'train')
for split, rows in fut.groupby('split'):
    print(
        f'{split}: {rows.member_id.nunique()} members · {len(rows)} future winters'
    )

fut_train = fut[fut['split'].eq('train')]
TF_FUTURE = fit_transforms(fut_train, years=FUTURE_YEARS)
fut = apply_transforms(fut, TF_FUTURE)
fut_train = fut[fut['split'].eq('train')]
print(
    f"future wet threshold (pooled training members) = {TF_FUTURE['wet_threshold']:+.4f} mm/day; "
    f"training-member future base rate {fut_train['wet'].mean():.3f}"
)
decade = (fut_train['winter'] // 10) * 10
for anom in ['nino34_anom', 'prect_anom']:
    g = fut_train.groupby(decade)[anom]
    z = g.mean() / (fut_train[anom].std() / np.sqrt(g.size()))
    print(
        f"{anom:12s} training members, decade-mean anomaly in standard errors: worst {z[z.abs().idxmax()]:+.2f} "
        f"(the {z.abs().idxmax()}s) · tolerance ±3 → {'PASS' if (z.abs() <= 3).all() else 'FAIL'}"
    )

# The shift in the input, training members only: no test-member future statistic is computed here.
sd_hist = S.train['nino34_anom'].std()
sd_fut = fut_train['nino34_anom'].std()
print(
    f'SON Niño 3.4 anomaly SD, training members: historical {sd_hist:.3f} K · '
    f'future {sd_fut:.3f} K · ratio {sd_fut / sd_hist:.3f}'
)

X_tr, y_tr = S.train[FEATURES].values, S.train['wet'].values
name5, make5, _ = TRIALS[4]
assert name5 == 'boosting@base-rate'
models = {
    'boosting@base-rate': make5(RNG_SEED).fit(X_tr, y_tr),
    'always-wet': ladder(RNG_SEED)['always-wet'].fit(X_tr, y_tr),
    'always-dry': ladder(RNG_SEED)['always-dry'].fit(X_tr, y_tr),
    'nino-sign': ladder(RNG_SEED)['nino-sign'].fit(X_tr, y_tr),
    'logistic@base-rate': at_base_rate('logistic')(RNG_SEED).fit(X_tr, y_tr),
}
edges = decile_edges(S.train)  # chapter 2's training edges
print('frozen:', ', '.join(models))

SCORED_RUNS = globals().get('SCORED_RUNS', 0) + 1
print(
    f'scored run #{SCORED_RUNS}'
    + ('' if SCORED_RUNS == 1 else '   — a re-run in this kernel')
)
fut_test = fut[fut['split'].eq('test')].reset_index(drop=True)
buf = io.StringIO()
with contextlib.redirect_stdout(buf):
    results = score_frozen(models, fut_test, edges)
pt, lo, hi = results['boosting@base-rate']['precision − base rate']
verdict = 'EDGE REPLICATES' if lo > 0 else 'DOES NOT REPLICATE'
print(
    f"PRIMARY OUTCOME — boosting@base-rate, precision − base rate on the test members' future winters: {pt:+.3f}, "
    f'95% interval over test members [{lo:+.3f}, {hi:+.3f}] → {verdict} '
    '(pre-registered rule: replicates only if the lower bound is above 0)'
)
print('\nSECONDARY OUTCOMES, in the protocol\'s order:')
print(buf.getvalue())

# The first scored run, as the author recorded it on 2026-10-29 (first
# committed 2026-10-30). A library upgrade that changes the scored result
# shows here as 'differs'; report a changed result as a re-run.
FIRST_RUN = (0.011, -0.009, 0.031, 'DOES NOT REPLICATE')
now = (round(pt, 3), round(lo, 3), round(hi, 3), verdict)
print(
    'scored result: '
    + ('matches' if now == FIRST_RUN else f'DIFFERS from {FIRST_RUN}:')
    + (' the first scored run' if now == FIRST_RUN else f' now {now}')
)

# Chapter 6's estimate, recomputed from the historical test members (already scored there), against this interval.
yhat_hist = models['boosting@base-rate'].predict(S.test[FEATURES].values)
w_hist = S.test['wet'].values
est6 = w_hist[yhat_hist == 1].mean() - w_hist.mean()
print(
    f"chapter 6's estimate (historical test members) {est6:+.3f}; inside this "
    f'interval [{lo:+.3f}, {hi:+.3f}]: {lo <= round(est6, 3) <= hi}'
)
# Outside the protocol: how often a test this size would confirm an edge of
# chapter 6's size (a normal approximation to the interval over members).
se7 = (hi - lo) / (2 * 1.96)
power7 = 0.5 * (1 - math.erf((1.96 - est6 / se7) / math.sqrt(2)))
print(
    f'outside the protocol: with this interval\'s width (SE ≈ {se7:.4f}), '
    f"a true edge of chapter 6's size clears a lower bound of 0 about "
    f'{power7:.0%} of the time'
)

names = [
    n for n in results if n != 'always-dry'
]  # undefined precision, not drawn
fig, ax = plt.subplots(figsize=(7, 3.2))
for i, name in enumerate(names):
    p, a, b = results[name]['precision − base rate']
    ax.plot([a, b], [i, i], color='#C9A23A', lw=9, solid_capstyle='butt')
    ax.plot(
        p, i, 'o', color='#8F7022', markeredgecolor='#13161A', markersize=9
    )
ax.axvline(0, color='#13161A', ls='--', lw=1)
ax.set_yticks(range(len(names)), names)
ax.invert_yaxis()
ax.set_xlabel(
    'precision − base rate, test members\' SSP3-7.0 winters (95% interval)'
)
ax.set_title(
    'The scored run: precision against no skill, frozen models\n'
    '(always-dry calls nothing wet: its precision is undefined, so it is not drawn)',
    fontsize=9,
)
cite(fig, FUTURE_ACCESS_DATE)  # the ssp370 table's access date
plt.tight_layout()
show_fig(
    fig,
    (
        "Horizontal interval plot: for each frozen model — boosting at "
        "the base-rate threshold, always-wet, the Niño-sign rule and "
        "logistic at the base-rate threshold — a dot for precision minus "
        "the base rate on the test members' SSP3-7.0 winters and a gold "
        "bar for its 95% interval over test members, against a dashed "
        "line at zero, no skill. Always-dry is not drawn: it calls "
        "nothing wet, so its precision is undefined."
    ),
)

# Sensitivity, outside the protocol: the decile chi-square on the test
# members' future winters under each reading of "the training edges".
for label, e in [
    ("chapter 2's historical training edges (the run)", edges),
    ("the training members' future-winter edges", decile_edges(fut_train)),
]:
    _, chi2, p, dof = wet_by_bin(fut_test, e)
    print(f'{label}: chi2 {chi2:.1f}, dof {dof}, p {p:.2f}')