"""
n18 — The interacting arena. Can a record-creating (dephasing) scatterer
attract in a NON-GAUSSIAN medium?

n13 executed the coherence-shadow in the LINEAR vacuum: Maxwell's cancellation
survives quantization — dephasing grains repel (pump radiation absorbed by the
medium), never attract, and the repulsion dies in the transparent limit.
But Maxwell's theorem is a linear-response statement. In an interacting medium
a recorder's unavoidable backaction DRESSES it: the pump's local heating
delta<x^2> shifts the local stiffness through the interaction (Hartree:
mu_eff^2 = mu^2 + 3 g delta<x^2>), i.e. KEEPING RECORDS MANUFACTURES A
COHERENT DEFECT — and coherent stiffening defects attract (n13's control row;
Kenneth–Klich guarantees the sign for mirror-symmetric scatterers).

Hypothesis under test (the resurrection): in an interacting vacuum the
record-shadow force flips sign — dephasing grains attract for g > g*,
and g* -> 0 in the transparent limit, because the pump repulsion is
proportional to the medium absorption eta while the dressing is not.

PRE-REGISTERED GATES (written before the first run)
 G0  g = 0 reproduces n13 (repulsion at d=24, gamma=0.01, eta=0.005).
 G1  Sign of F_att vs g at fixed (d, gamma, eta): a sign flip at finite g*
     is the resurrection; monotone repulsion at all g is the second death.
 G2  Transparent limit: at fixed g above g*, the attraction must survive
     (grow relative to) eta -> 0.
 G3  Scaling: the induced part of the force scales as gamma_A * gamma_B
     (each grain's dressing is linear in its own record rate).
 G4  Mechanism audit: freeze the converged single-grain dressing profiles
     into STATIC (noiseless) defects; their coherent Casimir force must
     account for the g-induced part of the noisy force (same sign,
     same order of magnitude; ratio reported).
 G5  Range: the induced attraction's d-dependence follows the static
     dressed-defect force, not the pump repulsion.

Method: self-consistent Hartree on n13's exact testbed. Exact NESS by
Lyapunov equations at every Hartree step; "vacuum continuation" boundary
baths rebuilt from the current K each step; force from the ABBA-subtracted
momentum-flux (T^xx) discontinuity, evaluated with each configuration's own
converged stiffness profile. Interaction H_int = (g/4) sum_i x_i^4, treated
in leading Wick order. Normal ordering: the counterterm subtracts the CURRENT
K's own vacuum <x^2>, so the dressing is sourced purely by the noise-induced
excess dx2 = diag(sigma_noisy - sigma_undriven). That excess is a Lyapunov
integral of a positive-semidefinite source, hence dx2 >= 0 exactly: the
dressing only ever stiffens (mu_eff^2 >= mu^2), the interacting vacuum stays
uniform, and no tachyonic runaway is possible. (A first scheme referencing
the BARE uniform vacuum diverged through window-edge reflections feeding
ground-state reshaping — kept on the record, chronicle item 1.)

Caveats, registered: Hartree is the LEADING interaction effect, not the full
non-Gaussian theory; 1D; white-noise records; T=0. If the sign flip appears,
the claim is "the cancellation is not protected beyond linear response",
demonstrated at leading order — not a solved interacting QFT.

Runtime: ~10-20 min (numpy + scipy, no Monte Carlo).
"""

import numpy as np
from scipy.linalg import solve_continuous_lyapunov, eigh

# macOS Accelerate raises spurious divide/overflow FP flags on healthy
# matmuls (verified: fires with min w2 = 5.5e-4, no division anywhere).
# Real divergence is caught by the Hartree delta monitor below.
np.seterr(divide='ignore', over='ignore', invalid='ignore')

N        = 256
MU2      = 4e-4
LABS     = 28
ETA_MAX  = 0.5
ETA_BULK = 0.005
CENTER   = N // 2
SMEAR    = 2.0
OFF      = 8

ALPHA    = 0.7      # Hartree damping
TOL      = 1e-9     # convergence: max |delta mu^2|
MAXIT    = 60

# ---------------------------------------------------------------- testbed ---

def K_matrix(mu2_vec):
    K = np.diag(2.0 + mu2_vec).astype(float)
    K -= np.diag(np.ones(N - 1), 1) + np.diag(np.ones(N - 1), -1)
    return K

