Merge pull request #1040 from JuliaReach/mforets/polyhedron_constraints

polyhedron redundant constraints
JuliaReach · Jan 22, 2019 · 4847490 · 4847490
2 parents 82b0c12 + 6e0f5f1
commit 4847490
Show file tree

Hide file tree

Showing 2 changed files with 83 additions and 57 deletions.
diff --git 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{N}} 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{N}} 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)
@@ -484,21 +489,25 @@ 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")
-        end
-        objval = -lp.objval
-        if objval <= b[j]
-            # the constraint is redundant
-            non_redundant_indices = setdiff(non_redundant_indices, j)
+        if lp.status == :Infeasible
+            # the polyhedron is empty
+            return false
+        elseif lp.status == :Optimal
+            objval = -lp.objval
+            if 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
-            # the constraint is not redundant
-            i = i+1
+            error("LP is not optimal; the status of the LP is $(lp.status)")
         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
@@ -311,79 +311,96 @@ 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
 
 """
     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
+in 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=default_polyhedra_backend(P1, N),
-                      prunefunc=removehredundancy!) where {N<:Real}
+                      backend=nothing,
+                      use_polyhedra_interface=false) 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)])
+    Q = HPOLY([constraints_list(P1); constraints_list(P2)])
 
     # remove redundant constraints
-    if prunefunc == removehredundancy!
-        # convert to polyhedron
+    if use_polyhedra_interface
+        if backend == nothing
+            backend = default_polyhedra_backend(P1, N)
+        end
+        # convert to an hrep, remove the redundancies and convert back to HPOLY
         ph = polyhedron(Q; backend=backend)
-        prunefunc(ph)
+        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