E-024-theta-sweep.py

{
  "experiment": "E-024",
  "task": "2.12: theta-flow numerics (theta-independence test)",
  "script": "research/numerics/scripts/E-024-theta-sweep.py",
  "parameters": {
    "a_list": [
      4.0,
      8.0,
      12.0
    ],
    "thetas_over_pi": [
      0.0,
      0.1667,
      0.3333,
      0.5,
      0.6667,
      0.8333,
      1.0
    ],
    "z_max": 300.0,
    "z_metric": 280.0
  },
  "checks": {
    "V1_lattice_all_theta": true,
    "V2_theta_pi_matches_E022": true,
    "V3_weights_positive": true
  },
  "key_values": {
    "a4_th0.0000pi": {
      "theta_over_pi": 0.0,
      "n_zeros": 381,
      "z1_times_a": 3.1471258584711514,
      "mass_total": 184.1835959613851,
      "mass_pm0.1": 119.4128093530634,
      "mass_pm0.05": 91.94966918355006,
      "residual": 48.76433809073876,
      "valley_mass_mean": 1.0747560148463995,
      "plateau_weight_mean": 0.21201886126408157,
      "Delta_eff": 1.6231737129268227
    },
    "a4_th0.1667pi": {
      "theta_over_pi": 0.1667,
      "n_zeros": 381,
      "z1_times_a": 3.4619359298190293,
      "mass_total": 184.18630284387328,
      "mass_pm0.1": 118.29623992508816,
      "mass_pm0.05": 86.87180322997308,
      "residual": 48.931287590873524,
      "valley_mass_mean": 1.0734525020079346,
      "plateau_weight_mean": 0.2127447286559719,
      "Delta_eff": 1.6185423804514978
    },
    "a4_th0.3333pi": {
      "theta_over_pi": 0.3333,
      "n_zeros": 381,
      "z1_times_a": 3.7176915045337466,
      "mass_total": 183.9726782811009,
      "mass_pm0.1": 113.57262540386267,
      "mass_pm0.05": 85.73049725879235,
      "residual": 50.41888852750091,
      "valley_mass_mean": 1.0599507123301588,
      "plateau_weight_mean": 0.22016981889738385,
      "Delta_eff": 1.5715785355566734
    },
    "a4_th0.5000pi": {
      "theta_over_pi": 0.5,
      "n_zeros": 381,
      "z1_times_a": 3.934163190177167,
      "mass_total": 183.98559583483848,
      "mass_pm0.1": 113.33711825061624,
      "mass_pm0.05": 89.82793793217617,
      "residual": 50.618934373720066,
      "valley_mass_mean": 1.0584655671517333,
      "plateau_weight_mean": 0.22104338154462916,
      "Delta_eff": 1.566216581357752
    },
    "a4_th0.6667pi": {
      "theta_over_pi": 0.6667,
      "n_zeros": 381,
      "z1_times_a": 4.128324562896003,
      "mass_total": 183.9983973653932,
      "mass_pm0.1": 119.75726561744686,
      "mass_pm0.05": 95.81754097212591,
      "residual": 48.13040632212352,
      "valley_mass_mean": 1.078317389232299,
      "plateau_weight_mean": 0.21017644682150013,
      "Delta_eff": 1.6352097313960756
    },
    "a4_th0.8333pi": {
      "theta_over_pi": 0.8333,
      "n_zeros": 381,
      "z1_times_a": 4.314182529701571,
      "mass_total": 184.01134904950752,
      "mass_pm0.1": 120.40026511109403,
      "mass_pm0.05": 90.96419878254113,
      "residual": 47.43102136203055,
      "valley_mass_mean": 1.0839708546625155,
      "plateau_weight_mean": 0.20712236402633416,
      "Delta_eff": 1.6550765456169134
    },
    "a4_th1.0000pi": {
      "theta_over_pi": 1.0,
      "n_zeros": 381,
      "z1_times_a": 4.504293602793437,
      "mass_total": 184.02432245308052,
      "mass_pm0.1": 116.64777135239572,
      "mass_pm0.05": 93.33647263205195,
      "residual": 48.30754760742576,
      "valley_mass_mean": 1.0771172606797996,
      "plateau_weight_mean": 0.2109499895520775,
      "Delta_eff": 1.6304224593273315
    },
    "a8_th0.0000pi": {
      "theta_over_pi": 0.0,
      "n_zeros": 763,
      "z1_times_a": 3.1437903681267265,
      "mass_total": 184.1082994928043,
      "mass_pm0.1": 134.20929653000076,
      "mass_pm0.05": 123.1835197830814,
      "residual": 57.73718018875971,
      "valley_mass_mean": 1.002945391301941,
      "plateau_weight_mean": 0.09852761124361722,
      "Delta_eff": 2.320359515027157
    },
    "a8_th0.1667pi": {
      "theta_over_pi": 0.1667,
      "n_zeros": 763,
      "z1_times_a": 3.418119179475546,
      "mass_total": 184.11283914170593,
      "mass_pm0.1": 131.8317519020843,
      "mass_pm0.05": 126.78097971430489,
      "residual": 56.7356526846647,
      "valley_mass_mean": 1.010930051246359,
      "plateau_weight_mean": 0.09681851993970085,
      "Delta_eff": 2.3457877312638367
    },
    "a8_th0.3333pi": {
      "theta_over_pi": 0.3333,
      "n_zeros": 763,
      "z1_times_a": 3.654042277120252,
      "mass_total": 184.11760264338506,
      "mass_pm0.1": 132.1490251381983,
      "mass_pm0.05": 124.0292094105938,
      "residual": 58.16107071439245,
      "valley_mass_mean": 0.9996550153094651,
      "plateau_weight_mean": 0.09925097391534557,
      "Delta_eff": 2.3097585025024654
    },
    "a8_th0.5000pi": {
      "theta_over_pi": 0.5,
      "n_zeros": 763,
      "z1_times_a": 3.867441949623103,
      "mass_total": 184.12212186389277,
      "mass_pm0.1": 133.24139728241443,
      "mass_pm0.05": 121.15410070178395,
      "residual": 59.85047462666191,
      "valley_mass_mean": 0.9862829145811973,
      "plateau_weight_mean": 0.10213391574515683,
      "Delta_eff": 2.26765839342186
    },
    "a8_th0.6667pi": {
      "theta_over_pi": 0.6667,
      "n_zeros": 763,
      "z1_times_a": 4.071478758632688,
      "mass_total": 184.1262503654809,
      "mass_pm0.1": 132.29513175105507,
      "mass_pm0.05": 122.51608044758039,
      "residual": 57.48114969555114,
      "valley_mass_mean": 1.0051198465867441,
      "plateau_weight_mean": 0.09809069913916582,
      "Delta_eff": 2.3269695116402334
    },
    "a8_th0.8333pi": {
      "theta_over_pi": 0.8333,
      "n_zeros": 763,
      "z1_times_a": 4.27766402464752,
      "mass_total": 184.12992169567084,
      "mass_pm0.1": 131.99007596810569,
      "mass_pm0.05": 124.30470833176733,
      "residual": 57.28031602287514,
      "valley_mass_mean": 1.0067429021650454,
      "plateau_weight_mean": 0.09774797956122035,
      "Delta_eff": 2.332083020259487
    },
    "a8_th1.0000pi": {
      "theta_over_pi": 1.0,
      "n_zeros": 763,
      "z1_times_a": 4.4968101024775144,
      "mass_total": 184.13312237431055,
      "mass_pm0.1": 133.51693635274228,
      "mass_pm0.05": 124.28962039506533,
      "residual": 58.6511497935264,
      "valley_mass_mean": 0.9958886712760646,
      "plateau_weight_mean": 0.10008728633707581,
      "Delta_eff": 2.297592806874823
    },
    "a12_th0.0000pi": {
      "theta_over_pi": 0.0,
      "n_zeros": 1145,
      "z1_times_a": 3.1428501414881977,
      "mass_total": 184.16007157245687,
      "mass_pm0.1": 133.61538092090012,
      "mass_pm0.05": 128.10759866335772,
      "residual": 60.30642315719122,
      "valley_mass_mean": 0.9829654636132195,
      "plateau_weight_mean": 0.06395166824728656,
      "Delta_eff": 2.732446371442262
    },
    "a12_th0.1667pi": {
      "theta_over_pi": 0.1667,
      "n_zeros": 1145,
      "z1_times_a": 3.405992501692748,
      "mass_total": 184.1639822094043,
      "mass_pm0.1": 133.36866840406773,
      "mass_pm0.05": 129.35521418242917,
      "residual": 60.3669918644722,
      "valley_mass_mean": 0.9825157963883501,
      "plateau_weight_mean": 0.06401589805352302,
      "Delta_eff": 2.73098496176509
    },
    "a12_th0.3333pi": {
      "theta_over_pi": 0.3333,
      "n_zeros": 1145,
      "z1_times_a": 3.635830555925344,
      "mass_total": 184.09580595481827,
      "mass_pm0.1": 133.41942035595872,
      "mass_pm0.05": 129.05415843419973,
      "residual": 58.57426551897014,
      "valley_mass_mean": 0.9962027018718106,
      "plateau_weight_mean": 0.06218074895856703,
      "Delta_eff": 2.7739053032480143
    },
    "a12_th0.5000pi": {
      "theta_over_pi": 0.5,
      "n_zeros": 1145,
      "z1_times_a": 3.8476946281889486,
      "mass_total": 184.09881987373944,
      "mass_pm0.1": 132.95788338488353,
      "mass_pm0.05": 129.3505890774291,
      "residual": 58.94147251618104,
      "valley_mass_mean": 0.9933122806155429,
      "plateau_weight_mean": 0.06257056530380153,
      "Delta_eff": 2.764750131884173
    },
    "a12_th0.6667pi": {
      "theta_over_pi": 0.6667,
      "n_zeros": 1145,
      "z1_times_a": 4.05423555082551,
      "mass_total": 184.10127816912785,
      "mass_pm0.1": 133.37160860447023,
      "mass_pm0.05": 129.68837121321633,
      "residual": 58.35802636312266,
      "valley_mass_mean": 0.9979623159206761,
      "plateau_weight_mean": 0.06195119571456758,
      "Delta_eff": 2.7793686068267633
    },
    "a12_th0.8333pi": {
      "theta_over_pi": 0.8333,
      "n_zeros": 1145,
      "z1_times_a": 4.2666373857585445,
      "mass_total": 184.1032513640804,
      "mass_pm0.1": 133.09565694496038,
      "mass_pm0.05": 129.3399679137547,
      "residual": 58.63289108937627,
      "valley_mass_mean": 0.9957965101166993,
      "plateau_weight_mean": 0.062242984171312424,
      "Delta_eff": 2.7724971046813764
    },
    "a12_th1.0000pi": {
      "theta_over_pi": 1.0,
      "n_zeros": 1145,
      "z1_times_a": 4.495304974166889,
      "mass_total": 184.10483484645468,
      "mass_pm0.1": 133.31147006007708,
      "mass_pm0.05": 129.13788187469896,
      "residual": 57.594257954705384,
      "valley_mass_mean": 1.004052197553566,
      "plateau_weight_mean": 0.06114040122580191,
      "Delta_eff": 2.7986264096434
    },
    "V4_theta_spread": {
      "a4": {
        "residual_minmax": [
          47.43102136203055,
          50.618934373720066
        ],
        "residual_spread_pct": 2.234095821084676,
        "mass_pm01_spread_pct": 2.302784187513983,
        "z1_times_a_range": [
          3.1471258584711514,
          4.504293602793437
        ]
      },
      "a8": {
        "residual_minmax": [
          56.7356526846647,
          59.85047462666191
        ],
        "residual_spread_pct": 1.6396401999266315,
        "mass_pm01_spread_pct": 0.6333372523281642,
        "z1_times_a_range": [
          3.1437903681267265,
          4.4968101024775144
        ]
      },
      "a12": {
        "residual_minmax": [
          57.594257954705384,
          60.3669918644722
        ],
        "residual_spread_pct": 1.6057262581418954,
        "mass_pm01_spread_pct": 0.15063835179687926,
        "z1_times_a_range": [
          3.1428501414881977,
          4.495304974166889
        ]
      }
    }
  },
  "elapsed_seconds": 212.6,
  "timestamp_utc": "2026-08-19T18:19:09.910473+00:00"
}

