Add getindex for DunklXuDisk and tests

Author: dlfivefifty

src/rectdisk.jl

@@ -7,8 +7,10 @@ struct DunklXuDisk{T, V} <: BivariateOrthogonalPolynomial{T}
     β::V
 end
 
+DunklXuDisk{T}(β) where T = DunklXuDisk{T, typeof(β)}(β)
 DunklXuDisk(β::T) where T = DunklXuDisk{float(T), T}(β)
 DunklXuDisk() = DunklXuDisk(0)
+DunklXuDisk{T}() where T = DunklXuDisk{T}(0)
 
 ==(D1::DunklXuDisk, D2::DunklXuDisk) = D1.β == D2.β
 
@@ -18,6 +20,17 @@ copy(A::DunklXuDisk) = A
 show(io::IO, P::DunklXuDisk) = summary(io, P)
 summary(io::IO, P::DunklXuDisk) = print(io, "DunklXuDisk($(P.β))")
 
+function getindex(P::DunklXuDisk{T}, 𝐱::StaticVector{2}, JR::BlockOneTo) where T
+    x,y = 𝐱
+    n = Int(last(JR))
+    ret = zeros(T, n, n)
+    β = P.β
+    ρ = sqrt(1-x^2)
+    for j = 1:n
+        ret[1:n-j+1,j] = Jacobi{T}(j+β-1/2, j+β-1/2)[x,1:n-j+1] * ρ^(j-1) * Jacobi{T}(β, β)[y/ρ,j]
+    end
+    DiagTrav(ret)
+end
 
 """
     DunklXuDiskWeight(β)
@@ -37,6 +50,11 @@ axes(P::DunklXuDiskWeight{T}) where T = (Inclusion(UnitDisk{T}()),)
 show(io::IO, P::DunklXuDiskWeight) = summary(io, P)
 summary(io::IO, P::DunklXuDiskWeight) = print(io, "(1-x^2-y^2)^$(P.β) on the unit disk")
 
+function getindex(P::DunklXuDiskWeight, 𝐱::StaticVector{2})
+    r = norm(𝐱)
+    (1-r^2)^P.β
+end
+
 const WeightedDunklXuDisk{T} = WeightedBasis{T,<:DunklXuDiskWeight,<:DunklXuDisk}
 
 WeightedDunklXuDisk(β) = DunklXuDiskWeight(β) .* DunklXuDisk(β)

test/test_rectdisk.jl

@@ -1,6 +1,7 @@
 using MultivariateOrthogonalPolynomials, StaticArrays, BlockArrays, BlockBandedMatrices, ArrayLayouts, Base64,
         QuasiArrays, Test, ClassicalOrthogonalPolynomials, BandedMatrices, FastTransforms, LinearAlgebra
 import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, AngularMomentum
+using ForwardDiff
 
 @testset "Dunkl-Xu disk" begin
     @testset "basics" begin
@@ -13,6 +14,11 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang
         @test xy[SVector(0.1,0.2)] == SVector(0.1,0.2)
         @test x[SVector(0.1,0.2)] == 0.1
         @test y[SVector(0.1,0.2)] == 0.2
+
+        ρ = sqrt(1-0.1^2)
+        @test P[SVector(0.1,0.2),1] ≈ 1
+        @test P[SVector(0.1,0.2),Block(2)] ≈ [0.15,0.2]
+        @test P[SVector(0.1,0.2),Block(3)] ≈ [jacobip(2,1/2,1/2,0.1),jacobip(1,3/2,3/2,0.1)*ρ*legendrep(1,0.2/ρ),ρ^2*legendrep(2,0.2/ρ)]
     end
 
     @testset "operators" begin
@@ -29,26 +35,39 @@ import MultivariateOrthogonalPolynomials: dunklxu_raising, dunklxu_lowering, Ang
 
         x, y = coordinates(P)
 
-        L = WP \ WQ
-        R = Q \ P
+        @testset "lowering/raising" begin
+            L = WP \ WQ
+            @test WP[SVector(0.1,0.2),Block.(1:6)]'L[Block.(1:6),Block.(1:4)] ≈ WQ[SVector(0.1,0.2),Block.(1:4)]'
+            R = Q \ P
+            @test Q[SVector(0.1,0.2),Block.(1:4)]'R[Block.(1:4),Block.(1:4)] ≈ P[SVector(0.1,0.2),Block.(1:4)]'
 
-        X = P \ (x .* P)
-        Y = P \ (y .* P)
+            @test (DunklXuDisk() \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)] ≈ (WeightedDunklXuDisk(0.0) \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)]
+        end
 
-        @test (L * R)[Block.(1:N), Block.(1:N)] ≈ (I - X^2 - Y^2)[Block.(1:N), Block.(1:N)]
 
-        @test (DunklXuDisk() \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)] ≈ (WeightedDunklXuDisk(0.0) \ WeightedDunklXuDisk(1.0))[Block.(1:N), Block.(1:N)]
+        @testset "jacobi" begin
+            X = P \ (x .* P)
+            Y = P \ (y .* P)
 
-        ∂x = Derivative(P, (1,0))
-        ∂y = Derivative(P, (0,1))
+            @test (L * R)[Block.(1:N), Block.(1:N)] ≈ (I - X^2 - Y^2)[Block.(1:N), Block.(1:N)]
+            @test P[SVector(0.1,0.2),Block.(1:5)]'X[Block.(1:5),Block.(1:4)] ≈ 0.1P[SVector(0.1,0.2),Block.(1:4)]'
+            @test P[SVector(0.1,0.2),Block.(1:5)]'Y[Block.(1:5),Block.(1:4)] ≈ 0.2P[SVector(0.1,0.2),Block.(1:4)]'
+        end
 
-        Dx = Q \ (∂x * P)
-        Dy = Q \ (∂y * P)
+        @testset "derivatives" begin
+            ∂x = Derivative(P, (1,0))
+            ∂y = Derivative(P, (0,1))
 
-        Mx = Q \ (x .* Q)
-        My = Q \ (y .* Q)
+            Dx = Q \ (∂x * P)
+            Dy = Q \ (∂y * P)
+
+            @test Q[SVector(0.1,0.2),Block.(1:3)]'Dx[Block.(1:3),Block.(1:4)] ≈ [ForwardDiff.gradient(𝐱 -> DunklXuDisk{eltype(𝐱)}(P.β)[𝐱,k], SVector(0.1,0.2))[1] for k=1:10]'
+            Mx = Q \ (x .* Q)
+            My = Q \ (y .* Q)
+    
+            A = (Mx * Dy - My * Dx)[Block.(1:N), Block.(1:N)]
+        end
 
-        A = (Mx * Dy - My * Dx)[Block.(1:N), Block.(1:N)]
 
         B = (Q \ P)[Block.(1:N), Block.(1:N)]