# 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', 'specparam')
_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
from scipy import signal as sps

try:
    from bycycle.features import compute_features as _bycycle_features
    HAVE_BYCYCLE = True
except ImportError:  # optional dependency; section 4 falls back and says so
    _bycycle_features = None
    HAVE_BYCYCLE = False

# 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")
print(f"bycycle available: {HAVE_BYCYCLE}" if HAVE_BYCYCLE else "bycycle not installed -- section 4 uses the documented NumPy fallback")

REF_SUBJECT = "sub-010002"
REF_BLOCK = 6
DURATION_S = 60.0
TARGET_SFREQ = 250.0
OCCIPITAL = ["O1", "Oz", "O2"]

dl = L4.Downloads("nb-4-6").start()
tmpdir = Path(helpers_l1.download_dir()) / "lemon-windows"
hdr = L4.lemon_headers(REF_SUBJECT)
blocks = L4.lemon_blocks(hdr)
print(f"{REF_SUBJECT} (mirror id {hdr['mirror_id']}): {hdr['n_channels']} channels, {hdr['sfreq']:g} Hz, "
      f"{hdr['binary_format']} {hdr['orientation']}, {len(hdr['markers'])} stimulus markers -> "
      f"{len(blocks)} blocks")
for b in blocks[:8]:
    print(f"    block {b['index']:2d}: code {b['code']} {b['t0_s']:8.3f}..{b['t1_s']:8.3f} s "
          f"({b['duration_s']:.1f} s)")
print(f"    ... {len(blocks)} blocks in total")
print(f"the whole .eeg is {hdr['n_channels'] * 2 * int(hdr['sfreq']) / 1e6:.2f} MB per second of recording; "
      f"{DURATION_S:g} s is {DURATION_S * hdr['n_channels'] * 2 * hdr['sfreq'] / 1e6:.1f} MB")

def load_window(subject, headers, t0_s, duration_s, downloads, resample_hz=TARGET_SFREQ):
    """Range-fetch one window and return it as a Raw resampled to `resample_hz` (uV in, uV out)."""
    vhdr = L4.fetch_lemon_window(subject, t0_s, duration_s, tmpdir, downloads, headers=headers,
                                 verbose=True)
    if vhdr is None:
        return None
    r = mne.io.read_raw_brainvision(vhdr, preload=True, verbose=False)
    if resample_hz and abs(r.info["sfreq"] - resample_hz) > 1e-6:
        r.filter(None, resample_hz * 0.4, picks="all", method="fir", phase="zero", fir_design="firwin",
                 verbose=False)
        r.resample(resample_hz, npad="auto", verbose=False)
    return r


def relative_alpha(raw_obj, picks=OCCIPITAL, band=(8.0, 12.0), broad=(1.0, 40.0)):
    """Median over 4-s windows of log10(alpha power / broadband power), averaged over `picks`."""
    sfq = float(raw_obj.info["sfreq"])
    d = raw_obj.get_data(picks=[c for c in picks if c in raw_obj.ch_names]) * 1e6
    f, P = sps.welch(d, fs=sfq, nperseg=int(4 * sfq), noverlap=int(2 * sfq), axis=-1)
    a = P[:, (f >= band[0]) & (f <= band[1])].mean(1)
    b = P[:, (f >= broad[0]) & (f <= broad[1])].mean(1)
    return float(np.median(np.log10(a / b)))


probe = {}
for code in ("S210", "S200"):
    b = next(bb for bb in blocks if bb["code"] == code)
    r = load_window(REF_SUBJECT, hdr, b["t0_s"] + 5.0, 12.0, dl)
    probe[code] = relative_alpha(r)
    print(f"  probe of the first {code} block ({b['t0_s'] + 5:.1f} s, 12 s): "
          f"median log10 relative occipital alpha = {probe[code]:+.4f}")
EC_CODE = max(probe, key=probe.get)
print(f"-> eyes-closed code inferred from the signal: {EC_CODE} "
      f"(label_source: algorithmic; TODO(confirm) against the dataset descriptor)")

