#!/usr/bin/env python3
"""
test_345 -- CLOSED HYPERSURFACES OF R^4: THE FARTHEST-POINT BOUND, ITS INTRINSIC FORM, THE DEGENERATE GEOMETRIES,
            AND WHERE DEFINITENESS IS LOST

For MB §3 (the companion of MA). Conventions as in MA: K_ab = -N . d_a d_b x, so the unit sphere with outward normal has
K = g; M = -G = adj K.

PRE-REGISTERED:
  P1 (Lemma 3.1, exact) For the ellipsoid with semi-axes a1 >= a2, a3, a4 > 0, at the farthest point from the centre,
     x = a1 e1, the principal curvatures are a1/a_j^2 (j = 2, 3, 4), each >= 1/a1, so det K >= a1^{-3} and
     det M >= a1^{-6}; equality throughout iff a2 = a3 = a4 = a1 (the round sphere). Symbolic.
  P2 (Lemma 3.1, numerical, non-convex) For closed star-shaped hypersurfaces |x| = r(x/|x|) with r = 1 + (a random
     polynomial of degree <= 4 in the coordinates of S^3), scaled so that the hypersurface is visibly non-convex, the
     second fundamental form at the maximum of r has every eigenvalue >= 1/R (R = max r), to 1e-20, with the
     second fundamental form computed from the level-set Hessian at 40 digits. Three random examples; each has a point
     where K is indefinite (a negative eigenvalue somewhere on a sample), so the bound is not a convexity artefact.
  P3 (Corollary 3.3(2), the constant) The circumradius of a regular simplex in R^4 of edge d is d sqrt(2/5) (exact), the
     constant in Jung's theorem at n = 4. Hence max det M >= (5/2)^3 / diam^6 = 125/(8 diam^6) for a realizable closed
     metric; on the unit S^3 (det M = 1, diam = pi) the bound holds strictly.
  P4 (Corollary 3.3(1)) det M = 0 identically for the product metrics on S^2(r) x S^1, H^2(-1/r^2) x S^1 and the flat
     T^3 (symbolic); det M = 1 for the round metric of curvature 1 (any space form S^3/Gamma).
  P5 (Theorem 3.4's mechanism, and §6) For the radial graphs r = 1 + eps cos 2chi (rotationally symmetric about the
     x1-axis, chi the polar angle):
       (a) the formulas k_mer = (r^2 + 2r'^2 - r r'')/(r^2 + r'^2)^{3/2} and k_sph = (r - r' cot chi)/(r sqrt(r^2 + r'^2))
           agree with the level-set computation at three points for eps = 0.3 (to 1e-25);
       (b) at the equator k_sph > 0 and r^2 + 2r'^2 - r r'' = (1 - eps)(1 - 5 eps) exactly;
       (c) for eps = 1/10 and 3/20, det K > 0 and K is definite at every grid point (chi = k pi/400); at eps = 1/5 the
           minimum of det K over the grid is 0, attained at the equator; at eps = 3/10, det K < 0 at the equator, and on
           the grid K stays definite at every point between the poles and the first zero of k_mer.

FIRST RUN, AND THE CHANGES (registered claims unchanged in substance):
  P2 failed on its own witness, not on the bound. At the farthest point every principal curvature was >= 1/R in all
  three examples (min kappa 0.9505 vs 1/R 0.9103; 1.1869 vs 0.9303; 1.0820 vs 0.9243), but two of the three random
  hypersurfaces were convex, so the registered non-convexity witness failed for them. Re-posed: each example is now
  r = 1 + (7/20)(2 u1^2 - 1) + (a random polynomial of degree <= 4 with coefficients in [-9, 9]/200), which is
  non-convex by construction (the first part is P5's family at eps = 7/20 > 1/5) and has no symmetry; the witness and
  the bound are checked as registered.
  P5(b) failed in simplification only: with eps carrying no sign assumption sympy left k_sph(equator) as
  1/sqrt((eps - 1)^2). Re-posed as k_sph(equator)^2 = 1/(1 - eps)^2 with k_sph(equator) > 0 at eps = 1/10, 3/10, 1/2.
"""
import sys
import random
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)


