Trac #28982: Use CombinatorialPolyhedron to obtain faces lattice of p…

…olyhedra

We use `CombinatorialPolyhedron` to compute the face lattice of a
polyhedron.

Along the way we implement `hasse_diagram` for `CombinatorialPolyhedron`
and `Polyhedron_base`.

Instead of caching `face_lattice`, we cache `hasse_diagram` now.

This is much slower, but removes a memory leak. As `hasse_diagram` is
hardly used with #28646, this seems to be reasonable.

Caching the face lattice does create a memory leak:

{{{
sage: import gc
sage: P = polytopes.cube()
sage: P.my_face_lattice = P.face_lattice()
sage: n = get_memory_usage()
sage: P = polytopes.cube()
sage: P.my_face_lattice = P.face_lattice()
sage: _ = gc.collect()
sage: n == get_memory_usage()
False
}}}

URL: https://trac.sagemath.org/28982
Reported by: gh-kliem
Ticket author(s): Jonathan Kliem
Reviewer(s): Jean-Philippe Labbé, Matthias Koeppe

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

@@ -127,6 +127,7 @@ List of Polyhedron methods
 
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.combinatorial_polyhedron` | the combinatorial polyhedron
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.face_lattice` | the face lattice
+    :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.hasse_diagram` | the hasse diagram
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.combinatorial_automorphism_group` | the automorphism group of the underlying combinatorial polytope
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.graph`, :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.vertex_graph` | underlying graph
     :meth:`~sage.geometry.polyhedron.base.Polyhedron_base.vertex_digraph` | digraph (orientation of edges determined by a linear form)

src/sage/geometry/polyhedron/base.py

@@ -782,7 +782,7 @@ def vertex_facet_graph(self, labels=True):
              An inequality (0, 0, 1) x + 1 >= 0,
              An inequality (0, 1, 0) x + 1 >= 0,
              An inequality (1, 0, 0) x + 1 >= 0]
-            sage: G.automorphism_group().is_isomorphic(P.face_lattice().hasse_diagram().automorphism_group())
+            sage: G.automorphism_group().is_isomorphic(P.hasse_diagram().automorphism_group())
             True
             sage: O = polytopes.octahedron(); O
             A 3-dimensional polyhedron in ZZ^3 defined as the convex hull of 6 vertices
@@ -6364,7 +6364,6 @@ def barycentric_subdivision(self, subdivision_frac=None):
 
         return (polar.polar(in_affine_span=True)) + barycenter
 
-    @cached_method
     def face_lattice(self):
         """
         Return the face-lattice poset.
@@ -6398,40 +6397,20 @@ def face_lattice(self):
 
         ALGORITHM:
 
-        For a full-dimensional polytope, the basic algorithm is
-        described in
-        :func:`~sage.geometry.hasse_diagram.lattice_from_incidences`.
-        There are three generalizations of [KP2002]_ necessary to deal
-        with more general polytopes, corresponding to the extra
-        H/V-representation objects:
-
-        * Lines are removed before calling
-          :func:`lattice_from_incidences`, and then added back
-          to each face V-representation except for the "empty face".
+        See :mod:`sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator`.
 
-        * Equations are removed before calling
-          :func:`lattice_from_incidences`, and then added back
-          to each face H-representation.
+        .. NOTE::
 
-        * Rays: Consider the half line as an example. The
-          V-representation objects are a point and a ray, which we can
-          think of as a point at infinity. However, the point at
-          infinity has no inequality associated to it, so there is
-          only one H-representation object alltogether. The face
-          lattice does not contain the "face at infinity". This means
-          that in :func:`lattice_from_incidences`, one needs to
-          drop faces with V-representations that have no matching
-          H-representation. In addition, one needs to ensure that
-          every non-empty face contains at least one vertex.
+            The face lattice is not cached, as long as this creates a memory leak, see :trac:`28982`.
 
         EXAMPLES::
 
             sage: square = polytopes.hypercube(2)
             sage: fl = square.face_lattice();fl
-            Finite lattice containing 10 elements with distinguished linear extension
+            Finite lattice containing 10 elements
             sage: list(f.ambient_V_indices() for f in fl)
-            [(), (0,), (1,), (2,), (3,), (0, 1), (1, 2), (2, 3), (0, 3), (0, 1, 2, 3)]
-            sage: poset_element = fl[6]
+            [(), (0,), (1,), (0, 1), (2,), (1, 2), (3,), (0, 3), (2, 3), (0, 1, 2, 3)]
+            sage: poset_element = fl[5]
             sage: a_face = poset_element
             sage: a_face
             A 1-dimensional face of a Polyhedron in ZZ^2 defined as the convex hull of 2 vertices
@@ -6497,41 +6476,50 @@ def face_lattice(self):
             [[()], [(0, 1)]]
             sage: [[ls.ambient_V_indices() for ls in lss] for lss in Polyhedron(lines=[(1,0)], vertices=[(0,0)]).face_lattice().level_sets()]
             [[()], [(0, 1)]]
+
+        Test that computing the face lattice does not lead to a memory leak::
+
+            sage: import gc
+            sage: _ = gc.collect()
+            sage: P = polytopes.cube()
+            sage: a = P.face_lattice()
+            sage: n = get_memory_usage()
+            sage: P = polytopes.cube()
+            sage: a = P.face_lattice()
+            sage: _ = gc.collect()
+            sage: n == get_memory_usage()
+            True
         """
