# SCRIPT: V10_CORRECTIONS_CHECK
# Each correction to v8.5.1 is checked before it is written. Kill line per item: the correction is
# NOT applied if its check fails. Each item states why the check decides it.
import numpy as np, sympy as sp
rng = np.random.default_rng(10)
ok = {}

# C1 Theorem 4 needs T and W in one frame. Why: g_W(D_T H, K) is symmetric only if (t_i + t_j) multiplies
# H_ij in the eigenframe of W; a T that is not diagonal there breaks it. Kill: the counterexample fails
# AND the diagonal case fails.
def gW(W, H, K):
    w, V = np.linalg.eigh(W); Hp, Kp = V.T@H@V, V.T@K@V
    return 0.5*np.sum(Hp*Kp/(w[:, None] + w[None, :]))
W = np.diag([0.5, 0.3, 0.2])
Tn = np.array([[1, 0.2, 0], [0.2, 1, 0], [0, 0, 1.0]])        # not diagonal in W's frame
A = lambda T, H: T@H + H@T - 2*np.trace(T@H)*W
H = np.diag([1.0, -1, 0]); K = np.array([[0, 1, 0], [1, 0, 0], [0, 0, 0.0]])
asym = gW(W, A(Tn, H), K) - gW(W, H, A(Tn, K))
Td = np.diag([2, 1, 1.3]); worst = 0.0
for _ in range(200):
    X = rng.normal(size=(3, 3)); H1 = X + X.T; H1 -= np.trace(H1)/3*np.eye(3)
    Y = rng.normal(size=(3, 3)); K1 = Y + Y.T; K1 -= np.trace(K1)/3*np.eye(3)
    worst = max(worst, abs(gW(W, A(Td, H1), K1) - gW(W, H1, A(Td, K1))))
ok['C1'] = abs(asym) > 1e-3 and worst < 1e-12
print(f"C1 non-commuting T: g(AH,K) - g(H,AK) = {asym:+.6f} (magnitude 2/15 = {2/15:.6f}); diagonal T, 200 random pairs: {worst:.1e}  {'APPLY' if ok['C1'] else 'DO NOT APPLY'}")
# the normalised scaling flow G -> e^{sT} G e^{sT}/Tr has Jacobian A(H) - 2 Tr(TG) H: a shift by a multiple of I
Gs = sp.Matrix(3, 3, lambda i, j: sp.Symbol(f'g{min(i,j)}{max(i,j)}')); Ts = sp.diag(*sp.symbols('t1:4')); s = sp.Symbol('s')
Hs = sp.Matrix(3, 3, lambda i, j: sp.Symbol(f'h{min(i,j)}{max(i,j)}')); e = sp.Symbol('e')
F = lambda G: (G + s*(Ts*G + G*Ts))/(1 + s*sp.trace(Ts*G + G*Ts))              # first order in s
dF = sp.diff(F(Gs + e*Hs), e).subs(e, 0); dF = sp.diff(dF, s).subs(s, 0)
claim = Ts*Hs + Hs*Ts - 2*(Ts*Hs).trace()*Gs - 2*(Ts*Gs).trace()*Hs
ok['C1b'] = sp.simplify(dF - claim) == sp.zeros(3, 3)
print(f"C1b Jacobian of the normalised flow = A(H) - 2Tr(TG)H symbolically: {'YES' if ok['C1b'] else 'NO'}")

# C2 sigma coordinates: the eigenvalue block is the round octant (curvature 1), not flat.
# Why: Gaussian curvature of sum dw_i^2/(4w_i) on sum w_i = 1, computed from the metric alone (Brioschi).
u, v = sp.symbols('u v', positive=True)             # w1 = u, w2 = v, w3 = 1 - u - v
w3 = 1 - u - v
E_ = 1/(4*u) + 1/(4*w3); F_ = 1/(4*w3); G_ = 1/(4*v) + 1/(4*w3)
def brioschi(E, F, G):
    Eu, Ev, Fu, Fv, Gu, Gv = [sp.diff(f, x) for f in (E, F, G) for x in (u, v)]
    M1 = sp.Matrix([[-sp.diff(E, v, 2)/2 + sp.diff(F, u, v) - sp.diff(G, u, 2)/2, Eu/2, Fu - Ev/2],
                    [Fv - Gu/2, E, F], [Gv/2, F, G]])
    M2 = sp.Matrix([[0, Ev/2, Gu/2], [Ev/2, E, F], [Gu/2, F, G]])
    return (M1.det() - M2.det())/(E*G - F**2)**2
