-
Notifications
You must be signed in to change notification settings - Fork 103
/
TimingAccuracyDH.py
executable file
·79 lines (64 loc) · 2.42 KB
/
TimingAccuracyDH.py
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
"""
This script is a python version of TimingAccuracy. We use some numpy functions
to simplify the creation of random coefficients.
"""
import time
import numpy as np
import pyshtools as pysh
# ==== MAIN FUNCTION ====
def main():
TimingAccuracyDH(2)
# ==== TEST FUNCTIONS ====
def TimingAccuracyDH(sampling=1):
# ---- input parameters ----
maxdeg = 2800
ls = np.arange(maxdeg + 1)
beta = 1.5
print('Driscoll-Healy (real) sampling =', sampling)
# ---- create mask to filter out m<=l ----
mask = np.zeros((2, maxdeg + 1, maxdeg + 1), dtype=bool)
mask[0, 0, 0] = True
for l in ls:
mask[:, l, :l + 1] = True
mask[1, :, 0] = False
# ---- create Gaussian powerlaw coefficients ----
print('creating {:d} random coefficients'.format(2 * (maxdeg + 1) *
(maxdeg + 1)))
np.random.seed(0)
cilm = np.random.normal(loc=0., scale=1., size=(2, maxdeg + 1, maxdeg + 1))
old_power = pysh.spectralanalysis.spectrum(cilm)
new_power = 1. / (1. + ls)**beta # initialize degrees > 0 to power-law
cilm[:, :, :] *= np.sqrt(new_power / old_power)[None, :, None]
cilm[~mask] = 0.
# ---- time spherical harmonics transform for lmax set to increasing
# ---- powers of 2.
lmax = 2
print('lmax maxerror rms tinverse tforward')
while lmax <= maxdeg:
# trim coefficients to lmax
cilm_trim = cilm[:, :lmax + 1, :lmax + 1]
mask_trim = mask[:, :lmax + 1, :lmax + 1]
# synthesis / inverse
tstart = time.time()
grid = pysh.expand.MakeGridDH(cilm_trim, sampling=sampling)
tend = time.time()
tinverse = tend - tstart
# analysis / forward
tstart = time.time()
cilm2_trim = pysh.expand.SHExpandDH(grid, sampling=sampling)
tend = time.time()
tforward = tend - tstart
# compute error
err = np.abs(cilm_trim[mask_trim] - cilm2_trim[mask_trim]) / \
np.abs(cilm_trim[mask_trim])
maxerr = err.max()
rmserr = np.mean(err**2)
print('{:4d} {:1.2e} {:1.2e} {:1.1e}s {:1.1e}s'.format(
lmax, maxerr, rmserr, tinverse, tforward))
if maxerr > 100.:
raise RuntimeError('Tests Failed. Maximum relative error = ',
maxerr)
lmax = lmax * 2
# ==== EXECUTE SCRIPT ====
if __name__ == "__main__":
main()