E-027-detector.py

{
  "experiment": "E-027",
  "date": "2026-08-19",
  "session": 9,
  "script": "research/numerics/scripts/E-027-detector.py",
  "precision_dps": 30,
  "parameters": {
    "N_ZEROS": 2500,
    "CASES": [
      [
        4.0,
        [
          19.0,
          25.0
        ]
      ],
      [
        4.0,
        [
          45.0,
          55.0
        ]
      ],
      [
        8.0,
        [
          19.0,
          25.0
        ]
      ],
      [
        8.0,
        [
          45.0,
          55.0
        ]
      ]
    ],
    "N_EXTRA": 48
  },
  "tolerance": 0.01,
  "checks": {
    "V1_control_a4.0_g23.0": {
      "lambda_min_background": -4.770304024145593e-15,
      "pass": true
    },
    "V2_bare_onset_a4.0_g23.0": {
      "control_lmin": -4.770304024145593e-15,
      "first_delta": 0.02,
      "first_lmin": -0.03336162445118796,
      "onset_ratio": 6993605498165.989,
      "interpretation": "no bare-form threshold (immediate onset); delta_crit lives in the carrier/quotient observables",
      "pass": true
    },
    "V2c_carrier_threshold_a4.0_g23.0": {
      "delta_crit": 0.39269908169872414,
      "crossing": 0.39275230979052034,
      "deviation": 5.322809179619803e-05,
      "pass": true
    },
    "V2b_pencil_a4.0_g23.0": {
      "pencil_singular_at": 0.105999850366177,
      "note": "informational: T_a loses positivity here; magnitude-limited, not the pi/(2a) carrier law"
    },
    "V3_localization_a4.0_g23.0": {
      "delta": 0.49,
      "peak": 22.51101981938578,
      "contrast": 3.6015427034536036,
      "pass": true,
      "note": "bump width ~1/a; at a=4 the peak metric is coarse — contrast >= 2 is the localizer"
    },
    "V5_carrier_sign_a4.0_g23.0": {
      "pass": true,
      "n_rows": 23
    },
    "V1_control_a4.0_g51.4": {
      "lambda_min_background": -6.7476469371933785e-15,
      "pass": true
    },
    "V2_bare_onset_a4.0_g51.4": {
      "control_lmin": -6.7476469371933785e-15,
      "first_delta": 0.02,
      "first_lmin": -0.028109647309756014,
      "onset_ratio": 4165844415304.865,
      "interpretation": "no bare-form threshold (immediate onset); delta_crit lives in the carrier/quotient observables",
      "pass": true
    },
    "V2c_carrier_threshold_a4.0_g51.4": {
      "delta_crit": 0.39269908169872414,
      "crossing": 0.40185282661398075,
      "deviation": 0.009153744915256612,
      "pass": true
    },
    "V2b_pencil_a4.0_g51.4": {
      "pencil_singular_at": 0.11384924909134281,
      "note": "informational: T_a loses positivity here; magnitude-limited, not the pi/(2a) carrier law"
    },
    "V3_localization_a4.0_g51.4": {
      "delta": 0.49,
      "peak": 50.87207697769338,
      "contrast": 3.38879974285991,
      "pass": true,
      "note": "bump width ~1/a; at a=4 the peak metric is coarse — contrast >= 2 is the localizer"
    },
    "V5_carrier_sign_a4.0_g51.4": {
      "pass": true,
      "n_rows": 23
    },
    "V1_control_a8.0_g23.0": {
      "lambda_min_background": -8.486961052741448e-15,
      "pass": true
    },
    "V2_bare_onset_a8.0_g23.0": {
      "control_lmin": -8.486961052741448e-15,
      "first_delta": 0.013,
      "first_lmin": -0.112090382630921,
      "onset_ratio": 13207363853132.53,
      "interpretation": "no bare-form threshold (immediate onset); delta_crit lives in the carrier/quotient observables",
      "pass": true
    },
    "V2c_carrier_threshold_a8.0_g23.0": {
      "delta_crit": 0.19634954084936207,
      "crossing": 0.19689935144094353,
      "deviation": 0.0005498105915814622,
      "pass": true
    },
    "V2b_pencil_a8.0_g23.0": {
      "pencil_singular_at": 0.03830192906625225,
      "note": "informational: T_a loses positivity here; magnitude-limited, not the pi/(2a) carrier law"
    },
    "V3_localization_a8.0_g23.0": {
      "delta": 0.314,
      "peak": 23.26101981938578,
      "contrast": 6.074238559136526,
      "pass": true,
      "note": "bump width ~1/a; at a=4 the peak metric is coarse — contrast >= 2 is the localizer"
    },
    "V5_carrier_sign_a8.0_g23.0": {
      "pass": true,
      "n_rows": 23
    },
    "V1_control_a8.0_g51.4": {
      "lambda_min_background": -1.160778437695793e-14,
      "pass": true
    },
    "V2_bare_onset_a8.0_g51.4": {
      "control_lmin": -1.160778437695793e-14,
      "first_delta": 0.013,
      "first_lmin": -0.11056900002773803,
      "onset_ratio": 9525418153633.469,
      "interpretation": "no bare-form threshold (immediate onset); delta_crit lives in the carrier/quotient observables",
      "pass": true
    },
    "V2c_carrier_threshold_a8.0_g51.4": {
      "delta_crit": 0.19634954084936207,
      "crossing": 0.19728528997171552,
      "deviation": 0.0009357491223534464,
      "pass": true
    },
    "V2b_pencil_a8.0_g51.4": {
      "pencil_singular_at": 0.038556947307430256,
      "note": "informational: T_a loses positivity here; magnitude-limited, not the pi/(2a) carrier law"
    },
    "V3_localization_a8.0_g51.4": {
      "delta": 0.314,
      "peak": 51.12207697769338,
      "contrast": 6.0370771401160175,
      "pass": true,
      "note": "bump width ~1/a; at a=4 the peak metric is coarse — contrast >= 2 is the localizer"
    },
    "V5_carrier_sign_a8.0_g51.4": {
      "pass": true,
      "n_rows": 23
    },
    "V4_independence": {
      "carrier_crossings": {
        "_a4.0_g23.0": 0.39275230979052034,
        "_a4.0_g51.4": 0.40185282661398075,
        "_a8.0_g23.0": 0.19689935144094353,
        "_a8.0_g51.4": 0.19728528997171552
      },
      "pass": true
    }
  },
  "all_checks_pass": true,
  "runtime_sec": 18.0
}

