#!/usr/bin/env python3
"""
test_349 -- EMBEDDED NON-CONVEX REALIZATIONS: SO(3)-SYMMETRIC LEAVES, THE RADIAL GRAPHS r = 1 + eps cos 2chi, AND A
            REALIZABLE PATH THROUGH THE DEGENERATE LOCUS

For MB §6.1 (the companion of MA). Conventions as in MA: K_ab = -N . d_a d_b x (the unit sphere with outward normal has
K = g); M = -G = adj K. chi is the polar angle on S^3, omega in S^2.

PRE-REGISTERED:
  P1 (SO(3)-symmetric hypersurfaces) For X = (rho(s) omega, z(s)) with |rho'| < 1 and z' = sqrt(1 - rho'^2) (s arclength
     of the profile), with the normal N = (sqrt(1 - rho'^2) omega, -rho'), the principal curvatures are
     k_mer = -rho''/sqrt(1 - rho'^2) and k_sph = sqrt(1 - rho'^2)/rho (twice), so adj K = diag((1 - rho'^2)/rho^2,
     -rho''/rho, -rho''/rho), the intrinsic M of test_347 P1, and sign k_mer = -sign rho''. (Symbolic, rho arbitrary.)
  P2 (radial graphs, intrinsic against extrinsic) For X = r(chi) n with r an arbitrary function, the metric is
     A^2 dchi^2 + Phi^2 (dtheta^2 + sin^2 theta dphi^2) with A = sqrt(r^2 + r'^2), Phi = r sin chi; the intrinsic chain
     gives M = diag(Q, P, P) with Q = k_sph^2, P = k_mer k_sph, for test_345's
     k_mer = (r^2 + 2r'^2 - r r'')/A^3, k_sph = (r - r' cot chi)/(r A). (Symbolic.)
  P3 (the family r = 1 + eps cos 2chi, exactly) With z = cos^2 chi:
     (a) r >= 1 - eps > 0, so X is an embedding for 0 <= eps < 1;
     (b) the numerator of k_sph, r - r' cot chi, is 1 - eps + 6 eps z > 0;
     (c) the numerator of k_mer, r^2 + 2r'^2 - r r'', is N(z) = 1 - 6eps + 5eps^2 + 12eps(1 + eps)z - 12eps^2 z^2, with
         N(0) = (1 - eps)(1 - 5eps), N(1) = (1 + eps)(1 + 5eps), and N = (12/25) z (6 - z) at eps = 1/5;
     (d) hence for 1/5 < eps < 1, N has exactly one root z* in (0, 1): K is indefinite exactly on the equatorial band
         cos^2 chi < z*, degenerate on the two latitude spheres cos^2 chi = z*, and definite on the polar caps, which
         contain the farthest points (the poles, r = 1 + eps). At eps = 1/5, K is semidefinite and degenerate only on
         the equator. Checked symbolically, and at eps = 3/10, 1/2, 9/10 numerically (z* and the signs).
  P4 (the band is a neck) Along the family, with Phi as a function of the profile's arclength s, the band
     {k_mer < 0} is exactly {Phi_ss > 0} (P1 with rho = Phi): at eps = 3/10 on a grid of 400 points, sign k_mer =
     -sign Phi_ss wherever |Phi_ss| > 1e-20.
  P5 (a realizable path through the degenerate locus) For eps = 3/10 the path r_t = 1 + t eps cos 2chi, t in [0, 1],
     is embedded throughout; K is definite for t < 2/3, degenerate exactly on the equator at t = 2/3, indefinite on a
     band for t > 2/3. Along it J = 2 pi^2 for t <= 2/3 and J > 2 pi^2 beyond (at t = 1/2, 2/3, 5/6, 1; 40 digits).

FIRST RUN: P1-P4 passed; P5 crashed before its check, in mpmath's bracketing root-finder. The fold point is now found
by plain bisection on the same bracket. No claim changed.
"""
import sys
import sympy as sp
import mpmath as mp

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)


R3 = range(3)

# ------------------------------------------------------------------------------------------------------------- P1
print("P1  SO(3)-symmetric hypersurfaces")
s, th, ph = sp.symbols("s theta phi", real=True)
rho = sp.Function("rho")(s)
r1, r2 = sp.diff(rho, s), sp.diff(rho, s, 2)
w = sp.sqrt(1 - r1**2)
om = sp.Matrix([sp.cos(th), sp.sin(th) * sp.cos(ph), sp.sin(th) * sp.sin(ph)])
zf = sp.Function("z")(s)
Xv = sp.Matrix([*(rho * om), zf])
Nv = sp.Matrix([*(w * om), -r1])
Y = [s, th, ph]
rule = {sp.diff(zf, s): w, sp.diff(zf, s, 2): sp.diff(w, s)}
tang = [sp.diff(Xv, v).subs(rule) for v in Y]
g = sp.Matrix(3, 3, lambda a, b: sp.simplify(tang[a].dot(tang[b])))
Kx = sp.Matrix(3, 3, lambda a, b: sp.simplify(-Nv.dot(sp.diff(Xv, Y[a], Y[b]).subs(rule))))
S_ = (g.inv() * Kx).applyfunc(sp.simplify)
ok1 = g == sp.diag(1, rho**2, rho**2 * sp.sin(th)**2) and all(sp.simplify(Nv.dot(T)) == 0 for T in tang)
ok1 &= sp.simplify(Nv.dot(Nv) - 1) == 0
ok1 &= (S_ - sp.diag(-r2 / w, w / rho, w / rho)).applyfunc(sp.simplify) == sp.zeros(3, 3)
ok1 &= (S_.adjugate() - sp.diag((1 - r1**2) / rho**2, -r2 / rho, -r2 / rho)).applyfunc(sp.simplify) == sp.zeros(3, 3)
check("k_mer = -rho''/sqrt(1-rho'^2), k_sph = sqrt(1-rho'^2)/rho; adj K = diag((1-rho'^2)/rho^2, -rho''/rho, -rho''/rho)",
      ok1)

