E-018-descent-profile.py

{
  "experiment": "E-018",
  "task": "2.8: descent profile of |M|/|R_0| on z = x + i y0 (D-027 sec.3)",
  "script": "research/numerics/scripts/E-018-descent-profile.py",
  "parameters": {
    "a_sweep": [
      4.0,
      6.0,
      8.0,
      10.0
    ],
    "thetas": [
      "pi",
      "zero"
    ],
    "h": 0.01,
    "lambda": -1.0,
    "x_grid": [
      -20.0,
      20.0,
      81
    ],
    "y0_list": [
      1.0,
      0.3,
      0.1
    ]
  },
  "checks": {
    "C1_W_offaxis_vs_mp": true,
    "C2_F_complex_vs_acb": true
  },
  "key_values": {
    "C1_worst": 6.910665131911271e-15,
    "C2_worst": 5.272337420798375e-10,
    "report": {
      "4.0_pi_y1.0": {
        "median": 0.0058782721256475925,
        "mean": 0.008291178962858501,
        "p90": 0.019481733636454057,
        "max": 0.03396554015684706
      },
      "4.0_pi_y0.3": {
        "median": 0.006564175181002058,
        "mean": 0.010435809591546516,
        "p90": 0.022551552466665042,
        "max": 0.07278660391646481
      },
      "4.0_pi_y0.1": {
        "median": 0.006707751880242824,
        "mean": 0.012108924528906726,
        "p90": 0.0228730299200152,
        "max": 0.14810301480988547
      },
      "4.0_zero_y1.0": {
        "median": 0.004555296424658213,
        "mean": 0.00664953683969302,
        "p90": 0.015072268250518295,
        "max": 0.03396554015684706
      },
      "4.0_zero_y0.3": {
        "median": 0.005025643749280622,
        "mean": 0.009101172413255655,
        "p90": 0.016734980735750398,
        "max": 0.11312570363787529
      },
      "4.0_zero_y0.1": {
        "median": 0.0050887171054875445,
        "mean": 0.012271728907487723,
        "p90": 0.022571197533653917,
        "max": 0.32130537611362625
      },
      "6.0_pi_y1.0": {
        "median": 0.005941393135128255,
        "mean": 0.008489600378158872,
        "p90": 0.019841105716330586,
        "max": 0.03790590678226328
      },
      "6.0_pi_y0.3": {
        "median": 0.006536290457633615,
        "mean": 0.010806263351104477,
        "p90": 0.022677423020264244,
        "max": 0.07232104437585961
      },
      "6.0_pi_y0.1": {
        "median": 0.006724513347452873,
        "mean": 0.011478890149315125,
        "p90": 0.023648740968485068,
        "max": 0.09142267471562603
      },
      "6.0_zero_y1.0": {
        "median": 0.005124172727938264,
        "mean": 0.00746876580462103,
        "p90": 0.017026245677604794,
        "max": 0.03790590678226328
      },
      "6.0_zero_y0.3": {
        "median": 0.005650623829821188,
        "mean": 0.01021180418503125,
        "p90": 0.018928933792391946,
        "max": 0.12699513605911442
      },
      "6.0_zero_y0.1": {
        "median": 0.005579365441718518,
        "mean": 0.014022393997392265,
        "p90": 0.02552089348096885,
        "max": 0.3808913219434005
      },
      "8.0_pi_y1.0": {
        "median": 0.005995466672509657,
        "mean": 0.008644983185117708,
        "p90": 0.020036607368895728,
        "max": 0.03989553973139515
      },
      "8.0_pi_y0.3": {
        "median": 0.00680712183885795,
        "mean": 0.011120526604333894,
        "p90": 0.022652043996305035,
        "max": 0.09017250858015384
      },
      "8.0_pi_y0.1": {
        "median": 0.006712953176200126,
        "mean": 0.011681318246964178,
        "p90": 0.02287777460583436,
        "max": 0.07957427900229018
      },
      "8.0_zero_y1.0": {
        "median": 0.005392890424355264,
        "mean": 0.007844722503562416,
        "p90": 0.017950018792760015,
        "max": 0.03989553973139515
      },
      "8.0_zero_y0.3": {
        "median": 0.005944031943257411,
        "mean": 0.010734555526343437,
        "p90": 0.019952286436256017,
        "max": 0.13369578063199172
      },
      "8.0_zero_y0.1": {
        "median": 0.005893840548583633,
        "mean": 0.014735232502103948,
        "p90": 0.02689214879887816,
        "max": 0.40147689610052184
      },
      "10.0_pi_y1.0": {
        "median": 0.006039113063356135,
        "mean": 0.008720513184988234,
        "p90": 0.02015010481452299,
        "max": 0.041077396588368544
      },
      "10.0_pi_y0.3": {
        "median": 0.006834329345532207,
        "mean": 0.011358219449539573,
        "p90": 0.02269955693499299,
        "max": 0.10221282132184713
      },
      "10.0_pi_y0.1": {
        "median": 0.006785856028914855,
        "mean": 0.012429028565700174,
        "p90": 0.026315834945999973,
        "max": 0.10947969136140419
      },
      "10.0_zero_y1.0": {
        "median": 0.005552510043266734,
        "mean": 0.008080195039525702,
        "p90": 0.018485879275465477,
        "max": 0.041077396588368544
      },
      "10.0_zero_y0.3": {
        "median": 0.006112777419801557,
        "mean": 0.011042350086136783,
        "p90": 0.02054129578864212,
        "max": 0.13759688657079275
      },
      "10.0_zero_y0.1": {
        "median": 0.006202458156384956,
        "mean": 0.015145435561662934,
        "p90": 0.0276428672378227,
        "max": 0.4131908062703263
      }
    }
  },
  "elapsed_seconds": 10.7,
  "timestamp_utc": "2026-08-18T18:27:37.463165+00:00"
}

