#!/usr/bin/env python3
"""
test_342 -- THE RANK-ONE BRANCH: THE BIANCHI SPLIT AND THE CODAZZI BLOCK SYSTEM

Where rank K = 2, M = adj K = mu v v^T with mu = k1 k2 and v spanning ker K. In an orthonormal frame (e1, e2) of v-perp
put B_ab = <e_b, nabla_{e_a} v>; B is intrinsic, since v is an eigenvector of the Einstein tensor.

PRE-REGISTERED:
  Q1 adj diag(k1, k2, 0) = diag(0, 0, k1 k2), of rank one; adj M = 0 there; the nonzero eigenvector of M spans ker K.
  Q2 If -G = mu v v^T with |v| = 1 (at a point, v = e3, d_a v_3 = 0), the contracted Bianchi identity splits: its
     component along v is div(mu v), and its components across v are mu times the acceleration nabla_v v.
  P1 Codazzi contracted with v, using K v = 0, gives (nabla_a K)(e_b, v) = -(B K)_ab, so B K = K B^T on v-perp:
     b21 k11 + k12 (b22 - b11) - b12 k22 = 0; in a principal frame k2 b12 = k1 b21.
  P2 Equal principal curvatures force b12 = b21 (no twist). A pure rotation of the rulings (b11 = b22,
     b21 = -b12 != 0) forces the block to be trace-free, so mu < 0 there.
  P3 Pointwise, {det K_perp = mu, B K = K B^T} has rank 2 on the three entries of K_perp: one free parameter.
  P4 The transport nabla_v K = -B K, with Jacobi's formula, gives v(ln mu) = -tr B: div(mu v) = 0 again.
  P5 The cone t T over a Clifford torus T in S^3: K kills the ruling, B = I/t, the block is diag(tan al, -cot al)/t,
     the frame is parallel along the ruling and d/dt of the block is -B K exactly; P1's equation is empty there.
  P6 Seams: where det M > 0, sqrt(det M) M^{-1} = sign(det K) K (20 random symmetric K). Across a seam, where one
     principal curvature passes through zero, K~ flips sign while sign(det K) K~ = K stays continuous.

This script replaces test_279 (its Q1-Q2 are that script's P2-P3) and test_291 (P1-P6) for the paper's
verification bundle; the computations are theirs, the descriptions are restated.
"""
import sys
import random
import sympy as sp

FAILURES = []


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


k1, k2, k3 = sp.symbols('k1 k2 k3', nonzero=True)
# ---------------------------------------------------------------- Q1 (rank one)
K2 = sp.diag(k1, k2, 0)
M2 = K2.adjugate()
check("Q1 adj diag(k1, k2, 0) = diag(0, 0, k1 k2): rank one", M2 == sp.diag(0, 0, k1 * k2) and M2.rank() == 1)
check("Q1 adj M = 0 on the branch (adj M vanishes)", M2.adjugate() == sp.zeros(3, 3))
check("Q1 M's nonzero eigenvector spans ker K (the ruling)",
      K2 * sp.Matrix([0, 0, 1]) == sp.zeros(3, 1) and M2 * sp.Matrix([0, 0, 1]) == k1 * k2 * sp.Matrix([0, 0, 1]))

# ---------------------------------------------------------------- Q2 (the Bianchi split)
mu = sp.Symbol('mu', nonzero=True)
dmu = sp.symbols('dmu1:4')
dv = sp.Matrix(3, 3, lambda a, b: sp.Symbol(f'dv_{a}{b}'))   # dv[a,b] = d_a v_b at the point
dv[:, 2] = sp.zeros(3, 1)                                      # |v| = 1 with v = e3: d_a v_3 = 0
v = [0, 0, 1]
E = [sum(dmu[a] * v[a] * v[b] + mu * dv[a, a] * v[b] + mu * v[a] * dv[a, b] for a in range(3)) for b in range(3)]
divmuv = sum(dmu[a] * v[a] + mu * dv[a, a] for a in range(3))
acc = [sum(v[a] * dv[a, b] for a in range(3)) for b in range(3)]
check("Q2 along v: the Bianchi component is exactly div(mu v)", sp.simplify(E[2] - divmuv) == 0)
check("Q2 across v: the components are mu times the acceleration of the ruling",
      all(sp.simplify(E[b] - mu * acc[b]) == 0 for b in (0, 1)))


