# 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
import helpers_l3 as L3
from scipy import signal as sps
from scipy import stats as sstats

# 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")

N_TRIALS = 40
rng = np.random.default_rng(20260918)
draws = np.abs(np.exp(1j * rng.uniform(-np.pi, np.pi, (20000, N_TRIALS))).mean(1))
print(f"ITC of {N_TRIALS} uniformly random phases")
print(f"  analytic expectation sqrt(pi)/(2 sqrt(n)) = {L4.itc_expectation(N_TRIALS):.4f}")
print(f"  analytic SD sqrt(1/n - pi/(4n))           = {L4.itc_null_sd(N_TRIALS):.4f}")
print(f"  simulated over 20000 draws                = {draws.mean():.4f} +/- {draws.std(ddof=1):.4f}")
print(f"  Rayleigh critical value at p = .05        = {L4.rayleigh_critical(N_TRIALS):.4f}")
print(f"  fraction of null draws above it           = {(draws > L4.rayleigh_critical(N_TRIALS)).mean():.4f} "
      f"(should be about 0.05)")
print()
print(f"{'n trials':>9s} {'null expectation':>17s} {'null SD':>9s} {'Rayleigh p=.05':>15s}")
for n in (10, 20, 40, 80, 160, 320):
    print(f"{n:9d} {L4.itc_expectation(n):17.4f} {L4.itc_null_sd(n):9.4f} {L4.rayleigh_critical(n):15.4f}")
print()
print("The null falls as 1/sqrt(n) and never reaches zero, so 'ITC was low' is not a statement until the "
      "trial count is given with it.")

# The synthetic pair the L4.5 exercise asks for, through the full pipeline rather than by formula:
# generate trials, Morlet transform, take the phase at the burst centre.
SYN = dict(sfreq=250.0, tmin=-1.0, tmax=1.5, f_hz=10.0, n_cycles=7.0, onset_s=0.3, amplitude_uv=20.0,
           noise_exponent=1.343, noise_amplitude_uv=19.8)
N_SEEDS = 200

def pipeline_itc(jitter_deg, seed, n_trials=N_TRIALS, evoked_uv=0.0):
    tr, tt = L4.synthetic_trials(n_trials, phase_jitter_deg=jitter_deg, evoked_amplitude_uv=evoked_uv,
                                 seed=seed, **SYN)
    r = L4.morlet_power(tr[:, None, :], sfreq=SYN["sfreq"], freqs=np.array([SYN["f_hz"]]),
                        n_cycles=np.array([SYN["n_cycles"]]), decim=1, output="complex", times=tt)
    z = r["complex"][:, 0, 0, int(np.argmin(np.abs(tt - SYN["onset_s"])))]
    return float(L4.itc(z))

vals = {}
for label, jit in (("random phase (360 deg full width)", 360.0), ("phase-locked (0 deg jitter)", 0.0)):
    v = np.array([pipeline_itc(jit, 20260918 + s) for s in range(N_SEEDS)])
    vals[label] = v
    print(f"{label:36s}: ITC = {v.mean():.4f} +/- {v.std(ddof=1):.4f} over {N_SEEDS} seeds "
          f"(5th-95th percentile {np.percentile(v, 5):.3f}-{np.percentile(v, 95):.3f})")
ideal = np.abs(np.exp(1j * np.zeros(N_TRIALS)).mean())
print(f"{'an ideal noiseless phasor':36s}: ITC = {ideal:.6f}")
print()
print("The phase-locked answer is NOT 1 because the wavelet sees the 1/f background as well as the burst; the "
      "random-phase answer is NOT 0 because the null of ITC is not zero. Both facts are the exercise.")

ERPCORE_SUBJECTS = [1, 2, 3]
dl = L4.Downloads("nb-4-5").start()
p3 = {}
for s in ERPCORE_SUBJECTS:
    L3.fetch_erpcore_subject(s, "P3", verbose=False)
    # A wider epoch than Level 3's -0.2..0.8 s, because a 3 Hz wavelet with 3 cycles is 1.6 s long and
    # MNE refuses a kernel longer than the signal.  The rejection criterion is still evaluated over
    # -0.2..0.8 s (helpers_l3.REJECT_WINDOW), so the same trials survive as in every Level-3 notebook.
    ep, nfo = L3.load_p3_epochs(s, tmin=-1.0, tmax=1.5, verbose=False)
    tgt = L3.condition_epochs(ep, "target")
    p3[s] = {"data": tgt.get_data(picks=[L3.P3_CHANNEL])[:, 0, :] * 1e6, "times": ep.times,
             "sfreq": float(ep.info["sfreq"]), "n": len(tgt), "info": nfo}
    print(f"  {L3.erpcore_subject_id(s)}: {len(tgt)} target trials kept at {L3.P3_CHANNEL}, "
          f"{ep.times[0]:g}..{ep.times[-1]:g} s, {ep.info['sfreq']:g} Hz")
