Add DiskSlice

Project.toml

@@ -1,9 +1,11 @@
 name = "AnnuliOrthogonalPolynomials"
 uuid = "de1797fd-24c3-4035-91a2-b52aecdcfb01"
 authors = ["john.papad "]
-version = "0.1"
+version = "0.1.0"
 
 [deps]
+ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
 ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd"
@@ -12,14 +14,18 @@ DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 HarmonicOrthogonalPolynomials = "e416a80e-9640-42f3-8df8-80a93ca01ea5"
+InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MultivariateOrthogonalPolynomials = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SemiclassicalOrthogonalPolynomials = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
+ArrayLayouts = "1.11.1"
+BandedMatrices = "1.9.1"
 BlockArrays = "1.1"
 BlockBandedMatrices = "0.13"
 ClassicalOrthogonalPolynomials = "0.15"
@@ -28,17 +34,21 @@ DomainSets = "0.7"
 FastTransforms = "0.17"
 FillArrays = "1"
 HarmonicOrthogonalPolynomials = "0.7"
+InfiniteArrays = "0.15.5"
 LazyArrays = "2.1"
+LazyBandedMatrices = "0.11.4"
 MultivariateOrthogonalPolynomials = "0.9"
+QuadGK = "2"
 QuasiArrays = "0.12"
 SemiclassicalOrthogonalPolynomials = "0.7"
 StaticArrays = "1"
 julia = "1.10"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ForwardDiff", "Random", "Test"]
+test = ["ForwardDiff", "QuadGK", "Random", "Test"]

src/AnnuliOrthogonalPolynomials.jl

@@ -1,20 +1,27 @@
 module AnnuliOrthogonalPolynomials
-using BlockArrays, BlockBandedMatrices, ClassicalOrthogonalPolynomials, ContinuumArrays, DomainSets, FastTransforms, FillArrays,
-    HarmonicOrthogonalPolynomials, LazyArrays, LinearAlgebra,
-    MultivariateOrthogonalPolynomials, QuasiArrays, SemiclassicalOrthogonalPolynomials, StaticArrays 
+using BlockArrays, BandedMatrices, BlockBandedMatrices, ClassicalOrthogonalPolynomials, ContinuumArrays, DomainSets, FastTransforms, FillArrays,
+    HarmonicOrthogonalPolynomials, LazyArrays, LinearAlgebra, LazyBandedMatrices,
+    MultivariateOrthogonalPolynomials, QuasiArrays, SemiclassicalOrthogonalPolynomials, StaticArrays, ArrayLayouts, InfiniteArrays
 
-import Base: in, axes, getindex, broadcasted, tail, +, -, *, /, \, convert, OneTo, show, summary, ==, oneto, diff
-import BlockArrays: block, blockindex, _BlockedUnitRange, blockcolsupport
+
+import ArrayLayouts: MemoryLayout, sublayout, sub_materialize
+import Base: in, axes, getindex, broadcasted, tail, +, -, *, /, \, convert, OneTo, show, summary, ==, oneto, diff, copy, view, BroadcastStyle
+import BandedMatrices: _BandedMatrix
+import BlockArrays: block, blockindex, _BlockedUnitRange, blockcolsupport, BlockSlice
+import BlockBandedMatrices: blockbandwidths, subblockbandwidths, AbstractBandedBlockBandedLayout, AbstractBandedBlockBandedMatrix
 import ClassicalOrthogonalPolynomials: checkpoints, ShuffledR2HC, TransformFactorization, ldiv, paddeddata, jacobimatrix, orthogonalityweight, SetindexInterlace
-import ContinuumArrays: Weight, weight, grid, ℵ₁, ℵ₀, unweighted, plan_transform
+import ContinuumArrays: Weight, weight, grid, ℵ₁, ℵ₀, unweighted, plan_transform, @simplify
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan
 import MultivariateOrthogonalPolynomials: BlockOneTo, ModalInterlace, laplacian, ModalTrav
-import SemiclassicalOrthogonalPolynomials: divdiff, HalfWeighted
-
-export Block, SVector, CircleCoordinate, ZernikeAnnulus, ComplexZernikeAnnulus, Laplacian
+import SemiclassicalOrthogonalPolynomials: divdiff, HalfWeighted, SemiclassicalJacobiFamily
+import LazyBandedMatrices: AbstractLazyBandedBlockBandedLayout
+import LazyArrays: resizedata!
+using InfiniteArrays: InfiniteCardinal
 
-abstract type AnnuliOrthogonalPolynomial{d,T} <: MultivariateOrthogonalPolynomial{d,T} end
+export Block, SVector, CircleCoordinate, ZernikeAnnulus, ComplexZernikeAnnulus, Laplacian, JacobiDiskSlice
 
