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

# 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
import pooch

mne.set_log_level("WARNING")
pooch.get_logger().setLevel("WARNING")   # no download chatter: it would print local paths
plt.rcParams["figure.dpi"] = 72
print(f"MNE {mne.__version__}; helpers_l5 imported from notebooks/_shared")
print("mne-connectivity available:", L5.have_mne_connectivity())
print(L5.disk_line("disk at the start"))

print("Licences, from data/directory.yaml (never from memory):")
L5.print_licences("ds-lemon", notes=True)

SUBJECT = "sub-010002"
MAX_MINUTES = 3.2
BAND = L5.ALPHA_BAND

print(L5.disk_line("disk before the download"))
print()
for k, v in L5.LEMON_PIPELINE.items():
    print(f"  {k:16s}: {v}")
print()

paths = []
try:
    epochs, meta = L5.lemon_condition_epochs(SUBJECT, "eyes-closed", max_minutes=MAX_MINUTES)
    paths = L5.lemon_files(SUBJECT)
    if epochs is None:
        raise RuntimeError(f"ds-lemon could not be fetched: {meta.get('error')}")
    data_uv = epochs.get_data(copy=True) * 1e6          # (n_epochs, n_channels, n_times), uV
    ch_names = list(epochs.ch_names)
    sfreq = float(epochs.info["sfreq"])
    ep_info = epochs.info.copy()
finally:
    print()
    freed = L5.delete_files(paths)
    print(L5.disk_line("disk after the download was deleted"))

print()
print(f"blocks found in the first {MAX_MINUTES:g} min of the recording:")
for b in meta["blocks"]:
    print(f"  {b['condition']:12s} {b['code']}  {b['start_s']:7.1f}-{b['end_s']:7.1f} s  "
          f"({b['n_ticks']} ticks)")
a = meta["alpha_check"]
print()
print(f"falsifiable check on the S210 reading -- {BAND[0]:g}-{BAND[1]:g} Hz power at {a['channels']}:")
print(f"  eyes-closed {a['eyes_closed_uv2_hz']:.3f} uV^2/Hz, eyes-open {a['eyes_open_uv2_hz']:.3f} uV^2/Hz, "
      f"ratio {a['ratio']:.2f}")
print("  a ratio above 1 is consistent with S210 = eyes closed; it does not prove it, and the catalog does "
      "not say. TODO(confirm)")
print()
print(f"{len(data_uv)} eyes-closed epochs of {meta['epoch_s']:g} s x {len(ch_names)} channels at "
      f"{sfreq:g} Hz ({meta['n_dropped']} dropped above {meta['ptp_criterion_uv']:g} uV peak-to-peak)")

METHODS = ("coh", "imcoh", "plv", "pli", "wpli")
for m in METHODS:
    print(f"  {m:6s}: {L5.CONNECTIVITY_NOTES[m]}")
print()

con, info_lib = L5.connectivity(data_uv, sfreq=sfreq, methods=METHODS, fmin=BAND[0], fmax=BAND[1])
print(f"backend used for every number below: {info_lib['backend']}")

con_np, info_np = L5.connectivity(data_uv, sfreq=sfreq, methods=METHODS, fmin=BAND[0], fmax=BAND[1],
                                  backend="numpy")
if info_lib["backend"] != info_np["backend"]:
    print(f"cross-check against {info_np['backend']}:")
    for m in METHODS:
        print(f"  {m:6s}: largest disagreement over all "
              f"{len(ch_names) * (len(ch_names) - 1) // 2} pairs = "
              f"{np.abs(np.abs(con[m]) - np.abs(con_np[m])).max():.3e}")
    print("  the two paths are the same estimator -- one Hann-tapered DFT per epoch, averaged over epochs -- "
          "so this is")
    print("  machine precision rather than agreement by luck.  A fallback that gave different numbers would "
          "not be a fallback.")
else:
    print("mne-connectivity is not installed; the NumPy path in helpers_l5 produced every number below")
CON = {m: np.abs(con[m]) for m in METHODS}