print()
print("licence statements shipped with ds-erpcore, all three verbatim (helpers_l3.ERPCORE_LICENCE_STATEMENTS):")
for k, v in L3.ERPCORE_LICENCE_STATEMENTS.items():
    print(f"    {k}: {str(v)[:150]}")

ITC_FREQS = np.arange(3.0, 20.5, 0.5)
ITC_CYC = np.maximum(3.0, ITC_FREQS / 2.0)
fig, axes = plt.subplots(2, len(ERPCORE_SUBJECTS), figsize=(6 * len(ERPCORE_SUBJECTS), 7.6))
itc_rows = []
for j, s in enumerate(ERPCORE_SUBJECTS):
    d = p3[s]
    r = L4.morlet_power(d["data"][:, None, :], sfreq=d["sfreq"], freqs=ITC_FREQS, n_cycles=ITC_CYC,
                        decim=2, output="complex", times=d["times"])
    I = L4.itc(r["complex"][:, 0])
    crit = L4.rayleigh_critical(d["n"])
    L4.plot_tfr(I, ITC_FREQS, r["times"], ax=axes[0, j], vlim=(0.0, 0.8), cmap="magma", symmetric=False,
                edge_s=L4.edge_seconds(ITC_FREQS, ITC_CYC), cbar_label="ITC (0-1, dimensionless)",
                title=f"{L3.erpcore_subject_id(s)} ITC at {L3.P3_CHANNEL}, target trials (n = {d['n']})\n"
                      f"Rayleigh p=.05 critical value {crit:.3f}; null expectation "
                      f"{L4.itc_expectation(d['n']):.3f}")
    axes[1, j].plot(d["times"] * 1000, d["data"].mean(0), lw=1.5, color="k")
    axes[1, j].axvline(0, color="gray", lw=0.8)
    axes[1, j].axhline(0, color="gray", lw=0.6)
    axes[1, j].set(xlabel="Time from stimulus (ms)", ylabel="Amplitude (uV)",
                   title=f"{L3.erpcore_subject_id(s)}: the ERP it came from ({L3.P3_CHANNEL}, uV)")
    axes[1, j].grid(alpha=0.3)
    w = (r["times"] >= 0.3) & (r["times"] <= 0.6)
    low = (ITC_FREQS >= 3) & (ITC_FREQS <= 5)
    itc_rows.append({"subject": L3.erpcore_subject_id(s), "n target trials": d["n"],
                     "peak ITC (3-20 Hz, 0-0.8 s)": float(I[:, (r["times"] >= 0) & (r["times"] <= 0.8)].max()),
                     "ITC 3-5 Hz, 300-600 ms": float(I[np.ix_(low, np.where(w)[0])].mean()),
                     "null expectation": L4.itc_expectation(d["n"]),
                     "Rayleigh p=.05": L4.rayleigh_critical(d["n"])})
fig.suptitle("Inter-trial phase coherence on ds-erpcore P3 target trials (dimensionless, 0-1)", y=1.01)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline
print(L4.fmt_table(itc_rows, floatfmt="{:.3f}"))

import re

idx = L4.aszed_index()
km = L4.aszed_extract(idx, "keymaps.kmp").decode("utf-8", "replace")
print("--- the archive's own protocol table (keymaps.kmp):")
print(km[km.index("PROTOCOL:"):km.index("ANNOTATION:")].rstrip())
assr_files = sorted([e["name"] for e in idx["entries"]
                     if re.search(r"subset_[23]/subject_\d+/\d+/Phase 2\.edf$", e["name"])],
                    key=lambda n: (int(re.search(r"subject_(\d+)", n).group(1)), n))
print(f"\nPhase 2 ('Fixed Auditory Stimulus') recordings in the archive: {len(assr_files)}; "
      f"{sum(idx['by_name'][n]['bytes'] for n in assr_files) / 1e6:.1f} MB uncompressed in total.")
print(f"This notebook reads {min(30, len(assr_files))} of them "
      f"({sum(idx['by_name'][n]['bytes'] for n in assr_files[:30]) / 1e6:.1f} MB).")

