Add coefficients for orthogonal bases

Project.toml

@@ -4,6 +4,7 @@ repo = "https://github.com/JuliaAlgebra/MultivariateBases.jl.git"
 version = "0.2.1"
 
 [deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 

src/MultivariateBases.jl

@@ -29,6 +29,8 @@ export generators,
     reccurence_third_coef,
     reccurence_deno_coef
 include("fixed.jl")
+
+import LinearAlgebra
 include("orthogonal.jl")
 include("hermite.jl")
 include("laguerre.jl")

src/chebyshev.jl

@@ -28,6 +28,17 @@ function degree_one_univariate_polynomial(
     MA.@rewrite(variable + 0)
 end
 
+function _scalar_product_function(::Type{<:ChebyshevBasisFirstKind}, i::Int)
+    if i == 0
+        return π
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (π / 2^i) * prod(n+1:i) / factorial(n)
+    end
+end
+
 """
     struct ChebyshevBasisSecondKind{P} <: AbstractChebyshevBasis{P}
         polynomials::Vector{P}
@@ -45,3 +56,14 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(2variable + 0)
 end
+
+function _scalar_product_function(::Type{<:ChebyshevBasisSecondKind}, i::Int)
+    if i == 0
+        return π / 2
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return π / (2^(i + 1)) * prod(n+2:i) / factorial(n)
+    end
+end

src/hermite.jl

@@ -30,6 +30,17 @@ function degree_one_univariate_polynomial(
     MA.@rewrite(1variable)
 end
 
+function _scalar_product_function(::Type{<:ProbabilistsHermiteBasis}, i::Int)
+    if i == 0
+        return √(2 * π)
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (√(2 * π) / (2^n)) * prod(n+1:2*n)
+    end
+end
+
 """
     struct PhysicistsHermiteBasis{P} <: AbstractHermiteBasis{P}
         polynomials::Vector{P}
@@ -48,3 +59,14 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(2variable)
 end
+
+function _scalar_product_function(::Type{<:PhysicistsHermiteBasis}, i::Int)
+    if i == 0
+        return √(π)
+    elseif isodd(i)
+        return 0
+    else
+        n = div(i, 2)
+        return (√(π) / (2^i)) * prod(n+1:i)
+    end
+end

src/laguerre.jl

@@ -26,3 +26,7 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(1 - variable)
 end
+
+function _scalar_product_function(::Type{<:LaguerreBasis}, i::Int)
+    return factorial(i)
+end

src/legendre.jl

@@ -35,3 +35,11 @@ function degree_one_univariate_polynomial(
 )
     MA.@rewrite(variable + 0)
 end
+
+function _scalar_product_function(::Type{<:LegendreBasis}, i::Int)
+    if isodd(i)
+        return 0
+    else
+        return 2 / (i + 1)
+    end
+end

src/monomial.jl

@@ -91,6 +91,10 @@ function MP.polynomial(Q::AbstractMatrix, mb::MonomialBasis, T::Type)
     return MP.polynomial(Q, mb.monomials, T)
 end
 
+function MP.coefficients(p, basis::MonomialBasis)
+    return MP.coefficients(p, basis.monomials)
+end
+
 function MP.coefficients(p, ::Type{<:MonomialBasis})
     return MP.coefficients(p)
 end

src/orthogonal.jl

@@ -158,3 +158,67 @@ function basis_covering_monomials(
         MP.monomial_vector(collect(m)),
     )
 end
+
+function _scalar_product_function(
+    ::Type{<:AbstractMultipleOrthogonalBasis},
+    i::Int,
+) end
+
+function LinearAlgebra.dot(
+    p,
+    q,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return _integral(p * q, basis_type)
+end
+
+function _integral(
+    p::Number,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return p * _scalar_product_function(basis_type, 0)
+end
+
+function _integral(
+    p::MP.AbstractVariable,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return _scalar_product_function(basis_type, 1)
+end
+
+function _integral(
+    p::MP.AbstractMonomial,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return prod([
+        _scalar_product_function(basis_type, i) for i in MP.exponents(p)
+    ])
+end
+
+function _integral(
+    p::MP.AbstractTerm,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return MP.coefficient(p) * _integral(MP.monomial(p), basis_type)
+end
+
+function _integral(
+    p::MP.AbstractPolynomial,
+    basis_type::Type{<:AbstractMultipleOrthogonalBasis},
+)
+    return sum([_integral(t, basis_type) for t in MP.terms(p)])
+end
+
+function MP.coefficients(
+    p,
+    basis::AbstractMultipleOrthogonalBasis;
+    check = true,
+)
+    B = typeof(basis)
+    coeffs = [
+        LinearAlgebra.dot(p, el, B) / LinearAlgebra.dot(el, el, B) for
+        el in basis
+    ]
+    idx = findall(c -> !isapprox(c, 0; atol = 1e-10), coeffs)
+    return coeffs[idx]
+end

src/scaled.jl

@@ -95,3 +95,6 @@ unscale_coef(t::MP.AbstractTerm) = MP.coefficient(t) / scaling(MP.monomial(t))
 function MP.coefficients(p, ::Type{<:ScaledMonomialBasis})
     return unscale_coef.(MP.terms(p))
 end
+function MP.coefficients(p, basis::ScaledMonomialBasis)
+    return MP.coefficients(p, basis.monomials) ./ scaling.(MP.monomials(p))
+end

test/Project.toml

@@ -1,4 +1,5 @@
 [deps]
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

test/chebyshev.jl

@@ -7,14 +7,50 @@ using MultivariateBases
         x -> [1, x, 2x^2 - 1, 4x^3 - 3x, 8x^4 - 8x^2 + 1],
         true,
     )
+    univ_orthogonal_test(
+        ChebyshevBasis,
+        i -> π * (1 / 2 + (0^i) / 2);
+        atol = 1e-12,
+    )
     orthogonal_test(
         ChebyshevBasisSecondKind,
         x -> [1, 2x, 4x^2 - 1, 8x^3 - 4x, 16x^4 - 12x^2 + 1],
         true,
     )
+
+    univ_orthogonal_test(ChebyshevBasisSecondKind, i -> π / 2; atol = 1e-12)
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(ChebyshevBasis, degree)
     api_test(ChebyshevBasisSecondKind, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        ChebyshevBasis,
+        reverse([
+            0.0625,
+            0.0625,
+            0.0625,
+            -0.25,
+            0.0625,
+            -0.3125,
+            -0.3125,
+            0.625,
+        ]),
+    )
+    coefficient_test(
+        ChebyshevBasisSecondKind,
+        reverse([
+            0.015625,
+            0.015625,
+            0.015625,
+            -0.09375,
+            0.015625,
+            -0.109375,
+            -0.109375,
+            0.875,
+        ]),
+    )
+end

test/fixed.jl

@@ -72,3 +72,11 @@ end
     p = @inferred polynomial(zeros(Int, 0, 0), basis, Int)
     @test iszero(p)
 end
+
+@testset "Enumerate" begin
+    monos = [1, x, y^2]
+    basis = FixedPolynomialBasis(monos)
+    for (i, e) in enumerate(basis)
+        @test e == monos[i]
+    end
+end

test/hermite.jl

@@ -7,14 +7,41 @@ using MultivariateBases
         x -> [1, x, x^2 - 1, x^3 - 3x, x^4 - 6x^2 + 3],
         true,
     )
+    univ_orthogonal_test(ProbabilistsHermiteBasis, i -> √(2 * π) * factorial(i))
     orthogonal_test(
         PhysicistsHermiteBasis,
         x -> [1, 2x, 4x^2 - 2, 8x^3 - 12x, 16x^4 - 48x^2 + 12],
         true,
     )
+    univ_orthogonal_test(
+        PhysicistsHermiteBasis,
+        i -> √(π) * factorial(i) * 2^i;
+        atol = 1e-12,
+    ) #precision issue
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(ProbabilistsHermiteBasis, degree)
     api_test(PhysicistsHermiteBasis, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        ProbabilistsHermiteBasis,
+        [4, 6, 6, 1, 9, 1, 1, 1];
+        atol = 1e-12,
+    )
+    coefficient_test(
+        PhysicistsHermiteBasis,
+        reverse([
+            0.015625,
+            0.015625,
+            0.03125,
+            0.1875,
+            0.03125,
+            0.1875,
+            0.1875,
+            1.0,
+        ]),
+    )
+end

test/laguerre.jl

@@ -13,8 +13,39 @@ using MultivariateBases
         ],
         false,
     )
+
+    univ_orthogonal_test(LaguerreBasis, i -> 1; atol = 1e-12) # there are precision issues
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(LaguerreBasis, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        LaguerreBasis,
+        reverse([
+            48.0,
+            48.0,
+            -96.0,
+            -192.0,
+            -192.0,
+            -96.0,
+            48.0,
+            384.0,
+            564.0,
+            384.0,
+            48.0,
+            -192.0,
+            -744.0,
+            -744.0,
+            -192.0,
+            324.0,
+            720.0,
+            324.0,
+            -264.0,
+            -264.0,
+            85.0,
+        ]),
+    )
+end

test/legendre.jl

@@ -7,8 +7,26 @@ using MultivariateBases
         x -> [1, x, (3x^2 - 1) / 2, (5x^3 - 3x) / 2, (35x^4 - 30x^2 + 3) / 8],
         true,
     )
+
+    univ_orthogonal_test(LegendreBasis, i -> 2 / (2 * i + 1))
 end
 
 @testset "API degree = $degree" for degree in 0:3
     api_test(LegendreBasis, degree)
 end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        LegendreBasis,
+        reverse([
+            0.1523809523809524,
+            0.1523809523809524,
+            0.0761904761904762,
+            -0.5714285714285714,
+            0.0761904761904762,
+            -0.3428571428571428,
+            -0.3428571428571428,
+            0.8,
+        ]),
+    )
+end

test/monomial.jl

@@ -66,3 +66,15 @@ end
     p = @inferred polynomial(zeros(Int, 0, 0), basis, Int)
     @test iszero(p)
 end
+
+@testset "Enumerate" begin
+    monos = [1, y, x]
+    basis = MonomialBasis(monos)
+    for (i, e) in enumerate(basis)
+        @test e == monos[i]
+    end
+end
+
+@testset "Coefficients" begin
+    coefficient_test(MonomialBasis, [1, -3, 1, 1])
+end