From 9d100071e8336acf4dd402a4ff3c0fae2bf0008f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 30 Jan 2020 13:54:46 +0100 Subject: [PATCH] Add classical multiple orthogonal basis --- .travis.yml | 4 +- docs/Project.toml | 1 + docs/src/index.md | 31 ++++++++++ src/MultivariateBases.jl | 10 ++- src/chebyshev.jl | 49 +++++++++------ src/hermite.jl | 36 +++++++++++ src/interface.jl | 31 ++++++++++ src/laguerre.jl | 21 +++++++ src/legendre.jl | 30 +++++++++ src/monomial.jl | 3 + src/orthogonal.jl | 128 +++++++++++++++++++++++++++++++++++++++ test/chebyshev.jl | 24 +++++--- test/hermite.jl | 24 ++++++++ test/laguerre.jl | 16 +++++ test/legendre.jl | 16 +++++ test/runtests.jl | 74 ++++++++++++++++++---- 16 files changed, 454 insertions(+), 44 deletions(-) create mode 100644 src/hermite.jl create mode 100644 src/laguerre.jl create mode 100644 src/legendre.jl create mode 100644 src/orthogonal.jl create mode 100644 test/hermite.jl create mode 100644 test/laguerre.jl create mode 100644 test/legendre.jl diff --git a/.travis.yml b/.travis.yml index 2970c74..e6d54e2 100644 --- a/.travis.yml +++ b/.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 diff --git a/docs/Project.toml b/docs/Project.toml index ce87d15..47be4ee 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,6 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" [compat] Documenter = "~0.21" diff --git a/docs/src/index.md b/docs/src/index.md index cb20259..ea41871 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,3 +1,9 @@ +```@meta +DocTestSetup = quote + using MultivariateBases +end +``` + # MultivariateBases [MultivariateBases.jl](https://github.com/JuliaAlgebra/MultivariateBases.jl) is a standardized API for multivariate polynomial bases @@ -5,7 +11,32 @@ based on the [MultivariatePolynomials](https://github.com/JuliaAlgebra/Multivari ```@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 +``` diff --git a/src/MultivariateBases.jl b/src/MultivariateBases.jl index 40383ea..b7342d7 100644 --- a/src/MultivariateBases.jl +++ b/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 diff --git a/src/chebyshev.jl b/src/chebyshev.jl index ee3e717..f7c37e2 100644 --- a/src/chebyshev.jl +++ b/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) diff --git a/src/hermite.jl b/src/hermite.jl new file mode 100644 index 0000000..ee26091 --- /dev/null +++ b/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) diff --git a/src/interface.jl b/src/interface.jl index 2c4ac8b..6ff188e 100644 --- a/src/interface.jl +++ b/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 diff --git a/src/laguerre.jl b/src/laguerre.jl new file mode 100644 index 0000000..6441375 --- /dev/null +++ b/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) diff --git a/src/legendre.jl b/src/legendre.jl new file mode 100644 index 0000000..27d1cc2 --- /dev/null +++ b/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) diff --git a/src/monomial.jl b/src/monomial.jl index 24de6c0..a26df65 100644 --- a/src/monomial.jl +++ b/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. diff --git a/src/orthogonal.jl b/src/orthogonal.jl new file mode 100644 index 0000000..1fe94c5 --- /dev/null +++ b/src/orthogonal.jl @@ -0,0 +1,128 @@ +""" + abstract type AbstractMultipleOrthogonalBasis{P} <: AbstractPolynomialVectorBasis{P, Vector{P}} end + +Polynomial basis such that ``\\langle p_i(x), p_j(x) \\rangle = 0`` if ``i \\neq j`` where +```math +\\langle p(x), q(x) \\rangle = \\int p(x)q(x) w(x) dx +``` +where the weight is a product of weight functions +``w(x) = w_1(x_1)w_2(x_2) \\cdots w_n(x_n)`` in each variable. +The polynomial of the basis are product of univariate polynomials: +``p(x) = p_1(x_1)p_2(x_2) \\cdots p_n(x_n)``. +where the univariate polynomials of variable `x_i` form an univariate +orthogonal basis for the weight function `w_i(x_i)`. +Therefore, they satisfy the recurrence relation +```math +d_k p_k(x_i) = (a_k x_i + b_k) p_{k-1}(x_i) + c_k p_{k-2}(x_i) +``` +where [`reccurence_first_coef`](@ref) gives `a_k`, +[`reccurence_second_coef`](@ref) gives `b_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 + +""" + reccurence_first_coef(B::Type{<:AbstractMultipleOrthogonalBasis}, degree::Integer) + +Return `a_{degree}` in recurrence equation +```math +d_k p_k(x_i) = (a_k x_i + b_k) p_{k-1}(x_i) + c_k p_{k-2}(x_i) +``` +""" +function reccurence_first_coef end + +""" + reccurence_second_coef(B::Type{<:AbstractMultipleOrthogonalBasis}, degree::Integer) + +Return `b_{degree}` in recurrence equation +```math +d_k p_k(x_i) = (a_k x_i + b_k) p_{k-1}(x_i) + c_k p_{k-2}(x_i) +``` +""" +function reccurence_second_coef end + +""" + reccurence_third_coef(B::Type{<:AbstractMultipleOrthogonalBasis}, degree::Integer) + +Return `c_{degree}` in recurrence equation +```math +d_k p_k(x_i) = (a_k x_i + b_k) p_{k-1}(x_i) + c_k p_{k-2}(x_i) +``` +""" +function reccurence_third_coef end + +""" + reccurence_deno_coef(B::Type{<:AbstractMultipleOrthogonalBasis}, degree::Integer) + +Return `d_{degree}` in recurrence equation +```math +d_k p_k(x_i) = (a_k x_i + b_k) p_{k-1}(x_i) + c_k p_{k-2}(x_i) +``` +""" +function reccurence_deno_coef end + +""" + univariate_orthogonal_basis(B::Type{<:AbstractMultipleOrthogonalBasis}, + variable::MP.AbstractVariable, degree::Integer) + +Return the vector of univariate polynomials of the basis `B` up to `degree` +with variable `variable`. +""" +function univariate_orthogonal_basis( + B::Type{<:AbstractMultipleOrthogonalBasis}, variable::MP.AbstractVariable, degree::Integer) + + @assert degree >= 0 + if degree == 0 + return polynomial_type(B, typeof(variable))[one(variable)] + elseif degree == 1 + return push!(univariate_orthogonal_basis(B, variable, 0), + degree_one_univariate_polynomial(B, variable)) + else + previous = univariate_orthogonal_basis(B, variable, degree - 1) + a = reccurence_first_coef(B, degree) + b = reccurence_second_coef(B, degree) + c = reccurence_third_coef(B, degree) + d = reccurence_deno_coef(B, degree) + next = MA.@rewrite((a * variable + b) * previous[degree] + c * previous[degree - 1]) + if !isone(d) + next = next / d + end + push!(previous, next) + return previous + end +end + +function _basis_from_monomials(B::Type{<:AbstractMultipleOrthogonalBasis}, variables, monos) + univariate = [univariate_orthogonal_basis( + B, variable, maximum(mono -> degree(mono, variable), monos)) + for variable in variables] + return B([ + prod(i -> univariate[i][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{<:AbstractMonomial}) + to_add = collect(monos) + m = Set(monos) + while !isempty(to_add) + mono = pop!(to_add) + for v in MP.variables(mono) + step = even_odd_separated(B) ? 2 : 1 + vstep = v^step + if MP.divides(vstep, mono) + new_mono = MP.mapexponents(-, mono, vstep) + if !(new_mono in m) + push!(m, new_mono) + push!(to_add, new_mono) + end + end + end + end + return _basis_from_monomials(B, variables(monos), MP.monovec(collect(m))) +end diff --git a/test/chebyshev.jl b/test/chebyshev.jl index 509bb7d..d56ad7e 100644 --- a/test/chebyshev.jl +++ b/test/chebyshev.jl @@ -1,16 +1,24 @@ using Test using MultivariateBases -using DynamicPolynomials -@polyvar x y -@testset "Univariate" begin - basis = maxdegree_basis(ChebyshevBasis, (x,), 3) - @test basis.polynomials[4] == 1 - @test basis.polynomials[3] == x - @test basis.polynomials[2] == 2x^2 - 1 - @test basis.polynomials[1] == 4x^3 - 3x +@testset "Orthogonal" begin + orthogonal_test(ChebyshevBasis, x -> [ + 1, + x, + 2x^2 - 1, + 4x^3 - 3x, + 8x^4 - 8x^2 + 1 + ], true) + orthogonal_test(ChebyshevBasisSecondKind, x -> [ + 1, + 2x, + 4x^2 - 1, + 8x^3 - 4x, + 16x^4 - 12x^2 + 1 + ], true) end @testset "API degree = $degree" for degree in 0:3 api_test(ChebyshevBasis, degree) + api_test(ChebyshevBasisSecondKind, degree) end diff --git a/test/hermite.jl b/test/hermite.jl new file mode 100644 index 0000000..8b02a36 --- /dev/null +++ b/test/hermite.jl @@ -0,0 +1,24 @@ +using Test +using MultivariateBases + +@testset "Orthogonal" begin + orthogonal_test(ProbabilistsHermiteBasis, x -> [ + 1, + x, + x^2 - 1, + x^3 - 3x, + x^4 - 6x^2 + 3 + ], true) + orthogonal_test(PhysicistsHermiteBasis, x -> [ + 1, + 2x, + 4x^2 - 2, + 8x^3 - 12x, + 16x^4 - 48x^2 + 12 + ], true) +end + +@testset "API degree = $degree" for degree in 0:3 + api_test(ProbabilistsHermiteBasis, degree) + api_test(PhysicistsHermiteBasis, degree) +end diff --git a/test/laguerre.jl b/test/laguerre.jl new file mode 100644 index 0000000..2c8bf5e --- /dev/null +++ b/test/laguerre.jl @@ -0,0 +1,16 @@ +using Test +using MultivariateBases + +@testset "Orthogonal" begin + orthogonal_test(LaguerreBasis, x -> [ + 1, + 1 - x, + (x^2 - 4x + 2) / 2, + (-x^3 + 9x^2 - 18x + 6) / 6, + (x^4 - 16x^3 + 72x^2 - 96x + 24) / 24 + ], false) +end + +@testset "API degree = $degree" for degree in 0:3 + api_test(LaguerreBasis, degree) +end diff --git a/test/legendre.jl b/test/legendre.jl new file mode 100644 index 0000000..05b51eb --- /dev/null +++ b/test/legendre.jl @@ -0,0 +1,16 @@ +using Test +using MultivariateBases + +@testset "Orthogonal" begin + orthogonal_test(LegendreBasis, x -> [ + 1, + x, + (3x^2 - 1) / 2, + (5x^3 - 3x) / 2, + (35x^4 - 30x^2 + 3) / 8 + ], true) +end + +@testset "API degree = $degree" for degree in 0:3 + api_test(LegendreBasis, degree) +end diff --git a/test/runtests.jl b/test/runtests.jl index 32fae5c..4ca135c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,19 +6,58 @@ using DynamicPolynomials function api_test(B::Type{<:AbstractPolynomialBasis}, degree) @polyvar x[1:2] - basis = maxdegree_basis(B, x, degree) - n = binomial(2 + degree, 2) - @test length(basis) == n - @test typeof(copy(basis)) == typeof(basis) - @test nvariables(basis) == 2 - @test variables(basis) == x - @test monomialtype(typeof(basis)) == monomialtype(x[1]) - @test typeof(empty_basis(typeof(basis))) == typeof(basis) - @test length(empty_basis(typeof(basis))) == 0 - @test polynomialtype(basis, Float64) == polynomialtype(x[1], Float64) - @test polynomial(i -> 0.0, basis) isa polynomialtype(basis, Float64) - @test polynomial(zeros(n, n), basis, Float64) isa polynomialtype(basis, Float64) - @test polynomial(ones(n, n), basis, Float64) isa polynomialtype(basis, Float64) + for basis in [ + maxdegree_basis(B, x, degree), + basis_covering_monomials(B, monomials(x, 0:degree)) + ] + n = binomial(2 + degree, 2) + @test length(basis) == n + @test typeof(copy(basis)) == typeof(basis) + @test nvariables(basis) == 2 + @test variables(basis) == x + @test monomialtype(typeof(basis)) == monomialtype(x[1]) + @test typeof(empty_basis(typeof(basis))) == typeof(basis) + @test length(empty_basis(typeof(basis))) == 0 + @test polynomialtype(basis, Float64) == polynomialtype(x[1], Float64) + @test polynomial(i -> 0.0, basis) isa polynomialtype(basis, Float64) + @test polynomial(zeros(n, n), basis, Float64) isa polynomialtype(basis, Float64) + @test polynomial(ones(n, n), basis, Float64) isa polynomialtype(basis, Float64) + end +end + +function orthogonal_test(B::Type{<:AbstractMultipleOrthogonalBasis}, univ::Function, even_odd_separated::Bool) + @test MultivariateBases.even_odd_separated(B) == even_odd_separated + @polyvar x y + univariate_x = univ(x) + univariate_y = univ(y) + + @testset "Univariate $var" for (var, univ) in [(x, univariate_x), (y, univariate_y)] + basis = maxdegree_basis(B, (var,), 4) + for i in 1:5 + @test basis.polynomials[length(basis) + 1 - i] == univ[i] + end + end + + @testset "basis_covering_monomials" begin + monos = basis_covering_monomials(B, monovec([x^2 * y, y^2])) + if even_odd_separated + exps = [(2, 1), (0, 2), (0, 1), (0, 0)] + else + exps = [(2, 1), (2, 0), (1, 1), (0, 2), (1, 0), (0, 1), (0, 0)] + end + for i in 1:length(monos) + @test monos.polynomials[i] == univariate_x[exps[i][1] + 1] * univariate_y[exps[i][2] + 1] + end + monos = basis_covering_monomials(B, monovec([x^4, x^2, x])) + if even_odd_separated + exps = [4, 2, 1, 0] + else + exps = [4, 3, 2, 1, 0] + end + for i in 1:length(monos) + @test monos.polynomials[i] == univariate_x[exps[i] + 1] + end + end end @testset "Monomial" begin @@ -30,6 +69,15 @@ end @testset "Fixed" begin include("fixed.jl") end +@testset "Hermite" begin + include("hermite.jl") +end +@testset "Laguerre" begin + include("laguerre.jl") +end +@testset "Legendre" begin + include("legendre.jl") +end @testset "Chebyshev" begin include("chebyshev.jl") end