E-022-spectral-measure.py

{
  "experiment": "E-022",
  "task": "2.10 (D-035 D2): spectral-measure convergence",
  "script": "research/numerics/scripts/E-022-spectral-measure.py",
  "parameters": {
    "a_list": [
      4.0,
      6.0,
      8.0,
      10.0,
      12.0
    ],
    "h": 0.01,
    "lam": -1.0,
    "z_max": 300.0,
    "z_metric": 280.0,
    "sigma_bump": 8.0,
    "n_grid": 40000
  },
  "checks": {
    "V0_schur_vs_solve_v": true,
    "V1_weights_positive": true,
    "V2_refinement": true,
    "V3_lattice_and_z1law": true
  },
  "key_values": {
    "4": {
      "n_zeros": 381,
      "z1": 1.1260734006983593,
      "z1_times_a": 4.504293602793437,
      "lattice_expected": 381.9718634205488,
      "mass_total_Zmetric": 184.02432245308052,
      "O1_cumulative_Linf": 0.030472775278319375,
      "O1_cumulative_L2": 0.015369322417422383,
      "O2_logw_slope": 0.004773175983096674,
      "O2_logw_R2": 0.17155831473911887,
      "O2_logw_std_resid": 0.8476672146023767,
      "O3_weak_metric": 0.010527875427328966,
      "O5_mass_near_gamma": 1.6069393990772483,
      "O5_mass_near_random": 1.239913792291073,
      "O5_ratio": 1.2960089718076262,
      "elapsed_s": 2.9
    },
    "6": {
      "n_zeros": 572,
      "z1": 0.7499016225074712,
      "z1_times_a": 4.499409735044827,
      "lattice_expected": 572.9577951308232,
      "mass_total_Zmetric": 184.01489106241632,
      "O1_cumulative_Linf": 0.02723380344436621,
      "O1_cumulative_L2": 0.015298102993183555,
      "O2_logw_slope": 0.004573680792971365,
      "O2_logw_R2": 0.13177706690264612,
      "O2_logw_std_resid": 0.9483241546589634,
      "O3_weak_metric": 0.010569164981132785,
      "O5_mass_near_gamma": 1.6173747450375149,
      "O5_mass_near_random": 1.257329769139859,
      "O5_ratio": 1.2863568371120038,
      "elapsed_s": 4.8
    },
    "8": {
      "n_zeros": 763,
      "z1": 0.5621012628096893,
      "z1_times_a": 4.4968101024775144,
      "lattice_expected": 763.9437268410976,
      "mass_total_Zmetric": 184.13312237431055,
      "O1_cumulative_Linf": 0.027258094320286474,
      "O1_cumulative_L2": 0.015486596743652112,
      "O2_logw_slope": 0.004185396789595135,
      "O2_logw_R2": 0.10532818940743971,
      "O2_logw_std_resid": 0.986518773772079,
      "O3_weak_metric": 0.010568959142405882,
      "O5_mass_near_gamma": 1.6196075350475192,
      "O5_mass_near_random": 1.2561234028070363,
      "O5_ratio": 1.2893697636937673,
      "elapsed_s": 8.2
    },
    "10": {
      "n_zeros": 954,
      "z1": 0.44959701245436284,
      "z1_times_a": 4.4959701245436285,
      "lattice_expected": 954.9296585513721,
      "mass_total_Zmetric": 184.1239991524169,
      "O1_cumulative_Linf": 0.024844169503423608,
      "O1_cumulative_L2": 0.015446208652370818,
      "O2_logw_slope": 0.0038097364559334847,
      "O2_logw_R2": 0.08653440341105123,
      "O2_logw_std_resid": 1.0006504172097648,
      "O3_weak_metric": 0.010581710686957842,
      "O5_mass_near_gamma": 1.6290317856739722,
      "O5_mass_near_random": 1.2562149964535079,
      "O5_ratio": 1.2967778527345915,
      "elapsed_s": 13.4
    },
    "12": {
      "n_zeros": 1145,
      "z1": 0.3746087478472408,
      "z1_times_a": 4.495304974166889,
      "lattice_expected": 1145.9155902616465,
      "mass_total_Zmetric": 184.10483484645468,
      "O1_cumulative_Linf": 0.028584131798192303,
      "O1_cumulative_L2": 0.015417413187446537,
      "O2_logw_slope": 0.003682406300531162,
      "O2_logw_R2": 0.08233450272296328,
      "O2_logw_std_resid": 0.9935774298804366,
      "O3_weak_metric": 0.010592268896428007,
      "O5_mass_near_gamma": 1.6165573705799765,
      "O5_mass_near_random": 1.2532365071394718,
      "O5_ratio": 1.2899060643148588,
      "elapsed_s": 20.9
    },
    "O4_trend": {
      "a": [
        4.0,
        6.0,
        8.0,
        10.0,
        12.0
      ],
      "d_a": [
        0.010527875427328966,
        0.010569164981132785,
        0.010568959142405882,
        0.010581710686957842,
        0.010592268896428007
      ],
      "monotone_decreasing": false
    },
    "V0_max_dev": 3.560323900663321e-14,
    "V2_rel": 0.0011243096862693027
  },
  "elapsed_seconds": 743.1,
  "timestamp_utc": "2026-08-19T17:19:38.162052+00:00"
}

