E-017-thinning-measurement.py

{
  "experiment": "E-017",
  "task": "2.7: direct measurement of Cor 1.6 thinning remainder R_tau = W/F - R_0 (D-026)",
  "script": "research/numerics/scripts/E-017-thinning-measurement.py",
  "parameters": {
    "a_sweep": [
      4.0,
      6.0,
      8.0,
      10.0
    ],
    "thetas": [
      "pi",
      "zero"
    ],
    "h": 0.01,
    "lambda": -1.0,
    "y_grid": [
      0.5,
      20.0,
      40
    ]
  },
  "checks": {
    "C1_W_rescaled_vs_mp": true,
    "C2_R0_skeleton_dominant": false
  },
  "key_values": {
    "xi_flint_rel": 8.821833278620677e-17,
    "C1_worst": 3.8224510130899045e-15,
    "report": {
      "4.0_pi": {
        "max_Rtau_over_R0": 1.0471206540091407,
        "arg_M_deg": 90.0,
        "C2_R0_vs_asymptotic": 0.03663281379835252,
        "ratio_profile": {
          "0.5": 1.0471206540091407,
          "4.5": 1.0089222529454642,
          "8.5": 1.004385471663986,
          "12.5": 1.002691172170338,
          "16.5": 1.0018343882474554
        }
      },
      "4.0_zero": {
        "max_Rtau_over_R0": 1.0679333353283942,
        "arg_M_deg": -90.0,
        "C2_R0_vs_asymptotic": 0.07326332357745524,
        "ratio_profile": {
          "0.5": 1.0679333353283942,
          "4.5": 1.007167312755912,
          "8.5": 1.0034687273606262,
          "12.5": 1.0021290544442027,
          "16.5": 1.001455414041486
        }
      },
      "6.0_pi": {
        "max_Rtau_over_R0": 1.0610941706671497,
        "arg_M_deg": 90.0,
        "C2_R0_vs_asymptotic": 0.004957508160826192,
        "ratio_profile": {
          "0.5": 1.0610941706671497,
          "4.5": 1.0092001036041862,
          "8.5": 1.0045129740120455,
          "12.5": 1.0027792548957009,
          "16.5": 1.0019027560815716
        }
      },
      "6.0_zero": {
        "max_Rtau_over_R0": 1.076103329487372,
        "arg_M_deg": -90.0,
        "C2_R0_vs_asymptotic": 0.009915010610412719,
        "ratio_profile": {
          "0.5": 1.076103329487372,
          "4.5": 1.0080288607530927,
          "8.5": 1.0038984561227535,
          "12.5": 1.0023976261332865,
          "16.5": 1.0016410970456837
        }
      },
      "8.0_pi": {
        "max_Rtau_over_R0": 1.068688355229517,
        "arg_M_deg": 90.0,
        "C2_R0_vs_asymptotic": 0.0006709252652428602,
        "ratio_profile": {
          "0.5": 1.068688355229517,
          "4.5": 1.0093700235310568,
          "8.5": 1.0046033456882202,
          "12.5": 1.0028429914959156,
          "16.5": 1.0019507330410926
        }
      },
      "8.0_zero": {
        "max_Rtau_over_R0": 1.0801095258799194,
        "arg_M_deg": -90.0,
        "C2_R0_vs_asymptotic": 0.0013418505163289656,
        "ratio_profile": {
          "0.5": 1.0801095258799194,
          "4.5": 1.0084481599920176,
          "8.5": 1.0041028348136058,
          "12.5": 1.0025253381993795,
          "16.5": 1.0017307942059315
        }
      },
      "10.0_pi": {
        "max_Rtau_over_R0": 1.0733197718087135,
        "arg_M_deg": 90.0,
        "C2_R0_vs_asymptotic": 9.079985954836377e-05,
        "ratio_profile": {
          "0.5": 1.0733197718087135,
          "4.5": 1.0094431817705503,
          "8.5": 1.0046315068937992,
          "12.5": 1.0028604241200307,
          "16.5": 1.0019637968007054
        }
      },
      "10.0_zero": {
        "max_Rtau_over_R0": 1.0824561061830875,
        "arg_M_deg": -90.0,
        "C2_R0_vs_asymptotic": 0.00018159971906163639,
        "ratio_profile": {
          "0.5": 1.0824561061830875,
          "4.5": 1.008718091685042,
          "8.5": 1.0042436582330878,
          "12.5": 1.0026166536755394,
          "16.5": 1.001795570011337
        }
      }
    }
  },
  "elapsed_seconds": 9.3,
  "timestamp_utc": "2026-08-18T18:05:37.058871+00:00"
}

