#!/usr/bin/env python3
"""verify_P318.py -- Verifier for Addendum 318 (OI-287-1 magnitude piece, the
first-order B^4 radial Dirac reduction).

Recomputes from scratch the radial Dirac system of A318 and asserts, by machine,
its defining numbers and the NULL verdict's defining condition. A318 tested the
entry A317 named: the family / phase entering through the FIRST-ORDER Dirac /
off-diagonal structure (the full B^4 spinor problem of P18 S2.2/S3.1) rather than
as the additive m^2 energy of A317. The construction:

  - the two-component radial Dirac operator on (f,g), self-adjoint in r^3 dr, with
    the off-diagonal first-derivative coupling -d/dr - (kappa+c)/r and its
    adjoint, c = 3/2 the 4D dimension constant (A312/A314: c = Delta_psi = R/4);
  - the canonical family map: family n -> n-th boundary-Dirac sector
    kappa_n = (n-1)+3/2 = 1.5, 2.5, 3.5 (A314 tower +-(j+3/2), the omega^k / Z_3
    spinor sector);
  - the density overlap as the spinor bilinear f^2+g^2 in r^3 dr, divided by |E|;
  - the central comparison: is the first-order coupling gentler than A317's m^2
    overshoot (mu = 2.29-2.43)?

Honest scope: this machine-verifies the radial Dirac construction's numbers and
the NULL verdict. It does NOT close OI-287-1; the magnitude piece remains open.

  S1  the radial Dirac system is stable and gives three distinct families  - 1-3
  S2  free-Dirac boundary sanity vs the A314 +-(j+3/2) tower ordering      - 4
  S3  the canonical map's diagonal / mu / tau / L2 / peak                  - 5-7
  S4  the central comparison and the NULL verdict's defining condition     - 8-9

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
C = (2.0, 3.0, 16.0)
TARGET = [c / C[0] for c in C]
N = 800
C_DIM = 1.5
CANONICAL_KAPPA = [1.5, 2.5, 3.5]

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


# --- recompute the radial Dirac construction from scratch -------------------
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 mass(x):
    return np.sqrt(Vpot(x))


def grids():
    dr = 1.0 / N
    rf = np.array([(i + 1) * dr for i in range(N - 1)])
    rg = np.array([(i + 0.5) * dr for i in range(N - 1)])
    return rf, rg, dr


def dirac_matrix(kappa, m_on_f, m_on_g):
    rf, rg, dr = grids()
    n = N - 1
    M = np.zeros((2 * n, 2 * n))
    for i in range(n):
        M[i, i] = m_on_f[i]
        M[n + i, n + i] = -m_on_g[i]
    for i in range(n):
        c_kp = (kappa + C_DIM) / rf[i]
        coupl_right = -(1.0 / dr) - 0.5 * c_kp
        coupl_left = (1.0 / dr) - 0.5 * c_kp
        M[i, n + i] += coupl_right
        if i - 1 >= 0:
            M[i, n + (i - 1)] += coupl_left
    diag_keep = np.diag(M).copy()
    M = M + M.T
    for i in range(2 * n):
        M[i, i] = diag_keep[i]
    return M, rf, rg, dr


def solve_dirac(kappa, free=False):
    rf, rg, dr = grids()
    if free:
        m_on_f = np.zeros(N - 1)
        m_on_g = np.zeros(N - 1)
    else:
        m_on_f = mass(rf)
        m_on_g = mass(rg)
    M, rf, rg, dr = dirac_matrix(kappa, m_on_f, m_on_g)
    w, v = np.linalg.eigh(M)
    idx = int(np.argmin(np.abs(w)))
    E = float(w[idx])
    vec = v[:, idx]
    n = N - 1
    f = vec[:n]
    g = vec[n:]
    g_on_f = np.empty(n)
    g_on_f[1:] = 0.5 * (g[1:] + g[:-1])
    g_on_f[0] = g[0]
    return f, g_on_f, E, rf, dr, w


def overlap(f, g, rf, dr, E):
    w3 = rf**3
    d = f**2 + g**2
    norm = float(np.sum(d * w3) * dr)
    return [float(np.sum(d * C[k - 1] * PI**k * rf**k * w3) * dr) / norm / abs(E)
            for k in (1, 2, 3)]


def family_matrix():
    Es = []
    rows = []
    for n in range(3):
        f, g, E, rf, dr, _ = solve_dirac(CANONICAL_KAPPA[n])
        Es.append(E)
        rows.append(overlap(f, g, rf, dr, E))
    M = np.array(rows)
    return M, Es


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(rf, dr):
    psi = np.column_stack([math.sqrt(2.0) * np.sin((n + 1) * PI * rf)
                           for n in range(3)])
    lam = np.array([(n + 1)**2 * PI**2 for n in range(3)])
    M = np.zeros((3, 3))
    w3 = rf**3
    for n in range(3):
        d = psi[:, n]**2
        nm = float(np.sum(d * w3) * dr)
        for k in (1, 2, 3):
            M[n, k - 1] = float(np.sum(d * C[k - 1] * PI**k * rf**k * w3)
                                * dr) / nm / lam[n]
    return M


# ===========================================================================
print("S1  Radial Dirac system: stable, three genuinely distinct families")

M, Es = family_matrix()
rf, rg, dr = grids()
ratio = norm_diag(M)
no_nan = bool(np.all(np.isfinite(M)) and all(math.isfinite(e) for e in Es))
check(1, "radial Dirac eigenproblem solves with no NaN/Inf (stable staggered "
      "discretization)", no_nan)

distinct = (abs(Es[0] - Es[1]) > 1e-4 and abs(Es[0] - Es[2]) > 1e-4
            and abs(Es[1] - Es[2]) > 1e-4)
check(2, "three families (kappa = 1.5, 2.5, 3.5) give three DISTINCT Dirac "
      "eigenvalues E = (%.3f, %.3f, %.3f): A316 degeneracy broken by the "
      "first-order structure" % (Es[0], Es[1], Es[2]), distinct)

diag_entries = [M[n, n] for n in range(3)]
three_distinct = (abs(diag_entries[0] - diag_entries[1]) > 1e-6
                  and abs(diag_entries[0] - diag_entries[2]) > 1e-6
                  and abs(diag_entries[1] - diag_entries[2]) > 1e-6)
check(3, "three DISTINCT diagonal overlap entries (no degenerate collapse)",
      three_distinct)

# ===========================================================================
print("S2  Free-Dirac boundary sanity vs the A314 +-(j+3/2) tower ordering")

# A314: the genuine S^3 boundary-Dirac spectrum is +-(j+3/2), monotone increasing
# in j. The free (m=0) radial Dirac operator's lowest |E| mode per sector must be
# MONOTONE INCREASING in kappa (higher kappa = higher j = higher boundary mode),
# the corpus-faithful signature of the +-(j+3/2) tower carried into the 4-ball
# radial spectrum (the radial kinetic tower adds on top of the boundary kappa).
free_lowest = []
for kap in CANONICAL_KAPPA:
    _, _, _, _, _, w = solve_dirac(kap, free=True)
    pos = np.sort(w[w > 1e-9])
    free_lowest.append(float(pos[0]))
mono = free_lowest[0] < free_lowest[1] < free_lowest[2]
check(4, "free radial Dirac lowest |E|+ is monotone increasing in kappa "
      "(%.3f < %.3f < %.3f): the +-(j+3/2) tower ordering of A314 carried into "
      "the 4-ball radial spectrum" % tuple(free_lowest), mono)

# ===========================================================================
print("S3  Canonical map: diagonal / mu / tau / L2 / peak")

mu = ratio[1]
tau = ratio[2]
L2 = l2(ratio)
check(5, "canonical diagonal normalized to (1, %.3f, %.3f) with mu=%.3f, "
      "tau=%.3f (overshoots the target 1.5/8)" % (mu, tau, mu, tau),
      abs(ratio[0] - 1.0) < 1e-9 and mu > 1.7 and tau > 8.0)

peak = diag_peaks(M)
box_peak = diag_peaks(box_matrix(rf, dr))
valid_peak = peak and not box_peak
check(6, "no valid diagonal peak restored (peak=%s, box_control_peak=%s, "
      "valid_peak=%s)" % (peak, box_peak, valid_peak), not valid_peak)

within20 = all(abs(ratio[i] - TARGET[i]) <= 0.20 * TARGET[i] for i in range(3))
check(7, "canonical map does NOT land within 20%% of (1,1.5,8) "
      "(L2=%.3f, within20=%s): no HIT" % (L2, within20),
      (not within20) and L2 > 0.5)

# ===========================================================================
print("S4  Central comparison (first-order Dirac vs A317 m^2) + NULL verdict")

A317_overshoot_mu = 2.29        # F-Hopf, the closest of A317's two m^2 forms
SCALAR_closest_L2 = 0.399
SCALAR_closest_mu = 1.222

gentler = mu < A317_overshoot_mu
# the central honest comparison: the first-order coupling did NOT come down
# gentler than A317's m^2; mu stays at/above the A317 overshoot.
check(8, "CENTRAL COMPARISON: first-order Dirac coupling is NOT gentler than "
      "A317's m^2 overshoot -- mu=%.3f vs A317's 2.29-2.43; mu did NOT come "
      "down toward 1.5" % mu, (not gentler) and mu >= A317_overshoot_mu)

improves_scalar = L2 < SCALAR_closest_L2 - 1e-9
moves_mu = abs(mu - 1.5) < abs(SCALAR_closest_mu - 1.5) - 1e-9
# NULL verdict's defining condition: does NOT improve on the scalar-radial
# closest (L2=0.399, mu=1.222) and is not gentler than A317 -> NULL.
null_holds = (not improves_scalar) and (not moves_mu) and (not gentler)
check(9, "NULL verdict's defining condition holds: the Dirac reduction does "
      "NOT improve on the scalar-radial closest (L2=%.3f vs 0.399; "
      "improves=%s) and does not move mu toward 1.5 (moves=%s) -- magnitude "
      "not in this reduction" % (L2, improves_scalar, moves_mu), null_holds)

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