#!/usr/bin/env python3
"""
attack11.py -- independent reproduction of the "Attack 11" live-neighbour
regeneration numbers (source note: 2026-09-07_ATTACK_11_LIVE_NEIGHBOR_REGENERATION.md,
model conventions from 2026-09-05_D40_COHERENT_DIRICHLET_CONTROL.md, Sec. 0).

Model (hbar = 1, J = 1 unless stated):
  spins sigma_i in {+1,-1}, computational (Z) basis, bit 0 <-> Z=+1 ("+").
  H0 = -J sum_<ij> Z_i Z_j - J sum_i e_i Z_i          (e_i = frozen exterior field)
  h_i = sum_{j~i, live} Z_j + e_i
  K  = sum_i g X_i 1_{h_i = 0}                         ([K, H0] = 0)
  c_i(sigma) = nu / (1 + exp(2 beta J sigma_i h_i(sigma)))   (heat bath)
  monitoring L_i = sqrt(kappa/2) Z_i
  L(rho) = -i[H0+K, rho] + sum_J D[J] rho,  D[J]rho = J rho J^+ - {J^+J, rho}/2

Two jump-resolution conventions for the thermal jumps:
  R1  fully resolved rank-one flips  J_{i,sigma} = sqrt(c_i(sigma)) |sigma^i><sigma|
  R2  neighbour-resolved (D40)       J_{i,s,eta} = sqrt(c_i(s,eta)) |-s><s|_i (x) Pi_eta(star) (x) 1

Geometry: 2x3 strip, order U0=(top,left) U1=(top,mid) U2=(bot,left) U3=(bot,mid)
D0=(top,right) D1=(bot,right).  Live bonds (U0,U1),(U2,U3),(U0,U2),(U1,U3),(U1,D0),
(U3,D1),(D0,D1).  Frozen-plus exterior fields e = (2,1,2,1,2,2).

Channel Phi_T^+ : B(U0,U2,D0,D1) (16-dim) -> B(D0,D1) (4-dim): inject U1=U3=+,
evolve the 6-spin Lindbladian for time T, trace U0..U3.

Coefficients (unnormalised Choi J(Phi) = sum_ij |i><j| (x) Phi(|i><j|)):
  p_tau   = lambda_min[(I (x) tau)^(-1/2) J (I (x) tau)^(-1/2)]           (full-rank tau)
  r_clean = max{ r : J - r I (x) |++><++| >= 0 }                          (rank-one target)
          = lambda_min of the Schur complement J/J_{S'S'} where S = range(I (x) |++><++|)
          (cross-checked by 1/lambda_max(P J^{-1} P) and by bisection on positivity).

Run:  python3 attack11.py            (~2 min unloaded; --dense adds a dense-Pade-expm cross-check, 3-20 min per convention depending on load)
No files are written.
"""
import os
for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"):
    os.environ.setdefault(_v, "1")          # small dense blocks: extra BLAS threads only add contention
import argparse
import time

import numpy as np
import scipy.linalg as sla
import scipy.sparse as sp
from scipy.sparse.linalg import expm_multiply

# ---------------------------------------------------------------- geometry
NAMES6 = ["U0", "U1", "U2", "U3", "D0", "D1"]
BONDS6 = [(0, 1), (2, 3), (0, 2), (1, 3), (1, 4), (3, 5), (4, 5)]
E_PLUS6 = (2, 1, 2, 1, 2, 2)
IN_SITES6 = [0, 2, 4, 5]          # U0, U2, D0, D1
FIXED6 = {1: 0, 3: 0}             # U1 = U3 = + (bit 0)
OUT_SITES6 = [4, 5]               # D0, D1

# frozen two-spin comparison cell: D0, D1 with U1 = U3 = + frozen
BONDS2 = [(0, 1)]
E_FROZEN2 = (3, 3)                # 2 exterior + 1 frozen interface spin each

