Skip to content

Commit

Permalink
Add support for constraint dual (#1129)
Browse files Browse the repository at this point in the history 
* Add support for constraint dual

* Use MOICON in constraint dict

* Update comment
  • Loading branch information
@blegat @mlubin
blegat authored and mlubin committed Nov 23, 2017
1 parent db33c82 commit d3e40d9
Show file tree
Hide file tree
Showing 5 changed files with 116 additions and 51 deletions.
25 changes: 24 additions & 1 deletion src/JuMP.jl
View file
Expand Up @@ -46,7 +46,7 @@ export
setlowerbound, setupperbound,
#getlowerbound, getupperbound,
#getvalue, setvalue,
getdual,
#getdual,
#setcategory, getcategory,
setstart,
linearindex,
Expand Down Expand Up @@ -100,6 +100,10 @@ mutable struct Model <: AbstractModel
# obj#::QuadExpr
# objSense::Symbol

# Mapping from the constraint reference in `instance`
# and the constraint reference in `solverinstance`
constrainttosolverconstraint::Dict{MOICON,MOICON}

# linconstr#::Vector{LinearConstraint}
# quadconstr
# sosconstr
Expand Down Expand Up @@ -524,6 +528,25 @@ end

# linearindex(x::ConstraintRef) = x.idx

function solverinstanceref(cr::ConstraintRef{Model, MOICON{F, S}}) where {F, S}
cr.m.constrainttosolverconstraint[cr.instanceref]::MOICON{F, S}
end

function hasresultdual(cr::ConstraintRef{Model, <:MOICON})
MOI.get(cr.m.solverinstance, MOI.ConstraintDual(), cr.instanceref)
end

"""
resultdual(cr::ConstraintRef)
Get the dual value of this constraint in the result returned by a solver.
Use `hasresultdual` to check if a result exists before asking for values.
Replaces `getdual` for most use cases.
"""
function resultdual(cr::ConstraintRef{Model, MOICON{F, S}}) where {F, S}
MOI.get(cr.m.solverinstance, MOI.ConstraintDual(), solverinstanceref(cr))
end

###############################################################################
# GenericAffineExpression, AffExpr, AffExprConstraint
include("affexpr.jl")
Expand Down
25 changes: 21 additions & 4 deletions src/solverinterface.jl
View file
Expand Up @@ -3,6 +3,22 @@
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.


"""
copyconstraints!(m::JuMP.Model, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
Transfer the constraints of type `F`-in-`S` to the solver instance.
"""
function copyconstraints!(m::Model, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
for cref in MOI.get(m.instance, MOI.ListOfConstraintReferences{F, S}())
f = MOI.get(m.instance, MOI.ConstraintFunction(), cref)
s = MOI.get(m.instance, MOI.ConstraintSet(), cref)
solvercref = MOI.addconstraint!(m.solverinstance, f, s)
@assert !haskey(m.constrainttosolverconstraint, cref)
m.constrainttosolverconstraint[cref] = solvercref
end
end

"""
attach(m::JuMP.Model)
Expand All @@ -16,6 +32,7 @@ function attach(m::Model, solverinstance::MOI.AbstractSolverInstance)
solvervariables = MOI.addvariables!(m.solverinstance, numvar(m)) # TODO numvar shouldn't return unsigned int

m.variabletosolvervariable = Dict{MOIVAR,MOIVAR}()
m.constrainttosolverconstraint = Dict{UInt64,UInt64}()
# TODO: replace with ListOfVariableReferences()
# Now we're assuming all instance variables are numbered sequentially
for i in 1:numvar(m)
Expand All @@ -25,10 +42,10 @@ function attach(m::Model, solverinstance::MOI.AbstractSolverInstance)
MOI.set!(m.solverinstance, MOI.ObjectiveSense(), MOI.get(m.instance, MOI.ObjectiveSense()))
MOI.set!(m.solverinstance, MOI.ObjectiveFunction(), MOI.get(m.instance, MOI.ObjectiveFunction()))

# TODO: replace with ListOfConstraintReferences()
# This would be a bit more transparent than the call below
# TODO: keep track of constraint references!
MOIU.broadcastcall( constrs -> for (cref, f, s) in constrs; MOI.addconstraint!(m.solverinstance, f, s) end, m.instance)
for (F, S) in MOI.get(m.instance, MOI.ListOfConstraints())
# do the rest in copyconstraints! which is type stable
copyconstraints!(m, F, S)
end

m.solverinstanceattached = true
nothing
Expand Down
3 changes: 2 additions & 1 deletion test/lp.jl
View file
Expand Up @@ -10,7 +10,7 @@
@variable(m, y >= 0.0)
@objective(m, Min, -x)

@constraint(m, x + y <= 1)
c = @constraint(m, x + y <= 1)

JuMP.attach(m, CSDPInstance(printlevel=0))
JuMP.solve(m)
Expand All @@ -27,5 +27,6 @@
@test JuMP.resultvalue(x + y) ≈ 1.0 atol=1e-6
@test JuMP.objectivevalue(m) ≈ -1.0 atol=1e-6

@test JuMP.resultdual(c) ≈ -1 atol=1e-6
end
end
100 changes: 59 additions & 41 deletions test/sdp.jl
View file
Expand Up @@ -58,6 +58,7 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
@test JuMP.hasvariableresult(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test JuMP.objectivevalue(m) ≈ 0.705710509 atol=1e-6

Expand Down Expand Up @@ -242,11 +243,12 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 3, atol=1e-5)
#@test isapprox(getdual(c1), 0, atol=1e-5)
#@test isapprox(getdual(c2), 1, atol=1e-5)
#@test isapprox(getdual(Y), [1 0 0; 0 0 0; 0 0 1], atol=1e-5)
@test isapprox(JuMP.resultdual(c1), 0, atol=1e-5)
@test isapprox(JuMP.resultdual(c2), 1, atol=1e-5)
#@test isapprox(JuMP.resultdual(Y), [1 0 0; 0 0 0; 0 0 1], atol=1e-5)
end

# min Y[1,2] max y
Expand All @@ -265,13 +267,14 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 1, atol=1e-4)
Yval = JuMP.resultvalue.(Y)
@test isapprox(Yval[1,2], 1, atol=1e-4)
@test isapprox(Yval[2,1], 1, atol=1e-4)
#@test isapprox(getdual(c), 1, atol=1e-5)
#@test isapprox(getdual(Y), zeros(3,3), atol=1e-4)
@test isapprox(JuMP.resultdual(c), 1, atol=1e-5)
#@test isapprox(JuMP.resultdual(Y), zeros(3,3), atol=1e-4)
end

# min x + Y[1,1] max y + z
Expand All @@ -294,14 +297,15 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 1, atol=1e-3)
@test isapprox(JuMP.resultvalue(x), 1, atol=1e-5)
@test isapprox(JuMP.resultvalue.(Y)[1,1], 0, atol=1e-4)
#@test isapprox(getdual(x), 0, atol=1e-5)
#@test isapprox(getdual(Y), [1 0 0; 0 0 0; 0 0 0], atol=1e-4)
#@test isapprox(getdual(c1), 1, atol=1e-5)
#@test isapprox(getdual(c2), 0, atol=1e-4)
#@test isapprox(JuMP.resultdual(x), 0, atol=1e-5)
#@test isapprox(JuMP.resultdual(Y), [1 0 0; 0 0 0; 0 0 0], atol=1e-4)
@test isapprox(JuMP.resultdual(c1), 1, atol=1e-5)
@test isapprox(JuMP.resultdual(c2), 0, atol=1e-4)
end

# function nuclear_norm(model, A)
Expand Down Expand Up @@ -346,6 +350,7 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 4, atol=1e-5)
end
Expand Down Expand Up @@ -389,6 +394,7 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 2, atol=1e-5)
end
Expand Down Expand Up @@ -464,21 +470,24 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
@variable(m, Q[1:2, 1:2], PSD)
c1 = @constraint(m, Q[1,1] - 1 == Q[2,2])
@variable(m, objective)
@SDconstraint(m, [1 Q[1,1]; Q[1,1] objective] ⪰ 0)
c2 = @SDconstraint(m, [1 Q[1,1]; Q[1,1] objective] ⪰ 0)
@objective(m, Min, objective)

JuMP.attach(m, CSDPInstance(printlevel=0))
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test JuMP.resultvalue.(Q) ≈ [1 0; 0 0] atol=1e-3
@test JuMP.objectivevalue(m) ≈ 1 atol=1e-4
@test JuMP.resultvalue(objective) ≈ 1 atol=1e-4
#@test getdual(objective) ≈ 0 atol=1e-5
#@test getdual(Q) ≈ [0 0; 0 2] atol=1e-3
#@test getdual(c1) ≈ 2 atol=1e-4 # y
#@test getdual(c2) ≈ [-1 1; 1 -1] atol=1e-3 # X
#@test JuMP.resultdual(objective) ≈ 0 atol=1e-5
#@test JuMP.resultdual(Q) ≈ [0 0; 0 2] atol=1e-3
@test JuMP.resultdual(c1) ≈ 2 atol=1e-4 # y
@test JuMP.resultdual(c2) ≈ [1, -1, 1] atol=1e-3 # X
# TODO
#@test JuMP.resultdual(c2) ≈ [1 -1; -1 1] atol=1e-3 # X
end

# The four following tests are from Example 2.11, Example 2.13 and Example 2.27 of:
Expand All @@ -492,17 +501,18 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
@variable(m, X[1:2,1:2], PSD)
c = @constraint(m, X[1,1]+X[2,2] == 1)
@objective(m, Min, 2*X[1,1]+2*X[1,2])
# @test all(isnan.(getdual(X)))
# @test all(isnan.(JuMP.resultdual(X)))

JuMP.attach(m, CSDPInstance(printlevel=0))
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test JuMP.objectivevalue(m) ≈ 1-sqrt(2) atol=1e-5
@test JuMP.resultvalue.(X) ≈ [(2-sqrt(2))/4 -1/(2*sqrt(2)); -1/(2*sqrt(2)) (2+sqrt(2))/4] atol=1e-4
# @test getdual(X) ≈ [1+sqrt(2) 1; 1 sqrt(2)-1] atol=1e-4
# @test getdual(c) ≈ 1-sqrt(2) atol=1e-5
# @test JuMP.resultdual(X) ≈ [1+sqrt(2) 1; 1 sqrt(2)-1] atol=1e-4
@test JuMP.resultdual(c) ≈ 1-sqrt(2) atol=1e-5
end

# Example 2.13
Expand All @@ -512,19 +522,22 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
@variable(m, y)
c = @SDconstraint(m, [2-y 1; 1 -y] >= 0)
@objective(m, Max, y)
#@test all(isnan, getdual(c))
#@test all(isnan, JuMP.resultdual(c))

JuMP.attach(m, CSDPInstance(printlevel=0))
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 1-sqrt(2), atol=1e-5)
@test isapprox(JuMP.resultvalue(y), 1-sqrt(2), atol=1e-5)

#X = getdual(c)
#@test isapprox(getdual(c), [(2-sqrt(2))/4 -1/(2*sqrt(2)); -1/(2*sqrt(2)) (2+sqrt(2))/4], atol=1e-4)
#@test isapprox(getdual(y), 0, atol=1e-5)
X = JuMP.resultdual(c)
@test isapprox(JuMP.resultdual(c), [(2-sqrt(2))/4, -1/(2*sqrt(2)), (2+sqrt(2))/4], atol=1e-4)
# TODO
#@test isapprox(JuMP.resultdual(c), [(2-sqrt(2))/4 -1/(2*sqrt(2)); -1/(2*sqrt(2)) (2+sqrt(2))/4], atol=1e-4)
#@test isapprox(JuMP.resultdual(y), 0, atol=1e-5)
end

# Example 2.27
Expand All @@ -541,21 +554,25 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
@variable(m, y)
c = @SDconstraint(m, [0 y; y 0] <= [1 0; 0 0])
@objective(m, Max, y)
#@test all(isnan, getdual(c))
#@test all(isnan, JuMP.resultdual(c))

JuMP.attach(m, CSDPInstance(printlevel=0))
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 0, atol=1e-5)
@test isapprox(JuMP.resultvalue(y), 0, atol=1e-5)

#X = getdual(c)
X = JuMP.resultdual(c)
@test isapprox(X[1], 0, atol=1e-5)
@test isapprox(X[2], 1/2, atol=1e-5)
# TODO
#@test isapprox(X[1,1], 0, atol=1e-5)
#@test isapprox(X[1,2], 1/2, atol=1e-5)
#@test isapprox(X[2,1], 1/2, atol=1e-5)
#@test isapprox(getdual(y), 0, atol=1e-5)
#@test isapprox(JuMP.resultdual(y), 0, atol=1e-5)
end

@testset "SDP with primal solution not attained" begin
Expand All @@ -564,24 +581,22 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
@variable(m, X[1:2,1:2], PSD)
c = @constraint(m, 2*X[1,2] == 1)
@objective(m, Min, X[1,1])
# @test all(isnan, getdual(X))
# @test all(isnan, JuMP.resultdual(X))

JuMP.attach(m, CSDPInstance(printlevel=0))
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint

@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test JuMP.objectivevalue(m) ≈ 0 atol=1e-5
Xval = JuMP.resultvalue.(X)
@test Xval[1,1] ≈ 0 atol=1e-5
@test Xval[1,2] ≈ 1/2 atol=1e-5
@test Xval[2,1] ≈ 1/2 atol=1e-5

# @test isapprox(getdual(X), [1 0; 0 0], atol=1e-4)
# @test isapprox(getdual(c), 0, atol=1e-5)
# @test isapprox(JuMP.resultdual(X), [1 0; 0 0], atol=1e-4)
@test isapprox(JuMP.resultdual(c), 0, atol=1e-4)
end

# # min X[1,1] max y/2+z/2
Expand All @@ -603,12 +618,12 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
# @test isapprox(JuMP.resultvalue(y), 0, atol=1e-5)
# @test isapprox(JuMP.resultvalue(z), 0, atol=1e-5)
#
# #X = getdual(c)
# #X = JuMP.resultdual(c)
# #@test isapprox(X[1,1], 0, atol=1e-5)
# #@test isapprox(X[1,2], 1/2, atol=1e-5)
# #@test isapprox(X[2,1], 1/2, atol=1e-5)
# #@test isapprox(getdual(y), 0, atol=1e-5)
# #@test isapprox(getdual(z), 0, atol=1e-5)
# #@test isapprox(JuMP.resultdual(y), 0, atol=1e-5)
# #@test isapprox(JuMP.resultdual(z), 0, atol=1e-5)
# end

# # min X[1,1] max y
Expand All @@ -631,12 +646,12 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
# @test isapprox(JuMP.resultvalue(y), 0, atol=1e-5)
# @test isapprox(JuMP.resultvalue(z), 0, atol=1e-5)
#
# #X = getdual(c)
# #X = JuMP.resultdual(c)
# #@test isapprox(X[1,1], 0, atol=1e-5)
# #@test isapprox(X[1,2], 1, atol=1e-5) # X is not symmetric !
# #@test isapprox(X[2,1], 0, atol=1e-5)
# #@test isapprox(getdual(y), 0, atol=1e-5)
# #@test isapprox(getdual(z), 0, atol=1e-5)
# #@test isapprox(JuMP.resultdual(y), 0, atol=1e-5)
# #@test isapprox(JuMP.resultdual(z), 0, atol=1e-5)
# end

@testset "Nonzero dual for a scalar variable with sdp solver" begin
Expand All @@ -653,16 +668,19 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3
JuMP.solve(m)
@test JuMP.terminationstatus(m) == MOI.Success
@test JuMP.primalstatus(m) == MOI.FeasiblePoint
@test JuMP.dualstatus(m) == MOI.FeasiblePoint

@test isapprox(JuMP.objectivevalue(m), 10/3, atol=1e-5)
@test isapprox(JuMP.resultvalue(x1), 5/3, atol=1e-5)
@test isapprox(JuMP.resultvalue(x2), 0, atol=1e-5)
@test isapprox(JuMP.resultvalue(x3), 1/3, atol=1e-5)

#@test isapprox(getdual(c), [-2/3 0; 0 -2/3], atol=1e-4)
#@test isapprox(getdual(x1), 0, atol=1e-5)
#@test isapprox(getdual(x2), 1/3, atol=1e-5)
#@test isapprox(getdual(x3), 0, atol=1e-5)
@test isapprox(JuMP.resultdual(c), [2/3, 0, 2/3], atol=1e-4)
# TODO
#@test isapprox(JuMP.resultdual(c), [-2/3 0; 0 -2/3], atol=1e-4)
#@test isapprox(JuMP.resultdual(x1), 0, atol=1e-5)
#@test isapprox(JuMP.resultdual(x2), 1/3, atol=1e-5)
#@test isapprox(JuMP.resultdual(x3), 0, atol=1e-5)
end


Expand All @@ -678,7 +696,7 @@ ispsd(x::Matrix) = minimum(eigvals(x)) ≥ -1e-3

@test abs(JuMP.objectivevalue(m)) < 1e-5
@test norm(JuMP.resultvalue.(X)) < 1e-5
#@test isapprox(getdual(X), eye(3), atol=1e-5)
#@test isapprox(JuMP.resultdual(X), eye(3), atol=1e-5)
end

end

0 comments on commit d3e40d9

Please sign in to comment.