lib/test_stats.py

#!/usr/bin/env python3
"""Validation: lab.lib must reproduce the published stat-hardening numbers.

Ground truth = results/verdict.json + mcnemar_matrix.csv from the 30-run
study (sandbox/stage2/stat_hardening_g). If these tests pass, the library is
a faithful generalization, not a rewrite.
"""
import csv
import json
import sys
from pathlib import Path

import numpy as np

sys.path.insert(0, str(Path.home() / "jupyterlab"))
from lab.lib.stats import (benjamini_hochberg, holm, lesion_cluster_bootstrap,
                           mean_confidence_vote, mcnemar_exact,
                           majority_vote, mean_std, pooled_mcnemar_exact,
                           wilson)
from lab.lib.metrics import (compute_metrics, lesion_frame_metrics,
                             lesion_level_metrics)

SH = Path.home() / "jupyterlab/sandbox/stage2/stat_hardening_g"
RUNS, RESULTS = SH / "runs", SH / "results"
ok = fail = 0


def check(name, cond, detail=""):
    global ok, fail
    status = "PASS" if cond else "FAIL"
    if cond:
        ok += 1
    else:
        fail += 1
    print(f"[{status}] {name}" + (f"  ({detail})" if detail else ""))


# ── unit sanity ───────────────────────────────────────────────────────────────
mu, sd, n = mean_std([0.1, 0.2, 0.3])
check("mean_std", abs(mu - 0.2) < 1e-12 and abs(sd - 0.1) < 1e-12 and n == 3)
lo, hi = wilson(3, 10)
check("wilson", lo < 0.3 < hi and lo > 0.08, f"{lo:.3f}..{hi:.3f}")
adj = holm([0.01, 0.04, 0.03])
# textbook Holm on these inputs: sorted 0.01,0.03,0.04 -> 0.03, 0.06, 0.06
check("holm textbook values", np.allclose(adj, [0.03, 0.06, 0.06]), f"{np.round(adj,4)}")
q = benjamini_hochberg([0.01, 0.02, 0.03])
# step-up monotone from largest: all saturate at p_max*m = 0.09
check("BH textbook values", np.allclose(q, [0.09, 0.09, 0.09]), f"{np.round(q,4)}")

rng = np.random.default_rng(0)
t = rng.integers(0, 3, 500)
pa = t.copy(); pa[rng.random(500) < .1] = 0          # A loses some frames
pb = t.copy(); pb[rng.random(500) < .05] = 1         # B loses fewer
ea, eb = pa != t, pb != t
r = mcnemar_exact(pa, pb, t)
b_manual = int((eb & ~ea).sum())                     # A right, B wrong
c_manual = int((ea & ~eb).sum())
check("mcnemar matches manual counts",
      r["b"] == c_manual and r["c"] == b_manual and 0 < r["p"] <= 1,
      f"b={r['b']} c={r['c']} p={r['p']:.3f}")

# majority_vote vs mean_confidence on a constructed lesion:
# one CONFIDENT NML frame + two WEAK NROI frames.
#   hard vote: argmax labels {NML, NROI, NROI} → NROI wins 2-1
#   soft mean: NML prob .98 dominates the average → NML
probs = np.array([[.01, .98, .01],   # → NML, confident
                  [.30, .30, .40],   # → NROI, weak
                  [.28, .30, .42]])  # → NROI, weak
lids = ["L1", "L1", "L1"]
hard = majority_vote(probs.argmax(1), lids)
soft = mean_confidence_vote(probs, lids)
check("soft vote rescues minority", hard[0] == 2 and soft[0] == 1,
      f"hard={hard} soft={soft} (mean row={probs.mean(0).round(2)})")

# ── reproduce study numbers ──────────────────────────────────────────────────
def preds(cfg, prep, seed):
    return np.load(RUNS / f"{cfg}_{prep}_seed{seed}/preds_test.npz",
                   allow_pickle=True)


