Skip to content

Commit e82cb05

Browse files
authored
Merge pull request #5 from JuliaApproximation/dl/approxfunbase
WIP Disk Transform / ApproxFunBase support
2 parents ff182bd + 35b773b commit e82cb05

22 files changed

+605
-83
lines changed

.gitignore

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
2+
Manifest.toml

Project.toml

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
name = "MultivariateOrthogonalPolynomials"
2+
uuid = "4f6956fd-4f93-5457-9149-7bfc4b2ce06d"
3+
4+
[deps]
5+
ApproxFun = "28f2ccd6-bb30-5033-b560-165f7b14dc2f"
6+
ApproxFunBase = "fbd15aa5-315a-5a7d-a8a4-24992e37be05"
7+
ApproxFunFourier = "59844689-9c9d-51bf-9583-5b794ec66d30"
8+
ApproxFunOrthogonalPolynomials = "b70543e2-c0d9-56b8-a290-0d4d6d4de211"
9+
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
10+
BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
11+
DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
12+
FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
13+
FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c"
14+
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
15+
InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
16+
InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
17+
LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
18+
Libdl = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
19+
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
20+
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
21+
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
22+
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
23+
SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
24+
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
25+
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
26+
27+
[compat]
28+
ApproxFun = "0.11"
29+
ApproxFunBase = "≥ 0.0.3"
30+
BandedMatrices = "0.9"
31+
BlockArrays = "0.8"
32+
DomainSets = "≥ 0.0.1"
33+
FastGaussQuadrature = "≥ 0.3.0"
34+
FastTransforms = "≥ 0.4.1"
35+
FillArrays = "≥ 0.5"
36+
InfiniteArrays = "0.0.3"
37+
LazyArrays = "0.8"
38+
RecipesBase = "≥ 0.5.0"
39+
SpecialFunctions = "≥ 0.6.0"
40+
StaticArrays = "≥ 0.3.0"
41+
julia = "≥ 0.7.0"

REQUIRE

Lines changed: 0 additions & 11 deletions
This file was deleted.

