diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py index eab62c280..98e8ca0c0 100644 --- a/clmm/dataops/__init__.py +++ b/clmm/dataops/__init__.py @@ -1,4 +1,5 @@ """Data operation for polar/azimuthal averages in radial bins and weights""" + import warnings import numpy as np from astropy.coordinates import SkyCoord @@ -25,7 +26,7 @@ def compute_tangential_and_cross_components( dec_source, shear1, shear2, - coordinate_system="euclidean", + coordinate_system=None, geometry="curve", is_deltasigma=False, sigma_c=None, @@ -95,7 +96,7 @@ def compute_tangential_and_cross_components( coordinate_system: str, optional Coordinate system of the ellipticity components. Must be either 'celestial' or euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details. - Default is 'euclidean'. + If not set, defaults to 'euclidean'. geometry: str, optional Sky geometry to compute angular separation. Options are curve (uses astropy) or flat. @@ -120,6 +121,10 @@ def compute_tangential_and_cross_components( # pylint: disable-msg=too-many-locals # Note: we make these quantities to be np.array so that a name is not passed from astropy # columns + if coordinate_system is None: + warnings.warn("The coordinate_system argument was not provided. Defaulting to 'euclidean'.") + coordinate_system = "euclidean" + if validate_input: _validate_ra(locals(), "ra_source", True) _validate_dec(locals(), "dec_source", True) @@ -129,7 +134,7 @@ def compute_tangential_and_cross_components( validate_argument(locals(), "shear2", "float_array") validate_argument(locals(), "geometry", str) validate_argument(locals(), "sigma_c", "float_array", none_ok=True) - _validate_coordinate_system(locals(), "coordinate_system", str) + _validate_coordinate_system(locals(), "coordinate_system") ra_source_, dec_source_, shear1_, shear2_ = arguments_consistency( [ra_source, dec_source, shear1, shear2], names=("Ra", "Dec", "Shear1", "Shear2"), @@ -482,6 +487,7 @@ def make_radial_profile( z_lens=None, validate_input=True, weights=None, + coordinate_system=None, ): r"""Compute the angular profile of given components @@ -535,6 +541,10 @@ def make_radial_profile( weights: array-like, optional Array of individual galaxy weights. If specified, the radial binned profile is computed using a weighted average + coordinate_system: str, optional + Coordinate system of the ellipticity components. Must be either 'celestial' or + euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details. + Default is 'euclidean'. Returns ------- @@ -553,6 +563,9 @@ def make_radial_profile( module. """ # pylint: disable-msg=too-many-locals + if not coordinate_system: + warnings.warn("coordinate_system not set, defaulting to 'euclidean'.") + coordinate_system = "euclidean" if validate_input: validate_argument(locals(), "angsep", "float_array") validate_argument(locals(), "angsep_units", str) @@ -564,6 +577,7 @@ def make_radial_profile( arguments_consistency(components, names=comp_dict.keys(), prefix="Input components") for component in comp_dict: validate_argument(comp_dict, component, "float_array") + _validate_coordinate_system(locals(), "coordinate_system") # Check to see if we need to do a unit conversion if angsep_units is not bin_units: source_seps = convert_units(angsep, angsep_units, bin_units, redshift=z_lens, cosmo=cosmo) @@ -576,7 +590,7 @@ def make_radial_profile( profile_table = GCData( [bins[:-1], np.zeros(len(bins) - 1), bins[1:]], names=("radius_min", "radius", "radius_max"), - meta={"bin_units": bin_units}, # Add metadata + meta={"bin_units": bin_units, "coordinate_system": coordinate_system}, # Add metadata ) # Compute the binned averages and associated errors for i, component in enumerate(components): diff --git a/clmm/galaxycluster.py b/clmm/galaxycluster.py index 35b4f44c6..57763e499 100644 --- a/clmm/galaxycluster.py +++ b/clmm/galaxycluster.py @@ -18,7 +18,6 @@ _validate_ra, _validate_dec, _draw_random_points_from_tab_distribution, - _validate_coordinate_system, ) @@ -37,10 +36,6 @@ class GalaxyCluster: Redshift of galaxy cluster center galcat : GCData Table of background galaxy data containing at least galaxy_id, ra, dec, e1, e2, z - coordinate_system : str, optional - Coordinate system of the ellipticity components. Must be either 'celestial' or - euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details. - Default is 'euclidean'. validate_input: bool Validade each input argument """ @@ -64,7 +59,6 @@ def _add_values( dec: float, z: float, galcat: GCData, - coordinate_system: str = "euclidean", ): """Add values for all attributes""" self.unique_id = unique_id @@ -72,7 +66,6 @@ def _add_values( self.dec = dec self.z = z self.galcat = galcat - self.coordinate_system = coordinate_system def _check_types(self): """Check types of all attributes""" @@ -81,7 +74,6 @@ def _check_types(self): _validate_dec(vars(self), "dec", False) validate_argument(vars(self), "z", (float, str), argmin=0, eqmin=True) validate_argument(vars(self), "galcat", GCData) - _validate_coordinate_system(vars(self), "coordinate_system", str) self.unique_id = str(self.unique_id) self.ra = float(self.ra) self.dec = float(self.dec) @@ -297,7 +289,7 @@ def compute_tangential_and_cross_components( dec_lens=self.dec, geometry=geometry, validate_input=self.validate_input, - coordinate_system=self.coordinate_system, + coordinate_system=self.galcat.meta["coordinate_system"], **cols, ) if add: @@ -607,6 +599,7 @@ def make_radial_profile( for n in (tan_component_in_err, cross_component_in_err, None) ], weights=self.galcat[weights_in].data if use_weights else None, + coordinate_system=self.galcat.meta["coordinate_system"], ) # Reaname table columns for i, name in enumerate([tan_component_out, cross_component_out, "z"]): @@ -722,3 +715,20 @@ def set_ra_lower(self, ra_low=0): if "ra" in self.galcat.columns: self.galcat["ra"][self.galcat["ra"] < ra_low] += 360.0 self.galcat["ra"][self.galcat["ra"] >= ra_low + 360.0] -= 360.0 + + def update_coordinate_system(self, coordinate_system, *args): + """Updates coordinate_system metadata of the galcat and converts ellipticity + + Parameters + ---------- + coordinate_system : str + Coordinate system of the ellipticity components. Must be either 'celestial' or + euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details. + *args : tuple + Components to be converted, must be in data. + + Returns + ------- + None + """ + self.galcat.update_coordinate_system(coordinate_system, *args) diff --git a/clmm/gcdata.py b/clmm/gcdata.py index 8ed47f44f..7b12b5663 100644 --- a/clmm/gcdata.py +++ b/clmm/gcdata.py @@ -1,16 +1,18 @@ """ Define the custom data type """ + import warnings from collections import OrderedDict from astropy.table import Table as APtable -import numpy as np - import qp +import numpy as np +from .utils import _validate_coordinate_system class GCMetaData(OrderedDict): - r"""Object to store metadata, it always has a cosmo key with protective changes + r"""Object to store metadata, it always has a cosmo and coordinate_system keys with protective + changes Attributes ---------- @@ -29,6 +31,10 @@ def __setitem__(self, item, value): raise ValueError( "cosmo must be changed via update_cosmo or update_cosmo_ext_valid method" ) + if item == "coordinate_system" and self.get("coordinate_system"): + raise ValueError( + "coordinate_system must be changed via update_coordinate_system method" + ) OrderedDict.__setitem__(self, item, value) def __getitem__(self, item): @@ -43,7 +49,7 @@ def __getitem__(self, item): class GCData(APtable): """ - GCData: A data objetc for gcdata. Right now it behaves as an astropy table, with the following + GCData: A data object for gcdata. Right now it behaves as an astropy table, with the following modifications: `__getitem__` is case independent; The attribute .meta['cosmo'] is protected and can only be changed via update_cosmo or update_cosmo_ext_valid methods; @@ -65,6 +71,16 @@ def __init__(self, *args, **kwargs): APtable.__init__(self, *args, **kwargs) metakwargs = kwargs["meta"] if "meta" in kwargs else {} metakwargs = {} if metakwargs is None else metakwargs + if ( + "coordinate_system" not in metakwargs + and "copy_indices" not in kwargs + and "masked" not in kwargs + and "copy" not in kwargs + ): + warnings.warn("coordinate_system not set, defaulting to 'euclidean'") + metakwargs["coordinate_system"] = "euclidean" + if "coordinate_system" in metakwargs: + _validate_coordinate_system(metakwargs, "coordinate_system") self.meta = GCMetaData(**metakwargs) # this attribute is set when source galaxies have p(z) self.pzpdf_info = { @@ -89,7 +105,7 @@ def _str_pzpdf_info(self): out += " " + str(np.round(self.pzpdf_info.get("zbins"), 2)) np.set_printoptions(**default_cfg) elif out == "quantiles": - np.set_printoptions(formatter={'float': "{0:g}".format}, edgeitems=3, threshold=6) + np.set_printoptions(formatter={"float": "{0:g}".format}, edgeitems=3, threshold=6) out += " " + str(self.pzpdf_info["quantiles"]) out += " - unpacked with zgrid : " + str( self.pzpdf_info["unpack_quantile_zbins_limits"] @@ -171,6 +187,8 @@ def update_info_ext_valid(self, key, gcdata, ext_value, overwrite=False): ------- None """ + if key == "coordinate_system": + _validate_coordinate_system(locals(), "ext_value") if ext_value: in_value = gcdata.meta[key] if in_value and in_value != ext_value: @@ -221,6 +239,30 @@ def update_cosmo(self, cosmo, overwrite=False): """ self.update_cosmo_ext_valid(self, cosmo, overwrite=overwrite) + def update_coordinate_system(self, coordinate_system, *args): + r"""Updates coordinate_system metadata and converts ellipticity + + Parameters + ---------- + coordinate_system : str + Coordinate system of the ellipticity components. Must be either 'celestial' or + euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details. + *args : tuple + Components to be converted, must be in data. + + Returns + ------- + None + """ + in_coordinate_system = self.meta["coordinate_system"] + self.update_info_ext_valid("coordinate_system", self, coordinate_system, overwrite=True) + if coordinate_system != in_coordinate_system: + for col in args: + if col in self.colnames: + self[col] = -self[col] + else: + raise ValueError(f"component {col} not found in data") + def has_pzpdfs(self): """Get pzbins and pzpdfs of galaxies diff --git a/clmm/support/mock_data.py b/clmm/support/mock_data.py index 27858b822..71a0b800d 100644 --- a/clmm/support/mock_data.py +++ b/clmm/support/mock_data.py @@ -44,7 +44,7 @@ def generate_galaxy_catalog( pz_bins=101, pz_quantiles_conf=(5, 31), pzpdf_type="shared_bins", - coordinate_system="euclidean", + coordinate_system=None, validate_input=True, ): r"""Generates a mock dataset of sheared background galaxies. @@ -170,7 +170,7 @@ def generate_galaxy_catalog( coordinate_system : str, optional Coordinate system of the ellipticity components. Must be either 'celestial' or euclidean'. See https://doi.org/10.48550/arXiv.1407.7676 section 5.1 for more details. - Default is 'euclidean'. + If not set, defaults to 'euclidean'. validate_input: bool Validade each input argument @@ -188,6 +188,10 @@ def generate_galaxy_catalog( # Too many local variables (25/15) # pylint: disable=R0914 + if not coordinate_system: + warnings.warn("The coordinate_system argument was not provided. Defaulting to 'euclidean'.") + coordinate_system = "euclidean" + if validate_input: validate_argument(locals(), "cluster_m", float, argmin=0, eqmin=True) validate_argument(locals(), "cluster_z", float, argmin=0, eqmin=True) @@ -224,7 +228,7 @@ def generate_galaxy_catalog( validate_argument(locals(), "ngals", float, none_ok=True) validate_argument(locals(), "ngal_density", float, none_ok=True) validate_argument(locals(), "pz_bins", (int, "array")) - _validate_coordinate_system(locals(), "coordinate_system", str) + _validate_coordinate_system(locals(), "coordinate_system") if zsrc_min is None: zsrc_min = cluster_z + 0.1 @@ -393,7 +397,7 @@ def _generate_galaxy_catalog( # pylint: disable=R0914 # Set the source galaxy redshifts - galaxy_catalog = _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals) + galaxy_catalog = _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals, coordinate_system) # Add photo-z errors and pdfs to source galaxy redshifts if photoz_sigma_unscaled is not None: @@ -486,7 +490,7 @@ def _generate_galaxy_catalog( return galaxy_catalog[cols] -def _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals): +def _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals, coordinate_system): """Set source galaxy redshifts either set to a fixed value or draw from a predefined distribution. Return a table (GCData) of the source galaxies @@ -542,7 +546,9 @@ def _draw_source_redshifts(zsrc, zsrc_min, zsrc_max, ngals): else: raise ValueError(f"zsrc must be a float, chang13 or desc_srd. You set: {zsrc}") - return GCData([zsrc_list, zsrc_list], names=("ztrue", "z")) + return GCData( + [zsrc_list, zsrc_list], names=("ztrue", "z"), meta={"coordinate_system": coordinate_system} + ) def _compute_photoz_pdfs( diff --git a/clmm/utils/validation.py b/clmm/utils/validation.py index 1002ed222..0e775b704 100644 --- a/clmm/utils/validation.py +++ b/clmm/utils/validation.py @@ -226,24 +226,21 @@ def _validate_is_deltasigma_sigma_c(is_deltasigma, sigma_c): if not is_deltasigma and sigma_c is not None: raise TypeError(f"sigma_c (={sigma_c}) must be None when is_deltasigma=False") -def _validate_coordinate_system(loc, coordinate_system, valid_type): + +def _validate_coordinate_system(loc, argname): r"""Validate the coordinate system. Parameters ---------- loc: dict Dictionary with all input arguments. Should be locals(). - coordinate_system: str - Coordinate system of the ellipticity components. Must be either 'celestial' or 'euclidean'. - valid_type: str, type - Valid types for argument, options are object types, list/tuple of types, or: - - * 'int_array' - interger, interger array - * 'float_array' - float, float array + argname: str + Name of argument to be tested. """ - validate_argument(loc, coordinate_system, valid_type) - if loc[coordinate_system] not in ["celestial", "euclidean"]: - raise ValueError(f"{coordinate_system} must be 'celestial' or 'euclidean'.") + validate_argument(loc, argname, str) + if loc[argname] not in ["celestial", "euclidean"]: + raise ValueError(f"{argname} must be 'celestial' or 'euclidean'.") + class DiffArray: """Array where arr1==arr2 is actually all(arr1==arr)""" diff --git a/docs/doc-config.ini b/docs/doc-config.ini index cb9d24c5c..278ba1063 100644 --- a/docs/doc-config.ini +++ b/docs/doc-config.ini @@ -33,6 +33,6 @@ OTHER ../examples/demo_compute_deltasigma_weights.ipynb ../examples/demo_mass_conversion.ipynb ../examples/demo_mock_ensemble_realistic.ipynb -../examples/demo_coordinate_system_datasets.ipynb +../examples/demo_coordinate_system.ipynb #../examples/other_compare_NFW_critical_massdef.ipynb #../examples/other_flat_sky_limitations.ipynb diff --git a/examples/demo_coordinate_system.ipynb b/examples/demo_coordinate_system.ipynb new file mode 100644 index 000000000..1770383fc --- /dev/null +++ b/examples/demo_coordinate_system.ipynb @@ -0,0 +1,1285 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0711213a", + "metadata": {}, + "source": [ + "# Tests of coordinate system effects on shear profiles" + ] + }, + { + "cell_type": "markdown", + "id": "3a2b9d35", + "metadata": {}, + "source": [ + "Authors: Marina Ricci, Tomomi Sunayama\n", + "\n", + "Tested, modified, and documented by: Camille Avestruz, Caio Lima de Oliveira\n", + "\n", + "In this notebook we illustrate the importance of setting the correct coordinate system for shear catalogs. We start by showing these effects on generated mock data, then on HSC Y3 data.\n", + "\n", + "Throughout this notebook we also show how to correctly set the coordinate system and how to update it and convert data." + ] + }, + { + "cell_type": "markdown", + "id": "1d67ea10", + "metadata": {}, + "source": [ + "### Shear coordinate system definition and conversion" + ] + }, + { + "cell_type": "markdown", + "id": "3b03619a", + "metadata": {}, + "source": [ + "* Celestial coordinate system: the declination $\\delta$ and the right ascension $\\alpha$ take the role of the spherical angles $\\theta$ and $\\varphi$.\n", + "\n", + "* Euclidean coordinate system: a cartesian coordinate system defined on the plane tanget to the celestial sphere at the point of observation. Here, the $y$-axis is parallel to the declination $\\delta$ and the $x$-axis is antiparallel to the right ascension $\\alpha$.\n", + "\n", + "In a small angles, planar approximation of the Celestial coordinates, both coordinate systems are related by a parity transformation of the $x$-axis. The conversion between the Euclidean ellipticity $\\epsilon^E = \\epsilon_1^E + i \\epsilon_2^E$ and the Celestial ellipticity $\\epsilon^C = \\epsilon_1^C + i \\epsilon_2^C$ is then given by the transformation $\\varphi \\rightarrow \\pi - \\varphi$:\n", + "\n", + "$$\\epsilon^C_1 + i \\epsilon_2^C = |\\epsilon| e^{2 i \\varphi^\\prime} = |\\epsilon| e^{2 i (\\pi - \\varphi)} = |\\epsilon| e^{- 2 i \\varphi} = \\epsilon^E_1 - i \\epsilon_2^E$$" + ] + }, + { + "cell_type": "markdown", + "id": "28a917d8-04b1-4a4d-9f3d-3bf0be1b5dba", + "metadata": {}, + "source": [ + "### Here we generate mock source catalogs with different coordinate system and explore how that must be accounted to measure the correct shear profiles. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "423b5f84-831a-4da7-af39-92e1b6dc3818", + "metadata": {}, + "outputs": [], + "source": [ + "import clmm\n", + "\n", + "clmm.__version__" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "9ed8f467-dbaa-413f-a2cf-f1ff7fd5513d", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import pandas as pd\n", + "from astropy.io import fits\n", + "from scipy import spatial\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "db760451-21e7-4c5b-8fdf-af066f2b8dd9", + "metadata": {}, + "outputs": [], + "source": [ + "from clmm.support import mock_data as mock\n", + "from clmm.utils import convert_units\n", + "from clmm import Cosmology" + ] + }, + { + "cell_type": "markdown", + "id": "c31afe8f-062a-42e2-bc79-442ac06579e4", + "metadata": {}, + "source": [ + "## Generate the mock catalog with different source galaxy options\n", + "In this example, the mock data includes: shape noise, galaxies drawn from redshift distribution and photoz errors." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "7c7b2a1a-e59e-4f8e-a415-d957af0b4b0f", + "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": "code", + "execution_count": 5, + "id": "11263d57-3e76-4fe1-adbc-9e1299e6900e", + "metadata": {}, + "outputs": [], + "source": [ + "cosmo = mock_cosmo\n", + "cluster_id = \"Awesome_cluster\"\n", + "\n", + "cluster_m = 1.0e15 # M200,m\n", + "cluster_z = 0.3\n", + "# Cluster centre coordinates\n", + "cluster_ra = 50.0\n", + "cluster_dec = 87.0\n", + "concentration = 4" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "20841454-8b80-4b8b-adb9-c1f189f03c78", + "metadata": {}, + "outputs": [], + "source": [ + "# let's put all these quantities in a single dictionary to facilitate clarity\n", + "cluster_kwargs = {\n", + " \"cluster_m\": cluster_m,\n", + " \"cluster_z\": cluster_z,\n", + " \"cluster_ra\": cluster_ra,\n", + " \"cluster_dec\": cluster_dec,\n", + " \"cluster_c\": concentration,\n", + " \"cosmo\": cosmo,\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "bb8e1555-565e-40c3-8d76-b5039bf6043d", + "metadata": {}, + "outputs": [], + "source": [ + "# let's put all these quantities in a single dictionary to facilitate clarity\n", + "source_kwargs = {\n", + " \"zsrc\": \"chang13\",\n", + " \"zsrc_min\": cluster_z + 0.1,\n", + " \"photoz_sigma_unscaled\": 0.05,\n", + " \"ngals\": 1000,\n", + " \"pz_bins\": np.linspace(0, 10, 1001),\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "2fd4f56e", + "metadata": {}, + "source": [ + "We must supply the coordinate system information when generating a mock galaxy catalog. If we don't, a warning is issued and problems may arise!" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "4db4311b-7947-4089-b092-d69de52bff69", + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(679)\n", + "\n", + "mock_sources_euclidean_coord = mock.generate_galaxy_catalog(\n", + " **cluster_kwargs, **source_kwargs, coordinate_system=\"euclidean\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "e51cb813-c077-4d94-878e-0387e1b21801", + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(679)\n", + "\n", + "mock_sources_celestial_coord = mock.generate_galaxy_catalog(\n", + " **cluster_kwargs, **source_kwargs, coordinate_system=\"celestial\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec9fab17", + "metadata": {}, + "outputs": [], + "source": [ + "np.random.seed(679)\n", + "\n", + "# In this case, we are going to generate a mock catalog without setting\n", + "# the coordinate system, which will raise a warning and default to \"euclidean\".\n", + "mock_sources_default_coord = mock.generate_galaxy_catalog(\n", + " **cluster_kwargs,\n", + " **source_kwargs,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "1b1cfb30", + "metadata": {}, + "source": [ + "Note that $e_1$ remains the same for all catalogs, while $e_2$ changes its sign between coordinate systems. Also, notice the `coordinate_system` metadata!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc2a6f85-d71e-497e-875e-358e0c68cfad", + "metadata": {}, + "outputs": [], + "source": [ + "mock_sources_euclidean_coord[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d5fbe561-c072-41cf-abdf-14845f28c0c2", + "metadata": {}, + "outputs": [], + "source": [ + "mock_sources_celestial_coord[:5]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4f04aa68", + "metadata": {}, + "outputs": [], + "source": [ + "mock_sources_default_coord[:5]" + ] + }, + { + "cell_type": "markdown", + "id": "a352677d", + "metadata": {}, + "source": [ + "Now, we can easily convert between the diferente coordinates with the method `update_coordinate_system()`. However, pay attention to the fact we must supply which columns we want to update!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2969f85", + "metadata": {}, + "outputs": [], + "source": [ + "mock_sources_default_coord.update_coordinate_system(\"celestial\", \"e2\")\n", + "\n", + "mock_sources_default_coord[:5]" + ] + }, + { + "cell_type": "markdown", + "id": "588701b8", + "metadata": {}, + "source": [ + "We can also update the coordinate system but not convert the data. This is useful if the original coordinate system was incorrectly set but may also lead to errors (as we intentionally do now)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "935e75ba", + "metadata": {}, + "outputs": [], + "source": [ + "mock_sources_default_coord.update_coordinate_system(\"euclidean\")\n", + "\n", + "mock_sources_wrong_coord = mock_sources_default_coord\n", + "\n", + "mock_sources_wrong_coord[:5]" + ] + }, + { + "cell_type": "markdown", + "id": "7e852ccb-a49a-441f-89b9-da314fecd5ce", + "metadata": {}, + "source": [ + "## Generate cluster objects from mock data" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "6d76bb85-2c93-4778-8135-1f1e6c139689", + "metadata": {}, + "outputs": [], + "source": [ + "import clmm.dataops\n", + "from clmm.dataops import compute_tangential_and_cross_components, make_radial_profile, make_bins\n", + "from clmm.galaxycluster import GalaxyCluster" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "353b4dbb-6239-4763-9125-d5d53b2f7d4e", + "metadata": {}, + "outputs": [], + "source": [ + "cl_euclidean = GalaxyCluster(\n", + " cluster_id,\n", + " cluster_ra,\n", + " cluster_dec,\n", + " cluster_z,\n", + " mock_sources_euclidean_coord,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "20c559a0-0ce2-411c-8d12-dbcb3fc42026", + "metadata": {}, + "outputs": [], + "source": [ + "cl_celestial = GalaxyCluster(\n", + " cluster_id,\n", + " cluster_ra,\n", + " cluster_dec,\n", + " cluster_z,\n", + " mock_sources_celestial_coord,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "a275df25", + "metadata": {}, + "outputs": [], + "source": [ + "cl_wrong = GalaxyCluster(\n", + " cluster_id,\n", + " cluster_ra,\n", + " cluster_dec,\n", + " cluster_z,\n", + " mock_sources_wrong_coord,\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "efefdd51-eb7a-4e26-a2dd-da574f213932", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax1 = plt.subplots(1, 1)\n", + "\n", + "ax1.scatter(cl_euclidean.galcat[\"e1\"], cl_euclidean.galcat[\"e2\"], s=2, alpha=0.5, label=\"euclidean\")\n", + "ax1.scatter(\n", + " cl_celestial.galcat[\"e1\"],\n", + " cl_celestial.galcat[\"e2\"],\n", + " s=2,\n", + " alpha=0.5,\n", + " color=\"red\",\n", + " label=\"celestial\",\n", + ")\n", + "\n", + "ax1.set_xlabel(\"$\\\\epsilon_1$\")\n", + "ax1.set_ylabel(\"$\\\\epsilon_2$\")\n", + "ax1.set_aspect(\"equal\", \"datalim\")\n", + "ax1.set_xlim(-0.125, 0.125)\n", + "ax1.set_ylim(-0.125, 0.125)\n", + "ax1.axvline(0, linestyle=\"dotted\", color=\"black\")\n", + "ax1.axhline(0, linestyle=\"dotted\", color=\"black\")\n", + "\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "f134fc6f-fc61-462a-b283-51eb6cb0cfd2", + "metadata": {}, + "source": [ + "## Compute and plot shear profiles" + ] + }, + { + "cell_type": "markdown", + "id": "789c6e25", + "metadata": {}, + "source": [ + "We will now compute the tangential and cross components of the ellipticity for all three clusters. We should get the same result for `cl_euclidean` and `cl_celestial`, but a completely different dataset for `cl_wrong`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f8b46476-6ecc-4afe-a9a6-8aebaf077f0c", + "metadata": {}, + "outputs": [], + "source": [ + "cl_euclidean.compute_tangential_and_cross_components(add=True)\n", + "cl_euclidean.galcat[\"et\", \"ex\"].pprint(max_width=-1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43a27c24-be4a-4380-bc46-b49810af780d", + "metadata": {}, + "outputs": [], + "source": [ + "cl_celestial.compute_tangential_and_cross_components(add=True)\n", + "cl_celestial.galcat[\"et\", \"ex\"].pprint(max_width=-1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9632531", + "metadata": {}, + "outputs": [], + "source": [ + "cl_wrong.compute_tangential_and_cross_components(add=True)\n", + "cl_wrong.galcat[\"et\", \"ex\"].pprint(max_width=-1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "137d6945-67cc-42ef-b93b-07f3c20cad5f", + "metadata": {}, + "outputs": [], + "source": [ + "f, ax = plt.subplots(1, 2, figsize=(10, 4))\n", + "\n", + "ax[0].hist(cl_euclidean.galcat[\"et\"], bins=50, color=\"tab:blue\", alpha=0.5)\n", + "ax[0].hist(cl_celestial.galcat[\"et\"], bins=50, color=\"tab:red\", alpha=0.5)\n", + "ax[0].hist(cl_wrong.galcat[\"et\"], bins=50, color=\"tab:orange\", alpha=0.5)\n", + "ax[0].set_xlabel(\"$\\\\epsilon_t$\", fontsize=\"xx-large\")\n", + "\n", + "ax[1].hist(cl_euclidean.galcat[\"ex\"], bins=50, color=\"tab:blue\", alpha=0.5)\n", + "ax[1].hist(cl_celestial.galcat[\"ex\"], bins=50, color=\"tab:red\", alpha=0.5)\n", + "ax[1].hist(cl_wrong.galcat[\"ex\"], bins=50, color=\"tab:orange\", alpha=0.5)\n", + "ax[1].set_xlabel(\"$\\\\epsilon_x$\", fontsize=\"xx-large\")\n", + "ax[1].set_yscale(\"log\")\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8700b684-279f-485c-98e3-a25ca5d9bd2f", + "metadata": {}, + "outputs": [], + "source": [ + "cl_euclidean.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "cl_euclidean.profile.show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f6774006-189b-4226-a4dd-f31722b7713d", + "metadata": {}, + "outputs": [], + "source": [ + "cl_celestial.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "cl_celestial.profile.show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3dbc425f", + "metadata": {}, + "outputs": [], + "source": [ + "cl_wrong.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "cl_wrong.profile.show_in_notebook()" + ] + }, + { + "cell_type": "markdown", + "id": "f52c7dc3", + "metadata": {}, + "source": [ + "### => When the correct coordinate system is specified, the profiles coming from the two catalogs are identical." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ab100d0-589f-42bd-bb90-eaefedf331cb", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, height_ratios=[3, 1], sharex=True)\n", + "\n", + "ax[0].errorbar(\n", + " cl_euclidean.profile[\"radius\"],\n", + " cl_euclidean.profile[\"gt\"],\n", + " yerr=cl_euclidean.profile[\"gt_err\"],\n", + " alpha=0.5,\n", + " marker=\".\",\n", + " color=\"tab:red\",\n", + " label=\"euclidean\",\n", + ")\n", + "ax[1].errorbar(\n", + " cl_euclidean.profile[\"radius\"],\n", + " cl_euclidean.profile[\"gx\"],\n", + " yerr=cl_euclidean.profile[\"gx_err\"],\n", + " marker=\".\",\n", + " alpha=0.5,\n", + " color=\"tab:red\",\n", + ")\n", + "\n", + "ax[0].errorbar(\n", + " cl_celestial.profile[\"radius\"] * 1.02,\n", + " cl_celestial.profile[\"gt\"],\n", + " yerr=cl_celestial.profile[\"gt_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + " label=\"celestial\",\n", + ")\n", + "ax[1].errorbar(\n", + " cl_celestial.profile[\"radius\"] * 1.02,\n", + " cl_celestial.profile[\"gx\"],\n", + " yerr=cl_celestial.profile[\"gx_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + ")\n", + "\n", + "ax[0].legend()\n", + "ax[0].set_xscale(\"log\")\n", + "ax[1].set_xlabel(\"R [kpc]\")\n", + "ax[0].set_ylabel(\"$g_t$\")\n", + "ax[1].set_ylabel(\"$g_x$\")\n", + "\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "43c7340f", + "metadata": {}, + "source": [ + "### => However, when the coordinate system is not correctly specified, the profiles are incorrect." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b098be6-8364-45c3-aa5a-ce28a84beb03", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, height_ratios=[3, 1], sharex=True)\n", + "\n", + "ax[0].errorbar(\n", + " cl_wrong.profile[\"radius\"],\n", + " cl_wrong.profile[\"gt\"],\n", + " yerr=cl_wrong.profile[\"gt_err\"],\n", + " alpha=0.5,\n", + " marker=\".\",\n", + " color=\"tab:red\",\n", + " label=\"incorrect coordinate system\",\n", + ")\n", + "ax[1].errorbar(\n", + " cl_wrong.profile[\"radius\"],\n", + " cl_wrong.profile[\"gx\"],\n", + " yerr=cl_wrong.profile[\"gx_err\"],\n", + " marker=\".\",\n", + " alpha=0.5,\n", + " color=\"tab:red\",\n", + ")\n", + "\n", + "ax[0].legend()\n", + "ax[0].set_xscale(\"log\")\n", + "ax[1].set_xlabel(\"R [kpc]\")\n", + "ax[0].set_ylabel(\"$g_t$\")\n", + "ax[1].set_ylabel(\"$g_x$\")\n", + "\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "f735bc6a", + "metadata": {}, + "source": [ + "### To recover the correct profile we need to update the coordinate system, compute the tangential and cross components again and recalculate the profile." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "acdfb0e0", + "metadata": {}, + "outputs": [], + "source": [ + "cl_wrong.update_coordinate_system(\"celestial\")\n", + "cl_wrong.compute_tangential_and_cross_components(add=True)\n", + "\n", + "cl_wrong.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "cl_wrong.profile.show_in_notebook()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2d034194", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, height_ratios=[3, 1], sharex=True)\n", + "\n", + "ax[0].errorbar(\n", + " cl_wrong.profile[\"radius\"],\n", + " cl_wrong.profile[\"gt\"],\n", + " yerr=cl_wrong.profile[\"gt_err\"],\n", + " alpha=0.5,\n", + " marker=\".\",\n", + " color=\"tab:red\",\n", + " label=\"correct coordinate system\",\n", + ")\n", + "ax[1].errorbar(\n", + " cl_wrong.profile[\"radius\"],\n", + " cl_wrong.profile[\"gx\"],\n", + " yerr=cl_wrong.profile[\"gx_err\"],\n", + " marker=\".\",\n", + " alpha=0.5,\n", + " color=\"tab:red\",\n", + ")\n", + "\n", + "ax[0].legend()\n", + "ax[0].set_xscale(\"log\")\n", + "ax[1].set_xlabel(\"R [kpc]\")\n", + "ax[0].set_ylabel(\"$g_t$\")\n", + "ax[1].set_ylabel(\"$g_x$\")\n", + "\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8f7908a7", + "metadata": {}, + "source": [ + "## Let's now do the same test on real data" + ] + }, + { + "cell_type": "markdown", + "id": "b1572733", + "metadata": {}, + "source": [ + "Here we present three datasets, each with different coordinate systems:\n", + "1. CosmoDC2 source galaxies with shears extracted from `TXPipe` for a single galaxy cluster (data is in `euclidean` coordinates);\n", + "2. Example source galaxies for galaxy clusters from a [Summer School](https://github.com/oguri/wlcluster_tutorial) taught by Masamune Oguri (data is also in `euclidean` coordinates);\n", + "3. HSC Y3 source galaxies with shears post processed by Tomomi Sunayama (data is in `celestial` coordinates)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instructions to download text data\n", + "\n", + "First, create a directory where you want to put the example data, e.g. for a given `data_coords_dir`:\n", + "\n", + "```\n", + "mkdir -p /data_coords\n", + "cd /data_coords\n", + "```\n", + "\n", + "Download all files from this [dropbox link](https://www.dropbox.com/scl/fo/dwsccslr5iwb7lnkf8jvx/AJkjgFeemUEHpHaZaHHqpAg?rlkey=efbtsr15mdrs3y6xsm7l48o0r&st=xb58ap0g&dl=0). This will be a zip file, `CLMM_data.zip` of size 242 Mb. `scp` or `mv` this to `data_coords_dir`. From the directory, you should be able to unzip:\n", + "\n", + "```\n", + "unzip data_CLMM.zip -d .\n", + "```\n", + "\n", + "You now have the necessary data files to run this notebook. **Make sure to change the `data_coords_dir` variable in the cell below to the appropriate location where you unzipped these files.**\n" + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "2f61c348", + "metadata": {}, + "outputs": [], + "source": [ + "# CHANGE TO YOUR LOCATION\n", + "data_coords_dir = \"/data_coords\"" + ] + }, + { + "cell_type": "markdown", + "id": "b1e6b530", + "metadata": {}, + "source": [ + "### Example galaxy cluster from CosmoDC2\n", + "\n", + "Here, we plot an example galaxy cluster shear profile using cluster and source galaxy files generated from CosmoDC2 and processed through TXPipe. These are originally in `euclidean` coordinates. Let's start by importing the data:" + ] + }, + { + "cell_type": "code", + "execution_count": 33, + "id": "62569cbd", + "metadata": {}, + "outputs": [], + "source": [ + "dc2_cluster = pd.read_pickle(data_coords_dir + \"/test_cluster.pkl\")\n", + "\n", + "cluster_z = dc2_cluster[\"redshift\"] # Cluster redshift\n", + "cluster_ra = dc2_cluster[\"ra\"] # Cluster Ra in deg\n", + "cluster_dec = dc2_cluster[\"dec\"] # Cluster Dec in deg\n", + "\n", + "source = np.genfromtxt(data_coords_dir + \"/test_source.txt\", names=True)\n", + "\n", + "gal_ra = source[\"ra\"]\n", + "gal_dec = source[\"dec\"]\n", + "gal_e1 = source[\"e1\"]\n", + "gal_e2 = source[\"e2\"]\n", + "gal_z = source[\"zmean\"]\n", + "gal_id = np.arange(len(source[\"ra\"]))" + ] + }, + { + "cell_type": "markdown", + "id": "b8991a69", + "metadata": {}, + "source": [ + "Now we generate a `GCData` object from the imported data. We are going to create two catalogs, each with a different choice of coordinate system, which we specify with `meta={\"coordinate_system\": \"coord\"}`" + ] + }, + { + "cell_type": "code", + "execution_count": 34, + "id": "fcb2ae6c", + "metadata": {}, + "outputs": [], + "source": [ + "dc2_galaxies_euclidean = clmm.GCData(\n", + " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"],\n", + " meta={\"coordinate_system\": \"euclidean\"},\n", + ")\n", + "\n", + "dc2_galaxies_celestial = clmm.GCData(\n", + " [gal_ra, gal_dec, gal_e1, gal_e2, gal_z, gal_id],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"],\n", + " meta={\"coordinate_system\": \"celestial\"},\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2245f754", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a GalaxyCluster.\n", + "dc2_cluster_euclidean = clmm.GalaxyCluster(\n", + " \"Euclidean cluster\",\n", + " cluster_ra,\n", + " cluster_dec,\n", + " cluster_z,\n", + " dc2_galaxies_euclidean,\n", + ")\n", + "\n", + "dc2_cluster_celestial = clmm.GalaxyCluster(\n", + " \"Celestial cluster\",\n", + " cluster_ra,\n", + " cluster_dec,\n", + " cluster_z,\n", + " dc2_galaxies_celestial,\n", + ")\n", + "\n", + "# Convert elipticities into shears for the members.\n", + "dc2_cluster_euclidean.compute_tangential_and_cross_components(add=True)\n", + "dc2_cluster_celestial.compute_tangential_and_cross_components(add=True)\n", + "print(dc2_cluster_euclidean.galcat.colnames)\n", + "print(dc2_cluster_celestial.galcat.colnames)\n", + "\n", + "# Calculate the radial profile of the cluster.\n", + "dc2_cluster_euclidean.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "dc2_cluster_celestial.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "print(dc2_cluster_euclidean.profile.colnames)\n", + "print(dc2_cluster_celestial.profile.colnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1764ff31", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, height_ratios=[3, 1], sharex=True)\n", + "\n", + "ax[0].errorbar(\n", + " dc2_cluster_euclidean.profile[\"radius\"],\n", + " dc2_cluster_euclidean.profile[\"gt\"],\n", + " yerr=dc2_cluster_euclidean.profile[\"gt_err\"],\n", + " alpha=0.5,\n", + " marker=\".\",\n", + " color=\"tab:red\",\n", + " label=\"euclidean\",\n", + ")\n", + "ax[1].errorbar(\n", + " dc2_cluster_euclidean.profile[\"radius\"],\n", + " dc2_cluster_euclidean.profile[\"gx\"],\n", + " yerr=dc2_cluster_euclidean.profile[\"gx_err\"],\n", + " marker=\".\",\n", + " alpha=0.5,\n", + " color=\"tab:red\",\n", + ")\n", + "\n", + "ax[0].errorbar(\n", + " dc2_cluster_celestial.profile[\"radius\"] * 1.02,\n", + " dc2_cluster_celestial.profile[\"gt\"],\n", + " yerr=dc2_cluster_celestial.profile[\"gt_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + " label=\"celestial\",\n", + ")\n", + "ax[1].errorbar(\n", + " dc2_cluster_celestial.profile[\"radius\"] * 1.02,\n", + " dc2_cluster_celestial.profile[\"gx\"],\n", + " yerr=dc2_cluster_celestial.profile[\"gx_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + ")\n", + "\n", + "ax[0].legend()\n", + "ax[0].set_xscale(\"log\")\n", + "ax[1].set_xlabel(\"R [kpc]\")\n", + "ax[0].set_ylabel(\"$g_t$\")\n", + "ax[1].set_ylabel(\"$g_x$\")\n", + "\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "028ab5a1", + "metadata": {}, + "source": [ + "### => There's no signal in either cluster! Is there a problem in DC2?" + ] + }, + { + "cell_type": "markdown", + "id": "8902db07", + "metadata": {}, + "source": [ + "### Example source galaxies from M. Oguri\n", + "\n", + "This dataset is a curated selection of cluster and source catalogs from Summer School lectures delivered by Masamune Oguri. There are eight galaxy clusters in this selection. \n", + "\n", + "More details on the corresponding tutorial can be found at this [GitHub link](https://github.com/oguri/wlcluster_tutorial). These are also in the `euclidean` coordinate system." + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "8eb0581e", + "metadata": {}, + "outputs": [], + "source": [ + "clusters = [\n", + " \"a1703\",\n", + " \"gho1320\",\n", + " \"sdss0851\",\n", + " \"sdss1050\",\n", + " \"sdss1138\",\n", + " \"sdss1226\",\n", + " \"sdss1329\",\n", + " \"sdss1531\",\n", + "]\n", + "\n", + "zl_all = {\n", + " \"a1703\": 0.277,\n", + " \"gho1320\": 0.308,\n", + " \"sdss0851\": 0.370,\n", + " \"sdss1050\": 0.60,\n", + " \"sdss1138\": 0.451,\n", + " \"sdss1226\": 0.435,\n", + " \"sdss1329\": 0.443,\n", + " \"sdss1531\": 0.335,\n", + "}\n", + "\n", + "ra_cl_all = {\n", + " \"a1703\": 198.771833,\n", + " \"gho1320\": 200.703208,\n", + " \"sdss0851\": 132.911917,\n", + " \"sdss1050\": 162.666250,\n", + " \"sdss1138\": 174.537292,\n", + " \"sdss1226\": 186.712958,\n", + " \"sdss1329\": 202.393708,\n", + " \"sdss1531\": 232.794167,\n", + "}\n", + "\n", + "dec_cl_all = {\n", + " \"a1703\": 51.817389,\n", + " \"gho1320\": 31.654944,\n", + " \"sdss0851\": 33.518361,\n", + " \"sdss1050\": 0.285306,\n", + " \"sdss1138\": 27.908528,\n", + " \"sdss1226\": 21.831194,\n", + " \"sdss1329\": 22.721167,\n", + " \"sdss1531\": 34.240278,\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "0f237ea0", + "metadata": {}, + "outputs": [], + "source": [ + "cname = \"a1703\"\n", + "\n", + "# cluster redshift\n", + "zl = zl_all.get(cname)\n", + "\n", + "# coordinates of the cluster center\n", + "ra_cl = ra_cl_all.get(cname)\n", + "dec_cl = dec_cl_all.get(cname)\n", + "\n", + "# fix source redshift to 1.0\n", + "zs = 1.0" + ] + }, + { + "cell_type": "markdown", + "id": "920a2011", + "metadata": {}, + "source": [ + "We inspect the first galaxy cluster, Abell 1703." + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "2a8705d3", + "metadata": {}, + "outputs": [], + "source": [ + "rfile = data_coords_dir + \"/data/shear_\" + cname + \".dat\"\n", + "data = np.loadtxt(rfile, comments=\"#\")\n", + "\n", + "ra = data[:, 0]\n", + "dec = data[:, 1]\n", + "e1 = data[:, 2]\n", + "e2 = data[:, 3]\n", + "wei = data[:, 4]\n", + "ids = np.arange(np.shape(data)[0])\n", + "redshifts = np.ones(np.shape(data)[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 40, + "id": "d76c9464", + "metadata": {}, + "outputs": [], + "source": [ + "oguri_galaxies_euclidean = clmm.GCData(\n", + " [ra, dec, e1, e2, redshifts, ids],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"],\n", + " meta={\"coordinate_system\": \"euclidean\"},\n", + ")\n", + "\n", + "oguri_galaxies_celestial = clmm.GCData(\n", + " [ra, dec, e1, e2, redshifts, ids],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"],\n", + " meta={\"coordinate_system\": \"celestial\"},\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04b71872", + "metadata": {}, + "outputs": [], + "source": [ + "oguri_cluster_euclidean = clmm.GalaxyCluster(cname, ra_cl, dec_cl, zl, oguri_galaxies_euclidean)\n", + "\n", + "oguri_cluster_celestial = clmm.GalaxyCluster(cname, ra_cl, dec_cl, zl, oguri_galaxies_celestial)\n", + "\n", + "# Convert elipticities into shears for the members.\n", + "oguri_cluster_euclidean.compute_tangential_and_cross_components(add=True)\n", + "oguri_cluster_celestial.compute_tangential_and_cross_components(add=True)\n", + "print(oguri_cluster_euclidean.galcat.colnames)\n", + "print(oguri_cluster_celestial.galcat.colnames)\n", + "\n", + "# Calculate the radial profile of the cluster.\n", + "oguri_cluster_euclidean.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "oguri_cluster_celestial.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "print(oguri_cluster_euclidean.profile.colnames)\n", + "print(oguri_cluster_celestial.profile.colnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22bae756", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, height_ratios=[3, 1], sharex=True)\n", + "\n", + "ax[0].errorbar(\n", + " oguri_cluster_euclidean.profile[\"radius\"],\n", + " oguri_cluster_euclidean.profile[\"gt\"],\n", + " yerr=oguri_cluster_euclidean.profile[\"gt_err\"],\n", + " alpha=0.5,\n", + " marker=\".\",\n", + " color=\"tab:red\",\n", + " label=\"euclidean\",\n", + ")\n", + "ax[1].errorbar(\n", + " oguri_cluster_euclidean.profile[\"radius\"],\n", + " oguri_cluster_euclidean.profile[\"gx\"],\n", + " yerr=oguri_cluster_euclidean.profile[\"gx_err\"],\n", + " marker=\".\",\n", + " alpha=0.5,\n", + " color=\"tab:red\",\n", + ")\n", + "\n", + "ax[0].errorbar(\n", + " oguri_cluster_celestial.profile[\"radius\"] * 1.02,\n", + " oguri_cluster_celestial.profile[\"gt\"],\n", + " yerr=oguri_cluster_celestial.profile[\"gt_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + " label=\"celestial\",\n", + ")\n", + "ax[1].errorbar(\n", + " oguri_cluster_celestial.profile[\"radius\"] * 1.02,\n", + " oguri_cluster_celestial.profile[\"gx\"],\n", + " yerr=oguri_cluster_celestial.profile[\"gx_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + ")\n", + "\n", + "ax[0].legend()\n", + "ax[0].set_xscale(\"log\")\n", + "ax[1].set_xlabel(\"R [kpc]\")\n", + "ax[0].set_ylabel(\"$g_t$\")\n", + "ax[1].set_ylabel(\"$g_x$\")\n", + "\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "180b4eb2", + "metadata": {}, + "source": [ + "### => As expected, we only see a clear signal in the `euclidean` cluster" + ] + }, + { + "cell_type": "markdown", + "id": "bd826a37", + "metadata": {}, + "source": [ + "### Example source galaxies from HSC Y3\n", + "\n", + "This dataset is a simplified version of HSC Y3 data (GAMA15H), post-processed by Tomomi Sunayama for testing purposes. The pre-processed data is already public. These catalogs assume a **celestial** coordinate system." + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "cba6862c", + "metadata": {}, + "outputs": [], + "source": [ + "hsc_cluster_cat = np.genfromtxt(\n", + " data_coords_dir + \"/GAMA15H/redmapper_dr8_GAMA15H.txt\",\n", + " dtype=np.dtype(\n", + " [(\"ra\", np.float64), (\"dec\", np.float64), (\"z\", np.float64), (\"richness\", np.float64)]\n", + " ),\n", + ")\n", + "\n", + "hsc_source_cat = fits.getdata(data_coords_dir + \"/GAMA15H/GAMA15H_tutorial.fits\")\n", + "\n", + "cl = hsc_cluster_cat[0]" + ] + }, + { + "cell_type": "markdown", + "id": "bbc58a20", + "metadata": {}, + "source": [ + "Here, we use a KDTree implementation in scipy to extract the background source galaxies for the first galaxy cluster in the dataset." + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "id": "2657871a", + "metadata": {}, + "outputs": [], + "source": [ + "source1 = hsc_source_cat[hsc_source_cat[\"photoz\"] > (cl[\"z\"] + 0.3)]\n", + "tree = spatial.cKDTree(np.array((source1[\"ra\"], source1[\"dec\"])).T)\n", + "sel = tree.query_ball_point([cl[\"ra\"], cl[\"dec\"]], 3)\n", + "bg = source1[sel]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "68670871", + "metadata": {}, + "outputs": [], + "source": [ + "hsc_galaxies_euclidean = clmm.GCData(\n", + " [bg[\"RA\"], bg[\"Dec\"], bg[\"e1\"], bg[\"e2\"], bg[\"photoz\"], bg[\"weight\"]],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"w_ls\"],\n", + " meta={\"coordinate_system\": \"euclidean\"},\n", + ")\n", + "\n", + "hsc_galaxies_celestial = clmm.GCData(\n", + " [bg[\"RA\"], bg[\"Dec\"], bg[\"e1\"], bg[\"e2\"], bg[\"photoz\"], bg[\"weight\"]],\n", + " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"w_ls\"],\n", + " meta={\"coordinate_system\": \"celestial\"},\n", + ")\n", + "\n", + "hsc_cluster_euclidean = clmm.GalaxyCluster(\n", + " \"Eucliden HSC cluster\",\n", + " cl[\"ra\"],\n", + " cl[\"dec\"],\n", + " cl[\"z\"],\n", + " hsc_galaxies_euclidean,\n", + ")\n", + "\n", + "hsc_cluster_celestial = clmm.GalaxyCluster(\n", + " \"Celestial HSC cluster\",\n", + " cl[\"ra\"],\n", + " cl[\"dec\"],\n", + " cl[\"z\"],\n", + " hsc_galaxies_celestial,\n", + ")\n", + "\n", + "# Convert elipticities into shears for the members.\n", + "hsc_cluster_euclidean.compute_tangential_and_cross_components(add=True)\n", + "hsc_cluster_celestial.compute_tangential_and_cross_components(add=True)\n", + "print(hsc_cluster_euclidean.galcat.colnames)\n", + "print(hsc_cluster_celestial.galcat.colnames)\n", + "\n", + "# Calculate the radial profile of the cluster.\n", + "hsc_cluster_euclidean.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "hsc_cluster_celestial.make_radial_profile(\"kpc\", cosmo=cosmo)\n", + "print(hsc_cluster_euclidean.profile.colnames)\n", + "print(hsc_cluster_celestial.profile.colnames)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e2abb47", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, height_ratios=[3, 1], sharex=True)\n", + "\n", + "ax[0].errorbar(\n", + " hsc_cluster_euclidean.profile[\"radius\"],\n", + " hsc_cluster_euclidean.profile[\"gt\"],\n", + " yerr=hsc_cluster_euclidean.profile[\"gt_err\"],\n", + " alpha=0.5,\n", + " marker=\".\",\n", + " color=\"tab:red\",\n", + " label=\"euclidean\",\n", + ")\n", + "ax[1].errorbar(\n", + " hsc_cluster_euclidean.profile[\"radius\"],\n", + " hsc_cluster_euclidean.profile[\"gx\"],\n", + " yerr=hsc_cluster_euclidean.profile[\"gx_err\"],\n", + " marker=\".\",\n", + " alpha=0.5,\n", + " color=\"tab:red\",\n", + ")\n", + "\n", + "ax[0].errorbar(\n", + " hsc_cluster_celestial.profile[\"radius\"] * 1.02,\n", + " hsc_cluster_celestial.profile[\"gt\"],\n", + " yerr=hsc_cluster_celestial.profile[\"gt_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + " label=\"celestial\",\n", + ")\n", + "ax[1].errorbar(\n", + " hsc_cluster_celestial.profile[\"radius\"] * 1.02,\n", + " hsc_cluster_celestial.profile[\"gx\"],\n", + " yerr=hsc_cluster_celestial.profile[\"gx_err\"],\n", + " alpha=0.3,\n", + " marker=\".\",\n", + " color=\"tab:blue\",\n", + ")\n", + "\n", + "ax[0].legend()\n", + "ax[0].set_xscale(\"log\")\n", + "ax[1].set_xlabel(\"R [kpc]\")\n", + "ax[0].set_ylabel(\"$g_t$\")\n", + "ax[1].set_ylabel(\"$g_x$\")\n", + "\n", + "plt.subplots_adjust(hspace=0)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "78fe5a8d", + "metadata": {}, + "source": [ + "### => As expected, we only see a clear signal in the `celestial` cluster" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/examples/demo_coordinate_system_datasets.ipynb b/examples/demo_coordinate_system_datasets.ipynb deleted file mode 100644 index acbedfb65..000000000 --- a/examples/demo_coordinate_system_datasets.ipynb +++ /dev/null @@ -1,1013 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Tests of coordinate system effects on shear profiles\n", - "\n", - "Author: Tomomi Sunayama, July 15, 2024\n", - "\n", - "Tested, modified, and documented by: Camille Avestruz, July 18, 2024\n", - "\n", - "Reviewed by: Caio Lima de Oliveira, July 19, 2024\n", - "\n", - "This notebook illustrates the impact of having an incorrect coordinate system when calculating shear profiles. The default coordinate system for `CLMM` is the euclidean/pixel coordinate system. This is consistent with DC2 catalogs. However, if we input a catalog that assumes a celestial (sky) coordinate system and use the default euclidean (pixel) coordinate system (or vice versa), the signal of shear around a cluster disappears because the signal essentially looks random.\n", - "\n", - "We test:\n", - "- CosmoDC2 source galaxies with shears extracted from `TXPipe` for a single galaxy cluster (euclidean/pixel coordinate system)\n", - "- Example source galaxies for galaxy clusters from a [Summer School](https://github.com/oguri/wlcluster_tutorial) taught by Masamune Oguri (euclidean/pixel coordinate system)\n", - "- HSC Y3 source galaxies with shears post processed by Tomomi Sunayama (celestial/sky coordinate system)\n", - "\n", - "We also \n", - "- Compare the explicit calculation of a shear profile on the HSC Y3 source galaxies against a profile produced from `CLMM`. \n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Instructions to download text data\n", - "\n", - "Before running this notebook, you will need to download some data. The data is available through a [dropbox link](https://www.dropbox.com/scl/fo/dwsccslr5iwb7lnkf8jvx/AJkjgFeemUEHpHaZaHHqpAg?rlkey=efbtsr15mdrs3y6xsm7l48o0r&st=xb58ap0g&dl=0)\n", - "\n", - "First, create a directory where you want to put the example data, e.g. for a given `data_coords_dir`:\n", - "```\n", - "mkdir -p /data_coords\n", - "cd /data_coords\n", - "```\n", - "Download all files from [dropbox link](https://www.dropbox.com/scl/fo/dwsccslr5iwb7lnkf8jvx/AJkjgFeemUEHpHaZaHHqpAg?rlkey=efbtsr15mdrs3y6xsm7l48o0r&st=xb58ap0g&dl=0). This will be a zip file, `CLMM_data.zip` of size 242 Mb. scp or move this to `data_coords_dir`.\n", - "\n", - "From the directory, you should be able to unzip:\n", - "```\n", - "unzip data_CLMM.zip -d .\n", - "```\n", - "You now have the necessary data files to run this notebook. \n", - "\n", - "**Make sure to change the `data_coords_dir` variable in the cell below to the appropriate location where you unzipped these files.**\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# CHANGE TO YOUR LOCATION\n", - "data_coords_dir = \"/data_coords/\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "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\n", - "import pandas\n", - "from astropy.io import fits" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [] - }, - "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)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example galaxy cluster from CosmoDC2\n", - "\n", - "Here, we plot an example galaxy cluster shear profile using `clmm`. The cluster and source galaxy files are generated from the CosmoDC2 processed through TXPipe. We test the coordinate system." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cluster = pandas.read_pickle(data_coords_dir + \"/test_cluster.pkl\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cluster_z = cluster[\"redshift\"] # Cluster redshift\n", - "cluster_ra = cluster[\"ra\"] # Cluster Ra in deg\n", - "cluster_dec = cluster[\"dec\"]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "source = np.genfromtxt(data_coords_dir + \"/test_source.txt\", names=True)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "gal_ra = source[\"ra\"]\n", - "gal_dec = source[\"dec\"]\n", - "gal_e1 = source[\"e1\"]\n", - "gal_e2 = source[\"e2\"]\n", - "gal_z = source[\"zmean\"]\n", - "gal_id = np.arange(len(source[\"ra\"]))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "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", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here we create a `GalaxyCluster` object, specifying an *incorrect* coordinate system. For source galaxies from CosmoDC2, these are in the **euclidean** coordinate system. We use the implemented kwarg when defining the galaxy cluster object to specify the **celestial** coordinate system." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Create a GalaxyCluster.\n", - "cluster = clmm.GalaxyCluster(\n", - " \"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies, coordinate_system=\"celestial\"\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.1, 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": "markdown", - "metadata": {}, - "source": [ - "Here, we plot the resulting profile when `clmm` uses assumes a coordinate system inconsistent with the catalogs provided. You will note that the signal is virtually zero at almost all radii." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "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", - "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", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here we create a `GalaxyCluster` object, specifying the *correct* coordinate system. For source galaxies from CosmoDC2, these are in the **euclidean** coordinate system. We use the implemented kwarg when defining the galaxy cluster object to also specify the **euclidean** coordinate system. However, with a single galaxy cluster, the signal is not significant enough to clearly see a difference. There is a possible excess excess with the correct coordinate system at larger radii. Note: First, the lensing signal in CosmoDC2 clusters at the inner radii is known to be weak due to a limitation in the resolution when the ray tracing was performed in generating the source galaxy shears. Second, this has been process through `TXPipe`, which is something else we will need to understand." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cluster2 = clmm.GalaxyCluster(\n", - " \"Name of cluster\", cluster_ra, cluster_dec, cluster_z, galaxies, coordinate_system=\"euclidean\"\n", - ")\n", - "cluster2.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(cluster2.galcat[\"theta\"], \"radians\", \"Mpc\", cluster2.z, cosmo)\n", - "\n", - "cluster2.make_radial_profile(\n", - " bins=da.make_bins(0.1, 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", - "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", - "ax.errorbar(\n", - " cluster.profile[\"radius\"], cluster.profile[\"gt\"], cluster.profile[\"gt_err\"], label=\"celestial\"\n", - ")\n", - "ax.errorbar(\n", - " cluster2.profile[\"radius\"],\n", - " cluster2.profile[\"gt\"],\n", - " cluster2.profile[\"gt_err\"],\n", - " label=\"euclidean\",\n", - ")\n", - "ax.set_xlabel(\"r [Mpc]\", fontsize=20)\n", - "ax.set_ylabel(r\"$g_t$\", fontsize=20)\n", - "ax.grid(lw=0.3)\n", - "ax.minorticks_on()\n", - "ax.grid(which=\"minor\", lw=0.1)\n", - "plt.legend(fontsize=20)\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example source galaxies from M. Oguri\n", - "\n", - "This dataset is a curated selection of cluster and source catalogs from Summer School lectures delivered by Masamune Oguri. There are eight galaxy clusters in this selection. \n", - "\n", - "More details on the corresponding tutorial can be found at this [GitHub link](https://github.com/oguri/wlcluster_tutorial). These are also in the **euclidean** coordinate system." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "clusters = [\n", - " \"a1703\",\n", - " \"gho1320\",\n", - " \"sdss0851\",\n", - " \"sdss1050\",\n", - " \"sdss1138\",\n", - " \"sdss1226\",\n", - " \"sdss1329\",\n", - " \"sdss1531\",\n", - "]\n", - "\n", - "zl_all = {\n", - " \"a1703\": 0.277,\n", - " \"gho1320\": 0.308,\n", - " \"sdss0851\": 0.370,\n", - " \"sdss1050\": 0.60,\n", - " \"sdss1138\": 0.451,\n", - " \"sdss1226\": 0.435,\n", - " \"sdss1329\": 0.443,\n", - " \"sdss1531\": 0.335,\n", - "}\n", - "\n", - "ra_cl_all = {\n", - " \"a1703\": 198.771833,\n", - " \"gho1320\": 200.703208,\n", - " \"sdss0851\": 132.911917,\n", - " \"sdss1050\": 162.666250,\n", - " \"sdss1138\": 174.537292,\n", - " \"sdss1226\": 186.712958,\n", - " \"sdss1329\": 202.393708,\n", - " \"sdss1531\": 232.794167,\n", - "}\n", - "\n", - "dec_cl_all = {\n", - " \"a1703\": 51.817389,\n", - " \"gho1320\": 31.654944,\n", - " \"sdss0851\": 33.518361,\n", - " \"sdss1050\": 0.285306,\n", - " \"sdss1138\": 27.908528,\n", - " \"sdss1226\": 21.831194,\n", - " \"sdss1329\": 22.721167,\n", - " \"sdss1531\": 34.240278,\n", - "}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cname = \"a1703\"\n", - "\n", - "# cluster redshift\n", - "zl = zl_all.get(cname)\n", - "\n", - "# coordinates of the cluster center\n", - "ra_cl = ra_cl_all.get(cname)\n", - "dec_cl = dec_cl_all.get(cname)\n", - "\n", - "# fix source redshift to 1.0\n", - "zs = 1.0" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We inspect the first galaxy cluster, Abell 1703." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "rfile = data_coords_dir + \"/data/shear_\" + cname + \".dat\"\n", - "data = np.loadtxt(rfile, comments=\"#\")\n", - "\n", - "ra = data[:, 0]\n", - "dec = data[:, 1]\n", - "e1 = data[:, 2]\n", - "e2 = data[:, 3]\n", - "wei = data[:, 4]\n", - "ids = np.arange(np.shape(data)[0])\n", - "redshifts = np.ones(np.shape(data)[0])\n", - "galaxies = clmm.GCData(\n", - " [ra, dec, e1, e2, redshifts, ids], names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"id\"]\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Here we create a GalaxyCluster object, specifying the correct coordinate system. For source galaxies from the Oguri curated dataset, these are in the euclidean coordinate system. We use the implemented kwarg when defining the galaxy cluster object to also **specify the euclidean coordinate system**." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cluster = clmm.GalaxyCluster(\"A1703euc\", ra_cl, dec_cl, zl, galaxies, coordinate_system=\"euclidean\")\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.2, 3.0, 8, 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", - "\n", - "\n", - "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", - "ax.errorbar(\n", - " cluster.profile[\"radius\"], cluster.profile[\"gt\"], cluster.profile[\"gt_err\"], label=\"euclidean\"\n", - ")\n", - "\n", - "# assume incorrect coordinates\n", - "cluster2 = clmm.GalaxyCluster(\n", - " \"A1703cel\", ra_cl, dec_cl, zl, galaxies, coordinate_system=\"celestial\"\n", - ")\n", - "\n", - "cluster2.compute_tangential_and_cross_components()\n", - "print(cluster2.galcat.colnames)\n", - "\n", - "# Measure profile and add profile table to the cluster.\n", - "seps = convert_units(cluster2.galcat[\"theta\"], \"radians\", \"Mpc\", cluster.z, cosmo)\n", - "\n", - "cluster2.make_radial_profile(\n", - " bins=da.make_bins(0.2, 3.0, 8, method=\"evenlog10width\"),\n", - " bin_units=\"Mpc\",\n", - " cosmo=cosmo,\n", - " include_empty_bins=False,\n", - " gal_ids_in_bins=True,\n", - ")\n", - "print(cluster2.profile.colnames)\n", - "\n", - "ax.errorbar(\n", - " cluster2.profile[\"radius\"],\n", - " cluster2.profile[\"gt\"],\n", - " cluster2.profile[\"gt_err\"],\n", - " label=\"celestial\",\n", - ")\n", - "\n", - "\n", - "ax.set_xlabel(\"r [Mpc]\", fontsize=20)\n", - "ax.set_ylabel(r\"$g_t$\", fontsize=20)\n", - "ax.grid(lw=0.3)\n", - "ax.minorticks_on()\n", - "ax.grid(which=\"minor\", lw=0.1)\n", - "plt.legend(fontsize=20)\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Example source galaxies from HSC Y3\n", - "\n", - "This dataset is a simplified version of HSC Y3 data (GAMA15H), post-processed by Tomomi Sunayama for testing purposes. The pre-processed data is already public. These catalogs assume a **celestial** coordinate system." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cluster_cat = np.genfromtxt(\n", - " data_coords_dir + \"/GAMA15H/redmapper_dr8_GAMA15H.txt\",\n", - " dtype=np.dtype(\n", - " [(\"ra\", np.float64), (\"dec\", np.float64), (\"z\", np.float64), (\"richness\", np.float64)]\n", - " ),\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "source_cat = fits.getdata(data_coords_dir + \"/GAMA15H/GAMA15H_tutorial.fits\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "cl = cluster_cat[0]" - ] - }, - { - "cell_type": "markdown", - "id": "8820c46a", - "metadata": {}, - "source": [ - "Here, we use a KDTree implementation in scipy to extract the background source galaxies for the first galaxy cluster in the dataset." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from scipy import spatial\n", - "\n", - "source1 = source_cat[source_cat[\"photoz\"] > (cl[\"z\"] + 0.3)]\n", - "tree = spatial.cKDTree(np.array((source1[\"ra\"], source1[\"dec\"])).T)\n", - "sel = tree.query_ball_point([cl[\"ra\"], cl[\"dec\"]], 3)\n", - "bg = source1[sel]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bf2df4dc", - "metadata": {}, - "outputs": [], - "source": [ - "# Inspect background source galaxy selection\n", - "bg" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a77e1bfe", - "metadata": {}, - "outputs": [], - "source": [ - "sources = clmm.GCData(\n", - " [bg[\"RA\"], bg[\"Dec\"], bg[\"e1\"], bg[\"e2\"], bg[\"photoz\"], bg[\"weight\"]],\n", - " names=[\"ra\", \"dec\", \"e1\", \"e2\", \"z\", \"w_ls\"],\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "45d76a13", - "metadata": {}, - "outputs": [], - "source": [ - "cluster = clmm.GalaxyCluster(\n", - " \"redmapper\", cl[\"ra\"], cl[\"dec\"], cl[\"z\"], sources, coordinate_system=\"celestial\"\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b1f88b26", - "metadata": {}, - "outputs": [], - "source": [ - "sigma_c = cosmo.eval_sigma_crit(cl[\"z\"], sources[\"z\"])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6de9c9de", - "metadata": {}, - "outputs": [], - "source": [ - "cluster.compute_tangential_and_cross_components(\n", - " shape_component1=\"e1\",\n", - " shape_component2=\"e2\",\n", - " tan_component=\"DS_t\",\n", - " cross_component=\"DS_x\",\n", - " cosmo=cosmo,\n", - " is_deltasigma=True,\n", - " use_pdz=False,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5fa13fb2", - "metadata": {}, - "outputs": [], - "source": [ - "cluster" - ] - }, - { - "cell_type": "markdown", - "id": "a4dcd4a7", - "metadata": {}, - "source": [ - "Now we construct a radial profile of the tangential and cross terms for the galaxy cluster" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "seps = convert_units(cluster.galcat[\"theta\"], \"radians\", \"Mpc\", cluster.z, cosmo)\n", - "\n", - "cluster.make_radial_profile(\n", - " tan_component_in=\"DS_t\",\n", - " cross_component_in=\"DS_x\",\n", - " tan_component_out=\"DS_t\",\n", - " cross_component_out=\"DS_x\",\n", - " weights_in=\"w_ls\",\n", - " bins=da.make_bins(0.1, 20.0, 15, method=\"evenlog10width\"),\n", - " bin_units=\"Mpc\",\n", - " cosmo=cosmo,\n", - " include_empty_bins=False,\n", - " gal_ids_in_bins=False,\n", - ")\n", - "print(cluster.profile.colnames)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "286dc0f9", - "metadata": {}, - "outputs": [], - "source": [ - "fig = plt.figure(figsize=(6, 6))\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", - "ax.errorbar(\n", - " cluster.profile[\"radius\"],\n", - " cluster.profile[\"DS_t\"] / 1e13,\n", - " cluster.profile[\"DS_t_err\"] / 1e13,\n", - " label=\"celestial\",\n", - ")\n", - "plt.loglog()\n", - "\n", - "ax.set_xlabel(\"r [Mpc]\", fontsize=20)\n", - "ax.set_ylabel(r\"$\\Delta\\Sigma(r)$\", fontsize=20)\n", - "ax.grid(lw=0.3)\n", - "ax.minorticks_on()\n", - "ax.grid(which=\"minor\", lw=0.1)\n", - "plt.legend(fontsize=20)\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "id": "05e4ee26", - "metadata": {}, - "source": [ - "## Example explicit lensing profile measurement comparison with CLMM profile\n", - "\n", - "Here, we use the example HSC Y3 dataset to explicitly measure the lensing signal without using CLMM for comparison. Note, we need to still define a cosmology to calculate comoving distances." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "651e0048", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "from astropy.io import fits\n", - "import matplotlib as mpl\n", - "import matplotlib.pyplot as plt\n", - "from scipy import spatial\n", - "from astropy.cosmology import WMAP5" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "dbefe49d", - "metadata": {}, - "outputs": [], - "source": [ - "# Read in the data\n", - "cluster_cat = np.genfromtxt(\n", - " data_coords_dir + \"GAMA15H/redmapper_dr8_GAMA15H.txt\",\n", - " dtype=np.dtype(\n", - " [(\"RA\", np.float64), (\"Dec\", np.float64), (\"z\", np.float64), (\"richness\", np.float64)]\n", - " ),\n", - ")\n", - "source_cat = fits.getdata(data_coords_dir + \"GAMA15H/GAMA15H_tutorial.fits\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a826c6e0", - "metadata": {}, - "outputs": [], - "source": [ - "cosmo = WMAP5" - ] - }, - { - "cell_type": "markdown", - "id": "ca4bd885", - "metadata": {}, - "source": [ - "Below, we measure lensing signals with simplified assumptions. We do not account for responsivity, multiplicative, nor additive biases. " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3fdca16e", - "metadata": {}, - "outputs": [], - "source": [ - "# Define functions for the explicit calculation\n", - "\n", - "\n", - "def calcDistanceAngle(a1, d1, a2, d2):\n", - " \"\"\"Compute the angle between lens and source galaxies from ra, dec in radians\n", - " a1 (float, ndarray) : lens RA in radians\n", - " d1 (float, ndarray) : lens DEC in radians\n", - " a2 (float, ndarray) : src RA in radians\n", - " d2 (float, ndarray) : src DEC in radians\n", - " \"\"\"\n", - " return np.arccos(np.cos(d1) * np.cos(d2) * np.cos(a1 - a2) + np.sin(d1) * np.sin(d2))\n", - "\n", - "\n", - "def cosPhi2(a1, a2, d1, d2, distanceAngle):\n", - " \"\"\"Compute $\\cos(2\\phi)$ for a given distance angle between lens and source galaxies from ra, dec in radians\n", - " a1 (float, ndarray) : lens RA in radians\n", - " a2 (float, ndarray) : src RA in radians\n", - " d1 (float, ndarray) : lens DEC in radians\n", - " d2 (float, ndarray) : src DEC in radians\n", - " distanceAngle (float, ndarray) : angular distance between lens and source in radians\n", - " \"\"\"\n", - " return np.sin(a2 - a1) * np.cos(d1) / np.sin(distanceAngle)\n", - "\n", - "\n", - "def sinPhi2(a1, a2, d1, d2, distanceAngle):\n", - " \"\"\"Compute $\\sin(2\\phi)$ for a given distance angle between lens and source galaxies from ra, dec in radians\n", - " a1 (float, ndarray) : lens RA in radians\n", - " a2 (float, ndarray) : src RA in radians\n", - " d1 (float, ndarray) : lens DEC in radians\n", - " d2 (float, ndarray) : src DEC in radians\n", - " distanceAngle (float, ndarray) : angular distance between lens and source in radians\n", - " \"\"\"\n", - " return (-np.cos(d2) * np.sin(d1) + np.sin(d2) * np.cos(d1) * np.cos(a1 - a2)) / np.sin(\n", - " distanceAngle\n", - " )\n", - "\n", - "\n", - "def compute_sin2phi_cos2phi(a1, a2, d1, d2, distanceAngle):\n", - " \"\"\"Compute necessary coefficients for the et and ex components, sin2phi and cos2phi\n", - " a1 (float, ndarray) : lens RA in radians\n", - " a2 (float, ndarray) : src RA in radians\n", - " d1 (float, ndarray) : lens DEC in radians\n", - " d2 (float, ndarray) : src DEC in radians\n", - " distanceAngle (float, ndarray) : angular distance between lens and source in radians\n", - " \"\"\"\n", - " cosp = cosPhi2(a1, a2, d1, d2, distanceAngle)\n", - " sinp = sinPhi2(a1, a2, d1, d2, distanceAngle)\n", - " cos2p = cosp**2 - sinp**2\n", - " sin2p = 2.0 * sinp * cosp\n", - " return cos2p, sin2p\n", - "\n", - "\n", - "def calc_et_ex(e1, e2, cos2p, sin2p):\n", - " \"\"\"Calculate the et and ex from the e1 e2 values of all sources and their sin2phi, cos2phi\"\"\"\n", - " et = -(e1 * cos2p + e2 * sin2p)\n", - " ex = -(-e1 * sin2p + e2 * cos2p)\n", - " return et, ex\n", - "\n", - "\n", - "def populate_profile_sums(ps, i_r, src_in_bin, Sigma_cr, sel, et, ex):\n", - " \"\"\"Populate the profile sums at a given radian bin from the calculated selection, sigma_crit, et, and ex\"\"\"\n", - " ps[\"n\"][i_r] += sel.sum()\n", - " ps[\"e_sq\"][i_r] += np.sum(et**2 + ex**2)\n", - "\n", - " wt = src_in_bin[\"weight\"] * Sigma_cr**-2 # multiply by the lens weights if it is not one\n", - " ps[\"w\"][i_r] += np.sum(wt)\n", - "\n", - " wetsigma = wt * Sigma_cr * et\n", - " ps[\"wetsigma\"][i_r] += np.sum(wetsigma)\n", - " ps[\"wetsigma_sq\"][i_r] += np.sum(wetsigma**2)\n", - "\n", - " wexsigma = wt * Sigma_cr * ex\n", - " ps[\"wexsigma\"][i_r] += np.sum(wexsigma)\n", - " ps[\"wexsigma_sq\"][i_r] += np.sum(wexsigma**2)\n", - "\n", - " wsigmainv = wt * 1.0 / Sigma_cr\n", - " ps[\"wsigmainv\"][i_r] += np.sum(wsigmainv)\n", - "\n", - " wzs = wt * src_in_bin[\"photoz\"]\n", - " ps[\"wzs\"][i_r] += np.sum(wzs)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f343ae7c", - "metadata": {}, - "outputs": [], - "source": [ - "# Relevant constants, radial binning, source photoz range and lens photoz range\n", - "\n", - "d2r = np.pi / 180.0\n", - "r2d = 180.0 / np.pi\n", - "Mpc = 3.08568025 * 10**19 # 1Mpc = 3.08568025*10**19 km\n", - "M_sun = 1.9884 * 10**33 # solar mass [g]\n", - "c_light = 2.99792458 * 10**5 # speed of light [km/s]\n", - "G = 6.67384 * 10 ** (-20) # gravitational constant [km^3s^-2kg^-1]\n", - "Sigma_cr_fact = c_light**2 / (4 * np.pi * G) * Mpc * 10**3 / M_sun\n", - "rbin_edges = np.logspace(-1, np.log10(20), 15) # Define your radial bins\n", - "\n", - "# Named numpy arrays for relevant profile values to explicitly compute and sum at each radii\n", - "profile_names = [\n", - " \"e_sq\",\n", - " \"w\",\n", - " \"wetsigma\",\n", - " \"wetsigma_sq\",\n", - " \"wexsigma\",\n", - " \"wexsigma_sq\",\n", - " \"wsigmainv\",\n", - " \"wzs\",\n", - " \"n\",\n", - "]\n", - "profile_sums = np.rec.fromarrays(\n", - " [np.zeros(len(rbin_edges) - 1) for i in profile_names], names=profile_names\n", - ")\n", - "\n", - "source_pz = {\"min\": 0.5, \"max\": 10}\n", - "lens_pz = {\"min\": 0.1, \"max\": 0.2}" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8da9ec37", - "metadata": {}, - "outputs": [], - "source": [ - "# Select lens clusters and source galaxies from catalogs using kdtree\n", - "\n", - "source_pz_criteria = (source_cat[\"photoz\"] < source_pz[\"max\"]) & (\n", - " source_cat[\"photoz\"] > source_pz[\"min\"]\n", - ")\n", - "selected_sources = source_cat[source_pz_criteria]\n", - "\n", - "tree = spatial.cKDTree(np.array((selected_sources[\"RA\"], selected_sources[\"Dec\"])).T)\n", - "\n", - "# We only select one,selecting many will take much more time to compute.\n", - "lens_pz_criteria = (cluster_cat[\"z\"] > lens_pz[\"min\"]) & (cluster_cat[\"z\"] < lens_pz[\"max\"])\n", - "lens_clusters = cluster_cat[lens_pz_criteria][:1]\n", - "\n", - "# Set weights for the cluster lenses to one\n", - "lens_weights = np.ones(lens_clusters.size)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "972bdfcb", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "# Calculate lensing profiles for each cluster lens\n", - "\n", - "for ilens in np.arange(lens_clusters.size):\n", - "\n", - " # Select source galaxies for this cluster lens\n", - " sel = tree.query_ball_point([lens_clusters[\"RA\"][ilens], lens_clusters[\"Dec\"][ilens]], 3)\n", - " sel_z = (\n", - " source_cat[sel][\"photoz\"] > lens_clusters[\"z\"][ilens]\n", - " ) # Try to change the source galaxy selection\n", - " source_bg = source_cat[sel][sel_z]\n", - "\n", - " # Compute an angle between the lens and the source\n", - " theta_ls = calcDistanceAngle(\n", - " lens_clusters[\"RA\"][ilens] * d2r,\n", - " lens_clusters[\"Dec\"][ilens] * d2r,\n", - " source_bg[\"RA\"] * d2r,\n", - " source_bg[\"Dec\"] * d2r,\n", - " )\n", - "\n", - " # Compute the comoving distance of the lens\n", - " l_chi = cosmo.comoving_distance((lens_clusters[\"z\"][ilens])).value\n", - " r = theta_ls * l_chi\n", - " assign_r = np.digitize(r, rbin_edges)\n", - "\n", - " for i_r in range(len(rbin_edges) - 1):\n", - " # Subselection mask of source galaxies in the radial bin\n", - " sel = assign_r == i_r + 1\n", - "\n", - " # Subselected source galaxies and their respective angle, theta, to lens\n", - " source_bg_inbin = source_bg[sel]\n", - " theta_sub = theta_ls[sel]\n", - "\n", - " # Compute the cos(2*phi) and sin(2*phi) for a given distance angle between lens and source galaxies\n", - " cos2p, sin2p = compute_sin2phi_cos2phi(\n", - " lens_clusters[\"RA\"][ilens] * d2r,\n", - " source_bg_inbin[\"RA\"] * d2r,\n", - " lens_clusters[\"Dec\"][ilens] * d2r,\n", - " source_bg_inbin[\"Dec\"] * d2r,\n", - " theta_sub,\n", - " )\n", - "\n", - " # Calculate tangential and cross terms from e1, e2 of all source galaxies in the rbin\n", - " et, ex = calc_et_ex(source_bg_inbin[\"e1\"], source_bg_inbin[\"e2\"], cos2p, sin2p)\n", - "\n", - " # Calculate critical surface mass density [M_sun/ comoving Mpc^2]. (1+zl)**-2 is for comoving coordinates.\n", - " comoving_distance = cosmo.comoving_distance((source_bg_inbin[\"photoz\"])).value\n", - " Sigma_cr = (\n", - " Sigma_cr_fact\n", - " / (1.0 - l_chi / comoving_distance)\n", - " / l_chi\n", - " / (1.0 + lens_clusters[\"z\"][ilens])\n", - " / 10**12\n", - " )\n", - "\n", - " # Populate the profile_sums at this radial bin for this cluster lens\n", - " populate_profile_sums(profile_sums, i_r, source_bg_inbin, Sigma_cr, sel, et, ex)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8078ef52", - "metadata": {}, - "outputs": [], - "source": [ - "# Compute the lensing signal and errors to plot\n", - "\n", - "radial_bin = (\n", - " 2.0\n", - " * (rbin_edges[1:] ** 3 - rbin_edges[:-1] ** 3)\n", - " / (3.0 * (rbin_edges[1:] ** 2 - rbin_edges[:-1] ** 2))\n", - ")\n", - "gt = 0.5 * profile_sums[\"wetsigma\"] / profile_sums[\"w\"]\n", - "gt_err = 0.5 * np.sqrt(profile_sums[\"wetsigma_sq\"]) / profile_sums[\"w\"]\n", - "gx = 0.5 * profile_sums[\"wexsigma\"] / profile_sums[\"w\"]\n", - "gx_err = 0.5 * np.sqrt(profile_sums[\"wexsigma_sq\"]) / profile_sums[\"w\"]\n", - "sigma_cr = 1.0 / (profile_sums[\"wsigmainv\"] / profile_sums[\"w\"])" - ] - }, - { - "cell_type": "markdown", - "id": "4f4756df", - "metadata": {}, - "source": [ - "Below, we compare the explicitly calculated lensing signal against the CLMM calculated signal. You will notice that the `CLMM` calculated profile is systematically lower than the one calculated using Tomomi's code. This is likely due to a combination of assumed weighting scheme and other factors that differ between HSC post-processing and what `CLMM` assumes or a \"little h\" problem, which we will need to understand and possibly address." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Figure for the lensing signal\n", - "plt.figure(figsize=(6, 6))\n", - "\n", - "# From explicit calculation\n", - "plt.errorbar(radial_bin, gt, yerr=gt_err, label=\"original\")\n", - "\n", - "# From CLMM\n", - "plt.errorbar(\n", - " cluster.profile[\"radius\"],\n", - " cluster.profile[\"DS_t\"] / 1e13,\n", - " cluster.profile[\"DS_t_err\"] / 1e13,\n", - " label=\"CLMM\",\n", - ")\n", - "plt.loglog()\n", - "plt.legend(fontsize=20)\n", - "plt.xlabel(r\"$R[h^{-1}{\\rm Mpc}]$\", fontsize=20)\n", - "plt.ylabel(r\"$\\Delta\\Sigma(R)$\", fontsize=20)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.6.9" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/tests/test_galaxycluster.py b/tests/test_galaxycluster.py index 874f2453d..732f7ae35 100644 --- a/tests/test_galaxycluster.py +++ b/tests/test_galaxycluster.py @@ -24,7 +24,6 @@ def test_initialization(): "dec": 34.0, "z": 0.3, "galcat": GCData(), - "coordinate_system": "euclidean", } cl1 = clmm.GalaxyCluster(**testdict1) @@ -33,7 +32,6 @@ def test_initialization(): assert_equal(testdict1["dec"], cl1.dec) assert_equal(testdict1["z"], cl1.z) assert isinstance(cl1.galcat, GCData) - assert_equal(testdict1["coordinate_system"], cl1.coordinate_system) def test_integrity(): # Converge on name @@ -53,7 +51,6 @@ def test_integrity(): # Converge on name dec=34.0, z=0.3, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( ValueError, @@ -63,7 +60,6 @@ def test_integrity(): # Converge on name dec=34.0, z=0.3, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( ValueError, @@ -73,7 +69,6 @@ def test_integrity(): # Converge on name dec=95.0, z=0.3, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( ValueError, @@ -83,7 +78,6 @@ def test_integrity(): # Converge on name dec=-95.0, z=0.3, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( ValueError, @@ -93,17 +87,6 @@ def test_integrity(): # Converge on name dec=34.0, z=-0.3, galcat=GCData(), - coordinate_system="euclidean", - ) - assert_raises( - ValueError, - clmm.GalaxyCluster, - unique_id=1, - ra=161.3, - dec=34.0, - z=0.3, - galcat=GCData(), - coordinate_system="blah", ) # Test that inputs are the correct type @@ -115,7 +98,6 @@ def test_integrity(): # Converge on name dec=34.0, z=0.3, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( TypeError, @@ -125,7 +107,6 @@ def test_integrity(): # Converge on name dec=34.0, z=0.3, galcat=1, - coordinate_system="euclidean", ) assert_raises( TypeError, @@ -135,7 +116,6 @@ def test_integrity(): # Converge on name dec=34.0, z=0.3, galcat=[], - coordinate_system="euclidean", ) assert_raises( TypeError, @@ -145,7 +125,6 @@ def test_integrity(): # Converge on name dec=34.0, z=0.3, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( TypeError, @@ -155,7 +134,6 @@ def test_integrity(): # Converge on name dec=None, z=0.3, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( TypeError, @@ -165,7 +143,6 @@ def test_integrity(): # Converge on name dec=34.0, z=None, galcat=GCData(), - coordinate_system="euclidean", ) assert_raises( TypeError, @@ -175,20 +152,15 @@ def test_integrity(): # Converge on name dec=34.0, z=0.3, galcat=None, - coordinate_system=2, ) # Test that id can support numbers and strings assert isinstance( - clmm.GalaxyCluster( - unique_id=1, ra=161.3, dec=34.0, z=0.3, galcat=GCData(), coordinate_system="euclidean" - ).unique_id, + clmm.GalaxyCluster(unique_id=1, ra=161.3, dec=34.0, z=0.3, galcat=GCData()).unique_id, str, ) assert isinstance( - clmm.GalaxyCluster( - unique_id="1", ra=161.3, dec=34.0, z=0.3, galcat=GCData(), coordinate_system="euclidean" - ).unique_id, + clmm.GalaxyCluster(unique_id="1", ra=161.3, dec=34.0, z=0.3, galcat=GCData()).unique_id, str, ) @@ -515,9 +487,7 @@ def test_coordinate_system(): shear1 = [0.2, 0.4] shear2_euclidean = [0.3, 0.5] shear2_celestial = [-0.3, -0.5] - # Set up radial values - bins_radians = [0.002, 0.003, 0.004] - bin_units = "radians" + # create cluster cl_euclidean = clmm.GalaxyCluster( unique_id="test", @@ -527,8 +497,8 @@ def test_coordinate_system(): galcat=GCData( [ra_source, dec_source, shear1, shear2_euclidean, z_src], names=("ra", "dec", "e1", "e2", "z"), + meta={"coordinate_system": "euclidean"}, ), - coordinate_system="euclidean", ) cl_celestial = clmm.GalaxyCluster( unique_id="test", @@ -538,22 +508,39 @@ def test_coordinate_system(): galcat=GCData( [ra_source, dec_source, shear1, shear2_celestial, z_src], names=("ra", "dec", "e1", "e2", "z"), + meta={"coordinate_system": "celestial"}, ), - coordinate_system="celestial", ) cl_euclidean.compute_tangential_and_cross_components() cl_celestial.compute_tangential_and_cross_components() + # assert tangential and cross components are computed consistently assert_allclose( cl_euclidean.galcat["et"], cl_celestial.galcat["et"], **TOLERANCE, - err_msg="Tangential component conversion between ellipticity coordinate systems failed", + err_msg="Tangential component computation different between coordinate systems", ) assert_allclose( cl_euclidean.galcat["ex"], -cl_celestial.galcat["ex"], **TOLERANCE, - err_msg="Cross component conversion between ellipticity coordinate systems failed", + err_msg="Cross component computation different between coordinate systems", + ) + + cl_euclidean.update_coordinate_system("celestial", "e2", "ex") + + # assert components are updated correctly + assert_allclose( + cl_euclidean.galcat["e2"], + cl_celestial.galcat["e2"], + **TOLERANCE, + err_msg="e2 column update failed", + ) + assert_allclose( + cl_euclidean.galcat["ex"], + cl_celestial.galcat["ex"], + **TOLERANCE, + err_msg="ex column update failed", ) diff --git a/tests/test_gcdata.py b/tests/test_gcdata.py index dbe2a0a84..18c6f35c5 100644 --- a/tests/test_gcdata.py +++ b/tests/test_gcdata.py @@ -1,7 +1,8 @@ """ Tests for datatype and galaxycluster """ -from numpy.testing import assert_raises, assert_equal + +from numpy.testing import assert_raises, assert_equal, assert_warns, assert_no_warnings from clmm import GCData from clmm import Cosmology @@ -11,6 +12,22 @@ def test_init(): """test init""" gcdata = GCData() assert_equal(None, gcdata.meta["cosmo"]) + assert_equal("euclidean", gcdata.meta["coordinate_system"]) + + # assert that default coordinate_system warning is raised + assert_warns(UserWarning, GCData) + + # assert that no warning is raised when coordinate_system is set + assert_no_warnings(GCData, meta={"coordinate_system": "celestial"}) + assert_no_warnings(GCData, copy_indices=False) + assert_no_warnings(GCData, masked=False) + assert_no_warnings(GCData, copy=False) + + # assert errros are raised when coordinate_system is not valid + assert_raises(ValueError, GCData, meta={"coordinate_system": "other"}) + assert_raises(TypeError, GCData, meta={"coordinate_system": 2}) + assert_raises(TypeError, GCData, meta={"coordinate_system": None}) + assert_raises(TypeError, GCData, meta={"coordinate_system": ["celestial", 2]}) def test_update_cosmo(): @@ -164,3 +181,68 @@ def test_pzfuncs(): assert isinstance(gcdata._repr_html_(), str) assert_raises(NotImplementedError, gcdata.has_pzpdfs) assert_raises(NotImplementedError, gcdata.get_pzpdfs) + + +def test_update_coordinate_system(): + ngals = 5 + e1 = [0.1] * ngals + e2 = [0.2] * ngals + ex = [0.3] * ngals + + # raise error if not updated from update_coordinate_system + gcdata = GCData(meta={"coordinate_system": "celestial"}) + gcdata["e1"] = e1 + gcdata["e2"] = e2 + gcdata["ex"] = ex + assert_raises(ValueError, gcdata.meta.__setitem__, "coordinate_system", "euclidean") + + # raise error if column is not present + gcdata = GCData(meta={"coordinate_system": "euclidean"}) + gcdata["e1"] = e1 + gcdata["e2"] = e2 + gcdata["ex"] = ex + assert_raises(ValueError, gcdata.update_coordinate_system, "celestial", "et") + + # update coordinate system with bogus system + gcdata = GCData(meta={"coordinate_system": "celestial"}) + gcdata["e1"] = e1 + gcdata["e2"] = e2 + gcdata["ex"] = ex + assert_raises(ValueError, gcdata.update_coordinate_system, "other") + assert_raises(TypeError, gcdata.update_coordinate_system, 2) + assert_raises(TypeError, gcdata.update_coordinate_system, None) + assert_raises(TypeError, gcdata.update_coordinate_system, ["celestial", 2]) + + # update coordinate system without conversion + gcdata = GCData(meta={"coordinate_system": "euclidean"}) + gcdata["e1"] = e1 + gcdata["e2"] = e2 + gcdata["ex"] = ex + assert_warns(UserWarning, gcdata.update_coordinate_system, "celestial") + assert_equal("celestial", gcdata.meta["coordinate_system"]) + assert_equal(e1, gcdata["e1"]) + assert_equal(e2, gcdata["e2"]) + assert_equal(ex, gcdata["ex"]) + + # update coordinate system to the same system + gcdata = GCData(meta={"coordinate_system": "celestial"}) + gcdata["e1"] = e1 + gcdata["e2"] = e2 + gcdata["ex"] = ex + assert_no_warnings(gcdata.update_coordinate_system, "celestial") + assert_equal("celestial", gcdata.meta["coordinate_system"]) + assert_no_warnings(gcdata.update_coordinate_system, "celestial", "e2", "ex") + assert_equal(e1, gcdata["e1"]) + assert_equal(e2, gcdata["e2"]) + assert_equal(ex, gcdata["ex"]) + + # update coordinate system with conversion + gcdata = GCData(meta={"coordinate_system": "celestial"}) + gcdata["e1"] = e1 + gcdata["e2"] = e2 + gcdata["ex"] = ex + assert_warns(UserWarning, gcdata.update_coordinate_system, "euclidean", "e2", "ex") + assert_equal("euclidean", gcdata.meta["coordinate_system"]) + assert_equal(e1, gcdata["e1"]) + assert_equal(e2, -gcdata["e2"]) + assert_equal(ex, -gcdata["ex"]) diff --git a/tests/test_mockdata.py b/tests/test_mockdata.py index edb202410..8078dc5aa 100644 --- a/tests/test_mockdata.py +++ b/tests/test_mockdata.py @@ -34,7 +34,9 @@ def test_mock_data(): warnings.simplefilter("always") # Trigger a warning. np.random.seed(314) - mock.generate_galaxy_catalog(1e15, 0.3, 4, cosmo, 0.30001, ngals=1000, nretry=0) + mock.generate_galaxy_catalog( + 1e15, 0.3, 4, cosmo, 0.30001, ngals=1000, nretry=0, coordinate_system="euclidean" + ) # Verify some things assert len(warn) == 1 # Simple test to check if option with ngal_density is working @@ -257,17 +259,13 @@ def test_shapenoise(): # Verify that the shape noise is Gaussian around 0 (for the very small shear here) sigma = 0.25 - data = mock.generate_galaxy_catalog( - 10**12.0, 0.3, 4, cosmo, 0.8, ngals=50000, shapenoise=sigma - ) + data = mock.generate_galaxy_catalog(10**12.0, 0.3, 4, cosmo, 0.8, ngals=50000, shapenoise=sigma) # Check that there are no galaxies with |e|>1 assert_equal(np.count_nonzero((data["e1"] > 1) | (data["e1"] < -1)), 0) assert_equal(np.count_nonzero((data["e2"] > 1) | (data["e2"] < -1)), 0) # Check that shape noise is Guassian with correct std dev bins = np.arange(-1, 1.1, 0.1) - gauss = ( - 5000 * np.exp(-0.5 * (bins[:-1] + 0.05) ** 2 / sigma**2) / (sigma * np.sqrt(2 * np.pi)) - ) + gauss = 5000 * np.exp(-0.5 * (bins[:-1] + 0.05) ** 2 / sigma**2) / (sigma * np.sqrt(2 * np.pi)) assert_allclose(np.histogram(data["e1"], bins=bins)[0], gauss, atol=50, rtol=0.05) assert_allclose(np.histogram(data["e2"], bins=bins)[0], gauss, atol=50, rtol=0.05)