E-025-lambda-flow.py

{
  "experiment": "E-025",
  "task": "2.13: lambda-flow numerics (ground-state channel)",
  "script": "research/numerics/scripts/E-025-lambda-flow.py",
  "parameters": {
    "a_list": [
      4.0,
      8.0,
      12.0
    ],
    "ladder_rel": [
      0.0,
      0.5,
      0.9,
      0.99,
      0.999
    ],
    "z_max": 300.0
  },
  "checks": {
    "V1_lambda_a_positive_law": false,
    "V2_baseline_matches_E022": true,
    "V3_solve_residuals": false,
    "V4_sister_reported": true
  },
  "key_values": {
    "4": {
      "lambda_a": 5.147476282756457e-11,
      "lambda_law_ratio": 3.783885052191424e-08,
      "sister_corr": -0.2997460393487147,
      "sister_crossings": 57,
      "sister_cross_err": 16.71385138849778,
      "ladder": {
        "-1": {
          "lambda": -1.0,
          "Delta": 2.1593140478568413,
          "floor": -0.106956845297883,
          "plateau": 2.0523572025589583,
          "z1_times_a": 4.504293602793437,
          "n_zeros": 381,
          "lock_metric": 0.10068828882668726,
          "residual": 48.3075476074236,
          "mass_pm0.1": 116.6477713523963,
          "solve_resid": 1.588314368200427e-16
        },
        "0": {
          "lambda": 0.0,
          "Delta": 6.988266234191322,
          "floor": 4.4061523727558125,
          "plateau": 11.394418606947134,
          "z1_times_a": 6.9816622297266715,
          "n_zeros": 381,
          "lock_metric": 7.623144158814413e-05,
          "residual": 10.718615266368445,
          "mass_pm0.1": 127.66723750632949,
          "solve_resid": 5.831210378871265e-10
        },
        "0.5": {
          "lambda": 2.5737381413782287e-11,
          "Delta": 6.988147390288639,
          "floor": 4.40643861356427,
          "plateau": 11.39458600385291,
          "z1_times_a": 7.529134813217906,
          "n_zeros": 381,
          "lock_metric": 7.612992319993416e-05,
          "residual": 10.708784969160092,
          "mass_pm0.1": 127.66507880386828,
          "solve_resid": 2.2679633796708461e-07
        },
        "0.9": {
          "lambda": 4.632728654480812e-11,
          "Delta": 6.987816911690756,
          "floor": 4.407352992396391,
          "plateau": 11.395169904087147,
          "z1_times_a": 8.199973337356482,
          "n_zeros": 381,
          "lock_metric": 7.609529464864893e-05,
          "residual": 10.693566680716287,
          "mass_pm0.1": 127.66305186152309,
          "solve_resid": 5.74001354100028e-08
        },
        "0.99": {
          "lambda": 5.096001519928893e-11,
          "Delta": 6.985923725722728,
          "floor": 4.413501470282933,
          "plateau": 11.39942519600566,
          "z1_times_a": 8.385279853786248,
          "n_zeros": 381,
          "lock_metric": 7.613585558613195e-05,
          "residual": 10.687862236134379,
          "mass_pm0.1": 127.66303788362372,
          "solve_resid": 7.885784848149343e-08
        },
        "0.999": {
          "lambda": 5.142328806473701e-11,
          "Delta": 6.981356415160433,
          "floor": 4.428696114419206,
          "plateau": 11.41005252957964,
          "z1_times_a": 8.40265274248596,
          "n_zeros": 381,
          "lock_metric": 7.614180518407809e-05,
          "residual": 10.687276136059936,
          "mass_pm0.1": 127.66306304014373,
          "solve_resid": 4.3024270231059637e-07
        }
      }
    },
    "8": {
      "lambda_a": -0.028701446327309428,
      "lambda_law_ratio": -15509.255243545844,
      "sister_corr": 0.00255772313467671,
      "sister_crossings": 118,
      "sister_cross_err": 18.633851388497806,
      "ladder": {
        "-1": {
          "lambda": -1.0,
          "Delta": 2.7352201581877043,
          "floor": 0.030554703235580436,
          "plateau": 2.7657748614232847,
          "z1_times_a": 4.496810102477515,
          "n_zeros": 763,
          "lock_metric": 0.057829797607063305,
          "residual": 58.651149793527864,
          "mass_pm0.1": 133.51693635274552,
          "solve_resid": 3.0889933849343203e-16
        },
        "0": {
          "lambda": -0.0,
          "Delta": -33.528034311582985,
          "floor": -423.8211920422485,
          "plateau": -457.3492263538315,
          "z1_times_a": 1.8093924166647792,
          "n_zeros": 548,
          "lock_metric": 0.0003714172395450807,
          "residual": 11.144052463580522,
          "mass_pm0.1": 127.82233169936302,
          "solve_resid": 0.00020060867988680158
        },
        "0.5": {
          "lambda": -0.014350723163654714,
          "Delta": 5.2473962964853165,
          "floor": 1.7336013999547513,
          "plateau": 6.980997696440068,
          "z1_times_a": 4.498199381863492,
          "n_zeros": 755,
          "lock_metric": 0.0013619017100949234,
          "residual": 12.052748449542477,
          "mass_pm0.1": 127.6918385778458,
          "solve_resid": 6.867617904061254e-13
        },
        "0.9": {
          "lambda": -0.025831301694578486,
          "Delta": 5.070963483079262,
          "floor": 1.333969497055245,
          "plateau": 6.4049329801345065,
          "z1_times_a": 4.497946694779682,
          "n_zeros": 761,
          "lock_metric": 0.0030658627341504153,
          "residual": 13.166729243282433,
          "mass_pm0.1": 127.88295027859101,
          "solve_resid": 1.2366600020318642e-13
        },
        "0.99": {
          "lambda": -0.028414431864036332,
          "Delta": 5.037993177808347,
          "floor": 1.272757105974604,
          "plateau": 6.3107502837829506,
          "z1_times_a": 4.5002419528288815,
          "n_zeros": 761,
          "lock_metric": 0.0018016901302855674,
          "residual": 12.566089279144094,
          "mass_pm0.1": 127.72782983544735,
          "solve_resid": 2.3651958714205305e-14
        },
        "0.999": {
          "lambda": -0.028672744880982117,
          "Delta": 5.034476634521893,
          "floor": 1.2679657032940865,
          "plateau": 6.30244233781598,
          "z1_times_a": 4.495276414486893,
          "n_zeros": 762,
          "lock_metric": 0.0026471749147504637,
          "residual": 12.710218865379431,
          "mass_pm0.1": 127.70661909585449,
          "solve_resid": 1.413266823421088e-14
        }
      }
    },
    "12": {
      "lambda_a": -0.0682789599667037,
      "lambda_law_ratio": -27121745.659492664,
      "sister_corr": -0.014045368106618201,
      "sister_crossings": 137,
      "sister_cross_err": 18.87885138849781,
      "ladder": {
        "-1": {
          "lambda": -1.0,
          "Delta": 2.9442229398281423,
          "floor": 0.22834844487190686,
          "plateau": 3.172571384700049,
          "z1_times_a": 4.495304974166888,
          "n_zeros": 1145,
          "lock_metric": 0.032106115361490564,
          "residual": 57.59425795470206,
          "mass_pm0.1": 133.31147006008132,
          "solve_resid": 5.565591099499078e-16
        },
        "0": {
          "lambda": -0.0,
          "Delta": 145.87233394077197,
          "floor": -369.5573776739587,
          "plateau": -223.68504373318672,
          "z1_times_a": 47.08983339916813,
          "n_zeros": 788,
          "lock_metric": 0.11988236920414796,
          "residual": 10.704038663232183,
          "mass_pm0.1": 126.78597572681957,
          "solve_resid": 0.00035380049050859374
        },
        "0.5": {
          "lambda": -0.03413947998335185,
          "Delta": 4.517541168837466,
          "floor": 2.012961451365246,
          "plateau": 6.530502620202712,
          "z1_times_a": 4.495602611843745,
          "n_zeros": 1127,
          "lock_metric": 0.0016805564393911263,
          "residual": 12.576655387504829,
          "mass_pm0.1": 127.80026124326591,
          "solve_resid": 2.361326934047527e-10
        },
        "0.9": {
          "lambda": -0.06145106397003333,
          "Delta": 4.373132705161369,
          "floor": 1.5778486637685551,
          "plateau": 5.9509813689299245,
          "z1_times_a": 4.4955360310205545,
          "n_zeros": 1141,
          "lock_metric": 0.002769088717484152,
          "residual": 14.752940107360544,
          "mass_pm0.1": 128.0191898119376,
          "solve_resid": 3.98328418223779e-14
        },
        "0.99": {
          "lambda": -0.06759617036703666,
          "Delta": 4.345367791047291,
          "floor": 1.511240245787921,
          "plateau": 5.856608036835212,
          "z1_times_a": 4.495571310035841,
          "n_zeros": 1143,
          "lock_metric": 0.0030508525265727083,
          "residual": 15.875746711965263,
          "mass_pm0.1": 128.06515694909066,
          "solve_resid": 2.4872731506462405e-14
        },
        "0.999": {
          "lambda": -0.068210681006737,
          "Delta": 4.342626921736709,
          "floor": 1.505160225931791,
          "plateau": 5.847787147668501,
          "z1_times_a": 4.495984857538869,
          "n_zeros": 1144,
          "lock_metric": 0.003716572174764971,
          "residual": 15.417387806155247,
          "mass_pm0.1": 128.0996751837004,
          "solve_resid": 6.800740851599709e-15
        }
      }
    }
  },
  "elapsed_seconds": 307.2,
  "timestamp_utc": "2026-08-19T18:39:23.796254+00:00"
}

