Implement StarAlgebras (#23)

* Implement StarAlgebras

* Fix format

* Fixes

* Fixes

* Fix doc build

.github/workflows/ci.yml

@@ -42,6 +42,13 @@ 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"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:

.github/workflows/documentation.yml

@@ -15,7 +15,14 @@ jobs:
           # Build documentation on latest 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(path=pwd()),
+          ])
+          Pkg.instantiate()
       - name: Build and deploy
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token

Project.toml

@@ -7,6 +7,7 @@ version = "0.2.2"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 
 [compat]
 MultivariatePolynomials = "0.5.3"

docs/src/index.md

@@ -14,31 +14,29 @@ based on the [MultivariatePolynomials](https://github.com/JuliaAlgebra/Multivari
 AbstractPolynomialBasis
 maxdegree_basis
 basis_covering_monomials
-FixedPolynomialBasis
-OrthonormalCoefficientsBasis
 ```
 
 ## Monomial basis
 
 ```@docs
-MonomialBasis
-ScaledMonomialBasis
+Monomial
+ScaledMonomial
 ```
 
 ## Orthogonal basis
 
 ```@docs
-AbstractMultipleOrthogonalBasis
+AbstractMultipleOrthogonal
 univariate_orthogonal_basis
 reccurence_first_coef
 reccurence_second_coef
 reccurence_third_coef
 reccurence_deno_coef
-ProbabilistsHermiteBasis
-PhysicistsHermiteBasis
-LaguerreBasis
-AbstractGegenbauerBasis
-LegendreBasis
-ChebyshevBasisFirstKind
-ChebyshevBasisSecondKind
+ProbabilistsHermite
+PhysicistsHermite
+Laguerre
+AbstractGegenbauer
+Legendre
+ChebyshevFirstKind
+ChebyshevSecondKind
 ```

src/MultivariateBases.jl

@@ -1,38 +1,32 @@
 module MultivariateBases
 
 import MutableArithmetics as MA
-
+import StarAlgebras as SA
 import MultivariatePolynomials as MP
 
-export AbstractPolynomialBasis
+export AbstractPolynomialBasis, FullBasis, SubBasis
 export maxdegree_basis, basis_covering_monomials, empty_basis
 include("interface.jl")
 
 export AbstractMonomialBasis, MonomialBasis, ScaledMonomialBasis
+include("polynomial.jl")
 include("monomial.jl")
 include("scaled.jl")
 
-export FixedPolynomialBasis,
-    OrthonormalCoefficientsBasis,
-    AbstractMultipleOrthogonalBasis,
-    ProbabilistsHermiteBasis,
-    PhysicistsHermiteBasis,
-    LaguerreBasis
-export AbstractGegenbauerBasis,
-    LegendreBasis,
-    ChebyshevBasis,
-    ChebyshevBasisFirstKind,
-    ChebyshevBasisSecondKind
+export AbstractMultipleOrthogonal,
+    ProbabilistsHermite, PhysicistsHermite, Laguerre
+export AbstractGegenbauer,
+    Legendre, Chebyshev, ChebyshevFirstKind, ChebyshevSecondKind
 export generators,
     univariate_orthogonal_basis,
     reccurence_first_coef,
     reccurence_second_coef,
     reccurence_third_coef,
     reccurence_deno_coef
-include("fixed.jl")
+#include("fixed.jl")
 
 import LinearAlgebra
-include("orthonormal.jl")
+#include("orthonormal.jl")
 include("orthogonal.jl")
 include("hermite.jl")
 include("laguerre.jl")

src/chebyshev.jl

@@ -1,34 +1,71 @@
-abstract type AbstractChebyshevBasis{P} <: AbstractGegenbauerBasis{P} end
+abstract type AbstractChebyshev <: AbstractGegenbauer end
 
-function MP.polynomial_type(::Type{<:AbstractChebyshevBasis}, V::Type)
-    return MP.polynomial_type(V, Float64)
+_promote_div(::Type{I}) where {I<:Integer} = Rational{I}
+_promote_div(::Type{T}) where {T} = MA.promote_operation(/, T, Int)
+
+function MP.polynomial_type(
+    ::Type{Polynomial{B,M}},
+    ::Type{T},
+) where {B<:AbstractChebyshev,M,T}
+    return MP.polynomial_type(M, _promote_div(T))
 end
 
-reccurence_first_coef(::Type{<:AbstractChebyshevBasis}, degree) = 2
-reccurence_third_coef(::Type{<:AbstractChebyshevBasis}, degree) = -1
-reccurence_deno_coef(::Type{<:AbstractChebyshevBasis}, degree) = 1
+reccurence_first_coef(::Type{<:AbstractChebyshev}, degree) = 2
+reccurence_third_coef(::Type{<:AbstractChebyshev}, degree) = -1
+reccurence_deno_coef(::Type{<:AbstractChebyshev}, degree) = 1
 
 """
-    struct ChebyshevBasisFirstKind{P} <: AbstractChebyshevBasis{P}
-        polynomials::Vector{P}
-    end
+    struct ChebyshevFirstKind <: AbstractChebyshev end
 
 Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\frac{1}{\\sqrt{1 - x^2}}`` over the interval ``[-1, 1]``.
 """
-struct ChebyshevBasisFirstKind{P} <: AbstractChebyshevBasis{P}
-    polynomials::Vector{P}
-end
+struct ChebyshevFirstKind <: AbstractChebyshev end
+const Chebyshev = ChebyshevFirstKind
 
-const ChebyshevBasis{P} = ChebyshevBasisFirstKind{P}
+# https://en.wikipedia.org/wiki/Chebyshev_polynomials#Properties
+# T_n * T_m = T_{n + m} / 2 + T_{|n - m|} / 2
+function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial)
+    terms = [MP.term(1 // 1, MP.constant_monomial(a * b))]
+    vars_b = MP.variables(b)
+    var_state_b = iterate(vars_b)
+    for var_a in MP.variables(a)
+        if isnothing(var_state_b)
+            break
+        end
+        var_b, state_b = var_state_b
+        if var_b > var_a
+            var_state_b = iterate(vars_b, state_b)
+            for i in eachindex(terms)
+                terms[i] = MA.mul!!(terms[i], var_b^MP.degree(b, var_b))
+            end
+        elseif var_a > var_b
+            for i in eachindex(terms)
+                terms[i] = MA.mul!!(terms[i], var_a^MP.degree(a, var_a))
+            end
+        else
+            d_a = MP.degree(a, var_a)
+            d_b = MP.degree(b, var_b)
+            I = eachindex(terms)
+            for i in I
+                mono = MP.monomial(terms[i]) * var_a^(d_a + d_b)
+                terms[i] = MA.mul!!(terms[i], var_a^abs(d_a - d_b))
+                terms[i] = MA.operate!!(/, terms[i], 2)
+                α = MA.copy_if_mutable(MP.coefficient(terms[i]))
+                push!(terms, MP.term(α, mono))
+            end
+        end
+    end
+    return MP.polynomial!(terms)
+end
 
 function degree_one_univariate_polynomial(
-    ::Type{<:ChebyshevBasisFirstKind},
+    ::Type{Chebyshev},
     variable::MP.AbstractVariable,
 )
     MA.@rewrite(variable + 0)
 end
 
-function _scalar_product_function(::Type{<:ChebyshevBasisFirstKind}, i::Int)
+function _scalar_product_function(::Type{Chebyshev}, i::Int)
     if i == 0
         return π
     elseif isodd(i)
@@ -40,24 +77,20 @@ function _scalar_product_function(::Type{<:ChebyshevBasisFirstKind}, i::Int)
 end
 
 """
-    struct ChebyshevBasisSecondKind{P} <: AbstractChebyshevBasis{P}
-        polynomials::Vector{P}
-    end
+    struct ChebyshevSecondKind <: AbstractChebyshevBasis end
 
 Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\sqrt{1 - x^2}`` over the interval ``[-1, 1]``.
 """
-struct ChebyshevBasisSecondKind{P} <: AbstractChebyshevBasis{P}
-    polynomials::Vector{P}
-end
+struct ChebyshevSecondKind <: AbstractChebyshev end
 
 function degree_one_univariate_polynomial(
-    ::Type{<:ChebyshevBasisSecondKind},
+    ::Type{ChebyshevSecondKind},
     variable::MP.AbstractVariable,
 )
     MA.@rewrite(2variable + 0)
 end
 
-function _scalar_product_function(::Type{<:ChebyshevBasisSecondKind}, i::Int)
+function _scalar_product_function(::Type{<:ChebyshevSecondKind}, i::Int)
     if i == 0
         return π / 2
     elseif isodd(i)

src/fixed.jl

@@ -14,6 +14,7 @@ function empty_basis(
 ) where {PT}
     return B(PT[])
 end
+
 function MP.polynomial_type(
     ::AbstractPolynomialVectorBasis{PT},
     T::Type,
@@ -23,22 +24,6 @@ function MP.polynomial_type(
     V = MA.promote_operation(+, U, U)
     return MP.polynomial_type(PT, V)
 end
-function MP.polynomial(
-    f::Function,
-    basis::AbstractPolynomialVectorBasis{P},
-) where {P}
-    if isempty(generators(basis))
-        return zero(P)
-    else
-        return MP.polynomial(
-            mapreduce(
-                ip -> f(ip[1]) * ip[2],
-                MA.add!!,
-                enumerate(basis.polynomials),
-            ),
-        )
-    end
-end
 
 function _poly(::MA.Zero, ::Type{P}, ::Type{T}) where {P,T}
     return zero(MP.polynomial_type(P, T))

src/hermite.jl

@@ -1,13 +1,16 @@
-abstract type AbstractHermiteBasis{P} <: AbstractMultipleOrthogonalBasis{P} end
+abstract type AbstractHermite <: AbstractMultipleOrthogonal end
 
-function MP.polynomial_type(::Type{<:AbstractHermiteBasis}, V::Type)
-    return MP.polynomial_type(V, Int)
+function MP.polynomial_type(
+    ::Type{Polynomial{B,M}},
+    ::Type{T},
+) where {B<:AbstractHermite,M,T}
+    return MP.polynomial_type(M, float(T))
 end
 
-even_odd_separated(::Type{<:AbstractHermiteBasis}) = true
+even_odd_separated(::Type{<:AbstractHermite}) = true
 
-reccurence_second_coef(::Type{<:AbstractHermiteBasis}, degree) = 0
-reccurence_deno_coef(::Type{<:AbstractHermiteBasis}, degree) = 1
+reccurence_second_coef(::Type{<:AbstractHermite}, degree) = 0
+reccurence_deno_coef(::Type{<:AbstractHermite}, degree) = 1
 
 """
     struct ProbabilistsHermiteBasis{P} <: AbstractHermiteBasis{P}
@@ -16,21 +19,19 @@ reccurence_deno_coef(::Type{<:AbstractHermiteBasis}, degree) = 1
 
 Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\exp(-x^2/2)`` over the interval ``[-\\infty, \\infty]``.
 """
-struct ProbabilistsHermiteBasis{P} <: AbstractHermiteBasis{P}
-    polynomials::Vector{P}
-end
-reccurence_first_coef(::Type{<:ProbabilistsHermiteBasis}, degree) = 1
-function reccurence_third_coef(::Type{<:ProbabilistsHermiteBasis}, degree)
+struct ProbabilistsHermite <: AbstractHermite end
+reccurence_first_coef(::Type{ProbabilistsHermite}, degree) = 1
+function reccurence_third_coef(::Type{ProbabilistsHermite}, degree)
     return -(degree - 1)
 end
 function degree_one_univariate_polynomial(
-    ::Type{<:ProbabilistsHermiteBasis},
+    ::Type{ProbabilistsHermite},
     variable::MP.AbstractVariable,
 )
     MA.@rewrite(1variable)
 end
 
-function _scalar_product_function(::Type{<:ProbabilistsHermiteBasis}, i::Int)
+function _scalar_product_function(::Type{ProbabilistsHermite}, i::Int)
     if i == 0
         return √(2 * π)
     elseif isodd(i)
@@ -42,25 +43,23 @@ function _scalar_product_function(::Type{<:ProbabilistsHermiteBasis}, i::Int)
 end
 
 """
-    struct PhysicistsHermiteBasis{P} <: AbstractHermiteBasis{P}
+    struct PhysicistsHermite{P} <: AbstractHermite{P}
         polynomials::Vector{P}
     end
 
 Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\exp(-x^2)`` over the interval ``[-\\infty, \\infty]``.
 """
-struct PhysicistsHermiteBasis{P} <: AbstractHermiteBasis{P}
-    polynomials::Vector{P}
-end
-reccurence_first_coef(::Type{<:PhysicistsHermiteBasis}, degree) = 2
-reccurence_third_coef(::Type{<:PhysicistsHermiteBasis}, degree) = -2(degree - 1)
+struct PhysicistsHermite <: AbstractHermite end
+reccurence_first_coef(::Type{PhysicistsHermite}, degree) = 2
+reccurence_third_coef(::Type{PhysicistsHermite}, degree) = -2(degree - 1)
 function degree_one_univariate_polynomial(
-    ::Type{<:PhysicistsHermiteBasis},
+    ::Type{PhysicistsHermite},
     variable::MP.AbstractVariable,
 )
     MA.@rewrite(2variable)
 end
 
-function _scalar_product_function(::Type{<:PhysicistsHermiteBasis}, i::Int)
+function _scalar_product_function(::Type{PhysicistsHermite}, i::Int)
     if i == 0
         return √(π)
     elseif isodd(i)