How to find integer points (coordinates) inside a polytope/polyhedra?

Question

I am using Python but I wouldn't mind changing language. All I have gotten from my research are tools to count the number of (lattice) points inside a region given the equations for the planes that enclose it. Other tools are made to optimize a given function inside the polytope (linear programming).

How about finding the lattice points alone? For example, a function of the kind

latticePoints( 'x < 5 & x > 0' ) = [ 1, 2, 3, 4]

Plus I am looking for something to work in the multivariate scenario (constrains on x, y, z, ...).

I am currently trying to solve this using ppl.

Answer 1

Using the Python package polytope, the integral points within a d-dimensional polytope can be computed as follows (this script is based on a test that I wrote: (polytope_test.py lines 415--455):

"""How to compute all points with integer coordinates inside a polytope."""
import numpy as np
import polytope.polytope as alg


def example():
    """Demonstrate the integral points computation."""
    # convex polytope
    vertices = np.array([[0.5, 1.5], [0.5, 1.5]])
    hull = alg.box2poly(vertices)
        # `hull` is an instance of the class `polytope.polytope.Polytope`,
        # which is for representing convex polytopes
    integral_points = alg.enumerate_integral_points(hull)
    print(hull)
    print('contains the integral points:')
    print(integral_points)
    #
    # nonconvex polytope
    vertices = np.array([[0.0, 0.0], [1.0, 1.0], [2.0, 1.0]])
    hull_1 = alg.qhull(vertices)  # convex hull of vertices in `vertices`
    hull_2 = alg.box2poly([[1.0, 2.0], [1.0, 2.0]])
    nonconvex = hull_1.union(hull_2)
        # `nonconvex` is an instance of the class `polytope.polytope.Region`,
        # which is for representing any polytope, including nonconvex ones,
        # and in this case can also be constructed with
        # `polytope.polytope.Region([hull_1, hull_2])`
    integral_points = alg.enumerate_integral_points(nonconvex)
    print('The polytope that is the union of the following polytopes:')
    print(nonconvex)
    print('contains the integral points:')
    print(integral_points)
    #
    # 3-dimensional polytope
    vertices = np.array([
        [0.0, 0.0, 0.0],
        [1.0, 0.0, 0.0],
        [0.0, 1.0, 0.0],
        [0.0, 0.0, 1.0]])
    hull = alg.qhull(vertices)
    integral_points = alg.enumerate_integral_points(hull)
    print(hull)
    print('contains the integral points:')
    print(integral_points)


if __name__ == '__main__':
    example()

Currently, the above Python code works with the development version of polytope, which can be installed either using the package installer pip:

pip install git+git://github.com/tulip-control/polytope.git

or by cloning the GitHub repository, and installing from the cloned repository:

git clone [email protected]:tulip-control/polytope
cd polytope
pip install .

The Python script above outputs:

Single polytope
  [[ 1.  0.] |    [[ 1.5]
   [ 0.  1.] x <=  [ 1.5]
   [-1. -0.] |     [-0.5]
   [-0. -1.]]|     [-0.5]]

contains the integral points:
[[1.]
 [1.]]
The polytope that is the union of the following polytopes:
     Polytope number 1:
     Single polytope
          [[-0.70711  0.70711] |    [[0.]
           [ 0.       1.     ] x <=  [1.]
           [ 0.44721 -0.89443]]|     [0.]]

     Polytope number 2:
     Single polytope
          [[ 1.  0.] |    [[ 2.]
          [ 0.  1.] x <=  [ 2.]
          [-1. -0.] |     [-1.]
          [-0. -1.]]|     [-1.]]



contains the integral points:
[[0. 1. 2. 1. 2.]
 [0. 1. 1. 2. 2.]]
Single polytope
  [[ 0.      -1.      -0.     ] |    [[0.     ]
   [-1.      -0.      -0.     ] x <=  [0.     ]
   [ 0.       0.      -1.     ] |     [0.     ]
   [ 0.57735  0.57735  0.57735]]|     [0.57735]]

contains the integral points:
[[0. 0. 1. 0.]
 [0. 0. 0. 1.]
 [0. 1. 0. 0.]]

Answer 2

There is a nice answer in Mathematica here:

points = {x, y} /. List@ToRules@ Reduce[x >= 4 y && x <= 4 y + 3 && 0 < x < 63 && 0 < y < 15, {x, y}, Integers]

---

         

---