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

import warnings

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"]
GRID_SPACING_M, GRID_MAX_FRACTION = 0.012, 0.85
SFREQ, NOISE_UV, SNR = 250.0, 1.0, 3.0
LAMBDA2 = 1.0 / SNR ** 2

info = L5.sphere_info(MONTAGE, sfreq=SFREQ)
# An average-reference projector, because an inverse operator needs to know what the reference did.
_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=GRID_SPACING_M, max_fraction=GRID_MAX_FRACTION, min_radius_m=0.012)
nn = rr / np.linalg.norm(rr, axis=1, keepdims=True)
fwd, G = L5.leadfield(info, bem, rr, fixed_normals=nn)
depth_mm = (L5.HEAD_RADIUS_M - np.linalg.norm(rr, axis=1)) * 1000
print(f"{len(MONTAGE)} electrodes, {len(rr)} radial sources on a {GRID_SPACING_M * 1000:.0f} mm grid; "
      f"leadfield {G.shape}")

rank = np.linalg.matrix_rank(G)
print(f"rank(G) = {rank}; null space dimension = {G.shape[1] - rank} of {G.shape[1]} source dimensions")
U, S, Vt = np.linalg.svd(G, full_matrices=True)
null_vec = Vt[rank]                                       # the first direction G annihilates
sensor = G @ null_vec
NULL_REL = float(np.abs(sensor).max() / np.abs(G).max())
print(f"a unit-norm null-space source configuration produces at most {NULL_REL:.2e} of the leadfield's "
      f"largest entry at any electrode -- zero, to double precision")
print(f"  (MNE stores a forward solution in SINGLE precision; helpers_l5.leadfield hands back a float64 copy. "
      f"Left in float32")
print(f"   the same residual reads about 6e-08, which is the arithmetic and not the physics.)")
print(f"singular values of G: largest {S[0]:.3e}, smallest {S[rank - 1]:.3e}, "
      f"ratio {S[0] / S[rank - 1]:.1f}")
print()
print("The ratio above is the honest version of 'EEG cannot see deep sources': the worst-seen source "
      "direction")
print(f"reaches the sensors {S[0] / S[rank - 1]:.0f} times more weakly than the best-seen one, so noise that "
      "is small compared")
print("with the best direction is large compared with the worst.  Regularization is the decision about where "
      "to stop.")

fig, axes = plt.subplots(1, 3, figsize=(15, 4.0))
sl = np.abs(rr[:, 0]) < GRID_SPACING_M / 2
v = np.abs(null_vec[sl]).max()
sc = axes[0].scatter(rr[sl, 1] * 1000, rr[sl, 2] * 1000, c=null_vec[sl], s=70, cmap="RdBu_r",
                     vmin=-v, vmax=v)
cb = fig.colorbar(sc, ax=axes[0]); cb.set_label("source amplitude (arbitrary, unit norm)")
axes[0].set_xlabel("y, anterior + (mm)"); axes[0].set_ylabel("z, superior + (mm)")
axes[0].set_title("a source pattern the sensors cannot see\n(mid-sagittal slice)", fontsize=9)
axes[0].set_aspect("equal")

im, _ = mne.viz.plot_topomap(sensor * 1e6, info, axes=axes[1], show=False, contours=6, sensors=True)
cb = fig.colorbar(im, ax=axes[1]); cb.set_label("uV")
axes[1].set_title(f"what it produces at the scalp\n(peak {np.abs(sensor).max() * 1e6:.1e} uV)", fontsize=9)

