E-029-convexity.py

{
  "experiment": "E-029",
  "task": "4.5b support: D-057 convexity proof verification",
  "script": "research/numerics/scripts/E-029-convexity.py",
  "checks": {
    "V1_sympy_matches_hand": false,
    "V2_psi_negative": true,
    "V3_margin_tilde_convex_positive": true,
    "V4_I1_second_derivative": true,
    "V5_majorant": true,
    "table": [
      {
        "block": "T log T",
        "max_rel_dev": 3.86713266584926e-16,
        "max_abs_dev": 1.1102230246251565e-16,
        "psi_max_on_grid": -9.278201041106209e-12,
        "negative_everywhere": true
      },
      {
        "block": "T",
        "max_rel_dev": 0.0,
        "max_abs_dev": 0.0,
        "psi_max_on_grid": -3.225806451613944e-13,
        "negative_everywhere": true
      },
      {
        "block": "log T",
        "max_rel_dev": 2.884096233990709e-16,
        "max_abs_dev": 1.3877787807814457e-17,
        "psi_max_on_grid": -2.0811654526538217e-25,
        "negative_everywhere": true
      },
      {
        "block": "log log T",
        "max_rel_dev": Infinity,
        "max_abs_dev": Infinity,
        "psi_max_on_grid": -7.361493795231685e-27,
        "negative_everywhere": true
      },
      {
        "block": "T L1",
        "max_rel_dev": 4.18537251328277e-16,
        "max_abs_dev": 1.1102230246251565e-16,
        "psi_max_on_grid": -9.278201041106209e-12,
        "negative_everywhere": true
      },
      {
        "block": "T L2",
        "max_rel_dev": 4.87560538773749e-16,
        "max_abs_dev": 4.163336342344337e-17,
        "psi_max_on_grid": -1.0839608294029891e-12,
        "negative_everywhere": true
      },
      {
        "block": "L1^2",
        "max_rel_dev": 3.8295728938905445e-16,
        "max_abs_dev": 1.3877787807814457e-17,
        "psi_max_on_grid": -1.1763755765826322e-23,
        "negative_everywhere": true
      },
      {
        "block": "L1 L2",
        "max_rel_dev": 4.506079009424732e-16,
        "max_abs_dev": 1.3877787807814457e-17,
        "psi_max_on_grid": -9.035766896079908e-25,
        "negative_everywhere": true
      },
      {
        "block": "L2^2",
        "max_rel_dev": Infinity,
        "max_abs_dev": Infinity,
        "psi_max_on_grid": -4.920397462312061e-26,
        "negative_everywhere": true
      },
      {
        "block": "L1",
        "max_rel_dev": 2.2639034786574087e-16,
        "max_abs_dev": 6.938893903907228e-18,
        "psi_max_on_grid": -2.081165452651808e-25,
        "negative_everywhere": true
      },
      {
        "block": "L2",
        "max_rel_dev": Infinity,
        "max_abs_dev": Infinity,
        "psi_max_on_grid": -7.361493795224352e-27,
        "negative_everywhere": true
      }
    ]
  },
  "key_values": {
    "N0_tilde": 75621,
    "margin_tilde_8e4": 16020.1971966516,
    "margin_tilde_prime_8e4": 3.6877762759686448,
    "max_2n_B4d2": 0.18932637317273204
  },
  "elapsed_seconds": 1.0
}

"""E-029 — task 4.5b support: D-057 convexity proof verification.

  V1. sympy derives psi_i(T) = phi''(T) - phi'(T)/(T-1) for the 11
      building blocks; matches D-057's hand-derived forms (rel dev).
  V2. psi_i(T) < 0 on T in [T0, 3.1e12] (log grid) — supports the
      hand sign arguments.
  V3. margin-tilde (Chat_0 closed-form majorant version) convex on
      [20, 4e5] (2nd differences > 0), positivity contiguous on
      [N0_tilde, 4e5], margin-tilde(8e4) > 0, margin-tilde'(8e4) > 0.
  V4. 2n|B4''(n)| < 1 on [20, 1e4] (I_I second-derivative bound).
  V5. B-tilde >= B pointwise on the E-028 grid (majorant sanity).

Outputs: research/numerics/data/E-029-convexity.csv (+ .json)
"""
import json
import sys
import time
from pathlib import Path

import mpmath as mp
import numpy as np
import sympy as sp

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

A_TR, B0_TR = 10.5, 48.0
M0 = 1.335180201468999
GAMMA = 0.577215664901532860606512090082402431
C1L = (GAMMA - 1 - np.log(2 * np.pi)) / 2 - GAMMA / 2
T0 = np.sqrt(20.0) + 1.0


