add script to produce relevant figure

project/data_analysis/python/paper_plots/results_plot_TT.py

@@ -0,0 +1,83 @@
+"""
+This script plot the combined TT spectra
+"""
+
+from pspy import so_dict, so_spectra, so_cov, pspy_utils
+from pspipe_utils import  log
+import numpy as np
+import pylab as plt
+import sys, os
+import scipy.stats as ss
+from matplotlib import rcParams, gridspec
+
+rcParams["xtick.labelsize"] = 16
+rcParams["ytick.labelsize"] = 16
+rcParams["axes.labelsize"] = 20
+rcParams["axes.titlesize"] = 20
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+show_220 = False
+
+bestfit_dir = "best_fits"
+combined_spec_dir = "combined_spectra"
+plot_dir = "plots/combined_spectra/"
+pspy_utils.create_directory(plot_dir)
+
+type = d["type"]
+
+########################################################################################
+spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+
+if show_220 == True:
+    case_list = ["150x150", "90x90", "90x150", "90x220", "150x220", "220x220"]
+else:
+    case_list = ["150x150", "90x90", "90x150", "90x220", "150x220"]
+
+########################################################################################
+
+lth, Dlth = so_spectra.read_ps(f"{bestfit_dir}/cmb.dat", spectra=spectra)
+
+fig = plt.figure(figsize=(18, 10))
+gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1])
+
+color_list = ["blue", "green", "orange", "red", "gray", "cyan"]
+
+count = 0
+for color, case in zip(color_list, case_list):
+
+    lb_ml, vec_ml, sigma_ml = np.loadtxt(f"{combined_spec_dir}/{type}_{case}_TT.dat", unpack=True)
+    lb_ml, vec_th_ml = np.loadtxt(f"{combined_spec_dir}/bestfit_{case}_TT.dat", unpack=True)
+    cov_ml = np.load(f"{combined_spec_dir}/cov_{case}_TT.npy")
+    corr_ml = so_cov.cov2corr(cov_ml, remove_diag=True)
+
+    chi2 = (vec_ml - vec_th_ml) @ np.linalg.inv(cov_ml) @ (vec_ml - vec_th_ml)
+    ndof = len(lb_ml)
+    pte = 1 - ss.chi2(ndof).cdf(chi2)
+
+    ax1 = plt.subplot(gs[0])
+    ax1.semilogy()
+    ax1.set_ylabel(r"$D_{\ell} [\mu K^{2}] $", fontsize=19)
+    ax1.plot(lth, Dlth["TT"], color="black", alpha=0.6, linestyle="--")
+    ax1.errorbar(lb_ml-10+count*3, vec_ml, sigma_ml, fmt=".", color=color, label=f"{case}, p = {pte:.3f}")
+    ax1.errorbar(lb_ml-10+count*3, vec_th_ml, color=color)
+    ax1.legend(fontsize=16)
+    ax1.set_xticks([])
+    ax1.set_ylim(20, 6000)
+    ax1.set_xlim(0, 8000)
+    ax2 = plt.subplot(gs[1])
+
+    ax2.errorbar(lb_ml-10+count*4, (vec_ml - vec_th_ml) * lb_ml , sigma_ml * lb_ml, fmt=".", color=color)
+    ax2.set_xlim(0, 8000)
+    ax2.set_xlabel(r"$\ell$", fontsize=19)
+    ax2.set_ylabel(r"$\ell (D_{\ell} - D^{\rm th}_{\ell}) [\mu K^{2}] $", fontsize=19)
+    ax2.plot(lth, lth*0, color="black", alpha=0.6, linestyle="--")
+
+    count += 1
+
+plt.subplots_adjust(wspace=0, hspace=0)
+plt.savefig(f"{plot_dir}/all_spectra_TT.png", bbox_inches="tight")
+plt.clf()
+plt.close()

project/data_analysis/python/paper_plots/results_plot_combined_spectra.py