Kc = [float(brioschi(E_, F_, G_).subs({u: a, v: b})) for a, b in [(0.2, 0.3), (0.5, 0.1), (0.05, 0.05)]]
ok['C2'] = all(abs(k - 1) < 1e-9 for k in Kc)
print(f"C2 Gaussian curvature of the eigenvalue block at three points: {Kc}  {'APPLY (round, not flat)' if ok['C2'] else 'DO NOT APPLY'}")

# C3 real qubit: Bures metric = hemisphere of radius 1/2 (curvature 4), boundary the circle of pure states.
# Why: Bures metric from the SLD on rho = (I + x sx + z sz)/2, compared with (1/4)(dchi^2 + sin^2 chi dphi^2), r = sin chi.
sx = np.array([[0, 1], [1, 0.0]]); sz = np.array([[1, 0], [0, -1.0]])
def bures(rho, d):
    l, V = np.linalg.eigh(rho); dp = V.T@d@V
    return 0.5*np.sum(dp**2/(l[:, None] + l[None, :]))   # (1/4) * SLD Fisher = (1/2) sum |d_ij|^2/(l_i+l_j)
dev = 0.0
for _ in range(100):
    chi, ph = rng.uniform(0.05, 1.5), rng.uniform(0, 2*np.pi); r = np.sin(chi)
    rho = 0.5*(np.eye(2) + r*np.cos(ph)*sx + r*np.sin(ph)*sz)
    dchi = 0.5*np.cos(chi)*(np.cos(ph)*sx + np.sin(ph)*sz); dph = 0.5*r*(-np.sin(ph)*sx + np.cos(ph)*sz)
    dev = max(dev, abs(bures(rho, dchi) - 0.25), abs(bures(rho, dph) - 0.25*np.sin(chi)**2))
ok['C3'] = dev < 1e-12
print(f"C3 real-qubit Bures = (1/4)(dchi^2 + sin^2 chi dphi^2), chi in [0, pi/2): max dev {dev:.1e}; boundary circle length pi  {'APPLY' if ok['C3'] else 'DO NOT APPLY'}")

# C4 hidden block of the transport has rank at most one, zero when T is a multiple of I.
def spine_rank(t):
    w = np.array([0.5, 0.3, 0.2]); Wm = np.diag(w); T = np.diag(t)
    basis = []
    for i in range(3):
        for j in range(i, 3):
            M = np.zeros((3, 3)); M[i, j] = M[j, i] = 1; basis.append(M)
    PW = lambda Z: Z - np.trace(Z)*Wm
    DT = lambda Z: T@Z + Z@T
    QW = lambda Z: np.trace(Z)*Wm
    cols = [ (PW(DT(QW(DT(PW(b)))))).ravel() for b in basis]
    return np.linalg.matrix_rank(np.array(cols).T, tol=1e-10)
r_an, r_iso = spine_rank([2, 1, 1.3]), spine_rank([1.5, 1.5, 1.5])
ok['C4'] = r_an == 1 and r_iso == 0
print(f"C4 rank of P_W D_T Q_W D_T P_W: anisotropic {r_an}, isotropic {r_iso}  {'APPLY (rank at most one)' if ok['C4'] else 'DO NOT APPLY'}")

# C5 the geometric-phase frame-change claim. Why: a single-valued invertible F(theta) = rotation, constant
# biorthogonal pair u, v; the transported phase picks up oint u.F^{-1}dF v/(u.v), which is not zero.
uu = np.array([1.0, 0]); vv = np.array([1.0, 1])
Jr = np.array([[0, -1], [1, 0.0]])                  # F^{-1} dF = J dtheta for F = rotation
extra = 2*np.pi*(uu@Jr@vv)/(uu@vv)
ok['C5'] = abs(extra) > 1e-6
print(f"C5 rotation frame change adds oint = {extra:+.6f} (exp = {np.exp(-extra):.4g}, not 1)  {'APPLY (claim false in general)' if ok['C5'] else 'DO NOT APPLY'}")
# diagonal F: the extra term is sum_i (u_i v_i/u.v) dlog F_i, the share-weighted log-derivative (the residue mechanism)
th = sp.Symbol('th'); f1, f2 = sp.exp(sp.sin(th)), 2 + sp.cos(th)
U2, V2 = sp.Matrix([1, 2]), sp.Matrix([3, 1])
Fm = sp.diag(f1, f2); ex = (U2.T*Fm.inv()*sp.diff(Fm, th)*V2)[0]/(U2.T*V2)[0]
ok['C5b'] = sp.simplify(ex - (3*sp.diff(sp.log(f1), th) + 2*sp.diff(sp.log(f2), th))/5) == 0
print(f"C5b diagonal F: extra term = share-weighted dlog F_i: {'YES' if ok['C5b'] else 'NO'}")

