Remove templating of types for cached vectors (#19)

* Remove templating of types for cached vectors

* Update Project.toml

* Add Laguerre

* Update test_laguerre.jl

* Update Project.toml

* Update Project.toml

* update op name

* Update test_stieltjes.jl

* Update test_lanczos.jl

* Laguerre tests

* increase coverage

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.3"
+version = "0.3.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -26,19 +26,19 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 [compat]
 ArrayLayouts = "0.6.2"
 BandedMatrices = "0.16.5"
-BlockArrays = "0.14.1"
+BlockArrays = "0.14.1, 0.15"
 BlockBandedMatrices = "0.10"
 ContinuumArrays = "0.6"
 DomainSets = "0.4, 0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
 FastTransforms = "0.11, 0.12"
 FillArrays = "0.11"
-InfiniteArrays = "0.10"
-InfiniteLinearAlgebra = "0.5.2"
+InfiniteArrays = "0.10.3"
+InfiniteLinearAlgebra = "0.5.3"
 IntervalSets = "0.5"
-LazyArrays = "0.20.8"
-LazyBandedMatrices = "0.5"
+LazyArrays = "0.21"
+LazyBandedMatrices = "0.5.1"
 QuasiArrays = "0.4.6"
 SpecialFunctions = "0.10, 1"
 julia = "1.5"

README.md

@@ -18,7 +18,7 @@ julia> chebyshevt.(0:5,0.1)
   0.9208
   0.48016
 ```
-Other examples include `chebyshevu`, `legendrep`, `jacobip`, and `ultrasphericalp`.
+Other examples include `chebyshevu`, `legendrep`, `jacobip`, `ultrasphericalc`, `hermiteh` and `laguerrel`.
 
 For expansion, it supports usafe as quasi-arrays where one one axes is continuous and the other axis is discrete (countably infinite), as implemented in [QuasiArrays.jl](https://github.com/JuliaApproximation/QuasiArrays.jl) and  [ContinuumArrays.jl](https://github.com/JuliaApproximation/ContinuumArrays.jl).  
 ```julia

src/ClassicalOrthogonalPolynomials.jl

@@ -39,10 +39,14 @@ import FastGaussQuadrature: jacobimoment
 import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray
 import BandedMatrices: bandwidths
 
-export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial, Hermite, Jacobi, Legendre, Chebyshev, ChebyshevT, ChebyshevU, ChebyshevInterval, Ultraspherical, Fourier,
-            HermiteWeight, JacobiWeight, ChebyshevWeight, ChebyshevGrid, ChebyshevTWeight, ChebyshevUWeight, UltrasphericalWeight, LegendreWeight,
+export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial, 
+            Hermite, Jacobi, Legendre, Chebyshev, ChebyshevT, ChebyshevU, ChebyshevInterval, Ultraspherical, Fourier, Laguerre,
+            HermiteWeight, JacobiWeight, ChebyshevWeight, ChebyshevGrid, ChebyshevTWeight, ChebyshevUWeight, UltrasphericalWeight, LegendreWeight, LaguerreWeight,
             WeightedUltraspherical, WeightedChebyshev, WeightedChebyshevT, WeightedChebyshevU, WeightedJacobi,
-            ∞, Derivative, .., Inclusion, chebyshevt, chebyshevu, legendre, jacobi, legendrep, jacobip, ultrasphericalp, jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight
+            ∞, Derivative, .., Inclusion, 
+            chebyshevt, chebyshevu, legendre, jacobi,
+            legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh,
+            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight
 
 if VERSION < v"1.6-"
     oneto(n) = Base.OneTo(n)
@@ -277,12 +281,14 @@ function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
     A\B
 end
 
+include("ratios.jl")
 include("normalized.jl")
 include("lanczos.jl")
 include("hermite.jl")
 include("jacobi.jl")
 include("chebyshev.jl")
 include("ultraspherical.jl")
+include("laguerre.jl")
 include("fourier.jl")
 include("stieltjes.jl")
 

src/hermite.jl

@@ -1,3 +1,8 @@
+"""
+   HermiteWeight()
+
+is a quasi-vector representing `exp(-x^2)` on ℝ.
+"""
 struct HermiteWeight{T} <: Weight{T} end
 
 HermiteWeight() = HermiteWeight{Float64}()
@@ -16,6 +21,14 @@ orthogonalityweight(::Hermite{T}) where T = HermiteWeight{T}()
 ==(::Hermite, ::Hermite) = true
 axes(::Hermite{T}) where T = (Inclusion(ℝ), oneto(∞))
 
+"""
+     hermiteh(n, z)
+
+computes the `n`-th Hermite polynomial, orthogonal with 
+respec to `exp(-x^2)`, at `z`.
+"""
+hermiteh(n::Integer, z::Number) = Base.unsafe_getindex(Hermite{typeof(z)}(), z, n+1)
+
 # H_{n+1} = 2x H_n - 2n H_{n-1}
 # 1/2 * H_{n+1} + n H_{n-1} = x H_n 
 # x*[H_0 H_1 H_2 …] = [H_0 H_1 H_2 …] * [0    1; 1/2  0     2; 1/2   0  3; …]   

src/laguerre.jl

@@ -0,0 +1,67 @@
+"""
+   LaguerreWeight(α)
+
+is a quasi-vector representing `x^α * exp(-x)` on `0..Inf`.
+"""
+struct LaguerreWeight{T} <: Weight{T}
+    α::T
+end
+
+LaguerreWeight{T}() where T = LaguerreWeight{T}(zero(T))
+LaguerreWeight() = LaguerreWeight{Float64}()
+axes(::LaguerreWeight{T}) where T = (Inclusion(ℝ),)
+function getindex(w::LaguerreWeight, x::Number)
+    x ∈ axes(w,1) || throw(BoundsError())
+    x^w.α * exp(-x)
+end
+
+sum(L::LaguerreWeight{T}) where T = gamma(L.α + 1)
+
+struct Laguerre{T} <: OrthogonalPolynomial{T} 
+    α::T
+    Laguerre{T}(α) where T = new{T}(convert(T, α))
+end
+Laguerre{T}() where T = Laguerre{T}(zero(T))
+Laguerre() = Laguerre{Float64}()
+Laguerre(α::T) where T = Laguerre{float(T)}(α)
+orthogonalityweight(L::Laguerre)= LaguerreWeight(L.α)
+
+==(L1::Laguerre, L2::Laguerre) = L1.α == L2.α
+axes(::Laguerre{T}) where T = (Inclusion(HalfLine{T}()), oneto(∞))
+
+"""
+     laguerrel(n, α, z)
+
+computes the `n`-th generalized Laguerre polynomial, orthogonal with 
+respec to `x^α * exp(-x)`, at `z`.
+"""
+laguerrel(n::Integer, α, z::Number) = Base.unsafe_getindex(Laguerre{promote_type(typeof(α), typeof(z))}(α), z, n+1)
+
+"""
+     laguerrel(n, z)
+
+computes the `n`-th Laguerre polynomial, orthogonal with 
+respec to `exp(-x)`, at `z`.
+"""
+laguerrel(n::Integer, z::Number) = laguerrel(n, 0, z)
+
+
+# L_{n+1} = (-1/(n+1) x + (2n+α+1)/(n+1)) L_n - (n+α)/(n+1) L_{n-1}
+# - (n+α) L_{n-1} + (2n+α+1)* L_n -(n+1) L_{n+1} = x  L_n
+# x*[L_0 L_1 L_2 …] = [L_0 L_1 L_2 …] * [(α+1)    -(α+1); -1  (α+3)     -(α+2);0  -2   (α+5) -(α+3); …]   
+function jacobimatrix(L::Laguerre{T}) where T
+    α = L.α
+    Tridiagonal(-(1:∞), (α+1):2:∞, -(α+1:∞))
+end
+
+recurrencecoefficients(L::Laguerre{T}) where T = ((-one(T)) ./ (1:∞), ((L.α+1):2:∞) ./ (1:∞), (L.α:∞) ./ (1:∞))
+
+##########
+# Derivatives
+##########
+
+@simplify function *(D::Derivative, L::Laguerre)
+    T = promote_type(eltype(D),eltype(L))
+    D = _BandedMatrix(Fill(-one(T),1,∞), ∞, -1,1)
+    Laguerre(L.α+1)*D
+end

src/lanczos.jl

@@ -25,24 +25,24 @@ end
 const PaddedVector{T} = CachedVector{T,Vector{T},Zeros{T,1,Tuple{OneToInf{Int}}}}
 const PaddedMatrix{T} = CachedMatrix{T,Matrix{T},Zeros{T,2,NTuple{2,OneToInf{Int}}}}
 
-mutable struct LanczosData{T,XX,WW}
-    X::XX
-    W::WW
+mutable struct LanczosData{T}
+    X
+    W
     γ::PaddedVector{T}
     β::PaddedVector{T}
     R::UpperTriangular{T,PaddedMatrix{T}}
     ncols::Int
 
-    function LanczosData{T,XX,WW}(X, W, γ, β, R) where {T,XX,WW}
+    function LanczosData{T}(X, W, γ, β, R) where {T}
         R[1,1] = 1;
         p0 = view(R,:,1);
         γ[1] = sqrt(dot(p0,W,p0))
         lmul!(inv(γ[1]), p0)
-        new{T,XX,WW}(X, W, γ, β, R, 1)
+        new{T}(X, W, γ, β, R, 1)
     end
 end
 
-LanczosData(X::XX, W::WW, γ::AbstractVector{T}, β, R) where {T,XX,WW} = LanczosData{T,XX,WW}(X, W, γ, β, R)
+LanczosData(X, W, γ::AbstractVector{T}, β, R) where {T} = LanczosData{T}(X, W, γ, β, R)
 LanczosData(X::AbstractMatrix{T}, W::AbstractMatrix{T}) where T = LanczosData(X, W, zeros(T,∞), zeros(T,∞), UpperTriangular(zeros(T,∞,∞)))
 
 function LanczosData(w::AbstractQuasiVector, P::AbstractQuasiMatrix)
@@ -63,8 +63,8 @@ function resizedata!(L::LanczosData, N)
     L
 end
 
-struct LanczosConversion{T,XX,WW} <: LayoutMatrix{T}
-    data::LanczosData{T,XX,WW}
+struct LanczosConversion{T} <: LayoutMatrix{T}
+    data::LanczosData{T}
 end
 
 size(::LanczosConversion) = (∞,∞)
@@ -115,12 +115,12 @@ function sub_paddeddata(::LanczosConversionLayout, S::SubArray{<:Any,1,<:Abstrac
     paddeddata(view(UpperTriangular(P.data.R.data.data), kr, j))
 end
 
-# struct LanczosJacobiMatrix{T,XX,WW} <: AbstractBandedMatrix{T}
-#     data::LanczosData{T,XX,WW}
+# struct LanczosJacobiMatrix{T} <: AbstractBandedMatrix{T}
+#     data::LanczosData{T}
 # end
 
-struct LanczosJacobiBand{T,XX,WW} <: LazyVector{T}
-    data::LanczosData{T,XX,WW}
+struct LanczosJacobiBand{T} <: LazyVector{T}
+    data::LanczosData{T}
     diag::Symbol
 end
 
@@ -145,46 +145,11 @@ getindex(K::SubArray{<:Any,1,<:LanczosJacobiBand}, k::AbstractInfUnitRange{<:Int
 
 copy(A::LanczosJacobiBand) = A # immutable
 
-struct LanczosRecurrence{ABC,T,XX,WW} <: LazyVector{T}
-    data::LanczosData{T,XX,WW}
-end
-
-LanczosRecurrence{ABC}(data::LanczosData{T,XX,WW}) where {ABC,T,XX,WW} = LanczosRecurrence{ABC,T,XX,WW}(data)
-
-size(P::LanczosRecurrence) = (∞,)
-
-resizedata!(A::LanczosRecurrence, n) = resizedata!(A.data, n)
-
-
-
-# Q[x,n] = (x/γ[n] - β[n-1]/γ[n])  * Q[x,n-1] - γ[n-1]/γ[n] * Q[x,n]
-
-
-function _lanczos_getindex(A::LanczosRecurrence{:A}, I)
-    resizedata!(A, maximum(I)+1)
-    inv.(A.data.γ.data[I .+ 1])
-end
-
-function _lanczos_getindex(B::LanczosRecurrence{:B}, I)
-    resizedata!(B, maximum(I)+1)
-    B.data.β.data[I] ./ B.data.γ.data[I .+ 1]
-end
-
-function _lanczos_getindex(C::LanczosRecurrence{:C}, I)
-    resizedata!(C, maximum(I)+1)
-    C.data.γ.data[I] ./ C.data.γ.data[I  .+ 1]
-end
-
-getindex(A::LanczosRecurrence, I::Integer) = _lanczos_getindex(A, I)
-getindex(A::LanczosRecurrence, I::AbstractVector) = _lanczos_getindex(A, I)
-getindex(K::LanczosRecurrence, k::InfRanges{<:Integer}) = view(K, k)
-getindex(K::SubArray{<:Any,1,<:LanczosRecurrence}, k::InfRanges{<:Integer}) = view(K, k)
-
 
-struct LanczosPolynomial{T,XX,WW,Weight,Basis} <: OrthogonalPolynomial{T}
+struct LanczosPolynomial{T,Weight,Basis} <: OrthogonalPolynomial{T}
     w::Weight # Weight of orthogonality
     P::Basis # Basis we use to represent the OPs
-    data::LanczosData{T,XX,WW}
+    data::LanczosData{T}
 end
 
 ==(A::LanczosPolynomial, B::LanczosPolynomial) = A.w == B.w
@@ -217,8 +182,8 @@ axes(Q::LanczosPolynomial) = (axes(Q.w,1),OneToInf())
 
 _p0(Q::LanczosPolynomial) = inv(Q.data.γ[1])*_p0(Q.P)
 
-recurrencecoefficients(Q::LanczosPolynomial) = LanczosRecurrence{:A}(Q.data),LanczosRecurrence{:B}(Q.data),LanczosRecurrence{:C}(Q.data)
 jacobimatrix(Q::LanczosPolynomial) = SymTridiagonal(LanczosJacobiBand(Q.data, :d), LanczosJacobiBand(Q.data, :du))
+recurrencecoefficients(Q::LanczosPolynomial) = normalized_recurrencecoefficients(Q)
 
 Base.summary(io::IO, Q::LanczosPolynomial{T}) where T = print(io, "LanczosPolynomial{$T} with weight with singularities $(singularities(Q.w))")
 Base.show(io::IO, Q::LanczosPolynomial{T}) where T = summary(io, Q)

src/normalized.jl

@@ -1,50 +1,26 @@
 """
-   NormalizationConstant
+   normalizationconstant
 
 gives the normalization constants so that the jacobi matrix is symmetric,
 that is, so we have orthonormal OPs:
 
-    Q == P*NormalizationConstant(P)
+    Q == P*normalizationconstant(P)
 """
-
-
-mutable struct NormalizationConstant{T, PP<:AbstractQuasiMatrix{T}} <: AbstractCachedVector{T}
-    P::PP # OPs
-    data::Vector{T}
-    datasize::Tuple{Int}
-
-    NormalizationConstant{T, PP}(μ, P::PP) where {T,PP<:AbstractQuasiMatrix{T}} = new{T, PP}(P, T[μ], (1,))
-end
-
-NormalizationConstant(μ, P::AbstractQuasiMatrix{T}) where T = NormalizationConstant{T,typeof(P)}(μ, P)
-NormalizationConstant(P::AbstractQuasiMatrix) = NormalizationConstant(inv(sqrt(sum(orthogonalityweight(P)))), P)
-
-size(K::NormalizationConstant) = (ℵ₀,)
-
-# How we populate the data
-# function _normalizationconstant_fill_data!(K::NormalizationConstant, J::Union{BandedMatrix,Symmetric{<:Any,BandedMatrix},Tridiagonal,SymTridiagonal}, inds)
-#     dl, _, du = bands(J)
-#     @inbounds for k in inds
-#         K.data[k] = sqrt(du[k-1]/dl[k]) * K.data[k-1]
-#     end
-# end
-
-function _normalizationconstant_fill_data!(K::NormalizationConstant, J, inds)
-    @inbounds for k in inds
-        K.data[k] = sqrt(J[k,k-1]/J[k-1,k]) * K.data[k-1]
-    end
+function normalizationconstant(μ, P::AbstractQuasiMatrix{T}) where T
+    X = jacobimatrix(P)
+    c,b = subdiagonaldata(X),supdiagonaldata(X)
+    # hide array type to avoid crazy compilation
+    Accumulate{T,1,typeof(*),Vector{T},AbstractVector{T}}(*, T[μ], Vcat(zero(T),sqrt.(c ./ b)), 1, (1,))
 end
 
-
-LazyArrays.cache_filldata!(K::NormalizationConstant, inds) = _normalizationconstant_fill_data!(K, jacobimatrix(K.P), inds)
+normalizationconstant(P::AbstractQuasiMatrix) = normalizationconstant(inv(sqrt(sum(orthogonalityweight(P)))), P)
 
 
 struct Normalized{T, OPs<:AbstractQuasiMatrix{T}, NL} <: OrthogonalPolynomial{T}
     P::OPs
     scaling::NL # Q = P * Diagonal(scaling)
 end
 
-normalizationconstant(P) = NormalizationConstant(P)
 Normalized(P::AbstractQuasiMatrix{T}) where T = Normalized(P, normalizationconstant(P))
 Normalized(Q::Normalized) = Q
 normalized(P) = Normalized(P)
@@ -83,24 +59,26 @@ _p0(Q::Normalized) = Q.scaling[1]
 # q_{n+1}/h[n+1] = (A_n * x + B_n) * q_n/h[n] - C_n * p_{n-1}/h[n-1]
 # q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}
 
-function recurrencecoefficients(Q::Normalized)
-    X = jacobimatrix(Q.P)
-    c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
-    inv.(sqrt.(b .* c)), -(a ./ sqrt.(b .* c)), Vcat(zero(eltype(Q)), sqrt.(b .* c) ./ sqrt.(b[2:end] .* c[2:end]))
+function normalized_recurrencecoefficients(Q::AbstractQuasiMatrix{T}) where T
+    X = jacobimatrix(Q)
+    a,b = diagonaldata(X), supdiagonaldata(X)
+    inv.(b), -(a ./ b), Vcat(zero(T), b) ./ b
 end
 
+recurrencecoefficients(Q::Normalized{T}) where T = normalized_recurrencecoefficients(Q)
+
 # x * p[n] = c[n-1] * p[n-1] + a[n] * p[n] + b[n] * p[n+1]
 # x * q[n]/h[n] = c[n-1] * q[n-1]/h[n-1] + a[n] * q[n]/h[n] + b[n] * q[n+1]/h[n+1]
 # x * q[n+1] = c[n-1] * h[n]/h[n-1] * q[n-1] + a[n] * q[n] + b[n] * h[n]/h[n+1] * q[n+1]
 
 # q_{n+1}/h[n+1] = (A_n * x + B_n) * q_n/h[n] - C_n * p_{n-1}/h[n-1]
 # q_{n+1} = (h[n+1]/h[n] * A_n * x + h[n+1]/h[n] * B_n) * q_n - h[n+1]/h[n-1] * C_n * p_{n-1}
 
-function symtridagonalize(X)
+function symtridiagonalize(X)
     c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
     SymTridiagonal(a, sqrt.(b .* c))
 end
-jacobimatrix(Q::Normalized) = symtridagonalize(jacobimatrix(Q.P))
+jacobimatrix(Q::Normalized) = symtridiagonalize(jacobimatrix(Q.P))
 
 orthogonalityweight(Q::Normalized) = orthogonalityweight(Q.P)
 singularities(Q::Normalized) = singularities(Q.P)
@@ -163,14 +141,13 @@ end
 ###
 # show
 ###
-Base.array_summary(io::IO, C::NormalizationConstant{T}, inds) where T = print(io, "NormalizationConstant{$T}")
 show(io::IO, Q::Normalized) = print(io, "Normalized($(Q.P))")
 show(io::IO, ::MIME"text/plain", Q::Normalized) = show(io, Q)
 
 
 
 struct OrthonormalWeighted{T, PP<:AbstractQuasiMatrix{T}} <: Basis{T}
-    P::Normalized{T, PP, NormalizationConstant{T, PP}}
+    P::Normalized{T, PP}
 end
 
 function OrthonormalWeighted(P)