Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 30 additions & 0 deletions crow/cluster_modules/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
"""Cluster modules package.

Keep this module lightweight: do not import the submodules or classes eagerly
here because several submodules import the top-level ``crow`` package and that
can produce circular imports when the top-level ``crow.__init__`` imports
things from ``cluster_modules``.

To import classes or modules use the explicit submodule path, for example::

from crow.cluster_modules.abundance import ClusterAbundance
from crow.cluster_modules.shear_profile import ClusterShearProfile

Sphinx/apidoc will still find and document the submodules even if they are not
imported here; this file only needs to exist so Python treats the directory as
a package.
"""

# Expose the expected subpackage/module names for convenience. Do NOT import
# the modules at package import time to avoid circular import problems.
__all__ = [
"abundance",
"completeness_models",
"kernel",
"parameters",
"purity_models",
"_clmm_patches",
"shear_profile",
"shear_profile_parallel",
"mass_proxy",
]
58 changes: 54 additions & 4 deletions crow/cluster_modules/abundance.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,18 @@ class ClusterAbundance:
an area on the sky, a halo mass function, as well as multiple kernels, where
each kernel represents a different distribution involved in the final cluster
abundance integrand.

Attributes
----------
cosmo : pyccl.cosmology.Cosmology or None
Cosmology object used for predictions. Set via the `cosmo` property.
halo_mass_function : callable
Halo mass function object compatible with PyCCL's MassFunc interface.
parameters : Parameters
Container for optional model parameters.
_hmf_cache : dict
Cache for previously computed halo mass function evaluations keyed by
(log_mass, scale_factor).
"""

@property
Expand All @@ -47,9 +59,28 @@ def __init__(
def comoving_volume(
self, z: npt.NDArray[np.float64], sky_area: float = 0
) -> npt.NDArray[np.float64]:
"""The differential comoving volume given area sky_area at redshift z.

:param sky_area: The area of the survey on the sky in square degrees.
"""Differential comoving volume for a given sky area.

Parameters
----------
z : array_like
Redshift or array of redshifts at which to compute the differential
comoving volume (dV/dz per steradian).
sky_area : float, optional
Survey area in square degrees. Default 0 returns per-steradian volume.

Returns
-------
numpy.ndarray
Differential comoving volume (same shape as `z`) multiplied by the
survey area (converted to steradians). Units: [Mpc/h]^3 (consistent
with the internal PyCCL conventions used here).

Notes
-----
This uses PyCCL background helpers to evaluate the angular diameter
distance and h(z) factor. If `sky_area` is zero the returned array is
the per-steradian dV/dz.
"""
assert self.cosmo is not None
scale_factor = 1.0 / (1.0 + z)
Expand All @@ -73,7 +104,26 @@ def mass_function(
log_mass: npt.NDArray[np.float64],
z: npt.NDArray[np.float64],
) -> npt.NDArray[np.float64]:
"""The mass function at z and mass."""
"""Evaluate the halo mass function at given log-mass and redshift.

Parameters
----------
log_mass : array_like
Array of log10 halo masses (M_sun).
z : array_like
Array of redshifts matching `log_mass`.

Returns
-------
numpy.ndarray
Array with mass function values (dn/dlnM or as provided by the
configured `halo_mass_function`) evaluated at each (mass, z).

Notes
-----
Results are cached in `_hmf_cache` keyed by (log_mass, scale_factor)
to avoid repeated expensive evaluations for identical inputs.
"""
scale_factor = 1.0 / (1.0 + z)
return_vals = []

Expand Down
103 changes: 92 additions & 11 deletions crow/cluster_modules/completeness_models.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""The cluster completeness module.

This module holds the classes that define the kernels that can be included
in the cluster abundance integrand.
This module holds the classes that define completeness kernels that can be included
in the cluster prediction integrand.
"""

import numpy as np
Expand All @@ -11,10 +11,15 @@


class Completeness:
"""The completeness kernel for the numcosmo simulated survey.
"""The completeness kernel base class.

This kernel will affect the integrand by accounting for the incompleteness
of a cluster selection.
This kernel affects the prediction integrand by accounting for the incompleteness
of a cluster selection. Subclasses should implement the ``distribution`` method.

Attributes
----------
parameters : Parameters, optional
Container for completeness model parameters (defined by subclasses).
"""

def __init__(self):
Expand All @@ -25,7 +30,22 @@ def distribution(
log_mass: npt.NDArray[np.float64],
z: npt.NDArray[np.float64],
) -> npt.NDArray[np.float64]:
"""Evaluates and returns the completeness contribution to the integrand."""
"""Evaluate the completeness kernel contribution.

Parameters
----------
log_mass : array_like
Array of log10 halo masses (units: Msun).
z : array_like
Array of redshifts matching ``log_mass``.

Returns
-------
numpy.ndarray
Array of completeness values in the range [0, 1] with the same
broadcastable shape as the inputs. Subclasses should guarantee the
output dtype is floating point.
"""
raise NotImplementedError


