Add MathProgBase solver interface.

REQUIRE

@@ -1,3 +1,4 @@
 BinDeps
 julia 0.6
 Compat
+MathProgBase 0.7 0.8

src/OSQP.jl

@@ -2,6 +2,8 @@ __precompile__()
 
 module OSQP
 
+export OSQPMathProgBaseInterface
+
 # Compatibility stuff
 using Compat
 using Compat.SparseArrays
@@ -34,5 +36,6 @@ end
 include("constants.jl")
 include("types.jl")
 include("interface.jl")
+include("mpbinterface.jl")
 
 end # module

src/interface.jl

@@ -467,7 +467,7 @@ end
 
 
 
-function warm_start!(model::OSQP.Model; x::Vector{Float64} = nothing, y::Vector{Float64} = nothing)
+function warm_start!(model::OSQP.Model; x::Union{Vector{Float64}, Nothing} = nothing, y::Union{Vector{Float64}, Nothing} = nothing)
     # Get problem dimensions
     (n, m) = OSQP.dimensions(model)
 

src/mpbinterface.jl

@@ -0,0 +1,197 @@
+module OSQPMathProgBaseInterface
+
+import MathProgBase
+import OSQP
+
+import MathProgBase: numvar, numconstr
+
+struct OSQPSolver <: MathProgBase.AbstractMathProgSolver
+    settings::Dict{Symbol,Any}
+end
+
+OSQPSolver(; kwargs...) = OSQPSolver(Dict{Symbol, Any}(k => v for (k, v) in kwargs))
+
+mutable struct OSQPModel <: MathProgBase.AbstractLinearQuadraticModel
+    settings::Dict{Symbol, Any}
+    inner::OSQP.Model
+
+    P::SparseMatrixCSC{Float64, Int}
+    q::Vector{Float64}
+    A::SparseMatrixCSC{Float64, Int}
+    l::Vector{Float64}
+    u::Vector{Float64}
+
+    xwarmstart::Vector{Float64}
+    dowarmstart::Bool
+
+    results::OSQP.Results
+
+    function OSQPModel(settings::Associative{Symbol, Any})
+        P = spzeros(Float64, 0, 0)
+        q = Float64[]
+        A = spzeros(Float64, 0, 0)
+        l = Float64[]
+        u = Float64[]
+        xwarmstart = Float64[]
+        dowarmstart = false
+        new(copy(settings), OSQP.Model(), P, q, A, l, u, xwarmstart, dowarmstart) # leave results undefined
+    end
+end
+
+function Base.resize!(model::OSQPModel, n, m)
+    nchange = n != numvar(model)
+    mchange = m != numconstr(model)
+    if nchange
+        model.P = spzeros(Float64, n, n)
+        resize!(model.q, n)
+        resize!(model.xwarmstart, n)
+        model.dowarmstart = false
+    end
+    if mchange
+        resize!(model.l, m)
+        resize!(model.u, m)
+    end
+    if nchange || mchange
+        model.A = spzeros(Float64, m, n)
+    end
+
+    model
+end
+
+checksolved(model::OSQPModel) = isdefined(model, :results) || error("Model has not been solved.")
+variablebounderror() = error("Variable bounds are not supported (use constraints instead).")
+senseerror() = error("Only objective sense :Min is currently supported to avoid confusing behavior.")
+
+# http://mathprogbasejl.readthedocs.io/en/latest/solverinterface.html
+MathProgBase.LinearQuadraticModel(solver::OSQPSolver) = OSQPModel(solver.settings)
+MathProgBase.getsolution(model::OSQPModel) = (checksolved(model); model.results.x)
+MathProgBase.getobjval(model::OSQPModel) = (checksolved(model); model.results.info.obj_val)
+
+function MathProgBase.optimize!(model::OSQPModel)
+    OSQP.setup!(model.inner; P = model.P, q = model.q, A = model.A, l = model.l, u = model.u, model.settings...)
+    model.dowarmstart && OSQP.warm_start!(model.inner, x = model.xwarmstart)
+    model.results = OSQP.solve!(model.inner)
+end
+
+function MathProgBase.status(model::OSQPModel)::Symbol
+    checksolved(model)
+    osqpstatus = model.results.info.status
+    status = osqpstatus # if OSQP status can't be mapped to a standard return value, just return as is
+
+    # map to standard return status values:
+    osqpstatus == :Solved && (status = :Optimal)
+    osqpstatus == :Max_iter_reached && (status = :UserLimit)
+    osqpstatus == :Interrupted && (status = :UserLimit) # following Gurobi.jl
+    osqpstatus == :Primal_infeasible && (status = :Infeasible)
+    osqpstatus == :Primal_infeasible_inaccurate && (status = :Infeasible)
+    osqpstatus == :Dual_infeasible && (status = :DualityFailure)
+
+    return status
+end
+
+# TODO: getobjbound
+# TODO: getobjgap
+MathProgBase.getrawsolver(model::OSQPModel) = model.inner
+MathProgBase.getsolvetime(model::OSQPModel) = (checksolved(model); model.results.info.run_time)
+MathProgBase.setsense!(model::OSQPModel, sense::Symbol) = sense == :Min || senseerror()
+MathProgBase.getsense(model::OSQPModel) = :Min
+MathProgBase.numvar(model::OSQPModel) = size(model.A, 2)
+MathProgBase.numconstr(model::OSQPModel) = size(model.A, 1)
+MathProgBase.freemodel!(model::OSQPModel) = OSQP.clean!(model.inner)
+
+function Base.copy(model::OSQPModel)
+    ret = OSQPModel(model.settings)
+    resize!(ret, numvar(model), numconstr(model))
+    copy!(ret.P, model.P)
+    copy!(ret.q, model.q)
+    copy!(ret.A, model.A)
+    copy!(ret.l, model.l)
+    copy!(ret.u, model.u)
+    ret
+end
+
+MathProgBase.setvartype!(model::OSQPModel, v::Vector{Symbol}) = any(x -> x != :Cont, v) && error("OSQP only supports continuous variables.")
+MathProgBase.getvartype(model::OSQPModel) = fill(:Cont, numvar(model))
+
+function MathProgBase.setparameters!(x::Union{OSQPSolver, OSQPModel}; Silent = nothing)
+    if Silent != nothing
+        Silent::Bool
+        x.settings[:verbose] = !Silent
+    end
+    x
+end
+
+MathProgBase.setwarmstart!(model::OSQPModel, v) = (copy!(model.xwarmstart, v); model.dowarmstart = true)
+
+
+# http://mathprogbasejl.readthedocs.io/en/latest/lpqcqp.html#linearquadratic-models
+# TODO: loadproblem!(m::AbstractLinearQuadraticModel, filename::String)
+
+function MathProgBase.loadproblem!(model::OSQPModel, A, l, u, c, lb, ub, sense)
+    (any(x -> x != -Inf, l) || any(x -> x != Inf, u)) && variablebounderror()
+    sense == :Min || senseerror()
+    m, n = size(A)
+    resize!(model, n, m)
+    copy!(model.q, c)
+    copy!(model.l, lb)
+    copy!(model.u, ub)
+    copy!(model.A, A)
+    model
+end
+
+# TODO: writeproblem
+MathProgBase.getvarLB(model::OSQPModel) = fill(-Inf, numvar(model))
+MathProgBase.setvarLB!(model::OSQPModel, l) = variablebounderror()
+MathProgBase.getvarUB(model::OSQPModel) = fill(Inf, numvar(model))
+MathProgBase.setvarUB!(model::OSQPModel, l) = variablebounderror()
+MathProgBase.getconstrLB(model::OSQPModel) = model.l
+MathProgBase.setconstrLB!(model::OSQPModel, lb) = (copy!(model.l, lb); model)
+MathProgBase.getconstrUB(model::OSQPModel) = model.u
+MathProgBase.setconstrUB!(model::OSQPModel, ub) = (copy!(model.u, ub); model)
+MathProgBase.setobj!(model::OSQPModel, c) = (copy!(model.q, c); model)
+MathProgBase.getconstrmatrix(model::OSQPModel) = model.A
+# TODO: addvar!
+# TODO: delvars!
+# TODO: addconstr!
+# TODO: delconstrs!
+# TODO: changecoeffs!
+MathProgBase.numlinconstr(model::OSQPModel) = numconstr(model)
+MathProgBase.getconstrsolution(model::OSQPModel) = model.A * getsolution(model)
+# TODO: getreducedcosts
+MathProgBase.getconstrduals(model::OSQPModel) = (checksolved(model); model.results.y)
+# TODO: getinfeasibilityray
+# TODO: getbasis
+# TODO: getunboundedray
+
+
+# http://mathprogbasejl.readthedocs.io/en/latest/lpqcqp.html#quadratic-programming
+MathProgBase.numquadconstr(model::OSQPModel) = 0
+MathProgBase.setquadobj!(model::OSQPModel, Q::Matrix) = (copy!(model.P, Q); model)
+
+function MathProgBase.setquadobj!(model::OSQPModel, rowidx::Vector, colidx::Vector, quadval::Vector)
+    nterms = length(quadval)
+    @boundscheck length(rowidx) == nterms || error()
+    @boundscheck length(colidx) == nterms || error()
+
+    # zero out coeffs
+    for i = 1 : nterms
+        @inbounds row = rowidx[i]
+        @inbounds col = colidx[i]
+        model.P[row, col] = 0
+    end
+
+    # add new coeffs
+    for i = 1 : nterms
+        @inbounds row = rowidx[i]
+        @inbounds col = colidx[i]
+        @inbounds val = quadval[i]
+        model.P[row, col] += val
+        model.P[col, row] = model.P[row, col]
+    end
+
+    model
+end
+
+# Note: skipping quadconstr methods
+
+end # module