src/Cone/Cone.jl

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
export Cone, DuffyCone, LegendreCone
2+
3+
struct Cone <: Domain{Vec{3,Float64}} end
4+
struct DuffyCone <: Space{Cone,Float64} end
5+
struct LegendreCone <: Space{Cone,Float64} end
6+
7+
@containsconstants DuffyCone
8+
@containsconstants LegendreCone
9+
10+
spacescompatible(::DuffyCone, ::DuffyCone) = true
11+
spacescompatible(::LegendreCone, ::LegendreCone) = true
12+
13+
domain(::DuffyCone) = Cone()
14+
domain(::LegendreCone) = Cone()
15+
16+
tensorizer(K::DuffyCone) = DiskTensorizer()
17+
tensorizer(K::LegendreCone) = DiskTensorizer()
18+
19+
rectspace(::DuffyCone) = NormalizedJacobi(0,1,Segment(1,0))*Fourier()
20+
21+
conepolar(t,θ) = Vec(t, t*cos(θ), t*sin(θ))
22+
conepolar(tθ) = ((t,θ) = tθ; conepolar(t,θ))
23+
24+
points(::Cone, n) = conepolar.(points(rectspace(DuffyCone()), n))
25+
checkpoints(::Cone) = conepolar.(checkpoints(rectspace(DuffyCone())))
26+
27+
28+
plan_transform(S::DuffyCone, n::AbstractVector) =
29+
TransformPlan(S, plan_transform(rectspace(S),n), Val{false})
30+
plan_itransform(S::DuffyCone, n::AbstractVector) =
31+
ITransformPlan(S, plan_itransform(rectspace(S),n), Val{false})
32+
33+
*(P::TransformPlan{<:Any,<:DuffyCone}, v::AbstractArray) = P.plan*v
34+
*(P::ITransformPlan{<:Any,<:DuffyCone}, v::AbstractArray) = P.plan*v
35+
36+
evaluate(cfs::AbstractVector, S::DuffyCone, txy) = evaluate(cfs, rectspace(S), ipolar(txy[Vec(2,3)]))
37+
evaluate(cfs::AbstractVector, S::LegendreCone, txy) = evaluate(coefficients(cfs, S, DuffyCone()), DuffyCone(), txy)
38+
39+
40+
# function transform(::DuffyCone, v)
41+
# n = (1 + isqrt(1+8*length(v))) ÷ 4
42+
# F = copy(reshape(v,n,2n-1))
43+
# Pt = plan_transform(NormalizedJacobi(0,1,Segment(1,0)), n)
44+
# Pθ = plan_transform(Fourier(), 2n-1)
45+
# for j = 1:size(F,2)
46+
# F[:,j] = Pt*F[:,j]
47+
# end
48+
# for k = 1:size(F,1)
49+
# F[k,:] = Pθ * F[k,:]
50+
# end
51+
# fromtensor(rectspace(DuffyCone()), F)
52+
# end
53+
54+
function evaluate(cfs::AbstractVector, S::DuffyCone, txy::Vec{3})
55+
t,x,y = txy
56+
@assert t sqrt(x^2+y^2)
57+
θ = atan(y,x)
58+
Fun(rectspace(S), cfs)(t,θ)
59+
end
60+
61+
62+
function _coefficients(triangleplan, v::AbstractVector{T}, ::DuffyCone, ::LegendreCone) where T
63+
F = totensor(rectspace(DuffyCone()), v)
64+
F = pad(F, :, 2size(F,1)-1)
65+
Fc = F[:,1:2:end]
66+
c_execute_tri_lo2hi(triangleplan, Fc)
67+
F[:,1:2:end] .= Fc
68+
# ignore first column
69+
Fc[:,2:end] .= F[:,2:2:end]
70+
c_execute_tri_lo2hi(triangleplan, Fc)
71+
F[:,2:2:end] .= Fc[:,2:end]
72+
fromtensor(LegendreCone(), F)
73+
end
74+
75+
function _coefficients(triangleplan, v::AbstractVector{T}, ::LegendreCone, ::DuffyCone) where T
76+
F = totensor(LegendreCone(), v)
77+
F = pad(F, :, 2size(F,1)-1)
78+
Fc = F[:,1:2:end]
79+
c_execute_tri_hi2lo(triangleplan, Fc)
80+
F[:,1:2:end] .= Fc
81+
# ignore first column
82+
Fc[:,2:end] .= F[:,2:2:end]
83+
c_execute_tri_hi2lo(triangleplan, Fc)
84+
F[:,2:2:end] .= Fc[:,2:end]
85+
fromtensor(DuffyCone(), F)
86+
end
87+
88+
function coefficients(cfs::AbstractVector{T}, CD::DuffyCone, ZD::LegendreCone) where T
89+
c = totensor(rectspace(DuffyCone()), cfs) # TODO: wasteful
90+
n = size(c,1)
91+
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), zero(T)), cfs, CD, ZD)
92+
end
93+
94+
function coefficients(cfs::AbstractVector{T}, ZD::LegendreCone, CD::DuffyCone) where T
95+
c = tridevec(cfs) # TODO: wasteful
96+
n = size(c,1)
97+
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), zero(T)), cfs, ZD, CD)
98+
end
99+
100+
struct LegendreConeTransformPlan{DUF,CHEB}
101+
duffyplan::DUF
102+
triangleplan::CHEB
103+
end
104+
105+
function LegendreConeTransformPlan(S::LegendreCone, v::AbstractVector{T}) where T
106+
D = plan_transform(DuffyCone(),v)
107+
F = totensor(rectspace(DuffyCone()), D*v) # wasteful, just use to figure out `size(F,1)`
108+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
109+
LegendreConeTransformPlan(D,P)
110+
end
111+
112+
113+
*(P::LegendreConeTransformPlan, v::AbstractVector) = _coefficients(P.triangleplan, P.duffyplan*v, DuffyCone(), LegendreCone())
114+
115+
plan_transform(K::LegendreCone, v::AbstractVector) = LegendreConeTransformPlan(K, v)
116+
117+
struct LegendreConeITransformPlan{DUF,CHEB}
118+
iduffyplan::DUF
119+
triangleplan::CHEB
120+
end
121+
122+
function LegendreConeITransformPlan(S::LegendreCone, v::AbstractVector{T}) where T
123+
F = fromtensor(LegendreCone(), v) # wasteful, just use to figure out `size(F,1)`
124+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
125+
D = plan_itransform(DuffyCone(), _coefficients(P, v, S, DuffyCone()))
126+
LegendreConeITransformPlan(D,P)
127+
end
128+
129+
130+
*(P::LegendreConeITransformPlan, v::AbstractVector) = P.iduffyplan*_coefficients(P.triangleplan, v, LegendreCone(), DuffyCone())
131+
132+
133+
plan_itransform(K::LegendreCone, v::AbstractVector) = LegendreConeITransformPlan(K, v)
134+