"""E-022 — task 2.10 (D-035 deliverable D2): spectral-measure
convergence of Suzuki's D_{a,pi} toward the gamma-target measure.

For each a in {4,6,8,10,12} (theta=pi, all zeros real, E-016 F1):
  - locate zeros z_k of W(a,pi;z) on (0, Z_MAX] (E-016/E-019 route:
    grid sign-changes of g(z) = z Im V+ - Re V+ + brentq);
  - compute the Clark/de Branges spectral weights
        w_k = 1 / K_a(z_k, z_k),
    K_a(z,z) = Re(b_z^H v),  v solves the BORDERED (mean-zero, L^2_0)
    Galerkin system  T v + mu e = b_z,  e.v = 0,
    T = Q - lam M (lam = -1),  b_z = h sinc^2(zh/2pi) e^{-i z x}
    (exact hat-basis quadrature), via Schur reduction on shared
    u1 = T^{-1} e: v = T^{-1}(b_z - mu e), mu = (e.u0)/(e.u1).
    Toeplitz solves (scipy solve_toeplitz); V0 validates the Schur
    solve against E16.solve_v (bordered dense) for b[+1].

Target measure: unit atoms at the Riemann zeros gamma_k (mpmath
zetazero, dps 30, float64 sufficient for statistics).

Observables:
  O1. normalized cumulative mass M_a(x)/M_a(Z') vs N_gamma(x)/N(Z'),
      Z' = 280 (Linfty + L2 discrepancy).
  O2. weight profile: linear fit log w_k vs z_k (slope, spread R^2).
  O3. weak metric: Gaussian bumps f_j (30 log-spaced centers in
      [15,250], sigma=8): d_a = max_j |mu~_a(f_j) - nu~(f_j)|,
      normalized by total mass in (0, Z'].
  O4. d_a vs a  -> convergence-rate estimate (headline).
  O5. gamma-concentration: mean mu_a-mass within +/-1 of each gamma
      vs within +/-1 of random points (does the WEIGHT carry gamma
      information even though positions don't lock? — E-019 B2 found
      chance-level POSITION proximity).

Checks:
  V0. Schur-Toeplitz bordered solve vs E16.solve_v (b[+1]): < 1e-8.
  V1. all weights > 0.
  V2. refinement a=4, h=0.005: weight rel dev on first 50 zeros < 5e-3.
  V3. z1*a in [4.3, 4.7] (stiff law, E-019) and count = maximal
      lattice (a*Z_MAX/pi +/- 2) at every a.
  V4. total mass M_a(Z') recorded vs a (scaling documented).

Outputs: research/numerics/data/E-022-spectral-measure.csv (+ .json)
LONG RUN: detached + log file (per STATE.md protocol).
"""
import importlib.util
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import mpmath as mp
import numpy as np
import scipy.linalg
import scipy.optimize

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

_spec = importlib.util.spec_from_file_location(
    "E16", Path(__file__).with_name("E-016-herglotz-fine-structure.py"))
E16 = importlib.util.module_from_spec(_spec)
_spec.loader.exec_module(E16)                     # runs the e^20 sieve

