E-022b-valley-analysis.py

{
  "experiment": "E-022b",
  "task": "2.10 analysis: K-diagonal valleys + mass concentration",
  "script": "research/numerics/scripts/E-022b-valley-analysis.py",
  "checks": {
    "V1_valleys_at_gamma": true,
    "V2_mass_concentrated": true,
    "V4_fit_quality": true
  },
  "key_values": {
    "4": {
      "valley_depth": 2.1593140478568413,
      "logK_at_gamma": -0.106956845297883,
      "logK_at_midpoints": 2.0523572025589583,
      "mass_frac_pm0.1": 0.6338714893637662,
      "uniform_expectation": 0.09000000000000001,
      "mass_ratio_vs_uniform": 7.043016548486291,
      "within_bin_meanlogw": {
        "z10-60_d0.0-0.1": 0.05321233524569066,
        "z10-60_d0.3-0.7": -1.7997896503200193,
        "z120-200_d0.0-0.1": 0.10247759965464488,
        "z120-200_d0.3-0.7": -1.7570525979064509
      },
      "elapsed_s": 2.9
    },
    "6": {
      "valley_depth": 2.5087271936533595,
      "logK_at_gamma": -0.04029578260580941,
      "logK_at_midpoints": 2.46843141104755,
      "mass_frac_pm0.1": 0.6986826134175042,
      "uniform_expectation": 0.09000000000000001,
      "mass_ratio_vs_uniform": 7.7631401490833785,
      "within_bin_meanlogw": {
        "z10-60_d0.0-0.1": -0.10145611105454076,
        "z10-60_d0.3-0.7": -2.4097253102207548,
        "z120-200_d0.0-0.1": 0.029192545882919238,
        "z120-200_d0.3-0.7": -2.314439077156247
      },
      "elapsed_s": 5.1
    },
    "8": {
      "valley_depth": 2.7352201581877043,
      "logK_at_gamma": 0.030554703235580436,
      "logK_at_midpoints": 2.7657748614232847,
      "mass_frac_pm0.1": 0.7251109123178877,
      "uniform_expectation": 0.09000000000000001,
      "mass_ratio_vs_uniform": 8.056787914643197,
      "within_bin_meanlogw": {
        "z10-60_d0.0-0.1": -0.1389960979715335,
        "z10-60_d0.3-0.7": -2.7392603482112916,
        "z120-200_d0.0-0.1": -0.20043183472111312,
        "z120-200_d0.3-0.7": -2.6451890375288283
      },
      "elapsed_s": 7.9
    },
    "10": {
      "valley_depth": 2.8706885758089333,
      "logK_at_gamma": 0.11948603922925684,
      "logK_at_midpoints": 2.9901746150381903,
      "mass_frac_pm0.1": 0.7251715332604941,
      "uniform_expectation": 0.09000000000000001,
      "mass_ratio_vs_uniform": 8.057461480672156,
      "within_bin_meanlogw": {
        "z10-60_d0.0-0.1": -0.07371951844373904,
        "z10-60_d0.3-0.7": -2.974563705268562,
        "z120-200_d0.0-0.1": -0.2233480705431757,
        "z120-200_d0.3-0.7": -2.8788847145178686
      },
      "elapsed_s": 11.4
    },
    "12": {
      "valley_depth": 2.9442229398281423,
      "logK_at_gamma": 0.22834844487190686,
      "logK_at_midpoints": 3.172571384700049,
      "mass_frac_pm0.1": 0.7241062961289328,
      "uniform_expectation": 0.09000000000000001,
      "mass_ratio_vs_uniform": 8.045625512543698,
      "within_bin_meanlogw": {
        "z10-60_d0.0-0.1": -0.1831905054381875,
        "z10-60_d0.3-0.7": -3.1605455546413452,
        "z120-200_d0.0-0.1": -0.35321115246295404,
        "z120-200_d0.3-0.7": -3.0663074224107776
      },
      "elapsed_s": 15.5
    },
    "V4_saturation_fit": {
      "Delta_inf": 3.3338697242432738,
      "c": 4.760242352939846,
      "max_resid": 0.03176880509994007
    },
    "valley_depths": {
      "4.0": 2.1593140478568413,
      "6.0": 2.5087271936533595,
      "8.0": 2.7352201581877043,
      "10.0": 2.8706885758089333,
      "12.0": 2.9442229398281423
    }
  },
  "elapsed_seconds": 59.7,
  "timestamp_utc": "2026-08-19T17:37:47.138364+00:00"
}