# ------------------------------------------------------------------------------------------------------------- P1
print("P1  the ellipsoid at its farthest point")
a1, a2, a3, a4 = sp.symbols("a1:5", positive=True)
y = sp.symbols("y2:5")                      # local graph x1 = a1 sqrt(1 - sum (y_j/a_j)^2) near the vertex
x1 = a1 * sp.sqrt(1 - sum((y[j] / [a2, a3, a4][j])**2 for j in range(3)))
# at the vertex the tangent space is span(e2, e3, e4), the outward normal is e1, and K = -N . Hess x = -Hess x1
Hx1 = sp.Matrix(3, 3, lambda i, j: sp.diff(x1, y[i], y[j]).subs({v: 0 for v in y}))
Kv = (-Hx1).applyfunc(sp.simplify)
want = sp.diag(a1 / a2**2, a1 / a3**2, a1 / a4**2)
ok1 = Kv == want
detK = sp.simplify(Kv.det())
# each a1/a_j^2 >= 1/a1 when a_j <= a1: a1^2 - a_j^2 >= 0
ok1 &= all(sp.simplify(a1 / aj**2 - 1 / a1 - (a1**2 - aj**2) / (a1 * aj**2)) == 0 for aj in (a2, a3, a4))
ok1 &= sp.simplify(detK - a1**3 / (a2 * a3 * a4)**2) == 0
ok1 &= sp.simplify(detK.subs({a2: a1, a3: a1, a4: a1}) - a1**-3) == 0
check("K(vertex) = diag(a1/a_j^2) >= 1/a1; det K = a1^3/(a2 a3 a4)^2 >= a1^-3, equality iff round", ok1, f"K = {list(Kv.diagonal())}")

# ------------------------------------------------------------------------------------------------------------- P2
print("P2  non-convex star-shaped hypersurfaces: the bound at the farthest point")
mp.mp.dps = 40
random.seed(345)


def make_r(coeffs):
    def r(u):
        val = mp.mpf(1)
        for (i, j, k, l), c in coeffs.items():
            val += c * u[0]**i * u[1]**j * u[2]**k * u[3]**l
        return val
    return r


def level_set_K(r, x):
    """Second fundamental form (eigenvalues) of {|x| = r(x/|x|)} at x, outward normal, K = P Hess(Phi) P / |grad Phi|."""
    Phi = lambda z: mp.sqrt(sum(t * t for t in z)) - r([t / mp.sqrt(sum(s * s for s in z)) for t in z])
    h = mp.mpf("1e-12")
    g = [mp.diff(lambda t, i=i: Phi([x[k] + (t if k == i else 0) for k in range(4)]), 0) for i in range(4)]
    H = mp.matrix(4, 4)
    for i in range(4):
        for j in range(4):
            H[i, j] = mp.diff(lambda s, t: Phi([x[k] + (s if k == i else 0) + (t if k == j else 0) for k in range(4)]),
                              (0, 0), (1, 1)) if i != j else mp.diff(lambda t: Phi([x[k] + (t if k == i else 0)
                                                                                   for k in range(4)]), 0, 2)
    gn = mp.sqrt(sum(t * t for t in g))
    n = [t / gn for t in g]
    # orthonormal basis of the tangent space by Gram-Schmidt
    basis = []
    for e in [[1 if k == m else 0 for k in range(4)] for m in range(4)]:
        v = [mp.mpf(e[k]) - sum(n[a] * e[a] for a in range(4)) * n[k] for k in range(4)]
        for b in basis:
            d = sum(v[k] * b[k] for k in range(4))
            v = [v[k] - d * b[k] for k in range(4)]
        nv = mp.sqrt(sum(t * t for t in v))
        if nv > mp.mpf("1e-10") and len(basis) < 3:
            basis.append([t / nv for t in v])
    Kt = mp.matrix(3, 3)
    for a in range(3):
        for b in range(3):
            Kt[a, b] = sum(basis[a][i] * H[i, j] * basis[b][j] for i in range(4) for j in range(4)) / gn
    ev = mp.eig(Kt)[0]
    return sorted([mp.re(e) for e in ev])


def unit(v):
    nv = mp.sqrt(sum(t * t for t in v))
    return [t / nv for t in v]


