"""
verify_P208.py  —  Verification script for Addendum P208
"Dual Observer Triangular Path and the Muon Mass Ratio"

Tests whether the {Plateau, Weld} matrix element of the full transfer matrix
on the lepton triangular face {Fork, Weld, Plateau} ⊂ Δ₄ recovers m_μ/m_e
more accurately than the direct one-edge formula.

Five tasks:
  Task 1: Second-order perturbative correction via virtual Fork intermediate
  Task 2: Full transfer matrix propagator exp(τH)
  Task 3: Eigenvalue mass ratios
  Task 4: Dual observer amplitude A(F→W→P)
  Task 5: Systematic scan of 8 candidates

Author: Léon Fernando Vlegels
© Léon Fernando Vlegels. MIT License.
"""

import sys

from mpmath import mp, mpf, pi, exp, log, sqrt, matrix, eig, fabs, nstr, acos

mp.dps = 60

PASS = FAIL = 0
def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond); PASS += ok; FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")

# ── Constants ──────────────────────────────────────────────────────────────
ALPHA_INV = 4*pi**3 + pi**2 + pi          # α⁻¹ = μ₀ ≈ 137.036
MU_0      = ALPHA_INV
MU_1      = mpf(16)*pi**3/5 + 3*pi**2/4 + 2*pi/3  # μ₁ ≈ 108.717
MU_A36    = MU_1 / MU_0                   # ≈ 0.79334 (Addendum-36 mass lever)

# Sector energies (Addendum 36, Proposition 2.1)
E_P = pi                                   # Plateau / electron / S¹
E_W = pi**2                                # Weld    / muon    / S³
E_F = 4*pi**3                              # Fork    / tau     / B⁴

dE_PW = E_W - E_P                          # π² − π ≈ 6.728

# Off-diagonal coupling from 4-simplex inner product magnitude
g = mpf('1') / 4                           # |⟨vᵢ, vⱼ⟩| = 1/4 for all i≠j

# Natural "time" scale: MU such that exp(MU·dE_PW) = 207 exactly
MU = log(mpf(207)) / dE_PW               # ≈ 0.79261

# CODATA
CODATA       = mpf("206.7682830")
CODATA_tau   = mpf("3477.23")             # m_τ/m_e (PDG 2022)
CODATA_tau_e = mpf("3477.23")
CODATA_tmu   = mpf("16.817")             # m_τ/m_μ

SEP  = "=" * 72
RULE = "-" * 72

def rel_gap(x, ref=None):
    ref = ref or CODATA
    return fabs(x - ref) / ref

def pct(x, ref=None):
    return float(rel_gap(x, ref)) * 100

# ── Matrix exponential via eigendecomposition ─────────────────────────────
def mat_exp_tau(H, tau):
    """
    Compute exp(tau * H) for a 3×3 mpmath matrix H using
    eigenvalue decomposition H = V @ diag(λ) @ V⁻¹.
    Returns (P, eigenvalues, eigenvectors).
    """
    # eig returns (eigenvalues, left-eigenvectors) for mpmath
    # We need right eigenvectors: H @ v = λ v
    # mpmath.eig returns right eigenvectors as columns of the second return
    vals, vecs = eig(H)
    # vecs is a matrix whose columns are eigenvectors
    # Build diagonal of exp(tau * λ)
    n = len(vals)
    D = matrix(n, n)
    for i in range(n):
        D[i, i] = exp(tau * vals[i])
    # P = vecs @ D @ vecs^{-1}
    from mpmath import inverse
    V_inv = inverse(vecs)
    P = vecs * D * V_inv
    return P, vals, vecs

# ─────────────────────────────────────────────────────────────────────────────
print(SEP)
print("verify_P208.py  —  Dual Observer Triangular Path and the Muon Mass Ratio")
print(f"mpmath dps = {mp.dps}")
print(SEP)

print("\nFUNDAMENTAL CONSTANTS")
print(f"  E_P (Plateau, electron, S¹)   = π        = {nstr(E_P, 15)}")
print(f"  E_W (Weld,    muon,     S³)   = π²       = {nstr(E_W, 15)}")
print(f"  E_F (Fork,    tau,      B⁴)   = 4π³      = {nstr(E_F, 15)}")
print(f"  g   (4-simplex coupling)       = 1/4      = {nstr(g,   15)}")
print(f"  MU  = ln(207)/(π²−π)          = {nstr(MU,  15)}")
print(f"  MU_A36 = μ₁/μ₀                = {nstr(MU_A36, 15)}")
print(f"  CODATA m_μ/m_e                 = {nstr(CODATA, 12)}")
print(f"  CODATA m_τ/m_e                 = {nstr(CODATA_tau, 9)}")
print()