"""E-022b — task 2.10 analysis stage: reproducing-kernel valley structure
and gamma-mass concentration (post-hoc analysis of E-022).

E-022 found the coarse weak metric flat (O3) but a persistent mass
elevation near the true zeros. This script resolves the mechanism:

  K_a(z,z) := Re(b_z^H v), v the bordered (L^2_0) Galerkin solution of
  T_a v = e^{-izx} — the diagonal of the model-space reproducing kernel.

Measurements (a in {4,6,8,10,12}, theta=pi):
  V1. valley depth: Delta(a) = mean log K at gamma-gap midpoints minus
      mean log K at gammas, z-grid step 0.05 in [14,60).
  V2. mass concentration: fraction of the E-022 spectral mass within
      +/-0.1 of a true zero (vs uniform-lattice expectation).
  V3. artifact control: mean log w vs distance-to-nearest-gamma WITHIN
      z-bins (controls the gamma-density/weight-trend confound that
      naive O5 fails — shifted-gamma null test showed the raw O5
      elevation is mostly distributional).
  V4. saturation model: Delta(a) ~ Delta_inf - c/a fit.

Outputs: research/numerics/data/E-022b-valley.csv (+ .json)
"""
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

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-022b-valley.csv"
JSON = DATA / "E-022b-valley.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 = E16.H_MAIN
LAM = E16.LAM_SHIFT
A_LIST = (4.0, 6.0, 8.0, 10.0, 12.0)
XN = E16.XN


