# SCRIPT: BH-SPINE-09
# The off-diagonal block B and its square BB^dagger, wherever a space is split into a kept part and a
# hidden part, and the contractive pullback of a passive boundary.
# Kill conditions, written before the run (each on random instances; tolerance 1e-10):
#  S1  P an orthogonal projector of rank r, split by a coordinate block: C - C^2 = B B^dagger, and
#      Tr(C^2 - C) = -Tr B B^dagger <= 0.
#  S2  a block Hamiltonian H = [[A,B],[B^dagger,D]]: the hidden-sector kernel K(t) = B e^{-iDt} B^dagger has
#      K(0) = B B^dagger, and P H^2 P - (P H P)^2 = P H Q H P = B B^dagger (Q = I - P).
#  S3  the transport operator is the compression A = P_W D_T P_W of the free scaling D_T(H) = TH + HT
#      onto trace-zero tangents, and P_W D_T^2 P_W - A^2 = P_W D_T Q_W D_T P_W is rank one and positive
#      in the kinetic metric (Q_W projects onto the scale direction W).
#  S4  passive interior L, boundary impedance Z = W^* L W by Schur, S = (Z - I)(Z + I)^{-1}. With the
#      graph map G = [I; S] and the indefinite form J = diag(I, -I): G^* J G = I - S^* S
#      = 4 (W V)^* Herm(L) (W V) with V = (Z + I)^{-1}; hence ||S|| <= 1.
#  S5  one tensor, two jobs: for a two-parameter family of pure states, Q = g + (i/2) Omega is positive
#      semidefinite, so |Omega_12| <= 2 sqrt(det g).
#  S6  one-way boundary: L_b = [[-1,0],[0.3,-b]] with exterior readout (1,0); the exterior signal is the
#      same for every hidden rate b and every hidden initial state. And P(X + QYQ)P = PXP.
#  S7  retention model x' = P - (k+a) x + b y, y' = a x - b y: total balance d(x+y)/dt = P - k x;
#      mean residence (1 + a/b)/k for b > 0; at b = 0 the retained fraction is a/(k+a).
import numpy as np, sympy as sp
from scipy.linalg import expm
rng=np.random.default_rng(9); ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
def crand(n,m): return rng.normal(size=(n,m))+1j*rng.normal(size=(n,m))
# S1
err=0
for _ in range(20):
    n,r,m=7,3,4
    X=np.linalg.qr(crand(n,r))[0]; P=X@X.conj().T
    C=P[:m,:m]; Bk=P[:m,m:]
    err=max(err,np.abs(C-C@C-Bk@Bk.conj().T).max(), abs(np.trace(C@C-C)+np.trace(Bk@Bk.conj().T)))
report('S1 C - C^2 = BB^dagger and Tr(C^2 - C) = -Tr BB^dagger', err<1e-10, f'max error {err:.1e}')
# S2
err=0
for _ in range(20):
    m,k=3,4; H=crand(m+k,m+k); H=H+H.conj().T
    A=H[:m,:m]; Bk=H[:m,m:]; D=H[m:,m:]
    K0=Bk@expm(-1j*D*0)@Bk.conj().T
    P=np.zeros((m+k,m+k)); P[:m,:m]=np.eye(m); Q=np.eye(m+k)-P
    lhs=(P@H@H@P-(P@H@P)@(P@H@P))[:m,:m]
    err=max(err,np.abs(K0-Bk@Bk.conj().T).max(),np.abs(lhs-Bk@Bk.conj().T).max(),np.abs((P@H@Q@H@P)[:m,:m]-Bk@Bk.conj().T).max())
report('S2 K(0) = BB^dagger = PH^2P - (PHP)^2 = PHQHP', err<1e-10, f'max error {err:.1e}')
# S3 in the 6-dim space of symmetric 3x3 matrices with the kinetic metric
def sym_basis():
    Bs=[]
    for i in range(3):
        M=np.zeros((3,3)); M[i,i]=1; Bs.append(M)
    for (i,j) in [(0,1),(0,2),(1,2)]:
        M=np.zeros((3,3)); M[i,j]=M[j,i]=1/np.sqrt(2); Bs.append(M)
    return Bs
Bs=sym_basis(); worst=0; ranks=[]; mineig=1
for _ in range(10):
    w=rng.dirichlet([1,1,1]); t=rng.uniform(0.5,2,3); Wm=np.diag(w); Tm=np.diag(t)
    vec=lambda M: np.array([np.sum(M*Bk) for Bk in Bs])
    G=np.array([[0.5*np.sum(Ba*Bb/(w[:,None]+w[None,:])) for Bb in Bs] for Ba in Bs])
    DT=np.column_stack([vec(Tm@Bk+Bk@Tm) for Bk in Bs])
    PW=np.column_stack([vec(Bk-np.trace(Bk)*Wm) for Bk in Bs]); QW=np.eye(6)-PW
    A=PW@DT@PW; Ad=PW@DT@DT@PW-A@A; R=PW@DT@QW@DT@PW
    worst=max(worst,np.abs(Ad-R).max(), np.abs(G@PW-PW.T@G).max())
    ranks.append(np.linalg.matrix_rank(R,tol=1e-9))
    L=np.linalg.cholesky(G); Rs=L.T@R@np.linalg.inv(L.T); mineig=min(mineig,np.linalg.eigvalsh((Rs+Rs.T)/2).min())
