Merge pull request #9 from ioannisPApapadopoulos/jp/via_jacobi

Jp/via jacobi

Manifest.toml

@@ -20,9 +20,7 @@ version = "1.5.0"
 
 [[deps.AnnuliOrthogonalPolynomials]]
 deps = ["BlockArrays", "BlockBandedMatrices", "ClassicalOrthogonalPolynomials", "ContinuumArrays", "DomainSets", "FastTransforms", "FillArrays", "HarmonicOrthogonalPolynomials", "LazyArrays", "LinearAlgebra", "MultivariateOrthogonalPolynomials", "QuasiArrays", "SemiclassicalOrthogonalPolynomials", "StaticArrays"]
-git-tree-sha1 = "bed642acc470e3e06f83673a145894f752495205"
-repo-rev = "main"
-repo-url = "https://github.com/JuliaApproximation/AnnuliOrthogonalPolynomials.jl.git"
+path = "C:\\Users\\john.papad\\.julia\\dev\\AnnuliOrthogonalPolynomials"
 uuid = "de1797fd-24c3-4035-91a2-b52aecdcfb01"
 version = "0.0.1"
 

Project.toml

@@ -24,19 +24,19 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-AnnuliOrthogonalPolynomials = "0.0.1"
+AnnuliOrthogonalPolynomials = "0.0.1, 0.0.2"
 BandedMatrices = "0.17"
 BlockArrays = "0.16"
 BlockBandedMatrices = "0.12"
 ClassicalOrthogonalPolynomials = "0.11"
-ContinuumArrays = "0.12, 0.13, 0.14"
+ContinuumArrays = "0.12, 0.13, 0.14, 0.16"
 DomainSets = "0.6"
 FastTransforms = "0.15"
 FillArrays = "1.5"
 HypergeometricFunctions = "0.3"
 LazyArrays = "1.5"
 Memoization = "0.2"
-MultivariateOrthogonalPolynomials = "0.5"
+MultivariateOrthogonalPolynomials = "0.5, 0.6"
 QuasiArrays = "0.9, 0.10, 0.11"
 SemiclassicalOrthogonalPolynomials = "0.5"
 SpecialFunctions = "2"

src/annuluselement.jl

@@ -15,11 +15,11 @@ annulus(ρ::T, r::T) where T = (r*UnitDisk{T}()) \ (ρ*UnitDisk{T}())
 # and cannot be used for L2 inner-product of FiniteZernikeBasis
 # and FiniteContinuousZernike
 function _ann2element_via_Jacobi(t::T, m::Int) where T
-    Q₁₁ = SemiclassicalJacobi{T}(t, 1, 1, m)
+    Q₁₁ = SemiclassicalJacobi(t, 1, 1, m)
 
     x = axes(Q₁₁, 1)
-    X = Q₁₁ \ (x.*Q₁₁)
-    L₁₁ = (X - ApplyArray(*,X,X))/t^2
+    X = jacobimatrix(Q₁₁)
+    L₁₁ = (X - X*X)/t^2
     L₀₁ = (I - X)/t
     L₁₀ = X/t
 
@@ -49,26 +49,31 @@ struct ContinuousZernikeAnnulusElementMode{T} <: Basis{T}
     L₀₁::AbstractMatrix{T}
     L₁₀::AbstractMatrix{T}
     D::AbstractMatrix{T}
+    via_Jacobi::Bool       # I use this as a flag to quickly switch between _ann2element_via_Jacobi or _ann2element_via_lowering
     b::Int # Should remove once adaptive expansion has been figured out.
 end
 
