Add tests for ScaledDiagonallyDominantBridge (#367)

* Add tests for ScaledDiagonallyDominantBridge

* Add dual

* typo fix

docs/src/tutorials/Polynomial Optimization/bilinear.jl

@@ -2,11 +2,11 @@
 
 #md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Polynomial Optimization/bilinear.ipynb)
 #md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Polynomial Optimization/bilinear.ipynb)
-# **Adapted from**: [Floudas1999; Section 3.1](@cite) and [Lasserre2009; Table 5.1](@cite)
+# **Adapted from**: [Floudas1999; Section 3.2](@cite) and [Lasserre2009; Table 5.1](@cite)
 
 # ## Introduction
 
-# Consider the polynomial optimization problem from [Floudas1999; Section 3.1](@cite).
+# Consider the polynomial optimization problem from [Floudas1999; Section 3.2](@cite).
 
 using Test #src
 using DynamicPolynomials

src/Bridges/Variable/copositive_inner.jl

@@ -95,22 +95,12 @@ end
 
 # TODO ConstraintPrimal, ConstraintDual
 
-# See https://jump.dev/MathOptInterface.jl/v0.9.1/apireference/#MathOptInterface.AbstractSymmetricMatrixSetTriangle
-function matrix_indices(k)
-    j = div(1 + isqrt(8k - 7), 2)
-    i = k - div((j - 1) * j, 2)
-    return i, j
-end
-# Vector index for the vectorization of the triangular part.
-function vector_index(i, j)
-    return div((j - 1) * j, 2) + i
-end
 # Vector index for the vectorization of the off-diagonal triangular part.
 function offdiag_vector_index(i, j)
     if i < j
-        return vector_index(i, j - 1)
+        return MOI.Utilities.trimap(i, j - 1)
     else
-        throw(ArgumentError())
+        throw(ArgumentError("Not off-diagonal"))
     end
 end
 
@@ -121,7 +111,7 @@ function MOI.get(
     i::MOI.Bridges.IndexInVector,
 )
     value = MOI.get(model, attr, bridge.matrix_variables[i.value])
-    row, col = matrix_indices(i.value)
+    row, col = MOI.Utilities.inverse_trimap(i.value)
     if row != col
         value += MOI.get(
             model,
@@ -138,13 +128,13 @@ function MOI.Bridges.bridged_function(
 ) where {T}
     func =
         convert(MOI.ScalarAffineFunction{T}, bridge.matrix_variables[i.value])
-    row, col = matrix_indices(i.value)
+    row, col = MOI.Utilities.inverse_trimap(i.value)
     if row != col
         func = MOI.Utilities.operate!(
             +,
             T,
             func,
-            bridge.nonneg_variables[vector_index(row, col - 1)],
+            bridge.nonneg_variables[MOI.Utilities.trimap(row, col - 1)],
         )
     end
     return func
@@ -157,11 +147,11 @@ function MOI.Bridges.Variable.unbridged_map(
     F = MOI.ScalarAffineFunction{T}
     func = convert(F, vi)
     map = bridge.matrix_variables[i.value] => func
-    row, col = matrix_indices(i.value)
+    row, col = MOI.Utilities.inverse_trimap(i.value)
     if row == col
         return (map,)
     else
-        nneg = bridge.nonneg_variables[vector_index(row, col - 1)]
+        nneg = bridge.nonneg_variables[MOI.Utilities.trimap(row, col - 1)]
         return (map, nneg => zero(F))
     end
 end

src/Bridges/Variable/scaled_diagonally_dominant.jl

@@ -67,7 +67,7 @@ end
 function MOI.Bridges.added_constrained_variable_types(
     ::Type{<:ScaledDiagonallyDominantBridge},
 )
-    return [(SOS.PositiveSemidefinite2x2ConeTriangle,)]
+    return Tuple{Type}[(SOS.PositiveSemidefinite2x2ConeTriangle,)]
 end
 function MOI.Bridges.added_constraint_types(
     ::Type{<:ScaledDiagonallyDominantBridge},
@@ -86,7 +86,7 @@ function MOI.get(
     bridge::ScaledDiagonallyDominantBridge,
     ::MOI.ListOfVariableIndices,
 )
-    return Iterators.flatten(bridge.variables)
+    return collect(Iterators.flatten(bridge.variables))
 end
 function MOI.get(
     bridge::ScaledDiagonallyDominantBridge,
@@ -127,17 +127,88 @@ function MOI.get(
     return SOS.ScaledDiagonallyDominantConeTriangle(bridge.side_dimension)
 end
 
-# TODO ConstraintPrimal, ConstraintDual
+# The map `A` is not injective because it maps the entry of several 2x2 matrices
+# into the same index so it is not invertible hence it's unclear how to implemented
+# `set` for `VariablePrimalStart`.
+# The adjoint `A'` is however injective and it is easy to invert.
 
-trimap(i, j) = div(j * (j - 1), 2) + i
+function MOI.supports(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDualStart,
+    ::Type{<:ScaledDiagonallyDominantBridge},
+)
+    return MOI.supports(
+        model,
+        attr,
+        MOI.ConstraintIndex{
+            MOI.VectorOfVariables,
+            SOS.PositiveSemidefinite2x2ConeTriangle,
+        },
+    )
+end
+
+function MOI.set(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDualStart,
+    bridge::ScaledDiagonallyDominantBridge,
+    value,
+)
+    n = bridge.side_dimension
+    k = 0
+    for j in 1:n
+        for i in 1:(j-1)
+            k += 1
+            # PSD constraints on 2x2 matrices are SOC representable
+            if isnothing(value)
+                dual = nothing
+            else
+                dual = [
+                    value[MOI.Utilities.trimap(i, i)],
+                    value[MOI.Utilities.trimap(i, j)],
+                    value[MOI.Utilities.trimap(j, j)],
+                ]
+            end
+            MOI.set(model, attr, bridge.constraints[k], dual)
+        end
+    end
+    return
+end
+
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
+    bridge::ScaledDiagonallyDominantBridge{T},
+) where {T}
+    n = bridge.side_dimension
+    value = zeros(T, MOI.Utilities.trimap(n, n))
+    k = 0
+    for j in 1:n
+        for i in 1:(j-1)
+            k += 1
+            dual = MOI.get(model, attr, bridge.constraints[k])
+            if isnothing(dual)
+                return nothing
+            end
+            # There are `bridge.side_dimension - 1` possible candidate that should all have
+            # the same `dual` so we take an arbitrary choice
+            if j == i + 1
+                value[MOI.Utilities.trimap(i, i)] = dual[1]
+            elseif i == 1 && j == n
+                value[MOI.Utilities.trimap(j, j)] = dual[3]
+            end
+            value[MOI.Utilities.trimap(i, j)] = dual[2]
+        end
+    end
+    return value
+end
 
 function MOI.get(
     model::MOI.ModelLike,
     attr::MOI.VariablePrimal,
     bridge::ScaledDiagonallyDominantBridge{T},
-    i::MOI.Bridges.IndexInVector,
+    index::MOI.Bridges.IndexInVector,
 ) where {T}
-    i, j = matrix_indices(i.value)
+    i, j = MOI.Utilities.inverse_trimap(index.value)
     if i == j
         value = zero(T)
         for k in 1:(i-1)
@@ -159,7 +230,7 @@ function MOI.Bridges.bridged_function(
     bridge::ScaledDiagonallyDominantBridge{T},
     i::MOI.Bridges.IndexInVector,
 ) where {T}
-    i, j = matrix_indices(i.value)
+    i, j = MOI.Utilities.inverse_trimap(i.value)
     if i == j
         func = zero(MOI.ScalarAffineFunction{T})
         for k in 1:(i-1)
@@ -188,19 +259,28 @@ function MOI.Bridges.Variable.unbridged_map(
     k = 0
     z = zero(SAF)
     saf(i) = convert(SAF, vis[i])
-    # vis[trimap(j, j)] is replaced by a sum of several variables.
+    # vis[MOI.Utilities.trimap(j, j)] is replaced by a sum of several variables.
     # The strategy is to replace all of them by zero except one.
     for j in 1:bridge.side_dimension
         for i in 1:(j-1)
             k += 1
             if i == 1 && j == 2
-                push!(umap, bridge.variables[k][1] => saf(trimap(1, 1)))
+                push!(
+                    umap,
+                    bridge.variables[k][1] => saf(MOI.Utilities.trimap(1, 1)),
+                )
             else
                 push!(umap, bridge.variables[k][1] => z)
             end
-            push!(umap, bridge.variables[k][2] => saf(trimap(i, j)))
+            push!(
+                umap,
+                bridge.variables[k][2] => saf(MOI.Utilities.trimap(i, j)),
+            )
             if i == 1
-                push!(umap, bridge.variables[k][3] => saf(trimap(j, j)))
+                push!(
+                    umap,
+                    bridge.variables[k][3] => saf(MOI.Utilities.trimap(j, j)),
+                )
             else
                 push!(umap, bridge.variables[k][3] => z)
             end

test/Bridges/Variable/scaled_diagonally_dominant.jl

@@ -0,0 +1,65 @@
+module TestVariableScaledDiagonallyDominant
+
+using Test
+import MultivariateBases as MB
+using DynamicPolynomials
+using SumOfSquares
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$(name)", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+function test_error_dim_1()
+    model = MOI.Utilities.Model{Float64}()
+    bridged = MOI.Bridges.Variable.SingleBridgeOptimizer{
+        SumOfSquares.Bridges.Variable.ScaledDiagonallyDominantBridge{Float64},
+    }(
+        model,
+    )
+    err = ErrorException(
+        "The bridges does not work with 1, `matrix_cone` should have returned `Nonnegatives` instead.",
+    )
+    @test_throws err MOI.add_constrained_variables(
+        bridged,
+        SumOfSquares.ScaledDiagonallyDominantConeTriangle(1),
+    )
+end
+
+function test_runtests()
+    @polyvar x y
+    MOI.Bridges.runtests(
+        SumOfSquares.Bridges.Variable.ScaledDiagonallyDominantBridge,
+        model -> begin
+            p, _ = MOI.add_constrained_variables(
+                model,
+                SumOfSquares.ScaledDiagonallyDominantConeTriangle(3),
+            )
+        end,
+        model -> begin
+            q1, _ = MOI.add_constrained_variables(
+                model,
+                SumOfSquares.PositiveSemidefinite2x2ConeTriangle(),
+            )
+            q2, _ = MOI.add_constrained_variables(
+                model,
+                SumOfSquares.PositiveSemidefinite2x2ConeTriangle(),
+            )
+            q3, _ = MOI.add_constrained_variables(
+                model,
+                SumOfSquares.PositiveSemidefinite2x2ConeTriangle(),
+            )
+        end,
+    )
+    return
+end
+
+end  # module
+
+TestVariableScaledDiagonallyDominant.runtests()