N_FILES = 30
CH_A = "EEG Cz-LE"
tmpdir = Path(helpers_l1.download_dir()) / "aszed-tmp"
rows_a = []
paired = []
pooled_assr, pooled_sil = [], []
for name in assr_files[:N_FILES]:
    raw_a = L4.aszed_read_edf(idx, name, tmpdir, dl)
    if CH_A not in raw_a.ch_names:
        continue
    sfa = float(raw_a.info["sfreq"])
    xa = raw_a.get_data(picks=[CH_A])[0] * 1e6
    seg_n = int(round(sfa))                     # one-second segments, so the spectrum has 1 Hz bins
    A, S = [], []
    for onset, desc in zip(raw_a.annotations.onset, raw_a.annotations.description):
        i = int(round(onset * sfa))
        for k in range(2 if desc == "ASSR" else 1):
            piece = xa[i + k * seg_n: i + (k + 1) * seg_n]
            if piece.size == seg_n:
                (A if desc == "ASSR" else S).append(piece)
    if len(A) < 4:
        continue
    A, S = np.array(A), np.array(S)
    pooled_assr.append(A)
    pooled_sil.append(S)
    w = np.hanning(seg_n)
    F = np.fft.rfft((A - A.mean(1, keepdims=True)) * w, axis=1)
    fr = np.fft.rfftfreq(seg_n, 1 / sfa)
    P = (np.abs(F) ** 2).mean(0)
    i40 = int(np.argmin(np.abs(fr - 40.0)))
    nb = [k for k in range(len(fr)) if 30 <= fr[k] <= 50 and abs(fr[k] - 40) >= 3 and abs(fr[k] - 50) >= 3]
    PS_i = None
    if len(S) >= 2:
        FS_i = np.fft.rfft((S - S.mean(1, keepdims=True)) * w, axis=1)
        PS_i = (np.abs(FS_i) ** 2).mean(0)
    paired.append(10 * np.log10((P[i40] / np.median(P[nb])) / (PS_i[i40] / np.median(PS_i[nb])))
                  if PS_i is not None else np.nan)
    rows_a.append({"file": "/".join(name.split("/")[-3:-1]), "sfreq": sfa, "ASSR 1-s epochs": len(A),
                   "40 Hz power / neighbourhood median": float(P[i40] / np.median(P[nb])),
                   "ITC at 40 Hz": float(L4.itc(F[:, i40])),
                   "Rayleigh p=.05": L4.rayleigh_critical(len(A))})
snr = np.array([r["40 Hz power / neighbourhood median"] for r in rows_a])
it40 = np.array([r["ITC at 40 Hz"] for r in rows_a])
crit = np.array([r["Rayleigh p=.05"] for r in rows_a])
print(f"{len(rows_a)} recordings analysed at {CH_A} (device {L4.DATASETS_L4['ds-aszed']['device']})")
print(L4.fmt_table(rows_a[:8], floatfmt="{:.3f}"))
print("  ... (first 8 rows)")
print()
print(f"40 Hz power as a multiple of the 30-50 Hz neighbourhood median: median {np.median(snr):.3f}, "
      f"IQR {np.percentile(snr, 25):.3f}-{np.percentile(snr, 75):.3f}, max {snr.max():.2f}; "
      f"recordings with a 2x peak: {(snr > 2).sum()} of {len(snr)}")
print(f"ITC at 40 Hz: median {np.median(it40):.3f}, max {it40.max():.3f}; above the Rayleigh critical value "
      f"in {(it40 > crit).sum()} of {len(it40)} recordings")

# The decisive test: WITHIN each recording, is 40 Hz louder during ASSR than during Silence, over and above
# whatever its own 30-50 Hz neighbourhood does?  Paired, so every per-recording amplitude, montage and filter
# setting cancels.
pr = np.array([v for v in paired if np.isfinite(v)])
t_pair, p_pair = sstats.ttest_1samp(pr, 0.0)
w_pair, pw_pair = sstats.wilcoxon(pr)
print()
print(f"PAIRED within-recording test ({len(pr)} recordings): 40 Hz during ASSR minus 40 Hz during Silence, "
      f"each normalised by that condition's own 30-50 Hz neighbourhood:")
print(f"    mean {pr.mean():+.3f} dB, median {np.median(pr):+.3f} dB, SD {pr.std(ddof=1):.3f}; "
      f"t({len(pr) - 1}) = {t_pair:+.3f}, p = {p_pair:.4f}; Wilcoxon p = {pw_pair:.4f}; "
      f"{int((pr > 0).sum())} of {len(pr)} recordings above zero.")

PA = np.concatenate(pooled_assr)
PS = np.concatenate(pooled_sil)
n_bin = PA.shape[1]
w = np.hanning(n_bin)
fr = np.fft.rfftfreq(n_bin, 1 / 256.0)
FA = np.fft.rfft((PA - PA.mean(1, keepdims=True)) * w, axis=1)
FS = np.fft.rfft((PS - PS.mean(1, keepdims=True)) * w, axis=1)
spec_a = np.median(np.abs(FA) ** 2, axis=0)
spec_s = np.median(np.abs(FS) ** 2, axis=0)
itc_a = L4.itc(FA)
crit_pooled = L4.rayleigh_critical(len(PA))

fig, axes = plt.subplots(1, 2, figsize=(14, 4.2))
keep = (fr >= 2) & (fr <= 120)
axes[0].semilogy(fr[keep], spec_a[keep], lw=1.3, label=f"ASSR blocks (n = {len(PA)} 1-s epochs)")
axes[0].semilogy(fr[keep], spec_s[keep], lw=1.3, label=f"Silence blocks (n = {len(PS)})")
axes[0].axvline(40, color="tab:red", lw=1.2, ls="--", label="40 Hz, where the response should be")
axes[0].axvline(50, color="0.5", lw=1.0, ls=":", label="50 Hz mains")
axes[0].set(xlabel="Frequency (Hz)", ylabel="Power (uV^2 per 1-Hz bin, median over epochs)",
            title=f"ds-aszed 'Fixed Auditory Stimulus' at {CH_A}\nmedian spectrum, ASSR vs Silence (uV^2)")
axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3, which="both")
axes[1].plot(fr[keep], itc_a[keep], lw=1.2, color="k")
axes[1].axhline(crit_pooled, color="tab:red", lw=1.0, ls="--",
                label=f"Rayleigh p=.05 for n = {len(PA)}: {crit_pooled:.3f}")
axes[1].axhline(L4.itc_expectation(len(PA)), color="tab:blue", lw=1.0, ls=":",
                label=f"null expectation {L4.itc_expectation(len(PA)):.3f}")
axes[1].axvline(40, color="tab:red", lw=1.2, ls="--")
axes[1].set(xlabel="Frequency (Hz)", ylabel="ITC (0-1, dimensionless)",
            title=f"ITC of the pooled ASSR epochs at {CH_A}\n(dimensionless; a 40 Hz ASSR would spike here)")
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

i40 = int(np.argmin(np.abs(fr - 40)))
diff_db = 10 * np.log10(spec_a / spec_s)
nb_mask = np.array([(30 <= f_ <= 50) and abs(f_ - 40) >= 3 and abs(f_ - 50) >= 3 for f_ in fr])
print(f"Pooled over {len(rows_a)} recordings ({len(PA)} one-second ASSR epochs, {len(PS)} Silence epochs):")
print(f"  ASSR-minus-Silence power at 40 Hz: {diff_db[i40]:+.2f} dB")
print(f"  ... but at its NEIGHBOURS (30-50 Hz, excluding 40 and the mains): median "
      f"{np.median(diff_db[nb_mask]):+.2f} dB, IQR {np.percentile(diff_db[nb_mask], 25):+.2f} to "
      f"{np.percentile(diff_db[nb_mask], 75):+.2f} dB. The ASSR blocks are broadband louder than the "
      f"Silence blocks -- 40 Hz is not singled out, and 40 Hz minus its own neighbourhood is "
      f"{diff_db[i40] - np.median(diff_db[nb_mask]):+.2f} dB.")
print(f"  ITC at 40 Hz: {itc_a[i40]:.4f}  (null expectation {L4.itc_expectation(len(PA)):.4f}, "
      f"Rayleigh p=.05 critical {crit_pooled:.4f})")
band = (fr >= 2) & (fr <= 120)
k = int(np.argmax(itc_a[band]))
print(f"  Largest ITC anywhere in 2-120 Hz: {itc_a[band][k]:.4f} at {fr[band][k]:.0f} Hz")
print()
print("FINDING -- this notebook cannot detect a 40 Hz auditory steady-state response in ds-aszed's Phase 2.")
print(f"  THE DECISIVE TEST is the paired one above, because it cancels everything that differs between "
      f"recordings: within each recording, 40 Hz during ASSR against 40 Hz during Silence, each normalised "
      f"by its own neighbourhood. It gives {pr.mean():+.3f} dB, t({len(pr) - 1}) = {t_pair:+.2f}, "
      f"p = {p_pair:.3f}, with {int((pr > 0).sum())} of {len(pr)} recordings on either side of zero. That is "
      f"a null result, not a small effect.")
print(f"  The pooled ASSR-minus-Silence line above reads {diff_db[i40] - np.median(diff_db[nb_mask]):+.2f} dB "
      f"at 40 Hz relative to its neighbours and looks like weak evidence FOR a response. It is not: pooling "
      f"epochs across recordings that differ by orders of magnitude in amplitude is not a paired comparison, "
      f"and the paired version of the same question is flat. The two are printed together rather than the "
      f"inconvenient one being dropped.")
print(f"  Within ASSR blocks alone, 40 Hz sits BELOW its own 30-50 Hz neighbourhood (median ratio "
      f"{np.median(snr):.2f}) -- a dip rather than a peak. That is a property of both conditions, so it says "
      f"something about the device's filtering rather than about the stimulus.")
print("  The ITC measure agrees and adds a caveat of its own: pooling epochs across recordings assumes the "
      "stimulus phase is reproducible against the marker, and the markers here are quantised to 10 ms, which "
      "is 0.4 of a cycle at 40 Hz. So the pooled ITC alone would not settle it either; the paired power test "
      "does.")
print("  What is NOT claimed: that the dataset contains no ASSR. This is one channel (Cz), one derivation "
      f"(the device's own linked-ear reference), {len(rows_a)} recordings, and a 50 Hz mains line that is "
      "enormous in this montage. TODO(confirm) with the data descriptor what the modulation rate actually "
      "was, and whether a different channel or reference is intended.")