"""E-017 — direct measurement of the Cor 1.6 thinning remainder R_tau
(D-026): M(a,theta;iy) = W_gal(iy)/F(iy), F = z^2 xi/xi', and
R_tau = M - R_0 with the closed-form skeleton
  R_0(a,theta;z) = -2i[sinh(a(1+iz)) + e^{i th} sinh(a(iz-1))]
                   + (2i sin(az)/z)[(z-i) - e^{i th}(z+i)]
             (~= 2 e^{a+i th/2} sin(az - th/2) on |Im z| < 1).

The (star) criterion of D-026 sec.4: Cor 1.6 <=> R_tau cancels the sine
skeleton up to a zero-free remainder. Progress metric: max_y |R_tau/R_0|
vs a (1 = full cancellation; << 1 = pre-transition, expected at a <= 10
per D-023).

W_gal via the E-016c stable rescaled evaluator (validated 1.75e-16 vs
mpmath dps 50). F(iy) = -y^2 xi(1/2+y)/xi'(1/2+y) via mpmath dps 40
(E-016 conventions; flint cross-check at endpoints).

Checks:
  C1. W_gal rescaled numpy vs mpmath dps-50 reference at y in {1,10}:
      rel dev < 1e-12 (machinery regression).
  C2. exact R_0 vs asymptotic 2 e^{a+i th/2} sin(a z - th/2) on z=iy,
      y in [1,20]: rel dev < 1% (skeleton dominance).
  C3. status check (NOT a pass/fail gate): max |R_tau/R_0| recorded per
      a — the thinning-progress curve.

Outputs: research/numerics/data/E-017-thinning-measurement.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

REPO = Path(__file__).resolve().parents[3]
DATA = REPO / "research" / "numerics" / "data"
CSV = DATA / "E-017-thinning-measurement.csv"
JSON = DATA / "E-017-thinning-measurement.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, 6.0, 8.0, 10.0]
THETAS = [(np.pi, "pi"), (0.0, "zero")]
H = E16.H_MAIN
LAM = E16.LAM_SHIFT
Y_GRID = np.linspace(0.5, 20.0, 40)
mp.mp.dps = 40


def F_imag(y):
    """F(iy) = (iy)^2 xi(1/2+y)/xi'(1/2+y) (real)."""
    s = mp.mpf(0.5) + mp.mpf(y)
    return complex(-mp.mpf(y) ** 2 * E16.xi_real(s)
                   / mp.diff(E16.xi_real, s))


def R0_exact(a, th, z):
    """Closed-form skeleton (D-026 sec.2), complex z."""
    z = np.asarray(z, dtype=complex)
    with np.errstate(divide="ignore", invalid="ignore"):
        t1 = -2j * (np.sinh(a * (1 + 1j * z))
                    + np.exp(1j * th) * np.sinh(a * (1j * z - 1)))
        t2 = (2j * np.sin(a * z) / z) * ((z - 1j) - np.exp(1j * th) * (z + 1j))
    t2 = np.where(z == 0, 2j * a * ((z - 1j) - np.exp(1j * th) * (z + 1j)), t2)
    return t1 + t2


def R0_asymptotic(a, th, z):
    return 2.0 * np.exp(a + 1j * th / 2) * np.sin(a * z - th / 2)


