# SCRIPT: BH-EP-SPINOR-04
# The sign the orientation cover keeps, at the paper's joint point t=(2,1,1.3), delta=0.4.
# WHY THIS PROVES IT: biorthogonal parallel transport (u_L.v_R = 1, u_old.v_new = 1) is the connection
# a = u.dv/(u.v); its holonomy around a loop is gauge invariant for nonvanishing single-valued frames.
# (H1) closed form: with v=(A12, lam-A11), u=(A21, lam-A11), u.v = 2 sigma3^2 sqrt(Q) p, p = sqrt(Q)+c sigma3/2,
#      so a has residue exactly 1/2 at sigma3=0 and the sigma3-loop holonomy is e^{-i pi} = -1.
# (H2) numerics: w3 loops give a monodromy M with M^2 = -I, M^4 = I; a sigma3 loop gives -I.
# KILL: residue != 1/2, M^2 != -I, or the control loop not returning I.
# CONTROL: off the joint condition the small sigma3 loop returns +I.
# Also: Ivanov exchange algebra (B^2)^2 = -I, and the trace-reverse handed part at the same point.
import os as _os
_HERE=_os.path.dirname(_os.path.abspath(__file__))
_FIG=_os.path.join(_HERE,'..','paper','fig') if _os.path.isdir(_os.path.join(_HERE,'..','paper','fig')) else _os.path.join(_HERE,'..','fig')
import numpy as np, sympy as sp, matplotlib
matplotlib.use('Agg'); import matplotlib.pyplot as plt
ok=True
T=(2.0,1.0,1.3); DL=0.4
def block(w3,t=T,d=DL):
    t1,t2,t3=t; w1=(1-w3+d)/2; w2=(1-w3-d)/2
    A11=(t1+t2)-(t1-t2)*(w1-w2); A22=2*t3+w3*(t1+t2-2*t3)
    A12=((t1-t2)-(t1+t2-2*t3)*(w1-w2))/np.sqrt(3); A21=np.sqrt(3)*(t1-t2)*w3
    return np.array([[A11,A12],[A21,A22]],dtype=complex)
# ---- H1 symbolic
s=sp.symbols('sigma3'); t1,t2,t3=sp.Integer(2),sp.Integer(1),sp.Rational(13,10); d=sp.Rational(2,5)
w3=s**2; w1=(1-w3+d)/2; w2=(1-w3-d)/2
A11=(t1+t2)-(t1-t2)*(w1-w2); A22=2*t3+w3*(t1+t2-2*t3)
A12=((t1-t2)-(t1+t2-2*t3)*(w1-w2))/sp.sqrt(3); A21=sp.sqrt(3)*(t1-t2)*w3
c=t1+t2-2*t3; Q=sp.simplify((((A22-A11)/2)**2+A12*A21)/s**2); p=sp.sqrt(Q)+c*s/2
lam=(A11+A22)/2+s*sp.sqrt(Q); v=sp.Matrix([A12,lam-A11]); u=sp.Matrix([A21,lam-A11])
uv=sp.simplify((u.T*v)[0]); closed=sp.simplify(uv-2*s**2*sp.sqrt(Q)*p)
res=sp.residue(sp.simplify((u.T*v.diff(s))[0]/uv),s,0)
print('H1  Q(sigma3) =',sp.factor(Q),'  Q(0) =',Q.subs(s,0),'  A12 =',sp.nsimplify(A12))
print('    u.v - 2 sigma3^2 sqrt(Q) p =',closed,'   residue of a at the wall =',res)
ok&= closed==0 and res==sp.Rational(1,2)
# ---- H2 transport
def eig_lr(M):
    e,R=np.linalg.eig(M); return e,R,np.linalg.inv(R)
