Count integral points without PALP

[vbraun]

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​

[vbraun]

Updated patch

[vbraun]

The order of points is different from PALP (unsurprisingly), so doctests needed to be fixed. This is done in the last patch.

[novoselt]

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

[vbraun]

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

[novoselt]

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

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

[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 ;-)

[vbraun]

3. I added a len(pts) to the len(set(pts)) to clarify - the doctests checks that all points are distinct.
4. 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.
5. 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.
6. 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).

[vbraun]

Updated patch fixes the remaining points of comment:11

[vbraun]

Rebased patch

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

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

[novoselt]

OK, positive review!

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

[jdemeyer]

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.

[jdemeyer]

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

[vbraun]

Updated patch

[vbraun]

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

[jdemeyer]

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.