Merge branch 'main' into issue_#635_coordinate_implementation

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 <a name="access_to_the_proper_environment_on_cori"></a>
 
 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:

clmm/__init__.py

@@ -26,4 +26,4 @@
 )
 from . import support
 
-__version__ = "1.12.5"
+__version__ = "1.13.2"

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"):

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(

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

clmm/utils/__init__.py

@@ -41,6 +41,7 @@
     _validate_dec,
     _validate_is_deltasigma_sigma_c,
     _validate_coordinate_system,
+    DiffArray,
 )
 
 from .units import (

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]

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

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,

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,

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",

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,

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)