# Speed of enumerating integer points in polytope

I'm trying to enumerate all integer points in some polytopes, but sage takes way too long.

This is a pretty minimal example of one of these polytopes.

q = MixedIntegerLinearProgram(maximization=False)
x = q.new_variable(integer=True, nonnegative=True)

q.add_constraint(88 * x[2] == -306 * x[0] + 1 * x[1])
q.add_constraint( 1 * x[3] == 42636 * x[0] + 42 * x[1])
q.add_constraint(11 * x[4] == 30804 * x[0] + -42 * x[1])

P = q.polyhedron()

P
# A 1-dimensional polyhedron in QQ^5 defined as the convex hull of 2 vertices

P.Vrepresentation()
# (A vertex at (1, 34/7, 5134/7, 73440, 0), A vertex at (1, 0, 306, 55488, 1632))


Now, I would like to enumerate the integral points in this (1-dimensional) polytope. Looking at the second coordinate of the two boundary points, it is clear there can be at most five such points (in fact, there are exactly five).

However, P.integral_points() seems to take forever.

Is there anything I can do to speed this up? I would like to do similar problems that are higher-dimensional, and have a lot more constraints (bounding hyperplanes).

Earlier, I used a MILP library (COIN) to solve the system, giving one integer point, then further constrain the system to exclude this solution, solve again, etc. This is fast enough for this particular example, but becomes unfeasible when the dimension of the polytope and/or the number of integer points increases.

edit retag close merge delete

Maybe you can find one point, shift it to zero and then work on the subspace spanned by the polytope ? That would make the ambient dimension smaller.

( 2014-10-15 13:48:42 -0500 )edit

Sort by » oldest newest most voted

In Sage 8.1, I can reproduce your code as follows (notice that I added "ppl" as a solver since the default behavior changed between your post and this answer).

sage: q = MixedIntegerLinearProgram(maximization=False, solver="ppl")
....: x = q.new_variable(integer=True, nonnegative=True)
....:
....: q.add_constraint(88 * x[2] == -306 * x[0] + 1 * x[1])
....: q.add_constraint( 1 * x[3] == 42636 * x[0] + 42 * x[1])
....: q.add_constraint(11 * x[4] == 30804 * x[0] + -42 * x[1])


In order to speed things up, I use normaliz.

(You can install it by typing the commands in a terminal:

sage -i normaliz
sage -i pynormaliz


)

and then use the backend normaliz for the polyhedron:

....: P = q.polyhedron(backend="normaliz")
....: P.integral_points()
....:
((1, 0, 306, 55488, 1632),
(1, 1, 394, 59184, 1296),
(1, 2, 482, 62880, 960),
(1, 3, 570, 66576, 624),
(1, 4, 658, 70272, 288))


which should appear instantaneously.

more

try using Latte Integrale (an experimental package in Sage), called latte_int; it implements Barvinok's algorithm for integer point counting. There is however no interface between Latte and Sage developed, so you need to create files to be real into Latte...

more

Unfortunately, it seems Latte Integrale does not enumerate the lattice points, it only counts them. I need the coordinates of all the integer points, for further calculations.

( 2014-10-13 06:27:45 -0500 )edit

it does count them, in a way; you have to call count with the option --multivariate-generating-function; such a function "simplifies" to a (Laurent) polynomial, with exponent vectors of the monomials corresponding to the integer points.

( 2014-10-13 10:07:10 -0500 )edit

I see. It does not sound like a very convenient way to represent the thousands of points in somewhat higher-dimensional spaces of my other cases. I'll give it a try, though. Thanks for the suggestion!

( 2014-10-13 11:29:33 -0500 )edit

it is the most compact way known, and this is what counts.IIRC, the rational function consists of summands, each of them a rational function corresponding to a simplex, simplices forming a kind of partition of the polytope.So one can probably only take summands corresponding to x^K with positive signs, if you don't mind getting the same point several times.

( 2014-10-14 02:16:57 -0500 )edit

I'm trying to do an intermediately big case, where the dimension is 8, and the polytope has 1804 vertices. Latte Integrale has already been computing for over 20 minutes (just counting the points, for now), and there is no indication how long it will take. It just says: "431000 unimodular cones done", this number increasing by a 1000 every few seconds, but I don't know how many unimodular cones it must do.

( 2014-10-14 06:40:05 -0500 )edit