Opened 10 years ago
Closed 10 years ago
#11429 closed enhancement (fixed)
Count integral points without PALP
Reported by: | vbraun | Owned by: | AlexGhitza |
---|---|---|---|
Priority: | major | Milestone: | sage-5.0 |
Component: | geometry | Keywords: | sd31 |
Cc: | novoselt, zaf | Merged in: | sage-5.0.beta2 |
Authors: | Volker Braun | Reviewers: | Andrey Novoseltsev, Jeroen Demeyer |
Report Upstream: | N/A | Work issues: | |
Branch: | Commit: | ||
Dependencies: | #11312, #9958 | Stopgaps: |
Description (last modified by )
We want our own code to enumerate lattice points in polyhedra because:
- Going through the PALP pexpect interface is annoyingly slow.
- no more compile-time bounds
- It seems like PALP uses a very unsophisticated algorithm:
sage: v = [(1,0,7,-1), (-2,-2,4,-3), (-1,-1,-1,4), (2,9,0,-5), (-2,-1,5,1)] sage: lp = LatticePolytope(matrix(v).transpose()); lp A lattice polytope: 4-dimensional, 5 vertices. sage: lp.npoints()
takes forever with PALP but only 500ms with my Python code.
Comparing timings, it seems that PALP always runs over the integral points of a rectangular bounding box. This is good for small polytopes (low overhead) but bad for large ones. To match PALP's speed for small polytopes, I implemented the same algorithm in Cython (the second patch) and use it for bounding boxes containing <50k points.
Apply:
Attachments (3)
Change History (25)
comment:1 Changed 10 years ago by
- Dependencies set to #11312
Changed 10 years ago by
comment:2 Changed 10 years ago by
- Cc novoselt added
- Description modified (diff)
- Status changed from new to needs_review
comment:3 Changed 10 years ago by
The order of points is different from PALP (unsurprisingly), so doctests needed to be fixed. This is done in the last patch.
comment:4 Changed 10 years ago by
- Description modified (diff)
comment:5 Changed 10 years ago by
- Description modified (diff)
- Keywords sd31 added
- Reviewers set to Andrey Novoseltsev
comment:6 Changed 10 years ago by
Regarding the reuse of objects, there was this ticket which may be related/useful: #8702.
comment:7 Changed 10 years ago by
- Cc zaf added
comment:8 Changed 10 years ago by
Updated patch because of Andrey's documentation changes to #11312.
comment:9 Changed 10 years ago by
- Component changed from algebraic geometry to geometry
Hi Volker,
Would it make any sense to somehow unite these functions/classes with your PPL wrapper, which already has stuff for inequalities and their systems? Do you know if they (PPL-people) have plans on algorithms for point counts and point listing?
There are also functions without documentation or without full description of input/output...
comment:10 Changed 10 years ago by
I don't know the PPL future plans, but I don't think that some Cython code would fit into their project. Its also not functionality that fits with their core misson objective, I think.
I haven't used the usual standards of documentation for cdef'd functions since they aren't reachable from Python, so you'll never see the help. Nor can you write Python doctests for them. Also, I think documenting input/output is implicit when you have typed arguments.
comment:11 Changed 10 years ago by
Line numbers are for the new module in the big patch:
- 224: should be
if A is not None:
for robustness (zero matrices evaluate toFalse
) - 231: should be
if A is None:
- 273-275: Why there is a conversion to set? Is it expected that the output contains repetitions?
- 293: shouldn't
translate_points
be used here? - 336: No description of input/output and what happens if the polyhedron is not specified.
- In the code of this function: why is it necessary to use this permutation? It seems like an extra operation on all of the points, although I believe that there is some reason ;-) Perhaps it can be explained in the documentation or comments.
- 421: Incomplete INPUT, missing OUTPUT, no real EXAMPLE.
- 440: While this function is inaccessible from Sage directly, it still would be nice to have a comment on what it does and what is d.
- 477: If I am correct, this condition is always true. In which case it can be removed, as well the next line, and 444 moved into the beginning of outer while. (Just for making code a little simpler.)
- I didn't yet looked carefully at the rest, but just one more pick at 575: Is it really necessary to have
DEF INEQ_INT_MAX_DIM = 20
? I mean can't it be without hard-coded limits?
Will look over the second half of the new module in the next few days, I am fine with the rest of the code.
Speedwise PALP is still more efficient for small polytopes, but of course only if one can call PALP at once for a lot of them, which is not always possible:
sage: len(ReflexivePolytopes(3)) 4319 sage: %time sage: ps = [Polyhedron(vertices=p.vertices().columns()) for p in ReflexivePolytopes(3)] CPU time: 83.87 s, Wall time: 218.85 s sage: %time sage: for p, lp in zip(ps, ReflexivePolytopes(3)): ... if len(p.integral_points()) != lp.npoints(): ... print lp CPU time: 10.10 s, Wall time: 10.30 s sage: %time sage: lattice_polytope.all_points(ReflexivePolytopes(3)) CPU time: 2.90 s, Wall time: 3.46 s
The last line is actually somewhat sad: PALP spends half a second to compute all the points and then Sage spends 3 seconds to read its output and convert into matrices. But hopefully objects of new generations will be more efficient ;-)
comment:12 Changed 10 years ago by
- I added a
len(pts)
to thelen(set(pts))
to clarify - the doctests checks that all points are distinct. - The coordinates are permuted such that the largest edge of the bounding box becomes the inner loop. The inner loop is optimized to run very fast, so we save wall time by doing the maximal amount of work there.
- No
inc
will usually be1
. Usually you'll only have to callprepare_inner_loop
between two inner loops. Only if the "odometer ticks over" you need to callprepare_next_to_inner_loop
. In terms of the coordinates(x1,x2,x3,...)
of the box, the inner loop runs overx1
. It suffices to only callprepare_inner_loop
as long as onlyx1
,x2
change. You need to callprepare_next_to_inner_loop
if more coordinates changed from one inner loop to the next. - The dimension limit is only used for the machine integer representation of inequalities. Its not a real limit, in higher dimensions the generic Python implementation is used.
There is some overhead to avoid integer overflows and the like, which is a good thing. Possible improvements are
- don't permute coordinates (overhead) for really really small polyhedra, and
- allow to just count points without enumeration (since we already find upper/lower bound in inner loop).
comment:13 Changed 10 years ago by
Updated patch fixes the remaining points of comment:11
comment:14 Changed 10 years ago by
Last small picks, lines again refer to the new module in the big patch:
- 483: Looks like something was not pasted.
- 720: If it is indeed necessary to check if the dimension is less than 1, perhaps a different error message is in order?
- 912: A perhaps naive question: wouldn't it be more efficient to work in an affine span of the polytope if there are equations? (If yes, I don't insist on doing it, but a "todo" note could be nice.)
comment:15 Changed 10 years ago by
- The
permuted_list()
function was just a leftover that I forgot to delete. - The
OverflowError
is caught internally and the generic version (as opposed to the optimized machine int version) is used. So the user never gets to see the error message. - There is definitely room for improvement for special lattice polytopes. But for really small lattice polytopes I think the naive thing is faster than doing linear algebra for the affine span.
comment:16 Changed 10 years ago by
- Status changed from needs_review to positive_review
OK, positive review!
With deep apollogies for the duration of "the next few days"...
comment:17 Changed 10 years ago by
- Milestone changed from sage-4.8 to sage-5.0
comment:18 Changed 10 years ago by
- Dependencies changed from #11312 to #11312, #9958
- Status changed from positive_review to needs_work
When running this with Python-2.7.2 (#9958), there is a doctest failure:
sage -t -force_lib devel/sage/sage/geometry/integral_points.pyx ********************************************************************** File "/mnt/usb1/scratch/jdemeyer/merger/sage-5.0.beta1/devel/sage-main/sage/geometry/integral_points.pyx", line 668: sage: Inequality_int([2,3,7], -5*10^50, [10]*3) Expected: Traceback (most recent call last): ... OverflowError: long int too large to convert to int Got: Traceback (most recent call last): File "/mnt/usb1/scratch/jdemeyer/merger/sage-5.0.beta1/local/bin/ncadoctest.py", line 1231, in run_one_test self.run_one_example(test, example, filename, compileflags) File "/mnt/usb1/scratch/jdemeyer/merger/sage-5.0.beta1/local/bin/sagedoctest.py", line 38, in run_one_example OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags) File "/mnt/usb1/scratch/jdemeyer/merger/sage-5.0.beta1/local/bin/ncadoctest.py", line 1172, in run_one_example compileflags, 1) in test.globs File "<doctest __main__.example_14[7]>", line 1, in <module> Inequality_int([Integer(2),Integer(3),Integer(7)], -Integer(5)*Integer(10)**Integer(50), [Integer(10)]*Integer(3))###line 668: sage: Inequality_int([2,3,7], -5*10^50, [10]*3) File "integral_points.pyx", line 705, in sage.geometry.integral_points.Inequality_int.__cinit__ (sage/geometry/integral_points.c:495 0) self.b = self._to_int(b) File "integral_points.pyx", line 685, in sage.geometry.integral_points.Inequality_int._to_int (sage/geometry/integral_points.c:4765) return int(x) # raises OverflowError in Cython if necessary OverflowError: Python int too large to convert to C long **********************************************************************
Since Python-2.7 will surely be merged in sage-5.0.beta0, you should fix the error message.
comment:19 Changed 10 years ago by
By the way: #9958 suggests the error message might be different on 32-bit and 64-bit systems.
comment:20 follow-up: ↓ 21 Changed 10 years ago by
- Status changed from needs_work to needs_review
I've modified the doctest to only expect OverflowError: ...
instead of a particular message.
comment:21 in reply to: ↑ 20 Changed 10 years ago by
- Reviewers changed from Andrey Novoseltsev to Andrey Novoseltsev, Jeroen Demeyer
- Status changed from needs_review to positive_review
Replying to vbraun:
I've modified the doctest to only expect
OverflowError: ...
instead of a particular message.
I haven't yet tested it, but this should fix the issue.
comment:22 Changed 10 years ago by
- Merged in set to sage-5.0.beta2
- Resolution set to fixed
- Status changed from positive_review to closed
Updated patch