TestEsti.py 5.19 KB
Newer Older
Matthieu Tristram's avatar
Matthieu Tristram committed
1
#!/usr/bin/env python
2 3 4
"""
Test script for xQML
"""
Syl's avatar
Syl committed
5

6 7
from __future__ import division

8

Matthieu Tristram's avatar
Matthieu Tristram committed
9 10
import numpy as np
import healpy as hp
11
from pylab import *
Matthieu Tristram's avatar
Matthieu Tristram committed
12 13
import timeit
import sys
Matthieu Tristram's avatar
Matthieu Tristram committed
14

Syl's avatar
Syl committed
15
import xqml
Matthieu Tristram's avatar
Matthieu Tristram committed
16
from xqml.xqml_utils import progress_bar, getstokes
17
from xqml.simulation import Karcmin2var
Matthieu Tristram's avatar
Matthieu Tristram committed
18
from xqml.simulation import extrapolpixwin
19 20
#ion()
#show()
Syl's avatar
Syl committed
21

Matthieu Tristram's avatar
Matthieu Tristram committed
22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
exp = "Big"
if len(sys.argv) > 1:
    if sys.argv[1].lower()[0] == "s":
        exp = "Small"

if exp == "Big":
    nside = 8
    dell = 1
    glat = 10
elif exp == "Small":
    nside = 64
    dell = 10
    glat = 80
else:
    print( "Need a patch !")

#lmax = nside
Syl's avatar
Syl committed
39
lmax = 2 * nside - 1
Matthieu Tristram's avatar
Matthieu Tristram committed
40 41
nsimu = 100
MODELFILE = 'planck_base_planck_2015_TTlowP.fits'
42
Slmax = 3*nside-1
Syl's avatar
Syl committed
43

Syl's avatar
Syl committed
44

Syl's avatar
Syl committed
45
# provide list of specs to be computed, and/or options
46 47
spec = ['EE','BB','EB']
#spec = ['TT','EE','BB','TE']#'EE','BB', 'EB']#, 'TE', 'TB']
48
pixwin = True
Syl's avatar
Syl committed
49
ellbins = np.arange(2, lmax + 2, dell)
Matthieu Tristram's avatar
Matthieu Tristram committed
50 51 52
ellbins[-1] = lmax+1

muKarcmin = 1.0
53
fwhm = 10
Matthieu Tristram's avatar
Matthieu Tristram committed
54 55 56 57 58 59



##############################
#input model
clth = np.array(hp.read_cl(MODELFILE))
60
clth = array( list(clth) + list(clth[0:2]*0.))
Matthieu Tristram's avatar
Matthieu Tristram committed
61 62 63
lth = arange(2, lmax+1)
##############################

Syl's avatar
Syl committed
64 65


Matthieu Tristram's avatar
Matthieu Tristram committed
66
##############################
Syl's avatar
Syl committed
67
# Create mask
Syl's avatar
Syl committed
68

Matthieu Tristram's avatar
Matthieu Tristram committed
69 70
t, p = hp.pix2ang(nside, range(hp.nside2npix(nside)))
mask = np.ones(hp.nside2npix(nside), bool)
Syl's avatar
Syl committed
71 72
# import random
# random.shuffle(mask)
Syl's avatar
Syl committed
73

Matthieu Tristram's avatar
Matthieu Tristram committed
74 75 76 77
if exp == "Big":
    mask[abs(90 - rad2deg(t)) < glat] = False
elif exp == "Small":
    mask[(90 - rad2deg(t)) < glat] = False
Syl's avatar
Syl committed
78 79

fsky = np.mean(mask)
Syl's avatar
Syl committed
80
npix = sum(mask)
81
print("%s patch: fsky=%.2g %% (npix=%d)" % (exp,100*fsky,npix))
Matthieu Tristram's avatar
Matthieu Tristram committed
82 83 84
toGB = 1024. * 1024. * 1024.
emem = 8.*(npix*2*npix*2) * ( len(lth)*2 ) / toGB
print("mem=%.2g Gb" % emem)
Matthieu Tristram's avatar
Matthieu Tristram committed
85 86
##############################

Syl's avatar
Syl committed
87 88


Matthieu Tristram's avatar
Matthieu Tristram committed
89
stokes, spec, istokes, ispecs = getstokes( spec=spec)
Syl's avatar
Syl committed
90 91 92 93 94
print(stokes, spec, istokes, ispecs)
nspec = len(spec)
nstoke = len(stokes)


95 96
# ############## Generate White Noise ###############
pixvar = Karcmin2var(muKarcmin*1e-6, nside)
Syl's avatar
Syl committed
97 98 99
varmap = ones((nstoke * npix)) * pixvar
NoiseVar = np.diag(varmap)

Matthieu Tristram's avatar
Matthieu Tristram committed
100 101 102 103


# ############## Initialise xqml class ###############
start = timeit.default_timer()
Matthieu Tristram's avatar
Matthieu Tristram committed
104
esti = xqml.xQML(mask, ellbins, clth, NA=NoiseVar, NB=NoiseVar, lmax=lmax, fwhm=fwhm, spec=spec)
Matthieu Tristram's avatar
Matthieu Tristram committed
105 106 107 108 109
s1 = timeit.default_timer()
print( "Init: %d sec" % (s1-start))
ellval = esti.lbin()

# ############## Compute Analytical variance ###############
110 111
#V  = esti.get_covariance(cross=True )
#Va = esti.get_covariance(cross=False)
Matthieu Tristram's avatar
Matthieu Tristram committed
112 113
s2 = timeit.default_timer()
print( "construct covariance: %d sec" % (s2-s1))
Syl's avatar
Syl committed
114

115

