diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 0a823c4..6c73b17 100644 --- a/.github/workflows/ci.yml +++ b/.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: diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index d5291d5..090c4d7 100644 --- a/.github/workflows/documentation.yml +++ b/.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 diff --git a/docs/Project.toml b/docs/Project.toml index dd849ca..f7414f9 100644 --- a/docs/Project.toml +++ b/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" diff --git a/docs/src/moments.md b/docs/src/moments.md index 6e15743..06331e0 100644 --- a/docs/src/moments.md +++ b/docs/src/moments.md @@ -12,16 +12,15 @@ 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. @@ -29,7 +28,6 @@ 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) diff --git a/src/MultivariateMoments.jl b/src/MultivariateMoments.jl index 90c1c1c..4d4bb29 100644 --- a/src/MultivariateMoments.jl +++ b/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 diff --git a/src/atomic.jl b/src/atomic.jl index bc8d112..2cbdcda 100644 --- a/src/atomic.jl +++ b/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 diff --git a/src/border.jl b/src/border.jl index 8da46ed..1781937 100644 --- a/src/border.jl +++ b/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)` diff --git a/src/dependence.jl b/src/dependence.jl index c705b17..d124dd9 100644 --- a/src/dependence.jl +++ b/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)) diff --git a/src/deprecate.jl b/src/deprecate.jl index 68a8a5b..26ab10d 100644 --- a/src/deprecate.jl +++ b/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 diff --git a/src/expectation.jl b/src/expectation.jl index 520e57b..c105e81 100644 --- a/src/expectation.jl +++ b/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 diff --git a/src/extract.jl b/src/extract.jl index 68c9f68..8816848 100644 --- a/src/extract.jl +++ b/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 """ diff --git a/src/flat.jl b/src/flat.jl index 2d994a9..af405a4 100644 --- a/src/flat.jl +++ b/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 diff --git a/src/moment.jl b/src/moment.jl index d84bb8a..2cb4340 100644 --- a/src/moment.jl +++ b/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 diff --git a/src/moment_matrix.jl b/src/moment_matrix.jl index 1cf209b..a123a10 100644 --- a/src/moment_matrix.jl +++ b/src/moment_matrix.jl @@ -1,10 +1,10 @@ using SemialgebraicSets -abstract type AbstractMomentMatrix{T,B<:MB.AbstractPolynomialBasis} <: +abstract type AbstractMomentMatrix{T,B<:SA.ExplicitBasis} <: AbstractMeasureLike{T} end """ - mutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis, MT<:AbstractMatrix{T}} <: AbstractMeasureLike{T} + mutable struct MomentMatrix{T,B<:SA.ExplicitBasis,MT<:AbstractMatrix{T}} <: AbstractMeasureLike{T} Q::MT basis::B support::Union{Nothing, AlgebraicSet} @@ -15,7 +15,7 @@ The set of points that are zeros of all the polynomials ``p`` such that ``\\math """ mutable struct MomentMatrix{ T, - B<:MB.AbstractPolynomialBasis, + B<:SA.ExplicitBasis, MT<:AbstractMatrix{T}, } <: AbstractMomentMatrix{T,B} Q::MT @@ -24,19 +24,19 @@ mutable struct MomentMatrix{ end function MomentMatrix{T,B,MT}( Q::MT, - basis::MB.AbstractPolynomialBasis, + basis::SA.ExplicitBasis, ) where {T,B,MT} return MomentMatrix{T,B,MT}(Q, basis, nothing) end function MomentMatrix{T,B}( Q::AbstractMatrix{T}, - basis::MB.AbstractPolynomialBasis, + basis::SA.ExplicitBasis, ) where {T,B} return MomentMatrix{T,B,typeof(Q)}(Q, basis) end function MomentMatrix( Q::SymMatrix{T}, - basis::MB.AbstractPolynomialBasis, + basis::SA.ExplicitBasis, ) where {T} return MomentMatrix{T,typeof(basis)}(Q, basis) end @@ -46,7 +46,7 @@ MP.nvariables(μ::MomentMatrix) = MP.nvariables(μ.basis) function MomentMatrix{T}( f::Function, - basis::MB.AbstractPolynomialBasis, + basis::SA.ExplicitBasis, σ = 1:length(basis), ) where {T} return MomentMatrix( @@ -56,7 +56,7 @@ function MomentMatrix{T}( end function MomentMatrix{T}(f::Function, monos::AbstractVector) where {T} σ, sorted_monos = MP.sort_monomial_vector(monos) - return MomentMatrix{T}(f, MB.MonomialBasis(sorted_monos), σ) + return MomentMatrix{T}(f, MB.SubBasis{MB.Monomial}(sorted_monos), σ) end function show_basis_indexed_matrix(io::IO, A, pre = "") @@ -79,39 +79,39 @@ end Creates a matrix the moment matrix for the moment matrix ``x x^\\top`` using the moments of `μ`. """ -function moment_matrix(μ::MomentVector{T}, X) where {T} - return MomentMatrix{T}((i, j) -> moment_value(μ, X[i] * X[j]), X) +function moment_matrix(μ::MomentVector{T}, basis) where {T} + return MomentMatrix{T}((i, j) -> moment_value(μ, basis[i] * basis[j]), basis) end function MomentMatrix( Q::AbstractMatrix{T}, - basis::MB.AbstractPolynomialBasis, + basis::SA.ExplicitBasis, σ, ) where {T} return MomentMatrix{T}((i, j) -> Q[σ[i], σ[j]], basis) end function MomentMatrix(Q::AbstractMatrix, monos::AbstractVector) σ, sorted_monos = MP.sort_monomial_vector(monos) - return MomentMatrix(Q, MB.MonomialBasis(sorted_monos), σ) + return MomentMatrix(Q, MB.SubBasis{MB.Monomial}(sorted_monos), σ) end moment_matrix(Q::AbstractMatrix, monos) = MomentMatrix(Q, monos) value_matrix(μ::MomentMatrix) = Matrix(μ.Q) -function vectorized_basis(ν::MomentMatrix{T,<:MB.MonomialBasis}) where {T} +function vectorized_basis(ν::MomentMatrix{T,<:MB.SubBasis{MB.Monomial}}) where {T} monos = ν.basis.monomials n = length(monos) - # We don't wrap in `MonomialBasis` as we don't want the monomials + # We don't wrap in `MB.SubBasis` as we don't want the monomials # to be `sort`ed and `uniq`ed. return [monos[i] * monos[j] for j in 1:n for i in 1:j] end -function measure(ν::MomentMatrix; kws...) +function moment_vector(ν::MomentMatrix; kws...) n = length(ν.basis) - return measure(ν.Q.Q, vectorized_basis(ν); kws...) + return moment_vector(ν.Q.Q, vectorized_basis(ν); kws...) end -struct BlockDiagonalMomentMatrix{T,B<:MB.AbstractPolynomialBasis,MT} <: +struct BlockDiagonalMomentMatrix{T,B<:SA.ExplicitBasis,MT} <: AbstractMomentMatrix{T,B} blocks::Vector{MomentMatrix{T,B,MT}} end diff --git a/src/moment_vector.jl b/src/moment_vector.jl index f2bc939..80981cb 100644 --- a/src/moment_vector.jl +++ b/src/moment_vector.jl @@ -1,25 +1,24 @@ function _check_length(values, basis) - if SA.basis(parent) <: SA.ExplicitBasis && length(values) != length(basis) + if length(values) != length(basis) throw(DimensionMismatch("dimension must match: `values` has length `$(length(values))` and `basis` has length `$(length(basis))`")) end end # If a monomial is not in x, it does not mean that the moment is zero, it means that it is unknown/undefined -struct MomentVector{T,A,V} <: AbstractMomentArray{T,A} +struct MomentVector{T,B<:SA.AbstractBasis,V} <: AbstractMeasure{T} values::V - parent::A + basis::B - function MomentVector{T,A,V}(values::V, parent::A) where {T,A,V} - if SA.basis(parent) <: SA.ExplicitBasis - _check_length(values, SA.basis(parent)) + function MomentVector{T,B,V}(values::V, basis::B) where {T,B,V} + if basis isa SA.ExplicitBasis + _check_length(values, basis) end - return new{T,A,V}(values, parent) + return new{T,B,V}(values, basis) end end -function moment_vector(values::AbstractVector{T}, basis::MB.SubBasis) where {T} - parent = MB.algebra(basis) - return MomentVector{T,typeof(parent),typeof(values)}(values, parent) +function moment_vector(values::AbstractVector{T}, basis::SA.AbstractBasis) where {T} + return MomentVector{T,typeof(basis),typeof(values)}(values, basis) end """ @@ -52,68 +51,67 @@ function moment_vector( end end end - return moment_vector(values, MB.SubBasis{MB.Monomial}(sorted_monos)) + return moment_vector(sorted_values, MB.SubBasis{MB.Monomial}(sorted_monos)) end -""" - variables(μ::AbstractMeasureLike) +SA.basis(μ::MomentVector) = μ.basis -Returns the variables of `μ` in decreasing order. Just like in MultivariatePolynomials, it could contain variables of zero degree in every monomial. """ -MP.variables(μ::MomentVector) = MP.variables(μ.x) + variables(μ::MomentVector) +Returns the variables of `μ` in decreasing order. Just like in MultivariatePolynomials, it could contain variables of zero degree in every monomial. """ - monomials(μ::AbstractMeasureLike) - -Returns an iterator over the monomials of `μ` sorted in the decreasing order. -""" -MP.monomials(μ::MomentVector) = μ.x +MP.variables(μ::MomentVector) = MP.variables(SA.basis(μ)) """ maxdegree(μ::AbstractMeasureLike) Returns the maximal degree of the monomials of `μ`. """ -MP.maxdegree(μ::MomentVector) = MP.maxdegree(MP.monomials(μ)) +MP.maxdegree(μ::MomentVector) = MP.maxdegree(SA.basis(μ)) """ - mindegree(μ::AbstractMeasureLike) + mindegree(μ::MomentVector) Returns the minimal degree of the monomials of `μ`. """ -MP.mindegree(μ::MomentVector) = MP.mindegree(MP.monomials(μ)) +MP.mindegree(μ::MomentVector) = MP.mindegree(SA.basis(μ)) """ - extdegree(μ::AbstractMeasureLike) + extdegree(μ::MomentVector) Returns the extremal degrees of the monomials of `μ`. """ -MP.extdegree(μ::MomentVector) = MP.extdegree(MP.monomials(μ)) +MP.extdegree(μ::MomentVector) = MP.extdegree(SA.basis(μ)) """ - moments(μ::AbstractMeasureLike) + moments(μ::MomentVector) Returns an iterator over the moments of `μ` sorted in decreasing order of monomial. """ -moments(μ::MomentVector) = map((α, x) -> moment(α, x), μ.a, μ.x) +moments(μ::MomentVector) = map((α, x) -> moment(α, x), μ.values, SA.basis(μ)) -Base.:(*)(α, μ::MomentVector) = measure(α * μ.a, μ.x) -Base.:(*)(μ::MomentVector, α) = measure(μ.a * α, μ.x) -Base.:(-)(μ::MomentVector) = measure(-μ.a, μ.x) +Base.:(*)(α, μ::MomentVector) = moment_vector(α * μ.values, SA.basis(μ)) +Base.:(*)(μ::MomentVector, α) = moment_vector(μ.values * α, SA.basis(μ)) +Base.:(-)(μ::MomentVector) = moment_vector(-μ.values, SA.basis(μ)) function Base.:(+)(μ::MomentVector, ν::MomentVector) - @assert μ.x == ν.x - return measure(μ.a + ν.a, μ.x) -end -function _index(basis::MB.SubBasis{B}, mono) where {B} - return get(basis, MB.Polynomial{B}(mono), nothing) + @assert SA.basis(μ) == SA.basis(ν) + return moment_vector(μ.values + ν.values, SA.basis(μ)) end -function moment_value(μ, mono) - i = _index(SA.basis(μ.parent), mono) +moment_value(μ::MomentVector, t::MP.AbstractTerm) = MP.coefficient(t) * moment_value(μ, MP.monomial(t)) +moment_value(μ::MomentVector, mono::MP.AbstractMonomial) = moment_value(μ, MB.Polynomial{MB.Monomial}(mono)) + +function moment_value(μ::MomentVector{T,<:MB.SubBasis{B}}, p::MB.Polynomial{B}) where {T,B} + i = MB.monomial_index(SA.basis(μ), p.monomial) if isnothing(i) - throw(ArgumentError("`$μ` does not have the moment `$mono`")) + throw(ArgumentError("`$μ` does not have the moment `$p`")) end - return μ.a[i] + return μ.values[i] +end + +function moment_value(μ::MomentVector, p::SA.AlgebraElement) + return sum(coef * moment_value(μ, SA.basis(parent(p))[mono]) for (mono, coef) in SA.nonzero_pairs(p.coeffs)) end """ @@ -129,5 +127,5 @@ function dirac( x::AbstractVector{MT}, s::MP.AbstractSubstitution..., ) where {MT<:MP.AbstractMonomial} - return Measure([m(s...) for m in x], x) + return moment_vector([m(s...) for m in x], x) end diff --git a/src/null.jl b/src/null.jl index 27a9632..ac4e684 100644 --- a/src/null.jl +++ b/src/null.jl @@ -11,7 +11,7 @@ the image space of a moment matrix with rows and columns indexed by `basis`). The value of `matrix[i, j]` should be interpreted as the value of the `i`th element of `basis` for the `j`th generator of the null space (resp. image) space. """ -struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT<:MB.AbstractPolynomialBasis} +struct MacaulayNullspace{T,MT<:AbstractMatrix{T},BT<:SA.ExplicitBasis} matrix::MT basis::BT accuracy::T diff --git a/src/shift.jl b/src/shift.jl index 5c3034e..d69712e 100644 --- a/src/shift.jl +++ b/src/shift.jl @@ -46,7 +46,7 @@ function BasisDependence{StaircaseDependence}( end function _indices(in::MB.SubBasis, from::MB.SubBasis) - return Int[_index(in, mono) for mono in from.monomials] + return Int[MB.monomial_index(in, mono) for mono in from.monomials] end function BorderBasis( diff --git a/test/atomic.jl b/test/atomic.jl index 6a85395..28c2e7c 100644 --- a/test/atomic.jl +++ b/test/atomic.jl @@ -18,6 +18,6 @@ p = x^2 + x + y + x * y^2 + 1 @test dot(η, p) == 2 * p(x => 1, y => 0) + 3 * p(x => 1 / 2, y => 1 / 2) @test dot(η, p) == dot(p, η) - μ = measure(η, monomials(p)) + μ = moment_vector(η, monomials(p)) @test dot(μ, p) == 2 * p(x => 1, y => 0) + 3 * p(x => 1 / 2, y => 1 / 2) end diff --git a/test/commutativetests.jl b/test/commutativetests.jl index a8d868b..47a6b9b 100644 --- a/test/commutativetests.jl +++ b/test/commutativetests.jl @@ -1,4 +1,4 @@ -include("measure.jl") +include("moment_vector.jl") include("expectation.jl") include("moment_matrix.jl") include("rank.jl") diff --git a/test/dependence.jl b/test/dependence.jl index 9d66fbe..b3851e1 100644 --- a/test/dependence.jl +++ b/test/dependence.jl @@ -124,17 +124,17 @@ function test_recipe(x, y, z) @test sprint(show, dep) == """ BasisDependence for bases: Standard: - MonomialBasis([z]) + SubBasis{Monomial}([z]) Trivial Standard: - MonomialBasis([1, y, z^2, y*z, z^3, y*z^2]) + SubBasis{Monomial}([1, y, z^2, y*z, z^3, y*z^2]) Corners: - MonomialBasis([x, y^2]) + SubBasis{Monomial}([x, y^2]) Independent Border: - MonomialBasis([x*z]) + SubBasis{Monomial}([x*z]) Trivial Independent Border: - MonomialBasis([y^2*z, x*z^2, x*y*z]) + SubBasis{Monomial}([y^2*z, x*z^2, x*y*z]) Dependent Border: - MonomialBasis([x*y]) + SubBasis{Monomial}([x*y]) """ @test MM.corners_basis(dep).monomials == c _test_recipe( diff --git a/test/extract.jl b/test/extract.jl index f6634a1..4cbb2a9 100644 --- a/test/extract.jl +++ b/test/extract.jl @@ -29,7 +29,7 @@ function _atoms(atoms, rank_check, solver) Mod.@polyvar x[1:2] η = AtomicMeasure(x, WeightedDiracMeasure.(atoms, ones(length(atoms)))) monos = monomials(x, 0:(length(atoms)+2)) - μ = measure(η, monos) + μ = moment_vector(η, monos) ν = moment_matrix(μ, monomials(x, 0:(div(length(atoms), 2)+1))) atoms = atomic_measure(ν, rank_check, solver) @test atoms !== nothing @@ -80,7 +80,7 @@ function hl05_2_3(rank_check, lrc, solver, perturb::Bool = true) y, 1, ] - μ = measure(η, monos) + μ = moment_vector(η, monos) ν = moment_matrix(μ, [1, x, y, x^2, x * y, y^2]) atoms = atomic_measure(ν, rank_check, lrc, Echelon(), solver) @test atoms !== nothing @@ -173,7 +173,7 @@ function hl05_4(rank_check, lrc) [0.25, 0.25, 0.25, 0.25], ), ) - μ = measure( + μ = moment_vector( [1 / 9, 0, 1 / 9, 0, 1 / 9, 0, 0, 0, 0, 1 / 3, 0, 1 / 3, 0, 0, 1], [ x^4, @@ -199,7 +199,7 @@ function hl05_4(rank_check, lrc) if lrc isa LowRankLDLTAlgorithm @test atoms ≈ η end - @test measure(atoms, μ.x).a ≈ μ.a rtol = 1e-3 + @test moment_vector(atoms, μ.basis).values ≈ μ.values rtol = 1e-3 end """ @@ -210,7 +210,7 @@ end function lpj20_3_8_0(rank_check, lrc, ok::Bool = true) Mod.@polyvar x[1:2] η = AtomicMeasure(x, [WeightedDiracMeasure([1.0, 0.0], 2.0)]) - μ = measure([1e-6, 0.0, 2], monomials(x, 2)) + μ = moment_vector([1e-6, 0.0, 2], monomials(x, 2)) ν = moment_matrix(μ, monomials(x, 1)) atoms = atomic_measure(ν, rank_check, lrc) if ok diff --git a/test/moment_matrix.jl b/test/moment_matrix.jl index 7a6989a..e02b119 100644 --- a/test/moment_matrix.jl +++ b/test/moment_matrix.jl @@ -2,17 +2,17 @@ using Test @testset "MomentMatrix" begin Mod.@polyvar x y - μ = measure([1], [x]) + μ = moment_vector([1], [x]) for mono in [x^0, x, y] - err = ArgumentError("`$μ` does not have the moment `$(mono^2)`") + err = ArgumentError("`$μ` does not have the moment `$(MB.Polynomial{MB.Monomial}(mono^2))`") @test_throws err moment_matrix(μ, [mono]) end - μ = measure(1:3, [x^2, x * y, y^2]) + μ = moment_vector(1:3, [x^2, x * y, y^2]) X = [x, y] ν1 = moment_matrix(μ, X) @test sprint(show, ν1) == """ MomentMatrix with row/column basis: - MonomialBasis([y, x]) + SubBasis{Monomial}([y, x]) And entries in a 2×2 SymMatrix{$Int}: 3 2 2 1""" @@ -27,18 +27,18 @@ And entries in a 2×2 SymMatrix{$Int}: @test variables(ν)[2] == y @test nvariables(ν) == 2 end - @test_throws ArgumentError moment_matrix(measure([1], [x]), [y]) + @test_throws ArgumentError moment_matrix(moment_vector([1], [x]), [y]) block_ν = BlockDiagonalMomentMatrix([ν1, ν2]) @test block_ν isa BlockDiagonalMomentMatrix{Int,typeof(ν1.basis)} @test sprint(show, block_ν) == """ BlockDiagonalMomentMatrix with 2 blocks: [1] Block with row/column basis: - MonomialBasis([y, x]) + SubBasis{Monomial}([y, x]) And entries in a 2×2 SymMatrix{$Int}: 3 2 2 1 [2] Block with row/column basis: - MonomialBasis([y, x]) + SubBasis{Monomial}([y, x]) And entries in a 2×2 SymMatrix{$Int}: 3 2 2 1""" diff --git a/test/moment_vector.jl b/test/moment_vector.jl index 03ee58c..a017833 100644 --- a/test/moment_vector.jl +++ b/test/moment_vector.jl @@ -1,28 +1,29 @@ -@testset "Measure" begin +import StarAlgebras as SA + +@testset "MomentVector" begin Mod.@polyvar x y - @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_throws DimensionMismatch moment_vector([1, 2], [x, x * y, y]) + @test_throws DimensionMismatch moment_vector([1, 2, 3, 4], [x, x * y, y]) + μ = moment_vector([1, 0, 2, 3], [x^2 * y^2, y * x^2, x * y * x^2, x * y^2]) @test MP.mindegree(μ) == 3 @test MP.maxdegree(μ) == 4 @test MP.extdegree(μ) == (3, 4) - @test monomials(μ) == + @test SA.basis(μ).monomials == monomial_vector([x^3 * y, x^2 * y^2, x^2 * y, x * y^2]) - @test monomial.(moments(μ)) == + @test map(m -> m.polynomial.monomial, moments(μ)) == monomial_vector([x^3 * y, x^2 * y^2, x^2 * y, x * y^2]) @test MultivariateMoments.moment_value.(moments(μ)) == reverse([2, 1, 0, 3]) + @test all(m -> variables(m) == variables(x * y), moments(μ)) @test all(nvariables.(moments(μ)) .== 2) - @test degree.(moments(μ)) == reverse([4, 4, 3, 3]) - @test μ.a == reverse([2, 1, 0, 3]) - @test (-μ).a == reverse([-2, -1, 0, -3]) - @test (2 * μ).a == reverse([4, 2, 0, 6]) - @test (μ * 3).a == reverse([6, 3, 0, 9]) - μ = measure([1, 1], [x, x]) - @test monomials(μ) == [x] + @test μ.values == reverse([2, 1, 0, 3]) + @test (-μ).values == reverse([-2, -1, 0, -3]) + @test (2 * μ).values == reverse([4, 2, 0, 6]) + @test (μ * 3).values == reverse([6, 3, 0, 9]) + μ = moment_vector([1, 1], [x, x]) @test MultivariateMoments.moment_value.(moments(μ)) == [1] err = ErrorException( "The monomial `x` occurs twice with different values: `1` and `2`", ) - @test_throws err measure([1, 2], [x, x]) - #@test_throws ArgumentError measure([1], [x]) + measure([1], [y]) + @test_throws err moment_vector([1, 2], [x, x]) + #@test_throws ArgumentError moment_vector([1], [x]) + moment_vector([1], [y]) end