#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
stt11_simulation.py
===================
Simulation Evidence Pack for Scale-Time Theory 11.0 — release 1.0

This pack is the constructive-evidence companion to the manuscript
"Scale-Time Theory 11.0 — The Pre-Geometry of Scale Space". Every claim
tagged [S:E1]..[S:E5] in the manuscript resolves to one of the five
experiments below; Appendix A of the manuscript summarizes them, and the
`verify` mode of this script machine-checks every number the manuscript
quotes against a fresh run.

Model (identical in every experiment):
  * a deterministic rotating phasor (or a superposed counter-rotating pair),
  * additive Gaussian readout noise, eta per quadrature,
  * stroboscopic sampling at R_A samples per standing-mode cycle.
No Hilbert space, hbar, Born rule, wavefunction, or quantum postulate is
inserted anywhere. Any quantum-like behavior shown in these tests arises
within this toy model from sampling and estimation alone.

Experiments and the manuscript claims they support:
  E1  Spin-direction statistics: statistical ambiguity zone around the
      structural edge, alias-reversal band (R_A < 2), stroboscopic freeze
      (phase advance 2*pi per sample; R_A = 1 demonstrated), and
      amplitude-squared-weighted outcomes for superposed modes.
      -> manuscript Sec. 8.1, 8.3, 8.5, 8.6, 9.
  E2  Resolution scaling and passive recurrence: per-cycle phase error is
      smooth and Cramer-Rao-limited in R_A (no privileged sampling values);
      the quantized-readout convention Delta_alpha = 2*pi/R_A; passive
      (kappa = 0) overlay locks only as measure-zero rational coincidences.
      -> manuscript Sec. 7.2, 8.2, 8.4.
  E3  Kinematic aliasing and resolvability thresholds: a detuned mode's
      residual precession turns from sign-flicker to coherent rotation;
      the 95%-stability threshold R*(O,A) moves with readout noise as the
      resolvability criterion predicts. -> manuscript Sec. 8.2, 8.3, 9.
  E4  Structural Nyquist identity: the sampled + and - spin states are the
      same state at R_A = 2 (overlap 1.000, discrimination error 0.500);
      frequency uncertainty tracks the N^(-3/2) bound at depth and blows up
      near the edge. -> manuscript Sec. 8.1, 8.6.
  E5  Baseline coupling: kappa = 0 gives no lock windows; kappa > 0 gives
      finite windows ordered by rational simplicity, noise-rigid inside and
      diffusive outside. -> manuscript Sec. 7.2, 7.3.

