Polyhedron

As seen in the previous section, a polyhedron can be described in 2 ways: either using the H-representation (intersection of halfspaces) or the V-representation (convex hull of points and rays). The problem of computing the H-representation from the V-representation (or vice versa) is called the representation conversion problem. It can be solved by the Double-Description method

Polyhedra.doubledescriptionFunction
doubledescription(h::HRepresentation)

Computes the V-representation of the polyhedron represented by h using the Double-Description algorithm [MRTT53, FP96]. It maintains a list of the points, rays and lines that are in the resulting representation and in the hyperplane corresponding to each halfspace and then add their intersection for each halfspace with Polyhedra.add_intersection!. The resulting V-representation has no redundancy.

doubledescription(V::VRepresentation)

Computes the H-representation of the polyhedron represented by v using the Double-Description algorithm [MRTT53, FP96]. Currently, this fallbacks to lifting it to the V-representation of a cone by homogenization, interpreting it as the H-representation of the dual cone so that it can rely on doubledescription(::HRepresentation).

[MRTT53] Motzkin, T. S., Raiffa, H., Thompson, G. L. and Thrall, R. M. The double description method Contribution to the Theory of Games, Princeton University Press, 1953

[FP96] Fukuda, K. and Prodon, A. Double description method revisited Combinatorics and computer science, Springer, 1996, 91-111

source
Polyhedra.add_intersection!Function
add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)

Add the intersection of the elements of indices idx1 and idx2 that belong to the hyperplane corresponding the halfspace of index hp_idx (see 3.2 (ii) of [FP96]) or have inherited adjacency if hp_idx === nothing (see 3.2 (i) of [FP96]). In case of inherited adjacency, hs_idx is the halfspace where idx1 was created. If idx1 and idx2 are not adjacent or if they are in the hyperplane corresponding to a halfspace of lower index then the intersection is not added to avoid adding redundant elements.

[FP96] Fukuda, K. and Prodon, A. Double description method revisited Combinatorics and computer science, Springer, 1996, 91-111

source

However, other methods exist such as the reverse search implemented by LRS and the quick hull algorithm implemented by qhull.

This motivates the creation of a type representing polyhedra, transparently handling the conversion from H-representation to V-representation when needed for some operation. Just like the abstract type AbstractArray{N,T} represents an N-dimensional array with elements of type T, the abstract type Polyhedron{N,T} represents an N-dimensional polyhedron with elements of coefficient type T.

There is typically one concrete subtype of Polyhedron by library. For instance, the CDD library defines CDDLib.Polyhedron and the LRS library defines LRSLib.Polyhedron. It must be said that the type T is not necessarily how the elements are stored internally by the library but the polyhedron will behave just like it is stored that way. For instance, when retreiving an H-(or V-)representation, the representation will be of type T. Therefore using Int for T may result in InexactError. For this reason, by default, the type T chosen is not a subtype of Integer.

A polyhedron can be created from a representation and a library using the polyhedron function.

Polyhedra.polyhedronFunction
polyhedron(rep::Representation{T})

Creates a polyhedron from the representation rep using the default library included in the Polyhedra package.

source

To illustrate the usage of the polyhedron function, consider the following representations:

