# 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')
_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_l4.py").exists()), None)
if _shared is None:
    raise FileNotFoundError("start the kernel in notebooks/L4/ (or notebooks/) so that "
                            "_shared/helpers_l4.py is found")
sys.path.insert(0, str(_shared))
import helpers
import helpers_l1
import helpers_l4 as L4

# 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

# 4. Quiet the downloader.  pooch, which MNE uses to fetch datasets, logs
#    "Downloading file '...' from '...' to '<cache directory>'" at INFO, and that last field is an
#    ABSOLUTE PATH from whichever machine executed the notebook.  Absolute paths are not allowed in a
#    stored notebook (scripts/scrub-notebooks.py is a CI gate) and re-executing would put them straight
#    back, so the message is suppressed at the source rather than cleaned up afterwards.  Nothing is
#    hidden by this: every cell below prints the file NAMES it fetched and helpers_l4.Downloads prints
#    free disk before and after.  Please do not delete this as noise.
try:
    import pooch

    pooch.get_logger().setLevel("WARNING")
except Exception:  # pooch absent or its API moved: the scrub script is the backstop
    pass
print(f"MNE {mne.__version__}; helpers_l4 imported from notebooks/_shared")
print(f"downloads go to {helpers_l1.download_dir().name}/ (resolved relative to the working directory, "
      "or $EEG_COURSE_DOWNLOADS) and are deleted at the end of this notebook")

SUBJECT, CH, COND = "S001", "C3", "T2"
ASSET_RUN, ASSET_T0, ASSET_DUR = 8, 61.3, 3.0     # w-wavelet-explorer/epoch.json

dl = L4.Downloads("nb-4-3").start()
raw_r08 = L4.load_imagery_raw(SUBJECT, [ASSET_RUN], dl)
sf = float(raw_r08.info["sfreq"])
i0 = int(round(ASSET_T0 * sf))
seg = raw_r08.get_data(picks=[CH])[0, i0:i0 + int(round(ASSET_DUR * sf))] * 1e6
t_seg = np.arange(seg.size) / sf + (ASSET_T0 - 62.3)      # 0 s = the cue
print(f"reproduced the shipped epoch: {SUBJECT} R{ASSET_RUN:02d}, t0 = {ASSET_T0:g} s, {ASSET_DUR:g} s, "
      f"{CH}, average reference -> {seg.size} samples at {sf:g} Hz")
print(f"  largest absolute sample {np.abs(seg).max():.1f} uV "
      f"(w-wavelet-explorer/epoch.json records max_abs_uv 53.8 for its chosen trial)")
print("downloaded so far: " + ", ".join(p.name for p in dl.paths))

from scipy import signal as sps

FREQS = L4.TF_FREQS
NCYC = L4.TF_CYCLES

mor = L4.morlet_power(seg[None, None, :], sfreq=sf, freqs=FREQS, n_cycles=NCYC, decim=1, times=t_seg)
mt = L4.multitaper_power(seg[None, None, :], sfreq=sf, freqs=FREQS, n_cycles=NCYC, time_bandwidth=4.0,
                         decim=1, times=t_seg)
print(mor["call"])
print(mt["call"])
print(f"multitaper: time_bandwidth {mt['time_bandwidth']:g} -> {mt['n_tapers']} tapers; half bandwidth "
      f"{mt['half_bandwidth_hz'][np.argmin(abs(FREQS - 10))]:.3f} Hz at 10 Hz and "
      f"{mt['half_bandwidth_hz'][-1]:.3f} Hz at {FREQS[-1]:g} Hz")