axes[2].semilogy(np.arange(1, len(S) + 1), S, "o-", ms=4)
axes[2].set_xlabel("singular value index")
axes[2].set_ylabel("singular value of G (V per A.m)")
axes[2].set_title("how unevenly the sensors see\nthe source space", fontsize=9)
axes[2].grid(alpha=0.25)
fig.suptitle("Ill-posedness, drawn: the null space, its (absent) scalp map, and the spectrum of the leadfield",
             y=1.02, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

from mne.beamformer import apply_lcmv, make_lcmv
from mne.minimum_norm import apply_inverse, make_inverse_operator, make_inverse_resolution_matrix

rng = np.random.default_rng(L5.SEED)
N_EPOCHS, N_TIMES, F_HZ, MOMENT_nAm = 40, 250, 10.0, 50.0
wave = np.sin(2 * np.pi * F_HZ * np.arange(N_TIMES) / SFREQ)


def simulate(j, moment_nAm=MOMENT_nAm, noise_uv=NOISE_UV, n_epochs=N_EPOCHS, rng=rng):
    """Epochs of one radial dipole at grid point ``j`` plus independent sensor noise."""
    clean = G[:, [j]] @ (moment_nAm * 1e-9 * wave)[None, :]
    X = np.stack([clean + rng.standard_normal(clean.shape) * noise_uv * 1e-6 for _ in range(n_epochs)])
    return mne.EpochsArray(X, info, tmin=0.0, verbose=False), clean


TRUE = int(np.argmin(np.abs(depth_mm - 25) + np.abs(rr[:, 0]) * 1000 + np.abs(rr[:, 1] + 0.05) * 1000))
epochs, clean = simulate(TRUE)
evoked = epochs.average()
print(f"simulated source: grid point {TRUE}, {depth_mm[TRUE]:.0f} mm below the scalp, "
      f"{MOMENT_nAm:.0f} nA.m radial at {F_HZ:g} Hz")
print(f"  peak scalp amplitude {np.abs(clean).max() * 1e6:.2f} uV; sensor noise {NOISE_UV:g} uV SD per "
      f"electrode; {N_EPOCHS} epochs averaged")
print(f"  regularization: lambda^2 = 1/SNR^2 = {LAMBDA2:.4f} (SNR {SNR:g})")

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)
    data_cov = mne.compute_covariance(epochs, verbose=False)
    lcmv = make_lcmv(info, fwd, data_cov, reg=0.05, noise_cov=noise_cov, verbose=False)

estimates = {}
for method in ("MNE", "dSPM", "sLORETA", "eLORETA"):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        stc = apply_inverse(evoked, inv, lambda2=LAMBDA2, method=method, verbose=False)
    estimates[method] = np.abs(stc.data).max(axis=1)
with warnings.catch_warnings():
    warnings.simplefilter("ignore")
    estimates["LCMV"] = np.abs(apply_lcmv(evoked, lcmv, verbose=False).data).max(axis=1)


def where(amp):
    peak = int(np.argmax(amp))
    w = amp ** 2
    d = np.linalg.norm(rr - rr[TRUE], axis=1) * 1000
    return peak, float(np.linalg.norm(rr[peak] - rr[TRUE]) * 1000), \
        float(np.sqrt((w * d ** 2).sum() / w.sum())), float(depth_mm[peak])


print()
print(f"{'method':>10s} {'peak error (mm)':>16s} {'spatial dispersion (mm)':>24s} {'peak depth (mm)':>16s} "
      f"{'true depth (mm)':>16s}")
est_rows = []
for method, amp in estimates.items():
    peak, err, disp, pd = where(amp)
    est_rows.append((method, err, disp, pd))
    print(f"{method:>10s} {err:16.1f} {disp:24.1f} {pd:16.1f} {depth_mm[TRUE]:16.1f}")
print()
print("Peak error alone is a poor summary: a method can put its maximum in the right place and still spread "
      "the")
print("estimate over half the head, which is what the dispersion column is for.  Both are reported "
      "everywhere below.")

fig, axes = plt.subplots(1, len(estimates), figsize=(3.1 * len(estimates), 3.6))
sl = np.abs(rr[:, 0] - rr[TRUE, 0]) < GRID_SPACING_M / 2
for ax, (method, amp) in zip(axes, estimates.items()):
    a = amp / amp.max()
    sc = ax.scatter(rr[sl, 1] * 1000, rr[sl, 2] * 1000, c=a[sl], s=60, cmap="magma", vmin=0, vmax=1)
    ax.plot(rr[TRUE, 1] * 1000, rr[TRUE, 2] * 1000, "o", mfc="none", mec="cyan", ms=14, mew=2)
    ax.set_title(f"{method}", fontsize=9); ax.set_aspect("equal")
    ax.set_xlabel("y, anterior + (mm)")
axes[0].set_ylabel("z, superior + (mm)")
cb = fig.colorbar(sc, ax=axes, shrink=0.8); cb.set_label("estimate, normalised to its own peak")
fig.suptitle(f"One simulated radial dipole {depth_mm[TRUE]:.0f} mm below the scalp (cyan ring), five "
             f"estimators, coronal-plane slice of the grid", y=1.02, fontsize=10)
plt.show()   # render the static figure(s) of this cell inline

