Trac #31245: Implement parallel f-vector for polytopes

This ticket parallelizes the f-vector for polytopes.

Each thread has its private structure with which it does partial jobs.
Depending on the parallelization depth, there is one job per face of
fixed codimension (usually 1,2 or 3). After everything has finished, the
partial f-vectors will be added.

Actually, every face is visited and thus the code could be modified in
the future, to explore other properties of faces then just the
dimension. The parallelization seems to work well with at least 40
threads (for computations taking long enough such that this pays off,
see https://arxiv.org/pdf/1905.01945.pdf).

Also the algorithm does work in other situations (simplicial complex,
lattice of flats of a matroid) and this parallel structure could be used
for this as well.

On the downside, `sig_on()/sig_off()` doesn't work with with multiple
threads and has to be replaced by a simple `sig_check()`. Also raising
errors in parallel code results in terrible slow down. Hence the errors
are replaced by returning the exception value. In case of a bug there
won't be any traceback anymore, but at least also no segmenation fault.

Before:

{{{
sage: P = P.base_extend(QQ)
sage: P = P.base_extend(QQ)
sage: Q = P.join(P.polar(in_affine_span=True))
sage: C = CombinatorialPolyhedron(Q)
sage: %time C.f_vector()
CPU times: user 111 ms, sys: 186 µs, total: 111 ms
Wall time: 111 ms
(1, 150, 3990, 25590, 69450, 95402, 69450, 25590, 3990, 150, 1)

sage: P = polytopes.Birkhoff_polytope(5)
sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector()
CPU times: user 584 ms, sys: 25 µs, total: 584 ms
Wall time: 583 ms
(1, 120, 5040, 50250, 233400, 631700, 1113700, 1367040, 1220550, 817150,
419225, 167200, 52120, 12600, 2300, 300, 25, 1)

# Using the <simple> version of the algorithm.
sage: P = polytopes.associahedron(['A', 11], backend='normaliz')
sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector()
CPU times: user 37.9 s, sys: 17.2 ms, total: 37.9 s
Wall time: 37.9 s
(1, 208012, 1144066, 2735810, 3730650, 3197700, 1790712, 659736, 157080,
23100, 1925, 77, 1)
}}}

After (machine has 4 cores):

{{{
sage: P = polytopes.permutahedron(5)
sage: P = P.base_extend(QQ)
sage: Q = P.join(P.polar(in_affine_span=True))
sage: C = CombinatorialPolyhedron(Q)
sage: %time C.f_vector(num_threads=1)
CPU times: user 107 ms, sys: 0 ns, total: 107 ms
Wall time: 107 ms
(1, 150, 3990, 25590, 69450, 95402, 69450, 25590, 3990, 150, 1)

sage: C = CombinatorialPolyhedron(Q)
sage: %time C.f_vector(num_threads=2)
CPU times: user 108 ms, sys: 0 ns, total: 108 ms
Wall time: 55.6 ms
(1, 150, 3990, 25590, 69450, 95402, 69450, 25590, 3990, 150, 1)

sage: C = CombinatorialPolyhedron(Q)
sage: %time C.f_vector(num_threads=4)
CPU times: user 147 ms, sys: 52 µs, total: 147 ms
Wall time: 38.6 ms
(1, 150, 3990, 25590, 69450, 95402, 69450, 25590, 3990, 150, 1)

sage: C = CombinatorialPolyhedron(Q)
sage: %time C.f_vector(num_threads=8)
CPU times: user 236 ms, sys: 0 ns, total: 236 ms
Wall time: 31.2 ms
(1, 150, 3990, 25590, 69450, 95402, 69450, 25590, 3990, 150, 1)

sage: P = polytopes.Birkhoff_polytope(5)
sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=1)
CPU times: user 354 ms, sys: 0 ns, total: 354 ms
Wall time: 354 ms
(1, 120, 5040, 50250, 233400, 631700, 1113700, 1367040, 1220550, 817150,
419225, 167200, 52120, 12600, 2300, 300, 25, 1)

sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=2)
CPU times: user 363 ms, sys: 0 ns, total: 363 ms
Wall time: 181 ms
(1, 120, 5040, 50250, 233400, 631700, 1113700, 1367040, 1220550, 817150,
419225, 167200, 52120, 12600, 2300, 300, 25, 1)

sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=4)
CPU times: user 459 ms, sys: 0 ns, total: 459 ms
Wall time: 117 ms
(1, 120, 5040, 50250, 233400, 631700, 1113700, 1367040, 1220550, 817150,
419225, 167200, 52120, 12600, 2300, 300, 25, 1)

sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=8)
CPU times: user 776 ms, sys: 154 µs, total: 776 ms
Wall time: 103 ms
(1, 120, 5040, 50250, 233400, 631700, 1113700, 1367040, 1220550, 817150,
419225, 167200, 52120, 12600, 2300, 300, 25, 1)

