1# This program is free software: you can redistribute it and/or modify
2# it under the terms of the GNU General Public License as published by
3# the Free Software Foundation, either version 3 of the License, or
4# (at your option) any later version.
5#
6# This program is distributed in the hope that it will be useful,
7# but WITHOUT ANY WARRANTY; without even the implied warranty of
8# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
9# GNU General Public License for more details.
10#
11# You should have received a copy of the GNU General Public License
12# along with this program.  If not, see <http://www.gnu.org/licenses/>.
13#
14# Copyright(C) 2019-2020 Max-Planck-Society
15
16
17import time
18import numpy as np
19import ducc0.healpix as ph
20
21rng = np.random.default_rng(42)
22
23
24def report(name, vlen, ntry, nside, isnest, perf):
25    print(name, ": ", perf*1e-6, "MOps/s", sep="")
26
27
28def random_ptg(vlen):
29    res = np.empty((vlen, 2), dtype=np.float64)
30    res[:, 0] = np.arccos((rng.random(vlen)-0.5)*2)
31    res[:, 1] = rng.random(vlen)*2*np.pi
32    return res
33
34
35def random_pix(nside, vlen):
36    return rng.integers(low=0, high=12*nside*nside-1, size=vlen, dtype=np.int64)
37
38
39def dummy(vlen):
40    _ = np.zeros(vlen, dtype=np.int64)
41
42
43def genperf(func, fname, inp, vlen, ntry, nside, isnest):
44    cnt = 0
45    t = time.time()
46    while cnt < ntry:
47        func(inp)
48        cnt += 1
49    t = time.time()-t
50    p = (vlen*ntry)/t
51    report(fname, vlen, ntry, nside, isnest, p)
52
53
54def perf_pix2ang(vlen, ntry, nside, isnest):
55    inp = random_pix(nside, vlen)
56    base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
57    genperf(base.pix2ang, "pix2ang", inp, vlen, ntry, nside, isnest)
58
59
60def perf_ang2pix(vlen, ntry, nside, isnest):
61    inp = random_ptg(vlen)
62    base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
63    genperf(base.ang2pix, "ang2pix", inp, vlen, ntry, nside, isnest)
64
65
66def perf_pix2vec(vlen, ntry, nside, isnest):
67    inp = random_pix(nside, vlen)
68    base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
69    genperf(base.pix2vec, "pix2vec", inp, vlen, ntry, nside, isnest)
70
71
72def perf_vec2pix(vlen, ntry, nside, isnest):
73    inp = ph.ang2vec(random_ptg(vlen))
74    base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
75    genperf(base.vec2pix, "vec2pix", inp, vlen, ntry, nside, isnest)
76
77
78def perf_ring2nest(vlen, ntry, nside, isnest):
79    inp = random_pix(nside, vlen)
80    base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
81    genperf(base.ring2nest, "ring2nest", inp, vlen, ntry, nside, isnest)
82
83
84def perf_nest2ring(vlen, ntry, nside, isnest):
85    inp = random_pix(nside, vlen)
86    base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
87    genperf(base.nest2ring, "nest2ring", inp, vlen, ntry, nside, isnest)
88
89
90def perf_neighbors(vlen, ntry, nside, isnest):
91    inp = random_pix(nside, vlen)
92    base = ph.Healpix_Base(nside, "NEST" if isnest else "RING")
93    genperf(base.neighbors, "neighbors", inp, vlen, ntry, nside, isnest)
94
95
96def suite(vlen, ntry, nside, isnest):
97    print("vlen=", vlen, ", ", "NEST" if isnest else "RING", sep="")
98    dummy(vlen)
99    perf_pix2ang(vlen, ntry, nside, isnest)
100    perf_ang2pix(vlen, ntry, nside, isnest)
101    perf_pix2vec(vlen, ntry, nside, isnest)
102    perf_vec2pix(vlen, ntry, nside, isnest)
103    perf_neighbors(vlen, ntry, nside, isnest)
104
105
106nside = 512
107ntry = 1000
108print("nside=", nside, sep="")
109for vlen in (1, 10, 100, 1000, 10000):
110    for isnest in (True, False):
111        suite(vlen, ntry, nside, isnest)
112        perf_ring2nest(vlen, ntry, nside, isnest)
113        perf_nest2ring(vlen, ntry, nside, isnest)
114        print()
115