#!/usr/bin/env python3
"""verify_P298.py — Verifier for Addendum 298 (the nodal direction).

Asserts: (S1) no fixed nodal/derivative pairing makes the coupling
peak on the diagonal; (S2) the derivative/nodal functionals recover
the increasing (2,3,16) direction the A297 position moments reversed.

  S1  Fifth class excluded     - checks 1-4
  S2  Direction recovered      - checks 5-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


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


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


def increasing(M):
    r = ratio(M)
    return r[0] < r[1] < r[2]


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, evals, psi = eigfuns(Vpot)
dx = x[1] - x[0]

# position moments (A297 reference)
Apos = np.array([[float(np.sum(psi[:, n]**2 * x**k) * dx)
                  for k in (1, 2, 3)] for n in range(3)])
# P1 derivative overlap
d = np.gradient(psi, dx, axis=0)
M1 = np.array([[float(np.sum(d[:, n]**2 * x**k) * dx)
                for k in (1, 2, 3)] for n in range(3)])
# P2 kinetic coupling
def kin(n, k):
    f = x**k * psi[:, n]
    d2 = np.gradient(np.gradient(f, dx), dx)
    return float(np.sum(psi[:, n] * (-d2)) * dx)
M2 = np.array([[kin(n, k) for k in (1, 2, 3)] for n in range(3)])
# P3 density-weighted overlap
M3 = np.array([[float(np.sum(psi[:, n]**2 * C[k - 1] * PI**k * x**k)
                      * dx) for k in (1, 2, 3)] for n in range(3)])

print("S1  Fifth class excluded")
check(1, "eigenvalues ordered (%.2f, %.2f, %.2f)" % tuple(evals),
      evals[0] < evals[1] < evals[2])
check(2, "derivative overlap: no diagonal peak", not diag_pref(M1))
check(3, "kinetic coupling: no diagonal peak", not diag_pref(M2))
check(4, "density overlap: no diagonal peak", not diag_pref(M3))

print("S2  Direction recovered")
check(5, "position moments DECREASE (1, %.2f, %.2f), wrong sign"
      % (ratio(Apos)[1], ratio(Apos)[2]), not increasing(Apos))
check(6, "derivative overlap INCREASES (1, %.2f, %.2f)"
      % (ratio(M1)[1], ratio(M1)[2]), increasing(M1))
check(7, "kinetic coupling INCREASES (1, %.2f, %.2f)"
      % (ratio(M2)[1], ratio(M2)[2]), increasing(M2))
check(8, "density overlap INCREASES with degree-3 jump (1, %.2f, %.2f); "
         "target (1,1.5,8)" % (ratio(M3)[1], ratio(M3)[2]),
      increasing(M3) and ratio(M3)[2] > 8 * ratio(M3)[1] / 3)

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