# SCRIPT: BH-CORE-02
# Audit of the transport operator on the four-body shape space at the coplanar wall,
# and the braid gate BH-N-01.
#
# Objects: G = W = diag(w1,w2,w3), sum w = 1 (Gram observer at a diagonal point).
# Tangent space = traceless symmetric 3x3, Frobenius-orthonormal basis E1..E5.
# Transport operator A(H) = T H + H T - 2 Tr(T H) W, T = diag(t1,t2,t3).
#
# Why each check proves what it claims (kill conditions written before the run):
#  C1 blocks: the matrix of A in E1..E5 is computed symbolically; the claimed closed forms
#     (A11, A22, A12, A21, shear diagonal t_i+t_j, zero shape-shear coupling) are FALSE if any
#     difference does not simplify to 0.
#  C2 wall triangularity: A21 = sqrt(3)(t1-t2) w3 is FALSE if A21 at w3=0 is not identically 0.
#  C3 interior reality: the shape-block discriminant D = (A11-A22)^2 + 4 A12 A21 must be >= 0 on the
#     simplex; FALSE if a random sample of 20000 interior points gives D < -1e-12.
#  C4 Jordan condition: at w3=0 the shape block is a Jordan block iff A11=A22 and A12 != 0;
#     the claimed joint condition t1+t2-2t3 = (t1-t2)(w1-w2) is FALSE if A11-A22 at w3=0 differs from it.
#  C5 braid gate BH-N-01 (at the joint point, complexified w3 = r e^{i theta}):
#     (a) pair swaps after ONE loop in w3 and returns after TWO   -> FALSE if not;
#     (b) discriminant winds +1 per loop in w3                      -> FALSE if not;
#     (c) ledger L_c = log det(c I + A_pair) returns to its start after one loop (single-valued) -> FALSE if drift > 1e-10;
#     (d) Riesz block projector of the pair returns after one loop -> FALSE if drift > 1e-10;
#     (e) in the FOLD coordinate sigma3 (w3 = sigma3^2) one loop in sigma3 does NOT swap the pair -> FALSE if it swaps.
#  (C4 note: the identity uses w1+w2 = 1 at the wall.)  C6 ledger table (contour radius 0.25 around 2.6; the shear eigenvalue 2.3 must stay outside, so w3 <= 1e-2)
#  C6 ledger table at t=(2,1,1.3), w1-w2=0.4: contour count of the pair = 2 and Tr P_pair = 2 for w3 down to 0,
#     ledger log det(c+A_pair) continuous; FALSE if count != 2 or Tr P off by > 1e-9.
#  C7 volume structure T = G^{-1}: A is self-adjoint for <H,K>_G = tr(G^-1 H G^-1 K) -> no Jordan block possible;
#     FALSE if the G-adjoint differs from A by > 1e-9 at w3 = 1e-6.
#  C8 transient sizing: ||exp(-J t)|| for J = [[b,mu],[0,b]] equals e^{-bt}(sqrt(4+mu^2 t^2)+|mu| t)/2; FALSE if off by > 1e-10.
import numpy as np, sympy as sp, cmath, math

t1,t2,t3,w1,w2,w3 = sp.symbols('t1 t2 t3 w1 w2 w3', real=True)
def basis():
    E1 = sp.diag(1,-1,0)/sp.sqrt(2); E2 = sp.diag(1,1,-2)/sp.sqrt(6)
    E=[E1,E2]
    for (i,j) in [(0,1),(0,2),(1,2)]:
        Z=sp.zeros(3,3); Z[i,j]=Z[j,i]=1/sp.sqrt(2); E.append(Z)
    return E
