Fix #889 (#891)

* Refactor accessors on containers/arrays

* Unify behavior of getdual and getvariable: NaN+warn

src/JuMP.jl

@@ -13,6 +13,7 @@ isdefined(Base, :__precompile__) && __precompile__()
 module JuMP
 
 importall Base.Operators
+import Base.map
 
 import MathProgBase
 
@@ -123,7 +124,7 @@ type Model <: AbstractModel
     conDict::Dict{Symbol,Any} # dictionary from constraint names to constraint objects
     varData::ObjectIdDict
 
-    getvalue_counter::Int # number of times we call getvalue on a JuMPContainer, so that we can print out a warning
+    map_counter::Int # number of times we call getvalue, getdual, getlowerbound and getupperbound on a JuMPContainer, so that we can print out a warning
     operator_counter::Int # number of times we add large expressions
 
     # Extension dictionary - e.g. for robust
@@ -177,7 +178,7 @@ function Model(;solver=UnsetSolver(), simplify_nonlinear_expressions::Bool=false
           Dict{Symbol,Any}(),          # varDict
           Dict{Symbol,Any}(),          # conDict
           ObjectIdDict(),              # varData
-          0,                           # getvalue_counter
+          0,                           # map_counter
           0,                           # operator_counter
           Dict{Symbol,Any}(),          # ext
     )
@@ -394,10 +395,12 @@ end
 # internal method that doesn't print a warning if the value is NaN
 _getValue(v::Variable) = v.m.colVal[v.col]
 
+getvaluewarn(v) = Base.warn("Variable value not defined for $(getname(v)). Check that the model was properly solved.")
+
 function getvalue(v::Variable)
     ret = _getValue(v)
     if isnan(ret)
-        Base.warn("Variable value not defined for $(getname(v)). Check that the model was properly solved.")
+        getvaluewarn(v)
     end
     ret
 end
@@ -434,11 +437,19 @@ function getvalue(arr::Array{Variable})
 end
 
 # Dual value (reduced cost) getter
+
+# internal method that doesn't print a warning if the value is NaN
+_getDual(v::Variable) = v.m.redCosts[v.col]
+
+getdualwarn(::Variable) = warn("Variable bound duals (reduced costs) not available. Check that the model was properly solved and no integer variables are present.")
+
 function getdual(v::Variable)
     if length(v.m.redCosts) < MathProgBase.numvar(v.m)
-        error("Variable bound duals (reduced costs) not available. Check that the model was properly solved and no integer variables are present.")
+        getdualwarn(v)
+        NaN
+    else
+        _getDual(v)
     end
-    return v.m.redCosts[v.col]
 end
 
 const var_cats = [:Cont, :Int, :Bin, :SemiCont, :SemiInt]
@@ -538,11 +549,18 @@ LinearConstraint(ref::LinConstrRef) = ref.m.linconstr[ref.idx]::LinearConstraint
 
 linearindex(x::ConstraintRef) = x.idx
 
+# internal method that doesn't print a warning if the value is NaN
+_getDual(c::LinConstrRef) = c.m.linconstrDuals[c.idx]
+
+getdualwarn{T<:Union{ConstraintRef, Int}}(::T) = warn("Dual solution not available. Check that the model was properly solved and no integer variables are present.")
+
 function getdual(c::LinConstrRef)
     if length(c.m.linconstrDuals) != MathProgBase.numlinconstr(c.m)
-        error("Dual solution not available. Check that the model was properly solved and no integer variables are present.")
+        getdualwarn(c)
+        NaN
+    else
+        _getDual(c)
     end
-    return c.m.linconstrDuals[c.idx]
 end
 
 # Returns the number of non-infinity and nonzero bounds on variables
@@ -572,24 +590,18 @@ function getNumBndRows(m::Model)
 end
 
 # Returns the number of second-order cone constraints
