Start PiecewiseInterlace (#26)

* Start PiecewiseInterlace

* copy for LanczosConversion

* tests pass

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -7,7 +7,7 @@ using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, Bloc
     InfiniteLinearAlgebra, InfiniteArrays, LinearAlgebra, FastGaussQuadrature, FastTransforms, FFTW,
     LazyBandedMatrices, HypergeometricFunctions
 
-import Base: @_inline_meta, axes, getindex, convert, prod, *, /, \, +, -,
+import Base: @_inline_meta, axes, getindex, unsafe_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
@@ -18,7 +18,7 @@ import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupp
                 AbstractCachedVector, AbstractCachedMatrix, paddeddata, cache_filldata!
 import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
                         subdiagonaldata, diagonaldata, supdiagonaldata
-import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, AbstractLazyBandedLayout
+import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, unitblocks, BlockRange1, AbstractLazyBandedLayout
 import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot
 import BandedMatrices: AbstractBandedLayout, AbstractBandedMatrix, _BandedMatrix, bandeddata
 import FillArrays: AbstractFill, getindex_value, SquareEye
@@ -39,7 +39,7 @@ import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurr
 
 import FastGaussQuadrature: jacobimoment
 
-import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray
+import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray, block, blockindex, BlockSlice
 import BandedMatrices: bandwidths
 
 export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial,
@@ -49,7 +49,7 @@ export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomia
             ∞, Derivative, .., Inclusion,
             chebyshevt, chebyshevu, legendre, jacobi, ultraspherical,
             legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
-            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted
+            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted, PiecewiseInterlace
 
 
 import Base: oneto
@@ -61,7 +61,7 @@ include("interlace.jl")
 cardinality(::FullSpace{<:AbstractFloat}) = ℵ₁
 cardinality(::EuclideanDomain) = ℵ₁
 
-transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V}  = 
+transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V}  =
     adaptivetransform_ldiv(convert(AbstractQuasiArray{promote_type(T,V)}, A), f)
 
 function chop!(c::AbstractVector, tol::Real)

src/clenshaw.jl

@@ -95,11 +95,16 @@ end
 # Clenshaw
 ###
 
-function getindex(f::Expansion{<:Any,<:OrthogonalPolynomial}, x::Number)
+function unsafe_getindex(f::Expansion{<:Any,<:OrthogonalPolynomial}, x::Number)
     P,c = arguments(f)
     _p0(P)*clenshaw(paddeddata(c), recurrencecoefficients(P)..., x)
 end
 
+Base.@propagate_inbounds function getindex(f::Expansion{<:Any,<:OrthogonalPolynomial}, x::Number)
+    @inbounds checkbounds(f, x)
+    unsafe_getindex(f, x)
+end
+
 getindex(f::Expansion{T,<:OrthogonalPolynomial}, x::AbstractVector{<:Number}) where T = 
     copyto!(Vector{T}(undef, length(x)), view(f, x))
 

src/interlace.jl

@@ -93,3 +93,79 @@ function interlace!(v::AbstractVector,offset::Int)
     end
     v
 end
