wrap LP status checks in functions

src/Approximations/Approximations.jl

@@ -9,7 +9,8 @@ module Approximations
 using LazySets, LazySets.Arrays, Requires, LinearAlgebra, SparseArrays
 
 using LazySets: _isapprox, _leq, _geq, _rtol, _normal_Vector, isapproxzero,
-                default_lp_solver, _isbounded_stiemke, require, dim, linprog
+                default_lp_solver, _isbounded_stiemke, require, dim, linprog,
+                is_lp_optimal
 
 import LazySets: project
 

src/Approximations/overapproximate.jl

@@ -1548,7 +1548,7 @@ function _overapproximate_zonotope_vrep(X::LazySet{N},
     lower_bounds = vcat(zeros(N, l), fill(N(-Inf), nvariables - l))
     lp = linprog(obj, A, sense, b, lower_bounds, Inf, solver)
 
-    if lp.status != OPTIMAL
+    if !is_lp_optimal(lp.status)
         error("got unexpected status from LP solver: $(lp.status)")
     end
     c = lp.sol[(col_offset_p+1):(col_offset_p+n)]

src/Initialization/init_JuMP.jl

@@ -1,11 +1,33 @@
 using JuMP.MathOptInterface: AbstractOptimizer, OptimizerWithAttributes
-using JuMP: Model, @variable, @objective, @constraint, optimize!, termination_status, objective_value, value 
+using JuMP: Model, @variable, @objective, @constraint, optimize!,
+            termination_status, objective_value, value
 
-# solver statuses
-const OPTIMAL = JuMP.MathOptInterface.OPTIMAL
-const INFEASIBLE = JuMP.MathOptInterface.INFEASIBLE
-const UNBOUNDED = JuMP.MathOptInterface.INFEASIBLE_OR_UNBOUNDED
+# solver status
+function is_lp_optimal(status)
+    return status == JuMP.MathOptInterface.OPTIMAL
+end
+
+function is_lp_infeasible(status; strict::Bool=false)
+    if status == JuMP.MathOptInterface.INFEASIBLE
+        return true
+    end
+    if strict
+        return false
+    end
+    return status == JuMP.MathOptInterface.INFEASIBLE_OR_UNBOUNDED
+end
+
+function is_lp_unbounded(status; strict::Bool=false)
+    if status == JuMP.MathOptInterface.DUAL_INFEASIBLE
+        return true
+    end
+    if strict
+        return false
+    end
+    return status == JuMP.MathOptInterface.INFEASIBLE_OR_UNBOUNDED
+end
 
+# solve a linear program (in the old MathProgBase interface)
 function linprog(c, A, sense::Char, b, l::Number, u::Number, solver)
     n = length(c)
     m = length(b)

src/Interfaces/AbstractPolyhedron_functions.jl

@@ -233,10 +233,10 @@ function remove_redundant_constraints!(constraints::AbstractVector{<:LinearConst
         br = b[non_redundant_indices]
         br[i] = b[j] + one(N)
         lp = linprog(-α, Ar, '<', br, -Inf, Inf, backend)
-        if lp.status == INFEASIBLE
+        if is_lp_infeasible(lp.status)
             # the polyhedron is empty
             return false
-        elseif lp.status == OPTIMAL
+        elseif is_lp_optimal(lp.status)
             objval = -lp.objval
             if _leq(objval, b[j])
                 # the constraint is redundant
@@ -931,9 +931,9 @@ function an_element(P::AbstractPolyhedron{N};
     obj = zeros(N, size(A, 2))
     lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
 
-    if lp.status == OPTIMAL
+    if is_lp_optimal(lp.status)
         return lp.sol
-    elseif lp.status == INFEASIBLE
+    elseif is_lp_infeasible(lp.status)
         error("can't return an element, the polyhedron is empty")
     else
         error("LP returned status $(lp.status) unexpectedly")

src/Interfaces/AbstractZonotope.jl

@@ -330,7 +330,7 @@ function ∈(x::AbstractVector, Z::AbstractZonotope; solver=nothing)
         solver = default_lp_solver(N)
     end
     lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
-    return (lp.status == OPTIMAL) # Infeasible or Unbounded => false
+    return is_lp_optimal(lp.status) # Infeasible or Unbounded => false
 end
 
 """

src/Sets/HPolyhedron.jl

@@ -188,13 +188,15 @@ function σ_helper(d::AbstractVector, P::HPoly, solver)
         l = -Inf
         u = Inf
         lp = linprog(c, A, sense, b, l, u, solver)
-        if lp.status == UNBOUNDED
-            unbounded = true
-        elseif lp.status == INFEASIBLE
+        if is_lp_infeasible(lp.status, strict=true)
             error("the support vector is undefined because the polyhedron is " *
                   "empty")
-        else
+        elseif is_lp_unbounded(lp.status)
+            unbounded = true
+        elseif is_lp_optimal(lp.status)
             unbounded = false
+        else
+            error("got unknown LP status $(lp.status)")
         end
     end
     return (lp, unbounded)
@@ -532,9 +534,9 @@ function isempty(P::HPoly{N},
             solver = default_lp_solver(N)
         end
         lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
-        if lp.status == OPTIMAL
+        if is_lp_optimal(lp.status)
             return witness ? (false, lp.sol) : false
-        elseif lp.status == INFEASIBLE
+        elseif is_lp_infeasible(lp.status)
             return witness ? (true, N[]) : true
         end
         error("LP returned status $(lp.status) unexpectedly")
@@ -687,7 +689,7 @@ function _isbounded_stiemke(P::HPolyhedron{N}; solver=LazySets.default_lp_solver
     At = copy(transpose(A))
     c = ones(N, m)
     lp = linprog(c, At, '=', zeros(n), one(N), Inf, solver)
-    return (lp.status == OPTIMAL)
+    return is_lp_optimal(lp.status)
 end
 
 function is_hyperplanar(P::HPolyhedron)

src/Sets/VPolytope.jl

@@ -220,9 +220,9 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N};
     sense = '='
     obj = zeros(N, m)
     lp = linprog(obj, A, sense, b, lbounds, ubounds, solver)
-    if lp.status == OPTIMAL
+    if is_lp_optimal(lp.status)
         return true
-    elseif lp.status == INFEASIBLE
+    elseif is_lp_infeasible(lp.status)
         return false
     end
     @assert false "LP returned status $(lp.status) unexpectedly"