Update to MultivariatePolynomials v0.5

Project.toml

@@ -19,15 +19,15 @@ SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 
 [compat]
-DynamicPolynomials = "0.4"
+DynamicPolynomials = "0.5"
 JuMP = "1"
 MathOptInterface = "1"
-MultivariateBases = "0.1"
-MultivariateMoments = "0.3"
-MultivariatePolynomials = "0.4"
+MultivariateBases = "0.2"
+MultivariateMoments = "0.4"
+MultivariatePolynomials = "0.5"
 MutableArithmetics = "1"
-PolyJuMP = "0.6"
+PolyJuMP = "0.7"
 Reexport = "1"
-SemialgebraicSets = "0.2"
-SumOfSquares = "0.6"
+SemialgebraicSets = "0.3"
+SumOfSquares = "0.7"
 julia = "1.6"

src/MBext/MBextra.jl

@@ -4,23 +4,23 @@ function MB.change_basis(p::MP.AbstractPolynomialLike, Basis::Type{<:MB.Abstract
 end
 
 function MB.change_basis(p::MP.AbstractPolynomialLike, basis::MB.AbstractPolynomialBasis)
-    coeffs = promote_type(Float64, coefficienttype(p))[]
+    coeffs = Vector{promote_type(Float64, MP.coefficient_type(p))}(undef, length(basis))
     rem = p
     mons = monomials(basis)
-    for i = 1:length(basis)
+    for i in reverse(eachindex(basis))
         c, rem = divrem(rem, mons[i])
         if iszero(c)
-            push!(coeffs, zero(eltype(coeffs)))
+            coeffs[i] = zero(eltype(coeffs))
         else
             @assert length(monomials(c)) == 1 && MP.isconstant(first(monomials(c)))
-            push!(coeffs, first(coefficients(c)))
+            coeffs[i] = first(coefficients(c))
         end
     end
     idx = findall(!iszero, coeffs)
     return coeffs[idx], mons[idx]
 end
 
-function MB.monomials(basis::MB.AbstractPolynomialBasis)
+function MP.monomials(basis::MB.AbstractPolynomialBasis)
     if basis isa MB.AbstractMonomialBasis
         return basis.monomials
     else

src/MomentOpt.jl

@@ -2,8 +2,7 @@ module MomentOpt
 
 using Reexport
 
-using MathOptInterface
-const MOI = MathOptInterface
+import MathOptInterface as MOI
 
 using MutableArithmetics
 const MA = MutableArithmetics

src/approximate.jl

@@ -32,6 +32,11 @@ function dual_variable(model::JuMP.Model, con::AbstractGMPConstraint, sense::MOI
     return variable
 end
 
+# FIXME sometimes we get `GenericAffExpr{BigFloat}` so we need this
+function _fix(p)
+    return convert(MP.similar_type(typeof(p), JuMP.AffExpr), p)
+end
+
 function approximate!(model::GMPModel, ::AbstractDualMode)
     # init dual variables
     meas_data = Dict(i => measure_data(model, i) for i in keys(gmp_variables(model)))
@@ -42,7 +47,7 @@ function approximate!(model::GMPModel, ::AbstractDualMode)
     end
 
     # init dual constraints
-    PT = polynomialtype(eltype(variables(gmp_variables(model)[1].v)), JuMP.AffExpr)
+    PT = polynomial_type(eltype(variables(gmp_variables(model)[1].v)), JuMP.AffExpr)
     dlhs = Dict{Int, PT}()
     for i in keys(gmp_variables(model))
         dlhs[i] = zero(PT)        
@@ -92,9 +97,14 @@ function approximate!(model::GMPModel, ::AbstractDualMode)
             p = MA.add!!(p, coefficients(obj_expr)[id])
         end
         scheme_parts[i] = approximation_scheme(model, v, meas_data[i])
-        just[i] = [dual_scheme_variable(approximation_model(model), sp)*sp.polynomial for sp in scheme_parts[i].schemeparts]
+        just[i] = [
+            _fix(dual_scheme_variable(approximation_model(model), sp)*sp.polynomial)
+            for sp in scheme_parts[i].schemeparts
+        ]
 
-        p = MA.sub!!(p, sum(just[i] )) 
+        for justi in just[i]
+            p = MA.sub!!(p, justi) 
+        end
         dcon[i] = @constraint approximation_model(model) scheme_parts[i].denominator*p == 0
     end
 
@@ -205,4 +215,3 @@ function approximate!(model::GMPModel, ::AbstractPrimalMode)
 
     return nothing
 end
-

src/approximationscheme.jl

@@ -135,7 +135,7 @@ function approximation_scheme(scheme::PutinarScheme, K::AbstractBasicSemialgebra
             end
         end
     end
-    return Scheme(schemeparts, one(polynomialtype(eltype(vars))), unique!(monos))
+    return Scheme(schemeparts, one(polynomial_type(eltype(vars))), unique!(monos))
 end
 
 
@@ -215,6 +215,5 @@ function approximation_scheme(scheme::SchmuedgenScheme, K::AbstractBasicSemialge
         push!(schemeparts, SchemePart(eq, mons, MOI.Zeros(length(mons))))
         append!(monos, monomials(mons))
     end
-    return Scheme(schemeparts, one(polynomialtype(eltype(vars))), unique!(monos))
+    return Scheme(schemeparts, one(polynomial_type(eltype(vars))), unique!(monos))
 end
-

src/gmppostproc.jl

@@ -38,7 +38,7 @@ export atomic
 Tries to exctract the atomic support of a measure. 
 """
 function atomic(vref::GMPVariableRef; tol = 1e-3)
-	optmeas = extractatoms(moment_matrix(vref), tol)
+	optmeas = atomic_measure(moment_matrix(vref), tol)
 	if typeof(optmeas) == Nothing
         return nothing
 	else
@@ -55,14 +55,14 @@ function min_val(x::Pair{<:Vector{<:MP.AbstractVariable}, <:Vector{<:Number}},
                  poly::AbstractPolynomialLike)
     p = subs(poly, x)
     t = first(variables(p))
-    set = algebraicset([differentiate(p, t)])
+    set = algebraic_set([differentiate(p, t)])
     Y = [sol[1] for sol in set]
     return Y[argmin([p(t => y) for y in Y])]
 end
 
 function christoffel(vref::GMPVariableRef; regpar = 1e-8)
     M = moment_matrix(vref)
-    eva, eve = eigen(MM.getmat(M))
+    eva, eve = eigen(MM.value_matrix(M))
     return sum( (dot(eve[:, i] / √(eva[i] + regpar), M.basis.monomials))^2 for i = 1:length(eva))
 end
 

src/momexpr.jl

@@ -92,7 +92,7 @@ function Base.sum(mev::Vector{MomentExpr{T}}) where T
     if T <: Number
         coefs = T[]
     else
-        coefs = polynomialtype(T)[]
+        coefs = polynomial_type(T)[]
     end
     vars = GMPVariableRef[]
     for me in mev

test/affexpr.jl

@@ -67,8 +67,8 @@
             @test MO.degree(m3) == 0
             @test MO.degree(m4) == 1
 
-            @test sprint(show, m5) == "⟨-x + 0.5, μ⟩ + ⟨0.5, ν⟩"
-            @test sprint(show, m4 - m5) == "⟨1.5x + 0.5, μ⟩ + ⟨0.5x - 0.5, ν⟩"
+            @test sprint(show, m5) == "⟨0.5 - x, μ⟩ + ⟨0.5, ν⟩"
+            @test sprint(show, m4 - m5) == "⟨0.5 + 1.5x, μ⟩ + ⟨-0.5 + 0.5x, ν⟩"
             @test sprint(show, MomentExpr(Int[], MO.GMPVariableRef[])) == "⟨0, 0⟩"
 
             @test eltype([Mom(1, μ + ν), μ]) == AbstractJuMPScalar

test/approximate.jl

@@ -58,8 +58,8 @@ end
         optimize!(m)
 
         @test value(mu) isa MultivariateMoments.Measure
-        @test primal_justification(mu) isa Array{Array{Float64,N} where N, 1}
-        @test dual_justification(mu) isa Polynomial{true, Float64}
+        @test primal_justification(mu) isa Vector{Array{Float64,N} where N}
+        @test dual_justification(mu) isa Polynomial
         @test all(isapprox.(residual(con), 0.0; atol = 1e-3))
         @test dual(con) isa Vector{Float64}
 
@@ -69,9 +69,9 @@ end
 
         @test value(mu) isa MultivariateMoments.Measure
         @test primal_justification(mu) isa Array{Array{Float64,2}, 1}
-        @test dual_justification(mu) isa Polynomial{true, Float64}
+        @test dual_justification(mu) isa Polynomial
         @test all(isapprox.(residual(con), 0.0; atol = 1e-3))
-        @test dual(con) isa Polynomial{true, Float64}
+        @test dual(con) isa Polynomial
 
         @test isapprox(objective_value(m), integrate(1, mu); atol = 1e-3)
         @test isapprox(dual_objective_value(m), integrate(1, mu); atol = 1e-3)
@@ -90,7 +90,7 @@ end
 
         @test value(mu) isa MultivariateMoments.Measure
         @test primal_justification(mu) isa Array{Array{Float64,N} where N, 1}
-        @test dual_justification(mu) isa Polynomial{true, Float64}
+        @test dual_justification(mu) isa Polynomial
         @test all(isapprox.(residual(con), 0.0; atol = 1e-3))
         @test dual(con) isa Vector{Float64}
 
@@ -101,9 +101,9 @@ end
 
         @test value(mu) isa MultivariateMoments.Measure
         @test primal_justification(mu) isa Array{Array{Float64,2}, 1}
-        @test dual_justification(mu) isa Polynomial{true, Float64}
+        @test dual_justification(mu) isa Polynomial
         @test all(isapprox.(residual(con), 0.0; atol = 1e-3))
-        @test dual(con) isa Polynomial{true, Float64}
+        @test dual(con) isa Polynomial
 
     end
 
@@ -204,9 +204,9 @@ end
     f = x^2*y^2 + x^4 - y^2 + 1
     gmp = GMPModel()
     mu = @variable gmp μ Meas([x, y]; support = @set(x^2-x^4>=0 && 1-x^2>=0&& x-y^2==0))
-    MO.approximation_info(gmp).approx_vrefs[1] = MO.VarApprox(MultivariateMoments.Measure([0.04301405019233425, 3.725430564600456e-17, 0.0944511405218325, -3.667541691651414e-17, 0.20739824193451498, 0.09445115553035512, -6.278937572688163e-17, 0.207398112278084, 2.2189511821478325e-17, 0.20739824605238844, 6.290839177153608e-18, 0.4554099572402294, 0.45540988993201603, 4.3223432627485784e-17, 1.0000002747597936], Monomial{true}[x^4, x^3*y, x^2*y^2, x*y^3, y^4, x^3, x^2*y, x*y^2, y^3, x^2, x*y, y^2, x, y, 1]), nothing, nothing)
+    MO.approximation_info(gmp).approx_vrefs[1] = MO.VarApprox(MultivariateMoments.Measure([0.04301405019233425, 3.725430564600456e-17, 0.0944511405218325, -3.667541691651414e-17, 0.20739824193451498, 0.09445115553035512, -6.278937572688163e-17, 0.207398112278084, 2.2189511821478325e-17, 0.20739824605238844, 6.290839177153608e-18, 0.4554099572402294, 0.45540988993201603, 4.3223432627485784e-17, 1.0000002747597936], [x^4, x^3*y, x^2*y^2, x*y^3, y^4, x^3, x^2*y, x*y^2, y^3, x^2, x*y, y^2, x, y, 1]), nothing, nothing)
     @test length(atomic(mu)) == 2
-    @test integrate(f, mu) == 0.682055508233731
+    @test integrate(f, mu) ≈ 0.682055508233731
 
     g =  graph(mu, [x]; regpar = 1e-6)
     @test g isa Function

test/defaultmeasures.jl

@@ -1,6 +1,6 @@
 function test_moments(μ::AnalyticMeasure, maxdegree, moments)
     mons = monomials(maxdegree_basis(μ, maxdegree))
-    for (m, mm) in zip(mons, moments)
+    for (m, mm) in zip(mons, reverse(moments))
         @test isapprox(integrate(m, μ), mm; atol = 1e-6)
     end
 end
@@ -46,7 +46,7 @@ end
         mons = monomials([x, y], 0:2);
         z = [1/100*(sum(XX[i]^e[1]*(2*XX[i]-1)^e[2] for i = 1:100)) for e in exponents.(mons)];
         μ = moment_measure(MultivariateMoments.Measure(z, mons))
-        test_moments(μ, 2, z)
+        test_moments(μ, 2, reverse(z))
     end
 
     @testset "(Partial) Integration" begin
@@ -70,4 +70,3 @@ end
         end 
     end
 end
-

test/objects.jl

@@ -5,7 +5,7 @@
 
         mone = MO._mono_one([x, y])
         @test isone(mone)
-        @test mone isa Monomial{true}
+        @test mone isa Monomial
         @test variables(mone) == [x, y]
 
     end
@@ -31,8 +31,8 @@
         @test integrate(1, λ) == 1
         μ = Meas([x, y]; support = S2, basis = ChebyshevBasis)
         @test_throws MethodError integrate(p2, μ)
-        @test monomials(maxdegree_basis(μ, 2)) == [2.0*x^2 - 1.0, x*y, 2.0*y^2 - 1.0, x, y, 1.0]
-        @test monomials(MO.covering_basis(μ, 2)) == Polynomial{true, typeof(1.0)}[1.0]
+        @test monomials(maxdegree_basis(μ, 2)) == reverse([2.0*x^2 - 1.0, x*y, 2.0*y^2 - 1.0, x, y, 1.0])
+        @test monomials(MO.covering_basis(μ, 2)) == polynomial_type(x, Float64)[1.0]
     end
 
     @testset "DefaultMeasures" begin

test/variables.jl

@@ -4,7 +4,7 @@
         @polyvar x y
         mu = Meas([x, y])
         v = GMPVariable(mu)
-        @test MO.object_type(v) == VariableMeasure{PolyVar{true}, FullSpace, MonomialBasis, PutinarScheme}
+        @test MO.object_type(v) == VariableMeasure{typeof(x), FullSpace, MonomialBasis, PutinarScheme}
 
         _error = _ -> ""
         info = VariableInfo(true, 1, true, 1, true, 1, true, 1, false, false)