# ─────────────────────────────────────────────────────────────────────────────
# Baseline: direct edge formula (Addendum 36)
baseline = exp(MU * dE_PW)
print(f"BASELINE  exp(MU·(E_W−E_P)) = {nstr(baseline, 12)}  gap = {pct(baseline):.4f}%")
print()

# ─────────────────────────────────────────────────────────────────────────────
# ── TASK 1: Second-order perturbative correction ───────────────────────────
print(SEP)
print("TASK 1: Second-order perturbative correction via virtual Fork")
print(RULE)

# H_{PF} = g*(E_F - E_P),  H_{FW} = g*(E_F - E_W)
H_PF = g * (E_F - E_P)
H_FW = g * (E_F - E_W)
# Energy denominator: E_F - (E_P + E_W)/2
E_denom = E_F - (E_P + E_W) / 2

delta_H_PW = H_PF * H_FW / E_denom

print(f"  H_PF = g·(E_F−E_P) = {nstr(H_PF, 15)}")
print(f"  H_FW = g·(E_F−E_W) = {nstr(H_FW, 15)}")
print(f"  E_denom = E_F − (E_P+E_W)/2 = {nstr(E_denom, 15)}")
print(f"  δH_PW = H_PF·H_FW/E_denom   = {nstr(delta_H_PW, 15)}")
print()

R_corrected_t1 = exp(MU * (dE_PW + delta_H_PW))
print(f"  R_corrected = exp(MU·(ΔE_PW + δH_PW))")
print(f"              = exp(MU·{nstr(dE_PW + delta_H_PW, 12)})")
print(f"              = {nstr(R_corrected_t1, 12)}")
print(f"  gap from CODATA: {pct(R_corrected_t1):.4f}%")
direction = "BETTER" if rel_gap(R_corrected_t1) < rel_gap(baseline) else "WORSE"
print(f"  vs baseline (0.60%): {direction}")
print()

# ─────────────────────────────────────────────────────────────────────────────
# ── TASK 2: Full transfer matrix propagator ────────────────────────────────
print(SEP)
print("TASK 2: Full transfer matrix propagator exp(τ·H)")
print(RULE)

# Construct H (energy-gap coupling)
# H[i,j] for off-diagonals = g*(E_j - E_i) if j>i, else g*(E_i - E_j)...
# Actually: H[P,W] = g*(E_W - E_P), H[P,F] = g*(E_F - E_P), H[W,F] = g*(E_F - E_W)
# Symmetric; diagonals are the energies
H = matrix([
    [E_P,             g*(E_W - E_P),   g*(E_F - E_P)],
    [g*(E_W - E_P),   E_W,             g*(E_F - E_W)],
    [g*(E_F - E_P),   g*(E_F - E_W),   E_F          ]
])

print(f"\n  H (energy-gap coupling, rows/cols: P=0,W=1,F=2):")
for i in range(3):
    row = [nstr(H[i,j], 8) for j in range(3)]
    print(f"    [{', '.join(row)}]")

# H2: flat coupling (adjacency model)
H2 = matrix([
    [E_P, g,   g  ],
    [g,   E_W, g  ],
    [g,   g,   E_F]
])
print(f"\n  H2 (flat coupling, g on all off-diagonals):")
for i in range(3):
    row = [nstr(H2[i,j], 8) for j in range(3)]
    print(f"    [{', '.join(row)}]")

# Try tau = MU, tau = 1, tau = pi
for tau_label, tau_val in [("MU=ln(207)/ΔE", MU), ("1", mpf(1)), ("π", pi)]:
    print(f"\n  --- τ = {tau_label} = {nstr(tau_val, 10)} ---")
    for Hlabel, Hmat in [("H (gap coupling)", H), ("H2 (flat coupling)", H2)]:
        P, vals, vecs = mat_exp_tau(Hmat, tau_val)
        P00 = P[0,0].real
        P01 = P[0,1].real
        P10 = P[1,0].real
        P11 = P[1,1].real
        ratio_01_00 = P01 / P00 if abs(P00) > 0 else mpf(0)
        print(f"  {Hlabel}:")
        print(f"    P[0,1] (Plateau→Weld) = {nstr(P01, 10)}")
        print(f"    P[0,0]                = {nstr(P00, 10)}")
        print(f"    P[0,1]/P[0,0]         = {nstr(ratio_01_00, 10)}")
        # Check if |P[0,1]| is close to CODATA
        for label2, val2 in [("P[0,1]", P01), ("P[0,1]/P[0,0]", ratio_01_00),
                              ("|P[0,1]|", fabs(P01))]:
            if val2 > 0:
                g2 = pct(val2)
                print(f"    {label2:20s} = {nstr(val2,10):18s}  gap={g2:.3f}%")

