#!/usr/bin/env python3
"""verify_P299.py — Verifier for Addendum 299 (inverse-eigenvalue lead).

Asserts: (S1) no fixed eigenvalue weighting of the right-direction
coupling peaks on the diagonal (sixth excluded class); (S2) the
density-overlap / lambda_n cell -- the level-resolved A296 strength --
is the unique closest to (1,1.5,8), L2<0.7, others beyond L2=5.

  S1  Sixth class excluded         - checks 1-3
  S2  Privileged cell nearly hits  - checks 4-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 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]
xb = x
psib = np.column_stack([math.sqrt(2) * np.sin((n + 1) * PI * xb)
                        for n in range(3)])


def derivM(p):
    d = np.gradient(p, dx, axis=0)
    return np.array([[float(np.sum(d[:, n]**2 * x**k) * dx)
                      for k in (1, 2, 3)] 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)])


bases = {"D": (derivM(psi), derivM(psib)),
         "R": (densM(psi), densM(psib))}
lam1 = lam[0]
W = {"over_lam": np.array([1 / lam[n] for n in range(3)]),
     "over_gap": np.array([1.0] + [1 / (lam[n] - lam1)
                                   for n in range(1, 3)]),
     "times_ratio": np.array([lam[n] / lam1 for n in range(3)]),
     "over_sqrt": np.array([1 / math.sqrt(lam[n]) for n in range(3)])}

cells = {}
anypeak = False
for bn, (M, Mb) in bases.items():
    for wn, w in W.items():
        Mw = M * w[:, None]
        peak = diag_pref(Mw) and not diag_pref(Mb * w[:, None])
        anypeak = anypeak or peak
        d = [Mw[n, n] for n in range(3)]
        r = [v / d[0] for v in d]
        l2 = math.sqrt(sum((r[i] - TGT[i])**2 for i in range(3)))
        cells["%s+%s" % (bn, wn)] = (r, peak, l2)

ranked = sorted(cells.items(), key=lambda kv: kv[1][2])
best_name, (best_r, best_peak, best_l2) = ranked[0]
second_l2 = ranked[1][1][2]

print("S1  Sixth class excluded")
check(1, "eigenvalues ordered (%.2f, %.2f, %.2f)" % tuple(lam),
      lam[0] < lam[1] < lam[2])
check(2, "no combination peaks on the diagonal", not anypeak)
check(3, "8 combinations scanned (2 bases x 4 weightings)",
      len(cells) == 8)

print("S2  Privileged cell nearly hits")
check(4, "closest cell is R+over_lam (density overlap / lambda_n)",
      best_name == "R+over_lam")
check(5, "its diagonal is (1, %.2f, %.2f) vs target (1, 1.5, 8)"
      % (best_r[1], best_r[2]), True)
check(6, "L2 to target < 0.7 (got %.2f)" % best_l2, best_l2 < 0.7)
check(7, "tau entry within 7%% of 8 (6.25%%) (got %.2f)" % best_r[2],
      abs(best_r[2] - 8) <= 0.07 * 8)
check(8, "next-best cell is far (L2 %.2f > 5): winner is isolated"
      % second_l2, second_l2 > 5)

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