Expand All @@ -38,10 +58,25 @@ def distribution(


class CompletenessAguena16(Completeness):
"""The completeness kernel for the numcosmo simulated survey.

This kernel will affect the integrand by accounting for the incompleteness
of a cluster selection.
"""Completeness model following Aguena et al. (2016) parametrisation.

The model uses a pivot mass and a redshift-dependent power-law index to
compute a sigmoid-like completeness as a function of mass and redshift.

Parameters
----------
(set during initialization)
a_n, b_n : float
Parameters controlling the redshift evolution of the power-law index.
a_logm_piv, b_logm_piv : float
Parameters controlling the pivot mass (in log10 units) and its
redshift evolution.

Attributes
----------
parameters : Parameters
Container holding the parameter values; defaults are defined in
``REDMAPPER_DEFAULT_PARAMETERS``.
"""

def __init__(
Expand All @@ -50,13 +85,37 @@ def __init__(
self.parameters = Parameters({**REDMAPPER_DEFAULT_PARAMETERS})

def _mpiv(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
"""Return the pivot mass (not in log) at redshift `z`.

Parameters
----------
z : array_like
Redshift or array of redshifts.

Returns
-------
numpy.ndarray
Pivot mass values in Msun (10**log_mpiv). Returned dtype is float64.
"""
log_mpiv = self.parameters["a_logm_piv"] + self.parameters["b_logm_piv"] * (
1.0 + z
)
mpiv = 10.0**log_mpiv
return mpiv.astype(np.float64)

def _nc(self, z: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
"""Return the redshift-dependent power-law index nc(z).

Parameters
----------
z : array_like
Redshift or array of redshifts.

Returns
-------
numpy.ndarray
The value of the power-law index at each provided redshift.
"""
nc = self.parameters["a_n"] + self.parameters["b_n"] * (1.0 + z)
assert isinstance(nc, np.ndarray)
return nc
Expand All @@ -66,7 +125,29 @@ def distribution(
log_mass: npt.NDArray[np.float64],
z: npt.NDArray[np.float64],
) -> npt.NDArray[np.float64]:
"""Evaluates and returns the completeness contribution to the integrand."""
r"""Compute the completeness fraction for given mass and redshift.
The completeness is given by

.. math::
c(M, z) = \frac{\left(M / M_{\rm piv}(z)\right)^{n_c(z)}}
{1 + \left(M / M_{\rm piv}(z)\right)^{n_c(z)}}

where M = 10^{\text{log\_mass}}, M_{\rm piv}(z) is returned by _mpiv(z),
and n_c(z) is returned by _nc(z).

Parameters
----------
log_mass : array_like
Array of log10 halo masses (Msun).
z : array_like
Array of redshifts matching ``log_mass``.

Returns
-------
numpy.ndarray
Completeness values in the interval [0, 1] with shape matching the
broadcasted inputs. dtype is float64.
"""

mass_norm_pow = (10.0**log_mass / self._mpiv(z)) ** self._nc(z)

Expand Down
34 changes: 25 additions & 9 deletions crow/cluster_modules/kernel.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,29 +9,45 @@


class TrueMass:
"""The true mass kernel.
"""True-mass kernel used in abundance integrals.

Assuming we measure the true mass, this will always be 1.
This kernel represents the case where the observed mass equals the true
halo mass. It therefore contributes a multiplicative factor of unity to
the abundance integrand and does not alter the mass distribution.

Notes
-----
The distribution method returns a NumPy array (dtype float) that can be
broadcast with other integrand factors.
"""

def distribution(self) -> npt.NDArray[np.float64]:
"""Evaluates and returns the mass distribution contribution to the integrand.
"""Evaluate and return the mass kernel contribution.

We have set this to 1.0 (i.e. it does not affect the mass distribution)
Returns
-------
numpy.ndarray
Array containing the value 1.0. This can be broadcast to the shape
required by the integrand and is provided as float64.
"""
return np.atleast_1d(1.0)


class SpectroscopicRedshift:
"""The spec-z kernel.
"""Spectroscopic-redshift kernel for abundance integrals.

Assuming the spectroscopic redshift has no uncertainties, this is akin to
multiplying by 1.
Represents the idealized case where cluster redshifts are known exactly
(spectroscopic precision). The kernel thus contributes a factor of unity
to the redshift part of the integrand.
"""

def distribution(self) -> npt.NDArray[np.float64]:
"""Evaluates and returns the z distribution contribution to the integrand.
"""Evaluate and return the redshift kernel contribution.

We have set this to 1.0 (i.e. it does not affect the redshift distribution)
Returns
-------
numpy.ndarray
Array containing the value 1.0 (dtype float64). This is intended
to be broadcast with other integrand arrays.
"""
return np.atleast_1d(1.0)
Loading
Loading