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

import json

RADIUS_M, SIGMA = L5.HEAD_RADIUS_M, 0.33
MONTAGE = ["Fp1", "Fpz", "Fp2", "F7", "F3", "Fz", "F4", "F8", "T7", "C3", "Cz", "C4", "T8",
           "P7", "P3", "Pz", "P4", "P8", "O1", "Oz", "O2"]

# The widget ships its own copy of these positions.  Use them when the repository is at hand so the
# notebook and the browser are demonstrably the same geometry; fall back to rebuilding them from MNE's
# montage otherwise.  Both paths are printed.
_fix = next((d / "site/public/data/widgets/w-volume-conduction-sim/fixtures.json"
             for d in (Path.cwd(), *Path.cwd().parents)
             if (d / "site/public/data/widgets/w-volume-conduction-sim/fixtures.json").exists()), None)
if _fix is not None:
    FIX = json.loads(_fix.read_text())
    MONTAGE = list(FIX["montage"]["channels"])
    POS = np.array(FIX["montage"]["positions_xyz_m"], float)
    print(f"montage: {len(MONTAGE)} electrodes from the widget's own fixtures.json "
          f"({FIX['montage']['source']})")
else:
    FIX, POS = None, None
    print(f"montage: {len(MONTAGE)} electrodes rebuilt from MNE's standard_1005 (widget fixtures not found)")

info = L5.sphere_info(MONTAGE, POS, sfreq=250.0)
models = L5.sphere_models(RADIUS_M)
for name, bem in models.items():
    print(f"  {name:11s}: {L5.SPHERE_CALLS[name]}")
    print(f"               {L5.SPHERE_NOTES[name]}")
    print(f"               layers (relative radius, sigma S/m): {L5.sphere_conductivities(bem)}")
sphere = models["one-shell"]

names, GC = L5.great_circle_cm(info, RADIUS_M)
iu = np.triu_indices(len(names), 1)
k = int(np.argmax(GC[iu]))
FAR = (names[iu[0][k]], names[iu[1][k]])
_, CH = L5.chord_cm(info)
print()
print(f"{len(names)} sensors, {iu[0].size} pairs; the farthest pair across the scalp is "
      f"{FAR[0]}-{FAR[1]} at {GC[iu][k]:.2f} cm (arc along the scalp), {CH[iu][k]:.2f} cm (straight line "
      f"through the head)")

# Check the forward model against the widget's stored MNE-computed potentials, if they are there.
if FIX is not None:
    print("Forward check against w-volume-conduction-sim/fixtures.json (values computed by MNE 1.10.2):")
    for d in FIX["verification"]["dipoles"]:
        for key, model in (("homogeneous", "one-shell"), ("four_layer", "four-shell")):
            got = L5.dipole_potentials(info, models[model], d["pos_m"], d["moment_nAm"])
            want = np.array(d["potentials_uV"][key], float)
            rel = np.abs(got - want).max() / np.abs(want).max()
            print(f"  {d['id']:24s} {model:11s}: max |difference| = {rel:.2e} of the largest potential")
    print("  (the fixtures are rounded to six decimals, which is where the residual comes from)")
else:
    print("widget fixtures not found; the forward model is not cross-checked in this run")

SFREQ, DURATION_S, SEGMENT_S = 250.0, 20.0, 2.0
SOURCE = dict(depth=0.35, azimuth_deg=180.0, elevation_deg=22.0, strength_nAm=50.0)   # the widget's default
F_HZ = 10.0
K_SEGMENTS = int(DURATION_S / SEGMENT_S)

pos, moment = L5.dipole_from_params(**SOURCE)
gain_uv = L5.dipole_potentials(info, sphere, pos, moment)      # uV per unit of the source waveform
print(f"source: depth {SOURCE['depth']} R, azimuth {SOURCE['azimuth_deg']:.0f} deg, elevation "
      f"{SOURCE['elevation_deg']:.0f} deg, radial, {SOURCE['strength_nAm']:.0f} nA.m at {F_HZ:g} Hz")
