Merge pull request #742 from JuliaOpt/bl/rsoc_get

Implement getter of ConstraintFunction/Set for RSOC

Author: blegat

src/Bridges/rsocbridge.jl

@@ -1,5 +1,5 @@
 """
-    RSOCBridge{T}
+    RSOCBridge{T, F, G}
 
 The `RotatedSecondOrderCone` is `SecondOrderCone` representable; see [1, p. 104].
 Indeed, we have ``2tu = (t/√2 + u/√2)^2 - (t/√2 - u/√2)^2`` hence
@@ -16,13 +16,11 @@ That means in particular that the norm is of constraint primal and duals are pre
 
 [1] Ben-Tal, Aharon, and Arkadi Nemirovski. *Lectures on modern convex optimization: analysis, algorithms, and engineering applications*. Society for Industrial and Applied Mathematics, 2001.
 """
-struct RSOCBridge{T, F} <: AbstractBridge
+struct RSOCBridge{T, F, G} <: AbstractBridge
     soc::CI{F, MOI.SecondOrderCone}
 end
-function bridge_constraint(::Type{RSOCBridge{T, F}}, model,
-                           f::MOI.AbstractVectorFunction,
-                           s::MOI.RotatedSecondOrderCone) where {T, F}
-    d = s.dimension
+function rotate_function(f::MOI.AbstractVectorFunction, T::Type)
+    d = MOI.output_dimension(f)
     f_scalars = MOIU.eachscalar(f)
     t = f_scalars[1]
     u = f_scalars[2]