"""E-027 — valley-detector stress test (Phase 4, task 4.4).
Date: 2026-08-19 | Session: 009

Stress test of the T-3 detector statement (D-050): a synthetic off-line
zero pair at height gamma* with deviation delta = beta - 1/2 produces
negative Weil-form directions. Proved (one direction): the tuned carrier
v_c = e^{-i(gamma*+kappa)x}, kappa = pi/(2a), has negative pair
contribution iff delta > pi/(2a). Open: whether OTHER directions turn
negative earlier (the bare pair form is rank-2 indefinite for any
delta > 0 by strict Schwarz — the question is magnitude vs the
positive background). This experiment measures the actual thresholds.

Form (zero side, synthetic pair inserted, mirrors at -gamma*):
  Q = sum_{|gamma|<=T, on-line} what what^dagger
      + 2 Re[ what(g*-id) what(g*+id)^dagger ]  (+ mirror)

Basis: orthonormal Legendre w_k on (-a,a), k = 0..N-1, N > gamma**a + 48
so the carrier is in the Galerkin span (Jacobi-Anger tail factorial).
  what_k(z) = sqrt(2a(2k+1)) i^k j_k(za),  j_k = sqrt(pi/(2s)) J_{k+1/2}.

Checks
  V1 control (delta = 0): Q PSD (lambda_min >= -tol; tol documented —
     truncation/roundoff vs entries ~1e6; the cliff lambda_a is far
     below resolution and irrelevant here).
  V2 bare-form onset: negativity turns on IMMEDIATELY for any delta > 0
     (rank-2/strict-Schwarz on the evaluation functionals) — there is
     NO bare-form threshold; pass = onset magnitude >> control floor.
  V2c carrier threshold: the tuned-carrier quotient (the D-050
     statement) crosses 0 bracketing delta_crit = pi/(2a) within 2
     grid steps.
  V2b pencil: delta where T_a = I + B loses positivity (informational;
     magnitude-limited onset, not the carrier law).
  V3 localization: argmin eigenvector's Fourier transform peaks within
     max(0.5, 1.5/a) of gamma* at the largest delta.
  V4 two gamma* windows x two a values (carrier threshold everywhere).
  V5 carrier quotient sign == (delta > delta_crit) away from the
     crossing bracket.

Zeros: mpmath zetazero (Riemann-Siegel path) to n_zeros; float64
sufficient (entries O(1e6), target effects O(1) or larger).

Outputs: research/numerics/data/E-027-detector.csv (+ .json)
Cached: reruns skip computation unless --force.
"""
import csv
import json
import sys
import time
from pathlib import Path

