last commit before christmas

Author: Merry Duparc

project/SO/pISO/paramfiles/dr6xlegacyxdeep56_20251119.dict

@@ -30,27 +30,29 @@ def get_proj_pattern(test, map_set, ref_map_set):
 # log calib infos from calib yaml file
 with open(d['calib_yaml'], "r") as f:
     calib_dict: dict = yaml.safe_load(f)
-calib_infos = calib_dict['get_calibs.py']
+calib_infos: dict = calib_dict['get_calibs.py']
 
 planck_corr = calib_infos['planck_corr']
 subtract_bf_fg = calib_infos['subtract_bf_fg']
 
+data_dir = d["data_dir"]
+plot_dir = d["plots_dir"]
 spec_dir = d['spec_dir']
 cov_dir = d['cov_dir']
-bestfit_dir = d['bestfits_dir']
+bestfit_dir = d['best_fits_dir']
 
 # Create output dirs
-output_dir = d['calib_dir']
-residual_output_dir = f"{output_dir}/residuals"
-plot_output_dir = f"{output_dir}/plots"
-chains_dir = f"{output_dir}/chains"
+calib_dir = d['calib_dir']
+residual_output_dir = f"{calib_dir}/residuals"
+chains_dir = f"{calib_dir}/chains"
+plot_output_dir = plot_dir + '/calib/'
 
 if planck_corr:
     spec_dir = "spectra_leak_corr_planck_bias_corr"
-    output_dir += "_planck_bias_corrected"
+    calib_dir += "_planck_bias_corrected"
 
 if subtract_bf_fg:
-    output_dir += "_fg_sub"
+    calib_dir += "_fg_sub"
 
 _, _, lb, _ = pspy_utils.read_binning_file(d["binning_file"], d["lmax"])
 n_bins = len(lb)
@@ -62,7 +64,7 @@ def get_proj_pattern(test, map_set, ref_map_set):
 else:
     modes = spectra
 
-pspy_utils.create_directory(output_dir)
+pspy_utils.create_directory(calib_dir)
 pspy_utils.create_directory(residual_output_dir)
 pspy_utils.create_directory(chains_dir)
 pspy_utils.create_directory(plot_output_dir)
@@ -78,110 +80,104 @@ def get_proj_pattern(test, map_set, ref_map_set):
 tests = ["AxA-AxB", "AxA-BxB", "BxB-AxB"]
 
 for test in tests:
-    for sv in calib_infos["surveys_to_calib"]:
-        for ar in d[f"arrays_{sv}"]:
+    # Loop over each map_set-ref_map_set combination 
+    # (so you can measure calib on a same map_set with different ref_map_set)
+    for map_set, ref_map_set in calib_infos[f"ref_map_sets"].items():  
+        name, proj_pattern = get_proj_pattern(test, map_set, ref_map_set)
 
-            map_set = f"{sv}_{ar}"
-            ref_map_set = calib_infos[f"ref_map_sets"][map_set]
-
-            name, proj_pattern = get_proj_pattern(test, map_set, ref_map_set)
-            
-            spectra_for_cal = [
-                (ref_map_set, ref_map_set, "TT"),
-                (ref_map_set, map_set, "TT"),
-                (map_set, map_set, "TT"),
-                # (map_set, map_set, "TT"),
-                # (map_set, ref_map_set, "TT"),
-                # (ref_map_set, ref_map_set, "TT"),
-            ]
-
-            # Load spectra and cov
-            ps_dict = {}
-            cov_dict = {}
-            for i, (ms1, ms2, m1) in enumerate(spectra_for_cal):
-                _, ps = so_spectra.read_ps(f"{spec_dir}/Dl_{ms1}x{ms2}_cross.dat",
-                                        spectra=spectra)
-                    
-                ps_dict[ms1, ms2, m1] = ps[m1]
+        spectra_for_cal = [
+            (ref_map_set, ref_map_set, "TT"),
+            (ref_map_set, map_set, "TT"),
+            (map_set, map_set, "TT"),
+            # (map_set, map_set, "TT"),
+            # (map_set, ref_map_set, "TT"),
+            # (ref_map_set, ref_map_set, "TT"),
+        ]
+
+        # Load spectra and cov
+        ps_dict = {}
+        cov_dict = {}
+        for i, (ms1, ms2, m1) in enumerate(spectra_for_cal):
+            _, ps = so_spectra.read_ps(f"{spec_dir}/Dl_{ms1}x{ms2}_cross.dat",
+                                    spectra=spectra)
 
