Skip to content

Commit 0884eb0

Browse files
dlfivefiftydlfivefifty
committed
Cone transform works
1 parent ee3b13c commit 0884eb0

File tree

2 files changed

+115
-39
lines changed

2 files changed

+115
-39
lines changed

‎src/Cone/Cone.jl

Lines changed: 63 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -179,10 +179,10 @@ rectspace(::DuffyCone) = NormalizedJacobi(0,1,Segment(1,0))*ZernikeDisk()
179179

180180
blocklengths(d::DuffyCone) = blocklengths(rectspace(d))
181181

182-
function points(sp::DuffyCone,n)
182+
function points(d::Cone,n)
183183
N = ceil(Int, n^(1/3))
184-
pts=Array{float(eltype(domain(sp)))}(undef,0)
185-
a,b = rectspace(sp).spaces
184+
pts=Array{float(eltype(d))}(undef,0)
185+
a,b = rectspace(DuffyCone()).spaces
186186
for y in points(b,N^2), x in points(a,N)
187187
push!(pts,Vec(x...,y...))
188188
end
@@ -258,52 +258,90 @@ function evaluate(cfs::AbstractVector, S::DuffyCone, txy::Vec{3})
258258
Fun(b, view(mat,1,:))(x/t,y/t)
259259
end
260260

261+
function duffy2legendrecone_column_J!(P, Fc, F, Jin)
262+
J = sum(1:Jin)
263+
for j = Jin:size(Fc,2)
264+
if 1 ≤ J ≤ size(F,2)
265+
Fc[:,j] .= view(F,:,J)
266+
else
267+
Fc[:,j] .= 0
268+
end
269+
J += j
270+
end
271+
c_execute_tri_lo2hi(P, Fc)
272+
J = sum(1:Jin)
273+
for j = Jin:size(Fc,2)
274+
J > size(F,2) && break
275+
F[:,J] .= view(Fc,:,j)
276+
J += j
277+
end
278+
F
279+
end
280+
281+
function legendre2duffycone_column_J!(P, Fc, F, Jin)
282+
J = sum(1:Jin)
283+
for j = Jin:size(Fc,2)
284+
if 1 ≤ J ≤ size(F,2)
285+
Fc[:,j] .= view(F,:,J)
286+
else
287+
Fc[:,j] .= 0
288+
end
289+
J += j
290+
end
291+
c_execute_tri_hi2lo(P, Fc)
292+
J = sum(1:Jin)
293+
for j = Jin:size(Fc,2)
294+
J > size(F,2) && break
295+
F[:,J] .= view(Fc,:,j)
296+
J += j
297+
end
298+
F
299+
end
300+
261301
function duffy2legendrecone!(triangleplan, F::AbstractMatrix)
262-
Fc = F[:,1:2:end]
263-
c_execute_tri_lo2hi(triangleplan, Fc)
264-
F[:,1:2:end] .= Fc
265-
# ignore first column
266-
Fc[:,2:end] .= F[:,2:2:end]
267-
c_execute_tri_lo2hi(triangleplan, Fc)
268-
F[:,2:2:end] .= Fc[:,2:end]
302+
Fc = Matrix{Float64}(undef,size(F,1),size(F,1))
303+
for J = 1:(isqrt(1 + 8size(F,2))-1)÷ 2 # actually div
304+
duffy2legendrecone_column_J!(triangleplan, Fc, F, J)
305+
end
269306
F
270307
end
271308

272309
function legendre2duffycone!(triangleplan, F::AbstractMatrix)
273-
Fc = F[:,1:2:end]
274-
c_execute_tri_hi2lo(triangleplan, Fc)
275-
F[:,1:2:end] .= Fc
276-
# ignore first column
277-
Fc[:,2:end] .= F[:,2:2:end]
278-
c_execute_tri_hi2lo(triangleplan, Fc)
279-
F[:,2:2:end] .= Fc[:,2:end]
310+
Fc = Matrix{Float64}(undef,size(F,1),size(F,1))
311+
for J = 1:(isqrt(1 + 8size(F,2))-1)÷ 2 # actually div
312+
legendre2duffycone_column_J!(triangleplan, Fc, F, J)
313+
end
280314
F
281315
end
282316