fig, axes = plt.subplots(1, 5, figsize=(21, 4.2))
for ax, m in zip(axes, METHODS):
    L5.plot_matrix(CON[m], None, title=f"|{m}|  ({BAND[0]:g}-{BAND[1]:g} Hz)", ax=ax, vmin=0,
                   vmax=1.0 if m in ("coh", "plv") else 0.5, cbar_label="dimensionless (0 to 1)")
fig.suptitle(f"ds-lemon {SUBJECT}, eyes closed, {len(data_uv)} epochs x {len(ch_names)} channels: "
             f"five connectivity measures on the same alpha-band cross-spectra (dimensionless)",
             y=1.03, fontsize=11)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

print("Summary over all pairs (dimensionless):")
print(f"{'measure':>8s} {'median':>9s} {'mean':>9s} {'min':>9s} {'max':>9s}  strongest pair")
for m in METHODS:
    s = L5.matrix_summary(CON[m], m, labels=ch_names)
    print(f"{m:>8s} {s['median']:9.4f} {s['mean']:9.4f} {s['min']:9.4f} {s['max']:9.4f}  "
          f"{s['extreme_pair'][0]}-{s['extreme_pair'][1]}")

d_names, GC_CM = L5.great_circle_cm(ep_info, L5.HEAD_RADIUS_M)
assert d_names == ch_names
dist = L5.upper_pairs(GC_CM)
coh_v, wpli_v = L5.upper_pairs(CON["coh"]), L5.upper_pairs(CON["wpli"])

fig, axes = plt.subplots(1, 3, figsize=(16, 4.2))
axes[0].plot(coh_v, wpli_v, ".", ms=3, alpha=0.4)
axes[0].set_xlabel("coherence (dimensionless)"); axes[0].set_ylabel("wPLI (dimensionless)")
axes[0].set_title("every sensor pair, alpha band", fontsize=9); axes[0].grid(alpha=0.25)
axes[0].set_xlim(0, 1); axes[0].set_ylim(0, 1)
L5.plot_versus_distance(dist, {"coherence": coh_v, "wPLI": wpli_v},
                        title="against distance across the scalp", ax=axes[1],
                        ylabel="dimensionless (0 to 1)")
axes[1].set_ylim(0, 1)
axes[2].hist(coh_v, bins=40, alpha=0.6, label="coherence")
axes[2].hist(wpli_v, bins=40, alpha=0.6, label="wPLI")
axes[2].set_xlabel("value (dimensionless)"); axes[2].set_ylabel("number of sensor pairs")
axes[2].set_title("distribution over all pairs", fontsize=9); axes[2].legend(fontsize=8)
fig.suptitle(f"ds-lemon {SUBJECT} eyes closed, {BAND[0]:g}-{BAND[1]:g} Hz: coherence and wPLI on identical "
             f"epochs ({len(coh_v)} sensor pairs)", y=1.03, fontsize=11)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

r_cw = float(np.corrcoef(coh_v, wpli_v)[0, 1])
print(f"correlation between the two measures across the {coh_v.size} pairs: r = {r_cw:+.3f}")
print(f"median coherence {np.median(coh_v):.4f}, median wPLI {np.median(wpli_v):.4f}; "
      f"pairs above 0.5: coherence {int((coh_v > 0.5).sum())}, wPLI {int((wpli_v > 0.5).sum())}")
print(f"coherence falls with distance: median {np.median(coh_v[dist < 5]):.3f} under 5 cm, "
      f"{np.median(coh_v[dist > 15]):.3f} over 15 cm")
print(f"wPLI:                          median {np.median(wpli_v[dist < 5]):.3f} under 5 cm, "
      f"{np.median(wpli_v[dist > 15]):.3f} over 15 cm")

N_SURROGATE = 200
rng = np.random.default_rng(L5.SEED)


