Commit 319a0d26 authored by Julien's avatar Julien
Browse files

Restructuration: split the code into sub-modules. Apply PEP8, but not everywhere...

parent 167c0431
from pylab import *
import astropy.io.fits as fits
import xqml
"""
Test script for xQML
Author: Vanneste
"""
from __future__ import division
import timeit
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
from pylab import *
import astropy.io.fits as fits
from xqml import xQML
from simulation import getstokes, muKarcmin2var, GetBinningMatrix
from libcov import compute_ds_dcb
from xqml_utils import progress_bar
nside = 8
lmax = 3*nside-1
Slmax = 3*nside-1
deltal = 1
clth = hp.read_cl( 'planck_base_planck_2015_TTlowP.fits')
lth = arange(2,lmax+1)
if __name__ == "__main__":
nside = 8
lmax = 3*nside-1
Slmax = 3*nside-1
deltal = 1
clth = hp.read_cl('planck_base_planck_2015_TTlowP.fits')
lth = arange(2,lmax+1)
ellbins = arange(2,lmax+2,deltal)
ellbins[-1] = lmax+1
ellbins = arange(2,lmax+2,deltal)
ellbins[-1] = lmax+1
P, Q, ell, ellval = xqml.GetBinningMatrix(ellbins, lmax)
nbins=len(ellbins)-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)
#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)
fwhm = 0.5
bl=hp.gauss_beam(deg2rad(fwhm), lmax=Slmax+1)
allStoke, der, ind = xqml.getstokes(polar=True,temp=False,EBTB=False)
nder=len(der)
allStoke, der, ind = getstokes(polar=True,temp=False,EBTB=False)
nder=len(der)
muKarcmin = 0.1
muKarcmin = 0.1
pixvar = 2*xqml.muKarcmin2var(muKarcmin, nside)
varmap = ones((2*npix))*pixvar
NoiseVar = np.diag(varmap)
pixvar = 2 * muKarcmin2var(muKarcmin, nside)
varmap = ones((2*npix))*pixvar
NoiseVar = np.diag(varmap)
cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
noise = (randn(len(varmap)) * varmap**.5).reshape(2,-1)
dm = cmb[1:,mask] + noise
cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
noise = (randn(len(varmap)) * varmap**.5).reshape(2,-1)
dm = cmb[1:,mask] + noise
###################################### Compute ds_dcb ######################################
ip=arange(hp.nside2npix(nside))
ipok=ip[mask]
###################################### Compute ds_dcb ######################################
ip=arange(hp.nside2npix(nside))
ipok=ip[mask]
Pl, S = xqml.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)
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.xQML(mask,ellbins, clth, Pl=Pl, fwhm=fwhm)
esti.construct_esti( NoiseVar, NoiseVar)
cl = esti.get_spectra( dm, dm)
V = esti.get_covariance()
###################################### Compute spectra ######################################
esti = xQML(mask,ellbins, clth, Pl=Pl, fwhm=fwhm)
esti.construct_esti( NoiseVar, NoiseVar)
cl = esti.get_spectra( dm, dm)
V = esti.get_covariance()
###################################### Construct MC ######################################
allcl = []
esti = xqml.xQML(mask,ellbins, clth, Pl=Pl, fwhm=fwhm)
esti.construct_esti( NoiseVar, NoiseVar)
start = timeit.default_timer()
for n in np.arange(100):
xqml.progress_bar( n, 100,timeit.default_timer()-start)
cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, new=True, verbose=False)
dm = cmb[1:,mask] + (randn(2*npix)*sqrt(varmap)).reshape( 2, npix)
allcl.append(esti.get_spectra( dm, dm))
###################################### Construct MC ######################################
allcl = []
esti = xQML(mask,ellbins, clth, Pl=Pl, fwhm=fwhm)
esti.construct_esti( NoiseVar, NoiseVar)
start = timeit.default_timer()
for n in np.arange(100):
progress_bar( n, 100,timeit.default_timer()-start)
cmb = hp.synfast( clth, nside, fwhm=deg2rad(fwhm), pixwin=True, 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( 2,1,1)
plot( lth, clth.transpose()[lth,1:3], 'k')
plot( ellval, mean( allcl,0).transpose(), 'r')
plot( ellval, mean( allcl,0).transpose() + std( allcl,0).transpose(), 'r--')
plot( ellval, mean( allcl,0).transpose() - std( allcl,0).transpose(), 'r--')
semilogy()
subplot( 2,1,2)
cosmic = sqrt( 2./(2*lth+1))/mean(mask) * clth[1:3,lth]
plot( lth, cosmic.transpose(), 'k')
plot( ellval, std( allcl,0).transpose(), 'r')
plot( ellval, sqrt(diag(V)).reshape(nder,-1).transpose(), 'b')
semilogy()
figure()
subplot( 2,1,1)
plot( lth, clth.transpose()[lth,1:3], 'k')
plot( ellval, mean( allcl,0).transpose(), 'r')
plot( ellval, mean( allcl,0).transpose() + std( allcl,0).transpose(), 'r--')
plot( ellval, mean( allcl,0).transpose() - std( allcl,0).transpose(), 'r--')
semilogy()
subplot( 2,1,2)
cosmic = sqrt( 2./(2*lth+1))/mean(mask) * clth[1:3,lth]
plot( lth, cosmic.transpose(), 'k')
plot( ellval, std( allcl,0).transpose(), 'r')
plot( ellval, sqrt(diag(V)).reshape(nder,-1).transpose(), 'b')
semilogy()
show()
show()
"""
Set of routines to ...
Comment: what is Expd()? It is defined nowhere in the code...
Author: Vanneste
"""
from __future__ import division
import numpy as np
def Pl(ds_dcb):
"""
Reshape ds_dcbin Pl
Parameters
----------
ds_dcb : ???
???
Returns
----------
???
"""
nnpix = np.shape(ds_dcb)[-1]
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):
"""
Compute correlation matrix S = sum_l Pl*Cl
Parameters
----------
Clth : ???
???
Returns
----------
???
"""
if EBTB:
xx = ['TT', 'EE', 'BB', 'TE', 'EB', 'TB']
ind = [0, 1, 2, 3, 4, 5]
else:
xx = ['TT', 'EE', 'BB', 'TE']
ind = [0, 1, 2, 3]
if not temp:
allStoke = ['Q', 'U']
if EBTB:
xx = ['EE', 'BB', 'EB']
ind = [1, 2, 5]
else:
xx = ['EE', 'BB']
ind = [1, 2]
if not polar:
allStoke = ['I']
xx = ['TT']
ind = [0]
clth = Clth[ind][:, 2: int(ellbins[-1])].flatten()
return np.sum(Pl * clth[:, None, None], 0)
def El(invCAA, invCBB, Pl, Bll=None, expend=False):
"""
Compute El = invCAA.Pl.invCBB
Parameters
----------
invCAA : ???
???
Returns
----------
???
"""
# Tegmark B-matrix useless so far)
if Bll == None:
Bll = np.diagflat((np.ones(len(Pl))))
lmax = len(Pl) * 2**int(expend)
lrange = np.arange(lmax)
npix = len(invCAA)
# triangular shape ds_dcb
if expend:
El = np.array([np.dot(np.dot(invCAA, Expd(Pl,l)), invCBB) for l in lrange]).reshape((lmax,npix,npix))
else:
El = np.array([np.dot(np.dot(invCAA, Pl[l]), invCBB) for l in lrange]).reshape((lmax,npix,npix))
return El
def ElLong(invCAA, invCBB, Pl, Bll=None, expend=False):
"""
Compute El = invCAA.Pl.invCBB
Parameters
----------
invCAA : ???
???
Returns
----------
???
"""
# Tegmark B-matrix useless so far)
if Bll == None:
Bll = np.diagflat((np.ones(len(Pl))))
lmax = len(Pl) * 2**int(expend)
lrange = np.arange(lmax)
npix = len(invCAA)
# triangular shape ds_dcb
if expend:
for l in lrange:
Pl[l] = np.dot(np.dot(invCAA, Expd(Pl, l)), invCBB).reshape((npix, npix ))
else:
for l in lrange:
Pl[l] = np.dot(np.dot(invCAA, Pl[l]), invCBB).reshape((npix, npix ))
def CrossFisherMatrix(El, Pl, expend=False):
"""
Compute cross (or auto) fisher matrix Fll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
Parameters
----------
El : ???
???
Returns
----------
???
"""
lmax = len(Pl) * 2**int(expend)
lrange = np.arange(lmax)
if expend:
# pas de transpose car symm
FAB = np.array([np.sum(El[il] * Expd(Pl,jl)) for il in lrange for jl in lrange]).reshape(lmax, lmax)
else:
# pas de transpose car symm
FAB = np.array([np.sum(El[il] * Pl[jl]) for il in lrange for jl in lrange]).reshape(lmax, lmax)
return FAB
def CrossWindowFunction(El, Pl):
"""
Compute Tegmark cross (or auto) window matrix Wll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
Equivalent to Fisher matrix Fll when Tegmark B-matrix = 1
Parameters
----------
El : ???
???
Returns
----------
???
"""
lmax = len(El)
lrange = np.arange((lmax))
# pas de transpose car symm
Wll = np.array([np.sum(El[il] * Pl[jl]) for il in lrange for jl in lrange]).reshape(lmax, lmax)
return Wll
def CrossWindowFunctionLong(invCAA, invCBB, Pl):
"""
Compute Tegmark cross (or auto) window matrix Wll = Trace[invCAA.Pl.invCBB.Pl] = Trace[El.Pl]
Equivalent to Fisher matrix Fll when Tegmark B-matrix = 1
Parameters
----------
invCAA : ???
???
Returns
----------
???
"""
lmax = len(Pl)
lrange = np.arange((lmax))
# Pas de transpose car symm
Wll = np.array([np.sum(np.dot(np.dot(invCAA, Pl[il]), invCBB) * Pl[jl]) for il in lrange for jl in lrange]).reshape(lmax, lmax)
return Wll
def CrossMisherMatrix(El, CAA, CBB):
"""
Compute matrix MAB = Trace[El.CAA.El.CBB] (see paper for definition)
Parameters
----------
El : ???
???
Returns
----------
???
"""
lmax = len(El)
lrange = np.arange(lmax)
El_CAA = np.array([np.dot(CAA,El[il]) for il in lrange])
El_CBB = np.array([np.dot(CBB,El[il]) for il in lrange])
MAB = np.array([np.sum(El_CAA[il] * El_CBB[jl].T) for il in lrange for jl in lrange]).reshape(lmax, lmax)
return MAB
def CrossGisherMatrix(El, CAB):
"""
Compute matrix GAB = Trace[El.CAB.El.CAB] (see paper for definition)
Parameters
----------
El : ???
???
Returns
----------
???
"""
lmax = len(El)
lrange = np.arange(lmax)
El_CAB = np.array([np.dot(CAB,El[il]) for il in lrange])
GAB = np.array([np.sum(El_CAB[il] * El_CAB[jl].T) for il in lrange for jl in lrange]).reshape(lmax, lmax)
return GAB
def CrossGisherMatrixLong(El, CAB):
"""
Compute matrix GAB = Trace[El.CAB.El.CAB] (see paper for definition)
Parameters
----------
El : ???
???
Returns
----------
???
"""
lmax = len(El)
lrange = np.arange(lmax)
GAB = np.array([np.sum(np.dot(CAB, El[il]) * np.dot(CAB,El[jl]).T) for il in lrange for jl in lrange]).reshape(lmax, lmax)
return GAB
def yQuadEstimator(dA, dB, El):
"""
Compute pre-estimator 'y' such that Cl = Fll^-1 . yl
Parameters
----------
dA : ???
???
Returns
----------
???
"""
npix = len(dA)
lrange = np.arange((len(El)))
y = np.array([dA.dot(El[l]).dot(dB) for l in lrange])
return y
def ClQuadEstimator(invW, y):
"""
Compute estimator 'Cl' such that Cl = Fll^-1 . yl
Parameters
----------
invW : ???
???
Returns
----------
???
"""
Cl = np.dot(invW, y)
return Cl
def biasQuadEstimator(NoiseN, El):
"""
Compute bias term bl such that Cl = Fll^-1 . ( yl + bias)
Parameters
----------
NoiseN : ???
???
Returns
----------
???
"""
lrange = np.arange((len(El)))
return np.array([np.sum(NoiseN * El[l]) for l in lrange])
def blEstimatorFlat(NoiseN, El):
"""
Compute bias term bl such that Cl = Fll^-1 . ( yl + bl)
Not to be confonded with beam function bl(fwhm)
Parameters
----------
NoiseN : ???
???
Returns
----------
???
"""
lrange = np.arange((len(El)))
return np.array([np.sum(NoiseN * np.diag(El[l])) for l in lrange])
def CovAB(invWll, GAB):
"""
Compute analytical covariance matrix Cov(Cl, Cl_prime)
Parameters
----------
invWll : ???
???
Returns
----------
???
"""
covAB = np.dot(np.dot(invWll, GAB), invWll.T) + invWll
return covAB
"""
Set of routines to ...
Author: Vanneste
"""
from __future__ import division
import numpy as np
def polrotangle(ri, rj):
"""
???
Parameters
----------
ri : ???
???
rj : ???
???
Returns
----------
cos2a : 1D array of floats
???
sin2a : 1D array of floats
???
"""
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
This diff is collapsed.
"""
Set of routines to ...
Author: Vanneste
"""
from __future__ import division
import numpy as np
import healpy as hp
from scipy import sparse
def muKarcmin2var(muKarcmin, nside):
"""
Return pixel variance for a given nside and noise level [1e-6 K . arcmin]
Parameters
----------
muKarcmin : ???
???
...
Returns
----------
??? : ???
???
"""
pixarea = hp.nside2pixarea(nside, degrees = True)
varperpix = (muKarcmin*1e-6/60.)**2/pixarea
return varperpix
def pixvar2nl(pixvar, nside):
"""
Return noise spectrum level for a given nside and pixel variance
Parameters
----------
pixvar : ???
???
...
Returns
----------
??? : ???
???
"""
return pixvar*4.*np.pi/(12*nside**2.)
def getNl(pixvar, nside, nbins):
"""
Return noise spectrum for a given nside and pixel variance
Parameters
----------
pixvar : ???
???
...