-        coatom_to_Hindex = [ h.index() for h in self.inequality_generator() ]
-        Hindex_to_coatom = [None] * self.n_Hrepresentation()
-        for i in range(len(coatom_to_Hindex)):
-            Hindex_to_coatom[ coatom_to_Hindex[i] ] = i
+        from sage.combinat.posets.lattices import FiniteLatticePoset
+        return FiniteLatticePoset(self.hasse_diagram())
 
-        atom_to_Vindex = [ v.index() for v in self.Vrep_generator() if not v.is_line() ]
-        Vindex_to_atom = [None] * self.n_Vrepresentation()
-        for i in range(len(atom_to_Vindex)):
-                        Vindex_to_atom[ atom_to_Vindex[i] ] = i
+    @cached_method
+    def hasse_diagram(self):
+        r"""
+        Return the Hasse diagram of the face lattice of ``self``.
 
-        atoms_incidences   = [ tuple([ Hindex_to_coatom[h.index()]
-                                       for h in v.incident() if h.is_inequality() ])
-                               for v in self.Vrepresentation() if not v.is_line() ]
+        This is the Hasse diagram of the poset of the faces of ``self``.
 
-        coatoms_incidences = [ tuple([ Vindex_to_atom[v.index()]
-                                       for v in h.incident() if not v.is_line() ])
-                               for h in self.Hrepresentation() if h.is_inequality() ]
+        OUTPUT: a directed graph
 
-        atoms_vertices = [ Vindex_to_atom[v.index()] for v in self.vertex_generator() ]
-        equations = [ e.index() for e in self.equation_generator() ]
-        lines     = [ l.index() for l in self.line_generator() ]
+        EXAMPLES::
 
-        def face_constructor(atoms, coatoms):
-            from sage.geometry.polyhedron.face import PolyhedronFace
-            if not atoms:
-                Vindices = ()
-            else:
-                Vindices = tuple(sorted([atom_to_Vindex[i] for i in atoms] + lines))
-            Hindices = tuple(sorted([coatom_to_Hindex[i] for i in coatoms] + equations))
-            return PolyhedronFace(self, Vindices, Hindices)
+            sage: P = polytopes.regular_polygon(4).pyramid()
+            sage: D = P.hasse_diagram(); D
+            Digraph on 20 vertices
+            sage: D.degree_polynomial()
+            x^5 + x^4*y + x*y^4 + y^5 + 4*x^3*y + 8*x^2*y^2 + 4*x*y^3
+        """
+
+        from sage.geometry.polyhedron.face import combinatorial_face_to_polyhedral_face
+        C = self.combinatorial_polyhedron()
+        D = C.hasse_diagram()
+
+        def index_to_polyhedron_face(n):
+            return combinatorial_face_to_polyhedral_face(
+                    self, C.face_by_face_lattice_index(n))
 