ref_block = blocks[REF_BLOCK]
print(f"\nAnalysing block {REF_BLOCK} (code {ref_block['code']}, {ref_block['t0_s']:.3f} s):")
raw_ref = load_window(REF_SUBJECT, hdr, ref_block["t0_s"], DURATION_S, dl)
CH = max(OCCIPITAL, key=lambda c: np.var(raw_ref.get_data(picks=[c])[0]) if c in raw_ref.ch_names else -1)
rel = {c: None for c in OCCIPITAL}
sfq = float(raw_ref.info["sfreq"])
f_ref, P_ref = sps.welch(raw_ref.get_data(picks=OCCIPITAL)[:, :] * 1e6, fs=sfq, nperseg=int(4 * sfq),
                         noverlap=int(2 * sfq), axis=-1)
for i, c in enumerate(OCCIPITAL):
    rel[c] = float(P_ref[i, (f_ref >= 8) & (f_ref <= 13)].mean() /
                   P_ref[i, (f_ref >= 1) & (f_ref <= 40)].mean())
CH = max(rel, key=rel.get)
print(f"  channel selection (the shipped asset's criterion: largest ratio of mean PSD in 8-13 Hz to "
      f"mean PSD in 1-40 Hz):")
for c in sorted(rel, key=rel.get, reverse=True):
    print(f"      {c}: relative alpha {rel[c]:.4f}")
print(f"  -> {CH}")
x = raw_ref.get_data(picks=[CH])[0] * 1e6
x = x - x.mean()
print(f"  {DURATION_S:g} s of {CH} at {sfq:g} Hz: SD {x.std():.2f} uV, largest absolute sample "
      f"{np.abs(x).max():.1f} uV, peak-to-peak of a 1-s moving mean "
      f"{np.ptp(np.convolve(x, np.ones(int(sfq)) / int(sfq), 'valid')):.1f} uV")
print(f"  (w-burst-detector/envelope.json records sd_uv 21.68, max_abs_uv 113.2, slow_drift_ptp_uv 97.2 for "
      f"this same segment, and relative_alpha O2 3.1027 > Oz 2.9291 > O1 2.8725)")

f_w, P_w = sps.welch(x, fs=sfq, nperseg=int(2 * sfq), noverlap=int(sfq), window="hann")
m = (f_w >= 0.5) & (f_w <= 60.0)
fit = helpers_l1.fit_specparam(f_w[m], P_w[m])
print(f"specparam fit over {helpers_l1.FIT_RANGE[0]:g}-{helpers_l1.FIT_RANGE[1]:g} Hz "
      f"({fit['tool']}, settings {fit['settings']}):")
print(f"    aperiodic offset {fit['offset']:.4f}, exponent {fit['exponent']:.4f}, "
      f"r^2 {fit['r_squared']:.4f}, error {fit['error']:.4f}")
print(f"    peaks (centre Hz, power, bandwidth Hz):")
for cf, pw_, bw_ in fit["peaks"]:
    print(f"        {cf:7.3f}  {pw_:.4f}  {bw_:.3f}")
print(f"    individual alpha frequency in {helpers_l1.IAF_RANGE[0]:g}-{helpers_l1.IAF_RANGE[1]:g} Hz: "
      f"{fit['iaf']:.4f} Hz" if np.isfinite(fit["iaf"]) else "    NO alpha peak fitted")

fig, ax = plt.subplots(figsize=(8.6, 4.3))
ax.semilogy(f_w[m], P_w[m], lw=1.2, color="k", label=f"{CH}, Welch 2-s Hann, 50 % overlap")
ax.semilogy(fit["freqs"] if "freqs" in fit else f_w[m][(f_w[m] >= 1) & (f_w[m] <= 40)],
            10 ** fit["aperiodic_log"], lw=1.4, color="tab:orange",
            label=f"aperiodic fit (exponent {fit['exponent']:.2f})")
ax.semilogy(f_w[m][(f_w[m] >= 1) & (f_w[m] <= 40)], 10 ** fit["model_log"], lw=1.2, color="tab:blue",
            ls="--", label=f"full model (r^2 {fit['r_squared']:.3f})")
