Opened 8 years ago

Closed 6 years ago

#18029 closed enhancement (fixed)

speed up integral point enumeration

Reported by: Vincent Delecroix Owned by:
Priority: major Milestone: sage-7.5
Component: geometry Keywords:
Cc: Andrey Novoseltsev, Travis Scrimshaw, Winfried Bruns Merged in:
Authors: Travis Scrimshaw Reviewers: Matthias Koeppe, Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: 83d3322 (Commits, GitHub, GitLab) Commit: 83d3322c4ba5f9a99be21e55e65489ec94862ef4
Dependencies: #17562 Stopgaps:

Status badges

Description (last modified by Vincent Delecroix)

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.

Attachments (1)

18029-example.in (302 bytes) - added by Matthias Köppe 6 years ago.
normaliz input

Download all attachments as: .zip

Change History (38)

comment:1 Changed 8 years ago by Vincent Delecroix

Description: modified (diff)

comment:2 Changed 6 years ago by Travis Scrimshaw

Cc: Travis Scrimshaw added

comment:3 Changed 6 years ago by Andrey Novoseltsev

Cc: Andrey Novoseltsev added

comment:4 Changed 6 years ago by Travis Scrimshaw

Authors: Travis Scrimshaw
Branch: public/geometry/speedup_integral_points-18029
Commit: 1e7e811a0ff797b2fb7665dc3a77a50a63a68583
Milestone: sage-6.6sage-7.3

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:

1e7e811Optimizations and improve cython code for integral_points.pyx.
Last edited 6 years ago by Travis Scrimshaw (previous) (diff)

comment:5 Changed 6 years ago by Travis Scrimshaw

Status: newneeds_review

Changed 6 years ago by Matthias Köppe

Attachment: 18029-example.in added

normaliz input

comment:6 Changed 6 years ago by Matthias Köppe

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.

comment:7 Changed 6 years ago by Matthias Köppe

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).

comment:8 Changed 6 years ago by git

Commit: 1e7e811a0ff797b2fb7665dc3a77a50a63a685838cd48fcb69a150e4901c257a29dd8764d25474ac

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

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

comment:9 Changed 6 years ago by Travis Scrimshaw

FYI - this merged cleanly with #21037.

comment:10 Changed 6 years ago by Matthias Köppe

Reviewers: Matthias Koeppe
Status: needs_reviewpositive_review

Looks good.

So... how about that normaliz interface?

comment:11 Changed 6 years ago by git

Commit: 8cd48fcb69a150e4901c257a29dd8764d25474ac801226044e890ac38db3bbe6367cd4692b3095d8
Status: positive_reviewneeds_review

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

8012260Don't manipulate the lengths of lists unless necessary.

comment:12 in reply to:  10 Changed 6 years ago by Travis Scrimshaw

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.

comment:13 Changed 6 years ago by Travis Scrimshaw

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

comment:14 Changed 6 years ago by Matthias Köppe

Status: needs_reviewpositive_review

comment:15 Changed 6 years ago by Travis Scrimshaw

Thank you.

comment:16 Changed 6 years ago by Volker Braun

Status: positive_reviewneeds_work

comment:17 Changed 6 years ago by Jeroen Demeyer

You need to use an except value when you have a c(p)def function with a return type: replace cpdef bint ...(...) by cpdef bint ...(...) except -1.

And since you seem to like Cython micro-optimizations a lot in this ticket, you missed an obvious one: setting the boundscheck and wraparound directives to False to speed up list indexing.

comment:18 Changed 6 years ago by Jeroen Demeyer

I should also add that declaring the type of a variable or return value may or may not speed up the code. You gain if Cython can use special optimizations for that type. But you lose if Cython needs to add pointless checks that the thing you are assigning to the variable has the correct type.

Deciding whether to declare the type of a variable should be done on a case-on-case basis, it's not a general rule that you should declare as much types as possible.

comment:19 Changed 6 years ago by git

Commit: 801226044e890ac38db3bbe6367cd4692b3095d8d8292231fcf4a43f9f8a9e0d19ed4e15322da3e9

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

76cf443Merge branch 'develop' into public/geometry/speedup_integral_points-18029
d829223Using Cython directives and forcing min/max_box to be lists.

comment:20 Changed 6 years ago by git

Commit: d8292231fcf4a43f9f8a9e0d19ed4e15322da3e99dd64fb8681bb4043caf4f108e5822559d684e83

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

9dd64fbUsing Cython directives and forcing box_min/max to be lists.

comment:21 Changed 6 years ago by Travis Scrimshaw

Status: needs_workneeds_review

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.

comment:22 Changed 6 years ago by Travis Scrimshaw

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

comment:23 Changed 6 years ago by Matthias Köppe

Reviewers: Matthias KoeppeMatthias Koeppe, Jeroen Demeyer

comment:24 Changed 6 years ago by Matthias Köppe

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

comment:25 Changed 6 years ago by Travis Scrimshaw

Milestone: sage-7.3sage-7.4

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

comment:26 Changed 6 years ago by Travis Scrimshaw

Ping?

comment:27 Changed 6 years ago by Matthias Köppe

Cc: Winfried Bruns added

comment:28 Changed 6 years ago by Matthias Köppe

#20885 now has an implementation of integer_points using normaliz.

comment:29 Changed 6 years ago by Matthias Köppe

This ticket would need rebasing also.

comment:30 Changed 6 years ago by git

Commit: 9dd64fb8681bb4043caf4f108e5822559d684e8383d3322c4ba5f9a99be21e55e65489ec94862ef4

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

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

comment:31 Changed 6 years ago by Travis Scrimshaw

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

comment:32 Changed 6 years ago by Matthias Köppe

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:

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

comment:33 Changed 6 years ago by Matthias Köppe

(same error also with this branch).

comment:34 Changed 6 years ago by Travis Scrimshaw

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

comment:35 Changed 6 years ago by Matthias Köppe

That's now #21993.

comment:36 Changed 6 years ago by Matthias Köppe

Milestone: sage-7.4sage-7.5
Status: needs_reviewpositive_review

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

comment:37 Changed 6 years ago by Volker Braun

Branch: public/geometry/speedup_integral_points-1802983d3322c4ba5f9a99be21e55e65489ec94862ef4
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.