From ca673ccdd3b315125bc11e1d3572a50cbdd88a06 Mon Sep 17 00:00:00 2001 From: Sacha Guerrini Date: Mon, 2 Nov 2026 20:36:37 +0100 Subject: [PATCH 1/2] Push scripts for harmonic space paperplot after aanda resubmission. --- .../harmonic/2025_09_11_psf_leakage_cell.py | 226 ++++++-- .../2025_09_11_validation_cov_cell.py | 89 ++- .../2025_09_26_plot_contours_IA_prior.ipynb | 258 +++++++++ .../2025_09_26_plot_contours_bias_prior.ipynb | 258 +++++++++ .../2025_09_26_plot_contours_cl_vs_xi.ipynb | 215 +++++--- .../2025_09_26_plot_contours_unblind.ipynb | 256 +++++++++ .../harmonic/2025_10_08_plot_data_vectors.py | 124 ++++- papers/harmonic/2025_10_28_plot_whisker.ipynb | 516 +++++++++--------- 8 files changed, 1522 insertions(+), 420 deletions(-) create mode 100644 papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb create mode 100644 papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb create mode 100644 papers/harmonic/2025_09_26_plot_contours_unblind.ipynb diff --git a/papers/harmonic/2025_09_11_psf_leakage_cell.py b/papers/harmonic/2025_09_11_psf_leakage_cell.py index 45f133620..74631ad48 100644 --- a/papers/harmonic/2025_09_11_psf_leakage_cell.py +++ b/papers/harmonic/2025_09_11_psf_leakage_cell.py @@ -23,8 +23,8 @@ from astropy.io import fits from matplotlib import scale as mscale from mpl_toolkits.axes_grid1 import make_axes_locatable - -from sp_validation.rho_tau import SquareRootScale +from sp_validation.utils_cosmo_val import SquareRootScale +from utils import get_chi2_and_pte mscale.register_scale(SquareRootScale) @@ -103,7 +103,7 @@ def get_shape_noise(e1, e2, w): f_psf = nmt.NmtFieldCatalog( positions=[cat_star["RA"], cat_star["DEC"]], weights=np.ones_like(cat_star["RA"]), - field=[cat_star["HSM_G1_PSF"], -cat_star["HSM_G2_PSF"]], + field=[cat_star["E1_PSF_HSM"], -cat_star["E2_PSF_HSM"]], lmax=b_lmax, lmax_mask=b_lmax, spin=2, @@ -235,7 +235,7 @@ def compute_iNKA_covariance(field1, field2, coupled=False): # %% -field_psf_map = get_field_catalog(cat_star, "HSM_G1_PSF", "HSM_G2_PSF", nside=1024) +field_psf_map = get_field_catalog(cat_star, "E1_PSF_HSM", "E2_PSF_HSM", nside=1024) field_gal_map = get_field_catalog(cat_gal, "e1", "e2", nside=1024, w_col="w_des") field_gal_corrected_map = get_field_catalog( cat_gal, "e1_leak_corrected", "e2_leak_corrected", nside=1024, w_col="w_des" @@ -370,80 +370,159 @@ def cov_to_corr(cov): plt.savefig("rho_0_cl.png", dpi=300) plt.show() + + +# %% +def sci_notation(x, sig=2): + if x == 0: + return "0" + exponent = int(np.floor(np.log10(abs(x)))) + mantissa = x / 10**exponent + return rf"{mantissa:.{sig}f} \times 10^{{{exponent}}}" + + # %% # Paper plot tau_0 -plt.figure(figsize=(10, 6)) +fig, (ax1, ax2) = plt.subplots( + 2, + 1, + figsize=(10, 6), + sharex=True, + gridspec_kw={"hspace": 0.0}, # tighten vertical space between panels +) -offset = 0.15 +offset = 0.2 ell_widths = np.diff(ell_eff) ell_widths = np.append(ell_widths, ell_widths[-1]) list_offset = [] -for i in range(4): - # Better jittering: symmetric around original ell values - jitter_fraction = (i - (4 - 1) / 2) * offset +for i in range(2): + jitter_fraction = (i - (2 - 1) / 2) * offset jittered_ell = ell_eff + jitter_fraction * ell_widths list_offset.append(jittered_ell) - -plt.errorbar( +# --- Top panel --- +chi2, reduced_chi2, pte = get_chi2_and_pte(tau_cl[0], cov_tau_0_ee, verbose=True) +chi2_corr, reduced_chi2_corr, pte_corr = get_chi2_and_pte( + tau_cl_corrected[0], cov_tau_0_ee, verbose=True +) +ax1.errorbar( list_offset[0], ell_eff * tau_cl[0], yerr=ell_eff * np.sqrt(cov_tau_0_ee.diagonal()), label=r"$\tau_0$ EE", fmt="o", capsize=2, + color="C1", + markersize=5, ) -plt.errorbar( +ax1.errorbar( list_offset[1], ell_eff * tau_cl_corrected[0], yerr=ell_eff * np.sqrt(cov_tau_0_ee.diagonal()), label=r"$\tau_0$ corrected EE", - fmt="o", + fmt="v", capsize=2, + color="C7", + markersize=5, ) -plt.errorbar( - list_offset[2], + +ax1.set_xscale("squareroot") +ax1.set_xticks(np.array([100, 400, 900, 1600])) +ax1.minorticks_on() +ax1.tick_params(axis="x", which="minor", length=2, width=0.8) +minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)] +ax1.set_xticks(minor_ticks, minor=True) + +ax1.axhline(0, c="k", ls="--", alpha=0.7) +ax1.tick_params(axis="both", which="major", labelsize=16) +ax1.yaxis.get_offset_text().set_visible(False) + +ax1.text( + 400, + -2.6e-7, + rf"$\chi^2_\nu={reduced_chi2:.3f}$, ${{\rm PTE}}={sci_notation(pte)}$", + color="C1", + fontsize=18, +) +ax1.text( + 400, + -2.6e-7 - 1.5e-7, + rf"$\chi^2_\nu={reduced_chi2_corr:.3f}$, ${{\rm PTE}}={sci_notation(pte_corr)}$", + color="C7", + fontsize=18, +) + +ax1.set_ylabel(r"$\ell C^{\tau_0}_\ell \times 10^{-7}$", fontsize=20) +ax1.set_xlim(1, 2048) +ax1.set_ylim(-5e-7, 5e-7) +ax1.legend(fontsize=16) + +# hide the shared x tick labels on the top panel +ax1.tick_params(axis="x", which="both", labelbottom=False) + +# --- Bottom panel --- +chi2, reduced_chi2, pte = get_chi2_and_pte(tau_cl[3], cov_tau_0_bb, verbose=True) +chi2_corr, reduced_chi2_corr, pte_corr = get_chi2_and_pte( + tau_cl_corrected[3], cov_tau_0_bb, verbose=True +) +ax2.errorbar( + list_offset[0], ell_eff * tau_cl[3], yerr=ell_eff * np.sqrt(cov_tau_0_bb.diagonal()), label=r"$\tau_0$ BB", fmt="o", capsize=2, + color="C2", + markersize=5, ) -plt.errorbar( - list_offset[3], +ax2.errorbar( + list_offset[1], ell_eff * tau_cl_corrected[3], yerr=ell_eff * np.sqrt(cov_tau_0_bb.diagonal()), label=r"$\tau_0$ corrected BB", - fmt="o", + fmt="v", capsize=2, + color="C3", + markersize=5, ) -plt.xscale("squareroot") -plt.xticks(np.array([100, 400, 900, 1600])) -plt.minorticks_on() -plt.tick_params(axis="x", which="minor", length=2, width=0.8) -minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)] -plt.xticks(minor_ticks, minor=True) +ax2.set_xscale("squareroot") +ax2.set_xticks(np.array([100, 400, 900, 1600])) +ax2.minorticks_on() +ax2.tick_params(axis="x", which="minor", length=2, width=0.8) +ax2.set_xticks(minor_ticks, minor=True) -plt.axhline(0, c="k", ls="--", alpha=0.7) +ax2.axhline(0, c="k", ls="--", alpha=0.7) +ax2.tick_params(axis="both", which="major", labelsize=16) +ax2.yaxis.get_offset_text().set_visible(False) -plt.tick_params(axis="both", which="major", labelsize=16) -plt.gca().yaxis.get_offset_text().set_visible(False) +ax2.text( + 400, + -2.6e-7, + rf"$\chi^2_\nu={reduced_chi2:.3f}$, ${{\rm PTE}}={sci_notation(pte)}$", + color="C2", + fontsize=18, +) +ax2.text( + 400, + -2.6e-7 - 1.5e-7, + rf"$\chi^2_\nu={reduced_chi2_corr:.3f}$, ${{\rm PTE}}={sci_notation(pte_corr)}$", + color="C3", + fontsize=18, +) +ax2.set_xlabel(r"Multipole $\ell$", fontsize=20) +ax2.set_ylabel(r"$\ell C^{\tau_0}_\ell \times 10^{-7}$", fontsize=20) +ax2.set_xlim(1, 2048) +ax2.set_ylim(-5e-7, 5e-7) +ax2.legend(fontsize=16) -plt.xlabel(r"Multipole $\ell$", fontsize=20) -plt.ylabel(r"$\ell C^{\tau_0}_\ell \times 10^{-7}$", fontsize=20) -plt.xlim(1, 2048) -plt.ylim(-5e-7, 5e-7) -plt.legend(fontsize=16) plt.savefig("./plots/tau_0_cl.png", dpi=300) -# Save pdf plt.savefig("./plots/tau_0_cl.pdf") - plt.show() # %% @@ -631,6 +710,7 @@ def cov_to_corr(cov): leakage_bias_corrected_std = np.std(leakage_bias_corrected_samples, axis=0) # %% +# Paper plot ratio w.r.t the data vector. plt.figure() leakage_bias = tau_cl[0] ** 2 / rho_cl[0] @@ -642,17 +722,19 @@ def cov_to_corr(cov): yerr=leakage_bias_std / cell_cl[0], label=r"Without object-wise leakage correction", fmt="o", - markersize=3, + markersize=5, capsize=2, + color="C1", ) plt.errorbar( list_offset[1], leakage_bias_corrected / cell_cl_corrected[0], yerr=leakage_bias_corrected_std / cell_cl_corrected[0], label="With object-wise leakage correction", - fmt="o", - markersize=3, + fmt="v", + markersize=5, capsize=2, + color="C7", ) threshold = 0.05 @@ -685,6 +767,72 @@ def cov_to_corr(cov): plt.savefig("./plots/leakage_bias_fraction_ee.pdf") plt.show() +# %% +# Paper plot +# Plot the ratio w.r.t the error bars + +# Load the covariance +path_cov = "/home/guerrini/sp_validation/notebooks/cosmo_val/output/pseudo_cl_cov_g_ng_iNKA_SP_v1.4.6.3_B.fits" +path_cov_corr = "/home/guerrini/sp_validation/notebooks/cosmo_val/output/pseudo_cl_cov_g_ng_iNKA_SP_v1.4.6.3_leak_corr_B.fits" +cov_cl = fits.open(path_cov)["COVAR_FULL"].data +cov_cl_corrected = fits.open(path_cov_corr)["COVAR_FULL"].data + +plt.figure() + +leakage_bias = tau_cl[0] ** 2 / rho_cl[0] +leakage_bias_corrected = tau_cl_corrected[0] ** 2 / rho_cl[0] + +plt.errorbar( + list_offset[0], + leakage_bias / np.sqrt(cov_cl.diagonal()), + yerr=leakage_bias_std / np.sqrt(cov_cl.diagonal()), + label=r"Without object-wise leakage correction", + fmt="o", + markersize=5, + capsize=2, + color="C1", +) +plt.errorbar( + list_offset[1], + leakage_bias_corrected / np.sqrt(cov_cl_corrected.diagonal()), + yerr=leakage_bias_corrected_std / np.sqrt(cov_cl_corrected.diagonal()), + label="With object-wise leakage correction", + fmt="v", + markersize=5, + capsize=2, + color="C7", +) + +threshold = 0.20 +plt.fill_between( + np.arange(7, 2050), + -threshold, + threshold, + color="lightgray", + alpha=0.5, + label=rf"{int(threshold * 100)}\% threshold", +) +plt.xscale("squareroot") + +plt.axhline(threshold, color="black", linestyle="--", alpha=0.6) +plt.axhline(-threshold, color="black", linestyle="--", alpha=0.6) +plt.xticks(np.array([100, 400, 900, 1600])) +plt.minorticks_on() +plt.tick_params(axis="x", which="minor", length=2, width=0.8) +minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)] +plt.xticks(minor_ticks, minor=True) +plt.xlim(ell_eff[0], ell_eff[-1]) +plt.xlabel(r"Multipole $\ell$", fontsize=16) +plt.xlim(7, 2050) +plt.ylabel(r"$C_\ell^{\rm sys} / \sigma(C_\ell)$", fontsize=16) +plt.ylim(-0.50, 0.50) +plt.gca().yaxis.set_major_formatter(mticker.PercentFormatter(xmax=1)) +plt.legend(loc="lower right", fontsize=12) +plt.savefig("./plots/leakage_bias_fraction_wrt_cov_ee.png", dpi=300) +# Save PDF +plt.savefig("./plots/leakage_bias_fraction_wrt_cov_ee.pdf") +plt.show() + # %% tau_samples_bb = np.random.multivariate_normal( mean=tau_cl[3], cov=cov_tau_0_bb, size=10000 @@ -787,10 +935,10 @@ def cov_to_corr(cov): # %% # Plot signal with errorbars pseudo_cl = fits.getdata( - "/home/guerrini/sp_validation/cosmo_val/output/pseudo_cl_SP_v1.4.5_leak_corr.fits" + "/home/guerrini/sp_validation/notebooks/cosmo_val/output/pseudo_cl_SP_v1.4.5_leak_corr.fits" ) cov_cl = fits.open( - "/home/guerrini/sp_validation/cosmo_val/output/pseudo_cl_cov_SP_v1.4.5_leak_corr.fits" + "/home/guerrini/sp_validation/notebooks/cosmo_val/output/pseudo_cl_cov_SP_v1.4.5_leak_corr.fits" ) plt.figure() diff --git a/papers/harmonic/2025_09_11_validation_cov_cell.py b/papers/harmonic/2025_09_11_validation_cov_cell.py index f235b4d78..af2367e2e 100644 --- a/papers/harmonic/2025_09_11_validation_cov_cell.py +++ b/papers/harmonic/2025_09_11_validation_cov_cell.py @@ -21,10 +21,9 @@ from matplotlib import scale as mscale from matplotlib.gridspec import GridSpec from mpl_toolkits.axes_grid1 import make_axes_locatable +from sp_validation.utils_cosmo_val import SquareRootScale from tqdm import tqdm -from sp_validation.rho_tau import SquareRootScale - mscale.register_scale(SquareRootScale) plt.style.use("./matplotlib_config/paper.mplstyle") @@ -36,8 +35,8 @@ # %% blind = "B" -path_cov_namaster = f"/home/guerrini/sp_validation/cosmo_val/output/pseudo_cl_cov_SP_v1.4.6.3_{blind}_leak_corr.fits" -path_glass_sims_output = "/n09data/guerrini/glass_mock_v1.4.6.3_v2/results/" +path_cov_namaster = f"/home/guerrini/sp_validation/notebooks/cosmo_val/output/pseudo_cl_cov_SP_v1.4.6.3_{blind}_leak_corr.fits" +path_glass_sims_output = "/n09data/guerrini/glass_mocks_v1.4.6.3_v4/results/" cov_namaster = fits.open(path_cov_namaster) @@ -86,7 +85,7 @@ def get_cov_from_one_cov(cov_one_cov, gaussian=True): # %% cov_one_cov = np.genfromtxt( - f"/home/guerrini/sp_validation/cosmo_val/output/pseudo_cl_cov_onecov_SP_v1.4.6.3_{blind}_leak_corr/covariance_list_3x2pt_pure_Cell.dat" + f"/home/guerrini/sp_validation/notebooks/cosmo_val/output/pseudo_cl_cov_onecov_SP_v1.4.6.3_{blind}_leak_corr/covariance_list_3x2pt_pure_Cell.dat" ) gaussian_one_cov = get_cov_from_one_cov(cov_one_cov, gaussian=True) @@ -489,7 +488,7 @@ def get_cov_from_one_cov(cov_one_cov, gaussian=True): # Set ticks and scales for the y-axis ax2.set_yticks([-0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3]) -ax2.set_ylim(-0.3, 0.3) +ax2.set_ylim(-0.2, 0.1) ax2.set_ylabel(r"$ \Delta \hat{\sigma} / \sigma$") ax2.set_xlabel(r"Multipole $\ell$") @@ -502,9 +501,7 @@ def get_cov_from_one_cov(cov_one_cov, gaussian=True): # Validation against Gaussian simulations -path_gaussian_sims = ( - "/n17data/sguerrini/sp_validation/cosmo_val/harmonic_covariance_gaussian_sims_v1463" -) +path_gaussian_sims = "/n17data/sguerrini/sp_validation/notebooks/cosmo_val/harmonic_covariance_gaussian_sims_v1463" n_sims = 10_000 covs = {} @@ -872,3 +869,77 @@ def get_cov_from_one_cov(cov_one_cov, gaussian=True): plt.show() # %% +# Save gaussian simulation covariance for bb +cov_gaussian_bb = cov_sim_gaussian[96:, 96:] +np.save( + "/home/guerrini/sp_validation/notebooks/cosmo_val/harmonic_covariance_gaussian_sims/cov_gaussian_sims_BB.npy", + cov_gaussian_bb, +) +cov_gaussian_eb = cov_sim_gaussian[64:96, 64:96] +np.save( + "/home/guerrini/sp_validation/notebooks/cosmo_val/harmonic_covariance_gaussian_sims/cov_gaussian_sims_EB.npy", + cov_gaussian_eb, +) + +# %% +# Compare map-based and catalogue-based Cl covariance +n_sims = 350 +cls_all_map = np.array([]).reshape((0, 32)) +for i in tqdm(range(n_sims)): + cls_map = np.load( + f"{path_glass_sims_output}/cl_map_glass_mock_{str(i + 1).zfill(5)}_4096.npy" + ) + cls_all_map = np.vstack((cls_all_map, cls_map[1])) + +cov_sim_map = np.cov(cls_all_map.T) + +# %% +fig = plt.figure() + +gs = GridSpec(2, 1, height_ratios=[3, 1], hspace=0.0) + +ax1 = fig.add_subplot(gs[0]) +ax2 = fig.add_subplot(gs[1], sharex=ax1) + +ax1.plot(ell, np.sqrt(np.diag(cov_sim)), label="Cov sim") +ax1.plot(ell, np.sqrt(np.diag(cov_sim_map)), label="Cov sim map") + +# Set ticks and scales for the x-axis +ax1.set_xscale("squareroot") +ax1.set_xticks(np.array([100, 400, 900, 1600])) +ax1.minorticks_on() +ax1.tick_params(axis="x", which="minor", length=2, width=0.8) +minor_ticks = [i * 10 for i in range(1, 10)] + [i * 100 for i in range(1, 21)] +ax1.set_xticks(minor_ticks, minor=True) + +# Set ticks and scale for the y-axis +ax1.set_yscale("log") +ax1.set_yticks([1e-11, 1e-11, 1e-10, 1e-9]) +ax1.tick_params(axis="y", which="minor", length=2, width=0.8) +minor_ticks = ( + [i * 1e-11 for i in range(1, 10)] + + [i * 1e-10 for i in range(1, 10)] + + [i * 1e-9 for i in range(1, 8)] +) +ax1.set_yticks(minor_ticks, minor=True) + +# Set labels and legend +ax1.set_ylabel(r"$\sigma(C_\ell^{EE})$") +ax1.legend(fontsize=10) + +baseline = np.sqrt(np.diag(cov_sim)) + +relative_error = np.sqrt(np.diag(cov_sim)) / baseline +ax2.plot(ell, relative_error - 1, label="cov sim") + +relative_error_map = np.sqrt(np.diag(cov_sim_map)) / baseline +ax2.plot(ell, relative_error_map - 1, label="cov sim map") + +# Set ticks and scales for the y-axis +ax2.set_yticks([-0.1, -0.05, 0.0, 0.05, 0.1]) +ax2.set_ylim(-0.1, 0.1) +ax2.set_ylabel(r"$ \Delta \hat{\sigma} / \sigma$") +ax2.set_xlabel(r"Multipole $\ell$") + +plt.show() +# %% diff --git a/papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb b/papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb new file mode 100644 index 000000000..1a6a2fdfa --- /dev/null +++ b/papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb @@ -0,0 +1,258 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "0", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "# Trick to plot with tex\n", + "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", + "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", + "\n", + "from getdist import plots, loadMCSamples\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "\n", + "plt.style.use(\n", + " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", + ")\n", + "\n", + "sns.set_palette(\"husl\")\n", + "\n", + "g = plots.get_subplot_plotter(width_inch=30)\n", + "g.settings.axes_fontsize=60\n", + "g.settings.axes_labelsize=60\n", + "g.settings.alpha_filled_add = 0.7\n", + "g.settings.legend_fontsize = 60\n", + "\n", + "%matplotlib inline\n", + "\n", + "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "root_dir = \"/n09data/guerrini/output_chains/\"\n", + "\n", + "roots = [\n", + " \"SP_v1.4.6.3_leak_corr_B\",\n", + " \"SP_v1.4.6.3_leak_corr_unif_A_IA_B\",\n", + " \"SP_v1.4.6.3_leak_corr_no_IA_B\",\n", + "]\n", + "\n", + "print(roots)\n", + "\n", + "os.makedirs(\"plots\", exist_ok=True)\n", + "os.makedirs(\"plots/cosmological_constraints_ia_prior\", exist_ok=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "# MAKE PARAMNAMES FILE\n", + "\n", + "for root in roots:\n", + " with open(root_dir + '{}/{}/samples_{}_cell.txt'.format('/'+root ,root, root), \"r\") as file:\n", + " params = file.readline()[1:].split('\\t')[:-4]\n", + " file.close()\n", + " \n", + " with open(root_dir + '{}/{}/getdist_{}_cell.paramnames'.format('/'+root, root,root), \"w\") as file:\n", + " for i in range(len(params)):\n", + " if len(params[i].split('--')) > 1:\n", + " file.write(params[i].split('--')[1] + '\\n')\n", + " else:\n", + " file.write(params[i].split('--')[0] + '\\n')\n", + " file.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "#READ CHAIN\n", + "\n", + "chains=[]\n", + "\n", + "for root in roots:\n", + "\n", + " samples = np.loadtxt(root_dir + '{}/{}/samples_{}_cell.txt'.format(root,root,root))\n", + " print(len(samples))\n", + " if 'nautilus' in root:\n", + " samples = np.column_stack((np.exp(samples[:,-3]),samples[:,-1]-samples[:,-2],samples[:,0:-3]))\n", + " else:\n", + " samples = np.column_stack((samples[:,-1],samples[:,-3],samples[:,0:-4]))\n", + " np.savetxt(root_dir + '{}/{}/getdist_{}_cell.txt'.format(root,root,root), samples)\n", + " \n", + " chain = g.samples_for_root(root_dir + '{}/{}/getdist_{}_cell'.format(root,root,root),\n", + " cache=False,\n", + " settings={'ignore_rows':0,\n", + " 'smooth_scale_2D':0.5,\n", + " 'smooth_scale_1D':0.5})\n", + "\n", + " chains.append(chain)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", + "label_list = label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "\n", + "for chain in chains:\n", + " param_names = chain.getParamNames()\n", + " for name, label in zip(name_list, label_list):\n", + " param_names.parWithName(name).label = label" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "legend_labels = [\n", + " r\"UNIONS $C_\\ell$ (This work)\",\n", + " r\"Uniform $A_{\\rm IA}$\",\n", + " r\"No $A_{\\rm IA}$\",\n", + "]\n", + "\n", + "colours = [\n", + " \"crimson\",\n", + " \"royalblue\",\n", + " \"forestgreen\"\n", + "]\n", + "\n", + "linestyle = [\n", + " 'solid',\n", + " 'solid',\n", + " 'solid'\n", + "]\n", + "\n", + "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", + "\n", + "#Plot all parameters\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " legend_labels=legend_labels,\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " legend_loc='upper right',\n", + " filled=True\n", + ")\n", + "\n", + "g.export('plots/cosmological_constraints_ia_prior/contours_all_cell.png')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "#Plot only cosmological parameters\n", + "\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','s_8_input', 'SIGMA_8', 'a'],\n", + " legend_labels=legend_labels,\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " legend_loc='upper right',\n", + " filled=True\n", + ")\n", + "\n", + "g.export('plots/cosmological_constraints_ia_prior/contours_cosmo_cell.png')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "#Plot S8 Omega_m only\n", + "#Plot S8 Omega_m only\n", + "g = plots.get_subplot_plotter(width_inch=20)\n", + "g.settings.axes_fontsize=35\n", + "g.settings.axes_labelsize=50\n", + "g.settings.alpha_filled_add = 0.7\n", + "g.settings.legend_fontsize = 40\n", + "\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','s_8_input'],\n", + " legend_labels=legend_labels,\n", + " legend_loc='upper right',\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " filled=True,\n", + ")\n", + "\n", + "\"\"\" plt.figtext(0.5, 0.5, 'PRELIMINARY',\n", + " fontsize=150, color='gray',\n", + " ha='center', va='center',\n", + " alpha=0.3, rotation=330) \"\"\"\n", + "\n", + "g.export('plots/cosmological_constraints_ia_prior/contours_s8_omegam_cell.png')\n", + "#Save PDF\n", + "g.export('plots/cosmological_constraints_ia_prior/contours_s8_omegam_cell.pdf')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "sp_validation_3.11", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb b/papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb new file mode 100644 index 000000000..8df25bed3 --- /dev/null +++ b/papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb @@ -0,0 +1,258 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "0", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "# Trick to plot with tex\n", + "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", + "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", + "\n", + "from getdist import plots, loadMCSamples\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "\n", + "plt.style.use(\n", + " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", + ")\n", + "\n", + "sns.set_palette(\"husl\")\n", + "\n", + "g = plots.get_subplot_plotter(width_inch=30)\n", + "g.settings.axes_fontsize=60\n", + "g.settings.axes_labelsize=60\n", + "g.settings.alpha_filled_add = 0.7\n", + "g.settings.legend_fontsize = 60\n", + "\n", + "%matplotlib inline\n", + "\n", + "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "root_dir = \"/n09data/guerrini/output_chains/\"\n", + "\n", + "roots = [\n", + " \"SP_v1.4.6.3_leak_corr_unif_dz_B\",\n", + " \"SP_v1.4.6.3_leak_corr_no_dz_B\",\n", + " \"SP_v1.4.6.3_leak_corr_no_m_bias_B\",\n", + "]\n", + "\n", + "print(roots)\n", + "\n", + "os.makedirs(\"plots\", exist_ok=True)\n", + "os.makedirs(\"plots/cosmological_constraints_bias_prior\", exist_ok=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "# MAKE PARAMNAMES FILE\n", + "\n", + "for root in roots:\n", + " with open(root_dir + '{}/{}/samples_{}_cell.txt'.format('/'+root ,root, root), \"r\") as file:\n", + " params = file.readline()[1:].split('\\t')[:-4]\n", + " file.close()\n", + " \n", + " with open(root_dir + '{}/{}/getdist_{}_cell.paramnames'.format('/'+root, root,root), \"w\") as file:\n", + " for i in range(len(params)):\n", + " if len(params[i].split('--')) > 1:\n", + " file.write(params[i].split('--')[1] + '\\n')\n", + " else:\n", + " file.write(params[i].split('--')[0] + '\\n')\n", + " file.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "#READ CHAIN\n", + "\n", + "chains=[]\n", + "\n", + "for root in roots:\n", + "\n", + " samples = np.loadtxt(root_dir + '{}/{}/samples_{}_cell.txt'.format(root,root,root))\n", + " print(len(samples))\n", + " if 'nautilus' in root:\n", + " samples = np.column_stack((np.exp(samples[:,-3]),samples[:,-1]-samples[:,-2],samples[:,0:-3]))\n", + " else:\n", + " samples = np.column_stack((samples[:,-1],samples[:,-3],samples[:,0:-4]))\n", + " np.savetxt(root_dir + '{}/{}/getdist_{}_cell.txt'.format(root,root,root), samples)\n", + " \n", + " chain = g.samples_for_root(root_dir + '{}/{}/getdist_{}_cell'.format(root,root,root),\n", + " cache=False,\n", + " settings={'ignore_rows':0,\n", + " 'smooth_scale_2D':0.5,\n", + " 'smooth_scale_1D':0.5})\n", + "\n", + " chains.append(chain)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", + "label_list = label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "\n", + "for chain in chains:\n", + " param_names = chain.getParamNames()\n", + " for name, label in zip(name_list, label_list):\n", + " param_names.parWithName(name).label = label" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "legend_labels = [\n", + " r\"Uniform $\\Delta z$ prior\",\n", + " r\"No $\\Delta z$ prior\",\n", + " r\"No $m$ bias prior\"\n", + "]\n", + "\n", + "colours = [\n", + " \"royalblue\",\n", + " \"crimson\",\n", + " \"forestgreen\"\n", + "]\n", + "\n", + "linestyle = [\n", + " 'solid',\n", + " 'solid',\n", + " 'solid'\n", + "]\n", + "\n", + "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", + "\n", + "#Plot all parameters\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " legend_labels=legend_labels,\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " legend_loc='upper right',\n", + " filled=True\n", + ")\n", + "\n", + "g.export('plots/cosmological_constraints_bias_prior/contours_all_cell.png')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "#Plot only cosmological parameters\n", + "\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','s_8_input', 'SIGMA_8', 'a'],\n", + " legend_labels=legend_labels,\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " legend_loc='upper right',\n", + " filled=True\n", + ")\n", + "\n", + "g.export('plots/cosmological_constraints_bias_prior/contours_cosmo_cell.png')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "#Plot S8 Omega_m only\n", + "#Plot S8 Omega_m only\n", + "g = plots.get_subplot_plotter(width_inch=20)\n", + "g.settings.axes_fontsize=35\n", + "g.settings.axes_labelsize=50\n", + "g.settings.alpha_filled_add = 0.7\n", + "g.settings.legend_fontsize = 40\n", + "\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','s_8_input'],\n", + " legend_labels=legend_labels,\n", + " legend_loc='upper right',\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " filled=True,\n", + ")\n", + "\n", + "\"\"\" plt.figtext(0.5, 0.5, 'PRELIMINARY',\n", + " fontsize=150, color='gray',\n", + " ha='center', va='center',\n", + " alpha=0.3, rotation=330) \"\"\"\n", + "\n", + "g.export('plots/cosmological_constraints_bias_prior/contours_s8_omegam_cell.png')\n", + "#Save PDF\n", + "g.export('plots/cosmological_constraints_bias_prior/contours_s8_omegam_cell.pdf')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "sp_validation_3.11", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb b/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb index ca2ec8f55..35c0180e7 100644 --- a/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb +++ b/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb @@ -16,42 +16,70 @@ "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", "\n", "# Append any useful folder in the path\n", - "sys.path.append(\"/home/guerrini/sp_validation/cosmo_inference/scripts/\")\n", + "sys.path.append(\n", + " \"/home/guerrini/sp_validation/cosmo_inference/scripts/\"\n", + ")\n", "\n", - "import chain_postprocessing as cp\n", + "from getdist import plots, loadMCSamples\n", + "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", - "from getdist import plots\n", "\n", - "plt.style.use(\"/home/guerrini/matplotlib_config/paper.mplstyle\")\n", + "import chain_postprocessing as cp\n", + "\n", + "plt.style.use(\n", + " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", + ")\n", "\n", "sns.set_palette(\"husl\")\n", "\n", "g = plots.get_subplot_plotter(width_inch=20)\n", - "g.settings.axes_fontsize = 60\n", - "g.settings.axes_labelsize = 60\n", + "g.settings.axes_fontsize=60\n", + "g.settings.axes_labelsize=60\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 50\n", "\n", "%matplotlib inline\n", "\n", - "# SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", "root_dir = \"/n09data/guerrini/output_chains/\"\n", "root_external = \"/n09data/guerrini/output_chains/ext_data/\"\n", "\n", "blind = \"B\"\n", "\n", - "colour_blind = {\"A\": \"royalblue\", \"B\": \"crimson\", \"C\": \"forestgreen\"}\n", + "colour_blind = {\n", + " \"A\": \"royalblue\",\n", + " \"B\": \"crimson\",\n", + " \"C\": \"forestgreen\"\n", + "}\n", "\n", "roots = [\n", + " f\"SP_v1.4.6.3_leak_corr_{blind}_prior_volume\",\n", " f\"SP_v1.4.6.3_leak_corr_{blind}\",\n", " f\"SP_v1.4.6.3_{blind}_fiducial_config\",\n", " \"Planck18\",\n", "]\n", "\n", - "categories = [\"harmonic\", \"configuration\", \"external\"]\n", + "categories = [\n", + " \"harmonic\",\n", + " \"harmonic\",\n", + " \"configuration\",\n", + " \"external\"\n", + "]\n", + "\n", + "need_compute_getdist = [\n", + " True,\n", + " True,\n", + " True,\n", + " False\n", + "]\n", "\n", - "need_compute_getdist = [True, True, False]\n", + "chain_types = [\n", + " \"apriori\",\n", + " \"polychord\",\n", + " \"polychord\",\n", + " \"polychord\",\n", + "]\n", "\n", "print(roots)\n", "\n", @@ -70,22 +98,48 @@ "for i, root in enumerate(roots):\n", " analysis_type_i = categories[i]\n", " if analysis_type_i == \"harmonic\":\n", - " path_samples = os.path.join(root_dir, root, root, f\"samples_{root}_cell.txt\")\n", - " path_gd = os.path.join(root_dir, root, root, f\"getdist_{root}_cell\")\n", + " path_samples = os.path.join(\n", + " root_dir,\n", + " root,\n", + " root,\n", + " f\"samples_{root}_cell.txt\"\n", + " )\n", + " path_gd = os.path.join(\n", + " root_dir,\n", + " root,\n", + " root,\n", + " f\"getdist_{root}_cell\"\n", + " )\n", " elif analysis_type_i == \"configuration\":\n", - " path_samples = os.path.join(root_dir, root, f\"samples_{root}.txt\")\n", - " path_gd = os.path.join(root_dir, root, f\"getdist_{root}\")\n", + " path_samples = os.path.join(\n", + " root_dir,\n", + " root,\n", + " f\"samples_{root}.txt\"\n", + " )\n", + " path_gd = os.path.join(\n", + " root_dir,\n", + " root,\n", + " f\"getdist_{root}\"\n", + " )\n", " elif analysis_type_i == \"external\":\n", - " path_samples = os.path.join(root_external, root, f\"samples_{root}.txt\")\n", - " path_gd = os.path.join(root_external, root, f\"getdist_{root}\")\n", + " path_samples = os.path.join(\n", + " root_external,\n", + " root,\n", + " f\"samples_{root}.txt\"\n", + " )\n", + " path_gd = os.path.join(\n", + " root_external,\n", + " root,\n", + " f\"getdist_{root}\"\n", + " )\n", " else:\n", " raise ValueError(f\"Invalid analysis type: {analysis_type_i}\")\n", - "\n", + " \n", " need_getdist = need_compute_getdist[i]\n", " if need_getdist:\n", " # Read files, format and load chains\n", - " cp.load_samples_and_write_paramnames(path_samples, path_gd + \".paramnames\")\n", - " cp.write_samples_getdist_format(path_samples, path_gd + \".txt\")\n", + " cp.load_samples_and_write_paramnames(path_samples, path_gd+\".paramnames\", chain_type=chain_types[i])\n", + " cp.write_samples_getdist_format(path_samples, path_gd+\".txt\", chain_type=chain_types[i])\n", " chain = cp.load_chain(path_gd, smoothing_scale=0.5)\n", " chains.append(chain)" ] @@ -97,39 +151,15 @@ "metadata": {}, "outputs": [], "source": [ - "name_list = [\n", - " \"OMEGA_M\",\n", - " \"ombh2\",\n", - " \"h0\",\n", - " \"n_s\",\n", - " \"SIGMA_8\",\n", - " \"S_8\",\n", - " \"s_8_input\",\n", - " \"logt_agn\",\n", - " \"a\",\n", - " \"m1\",\n", - " \"bias_1\",\n", - "]\n", - "label_list = [\n", - " r\"\\Omega_{\\rm m}\",\n", - " r\"\\omega_{\\rm b} h^2\",\n", - " r\"h_0\",\n", - " r\"n_{\\rm s}\",\n", - " r\"\\sigma_8\",\n", - " r\"S_8\",\n", - " r\"S_8\",\n", - " r\"\\log T_{\\rm AGN}\",\n", - " r\"A_{\\rm IA}\",\n", - " r\"m_1\",\n", - " r\"\\Delta z_1\",\n", - "]\n", + "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8', 'S_8', 's_8_input', 'logt_agn','a','m1','bias_1']\n", + "label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", "\n", "for chain in chains:\n", " param_names = chain.getParamNames()\n", " for name, label in zip(name_list, label_list):\n", " try:\n", " param_names.parWithName(name).label = label\n", - " except Exception:\n", + " except:\n", " warnings.warn(f\"Parameter {name} not found in chain {roots[i]}.\")" ] }, @@ -141,48 +171,50 @@ "outputs": [], "source": [ "legend_labels = [\n", + " \"Prior volume\",\n", " r\"UNIONS $C_\\ell$ (This work)\",\n", " r\"UNIONS $\\xi_\\pm$ (Goh et al., 2026)\",\n", " r\"\\textit{Planck} 2018\",\n", "]\n", "\n", - "colours = [colour_blind[\"B\"], \"orange\", \"violet\"]\n", + "colours = [\n", + " \"lightgrey\",\n", + " colour_blind[\"B\"],\n", + " \"orange\",\n", + " \"olive\",\n", + "]\n", "\n", - "linestyle = [\"solid\", \"solid\", \"solid\"]\n", + "linestyle = [\n", + " 'solid',\n", + " 'solid',\n", + " 'solid',\n", + " 'solid'\n", + "]\n", "\n", "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", "\n", + "idx = [0, 1, 2]\n", "g = plots.get_subplot_plotter(width_inch=30)\n", - "g.settings.axes_fontsize = 20\n", - "g.settings.axes_labelsize = 30\n", + "g.settings.axes_fontsize=20\n", + "g.settings.axes_labelsize=30\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 50\n", "\n", - "# Plot all parameters\n", + "#Plot all parameters\n", "g.triangle_plot(\n", - " chains[:-1], # Plot everything but Planck\n", - " [\n", - " \"OMEGA_M\",\n", - " \"ombh2\",\n", - " \"h0\",\n", - " \"n_s\",\n", - " \"SIGMA_8\",\n", - " \"s_8_input\",\n", - " \"logt_agn\",\n", - " \"a\",\n", - " \"m1\",\n", - " \"bias_1\",\n", - " ],\n", - " legend_labels=legend_labels,\n", - " line_args=line_args,\n", - " contour_colors=colours,\n", - " legend_loc=\"upper right\",\n", - " filled=True,\n", + " [chains[i] for i in idx], # Plot everything but Planck\n", + " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " legend_labels=[legend_labels[i] for i in idx],\n", + " line_args=[line_args[i] for i in idx],\n", + " contour_colors=[colours[i] for i in idx],\n", + " legend_loc='upper right',\n", + "\n", + " filled=True\n", ")\n", "\n", - "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_all_cell.png\")\n", - "# Save pdf\n", - "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_all_cell.pdf\")\n", + "g.export('plots/cosmological_constraints_cl_vs_xi/contours_all_cell.png')\n", + "#Save pdf\n", + "g.export('plots/cosmological_constraints_cl_vs_xi/contours_all_cell.pdf')\n", "plt.show()" ] }, @@ -193,24 +225,24 @@ "metadata": {}, "outputs": [], "source": [ - "# Plot only cosmological parameters\n", + "#Plot only cosmological parameters\n", "g = plots.get_subplot_plotter(width_inch=20)\n", - "g.settings.axes_fontsize = 35\n", - "g.settings.axes_labelsize = 50\n", + "g.settings.axes_fontsize=35\n", + "g.settings.axes_labelsize=50\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 40\n", "\n", "g.triangle_plot(\n", " chains,\n", - " [\"OMEGA_M\", \"S_8\", \"SIGMA_8\", \"a\"],\n", + " ['OMEGA_M','S_8', 'SIGMA_8', 'a'],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc=\"upper right\",\n", - " filled=True,\n", + " legend_loc='upper right',\n", + " filled=True\n", ")\n", "\n", - "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_cosmo_cell.png\")\n", + "g.export('plots/cosmological_constraints_cl_vs_xi/contours_cosmo_cell.png')\n", "plt.show()" ] }, @@ -221,30 +253,37 @@ "metadata": {}, "outputs": [], "source": [ - "g.settings.axes_fontsize = 35\n", - "g.settings.axes_labelsize = 50\n", + "g.settings.axes_fontsize=35\n", + "g.settings.axes_labelsize=50\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 40\n", "\n", - "# Plot S8 Omega_m only\n", + "#Plot S8 Omega_m only\n", "g.triangle_plot(\n", - " chains,\n", - " [\"OMEGA_M\", \"S_8\"],\n", + " chains[:], # Plot everything but the prior volume\n", + " ['OMEGA_M','S_8'],\n", " legend_labels=legend_labels,\n", - " legend_loc=\"upper right\",\n", + " legend_loc='upper right',\n", " line_args=line_args,\n", " contour_colors=colours,\n", " filled=True,\n", ")\n", "\n", + "# Set S_8 range\n", + "s8_min, s8_max = 0.58, 1.15 # adjust to what you want\n", + "\n", + "g.subplots[1, 0].set_ylim(s8_min, s8_max) # 2D contour panel (S_8 on y-axis)\n", + "g.subplots[1, 0].set_xlim(0.05, 0.65)\n", + "g.subplots[1, 1].set_xlim(s8_min, s8_max) # S_8 marginal panel (S_8 on x-axis)\n", + "\n", "\"\"\" plt.figtext(0.5, 0.5, 'PRELIMINARY',\n", " fontsize=150, color='gray',\n", " ha='center', va='center',\n", " alpha=0.3, rotation=330) \"\"\"\n", "\n", - "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.png\")\n", - "# Save pdf\n", - "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.pdf\")\n", + "g.export('plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.png')\n", + "#Save pdf\n", + "g.export('plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.pdf')\n", "plt.show()" ] }, diff --git a/papers/harmonic/2025_09_26_plot_contours_unblind.ipynb b/papers/harmonic/2025_09_26_plot_contours_unblind.ipynb new file mode 100644 index 000000000..7a85ea81c --- /dev/null +++ b/papers/harmonic/2025_09_26_plot_contours_unblind.ipynb @@ -0,0 +1,256 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "0", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "# Trick to plot with tex\n", + "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", + "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", + "\n", + "from getdist import plots, loadMCSamples\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "\n", + "plt.style.use(\n", + " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", + ")\n", + "\n", + "sns.set_palette(\"husl\")\n", + "\n", + "g = plots.get_subplot_plotter(width_inch=30)\n", + "g.settings.axes_fontsize=60\n", + "g.settings.axes_labelsize=60\n", + "g.settings.alpha_filled_add = 0.7\n", + "g.settings.legend_fontsize = 60\n", + "\n", + "%matplotlib inline\n", + "\n", + "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "root_dir = \"/n09data/guerrini/output_chains/\"\n", + "\n", + "roots = [\n", + " \"SP_v1.4.6.3_leak_corr_unblind_A\",\n", + " \"SP_v1.4.6.3_leak_corr_unblind_B\",\n", + " \"SP_v1.4.6.3_leak_corr_unblind_C\",\n", + "]\n", + "\n", + "print(roots)\n", + "\n", + "os.makedirs(\"plots\", exist_ok=True)\n", + "os.makedirs(\"plots/cosmological_constraints_unblind\", exist_ok=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1", + "metadata": {}, + "outputs": [], + "source": [ + "# MAKE PARAMNAMES FILE\n", + "\n", + "for root in roots:\n", + " with open(root_dir + '{}/{}/samples_{}_cell.txt'.format('/'+root ,root, root), \"r\") as file:\n", + " params = file.readline()[1:].split('\\t')[:-4]\n", + " file.close()\n", + " \n", + " with open(root_dir + '{}/{}/getdist_{}_cell.paramnames'.format('/'+root, root,root), \"w\") as file:\n", + " for i in range(len(params)):\n", + " if len(params[i].split('--')) > 1:\n", + " file.write(params[i].split('--')[1] + '\\n')\n", + " else:\n", + " file.write(params[i].split('--')[0] + '\\n')\n", + " file.close()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2", + "metadata": {}, + "outputs": [], + "source": [ + "#READ CHAIN\n", + "\n", + "chains=[]\n", + "\n", + "for root in roots:\n", + "\n", + " samples = np.loadtxt(root_dir + '{}/{}/samples_{}_cell.txt'.format(root,root,root))\n", + " print(len(samples))\n", + " if 'nautilus' in root:\n", + " samples = np.column_stack((np.exp(samples[:,-3]),samples[:,-1]-samples[:,-2],samples[:,0:-3]))\n", + " else:\n", + " samples = np.column_stack((samples[:,-1],samples[:,-3],samples[:,0:-4]))\n", + " np.savetxt(root_dir + '{}/{}/getdist_{}_cell.txt'.format(root,root,root), samples)\n", + " \n", + " chain = g.samples_for_root(root_dir + '{}/{}/getdist_{}_cell'.format(root,root,root),\n", + " cache=False,\n", + " settings={'ignore_rows':0,\n", + " 'smooth_scale_2D':0.5,\n", + " 'smooth_scale_1D':0.5})\n", + "\n", + " chains.append(chain)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", + "label_list = label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "\n", + "for chain in chains:\n", + " param_names = chain.getParamNames()\n", + " for name, label in zip(name_list, label_list):\n", + " param_names.parWithName(name).label = label" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", + "metadata": {}, + "outputs": [], + "source": [ + "legend_labels = [\n", + " fr\"$\\langle \\Delta z \\rangle = -0.030$, Blind {bl}\" for bl in ['A', 'B', 'C']\n", + "]\n", + "\n", + "colours = [\n", + " \"royalblue\",\n", + " \"crimson\",\n", + " \"forestgreen\"\n", + "]\n", + "\n", + "linestyle = [\n", + " 'solid',\n", + " 'solid',\n", + " 'solid'\n", + "]\n", + "\n", + "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", + "\n", + "#Plot all parameters\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " legend_labels=legend_labels,\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " legend_loc='upper right',\n", + " filled=True\n", + ")\n", + "\n", + "g.export('plots/cosmological_constraints_unblind/contours_all_cell.png')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5", + "metadata": {}, + "outputs": [], + "source": [ + "#Plot only cosmological parameters\n", + "\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','s_8_input', 'SIGMA_8', 'a'],\n", + " legend_labels=legend_labels,\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " legend_loc='upper right',\n", + " filled=True\n", + ")\n", + "\n", + "g.export('plots/cosmological_constraints_unblind/contours_cosmo_cell.png')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6", + "metadata": {}, + "outputs": [], + "source": [ + "#Plot S8 Omega_m only\n", + "#Plot S8 Omega_m only\n", + "g = plots.get_subplot_plotter(width_inch=20)\n", + "g.settings.axes_fontsize=35\n", + "g.settings.axes_labelsize=50\n", + "g.settings.alpha_filled_add = 0.7\n", + "g.settings.legend_fontsize = 40\n", + "\n", + "g.triangle_plot(\n", + " chains,\n", + " ['OMEGA_M','s_8_input'],\n", + " legend_labels=legend_labels,\n", + " legend_loc='upper right',\n", + " line_args=line_args,\n", + " contour_colors=colours,\n", + " filled=True,\n", + ")\n", + "\n", + "\"\"\" plt.figtext(0.5, 0.5, 'PRELIMINARY',\n", + " fontsize=150, color='gray',\n", + " ha='center', va='center',\n", + " alpha=0.3, rotation=330) \"\"\"\n", + "\n", + "g.export('plots/cosmological_constraints_unblind/contours_s8_omegam_cell.png')\n", + "#Save PDF\n", + "g.export('plots/cosmological_constraints_unblind/contours_s8_omegam_cell.pdf')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "sp_validation_3.11", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/papers/harmonic/2025_10_08_plot_data_vectors.py b/papers/harmonic/2025_10_08_plot_data_vectors.py index 24edc97f1..385cc696a 100644 --- a/papers/harmonic/2025_10_08_plot_data_vectors.py +++ b/papers/harmonic/2025_10_08_plot_data_vectors.py @@ -20,10 +20,9 @@ import seaborn as sns from astropy.io import fits from matplotlib import scale as mscale +from sp_validation.utils_cosmo_val import SquareRootScale from utils import get_chi2_and_pte -from sp_validation.rho_tau import SquareRootScale - mscale.register_scale(SquareRootScale) plt.style.use("./matplotlib_config/paper.mplstyle") @@ -38,7 +37,7 @@ # %% -base_dir = "/home/guerrini/sp_validation/cosmo_val/output/" +base_dir = "/home/guerrini/sp_validation/notebooks/cosmo_val/output/" versions = ["SP_v1.4.5_leak_corr", "SP_v1.4.6_leak_corr", "SP_v1.4.8_leak_corr"] labels = [ @@ -194,24 +193,61 @@ plt.show() # %% +# Load the Gaussian sims covariance +from tqdm import tqdm + +path_gaussian_sims = "/n17data/sguerrini/sp_validation/notebooks/cosmo_val/harmonic_covariance_gaussian_sims_v1463" + +n_sims = 10_000 +covs = {} + +for blind in ["A", "B", "C"]: + path_folder = path_gaussian_sims + "_" + blind + cls_all_gaussian = np.array([]).reshape((0, 4, 32)) + cls_noise_gaussian = np.array([]).reshape((0, 4, 32)) + for i in tqdm(range(n_sims)): + try: + cls = np.load(f"{path_folder}/sample_{i}.npz") + cls_all_gaussian = np.vstack((cls_all_gaussian, cls["cl_all"][None, ...])) + cls_noise_gaussian = np.vstack( + (cls_noise_gaussian, cls["cl_noise"][None, ...]) + ) + except Exception as e: + print(f"Error loading {i}: {e}") + cls_all_gaussian = cls_all_gaussian - np.mean(cls_noise_gaussian, axis=0) + cls_all_gaussian = cls_all_gaussian.reshape((n_sims, -1)) + cov_gaussian = np.cov(cls_all_gaussian.T) + covs[blind] = {} + covs[blind]["cov_gaussian"] = cov_gaussian + covs[blind]["diag_EE"] = np.sqrt(np.diag(cov_gaussian[:32, :32])) + covs[blind]["diag_BB"] = np.sqrt(np.diag(cov_gaussian[96:, 96:])) + covs[blind]["diag_EB"] = np.sqrt(np.diag(cov_gaussian[64:96, 64:96])) +# %% # Same plot with EB and BB modes only of v1.4.6 fig, (ax1, ax2) = plt.subplots(ncols=1, nrows=2, figsize=(8, 6), sharex=True) -version_paper_plot = ["SP_v1.4.6.3_A", "SP_v1.4.6.3_leak_corr_A"] +version_paper_plot = ["SP_v1.4.6.3_B", "SP_v1.4.6.3_B_leak_corr"] labels = ["W/o leakage correction", "With leakage correction"] -colors = ["C0", "C1"] -markers = ["o", "h"] -offset = 0.2 +colors = ["C1", "C7"] +markers = ["o", "v"] +offset = 0.30 for i, ver in enumerate(version_paper_plot): cell = fits.getdata(f"{base_dir}/pseudo_cl_{ver}.fits") cov_cell = fits.open(f"{base_dir}/pseudo_cl_cov_{ver}.fits") ell = cell["ell"] + n_bins = len(ell) cl_eb = cell["EB"] cov_cl_eb = cov_cell["COVAR_EB_EB"].data + cov_cl_eb_sims = covs["B"]["cov_gaussian"][ + 2 * n_bins : 3 * n_bins, 2 * n_bins : 3 * n_bins + ] cl_bb = cell["BB"] cov_cl_bb = cov_cell["COVAR_BB_BB"].data + cov_cl_bb_sims = covs["B"]["cov_gaussian"][ + 3 * n_bins : 4 * n_bins, 3 * n_bins : 4 * n_bins + ] ell_widths = np.diff(ell) ell_widths = np.append(ell_widths, ell_widths[-1]) @@ -219,39 +255,66 @@ jitter_fraction = (i - (len(version_paper_plot) - 1) / 2) * offset jiterred_ell = ell + jitter_fraction * ell_widths - chi2, reduced_chi2, pte = get_chi2_and_pte(cl_eb, cov_cl_eb, verbose=False) + chi2_inka, reduced_chi2_inka, pte_inka = get_chi2_and_pte( + cl_eb, cov_cl_eb, verbose=False + ) + print( + f"Version: {ver}, EB iNKA chi2: {chi2_inka:.3f}, reduced chi2: {reduced_chi2_inka:.3f}, PTE: {pte_inka:.3f}" + ) + chi2_sim, reduced_chi2_sim, pte_sim = get_chi2_and_pte( + cl_eb, cov_cl_eb_sims, verbose=False + ) + print( + f"Version: {ver}, EB Gaussian sims chi2: {chi2_sim:.3f}, reduced chi2: {reduced_chi2_sim:.3f}, PTE: {pte_sim:.3f}" + ) ax1.errorbar( jiterred_ell, cl_eb * ell, - yerr=np.sqrt(np.diag(cov_cl_eb)) * ell, + yerr=np.sqrt(np.diag(cov_cl_eb_sims)) * ell, label=labels[i], color=colors[i], fmt=markers[i], capsize=2, + markersize=5, ) ax1.text( - 100, - -0.5e-7 - i * 0.25e-7, - rf"$\chi^2_\nu={reduced_chi2:.3f}$ | PTE={pte:.3f}", + 10, + -0.9e-7 - i * 0.28e-7, + rf"$\chi^2_{{\nu, \mathrm{{iNKA}}}}={reduced_chi2_inka:.3f}$, ${{\rm PTE}}_{{\rm iNKA}}={pte_inka:.3f}$ | $\chi^2_{{\nu, \textrm{{sim.}}}}={reduced_chi2_sim:.3f}$, ${{\rm PTE}}_{{\rm sim.}}={pte_sim:.3f}$", color=colors[i], + fontsize=15, + ) + + chi2_inka, reduced_chi2_inka, pte_inka = get_chi2_and_pte( + cl_bb, cov_cl_bb, verbose=False + ) + print( + f"Version: {ver}, BB iNKA chi2: {chi2_inka:.3f}, reduced chi2: {reduced_chi2_inka:.3f}, PTE: {pte_inka:.3f}" + ) + chi2_sim, reduced_chi2_sim, pte_sim = get_chi2_and_pte( + cl_bb, cov_cl_bb_sims, verbose=False + ) + print( + f"Version: {ver}, BB Gaussian sims chi2: {chi2_sim:.3f}, reduced chi2: {reduced_chi2_sim:.3f}, PTE: {pte_sim:.3f}" ) - chi2, reduced_chi2, pte = get_chi2_and_pte(cl_bb, cov_cl_bb, verbose=False) ax2.errorbar( jiterred_ell, cl_bb * ell, - yerr=np.sqrt(np.diag(cov_cl_bb)) * ell, + yerr=np.sqrt(np.diag(cov_cl_bb_sims)) * ell, label=labels[i], color=colors[i], fmt=markers[i], capsize=2, + markersize=5, ) ax2.text( - 100, - 1.25e-7 - i * 0.25e-7, - rf"$\chi^2_\nu={reduced_chi2:.3f}$ | PTE={pte:.3f}", + 10, + 1.5e-7 - i * 0.28e-7, + rf"$\chi^2_{{\nu, \mathrm{{iNKA}}}}={reduced_chi2_inka:.3f}$, ${{\rm PTE}}_{{\rm iNKA}}={pte_inka:.3f}$ | $\chi^2_{{\nu, \mathrm{{sim.}}}}={reduced_chi2_sim:.3f}$, ${{\rm PTE}}_{{\rm sim.}}={pte_sim:.3f}$", color=colors[i], + fontsize=15, ) ax1.set_xscale("squareroot") @@ -261,6 +324,7 @@ ax1.set_xticks(minor_ticks, minor=True) ax1.tick_params(axis="x", which="minor", length=2, width=0.8) ax1.axhline(0, color="gray", linestyle="--", linewidth=1) +ax1.set_ylim(-1.35e-7, 7.32e-8) ax2.set_xscale("squareroot") ax2.set_xticks(np.array([100, 400, 900, 1600])) @@ -269,7 +333,14 @@ ax2.set_xticks(minor_ticks, minor=True) ax2.tick_params(axis="x", which="minor", length=2, width=0.8) ax2.axhline(0, color="gray", linestyle="--", linewidth=1) -ax2.legend(fontsize=14) +ax2.legend( + fontsize=18, + bbox_to_anchor=(0.5, -0.5), + loc="lower center", + ncol=2, + frameon=False, +) +ax2.set_ylim(-3.56e-8, 1.83e-7) ax2.set_xlabel(r"Multipole $\ell$", fontsize=16) ax1.set_ylabel(r"$\ell \, C_\ell^{EB} \times 10^{-7}$", fontsize=16) @@ -285,3 +356,20 @@ # %% +mask = (ell > 300) & (ell < 1600) + +chi2_inka, reduced_chi2_inka, pte_inka = get_chi2_and_pte( + cl_bb[:-4], cov_cl_bb[:-4, :-4], verbose=True +) +chi2_sim, reduced_chi2_sim, pte_sim = get_chi2_and_pte( + cl_bb[:-4], cov_cl_bb_sims[:-4, :-4], verbose=True +) +print( + f"Version: {ver}, BB iNKA chi2: {chi2_inka:.3f}, reduced chi2: {reduced_chi2_inka:.3f}, PTE: {pte_inka:.3f}" +) +print( + f"Version: {ver}, BB Gaussian sims chi2: {chi2_sim:.3f}, reduced chi2: {reduced_chi2_sim:.3f}, PTE: {pte_sim:.3f}" +) +# %% +ell +# %% diff --git a/papers/harmonic/2025_10_28_plot_whisker.ipynb b/papers/harmonic/2025_10_28_plot_whisker.ipynb index 4a6515b63..67f18be0d 100644 --- a/papers/harmonic/2025_10_28_plot_whisker.ipynb +++ b/papers/harmonic/2025_10_28_plot_whisker.ipynb @@ -17,36 +17,43 @@ "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", "\n", - "sys.path.append(\"/home/guerrini/sp_validation/cosmo_inference/scripts/\")\n", - "\n", - "import warnings\n", + "sys.path.append(\n", + " \"/home/guerrini/sp_validation/cosmo_inference/scripts/\"\n", + ")\n", "\n", - "import chain_postprocessing as cp\n", - "import matplotlib.pyplot as plt\n", + "from getdist import plots, loadMCSamples\n", "import numpy as np\n", + "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", - "from getdist import plots\n", + "import warnings\n", + "import chain_postprocessing as cp\n", "\n", - "plt.style.use(\"/home/guerrini/matplotlib_config/paper.mplstyle\")\n", + "plt.style.use(\n", + " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", + ")\n", "\n", - "plt.rc(\"text\", usetex=True)\n", + "plt.rc('text', usetex=True)\n", "\n", "sns.set_palette(\"husl\")\n", "\n", "g = plots.get_subplot_plotter(width_inch=30)\n", - "g.settings.axes_fontsize = 60\n", - "g.settings.axes_labelsize = 60\n", + "g.settings.axes_fontsize=60\n", + "g.settings.axes_labelsize=60\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 60\n", "\n", "%matplotlib inline\n", "\n", - "# SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", "root_dir = \"/n09data/guerrini/output_chains/\"\n", "root_external = \"/n09data/guerrini/output_chains/ext_data/\"\n", "blind = \"B\"\n", "\n", - "colour_blind = {\"A\": \"royalblue\", \"B\": \"crimson\", \"C\": \"forestgreen\"}\n", + "colour_blind = {\n", + " \"A\": \"royalblue\",\n", + " \"B\": \"crimson\",\n", + " \"C\": \"forestgreen\"\n", + "}\n", "\n", "roots = [\n", " f\"SP_v1.4.6.3_leak_corr_{blind}\",\n", @@ -70,11 +77,17 @@ " f\"SP_v1.4.6.3_leak_corr_HMCode_nobar_{blind}\",\n", " f\"SP_v1.4.6.3_leak_corr_OneCov_{blind}\",\n", " f\"SP_v1.4.6.3_{blind}\",\n", + " f\"SP_v1.4.6.3_leak_corr_unif_dz_{blind}\",\n", + " f\"SP_v1.4.6.3_leak_corr_no_dz_{blind}\",\n", + " f\"SP_v1.4.6.3_leak_corr_unif_A_IA_{blind}\",\n", + " f\"SP_v1.4.6.3_leak_corr_no_IA_{blind}\",\n", + " f\"SP_v1.4.6.3_leak_corr_no_m_bias_{blind}\",\n", + " f\"SP_v1.3.6.3_leak_corr\"\n", "]\n", "\n", "legend_labels = [\n", " r\"UNIONS $C_\\ell$ (This work)\",\n", - " r\"UNIONS $\\xi_\\pm(\\vartheta)$, (Goh et al., 2026)\",\n", + " r\"UNIONS $\\xi_\\pm(\\vartheta)$ (Goh et al., 2026)\",\n", " r\"\\textit{Planck} 2018\",\n", " r\"DES Y3 $\\xi_\\pm(\\vartheta)$\",\n", " r\"DES Y3 $C_\\ell$\",\n", @@ -94,12 +107,18 @@ " r\"\\texttt{HMCode} no baryons\",\n", " r\"\\texttt{OneCovariance} only\",\n", " r\"No leakage correction\",\n", + " r\"Uniform $\\Delta z$\",\n", + " r\"No $\\Delta z$\",\n", + " r\"Uniform $A_{\\rm IA}$\",\n", + " r\"No $A_{\\rm IA}$\",\n", + " r\"No $m$-bias correction\",\n", + " r\"MCCD PSF model\"\n", "]\n", "\n", "colours = [\n", " colour_blind[\"B\"],\n", " \"darkorange\",\n", - " \"violet\",\n", + " \"olive\",\n", " \"black\",\n", " \"black\",\n", " \"black\",\n", @@ -118,6 +137,12 @@ " colour_blind[\"B\"],\n", " colour_blind[\"B\"],\n", " colour_blind[\"B\"],\n", + " colour_blind[\"B\"],\n", + " colour_blind[\"B\"],\n", + " colour_blind[\"B\"],\n", + " colour_blind[\"B\"],\n", + " colour_blind[\"B\"],\n", + " colour_blind[\"B\"]\n", "]\n", "\n", "categories = [\n", @@ -142,15 +167,46 @@ " \"harmonic\",\n", " \"harmonic\",\n", " \"harmonic\",\n", + " \"harmonic\",\n", + " \"harmonic\",\n", + " \"harmonic\",\n", + " \"harmonic\",\n", + " \"harmonic\",\n", + " \"harmonic\"\n", "]\n", "\n", + "assert len(roots) == len(legend_labels) == len(colours) == len(categories), f\"Lists must have the same length\\n roots: {len(roots)}, legend_labels: {len(legend_labels)}, colours: {len(colours)}, categories: {len(categories)}\"\n", + "\n", + "\n", "# Add the blinds to each list\n", "for bl in [\"A\", \"B\", \"C\"]:\n", " if bl != blind:\n", - " roots.append(f\"SP_v1.4.6.3_leak_corr_{bl}\")\n", - " legend_labels.append(rf\"UNIONS $C_\\ell$, Blind {bl}\")\n", - " colours.append(colour_blind[bl])\n", - " categories.append(\"harmonic\")\n", + " roots.append(\n", + " f\"SP_v1.4.6.3_leak_corr_{bl}\"\n", + " )\n", + " legend_labels.append(\n", + " fr\"UNIONS $C_\\ell$, Blind {bl}\"\n", + " )\n", + " colours.append(\n", + " colour_blind[bl]\n", + " )\n", + " categories.append(\n", + " \"harmonic\"\n", + " )\n", + "\n", + "for bl in [\"A\", \"B\", \"C\"]:\n", + " roots.append(\n", + " f\"SP_v1.4.6.3_leak_corr_unblind_{bl}\"\n", + " )\n", + " legend_labels.append(\n", + " fr\"$\\langle \\Delta z \\rangle = -0.030$, Blind {bl}\"\n", + " )\n", + " colours.append(\n", + " colour_blind[bl]\n", + " )\n", + " categories.append(\n", + " \"harmonic\"\n", + " )\n", "\n", "print(roots)" ] @@ -166,34 +222,53 @@ "for i, root in enumerate(roots):\n", " category = categories[i]\n", " if category != \"external\":\n", - " if category == \"configuration\":\n", - " path_samples = os.path.join(root_dir, f\"{root}/samples_{root}.txt\")\n", - " path_getdist = os.path.join(root_dir, f\"{root}/getdist_{root}\")\n", - " elif category == \"harmonic\":\n", + " if category == 'configuration':\n", + " path_samples = os.path.join(\n", + " root_dir,\n", + " f\"{root}/samples_{root}.txt\"\n", + " )\n", + " path_getdist = os.path.join(\n", + " root_dir,\n", + " f\"{root}/getdist_{root}\"\n", + " )\n", + " elif category == 'harmonic':\n", " path_samples = os.path.join(\n", - " root_dir, f\"{root}/{root}/samples_{root}_cell.txt\"\n", + " root_dir,\n", + " f\"{root}/{root}/samples_{root}_cell.txt\"\n", + " )\n", + " path_getdist = os.path.join(\n", + " root_dir,\n", + " f\"{root}/{root}/getdist_{root}\"\n", " )\n", - " path_getdist = os.path.join(root_dir, f\"{root}/{root}/getdist_{root}\")\n", " elif category == \"external_compute_sample\":\n", - " path_samples = os.path.join(root_dir, f\"ext_data/{root}/samples_{root}.txt\")\n", - " path_getdist = os.path.join(root_dir, f\"ext_data/{root}/getdist_{root}\")\n", + " path_samples = os.path.join(\n", + " root_dir,\n", + " f\"ext_data/{root}/samples_{root}.txt\"\n", + " )\n", + " path_getdist = os.path.join(\n", + " root_dir,\n", + " f\"ext_data/{root}/getdist_{root}\"\n", + " )\n", " else:\n", " raise ValueError(f\"The category, {category}, of {root} is not correct\")\n", - "\n", + " \n", " if category == \"external_compute_sample\" and \"Legacy\" in root:\n", - " chain_type = \"nautilus\"\n", + " chain_type=\"nautilus\"\n", " else:\n", - " chain_type = \"polychord\"\n", - " cp.load_samples_and_write_paramnames(\n", - " path_samples, path_getdist + \".paramnames\", chain_type=chain_type\n", - " )\n", - " cp.write_samples_getdist_format(\n", - " path_samples, path_getdist + \".txt\", chain_type=chain_type\n", + " chain_type=\"polychord\"\n", + " cp.load_samples_and_write_paramnames(path_samples, path_getdist+\".paramnames\", chain_type=chain_type)\n", + " cp.write_samples_getdist_format(path_samples, path_getdist+\".txt\", chain_type=chain_type)\n", + " chains.append(\n", + " cp.load_chain(path_getdist, smoothing_scale=0.5)\n", " )\n", - " chains.append(cp.load_chain(path_getdist, smoothing_scale=0.5))\n", " else:\n", - " path_getdist = os.path.join(root_dir, f\"ext_data/{root}/getdist_{root}\")\n", - " chains.append(cp.load_chain(path_getdist))" + " path_getdist = os.path.join(\n", + " root_dir,\n", + " f\"ext_data/{root}/getdist_{root}\"\n", + " )\n", + " chains.append(\n", + " cp.load_chain(path_getdist)\n", + " )" ] }, { @@ -203,32 +278,8 @@ "metadata": {}, "outputs": [], "source": [ - "name_list = [\n", - " \"OMEGA_M\",\n", - " \"ombh2\",\n", - " \"h0\",\n", - " \"n_s\",\n", - " \"SIGMA_8\",\n", - " \"S_8\",\n", - " \"s_8_input\",\n", - " \"logt_agn\",\n", - " \"a\",\n", - " \"m1\",\n", - " \"bias_1\",\n", - "]\n", - "label_list = [\n", - " r\"\\Omega_{\\rm m}\",\n", - " r\"\\omega_b h^2\",\n", - " r\"h_0\",\n", - " r\"n_s\",\n", - " r\"\\sigma_8\",\n", - " r\"S_8\",\n", - " r\"S_8\",\n", - " r\"\\log T_{\\rm AGN}\",\n", - " r\"A_{\\rm IA}\",\n", - " r\"m_1\",\n", - " r\"\\Delta z_1\",\n", - "]\n", + "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','S_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", + "label_list = [r'\\Omega_{\\rm m}', r'\\omega_b h^2', r'h_0', r'n_s', r'\\sigma_8', r'S_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", "\n", "for i, chain in enumerate(chains):\n", " print(legend_labels[i])\n", @@ -236,7 +287,7 @@ " for name, label in zip(name_list, label_list):\n", " try:\n", " param_names.parWithName(name).label = label\n", - " except Exception:\n", + " except:\n", " warnings.warn(f\"Parameter {name} not found in chain {roots[i]}.\")" ] }, @@ -249,42 +300,42 @@ "source": [ "# Micro management of external chains\n", "# Account for the missing parameter conventions\n", - "# OMEGA_M not in DES_Y3_cell\n", - "idx = roots.index(\"DES_Y3_cell\")\n", - "cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")\n", + "#OMEGA_M not in DES_Y3_cell\n", + "idx = roots.index('DES_Y3_cell')\n", + "cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')\n", "cp.derive_parameter_S8(chains[idx])\n", "\n", - "# OMEGA_M not in KiDS-1000\n", + "#OMEGA_M not in KiDS-1000\n", "try:\n", - " idx = roots.index(\"KiDS-1000\")\n", - " cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")\n", - "except Exception:\n", + " idx = roots.index('KiDS-1000')\n", + " cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')\n", + "except:\n", " print(\"KiDS-1000 chain not found, skipping parameter name adjustment.\")\n", "\n", - "# S8 to derive in KiDS-Legacy chains\n", + "#S8 to derive in KiDS-Legacy chains\n", "try:\n", - " idx = roots.index(\"KiDS-Legacy_xipm\")\n", + " idx = roots.index('KiDS-Legacy_xipm')\n", " cp.derive_parameter_S8(chains[idx])\n", - "except Exception:\n", + "except:\n", " print(\"KiDS-Legacy_xipm chain not found, skipping S8 derivation.\")\n", "try:\n", - " idx = roots.index(\"KiDS-Legacy_bandpowers\")\n", + " idx = roots.index('KiDS-Legacy_bandpowers')\n", " cp.derive_parameter_S8(chains[idx])\n", - "except Exception:\n", + "except:\n", " print(\"KiDS-Legacy_bandpowers chain not found, skipping S8 derivation.\")\n", "try:\n", - " idx = roots.index(\"KiDS-Legacy_cosebis\")\n", + " idx = roots.index('KiDS-Legacy_cosebis')\n", " cp.derive_parameter_S8(chains[idx])\n", - "except Exception:\n", + "except:\n", " print(\"KiDS-Legacy_cosebis chain not found, skipping S8 derivation.\")\n", "\n", - "# OMEGA_M not in DES+KiDS\n", - "idx = roots.index(\"DES+KiDS\")\n", - "cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")\n", + "#OMEGA_M not in DES+KiDS\n", + "idx = roots.index('DES+KiDS')\n", + "cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')\n", "\n", - "# OMEGA_M not in HSC_Y3_cell\n", - "idx = roots.index(\"HSC_Y3_cell\")\n", - "cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")" + "#OMEGA_M not in HSC_Y3_cell\n", + "idx = roots.index('HSC_Y3_cell')\n", + "cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')" ] }, { @@ -294,60 +345,39 @@ "metadata": {}, "outputs": [], "source": [ - "best_fit_method = \"2Dkde\"\n", - "\n", - "param_values = np.array(\n", - " [\n", - " \"# Expt\",\n", - " \"Colour\",\n", - " \"S8_Mean\",\n", - " \"S8_low\",\n", - " \"S8_high\",\n", - " \"sigma_8_Mean\",\n", - " \"sigma_8_low\",\n", - " \"sigma_8_high\",\n", - " \"Omega_m_Mean\",\n", - " \"Omega_m_low\",\n", - " \"Omega_m_high\",\n", - " ]\n", - ")\n", - "escaped = np.char.replace(legend_labels, \"\\\\\", \"\\\\\\\\\")\n", + "best_fit_method = \"2Dkde\" \n", + "\n", + "param_values = np.array([\"# Expt\", \"Colour\", \"S8_Mean\", \"S8_low\", \"S8_high\", \"sigma_8_Mean\", \"sigma_8_low\", \"sigma_8_high\", \"Omega_m_Mean\", \"Omega_m_low\", \"Omega_m_high\"])\n", + "escaped = np.char.replace(legend_labels, '\\\\', '\\\\\\\\')\n", "for i, chain in enumerate(chains):\n", " print(chain.root)\n", " margestats = chain.getMargeStats()\n", " likestats = chain.getLikeStats()\n", "\n", - " s8_stats = margestats.parWithName(\"S_8\")\n", - " sigma8_stats = margestats.parWithName(\"SIGMA_8\")\n", - " omegam_stats = margestats.parWithName(\"OMEGA_M\")\n", + " s8_stats = margestats.parWithName('S_8')\n", + " sigma8_stats = margestats.parWithName('SIGMA_8')\n", + " omegam_stats = margestats.parWithName('OMEGA_M')\n", "\n", " best_fit = cp.extract_best_fit_params(chain, best_fit_method=best_fit_method)\n", "\n", - " param_values = np.vstack(\n", - " (\n", - " param_values,\n", - " [\n", - " escaped[i],\n", - " colours[i],\n", - " best_fit[\"S_8\"],\n", - " best_fit[\"S_8\"] - s8_stats.limits[0].lower,\n", - " s8_stats.limits[0].upper - best_fit[\"S_8\"],\n", - " best_fit[\"SIGMA_8\"],\n", - " best_fit[\"SIGMA_8\"] - sigma8_stats.limits[0].lower,\n", - " sigma8_stats.limits[0].upper - best_fit[\"SIGMA_8\"],\n", - " best_fit[\"OMEGA_M\"],\n", - " best_fit[\"OMEGA_M\"] - omegam_stats.limits[0].lower,\n", - " omegam_stats.limits[0].upper - best_fit[\"OMEGA_M\"],\n", - " ],\n", - " )\n", - " )\n", + " param_values = np.vstack((\n", + " param_values,\n", + " [\n", + " escaped[i],\n", + " colours[i],\n", + " best_fit['S_8'],\n", + " best_fit['S_8']-s8_stats.limits[0].lower,\n", + " s8_stats.limits[0].upper-best_fit['S_8'],\n", + " best_fit['SIGMA_8'],\n", + " best_fit['SIGMA_8'] - sigma8_stats.limits[0].lower,\n", + " sigma8_stats.limits[0].upper - best_fit['SIGMA_8'],\n", + " best_fit['OMEGA_M'],\n", + " best_fit['OMEGA_M'] - omegam_stats.limits[0].lower,\n", + " omegam_stats.limits[0].upper - best_fit['OMEGA_M'],\n", + " ]\n", + " ))\n", "print(param_values)\n", - "np.savetxt(\n", - " f\"{root_dir}/param_values.txt\",\n", - " param_values,\n", - " fmt=[\"%s\" for i in range(11)],\n", - " delimiter=\";\",\n", - ")" + "np.savetxt(f\"{root_dir}/param_values.txt\", param_values, fmt=['%s' for i in range(11)], delimiter=';')" ] }, { @@ -358,50 +388,20 @@ "outputs": [], "source": [ "# Load the value of the parameters\n", - "cosmo = np.loadtxt(\n", - " f\"{root_dir}/param_values.txt\",\n", - " dtype={\n", - " \"names\": (\n", - " \"Expt\",\n", - " \"colour\",\n", - " \"s8_mean\",\n", - " \"s8_low\",\n", - " \"s8_high\",\n", - " \"sigma8_mean\",\n", - " \"sigma8_low\",\n", - " \"sigma8_high\",\n", - " \"omegam_mean\",\n", - " \"omegam_low\",\n", - " \"omegam_high\",\n", - " ),\n", - " \"formats\": (\n", - " \"U250\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " \"U20\",\n", - " ),\n", - " },\n", - " skiprows=1,\n", - " delimiter=\";\",\n", - ")\n", - "expt = np.char.replace(cosmo[\"Expt\"], \"\\\\\\\\\", \"\\\\\")\n", - "colours = cosmo[\"colour\"]\n", - "s8_mean = cosmo[\"s8_mean\"].astype(np.float64)\n", - "s8_low = cosmo[\"s8_low\"].astype(np.float64)\n", - "s8_high = cosmo[\"s8_high\"].astype(np.float64)\n", - "sigma8_mean = cosmo[\"sigma8_mean\"].astype(np.float64)\n", - "sigma8_low = cosmo[\"sigma8_low\"].astype(np.float64)\n", - "sigma8_high = cosmo[\"sigma8_high\"].astype(np.float64)\n", - "omegam_mean = cosmo[\"omegam_mean\"].astype(np.float64)\n", - "omegam_low = cosmo[\"omegam_low\"].astype(np.float64)\n", - "omegam_high = cosmo[\"omegam_high\"].astype(np.float64)" + "cosmo = np.loadtxt(f\"{root_dir}/param_values.txt\",\n", + " dtype={'names': ('Expt', 'colour', 's8_mean', 's8_low', 's8_high', 'sigma8_mean', 'sigma8_low', 'sigma8_high', 'omegam_mean', 'omegam_low', 'omegam_high'),\n", + " 'formats': ('U250', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20')}, skiprows=1, delimiter=';')\n", + "expt = np.char.replace(cosmo['Expt'], '\\\\\\\\', '\\\\')\n", + "colours = cosmo['colour']\n", + "s8_mean = cosmo['s8_mean'].astype(np.float64)\n", + "s8_low = cosmo['s8_low'].astype(np.float64)\n", + "s8_high = cosmo['s8_high'].astype(np.float64)\n", + "sigma8_mean = cosmo['sigma8_mean'].astype(np.float64)\n", + "sigma8_low = cosmo['sigma8_low'].astype(np.float64)\n", + "sigma8_high = cosmo['sigma8_high'].astype(np.float64)\n", + "omegam_mean = cosmo['omegam_mean'].astype(np.float64)\n", + "omegam_low = cosmo['omegam_low'].astype(np.float64)\n", + "omegam_high = cosmo['omegam_high'].astype(np.float64)" ] }, { @@ -417,6 +417,19 @@ " delta_mean = np.abs(mean1 - mean2)\n", " sigma_tension = delta_mean / np.sqrt(sigma1**2 + sigma2**2)\n", " sign = 1 if mean1 > mean2 else -1\n", + " return sigma_tension * sign\n", + "\n", + "def get_sigma_tension_asymmetric(mean1, low1, high1, mean2, low2, high2):\n", + " delta_mean = mean1 - mean2\n", + " if delta_mean > 0:\n", + " # mean1 is higher: use its lower error, and mean2's upper error\n", + " sigma1 = low1\n", + " sigma2 = high2\n", + " else:\n", + " sigma1 = high1\n", + " sigma2 = low2\n", + " sigma_tension = np.abs(delta_mean) / np.sqrt(sigma1**2 + sigma2**2)\n", + " sign = 1 if delta_mean > 0 else -1\n", " return sigma_tension * sign" ] }, @@ -429,7 +442,9 @@ "source": [ "from matplotlib.gridspec import GridSpec\n", "\n", - "fig = plt.figure(figsize=(10, 6))\n", + "tension_function = get_sigma_tension_asymmetric\n", + "\n", + "fig = plt.figure(figsize=(10, 8))\n", "gs = GridSpec(1, 3, width_ratios=[1, 0.5, 0.5])\n", "ax1 = fig.add_subplot(gs[0])\n", "ax2 = fig.add_subplot(gs[1], sharey=ax1)\n", @@ -444,14 +459,27 @@ "]\n", "reference = r\"UNIONS $C_\\ell$ (This work)\"\n", "separation_after = [\n", - " r\"UNIONS $\\xi_\\pm(\\vartheta)$, (Goh et al., 2026)\",\n", + " r\"UNIONS $\\xi_\\pm(\\vartheta)$ (Goh et al., 2026)\",\n", " r\"HSC Y3 $C_\\ell$\",\n", " r\"Large scales only\",\n", " r\"\\texttt{HMCode} no baryons\",\n", " r\"\\texttt{OneCovariance} only\",\n", " r\"No leakage correction\",\n", + " r\"No $m$-bias correction\",\n", + " r\"MCCD PSF model\",\n", + " r\"UNIONS $C_\\ell$, Blind C\"\n", + "]\n", + "list_section_index = [\n", + " r\"(ii)\",\n", + " r\"(iii)\",\n", + " r\"(iv)\",\n", + " r\"(v)\",\n", + " r\"(vi)\",\n", + " r\"(vii)\",\n", + " r\"(viii)\",\n", + " r\"(ix)\",\n", + " r\"(x)\"\n", "]\n", - "list_section_index = [r\"(ii)\", r\"(iii)\", r\"(iv)\", r\"(v)\", r\"(vi)\", r\"(vii)\"]\n", "\n", "preliminary_watermark = False\n", "blind_axes = False\n", @@ -463,132 +491,88 @@ "for ax, param in zip(axs, params):\n", " means, lows, highs, label = param\n", " for i, mean, low, high, color in zip(y, means, lows, highs, colours):\n", - " ax.errorbar(\n", - " mean,\n", - " 0.05 + i * row_spacing,\n", - " xerr=np.array([low, high])[:, None],\n", - " fmt=\"o\",\n", - " color=color,\n", - " ecolor=color,\n", - " elinewidth=2,\n", - " capsize=3,\n", - " )\n", + " ax.errorbar(mean, 0.05+i*row_spacing, xerr=np.array([low, high])[:, None], fmt='o', color=color, ecolor=color, elinewidth=2, capsize=3)\n", " ax.set_xlabel(label, fontsize=14)\n", - "\n", + " \n", " ax.grid(False)\n", - " ax.tick_params(axis=\"y\", left=False, labelleft=False)\n", + " ax.tick_params(axis='y', left=False, labelleft=False)\n", " if label == r\"$S_8$\":\n", - " ax.axvspan(\n", - " s8_mean[index_ref] - s8_low[index_ref],\n", - " s8_mean[index_ref] + s8_high[index_ref],\n", - " color=colours[index_ref],\n", - " alpha=0.2,\n", - " )\n", - " ax.set_xlim(0.20, 1.1)\n", + " ax.axvspan(s8_mean[index_ref] - s8_low[index_ref], s8_mean[index_ref] + s8_high[index_ref], color=colours[index_ref], alpha=0.2)\n", + " ax.set_xlim(-0.1, 1.41)\n", " if blind_axes:\n", " ref_tick = np.mean(s8_mean[:4])\n", - " ax.set_xticks([ref_tick + i * 0.1 for i in range(-5, 5)], labels=[])\n", + " ax.set_xticks(\n", + " [ref_tick + i*0.1 for i in range(-5, 5)], labels=[]\n", + " )\n", " else:\n", - " ax.set_xticks([0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0])\n", + " ax.set_xticks(\n", + " [0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4],\n", + " labels=[\n", + " r\"$0.6$\", r\"\", r\"$0.8$\", r\"\", r\"$1.0$\", r\"\", r\"$1.2$\", r\"\", r\"$1.4$\"\n", + " ]\n", + " )\n", " elif label == r\"$\\sigma_8$\":\n", - " ax.axvspan(\n", - " sigma8_mean[index_ref] - sigma8_low[index_ref],\n", - " sigma8_mean[index_ref] + sigma8_high[index_ref],\n", - " color=colours[index_ref],\n", - " alpha=0.2,\n", - " )\n", - " ax.set_xlim(0.5, 1.2)\n", + " ax.axvspan(sigma8_mean[index_ref] - sigma8_low[index_ref], sigma8_mean[index_ref] + sigma8_high[index_ref], color=colours[index_ref], alpha=0.2)\n", + " ax.set_xlim(0.5, 1.25)\n", " if blind_axes:\n", " ref_tick = np.mean(sigma8_mean[:4])\n", - " ax.set_xticks([ref_tick + i * 0.2 for i in range(-2, 2)], labels=[])\n", + " ax.set_xticks(\n", + " [ref_tick + i*0.2 for i in range(-2, 2)], labels=[]\n", + " )\n", " elif label == r\"$\\Omega_{\\rm m}$\":\n", - " ax.axvspan(\n", - " omegam_mean[index_ref] - omegam_low[index_ref],\n", - " omegam_mean[index_ref] + omegam_high[index_ref],\n", - " color=colours[index_ref],\n", - " alpha=0.2,\n", - " )\n", + " ax.axvspan(omegam_mean[index_ref] - omegam_low[index_ref], omegam_mean[index_ref] + omegam_high[index_ref], color=colours[index_ref], alpha=0.2)\n", " ax.set_xlim(0.1, 0.5)\n", " if blind_axes:\n", " ref_tick = np.mean(omegam_mean[:4])\n", - " ax.set_xticks([ref_tick + i * 0.1 for i in range(-2, 3)], labels=[])\n", + " ax.set_xticks(\n", + " [ref_tick + i*0.1 for i in range(-2, 3)], labels=[]\n", + " )\n", "\n", "\n", - "axs[0].set_yticks(0.05 + y * row_spacing)\n", + "axs[0].set_yticks(0.05+y*row_spacing)\n", "axs[0].set_yticklabels([])\n", "for label, color in zip(expt, colours):\n", - " axs[0].text(\n", - " 0.21,\n", - " 0.05 + row_spacing * np.where(expt == label)[0][0],\n", - " label,\n", - " fontsize=12,\n", - " ha=\"left\",\n", - " va=\"center\",\n", - " color=color,\n", - " )\n", + " axs[0].text(-0.09, 0.05 + row_spacing * np.where(expt == label)[0][0], label, fontsize=11, ha='left', va='center', color=color)\n", " if label != reference:\n", " index = np.where(expt == label)[0][0]\n", - " s8_tension = get_sigma_tension(\n", - " s8_mean[index],\n", - " s8_low[index],\n", - " s8_high[index],\n", - " s8_mean[index_ref],\n", - " s8_low[index_ref],\n", - " s8_high[index_ref],\n", + " s8_tension = tension_function(\n", + " s8_mean[index], s8_low[index], s8_high[index],\n", + " s8_mean[index_ref], s8_low[index_ref], s8_high[index_ref]\n", " )\n", " sign_str = \"+\" if s8_tension > 0 else \"-\"\n", - " axs[0].text(\n", - " 1.0,\n", - " 0.05 + row_spacing * index,\n", - " rf\"${sign_str}{np.abs(s8_tension):.2f}\" + r\"\\, \\sigma$\",\n", - " fontsize=10,\n", - " ha=\"left\",\n", - " va=\"center\",\n", - " color=color,\n", - " )\n", + " text_to_plot = rf\"${sign_str}{np.abs(s8_tension):.2f}\" + r\"\\, \\sigma$\"\n", + " if label == r\"UNIONS $\\xi_\\pm(\\vartheta)$ (Goh et al., 2026)\":\n", + " text_to_plot += r\" (ind.)\"\n", + " axs[0].text(1.38, 0.05 + row_spacing * (index - 0.45), text_to_plot, fontsize=10, ha='right', va='center', color=color)\n", + " additional_text = r\"$-2.18 \\, \\sigma$ (corr.)\"\n", + " axs[0].text(1.4, 0.05 + row_spacing * (index + 0.15), additional_text, fontsize=10, ha='right', va='center', color=color)\n", + " else:\n", + " axs[0].text(1.38, 0.05 + row_spacing * index, text_to_plot, fontsize=10, ha='right', va='center', color=color)\n", "# Add separation lines\n", "for i, sep in enumerate(separation_after):\n", " index_sep = np.where(expt == sep)[0][0]\n", " for ax in axs:\n", - " ax.axhline(\n", - " row_spacing * (index_sep + 0.95),\n", - " color=\"black\",\n", - " linestyle=\"dotted\",\n", - " linewidth=1,\n", - " )\n", - " axs[0].text(\n", - " 0.20,\n", - " 0.05 + row_spacing * (index_sep + 1),\n", - " list_section_index[i],\n", - " fontsize=14,\n", - " fontweight=\"bold\",\n", - " va=\"center\",\n", - " ha=\"right\",\n", - " )\n", + " ax.axhline(row_spacing * (index_sep + 0.95), color='black', linestyle='dotted', linewidth=1)\n", + " axs[0].text(-0.1, 0.05 + row_spacing * (index_sep + 1), \n", + " list_section_index[i], fontsize=14, fontweight='bold', va='center', ha='right')\n", "\n", "\n", "# --- Add section labels (i), (ii)) ---\n", - "axs[0].text(0.20, 0.05, r\"(i)\", fontsize=14, fontweight=\"bold\", va=\"center\", ha=\"right\")\n", + "axs[0].text(-0.1, 0.05, \n", + " r\"(i)\", fontsize=14, fontweight='bold', va='center', ha='right')\n", "\n", "if preliminary_watermark:\n", - " plt.figtext(\n", - " 0.5,\n", - " 0.5,\n", - " \"PRELIMINARY\",\n", - " fontsize=50,\n", - " color=\"gray\",\n", - " ha=\"center\",\n", - " va=\"center\",\n", - " alpha=0.3,\n", - " rotation=330,\n", - " )\n", + " plt.figtext(0.5, 0.5, 'PRELIMINARY',\n", + " fontsize=50, color='gray',\n", + " ha='center', va='center',\n", + " alpha=0.3, rotation=330)\n", "\n", "plt.gca().invert_yaxis()\n", "\n", "plt.tight_layout()\n", "\n", "plt.savefig(\"./plots/whisker_plot.png\", dpi=300)\n", - "# Save pdf\n", + "#Save pdf\n", "plt.savefig(\"./plots/whisker_plot.pdf\")\n", "plt.show()" ] From cacdf0c1c8ae310d90cacfbacd43893a0176f419 Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Mon, 2 Nov 2026 19:38:46 +0000 Subject: [PATCH 2/2] ruff autofix (format + safe lint fixes) Pushed by the lint gate. --- .../2025_09_26_plot_contours_IA_prior.ipynb | 151 +++--- .../2025_09_26_plot_contours_bias_prior.ipynb | 153 +++--- .../2025_09_26_plot_contours_cl_vs_xi.ipynb | 191 ++++--- .../2025_09_26_plot_contours_unblind.ipynb | 153 +++--- papers/harmonic/2025_10_28_plot_whisker.ipynb | 472 +++++++++++------- 5 files changed, 658 insertions(+), 462 deletions(-) diff --git a/papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb b/papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb index 1a6a2fdfa..443d11f8f 100644 --- a/papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb +++ b/papers/harmonic/2025_09_26_plot_contours_IA_prior.ipynb @@ -13,26 +13,24 @@ "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", "\n", - "from getdist import plots, loadMCSamples\n", - "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "import seaborn as sns\n", + "from getdist import plots\n", "\n", - "plt.style.use(\n", - " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", - ")\n", + "plt.style.use(\"/home/guerrini/matplotlib_config/paper.mplstyle\")\n", "\n", "sns.set_palette(\"husl\")\n", "\n", "g = plots.get_subplot_plotter(width_inch=30)\n", - "g.settings.axes_fontsize=60\n", - "g.settings.axes_labelsize=60\n", + "g.settings.axes_fontsize = 60\n", + "g.settings.axes_labelsize = 60\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 60\n", "\n", "%matplotlib inline\n", "\n", - "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "# SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", "root_dir = \"/n09data/guerrini/output_chains/\"\n", "\n", "roots = [\n", @@ -57,16 +55,21 @@ "# MAKE PARAMNAMES FILE\n", "\n", "for root in roots:\n", - " with open(root_dir + '{}/{}/samples_{}_cell.txt'.format('/'+root ,root, root), \"r\") as file:\n", - " params = file.readline()[1:].split('\\t')[:-4]\n", + " with open(\n", + " root_dir + \"{}/{}/samples_{}_cell.txt\".format(\"/\" + root, root, root), \"r\"\n", + " ) as file:\n", + " params = file.readline()[1:].split(\"\\t\")[:-4]\n", " file.close()\n", - " \n", - " with open(root_dir + '{}/{}/getdist_{}_cell.paramnames'.format('/'+root, root,root), \"w\") as file:\n", + "\n", + " with open(\n", + " root_dir + \"{}/{}/getdist_{}_cell.paramnames\".format(\"/\" + root, root, root),\n", + " \"w\",\n", + " ) as file:\n", " for i in range(len(params)):\n", - " if len(params[i].split('--')) > 1:\n", - " file.write(params[i].split('--')[1] + '\\n')\n", + " if len(params[i].split(\"--\")) > 1:\n", + " file.write(params[i].split(\"--\")[1] + \"\\n\")\n", " else:\n", - " file.write(params[i].split('--')[0] + '\\n')\n", + " file.write(params[i].split(\"--\")[0] + \"\\n\")\n", " file.close()" ] }, @@ -77,25 +80,28 @@ "metadata": {}, "outputs": [], "source": [ - "#READ CHAIN\n", + "# READ CHAIN\n", "\n", - "chains=[]\n", + "chains = []\n", "\n", "for root in roots:\n", - "\n", - " samples = np.loadtxt(root_dir + '{}/{}/samples_{}_cell.txt'.format(root,root,root))\n", + " samples = np.loadtxt(\n", + " root_dir + \"{}/{}/samples_{}_cell.txt\".format(root, root, root)\n", + " )\n", " print(len(samples))\n", - " if 'nautilus' in root:\n", - " samples = np.column_stack((np.exp(samples[:,-3]),samples[:,-1]-samples[:,-2],samples[:,0:-3]))\n", + " if \"nautilus\" in root:\n", + " samples = np.column_stack(\n", + " (np.exp(samples[:, -3]), samples[:, -1] - samples[:, -2], samples[:, 0:-3])\n", + " )\n", " else:\n", - " samples = np.column_stack((samples[:,-1],samples[:,-3],samples[:,0:-4]))\n", - " np.savetxt(root_dir + '{}/{}/getdist_{}_cell.txt'.format(root,root,root), samples)\n", - " \n", - " chain = g.samples_for_root(root_dir + '{}/{}/getdist_{}_cell'.format(root,root,root),\n", - " cache=False,\n", - " settings={'ignore_rows':0,\n", - " 'smooth_scale_2D':0.5,\n", - " 'smooth_scale_1D':0.5})\n", + " samples = np.column_stack((samples[:, -1], samples[:, -3], samples[:, 0:-4]))\n", + " np.savetxt(root_dir + \"{}/{}/getdist_{}_cell.txt\".format(root, root, root), samples)\n", + "\n", + " chain = g.samples_for_root(\n", + " root_dir + \"{}/{}/getdist_{}_cell\".format(root, root, root),\n", + " cache=False,\n", + " settings={\"ignore_rows\": 0, \"smooth_scale_2D\": 0.5, \"smooth_scale_1D\": 0.5},\n", + " )\n", "\n", " chains.append(chain)" ] @@ -107,8 +113,30 @@ "metadata": {}, "outputs": [], "source": [ - "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", - "label_list = label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "name_list = [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + "]\n", + "label_list = label_list = [\n", + " r\"\\Omega_{\\rm m}\",\n", + " r\"\\omega_{\\rm b} h^2\",\n", + " r\"h_0\",\n", + " r\"n_{\\rm s}\",\n", + " r\"\\sigma_8\",\n", + " r\"S_8\",\n", + " r\"\\log T_{\\rm AGN}\",\n", + " r\"A_{\\rm IA}\",\n", + " r\"m_1\",\n", + " r\"\\Delta z_1\",\n", + "]\n", "\n", "for chain in chains:\n", " param_names = chain.getParamNames()\n", @@ -129,32 +157,35 @@ " r\"No $A_{\\rm IA}$\",\n", "]\n", "\n", - "colours = [\n", - " \"crimson\",\n", - " \"royalblue\",\n", - " \"forestgreen\"\n", - "]\n", + "colours = [\"crimson\", \"royalblue\", \"forestgreen\"]\n", "\n", - "linestyle = [\n", - " 'solid',\n", - " 'solid',\n", - " 'solid'\n", - "]\n", + "linestyle = [\"solid\", \"solid\", \"solid\"]\n", "\n", "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", "\n", - "#Plot all parameters\n", + "# Plot all parameters\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + " ],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc='upper right',\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_ia_prior/contours_all_cell.png')\n", + "g.export(\"plots/cosmological_constraints_ia_prior/contours_all_cell.png\")\n", "plt.show()" ] }, @@ -165,19 +196,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Plot only cosmological parameters\n", + "# Plot only cosmological parameters\n", "\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','s_8_input', 'SIGMA_8', 'a'],\n", + " [\"OMEGA_M\", \"s_8_input\", \"SIGMA_8\", \"a\"],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc='upper right',\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_ia_prior/contours_cosmo_cell.png')\n", + "g.export(\"plots/cosmological_constraints_ia_prior/contours_cosmo_cell.png\")\n", "plt.show()" ] }, @@ -188,19 +219,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Plot S8 Omega_m only\n", - "#Plot S8 Omega_m only\n", + "# Plot S8 Omega_m only\n", + "# Plot S8 Omega_m only\n", "g = plots.get_subplot_plotter(width_inch=20)\n", - "g.settings.axes_fontsize=35\n", - "g.settings.axes_labelsize=50\n", + "g.settings.axes_fontsize = 35\n", + "g.settings.axes_labelsize = 50\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 40\n", "\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','s_8_input'],\n", + " [\"OMEGA_M\", \"s_8_input\"],\n", " legend_labels=legend_labels,\n", - " legend_loc='upper right',\n", + " legend_loc=\"upper right\",\n", " line_args=line_args,\n", " contour_colors=colours,\n", " filled=True,\n", @@ -211,9 +242,9 @@ " ha='center', va='center',\n", " alpha=0.3, rotation=330) \"\"\"\n", "\n", - "g.export('plots/cosmological_constraints_ia_prior/contours_s8_omegam_cell.png')\n", - "#Save PDF\n", - "g.export('plots/cosmological_constraints_ia_prior/contours_s8_omegam_cell.pdf')\n", + "g.export(\"plots/cosmological_constraints_ia_prior/contours_s8_omegam_cell.png\")\n", + "# Save PDF\n", + "g.export(\"plots/cosmological_constraints_ia_prior/contours_s8_omegam_cell.pdf\")\n", "plt.show()" ] }, diff --git a/papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb b/papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb index 8df25bed3..94df62f8e 100644 --- a/papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb +++ b/papers/harmonic/2025_09_26_plot_contours_bias_prior.ipynb @@ -13,26 +13,24 @@ "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", "\n", - "from getdist import plots, loadMCSamples\n", - "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "import seaborn as sns\n", + "from getdist import plots\n", "\n", - "plt.style.use(\n", - " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", - ")\n", + "plt.style.use(\"/home/guerrini/matplotlib_config/paper.mplstyle\")\n", "\n", "sns.set_palette(\"husl\")\n", "\n", "g = plots.get_subplot_plotter(width_inch=30)\n", - "g.settings.axes_fontsize=60\n", - "g.settings.axes_labelsize=60\n", + "g.settings.axes_fontsize = 60\n", + "g.settings.axes_labelsize = 60\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 60\n", "\n", "%matplotlib inline\n", "\n", - "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "# SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", "root_dir = \"/n09data/guerrini/output_chains/\"\n", "\n", "roots = [\n", @@ -57,16 +55,21 @@ "# MAKE PARAMNAMES FILE\n", "\n", "for root in roots:\n", - " with open(root_dir + '{}/{}/samples_{}_cell.txt'.format('/'+root ,root, root), \"r\") as file:\n", - " params = file.readline()[1:].split('\\t')[:-4]\n", + " with open(\n", + " root_dir + \"{}/{}/samples_{}_cell.txt\".format(\"/\" + root, root, root), \"r\"\n", + " ) as file:\n", + " params = file.readline()[1:].split(\"\\t\")[:-4]\n", " file.close()\n", - " \n", - " with open(root_dir + '{}/{}/getdist_{}_cell.paramnames'.format('/'+root, root,root), \"w\") as file:\n", + "\n", + " with open(\n", + " root_dir + \"{}/{}/getdist_{}_cell.paramnames\".format(\"/\" + root, root, root),\n", + " \"w\",\n", + " ) as file:\n", " for i in range(len(params)):\n", - " if len(params[i].split('--')) > 1:\n", - " file.write(params[i].split('--')[1] + '\\n')\n", + " if len(params[i].split(\"--\")) > 1:\n", + " file.write(params[i].split(\"--\")[1] + \"\\n\")\n", " else:\n", - " file.write(params[i].split('--')[0] + '\\n')\n", + " file.write(params[i].split(\"--\")[0] + \"\\n\")\n", " file.close()" ] }, @@ -77,25 +80,28 @@ "metadata": {}, "outputs": [], "source": [ - "#READ CHAIN\n", + "# READ CHAIN\n", "\n", - "chains=[]\n", + "chains = []\n", "\n", "for root in roots:\n", - "\n", - " samples = np.loadtxt(root_dir + '{}/{}/samples_{}_cell.txt'.format(root,root,root))\n", + " samples = np.loadtxt(\n", + " root_dir + \"{}/{}/samples_{}_cell.txt\".format(root, root, root)\n", + " )\n", " print(len(samples))\n", - " if 'nautilus' in root:\n", - " samples = np.column_stack((np.exp(samples[:,-3]),samples[:,-1]-samples[:,-2],samples[:,0:-3]))\n", + " if \"nautilus\" in root:\n", + " samples = np.column_stack(\n", + " (np.exp(samples[:, -3]), samples[:, -1] - samples[:, -2], samples[:, 0:-3])\n", + " )\n", " else:\n", - " samples = np.column_stack((samples[:,-1],samples[:,-3],samples[:,0:-4]))\n", - " np.savetxt(root_dir + '{}/{}/getdist_{}_cell.txt'.format(root,root,root), samples)\n", - " \n", - " chain = g.samples_for_root(root_dir + '{}/{}/getdist_{}_cell'.format(root,root,root),\n", - " cache=False,\n", - " settings={'ignore_rows':0,\n", - " 'smooth_scale_2D':0.5,\n", - " 'smooth_scale_1D':0.5})\n", + " samples = np.column_stack((samples[:, -1], samples[:, -3], samples[:, 0:-4]))\n", + " np.savetxt(root_dir + \"{}/{}/getdist_{}_cell.txt\".format(root, root, root), samples)\n", + "\n", + " chain = g.samples_for_root(\n", + " root_dir + \"{}/{}/getdist_{}_cell\".format(root, root, root),\n", + " cache=False,\n", + " settings={\"ignore_rows\": 0, \"smooth_scale_2D\": 0.5, \"smooth_scale_1D\": 0.5},\n", + " )\n", "\n", " chains.append(chain)" ] @@ -107,8 +113,30 @@ "metadata": {}, "outputs": [], "source": [ - "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", - "label_list = label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "name_list = [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + "]\n", + "label_list = label_list = [\n", + " r\"\\Omega_{\\rm m}\",\n", + " r\"\\omega_{\\rm b} h^2\",\n", + " r\"h_0\",\n", + " r\"n_{\\rm s}\",\n", + " r\"\\sigma_8\",\n", + " r\"S_8\",\n", + " r\"\\log T_{\\rm AGN}\",\n", + " r\"A_{\\rm IA}\",\n", + " r\"m_1\",\n", + " r\"\\Delta z_1\",\n", + "]\n", "\n", "for chain in chains:\n", " param_names = chain.getParamNames()\n", @@ -126,35 +154,38 @@ "legend_labels = [\n", " r\"Uniform $\\Delta z$ prior\",\n", " r\"No $\\Delta z$ prior\",\n", - " r\"No $m$ bias prior\"\n", + " r\"No $m$ bias prior\",\n", "]\n", "\n", - "colours = [\n", - " \"royalblue\",\n", - " \"crimson\",\n", - " \"forestgreen\"\n", - "]\n", + "colours = [\"royalblue\", \"crimson\", \"forestgreen\"]\n", "\n", - "linestyle = [\n", - " 'solid',\n", - " 'solid',\n", - " 'solid'\n", - "]\n", + "linestyle = [\"solid\", \"solid\", \"solid\"]\n", "\n", "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", "\n", - "#Plot all parameters\n", + "# Plot all parameters\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + " ],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc='upper right',\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_bias_prior/contours_all_cell.png')\n", + "g.export(\"plots/cosmological_constraints_bias_prior/contours_all_cell.png\")\n", "plt.show()" ] }, @@ -165,19 +196,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Plot only cosmological parameters\n", + "# Plot only cosmological parameters\n", "\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','s_8_input', 'SIGMA_8', 'a'],\n", + " [\"OMEGA_M\", \"s_8_input\", \"SIGMA_8\", \"a\"],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc='upper right',\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_bias_prior/contours_cosmo_cell.png')\n", + "g.export(\"plots/cosmological_constraints_bias_prior/contours_cosmo_cell.png\")\n", "plt.show()" ] }, @@ -188,19 +219,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Plot S8 Omega_m only\n", - "#Plot S8 Omega_m only\n", + "# Plot S8 Omega_m only\n", + "# Plot S8 Omega_m only\n", "g = plots.get_subplot_plotter(width_inch=20)\n", - "g.settings.axes_fontsize=35\n", - "g.settings.axes_labelsize=50\n", + "g.settings.axes_fontsize = 35\n", + "g.settings.axes_labelsize = 50\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 40\n", "\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','s_8_input'],\n", + " [\"OMEGA_M\", \"s_8_input\"],\n", " legend_labels=legend_labels,\n", - " legend_loc='upper right',\n", + " legend_loc=\"upper right\",\n", " line_args=line_args,\n", " contour_colors=colours,\n", " filled=True,\n", @@ -211,9 +242,9 @@ " ha='center', va='center',\n", " alpha=0.3, rotation=330) \"\"\"\n", "\n", - "g.export('plots/cosmological_constraints_bias_prior/contours_s8_omegam_cell.png')\n", - "#Save PDF\n", - "g.export('plots/cosmological_constraints_bias_prior/contours_s8_omegam_cell.pdf')\n", + "g.export(\"plots/cosmological_constraints_bias_prior/contours_s8_omegam_cell.png\")\n", + "# Save PDF\n", + "g.export(\"plots/cosmological_constraints_bias_prior/contours_s8_omegam_cell.pdf\")\n", "plt.show()" ] }, diff --git a/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb b/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb index 35c0180e7..a2f254ceb 100644 --- a/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb +++ b/papers/harmonic/2025_09_26_plot_contours_cl_vs_xi.ipynb @@ -16,42 +16,32 @@ "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", "\n", "# Append any useful folder in the path\n", - "sys.path.append(\n", - " \"/home/guerrini/sp_validation/cosmo_inference/scripts/\"\n", - ")\n", + "sys.path.append(\"/home/guerrini/sp_validation/cosmo_inference/scripts/\")\n", "\n", - "from getdist import plots, loadMCSamples\n", - "import numpy as np\n", + "import chain_postprocessing as cp\n", "import matplotlib.pyplot as plt\n", "import seaborn as sns\n", + "from getdist import plots\n", "\n", - "import chain_postprocessing as cp\n", - "\n", - "plt.style.use(\n", - " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", - ")\n", + "plt.style.use(\"/home/guerrini/matplotlib_config/paper.mplstyle\")\n", "\n", "sns.set_palette(\"husl\")\n", "\n", "g = plots.get_subplot_plotter(width_inch=20)\n", - "g.settings.axes_fontsize=60\n", - "g.settings.axes_labelsize=60\n", + "g.settings.axes_fontsize = 60\n", + "g.settings.axes_labelsize = 60\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 50\n", "\n", "%matplotlib inline\n", "\n", - "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "# SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", "root_dir = \"/n09data/guerrini/output_chains/\"\n", "root_external = \"/n09data/guerrini/output_chains/ext_data/\"\n", "\n", "blind = \"B\"\n", "\n", - "colour_blind = {\n", - " \"A\": \"royalblue\",\n", - " \"B\": \"crimson\",\n", - " \"C\": \"forestgreen\"\n", - "}\n", + "colour_blind = {\"A\": \"royalblue\", \"B\": \"crimson\", \"C\": \"forestgreen\"}\n", "\n", "roots = [\n", " f\"SP_v1.4.6.3_leak_corr_{blind}_prior_volume\",\n", @@ -60,19 +50,9 @@ " \"Planck18\",\n", "]\n", "\n", - "categories = [\n", - " \"harmonic\",\n", - " \"harmonic\",\n", - " \"configuration\",\n", - " \"external\"\n", - "]\n", + "categories = [\"harmonic\", \"harmonic\", \"configuration\", \"external\"]\n", "\n", - "need_compute_getdist = [\n", - " True,\n", - " True,\n", - " True,\n", - " False\n", - "]\n", + "need_compute_getdist = [True, True, True, False]\n", "\n", "chain_types = [\n", " \"apriori\",\n", @@ -98,48 +78,26 @@ "for i, root in enumerate(roots):\n", " analysis_type_i = categories[i]\n", " if analysis_type_i == \"harmonic\":\n", - " path_samples = os.path.join(\n", - " root_dir,\n", - " root,\n", - " root,\n", - " f\"samples_{root}_cell.txt\"\n", - " )\n", - " path_gd = os.path.join(\n", - " root_dir,\n", - " root,\n", - " root,\n", - " f\"getdist_{root}_cell\"\n", - " )\n", + " path_samples = os.path.join(root_dir, root, root, f\"samples_{root}_cell.txt\")\n", + " path_gd = os.path.join(root_dir, root, root, f\"getdist_{root}_cell\")\n", " elif analysis_type_i == \"configuration\":\n", - " path_samples = os.path.join(\n", - " root_dir,\n", - " root,\n", - " f\"samples_{root}.txt\"\n", - " )\n", - " path_gd = os.path.join(\n", - " root_dir,\n", - " root,\n", - " f\"getdist_{root}\"\n", - " )\n", + " path_samples = os.path.join(root_dir, root, f\"samples_{root}.txt\")\n", + " path_gd = os.path.join(root_dir, root, f\"getdist_{root}\")\n", " elif analysis_type_i == \"external\":\n", - " path_samples = os.path.join(\n", - " root_external,\n", - " root,\n", - " f\"samples_{root}.txt\"\n", - " )\n", - " path_gd = os.path.join(\n", - " root_external,\n", - " root,\n", - " f\"getdist_{root}\"\n", - " )\n", + " path_samples = os.path.join(root_external, root, f\"samples_{root}.txt\")\n", + " path_gd = os.path.join(root_external, root, f\"getdist_{root}\")\n", " else:\n", " raise ValueError(f\"Invalid analysis type: {analysis_type_i}\")\n", - " \n", + "\n", " need_getdist = need_compute_getdist[i]\n", " if need_getdist:\n", " # Read files, format and load chains\n", - " cp.load_samples_and_write_paramnames(path_samples, path_gd+\".paramnames\", chain_type=chain_types[i])\n", - " cp.write_samples_getdist_format(path_samples, path_gd+\".txt\", chain_type=chain_types[i])\n", + " cp.load_samples_and_write_paramnames(\n", + " path_samples, path_gd + \".paramnames\", chain_type=chain_types[i]\n", + " )\n", + " cp.write_samples_getdist_format(\n", + " path_samples, path_gd + \".txt\", chain_type=chain_types[i]\n", + " )\n", " chain = cp.load_chain(path_gd, smoothing_scale=0.5)\n", " chains.append(chain)" ] @@ -151,8 +109,32 @@ "metadata": {}, "outputs": [], "source": [ - "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8', 'S_8', 's_8_input', 'logt_agn','a','m1','bias_1']\n", - "label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "name_list = [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"S_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + "]\n", + "label_list = [\n", + " r\"\\Omega_{\\rm m}\",\n", + " r\"\\omega_{\\rm b} h^2\",\n", + " r\"h_0\",\n", + " r\"n_{\\rm s}\",\n", + " r\"\\sigma_8\",\n", + " r\"S_8\",\n", + " r\"S_8\",\n", + " r\"\\log T_{\\rm AGN}\",\n", + " r\"A_{\\rm IA}\",\n", + " r\"m_1\",\n", + " r\"\\Delta z_1\",\n", + "]\n", "\n", "for chain in chains:\n", " param_names = chain.getParamNames()\n", @@ -184,37 +166,42 @@ " \"olive\",\n", "]\n", "\n", - "linestyle = [\n", - " 'solid',\n", - " 'solid',\n", - " 'solid',\n", - " 'solid'\n", - "]\n", + "linestyle = [\"solid\", \"solid\", \"solid\", \"solid\"]\n", "\n", "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", "\n", "idx = [0, 1, 2]\n", "g = plots.get_subplot_plotter(width_inch=30)\n", - "g.settings.axes_fontsize=20\n", - "g.settings.axes_labelsize=30\n", + "g.settings.axes_fontsize = 20\n", + "g.settings.axes_labelsize = 30\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 50\n", "\n", - "#Plot all parameters\n", + "# Plot all parameters\n", "g.triangle_plot(\n", - " [chains[i] for i in idx], # Plot everything but Planck\n", - " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " [chains[i] for i in idx], # Plot everything but Planck\n", + " [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + " ],\n", " legend_labels=[legend_labels[i] for i in idx],\n", " line_args=[line_args[i] for i in idx],\n", " contour_colors=[colours[i] for i in idx],\n", - " legend_loc='upper right',\n", - "\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_cl_vs_xi/contours_all_cell.png')\n", - "#Save pdf\n", - "g.export('plots/cosmological_constraints_cl_vs_xi/contours_all_cell.pdf')\n", + "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_all_cell.png\")\n", + "# Save pdf\n", + "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_all_cell.pdf\")\n", "plt.show()" ] }, @@ -225,24 +212,24 @@ "metadata": {}, "outputs": [], "source": [ - "#Plot only cosmological parameters\n", + "# Plot only cosmological parameters\n", "g = plots.get_subplot_plotter(width_inch=20)\n", - "g.settings.axes_fontsize=35\n", - "g.settings.axes_labelsize=50\n", + "g.settings.axes_fontsize = 35\n", + "g.settings.axes_labelsize = 50\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 40\n", "\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','S_8', 'SIGMA_8', 'a'],\n", + " [\"OMEGA_M\", \"S_8\", \"SIGMA_8\", \"a\"],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc='upper right',\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_cl_vs_xi/contours_cosmo_cell.png')\n", + "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_cosmo_cell.png\")\n", "plt.show()" ] }, @@ -253,17 +240,17 @@ "metadata": {}, "outputs": [], "source": [ - "g.settings.axes_fontsize=35\n", - "g.settings.axes_labelsize=50\n", + "g.settings.axes_fontsize = 35\n", + "g.settings.axes_labelsize = 50\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 40\n", "\n", - "#Plot S8 Omega_m only\n", + "# Plot S8 Omega_m only\n", "g.triangle_plot(\n", - " chains[:], # Plot everything but the prior volume\n", - " ['OMEGA_M','S_8'],\n", + " chains[:], # Plot everything but the prior volume\n", + " [\"OMEGA_M\", \"S_8\"],\n", " legend_labels=legend_labels,\n", - " legend_loc='upper right',\n", + " legend_loc=\"upper right\",\n", " line_args=line_args,\n", " contour_colors=colours,\n", " filled=True,\n", @@ -272,18 +259,18 @@ "# Set S_8 range\n", "s8_min, s8_max = 0.58, 1.15 # adjust to what you want\n", "\n", - "g.subplots[1, 0].set_ylim(s8_min, s8_max) # 2D contour panel (S_8 on y-axis)\n", + "g.subplots[1, 0].set_ylim(s8_min, s8_max) # 2D contour panel (S_8 on y-axis)\n", "g.subplots[1, 0].set_xlim(0.05, 0.65)\n", - "g.subplots[1, 1].set_xlim(s8_min, s8_max) # S_8 marginal panel (S_8 on x-axis)\n", + "g.subplots[1, 1].set_xlim(s8_min, s8_max) # S_8 marginal panel (S_8 on x-axis)\n", "\n", "\"\"\" plt.figtext(0.5, 0.5, 'PRELIMINARY',\n", " fontsize=150, color='gray',\n", " ha='center', va='center',\n", " alpha=0.3, rotation=330) \"\"\"\n", "\n", - "g.export('plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.png')\n", - "#Save pdf\n", - "g.export('plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.pdf')\n", + "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.png\")\n", + "# Save pdf\n", + "g.export(\"plots/cosmological_constraints_cl_vs_xi/contours_s8_omegam_cell.pdf\")\n", "plt.show()" ] }, diff --git a/papers/harmonic/2025_09_26_plot_contours_unblind.ipynb b/papers/harmonic/2025_09_26_plot_contours_unblind.ipynb index 7a85ea81c..416ad8301 100644 --- a/papers/harmonic/2025_09_26_plot_contours_unblind.ipynb +++ b/papers/harmonic/2025_09_26_plot_contours_unblind.ipynb @@ -13,26 +13,24 @@ "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", "\n", - "from getdist import plots, loadMCSamples\n", - "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "import seaborn as sns\n", + "from getdist import plots\n", "\n", - "plt.style.use(\n", - " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", - ")\n", + "plt.style.use(\"/home/guerrini/matplotlib_config/paper.mplstyle\")\n", "\n", "sns.set_palette(\"husl\")\n", "\n", "g = plots.get_subplot_plotter(width_inch=30)\n", - "g.settings.axes_fontsize=60\n", - "g.settings.axes_labelsize=60\n", + "g.settings.axes_fontsize = 60\n", + "g.settings.axes_labelsize = 60\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 60\n", "\n", "%matplotlib inline\n", "\n", - "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "# SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", "root_dir = \"/n09data/guerrini/output_chains/\"\n", "\n", "roots = [\n", @@ -57,16 +55,21 @@ "# MAKE PARAMNAMES FILE\n", "\n", "for root in roots:\n", - " with open(root_dir + '{}/{}/samples_{}_cell.txt'.format('/'+root ,root, root), \"r\") as file:\n", - " params = file.readline()[1:].split('\\t')[:-4]\n", + " with open(\n", + " root_dir + \"{}/{}/samples_{}_cell.txt\".format(\"/\" + root, root, root), \"r\"\n", + " ) as file:\n", + " params = file.readline()[1:].split(\"\\t\")[:-4]\n", " file.close()\n", - " \n", - " with open(root_dir + '{}/{}/getdist_{}_cell.paramnames'.format('/'+root, root,root), \"w\") as file:\n", + "\n", + " with open(\n", + " root_dir + \"{}/{}/getdist_{}_cell.paramnames\".format(\"/\" + root, root, root),\n", + " \"w\",\n", + " ) as file:\n", " for i in range(len(params)):\n", - " if len(params[i].split('--')) > 1:\n", - " file.write(params[i].split('--')[1] + '\\n')\n", + " if len(params[i].split(\"--\")) > 1:\n", + " file.write(params[i].split(\"--\")[1] + \"\\n\")\n", " else:\n", - " file.write(params[i].split('--')[0] + '\\n')\n", + " file.write(params[i].split(\"--\")[0] + \"\\n\")\n", " file.close()" ] }, @@ -77,25 +80,28 @@ "metadata": {}, "outputs": [], "source": [ - "#READ CHAIN\n", + "# READ CHAIN\n", "\n", - "chains=[]\n", + "chains = []\n", "\n", "for root in roots:\n", - "\n", - " samples = np.loadtxt(root_dir + '{}/{}/samples_{}_cell.txt'.format(root,root,root))\n", + " samples = np.loadtxt(\n", + " root_dir + \"{}/{}/samples_{}_cell.txt\".format(root, root, root)\n", + " )\n", " print(len(samples))\n", - " if 'nautilus' in root:\n", - " samples = np.column_stack((np.exp(samples[:,-3]),samples[:,-1]-samples[:,-2],samples[:,0:-3]))\n", + " if \"nautilus\" in root:\n", + " samples = np.column_stack(\n", + " (np.exp(samples[:, -3]), samples[:, -1] - samples[:, -2], samples[:, 0:-3])\n", + " )\n", " else:\n", - " samples = np.column_stack((samples[:,-1],samples[:,-3],samples[:,0:-4]))\n", - " np.savetxt(root_dir + '{}/{}/getdist_{}_cell.txt'.format(root,root,root), samples)\n", - " \n", - " chain = g.samples_for_root(root_dir + '{}/{}/getdist_{}_cell'.format(root,root,root),\n", - " cache=False,\n", - " settings={'ignore_rows':0,\n", - " 'smooth_scale_2D':0.5,\n", - " 'smooth_scale_1D':0.5})\n", + " samples = np.column_stack((samples[:, -1], samples[:, -3], samples[:, 0:-4]))\n", + " np.savetxt(root_dir + \"{}/{}/getdist_{}_cell.txt\".format(root, root, root), samples)\n", + "\n", + " chain = g.samples_for_root(\n", + " root_dir + \"{}/{}/getdist_{}_cell\".format(root, root, root),\n", + " cache=False,\n", + " settings={\"ignore_rows\": 0, \"smooth_scale_2D\": 0.5, \"smooth_scale_1D\": 0.5},\n", + " )\n", "\n", " chains.append(chain)" ] @@ -107,8 +113,30 @@ "metadata": {}, "outputs": [], "source": [ - "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", - "label_list = label_list = [r'\\Omega_{\\rm m}', r'\\omega_{\\rm b} h^2', r'h_0', r'n_{\\rm s}', r'\\sigma_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "name_list = [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + "]\n", + "label_list = label_list = [\n", + " r\"\\Omega_{\\rm m}\",\n", + " r\"\\omega_{\\rm b} h^2\",\n", + " r\"h_0\",\n", + " r\"n_{\\rm s}\",\n", + " r\"\\sigma_8\",\n", + " r\"S_8\",\n", + " r\"\\log T_{\\rm AGN}\",\n", + " r\"A_{\\rm IA}\",\n", + " r\"m_1\",\n", + " r\"\\Delta z_1\",\n", + "]\n", "\n", "for chain in chains:\n", " param_names = chain.getParamNames()\n", @@ -124,35 +152,38 @@ "outputs": [], "source": [ "legend_labels = [\n", - " fr\"$\\langle \\Delta z \\rangle = -0.030$, Blind {bl}\" for bl in ['A', 'B', 'C']\n", + " rf\"$\\langle \\Delta z \\rangle = -0.030$, Blind {bl}\" for bl in [\"A\", \"B\", \"C\"]\n", "]\n", "\n", - "colours = [\n", - " \"royalblue\",\n", - " \"crimson\",\n", - " \"forestgreen\"\n", - "]\n", + "colours = [\"royalblue\", \"crimson\", \"forestgreen\"]\n", "\n", - "linestyle = [\n", - " 'solid',\n", - " 'solid',\n", - " 'solid'\n", - "]\n", + "linestyle = [\"solid\", \"solid\", \"solid\"]\n", "\n", "line_args = [dict(color=col, ls=ls) for col, ls in zip(colours, linestyle)]\n", "\n", - "#Plot all parameters\n", + "# Plot all parameters\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','ombh2', 'h0','n_s','SIGMA_8','s_8_input', 'logt_agn','a','m1','bias_1'],\n", + " [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + " ],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc='upper right',\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_unblind/contours_all_cell.png')\n", + "g.export(\"plots/cosmological_constraints_unblind/contours_all_cell.png\")\n", "plt.show()" ] }, @@ -163,19 +194,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Plot only cosmological parameters\n", + "# Plot only cosmological parameters\n", "\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','s_8_input', 'SIGMA_8', 'a'],\n", + " [\"OMEGA_M\", \"s_8_input\", \"SIGMA_8\", \"a\"],\n", " legend_labels=legend_labels,\n", " line_args=line_args,\n", " contour_colors=colours,\n", - " legend_loc='upper right',\n", - " filled=True\n", + " legend_loc=\"upper right\",\n", + " filled=True,\n", ")\n", "\n", - "g.export('plots/cosmological_constraints_unblind/contours_cosmo_cell.png')\n", + "g.export(\"plots/cosmological_constraints_unblind/contours_cosmo_cell.png\")\n", "plt.show()" ] }, @@ -186,19 +217,19 @@ "metadata": {}, "outputs": [], "source": [ - "#Plot S8 Omega_m only\n", - "#Plot S8 Omega_m only\n", + "# Plot S8 Omega_m only\n", + "# Plot S8 Omega_m only\n", "g = plots.get_subplot_plotter(width_inch=20)\n", - "g.settings.axes_fontsize=35\n", - "g.settings.axes_labelsize=50\n", + "g.settings.axes_fontsize = 35\n", + "g.settings.axes_labelsize = 50\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 40\n", "\n", "g.triangle_plot(\n", " chains,\n", - " ['OMEGA_M','s_8_input'],\n", + " [\"OMEGA_M\", \"s_8_input\"],\n", " legend_labels=legend_labels,\n", - " legend_loc='upper right',\n", + " legend_loc=\"upper right\",\n", " line_args=line_args,\n", " contour_colors=colours,\n", " filled=True,\n", @@ -209,9 +240,9 @@ " ha='center', va='center',\n", " alpha=0.3, rotation=330) \"\"\"\n", "\n", - "g.export('plots/cosmological_constraints_unblind/contours_s8_omegam_cell.png')\n", - "#Save PDF\n", - "g.export('plots/cosmological_constraints_unblind/contours_s8_omegam_cell.pdf')\n", + "g.export(\"plots/cosmological_constraints_unblind/contours_s8_omegam_cell.png\")\n", + "# Save PDF\n", + "g.export(\"plots/cosmological_constraints_unblind/contours_s8_omegam_cell.pdf\")\n", "plt.show()" ] }, diff --git a/papers/harmonic/2025_10_28_plot_whisker.ipynb b/papers/harmonic/2025_10_28_plot_whisker.ipynb index 67f18be0d..baca9f876 100644 --- a/papers/harmonic/2025_10_28_plot_whisker.ipynb +++ b/papers/harmonic/2025_10_28_plot_whisker.ipynb @@ -17,43 +17,36 @@ "os.environ[\"LD_LIBRARY_PATH\"] = \"\"\n", "os.environ[\"CONDA_PREFIX\"] = \"/home/guerrini/.conda/envs/sp_validation_3.11\"\n", "\n", - "sys.path.append(\n", - " \"/home/guerrini/sp_validation/cosmo_inference/scripts/\"\n", - ")\n", + "sys.path.append(\"/home/guerrini/sp_validation/cosmo_inference/scripts/\")\n", "\n", - "from getdist import plots, loadMCSamples\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import seaborn as sns\n", "import warnings\n", + "\n", "import chain_postprocessing as cp\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import seaborn as sns\n", + "from getdist import plots\n", "\n", - "plt.style.use(\n", - " \"/home/guerrini/matplotlib_config/paper.mplstyle\"\n", - ")\n", + "plt.style.use(\"/home/guerrini/matplotlib_config/paper.mplstyle\")\n", "\n", - "plt.rc('text', usetex=True)\n", + "plt.rc(\"text\", usetex=True)\n", "\n", "sns.set_palette(\"husl\")\n", "\n", "g = plots.get_subplot_plotter(width_inch=30)\n", - "g.settings.axes_fontsize=60\n", - "g.settings.axes_labelsize=60\n", + "g.settings.axes_fontsize = 60\n", + "g.settings.axes_labelsize = 60\n", "g.settings.alpha_filled_add = 0.7\n", "g.settings.legend_fontsize = 60\n", "\n", "%matplotlib inline\n", "\n", - "#SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", + "# SPECIFY DATA DIRECTORY AND DESIRED CHAINS TO ANALYSE\n", "root_dir = \"/n09data/guerrini/output_chains/\"\n", "root_external = \"/n09data/guerrini/output_chains/ext_data/\"\n", "blind = \"B\"\n", "\n", - "colour_blind = {\n", - " \"A\": \"royalblue\",\n", - " \"B\": \"crimson\",\n", - " \"C\": \"forestgreen\"\n", - "}\n", + "colour_blind = {\"A\": \"royalblue\", \"B\": \"crimson\", \"C\": \"forestgreen\"}\n", "\n", "roots = [\n", " f\"SP_v1.4.6.3_leak_corr_{blind}\",\n", @@ -82,7 +75,7 @@ " f\"SP_v1.4.6.3_leak_corr_unif_A_IA_{blind}\",\n", " f\"SP_v1.4.6.3_leak_corr_no_IA_{blind}\",\n", " f\"SP_v1.4.6.3_leak_corr_no_m_bias_{blind}\",\n", - " f\"SP_v1.3.6.3_leak_corr\"\n", + " \"SP_v1.3.6.3_leak_corr\",\n", "]\n", "\n", "legend_labels = [\n", @@ -112,7 +105,7 @@ " r\"Uniform $A_{\\rm IA}$\",\n", " r\"No $A_{\\rm IA}$\",\n", " r\"No $m$-bias correction\",\n", - " r\"MCCD PSF model\"\n", + " r\"MCCD PSF model\",\n", "]\n", "\n", "colours = [\n", @@ -142,7 +135,7 @@ " colour_blind[\"B\"],\n", " colour_blind[\"B\"],\n", " colour_blind[\"B\"],\n", - " colour_blind[\"B\"]\n", + " colour_blind[\"B\"],\n", "]\n", "\n", "categories = [\n", @@ -172,41 +165,27 @@ " \"harmonic\",\n", " \"harmonic\",\n", " \"harmonic\",\n", - " \"harmonic\"\n", + " \"harmonic\",\n", "]\n", "\n", - "assert len(roots) == len(legend_labels) == len(colours) == len(categories), f\"Lists must have the same length\\n roots: {len(roots)}, legend_labels: {len(legend_labels)}, colours: {len(colours)}, categories: {len(categories)}\"\n", + "assert len(roots) == len(legend_labels) == len(colours) == len(categories), (\n", + " f\"Lists must have the same length\\n roots: {len(roots)}, legend_labels: {len(legend_labels)}, colours: {len(colours)}, categories: {len(categories)}\"\n", + ")\n", "\n", "\n", "# Add the blinds to each list\n", "for bl in [\"A\", \"B\", \"C\"]:\n", " if bl != blind:\n", - " roots.append(\n", - " f\"SP_v1.4.6.3_leak_corr_{bl}\"\n", - " )\n", - " legend_labels.append(\n", - " fr\"UNIONS $C_\\ell$, Blind {bl}\"\n", - " )\n", - " colours.append(\n", - " colour_blind[bl]\n", - " )\n", - " categories.append(\n", - " \"harmonic\"\n", - " )\n", + " roots.append(f\"SP_v1.4.6.3_leak_corr_{bl}\")\n", + " legend_labels.append(rf\"UNIONS $C_\\ell$, Blind {bl}\")\n", + " colours.append(colour_blind[bl])\n", + " categories.append(\"harmonic\")\n", "\n", "for bl in [\"A\", \"B\", \"C\"]:\n", - " roots.append(\n", - " f\"SP_v1.4.6.3_leak_corr_unblind_{bl}\"\n", - " )\n", - " legend_labels.append(\n", - " fr\"$\\langle \\Delta z \\rangle = -0.030$, Blind {bl}\"\n", - " )\n", - " colours.append(\n", - " colour_blind[bl]\n", - " )\n", - " categories.append(\n", - " \"harmonic\"\n", - " )\n", + " roots.append(f\"SP_v1.4.6.3_leak_corr_unblind_{bl}\")\n", + " legend_labels.append(rf\"$\\langle \\Delta z \\rangle = -0.030$, Blind {bl}\")\n", + " colours.append(colour_blind[bl])\n", + " categories.append(\"harmonic\")\n", "\n", "print(roots)" ] @@ -222,53 +201,34 @@ "for i, root in enumerate(roots):\n", " category = categories[i]\n", " if category != \"external\":\n", - " if category == 'configuration':\n", - " path_samples = os.path.join(\n", - " root_dir,\n", - " f\"{root}/samples_{root}.txt\"\n", - " )\n", - " path_getdist = os.path.join(\n", - " root_dir,\n", - " f\"{root}/getdist_{root}\"\n", - " )\n", - " elif category == 'harmonic':\n", + " if category == \"configuration\":\n", + " path_samples = os.path.join(root_dir, f\"{root}/samples_{root}.txt\")\n", + " path_getdist = os.path.join(root_dir, f\"{root}/getdist_{root}\")\n", + " elif category == \"harmonic\":\n", " path_samples = os.path.join(\n", - " root_dir,\n", - " f\"{root}/{root}/samples_{root}_cell.txt\"\n", - " )\n", - " path_getdist = os.path.join(\n", - " root_dir,\n", - " f\"{root}/{root}/getdist_{root}\"\n", + " root_dir, f\"{root}/{root}/samples_{root}_cell.txt\"\n", " )\n", + " path_getdist = os.path.join(root_dir, f\"{root}/{root}/getdist_{root}\")\n", " elif category == \"external_compute_sample\":\n", - " path_samples = os.path.join(\n", - " root_dir,\n", - " f\"ext_data/{root}/samples_{root}.txt\"\n", - " )\n", - " path_getdist = os.path.join(\n", - " root_dir,\n", - " f\"ext_data/{root}/getdist_{root}\"\n", - " )\n", + " path_samples = os.path.join(root_dir, f\"ext_data/{root}/samples_{root}.txt\")\n", + " path_getdist = os.path.join(root_dir, f\"ext_data/{root}/getdist_{root}\")\n", " else:\n", " raise ValueError(f\"The category, {category}, of {root} is not correct\")\n", - " \n", + "\n", " if category == \"external_compute_sample\" and \"Legacy\" in root:\n", - " chain_type=\"nautilus\"\n", + " chain_type = \"nautilus\"\n", " else:\n", - " chain_type=\"polychord\"\n", - " cp.load_samples_and_write_paramnames(path_samples, path_getdist+\".paramnames\", chain_type=chain_type)\n", - " cp.write_samples_getdist_format(path_samples, path_getdist+\".txt\", chain_type=chain_type)\n", - " chains.append(\n", - " cp.load_chain(path_getdist, smoothing_scale=0.5)\n", + " chain_type = \"polychord\"\n", + " cp.load_samples_and_write_paramnames(\n", + " path_samples, path_getdist + \".paramnames\", chain_type=chain_type\n", " )\n", - " else:\n", - " path_getdist = os.path.join(\n", - " root_dir,\n", - " f\"ext_data/{root}/getdist_{root}\"\n", + " cp.write_samples_getdist_format(\n", + " path_samples, path_getdist + \".txt\", chain_type=chain_type\n", " )\n", - " chains.append(\n", - " cp.load_chain(path_getdist)\n", - " )" + " chains.append(cp.load_chain(path_getdist, smoothing_scale=0.5))\n", + " else:\n", + " path_getdist = os.path.join(root_dir, f\"ext_data/{root}/getdist_{root}\")\n", + " chains.append(cp.load_chain(path_getdist))" ] }, { @@ -278,8 +238,32 @@ "metadata": {}, "outputs": [], "source": [ - "name_list = ['OMEGA_M','ombh2','h0','n_s','SIGMA_8','S_8','s_8_input', 'logt_agn','a','m1','bias_1']\n", - "label_list = [r'\\Omega_{\\rm m}', r'\\omega_b h^2', r'h_0', r'n_s', r'\\sigma_8', r'S_8', r'S_8', r'\\log T_{\\rm AGN}', r'A_{\\rm IA}', r'm_1', r'\\Delta z_1']\n", + "name_list = [\n", + " \"OMEGA_M\",\n", + " \"ombh2\",\n", + " \"h0\",\n", + " \"n_s\",\n", + " \"SIGMA_8\",\n", + " \"S_8\",\n", + " \"s_8_input\",\n", + " \"logt_agn\",\n", + " \"a\",\n", + " \"m1\",\n", + " \"bias_1\",\n", + "]\n", + "label_list = [\n", + " r\"\\Omega_{\\rm m}\",\n", + " r\"\\omega_b h^2\",\n", + " r\"h_0\",\n", + " r\"n_s\",\n", + " r\"\\sigma_8\",\n", + " r\"S_8\",\n", + " r\"S_8\",\n", + " r\"\\log T_{\\rm AGN}\",\n", + " r\"A_{\\rm IA}\",\n", + " r\"m_1\",\n", + " r\"\\Delta z_1\",\n", + "]\n", "\n", "for i, chain in enumerate(chains):\n", " print(legend_labels[i])\n", @@ -300,42 +284,42 @@ "source": [ "# Micro management of external chains\n", "# Account for the missing parameter conventions\n", - "#OMEGA_M not in DES_Y3_cell\n", - "idx = roots.index('DES_Y3_cell')\n", - "cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')\n", + "# OMEGA_M not in DES_Y3_cell\n", + "idx = roots.index(\"DES_Y3_cell\")\n", + "cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")\n", "cp.derive_parameter_S8(chains[idx])\n", "\n", - "#OMEGA_M not in KiDS-1000\n", + "# OMEGA_M not in KiDS-1000\n", "try:\n", - " idx = roots.index('KiDS-1000')\n", - " cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')\n", + " idx = roots.index(\"KiDS-1000\")\n", + " cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")\n", "except:\n", " print(\"KiDS-1000 chain not found, skipping parameter name adjustment.\")\n", "\n", - "#S8 to derive in KiDS-Legacy chains\n", + "# S8 to derive in KiDS-Legacy chains\n", "try:\n", - " idx = roots.index('KiDS-Legacy_xipm')\n", + " idx = roots.index(\"KiDS-Legacy_xipm\")\n", " cp.derive_parameter_S8(chains[idx])\n", "except:\n", " print(\"KiDS-Legacy_xipm chain not found, skipping S8 derivation.\")\n", "try:\n", - " idx = roots.index('KiDS-Legacy_bandpowers')\n", + " idx = roots.index(\"KiDS-Legacy_bandpowers\")\n", " cp.derive_parameter_S8(chains[idx])\n", "except:\n", " print(\"KiDS-Legacy_bandpowers chain not found, skipping S8 derivation.\")\n", "try:\n", - " idx = roots.index('KiDS-Legacy_cosebis')\n", + " idx = roots.index(\"KiDS-Legacy_cosebis\")\n", " cp.derive_parameter_S8(chains[idx])\n", "except:\n", " print(\"KiDS-Legacy_cosebis chain not found, skipping S8 derivation.\")\n", "\n", - "#OMEGA_M not in DES+KiDS\n", - "idx = roots.index('DES+KiDS')\n", - "cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')\n", + "# OMEGA_M not in DES+KiDS\n", + "idx = roots.index(\"DES+KiDS\")\n", + "cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")\n", "\n", - "#OMEGA_M not in HSC_Y3_cell\n", - "idx = roots.index('HSC_Y3_cell')\n", - "cp.adjust_paramname_chain(chains[idx], 'omega_m', 'OMEGA_M', r'\\Omega_{\\rm m}')" + "# OMEGA_M not in HSC_Y3_cell\n", + "idx = roots.index(\"HSC_Y3_cell\")\n", + "cp.adjust_paramname_chain(chains[idx], \"omega_m\", \"OMEGA_M\", r\"\\Omega_{\\rm m}\")" ] }, { @@ -345,39 +329,60 @@ "metadata": {}, "outputs": [], "source": [ - "best_fit_method = \"2Dkde\" \n", - "\n", - "param_values = np.array([\"# Expt\", \"Colour\", \"S8_Mean\", \"S8_low\", \"S8_high\", \"sigma_8_Mean\", \"sigma_8_low\", \"sigma_8_high\", \"Omega_m_Mean\", \"Omega_m_low\", \"Omega_m_high\"])\n", - "escaped = np.char.replace(legend_labels, '\\\\', '\\\\\\\\')\n", + "best_fit_method = \"2Dkde\"\n", + "\n", + "param_values = np.array(\n", + " [\n", + " \"# Expt\",\n", + " \"Colour\",\n", + " \"S8_Mean\",\n", + " \"S8_low\",\n", + " \"S8_high\",\n", + " \"sigma_8_Mean\",\n", + " \"sigma_8_low\",\n", + " \"sigma_8_high\",\n", + " \"Omega_m_Mean\",\n", + " \"Omega_m_low\",\n", + " \"Omega_m_high\",\n", + " ]\n", + ")\n", + "escaped = np.char.replace(legend_labels, \"\\\\\", \"\\\\\\\\\")\n", "for i, chain in enumerate(chains):\n", " print(chain.root)\n", " margestats = chain.getMargeStats()\n", " likestats = chain.getLikeStats()\n", "\n", - " s8_stats = margestats.parWithName('S_8')\n", - " sigma8_stats = margestats.parWithName('SIGMA_8')\n", - " omegam_stats = margestats.parWithName('OMEGA_M')\n", + " s8_stats = margestats.parWithName(\"S_8\")\n", + " sigma8_stats = margestats.parWithName(\"SIGMA_8\")\n", + " omegam_stats = margestats.parWithName(\"OMEGA_M\")\n", "\n", " best_fit = cp.extract_best_fit_params(chain, best_fit_method=best_fit_method)\n", "\n", - " param_values = np.vstack((\n", - " param_values,\n", - " [\n", - " escaped[i],\n", - " colours[i],\n", - " best_fit['S_8'],\n", - " best_fit['S_8']-s8_stats.limits[0].lower,\n", - " s8_stats.limits[0].upper-best_fit['S_8'],\n", - " best_fit['SIGMA_8'],\n", - " best_fit['SIGMA_8'] - sigma8_stats.limits[0].lower,\n", - " sigma8_stats.limits[0].upper - best_fit['SIGMA_8'],\n", - " best_fit['OMEGA_M'],\n", - " best_fit['OMEGA_M'] - omegam_stats.limits[0].lower,\n", - " omegam_stats.limits[0].upper - best_fit['OMEGA_M'],\n", - " ]\n", - " ))\n", + " param_values = np.vstack(\n", + " (\n", + " param_values,\n", + " [\n", + " escaped[i],\n", + " colours[i],\n", + " best_fit[\"S_8\"],\n", + " best_fit[\"S_8\"] - s8_stats.limits[0].lower,\n", + " s8_stats.limits[0].upper - best_fit[\"S_8\"],\n", + " best_fit[\"SIGMA_8\"],\n", + " best_fit[\"SIGMA_8\"] - sigma8_stats.limits[0].lower,\n", + " sigma8_stats.limits[0].upper - best_fit[\"SIGMA_8\"],\n", + " best_fit[\"OMEGA_M\"],\n", + " best_fit[\"OMEGA_M\"] - omegam_stats.limits[0].lower,\n", + " omegam_stats.limits[0].upper - best_fit[\"OMEGA_M\"],\n", + " ],\n", + " )\n", + " )\n", "print(param_values)\n", - "np.savetxt(f\"{root_dir}/param_values.txt\", param_values, fmt=['%s' for i in range(11)], delimiter=';')" + "np.savetxt(\n", + " f\"{root_dir}/param_values.txt\",\n", + " param_values,\n", + " fmt=[\"%s\" for i in range(11)],\n", + " delimiter=\";\",\n", + ")" ] }, { @@ -388,20 +393,50 @@ "outputs": [], "source": [ "# Load the value of the parameters\n", - "cosmo = np.loadtxt(f\"{root_dir}/param_values.txt\",\n", - " dtype={'names': ('Expt', 'colour', 's8_mean', 's8_low', 's8_high', 'sigma8_mean', 'sigma8_low', 'sigma8_high', 'omegam_mean', 'omegam_low', 'omegam_high'),\n", - " 'formats': ('U250', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20', 'U20')}, skiprows=1, delimiter=';')\n", - "expt = np.char.replace(cosmo['Expt'], '\\\\\\\\', '\\\\')\n", - "colours = cosmo['colour']\n", - "s8_mean = cosmo['s8_mean'].astype(np.float64)\n", - "s8_low = cosmo['s8_low'].astype(np.float64)\n", - "s8_high = cosmo['s8_high'].astype(np.float64)\n", - "sigma8_mean = cosmo['sigma8_mean'].astype(np.float64)\n", - "sigma8_low = cosmo['sigma8_low'].astype(np.float64)\n", - "sigma8_high = cosmo['sigma8_high'].astype(np.float64)\n", - "omegam_mean = cosmo['omegam_mean'].astype(np.float64)\n", - "omegam_low = cosmo['omegam_low'].astype(np.float64)\n", - "omegam_high = cosmo['omegam_high'].astype(np.float64)" + "cosmo = np.loadtxt(\n", + " f\"{root_dir}/param_values.txt\",\n", + " dtype={\n", + " \"names\": (\n", + " \"Expt\",\n", + " \"colour\",\n", + " \"s8_mean\",\n", + " \"s8_low\",\n", + " \"s8_high\",\n", + " \"sigma8_mean\",\n", + " \"sigma8_low\",\n", + " \"sigma8_high\",\n", + " \"omegam_mean\",\n", + " \"omegam_low\",\n", + " \"omegam_high\",\n", + " ),\n", + " \"formats\": (\n", + " \"U250\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " \"U20\",\n", + " ),\n", + " },\n", + " skiprows=1,\n", + " delimiter=\";\",\n", + ")\n", + "expt = np.char.replace(cosmo[\"Expt\"], \"\\\\\\\\\", \"\\\\\")\n", + "colours = cosmo[\"colour\"]\n", + "s8_mean = cosmo[\"s8_mean\"].astype(np.float64)\n", + "s8_low = cosmo[\"s8_low\"].astype(np.float64)\n", + "s8_high = cosmo[\"s8_high\"].astype(np.float64)\n", + "sigma8_mean = cosmo[\"sigma8_mean\"].astype(np.float64)\n", + "sigma8_low = cosmo[\"sigma8_low\"].astype(np.float64)\n", + "sigma8_high = cosmo[\"sigma8_high\"].astype(np.float64)\n", + "omegam_mean = cosmo[\"omegam_mean\"].astype(np.float64)\n", + "omegam_low = cosmo[\"omegam_low\"].astype(np.float64)\n", + "omegam_high = cosmo[\"omegam_high\"].astype(np.float64)" ] }, { @@ -419,6 +454,7 @@ " sign = 1 if mean1 > mean2 else -1\n", " return sigma_tension * sign\n", "\n", + "\n", "def get_sigma_tension_asymmetric(mean1, low1, high1, mean2, low2, high2):\n", " delta_mean = mean1 - mean2\n", " if delta_mean > 0:\n", @@ -467,7 +503,7 @@ " r\"No leakage correction\",\n", " r\"No $m$-bias correction\",\n", " r\"MCCD PSF model\",\n", - " r\"UNIONS $C_\\ell$, Blind C\"\n", + " r\"UNIONS $C_\\ell$, Blind C\",\n", "]\n", "list_section_index = [\n", " r\"(ii)\",\n", @@ -478,7 +514,7 @@ " r\"(vii)\",\n", " r\"(viii)\",\n", " r\"(ix)\",\n", - " r\"(x)\"\n", + " r\"(x)\",\n", "]\n", "\n", "preliminary_watermark = False\n", @@ -491,88 +527,168 @@ "for ax, param in zip(axs, params):\n", " means, lows, highs, label = param\n", " for i, mean, low, high, color in zip(y, means, lows, highs, colours):\n", - " ax.errorbar(mean, 0.05+i*row_spacing, xerr=np.array([low, high])[:, None], fmt='o', color=color, ecolor=color, elinewidth=2, capsize=3)\n", + " ax.errorbar(\n", + " mean,\n", + " 0.05 + i * row_spacing,\n", + " xerr=np.array([low, high])[:, None],\n", + " fmt=\"o\",\n", + " color=color,\n", + " ecolor=color,\n", + " elinewidth=2,\n", + " capsize=3,\n", + " )\n", " ax.set_xlabel(label, fontsize=14)\n", - " \n", + "\n", " ax.grid(False)\n", - " ax.tick_params(axis='y', left=False, labelleft=False)\n", + " ax.tick_params(axis=\"y\", left=False, labelleft=False)\n", " if label == r\"$S_8$\":\n", - " ax.axvspan(s8_mean[index_ref] - s8_low[index_ref], s8_mean[index_ref] + s8_high[index_ref], color=colours[index_ref], alpha=0.2)\n", + " ax.axvspan(\n", + " s8_mean[index_ref] - s8_low[index_ref],\n", + " s8_mean[index_ref] + s8_high[index_ref],\n", + " color=colours[index_ref],\n", + " alpha=0.2,\n", + " )\n", " ax.set_xlim(-0.1, 1.41)\n", " if blind_axes:\n", " ref_tick = np.mean(s8_mean[:4])\n", - " ax.set_xticks(\n", - " [ref_tick + i*0.1 for i in range(-5, 5)], labels=[]\n", - " )\n", + " ax.set_xticks([ref_tick + i * 0.1 for i in range(-5, 5)], labels=[])\n", " else:\n", " ax.set_xticks(\n", " [0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4],\n", " labels=[\n", - " r\"$0.6$\", r\"\", r\"$0.8$\", r\"\", r\"$1.0$\", r\"\", r\"$1.2$\", r\"\", r\"$1.4$\"\n", - " ]\n", + " r\"$0.6$\",\n", + " r\"\",\n", + " r\"$0.8$\",\n", + " r\"\",\n", + " r\"$1.0$\",\n", + " r\"\",\n", + " r\"$1.2$\",\n", + " r\"\",\n", + " r\"$1.4$\",\n", + " ],\n", " )\n", " elif label == r\"$\\sigma_8$\":\n", - " ax.axvspan(sigma8_mean[index_ref] - sigma8_low[index_ref], sigma8_mean[index_ref] + sigma8_high[index_ref], color=colours[index_ref], alpha=0.2)\n", + " ax.axvspan(\n", + " sigma8_mean[index_ref] - sigma8_low[index_ref],\n", + " sigma8_mean[index_ref] + sigma8_high[index_ref],\n", + " color=colours[index_ref],\n", + " alpha=0.2,\n", + " )\n", " ax.set_xlim(0.5, 1.25)\n", " if blind_axes:\n", " ref_tick = np.mean(sigma8_mean[:4])\n", - " ax.set_xticks(\n", - " [ref_tick + i*0.2 for i in range(-2, 2)], labels=[]\n", - " )\n", + " ax.set_xticks([ref_tick + i * 0.2 for i in range(-2, 2)], labels=[])\n", " elif label == r\"$\\Omega_{\\rm m}$\":\n", - " ax.axvspan(omegam_mean[index_ref] - omegam_low[index_ref], omegam_mean[index_ref] + omegam_high[index_ref], color=colours[index_ref], alpha=0.2)\n", + " ax.axvspan(\n", + " omegam_mean[index_ref] - omegam_low[index_ref],\n", + " omegam_mean[index_ref] + omegam_high[index_ref],\n", + " color=colours[index_ref],\n", + " alpha=0.2,\n", + " )\n", " ax.set_xlim(0.1, 0.5)\n", " if blind_axes:\n", " ref_tick = np.mean(omegam_mean[:4])\n", - " ax.set_xticks(\n", - " [ref_tick + i*0.1 for i in range(-2, 3)], labels=[]\n", - " )\n", + " ax.set_xticks([ref_tick + i * 0.1 for i in range(-2, 3)], labels=[])\n", "\n", "\n", - "axs[0].set_yticks(0.05+y*row_spacing)\n", + "axs[0].set_yticks(0.05 + y * row_spacing)\n", "axs[0].set_yticklabels([])\n", "for label, color in zip(expt, colours):\n", - " axs[0].text(-0.09, 0.05 + row_spacing * np.where(expt == label)[0][0], label, fontsize=11, ha='left', va='center', color=color)\n", + " axs[0].text(\n", + " -0.09,\n", + " 0.05 + row_spacing * np.where(expt == label)[0][0],\n", + " label,\n", + " fontsize=11,\n", + " ha=\"left\",\n", + " va=\"center\",\n", + " color=color,\n", + " )\n", " if label != reference:\n", " index = np.where(expt == label)[0][0]\n", " s8_tension = tension_function(\n", - " s8_mean[index], s8_low[index], s8_high[index],\n", - " s8_mean[index_ref], s8_low[index_ref], s8_high[index_ref]\n", + " s8_mean[index],\n", + " s8_low[index],\n", + " s8_high[index],\n", + " s8_mean[index_ref],\n", + " s8_low[index_ref],\n", + " s8_high[index_ref],\n", " )\n", " sign_str = \"+\" if s8_tension > 0 else \"-\"\n", " text_to_plot = rf\"${sign_str}{np.abs(s8_tension):.2f}\" + r\"\\, \\sigma$\"\n", " if label == r\"UNIONS $\\xi_\\pm(\\vartheta)$ (Goh et al., 2026)\":\n", " text_to_plot += r\" (ind.)\"\n", - " axs[0].text(1.38, 0.05 + row_spacing * (index - 0.45), text_to_plot, fontsize=10, ha='right', va='center', color=color)\n", + " axs[0].text(\n", + " 1.38,\n", + " 0.05 + row_spacing * (index - 0.45),\n", + " text_to_plot,\n", + " fontsize=10,\n", + " ha=\"right\",\n", + " va=\"center\",\n", + " color=color,\n", + " )\n", " additional_text = r\"$-2.18 \\, \\sigma$ (corr.)\"\n", - " axs[0].text(1.4, 0.05 + row_spacing * (index + 0.15), additional_text, fontsize=10, ha='right', va='center', color=color)\n", + " axs[0].text(\n", + " 1.4,\n", + " 0.05 + row_spacing * (index + 0.15),\n", + " additional_text,\n", + " fontsize=10,\n", + " ha=\"right\",\n", + " va=\"center\",\n", + " color=color,\n", + " )\n", " else:\n", - " axs[0].text(1.38, 0.05 + row_spacing * index, text_to_plot, fontsize=10, ha='right', va='center', color=color)\n", + " axs[0].text(\n", + " 1.38,\n", + " 0.05 + row_spacing * index,\n", + " text_to_plot,\n", + " fontsize=10,\n", + " ha=\"right\",\n", + " va=\"center\",\n", + " color=color,\n", + " )\n", "# Add separation lines\n", "for i, sep in enumerate(separation_after):\n", " index_sep = np.where(expt == sep)[0][0]\n", " for ax in axs:\n", - " ax.axhline(row_spacing * (index_sep + 0.95), color='black', linestyle='dotted', linewidth=1)\n", - " axs[0].text(-0.1, 0.05 + row_spacing * (index_sep + 1), \n", - " list_section_index[i], fontsize=14, fontweight='bold', va='center', ha='right')\n", + " ax.axhline(\n", + " row_spacing * (index_sep + 0.95),\n", + " color=\"black\",\n", + " linestyle=\"dotted\",\n", + " linewidth=1,\n", + " )\n", + " axs[0].text(\n", + " -0.1,\n", + " 0.05 + row_spacing * (index_sep + 1),\n", + " list_section_index[i],\n", + " fontsize=14,\n", + " fontweight=\"bold\",\n", + " va=\"center\",\n", + " ha=\"right\",\n", + " )\n", "\n", "\n", "# --- Add section labels (i), (ii)) ---\n", - "axs[0].text(-0.1, 0.05, \n", - " r\"(i)\", fontsize=14, fontweight='bold', va='center', ha='right')\n", + "axs[0].text(-0.1, 0.05, r\"(i)\", fontsize=14, fontweight=\"bold\", va=\"center\", ha=\"right\")\n", "\n", "if preliminary_watermark:\n", - " plt.figtext(0.5, 0.5, 'PRELIMINARY',\n", - " fontsize=50, color='gray',\n", - " ha='center', va='center',\n", - " alpha=0.3, rotation=330)\n", + " plt.figtext(\n", + " 0.5,\n", + " 0.5,\n", + " \"PRELIMINARY\",\n", + " fontsize=50,\n", + " color=\"gray\",\n", + " ha=\"center\",\n", + " va=\"center\",\n", + " alpha=0.3,\n", + " rotation=330,\n", + " )\n", "\n", "plt.gca().invert_yaxis()\n", "\n", "plt.tight_layout()\n", "\n", "plt.savefig(\"./plots/whisker_plot.png\", dpi=300)\n", - "#Save pdf\n", + "# Save pdf\n", "plt.savefig(\"./plots/whisker_plot.pdf\")\n", "plt.show()" ]