From 39bfb328b5739c04ff3ab286a8f8b632945f0a24 Mon Sep 17 00:00:00 2001 From: Michel Date: Fri, 20 Sep 2024 17:01:10 +0200 Subject: [PATCH 01/17] add Paper_v2.0/paper_formating.py --- examples/Paper_v2.0/paper_formating.py | 423 +++++++++++++++++++++++++ 1 file changed, 423 insertions(+) create mode 100644 examples/Paper_v2.0/paper_formating.py diff --git a/examples/Paper_v2.0/paper_formating.py b/examples/Paper_v2.0/paper_formating.py new file mode 100644 index 000000000..36314e1ee --- /dev/null +++ b/examples/Paper_v2.0/paper_formating.py @@ -0,0 +1,423 @@ +import matplotlib as mpl +import pylab as plt + + +def prep_plot( + figsize=(9, 9), + adjust=None, # {'left':0.2, 'right':0.95, 'top':0.94, 'bottom':0.2}, + subplots=False, +): + main_parameters = { + "figure.dpi": 300.0, + "legend.fontsize": 10, + "axes.labelsize": 12, + "axes.titlesize": 12, + "xtick.labelsize": 12, + "ytick.labelsize": 12, + "axes.linewidth": 0.8, + "xtick.major.width": 0.8, + "ytick.major.width": 0.8, + "errorbar.capsize": 0.8, + #'errorbar.capthick': 0.4, + #'errorbar.elinewidth': 0.4, + } + other_parameters = { + "legend": { + "fontsize": 10, + }, + "xlabel": {}, # {'fontsize': 2,}, + "ylabel": {}, # {'fontsize': 20,}, + "xtick": {"labelsize": 10}, + "ytick": {"labelsize": 10}, + } + plt.clf() + _figsize = list(x / 2.54 for x in figsize) + print(_figsize) + if subplots: + try: + fig, axes = plt.subplots(*subplots, figsize=_figsize) + except: + fig, axes = plt.subplots(subplots, figsize=_figsize) + else: + fig = plt.figure(figsize=_figsize) + if adjust is not None: + plt.subplots_adjust(**adjust) + mpl.rcParams.update(main_parameters) + # mpl.rc('xtick', **other_parameters['xtick']) + # mpl.rc('ytick', **other_parameters['ytick']) + + if subplots: + return fig, axes + return fig + + +default = { + "_internal.classic_mode": False, + "agg.path.chunksize": 0, + "animation.avconv_args": [], + "animation.avconv_path": "avconv", + "animation.bitrate": -1, + "animation.codec": "h264", + "animation.convert_args": [], + "animation.convert_path": "convert", + "animation.embed_limit": 20.0, + "animation.ffmpeg_args": [], + "animation.ffmpeg_path": "ffmpeg", + "animation.frame_format": "png", + "animation.html": "none", + "animation.html_args": [], + "animation.writer": "ffmpeg", + "axes.autolimit_mode": "data", + "axes.axisbelow": "line", + "axes.edgecolor": "black", + "axes.facecolor": "white", + "axes.formatter.limits": [-7, 7], + "axes.formatter.min_exponent": 0, + "axes.formatter.offset_threshold": 4, + "axes.formatter.use_locale": False, + "axes.formatter.use_mathtext": False, + "axes.formatter.useoffset": True, + "axes.grid": False, + "axes.grid.axis": "both", + "axes.grid.which": "major", + "axes.labelcolor": "black", + "axes.labelpad": 4.0, + "axes.labelsize": 20.0, + "axes.labelweight": "normal", + "axes.linewidth": 0.8, + "axes.prop_cycle": mpl.cycler( + "color", + [ + "#1f77b4", + "#ff7f0e", + "#2ca02c", + "#d62728", + "#9467bd", + "#8c564b", + "#e377c2", + "#7f7f7f", + "#bcbd22", + "#17becf", + ], + ), + "axes.spines.bottom": True, + "axes.spines.left": True, + "axes.spines.right": True, + "axes.spines.top": True, + "axes.titlepad": 6.0, + "axes.titlesize": 20.0, + "axes.titleweight": "normal", + "axes.unicode_minus": True, + "axes.xmargin": 0.05, + "axes.ymargin": 0.05, + "axes3d.grid": True, + "backend": "module://ipykernel.pylab.backend_inline", + "backend.qt4": None, + "backend.qt5": None, + "backend_fallback": True, + "boxplot.bootstrap": None, + "boxplot.boxprops.color": "black", + "boxplot.boxprops.linestyle": "-", + "boxplot.boxprops.linewidth": 1.0, + "boxplot.capprops.color": "black", + "boxplot.capprops.linestyle": "-", + "boxplot.capprops.linewidth": 1.0, + "boxplot.flierprops.color": "black", + "boxplot.flierprops.linestyle": "none", + "boxplot.flierprops.linewidth": 1.0, + "boxplot.flierprops.marker": "o", + "boxplot.flierprops.markeredgecolor": "black", + "boxplot.flierprops.markerfacecolor": "none", + "boxplot.flierprops.markersize": 6.0, + "boxplot.meanline": False, + "boxplot.meanprops.color": "C2", + "boxplot.meanprops.linestyle": "--", + "boxplot.meanprops.linewidth": 1.0, + "boxplot.meanprops.marker": "^", + "boxplot.meanprops.markeredgecolor": "C2", + "boxplot.meanprops.markerfacecolor": "C2", + "boxplot.meanprops.markersize": 6.0, + "boxplot.medianprops.color": "C1", + "boxplot.medianprops.linestyle": "-", + "boxplot.medianprops.linewidth": 1.0, + "boxplot.notch": False, + "boxplot.patchartist": False, + "boxplot.showbox": True, + "boxplot.showcaps": True, + "boxplot.showfliers": True, + "boxplot.showmeans": False, + "boxplot.vertical": True, + "boxplot.whiskerprops.color": "black", + "boxplot.whiskerprops.linestyle": "-", + "boxplot.whiskerprops.linewidth": 1.0, + "boxplot.whiskers": 1.5, + "contour.corner_mask": True, + "contour.negative_linestyle": "dashed", + "datapath": "/opt/conda/lib/python3.6/site-packages/matplotlib/mpl-data", + "date.autoformatter.day": "%Y-%m-%d", + "date.autoformatter.hour": "%m-%d %H", + "date.autoformatter.microsecond": "%M:%S.%f", + "date.autoformatter.minute": "%d %H:%M", + "date.autoformatter.month": "%Y-%m", + "date.autoformatter.second": "%H:%M:%S", + "date.autoformatter.year": "%Y", + "docstring.hardcopy": False, + "errorbar.capsize": 0.0, + "examples.directory": "", + "figure.autolayout": False, + "figure.constrained_layout.h_pad": 0.04167, + "figure.constrained_layout.hspace": 0.02, + "figure.constrained_layout.use": False, + "figure.constrained_layout.w_pad": 0.04167, + "figure.constrained_layout.wspace": 0.02, + "figure.dpi": 100.0, + "figure.edgecolor": "white", + "figure.facecolor": "white", + "figure.figsize": [6.4, 4.8], + "figure.frameon": True, + "figure.max_open_warning": 20, + "figure.subplot.bottom": 0.11, + "figure.subplot.hspace": 0.2, + "figure.subplot.left": 0.125, + "figure.subplot.right": 0.9, + "figure.subplot.top": 0.88, + "figure.subplot.wspace": 0.2, + "figure.titlesize": "large", + "figure.titleweight": "normal", + "font.cursive": [ + "Apple Chancery", + "Textile", + "Zapf Chancery", + "Sand", + "Script MT", + "Felipa", + "cursive", + ], + "font.family": ["sans-serif"], + "font.fantasy": [ + "Comic Sans MS", + "Chicago", + "Charcoal", + "Impact", + "Western", + "Humor Sans", + "xkcd", + "fantasy", + ], + "font.monospace": [ + "DejaVu Sans Mono", + "Bitstream Vera Sans Mono", + "Computer Modern Typewriter", + "Andale Mono", + "Nimbus Mono L", + "Courier New", + "Courier", + "Fixed", + "Terminal", + "monospace", + ], + "font.sans-serif": [ + "DejaVu Sans", + "Bitstream Vera Sans", + "Computer Modern Sans Serif", + "Lucida Grande", + "Verdana", + "Geneva", + "Lucid", + "Arial", + "Helvetica", + "Avant Garde", + "sans-serif", + ], + "font.serif": [ + "DejaVu Serif", + "Bitstream Vera Serif", + "Computer Modern Roman", + "New Century Schoolbook", + "Century Schoolbook L", + "Utopia", + "ITC Bookman", + "Bookman", + "Nimbus Roman No9 L", + "Times New Roman", + "Times", + "Palatino", + "Charter", + "serif", + ], + "font.size": 10.0, + "font.stretch": "normal", + "font.style": "normal", + "font.variant": "normal", + "font.weight": "normal", + "grid.alpha": 1.0, + "grid.color": "#b0b0b0", + "grid.linestyle": "-", + "grid.linewidth": 0.8, + "hatch.color": "black", + "hatch.linewidth": 1.0, + "hist.bins": 10, + "image.aspect": "equal", + "image.cmap": "viridis", + "image.composite_image": True, + "image.interpolation": "nearest", + "image.lut": 256, + "image.origin": "upper", + "image.resample": True, + "interactive": True, + "keymap.all_axes": ["a"], + "keymap.back": ["left", "c", "backspace"], + "keymap.copy": ["ctrl+c", "cmd+c"], + "keymap.forward": ["right", "v"], + "keymap.fullscreen": ["f", "ctrl+f"], + "keymap.grid": ["g"], + "keymap.grid_minor": ["G"], + "keymap.help": ["f1"], + "keymap.home": ["h", "r", "home"], + "keymap.pan": ["p"], + "keymap.quit": ["ctrl+w", "cmd+w", "q"], + "keymap.quit_all": ["W", "cmd+W", "Q"], + "keymap.save": ["s", "ctrl+s"], + "keymap.xscale": ["k", "L"], + "keymap.yscale": ["l"], + "keymap.zoom": ["o"], + "legend.borderaxespad": 0.5, + "legend.borderpad": 0.4, + "legend.columnspacing": 2.0, + "legend.edgecolor": "0.8", + "legend.facecolor": "inherit", + "legend.fancybox": True, + "legend.fontsize": 20.0, + "legend.framealpha": 0.8, + "legend.frameon": True, + "legend.handleheight": 0.7, + "legend.handlelength": 2.0, + "legend.handletextpad": 0.8, + "legend.labelspacing": 0.5, + "legend.loc": "best", + "legend.markerscale": 1.0, + "legend.numpoints": 1, + "legend.scatterpoints": 1, + "legend.shadow": False, + "legend.title_fontsize": None, + "lines.antialiased": True, + "lines.color": "C0", + "lines.dash_capstyle": "butt", + "lines.dash_joinstyle": "round", + "lines.dashdot_pattern": [6.4, 1.6, 1.0, 1.6], + "lines.dashed_pattern": [3.7, 1.6], + "lines.dotted_pattern": [1.0, 1.65], + "lines.linestyle": "-", + "lines.linewidth": 1.5, + "lines.marker": "None", + "lines.markeredgecolor": "auto", + "lines.markeredgewidth": 1.0, + "lines.markerfacecolor": "auto", + "lines.markersize": 6.0, + "lines.scale_dashes": True, + "lines.solid_capstyle": "projecting", + "lines.solid_joinstyle": "round", + "markers.fillstyle": "full", + "mathtext.bf": "sans:bold", + "mathtext.cal": "cursive", + "mathtext.default": "it", + "mathtext.fallback_to_cm": True, + "mathtext.fontset": "dejavusans", + "mathtext.it": "sans:italic", + "mathtext.rm": "sans", + "mathtext.sf": "sans", + "mathtext.tt": "monospace", + "patch.antialiased": True, + "patch.edgecolor": "black", + "patch.facecolor": "C0", + "patch.force_edgecolor": False, + "patch.linewidth": 1.0, + "path.effects": [], + "path.simplify": True, + "path.simplify_threshold": 0.1111111111111111, + "path.sketch": None, + "path.snap": True, + "pdf.compression": 6, + "pdf.fonttype": 3, + "pdf.inheritcolor": False, + "pdf.use14corefonts": False, + "pgf.preamble": [], + "pgf.rcfonts": True, + "pgf.texsystem": "xelatex", + "polaraxes.grid": True, + "ps.distiller.res": 6000, + "ps.fonttype": 3, + "ps.papersize": "letter", + "ps.useafm": False, + "ps.usedistiller": False, + "savefig.bbox": None, + "savefig.directory": "~", + "savefig.dpi": "figure", + "savefig.edgecolor": "white", + "savefig.facecolor": "white", + "savefig.format": "png", + "savefig.frameon": True, + "savefig.jpeg_quality": 95, + "savefig.orientation": "portrait", + "savefig.pad_inches": 0.1, + "savefig.transparent": False, + "scatter.marker": "o", + "svg.fonttype": "path", + "svg.hashsalt": None, + "svg.image_inline": True, + "text.antialiased": True, + "text.color": "black", + "text.hinting": "auto", + "text.hinting_factor": 8, + "text.latex.preamble": [], + "text.latex.preview": False, + "text.latex.unicode": True, + "text.usetex": False, + "timezone": "UTC", + "tk.window_focus": False, + "toolbar": "toolbar2", + "verbose.fileo": "sys.stdout", + "verbose.level": "silent", + "webagg.address": "127.0.0.1", + "webagg.open_in_browser": True, + "webagg.port": 8988, + "webagg.port_retries": 50, + "xtick.alignment": "center", + "xtick.bottom": True, + "xtick.color": "black", + "xtick.direction": "out", + "xtick.labelbottom": True, + "xtick.labelsize": 18.0, + "xtick.labeltop": False, + "xtick.major.bottom": True, + "xtick.major.pad": 3.5, + "xtick.major.size": 3.5, + "xtick.major.top": True, + "xtick.major.width": 0.8, + "xtick.minor.bottom": True, + "xtick.minor.pad": 3.4, + "xtick.minor.size": 2.0, + "xtick.minor.top": True, + "xtick.minor.visible": False, + "xtick.minor.width": 0.6, + "xtick.top": False, + "ytick.alignment": "center_baseline", + "ytick.color": "black", + "ytick.direction": "out", + "ytick.labelleft": True, + "ytick.labelright": False, + "ytick.labelsize": 18.0, + "ytick.left": True, + "ytick.major.left": True, + "ytick.major.pad": 3.5, + "ytick.major.right": True, + "ytick.major.size": 3.5, + "ytick.major.width": 0.8, + "ytick.minor.left": True, + "ytick.minor.pad": 3.4, + "ytick.minor.right": True, + "ytick.minor.size": 2.0, + "ytick.minor.visible": False, + "ytick.minor.width": 0.6, + "ytick.right": False, +} From cf18c6ca6ed5204849f0c5458940e6ea4e0a1b9c Mon Sep 17 00:00:00 2001 From: Celine Combet Date: Thu, 30 Jan 2025 12:12:42 +0100 Subject: [PATCH 02/17] NB for v2.0 paper figure showing miscentering, boost and 2-halo term --- .../fig_2halo_and_miscentering.ipynb | 204 ++++++++++++++++++ 1 file changed, 204 insertions(+) create mode 100644 examples/Paper_v2.0/fig_2halo_and_miscentering.ipynb diff --git a/examples/Paper_v2.0/fig_2halo_and_miscentering.ipynb b/examples/Paper_v2.0/fig_2halo_and_miscentering.ipynb new file mode 100644 index 000000000..9c63f2915 --- /dev/null +++ b/examples/Paper_v2.0/fig_2halo_and_miscentering.ipynb @@ -0,0 +1,204 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Figure to show the new 2-halo term and miscentering functionality" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Imports specific to clmm " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "os.environ[\n", + " \"CLMM_MODELING_BACKEND\"\n", + "] = \"ccl\" # here you may choose ccl, nc (NumCosmo) or ct (cluster_toolkit)\n", + "\n", + "import clmm\n", + "from clmm import Cosmology\n", + "import clmm.utils as u" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make sure we know which version we're using" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clmm.__version__" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define a cosmology using astropy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cosmo = Cosmology(H0=67.0, Omega_dm0=0.315 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the galaxy cluster model. Here, we choose parameters that describe the galaxy cluster model, including the mass definition, concentration, and mass distribution. For the mass distribution, we choose a distribution that follows an NFW profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "moo = clmm.Modeling(massdef=\"mean\", delta_mdef=200, halo_profile_model=\"nfw\")\n", + "\n", + "mass_cl = 1.e14\n", + "z_cl = 0.4\n", + "\n", + "conc_cl = 5.4# Duffy08 value for this halo mass and redshift (see last commented cell)\n", + "halo_bias = 2.4 # Tkinker10 value for this halo mass and redshift (see last commented cell)\n", + "\n", + "moo.set_cosmo(cosmo)\n", + "moo.set_concentration(conc_cl)\n", + "moo.set_mass(mass_cl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r_proj = np.logspace(-2, 2, 100)\n", + "\n", + "DeltaSigma = moo.eval_excess_surface_density(r_proj, z_cl)\n", + "\n", + "# Miscentered DeltaSigma\n", + "DeltaSigma_mis = moo.eval_excess_surface_density(r_proj, z_cl, r_mis=0.2)\n", + "\n", + "# 2halo DeltaSigma\n", + "DeltaSigma_2h = moo.eval_excess_surface_density_2h(r_proj, z_cl, halobias=6)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r_scale = 0.3\n", + "nfw_boost = u.compute_nfw_boost(r_proj, r_scale, boost0=0.2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the predicted profiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7,6))\n", + "\n", + "ax.loglog(r_proj, DeltaSigma,label='1-halo, centered, reference', color='k')\n", + "ax.loglog(r_proj, DeltaSigma/nfw_boost, ls=':', label='1-halo, corrected with NFW boost model', color='grey')\n", + "ax.loglog(r_proj, DeltaSigma_mis , ls ='--', label='1-halo, miscentered, R_mis = 0.2 Mpc')\n", + "ax.loglog(r_proj, DeltaSigma_2h, ls='dashdot', label='2-halo', color='green')\n", + "ax.legend(loc=1, fontsize=13)\n", + "\n", + "ax.set_xlabel('R [Mpc]', fontsize=15)\n", + "ax.set_ylabel(r'$\\Delta\\Sigma$ [M$_\\odot$ Mpc$^{-2}$]', fontsize=15)\n", + "ax.tick_params(axis='x', labelsize=15)\n", + "ax.tick_params(axis='y', labelsize=15)\n", + "ax.set_ylim([8.e10, 1e16])\n", + "\n", + "fig.tight_layout()\n", + "fig.savefig('2h_miscentering_boost.png')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Typical values for halo_bias and concentration to use in the the example, computed using CCL\n", + "\n", + "import pyccl as ccl\n", + "\n", + "cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.83, n_s=0.96)\n", + "\n", + "a = 1./(1+z_cl)\n", + "\n", + "bias = ccl.halos.HaloBiasTinker10(mass_def='200m', mass_def_strict=True)\n", + "conc = ccl.halos.concentration.ConcentrationDuffy08(mass_def='200m')\n", + "\n", + "bias(cosmo, mass_cl, a), conc(cosmo, mass_cl, a)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "wrk", + "language": "python", + "name": "wrk" + }, + "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.10.9" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 7e3669bbf6faa219989a550a99ab90f20b7b0db5 Mon Sep 17 00:00:00 2001 From: Celine Combet Date: Thu, 30 Jan 2025 12:14:46 +0100 Subject: [PATCH 03/17] Rename NB --- ...scentering.ipynb => fig_2halo_miscentering_boost_theory.ipynb} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/Paper_v2.0/{fig_2halo_and_miscentering.ipynb => fig_2halo_miscentering_boost_theory.ipynb} (100%) diff --git a/examples/Paper_v2.0/fig_2halo_and_miscentering.ipynb b/examples/Paper_v2.0/fig_2halo_miscentering_boost_theory.ipynb similarity index 100% rename from examples/Paper_v2.0/fig_2halo_and_miscentering.ipynb rename to examples/Paper_v2.0/fig_2halo_miscentering_boost_theory.ipynb From 9a0593ef36d7c404e34b91d75536d92ac517b999 Mon Sep 17 00:00:00 2001 From: Radhakrishnan Srinivasan <47119549+rkrishnan2912@users.noreply.github.com> Date: Fri, 30 May 2025 12:30:22 -0400 Subject: [PATCH 04/17] Notebook that makes Figure 7 --- ...xial_excess_surface_density_profiles.ipynb | 156 ++++++++++++++++++ 1 file changed, 156 insertions(+) create mode 100644 examples/Paper_v2.0/Figure_7_triaxial_excess_surface_density_profiles.ipynb diff --git a/examples/Paper_v2.0/Figure_7_triaxial_excess_surface_density_profiles.ipynb b/examples/Paper_v2.0/Figure_7_triaxial_excess_surface_density_profiles.ipynb new file mode 100644 index 000000000..ab55dc86a --- /dev/null +++ b/examples/Paper_v2.0/Figure_7_triaxial_excess_surface_density_profiles.ipynb @@ -0,0 +1,156 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "7c213181-875a-408b-a95f-19dac19a324e", + "metadata": {}, + "source": [ + "# To make a plot that compares the contribution to excess surface density profiles for weak lensing when the lens is triaxial" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "43f767fb-627e-4c15-aede-d7a7b051f9bb", + "metadata": {}, + "outputs": [], + "source": [ + "import clmm\n", + "from clmm import Cosmology\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import seaborn as sns\n", + "\n", + "cosmo = Cosmology(H0=70.0, Omega_dm0=0.2248, Omega_b0=0.3 - 0.2248, Omega_k0=0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "583e00e3-2f75-4a10-a773-ee816be0514a", + "metadata": {}, + "outputs": [], + "source": [ + "r=np.linspace(0.1, 5, 100)\n", + "\n", + "# A sample Mass, concentration for a halo:\n", + "mdelta_fit = 2.0E14 # in M_sun/h\n", + "cdelta_fit = 3.89 \n", + "\n", + "# Ellipticities:\n", + "ell_1_0 = 0.0 # q = 1.0, ellipticity is (1-q)/(1+q)\n", + "ell_0_7 = 0.176 # q = 0.7, ellipticity is (1-q)/(1+q)\n", + "ell_0_3 = 0.497 # q = 0.3, ellipticity is (1-q)/(1+q)\n", + "\n", + "z_cl = 0.47\n", + "\n", + "mdef = 'critical'\n", + "\n", + "'''\n", + "A Spherical Halo: (Axis ratio, q = 1.0 )\n", + "'''\n", + "ds_model_e_0 = clmm.compute_excess_surface_density_triaxial(ell=ell_1_0, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='mono')\n", + "ds_const_model_e_0 = clmm.compute_excess_surface_density_triaxial(ell=ell_1_0, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='quad_const')\n", + "ds_4theta_model_e_0 = clmm.compute_excess_surface_density_triaxial(ell=ell_1_0, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='quad_4theta')\n", + "\n", + "'''\n", + "An Elliptical Halo: (Axis ratio, q = 0.7)\n", + "'''\n", + "ds_model = clmm.compute_excess_surface_density_triaxial(ell=ell_0_7, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, z_cl=z_cl, \n", + " cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='mono')\n", + "ds_const_model = clmm.compute_excess_surface_density_triaxial(ell=ell_0_7, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='quad_const')\n", + "ds_4theta_model = clmm.compute_excess_surface_density_triaxial(ell=ell_0_7, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='quad_4theta')\n", + "\n", + "\n", + "'''\n", + "An Elliptical Halo: (Axis ratio, q = 0.3)\n", + "'''\n", + "ds_model_0_3 = clmm.compute_excess_surface_density_triaxial(ell=ell_0_3, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='mono')\n", + "ds_const_model_0_3 = clmm.compute_excess_surface_density_triaxial(ell=ell_0_3, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='quad_const')\n", + "ds_4theta_model_0_3 = clmm.compute_excess_surface_density_triaxial(ell=ell_0_3, r_proj=r, mdelta=mdelta_fit, cdelta=cdelta_fit, \n", + " z_cl=z_cl, cosmo=cosmo, halo_profile_model='nfw', \n", + " delta_mdef=200, massdef=mdef, term='quad_4theta')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b87d23bd-7432-42cc-aee4-addbdd6f90c0", + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(1,2, figsize=[10,5])\n", + "\n", + "sns.lineplot(x=r, y=ds_model_e_0/1E12, ax=ax[0], markers = True, ms=20.0, color='black', ls='dashdot', label='Spherical NFW')\n", + "sns.lineplot(x=r, y=ds_model/1E12, ax=ax[0], color='darkorange', label='Elliptical NFW (q = 2/3)', ls='solid')\n", + "sns.lineplot(x=r, y=ds_model_0_3/1E12, ax=ax[0], color='darkorange', label='Elliptical NFW (q = 1/3)', ls='solid', alpha=0.4)\n", + "ax[0].set_yscale('log')\n", + "ax[0].set_xscale('log')\n", + "ax[0].set_xlabel('r [Mpc]', fontsize=18)\n", + "ax[0].set_ylabel(r'Monopole $\\Delta\\Sigma \\,[M_{\\odot} h/\\rm{pc}^{2}]$', fontsize=18)\n", + "ax[0].tick_params(axis='both', which='both', labelsize=12)\n", + "ax[0].legend(fontsize=11)\n", + "plt.tight_layout()\n", + "axin1 = ax[0].inset_axes(\n", + " [0.18, 7, 1, 20], transform=ax[0].transData)\n", + "\n", + "sns.lineplot(x=r,y=(1-ds_model/ds_model_e_0)*100, ax=axin1, color='crimson', label='q=2/3')\n", + "sns.lineplot(x=r,y=(1-ds_model_0_3/ds_model_e_0)*100, ax=axin1, color='crimson', alpha=0.4, label='q=1/3')\n", + "axin1.set_xscale('log')\n", + "axin1.set_ylabel('Difference (%)', fontsize=10)\n", + "axin1.legend()\n", + "#axin1.set_xlabel('r', fontsize=15)\n", + "\n", + "sns.lineplot(x=r, y=ds_const_model_e_0/1E12, ax=ax[1], color='black', ls='dashdot', label='Spherical NFW')\n", + "sns.lineplot(x=r, y=ds_const_model/1E12, ax=ax[1], color='crimson', label=r'Elliptical NFW (q = 2/3) | $4\\theta$')\n", + "sns.lineplot(x=r, y=ds_4theta_model/1E12, ax=ax[1], color='darkorange', label='Elliptical NFW (q = 2/3) | const', ls='dashed')\n", + "\n", + "sns.lineplot(x=r, y=ds_const_model_0_3/1E12, ax=ax[1], color='crimson', label=r'Elliptical NFW (q = 1/3) | $4\\theta$', alpha=0.4)\n", + "sns.lineplot(x=r, y=ds_4theta_model_0_3/1E12, ax=ax[1], color='darkorange', label='Elliptical NFW (q = 1/3) | const', ls='dashed', alpha=0.4)\n", + "#ax[1].set_yscale('log')\n", + "ax[1].set_xscale('log')\n", + "ax[1].set_xlabel('r [Mpc]', fontsize=18)\n", + "ax[1].set_ylabel(r'Quadrupole $\\Delta\\Sigma \\,[M_{\\odot} h/\\rm{pc}^{2}]$', fontsize=18)\n", + "ax[1].tick_params(axis='both', which='both', labelsize=12)\n", + "ax[1].legend(fontsize=11)\n", + "plt.tight_layout()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "base3", + "language": "python", + "name": "base3" + }, + "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.12.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From dd2a76f2413f07ce6d5261767782229f5b58fb4f Mon Sep 17 00:00:00 2001 From: Michel Date: Thu, 5 Jun 2025 20:20:36 +0200 Subject: [PATCH 05/17] add mass_conversion.ipynb --- examples/Paper_v2.0/mass_conversion.ipynb | 230 ++++++++++++++++++++++ 1 file changed, 230 insertions(+) create mode 100644 examples/Paper_v2.0/mass_conversion.ipynb diff --git a/examples/Paper_v2.0/mass_conversion.ipynb b/examples/Paper_v2.0/mass_conversion.ipynb new file mode 100644 index 000000000..0b12d1165 --- /dev/null +++ b/examples/Paper_v2.0/mass_conversion.ipynb @@ -0,0 +1,230 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Mass conversion between different mass definitions\n", + "## Mass conversion between spherical overdensity mass definitions\n", + "\n", + "In this notebook, we demonstrates how to convert the mass and concentration between various mass definitions (going from $200m$ to $500c$ in this example), and related functionalities, using both the object-oriented and functional interfaces of the code." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "os.environ[\"CLMM_MODELING_BACKEND\"] = (\n", + " \"ccl\" # here you may choose ccl, nc (NumCosmo) or ct (cluster_toolkit)\n", + ")\n", + "import clmm\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We define a cosmology first:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cosmo = clmm.Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# first SOD definition\n", + "m1 = 1e14\n", + "c1 = 4\n", + "massdef1 = \"mean\"\n", + "delta_mdef1 = 200\n", + "halo_profile_model1 = \"nfw\"\n", + "# cluster redshift\n", + "z_cl = 0.4" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# get a value for alpha\n", + "ein_def1 = clmm.Modeling(\n", + " massdef=massdef1, delta_mdef=delta_mdef1, halo_profile_model=\"einasto\"\n", + ")\n", + "ein_def1.set_mass(m1)\n", + "ein_def1.set_concentration(c1)\n", + "\n", + "alpha = ein_def1.get_einasto_alpha(z_cl)\n", + "alpha" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "var_mod = {\n", + " \"massdef2\": (\"mean\", \"critical\"),\n", + " \"delta_mdef2\": (150, 200, 500),\n", + " \"halo_profile_model2\": (\"nfw\", \"einasto\", \"hernquist\"),\n", + "}\n", + "\n", + "\n", + "def loop_keys(keys_use, current_info, cfg_dict, list_cases):\n", + " if len(keys_use) == 0:\n", + " list_cases.append(\n", + " current_info\n", + " | {\n", + " \"data\": clmm.convert_profile_mass_concentration(\n", + " mdelta=m1,\n", + " cdelta=c1,\n", + " redshift=z_cl,\n", + " cosmo=cosmo,\n", + " massdef=massdef1,\n", + " delta_mdef=delta_mdef1,\n", + " halo_profile_model=halo_profile_model1,\n", + " alpha=alpha,\n", + " **current_info\n", + " )\n", + " }\n", + " )\n", + " return\n", + "\n", + " new_keys = list(keys_use)\n", + " key = new_keys.pop(0)\n", + " for value in cfg_dict[key]:\n", + " new_info = current_info | {key: value}\n", + " loop_keys(new_keys, new_info, cfg_dict, list_cases)\n", + "\n", + "\n", + "new_defs3 = []\n", + "loop_keys(list(var_mod.keys()), {}, var_mod, new_defs3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from paper_formating import prep_plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_config = {\n", + " \"massdef2\": {\"mean\": {\"s\": 20}, \"critical\": {\"s\": 100}},\n", + " \"delta_mdef2\": {150: {\"c\": \"C0\"}, 200: {\"c\": \"C1\"}, 500: {\"c\": \"C2\"}},\n", + " \"halo_profile_model2\": {\n", + " \"nfw\": {\"marker\": \".\", \"lw\": 0.5, \"edgecolors\": \"0\", \"zorder\": 10},\n", + " \"einasto\": {\"marker\": \"^\", \"lw\": 0.5, \"edgecolors\": \"0\", \"zorder\": 2},\n", + " \"hernquist\": {\"marker\": \"v\", \"lw\": 0.5, \"edgecolors\": \"0\", \"zorder\": 2},\n", + " },\n", + "}\n", + "_leg = {\n", + " \"massdef2\": r\"$\\rho$: \",\n", + " \"delta_mdef2\": r\"$\\Delta$: \",\n", + " \"halo_profile_model2\": \"\",\n", + "}\n", + "_leg_val = {\n", + " \"nfw\": \"NFW\",\n", + " \"einasto\": \"Einasto\",\n", + " \"hernquist\": \"Hernquist\",\n", + "}\n", + "fig = prep_plot(figsize=(9, 9))\n", + "ax = plt.axes()\n", + "for _case in new_defs3:\n", + " _kwargs = {}\n", + " label = []\n", + " for key, value in _case.items():\n", + " if key != \"data\":\n", + " _kwargs.update(_config[key][value])\n", + " label.append(_leg[key] + f\"{_leg_val.get(value, value)}\")\n", + " label = \", \".join(label)\n", + " ax.scatter(*_case[\"data\"], **_kwargs, label=label)\n", + "ax.legend(ncols=3)\n", + "# ax.set_xscale(\"log\")\n", + "# \"\"\"\n", + "handles, labels = ax.get_legend_handles_labels()\n", + "\n", + "leg1 = ax.legend(\n", + " handles[:3],\n", + " [lab.split(\", \")[-1] for lab in labels[:3]],\n", + " loc=(0.38, 0.8),\n", + " fontsize=8,\n", + ")\n", + "leg2 = ax.legend(\n", + " handles[:9:3],\n", + " [lab.split(\", \")[1] for lab in labels[:9:3]],\n", + " loc=(0.015, 0.6),\n", + " fontsize=8,\n", + ")\n", + "leg3 = ax.legend(\n", + " handles[::9],\n", + " [lab.split(\", \")[0] for lab in labels[::9]],\n", + " loc=(0.015, 0.83),\n", + " fontsize=8,\n", + ")\n", + "ax.add_artist(leg1)\n", + "ax.add_artist(leg2)\n", + "# \"\"\"\n", + "ax.set_xlabel(\"$M$ [$M_\\odot$]\")\n", + "ax.set_ylabel(\"concentration\")\n", + "fit_kwargs = {\"color\": \"r\", \"zorder\": 1, \"lw\": 1} # , \"ls\": \"--\"}\n", + "ax.axvline(1e14, **fit_kwargs)\n", + "ax.axhline(4, **fit_kwargs)\n", + "ax.grid(lw=0.5)\n", + "ax.minorticks_on()\n", + "ax.grid(which=\"minor\", lw=0.1)\n", + "ax.set_axisbelow(True)\n", + "plt.tight_layout()\n", + "plt.savefig(\"mass_conversion.png\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "clmm", + "language": "python", + "name": "clmm" + }, + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From d9753e8952831489316beb16b7f32a79462ad4ab Mon Sep 17 00:00:00 2001 From: Michel Date: Thu, 5 Jun 2025 20:20:44 +0200 Subject: [PATCH 06/17] add magnification_bias.ipynb --- examples/Paper_v2.0/magnification_bias.ipynb | 189 +++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 examples/Paper_v2.0/magnification_bias.ipynb diff --git a/examples/Paper_v2.0/magnification_bias.ipynb b/examples/Paper_v2.0/magnification_bias.ipynb new file mode 100644 index 000000000..f6695d3ab --- /dev/null +++ b/examples/Paper_v2.0/magnification_bias.ipynb @@ -0,0 +1,189 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot magnification bias" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Imports specific to clmm " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "os.environ[\n", + " \"CLMM_MODELING_BACKEND\"\n", + "] = \"ccl\" # here you may choose ccl, nc (NumCosmo) or ct (cluster_toolkit)\n", + "\n", + "import clmm\n", + "from clmm import Cosmology" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make sure we know which version we're using" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clmm.__version__" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define a cosmology using astropy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Define the galaxy cluster model. Here, we choose parameters that describe the galaxy cluster model, including the mass definition, concentration, and mass distribution. For the mass distribution, we choose a distribution that follows an NFW profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "moo = clmm.Modeling(massdef=\"mean\", delta_mdef=200, halo_profile_model=\"nfw\")\n", + "\n", + "moo.set_cosmo(cosmo)\n", + "moo.set_concentration(4)\n", + "moo.set_mass(1.0e15)\n", + "\n", + "z_cl = 1.0\n", + "\n", + "# source properties\n", + "z_src = 2.0 # all sources in the same plane\n", + "z_distrib_func = clmm.redshift.distributions.chang2013 # sources redshift following a distribution\n", + "alpha = [2, -0.5]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute quantity" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "r3d = np.logspace(-2, 2, 100)\n", + "mu_bias = moo.eval_magnification_bias(r3d, z_cl, z_src, alpha)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the predicted profiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from paper_formating import prep_plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = prep_plot(figsize=(9, 9))\n", + "ax = plt.axes()\n", + "\n", + "ax.plot(\n", + " r3d, mu_bias[0] - 1, label=\"$\\\\alpha$ =\" + str(alpha[0])\n", + ")\n", + "ax.plot(r3d, mu_bias[1] - 1, label=\"$\\\\alpha$ =\" + str(alpha[1]))\n", + "\n", + "ax.set_ylabel(\"$\\delta_{\\mu}$\")\n", + "ax.set_xlabel(\"R [Mpc]\")\n", + "ax.legend()#fontsize=\"xx-large\")\n", + "ax.set_xscale(\"log\")\n", + "\n", + "ax.set_ylim(-1.1, 3)\n", + "ax.set_xlim(.06, 1e1)\n", + "\n", + "ax.grid(lw=0.5)\n", + "ax.minorticks_on()\n", + "ax.grid(which=\"minor\", lw=0.1)\n", + "ax.set_axisbelow(True)\n", + "plt.tight_layout()\n", + "plt.savefig(\"delta_mu.png\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "clmm", + "language": "python", + "name": "clmm" + }, + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From f4f5fea8e527eb1023428bc85c60fb876450cf67 Mon Sep 17 00:00:00 2001 From: Michel Date: Fri, 6 Jun 2025 10:32:54 +0200 Subject: [PATCH 07/17] add add_grid in paper_formating.py --- examples/Paper_v2.0/paper_formating.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/examples/Paper_v2.0/paper_formating.py b/examples/Paper_v2.0/paper_formating.py index 36314e1ee..2e6046dcf 100644 --- a/examples/Paper_v2.0/paper_formating.py +++ b/examples/Paper_v2.0/paper_formating.py @@ -51,6 +51,14 @@ def prep_plot( return fig +def add_grid(ax, lw_major=0.5, lw_minor=0.1): + + ax.minorticks_on() + ax.grid(lw=lw_major) + ax.grid(which="minor", lw=lw_minor) + ax.set_axisbelow(True) + + default = { "_internal.classic_mode": False, "agg.path.chunksize": 0, From 10e974a192c828bc18e22057528667847c4a50b6 Mon Sep 17 00:00:00 2001 From: Michel Date: Fri, 6 Jun 2025 10:48:42 +0200 Subject: [PATCH 08/17] format fig_2halo_miscentering_boost_theory.ipynb --- .../fig_2halo_miscentering_boost_theory.ipynb | 93 ++++++++----------- examples/Paper_v2.0/magnification_bias.ipynb | 1 + 2 files changed, 41 insertions(+), 53 deletions(-) diff --git a/examples/Paper_v2.0/fig_2halo_miscentering_boost_theory.ipynb b/examples/Paper_v2.0/fig_2halo_miscentering_boost_theory.ipynb index 9c63f2915..8a6a2adda 100644 --- a/examples/Paper_v2.0/fig_2halo_miscentering_boost_theory.ipynb +++ b/examples/Paper_v2.0/fig_2halo_miscentering_boost_theory.ipynb @@ -13,12 +13,21 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "\n", "%matplotlib inline" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from paper_formating import add_grid, prep_plot" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -34,13 +43,13 @@ "source": [ "import os\n", "\n", - "os.environ[\n", - " \"CLMM_MODELING_BACKEND\"\n", - "] = \"ccl\" # here you may choose ccl, nc (NumCosmo) or ct (cluster_toolkit)\n", + "os.environ[\"CLMM_MODELING_BACKEND\"] = (\n", + " \"ccl\" # here you may choose ccl, nc (NumCosmo) or ct (cluster_toolkit)\n", + ")\n", "\n", "import clmm\n", - "from clmm import Cosmology\n", - "import clmm.utils as u" + "import clmm.utils as u\n", + "from clmm import Cosmology" ] }, { @@ -90,11 +99,13 @@ "source": [ "moo = clmm.Modeling(massdef=\"mean\", delta_mdef=200, halo_profile_model=\"nfw\")\n", "\n", - "mass_cl = 1.e14\n", + "mass_cl = 1.0e14\n", "z_cl = 0.4\n", "\n", - "conc_cl = 5.4# Duffy08 value for this halo mass and redshift (see last commented cell)\n", - "halo_bias = 2.4 # Tkinker10 value for this halo mass and redshift (see last commented cell)\n", + "conc_cl = 5.4 # Duffy08 value for this halo mass and redshift (see last commented cell)\n", + "halo_bias = (\n", + " 2.4 # Tkinker10 value for this halo mass and redshift (see last commented cell)\n", + ")\n", "\n", "moo.set_cosmo(cosmo)\n", "moo.set_concentration(conc_cl)\n", @@ -115,15 +126,9 @@ "DeltaSigma_mis = moo.eval_excess_surface_density(r_proj, z_cl, r_mis=0.2)\n", "\n", "# 2halo DeltaSigma\n", - "DeltaSigma_2h = moo.eval_excess_surface_density_2h(r_proj, z_cl, halobias=6)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + "DeltaSigma_2h = moo.eval_excess_surface_density_2h(r_proj, z_cl, halobias=6)\n", + "\n", + "# boost model\n", "r_scale = 0.3\n", "nfw_boost = u.compute_nfw_boost(r_proj, r_scale, boost0=0.2)" ] @@ -141,50 +146,32 @@ "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(7,6))\n", - "\n", - "ax.loglog(r_proj, DeltaSigma,label='1-halo, centered, reference', color='k')\n", - "ax.loglog(r_proj, DeltaSigma/nfw_boost, ls=':', label='1-halo, corrected with NFW boost model', color='grey')\n", - "ax.loglog(r_proj, DeltaSigma_mis , ls ='--', label='1-halo, miscentered, R_mis = 0.2 Mpc')\n", - "ax.loglog(r_proj, DeltaSigma_2h, ls='dashdot', label='2-halo', color='green')\n", - "ax.legend(loc=1, fontsize=13)\n", - "\n", - "ax.set_xlabel('R [Mpc]', fontsize=15)\n", - "ax.set_ylabel(r'$\\Delta\\Sigma$ [M$_\\odot$ Mpc$^{-2}$]', fontsize=15)\n", - "ax.tick_params(axis='x', labelsize=15)\n", - "ax.tick_params(axis='y', labelsize=15)\n", - "ax.set_ylim([8.e10, 1e16])\n", - "\n", - "fig.tight_layout()\n", - "fig.savefig('2h_miscentering_boost.png')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Typical values for halo_bias and concentration to use in the the example, computed using CCL\n", + "fig = prep_plot(figsize=(9, 9))\n", + "ax = plt.axes()\n", "\n", - "import pyccl as ccl\n", + "ax.loglog(r_proj, DeltaSigma, label=\"1-halo (reference)\", color=\"k\")\n", + "ax.loglog(r_proj, DeltaSigma_mis, ls=\"--\", label=r\"1-halo ($R_{\\rm mis} = 0.2$Mpc)\")\n", + "ax.loglog(r_proj, DeltaSigma / nfw_boost, ls=\":\", label=\"1-halo (boost correction)\")\n", + "ax.loglog(r_proj, DeltaSigma_2h, ls=\"dashdot\", label=\"2-halo\")\n", + "ax.legend(loc=1, fontsize=6.5)\n", "\n", - "cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.83, n_s=0.96)\n", + "ax.set_xlabel(\"R [Mpc]\")\n", + "ax.set_ylabel(r\"$\\Delta\\Sigma$ [M$_\\odot$ Mpc$^{-2}$]\")\n", + "ax.set_ylim(8.0e10, 1e15)\n", + "ax.set_xlim(1e-2, 1e2)\n", "\n", - "a = 1./(1+z_cl)\n", + "add_grid(ax)\n", "\n", - "bias = ccl.halos.HaloBiasTinker10(mass_def='200m', mass_def_strict=True)\n", - "conc = ccl.halos.concentration.ConcentrationDuffy08(mass_def='200m')\n", - "\n", - "bias(cosmo, mass_cl, a), conc(cosmo, mass_cl, a)\n" + "fig.tight_layout()\n", + "fig.savefig(\"2h_miscentering_boost.png\")" ] } ], "metadata": { "kernelspec": { - "display_name": "wrk", + "display_name": "clmm", "language": "python", - "name": "wrk" + "name": "clmm" }, "language_info": { "codemirror_mode": { @@ -196,7 +183,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.9" + "version": "3.11.8" } }, "nbformat": 4, diff --git a/examples/Paper_v2.0/magnification_bias.ipynb b/examples/Paper_v2.0/magnification_bias.ipynb index f6695d3ab..9c48f7ba5 100644 --- a/examples/Paper_v2.0/magnification_bias.ipynb +++ b/examples/Paper_v2.0/magnification_bias.ipynb @@ -160,6 +160,7 @@ "ax.minorticks_on()\n", "ax.grid(which=\"minor\", lw=0.1)\n", "ax.set_axisbelow(True)\n", + "\n", "plt.tight_layout()\n", "plt.savefig(\"delta_mu.png\")" ] From bee33c256ba63bc4fe89bbb43e8fe1e6e1d8b67e Mon Sep 17 00:00:00 2001 From: Michel Date: Fri, 6 Jun 2025 11:44:35 +0200 Subject: [PATCH 09/17] add examples/Paper_v2.0/theo_diff_z_types.ipynb --- examples/Paper_v2.0/theo_diff_z_types.ipynb | 294 ++++++++++++++++++++ 1 file changed, 294 insertions(+) create mode 100644 examples/Paper_v2.0/theo_diff_z_types.ipynb diff --git a/examples/Paper_v2.0/theo_diff_z_types.ipynb b/examples/Paper_v2.0/theo_diff_z_types.ipynb new file mode 100644 index 000000000..28f9c4abd --- /dev/null +++ b/examples/Paper_v2.0/theo_diff_z_types.ipynb @@ -0,0 +1,294 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "tags": [] + }, + "source": [ + "# Model WL Profiles (different redshift inputs)\n", + "## Model profiles using different type of source redshift information as input" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we model lensing profiles by giving as input either : \n", + "- discrete source redshifts, \n", + "- a redshift distribution function,\n", + "- the value of the mean beta parameters : \n", + "$\\langle \\beta_s \\rangle = \\left\\langle \\frac{D_{LS}}{D_S}\\frac{D_\\infty}{D_{L,\\infty}}\\right\\rangle$ ,\n", + "$\\langle \\beta_s^2 \\rangle = \\left\\langle \\left(\\frac{D_{LS}}{D_S}\\frac{D_\\infty}{D_{L,\\infty}}\\right)^2 \\right\\rangle$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "## Uncomment the following line if you want to use a specific modeling backend among 'ct' (cluster-toolkit), 'ccl' (CCL) or 'nc' (Numcosmo). Default is 'ccl'\n", + "# os.environ['CLMM_MODELING_BACKEND'] = 'nc'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import clmm\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from IPython.display import Math, display\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from paper_formating import add_grid, prep_plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Make sure we know which version we're using" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "clmm.__version__" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Import mock data module and setup the configuration " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from clmm import Cosmology\n", + "from clmm.redshift.distributions import *\n", + "from clmm.support import mock_data as mock" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Mock data generation requires a defined cosmology" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mock_cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Mock data generation requires some cluster information" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cosmo = mock_cosmo\n", + "\n", + "# cluster properties from https://arxiv.org/pdf/1611.03866.pdf\n", + "cluster_id = \"SPT-CL J0000−5748\"\n", + "cluster_m = 4.56e14 # M500,c\n", + "cluster_z = 0.702\n", + "cluster_ra = 0.2499\n", + "cluster_dec = -57.8064\n", + "concentration = 5 # (arbitrary value, not from the paper)\n", + "\n", + "# source redshift distribution properties\n", + "cluster_beta_s_mean = 0.466\n", + "cluster_beta_s2_mean = 0.243\n", + "ngal_density = (\n", + " 26.0 * 100\n", + ") # density of source galaxies per arcmin^2 # (arbitrary value, not from the paper)\n", + "model_z_distrib_dict = {\"func\": desc_srd, \"name\": \"desc_srd\"}\n", + "delta_z_cut = 0.1\n", + "zsrc_min = cluster_z + delta_z_cut\n", + "zsrc_max = 3.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1 - Defining the different inputs for the source redshifts" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Discrete redshifts" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Generate the mock source catalog\n", + "\n", + "The CLMM mock data generation will provide, among other things, a redshift value for each background galaxy that is draw from the redshift distribution given by `model_z_distrib_dict`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(0)\n", + "source_catalog = mock.generate_galaxy_catalog(\n", + " cluster_m,\n", + " cluster_z,\n", + " concentration,\n", + " cosmo,\n", + " model_z_distrib_dict[\"name\"],\n", + " delta_so=500,\n", + " massdef=\"critical\",\n", + " zsrc_min=zsrc_min,\n", + " zsrc_max=zsrc_max,\n", + " ngal_density=ngal_density,\n", + " cluster_ra=cluster_ra,\n", + " cluster_dec=cluster_dec,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Beta parameters\n", + "From this udnerlying redshift distribution, one may directly compute the average $\\langle\\beta_s\\rangle$ and $\\langle\\beta_s^2\\rangle$ quantities" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "beta_label = lambda beta: rf\"\\langle\\beta_s\\rangle = {beta:.2f}\"\n", + "beta_sq_label = lambda beta_sq: rf\"\\langle\\beta_s^2\\rangle = {beta_sq:.2f}\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "It is also possible to compute $\\langle\\beta_s\\rangle$ and $\\langle\\beta_s^2\\rangle$ using galaxy shape weights:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Visualisation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z = np.linspace(0, zsrc_max, 1000)\n", + "fig = prep_plot(figsize=(9, 9))\n", + "ax = plt.axes()\n", + "\n", + "ax.hist(source_catalog[\"z\"], bins=50, alpha=0.3, density=True, label=\"discrete values\")\n", + "ax.axvline(zsrc_min, color=\"red\")\n", + "# here we multiply by a constant for visualisation purposes\n", + "ax.plot(\n", + " z,\n", + " model_z_distrib_dict[\"func\"](z) * 25,\n", + " linestyle=\"dashed\",\n", + " label=\"distribution function\",\n", + ")\n", + "\n", + "ax.text(\n", + " 1.75,\n", + " 0.9,\n", + " f\"Resulting efficiency\\n${beta_label(beta_s_mean)}$\\n${beta_sq_label(beta_s_square_mean)}$\",\n", + " bbox={\"boxstyle\": \"round\", \"alpha\": 0.5, \"edgecolor\": \".7\", \"facecolor\": \"1\"},\n", + " fontsize=8,\n", + ")\n", + "ax.text(zsrc_min, .03, \"$z_{min}$\", color=\"r\", fontsize=8, rotation=30)\n", + "\n", + "\n", + "ax.set_xlim(0, 2.99)\n", + "ax.set_xlabel(\"redshift\")\n", + "ax.set_ylabel(\"distribution\")\n", + "ax.set_yticklabels([])\n", + "ax.legend(fontsize=8)\n", + "\n", + "add_grid(ax)\n", + "plt.tight_layout()\n", + "plt.savefig(\"theo_diff_z_types.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "clmm", + "language": "python", + "name": "clmm" + }, + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 01f4f6ad8e4f4058713e6cb14c2f4ca277eb1a6c Mon Sep 17 00:00:00 2001 From: Michel Date: Sat, 7 Jun 2025 20:48:47 +0200 Subject: [PATCH 10/17] add subplots_kwargs to prep_plot --- examples/Paper_v2.0/paper_formating.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/examples/Paper_v2.0/paper_formating.py b/examples/Paper_v2.0/paper_formating.py index 2e6046dcf..b95f31ef0 100644 --- a/examples/Paper_v2.0/paper_formating.py +++ b/examples/Paper_v2.0/paper_formating.py @@ -6,6 +6,7 @@ def prep_plot( figsize=(9, 9), adjust=None, # {'left':0.2, 'right':0.95, 'top':0.94, 'bottom':0.2}, subplots=False, + subplots_kwargs=None, ): main_parameters = { "figure.dpi": 300.0, @@ -34,10 +35,12 @@ def prep_plot( _figsize = list(x / 2.54 for x in figsize) print(_figsize) if subplots: + if subplots_kwargs is None: + subplots_kwargs = {} try: - fig, axes = plt.subplots(*subplots, figsize=_figsize) + fig, axes = plt.subplots(*subplots, figsize=_figsize, **subplots_kwargs) except: - fig, axes = plt.subplots(subplots, figsize=_figsize) + fig, axes = plt.subplots(subplots, figsize=_figsize, **subplots_kwargs) else: fig = plt.figure(figsize=_figsize) if adjust is not None: From 9574d142d3f9ace030f56be7a007dfcef24fc597 Mon Sep 17 00:00:00 2001 From: Michel Date: Sat, 7 Jun 2025 20:49:38 +0200 Subject: [PATCH 11/17] add profiles to examples/Paper_v2.0/theo_diff_z_types.ipynb --- examples/Paper_v2.0/theo_diff_z_types.ipynb | 314 +++++++++++++++++++- 1 file changed, 311 insertions(+), 3 deletions(-) diff --git a/examples/Paper_v2.0/theo_diff_z_types.ipynb b/examples/Paper_v2.0/theo_diff_z_types.ipynb index 28f9c4abd..afb5554d5 100644 --- a/examples/Paper_v2.0/theo_diff_z_types.ipynb +++ b/examples/Paper_v2.0/theo_diff_z_types.ipynb @@ -146,7 +146,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 1 - Defining the different inputs for the source redshifts" + "## Different inputs for the source redshifts" ] }, { @@ -247,7 +247,7 @@ " bbox={\"boxstyle\": \"round\", \"alpha\": 0.5, \"edgecolor\": \".7\", \"facecolor\": \"1\"},\n", " fontsize=8,\n", ")\n", - "ax.text(zsrc_min, .03, \"$z_{min}$\", color=\"r\", fontsize=8, rotation=30)\n", + "ax.text(zsrc_min, 0.03, \"$z_{min}$\", color=\"r\", fontsize=8, rotation=30)\n", "\n", "\n", "ax.set_xlim(0, 2.99)\n", @@ -261,12 +261,320 @@ "plt.savefig(\"theo_diff_z_types.png\")" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Profiles" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "z_inf = 1000\n", + "\n", + "beta_s_mean = clmm.utils.compute_beta_s_mean_from_distribution(\n", + " cluster_z,\n", + " z_inf,\n", + " cosmo,\n", + " zmax=zsrc_max,\n", + " delta_z_cut=delta_z_cut,\n", + " zmin=None,\n", + " z_distrib_func=model_z_distrib_dict[\"func\"],\n", + ")\n", + "beta_s_square_mean = clmm.utils.compute_beta_s_square_mean_from_distribution(\n", + " cluster_z,\n", + " z_inf,\n", + " cosmo,\n", + " zmax=zsrc_max,\n", + " delta_z_cut=delta_z_cut,\n", + " zmin=None,\n", + " z_distrib_func=model_z_distrib_dict[\"func\"],\n", + ")\n", + "\n", + "display(Math(beta_label(beta_s_mean)))\n", + "display(Math(beta_sq_label(beta_s_square_mean)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "prof_kwargs = dict(\n", + " r_proj=np.logspace(np.log10(0.2), np.log10(5), 10),\n", + " mdelta=cluster_m,\n", + " cdelta=concentration,\n", + " z_cluster=cluster_z,\n", + " cosmo=cosmo,\n", + " delta_mdef=500,\n", + " massdef=\"critical\",\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "### g_t" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "gt = {\n", + " key: clmm.theory.compute_reduced_tangential_shear(**prof_kwargs, **_kwargs)\n", + " for key, _kwargs in (\n", + " (\n", + " \"exact\",\n", + " dict(\n", + " z_src=model_z_distrib_dict[\"func\"],\n", + " z_src_info=\"distribution\",\n", + " approx=None,\n", + " ),\n", + " ),\n", + " (\n", + " \"order1\",\n", + " dict(\n", + " z_src=[beta_s_mean, beta_s_square_mean],\n", + " z_src_info=\"beta\",\n", + " approx=\"order1\",\n", + " ),\n", + " ),\n", + " (\n", + " \"order2\",\n", + " dict(\n", + " z_src=[beta_s_mean, beta_s_square_mean],\n", + " z_src_info=\"beta\",\n", + " approx=\"order2\",\n", + " ),\n", + " ),\n", + " )\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "### mu" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "mu = {\n", + " key: clmm.theory.compute_magnification(**prof_kwargs, **_kwargs)\n", + " for key, _kwargs in (\n", + " (\n", + " \"exact\",\n", + " dict(\n", + " z_src=model_z_distrib_dict[\"func\"],\n", + " z_src_info=\"distribution\",\n", + " approx=None,\n", + " ),\n", + " ),\n", + " (\n", + " \"order1\",\n", + " dict(\n", + " z_src=[beta_s_mean, beta_s_square_mean],\n", + " z_src_info=\"beta\",\n", + " approx=\"order1\",\n", + " ),\n", + " ),\n", + " (\n", + " \"order2\",\n", + " dict(\n", + " z_src=[beta_s_mean, beta_s_square_mean],\n", + " z_src_info=\"beta\",\n", + " approx=\"order2\",\n", + " ),\n", + " ),\n", + " )\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "### delta_mu" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "alpha = [2, -0.5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "dmu = {\n", + " key: [\n", + " clmm.theory.compute_magnification_bias(**prof_kwargs, **_kwargs, alpha=_alpha)\n", + " for _alpha in alpha\n", + " ]\n", + " for key, _kwargs in (\n", + " (\n", + " \"exact\",\n", + " dict(\n", + " z_src=model_z_distrib_dict[\"func\"],\n", + " z_src_info=\"distribution\",\n", + " approx=None,\n", + " ),\n", + " ),\n", + " (\n", + " \"order1\",\n", + " dict(\n", + " z_src=[beta_s_mean, beta_s_square_mean],\n", + " z_src_info=\"beta\",\n", + " approx=\"order1\",\n", + " ),\n", + " ),\n", + " (\n", + " \"order2\",\n", + " dict(\n", + " z_src=[beta_s_mean, beta_s_square_mean],\n", + " z_src_info=\"beta\",\n", + " approx=\"order2\",\n", + " ),\n", + " ),\n", + " )\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def plot_cases(axes, radius, profiles):\n", + " axes[0].loglog(radius, profiles[0][0], **profiles[0][-1])\n", + " for case in profiles[1:]:\n", + " axes[0].loglog(radius, case[0], **case[-1])\n", + " axes[1].plot(radius, 100 * (case[0] / profiles[0][0] - 1), **case[-1])\n", + "\n", + " axes[1].set_xlabel(\"R [Mpc]\")\n", + " for ax in axes:\n", + " add_grid(ax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.ticker import NullFormatter" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "fig, axes = prep_plot(\n", + " figsize=(18, 7),\n", + " subplots=(2, 3),\n", + " subplots_kwargs=dict(sharex=True, height_ratios=[4, 1]),\n", + ")\n", + "plot_cases(\n", + " axes[:, 0],\n", + " prof_kwargs[\"r_proj\"],\n", + " [\n", + " (gt[\"exact\"], {\"label\": \"No Approx.\", \"color\": \"0\", \"lw\": 0.8}),\n", + " (gt[\"order1\"], {\"label\": \"Order 1\", \"color\": \"C0\"}),\n", + " (gt[\"order2\"], {\"label\": \"Order 2\", \"color\": \"C1\"}),\n", + " ],\n", + ")\n", + "plot_cases(\n", + " axes[:, 1],\n", + " prof_kwargs[\"r_proj\"],\n", + " [\n", + " (mu[\"exact\"], {\"label\": \"No Approx.\", \"color\": \"0\", \"lw\": 0.8}),\n", + " (mu[\"order1\"], {\"label\": \"Order 1\", \"color\": \"C0\"}),\n", + " (mu[\"order2\"], {\"label\": \"Order 2\", \"color\": \"C1\"}),\n", + " ],\n", + ")\n", + "plot_cases(\n", + " axes[:, 2],\n", + " prof_kwargs[\"r_proj\"],\n", + " [\n", + " (dmu[\"exact\"][0], {\"label\": \"No Approx.\", \"color\": \"0\", \"lw\": 0.8}),\n", + " (dmu[\"order1\"][0], {\"label\": \"Order 1\", \"color\": \"C0\"}),\n", + " (dmu[\"order2\"][0], {\"label\": \"Order 2\", \"color\": \"C1\"}),\n", + " ],\n", + ")\n", + "plot_cases(\n", + " axes[:, 2],\n", + " prof_kwargs[\"r_proj\"],\n", + " [\n", + " (dmu[\"exact\"][1], {\"label\": \"No Approx.\", \"color\": \"0\", \"ls\": \"--\", \"lw\": 0.6}),\n", + " (dmu[\"order1\"][1], {\"label\": \"Order 1\", \"color\": \"C0\", \"ls\": \"--\", \"lw\":0.8}),\n", + " (dmu[\"order2\"][1], {\"label\": \"Order 2\", \"color\": \"C1\", \"ls\": \"--\", \"lw\":0.8}),\n", + " ],\n", + ")\n", + "axes[0, 0].legend(fontsize=8)\n", + "axes[1, 0].set_ylabel(\"%\")\n", + "\n", + "for ax, label in zip(axes[0], (\"$\\gamma_t$\", \"$\\mu_t$\", \"$\\delta_{\\mu_t}$\")):\n", + " ax.set_title(label, fontsize=10)\n", + "\n", + "for ax in axes.flatten():\n", + " ax.tick_params(axis=\"both\", which=\"major\", labelsize=8)\n", + " #ax.tick_params(axis=\"both\", which=\"minor\", labelsize=8)\n", + " ax.yaxis.set_minor_formatter(NullFormatter())\n", + " ax.set_xlim(.2, 5)\n", + "\n", + "handles, labels = axes[0, -1].get_legend_handles_labels()\n", + "axes[0, -1].legend(\n", + " handles[::3],\n", + " [fr\"$\\alpha={_alpha}$\" for _alpha in alpha],\n", + " #loc=(0.38, 0.8),\n", + " fontsize=8,\n", + ")\n", + "\n", + "\n", + "plt.tight_layout()\n", + "plt.subplots_adjust(hspace=0)\n", + "\n", + "plt.savefig(\"theo_diff_z_types_profiles.png\")" + ] } ], "metadata": { From fa464aee93e38369a57dd952556e1ddd2d5debfd Mon Sep 17 00:00:00 2001 From: Shenming Fu Date: Mon, 18 Aug 2025 06:51:09 -0700 Subject: [PATCH 12/17] add nb for hsc fig --- ...le4_Fit_Halo_mass_to_HSC_data_paper2.ipynb | 828 ++++++++++++++++++ .../Example4_Fit_Halo_mass_to_HSC_data.ipynb | 99 ++- 2 files changed, 896 insertions(+), 31 deletions(-) create mode 100644 examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb diff --git a/examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb b/examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb new file mode 100644 index 000000000..859094468 --- /dev/null +++ b/examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb @@ -0,0 +1,828 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "74785dcc", + "metadata": { + "tags": [] + }, + "source": [ + "# Example 4: Using a real dataset (HSC)\n", + "## Fit halo mass to shear profile using HSC data\n", + "\n", + "_the LSST-DESC CLMM team_\n", + "\n", + "This notebook can be run on NERSC.\n", + "\n", + "Here we demonstrate how to run CLMM on real observational datasets. As an example, we use the data from the Hyper Suprime-Cam Subaru Strategic Program (HSC SSP) public releases (Aihara+2018ab, 2019; Mandelbaum+2018ab) (Credit: NAOJ / HSC Collaboration), which have similar observation conditions and data formats to the Rubin LSST.\n", + "\n", + "The steps in this notebook includes:\n", + "- [Setting things up](#Setup)\n", + "- [Selecting a cluster](#Selecting_a_cluster)\n", + "- [Downloading the published catalog at the cluster field](#Downloading_the_catalog)\n", + "- [Loading the catalog into CLMM](#Loading_the_catalog)\n", + "- [Running CLMM on the dataset](#Running_CLMM)\n", + "\n", + "Links:\n", + "\n", + "The data access of the HSC SSP Public Data Release: \n", + "https://hsc-release.mtk.nao.ac.jp/doc/index.php/data-access__pdr3/\n", + "\n", + "Shape catalog: \n", + "https://hsc-release.mtk.nao.ac.jp/doc/index.php/s16a-shape-catalog-pdr2/\n", + "\n", + "FAQ: \n", + "https://hsc-release.mtk.nao.ac.jp/doc/index.php/faq__pdr3/\n", + "\n", + "Photometric redshifts:\n", + "https://hsc-release.mtk.nao.ac.jp/doc/index.php/photometric-redshifts/\n", + "\n", + "Cluster catalog:\n", + "https://hsc-release.mtk.nao.ac.jp/doc/index.php/camira_pdr2/" + ] + }, + { + "cell_type": "markdown", + "id": "8b2a5bbe", + "metadata": {}, + "source": [ + "\n", + "## 1. Setup\n", + " \n", + "We import packages." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "89a38de9", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# %matplotlib inline\n", + "from astropy.table import Table\n", + "import pickle as pkl\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "markdown", + "id": "cf169f3c", + "metadata": {}, + "source": [ + "\n", + "## 2. Selecting a cluster\n", + "\n", + "We use the HSC SSP publications (https://hsc.mtk.nao.ac.jp/ssp/publications/) to select a list of reported massive galaxy clusters that have been measured by weak lensing. In the table below, the coordinates are for lensing peaks unless otherwise specified, and we assume h=0.7.\n", + "\n", + "Name | $$z_{cl}$$ | RA (deg) | DEC (deg) | WL Mass | Reference | Note\n", + "- | - | - | - | - | - | -\n", + "HWL16a-094 | 0.592 | 223.0801 | 0.1689 | 15.3, 7.8 | [Hamana+2020](https://ui.adsabs.harvard.edu/abs/2020PASJ...72...78H/abstract) | CAMIRA ID 1417; Miyazaki+2018 rank 34 \n", + "HWL16a-026 | 0.424 | 130.5895 | 1.6473 | 8.7, 4.7 | [Hamana+2020](https://ui.adsabs.harvard.edu/abs/2020PASJ...72...78H/abstract) | --\n", + "HWL16a-034 | 0.315 | 139.0387 | −0.3966 | 8.1, 5.6 | [Hamana+2020](https://ui.adsabs.harvard.edu/abs/2020PASJ...72...78H/abstract) | Abell 776; MACS J0916.1−0023; Miyazaki+2018 rank 8; see also Medezinski+2018 \n", + "Rank 9 | 0.312 | 37.3951 | −3.6099 | --, 5.9 | [Miyazaki+2018](https://ui.adsabs.harvard.edu/abs/2018PASJ...70S..27M/abstract) | --\n", + "Rank 48 | 0.529 | 220.7900 | 1.0509 | --, 10.4 | [Miyazaki+2018](https://ui.adsabs.harvard.edu/abs/2018PASJ...70S..27M/abstract) | --\n", + "Rank 62 | 0.592 | 216.6510 | 0.7982 | --, 10.2 | [Miyazaki+2018](https://ui.adsabs.harvard.edu/abs/2018PASJ...70S..27M/abstract) | --\n", + "MaxBCG J140.53188+03.76632 | 0.2701 | 140.54565 | 3.77820 | 44.3, 25.1 | [Medezinski+2018](https://ui.adsabs.harvard.edu/abs/2018PASJ...70S..28M%2F/abstract) | BCG center (close to the X-ray center); PSZ2 G228.50+34.95; double BCGs\n", + "XLSSC006 | 0.429 | 35.439 | −3.772 | 9.6, 5.6 | [Umetsu+2020](https://ui.adsabs.harvard.edu/abs/2020ApJ...890..148U/abstract) | X-ray center\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "71f2b585", + "metadata": { + "jp-MarkdownHeadingCollapsed": true, + "tags": [] + }, + "source": [ + "\n", + "## 3. Downloading the catalog at the cluster field\n", + "\n", + "The 3 most massive cluster-candidates are MaxBCG.J140.53188+03.76632 (in the GAMA09H field), Miyazaki+2018 (M18 hearafter) rank 48 and 62 (in the GAMA15H field). We consider MaxBCG.J140.53188+03.76632 first.\n", + "The webpage for HSC SSP data access is here [link](https://hsc-release.mtk.nao.ac.jp/doc/index.php/data-access__pdr3/).\n", + "To download the catalogs, we need to first register for a user account ([link](https://hsc-release.mtk.nao.ac.jp/datasearch/new_user/new)).\n", + "Then we log into the system, query and download the catalogs at [CAS Search](https://hsc-release.mtk.nao.ac.jp/datasearch/helps/sql_search); we use `object_id` to cross match the shape catalog, photo-z catalog, and photometry catalog. \n", + "Since the clusters are at redshift about 0.4, a radius of 10 arcmin would be about 3 Mpc. However, we make a query for the whole field to save time.\n", + "The final catalog includes shape info, photo-z, and photometry. \n", + "Here is an example of the query SQL command (thank Calum Murray; [example command](https://hsc-release.mtk.nao.ac.jp/doc/index.php/s16a-shape-catalog-pdr2/); [schema](https://hsc-release.mtk.nao.ac.jp/schema/)); the query could take 1 hour and the size of the catalog could be 400 MB (.csv.gz). If you would like to test it, please copy from \"select\" to \"--LIMIT 5\". Also select \"PDR1\" or press \"Guess release from your SQL\" at the [CAS Search](https://hsc-release.mtk.nao.ac.jp/datasearch/helps/sql_search) webpage.\n", + "To unpress the file \".gz\", use \"gunzip\" or \"gzip -d\".\n", + "\n", + "```\n", + "select\n", + " b.*, \n", + " c.ira, c.idec, \n", + " a.ishape_hsm_regauss_e1, a.ishape_hsm_regauss_e2, \n", + " a.ishape_hsm_regauss_resolution, a.ishape_hsm_regauss_sigma, \n", + " d1.photoz_best as ephor_ab_photoz_best, d1.photoz_risk_best as ephor_ab_photoz_risk_best, \n", + " d2.photoz_best as frankenz_photoz_best, d2.photoz_risk_best as frankenz_photoz_risk_best, \n", + " d3.photoz_best as nnpz_photoz_best, d3.photoz_risk_best as nnpz_photoz_risk_best, \n", + " e.icmodel_mag, e.icmodel_mag_err, \n", + " e.detect_is_primary, \n", + " e.iclassification_extendedness, \n", + " e.icmodel_flux_flags, \n", + " e.icmodel_flux, e.icmodel_flux_err, \n", + " c.iblendedness_abs_flux\n", + "from\n", + " s16a_wide.meas2 a\n", + " inner join s16a_wide.weaklensing_hsm_regauss b using (object_id)\n", + " inner join s16a_wide.meas c using (object_id)\n", + " -- inner join s16a_wide.photoz_demp d using (object_id)\n", + " -- inner join s16a_wide.photoz_ephor d using (object_id)\n", + " inner join s16a_wide.photoz_ephor_ab d1 using (object_id)\n", + " inner join s16a_wide.photoz_frankenz d2 using (object_id)\n", + " -- inner join s16a_wide.photoz_mizuki d using (object_id)\n", + " -- inner join s16a_wide.photoz_mlz d using (object_id)\n", + " inner join s16a_wide.photoz_nnpz d3 using (object_id)\n", + " inner join s16a_wide.forced e using (object_id)\n", + "-- Uncomment the specific lines depending upon the field to be used\n", + " -- where s16a_wide.search_xmm(c.skymap_id)\n", + " -- where s16a_wide.search_wide01h(c.skymap_id)\n", + " -- where s16a_wide.search_vvds(c.skymap_id)\n", + " -- where s16a_wide.search_hectomap(c.skymap_id)\n", + " -- where s16a_wide.search_gama15h(c.skymap_id)\n", + " where s16a_wide.search_gama09h(c.skymap_id)\n", + " --AND e.detect_is_primary\n", + " --AND conesearch(c.icoord, 140.54565, 3.77820, 600) \n", + " --AND NOT e.icmodel_flux_flags\n", + " --AND e.iclassification_extendedness>0.5\n", + " --LIMIT 5\n", + "```\n" + ] + }, + { + "cell_type": "markdown", + "id": "eba97e69", + "metadata": {}, + "source": [ + "\n", + "## 4. Loading the catalog into CLMM\n", + "\n", + "Once we have the catalog, we read in the catalog, make cuts on the catalog, and adjust column names to prepare for the analysis in CLMM.\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fc38b936", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "# Assume the downloaded catalog is at this path:\n", + "filename = \"197376_GAMMA09H.csv\"\n", + "catalog = filename.replace(\".csv\", \".pkl\")\n", + "if not Path(catalog).is_file():\n", + " data_0 = Table.read(filename, format=\"ascii.csv\")\n", + " pkl.dump(data_0, open(catalog, \"wb\"))\n", + "else:\n", + " data_0 = pkl.load(open(catalog, \"rb\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "03e1ca9a", + "metadata": {}, + "outputs": [], + "source": [ + "print(data_0.colnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d8da2cdb", + "metadata": {}, + "outputs": [], + "source": [ + "# We select \"frankenz\" for the test, but there are other methods available.\n", + "photoz_type = \"frankenz\"\n", + "# photoz_type = \"nnpz\"\n", + "# photoz_type = \"ephor_ab\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb7d297a", + "metadata": {}, + "outputs": [], + "source": [ + "# Cuts\n", + "def make_cuts(catalog_in):\n", + " # We consider some cuts in Mandelbaum et al. 2018 (HSC SSP Y1 shear catalog).\n", + " select = catalog_in[\"detect_is_primary\"] == \"True\"\n", + " select &= catalog_in[\"icmodel_flux_flags\"] == \"False\"\n", + " select &= catalog_in[\"iclassification_extendedness\"] > 0.5\n", + " select &= catalog_in[\"icmodel_mag_err\"] <= 2.5 / np.log(10.0) / 10.0\n", + " select &= (\n", + " catalog_in[\"ishape_hsm_regauss_e1\"] ** 2 + catalog_in[\"ishape_hsm_regauss_e2\"] ** 2 < 4.0\n", + " )\n", + " select &= catalog_in[\"icmodel_mag\"] <= 24.5\n", + " select &= catalog_in[\"iblendedness_abs_flux\"] < (10 ** (-0.375))\n", + " select &= catalog_in[\"ishape_hsm_regauss_resolution\"] >= 0.3 # similar to extendedness\n", + " select &= catalog_in[\"ishape_hsm_regauss_sigma\"] <= 0.4\n", + " # Note \"zbest\" minimizes the risk of the photo-z point estimate being far away from the true value.\n", + " # Details: https://hsc-release.mtk.nao.ac.jp/doc/wp-content/uploads/2017/02/pdr1_photoz_release_note.pdf\n", + " select &= catalog_in[\"%s_photoz_risk_best\" % photoz_type] < 0.5\n", + "\n", + " catalog_out = catalog_in[select]\n", + "\n", + " return catalog_out\n", + "\n", + "\n", + "data_1 = make_cuts(data_0)\n", + "print(len(data_0), len(data_1), len(data_1) * 1.0 / len(data_0))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bac2fc15", + "metadata": {}, + "outputs": [], + "source": [ + "# Reference: Mandelbaum et al. 2018 \"The first-year shear catalog of the Subaru Hyper Suprime-Cam Subaru Strategic Program Survey\".\n", + "# Section A.3.2: \"per-object galaxy shear estimate\".\n", + "def apply_shear_calibration(catalog_in):\n", + " e1_0 = catalog_in[\"ishape_hsm_regauss_e1\"]\n", + " e2_0 = -catalog_in[\"ishape_hsm_regauss_e2\"]\n", + " e_rms = catalog_in[\"ishape_hsm_regauss_derived_rms_e\"]\n", + " m = catalog_in[\"ishape_hsm_regauss_derived_shear_bias_m\"]\n", + " c1 = catalog_in[\"ishape_hsm_regauss_derived_shear_bias_c1\"]\n", + " c2 = catalog_in[\"ishape_hsm_regauss_derived_shear_bias_c2\"]\n", + " # Note: in the mass fit we have not implemented the weight yet.\n", + " weight = catalog_in[\"ishape_hsm_regauss_derived_shape_weight\"]\n", + "\n", + " R = 1.0 - np.sum(weight * e_rms**2.0) / np.sum(weight)\n", + " m_mean = np.sum(weight * m) / np.sum(weight)\n", + " c1_mean = np.sum(weight * c1) / np.sum(weight)\n", + " c2_mean = np.sum(weight * c2) / np.sum(weight)\n", + " print(\"R, m_mean, c1_mean, c2_mean: \", R, m_mean, c1_mean, c2_mean)\n", + "\n", + " g1 = (e1_0 / (2.0 * R) - c1) / (1.0 + m_mean)\n", + " g2 = (e2_0 / (2.0 * R) - c2) / (1.0 + m_mean)\n", + "\n", + " return g1, g2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cff6ac8e", + "metadata": {}, + "outputs": [], + "source": [ + "# Adjust column names.\n", + "def adjust_column_names(catalog_in):\n", + " # We consider a map between new and old column names.\n", + " # Note we have considered shear calibration here.\n", + " column_name_map = {\n", + " \"ra\": \"ira\",\n", + " \"dec\": \"idec\",\n", + " \"z\": \"%s_photoz_best\" % photoz_type,\n", + " \"id\": \"# object_id\",\n", + " }\n", + "\n", + " catalog_out = Table()\n", + " for i in column_name_map:\n", + " catalog_out[i] = catalog_in[column_name_map[i]]\n", + "\n", + " g1, g2 = apply_shear_calibration(catalog_in)\n", + " # CLMM uses \"epsilon shape\" rather than \"chi shape\".\n", + " catalog_out[\"e1\"] = g1\n", + " catalog_out[\"e2\"] = g2\n", + "\n", + " return catalog_out\n", + "\n", + "\n", + "data_2 = adjust_column_names(data_1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5315cfc", + "metadata": {}, + "outputs": [], + "source": [ + "# Make some figures for visualization.\n", + "def make_plots(catalog_in):\n", + " # Scatter plot\n", + " plt.figure()\n", + " plt.scatter(catalog_in[\"ra\"], catalog_in[\"dec\"], c=catalog_in[\"z\"], s=1.0, alpha=0.2)\n", + " plt.colorbar()\n", + " plt.xlabel(\"ra\")\n", + " plt.ylabel(\"dec\")\n", + " plt.title(\"z\")\n", + "\n", + " # Histogram\n", + " plt.figure()\n", + " plt.hist(catalog_in[\"z\"], bins=20)\n", + " plt.xlabel(\"z\")\n", + " plt.ylabel(\"count\")\n", + "\n", + " # Relation\n", + " plt.figure()\n", + " plt.plot(catalog_in[\"e1\"], catalog_in[\"e2\"], \",\")\n", + " plt.xlabel(\"e1\")\n", + " plt.ylabel(\"e2\")\n", + "\n", + "\n", + "# make_plots(data_2)" + ] + }, + { + "cell_type": "markdown", + "id": "5dc4163b", + "metadata": {}, + "source": [ + "\n", + "## 5. Running CLMM on the dataset\n", + "We use the functions similar to `examples/Paper_v1.0/gt_and_use_case.ipynb`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e469558f", + "metadata": {}, + "outputs": [], + "source": [ + "from clmm import Cosmology\n", + "\n", + "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", + "\n", + "# We consider MaxBCG J140.53188+03.76632\n", + "#cluster_z = 0.2701 # Cluster redshift\n", + "#cluster_ra = 140.54565\n", + "#cluster_dec = 3.77820\n", + "\n", + "#HWL16a-026\n", + "#cluster_z = 0.424\n", + "#cluster_ra = 130.5895\n", + "#cluster_dec = 1.6473\n", + "\n", + "#HWL16a-034\n", + "cluster_z = 0.315\n", + "cluster_ra = 139.0387\n", + "cluster_dec = -0.3966\n", + "\n", + "obs_galaxies = data_2\n", + "\n", + "obs_galaxies = obs_galaxies[(obs_galaxies[\"z\"] > (cluster_z + 0.1))]\n", + "\n", + "# Area cut: the query is made for the whole field and this can simplify the processing.\n", + "select = obs_galaxies[\"ra\"] < cluster_ra + 12.0 / 60.0 / np.cos(cluster_dec / 180.0 * np.pi)\n", + "select &= obs_galaxies[\"ra\"] > cluster_ra - 12.0 / 60.0 / np.cos(cluster_dec / 180.0 * np.pi)\n", + "select &= obs_galaxies[\"dec\"] < cluster_dec + 12.0 / 60.0\n", + "select &= obs_galaxies[\"dec\"] > cluster_dec - 12.0 / 60.0\n", + "obs_galaxies = obs_galaxies[select]\n", + "\n", + "obs_galaxies[\"id\"] = np.arange(len(obs_galaxies))\n", + "\n", + "# Put galaxy values on arrays.\n", + "gal_ra = obs_galaxies[\"ra\"] # Galaxies Ra in deg\n", + "gal_dec = obs_galaxies[\"dec\"] # Galaxies Dec in deg\n", + "gal_e1 = obs_galaxies[\"e1\"] # Galaxies elipticipy 1\n", + "gal_e2 = obs_galaxies[\"e2\"] # Galaxies elipticipy 2\n", + "gal_z = obs_galaxies[\"z\"] # Galaxies observed redshift\n", + "gal_id = obs_galaxies[\"id\"] # Galaxies ID" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a12e56be", + "metadata": {}, + "outputs": [], + "source": [ + "# Using the GalaxyCluster object.\n", + "\n", + "import clmm\n", + "import clmm.dataops as da\n", + "from clmm.utils import convert_units\n", + "\n", + "# Create a GCData with the galaxies.\n", + "galaxies = clmm.GCData(\n", + " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id], names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"]\n", + ")\n", + "\n", + "# Create a GalaxyCluster.\n", + "cluster = clmm.GalaxyCluster(\"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies)\n", + "\n", + "# Convert elipticities into shears for the members.\n", + "cluster.compute_tangential_and_cross_components()\n", + "print(cluster.galcat.colnames)\n", + "\n", + "# Measure profile and add profile table to the cluster.\n", + "seps = convert_units(cluster.galcat[\"theta\"], \"radians\", \"Mpc\", cluster.z, cosmo)\n", + "\n", + "cluster.make_radial_profile(\n", + " bins=da.make_bins(0.3, 3.0, 15, method=\"evenlog10width\"),\n", + " bin_units=\"Mpc\",\n", + " cosmo=cosmo,\n", + " include_empty_bins=False,\n", + " gal_ids_in_bins=True,\n", + ")\n", + "print(cluster.profile.colnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "89050c01", + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(6, 4))\n", + "ax = fig.add_axes((0, 0, 1, 1))\n", + "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", + "ax.errorbar(\n", + " cluster.profile[\"radius\"],\n", + " cluster.profile[\"gt\"],\n", + " cluster.profile[\"gt_err\"],\n", + " c=\"k\",\n", + " **errorbar_kwargs\n", + ")\n", + "ax.set_xlabel(\"r [Mpc]\", fontsize=10)\n", + "ax.set_ylabel(r\"$g_t$\", fontsize=10)\n", + "ax.grid(lw=0.3)\n", + "ax.minorticks_on()\n", + "ax.grid(which=\"minor\", lw=0.1)\n", + "ax.set_xscale('log')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a806b0f1", + "metadata": {}, + "outputs": [], + "source": [ + "# Theoretical predictions\n", + "\n", + "# Model relying on the overall redshift distribution of the sources (WtG III Applegate et al. 2014).\n", + "# Note the concentration of MaxBCG J140.53188+03.76632 was not reported.\n", + "# The value from the stacked sample in the paper is ~7.\n", + "# For the mass scale, a typical c-M relation (e.g. Child et al. 2018) would give c~3 though.\n", + "# And we have not considered a c-M relation in the fitting.\n", + "z_inf = 1000\n", + "concentration = 4.0\n", + "\n", + "bs_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo))\n", + "bs2_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo) ** 2)\n", + "\n", + "\n", + "def predict_reduced_tangential_shear_redshift_distribution(profile, logm):\n", + " gt = clmm.compute_reduced_tangential_shear(\n", + " r_proj=profile[\"radius\"], # Radial component of the profile\n", + " mdelta=10**logm, # Mass of the cluster [M_sun]\n", + " cdelta=concentration, # Concentration of the cluster\n", + " z_cluster=cluster_z, # Redshift of the cluster\n", + " z_src=(bs_mean, bs2_mean), # tuple of (bs_mean, bs2_mean)\n", + " z_src_info=\"beta\",\n", + " approx=\"order1\",\n", + " cosmo=cosmo,\n", + " delta_mdef=200,\n", + " massdef=\"critical\",\n", + " halo_profile_model=\"nfw\",\n", + " )\n", + " return gt\n", + "\n", + "\n", + "# Model using individual redshift and radial information, to compute the averaged shear in each radial bin, based on the galaxies actually present in that bin.\n", + "cluster.galcat[\"theta_mpc\"] = convert_units(\n", + " cluster.galcat[\"theta\"], \"radians\", \"mpc\", cluster.z, cosmo\n", + ")\n", + "\n", + "\n", + "def predict_reduced_tangential_shear_individual_redshift(profile, logm):\n", + " return np.array(\n", + " [\n", + " np.mean(\n", + " clmm.compute_reduced_tangential_shear(\n", + " # Radial component of each source galaxy inside the radial bin\n", + " r_proj=cluster.galcat[radial_bin[\"gal_id\"]][\"theta_mpc\"],\n", + " mdelta=10**logm, # Mass of the cluster [M_sun]\n", + " cdelta=concentration, # Concentration of the cluster\n", + " z_cluster=cluster_z, # Redshift of the cluster\n", + " # Redshift value of each source galaxy inside the radial bin\n", + " z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n", + " cosmo=cosmo,\n", + " delta_mdef=200,\n", + " massdef=\"critical\",\n", + " halo_profile_model=\"nfw\",\n", + " )\n", + " )\n", + " for radial_bin in profile\n", + " ]\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48dd5752", + "metadata": {}, + "outputs": [], + "source": [ + "# Mass fitting\n", + "\n", + "mask_for_fit = cluster.profile[\"n_src\"] > 2\n", + "data_for_fit = cluster.profile[mask_for_fit]\n", + "\n", + "from clmm.support.sampler import fitters\n", + "\n", + "\n", + "def fit_mass(predict_function):\n", + " popt, pcov = fitters[\"curve_fit\"](\n", + " predict_function,\n", + " data_for_fit,\n", + " data_for_fit[\"gt\"],\n", + " data_for_fit[\"gt_err\"],\n", + " bounds=[10.0, 17.0],\n", + " )\n", + " logm, logm_err = popt[0], np.sqrt(pcov[0][0])\n", + " return {\n", + " \"logm\": logm,\n", + " \"logm_err\": logm_err,\n", + " \"m\": 10**logm,\n", + " \"m_err\": (10**logm) * logm_err * np.log(10),\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e6a4e04", + "metadata": {}, + "outputs": [], + "source": [ + "# In the paper, the measured mass is 44.3 {+ 30.3} {- 19.9} * 10^14 Msun (M200c,WL).\n", + "# For convenience, we consider a mean value for the errorbar.\n", + "# We build a dictionary based on that result.\n", + "m_paper = 44.3e14\n", + "m_err_paper = 25.1e14\n", + "logm_paper = np.log10(m_paper)\n", + "logm_err_paper = m_err_paper / (10**logm_paper) / np.log(10)\n", + "paper_value = {\n", + " \"logm\": logm_paper,\n", + " \"logm_err\": logm_err_paper,\n", + " \"m\": 10**logm_paper,\n", + " \"m_err\": (10**logm_paper) * logm_err_paper * np.log(10),\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc3a7934", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "fit_redshift_distribution = fit_mass(predict_reduced_tangential_shear_redshift_distribution)\n", + "fit_individual_redshift = fit_mass(predict_reduced_tangential_shear_individual_redshift)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "79e35734", + "metadata": {}, + "outputs": [], + "source": [ + "print(\n", + " f'Best fit mass for N(z) model = {fit_redshift_distribution[\"m\"]:.3e} +/- {fit_redshift_distribution[\"m_err\"]:.3e} Msun'\n", + ")\n", + "print(\n", + " f'Best fit mass for individual redshift and radius = {fit_individual_redshift[\"m\"]:.3e} +/- {fit_individual_redshift[\"m_err\"]:.3e} Msun'\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73202bf8", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualization of the results.\n", + "def get_predicted_shear(predict_function, fit_values):\n", + " gt_est = predict_function(data_for_fit, fit_values[\"logm\"])\n", + " gt_est_err = [\n", + " predict_function(data_for_fit, fit_values[\"logm\"] + i * fit_values[\"logm_err\"])\n", + " for i in (-3, 3)\n", + " ]\n", + " return gt_est, gt_est_err\n", + "\n", + "\n", + "gt_redshift_distribution, gt_err_redshift_distribution = get_predicted_shear(\n", + " predict_reduced_tangential_shear_redshift_distribution, fit_redshift_distribution\n", + ")\n", + "gt_individual_redshift, gt_err_individual_redshift = get_predicted_shear(\n", + " predict_reduced_tangential_shear_individual_redshift, fit_individual_redshift\n", + ")\n", + "\n", + "gt_paper1, gt_err_paper1 = get_predicted_shear(\n", + " predict_reduced_tangential_shear_redshift_distribution, paper_value\n", + ")\n", + "gt_paper2, gt_err_paper2 = get_predicted_shear(\n", + " predict_reduced_tangential_shear_individual_redshift, paper_value\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec85a1b1", + "metadata": {}, + "outputs": [], + "source": [ + "chi2_redshift_distribution_dof = np.sum(\n", + " (gt_redshift_distribution - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", + ") / (len(data_for_fit) - 1)\n", + "chi2_individual_redshift_dof = np.sum(\n", + " (gt_individual_redshift - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", + ") / (len(data_for_fit) - 1)\n", + "\n", + "print(f\"Reduced chi2 (N(z) model) = {chi2_redshift_distribution_dof}\")\n", + "print(f\"Reduced chi2 (individual (R,z) model) = {chi2_individual_redshift_dof}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a7ca7d4", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.ticker import MultipleLocator\n", + "\n", + "fig = plt.figure(figsize=(6, 6))\n", + "gt_ax = fig.add_axes((0.25, 0.42, 0.7, 0.55))\n", + "gt_ax.errorbar(\n", + " data_for_fit[\"radius\"], data_for_fit[\"gt\"], data_for_fit[\"gt_err\"], c=\"k\", **errorbar_kwargs\n", + ")\n", + "\n", + "# Points in grey have not been used for the fit.\n", + "gt_ax.errorbar(\n", + " cluster.profile[\"radius\"][~mask_for_fit],\n", + " cluster.profile[\"gt\"][~mask_for_fit],\n", + " cluster.profile[\"gt_err\"][~mask_for_fit],\n", + " c=\"grey\",\n", + " **errorbar_kwargs,\n", + ")\n", + "\n", + "pow10 = 15\n", + "mlabel = lambda name, fits: (\n", + " rf\"$M_{{fit}}^{{{name}}} = \"\n", + " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", + " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", + " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", + ")\n", + "\n", + "# The model for the 1st method.\n", + "gt_ax.loglog(\n", + " data_for_fit[\"radius\"],\n", + " gt_redshift_distribution,\n", + " \"-C1\",\n", + " label=mlabel(\"N(z)\", fit_redshift_distribution),\n", + " lw=0.5,\n", + ")\n", + "gt_ax.fill_between(\n", + " data_for_fit[\"radius\"], *gt_err_redshift_distribution, lw=0, color=\"C1\", alpha=0.2\n", + ")\n", + "\n", + "# The model for the 2nd method.\n", + "gt_ax.loglog(\n", + " data_for_fit[\"radius\"],\n", + " gt_individual_redshift,\n", + " \"-C2\",\n", + " label=mlabel(\"z,R\", fit_individual_redshift),\n", + " lw=0.5,\n", + ")\n", + "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_individual_redshift, lw=0, color=\"C2\", alpha=0.2)\n", + "\n", + "# The value in the reference paper.\n", + "#gt_ax.loglog(\n", + "# data_for_fit[\"radius\"], gt_paper1, \"-C3\", label=mlabel(\"paper; N(z)\", paper_value), lw=0.5\n", + "#)\n", + "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper1, lw=0, color=\"C3\", alpha=0.2)\n", + "#\n", + "#gt_ax.loglog(\n", + "# data_for_fit[\"radius\"], gt_paper2, \"-C4\", label=mlabel(\"paper; Z,R\", paper_value), lw=0.5\n", + "#)\n", + "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper2, lw=0, color=\"C4\", alpha=0.2)\n", + "\n", + "gt_ax.set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", + "gt_ax.legend(fontsize=\"x-large\")\n", + "gt_ax.set_xticklabels([])\n", + "gt_ax.tick_params(\"x\", labelsize=\"x-large\")\n", + "gt_ax.tick_params(\"y\", labelsize=\"x-large\")\n", + "\n", + "\n", + "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", + "errorbar_kwargs2[\"markersize\"] = 5\n", + "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "res_ax = fig.add_axes((0.25, 0.2, 0.7, 0.2))\n", + "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.25\n", + "\n", + "\n", + "res_ax.errorbar(\n", + " data_for_fit[\"radius\"],\n", + " data_for_fit[\"gt\"] / gt_redshift_distribution - 1,\n", + " yerr=data_for_fit[\"gt_err\"] / gt_redshift_distribution,\n", + " marker=\"s\",\n", + " c=\"C1\",\n", + " **errorbar_kwargs2,\n", + ")\n", + "#errorbar_kwargs2[\"markersize\"] = 3\n", + "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "\n", + "res_ax.errorbar(\n", + " data_for_fit[\"radius\"] * delta,\n", + " data_for_fit[\"gt\"] / gt_individual_redshift - 1,\n", + " yerr=data_for_fit[\"gt_err\"] / gt_individual_redshift,\n", + " marker=\"*\",\n", + " c=\"C2\",\n", + " **errorbar_kwargs2,\n", + ")\n", + "res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", + "\n", + "res_ax.set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", + "res_ax.set_xscale(\"log\")\n", + "\n", + "res_ax.set_ylim(-5, 5)\n", + "res_ax.yaxis.set_minor_locator(MultipleLocator(10))\n", + "\n", + "res_ax.tick_params(\"x\", labelsize=\"x-large\")\n", + "res_ax.tick_params(\"y\", labelsize=\"x-large\")\n", + "\n", + "for ax in (gt_ax, res_ax):\n", + " ax.grid(lw=0.3, ls=':')\n", + " ax.minorticks_on()\n", + " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "bc63da5c", + "metadata": {}, + "source": [ + "### Reference\n", + "\n", + "Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, Publications of the Astronomical Society of Japan, 70, S4\n", + "\n", + "Aihara, H., Armstrong, R., Bickerton, S., et al. 2018, Publications of the Astronomical Society of Japan, 70, S8\n", + "\n", + "Aihara, H., AlSayyad, Y., Ando, M., et al. 2019, Publications of the Astronomical Society of Japan, 71\n", + "\n", + "Hamana, T., Shirasaki, M., Lin, Y.-T., 2020, Publications of the Astronomical Society of Japan, 72, 78\n", + "\n", + "Mandelbaum, R., Miyatake, H., Hamana, T., et al. 2018, Publications of the Astronomical Society of Japan, 70, S25\n", + "\n", + "Mandelbaum, R., Lanusse, F., Leauthaud, A., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 3170\n", + "\n", + "Medezinski, E., Battaglia, N., Umetsu, K., et al., 2018, Publications of the Astronomical Society of Japan, 70, S28\n", + "\n", + "Miyazaki, S., Oguri, M., Hamana, T., et al., 2018, Publications of the Astronomical Society of Japan, 70, S27\n", + "\n", + "Umetsu, K., Sereno, M., Lieu, M., et al., 2020, Astrophysical Journal, 890, 148\n", + "\n", + "\n", + "\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "conda-clmmenv", + "language": "python", + "name": "conda-clmmenv" + }, + "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.9.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb b/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb index 74502f5c7..859094468 100644 --- a/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb +++ b/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb @@ -2,6 +2,7 @@ "cells": [ { "cell_type": "markdown", + "id": "74785dcc", "metadata": { "tags": [] }, @@ -42,6 +43,7 @@ }, { "cell_type": "markdown", + "id": "8b2a5bbe", "metadata": {}, "source": [ "\n", @@ -53,6 +55,7 @@ { "cell_type": "code", "execution_count": null, + "id": "89a38de9", "metadata": {}, "outputs": [], "source": [ @@ -67,6 +70,7 @@ }, { "cell_type": "markdown", + "id": "cf169f3c", "metadata": {}, "source": [ "\n", @@ -89,6 +93,7 @@ }, { "cell_type": "markdown", + "id": "71f2b585", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] @@ -150,6 +155,7 @@ }, { "cell_type": "markdown", + "id": "eba97e69", "metadata": {}, "source": [ "\n", @@ -163,6 +169,7 @@ { "cell_type": "code", "execution_count": null, + "id": "fc38b936", "metadata": {}, "outputs": [], "source": [ @@ -180,6 +187,7 @@ { "cell_type": "code", "execution_count": null, + "id": "03e1ca9a", "metadata": {}, "outputs": [], "source": [ @@ -189,6 +197,7 @@ { "cell_type": "code", "execution_count": null, + "id": "d8da2cdb", "metadata": {}, "outputs": [], "source": [ @@ -201,6 +210,7 @@ { "cell_type": "code", "execution_count": null, + "id": "bb7d297a", "metadata": {}, "outputs": [], "source": [ @@ -234,6 +244,7 @@ { "cell_type": "code", "execution_count": null, + "id": "bac2fc15", "metadata": {}, "outputs": [], "source": [ @@ -241,7 +252,7 @@ "# Section A.3.2: \"per-object galaxy shear estimate\".\n", "def apply_shear_calibration(catalog_in):\n", " e1_0 = catalog_in[\"ishape_hsm_regauss_e1\"]\n", - " e2_0 = catalog_in[\"ishape_hsm_regauss_e2\"]\n", + " e2_0 = -catalog_in[\"ishape_hsm_regauss_e2\"]\n", " e_rms = catalog_in[\"ishape_hsm_regauss_derived_rms_e\"]\n", " m = catalog_in[\"ishape_hsm_regauss_derived_shear_bias_m\"]\n", " c1 = catalog_in[\"ishape_hsm_regauss_derived_shear_bias_c1\"]\n", @@ -264,6 +275,7 @@ { "cell_type": "code", "execution_count": null, + "id": "cff6ac8e", "metadata": {}, "outputs": [], "source": [ @@ -296,6 +308,7 @@ { "cell_type": "code", "execution_count": null, + "id": "a5315cfc", "metadata": {}, "outputs": [], "source": [ @@ -327,6 +340,7 @@ }, { "cell_type": "markdown", + "id": "5dc4163b", "metadata": {}, "source": [ "\n", @@ -337,6 +351,7 @@ { "cell_type": "code", "execution_count": null, + "id": "e469558f", "metadata": {}, "outputs": [], "source": [ @@ -345,9 +360,19 @@ "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", "\n", "# We consider MaxBCG J140.53188+03.76632\n", - "cluster_z = 0.2701 # Cluster redshift\n", - "cluster_ra = 140.54565 # Cluster Ra in deg\n", - "cluster_dec = 3.77820 # Cluster Dec in deg\n", + "#cluster_z = 0.2701 # Cluster redshift\n", + "#cluster_ra = 140.54565\n", + "#cluster_dec = 3.77820\n", + "\n", + "#HWL16a-026\n", + "#cluster_z = 0.424\n", + "#cluster_ra = 130.5895\n", + "#cluster_dec = 1.6473\n", + "\n", + "#HWL16a-034\n", + "cluster_z = 0.315\n", + "cluster_ra = 139.0387\n", + "cluster_dec = -0.3966\n", "\n", "obs_galaxies = data_2\n", "\n", @@ -374,6 +399,7 @@ { "cell_type": "code", "execution_count": null, + "id": "a12e56be", "metadata": {}, "outputs": [], "source": [ @@ -411,12 +437,13 @@ { "cell_type": "code", "execution_count": null, + "id": "89050c01", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(6, 4))\n", "ax = fig.add_axes((0, 0, 1, 1))\n", - "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=1, elinewidth=0.5, capthick=0.5)\n", + "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", "ax.errorbar(\n", " cluster.profile[\"radius\"],\n", " cluster.profile[\"gt\"],\n", @@ -429,12 +456,14 @@ "ax.grid(lw=0.3)\n", "ax.minorticks_on()\n", "ax.grid(which=\"minor\", lw=0.1)\n", + "ax.set_xscale('log')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, + "id": "a806b0f1", "metadata": {}, "outputs": [], "source": [ @@ -501,6 +530,7 @@ { "cell_type": "code", "execution_count": null, + "id": "48dd5752", "metadata": {}, "outputs": [], "source": [ @@ -532,6 +562,7 @@ { "cell_type": "code", "execution_count": null, + "id": "9e6a4e04", "metadata": {}, "outputs": [], "source": [ @@ -553,6 +584,7 @@ { "cell_type": "code", "execution_count": null, + "id": "dc3a7934", "metadata": {}, "outputs": [], "source": [ @@ -564,6 +596,7 @@ { "cell_type": "code", "execution_count": null, + "id": "79e35734", "metadata": {}, "outputs": [], "source": [ @@ -578,6 +611,7 @@ { "cell_type": "code", "execution_count": null, + "id": "73202bf8", "metadata": {}, "outputs": [], "source": [ @@ -609,6 +643,7 @@ { "cell_type": "code", "execution_count": null, + "id": "ec85a1b1", "metadata": {}, "outputs": [], "source": [ @@ -626,6 +661,7 @@ { "cell_type": "code", "execution_count": null, + "id": "5a7ca7d4", "metadata": {}, "outputs": [], "source": [ @@ -677,25 +713,25 @@ "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_individual_redshift, lw=0, color=\"C2\", alpha=0.2)\n", "\n", "# The value in the reference paper.\n", - "gt_ax.loglog(\n", - " data_for_fit[\"radius\"], gt_paper1, \"-C3\", label=mlabel(\"paper; N(z)\", paper_value), lw=0.5\n", - ")\n", - "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper1, lw=0, color=\"C3\", alpha=0.2)\n", - "\n", - "gt_ax.loglog(\n", - " data_for_fit[\"radius\"], gt_paper2, \"-C4\", label=mlabel(\"paper; Z,R\", paper_value), lw=0.5\n", - ")\n", - "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper2, lw=0, color=\"C4\", alpha=0.2)\n", - "\n", - "gt_ax.set_ylabel(r\"$g_t$\", fontsize=8)\n", - "gt_ax.legend(fontsize=6)\n", + "#gt_ax.loglog(\n", + "# data_for_fit[\"radius\"], gt_paper1, \"-C3\", label=mlabel(\"paper; N(z)\", paper_value), lw=0.5\n", + "#)\n", + "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper1, lw=0, color=\"C3\", alpha=0.2)\n", + "#\n", + "#gt_ax.loglog(\n", + "# data_for_fit[\"radius\"], gt_paper2, \"-C4\", label=mlabel(\"paper; Z,R\", paper_value), lw=0.5\n", + "#)\n", + "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper2, lw=0, color=\"C4\", alpha=0.2)\n", + "\n", + "gt_ax.set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", + "gt_ax.legend(fontsize=\"x-large\")\n", "gt_ax.set_xticklabels([])\n", - "gt_ax.tick_params(\"x\", labelsize=8)\n", - "gt_ax.tick_params(\"y\", labelsize=8)\n", + "gt_ax.tick_params(\"x\", labelsize=\"x-large\")\n", + "gt_ax.tick_params(\"y\", labelsize=\"x-large\")\n", "\n", "\n", "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", - "errorbar_kwargs2[\"markersize\"] = 3\n", + "errorbar_kwargs2[\"markersize\"] = 5\n", "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "res_ax = fig.add_axes((0.25, 0.2, 0.7, 0.2))\n", "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.25\n", @@ -709,8 +745,8 @@ " c=\"C1\",\n", " **errorbar_kwargs2,\n", ")\n", - "errorbar_kwargs2[\"markersize\"] = 3\n", - "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "#errorbar_kwargs2[\"markersize\"] = 3\n", + "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "\n", "res_ax.errorbar(\n", " data_for_fit[\"radius\"] * delta,\n", @@ -720,26 +756,27 @@ " c=\"C2\",\n", " **errorbar_kwargs2,\n", ")\n", - "res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=8)\n", + "res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", "\n", - "res_ax.set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=8)\n", + "res_ax.set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", "res_ax.set_xscale(\"log\")\n", "\n", "res_ax.set_ylim(-5, 5)\n", "res_ax.yaxis.set_minor_locator(MultipleLocator(10))\n", "\n", - "res_ax.tick_params(\"x\", labelsize=8)\n", - "res_ax.tick_params(\"y\", labelsize=8)\n", + "res_ax.tick_params(\"x\", labelsize=\"x-large\")\n", + "res_ax.tick_params(\"y\", labelsize=\"x-large\")\n", "\n", "for ax in (gt_ax, res_ax):\n", - " ax.grid(lw=0.3)\n", + " ax.grid(lw=0.3, ls=':')\n", " ax.minorticks_on()\n", - " ax.grid(which=\"minor\", lw=0.1)\n", + " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", "plt.show()" ] }, { "cell_type": "markdown", + "id": "bc63da5c", "metadata": {}, "source": [ "### Reference\n", @@ -769,9 +806,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "conda-clmmenv", "language": "python", - "name": "python3" + "name": "conda-clmmenv" }, "language_info": { "codemirror_mode": { @@ -783,7 +820,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.9" + "version": "3.9.9" } }, "nbformat": 4, From b3fda3afb2c06ffc408b8a5e615319124f5889f8 Mon Sep 17 00:00:00 2001 From: Shenming Fu Date: Tue, 19 Aug 2025 01:55:21 -0700 Subject: [PATCH 13/17] add nb for des fig --- .../Example5_Fit_Halo_mass_to_DES_data.ipynb | 76 ++++++++++++------- 1 file changed, 50 insertions(+), 26 deletions(-) diff --git a/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb b/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb index 3e24d6f69..cf18a85a5 100644 --- a/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb +++ b/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb @@ -131,7 +131,9 @@ "source": [ "%%time\n", "# Assume the downloaded catalog is at this path:\n", - "filename = \"ACOS520_DES.csv\"\n", + "#filename = \"ACOS520_DES.csv\"\n", + "#filename = \"ID1.csv\"\n", + "filename = \"ID10.csv\"\n", "catalog = filename.replace(\".csv\", \".pkl\")\n", "if not Path(catalog).is_file():\n", " data_0 = Table.read(filename, format=\"ascii.csv\")\n", @@ -150,6 +152,17 @@ "print(data_0.colnames)" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5a40301-a3fd-4b29-a7c7-8bccf4c01c72", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.table import unique\n", + "data_0 = unique(data_0, [\"ra\",\"dec\"])" + ] + }, { "cell_type": "markdown", "id": "e66ea8ec-90f7-4df1-86a1-678a483da776", @@ -290,9 +303,19 @@ "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", "\n", "# We consider RMJ051637.4-543001.6 (ACO S520)\n", - "cluster_z = 0.30416065 # Cluster redshift\n", - "cluster_ra = 79.155704 # Cluster Ra in deg\n", - "cluster_dec = -54.500456 # Cluster Dec in deg\n", + "#cluster_z = 0.30416065 # Cluster redshift\n", + "#cluster_ra = 79.155704 # Cluster Ra in deg\n", + "#cluster_dec = -54.500456 # Cluster Dec in deg\n", + "\n", + "# ID1\n", + "#cluster_z = 0.4298039972782135 # Cluster redshift\n", + "#cluster_ra = 43.564574 # Cluster Ra in deg\n", + "#cluster_dec = -58.95297 # Cluster Dec in deg\n", + "\n", + "# ID10\n", + "cluster_z = 0.3264264166355133 # Cluster redshift\n", + "cluster_ra = 323.800394 # Cluster Ra in deg\n", + "cluster_dec = -1.04959 # Cluster Dec in deg\n", "\n", "# Select background galaxies\n", "obs_galaxies = obs_galaxies[(obs_galaxies[\"z\"] > (cluster_z + 0.1)) & (obs_galaxies[\"z\"] < 1.5)]\n", @@ -339,7 +362,7 @@ "\n", "# Measure profile and add profile table to the cluster.\n", "cluster.make_radial_profile(\n", - " bins=da.make_bins(0.2, 5.0, 7, method=\"evenlog10width\"),\n", + " bins=da.make_bins(0.2, 4.0, 6, method=\"evenlog10width\"),\n", " bin_units=\"Mpc\",\n", " cosmo=cosmo,\n", " include_empty_bins=False,\n", @@ -356,7 +379,7 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(7, 5), ncols=1, nrows=1)\n", - "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=1, elinewidth=0.5, capthick=0.5)\n", + "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", "ax.errorbar(\n", " cluster.profile[\"radius\"],\n", " cluster.profile[\"gt\"],\n", @@ -571,7 +594,8 @@ "fig, axes = plt.subplots(\n", " nrows=2,\n", " ncols=1,\n", - " figsize=(8, 6),\n", + " #figsize=(8, 6),\n", + " figsize=(8*0.8, 6*0.8),\n", " gridspec_kw={\"height_ratios\": [3, 1], \"wspace\": 0.4, \"hspace\": 0.03},\n", ")\n", "\n", @@ -588,11 +612,11 @@ " **errorbar_kwargs,\n", ")\n", "\n", - "pow10 = 14\n", + "pow10 = 15\n", "mlabel = lambda name, fits: (\n", " rf\"$M_{{fit}}^{{{name}}} = \"\n", - " rf'{fits[\"m\"]/10**pow10:.3f}\\pm'\n", - " rf'{fits[\"m_err\"]/10**pow10:.3f}'\n", + " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", + " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", ")\n", "\n", @@ -622,16 +646,16 @@ ")\n", "\n", "\n", - "axes[0].set_ylabel(r\"$g_t$\", fontsize=12)\n", - "axes[0].legend(fontsize=12, loc=4)\n", + "axes[0].set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", + "axes[0].legend(fontsize=\"x-large\", loc=\"lower left\")\n", "axes[0].set_xticklabels([])\n", - "axes[0].tick_params(\"x\", labelsize=12)\n", - "axes[0].tick_params(\"y\", labelsize=12)\n", + "axes[0].tick_params(\"x\", labelsize=\"x-large\")\n", + "axes[0].tick_params(\"y\", labelsize=\"x-large\")\n", "axes[1].set_ylim(1.0e-3, 0.5)\n", "\n", "\n", "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", - "errorbar_kwargs2[\"markersize\"] = 3\n", + "errorbar_kwargs2[\"markersize\"] = 5\n", "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "\n", "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.15\n", @@ -645,8 +669,8 @@ " c=\"C1\",\n", " **errorbar_kwargs2,\n", ")\n", - "errorbar_kwargs2[\"markersize\"] = 3\n", - "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "#errorbar_kwargs2[\"markersize\"] = 3\n", + "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "\n", "axes[1].errorbar(\n", " data_for_fit[\"radius\"] * delta,\n", @@ -656,20 +680,20 @@ " c=\"C2\",\n", " **errorbar_kwargs2,\n", ")\n", - "axes[1].set_xlabel(r\"$R$ [Mpc]\", fontsize=12)\n", + "axes[1].set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", "\n", - "axes[1].set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=12)\n", + "axes[1].set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", "axes[1].set_xscale(\"log\")\n", "\n", "axes[1].set_ylim(-5, 5)\n", "\n", - "axes[1].tick_params(\"x\", labelsize=12)\n", - "axes[1].tick_params(\"y\", labelsize=12)\n", + "axes[1].tick_params(\"x\", labelsize=\"x-large\")\n", + "axes[1].tick_params(\"y\", labelsize=\"x-large\")\n", "\n", "for ax in axes:\n", - " ax.grid(lw=0.3)\n", + " ax.grid(lw=0.3, ls=':')\n", " ax.minorticks_on()\n", - " ax.grid(which=\"minor\", lw=0.1)\n", + " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", "plt.show()\n", "\n", "# Note since we made cuts on the catalog, the redshift distribution of the remaining sources might not be representative." @@ -693,9 +717,9 @@ ], "metadata": { "kernelspec": { - "display_name": "wrk", + "display_name": "conda-clmmenv", "language": "python", - "name": "wrk" + "name": "conda-clmmenv" }, "language_info": { "codemirror_mode": { @@ -707,7 +731,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.6" + "version": "3.9.9" } }, "nbformat": 4, From 1c265183fcdda3606a0b20f6a3f19d6d75bea480 Mon Sep 17 00:00:00 2001 From: Shenming Fu Date: Tue, 19 Aug 2025 02:00:50 -0700 Subject: [PATCH 14/17] add nb for des fig --- ...le5_Fit_Halo_mass_to_DES_data_paper2.ipynb | 739 ++++++++++++++++++ 1 file changed, 739 insertions(+) create mode 100644 examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb diff --git a/examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb b/examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb new file mode 100644 index 000000000..cf18a85a5 --- /dev/null +++ b/examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb @@ -0,0 +1,739 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e320d2bf-dce6-4a5b-ab38-04b0c948155f", + "metadata": { + "tags": [] + }, + "source": [ + "# Example 5: Using a real datasets (DES)\n", + "## Fit halo mass to shear profile using DES data\n", + "\n", + "_the LSST-DESC CLMM team_\n", + "\n", + "This notebook can be run on NERSC.\n", + "\n", + "Here we demonstrate how to run CLMM on real observational datasets. As an example, we use the data from the Dark Energy Survey (DES) public releases. The catalogs can be accessed from the NOIRLab Astro Data Lab.\n", + "\n", + "The steps in this notebook includes:\n", + "- [Setting things up](#Setup)\n", + "- [Selecting a cluster](#Selecting_a_cluster)\n", + "- [Downloading the published catalog at the cluster field](#Downloading_the_catalog)\n", + "- [Loading the catalog into CLMM](#Loading_the_catalog)\n", + "- [Running CLMM on the dataset](#Running_CLMM)\n", + "\n", + "Acknowledgement\n", + "\n", + "DES data: https://des.ncsa.illinois.edu/thanks\n", + "\n", + "Astro Data Lab: https://datalab.noirlab.edu/acknowledgements.php\n" + ] + }, + { + "cell_type": "markdown", + "id": "4ef8f003-be65-44a3-8b0a-04749d1f7f13", + "metadata": {}, + "source": [ + "\n", + "## 1. Setup\n", + " \n", + "We import packages." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d3503a1a-103f-48aa-b600-ef2d72de82a0", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from astropy.table import Table\n", + "import pickle as pkl\n", + "from pathlib import Path" + ] + }, + { + "cell_type": "markdown", + "id": "5e228cf4-f5f3-4b55-97b7-f1e022a5b29c", + "metadata": {}, + "source": [ + "\n", + "## 2. Selecting a cluster\n", + "\n", + "We use the DES Y1 redMaPPer Catalogs (https://des.ncsa.illinois.edu/releases/y1a1/key-catalogs/key-redmapper) to select a list of high-richness (LAMBDA) galaxy clusters, which likely have high masses.\n", + "\n", + "Name | RA (deg) | DEC (deg) | Z_LAMBDA | LAMBDA | \n", + "- | - | - | - | - \n", + "RMJ025415.5 585710.7 | 43.564574 | -58.95297 | 0.429804 | 234.50368 \n", + "RMJ051637.4 543001.6 | 79.155704 | -54.500456 | 0.30416065 | 195.06956 \n", + "RMJ224851.8 443106.3 | 342.215897 | -44.518403 | 0.3514858 | 178.83827 \n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "c684515e-e176-4b92-9509-0217ade681a0", + "metadata": { + "jp-MarkdownHeadingCollapsed": true, + "tags": [] + }, + "source": [ + "\n", + "## 3. Downloading the catalog at the cluster field\n", + "\n", + "We consider RMJ051637.4-543001.6 (ACO S520) as an example.\n", + "We can access the DES catalog from NOIRLab Data Lab (https://datalab.noirlab.edu/query.php?name=des_dr1.shape_metacal_riz_unblind). No registration is required.\n", + "We make the query and download the catalogs in \"Query Interface\". \n", + "We use `coadd_objects_id` to cross match the shape catalog and photo-z catalog (https://datalab.noirlab.edu/query.php?name=des_dr1.photo_z). \n", + "Since the cluster is at redshift about 0.3, a radius of 0.3 deg would be about a radial distance of 5 Mpc. \n", + "The final catalog includes shape info and photo-z. \n", + "Here is an example of the query SQL command. \n", + "The query could take a few minutes and the size of the catalog is about 1.4 MB (.csv). \n", + "\n", + "```\n", + "SELECT P.mean_z, \n", + "C.ra, C.dec, C.e1, C.e2, C.r11, C.r12, C.r21, C.r22 \n", + "FROM des_dr1.photo_z as P\n", + "INNER JOIN des_dr1.shape_metacal_riz_unblind as C\n", + "ON P.coadd_objects_id=C.coadd_objects_id\n", + "WHERE 't' = Q3C_RADIAL_QUERY(C.ra, C.dec,79.155704, -54.500456, 0.3) \n", + "AND P.minchi2<1\n", + "AND P.z_sigma<0.1\n", + "AND C.flags_select=0\n", + "\n", + "```\n" + ] + }, + { + "cell_type": "markdown", + "id": "9f4535f6-23b6-46e8-b9f8-7aeb899bfe7c", + "metadata": {}, + "source": [ + "\n", + "## 4. Loading the catalog into CLMM\n", + "\n", + "Once we have the catalog, we read in the catalog, make cuts on the catalog, and adjust column names to prepare for the analysis in CLMM.\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a04d966e-e00b-4e07-94d9-c3cc13b6183c", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "# Assume the downloaded catalog is at this path:\n", + "#filename = \"ACOS520_DES.csv\"\n", + "#filename = \"ID1.csv\"\n", + "filename = \"ID10.csv\"\n", + "catalog = filename.replace(\".csv\", \".pkl\")\n", + "if not Path(catalog).is_file():\n", + " data_0 = Table.read(filename, format=\"ascii.csv\")\n", + " pkl.dump(data_0, open(catalog, \"wb\"))\n", + "else:\n", + " data_0 = pkl.load(open(catalog, \"rb\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88e10ef1-f035-462b-98c6-ceac1d32a7e1", + "metadata": {}, + "outputs": [], + "source": [ + "print(data_0.colnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c5a40301-a3fd-4b29-a7c7-8bccf4c01c72", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.table import unique\n", + "data_0 = unique(data_0, [\"ra\",\"dec\"])" + ] + }, + { + "cell_type": "markdown", + "id": "e66ea8ec-90f7-4df1-86a1-678a483da776", + "metadata": {}, + "source": [ + "### Shear response\n", + "Shears in the DES data have been measured using the `metacal` method and the catalog provides the shear response terms ($r11,r22, r12, r21$) required to calibrate the shear values. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5df10e3-771d-499e-b462-b72682902b98", + "metadata": {}, + "outputs": [], + "source": [ + "print(np.mean(data_0[\"r11\"]), np.mean(data_0[\"r22\"]))\n", + "print(np.mean(data_0[\"r12\"]), np.mean(data_0[\"r21\"]))\n", + "r_diag = np.mean([np.mean(data_0[\"r11\"]), np.mean(data_0[\"r22\"])])\n", + "r_off_diag = np.mean([np.mean(data_0[\"r12\"]), np.mean(data_0[\"r21\"])])\n", + "print(r_diag, r_off_diag, r_off_diag / r_diag)" + ] + }, + { + "cell_type": "markdown", + "id": "93c5d48a-3325-405e-b20a-e9e862d66975", + "metadata": {}, + "source": [ + "The diagonal terms are close to each other. The off-diagonal terms are much smaller (<1%).\n", + "We use the mean of the diagonal terms to reduce noise.\n", + "We also skip the selection bias since it is typically at percent level." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a98993f-2261-452b-8515-a2c64d52c744", + "metadata": {}, + "outputs": [], + "source": [ + "# Adjust column names.\n", + "def adjust_column_names(catalog_in):\n", + " # We consider a map between new and old column names.\n", + " # Note we have considered shear calibration here.\n", + " column_name_map = {\n", + " \"ra\": \"ra\",\n", + " \"dec\": \"dec\",\n", + " \"z\": \"mean_z\",\n", + " \"e1\": \"e1\",\n", + " \"e2\": \"e2\",\n", + " }\n", + "\n", + " catalog_out = Table()\n", + " for i in column_name_map:\n", + " catalog_out[i] = catalog_in[column_name_map[i]]\n", + "\n", + " catalog_out[\"e1\"] /= r_diag\n", + " catalog_out[\"e2\"] /= r_diag\n", + "\n", + " return catalog_out\n", + "\n", + "\n", + "obs_galaxies = adjust_column_names(data_0)\n", + "\n", + "select = obs_galaxies[\"e1\"] ** 2 + obs_galaxies[\"e2\"] ** 2 <= 1.0\n", + "print(np.sum(~select))\n", + "obs_galaxies = obs_galaxies[select]" + ] + }, + { + "cell_type": "markdown", + "id": "09a7df70-8d40-4106-9137-9a9036fac7fc", + "metadata": {}, + "source": [ + "### Basic visualization" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6d599e88-8a45-4b44-a46a-f2c00d86bd7b", + "metadata": {}, + "outputs": [], + "source": [ + "def make_plots(catalog_in):\n", + " # Scatter plot\n", + " plt.figure()\n", + " plt.scatter(catalog_in[\"ra\"], catalog_in[\"dec\"], c=catalog_in[\"z\"], s=1.0, alpha=1)\n", + " cb = plt.colorbar()\n", + " plt.xlabel(\"ra\")\n", + " plt.ylabel(\"dec\")\n", + " cb.ax.set_title(\"z\")\n", + "\n", + " # Histogram\n", + " plt.figure()\n", + " plt.hist(catalog_in[\"z\"], bins=20)\n", + " plt.xlabel(\"z\")\n", + " plt.ylabel(\"count\")\n", + "\n", + " # Relation\n", + " plt.figure()\n", + " plt.scatter(catalog_in[\"e1\"], catalog_in[\"e2\"], s=1.0, alpha=0.2)\n", + " plt.xlabel(\"e1\")\n", + " plt.ylabel(\"e2\")\n", + "\n", + "\n", + "make_plots(obs_galaxies)" + ] + }, + { + "cell_type": "markdown", + "id": "e71d6d00-f18e-4e27-969b-f02f13032713", + "metadata": {}, + "source": [ + "\n", + "## 5. Running CLMM on the dataset\n", + "We use the functions similar to `examples/Paper_v1.0/gt_and_use_case.ipynb`." + ] + }, + { + "cell_type": "markdown", + "id": "4435913d-4d81-4d7c-b2a1-d828204cd889", + "metadata": {}, + "source": [ + "### Make a galaxy cluster object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c4b465d-9e1c-4b63-a43e-eeea6139b462", + "metadata": {}, + "outputs": [], + "source": [ + "from clmm import Cosmology\n", + "import clmm\n", + "\n", + "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", + "\n", + "# We consider RMJ051637.4-543001.6 (ACO S520)\n", + "#cluster_z = 0.30416065 # Cluster redshift\n", + "#cluster_ra = 79.155704 # Cluster Ra in deg\n", + "#cluster_dec = -54.500456 # Cluster Dec in deg\n", + "\n", + "# ID1\n", + "#cluster_z = 0.4298039972782135 # Cluster redshift\n", + "#cluster_ra = 43.564574 # Cluster Ra in deg\n", + "#cluster_dec = -58.95297 # Cluster Dec in deg\n", + "\n", + "# ID10\n", + "cluster_z = 0.3264264166355133 # Cluster redshift\n", + "cluster_ra = 323.800394 # Cluster Ra in deg\n", + "cluster_dec = -1.04959 # Cluster Dec in deg\n", + "\n", + "# Select background galaxies\n", + "obs_galaxies = obs_galaxies[(obs_galaxies[\"z\"] > (cluster_z + 0.1)) & (obs_galaxies[\"z\"] < 1.5)]\n", + "\n", + "obs_galaxies[\"id\"] = np.arange(len(obs_galaxies))\n", + "\n", + "# Put galaxy values on arrays\n", + "gal_ra = obs_galaxies[\"ra\"] # Galaxies Ra in deg\n", + "gal_dec = obs_galaxies[\"dec\"] # Galaxies Dec in deg\n", + "gal_e1 = obs_galaxies[\"e1\"] # Galaxies elipticipy 1\n", + "gal_e2 = obs_galaxies[\"e2\"] # Galaxies elipticipy 2\n", + "gal_z = obs_galaxies[\"z\"] # Galaxies observed redshift\n", + "gal_id = obs_galaxies[\"id\"] # Galaxies ID\n", + "\n", + "# Create a GCData with the galaxies.\n", + "galaxies = clmm.GCData(\n", + " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id], names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"]\n", + ")\n", + "\n", + "# Create a GalaxyCluster.\n", + "cluster = clmm.GalaxyCluster(\"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies)" + ] + }, + { + "cell_type": "markdown", + "id": "4d526e90-3c58-4f4e-ad9b-58c95464d72d", + "metadata": {}, + "source": [ + "### Measure the shear profile" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "42573e50-85cc-4977-baf4-723711792c6b", + "metadata": {}, + "outputs": [], + "source": [ + "import clmm.dataops as da\n", + "\n", + "# Convert ellipticities into shears for the members.\n", + "cluster.compute_tangential_and_cross_components()\n", + "print(cluster.galcat.colnames)\n", + "\n", + "# Measure profile and add profile table to the cluster.\n", + "cluster.make_radial_profile(\n", + " bins=da.make_bins(0.2, 4.0, 6, method=\"evenlog10width\"),\n", + " bin_units=\"Mpc\",\n", + " cosmo=cosmo,\n", + " include_empty_bins=False,\n", + " gal_ids_in_bins=True,\n", + ")\n", + "print(cluster.profile.colnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db4ebc79-67e8-43da-ad89-ca96666f541d", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(figsize=(7, 5), ncols=1, nrows=1)\n", + "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", + "ax.errorbar(\n", + " cluster.profile[\"radius\"],\n", + " cluster.profile[\"gt\"],\n", + " cluster.profile[\"gt_err\"],\n", + " c=\"k\",\n", + " **errorbar_kwargs\n", + ")\n", + "ax.set_xlabel(\"R [Mpc]\", fontsize=10)\n", + "ax.set_ylabel(r\"$g_t$\", fontsize=10)\n", + "ax.set_xscale(\"log\")\n", + "ax.grid(lw=0.3)\n", + "ax.minorticks_on()\n", + "ax.grid(which=\"minor\", lw=0.1)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "9cccdd77-3312-4780-a202-0b88c3aa9802", + "metadata": {}, + "source": [ + "### Theoretical predictions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87e51a27-936a-48f5-abc7-120cff89b3ac", + "metadata": {}, + "outputs": [], + "source": [ + "from clmm.utils import convert_units\n", + "\n", + "# Model relying on the overall redshift distribution of the sources (WtG III Applegate et al. 2014).\n", + "z_inf = 1000\n", + "concentration = 4.0\n", + "\n", + "bs_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo))\n", + "bs2_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo) ** 2)\n", + "\n", + "\n", + "def predict_reduced_tangential_shear_redshift_distribution(profile, logm):\n", + " gt = clmm.compute_reduced_tangential_shear(\n", + " r_proj=profile[\"radius\"], # Radial component of the profile\n", + " mdelta=10**logm, # Mass of the cluster [M_sun]\n", + " cdelta=concentration, # Concentration of the cluster\n", + " z_cluster=cluster_z, # Redshift of the cluster\n", + " z_src=(bs_mean, bs2_mean), # tuple of (bs_mean, bs2_mean)\n", + " z_src_info=\"beta\",\n", + " approx=\"order1\",\n", + " cosmo=cosmo,\n", + " delta_mdef=200,\n", + " massdef=\"critical\",\n", + " halo_profile_model=\"nfw\",\n", + " )\n", + " return gt\n", + "\n", + "\n", + "# Model using individual redshift and radial information, to compute the averaged shear in each radial bin, based on the galaxies actually present in that bin.\n", + "cluster.galcat[\"theta_mpc\"] = convert_units(\n", + " cluster.galcat[\"theta\"], \"radians\", \"mpc\", cluster.z, cosmo\n", + ")\n", + "\n", + "\n", + "def predict_reduced_tangential_shear_individual_redshift(profile, logm):\n", + " return np.array(\n", + " [\n", + " np.mean(\n", + " clmm.compute_reduced_tangential_shear(\n", + " # Radial component of each source galaxy inside the radial bin\n", + " r_proj=cluster.galcat[radial_bin[\"gal_id\"]][\"theta_mpc\"],\n", + " mdelta=10**logm, # Mass of the cluster [M_sun]\n", + " cdelta=concentration, # Concentration of the cluster\n", + " z_cluster=cluster_z, # Redshift of the cluster\n", + " # Redshift value of each source galaxy inside the radial bin\n", + " z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n", + " cosmo=cosmo,\n", + " delta_mdef=200,\n", + " massdef=\"critical\",\n", + " halo_profile_model=\"nfw\",\n", + " )\n", + " )\n", + " for radial_bin in profile\n", + " ]\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "e2fd16d5-6d8b-4e9b-a0c3-5505af0291b8", + "metadata": {}, + "source": [ + "### Mass fitting" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c83799a1-5122-4e89-b67d-8f691ff9d541", + "metadata": {}, + "outputs": [], + "source": [ + "mask_for_fit = cluster.profile[\"n_src\"] > 2\n", + "data_for_fit = cluster.profile[mask_for_fit]\n", + "\n", + "from clmm.support.sampler import fitters\n", + "\n", + "\n", + "def fit_mass(predict_function):\n", + " popt, pcov = fitters[\"curve_fit\"](\n", + " predict_function,\n", + " data_for_fit,\n", + " data_for_fit[\"gt\"],\n", + " data_for_fit[\"gt_err\"],\n", + " bounds=[10.0, 17.0],\n", + " )\n", + " logm, logm_err = popt[0], np.sqrt(pcov[0][0])\n", + " return {\n", + " \"logm\": logm,\n", + " \"logm_err\": logm_err,\n", + " \"m\": 10**logm,\n", + " \"m_err\": (10**logm) * logm_err * np.log(10),\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc59846f-a380-449d-a30b-5261085a34f2", + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "fit_redshift_distribution = fit_mass(predict_reduced_tangential_shear_redshift_distribution)\n", + "fit_individual_redshift = fit_mass(predict_reduced_tangential_shear_individual_redshift)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f5b8ac6-7294-4e85-97ac-28b62512be60", + "metadata": {}, + "outputs": [], + "source": [ + "print(\n", + " \"Best fit mass for N(z) model =\"\n", + " f' {fit_redshift_distribution[\"m\"]:.3e} +/- {fit_redshift_distribution[\"m_err\"]:.3e} Msun'\n", + ")\n", + "print(\n", + " \"Best fit mass for individual redshift and radius =\"\n", + " f' {fit_individual_redshift[\"m\"]:.3e} +/- {fit_individual_redshift[\"m_err\"]:.3e} Msun'\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "abd490e7-9404-4daa-8d79-f41a638d86da", + "metadata": {}, + "source": [ + "### Visualization of the results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84e71bf0-5aa8-4c48-8042-bd461430291a", + "metadata": {}, + "outputs": [], + "source": [ + "def get_predicted_shear(predict_function, fit_values):\n", + " gt_est = predict_function(data_for_fit, fit_values[\"logm\"])\n", + " gt_est_err = [\n", + " predict_function(data_for_fit, fit_values[\"logm\"] + i * fit_values[\"logm_err\"])\n", + " for i in (-3, 3)\n", + " ]\n", + " return gt_est, gt_est_err\n", + "\n", + "\n", + "gt_redshift_distribution, gt_err_redshift_distribution = get_predicted_shear(\n", + " predict_reduced_tangential_shear_redshift_distribution, fit_redshift_distribution\n", + ")\n", + "gt_individual_redshift, gt_err_individual_redshift = get_predicted_shear(\n", + " predict_reduced_tangential_shear_individual_redshift, fit_individual_redshift\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d0e95d3-bc69-4cf7-8ac0-1788b74fe5d6", + "metadata": {}, + "outputs": [], + "source": [ + "chi2_redshift_distribution_dof = np.sum(\n", + " (gt_redshift_distribution - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", + ") / (len(data_for_fit) - 1)\n", + "chi2_individual_redshift_dof = np.sum(\n", + " (gt_individual_redshift - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", + ") / (len(data_for_fit) - 1)\n", + "\n", + "print(f\"Reduced chi2 (N(z) model) = {chi2_redshift_distribution_dof}\")\n", + "print(f\"Reduced chi2 (individual (R,z) model) = {chi2_individual_redshift_dof}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf89c28a-c19d-4022-83ce-32c113d2a615", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = plt.subplots(\n", + " nrows=2,\n", + " ncols=1,\n", + " #figsize=(8, 6),\n", + " figsize=(8*0.8, 6*0.8),\n", + " gridspec_kw={\"height_ratios\": [3, 1], \"wspace\": 0.4, \"hspace\": 0.03},\n", + ")\n", + "\n", + "axes[0].errorbar(\n", + " data_for_fit[\"radius\"], data_for_fit[\"gt\"], data_for_fit[\"gt_err\"], c=\"k\", **errorbar_kwargs\n", + ")\n", + "\n", + "# Points in grey have not been used for the fit.\n", + "axes[0].errorbar(\n", + " cluster.profile[\"radius\"][~mask_for_fit],\n", + " cluster.profile[\"gt\"][~mask_for_fit],\n", + " cluster.profile[\"gt_err\"][~mask_for_fit],\n", + " c=\"grey\",\n", + " **errorbar_kwargs,\n", + ")\n", + "\n", + "pow10 = 15\n", + "mlabel = lambda name, fits: (\n", + " rf\"$M_{{fit}}^{{{name}}} = \"\n", + " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", + " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", + " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", + ")\n", + "\n", + "# The model for the 1st method.\n", + "axes[0].loglog(\n", + " data_for_fit[\"radius\"],\n", + " gt_redshift_distribution,\n", + " \"-C1\",\n", + " label=mlabel(\"N(z)\", fit_redshift_distribution),\n", + " lw=0.5,\n", + ")\n", + "\n", + "axes[0].fill_between(\n", + " data_for_fit[\"radius\"], *gt_err_redshift_distribution, lw=0, color=\"C1\", alpha=0.2\n", + ")\n", + "\n", + "# The model for the 2nd method.\n", + "axes[0].loglog(\n", + " data_for_fit[\"radius\"],\n", + " gt_individual_redshift,\n", + " \"-C2\",\n", + " label=mlabel(\"z,R\", fit_individual_redshift),\n", + " lw=0.5,\n", + ")\n", + "axes[0].fill_between(\n", + " data_for_fit[\"radius\"], *gt_err_individual_redshift, lw=0, color=\"C2\", alpha=0.2\n", + ")\n", + "\n", + "\n", + "axes[0].set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", + "axes[0].legend(fontsize=\"x-large\", loc=\"lower left\")\n", + "axes[0].set_xticklabels([])\n", + "axes[0].tick_params(\"x\", labelsize=\"x-large\")\n", + "axes[0].tick_params(\"y\", labelsize=\"x-large\")\n", + "axes[1].set_ylim(1.0e-3, 0.5)\n", + "\n", + "\n", + "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", + "errorbar_kwargs2[\"markersize\"] = 5\n", + "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "\n", + "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.15\n", + "\n", + "\n", + "axes[1].errorbar(\n", + " data_for_fit[\"radius\"],\n", + " data_for_fit[\"gt\"] / gt_redshift_distribution - 1,\n", + " yerr=data_for_fit[\"gt_err\"] / gt_redshift_distribution,\n", + " marker=\"s\",\n", + " c=\"C1\",\n", + " **errorbar_kwargs2,\n", + ")\n", + "#errorbar_kwargs2[\"markersize\"] = 3\n", + "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "\n", + "axes[1].errorbar(\n", + " data_for_fit[\"radius\"] * delta,\n", + " data_for_fit[\"gt\"] / gt_individual_redshift - 1,\n", + " yerr=data_for_fit[\"gt_err\"] / gt_individual_redshift,\n", + " marker=\"*\",\n", + " c=\"C2\",\n", + " **errorbar_kwargs2,\n", + ")\n", + "axes[1].set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", + "\n", + "axes[1].set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", + "axes[1].set_xscale(\"log\")\n", + "\n", + "axes[1].set_ylim(-5, 5)\n", + "\n", + "axes[1].tick_params(\"x\", labelsize=\"x-large\")\n", + "axes[1].tick_params(\"y\", labelsize=\"x-large\")\n", + "\n", + "for ax in axes:\n", + " ax.grid(lw=0.3, ls=':')\n", + " ax.minorticks_on()\n", + " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", + "plt.show()\n", + "\n", + "# Note since we made cuts on the catalog, the redshift distribution of the remaining sources might not be representative." + ] + }, + { + "cell_type": "markdown", + "id": "4e50e37a-5c53-4420-9453-310aed730df0", + "metadata": {}, + "source": [ + "## References\n", + "\n", + "Zuntz J., Sheldon E., Samuroff S., Troxel M. A., Jarvis M., MacCrann N., Gruen D., et al., 2018, MNRAS, 481, 1149. [doi:10.1093/mnras/sty2219](\n", + "http://adsabs.harvard.edu/cgi-bin/bib_query?arXiv:1708.01533)\n", + "\n", + "Hoyle B., Gruen D., Bernstein G. M., Rau M. M., De Vicente J., Hartley W. G., Gaztanaga E., et al., 2018, MNRAS, 478, 592. [doi:10.1093/mnras/sty957](http://adsabs.harvard.edu/abs/2018MNRAS.478..592H)\n", + "\n", + "McClintock T., Varga T. N., Gruen D., Rozo E., Rykoff E. S., Shin T., Melchior P., et al., 2019, MNRAS, 482, 1352. [doi:10.1093/mnras/sty2711](https://ui.adsabs.harvard.edu/abs/2019MNRAS.482.1352M/abstract)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "conda-clmmenv", + "language": "python", + "name": "conda-clmmenv" + }, + "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.9.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 272ba947b8decccb0b9d1f1c7ddda3e3029218c1 Mon Sep 17 00:00:00 2001 From: Shenming Fu Date: Tue, 19 Aug 2025 02:09:01 -0700 Subject: [PATCH 15/17] Revert "add nb for des fig" This reverts commit b3fda3afb2c06ffc408b8a5e615319124f5889f8. --- .../Example5_Fit_Halo_mass_to_DES_data.ipynb | 76 +++++++------------ 1 file changed, 26 insertions(+), 50 deletions(-) diff --git a/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb b/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb index cf18a85a5..3e24d6f69 100644 --- a/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb +++ b/examples/mass_fitting/Example5_Fit_Halo_mass_to_DES_data.ipynb @@ -131,9 +131,7 @@ "source": [ "%%time\n", "# Assume the downloaded catalog is at this path:\n", - "#filename = \"ACOS520_DES.csv\"\n", - "#filename = \"ID1.csv\"\n", - "filename = \"ID10.csv\"\n", + "filename = \"ACOS520_DES.csv\"\n", "catalog = filename.replace(\".csv\", \".pkl\")\n", "if not Path(catalog).is_file():\n", " data_0 = Table.read(filename, format=\"ascii.csv\")\n", @@ -152,17 +150,6 @@ "print(data_0.colnames)" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "c5a40301-a3fd-4b29-a7c7-8bccf4c01c72", - "metadata": {}, - "outputs": [], - "source": [ - "from astropy.table import unique\n", - "data_0 = unique(data_0, [\"ra\",\"dec\"])" - ] - }, { "cell_type": "markdown", "id": "e66ea8ec-90f7-4df1-86a1-678a483da776", @@ -303,19 +290,9 @@ "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", "\n", "# We consider RMJ051637.4-543001.6 (ACO S520)\n", - "#cluster_z = 0.30416065 # Cluster redshift\n", - "#cluster_ra = 79.155704 # Cluster Ra in deg\n", - "#cluster_dec = -54.500456 # Cluster Dec in deg\n", - "\n", - "# ID1\n", - "#cluster_z = 0.4298039972782135 # Cluster redshift\n", - "#cluster_ra = 43.564574 # Cluster Ra in deg\n", - "#cluster_dec = -58.95297 # Cluster Dec in deg\n", - "\n", - "# ID10\n", - "cluster_z = 0.3264264166355133 # Cluster redshift\n", - "cluster_ra = 323.800394 # Cluster Ra in deg\n", - "cluster_dec = -1.04959 # Cluster Dec in deg\n", + "cluster_z = 0.30416065 # Cluster redshift\n", + "cluster_ra = 79.155704 # Cluster Ra in deg\n", + "cluster_dec = -54.500456 # Cluster Dec in deg\n", "\n", "# Select background galaxies\n", "obs_galaxies = obs_galaxies[(obs_galaxies[\"z\"] > (cluster_z + 0.1)) & (obs_galaxies[\"z\"] < 1.5)]\n", @@ -362,7 +339,7 @@ "\n", "# Measure profile and add profile table to the cluster.\n", "cluster.make_radial_profile(\n", - " bins=da.make_bins(0.2, 4.0, 6, method=\"evenlog10width\"),\n", + " bins=da.make_bins(0.2, 5.0, 7, method=\"evenlog10width\"),\n", " bin_units=\"Mpc\",\n", " cosmo=cosmo,\n", " include_empty_bins=False,\n", @@ -379,7 +356,7 @@ "outputs": [], "source": [ "fig, ax = plt.subplots(figsize=(7, 5), ncols=1, nrows=1)\n", - "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", + "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=1, elinewidth=0.5, capthick=0.5)\n", "ax.errorbar(\n", " cluster.profile[\"radius\"],\n", " cluster.profile[\"gt\"],\n", @@ -594,8 +571,7 @@ "fig, axes = plt.subplots(\n", " nrows=2,\n", " ncols=1,\n", - " #figsize=(8, 6),\n", - " figsize=(8*0.8, 6*0.8),\n", + " figsize=(8, 6),\n", " gridspec_kw={\"height_ratios\": [3, 1], \"wspace\": 0.4, \"hspace\": 0.03},\n", ")\n", "\n", @@ -612,11 +588,11 @@ " **errorbar_kwargs,\n", ")\n", "\n", - "pow10 = 15\n", + "pow10 = 14\n", "mlabel = lambda name, fits: (\n", " rf\"$M_{{fit}}^{{{name}}} = \"\n", - " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", - " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", + " rf'{fits[\"m\"]/10**pow10:.3f}\\pm'\n", + " rf'{fits[\"m_err\"]/10**pow10:.3f}'\n", " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", ")\n", "\n", @@ -646,16 +622,16 @@ ")\n", "\n", "\n", - "axes[0].set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", - "axes[0].legend(fontsize=\"x-large\", loc=\"lower left\")\n", + "axes[0].set_ylabel(r\"$g_t$\", fontsize=12)\n", + "axes[0].legend(fontsize=12, loc=4)\n", "axes[0].set_xticklabels([])\n", - "axes[0].tick_params(\"x\", labelsize=\"x-large\")\n", - "axes[0].tick_params(\"y\", labelsize=\"x-large\")\n", + "axes[0].tick_params(\"x\", labelsize=12)\n", + "axes[0].tick_params(\"y\", labelsize=12)\n", "axes[1].set_ylim(1.0e-3, 0.5)\n", "\n", "\n", "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", - "errorbar_kwargs2[\"markersize\"] = 5\n", + "errorbar_kwargs2[\"markersize\"] = 3\n", "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "\n", "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.15\n", @@ -669,8 +645,8 @@ " c=\"C1\",\n", " **errorbar_kwargs2,\n", ")\n", - "#errorbar_kwargs2[\"markersize\"] = 3\n", - "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "errorbar_kwargs2[\"markersize\"] = 3\n", + "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "\n", "axes[1].errorbar(\n", " data_for_fit[\"radius\"] * delta,\n", @@ -680,20 +656,20 @@ " c=\"C2\",\n", " **errorbar_kwargs2,\n", ")\n", - "axes[1].set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", + "axes[1].set_xlabel(r\"$R$ [Mpc]\", fontsize=12)\n", "\n", - "axes[1].set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", + "axes[1].set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=12)\n", "axes[1].set_xscale(\"log\")\n", "\n", "axes[1].set_ylim(-5, 5)\n", "\n", - "axes[1].tick_params(\"x\", labelsize=\"x-large\")\n", - "axes[1].tick_params(\"y\", labelsize=\"x-large\")\n", + "axes[1].tick_params(\"x\", labelsize=12)\n", + "axes[1].tick_params(\"y\", labelsize=12)\n", "\n", "for ax in axes:\n", - " ax.grid(lw=0.3, ls=':')\n", + " ax.grid(lw=0.3)\n", " ax.minorticks_on()\n", - " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", + " ax.grid(which=\"minor\", lw=0.1)\n", "plt.show()\n", "\n", "# Note since we made cuts on the catalog, the redshift distribution of the remaining sources might not be representative." @@ -717,9 +693,9 @@ ], "metadata": { "kernelspec": { - "display_name": "conda-clmmenv", + "display_name": "wrk", "language": "python", - "name": "conda-clmmenv" + "name": "wrk" }, "language_info": { "codemirror_mode": { @@ -731,7 +707,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.9" + "version": "3.10.6" } }, "nbformat": 4, From 193be065be09e0164818233923bb56aba03104da Mon Sep 17 00:00:00 2001 From: Shenming Fu Date: Tue, 19 Aug 2025 02:24:55 -0700 Subject: [PATCH 16/17] revert accidental change to mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb --- .../Example4_Fit_Halo_mass_to_HSC_data.ipynb | 99 ++++++------------- 1 file changed, 31 insertions(+), 68 deletions(-) diff --git a/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb b/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb index 859094468..74502f5c7 100644 --- a/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb +++ b/examples/mass_fitting/Example4_Fit_Halo_mass_to_HSC_data.ipynb @@ -2,7 +2,6 @@ "cells": [ { "cell_type": "markdown", - "id": "74785dcc", "metadata": { "tags": [] }, @@ -43,7 +42,6 @@ }, { "cell_type": "markdown", - "id": "8b2a5bbe", "metadata": {}, "source": [ "\n", @@ -55,7 +53,6 @@ { "cell_type": "code", "execution_count": null, - "id": "89a38de9", "metadata": {}, "outputs": [], "source": [ @@ -70,7 +67,6 @@ }, { "cell_type": "markdown", - "id": "cf169f3c", "metadata": {}, "source": [ "\n", @@ -93,7 +89,6 @@ }, { "cell_type": "markdown", - "id": "71f2b585", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] @@ -155,7 +150,6 @@ }, { "cell_type": "markdown", - "id": "eba97e69", "metadata": {}, "source": [ "\n", @@ -169,7 +163,6 @@ { "cell_type": "code", "execution_count": null, - "id": "fc38b936", "metadata": {}, "outputs": [], "source": [ @@ -187,7 +180,6 @@ { "cell_type": "code", "execution_count": null, - "id": "03e1ca9a", "metadata": {}, "outputs": [], "source": [ @@ -197,7 +189,6 @@ { "cell_type": "code", "execution_count": null, - "id": "d8da2cdb", "metadata": {}, "outputs": [], "source": [ @@ -210,7 +201,6 @@ { "cell_type": "code", "execution_count": null, - "id": "bb7d297a", "metadata": {}, "outputs": [], "source": [ @@ -244,7 +234,6 @@ { "cell_type": "code", "execution_count": null, - "id": "bac2fc15", "metadata": {}, "outputs": [], "source": [ @@ -252,7 +241,7 @@ "# Section A.3.2: \"per-object galaxy shear estimate\".\n", "def apply_shear_calibration(catalog_in):\n", " e1_0 = catalog_in[\"ishape_hsm_regauss_e1\"]\n", - " e2_0 = -catalog_in[\"ishape_hsm_regauss_e2\"]\n", + " e2_0 = catalog_in[\"ishape_hsm_regauss_e2\"]\n", " e_rms = catalog_in[\"ishape_hsm_regauss_derived_rms_e\"]\n", " m = catalog_in[\"ishape_hsm_regauss_derived_shear_bias_m\"]\n", " c1 = catalog_in[\"ishape_hsm_regauss_derived_shear_bias_c1\"]\n", @@ -275,7 +264,6 @@ { "cell_type": "code", "execution_count": null, - "id": "cff6ac8e", "metadata": {}, "outputs": [], "source": [ @@ -308,7 +296,6 @@ { "cell_type": "code", "execution_count": null, - "id": "a5315cfc", "metadata": {}, "outputs": [], "source": [ @@ -340,7 +327,6 @@ }, { "cell_type": "markdown", - "id": "5dc4163b", "metadata": {}, "source": [ "\n", @@ -351,7 +337,6 @@ { "cell_type": "code", "execution_count": null, - "id": "e469558f", "metadata": {}, "outputs": [], "source": [ @@ -360,19 +345,9 @@ "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", "\n", "# We consider MaxBCG J140.53188+03.76632\n", - "#cluster_z = 0.2701 # Cluster redshift\n", - "#cluster_ra = 140.54565\n", - "#cluster_dec = 3.77820\n", - "\n", - "#HWL16a-026\n", - "#cluster_z = 0.424\n", - "#cluster_ra = 130.5895\n", - "#cluster_dec = 1.6473\n", - "\n", - "#HWL16a-034\n", - "cluster_z = 0.315\n", - "cluster_ra = 139.0387\n", - "cluster_dec = -0.3966\n", + "cluster_z = 0.2701 # Cluster redshift\n", + "cluster_ra = 140.54565 # Cluster Ra in deg\n", + "cluster_dec = 3.77820 # Cluster Dec in deg\n", "\n", "obs_galaxies = data_2\n", "\n", @@ -399,7 +374,6 @@ { "cell_type": "code", "execution_count": null, - "id": "a12e56be", "metadata": {}, "outputs": [], "source": [ @@ -437,13 +411,12 @@ { "cell_type": "code", "execution_count": null, - "id": "89050c01", "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(6, 4))\n", "ax = fig.add_axes((0, 0, 1, 1))\n", - "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", + "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=1, elinewidth=0.5, capthick=0.5)\n", "ax.errorbar(\n", " cluster.profile[\"radius\"],\n", " cluster.profile[\"gt\"],\n", @@ -456,14 +429,12 @@ "ax.grid(lw=0.3)\n", "ax.minorticks_on()\n", "ax.grid(which=\"minor\", lw=0.1)\n", - "ax.set_xscale('log')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, - "id": "a806b0f1", "metadata": {}, "outputs": [], "source": [ @@ -530,7 +501,6 @@ { "cell_type": "code", "execution_count": null, - "id": "48dd5752", "metadata": {}, "outputs": [], "source": [ @@ -562,7 +532,6 @@ { "cell_type": "code", "execution_count": null, - "id": "9e6a4e04", "metadata": {}, "outputs": [], "source": [ @@ -584,7 +553,6 @@ { "cell_type": "code", "execution_count": null, - "id": "dc3a7934", "metadata": {}, "outputs": [], "source": [ @@ -596,7 +564,6 @@ { "cell_type": "code", "execution_count": null, - "id": "79e35734", "metadata": {}, "outputs": [], "source": [ @@ -611,7 +578,6 @@ { "cell_type": "code", "execution_count": null, - "id": "73202bf8", "metadata": {}, "outputs": [], "source": [ @@ -643,7 +609,6 @@ { "cell_type": "code", "execution_count": null, - "id": "ec85a1b1", "metadata": {}, "outputs": [], "source": [ @@ -661,7 +626,6 @@ { "cell_type": "code", "execution_count": null, - "id": "5a7ca7d4", "metadata": {}, "outputs": [], "source": [ @@ -713,25 +677,25 @@ "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_individual_redshift, lw=0, color=\"C2\", alpha=0.2)\n", "\n", "# The value in the reference paper.\n", - "#gt_ax.loglog(\n", - "# data_for_fit[\"radius\"], gt_paper1, \"-C3\", label=mlabel(\"paper; N(z)\", paper_value), lw=0.5\n", - "#)\n", - "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper1, lw=0, color=\"C3\", alpha=0.2)\n", - "#\n", - "#gt_ax.loglog(\n", - "# data_for_fit[\"radius\"], gt_paper2, \"-C4\", label=mlabel(\"paper; Z,R\", paper_value), lw=0.5\n", - "#)\n", - "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper2, lw=0, color=\"C4\", alpha=0.2)\n", - "\n", - "gt_ax.set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", - "gt_ax.legend(fontsize=\"x-large\")\n", + "gt_ax.loglog(\n", + " data_for_fit[\"radius\"], gt_paper1, \"-C3\", label=mlabel(\"paper; N(z)\", paper_value), lw=0.5\n", + ")\n", + "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper1, lw=0, color=\"C3\", alpha=0.2)\n", + "\n", + "gt_ax.loglog(\n", + " data_for_fit[\"radius\"], gt_paper2, \"-C4\", label=mlabel(\"paper; Z,R\", paper_value), lw=0.5\n", + ")\n", + "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper2, lw=0, color=\"C4\", alpha=0.2)\n", + "\n", + "gt_ax.set_ylabel(r\"$g_t$\", fontsize=8)\n", + "gt_ax.legend(fontsize=6)\n", "gt_ax.set_xticklabels([])\n", - "gt_ax.tick_params(\"x\", labelsize=\"x-large\")\n", - "gt_ax.tick_params(\"y\", labelsize=\"x-large\")\n", + "gt_ax.tick_params(\"x\", labelsize=8)\n", + "gt_ax.tick_params(\"y\", labelsize=8)\n", "\n", "\n", "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", - "errorbar_kwargs2[\"markersize\"] = 5\n", + "errorbar_kwargs2[\"markersize\"] = 3\n", "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "res_ax = fig.add_axes((0.25, 0.2, 0.7, 0.2))\n", "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.25\n", @@ -745,8 +709,8 @@ " c=\"C1\",\n", " **errorbar_kwargs2,\n", ")\n", - "#errorbar_kwargs2[\"markersize\"] = 3\n", - "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "errorbar_kwargs2[\"markersize\"] = 3\n", + "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", "\n", "res_ax.errorbar(\n", " data_for_fit[\"radius\"] * delta,\n", @@ -756,27 +720,26 @@ " c=\"C2\",\n", " **errorbar_kwargs2,\n", ")\n", - "res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", + "res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=8)\n", "\n", - "res_ax.set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", + "res_ax.set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=8)\n", "res_ax.set_xscale(\"log\")\n", "\n", "res_ax.set_ylim(-5, 5)\n", "res_ax.yaxis.set_minor_locator(MultipleLocator(10))\n", "\n", - "res_ax.tick_params(\"x\", labelsize=\"x-large\")\n", - "res_ax.tick_params(\"y\", labelsize=\"x-large\")\n", + "res_ax.tick_params(\"x\", labelsize=8)\n", + "res_ax.tick_params(\"y\", labelsize=8)\n", "\n", "for ax in (gt_ax, res_ax):\n", - " ax.grid(lw=0.3, ls=':')\n", + " ax.grid(lw=0.3)\n", " ax.minorticks_on()\n", - " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", + " ax.grid(which=\"minor\", lw=0.1)\n", "plt.show()" ] }, { "cell_type": "markdown", - "id": "bc63da5c", "metadata": {}, "source": [ "### Reference\n", @@ -806,9 +769,9 @@ ], "metadata": { "kernelspec": { - "display_name": "conda-clmmenv", + "display_name": "Python 3", "language": "python", - "name": "conda-clmmenv" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -820,7 +783,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.9" + "version": "3.6.9" } }, "nbformat": 4, From cb35123bbf325f0fe7a44f11db4bf949b54449a2 Mon Sep 17 00:00:00 2001 From: Michel Date: Wed, 17 Dec 2025 22:26:08 +0100 Subject: [PATCH 17/17] update data nbs --- ...le4_Fit_Halo_mass_to_HSC_data_paper2.ipynb | 588 +++++++++--------- ...le5_Fit_Halo_mass_to_DES_data_paper2.ipynb | 572 +++++++++-------- 2 files changed, 609 insertions(+), 551 deletions(-) diff --git a/examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb b/examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb index 859094468..ae291b6f7 100644 --- a/examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb +++ b/examples/Paper_v2.0/Example4_Fit_Halo_mass_to_HSC_data_paper2.ipynb @@ -59,13 +59,14 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", + "import pickle as pkl\n", + "from pathlib import Path\n", + "\n", "import matplotlib.pyplot as plt\n", + "import numpy as np\n", "\n", "# %matplotlib inline\n", - "from astropy.table import Table\n", - "import pickle as pkl\n", - "from pathlib import Path" + "from astropy.table import Table" ] }, { @@ -169,19 +170,14 @@ { "cell_type": "code", "execution_count": null, - "id": "fc38b936", + "id": "d8dd8322-f639-4500-8624-7e973ab9294d", "metadata": {}, "outputs": [], "source": [ "%%time\n", "# Assume the downloaded catalog is at this path:\n", - "filename = \"197376_GAMMA09H.csv\"\n", - "catalog = filename.replace(\".csv\", \".pkl\")\n", - "if not Path(catalog).is_file():\n", - " data_0 = Table.read(filename, format=\"ascii.csv\")\n", - " pkl.dump(data_0, open(catalog, \"wb\"))\n", - "else:\n", - " data_0 = pkl.load(open(catalog, \"rb\"))" + "filename = \"/sps/lsst/groups/clusters/CLMM/HSC_data_example/197376_GAMMA09H.csv\"\n", + "data_0 = Table.read(filename, format=\"ascii.csv\")" ] }, { @@ -222,11 +218,15 @@ " select &= catalog_in[\"iclassification_extendedness\"] > 0.5\n", " select &= catalog_in[\"icmodel_mag_err\"] <= 2.5 / np.log(10.0) / 10.0\n", " select &= (\n", - " catalog_in[\"ishape_hsm_regauss_e1\"] ** 2 + catalog_in[\"ishape_hsm_regauss_e2\"] ** 2 < 4.0\n", + " catalog_in[\"ishape_hsm_regauss_e1\"] ** 2\n", + " + catalog_in[\"ishape_hsm_regauss_e2\"] ** 2\n", + " < 4.0\n", " )\n", " select &= catalog_in[\"icmodel_mag\"] <= 24.5\n", " select &= catalog_in[\"iblendedness_abs_flux\"] < (10 ** (-0.375))\n", - " select &= catalog_in[\"ishape_hsm_regauss_resolution\"] >= 0.3 # similar to extendedness\n", + " select &= (\n", + " catalog_in[\"ishape_hsm_regauss_resolution\"] >= 0.3\n", + " ) # similar to extendedness\n", " select &= catalog_in[\"ishape_hsm_regauss_sigma\"] <= 0.4\n", " # Note \"zbest\" minimizes the risk of the photo-z point estimate being far away from the true value.\n", " # Details: https://hsc-release.mtk.nao.ac.jp/doc/wp-content/uploads/2017/02/pdr1_photoz_release_note.pdf\n", @@ -316,7 +316,9 @@ "def make_plots(catalog_in):\n", " # Scatter plot\n", " plt.figure()\n", - " plt.scatter(catalog_in[\"ra\"], catalog_in[\"dec\"], c=catalog_in[\"z\"], s=1.0, alpha=0.2)\n", + " plt.scatter(\n", + " catalog_in[\"ra\"], catalog_in[\"dec\"], c=catalog_in[\"z\"], s=1.0, alpha=0.2\n", + " )\n", " plt.colorbar()\n", " plt.xlabel(\"ra\")\n", " plt.ylabel(\"dec\")\n", @@ -351,55 +353,19 @@ { "cell_type": "code", "execution_count": null, - "id": "e469558f", + "id": "2ae7b3db-1246-493a-89af-61daf50d06d9", "metadata": {}, "outputs": [], "source": [ "from clmm import Cosmology\n", "\n", - "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", - "\n", - "# We consider MaxBCG J140.53188+03.76632\n", - "#cluster_z = 0.2701 # Cluster redshift\n", - "#cluster_ra = 140.54565\n", - "#cluster_dec = 3.77820\n", - "\n", - "#HWL16a-026\n", - "#cluster_z = 0.424\n", - "#cluster_ra = 130.5895\n", - "#cluster_dec = 1.6473\n", - "\n", - "#HWL16a-034\n", - "cluster_z = 0.315\n", - "cluster_ra = 139.0387\n", - "cluster_dec = -0.3966\n", - "\n", - "obs_galaxies = data_2\n", - "\n", - "obs_galaxies = obs_galaxies[(obs_galaxies[\"z\"] > (cluster_z + 0.1))]\n", - "\n", - "# Area cut: the query is made for the whole field and this can simplify the processing.\n", - "select = obs_galaxies[\"ra\"] < cluster_ra + 12.0 / 60.0 / np.cos(cluster_dec / 180.0 * np.pi)\n", - "select &= obs_galaxies[\"ra\"] > cluster_ra - 12.0 / 60.0 / np.cos(cluster_dec / 180.0 * np.pi)\n", - "select &= obs_galaxies[\"dec\"] < cluster_dec + 12.0 / 60.0\n", - "select &= obs_galaxies[\"dec\"] > cluster_dec - 12.0 / 60.0\n", - "obs_galaxies = obs_galaxies[select]\n", - "\n", - "obs_galaxies[\"id\"] = np.arange(len(obs_galaxies))\n", - "\n", - "# Put galaxy values on arrays.\n", - "gal_ra = obs_galaxies[\"ra\"] # Galaxies Ra in deg\n", - "gal_dec = obs_galaxies[\"dec\"] # Galaxies Dec in deg\n", - "gal_e1 = obs_galaxies[\"e1\"] # Galaxies elipticipy 1\n", - "gal_e2 = obs_galaxies[\"e2\"] # Galaxies elipticipy 2\n", - "gal_z = obs_galaxies[\"z\"] # Galaxies observed redshift\n", - "gal_id = obs_galaxies[\"id\"] # Galaxies ID" + "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" ] }, { "cell_type": "code", "execution_count": null, - "id": "a12e56be", + "id": "a5d455cf-9657-4154-ae4a-9983748522e8", "metadata": {}, "outputs": [], "source": [ @@ -407,57 +373,86 @@ "\n", "import clmm\n", "import clmm.dataops as da\n", - "from clmm.utils import convert_units\n", - "\n", - "# Create a GCData with the galaxies.\n", - "galaxies = clmm.GCData(\n", - " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id], names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"]\n", - ")\n", - "\n", - "# Create a GalaxyCluster.\n", - "cluster = clmm.GalaxyCluster(\"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies)\n", - "\n", - "# Convert elipticities into shears for the members.\n", - "cluster.compute_tangential_and_cross_components()\n", - "print(cluster.galcat.colnames)\n", - "\n", - "# Measure profile and add profile table to the cluster.\n", - "seps = convert_units(cluster.galcat[\"theta\"], \"radians\", \"Mpc\", cluster.z, cosmo)\n", - "\n", - "cluster.make_radial_profile(\n", - " bins=da.make_bins(0.3, 3.0, 15, method=\"evenlog10width\"),\n", - " bin_units=\"Mpc\",\n", - " cosmo=cosmo,\n", - " include_empty_bins=False,\n", - " gal_ids_in_bins=True,\n", - ")\n", - "print(cluster.profile.colnames)" + "from clmm.utils import convert_units" ] }, { "cell_type": "code", "execution_count": null, - "id": "89050c01", + "id": "e2474c36-4b0c-4bab-86a4-021b4a4ef54d", "metadata": {}, "outputs": [], "source": [ - "fig = plt.figure(figsize=(6, 4))\n", - "ax = fig.add_axes((0, 0, 1, 1))\n", - "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", - "ax.errorbar(\n", - " cluster.profile[\"radius\"],\n", - " cluster.profile[\"gt\"],\n", - " cluster.profile[\"gt_err\"],\n", - " c=\"k\",\n", - " **errorbar_kwargs\n", + "def make_cluster_with_profile(\n", + " cluster_ra,\n", + " cluster_dec,\n", + " cluster_z,\n", + " galaxies,\n", + "):\n", + " obs_galaxies = galaxies[(galaxies[\"z\"] > (cluster_z + 0.1))]\n", + "\n", + " # Area cut: the query is made for the whole field and this can simplify the processing.\n", + " select = obs_galaxies[\"ra\"] < cluster_ra + 12.0 / 60.0 / np.cos(\n", + " cluster_dec / 180.0 * np.pi\n", + " )\n", + " select &= obs_galaxies[\"ra\"] > cluster_ra - 12.0 / 60.0 / np.cos(\n", + " cluster_dec / 180.0 * np.pi\n", + " )\n", + " select &= obs_galaxies[\"dec\"] < cluster_dec + 12.0 / 60.0\n", + " select &= obs_galaxies[\"dec\"] > cluster_dec - 12.0 / 60.0\n", + " obs_galaxies = obs_galaxies[select]\n", + "\n", + " obs_galaxies[\"id\"] = np.arange(len(obs_galaxies))\n", + "\n", + " # Put galaxy values on arrays.\n", + " gal_ra = obs_galaxies[\"ra\"] # Galaxies Ra in deg\n", + " gal_dec = obs_galaxies[\"dec\"] # Galaxies Dec in deg\n", + " gal_e1 = obs_galaxies[\"e1\"] # Galaxies elipticipy 1\n", + " gal_e2 = obs_galaxies[\"e2\"] # Galaxies elipticipy 2\n", + " gal_z = obs_galaxies[\"z\"] # Galaxies observed redshift\n", + " gal_id = obs_galaxies[\"id\"] # Galaxies ID\n", + "\n", + " # Create a GCData with the galaxies.\n", + " galaxies = clmm.GCData(\n", + " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"],\n", + " )\n", + "\n", + " # Create a GalaxyCluster.\n", + " cluster = clmm.GalaxyCluster(\n", + " \"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies\n", + " )\n", + "\n", + " # Convert elipticities into shears for the members.\n", + " cluster.compute_tangential_and_cross_components()\n", + " print(cluster.galcat.colnames)\n", + "\n", + " # Measure profile and add profile table to the cluster.\n", + " seps = convert_units(cluster.galcat[\"theta\"], \"radians\", \"Mpc\", cluster.z, cosmo)\n", + "\n", + " cluster.make_radial_profile(\n", + " bins=da.make_bins(0.3, 3.0, 15, method=\"evenlog10width\"),\n", + " bin_units=\"Mpc\",\n", + " cosmo=cosmo,\n", + " include_empty_bins=False,\n", + " gal_ids_in_bins=True,\n", + " )\n", + " print(cluster.profile.colnames)\n", + " return cluster\n", + "\n", + "\n", + "# We consider MaxBCG J140.53188+03.76632\n", + "# cluster_z = 0.2701 # Cluster redshift\n", + "# cluster_ra = 140.54565\n", + "# cluster_dec = 3.77820\n", + "\n", + "clusters = {}\n", + "clusters[\"HWL16a-026\"] = make_cluster_with_profile(\n", + " cluster_z=0.424, cluster_ra=130.5895, cluster_dec=1.6473, galaxies=data_2\n", ")\n", - "ax.set_xlabel(\"r [Mpc]\", fontsize=10)\n", - "ax.set_ylabel(r\"$g_t$\", fontsize=10)\n", - "ax.grid(lw=0.3)\n", - "ax.minorticks_on()\n", - "ax.grid(which=\"minor\", lw=0.1)\n", - "ax.set_xscale('log')\n", - "plt.show()" + "clusters[\"HWL16a-034\"] = make_cluster_with_profile(\n", + " cluster_z=0.315, cluster_ra=139.0387, cluster_dec=-0.3966, galaxies=data_2\n", + ")" ] }, { @@ -477,17 +472,25 @@ "z_inf = 1000\n", "concentration = 4.0\n", "\n", - "bs_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo))\n", - "bs2_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo) ** 2)\n", + "for cluster in clusters.values():\n", "\n", + " bs_mean = np.mean(\n", + " clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster.z, z_inf, cosmo)\n", + " )\n", + " bs2_mean = np.mean(\n", + " clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster.z, z_inf, cosmo) ** 2\n", + " )\n", "\n", - "def predict_reduced_tangential_shear_redshift_distribution(profile, logm):\n", + " cluster.z_src_betas = (bs_mean, bs2_mean)\n", + "\n", + "\n", + "def predict_reduced_tangential_shear_redshift_distribution(cluster, logm):\n", " gt = clmm.compute_reduced_tangential_shear(\n", - " r_proj=profile[\"radius\"], # Radial component of the profile\n", + " r_proj=cluster.data_for_fit[\"radius\"], # Radial component of the profile\n", " mdelta=10**logm, # Mass of the cluster [M_sun]\n", " cdelta=concentration, # Concentration of the cluster\n", - " z_cluster=cluster_z, # Redshift of the cluster\n", - " z_src=(bs_mean, bs2_mean), # tuple of (bs_mean, bs2_mean)\n", + " z_cluster=cluster.z, # Redshift of the cluster\n", + " z_src=cluster.z_src_betas, # tuple of (bs_mean, bs2_mean)\n", " z_src_info=\"beta\",\n", " approx=\"order1\",\n", " cosmo=cosmo,\n", @@ -499,12 +502,14 @@ "\n", "\n", "# Model using individual redshift and radial information, to compute the averaged shear in each radial bin, based on the galaxies actually present in that bin.\n", - "cluster.galcat[\"theta_mpc\"] = convert_units(\n", - " cluster.galcat[\"theta\"], \"radians\", \"mpc\", cluster.z, cosmo\n", - ")\n", + "\n", + "for cluster in clusters.values():\n", + " cluster.galcat[\"theta_mpc\"] = convert_units(\n", + " cluster.galcat[\"theta\"], \"radians\", \"mpc\", cluster.z, cosmo\n", + " )\n", "\n", "\n", - "def predict_reduced_tangential_shear_individual_redshift(profile, logm):\n", + "def predict_reduced_tangential_shear_individual_redshift(cluster, logm):\n", " return np.array(\n", " [\n", " np.mean(\n", @@ -513,7 +518,7 @@ " r_proj=cluster.galcat[radial_bin[\"gal_id\"]][\"theta_mpc\"],\n", " mdelta=10**logm, # Mass of the cluster [M_sun]\n", " cdelta=concentration, # Concentration of the cluster\n", - " z_cluster=cluster_z, # Redshift of the cluster\n", + " z_cluster=cluster.z, # Redshift of the cluster\n", " # Redshift value of each source galaxy inside the radial bin\n", " z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n", " cosmo=cosmo,\n", @@ -522,7 +527,7 @@ " halo_profile_model=\"nfw\",\n", " )\n", " )\n", - " for radial_bin in profile\n", + " for radial_bin in cluster.data_for_fit\n", " ]\n", " )" ] @@ -536,18 +541,20 @@ "source": [ "# Mass fitting\n", "\n", - "mask_for_fit = cluster.profile[\"n_src\"] > 2\n", - "data_for_fit = cluster.profile[mask_for_fit]\n", "\n", "from clmm.support.sampler import fitters\n", "\n", "\n", - "def fit_mass(predict_function):\n", + "def fit_mass(predict_function, cluster):\n", + "\n", + " cluster.mask_for_fit = cluster.profile[\"n_src\"] > 2\n", + " cluster.data_for_fit = cluster.profile[cluster.mask_for_fit]\n", + "\n", " popt, pcov = fitters[\"curve_fit\"](\n", " predict_function,\n", - " data_for_fit,\n", - " data_for_fit[\"gt\"],\n", - " data_for_fit[\"gt_err\"],\n", + " cluster,\n", + " cluster.data_for_fit[\"gt\"],\n", + " cluster.data_for_fit[\"gt_err\"],\n", " bounds=[10.0, 17.0],\n", " )\n", " logm, logm_err = popt[0], np.sqrt(pcov[0][0])\n", @@ -559,28 +566,6 @@ " }" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "9e6a4e04", - "metadata": {}, - "outputs": [], - "source": [ - "# In the paper, the measured mass is 44.3 {+ 30.3} {- 19.9} * 10^14 Msun (M200c,WL).\n", - "# For convenience, we consider a mean value for the errorbar.\n", - "# We build a dictionary based on that result.\n", - "m_paper = 44.3e14\n", - "m_err_paper = 25.1e14\n", - "logm_paper = np.log10(m_paper)\n", - "logm_err_paper = m_err_paper / (10**logm_paper) / np.log(10)\n", - "paper_value = {\n", - " \"logm\": logm_paper,\n", - " \"logm_err\": logm_err_paper,\n", - " \"m\": 10**logm_paper,\n", - " \"m_err\": (10**logm_paper) * logm_err_paper * np.log(10),\n", - "}" - ] - }, { "cell_type": "code", "execution_count": null, @@ -589,23 +574,13 @@ "outputs": [], "source": [ "%%time\n", - "fit_redshift_distribution = fit_mass(predict_reduced_tangential_shear_redshift_distribution)\n", - "fit_individual_redshift = fit_mass(predict_reduced_tangential_shear_individual_redshift)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "79e35734", - "metadata": {}, - "outputs": [], - "source": [ - "print(\n", - " f'Best fit mass for N(z) model = {fit_redshift_distribution[\"m\"]:.3e} +/- {fit_redshift_distribution[\"m_err\"]:.3e} Msun'\n", - ")\n", - "print(\n", - " f'Best fit mass for individual redshift and radius = {fit_individual_redshift[\"m\"]:.3e} +/- {fit_individual_redshift[\"m_err\"]:.3e} Msun'\n", - ")" + "for cluster in clusters.values():\n", + " cluster.fit_redshift_distribution = fit_mass(\n", + " predict_reduced_tangential_shear_redshift_distribution, cluster\n", + " )\n", + " cluster.fit_individual_redshift = fit_mass(\n", + " predict_reduced_tangential_shear_individual_redshift, cluster\n", + " )" ] }, { @@ -616,46 +591,46 @@ "outputs": [], "source": [ "# Visualization of the results.\n", - "def get_predicted_shear(predict_function, fit_values):\n", - " gt_est = predict_function(data_for_fit, fit_values[\"logm\"])\n", + "def get_predicted_shear(cluster, predict_function, fit_values):\n", + " gt_est = predict_function(cluster, fit_values[\"logm\"])\n", " gt_est_err = [\n", - " predict_function(data_for_fit, fit_values[\"logm\"] + i * fit_values[\"logm_err\"])\n", + " predict_function(cluster, fit_values[\"logm\"] + i * fit_values[\"logm_err\"])\n", " for i in (-3, 3)\n", " ]\n", " return gt_est, gt_est_err\n", "\n", "\n", - "gt_redshift_distribution, gt_err_redshift_distribution = get_predicted_shear(\n", - " predict_reduced_tangential_shear_redshift_distribution, fit_redshift_distribution\n", - ")\n", - "gt_individual_redshift, gt_err_individual_redshift = get_predicted_shear(\n", - " predict_reduced_tangential_shear_individual_redshift, fit_individual_redshift\n", - ")\n", - "\n", - "gt_paper1, gt_err_paper1 = get_predicted_shear(\n", - " predict_reduced_tangential_shear_redshift_distribution, paper_value\n", - ")\n", - "gt_paper2, gt_err_paper2 = get_predicted_shear(\n", - " predict_reduced_tangential_shear_individual_redshift, paper_value\n", - ")" + "for cluster in clusters.values():\n", + " cluster.gt_redshift_distribution, cluster.gt_err_redshift_distribution = (\n", + " get_predicted_shear(\n", + " cluster,\n", + " predict_reduced_tangential_shear_redshift_distribution,\n", + " cluster.fit_redshift_distribution,\n", + " )\n", + " )\n", + " cluster.gt_individual_redshift, cluster.gt_err_individual_redshift = (\n", + " get_predicted_shear(\n", + " cluster,\n", + " predict_reduced_tangential_shear_individual_redshift,\n", + " cluster.fit_individual_redshift,\n", + " )\n", + " )\n", + "# gt_paper1, gt_err_paper1 = get_predicted_shear(\n", + "# predict_reduced_tangential_shear_redshift_distribution, paper_value\n", + "# )\n", + "# gt_paper2, gt_err_paper2 = get_predicted_shear(\n", + "# predict_reduced_tangential_shear_individual_redshift, paper_value\n", + "# )" ] }, { "cell_type": "code", "execution_count": null, - "id": "ec85a1b1", + "id": "60c39758-3f56-4eca-b6df-8fc9d67db42e", "metadata": {}, "outputs": [], "source": [ - "chi2_redshift_distribution_dof = np.sum(\n", - " (gt_redshift_distribution - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", - ") / (len(data_for_fit) - 1)\n", - "chi2_individual_redshift_dof = np.sum(\n", - " (gt_individual_redshift - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", - ") / (len(data_for_fit) - 1)\n", - "\n", - "print(f\"Reduced chi2 (N(z) model) = {chi2_redshift_distribution_dof}\")\n", - "print(f\"Reduced chi2 (individual (R,z) model) = {chi2_individual_redshift_dof}\")" + "from paper_formating import prep_plot" ] }, { @@ -667,111 +642,166 @@ "source": [ "from matplotlib.ticker import MultipleLocator\n", "\n", - "fig = plt.figure(figsize=(6, 6))\n", - "gt_ax = fig.add_axes((0.25, 0.42, 0.7, 0.55))\n", - "gt_ax.errorbar(\n", - " data_for_fit[\"radius\"], data_for_fit[\"gt\"], data_for_fit[\"gt_err\"], c=\"k\", **errorbar_kwargs\n", - ")\n", "\n", - "# Points in grey have not been used for the fit.\n", - "gt_ax.errorbar(\n", - " cluster.profile[\"radius\"][~mask_for_fit],\n", - " cluster.profile[\"gt\"][~mask_for_fit],\n", - " cluster.profile[\"gt_err\"][~mask_for_fit],\n", - " c=\"grey\",\n", - " **errorbar_kwargs,\n", - ")\n", + "def make_paper_plot(cluster):\n", + " fig, axes = prep_plot(\n", + " subplots=3,\n", + " figsize=(9, 9),\n", + " subplots_kwargs=dict(sharex=True, height_ratios=[1, 3, 1]),\n", + " )\n", "\n", - "pow10 = 15\n", - "mlabel = lambda name, fits: (\n", - " rf\"$M_{{fit}}^{{{name}}} = \"\n", - " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", - " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", - " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", - ")\n", + " gx_ax, gt_ax, res_ax = axes\n", "\n", - "# The model for the 1st method.\n", - "gt_ax.loglog(\n", - " data_for_fit[\"radius\"],\n", - " gt_redshift_distribution,\n", - " \"-C1\",\n", - " label=mlabel(\"N(z)\", fit_redshift_distribution),\n", - " lw=0.5,\n", - ")\n", - "gt_ax.fill_between(\n", - " data_for_fit[\"radius\"], *gt_err_redshift_distribution, lw=0, color=\"C1\", alpha=0.2\n", - ")\n", + " errorbar_kwargs = dict(\n", + " linestyle=\"\", marker=\"o\", markersize=0.8, elinewidth=0.5, capthick=0.5\n", + " )\n", "\n", - "# The model for the 2nd method.\n", - "gt_ax.loglog(\n", - " data_for_fit[\"radius\"],\n", - " gt_individual_redshift,\n", - " \"-C2\",\n", - " label=mlabel(\"z,R\", fit_individual_redshift),\n", - " lw=0.5,\n", - ")\n", - "gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_individual_redshift, lw=0, color=\"C2\", alpha=0.2)\n", - "\n", - "# The value in the reference paper.\n", - "#gt_ax.loglog(\n", - "# data_for_fit[\"radius\"], gt_paper1, \"-C3\", label=mlabel(\"paper; N(z)\", paper_value), lw=0.5\n", - "#)\n", - "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper1, lw=0, color=\"C3\", alpha=0.2)\n", - "#\n", - "#gt_ax.loglog(\n", - "# data_for_fit[\"radius\"], gt_paper2, \"-C4\", label=mlabel(\"paper; Z,R\", paper_value), lw=0.5\n", - "#)\n", - "#gt_ax.fill_between(data_for_fit[\"radius\"], *gt_err_paper2, lw=0, color=\"C4\", alpha=0.2)\n", - "\n", - "gt_ax.set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", - "gt_ax.legend(fontsize=\"x-large\")\n", - "gt_ax.set_xticklabels([])\n", - "gt_ax.tick_params(\"x\", labelsize=\"x-large\")\n", - "gt_ax.tick_params(\"y\", labelsize=\"x-large\")\n", - "\n", - "\n", - "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", - "errorbar_kwargs2[\"markersize\"] = 5\n", - "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", - "res_ax = fig.add_axes((0.25, 0.2, 0.7, 0.2))\n", - "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.25\n", - "\n", - "\n", - "res_ax.errorbar(\n", - " data_for_fit[\"radius\"],\n", - " data_for_fit[\"gt\"] / gt_redshift_distribution - 1,\n", - " yerr=data_for_fit[\"gt_err\"] / gt_redshift_distribution,\n", - " marker=\"s\",\n", - " c=\"C1\",\n", - " **errorbar_kwargs2,\n", - ")\n", - "#errorbar_kwargs2[\"markersize\"] = 3\n", - "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", - "\n", - "res_ax.errorbar(\n", - " data_for_fit[\"radius\"] * delta,\n", - " data_for_fit[\"gt\"] / gt_individual_redshift - 1,\n", - " yerr=data_for_fit[\"gt_err\"] / gt_individual_redshift,\n", - " marker=\"*\",\n", - " c=\"C2\",\n", - " **errorbar_kwargs2,\n", - ")\n", - "res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", + " gx_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.data_for_fit[\"gx\"],\n", + " cluster.data_for_fit[\"gx_err\"],\n", + " c=\"k\",\n", + " **errorbar_kwargs,\n", + " )\n", + "\n", + " # Points in grey have not been used for the fit.\n", + " gx_ax.errorbar(\n", + " cluster.profile[\"radius\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gx\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gx_err\"][~cluster.mask_for_fit],\n", + " c=\"grey\",\n", + " **errorbar_kwargs,\n", + " )\n", "\n", - "res_ax.set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", - "res_ax.set_xscale(\"log\")\n", + " gt_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.data_for_fit[\"gt\"],\n", + " cluster.data_for_fit[\"gt_err\"],\n", + " c=\"k\",\n", + " **errorbar_kwargs,\n", + " )\n", "\n", - "res_ax.set_ylim(-5, 5)\n", - "res_ax.yaxis.set_minor_locator(MultipleLocator(10))\n", + " # Points in grey have not been used for the fit.\n", + " gt_ax.errorbar(\n", + " cluster.profile[\"radius\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gt\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gt_err\"][~cluster.mask_for_fit],\n", + " c=\"grey\",\n", + " **errorbar_kwargs,\n", + " )\n", "\n", - "res_ax.tick_params(\"x\", labelsize=\"x-large\")\n", - "res_ax.tick_params(\"y\", labelsize=\"x-large\")\n", + " pow10 = 15\n", + " mlabel = lambda name, fits: (\n", + " rf\"$M_{{fit}}^{{{name}}} = \"\n", + " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", + " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", + " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", + " )\n", + "\n", + " # The model for the 1st method.\n", + " gt_ax.loglog(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.gt_redshift_distribution,\n", + " \"-C1\",\n", + " label=mlabel(\"N(z)\", cluster.fit_redshift_distribution),\n", + " lw=0.5,\n", + " )\n", + " gt_ax.fill_between(\n", + " cluster.data_for_fit[\"radius\"],\n", + " *cluster.gt_err_redshift_distribution,\n", + " lw=0,\n", + " color=\"C1\",\n", + " alpha=0.2,\n", + " )\n", "\n", - "for ax in (gt_ax, res_ax):\n", - " ax.grid(lw=0.3, ls=':')\n", - " ax.minorticks_on()\n", - " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", - "plt.show()" + " # The model for the 2nd method.\n", + " gt_ax.loglog(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.gt_individual_redshift,\n", + " \"-C2\",\n", + " label=mlabel(\"z,R\", cluster.fit_individual_redshift),\n", + " lw=0.5,\n", + " )\n", + " gt_ax.fill_between(\n", + " cluster.data_for_fit[\"radius\"],\n", + " *cluster.gt_err_individual_redshift,\n", + " lw=0,\n", + " color=\"C2\",\n", + " alpha=0.2,\n", + " )\n", + "\n", + " errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", + " errorbar_kwargs2[\"markersize\"] = 1\n", + " errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + " delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.25\n", + "\n", + " res_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.data_for_fit[\"gt\"] / cluster.gt_redshift_distribution - 1,\n", + " yerr=cluster.data_for_fit[\"gt_err\"] / cluster.gt_redshift_distribution,\n", + " marker=\"s\",\n", + " c=\"C1\",\n", + " **errorbar_kwargs2,\n", + " )\n", + " # errorbar_kwargs2[\"markersize\"] = 3\n", + " # errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "\n", + " res_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"] * delta,\n", + " cluster.data_for_fit[\"gt\"] / cluster.gt_individual_redshift - 1,\n", + " yerr=cluster.data_for_fit[\"gt_err\"] / cluster.gt_individual_redshift,\n", + " marker=\"*\",\n", + " c=\"C2\",\n", + " **errorbar_kwargs2,\n", + " )\n", + "\n", + " gx_ax.set_ylabel(r\"$g_x$\", fontsize=8)\n", + " gt_ax.set_ylabel(r\"$g_t$\", fontsize=8)\n", + " gt_ax.legend(fontsize=8)\n", + " res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=8)\n", + "\n", + " res_ax.set_ylabel(\"$g_t$\\n [rel. diff]\", fontsize=8)\n", + " res_ax.set_xscale(\"log\")\n", + "\n", + " res_ax.set_ylim(-1.05, 1.05)\n", + " res_ax.yaxis.set_minor_locator(MultipleLocator(10))\n", + "\n", + " for ax in axes:\n", + " ax.tick_params(\"x\", labelsize=8)\n", + " ax.tick_params(\"y\", labelsize=8)\n", + "\n", + "\n", + " for ax in axes:\n", + " ax.grid(lw=0.3, ls=\"-\")\n", + " ax.minorticks_on()\n", + " ax.grid(which=\"minor\", lw=0.1, ls=\":\")\n", + "\n", + " plt.tight_layout()\n", + " plt.subplots_adjust(hspace=0.0)\n", + " return fig, axes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1a0a57c4-a199-4d15-b21a-148e490977fc", + "metadata": {}, + "outputs": [], + "source": [ + "make_paper_plot(clusters[\"HWL16a-026\"])\n", + "plt.savefig(\"hsc_1_new.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f51f28e9-07cd-4624-b491-db6eb2c01f78", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = make_paper_plot(clusters[\"HWL16a-034\"])\n", + "axes[1].set_ylim(1.01e-2, .25)\n", + "plt.savefig(\"hsc_2_new.png\")" ] }, { @@ -806,9 +836,9 @@ ], "metadata": { "kernelspec": { - "display_name": "conda-clmmenv", + "display_name": "clmm", "language": "python", - "name": "conda-clmmenv" + "name": "clmm" }, "language_info": { "codemirror_mode": { @@ -820,7 +850,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.9" + "version": "3.11.8" } }, "nbformat": 4, diff --git a/examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb b/examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb index cf18a85a5..52dc496ea 100644 --- a/examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb +++ b/examples/Paper_v2.0/Example5_Fit_Halo_mass_to_DES_data_paper2.ipynb @@ -48,11 +48,12 @@ "metadata": {}, "outputs": [], "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "from astropy.table import Table\n", "import pickle as pkl\n", - "from pathlib import Path" + "from pathlib import Path\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from astropy.table import Table" ] }, { @@ -125,21 +126,14 @@ { "cell_type": "code", "execution_count": null, - "id": "a04d966e-e00b-4e07-94d9-c3cc13b6183c", + "id": "42f0cb43-e6c8-4c1a-984e-8f0449a5e490", "metadata": {}, "outputs": [], "source": [ "%%time\n", "# Assume the downloaded catalog is at this path:\n", - "#filename = \"ACOS520_DES.csv\"\n", - "#filename = \"ID1.csv\"\n", - "filename = \"ID10.csv\"\n", - "catalog = filename.replace(\".csv\", \".pkl\")\n", - "if not Path(catalog).is_file():\n", - " data_0 = Table.read(filename, format=\"ascii.csv\")\n", - " pkl.dump(data_0, open(catalog, \"wb\"))\n", - "else:\n", - " data_0 = pkl.load(open(catalog, \"rb\"))" + "filename = \"/sps/lsst/groups/clusters/CLMM/DES_data_example/ID1+10.csv\"\n", + "data_0 = Table.read(filename, format=\"ascii.csv\")" ] }, { @@ -160,7 +154,8 @@ "outputs": [], "source": [ "from astropy.table import unique\n", - "data_0 = unique(data_0, [\"ra\",\"dec\"])" + "\n", + "data_0 = unique(data_0, [\"ra\", \"dec\"])" ] }, { @@ -240,38 +235,6 @@ "### Basic visualization" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "6d599e88-8a45-4b44-a46a-f2c00d86bd7b", - "metadata": {}, - "outputs": [], - "source": [ - "def make_plots(catalog_in):\n", - " # Scatter plot\n", - " plt.figure()\n", - " plt.scatter(catalog_in[\"ra\"], catalog_in[\"dec\"], c=catalog_in[\"z\"], s=1.0, alpha=1)\n", - " cb = plt.colorbar()\n", - " plt.xlabel(\"ra\")\n", - " plt.ylabel(\"dec\")\n", - " cb.ax.set_title(\"z\")\n", - "\n", - " # Histogram\n", - " plt.figure()\n", - " plt.hist(catalog_in[\"z\"], bins=20)\n", - " plt.xlabel(\"z\")\n", - " plt.ylabel(\"count\")\n", - "\n", - " # Relation\n", - " plt.figure()\n", - " plt.scatter(catalog_in[\"e1\"], catalog_in[\"e2\"], s=1.0, alpha=0.2)\n", - " plt.xlabel(\"e1\")\n", - " plt.ylabel(\"e2\")\n", - "\n", - "\n", - "make_plots(obs_galaxies)" - ] - }, { "cell_type": "markdown", "id": "e71d6d00-f18e-4e27-969b-f02f13032713", @@ -293,107 +256,109 @@ { "cell_type": "code", "execution_count": null, - "id": "1c4b465d-9e1c-4b63-a43e-eeea6139b462", + "id": "3514950f-eddf-4acc-8cbf-bd7adda54901", "metadata": {}, "outputs": [], "source": [ - "from clmm import Cosmology\n", "import clmm\n", + "from clmm import Cosmology\n", "\n", - "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)\n", - "\n", - "# We consider RMJ051637.4-543001.6 (ACO S520)\n", - "#cluster_z = 0.30416065 # Cluster redshift\n", - "#cluster_ra = 79.155704 # Cluster Ra in deg\n", - "#cluster_dec = -54.500456 # Cluster Dec in deg\n", - "\n", - "# ID1\n", - "#cluster_z = 0.4298039972782135 # Cluster redshift\n", - "#cluster_ra = 43.564574 # Cluster Ra in deg\n", - "#cluster_dec = -58.95297 # Cluster Dec in deg\n", - "\n", - "# ID10\n", - "cluster_z = 0.3264264166355133 # Cluster redshift\n", - "cluster_ra = 323.800394 # Cluster Ra in deg\n", - "cluster_dec = -1.04959 # Cluster Dec in deg\n", - "\n", - "# Select background galaxies\n", - "obs_galaxies = obs_galaxies[(obs_galaxies[\"z\"] > (cluster_z + 0.1)) & (obs_galaxies[\"z\"] < 1.5)]\n", - "\n", - "obs_galaxies[\"id\"] = np.arange(len(obs_galaxies))\n", - "\n", - "# Put galaxy values on arrays\n", - "gal_ra = obs_galaxies[\"ra\"] # Galaxies Ra in deg\n", - "gal_dec = obs_galaxies[\"dec\"] # Galaxies Dec in deg\n", - "gal_e1 = obs_galaxies[\"e1\"] # Galaxies elipticipy 1\n", - "gal_e2 = obs_galaxies[\"e2\"] # Galaxies elipticipy 2\n", - "gal_z = obs_galaxies[\"z\"] # Galaxies observed redshift\n", - "gal_id = obs_galaxies[\"id\"] # Galaxies ID\n", - "\n", - "# Create a GCData with the galaxies.\n", - "galaxies = clmm.GCData(\n", - " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id], names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"]\n", - ")\n", - "\n", - "# Create a GalaxyCluster.\n", - "cluster = clmm.GalaxyCluster(\"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies)" + "cosmo = Cosmology(H0=70.0, Omega_dm0=0.27 - 0.045, Omega_b0=0.045, Omega_k0=0.0)" ] }, { - "cell_type": "markdown", - "id": "4d526e90-3c58-4f4e-ad9b-58c95464d72d", + "cell_type": "code", + "execution_count": null, + "id": "f9a8c215-c6d8-4544-a57f-973fd165e257", "metadata": {}, + "outputs": [], "source": [ - "### Measure the shear profile" + "import clmm.dataops as da\n", + "from clmm.utils import convert_units" ] }, { - "cell_type": "code", - "execution_count": null, - "id": "42573e50-85cc-4977-baf4-723711792c6b", + "cell_type": "markdown", + "id": "4d526e90-3c58-4f4e-ad9b-58c95464d72d", "metadata": {}, - "outputs": [], "source": [ - "import clmm.dataops as da\n", - "\n", - "# Convert ellipticities into shears for the members.\n", - "cluster.compute_tangential_and_cross_components()\n", - "print(cluster.galcat.colnames)\n", - "\n", - "# Measure profile and add profile table to the cluster.\n", - "cluster.make_radial_profile(\n", - " bins=da.make_bins(0.2, 4.0, 6, method=\"evenlog10width\"),\n", - " bin_units=\"Mpc\",\n", - " cosmo=cosmo,\n", - " include_empty_bins=False,\n", - " gal_ids_in_bins=True,\n", - ")\n", - "print(cluster.profile.colnames)" + "### Measure the shear profile" ] }, { "cell_type": "code", "execution_count": null, - "id": "db4ebc79-67e8-43da-ad89-ca96666f541d", + "id": "1591b322-3417-49d9-81db-7d868cee6f88", "metadata": {}, "outputs": [], "source": [ - "fig, ax = plt.subplots(figsize=(7, 5), ncols=1, nrows=1)\n", - "errorbar_kwargs = dict(linestyle=\"\", marker=\"o\", markersize=3, elinewidth=0.5, capthick=0.5)\n", - "ax.errorbar(\n", - " cluster.profile[\"radius\"],\n", - " cluster.profile[\"gt\"],\n", - " cluster.profile[\"gt_err\"],\n", - " c=\"k\",\n", - " **errorbar_kwargs\n", + "def make_cluster_with_profile(\n", + " cluster_ra,\n", + " cluster_dec,\n", + " cluster_z,\n", + " galaxies,\n", + "):\n", + "\n", + " # Select background galaxies\n", + " obs_galaxies = galaxies[(galaxies[\"z\"] > (cluster_z + 0.1)) & (galaxies[\"z\"] < 1.5)]\n", + "\n", + " obs_galaxies[\"id\"] = np.arange(len(obs_galaxies))\n", + "\n", + " # Put galaxy values on arrays\n", + " gal_ra = obs_galaxies[\"ra\"] # Galaxies Ra in deg\n", + " gal_dec = obs_galaxies[\"dec\"] # Galaxies Dec in deg\n", + " gal_e1 = obs_galaxies[\"e1\"] # Galaxies elipticipy 1\n", + " gal_e2 = obs_galaxies[\"e2\"] # Galaxies elipticipy 2\n", + " gal_z = obs_galaxies[\"z\"] # Galaxies observed redshift\n", + " gal_id = obs_galaxies[\"id\"] # Galaxies ID\n", + "\n", + " # Create a GCData with the galaxies.\n", + " galaxies = clmm.GCData(\n", + " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"],\n", + " )\n", + "\n", + " # Create a GalaxyCluster.\n", + " cluster = clmm.GalaxyCluster(\n", + " \"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies\n", + " )\n", + "\n", + " # Convert ellipticities into shears for the members.\n", + " cluster.compute_tangential_and_cross_components()\n", + " print(cluster.galcat.colnames)\n", + "\n", + " # Measure profile and add profile table to the cluster.\n", + " cluster.make_radial_profile(\n", + " bins=da.make_bins(0.2, 4.0, 6, method=\"evenlog10width\"),\n", + " bin_units=\"Mpc\",\n", + " cosmo=cosmo,\n", + " include_empty_bins=True,\n", + " gal_ids_in_bins=True,\n", + " )\n", + " print(cluster.profile.colnames)\n", + "\n", + " return cluster\n", + "\n", + "\n", + "# We consider RMJ051637.4-543001.6 (ACO S520)\n", + "# cluster_z = 0.30416065 # Cluster redshift\n", + "# cluster_ra = 79.155704 # Cluster Ra in deg\n", + "# cluster_dec = -54.500456 # Cluster Dec in deg\n", + "\n", + "clusters = {}\n", + "\n", + "clusters[\"ID1\"] = make_cluster_with_profile(\n", + " cluster_z=0.4298039972782135,\n", + " cluster_ra=43.564574,\n", + " cluster_dec=-58.95297,\n", + " galaxies=obs_galaxies,\n", ")\n", - "ax.set_xlabel(\"R [Mpc]\", fontsize=10)\n", - "ax.set_ylabel(r\"$g_t$\", fontsize=10)\n", - "ax.set_xscale(\"log\")\n", - "ax.grid(lw=0.3)\n", - "ax.minorticks_on()\n", - "ax.grid(which=\"minor\", lw=0.1)\n", - "plt.show()" + "clusters[\"ID10\"] = make_cluster_with_profile(\n", + " cluster_z=0.3264264166355133,\n", + " cluster_ra=323.800394,\n", + " cluster_dec=-1.04959,\n", + " galaxies=obs_galaxies,\n", + ")" ] }, { @@ -411,23 +376,35 @@ "metadata": {}, "outputs": [], "source": [ - "from clmm.utils import convert_units\n", + "# Theoretical predictions\n", "\n", "# Model relying on the overall redshift distribution of the sources (WtG III Applegate et al. 2014).\n", + "# Note the concentration of MaxBCG J140.53188+03.76632 was not reported.\n", + "# The value from the stacked sample in the paper is ~7.\n", + "# For the mass scale, a typical c-M relation (e.g. Child et al. 2018) would give c~3 though.\n", + "# And we have not considered a c-M relation in the fitting.\n", "z_inf = 1000\n", "concentration = 4.0\n", "\n", - "bs_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo))\n", - "bs2_mean = np.mean(clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster_z, z_inf, cosmo) ** 2)\n", + "for cluster in clusters.values():\n", + "\n", + " bs_mean = np.mean(\n", + " clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster.z, z_inf, cosmo)\n", + " )\n", + " bs2_mean = np.mean(\n", + " clmm.utils.compute_beta_s(cluster.galcat[\"z\"], cluster.z, z_inf, cosmo) ** 2\n", + " )\n", + "\n", + " cluster.z_src_betas = (bs_mean, bs2_mean)\n", "\n", "\n", - "def predict_reduced_tangential_shear_redshift_distribution(profile, logm):\n", + "def predict_reduced_tangential_shear_redshift_distribution(cluster, logm):\n", " gt = clmm.compute_reduced_tangential_shear(\n", - " r_proj=profile[\"radius\"], # Radial component of the profile\n", + " r_proj=cluster.data_for_fit[\"radius\"], # Radial component of the profile\n", " mdelta=10**logm, # Mass of the cluster [M_sun]\n", " cdelta=concentration, # Concentration of the cluster\n", - " z_cluster=cluster_z, # Redshift of the cluster\n", - " z_src=(bs_mean, bs2_mean), # tuple of (bs_mean, bs2_mean)\n", + " z_cluster=cluster.z, # Redshift of the cluster\n", + " z_src=cluster.z_src_betas, # tuple of (bs_mean, bs2_mean)\n", " z_src_info=\"beta\",\n", " approx=\"order1\",\n", " cosmo=cosmo,\n", @@ -439,12 +416,14 @@ "\n", "\n", "# Model using individual redshift and radial information, to compute the averaged shear in each radial bin, based on the galaxies actually present in that bin.\n", - "cluster.galcat[\"theta_mpc\"] = convert_units(\n", - " cluster.galcat[\"theta\"], \"radians\", \"mpc\", cluster.z, cosmo\n", - ")\n", + "\n", + "for cluster in clusters.values():\n", + " cluster.galcat[\"theta_mpc\"] = convert_units(\n", + " cluster.galcat[\"theta\"], \"radians\", \"mpc\", cluster.z, cosmo\n", + " )\n", "\n", "\n", - "def predict_reduced_tangential_shear_individual_redshift(profile, logm):\n", + "def predict_reduced_tangential_shear_individual_redshift(cluster, logm):\n", " return np.array(\n", " [\n", " np.mean(\n", @@ -453,7 +432,7 @@ " r_proj=cluster.galcat[radial_bin[\"gal_id\"]][\"theta_mpc\"],\n", " mdelta=10**logm, # Mass of the cluster [M_sun]\n", " cdelta=concentration, # Concentration of the cluster\n", - " z_cluster=cluster_z, # Redshift of the cluster\n", + " z_cluster=cluster.z, # Redshift of the cluster\n", " # Redshift value of each source galaxy inside the radial bin\n", " z_src=cluster.galcat[radial_bin[\"gal_id\"]][\"z\"],\n", " cosmo=cosmo,\n", @@ -462,7 +441,7 @@ " halo_profile_model=\"nfw\",\n", " )\n", " )\n", - " for radial_bin in profile\n", + " for radial_bin in cluster.data_for_fit\n", " ]\n", " )" ] @@ -482,18 +461,22 @@ "metadata": {}, "outputs": [], "source": [ - "mask_for_fit = cluster.profile[\"n_src\"] > 2\n", - "data_for_fit = cluster.profile[mask_for_fit]\n", + "# Mass fitting\n", + "\n", "\n", "from clmm.support.sampler import fitters\n", "\n", "\n", - "def fit_mass(predict_function):\n", + "def fit_mass(predict_function, cluster):\n", + "\n", + " cluster.mask_for_fit = cluster.profile[\"n_src\"] > 2\n", + " cluster.data_for_fit = cluster.profile[cluster.mask_for_fit]\n", + "\n", " popt, pcov = fitters[\"curve_fit\"](\n", " predict_function,\n", - " data_for_fit,\n", - " data_for_fit[\"gt\"],\n", - " data_for_fit[\"gt_err\"],\n", + " cluster,\n", + " cluster.data_for_fit[\"gt\"],\n", + " cluster.data_for_fit[\"gt_err\"],\n", " bounds=[10.0, 17.0],\n", " )\n", " logm, logm_err = popt[0], np.sqrt(pcov[0][0])\n", @@ -513,25 +496,13 @@ "outputs": [], "source": [ "%%time\n", - "fit_redshift_distribution = fit_mass(predict_reduced_tangential_shear_redshift_distribution)\n", - "fit_individual_redshift = fit_mass(predict_reduced_tangential_shear_individual_redshift)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2f5b8ac6-7294-4e85-97ac-28b62512be60", - "metadata": {}, - "outputs": [], - "source": [ - "print(\n", - " \"Best fit mass for N(z) model =\"\n", - " f' {fit_redshift_distribution[\"m\"]:.3e} +/- {fit_redshift_distribution[\"m_err\"]:.3e} Msun'\n", - ")\n", - "print(\n", - " \"Best fit mass for individual redshift and radius =\"\n", - " f' {fit_individual_redshift[\"m\"]:.3e} +/- {fit_individual_redshift[\"m_err\"]:.3e} Msun'\n", - ")" + "for cluster in clusters.values():\n", + " cluster.fit_redshift_distribution = fit_mass(\n", + " predict_reduced_tangential_shear_redshift_distribution, cluster\n", + " )\n", + " cluster.fit_individual_redshift = fit_mass(\n", + " predict_reduced_tangential_shear_individual_redshift, cluster\n", + " )" ] }, { @@ -549,39 +520,41 @@ "metadata": {}, "outputs": [], "source": [ - "def get_predicted_shear(predict_function, fit_values):\n", - " gt_est = predict_function(data_for_fit, fit_values[\"logm\"])\n", + "# Visualization of the results.\n", + "def get_predicted_shear(cluster, predict_function, fit_values):\n", + " gt_est = predict_function(cluster, fit_values[\"logm\"])\n", " gt_est_err = [\n", - " predict_function(data_for_fit, fit_values[\"logm\"] + i * fit_values[\"logm_err\"])\n", + " predict_function(cluster, fit_values[\"logm\"] + i * fit_values[\"logm_err\"])\n", " for i in (-3, 3)\n", " ]\n", " return gt_est, gt_est_err\n", "\n", "\n", - "gt_redshift_distribution, gt_err_redshift_distribution = get_predicted_shear(\n", - " predict_reduced_tangential_shear_redshift_distribution, fit_redshift_distribution\n", - ")\n", - "gt_individual_redshift, gt_err_individual_redshift = get_predicted_shear(\n", - " predict_reduced_tangential_shear_individual_redshift, fit_individual_redshift\n", - ")" + "for cluster in clusters.values():\n", + " cluster.gt_redshift_distribution, cluster.gt_err_redshift_distribution = (\n", + " get_predicted_shear(\n", + " cluster,\n", + " predict_reduced_tangential_shear_redshift_distribution,\n", + " cluster.fit_redshift_distribution,\n", + " )\n", + " )\n", + " cluster.gt_individual_redshift, cluster.gt_err_individual_redshift = (\n", + " get_predicted_shear(\n", + " cluster,\n", + " predict_reduced_tangential_shear_individual_redshift,\n", + " cluster.fit_individual_redshift,\n", + " )\n", + " )" ] }, { "cell_type": "code", "execution_count": null, - "id": "2d0e95d3-bc69-4cf7-8ac0-1788b74fe5d6", + "id": "e8074a76-23ce-4244-a73d-483a91f5e450", "metadata": {}, "outputs": [], "source": [ - "chi2_redshift_distribution_dof = np.sum(\n", - " (gt_redshift_distribution - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", - ") / (len(data_for_fit) - 1)\n", - "chi2_individual_redshift_dof = np.sum(\n", - " (gt_individual_redshift - data_for_fit[\"gt\"]) ** 2 / (data_for_fit[\"gt_err\"]) ** 2\n", - ") / (len(data_for_fit) - 1)\n", - "\n", - "print(f\"Reduced chi2 (N(z) model) = {chi2_redshift_distribution_dof}\")\n", - "print(f\"Reduced chi2 (individual (R,z) model) = {chi2_individual_redshift_dof}\")" + "from paper_formating import prep_plot" ] }, { @@ -591,112 +564,167 @@ "metadata": {}, "outputs": [], "source": [ - "fig, axes = plt.subplots(\n", - " nrows=2,\n", - " ncols=1,\n", - " #figsize=(8, 6),\n", - " figsize=(8*0.8, 6*0.8),\n", - " gridspec_kw={\"height_ratios\": [3, 1], \"wspace\": 0.4, \"hspace\": 0.03},\n", - ")\n", + "from matplotlib.ticker import MultipleLocator\n", "\n", - "axes[0].errorbar(\n", - " data_for_fit[\"radius\"], data_for_fit[\"gt\"], data_for_fit[\"gt_err\"], c=\"k\", **errorbar_kwargs\n", - ")\n", "\n", - "# Points in grey have not been used for the fit.\n", - "axes[0].errorbar(\n", - " cluster.profile[\"radius\"][~mask_for_fit],\n", - " cluster.profile[\"gt\"][~mask_for_fit],\n", - " cluster.profile[\"gt_err\"][~mask_for_fit],\n", - " c=\"grey\",\n", - " **errorbar_kwargs,\n", - ")\n", + "def make_paper_plot(cluster):\n", + " fig, axes = prep_plot(\n", + " subplots=3,\n", + " figsize=(9, 9),\n", + " subplots_kwargs=dict(sharex=True, height_ratios=[1, 3, 1]),\n", + " )\n", "\n", - "pow10 = 15\n", - "mlabel = lambda name, fits: (\n", - " rf\"$M_{{fit}}^{{{name}}} = \"\n", - " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", - " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", - " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", - ")\n", + " gx_ax, gt_ax, res_ax = axes\n", "\n", - "# The model for the 1st method.\n", - "axes[0].loglog(\n", - " data_for_fit[\"radius\"],\n", - " gt_redshift_distribution,\n", - " \"-C1\",\n", - " label=mlabel(\"N(z)\", fit_redshift_distribution),\n", - " lw=0.5,\n", - ")\n", + " errorbar_kwargs = dict(\n", + " linestyle=\"\", marker=\"o\", markersize=0.8, elinewidth=0.5, capthick=0.5\n", + " )\n", "\n", - "axes[0].fill_between(\n", - " data_for_fit[\"radius\"], *gt_err_redshift_distribution, lw=0, color=\"C1\", alpha=0.2\n", - ")\n", + " gx_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.data_for_fit[\"gx\"],\n", + " cluster.data_for_fit[\"gx_err\"],\n", + " c=\"k\",\n", + " **errorbar_kwargs,\n", + " )\n", "\n", - "# The model for the 2nd method.\n", - "axes[0].loglog(\n", - " data_for_fit[\"radius\"],\n", - " gt_individual_redshift,\n", - " \"-C2\",\n", - " label=mlabel(\"z,R\", fit_individual_redshift),\n", - " lw=0.5,\n", - ")\n", - "axes[0].fill_between(\n", - " data_for_fit[\"radius\"], *gt_err_individual_redshift, lw=0, color=\"C2\", alpha=0.2\n", - ")\n", + " # Points in grey have not been used for the fit.\n", + " gx_ax.errorbar(\n", + " cluster.profile[\"radius\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gx\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gx_err\"][~cluster.mask_for_fit],\n", + " c=\"grey\",\n", + " **errorbar_kwargs,\n", + " )\n", "\n", + " gt_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.data_for_fit[\"gt\"],\n", + " cluster.data_for_fit[\"gt_err\"],\n", + " c=\"k\",\n", + " **errorbar_kwargs,\n", + " )\n", "\n", - "axes[0].set_ylabel(r\"$g_t$\", fontsize=\"x-large\")\n", - "axes[0].legend(fontsize=\"x-large\", loc=\"lower left\")\n", - "axes[0].set_xticklabels([])\n", - "axes[0].tick_params(\"x\", labelsize=\"x-large\")\n", - "axes[0].tick_params(\"y\", labelsize=\"x-large\")\n", - "axes[1].set_ylim(1.0e-3, 0.5)\n", + " # Points in grey have not been used for the fit.\n", + " gt_ax.errorbar(\n", + " cluster.profile[\"radius\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gt\"][~cluster.mask_for_fit],\n", + " cluster.profile[\"gt_err\"][~cluster.mask_for_fit],\n", + " c=\"grey\",\n", + " **errorbar_kwargs,\n", + " )\n", "\n", + " pow10 = 15\n", + " mlabel = lambda name, fits: (\n", + " rf\"$M_{{fit}}^{{{name}}} = \"\n", + " rf'{fits[\"m\"]/10**pow10:.2f}\\pm'\n", + " rf'{fits[\"m_err\"]/10**pow10:.2f}'\n", + " rf\"\\times 10^{{{pow10}}} M_\\odot$\"\n", + " )\n", "\n", - "errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", - "errorbar_kwargs2[\"markersize\"] = 5\n", - "errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + " # The model for the 1st method.\n", + " gt_ax.loglog(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.gt_redshift_distribution,\n", + " \"-C1\",\n", + " label=mlabel(\"N(z)\", cluster.fit_redshift_distribution),\n", + " lw=0.5,\n", + " )\n", + " gt_ax.fill_between(\n", + " cluster.data_for_fit[\"radius\"],\n", + " *cluster.gt_err_redshift_distribution,\n", + " lw=0,\n", + " color=\"C1\",\n", + " alpha=0.2,\n", + " )\n", "\n", - "delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.15\n", + " # The model for the 2nd method.\n", + " gt_ax.loglog(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.gt_individual_redshift,\n", + " \"-C2\",\n", + " label=mlabel(\"z,R\", cluster.fit_individual_redshift),\n", + " lw=0.5,\n", + " )\n", + " gt_ax.fill_between(\n", + " cluster.data_for_fit[\"radius\"],\n", + " *cluster.gt_err_individual_redshift,\n", + " lw=0,\n", + " color=\"C2\",\n", + " alpha=0.2,\n", + " )\n", "\n", + " errorbar_kwargs2 = {k: v for k, v in errorbar_kwargs.items() if \"marker\" not in k}\n", + " errorbar_kwargs2[\"markersize\"] = 1\n", + " errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + " delta = (cluster.profile[\"radius\"][1] / cluster.profile[\"radius\"][0]) ** 0.25\n", + "\n", + " res_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"],\n", + " cluster.data_for_fit[\"gt\"] / cluster.gt_redshift_distribution - 1,\n", + " yerr=cluster.data_for_fit[\"gt_err\"] / cluster.gt_redshift_distribution,\n", + " marker=\"s\",\n", + " c=\"C1\",\n", + " **errorbar_kwargs2,\n", + " )\n", + " # errorbar_kwargs2[\"markersize\"] = 3\n", + " # errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", + "\n", + " res_ax.errorbar(\n", + " cluster.data_for_fit[\"radius\"] * delta,\n", + " cluster.data_for_fit[\"gt\"] / cluster.gt_individual_redshift - 1,\n", + " yerr=cluster.data_for_fit[\"gt_err\"] / cluster.gt_individual_redshift,\n", + " marker=\"*\",\n", + " c=\"C2\",\n", + " **errorbar_kwargs2,\n", + " )\n", "\n", - "axes[1].errorbar(\n", - " data_for_fit[\"radius\"],\n", - " data_for_fit[\"gt\"] / gt_redshift_distribution - 1,\n", - " yerr=data_for_fit[\"gt_err\"] / gt_redshift_distribution,\n", - " marker=\"s\",\n", - " c=\"C1\",\n", - " **errorbar_kwargs2,\n", - ")\n", - "#errorbar_kwargs2[\"markersize\"] = 3\n", - "#errorbar_kwargs2[\"markeredgewidth\"] = 0.5\n", - "\n", - "axes[1].errorbar(\n", - " data_for_fit[\"radius\"] * delta,\n", - " data_for_fit[\"gt\"] / gt_individual_redshift - 1,\n", - " yerr=data_for_fit[\"gt_err\"] / gt_individual_redshift,\n", - " marker=\"*\",\n", - " c=\"C2\",\n", - " **errorbar_kwargs2,\n", - ")\n", - "axes[1].set_xlabel(r\"$R$ [Mpc]\", fontsize=\"x-large\")\n", + " gx_ax.set_ylabel(r\"$g_x$\", fontsize=8)\n", + " gt_ax.set_ylabel(r\"$g_t$\", fontsize=8)\n", + " gt_ax.legend(fontsize=8)\n", + " res_ax.set_xlabel(r\"$R$ [Mpc]\", fontsize=8)\n", + "\n", + " res_ax.set_ylabel(\"$g_t$\\n [rel. diff]\", fontsize=8)\n", + " res_ax.set_xscale(\"log\")\n", "\n", - "axes[1].set_ylabel(r\"$g_t^{data}/g_t^{mod.}-1$\", fontsize=\"x-large\")\n", - "axes[1].set_xscale(\"log\")\n", + " res_ax.set_ylim(-1.05, 1.05)\n", + " res_ax.yaxis.set_minor_locator(MultipleLocator(10))\n", "\n", - "axes[1].set_ylim(-5, 5)\n", + " for ax in axes:\n", + " ax.tick_params(\"x\", labelsize=8)\n", + " ax.tick_params(\"y\", labelsize=8)\n", "\n", - "axes[1].tick_params(\"x\", labelsize=\"x-large\")\n", - "axes[1].tick_params(\"y\", labelsize=\"x-large\")\n", "\n", - "for ax in axes:\n", - " ax.grid(lw=0.3, ls=':')\n", - " ax.minorticks_on()\n", - " ax.grid(which=\"minor\", lw=0.1, ls=':')\n", - "plt.show()\n", + " for ax in axes:\n", + " ax.grid(lw=0.3, ls=\"-\")\n", + " ax.minorticks_on()\n", + " ax.grid(which=\"minor\", lw=0.1, ls=\":\")\n", "\n", - "# Note since we made cuts on the catalog, the redshift distribution of the remaining sources might not be representative." + " plt.tight_layout()\n", + " plt.subplots_adjust(hspace=0.0)\n", + " return fig, axes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "64838792-6fb1-4c4e-8547-1c0cdda0c7aa", + "metadata": {}, + "outputs": [], + "source": [ + "make_paper_plot(clusters[\"ID1\"])\n", + "plt.savefig(\"des_ID1_new.png\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fbbdeea-fc25-465a-a0a7-852a93b8a7a7", + "metadata": {}, + "outputs": [], + "source": [ + "make_paper_plot(clusters[\"ID10\"])\n", + "plt.savefig(\"des_ID10_new.png\")" ] }, { @@ -717,9 +745,9 @@ ], "metadata": { "kernelspec": { - "display_name": "conda-clmmenv", + "display_name": "clmm", "language": "python", - "name": "conda-clmmenv" + "name": "clmm" }, "language_info": { "codemirror_mode": { @@ -731,7 +759,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.9" + "version": "3.11.8" } }, "nbformat": 4,