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

L5.print_source_model_block()
print()
print("Licences of what IS used:")
L5.print_licences("ds-lemon")

import warnings
from scipy import stats

MONTAGE = ["Fp1", "Fp2", "F7", "F3", "Fz", "F4", "F8", "FC5", "FC1", "FC2", "FC6", "T7", "C3", "Cz", "C4",
           "T8", "CP5", "CP1", "CP2", "CP6", "AFz", "P7", "P3", "Pz", "P4", "P8", "PO9", "O1", "Oz", "O2",
           "PO10", "AF7", "AF3", "AF4", "AF8", "F5", "F1", "F2", "F6", "FT7", "FC3", "FC4", "FT8", "C5",
           "C1", "C2", "C6", "TP7", "CP3", "CPz", "CP4", "TP8", "P5", "P1", "P2", "P6", "PO7", "PO3",
           "POz", "PO4", "PO8"]
SFREQ, BAND, NOISE_UV, SNR = 250.0, L5.ALPHA_BAND, 1.0, 3.0
LAMBDA2 = 1.0 / SNR ** 2
N_REGIONS, MAX_DEPTH_MM = 12, 35.0

info = L5.sphere_info(MONTAGE, sfreq=SFREQ)
_raw = mne.io.RawArray(np.zeros((len(MONTAGE), 10)), info.copy(), verbose=False)
_raw.set_eeg_reference("average", projection=True, verbose=False)
info = _raw.info
bem = L5.sphere_models()["four-shell"]
rr = L5.source_grid(spacing_m=0.012, max_fraction=0.85, min_radius_m=0.012)
depth_mm = (L5.HEAD_RADIUS_M - np.linalg.norm(rr, axis=1)) * 1000
keep = depth_mm <= MAX_DEPTH_MM
rr, depth_mm = rr[keep], depth_mm[keep]
nn = rr / np.linalg.norm(rr, axis=1, keepdims=True)
fwd, G = L5.leadfield(info, bem, rr, fixed_normals=nn)
G = G - G.mean(axis=0, keepdims=True)                 # the measurement is average-referenced
print(f"{len(MONTAGE)} electrodes, {len(rr)} radial sources within {MAX_DEPTH_MM:.0f} mm of the scalp")

# Geometric regions: k-means on position, seeded, so the partition is reproducible and has no anatomy in it.
rng = np.random.default_rng(L5.SEED)
centres = rr[rng.choice(len(rr), N_REGIONS, replace=False)]
for _ in range(200):
    lab = np.argmin(((rr[:, None, :] - centres[None]) ** 2).sum(-1), axis=1)
    new = np.stack([rr[lab == k].mean(axis=0) if (lab == k).any() else centres[k]
                    for k in range(N_REGIONS)])
    if np.allclose(new, centres):
        break
    centres = new
labels = np.argmin(((rr[:, None, :] - centres[None]) ** 2).sum(-1), axis=1)


def region_name(c):
    """A coordinate description.  Deliberately NOT an anatomical label."""
    lr = "left" if c[0] < -0.015 else ("right" if c[0] > 0.015 else "midline")
    ap = "anterior" if c[1] > 0.015 else ("posterior" if c[1] < -0.015 else "central")
    iz = "superior" if c[2] > 0.02 else ("inferior" if c[2] < -0.005 else "mid")
    return f"{lr}-{ap}-{iz}"


REGIONS = [f"R{k + 1:02d} {region_name(centres[k])}" for k in range(N_REGIONS)]
print(f"{N_REGIONS} geometric regions (k-means on position, seed {L5.SEED}); "
      f"sizes {np.bincount(labels, minlength=N_REGIONS).tolist()}")
for k, name in enumerate(REGIONS):
    c = centres[k] * 1000
    print(f"  {name:28s} centre ({c[0]:+6.1f}, {c[1]:+6.1f}, {c[2]:+6.1f}) mm, "
          f"{int((labels == k).sum()):3d} sources")
print()
print("These names are coordinates.  No anatomical structure is named anywhere in this notebook, because a")
print("sphere has no anatomy and the template that would supply it does not clear the licence gate.")

