From c3435827f5d1e167751adf17e891fc73f238dd35 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Fri, 28 Jan 2022 12:17:14 -0300 Subject: [PATCH 01/15] remove MathProgBase, use JuMP --- Project.toml | 5 ++--- src/ConcreteOperations/intersection.jl | 2 -- src/Initialization/init_JuMP.jl | 18 ++++++++++++++++++ src/Initialization/init_Polyhedra.jl | 1 - src/Interfaces/AbstractPolyhedron_functions.jl | 4 ++-- src/LazySets.jl | 6 +++--- src/init.jl | 1 + 7 files changed, 26 insertions(+), 11 deletions(-) create mode 100644 src/Initialization/init_JuMP.jl diff --git a/Project.toml b/Project.toml index b23dedf9be..1b750d7183 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,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 +34,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 - 1" Polyhedra = "0.6" @@ -70,4 +68,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/ConcreteOperations/intersection.jl b/src/ConcreteOperations/intersection.jl index fe4e906fbc..9a00e50905 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}; diff --git a/src/Initialization/init_JuMP.jl b/src/Initialization/init_JuMP.jl new file mode 100644 index 0000000000..81d5677431 --- /dev/null +++ b/src/Initialization/init_JuMP.jl @@ -0,0 +1,18 @@ +using JuMP: Model, @variable, @objective, @constraint, optimize!, termination_status, objective_value, value + +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) + ) +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..62d1a06694 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) + GLPK.Optimizer(method=GLPK.SIMPLEX) end # default LP solver for rational numbers function default_lp_solver(N::Type{<:Rational}) - GLPKSolverLP(method=:Exact) + GLPK.Optimizer(method=GLPK.EXACT) end # helper function given two possibly different numeric types diff --git a/src/LazySets.jl b/src/LazySets.jl index 60e3c18025..1379d8b4ce 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -3,15 +3,15 @@ __precompile__(true) # main module for `LazySets.jl` module LazySets -using GLPKMathProgInterface, LinearAlgebra, MathProgBase, Reexport, Requires, - SparseArrays +using LinearAlgebra, MathProgBase, 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 +using JuMP.MathOptInterface: AbstractOptimizer import IntervalArithmetic import IntervalArithmetic: radius, ⊂ diff --git a/src/init.jl b/src/init.jl index 5589dcc2a8..32f9d4fa9e 100644 --- a/src/init.jl +++ b/src/init.jl @@ -4,6 +4,7 @@ function __init__() @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" include("Initialization/init_Expokit.jl") @require ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" include("Initialization/init_ExponentialUtilities.jl") @require Javis = "78b212ba-a7f9-42d4-b726-60726080707e" include("Initialization/init_Javis.jl") + include("Initialization/init_JuMP.jl") @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("Initialization/init_Makie.jl") @require MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca" include("Initialization/init_MiniQhull.jl") @require Optim = "429524aa-4258-5aef-a3af-852621145aeb" include("Initialization/init_Optim.jl") From dad80175f5c7cc114827322a9abcde22873a21f7 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Fri, 28 Jan 2022 13:29:16 -0300 Subject: [PATCH 02/15] revamp JuMP dep --- Project.toml | 2 -- src/Approximations/Approximations.jl | 3 +-- src/Approximations/overapproximate.jl | 2 +- src/ConcreteOperations/intersection.jl | 2 +- src/Initialization/init_JuMP.jl | 16 ++++++++++++++-- src/Interfaces/AbstractPolyhedron_functions.jl | 12 ++++++------ src/Interfaces/AbstractZonotope.jl | 2 +- src/LazySets.jl | 3 +-- src/Sets/HPolyhedron.jl | 10 +++++----- src/Sets/VPolytope.jl | 4 ++-- src/init.jl | 3 ++- 11 files changed, 34 insertions(+), 25 deletions(-) diff --git a/Project.toml b/Project.toml index 1b750d7183..7d74cab56a 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" diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index 11d65e63e9..f94833effa 100644 --- a/src/Approximations/Approximations.jl +++ b/src/Approximations/Approximations.jl @@ -6,8 +6,7 @@ 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 diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index acdc54832d..3dd03281da 100644 --- a/src/Approximations/overapproximate.jl +++ b/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 lp.status != OPTIMAL 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 9a00e50905..637910b769 100644 --- a/src/ConcreteOperations/intersection.jl +++ b/src/ConcreteOperations/intersection.jl @@ -642,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 index 81d5677431..c69d869a41 100644 --- a/src/Initialization/init_JuMP.jl +++ b/src/Initialization/init_JuMP.jl @@ -1,9 +1,21 @@ +using JuMP.MathOptInterface: AbstractOptimizer, OptimizerWithAttributes 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 + +function linprog(c, A, sense::Char, b, l::AbstractFloat, u::AbstractFloat, 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) + n = length(c) model = Model(solver) - @variable(model, l[i] <= x[i=1:N] <= u[i]) + @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]) diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl index 62d1a06694..6b39dd0c07 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}) - GLPK.Optimizer(method=GLPK.SIMPLEX) + JuMP.optimizer_with_attributes(() -> GLPK.Optimizer(method=GLPK.SIMPLEX)) end # default LP solver for rational numbers function default_lp_solver(N::Type{<:Rational}) - GLPK.Optimizer(method=GLPK.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 lp.status == INFEASIBLE # the polyhedron is empty return false - elseif lp.status == :Optimal + elseif lp.status == OPTIMAL 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 lp.status == OPTIMAL return lp.sol - elseif lp.status == :Infeasible + elseif lp.status == INFEASIBLE 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..74798c3344 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 (lp.status == OPTIMAL) # Infeasible or Unbounded => false end """ diff --git a/src/LazySets.jl b/src/LazySets.jl index 1379d8b4ce..e823f8f62d 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -3,7 +3,7 @@ __precompile__(true) # main module for `LazySets.jl` module LazySets -using LinearAlgebra, MathProgBase, Reexport, Requires, SparseArrays +using LinearAlgebra, Reexport, Requires, SparseArrays using LinearAlgebra: checksquare import LinearAlgebra: norm, ×, normalize, normalize! import SparseArrays: permute @@ -11,7 +11,6 @@ import Random using Random: AbstractRNG, GLOBAL_RNG, SamplerType, shuffle, randperm import InteractiveUtils: subtypes import JuMP, GLPK -using JuMP.MathOptInterface: AbstractOptimizer import IntervalArithmetic import IntervalArithmetic: radius, ⊂ diff --git a/src/Sets/HPolyhedron.jl b/src/Sets/HPolyhedron.jl index a1860ac933..60648f7b5e 100644 --- a/src/Sets/HPolyhedron.jl +++ b/src/Sets/HPolyhedron.jl @@ -188,9 +188,9 @@ function σ_helper(d::AbstractVector, P::HPoly, solver) l = -Inf u = Inf lp = linprog(c, A, sense, b, l, u, solver) - if lp.status == :Unbounded + if lp.status == UNBOUNDED unbounded = true - elseif lp.status == :Infeasible + elseif lp.status == INFEASIBLE error("the support vector is undefined because the polyhedron is " * "empty") else @@ -532,9 +532,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 lp.status == OPTIMAL return witness ? (false, lp.sol) : false - elseif lp.status == :Infeasible + elseif lp.status == INFEASIBLE return witness ? (true, N[]) : true end error("LP returned status $(lp.status) unexpectedly") @@ -687,7 +687,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 (lp.status == OPTIMAL) end function is_hyperplanar(P::HPolyhedron) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index c0c83b9dce..5639a02284 100644 --- a/src/Sets/VPolytope.jl +++ b/src/Sets/VPolytope.jl @@ -222,9 +222,9 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N}; obj = zeros(N, m) lp = linprog(obj, A, sense, b, lbounds, ubounds, solver) - if lp.status == :Optimal + if lp.status == OPTIMAL return true - elseif lp.status == :Infeasible + elseif lp.status == INFEASIBLE return false end @assert false "LP returned status $(lp.status) unexpectedly" diff --git a/src/init.jl b/src/init.jl index 32f9d4fa9e..778022f9ee 100644 --- a/src/init.jl +++ b/src/init.jl @@ -1,10 +1,11 @@ +include("Initialization/init_JuMP.jl") + function __init__() @require CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60" include("Initialization/init_CDDLib.jl") @require Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" include("Initialization/init_Distributions.jl") @require Expokit = "a1e7a1ef-7a5d-5822-a38c-be74e1bb89f4" include("Initialization/init_Expokit.jl") @require ExponentialUtilities = "d4d017d3-3776-5f7e-afef-a10c40355c18" include("Initialization/init_ExponentialUtilities.jl") @require Javis = "78b212ba-a7f9-42d4-b726-60726080707e" include("Initialization/init_Javis.jl") - include("Initialization/init_JuMP.jl") @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" include("Initialization/init_Makie.jl") @require MiniQhull = "978d7f02-9e05-4691-894f-ae31a51d76ca" include("Initialization/init_MiniQhull.jl") @require Optim = "429524aa-4258-5aef-a3af-852621145aeb" include("Initialization/init_Optim.jl") From 8f6d2b63bfeef228225277a37ef454615d89f8a6 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Fri, 28 Jan 2022 15:20:43 -0300 Subject: [PATCH 03/15] use char in sense --- src/Sets/VPolytope.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index 5639a02284..5c234dd710 100644 --- a/src/Sets/VPolytope.jl +++ b/src/Sets/VPolytope.jl @@ -218,7 +218,7 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N}; lbounds = zeros(N, m) ubounds = Inf - sense = fill('=', n + 1) + sense = '=' obj = zeros(N, m) lp = linprog(obj, A, sense, b, lbounds, ubounds, solver) From d4586e00120c3960d26065208c7b460023a9b4f3 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Fri, 28 Jan 2022 15:21:53 -0300 Subject: [PATCH 04/15] make lbounds a scalar --- src/Sets/VPolytope.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index 5c234dd710..adf62801f9 100644 --- a/src/Sets/VPolytope.jl +++ b/src/Sets/VPolytope.jl @@ -216,7 +216,7 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N}; end b = [x; one(N)] - lbounds = zeros(N, m) + lbounds = zero(N) ubounds = Inf sense = '=' obj = zeros(N, m) From f26c48c7ce68105fc1a01ebc9eb51dba50612d2e Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Fri, 28 Jan 2022 15:22:30 -0300 Subject: [PATCH 05/15] cleanup --- src/Sets/VPolytope.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index adf62801f9..9a8f154308 100644 --- a/src/Sets/VPolytope.jl +++ b/src/Sets/VPolytope.jl @@ -215,13 +215,8 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N}; A[n+1, j] = one(N) end b = [x; one(N)] - - lbounds = zero(N) - ubounds = Inf - sense = '=' obj = zeros(N, m) - - lp = linprog(obj, A, sense, b, lbounds, ubounds, solver) + lp = linprog(obj, A, '=', b, Inf, zero(N), solver) if lp.status == OPTIMAL return true elseif lp.status == INFEASIBLE From d9cd2a27f1ec8acbbba9a580d4e9f2acc2814a9a Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Fri, 28 Jan 2022 15:23:04 -0300 Subject: [PATCH 06/15] swap args --- src/Sets/VPolytope.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index 9a8f154308..4fd2f46bbf 100644 --- a/src/Sets/VPolytope.jl +++ b/src/Sets/VPolytope.jl @@ -216,7 +216,7 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N}; end b = [x; one(N)] obj = zeros(N, m) - lp = linprog(obj, A, '=', b, Inf, zero(N), solver) + lp = linprog(obj, A, '=', b, zero(N), Inf, solver) if lp.status == OPTIMAL return true elseif lp.status == INFEASIBLE From fa7f473c3d10a65556602c8f90cb864efca321b9 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 29 Jan 2022 11:07:11 +0100 Subject: [PATCH 07/15] fix linprog for non-floats --- src/Initialization/init_JuMP.jl | 2 +- src/Sets/VPolytope.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Initialization/init_JuMP.jl b/src/Initialization/init_JuMP.jl index c69d869a41..6c743e1292 100644 --- a/src/Initialization/init_JuMP.jl +++ b/src/Initialization/init_JuMP.jl @@ -6,7 +6,7 @@ const OPTIMAL = JuMP.MathOptInterface.OPTIMAL const INFEASIBLE = JuMP.MathOptInterface.INFEASIBLE const UNBOUNDED = JuMP.MathOptInterface.INFEASIBLE_OR_UNBOUNDED -function linprog(c, A, sense::Char, b, l::AbstractFloat, u::AbstractFloat, solver) +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) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index 4fd2f46bbf..710d6ea644 100644 --- a/src/Sets/VPolytope.jl +++ b/src/Sets/VPolytope.jl @@ -216,7 +216,7 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N}; end b = [x; one(N)] obj = zeros(N, m) - lp = linprog(obj, A, '=', b, zero(N), Inf, solver) + lp = linprog(obj, A, '=', b, zero(N), N(Inf), solver) if lp.status == OPTIMAL return true elseif lp.status == INFEASIBLE From c4759cdc55a23596da533e1ff9e66b8db41ea57e Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 29 Jan 2022 11:10:24 +0100 Subject: [PATCH 08/15] make linprog available in Approximations --- src/Approximations/Approximations.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index f94833effa..71563cd72e 100644 --- a/src/Approximations/Approximations.jl +++ b/src/Approximations/Approximations.jl @@ -9,7 +9,7 @@ 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 + default_lp_solver, _isbounded_stiemke, require, dim, linprog import LazySets: project From 6a45195a89696bd012235ef0f55598fb0864cf90 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 29 Jan 2022 11:11:59 +0100 Subject: [PATCH 09/15] assign to variable names --- src/Sets/VPolytope.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index 710d6ea644..3b67f2139d 100644 --- a/src/Sets/VPolytope.jl +++ b/src/Sets/VPolytope.jl @@ -215,8 +215,11 @@ function ∈(x::AbstractVector{N}, P::VPolytope{N}; A[n+1, j] = one(N) end b = [x; one(N)] + lbounds = zero(N) + ubounds = N(Inf) + sense = '=' obj = zeros(N, m) - lp = linprog(obj, A, '=', b, zero(N), N(Inf), solver) + lp = linprog(obj, A, sense, b, lbounds, ubounds, solver) if lp.status == OPTIMAL return true elseif lp.status == INFEASIBLE From b99a3e97f55613e39b15d6c729aec99288d23680 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 29 Jan 2022 11:32:55 +0100 Subject: [PATCH 10/15] load JuMP earlier --- src/LazySets.jl | 5 +++++ src/init.jl | 2 -- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/LazySets.jl b/src/LazySets.jl index e823f8f62d..9e3aa8c0c8 100644 --- a/src/LazySets.jl +++ b/src/LazySets.jl @@ -28,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/init.jl b/src/init.jl index 778022f9ee..5589dcc2a8 100644 --- a/src/init.jl +++ b/src/init.jl @@ -1,5 +1,3 @@ -include("Initialization/init_JuMP.jl") - function __init__() @require CDDLib = "3391f64e-dcde-5f30-b752-e11513730f60" include("Initialization/init_CDDLib.jl") @require Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" include("Initialization/init_Distributions.jl") From 2b16f29697bd28a5e26bc5110ffcd7983967be21 Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 29 Jan 2022 12:43:47 +0100 Subject: [PATCH 11/15] wrap LP status checks in functions --- src/Approximations/Approximations.jl | 3 +- src/Approximations/overapproximate.jl | 2 +- src/Initialization/init_JuMP.jl | 32 ++++++++++++++++--- .../AbstractPolyhedron_functions.jl | 8 ++--- src/Interfaces/AbstractZonotope.jl | 2 +- src/Sets/HPolyhedron.jl | 16 ++++++---- src/Sets/VPolytope.jl | 4 +-- 7 files changed, 46 insertions(+), 21 deletions(-) diff --git a/src/Approximations/Approximations.jl b/src/Approximations/Approximations.jl index 71563cd72e..cdc911fca3 100644 --- a/src/Approximations/Approximations.jl +++ b/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 diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 3dd03281da..8db27472f7 100644 --- a/src/Approximations/overapproximate.jl +++ b/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)] diff --git a/src/Initialization/init_JuMP.jl b/src/Initialization/init_JuMP.jl index 6c743e1292..7d30f08a0c 100644 --- a/src/Initialization/init_JuMP.jl +++ b/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) diff --git a/src/Interfaces/AbstractPolyhedron_functions.jl b/src/Interfaces/AbstractPolyhedron_functions.jl index 6b39dd0c07..cb30a75636 100644 --- a/src/Interfaces/AbstractPolyhedron_functions.jl +++ 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") diff --git a/src/Interfaces/AbstractZonotope.jl b/src/Interfaces/AbstractZonotope.jl index 74798c3344..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/Sets/HPolyhedron.jl b/src/Sets/HPolyhedron.jl index 60648f7b5e..8b3f217214 100644 --- a/src/Sets/HPolyhedron.jl +++ b/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) diff --git a/src/Sets/VPolytope.jl b/src/Sets/VPolytope.jl index 3b67f2139d..4ca7102f8c 100644 --- a/src/Sets/VPolytope.jl +++ b/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" From 99be8bddc034a5137915b7e88f04b3dafb0eeccd Mon Sep 17 00:00:00 2001 From: schillic Date: Sat, 29 Jan 2022 14:07:40 +0100 Subject: [PATCH 12/15] fix infeasibility ray --- src/Initialization/init_JuMP.jl | 9 +++++++-- src/Sets/HPolyhedron.jl | 4 ++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/Initialization/init_JuMP.jl b/src/Initialization/init_JuMP.jl index 7d30f08a0c..7db49e9148 100644 --- a/src/Initialization/init_JuMP.jl +++ b/src/Initialization/init_JuMP.jl @@ -1,6 +1,6 @@ using JuMP.MathOptInterface: AbstractOptimizer, OptimizerWithAttributes using JuMP: Model, @variable, @objective, @constraint, optimize!, - termination_status, objective_value, value + termination_status, objective_value, value, primal_status # solver status function is_lp_optimal(status) @@ -27,6 +27,10 @@ function is_lp_unbounded(status; strict::Bool=false) 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) @@ -47,6 +51,7 @@ function linprog(c, A, sense, b, l, u, solver) return ( status = termination_status(model), objval = objective_value(model), - sol = value.(x) + sol = value.(x), + model = model ) end diff --git a/src/Sets/HPolyhedron.jl b/src/Sets/HPolyhedron.jl index 8b3f217214..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 From 23366ffb58ba968197ed57895d58cf629845fac0 Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Sat, 29 Jan 2022 17:10:46 -0300 Subject: [PATCH 13/15] remove unnecesary import --- test/Sets/Polygon.jl | 2 -- 1 file changed, 2 deletions(-) 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) From 8132f11553b99f871efd1ff9e66fb3e7339a2daf Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Sat, 29 Jan 2022 17:11:05 -0300 Subject: [PATCH 14/15] fix input for linprog --- src/Approximations/overapproximate.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Approximations/overapproximate.jl b/src/Approximations/overapproximate.jl index 8db27472f7..462ed06770 100644 --- a/src/Approximations/overapproximate.jl +++ b/src/Approximations/overapproximate.jl @@ -1545,8 +1545,9 @@ 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 !is_lp_optimal(lp.status) error("got unexpected status from LP solver: $(lp.status)") From 9e311db04040276ee615e250599a452acc7684bc Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Sat, 29 Jan 2022 17:11:28 -0300 Subject: [PATCH 15/15] skip broken method w/rationals --- test/LazyOperations/Intersection.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) 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