#!/usr/bin/env python3
"""
verify_P063.py -- Addendum 63: neutrino mass differences from J3(O).

This verifier checks the numerical claims in 63_Addendum_NeutrinoMasses.tex.
The main negative conclusion is supported: G2 phase data encode PMNS directions,
not neutrino mass eigenvalues, and a seesaw/Jordan-invariant input remains open.
The flagged issues are local formula/proof problems: the printed TBM diagonal
formula for alpha2=alpha3 is not the formula that gives the quoted 27.7 meV,
the diagonal-dominated perturbation paragraph uses a degenerate alpha2=alpha3
case as if it were nondegenerate, and the systematic numerology search is not
reproducible from the TeX alone.
"""

from __future__ import annotations

import math
import sys
from pathlib import Path

import numpy as np

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


class Verifier(_BaseVerifier):
    """Output-layer normalization only: same checks, modern [PASS]/[FAIL] format."""

    def __init__(self, name: str):
        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)
        note, info = "", detail
        if not ok and "Expected" in detail:
            i = detail.find("Expected")
            note = " -- " + detail[i:]
            info = detail[:i].rstrip().rstrip(";")
        print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {label}{note}")
        if computed != "" or claimed != "":
            print(f"        computed: {computed}")
            print(f"        claimed : {claimed}")
        if info:
            print(f"        {info}")
        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 = Verifier("P063 -- Neutrino Mass-Difference Ratio Gap")

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

PI = math.pi
MU0 = 4 * PI**3 + PI**2 + PI
MU1 = 16 * PI**3 / 5 + 3 * PI**2 / 4 + 2 * PI / 3
MU = MU1 / MU0

DELTA21 = 7.53e-5
DELTA31 = 2.453e-3
R_PDG = DELTA21 / DELTA31
M2 = math.sqrt(DELTA21) * 1000.0  # meV
M3 = math.sqrt(DELTA31) * 1000.0  # meV


def pct(value: float, target: float) -> float:
    return 100.0 * (value - target) / abs(target)


v.check("PDG mass-squared ratio r", R_PDG, 0.0307, rel=2e-3)
v.check("m2 from Delta m21", M2, 8.68, rel=4e-4)
v.check("m3 from Delta m31", M3, 49.5, rel=6e-4)
v.check("m2/m3 equals sqrt(r)", M2 / M3, math.sqrt(R_PDG), rel=1e-12)

alpha1 = M2 / 3.0
alpha2_printed_formula = (M2**2 + M3**2 / 2.0) / (M2 + M3)
alpha2_tbm = M2 / 3.0 + M3 / 2.0
x12 = M2 / 3.0
x13 = M2 / 3.0
x23 = M3 / 2.0 - M2 / 3.0

v.check("TBM alpha1=m2/3", alpha1, 2.89, rel=1e-3)
v.check(
    "printed alpha2 formula gives quoted 27.7 meV",
    alpha2_printed_formula,
    27.7,
    rel=4e-3,
    detail="Expected fail: the displayed formula gives about 22.37 meV; the quoted value comes from m2/3 + m3/2.",
)
v.check("correct TBM alpha2=m2/3+m3/2", alpha2_tbm, 27.7, rel=2e-3)
v.check("TBM |x12|=|x13|", x12, 2.89, rel=1e-3)
v.check("TBM |x23|", x23, 21.9, rel=2e-3)

u_tbm = np.array(
    [
        [math.sqrt(2.0 / 3.0), math.sqrt(1.0 / 3.0), 0.0],
        [-math.sqrt(1.0 / 6.0), math.sqrt(1.0 / 3.0), math.sqrt(1.0 / 2.0)],
        [math.sqrt(1.0 / 6.0), -math.sqrt(1.0 / 3.0), math.sqrt(1.0 / 2.0)],
    ]
)
mass_diag = np.diag([0.0, M2, M3])
m_flavor = u_tbm @ mass_diag @ u_tbm.T
evals = np.linalg.eigvalsh(m_flavor)
v.check("direct TBM flavor matrix M11", m_flavor[0, 0], alpha1, rel=1e-12)
v.check("direct TBM flavor matrix M22", m_flavor[1, 1], alpha2_tbm, rel=1e-12)
v.check("direct TBM |M12|", abs(m_flavor[0, 1]), x12, rel=1e-12)
v.check("direct TBM |M23|", abs(m_flavor[1, 2]), x23, rel=1e-12)
v.check("direct TBM spectrum m2", evals[1], M2, rel=1e-12)
v.check("direct TBM spectrum m3", evals[2], M3, rel=1e-12)

sigma = x12**2 + x13**2 + x23**2
v.check("|x12|^2 fraction percent", 100.0 * x12**2 / sigma, 1.7, rel=7e-3)
v.check("|x13|^2 fraction percent", 100.0 * x13**2 / sigma, 1.7, rel=7e-3)
v.check("|x23|^2 fraction percent", 100.0 * x23**2 / sigma, 96.6, rel=1e-3)

