Improve singularities behaviour

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.15.1"
+version = "0.15.2"
 
 
 [deps]

src/ClassicalOrthogonalPolynomials.jl

@@ -37,7 +37,7 @@
 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,
+    checkpoints, weight, unweighted, MappedBasisLayouts, sum_layout, invmap, plan_ldiv, layout_broadcasted, MappedBasisLayout, MappedWeightLayout, SubBasisLayout, broadcastbasis_layout,
     plan_grid_transform, plan_transform, MAX_PLOT_POINTS, MulPlan, grammatrix, AdjointBasisLayout, grammatrix_layout, plan_transform_layout, _cumsum
 import FastTransforms: Λ, ChebyshevGrid, chebyshevpoints, Plan, ScaledPlan, th_cheb2leg, pochhammer
 import RecurrenceRelationships: forwardrecurrence, forwardrecurrence!, clenshaw, clenshaw!,
@@ -197,27 +197,45 @@
 
 
 
+singularities_layout(lay::BroadcastLayout, a) = singularitiesbroadcast(call(lay, a), map(singularities, arguments(lay, a))...)
+singularities_layout(::WeightedBasisLayouts, a) = singularities_layout(BroadcastLayout{typeof(*)}(), a)
+singularities_layout(::MappedWeightLayout, a) = view(singularities(demap(a)), basismap(a))
+singularities_layout(::WeightedOPLayout, a) = singularities(weight(a))
+singularities_layout(::ExpansionLayout, f) = singularities(basis(f))
+singularities_layout(lay, a) = NoSingularities() # assume no singularities
+
 """
     singularities(f)
 
 gives the singularity structure of an expansion, e.g.,
 `JacobiWeight`.
 """
-singularities(::AbstractWeightLayout, w) = w
-singularities(lay::BroadcastLayout, a) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, a))...)
-singularities(::WeightedBasisLayouts, a) = singularities(BroadcastLayout{typeof(*)}(), a)
-singularities(::WeightedOPLayout, a) = singularities(weight(a))
-singularities(w) = singularities(MemoryLayout(w), w)
-singularities(::ExpansionLayout, f) = singularities(basis(f))
+singularities(w) = singularities_layout(MemoryLayout(w), w)
+
 
-singularitiesview(w, ::Inclusion) = w # for now just assume it doesn't change
-singularitiesview(w, ind) = view(w, ind)
-singularities(S::SubQuasiArray) = singularitiesview(singularities(parent(S)), parentindices(S)[1])
 
 struct NoSingularities end
 
-basis_singularities(ax, ::NoSingularities) = basis(ax)
-basis_singularities(ax, sing) = basis_singularities(sing)
+## default is to just assume no singularities
+singularitiesbroadcast(_...) = NoSingularities()
+
+for op in (:+, :*)
+    @eval singularitiesbroadcast(::typeof($op), A, B, C, D...) = singularitiesbroadcast(*, singularitiesbroadcast(*, A, B), C, D...)
+    @eval singularitiesbroadcast(::typeof($op), ::NoSingularities, ::NoSingularities) = NoSingularities()
+end
+
+singularitiesbroadcast(::typeof(*), ::NoSingularities, b) = b
+singularitiesbroadcast(::typeof(*), a, ::NoSingularities) = a
+
+
+# for singularitiesbroadcast(literal_pow), ^, ...)
+# singularitiesbroadcast(F::Function, G::Function, V::SubQuasiArray, K) = singularitiesbroadcast(F, G, parent(V), K)[parentindices(V)...]
+# singularitiesbroadcast(F, V::SubQuasiArray...) = singularitiesbroadcast(F, map(parent,V)...)[parentindices(V...)...]
+# singularitiesbroadcast(F, V::NoSingularities...) = NoSingularities() # default is to assume smooth
+
+
+
+
 basis_axes(ax::Inclusion{<:Any,<:AbstractInterval}, v) = convert(AbstractQuasiMatrix{eltype(v)}, basis_singularities(ax, singularities(v)))
 
 

src/choleskyQR.jl

@@ -26,7 +26,7 @@
 # transform to P * inv(U) if needed for differentiation, etc.
 arguments(::ApplyLayout{typeof(*)}, Q::ConvertedOrthogonalPolynomial) = Q.P, ApplyArray(inv, Q.U)
 
-OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthogonalpolynomial(singularities(w)))
+OrthogonalPolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w, orthogonalpolynomial_axis(axes(w,1), singularities(w)))
 function OrthogonalPolynomial(w::AbstractQuasiVector, P::AbstractQuasiMatrix)
     Q = normalized(P)
     X,U = cholesky_jacobimatrix(w, Q)