def weight(site):
    """Grain = fluctuating index coupled to the local strain (n13)."""
    idx = np.arange(N)
    g = np.exp(-0.5 * ((idx - site) / SMEAR) ** 2)
    g[np.abs(idx - site) > 4 * SMEAR] = 0.0
    c = np.zeros(N)
    c[:-1] = np.diff(g)
    return c / np.linalg.norm(c)

def ground_cov(K, temp=0.0):
    w2, U = eigh(K)
    w = np.sqrt(np.maximum(w2, 1e-12))
    occ = 0.5 * np.ones(N) if temp == 0 else 0.5 / np.tanh(0.5 * w / temp)
    sig = np.zeros((2 * N, 2 * N))
    sig[:N, :N] = U @ np.diag(occ / w) @ U.T
    sig[N:, N:] = U @ np.diag(occ * w) @ U.T
    return sig

def build_A_D(K, eta_bulk, temp=0.0):
    eta = np.full(N, eta_bulk)
    for j in range(LABS):
        s = ((LABS - j) / LABS) ** 2
        eta[j] += ETA_MAX * s
        eta[N - 1 - j] += ETA_MAX * s
    E = np.diag(np.concatenate([eta, eta]))
    A_H = np.zeros((2 * N, 2 * N))
    A_H[:N, N:] = np.eye(N)
    A_H[N:, :N] = -K
    A = A_H - E / 2
    sgs = ground_cov(K, temp)
    D = 0.5 * (E @ sgs + sgs @ E)
    return A, D

def ness(K, noisy=(), eta_bulk=ETA_BULK, temp=0.0):
    """noisy = ((site, gamma), ...); self-consistent record noise (n13).
    Returns (sig, svals, sig0) — sig0 is the grain-free NESS of the same K,
    i.e. the current vacuum, used as the normal-ordering reference."""
    A, D = build_A_D(K, eta_bulk, temp)
    sig0 = solve_continuous_lyapunov(A, -D)
    if not noisy:
        return sig0, np.zeros(0), sig0
    gammas = np.array([g for _, g in noisy])
    uvecs, parts = [], []
    for s, _ in noisy:
        w = weight(s)
        u = np.concatenate([w, np.zeros(N)])
        v = np.concatenate([np.zeros(N), w])
        uvecs.append(u)
        parts.append(solve_continuous_lyapunov(A, -np.outer(v, v)))
    n = len(noisy)
    b = np.array([u @ sig0 @ u for u in uvecs])
    M = np.array([[gammas[j] * (uvecs[i] @ parts[j] @ uvecs[i])
                   for j in range(n)] for i in range(n)])
    ev = np.linalg.eigvals(M)
    if np.max(ev.real) >= 1.0:
        raise RuntimeError(f"parametric instability (max eig M = {np.max(ev.real):.3f})")
    svals = np.linalg.solve(np.eye(n) - M, b)
    sig = sig0 + sum(gammas[j] * svals[j] * parts[j] for j in range(n))
    return sig, svals, sig0

K_UNI  = K_matrix(np.full(N, MU2))
XX_VAC = np.diag(ground_cov(K_UNI)[:N, :N]).copy()   # counterterm reference

# ---------------------------------------------------------------- Hartree ---

def hartree(noisy, g, eta_bulk=ETA_BULK, mu2_init=None):
    """Self-consistent NESS: mu_i^2 = MU2 + 3 g dx2_i, with dx2 the
    noise-induced <x^2> excess over the current vacuum (dx2 >= 0)."""
    mu2 = np.full(N, MU2) if mu2_init is None else mu2_init.copy()
    it_done = 0
    for it in range(MAXIT):
        K = K_matrix(mu2)
        sig, s, sig0 = ness(K, noisy, eta_bulk)
        dx2 = np.diag(sig[:N, :N]) - np.diag(sig0[:N, :N])
        tgt = MU2 + 3.0 * g * dx2
        new = (1 - ALPHA) * mu2 + ALPHA * tgt
        delta = np.max(np.abs(new - mu2))
        mu2 = new
        it_done = it + 1
        if delta < TOL:
            break
    else:
        print(f"    [warn] Hartree not converged to {TOL} in {MAXIT} iters "
              f"(last delta = {delta:.2e})")
    K = K_matrix(mu2)
    sig, s, _ = ness(K, noisy, eta_bulk)
    return sig, s, mu2, it_done

