Skip to content

Commit

Permalink
update pol angle scripts
Browse files Browse the repository at this point in the history
  • Loading branch information
Louis Thibaut committed Jan 30, 2025
1 parent ae19ae2 commit dd8c19f
Show file tree
Hide file tree
Showing 3 changed files with 147 additions and 171 deletions.
202 changes: 54 additions & 148 deletions project/ACT_DR6/python/paper_plots/B_modes/EB_summary.py
View file
Original file line number Diff line number Diff line change
@@ -1,192 +1,98 @@
"""
For this to work you need to have created two folders results_EB_optimal and results_EB_uniform
using results_plot_pol_angle.py for both the uniform and weighted analysis.
Similarly you need to have created combined_spectra_paper_optimal and the corresponding simulations
"""

import numpy as np
import pylab as plt
import pickle
from matplotlib import rcParams
from pspy import so_spectra, pspy_utils, so_cov
from pspipe_utils import pol_angle
import scipy.stats as ss


binning_file = "BIN_ACTPOL_50_4_SC_large_bin_at_low_ell"
lmax = 8500


rcParams["xtick.labelsize"] = 16
rcParams["ytick.labelsize"] = 16
rcParams["axes.labelsize"] = 20
rcParams["axes.titlesize"] = 20
rcParams["xtick.labelsize"] = 30
rcParams["ytick.labelsize"] = 30
rcParams["axes.labelsize"] = 25
rcParams["axes.titlesize"] = 25

latex_par_name = {}
latex_par_name["alpha_pa5_f090"] = r"$\alpha_{pa5 f090}$"
latex_par_name["alpha_pa6_f090"] = r"$\alpha_{pa6 f090}$"
latex_par_name["alpha_pa5_f150"] = r"$\alpha_{pa5 f150}$"
latex_par_name["alpha_pa6_f150"] = r"$\alpha_{pa6 f150}$"
latex_par_name["alpha_pa5_f090"] = r"$\psi_{pa5 f090}$"
latex_par_name["alpha_pa6_f090"] = r"$\psi_{pa6 f090}$"
latex_par_name["alpha_pa5_f150"] = r"$\psi_{pa5 f150}$"
latex_par_name["alpha_pa6_f150"] = r"$\psi_{pa6 f150}$"
latex_par_name["beta_pa5"] = r"\beta_{\rm pa5}"
latex_par_name["beta_pa6"] = r"\beta_{\rm pa6}"
latex_par_name["beta_ACT"] = r"\beta_{\rm ACT}"
latex_par_name["beta_ACT+komatsu"] = r"\beta_{\rm ACT + Komatsu}"
latex_par_name["beta_komatsu"] = r"\beta_{\rm Komatsu}"

cut_list = [ "post_unblinding", "pre_unblinding"]
name_list = [ "baseline multipole cut", "extended multipole cut"]
weight_list = ["uniform", "optimal"]
cut_list = ["post_unblinding", "pre_unblinding"]
cut_list_name = ["baseline", "extended"]
folder_list = ["dr6", "dr6_weight"]
folder_list_name = ["uniform", "optimal"]

angle= {}
for my_w in weight_list:
for my_cut in cut_list:
with open(f"results_EB_{my_w}/angle_{my_cut}.pkl", "rb") as fp:
angle[my_w, my_cut] = pickle.load(fp)
for folder, folder_name in zip(folder_list, folder_list_name):
for cut, cut_name in zip(cut_list, cut_list_name):
with open(f"{folder}/plots/results_EB_paper/angle_{cut}.pkl", "rb") as fp:
angle[folder_name, cut_name] = pickle.load(fp)


alpha_list = ["alpha_pa5_f090", "alpha_pa5_f150", "alpha_pa6_f090", "alpha_pa6_f150"]


x = np.linspace(0,4.5,100)

plt.figure(figsize=(12,6))
plt.figure(figsize=(20, 12))
plt.subplot(1,2,1)
for count, alpha in enumerate(alpha_list):
color_list = ["red", "orange", "blue", "lightblue"]
color_list = ["red", "blue", "orange", "lightblue"]
count_all = 0
for shift2, (name, my_cut) in enumerate(zip(name_list, cut_list)):
for shift1, my_w in enumerate(weight_list):

