Skip to content

Commit 22d2395

Browse files
dlfivefiftydlfivefifty
committed
Cleanup disk and cone
1 parent f12f4b1 commit 22d2395

File tree

6 files changed

+449
-138
lines changed

6 files changed

+449
-138
lines changed

‎src/Cone/Cone.jl

Lines changed: 84 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,46 @@ spacescompatible(::LegendreConic, ::LegendreConic) = true
1717
domain(::DuffyConic) = Conic()
1818
domain(::LegendreConic) = Conic()
1919

20+
21+
# groups by Fourier
22+
struct ConicTensorizer end
23+
2024
tensorizer(K::DuffyConic) = DiskTensorizer()
21-
tensorizer(K::LegendreConic) = DiskTensorizer()
25+
tensorizer(K::LegendreConic) = ConicTensorizer()
26+
27+
function isqrt_roundup(n)
28+
N = isqrt(n)
29+
N^2 == n ? N : N+1
30+
end
31+
32+
function fromtensor(::ConicTensorizer, A::AbstractMatrix{T}) where T
33+
N = size(A,1)
34+
M = 2N-1
35+
@assert size(A,2) == M
36+
B = PseudoBlockArray(A, Ones{Int}(N), [1; Fill(2,N-1)])
37+
a = Vector{T}()
38+
for N = 1:nblocks(B,2), K=1:N
39+
append!(a, vec(view(B,Block(K,N-K+1))))
40+
end
41+
a
42+
end
43+
44+
function totensor(::ConicTensorizer, a::AbstractVector{T}) where T
45+
N = isqrt_roundup(length(a))
46+
M = 2N-1
47+
à = zeros(eltype(a), N, M)
48+
B = PseudoBlockArray(Ã, Ones{Int}(N), [1; Fill(2,N-1)])
49+
k = 1
50+
for N = 1:nblocks(B,2), K=1:N
51+
V = view(B, Block(K,N-K+1))
52+
for j = 1:length(V)
53+
V[j] = a[k]
54+
k += 1
55+
k > length(a) && return Ã
56+
end
57+
end
58+
Ã
59+
end
2260

2361
rectspace(::DuffyConic) = NormalizedJacobi(0,1,Segment(1,0))*Fourier()
2462

@@ -27,12 +65,32 @@ conemap(txy) = ((t,x,y) = txy; conemap(t,x,y))
2765
conicpolar(t,θ) = conemap(t, cos(θ), sin(θ))
2866
conicpolar(tθ) = ((t,θ) = tθ; conicpolar(t,θ))
2967

30-
points(::Conic, n) = conicpolar.(points(rectspace(DuffyConic()), n))
68+
69+
# M = 2N-1
70+
# N*M == 2N^2 -N == n
71+
# N = (1 + sqrt(1 + 8n)/4)
72+
function points(d::Conic, n)
73+
a,b = rectspace(DuffyConic()).spaces
74+
pts=Array{float(eltype(d))}(undef,0)
75+
N = (1 + isqrt(1+8n)) ÷ 4
76+
M = 2N-1
77+
78+
for y in points(b,M),
79+
x in points(a,N)
80+
push!(pts,conicpolar(Vec(x...,y...)))
81+
end
82+
pts
83+
end
3184
checkpoints(::Conic) = conicpolar.(checkpoints(rectspace(DuffyConic())))
3285

3386

34-
plan_transform(S::DuffyConic, n::AbstractVector) =
35-
TransformPlan(S, plan_transform(rectspace(S),n), Val{false})
87+
function plan_transform(S::DuffyConic, v::AbstractVector{T}) where T
88+
n = length(v)
89+
N = (1 + isqrt(1+8n)) ÷ 4
90+
M = 2N-1
91+
D = plan_transform!(rectspace(S), Array{T}(undef,N,M))
92+
TransformPlan(S, D, Val{false})
93+
end
3694
plan_itransform(S::DuffyConic, n::AbstractVector) =
3795
ITransformPlan(S, plan_itransform(rectspace(S),n), Val{false})
3896

