E-016c-v3-repair.py

{
  "experiment": "E-016c",
  "task": "2.6: V3 repair — stable rescaled analytic rho_W + pair-augmented Poisson reconstruction",
  "script": "research/numerics/scripts/E-016c-v3-repair.py",
  "parameters": {
    "a_sweep": [
      4.0,
      10.0
    ],
    "h": 0.01,
    "lambda": -1.0,
    "y_vals": [
      0.5,
      1.0,
      2.0,
      5.0,
      10.0,
      15.0,
      20.0,
      30.0,
      40.0
    ],
    "z_max": 300.0
  },
  "checks": {
    "V5_rescaled_stable": true,
    "V3R_repaired": true
  },
  "report": {
    "4.0": {
      "n_zeros": 381,
      "eps": 1e-12,
      "V5_rescaled_vs_mp_max_rel": 1.746533135566571e-16,
      "V6_direct_vs_mp_rel_by_y": {
        "0.5": 4.361403255037962e-15,
        "1.0": 6.459307529272192e-16,
        "2.0": 2.747798737159212e-15,
        "5.0": 1.394742189728899e-15,
        "10.0": 1.0478266069041646e-15,
        "15.0": 1.746533135566571e-16,
        "20.0": 1.7464113580770303e-16,
        "30.0": 3.492287178806053e-16,
        "40.0": 0.0
      },
      "V3R_shape_rms_rel": 1.9891458745035763e-05,
      "resid_after_repair": {
        "0.5": -1.0312442324522664e-06,
        "1.0": -2.062471045727321e-06,
        "2.0": -4.124802745142375e-06,
        "5.0": -1.030956893210977e-05,
        "10.0": -2.0601741559600484e-05,
        "15.0": -3.085921220113974e-05,
        "20.0": -4.106485540478211e-05,
        "30.0": -6.125358190423036e-05,
        "40.0": -8.103777361379016e-05
      }
    },
    "10.0": {
      "n_zeros": 954,
      "eps": 1e-12,
      "V5_rescaled_vs_mp_max_rel": 1.3957277747606402e-16,
      "V6_direct_vs_mp_rel_by_y": {
        "0.5": 0.0,
        "1.0": 0.0,
        "2.0": 4.1768373220595484e-16,
        "5.0": 8.372113810976072e-16,
        "10.0": 6.978540770655109e-16,
        "15.0": 0.0,
        "20.0": 8.374304062150007e-16,
        "30.0": 1.3956955567962724e-16,
        "40.0": 0.0
      },
      "V3R_shape_rms_rel": 1.4618285263569542e-05,
      "resid_after_repair": {
        "0.5": -1.8063902298415258e-06,
        "1.0": -3.6127501874538837e-06,
        "2.0": -7.225258211285279e-06,
        "5.0": -1.8058908774598592e-05,
        "10.0": -3.608758578899085e-05,
        "15.0": -5.4055959257670594e-05,
        "20.0": -7.193427555707288e-05,
        "30.0": -0.00010730423592120175,
        "40.0": -0.00014197153716599686
      }
    }
  },
  "elapsed_seconds": 11.5,
  "timestamp_utc": "2026-08-18T17:33:41.931287+00:00"
}