# Using the <simple> version of the algorithm.
sage: P = polytopes.associahedron(['A', 11], backend='normaliz')
sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=1)
CPU times: user 33.5 s, sys: 0 ns, total: 33.5 s
Wall time: 33.5 s
(1, 208012, 1144066, 2735810, 3730650, 3197700, 1790712, 659736, 157080,
23100, 1925, 77, 1)

sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=2)
CPU times: user 34.4 s, sys: 3.49 ms, total: 34.4 s
Wall time: 17.2 s
(1, 208012, 1144066, 2735810, 3730650, 3197700, 1790712, 659736, 157080,
23100, 1925, 77, 1)

sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=4)
CPU times: user 35.9 s, sys: 15.5 ms, total: 35.9 s
Wall time: 9 s
(1, 208012, 1144066, 2735810, 3730650, 3197700, 1790712, 659736, 157080,
23100, 1925, 77, 1)

sage: C = CombinatorialPolyhedron(P)
sage: %time C.f_vector(num_threads=8)
CPU times: user 1min 6s, sys: 31.3 ms, total: 1min 6s
Wall time: 8.44 s
(1, 208012, 1144066, 2735810, 3730650, 3197700, 1790712, 659736, 157080,
23100, 1925, 77, 1)
}}}

URL: https://trac.sagemath.org/31245
Reported by: gh-kliem
Ticket author(s): Jonathan Kliem
Reviewer(s): Travis Scrimshaw

src/sage/geometry/polyhedron/base.py