@@ -0,0 +1,111 @@
+"""
+This script plot the combined DR6 power spectra 
+"""
+
+from pspy import so_dict, pspy_utils, so_spectra, so_cov
+from pspipe_utils import pspipe_list, log, best_fits
+import numpy as np
+import pylab as plt
+import sys, os
+import scipy.stats as ss
+from matplotlib import rcParams
+
+
+rcParams["xtick.labelsize"] = 16
+rcParams["ytick.labelsize"] = 16
+rcParams["axes.labelsize"] = 20
+rcParams["axes.titlesize"] = 20
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+combined_spec_dir = "combined_spectra"
+bestfit_dir = "best_fits"
+plot_dir = "plots/combined_spectra/"
+pspy_utils.create_directory(plot_dir)
+
+type = d["type"]
+
+########################################################################################
+
+spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+selected_spectra_list = [["TT"], ["TE", "ET"], ["TB", "BT"], ["EB", "BE"], ["EE"], ["BB"]]
+
+case_list = ["all", "150x150", "90x90", "90x150", "90x220", "150x220", "220x220"]
+
+
+ylim = {}
+ylim["TT"] = [20, 6500]
+ylim["TE"] = [-150, 150]
+ylim["TB"] = [-10, 10]
+ylim["EE"] = [-5, 45]
+ylim["EB"] = [-1, 1]
+ylim["BB"] = [-1, 1]
+
+ylim_res = {}
+ylim_res["TT"] = [-40, 40]
+ylim_res["TE"] = [-10, 10]
+ylim_res["TB"] = [-10, 10]
+ylim_res["EE"] = [-2, 2]
+ylim_res["EB"] = [-1, 1]
+ylim_res["BB"] = [-1, 1]
+########################################################################################
+
+lth, Dlth = so_spectra.read_ps(f"{bestfit_dir}/cmb.dat", spectra=spectra)
+
+for spec_select in selected_spectra_list:
+    s_name = spec_select[0]
+    for case in case_list:
+
+        if (s_name == "TT") & (case == "all"):
+            continue
+
+        if ("220" in case) & (s_name != "TT"): continue
+
+        lb_ml, vec_ml, sigma_ml = np.loadtxt(f"{combined_spec_dir}/{type}_{case}_{s_name}.dat", unpack=True)
+        lb_ml, vec_th_ml = np.loadtxt(f"{combined_spec_dir}/bestfit_{case}_{s_name}.dat", unpack=True)
+        cov_ml = np.load(f"{combined_spec_dir}/cov_{case}_{s_name}.npy")
+        corr_ml = so_cov.cov2corr(cov_ml, remove_diag=True)
+
+        chi2 = (vec_ml - vec_th_ml) @ np.linalg.inv(cov_ml) @ (vec_ml - vec_th_ml)
+        ndof = len(lb_ml)
+        pte = 1 - ss.chi2(ndof).cdf(chi2)
+
+
+        plt.figure(figsize=(15, 8))
+        plt.suptitle(f"{s_name} {case}", fontsize=22)
+        plt.subplot(2,2,1)
+
+        if (s_name == "TT"): plt.semilogy()
+
+        plt.plot(lth, Dlth[s_name], color="gray", alpha=0.6, linestyle="--")
+        plt.errorbar(lb_ml, vec_ml, sigma_ml, fmt=".")
+        plt.errorbar(lb_ml, vec_th_ml)
+        plt.xlabel(r"$\ell$", fontsize=19)
+        plt.ylabel(r"$D_{\ell} [\mu K^{2}] $", fontsize=19)
+
+        plt.ylim(ylim[s_name])
+
+        if s_name != "TT":
+            plt.xlim(0,4000)
+        plt.subplot(2,2,2)
+        plt.title("Correlation matrix")
+        plt.imshow(corr_ml, cmap="seismic")
+        plt.xticks(np.arange(len(lb_ml))[::3], lb_ml[::3], rotation=90)
+        plt.yticks(np.arange(len(lb_ml))[::3], lb_ml[::3])
+        plt.tight_layout()
+        plt.colorbar()
+        plt.subplot(2,2,3)
+        plt.errorbar(lb_ml, (vec_ml - vec_th_ml) , sigma_ml , label=r"$\chi^{2}$=%.2f, pte = %.4f" % (chi2, pte), fmt=".")
+        plt.plot(lb_ml, lb_ml * 0 )
+        plt.plot(lth, Dlth[s_name] * 0, color="gray", alpha=0.6, linestyle="--")
+        plt.xlabel(r"$\ell$", fontsize=19)
+        plt.ylabel(r"$D_{\ell} - D^{\rm th}_{\ell} [\mu K^{2}]$", fontsize=19)
+        plt.ylim(ylim_res[s_name])
+        if s_name != "TT":
+            plt.xlim(0,4000)
+        plt.legend(fontsize=16)
+        plt.savefig(f"{plot_dir}/spectra_{case}_{s_name}.png", bbox_inches="tight")
+        plt.clf()
+        plt.close()

project/data_analysis/python/paper_plots/results_plot_with_DR4.py

