Improve face iterator for almost simple/simplicial polytopes

[kliem]

#30040 improved the face iterator for simple/simplicial polytopes.

Actually, the same trick can be applied for simple faces. (Note that the iterator treats simplicial polytopes dually).

If a face is simple, then any intersections of it's facets is one of it's ridges or empty.

As this helps only sporadically, we cannot afford to compute for each face if it is simple. However, a cheap certificate of a face being simple is that it only contains "simple vertices" (with vertex figure a simplex).

Now, at initialization the face iterator computes those "simple vertices" and checks inclusion of each face to determine if it is simple (if not at least 10 percent of the vertices we skip this check). If a face is simple, we can skip the costly inclusion checks.

Before this ticket:

sage: P = polytopes.permutahedron(7)                                                                                                                                                                                                                                                                                                                                       
sage: P1 = P.stack(next(P.face_generator(1)))                                                                                                                                                                                                                                                                                                                              
sage: C = CombinatorialPolyhedron(P1)                                                                                                                                                                                                                                                                                                                                      
sage: %time C.f_vector()                                                                                                                                                                                                                                                                                                                                                   
CPU times: user 257 ms, sys: 7 µs, total: 257 ms
Wall time: 257 ms
(1, 5041, 16251, 19761, 11144, 2860, 267, 1)

sage: P = polytopes.associahedron(['A',7])                                                                                                                                                                                                                                                                                                                                 
sage: P1 = P.stack(next(P.face_generator(1)))                                                                                                                                                                                                                                                                                                                              
sage: C = CombinatorialPolyhedron(P1)                                                                                                                                                                                                                                                                                                                                      
sage: %time C.f_vector()                                                                                                                                                                                                                                                                                                                                                   
CPU times: user 13.5 ms, sys: 0 ns, total: 13.5 ms
Wall time: 13.4 ms
(1, 1431, 5977, 10177, 9040, 4440, 1163, 134, 1)
sage: P = polytopes.hypercube(12)                                                                                                                                                                                                                                                                                                                                          
sage: P1 = P.stack(next(P.face_generator(1)))                                                                                                                                                                                                                                                                                                                              
sage: C = CombinatorialPolyhedron(P1)                                                                                                                                                                                                                                                                                                                                      
sage: %time C.f_vector()                                                                                                                                                                                                                                                                                                                                                   
CPU times: user 419 ms, sys: 17 µs, total: 419 ms
Wall time: 419 ms
(1, 4097, 28669, 92125, 180037, 238755, 226776, 158466, 82170, 31350, 8525, 1529, 145, 1)

sage: P = polytopes.permutahedron(7)                                                                                                                                                                                                                                                                                                                                       
sage: P1 = P.stack(next(P.face_generator(1)))                                                                                                                                                                                                                                                                                                                              
sage: C = CombinatorialPolyhedron(P1)                                                                                                                                                                                                                                                                                                                                      
sage: %time C.f_vector()                                                                                                                                                                                                                                                                                                                                                   
CPU times: user 79.5 ms, sys: 0 ns, total: 79.5 ms
Wall time: 79.3 ms
(1, 5041, 16251, 19761, 11144, 2860, 267, 1)

sage: P = polytopes.associahedron(['A',7])                                                                                                                                                                                                                                                                                                                                 
sage: P1 = P.stack(next(P.face_generator(1)))                                                                                                                                                                                                                                                                                                                              
sage: C = CombinatorialPolyhedron(P1)                                                                                                                                                                                                                                                                                                                                      
sage: %time C.f_vector()                                                                                                                                                                                                                                                                                                                                                   
CPU times: user 9.85 ms, sys: 0 ns, total: 9.85 ms
Wall time: 9.72 ms
(1, 1431, 5977, 10177, 9040, 4440, 1163, 134, 1)

sage: P = polytopes.hypercube(12)                                                                                                                                                                                                                                                                                                                                          
sage: P1 = P.stack(next(P.face_generator(1)))                                                                                                                                                                                                                                                                                                                              
sage: C = CombinatorialPolyhedron(P1)                                                                                                                                                                                                                                                                                                                                      
sage: %time C.f_vector()                                                                                                                                                                                                                                                                                                                                                   
CPU times: user 391 ms, sys: 0 ns, total: 391 ms
Wall time: 391 ms
(1, 4097, 28669, 92125, 180037, 238755, 226776, 158466, 82170, 31350, 8525, 1529, 145, 1)

Note that for the last instance, there isn't a significant amount of simple vertices so the old approach is chosen.

Depends on #30040
Depends on #30435

CC: @jplab @LaisRast @tscrim

