Skip to content

Store SemiclassicalJacobiFamily #14

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dlfivefifty wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Store SemiclassicalJacobiFamily #14

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
b4b3eff
Store SemiclassicalJacobiFamily
dlfivefifty Mar 21, 2025
3e0d148
Use SemiclassicalJacobiFamily
dlfivefifty Mar 21, 2025
65e0046
Work on faster conversion
dlfivefifty Mar 24, 2025
a465362
work on tests
dlfivefifty Mar 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/AnnuliOrthogonalPolynomials.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,10 @@ using BlockArrays, BlockBandedMatrices, ClassicalOrthogonalPolynomials, Continuu
import Base: in, axes, getindex, broadcasted, tail, +, -, *, /, \, convert, OneTo, show, summary, ==, oneto, diff
import BlockArrays: block, blockindex, _BlockedUnitRange, blockcolsupport
import ClassicalOrthogonalPolynomials: checkpoints, ShuffledR2HC, TransformFactorization, ldiv, paddeddata, jacobimatrix, orthogonalityweight, SetindexInterlace
import ContinuumArrays: Weight, weight, grid, ℵ₁, ℵ₀, unweighted, plan_transform
import ContinuumArrays: Weight, weight, grid, ℵ₁, ℵ₀, unweighted, plan_transform, plotgrid
import HarmonicOrthogonalPolynomials: BivariateOrthogonalPolynomial, MultivariateOrthogonalPolynomial, Plan
import MultivariateOrthogonalPolynomials: BlockOneTo, ModalInterlace, laplacian, ModalTrav
import SemiclassicalOrthogonalPolynomials: divdiff, HalfWeighted
import SemiclassicalOrthogonalPolynomials: divdiff, HalfWeighted, SemiclassicalJacobiFamily

export Block, SVector, CircleCoordinate, ZernikeAnnulus, ComplexZernikeAnnulus, Laplacian

Expand Down
62 changes: 37 additions & 25 deletions src/annulus.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,8 @@ struct ZernikeAnnulus{T} <: AbstractZernikeAnnulus{T}
ρ::T
a::T
b::T
ZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b)
P::SemiclassicalJacobiFamily{T}
ZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b, SemiclassicalJacobiFamily(inv(1-ρ^2),b,a,0:∞))
end

