Fixes for other packages (#82)

* Fixes for other packages

* More Weighted Ldivs

* Update normalized.jl

* fixes

* Update test_chebyshev.jl

* fix tests

* 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.5.0"
+version = "0.5.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/ClassicalOrthogonalPolynomials.jl

@@ -109,7 +109,7 @@ _equals(::AbstractWeightedBasisLayout, ::AbstractOPLayout, _, _) = false # Weigh
 _equals(::WeightedOPLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ)
 _equals(::WeightedOPLayout, ::WeightedBasisLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
 _equals(::WeightedBasisLayout, ::WeightedOPLayout, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
-_equals(::WeightedBasisLayout{<:AbstractOPLayout}, ::WeightedBasisLayout{<:AbstractOPLayout}, A, B) = A.f == B.f && all(A.args .== B.args)
+_equals(::WeightedBasisLayout{<:AbstractOPLayout}, ::WeightedBasisLayout{<:AbstractOPLayout}, wP, wQ) = unweighted(wP) == unweighted(wQ) && weight(wP) == weight(wQ)
 
 
 copy(L::Ldiv{MappedOPLayout,Lay}) where Lay<:MappedBasisLayouts = copy(Ldiv{MappedBasisLayout,Lay}(L.A,L.B))

src/adaptivetransform.jl

@@ -31,7 +31,8 @@ function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiVector
             return pad(similar(cfs,0), ax)
         end
 
-        if maximum(abs,@views(cfs[end-2:end])) < 10tol*maxabsc
+        m = length(cfs)
+        if maximum(abs,@views(cfs[max(1,m-2):m])) < 10tol*maxabsc
             n = length(cfs)
             c = padchop!(cfs, tol, ax)
             un = A * c

src/classical/chebyshev.jl

@@ -255,13 +255,13 @@ end
 # sum
 ####
 
-function _sum(A::WeightedBasis{T,<:ChebyshevUWeight,<:ChebyshevU}, dims) where T
-    w, U = A.args
+function _sum(::Weighted{T,<:ChebyshevU}, dims) where T
     @assert dims == 1
     Hcat(convert(T, π)/2, Zeros{T}(1,∞))
 end
 
-function _sum(A::WeightedBasis{T,<:ChebyshevWeight,<:Chebyshev}, dims) where T
+# Same normalization for T,V,W
+function _sum(::Weighted{T,<:Chebyshev}, dims) where T
     @assert dims == 1
     Hcat(convert(T, π), Zeros{T}(1,∞))
 end

src/classical/legendre.jl

@@ -44,6 +44,8 @@ Legendre() = Legendre{Float64}()
 
 weighted(P::Legendre) = P
 weighted(P::Normalized{<:Any,<:Legendre}) = P
+weighted(P::SubQuasiArray{<:Any,2,<:Legendre}) = P
+weighted(P::SubQuasiArray{<:Any,2,<:Normalized{<:Any,<:Legendre}}) = P
 
 legendre() = Legendre()
 legendre(d::AbstractInterval{T}) where T = Legendre{float(T)}()[affine(d,ChebyshevInterval{T}()), :]

src/lanczos.jl

@@ -233,6 +233,8 @@ MemoryLayout(::Type{<:LanczosPolynomial}) = LanczosLayout()
 arguments(::ApplyLayout{typeof(*)}, Q::LanczosPolynomial) = Q.P, LanczosConversion(Q.data)
 copy(L::Ldiv{LanczosLayout,Lay}) where Lay<:AbstractLazyLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
 copy(L::Ldiv{LanczosLayout,Lay}) where Lay<:BroadcastLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
+copy(L::Ldiv{LanczosLayout,BroadcastLayout{typeof(*)}}) = copy(Ldiv{ApplyLayout{typeof(*)},BroadcastLayout{typeof(*)}}(L.A,L.B))
+copy(L::Ldiv{LanczosLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{ApplyLayout{typeof(*)},BroadcastLayout{typeof(*)}}(L.A,L.B))
 copy(L::Ldiv{LanczosLayout,ApplyLayout{typeof(*)}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
 copy(L::Ldiv{LanczosLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
 copy(L::Ldiv{LanczosLayout,<:ExpansionLayout}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A, L.B))

src/normalized.jl

@@ -176,6 +176,14 @@ ContinuumArrays.unweighted(wP::AbstractWeighted) = wP.P
 # copy(L::Ldiv{WeightedOPLayout,ApplyLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{UnknownLayout,ApplyLayout{typeof(*)}}(L.A, L.B))
 
 copy(L::Ldiv{<:WeightedOPLayout{<:NormalizedOPLayout},Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
+copy(L::Ldiv{Lay,<:WeightedOPLayout{<:NormalizedOPLayout}}) where Lay<:AbstractBasisLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A,L.B))
+copy(L::Ldiv{<:WeightedOPLayout{<:NormalizedOPLayout},<:WeightedOPLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
+copy(L::Ldiv{<:WeightedBasisLayout{<:NormalizedOPLayout},Lay}) where Lay<:AbstractBasisLayout = copy(Ldiv{ApplyLayout{typeof(*)},Lay}(L.A,L.B))
+copy(L::Ldiv{Lay,<:WeightedBasisLayout{<:NormalizedOPLayout}}) where Lay<:AbstractBasisLayout = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A,L.B))
+copy(L::Ldiv{<:WeightedBasisLayout{<:NormalizedOPLayout},<:WeightedBasisLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
+copy(L::Ldiv{<:WeightedOPLayout{<:NormalizedOPLayout},<:WeightedBasisLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
+copy(L::Ldiv{<:WeightedBasisLayout{<:NormalizedOPLayout},<:WeightedOPLayout{<:NormalizedOPLayout}}) = copy(Ldiv{ApplyLayout{typeof(*)},ApplyLayout{typeof(*)}}(L.A,L.B))
+
 
 # function layout_broadcasted(::ExpansionLayout{WeightedOPLayout}, ::OPLayout, ::typeof(*), a, P)
 #     axes(a,1) == axes(P,1) || throw(DimensionMismatch())
@@ -226,8 +234,15 @@ copy(Q::Weighted) = Q
 
 weight(wP::Weighted) = orthogonalityweight(wP.P)
 
-MemoryLayout(::Type{<:Weighted{<:Any,PP}}) where PP = WeightedOPLayout{typeof(MemoryLayout(PP))}()
+_weightedmemorylayout(::PP) where PP<:AbstractOPLayout =  WeightedOPLayout{PP}()
+_weightedmemorylayout(::PP) where PP =  WeightedBasisLayout{PP}()
+MemoryLayout(::Type{<:Weighted{<:Any,PP}}) where PP = _weightedmemorylayout(MemoryLayout(PP))
 
+function arguments(::ApplyLayout{typeof(*)}, Q::BroadcastQuasiMatrix{<:Any,typeof(*),<:Tuple{Weight,Normalized}})
+    w,Q = Q.args
+    P,D = arguments(*,Q)
+    (w .* P),D
+end
 function arguments(::ApplyLayout{typeof(*)}, Q::Weighted{<:Any,<:Normalized})
     P,D = arguments(*, Q.P)
     Weighted(P),D

test/test_chebyshev.jl

@@ -209,6 +209,13 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map, WeightedBasisLayout
 
             WU2 = Weighted(chebyshevu(0..1))
             @test (Derivative(x) * WU2)[0.1,1:10] ≈ 2(Derivative(axes(T,1))*Weighted(ChebyshevU()))[2*0.1-1,1:10]
+
+            @testset "Derivative" begin
+                y = affine(0..1,-1..1)
+                D = Derivative(Inclusion(0..1))
+                @test isbanded(wT[y,:] \ D * wU[y,:])
+                @test isbanded((wU[y,:]' * D').args[1])
+            end
         end
     end
 

test/test_jacobi.jl

@@ -466,5 +466,28 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
             @test HalfWeighted{:b}(Jacobi(2,0)) == Jacobi(2,0)
             @test Jacobi(2,0) == HalfWeighted{:b}(Jacobi(2,0))
         end
+
+        @testset "Weighted(Normalized)" begin
+            w = JacobiWeight(1,1)
+            P = Jacobi(1,1)
+            Q = Normalized(P)
+
+            @test w .* Q == w .* Q
+
+            @test Weighted(Q) \ Weighted(P) isa Diagonal
+            @test Weighted(P) \ Weighted(Q) isa Diagonal
+            @test Weighted(Q) \ Weighted(Q) isa Diagonal
+
+            @test (w .* Q) \ (w .* P) isa Diagonal
+            @test (w .* P) \ (w .* Q) isa Diagonal
+            @test (w .* Q) \ (w .* Q) isa Diagonal
+
+            @test (w .* Q) \ Weighted(P) isa Diagonal
+            @test (w .* P) \ Weighted(Q) isa Diagonal
+            @test (w .* Q) \ Weighted(Q) isa Diagonal
+            @test Weighted(P) \ (w .* Q) isa Diagonal
+            @test Weighted(Q) \ (w .* P) isa Diagonal
+            @test Weighted(Q) \ (w .* Q) isa Diagonal
+        end
     end
 end

test/test_lanczos.jl

@@ -90,6 +90,8 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, ort
         @test R \ [1; 2; 3; zeros(∞)] ≈ [1; 2; 3; zeros(∞)]
         @test (Q * (Q \ (1 .- x.^2)))[0.1] ≈ (1-0.1^2)
 
+        @test Q \ (x .* x) ≈ Q \ x.^2
+
         ũ = Normalized(P)*[1; 2; 3; zeros(∞)]
         u = Q*[1; 2; 3; zeros(∞)]
         ū = P * (P\u)
@@ -150,6 +152,8 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, ort
         Q = LanczosPolynomial(w, P)
         @test Q[0.1,1:10] ≈ Normalized(P)[0.1,1:10]
 
+        @test (Q \ (exp.(x) .* P)) * [1; zeros(∞)] ≈ Q \ exp.(x)
+
         w = @. exp(x) * sqrt(1-x)
         Q = LanczosPolynomial(w, P)
         # emperical from Julia