hr = HalfSpace([1, 1], 1) ∩ HalfSpace([1, -1], 0) ∩ HalfSpace([-1, 0], 0)
vre = convexhull([0, 0], [0, 1], [1//2, 1//2])
vrf = convexhull([0, 0], [0, 1], [1/2, 1/2])

One can use the CDD library, to create an instance of a concrete subtype of Polyhedron as follows:

julia> using CDDLib
julia> polyf = polyhedron(hr, CDDLib.Library())
julia> typeof(polyhf)
CDDLib.CDDLib.Polyhedron{2,Float64}

We see that the library has choosen to deal with floating point arithmetic. This decision does not depend on the type of hr but only on the instance of CDDLib.Library given. CDDLib.Library creates CDDLib.Polyhedron of type either Float64 or Rational{BigInt}. One can choose the first one using CDDLib.Library(:float) and the second one using CDDLib.Library(:exact), by default it is :float.

julia> poly = polyhedron(hr, CDDLib.Library(:exact))
julia> typeof(poly)
CDDLib.Polyhedron{2,Rational{BigInt}}

The first polyhedron polyf can also be created from its V-representation using either of the 4 following lines:

julia> polyf = polyhedron(vrf, CDDLib.Library(:float))
julia> polyf = polyhedron(vrf, CDDLib.Library())
julia> polyf = polyhedron(vre,  CDDLib.Library(:float))
julia> polyf = polyhedron(vre,  CDDLib.Library())

and poly using either of those lines:

julia> poly = polyhedron(vrf, CDDLib.Library(:exact))
julia> poly = polyhedron(vre, CDDLib.Library(:exact))

Of course, creating a representation in floating points with exact arithmetic works here because we have 0.5 which is 0.1 in binary but in general, is not a good idea.

julia> Rational{BigInt}(1/2)
1//2
julia> Rational{BigInt}(1/3)
6004799503160661//18014398509481984
julia> Rational{BigInt}(1/5)
3602879701896397//18014398509481984

Retrieving a representation

One can retrieve an H-representation (resp. V-representation) from a polyhedron using hrep (resp. vrep). The concrete subtype of HRepresentation (resp. VRepresentation) returned is not necessarily the same that the one used to create the polyhedron. As a rule of thumb, it is the representation the closest to the internal representation used by the library.

julia> hr = hrep(poly)
julia> typeof(hr)
Polyhedra.LiftedHRepresentation{2,Rational{BigInt}}
julia> hr = MixedMatHRep(hr)
julia> typeof(hr)
Polyhedra.MixedMatHRep{2,Rational{BigInt}}
julia> hr.A
3x2 Array{Rational{BigInt},2}:
  1//1   1//1
  1//1  -1//1
 -1//1   0//1
julia> hr.b
3-element Array{Rational{BigInt},1}:
 1//1
 0//1
 0//1
julia> vr = vrep(poly)
julia> typeof(vr)
Polyhedra.LiftedVRepresentation{2,Rational{BigInt}}
julia> vr = MixedMatVRep(vrep)
julia> typeof(vr)
Polyhedra.MixedMatVRep{2,Rational{BigInt}}
julia> vr.V
3x2 Array{Rational{BigInt},2}:
 1//2  1//2
 0//1  1//1
 0//1  0//1

julia> vr.R
0x2 Array{Rational{BigInt},2}

Checking if a representation has been computed

Incidence

Elements can be accessed in a representation or polyhedron using indices and Base.get:

Polyhedra.IndicesType
Indices{T, ElemT, RepT<:Rep{T}}

Iterator over the indices of the elements of type ElemT of the field rep.

source

The list of indices can be obtained using, e.g., eachindex(points(rep)). For instance, the following prints all points using indices

for pi in eachindex(points(rep))
    @show get(rep, pi)
end

A point $p$ (resp. ray $r$) is incident to an halfspace $\langle a, x \rangle \le \beta$ if $\langle a, p \rangle = \beta$ (resp. $\langle a, r \rangle = \beta$).

In a polyhedron, all points and rays are incident to all hyperplanes and all halfspaces are incident to all lines. The following methods are therefore redundant, e.g. incidenthyperplanes(p, idx) is equivalent to hyperplanes(p) and incidenthyperplaneindices(p, idx) is equivalent to eachindex(hyperplanes(p)). The methods are hence only defined for consistency.

Default libraries

The following functions allows to select a default library:

Polyhedra.default_libraryFunction
default_library(d::FullDim, ::Type{T}) where {T}

Returns the default polyhedral library for d-dimensional polyhedron of coefficient type T.

Examples

To obtain the default library for 2-dimensional polyhedra of eltype Float64, do default_library(2, Float64).

Given an StaticArrays.SVector v, to obtain a default library for points of the type of v in a type stable way, do default_library(Polyhedra.FullDim(v), eltype(v)).

source
Polyhedra.similar_libraryFunction
similar_library(lib::Library, d::FullDim, T::Type)

Returns a library that supports polyhedra of full dimension T with coefficient type T. If lib does not support it, this commonly calls default_library(d, T).

source
Polyhedra.default_typeFunction
default_type(d::FullDim, ::Type{T}) where {T}

Returns the default polyhedron type for d-dimensional polyhedron of coefficient type T.

source

The following libraries serves as fallback:

Polyhedra.DefaultLibraryType
DefaultLibrary{T}

Default library for polyhedra of dimension larger than 1 (IntervalLibrary is the default for polyhedra of dimension 1). The library implements the bare minimum and uses the fallback implementation for all operations.

source
Polyhedra.IntervalLibraryType
IntervalLibrary{T}

Default library for polyhedra of dimension 1. Many aspect of polyhedral computation become trivial in one dimension. This library exploits this fact. The library is also used as a fallback for libraries that do not support 1-dimensional polyhedra (e.g. qhull). That is projecting a polyhedron using such library produces a polyhedron using IntervalLibrary.

source

The type and library of the polyhedron obtained after applying an operation of several polyhedra (of possibly different type and/or library) is determined by the similar function.

Base.similarFunction
similar(p::Tuple{Vararg{Polyhedra.Rep}}, d::Polyhedra.FullDim, ::Type{T}, it::Polyhedra.It{T}...)

Creates a representation with a type similar to p of a polyhedron of full dimension d, element type T and initialize it with the iterators it. The type of the result will be chosen closer to the type of p[1].

source