totalconvolve_usage.py 2.12 KB
Newer Older
Martin Reinecke's avatar
comment  
Martin Reinecke committed
1
2
# Short usage demo for the interpol_ng module

Martin Reinecke's avatar
Martin Reinecke committed
3
import ducc0.totalconvolve as totalconvolve
Martin Reinecke's avatar
Martin Reinecke committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import numpy as np


# establish a random number generator
rng = np.random.default_rng(np.random.SeedSequence(42))


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


def random_alm(lmax, mmax, ncomp):
    res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
     + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
    # make a_lm with m==0 real-valued
    res[0:lmax+1,:].imag = 0.
    return res


# simulation parameters
lmax=1024 # highest l and m moment for sky, highest l moment for beam
kmax=13 # highest m moment for beam
ncomp=3 # T, E and B

# get random sky a_lm
# the a_lm arrays follow the same conventions as those in healpy
slm = random_alm(lmax, lmax, ncomp)

# build beam a_lm
blm = random_alm(lmax, kmax, ncomp)

# build random pointings
npnt = 1000000 # play with this to measure performance
ptg = rng.uniform(0., 1., (npnt,3))
ptg[:,0]*=np.pi   # theta
ptg[:,1]*=2*np.pi # phi
ptg[:,2]*=2*np.pi # psi

# build the interpolator
# For a "classic" CMB experiment we want to combine the T, E and B signals,
# so we set `separate` to `False`

print("classic interpolator setup...")
Martin Reinecke's avatar
Martin Reinecke committed
47
inter_classic = totalconvolve.PyInterpolator(
Martin Reinecke's avatar
Martin Reinecke committed
48
49
50
51
52
53
54
55
56
57
58
59
    slm,blm,separate=False,lmax=lmax, kmax=kmax, epsilon=1e-4, nthreads=2)
print("...done")

# get interpolated values
print("interpolating...")
res = inter_classic.interpol(ptg)
print("...done")

# res is an array of shape (nptg, 1).
# If we had specified `separate=True`, it would be of shape(nptg, 3).
print(res.shape)

Martin Reinecke's avatar
Martin Reinecke committed
60
61
# Since the interpolator object holds large data structures, it should be
# deleted once it is no longer needed
Martin Reinecke's avatar
Martin Reinecke committed
62
63
64
65
66
67
68
del inter_classic

# Now the same thing for an experiment with HWP. In this case we need the
# interpolated T, E and B signals separate.
separate = True

print("HWP interpolator setup...")
Martin Reinecke's avatar
Martin Reinecke committed
69
inter_hwp = totalconvolve.PyInterpolator(
Martin Reinecke's avatar
Martin Reinecke committed
70
71
72
73
74
75
76
77
78
79
    slm,blm,separate=True,lmax=lmax, kmax=kmax, epsilon=1e-4, nthreads=2)
print("...done")

# get interpolated values
print("interpolating...")
res = inter_hwp.interpol(ptg)
print("...done")

# now res has shape(nptg,3)
print(res.shape)