def W_stable(ziy, cp, cm, x, h, th):
    """W(a,th;z) at z = i y via the rescaled (E-016c) evaluation:
    V = e^{i xm z} Vt, W = e^{i xm z}[(z-i)Vtp + e^{i th}(z+i)Vtm]."""
    xm = x.min()
    u = x - xm
    E = np.exp(1j * np.outer(ziy, u))
    hat = E16.hat_FT(ziy, h)
    Vtp = hat * (E @ cp)
    Vtm = hat * (E @ cm)
    return np.exp(1j * xm * ziy) * ((ziy - 1j) * Vtp
                                    + np.exp(1j * th) * (ziy + 1j) * Vtm)


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

    print("F side (mpmath dps 40, flint cross-check) ...", flush=True)
    import flint
    flint.ctx.prec = 200
    Fv = np.array([F_imag(y) for y in Y_GRID])

    def xi_flint(s):
        aa = flint.arb(str(s))
        return float(0.5 * aa * (aa - 1) * (aa / 2).gamma()
                     * flint.arb.pi() ** (-aa / 2) * aa.zeta())
    xi_err = float(max(abs(E16.xi_real(s) - xi_flint(s)) / abs(xi_flint(s))
                       for s in (1.5, 10.5, 20.5)))  # s=1: arb zeta pole
    print(f"  xi vs flint max rel: {xi_err:.2e}", flush=True)

    rows = []
    report = {}
    c1_worst = 0.0
    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)

        # C1: W at y=1,10 rescaled numpy vs mpmath dps-50 reference
        for yc in (1.0, 10.0):
            w_np = W_stable(np.array([1j * yc]), cp, cm, x, hh, np.pi)[0]
            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]
            x0 = xmm[0]
            um = [v - x0 for v in xmm]
            yy = mp.mpf(yc)
            hhmp = mp.mpf(hh)
            u_hat = yy * hhmp / 2
            hat = hhmp * (mp.sinh(u_hat) / u_hat) ** 2

            def Vt(c):
                return hat * sum(ci * mp.e ** (-yy * ui)
                                 for ci, ui in zip(c, um))
            w_mp = mp.e ** (1j * x0 * 1j * yy) * (
                (1j * yy - 1j) * Vt(cpm) - (1j * yy + 1j) * Vt(cmm))
            dev = abs(complex(w_np) - complex(w_mp)) / abs(complex(w_mp))
            c1_worst = max(c1_worst, float(dev))
        mp.mp.dps = 40

        ziy = 1j * Y_GRID
        for th, tag in THETAS:
            W = W_stable(ziy, cp, cm, x, hh, th)
            Mv = W / Fv
            R0v = R0_exact(a, th, ziy)
            R0a = R0_asymptotic(a, th, ziy)
            Rtv = Mv - R0v
            ratio = np.abs(Rtv) / np.abs(R0v)
            c2 = float(np.max(np.abs(R0v - R0a) / np.abs(R0a)))
            rep = {
                "max_Rtau_over_R0": float(np.max(ratio)),
                "arg_M_deg": float(np.degrees(np.angle(np.mean(
                    Mv / np.abs(Mv))))),
                "C2_R0_vs_asymptotic": c2,
                "ratio_profile": {f"{y:.1f}": float(r)
                                  for y, r in zip(Y_GRID[::8], ratio[::8])},
            }
            report[f"{a}_{tag}"] = rep
            for k, y in enumerate(Y_GRID):
                rows.append(f"{a},{tag},{y:.3f},{W[k].real:.10e},"
                            f"{W[k].imag:.10e},{Mv[k].real:.10e},"
                            f"{Mv[k].imag:.10e},{R0v[k].real:.10e},"
                            f"{R0v[k].imag:.10e},{Rtv[k].real:.10e},"
                            f"{Rtv[k].imag:.10e},{ratio[k]:.6e}")
            print(f"  th={tag}: max|R_tau/R_0|={rep['max_Rtau_over_R0']:.4f} "
                  f"argM={rep['arg_M_deg']:.1f}deg C2={c2:.2e}", flush=True)

    checks = {
        "C1_W_rescaled_vs_mp": bool(c1_worst < 1e-12),
        "C2_R0_skeleton_dominant": bool(all(
            report[k]["C2_R0_vs_asymptotic"] < 0.01
            for k in report)),
    }
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("a,theta,y,W_re,W_im,M_re,M_im,R0_re,R0_im,Rtau_re,Rtau_im,"
                "ratio_Rtau_R0\n")
        f.write("\n".join(rows) + "\n")
    meta = {
        "experiment": "E-017",
        "task": "2.7: direct measurement of Cor 1.6 thinning remainder "
                "R_tau = W/F - R_0 (D-026)",
        "script": "research/numerics/scripts/E-017-thinning-measurement.py",
        "parameters": {"a_sweep": A_SWEEP, "thetas": ["pi", "zero"],
                       "h": H, "lambda": LAM, "y_grid": [0.5, 20.0, 40]},
        "checks": checks,
        "key_values": {"xi_flint_rel": xi_err, "C1_worst": c1_worst,
                       "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())