E = basis()
T = sp.diag(t1,t2,t3); W = sp.diag(w1,w2,w3)
def A_of(H): return T*H + H*T - 2*(T*H).trace()*W
M = sp.Matrix(5,5, lambda a,b: sp.simplify((E[a].T*A_of(E[b])).trace()))
Msub = sp.simplify(M.subs(w2, 1-w1-w3))
ok = True
claims = {
 (0,0): (t1+t2)-(t1-t2)*(w1-w2),
 (1,1): 2*t3 + w3*(t1+t2-2*t3),
 (0,1): ((t1-t2)-(t1+t2-2*t3)*(w1-w2))/sp.sqrt(3),
 (1,0): sp.sqrt(3)*(t1-t2)*w3,
 (2,2): t1+t2, (3,3): t1+t3, (4,4): t2+t3,
}
print('C1 matrix of A in E1..E5 (rows a, cols b), with w2 = 1 - w1 - w3:')
for (a,b),val in claims.items():
    d = sp.simplify(Msub[a,b] - sp.simplify(sp.sympify(val).subs(w2,1-w1-w3)))
    print(f'   A[{a+1}{b+1}] claim residual: {d}'); ok &= (d==0)
coup = [sp.simplify(Msub[a,b]) for a in range(2) for b in range(2,5)] + [sp.simplify(Msub[a,b]) for a in range(2,5) for b in range(2)]
offsh = [sp.simplify(Msub[a,b]) for a in range(2,5) for b in range(2,5) if a!=b]
print('   shape-shear coupling entries:', set(coup), ' shear off-diagonal:', set(offsh)); ok &= set(coup)=={0} and set(offsh)=={0}
print('C2 A21 at w3=0:', sp.simplify(Msub[1,0].subs(w3,0))); ok &= sp.simplify(Msub[1,0].subs(w3,0))==0
D = sp.expand((M[0,0]-M[1,1])**2 + 4*M[0,1]*M[1,0])
rng = np.random.default_rng(1); Df = sp.lambdify((t1,t2,t3,w1,w2,w3), D, 'numpy'); mn = 1e9
for _ in range(20000):
    w = rng.dirichlet([1,1,1]); t = rng.uniform(0.1,5,3); mn = min(mn, Df(*t,*w))
print(f'C3 min discriminant over 20000 interior samples: {mn:.3e}'); ok &= mn > -1e-12
diff0 = sp.simplify((M[0,0]-M[1,1]).subs(w3,0).subs(w2,1-w1))
print('C4 (A11 - A22) at w3=0 =', sp.factor(diff0), ' ; claim residual:', sp.simplify(diff0 - ((t1+t2-2*t3)-(t1-t2)*(w1-w2)).subs(w2,1-w1)))
ok &= sp.simplify(diff0 - ((t1+t2-2*t3)-(t1-t2)*(w1-w2)).subs(w2,1-w1))==0

# numeric helpers
def Anum(t, w):
    f = sp.lambdify((t1,t2,t3,w1,w2,w3), M, 'numpy')
    return np.array(f(*t,*w), dtype=complex)
Mf = sp.lambdify((t1,t2,t3,w1,w2,w3), M, 'numpy')
def shape_block(t, dlt, z):  # w1-w2 = dlt fixed, w3 = z (complex allowed)
    wa = (1 - z + dlt)/2; wb = (1 - z - dlt)/2
    return np.array(Mf(*t, wa, wb, z), dtype=complex)[:2,:2]
t = (2.0, 1.0, 1.3)
dstar = (t[0]+t[1]-2*t[2])/(t[0]-t[1])     # joint condition at w3=0
print(f'C5 joint point: w1-w2 = {dstar}, A12 at wall = {shape_block(t,dstar,0)[0,1].real:.6f}')
def track(path_fn, n=4000, loops=1):
    th = np.linspace(0, 2*np.pi*loops, n*loops+1)
    B0 = path_fn(th[0]); ev = np.linalg.eigvals(B0); cur = ev.copy(); start = ev.copy()
    for k in range(1,len(th)):
        e = np.linalg.eigvals(path_fn(th[k]))
        if abs(e[0]-cur[0])+abs(e[1]-cur[1]) > abs(e[1]-cur[0])+abs(e[0]-cur[1]): e = e[::-1]
        cur = e
    return start, cur
