# 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", "pandas")
_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_l3 as L3
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-erpcore", "ds-lemon", notes=True)
print()
print("The three ERP CORE statements, verbatim, as helpers_l3 records them:")
for k, v in L3.ERPCORE_LICENCE_STATEMENTS.items():
    print(f"  {k}: {v}")

CSD_PARAMS = dict(lambda2=1e-5, stiffness=4, n_legendre_terms=50)
print("mne.preprocessing.compute_current_source_density(inst, sphere='auto', "
      + ", ".join(f"{k}={v}" for k, v in CSD_PARAMS.items()) + ")")
print()

# Reference independence, first on a simulation where the map is known exactly (no dataset).
SPHERE_MONTAGE = ["Fp1", "Fpz", "Fp2", "F7", "F3", "Fz", "F4", "F8", "T7", "C3", "Cz", "C4", "T8",
                  "P7", "P3", "Pz", "P4", "P8", "O1", "Oz", "O2"]
sim_info = L5.sphere_info(SPHERE_MONTAGE, sfreq=250.0)
sim_bem = L5.sphere_models()["four-shell"]
pos, moment = L5.dipole_from_params(0.35, 180.0, 22.0, strength_nAm=50.0)
sim_map = L5.dipole_potentials(sim_info, sim_bem, pos, moment)          # uV
sim_ev = mne.EvokedArray(np.repeat(sim_map[:, None] * 1e-6, 3, axis=1), sim_info, tmin=0.0, verbose=False)


def rereferenced(ev, mode):
    out = ev.copy()
    if mode == "average":
        out.set_eeg_reference("average", projection=False, verbose=False)
    elif mode == "none (as computed)":
        pass
    else:
        out.data = out.data - out.data[out.ch_names.index(mode)]
    return out


print("CSD of ONE simulated dipole map under four references (no dataset; sphere model):")
base, ref_independence = None, 0.0
for mode in ("none (as computed)", "average", "Cz", "Oz"):
    c = mne.preprocessing.compute_current_source_density(rereferenced(sim_ev, mode), **CSD_PARAMS,
                                                         verbose=False)
    if base is None:
        base = c.data
        print(f"  {mode:20s}: reference map, peak |CSD| = {np.abs(c.data).max():.4e} (MNE units)")
    else:
        rel = float(np.abs(c.data - base).max() / np.abs(base).max())
        ref_independence = max(ref_independence, rel)
        print(f"  {mode:20s}: max |difference| from the first = {rel:.2e} of the peak")
print("  Exactly zero up to floating point, for every reference.  That is the defining property.")

ERPCORE_SUBJECTS = (1, 2, 3)          # documented subset; nb-3-4 uses ten, this notebook needs three maps
W = L3.P3_WINDOW

print(L5.disk_line("disk before the ERP CORE downloads"))

# helpers_l3.fetch_erpcore_subject checks the free space on the paradigm folder BEFORE anything creates it,
# and helpers_l3.delete_erpcore_subject (through helpers_l1.cleanup) removes that folder again as soon as it
# is empty.  A loop that fetches and deletes subject by subject therefore fails on a cold cache -- on the
# first subject if the folder never existed, and on the second if it did.  Re-creating it before each fetch
# is the whole workaround; it is reported in site/notes/notebooks-L5.md for the helpers_l3 owner.
# TODO(confirm).
evokeds, erp_info, kept = {}, None, []
for s in ERPCORE_SUBJECTS:
    (L3.erpcore_root() / "P3").mkdir(parents=True, exist_ok=True)
    try:
        ep, nfo = L3.load_p3_epochs(s, verbose=False)
        evokeds[s] = {c: L3.condition_epochs(ep, c).average().pick("eeg") for c in ("target", "standard")}
        erp_info = evokeds[s]["target"].info.copy()
        kept.append((s, nfo["n_kept"]["target"], nfo["n_kept"]["standard"]))
    finally:
        L3.delete_erpcore_subject(s, "P3", keep_small=False, verbose=False)
