mapped Hilbert transform (#58)

* mapped Hilbert transform

* Mapped log kernel

* Update test_stieltjes.jl

Author: dlfivefifty

src/ClassicalOrthogonalPolynomials.jl

@@ -323,7 +323,6 @@ end
 
 
 function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}) where T
-    x,w = golubwelsch(L)
     Q = Normalized(parent(L))[parentindices(L)...]
     D = L \ Q
     F = factorize(Q)

src/stieltjes.jl

@@ -71,13 +71,13 @@ end
 
 @simplify function *(H::Hilbert, w::ChebyshevUWeight)
     T = promote_type(eltype(H), eltype(w))
-    fill(convert(T,π), axes(w,1))
+    convert(T,π) * axes(w,1)
 end
 
 @simplify function *(H::Hilbert, w::LegendreWeight)
     T = promote_type(eltype(H), eltype(w))
     x = axes(w,1)
-    log.(x .+ 1) .- log.(1 .- x)
+    log.(x .+ one(T)) .- log.(one(T) .- x)
 end
 
 @simplify function *(H::Hilbert, wT::Weighted{<:Any,<:ChebyshevT}) 
@@ -95,18 +95,48 @@ end
     P = wP.P
     w = orthogonalityweight(P)
     A = recurrencecoefficients(P)[1]
-    (-A[1]*sum(w))*[zero(axes(P,1)) associated(P)] + (H*w) .* P
+    Q = associated(P)
+    (-A[1]*sum(w))*[zero(axes(P,1)) Q] + (H*w) .* P
 end
 
 @simplify *(H::Hilbert, P::Legendre) = H * Weighted(P)
 
