Merge branch 'main' into pr/63

Author: dlfivefifty

Project.toml

@@ -1,7 +1,7 @@
 name = "ClassicalOrthogonalPolynomials"
 uuid = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.4.5"
+version = "0.4.6"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -28,7 +28,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ArrayLayouts = "0.7"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.16"
-BlockBandedMatrices = "0.10"
+BlockBandedMatrices = "0.10.9"
 ContinuumArrays = "0.9.1"
 DomainSets = "0.5"
 FFTW = "1.1"
@@ -39,8 +39,8 @@ HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.11.1"
 InfiniteLinearAlgebra = "0.5.8"
 IntervalSets = "0.5"
-LazyArrays = "0.21.15"
-LazyBandedMatrices = "0.6.5"
+LazyArrays = "0.21.16"
+LazyBandedMatrices = "0.6.6"
 QuasiArrays = "0.7"
 SpecialFunctions = "1.0"
 julia = "1.6"

src/stieltjes.jl

@@ -57,7 +57,8 @@ associated(::ChebyshevT{T}) where T = ChebyshevU{T}()
 associated(::ChebyshevU{T}) where T = ChebyshevU{T}()
 
 
-const StieltjesPoint{T,V,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{T,QuasiAdjoint{V,Inclusion{V,D}}}}}}
+const StieltjesPoint{T,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{BroadcastQuasiMatrix{T,typeof(-),Tuple{W,QuasiAdjoint{V,Inclusion{V,D}}}}}}
+const LogKernelPoint{T<:Real,C,W<:Number,V,D} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{BroadcastQuasiMatrix{C,typeof(-),Tuple{W,QuasiAdjoint{V,Inclusion{V,D}}}}}}}}
 const ConvKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(-),Tuple{D1,QuasiAdjoint{T,D2}}}
 const Hilbert{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}
 const LogKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}}}
@@ -79,18 +80,18 @@ end
     log.(x .+ one(T)) .- log.(one(T) .- x)
 end
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT}) 
+@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT})
     T = promote_type(eltype(H), eltype(wT))
     ChebyshevU{T}() * _BandedMatrix(Fill(-convert(T,π),1,∞), ℵ₀, -1, 1)
 end
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wU::Weighted{<:Any,<:ChebyshevU}) 
+@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wU::Weighted{<:Any,<:ChebyshevU})
     T = promote_type(eltype(H), eltype(wU))
     ChebyshevT{T}() * _BandedMatrix(Fill(convert(T,π),1,∞), ℵ₀, 1, -1)
 end
 
 
-@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wP::Weighted{<:Any,<:OrthogonalPolynomial}) 
+@simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wP::Weighted{<:Any,<:OrthogonalPolynomial})
     P = wP.P
     w = orthogonalityweight(P)
     A = recurrencecoefficients(P)[1]
@@ -150,11 +151,11 @@ end
     PiecewiseInterlace(c,d)  * BlockBroadcastArray{promote_type(eltype(H),eltype(S))}(hvcat, 2, A, B, C, D)
 end
 
-### 
+###
 # LogKernel
 ###
 