if np.isfinite(fit["iaf"]):
    ax.axvline(fit["iaf"], color="tab:red", lw=1.0, ls=":", label=f"peak {fit['iaf']:.3f} Hz")
ax.axvline(L4.DATASETS_L4["ds-lemon"]["mains_hz"], color="0.6", lw=0.8, ls=":",
           label=f"{L4.DATASETS_L4['ds-lemon']['mains_hz']} Hz mains (no notch at any stage)")
ax.set(xlabel="Frequency (Hz)", ylabel="Power spectral density (uV^2/Hz)",
       title=f"ds-lemon {REF_SUBJECT} {CH}, {DURATION_S:g} s eyes closed: is there a peak? (uV^2/Hz)")
ax.legend(fontsize=8)
ax.grid(alpha=0.3, which="both")
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

IAF = float(fit["iaf"])
BAND = (8.0, 13.0)
print(f"\nIndividualised band, peak +/- 2 Hz: {IAF - 2:.2f}-{IAF + 2:.2f} Hz. The +/- 2 Hz is a CONVENTION, "
      f"not a measurement from this recording; the analysis below uses the conventional "
      f"{BAND[0]:g}-{BAND[1]:g} Hz band so that its numbers can be compared with the widget's.")

fh = L4.filter_hilbert(x, sfq, BAND)
env = fh["envelope"]
E = fh["edge_samples"]
core = env[E:len(env) - E]
print(f"band-pass {BAND[0]:g}-{BAND[1]:g} Hz: {fh['numtaps']} taps "
      f"({fh['numtaps'] / sfq:.3f} s); contaminated edge {E} samples = {fh['edge_s']:.3f} s each end, "
      f"so {len(core) / sfq:.3f} s of the {DURATION_S:g} s is analysed")
print(f"envelope percentiles over the ANALYSED region (uV): "
      + ", ".join(f"p{p}={np.percentile(core, p):.3f}" for p in (10, 25, 50, 75, 90, 95, 99)))
print(f"envelope percentiles over the WHOLE trace (uV):    "
      + ", ".join(f"p{p}={np.percentile(env, p):.3f}" for p in (10, 25, 50, 75, 90, 95, 99)))
print("  (w-burst-detector/envelope.json records the whole-trace values p50 6.841, p75 10.373, p95 16.801, "
      "mean 7.762, sd 4.937; the analysed-region values are lower because the contaminated edges are "
      "the quietest part of the trace)")

MIN_CYCLES = 3.0
res = {}
rows = []
for pct in (50, 75, 95):
    r = L4.detect_bursts(env, sfq, percentile=pct, min_cycles=MIN_CYCLES, band=BAND, edge_samples=E)
    res[pct] = r
    rows.append({"threshold": f"{pct}th percentile", "uV": r["threshold_uv"], "bursts": r["n_bursts"],
                 "rate (/min)": r["rate_per_min"], "mean duration (ms)": r["mean_duration_s"] * 1000,
                 "median duration (ms)": r["median_duration_s"] * 1000,
                 "duty cycle (%)": 100 * r["duty_cycle"]})
print()
print(f"Bursts in {BAND[0]:g}-{BAND[1]:g} Hz, minimum {MIN_CYCLES:g} cycles of the band centre "
      f"({res[50]['min_duration_s'] * 1000:.1f} ms = {res[50]['min_samples']} samples, ceil), "
      f"edges excluded, {res[50]['analysed_s']:.3f} s analysed:")
print(L4.fmt_table(rows, floatfmt="{:.2f}"))
print()
print(f"The signal is identical in all three rows. The burst rate falls "
      f"{res[50]['n_bursts'] / max(res[95]['n_bursts'], 1):.0f}x and the duty cycle "
      f"{res[50]['duty_cycle'] / max(res[95]['duty_cycle'], 1e-9):.0f}x on the threshold alone.")

