diff --git a/clmm/__init__.py b/clmm/__init__.py index 0111be973..d998efc18 100644 --- a/clmm/__init__.py +++ b/clmm/__init__.py @@ -26,4 +26,4 @@ ) from . import support -__version__ = "1.12.5" +__version__ = "1.13.0" diff --git a/clmm/clusterensemble.py b/clmm/clusterensemble.py index 6a518279c..787a1de2e 100644 --- a/clmm/clusterensemble.py +++ b/clmm/clusterensemble.py @@ -7,6 +7,7 @@ from .gcdata import GCData from .dataops import make_stacked_radial_profile +from .utils import DiffArray class ClusterEnsemble: @@ -27,7 +28,7 @@ class ClusterEnsemble: * "tan_sc" : tangential component computed with sample covariance * "cross_sc" : cross component computed with sample covariance - * "tan_jk" : tangential component computed with bootstrap + * "tan_bs" : tangential component computed with bootstrap * "cross_bs" : cross component computed with bootstrap * "tan_jk" : tangential component computed with jackknife * "cross_jk" : cross component computed with jackknife @@ -46,7 +47,7 @@ def __init__(self, unique_id, gc_list=None, **kwargs): else: raise TypeError(f"unique_id incorrect type: {type(unique_id)}") self.unique_id = unique_id - self.data = GCData(meta={"bin_units": None}) + self.data = GCData(meta={"bin_units": None, "radius_min": None, "radius_max": None}) if gc_list is not None: self._add_values(gc_list, **kwargs) self.stacked_data = None @@ -198,6 +199,9 @@ def add_individual_radial_profile( """ cl_bin_units = profile_table.meta.get("bin_units", None) self.data.update_info_ext_valid("bin_units", self.data, cl_bin_units, overwrite=False) + for col in ("radius_min", "radius_max"): + value = DiffArray(profile_table[col]) + self.data.update_info_ext_valid(col, self.data, value, overwrite=False) cl_cosmo = profile_table.meta.get("cosmo", None) self.data.update_info_ext_valid("cosmo", self.data, cl_cosmo, overwrite=False) @@ -248,9 +252,14 @@ def make_stacked_radial_profile(self, tan_component="gt", cross_component="gx", [self.data[tan_component], self.data[cross_component]], ) self.stacked_data = GCData( - [radius, *components], - meta=self.data.meta, - names=("radius", tan_component, cross_component), + [ + self.data.meta["radius_min"].value, + self.data.meta["radius_max"].value, + radius, + *components, + ], + meta={k: v for k, v in self.data.meta.items() if k not in ("radius_min", "radius_max")}, + names=("radius_min", "radius_max", "radius", tan_component, cross_component), ) def compute_sample_covariance(self, tan_component="gt", cross_component="gx"): diff --git a/clmm/dataops/__init__.py b/clmm/dataops/__init__.py index 262c31dfa..eab62c280 100644 --- a/clmm/dataops/__init__.py +++ b/clmm/dataops/__init__.py @@ -1,8 +1,6 @@ -"""Functions to compute polar/azimuthal averages in radial bins""" - +"""Data operation for polar/azimuthal averages in radial bins and weights""" import warnings import numpy as np -import scipy from astropy.coordinates import SkyCoord from astropy import units as u from ..gcdata import GCData @@ -17,11 +15,7 @@ _validate_is_deltasigma_sigma_c, _validate_coordinate_system, ) -from ..redshift import ( - _integ_pzfuncs, - compute_for_good_redshifts, -) -from ..theory import compute_critical_surface_density_eff +from ..redshift import _integ_pzfuncs def compute_tangential_and_cross_components( diff --git a/clmm/gcdata.py b/clmm/gcdata.py index eae2c102c..29f9eb5d5 100644 --- a/clmm/gcdata.py +++ b/clmm/gcdata.py @@ -150,7 +150,7 @@ def update_info_ext_valid(self, key, gcdata, ext_value, overwrite=False): key: str Name of key to compare and update. gcdata: GCData - Table to check if same cosmology. + Table to check if same cosmology and ensemble bins. ext_value: Value to be compared to. overwrite: bool diff --git a/clmm/utils/__init__.py b/clmm/utils/__init__.py index 3bbf9b81c..bcd1806ac 100644 --- a/clmm/utils/__init__.py +++ b/clmm/utils/__init__.py @@ -41,6 +41,7 @@ _validate_dec, _validate_is_deltasigma_sigma_c, _validate_coordinate_system, + DiffArray, ) from .units import ( diff --git a/clmm/utils/validation.py b/clmm/utils/validation.py index 767ef67f5..1002ed222 100644 --- a/clmm/utils/validation.py +++ b/clmm/utils/validation.py @@ -226,7 +226,6 @@ 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): r"""Validate the coordinate system. @@ -245,3 +244,37 @@ def _validate_coordinate_system(loc, coordinate_system, valid_type): 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'.") + +class DiffArray: + """Array where arr1==arr2 is actually all(arr1==arr)""" + + def __init__(self, array): + self.value = np.array(array) + + def __eq__(self, other): + # pylint: disable=unidiomatic-typecheck + if type(other) != type(self): + return False + if self.value.size != other.value.size: + return False + return (self.value == other.value).all() + + def __repr__(self): + out = str(self.value) + if self.value.size > 4: + out = self._get_lim_str(out) + " ... " + self._get_lim_str(out[::-1])[::-1] + return out + + def _get_lim_str(self, out): + # pylint: disable=undefined-loop-variable + # get count starting point + for init_index, char in enumerate(out): + if all(char != _char for _char in "[]() "): + break + # get str + sep = 0 + for i, char in enumerate(out[init_index + 1 :]): + sep += int(char == " " and out[i + init_index] != " ") + if sep == 2: + break + return out[: i + init_index + 1] diff --git a/examples/demo_mock_ensemble.ipynb b/examples/demo_mock_ensemble.ipynb index fd3c734db..ed0233561 100644 --- a/examples/demo_mock_ensemble.ipynb +++ b/examples/demo_mock_ensemble.ipynb @@ -320,6 +320,43 @@ }, { "cell_type": "markdown", + "id": "84bd786f", + "metadata": {}, + "source": [ + "The individual cluster data and profiles are stored at the `.data` table of the `ClusterEnsemble`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "864f3acb", + "metadata": {}, + "outputs": [], + "source": [ + "clusterensemble.data[:3]" + ] + }, + { + "cell_type": "markdown", + "id": "6ea533b0", + "metadata": {}, + "source": [ + "The edges of the radial bins, their units, and the cosmology are stored on the metadata of this table:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39723fb1", + "metadata": {}, + "outputs": [], + "source": [ + "clusterensemble.data.meta" + ] + }, + { + "cell_type": "markdown", + "id": "99e3fe18", "metadata": {}, "source": [ "## Stacked profile of the cluster ensemble\n", @@ -335,6 +372,16 @@ "clusterensemble.make_stacked_radial_profile(tan_component=\"DS_t\", cross_component=\"DS_x\")" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "98b9fa63", + "metadata": {}, + "outputs": [], + "source": [ + "clusterensemble.stacked_data" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -601,7 +648,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.6.9" + "version": "3.10.8" } }, "nbformat": 4, diff --git a/examples/demo_mock_ensemble_realistic.ipynb b/examples/demo_mock_ensemble_realistic.ipynb index 42cae6d7d..b251ab103 100644 --- a/examples/demo_mock_ensemble_realistic.ipynb +++ b/examples/demo_mock_ensemble_realistic.ipynb @@ -175,7 +175,7 @@ "source": [ "Ncm.cfg_init()\n", "# cosmo_nc = Nc.HICosmoDEXcdm()\n", - "cosmo_nc = Nc.HICosmo.new_from_name(Nc.HICosmo, \"NcHICosmoDECpl{'massnu-length':<1>}\")\n", + "cosmo_nc = Nc.HICosmoDECpl(massnu_length=1)\n", "cosmo_nc.omega_x2omega_k()\n", "cosmo_nc.param_set_by_name(\"w0\", -1.0)\n", "cosmo_nc.param_set_by_name(\"w1\", 0.0)\n", @@ -201,7 +201,7 @@ "\n", "dist = Nc.Distance.new(2.0)\n", "dist.prepare_if_needed(cosmo_nc)\n", - "tf = Nc.TransferFunc.new_from_name(\"NcTransferFuncEH\")\n", + "tf = Nc.TransferFuncEH()\n", "\n", "psml = Nc.PowspecMLTransfer.new(tf)\n", "psml.require_kmin(1.0e-6)\n", @@ -903,6 +903,42 @@ " )" ] }, + { + "cell_type": "markdown", + "id": "b5e17ee0", + "metadata": {}, + "source": [ + "The individual cluster data and profiles are stored at the `.data` table of the `ClusterEnsemble`:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a970edb", + "metadata": {}, + "outputs": [], + "source": [ + "clusterensemble.data[:3]" + ] + }, + { + "cell_type": "markdown", + "id": "7320e899", + "metadata": {}, + "source": [ + "The edges of the radial bins, their units, and the cosmology are stored on the metadata of this table:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9ea8436", + "metadata": {}, + "outputs": [], + "source": [ + "clusterensemble.data.meta" + ] + }, { "cell_type": "markdown", "id": "9091de7c", @@ -922,6 +958,16 @@ "clusterensemble.make_stacked_radial_profile(tan_component=\"DS_t\", cross_component=\"DS_x\")" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4125570", + "metadata": {}, + "outputs": [], + "source": [ + "clusterensemble.stacked_data" + ] + }, { "cell_type": "markdown", "id": "771c1382", @@ -1193,7 +1239,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.0" + "version": "3.11.9" } }, "nbformat": 4, diff --git a/tests/test_dataops.py b/tests/test_dataops.py index 54942c9dd..679f56917 100644 --- a/tests/test_dataops.py +++ b/tests/test_dataops.py @@ -145,6 +145,13 @@ def test_compute_lensing_angles_flatsky(): ra_l, dec_l, ra_s, dec_s, coordinate_system="celestial" ) + assert_allclose( + da._compute_lensing_angles_flatsky(-180, dec_l, np.array([180.1, 179.7]), dec_s), + [[0.0012916551296819666, 0.003424250083245557], [-2.570568636904587, 0.31079754672944354]], + TOLERANCE["rtol"], + err_msg="Failure when ra_l and ra_s are the same but one is defined negative", + ) + assert_allclose( thetas_celestial, thetas_euclidean, diff --git a/tests/test_mockdata.py b/tests/test_mockdata.py index 9cd53e219..4b2943d6a 100644 --- a/tests/test_mockdata.py +++ b/tests/test_mockdata.py @@ -264,13 +264,17 @@ 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) diff --git a/tests/test_theory.py b/tests/test_theory.py index a17ab661d..6a5a7f96d 100644 --- a/tests/test_theory.py +++ b/tests/test_theory.py @@ -7,7 +7,11 @@ from clmm.constants import Constants as clc from clmm.galaxycluster import GalaxyCluster from clmm import GCData -from clmm.utils import compute_beta_s_square_mean_from_distribution, compute_beta_s_mean_from_distribution, compute_beta_s_func +from clmm.utils import ( + compute_beta_s_square_mean_from_distribution, + compute_beta_s_mean_from_distribution, + compute_beta_s_func, +) from clmm.redshift.distributions import chang2013, desc_srd TOLERANCE = {"rtol": 1.0e-8} @@ -213,7 +217,7 @@ def test_compute_reduced_shear(modeling_data): assert_allclose( theo.compute_reduced_shear_from_convergence(np.array(shear), np.array(convergence)), np.array(truth), - **TOLERANCE + **TOLERANCE, ) @@ -254,7 +258,7 @@ def helper_profiles(func): assert_allclose( func(r3d, mdelta, cdelta, z_cl, cclcosmo, halo_profile_model="nfw"), defaulttruth, - **TOLERANCE + **TOLERANCE, ) assert_allclose( func(r3d, mdelta, cdelta, z_cl, cclcosmo, massdef="mean"), defaulttruth, **TOLERANCE @@ -263,7 +267,7 @@ def helper_profiles(func): assert_allclose( func(r3d, mdelta, cdelta, z_cl, cclcosmo, halo_profile_model="NFW"), defaulttruth, - **TOLERANCE + **TOLERANCE, ) assert_allclose( func(r3d, mdelta, cdelta, z_cl, cclcosmo, massdef="MEAN"), defaulttruth, **TOLERANCE @@ -375,22 +379,25 @@ def test_profiles(modeling_data, profile_init): # Test use_projected_quad if mod.backend == "ccl" and profile_init == "einasto": - if hasattr(mod.hdpm, 'projected_quad'): + if hasattr(mod.hdpm, "projected_quad"): mod.set_projected_quad(True) assert_allclose( mod.eval_surface_density( cfg["SIGMA_PARAMS"]["r_proj"], cfg["SIGMA_PARAMS"]["z_cl"], verbose=True ), cfg["numcosmo_profiles"]["Sigma"], - reltol*1e-1, + reltol * 1e-1, ) assert_allclose( theo.compute_surface_density( - cosmo=cosmo, **cfg["SIGMA_PARAMS"], alpha_ein=alpha_ein, verbose=True, + cosmo=cosmo, + **cfg["SIGMA_PARAMS"], + alpha_ein=alpha_ein, + verbose=True, use_projected_quad=True, ), cfg["numcosmo_profiles"]["Sigma"], - reltol*1e-1, + reltol * 1e-1, ) delattr(mod.hdpm, "projected_quad") @@ -547,11 +554,13 @@ def test_shear_convergence_unittests(modeling_data, profile_init): cfg_inf = load_validation_config() # compute some values - cfg_inf['GAMMA_PARAMS']['z_src'] = 1000. + cfg_inf["GAMMA_PARAMS"]["z_src"] = 1000.0 beta_s_mean = compute_beta_s_mean_from_distribution( - cfg_inf['GAMMA_PARAMS']['z_cluster'], cfg_inf['GAMMA_PARAMS']['z_src'], cosmo) + cfg_inf["GAMMA_PARAMS"]["z_cluster"], cfg_inf["GAMMA_PARAMS"]["z_src"], cosmo + ) beta_s_square_mean = compute_beta_s_square_mean_from_distribution( - cfg_inf['GAMMA_PARAMS']['z_cluster'], cfg_inf['GAMMA_PARAMS']['z_src'], cosmo) + cfg_inf["GAMMA_PARAMS"]["z_cluster"], cfg_inf["GAMMA_PARAMS"]["z_src"], cosmo + ) gammat_inf = theo.compute_tangential_shear(cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"]) kappa_inf = theo.compute_convergence(cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"]) @@ -581,14 +590,14 @@ def test_shear_convergence_unittests(modeling_data, profile_init): theo.compute_reduced_tangential_shear, cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], - approx="notvalid" + approx="notvalid", ) assert_raises( ValueError, theo.compute_magnification, cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], - approx="notvalid" + approx="notvalid", ) assert_raises( ValueError, @@ -596,7 +605,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init): cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], alpha=alpha, - approx="notvalid" + approx="notvalid", ) # test KeyError from invalid key in integ_kwargs assert_raises( @@ -604,14 +613,14 @@ def test_shear_convergence_unittests(modeling_data, profile_init): theo.compute_reduced_tangential_shear, cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], - integ_kwargs={"notavalidkey": 0.0} + integ_kwargs={"notavalidkey": 0.0}, ) assert_raises( KeyError, theo.compute_magnification, cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], - integ_kwargs={"notavalidkey": 0.0} + integ_kwargs={"notavalidkey": 0.0}, ) assert_raises( KeyError, @@ -619,7 +628,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init): cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], alpha=alpha, - integ_kwargs={"notavalidkey": 0.0} + integ_kwargs={"notavalidkey": 0.0}, ) # test ValueError from unsupported z_src_info cfg_inf["GAMMA_PARAMS"]["z_src_info"] = "notvalid" @@ -632,14 +641,14 @@ def test_shear_convergence_unittests(modeling_data, profile_init): theo.compute_reduced_tangential_shear, cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], - approx="order1" + approx="order1", ) assert_raises( ValueError, theo.compute_magnification, cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], - approx="order1" + approx="order1", ) assert_raises( ValueError, @@ -647,7 +656,7 @@ def test_shear_convergence_unittests(modeling_data, profile_init): cosmo=cosmo, **cfg_inf["GAMMA_PARAMS"], alpha=2, - approx="order1" + approx="order1", ) # test z_src_info = 'beta' @@ -1065,7 +1074,7 @@ def test_compute_magnification_bias(modeling_data): assert_allclose( theo.compute_magnification_bias_from_magnification(magnification[0], alpha[0]), truth[0][0], - **TOLERANCE + **TOLERANCE, ) assert_allclose( theo.compute_magnification_bias_from_magnification(magnification, alpha), truth, **TOLERANCE @@ -1075,7 +1084,7 @@ def test_compute_magnification_bias(modeling_data): np.array(magnification), np.array(alpha) ), np.array(truth), - **TOLERANCE + **TOLERANCE, ) diff --git a/tests/test_utils.py b/tests/test_utils.py index 0b2db3591..a70524472 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -11,6 +11,7 @@ convert_shapes_to_epsilon, arguments_consistency, validate_argument, + DiffArray, ) from clmm.redshift import distributions as zdist @@ -535,6 +536,21 @@ def test_validate_argument(): ) +def test_diff_array(): + """test validate argument""" + # Validate diffs + assert DiffArray([1, 2]) == DiffArray([1, 2]) + assert DiffArray([1, 2]) == DiffArray(np.array([1, 2])) + assert DiffArray([1, 2]) != DiffArray(np.array([2, 2])) + assert DiffArray([1, 2]) != DiffArray(np.array([1, 2, 3])) + assert DiffArray([1, 2]) != None + # Validate prints + arr = DiffArray([1, 2]) + assert str(arr.value) == arr.__repr__() + arr = DiffArray(range(10)) + assert str(arr.value) != arr.__repr__() + + def test_beta_functions(modeling_data): z_cl = 1.0 z_src = [2.4, 2.1]