""" Common tools used by the component sepration pipeline.
"""
import healpy
import spherelib
import spherelib.astro as astro
import planckdata as data
import spherelib.astro as astro
import numpy as np
import glob
import os.path as op
import matplotlib.pyplot as plt
import smica
import os
import cPickle
from matplotlib import rc
[docs]def get_bins(lmin, lmax):
""" Return Xpol TT bins between lmin and lmax.
"""
bins = np.loadtxt(data.get_binfile())[:,1:3]
if lmin<bins[0,0]:
bins = np.vstack((np.array([lmin,bins[0,0]]), bins))
else:
imin = np.where(bins[:,0]>lmin)[0][0]
bins = bins[imin:,:]
if lmax>bins[-1,1]:
bins = np.vstack((bins, np.array([bins[-1,1], lmax]), bins))
else:
imax = np.where(bins[:,1]<lmax)[0][-1]
bins = bins[:imax+1,:]
return bins
[docs]def get_bins_plot():
""" Provide standard bins used for plots.
This is an arbitrary choice of bin size, only used for plots.
"""
fn_bins = data.get_binfile()
bins = np.loadtxt(fn_bins)[5:23,1:3]
bins = np.vstack((np.array([5,9]),bins))
bins = np.vstack((bins,spherelib.uniformbin(150,199,20)))
bins = np.vstack((bins,spherelib.uniformbin(200,299,40)))
bins = np.vstack((bins,spherelib.uniformbin(300,1499,80)))
bins = np.vstack((bins,np.loadtxt(fn_bins)[91:,1:3]))
bins = np.vstack((bins, np.array([3501,4000])))
lmean = np.squeeze(spherelib.get_lmean(bins))
return bins,lmean
[docs]def get_reference_cmb_powspec(name, cmb_unit="uK_CMB", ref_beam=5, year='2015'):
""" Return Planck 2013 best fit cmb power spectra.
"""
dname = dict({"TT":1, "EE":3, "BB":4})
conv = astro.convfact(freq=1, fr="uK_CMB", to=cmb_unit)
C = np.loadtxt(data.get_lambdaCDMfile(year=year))
EE = np.zeros((C[-1,0]+1))
EE[(C[:,0]).astype(np.int32)] = C[:,dname[name]]/(C[:,0]*(C[:,0]+1)/(2*np.pi))
b5 = np.squeeze(spherelib.gaussbeam(ref_beam, len(EE)-1))
EE *= b5**2
EE *= conv
return EE
[docs]def get_empirical_cmb_powspec(name, cmb_unit="uK_CMB", ref_beam=5, year='2015'):
""" Return Planck 2015 empirical cmb power spectra.
"""
dname = dict({"TT":3, "TE":5, "EE":7})
conv = astro.convfact(freq=1, fr="uK_CMB", to=cmb_unit)
C = spherelib.read_table(data.get_PlanckClfile(year=year), hdu=dname[name])
EE = np.zeros((C[-1,0]+1))
EE[(C[:,0]).astype(np.int32)] = C[:,1]/(C[:,0]*(C[:,0]+1)/(2*np.pi))
b5 = np.squeeze(spherelib.gaussbeam(ref_beam, len(EE)-1))
EE *= b5**2
EE *= conv
return EE
[docs]def make_psmask(preprocdir, nside, lstname='all'):
""" Make a common binary point source mask
"""
if lstname=='all':
msk_fn = sorted(glob.glob(op.join(preprocdir, "*/psmask*.fits.gz")))
else:
msk_fn = list()
for name in lstname:
msk_fn.append(glob.glob(op.join(preprocdir, "*/psmask*%s*.fits.gz"%name))[0])
psmask = np.ones((12*nside*nside))
for fn in msk_fn:
ppmask = healpy.read_map(fn)
psmask *= healpy.ud_grade(ppmask, nside)
psmask[psmask!=0] = 1
return psmask
[docs]def load_ibl(lstibl, lmax):
""" Read beam correction functions from files.
"""
nmap = len(lstibl)
ibl = np.zeros((nmap, lmax+1))
for m in range(nmap):
ibl[m,:] = healpy.read_cl(lstibl[m])
return ibl
[docs]def compute_cov(lstalm, lmax, mask=None, ibl=None):
""" Compute covariance from alm.
A correction for beam and mask effect is performed is options are
given.
"""
nmap = len(lstalm)
nalm = (lmax+1)*(lmax+2)/2
alm = np.zeros((nmap,nalm), dtype=np.complex128)
nmap = len(lstalm)
for m in range(nmap):
alm[m,:] = healpy.read_alm(lstalm[m])
lstats, lnmodes = spherelib.alm2cov(alm)
if ibl is not None:
for i in range(nmap):
for j in range(i, nmap):
iblij = np.sqrt(ibl[i,:]*ibl[j,:])
lstats[i,j,:] *= iblij**2
lstats[j,i,:] = lstats[i,j,:]
if mask is not None:
spherelib.mask_correction (lstats, mask, corr='fsky', nmode=lnmodes)
return lstats,lnmodes
[docs]def cross_beam_correction(lstats, ibl, names):
"""
"""
pass
[docs]def compute_noise_from_jk(lsthr1, lsthr2, lmax, fsky=1, ibl=None, yr=False):
""" Compute noise power spectra from alm.
A correction for beam and mask effect is performed is options are
given.
"""
nmap = len(lsthr1)
if ibl is None:
ibl = np.ones((nmap, lmax+1))
Nl= np.zeros((nmap, lmax+1))
for hr1,hr2,ind in zip(lsthr1,lsthr2,range(nmap)):
alm = healpy.read_alm(hr1) - healpy.read_alm(hr2)
if yr:
if '030' in hr1 or '044' in hr1 or '070' in hr1:
cl = healpy.alm2cl(alm)/4.0/fsky*ibl[ind,:]**2 ## yr1+yr3 / yr2+yr4
else:
cl = healpy.alm2cl(alm)/5.0/fsky*ibl[ind,:]**2
else:
cl = healpy.alm2cl(alm)/4.0/fsky*ibl[ind,:]**2
Nl[ind,:] = cl
return Nl
[docs]def plot_cov(stats, N=None, filename=None, lmean=None):
""" Make a plot of the auto spectra of the covariance matrice.
Noise is superimposed if provided.
"""
nmap = stats.shape[0]
lmax = stats.shape[2]-1
colors = ["b", "r", "g", "k", "c", "y", "m", "orange", "grey", "lightblue", "lightgreen", "b", "r", "g", "k" ,"c", "y", "m", "orange", "grey","b", "r", "g", "k", "c", "y", "m", "orange", "grey", "lightblue", "lightgreen", "b", "r", "g", "k" ,"c", "y", "m", "orange", "grey"]
asq = np.array([0,25,250,750,1500,2500,4000],dtype=int)
if lmean is None:
lmean = np.arange(2, lmax+1)
rg = slice(2,lmax+1)
else:
rg = slice(0,lmax+1)
imax = asq.searchsorted(lmean[-1])
imin = asq.searchsorted(lmean[0])
asq = np.hstack((asq[imin:imax], np.array([lmean[-1]])))
plt.figure()
for m in range(nmap):
plt.semilogy(np.sqrt(lmean), np.squeeze(stats[m,m,rg]), color=colors[m])
if N is not None:
for m in range(nmap):
plt.semilogy(np.sqrt(lmean), N[m,rg], linestyle="--",color=colors[m])
plt.xticks(np.sqrt(asq), [str(i) for i in asq])
if filename is not None:
plt.savefig(filename)
[docs]def plot_ev_decomp(stats, N, lmean=None, filename=None):
""" Plot eigen value decomposition
"""
nmap,Q = N.shape
EIG = np.zeros((Q, nmap))
for q in range(Q):
R12 = np.linalg.inv(np.diag(np.sqrt(N[:,q])))
Rnorm = np.dot(np.dot(R12, stats[:,:,q]), R12)
(EIG[q,:], jk) = np.linalg.eig(Rnorm)
EIG[q,:] = np.sort(EIG[q,:])
asq = np.array([0,25,250,750,1500,2500,4000],dtype=int)
if lmean is None:
lmean = np.arange(5, Q)
rg = slice(5,Q)
else:
rg = slice(0,Q)
imax = asq.searchsorted(lmean[-1])
imin = asq.searchsorted(lmean[0])
asq = np.hstack((asq[imin:imax], np.array([lmean[-1]],dtype=int)))
plt.figure()
plt.subplot(1,2,1)
plt.semilogy(np.sqrt(lmean), EIG[rg,:])
plt.xticks(np.sqrt(asq), [str(i) for i in asq])
plt.axis([None,None,1e-1,None])
plt.subplot(1,2,2)
plt.semilogy(np.sqrt(lmean), EIG[rg,:]-1)
plt.xticks(np.sqrt(asq), [str(i) for i in asq])
if filename is not None:
plt.savefig(filename)
[docs]def get_ilc_filter(Rl, A):
""" A^T R^-1 ( A^T R^-1 A )^-1
"""
nmap,nmap,lmax = Rl.shape
lmax = lmax-1
nc = A.shape[1]
Fh = np.zeros((nc,nmap,lmax+1))
for ell in range(lmax+1):
iRn = np.linalg.pinv(Rl[:,:,ell])
#WL[:,q] = dot(dot(linalg.inv(dot(dot(A.T,iRn),A)),A.T),iRn)
Fh[...,ell] = np.dot(np.dot(np.linalg.pinv(np.dot(np.dot(A.T,iRn),A)),A.T),iRn)
Fh[...,0:2] = 0
return Fh
[docs]def get_fg_filter(Rl, ics, irs, A, Pl):
""" (A P A^T) R^-1
"""
nmap,nmap,lmax = Rl.shape
lmax = lmax-1
n = A.shape[1]
nc = len(ics)
F = np.zeros((nc,nmap,lmax+1))
snr = np.zeros((nc,lmax+1))
for ell in range(lmax+1):
iRn = np.linalg.pinv(Rl[:,:,ell])
for ic,ir in zip(ics,irs):
Pli = Pl[ic,:,ell].reshape(1,n)
Ai = A[:,ic].reshape(nmap,1)
snr[...,ell] = np.dot(Pli, np.dot(A.T, np.dot(iRn, Ai)))
F[...,ell] = np.dot(Ai, np.dot(Pli, np.dot(A.T, iRn)))[ir,:]
return F,snr
[docs]def smooth_cl (cl, dloverl=0.1):
""" Return smoothed power spectrum.
Power spectrum is averaged over window of size dloverl (in %)
"""
lmax = len(cl) - 1
ell = np.arange(0,lmax+1)
clsmooth = np.zeros(cl.shape);
for l in range(lmax+1):
llo = max(int(np.round(ell[l]*(1-dloverl/2.0))), 0)
lhi = min(int(np.round(ell[l]*(1+dloverl/2.0))), lmax)
if llo==lhi:
lhi=llo+1
clsmooth[l] = np.mean(cl[llo:lhi])
return clsmooth
[docs]def unbin_R(R, Rl_high, bins, lcst=None, lilc=None):
""" Provide a set of covariance matrices which are smoothed over multipole range.
High ell regime is read from Rl_high starting from lilc multipole.
Low ell regime can be forced to be constant from 0 to lcst.
Matrices are supposed to be symmetric.
"""
ll = np.squeeze(spherelib.get_lmean(bins))
lmin = bins[0,0]
lmax = bins[-1,1]
if lcst is None:
lcst = lmin
if lilc is None:
lilc = lmax
nmap,nmap,Q = R.shape
lmax_sht = Rl_high.shape[2]-1
Rl = np.zeros((nmap,nmap,lmax_sht+1))
for i in range(nmap):
for j in range(i, nmap):
Rl[i,j,lmin:lmax] = np.interp(np.arange(lmin, lmax), ll, R[i,j,:])
Rl[i,j,0:lmin] = Rl[i,j,lmin+1]
Rl[i,j,0:lcst] = Rl[i,j,lcst+1]
rls = smooth_cl (Rl_high[i,j,:], dloverl=0.3)
Rl[i,j,lilc:] = rls[lilc:]
Rl[j,i,:] = Rl[i,j,:]
return Rl
[docs]def plot_filter(F, names, fields, filename=None):
""" Make a plot of detector weight as a function of the multipole.
"""
colors = ["c", "m", "y", "r", "b", "g", "k", "c", "m", "y", "r", "b", "g", "k"]
bins_plot,lmean = get_bins_plot()
names = [n.replace("_", "") for n in names]
nfield = len(fields)
nmap = len(names)
F = F.reshape(nfield, nmap, -1)
plt.figure()
for j in range(nfield):
ll = []
plt.subplot(1,nfield,j+1)
plt.title(fields[j])
for m in range(nmap):
data = spherelib.bin_cl(F[j,m,:], bins_plot)
if names[m][:-1] in ["100", "143", "217", "353"] and fields[j]==names[m][-1]:
ll.append(plt.plot(np.sqrt(lmean), data, color=colors[m])[0])
else:
plt.plot(np.sqrt(lmean), data, color=colors[m])
asq = [0,25,250,750,1500,2500,4000]
plt.xticks(np.sqrt(asq), [str(i) for i in asq])
if j==0:
leg = plt.legend(ll, ["100", "143", "217", "353"], loc="upper center")
for t in leg.get_texts():
t.set_fontsize('small') # the legend text fontsize
plt.ylabel(r"$\mathbf{w_i} \/ \left[ \mu \mathrm{K} / \mu \mathrm{K_i^{RJ}} \right]$",fontsize=20)
if filename is not None:
plt.savefig (filename)
[docs]def plot_noise_decomp(F, N, names, fields, filename=None):
""" Make a plot of the detector noise contribution as a function of the multipole.
"""
colors = ["c", "m", "y", "r", "b", "g", "k", "c", "m", "y", "r", "b", "g", "k"]
bins_plot,lmean = get_bins_plot()
names = [n.replace("_", "") for n in names]
nfield = len(fields)
nmap = len(names)
lmax = F.shape[-1]-1
plt.figure()
for j in range(nfield):
ll = []
sumn = np.zeros((lmax+1))
plt.subplot(1,nfield,j+1)
plt.title(fields[j])
for m in range(nmap):
data = F[j,m,:]**2*N[m,:]
sumn += data
data = spherelib.bin_cl(data, bins_plot)
if names[m][:-1] in ["100", "143", "217", "353"] and fields[j]==names[m][-1]:
ll.append(plt.semilogy(np.sqrt(lmean), data, color=colors[m])[0])
else:
plt.semilogy(np.sqrt(lmean), data, color=colors[m])
data = spherelib.bin_cl(sumn, bins_plot)
ll.append(plt.semilogy(np.sqrt(lmean), data, lw=2, color='purple')[0])
asq = [0,25,250,750,1500,2500,4000]
plt.xticks(np.sqrt(asq), [str(i) for i in asq])
if j==0:
leg = plt.legend(ll, ["100", "143", "217", "353", "sum"], loc="lower left")
for t in leg.get_texts():
t.set_fontsize('small') # the legend text fontsize
plt.ylabel(r"$\mathbf{w_i}^2 N_i \/ \left[ \mu \mathrm{K}^2 \right]$",fontsize=20)
plt.xlabel("$\mathcal{l}$",fontsize=20)
plt.axis([None,None,1e-12,None])
if filename is not None:
plt.savefig (filename)
if __name__ == '__main__':
pass