Skip to content

Commit 4d4751b

Browse files
dlfivefiftydlfivefifty
committed
empty transforms, increase coverage
1 parent c35ccd1 commit 4d4751b

File tree

2 files changed

+136
-76
lines changed

2 files changed

+136
-76
lines changed

‎src/chebyshevtransform.jl

Lines changed: 54 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -17,13 +17,10 @@ ChebyshevTransformPlan{k,inp}(plan) where {k,inp} =
1717

1818
function plan_chebyshevtransform!(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
1919
if kind == 1
20-
plan = plan_r2r!(x, REDFT10)
20+
plan = isempty(x) ? fill(one(T),1,length(x)) : plan_r2r!(x, REDFT10)
2121
ChebyshevTransformPlan{1,true}(plan)
2222
elseif kind == 2
23-
if length(x) ≤ 1
24-
error("Cannot create a length $(length(x)) chebyshev transform")
25-
end
26-
plan = plan_r2r!(x, REDFT00)
23+
plan = length(x) ≤ 1 ? fill(one(T),1,length(x)) : plan_r2r!(x, REDFT00)
2724
ChebyshevTransformPlan{2,true}(plan)
2825
end
2926
end
@@ -35,29 +32,20 @@ end
3532

3633
function *(P::ChebyshevTransformPlan{T,1,true},x::AbstractVector{T}) where T
3734
n = length(x)
38-
if n == 1
39-
x
40-
else
41-
x = P.plan*x
42-
x[1]/=2
43-
lmul!(inv(convert(T,n)), x)
44-
end
35+
n ≤ 1 && return x
36+
37+
x = P.plan*x
38+
x[1] /= 2
39+
lmul!(inv(convert(T,n)), x)
4540
end
4641

47-
function *(P::ChebyshevTransformPlan{T,2,true},x::AbstractVector{T}) where T
42+
function *(P::ChebyshevTransformPlan{T,2,true}, x::AbstractVector{T}) where T
4843
n = length(x)
49-
if n == 1
50-
x
51-
else
52-
n = length(x)
53-
if n == 1
54-
x
55-
else
56-
x = P.plan*x
57-
x[1] /= 2;x[end] /= 2
58-
lmul!(inv(convert(T,n-1)),x)
59-
end
60-
end
44+
n ≤ 1 && return x
45+
46+
x = P.plan*x
47+
x[1] /= 2; x[end] /= 2
48+
lmul!(inv(convert(T,n-1)),x)
6149
end
6250

6351
chebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
@@ -88,15 +76,9 @@ inv(P::IChebyshevTransformPlan{T,2,true}) where T = P.plan
8876

8977
function plan_ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
9078
if kind == 1
91-
if length(x) == 0
92-
error("Cannot create a length 0 inverse chebyshev transform")
93-
end
94-
plan = plan_r2r!(x, REDFT01)
79+
plan = isempty(x) ? fill(one(T),1,length(x)) : plan_r2r!(x, REDFT01)
9580
IChebyshevTransformPlan{T,1,true,typeof(plan)}(plan)
9681
elseif kind == 2
97-
if length(x) ≤ 1
98-
error("Cannot create a length $(length(x)) inverse chebyshev transform")
99-
end
10082
inv(plan_chebyshevtransform!(x;kind=2))
10183
end
10284
end
@@ -107,28 +89,25 @@ function plan_ichebyshevtransform(x::AbstractVector{T};kind::Integer=1) where T<
10789
end
10890

10991
function *(P::IChebyshevTransformPlan{T,1,true},x::AbstractVector{T}) where T<:fftwNumber
92+
isempty(x) && return x
11093
x[1] *=2
111-
x = lmul!(convert(T,0.5),P.plan*x)
94+
x = lmul!(convert(T,0.5), P.plan*x)
11295
x
11396
end
11497

11598
function *(P::IChebyshevTransformPlan{T,2,true},x::AbstractVector{T}) where T<:fftwNumber
11699
n = length(x)
117-
if n == 1
118-
x
119-
else
120-
##TODO: make thread safe
121-
x[1] *= 2;x[end] *= 2
122-
x = P.plan*x
123-
x[1] *= 2;x[end] *= 2
124-
lmul!(convert(T,.5(n-1)),x)
125-
end
100+
n ≤ 1 && return x
101+
x[1] *= 2; x[end] *= 2
102+
x = P.plan*x
103+
x[1] *= 2; x[end] *= 2
104+
lmul!(convert(T,0.5(n-1)),x)
126105
end
127106

128107
ichebyshevtransform!(x::AbstractVector{T};kind::Integer=1) where {T<:fftwNumber} =
129108
plan_ichebyshevtransform!(x;kind=kind)*x
130109

131-
ichebyshevtransform(x;kind::Integer=1) = ichebyshevtransform!(copy(x);kind=kind)
110+
ichebyshevtransform(x;kind::Integer=1) = ichebyshevtransform!(copy(x); kind=kind)
132111

133112
*(P::IChebyshevTransformPlan{T,k,false},x::AbstractVector{T}) where {T,k} = P.plan*copy(x)
134113

@@ -142,7 +121,7 @@ end
142121
# Matrix inputs
143122

144123

145-
function chebyshevtransform!(X::AbstractMatrix{T};kind::Integer=1) where T<:fftwNumber
124+
function chebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
146125
if kind == 1
147126
if size(X) == (1,1)
148127
X
@@ -164,7 +143,7 @@ function chebyshevtransform!(X::AbstractMatrix{T};kind::Integer=1) where T<:fftw
164143
end
165144
end
166145

167-
function ichebyshevtransform!(X::AbstractMatrix{T};kind::Integer=1) where T<:fftwNumber
146+
function ichebyshevtransform!(X::AbstractMatrix{T}; kind::Integer=1) where T<:fftwNumber
168147
if kind == 1
169148
if size(X) == (1,1)
170149
X
@@ -200,13 +179,10 @@ ChebyshevUTransformPlan{k,inp}(plan) where {k,inp} =
200179

201180
function plan_chebyshevutransform!(x::AbstractVector{T}; kind::Integer=1) where T<:fftwNumber
202181
if kind == 1
203-
plan = plan_r2r!(x, RODFT10)
182+
plan = isempty(x) ? fill(one(T),1,length(x)) : plan_r2r!(x, RODFT10)
204183
ChebyshevUTransformPlan{1,true}(plan)
205184
elseif kind == 2
206-
if length(x) ≤ 1
207-
error("Cannot create a length $(length(x)) chebyshevu transform")
208-
end
209-
plan = plan_r2r!(x, RODFT00)
185+
plan = length(x) ≤ 1 ? fill(one(T),1,length(x)) : plan_r2r!(x, RODFT00)
210186
ChebyshevUTransformPlan{2,true}(plan)
211187
end
212188
end
@@ -218,6 +194,8 @@ end
218194

219195
function *(P::ChebyshevUTransformPlan{T,1,true},x::AbstractVector{T}) where T
220196
n = length(x)
197+
n ≤ 1 && return x
198+
221199
for k=1:n # sqrt(1-x_j^2) weight
222200
x[k] *= sinpi(one(T)/(2n) + (k-one(T))/n)/n
223201
end
@@ -226,6 +204,8 @@ end
226204

227205
function *(P::ChebyshevUTransformPlan{T,2,true},x::AbstractVector{T}) where T
228206
n = length(x)
207+
n ≤ 1 && return x
208+
229209
c = one(T)/ (n+1)
230210
for k=1:n # sqrt(1-x_j^2) weight
231211
x[k] *= sinpi(k*c)
@@ -247,18 +227,13 @@ struct IChebyshevUTransformPlan{T,kind,inplace,P}
247227
plan::P
248228
end
249229

230+
250231
function plan_ichebyshevutransform!(x::AbstractVector{T};kind::Integer=1) where T<:fftwNumber
251232
if kind == 1
252-
if length(x) == 0
253-
error("Cannot create a length 0 inverse chebyshevu transform")
254-
end
255-
plan = plan_r2r!(x, RODFT01)
233+
plan = isempty(x) ? fill(one(T),1,length(x)) : plan_r2r!(x, RODFT01)
256234
IChebyshevUTransformPlan{T,1,true,typeof(plan)}(plan)
257235
elseif kind == 2
258-
if length(x) ≤ 1
259-
error("Cannot create a length $(length(x)) inverse chebyshevu transform")
260-
end
261-
plan = plan_r2r!(x, RODFT00)
236+
plan = length(x) ≤ 1 ? fill(one(T),1,length(x)) : plan_r2r!(x, RODFT00)
262237
IChebyshevUTransformPlan{T,2,true,typeof(plan)}(plan)
263238
end
264239
end
@@ -270,6 +245,8 @@ end
270245

271246
function *(P::IChebyshevUTransformPlan{T,1,true}, x::AbstractVector{T}) where T<:fftwNumber
272247
n = length(x)
248+
n ≤ 1 && return x
249+
273250
x = P.plan * x
274251
for k=1:n # sqrt(1-x_j^2) weight
275252
x[k] /= 2sinpi(one(T)/(2n) + (k-one(T))/n)
@@ -281,6 +258,8 @@ end
281258

282259
function *(P::IChebyshevUTransformPlan{T,2,true}, x::AbstractVector{T}) where T<:fftwNumber
283260
n = length(x)
261+
n ≤ 1 && return x
262+
284263
c = one(T)/ (n+1)
285264
lmul!((n+1)/(2n+2*one(T)), x)
286265
x = P.plan * x
@@ -312,7 +291,7 @@ function chebyshevpoints(::Type{T}, n::Integer; kind::Int=1) where T<:Number
312291
if kind == 1
313292
T[sinpi((n-2k-one(T))/2n) for k=0:n-1]
314293
elseif kind == 2
315-
if n==1
294+
if n == 1
316295
zeros(T,1)
317296
else
318297
T[cospi(k/(n-one(T))) for k=0:n-1]
@@ -323,19 +302,19 @@ chebyshevpoints(n::Integer; kind::Int=1) = chebyshevpoints(Float64, n; kind=kind
323302

324303

325304
# sin(nθ) coefficients to values at Clenshaw-Curtis nodes except ±1
326-
327-
struct DSTPlan{T,kind,inplace,P} <: Plan{T}
328-
plan::P
329-
end
330-
331-
DSTPlan{k,inp}(plan) where {k,inp} =
332-
DSTPlan{eltype(plan),k,inp,typeof(plan)}(plan)
333-
334-
335-
plan_DSTI!(x) = length(x) > 0 ? DSTPlan{1,true}(FFTW.plan_r2r!(x, FFTW.RODFT00)) :
336-
ones(x)'
337-
338-
function *(P::DSTPlan{T,1}, x::AbstractArray) where {T}
339-
x = P.plan*x
340-
rmul!(x,half(T))
341-
end
305+
#
306+
# struct DSTPlan{T,kind,inplace,P} <: Plan{T}
307+
# plan::P
308+
# end
309+
#
310+
# DSTPlan{k,inp}(plan) where {k,inp} =
311+
# DSTPlan{eltype(plan),k,inp,typeof(plan)}(plan)
312+
#
313+
#
314+
# plan_DSTI!(x) = length(x) > 0 ? DSTPlan{1,true}(FFTW.plan_r2r!(x, FFTW.RODFT00)) :
315+
# fill(one(T),1,length(x))
316+
#
317+
# function *(P::DSTPlan{T,1}, x::AbstractArray) where {T}
318+
# x = P.plan*x
319+
# rmul!(x,half(T))
320+
# end

‎test/chebyshevtests.jl

Lines changed: 82 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,26 @@ using FastTransforms, Compat, Compat.Test
1111
f̃ = x -> [cos(k*acos(x)) for k=0:n-1]' * f̌
1212
@test f̃(0.1) ≈ exp(T(0.1))
1313
@test ichebyshevtransform(f̌; kind=1) ≈ exp.(p_1)
14+
15+
f̃ = copy(f)
16+
f̄ = copy(f̌)
17+
P = plan_chebyshevtransform(f; kind=1)
18+
@test P*f == f̌
19+
@test f == f̃
20+
P = plan_chebyshevtransform!(f; kind=1)
21+
@test P*f == f̌
22+
@test f == f̌
23+
Pi = plan_ichebyshevtransform(f̌; kind=1)
24+
@test Pi*f̌ ≈ f̃
25+
@test f̌ == f̄
26+
Pi = plan_ichebyshevtransform!(f̌; kind=1)
27+
@test Pi*f̌ ≈ f̃
28+
@test f̌ ≈ f̃
29+
30+
@test chebyshevtransform(T[1]; kind=1) == T[1]
31+
@test ichebyshevtransform(T[1]; kind=1) == T[1]
32+
@test chebyshevtransform(T[]; kind=1) == T[]
33+
@test ichebyshevtransform(T[]; kind=1) == T[]
1434
end
1535
end
1636
@testset "Chebyshev second kind points <-> first kind coefficients" begin
@@ -27,7 +47,28 @@ using FastTransforms, Compat, Compat.Test
2747
P = plan_chebyshevtransform!(f; kind=2)
2848
Pi = plan_ichebyshevtransform!(f; kind=2)
2949
@test all((P \ copy(f)) .=== Pi * copy(f))
30-
@test f ≈ P \ (P*copy(f)) ≈ P * (P\copy(f))
50+
@test all((Pi \ copy(f̌)) .=== P * copy(f̌))
51+
@test f ≈ P \ (P*copy(f)) ≈ P * (P\copy(f)) ≈ Pi \ (Pi*copy(f)) ≈ Pi * (Pi \ copy(f))
52+
53+
f̃ = copy(f)
54+
f̄ = copy(f̌)
55+
P = plan_chebyshevtransform(f; kind=2)
56+
@test P*f == f̌
57+
@test f == f̃
58+
P = plan_chebyshevtransform!(f; kind=2)
59+
@test P*f == f̌
60+
@test f == f̌
61+
Pi = plan_ichebyshevtransform(f̌; kind=2)
62+
@test Pi*f̌ ≈ f̃
63+
@test f̌ == f̄
64+
Pi = plan_ichebyshevtransform!(f̌; kind=2)
65+
@test Pi*f̌ ≈ f̃
66+
@test f̌ ≈ f̃
67+
68+
@test chebyshevtransform(T[1]; kind=2) == T[1]
69+
@test ichebyshevtransform(T[1]; kind=2) == T[1]
70+
@test chebyshevtransform(T[]; kind=2) == T[]
71+
@test ichebyshevtransform(T[]; kind=2) == T[]
3172
end
3273
end
3374

@@ -41,6 +82,26 @@ using FastTransforms, Compat, Compat.Test
4182
f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-1]' * f̌
4283
@test f̃(0.1) ≈ exp(T(0.1))
4384
@test ichebyshevutransform(f̌; kind=1) ≈ exp.(p_1)
85+
86+
f̃ = copy(f)
87+
f̄ = copy(f̌)
88+
P = plan_chebyshevutransform(f; kind=1)
89+
@test P*f == f̌
90+
@test f == f̃
91+
P = plan_chebyshevutransform!(f; kind=1)
92+
@test P*f == f̌
93+
@test f == f̌
94+
Pi = plan_ichebyshevutransform(f̌; kind=1)
95+
@test Pi*f̌ ≈ f̃
96+
@test f̌ == f̄
97+
Pi = plan_ichebyshevutransform!(f̌; kind=1)
98+
@test Pi*f̌ ≈ f̃
99+
@test f̌ ≈ f̃
100+
101+
@test chebyshevutransform(T[1]; kind=1) == T[1]
102+
@test ichebyshevutransform(T[1]; kind=1) == T[1]
103+
@test chebyshevutransform(T[]; kind=1) == T[]
104+
@test ichebyshevutransform(T[]; kind=1) == T[]
44105
end
45106
end
46107

@@ -54,6 +115,26 @@ using FastTransforms, Compat, Compat.Test
54115
f̃ = x -> [sin((k+1)*acos(x))/sin(acos(x)) for k=0:n-3]' * f̌
55116
@test f̃(0.1) ≈ exp(T(0.1))
56117
@test ichebyshevutransform(f̌; kind=2) ≈ exp.(p_2)
118+
119+
f̃ = copy(f)
120+
f̄ = copy(f̌)
121+
P = plan_chebyshevutransform(f; kind=2)
122+
@test P*f == f̌
123+
@test f == f̃
124+
P = plan_chebyshevutransform!(f; kind=2)
125+
@test P*f == f̌
126+
@test f == f̌
127+
Pi = plan_ichebyshevutransform(f̌; kind=2)
128+
@test Pi*f̌ ≈ f̃
129+
@test f̌ == f̄
130+
Pi = plan_ichebyshevutransform!(f̌; kind=2)
131+
@test Pi*f̌ ≈ f̃
132+
@test f̌ ≈ f̃
133+
134+
@test chebyshevutransform(T[1]; kind=2) == T[1]
135+
@test ichebyshevutransform(T[1]; kind=2) == T[1]
136+
@test chebyshevutransform(T[]; kind=2) == T[]
137+
@test ichebyshevutransform(T[]; kind=2) == T[]
57138
end
58139
end
59140
end

0 commit comments

Comments
 (0)