#!/usr/bin/env python3
"""Supergeometry preprint, Sec. 3: a regular tetrahedron in an octahedral cavity.

Cavity: octahedron with vertices at +-e_i scaled so that its edge is a = 1 (|x|+|y|+|z| <= 1/sqrt(2)).
Checks:
  1. A tetrahedron of the SAME edge fits exactly (maximal uniform scale = 1.0000): base on one face,
     apex at the centre of the opposite face.
  2. It fills exactly 1/4 of the cavity volume.
  3. It fits in 8 positions (one per face), all equivalent under the octahedral group.
  4. Equal-edge fit in random orientations is < 1 except at these 8 orientations.
  5. A bipyramid (two face-sharing tetrahedra) has end-to-end length 2/sqrt(3) times the cavity diagonal.

Run:  python3 cavity_geometry.py       (needs numpy, scipy)
"""
import itertools
import numpy as np
from scipy.optimize import linprog
from scipy.spatial.transform import Rotation as R

RO = 1 / np.sqrt(2)  # octaradius for edge 1
SIGNS = np.array(list(itertools.product([1, -1], repeat=3)), float)


def max_scale(V):
    """Largest s, centre c with |s*v + c|_(octahedral norm) <= RO for all vertices v (linear program)."""
    A, b = [], []
    for v in V:
        for s in SIGNS:
            A.append([s @ v, *s])
            b.append(RO)
    r = linprog([-1, 0, 0, 0], A_ub=A, b_ub=b, bounds=[(0, None)] + [(None, None)] * 3)
    return -r.fun


def tet_on_face(sig):
    base = RO * np.diag(sig)                  # three cavity vertices of the face with signs sig
    apex = -RO / 3 * np.asarray(sig, float)   # centre of the opposite face
    T = np.vstack([base, apex])
    return T - T.mean(0)


if __name__ == "__main__":
    T = tet_on_face((1, 1, 1))
    print("tetrahedron edge:", round(np.linalg.norm(T[0] - T[1]), 6), "(cavity edge = 1)")
    print("max uniform scale that still fits:", round(max_scale(T), 6))
    vt, vo = np.sqrt(2) / 12, np.sqrt(2) / 3
    print("volume fraction of the cavity:", round(vt / vo, 6))
    print("scale for each of the 8 face positions:", [round(max_scale(tet_on_face(s)), 6) for s in SIGNS])
    rng = np.random.default_rng(0)
    vals = [max_scale(T @ R.random(random_state=rng).as_matrix().T) for _ in range(3000)]
    print("random orientations: max %.4f, median %.4f, min %.4f" % (max(vals), np.median(vals), min(vals)))
    apex2 = 2 * T[1:].mean(0) - T[0]
    print("bipyramid end-to-end / cavity diagonal:", round(np.linalg.norm(T[0] - apex2) / (2 * RO), 6),
          "(2/sqrt(3) = %.6f)" % (2 / np.sqrt(3)))