def collect_t_points(a, h):
    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 make_K1(a):
    x, h, Q, M, _b = E16.build(a, H)
    n = Q.shape[0]
    mc = np.zeros(n)
    mc[0] = 2 * h / 3
    mc[1] = h / 6
    tc = Q[:, 0] - LAM * mc
    u1 = scipy.linalg.solve_toeplitz((tc, tc), np.ones(n))
    s1 = u1.sum()

    def K1(z):
        bz = E16.hat_FT(np.array([z]), h)[0] * np.exp(-1j * z * x)
        u0 = scipy.linalg.solve_toeplitz((tc, tc), bz)
        v = u0 - (u0.sum() / s1) * u1
        return float(np.real(np.conj(bz) @ v))

    return K1


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

    tq = np.concatenate([collect_t_points(a, H) for a in A_LIST])
    be = PrimeBackendV2(tq)
    be.load_cache(DATA / "E-022-backend-cache.npz")
    E16._BUILD.clear()
    E16.prime_sum_vec = be.prime_sum

    mp.mp.dps = 30
    gam_all = np.array([float(mp.zetazero(k).imag)
                        for k in range(1, 201)])
    gam = gam_all[gam_all <= 280]

    # E-022 weights
    import csv as csvmod
    data = {}
    with open(DATA / "E-022-spectral-measure.csv") as f:
        for r in csvmod.DictReader(f):
            data.setdefault(float(r["a"]), []).append(
                (float(r["z_k"]), float(r["w_k"])))

    zgrid = np.arange(14.0, 60.0, 0.05)
    gsel = gam_all[(gam_all >= 14) & (gam_all < 60)]
    mids = (gsel[:-1] + gsel[1:]) / 2
    idxg = np.round((gsel - 14) / 0.05).astype(int)
    idxm = np.round((mids - 14) / 0.05).astype(int)

    per_a = {}
    grid_rows = []
    for a in A_LIST:
        ta = time.time()
        K1 = make_K1(a)
        lg = np.array([np.log(K1(z)) for z in zgrid])
        dip = float(lg[idxm].mean() - lg[idxg].mean())
        # mass fraction within +-0.1 of gamma (full window to 280)
        z = np.array([p[0] for p in data[a]])
        w = np.array([p[1] for p in data[a]])
        m = z <= 280
        zm, wm = z[m], w[m]
        idx = np.clip(np.searchsorted(gam, zm), 1, len(gam) - 1)
        d = np.minimum(np.abs(zm - gam[idx - 1]), np.abs(zm - gam[idx]))
        frac = float(wm[d <= 0.1].sum() / wm.sum())
        unif = float(0.2 * len(gam) / 280.0)
        # V3: within-z-bin mean log w by distance bin
        lw = np.log(wm)
        bins = {}
        for lo, hi in ((10, 60), (120, 200)):
            sel = (zm >= lo) & (zm < hi)
            for blo, bhi in ((0.0, 0.1), (0.3, 0.7)):
                s2 = sel & (d >= blo) & (d < bhi)
                if s2.sum() > 3:
                    bins[f"z{lo}-{hi}_d{blo}-{bhi}"] = float(lw[s2].mean())
        per_a[f"{a:g}"] = {
            "valley_depth": dip,
            "logK_at_gamma": float(lg[idxg].mean()),
            "logK_at_midpoints": float(lg[idxm].mean()),
            "mass_frac_pm0.1": frac,
            "uniform_expectation": unif,
            "mass_ratio_vs_uniform": frac / unif,
            "within_bin_meanlogw": bins,
            "elapsed_s": round(time.time() - ta, 1),
        }
        print(f"a={a:g}: dip={dip:.4f}, massfrac={frac:.4f} "
              f"({per_a[f'{a:g}']['elapsed_s']}s)", flush=True)
        for zz, ll in zip(zgrid, lg):
            grid_rows.append((a, zz, ll))

    # V4: saturation fit Delta(a) = Dinf - c/a
    av = np.array(A_LIST)
    dv = np.array([per_a[f"{a:g}"]["valley_depth"] for a in A_LIST])
    Afit = np.vstack([np.ones_like(av), -1.0 / av]).T
    Dinf, c = np.linalg.lstsq(Afit, dv, rcond=None)[0]
    resid = dv - (Dinf - c / av)
    print(f"V4: Delta_inf = {Dinf:.3f}, c = {c:.3f}, "
          f"max fit resid {np.abs(resid).max():.4f}", flush=True)

    checks = {
        "V1_valleys_at_gamma": bool(np.all(dv > 1.5)),
        "V2_mass_concentrated": bool(np.all(
            [per_a[f"{a:g}"]["mass_frac_pm0.1"] > 0.5 for a in A_LIST])),
        "V4_fit_quality": bool(np.abs(resid).max() < 0.1),
    }

    with open(CSV, "w", newline="") as f:
        wcsv = csvmod.writer(f)
        wcsv.writerow(["a", "z", "logK"])
        for a, zz, ll in grid_rows:
            wcsv.writerow([f"{a:g}", f"{zz:.3f}", f"{ll:.6f}"])

    meta = {
        "experiment": "E-022b",
        "task": "2.10 analysis: K-diagonal valleys + mass concentration",
        "script": "research/numerics/scripts/E-022b-valley-analysis.py",
        "checks": checks,
        "key_values": {**per_a,
                       "V4_saturation_fit": {"Delta_inf": float(Dinf),
                                             "c": float(c),
                                             "max_resid": float(
                                                 np.abs(resid).max())},
                       "valley_depths": dict(zip([str(a) for a in av],
                                                 dv.tolist()))},
        "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())