src/Disk/.DS_Store

6 KB
Binary file not shown.

src/Disk/Disk.jl

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
export ChebyshevDisk, ZernikeDisk
2+
3+
struct ZernikeDisk{V,T} <: Space{Disk{V},T}
4+
domain::Disk{V}
5+
end
6+
7+
ZernikeDisk(d::Disk{V}) where V = ZernikeDisk{V,complex(eltype(V))}(d)
8+
ZernikeDisk() = ZernikeDisk(Disk())
9+
10+
spacescompatible(::ZernikeDisk, ::ZernikeDisk) = true
11+
12+
@containsconstants ZernikeDisk
13+
14+
points(K::ZernikeDisk, n::Integer) =
15+
fromcanonical.(Ref(K), points(ChebyshevDisk(), n))
16+
17+
evaluate(cfs::AbstractVector, Z::ZernikeDisk, xy) =
18+
Fun(Fun(Z, cfs), ChebyshevDisk())(xy)
19+
20+
Space(d::Disk) = ZernikeDisk(d)
21+
22+
struct ChebyshevDisk{V,T} <: Space{Disk{V},T}
23+
domain::Disk{V}
24+
end
25+
26+
@containsconstants ChebyshevDisk
27+
28+
ChebyshevDisk(d::Disk{V}) where V = ChebyshevDisk{V,eltype(V)}(d)
29+
ChebyshevDisk() = ChebyshevDisk(Disk())
30+
31+
rectspace(_) = Chebyshev() * Fourier()
32+
33+
spacescompatible(K1::ChebyshevDisk, K2::ChebyshevDisk) = domainscompatible(K1, K2)
34+
35+
function points(S::ChebyshevDisk, N)
36+
pts = points(rectspace(S), N)
37+
polar.(pts)
38+
end
39+
40+
DiskTensorizer() = Tensorizer((Ones{Int}(∞),Ones{Int}(∞)))
41+
42+
tensorizer(K::ChebyshevDisk) = DiskTensorizer()
43+
44+
# we have each polynomial
45+
blocklengths(K::ChebyshevDisk) = 1:
46+
47+
for OP in (:block,:blockstart,:blockstop)
48+
@eval begin
49+
$OP(s::ChebyshevDisk,M::Block) = $OP(tensorizer(s),M)
50+
$OP(s::ChebyshevDisk,M) = $OP(tensorizer(s),M)
51+
end
52+
end
53+
54+
55+
plan_transform(S::ChebyshevDisk, n::AbstractVector) =
56+
TransformPlan(S, plan_transform(rectspace(S),n), Val{false})
57+
plan_itransform(S::ChebyshevDisk, n::AbstractVector) =
58+
ITransformPlan(S, plan_itransform(rectspace(S),n), Val{false})
59+
60+
# drop every other entry
61+
function checkerboard(A::AbstractMatrix{T}) where T
62+
m,n = size(A)
63+
C = zeros(T, (m+1)÷2, n)
64+
z = @view(A[1:2:end,1])
65+
e1,e2 = @view(A[1:2:end,4:4:end]),@view(A[1:2:end,5:4:end])
66+
o1,o2 = @view(A[2:2:end,2:4:end]),@view(A[2:2:end,3:4:end])
67+
C[1:size(z,1),1] .= z
68+
C[1:size(e1,1),4:4:n] .= e1
69+
C[1:size(e2,1),5:4:n] .= e2
70+
C[1:size(o1,1),2:4:n] .= o1
71+
C[1:size(o2,1),3:4:n] .= o2
72+
C
73+
end
74+
75+
function icheckerboard(C::AbstractMatrix{T}) where T
76+
m,n = size(C)
77+
A = zeros(T, 2*m, n)
78+
A[1:2:end,1] .= C[:,1]
79+
A[1:2:end,4:4:end] .= C[:,4:4:end]
80+
A[1:2:end,5:4:end] .= C[:,5:4:end]
81+
A[2:2:end,2:4:end] .= C[:,2:4:end]
82+
A[2:2:end,3:4:end] .= C[:,3:4:end]
83+
A
84+
end
85+
86+
torectcfs(S, v) = fromtensor(rectspace(S), icheckerboard(totensor(S, v)))
87+
fromrectcfs(S, v) = fromtensor(S, checkerboard(totensor(rectspace(S),v)))
88+
89+
*(P::TransformPlan{<:Any,<:ChebyshevDisk}, v::AbstractArray) = fromrectcfs(P.space, P.plan*v)
90+
*(P::ITransformPlan{<:Any,<:ChebyshevDisk}, v::AbstractArray) = P.plan*torectcfs(P.space, v)
91+
92+
evaluate(cfs::AbstractVector, S::ChebyshevDisk, x) = evaluate(torectcfs(S,cfs), rectspace(S), ipolar(x))
93+
94+
95+
function _coefficients(disk2cxf, v̂::AbstractVector{T}, ::ChebyshevDisk, ::ZernikeDisk) where T
96+
F = totensor(ChebyshevDisk(), v̂)
97+
F *= sqrt(convert(T,π))
98+
n = size(F,1)
99+
= disk2cxf \ pad(F,n,4n-3)
100+
trivec(F̌)
101+
end
102+
103+
function _coefficients(disk2cxf, v::AbstractVector{T}, ::ZernikeDisk, ::ChebyshevDisk) where T
104+
= tridevec(v)
105+
n = size(F̌,1)
106+
F = disk2cxf * pad(F̌,n,4n-3)
107+
F /= sqrt(convert(T,π))
108+
fromtensor(ChebyshevDisk(), F)
109+
end
110+
111+
function coefficients(cfs::AbstractVector, CD::ChebyshevDisk, ZD::ZernikeDisk)
112+
c = totensor(ChebyshevDisk(), cfs) # TODO: wasteful
113+
n = size(c,1)
114+
_coefficients(CDisk2CxfPlan(n), cfs, CD, ZD)
115+
end
116+
117+
function coefficients(cfs::AbstractVector, ZD::ZernikeDisk, CD::ChebyshevDisk)
118+
c = tridevec(cfs) # TODO: wasteful
119+
n = size(c,1)
120+
_coefficients(CDisk2CxfPlan(n), cfs, ZD, CD)
121+
end
122+
123+
124+
struct FastZernikeDiskTransformPlan{DUF,CHEB}
125+
cxfplan::DUF
126+
disk2cxf::CHEB
127+
end
128+
129+
function FastZernikeDiskTransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
130+
# n = floor(Integer,sqrt(2length(v)) + 1/2)
131+
# v = Array{T}(undef, sum(1:n))
132+
cxfplan = plan_transform(ChebyshevDisk(),v)
133+
c = totensor(ChebyshevDisk(), cxfplan*v) # TODO: wasteful
134+
n = size(c,1)
135+
FastZernikeDiskTransformPlan(cxfplan, CDisk2CxfPlan(n))
136+
end
137+
138+
*(P::FastZernikeDiskTransformPlan, v::AbstractVector) =
139+
_coefficients(P.disk2cxf, P.cxfplan*v, ChebyshevDisk(), ZernikeDisk())
140+
141+
plan_transform(K::ZernikeDisk, v::AbstractVector) = FastZernikeDiskTransformPlan(K, v)
142+
143+
struct FastZernikeDiskITransformPlan{DUF,CHEB}
144+
icxfplan::DUF
145+
disk2cxf::CHEB
146+
end
147+
148+
function FastZernikeDiskITransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
149+
# n = floor(Integer,sqrt(2length(v)) + 1/2)
150+
# v = Array{T}(undef, sum(1:n))
151+
FastZernikeDiskITransformPlan(plan_itransform(ChebyshevDisk(), v), CDisk2CxfPlan(n))
152+
end
153+
154+
function *(P::FastZernikeDiskITransformPlan, v)
155+
# n = floor(Integer,sqrt(2length(v)) + 1/2)
156+
# v = pad(v, sum(1:n))
157+
= trinormalize!(tridevec(v))
158+
F = P.disk2cxf *
159+
= trivec(transpose(F))
160+
P.icxfplan*
161+
end
162+
163+
164+
plan_itransform(K::ZernikeDisk, v::AbstractVector) = FastZernikeDiskITransformPlan(K, v)
165+
Base.sum(f::Fun{<:ZernikeDisk}) = Fun(f,ZernikeDisk()).coefficients[1]/π
File renamed without changes.
File renamed without changes.

0 commit comments

Comments
 (0)