@@ -34,17 +32,22 @@ function bridge_constraint(::Type{RSOCBridge{T, F}}, model,
     # line
     y  = ts - us
     z  = MOIU.operate!(+, T, ts, us)
-    g = MOIU.operate(vcat, T, z, y, x)
-    soc = MOI.add_constraint(model, g, MOI.SecondOrderCone(d))
-    RSOCBridge{T, F}(soc)
+    return MOIU.operate(vcat, T, z, y, x)
+end
+function bridge_constraint(::Type{RSOCBridge{T, F, G}}, model,
+                           f::MOI.AbstractVectorFunction,
+                           s::MOI.RotatedSecondOrderCone) where {T, F, G}
+    soc = MOI.add_constraint(model, rotate_function(f, T),
+                             MOI.SecondOrderCone(MOI.dimension(s)))
+    return RSOCBridge{T, F, G}(soc)
 end
 
 function MOI.supports_constraint(::Type{RSOCBridge{T}},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{MOI.RotatedSecondOrderCone}) where T
     return true
 end
-function added_constraint_types(::Type{RSOCBridge{T, F}}) where {T, F}
+function added_constraint_types(::Type{<:RSOCBridge{T, F}}) where {T, F}
     return [(F, MOI.SecondOrderCone)]
 end
 function concrete_bridge_type(::Type{<:RSOCBridge{T}},
@@ -54,7 +57,7 @@ function concrete_bridge_type(::Type{<:RSOCBridge{T}},
     Y = MOIU.promote_operation(-, T, S, S)
     Z = MOIU.promote_operation(+, T, S, S)
     F = MOIU.promote_operation(vcat, T, Z, Y, G)
-    RSOCBridge{T, F}
+    RSOCBridge{T, F, G}
 end
 
 # Attributes, Bridge acting as an model
@@ -73,6 +76,16 @@ function MOI.delete(model::MOI.ModelLike, c::RSOCBridge)
 end
 
 # Attributes, Bridge acting as a constraint
+function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
+                 bridge::RSOCBridge{T, F, G}) where {T, F, G}
+    # As it is an involution, we can just reapply the same transformation
+    func = MOI.get(model, attr, bridge.soc)
+    return MOIU.convert_approx(G, rotate_function(func, T))
+end
+function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet, bridge::RSOCBridge)
+    set = MOI.get(model, attr, bridge.soc)
+    return MOI.RotatedSecondOrderCone(MOI.dimension(set))
+end
 # As the linear transformation is a symmetric involution,
 # the constraint primal and dual both need to be processed by reapplying the same transformation
 function _get(model, attr::Union{MOI.ConstraintPrimal, MOI.ConstraintDual}, c::RSOCBridge)

src/Test/UnitTests/basic_constraint_tests.jl

@@ -65,10 +65,16 @@ const BasicConstraintTests = Dict(
     (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) => ( dummy_vector_affine, 2, MOI.Nonpositives(2) ),
     (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => ( dummy_vector_affine, 2, MOI.Nonnegatives(2) ),
 
+    (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone)        => ( dummy_vector_affine, 3, MOI.SecondOrderCone(3) ),
+    (MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone) => ( dummy_vector_affine, 3, MOI.RotatedSecondOrderCone(3) ),
+
     (MOI.VectorQuadraticFunction{Float64}, MOI.Reals)        => ( dummy_vector_quadratic, 2, MOI.Reals(2) ),
     (MOI.VectorQuadraticFunction{Float64}, MOI.Zeros)        => ( dummy_vector_quadratic, 2, MOI.Zeros(2) ),
     (MOI.VectorQuadraticFunction{Float64}, MOI.Nonpositives) => ( dummy_vector_quadratic, 2, MOI.Nonpositives(2) ),
-    (MOI.VectorQuadraticFunction{Float64}, MOI.Nonnegatives) => ( dummy_vector_quadratic, 2, MOI.Nonnegatives(2) )
+    (MOI.VectorQuadraticFunction{Float64}, MOI.Nonnegatives) => ( dummy_vector_quadratic, 2, MOI.Nonnegatives(2) ),
+
+    (MOI.VectorQuadraticFunction{Float64}, MOI.SecondOrderCone)        => ( dummy_vector_quadratic, 3, MOI.SecondOrderCone(3) ),
+    (MOI.VectorQuadraticFunction{Float64}, MOI.RotatedSecondOrderCone) => ( dummy_vector_quadratic, 3, MOI.RotatedSecondOrderCone(3) )
 )
 """
     basic_constraint_tests(model::MOI.ModelLike, config::TestConfig;

src/Utilities/functions.jl

@@ -152,7 +152,21 @@ function Base.getindex(it::ScalarFunctionIterator{VAF{T}}, I::AbstractVector) wh
         append!(terms, map(t -> MOI.VectorAffineTerm(i, t), g.terms))
         constant[i] = g.constant
     end
-    VAF(terms, constant)
+    return VAF(terms, constant)
+end
+function Base.getindex(it::ScalarFunctionIterator{VQF{T}}, I::AbstractVector) where T
+    affine_terms = MOI.VectorAffineTerm{T}[]
+    quadratic_terms = MOI.VectorQuadraticTerm{T}[]
+    constant = Vector{T}(undef, length(I))
+    for (i, j) in enumerate(I)
+        g = it[j]
+        append!(affine_terms, map(t -> MOI.VectorAffineTerm(i, t),
+                                  g.affine_terms))
+        append!(quadratic_terms, map(t -> MOI.VectorQuadraticTerm(i, t),
+                                     g.quadratic_terms))
+        constant[i] = g.constant
+    end
+    return VQF(affine_terms, quadratic_terms, constant)
 end
 
 """
@@ -1465,6 +1479,20 @@ function operate(::typeof(vcat), ::Type{T},
     fill_vector(constant, T, 0, 0, fill_constant, output_dim, funcs...)
     return VAF(terms, constant)
 end
+function operate(::typeof(vcat), ::Type{T},
+                 funcs::Union{ScalarQuadraticLike{T}, VVF, VAF{T}, VQF{T}}...) where T
+    num_affine_terms = sum(func -> number_of_affine_terms(T, func), funcs)
+    num_quadratic_terms = sum(func -> number_of_quadratic_terms(T, func), funcs)
+    out_dim = sum(func -> output_dim(T, func), funcs)
+    affine_terms = Vector{MOI.VectorAffineTerm{T}}(undef, num_affine_terms)
+    quadratic_terms = Vector{MOI.VectorQuadraticTerm{T}}(undef, num_quadratic_terms)
+    constant = zeros(T, out_dim)
+    fill_vector(affine_terms, T, 0, 0, fill_terms, number_of_affine_terms, funcs...)
+    fill_vector(quadratic_terms, T, 0, 0, fill_terms, number_of_quadratic_terms, funcs...)
+    fill_vector(constant, T, 0, 0, fill_constant, output_dim, funcs...)
+    return VQF(affine_terms, quadratic_terms, constant)
+end
+
 
 # Similar to `eachscalar` but faster, see
 # https://github.com/JuliaOpt/MathOptInterface.jl/issues/418
@@ -1518,3 +1546,20 @@ function count_terms(dimension::I, terms::Vector{T}) where {I,T}
     count_terms(counting, terms)
     return counting
 end
+
+convert_approx(::Type{T}, func::T; kws...) where {T} = func
+function convert_approx(::Type{MOI.SingleVariable}, func::MOI.ScalarAffineFunction{T};
+                        tol=sqrt(eps(T))) where {T}
+    f = MOIU.canonical(func)
+    i = findfirst(t -> isapprox(t.coefficient, one(T), atol=tol), f.terms)
+    if abs(f.constant) > tol || i === nothing ||
+        any(j -> j != i && abs(f.terms[i]) > tol, eachindex(f.terms))
+        throw(InexactError(:convert_approx, MOI.SingleVariable, func))
+    end
+    return MOI.SingleVariable(f.terms[i].variable_index)
+end
+function convert_approx(::Type{MOI.VectorOfVariables}, func::MOI.VectorAffineFunction{T};
+    tol=sqrt(eps(T))) where {T}
+    return MOI.VectorOfVariables([convert_approx(MOI.SingleVariable, f, tol=tol).variable
+                                  for f in scalarize(func)])
+end

test/Bridges/rsocbridge.jl

@@ -1,5 +1,28 @@
+using Test
+
+using MathOptInterface
+const MOI = MathOptInterface
+const MOIT = MathOptInterface.Test
+const MOIU = MathOptInterface.Utilities
+const MOIB = MathOptInterface.Bridges
+
+include("utilities.jl")
+
+include("simple_model.jl")
+
+mock = MOIU.MockOptimizer(SimpleModel{Float64}())
+config = MOIT.TestConfig()
+
 @testset "RSOC" begin
     bridged_mock = MOIB.RSOC{Float64}(mock)
+
+    MOIT.basic_constraint_tests(
+        bridged_mock, config,
+        include = [(F, S)
+                   for F in [MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64},
+                             MOI.VectorQuadraticFunction{Float64}]
+                   for S in [MOI.RotatedSecondOrderCone]])
+
     mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1/√2, 1/√2, 0.5, 1.0],
                           (MOI.SingleVariable,                MOI.EqualTo{Float64}) => [-√2, -1/√2],
                           (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone)  => [[3/2, 1/2, -1.0, -1.0]])