Maximal angles between convex cones

src/doc/en/reference/discrete_geometry/index.rst

@@ -80,6 +80,7 @@ Toric geometry
    sage/geometry/toric_lattice
    sage/geometry/cone
    sage/geometry/cone_catalog
+   sage/geometry/cone_critical_angles
    sage/geometry/fan
    sage/geometry/fan_morphism
    sage/geometry/point_collection

src/doc/en/reference/references/index.rst

@@ -5120,6 +5120,13 @@ REFERENCES:
              Electronic Journal of Linear Algebra, 34:444-458, 2018,
              :doi:`10.13001/1081-3810.3782`.
 
+.. [Or2020] \M. Orlitzky. When a maximal angle among cones is nonobtuse.
+            Computational and Applied Mathematics 39(2), 2020.
+            :doi:`10.1007/s40314-020-1115-y`.
+
+.. [Or2024] \M. Orlitzky. Continuity of the conic hull.
+            Journal of Convex Analysis 31(1):255-264, 2024.
+
 .. [ORS2013] Peter Ozsvath, Jacob Rasmussen, Zoltan Szabo,
              *Odd Khovanov homology*,
              Algebraic & Geometric Topology 13 (2013) 1465–1488,

src/sage/geometry/cone.py

@@ -6197,6 +6197,201 @@ def Z_operators_gens(self):
         """
         return [ -cp for cp in self.cross_positive_operators_gens() ]
 
+    def max_angle(self, other=None, exact=True, epsilon=0):
+        r"""
+        Return the maximal angle between ``self`` and ``other``.
+
+        The maximal angle between two closed convex cones is the
+        unique largest angle formed by any two unit-norm vectors in
+        those cones. In pathological cases, this computation can fail.
+
+        If it fails when ``exact`` is ``True`` and if each of the
+        cones :meth:`is_strictly_convex`, then a second attempt will
+        be made using inexact arithmetic. (This sometimes avoids the
+        problem noted in [Or2024]_). If the computation fails when the
+        cones are not strictly convex or when ``exact`` is ``False``,
+        a :class:`ValueError` is raised.
+
+        INPUT:
+
+        - ``other`` -- (default: ``None``) a rational, polyhedral
+          convex cone
+
+        - ``exact`` -- (default: ``True``) whether or not to use exact
+          rational arithmetic instead of floating point computations;
+          beware that ``True`` is not guaranteed to avoid floating
+          point computations if the algorithm runs into trouble in
+          rational arithmetic
+
+        - ``epsilon`` -- (default: ``0``) the tolerance to use when
+          making comparisons
+
+        .. WARNING::
+
+            Using inexact arithmetic (``exact=False``) is faster, but
+            this computation is only known to be stable when both of
+            the cones are strictly convex (or when one of them is the
+            entire space, but the maximal angle is obviously `\pi` in
+            that case).
+
+        OUTPUT:
+
+        A triple `\left( \theta_{\max}, u, v \right)`
+        containing:
+
+        - the maximal angle `\theta_{\max}` between ``self`` and
+          ``other``
+
+        - a vector `u` in ``self`` that achieves the maximal angle
+
+        - a vector `v` in ``other`` that achieves the maximal angle
+
+        If ``other`` is ``None`` (the default), then the maximal angle
+        within this cone (between this cone and itself) is returned.
+
+        If an eigenspace of dimension greater than one is encountered
+        and if the corresponding angle cannot be ruled out as a
+        maximum, the behavior of this function depends on
+        ``exact``:
+
+        - If ``exact`` is ``True`` and if both ``self`` and ``other``
+          are strictly convex, then the algorithm will fall back to
+          inexact arithmetic. In that case, the returned angle and
+          vectors will be over :class:`sage.rings.real_double.RDF`.
+
+        - If ``exact`` is ``False`` or if either cone is not strictly
+          convex, then a :class:`ValueError` is raised to indicate
+          that we have failed; i.e. we cannot say with certainty what
+          the maximal angle is.
+
+        REFERENCES:
+
+        - [IS2005]_
+        - [SS2016]_
+        - [Or2020]_
+        - [Or2024]_
+
+        ALGORITHM:
+
+        Algorithm 3 in [Or2020]_ is used. If a potentially-maximal
+        angle corresponds to an eigenspace of dimension two or more,
+        we sometimes fall back to inexact arithmetic which has the
+        effect of perturbing the cones. That this will not affect the
+        answer too much is one conclusion of [Or2024]_.
+
+        EXAMPLES:
+
+        The maximal angle in a single ray is zero::
+
+            sage: K = random_cone(min_rays=1, max_rays=1, max_ambient_dim=5)
+            sage: K.max_angle()[0]
+            0
+
+        The maximal angle in the nonnegative octant is `\pi/2`::
+
+            sage: K = cones.nonnegative_orthant(3)
+            sage: K.max_angle()[0]
+            1/2*pi
+
+        The maximal angle between the nonnegative quintant and the
+        Schur cone of dimension 5 is about `0.8524 \pi`. The same
+        result can be obtained faster using inexact arithmetic, but
+        only confidently so because we already know the answer::
+
+            sage: # long time
+            sage: P = cones.nonnegative_orthant(5)
+            sage: Q = cones.schur(5)
+            sage: actual = P.max_angle(Q)[0]
+            sage: expected = 0.8524*pi
+            sage: bool( (actual - expected).abs() < 0.0001 )
+            True
+            sage: actual = P.max_angle(Q,exact=False)[0]
+            sage: bool( (actual - expected).abs() < 0.0001 )
+            True
+
+        The maximal angle within the Schur cone is known explicitly
+        via Gourion and Seeger's Proposition 2 [GS2010]_::
+
+            sage: n = 3
+            sage: K = cones.schur(n)
+            sage: bool(K.max_angle()[0] == ((n-1)/n)*pi)
+            True
+
+        Sage can't prove that the actual and expected results are
+        equal in the next two cases without a little nudge in the
+        right direction, and, moreover, it's crashy about it::
+
+            sage: n = 4
+            sage: K = cones.schur(n)
+            sage: actual = K.max_angle()[0].simplify()._sympy_()._sage_()
+            sage: expected = ((n-1)/n)*pi
+            sage: bool( actual == expected )
+            True
+            sage: n = 5
+            sage: K = cones.schur(n)
+            sage: actual = K.max_angle()[0].simplify()._sympy_()._sage_()
+            sage: expected = ((n-1)/n)*pi
+            sage: bool( actual == expected )
+            True
+
+        When there's a unit norm vector in ``self`` whose negation is
+        in ``other``, they form a maximal angle of `\pi`::
+
+            sage: P = Cone([(5,1), (1,-1)])
+            sage: Q = Cone([(-1,0), (-1,0)])
+            sage: P.max_angle(Q)[0]
+            pi
+
+        TESTS:
+
+        When ``self`` and ``-other`` have nontrivial intersection, we
+        expect the maximal angle to be `\pi`::
+
+            sage: # long time
+            sage: from sage.geometry.cone_critical_angles import (
+            ....:   _random_admissible_cone )
+            sage: n = ZZ.random_element(1,3)
+            sage: P = _random_admissible_cone(ambient_dim=n)
+            sage: Q = _random_admissible_cone(ambient_dim=n)
+            sage: expected = P.intersection(-Q).is_trivial()
+            sage: actual = bool(P.max_angle(Q)[0] != pi)
+            sage: actual == expected
+            True
+
+        In particular, this should happen when either cone is the full
+        space::
+
+            sage: from sage.geometry.cone_critical_angles import (
+            ....:   _random_admissible_cone )
+            sage: n = ZZ.random_element(1,5)
+            sage: P = _random_admissible_cone(ambient_dim=n)
+            sage: Q = cones.trivial(n, P.dual().lattice()).dual()
+            sage: Q.is_full_space()
+            True
+            sage: P.max_angle(Q)[0]
+            pi
+            sage: Q.max_angle(P)[0]
+            pi
+        """
+        # We do the argument checking here, in the public cone method,
+        # because the error message should say something like "this
+        # cone" and not reference the argument of a function the user
+        # has never heard of in some internal module.
+        if self.is_trivial():
+            raise ValueError("cone should not be trivial")
+
+        if other is None:
+            other = self
+        else:
+            if (other.lattice_dim() != self.lattice_dim()):
+                raise ValueError("lattice dimensions of self and other "
+                                 "must agree")
+            if other.is_trivial():
+                raise ValueError("other cone cannot be trivial")
+
+        from sage.geometry.cone_critical_angles import max_angle
+        return max_angle(self, other, exact, epsilon)
+
 
 def random_cone(lattice=None, min_ambient_dim=0, max_ambient_dim=None,
                 min_rays=0, max_rays=None, strictly_convex=None, solid=None):