Component: geometry

Keywords: simple, simplicial, face iterator, polytope

Author: Jonathan Kliem

Branch/Commit: u/gh-kliem/face_iterator_simple_face_2 @ d3c9384

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

[kliem]

Last 10 new commits:

0297853	improvements in documentation
7ffac00	changes in 30429
01b4154	typo
34a74b8	simpler algorithm for many simple faces
c9165b0	sort the faces according to being simple
b9dc7e1	small improvements
5dbe710	improved count atoms
e88c28c	small improvements
a50be39	compile warning in doctest
e0d117a	two bugs

[kliem]

Branch: u/gh-kliem/face_iterator_simple_face

[kliem]

Commit: e0d117a

[sagetrac-git]

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

59f6a12

improved documentation

[sagetrac-git]

Changed commit from e0d117a to 59f6a12

[kliem]

Changed dependencies from #30040, #30045 to #30040, #30435

[kliem]

comment:7

A remark to cythonizing and compiling:

• Current version:

+            foo = self.atoms.data
+            atoms_face_length = self.atoms.face_length
+
+            counter = 0
+            for j in range(self.atoms.n_faces):
+                if count_atoms(foo[j], atoms_face_length) == self.structure.dimension:
+                    counter += 1
+                    self.structure.simple_vertices[j//64] |= vertex_to_bit_dictionary(j % 64)

• Much slower:

+            counter = 0
+            for j in range(self.atoms.n_faces):
+                if count_atoms(self.atoms.data[j], self.atoms.face_length) == self.structure.dimension:
+                    counter += 1
+                    self.structure.simple_vertices[j//64] |= vertex_to_bit_dictionary(j % 64)

• As fast as the current version:

+            counter = 0
+            for j in range(self.atoms.n_faces):
+                if count_atoms(self.atoms.data[j], self.atoms.face_length) == self.structure.dimension:
+                    counter += 1
+                    self.structure.simple_vertices[j//64] |= vertex_to_bit_dictionary(j % 64)
+
+            for j in range(self.atoms.n_faces):
+                if count_atoms(self.atoms.data[j], self.atoms.face_length) == self.structure.dimension:
+                    self.structure.simple_vertices[j//64] |= vertex_to_bit_dictionary(j % 64)

I guess looking up self.atoms takes significant amount of time. What is suprising is that when you add a useless loop, things speed up. Either by cythonizing or compiling (I didn't check) it is figured, that we need the properties of self.atoms more often and this is optimized then.

Perhaps even more suprising. When you add

# cython: profile=True
# cython: linetrace=True
# distutils: define_macros=CYTHON_TRACE_NOGIL=1
# distutils: define_macros=CYTHON_TRACE=1

to the beginning of base.pyx the second version is just as fast. Makes it hard to profile where the difference in performance is, when it goes away when profiling.

[sagetrac-git]

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

20a39b6	outsource inclusion maximal
b672fca	removed redundant function
f62e770	merge in #30458
accea52	missing }
56ba0e0	merge in #30040

[sagetrac-git]

Changed commit from 59f6a12 to 56ba0e0

[kliem]

Changed commit from 56ba0e0 to 79a8666

[kliem]

comment:9

Changes from #30040.

---

Last 10 new commits:

ed6e966	improvements in documentation
bf91483	improved count atoms
7bab490	Merge branch 'u/gh-kliem/improved_count_atoms' of git://trac.sagemath.org/sage into u/gh-kliem/face_iterator_simple_face_2
d9aca3b	simpler algorithm for many simple faces
3550857	sort the faces according to being simple
19b293c	small improvements
b6e4659	small improvements
d4a9fd2	compile warning in doctest
de36d26	two bugs
79a8666	improved documentation

[kliem]

Changed branch from u/gh-kliem/face_iterator_simple_face to u/gh-kliem/face_iterator_simple_face_2

[mkoeppe]

comment:10

needs rebase

[sagetrac-git]

Changed commit from 79a8666 to d3c9384

[sagetrac-git]

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

d3c9384

Merge branch 'develop' of git://trac.sagemath.org/sage into u/gh-kliem/face_iterator_simple_face_2

[kliem]

comment:13

I want to redesign some of the setup before making everything more complicated.

More precisely, creating strucutures face_struct and faces_list_struct that take care of the details. Then this overly long argument list of get_next_level will just reduce to three arguments or so. This would also make future changes easier including the transition to bitsets.pxi.

[mkoeppe]

comment:15

Setting new milestone based on a cursory review of ticket status, priority, and last modification date.

[mkoeppe]

comment:16

Setting a new milestone for this ticket based on a cursory review.