Merge pull request #3235 from JuliaReach/schillic/1078

#1078 - Generalize isempty from HPolyhedron to (polyhedral) LazySet

docs/src/lib/interfaces.md

@@ -77,6 +77,7 @@ is_polyhedral(::LazySet)
 norm(::LazySet, ::Real=Inf)
 radius(::LazySet, ::Real=Inf)
 diameter(::LazySet, ::Real=Inf)
+isempty(::LazySet{N}, ::Bool=false) where {N}
 affine_map(::AbstractMatrix, ::LazySet, ::AbstractVector)
 exponential_map(::AbstractMatrix, ::LazySet)
 an_element(::LazySet)

docs/src/lib/sets/HPolyhedron.md

@@ -19,7 +19,6 @@ constraints_list(::HPoly)
 tohrep(::HPoly)
 tovrep(::HPoly)
 normalize(::HPoly{N}, p=N(2)) where {N}
-isempty(::HPoly{N}, ::Bool=false) where {N}
 translate(::HPoly, ::AbstractVector)
 remove_redundant_constraints(::HPoly)
 remove_redundant_constraints!(::HPoly)
@@ -31,6 +30,7 @@ Inherited from [`LazySet`](@ref):
 * [`radius`](@ref radius(::LazySet, ::Real))
 * [`diameter`](@ref diameter(::LazySet, ::Real))
 * [`singleton_list`](@ref singleton_list(::LazySet))
+* [`isempty`](@ref isempty(::LazySet{N}, ::Bool=false) where {N})
 
 Inherited from [`AbstractPolyhedron`](@ref):
 * [`∈`](@ref ∈(::AbstractVector, ::AbstractPolyhedron))

src/ConcreteOperations/isdisjoint.jl

