Merge pull request #1056 from JuliaReach/schillic/1049

#1049 - Default implementation of isempty for polyhedra

docs/src/lib/representations.md

@@ -439,7 +439,7 @@ addconstraint!(::HPoly{N}, ::LinearConstraint{N}) where {N<:Real}
 constraints_list(::HPoly{N}) where {N<:Real}
 tosimplehrep(::HPoly{N}) where {N<:Real}
 tohrep(::HPoly{N}) where {N<:Real}
-isempty(::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}}
 tovrep(::HPoly{N}) where {N<:Real}

src/HPolyhedron.jl

@@ -706,38 +706,70 @@ function singleton_list(P::HPolyhedron{N}) where {N<:Real}
 end
 
 """
-   isempty(P::HPoly{N}; [solver]=GLPKSolverLP())::Bool where {N<:Real}
+   isempty(P::HPoly{N}, witness::Bool=false;
+           [use_polyhedra_interface]::Bool=false, [solver]=GLPKSolverLP(),
+           [backend]=nothing
+          )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
 
 Determine whether a polyhedron is empty.
 
 ### Input
 
 - `P`       -- polyhedron
-- `backend` -- (optional, default: `default_polyhedra_backend(P, N)`)
-               the polyhedral computations backend
-- `solver`  -- (optional, default: `GLPKSolverLP()`) LP solver backend
+- `witness` -- (optional, default: `false`) compute a witness if activated
+- `use_polyhedra_interface` -- (optional, default: `false`) if `true`, we use
+               the `Polyhedra` interface for the emptiness test
+- `solver`  -- (optional, default: `GLPKSolverLP()`) LP-solver backend
+- `backend` -- (optional, default: `default_polyhedra_backend(P, N)`) backend
+               for polyhedral computations in `Polyhedra`
 
 ### Output
 
-`true` if and only if the constraints are inconsistent.
+* If `witness` option is deactivated: `true` iff ``P = ∅``
+* If `witness` option is activated:
+  * `(true, [])` iff ``P = ∅``
+  * `(false, v)` iff ``P ≠ ∅`` and ``v ∈ P``
 
-### Algorithm
+### Notes
 
-This function uses `Polyhedra.isempty` which evaluates the feasibility of the
-LP whose feasible set is determined by the set of constraints and whose
-objective function is zero.
+Witness production is not supported if `use_polyhedra_interface` is `true`.
 
-### Notes
+### Algorithm
 
-This implementation uses `GLPKSolverLP` as linear programming backend by
-default.
-"""
-function isempty(P::HPoly{N};
-                 backend=default_polyhedra_backend(P, N),
-                 solver=GLPKSolverLP())::Bool where {N<:Real}
-    @assert isdefined(@__MODULE__, :Polyhedra) "the function `isempty` needs the " *
-                                        "package 'Polyhedra' to be loaded"
-    return Polyhedra.isempty(polyhedron(P; backend=backend), solver)
+The algorithm sets up a feasibility LP for the constraints of `P`.
+If `use_polyhedra_interface` is `true`, we call `Polyhedra.isempty`.
+Otherwise, we set up the LP internally.
+"""
+function isempty(P::HPoly{N},
+                 witness::Bool=false;
+                 use_polyhedra_interface::Bool=false,
+                 solver=GLPKSolverLP(),
+                 backend=default_polyhedra_backend(P, N)
+                )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
+    if use_polyhedra_interface
+        @assert isdefined(@__MODULE__, :Polyhedra) "the function `isempty` " *
+            "with the given options requires the package 'Polyhedra'"
+        result = Polyhedra.isempty(polyhedron(P; backend=backend), solver)
+        if result
+            return witness ? (true, N[]) : true
+        elseif witness
+            error("witness production is not supported yet")
+        else
+            return false
+        end
+    else
+        A, b = tosimplehrep(P)
+        lbounds, ubounds = -Inf, Inf
+        sense = '<'
+        obj = zeros(N, size(A, 2))
+        lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
+        if lp.status == :Optimal
+            return witness ? (false, lp.sol) : false
+        elseif lp.status == :Infeasible
+            return witness ? (true, N[]) : true
+        end
+        error("LP returned status $(lp.status) unexpectedly")
+    end
 end
 
 convert(::Type{HPolytope}, P::HPolyhedron{N}) where {N<:Real} =

src/is_intersection_empty.jl

@@ -810,7 +810,7 @@ end
                           X::LazySet{N},
                           witness::Bool=false;
                           solver=GLPKSolverLP(method=:Simplex)
-                         ) where {N<:Real}
+                         )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
 
 Check whether two polyhedra do not intersect.
 
@@ -839,8 +839,7 @@ For `algorithm == "sufficient"`, witness production is not supported.
 
 ### Algorithm
 
-For `algorithm == "exact"`, we set up a feasibility LP for the union of
-constraints of `P` and `X`.
+For `algorithm == "exact"`, see @ref(isempty(P::HPoly{N}, ::Bool)).
 
 For `algorithm == "sufficient"`, we rely on the intersection check between the
 set `X` and each constraint in `P`.
@@ -851,7 +850,7 @@ function is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
                                witness::Bool=false;
                                solver=GLPKSolverLP(method=:Simplex),
                                algorithm="exact"
-                              ) where {N<:Real}
+                              )::Union{Bool, Tuple{Bool, Vector{N}}} where {N<:Real}
     if algorithm == "sufficient"
         # sufficient check for empty intersection using half-space checks
         for Hi in constraints_list(P)
@@ -868,21 +867,8 @@ function is_intersection_empty(P::Union{HPolyhedron{N}, AbstractPolytope{N}},
         return false
     elseif algorithm == "exact"
         # exact check for empty intersection using a feasibility LP
-        A1, b1 = tosimplehrep(P)
-        A2, b2 = tosimplehrep(X)
-        A = [A1; A2]
-        b = [b1; b2]
-
-        lbounds, ubounds = -Inf, Inf
-        sense = '<'
-        obj = zeros(N, size(A, 2))
-        lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
-        if lp.status == :Optimal
-            return witness ? (false, lp.sol) : false
-        elseif lp.status == :Infeasible
-            return witness ? (true, N[]) : true
-        end
-        error("LP returned status $(lp.status) unexpectedly")
+        return isempty(HPolyhedron([constraints_list(P); constraints_list(X)]),
+                       witness; solver=solver)
     else
         error("algorithm $algorithm unknown")
     end