speed up integral point enumeration

[videlec]

Some of the cython code in sage/geometry/integral_points.pyx can be optimized a lot (properly defining the types of objects, use the set_unsafe feature from #17562) and more.

Depends on #17562

CC: @novoselt @tscrim @w-bruns

Component: geometry

Author: Travis Scrimshaw

Branch/Commit: 83d3322

Reviewer: Matthias Koeppe, Jeroen Demeyer

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

[tscrim]

Branch: public/geometry/speedup_integral_points-18029

[tscrim]

comment:4

I was able to get ~15% speedup:

sage: sage: ieqs = [(-1, -1, -1, -1, -1, -1, -1, -1, -1),
....:         (0, -1, 0, 0, 0, 0, 0, 0, 0),
....:         (0, -1, 0, 2, -1, 0, 0, 0, 0),
....:         (0, 0, -1, -1, 2, -1, 0, 0, 0),
....:         (0, 2, 0, -1, 0, 0, 0, 0, 0),
....:         (0, 0, 0, 0, 0, 0, 0, -1, 2),
....:         (1, 0, 2, 0, -1, 0, 0, 0, 0),
....:         (0, 0, 0, 0, -1, 2, -1, 0, 0),
....:         (0, 0, 0, 0, 0, 0, 0, 0, -1),
....:         (0, 0, 0, 0, 0, -1, 2, -1, 0),
....:         (0, 0, 0, 0, 0, 0, -1, 2, -1)]
sage: P = Polyhedron(ieqs=ieqs)
sage: %time len(P.integral_points())
CPU times: user 12.4 s, sys: 5 µs, total: 12.4 s
Wall time: 12.4 s
4

In the current version, this takes ~14s. There is also some improvement for the simplex approach:

sage: from sage.geometry.integral_points import simplex_points
sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)] 
sage: simplex = Polyhedron(v); simplexA 4-dimensional polyhedron in ZZ^4 defined as the convex hull of 5 vertices
sage: %timeit len(simplex_points(simplex.vertices()))
100 loops, best of 3: 15.1 ms per loop

vs

10 loops, best of 3: 24.1 ms per loop

I used all the tricks that I know to squeeze speed out with doing any rewrites. However, it would be better to implement/use a priority queue (or linked list) data structure instead of a straight list for the inequalities since we are often moving one to the front as part of the algorithm.
Yet I think this is a good gain and should be included in the present state (unless someone feels like doing something more invasive, which I would then review).

---

New commits:

1e7e811

Optimizations and improve cython code for integral_points.pyx.

[tscrim]

Author: Travis Scrimshaw

[tscrim]

Commit: 1e7e811

[mkoeppe]

normaliz input

[mkoeppe]

comment:6

Attachment: 18029-example.in.gz

This looks fine, but let me point out that the performance of this code is very far from the state of the art.
The first example in dimension 8 is solved by normaliz in 0.1s (I have attached an input file for normaliz). So a simple interface to normaliz (#20885) would bring huge speedups.

[mkoeppe]

comment:7

That said, of course it's nice to have a basic, generic implementation in Sage if one does not want to install normaliz; but it turns out that the implementation is actually not very generic (#21037).

[sagetrac-git]

Changed commit from 1e7e811 to 8cd48fc

[sagetrac-git]

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

25903fc	Merge branch 'public/geometry/speedup_integral_points-18029' of trac.sagemath.org:sage into public/geometry/speedup_integral_points-18029
8cd48fc	Trying to get a little bit more speed out.

[tscrim]

comment:9

FYI - this merged cleanly with #21037.

[mkoeppe]

Reviewer: Matthias Koeppe

[mkoeppe]

comment:10

Looks good.

So... how about that normaliz interface?

[sagetrac-git]

Changed commit from 8cd48fc to 8012260

[sagetrac-git]

Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. New commits:

8012260

Don't manipulate the lengths of lists unless necessary.

[tscrim]

comment:12

Replying to @mkoeppe:

Looks good.

Sorry, I found one other place where I got another .4s in the above example.

So... how about that normaliz interface?

That's quite a difference, and it would be great to have it available (and be extremely useful to me to have integral points that fast). I don't know how much I could help it writing the interface; I know somewhat how to write the cython bindings, but not what to use/expose from noraliz. I will definitely be happy to review it.

[tscrim]

comment:13

Just in case you didn't see it, could you check my last push?

[sagetrac-git]

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

76cf443	Merge branch 'develop' into public/geometry/speedup_integral_points-18029
d829223	Using Cython directives and forcing min/max_box to be lists.

[sagetrac-git]

Changed commit from d829223 to 9dd64fb

[sagetrac-git]

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

9dd64fb

Using Cython directives and forcing box_min/max to be lists.

[tscrim]

comment:21

I've added the Cython directives (I didn't know about them, thank you), and I've added a few more type declarations. I believe I have made it so once the type is known, it is kept that way so Cython should never have to recheck the type of things.