# Filter-Hilbert, one narrow band per frequency row, matched to the wavelet's own FWHM.
# A zero-phase filtfilt needs 3 * numtaps samples of signal, and numtaps grows as the pass band
# narrows -- so on a 3-s epoch at 160 Hz the low rows are simply not computable. That is not a
# bug to work around; it is the constraint filter-Hilbert carries and the reason the row is left
# empty and reported rather than quietly filled.
fh = np.full((len(FREQS), seg.size), np.nan)
bw = np.zeros(len(FREQS))
taps_n = np.zeros(len(FREQS), dtype=int)
feasible = np.zeros(len(FREQS), dtype=bool)
for i, (f0, nc) in enumerate(zip(FREQS, NCYC)):
    w = L4.wavelet_fwhm(f0, nc)["fwhm_f_hz"]
    lo, hi = max(1.0, f0 - w / 2), min(sf / 2 - 1.0, f0 + w / 2)
    bw[i] = hi - lo
    trans = max(1.0, 0.25 * lo)
    n_t = int(round(3.3 * sf / trans))
    n_t += (n_t + 1) % 2
    taps_n[i] = n_t
    if 3 * n_t < seg.size:
        feasible[i] = True
        fh[i] = L4.filter_hilbert(seg, sf, (lo, hi))["envelope"] ** 2
f_min = FREQS[feasible][0] if feasible.any() else float("nan")
print(f"filter-Hilbert: one zero-phase Hamming FIR per row, pass band = the wavelet's own FWHM "
      f"({bw[np.argmin(abs(FREQS - 10))]:.3f} Hz wide at 10 Hz), envelope squared; "
      f"{L4.BURST_CONVENTIONS['filter length']}")
print(f"  filtfilt needs 3 x numtaps = {3 * taps_n[0]} samples for the {FREQS[0]:g} Hz row and this epoch "
      f"has {seg.size}. {int((~feasible).sum())} of {len(FREQS)} rows are therefore NOT COMPUTABLE on a "
      f"{ASSET_DUR:g}-s epoch at {sf:g} Hz; the lowest row that is, is {f_min:g} Hz "
      f"({taps_n[int(np.argmax(feasible))]} taps = {taps_n[int(np.argmax(feasible))] / sf:.2f} s).")
print(f"  Morlet and multitaper have no such floor because they pad; they just return an edge-dominated "
      f"estimate there instead, which is worse in a different way (nb-4-2 section 3).")

fig, axes = plt.subplots(1, 3, figsize=(17.5, 3.9))
for ax, (name, P) in zip(axes, (("Morlet", mor["power"][0, 0]), ("Multitaper", mt["power"][0, 0]),
                                (f"Filter-Hilbert ({int(feasible.sum())}/{len(FREQS)} rows computable)", fh))):
    med = np.array([np.median(row) if np.isfinite(row).all() else np.nan for row in P])[:, None]
    L4.plot_tfr(10 * np.log10(P / med), FREQS, t_seg, ax=ax,
                vlim=(-12, 12), event_s=0.0, edge_s=L4.edge_seconds(FREQS, NCYC),
                cbar_label="Power re the row's own median (dB)",
                title=f"{name}\n{SUBJECT} R{ASSET_RUN:02d} {CH}, one trial (dB re each row's median)")
fig.suptitle("One 3-s epoch through three estimators (each row normalised to its own median so the three "
             "share a colour bar; the units are not the same, section 2)", y=1.04)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

cal_mor = L4.unit_cosine_power(sf, FREQS, NCYC, estimator="morlet")
cal_mt = L4.unit_cosine_power(sf, FREQS, NCYC, estimator="multitaper", time_bandwidth=4.0)
rows = [{"f (Hz)": f, "Morlet constant": a, "Multitaper constant": b, "ratio mt/morlet": b / a}
        for f, a, b in zip(FREQS, cal_mor, cal_mt) if f in (4, 8, 10, 13, 20, 30, 40)]
print("What each call returns for a unit-amplitude cosine (library units per uV^2 of mean square x 2):")
print(L4.fmt_table(rows, floatfmt="{:.4f}"))

P_mor = mor["power"][0, 0] / (2 * cal_mor[:, None])          # mean-square uV^2
P_mt = mt["power"][0, 0] / (2 * cal_mt[:, None])
P_fh = fh / 2.0                                              # |envelope|^2 / 2 is the mean square of a sinusoid
inner = np.abs(t_seg - t_seg.mean()) < (ASSET_DUR / 2 - L4.edge_seconds(FREQS, NCYC).max())
print(f"\nInterior window used for every comparison below: {t_seg[inner][0]:+.3f} .. {t_seg[inner][-1]:+.3f} s "
      f"({inner.sum()} of {seg.size} samples), so no edge-contaminated estimate enters a correlation.")

