# SCRIPT: BH-LIFT-05
# What it settles, with kill conditions written before the run:
#  K1  the transport operator is self-adjoint in the kinetic metric (shape block and full 5x5),
#      symbolically: g*A - A^T*g == 0. Kill if any entry is nonzero.
#  K2  A = P_W D_T with P_W(Z) = Z - Tr(Z) W the g-orthogonal projection onto trace-zero tangents.
#      Kill if g(W,H) != Tr(H)/4 or A != P_W D_T on random symmetric H.
#  K3  in the regular frame J = diag(1/sqrt2, -sqrt6*s) (coordinates delta, s = sigma3) the lift
#      B = J^-1 A J is analytic in s and equals lambda* I at s = 0 on the joint condition. Kill otherwise.
#  K4  complex loop of s around 0 at fixed joint data: Gram-frame biorthogonal holonomy -1,
#      regular-frame holonomy +1. Kill if not.
#  K5  the Gram residue 1/2 equals the folded-direction weight u2 v2 / (u.v) of the first-order
#      regular operator B1, which is exactly 1/2 because B1 is off-diagonal. Kill if not 1/2.
#  K6  eigenprojector norms: kinetic-metric operator norm exactly 1; Euclidean (Frobenius-frame)
#      norm times gap tends to |A12|. Kill if kinetic norm differs from 1 by more than 1e-9.
#  K7  PREDICTION: a real loop in (s, t3) around the joint point (it must pass through s<0, the
#      mirror side) returns each eigenvector with -1 in the regular frame (a conical, diabolic sign),
#      and +1 in the Gram frame (which sees the loop folded back on itself). A real loop that does
#      not enclose the joint point returns +1 in both. Kill the "cover keeps a sign" reading if the
#      regular-frame enclosing loop returns +1.
#  K8  along every real path the eigenvalues stay real on both sides of the wall (no real-to-complex
#      transition), unlike a spectral exceptional point. Kill if any complex pair appears.
import numpy as np, sympy as sp
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok = ok and cond

# ---------- K1 symbolic ----------
t1,t2,t3,w1,w2,w3=sp.symbols('t1 t2 t3 w1 w2 w3',positive=True)
E=[sp.diag(1,-1,0)/sp.sqrt(2), sp.diag(1,1,-2)/sp.sqrt(6)]
def Eij(i,j):
    M=sp.zeros(3); M[i,j]=M[j,i]=1/sp.sqrt(2); return M
E+= [Eij(0,1),Eij(0,2),Eij(1,2)]
T=sp.diag(t1,t2,t3); W=sp.diag(w1,w2,w3)
def Aop(H): return T*H+H*T-2*(T*H).trace()*W
def gW(H,K): return sp.Rational(1,2)*sum(H[i,j]*K[i,j]/(W[i,i]+W[j,j]) for i in range(3) for j in range(3))
def coords(M):
    # coordinates of traceless symmetric M in the Frobenius-orthonormal basis E
    return [sp.simplify(sum(M[i,j]*Eb[i,j] for i in range(3) for j in range(3))) for Eb in E]
sub={w3:1-w1-w2}
A5=sp.zeros(5); G5=sp.zeros(5)
for b in range(5):
    c=coords(Aop(E[b]).subs(sub))
    for a in range(5): A5[a,b]=c[a]
for a in range(5):
    for b in range(5): G5[a,b]=sp.simplify(gW(E[a],E[b]).subs(sub))
R=sp.simplify(G5*A5-A5.T*G5)
report('K1 g A - A^T g = 0 (5x5, symbolic)', R==sp.zeros(5))
# ---------- K2 ----------
rng=np.random.default_rng(5)
for trial in range(3):
    wv=rng.dirichlet([1,1,1]); tv=rng.uniform(0.5,2,3)
    Wn=np.diag(wv); Tn=np.diag(tv)
    H=rng.normal(size=(3,3)); H=H+H.T; H-=np.trace(H)/3*np.eye(3)
    gWH=0.5*sum(Wn[i,j]*H[i,j]/(wv[i]+wv[j]) for i in range(3) for j in range(3))
    DT=Tn@H+H@Tn; PD=DT-np.trace(DT)*Wn
    An=Tn@H+H@Tn-2*np.trace(Tn@H)*Wn
    report(f'K2 trial {trial}: g(W,H)=TrH/4 and A=P_W D_T', abs(gWH-np.trace(H)/4)<1e-12 and np.allclose(PD,An), f'{abs(gWH-np.trace(H)/4):.1e}')
