Issue/628/rad min max in ensemble (#631)

* add radius min/max in ensamble and validation for its consistency

* pass rmin/max to ClusterEnsemble.stacked_data

* Update version to 1.13.0

---------

Co-authored-by: Hsin Fan <[email protected]>

Author: m-aguena

clmm/__init__.py

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

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

@@ -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

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

@@ -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]

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,

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)