# ---------------------------------------------------------------- P1
k11, k12, k22, b11, b12, b21, b22 = sp.symbols('k11 k12 k22 b11 b12 b21 b22', real=True)
K = sp.Matrix([[k11, k12], [k12, k22]])
B = sp.Matrix([[b11, b12], [b21, b22]])
nablaK_v = -(B * K)                      # (nabla_a K)(e_b, v) = -sum_c B_ac K_cb
cod = sp.expand(nablaK_v[0, 1] - nablaK_v[1, 0])
expected = b21 * k11 + k12 * (b22 - b11) - b12 * k22
check("P1 Codazzi with v: b21 k11 + k12 (b22 - b11) - b12 k22 = 0",
      sp.expand(cod - expected) == 0, f"{cod}")
check("P1 equivalently B K = K B^T (the diagonal is identically zero)",
      sp.expand((B * K - K * B.T)[0, 0]) == 0 and sp.expand((B * K - K * B.T)[1, 1]) == 0
      and sp.expand((B * K - K * B.T)[0, 1] + cod) == 0)
k1, k2 = sp.symbols('k1 k2', real=True)
princ = sp.expand(cod.subs({k11: k1, k22: k2, k12: 0}))
check("P1 principal frame: k2 b12 = k1 b21", sp.expand(princ - (k1 * b21 - k2 * b12)) == 0, f"{princ}")

# ---------------------------------------------------------------- P2
k = sp.Symbol('k', nonzero=True)
check("P2 k1 = k2 = k forces b12 = b21 (no twist)",
      sp.solve(princ.subs({k1: k, k2: k}), b12) == [b21])
f, c = sp.symbols('f c', real=True)
rot = sp.factor(cod.subs({b11: f, b22: f, b12: c, b21: -c}))
check("P2 a pure rotation of the rulings forces tr K = 0 on the block (so mu = det < 0)",
      sp.expand(rot + c * (k11 + k22)) == 0, f"{rot}")

# ---------------------------------------------------------------- P3
random.seed(291)
mu_v = sp.Rational(-3, 2)
Bv = {b11: sp.Rational(2, 7), b12: sp.Rational(5, 3), b21: sp.Rational(-1, 4), b22: sp.Rational(-3, 5)}
F = sp.Matrix([k11 * k22 - k12 ** 2 - mu_v, cod.subs(Bv)])
J = F.jacobian([k11, k12, k22])
sols = []
for s in (sp.Rational(1, 2), sp.Rational(3, 2), sp.Integer(-2)):
    k22s = sp.solve(cod.subs(Bv).subs(k12, s), k22)[0]
    for r in sp.solve(sp.Eq(k11 * k22s - s ** 2, mu_v), k11):
        if r.is_real:
            pt = {k11: r, k12: s, k22: k22s.subs(k11, r)}
            sols.append(pt)
ranks = [J.subs(pt).rank(simplify=True) for pt in sols]
resid = [max(abs(sp.N(e)) for e in F.subs(pt)) for pt in sols]
check("P3 the pointwise system has a curve of solutions with Jacobian rank 2 (one free parameter)",
      len(sols) >= 3 and all(rk == 2 for rk in ranks) and max(resid) < 1e-20,
      f"{len(sols)} solutions on distinct k12, ranks {ranks}")

# ---------------------------------------------------------------- P4
dK = -(B * K)                                   # nabla_v K in a frame parallel along v
dlogdet = sp.simplify((K.inv() * dK).trace())
check("P4 transport + Jacobi: v(ln det K) = -tr B = -div v, i.e. nabla(mu v) = 0",
      sp.simplify(dlogdet + B.trace()) == 0, f"{dlogdet}")

