adding codes

Author: Serena Giardiello

project/SO/pISO/python/dust/global_dust.dict

@@ -10,6 +10,16 @@ kspace_tf_path  = release_dir + "pspipe/kspace_tf/dust_binning_50/"
 sim_spec_dir = release_dir + "pspipe/sim_spectra/dust_binning_50/"
 source_catalog = release_dir + "pspipe/source_catalog/cat_skn_multifreq_20220526_nightonly.txt"
 
+mcm_dir = "/scratch/c.spxsg6/iso/dust/mcms"
+alms_dir = "/scratch/c.spxsg6/iso/dust/alms"
+spec_dir = "/scratch/c.spxsg6/iso/dust/spectra"
+best_fits_dir = "/scratch/c.spxsg6/iso/dust/best_fits"
+plots_base_dir = "/scratch/c.spxsg6/iso/dust/plots"
+ps_model_dir = "/scratch/c.spxsg6/iso/dust/noise_model"
+sq_win_alms_dir = "/scratch/c.spxsg6/iso/dust/sq_win_alms"
+cov_dir = "/scratch/c.spxsg6/iso/dust/covariances"
+montecarlo_dir = "/scratch/c.spxsg6/iso/dust/montecarlo"
+
 deconvolve_pixwin = True
 pixwin_Planck = {"pix": 'HEALPIX', "nside": 2048}
 
@@ -31,8 +41,8 @@ n_splits_Planck = 2
 # get from running pipeline_dust.yml
 # then extract the patch given by window_so_i5_f220_baseline.fits
 # enmap.extract(existing_planck_maps, window.shape, window.wcs) 
-maps_Planck_f143 = ["planck_projected/npipe6v20%s_f143_map_srcfree_deep56_patch.fits" % split for split in ["A", "B"]] 
-maps_Planck_f353 = ["planck_projected/npipe6v20%s_f353_map_srcfree_deep56_patch.fits" % split for split in ["A", "B"]]
+maps_Planck_f143 = ["planck_projected/npipe6v20%s_f143_map_deep56_patch_srcfree.fits" % split for split in ["A", "B"]] 
+maps_Planck_f353 = ["planck_projected/npipe6v20%s_f353_map_deep56_patch_srcfree.fits" % split for split in ["A", "B"]]
 
 cal_Planck_f143 = 1.
 cal_Planck_f353 = 1.
@@ -64,8 +74,9 @@ cross_link_threshold = 0.97
 n_med_ivar = 3
 
 mask_dir = release_dir + "pspipe/masks/"
-ps_mask_Planck_f143 = mask_dir + "source_mask_15mJy_and_dust_rad12.fits"
-ps_mask_Planck_f353 = mask_dir + "source_mask_15mJy_and_dust_rad12.fits"
+# also these have to be projected over the patch
+ps_mask_Planck_f143 = mask_dir + "source_mask_15mJy_and_dust_rad12_deep56.fits"
+ps_mask_Planck_f353 = mask_dir + "source_mask_15mJy_and_dust_rad12_deep56.fits"
 
 window_dir = release_dir + "windows/"
 

project/SO/pISO/python/get_covariance_blocks.py

@@ -12,12 +12,12 @@
 
 log = log.get_logger(**d)
 
-mcms_dir = "mcms"
-spectra_dir = "spectra"
-noise_dir = "noise_model"
-cov_dir = "covariances"
-bestfit_dir = "best_fits"
-sq_win_alms_dir = "sq_win_alms"
+mcms_dir = d["mcm_dir"]
+spectra_dir = d["spec_dir"]
+noise_dir = d["ps_model_dir"]
+cov_dir = d["cov_dir"]
+bestfit_dir = d["best_fits_dir"]
+sq_win_alms_dir = d["sq_win_alms_dir"]
 
 pspy_utils.create_directory(cov_dir)
 surveys = d["surveys"]

project/SO/pISO/python/get_noise_model.py

