Better support for lowering conversions

.github/workflows/ci.yml

@@ -11,6 +11,7 @@ jobs:
       matrix:
         version:
           - '1.5'
+          - '^1.6.0-0'
         os:
           - ubuntu-latest
           - macOS-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -24,7 +24,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
 ArrayLayouts = "0.5.3"
-BandedMatrices = "0.16.2"
+BandedMatrices = "0.16.4"
 BlockArrays = "0.14.1"
 BlockBandedMatrices = "0.10"
 ContinuumArrays = "0.5"
@@ -36,8 +36,8 @@ FillArrays = "0.11"
 InfiniteArrays = "0.9"
 InfiniteLinearAlgebra = "0.4.6"
 IntervalSets = "0.3.1, 0.4, 0.5"
-LazyArrays = "0.20.4"
-QuasiArrays = "0.4.2"
+LazyArrays = "0.20.5"
+QuasiArrays = "0.4.5"
 SpecialFunctions = "0.10, 1"
 julia = "1.5"
 

src/ClassicalOrthogonalPolynomials.jl

@@ -40,9 +40,15 @@ export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomia
             WeightedUltraspherical, WeightedChebyshev, WeightedChebyshevT, WeightedChebyshevU, WeightedJacobi,
             ∞, Derivative, .., Inclusion, chebyshevt, chebyshevu, legendre, jacobi, legendrep, jacobip, ultrasphericalp, jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight
 
+if VERSION < v"1.6-"
+    oneto(n) = Base.OneTo(n)
+else
+    import Base: oneto
+end
+
 
 include("interlace.jl")
-    
+
 
 cardinality(::FullSpace{<:AbstractFloat}) = ℵ₁
 cardinality(::EuclideanDomain) = ℵ₁
@@ -75,7 +81,7 @@ function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiArray{
     tol = 20eps(real(T))
 
     for n = 2 .^ (4:∞)
-        An = A[:,OneTo(n)]
+        An = A[:,oneto(n)]
         cfs = An \ f
         maxabsc = maximum(abs, cfs)
         if maxabsc == 0 && maxabsfr == 0
@@ -121,9 +127,9 @@ satisfying for `x in axes(P,1)`
 P[x,2] == (A[1]*x + B[1])*P[x,1]
 P[x,n+1] == (A[n]*x + B[n])*P[x,n] - C[n]*P[x,n-1]
 ```
-Note that `C[1]`` is unused. 
+Note that `C[1]`` is unused.
 
-The relationship with the Jacobi matrix is: 
+The relationship with the Jacobi matrix is:
 ```julia
 1/A[n] == X[n+1,n]
 -B[n]/A[n] == X[n,n]
@@ -234,7 +240,7 @@ factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:OneT
 
 function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
     _,jr = parentindices(L)
-    ProjectionFactorization(factorize(parent(L)[:,Base.OneTo(maximum(jr))]), jr)
+    ProjectionFactorization(factorize(parent(L)[:,oneto(maximum(jr))]), jr)
 end
 
 function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial})

src/chebyshev.jl

@@ -173,6 +173,7 @@ function \(w_A::WeightedChebyshevT, w_B::WeightedChebyshevU)
     _BandedMatrix(Vcat(Fill(one(T)/2, 1, ∞), Zeros{T}(1, ∞), Fill(-one(T)/2, 1, ∞)), ∞, 2, 0)
 end
 
+\(T::ChebyshevT, U::ChebyshevU) = inv(U \ T)
 
 ####
 # interrelationships

src/clenshaw.jl

@@ -67,17 +67,17 @@ end
 getindex(P::OrthogonalPolynomial, x::Number, n::AbstractVector) = layout_getindex(P, x, n)
 getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractVector) = layout_getindex(P, x, n)
 getindex(P::SubArray{<:Any,1,<:OrthogonalPolynomial}, x::AbstractVector) = layout_getindex(P, x)
