E-016b-v3-rootcause.py

{
  "experiment": "E-016b",
  "task": "2.6: root-cause E-016 V3 Poisson-reconstruction failure",
  "script": "research/numerics/scripts/E-016b-v3-rootcause.py",
  "parameters": {
    "a_sweep": [
      4.0,
      10.0
    ],
    "h": 0.01,
    "lambda": -1.0,
    "y_grid": [
      0.1,
      40.0,
      400
    ],
    "z_max": 300.0
  },
  "report": {
    "4.0": {
      "T1_kappa_re": -0.9914379801328772,
      "T1_kappa_im": -1.4222222222222223,
      "T2_census": {
        "n_real": 381,
        "z1": 1.1260734006983593,
        "z_last": 299.8646747694667,
        "lattice_count_window": 381.9718634205488
      },
      "T2_winding": {
        "0.3927": {
          "winding": 1.9971283856655768,
          "real_zeros_inside": 0,
          "excess": 1.9971283856655768
        },
        "1.1781": {
          "winding": 3.957366770074305,
          "real_zeros_inside": 2,
          "excess": 1.9573667700743052
        },
        "0.0000": {
          "winding": 1.9960937479061762,
          "real_zeros_inside": 0,
          "excess": 1.9960937479061762
        },
        "0.0001": {
          "winding": 1.9960937500549862,
          "real_zeros_inside": 0,
          "excess": 1.9960937500549862
        },
        "0.0100": {
          "winding": 1.996094371929568,
          "real_zeros_inside": 0,
          "excess": 1.996094371929568
        },
        "0.1000": {
          "winding": 1.9961562249900062,
          "real_zeros_inside": 0,
          "excess": 1.9961562249900062
        },
        "2.0000": {
          "winding": 5.942871677547558,
          "real_zeros_inside": 4,
          "excess": 1.9428716775475579
        },
        "5.0000": {
          "winding": 12.077973882052314,
          "real_zeros_inside": 10,
          "excess": 2.077973882052314
        }
      },
      "T2_imag_pair": {
        "eps": 5.829287623947368e-09,
        "W_at_ieps": 2.825742938081395e-49,
        "scan_min": 5.749824935036285e-06,
        "scan_at_y": 0.0001
      },
      "T3_A_over_y_fit": 0.6343007263567857,
      "T3_two_over_pi": 0.6366197723675814,
      "T3_mean_shift": -1.2815872994129114,
      "T4_variants": {
        "base": {
          "resid_y0.5": 1.2739134084231456,
          "resid_y5": 0.07332699024962142,
          "resid_y20": -0.9925381870041659,
          "shape_rms_rel": 1.0497908515365875
        },
        "trunc60": {
          "resid_y0.5": 1.2766229564752873,
          "resid_y5": -0.14349117160523672,
          "resid_y20": -4.94373487936198,
          "shape_rms_rel": 4.338262071341213
        },
        "trunc100": {
          "resid_y0.5": 1.27526491761311,
          "resid_y5": -0.034968760716001634,
          "resid_y20": -3.0133337980338,
          "shape_rms_rel": 2.809762823982386
        },
        "notail": {
          "resid_y0.5": 1.2745883033553684,
          "resid_y5": 0.14081030856832366,
          "resid_y20": 0.08570192959282363,
          "shape_rms_rel": 0.2895380312354026
        },
        "base+pair": {
          "resid_y0.5": 0.0006738636879826121,
          "resid_y5": -0.05399696422389488,
          "resid_y20": -1.0243691756225448,
          "shape_rms_rel": 0.9421400227183758
        },
        "mp60": {
          "resid_y0.5": 1.276622956475289,
          "resid_y5": -0.14349117160523361,
          "resid_y20": -4.943734879361979,
          "shape_rms_rel": 4.338262071341214
        }
      },
      "T4_mp60": {
        "n_relocated": 75,
        "max_abs_dz": 1.2434497875801753e-12
      },
      "T5_sinc_est_max": 0.0002055923399439411
    },
    "10.0": {
      "T1_kappa_re": -3.6447054985051435,
      "T1_kappa_im": -4.397114393679149,
      "T2_census": {
        "n_real": 954,
        "z1": 0.44959700551214055,
        "z_last": 299.89357740917427,
        "lattice_count_window": 954.9296585513721
      },
      "T2_winding": {
        "0.1571": {
          "winding": 1.9971325652242644,
          "real_zeros_inside": 0,
          "excess": 1.9971325652242644
        },
        "0.4712": {
          "winding": 3.9589810973742314,
          "real_zeros_inside": 2,
          "excess": 1.9589810973742314
        },
        "0.0000": {
          "winding": 1.996093695657393,
          "real_zeros_inside": 0,
          "excess": 1.996093695657393
        },
        "0.0001": {
          "winding": 1.9960937503824123,
          "real_zeros_inside": 0,
          "excess": 1.9960937503824123
        },
        "0.0100": {
          "winding": 1.9960976525143628,
          "real_zeros_inside": 0,
          "excess": 1.9960976525143628
        },
        "0.1000": {
          "winding": 1.9964955397146036,
          "real_zeros_inside": 0,
          "excess": 1.9964955397146036
        },
        "2.0000": {
          "winding": 12.085419712205663,
          "real_zeros_inside": 10,
          "excess": 2.0854197122056632
        },
        "5.0000": {
          "winding": 31.975335381644364,
          "real_zeros_inside": 30,
          "excess": 1.9753353816443635
        }
      },
      "T2_imag_pair": {
        "eps": 2.0797966440236883e-09,
        "W_at_ieps": 1.6504059507421108e-26,
        "scan_min": 0.01459451355707131,
        "scan_at_y": 0.0001
      },
      "T3_A_over_y_fit": 0.6308206254641555,
      "T3_two_over_pi": 0.6366197723675814,
      "T3_mean_shift": -3.3631535473854246,
      "T4_variants": {
        "base": {
          "resid_y0.5": 1.2749261380781172,
          "resid_y5": -0.007753707743303639,
          "resid_y20": -2.530820269504034,
          "shape_rms_rel": 1.0143675095987894
        },
        "trunc60": {
          "resid_y0.5": 1.2816709492693397,
          "resid_y5": -0.5465597395886301,
          "resid_y20": -12.354286478739445,
          "shape_rms_rel": 4.43884573453234
        },
        "trunc100": {
          "resid_y0.5": 1.278305028486991,
          "resid_y5": -0.27789168396584607,
          "resid_y20": -7.57271682340175,
          "shape_rms_rel": 2.860573136578675
        },
        "notail": {
          "resid_y0.5": 1.2766145378113016,
          "resid_y5": 0.16107079639175792,
          "resid_y20": 0.16663165033476446,
          "shape_rms_rel": 0.11754357982145858
        },
        "base+pair": {
          "resid_y0.5": 0.0016865933429546232,
          "resid_y5": -0.13507766221681994,
          "resid_y20": -2.5626512581224135,
          "shape_rms_rel": 0.9793208614044097
        },
        "mp60": {
          "resid_y0.5": 1.2816709492693397,
          "resid_y5": -0.5465597395886301,
          "resid_y20": -12.354286478739445,
          "shape_rms_rel": 4.438845734532339
        }
      },
      "T4_mp60": {
        "n_relocated": 190,
        "max_abs_dz": 1.0587086762825493e-12
      },
      "T5_sinc_est_max": 0.0002055923399439411
    }
  },
  "elapsed_seconds": 702.8,
  "timestamp_utc": "2026-08-18T17:19:13.123417+00:00"
}

