E-006-positivity-horizon.py

{
  "experiment": "E-006",
  "task": "3.1 / Q-3.4b (positivity horizon)",
  "description": "Eigenvalue floor of Gabor-truncated Weil form vs bandwidth lambda; conditioning of positivity certification (cf. Groskin B_T budget, Kim et al. superexponential decay)",
  "script": "research/numerics/scripts/E-006-positivity-horizon.py",
  "parameters": {
    "T": 2000.0,
    "lambdas": [
      0.5,
      1.0,
      1.5,
      2.0
    ],
    "window": "psi0",
    "zeros_source": "E-002",
    "precision": "float64"
  },
  "caveats": [
    "all zeros on-line => PSD by construction; the floor measures conditioning, not failure",
    "d > N' => dimensional rank deficiency (overcomplete basis), not form degeneracy"
  ],
  "elapsed_seconds": 2.9,
  "timestamp_utc": "2026-08-17T22:20:46.622267+00:00",
  "row_count": 4
}

"""E-006 — Task 3.1 / Q-3.4b: positivity horizon of the truncated Weil form.

For bandwidth lambda in {0.5, 1, 1.5, 2} build the Gabor-truncated Weil
form G_lambda on [T, 2T), T = 2000 (construction of E-003/D-002 §3 with
L_lambda = lambda * log(T/2pi), grid at critical density h = 2pi/L_lambda):

  phi_lambda(u) = chi(L/2+u) chi(L/2-u) * 1   (flat window, L = L_lambda)
  G = (a L^2)^-1 sum_{gamma in I'} v_gamma v_gamma^T,  a = ||phi||^2/L

All zeros here are on the line, so G is PSD by construction; the
informative quantity is HOW CLOSE to singular the truncation gets as the
support grows — the conditioning floor. NOTE: when d = lambda*d_1 exceeds
the number N' of zeros in I' (crossover lambda ~ 1.11 at this T), the
matrix is rank-deficient for dimensional reasons (overcomplete Gabor
basis); the meaningful floor is then the smallest POSITIVE eigenvalue.

Metrics per lambda: d, N', numerical rank, eigenvalue percentiles,
smallest positive eigenvalue lambda_min+ and the ratios
lambda_min+/lambda_max, lambda_min+/mean(eig>0); log10 of the ratio =
conditioning exponent. A superexponential decay of the floor (cf. Kim et
al. 2026 for Suzuki's operator) would be a complexity wall for positivity
certification (Groskin's B_T budget: resolving lambda_min needs precision
1/lambda_min).

Outputs: research/numerics/data/E-006-positivity-horizon.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-006-positivity-horizon.csv"
JSON = DATA / "E-006-positivity-horizon.json"

T = 2000.0
LAMBDAS = (0.5, 1.0, 1.5, 2.0)


def chi(u: np.ndarray) -> np.ndarray:
    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, xi_max_needed: float, pad: int = 8, M0: int = 4096):
    M = int(2 ** np.ceil(np.log2(max(M0, L * xi_max_needed / np.pi))))
    du = L / M
    u = -L / 2 + du * np.arange(M)
    phi = chi(L / 2 + u) * chi(L / 2 - u)  # flat window psi0
    a = np.sum(phi**2) * du / L
    full = np.zeros(pad * M)
    full[: M // 2] = phi[M // 2:]
    full[-M // 2:] = phi[: M // 2]
    F = np.fft.fft(full) * du
    xi = 2 * np.pi * np.fft.fftfreq(pad * M, d=du)
    order = np.argsort(xi)
    return CubicSpline(xi[order], F[order]), a


def main() -> int:
    if not ZEROS_CSV.exists():
        print("E-002 zeros not found"); 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)
    ell = np.log(T / (2 * np.pi))
    sqT = np.sqrt(T)
    in_Ip = (gam >= T - sqT) & (gam < 2 * T + sqT)
    zeros_Ip = gam[in_Ip]
    Np = len(zeros_Ip)
    print(f"T={T:.0f}  ell={ell:.4f}  zeros in I'={Np}", flush=True)

    rows = []
    for lam in LAMBDAS:
        t1 = time.time()
        L = lam * ell
        h = 2 * np.pi / L
        d = int(L * T / (2 * np.pi))
        alpha = T + h * np.arange(d)
        phi_hat, a = build_phi_hat(L, xi_max_needed=2 * T + sqT)
        Z = zeros_Ip[:, None] - alpha[None, :]
        V = np.real(phi_hat(Z))
        G = (V.T @ V) / (a * L**2)
        ev = eigvalsh(G)
        lam_max = float(ev[-1])
        # numerical rank + smallest positive eigenvalue
        noise = 1e-10 * lam_max
        pos = ev[ev > noise]
        rank_num = int(len(pos))
        lam_min_pos = float(pos[0]) if len(pos) else 0.0
        mean_pos = float(pos.mean()) if len(pos) else 0.0
        ratio_max = lam_min_pos / lam_max if lam_max else 0.0
        ratio_mean = lam_min_pos / mean_pos if mean_pos else 0.0
        rows.append({
            "lambda": lam, "d": d, "N_primes_zeros": Np,
            "dim_deficit": max(0, d - Np),
            "rank_numeric": rank_num,
            "lambda_max": f"{lam_max:.4e}",
            "lambda_min_pos": f"{lam_min_pos:.4e}",
            "minpos_over_max": f"{ratio_max:.3e}",
            "minpos_over_mean": f"{ratio_mean:.3e}",
            "log10_minpos_over_mean": f"{np.log10(ratio_mean) if ratio_mean>0 else float('-inf'):.2f}",
            "pct01_eig": f"{float(np.percentile(ev, 1)):.3e}",
            "pct10_eig": f"{float(np.percentile(ev, 10)):.3e}",
            "pct50_eig": f"{float(np.percentile(ev, 50)):.3e}",
        })
        print(f"lam={lam}: d={d} rank={rank_num} lam_max={lam_max:.3e} "
              f"lam_min+={lam_min_pos:.3e} min+/mean={ratio_mean:.2e} "
              f"({time.time()-t1:.0f}s)", flush=True)
        del G, V, Z

    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-006",
        "task": "3.1 / Q-3.4b (positivity horizon)",
        "description": "Eigenvalue floor of Gabor-truncated Weil form vs "
                       "bandwidth lambda; conditioning of positivity "
                       "certification (cf. Groskin B_T budget, Kim et al. "
                       "superexponential decay)",
        "script": "research/numerics/scripts/E-006-positivity-horizon.py",
        "parameters": {"T": T, "lambdas": list(LAMBDAS), "window": "psi0",
                       "zeros_source": "E-002", "precision": "float64"},
        "caveats": ["all zeros on-line => PSD by construction; the floor "
                    "measures conditioning, not failure",
                    "d > N' => dimensional rank deficiency (overcomplete "
                    "basis), not form degeneracy"],
        "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())