sys.path.insert(0, str(Path(__file__).parent))
from prime_backend_v2 import PrimeBackendV2       # noqa: E402

H_MAIN = E16.H_MAIN          # 0.01
LAM = E16.LAM_SHIFT          # -1.0
A_LIST = [4.0, 6.0, 8.0, 10.0, 12.0]
Z_MAX = 300.0
Z_METRIC = 280.0
N_GRID = 40000
SIGMA_BUMP = 8.0
mp.mp.dps = 40

XN, XW = E16.XN, E16.XW
hat_FT = E16.hat_FT


def collect_t_points(a, h):
    """Copy of E-019's collector (E-016 S_full/quad segment structure)."""
    n_elem = int(round(2 * a / h))
    pts = []

    def add_quad(lo, hi, peak):
        n_lo = int(np.ceil(np.exp(lo)))
        n_hi = int(np.floor(np.exp(hi)))
        if n_hi - n_lo > 64:
            for aa, bb in ((lo, peak), (peak, hi)):
                if bb - aa > 0:
                    mid, half = 0.5 * (aa + bb), 0.5 * (bb - aa)
                    pts.extend(mid + half * XN)
            return
        segs = sorted({lo, peak, hi} |
                      {float(np.log(n))
                       for n in range(max(n_lo, 2), n_hi + 1)
                       if lo < np.log(n) < hi})
        for p0, p1 in zip(segs, segs[1:]):
            if p1 - p0 > 0:
                mid, half = 0.5 * (p0 + p1), 0.5 * (p1 - p0)
                pts.extend(mid + half * XN)

    for k in range(n_elem):
        if k == 0:
            add_quad(0.0, h, 0.0)
        else:
            add_quad((k - 1) * h, (k + 1) * h, k * h)
    return np.array(pts)


def locate_zeros(cp, x, h, z_max=Z_MAX, n_grid=N_GRID):
    zg = np.linspace(1e-3, z_max, n_grid)
    Vpr, _ = E16.V_and_dV(zg, cp, x, h)
    g = zg * Vpr.imag - Vpr.real

    def g_direct(zz):
        V = hat_FT(np.array([zz]), h)[0] * complex(
            np.exp(1j * zz * x) @ cp)
        return zz * V.imag - V.real

    idx = np.flatnonzero(np.sign(g[:-1]) * np.sign(g[1:]) < 0)
    zeros = []
    for i in idx:
        try:
            zeros.append(scipy.optimize.brentq(g_direct, zg[i], zg[i + 1],
                                               xtol=1e-13, rtol=1e-13))
        except ValueError:
            pass
    return np.array(zeros)


class ToeplitzBordered:
    """Schur-reduced bordered solve for T = toeplitz(tc), constraint
    e.v = 0 (e = ones):  v = u0 - mu u1, mu = sum(u0)/sum(u1)."""

    def __init__(self, tc):
        self.tc = tc
        n = len(tc)
        self.u1 = scipy.linalg.solve_toeplitz((tc, tc), np.ones(n))
        self.s1 = float(self.u1.sum())

    def solve(self, b):
        u0 = scipy.linalg.solve_toeplitz((self.tc, self.tc), b)
        return u0 - (u0.sum() / self.s1) * self.u1