print(f"  position {np.round(pos, 4)} m, moment {np.round(moment, 2)} nA.m")
order = np.argsort(-np.abs(gain_uv))
print("  scalp amplitude of a unit-amplitude source, largest first (uV):")
print("   " + "  ".join(f"{names[i]} {gain_uv[i]:+.2f}" for i in order[:6]))
print(f"  {FAR[0]} receives {abs(gain_uv[names.index(FAR[1])] / gain_uv[names.index(FAR[0])]):.2f}x "
      f"less of this source than {FAR[1]} does, and they are {GC[names.index(FAR[0]), names.index(FAR[1])]:.2f} cm apart")

rng = np.random.default_rng(L5.SEED)
t = np.arange(0, DURATION_S, 1 / SFREQ)
tone = np.cos(2 * np.pi * F_HZ * t)
sensors_clean = gain_uv[:, None] * tone[None, :]      # uV

fig, axes = plt.subplots(1, 2, figsize=(11, 3.4), gridspec_kw={"width_ratios": [1, 2.1]})
im, _ = mne.viz.plot_topomap(gain_uv, info, axes=axes[0], show=False, contours=6, sensors=True)
cb = fig.colorbar(im, ax=axes[0], shrink=0.85); cb.set_label("uV")
axes[0].set_title(f"one radial source, scalp map (uV per unit source)", fontsize=9)
ia, ib = names.index(FAR[0]), names.index(FAR[1])
m = t < 3.0
axes[1].plot(t[m], sensors_clean[ib][m], label=f"{FAR[1]} ({gain_uv[ib]:+.2f} uV)", lw=1.2)
axes[1].plot(t[m], sensors_clean[ia][m], label=f"{FAR[0]} ({gain_uv[ia]:+.2f} uV)", lw=1.2)
axes[1].set_xlabel("time (s)"); axes[1].set_ylabel("amplitude (uV)")
axes[1].set_title(f"the two most distant sensors, {GC[ia, ib]:.1f} cm apart -- one waveform, two scale factors",
                  fontsize=9)
