From 8f4a611082d3fec490fe0b19642620f6263070c0 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Sun, 30 Jan 2022 07:55:54 -0300 Subject: [PATCH] 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 --- Project.toml | 7 +-- src/Approximations/Approximations.jl | 6 +- src/Approximations/overapproximate.jl | 7 ++- src/ConcreteOperations/intersection.jl | 4 +- src/Initialization/init_JuMP.jl | 57 +++++++++++++++++++ src/Initialization/init_Polyhedra.jl | 1 - .../AbstractPolyhedron_functions.jl | 12 ++-- src/Interfaces/AbstractZonotope.jl | 2 +- src/LazySets.jl | 10 +++- src/Sets/HPolyhedron.jl | 20 ++++--- src/Sets/VPolytope.jl | 12 ++-- test/LazyOperations/Intersection.jl | 6 +- test/Sets/Polygon.jl | 2 - 13 files changed, 101 insertions(+), 45 deletions(-) create mode 100644 src/Initialization/init_JuMP.jl diff --git a/Project.toml b/Project.toml index b5aa44e9ca..d3f6268d27 100644 --- a/Project.toml +++ b/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,7 +25,6 @@ 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" @@ -35,7 +32,6 @@ 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"] diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index 11d65e63e9..cdc911fca3 100644 --- a/src/Approximations/Approximations.jl +++ b/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 diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index acdc54832d..462ed06770 100644 --- a/src/Approximations/overapproximate.jl +++ b/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)] diff --git a/src/ConcreteOperations/intersection.jl b/src/ConcreteOperations/intersection.jl index fe4e906fbc..637910b769 100644 --- a/src/ConcreteOperations/intersection.jl +++ b/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) diff --git a/src/Initialization/init_JuMP.jl b/src/Initialization/init_JuMP.jl new file mode 100644 index 0000000000..7db49e9148 --- /dev/null +++ b/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 diff --git a/src/Initialization/init_Polyhedra.jl b/src/Initialization/init_Polyhedra.jl index aa52379655..70fcb00e14 100644 --- a/src/Initialization/init_Polyhedra.jl +++ b/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 diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl index d691712ad3..cb30a75636 100644 --- a/src/Interfaces/AbstractPolyhedron_functions.jl +++ b/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") diff --git a/src/Interfaces/AbstractZonotope.jl b/src/Interfaces/AbstractZonotope.jl index a3bc3ea0ec..09dec543ad 100644 --- a/src/Interfaces/AbstractZonotope.jl +++ b/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 """ diff --git a/src/LazySets.jl b/src/LazySets.jl index 60e3c18025..9e3aa8c0c8 100644 --- a/src/LazySets.jl +++ b/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 # ===================== diff --git a/src/Sets/HPolyhedron.jl b/src/Sets/HPolyhedron.jl index a1860ac933..fa6bd67093 100644 --- a/src/Sets/HPolyhedron.jl +++ b/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) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index c0c83b9dce..4ca7102f8c 100644 --- a/src/Sets/VPolytope.jl +++ b/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" diff --git a/test/LazyOperations/Intersection.jl b/test/LazyOperations/Intersection.jl index 03460f0f8a..e5edae3a50 100644 --- a/test/LazyOperations/Intersection.jl +++ b/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 diff --git a/test/Sets/Polygon.jl b/test/Sets/Polygon.jl index 1870279be5..186755d86e 100644 --- a/test/Sets/Polygon.jl +++ b/test/Sets/Polygon.jl @@ -1,5 +1,3 @@ -using LazySets: _isapprox - for N in [Float64, Float32, Rational{Int}] # random polygons rand(HPolygon)