Optimize computing points of lattice polytopes

[Andrey Novoseltsev]

A follow up to #22391. My original problem was too slow computation of zero via

Delta3 = LatticePolytope([(1,0,0), (0,1,0), (0,0,1), (-1,-4,-6)])
sum(len(e.interior_points()) * len(e.dual().interior_points()) for e in Delta3.edges())

Sage 7.6.beta4 gives

• timeit: 5 loops, best of 3: 1.73 s per loop
• profiler: 71844 function calls (71620 primitive calls) in 1.841 seconds

With #22391

• timeit: 5 loops, best of 3: 451 ms per loop
• profiler: 46502 function calls (46434 primitive calls) in 0.570 seconds

With this one

• timeit: 25 loops, best of 3: 11.3 ms per loop
• profiler: 7793 function calls (7747 primitive calls) in 0.031 seconds

As a less contrived example, it now takes 0.41s to get ReflexivePolytopes(3) (cached afterwards) instead of 7.61s. Part of the speed up is due to omitting some precomputation, but it is possible to do so because of speed up on demand.

Next goal is not to rely on PALP for computing points of standalone polytopes at all, although for large sets of polytopes it still may be the fastest option.

[Andrey Novoseltsev]

Had to fix a doctest, results of sage -t --long src/sage/schemes/toric/fano_variety.py

BEFORE (7.6.beta4)
Total time for all tests: 20.5 seconds
    cpu time: 13.5 seconds
    cumulative wall time: 20.4 seconds

AFTER
Total time for all tests: 5.0 seconds
    cpu time: 3.8 seconds
    cumulative wall time: 4.9 seconds

---

Last 10 new commits:

​e30eeb4	Compute edge points directly
​ee46c5d	Optimize point creation using zero_vector
​f42ea43	Transpose databases of small reflexive polytopes
​0d734aa	Merge transpose_lattice_polytopes into points_without_PALP
​e56ec40	Read/write PointCollection in PALP format directly
​1b0467f	Do not precompute all polars of reflexive polytopes
​7cbc9de	Use zero_vector when constructing normals for speed
​fde45ab	Fix normal form computation
​7102115	Optimize getting points from PALP without embedding
​356cca4	Fix doctest depending on normals order and drop long time

[Matthias Köppe]

For large enough lattice point listing problems, considering using normaliz. Recent tickets have added it as a backend for Polyhedron.

[Andrey Novoseltsev]

Sure, but in toric geometry especially related to mirror symmetry it is important to handle very simple cases fast (e.g. all reflexive polytopes have precisely one interior lattice point, the origin).

Of course, it would be nice to handle all cases efficiently, and to be able to choose new methods appropriately I tried to first optimize existing one as much as possible, in particular not to take too much time just to convert data between different types.

Thanks for taking a look!

[git]

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

​d1eec1b

Fix print syntax

[git]

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

​dd6a616

Merge tag '8.1.rc0' into t/22524/points_without_PALP

[Travis Scrimshaw]

I think an even faster way than

n = N.zero_vector()
for i, coef in enumerate(c.coefficients()):
    n[i] = coef
n.set_immutable()
normals.append(n)

would be

n = N.element_class(N, c.coefficients(), coerce=False, copy=False)
n.set_immutable()

(although perhaps the coerce and copy are unnecessary). In general, if you know your input is good, you can bypass the extra overhead of __call__ and _element_constructor_ (there are a few places where you can do this I believe).

Also, it might be better to do

-                points = list(p.vertices())
-                for j in range(nv, m.ncols()):
-                    current = M.zero_vector()
-                    for i in range(M.rank()):
-                        current[i] = m[i, j]
-                    current.set_immutable()
-                    points.append(current)
+                points = list(p.vertices()) + m.columns(copy=False)
                 p._points = PointCollection(points, M)

Although I have not checked if the columns are longer than M.rank().

Otherwise LGTM.

[git]

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

​15e7043

Use element_class instead of zero_vector

[Andrey Novoseltsev]

Thank you for the pointer, element_class does lead to better timing and fewer function calls.

Regarding using matrix columns - that code deals with only some part of a matrix and, more importantly, is only used in "bad" cases: when dealing with non-full-dimensional polytopes or with a separate PALP call for a single polytope rather than sequence. My quick attempts to improve them now just broke things and I don't think they are worthy of more effort - the real solution is to rely on different backend. On the other hand, for computing points of many full-dimensional polytopes at once PALP interface still may be the fastest.

[Travis Scrimshaw]

I see. Well, you can do the same element_class trick for

                        current = M.zero_vector()
                        for i in range(M.rank()):
                            current[i] = m[i, j]

(there are 2 of these) and

            p = lattice.zero_vector()
            for j, e in enumerate(f.readline().split()):
                p[j] = int(e)

[git]

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

​e37e8f7

Use element_class instead of zero_vector

[Andrey Novoseltsev]

Travis - many many many thanks for reviewing all these tickets, making suggestions, and pushing me to implement them! I got 3-fold reduction in function calls for lattice_polytope.all_points(ReflexivePolytopes(3)) and the actual call to PALP is now responsible for a half of the time.

Regarding using columns of matrices - it is kind of incompatible with element_constructor since it expects lists or tuples, not vectors. Constructing a vector and then converting it to a list takes 3 times longer than just constructing this list directly from matrix elements.

I have also dropped copy=False, coerce=False since they do not seem to be used in the code and (while I cannot use vectors...) I can feed strings directly into element_constructor without converting them to int first.

[Travis Scrimshaw]

Replying to novoselt:

Travis - many many many thanks for reviewing all these tickets, making suggestions, and pushing me to implement them! I got 3-fold reduction in function calls for lattice_polytope.all_points(ReflexivePolytopes(3)) and the actual call to PALP is now responsible for a half of the time.

Not a problem. I'm sorry to took me a while to get through them all.

Regarding using columns of matrices - it is kind of incompatible with element_constructor since it expects lists or tuples, not vectors. Constructing a vector and then converting it to a list takes 3 times longer than just constructing this list directly from matrix elements.

I see.

I have also dropped copy=False, coerce=False since they do not seem to be used in the code and (while I cannot use vectors...) I can feed strings directly into element_constructor without converting them to int first.

Yea, that's not too surprising. I didn't think they were being used considering the datatypes, but it was more just-in-case.

LGTM. Positive review.