def build_table():
    """The 11 building blocks: (name, phi(T) sympy, psi hand form)."""
    T = sp.symbols('T', positive=True)
    L1 = sp.log(T + 2)
    L2 = sp.log(sp.log(T + 2))
    phis = [
        ("T log T", T * sp.log(T),
         lambda t: 1 / t - (np.log(t) + 1) / (t - 1)),
        ("T", T,
         lambda t: -1 / (t - 1)),
        ("log T", sp.log(T),
         lambda t: -1 / t ** 2 - 1 / (t * (t - 1))),
        ("log log T", sp.log(sp.log(T)),
         lambda t: -(np.log(t) + 1) / (t * np.log(t)) ** 2
         - 1 / (t * np.log(t) * (t - 1))),
        ("T L1", T * L1,
         lambda t: (t + 4) / (t + 2) ** 2
         - np.log(t + 2) / (t - 1) - t / ((t + 2) * (t - 1))),
        ("T L2", T * L2,
         lambda t: (t * (np.log(t + 2) - 1) + 4 * np.log(t + 2))
         / ((t + 2) * np.log(t + 2)) ** 2
         - np.log(np.log(t + 2)) / (t - 1)
         - t / ((t + 2) * np.log(t + 2) * (t - 1))),
        ("L1^2", L1 ** 2,
         lambda t: 2 * (1 - np.log(t + 2)) / (t + 2) ** 2
         - 2 * np.log(t + 2) / ((t + 2) * (t - 1))),
        ("L1 L2", L1 * L2,
         lambda t: (1 / np.log(t + 2) - np.log(np.log(t + 2)) - 1)
         / (t + 2) ** 2
         - (np.log(np.log(t + 2)) + 1) / ((t + 2) * (t - 1))),
        ("L2^2", L2 ** 2,
         lambda t: 2 * (1 - np.log(np.log(t + 2))
                        * (np.log(t + 2) + 1))
         / ((t + 2) * np.log(t + 2)) ** 2
         - 2 * np.log(np.log(t + 2))
         / ((t + 2) * np.log(t + 2) * (t - 1))),
        ("L1", L1,
         lambda t: -1 / (t + 2) ** 2 - 1 / ((t + 2) * (t - 1))),
        ("L2", L2,
         lambda t: -(np.log(t + 2) + 1) / ((t + 2) * np.log(t + 2)) ** 2
         - 1 / ((t + 2) * np.log(t + 2) * (t - 1))),
    ]
    rows = []
    tpts = np.geomspace(T0, 1e6, 300)      # symbolic-match grid
    tpts_full = np.geomspace(T0, 3.1e12, 400)  # negativity grid
    for name, phi, psi_hand in phis:
        psi_sym = sp.lambdify(
            T, sp.simplify(sp.diff(phi, T, 2)
                           - sp.diff(phi, T) / (T - 1)), "numpy")
        v_sym = psi_sym(tpts)
        v_hand = psi_hand(tpts)
        diff = np.abs(v_sym - v_hand)
        rel = float(np.max(np.where(np.abs(v_hand) > 1e-15,
                                    diff / np.abs(v_hand), 0.0)))
        absd = float(np.max(diff))
        rows.append({"block": name, "max_rel_dev": rel,
                     "max_abs_dev": absd,
                     "psi_max_on_grid":
                         float(max(np.max(v_sym),
                                   np.max(psi_hand(tpts_full)))),
                     "negative_everywhere":
                         bool(np.max(v_sym) < 0
                              and np.max(psi_hand(tpts_full)) < 0)})
    return rows


def Chat(t):
    return (0.433 * np.log(t + 2) + 0.886 * np.log(np.log(t + 2)) + 2.9)


def C0_honest(t):
    def F(u):
        return (u / (2 * np.pi)) * np.log(u / (2 * np.pi)) - u / (2 * np.pi) \
            + 7 / 8

    def Sb(u):
        return 0.137 * np.log(u) + 0.443 * np.log(np.log(u)) + 1.588
    return (F(t + 1) - F(t)) + Sb(t + 1) + Sb(t)


def B_tilde(n):
    n = np.asarray(n, dtype=np.float64)
    T = np.sqrt(n) + 1.0
    C = Chat(T)
    H = 3 * C ** 2 + A_TR * np.log(T + 2) + B0_TR
    twoN = 2 * ((T / (2 * np.pi)) * np.log(T / (2 * np.pi))
                - T / (2 * np.pi) + 7 / 8
                + 0.137 * np.log(T) + 0.443 * np.log(np.log(T)) + 1.588)
    I3 = 4 * T * ((2 / 3) * C + A_TR * np.log(T + 2) + B0_TR)
    I24 = 2 * (np.e ** 4 + 1) * (7 * H + 3.3)
    r = np.sqrt(n)
    logI1 = np.log(2 * (r + 1) * 0.5700) + r / 2 - (2 * r - 2) * np.log(2)
    I1 = np.where(logI1 > -700, np.exp(np.maximum(logI1, -745)), 0.0)
    return twoN + I3 + 40.0 * M0 + I24 + I1 + 1.0