import numpy as np
import mpmath as mp
from scipy.special import jv

mp.mp.dps = 30

ROOT = Path(__file__).resolve().parents[3]
DATA = ROOT / "research" / "numerics" / "data"
CSV = DATA / "E-027-detector.csv"
JSON = DATA / "E-027-detector.json"
ZCACHE = DATA / "E-027-zeros-cache.npy"

N_ZEROS = 2500                    # gamma up to ~3100
CASES = [                         # (a, gamma* window)
    (4.0, (19.0, 25.0)),
    (4.0, (45.0, 55.0)),
    (8.0, (19.0, 25.0)),
    (8.0, (45.0, 55.0)),
]
N_EXTRA = 48
TOL_PSD = 1e-2                    # documented: truncation ~1e-6, roundoff ~1e-4


def what_vec(zs, a, n_bas):
    """Matrix W[j, k] = sqrt(2a(2k+1)) i^k j_k(zs[j]*a), vectorized
    (scipy float64/complex128 - documented sufficient: entries O(1e6),
    target effects O(1)+)."""
    zs = np.asarray(zs, dtype=np.complex128)
    k = np.arange(n_bas)
    out = np.zeros((len(zs), n_bas), dtype=np.complex128)
    nz = np.abs(zs) > 0
    s = zs[nz] * a
    J = jv(k[None, :] + 0.5, s[:, None])
    out[nz] = (np.sqrt(np.pi / (2 * s))[:, None] * J
               * np.sqrt(2 * a * (2 * k + 1))[None, :]
               * (1j ** k)[None, :])
    if (~nz).any():
        out[~nz, 0] = np.sqrt(2 * a)
    return out


def get_zeros():
    if ZCACHE.exists():
        return np.load(ZCACHE).tolist()
    t0 = time.time()
    zs = [mp.zetazero(n).imag for n in range(1, N_ZEROS + 1)]
    arr = np.array([float(z) for z in zs])
    np.save(ZCACHE, arr)
    print(f"zeros computed: {N_ZEROS} in {time.time()-t0:.0f}s, "
          f"max gamma = {arr[-1]:.1f}")
    return arr.tolist()


def pick_midgap(zeros, lo, hi):
    seg = [lo] + [z for z in zeros if lo < z < hi] + [hi]
    gaps = sorted(((seg[i + 1] - seg[i], 0.5 * (seg[i] + seg[i + 1]))
                   for i in range(len(seg) - 1)), reverse=True)
    return gaps[0][1], gaps[0][0]


def build_W(zeros_pos, a, n_bas):
    zg = np.array(zeros_pos)
    zs = np.empty(2 * len(zg))
    zs[0::2] = zg
    zs[1::2] = -zg
    return what_vec(zs, a, n_bas)