-        from sage.geometry.hasse_diagram import lattice_from_incidences
-        return lattice_from_incidences(atoms_incidences, coatoms_incidences,
-             face_constructor=face_constructor, required_atoms=atoms_vertices)
+        return D.relabel(index_to_polyhedron_face, inplace=False, immutable=True)
 
     def face_generator(self, face_dimension=None, dual=None):
         r"""
@@ -9281,10 +9269,10 @@ def combinatorial_automorphism_group(self, vertex_graph_only=False):
 
         The automorphism group of the face lattice is isomorphic to the combinatorial automorphism group::
 
-            sage: CG = C.face_lattice().hasse_diagram().automorphism_group()
+            sage: CG = C.hasse_diagram().automorphism_group()
             sage: C.combinatorial_automorphism_group().is_isomorphic(CG)
             True
-            sage: QG = Q.face_lattice().hasse_diagram().automorphism_group()
+            sage: QG = Q.hasse_diagram().automorphism_group()
             sage: Q.combinatorial_automorphism_group().is_isomorphic(QG)
             True
 

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

@@ -90,8 +90,6 @@ AUTHOR:
 import numbers
 from sage.rings.integer             import Integer
 from sage.graphs.graph              import Graph
-from sage.graphs.digraph            import DiGraph
-from sage.combinat.posets.lattices  import FiniteLatticePoset
 from sage.geometry.polyhedron.base  import Polyhedron_base
 from sage.geometry.lattice_polytope import LatticePolytopeClass
 from sage.geometry.cone             import ConvexRationalPolyhedralCone
@@ -1543,6 +1541,7 @@ cdef class CombinatorialPolyhedron(SageObject):
             sage: Polyhedron([[0]]).vertex_facet_graph(False)
             Digraph on 1 vertex
         """
+        from sage.graphs.digraph import DiGraph
         if self.dimension() == -1:
             return DiGraph()
         if self.dimension() == 0:
@@ -2396,7 +2395,7 @@ cdef class CombinatorialPolyhedron(SageObject):
 
         OUTPUT:
 
-        - :class:'~sage.combinat.posets.lattices.FiniteLatticePoset'
+        - :class:`~sage.combinat.posets.lattices.FiniteLatticePoset`
 
         .. NOTE::
 
@@ -2406,10 +2405,9 @@ cdef class CombinatorialPolyhedron(SageObject):
         .. WARNING::
 
             The labeling of the face lattice might depend on architecture
-            and implementation.
-            Relabeling the face lattice with
+            and implementation. Relabeling the face lattice with
             :meth:`CombinatorialPolyhedron.face_by_face_lattice_index` or
-             the properties obtained from this face will be platform independent
+            the properties obtained from this face will be platform independent.
 
         EXAMPLES::
 