res_rows = []
res_metrics = {}
for method in ("MNE", "dSPM", "sLORETA", "eLORETA"):
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        R = make_inverse_resolution_matrix(fwd, inv, method=method, lambda2=LAMBDA2)
    met = L5.resolution_metrics_grid(R, rr)
    res_metrics[method] = met
    res_rows.append((method, float(np.median(met["peak_error_mm"])), float(met["peak_error_mm"].mean()),
                     float(np.median(met["spatial_dispersion_mm"]))))

# The beamformer has a kernel too, one row per source; build its resolution matrix the same way.
K_lcmv = lcmv["weights"]
R_lcmv = L5.resolution_matrix(G, K_lcmv)
met = L5.resolution_metrics_grid(R_lcmv, rr)
res_metrics["LCMV"] = met
res_rows.append(("LCMV", float(np.median(met["peak_error_mm"])), float(met["peak_error_mm"].mean()),
                 float(np.median(met["spatial_dispersion_mm"]))))

print(f"Resolution over all {len(rr)} grid sources, noiseless (the resolution matrix has no noise in it):")
print(f"{'method':>10s} {'peak error median (mm)':>24s} {'mean':>7s} {'dispersion median (mm)':>24s}")
for method, med, mean_e, disp in res_rows:
    print(f"{method:>10s} {med:24.1f} {mean_e:7.1f} {disp:24.1f}")
print()
print("sLORETA and eLORETA reach ZERO median peak error, which is not an accident and not an achievement "
      "either:")
print("both are constructed so that the resolution matrix's diagonal dominates its column, which is exactly")
print("'the peak lands in the right place'.  Their dispersion columns say what that construction costs -- "
      "the")
print("estimate is still spread over centimetres.  Zero localisation error for a point source is a property "
      "of")
print("the estimator, not evidence that the estimate is sharp.")

fig, axes = plt.subplots(1, 2, figsize=(12.5, 4.2))
for method in res_metrics:
    o = np.argsort(depth_mm)
    for ax, key, lab in ((axes[0], "peak_error_mm", "peak localisation error (mm)"),
                         (axes[1], "spatial_dispersion_mm", "spatial dispersion (mm)")):
        y = res_metrics[method][key][o]
        # median in 10 mm depth bins, so five curves stay readable
        edges = np.arange(depth_mm.min(), depth_mm.max() + 10, 10)
        xs = 0.5 * (edges[:-1] + edges[1:])
        ys = [np.median(res_metrics[method][key][(depth_mm >= a) & (depth_mm < b)])
              if ((depth_mm >= a) & (depth_mm < b)).any() else np.nan
              for a, b in zip(edges[:-1], edges[1:])]
        ax.plot(xs, ys, "-o", ms=4, label=method)
        ax.set_xlabel("depth below the scalp (mm)"); ax.set_ylabel(lab)
        ax.grid(alpha=0.25)
axes[0].set_title("where the peak lands", fontsize=9)
axes[1].set_title("how far the estimate spreads", fontsize=9)
axes[0].legend(fontsize=8)
fig.suptitle(f"Resolution against depth, {len(MONTAGE)} electrodes, four-shell sphere "
             f"(median per 10 mm depth bin)", y=1.02, fontsize=10)
fig.tight_layout()
plt.show()   # render the static figure(s) of this cell inline

# Both the leadfield and the data are average-referenced here, because that is the reference the simulated
# measurement carries; an inverse operator built on a differently referenced leadfield is solving a different
# problem.
G_avg = G - G.mean(axis=0, keepdims=True)
col_norm = np.linalg.norm(G_avg, axis=0)
print(f"column norms of the average-referenced leadfield span {col_norm.min():.1f} to {col_norm.max():.1f} "
      f"(a ratio of {col_norm.max() / col_norm.min():.1f}),")
print("which is all the leverage a depth weighting has.")


def mne_kernel(gamma, lambda2=LAMBDA2):
    """Depth-weighted minimum-norm kernel, written out: K = R G' (G R G' + lambda^2 tr/n I)^-1."""
    w = col_norm ** (-2.0 * gamma)                          # the source covariance R = diag(w)
    GR = G_avg * w
    A = GR @ G_avg.T
    A = A + lambda2 * np.trace(A) / len(MONTAGE) * np.eye(len(MONTAGE))
    return GR.T @ np.linalg.inv(A)


def avg_ref(x):
    return x - x.mean(axis=0, keepdims=True)


