E-026-inverse-spectral.py

{
  "experiment": "E-026",
  "task": "inverse-spectral reconstruction (Jacobi/Stieltjes coefficients of the Clark measures)",
  "script": "research/numerics/scripts/E-026-inverse-spectral.py",
  "parameters": {
    "a_list": [
      4.0,
      8.0,
      12.0
    ],
    "k_max": 240,
    "dps": 50
  },
  "checks": {
    "V1_roundtrip": false,
    "V2_symmetry": true,
    "V3_free_reference_documented": true
  },
  "key_values": {
    "4": {
      "n_atoms": 762,
      "b_first10": [
        197.53030725422184,
        142.92095536812724,
        159.15960898964033,
        144.02245566814906,
        156.54864075079917,
        145.58592069411515,
        155.01008099518884,
        146.33185226710341,
        154.10605961168793,
        146.84721248193748
      ],
      "b_tail_mean_last20": 142.99220059266153,
      "b_ref_tail_mean_last20": 143.37936432309326,
      "b_min": 136.4184875527748,
      "b_max": 197.53030725422184,
      "roundtrip_dev": 99.95090939634419
    },
    "8": {
      "n_atoms": 1526,
      "b_first10": [
        197.4687711112365,
        142.93932543312803,
        159.09877103726024,
        144.04943192616722,
        156.47961609870868,
        145.62114275673025,
        154.93543142819607,
        146.37506545353162,
        154.02642581802172,
        146.89798821598438
      ],
      "b_tail_mean_last20": 148.12846034636544,
      "b_ref_tail_mean_last20": 148.34613679083813,
      "b_min": 141.78183290120748,
      "b_max": 197.4687711112365,
      "roundtrip_dev": 102.09127232819776
    },
    "12": {
      "n_atoms": 2290,
      "b_first10": [
        197.4516510445814,
        142.94721881560383,
        159.08047821579584,
        144.06020454619352,
        156.45866720384134,
        145.6346431789002,
        154.91256454534792,
        146.3911624786654,
        154.00186642701698,
        146.9165327438543
      ],
      "b_tail_mean_last20": 149.2050889664049,
      "b_ref_tail_mean_last20": 149.18259058194445,
      "b_min": 141.8742264108275,
      "b_max": 197.4516510445814,
      "roundtrip_dev": 102.80351310637383
    }
  },
  "elapsed_seconds": 957.7,
  "timestamp_utc": "2026-08-19T20:21:10.797721+00:00"
}

"""E-026 — inverse-spectral reconstruction (sub-stages A+B).

From the E-022 Clark measures (atoms z_k = W(a,pi;.) zeros, masses
w_k = 1/K_a(z_k,z_k)) reconstruct the Jacobi/Stieltjes recurrence
coefficients of the symmetric spectral measure
    mu_a = sum_k w_k (delta_{z_k} + delta_{-z_k})  (+ central handling).
The coefficient sequence {b_k} (a_k = 0 by symmetry) is the discrete
shadow of the canonical-system Hamiltonian (Krein-string data) — the
inverse-spectral observable behind the spectral-purity question
(D-001 sweep entry, 2026-08-19).

Method: Stieltjes/Lanczos on the discrete inner product
    <p,q> = sum_atoms w p q,
with full reorthogonalization, mpmath dps 50. The measure is symmetric
about 0, so a_k = 0 identically (asserted numerically).

Checks:
  V1. Round-trip: eigendecomposition of the K x K Jacobi matrix
      reproduces the first K atoms/masses of the input measure
      (Gauss-quadrature exactness: a K-point Jacobi matrix matches
      the 2K moments; compare spectral atoms to the input's innermost
      atoms within tolerance, weights within 1%).
  V2. Symmetry: |a_k| < 1e-30 (dps-50 zero).
  V3. Free-reference: for the equal-mass lattice measure (same atoms,
      unit masses), b_k settle to a constant (documented).

Outputs: research/numerics/data/E-026-jacobi.csv (+ .json)
Columns: a, k, b_k (measured) ; b_k for the equal-mass reference.
"""
import csv
import json
import sys
import time
from datetime import datetime, timezone
from pathlib import Path

import mpmath as mp
import numpy as np

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

A_LIST = [4.0, 8.0, 12.0]
K_MAX = 240          # recurrence coefficients per side
DPS = 50


def load_measure(a):
    z, w = [], []
    with open(DATA / "E-022-spectral-measure.csv") as f:
        for r in csv.DictReader(f):
            if float(r["a"]) == a:
                z.append(float(r["z_k"]))
                w.append(float(r["w_k"]))
    z = np.array(z)
    w = np.array(w)
    # symmetric measure on +- atoms
    zz = np.concatenate([-z[::-1], z])
    ww = np.concatenate([w[::-1], w]) / 2.0
    return zz, ww