@@ -94,7 +152,7 @@ function _coefficients(triangleplan, v::AbstractVector{T}, ::DuffyConic, ::Legen
94152
fromtensor(LegendreConic(), F)
95153
end
96154

97-
function _coefficients(triangleplan, v::AbstractVector{T}, ::LegendreConic, ::DuffyConic) where T
155+
function _coefficients(triangleplan, v::AbstractVector{T}, ::LegendreConic, ::DuffyConic) where T
98156
F = totensor(LegendreConic(), v)
99157
F = pad(F, :, 2size(F,1)-1)
100158
legendre2duffyconic!(triangleplan, F)
@@ -108,8 +166,7 @@ function coefficients(cfs::AbstractVector{T}, CD::DuffyConic, ZD::LegendreConic)
108166
end
109167

110168
function coefficients(cfs::AbstractVector{T}, ZD::LegendreConic, CD::DuffyConic) where T
111-
c = tridevec(cfs) # TODO: wasteful
112-
n = size(c,1)
169+
n = isqrt_roundup(length(cfs))
113170
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), zero(T)), cfs, ZD, CD)
114171
end
115172

@@ -119,14 +176,21 @@ struct LegendreConicTransformPlan{DUF,CHEB}
119176
end
120177

121178
function LegendreConicTransformPlan(S::LegendreConic, v::AbstractVector{T}) where T
122-
D = plan_transform(DuffyConic(),v)
123-
F = totensor(rectspace(DuffyConic()), D*v) # wasteful, just use to figure out `size(F,1)`
124-
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
179+
n = length(v)
180+
N = (1 + isqrt(1+8n)) ÷ 4
181+
M = 2N-1
182+
D = plan_transform!(rectspace(DuffyConic()), Array{T}(undef,N,M))
183+
P = c_plan_rottriangle(N, zero(T), zero(T), zero(T))
125184
LegendreConicTransformPlan(D,P)
126185
end
127186

128187

129-
*(P::LegendreConicTransformPlan, v::AbstractVector) = _coefficients(P.triangleplan, P.duffyplan*v, DuffyConic(), LegendreConic())
188+
function *(P::LegendreConicTransformPlan, v::AbstractVector)
189+
N,M = P.duffyplan.plan[1][2],P.duffyplan.plan[2][2]
190+
V=reshape(v,N,M)
191+
fromtensor(LegendreConic(),
192+
duffy2legendreconic!(P.triangleplan,P.duffyplan*copy(V)))
193+
end
130194

131195
plan_transform(K::LegendreConic, v::AbstractVector) = LegendreConicTransformPlan(K, v)
132196

@@ -179,14 +243,18 @@ rectspace(::DuffyCone) = NormalizedJacobi(0,1,Segment(1,0))*ZernikeDisk()
179243

180244
blocklengths(d::DuffyCone) = blocklengths(rectspace(d))
181245

246+
# M = N(N+1)/2
247+
# N*M == N^2(N+1)/2 == n
248+
# N = 1/3 (-1 + 1/(-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27 n)))^(1/3) + (-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3))
249+
182250
function points(d::Cone,n)
183-
N = ceil(Int, n^(1/3))
184-
pts=Array{float(eltype(d))}(undef,0)
251+
N = round(Int,1/3*(-1 + 1/(-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3) + (-1 + 27n + 3sqrt(3)sqrt(n*(-2 + 27n)))^(1/3)),RoundUp)
252+
pts = Array{float(eltype(d))}(undef,0)
185253
a,b = rectspace(DuffyCone()).spaces
186-
for y in points(b,N^2), x in points(a,N)
187-
push!(pts,Vec(x...,y...))
254+
for y in points(b,M), x in points(a,N)
255+
push!(pts,conemap(Vec(x...,y...)))
188256
end
189-
conemap.(pts)
257+
pts
190258
end
191259

192260

‎src/Disk/Disk.jl