from mne.minimum_norm import make_inverse_operator, make_inverse_resolution_matrix
from mne.minimum_norm.inverse import _assemble_kernel, prepare_inverse_operator

noise_cov = mne.make_ad_hoc_cov(info, std=NOISE_UV * 1e-6, verbose=False)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    inv = make_inverse_operator(info, fwd, noise_cov, loose=0.0, depth=None, fixed=True, verbose=False)
    R_full = make_inverse_resolution_matrix(fwd, inv, method="MNE", lambda2=LAMBDA2)
# The inverse kernel itself, so region time courses can be computed without ever forming a source-space
# time series (61 channels x 30000 samples stays small; 1000 sources x 30000 samples would not).
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    inv_prep = prepare_inverse_operator(inv, nave=1, lambda2=LAMBDA2, method="MNE", verbose=False)
    K, noise_norm, vertno, source_nn = _assemble_kernel(inv_prep, None, "MNE", pick_ori=None, verbose=False)
print(f"inverse kernel K: {K.shape} (sources x channels); resolution matrix R = K G: {R_full.shape}")
print(f"  check: max |K G - R| / max |R| = {np.abs(K @ G - R_full).max() / np.abs(R_full).max():.2e}")

# Region-level leakage: what fraction of region A's estimated power comes from region B's true activity.
A = np.abs(R_full) ** 2
LEAK = np.zeros((N_REGIONS, N_REGIONS))
for a in range(N_REGIONS):
    for b in range(N_REGIONS):
        LEAK[a, b] = A[np.ix_(labels == a, labels == b)].sum()
LEAK = LEAK / LEAK.sum(axis=1, keepdims=True)
print()
print(f"Leakage matrix: row A, column B = the share of region A's estimated power that a unit source in")
print(f"region B would contribute.  Rows sum to 1.")
print(f"  median self share (the diagonal): {np.median(np.diag(LEAK)):.3f}")
print(f"  median off-diagonal share:        {np.median(L5.upper_pairs(LEAK)):.3f}")
print(f"  largest off-diagonal share:       {L5.upper_pairs(LEAK).max():.3f} "
      f"({REGIONS[int(np.unravel_index(np.argmax(LEAK - np.diag(np.diag(LEAK))), LEAK.shape)[0])].split()[0]}"
      f"-{REGIONS[int(np.unravel_index(np.argmax(LEAK - np.diag(np.diag(LEAK))), LEAK.shape)[1])].split()[0]})")
print()
print(f"Only {100 * np.median(np.diag(LEAK)):.0f} % of a region's estimated power is its own.  Nothing has "
      f"been recorded yet: this is the")
print("operator, the geometry and the electrode montage, and it is the floor under every source-space")
print("connectivity value computed below.")

fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.6))
L5.plot_matrix(LEAK, [r.split()[0] for r in REGIONS], title="region-to-region leakage (rows sum to 1)",
               ax=axes[0], vmin=0, vmax=float(np.diag(LEAK).max()),
               cbar_label="share of estimated power (dimensionless)")
