E-020b-window-optimum.py

{
  "experiment": "E-020b",
  "task": "1.6: window-class optimum for kernel D1",
  "script": "research/numerics/scripts/E-020b-window-optimum.py",
  "parameters": {
    "K_poly": 6,
    "n_grid": 801
  },
  "checks": {
    "C1_control_recovers_MT": true,
    "C2_D1_optimum_sane": true
  },
  "key_values": {
    "R_opt_x": 1.3274994741768826,
    "R_MT": 1.3274992963205885,
    "R_opt_D1": 1.131359211780498,
    "R_quartic_D1": 1.1313601952861119,
    "cert_opt_D1": 0.868640788219502,
    "opt_coeffs_D1": [
      1.4747096406032818,
      -0.11822069827325471,
      -0.3332111306955125,
      -0.05188711282919025,
      0.018718518018636963,
      0.014002590280189901,
      -0.014950796476394187
    ]
  },
  "elapsed_seconds": 5.9,
  "timestamp_utc": "2026-08-18T20:22:04.723162+00:00"
}

"""E-020b — window-class optimum for kernel D1 (task 1.6, attack-side
bound): is the paper's quartic window (0.86864) near-optimal for the
xi' certificate?

Minimize the Rayleigh quotient
  R_K(psi) = [∫ψ² + ∬ K(|u-v|) ψ(u)ψ(v)] / (∫ψ)²,  u,v in [-1/2, 1/2]
over even polynomial windows psi(x) = sum_{k=0}^{K} a_k (2x)^{2k}
(multi-start L-BFGS; K = 6). Certificate value = 2 - R.

CONTROL (kernel x): optimum should recover the Montgomery-Taylor value
R = 0.5 + (1/sqrt2) cot(1/sqrt2) = 1.32749930 (E-001) — validates the
optimizer. TARGET (kernel D1, E-013): optimum vs the paper's quartic
R = 2 - 0.86864152.

Checks:
  C1. control: |R_opt(x-kernel) - 1.32749930| < 2e-4.
  C2. D1: R_opt(D1) <= R_quartic + 1e-8 (optimizer does at least as
      well as the paper's window) and R_opt >= R_quartic - 5e-2
      (sanity: no runaway).
Outputs: research/numerics/data/E-020b-window-optimum.csv (+ .json)
"""
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import numpy as np
import scipy.optimize

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

K_TAIL = 40
K_POLY = 6
N_GRID = 801
R_MT = 0.5 + (1 / np.sqrt(2)) / np.tan(1 / np.sqrt(2))   # 1.32749930..
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


# quadrature grid
xg = np.linspace(-0.5, 0.5, N_GRID)
wgt = np.full(N_GRID, 1.0 / (N_GRID - 1))
wgt[0] = wgt[-1] = 0.5 / (N_GRID - 1)          # trapezoid
BASIS = np.stack([(2 * xg) ** (2 * k) for k in range(K_POLY + 1)], 1)


def make_R(kernel_vals):
    Km = kernel_vals[np.abs(np.subtract.outer(xg, xg)).astype(int)]
    # note: |u-v| in [0,1]; index into precomputed kernel table
    return Km


DX = np.abs(np.subtract.outer(xg, xg))
K_TAB_X = DX.copy()
K_TAB_D1 = d1(DX)


def R_of(a, Km):
    psi = BASIS @ a
    i1 = np.sum(wgt * psi**2)
    i2 = np.einsum("i,ij,j->", wgt * psi, Km, wgt * psi)
    n = np.sum(wgt * psi)
    return (i1 + i2) / n**2


def optimize(Km, starts=24):
    best = None
    seeds = [np.eye(K_POLY + 1)[k] for k in range(K_POLY + 1)]
    seeds += [rng.normal(size=K_POLY + 1) for _ in range(starts)]
    for a0 in seeds:
        res = scipy.optimize.minimize(
            lambda a: R_of(a, Km), a0, method="L-BFGS-B",
            options={"maxiter": 500})
        if best is None or res.fun < best.fun:
            best = res
    return best


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()

    print("control: kernel x ...", flush=True)
    b_x = optimize(K_TAB_X)
    r_x = float(b_x.fun)
    print(f"  R_opt = {r_x:.8f} (MT = {R_MT:.8f})", flush=True)

    print("target: kernel D1 ...", flush=True)
    b_d = optimize(K_TAB_D1)
    r_d = float(b_d.fun)
    # paper windows for reference
    r_flat = R_of(np.eye(K_POLY + 1)[0], K_TAB_D1)
    aq = np.zeros(K_POLY + 1)
    aq[0], aq[1], aq[2] = 1.0, -0.07, -0.255
    r_quart = R_of(aq, K_TAB_D1)
    print(f"  R_opt = {r_d:.8f}  R_flat = {r_flat:.8f}  "
          f"R_quartic = {r_quart:.8f}", flush=True)
    print(f"  certificate 2-R: opt {2-r_d:.5f} vs paper 0.86864",
          flush=True)

    c1 = bool(abs(r_x - R_MT) < 2e-4)
    c2 = bool(r_d <= r_quart + 1e-8 and r_d >= r_quart - 5e-2)
    checks = {"C1_control_recovers_MT": c1, "C2_D1_optimum_sane": c2}
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("kernel,R_opt,certificate_2_minus_R,reference\n")
        f.write(f"x,{r_x:.10f},{2-r_x:.10f},MT={R_MT:.10f}\n")
        f.write(f"D1,{r_d:.10f},{2-r_d:.10f},quartic={r_quart:.10f}\n")
    meta = {
        "experiment": "E-020b",
        "task": "1.6: window-class optimum for kernel D1",
        "script": "research/numerics/scripts/E-020b-window-optimum.py",
        "parameters": {"K_poly": K_POLY, "n_grid": N_GRID},
        "checks": checks,
        "key_values": {"R_opt_x": r_x, "R_MT": R_MT, "R_opt_D1": r_d,
                       "R_quartic_D1": r_quart,
                       "cert_opt_D1": 2 - r_d,
                       "opt_coeffs_D1": [float(v) for v in b_d.x]},
        "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())