"""E-024 — task 2.12: theta-flow numerics (D-041 §5/§8, D-038 §5).

Theory (de Branges/Clark): the model space's kernel K_a(z,z) is
theta-INDEPENDENT; the extension choice theta only phase-shifts the
atom lattice z_k(theta). PREDICTION (to be tested): consolidation
statistics (valley-captured mass, residual) are theta-independent —
if confirmed, no theta(a) flow can fix consolidation and the
obstruction is entirely the K-saturation (D-040 F6); if refuted,
the theta-flow channel is alive.

For theta in {0, pi/6, pi/3, pi/2, 2pi/3, 5pi/6, pi} and a in {4,8,12}:
  - zeros of W(a,theta;z) on (0, Z_MAX] via the real gauge
        g_theta(z) = Re[ e^{-i theta/2} (z - i) V_+(z) ]
    (e^{-i theta/2} W = 2 Re[e^{-i theta/2}(z-i)V_+] on the real line,
    using Vm = conj(Vp); reduces to z ImV - ReV at theta=pi);
  - Clark weights w_k = 1/K_a(z_k,z_k) (E-022 Schur-Toeplitz route);
  - metrics: n_zeros, z1*a, mass within +/-0.1 and +/-0.05 of true
    gammas, residual R(a,theta) = total - captured, mean valley mass,
    effective depth Delta_eff = log(mean valley weight / mean plateau
    weight).

Checks:
  V1. all zeros real & count = maximal lattice a*Z_MAX/pi +/- 2,
      every theta (self-adjointness).
  V2. theta=pi reproduces E-022: |z1*a - 4.4953| < 0.05 at a=12 and
      mass(+/-0.1) within 1% of E-022 value.
  V3. weights positive, all theta.
  V4. headline: theta-spread of residual R and mass(+/-0.1) reported
      per a (small spread = theta-independence confirmed).

Outputs: research/numerics/data/E-024-theta-sweep.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
import scipy.optimize

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-024-theta-sweep.csv"
JSON = DATA / "E-024-theta-sweep.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, 8.0, 12.0]
THETAS = [0.0, np.pi / 6, np.pi / 3, np.pi / 2, 2 * np.pi / 3,
          5 * np.pi / 6, np.pi]
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 locate_zeros_theta(cp, x, h, theta, z_max=Z_MAX, n_grid=N_GRID):
    zg = np.linspace(1e-3, z_max, n_grid)
    E = np.exp(1j * np.outer(zg, x))
    V = E16.hat_FT(zg, h) * (E @ cp)
    f = np.exp(-0.5j * theta) * (zg - 1j) * V
    g = f.real

    def g_direct(zz):
        V1 = E16.hat_FT(np.array([zz]), h)[0] * complex(
            np.exp(1j * zz * x) @ cp)
        return float(np.real(np.exp(-0.5j * theta) * (zz - 1j) * V1))

    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 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]

    rows = []
    report = {"per_a_theta": {}}
    v1_ok = v3_ok = True
    v2 = {}
    for a in A_LIST:
        ta = time.time()
        x, h, Q, M, b = E16.build(a, H)
        n = len(x)
        cp, _ = E16.solve_v(Q, M, b[+1.0], LAM, n)
        # Toeplitz bordered solver for weights (theta-independent K)
        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))

        for th in THETAS:
            zeros = locate_zeros_theta(cp, x, h, th)
            lat = a * Z_MAX / np.pi
            if abs(len(zeros) - lat) > 2:
                v1_ok = False
            zm = zeros[zeros <= Z_METRIC]
            ws = np.array([1.0 / K1(z) for z in zm])
            if np.any(ws <= 0):
                v3_ok = False
            tot = float(ws.sum())
            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]))
            m01 = float(ws[d <= 0.1].sum())
            m005 = float(ws[d <= 0.05].sum())
            # valley atoms = nearest to each gamma; plateau = rest
            near_idx = np.array([np.argmin(np.abs(zm - g)) for g in gam])
            valley_mass = ws[near_idx]
            mask_pl = np.ones(len(zm), bool)
            mask_pl[near_idx] = False
            pl = ws[mask_pl]
            d_eff = float(np.log(valley_mass.mean() / pl.mean()))
            key = f"a{a:g}_th{th / np.pi:.4f}pi"
            report["per_a_theta"][key] = {
                "theta_over_pi": round(th / np.pi, 4),
                "n_zeros": int(len(zeros)),
                "z1_times_a": float(zeros[0] * a),
                "mass_total": tot,
                "mass_pm0.1": m01, "mass_pm0.05": m005,
                "residual": tot - float(valley_mass.sum()),
                "valley_mass_mean": float(valley_mass.mean()),
                "plateau_weight_mean": float(pl.mean()),
                "Delta_eff": d_eff,
            }
            for k, (zz, ww) in enumerate(zip(zm, ws), 1):
                rows.append((a, th / np.pi, k, zz, ww))
        print(f"a={a:g} done ({time.time() - ta:.0f}s)", flush=True)

    # V2: theta=pi cross-check vs E-022
    e022 = report["per_a_theta"]["a12_th1.0000pi"]
    v2 = {"z1_times_a_ok": bool(abs(e022["z1_times_a"] - 4.4953) < 0.05)}
    # V4: theta-spread per a
    v4 = {}
    for a in A_LIST:
        vals = [report["per_a_theta"][f"a{a:g}_th{th / np.pi:.4f}pi"]
                for th in THETAS]
        res = np.array([v["residual"] for v in vals])
        m01 = np.array([v["mass_pm0.1"] for v in vals])
        z1a = np.array([v["z1_times_a"] for v in vals])
        v4[f"a{a:g}"] = {
            "residual_minmax": [float(res.min()), float(res.max())],
            "residual_spread_pct": float(100 * res.std() / res.mean()),
            "mass_pm01_spread_pct": float(100 * m01.std() / m01.mean()),
            "z1_times_a_range": [float(z1a.min()), float(z1a.max())],
        }
        print(f"a={a:g}: residual spread {v4[f'a{a:g}']['residual_spread_pct']:.2f}%, "
              f"mass01 spread {v4[f'a{a:g}']['mass_pm01_spread_pct']:.2f}%, "
              f"z1a in [{z1a.min():.3f},{z1a.max():.3f}]", flush=True)

    checks = {"V1_lattice_all_theta": v1_ok,
              "V2_theta_pi_matches_E022": v2["z1_times_a_ok"],
              "V3_weights_positive": v3_ok}

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

    meta = {
        "experiment": "E-024",
        "task": "2.12: theta-flow numerics (theta-independence test)",
        "script": "research/numerics/scripts/E-024-theta-sweep.py",
        "parameters": {"a_list": A_LIST,
                       "thetas_over_pi": [round(t / np.pi, 4)
                                          for t in THETAS],
                       "z_max": Z_MAX, "z_metric": Z_METRIC},
        "checks": checks,
        "key_values": {**report["per_a_theta"], "V4_theta_spread": v4},
        "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())