@@ -0,0 +1,117 @@
+"""
+The script compares the dr6 results with the DR4 (https://arxiv.org/abs/2007.07289) results
+"""
+from pspy import so_dict, pspy_utils
+from pspipe_utils import external_data
+import matplotlib.pyplot as plt
+import numpy as np
+import sys, os
+import scipy.stats as ss
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+
+spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+combined_spectra_dir = "combined_spectra"
+plot_dir = "plots/combined_spectra/"
+pspy_utils.create_directory(plot_dir)
+
+
+pow_res = {}
+pow_res["TT"] = 1.5
+pow_res["TE"] = 0
+pow_res["EE"] = -1
+
+# Choi best fits
+yp = {}
+yp["90"] = 0.9853
+yp["150"] = 0.9705
+
+ylim_res = {}
+ylim_res["TT", "90x90"] = [-2 * 10 ** 7, 2 * 10 ** 7]
+ylim_res["TT", "90x150"] =  [- 8.66 * 10 ** 6, 7 * 10 ** 6]
+ylim_res["TE", "90x90"] = [-28, 55]
+ylim_res["EE", "90x90"] = [-0.01, 0.015]
+
+ylim = {}
+ylim["EE", "90x90"] = [-20, 50]
+ylim["TT", "90x90"] = [20, 2000]
+ylim["TT", "90x150"] = [20, 2000]
+ylim["TT", "150x150"] = [20, 4000]
+
+for spec in ["TT", "TE", "EE"]:
+    fp_list, ell, cl_deep, err_deep = external_data.get_choi_spectra(spec, survey="deep", return_Dl=True)
+    _, _ , cl_wide, err_wide = external_data.get_choi_spectra(spec, survey="wide", return_Dl=True)
+
+    for fp in fp_list:
+        print(fp)
+        f0, f1 = fp.split("x")
+        if fp == "150x90": continue
+
+        l = ell[fp]
+        if spec == "TT":
+            # for TT only use deep because that's similar masking as DR6
+            dr4_label = "DR4 deep"
+            cl_choi = cl_deep[fp]
+            err_choi = err_deep[fp]
+        else:
+            dr4_label = "DR4 deep + wide"
+
+            # otherwise do an inverse variance combination
+            err_choi = np.sqrt(1 / ( 1 / err_deep[fp] ** 2 + 1 / err_wide[fp] ** 2))
+            cl_choi = (cl_deep[fp] / err_deep[fp] ** 2 + cl_wide[fp] / err_wide[fp] ** 2) * err_choi ** 2
+
+        if spec == "TE":
+            cl_choi, err_choi = cl_choi / yp[f1], err_choi / yp[f1]
+        elif spec == "EE":
+            cl_choi, err_choi = cl_choi / (yp[f0] * yp[f1]), err_choi / (yp[f0] * yp[f1])
+
+        lb, Db, sigmab = np.loadtxt(f"{combined_spectra_dir}/Dl_{fp}_{spec}.dat", unpack=True)
+
+        id = np.where((lb >= l[0]) & (lb <= l[-1]))
+        lb, Db, sigmab = lb[id], Db[id], sigmab[id]
+
+        id = np.where((l >= lb[0]) & (l <= lb[-1]))
+        l, cl_choi, err_choi = l[id], cl_choi[id], err_choi[id]
+
+        res = Db - cl_choi
+        err_res = np.sqrt(sigmab ** 2 + err_choi ** 2)
+        chi2 = np.sum(res ** 2/err_res ** 2)
+        ndof = len(l)
+        pte = 1 - ss.chi2(ndof).cdf(chi2)
+
+
+        plt.figure(figsize=(12,8))
+        plt.suptitle(f"{spec} {f0}x{f1}", fontsize=24)
+        plt.subplot(2, 1, 1)
+        if spec == "TT":
+            plt.semilogy()
+        plt.errorbar(lb, Db, sigmab, fmt=".-", label="DR6", alpha=0.5)
+        plt.errorbar(l,  cl_choi, err_choi, fmt=".-", label=dr4_label, alpha=0.5)
+        if spec == "EE":
+            plt.errorbar(l,  l*0, fmt="--",  alpha=0.5, color="gray")
+
+        plt.ylabel(r"$D_{\ell}$", fontsize=18)
+        plt.xlabel(r"$\ell$", fontsize=18)
+        try:
+            ylim0 = ylim[spec, fp][0]
+            ylim1 = ylim[spec, fp][1]
+            plt.ylim(ylim0, ylim1)
+        except:
+            pass
+        plt.legend(fontsize=15)
+        plt.subplot(2, 1, 2)
+        plt.ylabel(r"$ell^{%s} (D^{DR6}_{\ell} - D^{DR4}_{\ell}) \ $" % pow_res[spec], fontsize=15)
+        plt.xlabel(r"$\ell$", fontsize=18)
+        plt.errorbar(l,  res * l ** pow_res[spec], err_res * l ** pow_res[spec], fmt=".", label=r"$\chi^{2}$/Dof= %.02f/%d, pte:%0.3f"% (chi2, ndof, pte))
+        try:
+            ylim0 = ylim_res[spec, fp][0]
+            ylim1 = ylim_res[spec, fp][1]
+            plt.ylim(ylim0, ylim1)
+        except:
+            pass
+        plt.plot(l, l * 0, color="gray")
+        plt.legend(fontsize=15)
+        plt.savefig(f"{plot_dir}/dr6_choi_{spec}_{fp}.png", bbox_inches="tight")
+        plt.clf()
+        plt.close()