print("  Consequence for the lesson: L4.5 needs a different positive control. Section 2's ERP CORE P3 is a "
      "real one (measured above), and section 1's synthetic phasor is an exact one. Both are in this "
      "notebook; the 40 Hz claim is not supportable from these data as read here.")

sweep = []
for amp in (0.0, 2.5, 5.0, 7.5, 10.0, 15.0, 20.0):
    v = np.array([pipeline_itc(360.0, 777000 + s, evoked_uv=amp) for s in range(60)])
    # total power at the burst frequency, same trials
    tr, tt = L4.synthetic_trials(N_TRIALS, phase_jitter_deg=360.0, evoked_amplitude_uv=amp, seed=777000,
                                 **SYN)
    r = L4.morlet_power(tr[:, None, :], sfreq=SYN["sfreq"], freqs=np.array([SYN["f_hz"]]),
                        n_cycles=np.array([SYN["n_cycles"]]), decim=1, times=tt)
    at = int(np.argmin(np.abs(tt - SYN["onset_s"])))
    pre = (tt > -0.9) & (tt < -0.4)
    tot = r["power"][:, 0, 0]
    sweep.append({"added evoked (uV)": amp, "ITC mean": float(v.mean()), "ITC SD": float(v.std(ddof=1)),
                  "total power at the burst / pre-burst": float(tot[:, at].mean() / tot[:, pre].mean())})
print(f"Additive evoked component on trials whose oscillation is NOT phase-locked (jitter held at 360 deg), "
      f"{N_TRIALS} trials, 60 seeds each:")
print(L4.fmt_table(sweep, floatfmt="{:.4f}"))

fig, ax = plt.subplots(figsize=(8.6, 4.2))
a = np.array([s["added evoked (uV)"] for s in sweep])
m = np.array([s["ITC mean"] for s in sweep])
sd = np.array([s["ITC SD"] for s in sweep])
ax.errorbar(a, m, yerr=sd, fmt="-o", capsize=3, label="ITC (mean +/- SD over 60 seeds)")
ax.axhline(L4.itc_expectation(N_TRIALS), color="tab:blue", lw=1.0, ls=":",
           label=f"null expectation for n = {N_TRIALS}: {L4.itc_expectation(N_TRIALS):.3f}")
ax.axhline(L4.rayleigh_critical(N_TRIALS), color="tab:red", lw=1.0, ls="--",
           label=f"Rayleigh p=.05: {L4.rayleigh_critical(N_TRIALS):.3f}")
ax.set(xlabel="Amplitude of the added evoked component (uV)", ylabel="ITC (0-1, dimensionless)",
       title="SYNTHETIC: ITC rises with a purely additive component and no phase reset at all\n"
             "(dimensionless; the oscillation's phase jitter is held at 360 degrees throughout)")
ax.legend(fontsize=8)
ax.grid(alpha=0.3)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

raw_ec = L4.load_imagery_raw("S001", [2], dl)          # R02 = eyes-closed baseline, ~1 min
sf_ec = float(raw_ec.info["sfreq"])
x_ec = raw_ec.get_data(picks=["O1"])[0] * 1e6
print(f"ds-eegbci S001 R02 ({helpers.EEGBCI_RUNS['R02']}), channel O1: {x_ec.size} samples at {sf_ec:g} Hz "
      f"= {x_ec.size / sf_ec:.1f} s; average reference; mains {L4.DATASETS_L4['ds-eegbci']['mains_hz']} Hz")

PH_F = np.arange(4.0, 16.1, 1.0)
AM_F = np.arange(20.0, 50.1, 2.5)
MI = np.zeros((len(PH_F), len(AM_F)))
for i, pf in enumerate(PH_F):
    ph = np.angle(sps.hilbert(L4.filter_hilbert(x_ec, sf_ec, (pf - 1.0, pf + 1.0))["filtered"]))
    for j, af in enumerate(AM_F):
        half = max(pf + 1.0, 4.0)
        am = L4.filter_hilbert(x_ec, sf_ec, (af - half, af + half))["envelope"]
        MI[i, j] = L4.modulation_index(ph, am)["mi"]
fig, axes = plt.subplots(1, 2, figsize=(14, 4.2))
L4.plot_comodulogram(MI, PH_F, AM_F, ax=axes[0],
                     title="ds-eegbci S001 R02 (eyes closed), O1\nTort modulation index (dimensionless)")
i_max, j_max = np.unravel_index(int(np.argmax(MI)), MI.shape)
pf, af = PH_F[i_max], AM_F[j_max]
ph = np.angle(sps.hilbert(L4.filter_hilbert(x_ec, sf_ec, (pf - 1.0, pf + 1.0))["filtered"]))
am = L4.filter_hilbert(x_ec, sf_ec, (af - max(pf + 1.0, 4.0), af + max(pf + 1.0, 4.0)))["envelope"]
test = L4.pac_surrogate_test(ph, am, n_surrogates=200)
axes[1].hist(test["null"], bins=30, color="0.7", label=f"{test['n_surrogates']} block-shift surrogates")
axes[1].axvline(test["mi"], color="tab:red", lw=2,
                label=f"observed MI {test['mi']:.5f} (z = {test['z']:.2f}, p = {test['p_empirical']:.4f})")
