demo.py 2.94 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
import interpol_ng
import numpy as np
import pysharp
import time
Martin Reinecke's avatar
Martin Reinecke committed
5
import matplotlib.pyplot as plt
Martin Reinecke's avatar
Martin Reinecke committed
6

Martin Reinecke's avatar
Martin Reinecke committed
7
np.random.seed(20)
Martin Reinecke's avatar
Martin Reinecke committed
8
9
10
11

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

Martin Reinecke's avatar
Martin Reinecke committed
12

Martin Reinecke's avatar
Martin Reinecke committed
13
14
15
16
17
18
19
def random_alm(lmax, mmax):
    res = np.random.uniform(-1., 1., nalm(lmax, mmax)) \
     + 1j*np.random.uniform(-1., 1., nalm(lmax, mmax))
    # make a_lm with m==0 real-valued
    res[0:lmax+1].imag = 0.
    return res

Martin Reinecke's avatar
Martin Reinecke committed
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
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


def myalmdot(a1,a2,lmax,mmax,spin):
    return np.vdot(compress_alm(a1,lmax),compress_alm(a2,lmax))
def theta_extend(arr, spin):
    nth, nph = arr.shape
    arr2 = np.zeros(((nth-1)*2,nph))
    arr2[0:nth,:] = arr
    arr2[nth:,:] = np.roll(arr[nth-2:0:-1,:],nph//2,axis=1)
    if spin&1:
        arr2 = -arr2
    return arr2

def mydot(a1,a2,spin):
    return np.vdot(theta_extend(a1,spin),theta_extend(a2,spin))
Martin Reinecke's avatar
Martin Reinecke committed
41

Martin Reinecke's avatar
Martin Reinecke committed
42
def deltabeam(lmax,kmax):
Martin Reinecke's avatar
Martin Reinecke committed
43
    beam=np.zeros(nalm(lmax, kmax))+0j
Martin Reinecke's avatar
Martin Reinecke committed
44
45
46
47
    for l in range(lmax+1):
        beam[l] = np.sqrt((2*l+1.)/(4*np.pi))
    return beam

48
49
50
51
52
53
54
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
55

Martin Reinecke's avatar
Martin Reinecke committed
56
lmax=40
Martin Reinecke's avatar
Martin Reinecke committed
57
mmax=lmax
Martin Reinecke's avatar
Martin Reinecke committed
58
kmax=10
Martin Reinecke's avatar
Martin Reinecke committed
59
60
61
62


# get random sky a_lm
# the a_lm arrays follow the same conventions as those in healpy
Martin Reinecke's avatar
Martin Reinecke committed
63
slmT = random_alm(lmax, mmax)
Martin Reinecke's avatar
Martin Reinecke committed
64

Martin Reinecke's avatar
Martin Reinecke committed
65
# build beam a_lm
Martin Reinecke's avatar
Martin Reinecke committed
66
blmT = random_alm(lmax, kmax)
Martin Reinecke's avatar
Martin Reinecke committed
67
68

t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
69
# build interpolator object for slmT and blmT
Martin Reinecke's avatar
Martin Reinecke committed
70
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=1)
Martin Reinecke's avatar
Martin Reinecke committed
71
print("setup time: ",time.time()-t0)
72
nth = 2*lmax+1
Martin Reinecke's avatar
Martin Reinecke committed
73
nph = 2*mmax+1
Martin Reinecke's avatar
Martin Reinecke committed
74
#exit()
Martin Reinecke's avatar
Martin Reinecke committed
75
ptg = np.zeros((nth,nph,3))
76
ptg[:,:,0] = (np.pi*np.arange(nth)/(nth-1)).reshape((-1,1))
Martin Reinecke's avatar
Martin Reinecke committed
77
ptg[:,:,1] = (2*np.pi*np.arange(nph)/nph).reshape((1,-1))
Martin Reinecke's avatar
Martin Reinecke committed
78
ptg[:,:,2] = np.pi*0.2
Martin Reinecke's avatar
Martin Reinecke committed
79
t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
80
# do the actual interpolation
Martin Reinecke's avatar
Martin Reinecke committed
81
bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph))
Martin Reinecke's avatar
Martin Reinecke committed
82
print("interpolation time: ", time.time()-t0)
83
84
85
plt.subplot(2,2,1)
plt.imshow(bar.reshape((nth,nph)))
bar2 = np.zeros((nth,nph))
Martin Reinecke's avatar
Martin Reinecke committed
86
87
blmTfull = np.zeros(slmT.size)+0j
blmTfull[0:blmT.size] = blmT
88
for ith in range(nth):
Martin Reinecke's avatar
Martin Reinecke committed
89
    rbeamth=interpol_ng.rotate_alm(blmTfull, lmax, ptg[ith,0,2],ptg[ith,0,0],0)
90
    for iph in range(nph):
Martin Reinecke's avatar
Martin Reinecke committed
91
        rbeam=interpol_ng.rotate_alm(rbeamth, lmax, 0, 0, ptg[ith,iph,1])
Martin Reinecke's avatar
Martin Reinecke committed
92
        bar2[ith,iph] = convolve(slmT, rbeam, lmax).real
93
94
95
96
plt.subplot(2,2,2)
plt.imshow(bar2)
plt.subplot(2,2,3)
plt.imshow((bar2-bar.reshape((nth,nph))))
Martin Reinecke's avatar
Martin Reinecke committed
97
plt.show()
Martin Reinecke's avatar
Martin Reinecke committed
98
99
100
101
102

fake = np.random.uniform(-1.,1., bar.size)
foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
foo2.deinterpol(ptg.reshape((-1,3)), fake)
bla=foo2.getSlm(blmT)
Martin Reinecke's avatar
Martin Reinecke committed
103
print(myalmdot(slmT, bla, lmax, lmax, 0))
Martin Reinecke's avatar
Martin Reinecke committed
104
print(np.vdot(fake,bar))