"""E-025 — task 2.13: lambda-flow numerics (D-042 §3, Suzuki's lever).

At fixed lambda = -1 the family D_{a,theta} equilibrates at the
consolidation boundary (linear-K law, D-042). The remaining channel:
lambda -> lambda_a^- (bottom of the Weil-form spectrum, lambda_a ~
e^{-1.65a}), where T_a(lambda) = Q - lambda M becomes near-singular
and the ground state v_a dominates — the regime of Suzuki's sister
conjecture (1.2) (c_a v^_a -> xi(1/2+iz) = Xi(z)).

Mechanism being tested: as lambda -> lambda_a^-, the deficiency
vectors v_pm(lambda) ~ v_a/(lambda_a - lambda), so the zeros of
W(a,theta;.) flow toward the zeros of v^_a; if v^_a ~ Xi, those are
the gamma's => locking.

Per a in {4, 8, 12}:
  A. lambda_a and ground state v_a: smallest generalized eigenpair of
     the pencil (Q, M) restricted to the mean-zero subspace (QR
     projection, dense eigh; n <= 2400).
  B. sister-conjecture observable: v^_a(z) = hat_FT * sum e^{izx} v_a
     on z-grid; correlation of Re v^_a with Xi(z) (flint xi), zero
     crossings vs the first gammas.
  C. lambda ladder: lambda/lambda_a in {-1(abs), 0, 0.5, 0.9, 0.99,
     0.999}:  K_lambda(z,z) on z-grid in [14,60] (valley depth
     Delta(lambda), floor, plateau); zeros of W (theta=pi) ->
     LOCKING metric: mean over first 20 gammas of dist(gamma, nearest
     W-zero)/(pi/2a) (1.0 = chance level); residual mass R(lambda).

Checks:
  V1. lambda_a positive and consistent with the e^{-1.65a} law
      (order of magnitude; D-017/D-019).
  V2. lambda=-1 baseline reproduces E-022: z1*a in [4.3,4.7],
      Delta within 0.05 of E-022b valley depth.
  V3. Levinson residuals sane at all ladder rungs
      (|T v - b|/|b| < 1e-6 except possibly the 0.999 rung, reported).
  V4. sister observable reported (correlation v^_a vs Xi).

Outputs: research/numerics/data/E-025-lambda-flow.csv (+ .json)
"""
import importlib.util
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import flint
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-025-lambda-flow.csv"
JSON = DATA / "E-025-lambda-flow.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
LAM0 = E16.LAM_SHIFT          # -1.0 baseline
A_LIST = [4.0, 8.0, 12.0]
LADDER_REL = [0.0, 0.5, 0.9, 0.99, 0.999]   # x lambda_a (plus lam=-1)
Z_MAX = 300.0
Z_METRIC = 280.0
N_GRID = 40000
XN = E16.XN
mp.mp.dps = 40


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 xi_real(z):
    """Xi(z) = xi(1/2 - i z) via flint (real on R)."""
    flint.ctx.prec = 120
    out = np.empty(len(z))
    for i, zz in enumerate(z):
        s = flint.acb(0.5, -float(zz))
        val = (0.5 * s * (s - 1) * flint.acb.pi() ** (-s / 2)
               * (s / 2).gamma() * s.zeta())
        out[i] = float(val.real)
    return out