# package values (from the Attack 11 note / task statement)
PKG = {
    "p_live": {2: .6844846359, 4: .9284100882, 6: .9827314629, 8: .9947830953},
    "p_frozen": {4: .9542899597, 8: .9988628010},
    "r_clean": {4: .9127059959, 8: .9779563010},
    "tau_pp": .983084962,
    "stress": [("all plus exterior", (2, 1, 2, 1, 2, 2), .9284101),
               ("one far upstream exterior defect", (0, 1, 2, 1, 2, 2), .9219624),
               ("one interface-adjacent exterior defect", (2, -1, 2, 1, 2, 2), .8959041),
               ("both interface external neighbours flipped", (2, -1, 2, -1, 2, 2), .8495878),
               ("far+interface defects on both rows", (0, -1, 0, -1, 2, 2), .8143947)],
    "spot": [((.4, .3, .8), .8912585886, .8178528933),
             ((.6, .3, .8), .9284100882, .9127059959),
             ((1.0, 1.0, .8), .9516407118, .9509968209)],
    "iface_flip": {4: .03812, 8: .02145},
    "bracket": {4: (.02127, .08507), 8: (.001988, .007951)},
}


# ---------------------------------------------------------------- model
def configs(n):
    d = 2 ** n
    x = np.arange(d)
    bits = np.array([(x >> (n - 1 - k)) & 1 for k in range(n)]).T   # (d, n)
    return bits, 1 - 2 * bits                                        # bits, z


def build_model(n, bonds, e, beta, Jc=1.0, g=0.3, nu=1.0, kappa=0.8, conv="R1",
                rate_sign=+1.0, mon_rate=None):
    """Return dict with sparse Liouvillian (row-major vec), H0, K, pi."""
    d = 2 ** n
    bits, z = configs(n)
    nbr = [[] for _ in range(n)]
    for a, b in bonds:
        nbr[a].append(b)
        nbr[b].append(a)
    e = np.asarray(e, float)
    h = np.stack([z[:, nbr[i]].sum(axis=1) + e[i] for i in range(n)], axis=1)   # (d, n)
    E = -Jc * sum(z[:, a] * z[:, b] for a, b in bonds) - Jc * (z * e).sum(axis=1)
    c = nu / (1.0 + np.exp(rate_sign * 2.0 * beta * Jc * z * h))                  # (d, n)
    flip = [1 << (n - 1 - i) for i in range(n)]
    # tie Hamiltonian
    rows, cols = [], []
    for i in range(n):
        xs = np.nonzero(h[:, i] == 0)[0]
        rows.append(xs ^ flip[i])
        cols.append(xs)
    rows = np.concatenate(rows) if rows else np.zeros(0, int)
    cols = np.concatenate(cols) if cols else np.zeros(0, int)
    K = sp.csr_matrix((np.full(len(rows), g, float), (rows, cols)), shape=(d, d))
    H0 = sp.diags(E)
    H = (H0 + K).tocsr()
    I = sp.identity(d, format="csr")
    Lh = -1j * (sp.kron(H, I) - sp.kron(I, H.T))
    km = kappa if mon_rate is None else mon_rate
    G = c.sum(axis=1) + n * km / 2.0                       # sum_J J^+J (diagonal)
    diag_anti = -0.5 * (G[:, None] + G[None, :]).ravel()
    diag_mon = (km / 2.0) * (z @ z.T).astype(float).ravel()
    Ld = sp.diags(diag_anti + diag_mon)
    # thermal sandwich terms  sum_J J (x) J^*
    rr, cc, vv = [], [], []
    x = np.arange(d)
    for i in range(n):
        if conv == "R1":
            xi = x ^ flip[i]
            rr.append(xi * d + xi)
            cc.append(x * d + x)
            vv.append(c[:, i])
        elif conv == "R2":
            mask = flip[i]
            for j in nbr[i]:
                mask |= flip[j]
            X, Y = np.meshgrid(x, x, indexing="ij")
            sel = ((X ^ Y) & mask) == 0
            Xs, Ys = X[sel], Y[sel]
            rr.append((Xs ^ flip[i]) * d + (Ys ^ flip[i]))
            cc.append(Xs * d + Ys)
            vv.append(c[Xs, i])
        else:
            raise ValueError(conv)
    S = sp.csr_matrix((np.concatenate(vv), (np.concatenate(rr), np.concatenate(cc))),
                      shape=(d * d, d * d))
    L = (Lh + Ld + S).tocsr()
    w = np.exp(-beta * E)
    pi = w / w.sum()
    return dict(L=L, H0=H0, K=K, pi=pi, E=E, n=n, d=d, bits=bits, h=h, c=c)