+
+
+abstract type AbstractInterlaceBasis{T} <: Basis{T} end
+copy(A::AbstractInterlaceBasis) = interlacebasis(A, map(copy, A.args))
+
+"""
+    PiecewiseInterlace(args...)
+
+is an analogue of `Basis` that takes the union of the first axis,
+and the second axis is a blocked interlace of args.
+If there is overlap, it uses the first in order.
+"""
+struct PiecewiseInterlace{T, Args} <: AbstractInterlaceBasis{T}
+    args::Args
+end
+
+PiecewiseInterlace{T}(args...) where T = PiecewiseInterlace{T,typeof(args)}(args)
+PiecewiseInterlace(args...) = PiecewiseInterlace{mapreduce(eltype,promote_type,args)}(args...)
+PiecewiseInterlace{T}(args::AbstractVector) where T = PiecewiseInterlace{T,typeof(args)}(args)
+PiecewiseInterlace(args::AbstractVector) = PiecewiseInterlace{eltype(eltype(args))}(args)
+
+interlacebasis(::PiecewiseInterlace, args...) = PiecewiseInterlace(args...)
+
+
+axes(A::PiecewiseInterlace) = (union(axes.(A.args,1)...), LazyBandedMatrices._block_vcat_axes(unitblocks.(axes.(A.args,2))...))
+
+==(A::PiecewiseInterlace, B::PiecewiseInterlace) = all(A.args .== B.args)
+
+
+function QuasiArrays._getindex(::Type{IND}, A::PiecewiseInterlace{T}, (x,j)::IND) where {IND,T}
+    Jj = findblockindex(axes(A,2), j)
+    @boundscheck x in axes(A,1) || throw(BoundsError(A, (x,j)))
+    J = Int(block(Jj))
+    i = blockindex(Jj)
+    x in axes(A.args[i],1) && return A.args[i][x, J]
+    zero(T)
+end
+
+function \(A::AbstractInterlaceBasis, B::AbstractInterlaceBasis)
+    axes(A,1) == axes(B,1) || throw(DimensionMismatch())
+    T = promote_type(eltype(A),eltype(B))
+    A == B && return Eye{T}((axes(A,2),))
+    BlockBroadcastArray{T}(Diagonal, unitblocks.((\).(A.args, B.args))...)
+end
+
+@simplify function *(D::Derivative, S::AbstractInterlaceBasis)
+    axes(D,2) == axes(S,1) || throw(DimensionMismatch())
+    args = arguments.(Ref(ApplyLayout{typeof(*)}()), Derivative.(axes.(S.args,1)) .* S.args)
+    all(length.(args) .== 2) || error("Not implemented")
+    interlacebasis(S, map(first, args)...) * BlockBroadcastArray{promote_type(eltype(D),eltype(S))}(Diagonal, unitblocks.(last.(args))...)
+end
+
+struct PiecewiseFactorization{T,FF,Ax} <: Factorization{T}
+    factorizations::FF
+    axes::Ax
+end
+
+PiecewiseFactorization{T}(fac, ax) where T = PiecewiseFactorization{T,typeof(fac),typeof(ax)}(fac, ax)
+
+function \(F::PiecewiseFactorization{T}, v::AbstractQuasiVector) where {T}
+    BlockBroadcastArray{T}(vcat, unitblocks.((\).(F.factorizations, getindex.(Ref(v), F.axes)))...)
+end
+
+function factorize(V::SubQuasiArray{T,2,<:AbstractInterlaceBasis,<:Tuple{Inclusion,BlockSlice{BlockRange1{OneTo{Int}}}}}) where T
+    P = parent(V)
+    _,jr = parentindices(V)
+    N = Int(last(jr.block))
+    PiecewiseFactorization{T}(factorize.(view.(P.args, :, Ref(Base.OneTo(N)))), axes.(P.args,1))
+end
+
+function factorize(V::SubQuasiArray{<:Any,2,<:AbstractInterlaceBasis,<:Tuple{Inclusion,AbstractVector{Int}}})
+    P = parent(V)
+    _,jr = parentindices(V)
+    J = findblock(axes(P,2),maximum(jr))
+    ProjectionFactorization(factorize(P[:,Block.(OneTo(Int(J)))]), jr)
+end

src/lanczos.jl

@@ -87,7 +87,6 @@ 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
@@ -157,6 +156,15 @@ struct LanczosPolynomial{T,Weight,Basis} <: OrthogonalPolynomial{T}
     data::LanczosData{T}
 end
 
+function LanczosPolynomial(w_in::AbstractQuasiVector, P::AbstractQuasiMatrix)
+    Q = normalize(P)
+    wQ = weighted(Q)
+    w = wQ * (wQ \ w_in) # expand weight in basis
+    LanczosPolynomial(w, Q, LanczosData(w, Q))
+end
+
+LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, orthonormalpolynomial(singularities(w)))
+
 ==(A::LanczosPolynomial, B::LanczosPolynomial) = A.w == B.w
 ==(::LanczosPolynomial, ::OrthogonalPolynomial) = false # TODO: fix
 ==(::OrthogonalPolynomial, ::LanczosPolynomial) = false # TODO: fix