def surrogate_null(data, sfreq, methods, band, n=N_SURROGATE, rng=rng):
    """Max-statistic null: shuffle epochs per channel, recompute, keep the largest value over pairs."""
    n_ep, n_ch, _ = data.shape
    out = {m: np.empty(n) for m in methods}
    for k in range(n):
        shuffled = np.stack([data[rng.permutation(n_ep), c, :] for c in range(n_ch)], axis=1)
        c_k, _ = L5.connectivity(shuffled, sfreq=sfreq, methods=methods, fmin=band[0], fmax=band[1],
                                 backend="numpy")
        for m in methods:
            out[m][k] = np.abs(L5.upper_pairs(np.abs(c_k[m]))).max()
    return out


NULL_METHODS = METHODS
null = surrogate_null(data_uv, sfreq, NULL_METHODS, BAND)
print(f"{N_SURROGATE} epoch-shuffled surrogates; the null is the distribution of the LARGEST value over all "
      f"{coh_v.size} pairs,")
print("so a pair that beats the 95th percentile is significant with the multiple comparisons already paid for.")
print()
thr = {}
for m in NULL_METHODS:
    thr[m] = float(np.percentile(null[m], 95))
    obs = L5.upper_pairs(CON[m])
    print(f"  {m:6s}: null max-statistic median {np.median(null[m]):.4f}, 95th percentile {thr[m]:.4f}; "
          f"{int((obs > thr[m]).sum())} of {obs.size} pairs survive "
          f"({100 * (obs > thr[m]).mean():.1f} %)")
print()
print("Read the coherence row carefully.  This null destroys ALL dependence between channels, so what it tests "
      "is")
print("'is there any relationship at all', and volume conduction guarantees the answer is yes.  It is NOT a "
      "null for")
print("'connected beyond mixing', and no sensor-level surrogate can be: the mixing is in the data, not in the "
      "estimator.")

# The exercise: pairs that coherence calls connected and wPLI does not.
iu = np.triu_indices(len(ch_names), 1)
flag_coh = coh_v > thr["coh"]
flag_wpli = wpli_v > thr["wpli"]
only_coh = flag_coh & ~flag_wpli
print(f"ANSWER KEY -- ex-5-2 (which measure flags a connection that disappears under wPLI):")
print(f"  coherence flags {int(flag_coh.sum())} of {coh_v.size} pairs; wPLI flags {int(flag_wpli.sum())}.")
print(f"  {int(only_coh.sum())} pairs are flagged by coherence and NOT by wPLI "
      f"({100 * only_coh.mean():.1f} % of all pairs); "
      f"{int((flag_wpli & ~flag_coh).sum())} the other way round.")
print()
order = np.argsort(-(coh_v - wpli_v))
print(f"  the ten pairs where the two measures disagree most:")
print(f"  {'pair':>14s} {'distance (cm)':>14s} {'coherence':>10s} {'wPLI':>8s} {'|imcoh|':>9s} {'PLV':>7s}")
for k in order[:10]:
    i, j = iu[0][k], iu[1][k]
    print(f"  {ch_names[i] + '-' + ch_names[j]:>14s} {GC_CM[i, j]:14.1f} {CON['coh'][i, j]:10.4f} "
          f"{CON['wpli'][i, j]:8.4f} {CON['imcoh'][i, j]:9.4f} {CON['plv'][i, j]:7.4f}")
print()
print("  What it implies: a pair with high coherence and wPLI at the null is a pair whose relationship has no")
print("  measurable lag.  That is what one generator seen by two electrodes looks like (L5.1), and it is also")
print("  what a genuine zero-lag interaction would look like -- the measure cannot separate them, so the")
print("  honest reading is 'consistent with mixing', not 'no connection'.  wPLI buys robustness to mixing by")
print("  giving up the ability to see zero-lag coupling at all.")

# A single-electrode reference makes that electrode identically zero, and a channel with no power has no
# coherence with anything: the estimate is 0/0.  Every reference is therefore compared on the SAME reduced
# channel set -- the 59 channels that are not Cz, TP7 or TP8 -- so the numbers are comparable.
DROP = ["Cz", "TP7", "TP8"]
keep_idx = [i for i, c in enumerate(ch_names) if c not in DROP]
keep_names = [ch_names[i] for i in keep_idx]