z = preds("cascade_majvote", "norm", 44)
mv = majority_vote(z["pred_novote"], z["lesion_ids"])
agree = float((mv[np.array([l is not None for l in z["lesion_ids"]])] ==
               z["pred"][np.array([l is not None for l in z["lesion_ids"]])]).mean())
check("majority_vote reproduces harness voted preds", agree > 0.99, f"agreement={agree:.4f}")

lfm = lesion_frame_metrics(z["pred"], z["true"], z["lesion_ids"])
ref = json.loads((RUNS / "cascade_majvote_norm_seed44/metrics.json").read_text())["per_lesion_unified_test"]
check("lesion_frame_metrics == harness per_lesion_unified_test",
      abs(lfm["mcc"] - ref["mcc"]) < 1e-9 and
      abs(lfm["per_class_recall"]["non-malignant lesion"] -
          ref["per_class_recall"]["non-malignant lesion"]) < 1e-9,
      f"mcc {lfm['mcc']:.4f} vs {ref['mcc']:.4f}")

llm = lesion_level_metrics(z["pred"], z["true"], z["lesion_ids"])
check("lesion_level_metrics (one sample/lesion) computed, differs as expected",
      0 <= llm["mcc"] <= 1 and llm["accuracy"] != lfm["accuracy"],
      f"per-lesion mcc={llm['mcc']:.4f} n_lesion_samples≈{round(llm['confusion_matrix'][0][0]+llm['confusion_matrix'][0][1]+llm['confusion_matrix'][0][2] + sum(llm['confusion_matrix'][i][0] for i in range(3)))}")

fm = compute_metrics(z["pred"], z["true"], sources=z["sources"])
ref_f = json.loads((RUNS / "cascade_majvote_norm_seed44/metrics.json").read_text())["combined"]["test"]
check("compute_metrics == combined.test",
      abs(fm["mcc"] - ref_f["mcc"]) < 1e-9 and fm["confusion_matrix"] == ref_f["confusion_matrix"],
      f"mcc {fm['mcc']:.4f}")

rows = list(csv.DictReader(open(RESULTS / "mcnemar_matrix.csv")))
sel = [r_ for r_ in rows if r_["config_a"] == "cascade" and r_["config_b"] == "baseline"
       and r_["prep"] == "norm"]
recomputed = []
for r_ in sel:
    za = preds(r_["config_a"], "norm", int(r_["seed"]))
    zb = preds(r_["config_b"], "norm", int(r_["seed"]))
    res = mcnemar_exact(za["pred"], zb["pred"], za["true"])
    recomputed.append((res, r_))
match = all(abs(res["p"] - float(r_["p_raw"])) < 1e-4 and res["b"] == int(r_["b_disc"])
            for res, r_ in recomputed)
check("mcnemar reproduces matrix.csv (all seeds)", match,
      f"{len(recomputed)} pairs")
pool = pooled_mcnemar_exact([res for res, _ in recomputed])
check("pooled b/c sums match manual audit (296/354 was norm-vs-raw; here cascade-vs-baseline)",
      pool["b"] == sum(int(r_["b_disc"]) for r_ in sel), f"b={pool['b']} p={pool['p']:.5f}")

bs = lesion_cluster_bootstrap(preds("cascade", "norm", 42)["true"],
                              preds("cascade", "norm", 42)["pred"],
                              preds("cascade", "norm", 42)["lesion_ids"], "nml_recall",
                              n_boot=2000)
cis = json.loads((RESULTS / "bootstrap_cis.json").read_text())["cascade_norm"]["nml_recall_per_seed"][0]
close = abs(bs["point"] - cis["point"]) < 1e-9 and bs["n_lesions"] == cis["n_lesions"]
check("cluster bootstrap point+n_lesions identical", close,
      f"point {bs['point']:.4f}, L={bs['n_lesions']}")

print(f"\n{ok} passed, {fail} failed")
sys.exit(1 if fail else 0)