From ed03cca8126a0238d07819f93cd1821d3f245b4b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 14 May 2024 16:26:53 +0200 Subject: [PATCH] up --- .github/workflows/ci.yml | 7 ++ src/MultivariateBases.jl | 5 +- src/chebyshev.jl | 134 ++++++--------------------------------- src/fixed.jl | 1 + src/hermite.jl | 36 +++++------ src/laguerre.jl | 26 ++++---- src/legendre.jl | 22 +++---- src/monomial.jl | 81 ++++++++++++----------- src/orthogonal.jl | 70 +++++--------------- src/orthonormal.jl | 2 +- src/polynomial.jl | 74 +++++++++++++++++++++ test/chebyshev.jl | 7 +- 12 files changed, 205 insertions(+), 260 deletions(-) create mode 100644 src/polynomial.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4b564d9..1b7c314 100644 --- a/.github/workflows/ci.yml +++ b/.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: diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index a8ce671..2515320 100644 --- a/src/MultivariateBases.jl +++ b/src/MultivariateBases.jl @@ -9,6 +9,7 @@ export maxdegree_basis, basis_covering_monomials, empty_basis include("interface.jl") export AbstractMonomialBasis, MonomialBasis, ScaledMonomialBasis +include("polynomial.jl") include("monomial.jl") include("scaled.jl") @@ -29,10 +30,10 @@ export generators, 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") diff --git a/src/chebyshev.jl b/src/chebyshev.jl index 809b589..d2e3c39 100644 --- a/src/chebyshev.jl +++ b/src/chebyshev.jl @@ -1,95 +1,26 @@ -# TODO Add to MultivariatePolynomials -Base.keytype(p::MP.AbstractPolynomial) = MP.monomial_type(p) -Base.valtype(p::MP.AbstractPolynomial) = MP.coefficient_type(p) -#Base.keys(p::MP.AbstractPolynomial) = MP.monomials(p) -Base.pairs(p::MP.AbstractPolynomial) = MP.terms(p) -function Base.similar(p::PT, ::Type{T}) where {PT<:MP.AbstractPolynomial,T} - return convert(MP.similar_type(PT, T), copy(p)) # Missing the `copy` in MP -end -Base.iterate(t::MP.Term) = iterate(t, 1) -function Base.iterate(t::MP.Term, state) - if state == 1 - return MP.monomial(t), 2 - elseif state == 2 - return MP.coefficient(t), 3 - else - return nothing - end -end -function MA.operate!( - ::SA.UnsafeAddMul{typeof(*)}, - mc::MP.AbstractPolynomial, - val, - c::MP.AbstractPolynomial, -) - return MA.operate!(MA.add_mul, mc, val, c) -end - -function SA.__canonicalize!(::MP.AbstractPolynomial) end - -struct Polynomial{B,M<:MP.AbstractMonomial} - monomial::M - function Polynomial{B}(mono::MP.AbstractMonomial) where {B} - return new{B,typeof(mono)}(mono) - end -end +export ChebyshevFirstKind, Chebyshev -function Polynomial{B}(v::MP.AbstractVariable) where {B} - return Polynomial{B}(MP.monomial(v)) -end - -function _algebra_element(p, ::Type{B}) where {B} - basis = MonomialIndexedBasis{B,MP.monomial_type(p)}() - return SA.AlgebraElement( - p, - SA.StarAlgebra( - Polynomial{B}(MP.constant_monomial(p)), - basis, - ), - ) -end - -function Base.:*(a::Polynomial{B}, b::Polynomial{B}) where {B} - _algebra_element(Mul{B}()(a.monomial, b.monomial), B) -end +abstract type AbstractChebyshev <: AbstractGegenbauer end -function _show(io::IO, mime::MIME, p::Polynomial{B}) where {B} - print(io, B) - print(io, "(") - show(io, mime, p.monomial) - print(io, ")") -end -function Base.show(io::IO, mime::MIME"text/plain", p::Polynomial) - _show(io, mime, p) -end -function Base.show(io::IO, p::Polynomial) - show(io, MIME"text/plain"(), p) -end - -struct MonomialIndexedBasis{B,M<:MP.AbstractMonomial} <: SA.ImplicitBasis{Polynomial{B,M},M} end -function Base.getindex(::MonomialIndexedBasis{B,M}, mono::M) where {B,M} - return Polynomial{B}(mono) -end -function Base.getindex(::MonomialIndexedBasis{B,M}, p::Polynomial{B,M}) where {B,M} - return mono -end -SA.mstructure(::MonomialIndexedBasis{B}) where {B} = Mul{B}() +_promote_div(::Type{I}) where {I<:Integer} = Rational{I} +_promote_div(::Type{T}) where {T} = MA.promote_operation(/, T, Int) -function Base.zero(::Type{Polynomial{B,M}}) where {B,M} - return _algebra_element(zero(MP.polynomial_type(M, Rational{Int})), B) +function MP.polynomial_type(::Type{<:AbstractChebyshev}, V::Type) + return MP.polynomial_type(V, _promote_div(MP.coefficient_type(V))) end -Base.zero(p::Polynomial) = zero(typeof(p)) - -struct Mul{B} <: SA.MultiplicativeStructure end +reccurence_first_coef(::Type{<:AbstractChebyshev}, degree) = 2 +reccurence_third_coef(::Type{<:AbstractChebyshev}, degree) = -1 +reccurence_deno_coef(::Type{<:AbstractChebyshev}, degree) = 1 -export ChebyshevFirstKind, ChebyshevFirstKindBasis +""" + 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 ChebyshevFirstKind end const Chebyshev = ChebyshevFirstKind -export Chebyshev - # 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) @@ -126,37 +57,14 @@ function (::Mul{Chebyshev})(a::MP.AbstractMonomial, b::MP.AbstractMonomial) return MP.polynomial!(terms) end -abstract type AbstractChebyshevBasis{P} <: AbstractGegenbauerBasis{P} end - -function MP.polynomial_type(::Type{<:AbstractChebyshevBasis}, V::Type) - return MP.polynomial_type(V, Float64) -end - -reccurence_first_coef(::Type{<:AbstractChebyshevBasis}, degree) = 2 -reccurence_third_coef(::Type{<:AbstractChebyshevBasis}, degree) = -1 -reccurence_deno_coef(::Type{<:AbstractChebyshevBasis}, degree) = 1 - -""" - struct ChebyshevBasisFirstKind{P} <: AbstractChebyshevBasis{P} - polynomials::Vector{P} - 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 - -const ChebyshevBasis{P} = ChebyshevBasisFirstKind{P} - 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) @@ -168,24 +76,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) diff --git a/src/fixed.jl b/src/fixed.jl index d7a735b..bd6ca20 100644 --- a/src/fixed.jl +++ b/src/fixed.jl @@ -14,6 +14,7 @@ function empty_basis( ) where {PT} return B(PT[]) end + function MP.polynomial_type( ::AbstractPolynomialVectorBasis{PT}, T::Type, diff --git a/src/hermite.jl b/src/hermite.jl index cf64476..76a5a73 100644 --- a/src/hermite.jl +++ b/src/hermite.jl @@ -1,13 +1,13 @@ -abstract type AbstractHermiteBasis{P} <: AbstractMultipleOrthogonalBasis{P} end +abstract type AbstractHermite{P} <: AbstractMultipleOrthogonal{P} end -function MP.polynomial_type(::Type{<:AbstractHermiteBasis}, V::Type) +function MP.polynomial_type(::Type{<:AbstractHermite}, V::Type) return MP.polynomial_type(V, Int) 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 +16,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 +40,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) diff --git a/src/laguerre.jl b/src/laguerre.jl index fcec105..cfe9670 100644 --- a/src/laguerre.jl +++ b/src/laguerre.jl @@ -1,32 +1,30 @@ """ - struct LaguerreBasis{P} <: AbstractMultipleOrthogonalBasis{P} - polynomials::Vector{P} - end + struct LaguerreBasis <: AbstractMultipleOrthogonal end Orthogonal polynomial with respect to the univariate weight function ``w(x) = \\exp(-x)`` over the interval ``[0, \\infty]``. """ -struct LaguerreBasis{P} <: AbstractMultipleOrthogonalBasis{P} - polynomials::Vector{P} -end +struct Laguerre <: AbstractMultipleOrthogonal end + +# TODO implement multiplication with https://www.jstor.org/stable/2002985 -function MP.polynomial_type(::Type{<:LaguerreBasis}, V::Type) +function MP.polynomial_type(::Type{Laguerre}, V::Type) return MP.polynomial_type(V, Float64) end -even_odd_separated(::Type{<:LaguerreBasis}) = false +even_odd_separated(::Type{Laguerre}) = false -reccurence_first_coef(::Type{<:LaguerreBasis}, degree) = -1 -reccurence_second_coef(::Type{<:LaguerreBasis}, degree) = (2degree - 1) -reccurence_third_coef(::Type{<:LaguerreBasis}, degree) = -(degree - 1) -reccurence_deno_coef(::Type{<:LaguerreBasis}, degree) = degree +reccurence_first_coef(::Type{Laguerre}, degree) = -1 +reccurence_second_coef(::Type{Laguerre}, degree) = (2degree - 1) +reccurence_third_coef(::Type{Laguerre}, degree) = -(degree - 1) +reccurence_deno_coef(::Type{Laguerre}, degree) = degree function degree_one_univariate_polynomial( - ::Type{<:LaguerreBasis}, + ::Type{Laguerre}, variable::MP.AbstractVariable, ) MA.@rewrite(1 - variable) end -function _scalar_product_function(::Type{<:LaguerreBasis}, i::Int) +function _scalar_product_function(::Type{Laguerre}, i::Int) return factorial(i) end diff --git a/src/legendre.jl b/src/legendre.jl index bb9b8e1..cd0e014 100644 --- a/src/legendre.jl +++ b/src/legendre.jl @@ -5,10 +5,10 @@ Orthogonal polynomial with respect to the univariate weight function ``w(x) = (1 - x^2)^{\\alpha - 1/2}`` over the interval ``[-1, 1]``. """ -abstract type AbstractGegenbauerBasis{P} <: AbstractMultipleOrthogonalBasis{P} end +abstract type AbstractGegenbauer <: AbstractMultipleOrthogonal end -even_odd_separated(::Type{<:AbstractGegenbauerBasis}) = true -reccurence_second_coef(::Type{<:AbstractGegenbauerBasis}, degree) = 0 +even_odd_separated(::Type{<:AbstractGegenbauer}) = true +reccurence_second_coef(::Type{<:AbstractGegenbauer}, degree) = 0 """ struct LegendreBasis{P} <: AbstractGegenbauerBasis{P} @@ -17,26 +17,24 @@ reccurence_second_coef(::Type{<:AbstractGegenbauerBasis}, degree) = 0 Orthogonal polynomial with respect to the univariate weight function ``w(x) = 1`` over the interval ``[-1, 1]``. """ -struct LegendreBasis{P} <: AbstractGegenbauerBasis{P} - polynomials::Vector{P} -end +struct Legendre <: AbstractGegenbauer end -function MP.polynomial_type(::Type{<:LegendreBasis}, V::Type) +function MP.polynomial_type(::Type{<:Legendre}, V::Type) return MP.polynomial_type(V, Float64) end -reccurence_first_coef(::Type{<:LegendreBasis}, degree) = (2degree - 1) -reccurence_third_coef(::Type{<:LegendreBasis}, degree) = -(degree - 1) -reccurence_deno_coef(::Type{<:LegendreBasis}, degree) = degree +reccurence_first_coef(::Type{Legendre}, degree) = (2degree - 1) +reccurence_third_coef(::Type{Legendre}, degree) = -(degree - 1) +reccurence_deno_coef(::Type{Legendre}, degree) = degree function degree_one_univariate_polynomial( - ::Type{<:LegendreBasis}, + ::Type{Legendre}, variable::MP.AbstractVariable, ) MA.@rewrite(variable + 0) end -function _scalar_product_function(::Type{<:LegendreBasis}, i::Int) +function _scalar_product_function(::Type{Legendre}, i::Int) if isodd(i) return 0 else diff --git a/src/monomial.jl b/src/monomial.jl index ee2dfc6..870ab9b 100644 --- a/src/monomial.jl +++ b/src/monomial.jl @@ -1,29 +1,52 @@ -abstract type AbstractMonomialBasis{ - MT<:MP.AbstractMonomial, - MV<:AbstractVector{MT}, -} <: AbstractPolynomialBasis end +struct ImplicitMonomialIndexedBasis{B<:AbstractBasis,M<:MP.AbstractMonomial} <: SA.ImplicitBasis{Polynomial{B,M},M} end -generators(basis::AbstractMonomialBasis) = basis.monomials +SA.mstructure(::ImplicitMonomialIndexedBasis{B}) where {B} = Mul{B}() -MP.monomial_type(::Type{<:AbstractMonomialBasis{MT}}) where {MT} = MT +const ExplicitMonomialIndexedBasis{B<:AbstractBasis,M<:MP.AbstractMonomial,V<:AbstractVector{M}} = SA.FixedBasis{Polynomial{B,M},Int,V,MTable{Polynomial{B,M},Int,V,Matrix{Polynomial{B,M}},Mul{B}}} -function empty_basis(MB::Type{<:AbstractMonomialBasis{MT}}) where {MT} - return MB(MP.empty_monomial_vector(MT)) +function MonomialBasis{B,M,V}( + monomials::V, +) where {B<:AbstractBasis,M<:MP.AbstractMonomial,V<:AbstractVector{MT}} + return new{B,M,V}(monomials) +end +function MonomialBasis{B}(monomials::AbstractVector) where {B<:AbstractBasis} + sorted = MP.monomial_vector(monomials) + return MonomialBasis{B,eltype(sorted),typeof(sorted)}(sorted) +end + +generators(basis::ExplicitMonomialIndexedBasis) = basis.monomials + +const MonomialIndexedBasis{B,M} = Union{ImplicitMonomialIndexedBasis{B,M}, ExplicitMonomialIndexedBasis{B,M}} + +function Base.getindex(::MonomialIndexedBasis{B,M}, mono::M) where {B,M} + return Polynomial{B}(mono) +end +function Base.getindex(::MonomialIndexedBasis{B,M}, p::Polynomial{B,M}) where {B,M} + return mono +end + +abstract type AbstractBasis end +abstract type AbstractMonomial <: AbstractBasis end + +MP.monomial_type(::Type{<:MonomialIndexedBasis{B,M}}) where {B,M} = M + +function empty_basis(::Type{ImplicitMonomialIndexedBasis{B,M}}) where {B,M} + return ExplicitMonomialIndexedBasis{B,M}(MP.empty_monomial_vector(M)) end function maxdegree_basis( - B::Type{<:AbstractMonomialBasis}, + ::Type{ImplicitMonomialIndexedBasis{B,M}}, variables, maxdegree::Int, -) - return B(MP.monomials(variables, 0:maxdegree)) +) where {B<:AbstractMonomial,M} + return ExplicitMonomialIndexedBasis{B,M}(MP.monomials(variables, 0:maxdegree)) end function basis_covering_monomials( - B::Type{<:AbstractMonomialBasis}, + ::Type{ImplicitMonomialIndexedBasis{B,M}}, monos::AbstractVector{<:MP.AbstractMonomial}, -) - return B(monos) +) where {B<:AbstractMonomial,M} + return ExplicitMonomialIndexedBasis{B,M}(monos) end # The `i`th index of output is the index of occurence of `x[i]` in `y`, @@ -42,7 +65,7 @@ function multi_findsorted(x, y) return I end -function merge_bases(basis1::MB, basis2::MB) where {MB<:AbstractMonomialBasis} +function merge_bases(basis1::MB, basis2::MB) where {MB<:ExplicitMonomialIndexedBasis} monos = MP.merge_monomial_vectors([basis1.monomials, basis2.monomials]) I1 = multi_findsorted(monos, basis1.monomials) I2 = multi_findsorted(monos, basis2.monomials) @@ -50,9 +73,7 @@ function merge_bases(basis1::MB, basis2::MB) where {MB<:AbstractMonomialBasis} end """ - struct MonomialBasis{MT<:MP.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis - monomials::MV - end + struct Monomial <: AbstractBasis end Monomial basis with the monomials of the vector `monomials`. For instance, `MonomialBasis([1, x, y, x^2, x*y, y^2])` is the monomial basis @@ -60,29 +81,15 @@ for the subspace of quadratic polynomials in the variables `x`, `y`. This basis is orthogonal under a scalar product defined with the complex Gaussian measure as density. Once normalized so as to be orthonormal with this scalar product, -one get ths [`ScaledMonomialBasis`](@ref). +one get ths [`ScaledMonomial`](@ref). """ -struct MonomialBasis{MT<:MP.AbstractMonomial,MV<:AbstractVector{MT}} <: - AbstractMonomialBasis{MT,MV} - monomials::MV - function MonomialBasis{MT,MV}( - monomials::MV, - ) where {MT<:MP.AbstractMonomial,MV<:AbstractVector{MT}} - return new{MT,MV}(monomials) - end - function MonomialBasis(monomials::AbstractVector) - sorted = MP.monomial_vector(monomials) - return MonomialBasis{eltype(sorted),typeof(sorted)}(sorted) - end -end - -Base.getindex(basis::MonomialBasis, i::Int) = basis.monomials[i] +struct Monomial <: AbstractMonomial end function MP.polynomial_type( - ::Union{MonomialBasis{MT},Type{<:MonomialBasis{MT}}}, + ::Union{MonomialIndexedBasis{B,M},Type{<:MonomialIndexedBasis{B,M}}}, T::Type, -) where {MT} - return MP.polynomial_type(MT, T) +) where {B,M} + return MP.polynomial_type(M, T) end MP.polynomial(f::Function, mb::MonomialBasis) = MP.polynomial(f, mb.monomials) diff --git a/src/orthogonal.jl b/src/orthogonal.jl index 5085c4a..e2b525b 100644 --- a/src/orthogonal.jl +++ b/src/orthogonal.jl @@ -20,8 +20,7 @@ where [`reccurence_first_coef`](@ref) gives `a_k`, [`reccurence_third_coef`](@ref) gives `c_k` and [`reccurence_deno_coef`](@ref) gives `d_k`. """ -abstract type AbstractMultipleOrthogonalBasis{P} <: - AbstractPolynomialVectorBasis{P,Vector{P}} end +abstract type AbstractMultipleOrthogonal <: AbstractBasis end """ reccurence_first_coef(B::Type{<:AbstractMultipleOrthogonalBasis}, degree::Integer) @@ -71,7 +70,7 @@ Return the vector of univariate polynomials of the basis `B` up to `degree` with variable `variable`. """ function univariate_orthogonal_basis( - B::Type{<:AbstractMultipleOrthogonalBasis}, + B::Type{<:AbstractMultipleOrthogonal}, variable::MP.AbstractVariable, degree::Integer, ) @@ -100,44 +99,12 @@ function univariate_orthogonal_basis( end end -function _basis_from_monomials( - B::Type{<:AbstractMultipleOrthogonalBasis}, - variables, - monos, -) - univariate = [ - univariate_orthogonal_basis( - B, - variable, - maximum(mono -> MP.degree(mono, variable), monos), - ) for variable in variables - ] - return B([ - prod( - i -> univariate[i][MP.degree(mono, variables[i])+1], - eachindex(variables), - ) for mono in monos - ]) -end - -function maxdegree_basis( - B::Type{<:AbstractMultipleOrthogonalBasis}, - variables, - maxdegree::Int, -) - return _basis_from_monomials( - B, - variables, - MP.monomials(variables, 0:maxdegree), - ) -end - function basis_covering_monomials( - B::Type{<:AbstractMultipleOrthogonalBasis}, - monos::AbstractVector{<:MP.AbstractMonomial}, -) + ::Type{ImplicitMonomialIndexedBasis{B,M1}}, + monos::AbstractVector{M2}, +) where {B<:AbstractMultipleOrthogonal,M1,M2<:MP.AbstractMonomial} to_add = collect(monos) - m = Set(monos) + m = Set{promote_type(M1, M2)}(monos) while !isempty(to_add) mono = pop!(to_add) for v in MP.variables(mono) @@ -152,43 +119,39 @@ function basis_covering_monomials( end end end - return _basis_from_monomials( - B, - MP.variables(monos), - MP.monomial_vector(collect(m)), - ) + return ExplicitMonomialIndexedBasis{B}(MP.monomial_vector(collect(m))) end function _scalar_product_function( - ::Type{<:AbstractMultipleOrthogonalBasis}, + ::Type{<:AbstractMultipleOrthogonal}, i::Int, ) end function LinearAlgebra.dot( p, q, - basis_type::Type{<:AbstractMultipleOrthogonalBasis}, + basis_type::Type{<:AbstractMultipleOrthogonal}, ) return _integral(p * q, basis_type) end function _integral( p::Number, - basis_type::Type{<:AbstractMultipleOrthogonalBasis}, + basis_type::Type{<:AbstractMultipleOrthogonal}, ) return p * _scalar_product_function(basis_type, 0) end function _integral( - p::MP.AbstractVariable, - basis_type::Type{<:AbstractMultipleOrthogonalBasis}, + ::MP.AbstractVariable, + basis_type::Type{<:AbstractMultipleOrthogonal}, ) return _scalar_product_function(basis_type, 1) end function _integral( p::MP.AbstractMonomial, - basis_type::Type{<:AbstractMultipleOrthogonalBasis}, + basis_type::Type{<:AbstractMultipleOrthogonal}, ) return prod([ _scalar_product_function(basis_type, i) for i in MP.exponents(p) @@ -197,20 +160,19 @@ end function _integral( p::MP.AbstractTerm, - basis_type::Type{<:AbstractMultipleOrthogonalBasis}, + basis_type::Type{<:AbstractMultipleOrthogonal}, ) return MP.coefficient(p) * _integral(MP.monomial(p), basis_type) end function _integral( p::MP.AbstractPolynomial, - basis_type::Type{<:AbstractMultipleOrthogonalBasis}, + basis_type::Type{<:AbstractMultipleOrthogonal}, ) return sum([_integral(t, basis_type) for t in MP.terms(p)]) end -function MP.coefficients(p, basis::AbstractMultipleOrthogonalBasis) - B = typeof(basis) +function MP.coefficients(p, basis::ExplicitMonomialIndexedBasis{B,M}) where {B<:AbstractMultipleOrthogonal,M} return [ LinearAlgebra.dot(p, el, B) / LinearAlgebra.dot(el, el, B) for el in basis diff --git a/src/orthonormal.jl b/src/orthonormal.jl index 5475a7d..0827f76 100644 --- a/src/orthonormal.jl +++ b/src/orthonormal.jl @@ -12,7 +12,7 @@ polynomial basis for cubic polynomials in the variable `x`. struct OrthonormalCoefficientsBasis{ PT<:MP.AbstractPolynomialLike, PV<:AbstractVector{PT}, -} <: AbstractPolynomialVectorBasis{PT,PV} +} <: AbstractBasis{PT,PV} polynomials::PV end diff --git a/src/polynomial.jl b/src/polynomial.jl new file mode 100644 index 0000000..a4545c0 --- /dev/null +++ b/src/polynomial.jl @@ -0,0 +1,74 @@ +# TODO Add to MultivariatePolynomials +Base.keytype(p::MP.AbstractPolynomial) = MP.monomial_type(p) +Base.valtype(p::MP.AbstractPolynomial) = MP.coefficient_type(p) +#Base.keys(p::MP.AbstractPolynomial) = MP.monomials(p) +SA.nonzero_pairs(p::MP.AbstractPolynomial) = MP.terms(p) +function Base.similar(p::PT, ::Type{T}) where {PT<:MP.AbstractPolynomial,T} + return convert(MP.similar_type(PT, T), copy(p)) # Missing the `copy` in MP +end +Base.iterate(t::MP.Term) = iterate(t, 1) +function Base.iterate(t::MP.Term, state) + if state == 1 + return MP.monomial(t), 2 + elseif state == 2 + return MP.coefficient(t), 3 + else + return nothing + end +end +function MA.operate!( + ::SA.UnsafeAddMul{typeof(*)}, + mc::MP.AbstractPolynomial, + val, + c::MP.AbstractPolynomial, +) + return MA.operate!(MA.add_mul, mc, val, c) +end +function MA.operate!(::typeof(SA.canonical), ::MP.AbstractPolynomial) end + +struct Polynomial{B,M<:MP.AbstractMonomial} + monomial::M + function Polynomial{B}(mono::MP.AbstractMonomial) where {B} + return new{B,typeof(mono)}(mono) + end +end + +function Polynomial{B}(v::MP.AbstractVariable) where {B} + return Polynomial{B}(MP.monomial(v)) +end + +function _algebra_element(p, ::Type{B}) where {B} + basis = MonomialIndexedBasis{B,MP.monomial_type(p)}() + return SA.AlgebraElement( + p, + SA.StarAlgebra( + Polynomial{B}(MP.constant_monomial(p)), + basis, + ), + ) +end + +function Base.:*(a::Polynomial{B}, b::Polynomial{B}) where {B} + _algebra_element(Mul{B}()(a.monomial, b.monomial), B) +end + +function _show(io::IO, mime::MIME, p::Polynomial{B}) where {B} + print(io, B) + print(io, "(") + show(io, mime, p.monomial) + print(io, ")") +end +function Base.show(io::IO, mime::MIME"text/plain", p::Polynomial) + _show(io, mime, p) +end +function Base.show(io::IO, p::Polynomial) + show(io, MIME"text/plain"(), p) +end + +function Base.zero(::Type{Polynomial{B,M}}) where {B,M} + return _algebra_element(zero(MP.polynomial_type(M, Rational{Int})), B) +end + +Base.zero(p::Polynomial) = zero(typeof(p)) + +struct Mul{B} <: SA.MultiplicativeStructure end diff --git a/test/chebyshev.jl b/test/chebyshev.jl index cbb3f8c..fda6013 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -6,12 +6,9 @@ using DynamicPolynomials @polyvar x a = MultivariateBases.Polynomial{Chebyshev}(x) b = a * a - @show b.coeffs - #@test length(b.coeffs.basis_elements) == 2 + @test b.coeffs == 1//2 + 1//2 * x^2 c = b * b - @show c.coeffs - # TODO after __canonicalize, it should be smaller - #@test length(c.coeffs.basis_elements) == 8 + @test c.coeffs == 3//8 + 1//2 * x^2 + 1//8 * x^4 end @testset "Orthogonal" begin