def stieltjes(z, w, k_max, dps=DPS):
    """Lanczos/Stieltjes with full reorthogonalization at dps."""
    mp.mp.dps = dps
    n = len(z)
    zm = [mp.mpf(float(v)) for v in z]
    wm = [mp.mpf(float(v)) for v in w]
    tot = sum(wm)
    # p_0 = 1/sqrt(tot): work with normalized vectors u_k on the atoms
    # (u_k[i] = p_k(z_i) * sqrt(w_i)) — the L2(w) picture.
    sq = [mp.sqrt(v) for v in wm]
    U = []                       # orthonormal vectors (mp lists)
    u0 = [sq[i] / mp.sqrt(tot) for i in range(n)]
    U.append(u0)
    bs = []
    a_dev_max = mp.mpf(0)
    for k in range(k_max):
        # q = z u_k - b_k u_{k-1}   (a_k = <z u_k, u_k> ~ 0 by symmetry)
        uk = U[k]
        a_k = sum(zm[i] * uk[i] ** 2 for i in range(n))
        a_dev_max = max(a_dev_max, abs(a_k))
        q = [zm[i] * uk[i] - a_k * uk[i] for i in range(n)]
        if k > 0:
            q = [q[i] - bs[-1] * U[k - 1][i] for i in range(n)]
        # full reorthogonalization
        for j in range(k + 1):
            c = sum(q[i] * U[j][i] for i in range(n))
            if c:
                q = [q[i] - c * U[j][i] for i in range(n)]
        nq = mp.sqrt(sum(v * v for v in q))
        bs.append(nq)
        U.append([v / nq for v in q])
    return bs, float(a_dev_max)


def roundtrip_check(z, w, bs, k_use=40):
    """Jacobi matrix from bs[:k_use-1] -> Gauss atoms/weights vs input's
    innermost atoms."""
    k = k_use
    J = np.diag(np.array([float(b) for b in bs[:k - 1]]), 1)
    J = J + J.T
    ev, evec = np.linalg.eigh(J)
    wts = evec[0, :] ** 2 * float(np.sum(w))
    # compare: sorted evs near 0 vs input atoms near 0
    zin = np.sort(z)
    inner_in = zin[len(zin) // 2 - 5: len(zin) // 2 + 5]
    inner_ev = ev[len(ev) // 2 - 5: len(ev) // 2 + 5]
    dev = float(np.max(np.abs(np.sort(inner_in) - np.sort(inner_ev))))
    return dev


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()
    rows = []
    report = {"per_a": {}}
    rt_max = 0.0
    asym_max = 0.0

    for a in A_LIST:
        ta = time.time()
        z, w = load_measure(a)
        print(f"a={a:g}: {len(z)} atoms; Stieltjes to {K_MAX} ...",
              flush=True)
        bs, adev = stieltjes(z, w, K_MAX)
        asym_max = max(asym_max, adev)
        # equal-mass reference on the same atoms
        wref = np.ones_like(w)
        bs_ref, _ = stieltjes(z, wref, K_MAX)
        dev = roundtrip_check(z, w, bs)
        rt_max = max(rt_max, dev)
        bs_f = [float(b) for b in bs]
        bsr_f = [float(b) for b in bs_ref]
        report["per_a"][f"{a:g}"] = {
            "n_atoms": int(len(z)),
            "b_first10": bs_f[:10],
            "b_tail_mean_last20": float(np.mean(bs_f[-20:])),
            "b_ref_tail_mean_last20": float(np.mean(bsr_f[-20:])),
            "b_min": float(min(bs_f)), "b_max": float(max(bs_f)),
            "roundtrip_dev": dev,
        }
        for k, (b, br) in enumerate(zip(bs_f, bsr_f), 1):
            rows.append((a, k, b, br))
        print(f"  b_1..5 = {np.round(bs_f[:5], 5)}; "
              f"tail mean {np.mean(bs_f[-20:]):.5f} "
              f"(ref {np.mean(bsr_f[-20:]):.5f}); roundtrip {dev:.2e} "
              f"({time.time() - ta:.0f}s)", flush=True)

    checks = {
        "V1_roundtrip": bool(rt_max < 0.05),
        "V2_symmetry": bool(asym_max < 1e-30),
        "V3_free_reference_documented": True,
    }
    with open(CSV, "w", newline="") as f:
        wcsv = csv.writer(f)
        wcsv.writerow(["a", "k", "b_k", "b_k_ref_equalmass"])
        for a, k, b, br in rows:
            wcsv.writerow([f"{a:g}", k, f"{b:.12f}", f"{br:.12f}"])
    meta = {
        "experiment": "E-026",
        "task": "inverse-spectral reconstruction (Jacobi/Stieltjes "
                "coefficients of the Clark measures)",
        "script": "research/numerics/scripts/E-026-inverse-spectral.py",
        "parameters": {"a_list": A_LIST, "k_max": K_MAX, "dps": DPS},
        "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())