print(L5.disk_line("disk after the ERP CORE downloads were deleted"))
print()
print(f"ERP CORE P3, {len(evokeds)} subjects, pipeline helpers_l3.P3_PIPELINE "
      f"(the same one nb-3-1..nb-3-7 use):")
for s, nt, ns in kept:
    print(f"  sub-{s:03d}: {nt} target and {ns} standard epochs kept")

times = evokeds[ERPCORE_SUBJECTS[0]]["target"].times
erp_names = list(erp_info["ch_names"])
grand = {c: np.mean([evokeds[s][c].data for s in ERPCORE_SUBJECTS], axis=0) for c in ("target", "standard")}
grand["difference"] = grand["target"] - grand["standard"]
print(f"  {len(erp_names)} EEG channels, {len(times)} samples, window "
      f"{W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms")

def window_map(X, a, b):
    m = (times >= a) & (times <= b)
    return X[:, m].mean(axis=1)


def spread_cm(values, positions_cm):
    """Amplitude-weighted RMS distance of a map's energy from its own peak channel, in cm."""
    w = np.abs(values) ** 2
    peak = int(np.argmax(np.abs(values)))
    d = np.linalg.norm(positions_cm - positions_cm[peak], axis=1)
    return float(np.sqrt((w * d ** 2).sum() / w.sum())), peak


_, ERP_GC = L5.great_circle_cm(erp_info, L5.HEAD_RADIUS_M)
erp_pos_cm = np.array([erp_info["chs"][i]["loc"][:3] for i in range(len(erp_names))]) * 100.0

erp_ev = mne.EvokedArray(grand["difference"], erp_info, tmin=float(times[0]), verbose=False)
erp_csd = mne.preprocessing.compute_current_source_density(erp_ev, **CSD_PARAMS, verbose=False)

pot = window_map(erp_ev.data, *W) * 1e6                      # uV
csd = window_map(erp_csd.data, *W)                           # MNE units (V/m^2)
csd_disp = csd * 1e6                                         # uV/m^2 by the display convention above

rows = []
for label, v in (("potential (uV)", pot), ("CSD (uV/m^2)", csd_disp)):
    half = int((np.abs(v) >= 0.5 * np.abs(v).max()).sum())
    sp, peak = spread_cm(v, erp_pos_cm)
    rows.append((label, erp_names[peak], float(v[peak]), half, sp))
print(f"Grand-average P3 difference wave, {W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms, "
      f"{len(ERPCORE_SUBJECTS)} subjects, {len(erp_names)} channels")
print(f"{'map':>16s} {'peak channel':>13s} {'peak value':>12s} {'channels >= 50 % of peak':>25s} "
      f"{'weighted spread (cm)':>21s}")
for label, ch, val, half, sp in rows:
    print(f"{label:>16s} {ch:>13s} {val:12.3f} {half:25d} {sp:21.2f}")
r_maps = float(np.corrcoef(pot, csd_disp)[0, 1])
print()
print(f"across-channel correlation between the two maps: r = {r_maps:+.3f}")
print(f"the Laplacian narrows the map from {rows[0][4]:.2f} cm to {rows[1][4]:.2f} cm of weighted spread, "
      f"a factor of {rows[0][4] / rows[1][4]:.2f}")

fig, axes = plt.subplots(2, 5, figsize=(15, 6.2))
TIMES = [0.2, 0.3, 0.4, 0.5, 0.6]
v_pot = np.abs(np.array([window_map(erp_ev.data, t - 0.025, t + 0.025) for t in TIMES]) * 1e6).max()
v_csd = np.abs(np.array([window_map(erp_csd.data, t - 0.025, t + 0.025) for t in TIMES]) * 1e6).max()
for col, t_ in enumerate(TIMES):
    m_pot = window_map(erp_ev.data, t_ - 0.025, t_ + 0.025) * 1e6
    m_csd = window_map(erp_csd.data, t_ - 0.025, t_ + 0.025) * 1e6
    im0, _ = mne.viz.plot_topomap(m_pot, erp_info, axes=axes[0, col], show=False, contours=4,
                                  vlim=(-v_pot, v_pot), sensors=True)
    im1, _ = mne.viz.plot_topomap(m_csd, erp_info, axes=axes[1, col], show=False, contours=4,
                                  vlim=(-v_csd, v_csd), sensors=True)
    axes[0, col].set_title(f"{t_ * 1000:.0f} ms", fontsize=9)
