From e7164d27cecfe4bb0305f3a252f18bb34ae5f26a Mon Sep 17 00:00:00 2001 From: Celine Combet Date: Fri, 18 Oct 2024 17:53:39 +0200 Subject: [PATCH] Fix some minor problems in boost calculation and improve docstring (#646) * Make keyword non-optional, fix docstrings for boost models * update boost correction function names * Update version to 1.14.1 --------- Co-authored-by: m-aguena --- clmm/__init__.py | 2 +- clmm/utils/__init__.py | 4 +- clmm/utils/boost.py | 108 ++++++++++++++++++------------ examples/demo_boost_factors.ipynb | 68 ++++++++++++++----- tests/test_utils.py | 19 ++++-- 5 files changed, 131 insertions(+), 70 deletions(-) diff --git a/clmm/__init__.py b/clmm/__init__.py index 64fee9cf1..0b688448d 100644 --- a/clmm/__init__.py +++ b/clmm/__init__.py @@ -26,4 +26,4 @@ ) from . import support -__version__ = "1.14.0" +__version__ = "1.14.1" diff --git a/clmm/utils/__init__.py b/clmm/utils/__init__.py index bcd1806ac..a2c0854bb 100644 --- a/clmm/utils/__init__.py +++ b/clmm/utils/__init__.py @@ -13,8 +13,8 @@ from .boost import ( compute_nfw_boost, compute_powerlaw_boost, - correct_sigma_with_boost_values, - correct_sigma_with_boost_model, + correct_with_boost_values, + correct_with_boost_model, boost_models, ) diff --git a/clmm/utils/boost.py b/clmm/utils/boost.py index 0abaa8a1f..78e22c695 100644 --- a/clmm/utils/boost.py +++ b/clmm/utils/boost.py @@ -1,27 +1,42 @@ """General utility functions that are used in multiple modules""" + import numpy as np -def compute_nfw_boost(rvals, rscale=1000, boost0=0.1): - """Given a list of `rvals`, and optional `rscale` and `boost0`, return the corresponding - boost factor at each rval +def compute_nfw_boost(r_proj, r_scale, boost0=0.1): + r"""Computes the boost factor at radii r_proj using a parametric form + following that of the analytical NFW excess surface density. + + Setting :math:`x = r_{\rm proj}/r_{\rm scale}` + + .. math:: + B(x) = 1 + B_0\frac{1-F(x)}{(x)^2 -1} + + where + + .. math:: + F(x) = \begin{cases} \frac{\arctan \sqrt{x^2 -1 }}{\sqrt{x^2 -1 }} & \text{if $x>1$}, \\ + \text{1} & \text{if x = 1}, \\ + \frac{\text{arctanh} \sqrt{1-x^2 }}{\sqrt{1- x^2 }} & \text{if $x<1$}. + \end{cases} + Parameters ---------- - rvals : array_like + r_proj : array_like Radii - rscale : float, optional - scale radius for NFW in same units as rvals (default 2000 kpc) - boost0 : float, optional - Boost factor at each value of rvals + r_scale : float + Scale radius in same units as r_proj + boost0: float, optional + Boost factor normalisation Returns ------- array - Boost factor + Boost factor at each value of r_proj """ - r_norm = np.array(rvals) / rscale + r_norm = np.array(r_proj) / r_scale def _calc_finternal(r_norm): radicand = r_norm**2 - 1 @@ -34,33 +49,37 @@ def _calc_finternal(r_norm): / (2 * np.lib.scimath.sqrt(radicand)) ) - return np.nan_to_num(finternal, copy=False, nan=1.0).real + return np.nan_to_num((1 - finternal) / radicand, copy=False, nan=1.0 / 3.0).real + # return np.nan_to_num(finternal, copy=False, nan=1.0).real + + return 1.0 + boost0 * _calc_finternal(r_norm) + # return 1.0 + boost0 * (1 - _calc_finternal(r_norm)) / (r_norm**2 - 1) - return 1.0 + boost0 * (1 - _calc_finternal(r_norm)) / (r_norm**2 - 1) +def compute_powerlaw_boost(r_proj, r_scale, boost0=0.1, alpha=-1): + r"""Computes the boost factor at radii r_proj using a power-law parametric form -def compute_powerlaw_boost(rvals, rscale=1000, boost0=0.1, alpha=-1.0): - """Given a list of `rvals`, and optional `rscale` and `boost0`, and `alpha`, - return the corresponding boost factor at each `rval` + .. math:: + B(r_{\rm proj}) = 1 + B_0 \left(\frac{r_{\rm proj}}{r_{\rm scale}}\right)^\alpha Parameters ---------- - rvals : array_like + r_proj : array_like Radii - rscale : float, optional - Scale radius for NFW in same units as rvals (default 2000 kpc) + r_scale : float + Scale radius in same units as r_proj boost0 : float, optional - Boost factor at each value of rvals + Boost factor normalisation alpha : float, optional - Exponent from Melchior+16. Default: -1.0 + Exponent for the power-law parametrisation. NB: -1.0 (as in Melchior+16) Returns ------- array - Boost factor + Boost factor at each value of r_proj """ - r_norm = np.array(rvals) / rscale + r_norm = np.array(r_proj) / r_scale return 1.0 + boost0 * (r_norm) ** alpha @@ -71,51 +90,54 @@ def compute_powerlaw_boost(rvals, rscale=1000, boost0=0.1, alpha=-1.0): } -def correct_sigma_with_boost_values(sigma_vals, boost_factors): - """Given a list of boost values and sigma profile, compute corrected sigma +def correct_with_boost_values(profile_vals, boost_factors): + """Given a list of profile values (e.g., shear or DeltaSigma) and boost values, + computes corrected profile Parameters ---------- - sigma_vals : array_like - uncorrected sigma with cluster member dilution + profile_vals : array_like + Uncorrected profile values boost_factors : array_like - Boost values pre-computed + Boost values, pre-computed Returns ------- - sigma_corrected : numpy.ndarray - correted radial profile + profile_vals_corrected : numpy.ndarray + Corrected radial profile """ - sigma_corrected = np.array(sigma_vals) / np.array(boost_factors) - return sigma_corrected + profile_vals_corrected = np.array(profile_vals) / np.array(boost_factors) + return profile_vals_corrected -def correct_sigma_with_boost_model(rvals, sigma_vals, boost_model="nfw_boost", **boost_model_kw): +def correct_with_boost_model(r_proj, profile_vals, boost_model, boost_rscale, **boost_model_kw): """Given a boost model and sigma profile, compute corrected sigma Parameters ---------- - rvals : array_like - radii - sigma_vals : array_like - uncorrected sigma with cluster member dilution - boost_model : str, optional + r_proj : array_like + Radii + profile_vals : array_like + Uncorrected profile values + boost_model : str Boost model to use for correcting sigma - * 'nfw_boost' - NFW profile model (Default) * 'powerlaw_boost' - Powerlaw profile + boost_rscale : float + Scale radius to be used with the boost model + Returns ------- - sigma_corrected : numpy.ndarray - correted radial profile + profile_vals_corrected : numpy.ndarray + Correted radial profile """ boost_model_func = boost_models[boost_model] - boost_factors = boost_model_func(rvals, **boost_model_kw) + boost_factors = boost_model_func(r_proj, boost_rscale, **boost_model_kw) - sigma_corrected = np.array(sigma_vals) / boost_factors - return sigma_corrected + profile_vals_corrected = np.array(profile_vals) / boost_factors + return profile_vals_corrected boost_models = { diff --git a/examples/demo_boost_factors.ipynb b/examples/demo_boost_factors.ipynb index 6f236bd22..d7f8b320e 100644 --- a/examples/demo_boost_factors.ipynb +++ b/examples/demo_boost_factors.ipynb @@ -2,6 +2,7 @@ "cells": [ { "cell_type": "markdown", + "id": "6afbf4b7", "metadata": {}, "source": [ "# Boost factors" @@ -10,6 +11,7 @@ { "cell_type": "code", "execution_count": null, + "id": "e3cd679b", "metadata": {}, "outputs": [], "source": [ @@ -24,6 +26,7 @@ }, { "cell_type": "markdown", + "id": "257a1892", "metadata": {}, "source": [ "Make sure we know which version we're using" @@ -32,6 +35,7 @@ { "cell_type": "code", "execution_count": null, + "id": "2ebc9e45", "metadata": {}, "outputs": [], "source": [ @@ -40,6 +44,7 @@ }, { "cell_type": "markdown", + "id": "67544ec1", "metadata": {}, "source": [ "### Define cosmology object" @@ -48,6 +53,7 @@ { "cell_type": "code", "execution_count": null, + "id": "cb23750b", "metadata": {}, "outputs": [], "source": [ @@ -56,6 +62,7 @@ }, { "cell_type": "markdown", + "id": "7b50a10e", "metadata": {}, "source": [ "First, we want to generate a $\\Delta\\Sigma$ (excess surface density) profile from mock data, to which we can apply boost factors. The mock data is generated in the following cells." @@ -63,6 +70,7 @@ }, { "cell_type": "markdown", + "id": "f4b319de", "metadata": {}, "source": [ "Generate cluster object from mock data" @@ -71,6 +79,7 @@ { "cell_type": "code", "execution_count": null, + "id": "2055117e", "metadata": {}, "outputs": [], "source": [ @@ -98,6 +107,7 @@ }, { "cell_type": "markdown", + "id": "8e24cea1", "metadata": {}, "source": [ "Loading this into a CLMM cluster object centered on (0,0)" @@ -106,6 +116,7 @@ { "cell_type": "code", "execution_count": null, + "id": "a62553b4", "metadata": {}, "outputs": [], "source": [ @@ -116,6 +127,7 @@ }, { "cell_type": "markdown", + "id": "a9f5c58f", "metadata": {}, "source": [ "Compute cross and tangential excess surface density for each source galaxy" @@ -124,6 +136,7 @@ { "cell_type": "code", "execution_count": null, + "id": "c961c355", "metadata": {}, "outputs": [], "source": [ @@ -141,6 +154,7 @@ }, { "cell_type": "markdown", + "id": "2aef9237", "metadata": {}, "source": [ "Calculate the binned profile" @@ -149,6 +163,7 @@ { "cell_type": "code", "execution_count": null, + "id": "8dbb5041", "metadata": {}, "outputs": [], "source": [ @@ -175,6 +190,7 @@ }, { "cell_type": "markdown", + "id": "d1e5e1d8", "metadata": {}, "source": [ "Plot the $\\Delta\\Sigma$ profile" @@ -183,6 +199,7 @@ { "cell_type": "code", "execution_count": null, + "id": "afb483c9", "metadata": {}, "outputs": [], "source": [ @@ -202,6 +219,7 @@ }, { "cell_type": "markdown", + "id": "8323cb3b", "metadata": {}, "source": [ "## Boost Factors" @@ -209,11 +227,13 @@ }, { "cell_type": "markdown", + "id": "de848b74", "metadata": {}, "source": [ "CLMM offers two boost models, the NFW boost model, and a powerlaw boost model. \n", "\n", - "Note that `compute_nfw_boost` requires two parameters to be specified, `rs` and `b0`, and `compute_powerlaw_boost` requires three paramters, `rs`, `b0` and `alpha`. The default values are in kpc. \n", + "- `compute_nfw_boost` requires two parameters to be specified, `rs` (non-optional) and `b0` (optional)\n", + "- `compute_powerlaw_boost` requires three paramters, `rs`(non-optional), `b0`(optional) and `alpha`(optional). \n", "\n", "Details on these boost models can be found [here](https://cluster-toolkit.readthedocs.io/en/latest/source/boostfactors.html)\n", "\n", @@ -223,18 +243,21 @@ { "cell_type": "code", "execution_count": null, + "id": "fab7bd2a", "metadata": {}, "outputs": [], "source": [ - "nfw_boost = u.compute_nfw_boost(cl.DeltaSigma_profile[\"radius\"], rscale=1000, boost0=0.1)\n", + "r_proj = cl.DeltaSigma_profile[\"radius\"]\n", + "r_scale = 1000 # in kpc to match r_proj units\n", + "nfw_boost = u.compute_nfw_boost(r_proj, r_scale, boost0=0.1)\n", "\n", - "powerlaw_boost = u.compute_powerlaw_boost(\n", - " cl.DeltaSigma_profile[\"radius\"], rscale=1000, boost0=0.1, alpha=-1.0\n", + "powerlaw_boost = u.compute_powerlaw_boost(r_proj, r_scale, boost0=0.1, alpha=-1.0\n", ")" ] }, { "cell_type": "markdown", + "id": "2e356f74", "metadata": {}, "source": [ "Plot the two boost factors, $B(R)$" @@ -243,11 +266,12 @@ { "cell_type": "code", "execution_count": null, + "id": "05fd04e0", "metadata": {}, "outputs": [], "source": [ - "plt.plot(cl.DeltaSigma_profile[\"radius\"], nfw_boost, label=\"NFW boost factor\")\n", - "plt.plot(cl.DeltaSigma_profile[\"radius\"], powerlaw_boost, label=\"Powerlaw boost factor\")\n", + "plt.loglog(r_proj, nfw_boost, label=\"NFW boost factor\", marker='.')\n", + "plt.loglog(r_proj, powerlaw_boost, label=\"Powerlaw boost factor\", marker='.')\n", "plt.xlabel(\"Radius [kpc]\")\n", "plt.ylabel(\"$B(R)$\")\n", "plt.legend()\n", @@ -256,6 +280,7 @@ }, { "cell_type": "markdown", + "id": "e48db29d", "metadata": {}, "source": [ "The $\\Delta\\Sigma$ profiles can be corrected with the boost factor using `correct_sigma_with_boost_values` or `correct_sigma_with_boost_model`. \n", @@ -266,6 +291,7 @@ }, { "cell_type": "markdown", + "id": "f6e65d7b", "metadata": {}, "source": [ "\n", @@ -285,6 +311,7 @@ }, { "cell_type": "markdown", + "id": "e57a94d2", "metadata": {}, "source": [ "First we will apply the boost factor with `correct_sigma_with_boost_values`" @@ -293,19 +320,21 @@ { "cell_type": "code", "execution_count": null, + "id": "90b9c8d2", "metadata": {}, "outputs": [], "source": [ - "Sigma_corrected_powerlaw_boost = u.correct_sigma_with_boost_values(\n", + "Sigma_corrected_powerlaw_boost = u.correct_with_boost_values(\n", " cl.DeltaSigma_profile[\"DeltaSigma_tan\"], powerlaw_boost\n", ")\n", - "Sigma_corrected_nfw_boost = u.correct_sigma_with_boost_values(\n", + "Sigma_corrected_nfw_boost = u.correct_with_boost_values(\n", " cl.DeltaSigma_profile[\"DeltaSigma_tan\"], nfw_boost\n", ")" ] }, { "cell_type": "markdown", + "id": "28040f51", "metadata": {}, "source": [ "Plot the result" @@ -314,6 +343,7 @@ { "cell_type": "code", "execution_count": null, + "id": "65fcd8b7", "metadata": {}, "outputs": [], "source": [ @@ -342,7 +372,7 @@ " color=\"grey\",\n", ")\n", "\n", - "# plt.loglog()\n", + "#plt.loglog()\n", "plt.title(\"DeltaSigma profile\")\n", "plt.xlabel(\"Radius [kpc]\")\n", "plt.ylabel(\"$\\Delta\\Sigma [M_\\odot\\; Mpc^{-2}]$\")\n", @@ -352,26 +382,30 @@ }, { "cell_type": "markdown", + "id": "67c5e14a", "metadata": {}, "source": [ - "Now the same again but with `correct_sigma_with_boost_model`" + "Now the same again but with `correct_with_boost_model`" ] }, { "cell_type": "code", "execution_count": null, + "id": "115a8687", "metadata": {}, "outputs": [], "source": [ - "Sigma_corrected_powerlaw_boost = u.correct_sigma_with_boost_model(\n", + "Sigma_corrected_powerlaw_boost = u.correct_with_boost_model(\n", " cl.DeltaSigma_profile[\"radius\"],\n", " cl.DeltaSigma_profile[\"DeltaSigma_tan\"],\n", - " boost_model=\"powerlaw_boost\",\n", + " \"powerlaw_boost\", # boost_model\n", + " r_scale, # boost_rscale (in units of cl.DeltaSigma_profile[\"radius\"])\n", ")\n", - "Sigma_corrected_nfw_boost = u.correct_sigma_with_boost_model(\n", + "Sigma_corrected_nfw_boost = u.correct_with_boost_model(\n", " cl.DeltaSigma_profile[\"radius\"],\n", " cl.DeltaSigma_profile[\"DeltaSigma_tan\"],\n", - " boost_model=\"nfw_boost\",\n", + " \"nfw_boost\", # boost_model\n", + " r_scale, # boost_rscale (in units of cl.DeltaSigma_profile[\"radius\"])\n", ")\n", "\n", "plt.errorbar(\n", @@ -410,9 +444,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "wrk", "language": "python", - "name": "python3" + "name": "wrk" }, "language_info": { "codemirror_mode": { @@ -424,7 +458,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.9" + "version": "3.10.9" } }, "nbformat": 4, diff --git a/tests/test_utils.py b/tests/test_utils.py index a70524472..f4d41f2cb 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -23,8 +23,9 @@ def test_compute_nfw_boost(): """Test the nfw model for boost factor""" # Test data rvals = np.arange(1, 11).tolist() + rscale = 1000 - boost_factors = utils.compute_nfw_boost(rvals) + boost_factors = utils.compute_nfw_boost(rvals, rscale) test_boost_factors = np.array( [ @@ -49,8 +50,9 @@ def test_compute_powerlaw_boost(): """Test the powerlaw model for boost factor""" # Test data rvals = np.arange(1, 11).tolist() # Cannot contain 0 due to reciprocal term + rscale = 1000 - boost_factors = utils.compute_powerlaw_boost(rvals) + boost_factors = utils.compute_powerlaw_boost(rvals, rscale) test_boost_factors = np.array( [101.0, 51.0, 34.33333333, 26.0, 21.0, 17.66666667, 15.28571429, 13.5, 12.11111111, 11.0] @@ -60,7 +62,7 @@ def test_compute_powerlaw_boost(): assert_allclose(boost_factors, test_boost_factors) -def test_correct_sigma_with_boost_values(): +def test_correct_with_boost_values(): """ """ # Make test data rvals = np.arange(1, 11) @@ -68,26 +70,29 @@ def test_correct_sigma_with_boost_values(): test_unit_boost_factors = np.ones(rvals.shape).tolist() - corrected_sigma = utils.correct_sigma_with_boost_values(sigma_vals, test_unit_boost_factors) + corrected_sigma = utils.correct_with_boost_values(sigma_vals, test_unit_boost_factors) assert_allclose(sigma_vals, corrected_sigma) -def test_correct_sigma_with_boost_model(): +def test_correct_with_boost_model(): """ """ # Make test data rvals = np.arange(1, 11).tolist() sigma_vals = (2 ** np.arange(10)).tolist() + boost_rscale = 1000 for boost_model in utils.boost_models.keys(): # Check for no nans or inf with positive-definite rvals and sigma vals assert np.all( np.isfinite( - utils.correct_sigma_with_boost_model(rvals, sigma_vals, boost_model=boost_model) + utils.correct_with_boost_model(rvals, sigma_vals, boost_model, boost_rscale) ) ) # Test requesting unsupported boost model - assert_raises(KeyError, utils.correct_sigma_with_boost_model, rvals, sigma_vals, "glue") + assert_raises( + KeyError, utils.correct_with_boost_model, rvals, sigma_vals, "glue", boost_rscale + ) def test_compute_radial_averages():