report('S3 A = P_W D_T P_W; P D_T^2 P - A^2 = P D_T Q D_T P, rank one, kinetic-positive', worst<1e-10 and set(ranks)=={1} and mineig>-1e-10, f'error {worst:.1e}, ranks {set(ranks)}, min eigenvalue {mineig:.1e}')
# S4
err=0; nmax=0
for _ in range(50):
    nb,ni=2,3; X=crand(nb+ni,nb+ni); K=crand(nb+ni,nb+ni); L=X@X.conj().T*0.3+(K-K.conj().T)*0.5   # Herm L >= 0: passive
    Lbb,Lbi,Lib,Lii=L[:nb,:nb],L[:nb,nb:],L[nb:,:nb],L[nb:,nb:]
    Z=Lbb-Lbi@np.linalg.solve(Lii,Lib); Wm=np.vstack([np.eye(nb),-np.linalg.solve(Lii,Lib)])
    I=np.eye(nb); S=(Z-I)@np.linalg.inv(Z+I); V=np.linalg.inv(Z+I)
    Gm=np.vstack([I,S]); Jm=np.diag([1]*nb+[-1]*nb)
    HL=(L+L.conj().T)/2
    err=max(err,np.abs(Gm.conj().T@Jm@Gm-(I-S.conj().T@S)).max(),np.abs(I-S.conj().T@S-4*(Wm@V).conj().T@HL@(Wm@V)).max(),np.abs(Wm.conj().T@L@Wm-Z).max())
    nmax=max(nmax,np.linalg.norm(S,2))
report('S4 G*JG = I - S*S = 4(WV)* Herm L (WV), ||S|| <= 1', err<1e-10 and nmax<=1+1e-12, f'error {err:.1e}, max ||S|| {nmax:.4f}')
# S5
worst=1
for _ in range(200):
    M0,M1,M2=crand(4,1),crand(4,1),crand(4,1)
    M3=crand(4,1)
    def psi(x,y):
        v=M0+x*M1+y*M2+x*y*M3; return v/np.linalg.norm(v)
    x,y=rng.normal(size=2); h=1e-6
    p=psi(x,y); dx=(psi(x+h,y)-psi(x-h,y))/(2*h); dy=(psi(x,y+h)-psi(x,y-h))/(2*h)
    Pp=np.eye(4)-p@p.conj().T; D=[dx,dy]
    Q=np.array([[ (D[a].conj().T@Pp@D[b])[0,0] for b in range(2)] for a in range(2)])
    g=Q.real; Om=-2*Q.imag
    worst=min(worst, 2*np.sqrt(max(np.linalg.det(g),0))-abs(Om[0,1]), np.linalg.eigvalsh((Q+Q.conj().T)/2).min()+1e-9)
report('S5 Q = g + (i/2) Omega >= 0, |Omega_12| <= 2 sqrt(det g)', worst>-1e-6, f'min slack {worst:.2e}')
# S6
outs=[]
for b in [0.1,1.0,7.0]:
    for h0 in [0.0,1.0,-3.0]:
        L=np.array([[-1,0],[0.3,-b]]); ts=np.linspace(0,5,11)
        outs.append([ (np.array([1,0])@expm(L*t)@np.array([1.0,h0])) for t in ts])
outs=np.array(outs); spread=np.abs(outs-outs[0]).max()
P=np.diag([1,1,0,0]); Q=np.eye(4)-P; X=rng.normal(size=(4,4)); Y=rng.normal(size=(4,4))
report('S6 exterior signal independent of hidden rate and state; P(X+QYQ)P = PXP', spread<1e-12 and np.abs(P@(X+Q@Y@Q)@P-P@X@P).max()<1e-14, f'spread {spread:.1e}')
# S7
s=sp.symbols('s'); k,a,bb,E0=sp.symbols('kappa a b E0',positive=True)
Lm=sp.Matrix([[-(k+a),bb],[a,-bb]])
Xs=(s*sp.eye(2)-Lm).inv()*sp.Matrix([E0,0]); F=sp.simplify(k*Xs[0]/E0)
tot=sp.simplify(F.subs(s,0)); mean=sp.simplify(-sp.diff(F,s).subs(s,0))
y_inf=sp.simplify(sp.limit(s*((s*sp.eye(2)-Lm.subs(bb,0)).inv()*sp.Matrix([E0,0]))[1],s,0))
report('S7 retention: escaped fraction 1, mean residence (1+a/b)/kappa; b=0 keeps a/(kappa+a)', sp.simplify(tot-1)==0 and sp.simplify(mean-(1+a/bb)/k)==0 and sp.simplify(y_inf-a*E0/(k+a))==0, f'mean = {mean}, retained(b=0) = {y_inf}')
print('ALL PASS' if ok else 'SOME FAILED')