Optimize computing points of lattice polytopes

[novoselt]

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 Always use PPL for facet normals of lattice polytopes #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.

Depends on #22391
Depends on #22400

Component: geometry

Keywords: sd91

Author: Andrey Novoseltsev

Branch/Commit: e37e8f7

Reviewer: Travis Scrimshaw

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

[novoselt]

Branch: u/novoselt/points_without_PALP

[novoselt]

comment:2

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

[novoselt]

Commit: 356cca4

[mkoeppe]

comment:3

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

[novoselt]

comment:4

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!

[sagetrac-git]

Changed commit from 356cca4 to d1eec1b

[sagetrac-git]

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

d1eec1b

Fix print syntax

[simonbrandhorst]

Changed keywords from none to sd91

[sagetrac-git]

Changed commit from d1eec1b to dd6a616

[sagetrac-git]

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

dd6a616

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

[tscrim]

Reviewer: Travis Scrimshaw

[tscrim]

comment:9

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.

[sagetrac-git]

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

15e7043

Use element_class instead of zero_vector

[sagetrac-git]

Changed commit from dd6a616 to 15e7043

[novoselt]

comment:11

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.

[tscrim]

comment:12

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)

[sagetrac-git]

Changed commit from 15e7043 to e37e8f7

[sagetrac-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

[novoselt]

comment:14

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.

[tscrim]

comment:15

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.

[vbraun]

Changed branch from u/novoselt/points_without_PALP to e37e8f7