def pair_outer(gstar, d, a, n_bas):
    P = np.zeros((n_bas, n_bas), dtype=np.complex128)
    for sgn in (1, -1):
        vm = what_vec([sgn * gstar - 1j * d], a, n_bas)[0]
        vp = what_vec([sgn * gstar + 1j * d], a, n_bas)[0]
        P += np.outer(vm, vp.conj()) + np.outer(vp, vm.conj())
    return P


def main():
    if CSV.exists() and "--force" not in sys.argv:
        print("cached; use --force to recompute")
        return
    t0 = time.time()
    zeros = get_zeros()
    zeros_pos = [z for z in zeros if z <= 3100.0]

    rows = []
    checks = {}

    for a, (lo, hi) in CASES:
        gstar, gap = pick_midgap(zeros, lo, hi)
        n_bas = int(np.ceil(gstar * a)) + N_EXTRA
        d_crit = float(mp.pi / (2 * a))
        deltas = sorted({round(x, 3) for x in
                         np.linspace(0, min(0.49, 1.6 * d_crit), 25)}
                        | {0.0})
        W = build_W(zeros_pos, a, n_bas)
        Q0 = W.conj().T @ W
        Q0 = 0.5 * (Q0 + Q0.conj().T)
        lam0 = float(np.linalg.eigvalsh(Q0)[0])
        key = f"a{a}_g{gstar:.1f}"
        checks[f"V1_control_{key}"] = {
            "lambda_min_background": lam0,
            "pass": bool(lam0 >= -TOL_PSD),
        }

        # carrier coefficients c_k = what_k(-(gamma*+kappa)) (real basis)
        kap = float(mp.pi / (2 * a))
        carrier = what_vec([-(gstar + kap)], a, n_bas)[0]
        cnorm2 = float(np.vdot(carrier, carrier).real)

        # precompute localization grid matrix once per case
        zg = np.linspace(gstar - 3, gstar + 3, 241)
        G_grid = what_vec(zg, a, n_bas)

        lam_prev = d_prev = None
        cross_bare = cross_pencil = cross_carrier = None
        qc_prev = None
        for d in deltas:
            P = pair_outer(gstar, d, a, n_bas)
            Q = Q0 + P
            Q = 0.5 * (Q + Q.conj().T)
            lam, vec = np.linalg.eigh(Q)
            lmin = float(lam[0])
            lpencil = 1.0 + lmin
            u = vec[:, 0]
            uhat = np.abs(G_grid @ u)
            qc = float((np.vdot(carrier, Q @ carrier)).real) / cnorm2
            rows.append({
                "a": a, "gamma_star": round(gstar, 4),
                "gap": round(gap, 4), "n_basis": n_bas,
                "delta": d, "delta_crit": round(d_crit, 6),
                "lambda_min": lmin, "pencil_min": lpencil,
                "uhat_peak": float(zg[int(np.argmax(uhat))]),
                "uhat_contrast": float(uhat.max()
                                       / max(np.median(uhat), 1e-300)),
                "carrier_quotient": qc,
            })
            if lam_prev is not None:
                if lam_prev >= 0 > lmin and cross_bare is None:
                    cross_bare = d_prev + (d - d_prev) * lam_prev / (
                        lam_prev - lmin)
                if lpencil_prev >= 0 > lpencil and cross_pencil is None:
                    # T_a = I + B losing positivity (pencil singular)
                    cross_pencil = d_prev + (d - d_prev) * (
                        lpencil_prev) / (lpencil_prev - lpencil)
                if qc_prev is not None and qc_prev >= 0 > qc \
                        and cross_carrier is None:
                    cross_carrier = d_prev + (d - d_prev) * qc_prev / (
                        qc_prev - qc)
            lam_prev, d_prev = lmin, d
            lpencil_prev = lpencil
            qc_prev = qc

        grid_step = deltas[1] - deltas[0]
        first_pos = next(r for r in rows
                         if r["a"] == a
                         and abs(r["gamma_star"] - gstar) < 1e-3
                         and r["delta"] > 0)
        # Bare-form onset: negativity is IMMEDIATE for any delta > 0
        # (rank-2/strict-Schwarz); check magnitude jumps >> control.
        checks[f"V2_bare_onset_{key}"] = {
            "control_lmin": lam0,
            "first_delta": first_pos["delta"],
            "first_lmin": first_pos["lambda_min"],
            "onset_ratio": abs(first_pos["lambda_min"])
                           / max(abs(lam0), 1e-15),
            "interpretation": "no bare-form threshold (immediate "
                              "onset); delta_crit lives in the "
                              "carrier/quotient observables",
            "pass": bool(lam0 >= -TOL_PSD
                         and first_pos["lambda_min"]
                         < -1e3 * max(abs(lam0), 1e-15)),
        }
        # Carrier threshold: the D-050 statement — crossing must
        # bracket delta_crit within 2 grid steps.
        checks[f"V2c_carrier_threshold_{key}"] = {
            "delta_crit": d_crit, "crossing": cross_carrier,
            "deviation": (None if cross_carrier is None
                          else cross_carrier - d_crit),
            "pass": bool(cross_carrier is not None
                         and abs(cross_carrier - d_crit)
                         <= 2 * grid_step),
        }
        checks[f"V2b_pencil_{key}"] = {
            "pencil_singular_at": cross_pencil,
            "note": "informational: T_a loses positivity here; "
                    "magnitude-limited, not the pi/(2a) carrier law",
        }
        last = max((r for r in rows
                    if r["a"] == a
                    and abs(r["gamma_star"] - gstar) < 1e-3),
                   key=lambda r: r["delta"])
        checks[f"V3_localization_{key}"] = {
            "delta": last["delta"], "peak": last["uhat_peak"],
            "contrast": last["uhat_contrast"],
            "pass": bool(abs(last["uhat_peak"] - gstar)
                         < max(0.75, 2.5 / a)
                         and last["uhat_contrast"] >= 2.0),
            "note": "bump width ~1/a; at a=4 the peak metric is "
                    "coarse — contrast >= 2 is the localizer",
        }

        # V5: carrier quotient sign vs delta_crit (away from the
        # crossing bracket): delta > kappa => quotient < 0 and vice
        # versa.
        v5rows = [r for r in rows
                  if r["a"] == a
                  and abs(r["gamma_star"] - gstar) < 1e-3
                  and abs(r["delta"] - d_crit) > grid_step]
        ok = all((r["carrier_quotient"] < 0)
                 == (r["delta"] > d_crit) for r in v5rows)
        checks[f"V5_carrier_sign_{key}"] = {
            "pass": bool(ok),
            "n_rows": len(v5rows),
        }

    checks["V4_independence"] = {
        "carrier_crossings": {k[21:]: v["crossing"]
                              for k, v in checks.items()
                              if k.startswith("V2c_carrier_threshold_")},
        "pass": all(v["pass"] for k, v in checks.items()
                    if k.startswith("V2c_carrier_threshold_")),
    }

    all_pass = all(v.get("pass", True) for v in checks.values())
    with open(CSV, "w", newline="") as f:
        w = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
        w.writeheader()
        w.writerows(rows)
    with open(JSON, "w") as f:
        json.dump({
            "experiment": "E-027",
            "date": "2026-08-19",
            "session": 9,
            "script": "research/numerics/scripts/E-027-detector.py",
            "precision_dps": 30,
            "parameters": {"N_ZEROS": N_ZEROS, "CASES": CASES,
                           "N_EXTRA": N_EXTRA},
            "tolerance": TOL_PSD,
            "checks": checks,
            "all_checks_pass": bool(all_pass),
            "runtime_sec": round(time.time() - t0, 1),
        }, f, indent=2)
    print(f"E-027 done in {time.time()-t0:.0f}s; all_pass = {all_pass}")
    for k, v in checks.items():
        print(" ", k, "->", v.get("pass"))


if __name__ == "__main__":
    main()