#!/usr/bin/env python3
"""verify_P315.py -- Verifier for Addendum 315 (angular-coupling scan, OI-287-1).

Recomputes the one-parameter centrifugal-scaling families from scratch and
asserts the pre-registered conditions of A315. A315 asks whether the magnitude
residual to mu = 1.5 (A311/A313: the forced r^3 arena measure lifts mu to 1.34,
the full l=n angular identity overshoots to 2.009) is a ONE-PARAMETER angular
coupling question. It scales the A313 centrifugal term g*l(l+2)/r^2, l = n, by a
scalar g in [0,1] holding the construction fixed, in two mode-consistent readings:
  FAMILY A (nth)   : family n = n-th mode of O_{l=n}(g); g=0 -> A313 C_4Dself.
  FAMILY B (ground): family n = ground mode of O_{l=n}(g); g=1 -> A313 C_ang.

The honest finding is NULL with a structural reason: A313's two endpoints
(mu=1.222 C_4Dself, mu=2.009 C_ang) differ in MODE-SELECTION (modes 0,1,2 of one
operator vs the ground mode of each l-sector), not in centrifugal strength.
Holding the construction fixed and scaling only g, family A runs mu DOWN
1.222->0.975 and family B runs mu DOWN 2.532->2.009; NEITHER crosses 1.5. There
is no g*, so there is nothing to compare against the pre-declared corpus
constants. The verifier encodes this honest verdict: it asserts no single-knob
family crosses mu = 1.5, and reports the closest corpus constant to each family's
nearest-to-1.5 g for completeness while asserting (honestly) that no such g
exists as a genuine crossing.

  S1  A313 endpoint reproduction + family A behavior - checks 1-3
  S2  Family B behavior + the no-crossing structural fact - checks 4-6
  S3  Verdict NULL + corpus-constant comparison is vacuous + no NaN - checks 7-8
"""
import sys
import math

import numpy as np


def gen_eigh(A, Bm):
    L = np.linalg.cholesky(Bm)
    Linv = np.linalg.inv(L)
    Cmat = Linv @ A @ Linv.T
    Cmat = 0.5 * (Cmat + Cmat.T)
    w, y = np.linalg.eigh(Cmat)
    return w, Linv.T @ y


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}")


PI = math.pi
OM = 4 * PI**3 + PI**2 + PI
ES = 13.177
GAMMA = 0.75
C = (2.0, 3.0, 16.0)
TGT = [c / C[0] for c in C]
N = 1500

ALPHA = 1.0 / 137.036
CANDIDATE_CONSTANTS = {
    "FRAC_EDGE": PI / OM,
    "FRAC_BOUNDARY": PI**2 / OM,
    "FRAC_BULK": 4 * PI**3 / OM,
    "gamma_3_4": GAMMA,
    "one_half": 0.5,
    "one_third": 1.0 / 3.0,
    "two_thirds": 2.0 / 3.0,
    "alpha": ALPHA,
    "alpha_quarter": ALPHA**0.25,
    "inv_pi": 1.0 / PI,
}
HIT_TOL_REL = 0.02


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


def Vpot(x):
    return rho_p(x)**2 / (2 * OM**2) + ES * x**2 * (1 - x)**2


def grid():
    dx = 1.0 / N
    return np.linspace(dx, 1 - dx, N - 1), dx


def radial_matrices(x, dx, V_on_grid, weight):
    n = N - 1
    w_mid = 0.5 * (weight[:-1] + weight[1:])
    w_mid_left0 = max(0.5 * (0.0 + weight[0]), 0.0)
    A = np.zeros((n, n))
    for i in range(n):
        wl = w_mid[i - 1] if i > 0 else w_mid_left0
        wr = w_mid[i] if i + 1 < n else weight[i]
        A[i, i] = (wl + wr) / dx**2
        if i + 1 < n:
            A[i, i + 1] = -w_mid[i] / dx**2
            A[i + 1, i] = -w_mid[i] / dx**2
    A += np.diag(weight * V_on_grid)
    A = 0.5 * (A + A.T)
    return A, np.diag(weight)


def overlap_row(col, x, dx, w_arr, lam):
    """Normalized degree-k overlap row for one mode (n-th col), divided by lam."""
    norm = float(np.sum(col**2 * w_arr) * dx)
    return [float(np.sum(col**2 * C[k - 1] * PI**k * x**k * w_arr) * dx)
            / norm / lam for k in (1, 2, 3)]


def mu_both_families(g):
    """Solve O_{l=n}(g) ONCE per n; extract BOTH the n-th mode (family A) and the
    ground mode (family B) from the same eigendecomposition. Returns
    (mu_nth, mu_ground) for coupling g (overlap and states in r^3 dr)."""
    x, dx = grid()
    w_arr = x**3
    M_nth = np.zeros((3, 3))
    M_grd = np.zeros((3, 3))
    for n in range(3):
        cent = g * n * (n + 2) / x**2
        A, Wm = radial_matrices(x, dx, Vpot(x) + cent, w_arr)
        w, v = gen_eigh(A, Wm)
        o = np.argsort(w)
        col_nth = v[:, o[n]]
        col_grd = v[:, o[0]]
        M_nth[n, :] = overlap_row(col_nth, x, dx, w_arr, w[o[n]])
        M_grd[n, :] = overlap_row(col_grd, x, dx, w_arr, w[o[0]])
    mu_nth = M_nth[1, 1] / M_nth[0, 0]
    mu_grd = M_grd[1, 1] / M_grd[0, 0]
    return mu_nth, mu_grd


