From 6a15f929d3be603a16f3a0310abfe33fdb4d7615 Mon Sep 17 00:00:00 2001 From: Jonathan Kliem Date: Tue, 16 Nov 2021 09:57:52 +0100 Subject: [PATCH] move combinatorics to Polyhedron_base3 --- src/sage/geometry/polyhedron/base.py | 1920 +------------------------ src/sage/geometry/polyhedron/base0.py | 63 + src/sage/geometry/polyhedron/base3.py | 1902 ++++++++++++++++++++++++ 3 files changed, 1980 insertions(+), 1905 deletions(-) create mode 100644 src/sage/geometry/polyhedron/base3.py diff --git a/src/sage/geometry/polyhedron/base.py b/src/sage/geometry/polyhedron/base.py index 29f9a543f23..946f9daee8f 100644 --- a/src/sage/geometry/polyhedron/base.py +++ b/src/sage/geometry/polyhedron/base.py @@ -47,7 +47,7 @@ from sage.geometry.convex_set import AffineHullProjectionData from .constructor import Polyhedron -from .base2 import Polyhedron_base2 +from .base3 import Polyhedron_base3 ######################################################################### # Notes if you want to implement your own backend: @@ -92,7 +92,7 @@ def is_Polyhedron(X): ######################################################################### -class Polyhedron_base(Polyhedron_base2): +class Polyhedron_base(Polyhedron_base3): """ Base class for Polyhedron objects @@ -165,27 +165,6 @@ class Polyhedron_base(Polyhedron_base2): sage: TestSuite(P).run() """ - def _init_empty_polyhedron(self): - """ - Initializes an empty polyhedron. - - TESTS:: - - sage: Polyhedron().vertex_adjacency_matrix() # indirect doctest - [] - sage: Polyhedron().facet_adjacency_matrix() - [0] - """ - Polyhedron_base2._init_empty_polyhedron(self) - - V_matrix = matrix(ZZ, 0, 0, 0) - V_matrix.set_immutable() - self.vertex_adjacency_matrix.set_cache(V_matrix) - - H_matrix = matrix(ZZ, 1, 1, 0) - H_matrix.set_immutable() - self.facet_adjacency_matrix.set_cache(H_matrix) - def _sage_input_(self, sib, coerced): """ Return Sage command to reconstruct ``self``. @@ -219,48 +198,6 @@ def _sage_input_(self, sib, coerced): kwds['lines'] = [sib(tuple(l)) for l in self.lines()] return sib.name('Polyhedron')(**kwds) - def _facet_adjacency_matrix(self): - """ - Compute the facet adjacency matrix in case it has not been - computed during initialization. - - EXAMPLES:: - - sage: p = Polyhedron(vertices=[(0,0),(1,0),(0,1)]) - sage: p._facet_adjacency_matrix() - [0 1 1] - [1 0 1] - [1 1 0] - - Checks that :trac:`22455` is fixed:: - - sage: s = polytopes.simplex(2) - sage: s._facet_adjacency_matrix() - [0 1 1] - [1 0 1] - [1 1 0] - - """ - # TODO: This implementation computes the whole face lattice, - # which is much more information than necessary. - M = matrix(ZZ, self.n_facets(), self.n_facets(), 0) - codim = self.ambient_dim()-self.dim() - - def set_adjacent(h1, h2): - if h1 is h2: - return - i = h1.index() - codim - j = h2.index() - codim - M[i, j] = 1 - M[j, i] = 1 - - for face in self.faces(self.dim()-2): - Hrep = face.ambient_Hrepresentation() - assert(len(Hrep) == codim+2) - set_adjacent(Hrep[-2], Hrep[-1]) - M.set_immutable() - return M - def _delete(self): """ Delete this polyhedron. @@ -1079,223 +1016,6 @@ def to_linear_program(self, solver=None, return_variable=False, base_ring=None): else: return p - @cached_method - def vertices_matrix(self, base_ring=None): - """ - Return the coordinates of the vertices as the columns of a matrix. - - INPUT: - - - ``base_ring`` -- A ring or ``None`` (default). The base ring - of the returned matrix. If not specified, the base ring of - the polyhedron is used. - - OUTPUT: - - A matrix over ``base_ring`` whose columns are the coordinates - of the vertices. A ``TypeError`` is raised if the coordinates - cannot be converted to ``base_ring``. - - .. WARNING:: - - If the polyhedron has lines, return the coordinates of the vertices - of the ``Vrepresentation``. However, the represented polyhedron - has no 0-dimensional faces (i.e. vertices):: - - sage: P = Polyhedron(rays=[[1,0,0]],lines=[[0,1,0]]) - sage: P.vertices_matrix() - [0] - [0] - [0] - sage: P.faces(0) - () - - EXAMPLES:: - - sage: triangle = Polyhedron(vertices=[[1,0],[0,1],[1,1]]) - sage: triangle.vertices_matrix() - [0 1 1] - [1 0 1] - sage: (triangle/2).vertices_matrix() - [ 0 1/2 1/2] - [1/2 0 1/2] - sage: (triangle/2).vertices_matrix(ZZ) - Traceback (most recent call last): - ... - TypeError: no conversion of this rational to integer - - TESTS: - - Check that :trac:`28828` is fixed:: - - sage: P.vertices_matrix().is_immutable() - True - """ - if base_ring is None: - base_ring = self.base_ring() - m = matrix(base_ring, self.ambient_dim(), self.n_vertices()) - for i, v in enumerate(self.vertices()): - for j in range(self.ambient_dim()): - m[j, i] = v[j] - m.set_immutable() - return m - - def bounded_edges(self): - """ - Return the bounded edges (excluding rays and lines). - - OUTPUT: - - A generator for pairs of vertices, one pair per edge. - - EXAMPLES:: - - sage: p = Polyhedron(vertices=[[1,0],[0,1]], rays=[[1,0],[0,1]]) - sage: [ e for e in p.bounded_edges() ] - [(A vertex at (0, 1), A vertex at (1, 0))] - sage: for e in p.bounded_edges(): print(e) - (A vertex at (0, 1), A vertex at (1, 0)) - """ - obj = self.Vrepresentation() - for i in range(len(obj)): - if not obj[i].is_vertex(): - continue - for j in range(i+1, len(obj)): - if not obj[j].is_vertex(): - continue - if self.vertex_adjacency_matrix()[i, j] == 0: - continue - yield (obj[i], obj[j]) - - @cached_method - def vertex_adjacency_matrix(self): - """ - Return the binary matrix of vertex adjacencies. - - EXAMPLES:: - - sage: polytopes.simplex(4).vertex_adjacency_matrix() - [0 1 1 1 1] - [1 0 1 1 1] - [1 1 0 1 1] - [1 1 1 0 1] - [1 1 1 1 0] - - The rows and columns of the vertex adjacency matrix correspond - to the :meth:`Vrepresentation` objects: vertices, rays, and - lines. The `(i,j)` matrix entry equals `1` if the `i`-th and - `j`-th V-representation object are adjacent. - - Two vertices are adjacent if they are the endpoints of an - edge, that is, a one-dimensional face. For unbounded polyhedra - this clearly needs to be generalized and we define two - V-representation objects (see - :mod:`sage.geometry.polyhedron.constructor`) to be adjacent if - they together generate a one-face. There are three possible - combinations: - - * Two vertices can bound a finite-length edge. - - * A vertex and a ray can generate a half-infinite edge - starting at the vertex and with the direction given by the - ray. - - * A vertex and a line can generate an infinite edge. The - position of the vertex on the line is arbitrary in this - case, only its transverse position matters. The direction of - the edge is given by the line generator. - - For example, take the half-plane:: - - sage: half_plane = Polyhedron(ieqs=[(0,1,0)]) - sage: half_plane.Hrepresentation() - (An inequality (1, 0) x + 0 >= 0,) - - Its (non-unique) V-representation consists of a vertex, a ray, - and a line. The only edge is spanned by the vertex and the - line generator, so they are adjacent:: - - sage: half_plane.Vrepresentation() - (A line in the direction (0, 1), A ray in the direction (1, 0), A vertex at (0, 0)) - sage: half_plane.vertex_adjacency_matrix() - [0 0 1] - [0 0 0] - [1 0 0] - - In one dimension higher, that is for a half-space in 3 - dimensions, there is no one-dimensional face. Hence nothing is - adjacent:: - - sage: Polyhedron(ieqs=[(0,1,0,0)]).vertex_adjacency_matrix() - [0 0 0 0] - [0 0 0 0] - [0 0 0 0] - [0 0 0 0] - - EXAMPLES: - - In a bounded polygon, every vertex has precisely two adjacent ones:: - - sage: P = Polyhedron(vertices=[(0, 1), (1, 0), (3, 0), (4, 1)]) - sage: for v in P.Vrep_generator(): - ....: print("{} {}".format(P.adjacency_matrix().row(v.index()), v)) - (0, 1, 0, 1) A vertex at (0, 1) - (1, 0, 1, 0) A vertex at (1, 0) - (0, 1, 0, 1) A vertex at (3, 0) - (1, 0, 1, 0) A vertex at (4, 1) - - If the V-representation of the polygon contains vertices and - one ray, then each V-representation object is adjacent to two - V-representation objects:: - - sage: P = Polyhedron(vertices=[(0, 1), (1, 0), (3, 0), (4, 1)], - ....: rays=[(0,1)]) - sage: for v in P.Vrep_generator(): - ....: print("{} {}".format(P.adjacency_matrix().row(v.index()), v)) - (0, 1, 0, 0, 1) A ray in the direction (0, 1) - (1, 0, 1, 0, 0) A vertex at (0, 1) - (0, 1, 0, 1, 0) A vertex at (1, 0) - (0, 0, 1, 0, 1) A vertex at (3, 0) - (1, 0, 0, 1, 0) A vertex at (4, 1) - - If the V-representation of the polygon contains vertices and - two distinct rays, then each vertex is adjacent to two - V-representation objects (which can now be vertices or - rays). The two rays are not adjacent to each other:: - - sage: P = Polyhedron(vertices=[(0, 1), (1, 0), (3, 0), (4, 1)], - ....: rays=[(0,1), (1,1)]) - sage: for v in P.Vrep_generator(): - ....: print("{} {}".format(P.adjacency_matrix().row(v.index()), v)) - (0, 1, 0, 0, 0) A ray in the direction (0, 1) - (1, 0, 1, 0, 0) A vertex at (0, 1) - (0, 1, 0, 0, 1) A vertex at (1, 0) - (0, 0, 0, 0, 1) A ray in the direction (1, 1) - (0, 0, 1, 1, 0) A vertex at (3, 0) - - The vertex adjacency matrix has base ring integers. This way one can express various - counting questions:: - - sage: P = polytopes.cube() - sage: Q = P.stack(P.faces(2)[0]) - sage: M = Q.vertex_adjacency_matrix() - sage: sum(M) - (4, 4, 3, 3, 4, 4, 4, 3, 3) - sage: G = Q.vertex_graph() - sage: G.degree() - [4, 4, 3, 3, 4, 4, 4, 3, 3] - - TESTS: - - Check that :trac:`28828` is fixed:: - - sage: P.adjacency_matrix().is_immutable() - True - """ - return self.combinatorial_polyhedron().vertex_adjacency_matrix() - - adjacency_matrix = vertex_adjacency_matrix - def boundary_complex(self): """ Return the simplicial complex given by the boundary faces of ``self``, @@ -1351,280 +1071,6 @@ def boundary_complex(self): else: raise NotImplementedError("this function is only implemented for simplicial polytopes") - @cached_method - def facet_adjacency_matrix(self): - """ - Return the adjacency matrix for the facets and hyperplanes. - - EXAMPLES:: - - sage: s4 = polytopes.simplex(4, project=True) - sage: s4.facet_adjacency_matrix() - [0 1 1 1 1] - [1 0 1 1 1] - [1 1 0 1 1] - [1 1 1 0 1] - [1 1 1 1 0] - - The facet adjacency matrix has base ring integers. This way one can express various - counting questions:: - - sage: P = polytopes.cube() - sage: Q = P.stack(P.faces(2)[0]) - sage: M = Q.facet_adjacency_matrix() - sage: sum(M) - (4, 4, 4, 4, 3, 3, 3, 3, 4) - - TESTS: - - Check that :trac:`28828` is fixed:: - - sage: s4.facet_adjacency_matrix().is_immutable() - True - """ - return self._facet_adjacency_matrix() - - @cached_method - def incidence_matrix(self): - """ - Return the incidence matrix. - - .. NOTE:: - - The columns correspond to inequalities/equations in the - order :meth:`Hrepresentation`, the rows correspond to - vertices/rays/lines in the order - :meth:`Vrepresentation`. - - .. SEEALSO:: - - :meth:`slack_matrix`. - - EXAMPLES:: - - sage: p = polytopes.cuboctahedron() - sage: p.incidence_matrix() - [0 0 1 1 0 1 0 0 0 0 1 0 0 0] - [0 0 0 1 0 0 1 0 1 0 1 0 0 0] - [0 0 1 1 1 0 0 1 0 0 0 0 0 0] - [1 0 0 1 1 0 1 0 0 0 0 0 0 0] - [0 0 0 0 0 1 0 0 1 1 1 0 0 0] - [0 0 1 0 0 1 0 1 0 0 0 1 0 0] - [1 0 0 0 0 0 1 0 1 0 0 0 1 0] - [1 0 0 0 1 0 0 1 0 0 0 0 0 1] - [0 1 0 0 0 1 0 0 0 1 0 1 0 0] - [0 1 0 0 0 0 0 0 1 1 0 0 1 0] - [0 1 0 0 0 0 0 1 0 0 0 1 0 1] - [1 1 0 0 0 0 0 0 0 0 0 0 1 1] - sage: v = p.Vrepresentation(0) - sage: v - A vertex at (-1, -1, 0) - sage: h = p.Hrepresentation(2) - sage: h - An inequality (1, 1, -1) x + 2 >= 0 - sage: h.eval(v) # evaluation (1, 1, -1) * (-1/2, -1/2, 0) + 1 - 0 - sage: h*v # same as h.eval(v) - 0 - sage: p.incidence_matrix() [0,2] # this entry is (v,h) - 1 - sage: h.contains(v) - True - sage: p.incidence_matrix() [2,0] # note: not symmetric - 0 - - The incidence matrix depends on the ambient dimension:: - - sage: simplex = polytopes.simplex(); simplex - A 3-dimensional polyhedron in ZZ^4 defined as the convex hull of 4 vertices - sage: simplex.incidence_matrix() - [1 1 1 1 0] - [1 1 1 0 1] - [1 1 0 1 1] - [1 0 1 1 1] - sage: simplex = simplex.affine_hull_projection(); simplex - A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices - sage: simplex.incidence_matrix() - [1 1 1 0] - [1 1 0 1] - [1 0 1 1] - [0 1 1 1] - - An incidence matrix does not determine a unique - polyhedron:: - - sage: P = Polyhedron(vertices=[[0,1],[1,1],[1,0]]) - sage: P.incidence_matrix() - [1 1 0] - [1 0 1] - [0 1 1] - - sage: Q = Polyhedron(vertices=[[0,1], [1,0]], rays=[[1,1]]) - sage: Q.incidence_matrix() - [1 1 0] - [1 0 1] - [0 1 1] - - - An example of two polyhedra with isomorphic face lattices - but different incidence matrices:: - - sage: Q.incidence_matrix() - [1 1 0] - [1 0 1] - [0 1 1] - - sage: R = Polyhedron(vertices=[[0,1], [1,0]], rays=[[1,3/2], [3/2,1]]) - sage: R.incidence_matrix() - [1 1 0] - [1 0 1] - [0 1 0] - [0 0 1] - - The incidence matrix has base ring integers. This way one can express various - counting questions:: - - sage: P = polytopes.twenty_four_cell() - sage: M = P.incidence_matrix() - sage: sum(sum(x) for x in M) == P.flag_f_vector(0,3) - True - - TESTS: - - Check that :trac:`28828` is fixed:: - - sage: R.incidence_matrix().is_immutable() - True - - Test that this method works for inexact base ring - (``cdd`` sets the cache already):: - - sage: P = polytopes.dodecahedron(exact=False) - sage: M = P.incidence_matrix.cache - sage: P.incidence_matrix.clear_cache() - sage: M == P.incidence_matrix() - True - """ - if self.base_ring() in (ZZ, QQ): - # Much faster for integers or rationals. - incidence_matrix = self.slack_matrix().zero_pattern_matrix(ZZ) - incidence_matrix.set_immutable() - return incidence_matrix - - incidence_matrix = matrix(ZZ, self.n_Vrepresentation(), - self.n_Hrepresentation(), 0) - - Vvectors_vertices = tuple((v.vector(), v.index()) - for v in self.Vrep_generator() - if v.is_vertex()) - Vvectors_rays_lines = tuple((v.vector(), v.index()) - for v in self.Vrep_generator() - if not v.is_vertex()) - - # Determine ``is_zero`` to save lots of time. - if self.base_ring().is_exact(): - def is_zero(x): - return not x - else: - is_zero = self._is_zero - - for H in self.Hrep_generator(): - Hconst = H.b() - Hvec = H.A() - Hindex = H.index() - for Vvec, Vindex in Vvectors_vertices: - if is_zero(Hvec*Vvec + Hconst): - incidence_matrix[Vindex, Hindex] = 1 - - # A ray or line is considered incident with a hyperplane, - # if it is orthogonal to the normal vector of the hyperplane. - for Vvec, Vindex in Vvectors_rays_lines: - if is_zero(Hvec*Vvec): - incidence_matrix[Vindex, Hindex] = 1 - - incidence_matrix.set_immutable() - return incidence_matrix - - @cached_method - def slack_matrix(self): - r""" - Return the slack matrix. - - The entries correspond to the evaluation of the Hrepresentation - elements on the Vrepresentation elements. - - .. NOTE:: - - The columns correspond to inequalities/equations in the - order :meth:`Hrepresentation`, the rows correspond to - vertices/rays/lines in the order - :meth:`Vrepresentation`. - - .. SEEALSO:: - - :meth:`incidence_matrix`. - - EXAMPLES:: - - sage: P = polytopes.cube() - sage: P.slack_matrix() - [0 2 2 2 0 0] - [0 0 2 2 0 2] - [0 0 0 2 2 2] - [0 2 0 2 2 0] - [2 2 0 0 2 0] - [2 2 2 0 0 0] - [2 0 2 0 0 2] - [2 0 0 0 2 2] - - sage: P = polytopes.cube(intervals='zero_one') - sage: P.slack_matrix() - [0 1 1 1 0 0] - [0 0 1 1 0 1] - [0 0 0 1 1 1] - [0 1 0 1 1 0] - [1 1 0 0 1 0] - [1 1 1 0 0 0] - [1 0 1 0 0 1] - [1 0 0 0 1 1] - - sage: P = polytopes.dodecahedron().faces(2)[0].as_polyhedron() - sage: P.slack_matrix() - [1/2*sqrt5 - 1/2 0 0 1 1/2*sqrt5 - 1/2 0] - [ 0 0 1/2*sqrt5 - 1/2 1/2*sqrt5 - 1/2 1 0] - [ 0 1/2*sqrt5 - 1/2 1 0 1/2*sqrt5 - 1/2 0] - [ 1 1/2*sqrt5 - 1/2 0 1/2*sqrt5 - 1/2 0 0] - [1/2*sqrt5 - 1/2 1 1/2*sqrt5 - 1/2 0 0 0] - - sage: P = Polyhedron(rays=[[1, 0], [0, 1]]) - sage: P.slack_matrix() - [0 0] - [0 1] - [1 0] - - TESTS:: - - sage: Polyhedron().slack_matrix() - [] - sage: Polyhedron(base_ring=QuadraticField(2)).slack_matrix().base_ring() - Number Field in a with defining polynomial x^2 - 2 with a = 1.41... - """ - if not self.n_Vrepresentation() or not self.n_Hrepresentation(): - slack_matrix = matrix(self.base_ring(), self.n_Vrepresentation(), - self.n_Hrepresentation(), 0) - else: - Vrep_matrix = matrix(self.base_ring(), self.Vrepresentation()) - Hrep_matrix = matrix(self.base_ring(), self.Hrepresentation()) - - # Getting homogeneous coordinates of the Vrepresentation. - hom_helper = matrix(self.base_ring(), [1 if v.is_vertex() else 0 for v in self.Vrepresentation()]) - hom_Vrep = hom_helper.stack(Vrep_matrix.transpose()) - - slack_matrix = (Hrep_matrix * hom_Vrep).transpose() - - slack_matrix.set_immutable() - return slack_matrix - @cached_method def center(self): """ @@ -1748,51 +1194,6 @@ def centroid(self, engine='auto', **kwds): pass return centroid - def a_maximal_chain(self): - r""" - Return a maximal chain of the face lattice in increasing order. - - EXAMPLES:: - - sage: P = polytopes.cube() - sage: P.a_maximal_chain() - [A -1-dimensional face of a Polyhedron in ZZ^3, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices] - sage: P = polytopes.cube() - sage: chain = P.a_maximal_chain(); chain - [A -1-dimensional face of a Polyhedron in ZZ^3, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices] - sage: [face.ambient_V_indices() for face in chain] - [(), (5,), (0, 5), (0, 3, 4, 5), (0, 1, 2, 3, 4, 5, 6, 7)] - - TESTS:: - - Check output for the empty polyhedron:: - - sage: P = Polyhedron() - sage: P.a_maximal_chain() - [A -1-dimensional face of a Polyhedron in ZZ^0] - """ - comb_chain = self.combinatorial_polyhedron().a_maximal_chain() - - from sage.geometry.polyhedron.face import combinatorial_face_to_polyhedral_face - empty_face = self.faces(-1)[0] - universe = self.faces(self.dim())[0] - - if self.dim() == -1: - return [empty_face] - - return [empty_face] + \ - [combinatorial_face_to_polyhedral_face(self, face) - for face in comb_chain] + \ - [universe] - @cached_method def radius_square(self): """ @@ -2009,471 +1410,6 @@ def is_inscribed(self, certificate=False): else: return is_inscribed - @cached_method - def combinatorial_polyhedron(self): - r""" - Return the combinatorial type of ``self``. - - See :class:`sage.geometry.polyhedron.combinatorial_polyhedron.base.CombinatorialPolyhedron`. - - EXAMPLES:: - - sage: polytopes.cube().combinatorial_polyhedron() - A 3-dimensional combinatorial polyhedron with 6 facets - - sage: polytopes.cyclic_polytope(4,10).combinatorial_polyhedron() - A 4-dimensional combinatorial polyhedron with 35 facets - - sage: Polyhedron(rays=[[0,1], [1,0]]).combinatorial_polyhedron() - A 2-dimensional combinatorial polyhedron with 2 facets - """ - from sage.geometry.polyhedron.combinatorial_polyhedron.base import CombinatorialPolyhedron - return CombinatorialPolyhedron(self) - - def _test_combinatorial_polyhedron(self, tester=None, **options): - """ - Run test suite of combinatorial polyhedron. - - TESTS:: - - sage: polytopes.cross_polytope(3)._test_combinatorial_polyhedron() - """ - from sage.misc.sage_unittest import TestSuite - - tester = self._tester(tester=tester, **options) - tester.info("\n Running the test suite of self.combinatorial_polyhedron()") - TestSuite(self.combinatorial_polyhedron()).run(verbose=tester._verbose, - prefix=tester._prefix+" ") - tester.info(tester._prefix+" ", newline = False) - - def simplicity(self): - r""" - Return the largest integer `k` such that the polytope is `k`-simple. - - A polytope `P` is `k`-simple, if every `(d-1-k)`-face - is contained in exactly `k+1` facets of `P` for `1 \leq k \leq d-1`. - Equivalently it is `k`-simple if the polar/dual polytope is `k`-simplicial. - If `self` is a simplex, it returns its dimension. - - EXAMPLES:: - - sage: polytopes.hypersimplex(4,2).simplicity() - 1 - sage: polytopes.hypersimplex(5,2).simplicity() - 2 - sage: polytopes.hypersimplex(6,2).simplicity() - 3 - sage: polytopes.simplex(3).simplicity() - 3 - sage: polytopes.simplex(1).simplicity() - 1 - - The method is not implemented for unbounded polyhedra:: - - sage: p = Polyhedron(vertices=[(0,0)],rays=[(1,0),(0,1)]) - sage: p.simplicity() - Traceback (most recent call last): - ... - NotImplementedError: this function is implemented for polytopes only - """ - if not(self.is_compact()): - raise NotImplementedError("this function is implemented for polytopes only") - return self.combinatorial_polyhedron().simplicity() - - def is_simple(self): - """ - Test for simplicity of a polytope. - - See :wikipedia:`Simple_polytope` - - EXAMPLES:: - - sage: p = Polyhedron([[0,0,0],[1,0,0],[0,1,0],[0,0,1]]) - sage: p.is_simple() - True - sage: p = Polyhedron([[0,0,0],[4,4,0],[4,0,0],[0,4,0],[2,2,2]]) - sage: p.is_simple() - False - """ - if not self.is_compact(): - return False - return self.combinatorial_polyhedron().is_simple() - - def simpliciality(self): - r""" - Return the largest integer `k` such that the polytope is `k`-simplicial. - - A polytope is `k`-simplicial, if every `k`-face is a simplex. - If `self` is a simplex, returns its dimension. - - EXAMPLES:: - - sage: polytopes.cyclic_polytope(10,4).simpliciality() - 3 - sage: polytopes.hypersimplex(5,2).simpliciality() - 2 - sage: polytopes.cross_polytope(4).simpliciality() - 3 - sage: polytopes.simplex(3).simpliciality() - 3 - sage: polytopes.simplex(1).simpliciality() - 1 - - The method is not implemented for unbounded polyhedra:: - - sage: p = Polyhedron(vertices=[(0,0)],rays=[(1,0),(0,1)]) - sage: p.simpliciality() - Traceback (most recent call last): - ... - NotImplementedError: this function is implemented for polytopes only - """ - if not(self.is_compact()): - raise NotImplementedError("this function is implemented for polytopes only") - return self.combinatorial_polyhedron().simpliciality() - - def is_simplicial(self): - """ - Tests if the polytope is simplicial - - A polytope is simplicial if every facet is a simplex. - - See :wikipedia:`Simplicial_polytope` - - EXAMPLES:: - - sage: p = polytopes.hypercube(3) - sage: p.is_simplicial() - False - sage: q = polytopes.simplex(5, project=True) - sage: q.is_simplicial() - True - sage: p = Polyhedron([[0,0,0],[1,0,0],[0,1,0],[0,0,1]]) - sage: p.is_simplicial() - True - sage: q = Polyhedron([[1,1,1],[-1,1,1],[1,-1,1],[-1,-1,1],[1,1,-1]]) - sage: q.is_simplicial() - False - sage: P = polytopes.simplex(); P - A 3-dimensional polyhedron in ZZ^4 defined as the convex hull of 4 vertices - sage: P.is_simplicial() - True - - The method is not implemented for unbounded polyhedra:: - - sage: p = Polyhedron(vertices=[(0,0)],rays=[(1,0),(0,1)]) - sage: p.is_simplicial() - Traceback (most recent call last): - ... - NotImplementedError: this function is implemented for polytopes only - """ - if not(self.is_compact()): - raise NotImplementedError("this function is implemented for polytopes only") - return self.combinatorial_polyhedron().is_simplicial() - - 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 (1, 0, 0, 0)) - sage: egyptian_pyramid = polytopes.regular_polygon(4).pyramid() # optional - sage.rings.number_field - sage: egyptian_pyramid.is_pyramid() # optional - sage.rings.number_field - True - 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") - - return self.combinatorial_polyhedron().is_pyramid(certificate) - - 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, 6, 36, 216), - A vertex at (1, 0, 0, 0), - A vertex at (1, 7, 49, 343), - A vertex at (1, 5, 25, 125), - A vertex at (1, 1, 1, 1), - A vertex at (1, 2, 4, 8), - A vertex at (1, 4, 16, 64), - A vertex at (1, 3, 9, 27)], - [A vertex at (0, 3, 9, 27), - A vertex at (0, 6, 36, 216), - A vertex at (0, 0, 0, 0), - A vertex at (0, 7, 49, 343), - A vertex at (0, 5, 25, 125), - A vertex at (0, 1, 1, 1), - A vertex at (0, 2, 4, 8), - A vertex at (0, 4, 16, 64)]]) - - 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 @@ -4630,55 +3566,27 @@ def lawrence_polytope(self): sage: V = P.vertices_list() sage: i = 0 sage: for v in V: - ....: v = v + i*[0] - ....: P = P.lawrence_extension(v) - ....: i = i + 1 - sage: P == L - True - - REFERENCES: - - For more information, see Section 6.6 of [Zie2007]_. - """ - from sage.matrix.constructor import block_matrix - - if not self.is_compact(): - raise NotImplementedError("self must be a polytope") - - V = self.vertices_matrix().transpose() - n = self.n_vertices() - I_n = matrix.identity(n) - lambda_V = block_matrix([[V, I_n], [V, 2*I_n]]) - parent = self.parent().change_ring(self.base_ring(), ambient_dim=self.ambient_dim() + n) - return parent.element_class(parent, [lambda_V, [], []], None) - - 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. - - EXAMPLES:: - - sage: P = polytopes.hypersimplex(5,2) - sage: L = P.lawrence_polytope() - sage: L.is_lattice_polytope() - True - sage: egyptian_pyramid = polytopes.regular_polygon(4).pyramid() - sage: egyptian_pyramid.is_lawrence_polytope() + ....: v = v + i*[0] + ....: P = P.lawrence_extension(v) + ....: i = i + 1 + sage: P == L True - sage: polytopes.octahedron().is_lawrence_polytope() - False REFERENCES: - For more information, see [BaSt1990]_. + For more information, see Section 6.6 of [Zie2007]_. """ + from sage.matrix.constructor import block_matrix + if not self.is_compact(): raise NotImplementedError("self must be a polytope") - return self.combinatorial_polyhedron().is_lawrence_polytope() + V = self.vertices_matrix().transpose() + n = self.n_vertices() + I_n = matrix.identity(n) + lambda_V = block_matrix([[V, I_n], [V, 2*I_n]]) + parent = self.parent().change_ring(self.base_ring(), ambient_dim=self.ambient_dim() + n) + return parent.element_class(parent, [lambda_V, [], []], None) def _test_lawrence(self, tester=None, **options): """ @@ -5084,736 +3992,6 @@ def index_to_polyhedron_face(n): return D.relabel(index_to_polyhedron_face, inplace=False, immutable=True) - def face_generator(self, face_dimension=None, dual=None): - r""" - Return an iterator over the faces of given dimension. - - If dimension is not specified return an iterator over all faces. - - INPUT: - - - ``face_dimension`` -- integer (default ``None``), - yield only faces of this dimension if specified - - ``dual`` -- boolean (default ``None``); - if ``True``, generate the faces using the vertices; - if ``False``, generate the faces using the facets; - if ``None``, pick automatically - - OUTPUT: - - A :class:`~sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator.FaceIterator_geom`. - This class iterates over faces as - :class:`~sage.geometry.polyhedron.face.PolyhedronFace`. See - :mod:`~sage.geometry.polyhedron.face` for details. The order - is random but fixed. - - EXAMPLES:: - - sage: P = polytopes.cube() - sage: it = P.face_generator() - sage: it - Iterator over the faces of a 3-dimensional polyhedron in ZZ^3 - sage: list(it) - [A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices, - A -1-dimensional face of a Polyhedron in ZZ^3, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices] - - sage: P = polytopes.hypercube(4) - sage: list(P.face_generator(2))[:4] - [A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices, - A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices] - - If a polytope has more facets than vertices, the dual mode is chosen:: - - sage: P = polytopes.cross_polytope(3) - sage: list(P.face_generator()) - [A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 6 vertices, - A -1-dimensional face of a Polyhedron in ZZ^3, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices] - - The face iterator can also be slightly modified. - In non-dual mode we can skip subfaces of the current (proper) face:: - - sage: P = polytopes.cube() - sage: it = P.face_generator(dual=False) - sage: _ = next(it), next(it) - sage: face = next(it) - sage: face.ambient_H_indices() - (5,) - sage: it.ignore_subfaces() - sage: face = next(it) - sage: face.ambient_H_indices() - (4,) - sage: it.ignore_subfaces() - sage: [face.ambient_H_indices() for face in it] - [(3,), - (2,), - (1,), - (0,), - (2, 3), - (1, 3), - (1, 2, 3), - (1, 2), - (0, 2), - (0, 1, 2), - (0, 1)] - - In dual mode we can skip supfaces of the current (proper) face:: - - sage: P = polytopes.cube() - sage: it = P.face_generator(dual=True) - sage: _ = next(it), next(it) - sage: face = next(it) - sage: face.ambient_V_indices() - (7,) - sage: it.ignore_supfaces() - sage: next(it) - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex - sage: face = next(it) - sage: face.ambient_V_indices() - (5,) - sage: it.ignore_supfaces() - sage: [face.ambient_V_indices() for face in it] - [(4,), - (3,), - (2,), - (1,), - (0,), - (1, 6), - (3, 4), - (2, 3), - (0, 3), - (0, 1, 2, 3), - (1, 2), - (0, 1)] - - In non-dual mode, we cannot skip supfaces:: - - sage: it = P.face_generator(dual=False) - sage: _ = next(it), next(it) - sage: next(it) - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices - sage: it.ignore_supfaces() - Traceback (most recent call last): - ... - ValueError: only possible when in dual mode - - In dual mode, we cannot skip subfaces:: - - sage: it = P.face_generator(dual=True) - sage: _ = next(it), next(it) - sage: next(it) - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex - sage: it.ignore_subfaces() - Traceback (most recent call last): - ... - ValueError: only possible when not in dual mode - - We can only skip sub-/supfaces of proper faces:: - - sage: it = P.face_generator(dual=False) - sage: next(it) - A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices - sage: it.ignore_subfaces() - Traceback (most recent call last): - ... - ValueError: iterator not set to a face yet - - .. SEEALSO:: - - :class:`~sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator.FaceIterator_geom`. - - ALGORITHM: - - See :class:`~sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator.FaceIterator`. - - TESTS:: - - sage: P = polytopes.simplex() - sage: list(P.face_generator(-2)) - [] - sage: list(P.face_generator(-1)) - [A -1-dimensional face of a Polyhedron in ZZ^4] - sage: list(P.face_generator(3)) - [A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices] - - sage: list(Polyhedron().face_generator()) - [A -1-dimensional face of a Polyhedron in ZZ^0] - - Check that :trac:`29155` is fixed:: - - sage: P = polytopes.permutahedron(3) - sage: [f] = P.face_generator(2) - sage: f.ambient_Hrepresentation() - (An equation (1, 1, 1) x - 6 == 0,) - """ - from sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator import FaceIterator_geom - return FaceIterator_geom(self, output_dimension=face_dimension, dual=dual) - - def faces(self, face_dimension): - """ - Return the faces of given dimension - - INPUT: - - - ``face_dimension`` -- integer. The dimension of the faces - whose representation will be returned. - - OUTPUT: - - A tuple of - :class:`~sage.geometry.polyhedron.face.PolyhedronFace`. See - :mod:`~sage.geometry.polyhedron.face` for details. The order - is random but fixed. - - .. SEEALSO:: - - :meth:`face_generator`, - :meth:`facet`. - - EXAMPLES: - - Here we find the vertex and face indices of the eight three-dimensional - facets of the four-dimensional hypercube:: - - sage: p = polytopes.hypercube(4) - sage: list(f.ambient_V_indices() for f in p.faces(3)) - [(0, 5, 6, 7, 8, 9, 14, 15), - (1, 4, 5, 6, 10, 13, 14, 15), - (1, 2, 6, 7, 8, 10, 11, 15), - (8, 9, 10, 11, 12, 13, 14, 15), - (0, 3, 4, 5, 9, 12, 13, 14), - (0, 2, 3, 7, 8, 9, 11, 12), - (1, 2, 3, 4, 10, 11, 12, 13), - (0, 1, 2, 3, 4, 5, 6, 7)] - - sage: face = p.faces(3)[3] - sage: face.ambient_Hrepresentation() - (An inequality (1, 0, 0, 0) x + 1 >= 0,) - sage: face.vertices() - (A vertex at (-1, -1, 1, -1), - A vertex at (-1, -1, 1, 1), - A vertex at (-1, 1, -1, -1), - A vertex at (-1, 1, 1, -1), - A vertex at (-1, 1, 1, 1), - A vertex at (-1, 1, -1, 1), - A vertex at (-1, -1, -1, 1), - A vertex at (-1, -1, -1, -1)) - - You can use the - :meth:`~sage.geometry.polyhedron.representation.PolyhedronRepresentation.index` - method to enumerate vertices and inequalities:: - - sage: def get_idx(rep): return rep.index() - sage: [get_idx(_) for _ in face.ambient_Hrepresentation()] - [4] - sage: [get_idx(_) for _ in face.ambient_Vrepresentation()] - [8, 9, 10, 11, 12, 13, 14, 15] - - sage: [ ([get_idx(_) for _ in face.ambient_Vrepresentation()], - ....: [get_idx(_) for _ in face.ambient_Hrepresentation()]) - ....: for face in p.faces(3) ] - [([0, 5, 6, 7, 8, 9, 14, 15], [7]), - ([1, 4, 5, 6, 10, 13, 14, 15], [6]), - ([1, 2, 6, 7, 8, 10, 11, 15], [5]), - ([8, 9, 10, 11, 12, 13, 14, 15], [4]), - ([0, 3, 4, 5, 9, 12, 13, 14], [3]), - ([0, 2, 3, 7, 8, 9, 11, 12], [2]), - ([1, 2, 3, 4, 10, 11, 12, 13], [1]), - ([0, 1, 2, 3, 4, 5, 6, 7], [0])] - - TESTS:: - - sage: pr = Polyhedron(rays = [[1,0,0],[-1,0,0],[0,1,0]], vertices = [[-1,-1,-1]], lines=[(0,0,1)]) - sage: pr.faces(4) - () - sage: pr.faces(3)[0].ambient_V_indices() - (0, 1, 2, 3) - sage: pr.facets()[0].ambient_V_indices() - (0, 1, 2) - sage: pr.faces(1) - () - sage: pr.faces(0) - () - sage: pr.faces(-1) - (A -1-dimensional face of a Polyhedron in QQ^3,) - """ - return tuple(self.face_generator(face_dimension)) - - def facets(self): - r""" - Return the facets of the polyhedron. - - Facets are the maximal nontrivial faces of polyhedra. - The empty face and the polyhedron itself are trivial. - - A facet of a `d`-dimensional polyhedron is a face of dimension - `d-1`. For `d \neq 0` the converse is true as well. - - OUTPUT: - - A tuple of - :class:`~sage.geometry.polyhedron.face.PolyhedronFace`. See - :mod:`~sage.geometry.polyhedron.face` for details. The order - is random but fixed. - - .. SEEALSO:: :meth:`facets` - - EXAMPLES: - - Here we find the eight three-dimensional facets of the - four-dimensional hypercube:: - - sage: p = polytopes.hypercube(4) - sage: p.facets() - (A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices) - - This is the same result as explicitly finding the - three-dimensional faces:: - - sage: dim = p.dimension() - sage: p.faces(dim-1) - (A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, - A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices) - - The ``0``-dimensional polyhedron does not have facets:: - - sage: P = Polyhedron([[0]]) - sage: P.facets() - () - """ - if self.dimension() == 0: - return () - return self.faces(self.dimension()-1) - - def join_of_Vrep(self, *Vrepresentatives): - r""" - Return the smallest face that contains ``Vrepresentatives``. - - INPUT: - - - ``Vrepresentatives`` -- vertices/rays/lines of ``self`` or indices of such - - OUTPUT: a :class:`~sage.geometry.polyhedron.face.PolyhedronFace` - - .. NOTE:: - - In the case of unbounded polyhedra, the join of rays etc. may not be well-defined. - - EXAMPLES:: - - sage: P = polytopes.permutahedron(5) - sage: P.join_of_Vrep(1) - A 0-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 1 vertex - sage: P.join_of_Vrep() - A -1-dimensional face of a Polyhedron in ZZ^5 - sage: P.join_of_Vrep(0,12,13).ambient_V_indices() - (0, 12, 13, 68) - - The input is flexible:: - - sage: P.join_of_Vrep(2, P.vertices()[3], P.Vrepresentation(4)) - A 2-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 6 vertices - - :: - - sage: P = polytopes.cube() - sage: a, b = P.faces(0)[:2] - sage: P.join_of_Vrep(a, b) - A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices - - In the case of an unbounded polyhedron, the join may not be well-defined:: - - sage: P = Polyhedron(vertices=[[1,0], [0,1]], rays=[[1,1]]) - sage: P.join_of_Vrep(0) - A 0-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex - sage: P.join_of_Vrep(0,1) - A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 2 vertices - sage: P.join_of_Vrep(0,2) - A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray - sage: P.join_of_Vrep(1,2) - A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray - sage: P.join_of_Vrep(2) - Traceback (most recent call last): - ... - ValueError: the join is not well-defined - - The ``Vrepresentatives`` must be of ``self``:: - - sage: P = polytopes.cube(backend='ppl') - sage: Q = polytopes.cube(backend='field') - sage: v = P.vertices()[0] - sage: P.join_of_Vrep(v) - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex - sage: Q.join_of_Vrep(v) - Traceback (most recent call last): - ... - ValueError: not a Vrepresentative of ``self`` - sage: f = P.faces(0)[0] - sage: P.join_of_Vrep(v) - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex - sage: Q.join_of_Vrep(v) - Traceback (most recent call last): - ... - ValueError: not a Vrepresentative of ``self`` - - TESTS: - - ``least_common_superface_of_Vrep`` is an alias:: - - sage: P.least_common_superface_of_Vrep(v) - A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex - sage: P.least_common_superface_of_Vrep == P.join_of_Vrep - True - - Error message for invalid input:: - - sage: P.join_of_Vrep('foo') - Traceback (most recent call last): - ... - ValueError: foo is not a Vrepresentative - """ - from sage.geometry.polyhedron.representation import Vrepresentation - from sage.geometry.polyhedron.face import PolyhedronFace - - new_indices = [0]*len(Vrepresentatives) - for i, v in enumerate(Vrepresentatives): - if isinstance(v, PolyhedronFace) and v.dim() == 0: - if v.polyhedron() is not self: - raise ValueError("not a Vrepresentative of ``self``") - new_indices[i] = v.ambient_V_indices()[0] - elif v in ZZ: - new_indices[i] = v - elif isinstance(v, Vrepresentation): - if v.polyhedron() is not self: - raise ValueError("not a Vrepresentative of ``self``") - new_indices[i] = v.index() - else: - raise ValueError("{} is not a Vrepresentative".format(v)) - - return self.face_generator().join_of_Vrep(*new_indices) - - least_common_superface_of_Vrep = join_of_Vrep - - def meet_of_Hrep(self, *Hrepresentatives): - r""" - Return the largest face that is contained in ``Hrepresentatives``. - - INPUT: - - - ``Hrepresentatives`` -- facets or indices of Hrepresentatives; - the indices are assumed to be the indices of the Hrepresentation - - OUTPUT: a :class:`~sage.geometry.polyhedron.face.PolyhedronFace` - - EXAMPLES:: - - sage: P = polytopes.permutahedron(5) - sage: P.meet_of_Hrep() - A 4-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 120 vertices - sage: P.meet_of_Hrep(1) - A 3-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 24 vertices - sage: P.meet_of_Hrep(4) - A 3-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 12 vertices - sage: P.meet_of_Hrep(1,3,7) - A 1-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 2 vertices - sage: P.meet_of_Hrep(1,3,7).ambient_H_indices() - (0, 1, 3, 7) - - The indices are the indices of the Hrepresentation. - ``0`` corresponds to an equation and is ignored:: - - sage: P.meet_of_Hrep(0) - A 4-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 120 vertices - - The input is flexible:: - - sage: P.meet_of_Hrep(P.facets()[-1], P.inequalities()[2], 7) - A 1-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 2 vertices - - The ``Hrepresentatives`` must belong to ``self``:: - - sage: P = polytopes.cube(backend='ppl') - sage: Q = polytopes.cube(backend='field') - sage: f = P.facets()[0] - sage: P.meet_of_Hrep(f) - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices - sage: Q.meet_of_Hrep(f) - Traceback (most recent call last): - ... - ValueError: not a facet of ``self`` - sage: f = P.inequalities()[0] - sage: P.meet_of_Hrep(f) - A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices - sage: Q.meet_of_Hrep(f) - Traceback (most recent call last): - ... - ValueError: not a facet of ``self`` - - TESTS: - - Equations are not considered by the combinatorial polyhedron. - We check that the index corresponds to the Hrepresentation index:: - - sage: P = polytopes.permutahedron(3, backend='field') - sage: P.Hrepresentation() - (An inequality (0, 0, 1) x - 1 >= 0, - An inequality (0, 1, 0) x - 1 >= 0, - An inequality (0, 1, 1) x - 3 >= 0, - An inequality (1, 0, 0) x - 1 >= 0, - An inequality (1, 0, 1) x - 3 >= 0, - An inequality (1, 1, 0) x - 3 >= 0, - An equation (1, 1, 1) x - 6 == 0) - sage: P.meet_of_Hrep(0).ambient_Hrepresentation() - (An inequality (0, 0, 1) x - 1 >= 0, An equation (1, 1, 1) x - 6 == 0) - - sage: P = polytopes.permutahedron(3, backend='ppl') - sage: P.Hrepresentation() - (An equation (1, 1, 1) x - 6 == 0, - An inequality (1, 1, 0) x - 3 >= 0, - An inequality (-1, -1, 0) x + 5 >= 0, - An inequality (0, 1, 0) x - 1 >= 0, - An inequality (-1, 0, 0) x + 3 >= 0, - An inequality (1, 0, 0) x - 1 >= 0, - An inequality (0, -1, 0) x + 3 >= 0) - sage: P.meet_of_Hrep(1).ambient_Hrepresentation() - (An equation (1, 1, 1) x - 6 == 0, An inequality (1, 1, 0) x - 3 >= 0) - - ``greatest_common_subface_of_Hrep`` is an alias:: - - sage: P.greatest_common_subface_of_Hrep(1).ambient_Hrepresentation() - (An equation (1, 1, 1) x - 6 == 0, An inequality (1, 1, 0) x - 3 >= 0) - sage: P.greatest_common_subface_of_Hrep == P.meet_of_Hrep - True - - Error message for invalid input:: - - sage: P.meet_of_Hrep('foo') - Traceback (most recent call last): - ... - ValueError: foo is not a Hrepresentative - """ - from sage.geometry.polyhedron.representation import Inequality, Equation - from sage.geometry.polyhedron.face import PolyhedronFace - - # Equations are ignored by combinatorial polyhedron for indexing. - offset = 0 - if self.n_equations() and self.Hrepresentation(0).is_equation(): - offset = self.n_equations() - - new_indices = [] - for i, facet in enumerate(Hrepresentatives): - if isinstance(facet, PolyhedronFace) and facet.dim() + 1 == self.dim(): - if facet.polyhedron() is not self: - raise ValueError("not a facet of ``self``") - H_indices = facet.ambient_H_indices() - facet = H_indices[0] if H_indices[0] >= offset else H_indices[-1] - - if facet in ZZ and facet >= offset: - # Note that ``CombinatorialPolyhedron`` ignores indices of equations - # and has equations last. - new_indices.append(facet - offset) - elif isinstance(facet, Inequality): - if facet.polyhedron() is not self: - raise ValueError("not a facet of ``self``") - new_indices.append(facet.index() - offset) - elif isinstance(facet, Equation): - # Ignore equations. - continue - elif facet in ZZ and 0 <= facet < offset: - # Ignore equations. - continue - else: - raise ValueError("{} is not a Hrepresentative".format(facet)) - - return self.face_generator().meet_of_Hrep(*new_indices) - - greatest_common_subface_of_Hrep = meet_of_Hrep - - def _test_combinatorial_face_as_combinatorial_polyhedron(self, tester=None, **options): - """ - Run tests on obtaining the combinatorial face as combinatorial polyhedron. - - TESTS:: - - sage: polytopes.cross_polytope(3)._test_combinatorial_face_as_combinatorial_polyhedron() - """ - if not self.is_compact(): - return - if self.dim() > 7 or self.n_vertices() > 100 or self.n_facets() > 100: - # Avoid very long tests. - return - if self.dim() < 1: - # Avoid trivial cases. - return - - from sage.misc.prandom import random - - if tester is None: - tester = self._tester(**options) - - it = self.face_generator() - _ = next(it), next(it) # get rid of non_proper faces - C1 = self.combinatorial_polyhedron() - it1 = C1.face_iter() - C2 = C1.dual() - it2 = C2.face_iter(dual=not it1.dual) - - for f in it: - f1 = next(it1) - f2 = next(it2) - if random() < 0.95: - # Only test a random 5 percent of the faces. - continue - - P = f.as_polyhedron() - D1 = f1.as_combinatorial_polyhedron() - D2 = f2.as_combinatorial_polyhedron(quotient=True).dual() - D1._test_bitsets(tester, **options) - D2._test_bitsets(tester, **options) - try: - import sage.graphs.graph - except ImportError: - pass - else: - tester.assertTrue(P.combinatorial_polyhedron().vertex_facet_graph().is_isomorphic(D1.vertex_facet_graph())) - tester.assertTrue(P.combinatorial_polyhedron().vertex_facet_graph().is_isomorphic(D2.vertex_facet_graph())) - - @cached_method(do_pickle=True) - def f_vector(self, num_threads=None, parallelization_depth=None): - r""" - Return the f-vector. - - INPUT: - - - ``num_threads`` -- integer (optional); specify the number of threads; - otherwise determined by :func:`~sage.parallel.ncpus.ncpus` - - - ``parallelization_depth`` -- integer (optional); specify - how deep in the lattice the parallelization is done - - OUTPUT: - - Return a vector whose `i`-th entry is the number of - `i-2`-dimensional faces of the polytope. - - .. NOTE:: - - The ``vertices`` as given by :meth:`Polyhedron_base.vertices` - do not need to correspond to `0`-dimensional faces. If a polyhedron - contains `k` lines they correspond to `k`-dimensional faces. - See example below - - EXAMPLES:: - - sage: p = Polyhedron(vertices=[[1, 2, 3], [1, 3, 2], - ....: [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1], [0, 0, 0]]) - sage: p.f_vector() - (1, 7, 12, 7, 1) - - sage: polytopes.cyclic_polytope(4,10).f_vector() - (1, 10, 45, 70, 35, 1) - - sage: polytopes.hypercube(5).f_vector() - (1, 32, 80, 80, 40, 10, 1) - - Polyhedra with lines do not have `0`-faces:: - - sage: Polyhedron(ieqs=[[1,-1,0,0],[1,1,0,0]]).f_vector() - (1, 0, 0, 2, 1) - - However, the method :meth:`Polyhedron_base.vertices` returns - two points that belong to the ``Vrepresentation``:: - - sage: P = Polyhedron(ieqs=[[1,-1,0],[1,1,0]]) - sage: P.vertices() - (A vertex at (1, 0), A vertex at (-1, 0)) - sage: P.f_vector() - (1, 0, 2, 1) - - TESTS: - - Check that :trac:`28828` is fixed:: - - sage: P.f_vector().is_immutable() - True - - The cache of the f-vector is being pickled:: - - sage: P = polytopes.cube() - sage: P.f_vector() - (1, 8, 12, 6, 1) - sage: Q = loads(dumps(P)) - sage: Q.f_vector.is_in_cache() - True - """ - return self.combinatorial_polyhedron().f_vector(num_threads, parallelization_depth) - def flag_f_vector(self, *args): r""" Return the flag f-vector. @@ -7435,24 +5613,6 @@ def _integrate_latte_(self, polynomial, **kwds): polynomial, cdd=True, **kwds) - def is_simplex(self): - r""" - Return whether the polyhedron is a simplex. - - A simplex is a bounded polyhedron with `d+1` vertices, where - `d` is the dimension. - - EXAMPLES:: - - sage: Polyhedron([(0,0,0), (1,0,0), (0,1,0)]).is_simplex() - True - sage: polytopes.simplex(3).is_simplex() - True - sage: polytopes.hypercube(3).is_simplex() - False - """ - return self.is_compact() and (self.dim()+1 == self.n_vertices()) - def neighborliness(self): r""" Return the largest ``k``, such that the polyhedron is ``k``-neighborly. @@ -7497,56 +5657,6 @@ def neighborliness(self): """ return self.combinatorial_polyhedron().neighborliness() - def is_neighborly(self, k=None): - r""" - Return whether the polyhedron is neighborly. - - If the input ``k`` is provided, then return whether the polyhedron is ``k``-neighborly - - A polyhedron is neighborly if every set of `n` vertices forms a face - for `n` up to floor of half the dimension of the polyhedron. - It is `k`-neighborly if this is true for `n` up to `k`. - - INPUT: - - - ``k`` -- the dimension up to which to check if every set of ``k`` - vertices forms a face. If no ``k`` is provided, check up to floor - of half the dimension of the polyhedron. - - OUTPUT: - - - ``True`` if every set of up to ``k`` vertices forms a face, - - ``False`` otherwise - - .. SEEALSO:: - - :meth:`neighborliness` - - EXAMPLES:: - - sage: cube = polytopes.hypercube(3) - sage: cube.is_neighborly() - True - sage: cube = polytopes.hypercube(4) - sage: cube.is_neighborly() - False - - Cyclic polytopes are neighborly:: - - sage: all(polytopes.cyclic_polytope(i, i + 1 + j).is_neighborly() for i in range(5) for j in range(3)) - True - - The neighborliness of a polyhedron equals floor of dimension half - (or larger in case of a simplex) if and only if the polyhedron - is neighborly:: - - sage: testpolys = [polytopes.cube(), polytopes.cyclic_polytope(6, 9), polytopes.simplex(6)] - sage: [(P.neighborliness()>=floor(P.dim()/2)) == P.is_neighborly() for P in testpolys] - [True, True, True] - - """ - return self.combinatorial_polyhedron().is_neighborly() - @cached_method def bounding_box(self, integral=False, integral_hull=False): r""" diff --git a/src/sage/geometry/polyhedron/base0.py b/src/sage/geometry/polyhedron/base0.py index c9bbd20dd81..a16a4016124 100644 --- a/src/sage/geometry/polyhedron/base0.py +++ b/src/sage/geometry/polyhedron/base0.py @@ -1025,6 +1025,69 @@ def vertices(self): """ return tuple(self.vertex_generator()) + @cached_method + def vertices_matrix(self, base_ring=None): + """ + Return the coordinates of the vertices as the columns of a matrix. + + INPUT: + + - ``base_ring`` -- A ring or ``None`` (default). The base ring + of the returned matrix. If not specified, the base ring of + the polyhedron is used. + + OUTPUT: + + A matrix over ``base_ring`` whose columns are the coordinates + of the vertices. A ``TypeError`` is raised if the coordinates + cannot be converted to ``base_ring``. + + .. WARNING:: + + If the polyhedron has lines, return the coordinates of the vertices + of the ``Vrepresentation``. However, the represented polyhedron + has no 0-dimensional faces (i.e. vertices):: + + sage: P = Polyhedron(rays=[[1,0,0]],lines=[[0,1,0]]) + sage: P.vertices_matrix() + [0] + [0] + [0] + sage: P.faces(0) + () + + EXAMPLES:: + + sage: triangle = Polyhedron(vertices=[[1,0],[0,1],[1,1]]) + sage: triangle.vertices_matrix() + [0 1 1] + [1 0 1] + sage: (triangle/2).vertices_matrix() + [ 0 1/2 1/2] + [1/2 0 1/2] + sage: (triangle/2).vertices_matrix(ZZ) + Traceback (most recent call last): + ... + TypeError: no conversion of this rational to integer + + TESTS: + + Check that :trac:`28828` is fixed:: + + sage: P.vertices_matrix().is_immutable() + True + """ + from sage.matrix.constructor import matrix + + if base_ring is None: + base_ring = self.base_ring() + m = matrix(base_ring, self.ambient_dim(), self.n_vertices()) + for i, v in enumerate(self.vertices()): + for j in range(self.ambient_dim()): + m[j, i] = v[j] + m.set_immutable() + return m + def ray_generator(self): """ Return a generator for the rays of the polyhedron. diff --git a/src/sage/geometry/polyhedron/base3.py b/src/sage/geometry/polyhedron/base3.py new file mode 100644 index 00000000000..a8d3b8bbc9b --- /dev/null +++ b/src/sage/geometry/polyhedron/base3.py @@ -0,0 +1,1902 @@ +r""" +Base class for polyhedra, part 3 + +Define methods related to the combinatorics of a polyhedron +excluding methods relying on :mod:`sage.graphs`. +""" + +# **************************************************************************** +# Copyright (C) 2008-2012 Marshall Hampton +# Copyright (C) 2011-2015 Volker Braun +# Copyright (C) 2012-2018 Frederic Chapoton +# Copyright (C) 2013 Andrey Novoseltsev +# Copyright (C) 2014-2017 Moritz Firsching +# Copyright (C) 2014-2019 Thierry Monteil +# Copyright (C) 2015 Nathann Cohen +# Copyright (C) 2015-2017 Jeroen Demeyer +# Copyright (C) 2015-2017 Vincent Delecroix +# Copyright (C) 2015-2018 Dima Pasechnik +# Copyright (C) 2015-2020 Jean-Philippe Labbe +# Copyright (C) 2015-2021 Matthias Koeppe +# Copyright (C) 2016-2019 Daniel Krenn +# Copyright (C) 2017 Marcelo Forets +# Copyright (C) 2017-2018 Mark Bell +# Copyright (C) 2019 Julian Ritter +# Copyright (C) 2019-2020 Laith Rastanawi +# Copyright (C) 2019-2020 Sophia Elia +# Copyright (C) 2019-2021 Jonathan Kliem +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 2 of the License, or +# (at your option) any later version. +# https://www.gnu.org/licenses/ +# **************************************************************************** + +from sage.misc.cachefunc import cached_method +from sage.matrix.constructor import matrix +from sage.rings.integer_ring import ZZ +from sage.rings.rational_field import QQ +from .base2 import Polyhedron_base2 + +class Polyhedron_base3(Polyhedron_base2): + """ + Methods related to lattice points. + + See :class:`sage.geometry.polyhedron.base.Polyhedron_base`. + + TESTS:: + + sage: from sage.geometry.polyhedron.base3 import Polyhedron_base3 + sage: P = polytopes.cube() + sage: Polyhedron_base3.is_simple(P) + True + sage: Polyhedron_base3.is_simplicial(P) + False + sage: Polyhedron_base3.is_prism(P) + True + sage: Polyhedron_base3.is_pyramid(P) + False + sage: Polyhedron_base3.combinatorial_polyhedron.f(P) + A 3-dimensional combinatorial polyhedron with 6 facets + sage: Polyhedron_base3.facets(P) + (A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices) + sage: Polyhedron_base3.f_vector.f(P) + (1, 8, 12, 6, 1) + sage: next(Polyhedron_base3.face_generator(P)) + A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices + """ + + def _init_empty_polyhedron(self): + """ + Initializes an empty polyhedron. + + TESTS:: + + sage: Polyhedron().vertex_adjacency_matrix() # indirect doctest + [] + sage: Polyhedron().facet_adjacency_matrix() + [0] + """ + Polyhedron_base2._init_empty_polyhedron(self) + + V_matrix = matrix(ZZ, 0, 0, 0) + V_matrix.set_immutable() + self.vertex_adjacency_matrix.set_cache(V_matrix) + + H_matrix = matrix(ZZ, 1, 1, 0) + H_matrix.set_immutable() + self.facet_adjacency_matrix.set_cache(H_matrix) + + @cached_method + def slack_matrix(self): + r""" + Return the slack matrix. + + The entries correspond to the evaluation of the Hrepresentation + elements on the Vrepresentation elements. + + .. NOTE:: + + The columns correspond to inequalities/equations in the + order :meth:`Hrepresentation`, the rows correspond to + vertices/rays/lines in the order + :meth:`Vrepresentation`. + + .. SEEALSO:: + + :meth:`incidence_matrix`. + + EXAMPLES:: + + sage: P = polytopes.cube() + sage: P.slack_matrix() + [0 2 2 2 0 0] + [0 0 2 2 0 2] + [0 0 0 2 2 2] + [0 2 0 2 2 0] + [2 2 0 0 2 0] + [2 2 2 0 0 0] + [2 0 2 0 0 2] + [2 0 0 0 2 2] + + sage: P = polytopes.cube(intervals='zero_one') + sage: P.slack_matrix() + [0 1 1 1 0 0] + [0 0 1 1 0 1] + [0 0 0 1 1 1] + [0 1 0 1 1 0] + [1 1 0 0 1 0] + [1 1 1 0 0 0] + [1 0 1 0 0 1] + [1 0 0 0 1 1] + + sage: P = polytopes.dodecahedron().faces(2)[0].as_polyhedron() + sage: P.slack_matrix() + [1/2*sqrt5 - 1/2 0 0 1 1/2*sqrt5 - 1/2 0] + [ 0 0 1/2*sqrt5 - 1/2 1/2*sqrt5 - 1/2 1 0] + [ 0 1/2*sqrt5 - 1/2 1 0 1/2*sqrt5 - 1/2 0] + [ 1 1/2*sqrt5 - 1/2 0 1/2*sqrt5 - 1/2 0 0] + [1/2*sqrt5 - 1/2 1 1/2*sqrt5 - 1/2 0 0 0] + + sage: P = Polyhedron(rays=[[1, 0], [0, 1]]) + sage: P.slack_matrix() + [0 0] + [0 1] + [1 0] + + TESTS:: + + sage: Polyhedron().slack_matrix() + [] + sage: Polyhedron(base_ring=QuadraticField(2)).slack_matrix().base_ring() + Number Field in a with defining polynomial x^2 - 2 with a = 1.41... + """ + if not self.n_Vrepresentation() or not self.n_Hrepresentation(): + slack_matrix = matrix(self.base_ring(), self.n_Vrepresentation(), + self.n_Hrepresentation(), 0) + else: + Vrep_matrix = matrix(self.base_ring(), self.Vrepresentation()) + Hrep_matrix = matrix(self.base_ring(), self.Hrepresentation()) + + # Getting homogeneous coordinates of the Vrepresentation. + hom_helper = matrix(self.base_ring(), [1 if v.is_vertex() else 0 for v in self.Vrepresentation()]) + hom_Vrep = hom_helper.stack(Vrep_matrix.transpose()) + + slack_matrix = (Hrep_matrix * hom_Vrep).transpose() + + slack_matrix.set_immutable() + return slack_matrix + + @cached_method + def incidence_matrix(self): + """ + Return the incidence matrix. + + .. NOTE:: + + The columns correspond to inequalities/equations in the + order :meth:`Hrepresentation`, the rows correspond to + vertices/rays/lines in the order + :meth:`Vrepresentation`. + + .. SEEALSO:: + + :meth:`slack_matrix`. + + EXAMPLES:: + + sage: p = polytopes.cuboctahedron() + sage: p.incidence_matrix() + [0 0 1 1 0 1 0 0 0 0 1 0 0 0] + [0 0 0 1 0 0 1 0 1 0 1 0 0 0] + [0 0 1 1 1 0 0 1 0 0 0 0 0 0] + [1 0 0 1 1 0 1 0 0 0 0 0 0 0] + [0 0 0 0 0 1 0 0 1 1 1 0 0 0] + [0 0 1 0 0 1 0 1 0 0 0 1 0 0] + [1 0 0 0 0 0 1 0 1 0 0 0 1 0] + [1 0 0 0 1 0 0 1 0 0 0 0 0 1] + [0 1 0 0 0 1 0 0 0 1 0 1 0 0] + [0 1 0 0 0 0 0 0 1 1 0 0 1 0] + [0 1 0 0 0 0 0 1 0 0 0 1 0 1] + [1 1 0 0 0 0 0 0 0 0 0 0 1 1] + sage: v = p.Vrepresentation(0) + sage: v + A vertex at (-1, -1, 0) + sage: h = p.Hrepresentation(2) + sage: h + An inequality (1, 1, -1) x + 2 >= 0 + sage: h.eval(v) # evaluation (1, 1, -1) * (-1/2, -1/2, 0) + 1 + 0 + sage: h*v # same as h.eval(v) + 0 + sage: p.incidence_matrix() [0,2] # this entry is (v,h) + 1 + sage: h.contains(v) + True + sage: p.incidence_matrix() [2,0] # note: not symmetric + 0 + + The incidence matrix depends on the ambient dimension:: + + sage: simplex = polytopes.simplex(); simplex + A 3-dimensional polyhedron in ZZ^4 defined as the convex hull of 4 vertices + sage: simplex.incidence_matrix() + [1 1 1 1 0] + [1 1 1 0 1] + [1 1 0 1 1] + [1 0 1 1 1] + sage: simplex = simplex.affine_hull_projection(); simplex + A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 4 vertices + sage: simplex.incidence_matrix() + [1 1 1 0] + [1 1 0 1] + [1 0 1 1] + [0 1 1 1] + + An incidence matrix does not determine a unique + polyhedron:: + + sage: P = Polyhedron(vertices=[[0,1],[1,1],[1,0]]) + sage: P.incidence_matrix() + [1 1 0] + [1 0 1] + [0 1 1] + + sage: Q = Polyhedron(vertices=[[0,1], [1,0]], rays=[[1,1]]) + sage: Q.incidence_matrix() + [1 1 0] + [1 0 1] + [0 1 1] + + + An example of two polyhedra with isomorphic face lattices + but different incidence matrices:: + + sage: Q.incidence_matrix() + [1 1 0] + [1 0 1] + [0 1 1] + + sage: R = Polyhedron(vertices=[[0,1], [1,0]], rays=[[1,3/2], [3/2,1]]) + sage: R.incidence_matrix() + [1 1 0] + [1 0 1] + [0 1 0] + [0 0 1] + + The incidence matrix has base ring integers. This way one can express various + counting questions:: + + sage: P = polytopes.twenty_four_cell() + sage: M = P.incidence_matrix() + sage: sum(sum(x) for x in M) == P.flag_f_vector(0,3) + True + + TESTS: + + Check that :trac:`28828` is fixed:: + + sage: R.incidence_matrix().is_immutable() + True + + Test that this method works for inexact base ring + (``cdd`` sets the cache already):: + + sage: P = polytopes.dodecahedron(exact=False) + sage: M = P.incidence_matrix.cache + sage: P.incidence_matrix.clear_cache() + sage: M == P.incidence_matrix() + True + """ + if self.base_ring() in (ZZ, QQ): + # Much faster for integers or rationals. + incidence_matrix = self.slack_matrix().zero_pattern_matrix(ZZ) + incidence_matrix.set_immutable() + return incidence_matrix + + incidence_matrix = matrix(ZZ, self.n_Vrepresentation(), + self.n_Hrepresentation(), 0) + + Vvectors_vertices = tuple((v.vector(), v.index()) + for v in self.Vrep_generator() + if v.is_vertex()) + Vvectors_rays_lines = tuple((v.vector(), v.index()) + for v in self.Vrep_generator() + if not v.is_vertex()) + + # Determine ``is_zero`` to save lots of time. + if self.base_ring().is_exact(): + def is_zero(x): + return not x + else: + is_zero = self._is_zero + + for H in self.Hrep_generator(): + Hconst = H.b() + Hvec = H.A() + Hindex = H.index() + for Vvec, Vindex in Vvectors_vertices: + if is_zero(Hvec*Vvec + Hconst): + incidence_matrix[Vindex, Hindex] = 1 + + # A ray or line is considered incident with a hyperplane, + # if it is orthogonal to the normal vector of the hyperplane. + for Vvec, Vindex in Vvectors_rays_lines: + if is_zero(Hvec*Vvec): + incidence_matrix[Vindex, Hindex] = 1 + + incidence_matrix.set_immutable() + return incidence_matrix + + @cached_method + def combinatorial_polyhedron(self): + r""" + Return the combinatorial type of ``self``. + + See :class:`sage.geometry.polyhedron.combinatorial_polyhedron.base.CombinatorialPolyhedron`. + + EXAMPLES:: + + sage: polytopes.cube().combinatorial_polyhedron() + A 3-dimensional combinatorial polyhedron with 6 facets + + sage: polytopes.cyclic_polytope(4,10).combinatorial_polyhedron() + A 4-dimensional combinatorial polyhedron with 35 facets + + sage: Polyhedron(rays=[[0,1], [1,0]]).combinatorial_polyhedron() + A 2-dimensional combinatorial polyhedron with 2 facets + """ + from sage.geometry.polyhedron.combinatorial_polyhedron.base import CombinatorialPolyhedron + return CombinatorialPolyhedron(self) + + def _test_combinatorial_polyhedron(self, tester=None, **options): + """ + Run test suite of combinatorial polyhedron. + + TESTS:: + + sage: polytopes.cross_polytope(3)._test_combinatorial_polyhedron() + """ + from sage.misc.sage_unittest import TestSuite + + tester = self._tester(tester=tester, **options) + tester.info("\n Running the test suite of self.combinatorial_polyhedron()") + TestSuite(self.combinatorial_polyhedron()).run(verbose=tester._verbose, + prefix=tester._prefix+" ") + tester.info(tester._prefix+" ", newline = False) + + def face_generator(self, face_dimension=None, dual=None): + r""" + Return an iterator over the faces of given dimension. + + If dimension is not specified return an iterator over all faces. + + INPUT: + + - ``face_dimension`` -- integer (default ``None``), + yield only faces of this dimension if specified + - ``dual`` -- boolean (default ``None``); + if ``True``, generate the faces using the vertices; + if ``False``, generate the faces using the facets; + if ``None``, pick automatically + + OUTPUT: + + A :class:`~sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator.FaceIterator_geom`. + This class iterates over faces as + :class:`~sage.geometry.polyhedron.face.PolyhedronFace`. See + :mod:`~sage.geometry.polyhedron.face` for details. The order + is random but fixed. + + EXAMPLES:: + + sage: P = polytopes.cube() + sage: it = P.face_generator() + sage: it + Iterator over the faces of a 3-dimensional polyhedron in ZZ^3 + sage: list(it) + [A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices, + A -1-dimensional face of a Polyhedron in ZZ^3, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices] + + sage: P = polytopes.hypercube(4) + sage: list(P.face_generator(2))[:4] + [A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices, + A 2-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices] + + If a polytope has more facets than vertices, the dual mode is chosen:: + + sage: P = polytopes.cross_polytope(3) + sage: list(P.face_generator()) + [A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 6 vertices, + A -1-dimensional face of a Polyhedron in ZZ^3, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 3 vertices, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices] + + The face iterator can also be slightly modified. + In non-dual mode we can skip subfaces of the current (proper) face:: + + sage: P = polytopes.cube() + sage: it = P.face_generator(dual=False) + sage: _ = next(it), next(it) + sage: face = next(it) + sage: face.ambient_H_indices() + (5,) + sage: it.ignore_subfaces() + sage: face = next(it) + sage: face.ambient_H_indices() + (4,) + sage: it.ignore_subfaces() + sage: [face.ambient_H_indices() for face in it] + [(3,), + (2,), + (1,), + (0,), + (2, 3), + (1, 3), + (1, 2, 3), + (1, 2), + (0, 2), + (0, 1, 2), + (0, 1)] + + In dual mode we can skip supfaces of the current (proper) face:: + + sage: P = polytopes.cube() + sage: it = P.face_generator(dual=True) + sage: _ = next(it), next(it) + sage: face = next(it) + sage: face.ambient_V_indices() + (7,) + sage: it.ignore_supfaces() + sage: next(it) + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex + sage: face = next(it) + sage: face.ambient_V_indices() + (5,) + sage: it.ignore_supfaces() + sage: [face.ambient_V_indices() for face in it] + [(4,), + (3,), + (2,), + (1,), + (0,), + (1, 6), + (3, 4), + (2, 3), + (0, 3), + (0, 1, 2, 3), + (1, 2), + (0, 1)] + + In non-dual mode, we cannot skip supfaces:: + + sage: it = P.face_generator(dual=False) + sage: _ = next(it), next(it) + sage: next(it) + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices + sage: it.ignore_supfaces() + Traceback (most recent call last): + ... + ValueError: only possible when in dual mode + + In dual mode, we cannot skip subfaces:: + + sage: it = P.face_generator(dual=True) + sage: _ = next(it), next(it) + sage: next(it) + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex + sage: it.ignore_subfaces() + Traceback (most recent call last): + ... + ValueError: only possible when not in dual mode + + We can only skip sub-/supfaces of proper faces:: + + sage: it = P.face_generator(dual=False) + sage: next(it) + A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices + sage: it.ignore_subfaces() + Traceback (most recent call last): + ... + ValueError: iterator not set to a face yet + + .. SEEALSO:: + + :class:`~sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator.FaceIterator_geom`. + + ALGORITHM: + + See :class:`~sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator.FaceIterator`. + + TESTS:: + + sage: P = polytopes.simplex() + sage: list(P.face_generator(-2)) + [] + sage: list(P.face_generator(-1)) + [A -1-dimensional face of a Polyhedron in ZZ^4] + sage: list(P.face_generator(3)) + [A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 4 vertices] + + sage: list(Polyhedron().face_generator()) + [A -1-dimensional face of a Polyhedron in ZZ^0] + + Check that :trac:`29155` is fixed:: + + sage: P = polytopes.permutahedron(3) + sage: [f] = P.face_generator(2) + sage: f.ambient_Hrepresentation() + (An equation (1, 1, 1) x - 6 == 0,) + """ + from sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator import FaceIterator_geom + return FaceIterator_geom(self, output_dimension=face_dimension, dual=dual) + + def faces(self, face_dimension): + """ + Return the faces of given dimension + + INPUT: + + - ``face_dimension`` -- integer. The dimension of the faces + whose representation will be returned. + + OUTPUT: + + A tuple of + :class:`~sage.geometry.polyhedron.face.PolyhedronFace`. See + :mod:`~sage.geometry.polyhedron.face` for details. The order + is random but fixed. + + .. SEEALSO:: + + :meth:`face_generator`, + :meth:`facet`. + + EXAMPLES: + + Here we find the vertex and face indices of the eight three-dimensional + facets of the four-dimensional hypercube:: + + sage: p = polytopes.hypercube(4) + sage: list(f.ambient_V_indices() for f in p.faces(3)) + [(0, 5, 6, 7, 8, 9, 14, 15), + (1, 4, 5, 6, 10, 13, 14, 15), + (1, 2, 6, 7, 8, 10, 11, 15), + (8, 9, 10, 11, 12, 13, 14, 15), + (0, 3, 4, 5, 9, 12, 13, 14), + (0, 2, 3, 7, 8, 9, 11, 12), + (1, 2, 3, 4, 10, 11, 12, 13), + (0, 1, 2, 3, 4, 5, 6, 7)] + + sage: face = p.faces(3)[3] + sage: face.ambient_Hrepresentation() + (An inequality (1, 0, 0, 0) x + 1 >= 0,) + sage: face.vertices() + (A vertex at (-1, -1, 1, -1), + A vertex at (-1, -1, 1, 1), + A vertex at (-1, 1, -1, -1), + A vertex at (-1, 1, 1, -1), + A vertex at (-1, 1, 1, 1), + A vertex at (-1, 1, -1, 1), + A vertex at (-1, -1, -1, 1), + A vertex at (-1, -1, -1, -1)) + + You can use the + :meth:`~sage.geometry.polyhedron.representation.PolyhedronRepresentation.index` + method to enumerate vertices and inequalities:: + + sage: def get_idx(rep): return rep.index() + sage: [get_idx(_) for _ in face.ambient_Hrepresentation()] + [4] + sage: [get_idx(_) for _ in face.ambient_Vrepresentation()] + [8, 9, 10, 11, 12, 13, 14, 15] + + sage: [ ([get_idx(_) for _ in face.ambient_Vrepresentation()], + ....: [get_idx(_) for _ in face.ambient_Hrepresentation()]) + ....: for face in p.faces(3) ] + [([0, 5, 6, 7, 8, 9, 14, 15], [7]), + ([1, 4, 5, 6, 10, 13, 14, 15], [6]), + ([1, 2, 6, 7, 8, 10, 11, 15], [5]), + ([8, 9, 10, 11, 12, 13, 14, 15], [4]), + ([0, 3, 4, 5, 9, 12, 13, 14], [3]), + ([0, 2, 3, 7, 8, 9, 11, 12], [2]), + ([1, 2, 3, 4, 10, 11, 12, 13], [1]), + ([0, 1, 2, 3, 4, 5, 6, 7], [0])] + + TESTS:: + + sage: pr = Polyhedron(rays = [[1,0,0],[-1,0,0],[0,1,0]], vertices = [[-1,-1,-1]], lines=[(0,0,1)]) + sage: pr.faces(4) + () + sage: pr.faces(3)[0].ambient_V_indices() + (0, 1, 2, 3) + sage: pr.facets()[0].ambient_V_indices() + (0, 1, 2) + sage: pr.faces(1) + () + sage: pr.faces(0) + () + sage: pr.faces(-1) + (A -1-dimensional face of a Polyhedron in QQ^3,) + """ + return tuple(self.face_generator(face_dimension)) + + def facets(self): + r""" + Return the facets of the polyhedron. + + Facets are the maximal nontrivial faces of polyhedra. + The empty face and the polyhedron itself are trivial. + + A facet of a `d`-dimensional polyhedron is a face of dimension + `d-1`. For `d \neq 0` the converse is true as well. + + OUTPUT: + + A tuple of + :class:`~sage.geometry.polyhedron.face.PolyhedronFace`. See + :mod:`~sage.geometry.polyhedron.face` for details. The order + is random but fixed. + + .. SEEALSO:: :meth:`facets` + + EXAMPLES: + + Here we find the eight three-dimensional facets of the + four-dimensional hypercube:: + + sage: p = polytopes.hypercube(4) + sage: p.facets() + (A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices) + + This is the same result as explicitly finding the + three-dimensional faces:: + + sage: dim = p.dimension() + sage: p.faces(dim-1) + (A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices, + A 3-dimensional face of a Polyhedron in ZZ^4 defined as the convex hull of 8 vertices) + + The ``0``-dimensional polyhedron does not have facets:: + + sage: P = Polyhedron([[0]]) + sage: P.facets() + () + """ + if self.dimension() == 0: + return () + return self.faces(self.dimension()-1) + + @cached_method(do_pickle=True) + def f_vector(self, num_threads=None, parallelization_depth=None): + r""" + Return the f-vector. + + INPUT: + + - ``num_threads`` -- integer (optional); specify the number of threads; + otherwise determined by :func:`~sage.parallel.ncpus.ncpus` + + - ``parallelization_depth`` -- integer (optional); specify + how deep in the lattice the parallelization is done + + OUTPUT: + + Return a vector whose `i`-th entry is the number of + `i-2`-dimensional faces of the polytope. + + .. NOTE:: + + The ``vertices`` as given by :meth:`Polyhedron_base.vertices` + do not need to correspond to `0`-dimensional faces. If a polyhedron + contains `k` lines they correspond to `k`-dimensional faces. + See example below + + EXAMPLES:: + + sage: p = Polyhedron(vertices=[[1, 2, 3], [1, 3, 2], + ....: [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1], [0, 0, 0]]) + sage: p.f_vector() + (1, 7, 12, 7, 1) + + sage: polytopes.cyclic_polytope(4,10).f_vector() + (1, 10, 45, 70, 35, 1) + + sage: polytopes.hypercube(5).f_vector() + (1, 32, 80, 80, 40, 10, 1) + + Polyhedra with lines do not have `0`-faces:: + + sage: Polyhedron(ieqs=[[1,-1,0,0],[1,1,0,0]]).f_vector() + (1, 0, 0, 2, 1) + + However, the method :meth:`Polyhedron_base.vertices` returns + two points that belong to the ``Vrepresentation``:: + + sage: P = Polyhedron(ieqs=[[1,-1,0],[1,1,0]]) + sage: P.vertices() + (A vertex at (1, 0), A vertex at (-1, 0)) + sage: P.f_vector() + (1, 0, 2, 1) + + TESTS: + + Check that :trac:`28828` is fixed:: + + sage: P.f_vector().is_immutable() + True + + The cache of the f-vector is being pickled:: + + sage: P = polytopes.cube() + sage: P.f_vector() + (1, 8, 12, 6, 1) + sage: Q = loads(dumps(P)) + sage: Q.f_vector.is_in_cache() + True + """ + return self.combinatorial_polyhedron().f_vector(num_threads, parallelization_depth) + + def bounded_edges(self): + """ + Return the bounded edges (excluding rays and lines). + + OUTPUT: + + A generator for pairs of vertices, one pair per edge. + + EXAMPLES:: + + sage: p = Polyhedron(vertices=[[1,0],[0,1]], rays=[[1,0],[0,1]]) + sage: [ e for e in p.bounded_edges() ] + [(A vertex at (0, 1), A vertex at (1, 0))] + sage: for e in p.bounded_edges(): print(e) + (A vertex at (0, 1), A vertex at (1, 0)) + """ + obj = self.Vrepresentation() + for i in range(len(obj)): + if not obj[i].is_vertex(): + continue + for j in range(i+1, len(obj)): + if not obj[j].is_vertex(): + continue + if self.vertex_adjacency_matrix()[i, j] == 0: + continue + yield (obj[i], obj[j]) + + @cached_method + def vertex_adjacency_matrix(self): + """ + Return the binary matrix of vertex adjacencies. + + EXAMPLES:: + + sage: polytopes.simplex(4).vertex_adjacency_matrix() + [0 1 1 1 1] + [1 0 1 1 1] + [1 1 0 1 1] + [1 1 1 0 1] + [1 1 1 1 0] + + The rows and columns of the vertex adjacency matrix correspond + to the :meth:`Vrepresentation` objects: vertices, rays, and + lines. The `(i,j)` matrix entry equals `1` if the `i`-th and + `j`-th V-representation object are adjacent. + + Two vertices are adjacent if they are the endpoints of an + edge, that is, a one-dimensional face. For unbounded polyhedra + this clearly needs to be generalized and we define two + V-representation objects (see + :mod:`sage.geometry.polyhedron.constructor`) to be adjacent if + they together generate a one-face. There are three possible + combinations: + + * Two vertices can bound a finite-length edge. + + * A vertex and a ray can generate a half-infinite edge + starting at the vertex and with the direction given by the + ray. + + * A vertex and a line can generate an infinite edge. The + position of the vertex on the line is arbitrary in this + case, only its transverse position matters. The direction of + the edge is given by the line generator. + + For example, take the half-plane:: + + sage: half_plane = Polyhedron(ieqs=[(0,1,0)]) + sage: half_plane.Hrepresentation() + (An inequality (1, 0) x + 0 >= 0,) + + Its (non-unique) V-representation consists of a vertex, a ray, + and a line. The only edge is spanned by the vertex and the + line generator, so they are adjacent:: + + sage: half_plane.Vrepresentation() + (A line in the direction (0, 1), A ray in the direction (1, 0), A vertex at (0, 0)) + sage: half_plane.vertex_adjacency_matrix() + [0 0 1] + [0 0 0] + [1 0 0] + + In one dimension higher, that is for a half-space in 3 + dimensions, there is no one-dimensional face. Hence nothing is + adjacent:: + + sage: Polyhedron(ieqs=[(0,1,0,0)]).vertex_adjacency_matrix() + [0 0 0 0] + [0 0 0 0] + [0 0 0 0] + [0 0 0 0] + + EXAMPLES: + + In a bounded polygon, every vertex has precisely two adjacent ones:: + + sage: P = Polyhedron(vertices=[(0, 1), (1, 0), (3, 0), (4, 1)]) + sage: for v in P.Vrep_generator(): + ....: print("{} {}".format(P.adjacency_matrix().row(v.index()), v)) + (0, 1, 0, 1) A vertex at (0, 1) + (1, 0, 1, 0) A vertex at (1, 0) + (0, 1, 0, 1) A vertex at (3, 0) + (1, 0, 1, 0) A vertex at (4, 1) + + If the V-representation of the polygon contains vertices and + one ray, then each V-representation object is adjacent to two + V-representation objects:: + + sage: P = Polyhedron(vertices=[(0, 1), (1, 0), (3, 0), (4, 1)], + ....: rays=[(0,1)]) + sage: for v in P.Vrep_generator(): + ....: print("{} {}".format(P.adjacency_matrix().row(v.index()), v)) + (0, 1, 0, 0, 1) A ray in the direction (0, 1) + (1, 0, 1, 0, 0) A vertex at (0, 1) + (0, 1, 0, 1, 0) A vertex at (1, 0) + (0, 0, 1, 0, 1) A vertex at (3, 0) + (1, 0, 0, 1, 0) A vertex at (4, 1) + + If the V-representation of the polygon contains vertices and + two distinct rays, then each vertex is adjacent to two + V-representation objects (which can now be vertices or + rays). The two rays are not adjacent to each other:: + + sage: P = Polyhedron(vertices=[(0, 1), (1, 0), (3, 0), (4, 1)], + ....: rays=[(0,1), (1,1)]) + sage: for v in P.Vrep_generator(): + ....: print("{} {}".format(P.adjacency_matrix().row(v.index()), v)) + (0, 1, 0, 0, 0) A ray in the direction (0, 1) + (1, 0, 1, 0, 0) A vertex at (0, 1) + (0, 1, 0, 0, 1) A vertex at (1, 0) + (0, 0, 0, 0, 1) A ray in the direction (1, 1) + (0, 0, 1, 1, 0) A vertex at (3, 0) + + The vertex adjacency matrix has base ring integers. This way one can express various + counting questions:: + + sage: P = polytopes.cube() + sage: Q = P.stack(P.faces(2)[0]) + sage: M = Q.vertex_adjacency_matrix() + sage: sum(M) + (4, 4, 3, 3, 4, 4, 4, 3, 3) + sage: G = Q.vertex_graph() + sage: G.degree() + [4, 4, 3, 3, 4, 4, 4, 3, 3] + + TESTS: + + Check that :trac:`28828` is fixed:: + + sage: P.adjacency_matrix().is_immutable() + True + """ + return self.combinatorial_polyhedron().vertex_adjacency_matrix() + + adjacency_matrix = vertex_adjacency_matrix + + @cached_method + def facet_adjacency_matrix(self): + """ + Return the adjacency matrix for the facets and hyperplanes. + + EXAMPLES:: + + sage: s4 = polytopes.simplex(4, project=True) + sage: s4.facet_adjacency_matrix() + [0 1 1 1 1] + [1 0 1 1 1] + [1 1 0 1 1] + [1 1 1 0 1] + [1 1 1 1 0] + + The facet adjacency matrix has base ring integers. This way one can express various + counting questions:: + + sage: P = polytopes.cube() + sage: Q = P.stack(P.faces(2)[0]) + sage: M = Q.facet_adjacency_matrix() + sage: sum(M) + (4, 4, 4, 4, 3, 3, 3, 3, 4) + + TESTS: + + Check that :trac:`28828` is fixed:: + + sage: s4.facet_adjacency_matrix().is_immutable() + True + """ + return self._facet_adjacency_matrix() + + def _facet_adjacency_matrix(self): + """ + Compute the facet adjacency matrix in case it has not been + computed during initialization. + + EXAMPLES:: + + sage: p = Polyhedron(vertices=[(0,0),(1,0),(0,1)]) + sage: p._facet_adjacency_matrix() + [0 1 1] + [1 0 1] + [1 1 0] + + Checks that :trac:`22455` is fixed:: + + sage: s = polytopes.simplex(2) + sage: s._facet_adjacency_matrix() + [0 1 1] + [1 0 1] + [1 1 0] + + """ + # TODO: This implementation computes the whole face lattice, + # which is much more information than necessary. + M = matrix(ZZ, self.n_facets(), self.n_facets(), 0) + codim = self.ambient_dim()-self.dim() + + def set_adjacent(h1, h2): + if h1 is h2: + return + i = h1.index() - codim + j = h2.index() - codim + M[i, j] = 1 + M[j, i] = 1 + + for face in self.faces(self.dim()-2): + Hrep = face.ambient_Hrepresentation() + assert(len(Hrep) == codim+2) + set_adjacent(Hrep[-2], Hrep[-1]) + M.set_immutable() + return M + + def a_maximal_chain(self): + r""" + Return a maximal chain of the face lattice in increasing order. + + EXAMPLES:: + + sage: P = polytopes.cube() + sage: P.a_maximal_chain() + [A -1-dimensional face of a Polyhedron in ZZ^3, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices] + sage: P = polytopes.cube() + sage: chain = P.a_maximal_chain(); chain + [A -1-dimensional face of a Polyhedron in ZZ^3, + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex, + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices, + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices, + A 3-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 8 vertices] + sage: [face.ambient_V_indices() for face in chain] + [(), (5,), (0, 5), (0, 3, 4, 5), (0, 1, 2, 3, 4, 5, 6, 7)] + + TESTS:: + + Check output for the empty polyhedron:: + + sage: P = Polyhedron() + sage: P.a_maximal_chain() + [A -1-dimensional face of a Polyhedron in ZZ^0] + """ + comb_chain = self.combinatorial_polyhedron().a_maximal_chain() + + from sage.geometry.polyhedron.face import combinatorial_face_to_polyhedral_face + empty_face = self.faces(-1)[0] + universe = self.faces(self.dim())[0] + + if self.dim() == -1: + return [empty_face] + + return [empty_face] + \ + [combinatorial_face_to_polyhedral_face(self, face) + for face in comb_chain] + \ + [universe] + + def is_simplex(self): + r""" + Return whether the polyhedron is a simplex. + + A simplex is a bounded polyhedron with `d+1` vertices, where + `d` is the dimension. + + EXAMPLES:: + + sage: Polyhedron([(0,0,0), (1,0,0), (0,1,0)]).is_simplex() + True + sage: polytopes.simplex(3).is_simplex() + True + sage: polytopes.hypercube(3).is_simplex() + False + """ + return self.is_compact() and (self.dim()+1 == self.n_vertices()) + + def simplicity(self): + r""" + Return the largest integer `k` such that the polytope is `k`-simple. + + A polytope `P` is `k`-simple, if every `(d-1-k)`-face + is contained in exactly `k+1` facets of `P` for `1 \leq k \leq d-1`. + Equivalently it is `k`-simple if the polar/dual polytope is `k`-simplicial. + If `self` is a simplex, it returns its dimension. + + EXAMPLES:: + + sage: polytopes.hypersimplex(4,2).simplicity() + 1 + sage: polytopes.hypersimplex(5,2).simplicity() + 2 + sage: polytopes.hypersimplex(6,2).simplicity() + 3 + sage: polytopes.simplex(3).simplicity() + 3 + sage: polytopes.simplex(1).simplicity() + 1 + + The method is not implemented for unbounded polyhedra:: + + sage: p = Polyhedron(vertices=[(0,0)],rays=[(1,0),(0,1)]) + sage: p.simplicity() + Traceback (most recent call last): + ... + NotImplementedError: this function is implemented for polytopes only + """ + if not(self.is_compact()): + raise NotImplementedError("this function is implemented for polytopes only") + return self.combinatorial_polyhedron().simplicity() + + def is_simple(self): + """ + Test for simplicity of a polytope. + + See :wikipedia:`Simple_polytope` + + EXAMPLES:: + + sage: p = Polyhedron([[0,0,0],[1,0,0],[0,1,0],[0,0,1]]) + sage: p.is_simple() + True + sage: p = Polyhedron([[0,0,0],[4,4,0],[4,0,0],[0,4,0],[2,2,2]]) + sage: p.is_simple() + False + """ + if not self.is_compact(): + return False + return self.combinatorial_polyhedron().is_simple() + + def simpliciality(self): + r""" + Return the largest integer `k` such that the polytope is `k`-simplicial. + + A polytope is `k`-simplicial, if every `k`-face is a simplex. + If `self` is a simplex, returns its dimension. + + EXAMPLES:: + + sage: polytopes.cyclic_polytope(10,4).simpliciality() + 3 + sage: polytopes.hypersimplex(5,2).simpliciality() + 2 + sage: polytopes.cross_polytope(4).simpliciality() + 3 + sage: polytopes.simplex(3).simpliciality() + 3 + sage: polytopes.simplex(1).simpliciality() + 1 + + The method is not implemented for unbounded polyhedra:: + + sage: p = Polyhedron(vertices=[(0,0)],rays=[(1,0),(0,1)]) + sage: p.simpliciality() + Traceback (most recent call last): + ... + NotImplementedError: this function is implemented for polytopes only + """ + if not(self.is_compact()): + raise NotImplementedError("this function is implemented for polytopes only") + return self.combinatorial_polyhedron().simpliciality() + + def is_simplicial(self): + """ + Tests if the polytope is simplicial + + A polytope is simplicial if every facet is a simplex. + + See :wikipedia:`Simplicial_polytope` + + EXAMPLES:: + + sage: p = polytopes.hypercube(3) + sage: p.is_simplicial() + False + sage: q = polytopes.simplex(5, project=True) + sage: q.is_simplicial() + True + sage: p = Polyhedron([[0,0,0],[1,0,0],[0,1,0],[0,0,1]]) + sage: p.is_simplicial() + True + sage: q = Polyhedron([[1,1,1],[-1,1,1],[1,-1,1],[-1,-1,1],[1,1,-1]]) + sage: q.is_simplicial() + False + sage: P = polytopes.simplex(); P + A 3-dimensional polyhedron in ZZ^4 defined as the convex hull of 4 vertices + sage: P.is_simplicial() + True + + The method is not implemented for unbounded polyhedra:: + + sage: p = Polyhedron(vertices=[(0,0)],rays=[(1,0),(0,1)]) + sage: p.is_simplicial() + Traceback (most recent call last): + ... + NotImplementedError: this function is implemented for polytopes only + """ + if not(self.is_compact()): + raise NotImplementedError("this function is implemented for polytopes only") + return self.combinatorial_polyhedron().is_simplicial() + + 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 (1, 0, 0, 0)) + sage: egyptian_pyramid = polytopes.regular_polygon(4).pyramid() # optional - sage.rings.number_field + sage: egyptian_pyramid.is_pyramid() # optional - sage.rings.number_field + True + 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") + + return self.combinatorial_polyhedron().is_pyramid(certificate) + + 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, 6, 36, 216), + A vertex at (1, 0, 0, 0), + A vertex at (1, 7, 49, 343), + A vertex at (1, 5, 25, 125), + A vertex at (1, 1, 1, 1), + A vertex at (1, 2, 4, 8), + A vertex at (1, 4, 16, 64), + A vertex at (1, 3, 9, 27)], + [A vertex at (0, 3, 9, 27), + A vertex at (0, 6, 36, 216), + A vertex at (0, 0, 0, 0), + A vertex at (0, 7, 49, 343), + A vertex at (0, 5, 25, 125), + A vertex at (0, 1, 1, 1), + A vertex at (0, 2, 4, 8), + A vertex at (0, 4, 16, 64)]]) + + 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 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. + + EXAMPLES:: + + sage: P = polytopes.hypersimplex(5,2) + sage: L = P.lawrence_polytope() + sage: L.is_lattice_polytope() + True + sage: egyptian_pyramid = polytopes.regular_polygon(4).pyramid() + sage: egyptian_pyramid.is_lawrence_polytope() + True + sage: polytopes.octahedron().is_lawrence_polytope() + False + + REFERENCES: + + For more information, see [BaSt1990]_. + """ + if not self.is_compact(): + raise NotImplementedError("self must be a polytope") + + return self.combinatorial_polyhedron().is_lawrence_polytope() + + def is_neighborly(self, k=None): + r""" + Return whether the polyhedron is neighborly. + + If the input ``k`` is provided, then return whether the polyhedron is ``k``-neighborly + + A polyhedron is neighborly if every set of `n` vertices forms a face + for `n` up to floor of half the dimension of the polyhedron. + It is `k`-neighborly if this is true for `n` up to `k`. + + INPUT: + + - ``k`` -- the dimension up to which to check if every set of ``k`` + vertices forms a face. If no ``k`` is provided, check up to floor + of half the dimension of the polyhedron. + + OUTPUT: + + - ``True`` if every set of up to ``k`` vertices forms a face, + - ``False`` otherwise + + .. SEEALSO:: + + :meth:`neighborliness` + + EXAMPLES:: + + sage: cube = polytopes.hypercube(3) + sage: cube.is_neighborly() + True + sage: cube = polytopes.hypercube(4) + sage: cube.is_neighborly() + False + + Cyclic polytopes are neighborly:: + + sage: all(polytopes.cyclic_polytope(i, i + 1 + j).is_neighborly() for i in range(5) for j in range(3)) + True + + The neighborliness of a polyhedron equals floor of dimension half + (or larger in case of a simplex) if and only if the polyhedron + is neighborly:: + + sage: testpolys = [polytopes.cube(), polytopes.cyclic_polytope(6, 9), polytopes.simplex(6)] + sage: [(P.neighborliness()>=floor(P.dim()/2)) == P.is_neighborly() for P in testpolys] + [True, True, True] + + """ + return self.combinatorial_polyhedron().is_neighborly() + + def join_of_Vrep(self, *Vrepresentatives): + r""" + Return the smallest face that contains ``Vrepresentatives``. + + INPUT: + + - ``Vrepresentatives`` -- vertices/rays/lines of ``self`` or indices of such + + OUTPUT: a :class:`~sage.geometry.polyhedron.face.PolyhedronFace` + + .. NOTE:: + + In the case of unbounded polyhedra, the join of rays etc. may not be well-defined. + + EXAMPLES:: + + sage: P = polytopes.permutahedron(5) + sage: P.join_of_Vrep(1) + A 0-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 1 vertex + sage: P.join_of_Vrep() + A -1-dimensional face of a Polyhedron in ZZ^5 + sage: P.join_of_Vrep(0,12,13).ambient_V_indices() + (0, 12, 13, 68) + + The input is flexible:: + + sage: P.join_of_Vrep(2, P.vertices()[3], P.Vrepresentation(4)) + A 2-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 6 vertices + + :: + + sage: P = polytopes.cube() + sage: a, b = P.faces(0)[:2] + sage: P.join_of_Vrep(a, b) + A 1-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 2 vertices + + In the case of an unbounded polyhedron, the join may not be well-defined:: + + sage: P = Polyhedron(vertices=[[1,0], [0,1]], rays=[[1,1]]) + sage: P.join_of_Vrep(0) + A 0-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex + sage: P.join_of_Vrep(0,1) + A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 2 vertices + sage: P.join_of_Vrep(0,2) + A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray + sage: P.join_of_Vrep(1,2) + A 1-dimensional face of a Polyhedron in QQ^2 defined as the convex hull of 1 vertex and 1 ray + sage: P.join_of_Vrep(2) + Traceback (most recent call last): + ... + ValueError: the join is not well-defined + + The ``Vrepresentatives`` must be of ``self``:: + + sage: P = polytopes.cube(backend='ppl') + sage: Q = polytopes.cube(backend='field') + sage: v = P.vertices()[0] + sage: P.join_of_Vrep(v) + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex + sage: Q.join_of_Vrep(v) + Traceback (most recent call last): + ... + ValueError: not a Vrepresentative of ``self`` + sage: f = P.faces(0)[0] + sage: P.join_of_Vrep(v) + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex + sage: Q.join_of_Vrep(v) + Traceback (most recent call last): + ... + ValueError: not a Vrepresentative of ``self`` + + TESTS: + + ``least_common_superface_of_Vrep`` is an alias:: + + sage: P.least_common_superface_of_Vrep(v) + A 0-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 1 vertex + sage: P.least_common_superface_of_Vrep == P.join_of_Vrep + True + + Error message for invalid input:: + + sage: P.join_of_Vrep('foo') + Traceback (most recent call last): + ... + ValueError: foo is not a Vrepresentative + """ + from sage.geometry.polyhedron.representation import Vrepresentation + from sage.geometry.polyhedron.face import PolyhedronFace + + new_indices = [0]*len(Vrepresentatives) + for i, v in enumerate(Vrepresentatives): + if isinstance(v, PolyhedronFace) and v.dim() == 0: + if v.polyhedron() is not self: + raise ValueError("not a Vrepresentative of ``self``") + new_indices[i] = v.ambient_V_indices()[0] + elif v in ZZ: + new_indices[i] = v + elif isinstance(v, Vrepresentation): + if v.polyhedron() is not self: + raise ValueError("not a Vrepresentative of ``self``") + new_indices[i] = v.index() + else: + raise ValueError("{} is not a Vrepresentative".format(v)) + + return self.face_generator().join_of_Vrep(*new_indices) + + least_common_superface_of_Vrep = join_of_Vrep + + def meet_of_Hrep(self, *Hrepresentatives): + r""" + Return the largest face that is contained in ``Hrepresentatives``. + + INPUT: + + - ``Hrepresentatives`` -- facets or indices of Hrepresentatives; + the indices are assumed to be the indices of the Hrepresentation + + OUTPUT: a :class:`~sage.geometry.polyhedron.face.PolyhedronFace` + + EXAMPLES:: + + sage: P = polytopes.permutahedron(5) + sage: P.meet_of_Hrep() + A 4-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 120 vertices + sage: P.meet_of_Hrep(1) + A 3-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 24 vertices + sage: P.meet_of_Hrep(4) + A 3-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 12 vertices + sage: P.meet_of_Hrep(1,3,7) + A 1-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 2 vertices + sage: P.meet_of_Hrep(1,3,7).ambient_H_indices() + (0, 1, 3, 7) + + The indices are the indices of the Hrepresentation. + ``0`` corresponds to an equation and is ignored:: + + sage: P.meet_of_Hrep(0) + A 4-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 120 vertices + + The input is flexible:: + + sage: P.meet_of_Hrep(P.facets()[-1], P.inequalities()[2], 7) + A 1-dimensional face of a Polyhedron in ZZ^5 defined as the convex hull of 2 vertices + + The ``Hrepresentatives`` must belong to ``self``:: + + sage: P = polytopes.cube(backend='ppl') + sage: Q = polytopes.cube(backend='field') + sage: f = P.facets()[0] + sage: P.meet_of_Hrep(f) + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices + sage: Q.meet_of_Hrep(f) + Traceback (most recent call last): + ... + ValueError: not a facet of ``self`` + sage: f = P.inequalities()[0] + sage: P.meet_of_Hrep(f) + A 2-dimensional face of a Polyhedron in ZZ^3 defined as the convex hull of 4 vertices + sage: Q.meet_of_Hrep(f) + Traceback (most recent call last): + ... + ValueError: not a facet of ``self`` + + TESTS: + + Equations are not considered by the combinatorial polyhedron. + We check that the index corresponds to the Hrepresentation index:: + + sage: P = polytopes.permutahedron(3, backend='field') + sage: P.Hrepresentation() + (An inequality (0, 0, 1) x - 1 >= 0, + An inequality (0, 1, 0) x - 1 >= 0, + An inequality (0, 1, 1) x - 3 >= 0, + An inequality (1, 0, 0) x - 1 >= 0, + An inequality (1, 0, 1) x - 3 >= 0, + An inequality (1, 1, 0) x - 3 >= 0, + An equation (1, 1, 1) x - 6 == 0) + sage: P.meet_of_Hrep(0).ambient_Hrepresentation() + (An inequality (0, 0, 1) x - 1 >= 0, An equation (1, 1, 1) x - 6 == 0) + + sage: P = polytopes.permutahedron(3, backend='ppl') + sage: P.Hrepresentation() + (An equation (1, 1, 1) x - 6 == 0, + An inequality (1, 1, 0) x - 3 >= 0, + An inequality (-1, -1, 0) x + 5 >= 0, + An inequality (0, 1, 0) x - 1 >= 0, + An inequality (-1, 0, 0) x + 3 >= 0, + An inequality (1, 0, 0) x - 1 >= 0, + An inequality (0, -1, 0) x + 3 >= 0) + sage: P.meet_of_Hrep(1).ambient_Hrepresentation() + (An equation (1, 1, 1) x - 6 == 0, An inequality (1, 1, 0) x - 3 >= 0) + + ``greatest_common_subface_of_Hrep`` is an alias:: + + sage: P.greatest_common_subface_of_Hrep(1).ambient_Hrepresentation() + (An equation (1, 1, 1) x - 6 == 0, An inequality (1, 1, 0) x - 3 >= 0) + sage: P.greatest_common_subface_of_Hrep == P.meet_of_Hrep + True + + Error message for invalid input:: + + sage: P.meet_of_Hrep('foo') + Traceback (most recent call last): + ... + ValueError: foo is not a Hrepresentative + """ + from sage.geometry.polyhedron.representation import Inequality, Equation + from sage.geometry.polyhedron.face import PolyhedronFace + + # Equations are ignored by combinatorial polyhedron for indexing. + offset = 0 + if self.n_equations() and self.Hrepresentation(0).is_equation(): + offset = self.n_equations() + + new_indices = [] + for i, facet in enumerate(Hrepresentatives): + if isinstance(facet, PolyhedronFace) and facet.dim() + 1 == self.dim(): + if facet.polyhedron() is not self: + raise ValueError("not a facet of ``self``") + H_indices = facet.ambient_H_indices() + facet = H_indices[0] if H_indices[0] >= offset else H_indices[-1] + + if facet in ZZ and facet >= offset: + # Note that ``CombinatorialPolyhedron`` ignores indices of equations + # and has equations last. + new_indices.append(facet - offset) + elif isinstance(facet, Inequality): + if facet.polyhedron() is not self: + raise ValueError("not a facet of ``self``") + new_indices.append(facet.index() - offset) + elif isinstance(facet, Equation): + # Ignore equations. + continue + elif facet in ZZ and 0 <= facet < offset: + # Ignore equations. + continue + else: + raise ValueError("{} is not a Hrepresentative".format(facet)) + + return self.face_generator().meet_of_Hrep(*new_indices) + + greatest_common_subface_of_Hrep = meet_of_Hrep + + def _test_combinatorial_face_as_combinatorial_polyhedron(self, tester=None, **options): + """ + Run tests on obtaining the combinatorial face as combinatorial polyhedron. + + TESTS:: + + sage: polytopes.cross_polytope(3)._test_combinatorial_face_as_combinatorial_polyhedron() + """ + if not self.is_compact(): + return + if self.dim() > 7 or self.n_vertices() > 100 or self.n_facets() > 100: + # Avoid very long tests. + return + if self.dim() < 1: + # Avoid trivial cases. + return + + from sage.misc.prandom import random + + if tester is None: + tester = self._tester(**options) + + it = self.face_generator() + _ = next(it), next(it) # get rid of non_proper faces + C1 = self.combinatorial_polyhedron() + it1 = C1.face_iter() + C2 = C1.dual() + it2 = C2.face_iter(dual=not it1.dual) + + for f in it: + f1 = next(it1) + f2 = next(it2) + if random() < 0.95: + # Only test a random 5 percent of the faces. + continue + + P = f.as_polyhedron() + D1 = f1.as_combinatorial_polyhedron() + D2 = f2.as_combinatorial_polyhedron(quotient=True).dual() + D1._test_bitsets(tester, **options) + D2._test_bitsets(tester, **options) + try: + import sage.graphs.graph + except ImportError: + pass + else: + tester.assertTrue(P.combinatorial_polyhedron().vertex_facet_graph().is_isomorphic(D1.vertex_facet_graph())) + tester.assertTrue(P.combinatorial_polyhedron().vertex_facet_graph().is_isomorphic(D2.vertex_facet_graph()))