project/data_analysis/python/paper_plots/results_plot_with_planck.py

@@ -0,0 +1,79 @@
+"""
+This script plot the combined dr6 spectra together with planck
+"""
+
+from pspy import so_dict, so_spectra, pspy_utils
+from pspipe_utils import  log, best_fits
+import numpy as np
+import pylab as plt
+import sys, os
+from matplotlib import rcParams
+import pspipe_utils
+
+
+rcParams["xtick.labelsize"] = 16
+rcParams["ytick.labelsize"] = 16
+rcParams["axes.labelsize"] = 20
+rcParams["axes.titlesize"] = 20
+
+d = so_dict.so_dict()
+d.read_from_file(sys.argv[1])
+log = log.get_logger(**d)
+
+combined_spec_dir = "combined_spectra"
+bestfit_dir = "best_fits"
+
+plot_dir = "plots/combined_spectra/"
+pspy_utils.create_directory(plot_dir)
+
+type = d["type"]
+
+planck_data_path = os.path.join(os.path.dirname(os.path.abspath(pspipe_utils.__file__)), "data/spectra/planck")
+
+########################################################################################
+selected_spectra_list = [["TE", "ET"], ["EE"]]
+########################################################################################
+
+ylim = {}
+
+ylim["TE"] = [0, 60000]
+ylim["TE"] = [-105000, 75000]
+ylim["EE"] = [0, 45]
+fac = {}
+fac["TE"] = 1
+fac["EE"] = 0
+fac["TT"] = 0
+
+spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]
+
+lth, Dlth = so_spectra.read_ps(f"{bestfit_dir}/cmb.dat", spectra=spectra)
+
+for spec_select in selected_spectra_list:
+    s_name = spec_select[0]
+
+    lb_ml, vec_ml, sigma_ml = np.loadtxt(f"{combined_spec_dir}/{type}_all_{s_name}.dat", unpack=True)
+    lb_ml, vec_th_ml = np.loadtxt(f"{combined_spec_dir}/bestfit_all_{s_name}.dat", unpack=True)
+    cov_ml = np.load(f"{combined_spec_dir}/cov_all_{s_name}.npy")
+
+    plt.figure(figsize=(15, 10))
+    lp, Dlp, sigmap, _, _ = np.loadtxt(f"{planck_data_path}/COM_PowerSpect_CMB-{s_name}-binned_R3.02.txt", unpack=True)
+    if s_name == "TT": plt.semilogy()
+
+    plt.xlim(0,4000)
+    plt.ylim(ylim[s_name])
+    plt.errorbar(lp, Dlp * lp ** fac[s_name], sigmap * lp ** fac[s_name], fmt=".", color="royalblue", markersize=2, alpha=0.6, label="Planck PR3")
+    plt.errorbar(lb_ml, vec_ml *  lb_ml ** fac[s_name], sigma_ml * lb_ml ** fac[s_name] , fmt=".", color="red", markersize=2, label="ACT DR6")
+    plt.plot(lth, Dlth[s_name] * lth ** fac[s_name], color="gray", alpha=0.4)
+    plt.xlabel(r"$\ell$", fontsize=30)
+
+    if fac[s_name] == 0:
+        plt.ylabel(r"$D_{\ell}$", fontsize=30)
+    if fac[s_name] == 1:
+        plt.ylabel(r"$\ell D_{\ell}$", fontsize=30)
+    if fac[s_name] > 1:
+        plt.ylabel(r"$\ell^{%s}D_{\ell}$" % fac[s_name], fontsize=30)
+
+    plt.legend(fontsize=25)
+    plt.savefig(f"{plot_dir}/all_spectra_{s_name}_with_planck.png", bbox_inches="tight")
+    plt.clf()
+    plt.close()