Opened 9 years ago

Closed 9 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 novoselt)

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:

  1. trac_11429_native_enumeration_of_lattice_polytope_points.patch
  2. trac_11429_cythonize_lattice_points.patch
  3. trac_11429_fix_doctests.patch

Attachments (3)

trac_11429_native_enumeration_of_lattice_polytope_points.patch (16.0 KB) - added by vbraun 9 years ago.
Updated patch
trac_11429_fix_doctests.patch (4.8 KB) - added by vbraun 9 years ago.
Rebased patch
trac_11429_cythonize_lattice_points.patch (52.0 KB) - added by vbraun 9 years ago.
Updated patch

Download all attachments as: .zip

Change History (25)

comment:1 Changed 9 years ago by vbraun

  • Dependencies set to #11312

Changed 9 years ago by vbraun

Updated patch

comment:2 Changed 9 years ago by vbraun

  • Cc novoselt added
  • Description modified (diff)
  • Status changed from new to needs_review

comment:3 Changed 9 years ago by vbraun

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 9 years ago by vbraun

  • Description modified (diff)

comment:5 Changed 9 years ago by novoselt

  • Description modified (diff)
  • Keywords sd31 added
  • Reviewers set to Andrey Novoseltsev

comment:6 Changed 9 years ago by novoselt

Regarding the reuse of objects, there was this ticket which may be related/useful: #8702.

comment:7 Changed 9 years ago by zaf

  • Cc zaf added

comment:8 Changed 9 years ago by vbraun

Updated patch because of Andrey's documentation changes to #11312.

comment:9 Changed 9 years ago by novoselt

  • 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 9 years ago by vbraun

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 9 years ago by novoselt

Line numbers are for the new module in the big patch:

  1. 224: should be if A is not None: for robustness (zero matrices evaluate to False)
  2. 231: should be if A is None:
  3. 273-275: Why there is a conversion to set? Is it expected that the output contains repetitions?
  4. 293: shouldn't translate_points be used here?
  5. 336: No description of input/output and what happens if the polyhedron is not specified.
  6. 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.
  7. 421: Incomplete INPUT, missing OUTPUT, no real EXAMPLE.
  8. 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.
  9. 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.)
  10. 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 9 years ago by vbraun

  1. I added a len(pts) to the len(set(pts)) to clarify - the doctests checks that all points are distinct.
  2. 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.
  3. No inc will usually be 1. Usually you'll only have to call prepare_inner_loop between two inner loops. Only if the "odometer ticks over" you need to call prepare_next_to_inner_loop. In terms of the coordinates (x1,x2,x3,...) of the box, the inner loop runs over x1. It suffices to only call prepare_inner_loop as long as only x1, x2 change. You need to call prepare_next_to_inner_loop if more coordinates changed from one inner loop to the next.
  4. 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 9 years ago by vbraun

Updated patch fixes the remaining points of comment:11

Changed 9 years ago by vbraun

Rebased patch

comment:14 Changed 9 years ago by novoselt

Last small picks, lines again refer to the new module in the big patch:

  1. 483: Looks like something was not pasted.
  2. 720: If it is indeed necessary to check if the dimension is less than 1, perhaps a different error message is in order?
  3. 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 9 years ago by vbraun

  1. The permuted_list() function was just a leftover that I forgot to delete.
  2. 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.
  3. 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 9 years ago by novoselt

  • 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 9 years ago by jdemeyer

  • Milestone changed from sage-4.8 to sage-5.0

comment:18 Changed 9 years ago by jdemeyer

  • 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 9 years ago by jdemeyer

By the way: #9958 suggests the error message might be different on 32-bit and 64-bit systems.

Changed 9 years ago by vbraun

Updated patch

comment:20 follow-up: Changed 9 years ago by vbraun

  • 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 9 years ago by jdemeyer

  • 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 9 years ago by jdemeyer

  • Merged in set to sage-5.0.beta2
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.