283317
function _coefficients(triangleplan, v::AbstractVector{T}, ::DuffyCone, ::LegendreCone) where T
284318
F = totensor(rectspace(DuffyCone()), v)
285-
F = pad(F, :, 2size(F,1)-1)
319+
for j = 1:size(F,2)
320+
F[:,j] = coefficients(F[:,j], NormalizedJacobi(0,1,Segment(1,0)), NormalizedJacobi(0,2,Segment(1,0)))
321+
end
286322
duffy2legendrecone!(triangleplan, F)
287323
fromtensor(LegendreCone(), F)
288324
end
289325

290-
function _coefficients(triangleplan, v::AbstractVector{T}, ::LegendreCone, ::DuffyCone) where T
326+
function _coefficients(triangleplan, v::AbstractVector{T}, ::LegendreCone, ::DuffyCone) where T
291327
F = totensor(LegendreCone(), v)
292-
F = pad(F, :, 2size(F,1)-1)
293328
legendre2duffycone!(triangleplan, F)
329+
for j = 1:size(F,2)
330+
F[:,j] = coefficients(F[:,j], NormalizedJacobi(0,2,Segment(1,0)), NormalizedJacobi(0,1,Segment(1,0)))
331+
end
294332
fromtensor(rectspace(DuffyCone()), F)
295333
end
296334

297335
function coefficients(cfs::AbstractVector{T}, CD::DuffyCone, ZD::LegendreCone) where T
298336
c = totensor(rectspace(DuffyCone()), cfs) # TODO: wasteful
299337
n = size(c,1)
300-
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), zero(T)), cfs, CD, ZD)
338+
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), one(T)), cfs, CD, ZD)
301339
end
302340

303341
function coefficients(cfs::AbstractVector{T}, ZD::LegendreCone, CD::DuffyCone) where T
304-
c = tridevec(cfs) # TODO: wasteful
342+
c = totensor(LegendreCone(), cfs) # TODO: wasteful
305343
n = size(c,1)
306-
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), zero(T)), cfs, ZD, CD)
344+
_coefficients(c_plan_rottriangle(n, zero(T), zero(T), one(T)), cfs, ZD, CD)
307345
end
308346

309347
struct LegendreConeTransformPlan{DUF,CHEB}
@@ -314,7 +352,7 @@ end
314352
function LegendreConeTransformPlan(S::LegendreCone, v::AbstractVector{T}) where T
315353
D = plan_transform(DuffyCone(),v)
316354
F = totensor(rectspace(DuffyCone()), D*v) # wasteful, just use to figure out `size(F,1)`
317-
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
355+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))
318356
LegendreConeTransformPlan(D,P)
319357
end
320358

@@ -330,7 +368,7 @@ end
330368

331369
function LegendreConeITransformPlan(S::LegendreCone, v::AbstractVector{T}) where T
332370
F = fromtensor(LegendreCone(), v) # wasteful, just use to figure out `size(F,1)`
333-
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
371+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))
334372
D = plan_itransform(DuffyCone(), _coefficients(P, v, S, DuffyCone()))
335373
LegendreConeITransformPlan(D,P)
336374
end

‎test/test_cone.jl

Lines changed: 52 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
using ApproxFun, MultivariateOrthogonalPolynomials, StaticArrays, FillArrays, Test
22
import ApproxFunBase: checkpoints
3-
import MultivariateOrthogonalPolynomials: rectspace, totensor, duffy2legendreconic!, legendre2duffyconic!, c_plan_rottriangle, plan_transform
3+
import MultivariateOrthogonalPolynomials: rectspace, totensor, duffy2legendreconic!, legendre2duffyconic!, c_plan_rottriangle, plan_transform,
4+
c_execute_tri_hi2lo, c_execute_tri_lo2hi, duffy2legendrecone_column_J!, duffy2legendrecone!, legendre2duffycone!
45

56
@testset "Conic" begin
67
@testset "DuffyConic" begin
@@ -67,6 +68,7 @@ import MultivariateOrthogonalPolynomials: rectspace, totensor, duffy2legendrecon
6768
end
6869

6970

71+
7072
@testset "Cone" begin
7173
@testset "rectspace" begin
7274
rs = rectspace(DuffyCone())
@@ -102,30 +104,66 @@ end
102104
@test f(0.3,0.1,0.2) ≈ exp(cos(0.3*0.1)*0.2)
103105
end
104106

