Skip to content

Commit 8c98abd

Browse files
dlfivefiftydlfivefifty
committed
fix cone transform
1 parent 57624fc commit 8c98abd

File tree

3 files changed

+126
-98
lines changed

3 files changed

+126
-98
lines changed

‎src/Cone/Cone.jl

Lines changed: 67 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,8 @@ domain(::LegendreConic) = Conic()
2121
# groups by Fourier
2222
struct ConicTensorizer end
2323

24+
25+
2426
tensorizer(K::DuffyConic) = DiskTensorizer()
2527
tensorizer(K::LegendreConic) = ConicTensorizer()
2628

@@ -64,6 +66,8 @@ conemap(t,x,y) = Vec(t, t*x, t*y)
6466
conemap(txy) = ((t,x,y) = txy; conemap(t,x,y))
6567
conicpolar(t,θ) = conemap(t, cos(θ), sin(θ))
6668
conicpolar(tθ) = ((t,θ) = tθ; conicpolar(t,θ))
69+
conepolar(t,r,θ) = conemap(t, r*cos(θ), r*sin(θ))
70+
conepolar(v::Vec) = conepolar(v...)
6771

6872

6973
# M = 2N-1
@@ -74,11 +78,12 @@ function pointsize(::Conic, n)
7478
N,2N-1
7579
end
7680

81+
pointsize(s::Space, n) = pointsize(domain(s), n)
7782

7883
function points(d::Conic, n)
7984
a,b = rectspace(DuffyConic()).spaces
8085
pts=Array{float(eltype(d))}(undef,0)
81-
N,M = pointsize(d)
86+
N,M = pointsize(d,n)
8287

8388
for y in points(b,M),
8489
x in points(a,N)
@@ -254,57 +259,65 @@ function pointsize(::Cone, n)
254259
N, N, 2N-1
255260
end
256261

262+
function points(d::Cone, n)
263+
a,_ = rectspace(DuffyCone()).spaces
264+
b,c = rectspace(ZernikeDisk()).spaces
257265

258-
function points(d::Cone,n)
259-
N,M = pointsize(d,n)
260-
pts = Array{ApproxFunBase.float(eltype(d))}(undef,0)
261-
a,b = rectspace(DuffyCone()).spaces
262-
for y in points(b,M), x in points(a,N)
263-
push!(pts,conemap(Vec(x...,y...)))
264-
end
265-
pts
266+
M,_,N = pointsize(d, n)
267+
p_a = points(a,M)
268+
p_b = points(b,M)
269+
p_c = points(c,N)
270+
271+
conepolar.(Vec.(p_a,reshape(p_b,1,M),reshape(p_c,1,1,N)))
266272
end
267273

268274

275+
276+
# function points(d::Cone,n)
277+
# N,M = pointsize(d,n)
278+
# pts = Array{ApproxFunBase.float(eltype(d))}(undef,0)
279+
# a,b = rectspace(DuffyCone()).spaces
280+
# for y in points(b,M), x in points(a,N)
281+
# push!(pts,conemap(Vec(x...,y...)))
282+
# end
283+
# pts
284+
# end
285+
286+
269287
checkpoints(::Cone) = conemap.(checkpoints(rectspace(DuffyCone())))
270288

271289

272-
function plan_transform(sp::DuffyCone, n::AbstractVector{T}) where T
273-
rs = rectspace(sp)
274-
N = ceil(Int, length(n)^(1/3))
275-
M = length(points(rs.spaces[2], N^2))
276-
@assert N*M == length(n)
277-
TransformPlan(sp,((plan_transform(rs.spaces[1],T,N),N), (plan_transform(rs.spaces[2],T,M),M)), Val{false})
290+
function plan_transform(sp::DuffyCone, V::AbstractArray{T,3}) where T
291+
M,M̃,N = size(V)
292+
@assert M == M̃
293+
a,b = rectspace(sp).spaces
294+
TransformPlan(sp, (plan_transform(a, M), plan_transform(b, Array{T}(undef,M,N))), Val{false})
278295
end
296+
279297
plan_itransform(S::DuffyCone, n::AbstractVector) =
280298
ITransformPlan(S, plan_itransform(rectspace(S),n), Val{false})
281299