"""E-016b — root-cause the E-016 V3 (zero-Poisson reconstruction) failure.

E-016 V3: rho_W^analytic (Galerkin z-derivative) vs rho_W^recon (located
real zeros + Euler-Maclaurin lattice tail) mismatch rel ~0.3 at a=4.
Residual (D-024 sec.4): +1.274 (y=0.5), +0.637 (1), +0.313 (2), +0.073
(5), -0.179 (10), -0.993 (20). Note 2/(pi*y) = 1.2732, 0.6366, 0.3183 —
the low-y residual matches a MISSING zero pair at |z|~0 to 4 digits.

Tests (D-024 sec.4 candidates + new hypothesis):
  T1. kappa = W'/W(0): exact gauge arithmetic gives kappa = 0 by symmetry
      (x symmetric, c real). Numeric confirmation. => candidate (a) dead.
  T2. Zero census + ARGUMENT PRINCIPLE on circles |z|=r: the real-gauge
      brentq scan can only find REAL zeros. If Galerkin W has a complex
      pair near 0 (candidate: purely imaginary +/-i*eps, self-conjugate
      under the mirror symmetries), the winding number exceeds the real-
      zero count. If found: locate eps by scanning |W(iy)|.
  T3. Residual shape on a fine y-grid: test A/y law; report A vs 2/pi.
  T4. Recon variants: (i) base (zeros<=300 + EM tail); (ii) trunc60;
      (iii) trunc100; (iv) notail; (v) mp60 (zeros<=60 relocated at
      dps 50 + EM tail from 60 — float64 contamination test, D-024 (c));
      (vi) base+pair (recon + exact contribution of the T2 imaginary
      pair, (2/pi) y/(y^2-eps^2) — decisive if T2 fires).
  T5. Sinc-artifact zeros of hat h at z=2*pi*k/h (=628k, beyond window):
      analytic size of their rho contribution (expected ~1e-4, dismiss).

Outputs: research/numerics/data/E-016b-v3-rootcause.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-016b-v3-rootcause.csv"
JSON = DATA / "E-016b-v3-rootcause.json"

# ---- import E-016 machinery (runs the sieve once, ~15 s) ---------------
_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_GRID = np.linspace(0.1, 40.0, 400)
Z_MAX = 300.0
N_GRID = 40000
WIND_RADII_FRAC = [0.5, 1.5]   # in units of pi/a, plus absolute below
WIND_RADII_ABS = [1e-6, 1e-4, 1e-2, 0.1, 2.0, 5.0]


def locate_real_zeros(cp, x, h, z_max=Z_MAX, n_grid=N_GRID):
    """Real-axis sign-change + brentq on g(z) = z Im Vp - Re Vp (E-016)."""
    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 winding(cp, cm, x, h, r, n_pts=512):
    """Argument-principle zero count of W(theta=pi) on |z|=r."""
    phi = np.linspace(0, 2 * np.pi, n_pts, endpoint=False)
    z = r * np.exp(1j * phi)
    Vp, dVp = E16.V_and_dV(z, cp, x, h)
    Vm, dVm = E16.V_and_dV(z, cm, x, h)
    W, _ = E16.W_and_dW(z, Vp, dVp, Vm, dVm, np.pi)
    dphi = np.diff(np.unwrap(np.angle(W)))
    return float(dphi.sum() / (2 * np.pi))


def imag_pair_scan(cp, cm, x, h, y_max=3.0):
    """Scan |W(iy)| for a minimum -> candidate purely-imaginary pair."""
    ys = np.linspace(1e-4, y_max, 30000)
    z = 1j * ys
    Vp, _ = E16.V_and_dV(z, cp, x, h)
    Vm, _ = E16.V_and_dV(z, cm, x, h)
    W = (z - 1j) * Vp - (z + 1j) * Vm          # theta = pi
    aw = np.abs(W)
    i = int(np.argmin(aw))
    # polish with mpmath findroot on W(iy)=0 (real unknown y)
    mp.mp.dps = 50
    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]

    def Wimag(yy):
        yy = mp.mpf(yy)
        Sp = sum(c * mp.e ** (1j * yy * xx) for c, xx in zip(cpm, xm))
        Sm = sum(c * mp.e ** (1j * yy * xx) for c, xx in zip(cmm, xm))
        hh = mp.mpf(h)
        hat = hh * (mp.sin(yy * hh / 2) / (yy * hh / 2)) ** 2
        Vpv, Vmv = hat * Sp, hat * Sm
        W = 1j * (yy - 1) * Vpv - 1j * (yy + 1) * Vmv      # W(i*yy), th=pi
        return mp.im(W)          # mirror symmetry: W(iy) = i * (real)

    try:
        root = mp.findroot(Wimag, mp.mpf(float(ys[i])),
                           solver="secant", tol=mp.mpf("1e-40"))
        eps = float(root)
        resid = float(abs(Wimag(root)))
        return eps, resid, float(aw[i]), float(ys[i])
    except (ValueError, ZeroDivisionError):
        return None, None, float(aw[i]), float(ys[i])


def em_tail(y, z_start, d_lat):
    return (y / d_lat) * (np.pi / 2.0 - np.arctan(z_start / y))


def recon_rho(y, zeros_pos, d_lat, z_tail_start=None, extra_pairs=()):
    """Poisson density from positive zeros (doubled) + EM tail + explicit
    extra pairs (purely imaginary +/-i*eps: 2y/(y^2 - eps^2))."""
    y = np.asarray(y, dtype=float)
    zeros_pos = np.asarray(zeros_pos, dtype=float)
    if zeros_pos.size:
        s = 2.0 * np.sum(y[:, None] / (y[:, None] ** 2
                                       + zeros_pos[None, :] ** 2), axis=1)
    else:
        s = np.zeros_like(y)
    if z_tail_start is not None:
        s = s + 2.0 * em_tail(y, z_tail_start, d_lat)
    for eps in extra_pairs:
        s = s + 2.0 * y / (y**2 - eps**2)
    return s / np.pi


def mp_relocate(cp, x, h, zeros_pos, z_hi=60.0):
    """Re-locate zeros <= z_hi with mpmath dps 50 (float64 contamination
    test). Bisection on the real gauge, 80 iterations."""
    mp.mp.dps = 50
    cpm = [mp.mpf(float(v)) for v in cp]
    xm = [mp.mpf(float(v)) for v in x]
    hh = mp.mpf(h)

    def g_mp(zz):
        zz = mp.mpf(zz)
        S = sum(c * mp.e ** (1j * zz * xx) for c, xx in zip(cpm, xm))
        hat = hh * (mp.sin(zz * hh / 2) / (zz * hh / 2)) ** 2
        V = hat * S
        return zz * mp.im(V) - mp.re(V)

    out = []
    for z0 in zeros_pos:
        if z0 > z_hi:
            continue
        lo, hi = mp.mpf(float(z0)) - mp.mpf("0.05"), \
            mp.mpf(float(z0)) + mp.mpf("0.05")
        flo, fhi = g_mp(lo), g_mp(hi)
        if flo * fhi > 0:
            out.append((float(z0), None))
            continue
        for _ in range(80):
            mid = (lo + hi) / 2
            fm = g_mp(mid)
            if flo * fm <= 0:
                hi, fhi = mid, fm
            else:
                lo, flo = mid, fm
        out.append((float(z0), float((lo + hi) / 2)))
    return out


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 ...", 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)
        d_lat = np.pi / a

        # analytic rho_W at theta=pi on the fine grid
        ziy = 1j * Y_GRID
        Vp, dVp = E16.V_and_dV(ziy, cp, x, hh)
        Vm, dVm = E16.V_and_dV(ziy, cm, x, hh)
        W, dW = E16.W_and_dW(ziy, Vp, dVp, Vm, dVm, np.pi)
        rho_an = -np.imag(dW / W) / np.pi

        # T1: kappa = W'/W(0)
        z0 = np.array([0.0 + 0.0j])
        Vp0, dVp0 = E16.V_and_dV(z0, cp, x, hh)
        Vm0, dVm0 = E16.V_and_dV(z0, cm, x, hh)
        W0, dW0 = E16.W_and_dW(z0, Vp0, dVp0, Vm0, dVm0, np.pi)
        kappa = complex(dW0[0] / W0[0])
        print(f"  T1 kappa = {kappa:.3e}", flush=True)

        # T2: census + argument principle
        zeros = locate_real_zeros(cp, x, hh)
        census = {"n_real": int(len(zeros)), "z1": float(zeros[0]),
                  "z_last": float(zeros[-1]),
                  "lattice_count_window": a * Z_MAX / np.pi}
        radii = [f * d_lat for f in WIND_RADII_FRAC] + WIND_RADII_ABS
        wind = {}
        for r in radii:
            w = winding(cp, cm, x, hh, r)
            n_real_in = int(2 * np.sum(zeros < r))
            wind[f"{r:.4f}"] = {"winding": w, "real_zeros_inside": n_real_in,
                                "excess": w - n_real_in}
        print(f"  T2 winding: {wind}", flush=True)

        eps, eps_resid, scan_min, scan_at = None, None, None, None
        if any(abs(v["excess"]) > 0.5 for v in wind.values()):
            eps, eps_resid, scan_min, scan_at = imag_pair_scan(cp, cm, x, hh)
            print(f"  T2 imaginary pair: eps={eps} |W(i*eps)|={eps_resid:.2e} "
                  f"(scan min {scan_min:.3e} at y={scan_at:.4f})", flush=True)

        # T3: residual shape (base variant), A/y fit at small y
        z_tail = zeros[-1] + d_lat / 2.0
        rec_base = recon_rho(Y_GRID, zeros, d_lat, z_tail)
        resid = rho_an - rec_base
        sel_small = (Y_GRID >= 0.3) & (Y_GRID <= 2.0)
        A_fit = float(np.mean(resid[sel_small] * Y_GRID[sel_small]))
        # mean-shift removed version for shape comparison
        shift = float(np.mean(resid))

        # T4: variants
        z60 = zeros[zeros <= 60.0]
        z100 = zeros[zeros <= 100.0]
        variants = {
            "base": rec_base,
            "trunc60": recon_rho(Y_GRID, z60, d_lat, z60[-1] + d_lat / 2),
            "trunc100": recon_rho(Y_GRID, z100, d_lat, z100[-1] + d_lat / 2),
            "notail": recon_rho(Y_GRID, zeros, d_lat, None),
        }
        if eps is not None and eps_resid is not None and eps_resid < 1e-20:
            variants["base+pair"] = recon_rho(
                Y_GRID, zeros, d_lat, z_tail, extra_pairs=(eps,))
        print("  T4 mp60 relocation (dps 50) ...", flush=True)
        reloc = mp_relocate(cp, x, hh, zeros, 60.0)
        n_ok = sum(1 for _, r in reloc if r is not None)
        max_dz = max((abs(z0 - r) for z0, r in reloc if r is not None),
                     default=0.0)
        z60mp = np.array([r for _, r in reloc if r is not None])
        variants["mp60"] = recon_rho(Y_GRID, z60mp, d_lat,
                                     z60mp[-1] + d_lat / 2)

        var_summary = {}
        for name, rec in variants.items():
            res = rho_an - rec
            res_shifted = res - np.mean(res)
            var_summary[name] = {
                "resid_y0.5": float(res[np.argmin(abs(Y_GRID - 0.5))]),
                "resid_y5": float(res[np.argmin(abs(Y_GRID - 5.0))]),
                "resid_y20": float(res[np.argmin(abs(Y_GRID - 20.0))]),
                "shape_rms_rel": float(np.sqrt(np.mean(res_shifted**2))
                                       / np.sqrt(np.mean(rho_an**2))),
            }
            for y, rv, rsv in zip(Y_GRID, res, res_shifted):
                rows.append(f"resid,{name},{a},{y:.6f},{rv:.10e},{rsv:.10e}")
        for y, rv in zip(Y_GRID, rho_an):
            rows.append(f"rho_analytic,-,{a},{y:.6f},{rv:.10e},")

        # T5: sinc-artifact zero contribution estimate (analytic formula)
        k = np.arange(1, 20)
        zk = 2 * np.pi * k / hh
        sinc_est = float(np.max(
            (4.0 / np.pi) * np.max(Y_GRID) * np.sum(1.0 / zk**2)))

        report[str(a)] = {
            "T1_kappa_re": kappa.real, "T1_kappa_im": kappa.imag,
            "T2_census": census, "T2_winding": wind,
            "T2_imag_pair": {"eps": eps, "W_at_ieps": eps_resid,
                             "scan_min": scan_min, "scan_at_y": scan_at},
            "T3_A_over_y_fit": A_fit, "T3_two_over_pi": 2.0 / np.pi,
            "T3_mean_shift": shift,
            "T4_variants": var_summary,
            "T4_mp60": {"n_relocated": n_ok, "max_abs_dz": max_dz},
            "T5_sinc_est_max": sinc_est,
        }
        print(f"  T3 A/y fit: A={A_fit:.4f} vs 2/pi={2/np.pi:.4f} "
              f"shift={shift:.4f}", flush=True)
        print(f"  T4 variants: {json.dumps(var_summary, indent=1)}",
              flush=True)

    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("series,variant,a,y,resid,resid_shifted\n")
        f.write("\n".join(rows) + "\n")

    meta = {
        "experiment": "E-016b",
        "task": "2.6: root-cause E-016 V3 Poisson-reconstruction failure",
        "script": "research/numerics/scripts/E-016b-v3-rootcause.py",
        "parameters": {"a_sweep": A_SWEEP, "h": H, "lambda": LAM,
                       "y_grid": [0.1, 40.0, 400], "z_max": Z_MAX},
        "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")
    print(f"wrote {CSV.name}, {JSON.name}")
    return 0


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