@@ -39,6 +39,8 @@
 OrthogonalPolynomial(w::Function, P::AbstractQuasiMatrix) = OrthogonalPolynomial(w.(axes(P,1)), P)
 orthogonalpolynomial(w::Function, P::AbstractQuasiMatrix) = orthogonalpolynomial(w.(axes(P,1)), P)
 
+orthogonalpolynomial_axis(ax, ::NoSingularities) = legendre(ax)
+orthogonalpolynomial_axis(ax, w) = orthogonalpolynomial(w)
 resizedata!(P::ConvertedOrthogonalPolynomial, ::Colon, n::Int) = resizedata!(P.X.dv, n)
 
 

src/classical/chebyshev.jl

@@ -49,10 +49,12 @@
 
 chebyshevt() = ChebyshevT()
 chebyshevt(d::AbstractInterval{T}) where T = ChebyshevT{float(T)}()[affine(d, ChebyshevInterval{T}()), :]
+chebyshevt(d::ChebyshevInterval{T}) where T = ChebyshevT{float(T)}()
 chebyshevt(d::Inclusion) = chebyshevt(d.domain)
 chebyshevt(S::AbstractQuasiMatrix) = chebyshevt(axes(S,1))
 chebyshevu() = ChebyshevU()
 chebyshevu(d::AbstractInterval{T}) where T = ChebyshevU{float(T)}()[affine(d, ChebyshevInterval{T}()), :]
+chebyshevu(d::ChebyshevInterval{T}) where T = ChebyshevU{float(T)}()
 chebyshevu(d::Inclusion) = chebyshevu(d.domain)
 chebyshevu(S::AbstractQuasiMatrix) = chebyshevu(axes(S,1))
 

src/classical/jacobi.jl

@@ -92,23 +92,7 @@
 singularities(a::AbstractAffineQuasiVector) = singularities(a.x)
 
 