def weights_for(a, h, zeros, x):
    """K_k and w_k for the given zero list (build(a,h) cached)."""
    _, _, Q, M, _b = E16.build(a, h)
    n = Q.shape[0]
    mc = np.zeros(n)
    mc[0] = 2 * h / 3
    if n > 1:
        mc[1] = h / 6
    tc = Q[:, 0] - LAM * mc
    solver = ToeplitzBordered(tc)
    Ks = np.empty(len(zeros))
    for i, z in enumerate(zeros):
        bz = hat_FT(np.array([z]), h)[0] * np.exp(-1j * z * x)
        v = solver.solve(bz)
        Ks[i] = float(np.real(np.conj(bz) @ v))
    return Ks


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

    # ---- backend v2 (e^24) ------------------------------------------
    print("collecting quadrature t-points ...", flush=True)
    tq = np.concatenate([collect_t_points(a, H_MAIN) for a in A_LIST]
                        + [collect_t_points(4.0, 0.005)])
    print(f"  {tq.size} points; building/loading backend ...", flush=True)
    CACHE = DATA / "E-022-backend-cache.npz"
    be = PrimeBackendV2(tq)
    if CACHE.exists() and "--force" not in sys.argv:
        be.load_cache(CACHE)
        print("  backend loaded from cache", flush=True)
    else:
        be.build(verbose=True)
        be.save_cache(CACHE)
    E16._BUILD.clear()
    E16.prime_sum_vec = be.prime_sum

    # ---- gamma target ------------------------------------------------
    print("target: first 200 Riemann zeros (zetazero) ...", flush=True)
    mp.mp.dps = 30
    gam = np.array([float(mp.zetazero(k).imag) for k in range(1, 201)])
    gam = gam[gam <= Z_METRIC]

    # ---- main sweep --------------------------------------------------
    rows = []
    report = {"per_a": {}}
    v0_max = 0.0
    v1_ok = True
    v3_ok = True
    zc = np.linspace(10.0, Z_METRIC, 800)
    centers = np.geomspace(15.0, 250.0, 30)

    for a in A_LIST:
        ta = time.time()
        x, h, Q, M, b = E16.build(a, H_MAIN)
        n = len(x)
        # V0: Schur-Toeplitz vs E16.solve_v on b[+1]
        mc = np.zeros(n)
        mc[0] = 2 * h / 3
        mc[1] = h / 6
        tc = Q[:, 0] - LAM * mc
        solver = ToeplitzBordered(tc)
        cp_ref, _T = E16.solve_v(Q, M, b[+1.0], LAM, n)
        cp_schur = solver.solve(b[+1.0])
        dev = float(np.max(np.abs(cp_ref - cp_schur))
                    / max(1.0, float(np.max(np.abs(cp_ref)))))
        v0_max = max(v0_max, dev)

        cp, _ = E16.solve_v(Q, M, b[+1.0], LAM, n)
        cm, _ = E16.solve_v(Q, M, b[-1.0], LAM, n)
        zeros = locate_zeros(cp, x, h)

        Ks = weights_for(a, H_MAIN, zeros, x)
        ws = 1.0 / Ks
        if np.any(Ks <= 0):
            v1_ok = False

        z1 = float(zeros[0])
        lat = a * Z_MAX / np.pi
        if not (4.3 < z1 * a < 4.7) or abs(len(zeros) - lat) > 2:
            v3_ok = False

        # --- metrics (zeros in (0, Z_METRIC])
        zm = zeros[zeros <= Z_METRIC]
        wm = ws[zeros <= Z_METRIC]
        Mtot = float(wm.sum())
        # O1: cumulative shape vs unit-weight gamma measure
        Ma = np.array([wm[zm <= t].sum() for t in zc]) / Mtot
        Ng = np.array([np.sum(gam <= t) for t in zc]) / len(gam)
        o1_linf = float(np.max(np.abs(Ma - Ng)))
        o1_l2 = float(np.sqrt(np.mean((Ma - Ng) ** 2)))
        # O2: weight profile
        A_ = np.vstack([zm, np.ones_like(zm)]).T
        slope, intercept = np.linalg.lstsq(A_, np.log(wm), rcond=None)[0]
        pred = slope * zm + intercept
        ss_res = float(np.sum((np.log(wm) - pred) ** 2))
        ss_tot = float(np.sum((np.log(wm) - np.log(wm).mean()) ** 2))
        r2 = 1.0 - ss_res / max(ss_tot, 1e-300)
        # O3: Gaussian-bump weak metric
        def bump(z, c):
            return np.exp(-0.5 * ((z - c) / SIGMA_BUMP) ** 2)
        mu_f = np.array([np.dot(wm, bump(zm, c)) for c in centers]) / Mtot
        nu_f = np.array([bump(gam, c).sum() for c in centers]) / len(gam)
        o3 = float(np.max(np.abs(mu_f - nu_f)))
        # O5: gamma-concentration of mass
        mass_near_g = float(np.mean([wm[np.abs(zm - g) <= 1.0].sum()
                                     for g in gam]))
        rng = np.random.default_rng(7)
        rpts = rng.uniform(15.0, 250.0, len(gam))
        mass_near_r = float(np.mean([wm[np.abs(zm - r) <= 1.0].sum()
                                     for r in rpts]))

        for k, (zz, ww) in enumerate(zip(zeros, ws), 1):
            rows.append((a, k, zz, ww))
        report["per_a"][f"{a:g}"] = {
            "n_zeros": int(len(zeros)), "z1": z1, "z1_times_a": z1 * a,
            "lattice_expected": lat,
            "mass_total_Zmetric": Mtot,
            "O1_cumulative_Linf": o1_linf, "O1_cumulative_L2": o1_l2,
            "O2_logw_slope": float(slope), "O2_logw_R2": r2,
            "O2_logw_std_resid": float(np.sqrt(ss_res / max(len(zm) - 2, 1))),
            "O3_weak_metric": o3,
            "O5_mass_near_gamma": mass_near_g,
            "O5_mass_near_random": mass_near_r,
            "O5_ratio": mass_near_g / max(mass_near_r, 1e-300),
            "elapsed_s": round(time.time() - ta, 1),
        }
        print(f"a={a:g}: {len(zeros)} zeros, z1a={z1 * a:.3f}, "
              f"mass={Mtot:.3e}, O3={o3:.4f}, "
              f"O5={mass_near_g:.3e}/{mass_near_r:.3e} "
              f"({report['per_a'][f'{a:g}']['elapsed_s']}s)", flush=True)

    # ---- V2: refinement at a=4, h=0.005 ------------------------------
    print("V2: refinement a=4 h=0.005 (first 50 zeros) ...", flush=True)
    x5, h5, Q5, M5, b5 = E16.build(4.0, 0.005)
    cp5, _ = E16.solve_v(Q5, M5, b5[+1.0], LAM, len(x5))
    zeros5 = locate_zeros(cp5, x5, h5)
    K5 = weights_for(4.0, 0.005, zeros5[:50], x5)
    x01, _, _, _, _ = E16.build(4.0, H_MAIN)
    K4 = weights_for(4.0, H_MAIN, zeros5[:50], x01)  # same z's, h=0.01
    v2_rel = float(np.max(np.abs(K5 - K4) / np.abs(K5)))
    print(f"  V2 max rel dev K(h=0.005) vs K(h=0.01): {v2_rel:.2e}",
          flush=True)

    # ---- O4: convergence trend --------------------------------------
    avals = np.array(A_LIST)
    dvals = np.array([report["per_a"][f"{a:g}"]["O3_weak_metric"]
                      for a in A_LIST])
    o4 = {"a": list(avals), "d_a": list(dvals),
          "monotone_decreasing": bool(np.all(np.diff(dvals) < 0))}

    checks = {
        "V0_schur_vs_solve_v": bool(v0_max < 1e-8),
        "V1_weights_positive": bool(v1_ok),
        "V2_refinement": bool(v2_rel < 5e-3),
        "V3_lattice_and_z1law": bool(v3_ok),
    }
    print("checks:", checks, flush=True)

    import csv as csvmod
    with open(CSV, "w", newline="") as f:
        w = csvmod.writer(f)
        w.writerow(["a", "k", "z_k", "w_k"])
        for a, k, zz, ww in rows:
            w.writerow([f"{a:g}", k, f"{zz:.12f}", f"{ww:.10e}"])

    meta = {
        "experiment": "E-022",
        "task": "2.10 (D-035 D2): spectral-measure convergence",
        "script": "research/numerics/scripts/E-022-spectral-measure.py",
        "parameters": {"a_list": A_LIST, "h": H_MAIN, "lam": LAM,
                       "z_max": Z_MAX, "z_metric": Z_METRIC,
                       "sigma_bump": SIGMA_BUMP, "n_grid": N_GRID},
        "checks": checks,
        "key_values": {**report["per_a"], "O4_trend": o4,
                       "V0_max_dev": v0_max, "V2_rel": v2_rel},
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
    }
    JSON.write_text(json.dumps(meta, indent=2))
    print(f"wrote {CSV.name}, {JSON.name}; elapsed "
          f"{meta['elapsed_seconds']}s", flush=True)
    return 0 if all(checks.values()) else 1


if __name__ == "__main__":
    sys.exit(main())