axes[1].set(xlabel="Modulation index (dimensionless)", ylabel="Surrogate count",
            title=f"The strongest cell ({pf:g} Hz phase, {af:g} Hz amplitude)\nagainst its own surrogates "
                  f"(dimensionless)")
axes[1].legend(fontsize=8)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline
print(f"Strongest cell: {pf:g} Hz phase / {af:g} Hz amplitude, MI = {test['mi']:.5f}")
print(f"  surrogates: mean {test['surrogate_mean']:.5f}, SD {test['surrogate_sd']:.5f} -> z = {test['z']:.2f}, "
      f"empirical p = {test['p_empirical']:.4f}")
print(f"  method: {test['method']}")
print(f"  The raw MI is {test['mi']:.5f}. On its own that number says nothing at all: the surrogate mean is "
      f"{test['surrogate_mean']:.5f} for the SAME signals with the timing broken.")

# The confound: one non-sinusoidal rhythm, no coupling of any kind, seeded.
rng = np.random.default_rng(20260918)
SF_C, DUR_C = 250.0, 120.0
t_c = np.arange(int(SF_C * DUR_C)) / SF_C
sharp = sps.sawtooth(2 * np.pi * 10.0 * t_c, width=0.15) * 15.0      # 10 Hz, strongly non-sinusoidal
sine = 15.0 * np.cos(2 * np.pi * 10.0 * t_c)                          # 10 Hz, sinusoidal
noise_c = L4.pink_noise(t_c.size, SF_C, 1.343, 10.0, rng)

fig, axes = plt.subplots(1, 3, figsize=(17.5, 3.9))
axes[0].plot(t_c[:int(0.6 * SF_C)], (sharp + noise_c)[:int(0.6 * SF_C)], lw=1.1, label="non-sinusoidal 10 Hz")
axes[0].plot(t_c[:int(0.6 * SF_C)], (sine + noise_c)[:int(0.6 * SF_C)], lw=1.1, label="sinusoidal 10 Hz")
axes[0].set(xlabel="Time (s)", ylabel="Amplitude (uV)",
            title="SYNTHETIC: two 10 Hz rhythms, same amplitude (uV)\nneither is coupled to anything")
axes[0].legend(fontsize=8)
axes[0].grid(alpha=0.3)
res_c = {}
for ax, (nm, sig) in zip(axes[1:], (("non-sinusoidal", sharp + noise_c), ("sinusoidal", sine + noise_c))):
    M = np.zeros((len(PH_F), len(AM_F)))
    for i, pf_ in enumerate(PH_F):
        ph_ = np.angle(sps.hilbert(L4.filter_hilbert(sig, SF_C, (pf_ - 1.0, pf_ + 1.0))["filtered"]))
        for j, af_ in enumerate(AM_F):
            h = max(pf_ + 1.0, 4.0)
            am_ = L4.filter_hilbert(sig, SF_C, (af_ - h, af_ + h))["envelope"]
            M[i, j] = L4.modulation_index(ph_, am_)["mi"]
    i_, j_ = np.unravel_index(int(np.argmax(M)), M.shape)
    ph_ = np.angle(sps.hilbert(L4.filter_hilbert(sig, SF_C, (PH_F[i_] - 1.0, PH_F[i_] + 1.0))["filtered"]))
    h = max(PH_F[i_] + 1.0, 4.0)
    am_ = L4.filter_hilbert(sig, SF_C, (AM_F[j_] - h, AM_F[j_] + h))["envelope"]
    tt_ = L4.pac_surrogate_test(ph_, am_, n_surrogates=200)
    res_c[nm] = (float(M.max()), PH_F[i_], AM_F[j_], tt_)
    L4.plot_comodulogram(M, PH_F, AM_F, ax=ax,
                         title=f"SYNTHETIC {nm} 10 Hz, NO coupling\nmax MI {M.max():.5f} at "
                               f"{PH_F[i_]:g}/{AM_F[j_]:g} Hz, z = {tt_['z']:.1f}, p = {tt_['p_empirical']:.4f}")
fig.suptitle("pf-harmonics-as-pac: waveform shape alone produces significant phase-amplitude coupling", y=1.04)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline
for nm, (mx, pf_, af_, tt_) in res_c.items():
    print(f"  {nm:15s}: max MI {mx:.5f} at {pf_:g} Hz phase / {af_:g} Hz amplitude, "
          f"z = {tt_['z']:.2f}, empirical p = {tt_['p_empirical']:.4f}")
print(f"  ratio of the two maxima: {res_c['non-sinusoidal'][0] / res_c['sinusoidal'][0]:.1f}x")
print("  Neither signal contains any coupling: each is one rhythm plus independent noise. Read the two "
      "columns separately.")
