E-020-xiprime-ceiling.py

{
  "experiment": "E-020",
  "task": "1.6: xi'-world bandwidth-one ceiling (LawN256 analogue)",
  "script": "research/numerics/scripts/E-020-xiprime-ceiling.py",
  "parameters": {
    "N": 256,
    "n_cue": 500,
    "pool": 2304,
    "kappas": [
      1.0,
      1.5,
      3.0
    ]
  },
  "checks": {
    "C1_control_validates": false,
    "C2_D1_law_sane": true,
    "C3_LP_optimal": true
  },
  "report": {
    "xi_control": {
      "eps_star": 0.0,
      "stage1_status": 0,
      "ceiling_k1.0": 0.9958868145666365,
      "status_k1.0": 0,
      "ceiling_k1.5": 0.9958868145666365,
      "status_k1.5": 0,
      "ceiling_k3.0": 0.9958868145666365,
      "status_k3.0": 0
    },
    "xiprime": {
      "eps_star": 155.6324559280527,
      "stage1_status": 0,
      "ceiling_k1.0": 0.8312190947136757,
      "status_k1.0": 0,
      "ceiling_k1.5": 0.47884448965860277,
      "status_k1.5": 0,
      "ceiling_k3.0": 0.39148654100165314,
      "status_k3.0": 0
    }
  },
  "caveats": [
    "finite pool => ceiling estimates are UPPER bounds on the true LP ceiling (fewer adversary configs)",
    "eps* is the pool's best uniform row matching; the paper's exact-rational continuous-position LP reaches 1e-40",
    "ceiling(kappa) at tol = kappa*eps*; kappa=1 is the pool-reach value"
  ],
  "elapsed_seconds": 78.6,
  "timestamp_utc": "2026-08-18T20:18:07.930643+00:00"
}

"""E-020 — xi'-world bandwidth-one ceiling (task 1.6): LawN256 analogue
for kernel D1. TWO-STAGE LP version (v2).

Mirror of the paper's Sec.7.2 ceiling LP (D-009): the adversary mixes
256-point configurations (marks {1,2}, mass 256); row quantity is the
grid form factor |F_c(j)|^2, F_c(j) = sum_i m_i e^{2 pi i j x_i/256}.
  Stage 1: eps* = min ||A w - law||_inf   (best row matching the pool
                                          can do; LP with slack)
  Stage 2: ceiling = min p1.w  s.t. ||A w - law||_inf <= kappa*eps*,
           kappa in {1, 1.5, 3}
Laws: xi CONTROL law(j) = j (paper's ceiling 0.6818 validates the
machinery); xi' law(j) = 256*D1(j/256) (E-013/D-022 kernel).

Pool: CUE-256 eigenangle samples (E|F|^2 = j in expectation — the GUE
datum) + adversarial corruptions (doubles eating positions), lattice &
lattice+doubles (zero-row "dark matter"), random configs. Finite pool =>
UPPER bound on the true ceiling; eps* reports the pool's row reach.

Checks:
  C1. control: eps* <= 1.0 and ceiling(1x) in [0.66, 0.85].
  C2. D1 law sanity: D1 >= 0 on [0,1], node ratio < 0.1.
  C3. all stage-1/2 LPs solve to optimality.

Outputs: research/numerics/data/E-020-xiprime-ceiling.csv (+ .json)
"""
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import numpy as np
import scipy.optimize
import scipy.stats

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-020-xiprime-ceiling.csv"
JSON = DATA / "E-020-xiprime-ceiling.json"

N = 256
N_CUE = 500
DBL_FRAC = [(4, 500), (16, 500), (64, 300)]
K_TAIL = 40
KAPPAS = [1.0, 1.5, 3.0]
rng = np.random.default_rng(20)


def d1(x):
    from math import factorial
    x = np.asarray(x, dtype=float)
    s = x - 4 * x * x
    xp = x ** 3
    for k in range(K_TAIL):
        ck = 2 * 4 ** (k + 1) * factorial(k) / factorial(2 * k + 2)
        s = s + ck * xp
        xp = xp * x * x
    return s


def cue_positions():
    """256 eigenangles of a Haar-CUE matrix, mapped to [0,256)."""
    U = scipy.stats.unitary_group.rvs(N, random_state=rng)
    ev = np.linalg.eigvals(U)
    ang = np.angle(ev) % (2 * np.pi)
    return ang * (N / (2 * np.pi))


