# SCRIPT: BH-QUTRIT-18
# Tests the review's proposed closings of two open problems before anything is written into the paper.
# Shape: four equal masses in R^3, mass-weighted Jacobi matrix X (3x3), G = X^T X, Tr G = 1 (unit moment of inertia).
# Kinetic metric of the shape = the norm of the horizontal velocity: g(Gdot,Gdot) = min |Xdot|_F^2 subject to
# X^T Xdot + Xdot^T X = Gdot (computed here as a minimum-norm least-squares solution, no formula assumed).
#
# Claims under test, with kill conditions written first:
#  Q1  quantum Fisher information of the density matrix rho = G along any curve equals 4 g(Gdot,Gdot), so it is
#      8 times the shape kinetic energy (1/2)|Xdot|^2 at unit moment of inertia (the review said 4 times; tested).
#      SLD formula F_Q = 2 sum_ij |<i|Gdot|j>|^2/(p_i+p_j). Kill if F_Q/g differs from 4 by more than 1e-8.
#  Q2  the paper's wall distance arcsin sqrt(eps) (Theorem 2) is the Bures angle arccos sqrt(F) between the coplanar
#      configuration diag(1/2,1/2,0) and diag((1-eps)/2,(1-eps)/2,eps). Kill if they differ by more than 1e-12.
#  Q3  a configuration and its mirror image have fidelity exactly 1 (same G) while their distance on the
#      orientation cover is 2 arcsin|sigma3| for the fold path: fidelity is the Gram observer. Kill if F != 1.
#  P1  a thin-disk external potential V = (1/2) nu^2 sum_i m_i z_i^2 (COM in the disk plane), minimised over the
#      orientation of the configuration, equals (1/2) nu^2 w3 = (1/2) nu^2 sigma3^2: the transverse frequency in the
#      kinetic metric is exactly nu. Kill if min over orientations differs from w3 by more than 1e-10.
#  P2  any potential built from mutual distances is a function of G, hence even in sigma3; its transverse
#      stiffness at the wall is 2 dV/dw3. For Newtonian gravity and for repulsive Coulomb, equal masses/charges,
#      report the sign of dV/dw3 at the wall on random planar configurations (does pair gravity stiffen? measured).
#      Run note: the first run's kill condition predicted that neither stiffens; it failed for gravity (every sampled
#      configuration stiffens). The condition below is the corrected statement, with an exact control: for the
#      square, lifting alternate corners at fixed I gives U = U0 + (2 sqrt2 - 1) w3 + O(w3^2) (derived by hand).
#  C1  closed form of the cone slope on the joint condition: gap/|sigma3| -> 2|t1-t2| sqrt(1-delta^2)
#      (from D = 4 (t1-t2)^2 (1-delta^2) w3 + O(w3^2)); compare with the paper's 1.833030 at t=(2,1,1.3), delta=0.4.
import numpy as np
rng=np.random.default_rng(18)
ok=True
def report(name,cond,val=''):
    global ok
    print(('PASS ' if cond else 'FAIL ')+name, val); ok=ok and cond
def kin(X,Gd):
    # unknown Xdot (9 entries); constraint X^T Xdot + Xdot^T X = Gd (6 independent equations)
    rows=[];rhs=[]
    for a in range(3):
        for b in range(a,3):
            r=np.zeros((3,3))
            # (X^T Xd)_ab + (Xd^T X)_ab = sum_k X_ka Xd_kb + Xd_ka X_kb
            for k in range(3):
                r[k,b]+=X[k,a]; r[k,a]+=X[k,b]
            rows.append(r.ravel()); rhs.append(Gd[a,b])
    A=np.array(rows); xd=np.linalg.lstsq(A,np.array(rhs),rcond=None)[0]
    return float(xd@xd), xd.reshape(3,3)
def qfi(G,Gd):
    p,U=np.linalg.eigh(G); M=U.T@Gd@U; s=0
    for i in range(3):
        for j in range(3):
            if p[i]+p[j]>1e-14: s+=2*M[i,j]**2/(p[i]+p[j])
    return s
