#780 - Add AbstractPolyhedron interface

docs/src/lib/interfaces.md

@@ -100,10 +100,26 @@ an_element(::AbstractCentrallySymmetric{N}) where {N<:Real}
 isempty(::AbstractCentrallySymmetric)
 ```
 
-## Polytope
+## Polyhedron
 
-A polytope has finitely many vertices (*V-representation*) resp. facets
-(*H-representation*).
+A polyhedron has finitely many facets (*H-representation*) and is not
+necessarily bounded.
+
+```@docs
+AbstractPolyhedron
+```
+
+This interface defines the following functions:
+
+```@docs
+∈(::AbstractVector{N}, ::AbstractPolyhedron{N}) where {N<:Real}
+constrained_dimensions(::AbstractPolyhedron{N}) where {N<:Real}
+```
+
+### Polytope
+
+A polytope is a bounded set with finitely many vertices (*V-representation*)
+resp. facets (*H-representation*).
 Note that there is a special interface combination
 [Centrally symmetric polytope](@ref).
 
@@ -122,7 +138,7 @@ RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::AbstractPolytope)
 RecipesBase.apply_recipe(::Dict{Symbol,Any}, ::Vector{S}) where {S<:AbstractPolytope}
 ```
 
-### Polygon
+#### Polygon
 
 A polygon is a two-dimensional polytope.
 
@@ -137,7 +153,7 @@ dim(P::AbstractPolygon)
 linear_map(::AbstractMatrix{N}, P::AbstractPolygon{N}) where {N<:Real}
 ```
 
-#### HPolygon
+##### HPolygon
 
 An HPolygon is a polygon in H-representation (or constraint representation).
 
@@ -160,7 +176,7 @@ vertices_list(::AbstractHPolygon{N}, ::Bool=false, ::Bool=true) where {N<:Real}
 isbounded(::AbstractHPolygon, ::Bool=true)
 ```
 
-### Centrally symmetric polytope
+#### Centrally symmetric polytope
 
 A centrally symmetric polytope is a combination of two other interfaces:
 [Centrally symmetric set](@ref) and [Polytope](@ref).
@@ -177,7 +193,7 @@ an_element(::AbstractCentrallySymmetricPolytope{N}) where {N<:Real}
 isempty(::AbstractCentrallySymmetricPolytope)
 ```
 
-#### Hyperrectangle
+##### Hyperrectangle
 
 A hyperrectangle is a special centrally symmetric polytope with axis-aligned
 facets.
@@ -199,7 +215,7 @@ high(::AbstractHyperrectangle{N}) where {N<:Real}
 low(::AbstractHyperrectangle{N}) where {N<:Real}
 ```
 
-#### Singleton
+##### Singleton
 
 A singleton is a special hyperrectangle consisting of only one point.
 

docs/src/lib/representations.md

@@ -174,8 +174,10 @@ constraints_list(::HalfSpace{N}) where {N<:Real}
 constrained_dimensions(::HalfSpace{N}) where {N<:Real}
 halfspace_left(::AbstractVector{N}, ::AbstractVector{N}) where {N<:Real}
 halfspace_right(::AbstractVector{N}, ::AbstractVector{N}) where {N<:Real}
-tosimplehrep(::AbstractVector{HalfSpace{N}}) where {N<:Real}
 linear_map(::AbstractMatrix{N}, ::HalfSpace{N}) where {N}
+tosimplehrep(::AbstractVector{HalfSpace{N}}) where {N<:Real}
+remove_redundant_constraints(::AbstractVector{LinearConstraint{N}}) where {N<:Real}
+remove_redundant_constraints!(::AbstractVector{LinearConstraint{N}}) where {N<:Real}
 ```
 Inherited from [`LazySet`](@ref):
 * [`norm`](@ref norm(::LazySet, ::Real))
@@ -437,24 +439,27 @@ The following methods are shared between `HPolytope` and `HPolyhedron`.
 dim(::HPoly{N}) where {N<:Real}
 ρ(::AbstractVector{N}, ::HPoly{N}) where {N<:Real}
 σ(::AbstractVector{N}, ::HPoly{N}) where {N<:Real}
-∈(::AbstractVector{N}, ::HPoly{N}) where {N<:Real}
 addconstraint!(::HPoly{N}, ::LinearConstraint{N}) where {N<:Real}
 constraints_list(::HPoly{N}) where {N<:Real}
 tohrep(::HPoly{N}) where {N<:Real}
 isempty(::HPoly{N}, ::Bool=false) where {N<:Real}
 cartesian_product(::HPoly{N}, ::HPoly{N}) where {N<:Real}
-linear_map(M::AbstractMatrix{N}, P::PT) where {N<:Real, PT<:HPoly{N}}
+linear_map(::AbstractMatrix{N}, ::PT) where {N<:Real, PT<:HPoly{N}}
 tovrep(::HPoly{N}) where {N<:Real}
 polyhedron(::HPoly{N}) where {N<:Real}
-remove_redundant_constraints
-remove_redundant_constraints!
+remove_redundant_constraints(::PT) where {N<:Real, PT<:HPoly{N}}
+remove_redundant_constraints!(::HPoly{N}) where {N<:Real}
 ```
 
 Inherited from [`LazySet`](@ref):
 * [`norm`](@ref norm(::LazySet, ::Real))
 * [`radius`](@ref radius(::LazySet, ::Real))
 * [`diameter`](@ref diameter(::LazySet, ::Real))
 
