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}(η, 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,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

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)

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)
 

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

test/runtests.jl

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