@@ -2442,6 +2440,49 @@ cdef class CombinatorialPolyhedron(SageObject):
             sage: C.face_lattice().is_isomorphic(P.face_lattice())
             True
         """
+        from sage.combinat.posets.lattices import FiniteLatticePoset
+        return FiniteLatticePoset(self.hasse_diagram())
+
+    @cached_method
+    def hasse_diagram(self):
+        r"""
+        Return the Hasse diagram of ``self``.
+
+        This is the Hasse diagram of the poset of the faces of ``self``:
+        A directed graph consisting of a vertex for each face
+        and an edge for each minimal inclusion of faces.
+
+        .. NOTE::
+
+            The vertices of the Hasse diagram are given by indices.
+            Use :meth:`CombinatorialPolyhedron.face_by_face_lattice_index`
+            to relabel.
+
+        .. WARNING::
+
+            The indices of the Hasse diagram might depend on architecture
+            and implementation. Relabeling the face lattice with
+            :meth:`CombinatorialPolyhedron.face_by_face_lattice_index` or
+            the properties obtained from this face will be platform independent
+
+        EXAMPLES::
+
+            sage: P = polytopes.regular_polygon(4).pyramid()
+            sage: C = CombinatorialPolyhedron(P)
+            sage: D = C.hasse_diagram(); D
+            Digraph on 20 vertices
+            sage: D.average_degree()
+            21/5
+            sage: D.relabel(C.face_by_face_lattice_index)
+            sage: dim_0_vert = D.vertices()[1:6]; dim_0_vert
+            [A 0-dimensional face of a 3-dimensional combinatorial polyhedron,
+             A 0-dimensional face of a 3-dimensional combinatorial polyhedron,
+             A 0-dimensional face of a 3-dimensional combinatorial polyhedron,
+             A 0-dimensional face of a 3-dimensional combinatorial polyhedron,
+             A 0-dimensional face of a 3-dimensional combinatorial polyhedron]
+            sage: sorted(D.out_degree(vertices=dim_0_vert))
+            [3, 3, 3, 3, 4]
+        """
         if not self._face_lattice_incidences:
             # compute all incidences.
             self._compute_face_lattice_incidences()
@@ -2472,9 +2513,11 @@ cdef class CombinatorialPolyhedron(SageObject):
         edges = tuple((face_one(j), face_two(j))
                           for j in range(n_incidences))
 
-        V = tuple(range(sum(self._f_vector)))
-        D = DiGraph([V, edges], format='vertices_and_edges')
-        return FiniteLatticePoset(D)
+        V = tuple(smallInteger(i) for i in range(sum(self._f_vector)))
+
+        from sage.graphs.digraph import DiGraph
+        D = DiGraph([V, edges], format='vertices_and_edges', vertex_labels=False)
+        return D
 
     def _face_lattice_dimension(self, index):
         r"""
@@ -2572,7 +2615,7 @@ cdef class CombinatorialPolyhedron(SageObject):
             ((), (0,), (0, 1), (0, 2), (0, 1, 2))
         """
         self._record_all_faces()                            # Initalize ``_all_faces``, if not done yet.
-        dim = self._face_lattice_dimension(index)            # Determine dimension to that index.
+        dim = self._face_lattice_dimension(index)           # Determine dimension to that index.
         newindex = index - sum(self._f_vector[:dim + 1])    # Index in that level-set.
 
         # Let ``_all_faces`` determine Vrepresentation.

src/sage/geometry/polyhedron/combinatorial_polyhedron/combinatorial_face.pxd

@@ -17,7 +17,12 @@ cdef class CombinatorialFace(SageObject):
     cdef int _dimension             # dimension of current face, dual dimension if ``dual``
     cdef int _ambient_dimension     # dimension of the polyhedron
     cdef size_t face_length         # stores length of the faces in terms of uint64_t
-    cdef size_t _hash_index         # an index to give different hashes for all faces of a Polyhedron
+
+    # An index to give different hashes for all faces of a Polyhedron.
+    # The index must be chosen such that `F \subset G` implies ``hash(F) < hash(G)``.
+    cdef size_t _hash_index
+
+    cdef bint _initialized_from_face_lattice
 
     # some copies from ``CombinatorialPolyhedron``
     cdef tuple _ambient_Vrep, _ambient_facets, _equalities

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

@@ -198,6 +198,8 @@ cdef class CombinatorialFace(SageObject):
             self.coatoms            = it.coatoms
             self._hash_index        = it.structure._index
 
+            self._initialized_from_face_lattice = False
+
         elif isinstance(data, PolyhedronFaceLattice):
             all_faces = data
             assert isinstance(dimension, numbers.Integral), "dimension must be an integer"
@@ -220,12 +222,24 @@ cdef class CombinatorialFace(SageObject):
             self.atoms              = all_faces.atoms
             self.coatoms            = all_faces.coatoms
 
+            self._initialized_from_face_lattice = True
+
             self._hash_index = index
             for i in range(-1,dimension):
                 self._hash_index += all_faces.f_vector[i+1]
+
+            # Add the complete ``f-vector`` to the hash index,
+            # such that hash values obtained by an iterator or by the face lattice
+            # do not collide.
+            for i in range(-1,self._ambient_dimension+1):
+                self._hash_index += all_faces.f_vector[i+1]
         else:
             raise NotImplementedError("data must be face iterator or a list of all faces")
 
+        if self._dual:
+            # Reverse the hash index in dual mode to respect inclusion of faces.
+            self._hash_index = -self._hash_index - 1
+
     def _repr_(self):
         r"""
         Return a description of the combinatorial face.
