Merge pull request #40 from JuliaOpt/bl/basis

Update for polynomial basis

Author: blegat

src/SumOfSquares.jl

@@ -21,7 +21,6 @@ include("variable.jl")
 include("constraint.jl")
 
 function setdefaults!(data::PolyJuMP.Data)
-    PolyJuMP.setdefault!(data, PolyJuMP.Poly{true}, SOSPoly)
     PolyJuMP.setdefault!(data, PolyJuMP.NonNegPoly, SOSCone)
     PolyJuMP.setdefault!(data, PolyJuMP.NonNegPolyMatrix, SOSMatrixCone)
 end

src/constraint.jl

@@ -35,14 +35,14 @@ certificate_monomials(c::SOSConstraint) = c.slack.x
 PolyJuMP.getslack(c::SOSConstraint) = getvalue(c.slack)
 JuMP.getdual(c::SOSConstraint) = getdual(c.zero_constraint)
 
-function PolyJuMP.addpolyconstraint!(m::JuMP.Model, P::Matrix{PT}, ::SOSMatrixCone, domain::AbstractBasicSemialgebraicSet) where PT <: APL
+function PolyJuMP.addpolyconstraint!(m::JuMP.Model, P::Matrix{PT}, ::SOSMatrixCone, domain::AbstractBasicSemialgebraicSet, basis) where PT <: APL
     n = Base.LinAlg.checksquare(P)
     if !issymmetric(P)
         throw(ArgumentError("The polynomial matrix constrained to be SOS must be symmetric"))
     end
     y = [similarvariable(PT, gensym()) for i in 1:n]
     p = dot(y, P * y)
-    PolyJuMP.addpolyconstraint!(m, p, SOSCone(), domain)
+    PolyJuMP.addpolyconstraint!(m, p, SOSCone(), domain, basis)
 end
 
 function _createslack(m, x, set::SOSLikeCones)
@@ -57,16 +57,16 @@ function _createslack(m, x, set::CoSOSLikeCones)
     _matplus(_createslack(m, x, _nococone(set)), _matposynomial(m, x))
 end
 
-function PolyJuMP.addpolyconstraint!(m::JuMP.Model, p, set::SOSSubCones, domain::AbstractAlgebraicSet)
+function PolyJuMP.addpolyconstraint!(m::JuMP.Model, p, set::SOSSubCones, domain::AbstractAlgebraicSet, basis)
     r = rem(p, ideal(domain))
     X = getmonomialsforcertificate(monomials(r))
     slack = _createslack(m, X, set)
     q = r - slack
-    zero_constraint = PolyJuMP.addpolyconstraint!(m, q, ZeroPoly(), domain)
+    zero_constraint = PolyJuMP.addpolyconstraint!(m, q, ZeroPoly(), domain, basis)
     SOSConstraint(slack, zero_constraint)
 end
 
