Fix Lanczos with x^m*(1-x)^m*(2-x)^m (#52)

* Fix Lanczos with x^m*(1-x)^m*(2-x)^m

* Better mapped operator construction

* ∞-sub-views of Clenshaw

* Update Project.toml

* Update jacobi.jl

* work on bugs

* Update clenshaw.jl

* fixes for lanczos tests

* Update test_stieltjes.jl

* Update Project.toml

* increase coverage

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <solver@mac.com>"]
-version = "0.4.3"
+version = "0.4.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -27,9 +27,9 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 [compat]
 ArrayLayouts = "0.7"
 BandedMatrices = "0.16.5"
-BlockArrays = "0.15, 0.16"
+BlockArrays = "0.16"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.8.5, 0.9"
+ContinuumArrays = "0.9"
 DomainSets = "0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
@@ -39,7 +39,7 @@ HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.11.1"
 InfiniteLinearAlgebra = "0.5.8"
 IntervalSets = "0.5"
-LazyArrays = "0.21.9"
+LazyArrays = "0.21.14"
 LazyBandedMatrices = "0.6"
 QuasiArrays = "0.7"
 SpecialFunctions = "1.0"

src/ClassicalOrthogonalPolynomials.jl

@@ -1,4 +1,6 @@
 module ClassicalOrthogonalPolynomials
+using LazyBandedMatrices: LazyBandedLayout
+using InfiniteArrays: parentindices
 using IntervalSets: UnitRange
 using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, BlockArrays,
     IntervalSets, DomainSets, ArrayLayouts, SpecialFunctions,
@@ -9,7 +11,7 @@ import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
                 IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
                 first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum, cumsum,
                 to_indices, _maybetail, tail, getproperty, inv, show, isapprox, summary
-import Base.Broadcast: materialize, BroadcastStyle, broadcasted
+import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
 import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout,
                 sub_materialize, arguments, sub_paddeddata, paddeddata, PaddedLayout, resizedata!, LazyVector, ApplyLayout, call,
                 _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout,
@@ -31,7 +33,7 @@ import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, Infini
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     inbounds_getindex, grid, plotgrid, transform, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, Expansion,
     AffineQuasiVector, AffineMap, WeightLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
-    checkpoints, weight, unweightedbasis, MappedBasisLayouts, __sum
+    checkpoints, weight, unweightedbasis, MappedBasisLayouts, __sum, invmap
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints
 
@@ -45,7 +47,7 @@ export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomia
             HermiteWeight, JacobiWeight, ChebyshevWeight, ChebyshevGrid, ChebyshevTWeight, ChebyshevUWeight, UltrasphericalWeight, LegendreWeight, LaguerreWeight,
             WeightedUltraspherical, WeightedChebyshev, WeightedChebyshevT, WeightedChebyshevU, WeightedJacobi,
             ∞, Derivative, .., Inclusion,
-            chebyshevt, chebyshevu, legendre, jacobi,
+            chebyshevt, chebyshevu, legendre, jacobi, ultraspherical,
             legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
             jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted
 
@@ -241,6 +243,26 @@ function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::Sub
     P[kr, :] * X
 end
 
+
+# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
+#     T = promote_type(eltype(f), eltype(C))
+#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
+#     P = parent(C)
+#     kr,jr = parentindices(C)
+#     (f[invmap(kr)] .* P)[kr,jr]
+# end
+
+# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
+#     T = promote_type(eltype(f), eltype(C))
+#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
+#     P = parent(C)
+#     kr,jr = parentindices(C)
+#     (f[invmap(kr)] .* P)[kr,jr]
+# end
+
+broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::Broadcasted, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) =
+    broadcast(*, materialize(f), C)
+
 function jacobimatrix(C::SubQuasiArray{T,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) where T
     P = parent(C)
     kr,jr = parentindices(C)

src/classical/jacobi.jl

@@ -62,9 +62,14 @@ singularities(a::AbstractAffineQuasiVector) = singularities(a.x)
 singularitiesbroadcast(_, L::LegendreWeight) = L # Assume we stay smooth
 singularitiesbroadcast(::typeof(exp), L::LegendreWeight) = L
 singularitiesbroadcast(::typeof(Base.literal_pow), ::typeof(^), L::LegendreWeight, ::Val) = L
+
+
+for op in (:+, :*)
+    @eval singularitiesbroadcast(::typeof($op), A, B, C, D...) = singularitiesbroadcast(*, singularitiesbroadcast(*, A, B), C, D...)
+end
 for op in (:+, :-, :*)
     @eval begin
-        singularitiesbroadcast(::typeof($op), ::LegendreWeight{T}, ::LegendreWeight{V}) where {T,V} = LegendreWeight{promote_type(T,V)}()
+        singularitiesbroadcast(::typeof($op), A::LegendreWeight, B::LegendreWeight) = LegendreWeight{promote_type(eltype(A), eltype(B))}()
         singularitiesbroadcast(::typeof($op), L::LegendreWeight, ::NoSingularities) = L
         singularitiesbroadcast(::typeof($op), ::NoSingularities, L::LegendreWeight) = L
     end
@@ -76,8 +81,10 @@ _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(::typeof(*), V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(*, map(_parent,V)...)[_parentindices(V...)...]
 
 
 singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::AbstractJacobiWeight) = b

src/classical/ultraspherical.jl

@@ -42,6 +42,7 @@ WeightedUltraspherical{T}(λ) where T = UltrasphericalWeight{T}(λ) .* Ultrasphe
 orthogonalityweight(C::Ultraspherical) = UltrasphericalWeight(C.λ)
 
 ultrasphericalc(n::Integer, λ, z::Number) = Base.unsafe_getindex(Ultraspherical{promote_type(typeof(λ),typeof(z))}(λ), z, n+1)
+ultraspherical(λ, d::AbstractInterval{T}) where T = Ultraspherical{float(promote_type(eltype(λ),T))}(λ)[affine(d,ChebyshevInterval{T}()), :]
 
 ==(a::Ultraspherical, b::Ultraspherical) = a.λ == b.λ
 ==(::Ultraspherical, ::ChebyshevT) = false

src/clenshaw.jl

@@ -31,14 +31,14 @@ for (get, vie) in ((:getindex, :view), (:(Base.unsafe_getindex), :(Base.unsafe_v
     end
 end
 
-function copyto!(dest::AbstractArray, V::SubArray{<:Any,1,<:OrthogonalPolynomial,<:Tuple{<:Number,<:OneTo}})
+function copyto!(dest::AbstractVector, V::SubArray{<:Any,1,<:OrthogonalPolynomial,<:Tuple{<:Number,<:OneTo}})
     P = parent(V)
     x,n = parentindices(V)
     A,B,C = recurrencecoefficients(P)
     forwardrecurrence!(dest, A, B, C, x, _p0(P))
 end
 
-function copyto!(dest::AbstractArray, V::SubArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{<:AbstractVector,<:AbstractUnitRange}})
+function forwardrecurrence_copyto!(dest::AbstractMatrix, V)
     checkbounds(dest, axes(V)...)
     P = parent(V)
     xr,jr = parentindices(V)
@@ -51,8 +51,10 @@ function copyto!(dest::AbstractArray, V::SubArray{<:Any,2,<:OrthogonalPolynomial
     end
     dest
 end
+copyto!(dest::AbstractMatrix, V::SubArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{<:AbstractVector,<:AbstractUnitRange}}) = forwardrecurrence_copyto!(dest, V)
+copyto!(dest::LayoutMatrix, V::SubArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{<:AbstractVector,<:AbstractUnitRange}}) = forwardrecurrence_copyto!(dest, V)
 
-function copyto!(dest::AbstractArray, V::SubArray{<:Any,1,<:OrthogonalPolynomial,<:Tuple{<:Number,<:UnitRange}})
+function copyto!(dest::AbstractVector, V::SubArray{<:Any,1,<:OrthogonalPolynomial,<:Tuple{<:Number,<:UnitRange}})
     checkbounds(dest, axes(V)...)
     P = parent(V)
     x,jr = parentindices(V)
@@ -267,6 +269,9 @@ Base.array_summary(io::IO, C::Clenshaw{T}, inds::Tuple{Vararg{OneToInf{Int}}}) w
 struct ClenshawLayout <: AbstractLazyBandedLayout end
 MemoryLayout(::Type{<:Clenshaw}) = ClenshawLayout()
 sublayout(::ClenshawLayout, ::Type{<:NTuple{2,AbstractUnitRange{Int}}}) = ClenshawLayout()
+sublayout(::ClenshawLayout, ::Type{<:Tuple{AbstractUnitRange{Int},Union{Slice,AbstractInfUnitRange{Int}}}}) = LazyBandedLayout()
+sublayout(::ClenshawLayout, ::Type{<:Tuple{Union{Slice,AbstractInfUnitRange{Int}},AbstractUnitRange{Int}}}) = LazyBandedLayout()
+sublayout(::ClenshawLayout, ::Type{<:Tuple{Union{Slice,AbstractInfUnitRange{Int}},Union{Slice,AbstractInfUnitRange{Int}}}}) = LazyBandedLayout()
 sub_materialize(::ClenshawLayout, V) = BandedMatrix(V)
 
 function _BandedMatrix(::ClenshawLayout, V::SubArray{<:Any,2})

test/test_chebyshev.jl

@@ -36,6 +36,8 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
             @testset "inf-range-indexing" begin
                 @test T[[begin,end],2:∞][:,2:5] == T[[-1,1],3:6]
             end
+
+            @test copyto!(BandedMatrix{Float64}(undef,(2,5),(1,4)),view(T,[0.1,0.2],1:5)) == T[[0.1,0.2],1:5]
         end
 
         @testset "U" begin
@@ -124,6 +126,10 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
             @test T == T[x,:] == T[x,1:∞] == T[:,1:∞]
             @test T[x,:] == T
         end
+
+        @testset "broadcast" begin
+            @test (x.^2 .* T)[0.1,1:10] ≈ 0.1^2 * T[0.1,1:10]
+        end
     end
 
     @testset "weighted" begin
@@ -164,8 +170,8 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
             @test v[0.1] ≈ let x = 0.1; exp(x)/(sqrt(x)*sqrt(1-x)) end
 
             WT̃ = w[2x .- 1] .* T[2x .- 1, :]
-            @test MemoryLayout(wT̃) isa MappedWeightedBasisLayout
-            v = wT̃ * (wT̃ \ @.(exp(x)/(sqrt(x)*sqrt(1-x))))
+            @test MemoryLayout(WT̃) isa MappedWeightedBasisLayout
+            v = WT̃ * (WT̃ \ @.(exp(x)/(sqrt(x)*sqrt(1-x))))
             @test v[0.1] ≈ let x = 0.1; exp(x)/(sqrt(x)*sqrt(1-x)) end
         end
     end
@@ -291,6 +297,14 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
         @test M[2:10,2:10] == M[1:10,1:10][2:10,2:10]
         @test M[3:10,5:20] == M[1:10,1:20][3:10,5:20]
         @test M[100_000:101_000,100_000:101_000] == M[2:1002,2:1002]
+        @test M[1:10,2:∞][:,1:5] == M[1:10,2:6]
+        @test M[1:10,:][:,1:5] == M[1:10,1:5]
+        @test M[2:∞,1:5][1:5,:] == M[2:6,1:5]
+        @test M[:,1:5][1:5,:] == M[1:5,1:5]
+        @test M[2:∞,3:∞][1:5,1:5] == M[2:6,3:7]
+        @test M[:,3:∞][1:5,1:5] == M[1:5,3:7]
+        @test M[2:∞,:][1:5,1:5] == M[2:6,1:5]
+        @test M[:,:][1:5,1:5] == M[1:5,1:5]
 
         @test Eye(∞) * M isa Clenshaw
 

test/test_lanczos.jl

@@ -240,4 +240,11 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, ort
         P = Normalized(Legendre())
         Cϕ = Pϕ\P
     end
+
+    @testset "3-mul-singularity" begin
+        m = 1
+        x = Inclusion(0..1)
+        Q = LanczosPolynomial(@. x^m*(1-x)^m*(2-x)^m)
+        @test Q[0.1,:]'*(Q \exp.(x)) ≈ exp(0.1)
+    end
 end

test/test_stieltjes.jl

@@ -90,7 +90,7 @@ end
         wT2 = wT[y,:]
         wU2 = wU[y,:]
         @test (Ultraspherical(1)[y,:]\(H*wT2))[1:10,1:10] == diagm(1 => fill(-π,9))
-        @test_broken (Chebyshev()[y,:]\(H*wU2))[1:10,1:10] == diagm(-1 => fill(1.0π,9))
+        @test (Chebyshev()[y,:]\(H*wU2))[1:10,1:10] == diagm(-1 => fill(1.0π,9))
 
         @testset "Legendre" begin
             P = Legendre()

test/test_ultraspherical.jl

@@ -119,6 +119,19 @@ import LazyArrays: rowsupport, colsupport
         @test (JacobiWeight(0,0) .* Jacobi(0,0)) \ Ultraspherical(1/2) isa Diagonal
     end
 
+    @testset "D^2 * mapped" begin
+        T = chebyshevt(0..1)
+        C = ultraspherical(2,0..1)
+        r = axes(T,1)
+        D = Derivative(r)
+        D₂ = C \ (D^2 * T)
+        # r²D₂ = C \ (r.^2 .* (D^2 * T))
+
+        c = [randn(100); zeros(∞)]
+        @test C[0.1,:]'*(D₂ * c) ≈ 4*(Derivative(axes(ChebyshevT(),1))^2 * (ChebyshevT() * c))[2*0.1-1]
+        @test_broken C[0.1,:]'*(r²D₂ * c) ≈ 0.1^2 * C[0.1,:]'*(D₂ * c)
+    end
+  
     @testset "BigFloat" begin
         U = Ultraspherical{BigFloat}(1)
         T = ChebyshevT{BigFloat}()