def transport(path,n=20000,record=False):
    e,R,L=eig_lr(path(0.0)); V0=R.copy(); U=L.copy(); E=e.copy(); rec=[]
    for x in np.linspace(0,1,n+1)[1:]:
        e,R,L=eig_lr(path(x))
        o=[0,1] if abs(e[0]-E[0])+abs(e[1]-E[1])<=abs(e[1]-E[0])+abs(e[0]-E[1]) else [1,0]
        Rn=R[:,o]
        for k in range(2): Rn[:,k]/=U[k]@Rn[:,k]
        U=np.linalg.inv(Rn); E=e[o]
        if record: rec.append((x,np.linalg.solve(V0,Rn)))
    return np.linalg.solve(V0,Rn),rec
r=0.02
M1,_=transport(lambda x: block(r*np.exp(2j*np.pi*x)))
M2,_=transport(lambda x: block(r*np.exp(4j*np.pi*x)),40000)
M4,_=transport(lambda x: block(r*np.exp(8j*np.pi*x)),80000)
Ms,_=transport(lambda x: block((np.sqrt(r)*np.exp(2j*np.pi*x))**2),40000)
Mc,_=transport(lambda x: block(0.3+r*np.exp(2j*np.pi*x)))
print('H2  one w3 loop: M =',np.round(M1,4).tolist(),'  M[0,1]*M[1,0] =',np.round(M1[0,1]*M1[1,0],4))
print('    M1^2 =',np.round(M1@M1,4).tolist(),'  two loops =',np.round(M2,4).tolist(),'  four loops =',np.round(M4,4).tolist())
print('    one sigma3 loop =',np.round(Ms,4).tolist(),'  control (loop off the wall point) =',np.round(Mc,4).tolist())
ok&= np.allclose(M1@M1,-np.eye(2),atol=5e-3) and np.allclose(M2,-np.eye(2),atol=5e-3) and np.allclose(M4,np.eye(2),atol=5e-3) and np.allclose(Ms,-np.eye(2),atol=5e-3) and np.allclose(Mc,np.eye(2),atol=1e-6)
Toff=(2.0,1.0,1.25)
Mo,_=transport(lambda x: block((0.05*np.exp(2j*np.pi*x))**2,t=Toff),40000)
cc=Toff[0]+Toff[1]-2*Toff[2]; ws=-((Toff[0]+Toff[1])-(Toff[0]-Toff[1])*DL-2*Toff[2])**2
print('    off the joint condition t3=1.25: small sigma3 loop =',np.round(Mo,4).tolist())
ok&= np.allclose(Mo,np.eye(2),atol=5e-3)
# ---- Ivanov exchange algebra
X=np.array([[0,1],[1,0]]);Y=np.array([[0,-1j],[1j,0]]);Z=np.diag([1,-1]);I2=np.eye(2)
g=[np.kron(X,I2),np.kron(Y,I2),np.kron(Z,X),np.kron(Z,Y)]; B=(np.eye(4)+g[0]@g[1])/np.sqrt(2)
print('IV  B^2 = g1 g2:',np.allclose(B@B,g[0]@g[1]),'  (B^2)^2 = -I:',np.allclose(np.linalg.matrix_power(B,4),-np.eye(4)),'  B^8 = I:',np.allclose(np.linalg.matrix_power(B,8),np.eye(4)))
# ---- trace reverse handed part at the same point
def Pplus(M):
    e,R=np.linalg.eig(M); L=np.linalg.inv(R); i=np.argmax(e.real); return np.real(np.outer(R[:,i],L[i,:]))
def omega(sym,s3,dd,h=1e-5):
    def P(a,b):
        M=block(a*a,d=b).real; M=(M+M.T)/2 if sym else M; return Pplus(M)
    P0=P(s3,dd); V=(P(s3+h,dd)-P(s3-h,dd))/(2*h); W=(P(s3,dd+h)-P(s3,dd-h))/(2*h); return np.trace(P0@(V@W-W@V))
Toff2=(2.0,1.0,1.25)
print('TR  handed part Tr(P+[dP/dsigma3, dP/ddelta]) at t=(2,1,1.3), delta=0.4:')
for s3 in (0.1,0.2,0.4):
    a_,b_=omega(False,s3,DL),omega(False,-s3,DL); c_=omega(True,s3,DL)
    print(f'    sigma3=+-{s3}: {a_:+.6f} / {b_:+.6f}  sum {a_+b_:.1e}   symmetrised operator: {c_:.1e}')
    ok&= abs(a_+b_)<1e-8 and abs(c_)<1e-8 and abs(a_)>1e-3
