demo.py 1.55 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
28
29
30
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 mydot(a1,a2,spin):
    return np.vdot(theta_extend(a1,spin),theta_extend(a2,spin))
Martin Reinecke's avatar
Martin Reinecke committed
31

32
33
lmax=30
kmax=13
Martin Reinecke's avatar
Martin Reinecke committed
34
35
36
37


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

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