def rereference(x, mode):
    if mode == "average":
        y = x - x.mean(axis=1, keepdims=True)
    elif mode == "Cz":
        y = x - x[:, [ch_names.index("Cz")], :]
    elif mode == "linked TP7/TP8":
        y = x - 0.5 * (x[:, [ch_names.index("TP7")], :] + x[:, [ch_names.index("TP8")], :])
    else:
        raise ValueError(mode)
    return y[:, keep_idx, :]


REFS = ["average", "Cz", "linked TP7/TP8"]
ref_con = {}
print(f"the same {len(keep_names)} channels under every reference "
      f"({', '.join(DROP)} dropped, because a reference electrode is identically zero)")
print()
print(f"{'reference':>18s} {'median coh':>12s} {'median wPLI':>13s} {'r with average ref (coh)':>26s} "
      f"{'(wPLI)':>10s}")
for mode in REFS:
    c, _ = L5.connectivity(rereference(data_uv, mode), sfreq=sfreq, methods=("coh", "wpli"),
                           fmin=BAND[0], fmax=BAND[1])
    ref_con[mode] = {k: np.abs(v) for k, v in c.items()}
    rc = float(np.corrcoef(L5.upper_pairs(ref_con[mode]["coh"]), L5.upper_pairs(ref_con["average"]["coh"]))[0, 1])
    rw = float(np.corrcoef(L5.upper_pairs(ref_con[mode]["wpli"]), L5.upper_pairs(ref_con["average"]["wpli"]))[0, 1])
    print(f"{mode:>18s} {np.median(L5.upper_pairs(ref_con[mode]['coh'])):12.4f} "
          f"{np.median(L5.upper_pairs(ref_con[mode]['wpli'])):13.4f} {rc:26.3f} {rw:10.3f}")
print()
print("A reference change cannot change a topography's shape (L3.4) but it can and does change every "
      "connectivity value,")
print("because it changes what each channel IS.  The ranking of pairs survives better than the values do -- "
      "which is why a")
print("paper that reports connectivity without naming its reference has not reported its result.")

flat = data_uv.transpose(1, 0, 2).reshape(len(ch_names), -1)
U, S, Vt = np.linalg.svd(flat - flat.mean(axis=1, keepdims=True), full_matrices=False)
var = S ** 2 / (S ** 2).sum()
print(f"leading principal component carries {100 * var[0]:.1f} % of the variance across "
      f"{len(ch_names)} channels; the first five carry {100 * var[:5].sum():.1f} %")

proj = np.eye(len(ch_names)) - np.outer(U[:, 0], U[:, 0])
cleaned = np.einsum("ij,ejt->eit", proj, data_uv)
c_clean, _ = L5.connectivity(cleaned, sfreq=sfreq, methods=("coh", "wpli"), fmin=BAND[0], fmax=BAND[1])
c_clean = {k: np.abs(v) for k, v in c_clean.items()}
print()
print(f"{'measure':>8s} {'median before':>15s} {'median after':>14s} {'max |change|':>14s} "
      f"{'r before vs after':>19s}")
for m in ("coh", "wpli"):
    b, a_ = L5.upper_pairs(CON[m]), L5.upper_pairs(c_clean[m])
    print(f"{m:>8s} {np.median(b):15.4f} {np.median(a_):14.4f} {np.abs(a_ - b).max():14.4f} "
          f"{float(np.corrcoef(b, a_)[0, 1]):19.3f}")
print()
print("Every pair moves, including pairs whose channels barely load on the removed component: the projection")
print("is applied to the whole channel set at once.  Reporting how many components were removed is therefore "
      "part of")
print("reporting a connectivity result, not a preprocessing footnote.")

def granger(X, a, b, sfreq=250.0, fmin=5.0, fmax=60.0):
    """Mean spectral Granger causality from channel a to channel b, or None when the library is absent."""
    if not L5.have_mne_connectivity():
        return None
    from mne_connectivity import spectral_connectivity_epochs
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        c = spectral_connectivity_epochs(X, method=["gc"], indices=([[a]], [[b]]), mode="multitaper",
                                         sfreq=sfreq, fmin=fmin, fmax=fmax, verbose=False)
    return float(np.asarray(c.get_data()).squeeze().mean())