# ---- recompute -----------------------------------------------------------
gridg = [round(0.1 * i, 1) for i in range(11)]
both = [mu_both_families(g) for g in gridg]
mus_nth = [b[0] for b in both]
mus_grd = [b[1] for b in both]

nth_g0 = mus_nth[0]
nth_g1 = mus_nth[-1]
grd_g0 = mus_grd[0]
grd_g1 = mus_grd[-1]

nth_mono_dec = all(mus_nth[i + 1] <= mus_nth[i] + 1e-9
                   for i in range(len(mus_nth) - 1))
grd_mono_dec = all(mus_grd[i + 1] <= mus_grd[i] + 1e-9
                   for i in range(len(mus_grd) - 1))

nth_crosses = min(mus_nth) <= 1.5 <= max(mus_nth)
grd_crosses = min(mus_grd) <= 1.5 <= max(mus_grd)
any_cross = nth_crosses or grd_crosses

# nan guard
allmu = mus_nth + mus_grd
nan_found = any(m != m for m in allmu)

# For honesty completeness: the g in each family nearest to mu=1.5, and the
# corpus constant nearest to that g (this is NOT a crossing; reported only to
# state the comparison explicitly and assert it is vacuous as a derivation).
def nearest_g_to_15(mus):
    return gridg[min(range(len(mus)), key=lambda i: abs(mus[i] - 1.5))]


g_near_nth = nearest_g_to_15(mus_nth)   # family A best approach to 1.5
g_near_grd = nearest_g_to_15(mus_grd)   # family B best approach to 1.5


def closest_constant(g):
    name = min(CANDIDATE_CONSTANTS,
               key=lambda n: abs(g - CANDIDATE_CONSTANTS[n]))
    rel = abs(g - CANDIDATE_CONSTANTS[name]) / abs(CANDIDATE_CONSTANTS[name])
    return name, rel


cc_nth, rel_nth = closest_constant(g_near_nth)

verdict = "NULL" if not any_cross else "FIT_FLAGGED_PARTIAL"

# ---- assertions ----------------------------------------------------------
print("S1  A313 endpoint reproduction + family A (n-th mode) behavior")
check(1, "family A (n-th mode) at g=0 reproduces A313 C_4Dself mu=%.4f ~ 1.222"
      % nth_g0,
      abs(nth_g0 - 1.222) < 0.01)
check(2, "family B (ground mode) at g=1 reproduces A313 C_ang mu=%.4f ~ 2.009"
      % grd_g1,
      abs(grd_g1 - 2.009) < 0.01)
check(3, "family A runs mu monotonically DOWN from %.4f (g=0) to %.4f (g=1), "
      "never reaching 1.5 (it stays below the target)"
      % (nth_g0, nth_g1),
      nth_mono_dec and max(mus_nth) < 1.5 and not nth_crosses)

print("S2  Family B behavior + the no-crossing structural fact")
check(4, "family B runs mu monotonically DOWN from %.4f (g=0, degenerate) to "
      "%.4f (g=1), never reaching 1.5 (it stays above the target)"
      % (grd_g0, grd_g1),
      grd_mono_dec and min(mus_grd) > 1.5 and not grd_crosses)
check(5, "NO mode-consistent single-knob centrifugal family crosses mu=1.5 on "
      "[0,1]: family A max %.4f < 1.5 < family B min %.4f (mu=1.5 lies in the "
      "GAP between the two families, unreachable by either knob)"
      % (max(mus_nth), min(mus_grd)),
      (not any_cross) and max(mus_nth) < 1.5 < min(mus_grd))
check(6, "the A313 1.222->2.009 motion is MODE-REINDEXING, not coupling "
      "strength: family A (modes 0,1,2) and family B (ground of each l-sector) "
      "are disjoint mu-ranges [%.3f,%.3f] and [%.3f,%.3f] that do not overlap"
      % (min(mus_nth), max(mus_nth), min(mus_grd), max(mus_grd)),
      max(mus_nth) < min(mus_grd))

print("S3  Verdict NULL + corpus-constant comparison vacuous + no NaN")
check(7, "VERDICT NULL: no g* exists (no crossing), so there is NO g to test "
      "against the %d pre-declared corpus constants; the FIT/FORCED test is "
      "vacuous (closest constant to family A's best-approach g=%.2f would be "
      "%s at rel %.2f%%, but g=%.2f is NOT a mu=1.5 crossing)"
      % (len(CANDIDATE_CONSTANTS), g_near_nth, cc_nth, rel_nth * 100,
         g_near_nth),
      verdict == "NULL" and not any_cross)
check(8, "clean: no NaN across both family scans; A313 endpoints both "
      "reproduced (C_4Dself@g0=%.4f, C_ang@g1=%.4f); the honest verdict is "
      "NULL (the magnitude residual is structural, not a one-parameter "
      "coupling-strength question)"
      % (nth_g0, grd_g1),
      (not nan_found) and abs(nth_g0 - 1.222) < 0.01
      and abs(grd_g1 - 2.009) < 0.01 and verdict == "NULL")

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