# ---------- shape block numerics ----------
tt=np.array([2.0,1.0,1.3]); d0=0.4
def shapeA(w3c,d,t):
    t1,t2,t3=t; w1=(1-w3c+d)/2; w2=(1-w3c-d)/2
    A11=(t1+t2)-(t1-t2)*(w1-w2); A22=2*t3+w3c*(t1+t2-2*t3)
    A12=((t1-t2)-(t1+t2-2*t3)*(w1-w2))/np.sqrt(3); A21=np.sqrt(3)*(t1-t2)*w3c
    return np.array([[A11,A12],[A21,A22]],dtype=complex)
def gshape(w3c,d):
    w1=(1-w3c+d)/2; w2=(1-w3c-d)/2
    e1=np.array([1,-1,0])/np.sqrt(2); e2=np.array([1,1,-2])/np.sqrt(6); wv=np.array([w1,w2,w3c])
    f=lambda a,b: np.sum(a*b/(4*wv))
    return np.array([[f(e1,e1),f(e1,e2)],[f(e2,e1),f(e2,e2)]])
def Jf(s): return np.diag([1/np.sqrt(2), -np.sqrt(6)*s]).astype(complex)
def Breg(s,d,t):
    t1,t2,t3=t; c=t1+t2-2*t3
    # analytic closed form of J^-1 A J (checked against the product below for s != 0)
    return np.array([[-d*(t1-t2)+t1+t2, 2*s*(d*c-t1+t2)],[s*(t2-t1)/2, s*s*c+2*t3]],dtype=complex)
for s in [0.3,0.05,-0.2]:
    Bp=np.linalg.inv(Jf(s))@shapeA(s*s,d0,tt)@Jf(s)
    report(f'K3 closed form of B at s={s}', np.allclose(Bp,Breg(s,d0,tt)), f'{np.abs(Bp-Breg(s,d0,tt)).max():.1e}')
B0=Breg(0,d0,tt); lam=2.6
report('K3 B(0) = lambda* I on the joint condition', np.allclose(B0,lam*np.eye(2)), str(np.round(B0.real,12).tolist()))
# ---------- transport ----------
def transport(Mfun, path, branch0):
    # biorthogonal parallel transport of the eigenvector of the branch continued from branch0
    M=Mfun(path[0]); ev,V=np.linalg.eig(M); k=np.argmin(abs(ev-branch0))
    lamc=ev[k]; v=V[:,k]; evl,U=np.linalg.eig(M.T); u=U[:,np.argmin(abs(evl-lamc))]; u=u/(u@v)
    v0=v.copy()
    for p in path[1:]:
        M=Mfun(p); ev,V=np.linalg.eig(M); k=np.argmin(abs(ev-lamc)); lamc=ev[k]; vn=V[:,k]
        vn=vn*(u@v)/(u@vn)            # parallel: u.dv = 0
        evl,U=np.linalg.eig(M.T); un=U[:,np.argmin(abs(evl-lamc))]; un=un/(un@vn)
        v,u=vn,un
    return (v0.conj()@v)/(v0.conj()@v0)
n=40000
for frame in ['Gram','regular']:
    r=0.02; path=[r*np.exp(1j*th) for th in np.linspace(0,2*np.pi,n)]
    if frame=='Gram': Mf=lambda s: shapeA(s*s,d0,tt)
    else: Mf=lambda s: Breg(s,d0,tt)
    ev0=np.linalg.eigvals(Mf(path[0])); h=transport(Mf,path,ev0[0])
    want=-1 if frame=='Gram' else 1
    report(f'K4 complex s-loop, {frame} frame holonomy', abs(h-want)<5e-3, f'{h.real:+.5f}{h.imag:+.5f}i')
