Revamp optimization deps (#2923)

* remove MathProgBase, use JuMP

* revamp JuMP dep

* use char in sense

* make lbounds a scalar

* cleanup

* swap args

* fix linprog for non-floats

* make linprog available in Approximations

* assign to variable names

* load JuMP earlier

* wrap LP status checks in functions

* fix infeasibility ray

* remove unnecesary import

* fix input for linprog

* skip broken method w/rationals

Co-authored-by: schillic <[email protected]>

Project.toml

@@ -6,12 +6,10 @@ version = "1.54"
 Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 ExprTools = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
 GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
-GLPKMathProgInterface = "3c7084bd-78ad-589a-b5bb-dbd673274bea"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
@@ -27,15 +25,13 @@ Expokit = "0.2"
 ExponentialUtilities = "1"
 ExprTools = "0.1"
 GLPK = "0.11 - 0.15"
-GLPKMathProgInterface = "0.4 - 0.5"
 GR = "0"
 IntervalArithmetic = "0.15 - 0.20"
 IntervalConstraintProgramming = "0.9 - 0.12"
 IntervalMatrices = "0.8"
 Javis = "0.5 - 0.7"
 JuMP = "0.21 - 0.22"
 Makie = "0.9 - 0.16"
-MathProgBase = "0.7"
 MiniQhull = "0.1 - 0.3"
 Optim = "0.15 - 0.22, 1"
 Polyhedra = "0.6"
@@ -70,4 +66,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 
 [targets]
-test = ["CDDLib", "Distributions", "Documenter", "Expokit", "ExponentialUtilities", "GR", "IntervalConstraintProgramming", "IntervalMatrices", "Makie", "Optim", "Polyhedra", "RecipesBase", "StaticArrays", "Symbolics", "TaylorModels", "Test"]
+test = ["CDDLib", "Distributions", "Documenter", "Expokit", "ExponentialUtilities", "GR", "IntervalConstraintProgramming",
+        "IntervalMatrices", "Makie", "Optim", "Polyhedra", "RecipesBase", "StaticArrays", "Symbolics", "TaylorModels", "Test"]

src/Approximations/Approximations.jl

@@ -6,11 +6,11 @@ support vectors.
 """
 module Approximations
 
-using LazySets, LazySets.Arrays, Requires, LinearAlgebra, SparseArrays,
-      MathProgBase
+using LazySets, LazySets.Arrays, Requires, LinearAlgebra, SparseArrays
 
 using LazySets: _isapprox, _leq, _geq, _rtol, _normal_Vector, isapproxzero,
-                default_lp_solver, _isbounded_stiemke, require, dim
+                default_lp_solver, _isbounded_stiemke, require, dim, linprog,
+                is_lp_optimal
 
 import LazySets: project
 

src/Approximations/overapproximate.jl

@@ -1545,10 +1545,11 @@ function _overapproximate_zonotope_vrep(X::LazySet{N},
         end
         col_offset_b += m
     end
-    lower_bounds = vcat(zeros(N, l), fill(N(-Inf), nvariables - l))
-    lp = linprog(obj, A, sense, b, lower_bounds, Inf, solver)
+    lbounds = vcat(zeros(N, l), fill(N(-Inf), nvariables - l))
+    ubounds = fill(Inf, nvariables)
+    lp = linprog(obj, A, sense, b, lbounds, ubounds, 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/ConcreteOperations/intersection.jl

@@ -578,8 +578,6 @@ function intersection(P1::AbstractHPolygon, P2::AbstractHPolygon, prune::Bool=tr
     return P
 end
 
-using MathProgBase.SolverInterface: AbstractMathProgSolver
-
 """
     intersection(P1::AbstractPolyhedron{N},
                  P2::AbstractPolyhedron{N};
@@ -644,7 +642,7 @@ function _intersection_poly(P1::AbstractPolyhedron{N},
     Q = HPOLY([clist_P1; clist_P2])
 
     # remove redundant constraints
-    if backend isa AbstractMathProgSolver
+    if (backend isa AbstractOptimizer) || (backend isa OptimizerWithAttributes)
         # if Q is empty => the feasiblity LP for the list of constraints of Q
         # is infeasible and remove_redundant_constraints! returns `false`
         if !prune || remove_redundant_constraints!(Q, backend=backend)

src/Initialization/init_JuMP.jl

@@ -0,0 +1,57 @@
+using JuMP.MathOptInterface: AbstractOptimizer, OptimizerWithAttributes
+using JuMP: Model, @variable, @objective, @constraint, optimize!,
+            termination_status, objective_value, value, primal_status
+
+# 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
+
+function has_lp_infeasibility_ray(model)
+    return primal_status(model) == JuMP.MathOptInterface.INFEASIBILITY_CERTIFICATE
+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)
+    return linprog(c, A, fill(sense, m), b, fill(l, n), fill(u, n), solver)
+end
+
+function linprog(c, A, sense, b, l, u, solver)
+    n = length(c)
+    model = Model(solver)
+    @variable(model, l[i] <= x[i=1:n] <= u[i])
+    @objective(model, Min, c' * x)
+    eq_rows, ge_rows, le_rows = sense .== '=', sense .== '>', sense .== '<'
+    @constraint(model, A[eq_rows, :] * x .== b[eq_rows])
+    @constraint(model, A[ge_rows, :] * x .>= b[ge_rows])
+    @constraint(model, A[le_rows, :] * x .<= b[le_rows])
+    optimize!(model)
+    return (
+        status = termination_status(model),
+        objval = objective_value(model),
+        sol = value.(x),
+        model = model
+    )
+end

src/Initialization/init_Polyhedra.jl

@@ -3,7 +3,6 @@ eval(quote
     export polyhedron
     using .Polyhedra: HRep, VRep,
                       removehredundancy!, removevredundancy!
-    import JuMP, GLPK
 
     function default_polyhedra_backend(P::LazySet{N}) where {N}
         if LazySets.dim(P) == 1

src/Interfaces/AbstractPolyhedron_functions.jl

@@ -11,12 +11,12 @@ export constrained_dimensions,
 
 # default LP solver for floating-point numbers
 function default_lp_solver(N::Type{<:AbstractFloat})
-    GLPKSolverLP(method=:Simplex)
+    JuMP.optimizer_with_attributes(() -> GLPK.Optimizer(method=GLPK.SIMPLEX))
 end
 
 # default LP solver for rational numbers
 function default_lp_solver(N::Type{<:Rational})
-    GLPKSolverLP(method=:Exact)
+    JuMP.optimizer_with_attributes(() -> GLPK.Optimizer(method=GLPK.EXACT))
 end
 
 # helper function given two possibly different numeric types
@@ -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/LazySets.jl

@@ -3,15 +3,14 @@ __precompile__(true)
 # main module for `LazySets.jl`
 module LazySets
 
-using GLPKMathProgInterface, LinearAlgebra, MathProgBase, Reexport, Requires,
-      SparseArrays
+using LinearAlgebra, Reexport, Requires, SparseArrays
 using LinearAlgebra: checksquare
 import LinearAlgebra: norm, ×, normalize, normalize!
 import SparseArrays: permute
 import Random
 using Random: AbstractRNG, GLOBAL_RNG, SamplerType, shuffle, randperm
 import InteractiveUtils: subtypes
-
+import JuMP, GLPK
 import IntervalArithmetic
 import IntervalArithmetic: radius, ⊂
 
@@ -29,6 +28,11 @@ import .Assertions: activate_assertions, deactivate_assertions
 # activate assertions by default
 activate_assertions(LazySets)
 
+# ==================
+# Linear programming
+# ==================
+include("Initialization/init_JuMP.jl")
+
 # =====================
 # Numeric approximation
 # =====================

src/Sets/HPolyhedron.jl

@@ -153,8 +153,8 @@ function σ(d::AbstractVector{M}, P::HPoly{N};
         # construct the solution from the solver's ray result
         if lp == nothing
             ray = d
-        elseif haskey(lp.attrs, :unboundedray)
-            ray = lp.attrs[:unboundedray]
+        elseif has_lp_infeasibility_ray(lp.model)
+            ray = lp.sol  # infeasibility ray is stored as the solution
         else
             error("LP solver did not return an infeasibility ray")
         end
@@ -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

@@ -215,16 +215,14 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N};
         A[n+1, j] = one(N)
     end
     b = [x; one(N)]
-
-    lbounds = zeros(N, m)
-    ubounds = Inf
-    sense = fill('=', n + 1)
+    lbounds = zero(N)
+    ubounds = N(Inf)
+    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"

test/LazyOperations/Intersection.jl

@@ -24,8 +24,10 @@ for N in [Float64, Rational{Int}, Float32]
     I2 = Singleton(N[1]) ∩ HalfSpace(N[1], N(1))
     @test isbounded(I2) && isboundedtype(typeof(I2))
     I2 = Hyperplane(ones(N, 2), N(1)) ∩ HalfSpace(ones(N, 2), N(1))
-    @test !isbounded(I2) && !isboundedtype(typeof(I2))
-
+    if N != Rational{Int}
+        # Julia crashes if GLPK is used with rationals here, see https://github.com/jump-dev/GLPK.jl/issues/207
+        @test !isbounded(I2) && !isboundedtype(typeof(I2))
+    end
     # is_polyhedral
     @test is_polyhedral(I)
     if N isa AbstractFloat

test/Sets/Polygon.jl

@@ -1,5 +1,3 @@
-using LazySets: _isapprox
-
 for N in [Float64, Float32, Rational{Int}]
     # random polygons
     rand(HPolygon)