282-
function *(T::TransformPlan{<:Any,<:DuffyCone,true}, M::AbstractMatrix)
283-
n=size(M,1)
284-
285-
for k=1:size(M,2)
286-
M[:,k]=T.plan[1][1]*M[:,k]
300+
function tensortransform(P::TransformPlan{<:Any,<:DuffyCone,false}, V::AbstractArray{<:Any,3})
301+
M,_,N = size(V)
302+
R = Array{Float64}(undef,M,sum(1:M))
303+
for k = 1:M
304+
R[k,:] = P.plan[2]*V[k,:,:]
287305
end
288-
for k=1:n
289-
M[k,:] = (T.plan[2][1]*M[k,:])[1:size(M,2)]
306+
for j = 1:size(R,2)
307+
R[:,j] = P.plan[1]*R[:,j]
290308
end
291-
M
309+
R
292310
end
293311

294-
function *(T::TransformPlan{<:Any,<:DuffyCone,true},v::AbstractVector)
295-
N,M = T.plan[1][2],T.plan[2][2]
296-
V=reshape(v,N,M)
297-
fromtensor(T.space,T*V)
298-
end
299-
300-
function *(T::TransformPlan{<:Any,SS,false},v::AbstractVector) where {SS<:DuffyCone}
301-
P = TransformPlan(T.space,T.plan,Val{true})
302-
P*AbstractVector{rangetype(SS)}(v)
303-
end
312+
*(P::TransformPlan{<:Any,<:DuffyCone,false}, V::AbstractArray{<:Any,3}) =
313+
fromtensor(P.space,tensortransform(P, V))
304314

305315
*(P::ITransformPlan{<:Any,<:DuffyCone}, v::AbstractArray) = P.plan*v
306316

307-
evaluate(cfs::AbstractVector, S::DuffyCone, txy) = evaluate(cfs, rectspace(S), ipolar(txy[Vec(2,3)]))
317+
# function evaluate(cfs::AbstractVector, S::DuffyCone, txy)
318+
# t,x,y = txy
319+
# evaluate(cfs, rectspace(S), Vec(t,
320+
# end
308321
evaluate(cfs::AbstractVector, S::LegendreCone, txy) = evaluate(coefficients(cfs, S, DuffyCone()), DuffyCone(), txy)
309322

310323

@@ -375,28 +388,31 @@ function legendre2duffycone_column_J!(P, Fc, F, Jin)
375388
end
376389

377390
function duffy2legendrecone!(triangleplan, F::AbstractMatrix)
378-
Fc = Matrix{Float64}(undef,size(F,1),size(F,1))
391+
Fc = Matrix{eltype(F)}(undef,size(F,1),size(F,1))
379392
for J = 1:(isqrt(1 + 8size(F,2))-1)÷ 2 # actually div
380393
duffy2legendrecone_column_J!(triangleplan, Fc, F, J)
381394
end
382395
F
383396
end
384397

385398
function legendre2duffycone!(triangleplan, F::AbstractMatrix)
386-
Fc = Matrix{Float64}(undef,size(F,1),size(F,1))
399+
Fc = Matrix{eltype(F)}(undef,size(F,1),size(F,1))
387400
for J = 1:(isqrt(1 + 8size(F,2))-1)÷ 2 # actually div
388401
legendre2duffycone_column_J!(triangleplan, Fc, F, J)
389402
end
390403
F
391404
end
392405

393-
function _coefficients(triangleplan, v::AbstractVector{T}, ::DuffyCone, ::LegendreCone) where T
394-
F = totensor(rectspace(DuffyCone()), v)
406+
function _coefficients!(triangleplan, F::AbstractMatrix{T}, ::DuffyCone, ::LegendreCone) where T
395407
for j = 1:size(F,2)
396408
F[:,j] = coefficients(F[:,j], NormalizedJacobi(0,1,Segment(1,0)), NormalizedJacobi(0,2,Segment(1,0)))
397409
end
398410
duffy2legendrecone!(triangleplan, F)
399-
fromtensor(LegendreCone(), F)
411+
end
412+
413+
function _coefficients(triangleplan, v::AbstractVector{T}, D::DuffyCone, L::LegendreCone) where T
414+
F = totensor(rectspace(DuffyCone()), v)
415+
fromtensor(LegendreCone(), _coefficients!(triangleplan, F, D, L))
400416
end
401417

