Merge 44744ca into 9be00a2

REQUIRE

@@ -1,4 +1,5 @@
 julia 0.6
 MultivariatePolynomials 0.1.0
+MultivariateMoments 0.0.1
 SemialgebraicSets 0.0.1
 JuMP 0.17.1

src/PolyJuMP.jl

@@ -3,95 +3,65 @@ __precompile__()
 module PolyJuMP
 
 using MultivariatePolynomials
+using MultivariateMoments
 using SemialgebraicSets
 using JuMP
 export getslack, setpolymodule!
 
 # Polynomial Constraint
 
-export ZeroPoly, NonNegPoly
+export ZeroPoly, NonNegPoly, NonNegPolyMatrix
+abstract type PolynomialSet end
+struct ZeroPoly <: PolynomialSet end
+struct NonNegPoly <: PolynomialSet end
+struct NonNegPolyMatrix <: PolynomialSet end
 
-struct ZeroPoly end
-struct NonNegPoly end
-
-mutable struct PolyConstraint <: JuMP.AbstractConstraint
-    p # typically either be a polynomial or a Matrix of polynomials
-    set
-    kwargs::Vector{Any}
-    polymodule::Nullable{Module}
-    domain::AbstractSemialgebraicSet
-    delegate::Nullable
-    function PolyConstraint(p, s)
-        new(p, s, Any[], nothing, FullSpace(), nothing)
-    end
-end
-function setpolymodule!(c::PolyConstraint, pm::Module)
-    c.polymodule = pm
+struct PolyConstraint{PT, ST<:PolynomialSet} <: JuMP.AbstractConstraint
+    p::PT # typically either be a polynomial or a Matrix of polynomials
+    set::ST
 end
-getpolymodule(c::PolyConstraint) = get(c.polymodule)
-
 const PolyConstraintRef = ConstraintRef{Model, PolyConstraint}
 
-function JuMP.addconstraint(m::Model, c::PolyConstraint; domain::AbstractSemialgebraicSet=FullSpace(), kwargs...)
-    setpolymodule!(c, getpolymodule(m))
-    c.domain = domain
-    c.kwargs = kwargs
-    polyconstr = getpolyconstr(m)
-    push!(polyconstr, c)
+# Responsible for getting slack and dual values
+abstract type ConstraintDelegate end
+
+function JuMP.addconstraint(m::Model, pc::PolyConstraint; domain::AbstractSemialgebraicSet=FullSpace(), kwargs...)
+    delegates = getdelegates(m)
+    delegate = addpolyconstraint!(m, pc.p, pc.set, domain; kwargs...)
+    push!(delegates, delegate)
     m.internalModelLoaded = false
-    PolyConstraintRef(m, length(polyconstr))
+    PolyConstraintRef(m, length(delegates))
 end
 
-function getdelegate(c::PolyConstraintRef, s::Symbol)
-    delegate = getpolyconstr(c.m)[c.idx].delegate
-    if isnull(delegate)
-        error("$(string(s)) value not defined for constraint with index $c. Check that the model was properly solved.")
-    end
-    get(delegate)
-end
-
-function getslack(c::PolyConstraintRef)
-    getslack(getdelegate(c, :Slack))
-end
-function JuMP.getdual(c::PolyConstraintRef)
-    getdual(getdelegate(c, :Dual))
-end
+getdelegate(c::PolyConstraintRef) = getdelegates(c.m)[c.idx]
+getslack(c::PolyConstraintRef) = getslack(getdelegate(c))
+JuMP.getdual(c::PolyConstraintRef) = getdual(getdelegate(c))
 
 # PolyJuMP Data
-
-type PolyData
-    polyconstr::Vector{PolyConstraint}
-    polymodule::Nullable{Module}
-    function PolyData()
-        new(PolyConstraint[], nothing)
+type Data
+    # Delegates for polynomial constraints created
+    delegates::Vector{ConstraintDelegate}
+    # Default set for Poly{true}
+    nonnegpolyvardefault::Nullable
+    # Default set for NonNegPoly
+    nonnegpolydefault::Nullable
+    # Default set for NonNegPolyMatrix
+    nonnegpolymatrixdefault::Nullable
+    function Data()
+        new(ConstraintDelegate[], nothing, nothing, nothing)
     end
 end
 
 function getpolydata(m::JuMP.Model)
     if !haskey(m.ext, :Poly)
