demo.py 1.44 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

7
np.random.seed(2099)
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
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):
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
lmax=30
kmax=13
Martin Reinecke's avatar
Martin Reinecke committed
32
33
34
35


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

Martin Reinecke's avatar
Martin Reinecke committed
38
# build beam a_lm
Martin Reinecke's avatar
Martin Reinecke committed
39
blmT = random_alm(lmax, kmax)
40
41
42
43
ptg=np.random.uniform(0.,1.,3*1000000).reshape(1000000,3)
ptg[:,0]*=np.pi
ptg[:,1]*=2*np.pi
ptg[:,2]*=2*np.pi
Martin Reinecke's avatar
Martin Reinecke committed
44
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=1)
45
46
bar=foo.interpol(ptg)
fake = np.random.uniform(0.,1., ptg.shape[0])
Martin Reinecke's avatar
Martin Reinecke committed
47
48
49
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
cleanup    
Martin Reinecke committed
50
print(myalmdot(slmT, bla, lmax, lmax, 0))
Martin Reinecke's avatar
Martin Reinecke committed
51
print(np.vdot(fake,bar))
Martin Reinecke's avatar
cleanup    
Martin Reinecke committed
52
print(myalmdot(slmT, bla, lmax, lmax, 0)/np.vdot(fake,bar))