@@ -172,14 +180,6 @@ orthogonalpolynomial(w::AbstractQuasiVector) = OrthogonalPolynomial(w)
 orthogonalpolynomial(w::SubQuasiArray) = orthogonalpolynomial(parent(w))[parentindices(w)[1],:]
 orthonormalpolynomial(w::AbstractQuasiVector) = normalize(orthogonalpolynomial(w))
 
-function LanczosPolynomial(w_in::AbstractQuasiVector, P::AbstractQuasiMatrix)
-    Q = normalize(P)
-    wQ = weighted(Q)
-    w = wQ * (wQ \ w_in) # expand weight in basis
-    LanczosPolynomial(w, Q, LanczosData(w, Q))
-end
-
-LanczosPolynomial(w::AbstractQuasiVector) = LanczosPolynomial(w, orthonormalpolynomial(singularities(w)))
 
 
 

test/runtests.jl

@@ -39,6 +39,7 @@ include("test_ratios.jl")
 include("test_normalized.jl")
 include("test_lanczos.jl")
 include("test_stieltjes.jl")
+include("test_interlace.jl")
 
 @testset "Auto-diff" begin
     U = Ultraspherical(1)

test/test_interlace.jl

@@ -0,0 +1,23 @@
+using ClassicalOrthogonalPolynomials, BlockArrays, LazyBandedMatrices, Test
+import ClassicalOrthogonalPolynomials: PiecewiseInterlace
+
+@testset "Piecewise" begin
+    T1,T2 = chebyshevt((-1)..0), chebyshevt(0..1)
+    U1,U2 = chebyshevu((-1)..0), chebyshevu(0..1)
+    T = PiecewiseInterlace(T1, T2)
+    U = PiecewiseInterlace(U1, U2)
+    D = U \ (Derivative(axes(T,1))*T)
+    C = U \ T
+
+    A = BlockVcat(T[-1,:]',
+                    BlockBroadcastArray(vcat,unitblocks(T1[end,:]),-unitblocks(T2[begin,:]))',
+                    D-C)
+    N = 20
+    M = BlockArray(A[Block.(1:N+1), Block.(1:N)])
+
+    u = M \ [exp(-1); zeros(size(M,1)-1)]
+    x = axes(T,1)
+
+    F = factorize(T[:,Block.(Base.OneTo(N))])
+    @test F \ exp.(x) ≈ (T \ exp.(x))[Block.(1:N)] ≈ u
+end

test/test_lanczos.jl

@@ -1,5 +1,5 @@
-using ClassicalOrthogonalPolynomials, BandedMatrices, ArrayLayouts, Test
-import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, orthogonalityweight, golubwelsch
+using ClassicalOrthogonalPolynomials, BandedMatrices, ArrayLayouts, QuasiArrays, ContinuumArrays, Test
+import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, orthogonalityweight, golubwelsch, LanczosData
 
 @testset "Lanczos" begin
     @testset "Legendre" begin
@@ -247,4 +247,15 @@ import ClassicalOrthogonalPolynomials: recurrencecoefficients, PaddedLayout, ort
         Q = LanczosPolynomial(@. x^m*(1-x)^m*(2-x)^m)
         @test Q[0.1,:]'*(Q \exp.(x)) ≈ exp(0.1)
     end
+
+    @testset "1/sqrt(1-x^2) + δ₂" begin
+        U = ChebyshevU()
+        W = π/2*I + (Base.unsafe_getindex(U,2,:) * Base.unsafe_getindex(U,2,:)')
+        X = jacobimatrix(U)
+        dat = LanczosData(X, W);
+        w = QuasiArrays.UnionVcat(ChebyshevUWeight(), DiracDelta(2))
+        Q = LanczosPolynomial(w, U, dat);
+        R = U \ Q;
+        @test R[1:5,1:5] isa Matrix{Float64}
+    end
 end