SIM_FS, SIM_T, SIM_EP = 250.0, 500, 60
r = np.random.default_rng(L5.SEED)
sims = {}

X = np.zeros((SIM_EP, 2, SIM_T))                      # 1: x drives y at one sample's lag
for e in range(SIM_EP):
    x = r.standard_normal(SIM_T)
    y = np.zeros(SIM_T); y[1:] = 0.8 * x[:-1] + 0.3 * r.standard_normal(SIM_T - 1)
    X[e] = [x, y]
sims["A  x drives y, lag 1 sample, no shared dynamics"] = X

Z = np.zeros((SIM_EP, 2, SIM_T))                      # 2: both are the same oscillator, x drives y
for e in range(SIM_EP):
    x = np.zeros(SIM_T); y = np.zeros(SIM_T)
    ex, ey = r.standard_normal(SIM_T), r.standard_normal(SIM_T)
    for tt in range(2, SIM_T):
        x[tt] = 0.55 * x[tt - 1] - 0.8 * x[tt - 2] + ex[tt]
        y[tt] = 0.55 * y[tt - 1] - 0.8 * y[tt - 2] + 0.5 * x[tt - 1] + ey[tt]
    Z[e] = [x, y]
sims["B  x drives y, both strongly oscillatory"] = Z

for label, (g1, g2, n1, n2) in {
    "C  ONE source, zero lag, equal sensor noise": (1.0, 1.0, 0.2, 0.2),
    "D  ONE source, zero lag, UNEQUAL sensor noise": (1.0, 0.3, 0.2, 0.2),
}.items():
    W = np.zeros((SIM_EP, 2, SIM_T))
    for e in range(SIM_EP):
        s = mne.filter.filter_data(r.standard_normal(SIM_T), SIM_FS, BAND[0], BAND[1], verbose=False)
        s = s / s.std()
        W[e] = [g1 * s + n1 * r.standard_normal(SIM_T), g2 * s + n2 * r.standard_normal(SIM_T)]
    sims[label] = W

print(f"Spectral Granger causality, 5-60 Hz, {SIM_EP} epochs of {SIM_T / SIM_FS:g} s "
      f"(multitaper; mne-connectivity)")
print(f"{'simulation':>48s} {'GC 1->2':>9s} {'GC 2->1':>9s} {'coherence':>11s} {'wPLI':>8s}")
gc_rows = []
for label, Xs in sims.items():
    g12, g21 = granger(Xs, 0, 1, SIM_FS), granger(Xs, 1, 0, SIM_FS)
    c, _ = L5.connectivity(Xs, sfreq=SIM_FS, methods=("coh", "wpli"), fmin=BAND[0], fmax=BAND[1])
    row = (label, g12, g21, float(np.abs(c["coh"])[0, 1]), float(np.abs(c["wpli"])[0, 1]))
    gc_rows.append(row)
    if g12 is None:
        print(f"{label:>48s} {'n/a':>9s} {'n/a':>9s} {row[3]:11.4f} {row[4]:8.4f}")
    else:
        print(f"{label:>48s} {g12:9.4f} {g21:9.4f} {row[3]:11.4f} {row[4]:8.4f}")
print()
if gc_rows[0][1] is None:
    print("mne-connectivity is not installed, so the directed measures were not computed in this run; the "
          "symmetric")
    print("measures above came from the NumPy fallback and are unaffected.")
else:
    print("Read row D.  There is ONE generator, there is no lag, and Granger causality reports a strong and")
    print("asymmetric directed influence -- from the high-SNR sensor to the low-SNR one -- because unequal")
    print("measurement noise makes one channel's past a better predictor than the other's.  Row C shows the same")
    print("geometry with equal noise: symmetric, but not zero, which is the finite-sample floor.  wPLI is near")
    print("the floor for both, because neither has any lag to find.")