def build_pool():
    """Yields (positions array, marks array)."""
    cue_sets = []
    for _ in range(N_CUE):
        pos = cue_positions()
        cue_sets.append(pos)
        yield pos, np.ones(N)
        for n_d, cap in DBL_FRAC:
            if len(cue_sets) > cap:
                continue
            idx = rng.choice(N, size=2 * n_d, replace=False)
            keep = np.ones(N, dtype=bool)
            keep[idx[n_d:]] = False            # eaten positions
            marks = np.ones(N)
            marks[idx[:n_d]] = 2.0
            yield pos[keep], marks[keep]
    # lattice & lattice+doubles (zero-row dark matter + near-dark)
    lat = np.arange(N, dtype=float)
    yield lat, np.ones(N)
    for n_d in (8, 32, 64):
        marks = np.ones(N)
        idx = rng.choice(N, size=2 * n_d, replace=False)
        keep = np.ones(N, dtype=bool)
        keep[idx[n_d:]] = False
        marks[idx[:n_d]] = 2.0
        yield lat[keep], marks[keep]
    # random scatter configs
    for _ in range(500):
        n_d = int(rng.integers(4, 96))
        n_s = N - 2 * n_d
        pos = np.sort(rng.uniform(0, N, n_s + n_d))
        marks = np.concatenate([np.ones(n_s), np.full(n_d, 2.0)])
        yield pos, marks


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

    xs = np.linspace(0, 1, 1001)
    node_ratio = float(d1(np.array([0.5]))[0] / 0.5)
    c2 = bool(np.all(d1(xs) >= -1e-12) and node_ratio < 0.1)
    print(f"C2 D1 sanity: node ratio {node_ratio:.4f}", flush=True)

    print("building pool (CUE samples + corruptions) ...", flush=True)
    j_rows = np.arange(1, N)
    A_cols, p1 = [], []
    for pos, marks in build_pool():
        E = np.exp(2j * np.pi * np.outer(j_rows, pos) / N)
        F = E @ marks
        A_cols.append(np.abs(F) ** 2)
        p1.append(float(np.sum(marks == 1.0)) / N)
    A = np.array(A_cols).T              # 255 x pool
    p1 = np.array(p1)
    pool = len(p1)
    print(f"  pool = {pool}", flush=True)

    laws = {"xi_control": j_rows.astype(float),
            "xiprime": N * d1(j_rows / N)}

    rows_out, report = [], {}
    for lname, law in laws.items():
        # stage 1: min eps
        c = np.concatenate([np.zeros(pool), [1.0]])
        A_ub = np.hstack([np.vstack([A, -A]),
                          np.full((2 * A.shape[0], 1), -1.0)])
        b_ub = np.concatenate([law, -law])
        A_eq = np.concatenate([np.ones((1, pool)), np.zeros((1, 1))], 1)
        res1 = scipy.optimize.linprog(c, A_ub=A_ub, b_ub=b_ub,
                                      A_eq=A_eq, b_eq=[1.0],
                                      bounds=(0.0, None), method="highs")
        eps_star = float(res1.x[-1]) if res1.status == 0 else None
        report[lname] = {"eps_star": eps_star,
                         "stage1_status": int(res1.status)}
        print(f"  {lname}: eps* = {eps_star}", flush=True)
        if eps_star is None:
            continue
        for kap in KAPPAS:
            tol = kap * eps_star
            A_ub2 = np.vstack([A, -A])
            b_ub2 = np.concatenate([law + tol, -(law - tol)])
            res2 = scipy.optimize.linprog(
                p1, A_ub=A_ub2, b_ub=b_ub2,
                A_eq=np.ones((1, pool)), b_eq=[1.0],
                bounds=(0.0, None), method="highs")
            ceil = float(res2.fun) if res2.status == 0 else None
            report[lname][f"ceiling_k{kap}"] = ceil
            report[lname][f"status_k{kap}"] = int(res2.status)
            rows_out.append(f"{lname},{eps_star:.6e},{kap},{tol:.6e},"
                            f"{res2.status},{ceil if ceil else ''}")
            print(f"    kappa={kap}: ceiling={ceil}", flush=True)

    ctrl_eps = report["xi_control"]["eps_star"]
    ctrl_ceil = report["xi_control"].get("ceiling_k1.0")
    c1 = bool(ctrl_eps is not None and ctrl_eps <= 1.0
              and ctrl_ceil is not None and 0.66 <= ctrl_ceil <= 0.85)
    c3 = all(v == 0 for l in report.values() for k, v in l.items()
             if k.startswith("status"))
    checks = {"C1_control_validates": c1, "C2_D1_law_sane": c2,
              "C3_LP_optimal": c3}
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("law,eps_star,kappa,tol,lp_status,ceiling_estimate\n")
        f.write("\n".join(rows_out) + "\n")
    meta = {
        "experiment": "E-020",
        "task": "1.6: xi'-world bandwidth-one ceiling (LawN256 analogue)",
        "script": "research/numerics/scripts/E-020-xiprime-ceiling.py",
        "parameters": {"N": N, "n_cue": N_CUE, "pool": pool,
                       "kappas": KAPPAS},
        "checks": checks, "report": report,
        "caveats": [
            "finite pool => ceiling estimates are UPPER bounds on the "
            "true LP ceiling (fewer adversary configs)",
            "eps* is the pool's best uniform row matching; the paper's "
            "exact-rational continuous-position LP reaches 1e-40",
            "ceiling(kappa) at tol = kappa*eps*; kappa=1 is the "
            "pool-reach value",
        ],
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    ok = all(checks.values())
    print(f"wrote {CSV.name}, {JSON.name} — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE: ' + str(checks)}")
    return 0 if ok else 1


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