#!/usr/bin/env python3
"""verify_P320.py -- Verifier for Addendum 320 (OI-287-1, the coupled
bulk+boundary B^4 problem).

Recomputes from scratch and asserts, by machine:

  S1  The coupled operator O_{l,k} (P18 S4.3) is the genuine corpus object: the
      radial operator with BOTH the boundary-harmonic centrifugal term l(l+2)/r^2
      AND the layer-cycle family term omega^k gamma (gamma = 3/4), and the two
      controls reproduce the two limits the coupling should interpolate between --
      C3a pure radial recovers the A311 undershoot mu = 1.222, C3b pure boundary
      recovers the A319 boundary overshoot (mu >= 2.0). Sanity.

  S2  The full coupled spectrum lambda_{n,L} = (radial eigenvalue of the l=L
      sector) + L(L+2) (boundary Casimir) + Re(omega^k) gamma is positive, and the
      boundary Casimir L(L+2) = (0, 3, 8) is the genuine A314/A319 datum.

  S3  Each coupled map's diagonal / mu / tau / peak / L2, recomputed from scratch,
      matches the probe; whether the coupled problem reaches the GAP [1.34, 2.0] no
      single reduction could, and the HIT / PARTIAL / NULL verdict's defining
      condition, asserted honestly.

The coupled problem OVERSHOOTS: the canonical C1 (corpus O_{l,k}, family n ->
degree L=n-1 ground) gives mu = 7.77, the linked C2 gives mu = 3.90, both ABOVE the
[1.34, 2.0] gap. The coupling does NOT interpolate the radial undershoot and the
boundary overshoot into the gap; the boundary multiplicity factor compounds the
overshoot. NULL: the magnitude is in NONE of radial, boundary, or their coupling;
the framing is the likely issue (companion probe A321).

Copyright Léon Fernando Vlegels -- CC BY 4.0
"""
import sys
import math

import numpy as np

PI = math.pi
OMEGA = 4 * PI**3 + PI**2 + PI
E_SELF = 13.177
GAMMA = 0.75                              # P18 S4.2: gamma = 3/4
C = (2.0, 3.0, 16.0)
TARGET = [c / C[0] for c in C]           # (1, 1.5, 8)
RADIAL_BEST_L2 = 0.399
GAP_LO, GAP_HI = 1.34, 2.0
N = 1000
FAMILY_L = [0, 1, 2]
FAMILY_K = [0, 1, 2]
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}")


# --- independent rebuild of the coupled probe --------------------------------
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


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 * OMEGA**2) + E_SELF * x**2 * (1 - x)**2


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


def radial4d_matrices(x, dx, V_on_grid, weight):
    n = N - 1
    w_mid = 0.5 * (weight[:-1] + weight[1:])
    w_mid_left0 = max(0.5 * 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)
    return 0.5 * (A + A.T), np.diag(weight)


def solve_coupled_sector(L, n_modes=3):
    x, dx = grid()
    w_arr = x**3
    centrifugal = L * (L + 2) / x**2 if L > 0 else np.zeros_like(x)
    A, Wm = radial4d_matrices(x, dx, Vpot(x) + centrifugal, w_arr)
    w, v = gen_eigh(A, Wm)
    order = np.argsort(w)
    w, v = w[order], v[:, order]
    vv = v[:, :n_modes].copy()
    for k in range(n_modes):
        nrm = math.sqrt(float(np.sum(vv[:, k]**2 * w_arr) * dx))
        vv[:, k] = vv[:, k] / nrm
    return x, dx, w[:n_modes], vv


def harm_factor(L):
    return float((L + 1)**2)


def coupled_ev(radial_ev, L, k):
    return radial_ev + L * (L + 2) + math.cos(2 * PI * k / 3.0) * GAMMA


def overlap_row(col, x, dx, w_arr, lam, hf):
    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 * hf / lam for k in (1, 2, 3)]


def node_count(col):
    s = np.sign(col[np.abs(col) > 1e-9])
    return int(np.sum(s[1:] != s[:-1]))


def diag_peaks(M):
    Z = np.zeros_like(M)
    for k in range(3):
        c = M[:, k]
        Z[:, k] = (c - c.mean()) / c.std() if c.std() > 0 else 0.0
    return all(int(np.argmax(Z[:, k])) == k for k in range(3))


def norm_diag(M):
    d = [M[n, n] for n in range(3)]
    return [v / d[0] for v in d]


def l2(ratio):
    return math.sqrt(sum((ratio[i] - TARGET[i])**2 for i in range(3)))