-function getNumSOCRows(m::Model)
-    numSOCRows = 0
-    for con in m.socconstr
-        numSOCRows += length(con.normexpr.norm.terms) + 1
-    end
-    return numSOCRows
-end
+getNumRows(c::SOCConstraint) = length(c.normexpr.norm.terms) + 1
+getNumSOCRows(m::Model) = sum(getNumRows.(m.socconstr))
 
 # Returns the number of rows used by SDP constraints in the MPB conic representation
 # (excluding symmetry constraints)
-function getNumSDPRows(m::Model)
-    numSDPRows = 0
-    for con in m.sdpconstr
-        n = size(con.terms, 1)
-        numSDPRows += convert(Int, n*(n+1)/2)
-    end
-    return numSDPRows
+#   Julia seems to not be able to infer the return type (probably because c.terms is Any)
+#   so getNumSDPRows tries to call zero(Any)... Using ::Int solves this issue
+function getNumRows(c::SDConstraint)::Int
+    n = size(c.terms, 1)
+    (n * (n+1)) ÷ 2
 end
+getNumSDPRows(m::Model) = sum(getNumRows.(m.sdpconstr))
 
 # Returns the number of symmetry-enforcing constraints for SDP constraints
 function getNumSymRows(m::Model)
@@ -599,26 +611,36 @@ end
 # Returns the dual variables corresponding to
 # m.sdpconstr[idx] if issdp is true
 # m.socconstr[idx] if sdp is not true
-function getconicdualaux(m::Model, idx, issdp)
+function getconicdualaux(m::Model, idx::Int, issdp::Bool)
     numLinRows = MathProgBase.numlinconstr(m)
     numBndRows = getNumBndRows(m)
     numSOCRows = getNumSOCRows(m)
     numSDPRows = getNumSDPRows(m)
     numSymRows = getNumSymRows(m)
-    if length(m.conicconstrDuals) != (numLinRows + numBndRows + numSOCRows + numSDPRows + numSymRows)
-        error("Dual solution not available. Check that the model was properly solved and no integer variables are present.")
-    end
-    offset = numLinRows + numBndRows
-    if issdp
-        offset += length(m.socconstr)
-    end
-    dual = m.conicconstrDuals[m.constrDualMap[offset + idx]]
-    if issdp
-        offset += length(m.sdpconstr)
-        symdual = m.conicconstrDuals[m.constrDualMap[offset + idx]]
-        dual, symdual
+    numRows = numLinRows + numBndRows + numSOCRows + numSDPRows + numSymRows
+    if length(m.conicconstrDuals) != numRows
+        # solve might not have been called so m.constrDualMap might be empty
+        getdualwarn(idx)
+        c = issdp ? m.sdpconstr[idx] : m.socconstr[idx]
+        duals = fill(NaN, getNumRows(c))
+        if issdp
+            duals, Float64[]
+        else
+            duals
+        end
     else
-        dual
+        offset = numLinRows + numBndRows
+        if issdp
+            offset += length(m.socconstr)
+        end
+        dual = m.conicconstrDuals[m.constrDualMap[offset + idx]]
+        if issdp
+            offset += length(m.sdpconstr)
+            symdual = m.conicconstrDuals[m.constrDualMap[offset + idx]]
+            dual, symdual
+        else
+            dual
+        end
     end
 end
 
@@ -677,13 +699,15 @@ function getdual(c::ConstraintRef{Model,SDConstraint})
             end
         end
     end
-    @assert length(symdual) == length(c.m.sdpconstrSym[c.idx])
-    idx = 0
-    for (i,j) in c.m.sdpconstrSym[c.idx]
-        idx += 1
-        s = symdual[idx]
-        X[i,j] -= s
-        X[j,i] += s
+    if !isempty(symdual)
+        @assert length(symdual) == length(c.m.sdpconstrSym[c.idx])
+        idx = 0
+        for (i,j) in c.m.sdpconstrSym[c.idx]
+            idx += 1
+            s = symdual[idx]
+            X[i,j] -= s
+            X[j,i] += s
+        end
     end
     X
 end
@@ -796,13 +820,14 @@ function getconstraint(m::Model, conname::Symbol)
 end
 
 # usage warnings
