Trac #29189: Migrate is_lawrence_polytope and is_pyramid to combi…

…natorial polyhedron

This ticket migrates the methods `is_lawrence_polytope` and `is_pyramid`
from `Polyhedron_base` to `CombinatorialPolyhedron`.

Also, we change the output for the `0`-dimensional polyhedron. It is a
pyramid over the empty polyhedron, even if it is not constructable in
sage.

Along the way we fix a small bug, where the trivial combinatorial
polyhedron in dimension 0 was set up without vertices (and facets). The
bug fix is tested by
`CombinatorialPolyhedron(0).is_pyramid(certificate=True)`.

URL: https://trac.sagemath.org/29189
Reported by: gh-kliem
Ticket author(s): Jonathan Kliem
Reviewer(s): Laith Rastanawi

src/sage/geometry/polyhedron/base.py

@@ -3171,21 +3171,22 @@ def is_pyramid(self, certificate=False):
             sage: Q = polytopes.octahedron()
             sage: Q.is_pyramid()
             False
+
+        For the `0`-dimensional polyhedron, the output is ``True``,
+        but it cannot be constructed as a pyramid over the empty polyhedron::
+
+            sage: P = Polyhedron([[0]])
+            sage: P.is_pyramid()
+            True
+            sage: Polyhedron().pyramid()
+            Traceback (most recent call last):
+            ...
+            ZeroDivisionError: rational division by zero
         """
         if not self.is_compact():
             raise ValueError("polyhedron has to be compact")
 
-        # Find a vertex that is incident to all elements in Hrepresentation but one.
-        IM = self.incidence_matrix()
-        for index in range(self.n_vertices()):
-            vertex_incidences = IM.row(index)
-            if sum(vertex_incidences) == IM.ncols() - 1:
-                if certificate:
-                    return (True, self.vertices()[index])
-                return True
-        if certificate:
-            return (False, None)
-        return False
+        return self.combinatorial_polyhedron().is_pyramid(certificate)
 
     def is_bipyramid(self, certificate=False):
         r"""
@@ -5411,18 +5412,7 @@ def is_lawrence_polytope(self):
         if not self.is_compact():
             raise NotImplementedError("self must be a polytope")
 