# The comparison the L4.3 exercise keys is precisely specified, and every word of it matters:
#   Pearson r of LOG10 power, LIBRARY (uncalibrated) units, pooled over all frequency rows, restricted
#   to samples more than ONE edge margin (the global edge_s, 0.597 s here) from either end.
# All four combinations of {log10, linear} x {one global margin, a per-frequency margin} are computed
# below, because they are not the same number and a tight tolerance can only admit one of them.
def agreement(A, B, rows_ok=None, *, log10=True, per_freq=False):
    rows_ok = np.ones(len(FREQS), bool) if rows_ok is None else rows_ok
    a, b = [], []
    for i in np.where(rows_ok)[0]:
        m = (t_seg >= t_seg[0] + (edge_by_f[i] if per_freq else edge_by_f.max())) & \
            (t_seg <= t_seg[-1] - (edge_by_f[i] if per_freq else edge_by_f.max()))
        x, y = A[i][m], B[i][m]
        if log10:
            x, y = np.log10(x), np.log10(y)
        a.append(x)
        b.append(y)
    a, b = np.concatenate(a), np.concatenate(b)
    return float(np.corrcoef(a, b)[0, 1]), int(a.size)

edge_by_f = L4.edge_seconds(FREQS, NCYC)
MOR_LIB, MT_LIB = mor["power"][0, 0], mt["power"][0, 0]
variants = []
for log10_ in (True, False):
    for per_freq_ in (False, True):
        r_, n_ = agreement(MOR_LIB, MT_LIB, log10=log10_, per_freq=per_freq_)
        variants.append({"scale": "log10 power" if log10_ else "linear power",
                         "edge margin": "one global (edge_s)" if not per_freq_ else "per frequency row",
                         "units": "library", "r": r_, "cells": n_})
r_cal, _ = agreement(P_mor, P_mt, log10=True, per_freq=False)
variants.append({"scale": "log10 power", "edge margin": "one global (edge_s)",
                 "units": "calibrated mean-square uV^2", "r": r_cal, "cells": variants[0]["cells"]})
print()
print("Morlet against multitaper, the same two maps, four ways of asking the same question:")
print(L4.fmt_table(variants, floatfmt="{:.4f}"))

pairs = {"Morlet vs multitaper": (P_mor, P_mt, np.ones(len(FREQS), bool)),
         "Morlet vs filter-Hilbert": (P_mor, P_fh, feasible),
         "multitaper vs filter-Hilbert": (P_mt, P_fh, feasible)}
agree = []
for name, (A, B, rows_ok) in pairs.items():
    a = A[np.ix_(rows_ok, np.where(inner)[0])].ravel()
    b = B[np.ix_(rows_ok, np.where(inner)[0])].ravel()
    agree.append({"pair": name + ("" if rows_ok.all() else f" ({int(rows_ok.sum())} rows)"),
                  "r (calibrated uV^2)": float(np.corrcoef(a, b)[0, 1]),
                  "mean ratio": float(np.mean(b / a)), "median ratio": float(np.median(b / a))})
a_lib = MOR_LIB[:, inner].ravel()
b_lib = MT_LIB[:, inner].ravel()
agree.append({"pair": "Morlet vs multitaper, LIBRARY units",
              "r (calibrated uV^2)": float(np.corrcoef(a_lib, b_lib)[0, 1]),
              "mean ratio": float(np.mean(b_lib / a_lib)), "median ratio": float(np.median(b_lib / a_lib))})
print()
print("And on a linear scale, with the ratios beside the correlations:")
print(L4.fmt_table(agree, floatfmt="{:.4f}"))
print()
print("The last two rows are the same two maps. The correlation barely moves between library and calibrated "
      "units (the calibration is a per-row constant, and a correlation over all rows still feels it); the "
      "RATIO moves a great deal, because that is exactly what the calibration is.")