-getindex(P::OrthogonalPolynomial, x::Number, n::Number) = P[x,OneTo(n)][end]
+getindex(P::OrthogonalPolynomial, x::Number, n::Number) = P[x,oneto(n)][end]
 
 unsafe_layout_getindex(A...) = sub_materialize(Base.unsafe_view(A...))
 
 Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::AbstractUnitRange) = unsafe_layout_getindex(P, x, n)
 Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractUnitRange) = unsafe_layout_getindex(P, x, n)
-Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::AbstractVector) = Base.unsafe_getindex(P,x,OneTo(maximum(n)))[n]
-Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractVector) = Base.unsafe_getindex(P,x,OneTo(maximum(n)))[:,n]
+Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::AbstractVector) = Base.unsafe_getindex(P,x,oneto(maximum(n)))[n]
+Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractVector, n::AbstractVector) = Base.unsafe_getindex(P,x,oneto(maximum(n)))[:,n]
 Base.unsafe_getindex(P::OrthogonalPolynomial, x::AbstractVector, n::Number) = Base.unsafe_getindex(P, x, 1:n)[:,end]
 Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, ::Colon) = Base.unsafe_getindex(P, x, axes(P,2))
-Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::Number) = Base.unsafe_getindex(P,x,OneTo(n))[end]
+Base.unsafe_getindex(P::OrthogonalPolynomial, x::Number, n::Number) = Base.unsafe_getindex(P,x,oneto(n))[end]
 
 getindex(P::OrthogonalPolynomial, x::Number, jr::AbstractInfUnitRange{Int}) = view(P, x, jr)
 Base.unsafe_getindex(P::OrthogonalPolynomial{T}, x::Number, jr::AbstractInfUnitRange{Int}) where T = 

src/hermite.jl

@@ -12,7 +12,7 @@ struct Hermite{T} <: OrthogonalPolynomial{T} end
 Hermite() = Hermite{Float64}()
 
 ==(::Hermite, ::Hermite) = true
-axes(::Hermite{T}) where T = (Inclusion(ℝ), OneTo(∞))
+axes(::Hermite{T}) where T = (Inclusion(ℝ), oneto(∞))
 
 # H_{n+1} = 2x H_n - 2n H_{n-1}
 # 1/2 * H_{n+1} + n H_{n-1} = x H_n 

src/jacobi.jl

@@ -141,7 +141,7 @@ WeightedJacobi(a,b) = JacobiWeight(a,b) .* Jacobi(a,b)
 WeightedJacobi{T}(a,b) where T = JacobiWeight{T}(a,b) .* Jacobi{T}(a,b)
 
 
-axes(::AbstractJacobi{T}) where T = (Inclusion{T}(ChebyshevInterval{real(T)}()), OneTo(∞))
+axes(::AbstractJacobi{T}) where T = (Inclusion{T}(ChebyshevInterval{real(T)}()), oneto(∞))
 ==(P::Jacobi, Q::Jacobi) = P.a == Q.a && P.b == Q.b
 ==(P::Legendre, Q::Jacobi) = Jacobi(P) == Q
 ==(P::Jacobi, Q::Legendre) = P == Jacobi(Q)
@@ -276,6 +276,8 @@ function \(A::Jacobi, B::Jacobi)
     elseif A.b ≥ b+1
         J = Jacobi(a,b+1)
         (A \ J) * (J \ B)
+    elseif isinteger(A.a-a) && isinteger(A.b-b)
+        inv(B \ A)
     else
         error("not implemented for $A and $B")
     end

src/lanczos.jl

@@ -82,7 +82,12 @@ getindex(R::LanczosConversion, k::AbstractUnitRange, j::AbstractUnitRange) = _la
 
 inv(R::LanczosConversion) = ApplyArray(inv, R)
 
+
+Base.BroadcastStyle(::Type{<:LanczosConversion}) = LazyArrays.LazyArrayStyle{2}()
+
 struct LanczosConversionLayout <: AbstractLazyLayout end
+
+LazyArrays.simplifiable(::Mul{LanczosConversionLayout,<:PaddedLayout}) = Val(true)
 function copy(M::Mul{LanczosConversionLayout,<:PaddedLayout})
     resizedata!(M.A.data, maximum(colsupport(M.B)))
     M.A.data.R * M.B
