E-005-oversampled-third-moment.py

{
  "experiment": "E-005",
  "task": "1.3 decider (D-011 oversampling hatch)",
  "description": "3-moment one-sided bounds for (lambda,q) sinc-Gram matrices over the sine process; decides whether an unconditional third moment (lambda <= 2/3) can beat 2/3",
  "script": "research/numerics/scripts/E-005-oversampled-third-moment.py",
  "parameters": {
    "N_GUE": 400,
    "W": 500,
    "lambdas": [
      0.5,
      0.6666666666666666,
      0.8,
      1.0
    ],
    "qs": [
      1,
      2
    ],
    "n_samples": 30,
    "seed": 20260817
  },
  "decision": "close",
  "caveat": "sine-process model cannot see the D-H/Epstein obstruction; heuristic 8/9 conversion calibrated on paper's lambda=1 two-moment values",
  "elapsed_seconds": 1.8,
  "timestamp_utc": "2026-08-17T21:17:53.219249+00:00",
  "row_count": 8
}

"""E-005 — Task 1.3 decider: can a third moment at reduced bandwidth beat
2/3? (Oversampling escape hatch, D-011.)

Model: sine process (GUE bulk, unfolded spacing 1), window of W points.
The (lambda, q) sinc-Gram matrix (dimensionless, D-011 §4):
  beta_j = lambda * gamma_j          (bandwidth lambda <-> prime support)
  grid alpha_k = k/q, k = 0..q*lambda*W   (q-oversampled, EXACT no-aliasing:
  sum_k sinc^2(pi(x - k/q)) = q for integer q >= 1 — triangle FT vanishes)
  V[j,k] = sinc(pi (beta_j - alpha_k))  (np.sinc(x) = sin(pi x)/(pi x))
  G = (1/q) V^T V                     (=> tr G = N exactly)
  d = q*lambda*W rows; dimension cap on any certificate: d/N = q*lambda.

Moments per-d: mu_r = tr G^r / d (r = 1..4; trace powers without
eigendecomposition: tr G^2 = |G|_F^2, tr G^3 = sum(G * (G@G)),
tr G^4 = |G@G|_F^2).

One-sided bounds for n_+(G)/d (the certificate-relevant quantity):
  CMS2  = 1 - Lambda_1(0)  from (mu0..mu2)   [m=1 Christoffel]
  MARK3 = optimal degree-3 polynomial minorant of 1_{x>0} given
          (mu0..mu3) — discretized LP (scipy.linprog)
  CMS4  = 1 - Lambda_2(0)  from (mu0..mu4)   [reference; needs k=4]
Baselines: paper 2-moment chain at bandwidth lambda (flat window):
  C2(1) = 2/3; general R_lambda(psi0) = 1/lambda + lambda/3 (Remark 6.1);
  cap = q*lambda.
Calibrated conversion (heuristic): paper's sharp n_+ 2-moment bound 3/4
maps to certificate 2/3 -> factor 8/9; report MARK3*(8/9) as the
heuristic certificate estimate, clearly tagged.

Decision rule (recorded in JSON):
  pursue  iff exists (lambda <= 2/3, q) with MARK3 > CMS2 + 0.01 AND
          q*lambda > 2/3 AND heuristic cert > 2/3;
  else G-002 (third-moment route closed with exact obstruction).

Outputs: research/numerics/data/E-005-oversampled-third-moment.csv (+json)
"""
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import numpy as np
from scipy.optimize import linprog

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-005-oversampled-third-moment.csv"
JSON = DATA / "E-005-oversampled-third-moment.json"

RNG_SEED = 20260817
N_GUE = 400
W = 500
N_SAMPLES = 30
LAMBDAS = (0.5, 2.0 / 3.0, 0.8, 1.0)
QS = (1, 2)
XGRID = np.linspace(0.0, 6.0, 601)  # minorant enforcement grid