r = 1e-3
pw = lambda th: shape_block(t, dstar, r*np.exp(1j*th))
s1,e1 = track(pw,loops=1); s2,e2 = track(pw,loops=2)
swap1 = abs(e1[0]-s1[1]) < 1e-9 and abs(e1[1]-s1[0]) < 1e-9
ret2 = abs(e2[0]-s2[0]) < 1e-9 and abs(e2[1]-s2[1]) < 1e-9
print(f'C5a one loop in w3: swapped={swap1}; two loops: returned={ret2}'); ok &= swap1 and ret2
th = np.linspace(0,2*np.pi,8001)
Dv = np.array([ (lambda B: (B[0,0]-B[1,1])**2 + 4*B[0,1]*B[1,0])(pw(x)) for x in th])
wind = np.sum(np.diff(np.unwrap(np.angle(Dv))))/(2*np.pi)
print(f'C5b discriminant winding per loop: {wind:.6f}'); ok &= abs(wind-1) < 1e-6
c = 1.0
Lc = lambda B: np.log(np.linalg.det(c*np.eye(2)+B))
L0, L1 = Lc(pw(0)), Lc(pw(2*np.pi));
# continuity along the loop (no branch jump of the log because det(c+B) stays near (c+2.6)^2)
Ls = np.array([Lc(pw(x)) for x in th]); jump = np.max(np.abs(np.diff(Ls)))
print(f'C5c ledger log det(cI+A_pair): start {L0:.12f}, after loop {L1:.12f}, drift {abs(L1-L0):.2e}, max step {jump:.2e}')
ok &= abs(L1-L0) < 1e-10
# Riesz block projector of the pair (whole shape block -> identity on the 2-dim block of the 5x5)
def riesz(Afull, center, rad, n=2000):
    zs = center + rad*np.exp(1j*np.linspace(0,2*np.pi,n,endpoint=False)); P = np.zeros_like(Afull)
    for z in zs: P += np.linalg.inv(z*np.eye(Afull.shape[0]) - Afull)*(z-center)
    return P/n
def full(z, dlt=dstar):
    wa=(1-z+dlt)/2; wb=(1-z-dlt)/2; return np.array(Mf(*t,wa,wb,z),dtype=complex)
P0 = riesz(full(r), 2.6, 0.2); P1 = riesz(full(r*np.exp(1j*np.pi)), 2.6, 0.2); P2 = riesz(full(r*np.exp(2j*np.pi)), 2.6, 0.2)
print(f'C5d Riesz pair projector: Tr = {np.trace(P0).real:.12f}; drift after loop {np.max(np.abs(P2-P0)):.2e}; half-loop differs by {np.max(np.abs(P1-P0)):.2e} (moves, but returns)')
ok &= np.max(np.abs(P2-P0)) < 1e-10
ps = lambda th: shape_block(t, dstar, (np.sqrt(r)*np.exp(1j*th))**2)
s3,e3 = track(ps, loops=1)
swapS = abs(e3[0]-s3[1]) < 1e-9 and abs(e3[1]-s3[0]) < 1e-9
print(f'C5e one loop in fold coordinate sigma3 (w3 = sigma3^2, encircles w3 twice): swapped={swapS}'); ok &= not swapS
# real sigma3 line: eigenvalues analytic in sigma3, exchanged by sigma3 -> -sigma3, eigenvectors coalesce at 0
for s in [0.02, -0.02]:
    B = shape_block(t, dstar, s*s); ev, V = np.linalg.eig(B); o = np.argsort(ev.real)
    print(f'     sigma3={s:+.2f}: eigenvalues {np.round(ev[o].real,8)}, |cos angle between eigenvectors| = {abs(np.vdot(V[:,0],V[:,1]))/np.linalg.norm(V[:,0])/np.linalg.norm(V[:,1]):.6f}')
