Skip to content

Commit

Permalink
Fix format
Browse files Browse the repository at this point in the history
  • Loading branch information
blegat committed Nov 28, 2023
1 parent d4ce741 commit 954e83f
Show file tree
Hide file tree
Showing 5 changed files with 84 additions and 24 deletions.
43 changes: 33 additions & 10 deletions src/orthogonal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -159,12 +159,23 @@ function basis_covering_monomials(
)
end

function scalar_product_function(::Type{<:AbstractMultipleOrthogonalBasis}, i::Int) end

LinearAlgebra.dot(p, q, basis_type::Type{<:AbstractMultipleOrthogonalBasis}) =
_integral(p * q, basis_type)
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})
function _integral(
p::Number,
basis_type::Type{<:AbstractMultipleOrthogonalBasis},
)
return p * scalar_product_function(basis_type, 0)
end

Expand All @@ -179,10 +190,15 @@ function _integral(
p::MP.AbstractMonomial,
basis_type::Type{<:AbstractMultipleOrthogonalBasis},
)
return prod([scalar_product_function(basis_type, i) for i in MP.exponents(p)])
return prod([
scalar_product_function(basis_type, i) for i in MP.exponents(p)
])
end

function _integral(p::MP.AbstractTerm, basis_type::Type{<:AbstractMultipleOrthogonalBasis})
function _integral(
p::MP.AbstractTerm,
basis_type::Type{<:AbstractMultipleOrthogonalBasis},
)
return MP.coefficient(p) * _integral(MP.monomial(p), basis_type)
end

Expand All @@ -193,9 +209,16 @@ function _integral(
return sum([_integral(t, basis_type) for t in MP.terms(p)])
end

function MP.coefficients(p, basis::AbstractMultipleOrthogonalBasis; check = true)
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)
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
2 changes: 1 addition & 1 deletion src/scaled.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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 = 1:n, j = 1:n]
scaled_Q = [Q[i, j] * scalings[i] * scalings[j] for i in 1:n, j in 1:n]
return scaled_Q, MonomialBasis(basis.monomials)
end

Expand Down
31 changes: 26 additions & 5 deletions test/chebyshev.jl
Original file line number Diff line number Diff line change
Expand Up @@ -7,29 +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)
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)
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]),
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]),
reverse([
0.015625,
0.015625,
0.015625,
-0.09375,
0.015625,
-0.109375,
-0.109375,
0.875,
]),
)
end
15 changes: 12 additions & 3 deletions test/hermite.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ using MultivariateBases
)
univ_orthogonal_test(
PhysicistsHermiteBasis,
i -> (π) * factorial(i) * 2^i,
i -> (π) * factorial(i) * 2^i;
atol = 1e-12,
) #precision issue
end
Expand All @@ -28,11 +28,20 @@ end
@testset "Coefficients" begin
coefficient_test(
ProbabilistsHermiteBasis,
[4, 6, 6, 1, 9, 1, 1, 1],
[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]),
reverse([
0.015625,
0.015625,
0.03125,
0.1875,
0.03125,
0.1875,
0.1875,
1.0,
]),
)
end
17 changes: 12 additions & 5 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@ using MultivariateBases
using LinearAlgebra
using DynamicPolynomials


function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
@polyvar x[1:2]
for basis in [
Expand All @@ -28,12 +27,20 @@ function api_test(B::Type{<:AbstractPolynomialBasis}, degree)
end
end

function univ_orthogonal_test(B::Type{<:AbstractMultipleOrthogonalBasis}, univ::Function; kwargs...)
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
for i in eachindex(basis)
@test isapprox(
dot(basis[i], basis[i], B),
univ(maxdegree(basis[i]));
kwargs...,
)
for j in 1:i-1
@test isapprox(dot(basis[i], basis[j], B), 0.0; kwargs...)
end
end
Expand Down

0 comments on commit 954e83f

Please sign in to comment.