work on grid for Zernike (#135)

* work on grid for Zernike

* tests pass

Author: dlfivefifty

Project.toml

@@ -1,6 +1,6 @@
 name = "MultivariateOrthogonalPolynomials"
 uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
-version = "0.4"
+version = "0.4.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -28,12 +28,12 @@ ArrayLayouts = "0.8.6"
 BandedMatrices = "0.17"
 BlockArrays = "0.16.14"
 BlockBandedMatrices = "0.11.5"
-ClassicalOrthogonalPolynomials = "0.6.9"
-ContinuumArrays = "0.11"
+ClassicalOrthogonalPolynomials = "0.7"
+ContinuumArrays = "0.12"
 DomainSets = "0.5, 0.6"
 FastTransforms = "0.14"
 FillArrays = "0.13"
-HarmonicOrthogonalPolynomials = "0.3"
+HarmonicOrthogonalPolynomials = "0.4"
 InfiniteArrays = "0.12"
 InfiniteLinearAlgebra = "0.6.6"
 LazyArrays = "0.22.10"

src/MultivariateOrthogonalPolynomials.jl

@@ -24,7 +24,7 @@ import InfiniteArrays: InfiniteCardinal, OneToInf
 import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw
 import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan,
                                           PartialDerivative, AngularMomentum, BlockOneTo, BlockRange1, interlace,
-                                          MultivariateOPLayout
+                                          MultivariateOPLayout, MAX_PLOT_BLOCKS
 
 export MultivariateOrthogonalPolynomial, BivariateOrthogonalPolynomial,
        UnitTriangle, UnitDisk,

src/disk.jl

@@ -191,11 +191,9 @@ end
 # Transforms
 ###
 
-const FiniteZernike{T} = SubQuasiArray{T,2,Zernike{T},<:Tuple{<:Inclusion,<:BlockSlice{BlockRange1{OneTo{Int}}}}}
-
-function grid(S::FiniteZernike{T}) where T
-    N = blocksize(S,2) ÷ 2 + 1 # polynomial degree
-    M = 4N-3
+function grid(S::Zernike{T}, B::Block{1}) where T
+    N = Int(B) ÷ 2 + 1 # matrix rows
+    M = 4N-3 # matrix columns
 
     r = sinpi.((N .-(0:N-1) .- one(T)/2) ./ (2N))
 
@@ -206,26 +204,21 @@ end
 
 _angle(rθ::RadialCoordinate) = rθ.θ
 
-function plotgrid(S::FiniteZernike{T}) where T
-    N = blocksize(S,2) ÷ 2 + 1 # polynomial degree
-    g = grid(parent(S)[:,Block.(OneTo(2N))]) # double sampling
+function plotgrid(S::Zernike{T}, B::Block{1}) where T
+    N = Int(B) ÷ 2 + 1 # polynomial degree
+    g = grid(S, Block(min(2N, MAX_PLOT_BLOCKS))) # double sampling
     θ = [map(_angle,g[1,:]); 0]
-    [permutedims(RadialCoordinate.(1,θ)); g g[:,1]; permutedims(RadialCoordinate.(0,θ))]
-end
-
-function plotgrid(S::SubQuasiArray{<:Any,2,<:Zernike})
-    kr,jr = parentindices(S)
-    Z = parent(S)
-    plotgrid(Z[kr,Block.(OneTo(Int(findblock(axes(Z,2),maximum(jr)))))])
+    [permutedims(RadialCoordinate.(1,θ));
+     g g[:,1];
+     permutedims(RadialCoordinate.(0,θ))]
 end
 
-
 function plotvalues(u::ApplyQuasiVector{T,typeof(*),<:Tuple{Zernike, AbstractVector}}, x) where T
     Z,c = u.args
-    CS = blockcolsupport(c)
-    N = Int(last(CS)) ÷ 2 + 1 # polynomial degree
-    F = ZernikeITransform{T}(2N, Z.a, Z.b)
-    C = F * c[Block.(OneTo(2N))] # transform to grid
+    B = findblock(axes(Z,2), last(colsupport(c)))
+    N = Int(B) ÷ 2 + 1 # polynomial degree
+    F = ZernikeITransform{T}(min(2N, MAX_PLOT_BLOCKS), Z.a, Z.b)
+    C = F * c[Block.(OneTo(min(2N, MAX_PLOT_BLOCKS)))] # transform to grid
     [permutedims(u[x[1,:]]); # evaluate on edge of disk
      C C[:,1];
      fill(u[x[end,1]], 1, size(x,2))] # evaluate at origin and repeat
@@ -262,8 +255,11 @@ end
 *(P::ZernikeTransform{T}, f::Matrix{T}) where T = ModalTrav(P.disk2cxf \ (P.analysis * f))
 *(P::ZernikeITransform, f::AbstractVector) = P.synthesis * (P.disk2cxf * ModalTrav(f).matrix)
 
-factorize(S::FiniteZernike{T}) where T = TransformFactorization(grid(S), ZernikeTransform{T}(blocksize(S,2), parent(S).a, parent(S).b))
+plan_grid_transform(Z::Zernike{T}, B::Tuple{Block{1}}, dims=1:1) where T = grid(Z,B[1]), ZernikeTransform{T}(Int(B[1]), Z.a, Z.b)
 
+##
+# Laplacian
+###
 
 @simplify function *(Δ::Laplacian, WZ::Weighted{<:Any,<:Zernike})
     @assert WZ.P.a == 0 && WZ.P.b == 1

src/rect.jl

@@ -59,11 +59,11 @@ function checkpoints(P::RectPolynomial)
     SVector.(x, y')
 end
 
-function plan_grid_transform(P::KronPolynomial{d,<:Any,<:Fill}, B, dims=1:ndims(B)) where d
+function plan_grid_transform(P::KronPolynomial{d,<:Any,<:Fill}, B::Tuple{Block{1}}, dims=1:1) where d
     @assert dims == 1
 
     T = first(P.args)
-    x, F = plan_grid_transform(T, Array{eltype(B)}(undef, Fill(blocksize(B,1),d)...))
+    x, F = plan_grid_transform(T, Array{eltype(P)}(undef, Fill(Int(B[1]),d)...))
     @assert d == 2
     x̃ = Vector(x)
     SVector.(x̃, x̃'), ApplyPlan(DiagTrav, F)

src/triangle.jl

@@ -637,19 +637,13 @@ function tridenormalize!(F̌,a,b,c)
     F̌
 end
 
-function _trigrid(N::Integer)
+function trigrid(N::Integer)
     M = N
     x = [sinpi((2N-2n-1)/(4N))^2 for n in 0:N-1]
     w = [sinpi((2M-2m-1)/(4M))^2 for m in 0:M-1]
     [SVector(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]
-end
-trigrid(N::Block{1}) = _trigrid(Int(N))
-
-
-function grid(Pn::SubQuasiArray{T,2,<:JacobiTriangle,<:Tuple{Inclusion,BlockSlice{<:BlockRange{1}}}}) where T
-    kr,jr = parentindices(Pn)
-    trigrid(maximum(jr.block))
-end
+end 
+grid(P::JacobiTriangle, B::Block{1}) = trigrid(Int(B))
 
 struct TriPlan{T}
     tri2cheb::FastTransforms.FTPlan{T,2,FastTransforms.TRIANGLE}
@@ -666,8 +660,8 @@ TriPlan{T}(N::Block{1}, a, b, c) where T = TriPlan{T}(Matrix{T}(undef, Int(N), I
 
 *(T::TriPlan, F::AbstractMatrix) = DiagTrav(tridenormalize!(T.tri2cheb\(T.grid2cheb*F),T.a,T.b,T.c))
 
-function plan_grid_transform(P::JacobiTriangle, arr::AbstractVector, dims=1:ndims(arr))
-    T = promote_type(eltype(P), eltype(arr))
-    N = findblock(axes(P,2), length(arr))
-    grid(P[:,Block.(OneTo(Int(N)))]), TriPlan{T}(N, P.a, P.b, P.c)
+function plan_grid_transform(P::JacobiTriangle, Bs::Tuple{Block{1}}, dims=1:1)
+    T = eltype(P)
+    N = Bs[1]
+    grid(P, N), TriPlan{T}(N, P.a, P.b, P.c)
 end

test/test_disk.jl

@@ -309,7 +309,7 @@ import ForwardDiff: hessian
     @testset "plotting" begin
         Z = Zernike()
         u = Z * [1; 2; zeros(∞)];
-        rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u)
+        rep = RecipesBase.apply_recipe(Dict{Symbol, Any}(), u);
         g = MultivariateOrthogonalPolynomials.plotgrid(Z[:,1:3])
         @test all(rep[1].args .≈ (first.(g),last.(g),u[g]))