# 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')
_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"]
    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 scipy import signal

SUBJECTS = ["sub-001", "sub-002", "sub-003"]
RESAMPLE_HZ, LP_HZ = 256.0, 30.0
CUTOFFS = [0.01, 0.1, 0.3, 0.5, 1.0, 2.0]
TMIN, TMAX, BASELINE = -0.2, 0.8, (-0.2, 0.0)
CH, WINDOW = l2.P3_CHANNEL, l2.P3_WINDOW
PRE_WINDOW, LATE_WINDOW = (-0.2, 0.0), (0.6, 0.8)
SELECT_PTP_UV, SELECT_HP_HZ = 150.0, 1.0

print(f"cutoffs (Hz): {CUTOFFS}; common low-pass {LP_HZ:g} Hz; epochs {TMIN:g}..{TMAX:g} s, "
      f"baseline {BASELINE[0]:g}..{BASELINE[1]:g} s")
print(f"measure: target minus standard at {CH}, mean over {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms; "
      f"also {PRE_WINDOW[0] * 1000:.0f}..{PRE_WINDOW[1] * 1000:.0f} ms and "
      f"{LATE_WINDOW[0] * 1000:.0f}..{LATE_WINDOW[1] * 1000:.0f} ms")
print(f"trial selection: peak-to-peak at {CH} <= {SELECT_PTP_UV:g} uV on a {SELECT_HP_HZ:g}-{LP_HZ:g} Hz copy, "
      f"once, for every cutoff")

# What each cutoff costs in filter length, before a single epoch is cut.
rows = []
for lf in CUTOFFS:
    h = mne.filter.create_filter(None, RESAMPLE_HZ, l_freq=lf, h_freq=None, verbose=False)
    rows.append({"high-pass (Hz)": lf, "FIR taps": len(h), "impulse response (s)": len(h) / RESAMPLE_HZ,
                 "half length (s)": (len(h) - 1) / 2 / RESAMPLE_HZ,
                 "transition band (Hz)": min(max(0.25 * lf, 2.0), lf)})
print()
print(l2.fmt_table(rows, ["high-pass (Hz)", "FIR taps", "impulse response (s)", "half length (s)",
                          "transition band (Hz)"], floatfmt="{:.2f}"))
print("MNE's default transition bandwidth for a high-pass is min(max(0.25 x cutoff, 2), cutoff) Hz, so below "
      "2 Hz the transition equals the cutoff and the Hamming-window length scales with its inverse: a 0.01 Hz "
      "high-pass has an impulse response longer than a third of this recording.")

t_start = time.time()
results, per_subject = {}, []
for sid in SUBJECTS:
    raw, _ = l2.load_erpcore("P3", sid, resample_hz=RESAMPLE_HZ)
    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.info["bads"] = list(dec["bads"])
    if dec["bads"]:
        raw.interpolate_bads(reset_bads=True, verbose=False)
    raw_lp = raw.copy().filter(None, LP_HZ, picks=["eeg", "eog"], verbose=False)     # the common low-pass
    raw_lp.set_eeg_reference("average", verbose=False)
    events, _, _ = l2.p3_events(raw_lp)

    sel = raw_lp.copy().filter(SELECT_HP_HZ, None, picks=["eeg", "eog"], verbose=False)
    sel_ep = l2.epochs_p3(sel, events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)
    keep = np.ptp(sel_ep.get_data(picks=[CH]), axis=2)[:, 0] * 1e6 <= SELECT_PTP_UV
    row = {"subject": sid, "interpolated": dec["bads"] or ["-"], "n_kept": int(keep.sum()),
           "n_trials": int(len(keep)),
           "n_target": int((events[keep, 2] == l2.P3_EVENT_ID["target"]).sum())}
    for lf in CUTOFFS:
        r = raw_lp.copy().filter(lf, None, picks=["eeg", "eog"], verbose=False)       # MNE defaults
        ep = l2.epochs_p3(r, events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)[keep]
        d = l2.difference_wave(ep)
        results.setdefault(lf, {})[sid] = {"difference": d, "target": ep["target"].average(),
                                           "standard": ep["standard"].average()}
        row[f"{lf:g} Hz"] = l2.mean_amplitude(d, CH, WINDOW)
    per_subject.append(row)
    print(f"  {sid}: interpolated {dec['bads'] or 'none'}, {row['n_kept']}/{row['n_trials']} trials kept "
          f"({row['n_target']} targets)  |  " + "  ".join(f"{lf:g} Hz {row[f'{lf:g} Hz']:+5.2f}" for lf in CUTOFFS))