print(f"Unit-cosine ratio, multitaper over Morlet, median over rows: "
      f"{np.median(cal_mt / cal_mor):.4f} -- on a pure tone the two estimators return nearly the same "
      f"number, so the factor of ~2.4 on this epoch is not a normalisation difference.")

rng = np.random.default_rng(20260918)
SF_S = 500.0
N_S = int(20 * SF_S)
t_s = np.arange(N_S) / SF_S
truth = 6.0 * np.cos(2 * np.pi * 10.0 * t_s)
noise = L4.pink_noise(N_S, SF_S, 1.343, 20.0, rng)
sig = truth + noise
probe = np.arange(6.0, 15.1, 0.25)
probe_nc = np.full(len(probe), 7.0)
mo = L4.morlet_power(sig[None, None, :], sfreq=SF_S, freqs=probe, n_cycles=probe_nc, decim=5, times=t_s)
mu_ = L4.multitaper_power(sig[None, None, :], sfreq=SF_S, freqs=probe, n_cycles=probe_nc,
                          time_bandwidth=4.0, decim=5, times=t_s)
cal_a = L4.unit_cosine_power(SF_S, probe, probe_nc, estimator="morlet")
cal_b = L4.unit_cosine_power(SF_S, probe, probe_nc, estimator="multitaper", time_bandwidth=4.0)
A = mo["power"][0, 0] / (2 * cal_a[:, None])
B = mu_["power"][0, 0] / (2 * cal_b[:, None])
keep = (mo["times"] > 2.0) & (mo["times"] < 18.0)
i10 = int(np.argmin(abs(probe - 10.0)))

fig, axes = plt.subplots(1, 2, figsize=(13.5, 4.0))
axes[0].plot(mo["times"][keep], A[i10][keep], lw=1.0, label=f"Morlet, 7 cycles (CV {A[i10][keep].std()/A[i10][keep].mean():.3f})")
axes[0].plot(mu_["times"][keep], B[i10][keep], lw=1.0,
             label=f"Multitaper, 7 cycles, TW 4 ({mu_['n_tapers']} tapers; CV {B[i10][keep].std()/B[i10][keep].mean():.3f})")
axes[0].axhline(0.5 * 6.0 ** 2, color="k", lw=1.2, ls="--", label="true mean-square power of the 6 uV tone")
axes[0].set(xlabel="Time (s)", ylabel="Power (uV^2, mean square)",
            title="SYNTHETIC 10 Hz 6 uV tone on 1/f noise of 20 uV SD:\n10 Hz row of each estimate (uV^2)")
axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3)
axes[1].plot(probe, A[:, keep].mean(1), "-o", ms=3, label="Morlet")
axes[1].plot(probe, B[:, keep].mean(1), "-s", ms=3, label="Multitaper, TW 4")
axes[1].axvline(10.0, color="k", lw=1.0, ls="--")
axes[1].set(xlabel="Frequency (Hz)", ylabel="Mean power over the interior (uV^2, mean square)",
            title="The same estimates across frequency (uV^2):\nmultitaper is smoother and wider")
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(f"{'estimator':26s} {'peak f (Hz)':>12s} {'peak power (uV^2)':>18s} {'CV over time at 10 Hz':>22s} "
      f"{'-3 dB width (Hz)':>17s}")
for nm, M in (("Morlet, 7 cycles", A), (f"Multitaper, TW 4 ({mu_['n_tapers']} tapers)", B)):
    prof = M[:, keep].mean(1)
    half = prof.max() / 2
    over = probe[prof >= half]
    print(f"{nm:26s} {probe[int(prof.argmax())]:12.2f} {prof.max():18.3f} "
          f"{M[i10][keep].std() / M[i10][keep].mean():22.4f} {over[-1] - over[0]:17.3f}")
print(f"True mean-square power of the tone: {0.5 * 6.0 ** 2:.3f} uV^2. Both estimators sit above it because "
      f"the 1/f background contributes at 10 Hz as well; neither is measuring the tone alone.")

