#!/usr/bin/env python3
"""
verify_P142.py -- Addendum 142: GS lineshape log-symmetry failure.

This verifier reproduces the composite-Simpson numerical integration in
142_Addendum_GSLogSymmetry.tex for the Gounaris-Sakurai running-width
Breit-Wigner toy model. P142 is primarily a negative result: the proposed
log-symmetry fails at O(epsilon^2), and the correct reduction must use the full
external kernel K(s).

The numerical failure mode reproduces. The only flagged issue is a small wording
problem: NWA does not make <e^u> and <e^-u> wrong by roughly +/-10% individually;
with the displayed moments the errors are about +16.8% and -3.3%.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path

sys.path.insert(0, str(Path(__file__).resolve().parent))
from verify_common import CheckResult, Verifier


class ModernVerifier(Verifier):
    """Local adapter: tolerance logic inherited byte-identical from
    verify_common.Verifier; only the output layer is modernised."""

    def __init__(self, name: str) -> None:
        self.name = name
        self.results = []
        print(name)

    def record(self, label, ok, computed="", claimed="", detail=""):
        self.results.append(CheckResult(label, ok, computed, claimed, detail))
        n = len(self.results)
        desc = f"{label} -- {detail}" if detail else label
        print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")
        if computed != "" or claimed != "":
            print(f"        computed: {computed}")
            print(f"        claimed : {claimed}")
        return ok

    def summary(self):
        passed = sum(r.ok for r in self.results)
        failed = len(self.results) - passed
        print(f"\n{'='*60}\nRESULT: {passed} PASS / {failed} FAIL")
        return 1 if failed else 0


v = ModernVerifier("P142 -- GS Log-Symmetry Failure")

ROOT = Path(__file__).resolve().parents[1]
TEX = (ROOT / "142_Addendum_GSLogSymmetry.tex").read_text()

M_RHO = 0.77526
GAMMA_RHO = 0.1491
M_PI = 0.13957
S_MIN = 4.0 * M_PI**2
S_MAX = 4.0
N = 30000
EPS = GAMMA_RHO / M_RHO
XI = S_MIN / M_RHO**2
P0 = math.sqrt(M_RHO**2 / 4.0 - M_PI**2)


def gamma_s(s: float) -> float:
    p = math.sqrt(max(s / 4.0 - M_PI**2, 0.0))
    return GAMMA_RHO * (s / M_RHO**2) ** 1.5 * (p / P0) ** 3


def b_gs(s: float) -> float:
    width = gamma_s(s)
    return M_RHO**2 * width / ((M_RHO**2 - s) ** 2 + M_RHO**2 * width**2)


def simpson(func) -> float:
    n = N if N % 2 == 0 else N + 1
    h = (S_MAX - S_MIN) / n
    total = func(S_MIN) + func(S_MAX)
    for i in range(1, n):
        total += (4.0 if i % 2 else 2.0) * func(S_MIN + i * h)
    return total * h / 3.0


def expect(func) -> float:
    return simpson(lambda s: func(s) * b_gs(s) / s) / NORM


def u_of_s(s: float) -> float:
    return math.log(s / M_RHO**2)


NORM = simpson(lambda s: b_gs(s) / s)
EU = expect(lambda s: s / M_RHO**2)
EMU = expect(lambda s: M_RHO**2 / s)
U1 = expect(u_of_s)
U2 = expect(lambda s: u_of_s(s) ** 2)
U3 = expect(lambda s: u_of_s(s) ** 3)
DIFF = EU - EMU
GM = math.sqrt(EU * EMU)
VAR_U = U2 - U1**2
GM_CUM = 1.0 + 0.5 * VAR_U
TAYLOR_DIFF = 2.0 * U1 + U3 / 3.0


v.check("epsilon", EPS, 0.1923, rel=2e-4)
v.check("xi", XI, 0.1296, rel=4e-4)
v.check("u threshold", math.log(XI), -2.043, rel=2e-5)
f1 = 3.0 * (2.0 - XI) / (2.0 * (1.0 - XI))
v.check("f1 running-width coefficient", f1, 3.223, rel=2e-4)
v.check(
    "Gamma ratio at u=0.3",
    gamma_s(M_RHO**2 * math.exp(0.3)) / GAMMA_RHO,
    2.603,
    rel=5e-4,
)
v.check("e^u expectation", EU, 1.167581, rel=1e-6)
v.check("e^-u expectation", EMU, 0.966847, rel=1e-6)
v.check("expectation ratio", EU / EMU, 1.2076, rel=2e-5)
v.check("u mean", U1, 0.08514, rel=1e-5)
v.check("u2 mean", U2, 0.12432, rel=2e-5)
v.check("u3 mean", U3, 0.08450, rel=6e-5)
v.check("u mean in epsilon^2", U1 / EPS**2, 2.30, rel=9e-4)
v.check("u2 mean in epsilon^2", U2 / EPS**2, 3.36, rel=4e-4)
v.check("e^u minus e^-u", DIFF, 0.20073, rel=2e-5)
v.check("difference in epsilon^2", DIFF / EPS**2, 5.43, rel=7e-4)
v.check("Taylor reconstruction", TAYLOR_DIFF, 0.19845, rel=2e-5)
v.check("Taylor reconstruction percent gap", abs(100.0 * (TAYLOR_DIFF - DIFF) / DIFF), 1.14, rel=3e-3)
v.check("geometric mean direct", GM, 1.0625, rel=2e-5)
v.check("geometric mean cumulant approximation", GM_CUM, 1.0585, rel=4e-5)
v.check("geometric mean deviation percent", 100.0 * (GM - 1.0), 6.25, rel=5e-4)
v.check("cumulant approximation relative error percent", abs(100.0 * (GM_CUM - GM) / GM), 0.4, rel=8e-2)

v.record(
    "log-symmetry proof correctly fails",
    "The proof fails" in TEX and "Log-symmetry fails" in TEX,
    computed=f"<e^u>-<e^-u>={DIFF:.6f} = {DIFF / EPS**2:.3f} eps^2",
    claimed="symmetry fails at O(epsilon^2)",
)
v.record(
    "full kernel K(s) remains open",
    "full kernel $K(s)$" in TEX and "Open Items" in TEX,
    computed="pure dispersive log-symmetry is not sufficient",
    claimed="correct reduction must include K(s)",
)
v.record(
    "NWA individual errors are approximately +/-10%",
    False,
    computed=f"NWA=1 gives errors {100*(EU-1):+.2f}% for <e^u> and {100*(EMU-1):+.2f}% for <e^-u>",
    claimed="wrong by +/-10% individually",
    detail="Expected wording fail: the average scale is around 10%, but the individual errors are asymmetric.",
)

sys.exit(v.summary())