DONE = [r["subject"] for r in per_subject]
print(f"\n{len(DONE)} subjects x {len(CUTOFFS)} cutoffs in {time.time() - t_start:.0f} s")
print()
print(l2.fmt_table(per_subject, ["subject", "interpolated", "n_kept", *[f"{lf:g} Hz" for lf in CUTOFFS]],
                   floatfmt="{:+.2f}"))

GA = {lf: {k: mne.grand_average([results[lf][s][k] for s in DONE]) for k in ("difference", "target", "standard")}
      for lf in CUTOFFS}

summary = []
for lf in CUTOFFS:
    d = GA[lf]["difference"]
    summary.append({
        "high-pass (Hz)": lf,
        f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)": l2.mean_amplitude(d, CH, WINDOW),
        f"pre-stimulus {PRE_WINDOW[0] * 1000:.0f}..{PRE_WINDOW[1] * 1000:.0f} ms (uV)": l2.mean_amplitude(d, CH, PRE_WINDOW),
        f"late {LATE_WINDOW[0] * 1000:.0f}-{LATE_WINDOW[1] * 1000:.0f} ms (uV)": l2.mean_amplitude(d, CH, LATE_WINDOW),
        "peak (uV)": float(d.data[d.ch_names.index(CH)].max() * 1e6),
        "peak latency (ms)": float(d.times[d.data[d.ch_names.index(CH)].argmax()] * 1000),
    })
base = summary[0][f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)"]
for s in summary:
    s["% of the 0.01 Hz value"] = 100 * s[f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)"] / base
print(f"grand average over {len(DONE)} subjects, target minus standard at {CH} (uV):\n")
print(l2.fmt_table(summary, list(summary[0]), floatfmt="{:+.2f}"))

fig, axes = plt.subplots(1, 2, figsize=(14, 4.6))
colors = plt.cm.viridis(np.linspace(0, 0.9, len(CUTOFFS)))
for lf, c in zip(CUTOFFS, colors):
    d = GA[lf]["difference"]
    axes[0].plot(d.times * 1000, d.data[d.ch_names.index(CH)] * 1e6, lw=1.4, color=c, label=f"{lf:g} Hz")
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"Target minus standard at {CH} for six high-pass cutoffs "
                  f"(uV, positive up; common {LP_HZ:g} Hz low-pass, n = {len(DONE)})")
axes[0].grid(alpha=0.3); axes[0].legend(fontsize=8, title="high-pass", title_fontsize=8)

