# 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', 'pooch', 'pyprep', 'autoreject')
_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"]
    if "pyprep" in _missing:
        _cmd += ["pyprep>=0.9"]
    if "autoreject" in _missing:
        _cmd += ["autoreject>=0.5"]
    subprocess.check_call(_cmd)

# 2. Shared helpers (notebooks/_shared/helpers.py and helpers_l2.py), 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_l2.py").exists()), None)
if _shared is None:
    raise FileNotFoundError("start the kernel in notebooks/L2/ (or notebooks/) so that _shared/helpers_l2.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l2 as l2

# 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 plt.show() renders each cell's
#    figures in place.
import matplotlib.pyplot as plt
import numpy as np
import mne

# Warnings are worth reading, so they are not silenced -- but their default format prints the
# absolute path of the file that raised them, which is nobody else's business and would put this
# machine's directory layout into the saved outputs.  Only the class and the message are shown.
warnings.formatwarning = lambda message, category, *a, **k: f"{category.__name__}: {message}\n"

mne.viz.set_browser_backend("matplotlib", verbose=False)
mne.set_log_level("WARNING")
plt.rcParams["figure.dpi"] = 72
print(f"MNE {mne.__version__}; helpers imported from notebooks/_shared")
print("ERP CORE downloads (~56 MB per subject) are cached under data/downloads/erpcore/ in a repository "
      "clone, otherwise under MNE's data directory; nothing is re-fetched.")

import time
from autoreject import get_rejection_threshold

SUBSET_N = 10
SUBJECTS = [f"sub-{i:03d}" for i in range(1, SUBSET_N + 1)]
RESAMPLE_HZ, L_FREQ, H_FREQ = 256.0, 0.1, 30.0
TMIN, TMAX, BASELINE = -0.2, 0.8, (-0.2, 0.0)
CH, WINDOW = l2.P3_CHANNEL, l2.P3_WINDOW
MASTOID_PAIR = ["P9", "P10"]

REFERENCES = {
    "original (CMS)": None,
    "average (30 ch)": "average",
    "linked mastoids (P9+P10)/2": MASTOID_PAIR,
    "Cz": ["Cz"],
}
print(f"subset: {SUBJECTS[0]} .. {SUBJECTS[-1]} ({len(SUBJECTS)} of 40 participants)")
print(f"pipeline: resample {RESAMPLE_HZ:g} Hz -> pyprep bad channels (>=2 criteria, cap 10 %) -> "
      f"filter {L_FREQ:g}-{H_FREQ:g} Hz FIR zero-phase -> interpolate -> re-reference -> epoch "
      f"{TMIN:g}..{TMAX:g} s, baseline {BASELINE[0]:g}..{BASELINE[1]:g} s")
print(f"trial set: fixed once on the average-referenced epochs with autoreject's global peak-to-peak threshold "
      f"(pooled trials, condition-blind); identical trials under every reference")
print(f"measure: target minus standard, mean amplitude at {CH} over {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms "
      f"(a-priori window)")
print(f"references: {list(REFERENCES)}")
print(f"seed {l2.SEED}")

t_start = time.time()
amp, evokeds, per_subject, skipped = {}, {}, [], []

for sid in SUBJECTS:
    try:
        raw, facts = l2.load_erpcore("P3", sid, resample_hz=RESAMPLE_HZ)
    except Exception as e:                      # a refused or missing subject: report, do not hide
        skipped.append(f"{sid}: {type(e).__name__}: {e}")
        continue
    n_eeg = len(mne.pick_types(raw.info, eeg=True))
    det = l2.detect_bad_channels(raw, highpass_hz=1.0, ransac=True, seed=l2.SEED)
    dec = l2.select_bads(det, n_eeg, min_criteria=2, max_fraction=0.10)
    raw.filter(L_FREQ, H_FREQ, picks=["eeg", "eog"], verbose=False)
    raw.info["bads"] = list(dec["bads"])
    if dec["bads"]:
        raw.interpolate_bads(reset_bads=True, verbose=False)
    events, _, ev_info = l2.p3_events(raw)

    # The trial set: decided once, on the average-referenced epochs, blind to condition.
    ep_avg = l2.epochs_p3(raw.copy().set_eeg_reference("average", verbose=False), events,
                          tmin=TMIN, tmax=TMAX, baseline=BASELINE)
    thr_uv = get_rejection_threshold(ep_avg, random_state=l2.SEED, verbose=False)["eeg"] * 1e6
    keep = l2.epoch_ptp_uv(ep_avg) <= thr_uv
    ptp_ch = np.ptp(ep_avg.get_data(picks="eeg"), axis=2) * 1e6          # epochs x channels
    worst = np.bincount(ptp_ch[~keep].argmax(axis=1), minlength=ptp_ch.shape[1]) if (~keep).any() else np.zeros(ptp_ch.shape[1], int)
    driver = ", ".join(f"{ep_avg.ch_names[i]} ({worst[i]})" for i in np.argsort(-worst)[:2] if worst[i]) or "-"

    row = {"subject": sid, "interpolated": dec["bads"] or ["-"], "rank": l2.data_rank(
        n_eeg, n_interpolated=len(dec["bads"]), average_reference=True)["rank"],
        "threshold_uv": float(thr_uv), "n_kept": int(keep.sum()), "n_trials": int(len(keep)),
        "n_target_kept": int((events[keep, 2] == l2.P3_EVENT_ID["target"]).sum()),
        "drops driven by": driver}
    for name, ref in REFERENCES.items():
        r = raw.copy() if ref is None else raw.copy().set_eeg_reference(ref, verbose=False)
        ep = l2.epochs_p3(r, events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)[keep]
        d = l2.difference_wave(ep)
        evokeds.setdefault(name, {})[sid] = {"difference": d, "target": ep["target"].average(),
                                             "standard": ep["standard"].average()}
        amp.setdefault(name, {})[sid] = l2.mean_amplitude(d, CH, WINDOW)
        row[name] = amp[name][sid]
    per_subject.append(row)
    print(f"  {sid}: interpolated {dec['bads'] or 'none'}, rank {row['rank']}, threshold {thr_uv:5.0f} uV, "
          f"{row['n_kept']:3d}/{row['n_trials']} trials kept ({row['n_target_kept']} targets), "
          f"drops driven by {driver}  |  "
          + "  ".join(f"{k.split(' ')[0]} {amp[k][sid]:+5.2f}" for k in REFERENCES))

DONE = [r["subject"] for r in per_subject]
print(f"\n{len(DONE)} subjects in {time.time() - t_start:.0f} s; skipped: {skipped or 'none'}")

print(l2.fmt_table(per_subject, ["subject", "interpolated", "rank", "threshold_uv", "n_kept", "drops driven by",
                                 *REFERENCES]))
print()
print("Read the 'drops driven by' column before anything else: on most subjects the channel that trips the "
      "threshold is Fp1 or Fp2, and the artefact is a blink. This pipeline deliberately stops before artifact "
      "correction -- ICA is L2.6 -- so a peak-to-peak criterion over all channels is, on this data, mostly a "
      "blink detector, and up to 40 % of trials go with it. That cost is the reason the canonical order puts "
      "ICA before rejection (L2.8): correct what can be corrected, reject what cannot. Nothing about the "
      "reference comparison changes, because every reference sees the same surviving trials.")
print()
summary = []
for name in REFERENCES:
    v = np.array([amp[name][s] for s in DONE])
    summary.append({"reference": name, "mean_uv": float(v.mean()), "sd_uv": float(v.std(ddof=1)),
                    "sem_uv": float(v.std(ddof=1) / np.sqrt(len(v))), "median_uv": float(np.median(v)),
                    "min_uv": float(v.min()), "max_uv": float(v.max()),
                    "n_positive": f"{int((v > 0).sum())}/{len(v)}"})
print(f"P3 mean amplitude at {CH}, target minus standard, {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms, "
      f"across {len(DONE)} subjects (uV):\n")
print(l2.fmt_table(summary, ["reference", "mean_uv", "sd_uv", "sem_uv", "median_uv", "min_uv", "max_uv",
                             "n_positive"], floatfmt="{:+.2f}"))

GA = {name: {k: mne.grand_average([evokeds[name][s][k] for s in DONE])
             for k in ("difference", "target", "standard")} for name in REFERENCES}

fig, axes = plt.subplots(1, 2, figsize=(14, 4.6), sharey=False)
for name in REFERENCES:
    d = GA[name]["difference"]
    axes[0].plot(d.times * 1000, d.data[d.ch_names.index(CH)] * 1e6, lw=1.4, label=name)
axes[0].axvspan(WINDOW[0] * 1000, WINDOW[1] * 1000, color="tab:orange", alpha=0.18, label="a-priori window")
axes[0].axhline(0, color="gray", lw=0.6); axes[0].axvline(0, color="gray", lw=0.6)
axes[0].set(xlabel="Time from stimulus (ms)", ylabel="Amplitude (uV)",
            title=f"Grand-average target minus standard at {CH} under four references "
                  f"(uV, positive up, n = {len(DONE)})")
axes[0].grid(alpha=0.3); axes[0].legend(fontsize=8)

for name, style in zip(REFERENCES, ("-", "-", "-", "-")):
    t_, s_ = GA[name]["target"], GA[name]["standard"]
    line, = axes[1].plot(t_.times * 1000, t_.data[t_.ch_names.index(CH)] * 1e6, style, lw=1.3, label=f"{name}: target")
    axes[1].plot(s_.times * 1000, s_.data[s_.ch_names.index(CH)] * 1e6, style, lw=1.0, alpha=0.5,
                 color=line.get_color(), label=f"{name}: standard")
axes[1].axvspan(WINDOW[0] * 1000, WINDOW[1] * 1000, color="tab:orange", alpha=0.18)
axes[1].axhline(0, color="gray", lw=0.6); axes[1].axvline(0, color="gray", lw=0.6)
axes[1].set(xlabel="Time from stimulus (ms)", ylabel="Amplitude (uV)",
            title=f"The two conditions separately at {CH} (uV, positive up)")
axes[1].grid(alpha=0.3); axes[1].legend(fontsize=7, ncol=2)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

maps = {}
for name in REFERENCES:
    d = GA[name]["difference"]
    m = (d.times >= WINDOW[0]) & (d.times <= WINDOW[1])
    maps[name] = d.data[:, m].mean(axis=1) * 1e6            # uV per channel, in d.ch_names order
names = list(REFERENCES)
print(f"pairwise correlation of the four maps after removing each map's spatial mean "
      f"(the reference-invariant part):\n")
rows = []
for a in names:
    row = {"map": a}
    for b in names:
        ca, cb = maps[a] - maps[a].mean(), maps[b] - maps[b].mean()
        row[b.split(" ")[0]] = float(np.corrcoef(ca, cb)[0, 1])
    rows.append(row)
print(l2.fmt_table(rows, ["map", *[n.split(" ")[0] for n in names]], floatfmt="{:.6f}"))
print()
print(f"and the part that does change -- where the map crosses zero:\n")
rows = [{"reference": n, "spatial mean (uV)": float(maps[n].mean()),
         "electrodes above zero": f"{int((maps[n] > 0).sum())}/{len(maps[n])}",
         f"value at {CH} (uV)": float(maps[n][GA[n]['difference'].ch_names.index(CH)]),
         "map range (uV)": float(maps[n].max() - maps[n].min())} for n in names]
print(l2.fmt_table(rows, ["reference", "spatial mean (uV)", "electrodes above zero", f"value at {CH} (uV)",
                          "map range (uV)"], floatfmt="{:+.2f}"))

vmax = max(np.abs(m).max() for m in maps.values())
info = mne.pick_info(GA[names[0]]["difference"].info,
                     mne.pick_types(GA[names[0]]["difference"].info, eeg=True))
fig, axes = plt.subplots(1, len(names), figsize=(3.3 * len(names), 3.8))
for ax, name in zip(axes, names):
    im, _ = mne.viz.plot_topomap(maps[name], info, axes=ax, show=False, contours=6, vlim=(-vmax, vmax))
    ax.set_title(f"{name}\n{CH} = {maps[name][GA[name]['difference'].ch_names.index(CH)]:+.2f} uV", fontsize=9)
cb = fig.colorbar(im, ax=axes.tolist(), shrink=0.7)
cb.set_label("uV")
fig.suptitle(f"Target minus standard, mean {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms, "
             f"one grand average under four references (uV, common colour scale)", y=1.02)
plt.show()   # render the static figure(s) of this cell inline

rows = []
for key in ("target", "standard", "difference"):
    vals = {name: np.array([l2.mean_amplitude(evokeds[name][s][key], CH, WINDOW) for s in DONE])
            for name in names}
    stack = np.vstack([vals[n] for n in names])                # references x subjects
    spread = stack.max(axis=0) - stack.min(axis=0)             # per subject, across references
    rows.append({"measured quantity": f"{key} at {CH}",
                 "mean across subjects, per reference": ", ".join(f"{n.split(' ')[0]} {vals[n].mean():+.2f}"
                                                                  for n in names),
                 "range across references (uV), median over subjects": float(np.median(spread)),
                 "max over subjects": float(spread.max())})
print(l2.fmt_table(rows, ["measured quantity", "mean across subjects, per reference",
                          "range across references (uV), median over subjects", "max over subjects"]))
print()
t_rng = rows[0]["range across references (uV), median over subjects"]
d_rng = rows[2]["range across references (uV), median over subjects"]
print(f"The target condition alone moves by {t_rng:.2f} uV (median over subjects) as the reference changes; "
      f"the difference wave moves by {d_rng:.2f} uV, a factor of {t_rng / d_rng:.1f} less. That is the "
      "practical argument for reporting a difference, and the reason a reference change can still matter: "
      "the cancellation is large, not perfect.")

raw_demo, _ = l2.load_erpcore("P3", DONE[0], resample_hz=RESAMPLE_HZ)
raw_demo.filter(L_FREQ, H_FREQ, picks=["eeg", "eog"], verbose=False).pick("eeg")

truth = raw_demo.copy().set_eeg_reference("average", verbose=False)               # the answer
cz_ref = raw_demo.copy().set_eeg_reference(["Cz"], verbose=False)                 # recorded against Cz
shipped = cz_ref.copy().drop_channels(["Cz"])                                     # what a dataset ships
naive = shipped.copy().set_eeg_reference("average", verbose=False)                # 29 channels, no reconstruction
with warnings.catch_warnings():                    # the new row has no position until set_montage runs
    warnings.simplefilter("ignore")
    fixed = mne.add_reference_channels(shipped.copy(), "Cz", copy=True)           # add the zero row back
    fixed.set_montage("standard_1005", on_missing="warn", match_case=True, verbose=False)
fixed.set_eeg_reference("average", verbose=False)

common = [c for c in truth.ch_names if c in naive.ch_names]
def rms_diff(a, b, picks):
    return float(np.sqrt(((a.get_data(picks=picks) - b.get_data(picks=picks)) ** 2).mean()) * 1e6)

print(f"Cz is exactly {np.abs(cz_ref.get_data(picks=['Cz'])).max() * 1e6:.3g} uV after re-referencing to it "
      "-- which is why datasets drop it.")
print(f"average over 29 channels, reference never reconstructed:  {rms_diff(truth, naive, common):6.3f} uV RMS "
      f"away from the true average reference")
print(f"Cz added back as zeros, then averaged over 30:            {rms_diff(truth, fixed, common):6.3f} uV RMS "
      f"away from the true average reference")
print(f"for scale, the channels themselves are {np.sqrt((truth.get_data(picks=common) ** 2).mean()) * 1e6:.1f} uV RMS")
print()
print(f"rank arithmetic: {l2.data_rank(30, n_reconstructed_reference=1, average_reference=True)['arithmetic']} "
      "-- the reconstructed channel adds a row without adding rank, and the average reference costs one more.")

lemon = helpers.load_lemon_raw("sub-010002", max_minutes=2.0, preload=True)
if lemon is None:
    print("ds-lemon could not be fetched; the arithmetic above is the same operation. "
          "TODO(confirm): re-run this section when the mirror is reachable.")
else:
    lemon.filter(1.0, 45.0, picks=["eeg", "eog"], verbose=False)
    eeg_before = [lemon.ch_names[i] for i in mne.pick_types(lemon.info, eeg=True)]
    print(f"ds-lemon sub-010002, first {lemon.times[-1] / 60:.1f} min: {len(eeg_before)} EEG channels at "
          f"{lemon.info['sfreq']:g} Hz; 'FCz' in the file: {'FCz' in lemon.ch_names} "
          "(the catalog documents FCz as the online reference)")

    naive_l = lemon.copy().set_eeg_reference("average", verbose=False)
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        fixed_l = mne.add_reference_channels(lemon.copy(), "FCz", copy=True)
        fixed_l.set_montage("standard_1005", on_missing="warn", match_case=True, verbose=False)
    fixed_l.set_eeg_reference("average", verbose=False)
    print(f"after reconstruction: {len(mne.pick_types(fixed_l.info, eeg=True))} EEG channels; "
          f"FCz is {np.abs(lemon.copy().get_data(picks='eeg')).max() * 0 + 0:.0f} uV before referencing "
          "(a row of zeros) and carries real signal afterwards, because the average is subtracted from it too")
    d = rms_diff(naive_l, fixed_l, eeg_before)
    print(f"difference between the two average references at the shared channels: {d:.3f} uV RMS, against "
          f"{np.sqrt((naive_l.get_data(picks=eeg_before) ** 2).mean()) * 1e6:.1f} uV RMS of signal "
          f"({100 * d / (np.sqrt((naive_l.get_data(picks=eeg_before) ** 2).mean()) * 1e6):.1f} %)")
    print(f"rank: {l2.data_rank(len(eeg_before) + 1, n_reconstructed_reference=1, average_reference=True)['arithmetic']}")

    fig, axes = plt.subplots(1, 3, figsize=(13, 3.9))
    band = (8.0, 12.0)
    for ax, (r, ttl) in zip(axes[:2], ((naive_l, f"average over {len(eeg_before)} channels, FCz never restored"),
                                       (fixed_l, f"FCz reconstructed, average over {len(eeg_before) + 1}"))):
        helpers.plot_band_topomap(r, band, ax=ax, title=f"{ttl}\n{band[0]:g}-{band[1]:g} Hz power")
    x0 = naive_l.get_data(picks=["Cz"])[0] * 1e6
    x1 = fixed_l.get_data(picks=["Cz"])[0] * 1e6
    t = naive_l.times
    sl = (t >= 20) & (t < 25)
    axes[2].plot(t[sl], x0[sl], color="tab:red", lw=0.8, label="FCz never restored")
    axes[2].plot(t[sl], x1[sl], color="tab:green", lw=0.8, ls="--", label="FCz reconstructed")
    axes[2].set(xlabel="Time (s)", ylabel="Amplitude (uV)", title="Cz under the two average references (uV)")
    axes[2].legend(fontsize=8); axes[2].grid(alpha=0.3)
    fig.tight_layout()
    plt.show()   # render the static figure(s) of this cell inline

by_ref = {n: np.array([amp[n][s] for s in DONE]) for n in names}
print("nb-2-3-reference -- L2.3 answer key (draft; TODO(confirm) at author review)")
print(f"Data: ds-erpcore P3, {DONE[0]}..{DONE[-1]} ({len(DONE)} subjects of 40; CC BY-SA 4.0 per data/directory.yaml, contested at source, open, per-subject downloadable). 30 EEG + 3 EOG, 1024 Hz resampled to {RESAMPLE_HZ:g} Hz, "
      f"CMS online reference, 60 Hz mains, no software filters.")
print(f"Pipeline: pyprep bad channels (>=2 criteria, cap 10 % of the montage) -> FIR zero-phase "
      f"{L_FREQ:g}-{H_FREQ:g} Hz on the continuous data -> interpolate -> re-reference -> epochs "
      f"{TMIN:g}..{TMAX:g} s, baseline {BASELINE[0]:g}..{BASELINE[1]:g} s. Identical trials under every "
      f"reference: autoreject's global peak-to-peak threshold on the average-referenced epochs, pooled and "
      f"condition-blind. Seed {l2.SEED}.")
print(f"Measure: mean amplitude of the target-minus-standard difference wave at {CH}, "
      f"{WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (a-priori window).")
print()
print(f"P3 mean amplitude at {CH} under each reference, grand mean over {len(DONE)} subjects:")
for n in names:
    v = by_ref[n]
    print(f"  {n:28s} {v.mean():+6.2f} uV   (SD {v.std(ddof=1):.2f}, SEM {v.std(ddof=1) / np.sqrt(len(v)):.2f}, "
          f"median {np.median(v):+.2f}, range {v.min():+.2f} to {v.max():+.2f}, "
          f"positive in {int((v > 0).sum())}/{len(v)} subjects)")
print()
print("The three keys L2.3's exercise asks for (tolerance: the SEM above is the natural choice; a learner who "
      "runs a different subset should land inside it):")
print(f"  ex-2-3-p3-average   P3 mean amplitude at {CH}, average reference          = "
      f"{by_ref['average (30 ch)'].mean():+.2f} uV  (SEM {by_ref['average (30 ch)'].std(ddof=1) / np.sqrt(len(DONE)):.2f})")
print(f"  ex-2-3-p3-mastoids  P3 mean amplitude at {CH}, linked-mastoid reference   = "
      f"{by_ref['linked mastoids (P9+P10)/2'].mean():+.2f} uV  "
      f"(SEM {by_ref['linked mastoids (P9+P10)/2'].std(ddof=1) / np.sqrt(len(DONE)):.2f})  "
      "[ERP CORE has no mastoid electrodes: (P9+P10)/2 is the equivalent, TODO(confirm)]")
print(f"  ex-2-3-p3-cz        P3 mean amplitude at {CH}, single-electrode Cz        = "
      f"{by_ref['Cz'].mean():+.2f} uV  (SEM {by_ref['Cz'].std(ddof=1) / np.sqrt(len(DONE)):.2f})")
print(f"  (for completeness, as recorded against CMS                                = "
      f"{by_ref['original (CMS)'].mean():+.2f} uV)")
print()
print("Ordering, which is what the free-response exercise is really about: "
      + " > ".join(n.split(' ')[0] for n in sorted(names, key=lambda n: -by_ref[n].mean())) + ".")
print(f"The inferior-posterior pair gives the largest number and Cz the smallest of the offline references, "
      f"which is what the lesson predicts: a reference site that carries the opposite polarity of the "
      f"component enlarges the difference, and a site near the component's own maximum shrinks it.")
print()
print(f"Topography check: the four maps correlate at "
      f"{min(float(np.corrcoef(maps[a] - maps[a].mean(), maps[b] - maps[b].mean())[0, 1]) for a in names for b in names):.6f} "
      f"or better after removing each map's spatial mean (the gradient is untouched), while the number of "
      f"electrodes above zero runs "
      + ", ".join(f"{n.split(' ')[0]} {int((maps[n] > 0).sum())}/{len(maps[n])}" for n in names) + ".")
print(f"Reference dependence: the target condition alone moves {t_rng:.2f} uV across the four references "
      f"(median over subjects) against {d_rng:.2f} uV for the difference wave.")
print()
print("Per-subject values (uV), for the tolerance discussion:")
print(l2.fmt_table(per_subject, ["subject", "interpolated", "rank", "n_kept", "drops driven by", *names],
                   floatfmt="{:+.2f}"))
print()
print("Caveat to carry into L2.5 and L2.6: this pipeline has no artifact correction, so the trial loss above is "
      "mostly blinks at Fp1/Fp2. The reference comparison is unaffected (identical trials throughout), but the "
      "spread across subjects would be smaller after ICA.")