@@ -0,0 +1,117 @@
+#TO BE TESTED MORE
+
+import matplotlib
+
+matplotlib.use("Agg")
+import sys
+
+import numpy as np
+import pylab as plt
+import scipy.interpolate
+from pspipe_utils import log, misc
+from pspy import pspy_utils, so_dict, so_spectra
+
+
+def interpolate_dict(lb, cb, lth, spectra, force_positive=True, l_inf_lmin_equal_lmin=True, discard_cross=True):
+       cl_dict = {}
+       for spec in spectra:
+            cl = scipy.interpolate.interp1d(lb, cb[spec], fill_value = "extrapolate")
+            cl_dict[spec] = cl(lth)
+            if l_inf_lmin_equal_lmin:
+                id = np.where(lth <= np.min(lb))
+                cl_dict[spec][id]= cb[spec][0]
+            if force_positive:
+                cl_dict[spec] = np.abs(cl_dict[spec])
+            if discard_cross:
+                if spec not in ["TT", "EE", "BB"]:
+                    cl_dict[spec] = np.zeros(len(lth))
+       return cl_dict
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+spectra_dir = d['spec_dir']
+ps_model_dir = d['ps_model_dir'] 
+plot_dir = d['plot_base_dir'] + "/noise_model/"
+
+pspy_utils.create_directory(ps_model_dir)
+pspy_utils.create_directory(plot_dir)
+
+spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+surveys = d["surveys"]
+lmax = d["lmax"]
+type = d["type"]
+binning_file = d["binning_file"]
+
+lth = np.arange(2, lmax+2)
+
+for sv in surveys:
+    arrays = d[f"arrays_{sv}"]
+    for id_ar1, ar1 in enumerate(arrays):
+        for id_ar2, ar2 in enumerate(arrays):
+            if id_ar1 > id_ar2: continue
+
+            log.info(f"Computing noise for '{sv}' survey and '{ar1}x{ar2}' arrays")
+
+            l, bl_ar1 = misc.read_beams(d[f"beam_T_{sv}_{ar1}"], d[f"beam_pol_{sv}_{ar1}"])
+            l, bl_ar2 = misc.read_beams(d[f"beam_T_{sv}_{ar2}"], d[f"beam_pol_{sv}_{ar2}"])
+
+            lb, nbs_ar1xar1 = so_spectra.read_ps(f"{spectra_dir}/{type}_{sv}_{ar1}x{sv}_{ar1}_noise.dat", spectra=spectra)
+            lb, nbs_ar1xar2 = so_spectra.read_ps(f"{spectra_dir}/{type}_{sv}_{ar1}x{sv}_{ar2}_noise.dat", spectra=spectra)
+            lb, nbs_ar2xar2 = so_spectra.read_ps(f"{spectra_dir}/{type}_{sv}_{ar2}x{sv}_{ar2}_noise.dat", spectra=spectra)
+        
+            bb_ar1, bb_ar2  = {}, {}
+            for field in ["T", "E", "B"]:
+                lb, bb_ar1[field] = pspy_utils.naive_binning(l, bl_ar1[field], binning_file, lmax)
+                lb, bb_ar2[field] = pspy_utils.naive_binning(l, bl_ar2[field], binning_file, lmax)
+
+            Rb = {}
+            for spec in spectra:
+                X,Y = spec
+                nbs_ar1xar1[spec] *= bb_ar1[X] * bb_ar1[Y]
+                nbs_ar1xar2[spec] *= bb_ar1[X] * bb_ar2[Y]
+                nbs_ar2xar2[spec] *= bb_ar2[X] * bb_ar2[Y]
+                
+                Rb[spec] = nbs_ar1xar2[spec] / np.sqrt(np.abs(nbs_ar1xar1[spec] * nbs_ar2xar2[spec]))
+
+            if ar1 == ar2:
+                nlth = interpolate_dict(lb, nbs_ar1xar1, lth, spectra)
+            else:
+                Rlth = interpolate_dict(lb, Rb, lth, spectra)
+                nlth_ar1xar1 = interpolate_dict(lb, nbs_ar1xar1, lth, spectra)
+                nlth_ar2xar2 = interpolate_dict(lb, nbs_ar2xar2, lth, spectra)
+                nlth = {}
+                for spec in spectra:
+                    nlth[spec] = Rlth[spec] * np.sqrt(nlth_ar1xar1[spec] * nlth_ar2xar2[spec])
+
+
+            for spec in spectra:
+                plt.figure(figsize=(12,12))
+                plt.plot(lth,
+                        nlth[spec],
+                        label="interpolate",
+                        color="lightblue")
+                if ar1 == ar2:
+                    nbs= nbs_ar1xar1[spec]
+                else:
+                    nbs= nbs_ar1xar2[spec]
+
+                plt.plot(lb,
+                         nbs,
+                         ".",
+                         label = f"{sv} {ar1}x{ar2}",
+                         color="red")
+                plt.legend(fontsize=20)
+                plt.savefig(f"{plot_dir}/noise_interpolate_{ar1}x{ar2}_{sv}_{spec}.png", bbox_inches="tight")
+                plt.clf()
+                plt.close()
+
+            spec_name_noise_mean = f"mean_{ar1}x{ar2}_{sv}_noise"
+            so_spectra.write_ps(ps_model_dir + f"/{spec_name_noise_mean}.dat", lth, nlth, type, spectra=spectra)
+
+            if ar2 != ar1:
+                spec_name_noise_mean = f"mean_{ar2}x{ar1}_{sv}_noise"
+                TE, ET, TB, BT, EB, BE = nlth["ET"], nlth["TE"], nlth["BT"], nlth["TB"], nlth["BE"], nlth["EB"]
+                nlth["TE"], nlth["ET"], nlth["TB"], nlth["BT"], nlth["EB"], nlth["BE"] = TE, ET, TB, BT, EB, BE
+                so_spectra.write_ps(ps_model_dir + f"/{spec_name_noise_mean}.dat", lth, nlth, type, spectra=spectra)