"""E-016c — V3 repair: stable rescaled analytic rho_W + pair-augmented
reconstruction. Closes task 2.6 if V3 passes with the fixed machinery.

E-016b root cause (D-025): the V3 failure had TWO independent defects:
  (1) Galerkin W(a,pi;.) has a double zero at z ~ 0 (pair +/- i eps,
      eps ~ 2e-9..6e-9; winding excess exactly 2 at every radius). The
      real-axis brentq recon misses it -> residual 2/(pi*y) at small y.
  (2) float64 catastrophic cancellation in the ANALYTIC rho_W at large
      y: V(iy) sums span e^{+/- y(a-h)} of dynamic range (e^80 at a=4,
      y=20; e^200 at a=10). rho_analytic loses the entire zero tail
      (notail recon resid at y=20: +0.09 (a=4), +0.17 (a=10) vs the
      ~1.08 lattice tail the EM formula correctly expects).

Fix: rescale V(z) = e^{-i(a-h) z} Vtilde(z) with all exponents <= 0
(well-conditioned); W'/W = -i(a-h) + bracket'/bracket;
-Im W'/W(iy) = (a-h) - Im(bracket'/bracket)(iy).
Augment recon with the pair: + (2/pi) y/(y^2 - eps^2).

Checks:
  V5. rescaled float64 vs mpmath dps-50 rho_W: max rel dev < 1e-9.
  V6. direct float64 vs mpmath: degradation documented (expected to
      FAIL at large y — this is the measured defect).
  V3R. repaired V3: rescaled analytic vs recon+pair, mean-matched
      shape rms / rms(rho) < 0.02 (the original V3 threshold).

Outputs: research/numerics/data/E-016c-v3-repair.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.optimize

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-016c-v3-repair.csv"
JSON = DATA / "E-016c-v3-repair.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)

A_SWEEP = [4.0, 10.0]
H = E16.H_MAIN
LAM = E16.LAM_SHIFT
Y_VALS = [0.5, 1.0, 2.0, 5.0, 10.0, 15.0, 20.0, 30.0, 40.0]
Z_MAX = 300.0
N_GRID = 40000


def rho_direct(ziy, cp, cm, x, h):
    """Original E-016 float64 evaluation."""
    Vp, dVp = E16.V_and_dV(ziy, cp, x, h)
    Vm, dVm = E16.V_and_dV(ziy, cm, x, h)
    W, dW = E16.W_and_dW(ziy, Vp, dVp, Vm, dVm, np.pi)
    return -np.imag(dW / W) / np.pi


def rho_rescaled(ziy, cp, cm, x, h):
    """Stable evaluation for z = i y, y > 0.
    V(z) = e^{i xm z} Vt(z) with xm = min x node (u = x - xm >= 0, so
    e^{i z u} = e^{-y u} <= 1). W = e^{i xm z} B, B = (z-i)Vtp - (z+i)Vtm
    => W'/W = i xm + B'/B => -Im W'/W(iy) = -xm - Im(B'/B)(iy)."""
    xm = x.min()
    u = x - xm                                     # >= 0
    E = np.exp(1j * np.outer(ziy, u))              # |e^{i z u}| = e^{-y u}
    hat = E16.hat_FT(ziy, h)
    hatd = E16.hat_FT_d(ziy, h)
    S1p, S2p = E @ cp, (1j * u * E) @ cp
    S1m, S2m = E @ cm, (1j * u * E) @ cm
    # dVt = hat'*S1 + hat*S2  (Vt' w.r.t. z); dVtp below adds + i xm Vt
    # so that (dVtp - i xm Vtp) recovers Vt'.
    Vtp, dVtp = hat * S1p, hatd * S1p + hat * (S2p + 1j * xm * S1p)
    Vtm, dVtm = hat * S1m, hatd * S1m + hat * (S2m + 1j * xm * S1m)
    B = (ziy - 1j) * Vtp - (ziy + 1j) * Vtm
    dB = Vtp + (ziy - 1j) * (dVtp - 1j * xm * Vtp) \
        - (Vtm + (ziy + 1j) * (dVtm - 1j * xm * Vtm))
    return (-xm - np.imag(dB / B)) / np.pi


def rho_mpmath(y, cp, cm, x, h, dps=50):
    """Reference: mpmath dps-50 evaluation of -Im W'/W(iy)/pi."""
    mp.mp.dps = dps
    cpm = [mp.mpf(float(v)) for v in cp]
    cmm = [mp.mpf(float(v)) for v in cm]
    xm = [mp.mpf(float(v)) for v in x]
    hh = mp.mpf(h)
    yy = mp.mpf(y)
    # work with the rescaled (well-conditioned) sums directly
    x0 = xm[0]                                      # = min x
    um = [v - x0 for v in xm]                       # >= 0

    def Vt(c):
        S = sum(ci * mp.e ** (-yy * ui) for ci, ui in zip(c, um))
        u_hat = yy * hh / 2
        hat = hh * (mp.sinh(u_hat) / u_hat) ** 2    # hat(iy) real
        return hat * S

    def dVt(c):
        # d/dz [hat(z) S(z)] at z = i y, S(z) = sum c e^{i z u}
        u_hat = yy * hh / 2
        hat = hh * (mp.sinh(u_hat) / u_hat) ** 2
        # hat'(z) at z = i y: d/dz hat = i * d/dy hat(iy)
        # d/dy [hh (sinh(u)/u)^2], u = y hh/2
        d_hat_dy = hh**2 / 2 * 2 * (mp.sinh(u_hat) / u_hat) * (
            (u_hat * mp.cosh(u_hat) - mp.sinh(u_hat)) / u_hat**2)
        hatd = -1j * d_hat_dy      # d/dz f(iy) = -i * d/dy f(iy)
        S = sum(ci * mp.e ** (-yy * ui) for ci, ui in zip(c, um))
        dS = sum(1j * ui * ci * mp.e ** (-yy * ui)
                 for ci, ui in zip(c, um))
        return hatd * S + hat * dS

    Vtp, dVtp = Vt(cpm), dVt(cpm)
    Vtm, dVtm = Vt(cmm), dVt(cmm)
    z = 1j * yy
    xmn = float(x0)
    B = (z - 1j) * Vtp - (z + 1j) * Vtm
    dB = Vtp + (z - 1j) * dVtp - (Vtm + (z + 1j) * dVtm)  # dVt = Vt' already
    return float((-xmn - mp.im(dB / B)) / mp.pi)