Lines changed: 90 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -11,11 +11,27 @@ spacescompatible(::ZernikeDisk, ::ZernikeDisk) = true
1111

1212
@containsconstants ZernikeDisk
1313

14-
points(K::ZernikeDisk, n::Integer) =
15-
fromcanonical.(Ref(K), points(ChebyshevDisk(), n))
14+
# M = 2N-1
15+
# N*M == 2N^2 -N == n
16+
# N = (1 + sqrt(1 + 8n)/4)
17+
function points(d::ZernikeDisk, n)
18+
a,b = rectspace(ZernikeDisk()).spaces
19+
pts=Array{float(eltype(domain(d)))}(undef,0)
20+
N = (1 + isqrt(1+8n)) ÷ 4
21+
M = 2N-1
22+
23+
for y in points(b,M),
24+
x in points(a,N)
25+
push!(pts,polar(Vec(x...,y...)))
26+
end
27+
pts
28+
end
1629

17-
evaluate(cfs::AbstractVector, Z::ZernikeDisk, xy) =
18-
Fun(Fun(Z, cfs), ChebyshevDisk())(xy)
30+
function evaluate(cfs::AbstractVector, Z::ZernikeDisk, xy)
31+
C = totensor(ZernikeDisk(), cfs)
32+
F = _coefficients(CDisk2CxfPlan(size(C,1)), C, ZernikeDisk(), ChebyshevDisk())
33+
ProductFun(icheckerboard(F), rectspace(ChebyshevDisk()))(ipolar(xy...)...)
34+
end
1935

2036
Space(d::Disk) = ZernikeDisk(d)
2137

@@ -37,9 +53,42 @@ function points(S::ChebyshevDisk, N)
3753
polar.(pts)
3854
end
3955

40-
DiskTensorizer() = Tensorizer((Ones{Int}(∞),Ones{Int}(∞)))
56+
struct DiskTensorizer end
57+
58+
function fromtensor(::DiskTensorizer, A::AbstractMatrix{T}) where T
59+
N = size(A,1)
60+
M = 4(N-1)+1
61+
@assert size(A,2) == M
62+
B = PseudoBlockArray(A, Ones{Int}(N), [1; Fill(2,2*(N-1))])
63+
a = Vector{T}()
64+
for N = 1:nblocks(B,2), K=1:(N+1)÷2
65+
append!(a, vec(view(B,Block(K,N-2K+2))))
66+
end
67+
a
68+
end
69+
70+
# N + 4 * sum(1:N-1)
71+
# N + 4 * N (N-1)/2 == n
72+
function totensor(::DiskTensorizer, a::AbstractVector{T}) where T
73+
n = length(a)
74+
N = round(Int,1/4*(1 + sqrt(1 + 8n)),RoundUp)
75+
M = 4*(N-1)+1
76+
à = zeros(eltype(a), N, M)
77+
B = PseudoBlockArray(Ã, Ones{Int}(N), [1; Fill(2,2*(N-1))])
78+
k = 1
79+
for N = 1:nblocks(B,2), K=1:(N+1)÷2
80+
V = view(B, Block(K,N-2K+2))
81+
for j = 1:length(V)
82+
V[j] = a[k]
83+
k += 1
84+
k > length(a) && return Ã
85+
end
86+
end
87+
Ã
88+
end
4189

42-
tensorizer(K::ChebyshevDisk) = DiskTensorizer()
90+
tensorizer(K::ZernikeDisk) = DiskTensorizer()
91+
tensorizer(K::ChebyshevDisk) = Tensorizer((Ones{Int}(∞),Ones{Int}(∞)))
4392

