Simplify conic

src/MatrixOptInterface.jl

@@ -65,7 +65,7 @@ end
 
 @enum VariableType CONTINUOUS INTEGER BINARY
 
-include("sparse_matrix.jl")
+include("product_of_sets.jl")
 include("conic_form.jl")
 include("matrix_input.jl")
 #include("change_form.jl")

src/conic_form.jl

@@ -4,255 +4,48 @@
 # in the LICENSE.md file or at https://opensource.org/licenses/MIT.
 
 """
-    GeometricConicForm{T, AT, VT, C} <: MOI.ModelLike
+    empty_geometric_conic_form
 
 Represents an optimization model of the form:
 ```
 sense ⟨c, x⟩ + c0
-s.t.  b_i - A_i x ∈ C_i ∀ i
+s.t.  A_i x + b_i ∈ C_i ∀ i
 ```
 with each `C_i` a cone defined in MOI.
 """
-mutable struct GeometricConicForm{T,AT,VB,VC,C} <: MOI.ModelLike
-    num_rows::Vector{Int}
-    dimension::Dict{Int,Int}
-    sense::MOI.OptimizationSense
-    objective_constant::T # The objective
-    A::Union{Nothing,AT} # The constraints
-    b::VB          # `b - Ax in cones`
-    c::VC          # `sense c'x + objective_constant`
-    cone_types::C
-    cone_types_dict::Dict{DataType,Int}
-
-    function GeometricConicForm{T,AT,VB,VC}(cone_types) where {T,AT,VB,VC}
-        model = new{T,AT,VB,VC,typeof(cone_types)}()
-        model.cone_types = cone_types
-        model.cone_types_dict =
-            Dict{DataType,Int}(s => i for (i, s) in enumerate(cone_types))
-        model.num_rows = zeros(Int, length(cone_types))
-        model.dimension = Dict{Int,Int}()
-        model.A = nothing
-        return model
-    end
-end
-
-function GeometricConicForm{T,AT,VT}(cone_types) where {T,AT,VT}
-    return GeometricConicForm{T,AT,VT,VT}(cone_types)
-end
-
-_set_type(::MOI.ConstraintIndex{F,S}) where {F,S} = S
-
-MOI.is_empty(model::GeometricConicForm) = model.A === nothing
-
-function MOI.empty!(model::GeometricConicForm{T}) where {T}
-    empty!(model.dimension)
-    fill!(model.num_rows, 0)
-    model.A = nothing
-    model.sense = MOI.FEASIBILITY_SENSE
-    return model.objective_constant = zero(T)
-end
-
-function MOI.supports_constraint(
-    model::GeometricConicForm{T},
-    ::Type{MOI.VectorAffineFunction{T}},
-    ::Type{S},
-) where {T,S<:MOI.AbstractVectorSet}
-    return haskey(model.cone_types_dict, S)
-end
-
-function _allocate_variables(
-    model::GeometricConicForm{T,AT,VT},
-    vis_src,
-    idxmap,
-) where {T,AT,VT}
-    model.A = AT(length(vis_src))
-    for (i, vi) in enumerate(vis_src)
-        idxmap[vi] = MOI.VariableIndex(i)
-    end
-    return
-end
-
-function rows(
-    model::GeometricConicForm{T},
-    ci::CI{MOI.VectorAffineFunction{T}},
-) where {T}
-    return ci.value .+ (1:model.dimension[ci.value])
-end
-
-function MOI.set(
-    ::GeometricConicForm,
-    ::MOI.VariablePrimalStart,
-    ::MOI.VariableIndex,
-    ::Nothing,
-)
-    return
-end
-
-function MOI.set(
-    model::GeometricConicForm{T},
-    ::MOI.VariablePrimalStart,
-    vi::MOI.VariableIndex,
-    value::T,
-) where {T}
-    return model.primal[vi.value] = value
-end
-
-function MOI.set(
-    ::GeometricConicForm,
-    ::MOI.ConstraintPrimalStart,
-    ::MOI.ConstraintIndex,
-    ::Nothing,
-)
-    return
-end
-
-function MOI.set(
-    model::GeometricConicForm,
-    ::MOI.ConstraintPrimalStart,
-    ci::MOI.ConstraintIndex,
-    value,
-)
-    offset = constroffset(model, ci)
-    return model.slack[rows(model, ci)] .= value
-end
-
-function MOI.set(
-    ::GeometricConicForm,
-    ::MOI.ConstraintDualStart,
-    ::MOI.ConstraintIndex,
-    ::Nothing,
-)
-    return
-end
-
-function MOI.set(
-    model::GeometricConicForm,
-    ::MOI.ConstraintDualStart,
-    ci::MOI.ConstraintIndex,
-    value,
-)
-    offset = constroffset(model, ci)
-    return model.dual[rows(model, ci)] .= value
-end
-
-function MOI.set(
-    model::GeometricConicForm,
-    ::MOI.ObjectiveSense,
-    sense::MOI.OptimizationSense,
-)
-    return model.sense = sense
-end
-
-variable_index_value(t::MOI.ScalarAffineTerm) = t.variable.value
-
-function variable_index_value(t::MOI.VectorAffineTerm)
-    return variable_index_value(t.scalar_term)
-end
-
-function MOI.set(
-    model::GeometricConicForm{T},
-    ::MOI.ObjectiveFunction,
-    f::MOI.ScalarAffineFunction{T},
-) where {T}
-    c = Vector(
-        sparsevec(
-            variable_index_value.(f.terms),
-            MOI.coefficient.(f.terms),
-            model.A.n,
-        ),
+function empty_geometric_conic_form(
+    cones;
+    Tv = Float64,
+    Ti = Int,
+    I = MOI.Utilities.OneBasedIndexing,
+)
+    model = MOI.Utilities.GenericModel{Tv}(
+        MOI.Utilities.ObjectiveContainer{Tv}(),
+        MOI.Utilities.FreeVariables(),
+        MOI.Utilities.MatrixOfConstraints{
+            Tv,
+            MOI.Utilities.MutableSparseMatrixCSC{Tv,Ti,I},
+            Vector{Tv},
+            ProductOfSets{Tv},
+        }(),
     )
-    model.objective_constant = f.constant
-    model.c = c
-    return
+    set_set_types(model.constraints.sets, cones)
+    return model
 end
 
-function _allocate_constraint(
-    model::GeometricConicForm,
-    src,
-    indexmap,
-    cone_id,
-    ci,
-)
-    # TODO use `CanonicalConstraintFunction`
-    func = MOI.get(src, MOI.ConstraintFunction(), ci)
-    func = MOIU.is_canonical(func) ? func : MOI.Utilities.canonical(func)
-    allocate_terms(model.A, indexmap, func)
-    offset = model.num_rows[cone_id]
-    model.num_rows[cone_id] = offset + MOI.output_dimension(func)
-    return ci, offset, func
+function geometric_conic_form(model::MOI.ModelLike, cones; kws...)
+    form = empty_geometric_conic_form(cones; kws...)
+    index_map = MOI.copy_to(form, model)
+    return form, index_map
 end
 
-function _allocate_constraints(
-    model::GeometricConicForm{T},
-    src,
-    indexmap,
-    cone_id,
-    ::Type{S},
-) where {T,S}
-    cis = MOI.get(
-        src,
-        MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T},S}(),
-    )
-    return map(cis) do ci
-        return _allocate_constraint(model, src, indexmap, cone_id, ci)
-    end
-end
-
-function _load_variables(model::GeometricConicForm, nvars::Integer)
-    m = sum(model.num_rows)
-    model.A.m = m
-    model.b = zeros(m)
-    model.c = zeros(model.A.n)
-    return allocate_nonzeros(model.A)
-end
+_coef_type(::MOI.Utilities.AbstractModel{T}) where {T} = T
 
