pipe demo cmb debug + mc

test/cmb/clcmb.py

@@ -2,36 +2,48 @@
 
 Coadd cmb and noise map.
 Compute the power spectrum 
+Correct it from the noise bias
 Correct it from the mode coupling matrix
 """
 
+### Gather inputs by sim_id
+#multiplex cross_prod group_by 'cmb'
+sim_id   = seg_input 
+
 import healpy as hp
 import pylab as pl
 import spherelib as sp
 
+### Define some global parameters
+lst_par = ['sim_id']
+lst_tag = lst_par
+
 ### Retrieve some global parameters
 load_param("cmb", globals(), ["lmax", "cmb_unit"])
-load_param("noise", globals(), ["noise_unit", "noise_freq"])
-
-#multiplex cross_prod group_by 'cmb'
-seed   = seg_input["cmb"]
+load_param("noise", globals(), ["noise_unit", "noise_power"])
 
 ### Get inputs
-cmb_fn  = glob_seg("cmb", "map*.fits")[0]
-cmb_map = hp.read_map(cmb_fn)
-
+cmb_fn   = glob_seg("cmb", "map*.fits")[0]
+cmb_map  = hp.read_map(cmb_fn)
 noise_fn = glob_seg("noise", "map*.fits")[0]
 noise    = hp.read_map(noise_fn)
-cv_factor= sp.astro.convfact(freq=noise_freq, from=noise_unit, to=cmb_unit)
+
+### Coadd cmb and noise maps
+cv_factor= sp.astro.convfact(fr=noise_unit, to=cmb_unit)
 noise    = noise*cv_factor
+cmb_map  = cmb_map + noise
+
+### Make a plot
+coadd_fn = get_data_fn("noisy-cmb.png")
+hp.mollview(cmb_map, title="noisy cmb in %s"%cmb_unit)
+pl.savefig(coadd_fn)
 
-### Mask the input map
-cmb_map = cmb_map + noise
 
 ### Compute power spectrum
 cl = hp.anafast (cmb_map, lmax)
+cl = cl - (noise_power*cv_factor*cv_factor) ## debias
 
-### save it to disk
+### save unbiased cl to disk
 cl_fn = get_data_fn("cls.txt")
 pl.savetxt(cl_fn, cl)
 
@@ -39,14 +51,15 @@ pl.savetxt(cl_fn, cl)
 mll_fn = glob_seg('noise', 'mll*')[0]
 load_products (mll_fn , globals(), ["mll"])
 
-## correct cls
+## correct cl from coupling matrix
 cl = pl.linalg.solve(mll , cl)
+cl = cl * sum(mll[0,:])
 
-### save it to disk
+### save mode couplit corrected cl to disk
 cl_fn = get_data_fn("cls-mll.txt")
 pl.savetxt(cl_fn, cl)
 
-### Set output
-seg_output = [seed] 
+### Set output : forward as many childs as sim ids
+seg_output = [sim_id] 
 
 

test/cmb/clplot.py

@@ -7,26 +7,76 @@ import healpy as hp
 import pylab as pl
 import spherelib as sp
 
-load_param('cmb',  globals(), ["input_cl"])
+### Gather all sim_ids
+#multiplex cross_prod group_by '0'
 
-## load cls
-cl_cmb = pl.loadtxt(glob_seg('clcmb', 'cl*.txt')[0])
-cl_mll = pl.loadtxt(glob_seg('clcmb', 'cl*mll*.txt')[0])
+### Define some global parameters
+lst_par = ['deltal', 'lmin_ratio']
+lst_tag = lst_par
 
-## smoothed cls
-dloverl = 0.05
-bin_cl_cmb = sp.smooth_cl(cl_cmb, dloverl=dloverl)
-bin_cl_mll = sp.smooth_cl(cl_mll, dloverl=dloverl)
+### Retrieve some global parameters
+load_param('cmb',   globals(), ["input_cl", "lmax", "cmb_unit"])
+load_param("noise", globals(), ["noise_unit", "noise_power"])
 
-## input cls
+## load cls
 inp_cls = pl.loadtxt(input_cl)[0:lmax+1,0]
+try:
+    cl_n = pl.loadtxt(glob_seg('clnoise', 'cl*.txt')[0])
+except:
+    cl_n = pl.ones((lmax+1))*noise_power
+
+## Convert noise in cmb unit
+cv_factor= sp.astro.convfact(fr=noise_unit, to=cmb_unit)
+cl_n = cl_n*cv_factor**2
+
+## Get mean cl values and error bars
+cl_unbs  = glob_seg('clcmb', 'cls.txt')
+cl_mlls  = glob_seg('clcmb', 'cl*mll*.txt')
+
+## Binning
+deltal = 10
+bins   = sp.uniformbin(0, lmax, deltal)
+lmean  = pl.squeeze(sp.get_lmean(bins))
+nbins  = len(lmean)
+nsims    = len(cl_unbs)
+cl_unb   = pl.zeros((nbins, nsims))
+cl_mll   = pl.zeros((nbins, nsims))
+for sim_id in range(nsims):
+    cl_unb[:,sim_id] = sp.bin_cl(pl.loadtxt(cl_unbs[sim_id]), bins)
+    cl_mll[:,sim_id] = sp.bin_cl(pl.loadtxt(cl_mlls[sim_id]), bins)
+er_unb = pl.std(cl_unb, axis=1)
+er_mll = pl.std(cl_mll, axis=1)
+cl_unb = pl.mean(cl_unb, axis=1)
+cl_mll = pl.mean(cl_mll, axis=1)
+
 
 ## make a plot
+smooth_win = 0.05
 ell = pl.arange(lmax+1)
 pl.figure()
-pl.plot(ell, ell*(ell+1)*bin_cl_cmb/(2*pl.pi))
-pl.plot(ell, ell*(ell+1)*bin_cl_mll/(2*pl.pi))
-pl.plot(ell, ell*(ell+1)*inp_cls/(2*pl.pi), color="black")
-pl.legend([ "empirical", "mll", "input"], loc="upper right")
+pl.errorbar(lmean, lmean*(lmean+1)*cl_unb/(2*pl.pi), yerr=lmean*(lmean+1)*er_unb /(2*pl.pi),color="b")
+pl.errorbar(lmean, lmean*(lmean+1)*cl_mll/(2*pl.pi), yerr=lmean*(lmean+1)*er_mll /(2*pl.pi),color="c" )
+pl.plot(ell, ell*(ell+1)*cl_n /(2*pl.pi),color="r")
+pl.plot(ell, ell*(ell+1)*inp_cls /(2*pl.pi),color="k")
+pl.ylabel("$\ell (\ell+1) c_\ell / 2 \pi$", fontsize=18)
+pl.xlabel("$\ell$", fontsize=18)
+pl.title("Power spectra in %s$^2$"%cmb_unit)
+
+line1 = pl.Line2D(lmean, lmean, linestyle='-', color='b',linewidth=1)
+line2 = pl.Line2D(lmean, lmean, linestyle='-',color='c', linewidth=1)
+line3 = pl.Line2D(lmean, lmean, linestyle='-',color='r', linewidth=1)
+line4 = pl.Line2D(lmean, lmean, linestyle='-',color='k', linewidth=1)
+pl.legend((line1, line2,line3,line4),["unbiased", "mll", "noise", "input"], loc="upper center")
+
 pl.savefig(get_data_fn("cls.png"))
 pl.savefig(get_data_fn("cls.pdf"))
+
+# lmin_ratio = 200
+# pl.figure()
+# pl.plot(pl.arange(lmin_ratio,lmax+1), sp.smooth_cl(cl_unb[lmin_ratio:]/cl_mll[lmin_ratio:],dloverl=smooth_win))
+# pl.legend([ "unbiased/mll"], loc="upper center")
+# pl.ylabel("$c_\ell^{unb} /c_\ell^{mll} $", fontsize=18)
+# pl.xlabel("$\ell$", fontsize=18)
+# pl.title("Power spectra ratio")
+# pl.savefig(get_data_fn("ratio.png"))
+

test/cmb/cmb.py

@@ -7,35 +7,31 @@ import healpy as hp
 import pylab as pl
 
 ### Define some global parameters
-lst_par = ['lmax', 'nside', 'cmb_unit', 'seed', 'input_cl']
-lst_tag = lst_par[:-2]
-
-nside = seg_input.values()[0]
-seed  = seg_input.values()[1] ## actually not use by synfast
-lmax  = 2*nside
-cmb_unit = "uK_CMB"
+lst_par = ['lmax', 'nside', 'cmb_unit', 'sim_id', 'input_cl']
+lst_tag = lst_par[:-1]
 
+nside  = seg_input.values()[0][0] ## pushed from main
+sim_id = seg_input.values()[0][1] ## pushed from main
+lmax   = 2*nside
 
 ### Generate a cmb map
 input_cl = "lambda_best_fit.txt"
-cmb_cl  = pl.loadtxt(input_cl)[0:lmax+1,0]
-cmb_map = hp.synfast(cmb_cl, nside, lmax=lmax)
-
+cmb_unit = "uK_CMB"
+cmb_cl   = pl.loadtxt(input_cl)[0:lmax+1,0]     ## load cl
+cmb_map  = hp.synfast(cmb_cl, nside, lmax=lmax) ## make a map
 
-### Save it to disk
+### Save to disk
 cmb_cl_fn  = get_data_fn ('cls_cmb.txt')
+pl.savetxt  (cmb_cl_fn , cmb_cl)
 cmb_map_fn = get_data_fn ('map_cmb.fits')
 hp.write_map(cmb_map_fn, cmb_map)
-pl.savetxt  (cmb_cl_fn , cmb_cl)
-
 
 ### Make a plot
 cmb_map_fig = cmb_map_fn.replace('.fits', '.png')
-hp.mollview(cmb_map)
+hp.mollview(cmb_map, title="cmb in %s"%cmb_unit)
 pl.savefig (cmb_map_fig)
 
-
-### Set output
-seg_output = [seed] 
+### Set output : forward as many childs as sim ids
+seg_output = [sim_id] 
 
 

test/cmb/main.py

@@ -32,12 +32,11 @@ import os.path as path
 ###
 
 nside = 512
-seeds = [1]
+sim_ids = [2,4,5]
 
 pipe_dot = """ 