@@ -107,7 +112,7 @@ function sub_paddeddata(::LanczosConversionLayout, S::SubArray{<:Any,1,<:Abstrac
     P = parent(S)
     (kr,j) = parentindices(S)
     resizedata!(P.data, j)
-    paddeddata(view(P.data.R, kr, j))
+    paddeddata(view(UpperTriangular(P.data.R.data.data), kr, j))
 end
 
 # struct LanczosJacobiMatrix{T,XX,WW} <: AbstractBandedMatrix{T}

src/ultraspherical.jl

@@ -141,6 +141,8 @@ function \(U::Ultraspherical{<:Any,<:Integer}, C::ChebyshevU)
     U\Ultraspherical(C)
 end
 
+\(T::Chebyshev, C::Ultraspherical) = inv(C \ T)
+
 function \(C2::Ultraspherical{<:Any,<:Integer}, C1::Ultraspherical{<:Any,<:Integer})
     λ = C1.λ
     T = promote_type(eltype(C2), eltype(C1))
@@ -162,6 +164,11 @@ function \(C2::Ultraspherical, C1::Ultraspherical)
         _BandedMatrix( Vcat(-(λ ./ ((0:∞) .+ λ))', Zeros(1,∞), (λ ./ ((0:∞) .+ λ))'), ∞, 0, 2)
     elseif C2.λ == λ
         Eye{T}(∞)
+    elseif isinteger(C2.λ-λ) && C2.λ > λ
+        Cm = Ultraspherical{T}(λ+1)
+        (C2 \ Cm) * (Cm \ C1)
+    elseif isinteger(C2.λ-λ)
+        inv(C1 \ C2)
     else
         error("Not implemented")
     end

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 Base: OneTo
+import ClassicalOrthogonalPolynomials: oneto
 import InfiniteLinearAlgebra: KronTrav, Block
 import FastTransforms: clenshaw!
 

test/test_chebyshev.jl

@@ -20,27 +20,27 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
     @testset "Evaluation" begin
         @testset "T" begin
             T = Chebyshev()
-            @test @inferred(T[0.1,Base.OneTo(0)]) == Float64[]
-            @test @inferred(T[0.1,Base.OneTo(1)]) == [1.0]
-            @test @inferred(T[0.1,Base.OneTo(2)]) == [1.0,0.1]
+            @test @inferred(T[0.1,oneto(0)]) == Float64[]
+            @test @inferred(T[0.1,oneto(1)]) == [1.0]
+            @test @inferred(T[0.1,oneto(2)]) == [1.0,0.1]
             for N = 1:10
-                @test @inferred(T[0.1,Base.OneTo(N)]) ≈ @inferred(T[0.1,1:N]) ≈ cos.((0:N-1)*acos(0.1))
+                @test @inferred(T[0.1,oneto(N)]) ≈ @inferred(T[0.1,1:N]) ≈ cos.((0:N-1)*acos(0.1))
                 @test @inferred(T[0.1,N]) ≈ cos((N-1)*acos(0.1))
             end
             @test T[0.1,[2,5,10]] ≈ [0.1,cos(4acos(0.1)),cos(9acos(0.1))]
 
-            @test axes(T[1:1,:]) === (Base.OneTo(1), Base.OneTo(∞))
+            @test axes(T[1:1,:]) === (oneto(1), oneto(∞))
             @test T[1:1,:][:,1:5] == ones(1,5)
             @test T[0.1,:][1:10] ≈ T[0.1,1:10] ≈ (T')[1:10,0.1]
         end
 
         @testset "U" begin
             U = ChebyshevU()
-            @test @inferred(U[0.1,Base.OneTo(0)]) == Float64[]
-            @test @inferred(U[0.1,Base.OneTo(1)]) == [1.0]
-            @test @inferred(U[0.1,Base.OneTo(2)]) == [1.0,0.2]
+            @test @inferred(U[0.1,oneto(0)]) == Float64[]
+            @test @inferred(U[0.1,oneto(1)]) == [1.0]
+            @test @inferred(U[0.1,oneto(2)]) == [1.0,0.2]
             for N = 1:10
-                @test @inferred(U[0.1,Base.OneTo(N)]) ≈ @inferred(U[0.1,1:N]) ≈ [sin((n+1)*acos(0.1))/sin(acos(0.1)) for n = 0:N-1]
+                @test @inferred(U[0.1,oneto(N)]) ≈ @inferred(U[0.1,1:N]) ≈ [sin((n+1)*acos(0.1))/sin(acos(0.1)) for n = 0:N-1]
                 @test @inferred(U[0.1,N]) ≈ sin(N*acos(0.1))/sin(acos(0.1))
             end
             @test U[0.1,[2,5,10]] ≈ [0.2,sin(5acos(0.1))/sin(acos(0.1)),sin(10acos(0.1))/sin(acos(0.1))]
@@ -59,7 +59,7 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
         @testset "ChebyshevT" begin
             T = Chebyshev()
             @test T == ChebyshevT() == chebyshevt()
-            Tn = @inferred(T[:,OneTo(100)])
+            Tn = @inferred(T[:,oneto(100)])
             @test grid(Tn) == chebyshevpoints(100, Val(1))
             P = factorize(Tn)
             u = T*[P.plan * exp.(grid(Tn)); zeros(∞)]
@@ -89,7 +89,7 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
             @test U == chebyshevu()
             @test U ≠ ChebyshevT()
             x = axes(U,1)
-            F = factorize(U[:,Base.OneTo(5)])
+            F = factorize(U[:,oneto(5)])
             @test @inferred(F \ x) ≈ [0,0.5,0,0,0]
             v = (x -> (3/20 + x + (2/5) * x^2)*exp(x)).(x)
             @inferred(U[:,Base.OneTo(5)]\v)
@@ -140,7 +140,7 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
     @testset "operators" begin
         T = ChebyshevT()
         U = ChebyshevU()
-        @test axes(T) == axes(U) == (Inclusion(ChebyshevInterval()),Base.OneTo(∞))
+        @test axes(T) == axes(U) == (Inclusion(ChebyshevInterval()),oneto(∞))
         D = Derivative(axes(T,1))
 
         @test T\T === pinv(T)*T === Eye(∞)
@@ -161,6 +161,10 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
         J = T\(x.*T)
         @test J isa BandedMatrix
         @test J[1:10,1:10] == jacobimatrix(T)[1:10,1:10]
+
+        @testset "inv" begin
+            @test (T \ U)[1:10,1:10] ≈ inv((U \ T)[1:10,1:10])
+        end
     end
 
     @testset "test on functions" begin
@@ -261,7 +265,11 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
 
     @testset "show" begin
         T = Chebyshev()
-        @test stringmime("text/plain", T * [1; 2; Zeros(∞)]) == "Chebyshev{1,Float64} * [1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, …]"
+        if VERSION < v"1.6-"
+            @test stringmime("text/plain", T * [1; 2; Zeros(∞)]) == "Chebyshev{1,Float64} * [1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, …]"
+        else
+            @test stringmime("text/plain", T * [1; 2; Zeros(∞)]) == "ChebyshevT{Float64} * [1.0, 2.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, …]"
+        end
     end
 
     @testset "Complex eltype" begin

test/test_hermite.jl

@@ -3,7 +3,7 @@ import ClassicalOrthogonalPolynomials: jacobimatrix
 
 @testset "Hermite" begin
     H = Hermite()
-    @test axes(H) == (Inclusion(ℝ), Base.OneTo(∞))
+    @test axes(H) == (Inclusion(ℝ), oneto(∞))
     x = axes(H,1)
     X = jacobimatrix(H)
 

test/test_jacobi.jl

@@ -105,6 +105,16 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, basis, MulQuasiMa
             @test A[1:10,1:10] == I
         end
 
+        @testset "Conversion" begin
+            A,B = Jacobi(0.25,-0.7), Jacobi(3.25,1.3)
+            R = B \ A
+            c = [[1,2,3,4,5]; zeros(∞)]
+            @test B[0.1,:]' * (R * c) ≈ A[0.1,:]' * c
+            Ri = A \ B
+            @test Ri[1:10,1:10] ≈ inv(R[1:10,1:10])
+            @test A[0.1,:]' * (Ri * c) ≈ B[0.1,:]' * c
+        end
+
         @testset "Derivative" begin
             a,b,c = 0.1,0.2,0.3
             S = Jacobi(a,b)

test/test_legendre.jl

@@ -12,7 +12,7 @@ import QuasiArrays: MulQuasiArray
 
     @testset "basics" begin
         P = Legendre()
-        @test axes(P) == (Inclusion(ChebyshevInterval()),Base.OneTo(∞))
+        @test axes(P) == (Inclusion(ChebyshevInterval()),oneto(∞))
         @test P == P == Legendre{Float32}()
         A,B,C = recurrencecoefficients(P)
         @test B isa Zeros

test/test_odes.jl

@@ -56,9 +56,9 @@ import SemiseparableMatrices: VcatAlmostBandedLayout
         B = BroadcastArray(+, Δ, (P\WS)'*(P'P)*(P\WS))
         @test colsupport(B,1) == 1:3
 
-        @test axes(B.args[2].args[1]) == (Base.OneTo(∞),Base.OneTo(∞))
-        @test axes(B.args[2]) == (Base.OneTo(∞),Base.OneTo(∞))
-        @test axes(B) == (Base.OneTo(∞),Base.OneTo(∞))
+        @test axes(B.args[2].args[1]) == (oneto(∞),oneto(∞))
+        @test axes(B.args[2]) == (oneto(∞),oneto(∞))
+        @test axes(B) == (oneto(∞),oneto(∞))
 
         @test BandedMatrix(view(B,1:10,13:20)) == zeros(10,8)
 

test/test_ultraspherical.jl

@@ -40,7 +40,7 @@ using ClassicalOrthogonalPolynomials, BandedMatrices, LazyArrays, Test
 
         @testset "Interrelationships" begin
             @testset "Chebyshev–Ultrashperical" begin
-                T = Chebyshev()
+                T = ChebyshevT()
                 U = ChebyshevU()
                 C = Ultraspherical(2)
                 D = Derivative(axes(T,1))
@@ -56,12 +56,28 @@ using ClassicalOrthogonalPolynomials, BandedMatrices, LazyArrays, Test
                 S₁ = (C\U)[1:10,1:10]
                 @test S₁ isa BandedMatrix{Float64}
                 @test S₁ == diagm(0 => 1 ./ (1:10), 2=> -(1 ./ (3:10)))
+
+                @test (U\C)[1:10,1:10] ≈ inv((C\U)[1:10,1:10])
+                @test (T\C)[1:10,1:10] ≈ inv((C\T)[1:10,1:10])
+                @test bandwidths(U\C) == bandwidths(T\C) == (0,∞)
+                @test colsupport(U\C,5) == colsupport(T\C,5) == 1:5
+                @test rowsupport(U\C,5) == rowsupport(T\C,5) == 5:∞
             end
             @testset "Legendre" begin
                 @test Ultraspherical(0.5) \ (UltrasphericalWeight(0.0) .* Ultraspherical(0.5)) == Eye(∞)
                 @test Legendre() \ (UltrasphericalWeight(0.0) .* Ultraspherical(0.5)) == Eye(∞)
+                @test (Legendre() \ Ultraspherical(1.5))[1:10,1:10] ≈ inv(Ultraspherical(1.5) \ Legendre())[1:10,1:10]
             end
         end
+
+        @testset "Conversion" begin
+            R = Ultraspherical(3.5) \ Ultraspherical(0.5)
+            c = [[1,2,3,4,5]; zeros(∞)]
+            @test Ultraspherical(3.5)[0.1,:]' * (R * c) ≈ Ultraspherical(0.5)[0.1,:]' * c
+            Ri = Ultraspherical(0.5) \ Ultraspherical(3.5)
+            @test Ri[1:10,1:10] ≈ inv(R[1:10,1:10])
+            @test Ultraspherical(0.5)[0.1,:]' * (Ri * c) ≈ Ultraspherical(3.5)[0.1,:]' * c
+        end
     end
 
     @testset "test on functions" begin