#multiplex cross_prod group_by 'map2alm[0],powspec' where 'powspec in ["cmb", "sync", "dust"] and map2alm[0]=="full"'
import healpy
import spherelib
from pipetools import *
from smicatools import *
from maptools import *
import numpy as np
import matplotlib.pyplot as plt
import planckdata as data
def plot_snr(snr, beam, comp):
[docs] """
"""
plt.figure()
plt.title(comp)
plt.plot(snr.T)
plt.plot(beam )
plt.legend(["snrE", "snrB", "beam"])
plt.savefig(get_data_fn("snr.png"))
def plot_outcmbcl(Fl, hRl, Rl, acmb, N, name=""):
[docs] lmax = hRl.shape[2]-1
nmap = hRl.shape[0]
nfield = Fl.shape[0]
outcl = np.zeros((nfield, lmax+1))
outnl = np.zeros((nfield, lmax+1))
plt.figure()
ells = np.arange(lmax+1)
for j in range(1):#nfield):
for ell in range(lmax+1):
outcl[j, ell] = np.dot(Fl[j,:,ell].reshape(1,nmap), np.dot(hRl[:,:,ell], Fl[j,:,ell].reshape(nmap,1)))
outnl[j, ell] = np.sum(Fl[j,:,ell]**2*N[:,ell])
healpy.write_cl("cl_cmb_proj_%d.fits"%j, outcl[j,:])
plt.semilogy(ells**2*outcl[j,:], label="total", color="brown")
plt.semilogy(ells**2*outnl[j,:], label="noise", color="r")
plt.semilogy(ells**2*spherelib.smooth_cl((outcl[j,:]-outnl[j,:])), label="total-noise", color="c")
plt.semilogy(ells**2*np.squeeze(Rl[j*nmap/2+0,j*nmap/2+0]/acmb[j*nmap/2+0,j]**2), label="cmb fit", color="b")
if nfield==1:
ref = get_reference_cmb_powspec("TT")
else:
ref = get_reference_cmb_powspec("EE")
ref = np.hstack((ref, np.nan*np.ones((lmax+1-len(ref)))))
plt.semilogy(ells**2*ref[0:lmax+1], label="best fit", color="k")
plt.legend(loc="lower center")
if nfield==1:
plt.axis([0,2500,1e-1,7e4])
else:
plt.axis([0,2000,1e-1,5e3])
plt.savefig(get_data_fn("cl_cmb_proj_%s.png"%name))
bins = get_bins(0,lmax)
ll = np.squeeze(spherelib.get_lmean(bins))
plt.figure()
plt.plot(ll, ll**2*spherelib.bin_cl((outcl[j,:]-outnl[j,:])-ref, bins), color="c")
plt.plot(ll, ll**2*spherelib.bin_cl(np.squeeze(Rl[j*nmap/2+0,j*nmap/2+0]/acmb[j*nmap/2+0,j]**2)-ref, bins), color="b")
plt.plot(ll, np.ones((len(ll))), color="k")
plt.legend(loc="upper center")
if nfield==1:
plt.axis([0,2500,-100,150])
else:
plt.axis([0,2000,-30, 30])
plt.savefig(get_data_fn("deltacl_cmb_proj_%s.png"%name))
def alm2map(alms, nfield, nside=2048):
[docs] """ Run healpy.alm2map with correct input format, and return a list of map
T: alm array
E,B: list of alm
"""
nalm = len(alms[0])
alms = alms[0] if nfield==1 else [np.zeros((nalm),dtype=np.complex128)]+alms
maps = healpy.alm2map(alms, nside)
maps = [maps] if nfield==1 else maps[1:]
return maps
def get_minmax(comp, polar, cmb_unit, lowres=False, residual=False):
[docs] """
"""
if comp=='cmb':
if polar:
if lowres:
vmin,vmax = -2.5,2.5
else:
vmin,vmax = -50,50
else:
if lowres:
if residual:
vmin,vmax = -10, 10
else:
vmin,vmax = -200, 200
else:
vmin,vmax = -300, 300
cv = astro.convfact(freq=1e9, fr='uK_CMB', to=cmb_unit)
vmin *= cv
vmax *= cv
norm = None
else:
vmin,vmax,norm = None, None, "norm"
return vmin,vmax,norm
def diff_with_input(almFs, fields, bl_lowres, nside, vmin, vmax, cmb_unit):
[docs] """
Compute component map
"""
dset,comp = get_input()
load_param("map2alm", globals(), ["lmax", "unit_to"])
load_param("mixmat", globals(), ["names", "polar", "cmb_unit"])
load_param("powspec", globals(), ["bins"])
units = dict({"cmb":cmb_unit, "dust":unit_to, "sync":unit_to})
fields = ["E","B"] if polar else ["T"]
fieldsrs = ["Q","U"] if polar else ["I"] #real space
nfield = len(fields)
nmap = len(names)
# compute ilc filter
load_products (glob_parent("model_full.pkl")[0], globals(), ['model'])
load_products (glob_seg("mixmat", "covmat_unbinned_full.pkl")[0], globals(), ['lstats'])
N = np.loadtxt(glob_seg("mixmat", "noise_unbinned.txt")[0])
if dset!="full":
N/=2.
Acmb = model.get_comp_by_name("cmb").mixmat()
ibl = load_ibl(get_lstibl(names), lmax)
if comp=="cmb":
lcst = 10
lilc = 1500
Rl = unbin_R(model.covariance(), lstats, bins, lcst=lcst, lilc=lilc)
Fh = get_ilc_filter(Rl, Acmb)
plot_noise_decomp(Fh, N, names, fields, filename=get_data_fn("noise.png"))
np.savetxt (get_data_fn("raw_filter_full.txt"), Fh.reshape(nfield*nmap, lmax+1))
Rlcmb = unbin_R(model.get_comp_by_name("cmb").covariance(), lstats, bins)
plot_outcmbcl(Fh, lstats, Rlcmb, Acmb, N, name="full-yrhd")
# save filter for monte carlo
Fs = Fh.copy()
for m in range(nmap):
Fh[:,m,:]*= ibl[m,:] # include beam correction to filters
Fs[:,m,:] = Fh[:,m,:]*np.sum(Acmb[m,:]) # filters to be applied to maps in uK_CMB
np.savetxt (get_data_fn("mc_filter_full.txt"), Fs.reshape(nfield*nmap, lmax+1))
plot_filter(Fh*1e-6, names, fields, filename=get_data_fn("filters.png"))
elif comp in ["sync", "dust"] and polar:
lcst = 2
Rl = unbin_R(model.covariance(), lstats, bins, lcst=lcst)
A = model.get_comp_by_name("gal").mixmat()
d = A.shape[1]
P = model.get_comp_by_name("gal").powspec()
Pl = unbin_R(P, np.zeros((d,d,lmax+1)), bins, lcst=lcst)
dreco = dict({"sync":((2,3),(0,7),60,256), "dust":((0,1),(6,13),10,1024)})
units["sync"]+= "_%s"%(names[dreco["sync"][1][0]].split("_")[0])
units["dust"]+= "_%s"%(names[dreco["dust"][1][0]].split("_")[0])
Fh, snr = get_fg_filter(Rl, dreco[comp][0], dreco[comp][1], A, Pl)
beam = np.squeeze(spherelib.gaussbeam(dreco[comp][2], lmax))
Fs = Fh.copy()
for m in range(nmap):
Fh[:,m,:] /= snr
Fh[np.isinf(Fh)] = 0
Fh[np.isnan(Fh)] = 0
Fh[:,m,:] *= ibl[m,:] # include beam correction to filters
Fh[:,m,:] *= beam
Fs[:,m,:] = Fh[:,m,:]*np.sum(Acmb[m,:]) # filters to be applied to maps in uK_CMB
np.savetxt (get_data_fn("raw_filter_full.txt"), Fh.reshape(nfield*nmap, lmax+1))
plot_filter(Fh, names, fields, filename=get_data_fn("filters.png"))
plot_snr(snr, beam, comp)
# compute component map
option = 'join'
nalm = (lmax+1)*(lmax+2)/2
salm = np.zeros((nmap,nalm), dtype=np.complex128)
for name, m in zip(names, range(nmap)):
if polar and comp=='cmb':
inalm = glob_parent("alm_masked_%s_%s.fits"%(dset,name))[0]
else:
inalm = glob_parent("alm_full_sky_%s_%s.fits"%(dset,name))[0]
salm[m,:] = healpy.read_alm(inalm)
combalm = [np.zeros((nalm), dtype=np.complex128) for i in range(nfield)]
for m,name in zip(range(nmap),names):
for f,field in zip(range(nfield),fields):
if (option=='indep' and name[-1]==field) or option=='join':
combalm[f] += healpy.almxfl(salm[m,:]/np.sum(Acmb[m,:]), np.squeeze(Fs[f,m,:]))
del salm
write_alm(combalm, fields, "alm_%s_%s"%(dset,comp))
maps = alm2map(combalm, nfield)
vmin,vmax,norm = get_minmax(comp, polar, cmb_unit)
for field,mapF in zip(fieldsrs, maps):
healpy.write_map(get_data_fn("map_%s_%s_%s.fits"%(dset,comp,field)), mapF)
plot_map_in_pipe(mapF, "map", "%s-%s-%s"%(dset,comp,field), min=vmin, max=vmax, norm=norm)
if polar:
write_map_QU(get_data_fn("%s_smica_harmonic_%s_%s.fits"%(os.getenv("RELEASE"),comp,dset)), maps[0], maps[1], units=units[comp])
else:
write_map_I(get_data_fn("%s_smica_harmonic_%s_%s.fits"%(os.getenv("RELEASE"),comp,dset)), maps[0], units=units[comp])
# show low resolution map
bl_lowres = np.squeeze(spherelib.gaussbeam (80, lmax))
nside_lowres = 512
vmin,vmax,norm = get_minmax(comp, polar, cmb_unit, lowres=True)
for i in range(nfield):
combalm[i] = healpy.almxfl(combalm[i], bl_lowres)
diff_with_input(combalm, fields, bl_lowres, nside_lowres, vmin, vmax, cmb_unit)
maps = alm2map(combalm, nfield)
for field,mapF in zip(fieldsrs, maps):
healpy.write_map(get_data_fn("map_lowres_%d_%s_%s_%s.fits"%(nside_lowres,dset,comp,field)), mapF)
plot_map_in_pipe(mapF, "map_lowres", "%s-%s-%s"%(dset,comp,field), min=vmin, max=vmax, norm=norm)
set_output(get_input())
globals().update(locals())
save_param([])
expose(["option"])
if __name__=="__builtin__":
main()