rat=[]
for _ in range(200):
    X=rng.normal(size=(3,3)); X/=np.sqrt(np.trace(X.T@X))
    Xd=rng.normal(size=(3,3)); Xd-=np.trace(X.T@Xd)*X   # keep Tr G = 1 to first order
    G=X.T@X; Gd=X.T@Xd+Xd.T@X
    g,_=kin(X,Gd); rat.append(qfi(G,Gd)/g)
rat=np.array(rat)
report('Q1 F_Q = 4 g (so F_Q = 8 x shape kinetic energy at unit inertia)', np.abs(rat-4).max()<1e-8, f'F_Q/g over 200 random curves: {rat.min():.12f} .. {rat.max():.12f}')
def fid(A,B):
    w,V=np.linalg.eigh(A); sA=V@np.diag(np.sqrt(np.clip(w,0,None)))@V.T
    m=sA@B@sA; w2=np.linalg.eigvalsh(m); return np.sum(np.sqrt(np.clip(w2,0,None)))**2
err=0
for eps in [1e-6,1e-3,0.05,0.3,0.9]:
    G0=np.diag([.5,.5,0.]); G1=np.diag([(1-eps)/2,(1-eps)/2,eps])
    err=max(err,abs(np.arccos(np.sqrt(fid(G0,G1)))-np.arcsin(np.sqrt(eps))))
report('Q2 wall distance arcsin sqrt(eps) = Bures angle arccos sqrt(F)', err<1e-12, f'max difference {err:.1e}')
X=rng.normal(size=(3,3)); X/=np.sqrt(np.trace(X.T@X)); Xm=np.diag([1,1,-1])@X
F=fid(X.T@X,Xm.T@Xm); s3=np.sqrt(np.linalg.eigvalsh(X.T@X)[0])
report('Q3 mirror pair: fidelity 1, cover distance 2 arcsin|sigma3| > 0', abs(F-1)<1e-12 and s3>1e-3, f'F = {F:.15f}, sign det X = {np.sign(np.linalg.det(X)):+.0f} vs {np.sign(np.linalg.det(Xm)):+.0f}, cover distance on the fold path 2 arcsin|sigma3| = {2*np.arcsin(s3):.4f}')
# P1 disk potential
def rot(v):
    v=v/np.linalg.norm(v); a=np.array([0,0,1.]); c=np.cross(a,v); s=np.linalg.norm(c); cc=a@v
    if s<1e-12: return np.eye(3)
    K=np.array([[0,-c[2],c[1]],[c[2],0,-c[0]],[-c[1],c[0],0]])/s
    return np.eye(3)+s*K+(1-cc)*K@K
err=0
for _ in range(50):
    X=rng.normal(size=(3,3)); X/=np.sqrt(np.trace(X.T@X)); w3=np.linalg.eigvalsh(X.T@X)[0]
    # sum_i m_i z_i^2 over bodies with COM at z=0 equals the zz entry of X X^T (mass-weighted Jacobi identity)
    best=np.inf
    for _ in range(4000):
        R=rot(rng.normal(size=3)); Y=R@X; best=min(best,(Y@Y.T)[2,2])
    M=X@X.T; best=min(best,np.linalg.eigvalsh(M)[0])  # the aligned orientation (eigenvector of X X^T along z)
    err=max(err,abs(best-w3))
