# 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 AutoReject, 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
THRESHOLDS = np.arange(20.0, 305.0, 5.0)
GAP_PP = 20.0
FIXED_EXAMPLES = [75.0, 100.0, 150.0, 200.0]

print(f"subset {SUBJECTS[0]} .. {SUBJECTS[-1]} ({len(SUBJECTS)} of 40 participants)")
print(f"statistic: peak-to-peak within the epoch, maximised over all 30 EEG channels")
print(f"sweep: {THRESHOLDS[0]:g} to {THRESHOLDS[-1]:g} uV in {THRESHOLDS[1] - THRESHOLDS[0]:g} uV steps "
      f"({len(THRESHOLDS)} thresholds); the key is the largest threshold with a per-condition gap of at least "
      f"{GAP_PP:g} percentage points")
print(f"measure: target minus standard at {CH}, mean over {WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms, "
      f"and (for contrast) the peak in the same window")
print(f"seed {l2.SEED}")

t_start = time.time()
store, skipped = {}, []
for sid in SUBJECTS:
    try:
        raw, _ = l2.load_erpcore("P3", sid, resample_hz=RESAMPLE_HZ)
    except Exception as e:
        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)
    raw.set_eeg_reference("average", verbose=False)
    events, _, ev_info = l2.p3_events(raw)
    ep = l2.epochs_p3(raw, events, tmin=TMIN, tmax=TMAX, baseline=BASELINE)
    ptp = l2.epoch_ptp_uv(ep)
    is_target = ep.events[:, 2] == l2.P3_EVENT_ID["target"]
    store[sid] = {"epochs": ep, "ptp": ptp, "is_target": is_target, "bads": dec["bads"],
                  "raw": raw, "events": events,
                  "rank": l2.data_rank(n_eeg, n_interpolated=len(dec["bads"]), average_reference=True)["rank"]}
    print(f"  {sid}: {len(ep)} epochs ({int(is_target.sum())} target / {int((~is_target).sum())} standard), "
          f"interpolated {dec['bads'] or 'none'}, peak-to-peak median {np.median(ptp):5.0f} uV, "
          f"95th percentile {np.percentile(ptp, 95):5.0f} uV, max {ptp.max():6.0f} uV")
DONE = list(store)
print(f"\n{len(DONE)} subjects in {time.time() - t_start:.0f} s; skipped: {skipped or 'none'}")

sweeps, key_rows = {}, []
for sid in DONE:
    sw = l2.condition_bias_sweep(store[sid]["ptp"], store[sid]["is_target"], THRESHOLDS, gap_pp=GAP_PP)
    sweeps[sid] = sw
    k = sw["key"]
    key_rows.append({
        "subject": sid,
        "n target / standard": f"{int(store[sid]['is_target'].sum())} / {int((~store[sid]['is_target']).sum())}",
        "ptp median (uV)": float(np.median(store[sid]["ptp"])),
        f"largest threshold with a >= {GAP_PP:g} pp gap (uV)": k["threshold_uv"] if k else float("nan"),
        "target rejected (%)": k["percent_rejected_target"] if k else float("nan"),
        "standard rejected (%)": k["percent_rejected_standard"] if k else float("nan"),
        "gap (pp)": k["difference_pp"] if k else float("nan"),
        "which loses more": ("target" if k and k["difference_pp"] > 0 else "standard" if k else "-"),
    })
print(l2.fmt_table(key_rows, list(key_rows[0]), floatfmt="{:.1f}"))
with_key = [r for r in key_rows if not np.isnan(r[f"largest threshold with a >= {GAP_PP:g} pp gap (uV)"])]
print()
print(f"{len(with_key)} of {len(DONE)} subjects reach a {GAP_PP:g} percentage-point gap anywhere in the sweep.")

# The detail subject is chosen by a stated rule, not by eye: among the subjects that reach the gap, the one
# whose gap appears at the *largest* threshold -- the least total data loss for the clearest demonstration.
KEY_COL = f"largest threshold with a >= {GAP_PP:g} pp gap (uV)"
if with_key:
    DETAIL = max(with_key, key=lambda r: r[KEY_COL])["subject"]
    detail_key = sweeps[DETAIL]["key"]
