Opened 11 years ago

Closed 11 years ago

#11429 closed enhancement (fixed)

Count integral points without PALP

Reported by: Volker Braun Owned by: Alex Ghitza
Priority: major Milestone: sage-5.0
Component: geometry Keywords: sd31
Cc: Andrey Novoseltsev, Zafeirakis Zafeirakopoulos 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:

Status badges

Description (last modified by Andrey Novoseltsev)

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 Volker Braun 11 years ago.
Updated patch
trac_11429_fix_doctests.patch (4.8 KB) - added by Volker Braun 11 years ago.
Rebased patch
trac_11429_cythonize_lattice_points.patch (52.0 KB) - added by Volker Braun 11 years ago.
Updated patch

Download all attachments as: .zip

Change History (25)

comment:1 Changed 11 years ago by Volker Braun

Dependencies: #11312

Changed 11 years ago by Volker Braun

Updated patch

comment:2 Changed 11 years ago by Volker Braun

Cc: Andrey Novoseltsev added
Description: modified (diff)
Status: newneeds_review

comment:3 Changed 11 years ago by Volker Braun

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 11 years ago by Volker Braun

Description: modified (diff)

comment:5 Changed 11 years ago by Andrey Novoseltsev

Description: modified (diff)
Keywords: sd31 added
Reviewers: Andrey Novoseltsev

comment:6 Changed 11 years ago by Andrey Novoseltsev

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

comment:7 Changed 11 years ago by Zafeirakis Zafeirakopoulos

Cc: Zafeirakis Zafeirakopoulos added

comment:8 Changed 11 years ago by Volker Braun

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

comment:9 Changed 11 years ago by Andrey Novoseltsev

Component: algebraic geometrygeometry

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 11 years ago by Volker Braun

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 11 years ago by Andrey Novoseltsev

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 11 years ago by Volker Braun

  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 11 years ago by Volker Braun

Updated patch fixes the remaining points of comment:11

Changed 11 years ago by Volker Braun

Rebased patch

comment:14 Changed 11 years ago by Andrey Novoseltsev

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 11 years ago by Volker Braun

  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 11 years ago by Andrey Novoseltsev

Status: needs_reviewpositive_review

OK, positive review!

With deep apollogies for the duration of "the next few days"...

comment:17 Changed 11 years ago by Jeroen Demeyer

Milestone: sage-4.8sage-5.0

comment:18 Changed 11 years ago by Jeroen Demeyer

Dependencies: #11312#11312, #9958
Status: positive_reviewneeds_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 11 years ago by Jeroen Demeyer

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

Changed 11 years ago by Volker Braun

Updated patch

comment:20 Changed 11 years ago by Volker Braun

Status: needs_workneeds_review

I've modified the doctest to only expect OverflowError: ... instead of a particular message.

comment:21 in reply to:  20 Changed 11 years ago by Jeroen Demeyer

Reviewers: Andrey NovoseltsevAndrey Novoseltsev, Jeroen Demeyer
Status: needs_reviewpositive_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 11 years ago by Jeroen Demeyer

Merged in: sage-5.0.beta2
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.