Merge pull request #157 from JuliaPolyhedra/bl/jump0.19

Update to JuMP v0.19

Author: blegat

.travis.yml

@@ -4,7 +4,6 @@ os:
   - linux
   - osx
 julia:
-  - 0.7
   - 1.0
   - 1.1
 #addons:

REQUIRE

@@ -1,6 +1,6 @@
-julia 0.7
+julia 1.0
 StaticArrays 0.5
 MathProgBase 0.5.10 0.8
-JuMP 0.16 0.19
+JuMP 0.19
 GeometryTypes 0.4
 RecipesBase 0.2

appveyor.yml

@@ -1,6 +1,5 @@
 environment:
   matrix:
-  - julia_version: 0.7
   - julia_version: 1.0
 
 platform:

docs/src/optimization.md

@@ -10,7 +10,7 @@ Polyhedra.MathProgBase.linprog
 
 If the V-representation of the polyhedron has been computed, it can be used to solve the linear program.
 ```@docs
-VRepSolver
+VRepOptimizer
 ```
 
 Otherwise, any programming solver implementing the [MathProgBase](https://github.com/JuliaOpt/MathProgBase.jl) interface can be used. See [here](http://www.juliaopt.org/) for a list of available solvers.
@@ -46,5 +46,5 @@ In fact, the MathProgBase representation of the feasible set of a linear program
 \end{align*}
 ```
 
-has `LPHRepresentation` as a corresponding H-representation.
-A JuMP model can be converted to this representation using `LPHRepresentation(m)`.
+has `LPHRep` as a corresponding H-representation.
+A JuMP model can be converted to this representation using `LPHRep(m)`.

src/Polyhedra.jl

@@ -1,5 +1,3 @@
-__precompile__()
-
 module Polyhedra
 
 using SparseArrays
@@ -10,9 +8,9 @@ export Polyhedron
 abstract type Library end
 abstract type Polyhedron{T} end
 
-import MathProgBase
-const MPB = MathProgBase
-const MPBSI = MPB.SolverInterface
+import JuMP
+const Solver = JuMP.OptimizerFactory
+const SolverOrNot = Union{Nothing, Solver}
 
 coefficient_type(::Union{AbstractVector{T}, Type{<:AbstractVector{T}}}) where T = T
 similar_type(::Type{<:Vector}, ::Int, ::Type{T}) where T = Vector{T}
@@ -82,9 +80,8 @@ include("fulldim.jl")
 
 # Optimization
 include("opt.jl")
-include("lpqp_to_polyhedra.jl")
-include("polyhedra_to_lpqp.jl")
-include("vrepsolver.jl")
+include("polyhedra_to_lp_bridge.jl")
+include("vrep_optimizer.jl")
 include("default.jl")
 
 # Visualization

src/center.jl

@@ -1,55 +1,58 @@
 export hchebyshevcenter, vchebyshevcenter, chebyshevcenter
 using JuMP
+
 """
     hchebyshevcenter(p::HRep[, solver])
 
 Return a tuple with the center and radius of the largest euclidean ball contained in the polyhedron `p`.
 Throws an error if the polyhedron is empty or if the radius is infinite.
 """
-function hchebyshevcenter(p::HRep, solver=Polyhedra.solver(p))
-    m = Model(solver=solver)
-    c = @variable m [1:fulldim(p)]
+function hchebyshevcenter(p::HRep, solver=Polyhedra.default_solver(p))
+    model = JuMP.Model(solver)
+    c = JuMP.@variable(model, [1:fulldim(p)])
     for hp in hyperplanes(p)
-        a = Vector{Float64}(hp.a)
-        β = Float64(hp.β)
-        @constraint m dot(a, c) == β
+        a = convert(Vector{Float64}, hp.a)
+        β = convert(Float64, hp.β)
+        JuMP.@constraint(model, dot(a, c) == β)
     end
-    @variable m r[1:nhalfspaces(p)] >= 0
+    JuMP.@variable(model, r[1:nhalfspaces(p)] >= 0)
     for (i, hs) in enumerate(halfspaces(p))
-        a = Vector{Float64}(hs.a)
-        β = Float64(hs.β)
-        @constraint m dot(a, c) + r[i] * norm(a, 2) <= β
+        a = convert(Vector{Float64}, hs.a)
+        β = convert(Float64, hs.β)
+        JuMP.@constraint(model, dot(a, c) + r[i] * norm(a, 2) <= β)
     end
-    @variable m minr >= 0
-    @constraint m minr .<= r
-    @objective m Max minr
-    status = solve(m)
-    if status != :Optimal
-        if status == :Infeasible
-            error("An empty polyhedron has no H-Chebyshev center")
-        elseif status == :Unbounded
-            error("The polyhedron contains euclidean ball of arbitrary large radius")
+    JuMP.@variable(model, minr >= 0)
+    JuMP.@constraint(model, minr .<= r)
+    JuMP.@objective(model, Max, minr)
+    JuMP.optimize!(model)
+    term = JuMP.termination_status(model)
+    if term != MOI.OPTIMAL
+        if term == MOI.INFEASIBLE
+            error("An empty polyhedron has no H-Chebyshev center.")
+        elseif term == MOI.DUAL_INFEASIBLE
+            error("The polyhedron contains euclidean ball of arbitrary large radius.")
         else
-            error("Solver returned $status when computing the H-Chebyshev center")
+            error("Solver returned $term when computing the H-Chebyshev center.")
         end
     end
-    @constraint m minr == getvalue(minr)
-    @variable m maxr >= 0
-    @constraint m maxr .>= r
-    @objective m Min maxr
-    status = solve(m)
-    @assert status == :Optimal
-    (getvalue(c), getvalue(minr))
+    JuMP.@constraint(model, minr == JuMP.value(minr))
+    JuMP.@variable(model, maxr >= 0)
+    JuMP.@constraint(model, maxr .>= r)
+    JuMP.@objective(model, Min, maxr)
+    JuMP.optimize!(model)
+    term = JuMP.termination_status(model)
+    @assert term == MOI.OPTIMAL
+    return (JuMP.value.(c), JuMP.value(minr))
 end
 
-# TODO solver here should not be VRepSolver
+# TODO solver here should not be VRepOptimizer
 """
     vchebyshevcenter(p::VRep[, solver])
 
 Return a tuple with the center and radius of the smallest euclidean ball containing the polyhedron `p`.
 Throws an error if the polyhedron is empty or if the radius is infinite (i.e. `p` is not a polytope, it contains rays).
 """
-function vchebyshevcenter(p::VRep, solver=Polyhedra.solver(p))
+function vchebyshevcenter(p::VRep, solver=Polyhedra.default_solver(p))
     error("TODO")
 end
 
@@ -58,12 +61,12 @@ end
 
 If `p` is a H-representation or is a polyhedron for which the H-representation has already been computed, calls `hchebyshevcenter`, otherwise, call `vchebyshevcenter`.
 """
-function chebyshevcenter(p::Polyhedron, solver=Polyhedra.solver(p))
+function chebyshevcenter(p::Polyhedron, solver=Polyhedra.default_solver(p))
     if hrepiscomputed(p)
         hchebyshevcenter(p, solver)
     else
         vchebyshevcenter(p, solver)
     end
 end
-chebyshevcenter(p::HRepresentation, solver=Polyhedra.solver(p)) = hchebyshevcenter(p, solver)
-chebyshevcenter(p::VRepresentation, solver=Polyhedra.solver(p)) = vchebyshevcenter(p, solver)
+chebyshevcenter(p::HRepresentation, solver=Polyhedra.default_solver(p)) = hchebyshevcenter(p, solver)
+chebyshevcenter(p::VRepresentation, solver=Polyhedra.default_solver(p)) = vchebyshevcenter(p, solver)

src/default.jl

@@ -13,8 +13,8 @@ end
 function default_type(d::StaticArrays.Size{(1,)}, ::Type{T}) where T
     return Interval{T, StaticArrays.SVector{1, T}, typeof(d)}
 end
-function default_type(d::Int, ::Type{T}) where T
-    if d == 1
+function default_type(d::Integer, ::Type{T}) where T
+    if isone(d)
         return Interval{T, Vector{T}, typeof(d)}
     else
         return DefaultPolyhedron{T, Intersection{T, Vector{T}, typeof(d)}, Hull{T, Vector{T}, typeof(d)}}
@@ -42,8 +42,8 @@ _default_type(::Type{T}) where T = T
 _default_type(::Type{AbstractFloat}) = Float64
 default_library(::StaticArrays.Size, T::Type) = DefaultLibrary{_default_type(T)}()
 default_library(::StaticArrays.Size{(1,)}, T::Type) = IntervalLibrary{_default_type(T)}()
-function default_library(d::Int, ::Type{T}) where T
-    if d == 1
+function default_library(d::Integer, ::Type{T}) where T
+    if isone(d)
         return IntervalLibrary{_default_type(T)}()
     else
         return DefaultLibrary{_default_type(T)}()
@@ -78,35 +78,48 @@ default_solver(p::Rep) = nothing
 function default_solver(p::Rep, ps::Rep...)
     s = default_solver(p)
     if s === nothing
-        default_solver(ps...)
+        return default_solver(ps...)
     else
-        s
+        return s
     end
 end
 
 """
-    solver(p::Rep, solver::MathProgBase.AbstractMathProgSolver=default_solver(p))
-
-If the V-representation of `p` has been computed, returns `VRepSolver()`, otherwise, returns `solver`.
+    linear_objective_solver(p::Rep, solver::Union{Nothing, JuMP.OptimizerFactory}=default_solver(p))
+
+Return the solver to use for optimizing a linear objective over the polyhedron `p`, i.e.
+```julia
+model = Model(solver)
+x = @variable(model, [1:fulldim(p)])
+@constraint(model, x in p)
+@objective(model, c ⋅ x)
+```
+for some vector `c`.
+
+By default, if the V-representation of `p` has been computed, it returns
+`VRepOptimizer()`, otherwise, it returns `solver`.
+
+If the problem has constraints different to `x in p`, use `default_solver(p)` instead
+as the fact that the V-representation of `p` has been computed does not help.
 """
-function solver(p::Rep, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(p))
+function linear_objective_solver(p::Rep{T}, solver::SolverOrNot=default_solver(p)) where T
     if vrepiscomputed(p)
-        VRepSolver()
+        return with_optimizer(VRepOptimizer{T})
     else
-        solver
+        return solver
     end
 end
 
-solver(v::VRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(v)) = VRepSolver()
-solver(h::HRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}=default_solver(h)) = solver
+linear_objective_solver(v::VRepresentation{T}, solver::SolverOrNot=default_solver(v)) where {T} = with_optimizer(VRepOptimizer{T})
+linear_objective_solver(h::HRepresentation, solver::SolverOrNot=default_solver(h)) = solver
 
 _promote_reptype(P::Type{<:HRep}, ::Type{<:HRep}) = P
 _promote_reptype(P::Type{<:VRep}, ::Type{<:VRep}) = P
 # Breaks ambiguity with above two methods
 _promote_reptype(P::Type{<:Polyhedron}, ::Type{<:Polyhedron}) = P
 
 function promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep})
-    _promote_reptype(P1, P2)
+    return _promote_reptype(P1, P2)
 end
 promote_reptype(P1::Type{<:Rep}, P2::Type{<:Rep}, P::Type{<:Rep}...) = promote_reptype(promote_reptype(P1, P2), P...)
 promote_reptype(P::Type{<:Rep}) = P
@@ -121,13 +134,13 @@ function constructpolyhedron(RepT::Type{<:Rep{T}}, d::FullDim, p::Tuple{Vararg{R
             return RepT(d, it..., solver=solver)
         end
     end
-    RepT(d, it...)::RepT # FIXME without this type annotation even convexhull(::PointsHull{2,Int64,Array{Int64,1}}, ::PointsHull{2,Int64,Array{Int64,1}}) is not type stable, why ?
+    return RepT(d, it...)::RepT # FIXME without this type annotation even convexhull(::PointsHull{2,Int64,Array{Int64,1}}, ::PointsHull{2,Int64,Array{Int64,1}}) is not type stable, why ?
 end
 
 function default_similar(p::Tuple{Vararg{Rep}}, d::FullDim, ::Type{T}, it::It{T}...) where T
     # Some types in p may not support `d` or `T` so we call `similar_type` after `promote_reptype`
     RepT = similar_type(promote_reptype(typeof.(p)...), d, T)
-    constructpolyhedron(RepT, d, p, it...)
+    return constructpolyhedron(RepT, d, p, it...)
 end
 
 """
@@ -138,13 +151,13 @@ The type of the result will be chosen closer to the type of `p[1]`.
 """
 Base.similar(p::Tuple{Vararg{Rep}}, d::FullDim, ::Type{T}, it::It{T}...) where {T} = default_similar(p, d, T, it...)
 function promote_coefficient_type(p::Tuple{Vararg{Rep}})
-    promote_type(coefficient_type.(p)...)
+    return promote_type(coefficient_type.(p)...)
 end
 function Base.similar(p::Tuple{Vararg{Rep}}, d::FullDim, it::It...)
     T = promote_coefficient_type(p)
-    similar(p, d, T, it...)
+    return similar(p, d, T, it...)
 end
 function Base.similar(p::Tuple{Vararg{Rep}}, it::It...)
-    similar(p, FullDim(p[1]), it...)
+    return similar(p, FullDim(p[1]), it...)
 end
 Base.similar(p::Rep, args...) = similar((p,), args...)

src/defaultlibrary.jl

@@ -7,7 +7,7 @@ Default library for polyhedra of dimension larger than 1 ([`IntervalLibrary`](@r
 The library implements the bare minimum and uses the fallback implementation for all operations.
 """
 struct DefaultLibrary{T} <: Library
-    solver::Union{Nothing, MPB.AbstractMathProgSolver}
+    solver::SolverOrNot
     function DefaultLibrary{T}(solver=nothing) where T
         new{T}(solver)
     end
@@ -18,21 +18,21 @@ similar_library(lib::DefaultLibrary, d::FullDim, ::Type{T}) where T = default_li
 mutable struct DefaultPolyhedron{T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}} <: Polyhedron{T}
     hrep::Union{HRepT, Nothing}
     vrep::Union{VRepT, Nothing}
-    solver::Union{Nothing, MPB.AbstractMathProgSolver}
-    function DefaultPolyhedron{T, HRepT, VRepT}(hrep::Union{HRepT, Nothing}, vrep::Union{VRepT, Nothing}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
+    solver::SolverOrNot
+    function DefaultPolyhedron{T, HRepT, VRepT}(hrep::Union{HRepT, Nothing}, vrep::Union{VRepT, Nothing}, solver::SolverOrNot) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
         new{T, HRepT, VRepT}(hrep, vrep, solver)
     end
 end
-function DefaultPolyhedron{T, HRepT, VRepT}(hrep::HRepT, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
+function DefaultPolyhedron{T, HRepT, VRepT}(hrep::HRepT, solver::SolverOrNot) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
     DefaultPolyhedron{T, HRepT, VRepT}(hrep, nothing, solver)
 end
-function DefaultPolyhedron{T, HRepT, VRepT}(vrep::VRepT, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
+function DefaultPolyhedron{T, HRepT, VRepT}(vrep::VRepT, solver::SolverOrNot) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
     DefaultPolyhedron{T, HRepT, VRepT}(nothing, vrep, solver)
 end
-function DefaultPolyhedron{T, HRepT, VRepT}(hrep::HRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
+function DefaultPolyhedron{T, HRepT, VRepT}(hrep::HRepresentation, solver::SolverOrNot) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
     DefaultPolyhedron{T, HRepT, VRepT}(convert(HRepT, hrep), solver)
 end
-function DefaultPolyhedron{T, HRepT, VRepT}(vrep::VRepresentation, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
+function DefaultPolyhedron{T, HRepT, VRepT}(vrep::VRepresentation, solver::SolverOrNot) where {T, HRepT<:HRepresentation{T}, VRepT<:VRepresentation{T}}
     DefaultPolyhedron{T, HRepT, VRepT}(convert(VRepT, vrep), solver)
 end
 
@@ -54,13 +54,13 @@ function DefaultPolyhedron{T, HRepT, VRepT}(d::FullDim, vits::VIt...; solver=not
 end
 
 # Need fulltype in case the use does `intersect!` with another element
-DefaultPolyhedron{T}(rep::Representation, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T} = DefaultPolyhedron{T}(change_coefficient_type(rep, T), solver)
-function DefaultPolyhedron{T}(rep::HRepresentation{T}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T}
+DefaultPolyhedron{T}(rep::Representation, solver::SolverOrNot) where {T} = DefaultPolyhedron{T}(change_coefficient_type(rep, T), solver)
+function DefaultPolyhedron{T}(rep::HRepresentation{T}, solver::SolverOrNot) where {T}
     HRepT = fulltype(typeof(rep))
     VRepT = dualtype(HRepT)
     DefaultPolyhedron{T, HRepT, VRepT}(rep, solver)
 end
-function DefaultPolyhedron{T}(rep::VRepresentation{T}, solver::Union{Nothing, MPB.AbstractMathProgSolver}) where {T}
+function DefaultPolyhedron{T}(rep::VRepresentation{T}, solver::SolverOrNot) where {T}
     VRepT = fulltype(typeof(rep))
     HRepT = dualtype(VRepT)
     DefaultPolyhedron{T, HRepT, VRepT}(rep, solver)

src/dimension.jl

@@ -45,7 +45,7 @@ fulldim(p) = fulldim(FullDim(p))
 
 function fulldim end
 
-fulldim(N::Int) = N
+fulldim(N::Integer) = convert(Int, N)
 fulldim(::StaticArrays.Size{N}) where N = N[1]
 
 neg_fulldim(N::Int) = -N

src/jump.jl

@@ -6,24 +6,11 @@ using JuMP
 Builds an H-representation from the feasibility set of the JuMP model `model`.
 Note that if non-linear constraint are present in the model, they are ignored.
 """
-hrep(model::JuMP.Model) = LPHRepresentation(model)
+hrep(model::JuMP.Model) = LPHRep(model)
 
-function LPHRepresentation(model::JuMP.Model)
-    # Inspired from Joey Huchette's code in ConvexHull.jl
-    A = JuMP.prepConstrMatrix(model)
-    c = JuMP.prepAffObjective(model)
-    lb, ub = JuMP.prepConstrBounds(model)
-    l, u = model.colLower, model.colUpper
+LPHRep(model::JuMP.Model) = LPHRep(backend(model))
 
-    m, n = size(A)
-    @assert m == length(lb) == length(ub)
-    @assert model.nlpdata == nothing
-    @assert isempty(model.quadconstr)
-    @assert isempty(model.sosconstr)
-
-    LPHRepresentation(A, l, u, lb, ub)
-end
-
-function polyhedron(model::JuMP.Model, lib::Library=default_library(FullDim{model.numCols}(), Float64))
-    polyhedron(LPHRepresentation(model), lib)
+function polyhedron(model::JuMP.Model,
+                    lib::Library=default_library(FullDim(JuMP.num_variables(model)), Float64))
+    polyhedron(LPHRep(model), lib)
 end