E-003-finite-Gtilde.py

{
  "experiment": "E-003",
  "task": "0.6b",
  "description": "Finite compression G̃ on [T,2T]=[2000,4000): trace, HS norm, PSD-ness, rank–trace bound vs paper predictions (Props 4.1–4.3, Thm 5.7)",
  "script": "research/numerics/scripts/E-003-finite-Gtilde.py",
  "parameters": {
    "T_values": [
      1000,
      2000,
      4000
    ],
    "windows": [
      "psi0",
      "psiMT"
    ],
    "zeros_source": "E-002 (dps 25)",
    "precision": "float64"
  },
  "checks": {
    "psi0@T=1000": {
      "trace_within_1pct": true,
      "hs2_within_expected_finite_L": true,
      "frame_selftest_ok": true,
      "numerically_psd": true
    },
    "psiMT@T=1000": {
      "trace_within_1pct": true,
      "hs2_within_expected_finite_L": true,
      "frame_selftest_ok": true,
      "numerically_psd": true
    },
    "psi0@T=2000": {
      "trace_within_1pct": true,
      "hs2_within_expected_finite_L": true,
      "frame_selftest_ok": true,
      "numerically_psd": true
    },
    "psiMT@T=2000": {
      "trace_within_1pct": true,
      "hs2_within_expected_finite_L": true,
      "frame_selftest_ok": true,
      "numerically_psd": true
    },
    "psi0@T=4000": {
      "trace_within_1pct": true,
      "hs2_within_expected_finite_L": true,
      "frame_selftest_ok": true,
      "numerically_psd": true
    },
    "psiMT@T=4000": {
      "trace_within_1pct": true,
      "hs2_within_expected_finite_L": true,
      "frame_selftest_ok": true,
      "numerically_psd": true
    }
  },
  "elapsed_seconds": 6.0,
  "timestamp_utc": "2026-08-17T20:14:13.193268+00:00",
  "row_count": 6
}