-function PolyJuMP.addpolyconstraint!(m::JuMP.Model, p, set::SOSSubCones, domain::BasicSemialgebraicSet;
+function PolyJuMP.addpolyconstraint!(m::JuMP.Model, p, set::SOSSubCones, domain::BasicSemialgebraicSet, basis;
                             mindegree=MultivariatePolynomials.mindegree(p),
                             maxdegree=MultivariatePolynomials.maxdegree(p))
     for q in domain.p
@@ -85,5 +85,5 @@ function PolyJuMP.addpolyconstraint!(m::JuMP.Model, p, set::SOSSubCones, domain:
         s2 = createpoly(m, _varconetype(set)(monomials(variables(p), mindegree_s:maxdegree_s)), :Cont)
         p -= s2 * q
     end
-    PolyJuMP.addpolyconstraint!(m, p, set, domain.V)
+    PolyJuMP.addpolyconstraint!(m, p, set, domain.V, basis)
 end

src/variable.jl

@@ -6,24 +6,22 @@ end
 
 for poly in (:DSOSPoly, :SDSOSPoly, :SOSPoly)
     @eval begin
-        struct $poly{MS, MT<:MultivariatePolynomials.AbstractMonomial, MV<:AbstractVector{MT}} <: PolyJuMP.AbstractPoly
-            x::MV
+        struct $poly{PB<:PolyJuMP.AbstractPolynomialBasis} <: PolyJuMP.AbstractPoly
+            polynomial_basis::PB
         end
-        $poly{MS}(x::AbstractVector{MT}) where {MS, MT} = $poly{MS, MT, typeof(x)}(x)
-        $poly{MS}(x) where MS = $poly{MS}(monovec(x))
-        $poly(x) = $poly{:Default}(x)
+        $poly(x::AbstractVector{<:MultivariatePolynomials.AbstractPolynomialLike}) = $poly(MonomialBasis(x))
     end
 end
 
-const PosPoly{MT, MV} = Union{DSOSPoly{MT, MV}, SDSOSPoly{MT, MV}, SOSPoly{MT, MV}}
+const PosPoly{PB} = Union{DSOSPoly{PB}, SDSOSPoly{PB}, SOSPoly{PB}}
 
-JuMP.variabletype(m::JuMP.Model, p::PosPoly) = PolyJuMP.polytype(m, p, p.x)
-PolyJuMP.polytype(m::JuMP.Model, ::PosPoly, x::MVT) where {MT<:AbstractMonomial, MVT<:AbstractVector{MT}} = MatPolynomial{JuMP.Variable, MT, MVT}
+JuMP.variabletype(m::JuMP.Model, p::PosPoly) = PolyJuMP.polytype(m, p, p.polynomial_basis)
+PolyJuMP.polytype(m::JuMP.Model, ::PosPoly, basis::PolyJuMP.MonomialBasis{MT, MV}) where {MT<:AbstractMonomial, MV<:AbstractVector{MT}} = MatPolynomial{JuMP.Variable, MT, MV}
 
-function _constraintmatpoly!(m, p, ::Union{SOSPoly, Poly{true}})
+function _constraintmatpoly!(m, p::MatPolynomial, ::SOSPoly)
     push!(m.varCones, (:SDP, [p.Q[i, j].col for i in 1:size(p.Q, 1) for j in i:size(p.Q, 2)]))
 end
-function _constraintmatpoly!(m, p, ::DSOSPoly)
+function _constraintmatpoly!(m, p::MatPolynomial, ::DSOSPoly)
     n = length(p.x)
     Q = Matrix{JuMP.Variable}(n, n)
     for i in 1:n
@@ -58,21 +56,15 @@ function _matpolynomial(m, x::AbstractVector{<:AbstractMonomial}, category::Symb
         MatPolynomial{JuMP.Variable}((i, j) -> Variable(m, lb, Inf, category), x)
     end
 end
-function _createpoly(m::JuMP.Model, set::PosPoly, x::AbstractVector{<:AbstractMonomial}, category::Symbol)
-    p = _matpolynomial(m, x, category)
-    if length(x) > 1
+function _createpoly(m::JuMP.Model, set::PosPoly, basis::PolyJuMP.MonomialBasis, category::Symbol)
+    p = _matpolynomial(m, basis.monomials, category)
+    if length(basis.monomials) > 1
         _constraintmatpoly!(m, p, set)
     end
     p
 end
-function PolyJuMP.createpoly(m::JuMP.Model, p::Union{PosPoly{:Default}, PosPoly{:Gram}}, category::Symbol)
-    _createpoly(m, p, monovec(p.x), category)
-end
-function PolyJuMP.createpoly(m::JuMP.Model, pp::PosPoly{:Classic}, category::Symbol)
-    p = _createpoly(m, pp, getmonomialsforcertificate(pp.x), category)
-    # The coefficients of a monomial not in Z do not all have to be zero, only their sum
-    addpolyconstraint!(m, removemonomials(Polynomial(p), p.x), ZeroPoly(), FullSpace())
-    p
+function PolyJuMP.createpoly(m::JuMP.Model, p::PosPoly, category::Symbol)
+    _createpoly(m, p, p.polynomial_basis, category)
 end
 
 # Defer other methods to defaults in PolyJuMP

test/constraint.jl

@@ -2,5 +2,5 @@
     @polyvar x
     m = SOSModel()
     # Throws an Argument Error because it should be symmetric to be an SOS matrix
-    @test_throws ArgumentError PolyJuMP.addpolyconstraint!(m, [1 x; -x 0], SOSMatrixCone(), BasicSemialgebraicSet{Int, polynomialtype(x, Int)}())
+    @test_throws ArgumentError PolyJuMP.addpolyconstraint!(m, [1 x; -x 0], SOSMatrixCone(), BasicSemialgebraicSet{Int, polynomialtype(x, Int)}(), PolyJuMP.MonomialBasis)
 end

test/sosdemo5.jl

@@ -33,7 +33,7 @@
             # -- Q(x)'s -- : sums of squares
             # Monomial vector: [x1; ... x8]
             Q = Vector{MatPolynomial{JuMP.Variable, monomialtype(x[1]), monovectype(x[1])}}(4)
-            @variable m Q[1:4] >= 0 Poly{true}(Z)
+            @variable m Q[1:4] SOSPoly(Z)
 
             # -- r's -- : constant sum of squares
             Z = monomials(x, 0)

test/sosdemo6.jl

@@ -16,7 +16,7 @@
     m = SOSModel(solver = solver)
 
     Z = monomials(x, 0:1)
-    @variable m p1 >= 0 SOSPoly(Z)
+    @variable m p1 SOSPoly(Z)
 
     Z = monomials(x, 0:2)
     @variable m p[1:5] Poly(Z)

test/switchedsystems.jl

@@ -18,8 +18,7 @@
     expected_ub = [√2, 1]
     function testlyap(d, γ, expected_status)
         m = SOSModel(solver = solver)
-        X = monomials(x, d)
-        @variable m p Poly{false, :Gram}(X)
+        @variable m p Poly(monomials(x, 2d))
         # p strictly positive
         q = sum(x.^(2*d))
         @constraint m p >= q

test/variable.jl

@@ -2,8 +2,8 @@
     @polyvar x
     X = emptymonovec(typeof(x^2))
     m = SOSModel()
-    @inferred PolyJuMP.createpoly(m, SOSPoly{:Gram}(X), :Cont)
-    @test PolyJuMP.createpoly(m, SOSPoly{:Gram}(X), :Cont) == 0
+    @inferred PolyJuMP.createpoly(m, SOSPoly(X), :Cont)
+    @test PolyJuMP.createpoly(m, SOSPoly(X), :Cont) == 0
 end
 @testset "MatPolynomial getvalue" begin
     m = JuMP.Model()