From 38650587bd385423d058ffe6810c43c902b48e98 Mon Sep 17 00:00:00 2001 From: mforets Date: Sun, 20 Jan 2019 22:48:53 -0300 Subject: [PATCH 1/9] change default intersection algorithm --- src/concrete_intersection.jl | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl index 0ca1a20835..729e573de8 100644 --- a/src/concrete_intersection.jl +++ b/src/concrete_intersection.jl @@ -364,23 +364,33 @@ A new same polytope in H-representation with just one more constraint. """ function intersection(P1::HPoly{N}, P2::HPoly{N}; - backend=default_polyhedra_backend(P1, N), - prunefunc=removehredundancy!) where {N<:Real} + backend=nothing, + prunefunc=remove_redundant_constraints!) where {N<:Real} if typeof(P1) == typeof(P2) HPOLY = typeof(P1) else # one of them must be a polytope, so the intersection will be bounded HPOLY = HPolytope{N} end + # concatenate the linear constraints Q = HPOLY([constraints_list(P1); constraints_list(P2)]) # remove redundant constraints - if prunefunc == removehredundancy! + if prunefunc == remove_redundant_constraints! + if backend == nothing + backend = GLPKSolverLP() + end + # here, detection of empty intersection may be reported as an infeasible LP + remove_redundant_constraints!(Q, backend=backend) + elseif prunefunc == removehredundancy! + if backend == nothing + backend = default_polyhedra_backend(P1, N) + end # convert to polyhedron ph = polyhedron(Q; backend=backend) - prunefunc(ph) + removehredundancy!(ph) Q = convert(HPOLY, ph) else prunefunc(Q) From d991221b2e757a4cb6dc744937221f89250a239a Mon Sep 17 00:00:00 2001 From: mforets Date: Sun, 20 Jan 2019 22:52:11 -0300 Subject: [PATCH 2/9] error for empty polyhedron --- src/HPolyhedron.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 66e5a99f76..b9983a43cd 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -484,8 +484,10 @@ function remove_redundant_constraints!(P::PT; br = b[non_redundant_indices] br[i] = b[j] + one(N) lp = linprog(-α, Ar, '<', br, -Inf, Inf, backend) - if lp.status != :Optimal - error("LP is not optimal") + if lp.status == :Infeasible + throw(DomainError(P, "the polyhedron is empty")) + else + error("LP is not optimal; the status of the LP is $(lp.status)") end objval = -lp.objval if objval <= b[j] From a4cffed3fa7d477b908a4b44164f4d7f3c1ec6c2 Mon Sep 17 00:00:00 2001 From: mforets Date: Sun, 20 Jan 2019 23:01:57 -0300 Subject: [PATCH 3/9] use relaxed _leq --- src/HPolyhedron.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index b9983a43cd..0ff541ac9c 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -490,7 +490,7 @@ function remove_redundant_constraints!(P::PT; error("LP is not optimal; the status of the LP is $(lp.status)") end objval = -lp.objval - if objval <= b[j] + if _leq(objval, b[j]) # the constraint is redundant non_redundant_indices = setdiff(non_redundant_indices, j) else From 938a8988e42bcdccb603316b45b0c460f03aced8 Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 22 Jan 2019 06:36:00 -0300 Subject: [PATCH 4/9] update branch --- src/HPolyhedron.jl | 53 ++++++++++++++++++-------------- src/concrete_intersection.jl | 58 +++++++++++++++++++++--------------- 2 files changed, 64 insertions(+), 47 deletions(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 0ff541ac9c..354c9f13aa 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -418,10 +418,9 @@ end """ remove_redundant_constraints(P::PT; - backend=GLPKSolverLP()) where {N, PT<:HPoly{N}} + backend=GLPKSolverLP())::Union{PT, EmptySet} where {N, PT<:HPoly{N}} -Given a polyhedron in H-representation, return a new polyhedron with no reundant -constraints. +Remove the redundant constraints in a polyhedron in H-representation. ### Input @@ -430,25 +429,29 @@ constraints. ### Output -A new polyhedron obtained by removing the redundant constraints in `P`. +A polyhedron equivalent to `P` but with no redundant constraints, or an empty set +if `P` is detected to be empty (which happens if the constraints are infeasible). ### Algorithm -See `remove_redundant_constraints!`. +See [`remove_redundant_constraints!`](@ref) for details. """ function remove_redundant_constraints(P::PT; - backend=GLPKSolverLP() - ) where {N, PT<:HPoly{N}} - return remove_redundant_constraints!(PT(copy(constraints_list(P))), - backend=backend) + backend=GLPKSolverLP())::Union{PT, EmptySet} where {N, PT<:HPoly{N}} + Pred = copy(P) + if remove_redundant_constraints!(Pred, backend=backend) + return Pred + else # the polyhedron P is empty + return EmptySet{N}() + end end """ remove_redundant_constraints!(P::PT; - backend=GLPKSolverLP()) where {N, PT<:HPoly{N}} + backend=GLPKSolverLP())::Bool where {N, PT<:HPoly{N}} Remove the redundant constraints in a polyhedron in H-representation; the polyhedron -is updated inplace. +is updated in-place. ### Input @@ -457,7 +460,9 @@ is updated inplace. ### Output -The polyhedron obtained by removing the redundant constraints in `P`. +`true` if the method was successful and the polyhedron `P` is modified by +removing its redundant constraints, and `false` if `P` is detected to be empty +(which happens if the constraints are infeasible). ### Algorithm @@ -470,7 +475,7 @@ For details, see [Fukuda's Polyhedra FAQ](https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node24.html). """ function remove_redundant_constraints!(P::PT; - backend=GLPKSolverLP()) where {N, PT<:HPoly{N}} + backend=GLPKSolverLP())::Bool where {N, PT<:HPoly{N}} A, b = tosimplehrep(P) m, n = size(A) @@ -485,22 +490,24 @@ function remove_redundant_constraints!(P::PT; br[i] = b[j] + one(N) lp = linprog(-α, Ar, '<', br, -Inf, Inf, backend) if lp.status == :Infeasible - throw(DomainError(P, "the polyhedron is empty")) + # the polyhedron is empty + return false + elseif lp.status == :Optimal + objval = -lp.objval + if _leq(objval, b[j]) + # the constraint is redundant + non_redundant_indices = setdiff(non_redundant_indices, j) + else + # the constraint is not redundant + i = i+1 + end else error("LP is not optimal; the status of the LP is $(lp.status)") end - objval = -lp.objval - if _leq(objval, b[j]) - # the constraint is redundant - non_redundant_indices = setdiff(non_redundant_indices, j) - else - # the constraint is not redundant - i = i+1 - end end deleteat!(P.constraints, setdiff(1:m, non_redundant_indices)) - return P + return true end """ diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl index 729e573de8..d4206a35ae 100644 --- a/src/concrete_intersection.jl +++ b/src/concrete_intersection.jl @@ -341,31 +341,44 @@ end """ intersection(P1::HPoly{N}, P2::HPoly{N}; - backend=default_polyhedra_backend(P1, N), - prunefunc=removehredundancy!) where {N<:Real} + backend=nothing, + use_polyhedra_interface=false) where {N<:Real} Compute the intersection of two polyhedra in H-representation. ### Input -- `P1` -- polytope -- `P2` -- polytope -- `backend` -- (optional, default: `default_polyhedra_backend(P1, N)`) the - polyhedral computations backend, see - [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1) - for further information -- `prunefunc` -- (optional, default: `removehredundancy!`) function to - post-process the polytope after adding the additional - constraint +- `P1` -- polyhedron +- `P2` -- polyhedron +- `backend` -- (optional, default: `nothing`) the LP solver or the the + polyhedral computations backend; its value is set internally, + see see below in the Notes for details +- `use_polyhedra_interface` -- (optional, default: `false`) if `true`, use the + `Polyhedra` interface for the removal of constraints ### Output -A new same polytope in H-representation with just one more constraint. +A polyhedron resulting from the intersection of `P1` and `P2`, with the redundant +constraints removed, or an empty set if the intersection is empty. + +### Notes + +The default value of the backend is set internally and depends on whether the +Polyhedra backend is used or not. The default backends are `GLPKSolverLP()` +and `default_polyhedra_backend(P1, N)` respectively. + +Note that if `use_polyhedra_interface` is set to `true`, there is no guarantee +that the removal of constraints keep the set empty (see #1038 and Polyhedra#146), +so it is better to check for emptiness of intersection before using this function +if that case. + +The method implemented in this function can be used for any pair of sets that can +handle the `constraints_list` option. """ function intersection(P1::HPoly{N}, P2::HPoly{N}; backend=nothing, - prunefunc=remove_redundant_constraints!) where {N<:Real} + use_polyhedra_interface=false) where {N<:Real} if typeof(P1) == typeof(P2) HPOLY = typeof(P1) else @@ -374,26 +387,23 @@ function intersection(P1::HPoly{N}, end # concatenate the linear constraints - Q = HPOLY([constraints_list(P1); - constraints_list(P2)]) + Q = HPOLY([constraints_list(P1); constraints_list(P2)]) # remove redundant constraints - if prunefunc == remove_redundant_constraints! - if backend == nothing - backend = GLPKSolverLP() - end - # here, detection of empty intersection may be reported as an infeasible LP - remove_redundant_constraints!(Q, backend=backend) - elseif prunefunc == removehredundancy! + if use_polyhedra_interface if backend == nothing backend = default_polyhedra_backend(P1, N) end - # convert to polyhedron + # convert to an hrep, remove the redundancies and convert back to HPOLY ph = polyhedron(Q; backend=backend) removehredundancy!(ph) Q = convert(HPOLY, ph) else - prunefunc(Q) + if backend == nothing + backend = GLPKSolverLP() + end + # here, detection of empty intersection may be reported as an infeasible LP + remove_redundant_constraints!(Q, backend=backend) end return Q end From 5b101fad085eb766941f7278d8064b35aa1f5b5b Mon Sep 17 00:00:00 2001 From: Christian Schilling Date: Tue, 22 Jan 2019 09:33:18 -0300 Subject: [PATCH 5/9] Update src/concrete_intersection.jl Co-Authored-By: mforets --- src/concrete_intersection.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl index d4206a35ae..035483ec52 100644 --- a/src/concrete_intersection.jl +++ b/src/concrete_intersection.jl @@ -370,7 +370,7 @@ and `default_polyhedra_backend(P1, N)` respectively. Note that if `use_polyhedra_interface` is set to `true`, there is no guarantee that the removal of constraints keep the set empty (see #1038 and Polyhedra#146), so it is better to check for emptiness of intersection before using this function -if that case. +in that case. The method implemented in this function can be used for any pair of sets that can handle the `constraints_list` option. From e032ee5c7ce6b7c9f59e96d44c1e7cfc1885c2aa Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 22 Jan 2019 09:34:10 -0300 Subject: [PATCH 6/9] remove _leq for <= --- src/HPolyhedron.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 354c9f13aa..6550304a3f 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -494,7 +494,7 @@ function remove_redundant_constraints!(P::PT; return false elseif lp.status == :Optimal objval = -lp.objval - if _leq(objval, b[j]) + if objval <= b[j] # the constraint is redundant non_redundant_indices = setdiff(non_redundant_indices, j) else From 4dd9513d5921837a6668659bd09e2f51f4ea2602 Mon Sep 17 00:00:00 2001 From: Christian Schilling Date: Tue, 22 Jan 2019 09:50:32 -0300 Subject: [PATCH 7/9] Update src/HPolyhedron.jl Co-Authored-By: mforets --- src/HPolyhedron.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index 6550304a3f..cd7109b889 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -437,7 +437,7 @@ if `P` is detected to be empty (which happens if the constraints are infeasible) See [`remove_redundant_constraints!`](@ref) for details. """ function remove_redundant_constraints(P::PT; - backend=GLPKSolverLP())::Union{PT, EmptySet} where {N, PT<:HPoly{N}} + backend=GLPKSolverLP())::Union{PT, EmptySet{N}} where {N, PT<:HPoly{N}} Pred = copy(P) if remove_redundant_constraints!(Pred, backend=backend) return Pred From 88a59ac495d7c3a4581d316b853ba44df1ae9439 Mon Sep 17 00:00:00 2001 From: Christian Schilling Date: Tue, 22 Jan 2019 09:50:43 -0300 Subject: [PATCH 8/9] Update src/HPolyhedron.jl Co-Authored-By: mforets --- src/HPolyhedron.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HPolyhedron.jl b/src/HPolyhedron.jl index cd7109b889..b9aa080c53 100644 --- a/src/HPolyhedron.jl +++ b/src/HPolyhedron.jl @@ -418,7 +418,7 @@ end """ remove_redundant_constraints(P::PT; - backend=GLPKSolverLP())::Union{PT, EmptySet} where {N, PT<:HPoly{N}} + backend=GLPKSolverLP())::Union{PT, EmptySet{N}} where {N, PT<:HPoly{N}} Remove the redundant constraints in a polyhedron in H-representation. From 6e0f5f1c5f1a6894842ec2103d1a3a1643ef4eaa Mon Sep 17 00:00:00 2001 From: mforets Date: Tue, 22 Jan 2019 10:57:45 -0300 Subject: [PATCH 9/9] update intersections with halfspace --- src/concrete_intersection.jl | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/src/concrete_intersection.jl b/src/concrete_intersection.jl index 035483ec52..134cd3a1e2 100644 --- a/src/concrete_intersection.jl +++ b/src/concrete_intersection.jl @@ -311,31 +311,28 @@ Compute the intersection of a polytope in H-representation and a half-space. - `P` -- polytope - `hs` -- half-space -- `backend` -- (optional, default: `default_polyhedra_backend(P, N)`) the - polyhedral computations backend, see - [Polyhedra's documentation](https://juliapolyhedra.github.io/Polyhedra.jl/latest/installation.html#Getting-Libraries-1) - for further information -- `prunefunc` -- (optional, default: `removehredundancy!`) function to - post-process the polytope after adding the additional - constraint - +- `backend` -- (optional, default: `nothing`) the LP solver or the the + polyhedral computations backend; its value is set internally, + see see below in the Notes for details +- `use_polyhedra_interface` -- (optional, default: `false`) if `true`, use the + `Polyhedra` interface for the removal of constraints ### Output The same polytope in H-representation with just one more constraint. """ function intersection(P::HPoly{N}, hs::HalfSpace{N}; - backend=default_polyhedra_backend(P, N), - prunefunc=removehredundancy!) where {N<:Real} - return intersection(P, HPolyhedron([hs]), backend=backend, prunefunc=prunefunc) + backend=nothing, + use_polyhedra_interface=false) where {N<:Real} + return intersection(P, HPolyhedron([hs]), backend=backend, use_polyhedra_interface=use_polyhedra_interface) end # symmetric method function intersection(hs::HalfSpace{N}, P::HPoly{N}; - backend=default_polyhedra_backend(P, N), - prunefunc=removehredundancy!) where {N<:Real} - return intersection(P, hs, backend=backend, prunefunc=prunefunc) + backend=nothing, + use_polyhedra_interface=false) where {N<:Real} + return intersection(P, hs, backend=backend, use_polyhedra_interface=use_polyhedra_interface) end """