+
+
+@simplify function *(H::Hilbert, w::SubQuasiArray{<:Any,1})
+    T = promote_type(eltype(H), eltype(w))
+    m = parentindices(w)[1]
+    P = parent(w)
+    x = axes(P,1)
+    (inv.(x .- x') * P)[m]
+end
+
+@simplify function *(H::Hilbert, wP::Weighted{<:Any,<:SubQuasiArray{<:Any,2}})
+    T = promote_type(eltype(H), eltype(wP))
+    kr,jr = parentindices(wP.P)
+    P = parent(wP.P)
+    x = axes(P,1)
+    (inv.(x .- x') * Weighted(P))[kr,jr]
+end
+
+
 ### 
 # LogKernel
 ###
 
 @simplify function *(L::LogKernel, wT::Weighted{<:Any,<:ChebyshevT}) 
     T = promote_type(eltype(L), eltype(wT))
-    ChebyshevT{T}() * Diagonal(Vcat(-π*log(2*one(T)),-convert(T,π)./(1:∞)))
+    ChebyshevT{T}() * Diagonal(Vcat(-convert(T,π)*log(2*one(T)),-convert(T,π)./(1:∞)))
+end
+
+@simplify function *(H::LogKernel, wT::Weighted{<:Any,<:SubQuasiArray{<:Any,2,<:ChebyshevT,<:Tuple{AbstractAffineQuasiVector,Slice}}})
+    V = promote_type(eltype(H), eltype(wT))
+    kr,jr = parentindices(wT.P)
+    T = parent(wT.P)
+    x = axes(T,1)
+    W = Weighted(T)
+    A = kr.A
+    T[kr,:] * Diagonal(Vcat(-convert(V,π)*(log(2*one(V))+log(abs(A)))/A,-convert(V,π)./(A * (1:∞))))
 end
 
 

test/test_stieltjes.jl

@@ -1,4 +1,4 @@
-using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, Test
+using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, BandedMatrices, Test
 import ClassicalOrthogonalPolynomials: Hilbert, StieltjesPoint, ChebyshevInterval, associated, Associated, orthogonalityweight, Weighted, gennormalizedpower, *, dot, PowerLawMatrix, PowKernelPoint
 import InfiniteArrays: I
 
@@ -23,6 +23,26 @@ end
 
 
 @testset "Singular integrals" begin
+    @testset "weights" begin
+        w_T = ChebyshevTWeight()
+        w_U = ChebyshevUWeight()
+        w_P = LegendreWeight()
+        x = axes(w_T,1)
+        H = inv.(x .- x')
+        @test iszero(H*w_T)
+        @test (H*w_U)[0.1] ≈ π/10
+        @test (H*w_P)[0.1] ≈ log(1.1) - log(1-0.1)
+
+        w_T = orthogonalityweight(chebyshevt(0..1))
+        w_U = orthogonalityweight(chebyshevu(0..1))
+        w_P = orthogonalityweight(legendre(0..1))
+        x = axes(w_T,1)
+        H = inv.(x .- x')
+        @test iszero(H*w_T)
+        @test (H*w_U)[0.1] ≈ (2*0.1-1)*π
+        @test (H*w_P)[0.1] ≈ (log(1+(-0.8)) - log(1-(-0.8)))
+    end
+
     @testset "Stieltjes" begin
         T = Chebyshev()
         wT = ChebyshevWeight() .* T
@@ -36,7 +56,7 @@ end
 
         f = wT * [[1,2,3]; zeros(∞)];
         J = T \ (x .* T)
-    
+
         @test π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f)[1]
         @test π*((z*I-J) \ f.args[2])[1,1] ≈ (S*f.args[1]*f.args[2])[1]
 
@@ -55,7 +75,7 @@ end
             x = axes(T,1)
             @test inv.(t .- x') * Weighted(T) ≈ inv.((t+eps()im) .- x') * Weighted(T)
             @test (inv.(t .- x') * Weighted(U))[1:10] ≈ (inv.((t+eps()im) .- x') * Weighted(U))[1:10]
-            
+
             t = 0.5
             @test_broken inv.(t .- x') * Weighted(T)
             @test_broken inv.(t .- x') * Weighted(U)
@@ -71,7 +91,7 @@ end
 
         @testset "weights" begin
             @test H * ChebyshevTWeight() ≡ QuasiZeros{Float64}((x,))
-            @test H * ChebyshevUWeight() ≡ QuasiFill(1.0π,(x,))
+            @test H * ChebyshevUWeight() == π*x
             @test (H * LegendreWeight())[0.1] ≈ log((0.1+1)/(1-0.1))
         end
 
@@ -82,15 +102,16 @@ end
         @test (Ultraspherical(1) \ (H*wT) * (wT \ wU))[1:10,1:10] ==
                     ((Ultraspherical(1) \ Chebyshev()) * (Chebyshev() \ (H*wU)))[1:10,1:10]
 
-        # Other axes
-        x = Inclusion(0..1)
-        y = 2x .- 1
-        H = inv.(x .- x')
+        @testset "Other axes" begin
+            x = Inclusion(0..1)
+            y = 2x .- 1
+            H = inv.(x .- x')
 
-        wT2 = wT[y,:]
-        wU2 = wU[y,:]
-        @test (Ultraspherical(1)[y,:]\(H*wT2))[1:10,1:10] == diagm(1 => fill(-π,9))
-        @test (Chebyshev()[y,:]\(H*wU2))[1:10,1:10] == diagm(-1 => fill(1.0π,9))
+            wT2 = wT[y,:]
+            wU2 = wU[y,:]
+            @test (Ultraspherical(1)[y,:]\(H*wT2))[1:10,1:10] == diagm(1 => fill(-π,9))
+            @test (Chebyshev()[y,:]\(H*wU2))[1:10,1:10] == diagm(-1 => fill(1.0π,9))
+        end
 
         @testset "Legendre" begin
             P = Legendre()
@@ -99,6 +120,14 @@ end
             X = jacobimatrix(P)
             @test Q[0.1,1:11]'*X[1:11,1:10] ≈ (0.1 * Array(Q[0.1,1:10])' - [2 zeros(1,9)])
         end
+
+        @testset "mapped" begin
+            T = chebyshevt(0..1)
+            U = chebyshevu(0..1)
+            x = axes(T,1)
+            H = inv.(x .- x')
+            @test U\H*Weighted(T) isa BandedMatrix
+        end
     end
 
     @testset "Log kernel" begin
@@ -122,6 +151,13 @@ end
         u =  wT * (2 *(T \ exp.(x)))
         @test u[0.1] ≈ exp(0.1)/sqrt(0.1-0.1^2)
         @test (L * u)[0.5] ≈ -7.471469928754152 # Mathematica
+
+        @testset "mapped" begin
+            T = chebyshevt(0..1)
+            x = axes(T,1)
+            L = log.(abs.(x .- x'))
+            @test T[0.2,:]'*((T\L*Weighted(T)) * (T\exp.(x))) ≈ -2.9976362326874373 # Mathematica
+        end
     end
 
     @testset "pow kernel" begin
@@ -143,7 +179,7 @@ end
 
         u = Weighted(U) * ((H * Weighted(U)) \ imag(c * x))
 
-        ε  = eps(); 
+        ε  = eps();
         @test (inv.(0.1+ε*im .- x') * u + inv.(0.1-ε*im .- x') * u)/2 ≈ imag(c*0.1)
         @test real(inv.(0.1+ε*im .- x') * u ) ≈ imag(c*0.1)
 
@@ -284,4 +320,4 @@ end
         @time D = ClassicalOrthogonalPolynomials.LanczosData((1.001 .- x).^0.5, P);
         @time ClassicalOrthogonalPolynomials.resizedata!(D,100);
     end
-end
+end