"""n31 — relational clock records: exact Level-I judge.
Chain tested (statuses in output):
 1 [E,L_ij]=0 -> all moments of E conserved under records; local eps diffuse (kernel L_G)
 2 Fisher I = 4 gamma L_G ; exact zero mode ; lambda(q) ~ q^2 gate p in [1.95,2.05]
 3 pure records: NO static stiffness for tau (drift matrix has zero tau-restoring block)
 4 relational KMS partner L_e = sqrt(g w)[(B tau)_e + i eta (B eps)_e]:
   still [E,L]=0 ; derived drift -2*g*eta*L_G on tau AND eps
 5 static field: 2 g eta L_G <tau> = lam (s - sbar)  ->  theta = (lam/2geta) L_G^+ (s-sbar)
 CONTROLS: A absolute records -> [E,L] != 0 (heating) ; D pinning -> Yukawa screening.
All Gaussian: first moments exact linear ODEs; covariances exact Lyapunov flow."""
import numpy as np

N=64; GAM=0.3; ETA=0.5; LAM=0.2
B=np.zeros((N,N))                       # ring incidence
for e in range(N): B[e,e]=1.0; B[e,(e+1)%N]=-1.0
LG=B.T@B

print("== 1: E exactly QND under relational records; eps diffuses ==")
# covariance flow under records L_e = sqrt(g)(B tau)_e : dSig_ee'/dt adds g * LG in eps-block
# (derived: [L, r] c-number => diffusion (Omega u)(Omega u)^T, u in tau => lands in eps)
D_eps = GAM*LG
varE_rate = float(np.ones(N)@D_eps@np.ones(N))
print(f"  d Var(E)/dt = 1^T (gamma L_G) 1 = {varE_rate:.2e}   (exact 0: zero mode)")
print(f"  d<eps_k^2>/dt = gamma*deg = {GAM*LG[3,3]:.3f} per site (energy redistributes)")
uabs=np.zeros(N); uabs[0]=1.0
print(f"  CONTROL A (absolute record tau_0): d Var(E)/dt = {float(np.ones(N)@(GAM*np.outer(uabs,uabs))@np.ones(N)):.3f}  != 0 -> heats")

print("== 2: Fisher kernel ==")
I_an = 4*GAM*LG
ev=np.linalg.eigvalsh(I_an)
print(f"  zero mode |lambda_0| = {abs(ev[0]):.2e}   (gate < 1e-12)")
q=2*np.pi*np.arange(1,7)/N
lam_q=4*GAM*2*(1-np.cos(q))
p=np.polyfit(np.log(q),np.log(lam_q),1)[0]
print(f"  small-q fit p = {p:.4f}   (gate [1.95,2.05])")
# Monte-Carlo score check on N=8 ring
rng=np.random.default_rng(1); n8=8
B8=np.zeros((n8,n8))
for e in range(n8): B8[e,e]=1; B8[e,(e+1)%n8]=-1
tau0=rng.normal(size=n8); T=2000; dt=0.01
sc=np.zeros((4000,n8))
for m in range(4000):
    dW=rng.normal(size=(T,n8))*np.sqrt(dt)
    dY=2*np.sqrt(GAM)*(B8@tau0)*dt+dW
    # score for tau: grad log P = sum 2 sqrt(g) B^T (dY - 2 sqrt(g) B tau dt)/1
    sc[m]=2*np.sqrt(GAM)*B8.T@ (dY-2*np.sqrt(GAM)*(B8@tau0)*dt).sum(0)
I_mc=(sc.T@sc)/4000/(T*dt)
I_th=4*GAM*(B8.T@B8)
off=np.abs(I_mc-I_th).max()/np.abs(I_th).max()
print(f"  MC score covariance vs 4 gamma L_G: max rel dev = {off:.3f}  (MC ~ 1/sqrt(4000))")

print("== 3: pure records -> tau restoring drift is EXACTLY zero ==")
print("  drift(tau-sector) from records = 0 (derived: [tau_j, L]=0 for Hermitian tau-records)")
print("  => static susceptibility undefined; K_phys does NOT emerge. ARROW BREAKS HERE.")

print("== 4-5: relational KMS partner -> derived drift -2 g eta L_G ; static field ==")
# first-moment stationary: 2*g*eta*LG <tau> = lam*(s - sbar)   [uniform drift on zero mode]
s=np.zeros(N); s[N//4]=1.0; s[N//4+1]=1.0            # matter lump (energy density)
src=LAM*(s-s.mean())
Z=2*GAM*ETA
LGp=np.linalg.pinv(LG)
theta=LGp@src/Z
# verify by integrating the actual ODE dtau/dt = lam(s-sbar) - Z LG tau
tau=np.zeros(N)
for _ in range(400000):
    tau+= 1e-3*(src - Z*LG@tau)
tau-=tau.mean()
dev=np.abs(tau-theta).max()/np.abs(theta).max()
print(f"  ODE NESS vs (1/Z) L_G^+ src : max rel dev = {dev:.2e}   Z = 2*gamma*eta = {Z:.3f} (derived)")
# 1D Green: field of lump on ring = parabolic/tent profile; kernel check: LG theta = src/Z
kres=np.abs(LG@theta - src/Z).max()
print(f"  kernel identity |L_G theta - src/Z|_max = {kres:.2e}   -> K(q) = Z q^2, K(0)=0 exact")

print("== CONTROL D: pinning mu^2 tau^2 -> screening ==")
MU2=0.05
thetaD=np.linalg.solve(Z*LG+MU2*np.eye(N), src)
# decay length vs Yukawa prediction sqrt(Z/mu^2)
prof=np.abs(thetaD-thetaD.mean()); i0=N//4
tail=prof[i0+3:i0+13]
xi=-1/np.polyfit(np.arange(10),np.log(tail),1)[0]
print(f"  measured decay length = {xi:.2f}   Yukawa prediction sqrt(Z/mu^2) = {np.sqrt(Z/MU2):.2f}")
print("done.")