# ---------- K5 ----------
t1n,t2n,t3n=tt; B1=np.array([[0,2*(d0*(t1n+t2n-2*t3n)-t1n+t2n)],[(t2n-t1n)/2,0]])
ev,V=np.linalg.eig(B1); evl,U=np.linalg.eig(B1.T)
for k in range(2):
    v=V[:,k]; u=U[:,np.argmin(abs(evl-ev[k]))]
    report(f'K5 folded weight u2v2/(u.v), branch {k}', abs(u[1]*v[1]/(u@v)-0.5)<1e-12, f'{(u[1]*v[1]/(u@v)).real:.12f}')
# ---------- K6 ----------
for s in [0.1,0.01,0.001]:
    M=shapeA(s*s,d0,tt).real; g=gshape(s*s,d0); ev,V=np.linalg.eig(M)
    lp,lm=sorted(ev.real)[::-1]; P=(M-lm*np.eye(2))/(lp-lm)
    L=np.linalg.cholesky(g); Pk=L.T@P@np.linalg.inv(L.T)
    nk=np.linalg.norm(Pk,2); ne=np.linalg.norm(P,2)
    report(f'K6 s={s}: kinetic norm 1, Euclidean norm*gap', abs(nk-1)<1e-9, f'kinetic {nk:.12f}  Euclid {ne:.2f}  Euclid*gap {ne*(lp-lm):.6f}')
print('   |A12| =', abs(shapeA(0,d0,tt)[0,1]))
# ---------- K7 real loops in (s,t3) ----------
# On a real loop the eigenvectors are real, and the topological content is the orientation of the
# eigenline: continue a kinetic-unit eigenvector with continuous sign and compare its end with its
# start (+1 or -1 exactly). The biorthogonal transport is also printed; on real loops its connection
# is not closed (the metric varies along the loop), so only its sign is topological.
def gfull(s,frame):
    g=gshape(s*s,d0)
    if frame=='Gram': return g
    J=Jf(s).real; return J.T@g@J
def orient(frame, center, rad, n=20000):
    path=[(center[0]+rad*np.cos(th), center[1]+rad*np.sin(th)) for th in np.linspace(0,2*np.pi,n)]
    def Mf(p):
        s,t3v=p; t=(2.0,1.0,t3v)
        return (shapeA(s*s,d0,t) if frame=='Gram' else Breg(s,d0,t)).real
    M=Mf(path[0]); ev,V=np.linalg.eig(M); k=np.argmax(ev); lamc=ev[k]; v=V[:,k]
    g=gfull(path[0][0],frame); v=v/np.sqrt(v@g@v); v0=v.copy(); g0=g
    for p in path[1:]:
        M=Mf(p); ev,V=np.linalg.eig(M); k=np.argmin(abs(ev-lamc)); lamc=ev[k]; vn=V[:,k].real
        g=gfull(p[0],frame); vn=vn/np.sqrt(vn@g@vn)
        if vn@g@v<0: vn=-vn
        v=vn
    bi=transport(lambda q: Mf(q).astype(complex), path, max(np.linalg.eigvals(Mf(path[0])).real))
    return (v0@g0@v), bi
for frame,want in [('regular',-1),('Gram',1)]:
    h,bi=orient(frame,(0,1.3),0.05)
    report(f'K7 real (s,t3) loop around the joint point, {frame} frame, eigenline orientation', abs(h-want)<1e-6, f'{h:+.9f}  (biorthogonal {bi.real:+.6f})')
for frame in ['regular','Gram']:
    h,bi=orient(frame,(0.2,1.3),0.05)
    report(f'K7 real loop not enclosing it, {frame} frame, eigenline orientation', abs(h-1)<1e-6, f'{h:+.9f}  (biorthogonal {bi.real:+.6f})')
# gap along each axis of the cone
for s,dt in [(1e-3,0),(0,1e-3)]:
    ev=np.linalg.eigvals(Breg(s,d0,(2,1,1.3+dt))); print(f'   cone: s={s}, dt3={dt}: gap = {abs(ev[0]-ev[1]):.6e}')
# ---------- K8 ----------
mx=0
for t3v in np.linspace(1.1,1.5,41):
    for s in np.linspace(-0.6,0.6,241):
        ev=np.linalg.eigvals(Breg(s,d0,(2,1,t3v))); mx=max(mx,abs(ev.imag).max())
report('K8 real paths: eigenvalues real on both sides of the wall', mx<1e-12, f'max |Im| = {mx:.1e}')
print('ALL PASS' if ok else 'SOME FAILED')