-function ContinuousZernikeAnnulusElementMode(points::AbstractVector{T}, m::Int, j::Int, L₁₁::AbstractMatrix, L₀₁::AbstractMatrix, L₁₀::AbstractMatrix, D::AbstractMatrix, b::Int) where T
+function ContinuousZernikeAnnulusElementMode(points::AbstractVector{T}, m::Int, j::Int, L₁₁::AbstractMatrix, L₀₁::AbstractMatrix, L₁₀::AbstractMatrix, D::AbstractMatrix, via_Jacobi::Bool, b::Int) where T
     @assert length(points) == 2 && zero(T) < points[1] < points[2] ≤ one(T)
     @assert m ≥ 0
     @assert m == 0 ? j == 1 : 0 ≤ j ≤ 1
     # @assert b ≥ m
-    ContinuousZernikeAnnulusElementMode{T}(points, m, j, L₁₁, L₀₁, L₁₀, D, b)
+    ContinuousZernikeAnnulusElementMode{T}(points, m, j, L₁₁, L₀₁, L₁₀, D, via_Jacobi, b)
 end
 
-function ContinuousZernikeAnnulusElementMode(points::AbstractVector{T}, m::Int, j::Int) where T
+function ContinuousZernikeAnnulusElementMode(points::AbstractVector{T}, m::Int, j::Int, via_Jacobi::Bool=false) where T
     α, β = convert(T, first(points)), convert(T, last(points))
     ρ = α / β
     t = inv(one(T)-ρ^2)
-    (L₁₁, L₀₁, L₁₀) = _ann2element_via_lowering(t, m)
+    if via_Jacobi
+        (L₁₁, L₀₁, L₁₀) = _ann2element_via_Jacobi(t, m)
+    else
+        (L₁₁, L₀₁, L₁₀) = _ann2element_via_lowering(t, m)
+    end
     Z = ZernikeAnnulus{T}(ρ,1,1)
     D = (Z \ (Laplacian(axes(Z,1))*Weighted(Z))).ops[m+1]
 
-    ContinuousZernikeAnnulusElementMode(points, m, j, L₁₁, L₀₁, L₁₀, D, 200) 
+    ContinuousZernikeAnnulusElementMode(points, m, j, L₁₁, L₀₁, L₁₀, D, via_Jacobi, 200)
 end
 
 # ContinuousZernikeAnnulusElementMode(points::AbstractVector, m::Int, j::Int, L₁₁::AbstractMatrix, L₀₁::AbstractMatrix, L₁₀::AbstractMatrix, D::AbstractMatrix, b::Int) = ContinuousZernikeAnnulusElementMode{Float64}(points, m, j, L₁₁, L₀₁, L₁₀, D, b)
@@ -122,8 +127,10 @@ function ldiv(C::ContinuousZernikeAnnulusElementMode{T}, f::AbstractQuasiVector)
     α, β = convert(T, first(C.points)), convert(T, last(C.points))
     ρ = α / β
     m, j = C.m, C.j
-    Z = ZernikeAnnulus{T}(ρ, zero(T), zero(T)) # via lowering
-    # Z = ZernikeAnnulus{T}(ρ, one(T), one(T))   # via Jacobi
+
+    w_a = C.via_Jacobi ? one(T) : zero(T)
+    Z = ZernikeAnnulus{T}(ρ, w_a, w_a)
+
     x = axes(Z,1)
     # # Need to take into account different scalings
     if β ≉  1
@@ -145,7 +152,7 @@ function ldiv(C::ContinuousZernikeAnnulusElementMode{T}, f::AbstractQuasiVector)
     L₁₁, L₀₁, L₁₀ = C.L₁₁, C.L₀₁, C.L₁₀
     R̃ = Hcat(view(L₁₀,1:N,1), view(L₀₁,1:N,1), view(L₁₁,1:N, 1:N-2))
 
-    # convert from ZernikeAnnulus(ρ,0,0) to hats + Bubble
+    # convert from ZernikeAnnulus(ρ,w_a,w_a) to hats + Bubble
     dat = R̃[1:N,1:N] \ c̃
     cfs = T[]
     pad(append!(cfs, dat), axes(C,2))