-function getvalue_warn(x::JuMPContainer{Variable})
+function mapcontainer_warn(f, x::JuMPContainer{Variable})
     isempty(x) && return
     v = first(values(x))
     m = v.m
-    m.getvalue_counter += 1
-    if m.getvalue_counter > 400
-        Base.warn_once("getvalue has been called on a collection of variables a large number of times. For performance reasons, this should be avoided. Instead of getvalue(x)[a,b,c], use getvalue(x[a,b,c]) to avoid temporary allocations.")
+    m.map_counter += 1
+    if m.map_counter > 400
+        # It might not be f that was called the 400 first times but most probably it is f
+        Base.warn_once("$f has been called on a collection of variables a large number of times. For performance reasons, this should be avoided. Instead of $f(x)[a,b,c], use $f(x[a,b,c]) to avoid temporary allocations.")
     end
     return
 end

src/JuMPContainer.jl

@@ -42,6 +42,9 @@ type JuMPContainerData
     condition
 end
 
+# Needed by getvaluewarn when called by _mapInner
+getname(data::JuMPContainerData) = data.name
+
 #JuMPDict{T,N}(name::AbstractString) =
 #    JuMPDict{T,N}(Dict{NTuple{N},T}(), name)
 
@@ -98,25 +101,40 @@ pushmeta!(x::JuMPContainer, sym::Symbol, val) = (x.meta[sym] = val)
 getmeta(x::JuMPContainer, sym::Symbol) = x.meta[sym]
 
 # duck typing approach -- if eltype(innerArray) doesn't support accessor, will fail
-for accessor in (:getdual, :getlowerbound, :getupperbound)
-    @eval $accessor(x::Union{JuMPContainer,Array}) = map($accessor,x)
+for accessor in (:getdual, :getlowerbound, :getupperbound, :getvalue)
+    @eval $accessor(x::AbstractArray) = map($accessor,x)
+end
+# With JuMPContainer, we take care in _mapInner of the warning if NaN values are returned
+# by the accessor so we use the inner accessor that does not generate warnings
+for (accessor, inner) in ((:getdual, :_getDual), (:getlowerbound, :getlowerbound), (:getupperbound, :getupperbound), (:getvalue, :_getValue))
+    @eval $accessor(x::JuMPContainer) = _map($inner,x)
 end
 
+
 _similar(x::Array) = Array(Float64,size(x))
 _similar{T}(x::Dict{T}) = Dict{T,Float64}()
 
 _innercontainer(x::JuMPArray) = x.innerArray
 _innercontainer(x::JuMPDict)  = x.tupledict
 
-function _getValueInner(x)
+# Warning for getter returning NaN
+function _warnnan(f, data)
+    if f === _getValue
+        getvaluewarn(data)
+    elseif f === _getDual
+        getdualwarn(data)
+    end
+end
+
+function _mapInner(f, x::JuMPContainer)
     vars = _innercontainer(x)
     vals = _similar(vars)
     data = printdata(x)
     warnedyet = false
     for I in eachindex(vars)
-        tmp = _getValue(vars[I])
+        tmp = f(vars[I])
         if isnan(tmp) && !warnedyet
-            warn("Variable value not defined for entry of $(data.name). Check that the model was properly solved.")
+            _warnnan(f, data.name)
             warnedyet = true
         end
         vals[I] = tmp
@@ -127,9 +145,10 @@ end
 JuMPContainer_from(x::JuMPDict,inner) = JuMPDict(inner)
 JuMPContainer_from(x::JuMPArray,inner) = JuMPArray(inner, x.indexsets)
 
-function getvalue(x::JuMPContainer)
-    getvalue_warn(x)
-    ret = JuMPContainer_from(x,_getValueInner(x))
+# The name _map is used instead of map so that this function is called instead of map(::Function, ::JuMPArray)
+function _map(f, x::JuMPContainer)
+    mapcontainer_warn(f, x)
+    ret = JuMPContainer_from(x, _mapInner(f, x))
     # I guess copy!(::Dict, ::Dict) isn't defined, so...
     for (key,val) in x.meta
         ret.meta[key] = val

