diff --git a/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pxd b/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pxd index c44f41ff6ad..1cc816276f2 100644 --- a/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pxd +++ b/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pxd @@ -13,13 +13,14 @@ cdef class CombinatorialPolyhedron(SageObject): cdef tuple _H # the names of HRep, if they exist cdef tuple _equalities # stores equalities, given on input (might belong to Hrep) cdef int _dimension # stores dimension, -2 on init - cdef unsigned int _length_Hrep # Hrep might include equalities - cdef unsigned int _length_Vrepr # Vrep might include rays/lines + cdef unsigned int _length_Hrepr # Hrepr might include equalities + cdef unsigned int _length_Vrepr # Vrepr might include rays/lines cdef size_t _n_facets # length Hrep without equalities cdef bint _unbounded # ``True`` iff Polyhedron is unbounded - cdef int _n_lines # number of affinely independent lines in the Polyhedron cdef ListOfFaces bitrep_facets # facets in bit representation cdef ListOfFaces bitrep_Vrepr # vertices in bit representation + cdef ListOfFaces far_face # a 'face' containing all none-vertices of Vrepr + cdef tuple far_face_tuple cdef tuple _f_vector # Edges, ridges and incidences are stored in a pointer of pointers. diff --git a/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx b/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx index ef7aa76e997..078767c26b8 100644 --- a/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx +++ b/src/sage/geometry/polyhedron/combinatorial_polyhedron/base.pyx @@ -117,7 +117,6 @@ cdef class CombinatorialPolyhedron(SageObject): * or an ``incidence_matrix`` as in :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.incidence_matrix` In this case you should also specify the ``Vrepr`` and ``facets`` arguments - If the polyhedron is unbounded, then ``n_lines`` as well * or list of facets, each facet given as a list of ``[vertices, rays, lines]`` if the polyhedron is unbounded, then rays and lines and the extra argument ``nr_lines`` are required @@ -133,9 +132,12 @@ cdef class CombinatorialPolyhedron(SageObject): - ``facets`` -- (optional) when ``data`` is an incidence matrix or a list of facets, it should be a list of facets that would be used instead of indices (of the columns of the incidence matrix). - - ``n_lines`` -- (semi-optional) when ``data` is an incidence matrix or a + - ``unbounded`` -- value will be overwritten if ``data`` is a polyhedron; + if ``unbounded`` and ``data`` is incidence matrix or a list of facets, + need to specify ``far_face`` + - ``far_face`` -- (semi-optional) when ``data` is an incidence matrix or a list of facets and the polyhedron is unbounded this needs to be set to - the dimension of the maximal linear subspace contained in the polyhedron + the list of indices of the rays and lines EXAMPLES: @@ -155,11 +157,13 @@ cdef class CombinatorialPolyhedron(SageObject): an incidence matrix:: - sage: data = Polyhedron(rays=[[0,1]]).incidence_matrix() - sage: CombinatorialPolyhedron(data, n_lines=0) + sage: P = Polyhedron(rays=[[0,1]]) + sage: data = P.incidence_matrix() + sage: far_face = [i for i in range(2) if not P.Vrepresentation()[i].is_vertex()] + sage: CombinatorialPolyhedron(data, unbounded=True, far_face=far_face) A 1-dimensional combinatorial polyhedron with 1 facet sage: C = CombinatorialPolyhedron(data, Vrepr=['myvertex'], - ....: facets=['myfacet'], n_lines=0) + ....: facets=['myfacet'], unbounded=True, far_face=far_face) sage: C.Vrepresentation() ('myvertex',) sage: C.Hrepresentation() @@ -186,7 +190,7 @@ cdef class CombinatorialPolyhedron(SageObject): sage: CombinatorialPolyhedron(5).f_vector() (1, 0, 0, 0, 0, 0, 1) - Specifying the number of lines is important. The following with a + Specifying that a polyhedron is unbounded is important. The following with a polyhedron works fine:: sage: P = Polyhedron(ieqs=[[1,-1,0],[1,1,0]]) @@ -194,7 +198,7 @@ cdef class CombinatorialPolyhedron(SageObject): sage: C A 2-dimensional combinatorial polyhedron with 2 facets - But it becomes incorrect here due to the missing number of lines:: + The following is incorrect, as ``unbounded`` is implicitly set to ``False``:: sage: data = P.incidence_matrix() sage: vert = P.Vrepresentation() @@ -204,13 +208,12 @@ cdef class CombinatorialPolyhedron(SageObject): sage: C.f_vector() (1, 1, 2, 1) sage: C.vertices() - (A line in the direction (0, 1), - A line in the direction (0, 1), - A line in the direction (0, 1)) + (A line in the direction (0, 1), A vertex at (1, 0), A vertex at (-1, 0)) The correct usage is:: - sage: C = CombinatorialPolyhedron(data, Vrepr=vert, n_lines=1) + sage: far_face = [i for i in range(3) if not P.Vrepresentation()[i].is_vertex()] + sage: C = CombinatorialPolyhedron(data, Vrepr=vert, unbounded=True, far_face=far_face) sage: C A 2-dimensional combinatorial polyhedron with 2 facets sage: C.f_vector() @@ -218,7 +221,26 @@ cdef class CombinatorialPolyhedron(SageObject): sage: C.vertices() () - TESTS:: + TESTS: + + Checking that :trac:`27987` is fixed:: + + sage: P1 = Polyhedron(vertices=[[0,1],[1,0]], rays=[[1,1]]) + sage: P2 = Polyhedron(vertices=[[0,1],[1,0],[1,1]]) + sage: P1.incidence_matrix() == P2.incidence_matrix() + True + sage: CombinatorialPolyhedron(P1).f_vector() + (1, 2, 3, 1) + sage: CombinatorialPolyhedron(P2).f_vector() + (1, 3, 3, 1) + sage: P1 = Polyhedron(vertices=[[0,1],[1,0]], rays=[[1,1]]) + sage: P2 = Polyhedron(vertices=[[0,1],[1,0],[1,1]]) + sage: CombinatorialPolyhedron(P1).f_vector() + (1, 2, 3, 1) + sage: CombinatorialPolyhedron(P2).f_vector() + (1, 3, 3, 1) + + Some other tests regaring small polyhedra:: sage: P = Polyhedron(rays=[[1,0],[0,1]]) sage: C = CombinatorialPolyhedron(P) @@ -236,8 +258,11 @@ cdef class CombinatorialPolyhedron(SageObject): sage: C.f_vector() (1, 1, 2, 1) sage: C.vertices() - (A vertex at (0, 0), A vertex at (0, 0), A vertex at (0, 0)) - sage: C = CombinatorialPolyhedron(data, Vrepr=vert, n_lines=0) + (A vertex at (0, 0), + A ray in the direction (0, 1), + A ray in the direction (1, 0)) + sage: far_face = [i for i in range(3) if not P.Vrepresentation()[i].is_vertex()] + sage: C = CombinatorialPolyhedron(data, Vrepr=vert, unbounded=True, far_face=far_face) sage: C A 2-dimensional combinatorial polyhedron with 2 facets sage: C.f_vector() @@ -246,8 +271,10 @@ cdef class CombinatorialPolyhedron(SageObject): (A vertex at (0, 0),) sage: CombinatorialPolyhedron(3r) A 3-dimensional combinatorial polyhedron with 0 facets + + """ - def __init__(self, data, Vrepr=None, facets=None, n_lines=None): + def __init__(self, data, Vrepr=None, facets=None, unbounded=False, far_face=None): r""" Initialize :class:`CombinatorialPolyhedron`. @@ -285,36 +312,29 @@ cdef class CombinatorialPolyhedron(SageObject): if not data.is_compact(): self._unbounded = True - self._n_lines = int(data.n_lines()) + far_face = tuple(i for i in range(data.n_Vrepresentation()) if not data.Vrepresentation()[i].is_vertex()) else: self._unbounded = False - self._n_lines = 0 data = data.incidence_matrix() elif isinstance(data, LatticePolytopeClass): # input is ``LatticePolytope`` self._unbounded = False - self._n_lines = 0 Vrepr = data.vertices() self._length_Vrepr = len(Vrepr) facets = data.facets() - self._length_Hrep = len(facets) + self._length_Hrepr = len(facets) data = tuple(tuple(vert for vert in facet.vertices()) for facet in facets) else: # Input is different from ``Polyhedron`` and ``LatticePolytope``. - if n_lines is None: + if not unbounded: # bounded polyhedron self._unbounded = False - self._n_lines = 0 - elif n_lines >= 0: - # unbounded polyhedron - # will be slower but not incorrect if ``n_lines == 0`` - assert isinstance(n_lines, numbers.Integral), "n_lines need to be an integer" - self._unbounded = True - self._n_lines = int(n_lines) + elif not far_face: + raise ValueError("must specify far face for unbounded polyhedron") else: - raise ValueError("n_lines must be a non-negative integer") + self._unbounded = True if Vrepr: # store vertices names @@ -344,7 +364,7 @@ cdef class CombinatorialPolyhedron(SageObject): if isinstance(data, Matrix): # Input is incidence-matrix or was converted to it. - self._length_Hrep = data.ncols() + self._length_Hrepr = data.ncols() self._length_Vrepr = data.nrows() # Initializing the facets in their Bit-representation. @@ -355,6 +375,12 @@ cdef class CombinatorialPolyhedron(SageObject): self._n_facets = self.bitrep_facets.n_faces + # Initialize far_face if unbounded. + if self._unbounded: + self.far_face = facets_tuple_to_bit_repr_of_facets((tuple(far_face),), self._length_Vrepr) + else: + self.far_face = None + elif isinstance(data, numbers.Integral): # To construct a trivial polyhedron, equal to its affine hull, # one can give an Integer as Input. @@ -369,6 +395,8 @@ cdef class CombinatorialPolyhedron(SageObject): # Initializing the Vrepr as their Bit-representation. self.bitrep_Vrepr = facets_tuple_to_bit_repr_of_Vrepr((), 0) + self.far_face = None + else: # Input is a "list" of facets. # The facets given by its ``[vertices, rays, lines]``. @@ -398,7 +426,7 @@ cdef class CombinatorialPolyhedron(SageObject): facets = tuple(tuple(f(i) for i in j) for j in data) self._n_facets = len(facets) - self._length_Hrep = len(facets) + self._length_Hrepr = len(facets) # Initializing the facets in their Bit-representation. self.bitrep_facets = facets_tuple_to_bit_repr_of_facets(facets, length_Vrepr) @@ -406,6 +434,17 @@ cdef class CombinatorialPolyhedron(SageObject): # Initializing the Vrepr as their Bit-representation. self.bitrep_Vrepr = facets_tuple_to_bit_repr_of_Vrepr(facets, length_Vrepr) + # Initialize far_face if unbounded. + if self._unbounded: + self.far_face = facets_tuple_to_bit_repr_of_facets((tuple(far_face),), length_Vrepr) + else: + self.far_face = None + + if self._unbounded: + self.far_face_tuple = tuple(far_face) + else: + self.far_face_tuple = () + def _repr_(self): r""" Return a description of the combinatorial polyhedron. @@ -490,12 +529,14 @@ cdef class CombinatorialPolyhedron(SageObject): sage: tup == tup1 True """ - n_lines = None - if self._unbounded: - n_lines = smallInteger(self._n_lines) # Give a constructor by list of facets. - return (CombinatorialPolyhedron, (self.facets(), - self.Vrepresentation(), self.Hrepresentation(), n_lines)) + if self._unbounded: + return (CombinatorialPolyhedron, (self.facets(), + self.Vrepresentation(), self.Hrepresentation(), + True, self.far_face_tuple)) + else: + return (CombinatorialPolyhedron, (self.facets(), + self.Vrepresentation(), self.Hrepresentation())) def Vrepresentation(self): r""" @@ -537,7 +578,7 @@ cdef class CombinatorialPolyhedron(SageObject): if self._H is not None: return self._equalities + self._H else: - return tuple(smallInteger(i) for i in range(self._length_Hrep)) + return tuple(smallInteger(i) for i in range(self._length_Hrepr)) def dimension(self): r""" @@ -633,12 +674,11 @@ cdef class CombinatorialPolyhedron(SageObject): if self.dimension() == 0: # This specific trivial polyhedron needs special attention. return smallInteger(1) - elif not self._unbounded: - # In the unbounded case, we need to actually computed the vertices, - # the the Vrepr contains also ``lines`` and ``rays``. - return Integer(self._length_Vrepr) + if self._unbounded: + # Some elements in the ``Vrepr`` might not correspond to actual combinatorial vertices. + return len(self.vertices()) else: - return Integer(len(self.vertices())) + return smallInteger(self._length_Vrepr) def vertices(self, names=True): r""" @@ -664,24 +704,23 @@ cdef class CombinatorialPolyhedron(SageObject): sage: P = polytopes.cross_polytope(3) sage: C = CombinatorialPolyhedron(P) sage: C.vertices() - (A vertex at (1, 0, 0), - A vertex at (0, 1, 0), - A vertex at (0, 0, 1), - A vertex at (0, 0, -1), + (A vertex at (-1, 0, 0), A vertex at (0, -1, 0), - A vertex at (-1, 0, 0)) - + A vertex at (0, 0, -1), + A vertex at (0, 0, 1), + A vertex at (0, 1, 0), + A vertex at (1, 0, 0)) sage: C.vertices(names=False) - (5, 4, 3, 2, 1, 0) + (0, 1, 2, 3, 4, 5) sage: points = [(1,0,0), (0,1,0), (0,0,1), ....: (-1,0,0), (0,-1,0), (0,0,-1)] sage: L = LatticePolytope(points) sage: C = CombinatorialPolyhedron(L) sage: C.vertices() - (M(0, 0, -1), M(0, -1, 0), M(-1, 0, 0), M(0, 0, 1), M(0, 1, 0), M(1, 0, 0)) + (M(1, 0, 0), M(0, 1, 0), M(0, 0, 1), M(-1, 0, 0), M(0, -1, 0), M(0, 0, -1)) sage: C.vertices(names=False) - (5, 4, 3, 2, 1, 0) + (0, 1, 2, 3, 4, 5) sage: P = Polyhedron(vertices=[[0,0]]) sage: C = CombinatorialPolyhedron(P) @@ -694,16 +733,21 @@ cdef class CombinatorialPolyhedron(SageObject): return (self._V[0],) else: return (smallInteger(0),) - dual = False - if not self._unbounded: - # In the bounded case, we already know all the vertices. - # Whereas, in the unbounded case, some elements in Vrepr might - # be ``rays`` or ``lines``. - dual = True - - # Get all `0`-dimensional faces from :meth:`face_iter`. - face_iter = self.face_iter(0, dual=dual) - return tuple(face.Vrepr(names=names)[0] for face in face_iter) + if self._unbounded: + it = self.face_iter(0) + try: + # The Polyhedron has at least one vertex. + # In this case every element in the ``Vrepr`` + # that is not contained in the far face + # is a vertex. + next(it) + except StopIteration: + # The Polyhedron has no vertex. + return () + if names and self._V: + return tuple(self._V[i] for i in range(self._length_Vrepr) if not i in self.far_face_tuple) + else: + return tuple(smallInteger(i) for i in range(self._length_Vrepr) if not i in self.far_face_tuple) def n_facets(self): r""" @@ -1864,7 +1908,7 @@ cdef class CombinatorialPolyhedron(SageObject): if dim > -1: while (f_vector[dimension_one + 1] == 0): # Taking care of cases, where there might be no faces - # of dimension 0, 1, etc (``self._n_lines > 0``). + # of dimension 0, 1, etc (``n_lines > 0``). dimension_one += 1 dimension_two = -1 diff --git a/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pxd b/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pxd index 398bb2a6b45..1d3bacf095f 100644 --- a/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pxd +++ b/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pxd @@ -13,7 +13,6 @@ cdef class FaceIterator(SageObject): cdef size_t *coatom_repr # a place where coatom-representaion of face will be stored cdef int current_dimension # dimension of current face, dual dimension if ``dual`` cdef int dimension # dimension of the polyhedron - cdef int n_lines # ``_n_lines`` of ``CombinatorialPolyhedron`` cdef int output_dimension # only faces of this (dual?) dimension are considered cdef int lowest_dimension # don't consider faces below this (dual?) dimension cdef size_t _index # this counts the number of seen faces, useful for hasing the faces diff --git a/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pyx b/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pyx index 6fd8653e27c..98c94ca1395 100644 --- a/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pyx +++ b/src/sage/geometry/polyhedron/combinatorial_polyhedron/face_iterator.pyx @@ -421,14 +421,13 @@ cdef class FaceIterator(SageObject): self.face = NULL self.dimension = C.dimension() self.current_dimension = self.dimension -1 - self.n_lines = C._n_lines self._mem = MemoryAllocator() # We will not yield the empty face. # If there are `n` lines, than there # are no faces below dimension `n`. # The dimension of the level-sets in the face lattice jumps from `n` to `-1`. - self.lowest_dimension = self.n_lines + self.lowest_dimension = 0 if output_dimension is not None: if not output_dimension in range(0,self.dimension): @@ -438,7 +437,7 @@ cdef class FaceIterator(SageObject): self.output_dimension = self.dimension - 1 - output_dimension else: self.output_dimension = output_dimension - self.lowest_dimension = max(self.n_lines, self.output_dimension) + self.lowest_dimension = max(0, self.output_dimension) else: self.output_dimension = -2 @@ -478,6 +477,18 @@ cdef class FaceIterator(SageObject): self.visited_all = self._mem.allocarray(self.coatoms.n_faces, sizeof(uint64_t *)) self.n_visited_all = self._mem.allocarray(self.dimension, sizeof(size_t)) self.n_visited_all[self.dimension -1] = 0 + if C._unbounded: + # Treating the far face as if we had visited all its elements. + # Hence we will visit all intersections of facets unless contained in the far face. + + # Regarding the length of ``self.visited_all``: + # The last facet will not yield any new faces thus the length of ``visited_all`` + # needs to be at most ``n_facets - 1``. + # Hence it is fine to use the first entry already for the far face, + # as ``self.visited_all`` holds ``n_facets`` pointers. + some_list = C.far_face + self.visited_all[0] = some_list.data[0] + self.n_visited_all[self.dimension -1] = 1 # Initialize ``newfaces``, which will point to the new faces of codimension 1, # which have not been visited yet.