# ---------------------------------------------------------------- P5
t, u, w = sp.symbols('t u w', positive=True)
al = sp.pi / 5
T = sp.Matrix([sp.cos(al) * sp.cos(u), sp.cos(al) * sp.sin(u), sp.sin(al) * sp.cos(w), sp.sin(al) * sp.sin(w)])
N = sp.Matrix([-sp.sin(al) * sp.cos(u), -sp.sin(al) * sp.sin(u), sp.cos(al) * sp.cos(w), sp.cos(al) * sp.sin(w)])
X = t * T
q = [t, u, w]
Xd = [sp.diff(X, x) for x in q]
ok_normal = all(sp.simplify(N.dot(v_)) == 0 for v_ in Xd) and sp.simplify(N.dot(N) - 1) == 0
g = sp.Matrix(3, 3, lambda i, j: sp.simplify(Xd[i].dot(Xd[j])))
h = sp.Matrix(3, 3, lambda i, j: sp.simplify(sp.diff(X, q[i], q[j]).dot(N)))
S = sp.simplify(g.inv() * h)
check("P5 cone over a Clifford torus: N is the unit normal and K kills the ruling d_t",
      ok_normal and all(sp.simplify(S[i, 0]) == 0 for i in range(3)) and all(sp.simplify(h[0, j]) == 0 for j in range(3)))
e_u = Xd[1] / sp.sqrt(g[1, 1])
e_w = Xd[2] / sp.sqrt(g[2, 2])
blk = sp.Matrix([[h[1, 1] / g[1, 1], h[1, 2] / sp.sqrt(g[1, 1] * g[2, 2])],
                 [h[2, 1] / sp.sqrt(g[1, 1] * g[2, 2]), h[2, 2] / g[2, 2]]]).applyfunc(sp.simplify)
Bc = sp.Matrix(2, 2, lambda a, b_: sp.simplify([e_u, e_w][b_].dot(sp.diff(Xd[0], q[a + 1]) / sp.sqrt(g[a + 1, a + 1]))))
check("P5 the block is diag(tan al, -cot al)/t and B = I/t",
      max(abs(float((blk - sp.diag(sp.tan(al), -1 / sp.tan(al)) / t).subs({t: tv, u: 0.3, w: 1.1})[i]))
          for tv in (0.7, 1.9) for i in range(4)) < 1e-14
      and sp.simplify(Bc - sp.eye(2) / t) == sp.zeros(2, 2), f"B {list(Bc)}")
check("P5 frame parallel along the ruling, and d/dt block = -B K exactly",
      sp.simplify(sp.diff(e_u, t)) == sp.zeros(4, 1) and sp.simplify(sp.diff(e_w, t)) == sp.zeros(4, 1)
      and sp.simplify(sp.diff(blk, t) + Bc * blk) == sp.zeros(2, 2))
check("P5 cone class B = f I: P1's equation is empty there",
      sp.simplify(cod.subs({k11: blk[0, 0], k12: blk[0, 1], k22: blk[1, 1],
                            b11: Bc[0, 0], b12: Bc[0, 1], b21: Bc[1, 0], b22: Bc[1, 1]})) == 0)

# ---------------------------------------------------------------- P6
random.seed(6)
okP6 = True
for _ in range(20):
    A = sp.Matrix(3, 3, lambda i, j: sp.Rational(random.randint(-9, 9), random.randint(1, 5)))
    K3 = (A + A.T) / 2
    d = K3.det()
    if d == 0:
        continue
    M = K3.adjugate()
    Kt = sp.sqrt(M.det()) * M.inv()
    okP6 &= sp.simplify(Kt - sp.sign(d) * K3) == sp.zeros(3, 3)
check("P6 where det M > 0, sqrt(det M) M^-1 = sign(det K) K (20 random symmetric K)", okP6)
s_, a2, a3 = sp.symbols('s a2 a3', real=True)
Kseam = sp.diag(s_, 2, 3)
Mseam = Kseam.adjugate()
Ktil = sp.sqrt(Mseam.det()) * Mseam.inv()
lim_m = sp.limit(Ktil[1, 1], s_, 0, '-')
lim_p = sp.limit(Ktil[1, 1], s_, 0, '+')
rec = [sp.limit((sp.sign(Kseam.det()) * Ktil)[1, 1], s_, 0, d_) for d_ in ('-', '+')]
check("P6 across a seam K~ flips sign; sign(det K) K~ = K is continuous, so K is fixed there",
      lim_m == -lim_p and lim_p != 0 and rec[0] == rec[1] == 2, f"K~_22: {lim_m} | {lim_p}; K_22: {rec}")

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