4493
# we have each polynomial
4594
blocklengths(K::ChebyshevDisk) = Base.OneTo(∞)
@@ -61,7 +110,11 @@ plan_itransform(S::ChebyshevDisk, n::AbstractVector) =
61110
# drop every other entry
62111
function checkerboard(A::AbstractMatrix{T}) where T
63112
m,n = size(A)
64-
C = zeros(T, (m+1)÷2, n)
113+
#4(N-1) == n-1
114+
N = (n-1)÷4+1
115+
4(N-1)+1 == n || (N += 1) # round up
116+
M = 4(N-1)+1
117+
C = zeros(T, N, M)
65118
z = @view(A[1:2:end,1])
66119
e1,e2 = @view(A[1:2:end,4:4:end]),@view(A[1:2:end,5:4:end])
67120
o1,o2 = @view(A[2:2:end,2:4:end]),@view(A[2:2:end,3:4:end])
@@ -95,64 +148,68 @@ evaluate(cfs::AbstractVector, S::ChebyshevDisk, x) = evaluate(torectcfs(S,cfs),
95148

96149
function _coefficients(disk2cxf, v̂::AbstractVector{T}, ::ChebyshevDisk, ::ZernikeDisk) where T
97150
F = totensor(ChebyshevDisk(), v̂)
98-
F *= sqrt(convert(T,π))
99-
n = size(F,1)
151+
n = disk2cxf.n
100152
F̌ = disk2cxf \ pad(F,n,4n-3)
101-
trivec(F̌)
153+
fromtensor(ZernikeDisk(), F̌)
102154
end
103155

104-
function _coefficients(disk2cxf, v::AbstractVector{T}, ::ZernikeDisk, ::ChebyshevDisk) where T
105-
F̌ = tridevec(v)
106-
n = size(F̌,1)
107-
F = disk2cxf * pad(F̌,n,4n-3)
108-
F /= sqrt(convert(T,π))
109-
fromtensor(ChebyshevDisk(), F)
156+
function _coefficients(disk2cxf, F̌::AbstractMatrix{T}, ::ZernikeDisk, ::ChebyshevDisk) where T
157+
n = disk2cxf.n
158+
disk2cxf * pad(F̌,n,4n-3)
110159
end
111160

161+
_coefficients(disk2cxf, v::AbstractVector{T}, ::ZernikeDisk, ::ChebyshevDisk) where T =
162+
fromtensor(ChebyshevDisk(), _coefficients(disk2cxf, totensor(ZernikeDisk(), v), ZernikeDisk(), ChebyshevDisk()))
163+
112164
function coefficients(cfs::AbstractVector, CD::ChebyshevDisk, ZD::ZernikeDisk)
113165
c = totensor(ChebyshevDisk(), cfs) # TODO: wasteful
114166
n = size(c,1)
115167
_coefficients(CDisk2CxfPlan(n), cfs, CD, ZD)
116168
end
117169

118170
function coefficients(cfs::AbstractVector, ZD::ZernikeDisk, CD::ChebyshevDisk)
119-
c = tridevec(cfs) # TODO: wasteful
171+
c = totensor(ZernikeDisk(), cfs) # TODO: wasteful
120172
n = size(c,1)
121173
_coefficients(CDisk2CxfPlan(n), cfs, ZD, CD)
122174
end
123175

124176

125-
struct FastZernikeDiskTransformPlan{DUF,CHEB}
177+
struct ZernikeDiskTransformPlan{DUF,CHEB}
126178
cxfplan::DUF
127179
disk2cxf::CHEB
128180
end
129181

130-
function FastZernikeDiskTransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
131-
# n = floor(Integer,sqrt(2length(v)) + 1/2)
132-
# v = Array{T}(undef, sum(1:n))
133-
cxfplan = plan_transform(ChebyshevDisk(),v)
134-
c = totensor(ChebyshevDisk(), cxfplan*v) # TODO: wasteful
135-
n = size(c,1)
136-
FastZernikeDiskTransformPlan(cxfplan, CDisk2CxfPlan(n))
182+
function ZernikeDiskTransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
183+
n = length(v)
184+
N = (1 + isqrt(1+8n)) ÷ 4
185+
M = 2N-1
186+
D = plan_transform!(rectspace(ChebyshevDisk()), Array{T}(undef,N,M))
187+
N = (M-1)÷4+1
188+
4(N-1)+1 == M || (N += 1) # round up
189+
ZernikeDiskTransformPlan(D, CDisk2CxfPlan(N))
137190
end
138191

139-
*(P::FastZernikeDiskTransformPlan, v::AbstractVector) =
140-
_coefficients(P.disk2cxf, P.cxfplan*v, ChebyshevDisk(), ZernikeDisk())
192+
function *(P::ZernikeDiskTransformPlan, v::AbstractVector)
193+
N,M = P.cxfplan.plan[1][2],P.cxfplan.plan[2][2]
194+
V = P.cxfplan*reshape(copy(v),N,M)
195+
C = checkerboard(V)
196+
fromtensor(ZernikeDisk(), P.disk2cxf\C)
197+
end
141198

142-
plan_transform(K::ZernikeDisk, v::AbstractVector) = FastZernikeDiskTransformPlan(K, v)
199+
plan_transform(K::ZernikeDisk, v::AbstractVector) = ZernikeDiskTransformPlan(K, v)
143200

144-
struct FastZernikeDiskITransformPlan{DUF,CHEB}
201+
struct ZernikeDiskITransformPlan{DUF,CHEB}
145202
icxfplan::DUF
146203
disk2cxf::CHEB
147204
end
148205

149-
function FastZernikeDiskITransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
206+
function ZernikeDiskITransformPlan(S::ZernikeDisk, v::AbstractVector{T}) where T
150207
# n = floor(Integer,sqrt(2length(v)) + 1/2)
151208
# v = Array{T}(undef, sum(1:n))
152-
FastZernikeDiskITransformPlan(plan_itransform(ChebyshevDisk(), v), CDisk2CxfPlan(n))
209+
ZernikeDiskITransformPlan(plan_itransform(ChebyshevDisk(), v), CDisk2CxfPlan(n))
153210
end
154211

155-
function *(P::FastZernikeDiskITransformPlan, v)
212+
function *(P::ZernikeDiskITransformPlan, v)
156213
# n = floor(Integer,sqrt(2length(v)) + 1/2)
157214
# v = pad(v, sum(1:n))
158215
F̌ = trinormalize!(tridevec(v))
@@ -162,5 +219,5 @@ function *(P::FastZernikeDiskITransformPlan, v)
162219
end
163220

164221

165-
plan_itransform(K::ZernikeDisk, v::AbstractVector) = FastZernikeDiskITransformPlan(K, v)
222+
plan_itransform(K::ZernikeDisk, v::AbstractVector) = ZernikeDiskITransformPlan(K, v)
166223
Base.sum(f::Fun{<:ZernikeDisk}) = Fun(f,ZernikeDisk()).coefficients[1]/π

‎src/MultivariateOrthogonalPolynomials.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ import BlockArrays: blocksizes, BlockSizes, getblock, global2blockindex
1717
import ApproxFunBase: BandedMatrix, blocksize,
1818
linesum,complexlength, BandedBlockBandedMatrix,
1919
real, eps, isapproxinteger, FiniteRange, DFunction,
20-
TransformPlan, ITransformPlan
20+
TransformPlan, ITransformPlan, plan_transform!
2121

2222
# Domains import
2323
import ApproxFunBase: fromcanonical, tocanonical, domainscompatible

‎src/c_transforms.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -66,7 +66,7 @@ end
6666
CDisk2CxfPlan(n::Int) = CDisk2CxfPlan(c_plan_disk2cxf(n), n)
6767

6868
function *(C::CDisk2CxfPlan, A::Matrix{Float64})
69-
(size(A,1) == C.n && size(A,2) == 4C.n-3) || throw(ArgumentError(A))
69+
(size(A,1) == C.n && size(A,2) == 4C.n-3) || throw(DimensionMismatch(""))
7070
B = copy(A)
7171
c_disk2cxf(C.plan, B)
7272
B

0 commit comments

Comments
 (0)