"""n47b — the linear clock channel, per n47b_preregistration.md.
Two-channel decomposition across an F scan; parity-clock probe (L4)."""
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

w2, V = np.linalg.eigh(K_of(np.zeros(N)))
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)
G = 0.125*(B*kern.reshape(1, N*N)) @ B.T
A = np.zeros((N, N))
for i in range(N): A[i,(i-1)%N] = 1; A[i,i] += 1
Ginv_At = np.linalg.solve(G, A.T)
Wm = A @ Ginv_At

def sqrtm_sym(K):
    w2, Vv = np.linalg.eigh(K)
    return Vv @ np.diag(np.sqrt(np.maximum(w2, 1e-30))) @ Vv.T

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

s1, s2 = 40, 139
NC = N//2; cc = [s1//2, s2//2]
far = np.array([min([min(abs(i-j), NC-abs(i-j)) for j in cc]) > 2 for i in range(NC)])

def observables(F, swap=False):
    a, b = (s2, s1) if swap else (s1, s2)   # parity swap: consumers on odd sites
    c = np.zeros(N); c[a] = F; c[b] = F
    if swap: c = np.roll(c, 1)              # move both to opposite sublattice
    lam, *_ = np.linalg.lstsq(Wm, -c, rcond=None)
    dk = Ginv_At @ lam
    Om = sqrtm_sym(K_of(dk)); OmI = np.linalg.inv(Om)
    Om0 = sqrtm_sym(K_of(np.zeros(N))); OmI0 = np.linalg.inv(Om0)
    w1 = np.diag(Om)/np.diag(Om0) - 1        # C1 per site
    w2c = np.sqrt(2*(1+dk))/np.sqrt(2) - 1   # C2 per edge
    w3 = np.diag(OmI0)/np.diag(OmI) - 1      # C3 per site
    E = cell_avg(np.abs(dk)); Q = cell_avg(dk**2)
    sh = [cell_avg(w1), cell_avg(w2c), cell_avg(w3)]
    intra = w1[0::2] - w1[1::2]              # L4: intra-cell asymmetry of C1
    return sh, E, Q, intra

Fs = (0.075, 0.15, 0.3)
data = {F: observables(F) for F in Fs}

names = ["C1 site", "C2 bond", "C3 corr"]
print("n47b: joint two-channel fit s = a*E + b*Q over F scan, far cells")
res = {}
for i, nm in enumerate(names):
    X, y = [], []
    for F in Fs:
        sh, E, Q, _ = data[F]
        X.append(np.vstack([E[far], Q[far]]).T)
        y.append(sh[i][far])
    X = np.vstack(X); y = np.concatenate(y)
    cf, rss, *_ = np.linalg.lstsq(X, y, rcond=None)
    resid = y - X@cf
    dof = len(y) - 2
    cov = np.linalg.inv(X.T@X) * (resid@resid)/dof
    se = np.sqrt(np.diag(cov))
    fr = np.sqrt(np.mean(resid**2))/np.ptp(y)
    res[nm] = (cf, se, fr)
    print(f"  {nm}: a = {cf[0]:+.6f} (se {se[0]:.2e}, |a|/se = {abs(cf[0])/se[0]:.1f})  "
          f"b = {cf[1]:+.5f} (se {se[1]:.2e})  fit fracRMS = {fr:.3%}")

a1, a3 = res["C1 site"][0][0], res["C3 corr"][0][0]
se1, se3 = res["C1 site"][1][0], res["C3 corr"][1][0]
l1 = abs(a1) > 10*se1 and abs(a3) > 10*se3
print(f"L1 linear channel exists (both spectral clocks, |a|>10se): {'PASS' if l1 else 'FAIL — metric reading DIES'}")
if a3 != 0:
    sp = abs(a1/a3 - 1)
    print(f"L2 universality: a_C1/a_C3 - 1 = {sp:+.2%} -> "
          f"{'PASS' if sp < 0.10 else ('DEAD (clock-dependent)' if sp > 0.50 else 'PARTIAL')}")
a2, b2 = res["C2 bond"][0]
se2 = res["C2 bond"][1]
l3 = bool(abs(a2) < 10*se2[0]) and abs(b2 - (-0.125)) < 0.02*0.125
print(f"L3 control: a_C2 = {a2:+.2e} (se {se2[0]:.1e}), b_C2 = {b2:+.5f} vs -0.125 -> {'PASS' if l3 else 'FAIL'}")

print("L4 parity clocks: intra-cell asymmetry of C1")
for F in Fs:
    _, _, _, intra = data[F]
    print(f"  F={F}: max|delta| = {np.abs(intra[far]).max():.3e}, mean delta @far = {intra[far].mean():+.3e}")
i1 = data[0.15][3][far]; i2 = data[0.3][3][far]
mask = np.abs(i1) > 1e-12
expo = np.median(np.log(np.abs(i2[mask]/i1[mask]))/np.log(2))
_, _, _, intra_sw = observables(0.3, swap=True)
flip = np.corrcoef(data[0.3][3][far], intra_sw[far])[0,1]
print(f"  F-scaling exponent of delta = {expo:.3f} (1 => parity-sensitive linear; 2 => blind)")
print(f"  correlation under consumer parity swap = {flip:+.3f} (-1 => sign flip = anti-GR signature)")