Commit 58be8805 authored by Syl's avatar Syl
Browse files

few corrections

parent f1d0cc3b
......@@ -3,6 +3,7 @@ Test script for xQML
Author: Vanneste
"""
from __future__ import division
import timeit
......@@ -11,107 +12,135 @@ import healpy as hp
from pylab import *
import astropy.io.fits as fits
from xqml import xQML
from libcov import compute_ds_dcb
import xqml
import libcov
import simulation
from xqml_utils import progress_bar
from simulation import getstokes, muKarcmin2var, GetBinningMatrix
from simulation import muKarcmin2var, GetBinningMatrix
from simulation import extrapolpixwin
if __name__ == "__main__":
nside = 4
lmax = 3 * nside - 1
Slmax = 3 * nside - 1
deltal = 1
nsimu = 500
clth = np.array(hp.read_cl('planck_base_planck_2015_TTlowP.fits'))
lth = arange(2, lmax+1)
ellbins = arange(2, lmax + 2, deltal)
ellbins[-1] = lmax + 1
P, Q, ell, ellval = GetBinningMatrix(ellbins, lmax)
nbins = len(ellbins) - 1
# Create mask
t, p = hp.pix2ang(nside, range(hp.nside2npix(nside)))
mask = np.ones(hp.nside2npix(nside), bool)
mask[abs(90 - rad2deg(t)) < 10] = False
npix = sum(mask)
fwhm = 0.5
bl = hp.gauss_beam(deg2rad(fwhm), lmax=Slmax + 1)
allStoke, der, ind = getstokes(polar=True, temp=False, EBTB=False)
nder = len(der)
muKarcmin = 0.1
pixvar = muKarcmin2var(muKarcmin, nside)
varmap = ones((2 * npix)) * pixvar
NoiseVar = np.diag(varmap)
cmb = np.array(hp.synfast(
clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False))
noise = (randn(len(varmap)) * varmap**0.5).reshape(2, -1)
# dm = cmb[1:, mask] + noise
dm = cmb[1:][:, mask] + noise
# ############## Compute ds_dcb ###############
ip = arange(hp.nside2npix(nside))
ipok = ip[mask]
Pl, S = compute_ds_dcb(
ellbins, nside, ipok, bl, clth, Slmax,
polar=True, temp=False, EBTB=False,
pixwining=True, timing=True, MC=False)
# Pl = Pl.reshape((nder)*(np.shape(Pl)[1]), 2 * npix, 2 * npix)
# ############## Compute spectra ###############
esti = xQML(mask, ellbins, clth, Pl=Pl, S=S, fwhm=fwhm)
esti.construct_esti(NoiseVar, NoiseVar)
cl = esti.get_spectra(dm, dm)
V = esti.get_covariance()
# ############## Construct MC ###############
allcl = []
esti = xQML(mask, ellbins, clth, Pl=Pl, fwhm=fwhm)
esti.construct_esti(NoiseVar, NoiseVar)
fpixw = extrapolpixwin(nside, lmax+2, pixwining=True)
start = timeit.default_timer()
for n in np.arange(nsimu):
progress_bar(n, nsimu, timeit.default_timer() - start)
cmb = np.array(hp.synfast(clth[:, :len(fpixw)]*fpixw**2, nside,
fwhm=deg2rad(fwhm), new=True, verbose=False))
dm = cmb[1:, mask] + (randn(2 * npix) * sqrt(varmap)).reshape(2, npix)
allcl.append(esti.get_spectra(dm, dm))
figure()
subplot(3, 1, 1)
plot(lth, clth.transpose()[lth, 1: 3], '--k')
hcl = mean(allcl, 0).transpose()
scl = std(allcl, 0).transpose()
plot(ellval, hcl, 'b.')
plot(ellval, hcl + scl, 'r--', label=r"$\pm 1\sigma$")
plot(ellval, hcl - scl, 'r--')
ylabel(r"$C_\ell$")
semilogy()
legend(loc=4)
subplot(3, 1, 2)
cosmic = sqrt(2./(2 * lth + 1)) / mean(mask) * clth[1: 3, lth]
plot(lth, cosmic.transpose(), '--k')
plot(ellval, scl, 'r-', label=r"$\sigma_{\rm MC}$")
plot(ellval, sqrt(diag(V)).reshape(nder, -1).transpose(), 'b.')
ylabel(r"$\sigma(C_\ell)$")
semilogy()
legend(loc=4)
subplot(3, 1, 3)
plot(ellval, (hcl-clth[1: 3, lth].T)/(scl/sqrt(nsimu)), 'o')
ylabel(r"$R[C_\ell]$")
xlabel(r"$\ell$")
ylim(-3, 3)
show()
ion()
# if __name__ == "__main__":
nside = 4
lmax = 2 * nside - 1
Slmax = 2 * nside - 1
deltal = 1
nsimu = 10000
clth = np.array(hp.read_cl('planck_base_planck_2015_TTlowP.fits'))
Clthshape = zeros(((6,)+shape(clth)[1:]))
Clthshape[:4] = clth
clth = Clthshape
EB = 0.5
clth[4] = EB*sqrt(clth[1]*clth[2])
TB = 0.5
clth[5] = TB*sqrt(clth[0]*clth[2])
lth = arange(2, lmax+1)
spec = ['EB', 'TE', 'TB']
temp = True
polar = True
corr = False
pixwin = False
ellbins = arange(2, lmax + 2, deltal)
ellbins[-1] = lmax + 1
P, Q, ell, ellval = GetBinningMatrix(ellbins, lmax)
nbins = len(ellbins) - 1
# Create mask
t, p = hp.pix2ang(nside, range(hp.nside2npix(nside)))
mask = np.ones(hp.nside2npix(nside), bool)
mask[abs(90 - rad2deg(t)) < 60] = False
npix = sum(mask)
fwhm = 0.5
bl = hp.gauss_beam(deg2rad(fwhm), lmax=Slmax + 1)
stokes, spec, istokes, ispecs = simulation.getstokes(
spec=spec, polar=polar, temp=temp, corr=corr)
print(stokes, spec, istokes, ispecs)
nspec = len(spec)
nstoke = len(stokes)
# ############## Compute ds_dcb ###############
ip = arange(hp.nside2npix(nside))
ipok = ip[mask]
Pl, S = libcov.compute_ds_dcb(ellbins, nside, ipok, bl, clth, Slmax, spec=spec,
pixwin=pixwin, timing=True, MC=False)
# ############## Compute spectra ###############
muKarcmin = 1.0
pixvar = muKarcmin2var(muKarcmin, nside)
varmap = ones((nstoke * npix)) * pixvar
NoiseVar = np.diag(varmap)
cmb = np.array(hp.synfast(clth, nside, fwhm=deg2rad(fwhm), pixwin=pixwin,
new=True, verbose=False, lmax=Slmax))
noise = (randn(len(varmap)) * varmap**0.5).reshape(nstoke, -1)
dm = cmb[istokes][:, mask] + noise
esti = xqml.xQML(mask, ellbins, clth, Pl=Pl, S=S, fwhm=fwhm,
spec=spec, temp=temp, polar=polar, corr=corr)
esti.construct_esti(NoiseVar, NoiseVar)
cl = esti.get_spectra(dm, dm)
V = esti.get_covariance()
# ############## Construct MC ###############
allcl = []
# allcmb = []
# esti = xqml.xQML(mask, ellbins, clth, Pl=Pl, fwhm=fwhm, spec=spec, temp=temp,
# polar=polar, corr=corr)
esti.construct_esti(NoiseVar, NoiseVar)
fpixw = extrapolpixwin(nside, lmax+2, pixwin=pixwin)
start = timeit.default_timer()
for n in np.arange(nsimu):
progress_bar(n, nsimu, timeit.default_timer() - start)
cmb = np.array(hp.synfast(clth[:, :len(fpixw)]*(fpixw*bl)**2, nside,
lmax=Slmax, fwhm=deg2rad(fwhm), new=True, verbose=False))
cmbm = cmb[istokes][:, mask]
dmA = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
dmB = cmbm + (randn(nstoke * npix) * sqrt(varmap)).reshape(nstoke, npix)
# allcmb.append(cmbm)
allcl.append(esti.get_spectra(dmA, dmB))
figure(figsize=[10, 8])
clf()
subplot(3, 1, 1)
plot(lth, clth[ispecs][:, lth].T, '--k')
hcl = mean(allcl, 0)
scl = std(allcl, 0)
[plot(ellval, hcl[s], 'o', color='C%i' % s, label=r"$%s$" % spec[s])
for s in np.arange(nspec)]
[fill_between(ellval, (hcl - scl/sqrt(nsimu))[s], (hcl + scl/sqrt(nsimu))[s],
color='C%i' % s, alpha=0.2) for s in np.arange(nspec)]
ylabel(r"$C_\ell$")
semilogy()
legend(loc=4)
subplot(3, 1, 2)
cosmic = sqrt(2./(2 * lth + 1)) / mean(mask) * clth[ispecs][:, lth]
# plot(lth, cosmic.transpose(), '--k')
[plot(ellval, scl[s], '--', color='C%i' % s, label=r"$\sigma^{%s}_{\rm MC}$" %
spec[s]) for s in np.arange(nspec)]
[plot(ellval, sqrt(diag(V)).reshape(nspec, -1)[s], 'o', color='C%i' % s)
for s in np.arange(nspec)]
ylabel(r"$\sigma(C_\ell)$")
semilogy()
legend(loc=4)
subplot(3, 1, 3)
[plot(ellval, (hcl[s]-clth[ispecs[s]][lth])/(scl[s]/sqrt(nsimu)), '--o',
color='C%i' % s) for s in np.arange(nspec)]
ylabel(r"$R[C_\ell]$")
xlabel(r"$\ell$")
ylim(-3, 3)
grid()
show()
# savefig("../Plots/Git/"+"test0.png")
if __name__ == "__main__":
"""
......
......@@ -42,7 +42,7 @@ def Pl(ds_dcb):
return np.copy(ds_dcb).reshape(2 * (np.shape(ds_dcb)[1]), nnpix, nnpix)
def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, EBTB=False):
def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, corr=False):
"""
Compute correlation matrix S = sum_l Pl*Cl
......@@ -58,7 +58,7 @@ def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, EBTB=False):
If True, get Stokes parameters for polar (default: True)
temp : bool
If True, get Stokes parameters for temperature (default: False)
EBTB : bool
corr : bool
If True, get Stokes parameters for EB and TB (default: False)
Returns
......@@ -76,7 +76,7 @@ def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, EBTB=False):
[[ 0 280 560 840 1120]
[1400 1680 1960 2240 2520]]
"""
if EBTB:
if corr:
xx = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
ind = [0, 1, 2, 3, 4, 5]
else:
......@@ -85,7 +85,7 @@ def CorrelationMatrix(Clth, Pl, ellbins, polar=True, temp=False, EBTB=False):
if not temp:
allStoke = ['Q', 'U']
if EBTB:
if corr:
xx = ['EE', 'BB', 'EB']
ind = [1, 2, 5]
else:
......
"""
Set of routines to ...
Author: Vanneste
"""
from __future__ import division
import numpy as np
def polrotangle(ri, rj):
"""
Computes cosine and sine of twice the angle between pixels i and j.
Parameters
----------
ri : 3D array of floats
Coordinates of vector corresponding to input pixels i following
healpy.pix2vec(nside,ipix) output
rj : 3D array of floats
Coordinates of vector corresponding to input pixels j following
healpy.pix2vec(nside,jpix) output
Returns
----------
cos2a : 1D array of floats
Cosine of twice the angle between pixels i and j
sin2a : 1D array of floats
Sine of twice the angle between pixels i and j
Example
----------
>>> cos2a, sin2a = polrotangle([0.1,0.2,0.3], [0.4,0.5,0.6])
>>> print(round(cos2a,5),round(sin2a,5))
(0.06667, 0.37333)
"""
z = np.array([0.0, 0.0, 1.0])
# Compute ri^rj : unit vector for the great circle connecting i and j
rij = np.cross(ri, rj)
norm = np.sqrt(np.dot(rij, np.transpose(rij)))
# case where pixels are identical or diametrically opposed on the sky
if norm <= 1e-15:
cos2a = 1.0
sin2a = 0.0
return cos2a, sin2a
rij = rij / norm
# Compute z^ri : unit vector for the meridian passing through pixel i
ris = np.cross(z, ri)
norm = np.sqrt(np.dot(ris, np.transpose(ris)))
# case where pixels is at the pole
if norm <= 1e-15:
cos2a = 1.0
sin2a = 0.0
return cos2a, sin2a
ris = ris / norm
# Now, the angle we want is that
# between these two great circles: defined by
cosa = np.dot(rij, np.transpose(ris))
# the sign is more subtle : see tegmark et de oliveira costa 2000 eq. A6
rijris = np.cross(rij, ris)
sina = np.dot(rijris, np.transpose(ri))
# so now we have directly cos2a and sin2a
cos2a = 2.0 * cosa * cosa - 1.0
sin2a = 2.0 * cosa * sina
return cos2a, sin2a
if __name__ == "__main__":
"""
Run the doctest using
python simulation.py
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()
This diff is collapsed.
"""
Set of routines to generate basic simulations
Set of routines to generate basic simulations
Author: Vanneste
"""
......@@ -96,62 +96,81 @@ def getNl(pixvar, nside, nbins):
return Nl
def getstokes(polar=True, temp=False, EBTB=False):
def getstokes(spec=None, temp=False, polar=False, corr=False):
"""
Get the Stokes parameters number and name(s)
Parameters
----------
spec : bool
If True, get Stokes parameters for polar (default: True)
polar : bool
If True, get Stokes parameters for polar (default: True)
temp : bool
If True, get Stokes parameters for temperature (default: False)
EBTB : bool
corr : bool
If True, get Stokes parameters for EB and TB (default: False)
Returns
----------
allStoke : list of string
stokes : list of string
Stokes variables names
spec : int
Spectra names
ind : list
istokes : list
Indexes of power spectra
Example
----------
>>> getstokes(polar=True, temp=False, EBTB=False)
>>> getstokes(polar=True, temp=False, corr=False)
(['Q', 'U'], ['EE', 'BB'], [1, 2])
>>> getstokes(polar=True, temp=True, EBTB=False)
>>> getstokes(polar=True, temp=True, corr=False)
(['I', 'Q', 'U'], ['TT', 'EE', 'BB', 'TE'], [0, 1, 2, 3])
>>> getstokes(polar=True, temp=True, EBTB=True)
>>> getstokes(polar=True, temp=True, corr=True)
(['I', 'Q', 'U'], ['TT', 'EE', 'BB', 'TE', 'EB', 'TB'], [0, 1, 2, 3, 4, 5])
"""
allStoke = ['I', 'Q', 'U']
if EBTB:
spec = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
ind = [0, 1, 2, 3, 4, 5]
if spec is not None:
_temp = "TT" in spec or "TE" in spec or "TB" in spec or temp
_polar = "EE" in spec or "BB" in spec or "TE" in spec or "TB" in \
spec or "EB" in spec or polar
_corr = "TE" in spec or "TB" in spec or "EB" in spec or corr
if not _temp and not _polar and not _corr:
print("invalid spectra list")
else:
spec = ['TT', 'EE', 'BB', 'TE']
ind = [0, 1, 2, 3]
if not temp:
allStoke = ['Q', 'U']
if EBTB:
spec = ['EE', 'BB', 'EB']
ind = [1, 2, 4]
else:
spec = ['EE', 'BB']
ind = [1, 2]
if not polar:
allStoke = ['I']
spec = ['TT']
ind = [0]
return allStoke, spec, ind
_temp = temp
_polar = polar
_corr = corr
speclist = []
if _temp or (spec is None and corr):
speclist.extend(["TT"])
if _polar:
speclist.extend(["EE", "BB"])
if spec is not None and not corr:
if 'TE' in spec:
speclist.extend(["TE"])
if 'EB' in spec:
speclist.extend(["EB"])
if 'TB' in spec:
speclist.extend(["TB"])
elif _corr:
speclist.extend(["TE", "EB", "TB"])
stokes = []
if _temp:
stokes.extend(["I"])
if _polar:
stokes.extend(["Q", "U"])
ispecs = [['TT', 'EE', 'BB', 'TE', 'EB', 'TB'].index(s) for s in speclist]
istokes = [['I', 'Q', 'U'].index(s) for s in stokes]
return stokes, speclist, istokes, ispecs
def GetBinningMatrix(
ellbins, lmax, norm=False, polar=True,
temp=False, EBTB=False):
temp=False, corr=False):
"""
Return P (m,n) and Q (n,m) binning matrices such that
Cb = P.Cl and Vbb = P.Vll.Q with m the number of bins and
......@@ -171,7 +190,7 @@ def GetBinningMatrix(
If True, get Stokes parameters for polar (default: True)
temp : bool
If True, get Stokes parameters for temperature (default: False)
EBTB : bool
corr : bool
If True, get Stokes parameters for EB and TB (default: False)
Returns
......@@ -229,8 +248,7 @@ def GetBinningMatrix(
[ 3. 7.]
"""
# ### define Stokes
allStoke, der, ind = getstokes(polar, temp, EBTB)
nder = len(der)
nspec = 1
nbins = len(ellbins) - 1
ellmin = np.array(ellbins[0: nbins])
......@@ -248,10 +266,10 @@ def GetBinningMatrix(
for i in np.arange(nbins):
masklm.append(((ell[:-1] >= minell[i]) & (ell[:-1] <= maxell[i])))
allmasklm = nder*[list(masklm)]
allmasklm = nspec*[list(masklm)]
masklM = np.array(sparse.block_diag(allmasklm[:]).toarray())
binsnorm = np.array(
nder * [list(np.arange(minell[0], np.max(ellbins)))]).flatten()
nspec * [list(np.arange(minell[0], np.max(ellbins)))]).flatten()
binsnorm = binsnorm*(binsnorm+1)/2./np.pi
P = np.array(masklM)*1.
......@@ -263,7 +281,7 @@ def GetBinningMatrix(
return P, Q, ell, ellval
def extrapolpixwin(nside, Slmax, pixwining=True):
def extrapolpixwin(nside, Slmax, pixwin=True):
'''
Parameters
----------
......@@ -271,7 +289,7 @@ def extrapolpixwin(nside, Slmax, pixwining=True):
Healpix map resolution
Slmax : int
Maximum multipole value computed for the pixel covariance pixel matrix
pixwining : bool
pixwin : bool
If True, multiplies the beam window function by the pixel
window function. Default: True
......@@ -293,7 +311,7 @@ def extrapolpixwin(nside, Slmax, pixwining=True):
[ 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
1. 1.]
'''
if pixwining:
if pixwin:
prepixwin = np.array(hp.pixwin(nside))
poly = np.polyfit(np.arange(len(prepixwin)), np.log(prepixwin),
deg=3, w=np.sqrt(prepixwin))
......@@ -318,4 +336,4 @@ if __name__ == "__main__":
import doctest
if np.__version__ >= "1.14.0":
np.set_printoptions(legacy="1.13")
doctest.testmod()
\ No newline at end of file
doctest.testmod()
"""
Set of routines to compute the pixel covariance matrix using
Legendre polynomials
Author: Vanneste
"""
from __future__ import division
import math
import numpy as np
from scipy import special
def dlss(z, s1, s2, lmax):
"""
Computes the reduced Wigner D-function d^l_ss'
Parameters
----------
z : float
Cosine of the angle between two pixels
s1 : int
Spin number 1
s2 : int
Spin number 2
lmax : int
Maximum multipole
Returns
----------
d : 1D array of floats
???
Example
----------
>>> d = dlss(0.1, 2, 2, 5)
>>> print(round(sum(d),5))
0.24351
"""
d = np.zeros((lmax + 1))
if s1 < abs(s2):
print("error spins, s1<|s2|")
return
# Conv: sign = -1 if (s1 + s2) and 1 else 1
sign = (-1)**(s1 - s2)
fs1 = math.factorial(2.0 * s1)
fs1ps2 = math.factorial(1.0 * s1 + s2)
fs1ms2 = math.factorial(1.0 * s1 - s2)
num = (1.0 + z)**(0.5 * (s1 + s2)) * (1.0 - z)**(0.5 * (s1 - s2))
# Initialise the recursion (l = s1 + 1)
d[s1] = sign / 2.0**s1 * np.sqrt(fs1 / fs1ps2 / fs1ms2) * num
l1 = s1 + 1.0
rhoSSL1 = np.sqrt((l1 * l1 - s1 * s1) * (l1 * l1 - s2 * s2)) / l1
d[s1+1] = (2 * s1 + 1.0)*(z - s2 / (s1 + 1.0)) * d[s1] / rhoSSL1
# Build the recursion for l > s1 + 1
for l in np.arange(s1 + 1, lmax, 1):
l1 = l + 1.0
numSSL = (l * l * 1.0 - s1 * s1) * (l * l * 1.0 - s2 * s2)
rhoSSL = np.sqrt(numSSL) / (l * 1.0)
numSSL1 = (l1 * l1 - s1 * s1) * (l1 * l1 - s2 * s2)
rhoSSL1 = np.sqrt(numSSL1) / l1
numd = (2.0 * l + 1.0) * (z - s1 * s2 / (l * 1.0) / l1) * d[l]
d[l+1] = (numd - rhoSSL * d[l-1]) / rhoSSL1
return d