else:                                    # no subject reaches the gap: report the largest gap that exists
    DETAIL = max(DONE, key=lambda s: max(abs(r["difference_pp"]) for r in sweeps[s]["rows"]))
    detail_key = max(sweeps[DETAIL]["rows"], key=lambda r: abs(r["difference_pp"]))
    print(f"no subject reaches a {GAP_PP:g} pp gap anywhere in the sweep; reporting the largest gap found "
          f"({DETAIL}, {detail_key['difference_pp']:+.1f} pp) and marking the key TODO(confirm).")
print(f"detail subject: {DETAIL} (rule: among subjects that reach a {GAP_PP:g} pp gap, the one whose gap "
      f"appears at the largest threshold)")
print(f"  key threshold {detail_key['threshold_uv']:.0f} uV -> target {detail_key['percent_rejected_target']:.1f} % "
      f"rejected, standard {detail_key['percent_rejected_standard']:.1f} % rejected, gap "
      f"{detail_key['difference_pp']:+.1f} pp, {detail_key['percent_rejected_all']:.1f} % of all trials lost")

fig, axes = plt.subplots(1, 2, figsize=(14, 4.6))
sw = sweeps[DETAIL]["rows"]
th = [r["threshold_uv"] for r in sw]
axes[0].plot(th, [r["percent_rejected_target"] for r in sw], color="tab:blue", lw=1.6, label="target (rare)")
axes[0].plot(th, [r["percent_rejected_standard"] for r in sw], color="tab:orange", lw=1.6, label="standard (frequent)")
axes[0].plot(th, [r["percent_rejected_all"] for r in sw], color="0.5", lw=1.0, ls="--", label="all trials")
axes[0].axvline(detail_key["threshold_uv"], color="tab:red", lw=1.0, ls=":",
                label=f"key: {detail_key['threshold_uv']:.0f} uV")
axes[0].set(xlabel="Peak-to-peak rejection threshold (uV)", ylabel="Trials rejected (%)",
            title=f"{DETAIL}: rejection rate per condition against threshold (%)")
axes[0].grid(alpha=0.3); axes[0].legend(fontsize=8)

for sid in DONE:
    rowsw = sweeps[sid]["rows"]
    axes[1].plot([r["threshold_uv"] for r in rowsw], [r["difference_pp"] for r in rowsw],
                 lw=1.6 if sid == DETAIL else 0.8, color="tab:red" if sid == DETAIL else "0.6",
                 label=DETAIL if sid == DETAIL else None)
axes[1].axhline(GAP_PP, color="k", lw=0.8, ls="--", label=f"+/- {GAP_PP:g} pp")
axes[1].axhline(-GAP_PP, color="k", lw=0.8, ls="--")
axes[1].axhline(0, color="gray", lw=0.6)
axes[1].set(xlabel="Peak-to-peak rejection threshold (uV)",
            ylabel="target % rejected - standard % rejected (percentage points)",
            title=f"Condition bias against threshold, all {len(DONE)} subjects (percentage points)")
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

ep_d = store[DETAIL]["epochs"]
ptp_d, is_t = store[DETAIL]["ptp"], store[DETAIL]["is_target"]
t0 = time.time()

thr_global = get_rejection_threshold(ep_d, random_state=l2.SEED, verbose=False)["eeg"] * 1e6
ar = AutoReject(n_interpolate=[1, 2, 4], consensus=np.linspace(0.3, 1.0, 5), random_state=l2.SEED,
                n_jobs=1, verbose=False)
ar.fit(ep_d)
reject_log = ar.get_reject_log(ep_d)
print(f"autoreject: global threshold {thr_global:.0f} uV; local fit in {time.time() - t0:.0f} s "
      f"(n_interpolate {list(ar.n_interpolate_.values())[:1]}, consensus "
      f"{[round(float(v), 2) for v in list(ar.consensus_.values())[:1]]} per channel-type)")

strategies = {}
for v in FIXED_EXAMPLES:
    strategies[f"fixed {v:.0f} uV"] = ptp_d > v
strategies[f"autoreject global ({thr_global:.0f} uV)"] = ptp_d > thr_global
strategies["autoreject local (repair or reject)"] = np.asarray(reject_log.bad_epochs, bool)
strategies["no rejection"] = np.zeros(len(ptp_d), bool)