def B_honest(n):
    n = np.asarray(n, dtype=np.float64)
    T = np.sqrt(n) + 1.0
    C = C0_honest(T)
    H = 3 * C ** 2 + A_TR * np.log(T + 2) + B0_TR
    twoN = 2 * ((T / (2 * np.pi)) * np.log(T / (2 * np.pi))
                - T / (2 * np.pi) + 7 / 8
                + 0.137 * np.log(T) + 0.443 * np.log(np.log(T)) + 1.588)
    I3 = 4 * T * ((2 / 3) * C + A_TR * np.log(T + 2) + B0_TR)
    I24 = 2 * (np.e ** 4 + 1) * (7 * H + 3.3)
    r = np.sqrt(n)
    logI1 = np.log(2 * (r + 1) * 0.5700) + r / 2 - (2 * r - 2) * np.log(2)
    I1 = np.where(logI1 > -700, np.exp(np.maximum(logI1, -745)), 0.0)
    return twoN + I3 + 40.0 * M0 + I24 + I1 + 1.0


def margin_tilde(n):
    n = np.asarray(n, dtype=np.float64)
    return 0.5 * n * np.log(n) + C1L * n - 4.0 - B_tilde(n) + 1.0


def B4_d2(n):
    """Second derivative of I_I(n) = 2(sqrt(n)+1)(e^{sqrt(n)/2}+1)*0.57
    * 2^{-(2 sqrt(n)-2)} via complex-step-free central differences at
    dps 40."""
    mp.mp.dps = 40

    def f(x):
        r = mp.sqrt(x)
        return 2 * (r + 1) * (mp.e ** (r / 2) + 1) * mp.mpf("0.57") \
            * 2 ** (-(2 * r - 2))

    h = mp.mpf("0.01")
    return float((f(n + h) - 2 * f(n) + f(n - h)) / h ** 2)


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

    print("V1/V2: building-block table ...", flush=True)
    rows = build_table()
    v1 = all(r["max_rel_dev"] < 1e-10 and r["max_abs_dev"] < 1e-12
             for r in rows)
    v2 = all(r["negative_everywhere"] for r in rows)
    for r in rows:
        print(f"  {r['block']:>10}: rel_dev {r['max_rel_dev']:.2e}, "
              f"psi_max {r['psi_max_on_grid']:.3e}", flush=True)

    print("V3: margin-tilde convexity/positivity ...", flush=True)
    ngrid = np.arange(20, 400001)
    mt = margin_tilde(ngrid)
    conv = bool(np.all(np.diff(mt, 2) > 0))
    pos = ngrid[(mt > 0) & (np.gradient(mt, ngrid) > 0)]
    N0t = int(pos[0]) if len(pos) else -1
    cont = bool(np.all(mt[ngrid >= N0t] > 0)) if N0t > 0 else False
    m8 = float(margin_tilde(80000.0))
    d8 = float((margin_tilde(80000.5) - margin_tilde(79999.5)))
    print(f"  convex on [20,4e5]: {conv}; N0_tilde = {N0t} "
          f"(contiguous {cont}); margin~(8e4) = {m8:.1f}, "
          f"margin~'(8e4) = {d8:.1f}", flush=True)
    v3 = bool(conv and cont and m8 > 0 and d8 > 0)

    print("V4: I_I second-derivative bound ...", flush=True)
    n4 = np.arange(20, 10001, 5)
    ratios = [2 * n * abs(B4_d2(float(n))) for n in n4]
    v4 = bool(max(ratios) < 1)
    print(f"  max 2n|B4''| = {max(ratios):.4f} at n = "
          f"{n4[int(np.argmax(ratios))]} (< 1: {v4})", flush=True)

    print("V5: majorant sanity ...", flush=True)
    ng = np.arange(20, 400001, 20)
    v5 = bool(np.all(B_tilde(ng) >= B_honest(ng)))
    print(f"  B_tilde >= B_honest on grid: {v5}", flush=True)

    with open(CSV, "w", newline="") as f:
        import csv as csvmod
        w = csvmod.writer(f)
        w.writerow(["n", "margin_tilde"])
        for n in range(20, 400001, 20):
            w.writerow([n, f"{float(margin_tilde(n)):.6f}"])

    meta = {
        "experiment": "E-029",
        "task": "4.5b support: D-057 convexity proof verification",
        "script": "research/numerics/scripts/E-029-convexity.py",
        "checks": {
            "V1_sympy_matches_hand": v1,
            "V2_psi_negative": v2,
            "V3_margin_tilde_convex_positive": v3,
            "V4_I1_second_derivative": v4,
            "V5_majorant": v5,
            "table": rows,
        },
        "key_values": {"N0_tilde": N0t, "margin_tilde_8e4": m8,
                       "margin_tilde_prime_8e4": d8,
                       "max_2n_B4d2": float(max(ratios))},
        "elapsed_seconds": round(time.time() - t0, 1),
    }
    JSON.write_text(json.dumps(meta, indent=2))
    print(f"  wrote {CSV.name}, {JSON.name}; elapsed "
          f"{meta['elapsed_seconds']}s")
    return 0 if all([v1, v2, v3, v4, v5]) else 1


if __name__ == "__main__":
    sys.exit(main())