vals = [s[f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)"] for s in summary]
late = [s[f"late {LATE_WINDOW[0] * 1000:.0f}-{LATE_WINDOW[1] * 1000:.0f} ms (uV)"] for s in summary]
axes[1].plot(CUTOFFS, vals, "o-", color="tab:blue", label=f"{WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms mean")
axes[1].plot(CUTOFFS, late, "s--", color="tab:red",
             label=f"{LATE_WINDOW[0] * 1000:.0f}-{LATE_WINDOW[1] * 1000:.0f} ms mean (the opposite-polarity lobe)")
axes[1].axhline(0, color="gray", lw=0.6)
axes[1].set(xscale="log", xlabel="High-pass cutoff (Hz) [log]", ylabel="Amplitude (uV)",
            title=f"P3 mean amplitude against high-pass cutoff (uV, {CH}, target minus standard)")
axes[1].grid(alpha=0.3, which="both"); axes[1].legend(fontsize=8)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

SID = DONE[0]
raw, _ = l2.load_erpcore("P3", SID, resample_hz=RESAMPLE_HZ)
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.info["bads"] = list(dec["bads"])
if dec["bads"]:
    raw.interpolate_bads(reset_bads=True, verbose=False)
rank = l2.data_rank(n_eeg, n_interpolated=len(dec["bads"]), average_reference=True)["rank"]

analysis_01 = raw.copy().filter(0.1, LP_HZ, picks=["eeg", "eog"], verbose=False).set_eeg_reference("average", verbose=False)
analysis_1 = raw.copy().filter(1.0, LP_HZ, picks=["eeg", "eog"], verbose=False).set_eeg_reference("average", verbose=False)
# The fit copy is built from the same *unfiltered* recording, not from the 0.1-30 Hz analysis branch:
# high-passing data that has already been low-passed at 30 Hz leaves the fit with no content above 30 Hz,
# so muscle and line-noise components cannot be separated and ICLabel is used outside its documented regime.
ica_src = raw.copy().filter(1.0, 100.0, picks=["eeg", "eog"], verbose=False).set_eeg_reference("average", verbose=False)

t0 = time.time()
ica, raw_for_ica, ica_log = l2.two_pass_ica(analysis_01, n_components=rank, fit_on=ica_src,
                                            fit_l_freq=1.0, fit_h_freq=100.0, seed=l2.SEED)
cls = l2.classify_components(raw_for_ica, ica)
drop = [i for i, (lab, p) in enumerate(zip(cls["labels"], cls["probabilities"]))
        if lab in ("eye", "muscle", "heart", "line noise", "channel noise") and p >= 0.80][:6]
ica.exclude = sorted(drop)
print(f"{SID}: rank {rank}, ICA ({ica_log['method']}, seed {ica_log['seed']}) fitted on a "
      f"{ica_log['fit_filter_hz'][0]:g}-{ica_log['fit_lowpass_hz']:g} Hz branch of the same recording in "
      f"{ica_log['duration_s']:.0f} s; classifier: {cls['tool']}")
print(f"  components removed: {ica.exclude} -> "
      + ", ".join(f"IC{i} {cls['labels'][i]} (p {cls['probabilities'][i]:.2f})" for i in ica.exclude))

events, _, _ = l2.p3_events(analysis_01)
variants = {}
for name, src in (("two-pass: fit at 1 Hz, apply to the 0.1 Hz data", analysis_01),
                  ("shortcut: fit and analyse at 1 Hz", analysis_1)):
    cleaned = ica.apply(src.copy(), verbose=False)
    ep = l2.epochs_p3(cleaned, events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)
    ptp = l2.epoch_ptp_uv(ep)
    variants[name] = {"epochs": ep, "ptp": ptp}
keep2 = (variants["two-pass: fit at 1 Hz, apply to the 0.1 Hz data"]["ptp"] <= SELECT_PTP_UV)
rows = []
for name, v in variants.items():
    d = l2.difference_wave(v["epochs"][keep2])
    rows.append({"pipeline": name, f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)":
                 l2.mean_amplitude(d, CH, WINDOW),
                 f"late {LATE_WINDOW[0] * 1000:.0f}-{LATE_WINDOW[1] * 1000:.0f} ms (uV)":
                 l2.mean_amplitude(d, CH, LATE_WINDOW),
                 "trials": int(keep2.sum())})
    variants[name]["difference"] = d
print()
print(l2.fmt_table(rows, list(rows[0]), floatfmt="{:+.2f}"))
a = rows[0][f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)"]
b = rows[1][f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)"]
print(f"\nSame decomposition, same removed components, same trials: the only difference is which data the "
      f"unmixing was applied to. Analysing the 1 Hz data instead of the 0.1 Hz data changes the measured P3 "
      f"from {a:+.2f} to {b:+.2f} uV ({100 * (b - a) / abs(a):+.0f} %). The extra copy is the whole cost of "
      "avoiding that.")