@@ -161,35 +168,6 @@ function _sum_semiclassicaljacobiweight(t::T, a::Number, b::Number, c::Number) w
     return convert(T, t^c * beta(1+a,1+b) * _₂F₁general2(1+a,-c,2+a+b,1/t))
 end
 
-# @simplify function *(A::QuasiAdjoint{<:Any,<:ContinuousZernikeAnnulusElementMode}, B::ContinuousZernikeAnnulusElementMode)
-#     T = promote_type(eltype(A), eltype(B))
-#     @assert A' == B
-
-#     α, β = convert(T, first(B.points)), convert(T, last(B.points))
-#     ρ = α / β
-#     m = B.m
-
-#     t = inv(one(T)-ρ^2)
-#     L₁₁, L₀₁, L₁₀ = B.L₁₁, B.L₀₁, B.L₁₀
-
-#     # Contribution from the mass matrix of harmonic polynomial
-#     m₀ = convert(T,π) / ( t^(one(T) + m) )
-#     m₀ = m == 0 ? m₀ : m₀ / T(2)
-
-#     jw = _sum_semiclassicaljacobiweight(t,1,1,m) / t^2
-#     M = jw.*L₁₁
-#     a = jw.*L₁₀[1:2,1]
-#     b = jw.*L₀₁[1:2,1]
-
-#     a11 = _sum_semiclassicaljacobiweight(t,2,0,m) / t^2
-#     a12 = jw
-#     a22 = _sum_semiclassicaljacobiweight(t,0,2,m) / t^2
-
-#     C = [a11 a12 a'; a12 a22 b']
-#     M = Hcat(Vcat(C[1:2,3:4]', Zeros{T}(∞,2)), M)
-#     β^2*m₀*Vcat(Hcat(C, Zeros{T}(2,∞)), M)
-# end
-
 @simplify function *(A::QuasiAdjoint{<:Any,<:ContinuousZernikeAnnulusElementMode}, B::ContinuousZernikeAnnulusElementMode)
     T = promote_type(eltype(A), eltype(B))
     @assert A' == B
@@ -201,25 +179,43 @@ end
     t = inv(one(T)-ρ^2)
     L₁₁, L₀₁, L₁₀ = B.L₁₁, B.L₀₁, B.L₁₀
 
-    # Contribution from the mass matrix of harmonic polynomial
-    jw = _sum_semiclassicaljacobiweight(t,0,0,m)
-    m₀ = convert(T,π) / ( t^(one(T) + m) ) * jw
-    m₀ = m == 0 ? m₀ : m₀ / T(2)
-
-
-    M = ApplyArray(*, L₁₁', L₁₁)
+    if B.via_Jacobi
+        # Contribution from the mass matrix of harmonic polynomial
+        m₀ = convert(T,π) / ( t^(one(T) + m) )
+        m₀ = m == 0 ? m₀ : m₀ / T(2)
 
-    a = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₁₀,1:2,1))
-    b = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₀₁,1:2,1))
+        jw = _sum_semiclassicaljacobiweight(t,1,1,m) / t^2
+        M = ApplyArray(*,Diagonal(Fill(jw,∞)),L₁₁)
+        a = jw.*view(L₁₀,1:2,1)
+        b = jw.*view(L₀₁,1:2,1)
 
