#!/usr/bin/env python3
# verify_P321.py
# Copyright Léon Fernando Vlegels -- CC BY 4.0
#
# Verifier for A321: the selection-rule framing audit (OI-287-1).
# Recomputes from scratch (numpy only) and asserts the AUDITED FACTS.
# Verdict: FRAMING-DRIFT.
#
# Corpus citations encoded as comments at each check:
#   - (2,3,16) = Taylor coefficients of rho:
#       287_Addendum_LeptonCorrectionLaw.tex lines 40-54, 105-119
#       296_Addendum_IndexMatching.tex lines 82-95 (c_k = rho^(k)(0)/(pi^k k!))
#   - the rule is a lepton MASS correction law, not an overlap-diagonal:
#       287 lines 49-58 (beta_n = 6 mu1/mu0 - c_n alpha^(1/4)/lambda_1)
#   - the overlap-diagonal is the EXCLUDED bare Taylor probe:
#       296 lines 97-107 (Third exclusion); 297 lines 61-77 (position moments)
#       revived as a "lead": 299 lines 46-55
#   - the moments mu_n = int x^n rho dx are the corpus's hierarchy object,
#     monotone-decreasing, NOT equal to (2,3,16):
#       20_Paper_StandardPhysicsEmbedding.tex lines 120,133-142,451,459-461

import sys
from math import factorial
import numpy as np

PI = np.pi
PASS = 0
FAIL = 0


def check(n, desc, cond):
    global PASS, FAIL
    if cond:
        PASS += 1
        print(f"  [PASS] ({n}) {desc}")
    else:
        FAIL += 1
        print(f"  [FAIL] ({n}) {desc}")


# ---- recompute from scratch ----
def taylor_coeff(k):
    derivs = {1: 2 * PI, 2: 6 * PI**2, 3: 96 * PI**3}
    return derivs[k] / (PI**k * factorial(k))


def rho_prime(x):
    return 48 * PI**3 * x**2 + 6 * PI**2 * x + 2 * PI


def moment(n):
    return 2 * PI / (n + 2) + 3 * PI**2 / (n + 3) + 16 * PI**3 / (n + 4)


def build_eigs(N=1500):
    Omega = 4 * PI**3 + PI**2 + PI
    E_self = 13.177
    x = np.linspace(0.0, 1.0, N + 2)[1:-1]
    h = 1.0 / (N + 1)
    V = (rho_prime(x))**2 / (2 * Omega**2) + E_self * x**2 * (1 - x)**2
    main = 2.0 / h**2 + V
    off = -1.0 / h**2 * np.ones(N - 1)
    H = np.diag(main) + np.diag(off, 1) + np.diag(off, -1)
    w, v = np.linalg.eigh(H)
    for j in range(v.shape[1]):
        v[:, j] /= np.sqrt(np.sum(v[:, j]**2) * h)
    return x, h, w[:3], v[:, :3]


print("=" * 60)
print("verify_P321 :: selection-rule framing audit (OI-287-1)")
print("=" * 60)

print("\n[ORIGIN] (2,3,16) are rho's Taylor coefficients")
c = np.array([taylor_coeff(1), taylor_coeff(2), taylor_coeff(3)])
# 287 lines 40-54, 105-119; 296 lines 82-95
check(1, "Taylor c_k = rho^(k)(0)/(pi^k k!) gives (2,3,16)",
      np.allclose(c, [2.0, 3.0, 16.0], atol=1e-12))
# rho' = 48pi^3 x^2 + 6pi^2 x + 2pi integrates to rho = 16pi^3 x^3 + 3pi^2 x^2 + 2pi x
# coefficients in pi^k units at degree 1,2,3 are exactly (2,3,16) -- 287 line 42
check(2, "antiderivative arithmetic: rho' -> rho coeffs (2,3,16) in pi^k units",
      abs(48.0 / 3 - 16.0) < 1e-12 and abs(6.0 / 2 - 3.0) < 1e-12 and abs(2.0 - 2.0) < 1e-12)

print("\n[TARGET] direct coefficient readout = (1,1.5,8) trivially")
target_norm = c / c[0]
check(3, "normalized coefficient readout equals target (1,1.5,8)",
      np.allclose(target_norm, [1.0, 1.5, 8.0], atol=1e-12))

print("\n[MOMENTS] mu_n = int x^n rho dx (corpus hierarchy object, P20)")
mu = np.array([moment(n) for n in range(4)])
# 20_Paper_StandardPhysicsEmbedding.tex line 133: mu_0 = alpha^-1 = 137.036
check(4, "mu_0 = alpha^-1 = 4pi^3+pi^2+pi (137.036)",
      abs(mu[0] - (4 * PI**3 + PI**2 + PI)) < 1e-9 and abs(mu[0] - 137.0363) < 1e-3)
# mu_1 = 108.717 (P20 line 136); the moment sequence is monotone DECREASING
check(5, "moments monotone-decreasing, mu_1~108.72; do NOT equal (2,3,16)",
      mu[0] > mu[1] > mu[2] > mu[3] and abs(mu[1] - 108.7167) < 1e-2)
mom_ratio = np.array([mu[1] / mu[0], mu[2] / mu[0], mu[3] / mu[0]])
check(6, "moment ratios mu_n/mu_0 are L2 >> 1 from target (not the home of (2,3,16))",
      np.linalg.norm(mom_ratio / mom_ratio[0] - target_norm) > 1.0)

print("\n[DRIFT] overlap-diagonal = A296's EXCLUDED bare Taylor probe")
x, h, lam, psi = build_eigs()
# eigenvalues ~16.51,51.24,102.02 (OPEN_HEARTS_STATUS line 12)
check(7, "production eigenvalues ~ (16.51, 51.24, 102.02)",
      abs(lam[0] - 16.51) < 0.1 and abs(lam[1] - 51.24) < 0.1 and abs(lam[2] - 102.02) < 0.2)
B = np.array([
    np.sum(psi[:, n]**2 * (c[n] * PI**(n + 1) * x**(n + 1))) * h / lam[n]
    for n in range(3)
])
Bn = B / B[0]
# A299 lines 49-50: normalized diagonal (1, 1.12, 7.50), L2 = 0.63
# the overlap-diagonal MISSES the target; mu entry 1.12 not 1.5.
# This is the DRIFT: the corpus excluded this functional (296,297),
# yet ~18 probes chased it. The defining numerical condition of
# FRAMING-DRIFT: readout L2=0 to target while overlap-diagonal L2>0.1.
overlap_L2 = np.linalg.norm(Bn - np.array([1.0, 1.5, 8.0]))
check(8, "overlap-diagonal (1,1.12,7.50) MISSES target (L2~0.63) while readout L2=0 => FRAMING-DRIFT",
      abs(Bn[1] - 1.12) < 0.03 and abs(overlap_L2 - 0.63) < 0.05
      and np.linalg.norm(target_norm - [1.0, 1.5, 8.0]) < 1e-9)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)