@@ -6994,10 +6994,18 @@ def _test_combinatorial_face_as_combinatorial_polyhedron(self, tester=None, **op
             tester.assertTrue(P.combinatorial_polyhedron().vertex_facet_graph().is_isomorphic(D2.vertex_facet_graph()))
 
     @cached_method(do_pickle=True)
-    def f_vector(self):
+    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
@@ -7053,7 +7061,7 @@ def f_vector(self):
             sage: Q.f_vector.is_in_cache()
             True
         """
-        return self.combinatorial_polyhedron().f_vector()
+        return self.combinatorial_polyhedron().f_vector(num_threads, parallelization_depth)
 
     def flag_f_vector(self, *args):
         r"""

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

@@ -62,7 +62,7 @@ cdef class CombinatorialPolyhedron(SageObject):
     cdef tuple _mem_tuple
 
     cdef FaceIterator _face_iter(self, bint dual, int dimension)
-    cdef int _compute_f_vector(self, bint compute_edges=*, given_dual=*) except -1
+    cdef int _compute_f_vector(self, size_t num_threads, size_t parallelization_depth) except -1
 
     cdef inline int _compute_edges(self, dual) except -1:
         return self._compute_edges_or_ridges(dual, True)
@@ -71,7 +71,10 @@ cdef class CombinatorialPolyhedron(SageObject):
         return self._compute_edges_or_ridges(dual, False)
 
     cdef int _compute_edges_or_ridges(self, bint dual, bint do_edges) except -1
-    cdef size_t _compute_edges_or_ridges_with_iterator(self, FaceIterator face_iter, bint do_atom_rep, size_t ***edges_pt, size_t *counter_pt, size_t *current_length_pt, MemoryAllocator mem) except -1
+    cdef size_t _compute_edges_or_ridges_with_iterator(
+            self, FaceIterator face_iter, const bint do_atom_rep, const bint do_f_vector,
+            size_t ***edges_pt, size_t *counter_pt, size_t *current_length_pt,
+            size_t* f_vector, MemoryAllocator mem) except -1
     cdef int _compute_face_lattice_incidences(self) except -1
 
     cdef inline int _set_edge(self, size_t a, size_t b, size_t ***edges_pt, size_t *counter_pt, size_t *current_length_pt, MemoryAllocator mem) except -1

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

@@ -104,6 +104,7 @@ from cysignals.signals                 cimport sig_check, sig_block, sig_unblock
 from sage.matrix.matrix_integer_dense  cimport Matrix_integer_dense
 
 from .face_data_structure cimport face_len_atoms, face_init
+from .face_iterator cimport iter_t, parallel_f_vector
 
 cdef extern from "Python.h":
     int unlikely(int) nogil  # Defined by Cython
@@ -1572,13 +1573,20 @@ cdef class CombinatorialPolyhedron(SageObject):
         return DiGraph([vertices, edges], format='vertices_and_edges', immutable=True)
 
     @cached_method
-    def f_vector(self):
+    def f_vector(self, num_threads=None, parallelization_depth=None):
         r"""
         Compute the ``f_vector`` of the polyhedron.
 
         The ``f_vector`` contains the number of faces of dimension `k`
         for each `k` in ``range(-1, self.dimension() + 1)``.
 
+        INPUT:
+
+        - ``num_threads`` -- integer (optional); specify the number of threads
+
+        - ``parallelization_depth`` -- integer (optional); specify
+          how deep in the lattice the parallelization is done
+
         .. NOTE::
 
             To obtain edges and/or ridges as well, first do so. This might
@@ -1596,13 +1604,35 @@ cdef class CombinatorialPolyhedron(SageObject):
             sage: C.f_vector()
             (1, 10, 45, 120, 185, 150, 50, 1)
 
+        Using two threads::
+
+            sage: P = polytopes.permutahedron(5)
+            sage: C = CombinatorialPolyhedron(P)
+            sage: C.f_vector(num_threads=2)
+            (1, 120, 240, 150, 30, 1)
+
         TESTS::
 
             sage: type(C.f_vector())
             <type 'sage.modules.vector_integer_dense.Vector_integer_dense'>
         """
+        if num_threads is None:
+            from sage.parallel.ncpus import ncpus
+            num_threads = ncpus()
+
+        if parallelization_depth is None:
+            # Setting some reasonable defaults.
+            if num_threads == 0:
+                parallelization_depth = 0
+            elif num_threads <= 3:
+                parallelization_depth = 1
+            elif num_threads <= 8:
+                parallelization_depth = 2
+            else:
+                parallelization_depth = 3
+
         if not self._f_vector:
-            self._compute_f_vector()
+            self._compute_f_vector(num_threads, parallelization_depth)
         if not self._f_vector:
             raise ValueError("could not determine f_vector")
         from sage.modules.free_module_element import vector
@@ -2921,68 +2951,49 @@ cdef class CombinatorialPolyhedron(SageObject):
 
     # Internal methods.
 
-    cdef int _compute_f_vector(self, bint compute_edges=False, given_dual=None) except -1:
+    cdef int _compute_f_vector(self, size_t num_threads, size_t parallelization_depth) except -1:
         r"""
         Compute the ``f_vector`` of the polyhedron.
 
-        If ``compute_edges`` computes the edges in non-dual mode as well.
-        In dual mode the ridges.
-
-        See :meth:`f_vector` and :meth:`_compute_edges`.
+        See :meth:`f_vector`.
         """
         if self._f_vector:
             return 0  # There is no need to recompute the f_vector.
 
+        cdef int dim = self.dimension()
+        cdef int d  # dimension of the current face of the iterator
+        cdef MemoryAllocator mem = MemoryAllocator()
+
+        if num_threads == 0:
+            # No need to complain.
+            num_threads = 1
+
+        if parallelization_depth > dim - 1:
+            # Is a very bad choice anyway, but prevent segmenation faults.
+            parallelization_depth = dim - 1
+
         cdef bint dual
         if not self.is_bounded() or self.n_facets() <= self.n_Vrepresentation():
             # In this case the non-dual approach is faster..
             dual = False
         else:
             # In this case the dual approach is faster.
             dual = True
-        if given_dual is not None:
-            dual = given_dual
-
-
-        cdef FaceIterator face_iter = self._face_iter(dual, -2)
 
-        cdef int dim = self.dimension()
-        cdef int d  # dimension of the current face of the iterator
-        cdef MemoryAllocator mem = MemoryAllocator()
+        cdef FaceIterator face_iter
+        cdef iter_t* structs = <iter_t*> mem.allocarray(num_threads, sizeof(iter_t))
+        cdef size_t i
 
-        # In case we compute the edges as well.
-        cdef size_t **edges
-        cdef size_t counter = 0         # the number of edges so far
-        cdef size_t current_length      # dynamically enlarge **edges
-        cdef size_t a,b                # vertices of an edge
-        if compute_edges:
-            edges = <size_t**> mem.malloc(sizeof(size_t**))
-            current_length = 1
+        # For each thread an independent structure.
+        face_iters = [self._face_iter(dual, -2) for _ in range(num_threads)]
+        for i in range(num_threads):
+            face_iter = face_iters[i]
+            structs[i][0] = face_iter.structure[0]
 
         # Initialize ``f_vector``.
         cdef size_t *f_vector = <size_t *> mem.calloc((dim + 2), sizeof(size_t))
-        f_vector[0] = 1         # Face iterator will only visit proper faces.
-        f_vector[dim + 1] = 1   # Face iterator will only visit proper faces.
-
-        # For each face in the iterator, add `1` to the corresponding entry in
-        # ``f_vector``.
-        if self.n_facets() > 0 and dim > 0:
-            d = face_iter.next_dimension()
-            while (d < dim):
-                sig_check()
-                f_vector[d+1] += 1
-
-                if compute_edges and d == 1:
-                    # If it is an edge.
-
-                    # Set up face_iter.atom_rep
-                    face_iter.set_atom_rep()
 
-                    # Copy the information.
-                    a = face_iter.structure.atom_rep[0]
-                    b = face_iter.structure.atom_rep[1]
-                    self._set_edge(a, b, &edges, &counter, &current_length, mem)
-                d = face_iter.next_dimension()
+        parallel_f_vector(structs, num_threads, parallelization_depth, f_vector)
 
         # Copy ``f_vector``.
         if dual:
@@ -3003,21 +3014,6 @@ cdef class CombinatorialPolyhedron(SageObject):
 
             self._f_vector = tuple(smallInteger(f_vector[i]) for i in range(dim+2))
 
-        if compute_edges:
-            # Success, copy the data to ``CombinatorialPolyhedron``.
-            if dual:
-                sig_block()
-                self._n_ridges = counter
-                self._ridges = edges
-                self._mem_tuple += (mem,)
-                sig_unblock()
-            else:
-                sig_block()
-                self._n_edges = counter
-                self._edges = edges
-                self._mem_tuple += (mem,)
-                sig_unblock()
-
     cdef int _compute_edges_or_ridges(self, bint dual, bint do_edges) except -1:
         r"""
         Compute the edges of the polyhedron if ``edges`` else the ridges.
@@ -3041,6 +3037,9 @@ cdef class CombinatorialPolyhedron(SageObject):
         cdef size_t current_length = 1  # dynamically enlarge **edges
         cdef int output_dim_init = 1 if do_edges else dim - 2
 
+        cdef bint do_f_vector = False
+        cdef size_t* f_vector
+
         if dim == 1 and (do_edges or self.n_facets() > 1):
             # In this case there is an edge/ridge, but its not a proper face.
             self._set_edge(0, 1, &edges, &counter, &current_length, mem)
@@ -3050,18 +3049,46 @@ cdef class CombinatorialPolyhedron(SageObject):
             # Prevent an error when calling the face iterator.
             pass
 
-        elif not self._f_vector and ((dual ^ do_edges)):
-            # While doing edges in non-dual mode or ridges in dual-mode
-            # one might as well do the f-vector.
-            return self._compute_f_vector(compute_edges=True, given_dual=dual)
-
         else:
-            # Only compute the edges/ridges.
-            face_iter = self._face_iter(dual, output_dim_init)
-
-            self._compute_edges_or_ridges_with_iterator(face_iter, (dual ^ do_edges), &edges, &counter, &current_length, mem)
+            if not self._f_vector and ((dual ^ do_edges)):
+                # While doing edges in non-dual mode or ridges in dual-mode
+                # one might as well do the f-vector.
+                do_f_vector = True
+                # Initialize ``f_vector``.
+                f_vector = <size_t *> mem.calloc((dim + 2), sizeof(size_t))
+                f_vector[0] = 1
+                f_vector[dim + 1] = 1
+                face_iter = self._face_iter(dual, -2)
+            else:
+                do_f_vector = False
+                face_iter = self._face_iter(dual, output_dim_init)
+            self._compute_edges_or_ridges_with_iterator(face_iter, (dual ^ do_edges), do_f_vector,
+                                                        &edges, &counter, &current_length,
+                                                        f_vector, mem)
 
         # Success, copy the data to ``CombinatorialPolyhedron``.
+
+        # Copy ``f_vector``.
+        if do_f_vector:
+            if dual:
+                if dim > 1 and f_vector[1] < self.n_facets():
+                    # The input seemed to be wrong.
+                    raise ValueError("not all facets are joins of vertices")
+
+                # We have computed the ``f_vector`` of the dual.
+                # Reverse it:
+                self._f_vector = \
+                    tuple(smallInteger(f_vector[dim+1-i]) for i in range(dim+2))
+
+            else:
+                if self.is_bounded() and dim > 1 \
+                        and f_vector[1] < self.n_Vrepresentation() - len(self.far_face_tuple()):
+                    # The input seemed to be wrong.
+                    raise ValueError("not all vertices are intersections of facets")
+
+                self._f_vector = tuple(smallInteger(f_vector[i]) for i in range(dim+2))
+
+        # Copy the edge or ridges.
         if do_edges:
             sig_block()
             self._n_edges = counter
@@ -3076,32 +3103,45 @@ cdef class CombinatorialPolyhedron(SageObject):
             sig_unblock()
 
     cdef size_t _compute_edges_or_ridges_with_iterator(
-            self, FaceIterator face_iter, bint do_atom_rep, size_t ***edges_pt,
-            size_t *counter_pt, size_t *current_length_pt,
-            MemoryAllocator mem) except -1:
+            self, FaceIterator face_iter, const bint do_atom_rep, const bint do_f_vector,
+            size_t ***edges_pt, size_t *counter_pt, size_t *current_length_pt,
+            size_t* f_vector, MemoryAllocator mem) except -1:
         r"""
         See :meth:`CombinatorialPolyhedron._compute_edges`.
         """
         cdef size_t a,b                # facets of an edge
         cdef int dim = self.dimension()
+
+        # The dimension in which to record the edges or ridges.
         cdef output_dimension = 1 if do_atom_rep else dim - 2
 
-        while face_iter.next_dimension() == output_dimension:
-            if do_atom_rep:
-                # Set up face_iter.atom_rep
-                face_iter.set_atom_rep()
+        cdef int d = face_iter.next_dimension()
+        while d < dim:
+            sig_check()
+            if do_f_vector:
+                f_vector[d + 1] += 1
+
+            # If ``not do_f_vector`` the iterator is set up
+            # for ``output_dimension`` and
+            # ``d < dim`` implies
+            # ``d == ouput_dimension``.
+            if not do_f_vector or d == output_dimension:
+                if do_atom_rep:
+                    # Set up face_iter.atom_rep
+                    face_iter.set_atom_rep()
 
-                # Copy the information.
-                a = face_iter.structure.atom_rep[0]
-                b = face_iter.structure.atom_rep[1]
-            else:
-                # Set up face_iter.coatom_rep
-                face_iter.set_coatom_rep()
+                    # Copy the information.
+                    a = face_iter.structure.atom_rep[0]
+                    b = face_iter.structure.atom_rep[1]
+                else:
+                    # Set up face_iter.coatom_rep
+                    face_iter.set_coatom_rep()
 
-                # Copy the information.
-                a = face_iter.structure.coatom_rep[0]
-                b = face_iter.structure.coatom_rep[1]
-            self._set_edge(a, b, edges_pt, counter_pt, current_length_pt, mem)
+                    # Copy the information.
+                    a = face_iter.structure.coatom_rep[0]
+                    b = face_iter.structure.coatom_rep[1]
+                self._set_edge(a, b, edges_pt, counter_pt, current_length_pt, mem)
+            d = face_iter.next_dimension()
 
     cdef int _compute_face_lattice_incidences(self) except -1:
         r"""

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

@@ -35,7 +35,7 @@ cdef struct iter_s:
     face_list_t* new_faces
 
     # After having visited a face completely, we want to add it to ``visited_all``.
-    # ``first_dim[i]`` will indicate, whether there is one more face in
+    # ``first_time[i]`` will indicate, wether there is one more face in
     # ``newfaces[i]`` then ``n_newfaces[i]`` suggests
     # that has to be added to ``visited_all``.
     # If ``first_time[i] == False``, we still need to
@@ -45,6 +45,7 @@ cdef struct iter_s:
     # The number of elements in newfaces[current_dimension],
     # that have not been visited yet.
     size_t yet_to_visit
+    size_t n_coatoms
 
 ctypedef iter_s iter_t[1]
 
@@ -80,8 +81,10 @@ cdef class FaceIterator_geom(FaceIterator_base):
     cdef object _requested_dim  # Dimension requested on init.
     cdef readonly object P      # The original polyhedron.
 
+cdef int parallel_f_vector(iter_t* structures, size_t num_threads, size_t parallelization_depth, size_t *f_vector) except -1
+
 # Nogil definitions of crucial functions.
 
-cdef int next_dimension(iter_t structure) nogil except -1
+cdef int next_dimension(iter_t structure, size_t parallelization_depth=?) nogil except -1
 cdef int next_face_loop(iter_t structure) nogil except -1
 cdef size_t n_atom_rep(iter_t structure) nogil except -1