Use JuMP's functions to simplify conicconstraintdata and add SDP support (#50)

Author: blegat

.travis.yml

@@ -19,4 +19,4 @@ script:
     - julia -e 'Pkg.clone(pwd())'
     - julia --check-bounds=yes -e 'Pkg.test("StructJuMP", coverage=true)'
 after_success:
-  - julia -e 'cd(Pkg.dir("StructJuMP")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
+  - julia -e 'cd(Pkg.dir("StructJuMP")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'

REQUIRE

@@ -1,2 +1,2 @@
 julia 0.3-
-JuMP
+JuMP 0.17

src/BendersBridge.jl

@@ -12,81 +12,34 @@ function conicconstraintdata(m::Model)
     stoch = getStructure(m)
     parent = stoch.parent
     numMasterCols = 0
-    if parent != nothing
+    if parent !== nothing
         numMasterCols = parent.numCols
     end
     v = Symbol[]
-    obj_coeff = zeros(m.numCols)
-    for i in eachindex(m.obj.aff.vars)
-        var = m.obj.aff.vars[i]
-        coeff = m.obj.aff.coeffs[i]
-        obj_coeff[var.col] = coeff
-    end
 
     var_cones = Any[cone for cone in m.varCones]
     con_cones = Any[]
     nnz = 0
 
+    numSDPRows, numSymRows, nnz = JuMP.getSDrowsinfo(m)
+
     linconstr::Vector{LinearConstraint} = m.linconstr
     numLinRows = length(linconstr)
-    numBounds = 0
-    nonNeg  = Int[]
-    nonPos  = Int[]
-    free    = Int[]
-    zeroVar = Int[]
-    for i in 1:m.numCols
-        seen = false
-        lb, ub = m.colLower[i], m.colUpper[i]
-        for (_,cone) in m.varCones
-            if i in cone
-                seen = true
-                @assert lb == -Inf && ub == Inf
-                break
-            end
-        end
-
-        if !seen
-            if lb != -Inf && lb != 0
-                numBounds += 1
-            end
-            if ub != Inf && ub != 0
-                numBounds += 1
-            end
-            if lb == 0 && ub == 0
-                push!(zeroVar, i)
-            elseif lb == 0
-                push!(nonNeg, i)
-            elseif ub == 0
-                push!(nonPos, i)
-            else
-                push!(free, i)
-            end
-        end
-    end
 
-    if !isempty(zeroVar)
-        push!(var_cones, (:Zero,zeroVar))
-    end
-    if !isempty(nonNeg)
-        push!(var_cones, (:NonNeg,nonNeg))
-    end
-    if !isempty(nonPos)
-        push!(var_cones, (:NonPos,nonPos))
-    end
-    if !isempty(free)
-        push!(var_cones, (:Free,free))
-    end
+    numBounds = JuMP.variable_range_to_cone!(var_cones, m)
 
     nnz += numBounds
     for c in 1:numLinRows
         nnz += length(linconstr[c].terms.coeffs)
     end
 
-    numSOCRows = 0
-    for con in m.socconstr
-        numSOCRows += length(con.normexpr.norm.terms) + 1
-    end
-    numRows = numLinRows + numBounds + numSOCRows
+    numSOCRows = JuMP.getNumSOCRows(m)
+    numNormRows = length(m.socconstr)
+
+    numRows = numLinRows + numBounds + numSOCRows + numSDPRows + numSymRows
+
+    # constr_to_row is not used but fill_bounds_constr! and fillconstr! for SDP needs them
+    constr_to_row = Array{Vector{Int}}(numBounds + 2*length(m.sdpconstr))
 
     b = Array(Float64, numRows)
 
@@ -98,134 +51,56 @@ function conicconstraintdata(m::Model)
     V_s = Float64[]
 
     # Fill it up
-    nnz = 0
     tmprow = JuMP.IndexedVector(Float64,m.numCols)
-    tmpelts = tmprow.elts
-    tmpnzidx = tmprow.nzidx
-    nonneg_rows = Int[]
-    nonpos_rows = Int[]
-    eq_rows     = Int[]
-    for c in 1:numLinRows
-        if linconstr[c].lb == -Inf
-            b[c] = linconstr[c].ub
-            push!(nonneg_rows, c)
-        elseif linconstr[c].ub == Inf
-            b[c] = linconstr[c].lb
-            push!(nonpos_rows, c)
-        elseif linconstr[c].lb == linconstr[c].ub
-            b[c] = linconstr[c].lb
-            push!(eq_rows, c)
-        else
-            error("We currently do not support ranged constraints with conic solvers")
-        end
 
-        JuMP.assert_isfinite(linconstr[c].terms)
-        coeffs = linconstr[c].terms.coeffs
-        vars = linconstr[c].terms.vars
-        # eliminated collect duplicates
-        for ind in eachindex(coeffs)
-            if vars[ind].m === parent
-                push!(I_m, c)
-                push!(J_m, vars[ind].col)
-                push!(V_m, coeffs[ind])
-            else
-                push!(I_s, c)
-                push!(J_s, vars[ind].col)
-                push!(V_s, coeffs[ind])
-            end
-        end
+    JuMP.fillconstrRHS!(b, con_cones, 0, m.linconstr)
+    if numMasterCols > 0
+        JuMP.fillconstrLHS!(I_m, J_m, V_m, tmprow, 0, m.linconstr, parent, true)
     end
