diff --git a/INSTALL.md b/INSTALL.md index 68d1fbd7e..05ca78243 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -64,6 +64,17 @@ Now, you can install CLMM and its dependencies as python setup.py install # build from source ``` +### Local environment for CLMM + +Alternatively, you can create a new local environment by running + +```bash + conda env create -f environment.yml + conda activate clmm +``` + +You can now install CLMM in a local and stable environment with the usual procedure. + ## Access to the proper environment on cori.nersc.gov If you have access to NERSC, this will likely be the easiest to make sure you have the appropriate environment. After logging into cori.nersc.gov, you will need to execute the following. We recommend executing line-by-line to avoid errors: diff --git a/clmm/__init__.py b/clmm/__init__.py index 0111be973..a4b97c140 100644 --- a/clmm/__init__.py +++ b/clmm/__init__.py @@ -26,4 +26,4 @@ ) from . import support -__version__ = "1.12.5" +__version__ = "1.13.2" 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 e5e5e18b8..8f764a70a 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 c5f11c4a7..96ae740dd 100644 --- a/clmm/gcdata.py +++ b/clmm/gcdata.py @@ -167,7 +167,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 7aacde7b9..f85895b1b 100644 --- a/clmm/utils/validation.py +++ b/clmm/utils/validation.py @@ -240,3 +240,37 @@ def _validate_coordinate_system(loc, argname): 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)""" + + 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/environment.yml b/environment.yml new file mode 100644 index 000000000..86a876f6b --- /dev/null +++ b/environment.yml @@ -0,0 +1,20 @@ +# CLMM Conda environment +name: clmm +channels: + - conda-forge + - defaults +dependencies: + - python==3.10.2 + - pip>=21.0 + - jupyter>=1.0 + - numpy==1.25.2 + - scipy==1.11.2 + - astropy==5.3.2 + - healpy == 1.16.5 + - matplotlib==3.7.2 + - gsl==2.7 + - cmake==3.30.0 + - pyccl==3.0 + - numcosmo==0.18.2 + - pytest=7.4.0 + - sphinx==7.2.4 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/pyproject.toml b/pyproject.toml index 5bd0e0a77..efc4a2346 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -3,7 +3,7 @@ line-length = 100 [tool.pylint] max-line-length = 100 -disable = ["R0913"] +disable = ["R0913", "R0917"] good-names = ["ra", "z", "z1", "z2", "a", "a1", "a2", "gt", "gx", "mu", "da", "i", "H0", "Omega_b0", "Omega_dm0", "Omega_k0", 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 851bcc655..8d52ea16a 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]