totalconvolve_demo.py 4.57 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2020 Max-Planck-Society


Martin Reinecke's avatar
Martin Reinecke committed
17
import ducc0.totalconvolve as totalconvolve
Martin Reinecke's avatar
Martin Reinecke committed
18
19
20
import numpy as np
import time

21
rng = np.random.default_rng(48)
Martin Reinecke's avatar
Martin Reinecke committed
22

Martin Reinecke's avatar
Martin Reinecke committed
23

Martin Reinecke's avatar
Martin Reinecke committed
24
25
26
def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)

Martin Reinecke's avatar
Martin Reinecke committed
27

28
def random_alm(lmax, mmax, ncomp):
29
30
    res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
     + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
Martin Reinecke's avatar
Martin Reinecke committed
31
    # make a_lm with m==0 real-valued
Martin Reinecke's avatar
Martin Reinecke committed
32
    res[0:lmax+1, :].imag = 0.
Martin Reinecke's avatar
Martin Reinecke committed
33
34
    return res

35

Martin Reinecke's avatar
Martin Reinecke committed
36
def compress_alm(alm, lmax):
Martin Reinecke's avatar
Martin Reinecke committed
37
38
39
40
41
42
    res = np.empty(2*len(alm)-lmax-1, dtype=np.float64)
    res[0:lmax+1] = alm[0:lmax+1].real
    res[lmax+1::2] = np.sqrt(2)*alm[lmax+1:].real
    res[lmax+2::2] = np.sqrt(2)*alm[lmax+1:].imag
    return res

43

Martin Reinecke's avatar
Martin Reinecke committed
44
45
def myalmdot(a1, a2, lmax, mmax, spin):
    return np.vdot(compress_alm(a1, lmax), compress_alm(np.conj(a2), lmax))
Martin Reinecke's avatar
Martin Reinecke committed
46

47

Martin Reinecke's avatar
Martin Reinecke committed
48
49
50
51
lmax = 1024
kmax = 13
ncomp = 1
separate = True
Martin Reinecke's avatar
Martin Reinecke committed
52
nptg = 50000000
53
epsilon = 1e-4
Martin Reinecke's avatar
Martin Reinecke committed
54
55
56
57
ofactor = 1.5
nthreads = 0  # use as many threads as available

ncomp2 = ncomp if separate else 1
Martin Reinecke's avatar
Martin Reinecke committed
58
59
60

# get random sky a_lm
# the a_lm arrays follow the same conventions as those in healpy
61
slm = random_alm(lmax, lmax, ncomp)
Martin Reinecke's avatar
Martin Reinecke committed
62

Martin Reinecke's avatar
Martin Reinecke committed
63
# build beam a_lm
64
blm = random_alm(lmax, kmax, ncomp)
65
66


Martin Reinecke's avatar
Martin Reinecke committed
67
t0 = time.time()
68
# build interpolator object for slm and blm
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
69
70
71
foo = totalconvolve.Interpolator(slm, blm, separate, lmax, kmax,
                                 epsilon=epsilon, ofactor=ofactor,
                                 nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
72
73
t1 = time.time()-t0

Martin Reinecke's avatar
Martin Reinecke committed
74
print("Convolving sky and beam with lmax=mmax={}, kmax={}".format(lmax, kmax))
Martin Reinecke's avatar
Martin Reinecke committed
75
76
77
print("Interpolation taking place with a maximum error of {}\n"
      "and an oversampling factor of {}".format(epsilon, ofactor))
supp = foo.support()
Martin Reinecke's avatar
Martin Reinecke committed
78
print("(resulting in a kernel support size of {}x{})".format(supp, supp))
Martin Reinecke's avatar
Martin Reinecke committed
79
80
81
82
83
84
85
if ncomp == 1:
    print("One component")
else:
    print("{} components, which are {}coadded".format(ncomp, "not " if separate else ""))

print("\nDouble precision convolution/interpolation:")
print("preparation of interpolation grid: {}s".format(t1))
Martin Reinecke's avatar
Martin Reinecke committed
86
t0 = time.time()
87
88
89
nth = lmax+1
nph = 2*lmax+1

Martin Reinecke's avatar
Martin Reinecke committed
90
91
92
93
ptg = rng.uniform(0., 1., 3*nptg).reshape(nptg, 3)
ptg[:, 0] *= np.pi
ptg[:, 1] *= 2*np.pi
ptg[:, 2] *= 2*np.pi
Martin Reinecke's avatar
Martin Reinecke committed
94

Martin Reinecke's avatar
Martin Reinecke committed
95
96
t0 = time.time()
bar = foo.interpol(ptg)
Martin Reinecke's avatar
Martin Reinecke committed
97
del foo
Martin Reinecke's avatar
Martin Reinecke committed
98
99
100
print("Interpolating {} random angle triplets: {}s".format(nptg, time.time()-t0))
t0 = time.time()
fake = rng.uniform(0., 1., (ptg.shape[0], ncomp2))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
101
foo2 = totalconvolve.Interpolator(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
102
103
104
105
106
t0 = time.time()
foo2.deinterpol(ptg.reshape((-1, 3)), fake)
print("Adjoint interpolation: {}s".format(time.time()-t0))
t0 = time.time()
bla = foo2.getSlm(blm)
Martin Reinecke's avatar
Martin Reinecke committed
107
del foo2
Martin Reinecke's avatar
Martin Reinecke committed
108
109
110
print("Computing s_lm: {}s".format(time.time()-t0))
v1 = np.sum([myalmdot(slm[:, i], bla[:, i], lmax, lmax, 0) for i in range(ncomp)])
v2 = np.sum([np.vdot(fake[:, i], bar[:, i]) for i in range(ncomp2)])
Martin Reinecke's avatar
Martin Reinecke committed
111
112
113
print("Adjointness error: {}".format(v1/v2-1.))

# build interpolator object for slm and blm
Martin Reinecke's avatar
Martin Reinecke committed
114
t0 = time.time()
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
115
foo_f = totalconvolve.Interpolator_f(slm.astype(np.complex64), blm.astype(np.complex64), separate, lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
116
117
118
119
120
121
122
print("\nSingle precision convolution/interpolation:")
print("preparation of interpolation grid: {}s".format(time.time()-t0))

ptgf = ptg.astype(np.float32)
del ptg
fake_f = fake.astype(np.float32)
del fake
Martin Reinecke's avatar
Martin Reinecke committed
123
124
t0 = time.time()
bar_f = foo_f.interpol(ptgf)
Martin Reinecke's avatar
Martin Reinecke committed
125
del foo_f
Martin Reinecke's avatar
Martin Reinecke committed
126
print("Interpolating {} random angle triplets: {}s".format(nptg, time.time()-t0))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
127
foo2_f = totalconvolve.Interpolator_f(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
128
129
130
131
132
t0 = time.time()
foo2_f.deinterpol(ptgf.reshape((-1, 3)), fake_f)
print("Adjoint interpolation: {}s".format(time.time()-t0))
t0 = time.time()
bla_f = foo2_f.getSlm(blm.astype(np.complex64))
Martin Reinecke's avatar
Martin Reinecke committed
133
del foo2_f
Martin Reinecke's avatar
Martin Reinecke committed
134
135
136
print("Computing s_lm: {}s".format(time.time()-t0))
v1 = np.sum([myalmdot(slm[:, i], bla_f[:, i], lmax, lmax, 0) for i in range(ncomp)])
v2 = np.sum([np.vdot(fake_f[:, i], bar_f[:, i]) for i in range(ncomp2)])
Martin Reinecke's avatar
Martin Reinecke committed
137
print("Adjointness error: {}".format(v1/v2-1.))