-        m.solvehook = solvehook
-        m.ext[:Poly] = PolyData()
+        m.ext[:Poly] = Data()
     end
     m.ext[:Poly]
 end
-
-function getpolyconstr(m::JuMP.Model)
-    getpolydata(m).polyconstr
-end
-
-setpolymodule!(m::JuMP.Model, pm::Module) = setpolymodule!(getpolydata(m), pm)
-setpolymodule!(data::PolyData, pm::Module) = data.polymodule = pm
-
-getpolymodule(m::JuMP.Model) = getpolymodule(getpolydata(m))
-function getpolymodule(data::PolyData)
-    if isnull(data.polymodule)
-        return PolyJuMP
-    end
-    get(data.polymodule)
-end
+getdelegates(m::JuMP.Model) = getpolydata(m).delegates
 
 include("macros.jl")
-include("solve.jl")
+include("default.jl")
 include("default_methods.jl")
 
 end # module

src/default.jl

@@ -0,0 +1,26 @@
+setdefault!(m::JuMP.Model, S::Type, T::Type) = setdefault!(getpolydata(m), S, T)
+setdefault!(data::Data, ::Type{Poly{true}}, S::Type) = data.nonnegpolyvardefault = S
+setdefault!(data::Data, ::Type{NonNegPoly}, S::Type) = data.nonnegpolydefault = S
+setdefault!(data::Data, ::Type{NonNegPolyMatrix}, S::Type) = data.nonnegpolymatrixdefault = S
+
+getdefault(m::JuMP.Model, s) = getdefault(getpolydata(m), s)
+
+# includes Poly{false} but also custom types defined by other modules
+getdefault(data::Data, p) = p
+getdefault(data::Data, p::Poly{true, MS}) where MS = getdefault(data, Poly{true}){MS}(p.x)
+getdefault(data::Data, ::Type{Poly{true}}) = get_default(data.nonnegpolyvardefault)
+
+# includes Poly{false} but also custom types defined by other modules
+getdefault(data::Data, s::Union{NonNegPoly, NonNegPolyMatrix}) = getdefault(data, typeof(s))()
+getdefault(data::Data, ::Type{NonNegPoly}) = get_default(data.nonnegpolydefault)
+getdefault(data::Data, ::Type{NonNegPolyMatrix}) = get_default(data.nonnegpolymatrixdefault)
+
+setpolymodule!(m::JuMP.Model, pm::Module) = setpolymodule!(getpolydata(m), pm)
+setpolymodule!(data::Data, pm::Module) = pm.setdefaults!(data)
+
+function get_default(d::Nullable)
+    if isnull(d)
+        error("PolyJuMP is just a JuMP extension for modelling Polynomial Optimization: it does not implement any reformulation. To use automatic sums of squares (SOS) reformulations, install the SumOfSquares Julia package and try \`using SumOfSquares\` and \`setpolymodule!(SumOfSquares)\` or use \`SOSModel\` instead of \`Model\`.")
+    end
+    get(d)
+end

src/default_methods.jl

@@ -1,8 +1,5 @@
-polytype(m::JuMP.Model, p::Poly) = getpolymodule(m).polytype(m, p, p.x)
-polytype(m::JuMP.Model, p, X::AbstractVector) = getpolymodule(m).polytype(m, p, monovec(X))
-
 # Free polynomial
-
+JuMP.variabletype(m::JuMP.Model, p::Poly{false}) = polytype(m, p, p.x)
 polytype(m::JuMP.Model, ::Poly{false}, x::AbstractVector{MT}) where MT<:AbstractMonomial = MultivariatePolynomials.polynomialtype(MT, JuMP.Variable)
 
 # x should be sorted and without duplicates
@@ -14,27 +11,41 @@ function createpoly(m::JuMP.Model, p::Poly{false, :Gram}, category::Symbol)
     _createpoly(m, p, monomials(sum(p.x)^2), category)
 end
 
-function createpoly(m::JuMP.Model, p::Union{Poly{false, :Default,Array{T}}, Poly{false, :Classic,Array{T}}}, category::Symbol) where T <: MultivariatePolynomials.AbstractMonomial
-    _createpoly(m, p, monovec(p.x), category)
+function createpoly(m::JuMP.Model, p::Union{Poly{false, :Default}, Poly{false, :Classic}}, category::Symbol)
+    _createpoly(m, p, p.x, category)
 end
 
