import healpy
import spherelib
import numpy as np
import os.path as op
import smica
from pipetools import *
from smicatools import *
import spherelib
import sys
[docs]def main():
"""
Fit the power spectra from 2nd order statistics
"""
# set input channels
load_param("mixmat", globals(), ["cmb_unit", "names", "polar", "dim"])
load_param("map2alm", globals(), ["unit_to", "ref_beam"])
nmap = len(names)
# compute statistics from alm
load_products(glob_parent("covmat_unbinned_full.pkl")[0], globals(), ['lstats', 'lnmodes'])
# compute noise from jackknives
Nl = np.loadtxt(glob_parent("noise_unbinned.txt")[0])
# binning
bins = get_bins(5, 2000)
lmin_fit = bins[0,0]
lmax_fit = bins[-1,1]
lmean = np.squeeze(spherelib.get_lmean(bins))
stats, nmodes = spherelib.bin_stat (lstats, bins, nmode=np.squeeze(lnmodes))
nmodes = np.squeeze(nmodes)
Q = len(nmodes) #number of bins
save_products(get_data_fn("covmat_binned_full.pkl"), locals(),
['stats', 'nmodes', 'bins', 'lmean'])
write_stats_in_txt(lstats, lnmodes, suffix="_unbinned_full")
write_stats_in_txt(stats, nmodes, suffix="_binned_full")
N = np.zeros((nmap, Q))
for m in range(nmap):
N[m,:] = spherelib.bin_cl(Nl[m,:], bins)
np.savetxt(get_data_fn("noise_binned.txt"), N)
plot_cov(stats, N=N, lmean=lmean, filename=get_data_fn("cov_binned.png"))
plot_ev_decomp(stats, N, lmean=lmean, filename=get_data_fn("ev_binned.png"))
# build the model
if 0:
(of,oe) = logged_subprocess (['lspipe', 'calibdx11dr2@planck:plotc/data/*/calib-mean.txt'])
cfile = file(of).read().strip()
else:
cfile = None
cmb_unit = 'uK_CMB'
cmb = get_cmb(names, unit_to, bins, polar, calib=cfile, ref_beam=ref_beam, cmb_unit=cmb_unit)
gal,eA = get_gal_from(Q, filename=glob_parent("model_full.pkl")[0])
noise_fix = False if polar else True
noise = get_noise(N, fixed=noise_fix)
model = smica.Model(complist=[cmb, gal, noise])
# fit model to covariance
model = fit_model_to_cov(model, stats, nmodes, maxiter=150, noise_fix=noise_fix, logger=logger)
final_mismatch = int(np.real(model.mismatch(stats, nmodes, exact=True)))
gal.fix_powspec('all') # just to save time
em = model.error_model(nmodes)
gal.fix_powspec('null') # just to save time
save_products(get_data_fn("model_full.pkl"), locals(), ['model', 'em'])
em.get_comp_by_name("gal").set_mixmat(eA) # replace with value from mixmat step
write_model_in_txt(model, bins, em=em)
plot_mismatch(model, stats, nmodes, lmean, filename=get_data_fn("mismatch.png"))
model.plot_power(get_data_fn("power.png"), xaxis=lmean)
plot_relative_diff(model, stats, nmodes, names, lmean, filename=get_data_fn("relativediff.pdf"))
plot_decomp(model, stats, nmodes, lmean, names, filename=get_data_fn("decomp.pdf"))
plot_mismatch_pair(model, stats, nmodes, names, lmean, filename=get_data_fn("mismatch_pair.pdf"))
plot_noise(model, em, N, lmean, names, filename=[get_data_fn("noise.png"), get_data_fn("noise_ratio.png")])
plot_cmb_powspec(cmb, em.get_comp_by_name("cmb"), bins, polar=polar, cmb_unit=cmb_unit, ref_beam=ref_beam, filename=get_data_fn("cmb_powspec.png"))
try:
Kss,Knn,Ksn = plot_mismatch_KLD(model, stats, nmodes, names, lmean, filename=[get_data_fn("KLD0.png"), get_data_fn("KLD1.png")])
Kss_mismatch = int(np.sum(Kss))
except:
logger.info(sys.exc_info()[0])
Kss_mismatch = None
if dim==2 and polar:
set_output(["cmb", "sync", "dust"])
else:
set_output([c.name for c in model._complist])
save_param(["bins", "cmb_unit", "names"])
rg_fit = "%d-%d"%(lmin_fit, lmax_fit)
globals().update(locals())
expose(["final_mismatch", "Kss_mismatch", "dim", "rg_fit"])
if __name__=="__builtin__":
main()