sin2_theta13 = (4.0 / 9.0) * math.sin(PI / 14.0) ** 2
sin2_theta23 = 15.0 / 28.0
sin2_theta12 = math.cos(PI / 14.0) ** 2 / 3.0
v.check("corrected PMNS sin^2 theta13", sin2_theta13, 0.02201, rel=2e-4)
v.check("corrected PMNS sin^2 theta23", sin2_theta23, 0.5357, rel=3e-5)
v.check("corrected PMNS sin^2 theta12", sin2_theta12, 0.3168, rel=9e-5)

# Use the TBM sign convention from the direct flavor matrix above:
# x12 positive, x13 negative, x23 positive.
x_off = np.array([[0.0, x12, -x13], [x12, 0.0, x23], [-x13, x23, 0.0]])
off_evals = np.linalg.eigvalsh(x_off)
v.check("off-diagonal-only eigenvalue 1", off_evals[0], -22.6, rel=7e-4)
v.check("off-diagonal-only eigenvalue 2", off_evals[1], 0.7, rel=6e-2)
v.check("off-diagonal-only eigenvalue 3", off_evals[2], 21.9, rel=2e-3)
v.check("off-diagonal-only Delta21 meV^2", off_evals[1] ** 2 - off_evals[0] ** 2, -511.0, rel=6e-4)

v.record(
    "nondegenerate perturbation formula is valid for the TBM alpha2/alpha3 block",
    False,
    computed=f"alpha2-alpha3 = {alpha2_tbm - alpha2_tbm:.6g}; denominator alpha3-alpha2 vanishes",
    claimed="delta_i^(2) = -sum |x_ij|^2/(alpha_j-alpha_i) applies directly",
    detail="Expected proof-audit fail: the 2-3 block is degenerate at TBM, so ordinary nondegenerate second-order perturbation theory is not applicable.",
)
v.check(
    "largest off-diagonal perturbative size is about 1 percent",
    (x23 / alpha2_tbm) ** 2,
    0.01,
    rel=2e-1,
    detail="Expected fail: the 2-3 entry is about 79% of alpha2, so the squared ratio is about 62.5%, not 1%.",
)

e_50mev = PI + math.log((50.0e-3 * 1e-6) / 0.511) / MU
v.check("50 meV mass energy coordinate", e_50mev, -17.0, rel=2e-2)
v.record("50 meV energy lies below electron floor pi", e_50mev < PI, computed=e_50mev, claimed="< pi")

v.check("candidate 1/(4*pi^2)", 1.0 / (4.0 * PI**2), 0.02533, rel=2e-4)
v.check("candidate 1/(4*pi^2) residual percent", pct(1.0 / (4.0 * PI**2), R_PDG), -17.5, rel=2e-3)
v.check("candidate pi/mu0", PI / MU0, 0.02293, rel=3e-4)
v.check("candidate pi/mu0 residual percent", pct(PI / MU0, R_PDG), -25.3, rel=1e-3)
v.check("candidate (4/9)sin^2(pi/14)", sin2_theta13, 0.02201, rel=2e-4)
v.check("candidate sin^2(pi/14)", math.sin(PI / 14.0) ** 2, 0.04952, rel=1e-4)
v.check("candidate (pi/14)^2", (PI / 14.0) ** 2, 0.05036, rel=1e-4)
v.check("candidate 1/(2*pi^2)", 1.0 / (2.0 * PI**2), 0.05066, rel=1e-4)
v.check("candidate 1/(4*pi)", 1.0 / (4.0 * PI), 0.07958, rel=4e-5)
v.check(
    "candidate sin^2(2pi/14)-sin^2(pi/14)",
    math.sin(2.0 * PI / 14.0) ** 2 - math.sin(PI / 14.0) ** 2,
    0.13874,
    rel=4e-5,
)

hopf_m = [abs(math.sin(2.0 * i * PI / 14.0)) for i in (1, 2, 3)]
hopf_r = (hopf_m[1] ** 2 - hopf_m[0] ** 2) / (hopf_m[2] ** 2 - hopf_m[0] ** 2)
v.check("Hopf-lapse edge-spacing ratio", hopf_r, 0.55, rel=1e-2)
v.check("Hopf-lapse factor above PDG ratio", hopf_r / R_PDG, 18.0, rel=5e-3)

v.record(
    "phase-level negative conclusion is explicit",
    "not determined by the off-diagonal phase structure alone" in TEX
    and "require the Jordan invariants" in TEX
    and "Phase~5" in TEX,
    computed="phase data -> PMNS angles; Jordan invariants/seesaw remain needed",
    claimed="mass ratio not closed by G2 phases",
)
v.record(
    "systematic numerology search is reproducible from included code/candidate generator",
    False,
    computed="only the resulting candidate table is included in the TeX",
    claimed="systematic search over TOE-natural combinations",
    detail="Expected reproducibility fail: the table checks, but the candidate-generation code and full search space are not present.",
)

sys.exit(v.summary())