-@simplify function *(L::LogKernel{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT}) 
+@simplify function *(L::LogKernel{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, wT::Weighted{<:Any,<:ChebyshevT})
     T = promote_type(eltype(L), eltype(wT))
     ChebyshevT{T}() * Diagonal(Vcat(-convert(T,π)*log(2*one(T)),-convert(T,π)./(1:∞)))
 end
@@ -172,11 +173,11 @@ end
 
 
 
-### 
+###
 # PowKernel
 ###
 
-@simplify function *(K::PowKernel, wT::Weighted{<:Any,<:Jacobi}) 
+@simplify function *(K::PowKernel, wT::Weighted{<:Any,<:Jacobi})
     T = promote_type(eltype(K), eltype(wT))
     cnv,α = K.args
     x,y = K.args[1].args[1].args
@@ -203,8 +204,8 @@ end
     P = wP.P
     w = orthogonalityweight(P)
     X = jacobimatrix(P)
-    z, x = parent(S).args[1].args
-    z in axes(P,1) && transpose((inv.(x .- x') * wP)[z,:])
+    z, xc = parent(S).args[1].args
+    z in axes(P,1) && return transpose(view(inv.(xc' .- xc) * wP,z,:))
     transpose((X'-z*I) \ [-sum(w)*_p0(P); zeros(∞)])
 end
 
@@ -213,12 +214,35 @@ sqrtx2(x::Real) = sign(x)*sqrt(x^2-1)
 
 @simplify function *(S::StieltjesPoint, wP::Weighted{<:Any,<:ChebyshevU})
     T = promote_type(eltype(S), eltype(wP))
-    z, x = parent(S).args[1].args
-    z in axes(wP,1) && transpose((inv.(x .- x') * wP)[z,:])
+    z, xc = parent(S).args[1].args
+    z in axes(wP,1) && return  (convert(T,π)*ChebyshevT()[z,2:end])'
     ξ = inv(z + sqrtx2(z))
     transpose(convert(T,π) * ξ.^oneto(∞))
 end
 
+####
+# LogKernelPoint
+####
+
+@simplify function *(L::LogKernelPoint, wP::Weighted{<:Any,<:ChebyshevU})
+    T = promote_type(eltype(L), eltype(wP))
+    z, xc = parent(L).args[1].args[1].args
+    if z in axes(wP,1)
+        Tn = Vcat(convert(T,π)*log(2*one(T)), convert(T,π)*ChebyshevT()[z,2:end]./oneto(∞))
+        return transpose((Tn[3:end]-Tn[1:end])/2)
+    else
+        # for U_k where k>=1
+        ξ = inv(z + sqrtx2(z))
+        ζ = (convert(T,π)*ξ.^oneto(∞))./oneto(∞)
+        ζ = (ζ[3:end]- ζ[1:end])/2
+
+        # for U_0
+        ζ = Vcat(convert(T,π)*(ξ^2/4 - (log.(abs.(ξ)) + log(2*one(T)))/2), ζ)
+        return transpose(ζ)
+    end
+
+end
+
 """
    HilbertVandermonde(M, data)
 
@@ -231,7 +255,7 @@ mutable struct HilbertVandermonde{T,MM} <: AbstractCachedMatrix{T}
     colsupport::Vector{Int}
 end
 
-HilbertVandermonde(M, data::Matrix) = HilbertVandermonde(M, data, size(data), Int[])
+HilbertVandermonde(M, data::Matrix) = HilbertVandermonde(M, data, size(data), fill(size(data,1), size(data,2)))
 size(H::HilbertVandermonde) = (ℵ₀,ℵ₀)
 function colsupport(H::HilbertVandermonde, j)
     resizedata!(H, H.datasize[1], maximum(j))
@@ -244,7 +268,7 @@ function cache_filldata!(H::HilbertVandermonde{T}, kr, jr) where T
     n,m = H.datasize
     isempty(jr) && return
     resize!(H.colsupport, max(length(H.colsupport), maximum(jr)))
-    
+
     isempty(kr) || (H.data[(n+1):maximum(kr),1:m] .= zero(T))
     for j in (m+1):maximum(jr)
         u = H.M * [H.data[:,j-1]; Zeros{T}(∞)]
@@ -256,8 +280,9 @@ end
 @simplify function *(H::Hilbert{<:Any,<:Any,<:ChebyshevInterval}, W::Weighted{<:Any,<:ChebyshevU})
     x = axes(H,1)
     T̃ = chebyshevt(x)
-    ψ_1 = T̃ \ inv.(x .+ sqrtx2.(x))
+    ψ_1 = T̃ \ inv.(x .+ sqrtx2.(x)) # same ψ_1 = x .- sqrt(x^2 - 1) but with relative accuracy as x -> ∞
     data = convert(eltype(H),π) * Matrix(reshape(paddeddata(ψ_1),:,1))
+    # Operator has columns π * ψ_1^k
     T̃ * HilbertVandermonde(Clenshaw(T̃ * ψ_1, T̃), data)
 end
 
@@ -272,14 +297,14 @@ end
     (inv.(z̃ .- x̃') * P)[:,parentindices(wT)[2]]
 end
 
-@simplify function *(H::Hilbert, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) 
+@simplify function *(H::Hilbert, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
     P = parent(wT)
     x = axes(P,1)
     apply(*, inv.(x .- x'), P)[parentindices(wT)...]
 end
 
 
-@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}}) 
+@simplify function *(L::LogKernel, wT::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
     V = promote_type(eltype(L), eltype(wT))
     wP = parent(wT)
     kr, jr = parentindices(wT)
@@ -294,7 +319,7 @@ end
 
 ### generic fallback
 for Op in (:Hilbert, :StieltjesPoint, :LogKernel, :PowKernel)
-    @eval @simplify function *(H::$Op, wP::WeightedBasis{<:Any,<:Weight,<:Any}) 
+    @eval @simplify function *(H::$Op, wP::WeightedBasis{<:Any,<:Weight,<:Any})
         w,P = wP.args
         Q = OrthogonalPolynomial(w)
         (H * Weighted(Q)) * (Q \ P)
@@ -496,4 +521,4 @@ function dot(v::AbstractVector{T}, W::PowerLawMatrix, q::AbstractVector{T}) wher
     vpad, qpad = paddeddata(v), paddeddata(q)
     vl, ql = length(vpad), length(qpad)
     return dot(vpad,W[1:vl,1:ql]*qpad)
-end
+end

test/test_interlace.jl

@@ -1,4 +1,4 @@
-using ClassicalOrthogonalPolynomials, BlockArrays, LazyBandedMatrices, Test
+using ClassicalOrthogonalPolynomials, BlockArrays, LazyBandedMatrices, FillArrays, Test
 import ClassicalOrthogonalPolynomials: PiecewiseInterlace
 
 @testset "Piecewise" begin
@@ -36,12 +36,37 @@ import ClassicalOrthogonalPolynomials: PiecewiseInterlace
     @testset "two-interval Hilbert" begin
         T1,T2 = chebyshevt((-2)..(-1)), chebyshevt(0..2)
         U1,U2 = chebyshevu((-2)..(-1)), chebyshevu(0..2)
-        W = PiecewiseInterlace(Weighted(U1), Weighted(U2)) 
+        W = PiecewiseInterlace(Weighted(U1), Weighted(U2))
         T = PiecewiseInterlace(T1, T2)
+        U = PiecewiseInterlace(U1, U2)
         x = axes(W,1)
         H = T \ inv.(x .- x') * W;
 
+        @test maximum(BlockArrays.blockcolsupport(H,Block(5))) ≤ Block(50)
+
         c = W \ broadcast(x -> exp(x)* (0 ≤ x ≤ 2 ? sqrt(2-x)*sqrt(x) : sqrt(-1-x)*sqrt(x+2)), x)
         @test T[0.5,1:200]'*(H*c)[1:200] ≈ -6.064426633490422
+
+        @testset "inversion" begin
+            H̃ = BlockHcat(Eye((axes(H,1),))[:,Block(1)], H)
+            @test BlockArrays.blockcolsupport(H̃,1) == Block.(1:1)
+            @test BlockArrays.blockcolsupport(H̃,2) == Block.(1:22)
+
+            UT = U \ T
+            D = U \ Derivative(x) * T
+            V = x -> x^4 - 10x^2
+            V_cfs = T \ V.(x)
+            Vp_cfs_U = D * V_cfs
+
+            N = 100
+            Vp_cfs_N = UT[Block.(1:N),Block.(1:N)] \ Vp_cfs_U[Block.(1:N)]
+
+            cμ = H̃[Block.(1:N), Block.(1:N)] \ Vp_cfs_N;
+            c1,c2 = cμ[Block(1)]
+            μ = W[:,Block.(1:N-1)] * cμ[Block.(2:N)]/2;
+
+            # H * μ == Vp(x) + c1 on first interval
+            # H * μ == Vp(x) + c2 on second interval
+        end
     end
 end