demo.py 3.43 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 pysharp
import time
Martin Reinecke's avatar
Martin Reinecke committed
5
import matplotlib.pyplot as plt
Martin Reinecke's avatar
Martin Reinecke committed
6

7
np.random.seed(48)
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

13
14
15
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
16
    # make a_lm with m==0 real-valued
17
    res[0:lmax+1,:].imag = 0.
Martin Reinecke's avatar
Martin Reinecke committed
18
19
    return res

20

Martin Reinecke's avatar
Martin Reinecke committed
21
22
23
24
25
26
27
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

28

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

32
33
34
35
36
37
38
39
40
41

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
42
43
lmax=2048
kmax=8
44
ncomp=1
Martin Reinecke's avatar
Martin Reinecke committed
45
46
separate=False
nptg = 10000000
47
ncomp2 = ncomp if separate else 1
48
epsilon = 1e-4
Martin Reinecke's avatar
Martin Reinecke committed
49
50
51
52
ofactor = 1.5
nthreads = 0  # use as many threads as available

ncomp2 = ncomp if separate else 1
Martin Reinecke's avatar
Martin Reinecke committed
53
54
55

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

Martin Reinecke's avatar
Martin Reinecke committed
58
# build beam a_lm
59
blm = random_alm(lmax, kmax, ncomp)
60
61
62


t0=time.time()
63
# build interpolator object for slm and blm
64
foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
65
print("setup time: ",time.time()-t0)
66
print("support:",foo.support())
67
68
69
70
71
72
73
nth = lmax+1
nph = 2*lmax+1


# compute a convolved map at a fixed psi and compare it to a map convolved
# "by hand"

Martin Reinecke's avatar
Martin Reinecke committed
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
# ptg = np.zeros((nth,nph,3))
# ptg[:,:,0] = (np.pi*(0.5+np.arange(nth))/nth).reshape((-1,1))
# ptg[:,:,1] = (2*np.pi*(0.5+np.arange(nph))/nph).reshape((1,-1))
# ptg[:,:,2] = np.pi*0.2
# t0=time.time()
# # do the actual interpolation
# bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph,ncomp2))
# print("interpolation time: ", time.time()-t0)
# plt.subplot(2,2,1)
# plt.imshow(bar[:,:,0])
# bar2 = np.zeros((nth,nph))
# blmfull = np.zeros(slm.shape)+0j
# blmfull[0:blm.shape[0],:] = blm
# for ith in range(nth):
#     rbeamth=pyinterpol_ng.rotate_alm(blmfull[:,0], lmax, ptg[ith,0,2],ptg[ith,0,0],0)
#     for iph in range(nph):
#         rbeam=pyinterpol_ng.rotate_alm(rbeamth, lmax, 0, 0, ptg[ith,iph,1])
#         bar2[ith,iph] = convolve(slm[:,0], rbeam, lmax).real
# plt.subplot(2,2,2)
# plt.imshow(bar2)
# plt.subplot(2,2,3)
# plt.imshow(bar2-bar[:,:,0])
# plt.show()


ptg=np.random.uniform(0.,1.,3*nptg).reshape(nptg,3)
100
101
102
ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi
103
104
#foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=1e-6, nthreads=2)
t0=time.time()
105
bar=foo.interpol(ptg)
Martin Reinecke's avatar
Martin Reinecke committed
106
del foo
107
108
print("interpolation time: ", time.time()-t0)
fake = np.random.uniform(0.,1., (ptg.shape[0],ncomp2))
109
foo2 = pyinterpol_ng.PyInterpolator(lmax, kmax, ncomp2, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
110
111
112
113
114
115
116
117
118
t0=time.time()
foo2.deinterpol(ptg.reshape((-1,3)), fake)
print("deinterpolation time: ", time.time()-t0)
t0=time.time()
bla=foo2.getSlm(blm)
print("getSlm time: ", 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)])
print(v1/v2-1.)