print("nb-5-2-connectivity -- L5.2 numbers (draft; TODO(confirm) at author review)")
print(f"Data: ds-lemon {SUBJECT}, eyes closed, first {MAX_MINUTES:g} min of the raw recording "
      f"(HTTP Range request, {freed:.0f} MB, deleted)")
print(f"      {len(data_uv)} epochs of {meta['epoch_s']:g} s x {len(ch_names)} EEG channels at {sfreq:g} Hz; "
      f"{L5.licence_line('ds-lemon')}")
print(f"Pipeline: {L5.LEMON_PIPELINE['filter_hz']} Hz FIR, resampled to {L5.LEMON_PIPELINE['resample_hz']:g} Hz, "
      f"{L5.LEMON_PIPELINE['reference']}")
print(f"Band: {BAND[0]:g}-{BAND[1]:g} Hz.  Backend: {info_lib['backend']}")
print()
print(f"{'measure':>8s} {'median':>9s} {'mean':>9s} {'max':>9s} {'pairs above the 95th pct of the surrogate null':>48s}")
for m in METHODS:
    s = L5.matrix_summary(CON[m], m)
    flagged = f"{int((L5.upper_pairs(CON[m]) > thr[m]).sum())} of {coh_v.size} (threshold {thr[m]:.4f})"
    print(f"{m:>8s} {s['median']:9.4f} {s['mean']:9.4f} {s['max']:9.4f} {flagged:>48s}")
print()
print(f"ANSWER KEY -- ex-5-2 (free response): coherence flags {int(flag_coh.sum())} pairs, wPLI flags "
      f"{int(flag_wpli.sum())}, and")
print(f"    {int(only_coh.sum())} pairs are flagged by coherence alone.  The two measures correlate at "
      f"r = {r_cw:+.3f} across pairs.")
print(f"    The largest single disagreement is "
      f"{ch_names[iu[0][order[0]]]}-{ch_names[iu[1][order[0]]]} at {GC_CM[iu[0][order[0]], iu[1][order[0]]]:.1f} cm: "
      f"coherence {coh_v[order[0]]:.4f}, wPLI {wpli_v[order[0]]:.4f}.")
print("    Rubric points: (i) the pair has no measurable lag; (ii) that is what volume conduction produces "
      "and also what")
print("    a true zero-lag interaction would produce, so the two cannot be separated at the sensors; (iii) "
      "the honest")
print("    statement is 'consistent with mixing', not 'not connected'; (iv) wPLI's robustness is bought by "
      "being blind to")
print("    zero-lag coupling entirely.")
print()
print("ANSWER KEY -- reference dependence (same epochs, same band, same estimator):")
for mode in REFS:
    print(f"    {mode:>18s}: median coherence "
          f"{np.median(L5.upper_pairs(ref_con[mode]['coh'])):.4f}, median wPLI "
          f"{np.median(L5.upper_pairs(ref_con[mode]['wpli'])):.4f}")
print("    TP7/TP8 stand in for mastoids; this montage has none.  TODO(confirm).")
print()
print("ANSWER KEY -- removing one linear component (the leading PC, "
      f"{100 * var[0]:.1f} % of variance):")
for m in ("coh", "wpli"):
    b, a_ = L5.upper_pairs(CON[m]), L5.upper_pairs(c_clean[m])
    print(f"    {m:>5s}: median {np.median(b):.4f} -> {np.median(a_):.4f}, largest single change "
          f"{np.abs(a_ - b).max():.4f}")
print()
if gc_rows[0][1] is not None:
    print("ANSWER KEY -- directed measures on simulations where the truth is known:")
    for label, g12, g21, ccoh, cw in gc_rows:
        print(f"    {label:>48s}: GC 1->2 {g12:.4f}, GC 2->1 {g21:.4f}, coherence {ccoh:.4f}, wPLI {cw:.4f}")
    print("    Row D is the warning: one generator, zero lag, unequal sensor noise, and a confident direction.")
print()
print("Pitfalls: pf-volume-conduction-connectivity, pf-reference-changes-everything.")
print(L5.disk_line("disk at the end"))