def box_matrix(x, dx, w_arr):
    M = np.zeros((3, 3))
    for n in range(3):
        col = math.sqrt(2.0) * np.sin((n + 1) * PI * x)
        lam = (n + 1)**2 * PI**2
        M[n, :] = overlap_row(col, x, dx, w_arr, lam, 1.0)
    return M


def build_C1(x, dx, w_arr):
    M = np.zeros((3, 3))
    lams = np.zeros(3)
    nodes = [0, 0, 0]
    for n in range(3):
        L, k = FAMILY_L[n], FAMILY_K[n]
        _, _, w_L, v_L = solve_coupled_sector(L, 3)
        lams[n] = coupled_ev(w_L[0], L, k)
        nodes[n] = node_count(v_L[:, 0])
        M[n, :] = overlap_row(v_L[:, 0], x, dx, w_arr, lams[n], harm_factor(L))
    return M, lams, nodes


def build_C2(x, dx, w_arr):
    M = np.zeros((3, 3))
    lams = np.zeros(3)
    nodes = [0, 0, 0]
    for n in range(3):
        L, k = FAMILY_L[n], FAMILY_K[n]
        _, _, w_L, v_L = solve_coupled_sector(L, 3)
        lams[n] = coupled_ev(w_L[n], L, k)
        nodes[n] = node_count(v_L[:, n])
        M[n, :] = overlap_row(v_L[:, n], x, dx, w_arr, lams[n], harm_factor(L))
    return M, lams, nodes


def build_C3a(x, dx, w_arr):
    _, _, w0, v0 = solve_coupled_sector(0, 3)
    M = np.zeros((3, 3))
    for n in range(3):
        M[n, :] = overlap_row(v0[:, n], x, dx, w_arr, w0[n], 1.0)
    return M


def build_C3b(x, dx, w_arr):
    _, _, w0, v0 = solve_coupled_sector(0, 1)
    ground = v0[:, 0]
    M = np.zeros((3, 3))
    for n in range(3):
        L = FAMILY_L[n]
        cas = L * (L + 2)
        w = cas if cas > 1e-9 else 1.0
        M[n, :] = overlap_row(ground, x, dx, w_arr, w, harm_factor(L))
    return M


# ===========================================================================
x, dx = grid()
w_arr = x**3
M_box = box_matrix(x, dx, w_arr)
box_peak = diag_peaks(M_box)

M_c1, lams_c1, nodes_c1 = build_C1(x, dx, w_arr)
M_c2, lams_c2, nodes_c2 = build_C2(x, dx, w_arr)
M_c3a = build_C3a(x, dx, w_arr)
M_c3b = build_C3b(x, dx, w_arr)

r_c1 = norm_diag(M_c1)
r_c2 = norm_diag(M_c2)
r_c3a = norm_diag(M_c3a)
r_c3b = norm_diag(M_c3b)

print("S1  Coupled operator O_{l,k}: the two controls recover the two limits")

# C3a pure radial recovers the A311 undershoot mu = 1.222
check(1, "C3a pure-radial control (L=0 all families) recovers the A311 radial "
      "undershoot: mu = %.3f ~ 1.222 (the undershoot limit, no valid peak)"
      % r_c3a[1],
      abs(r_c3a[1] - 1.222) < 0.02
      and not (diag_peaks(M_c3a) and not box_peak))

# C3b pure boundary recovers the A319 overshoot (mu >= 2.0, above the gap)
check(2, "C3b pure-boundary control (radial ground fixed, boundary Casimir + "
      "multiplicity) recovers the A319 boundary overshoot: mu = %.3f >= 2.0 "
      "(the overshoot limit, above the gap, no valid peak)" % r_c3b[1],
      r_c3b[1] >= GAP_HI
      and not (diag_peaks(M_c3b) and not box_peak))

print("S2  The full coupled spectrum lambda_{n,L} = radial + L(L+2) + "
      "Re(omega^k) gamma is positive; boundary Casimir is the A314 datum")

# coupled eigenvalues positive
lam_pos = all(float(v) > 0 for v in list(lams_c1) + list(lams_c2))
check(3, "all coupled eigenvalues lambda_{n,L} > 0 (C1 %s, C2 %s); the genuine "
      "coupled spectrum, no division pathology"
      % ([round(float(v), 2) for v in lams_c1],
         [round(float(v), 2) for v in lams_c2]),
      lam_pos)

# the boundary Casimir L(L+2) entering the coupled spectrum = (0, 3, 8)
casimir_ok = all(abs(FAMILY_L[n] * (FAMILY_L[n] + 2) - [0, 3, 8][n]) < 1e-9
                 for n in range(3))
