# SCRIPT: BH-WALL-CHECK-01
# Purpose: test whether the coplanar-wall "blow-ups" of blkhole v2 are geometric or coordinate effects.
# Kill conditions written before the run:
#  K1: if the Kendall d=2 control does not return sectional curvature 4, the machinery is wrong -> stop.
#  K2: claim "metric blow-up is a coordinate effect" FAILS if, in amplitude coordinates sigma_i=sqrt(w_i),
#      the shape block is not the identity (flat) to 1e-8.
#  K3: claim "wall at finite distance" FAILS if the length of the path W(w)=diag((1-w)/2,(1-w)/2,w), w from 0 to eps,
#      differs from arcsin(sqrt(eps)) by more than 1e-6.
#  K4: claim "curvature divergence is a normalisation artefact" FAILS if the metric-normalised sectional curvatures
#      K(E4,V), K(E5,V), K(E2,V) with V metric-unit grow like 1/w3 (ratio between w3=1e-2 and 1e-4 above 10).
import numpy as np
def basis3():
    E1=np.diag([1,-1,0])/np.sqrt(2); E2=np.diag([1,1,-2])/np.sqrt(6)
    E=[E1,E2]
    for (i,j) in [(0,1),(0,2),(1,2)]:
        Z=np.zeros((3,3)); Z[i,j]=Z[j,i]=1/np.sqrt(2); E.append(Z)
    return E
def skew(d):
    out=[]
    for i in range(d):
        for j in range(i+1,d):
            O=np.zeros((d,d)); O[i,j]=1; O[j,i]=-1; out.append(O)
    return out
def lift(X,E):
    d=X.shape[0]; cols=[(k,l) for k in range(d) for l in range(d)]
    M=[];b=[]
    for i in range(d):
        for j in range(i,d):
            row=[]
            for (k,l) in cols:
                H=np.zeros((d,d));H[k,l]=1; row.append((X.T@H+H.T@X)[i,j])
            M.append(row); b.append(E[i,j])
    for i in range(d):
        for j in range(i+1,d):
            row=[]
            for (k,l) in cols:
                H=np.zeros((d,d));H[k,l]=1; row.append((X@H.T-H@X.T)[i,j])
            M.append(row); b.append(0.0)
    return np.linalg.solve(np.array(M),np.array(b)).reshape(d,d)
def vert(X,Z):
    Vb=[O@X for O in skew(X.shape[0])]
    G=np.array([[np.sum(a*b) for b in Vb] for a in Vb]); r=np.array([np.sum(a*Z) for a in Vb])
    c=np.linalg.solve(G,r); return sum(ci*vi for ci,vi in zip(c,Vb))
def Aten(X,H1,E2,eps):
    Xp=X+eps*H1; Xp/=np.linalg.norm(Xp); Xm=X-eps*H1; Xm/=np.linalg.norm(Xm)
    D=(lift(Xp,E2)-lift(Xm,E2))/(2*eps); D=D-np.sum(X*D)*X
    return vert(X,D)
def sectional(X,Ea,Eb,eps):
    Ha,Hb=lift(X,Ea),lift(X,Eb)
    gaa,gbb,gab=np.sum(Ha*Ha),np.sum(Hb*Hb),np.sum(Ha*Hb)
    A=Aten(X,Ha,Eb,eps)
    return (gaa*gbb-gab**2+3*np.sum(A*A))/(gaa*gbb-gab**2), gaa, gbb
ok=True
# K1 Kendall
for (a,b) in [(0.7,0.3),(0.5,0.5),(0.9,0.1)]:
    X=np.diag([np.sqrt(a),np.sqrt(b)]); F1=np.diag([1,-1])/np.sqrt(2); F2=np.array([[0,1],[1,0]])/np.sqrt(2)
    K,_,_=sectional(X,F1,F2,1e-6); print(f"K1 Kendall w=({a},{b}) sectional {K:.7f}"); ok&=abs(K-4)<1e-4
# K2 amplitude coordinates: shape metric g_ij = sum dw_i dw_j/(4 w_i) -> in sigma: dw=2 sigma dsigma -> identity
w=np.array([0.5,0.3,0.2]); s=np.sqrt(w)
J=np.diag(2*s)   # dw/dsigma
gw=np.diag(1/(4*w)); gs=J.T@gw@J; print("K2 shape metric in sigma coords:\n",gs); ok&=np.allclose(gs,np.eye(3),atol=1e-8)
# K3 finite length
from scipy.integrate import quad
for eps in [1e-1,1e-2,1e-4]:
    L,_=quad(lambda w: 1/np.sqrt(4*w*(1-w)),0,eps); print(f"K3 eps={eps:g} length {L:.8f} arcsin sqrt {np.arcsin(np.sqrt(eps)):.8f}"); ok&=abs(L-np.arcsin(np.sqrt(eps)))<1e-6
# K4 sectional curvature near the wall, metric-normalised
E=basis3()
Vc=np.array([0.3,-0.2,0.4,0.1,-0.5])
res={}
for w3 in [1e-2,1e-3,1e-4]:
    w=np.array([0.55*(1-w3),0.45*(1-w3),w3]); X=np.diag(np.sqrt(w))
    V=sum(c*e for c,e in zip(Vc,E))
    row=[]
    for idx in [1,3,4]:
        K,_,gV=sectional(X,E[idx],V,min(1e-6,w3*1e-2)); row.append(K)
    res[w3]=row; print(f"K4 w3={w3:g}: sectional K(E2,V)={row[0]:.5f} K(E4,V)={row[1]:.5f} K(E5,V)={row[2]:.5f}")
ratio=max(abs(res[1e-4][i])/abs(res[1e-2][i]) for i in range(3)); print("K4 growth ratio 1e-2 -> 1e-4:",ratio); ok&=ratio<10
print("RESULT", "PASS" if ok else "FAIL")