+Inherited from [`AbstractPolyhedron`](@ref):
+* [`∈`](@ref ∈(::AbstractVector{N}, ::AbstractPolyhedron{N}) where {N<:Real})
+* [`constrained_dimensions`](@ref constrained_dimensions(::AbstractPolyhedron{N}) where {N<:Real})
+
 #### Polytopes in constraint representation
 
 The following methods are specific for `HPolytope`.
@@ -477,7 +482,6 @@ rand(::Type{HPolyhedron})
 isbounded(::HPolyhedron)
 vertices_list(::HPolyhedron{N}) where {N<:Real}
 singleton_list(::HPolyhedron{N}) where {N<:Real}
-constrained_dimensions(::HPolyhedron{N}) where {N<:Real}
 ```
 
 ### Vertex representation

src/AbstractPolyhedron.jl

@@ -0,0 +1,32 @@
+export AbstractPolyhedron
+
+"""
+    AbstractPolyhedron{N<:Real} <: LazySet{N}
+
+Abstract type for polyhedral sets, i.e., sets with finitely many flat facets, or
+equivalently, sets defined as an intersection of a finite number of halfspaces.
+
+### Notes
+
+Every concrete `AbstractPolyhedron` must define the following functions:
+- `constraints_list(::AbstractPolyhedron{N})::Vector{LinearConstraint{N}}` --
+    return a list of all facet constraints
+
+```jldoctest
+julia> subtypes(AbstractPolyhedron)
+5-element Array{Any,1}:
+ AbstractPolytope
+ HPolyhedron
+ HalfSpace
+ Hyperplane
+ Line
+```
+"""
+abstract type AbstractPolyhedron{N<:Real} <: LazySet{N} end
+
+
+# --- common AbstractPolyhedron functions ---
+
+
+# To account for the compilation order, functions are defined in the file
+# AbstractPolyhedron_functions.jl

src/AbstractPolyhedron_functions.jl

@@ -0,0 +1,215 @@
+import Base.∈
+
+export constrained_dimensions,
+       tosimplehrep,
+       remove_redundant_constraints,
+       remove_redundant_constraints!
+
+"""
+    ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N})::Bool where {N<:Real}
+
+Check whether a given point is contained in a polyhedron.
+
+### Input
+
+- `x` -- point/vector
+- `P` -- polyhedron
+
+### Output
+
+`true` iff ``x ∈ P``.
+
+### Algorithm
+
+This implementation checks if the point lies inside each defining half-space.
+"""
+function ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N})::Bool where {N<:Real}
+    @assert length(x) == dim(P) "a $(length(x))-dimensional point cannot be " *
+        "an element of a $(dim(P))-dimensional set"
+
+    for c in constraints_list(P)
+        if dot(c.a, x) > c.b
+            return false
+        end
+    end
+    return true
+end
+
+"""
+    constrained_dimensions(P::AbstractPolyhedron{N})::Vector{Int}
+        where {N<:Real}
+
+Return the indices in which a polyhedron is constrained.
+
+### Input
+
+- `P` -- polyhedron
+
+### Output
+
+A vector of ascending indices `i` such that the polyhedron is constrained in
+dimension `i`.
+
+### Examples
+
+A 2D polyhedron with constraint ``x1 ≥ 0`` is constrained in dimension 1 only.
+"""
+function constrained_dimensions(P::AbstractPolyhedron{N}
+                               )::Vector{Int} where {N<:Real}
+    zero_indices = zeros(Int, dim(P))
+    for constraint in constraints_list(P)
+        for i in constrained_dimensions(constraint)
+            zero_indices[i] = i
+        end
+    end
+    return filter(x -> x != 0, zero_indices)
+end
+
+"""
+    tosimplehrep(constraints::AbstractVector{LinearConstraint{N}})
+        where {N<:Real}
+
+Return the simple H-representation ``Ax ≤ b`` from a list of linear constraints.
+
+### Input
+
+- `constraints` -- a list of linear constraints
+
+### Output
+
+The tuple `(A, b)` where `A` is the matrix of normal directions and `b` is the
+vector of offsets.
+"""
+function tosimplehrep(constraints::AbstractVector{LinearConstraint{N}}
+                     ) where {N<:Real}
+    n = length(constraints)
+    if n == 0
+        A = Matrix{N}(undef, 0, 0)
+        b = Vector{N}(undef, 0)
+        return (A, b)
+    end
+    A = zeros(N, n, dim(first(constraints)))
+    b = zeros(N, n)
+    @inbounds begin
+        for (i, Pi) in enumerate(constraints)
+            A[i, :] = Pi.a
+            b[i] = Pi.b
+        end
+    end
+    return (A, b)
+end
+
+"""
+     remove_redundant_constraints!(
+         constraints::AbstractVector{LinearConstraint{N}};
+         [backend]=GLPKSolverLP())::Bool where {N<:Real}
+
+Remove the redundant constraints of a given list of linear constraints; the list
+is updated in-place.
+
+### Input
+
+- `constraints` -- list of constraints
+- `backend`     -- (optional, default: `GLPKSolverLP`) numeric LP solver backend
+
+### Output
+
+`true` if the function was successful and the list of constraints `constraints`
+is modified by removing the redundant constraints, and `false` only if the
+constraints are infeasible.
+
+### Notes
+
+Note that the result may be `true` even if the constraints are infeasible.
+For example, ``x ≤ 0 && x ≥ 1`` will return `true` without removing any
+constraint.
+To check if the constraints are infeasible, use
+`isempty(HPolyhedron(constraints)`.
+
+### Algorithm
+
+If there are `m` constraints in `n` dimensions, this function checks one by one
+if each of the `m` constraints is implied by the remaining ones.
+
+To check if the `k`-th constraint is redundant, an LP is formulated using the
+constraints that have not yet being removed.
+If, at an intermediate step, it is detected that a subgroup of the constraints
+is infeasible, this function returns `false`.
+If the calculation finished successfully, this function returns `true`.
+
+For details, see [Fukuda's Polyhedra
+FAQ](https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node24.html).
+"""
+function remove_redundant_constraints!(constraints::AbstractVector{LinearConstraint{N}};
+                                       backend=GLPKSolverLP()
+                                      )::Bool where {N<:Real}
+
+    A, b = tosimplehrep(constraints)
+    m, n = size(A)
+    non_redundant_indices = 1:m
+
+    i = 1 # counter over reduced constraints
+
+    for j in 1:m    # loop over original constraints
+        α = A[j, :]
+        Ar = A[non_redundant_indices, :]
+        br = b[non_redundant_indices]
+        br[i] = b[j] + one(N)
+        lp = linprog(-α, Ar, '<', br, -Inf, Inf, backend)
+        if lp.status == :Infeasible
+            # the polyhedron is empty
+            return false
+        elseif lp.status == :Optimal
+            objval = -lp.objval
+            if objval <= b[j]
+                # the constraint is redundant
+                non_redundant_indices = setdiff(non_redundant_indices, j)
+            else
+                # the constraint is not redundant
+                i = i+1
+            end
+        else
+            error("LP is not optimal; the status of the LP is $(lp.status)")
+        end
+    end
+
+    deleteat!(constraints, setdiff(1:m, non_redundant_indices))
+    return true
+end
+
+"""
+    remove_redundant_constraints(
+        constraints::AbstractVector{LinearConstraint{N}};
+        backend=GLPKSolverLP())::Union{AbstractVector{LinearConstraint{N}},
+                                       EmptySet{N}} where {N<:Real}
+
+Remove the redundant constraints of a given list of linear constraints.
+
+### Input
+
+- `constraints` -- list of constraints
+- `backend`     -- (optional, default: `GLPKSolverLP`) numeric LP solver backend
+
+### Output
+
+The list of constraints with the redundant ones removed, or an empty set if the
+constraints are infeasible.
+
+### Algorithm
+
+See
+[`remove_redundant_constraints!(::AbstractVector{LinearConstraint{<:Real}})`](@ref)
+for details.
+"""
+function remove_redundant_constraints(constraints::AbstractVector{LinearConstraint{N}};
+                                      backend=GLPKSolverLP()
+                                     )::Union{AbstractVector{LinearConstraint{N}},
+                                                             EmptySet{N}
+                                                            } where {N<:Real}
+    constraints_copy = copy(constraints)
+    if remove_redundant_constraints!(constraints_copy, backend=backend)
+        return constraints_copy
+    else  # the constraints are infeasible
+        return EmptySet{N}()
+    end
+end