Merge pull request #2 from JuliaAlgebra/bl/orthogonal

Add classical multiple orthogonal basis

Author: blegat

.travis.yml

@@ -18,6 +18,8 @@ jobs:
       julia: 1.0
       os: linux
       script:
-        - julia --project=docs/ -e 'import Pkg; Pkg.instantiate(); Pkg.develop(Pkg.PackageSpec(path=pwd()))'
+        - julia --project=docs/ -e 'import Pkg; Pkg.instantiate();
+                                    Pkg.add("DynamicPolynomials");
+                                    Pkg.develop(Pkg.PackageSpec(path=pwd()))'
         - julia --project=docs/ docs/make.jl
       after_success: skip

docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
 
 [compat]
 Documenter = "~0.21"

docs/src/index.md

@@ -1,11 +1,42 @@
+```@meta
+DocTestSetup = quote
+    using MultivariateBases
+end
+```
+
 # MultivariateBases
 
 [MultivariateBases.jl](https://github.com/JuliaAlgebra/MultivariateBases.jl) is a standardized API for multivariate polynomial bases
 based on the [MultivariatePolynomials](https://github.com/JuliaAlgebra/MultivariatePolynomials.jl) API.
 
 ```@docs
 AbstractPolynomialBasis
+maxdegree_basis
+basis_covering_monomials
 FixedPolynomialBasis
+```
+
+## Monomial basis
+
+```@docs
 MonomialBasis
 ScaledMonomialBasis
 ```
+
+## Orthogonal basis
+
+```@docs
+AbstractMultipleOrthogonalBasis
+univariate_orthogonal_basis
+reccurence_first_coef
+reccurence_second_coef
+reccurence_third_coef
+reccurence_deno_coef
+ProbabilistsHermiteBasis
+PhysicistsHermiteBasis
+LaguerreBasis
+AbstractGegenbauerBasis
+LegendreBasis
+ChebyshevBasisFirstKind
+ChebyshevBasisSecondKind
+```

src/MultivariateBases.jl

@@ -7,15 +7,21 @@ using MultivariatePolynomials
 const MP = MultivariatePolynomials
 
 export AbstractPolynomialBasis
-export maxdegree_basis, empty_basis
+export maxdegree_basis, basis_covering_monomials, empty_basis
 include("interface.jl")
 
 export AbstractMonomialBasis, MonomialBasis, ScaledMonomialBasis
 include("monomial.jl")
 include("scaled.jl")
 
-export FixedPolynomialBasis, ChebyshevBasis
+export FixedPolynomialBasis, AbstractMultipleOrthogonalBasis, ProbabilistsHermiteBasis, PhysicistsHermiteBasis, LaguerreBasis
+export AbstractGegenbauerBasis, LegendreBasis, ChebyshevBasis, ChebyshevBasisFirstKind, ChebyshevBasisSecondKind
+export univariate_orthogonal_basis, reccurence_first_coef, reccurence_second_coef, reccurence_third_coef, reccurence_deno_coef
 include("fixed.jl")
+include("orthogonal.jl")
+include("hermite.jl")
+include("laguerre.jl")
+include("legendre.jl")
 include("chebyshev.jl")
 
 end # module

src/chebyshev.jl

@@ -1,26 +1,35 @@
-struct ChebyshevBasis{P} <: AbstractPolynomialVectorBasis{P, Vector{P}}
+abstract type AbstractChebyshevBasis{P} <: AbstractGegenbauerBasis{P} end
+
+polynomial_type(::Type{<:AbstractChebyshevBasis}, V::Type) = MP.polynomialtype(V, Float64)
+
+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
 
-function chebyshev_polynomial_first_kind(variable::MP.AbstractVariable, degree::Integer)
-    @assert degree >= 0
-    if degree == 0
-        return [MA.@rewrite(0 * variable + 1)]
-    elseif degree == 1
-        return push!(chebyshev_polynomial_first_kind(variable, 0),
-                     MA.@rewrite(1 * variable + 0))
-    else
-        previous = chebyshev_polynomial_first_kind(variable, degree - 1)
-        next = MA.@rewrite(2variable * previous[degree] - previous[degree - 1])
-        push!(previous, next)
-        return previous
+const ChebyshevBasis{P} = ChebyshevBasisFirstKind{P}
+
+degree_one_univariate_polynomial(::Type{<:ChebyshevBasisFirstKind}, variable::MP.AbstractVariable) = MA.@rewrite(variable + 0)
+
+"""
+    struct ChebyshevBasisSecondKind{P} <: AbstractChebyshevBasis{P}
+        polynomials::Vector{P}
     end
-end
 
-function maxdegree_basis(B::Type{ChebyshevBasis}, variables, maxdegree::Int)
-    univariate = [chebyshev_polynomial_first_kind(variable, maxdegree) for variable in variables]
-    return ChebyshevBasis([
-        prod(i -> univariate[i][degree(mono, variables[i]) + 1],
-             eachindex(variables))
-        for mono in MP.monomials(variables, 0:maxdegree)])
+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
+
+degree_one_univariate_polynomial(::Type{<:ChebyshevBasisSecondKind}, variable::MP.AbstractVariable) = MA.@rewrite(2variable + 0)

src/hermite.jl

@@ -0,0 +1,36 @@
+abstract type AbstractHermiteBasis{P} <: AbstractMultipleOrthogonalBasis{P} end
+
+polynomial_type(::Type{<:AbstractHermiteBasis}, V::Type) = MP.polynomialtype(V, Int)
+
+even_odd_separated(::Type{<:AbstractHermiteBasis}) = true
+
+reccurence_second_coef(::Type{<:AbstractHermiteBasis}, degree) = 0
+reccurence_deno_coef(::Type{<:AbstractHermiteBasis}, degree) = 1
+
+"""
+    struct ProbabilistsHermiteBasis{P} <: AbstractHermiteBasis{P}
+        polynomials::Vector{P}
+    end
+
+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
+reccurence_third_coef(::Type{<:ProbabilistsHermiteBasis}, degree) = -(degree - 1)
+degree_one_univariate_polynomial(::Type{<:ProbabilistsHermiteBasis}, variable::MP.AbstractVariable) = MA.@rewrite(1variable)
+
+"""
+    struct PhysicistsHermiteBasis{P} <: AbstractHermiteBasis{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)
+degree_one_univariate_polynomial(::Type{<:PhysicistsHermiteBasis}, variable::MP.AbstractVariable) = MA.@rewrite(2variable)

src/interface.jl

@@ -12,3 +12,34 @@ abstract type AbstractPolynomialBasis end
 function MP.polynomial(coefs::Vector, basis::AbstractPolynomialBasis)
     return MP.polynomial(i -> coefs[i], basis)
 end
+
+"""
+    maxdegree_basis(B::Type{<:AbstractPolynomialBasis}, variables, maxdegree::Int)
+
+Return the basis of type `B` generating all polynomials of degree up to
+`maxdegree` with variables `variables`.
+"""
+function maxdegree_basis end
+
+"""
+    basis_covering_monomials(B::Type{<:AbstractPolynomialBasis}, monos::AbstractVector{<:AbstractMonomial})
+
+Return the minimal basis of type `B` that can generate all polynomials of the
+monomial basis generated by `monos`.
+
+## Examples
+
+For example, to generate all the polynomials with nonzero coefficients for the
+monomials `x^4` and `x^2`, we need three polynomials as otherwise, we generate
+polynomials with nonzero constant term.
+```jldoctest
+julia> using DynamicPolynomials
+
+julia> @polyvar x
+(x,)
+
+julia> basis_covering_monomials(ChebyshevBasis, [x^4, x^2])
+ChebyshevBasisFirstKind{Polynomial{true,Float64}}(DynamicPolynomials.Polynomial{true,Float64}[8.0x⁴ - 8.0x² + 1.0, 2.0x² - 1.0, 1.0])
+```
+"""
+function basis_covering_monomials end

src/laguerre.jl

@@ -0,0 +1,21 @@
+"""
+    struct LaguerreBasis{P} <: AbstractMultipleOrthogonalBasis{P}
+        polynomials::Vector{P}
+    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
+
+polynomial_type(::Type{<:LaguerreBasis}, V::Type) = MP.polynomialtype(V, Float64)
+
+even_odd_separated(::Type{<:LaguerreBasis}) = 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
+
+degree_one_univariate_polynomial(::Type{<:LaguerreBasis}, variable::MP.AbstractVariable) = MA.@rewrite(1 - variable)

src/legendre.jl

@@ -0,0 +1,30 @@
+"""
+    struct AbstractGegenbauerBasis{P} <: AbstractMultipleOrthogonalBasis{P}
+        polynomials::Vector{P}
+    end
+
+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
+
+even_odd_separated(::Type{<:AbstractGegenbauerBasis}) = true
+reccurence_second_coef(::Type{<:AbstractGegenbauerBasis}, degree) = 0
+
+"""
+    struct LegendreBasis{P} <: AbstractGegenbauerBasis{P}
+        polynomials::Vector{P}
+    end
+
+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
+
+polynomial_type(::Type{<:LegendreBasis}, V::Type) = MP.polynomialtype(V, Float64)
+
+reccurence_first_coef(::Type{<:LegendreBasis}, degree) = (2degree - 1)
+reccurence_third_coef(::Type{<:LegendreBasis}, degree) = -(degree - 1)
+reccurence_deno_coef(::Type{<:LegendreBasis}, degree) = degree
+
+degree_one_univariate_polynomial(::Type{<:LegendreBasis}, variable::MP.AbstractVariable) = MA.@rewrite(variable + 0)

src/monomial.jl

@@ -11,6 +11,9 @@ empty_basis(MB::Type{<:AbstractMonomialBasis{MT}}) where {MT} = MB(MP.emptymonov
 function maxdegree_basis(B::Type{<:AbstractMonomialBasis}, variables, maxdegree::Int)
     return B(MP.monomials(variables, 0:maxdegree))
 end
+function basis_covering_monomials(B::Type{<:AbstractMonomialBasis}, monos::AbstractVector{<:AbstractMonomial})
+    return B(monos)
+end
 
 # The `i`th index of output is the index of occurence of `x[i]` in `y`,
 # or `0` if it does not occur.