d_reg = np.linalg.norm(centres[:, None, :] - centres[None], axis=-1) * 1000
axes[1].plot(L5.upper_pairs(d_reg), L5.upper_pairs(LEAK), "o", ms=5)
axes[1].set_xlabel("distance between region centres (mm)")
axes[1].set_ylabel("leakage share (dimensionless)")
axes[1].set_title("leakage against distance", fontsize=9)
axes[1].grid(alpha=0.25)
fig.suptitle(f"Leakage of a minimum-norm operator, {len(MONTAGE)} electrodes, {N_REGIONS} geometric regions "
             f"(no anatomy)", y=1.02, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

def region_kernel(K, labels, data_cov):
    """One row of sensor weights per region: the first PC of the region's estimated source time courses."""
    W = np.zeros((N_REGIONS, K.shape[1]))
    for k in range(N_REGIONS):
        Kk = K[labels == k]
        M = Kk @ data_cov @ Kk.T
        w, V = np.linalg.eigh(M)
        u = V[:, -1]
        row = u @ Kk
        W[k] = row * np.sign(row[np.argmax(np.abs(row))])      # fix the sign so a PC flip cannot flip a lag
    return W


def sensor_epochs(source_amp, n_epochs=60, n_times=500, noise_uv=NOISE_UV, rng=None):
    """Epochs of sensor data from a dict {grid index: waveform generator}, plus sensor noise."""
    rng = rng or np.random.default_rng(L5.SEED)
    out = np.empty((n_epochs, len(MONTAGE), n_times))
    for e in range(n_epochs):
        s = np.zeros((len(rr), n_times))
        for j, gen in source_amp.items():
            s[j] = gen(rng)
        out[e] = G @ s + rng.standard_normal((len(MONTAGE), n_times)) * noise_uv * 1e-6
    return out


def narrowband(rng, n_times=500, band=BAND, sfreq=SFREQ):
    x = mne.filter.filter_data(rng.standard_normal(n_times + 400), sfreq, band[0], band[1], verbose=False)
    x = x[200:200 + n_times]
    return x / x.std()


N_BACKGROUND, BACKGROUND_SHARE = 24, 0.7
SRC_A = int(np.argmin(np.linalg.norm(rr - centres[0], axis=1)))
SRC_B = int(np.argmin(np.linalg.norm(rr - centres[min(6, N_REGIONS - 1)], axis=1)))
REG_A, REG_B = int(labels[SRC_A]), int(labels[SRC_B])
MOMENT = 50e-9
print(f"true sources: grid point {SRC_A} in {REGIONS[REG_A]} and grid point {SRC_B} in {REGIONS[REG_B]}, "
      f"{np.linalg.norm(rr[SRC_A] - rr[SRC_B]) * 1000:.0f} mm apart")
print(f"  B lags A by a quarter cycle; both {MOMENT * 1e9:.0f} nA.m; sensor noise {NOISE_UV:g} uV; "
      f"band {BAND[0]:g}-{BAND[1]:g} Hz")
print(f"  plus {N_BACKGROUND} INDEPENDENT background sources at random grid points, scaled together to "
      f"{BACKGROUND_SHARE:g} x the")
print("  coupled pair's sensor RMS.  Without them the simulation would have exactly two generators, and with")
print("  two generators every sensor pair's imaginary cross-spectrum is the same quantity times a real")
print("  constant -- so wPLI takes ONE value for the whole head and the comparison below would be vacuous")
print("  (nb-5-1 section 4 proves that).  A background is what makes the question non-degenerate.")


def coupled_pair(rng):
    a = narrowband(rng)
    quarter = np.real(__import__("scipy.signal", fromlist=["hilbert"]).hilbert(a) * np.exp(-1j * np.pi / 2))
    return a, quarter / quarter.std()


rng_sim = np.random.default_rng(L5.SEED)
bg_idx = rng_sim.choice([j for j in range(len(rr)) if j not in (SRC_A, SRC_B)], N_BACKGROUND,
                        replace=False)
X = np.empty((60, len(MONTAGE), 500))
X_pair = np.empty_like(X)
X_bg = np.empty_like(X)
for e in range(60):
    a, b = coupled_pair(rng_sim)
    s_pair = np.zeros((len(rr), 500))
    s_pair[SRC_A], s_pair[SRC_B] = MOMENT * a, MOMENT * b
    s_bg = np.zeros((len(rr), 500))
    for j in bg_idx:
        s_bg[j] = MOMENT * narrowband(rng_sim)
    X_pair[e], X_bg[e] = G @ s_pair, G @ s_bg
scale = BACKGROUND_SHARE * np.sqrt((X_pair ** 2).mean()) / np.sqrt((X_bg ** 2).mean())
for e in range(60):
    X[e] = X_pair[e] + scale * X_bg[e] + rng_sim.standard_normal((len(MONTAGE), 500)) * NOISE_UV * 1e-6

data_cov = np.mean([np.cov(x) for x in X], axis=0)
W = region_kernel(K, labels, data_cov)
TC = np.einsum("rc,ect->ert", W, X)              # (epochs, regions, times)
print(f"region time courses: {TC.shape} (epochs x regions x samples)")

con_src, backend = L5.connectivity(TC, sfreq=SFREQ, methods=("coh", "wpli"), fmin=BAND[0], fmax=BAND[1])
con_src = {k: np.abs(v) for k, v in con_src.items()}
print(f"backend: {backend['backend']}")
print()
print(f"the ONE true connection is {REGIONS[REG_A].split()[0]}-{REGIONS[REG_B].split()[0]}:")
print(f"  coherence {con_src['coh'][REG_A, REG_B]:.4f}, wPLI {con_src['wpli'][REG_A, REG_B]:.4f}")
others = [(a, b) for a in range(N_REGIONS) for b in range(a + 1, N_REGIONS) if {a, b} != {REG_A, REG_B}]
oc = np.array([con_src["coh"][a, b] for a, b in others])
ow = np.array([con_src["wpli"][a, b] for a, b in others])
print(f"  the other {len(others)} region pairs, which contain NO COUPLED source:")
print(f"    coherence median {np.median(oc):.4f}, max {oc.max():.4f} "
      f"({REGIONS[others[int(np.argmax(oc))][0]].split()[0]}-"
      f"{REGIONS[others[int(np.argmax(oc))][1]].split()[0]})")
print(f"    wPLI      median {np.median(ow):.4f}, max {ow.max():.4f} "
      f"({REGIONS[others[int(np.argmax(ow))][0]].split()[0]}-"
      f"{REGIONS[others[int(np.argmax(ow))][1]].split()[0]})")
print(f"  the true pair is ranked {int((oc >= con_src['coh'][REG_A, REG_B]).sum()) + 1} of "
      f"{len(others) + 1} by coherence and "
      f"{int((ow >= con_src['wpli'][REG_A, REG_B]).sum()) + 1} of {len(others) + 1} by wPLI")

flat = TC.transpose(1, 0, 2).reshape(N_REGIONS, -1)          # regions x (epochs*times)
Y = L5.symmetric_orthogonalise(flat)
print(f"symmetric orthogonalisation of {N_REGIONS} regions over {flat.shape[1]} samples:")
print(f"  largest |zero-lag correlation| between regions, before {np.abs(np.corrcoef(flat) - np.eye(N_REGIONS)).max():.4f}"
      f", after {np.abs(np.corrcoef(Y) - np.eye(N_REGIONS)).max():.2e}")
print(f"  each region keeps {np.median(np.abs(np.diag(np.corrcoef(flat, Y)[:N_REGIONS, N_REGIONS:]))):.3f} "
      f"of itself (median |correlation| with its own original)")

TC_orth = Y.reshape(N_REGIONS, TC.shape[0], TC.shape[2]).transpose(1, 0, 2)
con_orth, _ = L5.connectivity(TC_orth, sfreq=SFREQ, methods=("coh", "wpli"), fmin=BAND[0], fmax=BAND[1])
con_orth = {k: np.abs(v) for k, v in con_orth.items()}

ENV = L5.orthogonalised_envelope_correlation(TC, SFREQ, BAND, orthogonalise=False)
ENV_ORTH = L5.orthogonalised_envelope_correlation(TC, SFREQ, BAND, orthogonalise=True)

print()
print(f"{'measure':>44s} {'true pair':>10s} {'median other':>13s} {'max other':>10s} {'margin':>8s} "
      f"{'rank':>10s}")


def summarise(M, label):
    t = float(M[REG_A, REG_B])
    o = np.array([M[a, b] for a, b in others])
    rank = int((o >= t).sum()) + 1
    print(f"{label:>44s} {t:10.4f} {np.median(o):13.4f} {o.max():10.4f} {t - o.max():8.4f} "
          f"{str(rank) + '/' + str(len(others) + 1):>10s}")
    return t, float(np.median(o)), float(o.max()), rank


rows = {}
for label, M in (("coherence", con_src["coh"]),
                 ("wPLI", con_src["wpli"]),
                 ("coherence, symmetric orthogonalisation", con_orth["coh"]),
                 ("wPLI, symmetric orthogonalisation", con_orth["wpli"]),
                 ("envelope correlation", np.abs(ENV)),
                 ("envelope correlation, pairwise orth.", np.abs(ENV_ORTH))):
    rows[label] = summarise(M, label)
print()
print("EVERY measure ranks the true pair first here, so 'does it find the connection' does not separate them.")
print("What separates them is the MARGIN -- how far the true pair sits above the highest pair that contains no")
print("coupled source at all -- and that spans an order of magnitude:")
for label, (t, med, mx, rank) in sorted(rows.items(), key=lambda kv: kv[1][0] - kv[1][2]):
    print(f"    {label:>44s}: margin {t - mx:.4f}")
print()
print("Both orthogonalisations buy margin, and they buy a lot of it: the two uncorrected measures leave a")
print("leaked pair within 0.05-0.06 of the true one, which on real data -- where the truth is unknown and the")
print("effect is far weaker than this simulation's -- is no margin at all.  A ranking that is correct with a")
print("margin of 0.05 is a ranking that a different noise draw can reverse.")
print()
print("What the table does NOT support is 'wPLI needs no correction'.  wPLI is robust to INSTANTANEOUS mixing,")
print("and leakage is instantaneous -- but a leaked signal is a mixture of both true sources, and those two")
print("are lagged with respect to each other, so the mixture carries a lag that wPLI is built to see.  Leakage")
print("between two genuinely interacting regions manufactures apparent interactions with every region near")
print("either of them, and no lag-based measure is immune to that.")

fig, axes = plt.subplots(2, 3, figsize=(15, 8.4))
panels = [("coherence", con_src["coh"], 1.0), ("wPLI", con_src["wpli"], 0.6),
          ("envelope correlation", np.abs(ENV), 1.0),
          ("coherence, sym. orth.", con_orth["coh"], 1.0), ("wPLI, sym. orth.", con_orth["wpli"], 0.6),
          ("envelope corr., pairwise orth.", np.abs(ENV_ORTH), 1.0)]
short = [r.split()[0] for r in REGIONS]
for ax, (title, M, vmax) in zip(axes.ravel(), panels):
    L5.plot_matrix(M, short, title=title, ax=ax, vmin=0, vmax=vmax,
                   cbar_label="dimensionless (0 to 1)")
    ax.plot([REG_B], [REG_A], "o", mfc="none", mec="cyan", ms=12, mew=2)
    ax.plot([REG_A], [REG_B], "o", mfc="none", mec="cyan", ms=12, mew=2)
fig.suptitle(f"Simulated: ONE true connection (cyan rings), {REGIONS[REG_A].split()[0]}-"
             f"{REGIONS[REG_B].split()[0]}, at a quarter-cycle lag.  Everything else is leakage.",
             y=1.0, fontsize=11)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

SUBJECT, MAX_MINUTES = "sub-010002", 3.2
print(L5.disk_line("disk before the download"))
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')}")
    order = [epochs.ch_names.index(c) for c in MONTAGE if c in epochs.ch_names]
    have = [c for c in MONTAGE if c in epochs.ch_names]
    lemon_X = epochs.get_data(copy=True)[:, order, :]
finally:
    print()
    freed = L5.delete_files(paths)
    print(L5.disk_line("disk after the download was deleted"))
print()
print(f"{len(lemon_X)} eyes-closed epochs; {len(have)} of the {len(MONTAGE)} modelled channels are present "
      f"in this recording")
if len(have) != len(MONTAGE):
    print(f"  missing: {[c for c in MONTAGE if c not in have]} -- the leadfield is rebuilt on the channels "
          f"that exist")

info_l = L5.sphere_info(have, sfreq=float(epochs.info["sfreq"]))
_r = mne.io.RawArray(np.zeros((len(have), 10)), info_l.copy(), verbose=False)
_r.set_eeg_reference("average", projection=True, verbose=False)
info_l = _r.info
fwd_l, G_l = L5.leadfield(info_l, bem, rr, fixed_normals=nn)
G_l = G_l - G_l.mean(axis=0, keepdims=True)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    inv_l = make_inverse_operator(info_l, fwd_l, mne.make_ad_hoc_cov(info_l, std=NOISE_UV * 1e-6,
                                                                     verbose=False),
                                  loose=0.0, depth=None, fixed=True, verbose=False)
    prep_l = prepare_inverse_operator(inv_l, nave=1, lambda2=LAMBDA2, method="MNE", verbose=False)
    K_l, _, _, _ = _assemble_kernel(prep_l, None, "MNE", pick_ori=None, verbose=False)

cov_l = np.mean([np.cov(x) for x in lemon_X], axis=0)
W_l = region_kernel(K_l, labels, cov_l)
TC_l = np.einsum("rc,ect->ert", W_l, lemon_X)
flat_l = TC_l.transpose(1, 0, 2).reshape(N_REGIONS, -1)
TC_l_orth = L5.symmetric_orthogonalise(flat_l).reshape(N_REGIONS, TC_l.shape[0],
                                                       TC_l.shape[2]).transpose(1, 0, 2)
print(f"region time courses: {TC_l.shape}")

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


def shuffle_epochs(tc, rng):
    return np.stack([tc[rng.permutation(tc.shape[0]), c, :] for c in range(tc.shape[1])], axis=1)


def max_null(tc, method, n=N_SURROGATE):
    out = np.empty(n)
    for k in range(n):
        c_k, _ = L5.connectivity(shuffle_epochs(tc, rng_s), sfreq=SFREQ, methods=(method,),
                                 fmin=BAND[0], fmax=BAND[1], backend="numpy")
        out[k] = np.abs(L5.upper_pairs(np.abs(c_k[method]))).max()
    return out


N_SURROGATE_ENV = 50


def max_null_env(tc, orthogonalise, n=N_SURROGATE_ENV):
    out = np.empty(n)
    for k in range(n):
        M = L5.orthogonalised_envelope_correlation(shuffle_epochs(tc, rng_s), SFREQ, BAND,
                                                   orthogonalise=orthogonalise)
        out[k] = np.abs(L5.upper_pairs(np.abs(M))).max()
    return out


real = {}
c_plain, _ = L5.connectivity(TC_l, sfreq=SFREQ, methods=("coh", "wpli"), fmin=BAND[0], fmax=BAND[1])
c_orth, _ = L5.connectivity(TC_l_orth, sfreq=SFREQ, methods=("coh", "wpli"), fmin=BAND[0], fmax=BAND[1])
real["coherence"] = np.abs(c_plain["coh"])
real["wPLI"] = np.abs(c_plain["wpli"])
real["coherence, sym. orth."] = np.abs(c_orth["coh"])
real["wPLI, sym. orth."] = np.abs(c_orth["wpli"])
real["envelope correlation"] = np.abs(L5.orthogonalised_envelope_correlation(TC_l, SFREQ, BAND,
                                                                             orthogonalise=False))
real["envelope corr., pairwise orth."] = np.abs(L5.orthogonalised_envelope_correlation(TC_l, SFREQ, BAND,
                                                                                        orthogonalise=True))

null_coh = max_null(TC_l, "coh")
null_wpli = max_null(TC_l, "wpli")
null_coh_o = max_null(TC_l_orth, "coh")
null_wpli_o = max_null(TC_l_orth, "wpli")
THRESH = {"coherence": float(np.percentile(null_coh, 95)),
          "wPLI": float(np.percentile(null_wpli, 95)),
          "coherence, sym. orth.": float(np.percentile(null_coh_o, 95)),
          "wPLI, sym. orth.": float(np.percentile(null_wpli_o, 95)),
          "envelope correlation": float(np.percentile(max_null_env(TC_l, False), 95)),
          "envelope corr., pairwise orth.": float(np.percentile(max_null_env(TC_l, True), 95))}
n_pairs = N_REGIONS * (N_REGIONS - 1) // 2
print(f"ds-lemon {SUBJECT}, eyes closed, {BAND[0]:g}-{BAND[1]:g} Hz, {N_REGIONS} geometric regions, "
      f"{n_pairs} region pairs")
print(f"{N_SURROGATE} epoch-shuffled surrogates for the spectral measures and {N_SURROGATE_ENV} for the "
      f"envelope ones (they cost more);")
print("the null is the largest value over all pairs, so the multiple comparisons are paid for.")
print()
print(f"{'measure':>36s} {'median':>9s} {'max':>9s} {'95th pct of the null':>22s} {'pairs surviving':>17s}")
survive = {}
for label, M in real.items():
    v = L5.upper_pairs(M)
    if label in THRESH:
        k = int((v > THRESH[label]).sum())
        survive[label] = k
        print(f"{label:>36s} {np.median(v):9.4f} {v.max():9.4f} {THRESH[label]:22.4f} "
              f"{str(k) + ' of ' + str(n_pairs):>17s}")
    else:
        print(f"{label:>36s} {np.median(v):9.4f} {v.max():9.4f} {'not tested':>22s} {'-':>17s}")
print()
print("The six answers disagree, and there is no ground truth to adjudicate them.  That IS the")
print("result: at the source level, 'how many connections are there' is a question whose answer depends on")
print("the measure, the leakage correction and the null, and a paper that reports one number has chosen one")
print("of these rows without saying so.")

print("nb-5-6-source-connectivity -- L5.6 numbers (draft; TODO(confirm) at author review)")
print(f"Model: four-shell concentric sphere, {len(MONTAGE)} electrodes, {len(rr)} radial sources within "
      f"{MAX_DEPTH_MM:.0f} mm of the scalp,")
print(f"       minimum-norm inverse (lambda^2 = {LAMBDA2:.4f}), {N_REGIONS} GEOMETRIC regions from k-means "
      f"on position, seed {L5.SEED}.")
print("       NO parcellation, NO cortical surface and NO anatomical names -- see section 0.")
print(f"Real data: ds-lemon {SUBJECT}, eyes closed, first {MAX_MINUTES:g} min ({freed:.0f} MB, deleted); "
      f"{L5.licence_line('ds-lemon')}")
print()
print("ANSWER KEY -- leakage, from the resolution matrix alone (no data):")
print(f"    a region's estimate is a median {np.median(np.diag(LEAK)):.3f} its own power; the median "
      f"off-diagonal share is {np.median(L5.upper_pairs(LEAK)):.3f}")
print(f"    and the largest is {L5.upper_pairs(LEAK).max():.3f}.  This is the operator and the geometry, "
      f"before any recording.")
print()
print(f"ANSWER KEY -- ex-5-6 (which region pairs survive orthogonalisation), simulation with ONE true "
      f"connection")
print(f"    ({REGIONS[REG_A]} to {REGIONS[REG_B]}, quarter-cycle lag, "
      f"{np.linalg.norm(rr[SRC_A] - rr[SRC_B]) * 1000:.0f} mm apart):")
print(f"{'':6s}{'measure':>44s} {'true pair':>10s} {'median other':>13s} {'max other':>10s} "
      f"{'margin':>8s} {'rank':>10s}")
for label, (t, med, mx, rank) in rows.items():
    print(f"{'':6s}{label:>44s} {t:10.4f} {med:13.4f} {mx:10.4f} {t - mx:8.4f} "
          f"{str(rank) + '/' + str(len(others) + 1):>10s}")
print("      Every measure ranks the true pair first; the MARGIN over the highest leaked pair is what")
print("      separates them, and orthogonalisation multiplies it by about ten.")
print(f"      The {len(others)} other pairs contain no COUPLED source; anything they show is leakage or the")
print(f"      independent background ({N_BACKGROUND} sources at {BACKGROUND_SHARE:g} x the pair's sensor RMS).")
print()
print(f"ANSWER KEY -- real data, {N_REGIONS} regions, {n_pairs} pairs, surrogate null at the 95th percentile "
      f"of the max statistic:")
for label in real:
    print(f"      {label:>32s}: {survive[label]:2d} of {n_pairs} pairs survive "
          f"(threshold {THRESH[label]:.4f}, median {np.median(L5.upper_pairs(real[label])):.4f})")
print()
print("WHAT A READER LOSES HERE, STATED PLAINLY: the regions are geometric, so no surviving pair can be named,")
print("and 'which parcels are connected' is not answerable from anything on this site.  What a reader keeps is")
print("every property of leakage and of the two orthogonalisations, which are properties of the inverse")
print("operator and the geometry rather than of the labels.  If the author settles the ds-mne-sample licence")
print("(spec section 13), a cortical parcellation becomes possible and this notebook's structure carries over")
print("unchanged: only the region definition would change.")
print()
print("Pitfall: pf-volume-conduction-connectivity.")
print(L5.disk_line("disk at the end"))