# ---- figure
BG='#0f1115'; INK='#e6e8ec'; MUT='#9aa3ae'; RULE='#2a313c'; C1='#8fb8ff'; C2='#7fd6a4'; C3='#e6c27a'; C4='#ff9b8f'
plt.rcParams.update({'figure.facecolor':BG,'axes.facecolor':BG,'savefig.facecolor':BG,'text.color':INK,'axes.labelcolor':INK,
 'xtick.color':MUT,'ytick.color':MUT,'axes.edgecolor':RULE,'font.family':'serif','font.serif':['TeX Gyre Pagella','URW Palladio L','DejaVu Serif'],'mathtext.fontset':'cm','font.size':11})
_,rec=transport(lambda x: block(r*np.exp(8j*np.pi*x)),80000,record=True)
xs=np.array([q[0] for q in rec])*4; m11=np.array([q[1][0,0] for q in rec]); m21=np.array([q[1][1,0] for q in rec])
fig,(ax,bx)=plt.subplots(1,2,figsize=(10,4.1))
for a in (ax,bx):
    for sp_ in ['top','right']: a.spines[sp_].set_visible(False)
ax.plot(xs,m11.real,color=C1,lw=2,label=r'component on $v_+$ (start)')
ax.plot(xs,m21.imag,color=C3,lw=2,label=r'Im of the component on $v_-$ (partner)')
for k in range(5): ax.axvline(k,color=RULE,lw=0.8)
ax.set_xlabel(r'loops of $w_3$ around the wall point'); ax.set_ylabel('transported eigenvector')
ax.set_xticks(range(5)); ax.set_yticks([-1,0,1])
for k,lab in zip(range(5),[r'$v_+$',r'$\propto v_-$',r'$-v_+$',r'$\propto v_-$',r'$v_+$']): ax.text(k+0.03,1.12,lab,color=INK,fontsize=10)
ax.legend(frameon=False,fontsize=9,loc='upper center',bbox_to_anchor=(0.5,-0.17),ncol=2); ax.set_ylim(-1.25,1.35)
ax.set_title('one loop swaps, two give $-1$, four return',fontsize=11,color=INK)
ss=np.linspace(0.01,0.6,200)
aj=sp.lambdify(s,sp.simplify(s*(u.T*v.diff(s))[0]/uv),'numpy')
bx.plot(ss,np.real(aj(ss)),color=C2,lw=2.2,label=r'joint condition: $\sigma_3\,a(\sigma_3)\to 1/2$')
# off joint, same closed forms numerically
def a_off(x,t=Toff2,h=1e-6):
    M=block(x*x,t=t); e,R,L=eig_lr(M); i=np.argmax(e.real); v0=np.array([M[0,1],e[i]-M[0,0]]); u0=np.array([M[1,0],e[i]-M[0,0]])
    M2=block((x+h)**2,t=t); e2,_,_=eig_lr(M2); j=np.argmin(abs(e2-e[i])); v1=np.array([M2[0,1],e2[j]-M2[0,0]])
    return x*(u0@((v1-v0)/h))/(u0@v0)
bx.plot(ss,[np.real(a_off(x)) for x in ss],color=C4,lw=2,ls='--',label=r'off the joint condition: $\to0$')
bx.axhline(0.5,color=MUT,lw=0.8,ls=':'); bx.set_xlabel(r'distance to the wall $\sigma_3$'); bx.set_ylabel(r'$\sigma_3\,a(\sigma_3)$')
bx.set_title('half a unit of connection flux at the wall',fontsize=11,color=INK); bx.legend(frameon=False,fontsize=9,loc='upper center',bbox_to_anchor=(0.5,-0.17),ncol=1)
fig.tight_layout(); fig.savefig(_os.path.join(_FIG,'f8_spinor.png'),dpi=260,bbox_inches='tight',pad_inches=0.12)
print('RESULT','PASS' if ok else 'FAIL')