Ask Your Question
3

Speed of enumerating integer points in polytope

asked 2014-10-10 18:53:00 +0100

Arie gravatar image

updated 2015-01-13 20:34:42 +0100

FrédéricC gravatar image

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(     x[0] == 1)
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 flag offensive close merge delete

Comments

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.

FrédéricC gravatar imageFrédéricC ( 2014-10-15 20:48:42 +0100 )edit

2 Answers

Sort by » oldest newest most voted
2

answered 2014-10-11 23:31:10 +0100

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

edit flag offensive delete link more

Comments

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.

Arie gravatar imageArie ( 2014-10-13 13:27:45 +0100 )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.

Dima gravatar imageDima ( 2014-10-13 17:07:10 +0100 )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!

Arie gravatar imageArie ( 2014-10-13 18:29:33 +0100 )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.

Dima gravatar imageDima ( 2014-10-14 09:16:57 +0100 )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.

Arie gravatar imageArie ( 2014-10-14 13:40:05 +0100 )edit
2

answered 2018-05-02 15:56:23 +0100

jipilab gravatar image

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(     x[0] == 1)
....: 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.

edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2014-10-10 18:53:00 +0100

Seen: 1,394 times

Last updated: May 02 '18