-## default is to just assume no singularities
-singularitiesbroadcast(_...) = NoSingularities()
-
-for op in (:+, :*)
-    @eval singularitiesbroadcast(::typeof($op), A, B, C, D...) = singularitiesbroadcast(*, singularitiesbroadcast(*, A, B), C, D...)
-end
-
-singularitiesbroadcast(::typeof(*), V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(*, map(_parent,V)...)[_parentindices(V...)...]
-
-
-_parent(::NoSingularities) = NoSingularities()
-_parent(a) = parent(a)
-_parentindices(a::NoSingularities, b...) = _parentindices(b...)
-_parentindices(a, b...) = parentindices(a)
-# for singularitiesbroadcast(literal_pow), ^, ...)
-singularitiesbroadcast(F::Function, G::Function, V::SubQuasiArray, K) = singularitiesbroadcast(F, G, parent(V), K)[parentindices(V)...]
-singularitiesbroadcast(F, V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(F, map(_parent,V)...)[_parentindices(V...)...]
+singularities(w::AbstractJacobiWeight) = w
 
 
 abstract type AbstractJacobi{T} <: OrthogonalPolynomial{T} end
@@ -208,10 +192,11 @@
 """
 jacobi(a,b) = Jacobi(a,b)
 jacobi(a,b, d::AbstractInterval{T}) where T = Jacobi{float(promote_type(eltype(a),eltype(b),T))}(a,b)[affine(d,ChebyshevInterval{T}()), :]
+jacobi(a,b, d::ChebyshevInterval{T}) where T = Jacobi{float(promote_type(eltype(a),eltype(b),T))}(a,b)
 
 Jacobi(P::Legendre{T}) where T = Jacobi(zero(T), zero(T))
 
-basis_singularities(w::JacobiWeight) = Weighted(Jacobi(w.a, w.b))
+
 
 
 """

src/classical/legendre.jl

@@ -73,13 +73,9 @@
 singularitiesbroadcast(::typeof(^), L::LegendreWeight, ::NoSingularities) = L
 singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L # can't find roots
 
-singularities(::AbstractJacobi{T}) where T = LegendreWeight{T}()
-singularities(::Inclusion{T,<:ChebyshevInterval}) where T = LegendreWeight{T}()
-singularities(d::Inclusion{T,<:AbstractInterval}) where T = LegendreWeight{T}()[affine(d,ChebyshevInterval{T}())]
-singularities(::AbstractFillLayout, P) = LegendreWeight{eltype(P)}()
+basis_singularities(ax::Inclusion, ::NoSingularities) = legendre(ax)
+basis_singularities(ax, sing) = basis_singularities(sing) # fallback for back compatibility
 
-basis_singularities(::LegendreWeight{T}) where T = Legendre{T}()
-basis_singularities(v::SubQuasiArray) = view(basis_singularities(parent(v)), parentindices(v)[1], :)
 
 """
     Legendre{T=Float64}(a,b)

src/classical/ultraspherical.jl

@@ -56,6 +56,7 @@
 
 ultrasphericalc(n::Integer, λ, z) = Base.unsafe_getindex(Ultraspherical{polynomialtype(typeof(λ),typeof(z))}(λ), z, n+1)
 ultraspherical(λ, d::AbstractInterval{T}) where T = Ultraspherical{float(promote_type(eltype(λ),T))}(λ)[affine(d,ChebyshevInterval{T}()), :]
+ultraspherical(λ, d::ChebyshevInterval{T}) where T = Ultraspherical{float(promote_type(eltype(λ),T))}(λ)
 
 ==(a::Ultraspherical, b::Ultraspherical) = a.λ == b.λ
 ==(::Ultraspherical, ::ChebyshevT) = false

src/lanczos.jl

@@ -165,7 +165,7 @@
     LanczosPolynomial(w, Q, LanczosData(w, Q))
 end
 
-LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, normalized(orthogonalpolynomial(singularities(w))))
+LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, normalized(orthogonalpolynomial_axis(axes(w,1), singularities(w))))
 
 ==(A::LanczosPolynomial, B::LanczosPolynomial) = A.w == B.w
 ==(::LanczosPolynomial, ::OrthogonalPolynomial) = false # TODO: fix

test/runtests.jl

@@ -8,7 +8,7 @@ import ClassicalOrthogonalPolynomials: jacobimatrix, ∞, ChebyshevInterval, Leg
 import LazyArrays: ApplyStyle, colsupport, MemoryLayout, arguments
 import SemiseparableMatrices: VcatAlmostBandedLayout
 import QuasiArrays: MulQuasiMatrix
-import ClassicalOrthogonalPolynomials: oneto
+import ClassicalOrthogonalPolynomials: oneto, NoSingularities
 import InfiniteLinearAlgebra: KronTrav, Block
 import FastTransforms: clenshaw!
 
@@ -18,15 +18,15 @@ Random.seed!(0)
     x = Inclusion(ChebyshevInterval())
     @test singularities(x) == singularities(exp.(x)) == singularities(x.^2) == 
         singularities(x .+ 1) == singularities(1 .+ x) == singularities(x .+ x) == 
-        LegendreWeight()
+        NoSingularities()
     @test singularities(exp.(x) .* JacobiWeight(0.1,0.2)) == 
         singularities(JacobiWeight(0.1,0.2) .* exp.(x)) ==
         JacobiWeight(0.1,0.2)
 
     x = Inclusion(0..1)
     @test singularities(x) == singularities(exp.(x)) == singularities(x.^2) == 
         singularities(x .+ 1) == singularities(1 .+ x) == singularities(x .+ x) == 
-        legendreweight(0..1)
+        NoSingularities()
 end
 
 include("test_chebyshev.jl")

test/test_chebyshev.jl

@@ -582,6 +582,11 @@ import BandedMatrices: isbanded
     @testset "diff of truncation" begin
         @test MemoryLayout(diff(ChebyshevT()[:,1:5]) * randn(5)) isa ExpansionLayout
     end
+
+    @testset "ChebyshevInterval constructior" begin
+        @test chebyshevt(ChebyshevInterval()) ≡ ChebyshevT()
+        @test chebyshevu(ChebyshevInterval()) ≡ ChebyshevU()
+    end
 end
 
 struct QuadraticMap{T} <: Map{T} end
@@ -612,7 +617,7 @@ ContinuumArrays.invmap(::InvQuadraticMap{T}) where T = QuadraticMap{T}()
     g = chebyshevt(0..1) * [1:3; zeros(∞)]
     @test_broken (f + g)[0.1] ≈ f[0.1] + g[0.1] # ContinuumArrays needs to check maps are equal
 
-    @test ClassicalOrthogonalPolynomials.singularities(T) === LegendreWeight()[QuadraticMap()]
+    @test ClassicalOrthogonalPolynomials.singularities(T) isa NoSingularities
 end
 
 @testset "block structure" begin

test/test_jacobi.jl

@@ -1,5 +1,5 @@
 using ClassicalOrthogonalPolynomials, FillArrays, BandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FastGaussQuadrature, Test
-import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMatrix, arguments, Weighted, HalfWeighted, grammatrix
+import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMatrix, arguments, Weighted, HalfWeighted, grammatrix, singularities
 
 @testset "Jacobi" begin
     @testset "JacobiWeight" begin
@@ -552,4 +552,13 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         @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
+
+    @testset "JacobiWeight singularities" begin
+        @test singularities(JacobiWeight(1,2) .* Jacobi(2,3)) == JacobiWeight(1,2)
+        @test singularities(jacobiweight(1,2,1..2) .* jacobi(2,3,1..2)) == singularities(view(JacobiWeight(1,2) .* Jacobi(2,3), affine(1..2, -1..1),:)) == jacobiweight(1,2,1..2) 
+    end
+
+    @testset "ChebyshevInterval constructior" begin
+        @test jacobi(1,2,ChebyshevInterval()) ≡ Jacobi(1,2)
+    end
 end

test/test_legendre.jl

@@ -217,5 +217,16 @@ import QuasiArrays: MulQuasiArray
         @test (P¹ \ diff(P,1))[1:10,1:10] == (P¹ \ diff(P))[1:10,1:10]
         @test (P² \ diff(P,2))[1:10,1:10] ≈ (P² \ diff(diff(P)))[1:10,1:10]
         @test (P³ \ diff(P,3))[1:10,1:10] ≈ (P³ \ diff(diff(diff(P))))[1:10,1:10]
+
+        # cover back compatibilty
+        @test_throws MethodError ClassicalOrthogonalPolynomials.basis_singularities(nothing, nothing)
+    end
+
+    @testset "singularities expand" begin
+        x = Inclusion(1..2)
+        @test basis(expand(fill(2, x))) == basis(expand(broadcast(+, exp.(x), cos.(x), x))) == legendre(x) 
+    end
+    @testset "ChebyshevInterval constructor" begin
+        @test legendre(ChebyshevInterval()) ≡ Legendre()
     end
 end

test/test_monic.jl

@@ -2,57 +2,59 @@ using ClassicalOrthogonalPolynomials
 using Test
 using ClassicalOrthogonalPolynomials: Monic, _p0, orthogonalityweight, recurrencecoefficients
 
-@testset "Basic definition" begin
-    P1 = Legendre()
-    P2 = Normalized(P1)
-    P3 = Monic(P1)
-    @test P3.P == P2
-    @test Monic(P3) === P3
-    @test axes(P3) == axes(Legendre())
-    @test Normalized(P3) === P3.P
-    @test _p0(P3) == 1
-    @test orthogonalityweight(P3) == orthogonalityweight(P1)
-    @test sprint(show, MIME"text/plain"(), P3) == "Monic(Legendre())"
-end
+@testset "Monic" begin
+    @testset "Basic definition" begin
+        P1 = Legendre()
+        P2 = Normalized(P1)
+        P3 = Monic(P1)
+        @test P3.P == P2
+        @test Monic(P3) === P3
+        @test axes(P3) == axes(Legendre())
+        @test Normalized(P3) === P3.P
+        @test _p0(P3) == 1
+        @test orthogonalityweight(P3) == orthogonalityweight(P1)
+        @test sprint(show, MIME"text/plain"(), P3) == "Monic(Legendre())"
+    end
 
-@testset "evaluation" begin
-    function _pochhammer(x, n)
-        y = one(x)
-        for i in 0:(n-1)
-            y *= (x + i)
+    @testset "evaluation" begin
+        function _pochhammer(x, n)
+            y = one(x)
+            for i in 0:(n-1)
+                y *= (x + i)
+            end
+            return y
+        end
+        jacobi_kn = (α, β, n) -> _pochhammer(n + α + β + 1, n) / (2.0^n * factorial(n))
+        ultra_kn = (λ, n) -> 2^n * _pochhammer(λ, n) / factorial(n)
+        chebt_kn = n -> n == 0 ? 1.0 : 2.0 .^ (n - 1)
+        chebu_kn = n -> 2.0^n
+        leg_kn = n -> 2.0^n * _pochhammer(1 / 2, n) / factorial(n)
+        lag_kn = n -> (-1)^n / factorial(n)
+        herm_kn = n -> 2.0^n
+        _Jacobi(α, β, x, n) = Jacobi(α, β)[x, n+1] / jacobi_kn(α, β, n)
+        _Ultraspherical(λ, x, n) = Ultraspherical(λ)[x, n+1] / ultra_kn(λ, n)
+        _ChebyshevT(x, n) = ChebyshevT()[x, n+1] / chebt_kn(n)
+        _ChebyshevU(x, n) = ChebyshevU()[x, n+1] / chebu_kn(n)
+        _Legendre(x, n) = Legendre()[x, n+1] / leg_kn(n)
+        _Laguerre(α, x, n) = Laguerre(α)[x, n+1] / lag_kn(n)
+        _Hermite(x, n) = Hermite()[x, n+1] / herm_kn(n)
+        Ps = [
+            Jacobi(2.0, 5.0) (x, n)->_Jacobi(2.0, 5.0, x, n)
+            Ultraspherical(1.7) (x, n)->_Ultraspherical(1.7, x, n)
+            ChebyshevT() _ChebyshevT
+            ChebyshevU() _ChebyshevU
+            Legendre() _Legendre
+            Laguerre(1.5) (x, n)->_Laguerre(1.5, x, n)
+            Hermite() _Hermite
+        ]
+        for (P, _P) in eachrow(Ps)
+            Q = Monic(P)
+            @test Q[0.2, 1] == 1.0
+            @test Q[0.25, 2] ≈ _P(0.25, 1)
+            @test Q[0.17, 3] ≈ _P(0.17, 2)
+            @test Q[0.4, 17] ≈ _P(0.4, 16)
+            @test Q[0.9, 21] ≈ _P(0.9, 20)
+            # @inferred Q[0.2, 5] # no longer inferred
         end
-        return y
-    end
-    jacobi_kn = (α, β, n) -> _pochhammer(n + α + β + 1, n) / (2.0^n * factorial(n))
-    ultra_kn = (λ, n) -> 2^n * _pochhammer(λ, n) / factorial(n)
-    chebt_kn = n -> n == 0 ? 1.0 : 2.0 .^ (n - 1)
-    chebu_kn = n -> 2.0^n
-    leg_kn = n -> 2.0^n * _pochhammer(1 / 2, n) / factorial(n)
-    lag_kn = n -> (-1)^n / factorial(n)
-    herm_kn = n -> 2.0^n
-    _Jacobi(α, β, x, n) = Jacobi(α, β)[x, n+1] / jacobi_kn(α, β, n)
-    _Ultraspherical(λ, x, n) = Ultraspherical(λ)[x, n+1] / ultra_kn(λ, n)
-    _ChebyshevT(x, n) = ChebyshevT()[x, n+1] / chebt_kn(n)
-    _ChebyshevU(x, n) = ChebyshevU()[x, n+1] / chebu_kn(n)
-    _Legendre(x, n) = Legendre()[x, n+1] / leg_kn(n)
-    _Laguerre(α, x, n) = Laguerre(α)[x, n+1] / lag_kn(n)
-    _Hermite(x, n) = Hermite()[x, n+1] / herm_kn(n)
-    Ps = [
-        Jacobi(2.0, 5.0) (x, n)->_Jacobi(2.0, 5.0, x, n)
-        Ultraspherical(1.7) (x, n)->_Ultraspherical(1.7, x, n)
-        ChebyshevT() _ChebyshevT
-        ChebyshevU() _ChebyshevU
-        Legendre() _Legendre
-        Laguerre(1.5) (x, n)->_Laguerre(1.5, x, n)
-        Hermite() _Hermite
-    ]
-    for (P, _P) in eachrow(Ps)
-        Q = Monic(P)
-        @test Q[0.2, 1] == 1.0
-        @test Q[0.25, 2] ≈ _P(0.25, 1)
-        @test Q[0.17, 3] ≈ _P(0.17, 2)
-        @test Q[0.4, 17] ≈ _P(0.4, 16)
-        @test Q[0.9, 21] ≈ _P(0.9, 20)
-        # @inferred Q[0.2, 5] # no longer inferred
     end
 end

test/test_ultraspherical.jl

@@ -204,4 +204,7 @@ using ClassicalOrthogonalPolynomials: grammatrix
         @test (C \ diff(U,1))[1:10,1:10] == (C \ diff(U))[1:10,1:10]
         @test (C³ \ diff(U,2))[1:10,1:10] == (C³ \ diff(diff(U)))[1:10,1:10]
     end
+    @testset "ChebyshevInterval constructior" begin
+        @test ultraspherical(2,ChebyshevInterval()) ≡ Ultraspherical(2)
+    end
 end