+    c = JuMP.fillconstrLHS!(I_s, J_s, V_s, tmprow, 0, m.linconstr, m, true)
 
-    c = numLinRows
-    bndidx = 0
     for idx in 1:m.numCols
-        lb = m.colLower[idx]
         # identify integrality information
         push!(v, m.colCat[idx])
-        if lb != -Inf && lb != 0
-            bndidx += 1
-            nnz += 1
-            c   += 1
-            push!(I_s, c)
-            push!(J_s, idx)
-            push!(V_s, 1.0)
-            b[c] = lb
-            push!(nonpos_rows, c)
-        end
-        ub = m.colUpper[idx]
-        if ub != Inf && ub != 0
-            bndidx += 1
-            c   += 1
-            push!(I_s, c)
-            push!(J_s, idx)
-            push!(V_s, 1.0)
-            b[c] = ub
-            push!(nonneg_rows, c)
-        end
     end
+    c, d = JuMP.fill_bounds_constr!(I_s, J_s, V_s, b, con_cones, constr_to_row, c, 0, m)
+
+    @assert c == numLinRows + numBounds
+    @assert d == numBounds
 
-    if !isempty(nonneg_rows)
-        push!(con_cones, (:NonNeg,nonneg_rows))
+    JuMP.fillconstrRHS!(b, con_cones, c, m.socconstr)
+    if numMasterCols > 0
+        JuMP.fillconstrLHS!(I_m, J_m, V_m, tmprow, c, m.socconstr, parent, true)
     end
-    if !isempty(nonpos_rows)
-        push!(con_cones, (:NonPos,nonpos_rows))
+    c = JuMP.fillconstrLHS!(I_s, J_s, V_s, tmprow, c, m.socconstr, m, true)
+
+    @assert c == numLinRows + numBounds + numSOCRows
+
+    if numMasterCols > 0
+        c, d = JuMP.fillconstr!(I_m, J_m, V_m, b, con_cones, tmprow, constr_to_row, c, d, m.sdpconstr, m, true)
     end
-    if !isempty(eq_rows)
-        push!(con_cones, (:Zero,eq_rows))
+    c, d = JuMP.fillconstr!(I_s, J_s, V_s, b, con_cones, tmprow, constr_to_row, c, d, m.sdpconstr, m, true)
+
+    if c < length(b)
+        # This happens for example when symmetry constraints are dropped with SDP
+        resize!(b, c)
     end
-    @assert c == numLinRows + numBounds
 
-    tmpelts = tmprow.elts
-    tmpnzidx = tmprow.nzidx
-    socidx = 0
-    for con in m.socconstr
-        socidx += 1
-        expr = con.normexpr
-        c += 1
-        soc_start = c
-        JuMP.collect_expr!(m, tmprow, expr.aff)
-        nnz = tmprow.nnz
-        indices = tmpnzidx[1:nnz]
-        vars = expr.aff.vars
-        for i in eachindex(vars)
-            if vars[i].m === parent
-                push!(I_m, c)
-                push!(J_m, indices[i])
-                push!(V_m, tmpelts[indices[i]])
-            else
-                push!(I_s, c)
-                push!(J_s, indices[i])
-                push!(V_s, tmpelts[indices[i]])
-            end
-        end
-        b[c] = -expr.aff.constant
-        for term in expr.norm.terms
-            c += 1
-            JuMP.collect_expr!(m, tmprow, term)
-            nnz = tmprow.nnz
-            indices = tmpnzidx[1:nnz]
-            vars = term.vars
-            for i = 1:length(vars)
-                if vars[i].m == parent
-                    push!(I_m, c)
-                    push!(J_m, indices[i])
-                    push!(V_m, -expr.coeff*tmpelts[indices[i]])
-                else
-                    push!(I_s, c)
-                    push!(J_s, indices[i])
-                    push!(V_s, -expr.coeff*tmpelts[indices[i]])
-                end
-            end
-            b[c] = expr.coeff*term.constant
-        end
-        push!(con_cones, (:SOC, soc_start:c))
+    f_s = JuMP.prepAffObjective(m)
+
+    # The conic MPB interface defines conic problems as
+    # always being minimization problems, so flip if needed
+    m.objSense == :Max && scale!(f_s, -1.0)
+
+    if numMasterCols > 0
+        JuMP.rescaleSDcols!(spzeros(numMasterCols), J_m, V_m, parent)
     end
-    @assert c == numLinRows + numBounds + numSOCRows
+    JuMP.rescaleSDcols!(f_s, J_s, V_s, m)
 
     A = sparse(I_m, J_m, V_m, numRows, numMasterCols)
     B = sparse(I_s, J_s, V_s, numRows, m.numCols)
-  
-    return obj_coeff, A, B, b, var_cones, con_cones, v
+
+    return f_s, A, B, b, var_cones, con_cones, v
 end
 
 function BendersBridge(m::Model, master_solver, sub_solver)