def locate_zeros_and_pair(cp, cm, x, h, a):
    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
    # pair at +/- i eps: root of w(y) = Im W(iy) near 0 (mpmath)
    mp.mp.dps = 50
    cpm = [mp.mpf(float(v)) for v in cp]
    cmm = [mp.mpf(float(v)) for v in cm]
    xmm = [mp.mpf(float(v)) for v in x]
    hh = mp.mpf(h)

    def w(yy):
        yy = mp.mpf(yy)
        Sp = sum(c * mp.e ** (-yy * xx) for c, xx in zip(cpm, xmm))
        Sm = sum(c * mp.e ** (-yy * xx) for c, xx in zip(cmm, xmm))
        u_hat = yy * hh / 2
        hat = hh * (mp.sinh(u_hat) / u_hat) ** 2 if yy != 0 else hh
        return mp.re((yy - 1) * hat * Sp - (yy + 1) * hat * Sm)

    ys_scan = [mp.mpf(10) ** k for k in range(-12, -1)]
    fs = [w(yy) for yy in ys_scan]
    eps = None
    for i in range(len(ys_scan) - 1):
        if fs[i] * fs[i + 1] < 0:
            lo, hi = ys_scan[i], ys_scan[i + 1]
            flo = fs[i]
            for _ in range(200):
                mid = (lo + hi) / 2
                fm = w(mid)
                if flo * fm <= 0:
                    hi = mid
                else:
                    lo, flo = mid, fm
            eps = float((lo + hi) / 2)
            break
    if eps is None:                     # fallback: smallest |w|
        eps = float(ys_scan[int(np.argmin([abs(f) for f in fs]))])
    return np.array(zeros), eps


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

    for a in A_SWEEP:
        print(f"a={a}: build + solve + zeros ...", flush=True)
        x, hh, Q, M, b = E16.build(a, H)
        n = len(x)
        cp, T = E16.solve_v(Q, M, b[+1.0], LAM, n)
        cm, _ = E16.solve_v(Q, M, b[-1.0], LAM, n)
        np_ = float(cp @ T @ cp)
        nm = float(cm @ T @ cm)
        cm = cm * np.sqrt(np_ / nm)
        zeros, eps = locate_zeros_and_pair(cp, cm, x, hh, a)
        d_lat = np.pi / a
        z_tail = zeros[-1] + d_lat / 2.0
        print(f"  zeros={len(zeros)} eps={eps:.3e}", flush=True)

        yg = np.array(Y_VALS)
        r_dir = rho_direct(1j * yg, cp, cm, x, hh)
        r_res = rho_rescaled(1j * yg, cp, cm, x, hh)
        r_mp = np.array([rho_mpmath(y, cp, cm, x, hh) for y in yg])

        # recon + pair (Poisson). NOTE: the EM tail is
        # (1/d)∫ y/(y²+z²) dz = (1/d)(π/2 − arctan(z0/y)) — NO factor y
        # (E-016/E-016b carried a spurious y — the y=20 deficit source).
        r_rec = np.empty_like(yg)
        for k, y in enumerate(yg):
            s = 2.0 * np.sum(y / (y**2 + zeros**2))
            s += 2.0 * (1.0 / d_lat) * (np.pi / 2.0 - np.arctan(z_tail / y))
            s += 2.0 * y / (y**2 - eps**2)
            r_rec[k] = s / np.pi

        dev_res_mp = np.abs(r_res - r_mp) / np.maximum(np.abs(r_mp), 1e-12)
        dev_dir_mp = np.abs(r_dir - r_mp) / np.maximum(np.abs(r_mp), 1e-12)
        res = r_mp - r_rec
        res_shift = res - np.mean(res)
        v3r = float(np.sqrt(np.mean(res_shift**2))
                    / np.sqrt(np.mean(r_mp**2)))

        for k, y in enumerate(yg):
            rows.append(f"{a},{y:.2f},{r_dir[k]:.10e},{r_res[k]:.10e},"
                        f"{r_mp[k]:.10e},{r_rec[k]:.10e},"
                        f"{dev_dir_mp[k]:.3e},{dev_res_mp[k]:.3e}")
        report[str(a)] = {
            "n_zeros": int(len(zeros)), "eps": eps,
            "V5_rescaled_vs_mp_max_rel": float(dev_res_mp.max()),
            "V6_direct_vs_mp_rel_by_y": {str(y): float(d) for y, d in
                                         zip(yg, dev_dir_mp)},
            "V3R_shape_rms_rel": v3r,
            "resid_after_repair": {str(y): float(r) for y, r in
                                   zip(yg, res)},
        }
        print(f"  V5 rescaled-vs-mp max rel: {dev_res_mp.max():.2e}",
              flush=True)
        print(f"  V6 direct-vs-mp rel by y: "
              f"{ {y: f'{d:.1e}' for y, d in zip(yg, dev_dir_mp)} }",
              flush=True)
        print(f"  V3R shape rms rel: {v3r:.3e}", flush=True)

    checks = {
        "V5_rescaled_stable": bool(all(
            report[str(a)]["V5_rescaled_vs_mp_max_rel"] < 1e-9
            for a in A_SWEEP)),
        "V3R_repaired": bool(all(
            report[str(a)]["V3R_shape_rms_rel"] < 0.02 for a in A_SWEEP)),
    }
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("a,y,rho_direct,rho_rescaled,rho_mpmath,rho_recon_pair,"
                "rel_direct_vs_mp,rel_rescaled_vs_mp\n")
        f.write("\n".join(rows) + "\n")
    meta = {
        "experiment": "E-016c",
        "task": "2.6: V3 repair — stable rescaled analytic rho_W + "
                "pair-augmented Poisson reconstruction",
        "script": "research/numerics/scripts/E-016c-v3-repair.py",
        "parameters": {"a_sweep": A_SWEEP, "h": H, "lambda": LAM,
                       "y_vals": Y_VALS, "z_max": Z_MAX},
        "checks": checks,
        "report": report,
        "elapsed_seconds": round(time.time() - t0, 1),
        "timestamp_utc": datetime.now(timezone.utc).isoformat(),
    }
    JSON.write_text(json.dumps(meta, indent=2), encoding="utf-8")
    ok = all(checks.values())
    print(f"wrote {CSV.name}, {JSON.name} — "
          f"{'ALL CHECKS PASS' if ok else 'CHECK FAILURE: ' + str(checks)}")
    return 0 if ok else 1


if __name__ == "__main__":
    raise SystemExit(main())