def model_checks(m):
    d = m["d"]
    L = m["L"]
    comm = (m["K"] @ m["H0"] - m["H0"] @ m["K"])
    comm = abs(comm).max() if comm.nnz else 0.0
    tvec = np.zeros(d * d)
    tvec[np.arange(d) * d + np.arange(d)] = 1.0
    tp = np.abs(L.T @ tvec).max()                  # trace preservation
    stat = np.abs(L @ np.diag(m["pi"]).ravel()).max()
    return comm, tp, stat


# ---------------------------------------------------------------- channel
def channel_setup(n, in_sites, fixed, out_sites):
    d = 2 ** n
    din, dout = 2 ** len(in_sites), 2 ** len(out_sites)
    rest = [k for k in range(n) if k not in out_sites]
    xin = np.zeros(din, int)
    for i in range(din):
        x = 0
        for pos, s in enumerate(in_sites):
            b = (i >> (len(in_sites) - 1 - pos)) & 1
            x |= b << (n - 1 - s)
        for s, b in fixed.items():
            x |= b << (n - 1 - s)
        xin[i] = x
    # input matrix units |i><j|, column k = i*din + j
    B = np.zeros((d * d, din * din), complex)
    for i in range(din):
        for j in range(din):
            B[xin[i] * d + xin[j], i * din + j] = 1.0

    def sub(x, sites):
        v = 0
        for s in sites:
            v = (v << 1) | ((x >> (n - 1 - s)) & 1)
        return v
    outi = np.array([sub(x, out_sites) for x in range(d)])
    resti = np.array([sub(x, rest) for x in range(d)])
    XA, XB = np.meshgrid(np.arange(d), np.arange(d), indexing="ij")
    sel = resti[XA] == resti[XB]
    XA, XB = XA[sel], XB[sel]
    PT = sp.csr_matrix((np.ones(len(XA)), (outi[XA] * dout + outi[XB], XA * d + XB)),
                       shape=(dout * dout, d * d))
    return dict(B=B, PT=PT, din=din, dout=dout, xin=xin, d=d)


def evolve(L, B, times):
    """Return {T: exp(T L) B} using the semigroup property between sorted times."""
    out, Y, t0 = {}, B, 0.0
    for T in sorted(times):
        if T > t0:
            Y = expm_multiply(L * (T - t0), Y)
        out[T] = Y
        t0 = T
    return out


def heisenberg_channel(L, ch, times, extra_obs=None):
    """Channel matrices via the Heisenberg picture (16 observables instead of 256 inputs).

    Phi(rho)_ab = Tr[O_ab^+ e^{TL}(rho)] with O_ab = 1_rest (x) |a><b|_out, whose row-major
    vectorisations are the rows of PT.  So PT e^{TL} B = (e^{T L^H} PT^H)^H B.
    Returns {T: (Phi matrix (dout^2 x din^2), extra expectation rows)}."""
    O = ch["PT"].T.toarray().astype(complex)             # (d^2, dout^2)
    if extra_obs is not None:
        O = np.hstack([O, extra_obs])
    Wt = evolve(L.conj().T.tocsr(), O, times)
    out = {}
    k = ch["PT"].shape[0]
    for T, W in Wt.items():
        R = W.conj().T @ ch["B"]                           # (dout^2 [+extra], din^2)
        out[T] = (R[:k], R[k:])
    return out


def choi_from(PT, Y, din, dout):
    """PT=None means Y already is the channel matrix (dout^2 x din^2)."""
    Phi = np.asarray(Y if PT is None else PT @ Y)          # (dout^2, din^2)
    Phi = Phi.reshape(dout, dout, din, din)                # [a,b,i,j]
    return Phi.transpose(2, 0, 3, 1).reshape(din * dout, din * dout)


def p_tau(Jm, tau, din):
    w, V = np.linalg.eigh(tau)
    tmh = V @ np.diag(w ** -0.5) @ V.conj().T
    M = np.kron(np.eye(din), tmh)
    A = M @ Jm @ M
    return np.linalg.eigvalsh((A + A.conj().T) / 2).min()


