# SCRIPT: V10_PUSH_OTHERS
# Each result headed for v10 is checked a second and third way, with an opposite control.
# Kill lines are stated in each block before the computation.
import numpy as np, sympy as sp
from scipy.integrate import solve_ivp
from scipy.linalg import expm
rng = np.random.default_rng(1910)
verdict = {}

# ---------- P1 photon ring: e^{-pi} per half orbit, delay 3 sqrt3 pi M ----------
# Way 1 (symbolic, sympy): the wobble generator of the null orbit, as in PHOTON_RING_WALL.
# Way 2 (orbit equation): u'' + u = 3 M u^2 for null rays (u = 1/r, ' = d/dphi). Linearise at u0 = 1/(3M):
#   delta'' = (6 M u0 - 1) delta = delta, so growth per radian is exactly 1 and per half orbit e^{pi}.
#   We integrate the NONLINEAR equation from a tiny offset and fit the growth per radian.
# Way 3 (opposite): Reissner-Nordstrom null orbit, u'' + u = 3 M u^2 - 2 Q^2 u^3. The growth per radian
#   is sqrt(6 M u0 - 6 Q^2 u0^2 - 1) at the photon sphere: not 1 when Q != 0, so e^{-pi} is Schwarzschild's.
# Kill: way 2 rate differs from 1 by > 1e-4; way 3 at Q = 0.5 gives rate 1 within 1e-3.
M = 1.0
def growth_rate(Q):
    u0 = (3*M + np.sqrt(9*M**2 - 8*Q**2))/(4*Q**2) if Q > 0 else 1/(3*M)
    u0 = 1/((3*M + np.sqrt(9*M**2 - 8*Q**2))/2)      # photon sphere radius r = (3M + sqrt(9M^2 - 8Q^2))/2
    f = lambda ph, s: [s[1], -s[0] + 3*M*s[0]**2 - 2*Q**2*s[0]**3]
    d0 = 1e-9
    sol = solve_ivp(f, [0, 12], [u0 + d0, 0.0], method='DOP853', rtol=1e-13, atol=1e-16, dense_output=True)
    ph = np.linspace(4, 11, 50); dv = np.abs(sol.sol(ph)[0] - u0)
    rate = np.polyfit(ph, np.log(dv), 1)[0]
    pred = np.sqrt(6*M*u0 - 6*Q**2*u0**2 - 1)
    return rate, pred
r0, p0 = growth_rate(0.0); r5, p5 = growth_rate(0.5)
verdict['P1'] = abs(r0 - 1) < 1e-4 and abs(r5 - 1) > 1e-3
print(f"P1 photon ring, nonlinear orbit: growth per radian {r0:.6f} (exact 1) -> per half orbit e^pi, echo factor {np.exp(-np.pi*r0):.6f}")
print(f"   opposite, charged hole Q = 0.5 M: growth per radian {r5:.6f} (formula {p5:.6f}), echo factor {np.exp(-np.pi*r5):.6f}  {'PASS' if verdict['P1'] else 'KILL'}")
# delay: coordinate time per half orbit at r = 3M is pi / Omega, Omega = sqrt(M/r^3)
print(f"   half-orbit delay pi/Omega = {np.pi/np.sqrt(M/27):.6f} M against 3 sqrt3 pi M = {3*np.sqrt(3)*np.pi:.6f} M")

# ---------- P2 sine guard: |Omega(V, d_x)| <= 2 sqrt(g_par - C) = 2 sqrt(g_par) cos(theta) ----------
# Why: the QGT positivity bound |Omega_12| <= 2 sqrt(det g) applied to the pair (V, d_x), whose Gram
# determinant is g_par - C. Way 1: two-level families (equality expected: rank-one QGT on CP^1).
# Way 2: rank-one projectors in C^3 from random 3x3 Hermitian families (strict inequality expected).
# Kill: any violation above 1e-9, or the Gram identity off by > 1e-8.
def proj_family(n):
    Hs = [rng.normal(size=(n, n)) + 1j*rng.normal(size=(n, n)) for _ in range(3)]
    Hs = [(h + h.conj().T)/2 for h in Hs]
    def P(x, y):
        H = Hs[0] + np.cos(2*x)*Hs[1] + np.sin(2*x)*y*Hs[2] + y**2*Hs[2]
        w, v = np.linalg.eigh(H); return np.outer(v[:, 0], v[:, 0].conj())
    return P