-                if (m1 == "TT") & (subtract_bf_fg):
-                    log.info(f"remove fg {m1}  {ms1} x {ms2}")
-                    l_fg, bf_fg = so_spectra.read_ps(f"{bestfit_dir}/fg_{ms1}x{ms2}.dat", spectra=spectra)
-                    _, bf_fg_TT_binned = pspy_utils.naive_binning(l_fg, bf_fg["TT"], d["binning_file"], d["lmax"])
-                    ps_dict[ms1, ms2, m1] -= bf_fg_TT_binned
+            ps_dict[ms1, ms2, m1] = ps[m1]
+            
+            if (m1 == "TT") & (subtract_bf_fg):
+                log.info(f"remove fg {m1}  {ms1} x {ms2}")
+                l_fg, bf_fg = so_spectra.read_ps(f"{bestfit_dir}/fg_{ms1}x{ms2}.dat", spectra=spectra)
+                _, bf_fg_TT_binned = pspy_utils.naive_binning(l_fg, bf_fg["TT"], d["binning_file"], d["lmax"])
+                ps_dict[ms1, ms2, m1] -= bf_fg_TT_binned
 
 
-                for j, (ms3, ms4, m2) in enumerate(spectra_for_cal):
-                    if j < i: continue
-                    cov = np.load(f"{cov_dir}/analytic_cov_{ms1}x{ms2}_{ms3}x{ms4}.npy")
-                    cov = so_cov.selectblock(cov, modes, n_bins = n_bins, block = m1 + m2)
-                    cov_dict[(ms1, ms2, m1), (ms3, ms4, m2)] = cov
+            for j, (ms3, ms4, m2) in enumerate(spectra_for_cal):
+                if j < i: continue
+                cov = np.load(f"{cov_dir}/analytic_cov_{ms1}x{ms2}_{ms3}x{ms4}.npy")
+                cov = so_cov.selectblock(cov, modes, n_bins = n_bins, block = m1 + m2)
+                cov_dict[(ms1, ms2, m1), (ms3, ms4, m2)] = cov
 
-            # Concatenation
-            spec_vec, full_cov = consistency.append_spectra_and_cov(ps_dict, cov_dict, spectra_for_cal)
+        # Concatenation
+        spec_vec, full_cov = consistency.append_spectra_and_cov(ps_dict, cov_dict, spectra_for_cal)
 
-            # Save and plot residuals before calibration
-            res_spectrum, res_cov = consistency.project_spectra_vec_and_cov(spec_vec, full_cov, proj_pattern)
-            np.savetxt(f"{residual_output_dir}/residual_{name}_before.dat", np.array([lb, res_spectrum]).T)
-            np.savetxt(f"{residual_output_dir}/residual_cov_{name}.dat", res_cov)
+        # Save and plot residuals before calibration
+        res_spectrum, res_cov = consistency.project_spectra_vec_and_cov(spec_vec, full_cov, proj_pattern)
+        np.savetxt(f"{residual_output_dir}/residual_{name}_before.dat", np.array([lb, res_spectrum]).T)
+        np.savetxt(f"{residual_output_dir}/residual_cov_{name}.dat", res_cov)
 
-            lmin, lmax = calib_infos['ell_ranges'][map_set]
-            id = np.where((lb >= lmin) & (lb <= lmax))[0]
-            consistency.plot_residual(lb, res_spectrum, {"analytical": res_cov}, "TT", f"{map_set} {test}",
-                                    f"{plot_output_dir}/residual_{name}_before",
-                                    lrange=id, l_pow=1, ylims=calib_infos['y_lims_plot_TT'])
+        lmin, lmax = calib_infos['ell_ranges'][map_set]
+        id = np.where((lb >= lmin) & (lb <= lmax))[0]
+        consistency.plot_residual(lb, res_spectrum, {"analytical": res_cov}, "TT", f"{map_set} {test}",
+                                f"{plot_output_dir}/residual_{name}_before",
+                                lrange=id, l_pow=1, ylims=calib_infos['y_lims_plot_TT'])
 
-            # Calibrate the spectra
-            cal_mean, cal_std = consistency.get_calibration_amplitudes(spec_vec, full_cov,
-                                                                    proj_pattern, "TT", id,
-                                                                    f"{chains_dir}/{name}")
+        # Calibrate the spectra
+        cal_mean, cal_std = consistency.get_calibration_amplitudes(spec_vec, full_cov,
+                                                                proj_pattern, "TT", id,
+                                                                f"{chains_dir}/{name}",)
 
 
-            results_dict[test, ar] = {"multipole_range": calib_infos['ell_ranges'][map_set],
-                                    "ref_map_set": ref_map_set,
-                                    "calibs": [cal_mean, cal_std]}
+        results_dict[test, map_set] = {"multipole_range": calib_infos['ell_ranges'][map_set],
+                                "ref_map_set": ref_map_set,
+                                "calibs": [cal_mean, cal_std]}
 
