Remove trivial constraints form split_zero and fix index mapping of psd (#3)

Author: blegat

src/Bridges/Constraint/split_zero.jl

@@ -14,17 +14,30 @@ end
 similar_type(::Type{<:MOI.VectorAffineFunction}, T::Type) = MOI.VectorAffineFunction{T}
 
 struct SplitZeroBridge{T, F<:MOI.Utilities.TypedLike{T}, G<:MOI.Utilities.TypedLike{Complex{T}}} <: MOI.Bridges.Constraint.AbstractBridge
-    real::MOI.ConstraintIndex{F, MOI.Zeros}
-    imag::MOI.ConstraintIndex{F, MOI.Zeros}
+    dimension::Int
+    constraint::MOI.ConstraintIndex{F, MOI.Zeros}
+    real_indices::Vector{Int}
+    imag_indices::Vector{Int}
+end
+function _nonzero_indices(func::MOI.AbstractVectorFunction)
+    return [i for (i, scalar_func) in enumerate(MOIU.scalarize(func)) if !iszero(scalar_func)]
 end
 function MOI.Bridges.Constraint.bridge_constraint(
     ::Type{SplitZeroBridge{T, F, G}}, model::MOI.ModelLike,
     f::G,
     set::MOI.Zeros
 ) where {T, F, G}
-    real_con = MOI.add_constraint(model, operate_coefficients(real, T, f), set)
-    imag_con = MOI.add_constraint(model, operate_coefficients(imag, T, f), set)
-    return SplitZeroBridge{T, F, G}(real_con, imag_con)
+    real_part = operate_coefficients(real, T, f)
+    imag_part = operate_coefficients(imag, T, f)
+    real_indices = _nonzero_indices(real_part)
+    imag_indices = _nonzero_indices(imag_part)
+    func = MOIU.operate(
+        vcat, T,
+        MOIU.eachscalar(real_part)[real_indices],
+        MOIU.eachscalar(imag_part)[imag_indices]
+    )
+    constraint = MOI.add_constraint(model, func, MOI.Zeros(length(real_indices) + length(imag_indices)))
+    return SplitZeroBridge{T, F, G}(MOI.dimension(set), constraint, real_indices, imag_indices)
 end
 
 function MOI.supports_constraint(
@@ -45,17 +58,16 @@ end
 # Attributes, Bridge acting as a model
 function MOI.get(::SplitZeroBridge{T, F},
                  ::MOI.NumberOfConstraints{F, MOI.Zeros}) where {T, F}
-    return 2
+    return 1
 end
 function MOI.get(bridge::SplitZeroBridge{T, F},
                  ::MOI.ListOfConstraintIndices{F, MOI.Zeros}) where {T, F}
-    return [bridge.real, bridge.imag]
+    return [bridge.constraint]
 end
 
 # Indices
 function MOI.delete(model::MOI.ModelLike, bridge::SplitZeroBridge)
-    MOI.delete(model, bridge.imag)
-    MOI.delete(model, bridge.real)
+    MOI.delete(model, bridge.constraint)
 end
 
 # Attributes, Bridge acting as a constraint
@@ -68,10 +80,24 @@ function MOI.supports(
 end
 function MOI.get(model::MOI.ModelLike, attr::Union{MOI.ConstraintPrimal, MOI.ConstraintPrimalStart, MOI.ConstraintDual, MOI.ConstraintDualStart},
                  bridge::SplitZeroBridge)
-    return MOI.get(model, attr, bridge.real) + im * MOI.get(model, attr, bridge.imag)
+    values = MOI.get(model, attr, bridge.constraint)
+    output = zeros(Complex{eltype(values)}, bridge.dimension)
+    for (i, idx) in enumerate(bridge.real_indices)
+        output[idx] = values[i]
+    end
+    for (i, idx) in enumerate(bridge.imag_indices)
+        output[idx] = values[length(bridge.real_indices) + i] * im
+    end
+    return output
 end
 function MOI.set(model::MOI.ModelLike, attr::Union{MOI.ConstraintPrimalStart, MOI.ConstraintDualStart},
-                 bridge::SplitZeroBridge, value)
-    MOI.set(model, attr, bridge.real, map(real, value))
-    MOI.set(model, attr, bridge.imag, map(imag, value))
+                 bridge::SplitZeroBridge{T}, value) where T
+    input = Vector{T}(undef, length(bridge.real_indices) + length(bridge.imag_indices))
+    for (i, idx) in enumerate(bridge.real_indices)
+        input[i] = real(value[idx])
+    end
+    for (i, idx) in enumerate(bridge.imag_indices)
+        input[length(bridge.real_indices) + i] = imag(value[idx])
+    end
+    MOI.set(model, attr, bridge.constraint, input)
 end

src/Bridges/Variable/psd.jl

@@ -56,28 +56,30 @@ function MOIB.Variable.bridge_constrained_variable(
     k22 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
     X11() = MOI.SingleVariable(variables[k11])
     X12() = MOI.SingleVariable(variables[k12])
-    X21() = MOI.SingleVariable(variables[k21])
+    function X21(i, j)
+        I = j
+        J = n + i
+        k21 = MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(J - 1)) + I
+        return MOI.SingleVariable(variables[k21])
+    end
     X22() = MOI.SingleVariable(variables[k22])
     con_11_22 = EQ{T}[]
     con12diag = EQ{T}[]
     con_12_21 = EQ{T}[]
     for j in 1:n
-        k21 -= n + 1 - j
         k22 += n
         for i in 1:j
             k11 += 1
             k12 += 1
-            k21 -= 1
             k22 += 1
             push!(con_11_22, MOI.add_constraint(model, MOI.Utilities.operate(-, T, X11(), X22()), MOI.EqualTo(zero(T))))
             if i == j
                 push!(con12diag, MOI.add_constraint(model, convert(MOI.ScalarAffineFunction{T}, X12()), MOI.EqualTo(zero(T))))
             else
-                push!(con_12_21, MOI.add_constraint(model, MOI.Utilities.operate(+, T, X21(), X12()), MOI.EqualTo(zero(T))))
+                push!(con_12_21, MOI.add_constraint(model, MOI.Utilities.operate(+, T, X21(i, j), X12()), MOI.EqualTo(zero(T))))
             end
         end
         k12 += n
-        k21 -= n - j
     end
 
     return HermitianToSymmetricPSDBridge(variables, psd_constraint, con_11_22, con12diag, con_12_21)

test/runtests.jl

@@ -49,7 +49,23 @@ function projection_test(optimizer, config)
     @test MOI.get(optimizer, MOI.ConstraintDual(), cx) ≈ dual atol=atol rtol=rtol
 end
 
-function zero_test(optimizer, config)
+function hermitian_psd_test(optimizer, config)
+    atol = config.atol
+    rtol = config.rtol
+
+    MOI.empty!(optimizer)
+    set = COI.HermitianPositiveSemidefiniteConeTriangle(3)
+    x, cx = MOI.add_constrained_variables(optimizer, set)
+    fx = MOI.SingleVariable.(x)
+    MOI.add_constraints(optimizer, fx, MOI.EqualTo.([1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0]))
+    MOI.optimize!(optimizer)
+    primal = [1.0, 0.0, 1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0]
+    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ primal atol=atol rtol=rtol
+    @test MOI.get(optimizer, MOI.ConstraintPrimal(), cx) ≈ primal atol=atol rtol=rtol
+    @test MOI.get(optimizer, MOI.ConstraintDual(), cx) ≈ zeros(9) atol=atol rtol=rtol
+end
+
+function zero_1_test(optimizer, config)
     atol = config.atol
     rtol = config.rtol
 
@@ -71,14 +87,42 @@ function zero_test(optimizer, config)
     @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ [0.0 + 0.0im] atol=atol rtol=rtol
 end
 
+function zero_2_test(optimizer, config)
+    atol = config.atol
+    rtol = config.rtol
+
+    MOI.empty!(optimizer)
+    x, cx = MOI.add_constrained_variables(optimizer, MOI.Nonnegatives(1))
+    fx = MOI.SingleVariable.(x)
+    func = (1.0 + 0.0im) * fx[1] + 1.0im * fx[1] - 2.0im - (1.0 + 0.0im) * fx[1]
+    c = MOI.add_constraint(
+        optimizer,
+        MOI.Utilities.operate(vcat, Complex{Float64}, func),
+        MOI.Zeros(1)
+    )
+    MOI.optimize!(optimizer)
+    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
+    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ [2.0] atol=atol rtol=rtol
+    @test MOI.get(optimizer, MOI.ConstraintPrimal(), cx) ≈ [2.0] atol=atol rtol=rtol
+    @test MOI.get(optimizer, MOI.ConstraintDual(), cx) ≈ zeros(1) atol=atol rtol=rtol
+    @test MOI.get(optimizer, MOI.ConstraintPrimal(), c) ≈ [0.0 + 0.0im] atol=atol rtol=rtol
+    @test MOI.get(optimizer, MOI.ConstraintDual(), c) ≈ [0.0 + 0.0im] atol=atol rtol=rtol
+end
+
 import CSDP
 @testset "CSDP" begin
     config = MOI.Test.TestConfig(atol=1e-4, rtol=1e-4)
     bridged = MOI.instantiate(CSDP.Optimizer, with_bridge_type=Float64)
     MOI.Bridges.add_bridge(bridged, COI.Bridges.Variable.HermitianToSymmetricPSDBridge{Float64})
     projection_test(bridged, config)
+    hermitian_psd_test(bridged, config)
     bridged = MOI.Bridges.LazyBridgeOptimizer(MOI.Utilities.CachingOptimizer(MOI.Utilities.Model{Float64}(), CSDP.Optimizer()))
     MOI.Bridges.add_bridge(bridged, MOI.Bridges.Constraint.ScalarizeBridge{Float64})
     MOI.Bridges.add_bridge(bridged, COI.Bridges.Constraint.SplitZeroBridge{Float64})
-    zero_test(bridged, config)
+    zero_1_test(bridged, config)
+    cis = MOI.get(bridged.model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}())
+    @test length(cis) == 2
+    zero_2_test(bridged, config)
+    cis = MOI.get(bridged.model, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}())
+    @test length(cis) == 1
 end