fig, axes = plt.subplots(3, 1, figsize=(13.5, 9.0), sharex=True)
for ax, pct in zip(axes, (50, 75, 95)):
    L4.plot_envelope_bursts(fh["filtered"], env, sfq, res[pct], ax=ax, window_s=(5.0, 25.0),
                            title=f"ds-lemon {REF_SUBJECT} {CH}, {BAND[0]:g}-{BAND[1]:g} Hz, threshold at the "
                                  f"{pct}th percentile ({res[pct]['threshold_uv']:.2f} uV): "
                                  f"{res[pct]['n_bursts']} bursts, duty cycle "
                                  f"{100 * res[pct]['duty_cycle']:.1f} % (uV)")
axes[-1].set_xlabel("Time from the start of the 60-s block (s)")
fig.suptitle("The same twenty seconds under three thresholds (uV)", y=1.01)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

conv = []
for label, kwargs in (
        ("edges excluded, ceil (this notebook and the widget)", dict(edge_samples=E, round_mode="ceil")),
        ("edges excluded, round (the shipped sidecar's rule)", dict(edge_samples=E, round_mode="round")),
        ("whole trace, ceil", dict(edge_samples=0, round_mode="ceil")),
        ("whole trace, round", dict(edge_samples=0, round_mode="round"))):
    row = {"convention": label}
    for pct in (50, 75, 95):
        r = L4.detect_bursts(env, sfq, percentile=pct, min_cycles=MIN_CYCLES, band=BAND, **kwargs)
        row[f"p{pct} bursts"] = r["n_bursts"]
        row[f"p{pct} duty (%)"] = 100 * r["duty_cycle"]
    conv.append(row)
print("How much the two conventions move the answer, same signal, same band:")
print(L4.fmt_table(conv, floatfmt="{:.1f}"))
print()
r_side = L4.detect_bursts(env, sfq, median_multiple=1.5, min_duration_s=0.284, edge_samples=0,
                          round_mode="round")
print(f"The shipped sidecar's own reference row (1.5 x the whole-trace median = "
      f"{r_side['threshold_uv']:.3f} uV, 0.284 s minimum, whole trace, round): "
      f"{r_side['n_bursts']} bursts, mean duration {r_side['mean_duration_s'] * 1000:.1f} ms, duty cycle "
      f"{100 * r_side['duty_cycle']:.2f} %.")
print("  w-burst-detector/envelope.json records 20 bursts, 572.0 ms and 19.07 % for that row; "
      "site/notes/integration-phase3.md records that the widget's ceil rule gives 19 where the sidecar's "
      "round rule gives 20, on one extra sample. This notebook reproduces both.")

def cycle_features_fallback(sig_bp, sig_raw, fs, burst_segments):
    """Rise/decay and peak/trough symmetry per cycle, segmented at rising zero crossings of `sig_bp`.

    Documented fallback for bycycle: for each pair of consecutive rising zero crossings, the peak is the
    largest sample and the trough the smallest inside it; rise time is peak minus the preceding trough,
    decay time the following trough minus the peak.  Only cycles whose whole extent lies inside one of
    `burst_segments` are returned, since the ratios are meaningless outside a burst.
    """
    zc = np.where((sig_bp[:-1] <= 0) & (sig_bp[1:] > 0))[0] + 1
    out = []
    for a, b in zip(zc[:-1], zc[1:]):
        if not any(s <= a and b <= e for s, e in burst_segments):
            continue
        seg_ = sig_raw[a:b]
        if seg_.size < 4:
            continue
        pk = int(np.argmax(seg_))
        tr = int(np.argmin(seg_))
        if pk == 0 or tr == 0 or pk == seg_.size - 1 or tr == seg_.size - 1:
            continue
        rise = abs(pk - tr)
        decay = seg_.size - 1 - max(pk, tr) + min(pk, tr)
        if rise + decay == 0:
            continue
        out.append({"period_s": (b - a) / fs, "time_rdsym": rise / (rise + decay),
                    "time_ptsym": (pk if pk < tr else seg_.size - pk) / seg_.size,
                    "volt_amp": float(seg_.max() - seg_.min())})
    return out