402418
function _coefficients(triangleplan, v::AbstractVector{T}, ::LegendreCone, ::DuffyCone) where T
@@ -425,24 +441,24 @@ struct LegendreConeTransformPlan{DUF,CHEB}
425441
triangleplan::CHEB
426442
end
427443

428-
function LegendreConeTransformPlan(S::LegendreCone, v::AbstractVector{T}) where T
429-
n = length(v)
430-
N,M = pointsize(Cone(), n)
431-
D = plan_transform!(rectspace(DuffyCone()), Array{T}(undef,N,M))
432-
P = c_plan_rottriangle(N, zero(T), zero(T), zero(T))
444+
function LegendreConeTransformPlan(S::LegendreCone, V::AbstractArray{T,3}) where T
445+
M,M̃,N = size(V)
446+
@assert M̃ == M
447+
D = plan_transform(DuffyCone(), V)
448+
P = c_plan_rottriangle(M, zero(T), zero(T), one(T))
433449
LegendreConeTransformPlan(D,P)
434450
end
435451

436452

437-
function *(P::LegendreConeTransformPlan, v::AbstractVector)
438-
N,M = P.duffyplan.plan[1][2],P.duffyplan.plan[2][2]
439-
V=reshape(v,N,M)
440-
cfs = P.duffyplan*copy(V)
441-
duffy2legendrecone!(P.triangleplan,cfs)
442-
fromtensor(LegendreCone(), cfs)
453+
function tensortransform(P::LegendreConeTransformPlan, V::AbstractArray{<:Any,3})
454+
C = tensortransform(P.duffyplan,V)
455+
_coefficients!(P.triangleplan,C, DuffyCone(), LegendreCone())
443456
end
444457

445-
plan_transform(K::LegendreCone, v::AbstractVector) = LegendreConeTransformPlan(K, v)
458+
*(P::LegendreConeTransformPlan, V::AbstractArray{<:Any,3}) =
459+
fromtensor(LegendreCone(), tensortransform(P,V))
460+
461+
plan_transform(K::LegendreCone, V::AbstractArray{<:Any,3}) = LegendreConeTransformPlan(K, V)
446462

447463
struct LegendreConeITransformPlan{DUF,CHEB}
448464
iduffyplan::DUF

‎src/Disk/Disk.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -191,7 +191,7 @@ end
191191
function *(P::ZernikeDiskTransformPlan, V::AbstractMatrix)
192192
V = P.cxfplan*copy(V)
193193
C = checkerboard(V)
194-
fromtensor(ZernikeDisk(), P.disk2cxf\C)
194+
fromtensor(ZernikeDisk(), P.disk2cxf\C)[1:sum(1:size(V,1))]
195195
end
196196

197197
plan_transform(K::ZernikeDisk, v::AbstractMatrix) = ZernikeDiskTransformPlan(K, v)

‎test/test_cone.jl

Lines changed: 58 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@ using ApproxFun, MultivariateOrthogonalPolynomials, StaticArrays, FillArrays, Bl
22
import ApproxFunBase: checkpoints, plan_transform!, fromtensor
33
import MultivariateOrthogonalPolynomials: rectspace, totensor, duffy2legendreconic!, legendre2duffyconic!, c_plan_rottriangle, plan_transform,
44
c_execute_tri_hi2lo, c_execute_tri_lo2hi, duffy2legendrecone_column_J!, duffy2legendrecone!, legendre2duffycone!,
5-
duffy2legendreconic!, ConicTensorizer, pointsize
5+
duffy2legendreconic!, ConicTensorizer, pointsize, tensortransform, _coefficients!
66

77

