diff --git a/Project.toml b/Project.toml index be1bf1b..4a65c44 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "ClassicalOrthogonalPolynomials" uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" authors = ["Sheehan Olver "] -version = "0.14.4" +version = "0.15" [deps] @@ -38,7 +38,7 @@ ArrayLayouts = "1.3.1" BandedMatrices = "1" BlockArrays = "1" BlockBandedMatrices = "0.13" -ContinuumArrays = "0.18.3" +ContinuumArrays = "0.19" DomainSets = "0.6, 0.7" DynamicPolynomials = "0.6" FFTW = "1.1" @@ -46,15 +46,15 @@ FastGaussQuadrature = "1" FastTransforms = "0.16.6, 0.17" FillArrays = "1" HypergeometricFunctions = "0.3.4" -InfiniteArrays = " 0.14, 0.15" -InfiniteLinearAlgebra = "0.8, 0.9" +InfiniteArrays = " 0.15" +InfiniteLinearAlgebra = "0.9" IntervalSets = "0.7" LazyArrays = "2.2" -LazyBandedMatrices = "0.10, 0.11" +LazyBandedMatrices = "0.11" MutableArithmetics = "1" -QuasiArrays = "0.11" +QuasiArrays = "0.12" RecurrenceRelationshipArrays = "0.1.2" -RecurrenceRelationships = "0.1.1, 0.2" +RecurrenceRelationships = "0.2" SpecialFunctions = "1.0, 2" julia = "1.10" diff --git a/src/ClassicalOrthogonalPolynomials.jl b/src/ClassicalOrthogonalPolynomials.jl index 380145d..f0b688b 100644 --- a/src/ClassicalOrthogonalPolynomials.jl +++ b/src/ClassicalOrthogonalPolynomials.jl @@ -34,7 +34,7 @@ import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclu import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges import InfiniteLinearAlgebra: chop!, chop, pad, choplength, compatible_resize!, partialcholesky! -import ContinuumArrays: Basis, Weight, basis_axes, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization, +import ContinuumArrays: Basis, Weight, basis_axes, @simplify, AbstractAffineQuasiVector, ProjectionFactorization, grid, plotgrid, plotgrid_layout, plotvalues_layout, grid_layout, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, ExpansionLayout, basismap, AffineQuasiVector, AffineMap, AbstractWeightLayout, AbstractWeightedBasisLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout, checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, SubBasisLayout, broadcastbasis_layout, @@ -90,6 +90,8 @@ represents an OP multiplied by its orthogonality weight. """ struct WeightedOPLayout{Lay<:AbstractOPLayout} <: AbstractWeightedBasisLayout end +grid_layout(::WeightedOPLayout, P, n) = grid(unweighted(P), n) + isorthogonalityweighted(::WeightedOPLayout, _) = true function isorthogonalityweighted(::AbstractWeightedBasisLayout, wS) w,S = arguments(wS) @@ -246,7 +248,7 @@ grammatrix_layout(::WeightedOPLayout{MappedOPLayout}, P) = grammatrix_layout(Map OrthogonalPolynomial(w::Weight) =error("Override for $(typeof(w))") -@simplify *(B::Identity, C::OrthogonalPolynomial) = ApplyQuasiMatrix(*, C, jacobimatrix(C)) +@simplify *(B::QuasiDiagonal{<:Any,<:Inclusion}, C::OrthogonalPolynomial) = ApplyQuasiMatrix(*, C, jacobimatrix(C)) function layout_broadcasted(::Tuple{PolynomialLayout,AbstractOPLayout}, ::typeof(*), x::Inclusion, C) x == axes(C,1) || throw(DimensionMismatch()) diff --git a/src/classical/jacobi.jl b/src/classical/jacobi.jl index b5ef4c7..7067312 100644 --- a/src/classical/jacobi.jl +++ b/src/classical/jacobi.jl @@ -38,13 +38,16 @@ julia> axes(J) (Inclusion(-1.0 .. 1.0 (Chebyshev)),) ``` """ -struct JacobiWeight{T} <: AbstractJacobiWeight{T} - a::T - b::T - JacobiWeight{T}(a, b) where T = new{T}(convert(T,a), convert(T,b)) +struct JacobiWeight{T,V} <: AbstractJacobiWeight{T} + a::V + b::V + JacobiWeight{T,V}(a, b) where {T,V} = new{T,V}(convert(V,a), convert(V,b)) end -JacobiWeight(a::V, b::T) where {T,V} = JacobiWeight{promote_type(T,V)}(a,b) +JacobiWeight{T}(a::V, b::V) where {T,V} = JacobiWeight{T,V}(a, b) +JacobiWeight{T}(a, b) where T = JacobiWeight{T}(promote(a,b)...) +JacobiWeight(a::V, b::T) where {T,V} = JacobiWeight{float(promote_type(T,V))}(a, b) + """ jacobiweight(a,b, d::AbstractInterval) @@ -71,9 +74,9 @@ AbstractQuasiVector{T}(w::JacobiWeight) where T = JacobiWeight{T}(w.a, w.b) ==(A::JacobiWeight, B::JacobiWeight) = A.b == B.b && A.a == B.a -function getindex(w::JacobiWeight, x::Number) +function getindex(w::JacobiWeight{T}, x::Number) where T x ∈ axes(w,1) || throw(BoundsError()) - (1-x)^w.a * (1+x)^w.b + convert(T, (1-x)^w.a * (1+x)^w.b) end show(io::IO, P::JacobiWeight) = summary(io, P) @@ -479,6 +482,8 @@ end broadcastbasis(::typeof(+), w_A::Weighted{<:Any,<:Jacobi}, w_B::Weighted{<:Any,<:Jacobi}) = broadcastbasis(+, convert(WeightedBasis,w_A), convert(WeightedBasis,w_B)) broadcastbasis(::typeof(+), w_A::Weighted{<:Any,<:Jacobi}, w_B::WeightedJacobi) = broadcastbasis(+, convert(WeightedBasis,w_A), w_B) broadcastbasis(::typeof(+), w_A::WeightedJacobi, w_B::Weighted{<:Any,<:Jacobi}) = broadcastbasis(+, w_A, convert(WeightedBasis,w_B)) +broadcastbasis(::typeof(+), A::Jacobi, B::Weighted{<:Any,<:Jacobi{<:Any,<:Integer}}) = A +broadcastbasis(::typeof(+), A::Weighted{<:Any,<:Jacobi{<:Any,<:Integer}}, B::Jacobi) = B function \(w_A::WeightedJacobi, w_B::WeightedJacobi) wA,A = w_A.args diff --git a/src/classical/ultraspherical.jl b/src/classical/ultraspherical.jl index 9cff369..7b75f0e 100644 --- a/src/classical/ultraspherical.jl +++ b/src/classical/ultraspherical.jl @@ -139,7 +139,9 @@ end # higher order -function diff(::ChebyshevT{T}, m::Integer; dims=1) where T +function diff(S::ChebyshevT{T}, m::Integer; dims=1) where T + iszero(m) && return S + isone(m) && return diff(S) μ = pochhammer(one(T),m-1)*convert(T,2)^(m-1) D = _BandedMatrix((μ * (0:∞))', ℵ₀, -m, m) ApplyQuasiMatrix(*, Ultraspherical{T}(m), D) diff --git a/src/lanczos.jl b/src/lanczos.jl index 9a5417b..5b6f401 100644 --- a/src/lanczos.jl +++ b/src/lanczos.jl @@ -233,4 +233,4 @@ broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q::Union{Nor broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q) = broadcastbasis(+, P.P, Q) broadcastbasis(::typeof(+), P, Q::Union{Normalized,LanczosPolynomial}) = broadcastbasis(+, P, Q.P) -diff_layout(::LanczosLayout, P::AbstractQuasiMatrix, dims...) = diff_layout(ApplyLayout{typeof(*)}(), P, dims...) \ No newline at end of file +diff_layout(::LanczosLayout, P::AbstractQuasiMatrix, order...; dims...) = diff_layout(ApplyLayout{typeof(*)}(), P, order...; dims...) \ No newline at end of file diff --git a/src/normalized.jl b/src/normalized.jl index b4d4419..36f2fc8 100644 --- a/src/normalized.jl +++ b/src/normalized.jl @@ -299,6 +299,7 @@ _sum(p::SubQuasiArray{T,1,<:Weighted,<:Tuple{Inclusion,Int}}, ::Colon) where T = demap(W::Weighted) = Weighted(demap(W.P)) basismap(W::Weighted) = basismap(W.P) const MappedOPLayouts = Union{MappedOPLayout,WeightedOPLayout{MappedOPLayout}} -diff_layout(::MappedOPLayouts, A, dims...) = diff_layout(MappedBasisLayout(), A, dims...) +diff_layout(::MappedOPLayouts, A, order::Int; dims...) = diff_layout(MappedBasisLayout(), A, order; dims...) +diff_layout(::MappedOPLayouts, A, order...; dims...) = diff_layout(MappedBasisLayout(), A, order...; dims...) -diff_layout(::AbstractNormalizedOPLayout, A, dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, dims...) \ No newline at end of file +diff_layout(::AbstractNormalizedOPLayout, A, order...; dims...) = diff_layout(ApplyLayout{typeof(*)}(), A, order...; dims...) diff --git a/test/test_jacobi.jl b/test/test_jacobi.jl index cdb8d9a..b2b85e3 100644 --- a/test/test_jacobi.jl +++ b/test/test_jacobi.jl @@ -537,4 +537,19 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa @testset "conversion not implemented" begin @test_throws ArgumentError Jacobi(0,0) \ Jacobi(1.1,2.1) end + + @testset "broadcastbasis" begin + a = Jacobi(1,1) * [1; zeros(∞)] + b = Weighted(Jacobi(1,1)) * [1; zeros(∞)] + c = Weighted(Jacobi(2,3)) * [1; zeros(∞)] + + @test basis(a+b) == basis(b+a) == basis(a+c) == basis(c+a) == Jacobi(1,1) + @test a[0.1]+b[0.1] ≈ (a+b)[0.1] ≈ (b+a)[0.1] + end + + @testset "Weighted expand" begin + W = Weighted(Jacobi(1,1)) + @test expand(W, x -> (1-x^2)*exp(x))[0.1] ≈ (1-0.1^2)*exp(0.1) + @test grid(W, 5) == grid(W.P, 5) + end end \ No newline at end of file diff --git a/test/test_lanczos.jl b/test/test_lanczos.jl index 95b5290..65faa05 100644 --- a/test/test_lanczos.jl +++ b/test/test_lanczos.jl @@ -303,16 +303,23 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedColumns, or @test jacobimatrix(Q)[1,1] ≈ 1/3 @test Q[0.5,1:3] ≈ [1, 1.369306393762913, 0.6469364618834543] end -end - -@testset "#197" begin - @test _emptymaximum(1:5) == 5 - @test _emptymaximum(1:0) == 0 - x = Inclusion(ChebyshevInterval()) - f = exp.(x) - QQ = LanczosPolynomial(f) - R = LanczosConversion(QQ.data) - v = cache(Zeros(∞)) - @test (R \ v)[1:500] == zeros(500) - @test (R * v)[1:500] == zeros(500) + + @testset "#197" begin + @test _emptymaximum(1:5) == 5 + @test _emptymaximum(1:0) == 0 + x = Inclusion(ChebyshevInterval()) + f = exp.(x) + QQ = LanczosPolynomial(f) + R = LanczosConversion(QQ.data) + v = cache(Zeros(∞)) + @test (R \ v)[1:500] == zeros(500) + @test (R * v)[1:500] == zeros(500) + end + + @testset "diff" begin + P = Normalized(Legendre()) + x = axes(P,1) + Q = LanczosPolynomial(exp.(x)) + @test diff(Q)[0.1,3] ≈ -0.1637907411174539 + end end \ No newline at end of file