#!/usr/bin/env python3
"""verify_P300.py — Verifier for Addendum 300 (node-count bracket).

Asserts: (S1) the eigenfunction node counts are (0,1,2), so level n is
distinguished by n-1 interior zeros; (S2) no declared node-count
operation on the A299 base peaks on the diagonal within 20%% (seventh
excluded class); (S3) the target IS bracketed by p=0 (A299, undershoot)
and p=1 (linear node gain, overshoot), but no single node-power (1+nu)^p
reconciles mu and tau -- mu needs p~0.42, tau needs p~0.06 -- so a
uniform node-power gain is excluded.

  S1  Node counts are (0,1,2)        - checks 1-2
  S2  Seventh class excluded         - checks 3-5
  S3  Bracketed, single power fails  - checks 6-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)
N = 1500
TGT = [c / C[0] for c in C]


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 diag_pref(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[n])) == n for n 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)]


def densM(p):
    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)])


B = densM(psi) / np.array(lam)[:, None]
g = np.array([1 + nu[n] for n in range(3)], dtype=float)
NC1 = B * g[:, None]
d0 = [B[n, n] for n in range(3)]
r0 = [v / d0[0] for v in d0]
d1 = [NC1[n, n] for n in range(3)]
r1 = [v / d1[0] for v in d1]
peak_any = diag_pref(NC1) or diag_pref(B * (g**2)[:, None])
p_mu = math.log(1.5 / r0[1]) / math.log(2)
p_tau = math.log(8.0 / r0[2]) / math.log(3)

print("S1  Node counts are (0,1,2)")
check(1, "eigenvalues ordered (%.2f, %.2f, %.2f)" % tuple(lam),
      lam[0] < lam[1] < lam[2])
check(2, "node counts are (0,1,2): level n has n-1 interior zeros",
      nu == [0, 1, 2])

print("S2  Seventh class excluded")
check(3, "A299 base diagonal is (1, %.2f, %.2f) (undershoots)"
      % (r0[1], r0[2]), abs(r0[1] - 1.12) < 0.02 and abs(r0[2] - 7.5) < 0.05)
check(4, "linear node gain diagonal (1, %.2f, %.2f) (overshoots)"
      % (r1[1], r1[2]), r1[1] > 1.5 and r1[2] > 8)
check(5, "no node-count op peaks on the diagonal", not peak_any)

print("S3  Bracketed, but single node-power fails")
check(6, "mu=1.5 bracketed by (%.2f, %.2f)" % (r0[1], r1[1]),
      r0[1] < 1.5 < r1[1])
check(7, "tau=8 bracketed by (%.2f, %.2f)" % (r0[2], r1[2]),
      r0[2] < 8.0 < r1[2])
check(8, "no single power reconciles: mu wants p=%.2f, tau wants p=%.2f"
      % (p_mu, p_tau), abs(p_mu - p_tau) > 0.2)

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