project/SO/pISO/python/get_sq_windows_alms.py

@@ -0,0 +1,64 @@
+"""
+This script compute all alms squared windows, it's a necessary step of covariance computation.
+"""
+import sys
+
+import numpy as np
+from pspipe_utils import log, pspipe_list
+from pspy import pspy_utils, so_dict, so_map, so_mpi, sph_tools
+
+
+def mult(map_a, map_b, log):
+    res_a = 1 / map_a.data.pixsize()
+    res_b = 1 / map_b.data.pixsize()
+
+    if res_a == res_b:
+        prod = map_a.copy()
+        prod.data *= map_b.data
+    elif res_a < res_b:
+        log.info("resample map a")
+        prod = map_b.copy()
+        map_a_proj = so_map.car2car(map_a, map_b)
+        prod.data *= map_a_proj.data
+    elif res_b < res_a:
+        log.info("resample map b")
+        prod = map_a.copy()
+        map_b_proj = so_map.car2car(map_b, map_a)
+        prod.data *= map_b_proj.data
+
+    return prod
+
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+surveys = d["surveys"]
+lmax = d["lmax"]
+niter = d["niter"]
+sq_win_alms_dir = d["sq_win_alms_dir"]
+
+pspy_utils.create_directory(sq_win_alms_dir)
+
+n_sq_alms, sv1_list, ar1_list, sv2_list, ar2_list = pspipe_list.get_spectra_list(d)
+
+
+log.info(f"number of sq win alms to compute : {n_sq_alms}")
+so_mpi.init(True)
+subtasks = so_mpi.taskrange(imin=0, imax=n_sq_alms - 1)
+print(subtasks)
+for task in subtasks:
+    task = int(task)
+    sv1, ar1, sv2, ar2 = sv1_list[task], ar1_list[task], sv2_list[task], ar2_list[task]
+
+    log.info(f"[{task:02d}] Computing map2alm for {sv1}_{ar1}x{sv2}_{ar2}...")
+
+    win_T1 = so_map.read_map(d["window_T_%s_%s" % (sv1, ar1)])
+    win_T2 = so_map.read_map(d["window_T_%s_%s" % (sv2, ar2)])
+
+    sq_win = mult(win_T1, win_T2, log)
+    # sq_win = win_T1.copy()
+    # sq_win.data[:] *= win_T2.data[:]
+    sqwin_alm = sph_tools.map2alm(sq_win, niter=niter, lmax=lmax)
+
+    np.save(f"{sq_win_alms_dir}/alms_{sv1}_{ar1}x{sv2}_{ar2}.npy", sqwin_alm)