test/mpbinterface.jl

@@ -0,0 +1,115 @@
+using MathProgBase
+import MathProgBase: LinearQuadraticModel, loadproblem!, setquadobj!, optimize!, status, numvar, numconstr, setwarmstart!, getsense,
+    getobjval, getsolution, setsense!, getvarLB, setvarLB!, getvarUB, setvarUB!, getconstrLB, setconstrLB!, getconstrUB, setconstrUB!,
+    getconstrmatrix, numlinconstr, getsolvetime, getrawsolver, getvartype, setvartype!, getconstrduals
+
+@testset "MathProgBase" begin
+    solver = OSQPMathProgBaseInterface.OSQPSolver(eps_abs = 1e-7, eps_rel = 1e-16)
+    MathProgBase.setparameters!(solver, Silent=true)
+
+    @testset "quadprog" begin
+        # modified from joinpath(Pkg.dir("MathProgBase"), "test", "quadprog.jl"):
+        sol = quadprog([0., 0., 0.],[2. 1. 0.; 1. 2. 1.; 0. 1. 2.],[1. 2. 3.; 1. 1. 0.],'>',[4., 1.],-Inf,Inf,solver)
+        @test sol.status == :Optimal
+        @test isapprox(sol.objval, 130/70, atol=1e-6)
+        @test isapprox(norm(sol.sol[1:3] - [0.5714285714285715,0.4285714285714285,0.8571428571428572]), 0.0, atol=1e-6)
+    end
+
+    @testset "QP1" begin
+        # modified from joinpath(Pkg.dir("MathProgBase"), "test", "quadprog.jl"):
+        m = LinearQuadraticModel(solver)
+        loadproblem!(m, [1. 2. 3.; 1. 1. 0.],[-Inf,-Inf,-Inf],[Inf,Inf,Inf],[0.,0.,0.],[4., 1.],[Inf,Inf], :Min)
+
+        setquadobj!(m,diagm([10.0,10.0,10.0]))
+        rows = [1, 2, 2, 2, 3, 3, 3]
+        cols = [1, 1, 1, 2, 2, 3, 3]
+        vals = Float64[2, 0.5, 0.5, 2, 1, 1, 1]
+        setquadobj!(m,rows,cols,vals)
+        m2 = copy(m)
+
+        verify_solution = function (m)
+            stat = status(m)
+            @test stat == :Optimal
+            @test isapprox(getobjval(m), 130/70, atol=1e-6)
+            @test isapprox(norm(getsolution(m) - [0.5714285714285715,0.4285714285714285,0.8571428571428572]), 0.0, atol=1e-6)
+            @test getsolvetime(m) > 0
+        end
+
+        optimize!(m)
+        verify_solution(m)
+
+        setwarmstart!(m2, getsolution(m) .+ 0.1)
+        optimize!(m2)
+        verify_solution(m2)
+    end
+
+    @testset "basic_QP" begin
+        # same QP as test/basic.jl: "basic_QP" but using MathProgBase
+        P = sparse([11. 0.; 0. 0.])
+        q = [3.; 4]
+        A = sparse([-1 0; 0 -1; -1 -3; 2 5; 3 4])
+        u = [0.; 0.; -15; 100; 80]
+        l = fill(-Inf, length(u))
+
+        verify_solution = function (m)
+            tol = 1e-5
+            @show getsolution(m)
+            @show getconstrduals(m)
+            @show getobjval(m)
+            println()
+
+            @test isapprox(norm(getsolution(m) - [0.; 5.]), 0., atol=tol)
+            @test isapprox(norm(getconstrduals(m) - [1.666666666666; 0.; 1.3333333; 0.; 0.]), 0., atol=tol)
+            @test isapprox(getobjval(m), 20., atol=tol)
+        end
+
+        m1 = LinearQuadraticModel(solver)
+        loadproblem!(m1, A, [-Inf, -Inf], [Inf, Inf], q, l, u, :Min)
+        m2 = copy(m1)
+        setquadobj!(m1, P)
+        optimize!(m1)
+        verify_solution(m1)
+
+        setquadobj!(m2, findnz(triu(P))...) # triu to avoid duplicate elements
+        optimize!(m2)
+        verify_solution(m2)
+    end
+
+    @testset "Unsupported behavior" begin
+        m = LinearQuadraticModel(solver)
+        @test_throws ErrorException setsense!(m, :Max)
+        @test_throws ErrorException loadproblem!(m, spzeros(0, 2), [-Inf, -Inf], [Inf, Inf], [1.; 1.], Float64[], Float64[], :Max) # maximization not supported
+        @test_throws ErrorException loadproblem!(m, spzeros(0, 2), [-1., -Inf], [Inf, Inf], [1.; 1.], Float64[], Float64[], :Min) # variable bounds not supported
+        @test_throws ErrorException loadproblem!(m, spzeros(0, 2), [-Inf, -Inf], [10., Inf], [1.; 1.], Float64[], Float64[], :Min) # variable bounds not supported
+
+        loadproblem!(m, [1. 2. 3.; 1. 1. 0.],[-Inf,-Inf,-Inf],[Inf,Inf,Inf],[0.,0.,0.],[4., 1.],[Inf,Inf], :Min)
+        @test_throws ErrorException setvarLB!(m, rand(numvar(m)))
+        @test_throws ErrorException setvarUB!(m, rand(numvar(m)))
+
+        @test_throws ErrorException setvartype!(m, [:Cont, :Bin, :Cont])
+    end
+
+    @testset "Getters/setters" begin
+        m = LinearQuadraticModel(solver)
+        loadproblem!(m, [1. 2. 3.; 1. 1. 0.],[-Inf,-Inf,-Inf],[Inf,Inf,Inf],[0.,0.,0.],[4., 1.],[Inf,Inf], :Min)
+
+        @test getrawsolver(m) isa OSQP.Model
+
+        @test getsense(m) == :Min
+
+        @test getvarLB(m) == [-Inf, -Inf, -Inf]
+        @test getvarUB(m) == [Inf, Inf, Inf]
+
+        @test getconstrLB(m) == [4., 1.]
+        setconstrLB!(m, [-4., Inf])
+        @test getconstrLB(m) == [-4., Inf]
+        @test getconstrUB(m) == [Inf, Inf]
+        setconstrUB!(m, [5., 8.])
+        @test getconstrUB(m) == [5., 8.]
+
+        @test getconstrmatrix(m) == [1. 2. 3.; 1. 1. 0.]
+        @test numlinconstr(m) == 2
+
+        @test getvartype(m) == [:Cont, :Cont, :Cont]
+    end
+end

test/runtests.jl

@@ -11,7 +11,8 @@ tests = [
     "primal_infeasibility.jl",
     "unconstrained.jl",
     "warm_start.jl",
-    "update_matrices.jl"
+    "update_matrices.jl",
+    "mpbinterface.jl"
     ]
 
 println("Running tests:")