SF_G = 500.0
N_G = int(12 * SF_G)
t_g = np.arange(N_G) / SF_G
rng = np.random.default_rng(4321)
# a broadband 60-90 Hz burst: white noise band-passed to 60-90 Hz, windowed, on a 1/f background
wide = L4.filter_hilbert(rng.standard_normal(N_G), SF_G, (60.0, 90.0))["filtered"]
wide = wide / wide.std() * 4.0
env = np.exp(-((t_g - 6.0) ** 2) / (2 * 0.6 ** 2))
gamma = wide * env
bg = L4.pink_noise(N_G, SF_G, 1.343, 20.0, rng)
x_g = gamma + bg

gf = np.arange(45.0, 120.1, 2.5)
res = {}
for label, nc, tw in (("Morlet, 7 cycles", 7.0, None), ("Morlet, 3 cycles", 3.0, None),
                      ("Multitaper, 7 cycles, TW 4", 7.0, 4.0), ("Multitaper, 7 cycles, TW 8", 7.0, 8.0)):
    ncv = np.full(len(gf), nc)
    if tw is None:
        r = L4.morlet_power(x_g[None, None, :], sfreq=SF_G, freqs=gf, n_cycles=ncv, decim=5, times=t_g)
        cal = L4.unit_cosine_power(SF_G, gf, ncv, estimator="morlet")
    else:
        r = L4.multitaper_power(x_g[None, None, :], sfreq=SF_G, freqs=gf, n_cycles=ncv, time_bandwidth=tw,
                                decim=5, times=t_g)
        cal = L4.unit_cosine_power(SF_G, gf, ncv, estimator="multitaper", time_bandwidth=tw)
    res[label] = (r["power"][0, 0] / (2 * cal[:, None]), r["times"], r.get("n_tapers", 1))

fig, axes = plt.subplots(1, 4, figsize=(21, 3.9))
stats = []
for ax, (label, (P, tt, ntap)) in zip(axes, res.items()):
    band = (gf >= 60) & (gf <= 90)
    base = (tt < 3.0)
    L4.plot_tfr(10 * np.log10(P / P[:, base].mean(1, keepdims=True)), gf, tt, ax=ax, vlim=(-10, 10),
                event_s=6.0, cbar_label="Power change from 0-3 s (dB)",
                title=f"{label}\nSYNTHETIC 60-90 Hz burst at 6 s (dB re 0-3 s)")
    trace = P[band].mean(0)
    quiet = trace[base]
    stats.append({"estimator": label, "tapers": ntap,
                  "peak 60-90 Hz power (uV^2)": float(trace.max()),
                  "background CV before the burst": float(quiet.std() / quiet.mean()),
                  "peak / background": float(trace.max() / quiet.mean())})
fig.suptitle("A broadband 60-90 Hz burst (SYNTHETIC, 500 Hz), four estimators", y=1.04)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline
print(L4.fmt_table(stats, floatfmt="{:.4f}"))
print()
print("What the table says: more tapers -> lower background variability (the CV column) at the cost of wider "
      "frequency smoothing, which for a band 30 Hz wide costs nothing that matters. The detectability column "
      "(peak / background) is the one a gamma analysis cares about.")

rng = np.random.default_rng(99)
SF_N = 250.0
noise_only = rng.standard_normal(int(6 * SF_N)) * 10.0
t_n = np.arange(noise_only.size) / SF_N
fig, axes = plt.subplots(1, 3, figsize=(17, 3.6))
axes[0].plot(t_n, noise_only, lw=0.6, color="0.4")
axes[0].set(xlabel="Time (s)", ylabel="Amplitude (uV)",
            title="SYNTHETIC white noise, 10 uV SD -- no oscillation of any kind (uV)")