def stress_profile(sig, mu2_vec):
    sxx = sig[:N, :N]; spp = sig[N:, N:]
    pp  = 0.5 * (np.diag(spp)[:-1] + np.diag(spp)[1:])
    dx2 = np.diag(sxx)[:-1] + np.diag(sxx)[1:] - 2 * np.diag(sxx, 1)
    mxx = 0.5 * (mu2_vec[:-1] * np.diag(sxx)[:-1] + mu2_vec[1:] * np.diag(sxx)[1:])
    return 0.5 * pp + 0.5 * dx2 - 0.5 * mxx

def force_att(T_by, a, b, off=OFF):
    dT = T_by["AB"] - T_by["A0"] - T_by["0B"] + T_by["00"]
    FA = dT[a - off] - dT[a + off - 1]
    FB = dT[b + off - 1] - dT[b - off]
    return 0.5 * (FA + FB), FA, FB

def run_pair(a, b, gA, gB, g, eta_bulk=ETA_BULK, warm=None):
    """Full ABBA with per-configuration Hartree. Returns force + diagnostics."""
    T_by, info = {}, {}
    for tag in ("00", "A0", "0B", "AB"):
        noisy = []
        if "A" in tag: noisy.append((a, gA))
        if "B" in tag: noisy.append((b, gB))
        init = warm.get(tag) if warm is not None else None
        sig, s, mu2, its = hartree(tuple(noisy), g, eta_bulk, init)
        if warm is not None:
            warm[tag] = mu2
        T_by[tag] = stress_profile(sig, mu2)
        info[tag] = dict(mu2=mu2, s=s, iters=its)
    F, FA, FB = force_att(T_by, a, b)
    return F, FA, FB, info

# ------------------------------------------------------- static references --

def static_force(profiles, eta_bulk=ETA_BULK):
    """Casimir force between STATIC stiffness profiles, same ABBA."""
    T_by = {}
    for tag, mu2 in profiles.items():
        sig, _, _ = ness(K_matrix(mu2), (), eta_bulk)
        T_by[tag] = stress_profile(sig, mu2)
    return force_att(T_by, A_SITE, B_SITE)

def coherent_control(a, b, strength, eta_bulk=ETA_BULK):
    """n13's control: the same grain operator, static. Calibrates the sign."""
    T_by = {}
    for tag, sites in (("00", ()), ("A0", (a,)), ("0B", (b,)), ("AB", (a, b))):
        K = K_UNI.copy()
        for s_ in sites:
            w = weight(s_)
            K += strength * np.outer(w, w)
        sig, _, _ = ness(K, (), eta_bulk)
        T_by[tag] = stress_profile(sig, np.full(N, MU2))
    return force_att(T_by, a, b)

# ------------------------------------------------------------------- main ---