-    a11 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₁₀,1:2,1))
-    a12 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₀₁,1:2,1))
-    a22 = ApplyArray(*, view(L₀₁',1,1:2)', view(L₀₁,1:2,1))
-
-    # C = Vcat(Hcat(a11[1], a12[1], [a;Zeros{T}(∞)]'), Hcat(a12[1], a22[1], [b;Zeros{T}(∞)]'))
-    C = Vcat(Hcat(a11[1], a12[1], a'), Hcat(a12[1], a22[1], b'))[1:2,1:4]
-    M = Hcat(Vcat(C[1:2,3:4]', Zeros{T}(∞,2)), M)
-    Vcat(Hcat(β^2*m₀*C, Zeros{T}(2,∞)), β^2*m₀*M)
+        a11 = _sum_semiclassicaljacobiweight(t,2,0,m) / t^2
+        a12 = jw
+        a22 = _sum_semiclassicaljacobiweight(t,0,2,m) / t^2
+
+        C = [a11 a12 a'; a12 a22 b']
+        M = Hcat(Vcat(C[1:2,3:4]', Zeros{T}(∞,2)), M)
+        return β^2*m₀*Vcat(Hcat(C, Zeros{T}(2,∞)), M)
+    else
+        # Contribution from the mass matrix of harmonic polynomial
+        jw = _sum_semiclassicaljacobiweight(t,0,0,m)
+        m₀ = convert(T,π) / ( t^(one(T) + m) ) * jw
+        m₀ = m == 0 ? m₀ : m₀ / T(2)
+
+        M = ApplyArray(*, L₁₁', L₁₁)
+
+        a = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₁₀,1:2,1))
+        b = ApplyArray(*, view(L₁₁',1:2,1:2), view(L₀₁,1:2,1))
+
+        a11 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₁₀,1:2,1))
+        a12 = ApplyArray(*, view(L₁₀',1,1:2)', view(L₀₁,1:2,1))
+        a22 = ApplyArray(*, view(L₀₁',1,1:2)', view(L₀₁,1:2,1))
+
+        # C = Vcat(Hcat(a11[1], a12[1], [a;Zeros{T}(∞)]'), Hcat(a12[1], a22[1], [b;Zeros{T}(∞)]'))
+        C = Vcat(Hcat(a11[1], a12[1], a'), Hcat(a12[1], a22[1], b'))[1:2,1:4]
+        M = Hcat(Vcat(C[1:2,3:4]', Zeros{T}(∞,2)), M)
+        return Vcat(Hcat(β^2*m₀*C, Zeros{T}(2,∞)), β^2*m₀*M)
+    end
 end
 
 ###
@@ -245,7 +241,7 @@ end
 # GradientContinuousZernikeAnnulusElementMode(m::Int, j::Int) =  GradientContinuousZernikeAnnulusElementMode([0.0; 1.0], m, j)
 
 axes(Z::GradientContinuousZernikeAnnulusElementMode) = (Inclusion(last(Z.C.points)*UnitDisk{eltype(Z)}()), oneto(∞))
-==(P::GradientContinuousZernikeAnnulusElementMode, Q::GradientContinuousZernikeAnnulusElementMode) = P.C.points == Q.C.points && P.C.m == Q.C.m && P.C.j == Q.C.j
+==(P::GradientContinuousZernikeAnnulusElementMode, Q::GradientContinuousZernikeAnnulusElementMode) = P.C == Q.C
 
 @simplify function *(D::Derivative, C::ContinuousZernikeAnnulusElementMode)
     GradientContinuousZernikeAnnulusElementMode(C)
@@ -323,7 +319,7 @@ function bubble2ann(C::ContinuousZernikeAnnulusElementMode, c::AbstractVector{T}
     N = length(c)
     R̃ = Hcat(view(L₁₀,1:N,1), view(L₀₁,1:N,1), view(L₁₁,1:N, 1:N-2))
 
-    R̃ * c # coefficients for ZernikeAnnulus(ρ,0,0)
+    R̃ * c # coefficients for ZernikeAnnulus(ρ,w_a,w_a)
 end
 
 function plotvalues(u::ApplyQuasiVector{T,typeof(*),<:Tuple{ContinuousZernikeAnnulusElementMode, AbstractVector}}) where T
@@ -349,8 +345,9 @@ function plotvalues(u::ApplyQuasiVector{T,typeof(*),<:Tuple{ContinuousZernikeAnn
     # Scale the grid
     g = scalegrid(grid(C, N), α, β)
 
+    w_a = C.via_Jacobi ? 1 : 0
     # Use fast transforms for synthesis
-    FT = ZernikeAnnulusITransform{T}(N+C.m, 0, 0, 0, ρ)
+    FT = ZernikeAnnulusITransform{T}(N+C.m, w_a, w_a, 0, ρ)
     g, FT * pad(F,blockedrange(oneto(∞)))[Block.(OneTo(N+C.m))]
 end
 

src/diskelement.jl

@@ -69,7 +69,7 @@ function ldiv(C::ContinuousZernikeElementMode{T}, f::AbstractQuasiVector) where
 
     # Seems to cache on its own, no need to memoize unlike annulus
     c = Z0 \ f̃ # Zernike transform
-    c̃ = paddeddata(c)
+    c̃ = ModalTrav(paddeddata(c))
     N = size(c̃.matrix, 1) # degree
 
     # Restrict to relevant mode and add a column corresponding to the hat function.

src/finitecontinuous.jl

@@ -1,15 +1,18 @@
+const VIA_JACOBI = false # toggle this to go via _ann2element_via_Jacobi or _ann2element_via_lowering
+
 struct FiniteContinuousZernike{T, N<:Int, P<:AbstractVector} <: Basis{T}
     N::N
     points::P
     Fs::Tuple{Vararg{FiniteContinuousZernikeMode}}
+    via_Jacobi::Bool
 end
 
-function FiniteContinuousZernike{T}(N::Int, points::AbstractVector, Fs::Tuple{Vararg{FiniteContinuousZernikeMode}}) where {T}
+function FiniteContinuousZernike{T}(N::Int, points::AbstractVector, Fs::Tuple{Vararg{FiniteContinuousZernikeMode}}, via_Jacobi::Bool) where {T}
     @assert length(points) > 1 && points == sort(points)
     @assert length(Fs) == 2N-1
-    FiniteContinuousZernike{T, Int, typeof(points)}(N, points, Fs)
+    FiniteContinuousZernike{T, Int, typeof(points)}(N, points, Fs, via_Jacobi)
 end
-FiniteContinuousZernike(N::Int, points::AbstractVector) = FiniteContinuousZernike{Float64}(N, points, Tuple(_getFs(N, points)))
+FiniteContinuousZernike(N::Int, points::AbstractVector, via_Jacobi::Bool=VIA_JACOBI) = FiniteContinuousZernike{Float64}(N, points, Tuple(_getFs(N, points, via_Jacobi)), via_Jacobi)
 
 function axes(Z::FiniteContinuousZernike{T}) where T
     first(Z.points) ≈ 0 && return (Inclusion(last(Z.points)*UnitDisk{T}()), blockedrange(Vcat(length(Z.points)-1, Fill(length(Z.points) - 1, Z.N-2))))
@@ -56,7 +59,7 @@ function _getMs_ms_js(N::Int)
     (Ms, ms, js)
 end
 
-function _getFs(N::Int, points::AbstractVector{T}) where T
+function _getFs(N::Int, points::AbstractVector{T}, via_Jacobi::Bool) where T
     # Ordered list of Fourier modes (ms, js) and correct length for each Fourier mode Ms.
     Ms, ms, js = _getMs_ms_js(N)
     K = length(points)-1
@@ -74,7 +77,12 @@ function _getFs(N::Int, points::AbstractVector{T}) where T
 
     # Use broadcast notation to compute all the lowering matrices across all
     # intervals and Fourier modes simultaneously.
-    Ls = _ann2element_via_lowering.(ts)
+
+    if via_Jacobi
+        Ls = _ann2element_via_Jacobi.(ts)
+    else
+        Ls = _ann2element_via_lowering.(ts)
+    end
 
     # Use broadcast notation to compute all the derivative matrices across
     # all the intervals and Fourier modes simultaneously
@@ -83,28 +91,31 @@ function _getFs(N::Int, points::AbstractVector{T}) where T
 
     # Loop over the Fourier modes
     Fs = []
+
+    # Pre-calling gives a big speedup.
+    Lss = [[Ls[i][j][1:N+1] for j in 1:3] for i in 1:K+1-κ]
+
     for (M, m, j) in zip(Ms, ms, js)
         # Extract the lowering and differentiation matrices associated
         # with each Fourier mode and store in the Tuples
-        L₁₁ = NTuple{K+1-κ, AbstractMatrix}([Ls[i][1][m+1] for i in 1:K+1-κ])
-        L₀₁ = NTuple{K+1-κ, AbstractMatrix}([Ls[i][2][m+1] for i in 1:K+1-κ])
-        L₁₀ = NTuple{K+1-κ, AbstractMatrix}([Ls[i][3][m+1] for i in 1:K+1-κ])
+        L₁₁ = NTuple{K+1-κ, AbstractMatrix}([Lss[i][1][m+1] for i in 1:K+1-κ])
+        L₀₁ = NTuple{K+1-κ, AbstractMatrix}([Lss[i][2][m+1] for i in 1:K+1-κ])
+        L₁₀ = NTuple{K+1-κ, AbstractMatrix}([Lss[i][3][m+1] for i in 1:K+1-κ])
 
         D = NTuple{K+1-κ, AbstractMatrix}([(Ds[i]).ops[m+1] for i in 1:K+1-κ])
 
         # Construct the structs for each Fourier mode seperately
-        append!(Fs, [FiniteContinuousZernikeMode(M, points, m, j, L₁₁, L₀₁, L₁₀, D, N)])
+        append!(Fs, [FiniteContinuousZernikeMode(M, points, m, j, L₁₁, L₀₁, L₁₀, D, via_Jacobi, N)])
     end
     return Fs
 end
 
-_getFs(F::FiniteContinuousZernike{T}) where T = _getFs(F.N, F.points)
+_getFs(F::FiniteContinuousZernike{T}) where T = _getFs(F.N, F.points, F.via_Jacobi)
 
 function ldiv(F::FiniteContinuousZernike{V}, f::AbstractQuasiVector) where V
     # T = promote_type(V, eltype(f))
     @warn "Expanding via FiniteContinuousZernike is ill-conditioned, please use FiniteZernikeBasis."
     T = V
-    N = F.N; points = T.(F.points)
 
     Fs = F.Fs
     [F \ f.f.(axes(F, 1)) for F in Fs]
@@ -208,7 +219,6 @@ function _bubble2disk_or_ann_all_modes(F::FiniteContinuousZernike, us::AbstractV
                 Ũs[1:Ms[i],i,k] = bubble2disk(m, Us[1:Ms[i],i,k])
             end
         else
-            α, β = points[k], points[k+1]
             for (Fm, i) in zip(Fs, 1:2N-1)
                 C = _getCs(Fm)[k]
                 Ũs[1:Ms[i],i,k] = bubble2ann(C, Us[1:Ms[i],i,k])
@@ -234,11 +244,12 @@ function finite_plotvalues(F::FiniteContinuousZernike, us::AbstractVector)
         else
             α, β = points[k], points[k+1]
             ρ = α / β
-            Z = ZernikeAnnulus{T}(ρ, zero(T), zero(T))
+            w_a = F.via_Jacobi ? one(T) : zero(T)
+            Z = ZernikeAnnulus{T}(ρ, w_a, w_a)
             # Scale the grid
             g = scalegrid(AnnuliOrthogonalPolynomials.grid(Z, Block(N)), α, β)
             # Use fast transforms for synthesis
-            FT = ZernikeAnnulusITransform{T}(N, 0, 0, 0, ρ)
+            FT = ZernikeAnnulusITransform{T}(N, w_a, w_a, 0, ρ)
             val = FT * pad(ModalTrav(Ũs[:,:,k]),axes(Z,2))[Block.(OneTo(N))]
         end
         (θ, r, val) = plot_helper(g, val)