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:  sage7.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: 
Description (last modified by )
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)
Change History (38)
comment:1 Changed 8 years ago by
Description:  modified (diff) 

comment:2 Changed 6 years ago by
Cc:  Travis Scrimshaw added 

comment:3 Changed 6 years ago by
Cc:  Andrey Novoseltsev added 

comment:4 Changed 6 years ago by
Authors:  → Travis Scrimshaw 

Branch:  → public/geometry/speedup_integral_points18029 
Commit:  → 1e7e811a0ff797b2fb7665dc3a77a50a63a68583 
Milestone:  sage6.6 → sage7.3 
comment:5 Changed 6 years ago by
Status:  new → needs_review 

comment:6 Changed 6 years ago by
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
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
Commit:  1e7e811a0ff797b2fb7665dc3a77a50a63a68583 → 8cd48fcb69a150e4901c257a29dd8764d25474ac 

comment:10 followup: 12 Changed 6 years ago by
Reviewers:  → Matthias Koeppe 

Status:  needs_review → positive_review 
Looks good.
So... how about that normaliz interface?
comment:11 Changed 6 years ago by
Commit:  8cd48fcb69a150e4901c257a29dd8764d25474ac → 801226044e890ac38db3bbe6367cd4692b3095d8 

Status:  positive_review → needs_review 
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.

comment:12 Changed 6 years ago by
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:14 Changed 6 years ago by
Status:  needs_review → positive_review 

comment:16 Changed 6 years ago by
Status:  positive_review → needs_work 

comment:17 Changed 6 years ago by
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 microoptimizations 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
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 caseoncase basis, it's not a general rule that you should declare as much types as possible.
comment:19 Changed 6 years ago by
Commit:  801226044e890ac38db3bbe6367cd4692b3095d8 → d8292231fcf4a43f9f8a9e0d19ed4e15322da3e9 

comment:20 Changed 6 years ago by
Commit:  d8292231fcf4a43f9f8a9e0d19ed4e15322da3e9 → 9dd64fb8681bb4043caf4f108e5822559d684e83 

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.

comment:21 Changed 6 years ago by
Status:  needs_work → needs_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:23 Changed 6 years ago by
Reviewers:  Matthias Koeppe → Matthias Koeppe, Jeroen Demeyer 

comment:25 Changed 6 years ago by
Milestone:  sage7.3 → sage7.4 

Jeroen, if you have a minute, could you check my most recent changes? Thanks.
comment:27 Changed 6 years ago by
Cc:  Winfried Bruns added 

comment:28 Changed 6 years ago by
#20885 now has an implementation of integer_points
using normaliz.
comment:30 Changed 6 years ago by
Commit:  9dd64fb8681bb4043caf4f108e5822559d684e83 → 83d3322c4ba5f9a99be21e55e65489ec94862ef4 

Branch pushed to git repo; I updated commit sha1. New commits:
83d3322  Merge branch 'public/geometry/speedup_integral_points18029' of trac.sagemath.org:sage into public/geometry/speedup_integral_points18029

comment:31 Changed 6 years ago by
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
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/sagerebasing/yet/another/local/lib/python2.7/sitepackages/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/sagerebasing/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/sagerebasing/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/sagerebasing/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/sagerebasing/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_points18029' of trac.sagemath.org:sage into public/geometry/speedup_integral_points18029

comment:34 Changed 6 years ago by
I think that should be a separate ticket. (It likely might just be changing some int
to long
...)
comment:36 Changed 6 years ago by
Milestone:  sage7.4 → sage7.5 

Status:  needs_review → positive_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
Branch:  public/geometry/speedup_integral_points18029 → 83d3322c4ba5f9a99be21e55e65489ec94862ef4 

Resolution:  → fixed 
Status:  positive_review → closed 
I was able to get ~15% speedup:
In the current version, this takes ~14s. There is also could be some improvement for the simplex approach:
vs
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:
Optimizations and improve cython code for integral_points.pyx.