segs = res[75]["segments"]
if HAVE_BYCYCLE:
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        df = _bycycle_features(x, sfq, BAND, center_extrema="peak", burst_method="cycles",
                               threshold_kwargs=dict(amp_fraction_threshold=0.3,
                                                     amp_consistency_threshold=0.4,
                                                     period_consistency_threshold=0.5,
                                                     monotonicity_threshold=0.6, min_n_cycles=3))
    burst = df[df["is_burst"]]
    rdsym = burst["time_rdsym"].to_numpy(float)
    ptsym = burst["time_ptsym"].to_numpy(float)
    periods = burst["period"].to_numpy(float) / sfq
    amps = burst["volt_amp"].to_numpy(float)
    TOOL = f"bycycle {__import__('bycycle').__version__}, burst_method='cycles'"
    print(f"cycle-by-cycle via {TOOL}: {len(df)} cycles found, {len(burst)} of them in bursts "
          f"({100 * len(burst) / len(df):.1f} %)")
else:
    feats = cycle_features_fallback(fh["filtered"], x, sfq, segs)
    rdsym = np.array([f["time_rdsym"] for f in feats])
    ptsym = np.array([f["time_ptsym"] for f in feats])
    periods = np.array([f["period_s"] for f in feats])
    amps = np.array([f["volt_amp"] for f in feats])
    TOOL = "the documented NumPy fallback in this notebook (bycycle not installed)"
    print(f"cycle-by-cycle via {TOOL}: {len(feats)} cycles inside the 75th-percentile bursts")

fig, axes = plt.subplots(1, 3, figsize=(17, 3.8))
for ax, v, name in ((axes[0], rdsym, "rise-decay symmetry"), (axes[1], ptsym, "peak-trough symmetry")):
    ax.hist(v, bins=25, color="0.7")
    ax.axvline(0.5, color="tab:red", lw=1.4, ls="--", label="0.5 = perfectly sinusoidal")
    ax.axvline(np.median(v), color="tab:blue", lw=1.4, label=f"median {np.median(v):.4f}")
    ax.set(xlabel=f"{name} (dimensionless, 0-1)", ylabel="Cycles",
           title=f"{name}\n{CH}, {BAND[0]:g}-{BAND[1]:g} Hz, {len(v)} cycles in bursts")
    ax.legend(fontsize=8)
axes[2].scatter(1 / periods, amps, s=8, alpha=0.6)
axes[2].axvline(IAF, color="tab:red", lw=1.0, ls=":", label=f"specparam peak {IAF:.2f} Hz")
axes[2].set(xlabel="Instantaneous cycle frequency (Hz)", ylabel="Cycle amplitude, peak to trough (uV)",
            title=f"Every burst cycle: frequency against amplitude\n{CH} (uV)")
axes[2].legend(fontsize=8)
axes[2].grid(alpha=0.3)
fig.suptitle(f"ds-lemon {REF_SUBJECT}, cycle-by-cycle waveform features ({TOOL})", y=1.04)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

from scipy import stats as sstats
for name, v in (("rise-decay symmetry", rdsym), ("peak-trough symmetry", ptsym)):
    t_, p_ = sstats.ttest_1samp(v, 0.5)
    print(f"  {name:22s}: median {np.median(v):.4f}, mean {v.mean():.4f} +/- {v.std(ddof=1):.4f}, "
          f"t({len(v) - 1}) = {t_:+.2f} against 0.5, p = {p_:.2e}")
print(f"  cycle frequency: median {np.median(1 / periods):.3f} Hz (specparam peak {IAF:.3f} Hz); "
      f"amplitude: median {np.median(amps):.2f} uV")
print("  A departure from 0.5 is a statement about waveform shape, not about coupling. It is also exactly "
      "what puts energy at harmonics of the fundamental, which is why L4.5's comodulogram of a "
      "non-sinusoidal rhythm lights up with no coupling present.")

