From c3180e1fd296000b01c6d7a33df381dc900ff617 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 17 Jun 2024 16:49:07 +0200 Subject: [PATCH] Support StarAlgebras in expectation (#81) * Support StarAlgebras in expectation * Fix format * Fixes * Fix format * Remove unused --- .github/workflows/ci.yml | 2 +- .github/workflows/documentation.yml | 2 +- src/MultivariateMoments.jl | 2 +- src/border.jl | 5 ++-- src/expectation.jl | 40 ++++++++++++++++++++++------- test/expectation.jl | 15 ++++++++--- 6 files changed, 49 insertions(+), 17 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 6c73b17..99abe63 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -47,7 +47,7 @@ jobs: run: | using Pkg Pkg.add([ - PackageSpec(name="StarAlgebras", rev="mk/non_monomial_basis"), + PackageSpec(name="StarAlgebras", rev="main"), PackageSpec(name="MultivariateBases", rev="master"), ]) - uses: julia-actions/julia-buildpkg@v1 diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml index 090c4d7..e720daf 100644 --- a/.github/workflows/documentation.yml +++ b/.github/workflows/documentation.yml @@ -19,7 +19,7 @@ jobs: run: | using Pkg Pkg.add([ - PackageSpec(name="StarAlgebras", rev="mk/non_monomial_basis"), + PackageSpec(name="StarAlgebras", rev="main"), PackageSpec(name="MultivariateBases", rev="master"), PackageSpec(path=pwd()), ]) diff --git a/src/MultivariateMoments.jl b/src/MultivariateMoments.jl index 4d4bb29..165b0ec 100644 --- a/src/MultivariateMoments.jl +++ b/src/MultivariateMoments.jl @@ -7,7 +7,7 @@ import MutableArithmetics as MA import StarAlgebras as SA import MultivariatePolynomials as MP -const _APL = MP.AbstractPolynomialLike +const _APL = Union{MP.AbstractPolynomialLike,SA.AlgebraElement} import MultivariateBases as MB diff --git a/src/border.jl b/src/border.jl index 2535d5a..76971bd 100644 --- a/src/border.jl +++ b/src/border.jl @@ -468,8 +468,9 @@ function solve(b::BorderBasis{E}, solver::AlgebraicBorderSolver{D}) where {D,E} # 2005 ind = independent_basis(b) dep = dependent_basis(b.dependence) - system = [ - dep.monomials[col] - MP.polynomial(b.matrix[:, col], ind) for + system = MP.polynomial_type(ind, eltype(b.matrix))[ + dep.monomials[col] - + MP.polynomial(MB.algebra_element(b.matrix[:, col], ind)) for col in eachindex(dep.monomials) ] filter!(!MP.isconstant, system) diff --git a/src/expectation.jl b/src/expectation.jl index d7a1ac6..06e0c11 100644 --- a/src/expectation.jl +++ b/src/expectation.jl @@ -1,25 +1,37 @@ function _expectation( - μ::MomentVector{S,<:MB.SubBasis{MB.Monomial}}, - p::_APL{T}, + μ::MomentVector{S,<:MB.SubBasis{B}}, + p::SA.AlgebraElement{<:MB.Algebra{BT,B},T}, f, -) where {S,T} +) where {S,T,BT,B} i = firstindex(μ.basis) s = zero(MA.promote_operation(*, S, T)) - for t in MP.terms(p) - while i in eachindex(μ.basis) && MP.monomial(t) != μ.basis.monomials[i] + for (k, v) in SA.nonzero_pairs(SA.coeffs(p)) + mono = SA.basis(p)[k] + while i in eachindex(μ.basis) && mono != μ.basis[i] i += 1 end if !(i in eachindex(μ.basis)) error( - "The polynomial $p has a nonzero term $t with monomial $(MP.monomial(t)) for which the expectation is not known in $μ", + "The polynomial $p has a nonzero term $mono with coefficient $v for which the expectation is not known in $μ", ) end - s += f(μ.values[i], MP.coefficient(t)) + s += f(μ.values[i], v) i += 1 end return s end +function _expectation(μ::MomentVector, p::MP.AbstractPolynomialLike, f) + return _expectation( + μ, + MB.algebra_element( + MB.sparse_coefficients(MP.polynomial(p)), + MB.FullBasis{MB.Monomial,MP.monomial_type(p)}(), + ), + f, + ) +end + """ MultivariateMoments.expectation(μ::AbstractMeasureLike, p::AbstractPolynomialLike) MultivariateMoments.expectation(p::AbstractPolynomialLike, μ::AbstractMeasureLike) @@ -37,5 +49,15 @@ expectation(p::_APL, μ::MomentVector) = _expectation(μ, p, (a, b) -> b * a) # See [`expectation`](@ref) """ -LinearAlgebra.dot(μ::AbstractMeasureLike, p::_APL) = expectation(μ, p) -LinearAlgebra.dot(p::_APL, μ::AbstractMeasureLike) = expectation(p, μ) +LinearAlgebra.dot(μ::AbstractMeasureLike, p::MP.AbstractPolynomialLike) = + expectation(μ, p) +function LinearAlgebra.dot(p::MP.AbstractPolynomialLike, μ::AbstractMeasureLike) + return expectation(p, μ) +end +# Need to split it and not use `_APL` to avoid ambiguity +function LinearAlgebra.dot(μ::AbstractMeasureLike, p::SA.AlgebraElement) + return expectation(μ, p) +end +function LinearAlgebra.dot(p::SA.AlgebraElement, μ::AbstractMeasureLike) + return expectation(p, μ) +end diff --git a/test/expectation.jl b/test/expectation.jl index cdc8e73..33d397d 100644 --- a/test/expectation.jl +++ b/test/expectation.jl @@ -3,10 +3,19 @@ p = x[3] - 2x[1] * x[2]^2 + 3x[3] * x[1] - 5x[1]^3 v = (1, 2, 3) m = dirac(monomials(p), x => v) - @test (@inferred MultivariateMoments.expectation(m, p)) == - p(x => v) == - (@inferred MultivariateMoments.expectation(p, m)) + a = MB.algebra_element(p) + expected = p(x => v) + @test (@inferred MultivariateMoments.expectation(m, p)) == expected + @test (@inferred MultivariateMoments.expectation(m, a)) == expected + @test (@inferred MultivariateMoments.expectation(p, m)) == expected + @test (@inferred MultivariateMoments.expectation(a, m)) == expected + @test (@inferred dot(m, p)) == expected + @test (@inferred dot(m, a)) == expected + @test (@inferred dot(p, m)) == expected + @test (@inferred dot(a, m)) == expected @test_throws ErrorException dot(x[1] * x[2] * x[3], m) @test (@inferred dot(0.5 * x[1] * x[2]^2, m)) == 2.0 + a = MB.algebra_element([0.5], MB.SubBasis{MB.Monomial}([x[1] * x[2]^2])) + @test (@inferred dot(a, m)) == 2.0 @test (@inferred dot(m, x[1] * x[3])) == 3 end