# C6 exact projector norm for a real 2x2 with real distinct eigenvalues: ||P||_2^2 = 1 + (b - c)^2/g^2.
dev6 = 0.0; n6 = 0
while n6 < 500:
    a, b, c, d = rng.normal(size=4)
    disc = (a - d)**2 + 4*b*c
    if disc <= 1e-6: continue
    M = np.array([[a, b], [c, d]]); g = np.sqrt(disc); lm = (a + d - g)/2
    P = (M - lm*np.eye(2))/g
    dev6 = max(dev6, abs(np.linalg.norm(P, 2)**2 - (1 + (b - c)**2/g**2))); n6 += 1
ok['C6'] = dev6 < 1e-9
print(f"C6 ||P+||^2 = 1 + (b-c)^2/g^2 on 500 random real matrices: max dev {dev6:.1e}  {'APPLY' if ok['C6'] else 'DO NOT APPLY'}")

# C7 F_Q = 8K only for motions that keep full rank. Why: on G(s) = diag((1-s^2)/2, (1-s^2)/2, s^2) the SLD
# Fisher information on the support is computed exactly at s = 0 and beside it; 8K from the sigma metric.
def sld_fisher(G, dG):
    l, V = np.linalg.eigh(G); dp = V.T@dG@V; tot = 0.0
    for i in range(3):
        for j in range(3):
            if l[i] + l[j] > 1e-15: tot += 2*dp[i, j]**2/(l[i] + l[j])
    return tot
rows7 = []
for s_ in [0.1, 0.01, 0.001, 0.0]:
    G = np.diag([(1 - s_**2)/2, (1 - s_**2)/2, s_**2]); dG = np.diag([-s_, -s_, 2*s_])
    K = 0.5*(2*s_**2/(2*(1 - s_**2)) + 1.0)         # (1/2) sum dw^2/(4w); the w3 term is (d sigma3)^2 = 1
    rows7.append((s_, sld_fisher(G, dG), 8*K))
ok['C7'] = abs(rows7[-1][1]) < 1e-12 and abs(rows7[-1][2] - 4) < 1e-12 and abs(rows7[0][1] - rows7[0][2]) < 1e-9
print("C7 s, F_Q (SLD), 8K: " + "; ".join(f"{a:g}: {b:.6f}, {c:.6f}" for a, b, c in rows7) + f"  {'APPLY (full rank only)' if ok['C7'] else 'DO NOT APPLY'}")

# C8 finite-amplitude claim for the square. Why: the exact potential on the lift path, U(s) with s = sigma3.
U = lambda s: -4*np.sqrt(2)/np.sqrt(1 + s**2) - 2/np.sqrt(1 - s**2)
c2 = 2*np.sqrt(2) - 1; S = sp.Symbol('S')
quart = sp.series(-4*sp.sqrt(2)/sp.sqrt(1 + S**2) - 2/sp.sqrt(1 - S**2), S, 0, 6).removeO().coeff(S, 4)
dU, harm = U(0.3) - U(0), c2*0.09
from scipy.optimize import brentq
E = dU; s_exact = brentq(lambda s: U(s) - U(0) - E, 1e-9, 0.5); s_harm = np.sqrt(E/c2)
ok['C8'] = float(quart) < 0 and s_exact > s_harm
print(f"C8 quartic coefficient {sp.nsimplify(quart)} = {float(quart):.4f}; at sigma3 = 0.3 exact rise {dU:.10f} vs harmonic {harm:.10f}; "
      f"energy {E:.4f} reaches sigma3 = {s_exact:.4f}, harmonic bound says {s_harm:.4f}  {'APPLY (at most is false)' if ok['C8'] else 'DO NOT APPLY'}")

# C9 local-symbol control at r = 2M: the symbol is the zero matrix, so no band projectors there; the limit exists.
k_, m_ = 0.7, 1.0
P_lim = 0.5*(np.eye(2) + (k_*sx + m_*sz)/np.hypot(k_, m_))
hs = 0.0*(k_*sx + m_*sz)
ok['C9'] = np.allclose(hs, 0) and np.allclose(P_lim@P_lim, P_lim)
print(f"C9 at N = 0 the symbol is zero (every vector an eigenvector); the N > 0 projector is constant, so its limit exists  {'APPLY' if ok['C9'] else 'DO NOT APPLY'}")

print("\nSUMMARY:", ", ".join(f"{k} {'apply' if v else 'HOLD'}" for k, v in ok.items()))