# Setup: dependencies, the shared helpers, non-interactive plotting.
import importlib.util
import subprocess
import sys
import time
import warnings
from pathlib import Path

_needed = ("mne", "scipy", "matplotlib", "pandas", "pooch", "sklearn")
_missing = [p for p in _needed if importlib.util.find_spec(p) is None]
if _missing:
    _req = next((d / "requirements.txt" for d in (Path.cwd(), *Path.cwd().parents)
                 if (d / "requirements.txt").exists()), None)
    _cmd = [sys.executable, "-m", "pip", "install", "-q"]
    _cmd += ["-r", str(_req)] if _req else ["mne==1.10.2", "pooch>=1.8", "scikit-learn"]
    subprocess.check_call(_cmd)

_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
                if (d / "_shared" / "helpers_l6.py").exists()), None)
if _shared is None:
    raise FileNotFoundError("start the kernel in notebooks/L6/ (or notebooks/) so that _shared/helpers_l6.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l3 as L3
import helpers_l6 as L6

import matplotlib.pyplot as plt
import numpy as np
import mne
from scipy import stats
from sklearn.model_selection import StratifiedKFold, cross_val_score
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import roc_auc_score

mne.set_log_level("WARNING")
plt.rcParams["figure.dpi"] = 72

# Iteration counts.  Time-resolved decoding is linear in the number of time points and the number of
# folds; the group permutation test is linear in its permutation count.  Both are budgeted here.
FULL_RUN = False
N_SUBJECTS = 10 if not FULL_RUN else 20           # ds-erpcore subjects for the time-resolved analysis
DECIM = 2 if not FULL_RUN else 1                  # 256 Hz -> 128 Hz for the sliding estimator
DECIM_GEN = 4 if not FULL_RUN else 2              # 256 Hz -> 64 Hz for the generalization matrix
N_PERM = 10000 if not FULL_RUN else 50000         # group-level cluster permutation
N_PERM_WITHIN = 200 if not FULL_RUN else 2000     # label permutations inside one subject (section 3)
CSP_CANDIDATES = [f"S{i:03d}" for i in range(1, 13)]
FOLDS = 5
SEED = L6.SEED
ALPHA = L6.ALPHA

print(f"MNE {mne.__version__}; helpers_l6 imported from notebooks/_shared")
print(f"FULL_RUN = {FULL_RUN}: {N_SUBJECTS} ERP CORE subjects, decimation {DECIM} for the sliding "
      f"estimator and {DECIM_GEN} for the generalization matrix, {N_PERM} group permutations, "
      f"{N_PERM_WITHIN} within-subject label permutations, {len(CSP_CANDIDATES)} ds-eegbci candidates")
print(f"ERP CORE cache: {L3.erpcore_root().name}/; each subject's EEGLAB pair is deleted as soon as its "
      f"epochs have been measured")

SUBJECTS = list(L6.SUBSET_DEFAULT)[:N_SUBJECTS]
L6.disk_report("before any download")
t0 = time.time()
data = {}
for s in SUBJECTS:
    sid = L3.erpcore_subject_id(s)
    t1 = time.time()
    try:
        cms, cond, tt, rt, acc, code = L6._subject_epochs(s, L6.HIGHPASS[L6.PRESPECIFIED["highpass"]],
                                                          verbose=False)
    finally:
        L3.delete_erpcore_subject(s, "P3", keep_small=True, verbose=False)
    ref = L6.rereference(cms, L6.PRESPECIFIED["reference"])
    b0, b1 = L6.BASELINE[L6.PRESPECIFIED["baseline"]]
    ref = ref - ref[:, :, (tt >= b0) & (tt <= b1)].mean(axis=2, keepdims=True)
    keep = np.ptp(ref, axis=2).max(axis=1) < L6.REJECT_UV[L6.PRESPECIFIED["rejection"]]
    y = (cond == "target").astype(int)[keep]
    if min(int(y.sum()), int((1 - y).sum())) < int(L6.FIXED_CHOICES["min_trials_per_condition"]):
        print(f"  {sid}: EXCLUDED, {int(y.sum())} target / {int((1 - y).sum())} standard trials survive")
        continue
    data[sid] = {"X": ref[keep][:, :, ::DECIM].astype(np.float64), "y": y}
    print(f"  {sid}: {data[sid]['X'].shape[0]} trials ({int(y.sum())} target), "
          f"{data[sid]['X'].shape[2]} time points in {time.time() - t1:.1f} s", flush=True)
times = tt[::DECIM]
print(f"{len(data)} subjects in {time.time() - t0:.0f} s; "
      f"sampling {1 / np.diff(times).mean():.0f} Hz after decimation by {DECIM}")
L6.disk_report("after the download loop")

sids = sorted(data)
t0 = time.time()
auc = {}
for s in sids:
    sc, desc = L6.sliding_auc(data[s]["X"], data[s]["y"], folds=FOLDS, seed=SEED)
    auc[s] = sc.mean(0)                                # mean over folds, one curve per subject
print(f"time-resolved decoding of {len(sids)} subjects x {len(times)} time points x {FOLDS} folds "
      f"in {time.time() - t0:.0f} s")
print(f"estimator: {desc}")
A = np.stack([auc[s] for s in sids])                   # subjects x times
peak_i = int(np.argmax(A.mean(0)))
print(f"\ngroup mean AUC: peak {A.mean(0)[peak_i]:.4f} at {times[peak_i] * 1000:.0f} ms")
print(f"{'subject':9s} {'peak AUC':>9s} {'at (ms)':>8s} {'AUC at the group peak':>22s} {'pre-stimulus max':>18s}")
pre = times < 0
for s in sids:
    i = int(np.argmax(auc[s]))
    print(f"{s:9s} {auc[s][i]:9.4f} {times[i] * 1000:8.0f} {auc[s][peak_i]:22.4f} "
          f"{auc[s][pre].max():18.4f}")

t0 = time.time()
thr = float(stats.t.ppf(1 - ALPHA, len(sids) - 1))          # one-sided cluster-forming threshold
t_obs, clusters, cluster_p, H0 = mne.stats.permutation_cluster_1samp_test(
    A - 0.5, threshold=thr, n_permutations=N_PERM, tail=1, out_type="indices", seed=SEED, verbose=False)
mask = np.zeros(len(times), bool)
print(f"mne.stats.permutation_cluster_1samp_test(AUC - 0.5, threshold={thr:.4f}, "
      f"n_permutations={N_PERM}, tail=1, seed={SEED})")
print(f"permutations used {len(H0)} (2^{len(sids)} = {2 ** len(sids):,} sign flips exist); smallest "
      f"attainable p = {1 / len(H0):.5f}; computed in {time.time() - t0:.1f} s\n")
print(f"{'#':>2s} {'range (ms)':>20s} {'samples':>8s} {'t-sum':>9s} {'p':>8s} {'peak AUC in it':>15s}")
sig_clusters = []
for rank, i in enumerate(np.argsort(cluster_p)):
    idx = np.asarray(clusters[i][0])
    ok = cluster_p[i] <= ALPHA
    if ok:
        mask[idx] = True
        sig_clusters.append({"t_start_ms": float(times[idx[0]] * 1000),
                             "t_end_ms": float(times[idx[-1]] * 1000),
                             "p": float(cluster_p[i]), "n": int(idx.size)})
    if ok or rank < 3:
        print(f"{rank:2d} {times[idx[0]] * 1000:8.0f} to {times[idx[-1]] * 1000:8.0f} {idx.size:8d} "
              f"{t_obs[idx].sum():+9.1f} {cluster_p[i]:8.4f} {A.mean(0)[idx].max():15.4f}"
              + ("   <- significant" if ok else ""))
print(f"\n{len(sig_clusters)} cluster(s) survive at alpha = {ALPHA}")
print(f"ex-6-5 ANSWER: the group peak is at {times[peak_i] * 1000:.0f} ms with AUC "
      f"{A.mean(0)[peak_i]:.4f}, and it "
      f"{'IS' if mask[peak_i] else 'is NOT'} inside a significant cluster.")

fig, ax = L6.plot_decoding_curve(
    times, A, chance=0.5, mask=mask,
    title=f"ERP CORE P3, target vs standard: time-resolved decoding, {len(sids)} subjects (ROC AUC)")
for s in sids:
    ax.plot(times * 1000, auc[s], lw=0.6, alpha=0.45, color="0.5", zorder=0)
ax.plot([], [], lw=0.6, color="0.5", label="individual subjects")
ax.legend(fontsize=8, loc="upper left")
plt.show()   # render the static figure(s) of this cell inline

n_tr = int(np.median([len(data[s]["y"]) for s in sids]))
n_tg = int(np.median([int(data[s]["y"].sum()) for s in sids]))
maj = 1 - n_tg / n_tr
lo, hi = L6.chance_band(n_tr, 0.5)
print(f"On a median subject: {n_tr} trials, {n_tg} target and {n_tr - n_tg} standard.")
print(f"   1. 'accuracy of the majority class'  {maj:.3f}  <- what a classifier that always answers "
      f"'standard' scores.  This is why accuracy is the wrong metric for a 1:4 design.")
print(f"   2. 'AUC of a coin flip'              0.500  <- the right centre for AUC, which is invariant to "
      f"the class ratio.")
print(f"   3. the SPREAD of chance              binomial 95 % interval on {n_tr} balanced trials runs "
      f"{lo:.3f} to {hi:.3f}; a single subject inside that band is not evidence of anything.")
print()
print(f"But the analytic band above assumes independent trials and a balanced split, and cross-validated AUC "
      f"on {FOLDS} folds is neither.  The honest yardstick is the EMPIRICAL null: shuffle the labels inside "
      f"the subject and rerun the same cross-validation.")

t0 = time.time()
s_demo = sids[int(np.argsort([auc[s].max() for s in sids])[len(sids) // 2])]     # the median subject
Xd, yd = data[s_demo]["X"], data[s_demo]["y"]
i_win = (times >= 0.30) & (times <= 0.60)
Xw = Xd[:, :, i_win].reshape(len(yd), -1)
clf = make_pipeline(StandardScaler(), LogisticRegression(class_weight="balanced", max_iter=2000,
                                                         solver="liblinear"))
cv = StratifiedKFold(n_splits=FOLDS, shuffle=True, random_state=SEED)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    obs = float(cross_val_score(clf, Xw, yd, cv=cv, scoring="roc_auc").mean())
    rng = np.random.default_rng(SEED)
    null = np.array([cross_val_score(clf, Xw, rng.permutation(yd), cv=cv, scoring="roc_auc").mean()
                     for _ in range(N_PERM_WITHIN)])
p_emp = (1 + int((null >= obs).sum())) / (1 + len(null))
print(f"\n{s_demo}, 300-600 ms window, {Xw.shape[1]} features ({Xd.shape[1]} channels x "
      f"{int(i_win.sum())} samples):")
print(f"   observed cross-validated AUC {obs:.4f}")
print(f"   empirical null over {N_PERM_WITHIN} label shuffles: mean {null.mean():.4f}, SD {null.std(ddof=1):.4f}, "
      f"2.5-97.5 % [{np.percentile(null, 2.5):.4f}, {np.percentile(null, 97.5):.4f}]")
print(f"   p = (1 + #{{null >= observed}}) / (1 + {N_PERM_WITHIN}) = {p_emp:.4f}")
print(f"   the null's centre is {null.mean():.4f}, not exactly 0.500 -- with {len(yd)} trials and {FOLDS} "
      f"folds the estimator has a small bias, which is precisely why the permutation null is worth having")
print(f"   ({time.time() - t0:.0f} s; a full run would use {2000} shuffles, "
      f"about {(time.time() - t0) * 2000 / N_PERM_WITHIN:.0f} s)")

fig, ax = plt.subplots(figsize=(8.6, 4.0))
ax.hist(null, bins=40, color="tab:blue", alpha=0.8, label=f"{N_PERM_WITHIN} label shuffles")
ax.axvline(0.5, color="gray", lw=1.2, ls="--", label="0.5 (nominal chance)")
ax.axvline(null.mean(), color="k", lw=1.2, ls=":", label=f"null mean {null.mean():.3f}")
ax.axvline(obs, color="tab:orange", lw=2.2, label=f"observed {obs:.3f} (p = {p_emp:.4f})")
ax.set(xlabel="Cross-validated ROC AUC (dimensionless)", ylabel="Shuffles",
       title=f"The empirical null for one subject ({s_demo}), 300–600 ms window")
ax.grid(alpha=0.3)
ax.legend(fontsize=8)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

t0 = time.time()
Xg = data[s_demo]["X"][:, :, ::(DECIM_GEN // DECIM)]
tg = times[::(DECIM_GEN // DECIM)]
gen = L6.generalization_auc(Xg, data[s_demo]["y"], folds=FOLDS, seed=SEED).mean(0)
print(f"{s_demo}: {len(tg)} x {len(tg)} = {len(tg) ** 2:,} classifiers x {FOLDS} folds in "
      f"{time.time() - t0:.0f} s ({1 / np.diff(tg).mean():.0f} Hz after decimation by {DECIM_GEN})")
print(f"   a full run at {1 / np.diff(times).mean():.0f} Hz would be "
      f"{(len(times) / len(tg)) ** 2:.0f}x this, about {(time.time() - t0) * (len(times) / len(tg)) ** 2 / 60:.0f} min")
d = np.diag(gen)
print(f"   diagonal peak {d.max():.4f} at {tg[int(np.argmax(d))] * 1000:.0f} ms")
print(f"   matrix maximum {gen.max():.4f} at train {tg[np.unravel_index(gen.argmax(), gen.shape)[0]] * 1000:.0f} ms, "
      f"test {tg[np.unravel_index(gen.argmax(), gen.shape)[1]] * 1000:.0f} ms")
post = tg >= 0
off = gen[np.ix_(post, post)].copy()
np.fill_diagonal(off, np.nan)
print(f"   mean off-diagonal AUC after 0 ms: {np.nanmean(off):.4f} against a diagonal mean of "
      f"{d[post].mean():.4f}")
_ratio = (np.nanmean(off) - 0.5) / (d[post].mean() - 0.5)
print(f"   ratio {np.nanmean(off) - 0.5:+.4f} / {d[post].mean() - 0.5:+.4f} = {_ratio:+.3f}")
print(f"   Read it as: 1 would mean the pattern persists unchanged (a square block), 0 that it does not "
      f"generalise off the diagonal at all (a narrow band), and a NEGATIVE value that training at one "
      f"moment predicts the opposite class at another -- which happens when the ERP changes sign between "
      f"the early and late parts of the epoch, and is a finding about the waveform rather than a defect.")

fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.4), gridspec_kw={"width_ratios": [1.15, 1]})
v = float(np.abs(gen - 0.5).max())
im = axes[0].imshow(gen, origin="lower", cmap="RdBu_r", vmin=0.5 - v, vmax=0.5 + v,
                    extent=[tg[0] * 1000, tg[-1] * 1000, tg[0] * 1000, tg[-1] * 1000])
axes[0].axhline(0, color="k", lw=0.6)
axes[0].axvline(0, color="k", lw=0.6)
axes[0].plot([tg[0] * 1000, tg[-1] * 1000], [tg[0] * 1000, tg[-1] * 1000], "k:", lw=0.7)
axes[0].set(xlabel="Test time (ms)", ylabel="Train time (ms)",
            title=f"Temporal generalization, {s_demo} (ROC AUC)")
cb = fig.colorbar(im, ax=axes[0])
cb.set_label("ROC AUC (dimensionless, 0.5 = chance)")

axes[1].plot(tg * 1000, d, color="tab:blue", lw=1.8, label="diagonal (train = test)")
axes[1].plot(tg * 1000, gen[int(np.argmax(d))], color="tab:orange", lw=1.5,
             label=f"trained at {tg[int(np.argmax(d))] * 1000:.0f} ms, tested everywhere")
axes[1].axhline(0.5, color="gray", ls="--", lw=0.9, label="chance")
axes[1].axvline(0, color="gray", lw=0.6)
axes[1].set(xlabel="Test time (ms)", ylabel="ROC AUC (dimensionless)",
            title="A slice through the matrix")
axes[1].grid(alpha=0.3)
axes[1].legend(fontsize=8)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

print(f"Pre-stimulus decoding ({int(pre.sum())} time points before 0 ms):")
print(f"   group mean AUC over the pre-stimulus interval : {A[:, pre].mean():.4f}")
print(f"   group mean AUC at its pre-stimulus maximum    : {A.mean(0)[pre].max():.4f} at "
      f"{times[pre][int(np.argmax(A.mean(0)[pre]))] * 1000:.0f} ms")
print(f"   subjects whose pre-stimulus maximum exceeds 0.55: "
      f"{int(sum(auc[s][pre].max() > 0.55 for s in sids))} of {len(sids)}")
pre_sig = mask[pre].sum()
print(f"   pre-stimulus points inside a significant cluster: {int(pre_sig)}")
t_pre, p_pre = stats.ttest_1samp(A[:, pre].mean(axis=1), 0.5)
print(f"   t-test of the per-subject pre-stimulus mean against 0.5: t({len(sids) - 1}) = {t_pre:.3f}, "
      f"p = {p_pre:.4f}")
print()
print("What to do about it, in order of preference:")
print("   1. report it.  A decoding figure without its pre-stimulus interval hides its own control condition.")
print("   2. high-pass less aggressively, or use a causal filter when the timing of the effect is the claim.")
print("   3. baseline-correct with a window that is not inside the epoch being decoded, or not at all -- a")
print("      classifier with a per-trial scaler does not need the baseline the way an ERP average does.")
print("   TODO(confirm): how much of the pre-stimulus signal above is filter smearing and how much is")
print("   baseline coupling is not separated here.  Doing so needs the same epochs built without each step,")
print("   which is a comparison this notebook does not run.")

K_CHANNELS = 5
t0 = time.time()
rows = []
for s in sids:
    X_, y_ = data[s]["X"], data[s]["y"]
    Xw_ = X_[:, :, i_win]                                     # 300-600 ms window, channels x samples
    feat = Xw_.mean(axis=2)                                   # mean amplitude per channel, one per trial

    # LEAKED: pick the K channels with the largest |t| using every trial, then cross-validate.
    t_all = np.abs(stats.ttest_ind(feat[y_ == 1], feat[y_ == 0], axis=0).statistic)
    sel_all = np.argsort(t_all)[-K_CHANNELS:]
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        leaked = float(cross_val_score(clf, feat[:, sel_all], y_, cv=cv, scoring="roc_auc").mean())

        # HONEST: the same selection, refitted inside every training fold.
        scores = []
        for tr, te in cv.split(feat, y_):
            t_tr = np.abs(stats.ttest_ind(feat[tr][y_[tr] == 1], feat[tr][y_[tr] == 0], axis=0).statistic)
            sel = np.argsort(t_tr)[-K_CHANNELS:]
            m = clf.fit(feat[tr][:, sel], y_[tr])
            scores.append(float(roc_auc_score(y_[te], m.predict_proba(feat[te][:, sel])[:, 1])))
        honest = float(np.mean(scores))
    rows.append({"subject": s, "honest": honest, "leaked": leaked, "gap": leaked - honest,
                 "selected_all": [L6.ERPCORE_30[i] for i in sorted(sel_all)]})
print(f"channel selection inside vs outside the fold, {len(rows)} subjects, top {K_CHANNELS} of "
      f"{len(L6.ERPCORE_30)} channels, {time.time() - t0:.0f} s\n")
print(f"{'subject':9s} {'honest AUC':>11s} {'leaked AUC':>11s} {'gap':>8s}  channels chosen on all the data")
for r in rows:
    print(f"{r['subject']:9s} {r['honest']:11.4f} {r['leaked']:11.4f} {r['gap']:+8.4f}  "
          f"{', '.join(r['selected_all'])}")
gaps = np.array([r["gap"] for r in rows])
print(f"\nmean honest {np.mean([r['honest'] for r in rows]):.4f}, mean leaked "
      f"{np.mean([r['leaked'] for r in rows]):.4f}, mean gap {gaps.mean():+.4f} "
      f"(range {gaps.min():+.4f} to {gaps.max():+.4f})")
print(f"subjects where the leak HELPED: {int((gaps > 0).sum())} of {len(gaps)}")
print(f"\nThe gap is modest here because {K_CHANNELS} of {len(L6.ERPCORE_30)} channels is a small selection "
      f"and the P3 effect is large and spatially broad, so the leaked and honest selections usually agree.  "
      f"It grows with the size of the search (selecting 5 of 500 voxels, or a time window as well as a "
      f"channel) and with the weakness of the effect -- which is exactly the regime where it does the most "
      f"damage and is the hardest to notice.")

print("the L5.7 pipeline, restated (helpers_l6.CSP_SPEC):")
for k, v in L6.CSP_SPEC.items():
    print(f"   {k:14s} : {v}")

SELECTED = "S002"
DELETE_EEGBCI = True      # spec section 11: delete the subset this notebook downloads.
# The ds-eegbci cache is SHARED -- Level 0 and Level 5 notebooks read the same EDF files -- so the
# listing below records what was on disk before anything was fetched, and only files this run
# actually downloads are deleted afterwards.  Deleting a file another lesson had already cached
# would cost that lesson a re-download for no benefit.
PRE_EXISTING = sorted({str(p) for sid in ([SELECTED] + CSP_CANDIDATES)
                       for p in L6.eegbci_run_files(sid, L6.CSP_SPEC["runs"])})
print(f"ds-eegbci EDFs already cached before this notebook runs: {len(PRE_EXISTING)} of "
      f"{len(set([SELECTED] + CSP_CANDIDATES)) * len(L6.CSP_SPEC['runs'])}; only files fetched by this "
      f"run will be deleted")

t0 = time.time()
Xc, yc, ch = L6.csp_epochs(SELECTED)
lk = L6.csp_leakage(Xc, yc, n_components=int(L6.CSP_SPEC["n_components"]), folds=FOLDS, seed=SEED)
print(f"\n{SELECTED}: {len(yc)} epochs ({int((yc == 0).sum())} left, {int((yc == 1).sum())} right), "
      f"{len(ch)} channels, {time.time() - t0:.0f} s")
print(f"   HONEST (CSP refitted in every training fold) : {lk['honest']:.4f}  "
      f"folds {', '.join(f'{v:.4f}' for v in lk['honest_scores'])}")
print(f"   LEAKED (CSP fitted once on every epoch)      : {lk['leaked']:.4f}  "
      f"folds {', '.join(f'{v:.4f}' for v in lk['leaked_scores'])}")
print(f"   THE COST OF THE LEAK                         : {lk['inflation']:+.4f} "
      f"({100 * lk['inflation']:+.1f} percentage points)")
clo, chi = L6.chance_band(len(yc), 0.5)
print(f"   chance is 0.5; the binomial 95 % band on {len(yc)} epochs is {clo:.4f} to {chi:.4f}")
print(f"\n   {lk['what_leaks']}")

t0 = time.time()
survey = []
for sid in CSP_CANDIDATES:
    t1 = time.time()
    try:
        Xs, ys, _ = L6.csp_epochs(sid)
        r = L6.csp_accuracy(Xs, ys, n_components=int(L6.CSP_SPEC["n_components"]), folds=FOLDS, seed=SEED)
        survey.append({"subject": sid, "mean": r["mean"], "std": r["std"], "n": r["n_epochs"],
                       "scores": r["scores"]})
        print(f"  {sid}: {r['mean']:.4f} (n = {r['n_epochs']}) in {time.time() - t1:.0f} s", flush=True)
    except Exception as exc:                                     # noqa: BLE001
        print(f"  {sid}: skipped ({type(exc).__name__}: {exc})")
    finally:
        if DELETE_EEGBCI:
            L6.drop_eegbci_runs(sid, L6.CSP_SPEC["runs"], keep=PRE_EXISTING)
if DELETE_EEGBCI:
    L6.drop_eegbci_runs(SELECTED, L6.CSP_SPEC["runs"], keep=PRE_EXISTING)
means = np.array([r["mean"] for r in survey])
print(f"\n{len(survey)} of {len(CSP_CANDIDATES)} candidates in {time.time() - t0:.0f} s")
if len(survey) < len(CSP_CANDIDATES):
    print(f"   WARNING: {len(CSP_CANDIDATES) - len(survey)} candidate(s) could not be loaded, so the "
          f"summary below is NOT the full distribution.  Rerun before quoting it.")
L6.disk_report("after the ds-eegbci loop",
               folders={"ds-eegbci cache": helpers.download_root() / "MNE-eegbci-data"})
print(f"   median {np.median(means):.4f}, mean {means.mean():.4f}, range {means.min():.4f} to "
      f"{means.max():.4f}, IQR {np.percentile(means, 75) - np.percentile(means, 25):.4f}")
print(f"   the median and the mean differ by {abs(means.mean() - np.median(means)):.4f} because one subject "
      f"sits far above the rest; for a distribution this skewed the median is the summary and the range is "
      f"the caveat.")
print(f"   the selected subject, {SELECTED}: {dict((r['subject'], r['mean']) for r in survey).get(SELECTED, float('nan')):.4f} "
      f"-- the {100 * float(np.mean(means <= dict((r['subject'], r['mean']) for r in survey).get(SELECTED, np.nan))):.0f}th percentile of its own candidate pool")
print(f"   candidates inside the chance band {clo:.3f}-{chi:.3f}: "
      f"{int(((means >= clo) & (means <= chi)).sum())} of {len(means)}")
print()
print(f"   Quoting {dict((r['subject'], r['mean']) for r in survey).get(SELECTED, float('nan')):.1%} as "
      f"'what motor-imagery decoding achieves' would be wrong by about "
      f"{100 * (dict((r['subject'], r['mean']) for r in survey).get(SELECTED, np.nan) - np.median(means)):.0f} "
      f"percentage points.  The median is the honest summary and the range is the honest caveat.")

fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.2))
order = np.argsort(means)
axes[0].bar(range(len(survey)), means[order],
            color=["tab:orange" if survey[i]["subject"] == SELECTED else "tab:blue" for i in order])
axes[0].axhline(0.5, color="gray", ls="--", lw=1, label="chance")
axes[0].axhspan(clo, chi, color="gray", alpha=0.18, label="binomial 95 % band around chance")
axes[0].axhline(float(np.median(means)), color="k", ls=":", lw=1.2,
                label=f"median {np.median(means):.3f}")
axes[0].set_xticks(range(len(survey)), [survey[i]["subject"] for i in order], rotation=60, fontsize=7)
axes[0].set(ylabel="Cross-validated accuracy (proportion correct)",
            title=f"The same pipeline on {len(survey)} candidates; the selected one in orange")
axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3, axis="y")

axes[1].bar(["honest\n(CSP inside the fold)", "leaked\n(CSP outside the fold)"],
            [lk["honest"], lk["leaked"]], color=["tab:blue", "tab:orange"], width=0.5)
axes[1].axhline(0.5, color="gray", ls="--", lw=1)
axes[1].axhspan(clo, chi, color="gray", alpha=0.18)
for i, v in enumerate([lk["honest"], lk["leaked"]]):
    axes[1].text(i, v + 0.015, f"{v:.3f}", ha="center", fontsize=11)
axes[1].annotate("", xy=(1, lk["leaked"]), xytext=(1, lk["honest"]),
                 arrowprops=dict(arrowstyle="<->", color="k"))
axes[1].text(1.06, (lk["honest"] + lk["leaked"]) / 2, f"{100 * lk['inflation']:+.1f} pp", fontsize=10)
axes[1].set(ylim=(0, 1.12), ylabel="Cross-validated accuracy (proportion correct)",
            title=f"{SELECTED}: the same data, the same folds, one line moved")
axes[1].grid(alpha=0.3, axis="y")
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

print("nb-6-5-decoding -- L6.5 numbers (draft; TODO(confirm) at author review)")
print(f"Data: ds-erpcore P3, {len(sids)} subjects (CC-BY-SA-4.0, contested at source) and ds-eegbci "
      f"R04+R08+R12 (ODC-By-1.0)")
print(f"Pipeline: helpers_l6.PRESPECIFIED for the ERP CORE epochs; helpers_l6.CSP_SPEC for ds-eegbci")
print()
print(f"1. TIME-RESOLVED DECODING ({len(sids)} subjects, {len(times)} time points at "
      f"{1 / np.diff(times).mean():.0f} Hz, {FOLDS}-fold StratifiedKFold within subject, ROC AUC):")
print(f"     group peak AUC {A.mean(0)[peak_i]:.4f} at {times[peak_i] * 1000:.0f} ms")
print(f"     per-subject peak AUC: median {np.median([auc[s].max() for s in sids]):.4f}, "
      f"range {min(auc[s].max() for s in sids):.4f} to {max(auc[s].max() for s in sids):.4f}")
print()
print(f"2. ex-6-5 ANSWER PAIR -- peak decoding time and whether it is significant:")
print(f"     peak at {times[peak_i] * 1000:.0f} ms, AUC {A.mean(0)[peak_i]:.4f}")
print(f"     cluster permutation (threshold {thr:.4f}, tail 1, {len(H0)} permutations, seed {SEED}): "
      f"{len(sig_clusters)} significant cluster(s)")
for c in sig_clusters:
    print(f"        {c['t_start_ms']:.0f} to {c['t_end_ms']:.0f} ms, {c['n']} points, p = {c['p']:.4f}")
print(f"     the peak IS{'' if mask[peak_i] else ' NOT'} inside a significant cluster")
print()
print(f"3. CHANCE: AUC chance 0.500; majority-class accuracy {maj:.3f} on a median subject; binomial 95 % "
      f"band {lo:.3f}-{hi:.3f} on {n_tr} trials")
print(f"     empirical null for {s_demo} ({N_PERM_WITHIN} label shuffles): mean {null.mean():.4f}, "
      f"SD {null.std(ddof=1):.4f}; observed {obs:.4f}, p = {p_emp:.4f}")
print()
print(f"4. TEMPORAL GENERALIZATION ({s_demo}, {len(tg)} x {len(tg)} at "
      f"{1 / np.diff(tg).mean():.0f} Hz): diagonal peak {d.max():.4f} at "
      f"{tg[int(np.argmax(d))] * 1000:.0f} ms; off-diagonal/diagonal ratio after 0 ms "
      f"{(np.nanmean(off) - 0.5) / (d[post].mean() - 0.5):.3f}")
print()
print(f"5. LEAKAGE:")
print(f"     (a) pre-stimulus group mean AUC {A[:, pre].mean():.4f}, maximum {A.mean(0)[pre].max():.4f}; "
      f"t({len(sids) - 1}) = {t_pre:.3f} against 0.5, p = {p_pre:.4f}")
print(f"     (b) channel selection outside the fold: honest {np.mean([r['honest'] for r in rows]):.4f} vs "
      f"leaked {np.mean([r['leaked'] for r in rows]):.4f}, mean gap {gaps.mean():+.4f} AUC "
      f"(top {K_CHANNELS} of {len(L6.ERPCORE_30)} channels)")
print(f"     (c) CSP outside the fold, {SELECTED}: honest {lk['honest']:.4f} vs leaked {lk['leaked']:.4f}, "
      f"gap {100 * lk['inflation']:+.1f} percentage points")
print()
print(f"6. THE SELECTION DISCLOSURE -- the same pipeline on {len(survey)} of {len(CSP_CANDIDATES)} "
      f"candidates:")
print(f"     median {np.median(means):.4f}, mean {means.mean():.4f}, range {means.min():.4f} to "
      f"{means.max():.4f}")
print(f"     " + ", ".join(f"{r['subject']} {r['mean']:.4f}" for r in survey))
print(f"     {SELECTED} was selected by the rule 'lowest-numbered subject reaching 0.75'.  Its score must "
      f"never be quoted without the median and the range.")
print()
print(f"CROSS-CHECK against w-csp-explorer (csp.json, ds-eegbci S002, 8-30 Hz, 45 epochs):")
print(f"     widget:   86.7 % honest vs 97.8 % leaked, +11.1 points; folds "
      f"0.8889, 0.6667, 1.0, 1.0, 0.7778")
print(f"     notebook: {100 * lk['honest']:.1f} % honest vs {100 * lk['leaked']:.1f} % leaked, "
      f"{100 * lk['inflation']:+.1f} points; folds "
      f"{', '.join(f'{v:.4f}' for v in lk['honest_scores'])}")
print(f"     integration note site/notes/integration-phase3.md quotes a candidate 'median of 60.0 %'.  "
      f"The survey above gives median {100 * np.median(means):.1f} % and mean {100 * means.mean():.1f} % "
      f"over {len(survey)} candidates, so the published figure is the MEAN rather than the median (the "
      f"notebooks-L5 track reports the same).  The range matches exactly either way.  Which statistic the "
      f"lessons quote is a decision for the author: the median is the right summary of a distribution this "
      f"skewed, and both belong beside the selected subject's score.")
print()
print(f"Pitfall: pf-decoding-leakage.  No widget for this lesson; w-csp-explorer (L5.7) is the nearest.")