-function createpoly(m::JuMP.Model, p::Union{PolyJuMP.Poly{false, :Default}, PolyJuMP.Poly{false, :Classic}}, category::Symbol)
-    MultivariatePolynomials.polynomial([JuMP.Variable(m, -Inf, Inf, category)*p.x[i] for i in 1:length(p.x)])
+# NonNegPoly and NonNegPolyMatrix
+addpolyconstraint!(m::JuMP.Model, p, s::Union{NonNegPoly, NonNegPolyMatrix}, domain; kwargs...) = addpolyconstraint!(m, p, getdefault(m, s), domain; kwargs...)
+
+# ZeroPoly
+struct ZeroConstraint{MT <: AbstractMonomial, MVT <: AbstractVector{MT}, JC <: JuMP.AbstractConstraint} <: ConstraintDelegate
+    zero_constraints::Vector{JuMP.ConstraintRef{JuMP.Model, JC}} # These are typically affine or quadratic equality constraints
+    x::MVT
+end
+function ZeroConstraint(zero_constraints::Vector{JuMP.ConstraintRef{JuMP.Model, JC}}, x::MVT) where {MT <: AbstractMonomial, MVT <: AbstractVector{MT}, JC <: JuMP.AbstractConstraint}
+    ZeroConstraint{MT, MVT, JC}(zero_constraints, x)
 end
 
+JuMP.getdual(c::ZeroConstraint) = MultivariateMoments.measure(getdual.(c.zero_constraints), c.x)
+
 function addpolyconstraint!(m::JuMP.Model, p, s::ZeroPoly, domain::FullSpace)
     constraints = JuMP.constructconstraint!.(coefficients(p), :(==))
-    JuMP.addVectorizedConstraint(m, constraints)
+    zero_constraints = JuMP.addVectorizedConstraint(m, constraints)
+    ZeroConstraint(zero_constraints, monomials(p))
 end
 
 function addpolyconstraint!(m::JuMP.Model, p, s::ZeroPoly, domain::AbstractAlgebraicSet)
     addpolyconstraint!(m, rem(p, ideal(domain)), s, FullSpace())
 end
 
-function addpolyconstraint!(m::JuMP.Model, p, s::ZeroPoly, domain::BasicSemialgebraicSet)
-    addpolyconstraint!(m,  p, NonNegPoly(), domain)
-    addpolyconstraint!(m, -p, NonNegPoly(), domain)
-    nothing
+struct ZeroConstraintWithDomain{DT<:ConstraintDelegate} <: ConstraintDelegate
+    lower::DT
+    upper::DT
 end
 
-addpolyconstraint!(m::JuMP.Model, args...) = error("PolyJuMP is just a JuMP extension for modelling Polynomial Optimization: it does not implement any reformulation. To use automatic sums of squares (SOS) reformulations, install the SumOfSquares Julia package and try \`using SumOfSquares\` and \`setpolymodule!(SumOfSquares)\` or use \`SOSModel\` instead of \`Model\`.")
+function addpolyconstraint!(m::JuMP.Model, p, s::ZeroPoly, domain::BasicSemialgebraicSet)
+    lower = addpolyconstraint!(m,  p, NonNegPoly(), domain)
+    upper = addpolyconstraint!(m, -p, NonNegPoly(), domain)
+    ZeroConstraintWithDomain(lower, upper)
+end

src/macros.jl

@@ -1,8 +1,4 @@
-using JuMP
-import JuMP: validmodel, addtoexpr_reorder
-using Base.Meta
-
-export Poly, @set
+export Poly
 
 function JuMP.getvalue(t::AbstractTerm{<:JuMP.AbstractJuMPScalar})
     getvalue(coefficient(t)) * monomial(t)
@@ -14,20 +10,20 @@ end
 abstract type AbstractPoly end
 
 # x is a vector of monomials to be used to construct a polynomial variable
-# if MT is Gram, x represents the monomials of the form x^T Q x
-# if MT is Classic, it represents the monomials of the form a^T x
-# if MT is Default, it depends on whether the polynomials is constructed as nonnegative or not:
+# if MS is Gram, x represents the monomials of the form x^T Q x
+# if MS is Classic, it represents the monomials of the form a^T x
+# if MS is Default, it depends on whether the polynomials is constructed as nonnegative or not:
 # For a nonnegative polynomial, it corresponds to Gram, otherwise it corresponds to Classic.