-function _load_constraints(
-    model::GeometricConicForm,
-    src,
-    indexmap,
-    cone_offset,
-    i,
-    cache,
-)
-    for (ci_src, offset_in_cone, func) in cache
-        offset = cone_offset + offset_in_cone
-        set = MOI.get(src, MOI.ConstraintSet(), ci_src)
-        load_terms(model.A, indexmap, func, offset)
-        copyto!(model.b, offset + 1, func.constants)
-        model.dimension[offset] = MOI.output_dimension(func)
-        indexmap[ci_src] = typeof(ci_src)(offset)
-    end
-end
-
-function MOI.copy_to(dest::GeometricConicForm{T}, src::MOI.ModelLike) where {T}
-    MOI.empty!(dest)
-    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
-    idxmap = MOIU.IndexMap()
-    has_constraints = BitSet()
-    for (F, S) in MOI.get(src, MOI.ListOfConstraintTypesPresent())
-        i = get(dest.cone_types_dict, S, nothing)
-        if i === nothing || F != MOI.VectorAffineFunction{T}
-            throw(MOI.UnsupportedConstraint{F,S}())
-        end
-        push!(has_constraints, i)
-    end
-    _allocate_variables(dest, vis_src, idxmap)
-    # Allocate constraints
-    caches = map(collect(has_constraints)) do i
-        return _allocate_constraints(dest, src, idxmap, i, dest.cone_types[i])
-    end
-    # Load variables
-    _load_variables(dest, length(vis_src))
-    # Set variable attributes
-    MOIU.pass_attributes(dest, src, idxmap, vis_src)
-    # Set model attributes
-    MOIU.pass_attributes(dest, src, idxmap)
-    # Load constraints
-    offset = 0
-    for (i, cache) in zip(has_constraints, caches)
-        _load_constraints(dest, src, idxmap, offset, i, cache)
-        offset += dest.num_rows[i]
+function objective_vector(model::MOI.ModelLike; T = _coef_type(model))
+    obj = MOI.get(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}())
+    c = zeros(MOI.get(model, MOI.NumberOfVariables()))
+    for term in obj.terms
+        c[term.variable.value] += term.coefficient
     end
