From 565e53a09cfd5afcc53679dbc9c240f5689f5fc5 Mon Sep 17 00:00:00 2001 From: Louis Thibaut Date: Sun, 15 Sep 2024 14:57:08 +0200 Subject: [PATCH] add script to produce relevant figure --- .../python/paper_plots/results_plot_TT.py | 83 +++++++++++++ .../results_plot_combined_spectra.py | 111 +++++++++++++++++ .../paper_plots/results_plot_with_DR4.py | 117 ++++++++++++++++++ .../paper_plots/results_plot_with_planck.py | 79 ++++++++++++ 4 files changed, 390 insertions(+) create mode 100644 project/data_analysis/python/paper_plots/results_plot_TT.py create mode 100644 project/data_analysis/python/paper_plots/results_plot_combined_spectra.py create mode 100644 project/data_analysis/python/paper_plots/results_plot_with_DR4.py create mode 100644 project/data_analysis/python/paper_plots/results_plot_with_planck.py diff --git a/project/data_analysis/python/paper_plots/results_plot_TT.py b/project/data_analysis/python/paper_plots/results_plot_TT.py new file mode 100644 index 00000000..b12f48d9 --- /dev/null +++ b/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() diff --git a/project/data_analysis/python/paper_plots/results_plot_combined_spectra.py b/project/data_analysis/python/paper_plots/results_plot_combined_spectra.py new file mode 100644 index 00000000..565404b2 --- /dev/null +++ b/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() diff --git a/project/data_analysis/python/paper_plots/results_plot_with_DR4.py b/project/data_analysis/python/paper_plots/results_plot_with_DR4.py new file mode 100644 index 00000000..5b378707 --- /dev/null +++ b/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() diff --git a/project/data_analysis/python/paper_plots/results_plot_with_planck.py b/project/data_analysis/python/paper_plots/results_plot_with_planck.py new file mode 100644 index 00000000..b379548b --- /dev/null +++ b/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()