# and it genuinely enters: C1 coupled lam = radial + L(L+2) + Re(omega^k)gamma,
# so lam[n] - lam[0]-style structure carries the Casimir. Check the L=1,2 lam
# exceed the L=0 lam by at least their Casimir (3, 8) minus the radial spread.
casimir_enters = (lams_c1[1] > lams_c1[0] and lams_c1[2] > lams_c1[1])
check(4, "boundary Casimir L(L+2) = (0, 3, 8) (A314 datum) enters the coupled "
      "spectrum: the l=1,2 coupled eigenvalues exceed the l=0 one (Casimir + "
      "radial centrifugal both raise the spectrum)",
      casimir_ok and casimir_enters)

print("S3  Coupled maps' diagonal / mu / tau / peak / L2; gap reach; verdict")

# C1 CANONICAL: overshoots HARD, mu ~ 7.77, no peak
mu1, tau1, L21 = r_c1[1], r_c1[2], l2(r_c1)
peak1 = diag_peaks(M_c1) and not box_peak
check(5, "C1 CANONICAL (corpus O_{l,k}, family n -> degree L=n-1 ground) "
      "diagonal (1, %.2f, %.1f): mu = %.2f (overshoot, above the gap), "
      "tau = %.1f, L2 = %.1f, no valid peak (%s)"
      % (mu1, tau1, mu1, tau1, L21, peak1),
      abs(mu1 - 7.7652) < 5e-2 and abs(tau1 - 152.33) < 1.0
      and not peak1 and abs(L21 - 144.46) < 1.0)

# C2 LINKED: overshoots, mu ~ 3.90, no peak
mu2, tau2, L22 = r_c2[1], r_c2[2], l2(r_c2)
peak2 = diag_peaks(M_c2) and not box_peak
check(6, "C2 LINKED (excited mode n-1 of l=n-1) diagonal (1, %.2f, %.1f): "
      "mu = %.2f (overshoot, above the gap), tau = %.1f, L2 = %.1f, no valid "
      "peak (%s)" % (mu2, tau2, mu2, tau2, L22, peak2),
      abs(mu2 - 3.903) < 5e-2 and abs(tau2 - 50.54) < 1.0
      and not peak2 and abs(L22 - 42.61) < 1.0)

# the GAP: neither coupled map reaches [1.34, 2.0]; both overshoot past it
reaches_gap = ((GAP_LO <= mu1 <= GAP_HI) or (GAP_LO <= mu2 <= GAP_HI))
check(7, "the coupled problem does NOT reach the [1.34, 2.0] gap no single "
      "reduction could (C1 mu = %.2f, C2 mu = %.2f both above 2.0): the coupling "
      "does NOT interpolate the radial undershoot and boundary overshoot into the "
      "gap; reaches_gap = %s" % (mu1, mu2, reaches_gap),
      not reaches_gap)

# the verdict's defining condition: NULL. Canonical C1 does not hit; coupled
# does not reach the gap and does not beat the standing-best L2 = 0.399.
canonical_hits = (all(abs(r_c1[i] - TARGET[i]) <= 0.20 * TARGET[i]
                      for i in range(3)) and peak1)
beats_radial = (L21 < RADIAL_BEST_L2 - 1e-9 or L22 < RADIAL_BEST_L2 - 1e-9)
nan_found = any(v != v for v in r_c1 + r_c2 + r_c3a + r_c3b)
verdict_null = (not canonical_hits) and (not reaches_gap) and (not beats_radial)
check(8, "NULL verdict's defining condition holds: canonical C1 (corpus O_{l,k}) "
      "does NOT hit, coupled does NOT reach the gap, coupled does NOT beat the "
      "standing-best L2 = 0.399; no NaN (magnitude in NONE of radial, boundary, "
      "coupling; framing is the likely issue)",
      verdict_null and not nan_found)

# extra honesty check: the controls genuinely bracket the gap (undershoot below,
# overshoot above) while no coupled map lands inside it -- the seven-arena
# pattern holds one level out.
brackets = (r_c3a[1] <= GAP_LO + 0.1 and r_c3b[1] >= GAP_HI - 0.1
            and not reaches_gap)
check(9, "the controls bracket the gap (C3a undershoots <= 1.34, C3b overshoots "
      ">= 2.0) yet NO coupled map lands inside [1.34, 2.0]: the undershoot/"
      "overshoot pattern holds for the coupled problem too",
      brackets)

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