def guard(n, trials):
    worst_id, worst_b, slack = 0.0, -np.inf, []
    for _ in range(trials):
        P = proj_family(n); x, y = rng.uniform(-1, 1, 2); h = 1e-5
        dx = (P(x + h, y) - P(x - h, y))/(2*h); dy = (P(x, y + h) - P(x, y - h))/(2*h); Pi = P(x, y)
        g = lambda A, B: 0.5*np.trace(A@B).real
        Om = lambda A, B: (1j*np.trace(Pi@(A@B - B@A))).real
        gxx = g(dx, dx)
        # a unit vector V at angle theta from d_x's normal, built in an orthonormal frame of g
        e1 = dx/np.sqrt(gxx); t = dy - g(dy, e1)*e1; e2 = t/np.sqrt(g(t, t))
        th = rng.uniform(-1.5, 1.5); V = np.sin(th)*e1 + np.cos(th)*e2
        C = gxx*np.sin(th)**2; det = g(V, V)*gxx - g(V, dx)**2
        worst_id = max(worst_id, abs(det - (gxx - C)))
        b = abs(Om(V, dx)) - 2*np.sqrt(max(det, 0)); worst_b = max(worst_b, b); slack.append(-b/(2*np.sqrt(det)))
    return worst_id, worst_b, np.median(slack)
i2, b2, s2 = guard(2, 300); i3, b3, s3 = guard(3, 300)
verdict['P2'] = max(i2, i3) < 1e-8 and max(b2, b3) < 1e-6
print(f"P2 sine guard: Gram identity max dev {max(i2, i3):.1e}; two-level max(|Omega| - bound) = {b2:.1e} (median slack {s2:.1e}: equality), "
      f"three-level {b3:.1e} (median slack {s3:.2f}: strict)  {'PASS' if verdict['P2'] else 'KILL'}")

# ---------- P3 fermion offset as a Smith coordinate ----------
# Claim: tau = log((1-C)/C) is a log-impedance, z = e^tau, and Gamma = (z-1)/(z+1) = 1 - 2C = tanh(tau/2) = tanh(omega/2T).
# Way 1 symbolic. Way 2 from the pair state by partial trace (Dirac pair cos r|00> + e^{i phi} sin r|11>).
# Opposite: C -> 1 - C (particle-hole) sends Gamma -> -Gamma; Gamma = 0 only at C = 1/2, i.e. omega = 0.
# Kill: way 2 off by > 1e-12 for any phase.
w_, T_ = sp.symbols('omega T', positive=True)
Cs = 1/(sp.exp(w_/T_) + 1); z = (1 - Cs)/Cs
ok_sym = sp.simplify(((z - 1)/(z + 1) - sp.tanh(w_/(2*T_))).rewrite(sp.exp)) == 0 and sp.simplify((1 - 2*Cs - sp.tanh(w_/(2*T_))).rewrite(sp.exp)) == 0
dev3 = 0.0
for wt in [0.1, 0.5, 1, 2, 5, 9]:
    r = np.arctan(np.exp(-wt/2))
    for phi in np.linspace(0, 2*np.pi, 7):
        psi = np.zeros(4, complex); psi[0] = np.cos(r); psi[3] = np.exp(1j*phi)*np.sin(r)
        rho = np.outer(psi, psi.conj()).reshape(2, 2, 2, 2); rho_out = np.einsum('ajbj->ab', rho)
        Cocc = rho_out[1, 1].real
        dev3 = max(dev3, abs((1 - 2*Cocc) - np.tanh(wt/2)))
verdict['P3'] = ok_sym and dev3 < 1e-12
print(f"P3 fermion offset is a Smith coordinate: symbolic {'yes' if ok_sym else 'NO'}; partial trace over 6 frequencies x 7 phases max dev {dev3:.1e}; "
      f"Gamma = 0 only at omega = 0  {'PASS' if verdict['P3'] else 'KILL'}")