-cmb->clcmb->clplot;
-noise->clcmb;
-noise->clnoise->clplot;
+noise->clcmb->clplot;
+cmb->clcmb;
 """
 
 ### code_dir
@@ -98,9 +97,13 @@ def main(print_info=print_info):
     
     ## Build pipeline instance
     P = Pipeline(pipe_dot, code_dir=code_dir, prefix=prefix, matplotlib=True, matplotlib_interactive=True)
-    for seed in seeds:
-        P.push(cmb=[(nside,seed)])
-        P.push(mask=[(nside,seed)])
+    cmbin = []
+    noisein = []
+    for sim_id in sim_ids:
+        cmbin.append((nside, sim_id))
+        noisein.append((nside, sim_id))
+    P.push(cmb=cmbin)
+    P.push(noise=noisein)
     
     ## Interactive mode
     if options.debug:

test/cmb/noise.py

 """ noise.py
 
 Generate a noise map from hitcount map.
-Compute coupling matrix
+Compute coupling matrix. 
+Compute noise power for debiasing
 """
 
 import healpy as hp
 import pylab as pl
+import numpy as npy
+import spherelib as sp
 
 ### Define some global parameters
-lst_par = ['lmax', 'nside', 'noise_unit', 'seed', 'input_hitcount', 'noise_level']
-lst_tag = lst_par[:-2]
+lst_par = ['sigma', 'input_hitcount', 'sim_id', 'lmax', 'nside', 'noise_unit', 'noise_power']
+lst_tag = lst_par[:-4]
 
-nside = seg_input.values()[0]
-seed  = seg_input.values()[1] ## actually not use by synfast
-lmax  = 2*nside
-noise_unit = "K_RJ"
+nside  = seg_input.values()[0][0] ## pushed from main
+sim_id = seg_input.values()[0][1] ## pushed from main
+lmax   = 2*nside
 
+### Noise parameters
+input_hitcount = "hitmap-wmap5-W.fits"
+noise_unit     = "mK_CMB"
+sigma          = 6.538
 
 ### Generate a noise map
-input_hitcount = ""
 hitcount = hp.read_map(input_hitcount)
 hitcount = hp.ud_grade(hitcount, nside)
 hitnside = pl.sqrt(len(hitcount)/12)
 hitcount = hitcount * (hitnside/nside)**2
-noise_map= pl.random.standard_normal(12*nside*nside)*noise_level / pl.sqrt(hitcount)
+noise_map= npy.random.standard_normal(12*nside*nside)*sigma / pl.sqrt(hitcount)
 
 ### Save it to disk
 noise_fn = get_data_fn("map_noise.fits")
@@ -31,13 +36,20 @@ hp.write_map(noise_fn, noise_map)
 
 ### Compute mll'
 mll    = sp.mask2mll (hitcount, lmax)
+mll    = mll[0:lmax+1,0:lmax+1]
 mll_fn = get_data_fn("mll.dat")
 save_products (mll_fn, globals(), ["mll"])
 
+### Compute noise power
+noise_power = 4*pl.pi/len(hitcount)**2 * sigma**2 * pl.sum( 1/hitcount) 
+  
 ### Make a plot
 noise_fig = noise_fn.replace('.fits', '.png')
-hp.mollview(noise_map)
+hp.mollview(noise_map, title="noise in %s"%noise_unit)
+pl.savefig (noise_fig)
+noise_fig = get_data_fn(input_hitcount).replace('.fits', '.png')
+hp.mollview(hitcount, title="hitcount")
 pl.savefig (noise_fig)
 
-### Set output
-seg_output = [seed]
+### Set output : forward as many childs as sim ids
+seg_output = [sim_id]