Add deformation cones and checking for regularity for Point Configurations and normal fans of Polyhedra

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

@@ -104,6 +104,10 @@ REFERENCES:
              graphs and isoperimetric inequalities*, The Annals of Probability
              32 (2004), no.  3A, 1727-1745.
 
+.. [ACEP2020] Federico Ardila, Federico Castillo, Christopher Eur, Alexander Postnikov,
+         *Coxeter submodular functions and deformations of Coxeter permutahedra*,
+         Advances in Mathematics, Volume 365, 13 May 2020.
+
 .. [ALL2002] P. Auger, G. Labelle and P. Leroux, *Combinatorial
              addition formulas and applications*, Advances in Applied
              Mathematics 28 (2002) 302-342.

src/sage/geometry/fan.py

@@ -2439,6 +2439,69 @@ def Gale_transform(self):
         m = m.augment(matrix(ZZ, m.nrows(), 1, [1] * m.nrows()))
         return matrix(ZZ, m.integer_kernel().matrix())
 
+    def is_polytopal(self) -> bool:
+        r"""
+        Check if ``self`` is the normal fan of a polytope.
+
+        A rational polyhedral fan is *polytopal* if it is the normal fan of a
+        polytope. This is also called *regular*, or provides a *coherent*
+        subdivision or leads to a *projective* toric variety.
+
+        OUTPUT: ``True`` if ``self`` is polytopal and ``False`` otherwise
+
+        EXAMPLES:
+
+        This is the mother of all examples (see Section 7.1.1 in
+        [DLRS2010]_)::
+
+            sage: def mother(epsilon=0):
+            ....:     rays = [(4-epsilon,epsilon,0),(0,4-epsilon,epsilon),(epsilon,0,4-epsilon),(2,1,1),(1,2,1),(1,1,2),(-1,-1,-1)]
+            ....:     L = [(0,1,4),(0,3,4),(1,2,5),(1,4,5),(0,2,3),(2,3,5),(3,4,5),(6,0,1),(6,1,2),(6,2,0)]
+            ....:     S1 = [Cone([rays[i] for i in indices]) for indices in L]
+            ....:     return Fan(S1)
+
+        When epsilon=0, it is not polytopal::
+
+            sage: epsilon = 0
+            sage: mother(epsilon).is_polytopal()
+            False
+
+        Doing a slight perturbation makes the same subdivision polytopal::
+
+            sage: epsilon = 1/2
+            sage: mother(epsilon).is_polytopal()
+            True
+
+        TESTS::
+
+            sage: cone = Cone([(1,1), (2,1)])
+            sage: F = Fan([cone])
+            sage: F.is_polytopal()
+            Traceback (most recent call last):
+            ...
+            ValueError: to be polytopal, the fan should be complete
+
+        .. SEEALSO::
+
+            :meth:`is_projective`.
+        """
+        if not self.is_complete():
+            raise ValueError('to be polytopal, the fan should be complete')
+        from sage.geometry.triangulation.point_configuration import PointConfiguration
+        from sage.geometry.polyhedron.constructor import Polyhedron
+        pc = PointConfiguration(self.rays())
+        v_pc = [tuple(p) for p in pc]
+        pc_to_indices = {tuple(p):i for i, p in enumerate(pc)}
+        indices_to_vr = (tuple(r) for r in self.rays())
+        cone_indices = (cone.ambient_ray_indices() for cone in self.generating_cones())
+        translator = [pc_to_indices[t] for t in indices_to_vr]
+        translated_cone_indices = [[translator[i] for i in ci] for ci in cone_indices]
+        dc_pc = pc.deformation_cone(translated_cone_indices)
+        lift = dc_pc.an_element()
+        ieqs = [(lift_i,) + v for (lift_i, v) in zip(lift, v_pc)]
+        poly = Polyhedron(ieqs=ieqs)
+        return self.is_equivalent(poly.normal_fan())
+
     def generating_cone(self, n):
         r"""
         Return the ``n``-th generating cone of ``self``.

src/sage/geometry/polyhedron/base5.py

@@ -659,6 +659,72 @@ def lawrence_polytope(self):
         parent = self.parent().change_ring(self.base_ring(), ambient_dim=self.ambient_dim() + n)
         return parent.element_class(parent, [lambda_V, [], []], None)
 
+    def deformation_cone(self):
+        r"""
+        Return the deformation cone of ``self``.
+
+        Let `P` be a `d`-polytope in `\RR^r` with `n` facets. The deformation
+        cone is a polyhedron in `\RR^n` whose points are the right-hand side `b`
+        in `Ax\leq b` where `A` is the matrix of facet normals of ``self``, so
+        that the resulting polytope has a normal fan which is a coarsening of
+        the normal fan of ``self``.
+
+        EXAMPLES:
+
+        Let's examine the deformation cone of the square with one truncated
+        vertex::
+
+            sage: tc = Polyhedron([(1, -1), (1/3, 1), (1, 1/3), (-1, 1), (-1, -1)])
+            sage: dc = tc.deformation_cone()
+            sage: dc.an_element()
+            (2, 1, 1, 0, 0)
+            sage: [_.A() for _ in tc.Hrepresentation()]
+            [(1, 0), (0, 1), (0, -1), (-3, -3), (-1, 0)]
+            sage: P = Polyhedron(rays=[(1, 0, 2), (0, 1, 1), (0, -1, 1), (-3, -3, 0), (-1, 0, 0)])
+            sage: P.rays()
+            (A ray in the direction (-1, -1, 0),
+             A ray in the direction (-1, 0, 0),
+             A ray in the direction (0, -1, 1),
+             A ray in the direction (0, 1, 1),
+             A ray in the direction (1, 0, 2))
+
+        Now, let's compute the deformation cone of the pyramid over a square
+        and verify that it is not full dimensional::
+
+            sage: py = Polyhedron([(0, -1, -1), (0, -1, 1), (0, 1, -1), (0, 1, 1), (1, 0, 0)])
+            sage: dc_py = py.deformation_cone(); dc_py
+            A 4-dimensional polyhedron in QQ^5 defined as the convex hull of 1 vertex, 1 ray, 3 lines
+            sage: [ineq.b() for ineq in py.Hrepresentation()]
+            [0, 1, 1, 1, 1]
+            sage: r = dc_py.rays()[0]
+            sage: l1,l2,l3 = dc_py.lines()
+            sage: r.vector()-l1.vector()/2-l2.vector()-l3.vector()/2
+            (0, 1, 1, 1, 1)
+
+        .. SEEALSO::
+
+            :meth:`~sage.schemes.toric.variety.Kaehler_cone`
+
+        REFERENCES:
+
+        For more information, see Section 5.4 of [DLRS2010]_ and Section
+        2.2 of [ACEP2020].
+        """
+        from .constructor import Polyhedron
+        m = matrix([ineq.A() for ineq in self.Hrepresentation()])
+        m = m.transpose()
+        m_ker = m.right_kernel_matrix(basis='computed')
+        gale = tuple(m_ker.columns())
+        collection = (f.ambient_H_indices() for f in self.faces(0))
+        n = len(gale)
+        c = None
+        for cone_indices in collection:
+            dual_cone = Polyhedron(rays=[gale[i] for i in range(n) if i not in
+                                         cone_indices])
+            c = c.intersection(dual_cone) if c is not None else dual_cone
+        preimages = [m_ker.solve_right(r.vector()) for r in c.rays()]
+        return Polyhedron(lines=m.rows(), rays=preimages)
+
     ###########################################################
     # Binary operations.
     ###########################################################