print(f"  THE MODULATION INDEX: the non-sinusoidal rhythm scores "
      f"{res_c['non-sinusoidal'][0] / res_c['sinusoidal'][0]:.0f}x the sinusoidal one. Any comparison that "
      f"puts a raw MI beside another raw MI -- between conditions, between groups, between subjects -- will "
      f"read that as coupling, and it is waveform shape.")
print(f"  THE SURROGATE TEST: it does not rescue the comparison either, and not in the direction one might "
      f"expect. Here the non-sinusoidal signal's own surrogates are inflated too (z = "
      f"{res_c['non-sinusoidal'][3]['z']:.2f}, p = {res_c['non-sinusoidal'][3]['p_empirical']:.4f}) while "
      f"the sinusoidal control's tiny MI comes out 'significant' (z = {res_c['sinusoidal'][3]['z']:.2f}, "
      f"p = {res_c['sinusoidal'][3]['p_empirical']:.4f}). A block shift of a strictly periodic signal by a "
      f"near-multiple of its own period reproduces the relation the test is trying to destroy, so the null "
      f"distribution is not the null the test assumes. Both readings are wrong in opposite directions.")
print("  The defence is to look at the waveform (L4.6's cycle-by-cycle analysis) and to report the shape "
      "measures beside the coupling measure, not to run a better test on the same two numbers.")

