diff --git a/Project.toml b/Project.toml index dcf2bea..6a023ad 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "2.1.0" [deps] MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" SCS_jll = "f4f2fc5b-1d94-523c-97ea-2ab488bedf4b" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -19,6 +20,7 @@ SCSSCS_MKL_jllExt = ["SCS_MKL_jll"] [compat] MathOptInterface = "1.20" +MutableArithmetics = "1" Pkg = "<0.0.1, ^1.6" PrecompileTools = "1" SCS_GPU_jll = "=3.2.8" diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl index 18e651c..a1e7671 100644 --- a/src/MOI_wrapper/MOI_wrapper.jl +++ b/src/MOI_wrapper/MOI_wrapper.jl @@ -4,6 +4,8 @@ # in the LICENSE.md file or at https://opensource.org/licenses/MIT. include("scaled_psd_cone_bridge.jl") +include("hermitian_complex_psd_cone_bridge.jl") +include("scaled_complex_psd_cone_bridge.jl") MOI.Utilities.@product_of_sets( _Cones, @@ -11,6 +13,7 @@ MOI.Utilities.@product_of_sets( MOI.Nonnegatives, MOI.SecondOrderCone, ScaledPSDCone, + ScaledComplexPSDCone, MOI.ExponentialCone, MOI.DualExponentialCone, MOI.PowerCone{T}, @@ -132,7 +135,11 @@ mutable struct Optimizer <: MOI.AbstractOptimizer end function MOI.get(::Optimizer, ::MOI.Bridges.ListOfNonstandardBridges) - return Type[ScaledPSDConeBridge{Cdouble}] + return Type[ + ScaledPSDConeBridge{Cdouble}, + ScaledComplexPSDConeBridge{Cdouble}, + HermitianComplexPSDConeBridge{Cdouble}, + ] end MOI.get(::Optimizer, ::MOI.SolverName) = "SCS" @@ -239,6 +246,7 @@ function MOI.supports_constraint( MOI.Nonnegatives, MOI.SecondOrderCone, ScaledPSDCone, + ScaledComplexPSDCone, MOI.ExponentialCone, MOI.DualExponentialCone, MOI.PowerCone{Cdouble}, @@ -390,9 +398,9 @@ function MOI.optimize!( Float64[], # # placeholder: bl _map_sets(MOI.dimension, T, Ab, MOI.SecondOrderCone), _map_sets(MOI.side_dimension, T, Ab, ScaledPSDCone), - Int[], # placeholder complex PSD - div(Ab.sets.num_rows[5] - Ab.sets.num_rows[4], 3), + _map_sets(MOI.side_dimension, T, Ab, ScaledComplexPSDCone), div(Ab.sets.num_rows[6] - Ab.sets.num_rows[5], 3), + div(Ab.sets.num_rows[7] - Ab.sets.num_rows[6], 3), vcat( _map_sets(set -> set.exponent, Float64, Ab, MOI.PowerCone{Cdouble}), _map_sets( diff --git a/src/MOI_wrapper/hermitian_complex_psd_cone_bridge.jl b/src/MOI_wrapper/hermitian_complex_psd_cone_bridge.jl new file mode 100644 index 0000000..7aa207b --- /dev/null +++ b/src/MOI_wrapper/hermitian_complex_psd_cone_bridge.jl @@ -0,0 +1,192 @@ +# Copyright (c) 2014: SCS.jl 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. + +""" + struct ComplexPositiveSemidefiniteConeTriangle <: MOI.AbstractVectorSet + side_dimension::Int + end + +Similar to `HermitianPositiveSemidefiniteConeTriangle` but its +vectorization interleaves real and imaginary parts of the off-diagonal. +""" +struct ComplexPositiveSemidefiniteConeTriangle <: MOI.AbstractVectorSet + side_dimension::Int +end + +function MOI.Utilities.set_with_dimension( + ::Type{ComplexPositiveSemidefiniteConeTriangle}, + dim, +) + return ComplexPositiveSemidefiniteConeTriangle(isqrt(dim)) +end + +function Base.copy(x::ComplexPositiveSemidefiniteConeTriangle) + return ComplexPositiveSemidefiniteConeTriangle(x.side_dimension) +end + +function MOI.side_dimension(x::ComplexPositiveSemidefiniteConeTriangle) + return x.side_dimension +end + +function MOI.dimension(x::ComplexPositiveSemidefiniteConeTriangle) + return x.side_dimension^2 +end + +import MutableArithmetics as MA + +function MOI.Utilities.set_dot( + x::AbstractVector{S}, + y::AbstractVector{T}, + set::ComplexPositiveSemidefiniteConeTriangle, +) where {S,T} + U = MA.promote_operation(MA.add_mul, S, T) + result = zero(U) + d = set.side_dimension + k = 0 + for j in 1:d + for i in 1:j-1 + k += 1 + result = MA.add_mul!!(result, 2, x[k], y[k]) + k += 1 + result = MA.add_mul!!(result, 2, x[k], y[k]) + end + k += 1 + result = MA.add_mul!!(result, x[k], y[k]) + end + return result +end + +function MOI.Utilities.dot_coefficients( + a::AbstractVector, + set::ComplexPositiveSemidefiniteConeTriangle, +) + d = set.side_dimension + b = copy(a) + k = 0 + for j in 1:d + for i in 1:j-1 + k += 1 + b[k] /= 2 + k += 1 + b[k] /= 2 + end + k += 1 + end + return b +end + +function MOI.is_set_dot_scaled(::Type{ComplexPositiveSemidefiniteConeTriangle}) + return true +end + +struct HermitianComplexPSDConeBridge{T,F} <: + MOI.Bridges.Constraint.SetMapBridge{ + T, + ComplexPositiveSemidefiniteConeTriangle, + MOI.HermitianPositiveSemidefiniteConeTriangle, + F, + F, +} + constraint::MOI.ConstraintIndex{F,ComplexPositiveSemidefiniteConeTriangle} +end + +MOI.Bridges.bridging_cost(::Type{<:HermitianComplexPSDConeBridge}) = 0.1 + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{HermitianComplexPSDConeBridge{T}}, + ::Type{F}, + ::Type{MOI.HermitianPositiveSemidefiniteConeTriangle}, +) where {T,F<:MOI.AbstractVectorFunction} + return HermitianComplexPSDConeBridge{T,F} +end + +function MOI.Bridges.map_set( + ::Type{<:HermitianComplexPSDConeBridge}, + set::MOI.HermitianPositiveSemidefiniteConeTriangle, +) + return ComplexPositiveSemidefiniteConeTriangle(MOI.side_dimension(set)) +end + +function MOI.Bridges.inverse_map_set( + ::Type{<:HermitianComplexPSDConeBridge}, + set::ComplexPositiveSemidefiniteConeTriangle, +) + return MOI.HermitianPositiveSemidefiniteConeTriangle(set.side_dimension) +end + +function _hermitian_to_complex(func) + vals = MOI.Utilities.eachscalar(func) + dim = length(vals) + side = isqrt(dim) + k_re = 1 + k_im = div(side * (side + 1), 2) + 1 + l = 1 + perm = zeros(Int, dim) + for i in 1:side + for j in 1:(i-1) + perm[l] = k_re + perm[l+1] = k_im + k_re += 1 + k_im += 1 + l += 2 + end + perm[l] = k_re + k_re += 1 + l += 1 + end + return vals[perm] +end + +function _complex_to_hermitian(func) + vals = MOI.Utilities.eachscalar(func) + dim = length(vals) + side = isqrt(dim) + k_re = 1 + k_im = div(side * (side + 1), 2) + 1 + l = 1 + perm = zeros(Int, dim) + for i in 1:side + for j in 1:(i-1) + perm[k_re] = l + perm[k_im] = l + 1 + k_re += 1 + k_im += 1 + l += 2 + end + perm[k_re] = l + k_re += 1 + l += 1 + end + return vals[perm] +end + +# Map ConstraintFunction from Hermitian -> Complex +function MOI.Bridges.map_function(::Type{<:HermitianComplexPSDConeBridge}, f) + return _hermitian_to_complex(f) +end + +# Used to map the ConstraintPrimal from Complex -> Hermitian +function MOI.Bridges.inverse_map_function( + ::Type{<:HermitianComplexPSDConeBridge}, + f, +) + return _complex_to_hermitian(f) +end + +# Used to map the ConstraintDual from Complex -> Hermitian +function MOI.Bridges.adjoint_map_function( + ::Type{<:HermitianComplexPSDConeBridge}, + f, +) + return _complex_to_hermitian(f) +end + +# Used to set ConstraintDualStart +function MOI.Bridges.inverse_adjoint_map_function( + ::Type{<:HermitianComplexPSDConeBridge}, + f, +) + return _hermitian_to_complex(f) +end diff --git a/src/MOI_wrapper/scaled_complex_psd_cone_bridge.jl b/src/MOI_wrapper/scaled_complex_psd_cone_bridge.jl new file mode 100644 index 0000000..8e60ddd --- /dev/null +++ b/src/MOI_wrapper/scaled_complex_psd_cone_bridge.jl @@ -0,0 +1,131 @@ +# Copyright (c) 2014: SCS.jl 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. + +""" + struct ScaledComplexPSDCone <: MOI.AbstractVectorSet + side_dimension::Int + end + +Similar to `MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle}` but its +vectorization is the lower triangular part column-wise. +""" +struct ScaledComplexPSDCone <: MOI.AbstractVectorSet + side_dimension::Int +end + +function MOI.Utilities.set_with_dimension(::Type{ScaledComplexPSDCone}, dim) + return ScaledComplexPSDCone(isqrt(dim)) +end + +Base.copy(x::ScaledComplexPSDCone) = ScaledComplexPSDCone(x.side_dimension) + +MOI.side_dimension(x::ScaledComplexPSDCone) = x.side_dimension + +function MOI.dimension(x::ScaledComplexPSDCone) + return x.side_dimension^2 +end + +struct ScaledComplexPSDConeBridge{T,F} <: MOI.Bridges.Constraint.SetMapBridge{ + T, + ScaledComplexPSDCone, + MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle}, + F, + F, +} + constraint::MOI.ConstraintIndex{F,ScaledComplexPSDCone} +end + +MOI.Bridges.bridging_cost(::Type{<:ScaledComplexPSDConeBridge}) = 0.1 + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{ScaledComplexPSDConeBridge{T}}, + ::Type{F}, + ::Type{MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle}}, +) where {T,F<:MOI.AbstractVectorFunction} + return ScaledComplexPSDConeBridge{T,F} +end + +function MOI.Bridges.map_set( + ::Type{<:ScaledComplexPSDConeBridge}, + set::MOI.Scaled{ComplexPositiveSemidefiniteConeTriangle}, +) + return ScaledComplexPSDCone(MOI.side_dimension(set)) +end + +function MOI.Bridges.inverse_map_set( + ::Type{<:ScaledComplexPSDConeBridge}, + set::ScaledComplexPSDCone, +) + return MOI.Scaled( + ComplexPositiveSemidefiniteConeTriangle(set.side_dimension), + ) +end + +function _complex_to_scs(func) + vals = MOI.Utilities.eachscalar(func) + d = isqrt(length(vals)) + c = 0 + perm = zeros(Int, length(vals)) + for i in 1:d + c += 1 + perm[c] = i^2 + for j in i+1:d + triidx = 2i - 1 + (j - 1)^2 + c += 1 + perm[c] = triidx + c += 1 + perm[c] = triidx + 1 + end + end + return vals[perm] +end + +function _scs_to_complex(func) + vals = MOI.Utilities.eachscalar(func) + d = isqrt(length(vals)) + c = 0 + perm = zeros(Int, length(vals)) + for i in 1:d + c += 1 + perm[i^2] = c + for j in i+1:d + triidx = 2i - 1 + (j - 1)^2 + c += 1 + perm[triidx] = c + c += 1 + perm[triidx+1] = c + end + end + return vals[perm] +end + +# Map ConstraintFunction from MOI -> SCS +function MOI.Bridges.map_function(::Type{<:ScaledComplexPSDConeBridge}, f) + return _complex_to_scs(f) +end + +# Used to map the ConstraintPrimal from SCS -> MOI +function MOI.Bridges.inverse_map_function( + ::Type{<:ScaledComplexPSDConeBridge}, + f, +) + return _scs_to_complex(f) +end + +# Used to map the ConstraintDual from SCS -> MOI +function MOI.Bridges.adjoint_map_function( + ::Type{<:ScaledComplexPSDConeBridge}, + f, +) + return _scs_to_complex(f) +end + +# Used to set ConstraintDualStart +function MOI.Bridges.inverse_adjoint_map_function( + ::Type{<:ScaledComplexPSDConeBridge}, + f, +) + return _complex_to_scs(f) +end diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index bbbc765..eec18fc 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -61,7 +61,6 @@ function _test_runtests(linear_solver) # Unexpected failures: # TODO(odow): looks like a tolerance issue? "test_linear_add_constraints", - "test_conic_HermitianPositiveSemidefiniteConeTriangle_2", # Expected test failures: # TODO(odow): get not supported for primal/dual starts "test_model_ModelFilter_AbstractConstraintAttribute",