"""
Expand All @@ -49,7 +50,8 @@ struct ComplexZernikeAnnulus{T} <: AbstractZernikeAnnulus{Complex{T}}
ρ::T
a::T
b::T
ComplexZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b)
P::SemiclassicalJacobiFamily{T}
ComplexZernikeAnnulus{T}(ρ::T, a::T, b::T) where T = new{T}(ρ, a, b, SemiclassicalJacobiFamily(inv(1-ρ^2),b,a,0:∞))
end


Expand All @@ -71,7 +73,9 @@ copy(A::AbstractZernikeAnnulus) = A

orthogonalityweight(Z::AbstractZernikeAnnulus) = AnnulusWeight(Z.ρ, Z.a, Z.b)

zernikeannulusr(ρ, ℓ, m, a, b, r::T) where T = r^m * SemiclassicalJacobi{T}(inv(1-ρ^2),b,a,m)[(r^2 - 1)/(ρ^2 - 1), (ℓ-m) ÷ 2 + 1]
zernikeannulusr(ρ, ℓ, m, a, b, r, P) = r^m * P[(r^2 - 1)/(ρ^2 - 1), (ℓ-m) ÷ 2 + 1]

zernikeannulusr(ρ, ℓ, m, a, b, r::T) where T = zernikeannulusr(ρ, ℓ, m, a, b, r, SemiclassicalJacobi{T}(inv(1-ρ^2),b,a,m))
function zernikeannulusz(ρ, ℓ, ms, a, b, rθ::RadialCoordinate{T}) where T
r,θ = rθ.r,rθ.θ
m = abs(ms)
Expand All @@ -97,15 +101,19 @@ function getindex(Z::ZernikeAnnulus{T}, rθ::RadialCoordinate, B::BlockIndex{1})
ℓ = Int(block(B))-1
k = blockindex(B)
m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
zernikeannulusz(Z.ρ, ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
(; r, θ) = rθ
(; ρ, a, b, P) = Z
(isodd(k+ℓ) ? cos(m*θ) : sin(m*θ)) * zernikeannulusr(ρ, ℓ, m, a, b, r, P[m+1])
end


function getindex(Z::ComplexZernikeAnnulus{T}, rθ::RadialCoordinate, B::BlockIndex{1}) where T
ℓ = Int(block(B))-1
k = blockindex(B)
m = iseven(ℓ) ? k-isodd(k) : k-iseven(k)
complexzernikeannulusz(Z.ρ, ℓ, (isodd(k+ℓ) ? 1 : -1) * m, Z.a, Z.b, rθ)
(; r, θ) = rθ
(; ρ, a, b, P) = Z
exp(im*(isodd(k+ℓ) ? m : -m)*θ) * zernikeannulusr(ρ, ℓ, m, a, b, r, P[m+1])
end


Expand All @@ -125,25 +133,29 @@ end
function \(A::ZernikeAnnulus{T}, B::Weighted{V,ZernikeAnnulus{V}}) where {T,V}
TV = promote_type(T,V)
(A.a == B.P.a == A.b == B.P.b == 0 && A.ρ == B.P.ρ) && return Eye{TV}(∞)
@assert A.a == A.b == 1
@assert B.P.a == B.P.b == 1
@assert A.ρ == B.P.ρ

ρ = convert(TV, A.ρ); t=inv(one(TV)-ρ^2)

# L₁ = Weighted.(SemiclassicalJacobi{real(TV)}.(t,zero(TV),zero(TV),zero(TV):∞)) .\ Weighted.(SemiclassicalJacobi{real(TV)}.(t,one(TV),one(TV),zero(TV):∞))
# L₂ = SemiclassicalJacobi{real(TV)}.(t,one(TV),one(TV),zero(TV):∞) .\ SemiclassicalJacobi{real(TV)}.(t,zero(TV),zero(TV),zero(TV):∞)

Q₀₀ = SemiclassicalJacobi{real(TV)}.(t,0,0,0:∞)
Q₁₁ = SemiclassicalJacobi{real(TV)}.(t,1,1,0:∞)

L₁ = Weighted.(Q₀₀) .\ Weighted.(Q₁₁)
L₂ = Q₁₁ .\ Q₀₀

# L = (one(TV)-ρ^2)^2 .* (L₂ .* L₁)
# Workaround for broken lazy multiplication
L = BroadcastVector{AbstractMatrix{TV}}((L2, L1) -> (one(TV)-ρ^2)^2 .* ApplyArray(*,L2,L1), L₂, L₁)
ModalInterlace{TV}(L, (ℵ₀,ℵ₀), (4, 4))
Q₁₁ = B.P.P # SemiclassicalJacobi{real(TV)}.(t,1,1,0:∞)
if A.a == A.b == 0
Q₀₀ = A.P # SemiclassicalJacobi{real(TV)}.(t,0,0,0:∞)
L = BroadcastVector{AbstractMatrix{TV}}(L1 -> (one(TV)-ρ^2)^2 * L1, Weighted.(Q₀₀) .\ Weighted.(Q₁₁))
ModalInterlace{TV}(L, (ℵ₀,ℵ₀), (4, 0))
else
@assert A.a == A.b == 1
Q₀₀ = SemiclassicalJacobi{real(TV)}.(t,0,0,0:∞)

# L₁ = Weighted.(SemiclassicalJacobi{real(TV)}.(t,zero(TV),zero(TV),zero(TV):∞)) .\ Weighted.(SemiclassicalJacobi{real(TV)}.(t,one(TV),one(TV),zero(TV):∞))
# L₂ = SemiclassicalJacobi{real(TV)}.(t,one(TV),one(TV),zero(TV):∞) .\ SemiclassicalJacobi{real(TV)}.(t,zero(TV),zero(TV),zero(TV):∞)

L₁ = Weighted.(Q₀₀) .\ Weighted.(Q₁₁)
L₂ = Q₁₁ .\ Q₀₀

# L = (one(TV)-ρ^2)^2 .* (L₂ .* L₁)
# Workaround for broken lazy multiplication
L = BroadcastVector{AbstractMatrix{TV}}((L2, L1) -> (one(TV)-ρ^2)^2 .* ApplyArray(*,L2,L1), L₂, L₁)
ModalInterlace{TV}(L, (ℵ₀,ℵ₀), (4, 4))
end
end


Expand All @@ -157,7 +169,7 @@ function laplacian(W::Weighted{<:Any,<:ZernikeAnnulus})
@assert P.a == P.b == 1
ρ = P.ρ; t = inv(1-ρ^2)
T = eltype(P)
Ps = SemiclassicalJacobi{T}.(t,1,1,0:∞)
Ps = P.P
Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 4t*(1-ρ^2)^2*divdiff(HalfWeighted{:c}(C), HalfWeighted{:c}(B))*divdiff(HalfWeighted{:ab}(B), HalfWeighted{:ab}(A)), Ps, SemiclassicalJacobi.(t,0,0,1:∞), Ps)
P * ModalInterlace(Δs, (ℵ₀,ℵ₀), (2,2))
end
Expand All @@ -166,7 +178,7 @@ function laplacian(P::ZernikeAnnulus)
ρ,a,b = P.ρ,P.a,P.b
t = inv(1-ρ^2)
T = eltype(P)
Ps = SemiclassicalJacobi.(t,b,a,0:∞)
Ps = P.P
Δs = BroadcastVector{AbstractMatrix{T}}((C,B,A) -> 4t*divdiff(HalfWeighted{:c}(C), HalfWeighted{:c}(B))*divdiff(B, A), SemiclassicalJacobi.(t,b+2,a+2,0:∞), SemiclassicalJacobi.(t,b+1,a+1,1:∞), Ps)
ZernikeAnnulus(ρ,a+2,b+2) * ModalInterlace(Δs, (ℵ₀,ℵ₀), (-2,6))
end
Expand Down Expand Up @@ -241,8 +253,8 @@ function denormalize_annulus(A::AbstractVector, a, b, c, ρ, analysis=true)
w = AnnulusWeight(ρ, a, b)
constants = normalize_mmodes(w)[1:l] # m-mode constants
d = [inv(constants[mm+1]*ss) for (mm, ss) in zip(m, s)] # multiply by relevant (-1)
analysis && return d.*A # multiply vector by denormalization if analysis
A ./ d # divide vector by denormalization if synthesis
analysis && return broadcast!(*, A, d, A) # multiply vector by denormalization if analysis
broadcast!(/, A, A, d) # divide vector by denormalization if synthesis
end

# # FastTransforms uses orthonormalized annulus OPs so we need to correct the normalization
Expand Down
26 changes: 17 additions & 9 deletions test/test_annulus.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -155,24 +155,32 @@ import LazyArrays: Ones
end

@testset "Weighted" begin
P = ZernikeAnnulus(ρ,1,1)
W = Weighted(P)
Δ = P \ (Laplacian(axes(P,1)) * W)
Q = ZernikeAnnulus(ρ,1,1)
W = Weighted(Q)
P = ZernikeAnnulus(ρ)

c = transform(Q, splat((x,y) -> exp(x*cos(y))))
x,y = 𝐱 = SVector(0.7,0.2)
@test (W * c)[𝐱] ≈ (1 - norm(𝐱)^2)*(norm(𝐱)^2-ρ^2) * exp(x*cos(y))
KR = Block.(Base.OneTo(100))
@test (P[:,KR] * ((P\W)[KR,KR] * c[KR]))[𝐱] ≈ (1 - norm(𝐱)^2)*(norm(𝐱)^2-ρ^2) * exp(x*cos(y))

Δ = Q \ (Laplacian(axes(Q,1)) * W)

xy = SVector(0.5,0.1); r = norm(xy); τ = (1-r^2)/(1-ρ^2)
@test W[xy, 1:3] ≈ [0.007400000000000006;0.0007400000000000007;0.003700000000000003]
@test Weighted(ZernikeAnnulus{eltype(xy)}(ρ,1,1))[xy,1] ≈ (1-r^2) * (r^2-ρ^2) * ZernikeAnnulus{eltype(xy)}(ρ,1,1)[xy,1] ≈ (1-ρ^2)^2 * HalfWeighted{:ab}(SemiclassicalJacobi.(t,1,1,0:∞)[1])[τ,1]
@test tr(hessian(xy -> Weighted(ZernikeAnnulus{eltype(xy)}(ρ,1,1))[xy,1], xy)) ≈ P[xy,1:4]'* Δ[1:4,1]
@test tr(hessian(xy -> Weighted(ZernikeAnnulus{eltype(xy)}(ρ,1,1))[xy,1], xy)) ≈ Q[xy,1:4]'* Δ[1:4,1]

# TODO get hessian working again
# c = [randn(100); zeros(∞)]
# @test tr(hessian(xy -> (Weighted(ZernikeAnnulus{eltype(xy)}(ρ,1,1))*c)[xy], xy)) ≈ (P*(Δ*c))[xy]
# @test tr(hessian(xy -> (Weighted(ZernikeAnnulus{eltype(xy)}(ρ,1,1))*c)[xy], xy)) ≈ (Q*(Δ*c))[xy]

c = [Ones(10); zeros(∞)]
@test (P*(Δ*c))[xy] ≈ -0.4023800762027685
@test (Q*(Δ*c))[xy] ≈ -0.4023800762027685

L = P \ W
@test W[xy, 1:10] ≈ (P[xy, 1:50]'*L[1:50,1:50])[1:10]
L = Q \ W
@test W[xy, 1:10] ≈ (Q[xy, 1:50]'*L[1:50,1:50])[1:10]
end
end

Expand Down Expand Up @@ -223,7 +231,7 @@ import LazyArrays: Ones
B = ComplexZernikeAnnulus(ρ,0,1)

xy = SVector(0.5,0.1)
@test A[xy, 1:3] ≈ [1.0 + 0.0im;0.5 + 0.10000000000000002im;0.5 + 0.10000000000000002im]
@test A[xy, 1:3] ≈ [1.0 + 0.0im;0.5 - 0.10000000000000002im;0.5 + 0.10000000000000002im]
R = B \ A
@test A[SVector(0.5,0.1), Block.(1:3)]' ≈ B[SVector(0.5,0.1), Block.(1:3)]' * R[Block.(1:3),Block.(1:3)]

Expand Down
Loading