project/SO/pISO/python/montecarlo/mc_analysis.py

@@ -0,0 +1,96 @@
+"""
+This script analyze the simulations generated by mc_mnms_get_spectra_from_nlms.py and stored in "sim_spec_dir" specified in the dictionary file.
+it estimates the mean and numerical std from the simulations
+"""
+
+
+from pspy import pspy_utils, so_dict, so_spectra
+from pspipe_utils import log
+import numpy as np
+import sys
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+type = d["type"]
+surveys = d["surveys"]
+iStart = d["iStart"]
+iStop = d["iStop"]
+lmax = d["lmax"]
+sim_alm_dtype = d["sim_alm_dtype"]
+
+if sim_alm_dtype == "complex64":
+    spec_dtype = np.float32
+elif sim_alm_dtype == "complex128":
+    spec_dtype = np.float64
+
+
+sim_spec_dir = d["sim_spec_dir"]
+mcm_dir = d["mcm_dir"]
+mc_dir = d["montecarlo_dir"]
+cov_dir = d["cov_dir"]
+
+pspy_utils.create_directory(mc_dir)
+pspy_utils.create_directory(cov_dir)
+
+spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+
+for kind in ["cross", "noise", "auto"]:
+
+    vec_list = []
+
+    for iii in range(iStart, iStop + 1):
+        log.info(f"reading sim spectra {kind} {iii}")
+
+        vec = []
+
+        for spec in spectra:
+            for id_sv1, sv1 in enumerate(surveys):
+                arrays_1 = d[f"arrays_{sv1}"]
+                for id_ar1, ar1 in enumerate(arrays_1):
+                    for id_sv2, sv2 in enumerate(surveys):
+                        arrays_2 = d[f"arrays_{sv2}"]
+                        for id_ar2, ar2 in enumerate(arrays_2):
+
+                            if  (id_sv1 == id_sv2) & (id_ar1 > id_ar2) : continue
+                            if  (id_sv1 > id_sv2) : continue
+                            if (sv1 != sv2) & (kind == "noise"): continue
+                            if (sv1 != sv2) & (kind == "auto"): continue
+
+                            spec_name = f"{type}_{sv1}_{ar1}x{sv2}_{ar2}_{kind}_%05d" %  iii
+
+                            lb, Db = so_spectra.read_ps(sim_spec_dir + f"/{spec_name}.dat", spectra=spectra)
+
+                            n_bins = len(lb)
+                            vec = np.append(vec, Db[spec])
+                            
+        vec_list += [vec.astype(spec_dtype)]
+
+    log.info(f"computing monte carlo mean and std")
+
+    mean_vec = np.mean(vec_list, axis=0)
+    std_vec = np.std(vec_list, axis=0)
+
+    log.info(f"saving files on disk")
+
+    id_spec = 0
+    for spec in spectra:
+        for id_sv1, sv1 in enumerate(surveys):
+            arrays_1 = d[f"arrays_{sv1}"]
+            for id_ar1, ar1 in enumerate(arrays_1):
+                for id_sv2, sv2 in enumerate(surveys):
+                    arrays_2 = d[f"arrays_{sv2}"]
+                    for id_ar2, ar2 in enumerate(arrays_2):
+
+                        if  (id_sv1 == id_sv2) & (id_ar1 > id_ar2) : continue
+                        if  (id_sv1 > id_sv2) : continue
+                        if (sv1 != sv2) & (kind == "noise"): continue
+                        if (sv1 != sv2) & (kind == "auto"): continue
+
+                        mean = mean_vec[id_spec * n_bins:(id_spec + 1) * n_bins]
+                        std = std_vec[id_spec * n_bins:(id_spec + 1) * n_bins]
+
+                        np.savetxt(f"{mc_dir}/spectra_{spec}_{sv1}_{ar1}x{sv2}_{ar2}_{kind}.dat", np.array([lb, mean, std]).T)
+
+                        id_spec += 1