work on singularities

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -205,14 +205,17 @@ singularities_layout(lay::BroadcastLayout, a) = singularitiesbroadcast(call(a),
 singularities_layout(::WeightedBasisLayouts, a) = singularities(BroadcastLayout{typeof(*)}(), a)
 singularities_layout(::WeightedOPLayout, a) = singularities(weight(a))
 singularities_layout(::ExpansionLayout, f) = singularities(basis(f))
-singularities_layout(::PolynomialLayout, a) = NoSingularities()
+singularities_layout(lay, a) = NoSingularities() # assume no singularities
 singularities(w) = singularities_layout(MemoryLayout(w), w)
 
+struct NoSingularities end
+
 singularitiesview(w, ::Inclusion) = w # for now just assume it doesn't change
 singularitiesview(w, ind) = view(w, ind)
+singularitiesview(::NoSingularities, ind) = NoSingularities()
+singularitiesview(::NoSingularities, ::Inclusion) = NoSingularities()
 singularities(S::SubQuasiArray) = singularitiesview(singularities(parent(S)), parentindices(S)[1])
 
-struct NoSingularities end
 
 basis_axes(ax::Inclusion{<:Any,<:AbstractInterval}, v) = convert(AbstractQuasiMatrix{eltype(v)}, basis_singularities(ax, singularities(v)))
 

src/choleskyQR.jl

@@ -19,7 +19,7 @@ orthogonalityweight(Q::ConvertedOrthogonalPolynomial) = Q.weight
 # transform to P * 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(axes(w,1), singularities(w)))
 function OrthogonalPolynomial(w::AbstractQuasiVector, P::AbstractQuasiMatrix)
     Q = normalized(P)
     X = cholesky_jacobimatrix(w, Q)
@@ -32,6 +32,9 @@ orthogonalpolynomial(w::SubQuasiArray) = orthogonalpolynomial(parent(w))[parenti
 OrthogonalPolynomial(w::Function, P::AbstractQuasiMatrix) = OrthogonalPolynomial(w.(axes(P,1)), P)
 orthogonalpolynomial(w::Function, P::AbstractQuasiMatrix) = orthogonalpolynomial(w.(axes(P,1)), P)
 
+orthogonalpolynomial(ax, ::NoSingularities) = legendre(ax)
+orthogonalpolynomial(ax, w) = orthogonalpolynomial(w)
+
 
 """
 cholesky_jacobimatrix(w, P)

src/classical/jacobi.jl

@@ -56,16 +56,18 @@ singularities(a::AbstractAffineQuasiVector) = singularities(a.x)
 
 singularities(w::JacobiWeight) = w
 
-
 ## 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()
+    @eval singularitiesbroadcast(::typeof($op), ::NoSingularities, ::NoSingularities, ::NoSingularities) = NoSingularities()
 end
 
 singularitiesbroadcast(::typeof(*), V::SubQuasiArray...) = singularitiesbroadcast(*, map(parent,V)...)[parentindices(V...)...]
-
+singularitiesbroadcast(::typeof(*), ::NoSingularities, b) = b
+singularitiesbroadcast(::typeof(*), a, ::NoSingularities) = a
 
 
 # for singularitiesbroadcast(literal_pow), ^, ...)
@@ -90,9 +92,6 @@ size(P::JacobiTransformPlan, k...) = size(P.chebtransform, k...)
 
 include("legendre.jl")
 
-singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::AbstractJacobiWeight) = b
-singularitiesbroadcast(::typeof(*), a::AbstractJacobiWeight, ::LegendreWeight) = a
-
 
 struct Jacobi{T} <: AbstractJacobi{T}
     a::T

src/classical/legendre.jl

@@ -40,7 +40,6 @@ singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L #
 
 
 singularities(::AbstractFillLayout, P) = LegendreWeight{eltype(P)}()
-singularities(::Legendre) = NoSingularities()
 basis_singularities(ax::Inclusion, ::NoSingularities) = legendre(ax)
 
 struct Legendre{T} <: AbstractJacobi{T} end

src/lanczos.jl

@@ -165,7 +165,7 @@ function LanczosPolynomial(w_in::AbstractQuasiVector, P::AbstractQuasiMatrix)
     LanczosPolynomial(w, Q, LanczosData(w, Q))
 end
 
-LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, normalized(orthogonalpolynomial(singularities(w))))
+LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, normalized(orthogonalpolynomial(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

@@ -612,7 +612,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