From a7db59dc91110e26d6cbcf4189824dc7f249470e Mon Sep 17 00:00:00 2001 From: tweisser Date: Thu, 4 Jun 2020 13:20:52 -0600 Subject: [PATCH 1/3] working with tw/MB --- src/atomic.jl | 4 ++-- src/expectation.jl | 10 +++++----- src/extract.jl | 4 ++-- src/measure.jl | 27 ++++++++++++++++++++------- src/moment.jl | 29 +++++++++++++++++++++-------- src/moment_matrix.jl | 2 +- test/expectation.jl | 11 +++++++++++ test/measure.jl | 21 +++++++++++++++++++-- test/runtests.jl | 2 +- 9 files changed, 82 insertions(+), 28 deletions(-) diff --git a/src/atomic.jl b/src/atomic.jl index da51502..3a36f9f 100644 --- a/src/atomic.jl +++ b/src/atomic.jl @@ -45,9 +45,9 @@ end measure(η::AtomicMeasure, x) = Measure(η, x) function Measure(η::AtomicMeasure{T}, x::AbstractVector{TT}) where {T, TT} - Measure{T, monomialtype(TT), monovectype(x)}(η, x) + Measure{T, MB.MonomialBasis}(η, x) end -function Measure{T, MT, MVT}(η::AtomicMeasure{T}, x) where {T, MT, MVT} +function Measure{T, BT}(η::AtomicMeasure{T}, x) where {T, BT} X = monovec(x) sum(atom.weight * dirac(X, η.variables => atom.center) for atom in η.atoms) end diff --git a/src/expectation.jl b/src/expectation.jl index 6f10de8..82aca40 100644 --- a/src/expectation.jl +++ b/src/expectation.jl @@ -1,16 +1,16 @@ const APL = AbstractPolynomialLike -function _expectation(μ::Measure, p::APL, f) +function _expectation(μ::Measure{T, BT}, p::APL, f) where {T, BT} i = 1 s = 0 - for t in terms(p) - while i <= length(μ.x) && monomial(t) != μ.x[i] + for (c, m) in zip(coefficients(p, BT), MB.basis_covering_monomials(BT, monomials(p))) + while i <= length(μ.x) && m != μ.x[i] i += 1 end if i > length(μ.x) - error("The polynomial $p has a nonzero term $t with monomial $(t.x) for which the expectation is not known in $μ") + error("The polynomial $p has a nonzero term $(c*m) with basis function $(m) for which the expectation is not known in $μ") end - s += f(μ.a[i], coefficient(t)) + s += f(μ.a[i], c) i += 1 end s diff --git a/src/extract.jl b/src/extract.jl index aed46e8..d6fdc3c 100644 --- a/src/extract.jl +++ b/src/extract.jl @@ -60,7 +60,7 @@ function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol) r = length(pivots) U = U[keep, :] end - monos = basis.monomials + monos = basis.elements β = monovec(monos[m + 1 .- pivots]) # monovec makes sure it stays sorted, TypedPolynomials wouldn't let it sorted function equation(i) if iszero(r) # sum throws ArgumentError: reducing over an empty collection is not allowed, if r is zero @@ -181,7 +181,7 @@ function extractatoms(ν::MomentMatrix{T}, ranktol, args...) where T vars = variables(μ) A = Matrix{T}(undef, length(μ.x), r) for i in 1:r - A[:, i] = dirac(μ.x, vars => centers[i]).a + A[:, i] = dirac(μ.x.elements, vars => centers[i]).a end weights = A \ μ.a isf = isfinite.(weights) diff --git a/src/measure.jl b/src/measure.jl index 807dbf7..91ad712 100644 --- a/src/measure.jl +++ b/src/measure.jl @@ -1,17 +1,22 @@ export measure, dirac -export variables, monomials, moments +export variables, base_functions, moments # If a monomial is not in x, it does not mean that the moment is zero, it means that it is unknown/undefined -struct Measure{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractMeasure{T} +struct Measure{T, B<:MB.AbstractPolynomialBasis} <: AbstractMeasure{T} a::Vector{T} - x::MVT - - function Measure{T, MT, MVT}(a::Vector{T}, x::MVT) where {T, MT, MVT} + x::B + function Measure(a::Vector{T}, x::B) where {T, B<:MB.AbstractPolynomialBasis} @assert length(a) == length(x) - new(a, x) + new{T, B}(a, x) end end -Measure(a::AbstractVector{T}, x::AbstractVector{TT}) where {T, TT <: AbstractTermLike} = Measure{T, monomialtype(TT), monovectype(x)}(monovec(a, x)...) + +basis_functions_type(::Measure{T, BT}) where {T,BT} = BT + +function Measure(a::AbstractVector{T}, x::AbstractVector{TT}) where {T, TT <: AbstractTermLike} + b, y = monovec(a, x) + return Measure(b, MB.MonomialBasis(y)) +end """ measure(a, X::AbstractVector{<:AbstractMonomial}) @@ -34,6 +39,14 @@ Returns an iterator over the monomials of `μ` sorted in the decreasing order. """ MP.monomials(μ::Measure) = μ.x +""" + base_functions(μ::AbstractMeasureLike) + +Returns an iterator over the base_functions of `μ` sorted in the decreasing order. +""" +base_functions(μ::Measure) = μ.x + + """ moments(μ::AbstractMeasureLike) diff --git a/src/moment.jl b/src/moment.jl index fda571e..99a01c2 100644 --- a/src/moment.jl +++ b/src/moment.jl @@ -1,16 +1,16 @@ -export moment, moment_value, monomial +export moment, moment_value, monomial, base_function -struct Moment{T, MT <: AbstractMonomial} <: AbstractMoment{T} +struct Moment{T, PT <: AbstractPolynomialLike} <: AbstractMoment{T} α::T - x::MT + x::PT end """ - moment(α, m::AbstractMonomial) + moment(α, m::AbstractPolynomialLike) -Creates the moment of the monomial `m` of value `α`. +Creates the moment of the base function `m` of value `α`. """ -moment(α, m::AbstractMonomial) = Moment(α, m) +moment(α, m::AbstractPolynomialLike) = Moment(α, m) """ moment_value(m::AbstractMomentLike) @@ -34,8 +34,21 @@ Calling `monomial(moment(3.1, x*y^2))` should return `x*y^2`. """ MP.monomial(m::Moment) = m.x -for f in [:variables, :nvariables, :exponents, :degree, :powers] +""" + base_function(m::AbstractMomentLike) + +Returns the base_function of the moment `m`. + +## Examples + +Calling `base_function(moment(3.1, x*y^2))` should return `x*y^2`. +""" +base_function(m::Moment) = m.x + +for f in [:variables, :nvariables, :exponents, :powers] @eval begin - MP.$f(m::AbstractMomentLike) = $f(monomial(m)) + MP.$f(m::AbstractMomentLike) = $f(base_function(m)) end end + +MP.degree(m::AbstractMomentLike) = MP.maxdegree(base_function(m)) diff --git a/src/moment_matrix.jl b/src/moment_matrix.jl index ce68b7c..7aa7a1d 100644 --- a/src/moment_matrix.jl +++ b/src/moment_matrix.jl @@ -69,7 +69,7 @@ end function measure(ν::MomentMatrix{T, <:MB.MonomialBasis}) where T n = length(ν.basis) - monos = ν.basis.monomials + monos = ν.basis.elements measure(ν.Q.Q, [monos[i] * monos[j] for i in 1:n for j in 1:i]) end diff --git a/test/expectation.jl b/test/expectation.jl index 94be370..f5a5d36 100644 --- a/test/expectation.jl +++ b/test/expectation.jl @@ -7,4 +7,15 @@ @test_throws ErrorException dot(x[1] * x[2] * x[3], m) @test dot(0.5 * x[1] * x[2]^2, m) == 2.0 @test dot(m, x[1] * x[3]) == 3 + + @testset "Expectation - MB" begin + Mod.@polyvar x[1:2] + p = x[2] - 2x[1]*x[2]^2 + 3x[2]*x[1] - 5x[1]^3 + m = measure([1, 0, 2, 3, 0, 4, 0, 2, 5, 0], maxdegree_basis(ChebyshevBasis, x, 3)) + @test MultivariateMoments.expectation(m, p) == MultivariateMoments.expectation(p, m) == 4.25 + @test_throws ErrorException dot((x[1] * x[2])^2, m) + @test dot(0.5 * x[1] * x[2]^2, m) == 1.0 + @test dot(m, x[1] * x[2]) == 4 + end + end diff --git a/test/measure.jl b/test/measure.jl index b7d63a4..e6929ac 100644 --- a/test/measure.jl +++ b/test/measure.jl @@ -3,8 +3,8 @@ @test_throws ArgumentError measure([1, 2], [x, x*y, y]) @test_throws ArgumentError measure([1, 2, 3, 4], [x, x*y, y]) μ = measure([1, 0, 2, 3], [x^2*y^2, y*x^2, x*y*x^2, x*y^2]) - @test monomials(μ) == [x^3*y, x^2*y^2, x^2*y, x*y^2] - @test monomial.(moments(μ)) == [x^3*y, x^2*y^2, x^2*y, x*y^2] + @test base_functions(μ) == MonomialBasis(monovec([x^3*y, x^2*y^2, x^2*y, x*y^2])) + @test base_function.(moments(μ)) == [x^3*y, x^2*y^2, x^2*y, x*y^2] @test MultivariateMoments.moment_value.(moments(μ)) == [2, 1, 0, 3] @test all(nvariables.(moments(μ)) .== 2) @test degree.(moments(μ)) == [4, 4, 3, 3] @@ -13,4 +13,21 @@ @test (2 * μ).a == [4, 2, 0, 6] @test (μ * 3).a == [6, 3, 0, 9] #@test_throws ArgumentError measure([1], [x]) + measure([1], [y]) + + @testset "Measure MB" begin + Mod.@polyvar x y + @test_throws AssertionError measure([1, 0, 2, 3, 0], maxdegree_basis(ChebyshevBasis, [x, y], 2)) + @test_throws AssertionError measure([1, 0, 2, 3, 0, 4], maxdegree_basis(ChebyshevBasis, [x, y], 3)) + μ = measure([1, 0, 2, 3, 0, 4], maxdegree_basis(ChebyshevBasis, [x, y], 2)) + @test base_functions(μ) == maxdegree_basis(ChebyshevBasis, [x, y], 2) + @test base_function.(moments(μ)) == [2.0x^2 - 1.0, x*y, 2.0y^2 - 1.0, x, y, 1.0] + @test MultivariateMoments.moment_value.(moments(μ)) == [1, 0, 2, 3, 0, 4] + @test all(nvariables.(moments(μ)) .== 2) + @test degree.(moments(μ)) == [2, 2, 2, 1, 1, 0] + @test μ.a == [1, 0, 2, 3, 0, 4] + @test (-μ).a == -[1, 0, 2, 3, 0, 4] + @test (2 * μ).a == [1, 0, 2, 3, 0, 4].*2 + @test (μ * 3).a == [1, 0, 2, 3, 0, 4].*3 + + end end diff --git a/test/runtests.jl b/test/runtests.jl index 4403f74..0fc222b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,5 +1,5 @@ using MultivariatePolynomials -import MultivariateBases +using MultivariateBases const MB = MultivariateBases using MultivariateMoments using SemialgebraicSets From 30ded84f36d94baa99ec227f622b2c3a8d2f9a04 Mon Sep 17 00:00:00 2001 From: tweisser Date: Thu, 11 Jun 2020 21:27:39 -0600 Subject: [PATCH 2/3] incorporate suggested changes --- src/atomic.jl | 2 +- src/expectation.jl | 6 ++++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/atomic.jl b/src/atomic.jl index 3a36f9f..6bc0374 100644 --- a/src/atomic.jl +++ b/src/atomic.jl @@ -45,7 +45,7 @@ end measure(η::AtomicMeasure, x) = Measure(η, x) function Measure(η::AtomicMeasure{T}, x::AbstractVector{TT}) where {T, TT} - Measure{T, MB.MonomialBasis}(η, x) + Measure{T, MB.MonomialBasis{monomialtype(TT), monovectype(x)}}(η, x) end function Measure{T, BT}(η::AtomicMeasure{T}, x) where {T, BT} X = monovec(x) diff --git a/src/expectation.jl b/src/expectation.jl index 82aca40..1551800 100644 --- a/src/expectation.jl +++ b/src/expectation.jl @@ -3,12 +3,14 @@ const APL = AbstractPolynomialLike function _expectation(μ::Measure{T, BT}, p::APL, f) where {T, BT} i = 1 s = 0 - for (c, m) in zip(coefficients(p, BT), MB.basis_covering_monomials(BT, monomials(p))) + basis = MB.basis_covering_monomials(BT, monomials(p)) + for (c, m) in zip(coefficients(p, basis), basis) while i <= length(μ.x) && m != μ.x[i] i += 1 end if i > length(μ.x) - error("The polynomial $p has a nonzero term $(c*m) with basis function $(m) for which the expectation is not known in $μ") + error("The base function $m has a non-zero multiplier $c in $p, + but its expectation is not known in $μ") end s += f(μ.a[i], c) i += 1 From 5df0b32562a244f36d1b928be66946b8f9336f87 Mon Sep 17 00:00:00 2001 From: tweisser Date: Thu, 16 Jul 2020 08:23:25 -0600 Subject: [PATCH 3/3] change field names of Measure --- src/expectation.jl | 6 +++--- src/extract.jl | 6 +++--- src/measure.jl | 36 ++++++++++++++++++------------------ test/measure.jl | 16 ++++++++-------- 4 files changed, 32 insertions(+), 32 deletions(-) diff --git a/src/expectation.jl b/src/expectation.jl index 1551800..a77dd6b 100644 --- a/src/expectation.jl +++ b/src/expectation.jl @@ -5,14 +5,14 @@ function _expectation(μ::Measure{T, BT}, p::APL, f) where {T, BT} s = 0 basis = MB.basis_covering_monomials(BT, monomials(p)) for (c, m) in zip(coefficients(p, basis), basis) - while i <= length(μ.x) && m != μ.x[i] + while i <= length(μ.basis) && m != μ.basis[i] i += 1 end - if i > length(μ.x) + if i > length(μ.basis) error("The base function $m has a non-zero multiplier $c in $p, but its expectation is not known in $μ") end - s += f(μ.a[i], c) + s += f(μ.values[i], c) i += 1 end s diff --git a/src/extract.jl b/src/extract.jl index d6fdc3c..27b3f04 100644 --- a/src/extract.jl +++ b/src/extract.jl @@ -179,11 +179,11 @@ function extractatoms(ν::MomentMatrix{T}, ranktol, args...) where T # Determine weights μ = measure(ν) vars = variables(μ) - A = Matrix{T}(undef, length(μ.x), r) + A = Matrix{T}(undef, length(μ.basis), r) for i in 1:r - A[:, i] = dirac(μ.x.elements, vars => centers[i]).a + A[:, i] = dirac(μ.basis.elements, vars => centers[i]).values end - weights = A \ μ.a + weights = A \ μ.values isf = isfinite.(weights) weights = weights[isf] centers = centers[isf] diff --git a/src/measure.jl b/src/measure.jl index 91ad712..5649a0b 100644 --- a/src/measure.jl +++ b/src/measure.jl @@ -1,20 +1,20 @@ export measure, dirac export variables, base_functions, moments -# If a monomial is not in x, it does not mean that the moment is zero, it means that it is unknown/undefined +# If a monomial is not in basis, it does not mean that the moment is zero, it means that it is unknown/undefined struct Measure{T, B<:MB.AbstractPolynomialBasis} <: AbstractMeasure{T} - a::Vector{T} - x::B - function Measure(a::Vector{T}, x::B) where {T, B<:MB.AbstractPolynomialBasis} - @assert length(a) == length(x) - new{T, B}(a, x) + values::Vector{T} + basis::B + function Measure(values::Vector{T}, basis::B) where {T, B<:MB.AbstractPolynomialBasis} + @assert length(values) == length(basis) + new{T, B}(values, basis) end end -basis_functions_type(::Measure{T, BT}) where {T,BT} = BT +basis_functions_type(::Measure{T, BT}) where {T, BT} = BT -function Measure(a::AbstractVector{T}, x::AbstractVector{TT}) where {T, TT <: AbstractTermLike} - b, y = monovec(a, x) +function Measure(values::AbstractVector{T}, basis::AbstractVector{TT}) where {T, TT <: AbstractTermLike} + b, y = monovec(values, basis) return Measure(b, MB.MonomialBasis(y)) end @@ -30,21 +30,21 @@ measure(a, X) = Measure(a, X) Returns the variables of `μ` in decreasing order. Just like in MultivariatePolynomials, it could contain variables of zero degree in every monomial. """ -MP.variables(μ::Measure) = variables(μ.x) +MP.variables(μ::Measure) = variables(μ.basis) """ monomials(μ::AbstractMeasureLike) Returns an iterator over the monomials of `μ` sorted in the decreasing order. """ -MP.monomials(μ::Measure) = μ.x +MP.monomials(μ::Measure) = μ.basis """ base_functions(μ::AbstractMeasureLike) Returns an iterator over the base_functions of `μ` sorted in the decreasing order. """ -base_functions(μ::Measure) = μ.x +base_functions(μ::Measure) = μ.basis """ @@ -52,14 +52,14 @@ base_functions(μ::Measure) = μ.x Returns an iterator over the moments of `μ` sorted in decreasing order of monomial. """ -moments(μ::Measure) = map((α, x) -> moment(α, x), μ.a, μ.x) +moments(μ::Measure) = map((α, x) -> moment(α, x), μ.values, μ.basis) -Base.:(*)(α, μ::Measure) = measure(α * μ.a, μ.x) -Base.:(*)(μ::Measure, α) = measure(μ.a * α, μ.x) -Base.:(-)(μ::Measure) = measure(-μ.a, μ.x) +Base.:(*)(α, μ::Measure) = measure(α * μ.values, μ.basis) +Base.:(*)(μ::Measure, α) = measure(μ.values * α, μ.basis) +Base.:(-)(μ::Measure) = measure(-μ.values, μ.basis) function Base.:(+)(μ::Measure, ν::Measure) - @assert μ.x == ν.x - measure(μ.a + ν.a, μ.x) + @assert μ.basis == ν.basis + measure(μ.values + ν.values, μ.basis) end """ diff --git a/test/measure.jl b/test/measure.jl index e6929ac..80088de 100644 --- a/test/measure.jl +++ b/test/measure.jl @@ -8,10 +8,10 @@ @test MultivariateMoments.moment_value.(moments(μ)) == [2, 1, 0, 3] @test all(nvariables.(moments(μ)) .== 2) @test degree.(moments(μ)) == [4, 4, 3, 3] - @test μ.a == [2, 1, 0, 3] - @test (-μ).a == [-2, -1, 0, -3] - @test (2 * μ).a == [4, 2, 0, 6] - @test (μ * 3).a == [6, 3, 0, 9] + @test μ.values == [2, 1, 0, 3] + @test (-μ).values == [-2, -1, 0, -3] + @test (2 * μ).values == [4, 2, 0, 6] + @test (μ * 3).values == [6, 3, 0, 9] #@test_throws ArgumentError measure([1], [x]) + measure([1], [y]) @testset "Measure MB" begin @@ -24,10 +24,10 @@ @test MultivariateMoments.moment_value.(moments(μ)) == [1, 0, 2, 3, 0, 4] @test all(nvariables.(moments(μ)) .== 2) @test degree.(moments(μ)) == [2, 2, 2, 1, 1, 0] - @test μ.a == [1, 0, 2, 3, 0, 4] - @test (-μ).a == -[1, 0, 2, 3, 0, 4] - @test (2 * μ).a == [1, 0, 2, 3, 0, 4].*2 - @test (μ * 3).a == [1, 0, 2, 3, 0, 4].*3 + @test μ.values == [1, 0, 2, 3, 0, 4] + @test (-μ).values == -[1, 0, 2, 3, 0, 4] + @test (2 * μ).values == [1, 0, 2, 3, 0, 4].*2 + @test (μ * 3).values == [1, 0, 2, 3, 0, 4].*3 end end