-            calib_vec = np.array([cal_mean**2, cal_mean, 1])
-            res_spectrum, res_cov = consistency.project_spectra_vec_and_cov(spec_vec, full_cov,
-                                                                            proj_pattern,
-                                                                            calib_vec = calib_vec)
-    
-            np.savetxt(f"{residual_output_dir}/residual_{name}_after.dat", np.array([lb, res_spectrum]).T)
-            consistency.plot_residual(lb, res_spectrum, {"analytical": res_cov}, "TT", f"{map_set} {test}",
-                                    f"{plot_output_dir}/residual_{name}_after",
-                                    lrange=id, l_pow=1, ylims=calib_infos[f"y_lims_plot_TT"])
+        calib_vec = np.array([cal_mean**2, cal_mean, 1])
+        res_spectrum, res_cov = consistency.project_spectra_vec_and_cov(spec_vec, full_cov,
+                                                                        proj_pattern,
+                                                                        calib_vec = calib_vec)
+
+        np.savetxt(f"{residual_output_dir}/residual_{name}_after.dat", np.array([lb, res_spectrum]).T)
+        consistency.plot_residual(lb, res_spectrum, {"analytical": res_cov}, "TT", f"{map_set} {test}",
+                                f"{plot_output_dir}/residual_{name}_after",
+                                lrange=id, l_pow=1, ylims=calib_infos[f"y_lims_plot_TT"])
 
 
 
 # plot the cal factors
 color_list =  ["blue", "red", "green"]
 
-for sv in calib_infos['surveys_to_calib']:
-    for i, ar in enumerate(d[f"arrays_{sv}"]):
-        map_set = f"{sv}_{ar}"
-        ref_map_set = calib_infos[f"ref_map_sets"][map_set]
-        print(f"**************")
-        print(f"calibration {map_set} with {ref_map_set}")
-
-        for j, test in enumerate(tests):
-            cal, std = results_dict[test, ar]["calibs"]
-            print(f"{test}, cal: {cal}, sigma cal: {std}")
-
-            plt.errorbar(i + 0.9 + j * 0.1, cal, std, label = test,
-                        color = color_list[j], marker = ".",
-                        ls = "None",
-                        markersize=6.5,
-                        markeredgewidth=2)
-
-        if i == 0:
-            plt.legend(fontsize = 15)
-
-    x = np.arange(1, len(d[f"arrays_{sv}"]) + 1)
-    plt.xticks(x, d[f"arrays_{sv}"])
-    plt.ylim(0.967, 1.06)
-    plt.tight_layout()
-    plt.savefig(f"{plot_output_dir}/calibs_summary.pdf", bbox_inches="tight")
-    plt.clf()
-    plt.close()
-
-pickle.dump(results_dict, open(f"{output_dir}/calibs_dict.pkl", "wb"))
+for i, (map_set, ref_map_set) in enumerate(calib_infos[f"ref_map_sets"].items()):
+    print(f"**************")
+    print(f"calibration {map_set} with {ref_map_set}")
+
+    for j, test in enumerate(tests):
+        cal, std = results_dict[test, map_set]["calibs"]
+        print(f"{test}, cal: {cal}, sigma cal: {std}")
+
+        plt.errorbar(i + 0.9 + j * 0.1, cal, std, label = test,
+                    color = color_list[j], marker = ".",
+                    ls = "None",
+                    markersize=6.5,
+                    markeredgewidth=2)
+
+    if i == 0:
+        plt.legend(fontsize = 15)
+
+x = np.arange(1, len(calib_infos[f"ref_map_sets"]) + 1)
+plt.xticks(x, calib_infos[f"ref_map_sets"].keys())
+    # plt.ylim(0.967, 1.06)
+plt.tight_layout()
+plt.savefig(f"{plot_output_dir}/calibs_summary.pdf", bbox_inches="tight")
+plt.clf()
+plt.close()
+
+pickle.dump(results_dict, open(f"{calib_dir}/calibs_dict.pkl", "wb"))

project/SO/pISO/paramfiles/dr6xlegacyxdeep56_sims_20251201_set0.dict

@@ -10,24 +10,31 @@
 from cobaya.run import run
 import numpy as np
 import sys
+import yaml
 
 d = so_dict.so_dict()
 d.read_from_file(sys.argv[1])
 
+# log calib infos from calib yaml file
+with open(d['poleff_yaml'], "r") as f:
+    calib_dict: dict = yaml.safe_load(f)
+calib_infos: dict = calib_dict['get_polar_eff_LCDM.py']
+
 # Set up directories