Notation follows the manuscript: eta is the per-quadrature readout-noise
scale (the archival report's sigma); a detuning delta corresponds to a
residual structure Delta_phi_res = 2*pi*delta per cycle.

Usage:
  python3 stt11_simulation.py all        # run E1..E5, build manuscript
                                         # figures M3-M5, write results.json,
                                         # then auto-verify
  python3 stt11_simulation.py e1 e3      # run a subset
  python3 stt11_simulation.py verify     # re-check results.json against the
                                         # manuscript-quoted values

Outputs: figs/figE1..figE5*.png, figs/figM3..figM5*.png, results.json.
Fixed seed 20599373; all numbers regenerate exactly.
"""

import json
import os
import sys
import time

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

SEED = 20599373
FIGDIR = "figs"
os.makedirs(FIGDIR, exist_ok=True)
PACK = {"name": "Scale-Time Theory 11.0 — Simulation Evidence Pack",
        "release": "1.0", "seed": SEED,
        "companion_to": "Scale-Time Theory 11.0 (manuscript)"}
RESULTS = {"pack": PACK}
CACHE = {}

plt.rcParams.update({
    "figure.dpi": 150, "savefig.dpi": 150, "font.size": 9,
    "axes.grid": True, "grid.alpha": 0.3, "axes.axisbelow": True,
    "legend.framealpha": 0.9,
})


# ----------------------------------------------------------------------------
# helpers
# ----------------------------------------------------------------------------
def wrap(x):
    """Wrap angle(s) to (-pi, pi]."""
    return (x + np.pi) % (2 * np.pi) - np.pi


def wilson_ci(p, n, z=1.96):
    """Wilson score interval for a binomial proportion."""
    denom = 1 + z * z / n
    centre = (p + z * z / (2 * n)) / denom
    half = z * np.sqrt(p * (1 - p) / n + z * z / (4 * n * n)) / denom
    return centre - half, centre + half


def gen_phasor(R, m_cycles, eta, trials, rng, direction=1.0, delta=0.0,
               chi=None, n_samples=None):
    """
    Stroboscopically sampled phasor (the framework's minimal observable).

    Underlying deterministic signal at sweep ticks k = 0..N-1, with R samples
    per nominal mode cycle:
        theta_k = direction * 2*pi*(1+delta)/R * k + phi0
        s_k     = exp(i theta_k)                         (traveling mode)
        s_k     = cos(chi) e^{+i theta} + sin(chi) e^{-i theta}
                                                         (superposed modes)
    plus additive complex Gaussian readout noise, eta per quadrature.
    phi0 is uniform per trial: the observer does not know the initial phase.
    Returns (z, phi0) with z of shape (trials, N).
    """
    if n_samples is None:
        n_samples = max(int(round(m_cycles * R)), 3)
    k = np.arange(n_samples)[None, :]
    phi0 = rng.uniform(0, 2 * np.pi, size=(trials, 1))
    step = direction * 2 * np.pi * (1 + delta) / R
    theta = k * step + phi0
    if chi is None:
        s = np.exp(1j * theta)
    else:
        s = np.cos(chi) * np.exp(1j * theta) + np.sin(chi) * np.exp(-1j * theta)
    noise = eta * (rng.standard_normal(s.shape)
                   + 1j * rng.standard_normal(s.shape))
    return s + noise, phi0


def direction_statistic(z):
    """Sum of lag-1 quadrature products; its sign is the observer's estimate
    of the rotation direction (the spin sign / normal direction n_hat)."""
    return np.imag(np.conj(z[:, :-1]) * z[:, 1:]).sum(axis=1)


def mark_edge(ax, label=True):
    """Mark the one structural node of the framework: R_A = 2."""
    lo, hi = ax.get_xlim()
    if lo <= 2 <= hi:
        ax.axvline(2, color="k", ls=":", lw=0.9, alpha=0.7)
        if label:
            ax.text(2, ax.get_ylim()[0] + 0.98 * np.diff(ax.get_ylim())[0],
                    "$R_A{=}2$\n(structural edge)", ha="center", va="top",
                    fontsize=7, alpha=0.8)


# ----------------------------------------------------------------------------
# E1  Spin-direction statistics (manuscript Sec. 8.1, 8.3, 8.5, 8.6, 9)
# ----------------------------------------------------------------------------
def experiment_1():
    """Statistical ambiguity zone, alias-reversal band, stroboscopic freeze,
    and amplitude-squared-weighted outcomes. Supports the manuscript tags
    [S:E1] in Sec. 8.1 (statistical zone, reversal, freeze), 8.5
    (Born-adjacent outcomes), 8.6 and 9 (regime chain)."""
    rng = np.random.default_rng(SEED + 1)
    t0 = time.time()
    res = {}

    # --- (a) P(correct spin sign) vs R_A for two readout-noise levels -------
    # NOTE: the R-grid below is a dense scan; no grid point is privileged.
    Rs = np.unique(np.concatenate([
        np.linspace(1.05, 3.0, 45), np.geomspace(3.0, 256.0, 28),
        np.array([2.0, 128.0, 256.0])]))
    etas = [0.35, 0.7]
    m_cycles, trials = 3, 2500
    curves = {}
    for et in etas:
        P, lo, hi = [], [], []
        for R in Rs:
            z, _ = gen_phasor(R, m_cycles, et, trials, rng)
            stat = direction_statistic(z)
            p = (np.sum(stat > 0) + 0.5 * np.sum(stat == 0)) / trials
            c = wilson_ci(p, trials)
            P.append(p); lo.append(c[0]); hi.append(c[1])
        curves[et] = (np.array(P), np.array(lo), np.array(hi))

    idx2 = int(np.argmin(np.abs(Rs - 2.0)))
    res["P_correct_at_R2"] = {f"eta={e}": float(curves[e][0][idx2])
                              for e in etas}
    r95 = {}
    for et in etas:
        P = curves[et][0]
        above = np.where((Rs > 2) & (P >= 0.95))[0]
        r95[f"eta={et}"] = float(Rs[above[0]]) if len(above) else None
    res["R_at_P95"] = r95

    # --- (b) aliased apparent frequency: reversal band and freeze -----------
    nu_true = np.linspace(0.02, 1.55, 65)          # mode cycles per sweep tick
    n_samp, tr_b, eta_b = 64, 600, 0.3
    mean_nu, std_nu = [], []
    for nu in nu_true:
        R = 1.0 / nu
        z, _ = gen_phasor(R, None, eta_b, tr_b, rng, n_samples=n_samp)
        lag = np.mean(np.conj(z[:, :-1]) * z[:, 1:], axis=1)
        nu_hat = np.angle(lag) / (2 * np.pi)
        mean_nu.append(np.mean(nu_hat)); std_nu.append(np.std(nu_hat))
    mean_nu, std_nu = np.array(mean_nu), np.array(std_nu)
    fold = ((nu_true + 0.5) % 1.0) - 0.5           # ideal alias fold

    # --- (c) superposed counter-rotating modes ------------------------------
    chis = np.linspace(0, np.pi / 2, 25)
    Rs_c, eta_c, tr_c = [2.2, 4.0, 32.0], 0.6, 4000
    pplus = {}
    born = np.cos(chis) ** 2
    for R in Rs_c:
        row = []
        for chi in chis:
            z, _ = gen_phasor(R, 3, eta_c, tr_c, rng, chi=chi)
            stat = direction_statistic(z)
            row.append((np.sum(stat > 0) + 0.5 * np.sum(stat == 0)) / tr_c)
        pplus[R] = np.array(row)
    res["max_dev_from_born"] = {f"R={R}": round(float(np.max(np.abs(pplus[R] - born))), 3)
                                for R in Rs_c}

    CACHE["e1"] = dict(Rs=Rs, curves=curves, etas=etas, nu_true=nu_true,
                       mean_nu=mean_nu, std_nu=std_nu, fold=fold,
                       chis=chis, pplus=pplus, born=born)

    # --- figure --------------------------------------------------------------
    fig, axes = plt.subplots(1, 3, figsize=(12.5, 3.7), constrained_layout=True)

    ax = axes[0]
    cols = {etas[0]: "tab:blue", etas[1]: "tab:red"}
    for et in etas:
        P, lo, hi = curves[et]
        ax.fill_between(Rs, lo, hi, color=cols[et], alpha=0.2)
        ax.plot(Rs, P, color=cols[et], lw=1.4, label=fr"readout noise $\eta={et}$")
    ax.axhline(0.5, color="gray", lw=0.8)
    ax.axvspan(Rs.min(), 2.0, color="orange", alpha=0.12)
    ax.text(1.35, 0.06, "alias-reversal\nband", fontsize=7,
            ha="center", color="darkorange")
    ax.set_xscale("log")
    ax.set_xlabel(r"sampling ratio $R_A$ (samples per mode cycle)")
    ax.set_ylabel("P(correct spin sign)")
    ax.set_ylim(-0.03, 1.05)
    ax.set_title("(a) statistical ambiguity zone around the structural edge")
    ax.legend(loc="center right", fontsize=7)
    mark_edge(ax)

    ax = axes[1]
    ax.plot(nu_true, fold, "k--", lw=1, label="ideal alias fold")
    ax.errorbar(nu_true, mean_nu, yerr=std_nu, fmt=".", ms=3, lw=0.7,
                color="tab:purple", label="estimated (MC)")
    ax.axvline(0.5, color="k", ls=":", lw=0.8)
    ax.axvline(1.0, color="k", ls=":", lw=0.8)
    ax.text(0.5, 0.55, "structural\nedge\n$R_A{=}2$", ha="center", fontsize=7)
    ax.text(1.0, 0.55, "stroboscopic\nfreeze\n$R_A{=}1$", ha="center", fontsize=7)
    ax.set_xlabel(r"true mode frequency $\nu = 1/R_A$ (cycles per tick)")
    ax.set_ylabel(r"apparent frequency $\hat\nu$")
    ax.set_title("(b) aliased appearance: reversal and freeze")
    ax.legend(fontsize=7, loc="lower left")

    ax = axes[2]
    for R, c in zip(Rs_c, ["tab:orange", "tab:green", "tab:blue"]):
        ax.plot(np.degrees(chis), pplus[R], "-o", ms=2.5, lw=1.2, color=c,
                label=fr"$R_A={R}$")
    ax.plot(np.degrees(chis), born, "k--", lw=1,
            label=r"$\cos^2\chi$ (Born-rule form, reference)")
    ax.set_xlabel(r"superposition angle $\chi$  ($0$: pure $+$, $45°$: standing)")
    ax.set_ylabel("P(outcome $+$)")
    ax.set_title("(c) Born-adjacent squared-amplitude ordering")
    ax.legend(fontsize=7)

    fig.savefig(f"{FIGDIR}/figE1_spin_direction.png")
    plt.close(fig)
    res["runtime_s"] = round(time.time() - t0, 1)
    RESULTS["E1_spin_direction"] = res
    print("E1 done", res)


# ----------------------------------------------------------------------------
# E2  Resolution scaling and passive recurrence (manuscript Sec. 7.2, 8.2, 8.4)
# ----------------------------------------------------------------------------
def experiment_2():
    """Per-cycle phase error is smooth and CRLB-limited in R_A, with no
    privileged sampling values; the quantized-readout convention
    Delta_alpha = 2*pi/R_A is illustrated as a convention; passive
    (kappa = 0) overlay recurrence locks only at measure-zero rational
    coincidences. Supports [S:E2] in Sec. 7.2 and the smoothness content
    of Sec. 8.2; illustrates the [D] convention of Sec. 8.4."""
    rng = np.random.default_rng(SEED + 2)
    t0 = time.time()
    res = {}

    # --- (a) phase RMSE vs R_A (continuous scan) -----------------------------
    Rs = np.geomspace(2, 512, 40)
    m_cycles, eta, trials = 2, 0.4, 2000
    rmse, rmse_q, crlb = [], [], []
    for R in Rs:
        z, phi0 = gen_phasor(R, m_cycles, eta, trials, rng)
        N = z.shape[1]
        k = np.arange(N)
        template = np.exp(-1j * 2 * np.pi * k / R)[None, :]
        phi_hat = np.angle((z * template).sum(axis=1))
        err = wrap(phi_hat - phi0[:, 0])
        rmse.append(np.sqrt(np.mean(err ** 2)))
        binw = 2 * np.pi / R                       # readout convention (Sec. 8.4)
        phi_q = np.round(phi_hat / binw) * binw
        err_q = wrap(phi_q - phi0[:, 0])
        rmse_q.append(np.sqrt(np.mean(err_q ** 2)))
        crlb.append(eta / np.sqrt(N))
    rmse, rmse_q, crlb = map(np.array, (rmse, rmse_q, crlb))

    slope = np.polyfit(np.log(Rs), np.log(rmse), 1)[0]
    res["intrinsic_loglog_slope"] = round(float(slope), 3)
    fit = np.exp(np.polyval(np.polyfit(np.log(Rs), np.log(rmse), 1), np.log(Rs)))
    resid = rmse / fit - 1
    res["max_abs_residual_vs_smooth_fit"] = round(float(np.max(np.abs(resid))), 3)

    # --- (b) passive overlay recurrence in a neutral window -------------------
    # kappa = 0: recurrence is exact only at rational R_A (measure-zero locks).
    Rg = np.linspace(12.0, 24.0, 12001)
    Q = 8
    drifts = np.full_like(Rg, np.inf)
    for q in range(1, Q + 1):
        d = np.abs(wrap(2 * np.pi * np.round(q * Rg) / Rg))
        drifts = np.minimum(drifts, d)
    res["rational_locks_int_and_half_in_window_12_24"] = int(np.sum(
        (drifts < 1e-3) & np.isclose(Rg % 0.5, 0, atol=1e-9)))

    # --- figure ----------------------------------------------------------------
    fig, axes = plt.subplots(1, 2, figsize=(10.5, 3.7), constrained_layout=True)

    ax = axes[0]
    ax.loglog(Rs, rmse, "o-", ms=3, lw=1.2, color="tab:blue",
              label="per-cycle phase error (MC)")
    ax.loglog(Rs, crlb, "k--", lw=1, label=r"CRLB $\eta/\sqrt{N}$")
    ax.loglog(Rs, rmse_q, "s-", ms=3, lw=1.2, color="tab:red",
              label=r"with quantized readout $\Delta\alpha=2\pi/R_A$")
    ax.loglog(Rs, (2 * np.pi / Rs) / np.sqrt(12), ":", color="tab:red", lw=1,
              label=r"bin floor $\Delta\alpha/\sqrt{12}$ (convention)")
    ax.set_xlabel(r"$R_A$")
    ax.set_ylabel("phase RMSE [rad]")
    ax.set_title("(a) resolution is smooth: no privileged sampling values")
    ax.legend(fontsize=7)
    mark_edge(ax)

    ax = axes[1]
    ax.plot(Rg, drifts, lw=0.6, color="tab:green")
    ax.set_yscale("log")
    ax.set_xlabel(r"$R_A$")
    ax.set_ylabel(r"best overlay drift $\min_{m\leq 8}|\Delta\psi_m|$ [rad]")
    ax.set_title(r"(b) passive recurrence ($\kappa=0$): measure-zero rational locks")
    fig.savefig(f"{FIGDIR}/figE2_resolution_passive.png")
    plt.close(fig)
    res["runtime_s"] = round(time.time() - t0, 1)
    RESULTS["E2_resolution_passive"] = res
    print("E2 done", res)


# ----------------------------------------------------------------------------
# E3  Kinematic aliasing and resolvability thresholds (Sec. 8.2, 8.3, 9)
# ----------------------------------------------------------------------------
def experiment_3():
    """A detuned mode (residual structure Delta_phi_res = 2*pi*delta per
    cycle) shows sign-flickering precession near the edge and coherent
    rotation at depth; the 95%-stability threshold R*(O,A) moves with
    readout noise eta as the resolvability criterion predicts. Supports
    [S:E3] in Sec. 8.2 (criterion and threshold), 8.3 (structure-relativity)
    and 9 (coherent aliasing -> stable n_hat). Also produces manuscript
    Figure 4 (the kinematic-aliasing triptych)."""
    rng = np.random.default_rng(SEED + 3)
    t0 = time.time()
    res = {}
    delta = 0.02
    true_inc = 2 * np.pi * delta                    # 0.1257 rad per cycle
    Rs = np.array([2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256,
                   384, 512])
    etas = [0.2, 0.45, 0.8]
    m_cycles, trials = 24, 250
    S = {et: [] for et in etas}                     # direction stability
    J = {et: [] for et in etas}                     # increment jitter / signal

    for et in etas:
        for R in Rs:
            n = int(R)
            z, _ = gen_phasor(R, m_cycles, et, trials, rng, delta=delta,
                              n_samples=m_cycles * n)
            zc = z.reshape(trials, m_cycles, n)
            template = np.exp(-1j * 2 * np.pi * np.arange(n) / R)[None, None, :]
            phi = np.angle((zc * template).sum(axis=2))          # (trials, m)
            d = wrap(np.diff(phi, axis=1))                       # increments
            S[et].append(np.abs(np.mean(np.sign(d))))
            J[et].append(np.std(d) / true_inc)

    rstar = {}
    for et in etas:
        s = np.array(S[et])
        above = np.where(s >= 0.95)[0]
        if len(above) and above[0] > 0:
            i = above[0]
            x0, x1 = np.log(Rs[i - 1]), np.log(Rs[i])
            y0, y1 = s[i - 1], s[i]
            rstar[f"eta={et}"] = round(float(np.exp(
                x0 + (0.95 - y0) * (x1 - x0) / (y1 - y0))), 1)
        elif len(above):
            rstar[f"eta={et}"] = float(Rs[0])
        else:
            rstar[f"eta={et}"] = None
    res["Rstar_S95_vs_noise"] = rstar
    res["residual_structure_rad_per_cycle"] = round(true_inc, 4)

    # --- manuscript Figure 4: the kinematic-aliasing triptych -------------------
    fig2, axs = plt.subplots(2, 3, figsize=(11.5, 6.6), constrained_layout=True)
    show = [2.0, 8.0, 64.0]
    eta_t, m_t, dl_t = 0.15, 8, 0.03
    for j, R in enumerate(show):
        n = max(int(round(R)), 2)
        z, _ = gen_phasor(R, m_t, eta_t, 1, rng, delta=dl_t,
                          n_samples=m_t * n)
        zz = z[0]
        ax = axs[0, j]
        pts = ax.scatter(zz.real, zz.imag, c=np.arange(len(zz)),
                         cmap="viridis", s=8)
        ax.plot(zz.real, zz.imag, lw=0.4, alpha=0.5, color="gray")
        th = np.linspace(0, 2 * np.pi, 200)
        ax.plot(np.cos(th), np.sin(th), "k:", lw=0.6)
        ax.set_aspect("equal")
        ax.set_title(fr"$R_A={R:g}$", fontsize=9)
        ax.set_xticks([]); ax.set_yticks([])
        if j == 2:
            fig2.colorbar(pts, ax=ax, label="sweep tick", shrink=0.8)

        zc = zz.reshape(m_t, n)
        template = np.exp(-1j * 2 * np.pi * np.arange(n) / R)[None, :]
        phi = np.unwrap(np.angle((zc * template).sum(axis=1)))
        ax2 = axs[1, j]
        ax2.plot(phi - phi[0], "o-", ms=3, lw=1, color="tab:blue",
                 label=r"estimated $\hat\varphi_m$")
        ax2.plot(np.arange(m_t) * 2 * np.pi * dl_t * (n / R), "k--", lw=1,
                 label="true residual precession")
        ax2.set_xlabel("cycle $m$")
        if j == 0:
            ax2.set_ylabel("accumulated residual phase [rad]")
            ax2.legend(fontsize=7)
    axs[0, 0].set_ylabel("sampled trajectory\n(complex plane)")
    fig2.suptitle("Kinematic aliasing: identical dynamics at three sampling "
                  f"depths ($\\delta={dl_t}$, $\\eta={eta_t}$)", fontsize=10)
    fig2.savefig(f"{FIGDIR}/figM4_kinematic_aliasing.png")
    plt.close(fig2)

    CACHE["e3"] = dict(Rs=Rs, S=S, J=J, etas=etas, rstar=rstar,
                       true_inc=true_inc)

    # --- main threshold figure ---------------------------------------------------
    fig, axes = plt.subplots(1, 2, figsize=(10.5, 3.7), constrained_layout=True)
    cols = plt.cm.plasma(np.linspace(0.15, 0.8, len(etas)))
    ax = axes[0]
    for et, c in zip(etas, cols):
        lab = fr"$\eta={et}$"
        if rstar[f"eta={et}"]:
            lab += fr", $R^*\approx${rstar[f'eta={et}']}"
        ax.semilogx(Rs, S[et], "o-", ms=3.5, lw=1.3, color=c, label=lab)
        if rstar[f"eta={et}"]:
            ax.axvline(rstar[f"eta={et}"], color=c, ls="--", lw=0.8, alpha=0.6)
    ax.axhline(0.95, color="gray", ls="--", lw=0.8)
    ax.set_xlabel(r"$R_A$")
    ax.set_ylabel(r"normal-direction stability $S=|\langle\mathrm{sign}\,"
                  r"\Delta\hat\varphi\rangle|$")
    ax.set_ylim(0, 1.05)
    ax.set_title(r"(a) resolvability thresholds $R^*(O,A)$ move with noise")
    ax.legend(fontsize=7, loc="lower right")
    mark_edge(ax)

    ax = axes[1]
    for et, c in zip(etas, cols):
        ax.loglog(Rs, J[et], "s-", ms=3.5, lw=1.3, color=c, label=fr"$\eta={et}$")
    ax.axhline(1.0, color="gray", ls="--", lw=0.8)
    ax.text(Rs[0], 1.1, "jitter = signal", fontsize=7, color="gray")
    ax.set_xlabel(r"$R_A$")
    ax.set_ylabel("precession jitter / residual structure")
    ax.set_title(r"(b) coherence follows $\varepsilon_\phi \sim \eta/\sqrt{R_A}$")
    ax.legend(fontsize=7)
    mark_edge(ax, label=False)
    fig.savefig(f"{FIGDIR}/figE3_resolvability_thresholds.png")
    plt.close(fig)
    res["runtime_s"] = round(time.time() - t0, 1)
    RESULTS["E3_resolvability"] = res
    print("E3 done", res)


# ----------------------------------------------------------------------------
# E4  Structural Nyquist identity and uncertainty (manuscript Sec. 8.1, 8.6)
# ----------------------------------------------------------------------------
def experiment_4():
    """(a) The sampled + and - spin states are the SAME state at R_A = 2:
    overlap 1.000, minimum discrimination error 0.500. Structural, not
    statistical, ambiguity. (b) Frequency estimation tracks the single-tone
    N^(-3/2) bound at depth and blows up near the edge where the +/- alias
    hypotheses merge. Supports [S:E4] in Sec. 8.1 and the regime chain."""
    rng = np.random.default_rng(SEED + 4)
    t0 = time.time()
    res = {}
    m_cycles, eta, trials = 3, 0.5, 3000
    Rs = np.unique(np.concatenate([np.linspace(2.0, 8.0, 30),
                                   np.geomspace(8.0, 32.0, 8)]))

    # --- (a) +/- discrimination error and analytic overlap ---------------------
    err, overlap = [], []
    for R in Rs:
        z, _ = gen_phasor(R, m_cycles, eta, trials, rng, direction=+1.0)
        N = z.shape[1]
        k = np.arange(N)
        tp = np.exp(-1j * 2 * np.pi * k / R)[None, :]
        tm = np.exp(+1j * 2 * np.pi * k / R)[None, :]
        Ep = np.abs((z * tp).sum(axis=1)) ** 2
        Em = np.abs((z * tm).sum(axis=1)) ** 2
        ties = np.isclose(Ep, Em)
        e = (np.sum((Em > Ep) & ~ties) + 0.5 * np.sum(ties)) / trials
        err.append(e)
        overlap.append(np.abs(np.exp(-1j * 4 * np.pi * k / R).sum()) / N)
    err, overlap = np.array(err), np.array(overlap)
    res["error_at_R2"] = round(float(err[0]), 3)
    res["overlap_at_R2"] = round(float(overlap[0]), 3)
    below5 = np.where(err < 0.05)[0]
    res["R_error_below_5pct"] = round(float(Rs[below5[0]]), 2) if len(below5) else None

    # --- (b) frequency RMSE via full-band ML search -----------------------------
    Rs_b = np.geomspace(2.05, 64, 16)
    tr_b = 400
    nu_grid = np.linspace(-0.5, 0.5, 801, endpoint=False)
    f_rmse, crlb_f = [], []
    for R in Rs_b:
        z, _ = gen_phasor(R, m_cycles, eta, tr_b, rng, direction=+1.0)
        N = z.shape[1]
        k = np.arange(N)
        E = np.exp(-2j * np.pi * np.outer(nu_grid, k))          # (F, N)
        P = np.abs(z @ E.T) ** 2                                 # (trials, F)
        i0 = np.argmax(P, axis=1)
        i0 = np.clip(i0, 1, len(nu_grid) - 2)
        pm, p0, pp = (P[np.arange(tr_b), i0 - 1],
                      P[np.arange(tr_b), i0],
                      P[np.arange(tr_b), i0 + 1])
        denom = pm - 2 * p0 + pp
        shift = np.where(np.abs(denom) > 0, 0.5 * (pm - pp) / denom, 0.0)
        shift = np.clip(shift, -0.5, 0.5)
        nu_hat = nu_grid[i0] + shift * (nu_grid[1] - nu_grid[0])
        f_rmse.append(np.sqrt(np.mean((nu_hat - 1.0 / R) ** 2)))
        crlb_f.append(np.sqrt(6 * (2 * eta ** 2)
                              / ((2 * np.pi) ** 2 * N * (N ** 2 - 1))))
    f_rmse, crlb_f = np.array(f_rmse), np.array(crlb_f)

    fig, axes = plt.subplots(1, 2, figsize=(10.5, 3.7), constrained_layout=True)
    ax = axes[0]
    ax.plot(Rs, err, "o-", ms=3, lw=1.2, color="tab:blue",
            label="MC discrimination error")
    ax.plot(Rs, overlap, "k--", lw=1.1,
            label=r"state overlap $|\langle s_+,s_-\rangle|/N$")
    ax.axhline(0.5, color="gray", lw=0.8)
    ax.set_xlabel(r"$R_A$")
    ax.set_ylabel("min. error / overlap")
    ax.set_title(r"(a) structural identity: $s_+\equiv s_-$ at $R_A{=}2$")
    ax.legend(fontsize=7)

    ax = axes[1]
    ax.loglog(Rs_b, f_rmse, "o-", ms=3.5, lw=1.2, color="tab:red",
              label=r"RMSE$(\hat\nu)$, full-band ML")
    ax.loglog(Rs_b, crlb_f, "k--", lw=1, label=r"single-tone bound $\propto N^{-3/2}$")
    ax.set_xlabel(r"$R_A$")
    ax.set_ylabel(r"frequency error [cycles/tick]")
    ax.set_title("(b) uncertainty blow-up near the structural edge")
    ax.legend(fontsize=7)
    fig.savefig(f"{FIGDIR}/figE4_structural_identity.png")
    plt.close(fig)
    res["runtime_s"] = round(time.time() - t0, 1)
    RESULTS["E4_structural_identity"] = res
    print("E4 done", res)


# ----------------------------------------------------------------------------
# E5  Baseline coupling and lock windows (manuscript Sec. 7.2, 7.3)
# ----------------------------------------------------------------------------
def experiment_5():
    """kappa = 0 gives no lock windows (winding follows the bare detuning);
    kappa > 0 gives finite lock windows ordered by rational simplicity,
    noise-rigid inside and diffusive outside. Minimal model: the circle map
    theta' = theta + Omega + (kappa/2pi) sin(2pi theta). Supports [S:E5]
    in Sec. 7.2 and 7.3; produces manuscript Figure 5."""
    rng = np.random.default_rng(SEED + 5)
    t0 = time.time()
    res = {}

    def winding(omega, K, n_iter=2500, n_trans=500):
        th = np.zeros_like(omega)
        for _ in range(n_trans):
            th = th + omega + (K / (2 * np.pi)) * np.sin(2 * np.pi * th)
        th0 = th.copy()
        for _ in range(n_iter):
            th = th + omega + (K / (2 * np.pi)) * np.sin(2 * np.pi * th)
        return (th - th0) / n_iter

    Om = np.linspace(0, 1, 2401)
    W_free = winding(Om.copy(), 0.0)               # kappa = 0: passive
    W_lock = winding(Om.copy(), 0.95)              # kappa > 0: coupled

    rationals = [(0, 1), (1, 4), (1, 3), (1, 2), (2, 3), (3, 4), (1, 1)]
    widths = {}
    for p, q in rationals:
        sel = np.abs(W_lock - p / q) < 1e-4
        widths[f"{p}/{q}"] = round(float(np.sum(sel) * (Om[1] - Om[0])), 4)
    res["lock_window_widths_kappa0.95"] = widths
    res["locked_fraction_kappa0.95"] = round(float(np.mean(
        np.min(np.abs(W_lock[:, None]
                      - np.array([p / q for p, q in rationals])[None, :]),
               axis=1) < 1e-4)), 3)

    def diffuse(omega0, K, noise, n_steps, n_ens):
        th = rng.uniform(0, 1, n_ens)
        traj = np.empty((n_steps, n_ens))
        for i in range(n_steps):
            th = th + omega0 + (K / (2 * np.pi)) * np.sin(2 * np.pi * th) \
                 + noise * rng.standard_normal(n_ens)
            traj[i] = th
        return traj

    n_steps, n_ens, drive_noise, K = 400, 500, 0.03, 0.95
    golden = (np.sqrt(5) - 1) / 2
    tr_lock = diffuse(0.5, K, drive_noise, n_steps, n_ens)   # inside 1/2 window
    tr_free = diffuse(golden, K, drive_noise, n_steps, n_ens)
    spread_lock = np.std(tr_lock - tr_lock.mean(axis=1, keepdims=True), axis=1)
    spread_free = np.std(tr_free - tr_free.mean(axis=1, keepdims=True), axis=1)
    res["final_phase_spread_locked"] = round(float(spread_lock[-1]), 3)
    res["final_phase_spread_unlocked"] = round(float(spread_free[-1]), 3)

    CACHE["e5"] = dict(Om=Om, W_free=W_free, W_lock=W_lock,
                       rationals=rationals, widths=widths)

    fig, axes = plt.subplots(1, 3, figsize=(12.5, 3.7), constrained_layout=True)
    ax = axes[0]
    ax.plot(Om, W_free, lw=0.7, color="gray", label=r"$\kappa=0$ (passive)")
    ax.plot(Om, W_lock, lw=0.9, color="tab:blue", label=r"$\kappa=0.95$ (coupled)")
    for p, q in rationals:
        ax.axhline(p / q, color="k", lw=0.3, alpha=0.3)
    ax.set_xlabel(r"bare detuning $\Omega$ (mode freq / sweep freq)")
    ax.set_ylabel("winding number $W$ (locked appearance)")
    ax.set_title("(a) baseline coupling creates lock windows")
    ax.legend(fontsize=7)

    ax = axes[1]
    qs = np.array([q for _, q in rationals if q > 0])
    ws = np.array([widths[f"{p}/{q}"] for p, q in rationals if q > 0])
    order = np.argsort(qs)
    ax.semilogy(qs[order], np.maximum(ws[order], 1e-5), "o-", color="tab:green")
    ax.set_xlabel("rational complexity $q$ (denominator of lock ratio $p/q$)")
    ax.set_ylabel("lock-window width")
    ax.set_title("(b) window width follows rational simplicity")

    ax = axes[2]
    ax.plot(spread_free, lw=1.2, color="tab:red",
            label=r"unlocked ($\Omega$ = golden mean)")
    ax.plot(spread_lock, lw=1.2, color="tab:blue",
            label=r"locked ($\Omega = 1/2$ window)")
    ax.set_xlabel("sweep cycles")
    ax.set_ylabel("unwrapped phase-drift spread")
    ax.set_title("(c) locked: noise-rigid; unlocked: diffusive")
    ax.legend(fontsize=7)
    fig.savefig(f"{FIGDIR}/figE5_lock_windows.png")
    plt.close(fig)
    res["runtime_s"] = round(time.time() - t0, 1)
    RESULTS["E5_lock_windows"] = res
    print("E5 done", res)


# ----------------------------------------------------------------------------
# Manuscript-ready figures M3 and M5 (composed from cached experiment data)
# ----------------------------------------------------------------------------
def manuscript_figures():
    """Compose the data-grounded manuscript figures from E1/E3/E5 caches
    (no additional random draws, so all numbers stay identical):
      figM3 -- sampling-regime diagram (manuscript Figure 3)
      figM5 -- lock windows, single panel (manuscript Figure 5)
    figM4 (kinematic-aliasing triptych) is produced inside E3."""
    if not all(k in CACHE for k in ("e1", "e3", "e5")):
        print("manuscript figures skipped (need e1, e3, e5 in this run)")
        return

    # --- Figure M3: sampling-regime diagram --------------------------------------
    e1, e3 = CACHE["e1"], CACHE["e3"]
    fig, axes = plt.subplots(1, 2, figsize=(11.0, 3.8), constrained_layout=True)

    ax = axes[0]
    ax.plot(e1["nu_true"], e1["fold"], "k--", lw=1, label="ideal alias fold")
    ax.errorbar(e1["nu_true"], e1["mean_nu"], yerr=e1["std_nu"], fmt=".",
                ms=3, lw=0.7, color="tab:purple", label="estimated (MC)")
    ax.axvspan(0.5, 1.0, color="orange", alpha=0.10)
    ax.axvline(0.5, color="k", ls=":", lw=0.9)
    ax.axvline(1.0, color="k", ls=":", lw=0.9)
    ax.text(0.25, 0.58, "resolvable\nside", ha="center", fontsize=7)
    ax.text(0.75, 0.58, "alias-reversal\nband ($R_A<2$)", ha="center",
            fontsize=7, color="darkorange")
    ax.text(0.5, -0.62, "structural edge\n$R_A{=}2$", ha="center", fontsize=7)
    ax.text(1.0, -0.62, "stroboscopic freeze\n$R_A{=}1$", ha="center", fontsize=7)
    ax.set_ylim(-0.72, 0.72)
    ax.set_xlabel(r"true mode frequency $\nu = 1/R_A$ (cycles per tick)")
    ax.set_ylabel(r"apparent frequency $\hat\nu$")
    ax.set_title("(a) below and at the structural edge")
    ax.legend(fontsize=7, loc="upper left")

    ax = axes[1]
    cols = plt.cm.plasma(np.linspace(0.15, 0.8, len(e3["etas"])))
    for et, c in zip(e3["etas"], cols):
        key = f"eta={et}"
        lab = fr"$\eta={et}$"
        if e3["rstar"][key]:
            lab += fr", $R^*\approx${e3['rstar'][key]}"
        ax.semilogx(e3["Rs"], e3["S"][et], "o-", ms=3.5, lw=1.3, color=c,
                    label=lab)
        if e3["rstar"][key]:
            ax.axvline(e3["rstar"][key], color=c, ls="--", lw=0.8, alpha=0.6)
    ax.axhline(0.95, color="gray", ls="--", lw=0.8)
    ax.set_ylim(0, 1.05)
    ax.set_xlabel(r"$R_A$")
    ax.set_ylabel("normal-direction stability $S$")
    ax.set_title(r"(b) above the edge: structure-relative thresholds $R^*(O,A)$")
    ax.legend(fontsize=7, loc="lower right")
    mark_edge(ax)

    fig.suptitle("Sampling regimes: reversal, freeze, structural edge, and "
                 "resolvability (manuscript Figure 3)", fontsize=10)
    fig.savefig(f"{FIGDIR}/figM3_sampling_regimes.png")
    plt.close(fig)

    # --- Figure M5: lock windows, single panel ------------------------------------
    e5 = CACHE["e5"]
    fig, ax = plt.subplots(figsize=(6.2, 4.0), constrained_layout=True)
    ax.plot(e5["Om"], e5["W_free"], lw=0.8, color="gray",
            label=r"$\kappa=0$: passive, no windows")
    ax.plot(e5["Om"], e5["W_lock"], lw=1.1, color="tab:blue",
            label=r"$\kappa=0.95$: finite lock windows")
    for p, q in e5["rationals"]:
        ax.axhline(p / q, color="k", lw=0.3, alpha=0.25)
    for lbl, x, y in [("0/1", 0.06, 0.02), ("1/2", 0.5, 0.52),
                      ("1/3", 0.31, 0.355), ("2/3", 0.69, 0.685)]:
        ax.annotate(lbl, (x, y), fontsize=7, ha="center", alpha=0.8)
    ax.set_xlabel(r"bare detuning $\Omega$ (mode freq / sweep freq)")
    ax.set_ylabel("winding number $W$ (locked appearance)")
    ax.set_title("Baseline coupling turns rational coincidences into\n"
                 "finite, noise-robust lock windows (manuscript Figure 5)",
                 fontsize=9)
    ax.legend(fontsize=7, loc="upper left")
    fig.savefig(f"{FIGDIR}/figM5_lock_windows.png")
    plt.close(fig)
    print("manuscript figures M3, M5 written (M4 written by E3)")


# ----------------------------------------------------------------------------
# Verification against the numbers quoted in the STT 11.0 manuscript
# ----------------------------------------------------------------------------
MANUSCRIPT_CHECKS = [
    # (results path, expected, abs tolerance, quoted where)
    (("E1_spin_direction", "P_correct_at_R2", "eta=0.35"), 0.50, 0.03,
     "Appendix A / Sec. 8.1: coin flip at the structural edge"),
    (("E1_spin_direction", "P_correct_at_R2", "eta=0.7"), 0.50, 0.03,
     "Appendix A / Sec. 8.1: coin flip, noise-independent"),
    (("E1_spin_direction", "R_at_P95", "eta=0.35"), 2.16, 0.05,
     "Appendix A: stable by R_A ~ 2.2 (low noise)"),
    (("E1_spin_direction", "R_at_P95", "eta=0.7"), 2.60, 0.05,
     "Appendix A: stable by R_A ~ 2.6 (higher noise)"),
    (("E1_spin_direction", "max_dev_from_born", "R=2.2"), 0.151, 0.01,
     "Sec. 8.5: Born-adjacent near the edge"),
    (("E2_resolution_passive", "intrinsic_loglog_slope"), -0.50, 0.02,
     "Appendix A: smooth CRLB-limited scaling"),
    (("E3_resolvability", "Rstar_S95_vs_noise", "eta=0.2"), 20.9, 1.0,
     "Appendix A: R* ~ 21"),
    (("E3_resolvability", "Rstar_S95_vs_noise", "eta=0.45"), 99.9, 5.0,
     "Appendix A: R* ~ 100"),
    (("E3_resolvability", "Rstar_S95_vs_noise", "eta=0.8"), 330.1, 15.0,
     "Appendix A: R* ~ 330"),
    (("E3_resolvability", "residual_structure_rad_per_cycle"), 0.1257, 0.0005,
     "Appendix A: residual structure 0.126 rad per cycle"),
    (("E4_structural_identity", "overlap_at_R2"), 1.000, 1e-6,
     "Sec. 8.1 / Appendix A: overlap 1.000"),
    (("E4_structural_identity", "error_at_R2"), 0.500, 0.01,
     "Sec. 8.1 / Appendix A: discrimination error 0.500"),
    (("E5_lock_windows", "lock_window_widths_kappa0.95", "0/1"), 0.151, 0.003,
     "Appendix A: width 0.151 at q=1"),
    (("E5_lock_windows", "lock_window_widths_kappa0.95", "1/2"), 0.067, 0.003,
     "Appendix A: width 0.067 at q=2"),
    (("E5_lock_windows", "lock_window_widths_kappa0.95", "1/3"), 0.026, 0.003,
     "Appendix A: width 0.026 at q=3"),
    (("E5_lock_windows", "lock_window_widths_kappa0.95", "1/4"), 0.013, 0.003,
     "Appendix A: width 0.013 at q=4"),
    (("E5_lock_windows", "final_phase_spread_locked"), 0.302, 0.05,
     "Sec. 7.3: locked states noise-rigid"),
    (("E5_lock_windows", "final_phase_spread_unlocked"), 0.913, 0.08,
     "Sec. 7.3: unlocked states diffusive"),
]


def verify(results=None):
    """Check every manuscript-quoted number against the (re)generated results.
    Exits nonzero on any failure, so this can gate a release."""
    if results is None:
        with open("results.json") as fh:
            results = json.load(fh)
    n_fail = 0
    print("\n--- verification against STT 11.0 manuscript quotes ---")
    for path, expected, tol, where in MANUSCRIPT_CHECKS:
        node = results
        try:
            for key in path:
                node = node[key]
            ok = node is not None and abs(float(node) - expected) <= tol
        except (KeyError, TypeError):
            node, ok = "MISSING", False
        status = "PASS" if ok else "FAIL"
        if not ok:
            n_fail += 1
        print(f"[{status}] {'/'.join(path):58s} = {node} "
              f"(expected {expected} ± {tol})  <- {where}")
    print(f"--- {len(MANUSCRIPT_CHECKS) - n_fail}/{len(MANUSCRIPT_CHECKS)} "
          f"checks passed ---")
    return n_fail


# ----------------------------------------------------------------------------
def main():
    args = sys.argv[1:] or ["all"]
    if args == ["verify"]:
        sys.exit(1 if verify() else 0)
    todo = {"e1": experiment_1, "e2": experiment_2, "e3": experiment_3,
            "e4": experiment_4, "e5": experiment_5}
    keys = list(todo) if "all" in args else [a for a in args if a in todo]
    t0 = time.time()
    for kk in keys:
        todo[kk]()
    manuscript_figures()
    RESULTS["total_runtime_s"] = round(time.time() - t0, 1)
    with open("results.json", "w") as fh:
        json.dump(RESULTS, fh, indent=2, default=float)
    print(f"\nAll done in {RESULTS['total_runtime_s']} s -> results.json, {FIGDIR}/")
    if "all" in args:
        sys.exit(1 if verify(RESULTS) else 0)


if __name__ == "__main__":
    main()
