From 7c5c6b3d7a7aa52d6fbd9376bb37915f48128f09 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Fri, 26 May 2023 20:46:31 +0100 Subject: [PATCH 1/7] implementation of partial derivative y from Weighted(Zernike(1)) to Zernike(0) --- src/disk.jl | 30 ++++++++++++++++++++++++++++++ test/test_disk.jl | 18 ++++++++++++++++++ 2 files changed, 48 insertions(+) diff --git a/src/disk.jl b/src/disk.jl index 4595766..97fe9dc 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -342,4 +342,34 @@ function Base._sum(P::Zernike{T}, dims) where T @assert dims == 1 @assert P.a == P.b == 0 Hcat(sqrt(convert(T, π)), Zeros{T}(1,∞)) +end + + +### +# Partial derivatives +### + +normal_jacobi(a::T,b::T,n::Int) where T = sqrt(2^(a+b+1) / (2n + a + b + 1) * gamma(n+a+1)*gamma(n+b+1) / (gamma(n+a+b+1) * factorial(n))) + +@simplify function *(∂ʸ::PartialDerivative{2}, WZ::Weighted{<:Any,<:Zernike}) + @assert WZ.P.a == 0 && WZ.P.b == 1 + T = eltype(eltype(WZ)) + + k = mortar(Base.OneTo.(oneto(∞))) # k counts the the angular mode (+1) + n = mortar(Fill.(oneto(∞),oneto(∞))) .- 1 # n counts the block number which corresponds to the order + m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) + + c=-(normal_jacobi.(T(0), T(1), (n .- m) .÷ 2) ./ normal_jacobi.(T(1), T(0), (n .- m) .÷ 2 .+ m)) + c1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(T(2)) + (1 .- iszero.(m))) + c2 = c1 .* iseven.(n .- k) + c3 = c1 - c2 + + d1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(isodd.(k.-n)) .* (isone.(m) * sqrt(T(2)) + (1 .- isone.(m))) + d2 = d1 .* iseven.(n .- k) + d3 = d1 - d2 + + + A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (2,0)) + B = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d3, Zeros((axes(n,1),)), d2)', axes(n,1), (1,-1), (0,2)) + Zernike{T}(0) * (A + B) end \ No newline at end of file diff --git a/test/test_disk.jl b/test/test_disk.jl index 2cc358f..b45e754 100644 --- a/test/test_disk.jl +++ b/test/test_disk.jl @@ -319,6 +319,24 @@ import ForwardDiff: hessian g = MultivariateOrthogonalPolynomials.plotgrid(W[:,1:3]) @test all(rep[1].args .≈ (first.(g),last.(g),u[g])) end + + @testset "weighted partial derivatives" begin + W = Weighted(Zernike(1)) + Z = Zernike(0) + 𝐱 = axes(W,1) + # ∂ˣ = PartialDerivative{1}(𝐱) + ∂ʸ = PartialDerivative{2}(𝐱) + + ∂Y = Z \ (∂ʸ * W) + B = Block.(1:10); xy = SVector(0.2,0.3) + + for (i,n,m) in zip((1,2,3,4,5,6,14,17), (0,1,1,2,2,2,4,5), (0,-1,1,0,-2,2,-4,1)) + g = 𝐱 -> ForwardDiff.gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n, m, 1, 𝐱), 𝐱)[2] + c = ModalTrav(transform(Z, g)[B])[B] + @test Z[xy,B]'*∂Y[B,i] ≈ Z[xy, B]' * c + end + + end end @testset "Fractional Laplacian on Unit Disk" begin From 2df16c2bd131102006d438eec060091e80194490 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Sat, 27 May 2023 15:31:56 +0100 Subject: [PATCH 2/7] implementation of partial y derivative from Zernike(b) to Zernike(b+1) --- src/disk.jl | 39 ++++++++++++++++++++++++++++++++++++++- test/test_disk.jl | 26 ++++++++++++++++++++------ 2 files changed, 58 insertions(+), 7 deletions(-) diff --git a/src/disk.jl b/src/disk.jl index 97fe9dc..07c807b 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -368,8 +368,45 @@ normal_jacobi(a::T,b::T,n::Int) where T = sqrt(2^(a+b+1) / (2n + a + b + 1) * ga d2 = d1 .* iseven.(n .- k) d3 = d1 - d2 - A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (2,0)) B = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d3, Zeros((axes(n,1),)), d2)', axes(n,1), (1,-1), (0,2)) Zernike{T}(0) * (A + B) +end + +@simplify function *(∂ʸ::PartialDerivative{2}, Z::Zernike) + @assert Z.a == 0 + T = eltype(eltype(Z)) + b = convert(T, Z.b) + + k = mortar(Base.OneTo.(oneto(∞))) # k counts the the angular mode (+1) + n = mortar(Fill.(oneto(∞),oneto(∞))) .- 1 # n counts the block number which corresponds to the order + m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) + + x=axes(Jacobi(0,0),1) + D = BroadcastVector(P->Derivative(x) * P, HalfWeighted{:b}.(Normalized.(Jacobi.(b,1:∞)))) + Ds = BroadcastVector{AbstractMatrix{Float64}}((P,D) -> P \ D, HalfWeighted{:b}.(Normalized.(Jacobi.(b+1,0:∞))) , D) + M = ModalInterlace(Ds, (ℵ₀,ℵ₀), (0,0)) + db = ones(axes(n)) .* (view(view(M, 1:∞, 1:∞),band(0))) + + D = BroadcastVector(P->Derivative(x) * P, Normalized.(Jacobi.(b,-1:∞))) + Ds = BroadcastVector{AbstractMatrix{Float64}}((P,D) -> P \ D, Normalized.(Jacobi.(b+1,0:∞)), D) + Dss = BroadcastVector{AbstractMatrix{Float64}}(P -> Diagonal(view(P, band(1))), Ds) + M = ModalInterlace(Dss, (ℵ₀,ℵ₀), (0,0)) + # d = ones(axes(n)) .* view(Vcat(0,view(view(M, 1:∞, 1:∞),band(0))),1:∞) + d = ones(axes(n)) .* view(view(M, 1:∞, 1:∞),band(0)) + + c1 = d .* (-1).^(isodd.(k.-n)) .* (isone.(m) .* sqrt(T(2)) + (1 .- isone.(m))) + c2 = c1 .* iseven.(n .- k) + c3 = c1 - c2 + + + d1 = db .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(T(2)) + (1 .- iszero.(m))) + d2 = d1 .* iseven.(n .- k) + d3 = d1 - d2 + + + A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (0,2)) + B = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d3, Zeros((axes(n,1),)), d2)', axes(n,1), (1,-1), (2,0)) + + Zernike{T}(Z.b+1) * (A+B)' end \ No newline at end of file diff --git a/test/test_disk.jl b/test/test_disk.jl index b45e754..b499d1e 100644 --- a/test/test_disk.jl +++ b/test/test_disk.jl @@ -320,20 +320,34 @@ import ForwardDiff: hessian @test all(rep[1].args .≈ (first.(g),last.(g),u[g])) end - @testset "weighted partial derivatives" begin + @testset "partial derivatives" begin W = Weighted(Zernike(1)) - Z = Zernike(0) + Z⁰ = Zernike(0) + Z¹ = Zernike(1) + Z² = Zernike(2) + 𝐱 = axes(W,1) # ∂ˣ = PartialDerivative{1}(𝐱) ∂ʸ = PartialDerivative{2}(𝐱) - ∂Y = Z \ (∂ʸ * W) + ∂Y⁰ = Z⁰ \ (∂ʸ * W) + ∂Y¹ = Z¹ \ (∂ʸ * Z⁰) + ∂Y² = Z² \ (∂ʸ * Z¹) + B = Block.(1:10); xy = SVector(0.2,0.3) for (i,n,m) in zip((1,2,3,4,5,6,14,17), (0,1,1,2,2,2,4,5), (0,-1,1,0,-2,2,-4,1)) - g = 𝐱 -> ForwardDiff.gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n, m, 1, 𝐱), 𝐱)[2] - c = ModalTrav(transform(Z, g)[B])[B] - @test Z[xy,B]'*∂Y[B,i] ≈ Z[xy, B]' * c + g⁰ = 𝐱 -> ForwardDiff.gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n, m, 1, 𝐱), 𝐱)[2] + c⁰ = ModalTrav(transform(Z⁰, g⁰)[B])[B] + @test Z⁰[xy,B]'*∂Y⁰[B,i] ≈ Z⁰[xy, B]' * c⁰ + + g¹ = 𝐱 -> ForwardDiff.gradient(𝐱 -> zernikez(n, m, 0, 𝐱), 𝐱)[2] + c¹ = ModalTrav(transform(Z¹, g¹)[B])[B] + @test Z¹[xy,B]'*∂Y¹[B,i] ≈ Z¹[xy, B]' * c¹ + + g² = 𝐱 -> ForwardDiff.gradient(𝐱 -> zernikez(n, m, 1, 𝐱), 𝐱)[2] + c² = ModalTrav(transform(Z², g²)[B])[B] + @test Z²[xy,B]'*∂Y²[B,i] ≈ Z²[xy, B]' * c² end end From f903b14b619126547069741364865771a617224b Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Sat, 27 May 2023 15:49:24 +0100 Subject: [PATCH 3/7] add some comments --- src/disk.jl | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/disk.jl b/src/disk.jl index 07c807b..2dbbda2 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -349,7 +349,8 @@ end # Partial derivatives ### -normal_jacobi(a::T,b::T,n::Int) where T = sqrt(2^(a+b+1) / (2n + a + b + 1) * gamma(n+a+1)*gamma(n+b+1) / (gamma(n+a+b+1) * factorial(n))) +# Helper function for the case weighted(1) to zernike(0) +_normal_jacobi(a::T,b::T,n::Int) where T = sqrt(2^(a+b+1) / (2n + a + b + 1) * gamma(n+a+1)*gamma(n+b+1) / (gamma(n+a+b+1) * gamma(n+1))) @simplify function *(∂ʸ::PartialDerivative{2}, WZ::Weighted{<:Any,<:Zernike}) @assert WZ.P.a == 0 && WZ.P.b == 1 @@ -357,17 +358,20 @@ normal_jacobi(a::T,b::T,n::Int) where T = sqrt(2^(a+b+1) / (2n + a + b + 1) * ga k = mortar(Base.OneTo.(oneto(∞))) # k counts the the angular mode (+1) n = mortar(Fill.(oneto(∞),oneto(∞))) .- 1 # n counts the block number which corresponds to the order - m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) + m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) # Fourier mode number - c=-(normal_jacobi.(T(0), T(1), (n .- m) .÷ 2) ./ normal_jacobi.(T(1), T(0), (n .- m) .÷ 2 .+ m)) + # Raising Fourier mode coefficients, split into sin and cos parts + m=0 change + c=-(_normal_jacobi.(T(0), T(1), (n .- m) .÷ 2) ./ _normal_jacobi.(T(1), T(0), (n .- m) .÷ 2 .+ m)) c1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(T(2)) + (1 .- iszero.(m))) c2 = c1 .* iseven.(n .- k) c3 = c1 - c2 + # Lowering Fourier mode coefficients, split into sin and cos parts + m=1 change d1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(isodd.(k.-n)) .* (isone.(m) * sqrt(T(2)) + (1 .- isone.(m))) d2 = d1 .* iseven.(n .- k) d3 = d1 - d2 + # Put coefficients together A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (2,0)) B = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d3, Zeros((axes(n,1),)), d2)', axes(n,1), (1,-1), (0,2)) Zernike{T}(0) * (A + B) @@ -382,29 +386,33 @@ end n = mortar(Fill.(oneto(∞),oneto(∞))) .- 1 # n counts the block number which corresponds to the order m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) + # Computes the entries for the component that lowers the Fourier mode x=axes(Jacobi(0,0),1) D = BroadcastVector(P->Derivative(x) * P, HalfWeighted{:b}.(Normalized.(Jacobi.(b,1:∞)))) Ds = BroadcastVector{AbstractMatrix{Float64}}((P,D) -> P \ D, HalfWeighted{:b}.(Normalized.(Jacobi.(b+1,0:∞))) , D) M = ModalInterlace(Ds, (ℵ₀,ℵ₀), (0,0)) db = ones(axes(n)) .* (view(view(M, 1:∞, 1:∞),band(0))) + # Computes the entries for the component that lowers the Fourier mode + # the -1 is a hack is the parameters is a trick to make ModalInterlace think that the + # 0-parameter is the second Fourier mode (we use the -1 as a dummy matrix in the ModalInterlace) D = BroadcastVector(P->Derivative(x) * P, Normalized.(Jacobi.(b,-1:∞))) Ds = BroadcastVector{AbstractMatrix{Float64}}((P,D) -> P \ D, Normalized.(Jacobi.(b+1,0:∞)), D) Dss = BroadcastVector{AbstractMatrix{Float64}}(P -> Diagonal(view(P, band(1))), Ds) M = ModalInterlace(Dss, (ℵ₀,ℵ₀), (0,0)) - # d = ones(axes(n)) .* view(Vcat(0,view(view(M, 1:∞, 1:∞),band(0))),1:∞) d = ones(axes(n)) .* view(view(M, 1:∞, 1:∞),band(0)) + # Lowering Fourier mode coefficients, split into sin and cos parts + m=1 change c1 = d .* (-1).^(isodd.(k.-n)) .* (isone.(m) .* sqrt(T(2)) + (1 .- isone.(m))) c2 = c1 .* iseven.(n .- k) c3 = c1 - c2 - + # Raising Fourier mode coefficients, split into sin and cos parts + m=0 change d1 = db .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(T(2)) + (1 .- iszero.(m))) d2 = d1 .* iseven.(n .- k) d3 = d1 - d2 - + # Put coefficients together A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (0,2)) B = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d3, Zeros((axes(n,1),)), d2)', axes(n,1), (1,-1), (2,0)) From da990b6b8790ff02e43f97bb7279654ae6daeb27 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Sat, 27 May 2023 16:08:30 +0100 Subject: [PATCH 4/7] fix tests --- test/test_disk.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_disk.jl b/test/test_disk.jl index b499d1e..f884047 100644 --- a/test/test_disk.jl +++ b/test/test_disk.jl @@ -1,7 +1,7 @@ using MultivariateOrthogonalPolynomials, ClassicalOrthogonalPolynomials, StaticArrays, BlockArrays, BandedMatrices, FastTransforms, LinearAlgebra, RecipesBase, Test, SpecialFunctions, LazyArrays, InfiniteArrays import MultivariateOrthogonalPolynomials: ModalTrav, grid, ZernikeTransform, ZernikeITransform, *, ModalInterlace import ClassicalOrthogonalPolynomials: HalfWeighted, expand -import ForwardDiff: hessian +import ForwardDiff: hessian, gradient @testset "Disk" begin @testset "Transform" begin @@ -337,15 +337,15 @@ import ForwardDiff: hessian B = Block.(1:10); xy = SVector(0.2,0.3) for (i,n,m) in zip((1,2,3,4,5,6,14,17), (0,1,1,2,2,2,4,5), (0,-1,1,0,-2,2,-4,1)) - g⁰ = 𝐱 -> ForwardDiff.gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n, m, 1, 𝐱), 𝐱)[2] + g⁰ = 𝐱 -> gradient(𝐱 -> (1-norm(𝐱)^2)*zernikez(n, m, 1, 𝐱), 𝐱)[2] c⁰ = ModalTrav(transform(Z⁰, g⁰)[B])[B] @test Z⁰[xy,B]'*∂Y⁰[B,i] ≈ Z⁰[xy, B]' * c⁰ - g¹ = 𝐱 -> ForwardDiff.gradient(𝐱 -> zernikez(n, m, 0, 𝐱), 𝐱)[2] + g¹ = 𝐱 -> gradient(𝐱 -> zernikez(n, m, 0, 𝐱), 𝐱)[2] c¹ = ModalTrav(transform(Z¹, g¹)[B])[B] @test Z¹[xy,B]'*∂Y¹[B,i] ≈ Z¹[xy, B]' * c¹ - g² = 𝐱 -> ForwardDiff.gradient(𝐱 -> zernikez(n, m, 1, 𝐱), 𝐱)[2] + g² = 𝐱 -> gradient(𝐱 -> zernikez(n, m, 1, 𝐱), 𝐱)[2] c² = ModalTrav(transform(Z², g²)[B])[B] @test Z²[xy,B]'*∂Y²[B,i] ≈ Z²[xy, B]' * c² end From 11c5f2cf61139a4a3d33a974bec56bf9ed5d5bb0 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Mon, 29 May 2023 11:45:19 +0100 Subject: [PATCH 5/7] performance adjustments --- src/MultivariateOrthogonalPolynomials.jl | 2 +- src/disk.jl | 34 +++++++++++------------- 2 files changed, 17 insertions(+), 19 deletions(-) diff --git a/src/MultivariateOrthogonalPolynomials.jl b/src/MultivariateOrthogonalPolynomials.jl index f3679fe..b012858 100644 --- a/src/MultivariateOrthogonalPolynomials.jl +++ b/src/MultivariateOrthogonalPolynomials.jl @@ -21,7 +21,7 @@ import LazyArrays: arguments, paddeddata, LazyArrayStyle, LazyLayout, PaddedLayo import LazyBandedMatrices: LazyBandedBlockBandedLayout, AbstractBandedBlockBandedLayout, AbstractLazyBandedBlockBandedLayout, _krontrav_axes, DiagTravLayout import InfiniteArrays: InfiniteCardinal, OneToInf -import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw +import ClassicalOrthogonalPolynomials: jacobimatrix, Weighted, orthogonalityweight, HalfWeighted, WeightedBasis, pad, recurrencecoefficients, clenshaw, massmatrix import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan, PartialDerivative, AngularMomentum, BlockOneTo, BlockRange1, interlace, MultivariateOPLayout, MAX_PLOT_BLOCKS diff --git a/src/disk.jl b/src/disk.jl index 2dbbda2..e53e309 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -349,9 +349,6 @@ end # Partial derivatives ### -# Helper function for the case weighted(1) to zernike(0) -_normal_jacobi(a::T,b::T,n::Int) where T = sqrt(2^(a+b+1) / (2n + a + b + 1) * gamma(n+a+1)*gamma(n+b+1) / (gamma(n+a+b+1) * gamma(n+1))) - @simplify function *(∂ʸ::PartialDerivative{2}, WZ::Weighted{<:Any,<:Zernike}) @assert WZ.P.a == 0 && WZ.P.b == 1 T = eltype(eltype(WZ)) @@ -361,15 +358,16 @@ _normal_jacobi(a::T,b::T,n::Int) where T = sqrt(2^(a+b+1) / (2n + a + b + 1) * g m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) # Fourier mode number # Raising Fourier mode coefficients, split into sin and cos parts + m=0 change - c=-(_normal_jacobi.(T(0), T(1), (n .- m) .÷ 2) ./ _normal_jacobi.(T(1), T(0), (n .- m) .÷ 2 .+ m)) - c1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(T(2)) + (1 .- iszero.(m))) + N = sqrt.(massmatrix(Jacobi{T}(0,1)).diag) + c=-BroadcastVector(n->N[n+1], (n .- m) .÷ 2) ./ BroadcastVector(n->N[n+1], (n .- m) .÷ 2 .+ m) + c1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(2*one(T)) + (1 .- iszero.(m))) c2 = c1 .* iseven.(n .- k) - c3 = c1 - c2 + c3 = c1 .* isodd.(n .- k) # Lowering Fourier mode coefficients, split into sin and cos parts + m=1 change - d1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(isodd.(k.-n)) .* (isone.(m) * sqrt(T(2)) + (1 .- isone.(m))) + d1 = c .* ((n .- m) .÷ 2 .+ 1) .* (-1).^(isodd.(k.-n)) .* (isone.(m) * sqrt(2*one(T)) + (1 .- isone.(m))) d2 = d1 .* iseven.(n .- k) - d3 = d1 - d2 + d3 = d1 .* isodd.(n .- k) # Put coefficients together A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (2,0)) @@ -384,33 +382,33 @@ end k = mortar(Base.OneTo.(oneto(∞))) # k counts the the angular mode (+1) n = mortar(Fill.(oneto(∞),oneto(∞))) .- 1 # n counts the block number which corresponds to the order - m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) + m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) # Fourier mode number # Computes the entries for the component that lowers the Fourier mode - x=axes(Jacobi(0,0),1) + x = Inclusion(ChebyshevInterval()) D = BroadcastVector(P->Derivative(x) * P, HalfWeighted{:b}.(Normalized.(Jacobi.(b,1:∞)))) - Ds = BroadcastVector{AbstractMatrix{Float64}}((P,D) -> P \ D, HalfWeighted{:b}.(Normalized.(Jacobi.(b+1,0:∞))) , D) + Ds = BroadcastVector{AbstractMatrix{T}}((P,D) -> P \ D, HalfWeighted{:b}.(Normalized.(Jacobi.(b+1,0:∞))) , D) M = ModalInterlace(Ds, (ℵ₀,ℵ₀), (0,0)) db = ones(axes(n)) .* (view(view(M, 1:∞, 1:∞),band(0))) # Computes the entries for the component that lowers the Fourier mode - # the -1 is a hack is the parameters is a trick to make ModalInterlace think that the + # the -1 in the parameters is a trick to make ModalInterlace think that the # 0-parameter is the second Fourier mode (we use the -1 as a dummy matrix in the ModalInterlace) D = BroadcastVector(P->Derivative(x) * P, Normalized.(Jacobi.(b,-1:∞))) - Ds = BroadcastVector{AbstractMatrix{Float64}}((P,D) -> P \ D, Normalized.(Jacobi.(b+1,0:∞)), D) - Dss = BroadcastVector{AbstractMatrix{Float64}}(P -> Diagonal(view(P, band(1))), Ds) + Ds = BroadcastVector{AbstractMatrix{T}}((P,D) -> P \ D, Normalized.(Jacobi.(b+1,0:∞)), D) + Dss = BroadcastVector{AbstractMatrix{T}}(P -> Diagonal(view(P, band(1))), Ds) M = ModalInterlace(Dss, (ℵ₀,ℵ₀), (0,0)) d = ones(axes(n)) .* view(view(M, 1:∞, 1:∞),band(0)) # Lowering Fourier mode coefficients, split into sin and cos parts + m=1 change - c1 = d .* (-1).^(isodd.(k.-n)) .* (isone.(m) .* sqrt(T(2)) + (1 .- isone.(m))) + c1 = d .* (-1).^(isodd.(k.-n)) .* (isone.(m) .* sqrt(2*one(T)) + (1 .- isone.(m))) c2 = c1 .* iseven.(n .- k) - c3 = c1 - c2 + c3 = c1 .* isodd.(n .- k) # Raising Fourier mode coefficients, split into sin and cos parts + m=0 change - d1 = db .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(T(2)) + (1 .- iszero.(m))) + d1 = db .* (-1).^(iseven.(k.-n)) .* (iszero.(m) * sqrt(2*one(T)) + (1 .- iszero.(m))) d2 = d1 .* iseven.(n .- k) - d3 = d1 - d2 + d3 = d1 .* isodd.(n .- k) # Put coefficients together A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (0,2)) From 794cd83ac1252cef4c5efde5804426056b23d98c Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Wed, 21 Jun 2023 12:59:11 +0100 Subject: [PATCH 6/7] Sheehan's requested changes --- src/disk.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/disk.jl b/src/disk.jl index e53e309..3218405 100644 --- a/src/disk.jl +++ b/src/disk.jl @@ -385,8 +385,7 @@ end m = k .- isodd.(k).*iseven.(n) .- iseven.(k).*isodd.(n) # Fourier mode number # Computes the entries for the component that lowers the Fourier mode - x = Inclusion(ChebyshevInterval()) - D = BroadcastVector(P->Derivative(x) * P, HalfWeighted{:b}.(Normalized.(Jacobi.(b,1:∞)))) + D = BroadcastVector(P->diff(P; dims=1), HalfWeighted{:b}.(Normalized.(Jacobi.(b,1:∞)))) Ds = BroadcastVector{AbstractMatrix{T}}((P,D) -> P \ D, HalfWeighted{:b}.(Normalized.(Jacobi.(b+1,0:∞))) , D) M = ModalInterlace(Ds, (ℵ₀,ℵ₀), (0,0)) db = ones(axes(n)) .* (view(view(M, 1:∞, 1:∞),band(0))) @@ -394,7 +393,7 @@ end # Computes the entries for the component that lowers the Fourier mode # the -1 in the parameters is a trick to make ModalInterlace think that the # 0-parameter is the second Fourier mode (we use the -1 as a dummy matrix in the ModalInterlace) - D = BroadcastVector(P->Derivative(x) * P, Normalized.(Jacobi.(b,-1:∞))) + D = BroadcastVector(P->diff(P; dims=1), Normalized.(Jacobi.(b,-1:∞))) Ds = BroadcastVector{AbstractMatrix{T}}((P,D) -> P \ D, Normalized.(Jacobi.(b+1,0:∞)), D) Dss = BroadcastVector{AbstractMatrix{T}}(P -> Diagonal(view(P, band(1))), Ds) M = ModalInterlace(Dss, (ℵ₀,ℵ₀), (0,0)) @@ -411,8 +410,8 @@ end d3 = d1 .* isodd.(n .- k) # Put coefficients together - A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros((axes(n,1),)), c2)', axes(n,1), (1,-1), (0,2)) - B = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d3, Zeros((axes(n,1),)), d2)', axes(n,1), (1,-1), (2,0)) + A = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, c3, Zeros{T}((axes(n,1),)), c2)', axes(n,1), (1,-1), (0,2)) + B = BlockBandedMatrices._BandedBlockBandedMatrix(BlockBroadcastArray(hcat, d3, Zeros{T}((axes(n,1),)), d2)', axes(n,1), (1,-1), (2,0)) Zernike{T}(Z.b+1) * (A+B)' end \ No newline at end of file From adf6f56c26eb36f3ad0eaa13aded861de8021350 Mon Sep 17 00:00:00 2001 From: ioannisPApapadopoulos Date: Thu, 22 Jun 2023 12:15:08 +0100 Subject: [PATCH 7/7] restart tests