+include("columninterlace.jl")
 include("annulus.jl")
+include("diskslice.jl")
 
 end # module AnnuliOrthogonalPolynomials

src/annulus.jl

@@ -37,7 +37,8 @@ struct ZernikeAnnulus{T} <: AbstractZernikeAnnulus{T}
     ρ::T
     a::T
     b::T
-    ZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b)
+    P::SemiclassicalJacobiFamily{T}
+    ZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b, SemiclassicalJacobiFamily(inv(1-ρ^2),b,a,0:∞))
 end
 
 """
@@ -49,7 +50,8 @@ struct ComplexZernikeAnnulus{T} <: AbstractZernikeAnnulus{Complex{T}}
     ρ::T
     a::T
     b::T
-    ComplexZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b)
+    P::SemiclassicalJacobiFamily{T}
+    ComplexZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b, SemiclassicalJacobiFamily(inv(1-ρ^2),b,a,0:∞))
 end
 
 
@@ -71,7 +73,9 @@ copy(A::AbstractZernikeAnnulus) = A
 
 orthogonalityweight(Z::AbstractZernikeAnnulus) = AnnulusWeight(Z.ρ, Z.a, Z.b)
 
-zernikeannulusr(ρ, ℓ, m, a, b, r::T) where T = r^m * SemiclassicalJacobi{T}(inv(1-ρ^2),b,a,m)[(r^2 - 1)/(ρ^2 - 1), (ℓ-m) ÷ 2 + 1]
+zernikeannulusr(ρ, ℓ, m, a, b, r, P) = r^m * P[(r^2 - 1)/(ρ^2 - 1), (ℓ-m) ÷ 2 + 1]
+
+zernikeannulusr(ρ, ℓ, m, a, b, r::T) where T = zernikeannulusr(ρ, ℓ, m, a, b, r, SemiclassicalJacobi{T}(inv(1-ρ^2),b,a,m))
 function zernikeannulusz(ρ, ℓ, ms, a, b, rθ::RadialCoordinate{T}) where T
     r,θ = rθ.r,rθ.θ
     m = abs(ms)
@@ -97,15 +101,19 @@ function getindex(Z::ZernikeAnnulus{T}, rθ::RadialCoordinate, B::BlockIndex{1})
     ℓ = Int(block(B))-1
     k = blockindex(B)
     m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
-    zernikeannulusz(Z.ρ, ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
+    (; r, θ) = rθ
+    (; ρ, a, b, P) = Z
+    (isodd(k+ℓ) ? cos(m*θ) : sin(m*θ)) * zernikeannulusr(ρ, ℓ, m, a, b, r, P[m+1])
 end
 
 
 function getindex(Z::ComplexZernikeAnnulus{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
     ℓ = Int(block(B))-1
     k = blockindex(B)
     m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
-    complexzernikeannulusz(Z.ρ, ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
+    (; r, θ) = rθ
+    (; ρ, a, b, P) = Z
+    exp(im*(isodd(k+ℓ) ? m : -m)*θ) * zernikeannulusr(ρ, ℓ, m, a, b, r, P[m+1])
 end
 
 

src/columninterlace.jl

@@ -0,0 +1,80 @@
+"""
+    ColumnInterlace(ops, (M,N), (l,u))
+
+interlaces the entries of a vector of banded matrices
+acting on different columns of the matrix underling a `DiagTrav`. 
+"""
+struct ColumnInterlace{T, MMNN<:Tuple} <: AbstractBandedBlockBandedMatrix{T}
+    ops
+    MN::MMNN
+    bandwidths::NTuple{2,Int}
+end
+
+ColumnInterlace{T}(ops, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) where T = ColumnInterlace{T,typeof(MN)}(ops, MN, bandwidths)
+ColumnInterlace(ops::AbstractVector{<:AbstractMatrix}, MN::NTuple{2,Integer}, bandwidths::NTuple{2,Int}) = ColumnInterlace{eltype(eltype(ops))}(ops, MN, bandwidths)
+
+axes(Z::ColumnInterlace) = blockedrange.(oneto.(Z.MN))
+
+blockbandwidths(R::ColumnInterlace) = R.bandwidths
+subblockbandwidths(::ColumnInterlace) = (0,0)
+copy(M::ColumnInterlace) = M
+
+
+function view(R::ColumnInterlace{T}, KJ::Block{2}) where T
+    K,J = KJ.n
+    dat = Matrix{T}(undef,1,J)
+    l,u = blockbandwidths(R)
+    if -l ≤ J - K ≤ u
+        sh = (J-K)
+        for j in 1:min(K,J)
+            dat[1,j] = R.ops[j][K-j+1,J-j+1]
+        end
+    else
+        fill!(dat, zero(T))
+    end
+    _BandedMatrix(dat, K, 0, 0)
+end
+
+getindex(R::ColumnInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
+
+struct ColumnInterlaceLayout <: AbstractBandedBlockBandedLayout end
+struct LazyColumnInterlaceLayout <: AbstractLazyBandedBlockBandedLayout end
+
+MemoryLayout(::Type{<:ColumnInterlace}) = ColumnInterlaceLayout()
+MemoryLayout(::Type{<:ColumnInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyColumnInterlaceLayout()
+sublayout(::Union{ColumnInterlaceLayout,LazyColumnInterlaceLayout}, ::Type{<:NTuple{2,BlockSlice{<:BlockOneTo}}}) = ColumnInterlaceLayout()
+
+
+function sub_materialize(::ColumnInterlaceLayout, V::AbstractMatrix{T}) where T
+    kr,jr = parentindices(V)
+    KR,JR = kr.block,jr.block
+    M,N = Int(last(KR)), Int(last(JR))
+    R = parent(V)
+    ColumnInterlace{T}([layout_getindex(R.ops[k],1:(M-k+1),1:(N-k+1)) for k=1:min(N,M)], (M,N), R.bandwidths)
+end
+
+# act like lazy array
+BroadcastStyle(::Type{<:ColumnInterlace{<:Any,NTuple{2,InfiniteCardinal{0}}}}) = LazyArrayStyle{2}()
+
+# TODO: overload muladd!
+function *(A::ColumnInterlace, b::DiagTrav)
+    M = b.array
+    ret = zeros(promote_type(eltype(A), eltype(b)), size(M)...)
+    m = maximum(colsupport(M))
+    n = maximum(rowsupport(M))
+    resizedata!(ret, m, n)
+    for j = 1:n
+        mul!(view(ret,1:m-j+1,j), view(A.ops[j], 1:m-j+1, 1:m-j+1), view(M,1:m-j+1,j))
+    end
+    DiagTrav(ret)
+end
+
+
+function \(A::ColumnInterlace, b::DiagTrav)
+    M = b.matrix
+    ret = similar(M, promote_type(eltype(A),eltype(b)))
+    for j = 1:size(ret,2)
+        ldiv!(@view(ret[1:end-j+1,j]), A.ops[j], @view(M[1:end-j+1,j]))
+    end
+    DiagTrav(ret)
+end

src/diskslice.jl

@@ -0,0 +1,92 @@
+"""
+    DiskSlice(α)
+
+represents α ≤ x ≤ 1, y^2 ≤ 1-x^2.
+"""
+struct DiskSlice{T} <: Domain{T}
+    α::T
+end
+
+
+function ClassicalOrthogonalPolynomials.checkpoints(d::DiskSlice)
+    α = d.α
+    x = α + (1-α)/3
+    [SVector(x, sqrt(1-x^2)/7)]
+end
+
+"""
+    DiskSliceWeight(α, a, b)
+
+is a quasi-vector representing `(x-α)^a * (1-x^2-y^2)^b`
+"""
+struct DiskSliceWeight{T} <: Weight{T}
+    α::T
+    a::T
+    b::T
+end
+
+copy(w::DiskSliceWeight) = w
+
+axes(w::DiskSliceWeight{T}) where T = (Inclusion(DiskSlice(w.α)),)
+
+==(w::DiskSliceWeight, v::DiskSliceWeight) = w.a == v.a && w.b == v.b && w.α == v.α
+
+function getindex(w::DiskSliceWeight, 𝐱::StaticVector{2})
+    x,y = 𝐱
+    (x-w.α)^w.a * (1-x^2-y^2)^w.b
+end
+
+"""
+    JacobiDiskSlice(α, a, b)
+
+is a quasi-matrix orthogonal to `(x-α)^a * (1-x^2-y^2)^b`.
+"""
+struct JacobiDiskSlice{T} <: BivariateOrthogonalPolynomial{T} 
+    α::T
+    a::T
+    b::T
+    P::SemiclassicalJacobiFamily{T}
+    JacobiDiskSlice{T}(α::T, a::T, b::T) where T = new{T}(α, a, b, SemiclassicalJacobiFamily(2/(1-α), (b+one(T)/2):∞, a, (b+one(T)/2):∞))
+end
+
+JacobiDiskSlice{T}(α, a, b) where T = JacobiDiskSlice{T}(convert(T,α), convert(T,a), convert(T,b))
+JacobiDiskSlice(α::R, a::T, b::V) where {R,T,V} = JacobiDiskSlice{float(promote_type(R,T,V))}(α, a, b)
+JacobiDiskSlice{T}(α) where T = JacobiDiskSlice{T}(α, zero(α), zero(α))
+JacobiDiskSlice(α) = JacobiDiskSlice(α, zero(α), zero(α))
+
+convert(::Type{AbstractQuasiMatrix{T}}, P::JacobiDiskSlice) where T = JacobiDiskSlice{T}(P.α, P.a, P.b)
+
+==(w::JacobiDiskSlice, v::JacobiDiskSlice) = w.α == v.α && w.a == v.a && w.b == v.b
+
+axes(P::JacobiDiskSlice{T}) where T = (Inclusion(DiskSlice(P.α)),blockedrange(oneto(∞)))
+copy(A::JacobiDiskSlice) = A
+
+orthogonalityweight(P::JacobiDiskSlice) = DiskSliceWeight(P.α, P.a, P.b)
+
+
+function getindex(P::JacobiDiskSlice{T}, 𝐱::StaticVector{2}, B::BlockIndex{1}) where T
+    n,k = Int(block(B)), blockindex(B)
+    x,y = 𝐱
+    ρ = sqrt(1-x^2)
+    P.P[k][(x-1)/(P.α-1),n-k+1] * ρ^(k-1) * ultrasphericalc(k-1, P.b+one(T)/2, y/ρ)
+end
+getindex(P::JacobiDiskSlice, 𝐱::StaticVector{2}, B::Block{1}) = [P[𝐱, B[j]] for j=1:Int(B)]
+getindex(P::JacobiDiskSlice, 𝐱::StaticVector{2}, JR::BlockOneTo) = mortar([P[𝐱,Block(J)] for J = 1:Int(JR[end])])
+
+
+###
+# jacobimatrix
+###
+
+jacobimatrix(::Val{1}, P::JacobiDiskSlice) = ColumnInterlace(BroadcastVector{AbstractMatrix{Float64}}((α,Pₖ) -> (α-1)*jacobimatrix(Pₖ) + I, P.α, P.P), (ℵ₀,ℵ₀), (1,1))
+
+
+###
+# raising
+####
+
+@simplify function \(Q::JacobiDiskSlice, P::JacobiDiskSlice)
+    T = promote_type(eltype(Q), eltype(P))
+    @assert Q.a == P.a + 1 && Q.b == P.b
+    ColumnInterlace(BroadcastVector{AbstractMatrix{T}}(\, Q.P, P.P), (ℵ₀,ℵ₀), (0,1))
+end

test/runtests.jl

@@ -3,4 +3,5 @@
 # using AlgebraicCurveOrthogonalPolynomials, LinearAlgebra, BlockBandedMatrices, BlockArrays, FillArrays, Test
 using AnnuliOrthogonalPolynomials, Test
 
-include("test_annulus.jl")
+include("test_annulus.jl")
+include("test_diskslice.jl")

test/test_annulus.jl

@@ -223,7 +223,7 @@ import LazyArrays: Ones
         B = ComplexZernikeAnnulus(ρ,0,1)
 
         xy = SVector(0.5,0.1)
-        @test A[xy, 1:3] ≈ [1.0 + 0.0im;0.5 + 0.10000000000000002im;0.5 + 0.10000000000000002im]
+        @test A[xy, 1:3] ≈ [1.0 + 0.0im;0.5 - 0.10000000000000002im;0.5 + 0.10000000000000002im]
         R = B \ A
         @test A[SVector(0.5,0.1), Block.(1:3)]' ≈ B[SVector(0.5,0.1), Block.(1:3)]' * R[Block.(1:3),Block.(1:3)]
 

test/test_columninterlace.jl

@@ -0,0 +1,19 @@
+using AnnuliOrthogonalPolynomials, ClassicalOrthogonalPolynomials, InfiniteArrays, LazyBandedMatrices, Test
+using AnnuliOrthogonalPolynomials: ColumnInterlace
+
+
+@testset "ColumnInterlace" begin
+    X = ColumnInterlace(jacobimatrix.(Jacobi.(0,1:∞)), (ℵ₀,ℵ₀), (1,1))
+    Xₙ = X[Block.(Base.oneto(5)), Block.(Base.oneto(5))]
+    @test Xₙ isa ColumnInterlace
+    Cₙ = DiagTrav(zeros(5,5))
+    Cₙ[1:9] = 1:9
+    @test Xₙ*Cₙ ≈ X[Block.(1:5), Block.(1:5)] * Vector(Cₙ)
+
+    for j = 1:5
+        @test (Xₙ*Cₙ).array[:,j] ≈ X.ops[j][1:5,1:5] * Cₙ.array[:,j]
+    end
+    C = DiagTrav(zeros(∞,∞))
+    C[1:9] = 1:9
+    @test (X*C)[Block.(1:5)] ≈ Xₙ*Cₙ
+end