Merge 87f41a9 into 830487b

src/PolyJuMP.jl

@@ -13,6 +13,7 @@ include("basis.jl")
 using MathOptInterface
 const MOI = MathOptInterface
 
+include("attributes.jl")
 include("zero_polynomial.jl")
 
 # Bridges

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

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})
 

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

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

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},

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)]