diff --git a/docs/about/lf_overview.rst b/docs/about/lf_overview.rst index 97d3d84..53b50eb 100644 --- a/docs/about/lf_overview.rst +++ b/docs/about/lf_overview.rst @@ -5,7 +5,7 @@ |lfkitlogo| Luminosity Functions ================================ -LFKit provides a public luminosity-function interface through +LFKit provides a public luminosity function interface through :class:`~lfkit.api.lumfunc.LuminosityFunction`. For most users, the recommended import is: @@ -14,14 +14,14 @@ For most users, the recommended import is: from lfkit import LuminosityFunction -The ``LuminosityFunction`` object defines and evaluates luminosity-function +The ``LuminosityFunction`` object defines and evaluates luminosity function models in rest-frame absolute-magnitude space. It can also evaluate luminosity functions from apparent magnitudes, convert between apparent and absolute magnitudes, and compute number-density quantities for magnitude-limited catalog selections. File reading is intentionally not handled by this API. Catalog-derived -luminosity-function parameters, magnitude limits, or correction models should +luminosity function parameters, magnitude limits, or correction models should be loaded elsewhere and passed in as scalars, arrays, or correction objects. @@ -139,7 +139,7 @@ the standard Schechter form. Built-in parameter evolution models ----------------------------------- -Redshift-dependent luminosity-function parameters are handled by parameter +Redshift-dependent luminosity function parameters are handled by parameter evolution models. The built-in options are: @@ -176,7 +176,7 @@ To inspect the available models: models = LuminosityFunction.available_parameter_models() -Custom parameter-evolution models can be registered through: +Custom parameter evolution models can be registered through: .. code-block:: python @@ -191,7 +191,7 @@ return NumPy-compatible parameter values. Evaluating from apparent magnitudes ----------------------------------- -The luminosity-function object can evaluate from apparent magnitudes. +The luminosity function object can evaluate from apparent magnitudes. In this case, LFKit first converts apparent magnitude :math:`m` into absolute magnitude :math:`M` using @@ -281,6 +281,54 @@ This computes the number density inside the requested absolute-magnitude interval. +LF-weighted redshift density +---------------------------- + +Use :meth:`~lfkit.api.lumfunc.LuminosityFunction.lf_redshift_density` to build +an LF-weighted redshift trend for a magnitude-limited sample. + +This computes the magnitude-integrated luminosity function as a function of +redshift. Optionally, it can also include a redshift-dependent volume weight +and normalize the result over the supplied redshift grid. + +.. code-block:: python + + density = lf.lf_redshift_density( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + corrections=corr, + normalize=True, + ) + +This is useful for constructing LF-dependent redshift trends. +Without a volume weight, the result is the luminosity-function +selection factor only, not a full survey redshift distribution. + +A volume-weighted trend can be computed by passing a callable: + +.. code-block:: python + + import pyccl as ccl + + def volume_weight_fn(z): + a = 1.0 / (1.0 + z) + chi = ccl.comoving_radial_distance(cosmo, a) + h_over_h0 = ccl.h_over_h0(cosmo, a) + return chi**2 / h_over_h0 + + density = lf.lf_redshift_density( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + corrections=corr, + volume_weight_fn=volume_weight_fn, + normalize=True, + ) + + Magnitude-limited catalog completeness -------------------------------------- @@ -398,7 +446,7 @@ There are two different places where luminosity evolution can enter an analysis: 1. the apparent-to-absolute magnitude conversion through :math:`E(z)`, -2. the luminosity-function model through redshift evolution of :math:`M_\star(z)`. +2. the luminosity function model through redshift evolution of :math:`M_\star(z)`. For this reason, be careful when using an evolving Schechter model with non-zero ``linear_q`` evolution together with an explicit evolution correction. @@ -410,14 +458,14 @@ separated. Typical workflow ---------------- -A typical luminosity-function workflow is: +A typical luminosity function workflow is: 1. define a cosmology, -2. define a luminosity-function model, +2. define a luminosity function model, 3. evaluate :math:`\phi(M, z)` directly or :math:`\phi(m, z)` from apparent magnitudes, -4. compute integrated, observed, or missing number densities if using a - magnitude-limited catalog. +4. compute integrated number densities, LF-weighted redshift densities, or + observed/missing number densities for magnitude-limited selections. For example: @@ -455,6 +503,15 @@ For example: corrections=corr, ) + density = lf.lf_redshift_density( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + corrections=corr, + normalize=True, + ) + Lower-level functions --------------------- @@ -465,8 +522,8 @@ downstream package interfaces. The lower-level functions in ``lfkit.photometry`` are useful when: - testing individual mathematical pieces, -- adding new luminosity-function models, -- adding new parameter-evolution models, +- adding new luminosity function models, +- adding new parameter evolution models, - debugging the magnitude convention, - building new public API objects, - integrating LFKit into specialized workflows. diff --git a/docs/api/lfkit.photometry.lf_integrals.rst b/docs/api/lfkit.photometry.lf_integrals.rst new file mode 100644 index 0000000..55f6f13 --- /dev/null +++ b/docs/api/lfkit.photometry.lf_integrals.rst @@ -0,0 +1,7 @@ +lfkit.photometry.lf\_integrals module +===================================== + +.. automodule:: lfkit.photometry.lf_integrals + :members: + :show-inheritance: + :undoc-members: diff --git a/docs/api/lfkit.photometry.lf_redshift_density.rst b/docs/api/lfkit.photometry.lf_redshift_density.rst new file mode 100644 index 0000000..576d8ab --- /dev/null +++ b/docs/api/lfkit.photometry.lf_redshift_density.rst @@ -0,0 +1,7 @@ +lfkit.photometry.lf\_redshift\_density module +============================================= + +.. automodule:: lfkit.photometry.lf_redshift_density + :members: + :show-inheritance: + :undoc-members: diff --git a/docs/api/lfkit.photometry.rst b/docs/api/lfkit.photometry.rst index 890f6a8..254d341 100644 --- a/docs/api/lfkit.photometry.rst +++ b/docs/api/lfkit.photometry.rst @@ -8,7 +8,9 @@ Submodules :maxdepth: 2 lfkit.photometry.catalog_completeness + lfkit.photometry.lf_integrals lfkit.photometry.lf_parameter_models + lfkit.photometry.lf_redshift_density lfkit.photometry.luminosities lfkit.photometry.luminosity_function lfkit.photometry.magnitudes diff --git a/docs/api/lfkit.utils.evaluators.rst b/docs/api/lfkit.utils.evaluators.rst new file mode 100644 index 0000000..8e8ffa8 --- /dev/null +++ b/docs/api/lfkit.utils.evaluators.rst @@ -0,0 +1,7 @@ +lfkit.utils.evaluators module +============================= + +.. automodule:: lfkit.utils.evaluators + :members: + :show-inheritance: + :undoc-members: diff --git a/docs/api/lfkit.utils.rst b/docs/api/lfkit.utils.rst index 80b4e19..6e06a20 100644 --- a/docs/api/lfkit.utils.rst +++ b/docs/api/lfkit.utils.rst @@ -8,6 +8,7 @@ Submodules :maxdepth: 2 lfkit.utils.download_poggianti97_data + lfkit.utils.evaluators lfkit.utils.interpolation lfkit.utils.io lfkit.utils.types diff --git a/src/lfkit/api/lumfunc.py b/src/lfkit/api/lumfunc.py index c4e6e11..43712ad 100644 --- a/src/lfkit/api/lumfunc.py +++ b/src/lfkit/api/lumfunc.py @@ -16,7 +16,7 @@ from __future__ import annotations -from collections.abc import Mapping +from collections.abc import Callable, Mapping from typing import TYPE_CHECKING import numpy as np @@ -31,7 +31,17 @@ ) from lfkit.photometry.magnitudes import ( absolute_magnitude, + absolute_magnitude_from_luminosity_distance, apparent_magnitude, + apparent_magnitude_from_luminosity_distance, +) +from lfkit.photometry.lf_integrals import ( + cumulative_number_density, + integrated_luminosity_density, + integrated_number_density, + lf_weighted_integral, + mean_luminosity, + selection_weighted_number_density, ) from lfkit.photometry.lf_parameter_models import ( available_lf_parameter_models, @@ -40,13 +50,16 @@ register_m_star_model, register_phi_star_model, ) +from lfkit.photometry.lf_redshift_density import ( + lf_integrated_number_density, + lf_weighted_redshift_density, +) from lfkit.photometry.catalog_completeness import ( absolute_magnitude_limit, catalog_completeness_fraction, out_of_catalog_fraction, observed_number_density, missing_number_density, - integrated_number_density, ) from lfkit.utils.types import ( Cosmology, @@ -152,6 +165,153 @@ def apparent_magnitude( e_correction=e_corr, ) + def absolute_magnitude_from_luminosity_distance( + self, + apparent_mag: FloatInput, + luminosity_distance_mpc: FloatInput, + *, + z: FloatInput, + corrections: Corrections | None = None, + ) -> FloatArray: + """Convert apparent magnitude to absolute magnitude from luminosity distance. + + Args: + apparent_mag: Apparent magnitude values. + luminosity_distance_mpc: Luminosity distance values in Mpc. + corrections: Optional object providing k-correction and e-correction values. + z: Redshift values where corrections are evaluated. + + Returns: + Absolute magnitudes using ``M = m - mu - K + E``. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return absolute_magnitude_from_luminosity_distance( + apparent_mag, + luminosity_distance_mpc, + k_correction=k_corr, + e_correction=e_corr, + ) + + def apparent_magnitude_from_luminosity_distance( + self, + absolute_mag: FloatInput, + luminosity_distance_mpc: FloatInput, + *, + z: FloatInput, + corrections: Corrections | None = None, + ) -> FloatArray: + """Convert absolute magnitude to apparent magnitude from luminosity distance. + + Args: + absolute_mag: Absolute magnitude values. + luminosity_distance_mpc: Luminosity distance values in Mpc. + z: Redshift values where corrections are evaluated. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Apparent magnitudes using ``m = M + mu + K - E``. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return apparent_magnitude_from_luminosity_distance( + absolute_mag, + luminosity_distance_mpc, + k_correction=k_corr, + e_correction=e_corr, + ) + + def lf_integrated_number_density( + self, + z: FloatInput, + *, + m_lim: float, + m_bright: float, + n_m: int = 512, + luminosity_distance_mpc_fn: Callable[[FloatArray], FloatArray], + corrections: Corrections | None = None, + ) -> FloatArray: + """Return LF-integrated number density for an apparent-magnitude limit. + + This integrates the luminosity function from ``m_bright`` to the + absolute-magnitude limit implied by ``m_lim`` and the supplied + luminosity-distance callable. + + Args: + z: Redshift values. + m_lim: Apparent-magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration bound. + n_m: Number of magnitude-grid points used in the integration. + luminosity_distance_mpc_fn: Callable returning luminosity distance + in Mpc as a function of redshift. + corrections: Optional object providing k-correction and + e-correction values. + + Returns: + LF-integrated number density evaluated at redshift. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return lf_integrated_number_density( + z, + self._as_callable(), + m_lim=m_lim, + m_bright=m_bright, + n_m=n_m, + luminosity_distance_mpc_fn=luminosity_distance_mpc_fn, + k_correction=k_corr, + evolution_correction=e_corr, + ) + + def lf_weighted_redshift_density( + self, + z: FloatInput, + *, + m_lim: float, + m_bright: float, + n_m: int = 512, + luminosity_distance_mpc_fn: Callable[[FloatArray], FloatArray], + volume_weight_fn: Callable[[FloatArray], FloatArray], + corrections: Corrections | None = None, + normalize: bool = True, + ) -> FloatArray: + """Return an LF-weighted redshift-density curve. + + This computes an LF-selected redshift distribution by integrating the + luminosity function up to the apparent-magnitude limit and multiplying + by a supplied redshift or volume weight. + + Args: + z: Redshift values. + m_lim: Apparent-magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration bound. + n_m: Number of magnitude-grid points used in the integration. + luminosity_distance_mpc_fn: Callable returning luminosity distance + in Mpc as a function of redshift. + volume_weight_fn: Callable returning the redshift or volume weight. + corrections: Optional object providing k-correction and + e-correction values. + normalize: If True, normalize the returned curve to integrate to + one over redshift. + + Returns: + LF-weighted redshift-density curve. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return lf_weighted_redshift_density( + z, + self._as_callable(), + m_lim=m_lim, + m_bright=m_bright, + n_m=n_m, + luminosity_distance_mpc_fn=luminosity_distance_mpc_fn, + volume_weight_fn=volume_weight_fn, + k_correction=k_corr, + evolution_correction=e_corr, + normalize=normalize, + ) + @classmethod def schechter( cls, @@ -463,6 +623,103 @@ def integrated_number_density( n_m=n_m, ) + def lf_weighted_integral( + self, + z: FloatInput, + *, + m_bright: ParameterValue, + m_faint: ParameterValue, + weight_fn: Callable[[FloatArray, FloatArray], FloatArray], + n_m: int = 512, + ) -> FloatArray: + """Integrate the LF with a user-supplied magnitude-redshift weight.""" + return lf_weighted_integral( + z, + self._as_callable(), + m_bright=m_bright, + m_faint=m_faint, + weight_fn=weight_fn, + n_m=n_m, + ) + + def selection_weighted_number_density( + self, + z: FloatInput, + *, + m_bright: ParameterValue, + m_faint: ParameterValue, + selection_fn: Callable[[FloatArray, FloatArray], FloatArray], + n_m: int = 512, + ) -> FloatArray: + """Return LF number density weighted by a selection function.""" + return selection_weighted_number_density( + z, + self._as_callable(), + m_bright=m_bright, + m_faint=m_faint, + selection_fn=selection_fn, + n_m=n_m, + ) + + def integrated_luminosity_density( + self, + z: FloatInput, + *, + m_bright: ParameterValue, + m_faint: ParameterValue, + m_reference: float = 0.0, + n_m: int = 512, + ) -> FloatArray: + """Return luminosity density over an absolute-magnitude range.""" + return integrated_luminosity_density( + z, + self._as_callable(), + m_bright=m_bright, + m_faint=m_faint, + m_reference=m_reference, + n_m=n_m, + ) + + def mean_luminosity( + self, + z: FloatInput, + *, + m_bright: ParameterValue, + m_faint: ParameterValue, + m_reference: float = 0.0, + n_m: int = 512, + ) -> FloatArray: + """Return mean luminosity over an absolute-magnitude range.""" + return mean_luminosity( + z, + self._as_callable(), + m_bright=m_bright, + m_faint=m_faint, + m_reference=m_reference, + n_m=n_m, + ) + + def cumulative_number_density( + self, + z: FloatInput, + *, + m_threshold: ParameterValue, + m_bright: ParameterValue, + m_faint: ParameterValue, + brighter_than: bool = True, + n_m: int = 512, + ) -> FloatArray: + """Return cumulative LF number density around a magnitude threshold.""" + return cumulative_number_density( + z, + self._as_callable(), + m_threshold=m_threshold, + m_bright=m_bright, + m_faint=m_faint, + brighter_than=brighter_than, + n_m=n_m, + ) + def observed_number_density( self, cosmo_obj: Cosmology, diff --git a/src/lfkit/photometry/catalog_completeness.py b/src/lfkit/photometry/catalog_completeness.py index 540fccf..5dc0ff1 100644 --- a/src/lfkit/photometry/catalog_completeness.py +++ b/src/lfkit/photometry/catalog_completeness.py @@ -6,8 +6,8 @@ priors for gravitational-wave cosmology. The utilities convert an apparent magnitude limit into an absolute-magnitude -limit, integrate a luminosity function over finite absolute-magnitude ranges, -and return number densities or fractions. +limit and call the generic LF integration helpers to return number densities +or fractions. The core API accepts a luminosity-function callable with signature @@ -22,6 +22,9 @@ import numpy as np +from lfkit.photometry.lf_integrals import ( + integrated_number_density as _integrated_number_density, +) from lfkit.photometry.magnitudes import absolute_magnitude from lfkit.utils.types import ( Cosmology, @@ -34,7 +37,6 @@ __all__ = [ "absolute_magnitude_limit", - "integrated_number_density", "observed_number_density", "missing_number_density", "catalog_completeness_fraction", @@ -65,7 +67,8 @@ def absolute_magnitude_limit( z: Redshift value or array-like of redshift values. m_lim: Apparent magnitude limit of the catalog. h: Optional dimensionless Hubble parameter used in the - distance-modulus convention. + distance-modulus convention. If not provided, this is read from + ``cosmo_obj["h"]``. k_correction: Optional k-correction term(s). e_correction: Optional evolution-correction term(s). @@ -80,54 +83,18 @@ def absolute_magnitude_limit( if not np.isfinite(m_lim): raise ValueError("m_lim must be finite.") + h_resolved = _resolve_h(cosmo_obj, h) + return absolute_magnitude( cosmo_obj, z_arr, m_lim, - h=h, + h=h_resolved, k_correction=k_correction, e_correction=e_correction, ) -def integrated_number_density( - z: FloatInput, - lf: LuminosityFunction, - *, - m_bright: FloatInput, - m_faint: FloatInput, - n_m: int = 512, -) -> FloatArray: - r"""Return finite-range number density from a luminosity function. - - This computes - - .. math:: - - n(z) = \int_{M_{\mathrm{bright}}(z)}^{M_{\mathrm{faint}}(z)} - \phi(M, z) \, dM. - - Magnitudes are ordered so that more negative values are brighter. - - Args: - z: Redshift value or array-like of redshift values. - lf: Luminosity-function callable with signature ``lf(M, z)``. - m_bright: Bright absolute-magnitude bound. May be scalar or array-like. - m_faint: Faint absolute-magnitude bound. May be scalar or array-like. - n_m: Number of magnitude-grid points used for the integral. - - Returns: - NumPy array of number densities evaluated at ``z``. - """ - return _integrated_number_density_between_bounds( - z, - lf, - m_lower=m_bright, - m_upper=m_faint, - n_m=n_m, - ) - - def observed_number_density( cosmo_obj: Cosmology, z: FloatInput, @@ -182,11 +149,11 @@ def observed_number_density( observed_upper = np.minimum(m_abs_lim, m_faint) - return _integrated_number_density_between_bounds( + return _integrated_number_density( z_arr, lf, - m_lower=m_bright, - m_upper=observed_upper, + m_bright=m_bright, + m_faint=observed_upper, n_m=n_m, ) @@ -245,11 +212,11 @@ def missing_number_density( missing_lower = np.maximum(m_abs_lim, m_bright) - return _integrated_number_density_between_bounds( + return _integrated_number_density( z_arr, lf, - m_lower=missing_lower, - m_upper=m_faint, + m_bright=missing_lower, + m_faint=m_faint, n_m=n_m, ) @@ -375,111 +342,27 @@ def out_of_catalog_fraction( return np.asarray(1.0 - completeness, dtype=float) -def _integrated_number_density_between_bounds( - z: FloatInput, - lf: LuminosityFunction, - *, - m_lower: FloatInput, - m_upper: FloatInput, - n_m: int, -) -> FloatArray: - r"""Return finite-range number density between magnitude bounds.""" - z_arr = validate_array(z, name="z") - m_lower_arr = validate_array(m_lower, name="m_lower") - m_upper_arr = validate_array(m_upper, name="m_upper") - - _validate_integration_inputs( - z=z_arr, - m_lower=m_lower_arr, - m_upper=m_upper_arr, - n_m=n_m, - ) - - z_b, m_lower_b, m_upper_b = np.broadcast_arrays( - z_arr, - m_lower_arr, - m_upper_arr, - ) - - density = np.zeros_like(z_b, dtype=float) - valid = m_upper_b > m_lower_b - - if not np.any(valid): - return density - - m_grid = _magnitude_grid( - m_lower=m_lower_b[valid], - m_upper=m_upper_b[valid], - n_m=n_m, - ) - z_grid = np.broadcast_to(z_b[valid][None, :], m_grid.shape) - phi = _evaluate_lf_on_grid(lf, m_grid=m_grid, z_grid=z_grid) - - density[valid] = np.trapezoid(phi, x=m_grid, axis=0) - - return np.asarray(density, dtype=float) - - -def _magnitude_grid( - *, - m_lower: FloatArray, - m_upper: FloatArray, - n_m: int, -) -> FloatArray: - r"""Return a magnitude grid for column-wise finite-range integration.""" - t = np.linspace(0.0, 1.0, int(n_m), dtype=float) - return np.asarray( - m_lower[None, :] + t[:, None] * (m_upper[None, :] - m_lower[None, :]), - dtype=float, - ) - - -def _evaluate_lf_on_grid( - lf: LuminosityFunction, - *, - m_grid: FloatArray, - z_grid: FloatArray, -) -> FloatArray: - r"""Return LF values evaluated on a magnitude-redshift grid.""" - phi = np.asarray(lf(m_grid, z_grid), dtype=float) - - if phi.shape != m_grid.shape: - try: - phi = np.broadcast_to(phi, m_grid.shape) - except ValueError as exc: - raise ValueError( - "lf(M, z) must return values broadcastable to the shape " - "of the magnitude-redshift integration grid." - ) from exc - - if np.any(~np.isfinite(phi)): - raise ValueError("lf(M, z) returned non-finite values.") - - if np.any(phi < 0): - raise ValueError("lf(M, z) must be non-negative.") - - return np.asarray(phi, dtype=float) - - -def _validate_integration_inputs( - *, - z: FloatArray, - m_lower: FloatArray, - m_upper: FloatArray, - n_m: int, -) -> None: - r"""Validate inputs for finite-range LF integration.""" - if np.any(z < 0): - raise ValueError("Redshift z must be >= 0.") - - if np.any(~np.isfinite(m_lower)): - raise ValueError("m_lower must contain only finite values.") - - if np.any(~np.isfinite(m_upper)): - raise ValueError("m_upper must contain only finite values.") - - if n_m < 2: - raise ValueError("n_m must be at least 2.") +def _resolve_h( + cosmo_obj: Cosmology, + h: float | None, +) -> float: + """Return explicit h or read it from a PyCCL-style cosmology object.""" + if h is not None: + if not np.isfinite(h): + raise ValueError("h must be finite.") + return float(h) + + try: + h_from_cosmo = cosmo_obj["h"] + except (KeyError, TypeError, AttributeError) as exc: + raise ValueError( + "h was not provided and could not be read from cosmo_obj['h']." + ) from exc + + if not np.isfinite(h_from_cosmo): + raise ValueError("cosmo_obj['h'] must be finite.") + + return float(h_from_cosmo) def _fraction( diff --git a/src/lfkit/photometry/lf_integrals.py b/src/lfkit/photometry/lf_integrals.py new file mode 100644 index 0000000..b045cde --- /dev/null +++ b/src/lfkit/photometry/lf_integrals.py @@ -0,0 +1,500 @@ +r"""Luminosity-function integration utilities. + +This module provides generic numerical integrals of luminosity-function +callables over finite absolute-magnitude ranges. + +The core API accepts a luminosity-function callable with signature + + lf(absolute_mag, z) + +where ``absolute_mag`` and ``z`` are NumPy arrays that can be broadcast +together. This keeps the integration machinery independent of any specific +luminosity-function parameterization, catalog selection, or cosmology backend. + +These helpers are intentionally generic. Catalog completeness, LF-dependent +redshift densities, luminosity-density calculations, and selection-weighted +integrals can all call this module instead of duplicating magnitude-grid logic. +""" + +from __future__ import annotations + +from collections.abc import Callable + +import numpy as np + +from lfkit.utils.types import FloatArray, FloatInput, LuminosityFunction +from lfkit.utils.validators import validate_array + + +__all__ = [ + "integrated_number_density", + "lf_weighted_integral", + "selection_weighted_number_density", + "integrated_luminosity_density", + "mean_luminosity", + "cumulative_number_density", +] + + +def integrated_number_density( + z: FloatInput, + lf: LuminosityFunction, + *, + m_bright: FloatInput, + m_faint: FloatInput, + n_m: int = 512, +) -> FloatArray: + r"""Return finite-range number density from a luminosity function. + + This computes + + .. math:: + + n(z) = \int_{M_{\mathrm{bright}}(z)}^{M_{\mathrm{faint}}(z)} + \phi(M, z) \, dM. + + Magnitudes are ordered so that more negative values are brighter. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_bright: Bright absolute-magnitude bound. May be scalar or array-like. + m_faint: Faint absolute-magnitude bound. May be scalar or array-like. + n_m: Number of magnitude-grid points used for the integral. + + Returns: + NumPy array of number densities evaluated at ``z``. + """ + return _integrated_lf_between_bounds( + z, + lf, + m_lower=m_bright, + m_upper=m_faint, + n_m=n_m, + ) + + +def lf_weighted_integral( + z: FloatInput, + lf: LuminosityFunction, + *, + m_bright: FloatInput, + m_faint: FloatInput, + weight_fn: Callable[[FloatArray, FloatArray], FloatArray], + n_m: int = 512, +) -> FloatArray: + r"""Return a weighted luminosity-function integral. + + This computes + + .. math:: + + I(z) = \int_{M_{\mathrm{bright}}(z)}^{M_{\mathrm{faint}}(z)} + w(M, z)\,\phi(M, z)\,dM. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_bright: Bright absolute-magnitude bound. May be scalar or array-like. + m_faint: Faint absolute-magnitude bound. May be scalar or array-like. + weight_fn: Weight callable with signature ``weight_fn(M, z)``. Its + return values must be broadcastable to the magnitude-redshift grid. + n_m: Number of magnitude-grid points used for the integral. + + Returns: + NumPy array of weighted LF integrals evaluated at ``z``. + """ + return _integrated_lf_between_bounds( + z, + lf, + m_lower=m_bright, + m_upper=m_faint, + n_m=n_m, + weight_fn=weight_fn, + ) + + +def selection_weighted_number_density( + z: FloatInput, + lf: LuminosityFunction, + *, + m_bright: FloatInput, + m_faint: FloatInput, + selection_fn: Callable[[FloatArray, FloatArray], FloatArray], + n_m: int = 512, +) -> FloatArray: + r"""Return number density weighted by a selection function. + + This computes + + .. math:: + + n_{\mathrm{sel}}(z) = + \int_{M_{\mathrm{bright}}(z)}^{M_{\mathrm{faint}}(z)} + S(M, z)\,\phi(M, z)\,dM. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_bright: Bright absolute-magnitude bound. May be scalar or array-like. + m_faint: Faint absolute-magnitude bound. May be scalar or array-like. + selection_fn: Selection callable with signature ``selection_fn(M, z)``. + Values should usually lie between 0 and 1, although this function + only requires finite non-negative values. + n_m: Number of magnitude-grid points used for the integral. + + Returns: + NumPy array of selection-weighted number densities evaluated at ``z``. + """ + return lf_weighted_integral( + z, + lf, + m_bright=m_bright, + m_faint=m_faint, + weight_fn=selection_fn, + n_m=n_m, + ) + + +def integrated_luminosity_density( + z: FloatInput, + lf: LuminosityFunction, + *, + m_bright: FloatInput, + m_faint: FloatInput, + m_reference: float = 0.0, + n_m: int = 512, +) -> FloatArray: + r"""Return luminosity density from a luminosity function. + + This computes + + .. math:: + + \rho_L(z) = + \int_{M_{\mathrm{bright}}(z)}^{M_{\mathrm{faint}}(z)} + L(M)\,\phi(M, z)\,dM, + + using relative luminosities + + .. math:: + + L(M) / L_{\mathrm{ref}} = + 10^{-0.4(M - M_{\mathrm{ref}})}. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_bright: Bright absolute-magnitude bound. May be scalar or array-like. + m_faint: Faint absolute-magnitude bound. May be scalar or array-like. + m_reference: Reference absolute magnitude defining the luminosity unit. + n_m: Number of magnitude-grid points used for the integral. + + Returns: + NumPy array of luminosity densities in units of the reference + luminosity. + """ + if not np.isfinite(m_reference): + raise ValueError("m_reference must be finite.") + + def luminosity_weight( + absolute_mag: FloatArray, + _redshift: FloatArray, + ) -> FloatArray: + """Return relative luminosity weights for absolute magnitudes.""" + return np.asarray(10.0 ** (-0.4 * (absolute_mag - m_reference)), dtype=float) + + return lf_weighted_integral( + z, + lf, + m_bright=m_bright, + m_faint=m_faint, + weight_fn=luminosity_weight, + n_m=n_m, + ) + + +def mean_luminosity( + z: FloatInput, + lf: LuminosityFunction, + *, + m_bright: FloatInput, + m_faint: FloatInput, + m_reference: float = 0.0, + n_m: int = 512, +) -> FloatArray: + r"""Return mean luminosity over a finite magnitude range. + + This computes + + .. math:: + + \langle L \rangle(z) = + \frac{\int L(M)\,\phi(M, z)\,dM} + {\int \phi(M, z)\,dM}. + + The luminosity is returned in units of the reference luminosity defined by + ``m_reference``. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_bright: Bright absolute-magnitude bound. May be scalar or array-like. + m_faint: Faint absolute-magnitude bound. May be scalar or array-like. + m_reference: Reference absolute magnitude defining the luminosity unit. + n_m: Number of magnitude-grid points used for the integral. + + Returns: + NumPy array of mean luminosities. Entries are zero where the integrated + number density is zero. + """ + luminosity_density = integrated_luminosity_density( + z, + lf, + m_bright=m_bright, + m_faint=m_faint, + m_reference=m_reference, + n_m=n_m, + ) + number_density = integrated_number_density( + z, + lf, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + ) + + return _safe_divide(luminosity_density, number_density) + + +def cumulative_number_density( + z: FloatInput, + lf: LuminosityFunction, + *, + m_threshold: FloatInput, + m_bright: FloatInput, + m_faint: FloatInput, + brighter_than: bool = True, + n_m: int = 512, +) -> FloatArray: + r"""Return cumulative LF number density around a magnitude threshold. + + If ``brighter_than`` is True, this computes + + .. math:: + + n(< M_{\mathrm{thr}}, z) = + \int_{M_{\mathrm{bright}}}^{\min(M_{\mathrm{thr}}, M_{\mathrm{faint}})} + \phi(M, z)\,dM. + + If ``brighter_than`` is False, this computes + + .. math:: + + n(> M_{\mathrm{thr}}, z) = + \int_{\max(M_{\mathrm{thr}}, M_{\mathrm{bright}})}^{M_{\mathrm{faint}}} + \phi(M, z)\,dM. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_threshold: Absolute-magnitude threshold. May be scalar or array-like. + m_bright: Bright absolute-magnitude bound. May be scalar or array-like. + m_faint: Faint absolute-magnitude bound. May be scalar or array-like. + brighter_than: If True, integrate galaxies brighter than the threshold. + If False, integrate galaxies fainter than the threshold. + n_m: Number of magnitude-grid points used for the integral. + + Returns: + NumPy array of cumulative number densities evaluated at ``z``. + """ + z_arr = validate_array(z, name="z") + m_threshold_arr = validate_array(m_threshold, name="m_threshold") + m_bright_arr = validate_array(m_bright, name="m_bright") + m_faint_arr = validate_array(m_faint, name="m_faint") + + z_b, threshold_b, bright_b, faint_b = np.broadcast_arrays( + z_arr, + m_threshold_arr, + m_bright_arr, + m_faint_arr, + ) + + if brighter_than: + m_lower = bright_b + m_upper = np.minimum(threshold_b, faint_b) + else: + m_lower = np.maximum(threshold_b, bright_b) + m_upper = faint_b + + return _integrated_lf_between_bounds( + z_b, + lf, + m_lower=m_lower, + m_upper=m_upper, + n_m=n_m, + ) + + +def _integrated_lf_between_bounds( + z: FloatInput, + lf: LuminosityFunction, + *, + m_lower: FloatInput, + m_upper: FloatInput, + n_m: int, + weight_fn: Callable[[FloatArray, FloatArray], FloatArray] | None = None, +) -> FloatArray: + r"""Return finite-range LF integral between magnitude bounds.""" + z_arr = validate_array(z, name="z") + m_lower_arr = validate_array(m_lower, name="m_lower") + m_upper_arr = validate_array(m_upper, name="m_upper") + + _validate_integration_inputs( + z=z_arr, + m_lower=m_lower_arr, + m_upper=m_upper_arr, + n_m=n_m, + ) + + z_b, m_lower_b, m_upper_b = np.broadcast_arrays( + z_arr, + m_lower_arr, + m_upper_arr, + ) + + integral = np.zeros_like(z_b, dtype=float) + valid = m_upper_b > m_lower_b + + if not np.any(valid): + return np.asarray(integral, dtype=float) + + m_grid = _magnitude_grid( + m_lower=m_lower_b[valid], + m_upper=m_upper_b[valid], + n_m=n_m, + ) + z_grid = np.broadcast_to(z_b[valid][None, :], m_grid.shape) + + phi = _evaluate_lf_on_grid(lf, m_grid=m_grid, z_grid=z_grid) + + if weight_fn is not None: + weight = _evaluate_weight_on_grid( + weight_fn, + m_grid=m_grid, + z_grid=z_grid, + ) + phi = phi * weight + + integral[valid] = np.trapezoid(phi, x=m_grid, axis=0) + + return np.asarray(integral, dtype=float) + + +def _magnitude_grid( + *, + m_lower: FloatArray, + m_upper: FloatArray, + n_m: int, +) -> FloatArray: + r"""Return a magnitude grid for column-wise finite-range integration.""" + t = np.linspace(0.0, 1.0, int(n_m), dtype=float) + + return np.asarray( + m_lower[None, :] + t[:, None] * (m_upper[None, :] - m_lower[None, :]), + dtype=float, + ) + + +def _evaluate_lf_on_grid( + lf: LuminosityFunction, + *, + m_grid: FloatArray, + z_grid: FloatArray, +) -> FloatArray: + r"""Return LF values evaluated on a magnitude-redshift grid.""" + phi = np.asarray(lf(m_grid, z_grid), dtype=float) + + if phi.shape != m_grid.shape: + try: + phi = np.broadcast_to(phi, m_grid.shape) + except ValueError as exc: + raise ValueError( + "lf(M, z) must return values broadcastable to the shape " + "of the magnitude-redshift integration grid." + ) from exc + + if np.any(~np.isfinite(phi)): + raise ValueError("lf(M, z) returned non-finite values.") + + if np.any(phi < 0.0): + raise ValueError("lf(M, z) must be non-negative.") + + return np.asarray(phi, dtype=float) + + +def _evaluate_weight_on_grid( + weight_fn: Callable[[FloatArray, FloatArray], FloatArray], + *, + m_grid: FloatArray, + z_grid: FloatArray, +) -> FloatArray: + r"""Return finite non-negative weight values on a magnitude-redshift grid.""" + weight = np.asarray(weight_fn(m_grid, z_grid), dtype=float) + + if weight.shape != m_grid.shape: + try: + weight = np.broadcast_to(weight, m_grid.shape) + except ValueError as exc: + raise ValueError( + "weight_fn(M, z) must return values broadcastable to the shape " + "of the magnitude-redshift integration grid." + ) from exc + + if np.any(~np.isfinite(weight)): + raise ValueError("weight_fn(M, z) returned non-finite values.") + + if np.any(weight < 0.0): + raise ValueError("weight_fn(M, z) must be non-negative.") + + return np.asarray(weight, dtype=float) + + +def _validate_integration_inputs( + *, + z: FloatArray, + m_lower: FloatArray, + m_upper: FloatArray, + n_m: int, +) -> None: + r"""Validate inputs for finite-range LF integration.""" + if np.any(z < 0.0): + raise ValueError("Redshift z must be >= 0.") + + if np.any(~np.isfinite(m_lower)): + raise ValueError("m_lower must contain only finite values.") + + if np.any(~np.isfinite(m_upper)): + raise ValueError("m_upper must contain only finite values.") + + if n_m < 2: + raise ValueError("n_m must be at least 2.") + + +def _safe_divide( + numerator: FloatArray, + denominator: FloatArray, +) -> FloatArray: + r"""Return numerator / denominator with zero output for zero denominator.""" + numerator_arr = np.asarray(numerator, dtype=float) + denominator_arr = np.asarray(denominator, dtype=float) + + with np.errstate(divide="ignore", invalid="ignore"): + result = np.divide( + numerator_arr, + denominator_arr, + out=np.zeros_like(numerator_arr, dtype=float), + where=denominator_arr > 0.0, + ) + + return np.asarray(result, dtype=float) diff --git a/src/lfkit/photometry/lf_redshift_density.py b/src/lfkit/photometry/lf_redshift_density.py new file mode 100644 index 0000000..9a619b1 --- /dev/null +++ b/src/lfkit/photometry/lf_redshift_density.py @@ -0,0 +1,238 @@ +r"""Luminosity-function redshift-density utilities. + +This module provides helpers for converting a luminosity-function callable into +an LF-integrated redshift-density curve. + +The core operation is + + n_lf(z) = int phi(M, z) dM + +over the observable absolute-magnitude range implied by an apparent-magnitude +limit. A second helper multiplies this LF-integrated density by a user-supplied +redshift or volume weight. + +The cosmology-dependent distance and volume pieces are supplied as callables. +This keeps the interface independent of CCL, Astropy, or any other cosmology +backend, which is useful for downstream packages such as Binny. + +Magnitude corrections are supplied as scalar or array-like values evaluated on +the same redshift grid. +""" + +from __future__ import annotations + +from collections.abc import Callable + +import numpy as np + +from lfkit.photometry.lf_integrals import integrated_number_density +from lfkit.photometry.magnitudes import absolute_magnitude_from_luminosity_distance +from lfkit.utils.evaluators import ( + evaluate_non_negative_redshift_callable, + evaluate_positive_redshift_callable, +) +from lfkit.utils.types import FloatArray, FloatInput, LuminosityFunction +from lfkit.utils.validators import validate_array + + +__all__ = [ + "lf_integrated_number_density", + "lf_weighted_redshift_density", +] + + +def lf_integrated_number_density( + z: FloatInput, + lf: LuminosityFunction, + *, + m_lim: float, + m_bright: float, + n_m: int = 512, + luminosity_distance_mpc_fn: Callable[[FloatArray], FloatArray], + k_correction: FloatInput | None = None, + evolution_correction: FloatInput | None = None, +) -> FloatArray: + r"""Return LF-integrated number density as a function of redshift. + + This computes + + .. math:: + + n_{\mathrm{LF}}(z) = + \int_{M_{\mathrm{bright}}}^{M_{\mathrm{lim}}(z)} + \phi(M, z)\,dM, + + where ``M_lim(z)`` is the absolute-magnitude limit implied by the apparent + magnitude cut ``m_lim``. + + The magnitude conversion follows + + .. math:: + + M_{\mathrm{lim}}(z) = + m_{\mathrm{lim}} - \mu(z) - K(z) + E(z), + + where ``mu`` is the distance modulus, ``K`` is the k-correction, and ``E`` + is the evolution correction. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_lim: Apparent magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration bound. + n_m: Number of magnitude-grid points used for the integral. + luminosity_distance_mpc_fn: Callable returning luminosity distance in + Mpc as a function of redshift. + k_correction: Optional scalar or array-like k-correction values. + evolution_correction: Optional scalar or array-like evolution-correction + values. + + Returns: + NumPy array of LF-integrated number densities evaluated at ``z``. + """ + z_arr = validate_array(z, name="z") + + if np.any(z_arr < 0.0): + raise ValueError("Redshift z must be >= 0.") + + if not np.isfinite(m_lim): + raise ValueError("m_lim must be finite.") + + if not np.isfinite(m_bright): + raise ValueError("m_bright must be finite.") + + luminosity_distance = evaluate_positive_redshift_callable( + luminosity_distance_mpc_fn, + z_arr, + name="luminosity_distance_mpc_fn", + ) + + k_correction_arr = _optional_correction_array( + k_correction, + z_arr, + name="k_correction", + ) + + evolution_correction_arr = _optional_correction_array( + evolution_correction, + z_arr, + name="evolution_correction", + ) + + m_faint = absolute_magnitude_from_luminosity_distance( + m_lim, + luminosity_distance, + k_correction=k_correction_arr, + e_correction=evolution_correction_arr, + ) + + return integrated_number_density( + z_arr, + lf, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + ) + + +def lf_weighted_redshift_density( + z: FloatInput, + lf: LuminosityFunction, + *, + m_lim: float, + m_bright: float, + n_m: int = 512, + luminosity_distance_mpc_fn: Callable[[FloatArray], FloatArray], + volume_weight_fn: Callable[[FloatArray], FloatArray], + k_correction: FloatInput | None = None, + evolution_correction: FloatInput | None = None, + normalize: bool = True, +) -> FloatArray: + r"""Return an LF-weighted redshift-density curve. + + This computes + + .. math:: + + n(z) \propto W(z) + \int_{M_{\mathrm{bright}}}^{M_{\mathrm{lim}}(z)} + \phi(M, z)\,dM, + + where ``W(z)`` is supplied by ``volume_weight_fn``. + + Args: + z: Redshift value or array-like of redshift values. + lf: Luminosity-function callable with signature ``lf(M, z)``. + m_lim: Apparent magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration bound. + n_m: Number of magnitude-grid points used for the integral. + luminosity_distance_mpc_fn: Callable returning luminosity distance in + Mpc as a function of redshift. + volume_weight_fn: Callable returning the redshift or volume weight. + k_correction: Optional scalar or array-like k-correction values. + evolution_correction: Optional scalar or array-like evolution-correction + values. + normalize: If True, normalize the returned curve to integrate to one + over ``z``. + + Returns: + NumPy array of LF-weighted redshift-density values. + """ + z_arr = validate_array(z, name="z") + + if np.any(z_arr < 0.0): + raise ValueError("Redshift z must be >= 0.") + + n_lf = lf_integrated_number_density( + z_arr, + lf, + m_lim=m_lim, + m_bright=m_bright, + n_m=n_m, + luminosity_distance_mpc_fn=luminosity_distance_mpc_fn, + k_correction=k_correction, + evolution_correction=evolution_correction, + ) + + volume_weight = evaluate_non_negative_redshift_callable( + volume_weight_fn, + z_arr, + name="volume_weight_fn", + ) + + nz = n_lf * volume_weight + + if normalize: + norm = np.trapezoid(nz, x=z_arr) + + if not np.isfinite(norm) or norm <= 0.0: + raise ValueError( + "Cannot normalize LF-weighted redshift density with " + "non-positive integral." + ) + + nz = nz / norm + + return np.asarray(nz, dtype=float) + + +def _optional_correction_array( + correction: FloatInput | None, + z: FloatArray, + *, + name: str, +) -> FloatArray: + """Return an optional scalar or array correction broadcast to redshift.""" + if correction is None: + return np.zeros_like(z, dtype=float) + + correction_arr = validate_array(correction, name=name) + + try: + correction_b = np.broadcast_to(correction_arr, z.shape) + except ValueError as exc: + raise ValueError( + f"{name} must be scalar or broadcastable to the shape of z." + ) from exc + + return np.asarray(correction_b, dtype=float) diff --git a/src/lfkit/photometry/magnitudes.py b/src/lfkit/photometry/magnitudes.py index e24cba6..7e46d2a 100644 --- a/src/lfkit/photometry/magnitudes.py +++ b/src/lfkit/photometry/magnitudes.py @@ -20,12 +20,15 @@ from lfkit.cosmo.cosmology import distance_modulus from lfkit.utils.types import Cosmology, FloatArray, FloatInput +from lfkit.utils.validators import validate_luminosity_distance -__all__ = ( +__all__ = [ "total_magnitude_correction", "absolute_magnitude", + "absolute_magnitude_from_luminosity_distance", "apparent_magnitude", -) + "apparent_magnitude_from_luminosity_distance", +] def total_magnitude_correction( @@ -133,3 +136,79 @@ def apparent_magnitude( ) return np.asarray(M + mu + correction, dtype=np.float64) + + +def absolute_magnitude_from_luminosity_distance( + apparent_mag: FloatInput, + luminosity_distance_mpc: FloatInput, + *, + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: + """Convert apparent magnitude to absolute magnitude from luminosity distance. + + The convention used is + + ``M = m - mu - K + E``, + + with + + ``mu = 5 log10(d_L / Mpc) + 25``. + + Args: + apparent_mag: Apparent magnitude value(s). + luminosity_distance_mpc: Luminosity distance value(s) in Mpc. + k_correction: Optional k-correction term(s). + e_correction: Optional evolution-correction term(s). + + Returns: + NumPy array of absolute magnitudes. + """ + m = np.asarray(apparent_mag, dtype=float) + d_l = validate_luminosity_distance(luminosity_distance_mpc) + + mu = 5.0 * np.log10(d_l) + 25.0 + correction = total_magnitude_correction( + k_correction=k_correction, + e_correction=e_correction, + ) + + return np.asarray(m - mu - correction, dtype=np.float64) + + +def apparent_magnitude_from_luminosity_distance( + absolute_mag: FloatInput, + luminosity_distance_mpc: FloatInput, + *, + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: + """Convert absolute magnitude to apparent magnitude from luminosity distance. + + The convention used is + + ``m = M + mu + K - E``, + + with + + ``mu = 5 log10(d_L / Mpc) + 25``. + + Args: + absolute_mag: Absolute magnitude value(s). + luminosity_distance_mpc: Luminosity distance value(s) in Mpc. + k_correction: Optional k-correction term(s). + e_correction: Optional evolution-correction term(s). + + Returns: + NumPy array of apparent magnitudes. + """ + M = np.asarray(absolute_mag, dtype=float) + d_l = validate_luminosity_distance(luminosity_distance_mpc) + + mu = 5.0 * np.log10(d_l) + 25.0 + correction = total_magnitude_correction( + k_correction=k_correction, + e_correction=e_correction, + ) + + return np.asarray(M + mu + correction, dtype=np.float64) diff --git a/src/lfkit/utils/evaluators.py b/src/lfkit/utils/evaluators.py new file mode 100644 index 0000000..72b3f3e --- /dev/null +++ b/src/lfkit/utils/evaluators.py @@ -0,0 +1,104 @@ +"""Callable evaluation utilities.""" + +from __future__ import annotations + +from collections.abc import Callable + +import numpy as np + +from lfkit.utils.types import FloatArray + + +__all__ = [ + "evaluate_non_negative_redshift_callable", + "evaluate_optional_redshift_callable", + "evaluate_positive_redshift_callable", +] + + +def evaluate_optional_redshift_callable( + fn: Callable[[FloatArray], FloatArray] | None, + z: FloatArray, + *, + name: str, +) -> FloatArray | None: + """Evaluate an optional redshift-dependent callable. + + Args: + fn: Callable to evaluate. If None, None is returned. + z: Redshift array passed to the callable. + name: Name used in error messages. + + Returns: + Callable values with the same shape as ``z``, or None. + """ + if fn is None: + return None + + return _evaluate_redshift_callable(fn, z, name=name) + + +def evaluate_positive_redshift_callable( + fn: Callable[[FloatArray], FloatArray], + z: FloatArray, + *, + name: str, +) -> FloatArray: + """Evaluate a redshift-dependent callable that must be positive. + + Args: + fn: Callable to evaluate. + z: Redshift array passed to the callable. + name: Name used in error messages. + + Returns: + Positive callable values with the same shape as ``z``. + """ + values = _evaluate_redshift_callable(fn, z, name=name) + + if np.any(values <= 0.0): + raise ValueError(f"{name}(z) must return positive values.") + + return values + + +def evaluate_non_negative_redshift_callable( + fn: Callable[[FloatArray], FloatArray], + z: FloatArray, + *, + name: str, +) -> FloatArray: + """Evaluate a redshift-dependent callable that must be non-negative. + + Args: + fn: Callable to evaluate. + z: Redshift array passed to the callable. + name: Name used in error messages. + + Returns: + Non-negative callable values with the same shape as ``z``. + """ + values = _evaluate_redshift_callable(fn, z, name=name) + + if np.any(values < 0.0): + raise ValueError(f"{name}(z) must return non-negative values.") + + return values + + +def _evaluate_redshift_callable( + fn: Callable[[FloatArray], FloatArray], + z: FloatArray, + *, + name: str, +) -> FloatArray: + """Evaluate a redshift callable and validate shape and finite values.""" + values = np.asarray(fn(z), dtype=float) + + if values.shape != z.shape: + raise ValueError(f"{name}(z) must return an array with the same shape as z.") + + if np.any(~np.isfinite(values)): + raise ValueError(f"{name}(z) returned non-finite values.") + + return np.asarray(values, dtype=float) diff --git a/src/lfkit/utils/validators.py b/src/lfkit/utils/validators.py index 874d87f..cb4c85e 100644 --- a/src/lfkit/utils/validators.py +++ b/src/lfkit/utils/validators.py @@ -8,6 +8,7 @@ __all__ = [ "validate_array", + "validate_luminosity_distance", "validate_magnitude_range", ] @@ -30,6 +31,22 @@ def validate_array( return np.asarray(arr, dtype=np.float64) +def validate_luminosity_distance( + luminosity_distance_mpc: FloatInput, +) -> FloatArray: + """Return finite positive luminosity distances in Mpc.""" + distance = validate_array( + luminosity_distance_mpc, + name="luminosity_distance_mpc", + allow_negative=False, + ) + + if np.any(distance <= 0.0): + raise ValueError("luminosity_distance_mpc must contain positive values.") + + return distance + + def validate_magnitude_range( *, m_bright: float, diff --git a/tests/test_api_lumfunc.py b/tests/test_api_lumfunc.py index 1f830ba..36e8f4a 100644 --- a/tests/test_api_lumfunc.py +++ b/tests/test_api_lumfunc.py @@ -530,3 +530,800 @@ def test_catalog_and_out_of_catalog_fractions_sum_to_one(monkeypatch): ) assert np.allclose(completeness + missing, 1.0) + + +def test_absolute_magnitude_forwards_corrections(monkeypatch): + """Tests that absolute_magnitude forwards correction arrays.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + apparent_mag = np.array([22.0, 24.0]) + corrections = cast(Corrections, DummyCorrections()) + + def fake_absolute_magnitude( + received_cosmo, + received_z, + received_apparent_mag, + *, + h, + k_correction, + e_correction, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert np.allclose(received_apparent_mag, apparent_mag) + assert h == 0.7 + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(e_correction, [-0.9, -0.5]) + return np.array([-19.0, -20.0]) + + monkeypatch.setattr(lf_api, "absolute_magnitude", fake_absolute_magnitude) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.absolute_magnitude( + cosmo_obj, + z, + apparent_mag, + h=0.7, + corrections=corrections, + ) + + assert np.allclose(result, [-19.0, -20.0]) + + +def test_apparent_magnitude_forwards_corrections(monkeypatch): + """Tests that apparent_magnitude forwards correction arrays.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + absolute_mag = np.array([-19.0, -20.0]) + corrections = cast(Corrections, DummyCorrections()) + + def fake_apparent_magnitude( + received_cosmo, + received_z, + received_absolute_mag, + *, + h, + k_correction, + e_correction, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert np.allclose(received_absolute_mag, absolute_mag) + assert h == 0.7 + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(e_correction, [-0.9, -0.5]) + return np.array([22.0, 24.0]) + + monkeypatch.setattr(lf_api, "apparent_magnitude", fake_apparent_magnitude) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.apparent_magnitude( + cosmo_obj, + z, + absolute_mag, + h=0.7, + corrections=corrections, + ) + + assert np.allclose(result, [22.0, 24.0]) + + +def test_absolute_magnitude_from_luminosity_distance_forwards_corrections(monkeypatch): + """Tests that absolute magnitude from distance forwards correction arrays.""" + z = np.array([0.1, 0.5]) + apparent_mag = np.array([22.0, 24.0]) + luminosity_distance_mpc = np.array([500.0, 1500.0]) + corrections = cast(Corrections, DummyCorrections()) + + def fake_absolute_magnitude_from_luminosity_distance( + received_apparent_mag, + received_luminosity_distance_mpc, + *, + k_correction, + e_correction, + ): + assert np.allclose(received_apparent_mag, apparent_mag) + assert np.allclose(received_luminosity_distance_mpc, luminosity_distance_mpc) + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(e_correction, [-0.9, -0.5]) + return np.array([-18.0, -21.0]) + + monkeypatch.setattr( + lf_api, + "absolute_magnitude_from_luminosity_distance", + fake_absolute_magnitude_from_luminosity_distance, + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.absolute_magnitude_from_luminosity_distance( + apparent_mag, + luminosity_distance_mpc, + z=z, + corrections=corrections, + ) + + assert np.allclose(result, [-18.0, -21.0]) + + +def test_apparent_magnitude_from_luminosity_distance_forwards_corrections(monkeypatch): + """Tests that apparent magnitude from distance forwards correction arrays.""" + z = np.array([0.1, 0.5]) + absolute_mag = np.array([-18.0, -21.0]) + luminosity_distance_mpc = np.array([500.0, 1500.0]) + corrections = cast(Corrections, DummyCorrections()) + + def fake_apparent_magnitude_from_luminosity_distance( + received_absolute_mag, + received_luminosity_distance_mpc, + *, + k_correction, + e_correction, + ): + assert np.allclose(received_absolute_mag, absolute_mag) + assert np.allclose(received_luminosity_distance_mpc, luminosity_distance_mpc) + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(e_correction, [-0.9, -0.5]) + return np.array([22.0, 24.0]) + + monkeypatch.setattr( + lf_api, + "apparent_magnitude_from_luminosity_distance", + fake_apparent_magnitude_from_luminosity_distance, + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.apparent_magnitude_from_luminosity_distance( + absolute_mag, + luminosity_distance_mpc, + z=z, + corrections=corrections, + ) + + assert np.allclose(result, [22.0, 24.0]) + + +def test_absolute_magnitude_limit_forwards_corrections(monkeypatch): + """Tests that absolute_magnitude_limit forwards correction arrays.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + corrections = cast(Corrections, DummyCorrections()) + + def fake_absolute_magnitude_limit( + received_cosmo, + received_z, + *, + m_lim, + h, + k_correction, + e_correction, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert m_lim == 24.5 + assert h == 0.7 + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(e_correction, [-0.9, -0.5]) + return np.array([-18.0, -20.0]) + + monkeypatch.setattr(lf_api, "absolute_magnitude_limit", fake_absolute_magnitude_limit) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.absolute_magnitude_limit( + cosmo_obj, + z, + m_lim=24.5, + h=0.7, + corrections=corrections, + ) + + assert np.allclose(result, [-18.0, -20.0]) + + +def test_phi_from_m_dispatches_evolving_schechter_from_m(monkeypatch): + """Tests that phi_from_m dispatches to evolving Schechter from apparent magnitude.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + apparent_mag = np.array([22.0, 24.0]) + + def fake_schechter_evolving_from_m( + received_cosmo, + received_z, + received_apparent_mag, + *, + h, + k_correction, + e_correction, + phi_model, + phi_kwargs, + m_star_model, + m_star_kwargs, + alpha_model, + alpha_kwargs, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert np.allclose(received_apparent_mag, apparent_mag) + assert h is None + assert k_correction is None + assert e_correction is None + assert phi_model == "constant" + assert phi_kwargs == {"value": 1.0e-3} + assert m_star_model == "constant" + assert m_star_kwargs == {"value": -20.5} + assert alpha_model == "constant" + assert alpha_kwargs == {"value": -1.2} + return np.array([9.0, 10.0]) + + monkeypatch.setattr( + lf_api, + "schechter_evolving_from_m", + fake_schechter_evolving_from_m, + ) + + lf = LuminosityFunction.evolving_schechter( + phi_model="constant", + phi_kwargs={"value": 1.0e-3}, + m_star_model="constant", + m_star_kwargs={"value": -20.5}, + alpha_model="constant", + alpha_kwargs={"value": -1.2}, + ) + + result = lf.phi_from_m(cosmo_obj, z, apparent_mag) + + assert np.allclose(result, [9.0, 10.0]) + + +def test_phi_from_m_dispatches_double_schechter_from_m(monkeypatch): + """Tests that phi_from_m dispatches to double Schechter from apparent magnitude.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + apparent_mag = np.array([22.0, 24.0]) + + def fake_schechter_double_from_m( + received_cosmo, + received_z, + received_apparent_mag, + *, + h, + k_correction, + e_correction, + phi_star, + m_star, + alpha, + beta, + m_transition, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert np.allclose(received_apparent_mag, apparent_mag) + assert h is None + assert k_correction is None + assert e_correction is None + assert phi_star == 1.0e-3 + assert m_star == -20.5 + assert alpha == -1.2 + assert beta == -2.0 + assert m_transition == -19.0 + return np.array([11.0, 12.0]) + + monkeypatch.setattr( + lf_api, + "schechter_double_from_m", + fake_schechter_double_from_m, + ) + + lf = LuminosityFunction.double_schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + beta=-2.0, + m_transition=-19.0, + ) + + result = lf.phi_from_m(cosmo_obj, z, apparent_mag) + + assert np.allclose(result, [11.0, 12.0]) + + +def test_phi_from_m_raises_for_unsupported_model(): + """Tests that phi_from_m rejects unsupported LF models.""" + lf = LuminosityFunction(model="bad_model", parameters={}) + + with pytest.raises(ValueError, match="Unsupported luminosity-function model"): + lf.phi_from_m(object(), np.array([0.1]), np.array([22.0])) + + +def test_lf_weighted_integral_forwards_to_module_function(monkeypatch): + """Tests that lf_weighted_integral forwards inputs to the module function.""" + z = np.array([0.1, 0.5]) + + def weight_fn(absolute_mag, redshift): + return np.ones_like(absolute_mag, dtype=float) + + def fake_lf_weighted_integral( + received_z, + lf_callable, + *, + m_bright, + m_faint, + weight_fn: object, + n_m, + ): + mag = np.array([-21.0, -20.0]) + assert np.allclose(received_z, z) + assert m_bright == -24.0 + assert m_faint == -18.0 + assert n_m == 64 + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([3.0, 4.0]) + + monkeypatch.setattr(lf_api, "lf_weighted_integral", fake_lf_weighted_integral) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + monkeypatch.setattr( + lf, + "phi", + lambda absolute_mag, z=None: np.full_like(absolute_mag, 2.0), + ) + + result = lf.lf_weighted_integral( + z, + m_bright=-24.0, + m_faint=-18.0, + weight_fn=weight_fn, + n_m=64, + ) + + assert np.allclose(result, [3.0, 4.0]) + + +def test_selection_weighted_number_density_forwards_to_module_function(monkeypatch): + """Tests that selection_weighted_number_density forwards inputs.""" + z = np.array([0.1, 0.5]) + + def selection_fn(absolute_mag, redshift): + return np.ones_like(absolute_mag, dtype=float) + + def fake_selection_weighted_number_density( + received_z, + lf_callable, + *, + m_bright, + m_faint, + selection_fn: object, + n_m, + ): + mag = np.array([-21.0, -20.0]) + assert np.allclose(received_z, z) + assert m_bright == -24.0 + assert m_faint == -18.0 + assert n_m == 64 + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([5.0, 6.0]) + + monkeypatch.setattr( + lf_api, + "selection_weighted_number_density", + fake_selection_weighted_number_density, + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + monkeypatch.setattr( + lf, + "phi", + lambda absolute_mag, z=None: np.full_like(absolute_mag, 2.0), + ) + + result = lf.selection_weighted_number_density( + z, + m_bright=-24.0, + m_faint=-18.0, + selection_fn=selection_fn, + n_m=64, + ) + + assert np.allclose(result, [5.0, 6.0]) + + +def test_integrated_luminosity_density_forwards_to_module_function(monkeypatch): + """Tests that integrated_luminosity_density forwards inputs.""" + z = np.array([0.1, 0.5]) + + def fake_integrated_luminosity_density( + received_z, + lf_callable, + *, + m_bright, + m_faint, + m_reference, + n_m, + ): + mag = np.array([-21.0, -20.0]) + assert np.allclose(received_z, z) + assert m_bright == -24.0 + assert m_faint == -18.0 + assert m_reference == -20.0 + assert n_m == 64 + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([7.0, 8.0]) + + monkeypatch.setattr( + lf_api, + "integrated_luminosity_density", + fake_integrated_luminosity_density, + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + monkeypatch.setattr( + lf, + "phi", + lambda absolute_mag, z=None: np.full_like(absolute_mag, 2.0), + ) + + result = lf.integrated_luminosity_density( + z, + m_bright=-24.0, + m_faint=-18.0, + m_reference=-20.0, + n_m=64, + ) + + assert np.allclose(result, [7.0, 8.0]) + + +def test_mean_luminosity_forwards_to_module_function(monkeypatch): + """Tests that mean_luminosity forwards inputs.""" + z = np.array([0.1, 0.5]) + + def fake_mean_luminosity( + received_z, + lf_callable, + *, + m_bright, + m_faint, + m_reference, + n_m, + ): + mag = np.array([-21.0, -20.0]) + assert np.allclose(received_z, z) + assert m_bright == -24.0 + assert m_faint == -18.0 + assert m_reference == -20.0 + assert n_m == 64 + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([9.0, 10.0]) + + monkeypatch.setattr(lf_api, "mean_luminosity", fake_mean_luminosity) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + monkeypatch.setattr( + lf, + "phi", + lambda absolute_mag, z=None: np.full_like(absolute_mag, 2.0), + ) + + result = lf.mean_luminosity( + z, + m_bright=-24.0, + m_faint=-18.0, + m_reference=-20.0, + n_m=64, + ) + + assert np.allclose(result, [9.0, 10.0]) + + +def test_cumulative_number_density_forwards_to_module_function(monkeypatch): + """Tests that cumulative_number_density forwards inputs.""" + z = np.array([0.1, 0.5]) + + def fake_cumulative_number_density( + received_z, + lf_callable, + *, + m_threshold, + m_bright, + m_faint, + brighter_than, + n_m, + ): + mag = np.array([-21.0, -20.0]) + assert np.allclose(received_z, z) + assert m_threshold == -20.0 + assert m_bright == -24.0 + assert m_faint == -18.0 + assert brighter_than is False + assert n_m == 64 + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([11.0, 12.0]) + + monkeypatch.setattr( + lf_api, + "cumulative_number_density", + fake_cumulative_number_density, + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + monkeypatch.setattr( + lf, + "phi", + lambda absolute_mag, z=None: np.full_like(absolute_mag, 2.0), + ) + + result = lf.cumulative_number_density( + z, + m_threshold=-20.0, + m_bright=-24.0, + m_faint=-18.0, + brighter_than=False, + n_m=64, + ) + + assert np.allclose(result, [11.0, 12.0]) + + +def test_lf_integrated_number_density_forwards_corrections(monkeypatch): + """Tests that apparent-limit LF number density forwards corrections.""" + z = np.array([0.1, 0.5]) + corrections = cast(Corrections, DummyCorrections()) + + def luminosity_distance_mpc_fn(redshift): + return 1000.0 * np.asarray(redshift, dtype=float) + + def fake_lf_integrated_number_density( + received_z, + lf_callable, + *, + m_lim, + m_bright, + n_m, + luminosity_distance_mpc_fn, + k_correction, + evolution_correction, + ): + mag = np.array([-21.0, -20.0]) + assert np.allclose(received_z, z) + assert m_lim == 24.5 + assert m_bright == -24.0 + assert n_m == 64 + assert np.allclose(luminosity_distance_mpc_fn(z), [100.0, 500.0]) + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(evolution_correction, [-0.9, -0.5]) + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([13.0, 14.0]) + + monkeypatch.setattr( + lf_api, + "lf_integrated_number_density", + fake_lf_integrated_number_density, + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + monkeypatch.setattr( + lf, + "phi", + lambda absolute_mag, z=None: np.full_like(absolute_mag, 2.0), + ) + + result = lf.lf_integrated_number_density( + z, + m_lim=24.5, + m_bright=-24.0, + n_m=64, + luminosity_distance_mpc_fn=luminosity_distance_mpc_fn, + corrections=corrections, + ) + + assert np.allclose(result, [13.0, 14.0]) + + +def test_lf_weighted_redshift_density_forwards_corrections(monkeypatch): + """Tests that LF-weighted redshift density forwards corrections.""" + z = np.array([0.1, 0.5]) + corrections = cast(Corrections, DummyCorrections()) + + def luminosity_distance_mpc_fn(redshift): + return 1000.0 * np.asarray(redshift, dtype=float) + + def volume_weight_fn(redshift): + return 2.0 * np.asarray(redshift, dtype=float) + + def fake_lf_weighted_redshift_density( + received_z, + lf_callable, + *, + m_lim, + m_bright, + n_m, + luminosity_distance_mpc_fn, + volume_weight_fn, + k_correction, + evolution_correction, + normalize, + ): + mag = np.array([-21.0, -20.0]) + assert np.allclose(received_z, z) + assert m_lim == 24.5 + assert m_bright == -24.0 + assert n_m == 64 + assert np.allclose(luminosity_distance_mpc_fn(z), [100.0, 500.0]) + assert np.allclose(volume_weight_fn(z), [0.2, 1.0]) + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(evolution_correction, [-0.9, -0.5]) + assert normalize is False + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([15.0, 16.0]) + + monkeypatch.setattr( + lf_api, + "lf_weighted_redshift_density", + fake_lf_weighted_redshift_density, + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + monkeypatch.setattr( + lf, + "phi", + lambda absolute_mag, z=None: np.full_like(absolute_mag, 2.0), + ) + + result = lf.lf_weighted_redshift_density( + z, + m_lim=24.5, + m_bright=-24.0, + n_m=64, + luminosity_distance_mpc_fn=luminosity_distance_mpc_fn, + volume_weight_fn=volume_weight_fn, + corrections=corrections, + normalize=False, + ) + + assert np.allclose(result, [15.0, 16.0]) + + +def test_available_parameter_models_delegates_to_module_function(monkeypatch): + """Tests that available_parameter_models delegates to the registry helper.""" + expected = { + "phi_star": ["constant", "linear_p"], + "m_star": ["constant", "linear_q"], + "alpha": ["constant"], + } + + monkeypatch.setattr( + lf_api, + "available_lf_parameter_models", + lambda: expected, + ) + + result = LuminosityFunction.available_parameter_models() + + assert result == expected + + +def test_register_phi_star_model_delegates_to_module_function(monkeypatch): + """Tests that register_phi_star_model delegates to the registry helper.""" + captured = {} + + def model(z, *, value): + return np.full_like(np.asarray(z, dtype=float), value) + + def fake_register(name, received_model, *, overwrite): + captured["name"] = name + captured["model"] = received_model + captured["overwrite"] = overwrite + + monkeypatch.setattr(lf_api, "register_phi_star_model", fake_register) + + LuminosityFunction.register_phi_star_model( + "test_phi", + model, + overwrite=True, + ) + + assert captured["name"] == "test_phi" + assert captured["model"] is model + assert captured["overwrite"] is True + + +def test_register_m_star_model_delegates_to_module_function(monkeypatch): + """Tests that register_m_star_model delegates to the registry helper.""" + captured = {} + + def model(z, *, value): + return np.full_like(np.asarray(z, dtype=float), value) + + def fake_register(name, received_model, *, overwrite): + captured["name"] = name + captured["model"] = received_model + captured["overwrite"] = overwrite + + monkeypatch.setattr(lf_api, "register_m_star_model", fake_register) + + LuminosityFunction.register_m_star_model( + "test_m_star", + model, + overwrite=True, + ) + + assert captured["name"] == "test_m_star" + assert captured["model"] is model + assert captured["overwrite"] is True + + +def test_register_alpha_model_delegates_to_module_function(monkeypatch): + """Tests that register_alpha_model delegates to the registry helper.""" + captured = {} + + def model(z, *, value): + return np.full_like(np.asarray(z, dtype=float), value) + + def fake_register(name, received_model, *, overwrite): + captured["name"] = name + captured["model"] = received_model + captured["overwrite"] = overwrite + + monkeypatch.setattr(lf_api, "register_alpha_model", fake_register) + + LuminosityFunction.register_alpha_model( + "test_alpha", + model, + overwrite=True, + ) + + assert captured["name"] == "test_alpha" + assert captured["model"] is model + assert captured["overwrite"] is True diff --git a/tests/test_photometry_completeness.py b/tests/test_photometry_completeness.py index 944b70d..4acc5d0 100644 --- a/tests/test_photometry_completeness.py +++ b/tests/test_photometry_completeness.py @@ -1,4 +1,4 @@ -"""Unit tests for the ``lfkit.photometry.catalog_completeness.py``""" +"""Unit tests for the ``lfkit.photometry.catalog_completeness.py`` module.""" import numpy as np import pytest @@ -16,7 +16,9 @@ def double_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: return 2.0 * np.ones_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) -def test_absolute_magnitude_limit_calls_absolute_magnitude(monkeypatch: pytest.MonkeyPatch) -> None: +def test_absolute_magnitude_limit_calls_absolute_magnitude( + monkeypatch: pytest.MonkeyPatch, +) -> None: """Tests that absolute magnitude limits are delegated to magnitude conversion.""" calls = {} @@ -70,129 +72,6 @@ def test_absolute_magnitude_limit_rejects_nonfinite_m_lim() -> None: cc.absolute_magnitude_limit(object(), [0.1, 0.2], m_lim=np.inf) -def test_integrated_number_density_integrates_constant_lf() -> None: - """Tests that finite-range integration returns the expected width.""" - result = cc.integrated_number_density( - [0.1, 0.2], - constant_lf, - m_bright=-24.0, - m_faint=-18.0, - n_m=64, - ) - - np.testing.assert_allclose(result, np.array([6.0, 6.0])) - - -def test_integrated_number_density_integrates_lf_amplitude() -> None: - """Tests that finite-range integration preserves LF amplitude.""" - result = cc.integrated_number_density( - [0.1, 0.2], - double_lf, - m_bright=-24.0, - m_faint=-18.0, - n_m=64, - ) - - np.testing.assert_allclose(result, np.array([12.0, 12.0])) - - -def test_integrated_number_density_supports_array_bounds() -> None: - """Tests that finite-range integration supports redshift-dependent bounds.""" - result = cc.integrated_number_density( - [0.1, 0.2, 0.3], - constant_lf, - m_bright=np.array([-24.0, -23.0, -22.0]), - m_faint=np.array([-18.0, -18.0, -18.0]), - n_m=64, - ) - - np.testing.assert_allclose(result, np.array([6.0, 5.0, 4.0])) - - -def test_integrated_number_density_returns_zero_for_empty_ranges() -> None: - """Tests that empty magnitude ranges return zero density.""" - result = cc.integrated_number_density( - [0.1, 0.2], - constant_lf, - m_bright=np.array([-18.0, -20.0]), - m_faint=np.array([-20.0, -20.0]), - n_m=64, - ) - - np.testing.assert_allclose(result, np.array([0.0, 0.0])) - - -def test_integrated_number_density_rejects_negative_redshift() -> None: - """Tests that finite-range integration rejects negative redshifts.""" - with pytest.raises(ValueError, match="Redshift z must be >= 0"): - cc.integrated_number_density( - [-0.1, 0.2], - constant_lf, - m_bright=-24.0, - m_faint=-18.0, - ) - - -def test_integrated_number_density_rejects_small_magnitude_grid() -> None: - """Tests that finite-range integration requires at least two grid points.""" - with pytest.raises(ValueError, match="n_m must be at least 2"): - cc.integrated_number_density( - [0.1, 0.2], - constant_lf, - m_bright=-24.0, - m_faint=-18.0, - n_m=1, - ) - - -def test_integrated_number_density_rejects_nonfinite_lf_values() -> None: - """Tests that non-finite luminosity-function values are rejected.""" - - def bad_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: - return np.full_like(m_abs, np.nan, dtype=float) - - with pytest.raises(ValueError, match="lf\\(M, z\\) returned non-finite values"): - cc.integrated_number_density( - [0.1, 0.2], - bad_lf, - m_bright=-24.0, - m_faint=-18.0, - ) - - -def test_integrated_number_density_rejects_negative_lf_values() -> None: - """Tests that negative luminosity-function values are rejected.""" - - def bad_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: - return -np.ones_like(m_abs, dtype=float) - - with pytest.raises(ValueError, match="lf\\(M, z\\) must be non-negative"): - cc.integrated_number_density( - [0.1, 0.2], - bad_lf, - m_bright=-24.0, - m_faint=-18.0, - ) - - -def test_integrated_number_density_rejects_unbroadcastable_lf_values() -> None: - """Tests that unbroadcastable luminosity-function outputs are rejected.""" - - def bad_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: - return np.ones((3, 3), dtype=float) - - with pytest.raises( - ValueError, - match="lf\\(M, z\\) must return values broadcastable", - ): - cc.integrated_number_density( - [0.1, 0.2], - bad_lf, - m_bright=-24.0, - m_faint=-18.0, - ) - - def test_observed_number_density_integrates_to_catalog_limit( monkeypatch: pytest.MonkeyPatch, ) -> None: @@ -239,6 +118,29 @@ def test_observed_number_density_clips_catalog_limit_to_faint_bound( np.testing.assert_allclose(result, np.array([6.0, 6.0])) +def test_observed_number_density_returns_zero_when_limit_is_too_bright( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that observed density is zero when the catalog limit is too bright.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-25.0, -26.0]), + ) + + result = cc.observed_number_density( + object(), + [0.1, 0.2], + constant_lf, + m_lim=10.0, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + def test_missing_number_density_integrates_from_catalog_limit( monkeypatch: pytest.MonkeyPatch, ) -> None: @@ -285,6 +187,29 @@ def test_missing_number_density_clips_catalog_limit_to_bright_bound( np.testing.assert_allclose(result, np.array([6.0, 6.0])) +def test_missing_number_density_returns_zero_when_limit_is_too_faint( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that missing density is zero when the catalog reaches the faint bound.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-17.0, -16.0]), + ) + + result = cc.missing_number_density( + object(), + [0.1, 0.2], + constant_lf, + m_lim=30.0, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + def test_observed_number_density_rejects_invalid_magnitude_range() -> None: """Tests that observed density rejects invalid LF magnitude bounds.""" with pytest.raises(ValueError, match="m_faint must be larger than m_bright"): @@ -357,46 +282,28 @@ def test_out_of_catalog_fraction_returns_missing_fraction( np.testing.assert_allclose(result, np.array([0.5, 1.0 / 6.0])) -def test_integrated_number_density_accepts_scalar_redshift() -> None: - """Tests that finite-range integration accepts scalar redshift input.""" - result = cc.integrated_number_density( +def test_observed_number_density_accepts_scalar_redshift( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that observed density accepts scalar redshift input.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.asarray(-21.0), + ) + + result = cc.observed_number_density( + object(), 0.1, constant_lf, + m_lim=24.5, m_bright=-24.0, m_faint=-18.0, n_m=64, ) assert result.shape == () - assert result == pytest.approx(6.0) - - -def test_integrated_number_density_accepts_broadcastable_scalar_lf_output() -> None: - """Tests that scalar luminosity-function outputs are broadcast.""" - - def scalar_lf(m_abs: np.ndarray, z: np.ndarray) -> float: - return 2.0 - - result = cc.integrated_number_density( - [0.1, 0.2], - scalar_lf, - m_bright=-24.0, - m_faint=-18.0, - n_m=64, - ) - - np.testing.assert_allclose(result, np.array([12.0, 12.0])) - - -def test_integrated_number_density_rejects_nonfinite_magnitude_bounds() -> None: - """Tests that non-finite magnitude bounds are rejected.""" - with pytest.raises(ValueError, match="m_lower contains NaN or infinite values"): - cc.integrated_number_density( - [0.1, 0.2], - constant_lf, - m_bright=np.nan, - m_faint=-18.0, - ) + assert result == pytest.approx(3.0) def test_catalog_completeness_fraction_returns_zero_for_zero_total_density( @@ -456,3 +363,170 @@ def test_completeness_and_out_of_catalog_fractions_sum_to_one( ) np.testing.assert_allclose(completeness + missing, np.ones(2)) + + +def test_observed_number_density_uses_internal_lf_integral( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that observed density calls the generic LF integral internally.""" + calls = {} + + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + def fake_integrated_number_density( + z: np.ndarray, + lf: object, + *, + m_bright: float | np.ndarray, + m_faint: float | np.ndarray, + n_m: int, + ) -> np.ndarray: + calls["z"] = z + calls["lf"] = lf + calls["m_bright"] = m_bright + calls["m_faint"] = m_faint + calls["n_m"] = n_m + return np.array([10.0, 20.0]) + + monkeypatch.setattr(cc, "_integrated_number_density", fake_integrated_number_density) + + result = cc.observed_number_density( + object(), + [0.1, 0.2], + double_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=123, + ) + + np.testing.assert_allclose(result, np.array([10.0, 20.0])) + np.testing.assert_allclose(calls["z"], np.array([0.1, 0.2])) + assert calls["lf"] is double_lf + assert calls["m_bright"] == pytest.approx(-24.0) + np.testing.assert_allclose(calls["m_faint"], np.array([-21.0, -19.0])) + assert calls["n_m"] == 123 + + +def test_missing_number_density_uses_internal_lf_integral( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that missing density calls the generic LF integral internally.""" + calls = {} + + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + def fake_integrated_number_density( + z: np.ndarray, + lf: object, + *, + m_bright: float | np.ndarray, + m_faint: float | np.ndarray, + n_m: int, + ) -> np.ndarray: + calls["z"] = z + calls["lf"] = lf + calls["m_bright"] = m_bright + calls["m_faint"] = m_faint + calls["n_m"] = n_m + return np.array([30.0, 40.0]) + + monkeypatch.setattr(cc, "_integrated_number_density", fake_integrated_number_density) + + result = cc.missing_number_density( + object(), + [0.1, 0.2], + double_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=321, + ) + + np.testing.assert_allclose(result, np.array([30.0, 40.0])) + np.testing.assert_allclose(calls["z"], np.array([0.1, 0.2])) + assert calls["lf"] is double_lf + np.testing.assert_allclose(calls["m_bright"], np.array([-21.0, -19.0])) + assert calls["m_faint"] == pytest.approx(-18.0) + assert calls["n_m"] == 321 + + +def test_absolute_magnitude_limit_reads_h_from_cosmology( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that h is read from the cosmology object when not provided.""" + calls = {} + + def fake_absolute_magnitude( + cosmo_obj: object, + z: np.ndarray, + apparent_mag: float, + *, + h: float | None = None, + k_correction: float | np.ndarray | None = None, + e_correction: float | np.ndarray | None = None, + ) -> np.ndarray: + calls["h"] = h + return np.array([-20.0, -19.0]) + + monkeypatch.setattr(cc, "absolute_magnitude", fake_absolute_magnitude) + + cosmo_obj = {"h": 0.68} + + result = cc.absolute_magnitude_limit( + cosmo_obj, + [0.1, 0.2], + m_lim=24.5, + ) + + np.testing.assert_allclose(result, np.array([-20.0, -19.0])) + assert calls["h"] == pytest.approx(0.68) + + +def test_absolute_magnitude_limit_explicit_h_overrides_cosmology( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that explicit h is used instead of cosmo_obj['h'].""" + calls = {} + + def fake_absolute_magnitude( + cosmo_obj: object, + z: np.ndarray, + apparent_mag: float, + *, + h: float | None = None, + k_correction: float | np.ndarray | None = None, + e_correction: float | np.ndarray | None = None, + ) -> np.ndarray: + calls["h"] = h + return np.array([-20.0, -19.0]) + + monkeypatch.setattr(cc, "absolute_magnitude", fake_absolute_magnitude) + + cosmo_obj = {"h": 0.68} + + cc.absolute_magnitude_limit( + cosmo_obj, + [0.1, 0.2], + m_lim=24.5, + h=0.72, + ) + + assert calls["h"] == pytest.approx(0.72) + + +def test_absolute_magnitude_limit_rejects_missing_h() -> None: + """Tests that missing h raises a clear error.""" + with pytest.raises( + ValueError, + match="h was not provided and could not be read from cosmo_obj", + ): + cc.absolute_magnitude_limit(object(), [0.1, 0.2], m_lim=24.5) diff --git a/tests/test_photometry_lf_integrals.py b/tests/test_photometry_lf_integrals.py new file mode 100644 index 0000000..888f7d7 --- /dev/null +++ b/tests/test_photometry_lf_integrals.py @@ -0,0 +1,513 @@ +"""Unit tests for the ``lfkit.photometry.lf_integrals.py`` module.""" + +import numpy as np +import pytest + +import lfkit.photometry.lf_integrals as li + + +def constant_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a constant luminosity function.""" + return np.ones_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + +def double_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a constant luminosity function with amplitude two.""" + return 2.0 * np.ones_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + +def zero_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a zero luminosity function.""" + return np.zeros_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + +def linear_magnitude_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a positive LF that varies linearly with absolute magnitude.""" + return np.asarray(m_abs + 25.0, dtype=float) + + +def test_integrated_number_density_integrates_constant_lf() -> None: + """Tests that finite-range integration returns the expected width.""" + result = li.integrated_number_density( + [0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([6.0, 6.0])) + + +def test_integrated_number_density_integrates_lf_amplitude() -> None: + """Tests that finite-range integration preserves LF amplitude.""" + result = li.integrated_number_density( + [0.1, 0.2], + double_lf, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([12.0, 12.0])) + + +def test_integrated_number_density_supports_array_bounds() -> None: + """Tests that finite-range integration supports redshift-dependent bounds.""" + result = li.integrated_number_density( + [0.1, 0.2, 0.3], + constant_lf, + m_bright=np.array([-24.0, -23.0, -22.0]), + m_faint=np.array([-18.0, -18.0, -18.0]), + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([6.0, 5.0, 4.0])) + + +def test_integrated_number_density_returns_zero_for_empty_ranges() -> None: + """Tests that empty magnitude ranges return zero density.""" + result = li.integrated_number_density( + [0.1, 0.2], + constant_lf, + m_bright=np.array([-18.0, -20.0]), + m_faint=np.array([-20.0, -20.0]), + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + +def test_integrated_number_density_accepts_scalar_redshift() -> None: + """Tests that finite-range integration accepts scalar redshift input.""" + result = li.integrated_number_density( + 0.1, + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + assert result.shape == () + assert result == pytest.approx(6.0) + + +def test_integrated_number_density_accepts_broadcastable_scalar_lf_output() -> None: + """Tests that scalar luminosity-function outputs are broadcast.""" + + def scalar_lf(m_abs: np.ndarray, z: np.ndarray) -> float: + """Return a scalar LF value.""" + return 2.0 + + result = li.integrated_number_density( + [0.1, 0.2], + scalar_lf, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([12.0, 12.0])) + + +def test_integrated_number_density_rejects_negative_redshift() -> None: + """Tests that finite-range integration rejects negative redshifts.""" + with pytest.raises(ValueError, match="Redshift z must be >= 0"): + li.integrated_number_density( + [-0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + ) + + +def test_integrated_number_density_rejects_small_magnitude_grid() -> None: + """Tests that finite-range integration requires at least two grid points.""" + with pytest.raises(ValueError, match="n_m must be at least 2"): + li.integrated_number_density( + [0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + n_m=1, + ) + + +def test_integrated_number_density_rejects_nonfinite_magnitude_lower_bound() -> None: + """Tests that non-finite lower magnitude bounds are rejected.""" + with pytest.raises(ValueError, match="m_lower contains NaN or infinite values"): + li.integrated_number_density( + [0.1, 0.2], + constant_lf, + m_bright=np.nan, + m_faint=-18.0, + ) + + +def test_integrated_number_density_rejects_nonfinite_magnitude_upper_bound() -> None: + """Tests that non-finite upper magnitude bounds are rejected.""" + with pytest.raises(ValueError, match="m_upper contains NaN or infinite values"): + li.integrated_number_density( + [0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=np.inf, + ) + + +def test_integrated_number_density_rejects_nonfinite_lf_values() -> None: + """Tests that non-finite luminosity-function values are rejected.""" + + def bad_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return non-finite LF values.""" + return np.full_like(m_abs, np.nan, dtype=float) + + with pytest.raises(ValueError, match="lf\\(M, z\\) returned non-finite values"): + li.integrated_number_density( + [0.1, 0.2], + bad_lf, + m_bright=-24.0, + m_faint=-18.0, + ) + + +def test_integrated_number_density_rejects_negative_lf_values() -> None: + """Tests that negative luminosity-function values are rejected.""" + + def bad_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return negative LF values.""" + return -np.ones_like(m_abs, dtype=float) + + with pytest.raises(ValueError, match="lf\\(M, z\\) must be non-negative"): + li.integrated_number_density( + [0.1, 0.2], + bad_lf, + m_bright=-24.0, + m_faint=-18.0, + ) + + +def test_integrated_number_density_rejects_unbroadcastable_lf_values() -> None: + """Tests that unbroadcastable luminosity-function outputs are rejected.""" + + def bad_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return LF values with an invalid shape.""" + return np.ones((3, 3), dtype=float) + + with pytest.raises( + ValueError, + match="lf\\(M, z\\) must return values broadcastable", + ): + li.integrated_number_density( + [0.1, 0.2], + bad_lf, + m_bright=-24.0, + m_faint=-18.0, + ) + + +def test_lf_weighted_integral_applies_constant_weight() -> None: + """Tests that weighted LF integration applies a constant weight.""" + + def weight_fn(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return constant weights.""" + return 3.0 * np.ones_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + result = li.lf_weighted_integral( + [0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + weight_fn=weight_fn, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([18.0, 18.0])) + + +def test_lf_weighted_integral_applies_magnitude_weight() -> None: + """Tests that weighted LF integration supports magnitude-dependent weights.""" + + def weight_fn(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a positive magnitude-dependent weight.""" + return np.asarray(m_abs + 25.0, dtype=float) + + result = li.lf_weighted_integral( + 0.1, + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + weight_fn=weight_fn, + n_m=128, + ) + + assert result == pytest.approx(24.0) + + +def test_lf_weighted_integral_rejects_nonfinite_weight_values() -> None: + """Tests that non-finite weight values are rejected.""" + + def bad_weight_fn(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return non-finite weights.""" + return np.full_like(m_abs, np.nan, dtype=float) + + with pytest.raises(ValueError, match="weight_fn\\(M, z\\) returned non-finite values"): + li.lf_weighted_integral( + [0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + weight_fn=bad_weight_fn, + ) + + +def test_lf_weighted_integral_rejects_negative_weight_values() -> None: + """Tests that negative weight values are rejected.""" + + def bad_weight_fn(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return negative weights.""" + return -np.ones_like(m_abs, dtype=float) + + with pytest.raises(ValueError, match="weight_fn\\(M, z\\) must be non-negative"): + li.lf_weighted_integral( + [0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + weight_fn=bad_weight_fn, + ) + + +def test_lf_weighted_integral_rejects_unbroadcastable_weight_values() -> None: + """Tests that unbroadcastable weight outputs are rejected.""" + + def bad_weight_fn(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return weights with an invalid shape.""" + return np.ones((3, 3), dtype=float) + + with pytest.raises( + ValueError, + match="weight_fn\\(M, z\\) must return values broadcastable", + ): + li.lf_weighted_integral( + [0.1, 0.2], + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + weight_fn=bad_weight_fn, + ) + + +def test_selection_weighted_number_density_matches_weighted_integral() -> None: + """Tests that selection-weighted density applies the selection function.""" + + def selection_fn(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a constant fifty-percent selection.""" + return 0.5 * np.ones_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + result = li.selection_weighted_number_density( + [0.1, 0.2], + double_lf, + m_bright=-24.0, + m_faint=-18.0, + selection_fn=selection_fn, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([6.0, 6.0])) + + +def test_integrated_luminosity_density_uses_reference_magnitude() -> None: + """Tests luminosity-density integration with relative luminosity weights.""" + result = li.integrated_luminosity_density( + 0.1, + constant_lf, + m_bright=0.0, + m_faint=1.0, + m_reference=0.0, + n_m=1024, + ) + + expected = (1.0 - 10.0**-0.4) / (0.4 * np.log(10.0)) + assert result == pytest.approx(expected, rel=1.0e-5) + + +def test_integrated_luminosity_density_rejects_nonfinite_reference() -> None: + """Tests that luminosity-density integration rejects non-finite references.""" + with pytest.raises(ValueError, match="m_reference must be finite"): + li.integrated_luminosity_density( + 0.1, + constant_lf, + m_bright=-24.0, + m_faint=-18.0, + m_reference=np.nan, + ) + + +def test_mean_luminosity_returns_luminosity_density_over_number_density() -> None: + """Tests that mean luminosity divides luminosity density by number density.""" + result = li.mean_luminosity( + 0.1, + constant_lf, + m_bright=0.0, + m_faint=1.0, + m_reference=0.0, + n_m=1024, + ) + + expected = (1.0 - 10.0**-0.4) / (0.4 * np.log(10.0)) + assert result == pytest.approx(expected, rel=1.0e-5) + + +def test_mean_luminosity_returns_zero_for_zero_number_density() -> None: + """Tests that mean luminosity is zero when number density is zero.""" + result = li.mean_luminosity( + [0.1, 0.2], + zero_lf, + m_bright=-24.0, + m_faint=-18.0, + m_reference=0.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + +def test_cumulative_number_density_brighter_than_threshold() -> None: + """Tests cumulative number density brighter than a threshold.""" + result = li.cumulative_number_density( + [0.1, 0.2], + constant_lf, + m_threshold=-21.0, + m_bright=-24.0, + m_faint=-18.0, + brighter_than=True, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([3.0, 3.0])) + + +def test_cumulative_number_density_fainter_than_threshold() -> None: + """Tests cumulative number density fainter than a threshold.""" + result = li.cumulative_number_density( + [0.1, 0.2], + constant_lf, + m_threshold=-21.0, + m_bright=-24.0, + m_faint=-18.0, + brighter_than=False, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([3.0, 3.0])) + + +def test_cumulative_number_density_clips_bright_threshold_to_faint_bound() -> None: + """Tests that brighter-than cumulative density clips to the faint bound.""" + result = li.cumulative_number_density( + [0.1, 0.2], + constant_lf, + m_threshold=-17.0, + m_bright=-24.0, + m_faint=-18.0, + brighter_than=True, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([6.0, 6.0])) + + +def test_cumulative_number_density_returns_zero_for_too_bright_threshold() -> None: + """Tests that brighter-than cumulative density is zero for too-bright cuts.""" + result = li.cumulative_number_density( + [0.1, 0.2], + constant_lf, + m_threshold=-25.0, + m_bright=-24.0, + m_faint=-18.0, + brighter_than=True, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + +def test_cumulative_number_density_clips_faint_threshold_to_bright_bound() -> None: + """Tests that fainter-than cumulative density clips to the bright bound.""" + result = li.cumulative_number_density( + [0.1, 0.2], + constant_lf, + m_threshold=-25.0, + m_bright=-24.0, + m_faint=-18.0, + brighter_than=False, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([6.0, 6.0])) + + +def test_cumulative_number_density_returns_zero_for_too_faint_threshold() -> None: + """Tests that fainter-than cumulative density is zero for too-faint cuts.""" + result = li.cumulative_number_density( + [0.1, 0.2], + constant_lf, + m_threshold=-17.0, + m_bright=-24.0, + m_faint=-18.0, + brighter_than=False, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + +def test_cumulative_number_density_supports_array_thresholds() -> None: + """Tests that cumulative number density supports redshift-dependent thresholds.""" + result = li.cumulative_number_density( + [0.1, 0.2, 0.3], + constant_lf, + m_threshold=np.array([-22.0, -21.0, -20.0]), + m_bright=-24.0, + m_faint=-18.0, + brighter_than=True, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([2.0, 3.0, 4.0])) + + +def test_cumulative_number_density_rejects_negative_redshift() -> None: + """Tests that cumulative number density rejects negative redshifts.""" + with pytest.raises(ValueError, match="Redshift z must be >= 0"): + li.cumulative_number_density( + [-0.1, 0.2], + constant_lf, + m_threshold=-21.0, + m_bright=-24.0, + m_faint=-18.0, + ) + + +def test_private_magnitude_grid_has_expected_shape_and_edges() -> None: + """Tests that the private magnitude-grid helper builds column-wise grids.""" + result = li._magnitude_grid( + m_lower=np.array([-24.0, -23.0]), + m_upper=np.array([-18.0, -20.0]), + n_m=4, + ) + + assert result.shape == (4, 2) + np.testing.assert_allclose(result[0], np.array([-24.0, -23.0])) + np.testing.assert_allclose(result[-1], np.array([-18.0, -20.0])) + + +def test_private_safe_divide_returns_zero_for_zero_denominator() -> None: + """Tests that the private safe-division helper handles zero denominators.""" + result = li._safe_divide( + np.array([1.0, 2.0, 3.0]), + np.array([1.0, 0.0, 2.0]), + ) + + np.testing.assert_allclose(result, np.array([1.0, 0.0, 1.5])) diff --git a/tests/test_photometry_lf_redshift_density.py b/tests/test_photometry_lf_redshift_density.py new file mode 100644 index 0000000..35a55e3 --- /dev/null +++ b/tests/test_photometry_lf_redshift_density.py @@ -0,0 +1,565 @@ +"""Unit tests for the ``lfkit.photometry.lf_redshift_density.py`` module.""" + +import numpy as np +import pytest + +import lfkit.photometry.lf_redshift_density as lfrd + + +M_LIM = 26.0 +M_BRIGHT = -5.0 + + +def constant_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a constant luminosity function.""" + return np.ones_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + +def double_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a constant luminosity function with amplitude two.""" + return 2.0 * np.ones_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + +def zero_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + """Return a zero luminosity function.""" + return np.zeros_like(np.broadcast_arrays(m_abs, z)[0], dtype=float) + + +def constant_luminosity_distance(z: np.ndarray) -> np.ndarray: + """Return a constant positive luminosity distance in Mpc.""" + return np.full_like(z, 10.0, dtype=float) + + +def linear_luminosity_distance(z: np.ndarray) -> np.ndarray: + """Return a simple positive luminosity distance in Mpc.""" + return 10.0 * (1.0 + np.asarray(z, dtype=float)) + + +def constant_volume_weight(z: np.ndarray) -> np.ndarray: + """Return a constant non-negative volume weight.""" + return np.ones_like(z, dtype=float) + + +def linear_volume_weight(z: np.ndarray) -> np.ndarray: + """Return a simple redshift-dependent volume weight.""" + return 1.0 + np.asarray(z, dtype=float) + + +def expected_absolute_magnitude_limit( + m_lim: float, + luminosity_distance_mpc: np.ndarray, + *, + k_correction: float | np.ndarray = 0.0, + evolution_correction: float | np.ndarray = 0.0, +) -> np.ndarray: + """Return the expected absolute-magnitude limit for test inputs.""" + return ( + m_lim + - 5.0 * np.log10(luminosity_distance_mpc) + - 25.0 + - k_correction + + evolution_correction + ) + + +def test_lf_integrated_number_density_integrates_to_absolute_magnitude_limit() -> None: + """Tests LF integration to the apparent-magnitude-implied absolute limit.""" + result = lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + np.testing.assert_allclose(result, np.array([1.0, 1.0])) + + +def test_lf_integrated_number_density_preserves_lf_amplitude() -> None: + """Tests that LF amplitude is preserved in the redshift-density integral.""" + result = lfrd.lf_integrated_number_density( + [0.1, 0.2], + double_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + np.testing.assert_allclose(result, np.array([2.0, 2.0])) + + +def test_lf_integrated_number_density_accepts_scalar_redshift() -> None: + """Tests that LF-integrated number density accepts scalar redshift input.""" + result = lfrd.lf_integrated_number_density( + 0.1, + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + assert result.shape == () + assert result == pytest.approx(1.0) + + +def test_lf_integrated_number_density_uses_luminosity_distance_function() -> None: + """Tests that the luminosity-distance callable changes the magnitude limit.""" + z = np.array([0.0, 0.5]) + luminosity_distance = linear_luminosity_distance(z) + m_faint = expected_absolute_magnitude_limit(M_LIM, luminosity_distance) + + result = lfrd.lf_integrated_number_density( + z, + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=linear_luminosity_distance, + ) + + np.testing.assert_allclose(result, m_faint - M_BRIGHT) + + +def test_lf_integrated_number_density_applies_k_correction() -> None: + """Tests that k-corrections shift the absolute-magnitude limit.""" + result = lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + k_correction=0.5, + ) + + np.testing.assert_allclose(result, np.array([0.5, 0.5])) + + +def test_lf_integrated_number_density_applies_evolution_correction() -> None: + """Tests that evolution corrections shift the absolute-magnitude limit.""" + result = lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + evolution_correction=0.5, + ) + + np.testing.assert_allclose(result, np.array([1.5, 1.5])) + + +def test_lf_integrated_number_density_accepts_array_corrections() -> None: + """Tests that redshift-dependent corrections are applied elementwise.""" + result = lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + k_correction=np.array([0.0, 0.5]), + evolution_correction=np.array([0.0, 1.0]), + ) + + np.testing.assert_allclose(result, np.array([1.0, 1.5])) + + +def test_lf_integrated_number_density_returns_zero_when_limit_is_too_bright() -> None: + """Tests that the LF integral is zero when the faint limit is too bright.""" + result = lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=20.0, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + +def test_lf_integrated_number_density_rejects_negative_redshift() -> None: + """Tests that negative redshifts are rejected.""" + with pytest.raises(ValueError, match="Redshift z must be >= 0"): + lfrd.lf_integrated_number_density( + [-0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + +def test_lf_integrated_number_density_rejects_nonfinite_m_lim() -> None: + """Tests that non-finite apparent magnitude limits are rejected.""" + with pytest.raises(ValueError, match="m_lim must be finite"): + lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=np.inf, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + +def test_lf_integrated_number_density_rejects_nonfinite_m_bright() -> None: + """Tests that non-finite bright magnitude bounds are rejected.""" + with pytest.raises(ValueError, match="m_bright must be finite"): + lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=np.nan, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + +def test_lf_integrated_number_density_rejects_nonpositive_luminosity_distance() -> None: + """Tests that luminosity-distance callables must return positive values.""" + + def bad_luminosity_distance(z: np.ndarray) -> np.ndarray: + """Return invalid luminosity distances.""" + return np.zeros_like(z, dtype=float) + + with pytest.raises(ValueError, match="positive"): + lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=bad_luminosity_distance, + ) + + +def test_lf_integrated_number_density_rejects_nonfinite_luminosity_distance() -> None: + """Tests that luminosity-distance callables must return finite values.""" + + def bad_luminosity_distance(z: np.ndarray) -> np.ndarray: + """Return non-finite luminosity distances.""" + return np.full_like(z, np.nan, dtype=float) + + with pytest.raises(ValueError, match="finite"): + lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=bad_luminosity_distance, + ) + + +def test_lf_integrated_number_density_rejects_bad_correction_shape() -> None: + """Tests that corrections must broadcast to the redshift shape.""" + with pytest.raises( + ValueError, + match="k_correction must be scalar or broadcastable to the shape of z", + ): + lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + k_correction=np.array([0.0, 0.1, 0.2]), + ) + + +def test_lf_integrated_number_density_rejects_nonfinite_correction() -> None: + """Tests that corrections must be finite.""" + with pytest.raises(ValueError, match="evolution_correction contains NaN"): + lfrd.lf_integrated_number_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + evolution_correction=np.array([0.0, np.nan]), + ) + + +def test_lf_integrated_number_density_calls_integral_with_expected_bounds( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that LF redshift density delegates to the generic LF integral.""" + calls = {} + + def fake_integrated_number_density( + z: np.ndarray, + lf: object, + *, + m_bright: float, + m_faint: np.ndarray, + n_m: int, + ) -> np.ndarray: + calls["z"] = z + calls["lf"] = lf + calls["m_bright"] = m_bright + calls["m_faint"] = m_faint + calls["n_m"] = n_m + return np.array([10.0, 20.0]) + + monkeypatch.setattr( + lfrd, + "integrated_number_density", + fake_integrated_number_density, + ) + + result = lfrd.lf_integrated_number_density( + [0.1, 0.2], + double_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=123, + luminosity_distance_mpc_fn=constant_luminosity_distance, + ) + + np.testing.assert_allclose(result, np.array([10.0, 20.0])) + np.testing.assert_allclose(calls["z"], np.array([0.1, 0.2])) + assert calls["lf"] is double_lf + assert calls["m_bright"] == pytest.approx(M_BRIGHT) + np.testing.assert_allclose(calls["m_faint"], np.array([-4.0, -4.0])) + assert calls["n_m"] == 123 + + +def test_lf_weighted_redshift_density_returns_unnormalized_density() -> None: + """Tests unnormalized LF-weighted redshift density.""" + result = lfrd.lf_weighted_redshift_density( + [0.0, 1.0], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=linear_volume_weight, + normalize=False, + ) + + np.testing.assert_allclose(result, np.array([1.0, 2.0])) + + +def test_lf_weighted_redshift_density_normalizes_density() -> None: + """Tests normalized LF-weighted redshift density.""" + result = lfrd.lf_weighted_redshift_density( + [0.0, 1.0], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=linear_volume_weight, + normalize=True, + ) + + np.testing.assert_allclose(result, np.array([2.0 / 3.0, 4.0 / 3.0])) + assert np.trapezoid(result, x=np.array([0.0, 1.0])) == pytest.approx(1.0) + + +def test_lf_weighted_redshift_density_normalizes_on_multi_point_grid() -> None: + """Tests normalization on a multi-point redshift grid.""" + z = np.array([0.0, 0.5, 1.0]) + + result = lfrd.lf_weighted_redshift_density( + z, + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=linear_volume_weight, + normalize=True, + ) + + assert np.trapezoid(result, x=z) == pytest.approx(1.0) + + +def test_lf_weighted_redshift_density_accepts_scalar_redshift_without_normalizing() -> None: + """Tests scalar redshift input when normalization is disabled.""" + result = lfrd.lf_weighted_redshift_density( + 0.1, + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=constant_volume_weight, + normalize=False, + ) + + assert result.shape == () + assert result == pytest.approx(1.0) + + +def test_lf_weighted_redshift_density_rejects_scalar_redshift_when_normalizing() -> None: + """Tests that scalar redshift input cannot be normalized.""" + with pytest.raises(ValueError, match="at least one dimensional"): + lfrd.lf_weighted_redshift_density( + 0.1, + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=constant_volume_weight, + normalize=True, + ) + + +def test_lf_weighted_redshift_density_rejects_negative_redshift() -> None: + """Tests that weighted redshift density rejects negative redshifts.""" + with pytest.raises(ValueError, match="Redshift z must be >= 0"): + lfrd.lf_weighted_redshift_density( + [-0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=constant_volume_weight, + ) + + +def test_lf_weighted_redshift_density_rejects_negative_volume_weight() -> None: + """Tests that volume weights must be non-negative.""" + + def bad_volume_weight(z: np.ndarray) -> np.ndarray: + """Return negative volume weights.""" + return -np.ones_like(z, dtype=float) + + with pytest.raises(ValueError, match="non-negative"): + lfrd.lf_weighted_redshift_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=bad_volume_weight, + normalize=False, + ) + + +def test_lf_weighted_redshift_density_rejects_nonfinite_volume_weight() -> None: + """Tests that volume weights must be finite.""" + + def bad_volume_weight(z: np.ndarray) -> np.ndarray: + """Return non-finite volume weights.""" + return np.full_like(z, np.nan, dtype=float) + + with pytest.raises(ValueError, match="finite"): + lfrd.lf_weighted_redshift_density( + [0.1, 0.2], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=bad_volume_weight, + normalize=False, + ) + + +def test_lf_weighted_redshift_density_rejects_zero_normalization() -> None: + """Tests that normalization rejects zero-integral redshift densities.""" + with pytest.raises( + ValueError, + match="Cannot normalize LF-weighted redshift density", + ): + lfrd.lf_weighted_redshift_density( + [0.0, 1.0], + zero_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=constant_volume_weight, + normalize=True, + ) + + +def test_lf_weighted_redshift_density_forwards_corrections_to_lf_density() -> None: + """Tests that weighted redshift density forwards magnitude corrections.""" + result = lfrd.lf_weighted_redshift_density( + [0.0, 1.0], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + n_m=64, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=constant_volume_weight, + k_correction=np.array([0.0, 0.5]), + evolution_correction=np.array([0.0, 1.0]), + normalize=False, + ) + + np.testing.assert_allclose(result, np.array([1.0, 1.5])) + + +def test_lf_weighted_redshift_density_multiplies_lf_density_by_volume_weight( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that LF density is multiplied by the volume-weight callable.""" + monkeypatch.setattr( + lfrd, + "lf_integrated_number_density", + lambda *args, **kwargs: np.array([2.0, 3.0]), + ) + + result = lfrd.lf_weighted_redshift_density( + [0.0, 1.0], + constant_lf, + m_lim=M_LIM, + m_bright=M_BRIGHT, + luminosity_distance_mpc_fn=constant_luminosity_distance, + volume_weight_fn=linear_volume_weight, + normalize=False, + ) + + np.testing.assert_allclose(result, np.array([2.0, 6.0])) + + +def test_optional_correction_array_returns_zeros_for_none() -> None: + """Tests that missing corrections are converted to zeros.""" + result = lfrd._optional_correction_array( + None, + np.array([0.1, 0.2]), + name="k_correction", + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + +def test_optional_correction_array_broadcasts_scalar() -> None: + """Tests that scalar corrections are broadcast to redshift shape.""" + result = lfrd._optional_correction_array( + 0.5, + np.array([0.1, 0.2]), + name="k_correction", + ) + + np.testing.assert_allclose(result, np.array([0.5, 0.5])) + + +def test_optional_correction_array_accepts_matching_array() -> None: + """Tests that matching correction arrays are preserved.""" + result = lfrd._optional_correction_array( + np.array([0.1, 0.2]), + np.array([0.3, 0.4]), + name="k_correction", + ) + + np.testing.assert_allclose(result, np.array([0.1, 0.2])) + + +def test_optional_correction_array_rejects_unbroadcastable_array() -> None: + """Tests that correction arrays must broadcast to redshift shape.""" + with pytest.raises( + ValueError, + match="k_correction must be scalar or broadcastable to the shape of z", + ): + lfrd._optional_correction_array( + np.array([0.1, 0.2, 0.3]), + np.array([0.1, 0.2]), + name="k_correction", + ) diff --git a/tests/test_photometry_magnitudes.py b/tests/test_photometry_magnitudes.py index 9b33456..28b5ea8 100644 --- a/tests/test_photometry_magnitudes.py +++ b/tests/test_photometry_magnitudes.py @@ -5,62 +5,80 @@ import numpy as np import pytest -from lfkit.photometry import magnitudes +import lfkit.photometry.magnitudes as magnitudes +from lfkit.photometry.magnitudes import ( + absolute_magnitude_from_luminosity_distance, + apparent_magnitude_from_luminosity_distance, + total_magnitude_correction, +) -def test_total_magnitude_correction_none_inputs(): +def test_total_magnitude_correction_none_inputs() -> None: """Tests that total_magnitude_correction returns zero when no corrections are given.""" - out = magnitudes.total_magnitude_correction() + out = total_magnitude_correction() + assert isinstance(out, np.ndarray) - assert out.dtype == float - assert np.allclose(out, 0.0) + assert out.dtype == np.float64 + assert out.shape == () + assert out == pytest.approx(0.0) + +def test_total_magnitude_correction_k_only_scalar() -> None: + """Tests that total_magnitude_correction returns only the k-correction.""" + out = total_magnitude_correction(k_correction=0.3) -def test_total_magnitude_correction_k_only_scalar(): - """Tests that total_magnitude_correction returns only the k-correction for scalar input.""" - out = magnitudes.total_magnitude_correction(k_correction=0.3) assert isinstance(out, np.ndarray) - assert out.dtype == float - assert np.allclose(out, 0.3) + assert out.dtype == np.float64 + assert out.shape == () + assert out == pytest.approx(0.3) -def test_total_magnitude_correction_e_only_scalar(): - """Tests that total_magnitude_correction returns minus the e-correction for scalar input.""" - out = magnitudes.total_magnitude_correction(e_correction=0.2) +def test_total_magnitude_correction_e_only_scalar() -> None: + """Tests that total_magnitude_correction returns minus the e-correction.""" + out = total_magnitude_correction(e_correction=0.2) + assert isinstance(out, np.ndarray) - assert out.dtype == float - assert np.allclose(out, -0.2) + assert out.dtype == np.float64 + assert out.shape == () + assert out == pytest.approx(-0.2) -def test_total_magnitude_correction_adds_array_terms(): - """Tests that total_magnitude_correction combines k- and e-corrections elementwise.""" +def test_total_magnitude_correction_adds_array_terms() -> None: + """Tests that total_magnitude_correction combines K and E corrections.""" k_corr = np.array([0.1, 0.2, 0.3]) e_corr = np.array([0.05, 0.10, 0.15]) - out = magnitudes.total_magnitude_correction( + out = total_magnitude_correction( k_correction=k_corr, e_correction=e_corr, ) expected = np.array([0.05, 0.10, 0.15]) assert isinstance(out, np.ndarray) - assert out.dtype == float - assert np.allclose(out, expected) + assert out.dtype == np.float64 + np.testing.assert_allclose(out, expected) -def test_total_magnitude_correction_broadcasts_scalar_and_array(): - """Tests that total_magnitude_correction broadcasts scalar and array inputs correctly.""" - out = magnitudes.total_magnitude_correction( +def test_total_magnitude_correction_broadcasts_scalar_and_array() -> None: + """Tests that total_magnitude_correction broadcasts scalar and array inputs.""" + out = total_magnitude_correction( k_correction=0.2, e_correction=np.array([0.0, 0.1, 0.2]), ) + expected = np.array([0.2, 0.1, 0.0]) - assert np.allclose(out, expected) + np.testing.assert_allclose(out, expected) -def test_absolute_magnitude_without_corrections(monkeypatch: pytest.MonkeyPatch): - """Tests that absolute_magnitude applies M = m - mu without extra corrections.""" - def mock_distance_modulus(cosmo_obj, z, h=None): +def test_absolute_magnitude_without_corrections(monkeypatch: pytest.MonkeyPatch) -> None: + """Tests that absolute_magnitude applies M = m - mu.""" + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a fixed distance modulus.""" _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) @@ -74,13 +92,19 @@ def mock_distance_modulus(cosmo_obj, z, h=None): expected = np.array([-20.0, -19.5]) assert isinstance(out, np.ndarray) - assert out.dtype == float - assert np.allclose(out, expected) + assert out.dtype == np.float64 + np.testing.assert_allclose(out, expected) -def test_absolute_magnitude_with_corrections(monkeypatch: pytest.MonkeyPatch): +def test_absolute_magnitude_with_corrections(monkeypatch: pytest.MonkeyPatch) -> None: """Tests that absolute_magnitude applies M = m - mu - K + E.""" - def mock_distance_modulus(cosmo_obj, z, h=None): + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a fixed distance modulus.""" _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) @@ -100,12 +124,18 @@ def mock_distance_modulus(cosmo_obj, z, h=None): - np.array([0.1, 0.2]) + np.array([0.3, 0.4]) ) - assert np.allclose(out, expected) + np.testing.assert_allclose(out, expected) -def test_apparent_magnitude_without_corrections(monkeypatch: pytest.MonkeyPatch): - """Tests that apparent_magnitude applies m = M + mu without extra corrections.""" - def mock_distance_modulus(cosmo_obj, z, h=None): +def test_apparent_magnitude_without_corrections(monkeypatch: pytest.MonkeyPatch) -> None: + """Tests that apparent_magnitude applies m = M + mu.""" + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a fixed distance modulus.""" _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) @@ -119,13 +149,19 @@ def mock_distance_modulus(cosmo_obj, z, h=None): expected = np.array([20.0, 21.5]) assert isinstance(out, np.ndarray) - assert out.dtype == float - assert np.allclose(out, expected) + assert out.dtype == np.float64 + np.testing.assert_allclose(out, expected) -def test_apparent_magnitude_with_corrections(monkeypatch: pytest.MonkeyPatch): +def test_apparent_magnitude_with_corrections(monkeypatch: pytest.MonkeyPatch) -> None: """Tests that apparent_magnitude applies m = M + mu + K - E.""" - def mock_distance_modulus(cosmo_obj, z, h=None): + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a fixed distance modulus.""" _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) @@ -145,15 +181,21 @@ def mock_distance_modulus(cosmo_obj, z, h=None): + np.array([0.1, 0.2]) - np.array([0.3, 0.4]) ) - assert np.allclose(out, expected) + np.testing.assert_allclose(out, expected) def test_apparent_and_absolute_are_inverse_without_corrections( monkeypatch: pytest.MonkeyPatch, -): - """Tests that apparent_magnitude and absolute_magnitude invert each other without corrections.""" - def mock_distance_modulus(cosmo_obj, z, h=None): - _, _, _ = cosmo_obj, z, h +) -> None: + """Tests that apparent_magnitude and absolute_magnitude invert each other.""" + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a redshift-dependent distance modulus.""" + _, _ = cosmo_obj, h z_arr = np.asarray(z, dtype=float) return 40.0 + z_arr @@ -173,14 +215,20 @@ def mock_distance_modulus(cosmo_obj, z, h=None): absolute_mag=absolute_mag, ) - assert np.allclose(recovered, apparent_mag) + np.testing.assert_allclose(recovered, apparent_mag) def test_apparent_and_absolute_are_inverse_with_corrections( monkeypatch: pytest.MonkeyPatch, -): - """Tests that apparent_magnitude and absolute_magnitude invert each other with corrections.""" - def mock_distance_modulus(cosmo_obj, z, h=None): +) -> None: + """Tests that magnitude conversions invert each other with corrections.""" + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a redshift-dependent distance modulus.""" _, _ = cosmo_obj, h z_arr = np.asarray(z, dtype=float) return 39.5 + 2.0 * z_arr @@ -207,16 +255,21 @@ def mock_distance_modulus(cosmo_obj, z, h=None): e_correction=e_corr, ) - assert np.allclose(recovered, apparent_mag) + np.testing.assert_allclose(recovered, apparent_mag) def test_absolute_magnitude_passes_h_to_distance_modulus( monkeypatch: pytest.MonkeyPatch, -): +) -> None: """Tests that absolute_magnitude forwards h to distance_modulus.""" captured: dict[str, float | None] = {"h": None} - def mock_distance_modulus(cosmo_obj, z, h=None): + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Capture h and return a fixed distance modulus.""" _, _ = cosmo_obj, z captured["h"] = h return np.array([40.0], dtype=float) @@ -230,16 +283,21 @@ def mock_distance_modulus(cosmo_obj, z, h=None): h=0.7, ) - assert captured["h"] == 0.7 + assert captured["h"] == pytest.approx(0.7) def test_apparent_magnitude_passes_h_to_distance_modulus( monkeypatch: pytest.MonkeyPatch, -): +) -> None: """Tests that apparent_magnitude forwards h to distance_modulus.""" captured: dict[str, float | None] = {"h": None} - def mock_distance_modulus(cosmo_obj, z, h=None): + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Capture h and return a fixed distance modulus.""" _, _ = cosmo_obj, z captured["h"] = h return np.array([40.0], dtype=float) @@ -253,14 +311,20 @@ def mock_distance_modulus(cosmo_obj, z, h=None): h=0.7, ) - assert captured["h"] == 0.7 + assert captured["h"] == pytest.approx(0.7) def test_absolute_magnitude_broadcasts_scalar_correction( monkeypatch: pytest.MonkeyPatch, -): - """Tests that absolute_magnitude broadcasts scalar corrections over array magnitudes.""" - def mock_distance_modulus(cosmo_obj, z, h=None): +) -> None: + """Tests that absolute_magnitude broadcasts scalar corrections.""" + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a fixed distance modulus array.""" _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0, 42.0], dtype=float) @@ -274,14 +338,20 @@ def mock_distance_modulus(cosmo_obj, z, h=None): ) expected = np.array([20.0, 21.0, 22.0]) - np.array([40.0, 41.0, 42.0]) - 0.2 - assert np.allclose(out, expected) + np.testing.assert_allclose(out, expected) def test_apparent_magnitude_broadcasts_scalar_correction( monkeypatch: pytest.MonkeyPatch, -): - """Tests that apparent_magnitude broadcasts scalar corrections over array magnitudes.""" - def mock_distance_modulus(cosmo_obj, z, h=None): +) -> None: + """Tests that apparent_magnitude broadcasts scalar corrections.""" + + def mock_distance_modulus( + cosmo_obj: object, + z: np.ndarray, + h: float | None = None, + ) -> np.ndarray: + """Return a fixed distance modulus array.""" _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0, 42.0], dtype=float) @@ -295,4 +365,188 @@ def mock_distance_modulus(cosmo_obj, z, h=None): ) expected = np.array([-20.0, -20.0, -20.0]) + np.array([40.0, 41.0, 42.0]) - 0.3 - assert np.allclose(out, expected) + np.testing.assert_allclose(out, expected) + + +def test_absolute_magnitude_from_luminosity_distance_scalar() -> None: + """Tests apparent-to-absolute conversion from scalar luminosity distance.""" + result = absolute_magnitude_from_luminosity_distance( + 26.0, + 10.0, + ) + + assert result.shape == () + assert result == pytest.approx(-4.0) + + +def test_apparent_magnitude_from_luminosity_distance_scalar() -> None: + """Tests absolute-to-apparent conversion from scalar luminosity distance.""" + result = apparent_magnitude_from_luminosity_distance( + -4.0, + 10.0, + ) + + assert result.shape == () + assert result == pytest.approx(26.0) + + +def test_luminosity_distance_magnitude_conversions_are_inverse() -> None: + """Tests that luminosity-distance magnitude conversions invert each other.""" + apparent_mag = np.array([24.0, 25.0, 26.0]) + luminosity_distance = np.array([10.0, 20.0, 30.0]) + + absolute_mag = absolute_magnitude_from_luminosity_distance( + apparent_mag, + luminosity_distance, + ) + recovered = apparent_magnitude_from_luminosity_distance( + absolute_mag, + luminosity_distance, + ) + + np.testing.assert_allclose(recovered, apparent_mag) + + +def test_luminosity_distance_magnitude_conversions_are_inverse_with_corrections() -> None: + """Tests luminosity-distance conversions invert each other with corrections.""" + apparent_mag = np.array([24.0, 25.0, 26.0]) + luminosity_distance = np.array([10.0, 20.0, 30.0]) + k_corr = np.array([0.0, 0.2, 0.4]) + e_corr = np.array([0.1, 0.1, 0.2]) + + absolute_mag = absolute_magnitude_from_luminosity_distance( + apparent_mag, + luminosity_distance, + k_correction=k_corr, + e_correction=e_corr, + ) + recovered = apparent_magnitude_from_luminosity_distance( + absolute_mag, + luminosity_distance, + k_correction=k_corr, + e_correction=e_corr, + ) + + np.testing.assert_allclose(recovered, apparent_mag) + + +def test_absolute_magnitude_from_luminosity_distance_with_corrections() -> None: + """Tests apparent-to-absolute conversion with K and evolution corrections.""" + result = absolute_magnitude_from_luminosity_distance( + 26.0, + 10.0, + k_correction=0.5, + e_correction=0.2, + ) + + assert result == pytest.approx(-4.3) + + +def test_apparent_magnitude_from_luminosity_distance_with_corrections() -> None: + """Tests absolute-to-apparent conversion with K and evolution corrections.""" + result = apparent_magnitude_from_luminosity_distance( + -4.3, + 10.0, + k_correction=0.5, + e_correction=0.2, + ) + + assert result == pytest.approx(26.0) + + +def test_absolute_magnitude_from_luminosity_distance_broadcasts_arrays() -> None: + """Tests apparent-to-absolute conversion with broadcastable inputs.""" + result = absolute_magnitude_from_luminosity_distance( + 26.0, + np.array([10.0, 100.0]), + ) + + np.testing.assert_allclose(result, np.array([-4.0, -9.0])) + + +def test_apparent_magnitude_from_luminosity_distance_broadcasts_arrays() -> None: + """Tests absolute-to-apparent conversion with broadcastable inputs.""" + result = apparent_magnitude_from_luminosity_distance( + np.array([-4.0, -9.0]), + np.array([10.0, 100.0]), + ) + + np.testing.assert_allclose(result, np.array([26.0, 26.0])) + + +def test_absolute_magnitude_from_luminosity_distance_broadcasts_corrections() -> None: + """Tests apparent-to-absolute conversion with array corrections.""" + result = absolute_magnitude_from_luminosity_distance( + 26.0, + 10.0, + k_correction=np.array([0.0, 0.5]), + e_correction=np.array([0.0, 1.0]), + ) + + np.testing.assert_allclose(result, np.array([-4.0, -3.5])) + + +def test_apparent_magnitude_from_luminosity_distance_broadcasts_corrections() -> None: + """Tests absolute-to-apparent conversion with array corrections.""" + result = apparent_magnitude_from_luminosity_distance( + -4.0, + 10.0, + k_correction=np.array([0.0, 0.5]), + e_correction=np.array([0.0, 1.0]), + ) + + np.testing.assert_allclose(result, np.array([26.0, 25.5])) + + +def test_absolute_magnitude_from_luminosity_distance_rejects_zero_distance() -> None: + """Tests that apparent-to-absolute conversion rejects zero distance.""" + with pytest.raises(ValueError, match="positive values"): + absolute_magnitude_from_luminosity_distance( + 26.0, + 0.0, + ) + + +def test_apparent_magnitude_from_luminosity_distance_rejects_zero_distance() -> None: + """Tests that absolute-to-apparent conversion rejects zero distance.""" + with pytest.raises(ValueError, match="positive values"): + apparent_magnitude_from_luminosity_distance( + -4.0, + 0.0, + ) + + +def test_absolute_magnitude_from_luminosity_distance_rejects_negative_distance() -> None: + """Tests that apparent-to-absolute conversion rejects negative distance.""" + with pytest.raises(ValueError, match="negative values"): + absolute_magnitude_from_luminosity_distance( + 26.0, + -10.0, + ) + + +def test_apparent_magnitude_from_luminosity_distance_rejects_negative_distance() -> None: + """Tests that absolute-to-apparent conversion rejects negative distance.""" + with pytest.raises(ValueError, match="negative values"): + apparent_magnitude_from_luminosity_distance( + -4.0, + -10.0, + ) + + +def test_absolute_magnitude_from_luminosity_distance_rejects_nonfinite_distance() -> None: + """Tests that apparent-to-absolute conversion rejects non-finite distance.""" + with pytest.raises(ValueError, match="luminosity_distance_mpc contains NaN"): + absolute_magnitude_from_luminosity_distance( + 26.0, + np.nan, + ) + + +def test_apparent_magnitude_from_luminosity_distance_rejects_nonfinite_distance() -> None: + """Tests that absolute-to-apparent conversion rejects non-finite distance.""" + with pytest.raises(ValueError, match="luminosity_distance_mpc contains NaN"): + apparent_magnitude_from_luminosity_distance( + -4.0, + np.nan, + ) diff --git a/tests/test_utils_evaluators.py b/tests/test_utils_evaluators.py new file mode 100644 index 0000000..069f35b --- /dev/null +++ b/tests/test_utils_evaluators.py @@ -0,0 +1,397 @@ +"""Unit tests for ``lfkit.utils.evaluators.py``.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from lfkit.utils.evaluators import ( + evaluate_non_negative_redshift_callable, + evaluate_optional_redshift_callable, + evaluate_positive_redshift_callable, +) + + +def test_evaluate_optional_redshift_callable_returns_none() -> None: + """Tests that an optional callable returns None when the callable is None.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + result = evaluate_optional_redshift_callable( + None, + z, + name="test_fn", + ) + + assert result is None + + +def test_evaluate_optional_redshift_callable_evaluates_callable() -> None: + """Tests that an optional callable is evaluated when provided.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return a simple redshift-dependent array.""" + return 1.0 + z_arr + + result = evaluate_optional_redshift_callable( + fn, + z, + name="test_fn", + ) + + expected = np.array([1.0, 1.5, 2.0]) + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + np.testing.assert_allclose(result, expected) + + +def test_evaluate_optional_redshift_callable_converts_list_output() -> None: + """Tests that callable list output is converted to a float array.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> list[float]: + """Return a list with the same length as the redshift array.""" + _ = z_arr + return [1.0, 2.0, 3.0] + + result = evaluate_optional_redshift_callable( + fn, + z, + name="test_fn", + ) + + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + np.testing.assert_allclose(result, np.array([1.0, 2.0, 3.0])) + + +def test_evaluate_optional_redshift_callable_rejects_wrong_shape() -> None: + """Tests that optional callable output must match the redshift shape.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return an array with the wrong shape.""" + _ = z_arr + return np.array([1.0, 2.0], dtype=float) + + with pytest.raises(ValueError, match="test_fn\\(z\\) must return an array"): + evaluate_optional_redshift_callable( + fn, + z, + name="test_fn", + ) + + +def test_evaluate_optional_redshift_callable_rejects_scalar_output() -> None: + """Tests that scalar callable output is rejected for array redshift input.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> float: + """Return a scalar instead of an array.""" + _ = z_arr + return 1.0 + + with pytest.raises(ValueError, match="test_fn\\(z\\) must return an array"): + evaluate_optional_redshift_callable( + fn, + z, + name="test_fn", + ) + + +def test_evaluate_optional_redshift_callable_accepts_scalar_redshift_shape() -> None: + """Tests that scalar-shaped redshift input accepts scalar-shaped output.""" + z = np.asarray(0.5, dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return a scalar-shaped array.""" + return np.asarray(1.0 + z_arr, dtype=float) + + result = evaluate_optional_redshift_callable( + fn, + z, + name="test_fn", + ) + + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + assert result.shape == () + assert result == pytest.approx(1.5) + + +def test_evaluate_optional_redshift_callable_rejects_nan() -> None: + """Tests that optional callable output cannot contain NaN values.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing NaN.""" + _ = z_arr + return np.array([1.0, np.nan, 2.0], dtype=float) + + with pytest.raises(ValueError, match="test_fn\\(z\\) returned non-finite values"): + evaluate_optional_redshift_callable( + fn, + z, + name="test_fn", + ) + + +def test_evaluate_optional_redshift_callable_rejects_infinite_value() -> None: + """Tests that optional callable output cannot contain infinite values.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing infinity.""" + _ = z_arr + return np.array([1.0, np.inf, 2.0], dtype=float) + + with pytest.raises(ValueError, match="test_fn\\(z\\) returned non-finite values"): + evaluate_optional_redshift_callable( + fn, + z, + name="test_fn", + ) + + +def test_evaluate_positive_redshift_callable_accepts_positive_values() -> None: + """Tests that positive callable output is accepted.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return positive redshift-dependent values.""" + return 1.0 + z_arr + + result = evaluate_positive_redshift_callable( + fn, + z, + name="distance", + ) + + expected = np.array([1.0, 1.5, 2.0]) + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + np.testing.assert_allclose(result, expected) + + +def test_evaluate_positive_redshift_callable_rejects_zero_value() -> None: + """Tests that positive callable output cannot contain zero values.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing zero.""" + _ = z_arr + return np.array([1.0, 0.0, 2.0], dtype=float) + + with pytest.raises(ValueError, match="distance\\(z\\) must return positive values"): + evaluate_positive_redshift_callable( + fn, + z, + name="distance", + ) + + +def test_evaluate_positive_redshift_callable_rejects_negative_value() -> None: + """Tests that positive callable output cannot contain negative values.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing a negative value.""" + _ = z_arr + return np.array([1.0, -0.5, 2.0], dtype=float) + + with pytest.raises(ValueError, match="distance\\(z\\) must return positive values"): + evaluate_positive_redshift_callable( + fn, + z, + name="distance", + ) + + +def test_evaluate_positive_redshift_callable_rejects_wrong_shape() -> None: + """Tests that positive callable output must match the redshift shape.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return an array with the wrong shape.""" + _ = z_arr + return np.array([1.0, 2.0], dtype=float) + + with pytest.raises(ValueError, match="distance\\(z\\) must return an array"): + evaluate_positive_redshift_callable( + fn, + z, + name="distance", + ) + + +def test_evaluate_positive_redshift_callable_rejects_nonfinite_value() -> None: + """Tests that positive callable output must be finite.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing infinity.""" + _ = z_arr + return np.array([1.0, np.inf, 2.0], dtype=float) + + with pytest.raises(ValueError, match="distance\\(z\\) returned non-finite values"): + evaluate_positive_redshift_callable( + fn, + z, + name="distance", + ) + + +def test_evaluate_non_negative_redshift_callable_accepts_positive_values() -> None: + """Tests that non-negative callable output accepts positive values.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return positive redshift-dependent values.""" + return 1.0 + z_arr + + result = evaluate_non_negative_redshift_callable( + fn, + z, + name="weight", + ) + + expected = np.array([1.0, 1.5, 2.0]) + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + np.testing.assert_allclose(result, expected) + + +def test_evaluate_non_negative_redshift_callable_accepts_zero_value() -> None: + """Tests that non-negative callable output accepts zero values.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing zero.""" + _ = z_arr + return np.array([0.0, 1.0, 2.0], dtype=float) + + result = evaluate_non_negative_redshift_callable( + fn, + z, + name="weight", + ) + + np.testing.assert_allclose(result, np.array([0.0, 1.0, 2.0])) + + +def test_evaluate_non_negative_redshift_callable_rejects_negative_value() -> None: + """Tests that non-negative callable output cannot contain negative values.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing a negative value.""" + _ = z_arr + return np.array([1.0, -0.5, 2.0], dtype=float) + + with pytest.raises( + ValueError, + match="weight\\(z\\) must return non-negative values", + ): + evaluate_non_negative_redshift_callable( + fn, + z, + name="weight", + ) + + +def test_evaluate_non_negative_redshift_callable_rejects_wrong_shape() -> None: + """Tests that non-negative callable output must match the redshift shape.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return an array with the wrong shape.""" + _ = z_arr + return np.array([1.0, 2.0], dtype=float) + + with pytest.raises(ValueError, match="weight\\(z\\) must return an array"): + evaluate_non_negative_redshift_callable( + fn, + z, + name="weight", + ) + + +def test_evaluate_non_negative_redshift_callable_rejects_nonfinite_value() -> None: + """Tests that non-negative callable output must be finite.""" + z = np.array([0.0, 0.5, 1.0], dtype=float) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return output containing NaN.""" + _ = z_arr + return np.array([1.0, np.nan, 2.0], dtype=float) + + with pytest.raises(ValueError, match="weight\\(z\\) returned non-finite values"): + evaluate_non_negative_redshift_callable( + fn, + z, + name="weight", + ) + + +def test_evaluate_non_negative_redshift_callable_preserves_input_shape() -> None: + """Tests that callable output preserves multi-dimensional redshift shape.""" + z = np.array( + [ + [0.0, 0.5], + [1.0, 1.5], + ], + dtype=float, + ) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return values with the same two-dimensional shape as redshift.""" + return 1.0 + z_arr + + result = evaluate_non_negative_redshift_callable( + fn, + z, + name="weight", + ) + + expected = np.array( + [ + [1.0, 1.5], + [2.0, 2.5], + ], + dtype=float, + ) + assert result.shape == z.shape + assert result.dtype == np.float64 + np.testing.assert_allclose(result, expected) + + +def test_evaluate_positive_redshift_callable_preserves_input_shape() -> None: + """Tests that positive callable output preserves multi-dimensional shape.""" + z = np.array( + [ + [0.0, 0.5], + [1.0, 1.5], + ], + dtype=float, + ) + + def fn(z_arr: np.ndarray) -> np.ndarray: + """Return positive values with the same shape as redshift.""" + return 2.0 + z_arr + + result = evaluate_positive_redshift_callable( + fn, + z, + name="distance", + ) + + expected = np.array( + [ + [2.0, 2.5], + [3.0, 3.5], + ], + dtype=float, + ) + assert result.shape == z.shape + assert result.dtype == np.float64 + np.testing.assert_allclose(result, expected) diff --git a/tests/test_utils_validators.py b/tests/test_utils_validators.py index 52b7ff5..f0cf869 100644 --- a/tests/test_utils_validators.py +++ b/tests/test_utils_validators.py @@ -3,7 +3,11 @@ import numpy as np import pytest -from lfkit.utils.validators import validate_array, validate_magnitude_range +from lfkit.utils.validators import ( + validate_array, + validate_luminosity_distance, + validate_magnitude_range, +) def test_validate_array_accepts_scalar() -> None: @@ -104,4 +108,77 @@ def test_validate_magnitude_range_rejects_equal_bounds() -> None: def test_validate_magnitude_range_rejects_reversed_bounds() -> None: """Tests that reversed magnitude bounds are rejected.""" with pytest.raises(ValueError, match="m_faint must be larger than m_bright"): - validate_magnitude_range(m_bright=-18.0, m_faint=-24.0) \ No newline at end of file + validate_magnitude_range(m_bright=-18.0, m_faint=-24.0) + + +def test_validate_luminosity_distance_accepts_scalar() -> None: + """Tests that a finite positive scalar distance is accepted.""" + result = validate_luminosity_distance(100.0) + + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + assert result.shape == () + assert result == pytest.approx(100.0) + + +def test_validate_luminosity_distance_accepts_list() -> None: + """Tests that finite positive distance values are accepted.""" + result = validate_luminosity_distance([10.0, 100.0, 1000.0]) + + np.testing.assert_allclose(result, np.array([10.0, 100.0, 1000.0])) + assert result.dtype == np.float64 + + +def test_validate_luminosity_distance_accepts_numpy_array() -> None: + """Tests that a finite positive numpy array is accepted.""" + distance = np.array([10, 100, 1000]) + + result = validate_luminosity_distance(distance) + + np.testing.assert_allclose(result, np.array([10.0, 100.0, 1000.0])) + assert result.dtype == np.float64 + + +def test_validate_luminosity_distance_rejects_zero() -> None: + """Tests that zero luminosity distance is rejected.""" + with pytest.raises( + ValueError, + match="luminosity_distance_mpc must contain positive values", + ): + validate_luminosity_distance([0.0, 10.0]) + + +def test_validate_luminosity_distance_rejects_negative_values() -> None: + """Tests that negative luminosity distances are rejected.""" + with pytest.raises( + ValueError, + match="luminosity_distance_mpc contains negative values", + ): + validate_luminosity_distance([-1.0, 10.0]) + + +def test_validate_luminosity_distance_rejects_nan() -> None: + """Tests that NaN luminosity distances are rejected.""" + with pytest.raises( + ValueError, + match="luminosity_distance_mpc contains NaN or infinite values", + ): + validate_luminosity_distance([10.0, np.nan]) + + +def test_validate_luminosity_distance_rejects_positive_infinity() -> None: + """Tests that infinite luminosity distances are rejected.""" + with pytest.raises( + ValueError, + match="luminosity_distance_mpc contains NaN or infinite values", + ): + validate_luminosity_distance([10.0, np.inf]) + + +def test_validate_luminosity_distance_rejects_negative_infinity() -> None: + """Tests that negative infinite luminosity distances are rejected.""" + with pytest.raises( + ValueError, + match="luminosity_distance_mpc contains NaN or infinite values", + ): + validate_luminosity_distance([10.0, -np.inf])