DEEP = int(np.argmin(np.abs(depth_mm - 60) + np.abs(rr[:, 0]) * 1000))
deep_epochs, deep_clean = simulate(DEEP)
deep_v = deep_epochs.average().data
print(f"a second, deeper source for this section: grid point {DEEP}, {depth_mm[DEEP]:.0f} mm below the "
      f"scalp, peak {np.abs(deep_clean).max() * 1e6:.2f} uV")
print()
print(f"{'gamma':>7s} {'peak depth, shallow source (mm)':>32s} {'peak depth, deep source (mm)':>30s} "
      f"{'median peak error over the grid (mm)':>38s}")
gamma_rows = []
for gamma in (0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2.0):
    K = mne_kernel(gamma)
    a_shallow = np.abs(K @ avg_ref(evoked.data)).max(axis=1)
    a_deep = np.abs(K @ avg_ref(deep_v)).max(axis=1)
    R = L5.resolution_matrix(G_avg, K)
    met = L5.resolution_metrics_grid(R, rr)
    row = (gamma, float(depth_mm[int(np.argmax(a_shallow))]), float(depth_mm[int(np.argmax(a_deep))]),
           float(np.median(met["peak_error_mm"])), float(np.median(met["spatial_dispersion_mm"])))
    gamma_rows.append(row)
    print(f"{gamma:7.1f} {row[1]:32.1f} {row[2]:30.1f} {row[3]:38.1f}")
print()
print(f"true depths: shallow source {depth_mm[TRUE]:.0f} mm, deep source {depth_mm[DEEP]:.0f} mm")
print()
print("Depth weighting moves the estimate of a deep source deeper -- and then STOPS.  On this geometry it")
print(f"takes the {depth_mm[DEEP]:.0f} mm source's estimate from {gamma_rows[0][2]:.0f} mm to about "
      f"{max(r[2] for r in gamma_rows):.0f} mm and no further, however")
print("large gamma is made, because the column norms of the leadfield span only a factor of "
      f"{col_norm.max() / col_norm.min():.0f} across the")
print("whole grid: that is all the leverage the weighting has.  A source 60 mm below the scalp is never")
print("recovered at 60 mm by any gamma tried here.  Depth weighting is also not free -- the dispersion column")
print("grows with gamma -- and it is not a measurement: gamma is a knob, and a source estimate is only as deep")
print("as the knob was set.")
print(f"{'gamma':>7s} {'median dispersion over the grid (mm)':>38s}")
for gamma, ds, dd, pe, disp in gamma_rows:
    print(f"{gamma:7.1f} {disp:38.1f}")

def patch_dirs(centre_dir, radius_deg, n=61):
    """``n`` directions spread evenly over a cap of angular radius ``radius_deg`` about ``centre_dir``."""
    if radius_deg == 0:
        return centre_dir[None, :]
    e1 = np.cross(centre_dir, [0, 0, 1.0])
    e1 = e1 / np.linalg.norm(e1) if np.linalg.norm(e1) > 1e-9 else np.array([1.0, 0.0, 0.0])
    e2 = np.cross(centre_dir, e1)
    k = np.arange(n)
    th = np.deg2rad(radius_deg) * np.sqrt((k + 0.5) / n)
    ph = 2 * np.pi * ((k * 0.61803398875) % 1.0)
    d = (centre_dir[None, :] * np.cos(th)[:, None]
         + (e1[None, :] * np.cos(ph)[:, None] + e2[None, :] * np.sin(ph)[:, None]) * np.sin(th)[:, None])
    return d / np.linalg.norm(d, axis=1, keepdims=True)


def patch_map(centre_dir, radius_deg, depth_fraction, bem, info, n=61, moment_nAm=20.0):
    """Scalp map (uV) of a cap of radial dipoles, total moment held at ``moment_nAm``.

    One forward solution for the whole cap rather than one per dipole -- the search below tries several
    hundred caps and the difference is a minute of wall clock.
    """
    d = patch_dirs(centre_dir, radius_deg, n)
    pos = d * (1 - depth_fraction) * L5.HEAD_RADIUS_M
    _, Gp = L5.leadfield(info, bem, pos, fixed_normals=d)
    return Gp.sum(axis=1) * (moment_nAm / len(d)) * 1e-9 * 1e6


PATCH_DEPTH = 0.15
RADII = np.arange(0.0, 75.1, 2.5)
CENTRE = np.array([0.0, -np.cos(np.deg2rad(22.0)), np.sin(np.deg2rad(22.0))])
info21 = L5.sphere_info(["Fp1", "Fpz", "Fp2", "F7", "F3", "Fz", "F4", "F8", "T7", "C3", "Cz", "C4", "T8",
                         "P7", "P3", "Pz", "P4", "P8", "O1", "Oz", "O2"], sfreq=SFREQ)