HP_DEMO = 0.5
h_demo = mne.filter.create_filter(None, RESAMPLE_HZ, l_freq=HP_DEMO, h_freq=None, verbose=False)
epoch_len_s = TMAX - TMIN
print(f"{HP_DEMO:g} Hz high-pass at {RESAMPLE_HZ:g} Hz: {len(h_demo)} taps = {len(h_demo) / RESAMPLE_HZ:.1f} s "
      f"of impulse response, against an epoch of {epoch_len_s:.1f} s "
      f"({len(h_demo) / RESAMPLE_HZ / epoch_len_s:.0f} times longer than the epoch)")

cont = raw_lp_demo = raw.copy().filter(None, LP_HZ, picks=["eeg", "eog"], verbose=False)
cont.set_eeg_reference("average", verbose=False)
ev_demo, _, _ = l2.p3_events(cont)

ep_raw = l2.epochs_p3(cont, ev_demo, tmin=TMIN, tmax=TMAX, baseline=None)          # no baseline yet
before = l2.epochs_p3(cont.copy().filter(HP_DEMO, None, picks=["eeg", "eog"], verbose=False), ev_demo,
                      tmin=TMIN, tmax=TMAX, baseline=None)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")                     # MNE warns that the filter is longer than the epoch
    after = ep_raw.copy().filter(HP_DEMO, None, picks="eeg", verbose=False)
before.apply_baseline(BASELINE, verbose=False)
after.apply_baseline(BASELINE, verbose=False)

d_before = l2.difference_wave(before[keep2])
d_after = l2.difference_wave(after[keep2])
i = d_before.ch_names.index(CH)
print(f"\ntarget minus standard at {CH}, {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms:")
print(f"  filtered continuous, then epoched: {l2.mean_amplitude(d_before, CH, WINDOW):+.2f} uV")
print(f"  epoched, then filtered:            {l2.mean_amplitude(d_after, CH, WINDOW):+.2f} uV")
print(f"  RMS difference between the two waveforms over the whole epoch: "
      f"{np.sqrt(((d_before.data[i] - d_after.data[i]) ** 2).mean()) * 1e6:.2f} uV")

fig, axes = plt.subplots(1, 2, figsize=(13, 4.2))
axes[0].plot(d_before.times * 1000, d_before.data[i] * 1e6, color="tab:blue", lw=1.4,
             label="filter the continuous data, then epoch")
axes[0].plot(d_after.times * 1000, d_after.data[i] * 1e6, color="tab:red", lw=1.4, ls="--",
             label="epoch, then filter each epoch")
axes[0].axvspan(WINDOW[0] * 1000, WINDOW[1] * 1000, color="tab:orange", alpha=0.18)
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"Same {HP_DEMO:g} Hz high-pass, two placements ({CH}, uV, positive up)")
axes[0].grid(alpha=0.3); axes[0].legend(fontsize=8)

x_b = before.get_data(picks=[CH])[:, 0, :] * 1e6
x_a = after.get_data(picks=[CH])[:, 0, :] * 1e6
axes[1].plot(before.times * 1000, x_b.std(axis=0), color="tab:blue", lw=1.4, label="continuous, then epoched")
axes[1].plot(after.times * 1000, x_a.std(axis=0), color="tab:red", lw=1.4, ls="--", label="epoched, then filtered")
axes[1].set(xlabel="Time from stimulus (ms)", ylabel="Across-trial standard deviation (uV)",
            title=f"Where the variance sits inside the epoch ({CH}, uV)")
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

CUT_A, GAP_S = 150.0, 20.0
src = raw.copy().pick("eeg").filter(None, LP_HZ, verbose=False)
first = src.copy().crop(0, CUT_A)
second = src.copy().crop(CUT_A + GAP_S, min(CUT_A + GAP_S + 120.0, src.times[-1]))
joined = mne.concatenate_raws([first, second], verbose=False)
boundary = [(float(a["onset"]), str(a["description"])) for a in joined.annotations
            if "boundary" in str(a["description"]).lower()]
print(f"joined {first.times[-1]:.0f} s + {second.times[-1]:.0f} s with {GAP_S:g} s removed between them; "
      f"annotations inserted by concatenate_raws: {boundary}")