"""E-018 — descent profile: |M|/|R_0| on lines z = x + i y0, y0 -> 0
(D-027 sec.3): the (star) cancellation is ~99.8% complete on the
imaginary axis at a<=10 but the lattice survives on the real axis. The
Cor 1.6 transition = the cancellation DESCENDING to y=0. This experiment
measures the profile: M(a,th;z) = W_gal(z)/F(z), R_0 closed form
(D-026), leftover-after-cancellation ratio |M|/|R_0| on
z = x + i y0, x in [-20,20], y0 in {1, 0.3, 0.1}, a = 4..10,
theta = pi (theta = 0 control).

Interpretation: |M|/|R_0| ~ 2e-3..8e-2 on the imaginary axis (E-017);
skeleton dominance near the real axis (pre-transition) shows as
O(1)-or-larger ratios and spiking at the sine skeleton zeros
x = th/(2a) + pi k/a. A descending cancellation (transition onset) would
show the small-y0 ratios falling toward the imaginary-axis values.

Checks:
  C1. W_gal rescaled vs mpmath dps-50 at z = 5 + 0.3i, 12 + i (a=4,10):
      rel dev < 1e-10 (machinery regression off-axis).
  C2. F(z) at complex z vs flint acb xi/xi' (10 sample points):
      rel dev < 1e-10.

Outputs: research/numerics/data/E-018-descent-profile.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-018-descent-profile.csv"
JSON = DATA / "E-018-descent-profile.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
X_GRID = np.linspace(-20.0, 20.0, 81)
Y0_LIST = [1.0, 0.3, 0.1]
mp.mp.dps = 40


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 xi_complex(s):
    return 0.5 * s * (s - 1) * mp.pi ** (-s / 2) * mp.gamma(s / 2) * mp.zeta(s)


def F_complex(z):
    """F(z) = z^2 xi(1/2 - i z)/xi'(1/2 - i z), mpmath dps 40."""
    zc = complex(z)
    s = mp.mpc(0.5 + zc.imag, -zc.real)          # s = 1/2 - i z
    val = xi_complex(s)
    dval = mp.diff(xi_complex, s)
    return complex(zc**2 * val / dval)