88
@testset "Conic" begin
@@ -119,32 +119,43 @@ import MultivariateOrthogonalPolynomials: rectspace, totensor, duffy2legendrecon
119119
end
120120

121121
@testset "Cone" begin
122-
@testset "rectspace" begin
123-
rs = rectspace(DuffyCone())
124-
@test points(rs,10) isa Vector{SVector{3,Float64}}
125-
@test_broken @inferred(checkpoints(rs))
126-
@test checkpoints(rs) isa Vector{SVector{3,Float64}}
127-
end
128-
129122
@testset "ConeTensorizer" begin
130123
ts = MultivariateOrthogonalPolynomials.ConeTensorizer()
131124
@test totensor(ts,1:10) == [1 2 3 5 6 7; 4 8 9 0 0 0; 10 0 0 0 0 0]
132125
end
133126

134127
@testset "DuffyCone" begin
135128
p = points(DuffyCone(), 10)
136-
@test p isa Vector{SVector{3,Float64}}
137-
P = plan_transform(DuffyCone(), Vector{Float64}(undef, length(p)))
138-
139-
@test P * fill(1.0, length(p)) ≈ [1/sqrt(2); Zeros(55)] ≈ [Fun((x,y) -> 1, ZernikeDisk()).coefficients; Zeros(55)]
129+
@test p isa Array{SVector{3,Float64},3}
130+
V = fill(1.0, size(p))
131+
M,_,N = size(V)
132+
R = Array{Float64}(undef,M,sum(1:M))
133+
P = plan_transform(DuffyCone(), Array{Float64}(undef, size(p)))
134+
for k = 1:M
135+
R[k,:] = P.plan[2]*V[k,:,:]
136+
end
137+
for j = 1:size(M,2)
138+
R[:,j] = P.plan[1]*R[:,j]
139+
end
140+
@test R ≈ [1/sqrt(2) 0 0; 0 0 0]
141+
142+
@test P * V ≈ [1/sqrt(2); Zeros(9)] ≈ [Fun((x,y) -> 1, ZernikeDisk()).coefficients; Zeros(9)]
140143

141144
f = Fun((t,x,y) -> 1, DuffyCone(), 10)
142-
@test f.coefficients ≈ [1/sqrt(2); Zeros(55)]
145+
@test f.coefficients ≈ [1/sqrt(2); Zeros(9)]
143146
@test f(0.3,0.1,0.2) ≈ 1
144147

145148
f = Fun((t,x,y) -> t, DuffyCone(), 10)
146149
@test f(0.3,0.1,0.2) ≈ 0.3
147150

151+
g = Fun(ZernikeDisk(),[0,1.0])
152+
p = points(DuffyCone(), 10)
153+
154+
gg = txy -> ((t,x,y) = txy; g(x/t,y/t))
155+
V = gg.(p)
156+
P = plan_transform(DuffyCone(), V)
157+
@test P*V ≈ [0;1;Zeros(8)] ≈ Fun((t,x,y) -> g(x/t,y/t), DuffyCone(), 10).coefficients
158+
148159
f = Fun((t,x,y) -> x, DuffyCone(), 10)
149160
@test f(0.3,0.1,0.2) ≈ 0.1
150161

@@ -208,49 +219,40 @@ end
208219

209220