rows = []
for name, rejected in strategies.items():
    tab = {r["condition"]: r for r in l2.rejection_table(ep_d, rejected)}
    kept = ep_d[~rejected] if rejected.any() else ep_d
    d = l2.difference_wave(kept)
    snr = l2.erp_snr(d, CH, WINDOW, BASELINE)
    i = d.ch_names.index(CH)
    m = (d.times >= WINDOW[0]) & (d.times <= WINDOW[1])
    rows.append({
        "strategy": name,
        "kept": tab["all"]["n_kept"],
        "target % rejected": tab["target"]["percent_rejected"],
        "standard % rejected": tab["standard"]["percent_rejected"],
        "gap (pp)": tab["target"]["percent_rejected"] - tab["standard"]["percent_rejected"],
        f"{CH} mean (uV)": l2.mean_amplitude(d, CH, WINDOW),
        f"{CH} peak (uV)": float(d.data[i, m].max() * 1e6),
        "baseline noise (uV)": snr["noise_uv"],
        "SNR": snr["snr"],
    })
print()
print(l2.fmt_table(rows, list(rows[0]), floatfmt="{:.2f}"))
print()
print("The per-condition columns are the point. A number chosen because it is round has no reason to treat the "
      "two conditions alike; autoreject's is estimated on pooled trials and so is condition-blind by "
      "construction -- which does not guarantee equal rejection, only that the condition played no part in "
      "setting the number.")

balanced = min(sweeps[DETAIL]["rows"],
               key=lambda r: (abs(r["difference_pp"]) + 100 * (r["percent_rejected_all"] > 25)))
print(f"balanced threshold (smallest per-condition gap with under 25 % total loss): "
      f"{balanced['threshold_uv']:.0f} uV -> gap {balanced['difference_pp']:+.1f} pp, "
      f"{balanced['percent_rejected_all']:.1f} % of trials lost")
print(f"biased threshold (the key):        {detail_key['threshold_uv']:.0f} uV -> gap "
      f"{detail_key['difference_pp']:+.1f} pp, {detail_key['percent_rejected_all']:.1f} % of trials lost")

variants = {
    "no rejection": np.zeros(len(ptp_d), bool),
    f"balanced ({balanced['threshold_uv']:.0f} uV)": ptp_d > balanced["threshold_uv"],
    f"biased ({detail_key['threshold_uv']:.0f} uV)": ptp_d > detail_key["threshold_uv"],
}
comparison, waves = [], {}
for name, rejected in variants.items():
    kept = ep_d[~rejected] if rejected.any() else ep_d
    d = l2.difference_wave(kept)
    waves[name] = d
    tab = {r["condition"]: r for r in l2.rejection_table(ep_d, rejected)}
    i = d.ch_names.index(CH)
    m = (d.times >= WINDOW[0]) & (d.times <= WINDOW[1])
    snr = l2.erp_snr(d, CH, WINDOW, BASELINE)
    comparison.append({"trial set": name,
                       "target kept": tab["target"]["n_kept"], "standard kept": tab["standard"]["n_kept"],
                       "gap (pp)": tab["target"]["percent_rejected"] - tab["standard"]["percent_rejected"],
                       f"{CH} mean (uV)": l2.mean_amplitude(d, CH, WINDOW),
                       f"{CH} peak (uV)": float(d.data[i, m].max() * 1e6),
                       "peak - mean (uV)": float(d.data[i, m].max() * 1e6) - l2.mean_amplitude(d, CH, WINDOW),
                       "baseline noise (uV)": snr["noise_uv"]})
print()
print(l2.fmt_table(comparison, list(comparison[0]), floatfmt="{:+.2f}"))

fig, axes = plt.subplots(1, 2, figsize=(14, 4.4))
for (name, d), color in zip(waves.items(), ("0.4", "tab:blue", "tab:red")):
    i = d.ch_names.index(CH)
    axes[0].plot(d.times * 1000, d.data[i] * 1e6, lw=1.5, color=color, label=f"{name}, n = {d.nave}")
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"{DETAIL}: target minus standard at {CH} under three trial sets (uV, positive up)")
axes[0].grid(alpha=0.3); axes[0].legend(fontsize=8)

bins = np.linspace(0, min(400, ptp_d.max()), 40)
axes[1].hist(ptp_d[is_t], bins=bins, alpha=0.6, color="tab:blue", label=f"target (n = {int(is_t.sum())})", density=True)
axes[1].hist(ptp_d[~is_t], bins=bins, alpha=0.6, color="tab:orange",
             label=f"standard (n = {int((~is_t).sum())})", density=True)
