Improve singularities behaviour

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -201,21 +201,19 @@ recurrencecoefficients(Q) = recurrencecoefficients_layout(MemoryLayout(Q), Q)
 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_layout(lay::BroadcastLayout, a) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, 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(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)
 basis_axes(ax::Inclusion{<:Any,<:AbstractInterval}, v) = convert(AbstractQuasiMatrix{eltype(v)}, basis_singularities(ax, singularities(v)))
 
 

src/classical/jacobi.jl

@@ -54,23 +54,24 @@ hasboundedendpoints(w::AbstractJacobiWeight) = w.a ≥ 0 && w.b ≥ 0
 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...)
 end
 
-singularitiesbroadcast(::typeof(*), V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(*, map(_parent,V)...)[_parentindices(V...)...]
+singularitiesbroadcast(::typeof(*), V::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...)...]
+singularitiesbroadcast(F, V::SubQuasiArray...) = singularitiesbroadcast(F, map(parent,V)...)[parentindices(V...)...]
+singularitiesbroadcast(F, V::NoSingularities...) = NoSingularities() # default is to assume smooth
 
 
 abstract type AbstractJacobi{T} <: OrthogonalPolynomial{T} end
@@ -107,10 +108,11 @@ AbstractQuasiMatrix{T}(w::Jacobi) where T = Jacobi{T}(w.a, w.b)
 
 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

@@ -38,13 +38,10 @@ end
 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(::LegendreWeight{T}) where T = Legendre{T}()
-basis_singularities(v::SubQuasiArray) = view(basis_singularities(parent(v)), parentindices(v)[1], :)
+singularities(::AbstractFillLayout, P) = LegendreWeight{eltype(P)}()
+singularities(::Legendre) = NoSingularities()
+basis_singularities(ax::Inclusion, ::NoSingularities) = legendre(ax)
 
 struct Legendre{T} <: AbstractJacobi{T} end
 Legendre() = Legendre{Float64}()