210221
@testset "LegendreConePlan" begin
211-
for N = 1:10
212-
n = N*sum(1:N)
213-
@test round(Int, 1/6*(1 + 1/(1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3) + (1 + 54n + 6sqrt(3)sqrt(n + 27n^2))^(1/3)), RoundUp)
214-
end
215222
p = points(LegendreCone(),10)
216223
@test length(p) == 12 == 2*2*3
217-
v = fill(1.0,length(p))
218-
n = length(v)
219-
M,N,O = pointsize(Cone(),n)
220-
221-
222-
D = plan_transform!(rectspace(DuffyCone()), reshape(v,N,M))
223-
N,M = D.plan[1][2],D.plan[2][2]
224-
V=reshape(v,N,M)
225-
@test D*V ≈ [[1/sqrt(2) zeros(1,5)]; zeros(2,6)]
224+
@test size(p) == (2,2,3)
225+
V = fill(1.0,size(p))
226+
M,_,N = size(V)
227+
D = plan_transform(DuffyCone(), V)
228+
C = tensortransform(D,V)
229+
@test C ≈ [[1/sqrt(2) zeros(1,2)]; zeros(1,3)]
226230
T = Float64
227-
P = c_plan_rottriangle(N, zero(T), zero(T), zero(T))
228-
duffy2legendrecone!(P,V)
229-
@test V ≈ [[1/sqrt(2) zeros(1,5)]; zeros(2,6)]
230-
@test fromtensor(LegendreCone(), V) ≈ [1/sqrt(2); Zeros(55)]
231+
P = c_plan_rottriangle(M, zero(T), zero(T), one(T))
232+
_coefficients!(P,C,DuffyCone(),LegendreCone())
233+
@test C ≈ [[0.8164965809277259 zeros(1,2)]; zeros(1,3)]
234+
@test fromtensor(LegendreCone(), C) ≈ [0.8164965809277259; Zeros(9)]
231235

232-
v = fill(1.0,length(p))
233-
P = plan_transform(LegendreCone(),v)
234-
@test P*v ≈ [1/sqrt(2); Zeros(55)]
235-
@test v == fill(1.0,length(p))
236+
P = plan_transform(LegendreCone(),V)
237+
@test P*V ≈ [0.8164965809277259; Zeros(9)]
238+
@test V == fill(1.0,size(p))
236239

237240
p = points(LegendreCone(),20)
238-
v = fill(1.0,length(p))
239-
P = plan_transform(LegendreCone(),v)
240-
241-
242-
for (m,k,ℓ) in ((0,0,0), )
243-
Y = Fun(ZernikeDisk(), [Zeros(sum(1:m)+k); 1])
244-
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))
245-
g = Fun(f, LegendreCone(),20)
246-
t,x,y = 0.3,0.1,0.2
247-
@test g(t,x,y) ≈ f((t,x,y))
248-
end
241+
V = fill(1.0,size(p))
242+
P = plan_transform(LegendreCone(),V)
243+
@test P*V ≈ [0.8164965809277259; Zeros(55)]
244+
end
245+
@testset "iduffy" begin
246+
p = points(LegendreCone(),20)
247+
V = fill(1.0,size(p))
248+
D = plan_transform(DuffyCone(),V)
249+
P = plan_transform(LegendreCone(),V)
250+
249251
end
250252

251253
@testset "LegendreCone" begin
252254
f = Fun((t,x,y) -> 1, LegendreCone(), 10)
253-
@test f.coefficients ≈ [1; zeros(ncoefficients(f)-1)]
255+
@test f.coefficients ≈ [0.8164965809277259; zeros(ncoefficients(f)-1)]
254256
@test f(0.3,0.1,0.2) ≈ 1
255257

256258
f = Fun((t,x,y) -> t, LegendreCone(), 10)
@@ -264,6 +266,16 @@ end
264266
f = Fun((t,x,y) -> 1+t+x+y, LegendreConic(), 1000)
265267
@test f(0.3,0.1,0.2) ≈ 1+0.3+0.1+0.2
266268
@test ncoefficients(f) == 1771
269+
270+
271+
for (m,k,ℓ) in ((0,0,0), (1,0,0), (0,0,1), (1,1,0), (2,1,1))
272+
Y = Fun(ZernikeDisk(), [Zeros(sum(1:m)+k); 1])
273+
f = (txy) -> ((t,x,y) = txy; Fun(NormalizedJacobi(0,2m+2,Segment(1,0)),[zeros(ℓ);1])(t) * 2^m * t^m * Y(x/t,y/t))
274+
g = Fun(f, LegendreCone(),100)
275+
@test sum(g.coefficients) ≈ 1
276+
t,x,y = 0.3,0.1,0.2
277+
@test g(t,x,y) ≈ f((t,x,y))
278+
end
267279
end
268280
end
269281

0 commit comments

Comments
 (0)