SUBJECTS = ["sub-010002", "sub-010003", "sub-010004", "sub-010005", "sub-010006"]
PROBE_S = 12.0
cohort = []
for sid in SUBJECTS:
    sub_dl = L4.Downloads(f"{sid}", verbose=False).start() if False else None
    try:
        h = L4.lemon_headers(sid)
        bl = L4.lemon_blocks(h)
        if len(bl) < 4:
            print(f"  {sid}: only {len(bl)} blocks; skipped")
            continue
        scores = {}
        for code in sorted({b["code"] for b in bl}):
            b0 = next(bb for bb in bl if bb["code"] == code)
            rr = load_window(sid, h, b0["t0_s"] + 5.0, PROBE_S, dl)
            scores[code] = relative_alpha(rr) if rr is not None else -np.inf
        ec = max(scores, key=scores.get)
        # the first block of the eyes-closed code that has 60 s of room before the next block
        cand = [b for b in bl if b["code"] == ec and
                (b["next_block_t0_s"] is None or b["next_block_t0_s"] - b["t0_s"] >= DURATION_S)]
        b = cand[min(REF_BLOCK // 2, len(cand) - 1)] if cand else None
        if b is None:
            print(f"  {sid}: no {ec} block with {DURATION_S:g} s of room; skipped")
            continue
        rr = load_window(sid, h, b["t0_s"], DURATION_S, dl)
        if rr is None:
            continue
        sfq_i = float(rr.info["sfreq"])
        pk_ch, pk_rel = None, -np.inf
        f_i, P_i = sps.welch(rr.get_data(picks=[c for c in OCCIPITAL if c in rr.ch_names])[:, :] * 1e6,
                             fs=sfq_i, nperseg=int(4 * sfq_i), noverlap=int(2 * sfq_i), axis=-1)
        for k, c in enumerate([c for c in OCCIPITAL if c in rr.ch_names]):
            v = P_i[k, (f_i >= 8) & (f_i <= 13)].mean() / P_i[k, (f_i >= 1) & (f_i <= 40)].mean()
            if v > pk_rel:
                pk_ch, pk_rel = c, v
        xi = rr.get_data(picks=[pk_ch])[0] * 1e6
        xi = xi - xi.mean()
        fi, Pi = sps.welch(xi, fs=sfq_i, nperseg=int(2 * sfq_i), noverlap=int(sfq_i), window="hann")
        mi = (fi >= 0.5) & (fi <= 60.0)
        fit_i = helpers_l1.fit_specparam(fi[mi], Pi[mi])
        fhi = L4.filter_hilbert(xi, sfq_i, BAND)
        row = {"subject": sid, "eyes-closed code": ec, "block t0 (s)": b["t0_s"], "channel": pk_ch,
               "relative alpha": float(pk_rel), "SD (uV)": float(xi.std()),
               "exponent": fit_i["exponent"], "r^2": fit_i["r_squared"],
               "alpha peak (Hz)": fit_i["iaf"] if np.isfinite(fit_i["iaf"]) else float("nan")}
        for pct in (50, 75, 95):
            ri = L4.detect_bursts(fhi["envelope"], sfq_i, percentile=pct, min_cycles=MIN_CYCLES, band=BAND,
                                  edge_samples=fhi["edge_samples"])
            row[f"p{pct} bursts"] = ri["n_bursts"]
            row[f"p{pct} duty (%)"] = 100 * ri["duty_cycle"]
        cohort.append(row)
        print(f"  {sid}: {ec} block at {b['t0_s']:.0f} s, {pk_ch}, peak "
              + (f"{fit_i['iaf']:.3f} Hz" if np.isfinite(fit_i["iaf"]) else "NONE in 7-13 Hz")
              + f", {row['p75 bursts']} bursts at p75 ({row['p75 duty (%)']:.1f} % duty)")
    finally:
        dl.cleanup()      # every window of this subject goes before the next subject starts
print()
print(L4.fmt_table(cohort, ["subject", "channel", "relative alpha", "SD (uV)", "exponent", "r^2",
                            "alpha peak (Hz)", "p50 bursts", "p75 bursts", "p95 bursts",
                            "p50 duty (%)", "p75 duty (%)", "p95 duty (%)"], floatfmt="{:.3f}"))
no_peak = [r["subject"] for r in cohort if not np.isfinite(r["alpha peak (Hz)"])]
print()
print(f"Subjects with NO fitted alpha peak in {helpers_l1.IAF_RANGE[0]:g}-{helpers_l1.IAF_RANGE[1]:g} Hz: "
      + (", ".join(no_peak) if no_peak else "none of them"))
print(f"Burst rate at the 75th percentile across {len(cohort)} subjects: median "
      f"{np.median([r['p75 bursts'] for r in cohort]):.0f}, range "
      f"{min(r['p75 bursts'] for r in cohort)}-{max(r['p75 bursts'] for r in cohort)}; duty cycle median "
      f"{np.median([r['p75 duty (%)'] for r in cohort]):.1f} %, range "
      f"{min(r['p75 duty (%)'] for r in cohort):.1f}-{max(r['p75 duty (%)'] for r in cohort):.1f} %.")
print("A percentile threshold is defined on each subject's own envelope, so it adapts to their amplitude; "
      "the spread that remains is a spread in burstiness, not in microvolts.")

try:
    KEY = {50: (33, 40.3), 75: (19, 17.8), 95: (3, 2.5)}   # site/notes/integration-phase3.md, widgets-K
    KEY_PEAK_HZ = 10.272
    print("nb-4-6-bursts -- L4.6 numbers (draft; TODO(confirm) at author review)")
    print(f"Data: ds-lemon {REF_SUBJECT}, block {REF_BLOCK} (code {ref_block['code']}, "
          f"{ref_block['t0_s']:.3f} s), {DURATION_S:g} s, channel {CH}, resampled "
          f"{hdr['sfreq']:g} -> {sfq:g} Hz with a {sfq * 0.4:g} Hz anti-alias low-pass; mean removed; "
          f"nothing else filtered. Licence {L4.DATASETS_L4['ds-lemon']['license']}")
    print(f"Fetched as a byte range: {DURATION_S * hdr['n_channels'] * 2 * hdr['sfreq'] / 1e6:.1f} MB of a "
          f"317 MB file, plus two {PROBE_S:g}-s probes to identify the eyes-closed marker code.")
    print(f"Eyes-closed code inferred from the signal: {EC_CODE} (relative occipital alpha "
          + ", ".join(f"{k} {v:+.3f}" for k, v in probe.items())
          + "). label_source: algorithmic; TODO(confirm) against the dataset descriptor.")
    print()
    print("ANSWER KEY -- ex-4-6: does this channel have a spectral peak, and what do bursts do with the "
          "threshold?")
    print(f"    Spectral peak: {'YES' if np.isfinite(IAF) else 'NO'}, at {IAF:.3f} Hz "
          f"(specparam {fit['tool']}, aperiodic exponent {fit['exponent']:.3f}, r^2 "
          f"{fit['r_squared']:.3f}). Band power in {BAND[0]:g}-{BAND[1]:g} Hz is therefore about a rhythm "
          f"for this channel, which is what makes it a fair subject for the rest of the lesson.")
    print(f"    Cross-check: w-burst-detector/envelope.json records {KEY_PEAK_HZ:g} Hz, exponent 1.114 and "
          f"r^2 0.975. This notebook "
          f"{'AGREES' if abs(IAF - KEY_PEAK_HZ) < 0.005 else 'DISAGREES'} "
          f"(difference {IAF - KEY_PEAK_HZ:+.5f} Hz).")
    print()
    print(f"    Bursts in {BAND[0]:g}-{BAND[1]:g} Hz, minimum {MIN_CYCLES:g} cycles of the band centre "
          f"({res[50]['min_duration_s'] * 1000:.1f} ms -> {res[50]['min_samples']} samples, ceil), edges "
          f"excluded ({res[50]['edge_s']:.3f} s each end, {res[50]['analysed_s']:.3f} s analysed):")
    ok = True
    for pct in (50, 75, 95):
        r = res[pct]
        kb, kd = KEY[pct]
        same = (r["n_bursts"] == kb) and abs(100 * r["duty_cycle"] - kd) < 0.05
        ok &= same
        print(f"        {pct}th percentile = {r['threshold_uv']:6.2f} uV -> {r['n_bursts']:3d} bursts, "
              f"{r['rate_per_min']:5.1f}/min, mean {r['mean_duration_s'] * 1000:5.0f} ms, duty "
              f"{100 * r['duty_cycle']:5.1f} %   [key {kb} / {kd} % -> {'agrees' if same else 'DIFFERS'}]")
    print(f"    Cross-check against site/notes/integration-phase3.md (widgets-K): "
          f"{'ALL THREE ROWS AGREE' if ok else 'AT LEAST ONE ROW DIFFERS -- report both'}.")
    print(f"    The lesson: the rate falls {res[50]['n_bursts'] / max(res[95]['n_bursts'], 1):.0f}x and the "
          f"duty cycle {res[50]['duty_cycle'] / max(res[95]['duty_cycle'], 1e-9):.0f}x on one free parameter, "
          f"with the signal untouched.")
    print()
    print("    The conventions are choices, and two of them move the count:")
    for c in conv:
        print(f"        {c['convention']:52s} "
              + "  ".join(f"p{p} {c[f'p{p} bursts']:3d} ({c[f'p{p} duty (%)']:.1f} %)" for p in (50, 75, 95)))
    print(f"        the shipped sidecar's own row (1.5 x whole-trace median, 0.284 s, round, whole trace): "
          f"{r_side['n_bursts']} bursts, {100 * r_side['duty_cycle']:.2f} % -- the sidecar records 20 and "
          f"19.07 %.")
    print()
    print(f"Supporting -- cycle-by-cycle ({TOOL}):")
    print(f"    rise-decay symmetry median {np.median(rdsym):.4f}, peak-trough symmetry median "
          f"{np.median(ptsym):.4f} over {len(rdsym)} burst cycles (0.5 would be sinusoidal); median cycle "
          f"frequency {np.median(1 / periods):.3f} Hz, median amplitude {np.median(amps):.2f} uV.")
    print()
    print(f"Supporting -- burst rates per subject ({len(cohort)} LEMON subjects, one 60-s eyes-closed block "
          f"each, one occipital channel each chosen by relative alpha):")
    for r in cohort:
        print(f"    {r['subject']} {r['channel']:>3s}: peak "
              + (f"{r['alpha peak (Hz)']:6.3f} Hz" if np.isfinite(r["alpha peak (Hz)"]) else "   NONE   ")
              + f", exponent {r['exponent']:.3f}, SD {r['SD (uV)']:5.2f} uV | bursts "
              + "/".join(f"{r[f'p{p} bursts']:3d}" for p in (50, 75, 95))
              + " | duty " + "/".join(f"{r[f'p{p} duty (%)']:5.1f}" for p in (50, 75, 95)) + " %")
    print(f"    Subjects with no fitted alpha peak in {helpers_l1.IAF_RANGE[0]:g}-"
          f"{helpers_l1.IAF_RANGE[1]:g} Hz: " + (", ".join(no_peak) if no_peak else "none of them"))
    print()
    print("Pitfalls: pf-band-power-slope, pf-harmonics-as-pac, pf-narrowband-filter-oscillation, "
          "pf-notch-hole-in-band (this recording has NO notch at any stage, so its 50 Hz mains line is "
          "intact and visible in section 2's spectrum -- the pitfall is what a notch would have done to a "
          "band that touches it). Widgets: w-burst-detector, w-aperiodic-explorer (mode compare).")
    print("TODO(confirm): what LEMON's marker codes mean (inferred here from relative occipital alpha); the "
          "exact dataset terms behind the CC BY 4.0 record (spec section 13 item 5); and the +/- 2 Hz "
          "individualised band, which is a convention rather than a measurement.")
finally:
    dl.finish()