# ---------- P4 the compression D(X) = P X P ----------
# Claim: D is idempotent for every projector P; it is an orthogonal projection (self-adjoint in the
# Frobenius inner product) exactly when P is Hermitian. Way 1 Hermitian P; way 2 oblique Riesz projector
# of a non-normal matrix (the opposite). Kill: idempotence fails, or self-adjointness holds for oblique P.
def Dmat(P):
    n = P.shape[0]; cols = []
    for k in range(n*n):
        E = np.zeros(n*n, complex); E[k] = 1; X = E.reshape(n, n); cols.append((P@X@P).ravel())
    return np.array(cols).T
n = 4
Q, _ = np.linalg.qr(rng.normal(size=(n, n)) + 1j*rng.normal(size=(n, n))); Ph = Q[:, :2]@Q[:, :2].conj().T
A = np.triu(rng.normal(size=(n, n)), 1)*3 + np.diag([1, 1.2, 5, 6.0])
ev, R = np.linalg.eig(A); L = np.linalg.inv(R); idx = np.argsort(ev.real)[:2]
Pr = sum(np.outer(R[:, i], L[i, :]) for i in idx)
res4 = []
for P in [Ph, Pr]:
    Dm = Dmat(P); res4.append((np.abs(Dm@Dm - Dm).max(), np.abs(Dm - Dm.conj().T).max(), np.linalg.norm(P, 2)))
verdict['P4'] = res4[0][0] < 1e-12 and res4[1][0] < 1e-10 and res4[0][1] < 1e-12 and res4[1][1] > 1e-3
print(f"P4 compression D(X) = PXP: Hermitian P idempotent {res4[0][0]:.1e}, self-adjoint {res4[0][1]:.1e}; "
      f"oblique Riesz P (norm {res4[1][2]:.2f}) idempotent {res4[1][0]:.1e}, NOT self-adjoint ({res4[1][1]:.2f})  {'PASS' if verdict['P4'] else 'KILL'}")

# ---------- P5 tip and tail: vectors swing, the cluster subspace holds ----------
# Claim (Davis-Kahan): sin of the largest principal angle between cluster subspaces <= ||E|| / delta,
# delta the gap from the cluster to the rest; single eigenvectors inside the cluster can turn by O(1).
# Way: 2000 random symmetric 6x6 with a tight 3-cluster; opposite: a wide internal split (vectors hold too).
# Kill: any ratio > 1, or the median single-vector turn not exceeding 10x the subspace angle in the tight case.
def trial(spread):
    Qm, _ = np.linalg.qr(rng.normal(size=(6, 6)))
    lam = np.r_[1 + spread*rng.uniform(-1, 1, 3), 3 + rng.uniform(0, 1, 3)]
    A0 = Qm@np.diag(lam)@Qm.T
    Ep = rng.normal(size=(6, 6)); Ep = (Ep + Ep.T)/2; Ep *= 0.02/np.linalg.norm(Ep, 2)
    w0, V0 = np.linalg.eigh(A0); w1, V1 = np.linalg.eigh(A0 + Ep)
    S0, S1 = V0[:, :3], V1[:, :3]
    sv = np.linalg.svd(S0.T@S1, compute_uv=False); sin_th = np.sqrt(max(0, 1 - sv.min()**2))
    delta = w0[3] - w0[2]
    turn = np.degrees(np.arccos(np.clip(np.abs(np.sum(V0[:, :3]*V1[:, :3], 0)), 0, 1))).max()
    return sin_th/(0.02/delta), turn, np.degrees(np.arcsin(sin_th))
tight = np.array([trial(1e-3) for _ in range(2000)]); wide = np.array([trial(0.4) for _ in range(2000)])
verdict['P5'] = tight[:, 0].max() <= 1 and wide[:, 0].max() <= 1 and np.median(tight[:, 1]) > 10*np.median(tight[:, 2])
print(f"P5 tip and tail, 2000 trials each: max DK ratio {max(tight[:,0].max(), wide[:,0].max()):.3f} (must be <= 1); "
      f"tight cluster: median vector turn {np.median(tight[:,1]):.2f} deg vs subspace {np.median(tight[:,2]):.3f} deg; "
      f"wide split (opposite): vector turn {np.median(wide[:,1]):.3f} deg  {'PASS' if verdict['P5'] else 'KILL'}")

print("\nSUMMARY:", ", ".join(f"{k} {'pass' if v else 'KILL'}" for k, v in verdict.items()))