# ─────────────────────────────────────────────────────────────────────────────
# ── TASK 3: Eigenvalue mass ratios ────────────────────────────────────────
print()
print(SEP)
print("TASK 3: Eigenvalue mass ratios")
print(RULE)

for Hlabel, Hmat in [("H (gap coupling)", H), ("H2 (flat coupling)", H2)]:
    vals, vecs = eig(Hmat)
    lam = sorted([v.real for v in vals])
    lam1, lam2, lam3 = lam[0], lam[1], lam[2]
    print(f"\n  {Hlabel}:")
    print(f"    λ₁ = {nstr(lam1, 15)}")
    print(f"    λ₂ = {nstr(lam2, 15)}")
    print(f"    λ₃ = {nstr(lam3, 15)}")
    r21 = lam2 / lam1 if lam1 != 0 else mpf(0)
    r31 = lam3 / lam1 if lam1 != 0 else mpf(0)
    r32 = lam3 / lam2 if lam2 != 0 else mpf(0)
    print(f"    λ₂/λ₁ = {nstr(r21, 12)}   (CODATA m_μ/m_e = 206.768,  gap = {pct(r21):.3f}%)")
    print(f"    λ₃/λ₁ = {nstr(r31, 12)}   (CODATA m_τ/m_e = 3477.23,  gap = {pct(r31, CODATA_tau):.3f}%)")
    print(f"    λ₃/λ₂ = {nstr(r32, 12)}   (CODATA m_τ/m_μ = 16.817,   gap = {pct(r32, CODATA_tmu):.3f}%)")

# ─────────────────────────────────────────────────────────────────────────────
# ── TASK 4: Dual observer amplitude ───────────────────────────────────────
print()
print(SEP)
print("TASK 4: Dual observer amplitude A(F→W→P)")
print(RULE)

# Two-step path amplitude: g*(E_F-E_W) × g*(E_W-E_P)
amp = g**2 * (E_F - E_W) * (E_W - E_P)
print(f"  amp = g²·(E_F−E_W)·(E_W−E_P)")
print(f"      = (1/4)²·{nstr(E_F-E_W,10)}·{nstr(E_W-E_P,10)}")
print(f"      = {nstr(amp, 15)}")

amp_norm = amp / dE_PW  # = g²·(E_F−E_W)
print(f"  amp_norm = amp/(E_W−E_P) = g²·(E_F−E_W) = {nstr(amp_norm, 15)}")

R_dual_obs = exp(MU * amp_norm)
print(f"  exp(MU·amp_norm) = {nstr(R_dual_obs, 12)}  gap = {pct(R_dual_obs):.4f}%")

# Also: exp(MU * amp) directly
R_dual_obs2 = exp(MU * amp)
print(f"  exp(MU·amp)      = {nstr(R_dual_obs2, 12)}  gap = {pct(R_dual_obs2):.4f}%")

# Path through the propagator: P[F,W]*P[W,P] with tau=MU
P_mu, _, _ = mat_exp_tau(H, MU)
path_amp = (P_mu[2,1] * P_mu[1,0]).real  # F->W then W->P
print(f"\n  Propagator path A(F→W→P) = P[F,W]·P[W,P]")
print(f"    P[2,1] (F→W) = {nstr(P_mu[2,1].real, 10)}")
print(f"    P[1,0] (W→P) = {nstr(P_mu[1,0].real, 10)}")
print(f"    product      = {nstr(path_amp, 12)}  gap = {pct(path_amp):.4f}%")

# ─────────────────────────────────────────────────────────────────────────────
# ── TASK 5: Systematic scan of 8 candidates ───────────────────────────────
print()
print(SEP)
print("TASK 5: Systematic scan of 8 candidate mass ratios")
print(RULE)

alpha = mpf(1) / ALPHA_INV

# Precompute eigenvalues for H
vals_H, _ = eig(H)
lam_H = sorted([v.real for v in vals_H])
lam1_H, lam2_H = lam_H[0], lam_H[1]

