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,44 @@ function basis_covering_monomials(
         MP.monomial_vector(collect(m)),
     )
 end
+
+function scalar_product_function(::Type{<:AbstractMultipleOrthogonalBasis}, i::Int) end
+
+LinearAlgebra.dot(p, q, basis_type::Type{<:AbstractMultipleOrthogonalBasis}) =
+    _integral(p * q, basis_type)
+
+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

@@ -76,7 +76,7 @@ function change_basis(
 ) where {MT,MV}
     n = length(basis)
     scalings = map(scaling, basis.monomials)
-    scaled_Q = [Q[i, j] * scalings[i] * scalings[j] for i in 1:n, j in 1:n]
+    scaled_Q = [Q[i, j] * scalings[i] * scalings[j] for i = 1:n, j = 1:n]
     return scaled_Q, MonomialBasis(basis.monomials)
 end
 
@@ -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/chebyshev.jl

@@ -7,14 +7,29 @@ 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,32 @@ 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

test/runtests.jl

@@ -1,9 +1,10 @@
 using Test
 
 using MultivariateBases
-
+using LinearAlgebra
 using DynamicPolynomials
 
+
 function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
     @polyvar x[1:2]
     for basis in [
@@ -27,6 +28,17 @@ function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
     end
 end
 
+function univ_orthogonal_test(B::Type{<:AbstractMultipleOrthogonalBasis}, univ::Function; kwargs...)
+    @polyvar x
+    basis = maxdegree_basis(B, [x], 4)
+    for i = eachindex(basis)
+        @test isapprox(dot(basis[i], basis[i], B), univ(maxdegree(basis[i])); kwargs...)
+        for j = 1:i-1
+            @test isapprox(dot(basis[i], basis[j], B), 0.0; kwargs...)
+        end
+    end
+end
+
 function orthogonal_test(
     B::Type{<:AbstractMultipleOrthogonalBasis},
     univ::Function,
@@ -68,6 +80,15 @@ function orthogonal_test(
     end
 end
 
+function coefficient_test(B::Type{<:AbstractPolynomialBasis}, coefs; kwargs...)
+    @polyvar x y
+    p = x^4 * y^2 + x^2 * y^4 - 3 * x^2 * y^2 + 1
+    basis = basis_covering_monomials(B, monomials(p))
+    cc = coefficients(p, basis)
+    @test isapprox(coefs, cc; kwargs...)
+    @test isapprox(p, polynomial(cc, basis); kwargs...)
+end
+
 @testset "Monomial" begin
     include("monomial.jl")
 end

test/scaled.jl

@@ -29,3 +29,18 @@ end
 @testset "API degree = $degree" for degree in 0:3
     api_test(ScaledMonomialBasis, degree)
 end
+
+@testset "Enumerate" begin
+    monos = [1, y, x]
+    basis = ScaledMonomialBasis(monos)
+    for (i, e) in enumerate(basis)
+        @test e == monos[i]
+    end
+end
+
+@testset "Coefficients" begin
+    coefficient_test(
+        ScaledMonomialBasis,
+        [1, -√3/√2, 1/√15, 1/√15],
+    )
+end