print(f"{alpha} {name} ({my_w} weighting)", angle[my_w, my_cut][alpha, "mean"], angle[my_w, my_cut][alpha, "std"] )
sh1 = shift1*0.1
sh2 = shift2*0.2
plt.errorbar(count+sh1+sh2,
angle[my_w, my_cut][alpha, "mean"],
angle[my_w, my_cut][alpha, "std"],
label=f"{name} ({my_w} weighting)",

for shift1, folder_name in enumerate(folder_list_name):
for shift2, cut_name in enumerate(cut_list_name):


print(f"{alpha} {cut_name} ({folder_name} weighting)", angle[folder_name, cut_name][alpha, "mean"], angle[folder_name, cut_name][alpha, "std"])
sh1 = shift1 * 0.1
sh2 = shift2 * 0.2
plt.errorbar(count + sh1 + sh2,
angle[folder_name, cut_name][alpha, "mean"],
angle[folder_name, cut_name][alpha, "std"],
label=f"{cut_name} ({folder_name} weight)",
color=color_list[count_all],
fmt="o")

count_all += 1
if count == 0:
plt.legend(fontsize=12)
plt.legend(fontsize=26, loc="upper right")

xticks_list = ["pa5 f090", "pa5 f150", "pa6 f090", "pa6 f150"]
plt.xticks([0.15,1.15,2.15,3.15], xticks_list, rotation=90, fontsize=22)
plt.ylabel(r"$\alpha$", fontsize=32)
plt.tight_layout()
plt.savefig("all_alpha.png")
plt.clf()
plt.close()


plt.xticks([0.15, 1.15, 2.15, 3.15], xticks_list, rotation=90, fontsize=32)
plt.ylabel(r"$\psi \ [\rm deg]$", fontsize=42)
plt.ylim(0,0.5)

x = np.linspace(0,4.5,100)
beta_list = ["beta_pa5", "beta_pa6", "beta_ACT", "beta_ACT+komatsu"]
plt.figure(figsize=(12,8))
beta_list = ["beta_pa5", "beta_pa6", "beta_ACT"]
plt.subplot(1,2,2)

for count, beta in enumerate(beta_list):
count_all = 0
for shift2, (name, my_cut) in enumerate(zip(name_list, cut_list)):
for shift1, my_w in enumerate(weight_list):
print(f"{beta} {name} ({my_w} weighting)", angle[my_w, my_cut][beta, "mean"], angle[my_w, my_cut][beta, "std"] )

sh1 = shift1*0.1
sh2 = shift2*0.2
plt.errorbar(count+sh1+sh2,
angle[my_w, my_cut][beta, "mean"],
angle[my_w, my_cut][beta, "std"],
label=f"{name} ({my_w} weighting)",
for shift1, folder_name in enumerate(folder_list_name):
for shift2, cut_name in enumerate(cut_list_name):

print(f"{beta} {cut_name} ({folder_name} weighting)", angle[folder_name, cut_name][beta, "mean"], angle[folder_name, cut_name][beta, "std"] )

sh1 = shift1 * 0.1
sh2 = shift2 * 0.2
plt.errorbar(count + sh1 + sh2,
angle[folder_name, cut_name][beta, "mean"],
angle[folder_name, cut_name][beta, "std"],
label=f"{cut_name} ({folder_name} weighting)",
color=color_list[count_all],
fmt="o")

count_all += 1

if count == 0:
plt.legend(fontsize=12)
plt.ylim(0,0.5)

xticks_list = ["pa5", "pa6", "ACT"]#, "ACT+Planck"]
plt.xticks([0.15,1.15,2.15], xticks_list, rotation=90, fontsize=32)
plt.ylabel(r"$\hat{\psi} \ [\rm deg]$", fontsize=42)

xticks_list = ["pa5", "pa6", "ACT", "ACT+Planck"]
plt.xticks([0.15,1.15,2.15,3.15], xticks_list, rotation=90, fontsize=22)
plt.ylabel(r"$\beta$", fontsize=32)
plt.tight_layout()
plt.savefig("all_beta.png")
plt.savefig("all_alpha_and_beta.pdf")
plt.clf()
plt.close()

beta_ACT = angle["optimal", "post_unblinding"]["beta_ACT", "mean"]
print(beta_ACT)


spectra = ["TT", "TE", "TB", "ET", "BT", "EE", "EB", "BE", "BB"]

for spec in ["EB", "TB"]:
l, ps_th = so_spectra.read_ps("cmb.dat", spectra=spectra)
l, psth_rot = pol_angle.rot_theory_spectrum(l, ps_th, beta_ACT, beta_ACT)

lb, ps_th_b = pspy_utils.naive_binning(l, psth_rot[spec], binning_file, lmax)

lb_dat, ps, error = np.loadtxt(f"combined_spectra_paper_optimal/Dl_all_{spec}.dat", unpack=True)

id = np.where(lb>=lb_dat[0])

lb, ps_th_b = lb[id], ps_th_b[id]

cov = np.load(f"combined_spectra_paper_optimal/cov_all_{spec}.npy")

chi2 = (ps - ps_th_b) @ np.linalg.inv(cov) @ (ps - ps_th_b)
chi2_null = (ps) @ np.linalg.inv(cov) @ (ps)

pte = 1 - ss.chi2(len(lb)).cdf(chi2)
pte_null = 1 - ss.chi2(len(lb)).cdf(chi2_null)

plt.figure(figsize=(12,6))
plt.errorbar(l, ps_th["TT"]*0, label=f"pte (null): {pte_null*100:.5f} %", color="gray")
plt.errorbar(l, psth_rot[spec], label=f"pte (model): {pte*100:.2f} %", color="red")
plt.errorbar(lb_dat, ps, error, fmt="o", color="royalblue")
plt.xlabel(r"$\ell$", fontsize=25)
plt.ylabel(r"$D^{%s}_{\ell}$" % spec, fontsize=25)
plt.legend(fontsize=14)
if spec == "EB":
plt.ylim(-0.3,0.6)
if spec == "TB":
plt.ylim(-4,6)

plt.xlim(0,3500)
plt.savefig(f"{spec}.png")
plt.clf()
plt.close()

for spec in ["EB", "TB"]:
l, ps_th = so_spectra.read_ps("cmb.dat", spectra=spectra)
l, psth_rot = pol_angle.rot_theory_spectrum(l, ps_th, beta_ACT, beta_ACT)

for count, freq_pair in enumerate(["90x90", "90x150", "150x150"]):
lb, ps_th_b = pspy_utils.naive_binning(l, psth_rot[spec], binning_file, lmax)
lb_dat, ps, error = np.loadtxt(f"combined_spectra_paper_optimal/Dl_{freq_pair}_{spec}.dat", unpack=True)
id = np.where(lb>=lb_dat[0])

lb, ps_th_b = lb[id], ps_th_b[id]

cov = np.load(f"combined_spectra_paper_optimal/cov_{freq_pair}_{spec}.npy")

chi2 = (ps - ps_th_b) @ np.linalg.inv(cov) @ (ps - ps_th_b)
chi2_null = (ps) @ np.linalg.inv(cov) @ (ps)

pte = 1 - ss.chi2(len(lb)).cdf(chi2)
pte_null = 1 - ss.chi2(len(lb)).cdf(chi2_null)

plt.figure(figsize=(12,6))
plt.errorbar(l, ps_th["TT"]*0, label=f"pte (null): {pte_null*100:.5f} %", color="gray")
plt.errorbar(l, psth_rot[spec], label=f"pte (model): {pte*100:.2f} %", color="red")
plt.errorbar(lb_dat, ps, error, fmt=".")

plt.xlabel(r"$\ell$", fontsize=25)
plt.ylabel(r"$D^{%s}_{\ell}$ (%s)" % (spec, freq_pair), fontsize=25)

plt.legend(fontsize=14)
if spec == "EB":
plt.ylim(-0.3,0.6)
if spec == "TB":
plt.ylim(-4,6)

plt.show()
29 changes: 9 additions & 20 deletions project/ACT_DR6/python/paper_plots/B_modes/results_EB_likelihood.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -206,33 +206,22 @@ def loglike(alpha_pa5_f090, alpha_pa5_f150, alpha_pa6_f090, alpha_pa6_f150):


lb_ml = B_modes_utils.get_ml_bins(bin_out_dict, lb)
P_mat = B_modes_utils.get_P_mat(len(vec_EB), lb_ml, bin_out_dict, fig_name=f"{result_dir}/P_mat_all_EB.png")
P_mat = B_modes_utils.get_P_mat(len(vec_EB), lb_ml, bin_out_dict, fig_name=f"{result_dir}/P_mat_all_EB{cut}.png")

cov_ml = covariance.get_max_likelihood_cov(P_mat,
i_cov,
force_sim = True,
check_pos_def = True)

vec_ml_corr = covariance.max_likelihood_spectra(cov_ml,
i_cov,
P_mat,
vec_EB_corr)
vec_ml = covariance.max_likelihood_spectra(cov_ml,
i_cov,
P_mat,
vec_EB)


i_cov_ml = np.linalg.inv(cov_ml)
error_ml = np.sqrt(cov_ml.diagonal())

chi2_ml = vec_ml_corr @ i_cov_ml @ vec_ml_corr
PTE_ml = 1 - ss.chi2(len(lb_ml)).cdf(chi2_ml)

plt.figure(figsize=(12,8))
plt.ylabel(r"$D_{\ell}^{EB}$", fontsize=22)
plt.xlabel(r"$\ell$", fontsize=22)
plt.plot(lb_ml, lb_ml*0)
plt.errorbar(lb_ml, vec_ml_corr, error_ml, fmt=".", label=f"p-value {PTE_ml:.3f}")
plt.legend()
plt.tight_layout()
plt.savefig(f"{paper_plot_dir}/vec_ml_{cut}{tag}.pdf")
plt.clf()
plt.close()

np.savetxt(f"{result_dir}/combined_EB_{cut}.dat", np.transpose([lb_ml, vec_ml, error_ml]))
np.save(f"{result_dir}/combined_cov_EB_{cut}.npy", cov_ml)


Loading

0 comments on commit dd8c19f

Please sign in to comment.