-spec_dir = "/global/cfs/cdirs/sobs/users/merrydup/deep56/spectra_0925"
-cov_dir = "/pscratch/sd/m/merrydup/PSpipe_SO/covariances_d56_0925_i1"
-bestfit_dir = "best_fits"
-mcm_dir = "mcms"
+spec_dir = d['spec_dir']
+cov_dir = d['cov_dir']
+bestfit_dir = d['best_fits_dir']
+mcm_dir = d["mcm_dir"]
 planck_corr = False
 use_leakage_cov = False
 
 if planck_corr:
     spec_dir = "spectra_leak_corr_planck_bias_corr"
 
-
-output_dir = "pol_eff_results"
+output_dir = d['poleff_dir']
 pspy_utils.create_directory(output_dir)
+plots_dir = d['plots_dir'] + '/poleff/'
+pspy_utils.create_directory(plots_dir)
 
 ps_filename = f"{bestfit_dir}/cmb.dat"
 
@@ -60,6 +67,8 @@
 
 print(dust_priors)
 
+measure_list = calib_infos['sets_to_measure'] or [f'{sv}_{ar}' for sv in surveys for ar in arrays[sv]]
+
 passbands = {}
 for sv in surveys:
     for ar in arrays[sv]:
@@ -84,11 +93,10 @@ def get_model(cmb_th, fg_th, Bbl, dust_amp, pol_eff, mode):
     return Bbl @ ps_theory
 
 pol_eff_mean, pol_eff_std, dust_mean, dust_std = {}, {}, {}, {}
-for sv in surveys:
-    for ar in arrays[sv]:
+for sv_ar in measure_list:
         for spectrum in ["EE", "TE"]:
 
-            print(spectrum, sv, ar, d["lmin_cal"], d[f"lmax_cal_{sv}_{ar}"])
+            print(spectrum, sv, ar, calib_infos["lmin_cal"], calib_infos[f"lmax_cal_{sv}_{ar}"])
 
             # Load ps and cov
             spec_name = f"{sv}_{ar}x{sv}_{ar}"
@@ -104,7 +112,7 @@ def get_model(cmb_th, fg_th, Bbl, dust_amp, pol_eff, mode):
             cov = so_cov.selectblock(cov, spectra, n_bins=n_bins, block=spectrum+spectrum)
 
             # Multipole cuts
-            id = np.where((lb >= d["lmin_cal"]) & (lb <= d[f"lmax_cal_{sv}_{ar}"]))[0]
+            id = np.where((lb >= calib_infos["lmin_cal"]) & (lb <= calib_infos[f"lmax_cal_{sv}_{ar}"]))[0]
             ps = ps[id]
             cov = cov[np.ix_(id, id)]
             invcov = np.linalg.inv(cov)
@@ -174,7 +182,7 @@ def loglike(pol_eff, dust_amp):
 
             gdplot = gdplt.get_subplot_plotter()
             gdplot.triangle_plot(samples, ["pol_eff", "dust_amp"], filled=True, title_limit=1)
-            plt.savefig(f"{output_dir}/posterior_dist_{spectrum}_{sv}_{ar}.png", dpi=300, bbox_inches="tight")
+            plt.savefig(f"{plots_dir}/posterior_dist_{spectrum}_{sv}_{ar}.png", dpi=300, bbox_inches="tight")
 
 
 for sv in surveys:

project/SO/pISO/python/calibration/get_calibs.py

@@ -16,7 +16,7 @@
 spectra_dir = d['spec_dir']
 noise_dir = d['noise_model_dir']
 cov_dir = d['cov_dir']
-bestfit_dir = d["bestfits_dir"]
+bestfit_dir = d["best_fits_dir"]
 sq_win_alms_dir = d["sq_win_alms_dir"]
 
 pspy_utils.create_directory(cov_dir)

project/SO/pISO/python/calibration/get_polar_eff_LCDM.py

@@ -18,7 +18,9 @@
 
 spec_dir = d['spec_dir']
 cov_dir = d['cov_dir']
-bestfir_dir = d['bestfits_dir']
+bestfir_dir = d['best_fits_dir']
+plots_dir = d['plots_dir'] + '/tf/'
+os.makedirs(plots_dir, exist_ok=True)
 
 tf_dir = d['tf_dir']
 
@@ -88,4 +90,4 @@
             np.savetxt(f"{output_dir}/tf_cov_{sv_ar}.dat", tf_cov)
 
 
-        tf_tools.plot_tf(lb_list, tf_list, tf_err_list, d[f"arrays_{sv}"], f"{output_dir}/tf_estimation.png")
+        tf_tools.plot_tf(lb_list, tf_list, tf_err_list, d[f"arrays_{sv}"], f"{plots_dir}/tf_estimation_{combin}.png")