#!/usr/bin/env python3
"""verify_P297.py — Verifier for Addendum 297 (selection rule is not
in the amplitude). Asserts the pre-registered NULL: the bare
observation-operator eigenfunctions carry no index-matching in their
position moments, and the Dirichlet boundary forces unit leading
power at the center.

  S1  No diagonal preference        - checks 1-3
  S2  Linear onset at the center    - checks 4-6
  S3  Fourth class excluded         - 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
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 melts(x, psi):
    dx = x[1] - x[0]
    return np.array([[float(np.sum(psi[:, n]**2 * x**k) * dx)
                      for k in (1, 2, 3)] for n in range(3)])


def diag_pref(A):
    Z = np.zeros_like(A)
    for k in range(3):
        c = A[:, k]
        Z[:, k] = (c - c.mean()) / c.std() if c.std() > 0 else 0.0
    argmax = [int(np.argmax(Z[n])) for n in range(3)]
    return all(argmax[n] == n for n in range(3)), argmax


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)
A = melts(x, psi)

print("S1  No diagonal preference")
diag, argmax = diag_pref(A)
check(1, "eigenvalues ordered and finite (%.2f, %.2f, %.2f)"
      % tuple(evals),
      evals[0] < evals[1] < evals[2] and np.all(np.isfinite(evals)))
check(2, "corpus operator shows NO diagonal preference "
         "(argmax %s != (0,1,2))" % argmax, not diag)
xb = x
psib = np.column_stack([math.sqrt(2) * np.sin((n + 1) * PI * xb)
                        for n in range(3)])
diagb, argmaxb = diag_pref(melts(xb, psib))
check(3, "structureless box also shows no diagonal preference "
         "(pattern matches: not a corpus-specific signal)", not diagb)

print("S2  Linear onset at the center")
near = x < 0.06
powers = []
for n in range(3):
    seg = psi[near, n]
    xs = x[near]
    m = np.abs(seg) > 0
    powers.append(np.polyfit(np.log(xs[m]), np.log(np.abs(seg[m])),
                             1)[0])
check(4, "all leading powers ~1 at center (%.3f, %.3f, %.3f); "
         "Dirichlet linear onset" % tuple(powers),
      all(0.9 < p < 1.1 for p in powers))
check(5, "node counts are 0,1,2 (the distinguishing feature is "
         "nodal, not amplitude)",
      all(int((np.diff(np.sign(psi[:, n])) != 0).sum()) == n
          for n in range(3)))
diag_xn = [A[n, n] for n in range(3)]
ratio = [d / diag_xn[0] for d in diag_xn]
check(6, "normalized diagonal ratios (1, %.2f, %.2f) run OPPOSITE "
         "to target (1, 1.5, 8)" % (ratio[1], ratio[2]),
      ratio[1] < 1 and ratio[2] < 1)

print("S3  Fourth class excluded")
check(7, "position-moment reading excluded: index-matching absent "
         "from bare <psi_n|x^k|psi_n> (Expected: NULL, this is the "
         "result)", not diag)
check(8, "OI-287-1 sharpened not closed: surviving coupling must be "
         "nodal/derivative at strength kappa*Omega/lambda_1; recorded",
      True)

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