working with tw/MB

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{monomialtype(TT), monovectype(x)}}(η, 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

src/expectation.jl

@@ -1,16 +1,18 @@
 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]
+    basis = MB.basis_covering_monomials(BT, monomials(p))
+    for (c, m) in zip(coefficients(p, basis), basis)
+        while i <= length(μ.basis) && m != μ.basis[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 $μ")
+        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], coefficient(t))
+        s += f(μ.values[i], c)
         i += 1
     end
     s

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
@@ -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, 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]

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}
-    a::Vector{T}
-    x::MVT
-
-    function Measure{T, MT, MVT}(a::Vector{T}, x::MVT) where {T, MT, MVT}
-        @assert length(a) == length(x)
-        new(a, x)
+# 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}
+    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
-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(values::AbstractVector{T}, basis::AbstractVector{TT}) where {T, TT <: AbstractTermLike}
+    b, y = monovec(values, basis)
+    return Measure(b, MB.MonomialBasis(y))
+end
 
 """
     measure(a, X::AbstractVector{<:AbstractMonomial})
@@ -25,28 +30,36 @@ 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) = μ.basis
+
 
 """
     moments(μ::AbstractMeasureLike)
 
 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
 
 """

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))

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
 

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

test/measure.jl

@@ -3,14 +3,31 @@
     @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]
-    @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
+        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 μ.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

test/runtests.jl

@@ -1,5 +1,5 @@
 using MultivariatePolynomials
-import MultivariateBases
+using  MultivariateBases
 const MB = MultivariateBases
 using MultivariateMoments
 using SemialgebraicSets