axes[1].legend(fontsize=7); axes[1].grid(alpha=0.25)
fig.suptitle("A single source in a homogeneous sphere: topography (uV) and sensor traces (uV)", y=1.02, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

from scipy import signal as sp_signal


def coherence_matrix(x, sfreq=SFREQ, segment_s=SEGMENT_S, band=(F_HZ, F_HZ)):
    """Magnitude-squared coherence, Welch cross-spectra over non-overlapping Hann segments."""
    nper = int(round(segment_s * sfreq))
    f, S = sp_signal.csd(x[:, None, :], x[None, :, :], fs=sfreq, nperseg=nper, noverlap=0,
                         window="hann", detrend="constant", axis=-1)
    m = (f >= band[0] - 1e-9) & (f <= band[1] + 1e-9)
    P = np.real(np.einsum("iif->if", S))
    with np.errstate(divide="ignore", invalid="ignore"):
        g2 = (np.abs(S[..., m]) ** 2 / (P[:, None, m] * P[None, :, m])).mean(axis=-1)
    return np.nan_to_num(g2)


def epoch(x, sfreq=SFREQ, segment_s=SEGMENT_S):
    """Split a continuous record into non-overlapping epochs, for the helpers' estimator."""
    nper = int(round(segment_s * sfreq))
    n = x.shape[-1] // nper
    return x[:, : n * nper].reshape(x.shape[0], n, nper).transpose(1, 0, 2)


G2 = coherence_matrix(sensors_clean)
lib, backend = L5.connectivity(epoch(sensors_clean), sfreq=SFREQ, methods=("coh", "imcoh", "wpli"),
                               fmin=F_HZ, fmax=F_HZ)
v = L5.upper_pairs(G2)
print(f"ONE SOURCE, NO NOISE -- magnitude-squared coherence at {F_HZ:g} Hz over all {v.size} sensor pairs")
print(f"  Welch cross-spectra (scipy.signal.csd, {K_SEGMENTS} non-overlapping {SEGMENT_S:g}-s Hann segments):")
print(f"      minimum {v.min():.15f}   maximum {v.max():.15f}   largest deviation from 1: {np.abs(v - 1).max():.2e}")
print(f"  the most distant pair {FAR[0]}-{FAR[1]} ({GC[ia, ib]:.2f} cm): {G2[ia, ib]:.15f}")
c = L5.upper_pairs(np.abs(lib["coh"]))
print(f"  independent estimator ({backend['backend'].split('(')[0].strip()}), one DFT per epoch, "
      f"|coherency| not its square:")
print(f"      minimum {c.min():.15f}   maximum {c.max():.15f}")
print(f"  |imaginary coherency|, same data: maximum {np.abs(L5.upper_pairs(lib['imcoh'])).max():.2e} -- zero, "
      f"because instantaneous mixing has no lag")
print()
print("A WARNING THIS CASE PRINTS FOR FREE.  wPLI is |mean(Im S)| / mean(|Im S|), and with one noiseless "
      "source both")
print("halves are exactly zero: every Im S here is around 1e-17, which is rounding error.  The ratio is "
      "therefore 0/0 and")
print(f"  comes out at {np.abs(L5.upper_pairs(lib['wpli'])).max():.4f} -- a number with no meaning at all.  "
      "A lagged measure needs something lagged to")
print("  measure; it does not degrade gracefully towards zero when there is nothing there.  Section 4 uses "
      "it where it applies.")
print()
print("ANSWER KEY (ex-5-1): with a single source the coherence between two distant sensors is 1, exactly, "
      "with tolerance 0.")

# The result does not depend on the waveform, the frequency, the segment count or the distance.
print("The same measurement under four changes, to show which of them the result depends on:")
broadband = rng.standard_normal(t.size)
narrow = mne.filter.filter_data(rng.standard_normal(t.size), SFREQ, 8.0, 12.0, verbose=False)
cases = {
    "10 Hz tone, 10 segments": (gain_uv[:, None] * tone[None, :], SEGMENT_S, (F_HZ, F_HZ)),
    "broadband white source": (gain_uv[:, None] * broadband[None, :], SEGMENT_S, (2.0, 40.0)),
    "narrowband 8-12 Hz source": (gain_uv[:, None] * narrow[None, :], SEGMENT_S, (8.0, 12.0)),
    "tone, 4-s segments (5 of them)": (gain_uv[:, None] * tone[None, :], 4.0, (F_HZ, F_HZ)),
}
print(f"{'case':34s} {'min over pairs':>16s} {'max over pairs':>16s} {'farthest pair':>16s}")
for label, (x, seg, band) in cases.items():
    M = coherence_matrix(x, segment_s=seg, band=band)
    vv = L5.upper_pairs(M)
    print(f"  {label:32s} {vv.min():16.12f} {vv.max():16.12f} {M[ia, ib]:16.12f}")
print()
print("Nothing moves.  The cancellation is algebraic, so distance, waveform, band and segment count are all")
print("irrelevant -- which is the point: a coherence of 1 between two sensors 28 cm apart is not a finding.")

nper = int(round(SEGMENT_S * SFREQ))
w = sp_signal.get_window("hann", nper, fftbins=True)
seg_tone = tone[:nper] * w
bin_index = int(round(F_HZ * nper / SFREQ))
S_bin = np.abs(np.fft.rfft(seg_tone)[bin_index]) ** 2      # |S|^2 of a unit-amplitude source in that bin
nu_per_sigma2 = float((w ** 2).sum())                      # noise bin power per unit variance
print(f"bin arithmetic at {F_HZ:g} Hz, {SEGMENT_S:g}-s Hann segments ({nper} samples, bin {bin_index}):")
print(f"  |S|^2 of a unit source = {S_bin:.3f};  noise bin power = sigma^2 * {nu_per_sigma2:.1f}")

NOISE_UV = [0.0, 5.0, 10.0, 20.0]
rows = []
for sigma in NOISE_UV:
    noise = rng.standard_normal(sensors_clean.shape) * sigma
    M = coherence_matrix(sensors_clean + noise)
    if sigma > 0:
        snr = gain_uv ** 2 * S_bin / (sigma ** 2 * nu_per_sigma2)
        rho = snr / (1 + snr)
    else:
        rho = np.ones(len(names))          # no noise: every sensor's power is the source's, so rho = 1
    pred = rho[:, None] * rho[None, :]
    rows.append((sigma, M[ia, ib], pred[ia, ib], float(np.median(L5.upper_pairs(M))),
                 float(L5.upper_pairs(M).min())))
print()
print(f"{'sensor noise (uV SD)':>21s} {FAR[0] + '-' + FAR[1] + ' estimated':>22s} {'predicted rho_i*rho_j':>22s} "
      f"{'median over pairs':>18s} {'smallest pair':>15s}")
for sigma, est, pred_v, med, mn in rows:
    print(f"{sigma:21.0f} {est:22.4f} {pred_v:22.4f} {med:18.4f} {mn:15.4f}")
print()
print("w-volume-conduction-sim reports, for the same montage, source and segmentation, 1.0000 / 0.9340 / "
      "0.8038 / 0.5585 estimated")
print("and 1.0000 / 0.9258 / 0.7556 / 0.4277 predicted.  The predictions are deterministic and must agree "
      "exactly; the")
print("estimates are one draw of a random variable and the two implementations do not share a random number "
      "generator.")

# How far apart are the two estimates, and is the gap inside the estimator's own scatter?
WIDGET_EST = {0.0: 1.0000, 5.0: 0.9340, 10.0: 0.8038, 20.0: 0.5585}
WIDGET_PRED = {0.0: 1.0000, 5.0: 0.9258, 10.0: 0.7556, 20.0: 0.4277}
N_SEEDS = 200
print(f"One draw is not a result.  The same measurement over {N_SEEDS} noise realisations, against both "
      f"predictions:")
print()
print(f"{'sigma (uV)':>10s} {'predicted here':>15s} {'widget pred.':>13s} {'this draw':>10s} "
      f"{'widget draw':>12s} {'mean of ' + str(N_SEEDS):>13s} {'SD':>8s} {'pred + bias':>12s}")
seed_stats = {}
for sigma, est, pred_v, _, _ in rows:
    if sigma > 0:
        r2 = np.random.default_rng(L5.SEED + 1)
        draws = np.array([coherence_matrix(sensors_clean + r2.standard_normal(sensors_clean.shape) * sigma)[ia, ib]
                          for _ in range(N_SEEDS)])
        mean_draw, spread = float(draws.mean()), float(draws.std())
    else:
        mean_draw, spread = 1.0, 0.0
    seed_stats[sigma] = (mean_draw, spread)
    bias = pred_v + (1 - pred_v) / K_SEGMENTS
    print(f"{sigma:10.0f} {pred_v:15.4f} {WIDGET_PRED[sigma]:13.4f} {est:10.4f} "
          f"{WIDGET_EST[sigma]:12.4f} {mean_draw:13.4f} {spread:8.4f} {bias:12.4f}")
print()
print(f"The last column is the prediction plus (1 - gamma^2)/K, the finite-sample bias of a coherence "
      f"estimate from K = {K_SEGMENTS} segments.")
print("It tracks the mean over seeds closely, which is what says the two numbers in the middle are the same "
      "estimator seen twice")
print("rather than two different pipelines: the widget's draw and this notebook's draw both sit inside one SD "
      "of that mean.")

LONG_S = 120.0                                    # 60 non-overlapping 2-s epochs, so the estimates settle
COUPLING = 0.6                                    # share of source B's amplitude that is a delayed copy of A
t_long = np.arange(0, LONG_S, 1 / SFREQ)
POS_A = dict(depth=0.35, azimuth_deg=-90.0, elevation_deg=54.0, strength_nAm=50.0)   # left central
POS_B = dict(depth=0.35, azimuth_deg=90.0, elevation_deg=54.0, strength_nAm=50.0)    # right central
gA = L5.dipole_potentials(info, sphere, *L5.dipole_from_params(**POS_A))
gB = L5.dipole_potentials(info, sphere, *L5.dipole_from_params(**POS_B))

r3 = np.random.default_rng(L5.SEED + 7)
band = lambda x: mne.filter.filter_data(x, SFREQ, 8.0, 12.0, verbose=False)
sA = band(r3.standard_normal(t_long.size))
sB_ind = band(r3.standard_normal(t_long.size))
quarter = np.real(sp_signal.hilbert(sA) * np.exp(-1j * np.pi / 2))     # A shifted by a quarter cycle
sA, sB_ind, quarter = [z / z.std() for z in (sA, sB_ind, quarter)]
# A genuine lagged interaction is not a copy: source B is part delayed-A and part its own process, so the
# two are neither independent nor identical and the sensor mixture really does have rank 2.
sB_lag = COUPLING * quarter + np.sqrt(1 - COUPLING ** 2) * sB_ind

CASES = {
    "two independent sources": gA[:, None] * sA[None, :] + gB[:, None] * sB_ind[None, :],
    f"two sources, B = {COUPLING:g} x (A lagged 90 deg) + independent":
        gA[:, None] * sA[None, :] + gB[:, None] * sB_lag[None, :],
}
summary = {}
print(f"{LONG_S:g} s at {SFREQ:g} Hz = {int(LONG_S / SEGMENT_S)} non-overlapping {SEGMENT_S:g}-s epochs, "
      f"alpha band 8-12 Hz")
print(f"true source-level coherence between A and B:")
for label, sB in (("independent", sB_ind), ("coupled", sB_lag)):
    src_pair = np.stack([sA, sB])[None]
    c, _ = L5.connectivity(epoch(np.stack([sA, sB])), sfreq=SFREQ, methods=("coh", "imcoh", "wpli"),
                           fmin=8.0, fmax=12.0)
    print(f"   {label:12s}: coh {np.abs(c['coh'])[0, 1]:.4f}   |imcoh| {np.abs(c['imcoh'])[0, 1]:.4f}   "
          f"wPLI {np.abs(c['wpli'])[0, 1]:.4f}")
print()
for label, x in CASES.items():
    con, _ = L5.connectivity(epoch(x), sfreq=SFREQ, methods=("coh", "imcoh", "wpli"), fmin=8.0, fmax=12.0)
    summary[label] = {k: np.abs(v) for k, v in con.items()}
    print(f"AT THE SENSORS -- {label}:")
    for k in ("coh", "imcoh", "wpli"):
        vv = L5.upper_pairs(summary[label][k])
        print(f"   |{k:5s}|  median {np.median(vv):.4f}   min {vv.min():.4f}   max {vv.max():.4f}   "
              f"{FAR[0]}-{FAR[1]} {summary[label][k][ia, ib]:.4f}")
print()
print("Coherence is high in both cases -- including the one where the two generators share nothing -- so on "
      "its own it")
print("cannot distinguish an interaction from a mixture.  The imaginary measures are near zero for the "
      "independent pair and")
print("clearly non-zero for the coupled one.")

fig, axes = plt.subplots(2, 3, figsize=(12.5, 7.2))
for row, (label, mats) in enumerate(summary.items()):
    for col, key in enumerate(("coh", "imcoh", "wpli")):
        L5.plot_matrix(mats[key], names, title=f"|{key}| -- {label}", ax=axes[row, col],
                       vmin=0, vmax=1 if key == "coh" else 0.6,
                       cbar_label="dimensionless (0 to 1)")
fig.suptitle("Alpha-band (8-12 Hz) sensor connectivity, 21 electrodes, simulated: "
             "coherence cannot separate mixing from lag; the imaginary measures can", y=1.0, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

# One consequence worth stating, because it is exact rather than approximate.  With exactly TWO generators,
# Im(S_ij) = (a_i b_j - a_j b_i) * Im(S_AB): every sensor pair's imaginary cross-spectrum is the SAME quantity
# times a real constant.  wPLI divides one by the mean of the other, so the constant cancels and wPLI takes
# the same value for every pair.  That is what the min-max spread above is saying, and it is worth seeing.
for label, mats in summary.items():
    vv = L5.upper_pairs(mats["wpli"])
    print(f"wPLI spread over all {vv.size} sensor pairs, {label[:38]}: "
          f"max - min = {vv.max() - vv.min():.2e}")
print("With two sources wPLI is one number for the whole head, however many electrodes there are.  A real "
      "recording has")
print("many generators, so this degeneracy is a property of the simulation -- but it is also a reminder that "
      "a connectivity")
print("matrix with far fewer independent generators than sensors has far fewer degrees of freedom than it "
      "has entries.")


# Coherence and wPLI against distance, single source and two independent sources.
single, _ = L5.connectivity(epoch(sensors_clean + rng.standard_normal(sensors_clean.shape) * 5.0),
                            sfreq=SFREQ, methods=("coh", "wpli"), fmin=F_HZ, fmax=F_HZ)
two = summary["two independent sources"]
d = L5.upper_pairs(GC)
fig, axes = plt.subplots(1, 2, figsize=(12, 4.0))
L5.plot_versus_distance(d, {"coherence": L5.upper_pairs(np.abs(single["coh"])),
                            "wPLI": L5.upper_pairs(np.abs(single["wpli"]))},
                        title="one source + 5 uV sensor noise", ax=axes[0], ylabel="dimensionless (0 to 1)")
L5.plot_versus_distance(d, {"coherence": L5.upper_pairs(two["coh"]),
                            "wPLI": L5.upper_pairs(two["wpli"])},
                        title="two independent sources, no noise", ax=axes[1], ylabel="dimensionless (0 to 1)")
for ax in axes:
    ax.set_ylim(-0.02, 1.02)
fig.suptitle("Sensor coherence and wPLI against distance across the scalp (cm), 8-12 Hz -- "
             "simulated, so every value shown is mixing", y=1.02, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

_side = next((d / "site/public/figures/lessons/volume-conduction-intracranial-counterpoint.json"
              for d in (Path.cwd(), *Path.cwd().parents)
              if (d / "site/public/figures/lessons/volume-conduction-intracranial-counterpoint.json").exists()),
             None)
if _side is None:
    print("the counterpoint figure's sidecar was not found in this checkout; its recorded numbers are:")
    print("  ds-hup sub-HUP060, alpha coherence 0.7745 between adjacent SEEG contacts, 0.2015 three apart;")
    print("  ds-eegbci S001R02, 0.7764 for scalp pairs under 5.9 cm and 0.1007 for pairs over 17.0 cm.")
    CP = None
else:
    CP = json.loads(_side.read_text())
    n = CP["numbers"]
    print(f"Figure: site/public/figures/lessons/volume-conduction-intracranial-counterpoint.png")
    print(f"  datasets {CP['datasets']}, licence {CP['license']}, subject {CP['subject']}, "
          f"generated by {CP['generated_by']}")
    print(f"  estimator: {n['estimator']}, band {n['band_hz']} Hz")
    print(f"  intracranial ({n['n_good_contacts']} contacts in {n['n_electrode_blocks']} blocks, "
          f"{n['ieeg_sfreq_hz']:g} Hz): median coherence by contact separation {n['ieeg_median_by_separation']}")
    print(f"  scalp ({n['scalp_sfreq_hz']:g} Hz, eyes closed): median {n['scalp_median_near']:.4f} under "
          f"{n['scalp_near_cutoff_cm']:g} cm, {n['scalp_median_far']:.4f} over {n['scalp_far_cutoff_cm']:g} cm")
    print(f"  scalp distance at which coherence reaches the SEEG three-contact value: "
          f"{n['scalp_distance_reaching_ieeg_three_apart_cm']:g} cm")
    print(f"  caveats recorded with the figure: {n['caveat']}")
    print(f"  {n['contact_spacing_mm']}")

if CP is not None:
    n = CP["numbers"]
    fig, axes = plt.subplots(1, 2, figsize=(11.5, 3.8))
    sep = sorted(int(k) for k in n["ieeg_median_by_separation"])
    axes[0].plot(sep, [n["ieeg_median_by_separation"][str(s)] for s in sep], "-o", color="C1",
                 label="intracranial (contacts apart)")
    axes[0].set_xlabel("separation (contacts)"); axes[0].set_ylabel("median coherence (dimensionless)")
    axes[0].set_title("ds-hup SEEG, 8-12 Hz", fontsize=9); axes[0].set_xticks(sep); axes[0].grid(alpha=0.25)
    cm = {float(k): v for k, v in n["scalp_median_by_cm"].items() if v is not None}
    axes[1].plot(sorted(cm), [cm[k] for k in sorted(cm)], "-o", color="C0", label="scalp (cm apart)")
    axes[1].axhline(n["ieeg_median_three_apart"], ls="--", color="C1",
                    label=f"SEEG three contacts apart ({n['ieeg_median_three_apart']:.2f})")
    axes[1].set_xlabel("separation across the scalp (cm)"); axes[1].set_ylabel("median coherence (dimensionless)")
    axes[1].set_title("ds-eegbci scalp, 8-12 Hz, eyes closed", fontsize=9)
    axes[1].legend(fontsize=7); axes[1].grid(alpha=0.25)
    for ax in axes:
        ax.set_ylim(0, 1)
    fig.suptitle("'Nearby' means millimetres inside the head and something like ten centimetres on it "
                 "(recorded numbers, redrawn -- not recomputed)", y=1.03, fontsize=10)
    fig.tight_layout()
    plt.show()   # render the static figure(s) of this cell inline
else:
    print("nothing to redraw: the figure sidecar is not in this checkout")

print("nb-5-1-vc-sim -- L5.1 numbers (draft; TODO(confirm) at author review)")
print(f"Model: {L5.SPHERE_CALLS['one-shell']}")
print(f"       {L5.SPHERE_NOTES['one-shell']}; {len(names)} electrodes on a {RADIUS_M * 100:.0f} cm sphere")
print(f"Simulation: {DURATION_S:g} s at {SFREQ:g} Hz, {K_SEGMENTS} non-overlapping {SEGMENT_S:g}-s Hann "
      f"segments, seed {L5.SEED}; NO DATASET")
print(f"Source: depth {SOURCE['depth']} R, azimuth {SOURCE['azimuth_deg']:.0f} deg, elevation "
      f"{SOURCE['elevation_deg']:.0f} deg, radial, {SOURCE['strength_nAm']:.0f} nA.m, {F_HZ:g} Hz")
print()
print("ANSWER KEY -- ex-5-1 (coherence between two distant sensors with a single source):")
print(f"    1, EXACTLY.  Tolerance 0.")
print(f"    Measured here on the farthest pair {FAR[0]}-{FAR[1]} ({GC[ia, ib]:.2f} cm apart across the scalp): "
      f"{G2[ia, ib]:.15f}")
print(f"    Over all {L5.upper_pairs(G2).size} pairs the largest deviation from 1 is "
      f"{np.abs(L5.upper_pairs(G2) - 1).max():.2e}, which is floating-point arithmetic and nothing else.")
print(f"    It holds for a tone, a broadband source and a narrowband source, for 5 and for 10 segments, at "
      f"every distance (section 2).")
print()
print("ANSWER KEY -- what independent sensor noise does to that 1 (the same pair):")
print(f"{'sigma (uV)':>12s} {'estimated here':>16s} {'predicted here':>16s} {'widget estimate':>17s} "
      f"{'widget prediction':>19s}")
for sigma, est, pred_v, _, _ in rows:
    print(f"{sigma:12.0f} {est:16.4f} {pred_v:16.4f} {WIDGET_EST.get(sigma, 1.0):17.4f} "
          f"{WIDGET_PRED[sigma]:19.4f}")
print("    The predictions match the widget's to four decimals, which is the check that matters: they are")
print("    deterministic.  The estimates differ by one draw of the noise -- this notebook and the widget do")
print("    not share a random number generator -- and the difference is inside the estimator's own SD over")
print("    seeds, printed in section 3.")
print()
print(f"ANSWER KEY -- two sources, alpha band, no sensor noise, {int(LONG_S / SEGMENT_S)} epochs (section 4):")
for label, mats in summary.items():
    print(f"    {label}:")
    for k in ("coh", "imcoh", "wpli"):
        vv = L5.upper_pairs(mats[k])
        print(f"        |{k:5s}| median {np.median(vv):.4f}, min {vv.min():.4f}, max {vv.max():.4f}")
print("    With two generators wPLI is identical for every sensor pair (proved and checked in section 4); "
      "that is a property")
print("    of a two-source simulation, not of wPLI.")
print()
print("Pitfall: pf-volume-conduction-connectivity.  Widget: w-volume-conduction-sim.")
print("Counterpoint figure (not recomputed here): "
      "site/public/figures/lessons/volume-conduction-intracranial-counterpoint.png, ds-hup (CC0) + ds-eegbci.")
print()
print(L5.disk_line("disk at the end"))
print("Downloads in this notebook: none.  Sections 1-4 use no dataset; section 5 reads a figure sidecar "
      "already in the repository.")