gap = lambda s: (lambda e: abs(e[0]-e[1]))(np.linalg.eigvals(shape_block(t,dstar,s*s)))
print('     gap / |sigma3| at sigma3 = 1e-2, 1e-3, 1e-4:', [round(gap(s)/s,8) for s in [1e-2,1e-3,1e-4]], ' (constant => gap linear in sigma3, square-root in w3)')

print('C6 ledger table t=(2,1,1.3), w1-w2=0.4 (joint point):')
for z in [1e-2,1e-3,1e-4,1e-6,0.0]:
    Af = full(z, 0.4); P = riesz(Af, 2.6, 0.25); ev = np.sort(np.linalg.eigvals(Af[:2,:2]).real)
    cnt = sum(1 for e in np.linalg.eigvals(Af) if abs(e-2.6)<0.25)
    Ldg = np.log(np.linalg.det(np.eye(2)*1.0 + Af[:2,:2])).real
    Pd = riesz(Af[:2,:2], ev[0], abs(ev[1]-ev[0])/2) if z>0 else None
    nP1 = np.linalg.norm(Pd,2) if Pd is not None else float('inf')
    print(f'   w3={z:7.0e}: pair {ev}, count {cnt}, Tr P_pair {np.trace(P).real:.12f}, ledger log det(1+A) {Ldg:.12f}, ||P_single||*gap = {nP1*abs(ev[1]-ev[0]) if z>0 else float("nan"):.6f}')
    ok &= cnt==2 and abs(np.trace(P).real-2)<1e-9
print(f'   |A12| at wall = {abs(full(0,0.4)[0,1]):.6f}  (individual projector norm x gap -> |A12|)')

# C7 volume structure
w = np.array([0.6, 0.4-1e-6, 1e-6]); Gm = np.diag(w); Gi = np.diag(1/w)
En = [np.array(e.tolist(), dtype=float) for e in E]
def A_T(H, Tm): return Tm@H + H@Tm - 2*np.trace(Tm@H)*Gm
Amat = np.array([[np.trace(En[a].T @ A_T(En[b], Gi)) for b in range(5)] for a in range(5)])
Gmet = np.array([[np.trace(Gi@En[a]@Gi@En[b]) for b in range(5)] for a in range(5)])
# self-adjoint for Gmet  <=>  Gmet @ Amat symmetric
S = Gmet@Amat; asym = np.max(np.abs(S-S.T))/np.max(np.abs(S))
print(f'C7 T=G^-1 at w3=1e-6: relative asymmetry of G-metric * A = {asym:.2e}  (0 => self-adjoint => diagonalisable, no Jordan block)')
ok &= asym < 1e-9
for w3v in [1e-2,1e-4,1e-6]:
    w = np.array([0.6, 0.4-w3v, w3v]); Gm = np.diag(w); Gi=np.diag(1/w)
    Am = np.array([[np.trace(En[a].T @ (Gi@En[b] + En[b]@Gi - 2*np.trace(Gi@En[b])*Gm)) for b in range(5)] for a in range(5)])
    ev = np.sort(np.abs(np.linalg.eigvals(Am)))
    print(f'     w3={w3v:.0e}: |eigenvalues| = {np.round(ev,4)}  (x w3: {np.round(ev*w3v,4)})')

# C8 transient formula
b, mu = 2.6, 0.484974
for tt in [0.1, 1.0, 3.0]:
    J = np.array([[b,mu],[0,b]]); from scipy.linalg import expm
    nrm = np.linalg.norm(expm(-J*tt),2); frm = math.exp(-b*tt)*(math.sqrt(4+mu*mu*tt*tt)+abs(mu)*tt)/2
    print(f'C8 t={tt}: ||exp(-Jt)|| = {nrm:.12f}, formula {frm:.12f}'); ok &= abs(nrm-frm) < 1e-10
print(f'   growth ratio kappa = |mu|/b = {mu/b:.4f} (transient growth above 1 needs kappa > 2)')
print('RESULT', 'PASS' if ok else 'FAIL')