demo.py 3.35 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=60
kmax=13
44
ncomp=1
Martin Reinecke's avatar
Martin Reinecke committed
45
separate=False
Martin Reinecke's avatar
Martin Reinecke committed
46
nptg = 1000000
47
epsilon = 1e-4
Martin Reinecke's avatar
Martin Reinecke committed
48
49
50
51
ofactor = 1.5
nthreads = 0  # use as many threads as available

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

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

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


t0=time.time()
62
# build interpolator object for slm and blm
63
foo = pyinterpol_ng.PyInterpolator(slm,blm,separate,lmax, kmax, epsilon=epsilon, ofactor=ofactor, nthreads=nthreads)
64
print("setup time: ",time.time()-t0)
65
print("support:",foo.support())
66
67
68
69
70
71
72
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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()
Martin Reinecke's avatar
Martin Reinecke committed
96
97
98


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