Unify logkernel and stieltjes expansions (#34)

Project.toml

@@ -1,7 +1,7 @@
 name = "SingularIntegrals"
 uuid = "d7440221-8b5e-42fc-909c-0567823f424a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.3"
+version = "0.2.4"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/SingularIntegrals.jl

@@ -1,6 +1,6 @@
 module SingularIntegrals
 using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FillArrays, BandedMatrices, LinearAlgebra, SpecialFunctions, HypergeometricFunctions, InfiniteArrays
-using ContinuumArrays: @simplify, Weight, AbstractAffineQuasiVector, inbounds_getindex, broadcastbasis, MappedBasisLayouts, MemoryLayout, MappedWeightLayout, ExpansionLayout, demap, basismap, AbstractBasisLayout, SubBasisLayout
+using ContinuumArrays: @simplify, Weight, AbstractAffineQuasiVector, inbounds_getindex, broadcastbasis, MappedBasisLayouts, MemoryLayout, MappedWeightLayout, AbstractWeightLayout, ExpansionLayout, demap, basismap, AbstractBasisLayout, SubBasisLayout
 using QuasiArrays: AbstractQuasiMatrix, BroadcastQuasiMatrix, LazyQuasiArrayStyle, AbstractQuasiVecOrMat
 import ClassicalOrthogonalPolynomials: AbstractJacobiWeight, WeightedBasis, jacobimatrix, orthogonalityweight, recurrencecoefficients, _p0, Clenshaw, chop, initiateforwardrecurrence, MappedOPLayouts, unweighted
 using LazyBandedMatrices: Tridiagonal, SymTridiagonal, subdiagonaldata, supdiagonaldata, diagonaldata, ApplyLayout
@@ -32,4 +32,32 @@ for Op in (:Stieltjes, :StieltjesPoint, :LogKernelPoint, :PowerKernelPoint, :Log
 end
 
 
+
+# general routines
+for lk in (:logkernel, :complexlogkernel, :stieltjes)
+    lk_layout = Symbol(lk, :_layout)
+    @eval begin
+        $lk_layout(::AbstractBasisLayout, P, z...) = error("not implemented")
+        $lk_layout(::AbstractWeightLayout, w, zs::AbstractVector) = [stieltjes(w, z) for z in zs]
+        function $lk_layout(::AbstractWeightLayout, w, z::Inclusion)
+            axes(w,1) == z || error("Not implemented")
+            $lk(w)
+        end
+        function $lk_layout(::AbstractBasisLayout, w, z::Inclusion)
+            axes(w,1) == z || error("Not implemented")
+            $lk(w)
+        end
+        $lk_layout(lay, P, z...) = $lk(expand(P), z...)
+
+        function $lk_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, y...)
+            a = arguments(LAY, V)
+            *($lk(a[1], y...), tail(a)...)
+        end
+
+        $lk_layout(::ExpansionLayout, A, dims...) = $lk_layout(ApplyLayout{typeof(*)}(), A, dims...)
+        $lk_layout(::SubBasisLayout, A, dims...) = $lk(parent(A), dims...)[:, parentindices(A)[2]]
+    end
+end
+
+
 end # module SingularIntegrals

src/logkernel.jl

@@ -40,22 +40,6 @@ applies the log kernel log(x-t) to the columns of a quasi matrix, i.e., `(log.(x
 """
 complexlogkernel(P, z...) = complexlogkernel_layout(MemoryLayout(P), P, z...)
 
-for lk in (:logkernel, :complexlogkernel)
-    lk_layout = Symbol(lk, :_layout)
-    @eval begin
-        $lk_layout(::AbstractBasisLayout, P, z...) = error("not implemented")
-        $lk_layout(lay, P, z...) = $lk(expand(P), z...)
-
-        function $lk_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, y...)
-            a = arguments(LAY, V)
-            *($lk(a[1], y...), tail(a)...)
-        end
-
-        $lk_layout(::ExpansionLayout, A, dims...) = $lk_layout(ApplyLayout{typeof(*)}(), A, dims...)
-        $lk_layout(::SubBasisLayout, A, dims...) = $lk(parent(A), dims...)[:, parentindices(A)[2]]
-    end
-end
-
 logkernel(wT::Weighted{T,<:ChebyshevT}) where T = ChebyshevT{T}() * Diagonal(Vcat(-convert(T,π)*log(2*one(T)),-convert(T,π)./(1:∞)))
 function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::AbstractQuasiMatrix{V}) where V
     kr = basismap(wT)

src/stieltjes.jl

@@ -87,21 +87,6 @@ end
 computes inv.(y - x') * P understood in a principle value sense.
 """
 stieltjes(P, y...) = stieltjes_layout(MemoryLayout(P), P, y...)
-stieltjes_layout(lay, P, y) = stieltjes_size(size(P), P, y)
-function stieltjes_size(sz, P, y)
-    axes(P,1) == y && return stieltjes(P)
-    error("Not implemented")
-end
-
-stieltjes_size(::Tuple{Any}, P, zs::AbstractVector) = [stieltjes(P, z) for z in zs]
-
-function stieltjes_layout(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVecOrMat, y...)
-    a = arguments(LAY, V)
-    *(stieltjes(a[1], y...), tail(a)...)
-end
-
-stieltjes_layout(::ExpansionLayout, A, y...) = stieltjes_layout(ApplyLayout{typeof(*)}(), A, y...)
-
 
 """
     stieltjes(P)

test/test_stieltjes.jl

@@ -223,4 +223,10 @@ end
         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
+
+    @testset "expand" begin
+        P = Legendre()
+        @test stieltjes(P[:,1], 3) ≈ stieltjes(P,3)[1]
+        @test stieltjes(P[:,1:3], 3) ≈ stieltjes(P,3)[:,1:3]
+    end
 end