Add LanczosLayout (#32)

* Add LanczosLayout

* ambiguity

* increase coverage

* v0.3.6

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.3.5"
+version = "0.3.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -29,7 +29,7 @@ ArrayLayouts = "0.6.2"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.15"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.7"
+ContinuumArrays = "0.7.1"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"

src/ClassicalOrthogonalPolynomials.jl

@@ -30,7 +30,7 @@ import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, Infini
 import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
     inbounds_getindex, grid, transform, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, Expansion,
     AffineQuasiVector, AffineMap, WeightLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
-    checkpoints, weight, unweightedbasis
+    checkpoints, weight, unweightedbasis, MappedBasisLayouts
 import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                         _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints
 

src/classical/jacobi.jl

@@ -121,7 +121,7 @@ end
 Jacobi(a::V, b::T) where {T,V} = Jacobi{float(promote_type(T,V))}(a, b)
 
 jacobi(a,b) = Jacobi(a,b)
-jacobi(a,b, d::AbstractInterval{T}) where T = Jacobi(a,b)[affine(d,ChebyshevInterval{T}()), :]
+jacobi(a,b, d::AbstractInterval{T}) where T = Jacobi{float(promote_type(eltype(a),eltype(b),T))}(a,b)[affine(d,ChebyshevInterval{T}()), :]
 
 Jacobi(P::Legendre{T}) where T = Jacobi(zero(T), zero(T))
 

src/lanczos.jl

@@ -70,6 +70,8 @@ end
 size(::LanczosConversion) = (∞,∞)
 copy(R::LanczosConversion) = R
 
+Base.permutedims(R::LanczosConversion{<:Number}) = transpose(R)
+
 bandwidths(::LanczosConversion) = (0,∞)
 colsupport(L::LanczosConversion, j) = 1:maximum(j)
 rowsupport(L::LanczosConversion, j) = minimum(j):∞
@@ -178,6 +180,8 @@ end
 
 LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, orthonormalpolynomial(singularities(w)))
 
+
+
 orthogonalityweight(Q::LanczosPolynomial) = Q.w
 
 axes(Q::LanczosPolynomial) = (axes(Q.w,1),OneToInf())
@@ -221,8 +225,18 @@ function ldiv(Qn::SubQuasiArray{<:Any,2,<:LanczosPolynomial,<:Tuple{<:Inclusion,
     Q = parent(Qn)
     LanczosConversion(Q.data)[jr,jr] \ (Q.P[:,jr] \ C)
 end
+
+struct LanczosLayout <: AbstractBasisLayout end
+
+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,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))
 LazyArrays._mul_arguments(Q::LanczosPolynomial) = arguments(ApplyLayout{typeof(*)}(), Q)
 LazyArrays._mul_arguments(Q::QuasiAdjoint{<:Any,<:LanczosPolynomial}) = arguments(ApplyLayout{typeof(*)}(), Q)
 
-\(A::LanczosPolynomial, x::AbstractQuasiVector) = ApplyQuasiArray(A) \ x
+
+broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q::Union{Normalized,LanczosPolynomial}) = broadcastbasis(+, P.P, Q.P)
+broadcastbasis(::typeof(+), P::Union{Normalized,LanczosPolynomial}, Q) = broadcastbasis(+, P.P, Q)
+broadcastbasis(::typeof(+), P, Q::Union{Normalized,LanczosPolynomial}) = broadcastbasis(+, P, Q.P)

src/normalized.jl

@@ -123,6 +123,8 @@ for Lay in (:(ApplyLayout{typeof(*)}),:(BroadcastLayout{typeof(+)}),:(BroadcastL
     end
 end
 
+copy(L::Ldiv{Lay,<:NormalizedBasisLayout}) where Lay<:MappedBasisLayouts = copy(Ldiv{Lay,ApplyLayout{typeof(*)}}(L.A,L.B))
+
 # want to use special re-expansion routines without expanding Normalized basis
 @inline copy(L::Ldiv{<:NormalizedBasisLayout,BroadcastLayout{typeof(*)}}) = copy(Ldiv{BasisLayout,BroadcastLayout{typeof(*)}}(L.A, L.B))
 @inline copy(L::Ldiv{<:NormalizedBasisLayout,BroadcastLayout{typeof(*)},<:Any,<:AbstractQuasiVector}) = copy(Ldiv{BasisLayout,BroadcastLayout{typeof(*)}}(L.A, L.B))
@@ -196,4 +198,4 @@ broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, Q::Weighted) =
 @simplify *(Ac::QuasiAdjoint{<:Any,<:Weighted}, wB::Weighted) = 
     convert(WeightedOrthogonalPolynomial, parent(Ac))' * convert(WeightedOrthogonalPolynomial, wB)
 
-summary(io::IO, Q::Weighted) = print(io, "Weighted($(Q.P))")
+summary(io::IO, Q::Weighted) = print(io, "Weighted($(Q.P))")

test/test_jacobi.jl

@@ -163,6 +163,13 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
         f̃ = @.(sqrt(1 - x̃) * exp(x̃))
         @test wP̃[0.1,1:100]'*(wP̃[:,1:100] \ f̃) ≈ sqrt(1-0.1) * exp(0.1)
         @test (wP̃ * (wP̃ \ f̃))[0.1] ≈ sqrt(1-0.1) * exp(0.1)
+
+        @testset "bug" begin
+            P = jacobi(0,-1/2,0..1)
+            x = axes(P,1)
+            u = P * (P \ exp.(x))
+            @test u[0.1] ≈ exp(0.1)
+        end
     end
 
     @testset "trivial weight" begin

test/test_lanczos.jl

@@ -1,5 +1,5 @@
 using ClassicalOrthogonalPolynomials, BandedMatrices, ArrayLayouts, Test
-import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout
+import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, orthogonalityweight
 
 @testset "Lanczos" begin
     @testset "Legendre" begin
@@ -84,9 +84,19 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout
         w = P * [1; zeros(∞)];
         Q = LanczosPolynomial(w);
         R = Normalized(P) \ Q
+        @test bandwidths(R) == (0,∞)
+        @test orthogonalityweight(Q) == w
+        @test permutedims(R) === transpose(R)
         @test R * [1; 2; 3; zeros(∞)] ≈ [R[1:3,1:3] * [1,2,3]; zeros(∞)]
         @test R \ [1; 2; 3; zeros(∞)] ≈ [1; 2; 3; zeros(∞)]
         @test (Q * (Q \ (1 .- x.^2)))[0.1] ≈ (1-0.1^2)
+
+        ũ = Normalized(P)*[1; 2; 3; zeros(∞)]
+        u = Q*[1; 2; 3; zeros(∞)]
+        ū = P * (P\u)
+        @test (u + u)[0.1] ≈ (ũ + u)[0.1] ≈ (u + ũ)[0.1] ≈ (ũ + ũ)[0.1] ≈ (ū + u)[0.1] ≈ (u + ū)[0.1] ≈ (ū + ū)[0.1] ≈ 2u[0.1]
+
+        @test Q \ u ≈ Q \ ũ ≈ Q \ ū
     end
 
     @testset "Jacobi via Lanczos" begin