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

# 1. Dependencies are pinned in notebooks/requirements.txt.  Nothing is installed when the
#    pinned stack is already present (local runs, CI); a fresh Colab or Binder kernel installs
#    it once.  On Colab, run from a clone of the repository so that notebooks/_shared/ is
#    available (repository URL: TODO(confirm), spec section 13 item 3).
_needed = ("mne", "scipy", "matplotlib", "pandas", "pooch")
_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"]
    subprocess.check_call(_cmd)

# 2. Shared helpers, located relative to the working directory -- notebooks/<level>/ or
#    notebooks/ -- never through an absolute path.
_shared = next((d / "_shared" for d in (Path.cwd(), *Path.cwd().parents)
                if (d / "_shared" / "helpers_l3.py").exists()), None)
if _shared is None:
    raise FileNotFoundError("start the kernel in notebooks/L3/ (or notebooks/) so that _shared/helpers_l3.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l3 as L3

# 3. Plotting: Jupyter's default inline backend renders static PNGs through Agg (no windows,
#    nothing blocks); outside Jupyter the helpers select Agg.  Every MNE figure is requested
#    with show=False and each figure cell ends with plt.show().
import matplotlib.pyplot as plt
import numpy as np
import mne

mne.set_log_level("WARNING")
plt.rcParams["figure.dpi"] = 72
print(f"MNE {mne.__version__}; helpers_l3 imported from notebooks/_shared")
print(f"ERP CORE cache: {L3.erpcore_root().name}/ (resolved relative to the working directory, "
      "or $EEG_COURSE_ERPCORE); only the subjects this notebook names are fetched")

FULL_COHORT = False               # True -> all 40 subjects; see the disk check below

SUBJECTS = list(range(1, 41)) if FULL_COHORT else list(range(1, 21))
W = L3.P3_WINDOW
CH = L3.P3_CHANNEL
ALPHA = 0.05
SEED = 20260918
SME_THRESHOLD_UV = 1.0            # the threshold nb-3-5 states and uses

root = L3.erpcore_root()
have = [s for s in SUBJECTS if (root / "p3" / L3.erpcore_subject_id(s)).is_dir()
        or (root / "P3" / L3.erpcore_subject_id(s)).is_dir()]
missing = [s for s in SUBJECTS if s not in have]
free_mb = helpers.free_disk_mb(root)
need_mb = 58.0 * len(missing)
print(f"cohort: {'FULL (40 subjects)' if FULL_COHORT else f'documented subset ({len(SUBJECTS)} subjects)'}")
print(f"  already cached: {len(have)}; to download: {len(missing)} (~{need_mb:.0f} MB at ~58 MB each)")
print(f"  free disk where the cache lives: {free_mb:.0f} MB")
if need_mb > 0 and need_mb + 500 > free_mb:
    raise SystemExit(f"refusing to download {need_mb:.0f} MB with only {free_mb:.0f} MB free; "
                     f"set FULL_COHORT = False, or free space, or point $EEG_COURSE_ERPCORE elsewhere")
print("  proceeding" + ("" if missing else " -- nothing to download"))

# The one Level-3 pipeline, printed rather than described.  Every Level-3 notebook and the C3
# capstone call the same helpers_l3.load_p3_epochs, so their numbers are comparable.
for key, value in L3.P3_PIPELINE.items():
    print(f"{key:15s} : {value}")
print()
print(f"a-priori measurement window : {L3.P3_WINDOW[0] * 1000:.0f}-{L3.P3_WINDOW[1] * 1000:.0f} ms "
      f"at {L3.P3_CHANNEL}, fixed in helpers_l3.P3_WINDOW")

import time as _time

t0 = _time.time()
store, qc = {}, []
for s in SUBJECTS:
    ep, nfo = L3.load_p3_epochs(s, verbose=False)
    store[s] = {"target": L3.condition_epochs(ep, "target"), "standard": L3.condition_epochs(ep, "standard"),
                "info": nfo}
    qc.append(nfo)
    print(f"  {nfo['subject']}: {nfo['n_kept']['target']:3d} target / {nfo['n_kept']['standard']:3d} standard "
          f"kept, {nfo['n_rejected']:3d} rejected; bads {nfo['bad_channels'] or '[]'}; "
          f"ICA excluded {nfo['ica_excluded']}; RT target {nfo['mean_rt_ms']['target']:.0f} ms; "
          f"accuracy {100 * nfo['accuracy']:.1f} %", flush=True)
times = ep.times
eeg = qc[0]["eeg_channels"]
i_ch = eeg.index(CH)
proto = store[SUBJECTS[0]]["target"].average().copy().pick("eeg")
print(f"\n{len(SUBJECTS)} subjects processed from raw in {_time.time() - t0:.1f} s")

n_t = np.array([q["n_kept"]["target"] for q in qc])
n_s = np.array([q["n_kept"]["standard"] for q in qc])
n_rej = np.array([q["n_rejected"] for q in qc])
n_ica = np.array([len(q["ica_excluded"]) for q in qc])
n_bad = np.array([len(q["bad_channels"]) for q in qc])
acc = np.array([q["accuracy"] for q in qc])
rt_t = np.array([q["mean_rt_ms"]["target"] for q in qc])
rt_s = np.array([q["mean_rt_ms"]["standard"] for q in qc])

print(f"cohort summary ({len(SUBJECTS)} subjects, 200 stimulus events each by design: 40 target, 160 standard)")
print(f"  epochs rejected at {L3.REJECT_PTP_UV:g} uV peak-to-peak: median {np.median(n_rej):.0f} "
      f"({100 * np.median(n_rej) / 200:.1f} %), range {n_rej.min()}-{n_rej.max()}")
print(f"  target trials kept   : median {np.median(n_t):.0f}, range {n_t.min()}-{n_t.max()}")
print(f"  standard trials kept : median {np.median(n_s):.0f}, range {n_s.min()}-{n_s.max()}")
print(f"  condition-biased rejection (|% target rejected - % standard rejected|): median "
      f"{np.median(np.abs(100 * (1 - n_t / 40) - 100 * (1 - n_s / 160))):.1f} points, "
      f"max {np.max(np.abs(100 * (1 - n_t / 40) - 100 * (1 - n_s / 160))):.1f} points")
print(f"  bad channels interpolated: {int(n_bad.sum())} across the cohort "
      f"({int((n_bad > 0).sum())} subjects); ICA components removed as ocular: median {np.median(n_ica):.0f}, "
      f"range {n_ica.min()}-{n_ica.max()}")
print(f"  behaviour: accuracy {100 * acc.mean():.1f} % (range {100 * acc.min():.1f}-{100 * acc.max():.1f}); "
      f"reaction time target {rt_t.mean():.0f} ms, standard {rt_s.mean():.0f} ms "
      f"(target - standard {rt_t.mean() - rt_s.mean():+.0f} ms)")
from scipy import stats
t_rt, p_rt = stats.ttest_rel(rt_t, rt_s)
print(f"    paired t test on reaction time, target vs standard: t({len(SUBJECTS) - 1}) = {t_rt:.3f}, "
      f"p = {p_rt:.5f} -- the behavioural oddball effect, independent of the EEG")

fig, axes = plt.subplots(1, 3, figsize=(14, 3.6))
axes[0].bar(np.arange(len(SUBJECTS)), n_rej, color="tab:blue")
axes[0].set(xlabel="Subject (index in the subset)", ylabel="Epochs rejected (of 200)",
            title=f"Artifact rejection at {L3.REJECT_PTP_UV:g} uV peak-to-peak")
axes[1].bar(np.arange(len(SUBJECTS)) - 0.2, n_t, 0.4, label="target")
axes[1].bar(np.arange(len(SUBJECTS)) + 0.2, n_s / 4, 0.4, label="standard / 4")
axes[1].set(xlabel="Subject (index in the subset)", ylabel="Trials kept",
            title="Trials kept per condition (standard scaled by 1/4)")
axes[1].legend(fontsize=8)
axes[2].plot(rt_s, rt_t, "o")
lims = [min(rt_s.min(), rt_t.min()) - 20, max(rt_s.max(), rt_t.max()) + 20]
axes[2].plot(lims, lims, "k--", lw=0.8)
axes[2].set(xlabel="Reaction time, standard (ms)", ylabel="Reaction time, target (ms)",
            title="Behaviour: every point above the line is a slower target")
for ax in axes:
    ax.grid(alpha=0.3)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

X_all = np.stack([(store[s]["target"].average().copy().pick("eeg").data
                   - store[s]["standard"].average().copy().pick("eeg").data) * 1e6 for s in SUBJECTS])
X = X_all[:, i_ch, :]
grand_t = np.mean([store[s]["target"].average().copy().pick("eeg").data for s in SUBJECTS], axis=0) * 1e6
grand_s = np.mean([store[s]["standard"].average().copy().pick("eeg").data for s in SUBJECTS], axis=0) * 1e6
n = len(SUBJECTS)

amp_mean = np.array([L3.mean_amplitude(X[k], None, W, times=times) for k in range(n)])
amp_peak = np.array([L3.peak_amplitude(X[k], None, W, times=times) for k in range(n)])
t_stat, p_val = stats.ttest_1samp(amp_mean, 0)
dz = amp_mean.mean() / amp_mean.std(ddof=1)
ci = stats.t.ppf(1 - ALPHA / 2, n - 1) * amp_mean.std(ddof=1) / np.sqrt(n)

print(f"P3 amplitude, {CH}, {W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms, target minus standard, {n} subjects")
print(f"  mean amplitude : {amp_mean.mean():+.3f} uV  (SD {amp_mean.std(ddof=1):.3f}, "
      f"SEM {amp_mean.std(ddof=1) / np.sqrt(n):.3f}, 95% CI [{amp_mean.mean() - ci:+.3f}, "
      f"{amp_mean.mean() + ci:+.3f}])")
print(f"  peak amplitude : {amp_peak.mean():+.3f} uV  (SD {amp_peak.std(ddof=1):.3f}) -- reported for "
      f"comparison only; it is biased upward by noise")
print(f"  grand-average difference wave at {CH}: mean "
      f"{L3.mean_amplitude(grand_t[i_ch] - grand_s[i_ch], None, W, times=times):+.3f} uV, peak "
      f"{L3.peak_amplitude(grand_t[i_ch] - grand_s[i_ch], None, W, times=times):+.3f} uV")
print(f"  per condition at {CH}: target {L3.mean_amplitude(grand_t[i_ch], None, W, times=times):+.3f} uV, "
      f"standard {L3.mean_amplitude(grand_s[i_ch], None, W, times=times):+.3f} uV")
print(f"  {int((amp_mean > 0).sum())} of {n} subjects show a positive effect")

fig, axes = plt.subplots(1, 2, figsize=(13, 4.4))
L3.plot_erp({f"target (n = {int(n_t.sum())} trials)": grand_t[i_ch],
             f"standard (n = {int(n_s.sum())} trials)": grand_s[i_ch],
             "difference": ((grand_t - grand_s)[i_ch], {"color": "k", "lw": 1.8})},
            times, window=W, ax=axes[0],
            title=f"C3: grand-average P3 at {CH}, {n} ERP CORE subjects (uV, positive up)")
axes[1].hist(amp_mean, bins=max(6, n // 2), color="tab:blue", alpha=0.85)
axes[1].axvline(0, color="gray", lw=1.0)
axes[1].axvline(amp_mean.mean(), color="tab:orange", lw=2,
                label=f"mean {amp_mean.mean():+.2f} uV, 95% CI +/-{ci:.2f}")
axes[1].set(xlabel=f"P3 mean amplitude at {CH} (uV)", ylabel="Subjects",
            title=f"Distribution across {n} subjects (uV)")
axes[1].grid(alpha=0.3, axis="y")
axes[1].legend(fontsize=8)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

lat_peak = np.array([L3.peak_latency(X[k], None, W, times=times) for k in range(n)])
lat_fal = np.array([L3.fractional_area_latency(X[k], None, W, times=times) for k in range(n)])
g_diff = (grand_t - grand_s)[i_ch]

print(f"P3 latency, {CH}, measured inside the same a-priori window "
      f"({W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms)")
print(f"  peak latency, per subject then averaged      : {1000 * lat_peak.mean():.1f} ms "
      f"(SD {1000 * lat_peak.std(ddof=1):.1f}, range {1000 * lat_peak.min():.0f}-{1000 * lat_peak.max():.0f})")
print(f"  50% fractional-area latency, per subject      : {1000 * lat_fal.mean():.1f} ms "
      f"(SD {1000 * lat_fal.std(ddof=1):.1f}, range {1000 * lat_fal.min():.0f}-{1000 * lat_fal.max():.0f})")
print(f"  peak latency of the GRAND AVERAGE            : {1000 * L3.peak_latency(g_diff, None, W, times=times):.1f} ms")
print(f"  50% area latency of the GRAND AVERAGE        : "
      f"{1000 * L3.fractional_area_latency(g_diff, None, W, times=times):.1f} ms")
print(f"  the grand-average latency and the average of the per-subject latencies are different quantities "
      f"and neither is 'the' latency; both are reported because a paper that quotes one and a replication "
      f"that computes the other will not agree.")
print(f"  correlation between a subject's P3 latency and its mean reaction time: "
      f"r = {np.corrcoef(lat_fal, rt_t)[0, 1]:+.3f} (peak latency: "
      f"{np.corrcoef(lat_peak, rt_t)[0, 1]:+.3f}); nb-3-6 does this within subject, trial by trial")

m_win = (times >= W[0]) & (times <= W[1])
# get_data returns volts; every score in this notebook is in microvolts
scores = {s: {c: store[s][c].get_data(picks="eeg")[:, i_ch, :][:, m_win].mean(1) * 1e6
              for c in ("target", "standard")} for s in SUBJECTS}
print(f"single-trial mean amplitudes at {CH} over {W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms: "
      f"cohort trial SD {np.mean([scores[s]['target'].std(ddof=1) for s in SUBJECTS]):.2f} uV (target)")
sme_t = np.array([L3.sme(scores[s]["target"]) for s in SUBJECTS])
sme_s = np.array([L3.sme(scores[s]["standard"]) for s in SUBJECTS])
sme_d = np.sqrt(sme_t ** 2 + sme_s ** 2)
need = np.array([int(np.ceil((scores[s]["target"].std(ddof=1) / SME_THRESHOLD_UV) ** 2)) for s in SUBJECTS])
reach = int(np.sum(need <= n_t))

print(f"effect size ({n} subjects, within-subject contrast)")
print(f"  Cohen dz (mean / SD of the per-subject differences) : {dz:.3f}")
print(f"  t({n - 1}) = {t_stat:.3f}, p = {p_val:.3e}")
print(f"  95% CI on the mean difference                        : "
      f"[{amp_mean.mean() - ci:+.3f}, {amp_mean.mean() + ci:+.3f}] uV")
print()
print(f"standardized measurement error (mean amplitude at {CH}, {W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms)")
print(f"  target condition   : median {np.median(sme_t):.3f} uV, range {sme_t.min():.3f}-{sme_t.max():.3f}")
print(f"  standard condition : median {np.median(sme_s):.3f} uV, range {sme_s.min():.3f}-{sme_s.max():.3f}")
print(f"  the difference     : median {np.median(sme_d):.3f} uV (the two errors add in quadrature)")
print(f"  at the stated threshold SME <= {SME_THRESHOLD_UV:g} uV: median {int(np.median(need))} target trials "
      f"needed, range {need.min()}-{need.max()}; {reach}/{n} subjects reach it with the trials they have")
print(f"  measurement error accounts for about "
      f"{100 * np.mean(sme_d ** 2) / amp_mean.var(ddof=1):.0f} % of the between-subject variance in the "
      f"effect; the remaining {100 - 100 * np.mean(sme_d ** 2) / amp_mean.var(ddof=1):.0f} % is real "
      f"between-subject difference")

fig, ax = plt.subplots(figsize=(10, 4.4))
order = np.argsort(amp_mean)
ax.errorbar(np.arange(n), amp_mean[order], yerr=sme_d[order], fmt="o", capsize=3, ms=5)
ax.axhline(0, color="gray", lw=0.8)
ax.axhline(amp_mean.mean(), color="tab:orange", lw=1.4, ls="--",
           label=f"cohort mean {amp_mean.mean():+.2f} uV (95% CI +/-{ci:.2f})")
ax.fill_between([-0.5, n - 0.5], amp_mean.mean() - ci, amp_mean.mean() + ci, color="tab:orange", alpha=0.15)
ax.set_xticks(np.arange(n))
ax.set_xticklabels([qc[i]["subject"][-3:] for i in order], fontsize=7, rotation=90)
ax.set(xlabel="Subject (sorted by effect)", ylabel="P3, target - standard (uV)", xlim=(-0.5, n - 0.5),
       title=f"C3: every subject's P3 with its standardized measurement error "
             f"({CH}, {W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms, uV)")
ax.grid(alpha=0.3, axis="y")
ax.legend(fontsize=8)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

threshold = float(stats.t.ppf(1 - ALPHA / 2, n - 1))
N_PERM = 10000
t_obs, clusters, cluster_p, H0 = mne.stats.permutation_cluster_1samp_test(
    X, threshold=threshold, n_permutations=N_PERM, tail=0, out_type="indices", seed=SEED, verbose=False)

print(f"mne.stats.permutation_cluster_1samp_test(X, threshold={threshold:.4f}, n_permutations={N_PERM}, "
      f"tail=0, out_type='indices', seed={SEED})")
print(f"  X: {X.shape[0]} subjects x {X.shape[1]} time points of the {CH} difference wave (uV)")
print(f"  permutations used {len(H0)} of the {2 ** n if n <= 25 else 'many'} possible sign flips; "
      f"smallest attainable p-value {1 / len(H0):.5f}")
rows = []
for i in np.argsort(cluster_p):
    idx = np.asarray(clusters[i][0])
    rows.append({"start_ms": float(times[idx[0]] * 1000), "end_ms": float(times[idx[-1]] * 1000),
                 "n": int(len(idx)), "t_sum": float(t_obs[idx].sum()), "p": float(cluster_p[i]),
                 "sign": "positive" if t_obs[idx].sum() > 0 else "negative"})
for r in rows:
    print(f"  {r['sign']:8s} cluster {r['start_ms']:7.1f} to {r['end_ms']:7.1f} ms ({r['n']:3d} samples): "
          f"t-sum {r['t_sum']:+9.2f}, p = {r['p']:.4f}"
          + ("   <- significant" if r["p"] <= ALPHA else ""))
sig_clusters = [r for r in rows if r["p"] <= ALPHA]

adjacency, _ = mne.channels.find_ch_adjacency(proto.info, ch_type="eeg")
t_st, cl_st, p_st, H0_st = mne.stats.spatio_temporal_cluster_1samp_test(
    np.transpose(X_all, (0, 2, 1)), threshold=threshold, n_permutations=N_PERM, tail=0,
    adjacency=adjacency, out_type="mask", seed=SEED, verbose=False)
st_rows = []
for i in np.argsort(p_st):
    mask = cl_st[i]
    ti, ci_ = np.where(mask.any(axis=1))[0], np.where(mask.any(axis=0))[0]
    st_rows.append({"start_ms": float(times[ti[0]] * 1000), "end_ms": float(times[ti[-1]] * 1000),
                    "n_ch": int(len(ci_)), "t_sum": float(t_st[mask].sum()), "p": float(p_st[i]),
                    "channels": [eeg[j] for j in ci_]})
print(f"\nspatio-temporal over {len(eeg)} channels "
      f"(mne.stats.spatio_temporal_cluster_1samp_test, same threshold/seed/permutations, adjacency from "
      f"mne.channels.find_ch_adjacency):")
for r in st_rows[:4]:
    print(f"  {r['start_ms']:7.1f} to {r['end_ms']:7.1f} ms, {r['n_ch']:2d} channels: t-sum "
          f"{r['t_sum']:+9.1f}, p = {r['p']:.4f}"
          + ("   <- significant" if r["p"] <= ALPHA else ""))
    print(f"      {', '.join(r['channels'])}")
st_sig = [r for r in st_rows if r["p"] <= ALPHA]

fig, ax = plt.subplots(figsize=(10, 4.2))
L3.plot_cluster_test(times, t_obs, clusters, cluster_p, alpha=ALPHA, threshold=threshold, ax=ax,
                     title=f"C3 cluster permutation test at {CH} ({n} subjects, {len(H0)} permutations, "
                           f"seed {SEED})")
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

published = {
    # Fill each value from the ERP CORE paper (Kappenman et al. 2020, DOI 10.31234/osf.io/4azqm; dataset DOI
    # 10.18112/openneuro.ds003069.v1.0.0, both now in data/directory.yaml).  The catalog holds the citation but
    # no result values, so these stay None -- and None keeps the TODO(confirm) marker -- until read from the
    # paper itself (spec section 0.3: no scientific number from memory).
    "mean amplitude, target - standard (uV)": None,
    "peak latency (ms)": None,
    "effect size (Cohen dz)": None,
    "SME of the mean amplitude (uV)": None,
    "n subjects": None,
    "measurement window (ms)": None,
    "measurement channel": None,
    "reference": None,
}
ours = {
    "mean amplitude, target - standard (uV)": amp_mean.mean(),
    "peak latency (ms)": 1000 * lat_peak.mean(),
    "effect size (Cohen dz)": dz,
    "SME of the mean amplitude (uV)": float(np.median(sme_t)),
    "n subjects": float(n),
    "measurement window (ms)": f"{W[0] * 1000:.0f}-{W[1] * 1000:.0f}",
    "measurement channel": CH,
    "reference": "average of the 30 EEG channels",
}

print(f"{'quantity':42s} {'this replication':>22s} {'published':>18s} {'divergence':>14s}")
for key in ours:
    pub = published[key]
    mine = ours[key]
    mine_s = f"{mine:.3f}" if isinstance(mine, float) else str(mine)
    if pub is None:
        print(f"{key:42s} {mine_s:>22s} {'TODO(confirm)':>18s} {'TODO(confirm)':>14s}")
    elif isinstance(pub, (int, float)) and isinstance(mine, float):
        print(f"{key:42s} {mine_s:>22s} {pub:18.3f} {mine - pub:+14.3f}")
    else:
        print(f"{key:42s} {mine_s:>22s} {str(pub):>18s} {'(match)' if str(pub) == str(mine) else '(differs)':>14s}")
print()
print("Every published cell is TODO(confirm) by design.  The catalog (data/directory.yaml) now carries the "
      "citation and both DOIs -- paper 10.31234/osf.io/4azqm, dataset "
      "10.18112/openneuro.ds003069.v1.0.0 -- but it carries no published amplitudes, latencies or effect "
      "sizes, and spec section 0.3 forbids supplying those from memory.  Fill `published` above from the "
      "paper; nothing else in this notebook changes.")

i9, i10 = eeg.index("P9"), eeg.index("P10")
variants = {
    "this notebook (average reference, 300-600 ms, mean amplitude)": amp_mean,
    "linked-mastoid stand-in (P9/P10), same window":
        np.array([L3.mean_amplitude((X_all[k] - 0.5 * (X_all[k, [i9]] + X_all[k, [i10]]))[i_ch],
                                    None, W, times=times) for k in range(n)]),
    "window 300-500 ms (the window nb-1-5-filters used)":
        np.array([L3.mean_amplitude(X[k], None, (0.30, 0.50), times=times) for k in range(n)]),
    "window 400-700 ms": np.array([L3.mean_amplitude(X[k], None, (0.40, 0.70), times=times) for k in range(n)]),
    "channel CPz instead of Pz":
        np.array([L3.mean_amplitude(X_all[k, eeg.index("CPz")], None, W, times=times) for k in range(n)]),
    "peak amplitude instead of mean amplitude": amp_peak,
}
print(f"the same {n} subjects measured six ways")
print(f"{'variant':62s} {'mean (uV)':>10s} {'SD':>8s} {'dz':>7s} {'vs this notebook':>18s}")
for label, v in variants.items():
    t_, p_ = stats.ttest_1samp(v, 0)
    delta = v.mean() - amp_mean.mean()
    print(f"{label:62s} {v.mean():+10.3f} {v.std(ddof=1):8.3f} {v.mean() / v.std(ddof=1):7.3f} "
          f"{delta:+18.3f}")
print(f"\nthe spread across these six choices is "
      f"{max(v.mean() for v in variants.values()) - min(v.mean() for v in variants.values()):.3f} uV -- "
      f"which is the size of the divergence the analysis alone can produce, and the number any discussion "
      f"of divergence has to beat before invoking anything about the sample.")

alt_ref = variants["linked-mastoid stand-in (P9/P10), same window"].mean()
alt_win = variants["window 300-500 ms (the window nb-1-5-filters used)"].mean()
alt_pk = variants["peak amplitude instead of mean amplitude"].mean()
spread = max(v.mean() for v in variants.values()) - min(v.mean() for v in variants.values())
cluster_sentence = (
    f"A cluster-based permutation test on the subject-level difference waves at {CH} returned "
    f"{len(sig_clusters)} significant cluster(s) "
    + "; ".join(f"({r['start_ms']:.0f}-{r['end_ms']:.0f} ms, t-sum {r['t_sum']:+.0f}, p = {r['p']:.4f})"
                for r in sig_clusters)
    + f", which establishes that the conditions differ somewhere in the tested window and says nothing about "
      f"when the difference begins or ends."
) if sig_clusters else (
    f"A cluster-based permutation test on the subject-level difference waves at {CH} returned no significant "
    f"cluster, which is not evidence that the conditions are the same.")

print("nb-c3-replicate-erp-core -- C3 capstone numbers (draft; TODO(confirm) at author review)")
print(f"Cohort: ds-erpcore P3, sub-001 to sub-{SUBJECTS[-1]:03d} ({n} subjects; FULL_COHORT = {FULL_COHORT}; "
      f"set it True for all 40); CC-BY-SA-4.0 per data/directory.yaml (contested at source; spec 10.7 most-restrictive rule)")
print(f"Pipeline: helpers_l3.P3_PIPELINE, unchanged (printed in section 2)")
print(f"Measurement: mean amplitude, {CH}, {W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms, a-priori "
      f"(helpers_l3.P3_WINDOW)")
print()
print(f"  AMPLITUDE       {amp_mean.mean():+.3f} uV  (SD {amp_mean.std(ddof=1):.3f}, 95% CI "
      f"[{amp_mean.mean() - ci:+.3f}, {amp_mean.mean() + ci:+.3f}]); peak amplitude "
      f"{amp_peak.mean():+.3f} uV")
print(f"  LATENCY         peak {1000 * lat_peak.mean():.1f} ms (SD {1000 * lat_peak.std(ddof=1):.1f}); "
      f"50% fractional-area {1000 * lat_fal.mean():.1f} ms (SD {1000 * lat_fal.std(ddof=1):.1f}); "
      f"grand-average peak {1000 * L3.peak_latency(g_diff, None, W, times=times):.1f} ms")
print(f"  EFFECT SIZE     Cohen dz = {dz:.3f}; t({n - 1}) = {t_stat:.3f}, p = {p_val:.3e}")
print(f"  SME             median {np.median(sme_t):.3f} uV (target), {np.median(sme_s):.3f} uV (standard), "
      f"{np.median(sme_d):.3f} uV (difference); {reach}/{n} subjects reach SME <= "
      f"{SME_THRESHOLD_UV:g} uV with the trials they have")
print(f"  CLUSTER TEST    {CH}, threshold t = {threshold:.4f}, {len(H0)} permutations, seed {SEED}:")
for r in rows:
    print(f"                    {r['sign']:8s} {r['start_ms']:7.1f} to {r['end_ms']:7.1f} ms, t-sum "
          f"{r['t_sum']:+9.2f}, p = {r['p']:.4f}"
          + ("   <- significant" if r["p"] <= ALPHA else ""))
print(f"                  spatio-temporal, {len(eeg)} channels:")
for r in st_rows[:3]:
    print(f"                    {r['start_ms']:7.1f} to {r['end_ms']:7.1f} ms, {r['n_ch']:2d} channels, "
          f"t-sum {r['t_sum']:+9.1f}, p = {r['p']:.4f}"
          + ("   <- significant" if r["p"] <= ALPHA else ""))
print(f"  BEHAVIOUR       reaction time target {rt_t.mean():.0f} ms vs standard {rt_s.mean():.0f} ms "
      f"(t({n - 1}) = {t_rt:.3f}, p = {p_rt:.3e}); accuracy {100 * acc.mean():.1f} %")
print(f"  QC              median {np.median(n_rej):.0f}/200 epochs rejected; target trials kept median "
      f"{np.median(n_t):.0f} (range {n_t.min()}-{n_t.max()}); {int(n_bad.sum())} channels interpolated; "
      f"median {np.median(n_ica):.0f} ICA components removed")
print(f"  PUBLISHED       TODO(confirm) -- not in data/directory.yaml or data/catalog/; every cell of the "
      f"comparison table in section 7 is TODO(confirm) and none was supplied from memory")
print()
print("PARAGRAPH TEMPLATE, filled with this run's values (published values stay TODO(confirm)):")
print()
print(f"  We measured the P3 oddball effect as {amp_mean.mean():+.3f} uV (95% CI "
      f"[{amp_mean.mean() - ci:+.3f}, {amp_mean.mean() + ci:+.3f}], dz = {dz:.3f}, n = {n}) as the mean "
      f"amplitude of the target-minus-standard difference wave over "
      f"{W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms at {CH}, against a published value of TODO(confirm); the "
      f"difference is TODO(confirm) uV. Before attributing it to the sample, three analysis choices were "
      f"checked on our own data: the reference (a linked-mastoid stand-in gives {alt_ref:+.3f} uV, a "
      f"difference of {alt_ref - amp_mean.mean():+.3f} uV), the measurement window ({alt_win:+.3f} uV over "
      f"300-500 ms) and the measure itself (peak amplitude gives {alt_pk:+.3f} uV). Across the six analysis "
      f"variants we tried the measured effect spans {spread:.3f} uV, so any divergence smaller than that is "
      f"explained by analysis choices alone. Our cohort was {n} subjects of the 40 the dataset ships, with a "
      f"median of {np.median(n_rej):.0f} of 200 epochs rejected by a {L3.REJECT_PTP_UV:g} uV peak-to-peak "
      f"criterion and a median of {np.median(n_t):.0f} target trials retained per subject; the median "
      f"standardized measurement error of the target mean amplitude was {np.median(sme_t):.3f} uV, so "
      f"{reach} of {n} subjects met our stated precision threshold of {SME_THRESHOLD_UV:g} uV. "
      f"{cluster_sentence} Remaining candidate explanations, in the order we would test them: measurement "
      f"window and reference (quantified above); artifact handling and trial counts (ours documented, the "
      f"published criterion TODO(confirm)); the subject set ({n} of 40 here, TODO(confirm) published); and "
      f"only then the population or the recording itself.")
print()
print("Rubric: (1) pipeline unchanged from C2 -- helpers_l3.P3_PIPELINE, printed in section 2, defaults "
      "only; (2) windows fixed a priori -- helpers_l3.P3_WINDOW / P3_CHANNEL, module constants shared with "
      "every Level-3 notebook; (3) cluster interpretation -- section 8 and nb-3-7.")