def locate_zeros_pi(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 = E16.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)


def ground_state(Q, M):
    """Smallest generalized eigenpair of (Q, M) on the mean-zero
    subspace (E-009 route: QR projection)."""
    n = Q.shape[0]
    Z = np.linalg.qr(np.ones((n, 1)), mode="complete")[0][:, 1:]
    Qp = Z.T @ Q @ Z
    Mp = Z.T @ M @ Z
    evals, evecs = scipy.linalg.eigh(Qp, Mp)
    lam_a = float(evals[0])
    v = Z @ evecs[:, 0]
    return lam_a, v


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

    print("backend: load E-022 cache ...", flush=True)
    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 = np.array([float(mp.zetazero(k).imag) for k in range(1, 201)])
    gam = gam[gam <= Z_METRIC]

    zgrid = np.arange(14.0, 60.0, 0.05)
    gsel = gam[(gam >= 14) & (gam < 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)
    XI = xi_real(zgrid)

    rows = []
    report = {"per_a": {}}
    v1_ok = v2_ok = v3_ok = True
    e022b_dip = {4.0: 2.1593, 8.0: 2.7352, 12.0: 2.9442}

    for a in A_LIST:
        ta = time.time()
        x, h, Q, M, b = E16.build(a, H)
        n = len(x)

        # A: ground state
        lam_a, v_a = ground_state(Q, M)
        if not (lam_a > 0):
            v1_ok = False
        law = np.exp(-1.65 * a)
        print(f"a={a:g}: lambda_a = {lam_a:.3e} (e^-1.65a = {law:.1e})",
              flush=True)

        # B: sister-conjecture observable
        E = np.exp(1j * np.outer(zgrid, x))
        vhat = E16.hat_FT(zgrid, h) * (E @ v_a)
        rv = vhat.real
        ca = float(np.dot(rv, XI) / np.dot(rv, rv))
        corr = float(np.corrcoef(rv, XI)[0, 1])
        # zero crossings of Re v^_a on grid
        sc = np.flatnonzero(np.sign(rv[:-1]) * np.sign(rv[1:]) < 0)
        cross = zgrid[sc]
        ncr = min(10, len(cross), len(gsel))
        xerr = float(np.mean(np.abs(cross[:ncr] - gsel[:ncr]))) \
            if ncr else float("nan")
        print(f"  sister: corr(v^,Xi) = {corr:+.4f}, {len(cross)} "
              f"crossings, mean |cross-gamma| (first {ncr}) = {xerr:.3f}",
              flush=True)

        per = {"lambda_a": lam_a, "lambda_law_ratio": lam_a / law,
               "sister_corr": corr, "sister_crossings": int(len(cross)),
               "sister_cross_err": xerr, "ladder": {}}

        # C: lambda ladder
        for tag, lam in ([("-1", LAM0)]
                         + [(f"{r:g}", r * lam_a) for r in LADDER_REL]):
            tc = Q[:, 0] - lam * (np.r_[2 * h / 3, h / 6,
                                        np.zeros(max(n - 2, 0))])
            u1 = scipy.linalg.solve_toeplitz((tc, tc), np.ones(n))
            s1 = u1.sum()
            T_dense = None  # built lazily for residual spot-checks

            def K1(z, _tc=tc, _u1=u1, _s1=s1, _chk=False):
                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
                kv = float(np.real(np.conj(bz) @ v))
                if _chk:
                    Tm = scipy.linalg.toeplitz(_tc)
                    rs = float(np.max(np.abs(
                        Tm @ v + (u0.sum() / _s1) * np.ones(n) - bz))
                        / max(1.0, float(np.max(np.abs(bz)))))
                    return kv, rs
                return kv, 0.0

            lg = np.empty(len(zgrid))
            resmax = 0.0
            for i, zz in enumerate(zgrid):
                kv, rs = K1(zz, _chk=(i % 100 == 0))
                lg[i] = np.log(max(kv, 1e-300))
                resmax = max(resmax, rs)
            if resmax > 1e-6:
                v3_ok = False
            dip = float(lg[idxm].mean() - lg[idxg].mean())
            floor = float(lg[idxg].mean())
            plate = float(lg[idxm].mean())

            # zeros of W (theta=pi) at this lambda (Schur form,
            # validated 3.6e-14 in E-022)
            u0p = scipy.linalg.solve_toeplitz((tc, tc), b[+1.0])
            cp = u0p - (u0p.sum() / s1) * u1
            zeros = locate_zeros_pi(cp, x, h)
            zm = zeros[zeros <= Z_METRIC]
            # locking: mean normalized distance of first 20 gammas
            d20 = []
            for g in gam[:20]:
                d20.append(np.min(np.abs(zm - g)))
            lock = float(np.mean(d20) / (np.pi / (2 * a)))
            # weights & residual
            ws = np.array([1.0 / K1(z)[0] for z in zm])
            near_idx = np.array([np.argmin(np.abs(zm - g)) for g in gam])
            resid = float(ws.sum() - ws[near_idx].sum())
            m01 = float(ws[np.minimum.reduce([
                np.abs(zm - g) for g in gam]).reshape(-1) <= 0.1].sum()) \
                if len(zm) else 0.0

            per["ladder"][tag] = {
                "lambda": lam, "Delta": dip, "floor": floor,
                "plateau": plate, "z1_times_a": float(zeros[0] * a),
                "n_zeros": int(len(zeros)),
                "lock_metric": lock, "residual": resid,
                "mass_pm0.1": m01, "solve_resid": resmax,
            }
            print(f"  lam={tag}*lam_a ({lam:.3e}): Delta={dip:.3f}, "
                  f"z1a={zeros[0] * a:.3f}, lock={lock:.3f}, "
                  f"resid={resid:.1f}, K-resid={resmax:.1e}", flush=True)
            for zz, ll in zip(zgrid, lg):
                rows.append((a, tag, zz, ll))

        # V2: baseline checks (lambda = -1)
        base = per["ladder"]["-1"]
        if not (4.3 < base["z1_times_a"] < 4.7):
            v2_ok = False
        if abs(base["Delta"] - e022b_dip[a]) > 0.05:
            v2_ok = False
        report["per_a"][f"{a:g}"] = per
        print(f"a={a:g} done ({time.time() - ta:.0f}s)", flush=True)

    checks = {"V1_lambda_a_positive_law": v1_ok,
              "V2_baseline_matches_E022": v2_ok,
              "V3_solve_residuals": v3_ok,
              "V4_sister_reported": True}

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

    meta = {
        "experiment": "E-025",
        "task": "2.13: lambda-flow numerics (ground-state channel)",
        "script": "research/numerics/scripts/E-025-lambda-flow.py",
        "parameters": {"a_list": A_LIST, "ladder_rel": LADDER_REL,
                       "z_max": Z_MAX},
        "checks": checks,
        "key_values": report["per_a"],
        "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())