Check containment for combinatorial faces

[kliem]

We add a method is_subface to combinatorial faces.

Depends on #31819

CC: @jplab @mkoeppe @yuan-zhou

Component: geometry

Keywords: combinatorial polyhedron, subface

Author: Jonathan Kliem

Branch/Commit: ac5c327

Reviewer: Yuan Zhou, Matthias Koeppe

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

[kliem]

comment:3

It works fine in the intended case.

However, it should raise an error, when self and other or not from the same CombinatorialPolyhedron (NotImplementedError?) and it should consider the case that self is from a dual face iterator and `other isn't.

[sagetrac-git]

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

b4bfca6

only allow subface check for faces of identical combinatorial polyhedron

[sagetrac-git]

Changed commit from ce5ab6d to b4bfca6

[slel]

comment:6

The pyflakes plugin complains that facet is not defined in this line:

            if isinstance(v, PolyhedronFace) and facet.dim() == 0:

Did you mean v.dim()?

[sagetrac-git]

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

343b340

fix a variable name and add a doctest

[sagetrac-git]

Changed commit from b4bfca6 to 343b340

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

9a66e32	expose in Polyhedron_base
977c088	fix doctests
15e2c58	improved documentation
deca8e0	fix a variable name and add a doctest
be121ad	implement only subfaces/supfaces for face iterator
ee42e98	fix bug revealed by a_maximal_chain
05ec5b2	fix error message in doctest
6ddf228	fix the order of the coatoms when resetting the face iterator
f0a0afa	implement is_subface for combinatorial face
c632f56	only allow subface check for faces of identical combinatorial polyhedron

[sagetrac-git]

Changed commit from 343b340 to c632f56

[kliem]

comment:10

Replying to @slel:

The pyflakes plugin complains that facet is not defined in this line:

            if isinstance(v, PolyhedronFace) and facet.dim() == 0:

Did you mean v.dim()?

Fixed in #29683 and rebased.

[kliem]

comment:11

Thanks.

[kliem]

comment:12

Nees rebase.

[yuan-zhou]

comment:13

is_subface looks correct. However I don't understand why you distinguished two cases on whetherself._dual == other_face._dual. It seems to me that the ambient_V_indices containment check always works. Is it because using face_issubset in the first case is faster? Then, how about using only_subfaces?

sage: P = polytopes.cube()
sage: C = P.combinatorial_polyhedron()

is repeated in the doctest.

[sagetrac-git]

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

17720f3

delete duplication in doctest

[sagetrac-git]

Changed commit from c632f56 to 17720f3

[kliem]

comment:15

Replying to @yuan-zhou:

is_subface looks correct. However I don't understand why you distinguished two cases on whetherself._dual == other_face._dual. It seems to me that the ambient_V_indices containment check always works. Is it because using face_issubset in the first case is faster? Then, how about using only_subfaces?

Using the face as a bitset is much much faster than creating a list. Well there would be a faster method for the second case by comparing directly the array instead of the list, but I don't think this case is very important anyway.

Here is some speed comparison:

sage: P = polytopes.hypercube(14, backend='field')                                                                                                                                  
sage: C = CombinatorialPolyhedron(P)                                                                                                                                                
sage: it = C.face_iter(dual=False)                                                                                                                                                  
sage: f = next(it)                                                                                                                                                                  
sage: it.only_subfaces()                                                                                                                                                            
sage: f2 = next(it)                                                                                                                                                                 
sage: %timeit f2.is_subface(f)                                                                                                                                                      
80.2 ns ± 0.0408 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: %timeit f.is_subface(f2)                                                                                                                                                      
50.1 ns ± 0.0767 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: it2 = C.face_iter(dual=True)                                                                                                                                                  
sage: f3 = next(it2)                                                                                                                                                                
sage: %timeit f3.is_subface(f)                                                                                                                                                      
781 µs ± 444 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
sage: %timeit f.is_subface(f3)                                                                                                                                                      
708 µs ± 239 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
sage: while f3.dimension() < 12: 
....:     while not f3.is_subface(f): 
....:         f3 = next(it2) 
....:     it2.only_supfaces() 
....:     f3 = next(it2) 
....:                                                                                                                                                                               
sage: f3.dimension()                                                                                                                                                                
12
sage: %timeit f3.is_subface(f)                                                                                                                                                      
1.16 ms ± 619 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
sage: %timeit f.is_subface(f3)                                                                                                                                                      
1.06 ms ± 256 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Even the very difficult first case (f2 is indeed a subface of f) only takes a few nano seconds to check all 2^14 bits. In comparison the other case is really slow (but it is also not optimized at all).

sage: P = polytopes.cube()
sage: C = P.combinatorial_polyhedron()

is repeated in the doctest.

Fixed.

[yuan-zhou]

comment:16

Replying to @kliem:

Replying to @yuan-zhou:
In comparison the other case is really slow (but it is also not optimized at all).

If speed is a concern, then I'd recommend to check self.dimension() <= other.dimension() at least.

Does a face know about the original combinatorial polyhedron C? Is f2.is_subface(f1) through something like C.face_iter().find_face(f1).only_subfaces().find_face(f2) faster than comparing ambient_V_indices()?

[kliem]

comment:17

If f1 and f2 are in the same mode, then the current check is fast enough I think.

The other case is not really intended and I just have something that will work. However, if you want, I can do something more efficient. find_face is definitely not the best way of doing it, because it does way to much work. After all, I just have to get the coatom representation of f1 and compare it with the atom representation of f2 (the atoms of f1 are the coatoms of f2).

[yuan-zhou]

comment:18

Replying to @kliem:

After all, I just have to get the coatom representation of f1 and compare it with the atom representation of f2 (the atoms of f1 are the coatoms of f2).

If it's not too complicated to implement, this sounds like a great solution! I don't have a direct use case of "f1 and f2 not in the same mode", so I feel good about the ticket even if you won't implement the additional speed-up. However, I'd like to include the self.dimension() <= other.dimension() check, as mentioned above.

[yuan-zhou]

Reviewer: Yuan Zhou

[sagetrac-git]

Changed commit from 17720f3 to c0c8295

[sagetrac-git]

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

8ed706f	assume type to safe some time
c0c8295	optimized code

[kliem]

comment:20

Replying to @yuan-zhou:

Replying to @kliem:

After all, I just have to get the coatom representation of f1 and compare it with the atom representation of f2 (the atoms of f1 are the coatoms of f2).

If it's not too complicated to implement, this sounds like a great solution! I don't have a direct use case of "f1 and f2 not in the same mode", so I feel good about the ticket even if you won't implement the additional speed-up. However, I'd like to include the self.dimension() <= other.dimension() check, as mentioned above.

The check self.dimension() <= other.dimension() is no improvement either way, when self and other are in same mode.

Here are the new timings:

sage: P = polytopes.hypercube(14, backend='field')                                                                                                                                  
sage: C = CombinatorialPolyhedron(P)                                                                                                                                                
sage: it = C.face_iter(dual=False)                                                                                                                                                  
sage: f = next(it)                                                                                                                                                                  
sage: it.only_subfaces()                                                                                                                                                            
sage: f2 = next(it)                                                                                                                                                                 
sage: %timeit f.is_subface(f2)                                                                                                                                                      
42.2 ns ± 0.0401 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: %timeit f2.is_subface(f)                                                                                                                                                      
72.4 ns ± 0.0662 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: it2 = C.face_iter(dual=True)                                                                                                                                                  
sage: f3 = next(it2)                                                                                                                                                                
sage: %timeit f.is_subface(f3)                                                                                                                                                      
56.8 ns ± 0.0814 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: %timeit f3.is_subface(f)                                                                                                                                                      
78.1 µs ± 92.8 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
sage: while True: 
....:     while not f3.is_subface(f): 
....:         f3 = next(it2) 
....:     it2.only_supfaces() 
....:     if f3.dimension() == 12: 
....:         break 
....:     f3 = next(it2) 
....:                                                                                                                                                                               
sage: %timeit f.is_subface(f3)                                                                                                                                                      
58.5 ns ± 0.0699 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: %timeit f3.is_subface(f)                                                                                                                                                      
78.3 µs ± 46.5 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
sage: f3.dimension()                                                                                                                                                                
12
sage: f.is_subface(f3)                                                                                                                                                              
False
sage: f3.is_subface(f)                                                                                                                                                              
True

[kliem]

comment:21

Btw, currently face.coatom_rep etc. is reinitialized. I expect things to be much faster on repeated calls, if I change the code to realize this.

However, this might be the more honest timing, as it reflects more realistic the amount of time the first call needs.

[kliem]

comment:22

Caching atom_rep and coatom_rep (actually not cashing, but only remembering that it has already been computed), I get the following timings now:

sage: P = polytopes.hypercube(14, backend='field')                                                                                                                                  
sage: C = CombinatorialPolyhedron(P)                                                                                                                                                
sage: it = C.face_iter(dual=False)                                                                                                                                                  
sage: f = next(it)                                                                                                                                                                   
sage: it2 = C.face_iter(dual=True)                                                                                                                                                  
sage: f3 = next(it2)                                                                                                                                                                
sage: %timeit f.is_subface(f3)                                                                                                                                                      
56.7 ns ± 0.0563 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: %timeit f3.is_subface(f)                                                                                                                                                      
7.12 µs ± 68.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
sage: while True: 
....:     while not f3.is_subface(f): 
....:         f3 = next(it2) 
....:     it2.only_supfaces() 
....:     if f3.dimension() == 12: 
....:         break 
....:     f3 = next(it2) 
....:                                                                                                                                                                               
sage: %timeit f.is_subface(f3)                                                                                                                                                      
58.1 ns ± 0.0822 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)
sage: %timeit f3.is_subface(f)                                                                                                                                                      
5.34 µs ± 2.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

[sagetrac-git]

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

5de63be

cache atom_rep and coatom_rep of combinatorial face

[sagetrac-git]

Changed commit from c0c8295 to 5de63be

[mkoeppe]

Changed branch from u/gh-kliem/is_subface_for_combinatorial_face to u/mkoeppe/is_subface_for_combinatorial_face

[mkoeppe]

Changed commit from 5de63be to ac5c327

[mkoeppe]

comment:25

I am guessing that this ticket "needs review"

---

Last 10 new commits:

70a8791	Merge branch 'u/gh-kliem/check_openmp_at_configuration' of git://trac.sagemath.org/sage into u/gh-kliem/first_parallel_version_of_face_iterator_reb2
fdbe95f	use cython.parallel.threadid instead of omp_get_thread_num
4ae6966	Merge tag '9.4.beta0' into t/31245/first_parallel_version_of_face_iterator_reb2
bfb4efb	Merge branch 'u/mkoeppe/first_parallel_version_of_face_iterator_reb2' of git://trac.sagemath.org/sage into u/mkoeppe/first_parallel_version_of_face_iterator_reb3
4c0a4ae	initialize do_f_vector
542baee	merge in #31245
3e0229f	Merge tag '9.4.beta3' into t/31834/public/31834-reb2
8190917	Merge #31834
7f7e630	Merge #29683
ac5c327	Merge #31819

[mkoeppe]

Changed reviewer from Yuan Zhou to Yuan Zhou, Matthias Koeppe

[mkoeppe]

comment:26

I have not followed the discussion on timing, but the code looks fine to me, and passes tests.

[kliem]

comment:27

Thank you.

[vbraun]

Changed branch from u/mkoeppe/is_subface_for_combinatorial_face to ac5c327