det2 = []
ok2 = True
monos = [(i, j, k, l) for i in range(5) for j in range(5) for k in range(5) for l in range(5) if 1 <= i + j + k + l <= 4]
for ex in range(3):
    coeffs = {m: mp.mpf(random.randint(-9, 9)) / 200 for m in random.sample(monos, 8)}
    coeffs[(2, 0, 0, 0)] = coeffs.get((2, 0, 0, 0), 0) + mp.mpf(7) / 10      # (7/20)(2 u1^2 - 1)
    coeffs[(0, 0, 0, 0)] = -mp.mpf(7) / 20
    r = make_r(coeffs)
    # the farthest point: maximise r on S^3 (multistart, then refine)
    best = None
    for _ in range(60):
        u0 = unit([mp.mpf(random.gauss(0, 1)) for _ in range(4)])
        val = r(u0)
        if best is None or val > best[0]:
            best = (val, u0)
    u = best[1]
    for it in range(200):                      # projected gradient ascent on S^3
        gr = [mp.diff(lambda t, i=i: r(unit([u[k] + (t if k == i else 0) for k in range(4)])), 0) for i in range(4)]
        u = unit([u[k] + mp.mpf("0.05") * gr[k] for k in range(4)])
    R = r(u)
    xmax = [R * t for t in u]
    kap = level_set_K(r, xmax)
    # a non-convexity witness: some sample point with a negative principal curvature
    neg = False
    for _ in range(400):
        v = unit([mp.mpf(random.gauss(0, 1)) for _ in range(4)])
        if r(v) <= 0:
            continue
        if min(level_set_K(r, [r(v) * t for t in v])) < 0:
            neg = True
            break
    good = min(kap) >= 1 / R - mp.mpf("1e-20") and neg and min(r(unit([mp.mpf(random.gauss(0, 1)) for _ in range(4)]))
                                                             for _ in range(200)) > 0
    ok2 &= good
    det2.append(f"ex{ex}: R = {mp.nstr(R, 6)}, min kappa = {mp.nstr(min(kap), 6)} vs 1/R = {mp.nstr(1 / R, 6)}, "
                f"indefinite somewhere: {neg}")
check("at the maximum of r every principal curvature is >= 1/R, on three non-convex examples", ok2, "; ".join(det2))

# ------------------------------------------------------------------------------------------------------------- P3
print("P3  Jung's constant in R^4")
# regular simplex: the 5 standard basis vectors of R^5 lie in a hyperplane, pairwise distance sqrt(2)
E5 = [sp.Matrix([1 if i == k else 0 for i in range(5)]) for k in range(5)]
cen = sum(E5, sp.zeros(5, 1)) / 5
circ = sp.sqrt(((E5[0] - cen).T * (E5[0] - cen))[0])
edge = sp.sqrt(2)
ratio = sp.simplify(circ / edge)
ok3 = sp.simplify(ratio - sp.sqrt(sp.Rational(2, 5))) == 0
bound = sp.Rational(125, 8) / sp.pi**6
ok3 &= sp.N(bound) < 1
check("circumradius / edge = sqrt(2/5); on the unit S^3 the bound 125/(8 pi^6) < 1 = det M", ok3,
      f"ratio {ratio}; bound {sp.N(bound, 6)}")

# ------------------------------------------------------------------------------------------------------------- P4
print("P4  the degenerate products, and the space forms")


def M_of(g, X):
    n = len(X)
    gi = g.inv()
    Gam = [[[sum(gi[a, d] * (sp.diff(g[d, b], X[c]) + sp.diff(g[d, c], X[b]) - sp.diff(g[b, c], X[d])) for d in range(n)) / 2
             for c in range(n)] for b in range(n)] for a in range(n)]
    Ric = sp.Matrix(n, n, lambda i, j: sp.simplify(sum(sp.diff(Gam[a][i][j], X[a]) - sp.diff(Gam[a][i][a], X[j])
                                                       + sum(Gam[a][a][b] * Gam[b][i][j] - Gam[a][j][b] * Gam[b][i][a]
                                                             for b in range(n)) for a in range(n))))
    Rs = sp.simplify(sum(gi[i, j] * Ric[i, j] for i in range(n) for j in range(n)))
    return (-(gi * (Ric - Rs * g / 2))).applyfunc(sp.simplify)


t, ph, z, rr = sp.symbols("t phi z r", positive=True)
cases = {
    "S^2(r) x S^1": (sp.diag(rr**2, rr**2 * sp.sin(t)**2, 1), [t, ph, z]),
    "H^2 x S^1": (sp.diag(rr**2, rr**2 * sp.sinh(t)**2, 1), [t, ph, z]),
    "flat T^3": (sp.eye(3), [t, ph, z]),
}
det4 = []
ok4 = True
for name, (g, X) in cases.items():
    d = sp.simplify(M_of(g, X).det())
    ok4 &= d == 0
    det4.append(f"{name}: det M = {d}")
chi = sp.symbols("chi", positive=True)
gS3 = sp.diag(1, sp.sin(chi)**2, sp.sin(chi)**2 * sp.sin(t)**2)
MS3 = M_of(gS3, [chi, t, ph])
ok4 &= MS3 == sp.eye(3)
det4.append(f"round S^3: M = Id: {MS3 == sp.eye(3)}")
check("det M = 0 identically on S^2 x S^1, H^2 x S^1, T^3; M = Id on the round S^3 (so on every S^3/Gamma)", ok4,
      "; ".join(det4))