# ------------------------------------------------------------------------------------------------------------- P2
print("P2  radial graphs: the intrinsic M against the embedding")
ch = sp.Symbol("chi", real=True)
rf = sp.Function("r")(ch)
rp, rpp = sp.diff(rf, ch), sp.diff(rf, ch, 2)
A = sp.sqrt(rf**2 + rp**2)
Phi = rf * sp.sin(ch)
gR = sp.diag(A**2, Phi**2, Phi**2 * sp.sin(th)**2)
X3 = [ch, th, ph]
gi = gR.inv()
Gam = [[[sp.simplify(sum(gi[c, d] * (sp.diff(gR[d, b], X3[a]) + sp.diff(gR[d, a], X3[b]) - sp.diff(gR[a, b], X3[d]))
                         for d in R3) / 2) for b in R3] for a in R3] for c in R3]


def Riem(c, d, a, b):
    return (sp.diff(Gam[c][b][d], X3[a]) - sp.diff(Gam[c][a][d], X3[b])
            + sum(Gam[c][a][e] * Gam[e][b][d] - Gam[c][b][e] * Gam[e][a][d] for e in R3))


Ric = sp.Matrix(3, 3, lambda d, b: sp.simplify(sum(Riem(c, d, c, b) for c in R3)))
Rs = sp.simplify(sum(gi[i, j] * Ric[i, j] for i in R3 for j in R3))
M = (Rs / 2 * sp.eye(3) - gi * Ric).applyfunc(sp.simplify)
kmer = (rf**2 + 2 * rp**2 - rf * rpp) / A**3
ksph = (rf - rp * sp.cot(ch)) / (rf * A)
ok2 = all(M[i, j] == 0 for i in R3 for j in R3 if i != j)
ok2 &= sp.simplify(M[0, 0] - ksph**2) == 0
ok2 &= sp.simplify(M[1, 1] - kmer * ksph) == 0 and sp.simplify(M[2, 2] - kmer * ksph) == 0
check("intrinsic M = diag(k_sph^2, k_mer k_sph, k_mer k_sph) for every r(chi)", ok2)

# ------------------------------------------------------------------------------------------------------------- P3
print("P3  the family r = 1 + eps cos 2chi")
e, zz = sp.symbols("epsilon z", real=True)
r = 1 + e * sp.cos(2 * ch)
r_1, r_2 = sp.diff(r, ch), sp.diff(r, ch, 2)
to_z = lambda expr: sp.expand(sp.expand_trig(sp.expand(expr)).subs(sp.sin(ch)**2, 1 - sp.cos(ch)**2)
                              ).subs(sp.cos(ch), sp.sqrt(zz))
Nsph = sp.simplify(to_z(sp.expand(sp.expand_trig(r - r_1 * sp.cos(ch) / sp.sin(ch)))))
Nmer = sp.expand(to_z(r**2 + 2 * r_1**2 - r * r_2))
want = 1 - 6 * e + 5 * e**2 + 12 * e * (1 + e) * zz - 12 * e**2 * zz**2
ok3 = sp.simplify(Nsph - (1 - e + 6 * e * zz)) == 0 and sp.expand(Nmer - want) == 0
ok3 &= sp.factor(want.subs(zz, 0)) == sp.factor((1 - e) * (1 - 5 * e)) and sp.factor(want.subs(zz, 1)) == sp.factor((1 + e) * (1 + 5 * e))
ok3 &= sp.expand(want.subs(e, sp.Rational(1, 5)) - sp.Rational(12, 25) * zz * (6 - zz)) == 0
check("(a)-(c) embedded; k_sph numerator 1 - eps + 6eps z; k_mer numerator N(z) with N(0) = (1-eps)(1-5eps)", ok3,
      f"N = {sp.factor(Nmer)}")
