Merge pull request #159 from JuliaOpt/ml/rsoc

enable SOCRotated through conic_to_lpqp

Author: mlubin

src/SolverInterface/conic_to_lpqp.jl

@@ -14,9 +14,11 @@ type ConicToLPQPBridge <: AbstractLinearQuadraticModel
     sense::Symbol
     varboundmap::Vector{Int}
     SOCconstrs::Vector{Vector{Int}} # x'x <= y^2, y is first index in vector
+    RSOCconstrs::Vector{Vector{Int}} # x'x <= y*z, y is first index in vector, z is second
+    vartypes::Vector{Symbol}
 end
 
-ConicToLPQPBridge(s::AbstractConicModel) = ConicToLPQPBridge(s, sparse(Int[],Int[],Float64[]), Float64[], Float64[], Float64[], Float64[], Float64[], :Min, Int[], Array{Vector{Int}}(0))
+ConicToLPQPBridge(s::AbstractConicModel) = ConicToLPQPBridge(s, sparse(Int[],Int[],Float64[]), Float64[], Float64[], Float64[], Float64[], Float64[], :Min, Int[], Array{Vector{Int}}(0),Array{Vector{Int}}(0), Symbol[])
 
 export ConicToLPQPBridge
 
@@ -35,25 +37,66 @@ function loadproblem!(wrap::ConicToLPQPBridge, A, collb, colub, obj, rowlb, rowu
     wrap.sense = sense
 end
 
-const notsoc_error = "For conic solvers, only quadratic constraints in second-order conic format (x'x <= y^2) are supported"
+const notsoc_error = "For conic solvers, only quadratic constraints in second-order cone format (x'x <= y^2) or rotated second-order cone format (x'x <= yz) are supported"
 
 function addquadconstr!(wrap::ConicToLPQPBridge, linearidx, linearval, quadrowidx, quadcolidx, quadval, sense, rhs)
     if length(linearidx) > 0 || length(linearval) > 0 || sense != '<' || rhs != 0 || mapreduce(v-> (v == 0.0 || v == 1.0), +, 0, quadval) != length(quadval) - 1 || mapreduce(v->v == -1.0, +, 0, quadval) != 1
         error(notsoc_error)
     end
     length(quadrowidx) == length(quadcolidx) == length(quadval) || error("Inconsistent dimensions")
-    SOCconstr = Int[]
-    for i in 1:length(quadval)
-        quadrowidx[i] == quadcolidx[i] || quadval[i] == 0 || error(notsoc_error)
-        if quadval[i] == 1.0
-            push!(SOCconstr, quadcolidx[i])
-        elseif quadval[i] == -1.0
-            unshift!(SOCconstr, quadcolidx[i])
+    # check if SOC or RSOC or neither
+    n_pos_on_diag = 0
+    off_diag_idx  = 0
+    neg_diag_idx  = 0
+    n = length(quadval)
+    nz = 0
+    for i in 1:n
+        quadval[i] == 0 && continue
+        nz += 1
+        if quadrowidx[i] == quadcolidx[i]
+            if quadval[i] == 1
+                n_pos_on_diag += 1
+            elseif quadval[i] == -1
+                neg_diag_idx = i
+            else
+                error(notsoc_error)
+            end
         else
-            @assert quadval[i] == 0.0
+            if quadval[i] == -1
+                if !(neg_diag_idx == off_diag_idx == 0)
+                    error(notsoc_error)
+                end
+                off_diag_idx = i
+            else
+                error(notsoc_error)
+            end
         end
     end
-    push!(wrap.SOCconstrs, SOCconstr)
+    if n_pos_on_diag == nz-1 && neg_diag_idx > 0
+        # SOC
+        SOCconstr = Int[]
+        push!(SOCconstr, quadcolidx[neg_diag_idx])
+        for i in 1:n
+            if quadval[i] == 1
+                push!(SOCconstr,quadcolidx[i])
+            end
+        end
+        push!(wrap.SOCconstrs, SOCconstr)
+    elseif n_pos_on_diag == nz-1 && off_diag_idx > 0
+        # Rotated SOC
+        RSOCconstr = Int[]
+        push!(RSOCconstr, quadcolidx[off_diag_idx])
+        push!(RSOCconstr, quadrowidx[off_diag_idx])
+        for i in 1:n
+            if quadval[i] == 1
+                push!(RSOCconstr,quadcolidx[i])
+            end
+        end
+        push!(wrap.RSOCconstrs, RSOCconstr)
+
+    else
+        error(notsoc_error)
+    end
 end
 
 
@@ -138,6 +181,7 @@ function optimize!(wrap::ConicToLPQPBridge)
     # Instead, append to A.
     ISOC = Int[]
     JSOC = Int[]
+    VSOC = Float64[]
     SOCconstrs = wrap.SOCconstrs
     soc_row = 0
     for it in 1:length(SOCconstrs)
@@ -148,23 +192,46 @@ function optimize!(wrap::ConicToLPQPBridge)
         for k in 1:n
             push!(ISOC, soc_row+k)
             push!(JSOC, SOCconstrs[it][k])
+            push!(VSOC, -1.0)
             push!(b, 0.0)
         end
         push!(constr_cones, (:SOC, (nrow+extrarows+soc_row+1):(nrow+extrarows+soc_row+n)))
         soc_row += n
     end
+    RSOCconstrs = wrap.RSOCconstrs
+    for it in 1:length(RSOCconstrs)
+        n = length(RSOCconstrs[it])
+        if collb[RSOCconstrs[it][1]] < 0 || collb[RSOCconstrs[it][2]] < 0
+            error("Invalid rotated second-order conic constraint: x'x <= yz requires y >= 0 and z >= 0")
+        end
+        for k in 1:n
+            push!(ISOC, soc_row+k)
+            push!(JSOC, RSOCconstrs[it][k])
+            if k <= 2
+                push!(VSOC, -1.0/sqrt(2))
+            else
+                push!(VSOC, -1.0)
+            end
+            push!(b, 0.0)
+        end
+        push!(constr_cones, (:SOCRotated, (nrow+extrarows+soc_row+1):(nrow+extrarows+soc_row+n)))
+        soc_row += n
+    end
 
 
     A = vcat(sparse(I,J,ones(extrarows),extrarows, nvar),
              A,
-             sparse(ISOC,JSOC,fill(-1.0, length(ISOC)),soc_row,nvar))
+             sparse(ISOC,JSOC,VSOC,soc_row,nvar))
 
     if wrap.sense == :Max
         obj = -obj
     end
 
     #@show obj, full(A), b, constr_cones, var_cones
     loadproblem!(wrap.m, obj, A, b, constr_cones, var_cones)
+    if !isempty(wrap.vartypes)
+        setvartype!(wrap.m, wrap.vartypes)
+    end
     optimize!(wrap.m)
 
 end
@@ -245,3 +312,6 @@ function addconstr!{T<:Integer}(wrap::ConicToLPQPBridge, varidx::AbstractArray{T
     push!(wrap.rowlb, lb)
     push!(wrap.rowub, ub)
 end
+function setvartype!(wrap::ConicToLPQPBridge, v)
+    wrap.vartypes = copy(v)
+end