[tscrim]

comment:22

I've also added the except -1 to the necessary c(p)def methods.

[mkoeppe]

Changed reviewer from Matthias Koeppe to Matthias Koeppe, Jeroen Demeyer

[mkoeppe]

comment:24

It builds and works OK.
Jeroen, could you take another look?

[tscrim]

comment:25

Jeroen, if you have a minute, could you check my most recent changes? Thanks.

[tscrim]

comment:26

Ping?

[mkoeppe]

comment:28

#20885 now has an implementation of integer_points using normaliz.

[mkoeppe]

comment:29

This ticket would need rebasing also.

[sagetrac-git]

Changed commit from 9dd64fb to 83d3322

[sagetrac-git]

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

83d3322

Merge branch 'public/geometry/speedup_integral_points-18029' of trac.sagemath.org:sage into public/geometry/speedup_integral_points-18029

[tscrim]

comment:31

Rebased. It would be nice to finalize this and get this in until #20885 gets merged and becomes the default.

[mkoeppe]

comment:32

By the way, I noticed the following bug in Sage (without this branch -- haven't tested yet with this branch!):

sage: %timeit (Polyhedron(vertices=((0, 0), (1789345,37121))) + 1/1000*polytopes.hypercube(2)).integral_points()
...
/Users/mkoeppe/s/sage/sage-rebasing/yet/another/local/lib/python2.7/site-packages/sage/geometry/polyhedron/base.pyc in integral_points(self, threshold)
   4357                 box_points<threshold:
   4358             from sage.geometry.integral_points import rectangular_box_points
-> 4359             return rectangular_box_points(box_min, box_max, self)
   4360 
   4361         # for more complicate polytopes, triangulate & use smith normal form

/Users/mkoeppe/s/sage/sage-rebasing/src/sage/geometry/integral_points.pyx in sage.geometry.integral_points.rectangular_box_points (build/cythonized/sage/geometry/integral_points.c:6809)()
    552     v = vector(ZZ, d)
    553     if not return_saturated:
--> 554         for p in loop_over_rectangular_box_points(box_min, box_max, inequalities, d, count_only):
    555             #  v = vector(ZZ, orig_perm.action(p))   # too slow
    556             for i in range(0,d):

/Users/mkoeppe/s/sage/sage-rebasing/src/sage/geometry/integral_points.pyx in sage.geometry.integral_points.loop_over_rectangular_box_points (build/cythonized/sage/geometry/integral_points.c:7772)()
    601     while True:
    602         sig_check()
--> 603         inequalities.prepare_inner_loop(p)
    604         i_min = box_min[0]
    605         i_max = box_max[0]

/Users/mkoeppe/s/sage/sage-rebasing/src/sage/geometry/integral_points.pyx in sage.geometry.integral_points.InequalityCollection.prepare_inner_loop (build/cythonized/sage/geometry/integral_points.c:14134)()
   1256         """
   1257         for ineq in self.ineqs_int:
-> 1258             (<Inequality_int>ineq).prepare_inner_loop(p)
   1259         for ineq in self.ineqs_generic:
   1260             (<Inequality_generic>ineq).prepare_inner_loop(p)

/Users/mkoeppe/s/sage/sage-rebasing/src/sage/geometry/integral_points.pyx in sage.geometry.integral_points.Inequality_int.prepare_inner_loop (build/cythonized/sage/geometry/integral_points.c:10574)()
    938         cdef int j
    939         if self.dim>1:
--> 940             self.cache = self.cache_next + self.coeff_next*p[1]
    941         else:
    942             self.cache = self.cache_next

OverflowError: value too large to convert to int

(I can make this a separate ticket if that's better.)

---

New commits:

83d3322

Merge branch 'public/geometry/speedup_integral_points-18029' of trac.sagemath.org:sage into public/geometry/speedup_integral_points-18029

[mkoeppe]

comment:33

(same error also with this branch).

[tscrim]

comment:34

I think that should be a separate ticket. (It likely might just be changing some int to long...)

[mkoeppe]

comment:35

That's now #21993.

[mkoeppe]

comment:36

Looks good to me.
I'm setting this to positive review -- but Jeroen may have more comments.

[vbraun]

Changed branch from public/geometry/speedup_integral_points-18029 to 83d3322