def r_rank_one(Jm, din, dout, k=0):
    """max r with J - r I(x)|k><k| >= 0. Returns (schur, inverse, bisection)."""
    Jh = (Jm + Jm.conj().T) / 2
    S = np.array([i * dout + k for i in range(din)])
    Sc = np.array([q for q in range(din * dout) if q not in set(S)])
    A, Bm, C = Jh[np.ix_(S, S)], Jh[np.ix_(S, Sc)], Jh[np.ix_(Sc, Sc)]
    schur = A - Bm @ np.linalg.solve(C, Bm.conj().T)
    r_schur = np.linalg.eigvalsh((schur + schur.conj().T) / 2).min()
    Ji = np.linalg.inv(Jh)
    Pi = Ji[np.ix_(S, S)]
    r_inv = 1.0 / np.linalg.eigvalsh((Pi + Pi.conj().T) / 2).max()
    P = np.zeros(din * dout)
    P[S] = 1.0
    lo, hi = 0.0, 1.0
    for _ in range(60):
        mid = (lo + hi) / 2
        if np.linalg.eigvalsh(Jh - mid * np.diag(P)).min() >= 0:
            lo = mid
        else:
            hi = mid
    return r_schur, r_inv, lo, np.linalg.eigvalsh(C).min()


def tau_D_plus(beta, Jc=1.0):
    # Gibbs of D0,D1 with frozen plus interface + plus exterior: exp(bJ(Z0Z1 + 3Z0 + 3Z1))
    _, z = configs(2)
    w = np.exp(beta * Jc * (z[:, 0] * z[:, 1] + 3 * z[:, 0] + 3 * z[:, 1]))
    return np.diag(w / w.sum())


# ---------------------------------------------------------------- helpers
def agree(a, b):
    if a == b:
        return 16.0
    return -np.log10(abs(a - b) / abs(b))


def fmt_row(label, val, pkg, width=46):
    if pkg is None:
        return f"  {label:<{width}} {val:.10f}"
    return f"  {label:<{width}} {val:.10f}   pkg {pkg:.10f}   diff {val - pkg:+.2e}   ~{agree(val, pkg):4.1f} digits"


