[Haskell-cafe] Finding points contained within a convex hull.

[haskell at list.mightyreason.com]

Daniel McAllansmith wrote:
> Hello.
> 
> I've got a system of linear inequalities representing half-spaces.  The 
> half-spaces may or may not form a convex hull.

They could only fail to define a convex volume if they are inconsistent and
define an empty set.  Though they might define a convex volume of lesser
dimension than the space.

> 
> I need to find the integral coordinates that are contained within the convex 
> hull, if there is one.
> 
> For example, given
> 0 <= x <= 4
> 0 <= y <= 3
> 0 <= 2x - y
> 0 <= 1.2y - x
> 
> I want the following (x,y) coordinates
> [(0,0),(1,1),(1,2),(2,2),(2,3),(3,3)]
> 
> 
> Anybody have any suggestions, or libraries, for solving this in many 
> dimensions and equations?
> 
> 
> Thanks
> Daniel

This FAQ looks useful: http://www.cs.mcgill.ca/~fukuda/soft/polyfaq/polyfaq.html

You want to find something like the Volume, so Quoting from:
http://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node26.html

> 
> Is there any efficient algorithm to compute the volume of a convex polytope
> in $ R^d$?
> 
> It is known that computing the volume of a $ V$-polytope (or H-polytope) is
> #P-hard, see [DF88] and [Kha93]. There are theoretically efficient randomized
> algorithms to approximate the volume of a convex body [LS93] but no
> implementation seems to be available.
> 
> There is a comparative study [BEF00] of various volume computation algorithms
> for convex polytopes. It indicates that there is no single algorithm that
> works well for many different types of polytopes. For ``near'' simple
> polytopes, triangulation-based algorithms are more efficient. For ``near''
> simplicial polytopes, sign-decomposition-based algorithms are better. See the
> paper for the justification of these claims.
> 

I have never touched a problem like this, but it seems that enumerating the
point could be attacked by recursive brute force:

* Find an extreme vertex and any other point in the polyhedra (e.g. a different
extreme vertex).  If the polyhedra is empty this will fail.  If the polyhedra is
a single point this will find one vertex but no other point, so just check this
vertex to see if it has all integer coordinates.

* Pick a dimension for which the vertex and other point have different
coordinates. Call this dimension's coordinate x.  The vertex has value vx and
the other point ix.  Neither vx nor ix are required to be integers.

* Now you make several recursive problems that have different integer values of
coordinate x.  These subproblems are in two series:
  ** Increasing series x = a, a+1, a+2, a+3, ... where a = ceil of vx.
  ** Decreasing series x = b, b-1, b-2, b-3, ... where b = floor of vx.
  ** Obviously a==b is a special case.
  ** If you are very clever you can always find a dimensions for which it is
obvious that one of the two series is always outside the polygon and need not be
checked.  But this will also be noticed by the brute force algorithm once that
series is run for the first x not equal to vx (which will be either the first or
second element of that series depending on whether vx is an integer).

* For the two series plug 'x' back into your original set of equations.
  ** If any equation is False then stop this series of subproblems
  ** If an equation is True then you can remove it
  ** If all equations are True, then you have found a integer point in the
volume and you can return it and continue with the series of subproblems.
  ** Otherwise, take the remaining equation compute the vertices of the convex
hull (If the new polyhedra is empty then stop this series of subproblems) and
iterate this enumeration procedure.
  ** Note: A subproblem may find no integer points in its convex hull, but that
does not stop the algorithm from checking the next value of x in the series.

The above looks like a polynomial time complexity algorithm to me.


Cheers,
  Chris