Using ProximalOperators.SlicedSeperableSum instead of ConeProduct

Author: mfalt

src/FOSSolverInterface.jl

@@ -39,13 +39,13 @@ function loadproblem!(model::FOSMathProgModel, c, A::SparseMatrixCSC, b, constr_
         cone_vars[1] in badcones && error("Cone type $(cone_vars[1]) not supported")
     end
 
-    conesK1 = tuple([conemap[t[1]] for t in constr_cones]...)
-    indexK1 = tuple([t[2] for t in constr_cones]...)
-    model.K1 = ConeProduct(indexK1, conesK1)
+    conesK1 = [conemap[t[1]] for t in constr_cones]
+    indexK1 = [(t[2],) for t in constr_cones]
+    model.K1 = SlicedSeparableSum(conesK1, indexK1)
 
-    conesK2 = tuple([conemap[t[1]] for t in var_cones]...)
-    indexK2 = tuple([t[2] for t in var_cones]...)
-    model.K2 = ConeProduct(indexK2, conesK2)
+    conesK2 = [conemap[t[1]] for t in var_cones]
+    indexK2 = [(t[2],) for t in var_cones]
+    model.K2 = SlicedSeparableSum(conesK2, indexK2)
 
     model.A = A
     # TODO figure out :Min/:Max

src/cones.jl

@@ -14,67 +14,6 @@ const conemap = Dict{Symbol, ProximableFunction}(
 )
 const badcones = ProximableFunction[]
 
-# type DualCone{T<:ProximableFunction} <: ProximableFunction
-#     C::T
-# end
-
-# dual{T<:ProximableFunction}(C::T) = DualCone{T}(C)
-# dual(C::DualCone) = C.C
-
-# function prox!(cone::DualCone, x, y)
-#     prox!(cone.C, -x, y)
-#     for i = 1:length(x)
-#         y[i] = x[i] + y[i]
-#     end
-# end
-
-type ConeProduct{N,T} <: ProximableFunction
-    ranges::NTuple{N, UnitRange{Int64}}
-    cones::T
-end
-#
-ConeProduct() = ConeProduct{0,Any}((),())
-
-# function ConeProduct(ranges::NTuple{N, UnitRange{Int64}}, cones::T) where {N,T}
-#     return ConeProduct{N,T}(ranges, cones)
-# end
-
-toRanges{N}(rangesIn::NTuple{N, UnitRange{Int64}}) = rangesIn
-
-function toRanges{N, T<:Integer}(rangesIn::NTuple{N,Array{T,1}})
-    ranges = Array{UnitRange{Int64},1}(N)
-    for j in 1:N
-        range = rangesIn[j]
-        for (i,el) in enumerate(range[1]:range[end])
-            if el != range[i]
-                error("Invalid range in input")
-            end
-        end
-        ranges[j] = range[1]:range[end]
-    end
-    return tuple(ranges...)
-end
-
-#Should accept tro tuples with UnitRange{Int64} and ProximableFunction
-#Cant enforce types or ambigous with default constructor?
-function ConeProduct(rangesIn, cones)
-    ranges = toRanges(rangesIn)
-    N = length(cones)
-    @assert typeof(ranges) == NTuple{N, UnitRange{Int64}}
-    @assert length(ranges) == N
-    @assert typeof(cones) <: Tuple
-    #Verify that ranges covers the whole range
-    prevRangeEnd = 0
-    for i = 1:N
-      @assert ranges[i][1] == prevRangeEnd + 1
-      prevRangeEnd = ranges[i][end]
-      @assert typeof(cones[i]) <: ProximableFunction
-    end
-    T = typeof(cones)
-    #println("Dumping cones input")
-    #dump(cones)
-    ConeProduct(ranges, cones)
-end
 
 #Wrapper for dual prox to avoid double duals
 function proxDual!(y::AbstractArray, C::ProximableFunction, x::AbstractArray)
@@ -83,15 +22,6 @@ function proxDual!(y::AbstractArray, C::ProximableFunction, x::AbstractArray)
         y[i] = x[i] + y[i]
     end
 end
-# proxDual!(y, C::DualCone, x) = prox!(y, C.C, x)
-
-
-function prox!{N,T}(y::AbstractArray, C::ConeProduct{N,T}, x::AbstractArray)
-    #TODO Paralell implementation
-    for i = 1:N
-        prox!(view(y, C.ranges[i]), C.cones[i], view(x, C.ranges[i]))
-    end
-end
 
 ## Some better dual projections
 const myIndFree = IndFree()
@@ -103,22 +33,38 @@ proxDual!(y::AbstractArray, C::IndNonpositive, x::AbstractArray)   = prox!(y, C,
 # TODO figure out if self dual PSD
 #proxDual!(y::AbstractArray, C::IndPSD, x::AbstractArray)           = prox!(y, C, x)
 
-function proxDual!{N,T}(y::AbstractArray, C::ConeProduct{N,T}, x::AbstractArray)
-    #TODO Paralell implementation
-    for i = 1:N
-        proxDual!(view(y, C.ranges[i]), C.cones[i], view(x, C.ranges[i]))
+# Unroll the loop over the different types of functions to prox on, same as for prox! in ProximalOperators
+# Ignores function value
+@generated function proxDual!{T, A, B, N}(y::AbstractArray{T}, f::SlicedSeparableSum{A, B, N}, x::AbstractArray{T})
+  ex = :()
+  for i = 1:N # For each function type
+    ex = quote $ex;
+      for k in eachindex(f.fs[$i]) # For each function of that type
+        proxDual!(view(y, f.idxs[$i][k]...), f.fs[$i][k], view(x,f.idxs[$i][k]...))
+      end
     end
+  end
+  ex = :($ex; return)
 end
 
 """ Indicator of K2×K1*×R+ × K2*×K1×R+ ∈ (R^n,R^m,R)^2"""
-type DualConeProduct{T1<:ConeProduct,T2<:ConeProduct} <: ProximableFunction
+type DualConeProduct{T1<:SlicedSeparableSum,T2<:SlicedSeparableSum} <: ProximableFunction
     K1::T1
     K2::T2
     m::Int64
     n::Int64
 end
 
-DualConeProduct{T1,T2}(K1::T1,K2::T2) = DualConeProduct{T1,T2}(K1, K2, K1.ranges[end][end], K2.ranges[end][end])
+"""
+Finds the largest index in a SlicedSeparableSum
+"""
+largestindex(e) = error("Something unexpected happend with the product cones with index of type $(typeof(e))")
+largestindex(K::SlicedSeparableSum) = largestindex(K.idxs)                      # Recursive call on indices
+largestindex(l::Union{Tuple, AbstractArray})     = maximum(largestindex.(l))    # Recursive call, apply on each element
+largestindex(l::Union{AbstractUnitRange,Number}) = maximum(l)                   # Possible elements
+
+DualConeProduct{T1,T2}(K1::T1,K2::T2) = DualConeProduct{T1,T2}(K1, K2, largestindex(K1), largestindex(K2))
+
 function prox!(y::AbstractArray, K::DualConeProduct, x::AbstractArray)
     m, n = K.m, K.n
     K1, K2 = K.K1, K.K2

src/types.jl

@@ -26,8 +26,8 @@ end
 type FOSMathProgModel <: AbstractConicModel
     input_numconstr::Int64            # Only needed for interface?
     input_numvar::Int64               # Only needed for interface?
-    K1::ConeProduct
-    K2::ConeProduct
+    K1::SlicedSeparableSum
+    K2::SlicedSeparableSum
     A::SparseMatrixCSC{Float64,Int}   # The A matrix (equalities)
     b::Vector{Float64}                # RHS
     c::Vector{Float64}                # The objective coeffs (always min)
@@ -48,7 +48,7 @@ type FOSMathProgModel <: AbstractConicModel
 end
 
 function FOSMathProgModel(s::FOSAlgorithm; kwargs...)
-    FOSMathProgModel(0, 0, ConeProduct(), ConeProduct(), spzeros(0, 0),
+    FOSMathProgModel(0, 0, SlicedSeparableSum([],[]), SlicedSeparableSum([],[]), spzeros(0, 0),
                      Float64[], Float64[], s, FOSSolverDataPlaceholder(), :NotSolved,
                      0.0, Float64[], Float64[], Float64[], Dict{Symbol,Any}(kwargs), -1,
                      x -> error("No status generator defined"), UInt64(1), ValueHistories.MVHistory())