def W_stable_gen(z, cp, cm, x, h, th):
    """W(a,th;z) for Im z > 0 via the rescaled evaluation (E-016c)."""
    xm = x.min()
    u = x - xm
    E = np.exp(1j * np.outer(z, u))
    hat = E16.hat_FT(z, h)
    Vtp = hat * (E @ cp)
    Vtm = hat * (E @ cm)
    return np.exp(1j * xm * z) * ((z - 1j) * Vtp
                                  + np.exp(1j * th) * (z + 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()

    # ---- F on the grid (z only, a-independent)
    print("F on grid (mpmath dps 40) + flint acb cross-check ...",
          flush=True)
    zpts = np.array([x + 1j * y0 for y0 in Y0_LIST for x in X_GRID])
    Fv = np.array([F_complex(z) for z in zpts])

    import flint
    flint.ctx.prec = 200

    def xi_acb(s):
        aa = flint.acb(s.real, s.imag)
        return complex(0.5 * aa * (aa - 1) * (aa / 2).gamma()
                       * flint.arb.pi() ** (-aa / 2) * aa.zeta())
    c2 = 0.0
    eps_fd = 1e-6
    for z in zpts[::25]:
        zc = complex(z)
        s = 0.5 - 1j * zc
        xi0 = xi_acb(s)
        xid = complex(mp.diff(xi_complex, mp.mpc(s.real, s.imag)))
        xid_acb = (xi_acb(s + eps_fd) - xi_acb(s - eps_fd)) / (2 * eps_fd)
        F_ref = zc**2 * xi0 / xid_acb
        F_mp = zc**2 * xi0 / xid          # acb xi, mp xi' — isolate xi'
        c2 = max(c2, abs(F_mp - F_ref) / abs(F_ref))
    print(f"  C2 xi' mp-vs-acb (fd) max rel: {c2:.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: off-axis machinery regression vs mpmath dps 50
        for zc in (5 + 0.3j, 12 + 1.0j):
            w_np = W_stable_gen(np.array([zc]), 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]
            zz = mp.mpc(zc.real, zc.imag)
            hhmp = mp.mpf(hh)
            u_hat = zz * hhmp / 2
            hat = hhmp * (mp.sin(u_hat) / u_hat) ** 2

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

        for th, tag in THETAS:
            for iy, y0 in enumerate(Y0_LIST):
                zg = X_GRID + 1j * y0
                W = W_stable_gen(zg, cp, cm, x, hh, th)
                Fg = Fv[iy * len(X_GRID):(iy + 1) * len(X_GRID)]
                Mv = W / Fg
                R0v = R0_exact(a, th, zg)
                ratio = np.abs(Mv) / np.abs(R0v)
                rep = {
                    "median": float(np.median(ratio)),
                    "mean": float(np.mean(ratio)),
                    "p90": float(np.quantile(ratio, 0.9)),
                    "max": float(np.max(ratio)),
                }
                report[f"{a}_{tag}_y{y0}"] = rep
                for k, xx in enumerate(X_GRID):
                    rows.append(f"{a},{tag},{y0},{xx:.3f},{Mv[k].real:.10e},"
                                f"{Mv[k].imag:.10e},{R0v[k].real:.10e},"
                                f"{R0v[k].imag:.10e},{ratio[k]:.6e}")
            r1 = report[f"{a}_{tag}_y1.0"]["median"]
            r3 = report[f"{a}_{tag}_y0.3"]["median"]
            r01 = report[f"{a}_{tag}_y0.1"]["median"]
            print(f"  th={tag}: median|M/R0| y0=1: {r1:.3e}  "
                  f"0.3: {r3:.3e}  0.1: {r01:.3e}", flush=True)

    checks = {
        "C1_W_offaxis_vs_mp": bool(c1_worst < 1e-10),
        "C2_F_complex_vs_acb": bool(c2 < 1e-8),
    }
    with CSV.open("w", encoding="utf-8", newline="") as f:
        f.write("a,theta,y0,x,M_re,M_im,R0_re,R0_im,ratio_M_R0\n")
        f.write("\n".join(rows) + "\n")
    meta = {
        "experiment": "E-018",
        "task": "2.8: descent profile of |M|/|R_0| on z = x + i y0 "
                "(D-027 sec.3)",
        "script": "research/numerics/scripts/E-018-descent-profile.py",
        "parameters": {"a_sweep": A_SWEEP, "thetas": ["pi", "zero"],
                       "h": H, "lambda": LAM,
                       "x_grid": [-20.0, 20.0, 81], "y0_list": Y0_LIST},
        "checks": checks,
        "key_values": {"C1_worst": c1_worst, "C2_worst": c2,
                       "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())