Support two-interval Hilbert (#62)

* Support two-interval Hilbert

* store each colsupport

* tests pass

* Update Project.toml

* Update Project.toml

* improve coverage

* Update test_stieltjes.jl

Author: dlfivefifty

Project.toml

@@ -29,7 +29,7 @@ ArrayLayouts = "0.7"
 BandedMatrices = "0.16.5"
 BlockArrays = "0.16"
 BlockBandedMatrices = "0.10"
-ContinuumArrays = "0.9"
+ContinuumArrays = "0.9.1"
 DomainSets = "0.5"
 FFTW = "1.1"
 FastGaussQuadrature = "0.4.3"
@@ -39,8 +39,8 @@ HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.11.1"
 InfiniteLinearAlgebra = "0.5.8"
 IntervalSets = "0.5"
-LazyArrays = "0.21.14"
-LazyBandedMatrices = "0.6"
+LazyArrays = "0.21.15"
+LazyBandedMatrices = "0.6.5"
 QuasiArrays = "0.7"
 SpecialFunctions = "1.0"
 julia = "1.6"

src/classical/chebyshev.jl

@@ -44,9 +44,11 @@ WeightedChebyshev() = WeightedChebyshevT()
 chebyshevt() = ChebyshevT()
 chebyshevt(d::AbstractInterval{T}) where T = ChebyshevT{float(T)}()[affine(d, ChebyshevInterval{T}()), :]
 chebyshevt(d::Inclusion) = chebyshevt(d.domain)
+chebyshevt(S::AbstractQuasiMatrix) = chebyshevt(axes(S,1))
 chebyshevu() = ChebyshevU()
 chebyshevu(d::AbstractInterval{T}) where T = ChebyshevU{float(T)}()[affine(d, ChebyshevInterval{T}()), :]
 chebyshevu(d::Inclusion) = chebyshevu(d.domain)
+chebyshevu(S::AbstractQuasiMatrix) = chebyshevu(axes(S,1))
 
 """
      chebyshevt(n, z)

src/interlace.jl

@@ -140,7 +140,7 @@ 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)
+    args = arguments.(*, 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

src/stieltjes.jl

@@ -63,7 +63,6 @@ const Hilbert{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(inv),Tuple{ConvKernel{T,I
 const LogKernel{T,D1,D2} = BroadcastQuasiMatrix{T,typeof(log),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}}}}
 const PowKernel{T,D1,D2,F<:Real} = BroadcastQuasiMatrix{T,typeof(^),Tuple{BroadcastQuasiMatrix{T,typeof(abs),Tuple{ConvKernel{T,Inclusion{T,D1},Inclusion{T,D2}}}},F}}
 
-
 @simplify function *(H::Hilbert{<:Any,<:ChebyshevInterval,<:ChebyshevInterval}, w::ChebyshevTWeight)
     T = promote_type(eltype(H), eltype(w))
     zeros(T, axes(w,1))
@@ -137,6 +136,19 @@ end
     end
 end
 
+@simplify function *(H::Hilbert, S::PiecewiseInterlace)
+    axes(H,2) == axes(S,1) || throw(DimensionMismatch())
+    @assert length(S.args) == 2
+    a,b = S.args
+    xa,xb = axes(a,1),axes(b,1)
+    Ha_a = inv.(xa .- xa') * a
+    Ha_b = inv.(xb .- xa') * a
+    Hb_a = inv.(xa .- xb') * b
+    Hb_b = inv.(xb .- xb') * b
+    c,d = Hb_a.args[1], Ha_b.args[1]
+    A,B,C,D = unitblocks(c \ Ha_a), unitblocks(c \ Hb_a), unitblocks(d \ Ha_b), unitblocks(d \ Hb_b)
+    PiecewiseInterlace(c,d)  * BlockBroadcastArray{promote_type(eltype(H),eltype(S))}(hvcat, 2, A, B, C, D)
+end
 
 ### 
 # LogKernel
@@ -216,18 +228,28 @@ mutable struct HilbertVandermonde{T,MM} <: AbstractCachedMatrix{T}
     M::MM
     data::Matrix{T}
     datasize::NTuple{2,Int}
+    colsupport::Vector{Int}
 end
 
-HilbertVandermonde(M, data::Matrix) = HilbertVandermonde(M, data, size(data))
+HilbertVandermonde(M, data::Matrix) = HilbertVandermonde(M, data, size(data), Int[])
 size(H::HilbertVandermonde) = (ℵ₀,ℵ₀)
+function colsupport(H::HilbertVandermonde, j)
+    resizedata!(H, H.datasize[1], maximum(j))
+    1:maximum(H.colsupport[j])
+end
+
+copy(H::HilbertVandermonde) = HilbertVandermonde(H.M, copy(H.data), H.datasize, H.colsupport)
 
 function cache_filldata!(H::HilbertVandermonde{T}, kr, jr) where T
     n,m = H.datasize
-    isempty(kr) && return
     isempty(jr) && return
-    H.data[(n+1):maximum(kr),1:m] .= zero(T)
+    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)
-        H.data[kr,j] .= (H.M * [H.data[:,j-1]; Zeros{T}(∞)])[kr]
+        u = H.M * [H.data[:,j-1]; Zeros{T}(∞)]
+        H.colsupport[j] = maximum(colsupport(u,1))
+        isempty(kr) || (H.data[kr,j] .= u[kr])
     end
 end
 

test/test_chebyshev.jl

@@ -111,12 +111,13 @@ import ContinuumArrays: MappedWeightedBasisLayout, Map
         T = Chebyshev()[2x .- 1,:]
         @test (T*(T\x))[0.1] ≈ 0.1
         @test (T* (T \ exp.(x)))[0.1] ≈ exp(0.1)
-        @test chebyshevt(0..1) == T
+        @test chebyshevt(0..1) == chebyshevt(Inclusion(0..1)) == chebyshevt(T) == T
 
         Tn = Chebyshev()[2x .- 1, [1,3,4]]
         @test (axes(Tn,1) .* Tn).args[2][1:5,:] ≈ (axes(T,1) .* T).args[2][1:5,[1,3,4]]
 
         U = chebyshevu(0..1)
+        @test chebyshevu(0..1) == chebyshevu(Inclusion(0..1)) == chebyshevu(T) == U
         @test (U*(U\x))[0.1] ≈ 0.1
         @test (U* (U \ exp.(x)))[0.1] ≈ exp(0.1)
 

test/test_interlace.jl

@@ -2,22 +2,46 @@ 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
+    @testset "two-interval ODE" 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)])
+        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)
+        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
+        F = factorize(T[:,Block.(Base.OneTo(N))])
+        @test F \ exp.(x) ≈ (T \ exp.(x))[Block.(1:N)] ≈ u
+    end
+
+    @testset "two-interval Weighted Derivative" begin
+        T1,T2 = chebyshevt((-2)..(-1)), chebyshevt(0..1)
+        U1,U2 = chebyshevu((-2)..(-1)), chebyshevu(0..1)
+        W = PiecewiseInterlace(Weighted(T1), Weighted(T2))
+        U = PiecewiseInterlace(U1, U2)
+        x = axes(W,1)
+        D = Derivative(x)
+        @test_broken U\D*W isa BlockBroadcastArray
+    end
+
+    @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)) 
+        T = PiecewiseInterlace(T1, T2)
+        x = axes(W,1)
+        H = T \ inv.(x .- x') * W;
+
+        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
+    end
 end

test/test_stieltjes.jl

@@ -196,7 +196,9 @@ end
             t = axes(U,1)
             x = Inclusion(2..3)
             T = chebyshevt(2..3)
-            H = T \ inv.(x .- t') * W
+            H = T \ inv.(x .- t') * W;
+            @test last(colsupport(H,6)) ≤ 40
+            @test T[2.3,1:100]'*(H * (W \ @.(sqrt(1-t^2)exp(t))))[1:100] ≈ 0.9068295340935111
             @test T[2.3,1:100]' * H[1:100,1:100] ≈ (inv.(2.3 .- t') * W)[:,1:100]
         end
 
@@ -207,7 +209,16 @@ end
             x = Inclusion(2..3)
             T = chebyshevt(2..3)
             H = T \ inv.(x .- t') * W
-            @test T[2.3,1:100]' * H[1:100,1:100] ≈ (inv.(2.3 .- t') * W)[:,1:100]
+            N = 100
+            @test T[2.3,1:N]' * H[1:N,1:N] ≈ (inv.(2.3 .- t') * W)[:,1:N]
+
+            U = chebyshevu((-2)..(-1))
+            W = Weighted(U)
+            T = chebyshevt(0..2)
+            x = axes(T,1)
+            t = axes(W,1)
+            H = T \ inv.(x .- t') * W
+            @test T[0.5,1:N]'*(H * (W \ @.(sqrt(-1-t)*sqrt(t+2)*exp(t))))[1:N] ≈ 0.047390454610749054
         end
     end
 end