#!/usr/bin/env python3
"""verify_P303.py -- Verifier for Addendum 303 (degree-dependent factor).

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 degree-dependent candidate (D1 factorial, D2 inverse
degree, D3 pi-power strength, D5 linear strength) yields the diagonal,
peak status, and L2 the probe found, and none 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 is
D2 with L2~5.58, and its mu entry did NOT move toward 1.5 (moved away),
so the degree-dependent weight is an eighth excluded mechanism class.

  S1  Eigenstructure + A299 base reproduce  - checks 1-3
  S2  Candidate 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
KAPPA = OM ** (-1.25)
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 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)

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]

cand = {
    "D1": [1.0 / math.factorial(k) for k in (1, 2, 3)],
    "D2": [1.0 / k for k in (1, 2, 3)],
    "D3": [PI**(-k) * KAPPA * OM for k in (1, 2, 3)],
    "D5": [KAPPA * OM * k for k in (1, 2, 3)],
}
res = {}
for name, w in cand.items():
    wv = np.array(w)
    M = B * wv[None, :]
    Mb = B_box * wv[None, :]
    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 = 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)

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  Candidate diagonals and no valid peak")
# D1 factorial: diag (1, 0.56, 1.25); D2 inverse: (1, 0.56, 2.50)
check(4, "D1 factorial diag (1, %.2f, %.2f) within 0.02 of (1,0.56,1.25)"
      % (res["D1"]["ratio"][1], res["D1"]["ratio"][2]),
      abs(res["D1"]["ratio"][1] - 0.56) < 0.02
      and abs(res["D1"]["ratio"][2] - 1.25) < 0.02)
check(5, "D2 inverse-degree diag (1, %.2f, %.2f) within 0.02 of (1,0.56,2.50)"
      % (res["D2"]["ratio"][1], res["D2"]["ratio"][2]),
      abs(res["D2"]["ratio"][1] - 0.56) < 0.02
      and abs(res["D2"]["ratio"][2] - 2.50) < 0.02)
check(6, "no declared candidate yields a valid diagonal peak",
      not any_valid)

print("S3  NULL verdict conditions hold")
check(7, "closest non-control is D2 with L2=%.2f (> A299 L2=0.63)"
      % res["D2"]["L2"],
      closest == "D2" and res["D2"]["L2"] > 0.63)
check(8, "no HIT and mu did NOT move toward 1.5 (closest mu=%.2f vs A299 %.2f)"
      % (closest_mu, r_base[1]),
      (not hit) and (not mu_moved_toward))

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