implemented is_pyramid, is_bipyramid, is_prism

Author: Jonathan Kliem

src/doc/en/thematic_tutorials/geometry/polyhedra_quickref.rst

@@ -55,6 +55,9 @@ List of Polyhedron methods
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.is_simplicial` |  checks whether all faces of the polyhedron are simplices
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.is_lawrence_polytope` |  tests whether self is a Lawrence polytope
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.is_self_dual` |  tests whether the polytope is self-dual
+    :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.is_pyramid` | test whether the polytope is a pyramid over one of its facets
+    :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.is_bipyramid` | test whether the polytope is combinatorially equivalent to a bipyramid over some polytope
+    :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.is_prism` | test whether the polytope is combinatorially equivalent to a prism of some polytope
 
 **Enumerative properties**
 

src/sage/geometry/polyhedron/base.py

@@ -2685,6 +2685,309 @@ def is_simplicial(self):
                    for facet in self.Hrepresentation()
                    if not facet.is_equation())
 
+    def is_pyramid(self, certificate=False):
+        """
+        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.
+
+        EXAMPLES::
+
+            sage: P = polytopes.simplex(3)
+            sage: P.is_pyramid()
+            True
+            sage: P.is_pyramid(certificate=True)
+            (True, A vertex at (0, 0, 0, 1))
+            sage: egyptian_pyramid = polytopes.regular_polygon(4).pyramid()
+            sage: egyptian_pyramid.is_pyramid()
+            True
+            sage: Q = polytopes.octahedron()
+            sage: Q.is_pyramid()
+            False
+        """
+        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
+
+    def is_bipyramid(self, certificate=False):
+        r"""
+        Test whether the polytope is combinatorially equivalent to a
+        bipyramid over some polytope.
+
+        INPUT:
+
+        - ``certificate`` -- boolean (default: ``False``); specifies whether
+          to return two vertices of the polytope which are the apices of a
+          bipyramid, if found
+
+        OUTPUT:
+
+        If ``certificate`` is ``True``, returns a tuple containing:
+
+        1. Boolean.
+        2. ``None`` or a tuple containing:
+            a. The first apex.
+            b. The second apex.
+
+        If ``certificate`` is ``False`` returns a boolean.
+
+        EXAMPLES::
+
+            sage: P = polytopes.octahedron()
+            sage: P.is_bipyramid()
+            True
+            sage: P.is_bipyramid(certificate=True)
+            (True, [A vertex at (-1, 0, 0), A vertex at (1, 0, 0)])
+            sage: Q = polytopes.cyclic_polytope(3,7)
+            sage: Q.is_bipyramid()
+            False
+            sage: R = Q.bipyramid()
+            sage: R.is_bipyramid(certificate=True)
+            (True, [A vertex at (-1, 3, 13, 63), A vertex at (1, 3, 13, 63)])
+
+        TESTS::
+
+            sage: P = polytopes.permutahedron(4).bipyramid()
+            sage: P.is_bipyramid()
+            True
+
+            sage: P = polytopes.cube()
+            sage: P.is_bipyramid()
+            False
+
+            sage: P = Polyhedron(vertices=[[0,1], [1,0]], rays=[[1,1]])
+            sage: P.is_bipyramid()
+            Traceback (most recent call last):
+            ...
+            ValueError: polyhedron has to be compact
+
+        ALGORITHM:
+
+        Assume all faces of a polyhedron to be given as lists of vertices.
+
+        A polytope is a bipyramid with apexes `v`, `w` if and only if for each
+        proper face `v \in F` there exists a face `G` with
+        `G \setminus \{w\} = F \setminus \{v\}`
+        and vice versa (for each proper face
+        `w \in F` there exists ...).
+
+        To check this property it suffices to check for all facets of the polyhedron.
+        """
+        if not self.is_compact():
+            raise ValueError("polyhedron has to be compact")
+
+        from sage.misc.functional import is_odd
+        n_verts = self.n_vertices()
+        n_facets = self.n_facets()
+        if is_odd(n_facets):
+            if certificate:
+                return (False, None)
+            return False
+
+        IM = self.incidence_matrix()
+        if self.n_equations():
+            # Remove equations from the incidence matrix,
+            # such that this is the vertex-facet incidences matrix.
+            I1 = IM.transpose()
+            I2 = I1[[i for i in range(self.n_Hrepresentation())
+                     if not self.Hrepresentation()[i].is_equation()]]
+            IM = I2.transpose()
+
+        facets_incidences = [set(column.nonzero_positions()) for column in IM.columns()]
+        verts_incidences = dict()
+        for i in range(n_verts):
+            v_i = set(IM.row(i).nonzero_positions())
+            if len(v_i) == n_facets/2:
+                verts_incidences[i] = v_i
+
+        # Find two vertices ``vert1`` and ``vert2`` such that one of them
+        # lies on exactly half of the facets, and the other one lies on
+        # exactly the other half.
+        from itertools import combinations
+        for index1, index2 in combinations(verts_incidences, 2):
+            vert1_incidences = verts_incidences[index1]
+            vert2_incidences = verts_incidences[index2]
+            vert1and2 = vert1_incidences.union(vert2_incidences)
+            if len(vert1and2) == n_facets:
+                # We have found two candidates for apexes.
+                # Remove from each facet ``index1`` resp. ``index2``.
+                test_facets = set(frozenset(facet_inc.difference({index1, index2}))
+                                  for facet_inc in facets_incidences)
+                if len(test_facets) == n_facets/2:
+                    # For each `F` containing `index1` there is
+                    # `G` containing `index2` such that
+                    # `F \setminus \{index1\} =  G \setminus \{index2\}
+                    # and vice versa.
+                    if certificate:
+                        V = self.vertices()
+                        return (True, [V[index1], V[index2]])
+                    return True
+
+        if certificate:
+            return (False, None)
+        return False
+
+    def is_prism(self, certificate=False):
+        """
+        Test whether the polytope is combinatorially equivalent to a prism of
+        some polytope.
+
+        INPUT:
+
+        - ``certificate`` -- boolean (default: ``False``); specifies whether
+          to return two facets of the polytope which are the bases of a prism,
+          if found
+
+        OUTPUT:
+
+        If ``certificate`` is ``True``, returns a tuple containing:
+
+        1. Boolean.
+        2. ``None`` or a tuple containing:
+            a. List of the vertices of the first base facet.
+            b. List of the vertices of the second base facet.
+
+        If ``certificate`` is ``False`` returns a boolean.
+
+        EXAMPLES::
+
+            sage: P = polytopes.cube()
+            sage: P.is_prism()
+            True
+            sage: P.is_prism(certificate=True)
+            (True,
+             [[A vertex at (-1, -1, 1),
+               A vertex at (-1, 1, 1),
+               A vertex at (1, -1, 1),
+               A vertex at (1, 1, 1)],
+              [A vertex at (-1, -1, -1),
+               A vertex at (-1, 1, -1),
+               A vertex at (1, -1, -1),
+               A vertex at (1, 1, -1)]])
+            sage: Q = polytopes.cyclic_polytope(3,8)
+            sage: Q.is_prism()
+            False
+            sage: R = Q.prism()
+            sage: R.is_prism(certificate=True)
+            (True,
+             [[A vertex at (1, 0, 0, 0),
+               A vertex at (1, 1, 1, 1),
+               A vertex at (1, 2, 4, 8),
+               A vertex at (1, 3, 9, 27),
+               A vertex at (1, 4, 16, 64),
+               A vertex at (1, 5, 25, 125),
+               A vertex at (1, 6, 36, 216),
+               A vertex at (1, 7, 49, 343)],
+              [A vertex at (0, 0, 0, 0),
+               A vertex at (0, 1, 1, 1),
+               A vertex at (0, 2, 4, 8),
+               A vertex at (0, 3, 9, 27),
+               A vertex at (0, 4, 16, 64),
+               A vertex at (0, 5, 25, 125),
+               A vertex at (0, 6, 36, 216),
+               A vertex at (0, 7, 49, 343)]])
+
+        TESTS::
+
+            sage: P = polytopes.cross_polytope(5)
+            sage: P.is_prism()
+            False
+
+            sage: P = polytopes.permutahedron(4).prism()
+            sage: P.is_prism()
+            True
+
+            sage: P = Polyhedron(vertices=[[0,1], [1,0]], rays=[[1,1]])
+            sage: P.is_prism()
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: polyhedron has to be compact
+
+        ALGORITHM:
+
+        See :meth:`Polyhedron_base.is_bipyramid`.
+        """
+        if not self.is_compact():
+            raise NotImplementedError("polyhedron has to be compact")
+
+        from sage.misc.functional import is_odd
+        n_verts = self.n_vertices()
+        n_facets = self.n_facets()
+        if is_odd(n_verts):
+            if certificate:
+                return (False, None)
+            return False
+
+        IM = self.incidence_matrix()
+        if self.n_equations():
+            # Remove equations from the incidence matrix,
+            # such that this is the vertex-facet incidences matrix.
+            I1 = IM.transpose()
+            I2 = I1[[i for i in range(self.n_Hrepresentation())
+                     if not self.Hrepresentation()[i].is_equation()]]
+            IM = I2.transpose()
+
+        verts_incidences = [set(row.nonzero_positions()) for row in IM.rows()]
+        facets_incidences = dict()
+        for j in range(n_facets):
+            F_j = set(IM.column(j).nonzero_positions())
+            if len(F_j) == n_verts/2:
+                facets_incidences[j] = F_j
+
+        # Find two vertices ``facet1`` and ``facet2`` such that one of them
+        # contains exactly half of the vertices, and the other one contains
+        # exactly the other half.
+        from itertools import combinations
+        for index1, index2 in combinations(facets_incidences, 2):
+            facet1_incidences = facets_incidences[index1]
+            facet2_incidences = facets_incidences[index2]
+            facet1and2 = facet1_incidences.union(facet2_incidences)
+            if len(facet1and2) == n_verts:
+                # We have found two candidates for base faces.
+                # Remove from each vertex ``index1`` resp. ``index2``.
+                test_verts = set(frozenset(vert_inc.difference({index1, index2}))
+                                  for vert_inc in verts_incidences)
+                if len(test_verts) == n_verts/2:
+                    # For each vertex containing `index1` there is
+                    # another one contained in `index2`
+                    # and vice versa.
+                    # Other than `index1` and `index2` both are contained in
+                    # exactly the same facets.
+                    if certificate:
+                        V = self.vertices()
+                        facet1_vertices = [V[i] for i in facet1_incidences]
+                        facet2_vertices = [V[i] for i in facet2_incidences]
+                        return (True, [facet1_vertices, facet2_vertices])
+                    return True
+
+        if certificate:
+            return (False, None)
+        return False
+
     def hyperplane_arrangement(self):
         """
         Return the hyperplane arrangement defined by the equations and