axes[1].axvline(detail_key["threshold_uv"], color="tab:red", lw=1.2, ls=":",
                label=f"key threshold {detail_key['threshold_uv']:.0f} uV")
axes[1].axvline(balanced["threshold_uv"], color="tab:blue", lw=1.2, ls=":",
                label=f"balanced {balanced['threshold_uv']:.0f} uV")
axes[1].set(xlabel="Epoch peak-to-peak, maximum over channels (uV)", ylabel="Density (1/uV)",
            title=f"{DETAIL}: where the two conditions' peak-to-peak distributions differ (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

raw_d = store[DETAIL]["raw"]                     # the 0.1-30 Hz analysis branch built in section 2
rank_d = store[DETAIL]["rank"]
bads_d = store[DETAIL]["bads"]

# The ICA fit copy must be built from data that has NOT been low-passed at 30 Hz: the two-pass split of
# L2.4 happens at the filter step, and both branches then get the same interpolation and the same reference.
raw_pre, _ = l2.load_erpcore("P3", DETAIL, resample_hz=RESAMPLE_HZ)
raw_ica_src = raw_pre.filter(1.0, 100.0, picks=["eeg", "eog"], verbose=False)
raw_ica_src.info["bads"] = list(bads_d)
if bads_d:
    raw_ica_src.interpolate_bads(reset_bads=True, verbose=False)
raw_ica_src.set_eeg_reference("average", verbose=False)
ica, raw_for_ica, ica_log = l2.two_pass_ica(raw_d, n_components=rank_d, fit_on=raw_ica_src, 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"{DETAIL}: rank {rank_d}, ICA fitted on a {ica_log['fit_filter_hz'][0]:g}-"
      f"{ica_log['fit_filter_hz'][1]:g} Hz branch in {ica_log['duration_s']:.0f} s ({cls['tool']}); removed "
      + (", ".join(f"IC{i} {cls['labels'][i]} (p {cls['probabilities'][i]:.2f})" for i in ica.exclude) or "nothing"))

clean = ica.apply(raw_d.copy(), verbose=False)
ep_clean = l2.epochs_p3(clean, store[DETAIL]["events"], tmin=TMIN, tmax=TMAX, baseline=BASELINE)
ptp_clean = l2.epoch_ptp_uv(ep_clean)
sw_clean = l2.condition_bias_sweep(ptp_clean, is_t, THRESHOLDS, gap_pp=GAP_PP)

print(f"\npeak-to-peak, maximum over channels: median {np.median(ptp_d):.0f} uV before ICA -> "
      f"{np.median(ptp_clean):.0f} uV after; 95th percentile {np.percentile(ptp_d, 95):.0f} -> "
      f"{np.percentile(ptp_clean, 95):.0f} uV")
for name, ptp_v, sw_v in (("before ICA", ptp_d, sweeps[DETAIL]), ("after ICA", ptp_clean, sw_clean)):
    k = sw_v["key"]
    at100 = next(r for r in sw_v["rows"] if r["threshold_uv"] == 100.0)
    print(f"  {name:10s}: at a 100 uV threshold {at100['percent_rejected_all']:5.1f} % of trials are rejected "
          f"(target {at100['percent_rejected_target']:5.1f} %, standard {at100['percent_rejected_standard']:5.1f} %); "
          f"largest threshold with a {GAP_PP:g} pp gap: "
          + (f"{k['threshold_uv']:.0f} uV" if k else "none in the sweep"))
d_clean = l2.difference_wave(ep_clean[ptp_clean <= 100.0])
d_dirty = l2.difference_wave(ep_d[ptp_d <= 100.0])
print(f"  at that 100 uV threshold the difference wave at {CH} measures "
      f"{l2.mean_amplitude(d_dirty, CH, WINDOW):+.2f} uV from {d_dirty.nave} trials before ICA and "
      f"{l2.mean_amplitude(d_clean, CH, WINDOW):+.2f} uV from {d_clean.nave} trials after")

print("nb-2-5-rejection -- L2.5 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.")
print(f"Pipeline (canonical order, stopping before ICA): pyprep bad channels (>=2 criteria, cap 10 %) -> FIR "
      f"zero-phase {L_FREQ:g}-{H_FREQ:g} Hz on the continuous data -> interpolate -> average reference -> "
      f"epochs {TMIN:g}..{TMAX:g} s, baseline {BASELINE[0]:g}..{BASELINE[1]:g} s. Seed {l2.SEED}.")
print(f"Statistic: peak-to-peak within the epoch, maximum over all 30 EEG channels. Sweep {THRESHOLDS[0]:g}-"
      f"{THRESHOLDS[-1]:g} uV in {THRESHOLDS[1] - THRESHOLDS[0]:g} uV steps.")
print()
print(f"ex-2-5-biased-threshold -- the threshold at which one condition loses at least {GAP_PP:g} percentage "
      f"points more of its trials than the other (largest such threshold; tolerance: the sweep's step size, "
      f"{THRESHOLDS[1] - THRESHOLDS[0]:g} uV):")
print(f"  ANSWER, detail subject {DETAIL}: {detail_key['threshold_uv']:.0f} uV")
print(f"    target {detail_key['percent_rejected_target']:.1f} % rejected, standard "
      f"{detail_key['percent_rejected_standard']:.1f} % rejected, gap {detail_key['difference_pp']:+.1f} pp, "
      f"{detail_key['percent_rejected_all']:.1f} % of all trials lost")
print(f"    ({DETAIL} was selected by the stated rule: among subjects reaching the gap, the one whose gap "
      f"appears at the largest threshold.)")
_dir = "target" if detail_key["difference_pp"] > 0 else "standard"
print(f"    Direction: on this subject it is the {_dir} condition that loses more. The lesson's worked story "
      "runs the other way -- the rare, attended condition provokes more blinks -- and both directions occur "
      f"in this subset ({', '.join(r['subject'] + ': ' + r['which loses more'] for r in with_key)}). The "
      "direction is not the finding; the existence of a gap the experimenter did not choose is.")
print()
print(f"Every subject in the subset, so the key can be matched to whichever subject w-threshold-tuner ships:")
print(l2.fmt_table(key_rows, list(key_rows[0]), floatfmt="{:.1f}"))
print()
print(f"What the biased threshold does to the difference wave ({DETAIL}, {CH}, "
      f"{WINDOW[0] * 1000:.0f}-{WINDOW[1] * 1000:.0f} ms):")
print(l2.fmt_table(comparison, list(comparison[0]), floatfmt="{:+.2f}"))
nb_no = comparison[0]
nb_bias = comparison[-1]
print(f"  mean amplitude moves {nb_no[f'{CH} mean (uV)']:+.2f} -> {nb_bias[f'{CH} mean (uV)']:+.2f} uV "
      f"({nb_bias[f'{CH} mean (uV)'] - nb_no[f'{CH} mean (uV)']:+.2f} uV) between no rejection and the biased "
      f"threshold, while the peak moves {nb_no[f'{CH} peak (uV)']:+.2f} -> {nb_bias[f'{CH} peak (uV)']:+.2f} uV "
      f"({nb_bias[f'{CH} peak (uV)'] - nb_no[f'{CH} peak (uV)']:+.2f} uV), and the averaged baseline noise moves "
      f"{nb_no['baseline noise (uV)']:.2f} -> {nb_bias['baseline noise (uV)']:.2f} uV.")
print("  The trial-count imbalance is what the free-response exercise (ex-2-5-consequence) is about: the two "
      "averages are no longer estimated with the same precision, and a peak measurement reads a noisier "
      "waveform as a larger one.")
print()
print("autoreject against fixed thresholds, same subject:")
print(l2.fmt_table(rows, list(rows[0]), floatfmt="{:.2f}"))
print()
print(f"Correction before rejection ({DETAIL}, ICA on a 1 Hz copy, {len(ica.exclude)} components removed): "
      f"peak-to-peak median {np.median(ptp_d):.0f} -> {np.median(ptp_clean):.0f} uV, and at a 100 uV threshold "
      f"the loss falls from "
      f"{next(r for r in sweeps[DETAIL]['rows'] if r['threshold_uv'] == 100.0)['percent_rejected_all']:.1f} % to "
      f"{next(r for r in sw_clean['rows'] if r['threshold_uv'] == 100.0)['percent_rejected_all']:.1f} % of trials.")