print(f"For each source depth, the best-matching superficial patch (centred at {PATCH_DEPTH:g} R depth):")
print(f"{'source depth (mm)':>19s} {'peak (uV)':>11s} | {'21 electrodes: best r':>23s} {'radius (deg)':>13s} "
      f"| {'61 electrodes: best r':>23s} {'radius (deg)':>13s}")
deep_rows = []
for frac in (0.2, 0.3, 0.4, 0.5, 0.6, 0.7):
    depth = frac * L5.HEAD_RADIUS_M * 1000
    row = [depth]
    for label, inf in (("21", info21), ("61", info)):
        target = patch_map(CENTRE, 0.0, frac, bem, inf)
        best_r, best_rad = -2.0, np.nan
        for rad in RADII:
            m = patch_map(CENTRE, float(rad), PATCH_DEPTH, bem, inf)
            r = float(np.corrcoef(target - target.mean(), m - m.mean())[0, 1])
            if r > best_r:
                best_r, best_rad = r, float(rad)
        row += [best_r, best_rad]
        if label == "21":
            row.append(float(np.abs(target).max()))
    deep_rows.append(row)
    print(f"{row[0]:19.1f} {row[3]:11.3f} | {row[1]:23.4f} {row[2]:13.1f} | {row[4]:23.4f} {row[5]:13.1f}")
print()
print("Every depth has a superficial patch whose map is nearly the same.  More electrodes help -- the")
print("correlations at 61 are below those at 21 -- but they do not solve it, and neither number is anywhere")
print("near low enough to call the two configurations distinguishable in the presence of real noise.")

fig, axes = plt.subplots(1, 3, figsize=(14.5, 4.0))
frac = 0.5
target21 = patch_map(CENTRE, 0.0, frac, bem, info21)
best = max(RADII, key=lambda rad: float(np.corrcoef(target21 - target21.mean(),
                                                    patch_map(CENTRE, float(rad), PATCH_DEPTH, bem, info21) -
                                                    patch_map(CENTRE, float(rad), PATCH_DEPTH, bem,
                                                              info21).mean())[0, 1]))
match21 = patch_map(CENTRE, float(best), PATCH_DEPTH, bem, info21)
scale = np.abs(target21).max() / np.abs(match21).max()
v = np.abs(target21).max()
for ax, (vals, title) in zip(axes[:2], ((target21, f"one dipole {frac * L5.HEAD_RADIUS_M * 1000:.0f} mm deep"),
                                        (match21 * scale,
                                         f"{best:.0f}-degree superficial patch, rescaled"))):
    im, _ = mne.viz.plot_topomap(vals, info21, axes=ax, show=False, contours=6, sensors=True,
                                 vlim=(-v, v))
    ax.set_title(title, fontsize=9)
cb = fig.colorbar(im, ax=axes[:2], shrink=0.8); cb.set_label("uV")
for k, label in ((1, "21 electrodes"), (4, "61 electrodes")):
    axes[2].plot([r[0] for r in deep_rows], [r[k] for r in deep_rows], "-o", ms=5, label=label)
axes[2].set_xlabel("source depth below the scalp (mm)")
axes[2].set_ylabel("best correlation with a superficial patch (dimensionless)")
axes[2].set_ylim(0.8, 1.005); axes[2].grid(alpha=0.25); axes[2].legend(fontsize=8)
axes[2].set_title("how well a superficial patch can imitate a deep source", fontsize=9)
fig.suptitle("The measurement behind the exercise: depth and spatial extent trade off against each other, "
             "and the scalp cannot separate them", y=1.03, fontsize=10)
plt.show()   # render the static figure(s) of this cell inline

print("nb-5-5-inverse -- L5.5 numbers (draft; TODO(confirm) at author review)")
print(f"Model: four-shell concentric sphere ({L5.SPHERE_CALLS['four-shell']}), {len(MONTAGE)} electrodes at "
      f"standard_1005 template positions,")
print(f"       {len(rr)} radial sources on a {GRID_SPACING_M * 1000:.0f} mm grid.  NO dataset, NO anatomy, "
      f"NO cortical surface -- see section 0.")