# ---------------------------------------------------------------- main
def live_channel(conv, times, e=E_PLUS6, beta=.6, g=.3, nu=1.0, kappa=.8, ch=None, **kw):
    """Choi matrices {T: J(Phi_T^+)} (64x64) plus interface-flip rows, Heisenberg route."""
    m = build_model(6, BONDS6, e, beta, g=g, nu=nu, kappa=kappa, conv=conv, **kw)
    d = 64
    bits = configs(6)[0]
    not_pp = ((bits[:, 1] != 0) | (bits[:, 3] != 0)).astype(complex)
    obs = np.zeros((d * d, 1), complex)
    obs[np.arange(d) * d + np.arange(d), 0] = not_pp          # 1 - projector(U1=U3=+)
    out = heisenberg_channel(m["L"], ch, times, extra_obs=obs)
    chois = {T: choi_from(None, out[T][0], 16, 4) for T in times}
    flips = {T: np.real(out[T][1][0]) for T in times}           # indexed by column i*16+j
    return m, chois, flips


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--dense", action="store_true",
                    help="also cross-check with a dense scipy.linalg.expm of the 4096x4096 "
                         "Liouvillian (about 3 min per convention)")
    args = ap.parse_args()
    np.seterr(all="ignore")   # numpy 2.0 + Apple Accelerate emits spurious matmul FP warnings;
    import warnings           # every result below is explicitly checked to be finite instead.
    warnings.filterwarnings("ignore", category=RuntimeWarning)
    t_start = time.time()
    beta, g, nu, kappa = .6, .3, 1.0, .8
    TS = (2, 4, 6, 8)
    tau = tau_D_plus(beta)
    ch6 = channel_setup(6, IN_SITES6, FIXED6, OUT_SITES6)
    ch2 = channel_setup(2, [0, 1], {}, [0, 1])
    W = 118
    print("=" * W)
    print("ATTACK 11 REPRODUCTION -- 2x3 strip, live shared boundary.  Baseline beta J=.6, g=.3, nu=1, kappa=.8, J=1")
    print("=" * W)
    print(fmt_row("<++|tau_D^+|++>", tau[0, 0], PKG["tau_pp"]))

    # ---------------- frozen two-spin cell (R1 == R2 here: the star is the whole cell)
    mf = build_model(2, BONDS2, E_FROZEN2, beta, g=g, nu=nu, kappa=kappa, conv="R1")
    mf2 = build_model(2, BONDS2, E_FROZEN2, beta, g=g, nu=nu, kappa=kappa, conv="R2")
    same_frozen = abs(mf["L"] - mf2["L"]).max()
    outf = heisenberg_channel(mf["L"], ch2, [4, 8])
    pf = {T: p_tau(choi_from(None, outf[T][0], 4, 4), tau, 4) for T in (4, 8)}

    res = {}
    for conv in ("R1", "R2"):
        t0 = time.time()
        m, chois, flips = live_channel(conv, TS, ch=ch6)
        r = dict(m=m, chois=chois, flips=flips, checks=model_checks(m))
        r["p"] = {T: p_tau(chois[T], tau, 16) for T in TS}
        r["r"] = {T: r_rank_one(chois[T], 16, 4, 0) for T in TS}
        r["herm"] = max(np.abs(Jm - Jm.conj().T).max() for Jm in chois.values())
        r["tp"] = max(np.abs(np.einsum("iaja->ij", Jm.reshape(16, 4, 16, 4)) - np.eye(16)).max()
                      for Jm in chois.values())
        r["lmin"] = min(np.linalg.eigvalsh((Jm + Jm.conj().T) / 2).min() for Jm in chois.values())
        r["finite"] = all(np.isfinite(Jm).all() for Jm in chois.values())
        r["time"] = time.time() - t0
        res[conv] = r

    print("\n[1] Model / channel sanity checks")
    print(f"  frozen 2-spin cell: K nnz = {mf['K'].nnz} (h_D = Z_other + 3 is never 0), "
          f"max|L_R1 - L_R2| = {same_frozen:.1e}")
    for conv in ("R1", "R2"):
        comm, tp, stat = res[conv]["checks"]
        r = res[conv]
        print(f"  {conv}: nnz(L)={r['m']['L'].nnz:6d}  K nnz={r['m']['K'].nnz:3d}  |[K,H0]|={comm:.0e}  "
              f"|Tr L(.)|={tp:.0e}  |L(pi)|={stat:.0e}  Choi: herm {r['herm']:.0e}, TP {r['tp']:.0e}, "
              f"min eig {r['lmin']:.2e}, finite {r['finite']}  ({r['time']:.1f}s)")

    print("\n[2] Headline coefficients vs package (digits = -log10 |rel. diff|)")
    for conv in ("R1", "R2"):
        print(f"  --- thermal jumps {conv} "
              + ("(fully resolved rank-one flips |sigma^i><sigma|)" if conv == "R1"
                 else "(D40 neighbour-resolved |-s><s|_i (x) Pi_eta (x) 1)") + " ---")
        for T in TS:
            print(fmt_row(f"p_live  T={T}", res[conv]["p"][T], PKG["p_live"][T]))
        for T in TS:
            print(fmt_row(f"r_clean T={T}", res[conv]["r"][T][0], PKG["r_clean"].get(T)))
    print("  --- literally frozen-plus two-spin cell (L identical under R1/R2) ---")
    for T in (4, 8):
        print(fmt_row(f"p_frozen T={T}", pf[T], PKG["p_frozen"][T]))
    for conv in ("R1", "R2"):
        print(f"  live-neighbour penalty p_frozen - p_live [{conv}]:  "
              + "   ".join(f"T={T}: {pf[T] - res[conv]['p'][T]:.5f}" for T in (4, 8))
              + "     (package: .02588, .00408)")
    for conv in ("R1", "R2"):
        print(f"  q_clean = 1 - r_clean [{conv}]:  T=4 {1 - res[conv]['r'][4][0]:.6f}   T=8 {1 - res[conv]['r'][8][0]:.6f}"
              f"     4*q_clean: {4 * (1 - res[conv]['r'][4][0]):.5f}, {4 * (1 - res[conv]['r'][8][0]):.5f}"
              "   (package .087294/.022044; .34918/.08817)")

    print("\n[3] Rank-one-target coefficient: three methods, and the product bound p_live * <++|tau|++>")
    for conv in ("R1", "R2"):
        for T in TS:
            rs, ri, rb, cmin = res[conv]["r"][T]
            prod = res[conv]["p"][T] * tau[0, 0]
            # r = tau_++ * p  <=>  the lambda_min eigenvector of (I x tau^-1/2) J (I x tau^-1/2)
            # lies in range(I x |++><++|); report its weight there.
            tmh = np.kron(np.eye(16), np.diag(np.real(np.diag(tau)) ** -0.5))
            Mm = tmh @ res[conv]["chois"][T] @ tmh
            wv, Vv = np.linalg.eigh((Mm + Mm.conj().T) / 2)
            wS = np.sum(np.abs(Vv[np.arange(16) * 4, 0]) ** 2)
            print(f"  {conv} T={T}: Schur {rs:.12f} | 1/lmax(P J^-1 P) {ri:.12f} | bisection {rb:.12f} | "
                  f"product {prod:.12f} (r-prod {rs - prod:+.0e}) | weight of min-eigvec on I(x)|++> {wS:.12f}")

    print("\n[4] Independent propagation route at T=4: Schroedinger picture, 256 input matrix units")
    for conv in ("R1", "R2"):
        t0 = time.time()
        Y = evolve(res[conv]["m"]["L"], ch6["B"], [4])[4]
        Jm = choi_from(ch6["PT"], Y, 16, 4)
        print(f"  {conv}: max|J_schr - J_heis| = {np.abs(Jm - res[conv]['chois'][4]).max():.1e}   "
              f"p_live(schr) - p_live(heis) = {p_tau(Jm, tau, 16) - res[conv]['p'][4]:+.1e}   ({time.time() - t0:.1f}s)")
    if args.dense:
        print("  dense Pade expm (scipy.linalg.expm of 2L, then E2^k):")
        for conv in ("R1", "R2"):
            t0 = time.time()
            E2 = sla.expm(2.0 * res[conv]["m"]["L"].toarray())
            Y = ch6["B"]
            devs = []
            for T in TS:
                Y = E2 @ Y
                Jm = choi_from(ch6["PT"], Y, 16, 4)
                devs.append((T, p_tau(Jm, tau, 16) - res[conv]["p"][T]))
            del E2
            print(f"  {conv}: " + "  ".join(f"T={T}: dp={dp:+.1e}" for T, dp in devs)
                  + f"   ({time.time() - t0:.0f}s)")
    else:
        print("  (dense Pade expm cross-check available with --dense)")

    print("\n[5] Classical deep-upstream diagnostics (u = basis state of U0,U2; U1=U3=+ at entrance)")
    for conv in ("R1", "R2"):
        for T in (4, 8):
            J4 = res[conv]["chois"][T].reshape(16, 4, 16, 4)
            Ju = [J4[u * 4:(u + 1) * 4, :, u * 4:(u + 1) * 4, :].reshape(16, 16) for u in range(4)]
            pu = [p_tau(Ju[u], tau, 4) for u in range(4)]
            nrm = max(np.abs(np.linalg.eigvalsh((Ju[a] - Ju[b] + (Ju[a] - Ju[b]).conj().T) / 2)).sum()
                      for a in range(4) for b in range(a + 1, 4))
            fl = res[conv]["flips"][T]
            pflip = np.array([[fl[(u * 4 + dd) * 16 + (u * 4 + dd)] for dd in range(4)] for u in range(4)])
            print(f"  {conv} T={T}: min_u p(Phi_u -> tau) = {min(pu):.10f}  vs complete p_live {res[conv]['p'][T]:.10f}"
                  f"  (diff {min(pu) - res[conv]['p'][T]:+.1e})")
            print(f"           max pairwise ||J(Phi_u - Phi_u')||_1 bracket [{nrm / 4:.6f}, {nrm:.6f}]"
                  f"   package [{PKG['bracket'][T][0]}, {PKG['bracket'][T][1]}]")
            print(f"           max_u P(interface != ++ at time T):  D init ++: {pflip[:, 0].max():.5f}   "
                  f"D init tau_D^+: {(pflip @ np.real(np.diag(tau))).max():.5f}   worst D basis: {pflip.max():.5f}"
                  f"   package {PKG['iface_flip'][T]}")

    print("\n[6] Exterior-defect stress tests, T=4 (target tau_D^+ unchanged, entrance U1=U3=+)")
    print("    Interpretation: one flipped frozen exterior neighbour lowers that spin's e_i by 2.  'far' = one")
    print("    exterior neighbour of U0 (left column); 'interface-adjacent' = the exterior neighbour of U1;")
    print("    'both interface' = exterior neighbours of U1 and U3; 'far+interface both rows' = one flip at each")
    print("    of U0, U1, U2, U3.  (Top/bottom mirror symmetry makes the U0-vs-U2 choice immaterial.)")
    stress = {}
    for label, e, pkgv in PKG["stress"]:
        row = []
        for conv in ("R1", "R2"):
            _, chois, _ = live_channel(conv, [4], e=e, ch=ch6)
            row.append(p_tau(chois[4], tau, 16))
        stress[label] = row
        print(f"  {label:<43} e={str(e):<19} R1 {row[0]:.7f} ({agree(row[0], pkgv):4.1f}d)   "
              f"R2 {row[1]:.7f} ({agree(row[1], pkgv):4.1f}d)   pkg {pkgv:.7f}")
    # mirror check for the far defect
    _, chm, _ = live_channel("R1", [4], e=(2, 1, 0, 1, 2, 2), ch=ch6)
    print(f"  mirror check, far defect on bottom row (e=(2,1,0,1,2,2)), R1: {p_tau(chm[4], tau, 16):.7f}")

    print("\n[7] Parameter spot checks at T=4 (nu=1; tau_D^+ recomputed at each beta)")
    spot = []
    for (b, gg, kk), pl, rc in PKG["spot"]:
        taub = tau_D_plus(b)
        row = []
        for conv in ("R1", "R2"):
            _, chois, _ = live_channel(conv, [4], beta=b, g=gg, kappa=kk, ch=ch6)
            row.append((p_tau(chois[4], taub, 16), r_rank_one(chois[4], 16, 4, 0)[0]))
        spot.append(((b, gg, kk), row, pl, rc))
        print(f"  betaJ={b:<3} g={gg:<3} kappa={kk:<3} | R1 p={row[0][0]:.10f} ({agree(row[0][0], pl):4.1f}d) "
              f"r={row[0][1]:.10f} ({agree(row[0][1], rc):4.1f}d) | R2 p={row[1][0]:.10f} ({agree(row[1][0], pl):4.1f}d) "
              f"r={row[1][1]:.10f} ({agree(row[1][1], rc):4.1f}d)")
        print(f"  {'':<27}| package p={pl:.10f}  r={rc:.10f}")

    print("\n[8] Convention sensitivity, baseline, T=4 (all rows R1 thermal jumps unless labelled R2)")
    variants = [("R1 baseline", dict()),
                ("R2 baseline", dict(conv="R2")),
                ("monitoring coherence decay 2*kappa (L=sqrt(kappa) Z)", dict(mon_rate=2 * kappa)),
                ("monitoring coherence decay kappa/2", dict(mon_rate=kappa / 2)),
                ("rates with flipped field sign (anti-detailed-balance)", dict(rate_sign=-1.0)),
                ("no tie Hamiltonian, g=0", dict(g=0.0))]
    for label, kw in variants:
        kw = dict(kw)
        conv = kw.pop("conv", "R1")
        gv = kw.pop("g", g)
        _, chois, _ = live_channel(conv, [4], g=gv, ch=ch6, **kw)
        v = p_tau(chois[4], tau, 16)
        print(fmt_row(f"p_live T=4  {label}", v, PKG["p_live"][4], width=60))

    print("\n[9] Verdict")
    for conv in ("R1", "R2"):
        vals = [(res[conv]["p"][T], PKG["p_live"][T]) for T in TS] + \
               [(res[conv]["r"][T][0], PKG["r_clean"][T]) for T in (4, 8)] + \
               [(v[0 if conv == "R1" else 1][k], pv) for _, v, pl, rc in spot for k, pv in ((0, pl), (1, rc))]
        worst = min(agree(a, b) for a, b in vals)
        mx = max(abs(a - b) for a, b in vals)
        sv = [(stress[lab][0 if conv == "R1" else 1], pk) for lab, _, pk in PKG["stress"]]
        ws = min(agree(a, b) for a, b in sv)
        print(f"  {conv}: 12 ten-digit package values (p_live x4, r_clean x2, spot p/r x6): worst agreement "
              f"~{worst:.1f} sig. digits, max |diff| {mx:.1e};  stress tests (7-digit values): worst ~{ws:.1f} digits")
    wf = min(agree(pf[T], PKG["p_frozen"][T]) for T in (4, 8))
    print(f"  p_frozen (T=4,8; convention-independent): ~{wf:.1f} sig. digits")
    print(f"\n  total wall time {time.time() - t_start:.1f}s")


if __name__ == "__main__":
    main()