axes[0].grid(alpha=0.3)
for ax, band in zip(axes[1:], ((8.0, 13.0), (9.5, 10.5))):
    r = L4.filter_hilbert(noise_only, SF_N, band)
    ax.plot(t_n, r["filtered"], lw=0.9, color="tab:blue", label=f"{band[0]:g}-{band[1]:g} Hz, "
                                                                f"{r['numtaps']} taps")
    ax.plot(t_n, r["envelope"], lw=1.2, color="tab:red", label="Hilbert envelope")
    cyc = np.diff(np.unwrap(np.angle(sps.hilbert(r["filtered"])))) * SF_N / (2 * np.pi)
    ax.set(xlabel="Time (s)", ylabel="Amplitude (uV)",
           title=f"The same noise, band-passed {band[0]:g}-{band[1]:g} Hz (uV)\n"
                 f"median instantaneous frequency {np.median(cyc):.2f} Hz")
    ax.legend(fontsize=8)
    ax.grid(alpha=0.3)
fig.suptitle("A narrow filter manufactures a rhythm (pf-narrowband-filter-oscillation)", y=1.04)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline
print("Neither band-passed trace contains an oscillation: the input was white noise. Both have a smooth "
      "envelope and a well-behaved instantaneous frequency close to the centre of the pass band, because "
      "that is what a band-pass filter does to anything.")

