E-004-sine-moments.py

{
  "experiment": "E-004",
  "task": "1.3 warm-up (thread T1)",
  "description": "Monte Carlo reproduction of sine-kernel Gram matrix spectral moments m_k(1)=1,4/3,2,13/4 and Christoffel Lambda_2(0)=5/36 (paper §7.2(f))",
  "script": "research/numerics/scripts/E-004-sine-moments.py",
  "parameters": {
    "N_GUE": 400,
    "W": 600,
    "construction": "integer grid, sinc kernel (dimensionless limit object)",
    "n_samples": 60,
    "seed": 20260817
  },
  "christoffel_simulated": 0.14529137587214944,
  "checks": {
    "moments": {
      "m1_within_6pct": true,
      "m2_within_6pct": true,
      "m3_within_6pct": true,
      "m4_within_6pct": true
    },
    "christoffel": {
      "within_15pct": true
    }
  },
  "elapsed_seconds": 3.7,
  "timestamp_utc": "2026-08-17T20:56:30.918985+00:00",
  "row_count": 1
}

"""E-004 — Thread T1 / task 1.3 warm-up: reproduce the limiting spectral
moments m_k(1) = 1, 4/3, 2, 13/4 (k <= 4) of the sine-kernel Gram matrix
over the sine process (paper §7.2(f)), and the Christoffel value
Lambda_2(0) = 5/36.

Construction (flat window; "sine-kernel Gram matrix"):
  points gamma: sine-process sample (GUE bulk eigenvalues, unfolded to
  mean spacing 1) in a window of W points;
  grid alpha_k = 2 pi k / L at critical density h = 2pi/L, k = 0..d-1,
  d = floor(L W / 2 pi);
  G_tilde = V^T V with V[j,k] = sinc(L(gamma_j - alpha_k)/2),
  sinc(x) = sin(x)/x   (sharp cutoff; a = ||phi||^2/L = 1).
Moments: m_k = (1/d) tr G_tilde^k  = mean of eigenvalues^k.

Hand check (already verified, D-010): from m_0..m_4 = 1,1,4/3,2,13/4 the
moment matrix has det = 5/108 and (M^-1)_00 = 36/5, so Lambda_2(0) = 5/36.

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

import numpy as np

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

RNG_SEED = 20260817
N_GUE = 400          # GUE size per sample (bulk used)
W = 600              # window length in unfolded points
N_SAMPLES = 60
TARGET = {1: 1.0, 2: 4.0 / 3.0, 3: 2.0, 4: 13.0 / 4.0}


def gue_bulk_points(rng) -> np.ndarray:
    """Unfolded GUE bulk eigenvalues (mean spacing 1), W consecutive."""
    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)  # semicircle on [-2, 2]
    # take middle 60% (bulk) and unfold: rescale to mean spacing exactly 1
    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 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()

    # Dimensionless limit object (spacing-1 units): grid = integers,
    # kernel S(x) = sin(pi x)/(pi x) = np.sinc(x); G = sum_j S S^T.
    # (Paper's L cancels exactly; grid spacing = mean spacing = 1.)
    # Internal identity check: sum_k S(x-k)^2 = 1 exactly (sharp-window
    # Lemma 2.1) => m_1 = 1 up to window-edge effects.
    acc = {k: [] for k in (1, 2, 3, 4)}
    for s in range(N_SAMPLES):
        pts = gue_bulk_points(rng)
        pts = pts - pts[0]
        Wpts = int(pts[-1])
        alpha = np.arange(Wpts + 1)
        Z = pts[:, None] - alpha[None, :]
        V = np.sinc(Z)
        G = V.T @ V
        ev = np.linalg.eigvalsh(G)
        npts = len(pts)
        for k in (1, 2, 3, 4):
            acc[k].append(float(np.sum(ev**k) / npts))
        if (s + 1) % 20 == 0:
            print(f"  {s+1}/{N_SAMPLES} samples ({time.time()-t0:.0f}s)",
                  flush=True)

    rows = []
    checks = {}
    row = {"construction": "sinc-integer-grid (L->infinity limit object)"}
    for k in (1, 2, 3, 4):
        arr = np.array(acc[k])
        m, sd = arr.mean(), arr.std()
        row[f"m{k}"] = f"{m:.5f}"
        row[f"m{k}_std"] = f"{sd:.5f}"
        row[f"m{k}_target"] = f"{TARGET[k]:.5f}"
        row[f"m{k}_relerr"] = f"{abs(m / TARGET[k] - 1):.5f}"
    rows.append(row)
    # generous tolerances: finite-window edge effects + sampling noise
    checks["moments"] = {
        f"m{k}_within_6pct": bool(abs(np.mean(acc[k]) / TARGET[k] - 1) < 0.06)
        for k in (1, 2, 3, 4)
    }
    print("sine-kernel Gram (limit object): " + "  ".join(
        f"m{k}={np.mean(acc[k]):.4f}+-{np.std(acc[k]):.4f}"
        f" (target {TARGET[k]:.4f})" for k in (1, 2, 3, 4)), flush=True)

    # Christoffel from simulated moments: Lambda_2(0) = 1/(M^-1)_00
    m = [1.0] + [float(np.mean(acc[k])) for k in (1, 2, 3, 4)]
    M = np.array([[m[i + j] for j in range(3)] for i in range(3)])
    lam2 = 1.0 / np.linalg.inv(M)[0, 0]
    print(f"Christoffel Lambda_2(0) from simulated moments: "
          f"{lam2:.5f}  (target 5/36 = {5/36:.5f})")
    checks["christoffel"] = {"within_15pct": bool(abs(lam2 / (5 / 36) - 1) < 0.15)}

    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-004",
        "task": "1.3 warm-up (thread T1)",
        "description": "Monte Carlo reproduction of sine-kernel Gram matrix "
                       "spectral moments m_k(1)=1,4/3,2,13/4 and Christoffel "
                       "Lambda_2(0)=5/36 (paper §7.2(f))",
        "script": "research/numerics/scripts/E-004-sine-moments.py",
        "parameters": {"N_GUE": N_GUE, "W": W, "construction": "integer grid, "
                       "sinc kernel (dimensionless limit object)",
                       "n_samples": N_SAMPLES, "seed": RNG_SEED},
        "christoffel_simulated": lam2,
        "checks": checks,
        "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")
    ok = all(all(c.values()) for c in checks.values())
    print(f"wrote {CSV.name}, {JSON.name} — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE (inspect)'}")
    return 0 if ok else 1


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