demo cmb ok

README.org

@@ -142,6 +142,15 @@ P = Pipeline(pipedot, codedir=, prefix=)
 - codedir is the path where the segment scripts can be found
 - prefix  is the path to the data repository (see below Hierarchical data storage)
 
+It is possible to output the graphviz representation of the pipeline. 
+First, save the graph string into a .dot file with the pipelet function: 
+
+P.to_dot('pipeline.dot')
+
+Then, convert it to an image file with the dot command: 
+
+dot -Tpng -o pipeline.png pipeline.dot
+
 *** Dependencies between segments
 
 The modification of the code of one segment will trigger its

test/cmb/clcmb.py

@@ -6,54 +6,65 @@ 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
+import numpy as npy
 
 ### Define some global parameters
 lst_par = ['sim_id']
 lst_tag = lst_par
 
 ### Retrieve some global parameters
+sim_id   = seg_input['cmb']
+
 load_param("cmb", globals(), ["lmax", "cmb_unit"])
-load_param("noise", globals(), ["noise_unit", "noise_power"])
+load_param("noise", globals(), ["noise_unit", "noise_power", "nside", "sigma"])
 
 ### Get inputs
 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)
+hitcount_fn= glob_seg("noise", "hit*.fits")[0]
+hitcount = hp.read_map(hitcount_fn)
 
 ### 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
+nalm     = (lmax+1)*(lmax+2)/2
+maps     = []
+alms     = pl.zeros((2,nalm), dtype=complex)
+for m in range(2):
+    noise_map= npy.random.randn(12*nside*nside)*sigma / pl.sqrt(hitcount)
+    maps.append( cv_factor*noise_map + cmb_map)
+for m in range(2):
+    alms[m,:]= hp.map2alm(maps[m], lmax) 
 
 ### Make a plot
 coadd_fn = get_data_fn("noisy-cmb.png")
-hp.mollview(cmb_map, title="noisy cmb in %s"%cmb_unit)
+hp.mollview(maps[0], title="noisy cmb in %s"%cmb_unit)
 pl.savefig(coadd_fn)
 
-
-### Compute power spectrum
-cl = hp.anafast (cmb_map, lmax)
-cl = cl - (noise_power*cv_factor*cv_factor) ## debias
+### Compute cross power spectrum
+(cov,jk) = sp.alm2cov(alms)
+cl = cov[1,0,:]
 
 ### save unbiased cl to disk
 cl_fn = get_data_fn("cls.txt")
 pl.savetxt(cl_fn, cl)
 
+### Apply inverse noise weighting
+for m in range(2):
+    alms[m,:]= hp.map2alm(hitcount*maps[m], lmax) 
+
+### Compute cross power spectrum
+(cov,jk) = sp.alm2cov (alms)
+cl = cov[0,1,:]
+
 ## load mll'
 mll_fn = glob_seg('noise', 'mll*')[0]
 load_products (mll_fn , globals(), ["mll"])
 
 ## correct cl from coupling matrix
 cl = pl.linalg.solve(mll , cl)
-cl = cl * sum(mll[0,:])
 
 ### save mode couplit corrected cl to disk
 cl_fn = get_data_fn("cls-mll.txt")

test/cmb/clplot.py

@@ -34,8 +34,8 @@ cl_unbs  = glob_seg('clcmb', 'cls.txt')
 cl_mlls  = glob_seg('clcmb', 'cl*mll*.txt')
 
 ## Binning
-deltal = 10
-bins   = sp.uniformbin(0, lmax, deltal)
+deltal = 2
+bins   = sp.wg2bin(0, lmax)
 lmean  = pl.squeeze(sp.get_lmean(bins))
 nbins  = len(lmean)
 nsims    = len(cl_unbs)
@@ -49,34 +49,35 @@ 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.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)
+for fn_plot in [pl.plot, pl.semilogx]:
+    ell = pl.arange(lmax+1)
+    pl.figure()
+    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" )
+    fn_plot(ell, ell*(ell+1)*cl_n /(2*pl.pi),color="r")
+    fn_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),["pseudo $c_\ell$ "," pseudo $m_{\ell \ell}^{-1} c_\ell$", "noise", "input"], loc="upper center")
+    [xa, xb, ya, yb] = pl.axis()
+    pl.axis([0,lmax,ya,yb])
+    pl.savefig(get_data_fn("cls-%s.png"%fn_plot.func_name))
+    pl.savefig(get_data_fn("cls-%s.pdf"%fn_plot.func_name))
 
-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"))
+pl.figure()
+pl.plot(lmean, er_unb/er_mll)
+pl.plot(lmean, pl.ones(pl.shape(lmean)), color="black")
+pl.legend([ "pseudo $c_\ell$ / pseudo $m_{\ell \ell}^{-1} c_\ell$"], loc="upper center")
+pl.ylabel("$\sigma_\ell^{unb} / \sigma_\ell^{mll} $", fontsize=18)
+pl.xlabel("$\ell$", fontsize=18)
+pl.title("Error bars ratio")
+pl.savefig(get_data_fn("ratio.png"))
 

test/cmb/main.py

@@ -32,7 +32,7 @@ import os.path as path
 ###
 
 nside = 512
-sim_ids = [2,4,5]
+sim_ids = [1,2]
 
 pipe_dot = """ 
 noise->clcmb->clplot;
@@ -98,12 +98,10 @@ def main(print_info=print_info):
     ## Build pipeline instance
     P = Pipeline(pipe_dot, code_dir=code_dir, prefix=prefix, matplotlib=True, matplotlib_interactive=True)
     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)
+    P.push(noise=[nside])
     
     ## Interactive mode
     if options.debug:

test/cmb/noise.py

@@ -14,8 +14,7 @@ import spherelib as sp
 lst_par = ['sigma', 'input_hitcount', 'sim_id', 'lmax', 'nside', 'noise_unit', 'noise_power']
 lst_tag = lst_par[:-4]
 
-nside  = seg_input.values()[0][0] ## pushed from main
-sim_id = seg_input.values()[0][1] ## pushed from main
+nside  = seg_input.values()[0] ## pushed from main
 lmax   = 2*nside
 
 ### Noise parameters
@@ -28,14 +27,13 @@ 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= npy.random.standard_normal(12*nside*nside)*sigma / pl.sqrt(hitcount)
 
 ### Save it to disk
-noise_fn = get_data_fn("map_noise.fits")
-hp.write_map(noise_fn, noise_map)
+noise_fn = get_data_fn("hitmap_noise.fits")
+hp.write_map(noise_fn, hitcount)
 
 ### Compute mll'
-mll    = sp.mask2mll (hitcount, lmax)
+mll    = sp.mask2mll (hitcount, 2*lmax)
 mll    = mll[0:lmax+1,0:lmax+1]
 mll_fn = get_data_fn("mll.dat")
 save_products (mll_fn, globals(), ["mll"])
@@ -44,12 +42,9 @@ save_products (mll_fn, globals(), ["mll"])
 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, 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 : forward as many childs as sim ids
-seg_output = [sim_id]
+seg_output = [input_hitcount]