# ------------------------------------------------------------------------------------------------------------- P5
print("P5  the radial graphs r = 1 + eps cos 2chi: where definiteness is lost")
c = sp.symbols("chi")
eps = sp.symbols("epsilon")
rf = 1 + eps * sp.cos(2 * c)
rp, rpp = sp.diff(rf, c), sp.diff(rf, c, 2)
kmer = (rf**2 + 2 * rp**2 - rf * rpp) / (rf**2 + rp**2)**sp.Rational(3, 2)
ksph = (rf - rp * sp.cot(c)) / (rf * sp.sqrt(rf**2 + rp**2))
# (a) against the level-set computation, at eps = 3/10
ok5a = True
det5 = []
e0 = mp.mpf(3) / 10
rfun = lambda u: 1 + e0 * (2 * u[0]**2 - 1)          # cos 2chi = 2 cos^2 chi - 1, cos chi = u1
for cv in (mp.mpf("0.4"), mp.mpf("1.1"), mp.mpf("1.4")):
    u = [mp.cos(cv), mp.sin(cv) * mp.mpf("0.6"), mp.sin(cv) * mp.mpf("0.8"), mp.mpf(0)]
    xp = [rfun(u) * t for t in u]
    kap = level_set_K(rfun, xp)
    km = sp.N(kmer.subs({eps: sp.Rational(3, 10), c: sp.Float(str(cv), 40)}), 40)
    ks = sp.N(ksph.subs({eps: sp.Rational(3, 10), c: sp.Float(str(cv), 40)}), 40)
    want5 = sorted([mp.mpf(str(km)), mp.mpf(str(ks)), mp.mpf(str(ks))])
    err = max(abs(kap[i] - want5[i]) for i in range(3))
    ok5a &= err < mp.mpf("1e-25")
    det5.append(f"chi={mp.nstr(cv, 3)}: err {mp.nstr(err, 3)}")
check("(a) k_mer and k_sph formulas agree with the level-set second fundamental form", ok5a, "; ".join(det5))
# (b) the equator
num_eq = sp.factor(sp.simplify((rf**2 + 2 * rp**2 - rf * rpp).subs(c, sp.pi / 2)))
ks_eq = sp.simplify(ksph.subs(c, sp.pi / 2))
ok5b = (sp.simplify(num_eq - (1 - eps) * (1 - 5 * eps)) == 0 and sp.simplify(ks_eq**2 - 1 / (1 - eps)**2) == 0
        and all(ks_eq.subs(eps, e) > 0 for e in (sp.Rational(1, 10), sp.Rational(3, 10), sp.Rational(1, 2))))
check("(b) at the equator k_sph = 1/(1 - eps) > 0 and the k_mer numerator is (1 - eps)(1 - 5 eps)", ok5b,
      f"numerator {num_eq}, k_sph {ks_eq}")
# (c) on a grid
grid = [sp.Rational(k, 400) * sp.pi for k in range(1, 400)]
fk = sp.lambdify((c, eps), kmer, "mpmath")
fs = sp.lambdify((c, eps), ksph, "mpmath")
mp.mp.dps = 30


def scan(ev):
    vals = []
    for gpt in grid:
        cv = mp.mpf(sp.N(gpt, 30))
        km, ks = fk(cv, ev), fs(cv, ev)
        vals.append((cv, km, ks, km * ks * ks))
    return vals


ok5c = True
det5c = []
for ev in (mp.mpf(1) / 10, mp.mpf(3) / 20):
    vals = scan(ev)
    ok5c &= all(d > 0 and k1 > 0 and k2 > 0 for _, k1, k2, d in vals)
    det5c.append(f"eps={mp.nstr(ev, 3)}: min det K {mp.nstr(min(d for *_, d in vals), 5)}")
vals = scan(mp.mpf(1) / 5)
mind = min(vals, key=lambda v: v[3])
ok5c &= abs(mind[3]) < mp.mpf("1e-25") and abs(mind[0] - mp.pi / 2) < mp.mpf("1e-25")
det5c.append(f"eps=1/5: min det K {mp.nstr(mind[3], 3)} at chi = {mp.nstr(mind[0], 6)}")
vals = scan(mp.mpf(3) / 10)
eq = [v for v in vals if abs(v[0] - mp.pi / 2) < mp.mpf("1e-25")][0]
ok5c &= eq[3] < 0
first_zero = next(v[0] for v in vals if v[1] <= 0)
ok5c &= all(v[1] > 0 and v[2] > 0 for v in vals if v[0] < first_zero)
det5c.append(f"eps=3/10: det K(equator) {mp.nstr(eq[3], 5)}; definite for chi < {mp.nstr(first_zero, 5)}")
check("(c) definite with det K > 0 below eps = 1/5; det K = 0 first at the equator at eps = 1/5; indefinite only past "
      "a zero of det K", ok5c, "; ".join(det5c))

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