# SCRIPT: BH-STIFF-19
# The kill test BH-QUTRIT-18 lacked: does equal-mass gravity stiffen toward coplanarity at EVERY flat shape?
# Random Gaussian shapes are never very elongated, so a sign change there cannot show. Here the stiffness is
# derived in closed form and then pushed toward its minimum.
# Derivation (unit masses, I = 1). A flat shape q_i (COM 0) has one out-of-plane shape direction c, with c orthogonal
# to (1,1,1,1), to the x-coordinates and to the y-coordinates. Lift q_i -> q_i + h c_i z, rescale by 1/sqrt(1+h^2|c|^2).
# Then w3 = h^2|c|^2 + O(h^4) and r_ij -> r_ij (1 + h^2 (c_i-c_j)^2/(2 r_ij^2) - h^2|c|^2/2), so for V = sum phi(r_ij)
#   dV/dw3 = (1/2) [ sum phi'(r)(c_i-c_j)^2/(r |c|^2) - sum r phi'(r) ].
# Gravity phi = -1/r:  dU/dw3 = (1/2) [ sum (c_i-c_j)^2/(r^3 |c|^2) - sum 1/r ].
# Kill conditions, written first:
#  S1  closed form against the finite-difference lift of BH-QUTRIT-18 on 200 random flat shapes: kill if rel. diff > 1e-4.
#  S2  square: closed form = 2 sqrt2 - 1: kill if diff > 1e-12.
#  S3  rhombus with diagonal ratio k = short/long: the review's claim is a sign change near k = 1/7. Kill if the
#      stiffness has no zero on (0,1); report the root k*.
#  S4  minimise the stiffness over all flat shapes (Nelder-Mead from 50 starts): kill (claim 'not general' fails) if no
#      flat shape with negative stiffness is found away from collisions (min pair distance > 0.05).
#      Run note: the first run used a soft penalty and the minimiser ran into a collision (stiffness ~ -1e8, distance 0);
#      that run is void. The barrier is now hard (shapes with a pair closer than 0.05 are rejected).
#  S5  rhombus closed form: the edge is 1/sqrt2 for every k at I = 1, so dU/dw3 = (1/2)[4 sqrt2 - sqrt((1+k^2)/2)(1+1/k)],
#      zero where (1+k^2)(1+k)^2 = 64 k^2. Kill if this root differs from S3's k* by more than 1e-10.
import numpy as np
from scipy.optimize import brentq, minimize
ok=True
def report(n,c,v=''):
    global ok; print(('PASS ' if c else 'FAIL ')+n, v); ok=ok and c
def norm(q):
    q=q-q.mean(0); return q/np.sqrt((q**2).sum())
def cdir(q):
    M=np.vstack([np.ones(4),q[:,0],q[:,1]]); _,_,vt=np.linalg.svd(M); return vt[3]
def stiff(q,sign=-1):
    q=norm(q); c=cdir(q); a=b=0
    for i in range(4):
        for j in range(i+1,4):
            r=np.linalg.norm(q[i]-q[j]); a+=(c[i]-c[j])**2/(r**3); b+=1/r
    a/= (c@c)
    return 0.5*(a-b) if sign<0 else -0.5*(a-b)
def fd(q,h=1e-4):
    q=norm(q); c=cdir(q); Q=np.c_[q,np.zeros(4)]; L=Q.copy(); L[:,2]=h*c; L/=np.sqrt((L**2).sum())
    U=lambda P:-sum(1/np.linalg.norm(P[i]-P[j]) for i in range(4) for j in range(i+1,4))
    G=L.T@L; w3=np.linalg.eigvalsh(G)[0]; return (U(L)-U(Q))/w3
rng=np.random.default_rng(19); err=0
for _ in range(200):
    q=rng.normal(size=(4,2)); err=max(err,abs(stiff(q)-fd(q))/abs(fd(q)))
report('S1 closed form = finite-difference lift',err<1e-4,f'max rel diff {err:.1e}')
sq=np.array([[.5,0],[-.5,0],[0,.5],[0,-.5]])
report('S2 square',abs(stiff(sq)-(2*np.sqrt(2)-1))<1e-12,f'{stiff(sq):.15f} vs {2*np.sqrt(2)-1:.15f}')
rh=lambda k:np.array([[1,0],[-1,0],[0,k],[0,-k]])
ks=brentq(lambda k:stiff(rh(k)),0.01,0.99)
report('S3 rhombus sign change',0<ks<1,f'k* = {ks:.6f} (1/k* = {1/ks:.4f}); stiffness at k=1 {stiff(rh(1)):+.6f}, k=0.2 {stiff(rh(0.2)):+.4f}, k=0.1 {stiff(rh(0.1)):+.4f}')
# exact root: rhombus (+-1,0),(0,+-k), c=(1,1,-1,-1)/2; I-normalisation is a common scale s=1/sqrt(2+2k^2)
def exact(k):
    s=1/np.sqrt(2+2*k*k); re=s*np.sqrt(1+k*k)
    return 0.5*(4/re**3 - 4/re - 1/(2*s) - 1/(2*s*k))
print('   exact rhombus formula at k* :', exact(ks), ' and at k=0.5:', exact(0.5), stiff(rh(0.5)))
best=(np.inf,None)
for _ in range(50):
    x0=rng.normal(size=8)
    def f(x):
        q=norm(x.reshape(4,2)); d=min(np.linalg.norm(q[i]-q[j]) for i in range(4) for j in range(i+1,4))
        return 1e6 if d<0.05 else stiff(q)
    r=minimize(f,x0,method='Nelder-Mead',options={'maxiter':4000,'xatol':1e-9,'fatol':1e-12})
    if r.fun<best[0]: best=(r.fun,norm(r.x.reshape(4,2)))
q=best[1]; d=min(np.linalg.norm(q[i]-q[j]) for i in range(4) for j in range(i+1,4))
report('S4 negative stiffness exists away from collisions',stiff(q)<0 and d>0.049,f'min found {stiff(q):+.4f} with min pair distance {d:.4f}')
import numpy.polynomial.polynomial as P
rt=[x.real for x in np.roots([1,2,-62,2,1]) if abs(x.imag)<1e-12 and 0<x.real<1]
report('S5 rhombus root from (1+k^2)(1+k)^2 = 64k^2',len(rt)==1 and abs(rt[0]-ks)<1e-10,f'k* = {rt[0]:.12f}')
fr=np.mean([stiff(rng.normal(size=(4,2)))>0 for _ in range(20000)])
print(f'   fraction of 20000 Gaussian flat shapes that stiffen: {fr:.4f}')
print('ALL PASS' if ok else 'SOME FAILED')