-    final_touch(dest.A)
-    return idxmap
+    return c
 end

src/product_of_sets.jl

@@ -0,0 +1,75 @@
+# Copyright (c) 2019: Joaquim Dias Garcia, and contributors
+#
+# Use of this source code is governed by an MIT-style license that can be found
+# in the LICENSE.md file or at https://opensource.org/licenses/MIT.
+# Copied from DiffOpt
+
+"""
+    ProductOfSets{T} <: MOI.Utilities.OrderedProductOfSets{T}
+
+The `MOI.Utilities.@product_of_sets` macro requires to know the list of sets
+at compile time. In DiffOpt however, the list depends on what the user is going
+to use as set as DiffOpt supports any set as long as it implements the
+required function of MathOptSetDistances.
+For this type, the list of sets can be given a run-time.
+"""
+mutable struct ProductOfSets{T} <: MOI.Utilities.OrderedProductOfSets{T}
+    """
+    During the copy, this counts the number of rows corresponding to
+    each set. At the end of copy, `final_touch` is called, which
+    converts this list into a cumulative ordering.
+    """
+    num_rows::Vector{Int}
+
+    """
+    A dictionary which maps the `set_index` and `offset` of a set to the
+    dimension, i.e., `dimension[(set_index,offset)] → dim`.
+    """
+    dimension::Dict{Tuple{Int,Int},Int}
+
+    """
+    A sanity bit to check that we don't call functions out-of-order.
+    """
+    final_touch::Bool
+
+    set_types::Vector{Type}
+    set_types_dict::Dict{Type,Int}
+
+    function ProductOfSets{T}() where {T}
+        return new(
+            Int[],
+            Dict{Tuple{Int,Int},Int}(),
+            false,
+            Type[],
+            Dict{Type,Int}(),
+        )
+    end
+end
+
+function MOI.Utilities.set_index(set::ProductOfSets, S::Type{<:MOI.AbstractSet})
+    return get(set.set_types_dict, S, nothing)
+end
+
+MOI.Utilities.set_types(set::ProductOfSets) = set.set_types
+
+function set_set_types(set::ProductOfSets, set_types)
+    resize!(set.num_rows, length(set_types))
+    fill!(set.num_rows, 0)
+    resize!(set.set_types, length(set_types))
+    copy!(set.set_types, set_types)
+    empty!(set.set_types_dict)
+    for i in eachindex(set_types)
+        set.set_types_dict[set_types[i]] = i
+    end
+    return
+end
+
+function add_set_types(set::ProductOfSets, S::Type)
+    if !haskey(set.set_types_dict, S)
+        push!(set.num_rows, 0)
+        push!(set.set_types, S)
+        set.set_types_dict[S] = length(set.set_types)
+        return true
+    end
+    return false
+end