try:
    KEY_RANDOM, KEY_LOCKED = 0.140, 0.97       # site/notes/integration-phase3.md, data-p3a and widgets-I
    vr = vals["random phase (360 deg full width)"]
    vl = vals["phase-locked (0 deg jitter)"]
    print("nb-4-5-itc -- L4.5 exercise numbers (draft; TODO(confirm) at author review)")
    print(f"Datasets: ds-erpcore P3 (CC BY-SA 4.0 per data/directory.yaml, CONTESTED at source -- LICENSE "
          f"says BY-SA, dataset_description.json says CC0, the OSF node says BY; spec 10.7 makes the most "
          f"restrictive reading govern, so share-alike binds anything derived from it), subjects "
          f"{', '.join(L3.erpcore_subject_id(s) for s in ERPCORE_SUBJECTS)}; "
          f"ds-aszed (CC BY 4.0, Zenodo DOI {L4.DATASETS_L4['ds-aszed']['dataset_doi']}), "
          f"{len(rows_a)} Phase 2 recordings; ds-eegbci S001 R02 (ODC-By 1.0) for the coupling section.")
    print()
    print("ANSWER KEY -- ex-4-5 (numeric pair): ITC for 40 trials, uniformly random phase vs phase-locked")
    print(f"    random phase : {vr.mean():.4f}   (ask for the EXPECTED value, not one draw: the SD over "
          f"seeds is {vr.std(ddof=1):.4f} and the 5th-95th percentile is "
          f"{np.percentile(vr, 5):.3f}-{np.percentile(vr, 95):.3f})")
    print(f"    phase-locked : {vl.mean():.4f}   (SD {vl.std(ddof=1):.4f})")
    print(f"    The random-phase value is the analytic Rayleigh mean sqrt(pi)/(2 sqrt(40)) = "
          f"{L4.itc_expectation(N_TRIALS):.4f}, confirmed here by 20000 direct phase draws "
          f"({draws.mean():.4f} +/- {draws.std(ddof=1):.4f}).")
    print(f"    IT IS NOT ZERO and the phase-locked value IS NOT ONE. A tolerance reaching either round "
          f"number hides the point of the exercise. The Rayleigh critical value at p = .05 for 40 trials is "
          f"{L4.rayleigh_critical(N_TRIALS):.4f}, so 0.14 is not even close to significant.")
    print(f"    Cross-check: site/notes/integration-phase3.md records {KEY_RANDOM:g} and {KEY_LOCKED:g} from "
          f"data-p3a and widgets-I. This notebook "
          f"{'AGREES' if (abs(vr.mean() - KEY_RANDOM) < 0.02 and abs(vl.mean() - KEY_LOCKED) < 0.02) else 'DISAGREES'} "
          f"(differences {vr.mean() - KEY_RANDOM:+.4f} and {vl.mean() - KEY_LOCKED:+.4f}).")
    print()
    print("Supporting -- ITC on real trials (ds-erpcore P3 target trials at Pz):")
    for r in itc_rows:
        print(f"    {r['subject']}: n = {r['n target trials']:3d}, peak ITC "
              f"{r['peak ITC (3-20 Hz, 0-0.8 s)']:.3f}, 3-5 Hz over 300-600 ms "
              f"{r['ITC 3-5 Hz, 300-600 ms']:.3f}; null expectation {r['null expectation']:.3f}, "
              f"Rayleigh p=.05 {r['Rayleigh p=.05']:.3f}")
    print()
    print("FINDING -- the 40 Hz ASSR positive control spec section 6 L4.5 names does NOT work on these data.")
    print(f"    ds-aszed Phase 2 ('Fixed Auditory Stimulus', annotated ASSR/Silence in the files themselves), "
          f"{len(rows_a)} recordings at {CH_A}, 1-s epochs:")
    print(f"    40 Hz power / 30-50 Hz neighbourhood median: {np.median(snr):.3f} (IQR "
          f"{np.percentile(snr, 25):.3f}-{np.percentile(snr, 75):.3f}); a response would give a value above 1 "
          f"and this is a DIP, not a peak. {(snr > 2).sum()} of {len(snr)} recordings reach 2x.")
    print(f"    ITC at 40 Hz above the per-recording Rayleigh critical value: {(it40 > crit).sum()} of "
          f"{len(it40)}. Pooled over all epochs (n = {len(PA)}): ITC {itc_a[i40]:.4f} against a null "
          f"expectation of {L4.itc_expectation(len(PA)):.4f} and a critical value of {crit_pooled:.4f}.")
    print(f"    PAIRED within-recording test (the decisive one): 40 Hz during ASSR minus 40 Hz during "
          f"Silence, each normalised by its own 30-50 Hz neighbourhood, {pr.mean():+.3f} dB, "
          f"t({len(pr) - 1}) = {t_pair:+.3f}, p = {p_pair:.4f}, Wilcoxon p = {pw_pair:.4f}, "
          f"{int((pr > 0).sum())}/{len(pr)} recordings above zero.")
    print(f"    Pooled (NOT paired) ASSR-minus-Silence at 40 Hz: {diff_db[i40]:+.2f} dB against "
          f"{np.median(diff_db[nb_mask]):+.2f} dB at its neighbours. Reported because it points the other "
          f"way and a reader should see both; it is not a paired comparison and the paired one is flat.")
    print(f"    TODO(confirm): the archive's own protocol table calls this phase 'Fixed Auditory Stimulus' and "
          f"never states a modulation rate; 'assr-40hz' comes from data/directory.yaml. The rate, the "
          f"intended channel and the intended reference all need checking against the data descriptor before "
          f"the lesson claims a 40 Hz control. Until then L4.5's positive control should be the ERP CORE P3 "
          f"above (real) or the synthetic phasor of section 1 (exact).")
    print()
    print("Supporting -- additive evoked component with no phase reset (synthetic):")
    for s_ in sweep:
        print(f"    +{s_['added evoked (uV)']:4.1f} uV evoked -> ITC {s_['ITC mean']:.3f} "
              f"+/- {s_['ITC SD']:.3f}, total power at the burst x{s_['total power at the burst / pre-burst']:.2f}")
    print("    Phase jitter was held at 360 degrees throughout: nothing was reset. TODO(confirm): reading "
          "'total power rose' as evidence for the additive account is contested in the literature; the "
          "arithmetic here is exact, the inference is not.")
    print()
    print("Supporting -- phase-amplitude coupling:")
    print(f"    ds-eegbci S001 R02 eyes closed, O1: strongest cell {pf:g} Hz phase / {af:g} Hz amplitude, "
          f"MI {test['mi']:.5f}, z {test['z']:.2f}, empirical p {test['p_empirical']:.4f} against "
          f"{test['n_surrogates']} block-shift surrogates.")
    print(f"    pf-harmonics-as-pac: a non-sinusoidal 10 Hz rhythm with NO coupling gives MI "
          f"{res_c['non-sinusoidal'][0]:.5f}, "
          f"{res_c['non-sinusoidal'][0] / res_c['sinusoidal'][0]:.0f}x the sinusoidal control's "
          f"{res_c['sinusoidal'][0]:.5f}. Block-shift surrogates do not sort the two out: the "
          f"non-sinusoidal signal scores z {res_c['non-sinusoidal'][3]['z']:.2f} "
          f"(p {res_c['non-sinusoidal'][3]['p_empirical']:.4f}) and the sinusoidal control z "
          f"{res_c['sinusoidal'][3]['z']:.2f} (p {res_c['sinusoidal'][3]['p_empirical']:.4f}) -- shifting a "
          f"strictly periodic signal by a near-multiple of its own period reproduces the relation the "
          f"surrogate is meant to break. TODO(confirm) whether a surrogate better suited to periodic "
          f"signals (e.g. trial shuffling on epoched data) changes this; the point that MI alone cannot "
          f"separate waveform shape from coupling does not depend on it.")
    print()
    print("Pitfall: pf-harmonics-as-pac. Widget: w-phase-clock.")
finally:
    dl.finish()
    for s in ERPCORE_SUBJECTS:
        L3.delete_erpcore_subject(s, "P3", keep_small=True, verbose=False)
    print("ERP CORE heavy files (.set/.fdt) deleted; the small BIDS sidecars are left in the cache.")
    L4.free_disk_line("at the very end of nb-4-5")