def gue_bulk_points(rng) -> np.ndarray:
    A = rng.standard_normal((N_GUE, N_GUE)) + 1j * rng.standard_normal((N_GUE, N_GUE))
    H = (A + A.conj().T) / 2.0
    ev = np.linalg.eigvalsh(H) / np.sqrt(N_GUE)
    lo, hi = int(0.2 * N_GUE), int(0.8 * N_GUE)
    bulk = ev[lo:hi]
    bulk = bulk / np.mean(np.diff(bulk))
    mid = len(bulk) // 2
    return bulk[mid - W // 2: mid + W // 2]


def christoffel_1(mus):
    """Lambda_1(0) = 1/(M^-1)_00, M = [[1, mu1],[mu1, mu2]]."""
    M = np.array([[1.0, mus[1]], [mus[1], mus[2]]])
    return 1.0 / np.linalg.inv(M)[0, 0]


def christoffel_2(mus):
    M = np.array([[mus[i + j] for j in range(3)] for i in range(3)])
    return 1.0 / np.linalg.inv(M)[0, 0]


def markov3(mus):
    """Maximize sum c_r mu_r over deg-3 polys p <= 1_{x>0} on XGRID."""
    # linprog minimizes; negate objective
    c = -np.array([mus[0], mus[1], mus[2], mus[3]])
    # constraints: for x <= 0: p(x) <= 0 ; for x > 0: p(x) <= 1
    Arows, bvec = [], []
    for x in XGRID:
        powers = np.array([1.0, x, x**2, x**3])
        if x <= 0:
            Arows.append(powers)
            bvec.append(0.0)
        else:
            Arows.append(powers)
            bvec.append(1.0)
    res = linprog(c, A_ub=np.array(Arows), b_ub=np.array(bvec),
                  bounds=[(None, None)] * 4, method="highs")
    if not res.success:
        return float("nan")
    return float(-res.fun)


def main() -> int:
    if CSV.exists() and "--force" not in sys.argv:
        print(f"cached: {CSV} exists (use --force to recompute)")
        return 0
    rng = np.random.default_rng(RNG_SEED)
    t0 = time.time()

    acc = {(lam, q): {r: [] for r in (1, 2, 3, 4)}
           for lam in LAMBDAS for q in QS}
    for s in range(N_SAMPLES):
        pts = gue_bulk_points(rng)
        pts = pts - pts[0]
        npts = len(pts)
        for lam in LAMBDAS:
            beta = lam * pts
            for q in QS:
                d = int(round(q * lam * pts[-1]))
                alpha = np.arange(d + 1) / q
                V = np.sinc(beta[:, None] - alpha[None, :])
                G = (V.T @ V) / q
                G2 = G @ G
                mus_here = {
                    1: np.trace(G) / d,
                    2: np.sum(G * G) / d,
                    3: np.sum(G * G2) / d,
                    4: np.sum(G2 * G2) / d,
                }
                for r in (1, 2, 3, 4):
                    acc[(lam, q)][r].append(mus_here[r])
        if (s + 1) % 10 == 0:
            print(f"  {s+1}/{N_SAMPLES} samples ({time.time()-t0:.0f}s)",
                  flush=True)

    rows = []
    decisions = []
    for lam in LAMBDAS:
        for q in QS:
            mus = [1.0] + [float(np.mean(acc[(lam, q)][r]))
                           for r in (1, 2, 3, 4)]
            cms2 = 1.0 - christoffel_1(mus)
            cms4 = 1.0 - christoffel_2(mus)
            mk3 = markov3(mus)
            cap = q * lam
            base2 = 2.0 - (1.0 / lam + lam / 3.0)  # 2-moment baseline
            heur = mk3 * (8.0 / 9.0)
            rows.append({
                "lambda": f"{lam:.4f}", "q": q, "d_over_N": f"{cap:.4f}",
                "mu1": f"{mus[1]:.4f}", "mu2": f"{mus[2]:.4f}",
                "mu3": f"{mus[3]:.4f}", "mu4": f"{mus[4]:.4f}",
                "CMS2": f"{cms2:.4f}", "MARKOV3": f"{mk3:.4f}",
                "CMS4": f"{cms4:.4f}",
                "baseline_2moment": f"{base2:.4f}",
                "cap_q*lambda": f"{cap:.4f}",
                "heuristic_cert_MARK3*8/9": f"{heur:.4f}",
            })
            print(f"lam={lam:.3f} q={q}: mu=({mus[1]:.3f},{mus[2]:.3f},"
                  f"{mus[3]:.3f},{mus[4]:.3f})  CMS2={cms2:.4f}  "
                  f"MARK3={mk3:.4f}  CMS4={cms4:.4f}  cap={cap:.3f}  "
                  f"heur={heur:.4f}", flush=True)
            if lam <= 2.0 / 3.0 + 1e-9:
                decisions.append(
                    (mk3 > cms2 + 0.01) and (cap > 2.0 / 3.0)
                    and (heur > 2.0 / 3.0))

    pursue = any(decisions)
    print(f"\nDECISION: {'PURSUE — a config beats 2/3' if pursue else 'CLOSE — no config beats 2/3 (G-002)'}")

    DATA.mkdir(parents=True, exist_ok=True)
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write(",".join(rows[0].keys()) + "\n")
        for r in rows:
            f.write(",".join(str(v) for v in r.values()) + "\n")
    meta = {
        "experiment": "E-005",
        "task": "1.3 decider (D-011 oversampling hatch)",
        "description": "3-moment one-sided bounds for (lambda,q) sinc-Gram "
                       "matrices over the sine process; decides whether an "
                       "unconditional third moment (lambda <= 2/3) can beat 2/3",
        "script": "research/numerics/scripts/E-005-oversampled-third-moment.py",
        "parameters": {"N_GUE": N_GUE, "W": W, "lambdas": list(LAMBDAS),
                       "qs": list(QS), "n_samples": N_SAMPLES,
                       "seed": RNG_SEED},
        "decision": "pursue" if pursue else "close",
        "caveat": "sine-process model cannot see the D-H/Epstein "
                  "obstruction; heuristic 8/9 conversion calibrated on "
                  "paper's lambda=1 two-moment values",
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
        "row_count": len(rows),
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    print(f"wrote {CSV.name}, {JSON.name}")
    return 0


if __name__ == "__main__":
    raise SystemExit(main())