demo.py 4.04 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
import pyinterpol_ng
Martin Reinecke's avatar
Martin Reinecke committed
2
3
4
import numpy as np
import time

5
np.random.seed(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
12
13
def random_alm(lmax, mmax, ncomp):
    res = np.random.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
     + 1j*np.random.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
31
32
33
34
35
36
37
38
39

def convolve(alm1, alm2, lmax):
    job = pysharp.sharpjob_d()
    job.set_triangular_alm_info(lmax, lmax)
    job.set_gauss_geometry(lmax+1, 2*lmax+1)
    map = job.alm2map(alm1)*job.alm2map(alm2)
    job.set_triangular_alm_info(0,0)
    return job.map2alm(map)[0]*np.sqrt(4*np.pi)


Martin Reinecke's avatar
Martin Reinecke committed
40
lmax=1024
Martin Reinecke's avatar
Martin Reinecke committed
41
kmax=13
42
ncomp=1
Martin Reinecke's avatar
Martin Reinecke committed
43
44
separate=True
nptg = 50000000
45
epsilon = 1e-4
Martin Reinecke's avatar
Martin Reinecke committed
46
47
48
49
ofactor = 1.5
nthreads = 0  # use as many threads as available

ncomp2 = ncomp if separate else 1
Martin Reinecke's avatar
Martin Reinecke committed
50
51
52

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

Martin Reinecke's avatar
Martin Reinecke committed
55
# build beam a_lm
56
blm = random_alm(lmax, kmax, ncomp)
57
58
59


t0=time.time()
60
# build interpolator object for slm and blm
61
foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
Martin Reinecke's avatar
Martin Reinecke committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
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()
77
78
79
nth = lmax+1
nph = 2*lmax+1

Martin Reinecke's avatar
Martin Reinecke committed
80
ptg=np.random.uniform(0.,1.,3*nptg).reshape(nptg,3)
81
82
83
ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi
Martin Reinecke's avatar
Martin Reinecke committed
84

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

# build interpolator object for slm and blm
t0=time.time()
foo_f = pyinterpol_ng.PyInterpolator_f(slm.astype(np.complex64),blm.astype(np.complex64),separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
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))
foo2_f = pyinterpol_ng.PyInterpolator_f(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
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.))