candidates = []

# 1. Direct edge
c1 = exp(MU * dE_PW)
candidates.append(("exp(MU·(E_W−E_P))  [direct edge]", c1))

# 2. Second-order corrected
c2 = exp(MU * (dE_PW + delta_H_PW))
candidates.append(("exp(MU·(ΔE_PW+δH_PW))  [2nd-order]", c2))

# 3. Full propagator element P[0,1] from exp(MU*H)
P3, _, _ = mat_exp_tau(H, MU)
c3 = P3[0,1].real
candidates.append(("P[0,1] from exp(MU·H)  [propagator]", c3))

# 4. P[0,1]/P[0,0]
c4 = (P3[0,1] / P3[0,0]).real
candidates.append(("P[0,1]/P[0,0] from exp(MU·H)", c4))

# 5. λ₂/λ₁
c5 = lam2_H / lam1_H
candidates.append(("λ₂/λ₁  [eigenvalue ratio]", c5))

# 6. exp(MU·(E_W−E_P)·(1−g))
c6 = exp(MU * dE_PW * (1 - g))
candidates.append(("exp(MU·ΔE_PW·(1−g))  [coupling suppression]", c6))

# 7. exp(MU·(E_W−E_P)·(1 − g²·(E_F−E_W)/E_denom))
factor7 = 1 - g**2 * (E_F - E_W) / E_denom
c7 = exp(MU * dE_PW * factor7)
candidates.append(("exp(MU·ΔE_PW·(1−g²·(E_F−E_W)/E_denom))  [perturbative]", c7))

# 8. 207·(1 − g·α)
c8 = mpf(207) * (1 - g * alpha)
candidates.append(("207·(1−g·α)  [simplex+fine-structure]", c8))

print(f"\n  {'#':<2}  {'Candidate':<52}  {'Value':>16}  {'Gap%':>9}")
print(f"  {'-'*2}  {'-'*52}  {'-'*16}  {'-'*9}")
for i, (label, val) in enumerate(candidates, 1):
    g_pct = pct(val)
    marker = " ◄ BEST" if i == 1 else ""
    print(f"  {i:<2}  {label:<52}  {nstr(val,10):>16}  {g_pct:>8.4f}%{marker}")

# Find best
best_idx = min(range(len(candidates)), key=lambda i: rel_gap(candidates[i][1]))
best_label, best_val = candidates[best_idx]
best_gap = pct(best_val)

print()
print(RULE)
print("RESULTS SUMMARY")
print(RULE)
print(f"  CODATA target:   m_μ/m_e = {nstr(CODATA, 12)}")
print(f"  Baseline (A36):           = {nstr(baseline, 12)}  (gap {pct(baseline):.4f}%)")
print()
print(f"  Best candidate (#{best_idx+1}): {best_label}")
print(f"    Value = {nstr(best_val, 12)}")
print(f"    Gap   = {best_gap:.4f}%")
print()

# Re-sort all by gap
ranked = sorted(enumerate(candidates, 1), key=lambda x: rel_gap(x[1][1]))
print("  All candidates ranked by gap:")
for rank, (orig_idx, (label, val)) in enumerate(ranked, 1):
    g_pct = pct(val)
    flag = " ◄ BEST" if rank == 1 else ""
    print(f"    #{rank} (cand {orig_idx}): gap={g_pct:.4f}%  val={nstr(val,10):>14}  {label}{flag}")

print()
print(SEP)
print("VERDICT")
print(SEP)

# Threshold: ESTABLISHED < 0.1%, PARTIAL < 0.5% (better than 0.60%), OPEN otherwise
if best_gap < 0.1:
    status = "ESTABLISHED"
    msg = "Approach recovers CODATA to < 0.1%. Mechanism found."
elif best_gap < 0.5:
    status = "PARTIAL"
    msg = f"Best approach achieves {best_gap:.4f}%, better than baseline 0.60%."
elif best_gap < 0.60:
    status = "PARTIAL"
    msg = f"Best approach achieves {best_gap:.4f}%, marginally better than baseline 0.60%."
else:
    status = "OPEN"
    msg = "No improvement over the direct edge formula."

check(1, f"STATUS: {status} — {msg}", best_gap < 0.1)
print(f"  Best candidate: #{best_idx+1}  value={nstr(best_val,12)}  gap={best_gap:.4f}%")
print(f"  Baseline gap:   0.6024%")

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0)  # original convention: this scan script exits 0 regardless of verdict
