From 1b55e9c81ab73ab08d315b165ab43fdde67d1c2b Mon Sep 17 00:00:00 2001 From: Sheehan Olver Date: Fri, 19 Apr 2024 17:37:10 +0100 Subject: [PATCH] Fix #38 (#40) * logkernel defaults to complexlogkernel * Update Project.toml * better broadcasting with real * Weighted(ChebyshevU) now uses RecurrenceArray * Tests pass * Weighted Chebyshev T * split out recurrence from IC code * Fix recurrence forward mode bug --- .github/workflows/ci.yml | 2 +- Project.toml | 14 ++--- src/SingularIntegrals.jl | 9 ++- src/logkernel.jl | 133 +++++++++++++++++++++++++-------------- src/recurrence.jl | 13 +++- test/runtests.jl | 2 +- test/test_logkernel.jl | 18 +++++- 7 files changed, 130 insertions(+), 61 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2f7697c..7229ef4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -20,7 +20,7 @@ jobs: fail-fast: false matrix: version: - - '1.9' + - '1.10' os: - ubuntu-latest - macOS-latest diff --git a/Project.toml b/Project.toml index 2c47563..175debf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SingularIntegrals" uuid = "d7440221-8b5e-42fc-909c-0567823f424a" authors = ["Sheehan Olver "] -version = "0.2.6" +version = "0.3" [deps] ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a" @@ -20,18 +20,18 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" [compat] ArrayLayouts = "1.4" -BandedMatrices = "0.17.17, 1" -ClassicalOrthogonalPolynomials = "0.11, 0.12" -ContinuumArrays = "0.16, 0.17" -FastTransforms = "0.15" +BandedMatrices = "1" +ClassicalOrthogonalPolynomials = "0.12" +ContinuumArrays = "0.17" +FastTransforms = "0.15, 0.16" FillArrays = "1" HypergeometricFunctions = "0.3.4" InfiniteArrays = "0.13" -LazyArrays = "1.8" +LazyArrays = "1.10" LazyBandedMatrices = "0.9" QuasiArrays = "0.11" SpecialFunctions = "1, 2" -julia = "1.9" +julia = "1.10" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/SingularIntegrals.jl b/src/SingularIntegrals.jl index 7291ace..4e4d6c4 100644 --- a/src/SingularIntegrals.jl +++ b/src/SingularIntegrals.jl @@ -2,9 +2,9 @@ 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, 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 +import ClassicalOrthogonalPolynomials: AbstractJacobiWeight, WeightedBasis, jacobimatrix, orthogonalityweight, recurrencecoefficients, _p0, Clenshaw, chop, initiateforwardrecurrence, MappedOPLayouts, unweighted, WeightedOPLayout, MappedOPLayout using LazyBandedMatrices: Tridiagonal, SymTridiagonal, subdiagonaldata, supdiagonaldata, diagonaldata, ApplyLayout -import LazyArrays: AbstractCachedMatrix, AbstractCachedArray, paddeddata, arguments, resizedata!, cache_filldata!, zero!, cacheddata +import LazyArrays: AbstractCachedMatrix, AbstractCachedArray, paddeddata, arguments, resizedata!, cache_filldata!, zero!, cacheddata, LazyArrayStyle import Base: *, +, -, /, \, Slice, axes, getindex, sum, ==, oneto, size, broadcasted, copy, tail, view import LinearAlgebra: dot using BandedMatrices: _BandedMatrix @@ -32,12 +32,15 @@ for Op in (:Stieltjes, :StieltjesPoint, :LogKernelPoint, :PowerKernelPoint, :Log end +for lk in (:complexlogkernel, :stieltjes) + lk_layout = Symbol(lk, :_layout) + @eval $lk_layout(::AbstractBasisLayout, P, z...) = error("not implemented") +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") diff --git a/src/logkernel.jl b/src/logkernel.jl index 0243348..aab959c 100644 --- a/src/logkernel.jl +++ b/src/logkernel.jl @@ -41,7 +41,8 @@ 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...) 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 +function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT) + V = eltype(wT) kr = basismap(wT) @assert kr isa AbstractAffineQuasiVector W = demap(wT) @@ -53,8 +54,7 @@ function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::Abst unweighted(wT) * Diagonal(Vcat(-convert(V,π)*(log(2*one(V))+log(abs(A)))/A,-convert(V,π) ./ (A * (1:∞)))) end - -function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::AbstractQuasiMatrix{V}, z::Number) where V +function logkernel_demap(wT, z) P = demap(wT) kr = basismap(wT) z̃ = inbounds_getindex(kr, z) @@ -65,6 +65,9 @@ function logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT::Abst end +logkernel_layout(::Union{MappedBasisLayouts, MappedOPLayouts}, wT, z::Number) = logkernel_demap(wT, z) +logkernel_layout(::WeightedOPLayout{MappedOPLayout}, wT, z::Real) = logkernel_demap(wT, z) + @@ -74,71 +77,109 @@ end #### -function logkernel(wP::Weighted{T,<:ChebyshevU}, z::Number) where T - if z in axes(wP,1) - Tn = Vcat(convert(T,π)*log(2*one(T)), convert(T,π)*ChebyshevT{T}()[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 +function complexlogkernel_ics(wP::Weighted{<:Any,<:ChebyshevT}, z::Number) + V = promote_type(eltype(wP), typeof(z)) + ξ = inv(z + sqrtx2(z)) + r0 = convert(V,π)*(-log(ξ)-log(convert(V,2))) + r1 = -convert(V,π)*ξ + r2 = convert(V,π)/2 + z*r1 + r0,r1,r2 +end + +function complexlogkernel_recurrence(wP::Weighted{<:Any,<:ChebyshevT}) + # We have for n ≠ 0 L_n(z) = stieltjes(U_n, z) + # where U_n(z) = ∫_1^x T_n(x)/sqrt(1-x^2) dx = -(1-x^2)^(-1/2)T_{n-1}(x)/n + # We have the 3-term recurrence + # T_{n+1}(x) == 2 x T_n(x) - T_{n-1}(x) + # Thus + # U_{n+1}(x) == -(1-x^2)^(-1/2) T_n(x)/(n+1) + # == -(1-x^2)^(-1/2) * ( 2x T_{n-1}(x) - T_{n-2}(x))/(n+1) + # == -(1-x^2)^(-1/2) * ( 2n/(n+1) x T_{n-1}(x)/n - (n-1)/(n+1) T_{n-2}(x)/(n-1)) + # == (2n/(n+1) * x U_n(x) - (n - 1)/(n+1) * U_{n-1}(x) + + R = real(eltype(wP)) + n = zero(R):∞ + (2n) ./ (n .+ 1), Zeros{R}(∞), (n .- 1) ./ (n .+ 1) +end - # for U_0 - ζ = Vcat(convert(T,π)*(ξ^2/4 - (log.(abs.(ξ)) + log(2*one(T)))/2), ζ) - return transpose(ζ) - end +function complexlogkernel_ics(wP::Weighted{<:Any,<:ChebyshevU}, z::Number) + T = promote_type(eltype(wP), typeof(z)) + ξ = inv(z + sqrtx2(z)) + r0 = convert(T,π)*(ξ^2/4 - (log.(abs.(ξ)) + log(2*one(T)))/2) + r1 = convert(T,π)*(ξ^3/3 - ξ)/2 + r2 = convert(T,π)*(ξ^4/4 - ξ^2/2)/2 + r0,r1,r2 end -function complexlogkernel(P::Legendre, z::Number) +function complexlogkernel_recurrence(wP::Weighted{<:Any,<:ChebyshevU}) + # We have for n ≠ 0 L_n(z) = stieltjes(U_n, z) + # where U_n(z) = ∫_1^x sqrt(1-x^2) U_n(x) dx = -(1-x^2)^(3/2)C_{n-1}(x) * 2/(n * (n + 2)) + # where C_n(x) = C_n^{(2)}(x) + # We have the 3-term recurrence + # C_{n+1}(x) == 2(n + 2) / (n + 1) * x C_n(x) - (n + 3) / (n + 1) C_{n-1}(x) + # Thus + # U_{n+1}(x) == -(1-x^2)^(3/2) C_n(x) * 2/((n+1) * (n + 3)) + # == -(1-x^2)^(3/2) * 2/((n+1) * (n + 3)) ( 2(n + 1) / n * x C_{n-1}(x) - (n + 2) / n C_{n-2}(x)) + # == -(1-x^2)^(3/2) * (4 / (n*(n+3)) * x C_{n-1}(x) - 2 (n + 2) /(n*(n+1)*(n+3)) C_{n-2}(x)) + # == -(1-x^2)^(3/2) * (2 (n-1)*(n+1) / (n*(n+3)) * x 2/((n-1)*(n+1)) C_{n-1}(x) - (n + 2)*(n-1) /(n*(n+3)) * 2/((n-1)*(n+1)) C_{n-2}(x)) + # == (2 (n+2)/(n+3) * x U_n(x) - (n + 2)*(n-1) /(n*(n+3)) * U_{n-1}(x) + + R = real(eltype(wP)) + n = zero(R):∞ + (2*(n .+ 2)) ./ (n .+ 3), Zeros{R}(∞), (n .+ 2) .* (n .- 1) ./ (n .* (n .+ 3)) +end + +function complexlogkernel_ics(P::Weighted{<:Any,<:Legendre}, z::Number) T = promote_type(eltype(P), typeof(z)) r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T) r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1) r2 = z*r1 + 2*one(T)/3 - transpose(RecurrenceArray(z, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), [r0,r1,r2])) + r0,r1,r2 end -function complexlogkernel(P::Legendre, zs::AbstractVector) - T = promote_type(eltype(P), eltype(zs)) +function complexlogkernel_recurrence(wP::Weighted{<:Any,<:Legendre}) + # We have for n ≠ 0 L_n(z) = stieltjes(U_n, z) + # where U_n(z) = ∫_1^x P_n(x) dx = C_{n+1}^{(-1/2)}(x) + # Since these are equivalent to weihted OPs (1-x^2)C_{n-1}^(3/2)(x) + # we know they satisfy the same recurrence coefficients. + # Thus the following could also be written: + # A,B,C = recurrencecoefficients(Ultraspherical(-1/2)) + # A[2:end],B[2:end],C[2:end] + R = real(eltype(wP)) + ((one(R):2:∞)./(2:∞), Zeros{R}(∞), (-one(R):∞)./(2:∞)) +end +complexlogkernel(P::Legendre, z...) = complexlogkernel(Weighted(P), z...) +logkernel(P::Legendre, x...) = logkernel(Weighted(P), x...) + + + +complexlogkernel_layout(::WeightedOPLayout, wP, z::Number) = transpose(RecurrenceArray(z, complexlogkernel_recurrence(wP), [complexlogkernel_ics(wP,z)...])) + +function complexlogkernel_layout(::WeightedOPLayout, wP, zs::AbstractVector) + T = promote_type(eltype(wP), eltype(zs)) m = length(zs) data = Matrix{T}(undef, 3, m) for j = 1:m z = zs[j] - r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T) - r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1) - r2 = z*r1 + 2*one(T)/3 - data[1,j] = r0; data[2,j] = r1; data[3,j] = r2; + data[1:3,j] .= complexlogkernel_ics(wP, z) end - transpose(RecurrenceArray(zs, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), data)) + transpose(RecurrenceArray(zs, complexlogkernel_recurrence(wP), data)) end -logkernel(P::Legendre, z) = real.(complexlogkernel(P, complex(z))) -function logkernel(P::Legendre, x::Real) - T = promote_type(eltype(P), typeof(x)) - z = complex(x) - r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T) - r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1) - r2 = z*r1 + 2*one(T)/3 - transpose(RecurrenceArray(x, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), [real(r0),real(r1),real(r2)])) -end +logkernel_layout(::AbstractBasisLayout, P, z...) = real.(complexlogkernel(P, z...)) -function logkernel(P::Legendre, x::AbstractVector{<:Real}) - T = promote_type(eltype(P), eltype(x)) - m = length(x) - data = Matrix{T}(undef, 3, m) - for j = 1:m - z = complex(x[j]) - r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T) - r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1) - r2 = z*r1 + 2*one(T)/3 - data[1,j] = real(r0); data[2,j] = real(r1); data[3,j] = real(r2); - end +function logkernel_layout(::WeightedOPLayout, P, x::Real) + L = transpose(complexlogkernel(P, complex(x))) + transpose(RecurrenceArray(x, (L.A, L.B, L.C), real.(L.data))) +end - transpose(RecurrenceArray(x, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), data)) +function logkernel_layout(::WeightedOPLayout, P, x::AbstractVector{<:Real}) + L = transpose(complexlogkernel(P, complex(x))) + transpose(RecurrenceArray(x, (L.A, L.B, L.C), real.(L.data))) end diff --git a/src/recurrence.jl b/src/recurrence.jl index 8bc7880..b718e48 100644 --- a/src/recurrence.jl +++ b/src/recurrence.jl @@ -24,6 +24,9 @@ const RecurrenceMatrix{T, Z<:AbstractVector, A<:AbstractVector, B<:AbstractVecto function RecurrenceArray(z::Number, (A,B,C), data::AbstractVector{T}) where T N = length(data) p0, p1 = initiateforwardrecurrence(N, A, B, C, z, one(z)) + if iszero(p1) + p1 = one(p1) # avoid degeneracy in recurrence. Probably needs more thought + end RecurrenceVector{T,typeof(A),typeof(B),typeof(C)}(z, A, B, C, data, size(data), T[p0], T[p1], T[]) end @@ -33,6 +36,9 @@ function RecurrenceArray(z::AbstractVector, (A,B,C), data::AbstractMatrix{T}) wh p1 = Vector{T}(undef, N) for j = axes(z,1) p0[j], p1[j] = initiateforwardrecurrence(M, A, B, C, z[j], one(T)) + if iszero(p1[j]) + p1[j] = one(p1[j]) # avoid degeneracy in recurrence. Probably needs more thought + end end RecurrenceMatrix{T,typeof(z),typeof(A),typeof(B),typeof(C)}(z, A, B, C, data, size(data), p0, p1, T[]) end @@ -165,4 +171,9 @@ _getindex_resize_iffinite!(A, kr, jr, _) = layout_getindex(A, kr, jr) function view(A::RecurrenceVector, kr::AbstractVector) resizedata!(A, maximum(kr)) view(A.data, kr) -end \ No newline at end of file +end + +### +# broadcasted +### +broadcasted(::LazyArrayStyle, op, A::Transpose{<:Any,<:RecurrenceArray}) = transpose(op.(parent(A))) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 4f5163a..be536b6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,6 @@ using SingularIntegrals, ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, BandedMatrices, LinearAlgebra, Test using SingularIntegrals: Stieltjes, StieltjesPoint, ChebyshevInterval, associated, Associated, - orthogonalityweight, Weighted, *, dot, LogKernelPoint + orthogonalityweight, Weighted, *, dot using LazyArrays: MemoryLayout, PaddedLayout, colsupport, rowsupport, paddeddata using LazyBandedMatrices: blockcolsupport, Block, BlockHcat, blockbandwidths diff --git a/test/test_logkernel.jl b/test/test_logkernel.jl index 1876ee9..85a5861 100644 --- a/test/test_logkernel.jl +++ b/test/test_logkernel.jl @@ -1,5 +1,5 @@ using SingularIntegrals, ClassicalOrthogonalPolynomials, FillArrays, Test -using SingularIntegrals: RecurrenceArray +using SingularIntegrals: RecurrenceArray, LogKernelPoint using ClassicalOrthogonalPolynomials: affine @testset "ComplexLogKernelPoint" begin @@ -29,22 +29,36 @@ end wU = Weighted(ChebyshevU()) x = axes(wU,1) z = 0.1+0.2im - L = log.(abs.(z.-x')) + L = log.(abs.(z .- x')) @test L isa LogKernelPoint{Float64,ComplexF64,ComplexF64,Float64,ChebyshevInterval{Float64}} + @test (L * wU)[1:5] ≈ [ -1.2919202947616695, -0.20965486677056738, 0.6799687631764493, 0.13811497572177128, -0.2289481463304956] + + @test L * (wU / wU \ @.(exp(x) * sqrt(1 - x^2))) ≈ -1.4812979070884382 + + wT = Weighted(ChebyshevT()) + @test L * (wT / wT \ @.(exp(x) / sqrt(1 - x^2))) ≈ -1.9619040529776954 #mathematica end @testset "Real point" begin + T = ChebyshevT() U = ChebyshevU() x = axes(U,1) t = 2.0 + @test (log.(abs.(t .- x') )* Weighted(T))[1,1:3] ≈ [1.9597591637624774, -0.8417872144769223, -0.11277810215896047] #mathematica @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [1.0362686329607178,-0.4108206734393296, -0.054364775221816465] #mathematica t = 0.5 + @test (log.(abs.(t .- x') )* Weighted(T))[1,1:3] ≈ [-2.1775860903036017, -1.5707963267948832, 0.7853981633974272] #mathematica @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica t = 0.5+0im + @test (log.(abs.(t .- x') )* Weighted(T))[1,1:3] ≈ [-2.1775860903036017, -1.5707963267948832, 0.7853981633974272] #mathematica @test (log.(abs.(t .- x') )* Weighted(U))[1,1:3] ≈ [-1.4814921268505252, -1.308996938995747, 0.19634954084936207] #mathematica + + + # checks for degeneracy in recurrence initial conditions + @test iszero(logkernel(Weighted(T), 10.3)[100]) end @testset "mapped" begin