HP_BOUND = 0.5
respected = joined.copy().filter(HP_BOUND, None, verbose=False)                       # MNE's default
ignored = joined.copy().filter(HP_BOUND, None, skip_by_annotation=(), verbose=False)  # boundary information lost
t = joined.times
x_r = respected.get_data(picks=[CH])[0] * 1e6
x_i = ignored.get_data(picks=[CH])[0] * 1e6
diff = np.abs(x_r - x_i)
thr = 0.05 * float(np.std(x_r))
b_t = float(boundary[0][0]) if boundary else CUT_A
near = np.where(diff > thr)[0]
if len(near):
    ring_s = float(t[near[-1]] - t[near[0]])
    lo, hi = float(t[near[0]] - b_t), float(t[near[-1]] - b_t)
else:
    ring_s, lo, hi = 0.0, 0.0, 0.0
h_b = mne.filter.create_filter(None, RESAMPLE_HZ, l_freq=HP_BOUND, h_freq=None, verbose=False)
print(f"\nboundary at {b_t:.1f} s; {HP_BOUND:g} Hz high-pass, {len(h_b)} taps = {len(h_b) / RESAMPLE_HZ:.1f} s "
      f"(half length {(len(h_b) - 1) / 2 / RESAMPLE_HZ:.1f} s)")
print(f"the two outputs differ by more than 5 % of the signal SD ({thr:.2f} uV) from "
      f"{lo:+.1f} s to {hi:+.1f} s around the join: a ringing duration of {ring_s:.1f} s")
print(f"largest excursion introduced by ignoring the boundary: {diff.max():.1f} uV at "
      f"{t[int(diff.argmax())] - b_t:+.2f} s from the join, against {np.std(x_r):.1f} uV of signal SD")
print(f"how that compares with the filter: the ringing spans {ring_s / (len(h_b) / RESAMPLE_HZ):.2f} of one "
      "impulse response, which is what a discontinuity convolved with a long symmetric kernel should produce.")

fig, axes = plt.subplots(2, 1, figsize=(12, 6.5), sharex=True)
w = (t > b_t - 25) & (t < b_t + 25)
axes[0].plot(t[w], x_r[w], color="tab:blue", lw=0.8, label="boundary respected (MNE default)")
axes[0].plot(t[w], x_i[w], color="tab:red", lw=0.8, label="boundary ignored (skip_by_annotation=())")
axes[0].axvline(b_t, color="k", lw=1.0, ls=":", label="join")
axes[0].set(ylabel="Amplitude (uV)",
            title=f"{CH} around a discontinuity, {HP_BOUND:g} Hz high-pass (uV)")
axes[0].legend(fontsize=8); axes[0].grid(alpha=0.3)
axes[1].plot(t[w], diff[w], color="k", lw=0.8)
axes[1].axhline(thr, color="tab:orange", lw=0.9, label=f"5 % of the signal SD ({thr:.2f} uV)")
axes[1].axvline(b_t, color="k", lw=1.0, ls=":")
axes[1].set(xlabel="Time (s)", ylabel="|difference| (uV)",
            title=f"What ignoring the boundary costs (uV); ringing lasts {ring_s:.1f} s")
axes[1].legend(fontsize=8); axes[1].grid(alpha=0.3)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

print("nb-2-4-filter-choices -- L2.4 numbers (draft; TODO(confirm) at author review)")
print(f"Data: ds-erpcore P3, {', '.join(DONE)} ({len(DONE)} subjects; CC BY-SA 4.0 per data/directory.yaml, contested at source, open, "
      f"per-subject downloadable). 30 EEG + 3 EOG, 1024 Hz resampled to {RESAMPLE_HZ:g} Hz, CMS online "
      f"reference re-referenced to the average of the 30 EEG channels, 60 Hz mains, no software filters.")
