#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)
 ```
 
 ### 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

@@ -213,7 +213,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}

src/is_intersection_empty.jl

@@ -529,3 +529,48 @@ 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 a pair of 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

test/unit_Polytope.jl

@@ -88,8 +88,15 @@ 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}
+
+# 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)
 
 # Polyhedra tests that only work with Float64
 if test_suite_polyhedra