mp.mp.dps = 40
ok3d = True
rows = []
for ev in (sp.Rational(3, 10), sp.Rational(1, 2), sp.Rational(9, 10)):
    Nz = want.subs(e, ev)
    roots = [rt for rt in sp.solve(Nz, zz) if rt.is_real and 0 < rt < 1]
    ok3d &= len(roots) == 1 and Nz.subs(zz, 0) < 0 < Nz.subs(zz, 1)
    zs = roots[0]
    ok3d &= sp.diff(Nz, zz, 2) < 0                                    # concave: one crossing on (0, 1)
    for zv in (zs / 2, (zs + 1) / 2):
        ok3d &= (Nz.subs(zz, zv) < 0) == (zv < zs)
    rows.append(f"eps={ev}: z* = {sp.N(zs, 8)}, band |chi - pi/2| < {sp.N(sp.pi / 2 - sp.acos(sp.sqrt(zs)), 6)}")
N15 = want.subs(e, sp.Rational(1, 5))
ok3d &= all(N15.subs(zz, sp.Rational(k, 10)) > 0 for k in range(1, 11)) and N15.subs(zz, 0) == 0
check("(d) exactly one root z* in (0,1): indefinite band, two degenerate spheres, definite caps; eps = 1/5 semidefinite",
      ok3d, "; ".join(rows))

# ------------------------------------------------------------------------------------------------------------- P4
print("P4  the band is a neck of the orbit radius")
ev = mp.mpf(3) / 10
rr = lambda c: 1 + ev * mp.cos(2 * c)
r1n = lambda c: -2 * ev * mp.sin(2 * c)
r2n = lambda c: -4 * ev * mp.cos(2 * c)
An = lambda c: mp.sqrt(rr(c)**2 + r1n(c)**2)
Phin = lambda c: rr(c) * mp.sin(c)


def Phi_ss(c):                                   # d^2 Phi / ds^2 with ds = A dchi
    dPhi = lambda u: mp.diff(Phin, u) / An(u)
    return mp.diff(dPhi, c) / An(c)


ok4 = True
cnt = 0
for k in range(1, 400):
    c = mp.pi * k / 400
    km = (rr(c)**2 + 2 * r1n(c)**2 - rr(c) * r2n(c)) / An(c)**3
    pss = Phi_ss(c)
    if abs(pss) > mp.mpf("1e-20"):
        ok4 &= mp.sign(km) == -mp.sign(pss)
        cnt += 1
check("sign k_mer = -sign Phi_ss at eps = 3/10 (the band is where the orbit radius is strictly convex)", ok4 and cnt > 390,
      f"{cnt} points")

# ------------------------------------------------------------------------------------------------------------- P5
print("P5  a realizable path from the round sphere through the degenerate locus")
eps0 = mp.mpf(3) / 10


def Jpath(t):
    ee = t * eps0
    rr = lambda c: 1 + ee * mp.cos(2 * c)
    r1 = lambda c: -2 * ee * mp.sin(2 * c)
    r2 = lambda c: -4 * ee * mp.cos(2 * c)
    beta_p = lambda c: (rr(c)**2 + 2 * r1(c)**2 - rr(c) * r2(c)) / (rr(c)**2 + r1(c)**2)
    beta = lambda c: c - mp.atan(r1(c) / rr(c))
    nodes = [0, mp.pi / 2, mp.pi]
    if ee > mp.mpf(1) / 5:
        lo, hi = mp.mpf("0.6"), mp.pi / 2 - mp.mpf("1e-12")       # beta_p(lo) > 0 > beta_p(hi): bisect
        for _ in range(140):
            mid = (lo + hi) / 2
            lo, hi = (mid, hi) if beta_p(mid) > 0 else (lo, mid)
        a0 = (lo + hi) / 2
        nodes = [0, a0, mp.pi / 2, mp.pi - a0, mp.pi]
    J = 4 * mp.pi * mp.quad(lambda c: mp.sin(beta(c))**2 * abs(beta_p(c)), nodes)
    eq = beta_p(mp.pi / 2)                                       # sign of k_mer at the equator
    return J, eq, 1 - ee


ok5 = True
rows = []
for t in (mp.mpf(1) / 2, mp.mpf(2) / 3, mp.mpf(5) / 6, mp.mpf(1)):
    J, eq, rmin = Jpath(t)
    ok5 &= rmin > 0
    if t < mp.mpf(2) / 3:
        ok5 &= eq > 0 and abs(J - 2 * mp.pi**2) < mp.mpf("1e-25")
    elif t == mp.mpf(2) / 3:
        ok5 &= abs(eq) < mp.mpf("1e-30") and abs(J - 2 * mp.pi**2) < mp.mpf("1e-25")
    else:
        ok5 &= eq < 0 and J - 2 * mp.pi**2 > mp.mpf("1e-3")
    rows.append(f"t={mp.nstr(t, 3)}: J/2pi^2 = {mp.nstr(J / (2 * mp.pi**2), 10)}, beta'(equator) = {mp.nstr(eq, 3)}")
check("embedded throughout; definite, then degenerate at t = 2/3, then indefinite; J = 2 pi^2 until 2/3, then above",
      ok5, "; ".join(rows))

print()
if FAILURES:
    print(f"test_349: {len(FAILURES)} FAILED")
    sys.exit(1)
print("test_349: ALL PASS")