Matthieu Tristram's avatar
Matthieu Tristram committed
116
# ############## Construct MC ###############
117
allcla = []
Syl's avatar
Syl committed
118
allcl = []
Matthieu Tristram's avatar
Matthieu Tristram committed
119 120 121
t = []
bl = hp.gauss_beam(deg2rad(fwhm), lmax=Slmax)
fpixw = extrapolpixwin(nside, Slmax, pixwin=pixwin)
Syl's avatar
Syl committed
122
for n in np.arange(nsimu):
Matthieu Tristram's avatar
Matthieu Tristram committed
123
    progress_bar(n, nsimu)
Syl's avatar
Syl committed
124
    cmb = np.array(hp.synfast(clth[:, :len(fpixw)]*(fpixw*bl)**2, nside,
Matthieu Tristram's avatar
Matthieu Tristram committed
125
                   pixwin=False, lmax=Slmax, fwhm=0.0, new=True, verbose=False))
Syl's avatar
Syl committed
126 127 128
    cmbm = cmb[istokes][:, mask]
    dmA = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
    dmB = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
Matthieu Tristram's avatar
Matthieu Tristram committed
129
    s1 = timeit.default_timer()
Syl's avatar
Syl committed
130
    allcl.append(esti.get_spectra(dmA, dmB))
Matthieu Tristram's avatar
Matthieu Tristram committed
131
    t.append( timeit.default_timer() - s1)
132
    allcla.append(esti.get_spectra(dmA))
Syl's avatar
Syl committed
133

134
print( "get_spectra: %f sec" % mean(t))
Syl's avatar
Syl committed
135 136
hcl = mean(allcl, 0)
scl = std(allcl, 0)
137 138 139
hcla = mean(allcla, 0)
scla = std(allcla, 0)

Matthieu Tristram's avatar
Matthieu Tristram committed
140 141 142 143


# ############## Plot results ###############

144 145
figure(figsize=[12, 8])
clf()
Syl's avatar
Syl committed
146
Delta = (ellbins[1:] - ellbins[:-1])/2.
Matthieu Tristram's avatar
Matthieu Tristram committed
147

148 149
subplot(3, 2, 1)
title("Cross")
Syl's avatar
Syl committed
150
plot(lth, (lth*(lth+1)/2./np.pi)[:, None]*clth[ispecs][:, lth].T, '--k')
Matthieu Tristram's avatar
Matthieu Tristram committed
151 152
for s in np.arange(nspec):
    errorbar(ellval, ellval*(ellval+1)/2./np.pi*hcl[s], yerr=scl[s], xerr=Delta, fmt='o', color='C%i' % ispecs[s], label=r"$%s$" % spec[s])
153
semilogy()
Syl's avatar
Syl committed
154
ylabel(r"$D_\ell$")
155 156 157 158
legend(loc=4, frameon=True)

subplot(3, 2, 2)
title("Auto")
Syl's avatar
Syl committed
159
plot(lth,(lth*(lth+1)/2./np.pi)[:, None]*clth[ispecs][:, lth].T, '--k')
Matthieu Tristram's avatar
Matthieu Tristram committed
160 161
for s in np.arange(nspec):
    errorbar(ellval, ellval*(ellval+1)/2./np.pi*hcla[s], yerr=scla[s], xerr=Delta, fmt='o', color='C%i' % ispecs[s], label=r"$%s$" % spec[s])
Syl's avatar
Syl committed
162 163
semilogy()

164
subplot(3, 2, 3)
Matthieu Tristram's avatar
Matthieu Tristram committed
165
for s in np.arange(nspec):
166 167
    plot(ellval, scl[s], 'o', color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s])
#    plot(ellval, sqrt(diag(V)).reshape(nspec, -1)[s], '-', color='C%i' % ispecs[s])
Syl's avatar
Syl committed
168 169
ylabel(r"$\sigma(C_\ell)$")
semilogy()
170 171

subplot(3, 2, 4)
Matthieu Tristram's avatar
Matthieu Tristram committed
172
for s in np.arange(nspec):
173 174
    plot(ellval, scla[s], 'o', color='C%i' % ispecs[s], label=r"$\sigma^{%s}_{\rm MC}$" % spec[s])
#    plot(ellval, sqrt(diag(Va)).reshape(nspec, -1)[s], '-', color='C%i' % ispecs[s])
175
semilogy()
Syl's avatar
Syl committed
176

177
subplot(3, 2, 5)
Matthieu Tristram's avatar
Matthieu Tristram committed
178 179
for s in np.arange(nspec):
    plot(ellval, (hcl[s]-esti.BinSpectra(clth)[s])/(scl[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s])
Syl's avatar
Syl committed
180 181 182 183
ylabel(r"$R[C_\ell]$")
xlabel(r"$\ell$")
ylim(-3, 3)
grid()
184 185

subplot(3, 2, 6)
Matthieu Tristram's avatar
Matthieu Tristram committed
186 187
for s in np.arange(nspec):
    plot(ellval, (hcla[s]-esti.BinSpectra(clth)[s])/(scla[s]/sqrt(nsimu)), '--o', color='C%i' % ispecs[s])
188 189 190 191
xlabel(r"$\ell$")
ylim(-3, 3)
grid()

Syl's avatar
Syl committed
192 193
show()

Syl's avatar
Syl committed
194

195 196 197
if __name__ == "__main__":
    """
    Run the doctest using
Syl's avatar
Syl committed
198

199
    python simulation.py
Syl's avatar
Syl committed
200

201 202 203 204 205 206 207
    If the tests are OK, the script should exit gracefuly, otherwise the
    failure(s) will be printed out.
    """
    import doctest
    if np.__version__ >= "1.14.0":
        np.set_printoptions(legacy="1.13")
    doctest.testmod()