# the identity sum m z^2 = (X X^T)_zz checked on explicit bodies
q=rng.normal(size=(4,3)); q-=q.mean(0)
J=np.array([(q[1]-q[0])/np.sqrt(2),(q[2]-(q[0]+q[1])/2)*np.sqrt(2/3),(q[3]-(q[0]+q[1]+q[2])/3)*np.sqrt(3/4)]).T
idz=abs((J@J.T)[2,2]-np.sum(q[:,2]**2))
report('P1 disk potential: min over orientation of sum m z^2 = w3 = sigma3^2 (frequency nu exactly)', err<1e-10 and idz<1e-12, f'max |min - w3| {err:.1e}; Jacobi identity {idz:.1e}')
# P2 pair potentials near the wall
def bodies(X):
    r1,r2,r3=X[:,0],X[:,1],X[:,2]
    q1=-r1/np.sqrt(2)-r2/np.sqrt(6)-r3/np.sqrt(12); q2=r1/np.sqrt(2)-r2/np.sqrt(6)-r3/np.sqrt(12)
    q3=2*r2/np.sqrt(6)-r3/np.sqrt(12); q4=3*r3/np.sqrt(12)
    return [q1,q2,q3,q4]
def U(X,kind):
    q=bodies(X); s=0
    for i in range(4):
        for j in range(i+1,4):
            r=np.linalg.norm(q[i]-q[j]); s+= -1/r if kind=='grav' else 1/r
    return s
signs={'grav':[],'coul':[]}
for _ in range(400):
    P=rng.normal(size=(3,3)); P[2,:]=0                 # planar configuration (third row zero): sigma3 = 0
    P/=np.sqrt(np.trace(P.T@P))
    u,sv,vt=np.linalg.svd(P); n=u[:,2]                 # out-of-plane direction
    h=1e-4
    for kind in signs:
        # move out of the plane along the third singular direction, renormalise to Tr G = 1
        Xp=P+h*np.outer(n,vt[2]); Xp/=np.sqrt(np.trace(Xp.T@Xp))
        w3=np.linalg.eigvalsh(Xp.T@Xp)[0]
        signs[kind].append((U(Xp,kind)-U(P,kind))/w3)
for kind in signs:
    a=np.array(signs[kind]); print(f'   {kind}: dV/dw3 at the wall over 400 planar configurations: min {a.min():+.4f}, max {a.max():+.4f}, fraction stiffening (>0) {np.mean(a>0):.3f}')
dg=np.array(signs['grav']); dc=np.array(signs['coul'])
# square control
def Usq(h):
    a=np.sqrt(0.25-h*h); q=[np.array([a,0,h]),np.array([-a,0,h]),np.array([0,a,-h]),np.array([0,-a,-h])]
    return -sum(1/np.linalg.norm(q[i]-q[j]) for i in range(4) for j in range(i+1,4)), 4*h*h
u0,_=Usq(0); hh=1e-4; u1,w3=Usq(hh); slope=(u1-u0)/w3
report('P2 measured: equal-mass gravity stiffens toward coplanarity at every sampled planar configuration; like charges do the opposite', np.all(dg>0) and np.all(dc<0) and abs(slope-(2*np.sqrt(2)-1))<1e-6,
       f'square control dU/dw3 = {slope:.8f} vs 2 sqrt2 - 1 = {2*np.sqrt(2)-1:.8f}; transverse omega = sqrt(2 dU/dw3) = {np.sqrt(2*(2*np.sqrt(2)-1)):.6f} (G = m = I = 1)')
t1,t2,t3,dl=2,1,1.3,0.4
def gap(s3):
    w3=s3*s3; A=np.array([[(t1+t2)-(t1-t2)*dl,(1/np.sqrt(3))*((t1-t2)-(t1+t2-2*t3)*dl)],[np.sqrt(3)*(t1-t2)*w3,2*t3+w3*(t1+t2-2*t3)]])
    e=np.linalg.eigvals(A); return abs(e[0]-e[1])
cf=2*abs(t1-t2)*np.sqrt(1-dl*dl)
report('C1 cone slope closed form 2|t1-t2| sqrt(1-delta^2)', abs(gap(1e-5)/1e-5-cf)<1e-6 and abs(cf-1.833030)<1e-6, f'{cf:.9f} vs measured {gap(1e-5)/1e-5:.9f}')
print('ALL PASS' if ok else 'SOME FAILED')