print(f"Simulation: one radial dipole, {MOMENT_nAm:.0f} nA.m at {F_HZ:g} Hz, {N_EPOCHS} epochs, "
      f"{NOISE_UV:g} uV sensor noise, seed {L5.SEED}")
print(f"Regularization: lambda^2 = {LAMBDA2:.4f} (SNR {SNR:g}); noise covariance ad hoc, "
      f"{NOISE_UV:g} uV per electrode")
print()
print(f"ANSWER KEY -- ill-posedness: rank(G) = {rank} with {G.shape[1]} source dimensions, so the null space "
      f"has {G.shape[1] - rank}")
print(f"    dimensions.  A unit-norm null-space source pattern produces {NULL_REL:.1e} of the leadfield's "
      f"largest entry at any")
print(f"    electrode -- nothing.  Largest / smallest singular value of G = {S[0] / S[rank - 1]:.1f}.")
print()
print(f"ANSWER KEY -- five estimators on one simulated source {depth_mm[TRUE]:.0f} mm below the scalp:")
print(f"{'':6s}{'method':>10s} {'peak error (mm)':>17s} {'dispersion (mm)':>17s} {'estimated depth (mm)':>21s}")
for method, err, disp, pd in est_rows:
    print(f"{'':6s}{method:>10s} {err:17.1f} {disp:17.1f} {pd:21.1f}")
print()
print("ANSWER KEY -- resolution over the whole grid (noiseless):")
print(f"{'':6s}{'method':>10s} {'peak error median (mm)':>24s} {'mean':>7s} {'dispersion median (mm)':>24s}")
for method, med, mean_e, disp in res_rows:
    print(f"{'':6s}{method:>10s} {med:24.1f} {mean_e:7.1f} {disp:24.1f}")
print("      sLORETA and eLORETA reach zero median peak error BY CONSTRUCTION.  It is a property of the")
print("      estimator, not evidence about the data, and their dispersion is no better than MNE's.")
print()
print("ANSWER KEY -- depth weighting (minimum-norm, gamma sweep):")
print(f"{'':6s}{'gamma':>7s} {'peak depth, shallow (mm)':>26s} {'peak depth, deep (mm)':>23s} "
      f"{'median peak error (mm)':>24s} {'median dispersion (mm)':>24s}")
for gamma, ds, dd, pe, disp in gamma_rows:
    print(f"{'':6s}{gamma:7.1f} {ds:26.1f} {dd:23.1f} {pe:24.1f} {disp:24.1f}")
print(f"      true depths {depth_mm[TRUE]:.0f} mm and {depth_mm[DEEP]:.0f} mm.  gamma is a choice, not a "
      f"measurement.")
print()
print("ANSWER KEY -- ex-5-5 (which claim is defensible: 'lateral parietal' or 'hippocampal'):")
print("    THE ANSWER IS 'LATERAL PARIETAL', AND THE REASON IS NOT ABOUT PARIETAL CORTEX.")
print("    For every source depth tested, a superficial patch exists whose scalp map is nearly identical:")
print(f"{'':6s}{'source depth (mm)':>19s} {'best r, 21 electrodes':>23s} {'best r, 61 electrodes':>23s} "
      f"{'peak (uV)':>11s}")
for row in deep_rows:
    print(f"{'':6s}{row[0]:19.1f} {row[1]:23.4f} {row[4]:23.4f} {row[3]:11.3f}")
print("    A superficial claim is defensible because a superficial generator is the configuration the data")
print("    ALREADY support -- every estimator above places its peak superficially unless a depth weighting is")
print("    turned up by hand.  A deep claim needs the data to EXCLUDE the superficial alternative, and the")
print("    table shows they cannot: the alternative reproduces the map to r > 0.9 at every depth, on both")
print("    montages.  More electrodes lower the correlation but never far enough.")
print("    A hippocampal claim from scalp EEG needs a constraint from outside the EEG -- an anatomical prior,")
print("    a simultaneous intracranial recording, or a lesion -- and a source map on its own is not one.")
print()
print("WHAT THIS NOTEBOOK CANNOT DO, AND WHY: no cortical surface, no anatomical parcellation and no region")
print("names, because ds-fsaverage and ds-mne-sample both fail the spec 10.7 licence gate (section 0).  Every")
print("location above is a coordinate in a sphere, and the notebook never calls one by an anatomical name.")
print("That is not a workaround: it is the same discipline the exercise is testing.")
print()
print("Pitfall: pf-deep-source-claims.")
print(L5.disk_line("disk at the end"))
print("Downloads in this notebook: none.")