From bdba36ea98e175c8c4b4960e79af65871ddecfe7 Mon Sep 17 00:00:00 2001 From: schillic Date: Tue, 22 Jan 2019 21:02:00 +0100 Subject: [PATCH] refactoring: move code from isdisjoint to isempty --- docs/src/lib/representations.md | 2 +- src/HPolyhedron.jl | 70 ++++++++++++++++++++++++--------- src/is_intersection_empty.jl | 24 +++-------- 3 files changed, 57 insertions(+), 39 deletions(-) diff --git a/docs/src/lib/representations.md b/docs/src/lib/representations.md index 817d1c435e..415699181e 100644 --- a/docs/src/lib/representations.md +++ b/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} diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index b9aa080c53..b375dba3c2 100644 --- a/src/HPolyhedron.jl +++ b/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} = diff --git a/src/is_intersection_empty.jl b/src/is_intersection_empty.jl index 4ed97bc141..db892f22e2 100644 --- a/src/is_intersection_empty.jl +++ b/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