if __name__ == "__main__":
    D0 = 24
    A_SITE, B_SITE = CENTER - D0 // 2, CENTER - D0 // 2 + D0
    G0GAMMA = 0.01

    print("=== sign calibration: static coherent grains (d=24) ===")
    Fc, _, _ = coherent_control(A_SITE, B_SITE, 0.05)
    print(f"  F_att(static, strength 0.05) = {Fc:+.4e}   "
          f"(n13: coherent grains ATTRACT -> this sign = attraction)\n")

    print("=== G0: g = 0 must reproduce n13 (d=24, gamma=0.01, eta=0.005) ===")
    F0, FA0, FB0, info0 = run_pair(A_SITE, B_SITE, G0GAMMA, G0GAMMA, 0.0)
    sAB = info0["AB"]["s"]
    print(f"  F_att = {F0:+.4e}   (n13 reference: -1.3152e-05)")
    print(f"  pump <x~^2> per grain: {sAB}\n")

    # pilot diagnostic: how big is the heating the dressing feeds on?
    KU = K_matrix(np.full(N, MU2))
    sigA, sA, sig0A = ness(KU, ((A_SITE, G0GAMMA),))
    dx2A = np.diag(sigA[:N, :N]) - np.diag(sig0A[:N, :N])
    print(f"  single-grain noise excess: max dx2 = {dx2A.max():.3e}, "
          f"channel level ~ {dx2A[CENTER + 40]:.3e}\n")

    print("=== G1: force vs interaction g (d=24, gamma=0.01, eta=0.005) ===")
    print(f"{'g':>8} {'F_att':>12} {'max d(mu^2)':>12} {'iters(AB)':>10}")
    warm = {}
    G_SCAN = (0.0, 0.01, 0.03, 0.1, 0.3, 1.0)
    F_of_g = {}
    for g in G_SCAN:
        try:
            F, FA, FB, info = run_pair(A_SITE, B_SITE, G0GAMMA, G0GAMMA, g, warm=warm)
            dmu = np.max(np.abs(info["AB"]["mu2"] - MU2))
            F_of_g[g] = F
            print(f"{g:>8} {F:>+12.4e} {dmu:>12.3e} {info['AB']['iters']:>10}")
        except RuntimeError as e:
            print(f"{g:>8}  FAILED: {e}")
    print()

    # pick the smallest scanned g with attraction (if any) for the later gates
    g_res = next((g for g in G_SCAN if g > 0 and F_of_g.get(g, -1) > 0), None)
    if g_res is None:
        print(">>> G1 VERDICT: no sign flip in the scanned window — "
              "the record-shadow does not attract at leading interacting order.")
        g_res = G_SCAN[-1]   # still run the remaining gates for the record
    else:
        print(f">>> G1 VERDICT: sign flip — attraction from g = {g_res}.")
    print()

    print(f"=== G2: transparent limit (g = {g_res}, d=24, gamma=0.01) ===")
    print(f"{'eta_bulk':>9} {'F_att(g)':>12} {'F_att(0)':>12}")
    for eb in (0.02, 0.005, 0.00125):
        Fg, *_ = run_pair(A_SITE, B_SITE, G0GAMMA, G0GAMMA, g_res, eta_bulk=eb)
        F0e, *_ = run_pair(A_SITE, B_SITE, G0GAMMA, G0GAMMA, 0.0, eta_bulk=eb)
        print(f"{eb:>9} {Fg:>+12.4e} {F0e:>+12.4e}")
    print()

    print(f"=== G3: gamma scaling of the induced part (g = {g_res}, d=24) ===")
    print(f"{'gA':>7} {'gB':>7} {'dF = F(g)-F(0)':>15} {'dF/(gA*gB)':>12}")
    for gA, gB in ((0.005, 0.005), (0.01, 0.01), (0.02, 0.02), (0.02, 0.005)):
        Fg, *_ = run_pair(A_SITE, B_SITE, gA, gB, g_res)
        F0e, *_ = run_pair(A_SITE, B_SITE, gA, gB, 0.0)
        dF = Fg - F0e
        print(f"{gA:>7} {gB:>7} {dF:>+15.4e} {dF / (gA * gB):>12.4e}")
    print()

    print(f"=== G4: mechanism audit — static dressed defects (g = {g_res}) ===")
    _, _, mu2_A, _ = hartree(((A_SITE, G0GAMMA),), g_res)
    _, _, mu2_B, _ = hartree(((B_SITE, G0GAMMA),), g_res)
    exA, exB = mu2_A - MU2, mu2_B - MU2
    profiles = {"00": np.full(N, MU2), "A0": MU2 + exA,
                "0B": MU2 + exB, "AB": MU2 + exA + exB}
    Fs, _, _ = static_force(profiles)
    Fg, *_ = run_pair(A_SITE, B_SITE, G0GAMMA, G0GAMMA, g_res)
    F0e = F_of_g.get(0.0)
    dF = Fg - F0e
    print(f"  static Casimir of the frozen dressings: F = {Fs:+.4e}")
    print(f"  induced part of the noisy force  dF    = {dF:+.4e}")
    print(f"  ratio static/induced = {Fs / dF if dF else float('nan'):.3f}\n")

    print(f"=== G5: range (g = {g_res}, gamma=0.01, eta=0.005) ===")
    print(f"{'d':>4} {'F_att(g)':>12} {'F_att(0)':>12} {'dF':>12}")
    for d in (16, 24, 32, 48):
        aa, bb = CENTER - d // 2, CENTER - d // 2 + d
        Fg, *_ = run_pair(aa, bb, G0GAMMA, G0GAMMA, g_res)
        F0e, *_ = run_pair(aa, bb, G0GAMMA, G0GAMMA, 0.0)
        print(f"{d:>4} {Fg:>+12.4e} {F0e:>+12.4e} {Fg - F0e:>+12.4e}")