"""E-003 — Task 0.6b: finite compression G̃ — empirical check of the paper's
zero-side/trace/HS-norm predictions against the first 10^4 zeros (E-002).

Construction follows D-002 §3 exactly (paper §2):
  L := log(T/2π), X = T/2π, grid α_k = T + 2πk/L, d = ⌊LT/2π⌋
  ϕ(u) = χ(L/2+u)χ(L/2−u)·√ψ(u/L),   supp ϕ = [−L/2, L/2]
  a := ‖ϕ‖²₂/L,  v_ρ = (ϕ̂(γ_ρ − α_k))_{0≤k<d}
  G̃ := (aL²)⁻¹ Σ_{Re γ_ρ ∈ I′} m_ρ v_ρ v_ρᵀ,  I′ = [T−√T, 2T+√T)

Predictions checked (all zeros here are on the line and simple, so
P = G̃ up to the tail Ẽ, Q ≈ 0):
  (i)   tr G̃ ≈ N(I)                       (Prop 4.2)
  (ii)  ‖G̃‖²_HS ≈ R(ψ)·N(I)               (Thm 5.7; R from E-001)
  (iii) G̃ numerically PSD; rank ≈ N(I)    (Prop 4.1; n₊(Q) = 0)
  (iv)  rank–trace lower bound 4tr G̃ − 2N − ‖G̃‖² ≈ (2−R)N ≤ s₁
Windows: ψ₀ (indicator) and ψ_MT.

NOTE on scale: at T = 2000 the asymptotic error terms (O(√T L²) on the
trace, O(1/L) relative on the HS norm) are only moderately small; we expect
agreement at the few-percent level, not machine precision. This is a
structural check, not a precision verification.

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

import numpy as np
from scipy.interpolate import CubicSpline
from scipy.linalg import eigvalsh

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
ZEROS_CSV = DATA / "E-002-zeros-first-10000.csv"
CSV = DATA / "E-003-finite-Gtilde.csv"
JSON = DATA / "E-003-finite-Gtilde.json"

T = 2000.0
R_PSI0 = 4.0 / 3.0
R_PSIMT = 0.5 + (1.0 / np.sqrt(2.0)) / np.tan(1.0 / np.sqrt(2.0))  # cot


def chi(u: np.ndarray) -> np.ndarray:
    """Smooth 0→1 transition on [0,1] (C^∞ exp mollifier)."""
    out = np.zeros_like(u)
    mid = (u > 0) & (u < 1)
    e1 = np.exp(-1.0 / u[mid])
    e2 = np.exp(-1.0 / (1.0 - u[mid]))
    out[mid] = e1 / (e1 + e2)
    out[u >= 1] = 1.0
    return out


def build_phi_hat(L: float, psi, pad: int = 8, M: int = 4096,
                  xi_max_needed: float = 0.0):
    """Return callable ϕ̂(ξ) via zero-padded FFT + cubic spline.

    ϕ sampled on [−L/2, L/2] with M points; zero-padded to pad·M so the
    DFT frequency grid has spacing 2π/(pad·L). The spline domain is
    |ξ| ≤ πM/L — M is raised to the next power of two so the domain
    covers xi_max_needed (no extrapolation allowed).
    """
    M_need = int(2 ** np.ceil(np.log2(max(M, L * xi_max_needed / np.pi))))
    if M_need != M:
        M = M_need
    du = L / M
    u = -L / 2 + du * np.arange(M)
    phi = chi(L / 2 + u) * chi(L / 2 - u) * np.sqrt(psi(u / L))
    a = np.sum(phi**2) * du / L  # ‖ϕ‖²₂/L

    # Zero-pad to pad*M samples, rotating so index 0 ↔ u=0 (DFT convention);
    # then FFT[m]·du = ϕ̂(ξ_m) directly on the signed grid fftfreq.
    full = np.zeros(pad * M)
    full[: M // 2] = phi[M // 2:]          # u in [0, L/2)
    full[-M // 2:] = phi[: M // 2]         # u in [−L/2, 0)
    F = np.fft.fft(full) * du
    xi_signed = 2 * np.pi * np.fft.fftfreq(pad * M, d=du)
    order = np.argsort(xi_signed)
    spline = CubicSpline(xi_signed[order], F[order])
    return spline, a


def main() -> int:
    if not ZEROS_CSV.exists():
        print("E-002 zeros not found — run E-002-first-zeros.py first")
        return 1
    if CSV.exists() and "--force" not in sys.argv:
        print(f"cached: {CSV} exists (use --force to recompute)")
        return 0

    t0 = time.time()
    gam = np.loadtxt(ZEROS_CSV, delimiter=",", skiprows=1, usecols=1)

    rows = []
    checks = {}
    for T in (1000.0, 2000.0, 4000.0):
        L = np.log(T / (2 * np.pi))
        h = 2 * np.pi / L
        d = int(L * T / (2 * np.pi))
        alpha = T + h * np.arange(d)
        sqT = np.sqrt(T)
        in_I = (gam >= T) & (gam < 2 * T)
        in_Ip = (gam >= T - sqT) & (gam < 2 * T + sqT)
        N_I = int(in_I.sum())
        zeros_Ip = gam[in_Ip]
        print(f"\nT={T:.0f}: L={L:.4f}  d={d}  N(I)={N_I}  "
              f"zeros in I'={len(zeros_Ip)}", flush=True)

        psi_fns = {
            "psi0": lambda s: np.ones_like(s),
            "psiMT": lambda s: np.cos(np.sqrt(2.0) * s),
        }
        R_expected = {"psi0": R_PSI0, "psiMT": R_PSIMT}

        for name, psi in psi_fns.items():
            t1 = time.time()
            phi_hat, a = build_phi_hat(L, psi,
                                       xi_max_needed=2 * T + sqT)

            # Lemma 2.1 self-test: sum_k phi_hat(x - alpha_k)^2 == a*L**2
            # (validates the FFT+spline phi_hat pipeline independent of G)
            k_ext = np.arange(-2 * d, 3 * d)
            alpha_ext = T + h * k_ext
            x_tests = np.array([T + 0.3 * T, T + 0.5 * T, T + 0.8 * T])
            frame_err = 0.0
            for x in x_tests:
                s = np.sum(np.real(phi_hat(x - alpha_ext))**2)
                frame_err = max(frame_err, abs(s / (a * L**2) - 1))
            print(f"  [{name}] Lemma 2.1 frame self-test max rel err: "
                  f"{frame_err:.2e}", flush=True)

            Z = zeros_Ip[:, None] - alpha[None, :]
            V = np.real(phi_hat(Z))
            G = (V.T @ V) / (a * L**2)          # all m_rho = 1
            trG = np.trace(G)
            hs2 = np.sum(G**2)
            ev = eigvalsh(G)
            tol = 1e-8 * ev[-1]
            rank_num = int(np.sum(ev > tol))
            min_ev = float(ev[0])
            Nf = float(N_I)
            lb = 4 * trG - 2 * Nf - hs2
            rows.append({
                "T": int(T),
                "window": name,
                "L": f"{L:.4f}",
                "N_I": N_I,
                "tr_G_over_N": f"{trG / Nf:.5f}",
                "HS2_over_N": f"{hs2 / Nf:.5f}",
                "R_psi_expected": f"{R_expected[name]:.5f}",
                "HS2_rel_err": f"{abs(hs2 / Nf / R_expected[name] - 1):.5f}",
                "expected_O(1/L)": f"{1 / L:.4f}",
                "frame_selftest_rel_err": f"{frame_err:.2e}",
                "rank_numeric": rank_num,
                "d": d,
                "min_eigenvalue": f"{min_ev:.3e}",
                "rank_trace_lower_bound_over_N": f"{lb / Nf:.5f}",
                "two_minus_R": f"{2 - R_expected[name]:.5f}",
            })
            key = f"{name}@T={int(T)}"
            checks[key] = {
                "trace_within_1pct": bool(abs(trG / Nf - 1) < 0.01),
                # theory error is O(1/L) ~ 17% at T=2000; allow 3.5/L
                "hs2_within_expected_finite_L": bool(
                    abs(hs2 / Nf / R_expected[name] - 1) < 3.5 / L),
                "frame_selftest_ok": bool(frame_err < 5e-5),
                "numerically_psd": bool(min_ev > -1e-6 * float(ev[-1])),
            }
            print(f"  [{name}] tr/N={trG/Nf:.4f}  HS2/N={hs2/Nf:.4f} "
                  f"(R={R_expected[name]:.4f}, rel err "
                  f"{abs(hs2/Nf/R_expected[name]-1):.4f} vs 1/L={1/L:.4f})  "
                  f"rank={rank_num}/{d}  min_ev={min_ev:.2e}  "
                  f"(4tr-2N-HS2)/N={lb/Nf:.4f} "
                  f"({time.time()-t1:.1f}s)", flush=True)

    CSV.parent.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-003",
        "task": "0.6b",
        "description": "Finite compression G̃ on [T,2T]=[2000,4000): trace, "
                       "HS norm, PSD-ness, rank–trace bound vs paper "
                       "predictions (Props 4.1–4.3, Thm 5.7)",
        "script": "research/numerics/scripts/E-003-finite-Gtilde.py",
        "parameters": {"T_values": [1000, 2000, 4000],
                       "windows": ["psi0", "psiMT"],
                       "zeros_source": "E-002 (dps 25)",
                       "precision": "float64"},
        "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 tolerances)'}")
    return 0 if ok else 1


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