import healpy
import spherelib
import spherelib.astro as astro
import planckdata as data
import numpy as np
import glob
import os.path as op
import matplotlib.pyplot as plt
import smica
import os
import cPickle
from pipetools import *
from matplotlib import rc
[docs]def load_from_dir(dirname):
"""
"""
sfile = op.join(dirname, "hR_binned_full.txt")
s = np.loadtxt(sfile)
m2,Q = s.shape
w = np.loadtxt(op.join(dirname, "w_binned_full.txt"))
s = s.reshape(np.sqrt(m2),np.sqrt(m2), Q)
model = cPickle.load(file(op.join(dirname, "model_full.pkl")))["model"]
return s, w, model
[docs]def get_cmb(names, unit_to, bins, polar=False, calib=None, cmb_unit="uK_CMB", ref_beam=5, fixed=True):
""" Return a smica cmb component of known mixing matrix.
"""
# calibration
nmap = len(names)
Q = bins.shape[0]
a = np.loadtxt(calib) if calib is not None else np.ones((8))
dcalib = dict(zip([30,44,70,100,143,217,353,545],a))
dcalib[30] = 1 # no recalibration of the 30GHz
dcalib[857] = 1 # no recalibration of the 857GHz
# unit convert
if polar:
acmb = np.zeros((nmap,2))
for m in range(nmap/2): #assume E B for all channel
freq = data.get_freq(names[m])
acmb[m,0] = astro.convfact(freq=freq*1e9, fr=cmb_unit, to=unit_to)*dcalib[freq]
acmb[m+nmap/2,1] = acmb[m,0]
cmb = smica.SourceND (nmap, Q, 2, name='cmb')
else:
acmb = np.zeros((nmap,1))
for m in range(nmap):
freq = data.get_freq(names[m])
acmb[m,0] = astro.convfact(freq=freq*1e9, fr=cmb_unit, to=unit_to)*dcalib[freq]
cmb = smica.Source1D (nmap, Q, name='cmb')
if fixed:
cmb.set_mixmat(acmb, fixed='all')
elif not polar:
if '143_T' in names:
i143 = names.index('143_T')
elif '143_5_T' in names:
i143 = names.index('143_5_T')
i857 = names.index('857_T')
afix = np.zeros((nmap,1))
afix[i143] = 1
afix[i857] = 1
cmb.set_mixmat(acmb, fixed=afix)
else:
print "warning: calib in polar is not available"
#powspec
if polar:
EE = get_reference_cmb_powspec('EE', cmb_unit=cmb_unit, ref_beam=ref_beam)
BB = get_reference_cmb_powspec('BB', cmb_unit=cmb_unit, ref_beam=ref_beam)
cqcmb = np.zeros((2,2,Q))
cqcmb[0,0,:] = spherelib.bin_cl(EE,bins)
cqcmb[1,1,:] = spherelib.bin_cl(BB,bins)
cmb.set_powspec (cqcmb)
else:
TT = get_reference_cmb_powspec('TT', cmb_unit=cmb_unit, ref_beam=ref_beam)
cmbcq = spherelib.bin_cl(TT, bins)
cmb.set_powspec(cmbcq)
return cmb
[docs]def get_gal(nmap, dim, Q, polar=False, polar2D=False):
""" Return a smica N dimensional galactic component
"""
if polar:
nc = 2 #number of column per dim=2 for E and B
agfix = np.zeros((nmap, nc*dim))
for d in range(dim):
oc = d*nc
agfix[nmap/2:,oc] = 1 # E only
agfix[0:nmap/2,oc+1] = 1 # B only
# dust+sync constraint
if polar2D:
if d==0 and dim==2:
agfix[0,oc] = 1 # set 30 to 0 for dust
agfix[nmap/2,oc+1] = 1
elif d==0 and dim==3:
agfix[0:2,oc] = 1
agfix[nmap/2:nmap/2+2,oc+1] = 1
else:
agfix[nmap/2-1,oc] = 1 # set 353 to 0 for sync
agfix[-1,oc+1] = 1
dim = dim*nc
gal = smica.SourceND(nmap, Q, dim, name='gal')
gal.fix_mixmat(agfix)
else:
gal = smica.SourceND(nmap, Q, dim, name='gal')
return gal
[docs]def get_gal_from(Q, filename=None, model=None, em=None):
""" Return a smica N dimensional galactic component
of fixed mixmat read from file or model.
"""
if model is None:
model = cPickle.load(file(filename))["model"]
if em is None:
em = cPickle.load(file(filename))["em"]
An = model.get_comp_by_name("gal").mixmat()
eAn = em.get_comp_by_name("gal").mixmat()
nmap,dim= An.shape
gal= smica.SourceND (nmap, Q, dim, name='gal')
gal.set_mixmat(An, fixed='all')
return gal,eAn
[docs]def get_noise(N, fixed=False):
""" Return a smica noise component
"""
nmap,Q = N.shape
noise = smica.NoiseAmpl(nmap, Q, name='noise')
noise.set_powspec(N, fixed="all")
if fixed:
noise.set_ampl(np.ones((nmap,1)), fixed="all")
return noise
[docs]def is_mixmat_fixed(model):
return np.sum(model.get_comp_by_name("gal")._mixmat.get_mask())==0
[docs]def fit_model_to_cov(model, stats, nmodes, maxiter=50, noise_fix=False, afix=None, logger=None):
""" Fit the model to empirical covariance.
"""
cg_maxiter = 1
cg_eps = 1e-12
Q = len(nmodes)
nmap = stats.shape[0]
cmb,gal,noise = model._complist
# find a starting point
acmb = model.get_comp_by_name("cmb").mixmat()
cmbcq = np.squeeze(model.get_comp_by_name("cmb").powspec())
polar = True if acmb.shape[1]==2 else False
if not is_mixmat_fixed(model):
afix = 1-cmb._mixmat.get_mask() #0:free 1:fixed
model.ortho_subspace(stats, nmodes, acmb, qmin=0, qmax=Q)
cmb.set_mixmat(acmb, fixed=afix)
model.quasi_newton(stats, nmodes)
cmb.set_powspec (cmbcq)
model.close_form(stats)
cmb.set_powspec (cmbcq)
mm = model.mismatch(stats, nmodes, exact=True)
mmG = model.mismatch(stats, nmodes, exact=True)
# start CG/close form
for i in range(maxiter):
# fit mixing matrix
gal.fix_powspec("all")
cmbfix = "all" if polar else "null"
cmb.fix_powspec(cmbfix)
model.conjugate_gradient (stats, nmodes,maxiter=cg_maxiter, avextol=cg_eps)
# fit power spectra
gal.fix_powspec("null")
cmbfix = "null" if polar else "all"
cmb.fix_powspec(cmbfix)
if i==maxiter/2 and not noise_fix: # fit also noise at some point
Nt = noise.powspec()
noise = smica.NoiseDiag(nmap, Q, name='noise')
model = smica.Model([cmb, gal, noise])
noise.set_powspec(Nt)
model.close_form(stats)
# compute new mismatch
mm2 = model.mismatch(stats, nmodes, exact=True)
mm2G = model.mismatch(stats, nmodes)
gain = np.real(mmG-mm2G)
if gain==0 and i>maxiter/2.0:
break
strtoprint = "iter= % 4i mismatch = %10.5f gain= %7.5f " % (i, np.real(mm2), gain)
if logger is not None:
logger.info(strtoprint)
else:
print strtoprint
mm = mm2
mmG = mm2G
cmb.fix_powspec("null")
gal.fix_powspec("null")
return model
[docs]def plot_mismatch(model, stats, nmodes, ll, filename=None):
""" Plot spectral mismatch
"""
Q = len(nmodes)
ellmax = ll[-1]
nmap = stats.shape[0]
mm = np.zeros((Q))
for z in range(Q):
mm[z] = int(np.real(model.mismatch(stats[:,:,z], nmodes[z], bin=z, exact=True)))
dof = 0.5 * ( (nmap*(nmap+1))/2.0 - model.get_dim()/model.nbin)
smm = np.sum(mm)
plt.figure()
plt.step(ll, mm, color='r')
plt.plot(np.ones((ellmax+1))*dof, color='k')
plt.title("Total mismatch %2.3e" % smm)
if filename is not None:
plt.savefig(filename)
[docs]def plot_relative_diff(model, Stat, NModes, channelnames, lmean, filename=None):
""" Plot relative difference between model and data covariance matrices.
"""
channelnames = [n.replace("_", "") for n in channelnames]
plt.figure()
nchannels = Stat.shape[0]
SumCovMat = model.covariance()
mode = 'model'
for i in range(nchannels):
for j in range(i,nchannels):
plt.subplot(nchannels,nchannels,i*nchannels+j+1)
if (mode=='model'):
dcl_over_cl = (Stat[i,j,:]-SumCovMat[i,j,:])/SumCovMat[i,j,:]
else:
dcl_over_cl = (Stat[i,j,:]-SumCovMat[i,j,:])/Stat[i,j,:]
if (i==j):
cov = np.sqrt(2./NModes)
else:
if (mode == 'model'):
cov = np.sqrt(1./NModes * (1+SumCovMat[i,i,:]*SumCovMat[j,j,:]/SumCovMat[i,j,:]**2))
else:
cov = np.sqrt(1./NModes * (1+Stat[i,i,:]*Stat[j,j,:]/Stat[i,j,:]**2))
plt.ylim(-0.3,0.3)
plt.errorbar(lmean,dcl_over_cl,yerr=cov, linewidth=0.5,fmt='+')
# Horizontal line at zero
xmin,xmax=plt.xlim()
plt.hlines(0,xmin,xmax)
xticklabels = plt.getp(plt.gca(), 'xticklabels')
yticklabels = plt.getp(plt.gca(), 'yticklabels')
if (i==0):
plt.title(channelnames[j],fontsize=8)
if (j==nchannels-1):
ax=plt.gca()
ax.yaxis.set_label_position('right')
plt.ylabel(channelnames[i])
plt.setp(yticklabels, visible=True)
for tick in ax.yaxis.get_major_ticks():
tick.label1.set_fontsize(4)
tick.label2.set_fontsize(4)
ax.yaxis.set_ticks_position("right") # ticks and ticklabels
else:
plt.setp(yticklabels, visible=False)
plt.setp(xticklabels, fontsize=4)
plt.suptitle(r"$(\hat{C}_\ell - C_\ell)/C_\ell$ with naive error bars",fontsize=12)
if filename is not None:
plt.savefig(filename)
[docs]def plot_mismatch_pair(model, stats, nmodes, names, ll, filename=None):
""" Plot spectral mismatch per detector pair.
"""
names = [n.replace("_", "") for n in names]
R = model.covariance()
nmaps,nmaps,Q = stats.shape
fig = plt.figure()
s = np.zeros((2,2,Q))
r = np.zeros((2,2,Q))
for m1 in range(nmaps):
s[0,0,:] = stats[m1,m1,:]
r[0,0,:] = R[m1,m1,:]
for m2 in range(m1+1, nmaps):
s[0,1,:] = stats[m1,m2,:]
s[1,0,:] = stats[m2,m1,:]
s[1,1,:] = stats[m2,m2,:]
r[0,1,:] = R[m1,m2,:]
r[1,0,:] = R[m2,m1,:]
r[1,1,:] = R[m2,m2,:]
mm = np.zeros((Q))
for q in range(Q):
mm[q] = np.real(nmodes[q]*smica.kullback(s[:,:,q],r[:,:,q], exact=True))
smm = sum(mm)
plt.subplot(nmaps-1, nmaps-1,(m1)*(nmaps-1)+(m2-1)+1)
if m2==m1+1:
plt.ylabel(names[m1], fontsize=8)
if m1==0:
plt.title(names[m2], fontsize=8)
plt.step(ll, mm)
plt.plot(ll, np.ones((Q))*(model.get_dim()/(Q*1.0)), color='r')
plt.axis([ll[0],ll[-1], 0, 25])
xticklabels = plt.getp(plt.gca(), 'xticklabels')
yticklabels = plt.getp(plt.gca(), 'yticklabels')
if m2==nmaps-1:
plt.setp(yticklabels, visible=True)
ax = plt.gca()
for tick in ax.yaxis.get_major_ticks():
tick.label1.set_fontsize(4)
tick.label2.set_fontsize(4)
ax.yaxis.set_ticks_position("right") # ticks and ticklabels
else:
plt.setp(yticklabels, visible=False)
plt.setp(xticklabels, fontsize=4)
if filename is not None:
plt.savefig(filename)
[docs]def plot_mismatch_KLD(model, stats, nmodes, names, ll, filename=None):
""" Plot KLD decomposition of the spectral mismatch.
Not working yet.
"""
Q = len(ll)
nmap = stats.shape[0]
Kss = np.zeros((Q))
Knn = np.zeros((Q))
Ksn = np.zeros((Q))
K = np.zeros((Q))
cmb = model.get_comp_by_name("cmb")
gal = model.get_comp_by_name("gal")
dim = gal.dim+cmb.dim
for q in range(Q):
sc = cmb.covariance(bin=q)+gal.covariance(bin=q)
Kss[q], Knn[q], Ksn[q] = smica.KLDdecomp (stats[:,:,q], sc, model.covariance(bin=q)-sc, dimsignal=dim)
K[q] = smica.kullback (stats[:,:,q], model.covariance(bin=q), exact=True)
dof = 0.5 * ( (nmap*(nmap+1))/2.0 - model.get_dim()/model.nbin)
fig = plt.figure()
ind = np.arange(Q) # the x locations for the groups
width = 0.5 # the width of the bars
ax = fig.add_subplot(111)
rects1 = ax.bar(ind, nmodes*Knn, width, color='g')
rects2 = ax.bar(ind, nmodes*Ksn, width, color='r', bottom=nmodes*Knn)
rects3 = ax.bar(ind, nmodes*Kss, width, color='b', bottom=nmodes*(Knn+Ksn))
ax.plot(ind, np.ones((Q))*dof, color='k')
ax.set_ylabel('Mismatch')
ax.set_title('KLD decomposition')
ax.set_xlabel('bin number')
ax.legend( (rects1[0], rects2[0], rects3[0]), ('nn', 'sn', 'ss') )
if filename is not None:
plt.savefig(filename[0])
plt.figure()
i = 1
dg = np.shape(model.mixmat())[1]
titles = ["Total mismatch", "NN", "SS", "NS"]
for k in [K, Knn, Ksn, Kss]:
plt.subplot(2,2,i)
plt.step (ll, nmodes*k)
plt.ylabel('Mismatch')
plt.xlabel('Multipole $\ell$')
if i==1:
plt.plot(ll, dof*np.ones((Q)))
plt.title(titles[i-1])
i = i+1
(x1,x2,y1,y2) = plt.axis()
plt.axis([x1,x2,0,y2])
if filename is not None:
plt.savefig (filename[1])
return nmodes*Kss,nmodes*Knn,nmodes*Ksn
[docs]def plot_noise(model, em, N, ll, names, Nfix=None, filename=None):
""" Plot the fitted noise versus expected noise as predicted by jackknives.
"""
nmaps,Q = N.shape
plt.figure()
noise = model.get_comp_by_name("noise")
for i in range(nmaps):
plt.subplot(int(np.sqrt(nmaps))+1,int(np.ceil(np.sqrt(nmaps))),i+1)
ye = em.get_comp_by_name("noise").powspec()[i,:]
if Nfix is not None:
ye[Nfix[i,:]==1] = 0
plt.plot(np.sqrt(ll), np.ones((Q)), color='k')
plt.errorbar(np.sqrt(ll), noise.powspec()[i,:]/N[i,:], yerr=ye/N[i,:])
asq = [0,25,250,750,1500,2500]
plt.xticks(np.sqrt(asq), [str(j) for j in asq])
plt.title(names[i],fontsize=8)
if filename is not None:
plt.savefig(filename[0])
plt.figure()
for i in range(nmaps):
plt.subplot(int(np.sqrt(nmaps))+1,int(np.ceil(np.sqrt(nmaps))),i+1)
plt.semilogy(np.sqrt(ll), np.squeeze(N[i,:]))
ye = em.get_comp_by_name("noise").powspec()[i,:]
if Nfix is not None:
ye[Nfix[i,:]==1] = 0
plt.errorbar(np.sqrt(ll), np.squeeze(noise.powspec()[i,:]), yerr=ye)
asq = [0,25,250,750,1500,2500]
plt.xticks(np.sqrt(asq), [str(j) for j in asq])
plt.title(names[i],fontsize=8)
if filename is not None:
plt.savefig(filename[1])
[docs]def plot_decomp(model, stats, nmodes, ll, names, filename=None):
""" Plot a decomposition of the data into cmb + fg + noise for each cross spectra.
"""
Rc = model.covariance4D()
R = model.covariance()
aa = model.mixmat()
nmaps = stats.shape[0]
fig = plt.figure()
cols = ['b', 'r', 'c', 'g', 'y']
names = [n.replace("_", "") for n in names]
for m1 in range(nmaps):
for m2 in range(m1, nmaps):
plt.subplot(nmaps,nmaps,m1*nmaps+m2+1)
if m2==m1:
plt.ylabel(names[m1], fontsize=10)
if m1==0:
plt.title(names[m2], fontsize=10)
n = []
lines = []
for c in range(model.ncomp):
lines.append(plt.semilogy(ll, np.squeeze(np.abs(Rc[m1,m2,:,c])), linewidth=0.5, color=cols[c])[0])
n.append(model.get_comp_by_number(c).name)
lines.append(plt.semilogy(ll, np.squeeze(R[m1,m2,:]), linewidth=0.5, color='k')[0])
n.append("sum")
lines.append(plt.semilogy(ll, np.squeeze(stats[m1,m2,:]), linewidth=0.5, color='m')[0])
n.append("obs")
xticklabels = plt.getp(plt.gca(), 'xticklabels')
yticklabels = plt.getp(plt.gca(), 'yticklabels')
if m2==nmaps-1:
plt.setp(yticklabels, visible=True)
ax = plt.gca()
for tick in ax.yaxis.get_major_ticks():
tick.label1.set_fontsize(4)
tick.label2.set_fontsize(4)
ax.yaxis.set_ticks_position("right") # ticks and ticklabels
else:
plt.setp(yticklabels, visible=False)
plt.setp(xticklabels, fontsize=4)
leg = fig.legend(lines, n, loc="lower left")
for t in leg.get_texts():
t.set_fontsize('small') # the legend text fontsize
if filename is not None:
plt.savefig(filename)
[docs]def plot_cmb_powspec(cmb, ecmb, bins, polar=False, cmb_unit="uK_CMB", ref_beam=5, filename=None):
""" Plot the fitted CMB power spectra versus the 2013 lambdaCDM best fit.
"""
if not isinstance(cmb, list):
ecmb = [ecmb]
cmb = [cmb]
if polar:
BBcq= [np.squeeze(icmb.powspec()[1,1,:]) for icmb in cmb]
BBecq = [np.squeeze(icmb.powspec()[1,1,:]) for icmb in ecmb]
EBcq = [np.squeeze(icmb.powspec()[0,1,:]) for icmb in cmb]
EBecq = [np.squeeze(icmb.powspec()[0,1,:]) for icmb in ecmb]
EEcq = [np.squeeze(icmb.powspec()[0,0,:]) for icmb in cmb]
EEecq = [np.squeeze(icmb.powspec()[0,0,:]) for icmb in ecmb]
EE = get_reference_cmb_powspec('EE', cmb_unit=cmb_unit, ref_beam=ref_beam)
BB = get_reference_cmb_powspec('BB', cmb_unit=cmb_unit, ref_beam=ref_beam)
fits = [(EEcq, EEecq, EE, 'EE'),
(BBcq, BBecq, BB, 'BB'),
(EBcq, EBecq, np.zeros((len(EE))), "EB"),
]
else:
TTcq = [np.squeeze(icmb.powspec()[0,0,:]) for icmb in cmb]
TTecq = [np.squeeze(icmb.powspec()[0,0,:]) for icmb in ecmb]
TT = get_reference_cmb_powspec('TT', cmb_unit=cmb_unit, ref_beam=ref_beam)
TTd = [iTTcq-spherelib.bin_cl(TT, bins) for iTTcq in TTcq]
TTde = [iTTecq for iTTecq in TTecq]
TTr = [iTTcq/spherelib.bin_cl(TT, bins) for iTTcq in TTcq]
TTre = [iTTecq/spherelib.bin_cl(TT, bins) for iTTecq in TTecq]
fits = [(TTcq, TTecq, TT, "$D_\ell$"),
(TTr, TTre, np.ones(TT.shape), "$D_\ell / D_\ell^{ref}$"),
(TTd, TTde, np.zeros(TT.shape), "$D_\ell - D_\ell^{ref}$"),
]
plt.figure()
nf = len(fits)
ll = np.squeeze(spherelib.get_lmean(bins))
for (cq,ecq,ref,name),i_f in zip(fits[0:1],range(nf)[0:1]):
#plt.subplot(nf,1,i_f+1)
plt.ylabel(name)
ell = np.arange(len(ref))
if np.sum(ref)!=len(ref):
cq = [ll*(ll+1)*icq/(2*np.pi) for icq in cq]
ecq = [ll*(ll+1)*(iecq)/(2*np.pi) for iecq in ecq]
ref = ell*(ell+1)*ref/(2*np.pi)
for icq, iecq in zip(cq, ecq):
plt.errorbar(ll, icq, yerr=iecq, marker="o", markersize=5)
a0,a1,a2,a3 = plt.axis()
plt.plot(ell, ref,color='k')
plt.axis([ll[0],ll[-1],a2,a3])
#if i_f==0:
# plt.axis([ll[0],100,-0.05,1.1])
#else:
# plt.axis([ll[0],100,-0.1,0.3])
if filename is not None:
plt.savefig(filename)
[docs]def save_calib(filename, ref, model, em):
""" Save calibration result in txt file.
"""
nmap = model.ndet
CCalib = np.zeros((nmap, 4))
ref = np.squeeze(ref)
fit = np.squeeze(model.get_comp_by_name("cmb").mixmat())
sigma = np.squeeze(em.get_comp_by_name("cmb").mixmat())
CCalib[:,0] = ref
CCalib[:,1] = fit
CCalib[:,2] = fit/ref
CCalib[:,3] = sigma/ref
np.savetxt(filename, CCalib)
[docs]def plot_calib(ref, fit, efit, names, filename=None):
""" Plot calibration result.
"""
ref = np.squeeze(ref)
fit = np.squeeze(fit)
sigma = np.squeeze(efit)
sigma[sigma==fit] = 0
nmap = len(names)
plt.figure()
ax = plt.subplot(111)
majorFormatter = plt.FormatStrFormatter('%g')
ax.yaxis.set_major_formatter(majorFormatter)
plt.plot(np.arange(nmap), np.ones((nmap)),color="k")
colors = ["c", "m", "g", "r", "b", "grey", "orange", "purple", "yellow"]
freqs_all = [data.get_freq(name) for name in names]
freqs = sorted(list(set(freqs_all)))
print freqs_all
print freqs
for freq,j in zip(freqs, range(len(freqs))):
idx = np.where(np.array(freqs_all)==freq)
print idx
plt.errorbar(np.arange(nmap)[idx], fit[idx]/ref[idx], yerr=sigma[idx]/ref[idx], color=colors[j], linewidth=2, ls='.', marker='s')
names_show = [name.split("_T")[0] for name in names]
plt.xticks(np.arange(nmap), tuple(names_show),rotation=45)
[xi, xa, yi, ya] = plt.axis()
plt.axis([-0.5, nmap+0.5-1, 0.99, 1.06])
#plt.ylabel("$\hat a_i$", fontsize=18)
plt.grid(True)
if filename is not None:
plt.savefig(filename)