diff --git a/src/PolyJuMP.jl b/src/PolyJuMP.jl index 019e920..f4e9177 100644 --- a/src/PolyJuMP.jl +++ b/src/PolyJuMP.jl @@ -13,6 +13,7 @@ include("basis.jl") using MathOptInterface const MOI = MathOptInterface +include("attributes.jl") include("zero_polynomial.jl") # Bridges diff --git a/src/attributes.jl b/src/attributes.jl new file mode 100644 index 0000000..704d421 --- /dev/null +++ b/src/attributes.jl @@ -0,0 +1,38 @@ +""" + ConstraintDual(N) + ConstraintDual() + +A constraint attribute for the vector of moments corresponding to the +constraint that a polynomial is zero in the full space in result `N`. If the +polynomial is constrained to be zero in an algebraic set, it is the moments +for the constraint once it is rewritten into an constraint on the full space. +If `N` is omitted, it is 1 by default. + +## Examples + +Consider the following program: +```julia +@variable(model, α) +@variable(model, β ≤ 1) + +using DynamicPolynomials +@polyvar x y +cref = @constraint(model, α * x - β * y == 0, domain = @set x == y) +``` +The constraint is equivalent to +```julia +@constraint(model, (α - β) * y == 0) +``` +for which the dual is the 1-element vector with the moment of `y` of value `-1`. +This is the result of `moments(cref)`. However, the dual of `cref` obtained by +`dual(cref)` is the 2-elements vector with both the moments of `x` and `y` +of value `-1`. +""" +struct MomentsAttribute <: MOI.AbstractConstraintAttribute + N::Int +end +MomentsAttribute() = MomentsAttribute(1) + +function MOI.is_set_by_optimize(::MomentsAttribute) + return true +end diff --git a/src/constraint.jl b/src/constraint.jl index 3431871..0818138 100644 --- a/src/constraint.jl +++ b/src/constraint.jl @@ -1,3 +1,5 @@ +export moments + abstract type PolynomialSet end function JuMP.in_set_string(print_mode, set::PolynomialSet) @@ -45,6 +47,15 @@ function JuMP.moi_set(::ZeroPoly, monos::AbstractVector{<:AbstractMonomial}; return ZeroPolynomialSet(domain, basis, monos) end +""" + moments(cref::JuMP.ConstraintRef) + +Return the [`MomentsAttribute`](@ref) of `cref`. +""" +function MultivariateMoments.moments(cref::JuMP.ConstraintRef) + return MOI.get(cref.model, MomentsAttribute(), cref) +end + """ bridges(F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet}) diff --git a/src/zero_polynomial_bridge.jl b/src/zero_polynomial_bridge.jl index 2dfd4dc..d86633f 100644 --- a/src/zero_polynomial_bridge.jl +++ b/src/zero_polynomial_bridge.jl @@ -39,13 +39,18 @@ function MOI.get(b::ZeroPolynomialBridge{T, F}, end # Indices -function MOI.delete(model::MOI.ModelLike, c::ZeroPolynomialBridge) - MOI.delete(model, c.zero_constraint) +function MOI.delete(model::MOI.ModelLike, bridge::ZeroPolynomialBridge) + MOI.delete(model, bridge.zero_constraint) end # Attributes, Bridge acting as a constraint function MOI.get(model::MOI.ModelLike, attr::Union{MOI.ConstraintPrimal, MOI.ConstraintDual}, - c::ZeroPolynomialBridge) - return MOI.get(model, attr, c.zero_constraint) + bridge::ZeroPolynomialBridge) + return MOI.get(model, attr, bridge.zero_constraint) +end +function MOI.get(model::MOI.ModelLike, attr::MomentsAttribute, + bridge::ZeroPolynomialBridge) + values = MOI.get(model, MOI.ConstraintDual(attr.N), bridge) + return measure(values, bridge.monomials) end diff --git a/src/zero_polynomial_in_algebraic_set_bridge.jl b/src/zero_polynomial_in_algebraic_set_bridge.jl index a1f031e..121a9ff 100644 --- a/src/zero_polynomial_in_algebraic_set_bridge.jl +++ b/src/zero_polynomial_in_algebraic_set_bridge.jl @@ -57,6 +57,8 @@ function MOI.delete(model::MOI.ModelLike, c::ZeroPolynomialInAlgebraicSetBridge) MOI.delete(model, c.zero_constraint) end +# Attributes, Bridge acting as a constraint + # TODO ConstraintPrimal # Let A be the linear map corresponding to A(p) = rem(p, ideal(set.domain)) @@ -73,3 +75,7 @@ function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual, I = ideal(bridge.domain) return [dot(rem(mono, I), μ) for mono in bridge.monomials] end +function MOI.get(model::MOI.ModelLike, attr::MomentsAttribute, + bridge::ZeroPolynomialInAlgebraicSetBridge) + return MOI.get(model, attr, bridge.zero_constraint) +end diff --git a/test/Tests/zero_polynomial.jl b/test/Tests/zero_polynomial.jl index c9d4b62..24f4c3a 100644 --- a/test/Tests/zero_polynomial.jl +++ b/test/Tests/zero_polynomial.jl @@ -32,11 +32,13 @@ function zero_polynomial_test(optimizer::MOI.AbstractOptimizer, @test dual_status(model) == MOI.FEASIBLE_POINT @test dual(UpperBoundRef(β)) ≈ -1.0 atol=atol rtol=rtol - μ = dual(cref) - @test μ isa AbstractMeasure{Float64} - @test length(moments(μ)) == 1 - @test moment_value(moments(μ)[1]) ≈ -1.0 atol=atol rtol=rtol - @test monomial(moments(μ)[1]) == x*y + + for μ in [dual(cref), moments(cref)] + @test μ isa AbstractMeasure{Float64} + @test length(moments(μ)) == 1 + @test moment_value(moments(μ)[1]) ≈ -1.0 atol=atol rtol=rtol + @test monomial(moments(μ)[1]) == x*y + end F = MOI.VectorAffineFunction{Float64} S = PolyJuMP.ZeroPolynomialSet{FullSpace,MonomialBasis,Monomial{true}, diff --git a/test/Tests/zero_polynomial_in_algebraic_set.jl b/test/Tests/zero_polynomial_in_algebraic_set.jl index 86152c0..5f2c2b3 100644 --- a/test/Tests/zero_polynomial_in_algebraic_set.jl +++ b/test/Tests/zero_polynomial_in_algebraic_set.jl @@ -31,6 +31,7 @@ function zero_polynomial_in_algebraic_set_test(optimizer, @test dual_status(model) == MOI.FEASIBLE_POINT @test dual(UpperBoundRef(β)) ≈ -1.0 atol=atol rtol=rtol + μ = dual(cref) @test μ isa AbstractMeasure{Float64} @test length(moments(μ)) == 2 @@ -39,8 +40,14 @@ function zero_polynomial_in_algebraic_set_test(optimizer, @test moment_value(moments(μ)[2]) ≈ -1.0 atol=atol rtol=rtol @test monomial(moments(μ)[2]) == y + μ = moments(cref) + @test μ isa AbstractMeasure{Float64} + @test length(moments(μ)) == 1 + @test moment_value(moments(μ)[1]) ≈ -1.0 atol=atol rtol=rtol + @test monomial(moments(μ)[1]) == y + F = MOI.VectorAffineFunction{Float64} - S = PolyJuMP.ZeroPolynomialSet{typeof(@set x == 1),MonomialBasis, + S = PolyJuMP.ZeroPolynomialSet{typeof(@set x == 1), MonomialBasis, monomialtype(x), monovectype(x)} @test MOI.get(model, MOI.ListOfConstraints()) == [ (MOI.SingleVariable, MOI.LessThan{Float64}), (F, S)]