105-
@testset "Legendre<>DuffyConic" begin
107+
@testset "TriTransform" begin
108+
m = 1
109+
p = Fun(NormalizedJacobi(0,2m+2,0..1), [1])
110+
q = Fun(x -> p(x)*(1-x)^m * 2^m, NormalizedJacobi(0,2,0..1))
111+
F = [zeros(2) q.coefficients]
112+
F = pad(F, :, 2size(F,1)-1)
113+
T = Float64
114+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))
115+
c_execute_tri_lo2hi(P, F)
116+
@test F ≈ (B = zeros(size(F)); B[1,m+1] = 1; B)
117+
end
118+
119+
@testset "Legendre<>DuffyCone" begin
120+
for (m,k,ℓ) in ((0,0,0), (0,0,1), (0,0,2), (1,0,0), (1,0,1), (1,1,0), (1,1,1), (1,1,2))
121+
Y = Fun(ZernikeDisk(), [Zeros(sum(1:m)+k); 1])
122+
f = (txy) -> ((t,x,y) = txy; θ = atan(y,x); Fun(NormalizedJacobi(0,2m+2,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * Y(x/t,y/t))
123+
g = Fun(f, DuffyCone())
124+
a = g.coefficients
125+
F = totensor(rectspace(DuffyCone()), a)
126+
for j = 1:size(F,2)
127+
F[:,j] = coefficients(F[:,j], NormalizedJacobi(0,1,Segment(1,0)), NormalizedJacobi(0,2,Segment(1,0)))
128+
end
129+
130+
T = eltype(a)
131+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))
132+
133+
Fc = Matrix{Float64}(undef,size(F,1),size(F,1))
134+
for J = 1:size(F,2)
135+
duffy2legendrecone_column_J!(P, Fc, F, J)
136+
end
137+
138+
@test F ≈ (B = zeros(size(F)); B[ℓ+1,sum(1:m)+k+1] = 1; B)
139+
end
106140
for k = 0:10
107141
a = [zeros(k); 1.0; zeros(5)]
108-
F = totensor(rectspace(DuffyConic()), a)
142+
F = totensor(rectspace(DuffyCone()), a)
109143
F = pad(F, :, 2size(F,1)-1)
144+
110145
T = eltype(a)
111-
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), zero(T))
112-
@test legendre2duffyconic!(P, duffy2legendreconic!(P, copy(F))) ≈ F
146+
P = c_plan_rottriangle(size(F,1), zero(T), zero(T), one(T))
147+
148+
@test legendre2duffycone!(P, duffy2legendrecone!(P, copy(F))) ≈ F
113149

114-
b = coefficients(a, LegendreConic(), DuffyConic())
115-
@test a ≈ coefficients(b, DuffyConic(), LegendreConic())[1:length(a)]
150+
b = coefficients(a, LegendreCone(), DuffyCone())
151+
@test a ≈ coefficients(b, DuffyCone(), LegendreCone())[1:length(a)]
116152
end
117153
end
118-
119-
@testset "LegendreConicPlan" begin
120-
for (m,ℓ) in ((0,0), (0,1), (0,2), (1,1), (1,2), (2,2))
121-
f = (txy) -> ((t,x,y) = txy; θ = atan(y,x); Fun(NormalizedJacobi(0,2m+1,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * cos(m*θ))
122-
g = Fun(f, LegendreConic())
123-
t,x,y = sqrt(0.1^2+0.2^2),0.1,0.2
154+
155+
156+
@testset "LegendreConePlan" begin
157+
for (m,k,ℓ) in ((0,0,0), )
158+
Y = Fun(ZernikeDisk(), [Zeros(sum(1:m)+k); 1])
159+
f = (txy) -> ((t,x,y) = txy; θ = atan(y,x); Fun(NormalizedJacobi(0,2m+2,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * Y(x,y))
160+
g = Fun(f, LegendreCone(),20)
161+
t,x,y = 0.3,0.1,0.2
124162
@test g(t,x,y) ≈ f((t,x,y))
125163
end
126164
end
127165

128-
@testset "LegendreConic" begin
166+
@testset "LegendreCone" begin
129167
f = Fun((t,x,y) -> 1, LegendreConic(), 10)
130168
@test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
131169
@test f(sqrt(0.1^2+0.2^2),0.1,0.2) ≈ 1

0 commit comments

Comments
 (0)