axes[0, 0].set_ylabel("potential (uV)", fontsize=9)
axes[1, 0].set_ylabel("CSD (uV/m^2)", fontsize=9)
cb0 = fig.colorbar(im0, ax=axes[0, :], shrink=0.8); cb0.set_label("uV")
cb1 = fig.colorbar(im1, ax=axes[1, :], shrink=0.8); cb1.set_label("uV/m^2")
fig.suptitle(f"ERP CORE P3 difference wave (target - standard), {len(ERPCORE_SUBJECTS)} subjects: "
             f"potential in uV above, current source density in uV/m^2 below", y=1.0, fontsize=11)
plt.show()   # render the static figure(s) of this cell inline

def csd_of_map(values_uv, info):
    ev = mne.EvokedArray(np.repeat(np.asarray(values_uv, float)[:, None] * 1e-6, 3, axis=1), info,
                         tmin=0.0, verbose=False)
    return mne.preprocessing.compute_current_source_density(ev, **CSD_PARAMS, verbose=False).data[:, 0] * 1e6


DEPTHS = np.arange(0.20, 0.81, 0.05)
depth_rows = []
for d in DEPTHS:
    p, m = L5.dipole_from_params(float(d), 180.0, 22.0, strength_nAm=50.0)
    v = L5.dipole_potentials(sim_info, sim_bem, p, m)
    c = csd_of_map(v, sim_info)
    _, pk = spread_cm(v, np.array([sim_info["chs"][i]["loc"][:3] for i in range(len(SPHERE_MONTAGE))]) * 100)
    sp, _ = spread_cm(v, np.array([sim_info["chs"][i]["loc"][:3] for i in range(len(SPHERE_MONTAGE))]) * 100)
    depth_rows.append((float(d), d * L5.HEAD_RADIUS_M * 1000, float(np.abs(v).max()),
                       float(np.abs(c).max()), sp))

RADII = np.arange(0.0, 60.1, 7.5)
patch_rows = []
for rad in RADII:
    n = 1 if rad == 0 else 41
    total = np.zeros(len(SPHERE_MONTAGE))
    for k in range(n):
        if n == 1:
            az, el = 180.0, 22.0
        else:
            u = (k + 0.5) / n
            th = np.deg2rad(rad) * np.sqrt(u)
            ph = 2 * np.pi * ((k * 0.61803398875) % 1.0)
            az = 180.0 + np.rad2deg(th * np.cos(ph)) / max(np.cos(np.deg2rad(22.0)), 1e-6)
            el = 22.0 + np.rad2deg(th * np.sin(ph))
        p, m = L5.dipole_from_params(0.2, float(az), float(el), strength_nAm=50.0 / n)
        total += L5.dipole_potentials(sim_info, sim_bem, p, m)
    c = csd_of_map(total, sim_info)
    sp, _ = spread_cm(total, np.array([sim_info["chs"][i]["loc"][:3] for i in range(len(SPHERE_MONTAGE))]) * 100)
    patch_rows.append((float(rad), n, float(np.abs(total).max()), float(np.abs(c).max()), sp))

print("A single radial dipole moved deeper (four-shell sphere, 21 electrodes, no noise):")
print(f"{'depth (fraction of R)':>22s} {'depth below scalp (mm)':>23s} {'peak potential (uV)':>20s} "
      f"{'peak CSD (uV/m^2)':>18s} {'CSD/potential, relative':>24s}")
