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

_needed = ("mne", "scipy", "matplotlib", "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)

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

import matplotlib.pyplot as plt
import numpy as np
import mne
import pooch
from scipy import stats
from scipy.signal import lfilter, sosfiltfilt

mne.set_log_level("WARNING")
# Silence the downloader at the source rather than scrubbing cache paths out of stored outputs
# afterwards: a re-execution would put them straight back (Phase 3 note in notebooks/README.md).
pooch.get_logger().setLevel("WARNING")
plt.rcParams["figure.dpi"] = 72

FS = float(L7.LATENCY_SPEC["sfreq_hz"])
BAND = tuple(L7.LATENCY_SPEC["band_hz"])
N_TAPS = int(L7.LATENCY_SPEC["fir_taps"])
BLOCK = int(L7.LATENCY_SPEC["block_samples"])
PROCESSING_MS = float(L7.LATENCY_SPEC["processing_ms"])
SUBJECTS = [f"S{i:03d}" for i in range(1, 13)]     # ds-eegbci baseline runs, 12 subjects
SEGMENT_S = 60.0
SEED = L7.SEED

print(f"MNE {mne.__version__}; helpers_l7 imported from notebooks/_shared")
print(f"pipeline: {FS:.0f} Hz, {BAND[0]:g}-{BAND[1]:g} Hz, {N_TAPS}-tap FIR, {BLOCK}-sample block, "
      f"{PROCESSING_MS:g} ms processing, causal")

taps = L7.fir_taps(N_TAPS, BAND, FS)
algebra = L7.fir_group_delay_samples(N_TAPS)

from scipy.signal import group_delay
_, gd = group_delay((taps, [1.0]), w=[2 * np.pi * f / FS for f in (8.0, 10.0, 12.0)])

imp = np.zeros(2048)
imp[512] = 1.0
resp = lfilter(taps, [1.0], imp)
peak = int(np.argmax(np.abs(resp))) - 512
centroid = L7.energy_centroid(np.arange(len(resp)), resp) - 512

print(f"{'method':52s} {'samples':>10s} {'ms':>9s}")
print(f"{'(N - 1) / 2, the algebra':52s} {algebra:10.4f} {algebra / FS * 1000:9.4f}")
for f, g in zip((8.0, 10.0, 12.0), gd):
    print(f"{f'scipy.signal.group_delay at {f:g} Hz':52s} {g:10.4f} {g / FS * 1000:9.4f}")
print(f"{'position of a filtered impulse':52s} {peak:10.4f} {peak / FS * 1000:9.4f}")
print(f"{'energy centroid of the impulse response':52s} {centroid:10.4f} "
      f"{centroid / FS * 1000:9.4f}")
print(f"\ntransition band for {N_TAPS} taps (Hamming, 3.3 fs / N): {3.3 * FS / N_TAPS:.2f} Hz — "
      f"this is what the {N_TAPS} taps are buying, and the 400 ms is what they cost")

for k, v in L7.BUFFER_CONVENTIONS.items():
    print(f"[{k}]\n  {v}\n")
print(f"at {FS:.0f} Hz with a {BLOCK}-sample block:")
print(f"  block-period   B / fs       = {BLOCK} / {FS:.0f} = {BLOCK / FS * 1000:.2f} ms")
print(f"  first-to-last  (B - 1) / fs = {BLOCK - 1} / {FS:.0f} = {(BLOCK - 1) / FS * 1000:.2f} ms")
print(f"  they differ by exactly one sample = {1 / FS * 1000:.2f} ms, for any block size")
print(f"  loop update rate = fs / B = {FS / BLOCK:.2f} Hz")

budget = L7.latency_budget()
L7.print_budget(budget, title=f"{FS:.0f} Hz, {BAND[0]:g}-{BAND[1]:g} Hz, {N_TAPS}-tap FIR, "
                              f"{BLOCK}-sample block, {PROCESSING_MS:g} ms processing, causal")

widget = {                 # site/notes/widget-latency.md, sections 2.1 and 2.2
    "FIR group delay, 129 taps (ms)": 400.00,
    "buffer, 32 samples, block period (ms)": 200.00,
    "buffer, 32 samples, first-to-last (ms)": 193.75,
    "total, worst case, block period (ms)": 610.00,
    "total, worst case, first-to-last (ms)": 603.75,
    "total, mean over block position (ms)": 510.00,
    "loop update rate (Hz)": 5.00,
    "zero-phase FIR look-ahead, 129 taps (ms)": 800.00,
    "FIR group delay, 65 taps (ms)": 200.00,
    "FIR group delay, 257 taps (ms)": 800.00,
    "FIR group delay, 513 taps (ms)": 1600.00,
    "transition band, 129 taps (Hz)": 4.09,
    "Butterworth group delay at 10 Hz, order 2 (samples)": 17.802755,
    "Butterworth group delay at 10 Hz, order 4 (samples)": 32.727137,
    "Butterworth group delay at 10 Hz, order 6 (samples)": 48.368038,
    "Butterworth group delay at 10 Hz, order 8 (samples)": 64.159335,
    "zero-phase Butterworth look-ahead, order 2 (samples)": 83,
    "zero-phase Butterworth look-ahead, order 4 (samples)": 145,
    "zero-phase Butterworth look-ahead, order 6 (samples)": 186,
    "zero-phase Butterworth look-ahead, order 8 (samples)": 226,
}
other = L7.latency_budget(convention="first-to-last")
mine = {
    "FIR group delay, 129 taps (ms)": L7.fir_group_delay_samples(129) / FS * 1000,
    "buffer, 32 samples, block period (ms)": BLOCK / FS * 1000,
    "buffer, 32 samples, first-to-last (ms)": (BLOCK - 1) / FS * 1000,
    "total, worst case, block period (ms)": budget["total_ms"],
    "total, worst case, first-to-last (ms)": other["total_ms"],
    "total, mean over block position (ms)": budget["mean_ms"],
    "loop update rate (Hz)": budget["update_rate_hz"],
    "zero-phase FIR look-ahead, 129 taps (ms)":
        L7.zero_phase_lookahead_samples(taps, fir=True) / FS * 1000,
    "FIR group delay, 65 taps (ms)": L7.fir_group_delay_samples(65) / FS * 1000,
    "FIR group delay, 257 taps (ms)": L7.fir_group_delay_samples(257) / FS * 1000,
    "FIR group delay, 513 taps (ms)": L7.fir_group_delay_samples(513) / FS * 1000,
    "transition band, 129 taps (Hz)": 3.3 * FS / 129,
}
for o in (2, 4, 6, 8):
    sos = L7.butter_sos(o, BAND, FS)
    mine[f"Butterworth group delay at 10 Hz, order {o} (samples)"] = \
        L7.iir_group_delay_samples(sos, 10.0, FS)
    mine[f"zero-phase Butterworth look-ahead, order {o} (samples)"] = \
        L7.zero_phase_lookahead_samples(sos, fir=False, fs=FS)

print(f"{'quantity':54s} {'this notebook':>14s} {'the widget':>12s} {'difference':>12s}")
print("-" * 96)
worst = 0.0
for k, v in widget.items():
    m = mine[k]
    d = m - v
    worst = max(worst, abs(d))
    print(f"{k:54s} {m:14.6f} {v:12.6f} {d:12.3e}")
print("-" * 96)
print(f"largest disagreement anywhere: {worst:.3e}")
assert worst < 5e-3, "a number disagrees with the widget: report both, change nothing"
print("Every quantity agrees.  The independent derivation and the widget's TypeScript, checked "
      "against SciPy, give the same answers.")
print("The one non-zero row is the transition band: site/notes/widget-latency.md quotes 4.09 Hz to "
      f"two decimals and the exact value is {3.3 * FS / 129:.4f} Hz.  That is the published figure "
      "being rounded, not a disagreement.")

from scipy.signal import sos2tf

print(f"{'order':>6s} {'sum of sections':>16s} {'phase derivative':>18s} {'expanded (b, a)':>17s} "
      f"{'ms':>9s}")
complaints = set()
for o in (2, 4, 6, 8):
    sos = L7.butter_sos(o, BAND, FS)
    per_section = L7.iir_group_delay_samples(sos, 10.0, FS)
    by_phase = L7.iir_group_delay_by_phase(sos, 10.0, FS)
    b, a = sos2tf(sos)
    with warnings.catch_warnings(record=True) as caught:
        warnings.simplefilter("always")
        _, expanded = group_delay((b, a), w=[2 * np.pi * 10.0 / FS])
    complaints.update(str(w.message).split(";")[0] for w in caught)
    print(f"{o:6d} {per_section:16.6f} {by_phase:18.6f} {float(expanded[0]):17.6f} "
          f"{per_section / FS * 1000:9.4f}")
print("\nThe first two columns agree to six decimals at every order.  The third is the expansion, and "
      "it is the one that says an order-8 filter is faster than an order-6 one.")
for c in sorted(complaints):
    print(f"  SciPy says so itself while computing that column: \"{c}\"")

print(f"{'filter':34s} {'causal group delay':>20s} {'zero-phase look-ahead':>23s} {'ratio':>7s}")
rows = [(f"FIR, {N_TAPS} taps",
         L7.fir_group_delay_samples(N_TAPS),
         L7.zero_phase_lookahead_samples(taps, fir=True))]
for o in (2, 4, 6, 8):
    sos = L7.butter_sos(o, BAND, FS)
    rows.append((f"Butterworth, order {o}",
                 L7.iir_group_delay_samples(sos, 10.0, FS),
                 L7.zero_phase_lookahead_samples(sos, fir=False, fs=FS)))
for name, causal, ahead in rows:
    print(f"{name:34s} {causal / FS * 1000:17.2f} ms {ahead / FS * 1000:20.2f} ms "
          f"{ahead / causal:7.2f}")
print("\nFor the FIR the look-ahead is exact: N - 1 samples, the support of the kernel, which is "
      "twice the causal group delay of the same filter.  For an IIR there is no exact answer, only "
      "a truncation criterion — how far before an impulse the forward-backward response still "
      "exceeds 1 % of its peak.  A tighter tolerance demands a longer wait, so the criterion is "
      "part of the number and is quoted with it.")
zp = L7.latency_budget(phase="zero-phase-live")
print(f"\nA 'zero-phase' loop that waits for its kernel rather than cheating: "
      f"{zp['total_ms']:.2f} ms end to end against {budget['total_ms']:.2f} ms causal — "
      f"{zp['total_ms'] / budget['total_ms']:.2f}x, for a filter that is usually described as "
      f"having no delay at all.")

budgets = {
    "FIR 129, causal": budget,
    "FIR 129, zero-phase live": zp,
    "Butterworth 4, causal": L7.latency_budget(sos=L7.butter_sos(4, BAND, FS)),
    "Butterworth 8, causal": L7.latency_budget(sos=L7.butter_sos(8, BAND, FS)),
}
fig, _ = L7.plot_latency_budget(budgets)
plt.show()

rng = np.random.default_rng(SEED)
x = rng.standard_normal(16000)

whole = lfilter(taps, [1.0], x)
online = L7.OnlineFIR(taps)
blockwise = np.concatenate([online.process(b)[0] for _, b in L7.stream_blocks(x, BLOCK)])
err = np.abs(blockwise - whole[: len(blockwise)]).max() / np.std(x)
print(f"causal FIR, carried state:  max |error| / signal RMS = {err:.3e}")

sos = L7.butter_sos(4, BAND, FS)
ref = sosfiltfilt(sos, x)
per_block = np.concatenate([sosfiltfilt(sos, b[0])      # scipy's own default padding
                            for _, b in L7.stream_blocks(x, BLOCK)])
err2 = np.std(per_block - ref[: len(per_block)]) / np.std(ref)
print(f"zero-phase, each block on its own: error RMS / signal RMS = {err2:.3f}")
print(f"\nThe first is machine precision.  The second is {100 * err2:.0f} % of the signal — the same "
      f"filter, applied the same number of times, to blocks instead of to a recording.  That is "
      f"pf-filter-across-boundaries, and it is why the loop carries state.")

ring = L7.RingBuffer(1, capacity=4 * BLOCK)
for _, b in L7.stream_blocks(x[:400], BLOCK):
    ring.push(b)
print(f"\nring buffer: capacity {ring.capacity} samples, {ring.n_pushed} pushed, "
      f"latest({BLOCK}) returns {ring.latest(BLOCK).shape} and matches the last block: "
      f"{np.allclose(ring.latest(BLOCK)[0], x[384 - BLOCK:384])}")

probe = L7.measure_loop_delay()
print(f"{'quantity':44s} {'measured':>12s} {'predicted':>12s} {'difference':>12s}")
print("-" * 84)
print(f"{'filter stream, energy-centroid lag (ms)':44s} {probe['filter_ms']:12.4f} "
      f"{probe['expected_filter_ms']:12.4f} "
      f"{probe['filter_ms'] - probe['expected_filter_ms']:12.2e}")
print(f"{'feedback value, energy-centroid lag (ms)':44s} {probe['value_ms']:12.4f} "
      f"{probe['expected_value_ms']:12.4f} "
      f"{probe['value_ms'] - probe['expected_value_ms']:12.2e}")
print("-" * 84)
print(f"spread over the {probe['alignments']} phases of the block grid: "
      f"{probe['value_min_ms']:.2f}-{probe['value_max_ms']:.2f} ms "
      f"(SD {probe['value_sd_ms']:.2f}); that spread IS the buffer term")
print(f"criterion: {probe['criterion']}")
print(f"\nThe budget's worst case is {budget['total_ms']:.2f} ms and its mean is "
      f"{budget['mean_ms']:.2f} ms; the probe measures the centre of mass, not the worst case, and "
      f"the simulation has no processing stage, so the number to compare it with is "
      f"{budget['mean_ms'] - PROCESSING_MS:.2f} ms plus the half sample above.")

fig, _ = L7.plot_feedback(probe, budget_ms=budget["total_ms"])
plt.show()

print(f"{'smoothing time constant':>24s} {'filter lag (ms)':>16s} {'feedback lag (ms)':>18s} "
      f"{'added (ms)':>11s}")
base = None
for smooth in (0.0, 0.25, 0.5, 1.0):
    pr = L7.measure_loop_delay(smooth_s=smooth, alignments=8)
    base = pr["value_ms"] if base is None else base
    print(f"{smooth:21.2f} s {pr['filter_ms']:16.2f} {pr['value_ms']:18.2f} "
          f"{pr['value_ms'] - base:11.2f}")
print("\nThe filter row does not move; the smoother is a fourth row the budget above does not have. "
      "A neurofeedback protocol that reports 'a 0.5-second smoothing window' has reported a quarter "
      "of a second of latency without saying so.")

L7.disk_report("before any download")
runs = ("R01", "R02")
already_on_disk = [f for s in SUBJECTS for f in L6.eegbci_run_files(s, runs)]
print(f"{len(already_on_disk)} of the {2 * len(SUBJECTS)} runs were already cached "
      f"(the ds-eegbci cache is shared with Levels 0, 4, 5 and 6, so only what this run fetches "
      f"is deleted)")

t0 = time.time()
eyes_open, eyes_closed = {}, {}
n = int(SEGMENT_S * FS)
for s in SUBJECTS:
    for run, store in (("R01", eyes_open), ("R02", eyes_closed)):
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            raw = helpers.load_spine("ds-eegbci", s, [run], verbose=False).pick(["O1"])
        assert raw.info["sfreq"] == FS, f"{s} {run}: {raw.info['sfreq']} Hz, not {FS}"
        store[s] = raw.get_data()[0][:n] * 1e6
print(f"{len(SUBJECTS)} subjects x 2 runs in {time.time() - t0:.1f} s; "
      f"{SEGMENT_S:.0f} s of O1 each at {FS:.0f} Hz")
assert len(eyes_open) == len(eyes_closed) == len(SUBJECTS)
L7.disk_report("after the downloads")

ec = np.array([L7.alpha_feedback(eyes_closed[s])["values"].mean() for s in SUBJECTS])
eo = np.array([L7.alpha_feedback(eyes_open[s])["values"].mean() for s in SUBJECTS])
t_ec, p_ec = stats.ttest_rel(ec, eo)
print(f"{'subject':>8s} {'eyes closed':>12s} {'eyes open':>11s} {'ratio':>7s}   (µV RMS per block)")
for s, a, b in zip(SUBJECTS, ec, eo):
    print(f"{s:>8s} {a:12.3f} {b:11.3f} {a / b:7.2f}")
print(f"\nmean {ec.mean():.3f} against {eo.mean():.3f} µV, {ec.mean() / eo.mean():.2f}x; "
      f"{int((ec > eo).sum())} of {len(SUBJECTS)} subjects higher with the eyes closed; "
      f"t({len(ec) - 1}) = {t_ec:.2f}, p = {p_ec:.2e}")
print("The loop measures something real.  Everything in section 6 is about whether a CHANGE in it "
      "during a session means anything.")

one = SUBJECTS[0]
fb = L7.alpha_feedback(eyes_closed[one])
truth = L7.offline_envelope(eyes_closed[one], FS, BAND)
t_s = np.arange(len(eyes_closed[one])) / FS
fig, axes = plt.subplots(2, 1, figsize=(11, 4.6), sharex=True)
axes[0].plot(t_s, eyes_closed[one], lw=0.5)
axes[0].set_ylabel("O1 (µV)")
axes[0].set_title(f"{one} R02, eyes closed, O1 — the replayed signal (µV)", fontsize=10)
axes[1].plot(t_s, truth, lw=0.8, label="offline zero-phase envelope (needs the future)")
axes[1].step(fb["times_s"], fb["values"], where="post", lw=1.1,
             label=f"live feedback, {fb['update_rate_hz']:.0f} Hz update (causal)")
axes[1].set_ylabel("Alpha amplitude\n(µV)")
axes[1].set_xlabel("Time (s)")
axes[1].legend(fontsize=8)
for a in axes:
    a.grid(alpha=0.3)
fig.suptitle(f"What the loop can see against what an offline analysis would report "
             f"(µV; {BAND[0]:g}-{BAND[1]:g} Hz)", fontsize=10)
fig.tight_layout()
plt.show()

for k, v in L7.SHAM_MODES.items():
    print(f"[{k}]\n  {v}\n")

one = SUBJECTS[0]
donor = eyes_open[SUBJECTS[1]]
runs_shown = {m: L7.sham_feedback(eyes_closed[one], donor=donor, mode=m)
              for m in ("veridical", "sham-band", "sham-yoked")}
for m, r in runs_shown.items():
    print(f"{m:14s} source: {r['source']}")
fig, _ = L7.plot_sham_comparison(runs_shown,
                                 title=f"{one}: what the participant would see, by condition "
                                       f"(µV RMS per {BLOCK / FS:.2f}-s block)")
plt.show()

modes = ("veridical", "sham-band", "sham-yoked", "sham-inverted")
trend = {}
print(f"{'condition':16s} {'p < .05':>9s} {'and rising':>11s} {'median |slope|':>15s} "
      f"{'median lag-1 r':>15s}")
for m in modes:
    ps, slopes, lags = [], [], []
    for i, s in enumerate(SUBJECTS):
        r = L7.sham_feedback(eyes_closed[s], donor=eyes_open[SUBJECTS[(i + 1) % len(SUBJECTS)]],
                             mode=m)
        tr = L7.within_session_trend(r["values"], r["times_s"])
        ps.append(tr["p"])
        slopes.append(tr["slope_per_min"])
        lags.append(tr["lag1_autocorrelation"])
    ps, slopes, lags = np.asarray(ps), np.asarray(slopes), np.asarray(lags)
    trend[m] = (ps, slopes, lags)
    print(f"{m:16s} {int((ps < 0.05).sum()):5d}/{len(ps):<3d} "
          f"{int(((ps < 0.05) & (slopes > 0)).sum()):7d}/{len(ps):<3d} "
          f"{np.median(np.abs(slopes)):15.4f} {np.median(lags):15.3f}")
print(f"\nnominal false-positive rate at alpha = .05: {0.05 * len(SUBJECTS):.1f} of {len(SUBJECTS)}")
print("The cause is in the last column.  The per-block feedback values are strongly "
      "autocorrelated, so the ordinary-least-squares standard error is far too small and the "
      "p-value is anti-conservative.  The test is invalid, not unlucky.")
print("\nAnd the sham conditions flag about as often as the veridical one, which is the whole point: "
      "'alpha went up during the session' carries no information about whether the feedback was "
      "real, because it does the same thing when the feedback comes from somebody else.")

print(f"{'condition':16s} {'second half - first half':>26s} {'t':>7s} {'p':>9s}")
session_level = {}
for m in modes:
    d = []
    for i, s in enumerate(SUBJECTS):
        r = L7.sham_feedback(eyes_closed[s], donor=eyes_open[SUBJECTS[(i + 1) % len(SUBJECTS)]],
                             mode=m)
        v = r["values"]
        h = len(v) // 2
        d.append(float(v[h:].mean() - v[:h].mean()))
    d = np.asarray(d)
    t_s, p_s = stats.ttest_1samp(d, 0.0)
    session_level[m] = (d, t_s, p_s)
    print(f"{m:16s} {d.mean():+21.4f} µV {t_s:7.2f} {p_s:9.4f}")

dv = session_level["veridical"][0]
ds = session_level["sham-yoked"][0]
t_c, p_c = stats.ttest_rel(dv, ds)
print(f"\nthe controlled contrast — veridical change minus yoked-sham change: "
      f"{dv.mean() - ds.mean():+.4f} µV, t({len(dv) - 1}) = {t_c:.2f}, p = {p_c:.4f}")
print("Correctly null, which is what a valid test does on data where the null is true.")

print("What a neurofeedback claim has to carry (CRED-nf-style reporting items):")
for i, item in enumerate(L7.CREDNF_ITEMS, 1):
    print(f"  {i}. {item}")
print(f"\nItem 4 in full, for the pipeline this notebook budgets:")
print(f"  band              {BAND[0]:g}-{BAND[1]:g} Hz")
print(f"  derivation        {L7.LATENCY_SPEC['channel']} "
      f"({L7.DATASETS_L7['ds-eegbci']['name']}, reference TODO(confirm) in the directory)")
print(f"  online filter     {N_TAPS}-tap linear-phase FIR, causal, applied block by block with "
      f"carried state")
print(f"  update rate       {budget['update_rate_hz']:.2f} Hz ({BLOCK} samples at {FS:.0f} Hz)")
print(f"  end-to-end latency {budget['total_ms']:.2f} ms worst case, {budget['mean_ms']:.2f} ms mean, "
      f"buffer quoted on the {budget['convention']} convention")
print("\nTODO(confirm): typical real-world processing delays.  The 10 ms row here is a stated number, "
      "not a measured one; this notebook cannot measure compute, transport or display and does not "
      "pretend to.  An author with a real rig should replace it with a measurement.")

print("=" * 96)
print("nb-7-3-realtime — the numbers the L7.3 exercises ask for")
print("=" * 96)
print(f"\n[1] END-TO-END LATENCY for the stated pipeline "
      f"({FS:.0f} Hz, {BAND[0]:g}-{BAND[1]:g} Hz, {N_TAPS}-tap FIR, {BLOCK}-sample block, "
      f"{PROCESSING_MS:g} ms processing, causal):")
for r in budget["rows"]:
    print(f"      {r['row']:12s} {r['detail']:44s} {r['ms']:8.2f} ms")
print(f"      {'TOTAL':12s} {'worst case, block-period convention':44s} "
      f"{budget['total_ms']:8.2f} ms   <- the answer")
print(f"      {'':12s} {'worst case, first-to-last convention':44s} "
      f"{budget['total_other_convention_ms']:8.2f} ms")
print(f"      {'':12s} {'mean over block position':44s} {budget['mean_ms']:8.2f} ms")
print(f"    The convention has to be named with the number: the three readings sit within 100 ms of "
      f"each other and a tolerance wide enough to admit {budget['mean_ms']:.0f} would accept an "
      f"answer that budgeted nothing.")

print(f"\n[2] AGREEMENT WITH w-latency-budget: every one of the {len(widget)} quantities in "
      f"site/notes/widget-latency.md reproduces, largest disagreement {worst:.1e}.  Derived here "
      f"independently; no number was taken from the widget and checked against itself.")

print(f"\n[3] MEASURED, not assumed: the loopback probe reads the filter row at "
      f"{probe['filter_ms']:.4f} ms against {probe['expected_filter_ms']:.4f} predicted, and the "
      f"feedback value at {probe['value_ms']:.4f} ms against {probe['expected_value_ms']:.4f}; "
      f"the {probe['value_min_ms']:.2f}-{probe['value_max_ms']:.2f} ms spread over the block grid "
      f"is the buffer term.")

print(f"\n[4] A zero-phase filter cannot run live: the same {N_TAPS}-tap kernel needs "
      f"{L7.zero_phase_lookahead_samples(taps, fir=True) / FS * 1000:.0f} ms of FUTURE signal, so "
      f"the honest budget is {zp['total_ms']:.2f} ms — {zp['total_ms'] / budget['total_ms']:.2f}x "
      f"the causal one.  Block-wise causal filtering is machine-precision identical to whole-signal "
      f"filtering ({err:.1e}); filtering each block on its own is not ({100 * err2:.0f} % error, with "
      f"scipy's own default padding).")

print(f"\n[5] POSITIVE CONTROL: eyes closed {ec.mean():.3f} µV against eyes open {eo.mean():.3f} µV "
      f"through the same loop, {ec.mean() / eo.mean():.2f}x, "
      f"{int((ec > eo).sum())}/{len(SUBJECTS)} subjects, t({len(ec) - 1}) = {t_ec:.2f}, "
      f"p = {p_ec:.2e}.")

print(f"\n[6] THE SHAM RESULT.  On {len(SUBJECTS)} replayed recordings where learning is impossible "
      f"by construction, the naive within-session trend test called:")
for m in modes:
    ps, slopes, _ = trend[m]
    print(f"      {m:16s} {int((ps < 0.05).sum())}/{len(ps)} sessions significant, "
          f"{int(((ps < 0.05) & (slopes > 0)).sum())}/{len(ps)} of them 'an increase'")
print(f"      nominal rate at alpha = .05: {0.05 * len(SUBJECTS):.1f}/{len(SUBJECTS)}.  "
      f"Median lag-1 autocorrelation of the feedback values: "
      f"{np.median(trend['veridical'][2]):.3f}.")
print(f"      With the session as the unit instead of the block, every condition is null "
      f"(p = " + ", ".join(f"{session_level[m][2]:.3f}" for m in modes) + "), and the controlled "
      f"contrast is null too (p = {p_c:.4f}).")
print(f"      The design flaw the exercise asks for has a shape: a protocol with no control "
      f"condition, whose outcome is a within-session trend in the feedback signal itself, reports "
      f"the same thing whether the feedback was the participant's own brain or somebody else's.")
print("=" * 96)

# Only what this run downloaded is deleted: the ds-eegbci cache is shared with four other levels.
freed = 0
for s in SUBJECTS:
    freed += L6.drop_eegbci_runs(s, runs, keep=already_on_disk)
print(f"deleted {freed / 1e6:.1f} MB of ds-eegbci runs fetched by this notebook")
L7.disk_report("at the end")
left = {f.name for s in SUBJECTS for f in L6.eegbci_run_files(s, runs)}
assert left <= {f.name for f in already_on_disk}, "this run left a download behind"
print(f"ds-eegbci runs still on disk: {len(left)} (all pre-existing); "
      f"no download made by this run was left behind")