-        n = self.n_vertices()
-        vertices = list(range(n))
-        facets = self.incidence_matrix().columns()
-
-        for facet in facets:
-            facet_vertices = facet.nonzero_positions()
-            if len(facet_vertices) == n-1 or len(facet_vertices) == n-2:
-                facet_non_vertices = [i for i in range(n) if i not in facet_vertices]
-                if all(vertex in vertices for vertex in facet_non_vertices):
-                    for vertex in facet_non_vertices:
-                        vertices.remove(vertex)
-        return not vertices
+        return self.combinatorial_polyhedron().is_lawrence_polytope()
 
     def barycentric_subdivision(self, subdivision_frac=None):
         r"""

src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx

@@ -446,9 +446,15 @@ cdef class CombinatorialPolyhedron(SageObject):
             # one can give an Integer as Input.
             if data < -1:
                 ValueError("any polyhedron must have dimension at least -1")
-            self._n_facets = 0
             self._dimension = data
 
+            if self._dimension == 0:
+                self._n_facets = 1
+                self._n_Vrepresentation = 1
+            else:
+                self._n_facets = 0
+                self._n_Vrepresentation = 0
+
             # Initializing the facets in their Bit-representation.
             self._bitrep_facets = facets_tuple_to_bit_rep_of_facets((), 0)
 
@@ -836,7 +842,7 @@ cdef class CombinatorialPolyhedron(SageObject):
                 # The Polyhedron has no vertex.
                 return ()
         if names and self.Vrep():
-            return tuple(self.Vrep()[i]     for i in range(self.n_Vrepresentation()) if not i in self.far_face_tuple())
+            return tuple(self.Vrep()[i]  for i in range(self.n_Vrepresentation()) if not i in self.far_face_tuple())
         else:
             return tuple(smallInteger(i) for i in range(self.n_Vrepresentation()) if not i in self.far_face_tuple())
 
@@ -1764,6 +1770,191 @@ cdef class CombinatorialPolyhedron(SageObject):
                 d = dim
         return smallInteger(simplicity)
 
+    @cached_method
+    def is_lawrence_polytope(self):
+        """
+        Return ``True`` if ``self`` is a Lawrence polytope.
+
+        A polytope is called a Lawrence polytope if it has a centrally
+        symmetric (normalized) Gale diagram.
+
+        Equivalently, there exists a partition `P_1,\dots,P_k`
+        of the vertices `V` such that each part
+        `P_i` has size `2` or `1` and for each part there exists
+        a facet with vertices exactly `V \setminus P_i`.
+
+        EXAMPLES::
+
+            sage: C = polytopes.simplex(5).combinatorial_polyhedron()
+            sage: C.is_lawrence_polytope()
+            True
+            sage: P = polytopes.hypercube(4).lawrence_polytope()
+            sage: C = P.combinatorial_polyhedron()
+            sage: C.is_lawrence_polytope()
+            True
+            sage: P = polytopes.hypercube(4)
+            sage: C = P.combinatorial_polyhedron()
+            sage: C.is_lawrence_polytope()
+            False
+
+        For unbounded polyhedra, an error is raised::
+
+            sage: C = CombinatorialPolyhedron([[0,1], [0,2]], far_face=[1,2], unbounded=True)
+            sage: C.is_lawrence_polytope()
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: this function is implemented for polytopes only
+
+        AUTHORS:
+
+        - Laith Rastanawi
+        - Jonathan Kliem
+
+        REFERENCES:
+
+            For more information, see [BaSt1990]_.
+        """
+        if not self.is_compact():
+            raise NotImplementedError("this function is implemented for polytopes only")
+        if self.n_Vrepresentation() <= 2:
+            return True
+
+        cdef FaceIterator facet_iterator = self._face_iter(False, self.dimension()-1)
+        cdef CombinatorialFace facet
+        cdef size_t n_vertices = self.n_Vrepresentation()
+        cdef size_t one, two, length, counter
+        cdef list vertices = [1 for _ in range(n_vertices)]
+
+        for facet in facet_iterator:
+            length = facet.n_atom_rep()
+            if length >= n_vertices - 2:
+                # The facet has at most two non-vertices and corresponds to
+                # two symmetric vertices or a vertex at the origin
+                # in the Gale transform.
+                facet.set_atom_rep()
+                counter = 0
+                while counter < length:
+                    if facet.atom_rep[counter] != counter:
+                        # We have found our first non-vertex.
+                        one = counter
+                        break
+                    counter += 1
+                else:
+                    # The facet contains the first ``length`` vertices.
+                    one = length
+
+                if length == n_vertices - 1:
+                    # The facet corresponds to a vertex at the origin
+                    # of the Gale transform.
+                    vertices[one] = 0
+                else:
+                    # The facet corresponds to two symmetric vertices
+                    # of the Gale transform.
+                    while counter < length:
+                        if facet.atom_rep[counter] != counter + 1:
+                            # We have found our second non-vertex.
+                            two = counter + 1
+                            break
+                        counter += 1
+                    else:
+                        # The second non-vertex is the very last vertex.
+                        two = length + 1
+
+                    if vertices[one] == vertices[two]:
+                        # Possibly the Gale transform contains duplicates,
+                        # we must make sure that the mulitplicites are symmetric as well.
+                        # (And not two vertices are symmetric to just one).
+                        vertices[one] = 0
+                        vertices[two] = 0
+
+        return not any(vertices)
+
+    @cached_method
+    def is_pyramid(self, certificate=False):
+        r"""
+        Test whether the polytope is a pyramid over one of its facets.
+
+        INPUT:
+
+        - ``certificate`` -- boolean (default: ``False``); specifies whether
+          to return a vertex of the polytope which is the apex of a pyramid,
+          if found
+
+        OUTPUT:
+
+        If ``certificate`` is ``True``, returns a tuple containing:
+
+        1. Boolean.
+        2. The apex of the pyramid or ``None``.
+
+        If ``certificate`` is ``False`` returns a boolean.
+
+        AUTHORS:
+
+        - Laith Rastanawi
+        - Jonathan Kliem
+
+        EXAMPLES::
+
+            sage: C = polytopes.cross_polytope(4).combinatorial_polyhedron()
+            sage: C.is_pyramid()
+            False
+            sage: C.is_pyramid(certificate=True)
+            (False, None)
+            sage: C = polytopes.cross_polytope(4).pyramid().combinatorial_polyhedron()
+            sage: C.is_pyramid()
+            True
+            sage: C.is_pyramid(certificate=True)
+            (True, A vertex at (0, -1, 0, 0, 0))
+            sage: C = polytopes.simplex(5).combinatorial_polyhedron()
+            sage: C.is_pyramid(certificate=True)
+            (True, A vertex at (0, 0, 0, 0, 0, 1))
+
+        For unbounded polyhedra, an error is raised::
+
+            sage: C = CombinatorialPolyhedron([[0,1], [0,2]], far_face=[1,2], unbounded=True)
+            sage: C.is_pyramid()
+            Traceback (most recent call last):
+            ...
+            ValueError: polyhedron has to be compact
+
+        TESTS::
+
+            sage: CombinatorialPolyhedron(-1).is_pyramid()
+            False
+            sage: CombinatorialPolyhedron(-1).is_pyramid(True)
+            (False, None)
+            sage: CombinatorialPolyhedron(0).is_pyramid()
+            True
+            sage: CombinatorialPolyhedron(0).is_pyramid(True)
+            (True, 0)
+        """
+        if not self.is_bounded():
+            raise ValueError("polyhedron has to be compact")
+
+        if self.dim() == -1:
+            if certificate:
+                return (False, None)
+            return False
+
+        if self.dim() == 0:
+            if certificate:
+                return (True, self.Vrepresentation()[0])
+            return True
+
+        # Find a vertex that is incident to all elements in Hrepresentation but one.
+        vertex_iter = self._face_iter(True, 0)
+        n_facets = self.n_facets()
+        for index, vertex in enumerate(vertex_iter):
+            if vertex.n_ambient_Hrepresentation() == n_facets - 1:
+                if certificate:
+                    return (True, self.Vrepresentation()[index])
+                return True
+
+        if certificate:
+            return (False, None)
+        return False
+
     def face_iter(self, dimension=None, dual=None):
         r"""
         Iterator over all proper faces of specified dimension.