diff --git a/Project.toml b/Project.toml index 8bcb12b..88dd471 100644 --- a/Project.toml +++ b/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" diff --git a/src/MBext/MBextra.jl b/src/MBext/MBextra.jl index 67a1516..e7d69a7 100644 --- a/src/MBext/MBextra.jl +++ b/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, coefficienttype(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 diff --git a/src/MomentOpt.jl b/src/MomentOpt.jl index 7546558..eb5430f 100644 --- a/src/MomentOpt.jl +++ b/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 diff --git a/src/approximate.jl b/src/approximate.jl index 5dc63d1..2fc58aa 100644 --- a/src/approximate.jl +++ b/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))) @@ -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 - diff --git a/src/approximationscheme.jl b/src/approximationscheme.jl index cc32bac..384cd35 100644 --- a/src/approximationscheme.jl +++ b/src/approximationscheme.jl @@ -217,4 +217,3 @@ function approximation_scheme(scheme::SchmuedgenScheme, K::AbstractBasicSemialge end return Scheme(schemeparts, one(polynomialtype(eltype(vars))), unique!(monos)) end - diff --git a/test/affexpr.jl b/test/affexpr.jl index d2c19f5..5542eba 100644 --- a/test/affexpr.jl +++ b/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 diff --git a/test/approximate.jl b/test/approximate.jl index 220bbfc..60edd3b 100644 --- a/test/approximate.jl +++ b/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 diff --git a/test/defaultmeasures.jl b/test/defaultmeasures.jl index 4a77bb7..5533f17 100644 --- a/test/defaultmeasures.jl +++ b/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 - diff --git a/test/objects.jl b/test/objects.jl index 539f874..5d0a654 100644 --- a/test/objects.jl +++ b/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 diff --git a/test/variables.jl b/test/variables.jl index 8c694d5..a099a95 100644 --- a/test/variables.jl +++ b/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)