-struct Poly{P, MT, MV} <: AbstractPoly
+struct Poly{P, MS, MT<:MultivariatePolynomials.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPoly
     x::MV
 end
-Poly{P, MT}(x::MV) where {P, MT, MV} = Poly{P, MT, MV}(x)
+Poly{P, MS}(x::AbstractVector{MT}) where {P, MS, MT<:MultivariatePolynomials.AbstractMonomial} = Poly{P, MS, MT, typeof(x)}(x)
+Poly{P, MS}(x) where {P, MS} = Poly{P, MS}(monovec(x))
 Poly{P}(x) where P = Poly{P, :Default}(x)
 Poly(x) = Poly{false}(x)
 
-function JuMP.variabletype(m::Model, p::AbstractPoly)
-    getpolymodule(m).polytype(m, p)
-end
+JuMP.variabletype(m::JuMP.Model, p::Poly{true}) = JuMP.variabletype(m, getdefault(m, p))
+
 function cvarchecks(_error::Function, lowerbound::Number, upperbound::Number, start::Number; extra_kwargs...)
     for (kwarg, _) in extra_kwargs
         _error("Unrecognized keyword argument $kwarg")
@@ -58,17 +54,17 @@ end
 function JuMP.constructvariable!(m::Model, p::AbstractPoly, _error::Function, lowerbound::Number, upperbound::Number, category::Symbol, basename::AbstractString, start::Number; extra_kwargs...)
     cvarchecks(_error, lowerbound, upperbound, start; extra_kwargs...)
     _warnbounds(_error, p, lowerbound, upperbound)
-    getpolymodule(m).createpoly(m, p, category == :Default ? :Cont : category)
+    createpoly(m, getdefault(m, p), category == :Default ? :Cont : category)
 end
 
 function JuMP.constructconstraint!(p::AbstractPolynomialLike, sense::Symbol)
-    PolyConstraint(sense == :(<=) ? -p : p, sense == :(==) ? ZeroPoly() : NonNegPoly())
+    JuMP.constructconstraint!(sense == :(<=) ? -p : p, sense == :(==) ? ZeroPoly() : NonNegPoly())
 end
 
-function JuMP.constructconstraint!{PolyT<:AbstractPolynomialLike}(p::Union{PolyT, AbstractMatrix{PolyT}}, s)
+function JuMP.constructconstraint!(p::Union{AbstractPolynomialLike, AbstractMatrix{<:AbstractPolynomialLike}}, s)
     PolyConstraint(p, s)
 end
 # there is already a method for AbstractMatrix in PSDCone in JuMP so we need a more specific here to avoid ambiguity
-function JuMP.constructconstraint!{PolyT<:AbstractPolynomialLike}(p::AbstractMatrix{PolyT}, s::PSDCone)
-    PolyConstraint(p, s)
+function JuMP.constructconstraint!(p::AbstractMatrix{<:AbstractPolynomialLike}, s::PSDCone)
+    PolyConstraint(p, NonNegPolyMatrix())
 end

test/constraint.jl

@@ -24,9 +24,18 @@ _iszero(m, p::AbstractArray) = all(q -> _iszero(m, q), p)
     #@test macroexpand(:(@constraint(m, p >= 0, domain = (@set x >= -1 && x <= 1, domain = y >= -1 && y <= 1)))).head == :error
     @test macroexpand(:(@constraint(m, p + 0, domain = (@set x >= -1 && x <= 1)))).head == :error
 
+    function testcon(m, cref, set::ZeroPoly, p, ineqs, eqs, kwargs=[])
+        @test isa(cref, ConstraintRef{Model, PolyJuMP.PolyConstraint})
+        c = PolyJuMP.getdelegate(cref)
+        if isempty(ineqs)
+            @test c isa PolyJuMP.ZeroConstraint
+        else
+            @test c isa PolyJuMP.ZeroConstraintWithDomain
+        end
+    end
     function testcon(m, cref, set, p, ineqs, eqs, kwargs=[])
         @test isa(cref, ConstraintRef{Model, PolyJuMP.PolyConstraint})
-        c = PolyJuMP.getpolyconstr(m)[cref.idx]
+        c = PolyJuMP.getdelegate(cref)
         @test c.set == set
         @test c.kwargs == kwargs
         # == between JuMP affine expression is not accurate, e.g. β + α != α + β
@@ -49,15 +58,13 @@ _iszero(m, p::AbstractArray) = all(q -> _iszero(m, q), p)
 
     f(x, y) = @set x + y == 2
     dom = @set x^2 + y^2 == 1 && x^3 + x*y^2 + y >= 1
-    testcon(m, @constraint(m, p >= q + 1, domain = @set y >= 1 && dom), NonNegPoly(), p - q - 1, [y-1, x^3 + x*y^2 + y - 1], [x^2 + y^2 - 1])
-    testcon(m, @constraint(m, p <= q), NonNegPoly(), q - p, [], [])
-    testcon(m, @constraint(m, q - p in NonNegPoly()), NonNegPoly(), q - p, [], [])
-    testcon(m, @constraint(m, p + q >= 0, domain = @set x == y^3), NonNegPoly(), p + q, [], [x - y^3])
+    testcon(m, @constraint(m, p >= q + 1, domain = @set y >= 1 && dom), TestPolyModule.TestNonNegConstraint(), p - q - 1, [y-1, x^3 + x*y^2 + y - 1], [x^2 + y^2 - 1])
+    testcon(m, @constraint(m, p <= q), TestPolyModule.TestNonNegConstraint(), q - p, [], [])
+    testcon(m, @constraint(m, q - p in NonNegPoly()), TestPolyModule.TestNonNegConstraint(), q - p, [], [])
+    testcon(m, @constraint(m, p + q >= 0, domain = @set x == y^3), TestPolyModule.TestNonNegConstraint(), p + q, [], [x - y^3])
     testcon(m, @constraint(m, p == q, domain = @set x == 1 && f(x, y)), ZeroPoly(), p - q, [], [x - 1, x + y - 2])
+    testcon(m, @constraint(m, p == q, domain = dom), ZeroPoly(), p - q, [x^3 + x*y^2 + y - 1], [x^2 + y^2 - 1])
     testcon(m, @constraint(m, p - q in ZeroPoly(), domain = @set x == 1 && f(x, y)), ZeroPoly(), p - q, [], [x - 1, x + y - 2])
-    testcon(m, @SDconstraint(m, [p q; q 0] ⪰ [0 0; 0 p]), PSDCone(), [p q; q -p], [], [])
-    testcon(m, @constraint(m, p <= q, maxdegree=1), NonNegPoly(), q - p, [], [], [(:maxdegree, 1)])
-
-    cref = @constraint(m, p >= 0)
-    @test_throws ErrorException PolyJuMP.getdelegate(cref, :Test)
+    testcon(m, @SDconstraint(m, [p q; q 0] ⪰ [0 0; 0 p]), TestPolyModule.TestNonNegMatrixConstraint(), [p q; q -p], [], [])
+    testcon(m, @constraint(m, p <= q, maxdegree=1), TestPolyModule.TestNonNegConstraint(), q - p, [], [], [(:maxdegree, 1)])
 end

test/polymodule.jl

@@ -1,8 +1,16 @@
 @testset "PolyModule" begin
     m = Model()
     # Triggers the creation of polydata
-    @test isnull(PolyJuMP.getpolydata(m).polymodule)
-    @test PolyJuMP.getpolymodule(m) == PolyJuMP
+    @test isnull(PolyJuMP.getpolydata(m).nonnegpolyvardefault)
+    @test_throws ErrorException PolyJuMP.getdefault(m, Poly{true})
+    @test isnull(PolyJuMP.getpolydata(m).nonnegpolydefault)
+    @test_throws ErrorException PolyJuMP.getdefault(m, NonNegPoly)
+    @test isnull(PolyJuMP.getpolydata(m).nonnegpolymatrixdefault)
+    @test_throws ErrorException PolyJuMP.getdefault(m, NonNegPolyMatrix)
     setpolymodule!(m, TestPolyModule)
-    @test PolyJuMP.getpolymodule(m) == TestPolyModule
+    @test PolyJuMP.getdefault(m, Poly{true}) == TestPolyModule.TestPoly
+    @test PolyJuMP.getdefault(m, NonNegPoly) == TestPolyModule.TestNonNegConstraint
+    @test PolyJuMP.getdefault(m, NonNegPolyMatrix) == TestPolyModule.TestNonNegMatrixConstraint
+    PolyJuMP.setdefault!(m, NonNegPolyMatrix, TestPolyModule.TestNonNegConstraint)
+    @test PolyJuMP.getdefault(m, NonNegPolyMatrix) == TestPolyModule.TestNonNegConstraint
 end