#619 Is intersection empty for polytopes

docs/src/lib/binary_functions.md

@@ -44,6 +44,7 @@ is_intersection_empty(::Zonotope{N}, ::Hyperplane{N}, ::Bool=false) where {N<:Re
 is_intersection_empty(::Hyperplane{N}, ::Zonotope{N}, ::Bool=false) where {N<:Real}
 is_intersection_empty(::Ball2{N}, ::Ball2{N}, ::Bool=false) where {N<:AbstractFloat}
 is_intersection_empty(::LineSegment{N}, ::LineSegment{N}, ::Bool=false) where {N<:Real}
+is_intersection_empty(::AbstractPolytope{N}, ::AbstractPolytope{N}, ::Bool=false) where {N<:Real}
 ```
 
 ## Intersection of two sets

docs/src/lib/interfaces.md

@@ -93,6 +93,7 @@ This interface defines the following functions:
 ```@docs
 singleton_list(::AbstractPolytope)
 linear_map(::AbstractMatrix, ::AbstractPolytope)
+isempty(::AbstractPolytope{N}) where {N<:Real}
 ```
 
 ### Polygon

src/AbstractPolytope.jl

@@ -1,7 +1,10 @@
+import Base.isempty
+
 export AbstractPolytope,
        vertices_list,
        singleton_list,
-       linear_map
+       linear_map,
+       isempty
 
 """
     AbstractPolytope{N<:Real} <: LazySet{N}
@@ -81,3 +84,25 @@ function linear_map(M::AbstractMatrix, P::AbstractPolytope{N}) where {N<:Real}
     end
     return T(new_vlist)
 end
+
+"""
+    isempty(P::AbstractPolytope{N})::Bool where {N<:Real}
+
+Determine whether a polytope is empty.
+
+### Input
+
+- `P` -- abstract polytope
+
+### Output
+
+`true` if the given polytope contains no vertices, and `false` otherwise.
+
+### Algorithm
+
+This algorithm checks whether the `vertices_list` of the given polytope is empty
+or not.
+"""
+function isempty(P::AbstractPolytope{N})::Bool where {N<:Real}
+    return isempty(vertices_list(P))
+end

src/HPolytope.jl

@@ -195,6 +195,12 @@ Return the simple H-representation ``Ax ≤ b`` of a polytope.
 The tuple `(A, b)` where `A` is the matrix of normal directions and `b` are the offsets.
 """
 function tosimplehrep(P::HPolytope{N}) where {N}
+    if length(P.constraints) == 0
+        A = Matrix{N}(undef, 0, 0)
+        b = Vector{N}(undef, 0)
+        return (A, b)
+    end
+
     A = hcat([ci.a for ci in P.constraints]...)'
     b = [ci.b for ci in P.constraints]
     return (A, b)
@@ -213,7 +219,13 @@ import Polyhedra:polyhedron, SimpleHRepresentation, SimpleHRepresentation,
                  hcartesianproduct,
                  points
 
-export intersection, convex_hull, cartesian_product, vertices_list, tovrep, tohrep
+export intersection,
+       convex_hull,
+       cartesian_product,
+       vertices_list,
+       tovrep,
+       tohrep,
+       is_intersection_empty
 
 # HPolytope from an HRep
 function HPolytope(P::HRep{N, T}, backend=CDDLib.CDDLibrary()) where {N, T}
@@ -366,6 +378,9 @@ julia> vertices_list(P)
 function vertices_list(P::HPolytope{N};
                        backend=CDDLib.CDDLibrary(),
                        prunefunc=removevredundancy!)::Vector{Vector{N}} where {N<:Real}
+    if length(P.constraints) == 0
+        return Vector{N}(undef, Vector{N}(undef, 0))
+    end
     P = polyhedron(P, backend)
     prunefunc(P)
     return collect(points(P))

src/is_intersection_empty.jl

@@ -529,3 +529,69 @@ function is_intersection_empty(ls1::LineSegment{N},
         return empty_intersection
     end
 end
+
+# --- AbstractPolytope ---
+
+"""
+    is_intersection_empty(P::AbstractPolytope{N},
+                          Q::AbstractPolytope{N},
+                          witness::Bool=false
+                          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+
+Check whether two polytopes do not intersect, and otherwise optionally
+compute a witness.
+
+### Input
+
+- `P`       -- polytope
+- `Q`       -- another polytope
+- `witness` -- (optional, default: `false`) compute a witness if activated
+
+### Output
+
+* If `witness` option is deactivated: `true` iff ``P ∩ Q = ∅``
+* If `witness` option is activated:
+  * `(true, [])` iff ``P ∩ Q = ∅``
+  * `(false, v)` iff ``P ∩ Q ≠ ∅`` and ``v ∈ P ∩ Q``
+
+### Algorithm
+
+This is a fallback implementation of the `AbstractPolytope` interface that computes
+the concrete intersection, `intersection`, of the given pair of polytopes. If a
+witness is required, the first vertex of the resulting intersection polytope is
+returned.
+"""
+function is_intersection_empty(P::AbstractPolytope{N},
+                               Q::AbstractPolytope{N},
+                               witness::Bool=false
+                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    X = intersection(P, Q)
+    isempty_flag = isempty(X)
+    if witness
+        witness_vertex = isempty_flag ? N[] : vertices_list(X)[1]
+        return (isempty_flag, witness_vertex)
+    else
+        return isempty_flag
+    end
+end
+
+# method disambiguation
+function is_intersection_empty(point::AbstractSingleton{N},
+                               set::AbstractPolytope{N},
+                               witness::Bool=false
+                              )::Union{Bool,Tuple{Bool,Vector{N}}} where {N<:Real}
+    empty_intersection = !∈(element(point), set)
+    if witness
+        return (empty_intersection, empty_intersection ? N[] : element(S))
+    else
+        return empty_intersection
+    end
+end
+
+# symmetric function
+function is_intersection_empty(set::AbstractPolytope{N},
+                               point::AbstractSingleton{N},
+                               witness::Bool=false
+                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    return is_intersection_empty(point, set, witness)
+end

test/unit_Polytope.jl

@@ -88,8 +88,17 @@ for N in [Float64, Rational{Int}, Float32]
 end
 
 # default Float64 constructors
-@test HPolytope() isa HPolytope{Float64}
-@test VPolytope() isa VPolytope{Float64}
+empty_H_polytope = HPolytope()
+empty_V_polytope = VPolytope()
+@test empty_H_polytope isa HPolytope{Float64}
+@test empty_V_polytope isa VPolytope{Float64}
+
+if test_suite_polyhedra
+    # check isempty function of the AbstractPolytope interface
+    @test isempty(empty_H_polytope)
+    @test isempty(empty_V_polytope)
+    is_intersection_empty(empty_H_polytope, empty_V_polytope)
+end
 
 # Polyhedra tests that only work with Float64
 if test_suite_polyhedra