src/nlp.jl

@@ -54,9 +54,11 @@ function getdual(c::ConstraintRef{Model,NonlinearConstraint})
     initNLP(c.m)
     nldata::NLPData = c.m.nlpdata
     if length(nldata.nlconstrDuals) != length(nldata.nlconstr)
-        error("Dual solution not available. Check that the model was properly solved.")
+        getdualwarn(c)
+        NaN
+    else
+        nldata.nlconstrDuals[c.idx]
     end
-    return nldata.nlconstrDuals[c.idx]
 end
 
 type FunctionStorage

src/solvers.jl

@@ -121,7 +121,7 @@ function solve(m::Model; suppress_warnings=false,
     isempty(kwargs) || error("Unrecognized keyword arguments: $(join([k[1] for k in kwargs], ", "))")
 
     # Clear warning counters
-    m.getvalue_counter = 0
+    m.map_counter = 0
     m.operator_counter = 0
 
     # Remember if the solver was initially unset so we can restore

test/model.jl

@@ -39,14 +39,14 @@ end
         @test_throws ErrorException setobjectivesense(modErr, :Maximum)
         @variable(modErr, errVar)
         @test isnan(getvalue(errVar))
-        @test_throws ErrorException getdual(errVar)
+        @test isnan(getdual(errVar))
 
         modErr = Model()
         @variable(modErr, x, Bin)
         @objective(modErr, Max, x)
         con = @constraint(modErr, x <= 0.5)
         solve(modErr)
-        @test_throws ErrorException getdual(con)
+        @test isnan(getdual(con))
 
         modErr = Model()
         @variable(modErr, 0 <= x <= 1)

test/sdp.jl

@@ -582,7 +582,7 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
         @variable(m, X[1:2,1:2], SDP)
         c = @constraint(m, X[1,1]+X[2,2] == 1)
         @objective(m, Min, 2*X[1,1]+2*X[1,2])
-        @test_throws ErrorException getdual(X)
+        @test all(isnan.(getdual(X)))
         status = solve(m)
 
         @test status == :Optimal
@@ -599,7 +599,7 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
         @variable(m, y)
         c = @SDconstraint(m, [2-y 1; 1 -y] >= 0)
         @objective(m, Max, y)
-        @test_throws ErrorException getdual(c)
+        @test all(isnan(getdual(c)))
         status = solve(m)
 
         @test status == :Optimal
@@ -625,7 +625,7 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
         @variable(m, y)
         c = @SDconstraint(m, [0 y; y 0] <= [1 0; 0 0])
         @objective(m, Max, y)
-        @test_throws ErrorException getdual(c)
+        @test all(isnan(getdual(c)))
         status = solve(m)
 
         @test status == :Optimal
@@ -645,7 +645,7 @@ ispsd(x::JuMP.JuMPArray) = ispsd(x.innerArray)
         @variable(m, X[1:2,1:2], SDP)
         c = @constraint(m, 2*X[1,2] == 1)
         @objective(m, Min, X[1,1])
-        @test_throws ErrorException getdual(X)
+        @test all(isnan(getdual(X)))
         status = solve(m)
 
         @test status == :Optimal

test/variable.jl

@@ -224,3 +224,19 @@ end
         @test MathProgBase.numvar(m) == 4
     end
 end
+
+@testset "[variable] getvalue on sparse array #889" begin
+    m = Model()
+    @variable(m, x)
+    @variable(m, y)
+    X = sparse([1, 3], [2, 3], [x, y])
+
+    @test typeof(X) == SparseMatrixCSC{Variable, Int}
+
+    setvalue(x, 1)
+    setvalue(y, 2)
+
+    Y = getvalue(X)
+    @test typeof(Y) == SparseMatrixCSC{Float64, Int}
+    @test Y == sparse([1, 3], [2, 3], [1, 2])
+end