#!/usr/bin/env python3
"""
N12 — THE WELD: one metric as the zero-reflection fixed point
==============================================================
THE IDEA. The welding condition (time-of-flight metric ~ resistance
metric) on a chain with springs c(x) and masses m(x) is m*c = const —
UNIFORM IMPEDANCE Z = sqrt(m c). Uniform impedance is the classic
reflectionless condition of wave physics. Reflections = records
radiated backward = the FEE (Folio VII). Hence the thesis:
    the single metric of GR is the no-fee fixed point;
    bimetric substrates bleed; geometry is selected by dissipation.

PRE-REGISTERED (before the runs):
  A. The resistance dictionary holds in ALL mass-substrates: MI-inferred
     distance tracks d_R (= static Green metric) even when masses make
     the causal metric maximally different. KILL: if in the anti-welded
     substrate (m = c) MI tracks the causal metric instead, Folio XV's
     verdict was confounded by uniform masses and falls.
  B. Wave packet through the lens: reflection finite for m uniform,
     ~zero for the welded substrate m = 1/c.
  C. Reflection vs family m = c^(-s): minimized at s = 1 (the weld).
"""
import numpy as np, json, warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

# ---------------- Test A: MI distance with dynamical masses ----------------
def covs_mass(c, msite, N, mreg=3e-4):
    K = np.zeros((N, N))
    for i in range(N - 1):
        K[i, i] += c[i]; K[i+1, i+1] += c[i]
        K[i, i+1] -= c[i]; K[i+1, i] -= c[i]
    K += np.eye(N) * mreg * mreg
    Mi = 1.0/np.sqrt(msite)
    Kt = (K * Mi).T * Mi                  # M^-1/2 K M^-1/2
    w2, U = np.linalg.eigh(Kt)
    w = np.sqrt(np.clip(w2, 1e-30, None))
    Xt = (U*(0.5/w)) @ U.T; Pt = (U*(0.5*w)) @ U.T
    X = (Xt * Mi).T * Mi                  # back to physical q,p
    P = (Pt * np.sqrt(msite)).T * np.sqrt(msite)
    return X, P

def S(X, P, R):
    ix = np.ix_(R, R)
    nu = np.sqrt(np.clip(np.linalg.eigvals(X[ix] @ P[ix]).real, 0.25, None))
    up, dn = nu+.5, nu-.5
    return float(np.sum(up*np.log(up)) -
                 np.sum(np.where(dn > 1e-12, dn*np.log(np.clip(dn, 1e-300, None)), 0)))

def MI(X, P, i, j):
    return S(X, P, [i]) + S(X, P, [j]) - S(X, P, [i, j])

N = 800; ctr = N//2
x = np.arange(N - 1) + 0.5
lens = 1.0 - (1.0 - 0.0625) * np.exp(-((x - ctr)/30)**2)
ones = np.ones(N)

X0, P0 = covs_mass(np.ones(N-1), ones, N)
bd = np.array([16, 24, 32, 48, 64, 96, 128, 192, 256, 320, 400], float)
bI = np.array([MI(X0, P0, int(ctr-d//2), int(ctr+(d+1)//2)) for d in bd])
def dinf(I): return float(np.exp(np.interp(np.log(I), np.log(bI)[::-1], np.log(bd)[::-1])))

def mass_profile(s):
    m = np.ones(N)
    m[1:] *= lens**(-s)   # node mass tracks adjacent spring; ends at c=1 anyway
    return m

print("TEST A — MI distance in three substrates (pre-registered: d_R wins in all)")
print(f"{'substrate':>22} {'d_eucl':>7} {'d_R':>7} {'d_T':>7} {'d_inf':>7} {'verdict':>14}")
resA = []
for label, s in (("m uniform (XV)", 0.0), ("WELDED m=1/c", 1.0), ("ANTI-WELD m=c", -1.0)):
    msite = mass_profile(s)
    X, P = covs_mass(lens, msite, N)
    for d in (96, 160):
        i, j = int(ctr-d//2), int(ctr+(d+1)//2)
        dR = float(np.sum(1.0/lens[i:j]))
        # local speed v = sqrt(c/m_local); flight = sum 1/v
        mloc = 0.5*(msite[i:j] + msite[i+1:j+1])
        dT = float(np.sum(np.sqrt(mloc/lens[i:j])))
        I = MI(X, P, i, j); di = dinf(I)
        lean = "RESISTANCE" if abs(di-dR) < abs(di-dT) else "CAUSAL"
        resA.append(dict(sub=label, d=d, dR=dR, dT=dT, dinf=di, lean=lean))
        print(f"{label:>22} {d:7d} {dR:7.1f} {dT:7.1f} {di:7.1f} {lean:>14}")

# ---------------- Test B/C: classical reflection off the lens ----------------
def reflection(s, Nw=2400, T=1500.0, dt=0.04):
    xw = np.arange(Nw - 1) + 0.5
    cw = 1.0 - (1.0 - 0.0625) * np.exp(-((xw - Nw/2)/30)**2)
    mw = np.ones(Nw); mw[1:] *= cw**(-s)
    q = np.zeros(Nw); p = np.zeros(Nw)
    x0, sig, k0 = Nw/2 - 500, 40.0, 0.35
    xs = np.arange(Nw)
    q = np.exp(-((xs-x0)/sig)**2) * np.cos(k0*(xs-x0))
    p = mw * np.exp(-((xs-x0)/sig)**2) * np.sin(k0*(xs-x0))  # right-mover-ish
    nst = int(T/dt)
    for _ in range(nst):
        f = np.zeros(Nw)
        dq = np.diff(q)
        f[:-1] += cw*dq; f[1:] -= cw*dq
        p += dt*f
        q += dt*p/mw
    eloc = 0.5*p**2/mw
    eloc[:-1] += 0.25*cw*np.diff(q)**2; eloc[1:] += 0.25*cw*np.diff(q)**2
    left = float(np.sum(eloc[:Nw//2 - 60])); right = float(np.sum(eloc[Nw//2 + 60:]))
    return left/(left+right)

print("\nTEST B/C — backscatter R(s) for m = c^(-s) (pre-registered: min at s=1)")
resC = []
for s in (-1.0, 0.0, 0.5, 1.0, 1.5, 2.0):
    R = reflection(s)
    resC.append(dict(s=s, R=R))
    tag = "  <- THE WELD" if s == 1.0 else ""
    print(f"  s = {s:+4.1f}   R = {R:.5f}{tag}")
json.dump(dict(A=resA, C=resC), open('/Users/antoine/agi/ledger/n12_results.json','w'), indent=1)
print("-> n12_results.json")