Fixes

.github/workflows/ci.yml

@@ -42,6 +42,14 @@ jobs:
             ${{ runner.os }}-test-${{ env.cache-name }}-
             ${{ runner.os }}-test-
             ${{ runner.os }}-
+      - name: dev
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="StarAlgebras", rev="mk/non_monomial_basis"),
+              PackageSpec(name="MultivariateBases", rev="master"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:

.github/workflows/documentation.yml

@@ -15,7 +15,15 @@ jobs:
           # Build documentation on latest Julia release
           version: '1'
       - name: Install dependencies
-        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+        shell: julia --project=docs/ {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="StarAlgebras", rev="mk/non_monomial_basis"),
+              PackageSpec(name="MultivariateBases", rev="master"),
+              PackageSpec(path=pwd()),
+          ])
+          Pkg.instantiate()
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token

docs/Project.toml

@@ -1,6 +1,9 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
+MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
+StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 
 [compat]
 Documenter = "1"

docs/src/moments.md

@@ -12,24 +12,22 @@ An `AbstractMomentLike` is a type that can act like an `AbstractMoment` (it is s
 that is, it implements the following two functions
 ```@docs
 moment_value
-MultivariatePolynomials.monomial(::MultivariateMoments.Moment)
 ```
 
 ## Measure
 
-Given a monomials and a values for the moments, a "measure" can be created using the `measure` function
+Given a monomials and a values for the moments, a "moment vector" can be created using the `moment_vector` function
 ```@docs
-measure
+moment_vector
 ```
-The `measure` function returns an `AbstractMeasure` which is a subtype of `AbstractMeasureLike`.
+The `moment_vector` function returns an `AbstractMeasure` which is a subtype of `AbstractMeasureLike`.
 Note that it does not actually compute the probability density function of a measure having these moments, it simply stores a vector of moments belonging to a hypothetical measure.
 However, it acts like a measure when taking its scalar product with a polynomial.
 
 An `AbstractMeasureLike` is a type that can act like an `AbstractMeasure`,
 that is, it implements the following two functions
 ```@docs
 MultivariatePolynomials.variables(::MultivariateMoments.MomentVector)
-MultivariatePolynomials.monomials(::MultivariateMoments.MomentVector)
 MultivariatePolynomials.maxdegree(::MultivariateMoments.MomentVector)
 MultivariatePolynomials.mindegree(::MultivariateMoments.MomentVector)
 MultivariatePolynomials.extdegree(::MultivariateMoments.MomentVector)

src/MultivariateMoments.jl

@@ -4,6 +4,7 @@ using LinearAlgebra
 import SparseArrays
 
 import MutableArithmetics as MA
+import StarAlgebras as SA
 
 import MultivariatePolynomials as MP
 const _APL = MP.AbstractPolynomialLike

src/atomic.jl

@@ -46,17 +46,16 @@ function Base.show(io::IO, η::AtomicMeasure)
     end
 end
 
-measure(η::AtomicMeasure, x) = moment_vector(η, x)
-function Measure(η::AtomicMeasure{T}, x::AbstractVector{TT}) where {T,TT}
-    return Measure{T,MP.monomial_type(TT),MP.monomial_vector_type(x)}(η, x)
-end
-function Measure{T,MT,MVT}(η::AtomicMeasure{T}, x) where {T,MT,MVT}
-    X = MP.monomial_vector(x)
+function moment_vector(η::AtomicMeasure, basis::MB.SubBasis{MB.Monomial})
     return sum(
-        atom.weight * dirac(X, η.variables => atom.center) for atom in η.atoms
+        atom.weight * dirac(basis.monomials, η.variables => atom.center) for atom in η.atoms
     )
 end
 
+function moment_vector(η::AtomicMeasure, monos::AbstractVector{<:MP.AbstractTermLike})
+    return moment_vector(η, MB.SubBasis{MB.Monomial}(monos))
+end
+
 function expectation(η::AtomicMeasure, p::_APL)
     return sum(δ -> δ.weight * p(η.variables => δ.center), η.atoms)
 end

src/border.jl

@@ -121,21 +121,21 @@ function solve(
     mult = Matrix{T}[zeros(T, m, m) for _ in eachindex(vars)]
     completed_border = Dict{eltype(dependent.monomials),Vector{T}}()
     function known_border_coefficients(border)
-        return !isnothing(_index(standard, border)) ||
-               !isnothing(_index(dependent, border)) ||
+        return !isnothing(MB.monomial_index(standard, border)) ||
+               !isnothing(MB.monomial_index(dependent, border)) ||
                haskey(completed_border, border)
     end
     function border_coefficients(border)
-        k = _index(standard, border)
+        k = MB.monomial_index(standard, border)
         if !isnothing(k)
             return SparseArrays.sparsevec([k], [one(T)], m)
         end
-        k = _index(dependent, border)
+        k = MB.monomial_index(dependent, border)
         if !isnothing(k)
             v = zeros(T, m)
             row = 0
             for (i, std) in enumerate(standard.monomials)
-                j = _index(d.basis, std)
+                j = MB.monomial_index(d.basis, std)
                 if !is_trivial(d.dependence[j])
                     row += 1
                     v[i] = b.matrix[row, k]
@@ -243,7 +243,7 @@ function solve(
         new_basis, I1, I2 = MB.merge_bases(Ubasis, dependent)
         new_matrix = Matrix{T}(undef, length(new_basis), size(Uperp, 2))
         I_nontrivial_standard = [
-            _index(Ubasis, std) for
+            MB.monomial_index(Ubasis, std) for
             std in standard_basis(b.dependence, trivial = false).monomials
         ]
         Uperpstd = Uperp[I_nontrivial_standard, :]
@@ -322,7 +322,7 @@ function partial_commutation_fix(
                 continue
             end
             shifted = shift * standard.monomials[i]
-            j = _index(standard, shifted)
+            j = MB.monomial_index(standard, shifted)
             if !isnothing(j)
                 ret[j] += coef[i]
             elseif known_border_coefficients(shifted)
@@ -350,8 +350,8 @@ function partial_commutation_fix(
                     # but the other one is known ?
                     continue
                 end
-                if isnothing(_index(standard, mono_x))
-                    if isnothing(_index(standard, mono_y))
+                if isnothing(MB.monomial_index(standard, mono_x))
+                    if isnothing(MB.monomial_index(standard, mono_y))
                         coef_xy, unknowns_xy =
                             shifted_border_coefficients(mono_y, x)
                     else
@@ -367,7 +367,7 @@ function partial_commutation_fix(
                     coef_yx, unknowns_yx =
                         shifted_border_coefficients(mono_x, y)
                 else
-                    if !isnothing(_index(standard, mono_y))
+                    if !isnothing(MB.monomial_index(standard, mono_y))
                         # Let `f` be `known_border_coefficients`.
                         # They are both standard so we'll get
                         # `f(y * mono_x) - f(x * mono_y)`

src/dependence.jl

@@ -125,7 +125,7 @@ function category_markercolor(d::StaircaseDependence)
 end
 
 """
-    struct BasisDependence{D,B<:MB.AbstractPolynomialBasis}
+    struct BasisDependence{D,B<:SA.ExplicitBasis}
         basis::B
         dependence::Vector{D}
     end
@@ -135,7 +135,7 @@ The dependence between the elements of a basis.
 !!! tip
     If the number of variables is 2 or 3, it can be visualized with Plots.jl.
 """
-struct BasisDependence{D,B<:MB.AbstractPolynomialBasis}
+struct BasisDependence{D,B<:SA.ExplicitBasis}
     basis::B
     dependence::Vector{D}
 end
@@ -146,7 +146,7 @@ end
 
 _first_arg(cat, _) = cat
 
-function Base.show(io::IO, d::BasisDependence)
+function Base.show(io::IO, d::BasisDependence{D,<:MB.SubBasis{B}}) where {D,B}
     print(io, "BasisDependence for ")
     if isempty(d)
         println(io, "an empty basis")
@@ -155,7 +155,7 @@ function Base.show(io::IO, d::BasisDependence)
     end
     for (cat, monos) in basis_categories(d)
         println(io, " ", category_label(cat), ":")
-        println(io, " ", MB.MonomialBasis(monos))
+        println(io, " ", MB.SubBasis{B}(monos))
     end
     return
 end
@@ -197,8 +197,8 @@ end
 
 function BasisDependence{LinearDependence}(
     r,
-    basis::MB.MonomialBasis{M},
-) where {M}
+    basis::MB.SubBasis,
+)
     return BasisDependence(
         basis,
         LinearDependence[
@@ -222,7 +222,7 @@ end
 """
     function BasisDependence{StaircaseDependence}(
         is_dependent::Function,
-        basis::MB.AbstractPolynomialBasis,
+        basis::SA.ExplicitBasis,
     )
 
 Computes the set of standard monomials using the *greedy sieve* algorithm
@@ -236,13 +236,13 @@ Foundations of Computational Mathematics 8 (2008): 607-647.
 """
 function BasisDependence{StaircaseDependence}(
     r,
-    basis::MB.MonomialBasis{M},
+    basis::MB.SubBasis{MB.Monomial,M},
 ) where {M}
     if isempty(basis.monomials)
         return BasisDependence(basis, StaircaseDependence[])
     end
     function dependence(mono)
-        i = get(basis, mono, nothing)
+        i = MB.monomial_index(basis, mono)
         return if isnothing(i)
             TRIVIAL
         else
@@ -251,7 +251,7 @@ function BasisDependence{StaircaseDependence}(
     end
     vars = MP.variables(basis)
     full_basis =
-        MB.maxdegree_basis(typeof(basis), vars, MP.maxdegree(basis.monomials))
+        MB.maxdegree_basis(MB.Monomial, vars, MP.maxdegree(basis.monomials))
     d = StaircaseDependence[]
     # This sieve of [LLR08, Algorithm 1] is a performance improvement but not only.
     # It also ensures that the standard monomials have the "staircase structure".
@@ -273,7 +273,7 @@ function BasisDependence{StaircaseDependence}(
         end
     end
     full_basis = typeof(full_basis)(full_basis.monomials[keep])
-    new_basis = MB.MonomialBasis(
+    new_basis = MB.SubBasis{MB.Monomial}(
         eltype(basis.monomials)[
             full_basis.monomials[i] * shift for
             i in eachindex(d) if !is_dependent(d[i]) for shift in vars
@@ -289,7 +289,7 @@ function BasisDependence{StaircaseDependence}(
             else
                 # If it was not seen before, it means it is outside the basis
                 # so it is trivial standard
-                @assert isnothing(get(basis, mono, nothing))
+                @assert isnothing(MB.monomial_index(basis, mono))
                 std = true
             end
             deps[i] = StaircaseDependence(std, dependence(mono))

src/deprecate.jl

@@ -1,3 +1,2 @@
-@deprecate extractatoms atomic_measure
-@deprecate computesupport! compute_support!
-@deprecate getmat value_matrix
+@deprecate Measure MomentVector
+@deprecate measure moment_vector

src/expectation.jl

@@ -1,16 +1,16 @@
-function _expectation(μ::MomentVector{S}, p::_APL{T}, f) where {S,T}
-    i = 1
+function _expectation(μ::MomentVector{S,<:MB.SubBasis{MB.Monomial}}, p::_APL{T}, f) where {S,T}
+    i = firstindex(μ.basis)
     s = zero(MA.promote_operation(*, S, T))
     for t in MP.terms(p)
-        while i <= length(μ.x) && MP.monomial(t) != μ.x[i]
+        while i in eachindex(μ.basis) && MP.monomial(t) != μ.basis.monomials[i]
             i += 1
         end
-        if i > length(μ.x)
+        if !(i in eachindex(μ.basis))
             error(
-                "The polynomial $p has a nonzero term $t with monomial $(t.x) for which the expectation is not known in $μ",
+                "The polynomial $p has a nonzero term $t with monomial $(MP.monomial(t)) for which the expectation is not known in $μ",
             )
         end
-        s += f(μ.a[i], MP.coefficient(t))
+        s += f(μ.values[i], MP.coefficient(t))
         i += 1
     end
     return s

src/extract.jl

@@ -47,11 +47,11 @@ end
     end
 
 Given a moment matrix `ν` and the atom centers, first convert the moment matrix
-to a vector of moments, using [`measure(ν; rtol=rtol, atol=atol)`](@ref measure)
+to a vector of moments, using [`moment_vector(ν; rtol=rtol, atol=atol)`](@ref moment_vector)
 and then determine the weights by solving a linear system over the monomials obtained.
 
 If the moment values corresponding to the same monomials can have small differences,
-[`measure`](@ref) can throw an error if `rtol` and `atol` are not small enough.
+[`moment_vector`](@ref) can throw an error if `rtol` and `atol` are not small enough.
 Alternatively to tuning these tolerances [`MomentVectorWeightSolver`](@ref) can be used instead.
 """
 struct MomentVectorWeightSolver{T}
@@ -82,17 +82,17 @@ function MomentVectorWeightSolver(; rtol = nothing, atol = nothing)
 end
 
 function solve_weight(
-    ν::MomentMatrix{T},
+    ν::MomentMatrix{T,<:MB.SubBasis{MB.Monomial}},
     centers,
     solver::MomentVectorWeightSolver,
 ) where {T}
-    μ = measure(ν; rtol = solver.rtol, atol = solver.atol)
+    μ = moment_vector(ν; rtol = solver.rtol, atol = solver.atol)
     vars = MP.variables(μ)
-    A = Matrix{T}(undef, length(μ.x), length(centers))
+    A = Matrix{T}(undef, length(μ.basis), length(centers))
     for i in eachindex(centers)
-        A[:, i] = dirac(μ.x, vars => centers[i]).a
+        A[:, i] = dirac(μ.basis.monomials, vars => centers[i]).values
     end
-    return A \ μ.a
+    return A \ μ.values
 end
 
 """

src/flat.jl

@@ -100,6 +100,6 @@ function compute_support!(
     rank_check::RankCheck,
     solver::FlatExtension,
 )
-    ν.support = support(measure(ν), rank_check, solver)
+    ν.support = support(moment_vector(ν), rank_check, solver)
     return
 end

src/moment.jl

@@ -1,14 +1,14 @@
-struct Moment{T,MT<:MP.AbstractMonomial} <: AbstractMoment{T}
+struct Moment{T,P} <: AbstractMoment{T}
     α::T
-    x::MT
+    polynomial::P
 end
 
 """
-    moment(α, m::AbstractMonomial)
+    moment(α, p)
 
-Creates the moment of the monomial `m` of value `α`.
+Creates the moment of the polynomial `p` of value `α`.
 """
-moment(α, m::MP.AbstractMonomial) = Moment(α, m)
+moment(α, p) = Moment(α, p)
 
 """
     moment_value(m::AbstractMomentLike)
@@ -21,19 +21,8 @@ Calling `moment_value(moment(3.1, x*y^2))` should return `3.1`.
 """
 moment_value(m::Moment) = m.α
 
-"""
-    monomial(m::AbstractMomentLike)
-
-Returns the monomial of the moment `m`.
-
-## Examples
-
-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]
+for f in [:variables, :nvariables]
     @eval begin
-        MP.$f(m::AbstractMomentLike) = MP.$f(MP.monomial(m))
+        MP.$f(m::AbstractMomentLike) = MP.$f(m.polynomial)
     end
 end