totalconvolve_demo.py 3.78 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
import ducc_0_1.totalconvolve as totalconvolve
Martin Reinecke's avatar
Martin Reinecke committed
2
3
4
import numpy as np
import time

5
rng = np.random.default_rng(48)
Martin Reinecke's avatar
Martin Reinecke committed
6
7
8
9

def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)

Martin Reinecke's avatar
Martin Reinecke committed
10

11
def random_alm(lmax, mmax, ncomp):
12
13
    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
14
    # make a_lm with m==0 real-valued
15
    res[0:lmax+1,:].imag = 0.
Martin Reinecke's avatar
Martin Reinecke committed
16
17
    return res

18

Martin Reinecke's avatar
Martin Reinecke committed
19
20
21
22
23
24
25
def compress_alm(alm,lmax):
    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

26

Martin Reinecke's avatar
Martin Reinecke committed
27
def myalmdot(a1,a2,lmax,mmax,spin):
Martin Reinecke's avatar
cleanup  
Martin Reinecke committed
28
    return np.vdot(compress_alm(a1,lmax),compress_alm(np.conj(a2),lmax))
Martin Reinecke's avatar
Martin Reinecke committed
29

30

Martin Reinecke's avatar
Martin Reinecke committed
31
lmax=1024
Martin Reinecke's avatar
Martin Reinecke committed
32
kmax=13
33
ncomp=1
Martin Reinecke's avatar
Martin Reinecke committed
34
35
separate=True
nptg = 50000000
36
epsilon = 1e-4
Martin Reinecke's avatar
Martin Reinecke committed
37
38
39
40
ofactor = 1.5
nthreads = 0  # use as many threads as available

ncomp2 = ncomp if separate else 1
Martin Reinecke's avatar
Martin Reinecke committed
41
42
43

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

Martin Reinecke's avatar
Martin Reinecke committed
46
# build beam a_lm
47
blm = random_alm(lmax, kmax, ncomp)
48
49
50


t0=time.time()
51
# build interpolator object for slm and blm
Martin Reinecke's avatar
Martin Reinecke committed
52
foo = totalconvolve.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
t1 = time.time()-t0

print("Convolving sky and beam with lmax=mmax={}, kmax={}".format(lmax,kmax))
print("Interpolation taking place with a maximum error of {}\n"
      "and an oversampling factor of {}".format(epsilon, ofactor))
supp = foo.support()
print("(resulting in a kernel support size of {}x{})".format(supp,supp))
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))
t0=time.time()
68
69
70
nth = lmax+1
nph = 2*lmax+1

71
ptg=rng.uniform(0.,1.,3*nptg).reshape(nptg,3)
72
73
74
ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi
Martin Reinecke's avatar
Martin Reinecke committed
75

76
t0=time.time()
77
bar=foo.interpol(ptg)
Martin Reinecke's avatar
Martin Reinecke committed
78
del foo
Martin Reinecke's avatar
Martin Reinecke committed
79
80
print("Interpolating {} random angle triplets: {}s".format(nptg, time.time() -t0))
t0=time.time()
81
fake = rng.uniform(0.,1., (ptg.shape[0],ncomp2))
Martin Reinecke's avatar
Martin Reinecke committed
82
foo2 = totalconvolve.PyInterpolator(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
83
84
t0=time.time()
foo2.deinterpol(ptg.reshape((-1,3)), fake)
Martin Reinecke's avatar
Martin Reinecke committed
85
print("Adjoint interpolation: {}s".format(time.time() -t0))
86
87
t0=time.time()
bla=foo2.getSlm(blm)
Martin Reinecke's avatar
Martin Reinecke committed
88
89
del foo2
print("Computing s_lm: {}s".format(time.time() -t0))
90
91
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
92
93
94
95
print("Adjointness error: {}".format(v1/v2-1.))

# build interpolator object for slm and blm
t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
96
foo_f = totalconvolve.PyInterpolator_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
97
98
99
100
101
102
103
104
105
106
107
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
t0=time.time()
bar_f=foo_f.interpol(ptgf)
del foo_f
print("Interpolating {} random angle triplets: {}s".format(nptg, time.time() -t0))
Martin Reinecke's avatar
Martin Reinecke committed
108
foo2_f = totalconvolve.PyInterpolator_f(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
109
110
111
112
113
114
115
116
117
118
119
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))
del foo2_f
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)])
print("Adjointness error: {}".format(v1/v2-1.))