#!/usr/bin/env python3
"""verify_P305.py -- Verifier for Addendum 305 (non-separable joint kernel).

Asserts: (S1) the eigenstructure reproduces -- eigenvalues ordered, node
counts (0,1,2), and the A299 base diagonal reproduces (1, 1.12, 7.50);
(S2) each declared joint kernel (J1 Gaussian, J2 band, J3a self-structure,
J4 hard projection, J5 separable control) leaves the normalized diagonal
unchanged at (1, 1.12, 7.50) -- kernels equal 1 on the matching pair, so
off-diagonal suppression cannot move a diagonal entry -- and no kernel
yields a valid diagonal peak; (S3) the NULL verdict's defining condition
holds -- no candidate both peaks and lands within 20% of (1,1.5,8), the
closest non-control sits at L2~0.63 (the A299 control distance), mu did
NOT move toward 1.5, and the separable control J5 does not beat the joint
kernels, so the non-separable joint kernel is a ninth excluded class.

  S1  Eigenstructure + A299 base reproduce  - checks 1-3
  S2  Kernel diagonals / no valid peak      - checks 4-6
  S3  NULL verdict conditions hold          - checks 7-8
"""
import sys
import math

import numpy as np

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
C = (2.0, 3.0, 16.0)
TGT = [c / C[0] for c in C]
N = 1500


def eigfuns(potential):
    dx = 1.0 / N
    x = np.linspace(dx, 1 - dx, N - 1)
    off = -1.0 / dx**2 * np.ones(N - 2)
    H = (np.diag(2.0 / dx**2 + potential(x))
         + np.diag(off, 1) + np.diag(off, -1))
    w, v = np.linalg.eigh(H)
    return x, w[:3], v[:, :3] / math.sqrt(dx)


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


def densM(p, x, dx):
    return np.array([[float(np.sum(p[:, n]**2 * C[k - 1] * PI**k * x**k) * dx)
                      for k in (1, 2, 3)] for n in range(3)])


def overlapM(p, x, dx):
    return np.array([[float(np.sum(p[:, n]**2 * PI**k * x**k) * dx)
                      for k in (1, 2, 3)] for n in range(3)])


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(r):
    return math.sqrt(sum((r[i] - TGT[i])**2 for i in range(3)))


rho_p = lambda x: 48 * PI**3 * x**2 + 6 * PI**2 * x + 2 * PI
Vpot = lambda x: rho_p(x)**2 / (2 * OM**2) + ES * x**2 * (1 - x)**2
x, lam, psi = eigfuns(Vpot)
dx = x[1] - x[0]
nu = [node_count(psi[:, n]) for n in range(3)]

B = densM(psi, x, dx) / np.array(lam)[:, None]
r_base = norm_diag(B)

OV = overlapM(psi, x, dx)
OV_hat = OV / np.array([OV[n, n] for n in range(3)])[:, None]

psi_box = np.column_stack([math.sqrt(2.0) * np.sin((n + 1) * PI * x)
                           for n in range(3)])
lam_box = np.array([(n + 1)**2 * PI**2 for n in range(3)])
B_box = densM(psi_box, x, dx) / lam_box[:, None]
OV_box = overlapM(psi_box, x, dx)
OV_box_hat = OV_box / np.array([OV_box[n, n] for n in range(3)])[:, None]

Kg = np.array([[math.exp(-((k - (n + 1))**2) / 2.0)
                for k in (1, 2, 3)] for n in range(3)])
Kb = np.array([[0.5**abs(k - (n + 1))
                for k in (1, 2, 3)] for n in range(3)])
Khard = np.eye(3)
f_sep = np.array([1 + nu[n] for n in range(3)], dtype=float)
g_sep = np.array([1.0 / k for k in (1, 2, 3)])
Ksep = f_sep[:, None] * g_sep[None, :]

kernels = {
    "J1": (Kg, Kg),
    "J2": (Kb, Kb),
    "J3a": (OV_hat, OV_box_hat),
    "J4": (Khard, Khard),
    "J5": (Ksep, Ksep),
}
res = {}
for name, (K, Kbx) in kernels.items():
    M = B * K
    Mb = B_box * Kbx
    r = norm_diag(M)
    valid = diag_peaks(M) and not diag_peaks(Mb)
    within = all(abs(r[i] - TGT[i]) <= 0.20 * TGT[i] for i in range(3))
    res[name] = {"ratio": r, "valid_peak": valid, "L2": l2(r),
                 "within20": within}

# closest non-control = all of J1..J5 here (J0 control not in dict)
closest = min(res, key=lambda k: res[k]["L2"])
closest_mu = res[closest]["ratio"][1]
mu_moved_toward = abs(closest_mu - 1.5) < abs(r_base[1] - 1.5)
any_valid = any(res[k]["valid_peak"] for k in res)
hit = any(res[k]["valid_peak"] and res[k]["within20"] for k in res)

joint = {k: res[k] for k in ("J1", "J2", "J3a", "J4")}
best_joint_L2 = min(v["L2"] for v in joint.values())
sep_L2 = res["J5"]["L2"]
any_joint_valid = any(v["valid_peak"] for v in joint.values())
sep_valid = res["J5"]["valid_peak"]
sep_underperforms = (sep_L2 >= best_joint_L2) and not (
    sep_valid and not any_joint_valid)

print("S1  Eigenstructure + A299 base reproduce")
check(1, "eigenvalues ordered (%.2f, %.2f, %.2f)" % tuple(lam),
      lam[0] < lam[1] < lam[2])
check(2, "node counts are (0,1,2)", nu == [0, 1, 2])
check(3, "A299 base diagonal is (1, %.2f, %.2f), within 0.02 of (1,1.12,7.50)"
      % (r_base[1], r_base[2]),
      abs(r_base[1] - 1.12) < 0.02 and abs(r_base[2] - 7.50) < 0.02)

print("S2  Kernel diagonals and no valid peak")
check(4, "J1 Gaussian diag (1, %.2f, %.2f) equals A299 base (kernel=1 on pair)"
      % (res["J1"]["ratio"][1], res["J1"]["ratio"][2]),
      abs(res["J1"]["ratio"][1] - r_base[1]) < 1e-9
      and abs(res["J1"]["ratio"][2] - r_base[2]) < 1e-9)
check(5, "J2 band diag (1, %.2f, %.2f) equals A299 base (off-diagonal only)"
      % (res["J2"]["ratio"][1], res["J2"]["ratio"][2]),
      abs(res["J2"]["ratio"][1] - r_base[1]) < 1e-9
      and abs(res["J2"]["ratio"][2] - r_base[2]) < 1e-9)
check(6, "no declared joint kernel yields a valid diagonal peak",
      not any_valid)

print("S3  NULL verdict conditions hold")
check(7, "closest is %s at L2=%.2f (= A299 control L2=0.63); mu did NOT move "
      "toward 1.5 (closest mu=%.2f vs base %.2f)"
      % (closest, res[closest]["L2"], closest_mu, r_base[1]),
      abs(res[closest]["L2"] - 0.63) < 0.02 and not mu_moved_toward)
check(8, "no HIT and separable control J5 does not beat the joint kernels",
      (not hit) and sep_underperforms)

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