@@ -267,12 +281,18 @@ cdef class CombinatorialFace(SageObject):
         r"""
         Return an index for the face.
 
+        This is constructed such that for faces `F,G` constructed in the same manner (same face iterator or face lattice)
+        it holds that `F` contained in `G` implies ``hash(F) < hash(G)``.
+
         If the face was constructed from a :class:`sage.geometry.polyhedron.combinatorial_polyhedron.face_iterator.FaceIterator`,
         then this is the index of the occurence in the iterator.
+        In dual mode this value is then deducted from the maximal value of ``size_t``.
 
         If the face was constructed from
-        :meth:`sage:geometry.polyhedron.combinatorial_polyhedronn.base.CombinatorialPolyhedron.face_by_face_lattice_index`,
-        then this the index in the level set plus the number of lower dimension (or higher dimension).
+        :meth:`sage:geometry.polyhedron.combinatorial_polyhedron.base.CombinatorialPolyhedron.face_by_face_lattice_index`,
+        then this is the total number of faces plus the index in the level set plus the number of lower dimensional faces
+        (or higher dimensional faces in dual mode).
+        In dual mode this value is then deducted from the maximal value of ``size_t``.
 
         EXAMPLES::
 
@@ -296,6 +316,51 @@ cdef class CombinatorialFace(SageObject):
         """
         return self._hash_index
 
+    def __lt__(self, other):
+        r"""
+        Compare faces of the same polyhedron.
+
+        This is a helper function.
+        In order to construct a Hasse diagram (a digraph) with combinatorial faces,
+        we must define some order relation that is compatible with the Hasse diagram.
+
+        Any order relation compatible with ordering by dimension is suitable.
+        We use :meth:`__hash__` to define the relation.
+
+        EXAMPLES::
+
+            sage: P = polytopes.cube()
+            sage: C = CombinatorialPolyhedron(P)
+            sage: F1 = C.face_by_face_lattice_index(0)
+            sage: F2 = C.face_by_face_lattice_index(1)
+            sage: F1 < F2
+            True
+            sage: for i,j in Combinations(range(28), 2):
+            ....:     F1 = C.face_by_face_lattice_index(i)
+            ....:     F2 = C.face_by_face_lattice_index(j)
+            ....:     if F1.dim() != F2.dim():
+            ....:          assert (F1.dim() < F2.dim()) == (F1 < F2)
+
+            sage: P = polytopes.cross_polytope(3)
+            sage: C = CombinatorialPolyhedron(P)
+            sage: F1 = C.face_by_face_lattice_index(0)
+            sage: F2 = C.face_by_face_lattice_index(1)
+            sage: F1 < F2
+            True
+            sage: for i,j in Combinations(range(28), 2):
+            ....:     F1 = C.face_by_face_lattice_index(i)
+            ....:     F2 = C.face_by_face_lattice_index(j)
+            ....:     if F1.dim() != F2.dim():
+            ....:          assert (F1.dim() < F2.dim()) == (F1 < F2)
+        """
+        cdef CombinatorialFace other_face
+        if isinstance(other, CombinatorialFace):
+            other_face = other
+            if (self._initialized_from_face_lattice == other_face._initialized_from_face_lattice and
+                    self.atoms is other_face.atoms):
+                # They are faces of the same polyhedron obtained in the same way.
+                return hash(self) < hash(other)
+
     def dimension(self):
         r"""
         Return the dimension of the face.