#!/usr/bin/env python3
"""
test_329 -- THE POLYNOMIAL (THOMAS) FORM OF REALIZABILITY AT THE DEGENERATE LOCUS: NOT VACUOUS WHERE rank M = 2,
            IDENTICALLY ZERO ON REALIZABLE rank M <= 1, AND THE RANK-ONE EINSTEIN TENSOR'S BIANCHI SPLIT

For MA section 8. E(g) = 2 (det M) Cod(P) - d(det M) ^ P with P = adj M, computed as exact rational jets at a point,
independent of the audit's polynomial-check script and of tests 246, 279 and 291.

PRE-REGISTERED:
  P1 Algebra. For 3x3 symmetric matrices: rank adj(K) is 3, 1 or 0 as rank K is 3, 2 or <= 1, so a rank-2 M is
         never adj K; and for rank-2 M with kernel spanned by unit v, adj M = (m1 m2) v v^T. (Symbolic, diagonal.)
  P2 Not vacuous where rank M = 2 (the audit's example, P4-C04/P4-X04, rebuilt here): for
         g = dr^2 + (1 + r + r^3)^2 dx^2 + (1 + r + r^2)^2 dy^2 at r = 0, M = diag(-1, -2, 0), det M = 0,
         d_r det M = -12, P = diag(0, 0, 2), and E_ryy = 24. E there is exactly -(d det M ^ P), i.e. the condition
         d(det M) ^ v = 0 with v = d_y, which fails.
  P3 Identically zero on a realizable rank-one point: the cone over a Clifford torus,
         g = dt^2 + t^2 (cos^2 a dth1^2 + sin^2 a dth2^2), has M = diag(-1/t^2, 0, 0), rank 1, det M = 0, P = 0, and
         E = 0 at the point (all components), while F (built on sqrt(det M) M^{-1}) is undefined there.
  P4 The Bianchi split on that cone: v = d_t is geodesic and div(mu v) = 0 with mu = -1/t^2.

BUNDLE EDIT (29 October): the archive's [M] labels were removed from the claims; nothing else changed.
"""
import sys
import sympy as sp

FAILURES = []


def check(name, ok, detail=""):
    print(f"  [{'PASS' if ok else 'FAIL'}] {name}" + (f"  -- {detail}" if detail else ""))
    if not ok:
        FAILURES.append(name)


X = sp.symbols("u0:3")


def thomas(g, pt):
    """E, det M, M, P at a point, exactly (symbolic metric, then evaluated)."""
    gi = sp.simplify(g.inv())
    Gam = [[[sp.simplify(sum(gi[a, e] * (sp.diff(g[e, b], X[c]) + sp.diff(g[e, c], X[b]) - sp.diff(g[b, c], X[e]))
                             for e in range(3)) / 2) for c in range(3)] for b in range(3)] for a in range(3)]
    Riem = lambda a, b, c, d: (sp.diff(Gam[a][d][b], X[c]) - sp.diff(Gam[a][c][b], X[d])
                               + sum(Gam[a][c][e] * Gam[e][d][b] - Gam[a][d][e] * Gam[e][c][b] for e in range(3)))
    Ric = sp.Matrix(3, 3, lambda b, d: sp.simplify(sum(Riem(a, b, a, d) for a in range(3))))
    R = sp.simplify(sum(gi[b, d] * Ric[b, d] for b in range(3) for d in range(3)))
    G = Ric - R * g / 2
    M = sp.simplify(-gi * G)
    detM = sp.simplify(M.det())
    Pl = sp.simplify(g * M.adjugate())
    nab = lambda a, b, c: sp.diff(Pl[b, c], X[a]) - sum(Gam[e][a][b] * Pl[e, c] + Gam[e][a][c] * Pl[b, e] for e in range(3))
    sub = dict(zip(X, pt))
    E = {(a, b, c): sp.simplify((2 * detM * (nab(a, b, c) - nab(b, a, c))
                                 - (sp.diff(detM, X[a]) * Pl[b, c] - sp.diff(detM, X[b]) * Pl[a, c])).subs(sub))
         for (a, b) in ((0, 1), (0, 2), (1, 2)) for c in range(3)}
    return E, detM, M, Pl, Gam, sub


# ---------------------------------------------------------------- P1
print("P1  the adjugate's ranks; a rank-2 M is never adj K")
k1, k2, k3 = sp.symbols("k1 k2 k3", nonzero=True)
adj = lambda d: sp.diag(*d).adjugate()
ok1 = (adj((k1, k2, k3)).rank() == 3 and adj((k1, k2, 0)).rank() == 1 and adj((k1, 0, 0)).rank() == 0
       and adj((k1, k2, 0)) == sp.diag(0, 0, k1 * k2))
check("rank adj K = 3, 1, 0 for rank K = 3, 2, <= 1; adj of rank-2 diag(m1, m2, 0) is m1 m2 e3 e3^T", ok1)

# ---------------------------------------------------------------- P2
print("P2  not vacuous where rank M = 2 (the audit's example)")
r, xx, yy = X
g2 = sp.diag(1, (1 + r + r**3)**2, (1 + r + r**2)**2)
E2, detM2, M2, P2l, _, sub2 = thomas(g2, (0, 0, 0))
M2p, dM2 = M2.subs(sub2), sp.diff(detM2, r).subs(sub2)
nonzero = {k: v for k, v in E2.items() if v != 0}
ok2 = (M2p == sp.diag(-1, -2, 0) and detM2.subs(sub2) == 0 and dM2 == -12 and P2l.subs(sub2) == sp.diag(0, 0, 2)
       and nonzero == {(0, 2, 2): 24})
check("M = diag(-1,-2,0), det M = 0, d_r det M = -12, P = diag(0,0,2), E_ryy = 24 and nothing else", ok2,
      f"nonzero E components {nonzero}")

# ---------------------------------------------------------------- P3, P4
print("P3  identically zero on a realizable rank-one point: the cone over a Clifford torus")
t, th1, th2 = X
a = sp.Rational(1, 3)        # stand-in for the torus angle: cos^2 a, sin^2 a replaced by rationals summing to 1
c2, s2 = sp.Rational(16, 25), sp.Rational(9, 25)
g3 = sp.diag(1, t**2 * c2, t**2 * s2)
E3, detM3, M3, P3l, Gam3, sub3 = thomas(g3, (sp.Rational(3, 2), sp.Rational(1, 5), sp.Rational(2, 7)))
ok3 = (sp.simplify(M3 - sp.diag(-1 / t**2, 0, 0)) == sp.zeros(3, 3) and sp.simplify(detM3) == 0
       and sp.simplify(P3l) == sp.zeros(3, 3) and all(v == 0 for v in E3.values()))
check("M = diag(-1/t^2, 0, 0), det M = 0, P = 0 identically, E = 0 at the point", ok3)
print("P4  the Bianchi split on the cone")
mu = -1 / t**2
sqrtg = sp.sqrt(g3.det())
div = sp.simplify(sp.diff(sqrtg * mu, t) / sqrtg)
geod = all(sp.simplify(Gam3[i][0][0]) == 0 for i in range(3))
check("v = d_t geodesic (Gamma^i_tt = 0); div(mu v) = 0", geod and div == 0, f"div = {div}")

print()
print(f"test_329: {'ALL PASS' if not FAILURES else 'FAILURES: ' + ', '.join(FAILURES)}")
sys.exit(1 if FAILURES else 0)