-
Notifications
You must be signed in to change notification settings - Fork 1
/
runtests.jl
315 lines (255 loc) · 13.3 KB
/
runtests.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
using InverseHankelFunction
import Test: @test, @testset, @test_throws
import SpecialFunctions: hankelh1
import BenchmarkTools: @benchmark
include("utils.jl")
@testset "Hankel function asymptotics " begin
@testset "Small argument" begin
# Some small values to test our small argument asymptotics
small_values = 1e-4 .*[1+0im, 1im, -1+0im, -1im]
@testset "H ≈ H_asymptotic" begin
@test small_arg_hankelh1.(0, small_values) ≈ hankelh1.(0, small_values) rtol=0.1
@test small_arg_hankelh1.(1, small_values) ≈ hankelh1.(1, small_values) rtol=1e-5
@test small_arg_hankelh1.(2, small_values) ≈ hankelh1.(2, small_values) rtol=1e-5
@test small_arg_hankelh1.(4, small_values) ≈ hankelh1.(4, small_values) rtol=1e-5
@test small_arg_hankelh1.(8, small_values) ≈ hankelh1.(8, small_values) rtol=1e-5
@test small_arg_hankelh1.(11, small_values) ≈ hankelh1.(11, small_values) rtol=1e-5
end
# Test whether f f⁻¹ = identity for different Hankel order (n) and solution branches (b)
# We apply f⁻¹ first because action of f elliminates effect of branch
@testset "f f⁻¹(z) = z" begin
for n in [0,1,2,4,8,11]
for b in [0,1,2,3,4,5,10,11]
f = z-> small_arg_hankelh1(n, z)
invf = z->inv_small_arg_hankelh1(n, z, b)
# f is small argument asymptotic, so f(small_values) should be in domain of invf
testvec = f.(small_values)
if !all(isapprox.((f ∘ invf).(testvec), testvec; rtol=1e-5))
@show n, b, abs.(((f ∘ invf).(testvec) .- testvec)./testvec)
end
@test (f ∘ invf).(testvec) ≈ testvec rtol=1e-5
end
end
end
@testset "Domain enforcement for ν=0" begin
@test_throws DomainError inv_small_arg_hankelh1(0, 3.1+0.1im, 6)
@test_throws DomainError inv_small_arg_hankelh1(0, -1.01-0.7im, 9)
end
end
@testset "Large argument" begin
# Some large values to test our large argument asymptotics
large_values = [1.0e8+0im, 1.0e8+0.5im, 1.0e8-0.5im]
@testset "H ≈ H_asymptotic" begin
@test large_arg_hankelh1.(0, large_values) ≈ hankelh1.(0, large_values) rtol=1e-4
@test large_arg_hankelh1.(1, large_values) ≈ hankelh1.(1, large_values) rtol=1e-4
@test large_arg_hankelh1.(2, large_values) ≈ hankelh1.(2, large_values) rtol=1e-4
@test large_arg_hankelh1.(3, large_values) ≈ hankelh1.(3, large_values) rtol=1e-4
@test large_arg_hankelh1.(4, large_values) ≈ hankelh1.(4, large_values) rtol=1e-4
@test large_arg_hankelh1.(8, large_values) ≈ hankelh1.(8, large_values) rtol=1e-4
end
# Test whether f f⁻¹ = identity for different Hankel order (n) and solution branches (b)
# We apply f⁻¹ first because action of f elliminates effect of branch
@testset "f f⁻¹(z) = z" begin
for n in [0,1,2,4,8,11]
for b in [0,1,2,3]
f = z-> large_arg_hankelh1(n, z)
invf = z->inv_large_arg_hankelh1(n, z, b)
# f is large argument asymptotic, so f(large_values) should be in domain of invf
testvec = f.(large_values)
diffvec = abs.((f ∘ invf).(testvec) .- testvec)
if !all(isapprox.((f ∘ invf).(testvec), testvec; rtol=1e-5))
@show n, b, testvec, invf.(testvec), (f ∘ invf).(testvec)
end
@test (f ∘ invf).(testvec) ≈ testvec rtol=1e-5
end
end
end
end
@testset "Asymptotic scale" begin
# Small |z| is small
@test hankel_arg_asymptotic_scale(0, 0.001 ) < 1/2
@test hankel_arg_asymptotic_scale(4, -0.001 ) < 1/2
@test hankel_arg_asymptotic_scale(5, 0.001im) < 1/2
# Big |z| is big
@test hankel_arg_asymptotic_scale(0, 10000 ) > 2
@test hankel_arg_asymptotic_scale(9, -10000 ) > 2
@test hankel_arg_asymptotic_scale(3, -1000im) > 2
# Decreasing in |ν|
@test 2 < hankel_arg_asymptotic_scale(0 , 3)
@test 1/2 < hankel_arg_asymptotic_scale( 20, 3) < 2
@test hankel_arg_asymptotic_scale(450, 3) < 1/2
@test 2 < hankel_arg_asymptotic_scale(0 , -3)
@test 1/2 < hankel_arg_asymptotic_scale(-20, 3) < 2
@test hankel_arg_asymptotic_scale(450im, 3) < 1/2
# Although |z| dominates for large |z|
@test hankel_arg_asymptotic_scale( 100, 100 ) > 2
@test hankel_arg_asymptotic_scale( 100, -100 ) > 2
@test hankel_arg_asymptotic_scale(-100, -100im) > 2
@test hankel_arg_asymptotic_scale( 100, -100im) > 2
end
end
@testset "Normalised Hankel function " begin
ν = 1
z₀ = 2.0 - 0.3im
z = 2.0 + 0.1im
# Construct the partially evaluated version, which caches the value of hankel at z₀
hn1_partial = hankelh1n(1, z₀)
# Should be equal when we evaluate in stages or at once
@test hn1_partial(z) ≈ hankelh1n(ν, z₀, z)
# Check it isn't trivial
@test hn1_partial(z) != 1
@test hankelh1n(ν, z₀, z) != 1
# Normalised function evaluated at z₀ will always be 1.0
@test hn1_partial(z₀) ≈ 1
@test hankelh1n(ν, z₀, z₀) ≈ 1
# Cached version should be faster
b1 = @benchmark hankelh1n($ν, $z₀, $z) samples=20 evals=1
b2 = @benchmark $hn1_partial($z) samples=20 evals=1
@test minimum(b2.times) <= minimum(b1.times)
end
@testset "Inverse Normalised Hankel function" begin
@testset "Vector optimisation" begin
hns = 1.0:-0.01:0.01
hns = 1.0:-0.1:0.1
for z₀ in [1.0, 2.0, 4.0]
for ν in [0, 1, 4]
begin
# Test that the optimised and unoptimised versions give the same result
naive_method = invhankelh1n.(ν, z₀, hns)
@test naive_method ≈ invhankelh1n_sortedvec(ν, z₀, hns)
@test naive_method ≈ invhankelh1n(ν, z₀, hns)
# Test that vector optimised versions are faster (note no broadcasting on b3)
b1 = @benchmark invhankelh1n.($ν, $z₀, $hns) samples=10 evals=1
b2 = @benchmark invhankelh1n_sortedvec($ν, $z₀, $hns) samples=10 evals=1
b3 = @benchmark invhankelh1n($ν, $z₀, $hns) samples=10 evals=1
@test minimum(b2.times) <= minimum(b1.times)
@test minimum(b3.times) <= minimum(b1.times)
end
begin # diff version
# Test that the optimised and unoptimised versions give the same result
unzip(arr) = ([a[1] for a in arr], [a[2] for a in arr])
inv, diff = unzip(diffinvhankelh1n.(ν, z₀, hns))
inv_vec, diff_vec = diffinvhankelh1n_sortedvec(ν, z₀, hns)
@test inv ≈ inv_vec
@test diff ≈ diff_vec
# Test that optimised version is faster
b1 = @benchmark diffinvhankelh1n.($ν, $z₀, $hns) samples=10 evals=1
b2 = @benchmark diffinvhankelh1n_sortedvec($ν, $z₀, $hns) samples=10 evals=1
@test minimum(b2.times) <= minimum(b1.times)
end
end
end
end
@testset "Unsorted Vec function" begin
# Deliberately unsorted vec
hns = [0.4, 0.1, 1.1, 0.9, 1.7]
z₀ = 1.3
ν = 4
naive_method = invhankelh1n.(ν, z₀, hns)
@test naive_method ≈ invhankelh1n(ν, z₀, hns)
end
@testset "Householder corrector orders" begin
# Collect the total number of predictor and corrector steps so that we can assess performance
total_predictor_steps = 0
total_corrector_steps1 = 0
total_corrector_steps2 = 0
total_corrector_steps3 = 0
for ν in [0, 1, 4]
for z₀ in [1.0, 2.0, 4.0]
for ξ_target in [0.9, 0.34, 0.01]
iters1 = []; iters2 = []; iters3 =[];
ans1 = invhankelh1n_adaptive_solve(ν, z₀, ξ_target; householder_order=1, iter_vec=iters1)
ans2 = invhankelh1n_adaptive_solve(ν, z₀, ξ_target; householder_order=2, iter_vec=iters2)
ans3 = invhankelh1n_adaptive_solve(ν, z₀, ξ_target; householder_order=3, iter_vec=iters3)
# All methods should give similar answers
@test all(ans1 .≈ ans2)
@test all(ans2 .≈ ans3)
@test all(ans3 .≈ ans1)
# Higher order methods should require no more steps
@test all(iters2 .<= iters1)
@test all(iters3 .<= iters2)
# All methods should use the same predictor
@test length(iters1) == length(iters2) && length(iters2) == length(iters3)
total_predictor_steps += length(iters1)
total_corrector_steps1 += sum(iters1)
total_corrector_steps2 += sum(iters2)
total_corrector_steps3 += sum(iters3)
# Unimplemented/invalid householder corrector order
@test_throws Exception invhankelh1n_adaptive_solve(ν, z₀, ξ_target; householder_order=0)
@test_throws Exception invhankelh1n_adaptive_solve(ν, z₀, ξ_target; householder_order=4)
end
end
end
# Ratios of corrector to predictor steps for each method order
cp_ratio1 = total_corrector_steps1/total_predictor_steps
cp_ratio2 = total_corrector_steps2/total_predictor_steps
cp_ratio3 = total_corrector_steps3/total_predictor_steps
# Higher order methods should require less steps
@test cp_ratio3 < cp_ratio2 < cp_ratio1
# Empirically derived differences in steps, these may change if the examples change
@test cp_ratio2 <= cp_ratio1*0.7
@test cp_ratio3 <= cp_ratio2*0.85
end
@testset "Too many steps kill switch" begin
@test typeof(invhankelh1n_adaptive_solve(0, 1.0, 0; silent_failure=true)) <: Tuple
@test_throws Exception invhankelh1n_adaptive_solve(0, 1.0, 0; silent_failure=false)
end
end
@testset "Passing through " begin
for ν in [0, 1, 4]
for p in PassingThrough.([1.0, 2.0, 4.0, 1.0+0.5im, 2.0-0.2im])
h(z) = hankelh1(ν, z)
h⁻¹(z) = invhankelh1(ν, z, p)
# Check that invhankelh1 does indeed pass through this point
@test h⁻¹(h(p.point)) ≈ p.point
for hz in [0.9, 0.34, 0.01, 1.4, 1.0-0.5im]
h⁻¹hz = h⁻¹(hz)
@test h(h⁻¹hz) ≈ hz
end
end
end
@testset "Domain errors" begin
ν = 4
h = 3.4 + 0.1im
# NaNs
@test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN - 1.0*im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN - Inf*im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN - NaN*im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough(10 - NaN*im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough(NaN))
# Any kind of infinity
@test_throws DomainError invhankelh1(ν, h, PassingThrough(Inf + 0im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough(1.0 + Inf*im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough(Inf + -Inf*im))
# Zeros
@test_throws DomainError invhankelh1(ν, h, PassingThrough( 0.0 + 0.0*im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough( 0.0 - 0.0*im))
@test_throws DomainError invhankelh1(ν, h, PassingThrough(-0.0 - 0.0*im))
end
end
@testset "Derivatives of Hankel functions " begin
for ν in [0, 1, 4, 10]
for z in [1.0+0.0im, -1.0+0.1im, 0.0+2.0im, 0.0-2.0im, 100.0+5.0im]
findiff(f, δ=1e-8) = (f(ν, z+δ) .- f(ν, z))./δ
# Cache values
h = hankelh1(ν, z)
hm1 = hankelh1(ν-1, z)
hp1 = diffhankelh1( ν, z, h, hm1)
(hp2, hpp2) = diff2hankelh1(ν, z, h, hm1)
(hp3, hpp3, hppp3) = diff3hankelh1(ν, z, h, hm1)
# Test they are consistent
@test hp1 ≈ hp2
@test hp2 ≈ hp3
@test hp3 ≈ hp1
@test hpp2 ≈ hpp3
# Test optional arguments don't matter
@test diffhankelh1(ν, z, h, hm1) == diffhankelh1(ν, z)
@test diff2hankelh1(ν, z, h, hm1) == diff2hankelh1(ν, z)
@test diff3hankelh1(ν, z, h, hm1) == diff3hankelh1(ν, z)
# Test the derivatives are similar to naive finite difference
@test diffhankelh1(ν,z)[1] ≈ findiff(hankelh1) rtol=1e-6
@test diff2hankelh1(ν,z)[2] ≈ findiff(diffhankelh1)[1] rtol=1e-6
@test diff3hankelh1(ν,z)[3] ≈ findiff(diff2hankelh1)[2] rtol=1e-6
end
end
end