@@ -868,7 +868,7 @@ function _isdisjoint_polyhedron(P::AbstractPolyhedron, X::LazySet,
         if solver == nothing
             solver = default_lp_solver(N)
         end
-        return isempty(HPolyhedron([clist_P; clist_X]), witness; solver=solver)
+        return _isempty_polyhedron_lp([clist_P; clist_X], witness; solver=solver)
     else
         error("algorithm $algorithm unknown")
     end

src/Interfaces/AbstractPolyhedron_functions.jl

@@ -996,7 +996,7 @@ feasible solution: ``\\min∥y∥_1`` subject to ``A^Ty=0`` and ``y≥1``.
 function _isbounded_stiemke(constraints::AbstractVector{<:HalfSpace{N}};
                             solver=LazySets.default_lp_solver(N),
                             check_nonempty::Bool=true) where {N}
-    if check_nonempty && isempty(HPolyhedron(constraints))
+    if check_nonempty && _isempty_polyhedron_lp(constraints; solver=solver)
         return true
     end
 
@@ -1199,3 +1199,28 @@ function _project_polyhedron(P::LazySet{N}, block; kwargs...) where {N}
     πP = linear_map(M, P; kwargs...)
     return constraints_list(πP)
 end
+
+function _isempty_polyhedron_lp(constraints::AbstractVector{<:HalfSpace},
+                                witness::Bool=false; solver=nothing)
+    # the caller should verify that there is at least one constraint
+    A, b = tosimplehrep(constraints)
+    return _isempty_polyhedron_lp(A, b, witness; solver=solver)
+end
+
+function _isempty_polyhedron_lp(A::AbstractMatrix{N}, b::AbstractVector{N},
+                                witness::Bool=false; solver=nothing) where {N}
+    # feasibility LP
+    lbounds, ubounds = -Inf, Inf
+    sense = '<'
+    obj = zeros(N, size(A, 2))
+    if isnothing(solver)
+        solver = default_lp_solver(N)
+    end
+    lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
+    if is_lp_optimal(lp.status)
+        return witness ? (false, lp.sol) : false
+    elseif is_lp_infeasible(lp.status)
+        return witness ? (true, N[]) : true
+    end
+    error("LP returned status $(lp.status) unexpectedly")
+end

src/Interfaces/LazySet.jl

@@ -1,4 +1,4 @@
-import Base: extrema, ==, ≈, copy, eltype, rationalize
+import Base: isempty, extrema, ==, ≈, copy, eltype, rationalize
 import Random.rand
 
 export LazySet,
@@ -24,6 +24,7 @@ export LazySet,
        exponential_map,
        reflect,
        is_interior_point,
+       isempty,
        isoperation,
        isoperationtype,
        isequivalent,
@@ -1741,3 +1742,98 @@ function triangulate(X::LazySet)
 end
 
 end end  # quote / load_polyhedra_lazyset()
+
+"""
+    isempty(P::LazySet{N}, witness::Bool=false;
+            [use_polyhedra_interface]::Bool=false, [solver]=nothing,
+            [backend]=nothing) where {N}
+
+Check whether a polyhedral set is empty.
+
+### Input
+
+- `P`       -- polyhedral set
+- `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: `nothing`) LP-solver backend; uses
+               `default_lp_solver(N)` if not provided
+- `backend` -- (optional, default: `nothing`) backend for polyhedral
+               computations in `Polyhedra`; uses `default_polyhedra_backend(P)`
+               if not provided
+
+### Output
+
+* If `witness` option is deactivated: `true` iff ``P = ∅``
+* If `witness` option is activated:
+  * `(true, [])` iff ``P = ∅``
+  * `(false, v)` iff ``P ≠ ∅`` and ``v ∈ P``
+
+### Notes
+
+The default value of the `backend` is set internally and depends on whether the
+`use_polyhedra_interface` option is set or not.
+If the option is set, we use `default_polyhedra_backend(P)`.
+
+Witness production is not supported if `use_polyhedra_interface` is `true`.
+
+### Algorithm
+
+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::LazySet{N},
+                 witness::Bool=false;
+                 use_polyhedra_interface::Bool=false,
+                 solver=nothing,
+                 backend=nothing) where {N}
+    @assert is_polyhedral(P) "this algorithm requires a polyhedral set"
+    clist = constraints_list(P)
+    if length(clist) < 2
+        # catch corner case because of problems in LP solver for Rationals
+        return witness ? (false, an_element(P)) : false
+    end
+    if use_polyhedra_interface
+        return _isempty_polyhedron_polyhedra(P, witness; solver=solver,
+                                             backend=backend)
+    else
+        return _isempty_polyhedron_lp(clist, witness; solver=solver)
+    end
+end
+
+function _isempty_polyhedron(P::LazySet{N}, witness::Bool=false;
+                             use_polyhedra_interface::Bool=false,
+                             solver=nothing, backend=nothing) where {N}
+    if use_polyhedra_interface
+        return _isempty_polyhedron_polyhedra(P, witness; solver=solver,
+                                             backend=backend)
+    else
+        return _isempty_polyhedron_lp(constraints_list(P), witness;
+                                      solver=solver)
+    end
+end
+
+function _isempty_polyhedron_polyhedra(P::LazySet{N}, witness::Bool=false;
+                                       solver=nothing, backend=nothing) where {N}
+    require(@__MODULE__, :Polyhedra; fun_name="isempty",
+            explanation="with the active option `use_polyhedra_interface`")
+
+    if backend == nothing
+        backend = default_polyhedra_backend(P)
+    end
+
+    if isnothing(solver)
+        result = Polyhedra.isempty(polyhedron(P; backend=backend))
+    else
+        result = Polyhedra.isempty(polyhedron(P; backend=backend), solver)
+    end
+
+    if result
+        return witness ? (true, N[]) : true
+    elseif witness
+        error("witness production is not supported yet")
+    else
+        return false
+    end
+end

src/LazyOperations/IntersectionArray.jl

@@ -139,6 +139,10 @@ function isboundedtype(::Type{<:IntersectionArray{N, S}}) where {N, S}
     return isboundedtype(S)
 end
 
+function is_polyhedral(ia::IntersectionArray)
+    return all(is_polyhedral, array(ia))
+end
+
 """
     ∈(x::AbstractVector, ia::IntersectionArray)
 

src/Sets/HPolyhedron.jl

@@ -1,5 +1,4 @@
-import Base: isempty,
-             rand,
+import Base: rand,
              convert
 
 export HPolyhedron,
@@ -9,7 +8,6 @@ export HPolyhedron,
        tohrep, tovrep,
        convex_hull,
        cartesian_product,
-       isempty,
        remove_redundant_constraints,
        remove_redundant_constraints!,
        constrained_dimensions,
@@ -454,92 +452,21 @@ function tovrep(P::HPoly; backend=default_polyhedra_backend(P))
     return VPolytope(P)
 end
 
-"""
-    isempty(P::HPoly{N}, witness::Bool=false;
-            [use_polyhedra_interface]::Bool=false, [solver]=nothing,
-            [backend]=nothing) where {N}
-
-Check whether a polyhedron in constraint representation is empty.
-
-### Input
-
-- `P`       -- polyhedron in constraint representation
-- `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: `nothing`) LP-solver backend; uses
-               `default_lp_solver(N)` if not provided
-- `backend` -- (optional, default: `nothing`) backend for polyhedral
-               computations in `Polyhedra`; uses `default_polyhedra_backend(P)`
-               if not provided
-
-### Output
-
-* If `witness` option is deactivated: `true` iff ``P = ∅``
-* If `witness` option is activated:
-  * `(true, [])` iff ``P = ∅``
-  * `(false, v)` iff ``P ≠ ∅`` and ``v ∈ P``
-
-### Notes
-
-The default value of the `backend` is set internally and depends on whether the
-`use_polyhedra_interface` option is set or not.
-If the option is set, we use `default_polyhedra_backend(P)`.
-
-Witness production is not supported if `use_polyhedra_interface` is `true`.
-
-### Algorithm
-
-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},
+# this method is required mainly for HPolytope (because the fallback for
+# AbstractPolytope is incorrect with no constraints)
+#
+# the method also treats a corner case for problems with Rationals in LP solver
+function isempty(P::HPoly,
                  witness::Bool=false;
                  use_polyhedra_interface::Bool=false,
                  solver=nothing,
-                 backend=nothing) where {N}
+                 backend=nothing)
     if length(constraints_list(P)) < 2
-        # catch corner case because of problems in LP solver for Rationals
         return witness ? (false, an_element(P)) : false
     end
-    if use_polyhedra_interface
-        require(@__MODULE__, :Polyhedra; fun_name="isempty",
-                explanation="with the active option `use_polyhedra_interface`")
-
-        if backend == nothing
-            backend = default_polyhedra_backend(P)
-        end
-
-        if isnothing(solver)
-            result = Polyhedra.isempty(polyhedron(P; backend=backend))
-        else
-            result = Polyhedra.isempty(polyhedron(P; backend=backend), solver)
-        end
-
-        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))
-        if isnothing(solver)
-            solver = default_lp_solver(N)
-        end
-        lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
-        if is_lp_optimal(lp.status)
-            return witness ? (false, lp.sol) : false
-        elseif is_lp_infeasible(lp.status)
-            return witness ? (true, N[]) : true
-        end
-        error("LP returned status $(lp.status) unexpectedly")
-    end
+    return _isempty_polyhedron(P, witness;
+                               use_polyhedra_interface=use_polyhedra_interface,
+                               solver=solver, backend=backend)
 end
 
 function load_polyhedra_hpolyhedron() # function to be loaded by Requires

test/LazyOperations/Intersection.jl

@@ -98,8 +98,15 @@ for N in [Float64, Rational{Int}, Float32]
     IArr2 = IntersectionArray([HalfSpace(N[1], N(1))])
     @test !isboundedtype(typeof(IArr2))
 
+    # is_polyhedral
+    @test is_polyhedral(IArr)
+    if N isa AbstractFloat
+        IArr2 = IntersectionArray([B, Ball2(N[0, 0], N(1))])
+        @test !is_polyhedral(IArr2)
+    end
+
     # isempty
-    @test_throws MethodError isempty(IArr)
+    @test !isempty(IArr)
 
     # membership
     @test ones(N, 2) ∈ IArr && N[5, 5] ∉ IArr