Lowering for Zernike polynomials (#68)

* Generalize ModalInterlace

* Zernike lowering

* Update test_disk.jl

* Update Project.toml

* add tests

Author: dlfivefifty

Project.toml

@@ -28,8 +28,8 @@ ArrayLayouts = "0.6"
 BandedMatrices = "0.16"
 BlockArrays = "0.14.1, 0.15"
 BlockBandedMatrices = "0.10.1"
-ClassicalOrthogonalPolynomials = "0.3.1"
-ContinuumArrays = "0.6"
+ClassicalOrthogonalPolynomials = "0.3.2"
+ContinuumArrays = "0.6.3"
 DomainSets = "0.4"
 FastTransforms = "0.11, 0.12"
 FillArrays = "0.11"

src/MultivariateOrthogonalPolynomials.jl

@@ -15,7 +15,7 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, AbstractBandedBlockBandedM
 import LinearAlgebra: factorize
 import LazyArrays: arguments, paddeddata
 
-import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight
+import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
                                           PartialDerivative, BlockOneTo, BlockRange1, interlace
 

src/disk.jl

@@ -53,13 +53,19 @@ struct ZernikeWeight{T} <: Weight{T}
     b::T
 end
 
+
 """
     ZernikeWeight(b)
 
 is a quasi-vector representing `(1-r^2)^b`
 """
 
 ZernikeWeight(b) = ZernikeWeight(zero(b), b)
+ZernikeWeight{T}(b) where T = ZernikeWeight{T}(zero(T), b)
+ZernikeWeight{T}() where T = ZernikeWeight{T}(zero(T))
+ZernikeWeight() = ZernikeWeight{Float64}()
+
+copy(w::ZernikeWeight) = w
 
 axes(::ZernikeWeight{T}) where T = (Inclusion(UnitDisk{T}()),)
 
@@ -186,50 +192,55 @@ getindex(W::WeightedZernikeLaplacianDiag, k::Integer) = W[findblockindex(axes(W,
     WZ.P * Diagonal(WeightedZernikeLaplacianDiag{eltype(eltype(WZ))}())
 end
 
-struct ZernikeConversion{T} <: AbstractBandedBlockBandedMatrix{T} end
 
-axes(Z::ZernikeConversion) = (blockedrange(oneto(∞)), blockedrange(oneto(∞)))
+"""
+    ModalInterlace
+"""
+struct ModalInterlace{T} <: AbstractBandedBlockBandedMatrix{T}
+    ops
+    bandwidths::NTuple{2,Int}
+end
+
+axes(Z::ModalInterlace) = (blockedrange(oneto(∞)), blockedrange(oneto(∞)))
 
-blockbandwidths(::ZernikeConversion) = (0,2)
-subblockbandwidths(::ZernikeConversion) = (0,0)
+blockbandwidths(R::ModalInterlace) = R.bandwidths
+subblockbandwidths(::ModalInterlace) = (0,0)
 
 
-function Base.view(W::ZernikeConversion{T}, KJ::Block{2}) where T
+function Base.view(R::ModalInterlace{T}, KJ::Block{2}) where T
     K,J = KJ.n
     dat = Matrix{T}(undef,1,J)
-    if J == K
+    l,u = blockbandwidths(R)
+    if iseven(J-K) && -l ≤ J - K ≤ u
+        sh = (J-K)÷2
         if isodd(K)
-            R0 = Normalized(Jacobi(1,0)) \ Normalized(Jacobi(0,0))
-            dat[1,1] = R0[K÷2+1,K÷2+1]
+            k = K÷2+1
+            dat[1,1] = R.ops[1][k,k+sh]
         end
-        for m in range(2-iseven(K); step=2, length=J÷2)
-            Rm = Normalized(Jacobi(1,m)) \ Normalized(Jacobi(0,m))
-            j = K÷2-m÷2+isodd(K)
-            dat[1,m] = dat[1,m+1] = Rm[j,j]
-        end
-    elseif J == K + 2
-        if isodd(K)
-            R0 = Normalized(Jacobi(1,0)) \ Normalized(Jacobi(0,0))
-            j = K÷2+1
-            dat[1,1] = R0[j,j+1]
-        end
-        for m in range(2-iseven(K); step=2, length=K÷2)
-            Rm = Normalized(Jacobi(1,m)) \ Normalized(Jacobi(0,m))
-            j = K÷2-m÷2+isodd(K)
-            dat[1,m] = dat[1,m+1] = Rm[j,j+1]
+        for m in range(2-iseven(K); step=2, length=J÷2-max(0,sh))
+            k = K÷2-m÷2+isodd(K)
+            dat[1,m] = dat[1,m+1] = R.ops[m+1][k,k+sh]
         end
     else
         fill!(dat, zero(T))
     end
-    dat ./= sqrt(2one(T))
     _BandedMatrix(dat, K, 0, 0)
 end
 
-getindex(R::ZernikeConversion, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
+getindex(R::ModalInterlace, k::Integer, j::Integer) = R[findblockindex.(axes(R),(k,j))...]
 
 function \(A::Zernike{T}, B::Zernike{V}) where {T,V}
-    A.a == B.a && A.b == B.b && return Eye{promote_type(T,V)}(∞)
+    TV = promote_type(T,V)
+    A.a == B.a && A.b == B.b && return Eye{TV}(∞)
     @assert A.a == 0 && A.b == 1
     @assert B.a == 0 && B.b == 0
-    ZernikeConversion{promote_type(T,V)}()
+    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(1,0:∞)) .\ Normalized.(Jacobi{TV}.(0,0:∞))) ./ sqrt(convert(TV, 2)), (0,2))
+end
+
+function \(A::Zernike{T}, B::Weighted{V,Zernike{V}}) where {T,V}
+    TV = promote_type(T,V)
+    A.a == B.P.a == A.b == B.P.b == 0 && return Eye{TV}(∞)
+    @assert A.a == A.b == 0
+    @assert B.P.a == 0 && B.P.b == 1
+    ModalInterlace{TV}((Normalized.(Jacobi{TV}.(0, 0:∞)) .\ HalfWeighted{:a}.(Normalized.(Jacobi{TV}.(1, 0:∞)))) ./ sqrt(convert(TV, 2)), (2,0))
 end

test/test_disk.jl

@@ -1,5 +1,5 @@
 using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, Test
-import MultivariateOrthogonalPolynomials: DiskTrav, grid, ZernikeConversion
+import MultivariateOrthogonalPolynomials: DiskTrav, grid
 import ClassicalOrthogonalPolynomials: HalfWeighted
 
 
@@ -10,6 +10,13 @@ import ClassicalOrthogonalPolynomials: HalfWeighted
 
         @test Zernike() == Zernike()
         @test Zernike(1) ≠ Zernike()
+        @test Zernike() ≡ copy(Zernike())
+
+        @test ZernikeWeight() == ZernikeWeight() == ZernikeWeight(0,0) == 
+                ZernikeWeight(0) == ZernikeWeight{Float64}() == 
+                ZernikeWeight{Float64}(0) == ZernikeWeight{Float64}(0, 0)
+        @test ZernikeWeight(1) ≠ ZernikeWeight()
+        @test ZernikeWeight() ≡ copy(ZernikeWeight())
     end
 
     @testset "Evaluation" begin
@@ -176,6 +183,29 @@ import ClassicalOrthogonalPolynomials: HalfWeighted
 
 
         R = Zernike(1) \ Zernike()
-        @test Zernike()[xy,Block.(1:5)]' ≈ Zernike(1)[xy,Block.(1:5)]'*R[Block.(1:5),Block.(1:5)] 
+        @test Zernike()[xy,Block.(1:6)]' ≈ Zernike(1)[xy,Block.(1:6)]'*R[Block.(1:6),Block.(1:6)] 
+    end
+
+    @testset "Lowering" begin
+        L0 = Normalized(Jacobi(0, 0)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 0)))
+        L1 = Normalized(Jacobi(0, 1)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 1)))
+        L2 = Normalized(Jacobi(0, 2)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 2)))
+        L3 = Normalized(Jacobi(0, 3)) \ HalfWeighted{:a}(Normalized(Jacobi(1, 3)))
+
+        xy = SVector(0.1,0.2)
+        r = norm(xy)
+        w = 1 - r^2
+
+        @test w*Zernike(1)[xy,Block(1)[1]] ≈ L0[1:2,1]'*Zernike()[xy,getindex.(Block.(1:2:3),1)] / sqrt(2)
+
+        @test w*Zernike(1)[xy,Block(2)[1]] ≈ L1[1:2,1]'*Zernike()[xy,getindex.(Block.(2:2:4),1)]/sqrt(2)
+        @test w*Zernike(1)[xy,Block(2)[2]] ≈ L1[1:2,1]'*Zernike()[xy,getindex.(Block.(2:2:4),2)]/sqrt(2)
+
+        @test w*Zernike(1)[xy,Block(3)[1]] ≈ L0[2:3,2]'*Zernike()[xy,getindex.(Block.(3:2:5),1)]/sqrt(2)
+        @test w*Zernike(1)[xy,Block(3)[2]] ≈ L2[1:2,1]'*Zernike()[xy,getindex.(Block.(3:2:5),2)]/sqrt(2)
+        @test w*Zernike(1)[xy,Block(3)[3]] ≈ L2[1:2,1]'*Zernike()[xy,getindex.(Block.(3:2:5),3)]/sqrt(2)
+
+        L = Zernike() \ Weighted(Zernike(1))
+        @test w*Zernike(1)[xy,Block.(1:5)]' ≈ Zernike()[xy,Block.(1:7)]'*L[Block.(1:7),Block.(1:5)] 
     end
 end