"""n51a_verify — numerical gates for the dichotomy proofs, as declared in
n51a_dichotomy_theorem.md (committed before this file)."""
import numpy as np

N = 200; MU2 = 1e-6

def K_of(dk):
    k = 1.0 + dk
    K = MU2*np.eye(N)
    for e in range(N):
        i, j = e, (e+1) % N
        K[i,i] += k[e]; K[j,j] += k[e]; K[i,j] -= k[e]; K[j,i] -= k[e]
    return K

def info_metric(dk):
    w2, V = np.linalg.eigh(K_of(dk))
    wv = np.sqrt(np.maximum(w2, 1e-30))
    U = np.zeros((N, N))
    for e in range(N): U[e] = V[e] - V[(e+1) % N]
    W1, W2 = np.meshgrid(wv, wv, indexing='ij')
    kern = 1.0/((W1+W2)**2 * W1*W2)
    B = np.einsum('en,em->enm', U, U).reshape(N, N*N)
    return 0.125*(B*kern.reshape(1, N*N)) @ B.T

G = info_metric(np.zeros(N))

# ---------------- T1: bipartite linear response is tidal -------------------
Au = np.zeros((N, N))
for i in range(N): Au[i,(i-1)%N] = 1; Au[i,i] += 1
GiAu = np.linalg.solve(G, Au.T)
Wu = Au @ GiAu

def dk_of(c):
    lam, *_ = np.linalg.lstsq(Wu, -c, rcond=None)
    return GiAu @ lam

def cell_avg(x): return 0.5*(x[0::2] + x[1::2])

def c2_shift(F):
    c = np.zeros(N); c[40] = F; c[139] = F
    dk = dk_of(c)
    return cell_avg(np.sqrt(1+dk) - 1), dk

F = 0.1
sp, dkp = c2_shift(F)
sm, _ = c2_shift(-F)
odd = 0.5*(sp - sm)        # linear-in-F part
even = 0.5*(sp + sm)       # quadratic part
cells = np.arange(N//2); cc = [20, 69]
win = np.array([5 <= min([min(abs(i-j), N//2-abs(i-j)) for j in cc]) <= 30 for i in cells])
ratio = np.abs(odd[win]).max()/np.abs(even[win]).max()
print(f"T1a: max|linear|/max|quadratic| at cells 5..30 = {ratio:.4f} -> {'PASS' if ratio < 0.05 else 'FAIL'}")
env = cell_avg(np.abs(dkp))
grad = np.abs(np.gradient(env))
def fitfr(y, x, m):
    Ad = np.vstack([np.ones(m.sum()), x[m]]).T
    cf, *_ = np.linalg.lstsq(Ad, y[m], rcond=None)
    r = y[m] - Ad@cf
    return np.sqrt(np.mean(r**2))/max(np.ptp(y[m]), 1e-300)
fr_env = fitfr(odd, env, win); fr_grad = fitfr(odd, grad, win)
print(f"T1b: linear part vs envelope fracRMS = {fr_env:.3%}, vs |grad envelope| = {fr_grad:.3%} -> "
      f"{'PASS (tidal)' if fr_grad < fr_env else 'FAIL'}")

# ---------------- T2: signed identity and Coulomb slope --------------------
As = np.zeros((N, N))
for i in range(N): As[i, i] = 1.0; As[i, (i-1) % N] = -1.0
GiAs = np.linalg.solve(G, As.T)
Ws = As @ GiAs
Wp = np.linalg.pinv(Ws, rcond=1e-10)
Kk = Wp[0]                          # kernel K(d) = W^+(0, d)
f = 0.3
def Eint(d):
    c = np.zeros(N); c[40] = f; c[(40+d) % N] = f; c -= c.mean()
    lam, *_ = np.linalg.lstsq(Ws, -c, rcond=None)
    dk = GiAs @ lam
    E_pair = 0.5*dk @ G @ dk
    c1 = np.zeros(N); c1[40] = f; c1 -= c1.mean()
    lam1, *_ = np.linalg.lstsq(Ws, -c1, rcond=None)
    dk1 = GiAs @ lam1
    return E_pair - 2*(0.5*dk1 @ G @ dk1)
ds = [11, 21, 41, 61, 81]
errs = []
for d in ds:
    lhs = Eint(d)
    rhs = f*f*(Kk[(40+d) % N - 40 + 0] if False else Wp[40, (40+d) % N])
    errs.append(abs(lhs - rhs)/max(abs(lhs), 1e-300))
print(f"T2a: |E_int(d) - f^2 K(d)|/|E_int| max = {max(errs):.2e} -> {'PASS' if max(errs) < 1e-8 else 'FAIL'}")
S = 32.4
dd = np.arange(5, 60)
kd = np.array([Wp[40, 40+d] for d in dd])
slope = np.polyfit(dd, Wp[40, 40] - kd - dd*(N-dd)/(2*S*N)*0, 1)[0]  # raw slope of K(0)-K(d)
# compare with d(N-d)/(2 S N)-form: fit K(0)-K(d) against d(N-d)/N
xf = dd*(N-dd)/N
cf = np.polyfit(xf, Wp[40, 40] - kd, 1)[0]
print(f"T2b: fitted 1/(2S) = {cf:.5f} vs expected {1/(2*S):.5f} -> "
      f"{'PASS' if abs(cf*2*S - 1) < 0.02 else 'FAIL'} (dev {abs(cf*2*S-1):.3%})")