Implement parallel f-vector for polytopes

[kliem]

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)

CC: @jplab @LaisRast @stumpc5 @tscrim

Component: geometry

Keywords: parallel f-vector

Author: Jonathan Kliem

Branch/Commit: 4c0a4ae

Reviewer: Travis Scrimshaw

Issue created by migration from https://trac.sagemath.org/ticket/31245

[kliem]

Description changed:

--- 
+++ 
@@ -2,7 +2,7 @@
 
 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.
+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.
 

[tscrim]

comment:3

The sig_check() outside the loop is unnecessary:

     cdef size_t j
+    sig_check()
     for j in range(n_faces):
+        sig_check()
         if (is_not_maximal_fused(new_faces, j, algorithm) or  # Step 2
                 is_contained_in_one_fused(new_faces.faces[j], visited_all, algorithm)):  # Step 3
             is_not_new_face[j] = True

Likewise here:

         face_list_t visited_all) nogil except -1:
 
     cdef size_t output
-    sig_on()
+    sig_check()
     if faces.polyhedron_is_simple:
         output = get_next_level_fused(faces, new_faces, visited_all, <simple> 0)
     else:
         output = get_next_level_fused(faces, new_faces, visited_all, <standard> 0)
-    sig_off()
     return output
 
 cdef inline size_t bit_rep_to_coatom_rep(face_t face, face_list_t coatoms, size_t *output):

Here you either need the sig_on/sig_off or some more targeted checking. I don't see the point in performing a single check here.

Otherwise LGTM.

[tscrim]

Reviewer: Travis Scrimshaw

[sagetrac-git]

Changed commit from d051f28 to 0dae7a0

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

0dae7a0

remove redundant sig_checks

[sagetrac-git]

Changed commit from 0dae7a0 to e3bdc61

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. New commits:

609bd92	do not delete all duplicates
5c6ccee	fix index
e3bdc61	pull in 31226 and fix merge conflict

[kliem]

Dependencies: #31226

[kliem]

comment:7

Replying to @tscrim:

The sig_check() outside the loop is unnecessary:

     cdef size_t j
+    sig_check()
     for j in range(n_faces):
+        sig_check()
         if (is_not_maximal_fused(new_faces, j, algorithm) or  # Step 2
                 is_contained_in_one_fused(new_faces.faces[j], visited_all, algorithm)):  # Step 3
             is_not_new_face[j] = True

Likewise here:

         face_list_t visited_all) nogil except -1:
 
     cdef size_t output
-    sig_on()
+    sig_check()
     if faces.polyhedron_is_simple:
         output = get_next_level_fused(faces, new_faces, visited_all, <simple> 0)
     else:
         output = get_next_level_fused(faces, new_faces, visited_all, <standard> 0)
-    sig_off()
     return output
 
 cdef inline size_t bit_rep_to_coatom_rep(face_t face, face_list_t coatoms, size_t *output):

Here you either need the sig_on/sig_off or some more targeted checking. I don't see the point in performing a single check here.

Otherwise LGTM.

I removed the redundant sig_checks.

As mentioned, sig_on, sig_off doesn't work well with parallelization. It crashes hard because sig_on, sig_off doesn't appear in correct order anymore.

Wrapping prange in sig_on/sig_off has a large performance penalty (8.45 s -> 9 s).

Currently there is one sig_check before checking that a face is inclusion maximal. One can add more sig checks, which wouldn't affect performance much for large examples, but for the small examples there is a penalty of about 50 percent.

I think the current solution should be fine: The largest thing I ever computed with this algorithm, which took about a month on 40 cores, had 11,665,781 vertices and 162 facets. In this case between the sig_checks A &~ B == 0 is performed for about 2 billion bits. Currently, without intrinsics, I get:

sage: B = Bitset((randint(0,2**20) for _ in range(2**19)))                                                                                                    
sage: %timeit B.issubset(B)                                                                                                                                   
6.18 µs ± 9.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
sage: %timeit B.intersection_update(B)                                                                                                                        
6.02 µs ± 2.72 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

On this rate performing this for 2 billion bits should take about 12 ms. And even if the example has many more coatoms, in the worst case (completely unrealistic, because the dimension is not 1 and more memory is needed for other dimensions and probably for other threads), there are bitwise operations performed twice for the entire RAM. Say, we have 1 TB RAM, this would then take 96 seconds. Okay, that is long, but in any realistic sceanrio it would only take about a second, even with 1 TB RAM.

To sum up, I think the current amount of sig_checks suffices for the next years until more than 1 TB RAM is common (and even then, it doubt that there is an application for such a large example, as it would take forever to compute).

[kliem]

comment:8

Rebased.

I didn't have this ticket here ready when writing #30445.

Eventually, one could also store edges/ridges or any sort of face in parallel.
This would also help with applications for simplicial complexes and lattice of flats of manifolds.

But that should be left for another time.

---

New commits:

6a2a9b4	fix merge conflict with develop in particual #30445
928ea60	remove redundant allocation