ref_ratio = depth_rows[0][3] / depth_rows[0][2]
for d, mm, vp, vc, sp in depth_rows:
    print(f"{d:22.2f} {mm:23.1f} {vp:20.3f} {vc:18.3f} {(vc / vp) / ref_ratio:24.3f}")
print()
print("A patch of superficial dipoles of growing angular radius, total moment held constant:")
print(f"{'patch radius (deg)':>19s} {'n dipoles':>10s} {'peak potential (uV)':>20s} {'peak CSD (uV/m^2)':>18s} "
      f"{'CSD/potential, relative':>24s} {'map spread (cm)':>16s}")
ref_ratio_p = patch_rows[0][3] / patch_rows[0][2]
for rad, n, vp, vc, sp in patch_rows:
    print(f"{rad:19.1f} {n:10d} {vp:20.3f} {vc:18.3f} {(vc / vp) / ref_ratio_p:24.3f} {sp:16.2f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4.0))
d_mm = [r[1] for r in depth_rows]
axes[0].plot(d_mm, [(r[3] / r[2]) / ref_ratio for r in depth_rows], "-o", color="C0")
axes[0].set_xlabel("depth below the scalp (mm)")
axes[0].set_ylabel("CSD / potential, relative to the shallowest (dimensionless)")
axes[0].set_title("a single radial dipole, moved deeper", fontsize=9)
axes[0].grid(alpha=0.25)
ax2 = axes[0].twinx()
ax2.plot(d_mm, [r[2] for r in depth_rows], "--", color="C1", lw=1)
ax2.set_ylabel("peak potential (uV)", color="C1")
axes[1].plot([r[0] for r in patch_rows], [(r[3] / r[2]) / ref_ratio_p for r in patch_rows], "-o", color="C0")
axes[1].set_xlabel("patch radius (degrees of arc)")
axes[1].set_ylabel("CSD / potential, relative to a point source (dimensionless)")
axes[1].set_title("a superficial patch, spread wider (total moment fixed)", fontsize=9)
axes[1].grid(alpha=0.25)
fig.suptitle("What the surface Laplacian discounts: depth (left) and spatial extent (right), "
             "four-shell sphere, 21 electrodes", y=1.02, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

print("Both curves fall, and they fall for the same reason: a deep source and a broad patch both produce a")
print("SPATIALLY BROAD scalp map, and the Laplacian is a measure of spatial curvature.  The patch curve is not")
print("monotone past about 45 degrees: a patch that large wraps far enough around the sphere to develop new")
print("structure of its own, so the map stops simply getting broader.  That is a property of a sphere with a")
print(f"{L5.HEAD_RADIUS_M * 100:.0f} cm radius and 21 electrodes, and it is reported rather than trimmed off.")
print("CSD therefore cannot")
print("distinguish 'deep' from 'large' -- it discounts both, and it discounts a distant superficial generator")
print("for the same reason.  That is the whole trade: the sharpening you want and the attenuation you may not.")

SUBJECT, MAX_MINUTES, BAND = "sub-010002", 3.2, L5.ALPHA_BAND
print(L5.disk_line("disk before the ds-lemon 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')}")
    lemon_info = epochs.info.copy()
    lemon_names = list(epochs.ch_names)
    lemon_sfreq = float(epochs.info["sfreq"])
    pot_data = epochs.get_data(copy=True) * 1e6
    csd_epochs = mne.preprocessing.compute_current_source_density(epochs.copy(), **CSD_PARAMS, verbose=False)
    csd_data = csd_epochs.get_data(copy=True) * 1e6
finally:
    print()
    freed_lemon = L5.delete_files(paths)
    print(L5.disk_line("disk after the ds-lemon download was deleted"))

print()
print(f"{len(pot_data)} eyes-closed epochs x {len(lemon_names)} channels at {lemon_sfreq:g} Hz")
METHODS = ("coh", "imcoh", "wpli")
con_pot, backend = L5.connectivity(pot_data, sfreq=lemon_sfreq, methods=METHODS, fmin=BAND[0], fmax=BAND[1])
con_csd, _ = L5.connectivity(csd_data, sfreq=lemon_sfreq, methods=METHODS, fmin=BAND[0], fmax=BAND[1])
con_pot = {k: np.abs(v) for k, v in con_pot.items()}
con_csd = {k: np.abs(v) for k, v in con_csd.items()}
print(f"backend: {backend['backend']}")
print()
_, LEM_GC = L5.great_circle_cm(lemon_info, L5.HEAD_RADIUS_M)
dist = L5.upper_pairs(LEM_GC)
print(f"{'measure':>8s} {'median before':>14s} {'median after':>13s} {'r(before, after)':>17s} "
      f"{'median < 5 cm':>14s} {'median > 15 cm':>15s}")
for m in METHODS:
    b, a_ = L5.upper_pairs(con_pot[m]), L5.upper_pairs(con_csd[m])
    print(f"{m:>8s} {np.median(b):14.4f} {np.median(a_):13.4f} "
          f"{float(np.corrcoef(b, a_)[0, 1]):17.3f} "
          f"{np.median(b[dist < 5]):.3f} -> {np.median(a_[dist < 5]):.3f}   "
          f"{np.median(b[dist > 15]):.3f} -> {np.median(a_[dist > 15]):.3f}")

fig, axes = plt.subplots(2, 3, figsize=(14, 8))
for row, (label, mats) in enumerate((("potential", con_pot), ("CSD", con_csd))):
    for col, m in enumerate(METHODS):
        L5.plot_matrix(mats[m], None, title=f"|{m}| -- {label}", ax=axes[row, col], vmin=0,
                       vmax=1.0 if m == "coh" else 0.5, cbar_label="dimensionless (0 to 1)")
fig.suptitle(f"ds-lemon {SUBJECT} eyes closed, {BAND[0]:g}-{BAND[1]:g} Hz: connectivity before (top) and "
             f"after (bottom) the surface Laplacian", y=1.0, fontsize=11)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

fig, axes = plt.subplots(1, 2, figsize=(12, 4.0))
L5.plot_versus_distance(dist, {"coherence": L5.upper_pairs(con_pot["coh"]),
                               "wPLI": L5.upper_pairs(con_pot["wpli"])},
                        title="potential", ax=axes[0], ylabel="dimensionless (0 to 1)")
L5.plot_versus_distance(dist, {"coherence": L5.upper_pairs(con_csd["coh"]),
                               "wPLI": L5.upper_pairs(con_csd["wpli"])},
                        title="after CSD", ax=axes[1], ylabel="dimensionless (0 to 1)")
for ax in axes:
    ax.set_ylim(0, 1)
fig.suptitle("Connectivity against distance across the scalp (cm), before and after the surface Laplacian, "
             "8-12 Hz", y=1.02, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

COMPONENTS = {
    "P3 difference, 300-600 ms": window_map(grand["difference"], *W),
    "early response, 100-200 ms, target": window_map(grand["target"], 0.10, 0.20),
    "early response, 100-200 ms, standard": window_map(grand["standard"], 0.10, 0.20),
    "late slow wave, 600-800 ms, difference": window_map(grand["difference"], 0.60,
                                                         min(0.80, float(times[-1]))),
    "pre-stimulus window, -200-0 ms, difference": window_map(grand["difference"], -0.20, 0.0),
}
# The descriptions are MEASURED, not asserted: "broad" and "focal" are the spread column, and the survival
# of amplitude is an RMS ratio (the peak channel of the CSD map need not be the peak channel of the
# potential map, which makes a peak-to-peak ratio the wrong statistic).
comp_rows = []
base_ratio = None
for label, v_uv in COMPONENTS.items():
    v = np.asarray(v_uv, float) * 1e6
    c = csd_of_map(v, erp_info)
    sp_b, _ = spread_cm(v, erp_pos_cm)
    sp_a, _ = spread_cm(c, erp_pos_cm)
    rms_v, rms_c = float(np.sqrt((v ** 2).mean())), float(np.sqrt((c ** 2).mean()))
    ratio = rms_c / rms_v
    if base_ratio is None:
        base_ratio = ratio
    comp_rows.append((label, float(np.abs(v).max()), rms_v, rms_c, ratio, ratio / base_ratio, sp_b, sp_a))
print(f"{'component':>44s} {'peak (uV)':>10s} {'RMS (uV)':>9s} {'relative RMS ratio':>19s} "
      f"{'spread (cm)':>12s} {'-> after':>9s} {'narrowing':>10s}")
for label, pk, rms_v, rms_c, ratio, rel, sp_b, sp_a in comp_rows:
    print(f"{label:>44s} {pk:10.3f} {rms_v:9.3f} {rel:19.2f} {sp_b:12.2f} {sp_a:9.2f} "
          f"{sp_b / sp_a:9.2f}x")
print()
print("The 'relative RMS ratio' column is each component's CSD/potential RMS ratio divided by the P3's, so the")
print("P3 is 1 by construction: the absolute ratio carries the units convention and the spline parameters and")
print("means nothing on its own.")
print()
print("TWO RESULTS, AND THE SECOND IS A NEGATIVE ONE.")
print("  (a) Every map narrows.  The narrowing factor separates the components clearly: the P3 narrows most")
print(f"      ({comp_rows[0][6] / comp_rows[0][7]:.2f}x) and the early responses least "
      f"({comp_rows[1][6] / comp_rows[1][7]:.2f}x).")
print("  (b) The AMPLITUDE that survives does NOT separate them.  All five relative ratios sit inside")
print(f"      {min(r[5] for r in comp_rows):.2f}-{max(r[5] for r in comp_rows):.2f}, a spread of well under a")
print("      factor of two, where the controlled sphere sweep above moved by a factor of two on depth alone.")
print("      Spread does not predict survival here either: the P3 has the smallest spread of the five and a")
print("      middling ratio.  On 30 electrodes, over these windows, the attenuation side of the exercise is")
print("      not decided by the data.  TODO(confirm): the exercise should either ask about narrowing, which")
print("      these data answer, or about attenuation, which the sphere answers -- not about both at once.")
print()
print(f"The pre-stimulus row is a scale check, not a component: its amplitude is {comp_rows[-1][1]:.3f} uV, two")
print("orders below the P3, so its ratio measures the noise floor rather than a field.")

print("nb-5-3-csd -- L5.3 numbers (draft; TODO(confirm) at author review)")
print(f"Transform: mne.preprocessing.compute_current_source_density(sphere='auto', "
      + ", ".join(f"{k}={v}" for k, v in CSD_PARAMS.items()) + "), MNE " + mne.__version__)
print(f"ERP data: ds-erpcore P3, sub-001..sub-{ERPCORE_SUBJECTS[-1]:03d} ({len(ERPCORE_SUBJECTS)} subjects), "
      f"helpers_l3.P3_PIPELINE; {L5.licence_line('ds-erpcore')}")
print(f"Resting data: ds-lemon {SUBJECT}, eyes closed, first {MAX_MINUTES:g} min "
      f"({freed_lemon:.0f} MB, deleted); {L5.licence_line('ds-lemon')}")
print("Sphere sweep: no dataset (MNE four-shell concentric sphere, 21 electrodes)")
print()
print("ANSWER KEY -- reference independence: CSD computed after four different references differs by at most")
print(f"    {ref_independence:.1e} of the peak -- the Laplacian of a constant is zero, so this is exact, "
      "not approximate.")
print()
print(f"ANSWER KEY -- sharpening of the P3 ({W[0] * 1000:.0f}-{W[1] * 1000:.0f} ms grand-average difference "
      f"wave):")
for label, ch, val, half, sp in rows:
    print(f"    {label:>16s}: peak {ch} = {val:+.3f}; {half} of {len(erp_names)} channels above half the "
          f"peak; weighted spread {sp:.2f} cm")
print(f"    across-channel correlation between the two maps r = {r_maps:+.3f}; spread falls by a factor of "
      f"{rows[0][4] / rows[1][4]:.2f}")
print()
print("ANSWER KEY -- what CSD attenuates (four-shell sphere, relative CSD/potential ratio, no dataset):")
print(f"    depth   {depth_rows[0][1]:.0f} mm below the scalp -> 1.000 (reference);  "
      f"{depth_rows[len(depth_rows) // 2][1]:.0f} mm -> "
      f"{(depth_rows[len(depth_rows) // 2][3] / depth_rows[len(depth_rows) // 2][2]) / ref_ratio:.3f};  "
      f"{depth_rows[-1][1]:.0f} mm -> {(depth_rows[-1][3] / depth_rows[-1][2]) / ref_ratio:.3f}")
print(f"    extent  point source -> 1.000 (reference);  "
      f"{patch_rows[len(patch_rows) // 2][0]:.0f} deg patch -> "
      f"{(patch_rows[len(patch_rows) // 2][3] / patch_rows[len(patch_rows) // 2][2]) / ref_ratio_p:.3f};  "
      f"{patch_rows[-1][0]:.0f} deg patch -> {(patch_rows[-1][3] / patch_rows[-1][2]) / ref_ratio_p:.3f}")
print("    CSD cannot separate 'deep' from 'large': both make the scalp map broad and both are discounted.")
print()
print("ANSWER KEY -- ex-5-3 (which components sharpen and which shrink):")
print("    SHARPENING, which these data do answer -- weighted spread before -> after, and the narrowing factor:")
for label, pk, rms_v, rms_c, ratio, rel, sp_b, sp_a in comp_rows:
    print(f"      {label:>44s}: {sp_b:5.2f} -> {sp_a:5.2f} cm ({sp_b / sp_a:.2f}x), "
          f"relative RMS ratio {rel:.2f}")
print("    ATTENUATION, which these data do NOT answer: all five relative RMS ratios lie inside "
      f"{min(r[5] for r in comp_rows):.2f}-{max(r[5] for r in comp_rows):.2f},")
print("      which is not a separation.  The controlled answer is the sphere sweep above: depth costs a "
      "factor of 2.1")
print("      from 18 to 72 mm, and extent a factor of 1.4 from a point to a 60-degree patch.  TODO(confirm): "
      "the")
print("      exercise's option list should be built on one of those two measurements, not on both.")
print()
print(f"ANSWER KEY -- connectivity before and after CSD (ds-lemon {SUBJECT}, {BAND[0]:g}-{BAND[1]:g} Hz, "
      f"{len(pot_data)} epochs):")
for m in METHODS:
    b, a_ = L5.upper_pairs(con_pot[m]), L5.upper_pairs(con_csd[m])
    print(f"    {m:>6s}: median {np.median(b):.4f} -> {np.median(a_):.4f}; "
          f"under 5 cm {np.median(b[dist < 5]):.3f} -> {np.median(a_[dist < 5]):.3f}; "
          f"over 15 cm {np.median(b[dist > 15]):.3f} -> {np.median(a_[dist > 15]):.3f}; "
          f"r(before, after) = {float(np.corrcoef(b, a_)[0, 1]):+.3f}")
print()
print("CSD is not an inverse solution and does not remove mixing.  It changes the spatial scale the channels")
print("are sensitive to, which reduces the contribution of distant generators and also removes any genuinely")
print("broad or deep source.  A connectivity result that appears only after CSD, or only before it, is a")
print("result about spatial scale.")
print()
print("Widget: w-reference-explorer (mode csd).  Pitfalls: none assigned to L5.3 in the spec.")
print(L5.disk_line("disk at the end"))