try:
    WIDGET_R, WIDGET_RATIO = 0.7099, 2.398        # site/notes/integration-phase3.md, widgets-J
    lib_row = [a for a in agree if "LIBRARY" in a["pair"]][0]
    print("nb-4-3-multitaper-hilbert -- L4.3 numbers (draft; TODO(confirm) at author review)")
    print(f"Data: ds-eegbci {SUBJECT} run R{ASSET_RUN:02d}, the 3-s epoch at t0 = {ASSET_T0:g} s that "
          f"w-wavelet-explorer/epoch.json ships (channel {CH}, average reference, no filtering, "
          f"{sf:g} Hz); PhysioNet DOI {L4.DATASETS_L4['ds-eegbci']['dataset_doi']}, ODC-By 1.0.")
    print(f"Estimators: Morlet (max(3, f/2) cycles), multitaper (same cycles, time_bandwidth 4 -> "
          f"{mt['n_tapers']} tapers), filter-Hilbert (zero-phase Hamming FIR per row, pass band = the "
          f"wavelet's own FWHM). MNE {mne.__version__}.")
    print()
    print("L4.3's exercise is a multiple choice with reasoning, so there is no single numeric key. What this "
          "notebook supplies for it:")
    print(f"    ANSWER -- for a 60-90 Hz analysis, multitaper. The band is about 30 Hz wide and the activity "
          f"in it is broadband, so resolving inside it is not the goal; averaging tapers buys variance "
          f"reduction, and the frequency smoothing it costs is smaller than the band.")
    print(f"    Measured on this notebook's synthetic 60-90 Hz burst at 500 Hz:")
    for s in stats:
        print(f"        {s['estimator']:28s} {s['tapers']} taper(s): background CV "
              f"{s['background CV before the burst']:.4f}, peak/background {s['peak / background']:.2f}")
    print(f"    AND: ds-eegbci at {sf:g} Hz CANNOT answer this question -- Nyquist is {sf / 2:g} Hz, so 90 Hz "
          f"is absent, and 60 Hz is this recording's unfiltered mains line. A lesson that runs a 60-90 Hz "
          f"analysis on it is analysing padding and line noise.")
    print()
    print("Supporting -- how well the estimators agree on the real epoch (interior only, edges excluded):")
    for a in agree:
        print(f"    {a['pair']:38s} r = {a['r (calibrated uV^2)']:.4f}, "
              f"median ratio {a['median ratio']:.4f}, mean ratio {a['mean ratio']:.4f}")
    key_row = variants[0]      # log10 power, library units, one global edge margin
    print(f"    ANSWER KEY -- ex-4-3-agreement (numeric): Pearson r of LOG10 power between the Morlet and "
          f"multitaper maps of this epoch, LIBRARY units, pooled over all {len(FREQS)} rows, restricted to "
          f"samples more than ONE edge margin ({edge_by_f.max():.3f} s) from either end.")
    print(f"        r = {key_row['r']:.4f} over {key_row['cells']} cells.")
    print(f"    Cross-check: the lesson's key is r = {WIDGET_R:g} with tolerance +/- 0.004 "
          f"(widgets-J via site/notes/integration-phase3.md, recorded in "
          f"w-wavelet-explorer/multitaper.json as morlet_agreement.pearson_r_log10_power). This notebook "
          f"{'AGREES' if abs(key_row['r'] - WIDGET_R) < 0.004 else 'DISAGREES'} "
          f"(difference {key_row['r'] - WIDGET_R:+.4f}).")
    print(f"    Every word of the definition is load-bearing at that tolerance. The same two maps give:")
    for v in variants:
        inside = abs(v["r"] - WIDGET_R) < 0.004
        print(f"        r = {v['r']:.4f}  {v['scale']:13s} {v['units']:26s} edge margin "
              f"{v['edge margin']:20s} -> {'inside' if inside else 'OUTSIDE'} the +/- 0.004 band")
    print(f"    FINDING for content-L4 and the integrator: the key's own comment says the correlation is "
          f"'restricted to samples more than edge_s_by_freq from either end', i.e. a PER-FREQUENCY margin, "
          f"while the asset's morlet_agreement.measured_over says 'more than edge_s = 0.597 s from either "
          f"end', i.e. ONE GLOBAL margin, and the prompt says 'more than one edge margin'. The global "
          f"reading gives {variants[0]['r']:.4f} and reproduces the key; the per-frequency reading gives "
          f"{variants[1]['r']:.4f}, which the +/- 0.004 tolerance EXCLUDES. Two of the three sentences "
          f"describing the key agree with each other and the third does not; the key itself is right and "
          f"the comment's phrase should be corrected to edge_s. Nothing here was tuned -- the global "
          f"reading was tried because the asset named it.")
    print(f"    The calibrated counterpart the lesson's tolerance is set to exclude: r = {r_cal:.4f} "
          f"(the lesson's note records 0.7001, difference {r_cal - 0.7001:+.4f}).")
    print()
    print(f"    The RATIO, also keyed at {WIDGET_RATIO:g}, agrees only when read as a MEDIAN: "
          f"{lib_row['median ratio']:.4f} against {WIDGET_RATIO:g} "
          f"(difference {lib_row['median ratio'] - WIDGET_RATIO:+.4f}), while the MEAN ratio over the "
          f"same cells is {lib_row['mean ratio']:.4f}. The per-cell ratio of two power maps is heavy-tailed "
          f"-- a handful of cells where the Morlet estimate is near zero dominate any average of it -- so "
          f"'the ratio is 2.4' is a statement about the median and has to say so.")
    print(f"    The asset's other three numbers reproduce too: unit-cosine ratio "
          f"{np.median(cal_mt / cal_mor):.4f} (asset 1.1255), median ratio after calibration "
          f"{np.median((P_mt / P_mor)[:, inner]):.4f} (asset 2.132).")
    print()
    print("Supporting -- the calibration constants, without which no ratio means anything:")
    print(f"    Morlet at 10 Hz: {cal_mor[np.argmin(abs(FREQS - 10))]:.4f}; multitaper at 10 Hz: "
          f"{cal_mt[np.argmin(abs(FREQS - 10))]:.4f} (library units for a unit-amplitude cosine).")
    print()
    print("Supporting -- a constraint filter-Hilbert has and the other two do not:")
    print(f"    On a {ASSET_DUR:g}-s epoch at {sf:g} Hz, {int((~feasible).sum())} of {len(FREQS)} frequency "
          f"rows cannot be filtered at all: a zero-phase filtfilt needs 3 x numtaps samples and the "
          f"{FREQS[0]:g} Hz row's filter is {taps_n[0]} taps ({taps_n[0] / sf:.2f} s). The lowest "
          f"computable row is {f_min:g} Hz. Morlet and multitaper return a number there because they pad, "
          f"which is not the same as having data.")
    print()
    print("Pitfall: pf-narrowband-filter-oscillation -- section 5 band-passes white noise and gets a rhythm. "
          "Widget: w-wavelet-explorer (mode compare).")
finally:
    dl.finish()