diff --git a/src/Bridges/Constraint/det.jl b/src/Bridges/Constraint/det.jl index b0a5521391..40c4f7184e 100644 --- a/src/Bridges/Constraint/det.jl +++ b/src/Bridges/Constraint/det.jl @@ -1,12 +1,3 @@ -# TODO(odow): this appears elsewhere -function trimap(i::Integer, j::Integer) - if i < j - trimap(j, i) - else - div((i - 1) * i, 2) + j - end -end - """ extract_eigenvalues( model, @@ -32,14 +23,14 @@ function extract_eigenvalues( f_scalars = MOIU.eachscalar(f) tu = [f_scalars[i] for i in 1:offset] - n = trimap(d, d) + n = MOIU.trimap(d, d) X = f_scalars[offset.+(1:n)] m = length(X.terms) M = m + n + d terms = Vector{MOI.VectorAffineTerm{T}}(undef, M) terms[1:m] = X.terms - N = trimap(2d, 2d) + N = MOIU.trimap(2d, 2d) constant = zeros(T, N) constant[1:n] = X.constants @@ -50,14 +41,14 @@ function extract_eigenvalues( for i in j:d cur += 1 terms[cur] = MOI.VectorAffineTerm( - trimap(i, d + j), - MOI.ScalarAffineTerm(one(T), Δ[trimap(i, j)]), + MOIU.trimap(i, d + j), + MOI.ScalarAffineTerm(one(T), Δ[MOIU.trimap(i, j)]), ) end cur += 1 terms[cur] = MOI.VectorAffineTerm( - trimap(d + j, d + j), - MOI.ScalarAffineTerm(one(T), Δ[trimap(j, j)]), + MOIU.trimap(d + j, d + j), + MOI.ScalarAffineTerm(one(T), Δ[MOIU.trimap(j, j)]), ) end @assert cur == M @@ -65,7 +56,7 @@ function extract_eigenvalues( sdindex = MOI.add_constraint(model, Y, MOI.PositiveSemidefiniteConeTriangle(2d)) - D = Δ[trimap.(1:d, 1:d)] + D = Δ[MOIU.trimap.(1:d, 1:d)] return tu, D, Δ, sdindex end diff --git a/src/Bridges/Constraint/norm_spec_nuc_to_psd.jl b/src/Bridges/Constraint/norm_spec_nuc_to_psd.jl index be3a8764fd..d91550abda 100644 --- a/src/Bridges/Constraint/norm_spec_nuc_to_psd.jl +++ b/src/Bridges/Constraint/norm_spec_nuc_to_psd.jl @@ -24,14 +24,14 @@ function bridge_constraint( psd_set = MOI.PositiveSemidefiniteConeTriangle(side_dim) psd_func = MOIU.zero_with_output_dimension(F, MOI.dimension(psd_set)) for i in 1:side_dim - MOIU.operate_output_index!(+, T, trimap(i, i), psd_func, t) + MOIU.operate_output_index!(+, T, MOIU.trimap(i, i), psd_func, t) end X_idx = 2 for j in 1:column_dim, i in (column_dim+1):side_dim MOIU.operate_output_index!( +, T, - trimap(i, j), + MOIU.trimap(i, j), psd_func, f_scalars[X_idx], ) @@ -101,7 +101,7 @@ function MOI.get( t = psd_func[1] side_dim = bridge.row_dim + bridge.column_dim X = psd_func[[ - trimap(i, j) for j in 1:bridge.column_dim for + MOIU.trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim+1):side_dim ]] return MOIU.convert_approx(G, MOIU.operate(vcat, T, t, X)) @@ -132,7 +132,7 @@ function MOI.get( t = primal[1] side_dim = bridge.row_dim + bridge.column_dim X = primal[[ - trimap(i, j) for j in 1:bridge.column_dim for + MOIU.trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim+1):side_dim ]] return vcat(t, X) @@ -148,11 +148,11 @@ function MOI.set( side_dim = bridge.row_dim + column_dim primal = zeros(T, div(side_dim * (side_dim + 1), 2)) for i in 1:side_dim - primal[trimap(i, i)] = value[1] + primal[MOIU.trimap(i, i)] = value[1] end X_idx = 2 for j in 1:column_dim, i in (column_dim+1):side_dim - primal[trimap(i, j)] = value[X_idx] + primal[MOIU.trimap(i, j)] = value[X_idx] X_idx += 1 end MOI.set(model, MOI.ConstraintPrimalStart(), bridge.psd_index, primal) @@ -168,10 +168,11 @@ function MOI.get( dual = MOI.get(model, MOI.ConstraintDual(), bridge.psd_index) column_dim = bridge.column_dim side_dim = bridge.row_dim + column_dim - t = sum(dual[trimap(i, i)] for i in 1:side_dim) + t = sum(dual[MOIU.trimap(i, i)] for i in 1:side_dim) X = 2 * dual[[ - trimap(i, j) for j in 1:column_dim for i in (column_dim+1):side_dim + MOIU.trimap(i, j) for j in 1:column_dim for + i in (column_dim+1):side_dim ]] return vcat(t, X) end @@ -207,8 +208,8 @@ function bridge_constraint( U = MOI.add_variables(model, U_dim) V = MOI.add_variables(model, V_dim) diag_vars = vcat( - [U[trimap(i, i)] for i in 1:column_dim], - [V[trimap(i, i)] for i in 1:row_dim], + [U[MOIU.trimap(i, i)] for i in 1:column_dim], + [V[MOIU.trimap(i, i)] for i in 1:row_dim], ) rhs = MOI.ScalarAffineFunction( MOI.ScalarAffineTerm.(one(T), diag_vars), @@ -230,7 +231,7 @@ function bridge_constraint( vcat, T, psd_func, - MOI.VectorOfVariables(V[[trimap(i, j) for j in 1:i]]), + MOI.VectorOfVariables(V[[MOIU.trimap(i, j) for j in 1:i]]), ) end psd_index = MOI.add_constraint(model, psd_func, psd_set) @@ -346,13 +347,17 @@ function MOI.get( MOIU.operate( /, T, - MOIU.operate(+, T, [psd_func[trimap(i, i)] for i in 1:side_dim]...), + MOIU.operate( + +, + T, + [psd_func[MOIU.trimap(i, i)] for i in 1:side_dim]..., + ), T(2), ), ) t = MOIU.remove_variable(MOIU.remove_variable(t, bridge.U), bridge.V) X = psd_func[[ - trimap(i, j) for j in 1:bridge.column_dim for + MOIU.trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim+1):side_dim ]] return MOIU.convert_approx(H, MOIU.operate(vcat, T, t, X)) @@ -382,9 +387,9 @@ function MOI.get( ge_primal = MOI.get(model, MOI.ConstraintPrimal(), bridge.ge_index) psd_primal = MOI.get(model, MOI.ConstraintPrimal(), bridge.psd_index) side_dim = bridge.row_dim + bridge.column_dim - t = ge_primal + sum(psd_primal[trimap(i, i)] for i in 1:side_dim) / 2 + t = ge_primal + sum(psd_primal[MOIU.trimap(i, i)] for i in 1:side_dim) / 2 X = psd_primal[[ - trimap(i, j) for j in 1:bridge.column_dim for + MOIU.trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim+1):side_dim ]] return vcat(t, X) @@ -401,7 +406,7 @@ function MOI.get( side_dim = bridge.row_dim + bridge.column_dim X = 2 * psd_dual[[ - trimap(i, j) for j in 1:bridge.column_dim for + MOIU.trimap(i, j) for j in 1:bridge.column_dim for i in (bridge.column_dim+1):side_dim ]] return vcat(t, X) @@ -418,11 +423,11 @@ function MOI.set( side_dim = bridge.row_dim + column_dim dual = zeros(T, div(side_dim * (side_dim + 1), 2)) for i in 1:side_dim - dual[trimap(i, i)] = value[1] + dual[MOIU.trimap(i, i)] = value[1] end X_idx = 2 for j in 1:column_dim, i in (column_dim+1):side_dim - dual[trimap(i, j)] = value[X_idx] / 2 + dual[MOIU.trimap(i, j)] = value[X_idx] / 2 X_idx += 1 end MOI.set(model, MOI.ConstraintDualStart(), bridge.psd_index, dual) diff --git a/src/Bridges/Constraint/soc_to_psd.jl b/src/Bridges/Constraint/soc_to_psd.jl index ffd6de6036..3c5ef137b9 100644 --- a/src/Bridges/Constraint/soc_to_psd.jl +++ b/src/Bridges/Constraint/soc_to_psd.jl @@ -15,10 +15,10 @@ function _SOCtoPSDaff( dim = length(f_scalars) n = div(dim * (dim + 1), 2) h = MOIU.zero_with_output_dimension(F, n) - MOIU.operate_output_index!(+, T, trimap(1, 1), h, f_scalars[1]) + MOIU.operate_output_index!(+, T, MOIU.trimap(1, 1), h, f_scalars[1]) for i in 2:dim - MOIU.operate_output_index!(+, T, trimap(1, i), h, f_scalars[i]) - MOIU.operate_output_index!(+, T, trimap(i, i), h, g) + MOIU.operate_output_index!(+, T, MOIU.trimap(1, i), h, f_scalars[i]) + MOIU.operate_output_index!(+, T, MOIU.trimap(i, i), h, g) end return h end @@ -89,14 +89,14 @@ end function MOIB.inverse_map_function(::Type{<:SOCtoPSDBridge}, func) scalars = MOIU.eachscalar(func) dim = MOIU.side_dimension_for_vectorized_dimension(length(scalars)) - return scalars[trimap.(1, 1:dim)] + return scalars[MOIU.trimap.(1, 1:dim)] end function MOIB.adjoint_map_function(::Type{<:SOCtoPSDBridge{T}}, func) where {T} scalars = MOIU.eachscalar(func) dim = MOIU.side_dimension_for_vectorized_dimension(length(scalars)) - tdual = sum(i -> func[trimap(i, i)], 1:dim) - return MOIU.operate(vcat, T, tdual, func[trimap.(2:dim, 1)] * 2) + tdual = sum(i -> func[MOIU.trimap(i, i)], 1:dim) + return MOIU.operate(vcat, T, tdual, func[MOIU.trimap.(2:dim, 1)] * 2) end function MOIB.inverse_adjoint_map_function( @@ -189,14 +189,26 @@ function MOIB.inverse_map_function(::Type{<:RSOCtoPSDBridge{T}}, func) where {T} t = scalars[1] # It is (2u*I)[1,1] so it needs to be divided by 2 to get u u = MOIU.operate!(/, T, scalars[3], convert(T, 2)) - return MOIU.operate(vcat, T, t, u, scalars[[trimap(1, i) for i in 2:dim]]) + return MOIU.operate( + vcat, + T, + t, + u, + scalars[[MOIU.trimap(1, i) for i in 2:dim]], + ) end function MOIB.adjoint_map_function(::Type{<:RSOCtoPSDBridge{T}}, func) where {T} scalars = MOIU.eachscalar(func) dim = MOIU.side_dimension_for_vectorized_dimension(length(scalars)) - udual = sum(i -> func[trimap(i, i)], 2:dim) - return MOIU.operate(vcat, T, func[1], 2udual, func[trimap.(2:dim, 1)] * 2) + udual = sum(i -> func[MOIU.trimap(i, i)], 2:dim) + return MOIU.operate( + vcat, + T, + func[1], + 2udual, + func[MOIU.trimap.(2:dim, 1)] * 2, + ) end function MOIB.inverse_adjoint_map_function( diff --git a/src/Bridges/Variable/rsoc_to_psd.jl b/src/Bridges/Variable/rsoc_to_psd.jl index cc1646cab8..cd774451ba 100644 --- a/src/Bridges/Variable/rsoc_to_psd.jl +++ b/src/Bridges/Variable/rsoc_to_psd.jl @@ -62,7 +62,7 @@ function bridge_constrained_variable( func = MOIU.operate(-, T, u2, MOI.SingleVariable(variables[k])) push!(diag, MOI.add_constraint(model, func, MOI.EqualTo(zero(T)))) end - @assert k == trimap(dim, dim) + @assert k == MOIU.trimap(dim, dim) end return RSOCtoPSDBridge{T}(variables, psd, off_diag, diag) end @@ -168,14 +168,6 @@ function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::RSOCtoPSDBridge) return MOI.RotatedSecondOrderCone(dim) end -function trimap(i::Integer, j::Integer) - if i < j - return trimap(j, i) - else - return div((i - 1) * i, 2) + j - end -end - function _variable_map(bridge::RSOCtoPSDBridge, i::MOIB.IndexInVector) if bridge.psd isa MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives} @@ -185,7 +177,7 @@ function _variable_map(bridge::RSOCtoPSDBridge, i::MOIB.IndexInVector) elseif i.value == 2 return 3 else - return trimap(1, i.value - 1) + return MOIU.trimap(1, i.value - 1) end end diff --git a/src/Utilities/sets.jl b/src/Utilities/sets.jl index faf225c619..9577c1ff0e 100644 --- a/src/Utilities/sets.jl +++ b/src/Utilities/sets.jl @@ -122,3 +122,19 @@ Return the dimension `d` such that function side_dimension_for_vectorized_dimension(n::Base.Integer) return div(isqrt(1 + 8n), 2) end + +""" + trimap(row::Integer, column::Integer) + +Convert between the row and column indices of a matrix, to the linear index of +the corresponding element in the triangular representation. + +This is most useful when mapping between `ConeSquare` and `ConeTriangle` sets, +e.g., as part of an [`MOI.AbstractSymmetricMatrixSetTriangle`](@ref) set. +""" +function trimap(row::Integer, column::Integer) + if row < column + return trimap(column, row) + end + return div((row - 1) * row, 2) + column +end diff --git a/test/Utilities/sets.jl b/test/Utilities/sets.jl index c9a925ec80..b8a93f5af3 100644 --- a/test/Utilities/sets.jl +++ b/test/Utilities/sets.jl @@ -175,6 +175,19 @@ function test_dot_coefficients() @test MOIU.dot_coefficients(sp_vec, MOI.LogDetConeTriangle(3)) == sp_vec end +function test_trimap() + @test MOIU.trimap(1, 1) == 1 + @test MOIU.trimap(1, 2) == 2 + @test MOIU.trimap(2, 1) == 2 + @test MOIU.trimap(2, 2) == 3 + @test MOIU.trimap(3, 1) == 4 + @test MOIU.trimap(1, 3) == 4 + @test MOIU.trimap(3, 2) == 5 + @test MOIU.trimap(2, 3) == 5 + @test MOIU.trimap(3, 3) == 6 + return +end + end # module TestSets.runtests()