From 9de6dacccd6d2e60d4470cd80dd9780e06cb36d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 30 Jun 2023 21:18:01 +0200 Subject: [PATCH] Fix automatic adding complex bridge (#313) * Fix automatic adding complex bridge * Fix format * Fix doc ref --- docs/src/constraints.md | 2 +- .../sums_of_hermitian_squares.jl | 4 +- docs/src/tutorials/Symmetry/cyclic.jl | 5 +- src/constraint.jl | 60 +++++++++---------- src/variable.jl | 6 +- 5 files changed, 34 insertions(+), 43 deletions(-) diff --git a/docs/src/constraints.md b/docs/src/constraints.md index 621af67bc..19eeec63f 100644 --- a/docs/src/constraints.md +++ b/docs/src/constraints.md @@ -206,7 +206,7 @@ So `p(x)` is constrained to be equal to `s(x) = s_0(x) + s_1(x) * q_1(x) + s_2(x) * q_2(x) + ...` where the `s_i(x)` polynomials are Sum-of-Squares. The dual of the equality constraint between `p(x)` and `s(x)` is given -by [`SumOfSquares.PolyJuMP.moments`](@ref). +by [`SumOfSquares.MultivariateMoments.moments`](@ref). ```julia μ = moments(cref) ``` diff --git a/docs/src/tutorials/Noncommutative and Hermitian/sums_of_hermitian_squares.jl b/docs/src/tutorials/Noncommutative and Hermitian/sums_of_hermitian_squares.jl index fd7d544a9..f829f3c80 100644 --- a/docs/src/tutorials/Noncommutative and Hermitian/sums_of_hermitian_squares.jl +++ b/docs/src/tutorials/Noncommutative and Hermitian/sums_of_hermitian_squares.jl @@ -10,12 +10,10 @@ using SumOfSquares using DynamicPolynomials @ncpolyvar x y -p = (x + 1.0im * y) * (x - im * y) +p = (x + im * y) * (x - im * y) import CSDP model = Model(CSDP.Optimizer) -MOI.Bridges.add_bridge(backend(model).optimizer, PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Complex{Float64}}) -MOI.Bridges.add_bridge(backend(model).optimizer, SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Complex{Float64}}) cone = NonnegPolyInnerCone{MOI.HermitianPositiveSemidefiniteConeTriangle}() con_ref = @constraint(model, p in cone) optimize!(model) diff --git a/docs/src/tutorials/Symmetry/cyclic.jl b/docs/src/tutorials/Symmetry/cyclic.jl index def2b773f..f7088f4b0 100644 --- a/docs/src/tutorials/Symmetry/cyclic.jl +++ b/docs/src/tutorials/Symmetry/cyclic.jl @@ -140,14 +140,11 @@ b = √3/2 import CSDP solver = CSDP.Optimizer model = Model(solver) -MOI.Bridges.add_bridge(backend(model).optimizer, PolyJuMP.Bridges.Constraint.ZeroPolynomialBridge{Complex{Float64}}) -MOI.Bridges.add_bridge(backend(model).optimizer, SumOfSquares.Bridges.Constraint.SOSPolynomialBridge{Complex{Float64}}) @variable(model, t) @objective(model, Max, t) pattern = Symmetry.Pattern(G, action) cone = SumOfSquares.NonnegPolyInnerCone{MOI.HermitianPositiveSemidefiniteConeTriangle}() -pp = (1.0 + 0.0im) * poly - (1.0 + 0.0im) * t -con_ref = @constraint(model, pp in cone, symmetry = pattern) +con_ref = @constraint(model, poly - t in cone, symmetry = pattern) optimize!(model) @test termination_status(model) == MOI.OPTIMAL #src @test objective_value(model) ≈ 0.0 atol=1e-4 #src diff --git a/src/constraint.jl b/src/constraint.jl index 9b5c32d64..21a1e5413 100644 --- a/src/constraint.jl +++ b/src/constraint.jl @@ -290,42 +290,54 @@ function PolyJuMP.bridges( ::Type{<:MOI.AbstractVectorFunction}, ::Type{EmptyCone}, ) - return [Bridges.Constraint.EmptyBridge] + return [(Bridges.Constraint.EmptyBridge, Float64)] end function PolyJuMP.bridges( ::Type{<:MOI.AbstractVectorFunction}, ::Type{PositiveSemidefinite2x2ConeTriangle}, ) - return [Bridges.Constraint.PositiveSemidefinite2x2Bridge] + return [(Bridges.Constraint.PositiveSemidefinite2x2Bridge, Float64)] end function PolyJuMP.bridges( ::Type{<:MOI.AbstractVectorFunction}, ::Type{<:DiagonallyDominantConeTriangle}, ) - return [Bridges.Constraint.DiagonallyDominantBridge] + return [(Bridges.Constraint.DiagonallyDominantBridge, Float64)] end function PolyJuMP.bridges( ::Type{<:MOI.AbstractVectorFunction}, ::Type{<:ScaledDiagonallyDominantConeTriangle}, ) - return [Bridges.Constraint.ScaledDiagonallyDominantBridge] + return [(Bridges.Constraint.ScaledDiagonallyDominantBridge, Float64)] +end + +function _bridge_coefficient_type( + ::Type{SOSPolynomialSet{S,M,MV,C}}, +) where {S,M,MV,C} + return _complex(Float64, matrix_cone_type(C)) end function PolyJuMP.bridges( ::Type{<:MOI.AbstractVectorFunction}, - ::Type{<:SOSPolynomialSet{<:AbstractAlgebraicSet}}, + S::Type{<:SOSPolynomialSet{<:AbstractAlgebraicSet}}, ) - return [Bridges.Constraint.SOSPolynomialBridge] + return [( + Bridges.Constraint.SOSPolynomialBridge, + _bridge_coefficient_type(S), + )] end function PolyJuMP.bridges( ::Type{<:MOI.AbstractVectorFunction}, - ::Type{<:SOSPolynomialSet{<:BasicSemialgebraicSet}}, + S::Type{<:SOSPolynomialSet{<:BasicSemialgebraicSet}}, ) - return [Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge] + return [( + Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge, + _bridge_coefficient_type(S), + )] end # Syntax: `@constraint(model, a >= b, SOSCone())` @@ -339,10 +351,17 @@ function JuMP.build_constraint( return build_constraint(_error, f, extra; kws...) end +_promote_coef_type(::Type{V}, ::Type) where {V<:JuMP.AbstractVariableRef} = V +_promote_coef_type(::Type{F}, ::Type{T}) where {F,T} = promote_type(F, T) + function JuMP.build_constraint(_error::Function, p, cone::SOSLikeCone; kws...) - coefs = PolyJuMP.non_constant_coefficients(p) monos = MP.monomials(p) set = JuMP.moi_set(cone, monos; kws...) + _coefs = PolyJuMP.non_constant_coefficients(p) + # If a polynomial with real coefficients is used with the Hermitian SOS + # cone, we want to promote the coefficients to complex + T = _bridge_coefficient_type(typeof(set)) + coefs = convert(Vector{_promote_coef_type(eltype(_coefs), T)}, _coefs) shape = PolyJuMP.PolynomialShape(monos) return PolyJuMP.bridgeable( JuMP.VectorConstraint(coefs, set, shape), @@ -351,29 +370,6 @@ function JuMP.build_constraint(_error::Function, p, cone::SOSLikeCone; kws...) ) end -_non_constant(a::Vector{T}) where {T} = convert.(MOI.ScalarAffineFunction{T}, a) -_non_constant(a::Vector{<:JuMP.AbstractJuMPScalar}) = moi_function.(a) -_non_constant(a::Vector{<:MOI.AbstractFunction}) = a - -# Add constraint with `p` having coefficients being MOI functions. -# This is needed as a workaround as JuMP does not support complex numbers yet. -# We can remove it when https://github.com/jump-dev/JuMP.jl/pull/2391 is done -# We overload `JuMP.add_constraint` to avoid clash with the name. -function JuMP.add_constraint(model::MOI.ModelLike, p, cone::SOSLikeCone; kws...) - coefs = MOI.Utilities.vectorize(_non_constant(MP.coefficients(p))) - monos = MP.monomials(p) - set = JuMP.moi_set(cone, monos; kws...) - return MOI.add_constraint(model, coefs, set) -end -function JuMP.add_constraint(model::JuMP.Model, p, cone::SOSLikeCone; kws...) - ci = JuMP.add_constraint(JuMP.backend(model), p, cone; kws...) - return JuMP.ConstraintRef( - model, - ci, - PolyJuMP.PolynomialShape(MP.monomials(p)), - ) -end - struct ValueNotSupported <: Exception end function Base.showerror(io::IO, ::ValueNotSupported) return print( diff --git a/src/variable.jl b/src/variable.jl index af01309bc..44ec94087 100644 --- a/src/variable.jl +++ b/src/variable.jl @@ -1,13 +1,13 @@ export DSOSPoly, SDSOSPoly, SOSPoly function PolyJuMP.bridges(::Type{<:PositiveSemidefinite2x2ConeTriangle}) - return [Bridges.Variable.PositiveSemidefinite2x2Bridge] + return [(Bridges.Variable.PositiveSemidefinite2x2Bridge, Float64)] end function PolyJuMP.bridges(::Type{<:ScaledDiagonallyDominantConeTriangle}) - return [Bridges.Variable.ScaledDiagonallyDominantBridge] + return [(Bridges.Variable.ScaledDiagonallyDominantBridge, Float64)] end function PolyJuMP.bridges(::Type{<:CopositiveInnerCone}) - return [Bridges.Variable.CopositiveInnerBridge] + return [(Bridges.Variable.CopositiveInnerBridge, Float64)] end function JuMP.value(p::GramMatrix{<:JuMP.AbstractJuMPScalar})