Clean up old functions (#8)

* update nbs

* rm unused functions

* simplify loading

* optionally return counts as well

* update TODO

* update TODO

README.md

@@ -1,17 +1,15 @@
 # CSiBORG tools
 
 ## :scroll: Short-term TODO
-- [ ] Make a nice plot comparing the SZ clusters and their number density
-- [x] Compare empirical $M_{500c}$ to the NFW expectation.
+- [ ] Add half-mass radius and its corresponding concentration.
+- [ ] Model the Planck SZ selection function.
 - [ ] Calculate catalogues for all realisations.
-- [x] Add shortcut function for loading a catalogue
 
 
 ## :hourglass: Long-term TODO
-- [ ] Calculate the halo spin
-- [ ] Model the Planck SZ selection function.
+- [ ] Calculate the halo spin.
 - [ ] Calculate the cross-correlation in CSiBORG. Should see the scale of the constraints?
-- [ ] Improve file naming system
+- [ ] Calculate DM environmental properties.
 
 
 ## :bulb: Open questions

csiborgtools/match/num_density.py

@@ -51,7 +51,7 @@ def binned_counts(x, bins):
     return centres, out
 
 
-def number_density(data, feat, bins, max_dist, to_log10):
+def number_density(data, feat, bins, max_dist, to_log10, return_counts=False):
     """
     Calculate volume-limited number density of a feature `feat` from array
     `data`, normalised also by the bin width.
@@ -70,6 +70,9 @@ def number_density(data, feat, bins, max_dist, to_log10):
     to_log10 : bool
         Whether to take a logarithm of base 10 of the feature. If so, then the
         bins must also be logarithmic.
+    return_counts : bool, optional
+        Whether to also return number counts in each bin. By default `False`.
+
 
     Returns
     -------
@@ -80,6 +83,9 @@ def number_density(data, feat, bins, max_dist, to_log10):
         Number density of shape `(n_edges - 1, )`.
     nd_err : 1-dimensional array
         Poissonian uncertainty of `nd` of shape `(n_edges - 1, )`.
+    counts: 1-dimensional array, optional
+        Counts in each bin of shape `(n_edges - 1, )`. Returned only if
+        `return_counts`.
     """
     # Extract the param and optionally convert to log10
     x = data[feat]
@@ -104,4 +110,8 @@ def number_density(data, feat, bins, max_dist, to_log10):
     # Convert bins to linear space if log10
     if to_log10:
         bin_centres = 10**bin_centres
-    return bin_centres, nd, nd_err
+
+    out = (bin_centres, nd, nd_err)
+    if return_counts:
+        out += counts
+    return out

csiborgtools/units/__init__.py

@@ -13,6 +13,5 @@
 # with this program; if not, write to the Free Software Foundation, Inc.,
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
-from .transforms import (cartesian_to_radec, convert_mass_cols,  # noqa
-                         convert_position_cols)  # noqa
+from .transforms import cartesian_to_radec  # noqa
 from .box_units import (BoxUnits, convert_from_boxunits)  # noqa

csiborgtools/units/transforms.py

@@ -16,15 +16,9 @@
 Various coordinate transformations.
 """
 
-
 import numpy
 
 
-little_h = 0.705
-BOXSIZE = 677.7 / little_h  # Mpc. Otherwise positions in [0, 1].
-BOXMASS = 3.749e19  # Msun
-
-
 def cartesian_to_radec(arr, xpar="peak_x", ypar="peak_y", zpar="peak_z"):
     r"""
     Extract `x`, `y`, and `z` coordinates from a record array `arr` and
@@ -61,50 +55,3 @@ def cartesian_to_radec(arr, xpar="peak_x", ypar="peak_y", zpar="peak_z"):
     ra[ra < 0] += 360
 
     return dist, ra, dec
-
-
-def convert_mass_cols(arr, cols):
-    r"""
-    Convert mass columns from box units to :math:`M_{\odot}`. `arr` is passed
-    by reference and is not explicitly returned back.
-
-    Parameters
-    ----------
-    arr : structured array
-        The array whose columns are to be converted.
-    cols : str or list of str
-        The mass columns to be converted.
-
-    Returns
-    -------
-    None
-    """
-    cols = [cols] if isinstance(cols, str) else cols
-    for col in cols:
-        arr[col] *= BOXMASS
-
-
-def convert_position_cols(arr, cols, zero_centered=True):
-    r"""
-    Convert position columns from box units to :math:`\mathrm{Mpc}`. `arr` is
-    passed by reference and is not explicitly returned back.
-
-    Parameters
-    ----------
-    arr : structured array
-        The array whose columns are to be converted.
-    cols : str or list of str
-        The mass columns to be converted.
-    zero_centered : bool, optional
-        Whether to translate the well-resolved origin in the centre of the
-        simulation to the :math:`(0, 0 , 0)` point. By default `True`.
-
-    Returns
-    -------
-    None
-    """
-    cols = [cols] if isinstance(cols, str) else cols
-    for col in cols:
-        arr[col] *= BOXSIZE
-        if zero_centered:
-            arr[col] -= BOXSIZE / 2

scripts/utils.py

@@ -19,7 +19,6 @@ Notebook utility functions.
 
 import numpy
 from os.path import join
-from tqdm import trange
 from astropy.cosmology import FlatLambdaCDM
 
 try:
@@ -33,40 +32,6 @@ Nsplits = 200
 dumpdir = "/mnt/extraspace/rstiskalek/csiborg/"
 
 
-def load_mmain_convert(n):
-    srcdir = "/users/hdesmond/Mmain"
-    arr = csiborgtools.io.read_mmain(n, srcdir)
-
-    csiborgtools.utils.convert_mass_cols(arr, "mass_cl")
-    csiborgtools.utils.convert_position_cols(
-        arr, ["peak_x", "peak_y", "peak_z"])
-    csiborgtools.utils.flip_cols(arr, "peak_x", "peak_z")
-
-    d, ra, dec = csiborgtools.utils.cartesian_to_radec(arr)
-    arr = csiborgtools.utils.add_columns(
-        arr, [d, ra, dec], ["dist", "ra", "dec"])
-    return arr
-
-
-def load_mmains(N=None, verbose=True):
-    ids = csiborgtools.io.get_csiborg_ids("/mnt/extraspace/hdesmond")
-    N = ids.size if N is None else N
-    if N > ids.size:
-        raise ValueError("`N` cannot be larger than 101.")
-    # If N less than num of CSiBORG, then radomly choose
-    if N == ids.size:
-        choices = numpy.arange(N)
-    else:
-        choices = numpy.random.choice(ids.size, N, replace=False)
-
-    out = [None] * N
-    iters = trange(N) if verbose else range(N)
-    for i in iters:
-        j = choices[i]
-        out[i] = load_mmain_convert(ids[j])
-    return out
-
-
 def load_processed(Nsim, Nsnap):
     simpath = csiborgtools.io.get_sim_path(Nsim)
     outfname = join(
@@ -76,6 +41,7 @@ def load_processed(Nsim, Nsnap):
     # Add mmain
     mmain = csiborgtools.io.read_mmain(Nsim, "/mnt/zfsusers/hdesmond/Mmain")
     data = csiborgtools.io.merge_mmain_to_clumps(data, mmain)
+    csiborgtools.utils.flip_cols(data, "peak_x", "peak_z")
     # Cut on numbre of particles and finite m200
     data = data[(data["npart"] > 100) & numpy.isfinite(data["m200"])]
 
@@ -85,6 +51,10 @@ def load_processed(Nsim, Nsnap):
                     "r200", "r500", "Rs", "rho0", "peak_x", "peak_y", "peak_z"]
     data = csiborgtools.units.convert_from_boxunits(
         data, convert_cols, boxunits)
+    # Now calculate spherical coordinates
+    d, ra, dec = csiborgtools.units.cartesian_to_radec(data)
+    data = csiborgtools.utils.add_columns(
+        data, [d, ra, dec], ["dist", "ra", "dec"])
     return data