Support StarAlgebras in expectation (#81)

* Support StarAlgebras in expectation * Fix format * Fixes * Fix format * Remove unused
JuliaAlgebra · Jun 17, 2024 · c3180e1 · c3180e1
1 parent 71944e4
commit c3180e1
Show file tree

Hide file tree

Showing 6 changed files with 49 additions and 17 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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