print(f"Held fixed: common {LP_HZ:g} Hz FIR low-pass; pyprep bad channels (>=2 criteria, cap 10 %) and "
      f"interpolation; epochs {TMIN:g}..{TMAX:g} s, baseline {BASELINE[0]:g}..{BASELINE[1]:g} s; one trial set "
      f"for every cutoff (peak-to-peak at {CH} <= {SELECT_PTP_UV:g} uV on a {SELECT_HP_HZ:g}-{LP_HZ:g} Hz copy).")
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("Amplitude per high-pass cutoff (grand average over the three subjects):")
key = f"{CH} {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms (uV)"
late_key = f"late {LATE_WINDOW[0] * 1000:.0f}-{LATE_WINDOW[1] * 1000:.0f} ms (uV)"
for s in summary:
    print(f"  high-pass {s['high-pass (Hz)']:5.2f} Hz: {CH} mean amplitude = {s[key]:+5.2f} uV "
          f"({s['% of the 0.01 Hz value']:5.1f} % of the 0.01 Hz value); "
          f"late-window mean {s[late_key]:+5.2f} uV; peak {s['peak (uV)']:+5.2f} uV at "
          f"{s['peak latency (ms)']:.0f} ms")
print()
print("Per subject (uV):")
print(l2.fmt_table(per_subject, ["subject", "n_kept", *[f"{lf:g} Hz" for lf in CUTOFFS]], floatfmt="{:+.2f}"))
print()
v0, v1 = summary[0][key], summary[-1][key]
print(f"Trend: the P3 mean amplitude falls from {v0:+.2f} uV at 0.01 Hz to {v1:+.2f} uV at "
      f"{CUTOFFS[-1]:g} Hz ({100 * v1 / v0:.0f} % of the value at the lowest cutoff), and the "
      f"{LATE_WINDOW[0] * 1000:.0f}-{LATE_WINDOW[1] * 1000:.0f} ms window moves from "
      f"{summary[0][late_key]:+.2f} to {summary[-1][late_key]:+.2f} uV -- the opposite-polarity lobe a "
      "zero-phase filter creates out of the energy it removed (pf-hp-cutoff-erp).")
print()
print("Two-pass strategy (single subject " + SID + f", rank {rank}, {len(ica.exclude)} ICA components removed, "
      f"same decomposition and same trials in both rows):")
for r in rows:
    print(f"  {r['pipeline']:52s} {CH} mean amplitude = {r[key]:+5.2f} uV")
print(f"  cost of the shortcut: {b - a:+.2f} uV ({100 * (b - a) / abs(a):+.0f} %)")
print()
print("Filtering epoched data (pf-filter-epoched-data):")
print(f"  {HP_DEMO:g} Hz high-pass = {len(h_demo)} taps = {len(h_demo) / RESAMPLE_HZ:.1f} s of impulse response, "
      f"on a {epoch_len_s:.1f}-s epoch")
print(f"  continuous then epoched: {l2.mean_amplitude(d_before, CH, WINDOW):+.2f} uV; "
      f"epoched then filtered: {l2.mean_amplitude(d_after, CH, WINDOW):+.2f} uV; "
      f"RMS difference over the epoch {np.sqrt(((d_before.data[i] - d_after.data[i]) ** 2).mean()) * 1e6:.2f} uV")
print()
print("Boundary artifact (pf-filter-across-boundaries):")
print(f"  {GAP_S:g} s removed from the middle of the recording and the halves concatenated; "
      f"{HP_BOUND:g} Hz high-pass ({len(h_b)} taps = {len(h_b) / RESAMPLE_HZ:.1f} s)")
print(f"  ringing duration with the boundary ignored: {ring_s:.1f} s "
      f"(from {lo:+.1f} s to {hi:+.1f} s around the join), largest excursion {diff.max():.1f} uV against "
      f"{np.std(x_r):.1f} uV of signal SD")
print(f"  with MNE's default skip_by_annotation the two spans are filtered independently and the artifact "
      "does not occur at all.")
print()
print("L2.4's exercises are multiple-choice (ex-2-4-settings-erp / -ica / -oscillation) plus a free response; "
      "this notebook supplies the evidence behind the keys rather than a numeric answer.")