From 72c54d3f107fb549f11ae7edf595da1140082f72 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 11:44:38 -0400 Subject: [PATCH 01/10] Add LF catalog completeness utilities --- src/lfkit/photometry/completeness.py | 509 +++++++++++++++++++++++++++ src/lfkit/utils/validators.py | 17 + 2 files changed, 526 insertions(+) create mode 100644 src/lfkit/photometry/completeness.py diff --git a/src/lfkit/photometry/completeness.py b/src/lfkit/photometry/completeness.py new file mode 100644 index 0000000..b5a91a6 --- /dev/null +++ b/src/lfkit/photometry/completeness.py @@ -0,0 +1,509 @@ +r"""Catalog completeness utilities for luminosity-function models. + +This module provides helpers for estimating the observed and missing galaxy +population implied by a magnitude-limited catalog. These functions are useful +for applications that need an out-of-catalog correction, such as galaxy-catalog +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. + +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 completeness machinery independent of any specific +luminosity-function parameterization. +""" + +from __future__ import annotations + +from collections.abc import Callable +from typing import TYPE_CHECKING, TypeAlias + +import numpy as np +from numpy.typing import ArrayLike, NDArray + +from lfkit.photometry.lf_parameter_models import ParameterValue +from lfkit.photometry.magnitudes import absolute_magnitude +from lfkit.utils.validators import validate_array, validate_magnitude_range + +if TYPE_CHECKING: + import pyccl as ccl + + Cosmology = ccl.Cosmology +else: + Cosmology = object + + +FloatArray: TypeAlias = NDArray[np.float64] +LuminosityFunction: TypeAlias = Callable[[FloatArray, FloatArray], FloatArray] + + +__all__ = [ + "LuminosityFunction", + "absolute_magnitude_limit", + "integrated_number_density", + "observed_number_density", + "missing_number_density", + "catalog_completeness_fraction", + "out_of_catalog_fraction", +] + + +def absolute_magnitude_limit( + cosmo_obj: Cosmology, + z: ArrayLike, + *, + m_lim: float, + h: float | None = None, + k_correction: ParameterValue | None = None, + e_correction: ParameterValue | None = None, +) -> NDArray[np.float64]: + r"""Return the absolute-magnitude limit of an apparent-magnitude catalog cut. + + This converts an apparent magnitude limit into the corresponding limiting + absolute magnitude at each redshift, + + .. math:: + + M_{\mathrm{lim}}(z) = m_{\mathrm{lim}} - \mu(z) - K(z) - E(z). + + Args: + cosmo_obj: A PyCCL cosmology object. + 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. + k_correction: Optional k-correction term(s). + e_correction: Optional evolution-correction term(s). + + Returns: + NumPy array of limiting absolute magnitudes. + """ + z_arr = validate_array(z, name="z") + + if np.any(z_arr < 0): + raise ValueError("Redshift z must be >= 0.") + + if not np.isfinite(m_lim): + raise ValueError("m_lim must be finite.") + + return absolute_magnitude( + cosmo_obj, + z_arr, + m_lim, + h=h, + k_correction=k_correction, + e_correction=e_correction, + ) + + +def integrated_number_density( + z: ArrayLike, + lf: LuminosityFunction, + *, + m_bright: ArrayLike, + m_faint: ArrayLike, + n_m: int = 512, +) -> NDArray[np.float64]: + 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: ArrayLike, + lf: LuminosityFunction, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + k_correction: ParameterValue | None = None, + e_correction: ParameterValue | None = None, +) -> NDArray[np.float64]: + r"""Return number density observable in a magnitude-limited catalog. + + This integrates the luminosity function over galaxies brighter than the + catalog limit, + + .. math:: + + n_{\mathrm{obs}}(z) = + \int_{M_{\mathrm{bright}}}^{\min[M_{\lim}(z), M_{\mathrm{faint}}]} + \phi(M, z) \, dM. + + Args: + cosmo_obj: A PyCCL cosmology object. + 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 bound of the LF model. + m_faint: Faint absolute-magnitude bound of the LF model. + n_m: Number of magnitude-grid points used for the integral. + h: Optional dimensionless Hubble parameter used in the + distance-modulus convention. + k_correction: Optional k-correction term(s). + e_correction: Optional evolution-correction term(s). + + Returns: + NumPy array of observed number densities. + """ + validate_magnitude_range(m_bright=m_bright, m_faint=m_faint) + + z_arr = validate_array(z, name="z") + m_abs_lim = absolute_magnitude_limit( + cosmo_obj, + z_arr, + m_lim=m_lim, + h=h, + k_correction=k_correction, + e_correction=e_correction, + ) + + observed_upper = np.minimum(m_abs_lim, m_faint) + + return _integrated_number_density_between_bounds( + z_arr, + lf, + m_lower=m_bright, + m_upper=observed_upper, + n_m=n_m, + ) + + +def missing_number_density( + cosmo_obj: Cosmology, + z: ArrayLike, + lf: LuminosityFunction, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + k_correction: ParameterValue | None = None, + e_correction: ParameterValue | None = None, +) -> NDArray[np.float64]: + r"""Return number density missing from a magnitude-limited catalog. + + This integrates the luminosity function over galaxies fainter than the + catalog limit, + + .. math:: + + n_{\mathrm{miss}}(z) = + \int_{\max[M_{\lim}(z), M_{\mathrm{bright}}]}^{M_{\mathrm{faint}}} + \phi(M, z) \, dM. + + Args: + cosmo_obj: A PyCCL cosmology object. + 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 bound of the LF model. + m_faint: Faint absolute-magnitude bound of the LF model. + n_m: Number of magnitude-grid points used for the integral. + h: Optional dimensionless Hubble parameter used in the + distance-modulus convention. + k_correction: Optional k-correction term(s). + e_correction: Optional evolution-correction term(s). + + Returns: + NumPy array of missing number densities. + """ + validate_magnitude_range(m_bright=m_bright, m_faint=m_faint) + + z_arr = validate_array(z, name="z") + m_abs_lim = absolute_magnitude_limit( + cosmo_obj, + z_arr, + m_lim=m_lim, + h=h, + k_correction=k_correction, + e_correction=e_correction, + ) + + missing_lower = np.maximum(m_abs_lim, m_bright) + + return _integrated_number_density_between_bounds( + z_arr, + lf, + m_lower=missing_lower, + m_upper=m_faint, + n_m=n_m, + ) + + +def catalog_completeness_fraction( + cosmo_obj: Cosmology, + z: ArrayLike, + lf: LuminosityFunction, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + k_correction: ParameterValue | None = None, + e_correction: ParameterValue | None = None, +) -> NDArray[np.float64]: + r"""Return the LF fraction observable in a magnitude-limited catalog. + + This returns + + .. math:: + + f_{\mathrm{obs}}(z) = + \frac{n_{\mathrm{obs}}(z)} + {n_{\mathrm{obs}}(z) + n_{\mathrm{miss}}(z)}. + + Args: + cosmo_obj: A PyCCL cosmology object. + 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 bound of the LF model. + m_faint: Faint absolute-magnitude bound of the LF model. + n_m: Number of magnitude-grid points used for the integral. + h: Optional dimensionless Hubble parameter used in the + distance-modulus convention. + k_correction: Optional k-correction term(s). + e_correction: Optional evolution-correction term(s). + + Returns: + NumPy array of catalog completeness fractions. + """ + observed = observed_number_density( + cosmo_obj, + z, + lf, + m_lim=m_lim, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + h=h, + k_correction=k_correction, + e_correction=e_correction, + ) + + missing = missing_number_density( + cosmo_obj, + z, + lf, + m_lim=m_lim, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + h=h, + k_correction=k_correction, + e_correction=e_correction, + ) + + return _fraction(observed, observed + missing) + + +def out_of_catalog_fraction( + cosmo_obj: Cosmology, + z: ArrayLike, + lf: LuminosityFunction, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + k_correction: ParameterValue | None = None, + e_correction: ParameterValue | None = None, +) -> NDArray[np.float64]: + r"""Return the LF fraction missing from a magnitude-limited catalog. + + This returns + + .. math:: + + f_{\mathrm{miss}}(z) = 1 - f_{\mathrm{obs}}(z). + + Args: + cosmo_obj: A PyCCL cosmology object. + 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 bound of the LF model. + m_faint: Faint absolute-magnitude bound of the LF model. + n_m: Number of magnitude-grid points used for the integral. + h: Optional dimensionless Hubble parameter used in the + distance-modulus convention. + k_correction: Optional k-correction term(s). + e_correction: Optional evolution-correction term(s). + + Returns: + NumPy array of out-of-catalog fractions. + """ + completeness = catalog_completeness_fraction( + cosmo_obj, + z, + lf, + m_lim=m_lim, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + h=h, + k_correction=k_correction, + e_correction=e_correction, + ) + + return np.asarray(1.0 - completeness, dtype=float) + + +def _integrated_number_density_between_bounds( + z: ArrayLike, + lf: LuminosityFunction, + *, + m_lower: ArrayLike, + m_upper: ArrayLike, + n_m: int, +) -> NDArray[np.float64]: + 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, +) -> NDArray[np.float64]: + 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, +) -> NDArray[np.float64]: + 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 _fraction( + numerator: FloatArray, + denominator: FloatArray, +) -> NDArray[np.float64]: + r"""Return a clipped fraction with safe zero-denominator handling.""" + with np.errstate(divide="ignore", invalid="ignore"): + result = np.divide( + numerator, + denominator, + out=np.zeros_like(numerator, dtype=float), + where=denominator > 0.0, + ) + + return np.asarray(np.clip(result, 0.0, 1.0), dtype=float) diff --git a/src/lfkit/utils/validators.py b/src/lfkit/utils/validators.py index 8233ff4..9dbbca4 100644 --- a/src/lfkit/utils/validators.py +++ b/src/lfkit/utils/validators.py @@ -11,6 +11,7 @@ __all__ = [ "validate_array", + "validate_magnitude_range", ] @@ -30,3 +31,19 @@ def validate_array( raise ValueError(f"{name} contains negative values, which are not allowed.") return np.asarray(arr, dtype=np.float64) + + +def validate_magnitude_range( + *, + m_bright: float, + m_faint: float, +) -> None: + """Validate bright and faint magnitude bounds.""" + if not np.isfinite(m_bright): + raise ValueError("m_bright must be finite.") + + if not np.isfinite(m_faint): + raise ValueError("m_faint must be finite.") + + if m_faint <= m_bright: + raise ValueError("m_faint must be larger than m_bright.") From 17b638e8ddf87fd77a530db056a296e44dc29f27 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:25:32 -0400 Subject: [PATCH 02/10] Ignore local unused directory --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index e89d5c3..8ae807d 100644 --- a/.gitignore +++ b/.gitignore @@ -18,3 +18,4 @@ notebooks/plots/*.png notebooks/plots/*.pdf notebooks/logs/ docs/_generated/ +unused/ From 4e7ec88f104f479b3ea2965712751b7f96f38de2 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:25:40 -0400 Subject: [PATCH 03/10] Add shared typing utilities --- src/lfkit/utils/types.py | 27 +++++++++ src/lfkit/utils/validators.py | 11 ++-- tests/test_utils_validators.py | 107 +++++++++++++++++++++++++++++++++ 3 files changed, 138 insertions(+), 7 deletions(-) create mode 100644 src/lfkit/utils/types.py create mode 100644 tests/test_utils_validators.py diff --git a/src/lfkit/utils/types.py b/src/lfkit/utils/types.py new file mode 100644 index 0000000..cb426d7 --- /dev/null +++ b/src/lfkit/utils/types.py @@ -0,0 +1,27 @@ +"""Shared type aliases for LFKit.""" + +from __future__ import annotations + +from collections.abc import Callable, Sequence +from typing import Any, TypeAlias + +import numpy as np +from numpy.typing import NDArray + +FloatArray: TypeAlias = NDArray[np.float64] +FloatInput: TypeAlias = float | Sequence[float] | FloatArray + +ParameterValue: TypeAlias = FloatInput +ParameterModel: TypeAlias = Callable[..., FloatArray] + +Cosmology: TypeAlias = Any +LuminosityFunction: TypeAlias = Callable[[FloatArray, FloatArray], FloatArray] + +__all__ = [ + "Cosmology", + "FloatArray", + "FloatInput", + "LuminosityFunction", + "ParameterModel", + "ParameterValue", +] diff --git a/src/lfkit/utils/validators.py b/src/lfkit/utils/validators.py index 9dbbca4..874d87f 100644 --- a/src/lfkit/utils/validators.py +++ b/src/lfkit/utils/validators.py @@ -2,12 +2,9 @@ from __future__ import annotations -from numpy.typing import NDArray import numpy as np -from typing import TypeAlias -FloatArray: TypeAlias = NDArray[np.float64] -ParameterValue: TypeAlias = float | FloatArray +from lfkit.utils.types import FloatArray, FloatInput __all__ = [ "validate_array", @@ -16,12 +13,12 @@ def validate_array( - x: ParameterValue, + x: FloatInput, *, name: str, allow_negative: bool = True, -) -> NDArray[np.float64]: - """Basic validation for numeric arrays.""" +) -> FloatArray: + """Return a finite float array.""" arr = np.asarray(x, dtype=float) if not np.all(np.isfinite(arr)): diff --git a/tests/test_utils_validators.py b/tests/test_utils_validators.py new file mode 100644 index 0000000..52b7ff5 --- /dev/null +++ b/tests/test_utils_validators.py @@ -0,0 +1,107 @@ +"""Unit tests for the ``lfkit.utils.validators.py``""" + +import numpy as np +import pytest + +from lfkit.utils.validators import validate_array, validate_magnitude_range + + +def test_validate_array_accepts_scalar() -> None: + """Tests that a finite scalar is converted to a float array.""" + result = validate_array(1.5, name="x") + + assert isinstance(result, np.ndarray) + assert result.dtype == np.float64 + assert result.shape == () + assert result == pytest.approx(1.5) + + +def test_validate_array_accepts_list() -> None: + """Tests that a finite list is converted to a float array.""" + result = validate_array([1.0, 2.0, 3.0], name="x") + + np.testing.assert_allclose(result, np.array([1.0, 2.0, 3.0])) + assert result.dtype == np.float64 + + +def test_validate_array_accepts_numpy_array() -> None: + """Tests that a finite numpy array is preserved as a float array.""" + x = np.array([1, 2, 3]) + + result = validate_array(x, name="x") + + np.testing.assert_allclose(result, np.array([1.0, 2.0, 3.0])) + assert result.dtype == np.float64 + + +def test_validate_array_accepts_negative_values_by_default() -> None: + """Tests that negative values are allowed by default.""" + result = validate_array([-1.0, 0.0, 1.0], name="x") + + np.testing.assert_allclose(result, np.array([-1.0, 0.0, 1.0])) + + +def test_validate_array_rejects_negative_values_when_disallowed() -> None: + """Tests that negative values are rejected when requested.""" + with pytest.raises(ValueError, match="x contains negative values"): + validate_array([-1.0, 0.0, 1.0], name="x", allow_negative=False) + + +def test_validate_array_accepts_zero_when_negative_values_are_disallowed() -> None: + """Tests that zero is allowed when only negative values are disallowed.""" + result = validate_array([0.0, 1.0, 2.0], name="x", allow_negative=False) + + np.testing.assert_allclose(result, np.array([0.0, 1.0, 2.0])) + + +def test_validate_array_rejects_nan() -> None: + """Tests that NaN values are rejected.""" + with pytest.raises(ValueError, match="x contains NaN or infinite values"): + validate_array([1.0, np.nan], name="x") + + +def test_validate_array_rejects_positive_infinity() -> None: + """Tests that positive infinite values are rejected.""" + with pytest.raises(ValueError, match="x contains NaN or infinite values"): + validate_array([1.0, np.inf], name="x") + + +def test_validate_array_rejects_negative_infinity() -> None: + """Tests that negative infinite values are rejected.""" + with pytest.raises(ValueError, match="x contains NaN or infinite values"): + validate_array([1.0, -np.inf], name="x") + + +def test_validate_array_uses_name_in_error_message() -> None: + """Tests that the provided parameter name appears in errors.""" + with pytest.raises(ValueError, match="redshift contains NaN or infinite values"): + validate_array([1.0, np.nan], name="redshift") + + +def test_validate_magnitude_range_accepts_valid_bounds() -> None: + """Tests that valid magnitude bounds are accepted.""" + validate_magnitude_range(m_bright=-24.0, m_faint=-18.0) + + +def test_validate_magnitude_range_rejects_nonfinite_bright_bound() -> None: + """Tests that non-finite bright magnitude bounds are rejected.""" + with pytest.raises(ValueError, match="m_bright must be finite"): + validate_magnitude_range(m_bright=np.nan, m_faint=-18.0) + + +def test_validate_magnitude_range_rejects_nonfinite_faint_bound() -> None: + """Tests that non-finite faint magnitude bounds are rejected.""" + with pytest.raises(ValueError, match="m_faint must be finite"): + validate_magnitude_range(m_bright=-24.0, m_faint=np.inf) + + +def test_validate_magnitude_range_rejects_equal_bounds() -> None: + """Tests that equal magnitude bounds are rejected.""" + with pytest.raises(ValueError, match="m_faint must be larger than m_bright"): + validate_magnitude_range(m_bright=-20.0, m_faint=-20.0) + + +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 From 8dcade7e94b76789beeb90e331b8b9e0ca13f003 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:25:47 -0400 Subject: [PATCH 04/10] Add catalog completeness photometry utilities --- ...ompleteness.py => catalog_completeness.py} | 87 ++-- src/lfkit/photometry/lf_parameter_models.py | 117 ++++- src/lfkit/photometry/luminosities.py | 44 +- src/lfkit/photometry/luminosity_function.py | 73 ++- src/lfkit/photometry/magnitudes.py | 95 ++-- tests/test_photometry_completeness.py | 458 ++++++++++++++++++ ...st_photometry_completeness_fake_catalog.py | 251 ++++++++++ tests/test_photometry_lf_parameter_models.py | 194 +++++++- tests/test_photometry_magnitudes.py | 40 +- 9 files changed, 1155 insertions(+), 204 deletions(-) rename src/lfkit/photometry/{completeness.py => catalog_completeness.py} (89%) create mode 100644 tests/test_photometry_completeness.py create mode 100644 tests/test_photometry_completeness_fake_catalog.py diff --git a/src/lfkit/photometry/completeness.py b/src/lfkit/photometry/catalog_completeness.py similarity index 89% rename from src/lfkit/photometry/completeness.py rename to src/lfkit/photometry/catalog_completeness.py index b5a91a6..540fccf 100644 --- a/src/lfkit/photometry/completeness.py +++ b/src/lfkit/photometry/catalog_completeness.py @@ -20,30 +20,19 @@ from __future__ import annotations -from collections.abc import Callable -from typing import TYPE_CHECKING, TypeAlias - import numpy as np -from numpy.typing import ArrayLike, NDArray -from lfkit.photometry.lf_parameter_models import ParameterValue from lfkit.photometry.magnitudes import absolute_magnitude +from lfkit.utils.types import ( + Cosmology, + FloatArray, + FloatInput, + LuminosityFunction, +) from lfkit.utils.validators import validate_array, validate_magnitude_range -if TYPE_CHECKING: - import pyccl as ccl - - Cosmology = ccl.Cosmology -else: - Cosmology = object - - -FloatArray: TypeAlias = NDArray[np.float64] -LuminosityFunction: TypeAlias = Callable[[FloatArray, FloatArray], FloatArray] - __all__ = [ - "LuminosityFunction", "absolute_magnitude_limit", "integrated_number_density", "observed_number_density", @@ -55,13 +44,13 @@ def absolute_magnitude_limit( cosmo_obj: Cosmology, - z: ArrayLike, + z: FloatInput, *, m_lim: float, h: float | None = None, - k_correction: ParameterValue | None = None, - e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: r"""Return the absolute-magnitude limit of an apparent-magnitude catalog cut. This converts an apparent magnitude limit into the corresponding limiting @@ -69,7 +58,7 @@ def absolute_magnitude_limit( .. math:: - M_{\mathrm{lim}}(z) = m_{\mathrm{lim}} - \mu(z) - K(z) - E(z). + M_{\mathrm{lim}}(z) = m_{\mathrm{lim}} - \mu(z) - K(z) + E(z). Args: cosmo_obj: A PyCCL cosmology object. @@ -102,13 +91,13 @@ def absolute_magnitude_limit( def integrated_number_density( - z: ArrayLike, + z: FloatInput, lf: LuminosityFunction, *, - m_bright: ArrayLike, - m_faint: ArrayLike, + m_bright: FloatInput, + m_faint: FloatInput, n_m: int = 512, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return finite-range number density from a luminosity function. This computes @@ -141,7 +130,7 @@ def integrated_number_density( def observed_number_density( cosmo_obj: Cosmology, - z: ArrayLike, + z: FloatInput, lf: LuminosityFunction, *, m_lim: float, @@ -149,9 +138,9 @@ def observed_number_density( m_faint: float, n_m: int = 512, h: float | None = None, - k_correction: ParameterValue | None = None, - e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: r"""Return number density observable in a magnitude-limited catalog. This integrates the luminosity function over galaxies brighter than the @@ -204,7 +193,7 @@ def observed_number_density( def missing_number_density( cosmo_obj: Cosmology, - z: ArrayLike, + z: FloatInput, lf: LuminosityFunction, *, m_lim: float, @@ -212,9 +201,9 @@ def missing_number_density( m_faint: float, n_m: int = 512, h: float | None = None, - k_correction: ParameterValue | None = None, - e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: r"""Return number density missing from a magnitude-limited catalog. This integrates the luminosity function over galaxies fainter than the @@ -267,7 +256,7 @@ def missing_number_density( def catalog_completeness_fraction( cosmo_obj: Cosmology, - z: ArrayLike, + z: FloatInput, lf: LuminosityFunction, *, m_lim: float, @@ -275,9 +264,9 @@ def catalog_completeness_fraction( m_faint: float, n_m: int = 512, h: float | None = None, - k_correction: ParameterValue | None = None, - e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: r"""Return the LF fraction observable in a magnitude-limited catalog. This returns @@ -335,7 +324,7 @@ def catalog_completeness_fraction( def out_of_catalog_fraction( cosmo_obj: Cosmology, - z: ArrayLike, + z: FloatInput, lf: LuminosityFunction, *, m_lim: float, @@ -343,9 +332,9 @@ def out_of_catalog_fraction( m_faint: float, n_m: int = 512, h: float | None = None, - k_correction: ParameterValue | None = None, - e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: r"""Return the LF fraction missing from a magnitude-limited catalog. This returns @@ -387,13 +376,13 @@ def out_of_catalog_fraction( def _integrated_number_density_between_bounds( - z: ArrayLike, + z: FloatInput, lf: LuminosityFunction, *, - m_lower: ArrayLike, - m_upper: ArrayLike, + m_lower: FloatInput, + m_upper: FloatInput, n_m: int, -) -> NDArray[np.float64]: +) -> 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") @@ -436,7 +425,7 @@ def _magnitude_grid( m_lower: FloatArray, m_upper: FloatArray, n_m: int, -) -> NDArray[np.float64]: +) -> 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( @@ -450,7 +439,7 @@ def _evaluate_lf_on_grid( *, m_grid: FloatArray, z_grid: FloatArray, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return LF values evaluated on a magnitude-redshift grid.""" phi = np.asarray(lf(m_grid, z_grid), dtype=float) @@ -496,7 +485,7 @@ def _validate_integration_inputs( def _fraction( numerator: FloatArray, denominator: FloatArray, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a clipped fraction with safe zero-denominator handling.""" with np.errstate(divide="ignore", invalid="ignore"): result = np.divide( diff --git a/src/lfkit/photometry/lf_parameter_models.py b/src/lfkit/photometry/lf_parameter_models.py index 2c25a16..b947423 100644 --- a/src/lfkit/photometry/lf_parameter_models.py +++ b/src/lfkit/photometry/lf_parameter_models.py @@ -14,30 +14,14 @@ from __future__ import annotations from collections.abc import Mapping -from typing import TYPE_CHECKING, Callable import numpy as np -from numpy.typing import NDArray -from typing import TypeAlias +from lfkit.utils.types import FloatArray, ParameterModel, ParameterValue from lfkit.utils.validators import validate_array -if TYPE_CHECKING: - import pyccl as ccl - - Cosmology = ccl.Cosmology -else: - Cosmology = object - -ParameterModel = Callable[..., NDArray[np.float64]] -FloatArray: TypeAlias = NDArray[np.float64] -ParameterValue: TypeAlias = float | FloatArray - __all__ = [ - "ParameterModel", - "FloatArray", - "ParameterValue", "phi_star_constant", "phi_star_linear_p", "m_star_constant", @@ -47,6 +31,10 @@ "PHI_STAR_MODELS", "M_STAR_MODELS", "ALPHA_MODELS", + "available_lf_parameter_models", + "register_phi_star_model", + "register_m_star_model", + "register_alpha_model", "get_parameter_model", "evaluate_lf_parameters", ] @@ -56,7 +44,7 @@ def phi_star_constant( z: FloatArray, *, phi_star: float, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a constant Schechter normalization over redshift. This uses @@ -82,7 +70,7 @@ def phi_star_linear_p( *, phi_0_star: float, p: float, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a Schechter normalization with density evolution. This uses the common density-evolution form @@ -105,7 +93,7 @@ def m_star_constant( z: FloatArray, *, m_star: float, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a constant characteristic magnitude over redshift This uses @@ -132,7 +120,7 @@ def m_star_linear_q( m_0_star: float, q: float, z_ref: float = 0.1, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a characteristic magnitude with luminosity evolution. This uses the common luminosity-evolution form @@ -158,7 +146,7 @@ def alpha_constant( z: FloatArray, *, alpha: float, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a constant faint-end slope over redshift. This uses @@ -185,7 +173,7 @@ def alpha_linear( alpha_0: float, alpha_1: float, z_ref: float = 0.1, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a faint-end slope that varies linearly with redshift. This uses @@ -236,6 +224,87 @@ def get_parameter_model( ) from exc +def available_lf_parameter_models() -> dict[str, list[str]]: + """Return available luminosity-function parameter evolution models.""" + return { + "phi_star": sorted(PHI_STAR_MODELS), + "m_star": sorted(M_STAR_MODELS), + "alpha": sorted(ALPHA_MODELS), + } + + +def _register_parameter_model( + name: str, + model: ParameterModel, + registry: dict[str, ParameterModel], + *, + model_kind: str, + overwrite: bool = False, +) -> None: + """Register a luminosity-function parameter evolution model.""" + if not name: + raise ValueError(f"{model_kind} model name cannot be empty.") + + if not callable(model): + raise TypeError(f"{model_kind} model must be callable.") + + if name in registry and not overwrite: + raise ValueError( + f"{model_kind} model '{name}' is already registered. " + "Use overwrite=True to replace it." + ) + + registry[name] = model + + +def register_phi_star_model( + name: str, + model: ParameterModel, + *, + overwrite: bool = False, +) -> None: + """Register a phi_star evolution model.""" + _register_parameter_model( + name, + model, + PHI_STAR_MODELS, + model_kind="phi_star", + overwrite=overwrite, + ) + + +def register_m_star_model( + name: str, + model: ParameterModel, + *, + overwrite: bool = False, +) -> None: + """Register an M_star evolution model.""" + _register_parameter_model( + name, + model, + M_STAR_MODELS, + model_kind="m_star", + overwrite=overwrite, + ) + + +def register_alpha_model( + name: str, + model: ParameterModel, + *, + overwrite: bool = False, +) -> None: + """Register an alpha evolution model.""" + _register_parameter_model( + name, + model, + ALPHA_MODELS, + model_kind="alpha", + overwrite=overwrite, + ) + + def evaluate_lf_parameters( z: FloatArray, *, @@ -245,7 +314,7 @@ def evaluate_lf_parameters( m_star_kwargs: Mapping[str, ParameterValue] | None = None, alpha_model: str = "constant", alpha_kwargs: Mapping[str, ParameterValue] | None = None, -) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: +) -> tuple[FloatArray, FloatArray, FloatArray]: r"""Evaluate evolving luminosity-function parameters at redshift ``z``. Args: diff --git a/src/lfkit/photometry/luminosities.py b/src/lfkit/photometry/luminosities.py index bd38e48..30dd495 100644 --- a/src/lfkit/photometry/luminosities.py +++ b/src/lfkit/photometry/luminosities.py @@ -20,9 +20,9 @@ import numpy as np from numpy.random import Generator, default_rng -from numpy.typing import ArrayLike, NDArray from scipy.special import gamma, gammaincc +from lfkit.utils.types import FloatArray, FloatInput from lfkit.utils.validators import validate_array __all__ = ( @@ -40,9 +40,9 @@ def luminosity_ratio( - absolute_mag: ArrayLike, - m_star: ArrayLike, -) -> NDArray[np.float64]: + absolute_mag: FloatInput, + m_star: FloatInput, +) -> FloatArray: r"""Return the luminosity ratio relative to the characteristic luminosity. For magnitudes, @@ -69,9 +69,9 @@ def luminosity_ratio( def luminosity_ratio_from_magnitudes( - magnitude: ArrayLike, - ref_magnitude: ArrayLike, -) -> NDArray[np.float64]: + magnitude: FloatInput, + ref_magnitude: FloatInput, +) -> FloatArray: r"""Return luminosity relative to a reference magnitude. This uses :math:`L / L_{\mathrm{ref}} = 10^{-0.4 (m - m_{\mathrm{ref}})}`. @@ -93,8 +93,8 @@ def luminosity_ratio_from_magnitudes( def magnitude_difference_from_luminosity_ratio( - luminosity_ratio: ArrayLike, -) -> NDArray[np.float64]: + ratio: FloatInput, +) -> FloatArray: r"""Return the magnitude difference corresponding to a luminosity ratio. This uses @@ -104,7 +104,7 @@ def magnitude_difference_from_luminosity_ratio( m_1 - m_2 = -2.5 \log_{10}(L_1 / L_2). Args: - luminosity_ratio: Luminosity ratio value(s) ``L1 / L2``. + ratio: Luminosity ratio value(s) ``L1 / L2``. Returns: NumPy array of magnitude differences ``m1 - m2``. @@ -112,18 +112,18 @@ def magnitude_difference_from_luminosity_ratio( Raises: ValueError: If any luminosity ratio is not strictly positive. """ - ratio = validate_array(luminosity_ratio, name="luminosity_ratio") + ratio_arr = validate_array(ratio, name="luminosity_ratio") - if np.any(ratio <= 0): + if np.any(ratio_arr <= 0): raise ValueError("luminosity_ratio must be strictly positive.") - return np.asarray(-2.5 * np.log10(ratio), dtype=float) + return np.asarray(-2.5 * np.log10(ratio_arr), dtype=float) def luminosity_weight_from_magnitude( - magnitude: ArrayLike, + magnitude: FloatInput, reference_magnitude: float = 0.0, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return an unnormalized luminosity weight from a magnitude. This uses @@ -147,11 +147,11 @@ def luminosity_weight_from_magnitude( def luminosity_from_magnitude( - magnitude: ArrayLike, + magnitude: FloatInput, *, reference_magnitude: float = 0.0, reference_luminosity: float = 1.0, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return luminosity corresponding to a magnitude relative to a reference. This uses @@ -183,12 +183,12 @@ def luminosity_from_magnitude( def schechter_cumulative_number_density_luminosity( - luminosity_min: ArrayLike, + luminosity_min: FloatInput, *, phi_star: float, l_star: float, alpha: float, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return cumulative Schechter number density above a luminosity threshold. This evaluates @@ -327,7 +327,7 @@ def sample_schechter_luminosity( l_star: float, alpha: float, rng: Generator | None = None, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Sample luminosities from a normalized Schechter distribution. This samples from @@ -375,11 +375,11 @@ def sample_schechter_luminosity( def schechter_selection_function( - luminosity_min: ArrayLike, + luminosity_min: FloatInput, *, l_star: float, alpha: float, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return the Schechter selection fraction above a luminosity threshold. This evaluates diff --git a/src/lfkit/photometry/luminosity_function.py b/src/lfkit/photometry/luminosity_function.py index d571caf..c4b5f78 100644 --- a/src/lfkit/photometry/luminosity_function.py +++ b/src/lfkit/photometry/luminosity_function.py @@ -1,4 +1,4 @@ -r"""Luminosity function utilities for LFKit. +"""Luminosity function utilities for LFKit. This module provides simple standalone functions for evaluating common galaxy luminosity-function parameterization. @@ -49,31 +49,16 @@ from __future__ import annotations from collections.abc import Mapping -from typing import TYPE_CHECKING, Callable import warnings import numpy as np -from numpy.typing import NDArray -from typing import TypeAlias from scipy.special import gammaincc, gamma -from lfkit.photometry.magnitudes import absolute_magnitude +from lfkit.photometry.lf_parameter_models import evaluate_lf_parameters from lfkit.photometry.luminosities import luminosity_ratio +from lfkit.photometry.magnitudes import absolute_magnitude +from lfkit.utils.types import Cosmology, FloatArray, FloatInput, ParameterValue from lfkit.utils.validators import validate_array -from lfkit.photometry.lf_parameter_models import ( - ParameterValue, - evaluate_lf_parameters, -) - -if TYPE_CHECKING: - import pyccl as ccl - - Cosmology = ccl.Cosmology -else: - Cosmology = object - -ParameterModel = Callable[..., NDArray[np.float64]] -FloatArray: TypeAlias = NDArray[np.float64] __all__ = [ @@ -89,12 +74,12 @@ def schechter( - absolute_mag: FloatArray, + absolute_mag: FloatInput, *, phi_star: ParameterValue, m_star: ParameterValue, alpha: ParameterValue, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return the standard Schechter luminosity function in magnitude space. This computes @@ -140,8 +125,8 @@ def schechter( def schechter_evolving( - absolute_mag: FloatArray, - z: FloatArray, + absolute_mag: FloatInput, + z: FloatInput, *, phi_model: str = "linear_p", phi_kwargs: Mapping[str, ParameterValue] | None = None, @@ -149,7 +134,7 @@ def schechter_evolving( m_star_kwargs: Mapping[str, ParameterValue] | None = None, alpha_model: str = "constant", alpha_kwargs: Mapping[str, ParameterValue] | None = None, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return an evolving Schechter luminosity function in magnitude space. This evaluates @@ -203,14 +188,14 @@ def schechter_evolving( def schechter_double( - absolute_mag: FloatArray, + absolute_mag: FloatInput, *, phi_star: ParameterValue, m_star: ParameterValue, alpha: float, beta: float, m_transition: ParameterValue, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a double-power-law Schechter luminosity function in magnitude space. This implements the Loveday/GAMA-style faint-end extension @@ -292,23 +277,23 @@ def schechter_double( def schechter_from_m( cosmo_obj: Cosmology, - z: FloatArray, - apparent_mag: FloatArray, + z: FloatInput, + apparent_mag: FloatInput, *, phi_star: ParameterValue, m_star: ParameterValue, - alpha: float, + alpha: ParameterValue, h: float | None = None, k_correction: ParameterValue | None = None, e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return the standard Schechter luminosity function from apparent magnitudes. This uses .. math:: - M = m - \mu - K - E, + M = m - \mu - K + E, followed by @@ -368,8 +353,8 @@ def schechter_from_m( def schechter_evolving_from_m( cosmo_obj: Cosmology, - z: FloatArray, - apparent_mag: FloatArray, + z: FloatInput, + apparent_mag: FloatInput, *, phi_model: str = "linear_p", phi_kwargs: Mapping[str, ParameterValue] | None = None, @@ -380,14 +365,14 @@ def schechter_evolving_from_m( h: float | None = None, k_correction: ParameterValue | None = None, e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return an evolving Schechter luminosity function from apparent magnitudes. This uses .. math:: - M = m - \mu - K - E, + M = m - \mu - K + E, followed by @@ -468,8 +453,8 @@ def schechter_evolving_from_m( def schechter_double_from_m( cosmo_obj: Cosmology, - z: FloatArray, - apparent_mag: FloatArray, + z: FloatInput, + apparent_mag: FloatInput, *, phi_star: ParameterValue, m_star: ParameterValue, @@ -479,14 +464,14 @@ def schechter_double_from_m( h: float | None = None, k_correction: ParameterValue | None = None, e_correction: ParameterValue | None = None, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return a double-power-law Schechter luminosity function from apparent magnitudes. This uses .. math:: - M = m - \mu - K - E, + M = m - \mu - K + E, followed by @@ -550,13 +535,13 @@ def schechter_double_from_m( def schechter_cumulative( - magnitude_limit: FloatArray, + magnitude_limit: FloatInput, *, phi_star: float, m_star: float, alpha: float, brighter_than: bool = True, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return the cumulative number density for a standard Schechter LF. For a magnitude threshold :math:`M_{\mathrm{lim}}`, define @@ -627,8 +612,8 @@ def schechter_cumulative( def schechter_cumulative_evolving( - magnitude_limit: FloatArray, - z: FloatArray, + magnitude_limit: FloatInput, + z: FloatInput, *, phi_model: str = "linear_p", phi_kwargs: Mapping[str, ParameterValue] | None = None, @@ -637,7 +622,7 @@ def schechter_cumulative_evolving( alpha_model: str = "constant", alpha_kwargs: Mapping[str, ParameterValue] | None = None, brighter_than: bool = True, -) -> NDArray[np.float64]: +) -> FloatArray: r"""Return the cumulative number density for an evolving Schechter LF. For a magnitude threshold :math:`M_{\mathrm{lim}}`, define diff --git a/src/lfkit/photometry/magnitudes.py b/src/lfkit/photometry/magnitudes.py index 7bc2aa7..e24cba6 100644 --- a/src/lfkit/photometry/magnitudes.py +++ b/src/lfkit/photometry/magnitudes.py @@ -1,12 +1,12 @@ """Magnitude conversion utilities for LFKit. -This module provides helpers to translate between apparent and -absolute magnitudes using the cosmology utilities defined in LFKit. +This module provides helpers to translate between apparent ``m`` and +absolute magnitudes ``M`` using the cosmology utilities defined in LFKit. The conversions follow the convention - M = m - mu - K - E - m = M + mu + K + E + M = m - mu - K + E + m = M + mu + K - E where ``mu`` is the distance modulus, ``K`` is the k-correction, and ``E`` is the evolution correction. @@ -16,20 +16,10 @@ from __future__ import annotations -from typing import TYPE_CHECKING - import numpy as np -from numpy.typing import ArrayLike, NDArray from lfkit.cosmo.cosmology import distance_modulus - -if TYPE_CHECKING: - import pyccl as ccl - - Cosmology = ccl.Cosmology -else: - Cosmology = object - +from lfkit.utils.types import Cosmology, FloatArray, FloatInput __all__ = ( "total_magnitude_correction", @@ -40,59 +30,58 @@ def total_magnitude_correction( *, - k_correction: ArrayLike | float | None = None, - e_correction: ArrayLike | float | None = None, -) -> NDArray[np.float64]: - """Return the total additive magnitude correction. + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: + """Return the net correction added to apparent absolute conversion. + + This combines optional k-correction and evolution-correction terms as - This combines optional k-correction and evolution-correction - terms into a single array: + ``correction = K - E`` - ``correction = K + E`` + so that + + ``M = m - mu - correction``. Args: k_correction: Optional k-correction term(s). e_correction: Optional evolution-correction term(s). Returns: - NumPy array of total correction values. + NumPy array of net correction values. """ - correction = 0.0 + correction = np.asarray(0.0, dtype=float) if k_correction is not None: correction = correction + np.asarray(k_correction, dtype=float) if e_correction is not None: - correction = correction + np.asarray(e_correction, dtype=float) + correction = correction - np.asarray(e_correction, dtype=float) - return np.asarray(correction, dtype=float) + return np.asarray(correction, dtype=np.float64) def absolute_magnitude( cosmo_obj: Cosmology, - z: ArrayLike, - apparent_mag: ArrayLike, + z: FloatInput, + apparent_mag: FloatInput, *, h: float | None = None, - k_correction: ArrayLike | float | None = None, - e_correction: ArrayLike | float | None = None, -) -> NDArray[np.float64]: - """Convert apparent magnitude to absolute magnitude. + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: + """Convert apparent magnitude m to absolute magnitude M. The convention used is - ``M = m - mu - K - E`` - - where ``mu`` is the distance modulus, ``K`` is the k-correction, - and ``E`` is the evolution correction. + ``M = m - mu - K + E``. Args: - cosmo_obj: A PyCCL cosmology object. + cosmo_obj: Cosmology object passed to LFKit distance utilities. z: Redshift value or array-like of redshift values. apparent_mag: Apparent magnitude value(s). h: Optional dimensionless Hubble parameter. If provided, - the distance modulus uses the convention - ``mu = 5 log10(d_L * h) + 25``. + the distance modulus uses ``mu = 5 log10(d_L * h) + 25``. k_correction: Optional k-correction term(s). e_correction: Optional evolution-correction term(s). @@ -105,34 +94,31 @@ def absolute_magnitude( k_correction=k_correction, e_correction=e_correction, ) - return np.asarray(m - mu - correction, dtype=float) + + return np.asarray(m - mu - correction, dtype=np.float64) def apparent_magnitude( cosmo_obj: Cosmology, - z: ArrayLike, - absolute_mag: ArrayLike, + z: FloatInput, + absolute_mag: FloatInput, *, h: float | None = None, - k_correction: ArrayLike | float | None = None, - e_correction: ArrayLike | float | None = None, -) -> NDArray[np.float64]: - """Convert absolute magnitude to apparent magnitude. + k_correction: FloatInput | None = None, + e_correction: FloatInput | None = None, +) -> FloatArray: + """Convert absolute magnitude M to apparent magnitude m. The convention used is - ``m = M + mu + K + E`` - - where ``mu`` is the distance modulus, ``K`` is the k-correction, - and ``E`` is the evolution correction. + ``m = M + mu + K - E``. Args: - cosmo_obj: A PyCCL cosmology object. + cosmo_obj: Cosmology object passed to LFKit distance utilities. z: Redshift value or array-like of redshift values. absolute_mag: Absolute magnitude value(s). h: Optional dimensionless Hubble parameter. If provided, - the distance modulus uses the convention - ``mu = 5 log10(d_L * h) + 25``. + the distance modulus uses ``mu = 5 log10(d_L * h) + 25``. k_correction: Optional k-correction term(s). e_correction: Optional evolution-correction term(s). @@ -145,4 +131,5 @@ def apparent_magnitude( k_correction=k_correction, e_correction=e_correction, ) - return np.asarray(M + mu + correction, dtype=float) + + return np.asarray(M + mu + correction, dtype=np.float64) diff --git a/tests/test_photometry_completeness.py b/tests/test_photometry_completeness.py new file mode 100644 index 0000000..944b70d --- /dev/null +++ b/tests/test_photometry_completeness.py @@ -0,0 +1,458 @@ +"""Unit tests for the ``lfkit.photometry.catalog_completeness.py``""" + +import numpy as np +import pytest + +import lfkit.photometry.catalog_completeness as cc + + +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 test_absolute_magnitude_limit_calls_absolute_magnitude(monkeypatch: pytest.MonkeyPatch) -> None: + """Tests that absolute magnitude limits are delegated to magnitude conversion.""" + 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["cosmo_obj"] = cosmo_obj + calls["z"] = z + calls["apparent_mag"] = apparent_mag + calls["h"] = h + calls["k_correction"] = k_correction + calls["e_correction"] = e_correction + return np.array([-20.0, -19.0]) + + monkeypatch.setattr(cc, "absolute_magnitude", fake_absolute_magnitude) + + cosmo_obj = object() + result = cc.absolute_magnitude_limit( + cosmo_obj, + [0.1, 0.2], + m_lim=24.5, + h=0.7, + k_correction=0.1, + e_correction=0.2, + ) + + np.testing.assert_allclose(result, np.array([-20.0, -19.0])) + assert calls["cosmo_obj"] is cosmo_obj + np.testing.assert_allclose(calls["z"], np.array([0.1, 0.2])) + assert calls["apparent_mag"] == pytest.approx(24.5) + assert calls["h"] == pytest.approx(0.7) + assert calls["k_correction"] == pytest.approx(0.1) + assert calls["e_correction"] == pytest.approx(0.2) + + +def test_absolute_magnitude_limit_rejects_negative_redshift() -> None: + """Tests that negative redshifts are rejected.""" + with pytest.raises(ValueError, match="Redshift z must be >= 0"): + cc.absolute_magnitude_limit(object(), [-0.1, 0.2], m_lim=24.5) + + +def test_absolute_magnitude_limit_rejects_nonfinite_m_lim() -> None: + """Tests that non-finite apparent magnitude limits are rejected.""" + with pytest.raises(ValueError, match="m_lim must be finite"): + 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: + """Tests that observed density integrates from bright bound to catalog limit.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + result = cc.observed_number_density( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([3.0, 5.0])) + + +def test_observed_number_density_clips_catalog_limit_to_faint_bound( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that observed density clips limits fainter than the LF range.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-17.0, -16.0]), + ) + + result = cc.observed_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([6.0, 6.0])) + + +def test_missing_number_density_integrates_from_catalog_limit( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that missing density integrates from catalog limit to faint bound.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + result = cc.missing_number_density( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([3.0, 1.0])) + + +def test_missing_number_density_clips_catalog_limit_to_bright_bound( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that missing density clips limits brighter than the LF range.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-25.0, -26.0]), + ) + + result = cc.missing_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([6.0, 6.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"): + cc.observed_number_density( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-18.0, + m_faint=-24.0, + ) + + +def test_missing_number_density_rejects_invalid_magnitude_range() -> None: + """Tests that missing density rejects invalid LF magnitude bounds.""" + with pytest.raises(ValueError, match="m_faint must be larger than m_bright"): + cc.missing_number_density( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-18.0, + m_faint=-24.0, + ) + + +def test_catalog_completeness_fraction_returns_observed_fraction( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that catalog completeness returns observed over total density.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + result = cc.catalog_completeness_fraction( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.5, 5.0 / 6.0])) + + +def test_out_of_catalog_fraction_returns_missing_fraction( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that out-of-catalog fraction is one minus completeness.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + result = cc.out_of_catalog_fraction( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + 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( + 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 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, + ) + + +def test_catalog_completeness_fraction_returns_zero_for_zero_total_density( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that completeness is zero when total density is zero.""" + + def zero_lf(m_abs: np.ndarray, z: np.ndarray) -> np.ndarray: + return np.zeros_like(m_abs, dtype=float) + + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + result = cc.catalog_completeness_fraction( + object(), + [0.1, 0.2], + zero_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(result, np.array([0.0, 0.0])) + + +def test_completeness_and_out_of_catalog_fractions_sum_to_one( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """Tests that completeness and out-of-catalog fractions are complementary.""" + monkeypatch.setattr( + cc, + "absolute_magnitude_limit", + lambda *args, **kwargs: np.array([-21.0, -19.0]), + ) + + completeness = cc.catalog_completeness_fraction( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + missing = cc.out_of_catalog_fraction( + object(), + [0.1, 0.2], + constant_lf, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + ) + + np.testing.assert_allclose(completeness + missing, np.ones(2)) diff --git a/tests/test_photometry_completeness_fake_catalog.py b/tests/test_photometry_completeness_fake_catalog.py new file mode 100644 index 0000000..bbac27f --- /dev/null +++ b/tests/test_photometry_completeness_fake_catalog.py @@ -0,0 +1,251 @@ +"""Tests for LF completeness utilities using the packaged fake catalog.""" + +from __future__ import annotations + +from importlib.resources import files + +import numpy as np +import pandas as pd +import pyccl as ccl + +from lfkit.photometry.catalog_completeness import ( + catalog_completeness_fraction, + missing_number_density, + observed_number_density, + out_of_catalog_fraction, +) + + +def toy_luminosity_function( + absolute_mag: np.ndarray, + z: np.ndarray, +) -> np.ndarray: + """Return a smooth non-negative toy luminosity function.""" + phi0 = 1.0e-3 + m_star = -20.5 - 0.3 * z + width = 1.2 + + return phi0 * np.exp(-0.5 * ((absolute_mag - m_star) / width) ** 2) + + +def make_cosmology() -> ccl.Cosmology: + """Return a small test cosmology.""" + return ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + +def load_fake_catalog() -> pd.DataFrame: + """Return the packaged fake magnitude-limited catalog.""" + path = ( + files("lfkit") + / "data" + / "demo_catalogs" + / "fake_magnitude_limited_catalog.csv" + ) + return pd.read_csv(path) + + +def test_fake_catalog_is_packaged() -> None: + """Tests that the fake catalog is available as packaged data.""" + catalog = load_fake_catalog() + + assert len(catalog) == 200 + assert {"galaxy_id", "ra_deg", "dec_deg", "z", "m_app", "band"} <= set( + catalog.columns + ) + + +def test_completeness_fraction_on_fake_catalog_redshifts() -> None: + """Tests that completeness fractions are valid on fake catalog redshifts.""" + catalog = load_fake_catalog() + cosmo = make_cosmology() + + z = catalog["z"].to_numpy() + + completeness = catalog_completeness_fraction( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=256, + ) + + assert completeness.shape == z.shape + assert np.all(np.isfinite(completeness)) + assert np.all(completeness >= 0.0) + assert np.all(completeness <= 1.0) + + +def test_observed_and_missing_fractions_sum_to_one() -> None: + """Tests that observed and out-of-catalog fractions are complementary.""" + catalog = load_fake_catalog() + cosmo = make_cosmology() + + z = catalog["z"].to_numpy() + + completeness = catalog_completeness_fraction( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=256, + ) + missing_fraction = out_of_catalog_fraction( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=256, + ) + + np.testing.assert_allclose( + completeness + missing_fraction, + np.ones_like(z), + rtol=1.0e-12, + atol=1.0e-12, + ) + + +def test_deeper_catalog_limit_increases_completeness() -> None: + """Tests that a fainter apparent-magnitude limit increases completeness.""" + catalog = load_fake_catalog() + cosmo = make_cosmology() + + z = catalog["z"] + + shallow = catalog_completeness_fraction( + cosmo, + z, + toy_luminosity_function, + m_lim=22.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + deep = catalog_completeness_fraction( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + + np.testing.assert_array_less(shallow, deep + 1.0e-8) + assert np.mean(deep) > np.mean(shallow) + + +def test_deeper_catalog_limit_decreases_missing_density() -> None: + """Tests that a fainter apparent-magnitude limit misses no more galaxies.""" + catalog = load_fake_catalog() + cosmo = make_cosmology() + + z = catalog["z"] + + shallow = missing_number_density( + cosmo, + z, + toy_luminosity_function, + m_lim=22.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + deep = missing_number_density( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + + np.testing.assert_array_less(deep, shallow + 1.0e-8) + + +def test_deeper_catalog_limit_increases_observed_density() -> None: + """Tests that a fainter apparent-magnitude limit observes no fewer galaxies.""" + catalog = load_fake_catalog() + cosmo = make_cosmology() + + z = catalog["z"] + + shallow = observed_number_density( + cosmo, + z, + toy_luminosity_function, + m_lim=22.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + deep = observed_number_density( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + + np.testing.assert_array_less(shallow, deep + 1.0e-8) + assert np.mean(deep) > np.mean(shallow) + + +def test_observed_and_missing_number_densities_sum_to_total_density() -> None: + """Tests that observed and missing densities recover the total LF density.""" + catalog = load_fake_catalog() + cosmo = make_cosmology() + + z = catalog["z"].to_numpy() + + total = observed_number_density( + cosmo, + z, + toy_luminosity_function, + m_lim=1.0e6, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + observed = observed_number_density( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + missing = missing_number_density( + cosmo, + z, + toy_luminosity_function, + m_lim=24.0, + m_bright=-24.0, + m_faint=-14.0, + n_m=1024, + ) + + np.testing.assert_allclose( + observed + missing, + total, + rtol=1.0e-6, + atol=1.0e-12, + ) diff --git a/tests/test_photometry_lf_parameter_models.py b/tests/test_photometry_lf_parameter_models.py index 32658b1..b7017f8 100644 --- a/tests/test_photometry_lf_parameter_models.py +++ b/tests/test_photometry_lf_parameter_models.py @@ -5,7 +5,6 @@ import numpy as np import pytest -# Update this import path to match your actual module location. from lfkit.photometry.lf_parameter_models import ( ALPHA_MODELS, M_STAR_MODELS, @@ -18,6 +17,10 @@ m_star_linear_q, phi_star_constant, phi_star_linear_p, + available_lf_parameter_models, + register_alpha_model, + register_m_star_model, + register_phi_star_model, ) @@ -273,3 +276,192 @@ def test_scalar_input_returns_numpy_array() -> None: assert phi.dtype == np.float64 assert m_star.dtype == np.float64 assert alpha.dtype == np.float64 + + +def test_available_lf_parameter_models_returns_sorted_builtin_names() -> None: + """Tests that available_lf_parameter_models returns sorted built-in model names.""" + result = available_lf_parameter_models() + + assert result["phi_star"] == ["constant", "linear_p"] + assert result["m_star"] == ["constant", "linear_q"] + assert result["alpha"] == ["constant", "linear"] + + +def test_register_phi_star_model_adds_custom_model() -> None: + """Tests that register_phi_star_model adds a custom phi_star model.""" + def custom_phi(z: np.ndarray, *, amplitude: float) -> np.ndarray: + return np.full_like(z, amplitude, dtype=float) + + register_phi_star_model("custom_phi_test", custom_phi) + + try: + z = np.array([0.0, 0.5, 1.0]) + model = get_parameter_model( + "custom_phi_test", + PHI_STAR_MODELS, + model_kind="phi_model", + ) + result = model(z, amplitude=3.0e-3) + + np.testing.assert_allclose(result, np.array([3.0e-3, 3.0e-3, 3.0e-3])) + finally: + PHI_STAR_MODELS.pop("custom_phi_test", None) + + +def test_register_m_star_model_adds_custom_model() -> None: + """Tests that register_m_star_model adds a custom M_star model.""" + def custom_m_star(z: np.ndarray, *, base: float, slope: float) -> np.ndarray: + return np.asarray(base + slope * z, dtype=float) + + register_m_star_model("custom_m_star_test", custom_m_star) + + try: + z = np.array([0.0, 0.5, 1.0]) + model = get_parameter_model( + "custom_m_star_test", + M_STAR_MODELS, + model_kind="m_star_model", + ) + result = model(z, base=-20.0, slope=-1.0) + + np.testing.assert_allclose(result, np.array([-20.0, -20.5, -21.0])) + finally: + M_STAR_MODELS.pop("custom_m_star_test", None) + + +def test_register_alpha_model_adds_custom_model() -> None: + """Tests that register_alpha_model adds a custom alpha model.""" + def custom_alpha(z: np.ndarray, *, alpha: float) -> np.ndarray: + return np.full_like(z, alpha, dtype=float) + + register_alpha_model("custom_alpha_test", custom_alpha) + + try: + z = np.array([0.0, 0.5, 1.0]) + model = get_parameter_model( + "custom_alpha_test", + ALPHA_MODELS, + model_kind="alpha_model", + ) + result = model(z, alpha=-1.4) + + np.testing.assert_allclose(result, np.array([-1.4, -1.4, -1.4])) + finally: + ALPHA_MODELS.pop("custom_alpha_test", None) + + +def test_evaluate_lf_parameters_uses_registered_custom_models() -> None: + """Tests that evaluate_lf_parameters can use registered custom models.""" + def custom_phi(z: np.ndarray, *, amplitude: float) -> np.ndarray: + return np.full_like(z, amplitude, dtype=float) + + def custom_m_star(z: np.ndarray, *, base: float) -> np.ndarray: + return np.asarray(base - z, dtype=float) + + def custom_alpha(z: np.ndarray, *, alpha: float) -> np.ndarray: + return np.full_like(z, alpha, dtype=float) + + register_phi_star_model("eval_custom_phi_test", custom_phi) + register_m_star_model("eval_custom_m_star_test", custom_m_star) + register_alpha_model("eval_custom_alpha_test", custom_alpha) + + try: + z = np.array([0.0, 0.5, 1.0]) + + phi_star, m_star, alpha = evaluate_lf_parameters( + z, + phi_model="eval_custom_phi_test", + phi_kwargs={"amplitude": 2.0e-3}, + m_star_model="eval_custom_m_star_test", + m_star_kwargs={"base": -20.0}, + alpha_model="eval_custom_alpha_test", + alpha_kwargs={"alpha": -1.2}, + ) + + np.testing.assert_allclose(phi_star, np.array([2.0e-3, 2.0e-3, 2.0e-3])) + np.testing.assert_allclose(m_star, np.array([-20.0, -20.5, -21.0])) + np.testing.assert_allclose(alpha, np.array([-1.2, -1.2, -1.2])) + finally: + PHI_STAR_MODELS.pop("eval_custom_phi_test", None) + M_STAR_MODELS.pop("eval_custom_m_star_test", None) + ALPHA_MODELS.pop("eval_custom_alpha_test", None) + + +def test_register_parameter_model_raises_for_empty_name() -> None: + """Tests that registering a model with an empty name raises an error.""" + def custom_model(z: np.ndarray) -> np.ndarray: + return np.asarray(z, dtype=float) + + with pytest.raises(ValueError, match="phi_star model name cannot be empty"): + register_phi_star_model("", custom_model) + + +def test_register_parameter_model_raises_for_non_callable_model() -> None: + """Tests that registering a non-callable model raises an error.""" + with pytest.raises(TypeError, match="phi_star model must be callable"): + register_phi_star_model("not_callable_test", 3.0) # type: ignore[arg-type] + + +def test_register_parameter_model_raises_for_duplicate_without_overwrite() -> None: + """Tests that registering a duplicate model without overwrite raises an error.""" + def custom_model(z: np.ndarray) -> np.ndarray: + return np.asarray(z, dtype=float) + + register_phi_star_model("duplicate_phi_test", custom_model) + + try: + with pytest.raises(ValueError, match="already registered"): + register_phi_star_model("duplicate_phi_test", custom_model) + finally: + PHI_STAR_MODELS.pop("duplicate_phi_test", None) + + +def test_register_parameter_model_allows_duplicate_with_overwrite() -> None: + """Tests that registering a duplicate model with overwrite replaces the model.""" + def first_model(z: np.ndarray) -> np.ndarray: + return np.full_like(z, 1.0, dtype=float) + + def second_model(z: np.ndarray) -> np.ndarray: + return np.full_like(z, 2.0, dtype=float) + + register_phi_star_model("overwrite_phi_test", first_model) + + try: + register_phi_star_model("overwrite_phi_test", second_model, overwrite=True) + + z = np.array([0.0, 0.5]) + result = PHI_STAR_MODELS["overwrite_phi_test"](z) + + np.testing.assert_allclose(result, np.array([2.0, 2.0])) + finally: + PHI_STAR_MODELS.pop("overwrite_phi_test", None) + + +def test_parameter_models_raise_for_non_finite_redshift() -> None: + """Tests that parameter models reject non-finite redshift values.""" + z = np.array([0.0, np.nan, 1.0]) + + with pytest.raises(ValueError, match="z contains NaN or infinite values"): + phi_star_constant(z, phi_star=1.0e-3) + + with pytest.raises(ValueError, match="z contains NaN or infinite values"): + m_star_constant(z, m_star=-20.0) + + with pytest.raises(ValueError, match="z contains NaN or infinite values"): + alpha_constant(z, alpha=-1.0) + + +def test_evaluate_lf_parameters_rejects_non_finite_redshift() -> None: + """Tests that evaluate_lf_parameters rejects non-finite redshift values.""" + z = np.array([0.0, np.inf]) + + with pytest.raises(ValueError, match="z contains NaN or infinite values"): + evaluate_lf_parameters( + z, + phi_model="constant", + phi_kwargs={"phi_star": 1.0e-3}, + m_star_model="constant", + m_star_kwargs={"m_star": -20.0}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.0}, + ) diff --git a/tests/test_photometry_magnitudes.py b/tests/test_photometry_magnitudes.py index f214475..9b33456 100644 --- a/tests/test_photometry_magnitudes.py +++ b/tests/test_photometry_magnitudes.py @@ -25,15 +25,15 @@ def test_total_magnitude_correction_k_only_scalar(): def test_total_magnitude_correction_e_only_scalar(): - """Tests that total_magnitude_correction returns only the e-correction for scalar input.""" + """Tests that total_magnitude_correction returns minus the e-correction for scalar input.""" out = magnitudes.total_magnitude_correction(e_correction=0.2) assert isinstance(out, np.ndarray) assert out.dtype == float - assert np.allclose(out, 0.2) + assert np.allclose(out, -0.2) def test_total_magnitude_correction_adds_array_terms(): - """Tests that total_magnitude_correction adds k- and e-corrections elementwise.""" + """Tests that total_magnitude_correction combines k- and e-corrections elementwise.""" k_corr = np.array([0.1, 0.2, 0.3]) e_corr = np.array([0.05, 0.10, 0.15]) @@ -42,7 +42,7 @@ def test_total_magnitude_correction_adds_array_terms(): e_correction=e_corr, ) - expected = np.array([0.15, 0.30, 0.45]) + expected = np.array([0.05, 0.10, 0.15]) assert isinstance(out, np.ndarray) assert out.dtype == float assert np.allclose(out, expected) @@ -54,13 +54,14 @@ def test_total_magnitude_correction_broadcasts_scalar_and_array(): k_correction=0.2, e_correction=np.array([0.0, 0.1, 0.2]), ) - expected = np.array([0.2, 0.3, 0.4]) + expected = np.array([0.2, 0.1, 0.0]) assert np.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): + _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) monkeypatch.setattr(magnitudes, "distance_modulus", mock_distance_modulus) @@ -78,8 +79,9 @@ def mock_distance_modulus(cosmo_obj, z, h=None): def test_absolute_magnitude_with_corrections(monkeypatch: pytest.MonkeyPatch): - """Tests that absolute_magnitude subtracts distance modulus and total correction.""" + """Tests that absolute_magnitude applies M = m - mu - K + E.""" def mock_distance_modulus(cosmo_obj, z, h=None): + _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) monkeypatch.setattr(magnitudes, "distance_modulus", mock_distance_modulus) @@ -92,13 +94,19 @@ def mock_distance_modulus(cosmo_obj, z, h=None): e_correction=np.array([0.3, 0.4]), ) - expected = np.array([20.0, 21.5]) - np.array([40.0, 41.0]) - np.array([0.4, 0.6]) + expected = ( + np.array([20.0, 21.5]) + - np.array([40.0, 41.0]) + - np.array([0.1, 0.2]) + + np.array([0.3, 0.4]) + ) assert np.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): + _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) monkeypatch.setattr(magnitudes, "distance_modulus", mock_distance_modulus) @@ -116,8 +124,9 @@ def mock_distance_modulus(cosmo_obj, z, h=None): def test_apparent_magnitude_with_corrections(monkeypatch: pytest.MonkeyPatch): - """Tests that apparent_magnitude adds distance modulus and total correction.""" + """Tests that apparent_magnitude applies m = M + mu + K - E.""" def mock_distance_modulus(cosmo_obj, z, h=None): + _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0], dtype=float) monkeypatch.setattr(magnitudes, "distance_modulus", mock_distance_modulus) @@ -130,7 +139,12 @@ def mock_distance_modulus(cosmo_obj, z, h=None): e_correction=np.array([0.3, 0.4]), ) - expected = np.array([-20.0, -19.5]) + np.array([40.0, 41.0]) + np.array([0.4, 0.6]) + expected = ( + np.array([-20.0, -19.5]) + + np.array([40.0, 41.0]) + + np.array([0.1, 0.2]) + - np.array([0.3, 0.4]) + ) assert np.allclose(out, expected) @@ -139,6 +153,7 @@ def test_apparent_and_absolute_are_inverse_without_corrections( ): """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 z_arr = np.asarray(z, dtype=float) return 40.0 + z_arr @@ -166,6 +181,7 @@ def test_apparent_and_absolute_are_inverse_with_corrections( ): """Tests that apparent_magnitude and absolute_magnitude invert each other with corrections.""" def mock_distance_modulus(cosmo_obj, z, h=None): + _, _ = cosmo_obj, h z_arr = np.asarray(z, dtype=float) return 39.5 + 2.0 * z_arr @@ -201,6 +217,7 @@ def test_absolute_magnitude_passes_h_to_distance_modulus( captured: dict[str, float | None] = {"h": None} def mock_distance_modulus(cosmo_obj, z, h=None): + _, _ = cosmo_obj, z captured["h"] = h return np.array([40.0], dtype=float) @@ -223,6 +240,7 @@ def test_apparent_magnitude_passes_h_to_distance_modulus( captured: dict[str, float | None] = {"h": None} def mock_distance_modulus(cosmo_obj, z, h=None): + _, _ = cosmo_obj, z captured["h"] = h return np.array([40.0], dtype=float) @@ -243,6 +261,7 @@ def test_absolute_magnitude_broadcasts_scalar_correction( ): """Tests that absolute_magnitude broadcasts scalar corrections over array magnitudes.""" def mock_distance_modulus(cosmo_obj, z, h=None): + _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0, 42.0], dtype=float) monkeypatch.setattr(magnitudes, "distance_modulus", mock_distance_modulus) @@ -263,6 +282,7 @@ def test_apparent_magnitude_broadcasts_scalar_correction( ): """Tests that apparent_magnitude broadcasts scalar corrections over array magnitudes.""" def mock_distance_modulus(cosmo_obj, z, h=None): + _, _, _ = cosmo_obj, z, h return np.array([40.0, 41.0, 42.0], dtype=float) monkeypatch.setattr(magnitudes, "distance_modulus", mock_distance_modulus) @@ -274,5 +294,5 @@ def mock_distance_modulus(cosmo_obj, z, h=None): e_correction=0.3, ) - expected = np.array([-20.0, -20.0, -20.0]) + np.array([40.0, 41.0, 42.0]) + 0.3 + expected = np.array([-20.0, -20.0, -20.0]) + np.array([40.0, 41.0, 42.0]) - 0.3 assert np.allclose(out, expected) From 977ef9a4531e7bdd7b7ed55ac3c46c830ec50aa1 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:25:57 -0400 Subject: [PATCH 05/10] Expose luminosity-function completeness API --- src/lfkit/__init__.py | 10 +- src/lfkit/api/lumfunc.py | 657 ++++++++++++++++++++++++++++++++++++++ tests/test_api_lumfunc.py | 532 ++++++++++++++++++++++++++++++ 3 files changed, 1198 insertions(+), 1 deletion(-) create mode 100644 tests/test_api_lumfunc.py diff --git a/src/lfkit/__init__.py b/src/lfkit/__init__.py index c9dd3d1..00b80a6 100644 --- a/src/lfkit/__init__.py +++ b/src/lfkit/__init__.py @@ -1,5 +1,13 @@ +"""Top-level package for LFKit.""" + from __future__ import annotations from lfkit.api.corrections import Corrections +from lfkit.api.lumfunc import LuminosityFunction + +__all__ = [ + "Corrections", + "LuminosityFunction", +] -__all__ = ["Corrections"] +__author__ = """Niko Sarcevic and collaborators.""" diff --git a/src/lfkit/api/lumfunc.py b/src/lfkit/api/lumfunc.py index e69de29..c4e6e11 100644 --- a/src/lfkit/api/lumfunc.py +++ b/src/lfkit/api/lumfunc.py @@ -0,0 +1,657 @@ +r"""Public luminosity-function interface. + +This module provides the user-facing :class:`LuminosityFunction` API for +evaluating luminosity functions in absolute- or apparent-magnitude space. + +The class wraps the lower-level LFKit photometry functions behind a small, +stable interface. Users can construct standard, evolving, or double Schechter +models, evaluate :math:`\phi(M, z)`, evaluate from apparent magnitudes, and +compute integrated, observed, and missing number densities for +magnitude-limited catalog selections. + +File reading is intentionally not handled here. Catalog-derived LF parameters, +magnitude limits, or correction models should be loaded elsewhere and passed +into this API as scalars, arrays, or correction objects. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from typing import TYPE_CHECKING + +import numpy as np + +from lfkit.photometry.luminosity_function import ( + schechter, + schechter_evolving, + schechter_double, + schechter_from_m, + schechter_evolving_from_m, + schechter_double_from_m, +) +from lfkit.photometry.magnitudes import ( + absolute_magnitude, + apparent_magnitude, +) +from lfkit.photometry.lf_parameter_models import ( + available_lf_parameter_models, + evaluate_lf_parameters, + register_alpha_model, + register_m_star_model, + register_phi_star_model, +) +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, + FloatArray, + FloatInput, + ParameterModel, + ParameterValue, +) + +if TYPE_CHECKING: + from lfkit.api.corrections import Corrections +else: + Corrections = object + + +__all__ = ["LuminosityFunction"] + + +class LuminosityFunction: + """User-facing wrapper for luminosity-function evaluation.""" + + def __init__( + self, + *, + model: str, + parameters: Mapping[str, object], + meta: Mapping[str, object] | None = None, + ) -> None: + """Store a luminosity-function model and its parameters. + + Args: + model: Name of the luminosity-function model. + parameters: Model parameters passed to the underlying LF function. + meta: Optional metadata describing the LF source or calibration. + """ + self.model = str(model) + self.parameters_dict = dict(parameters) + self.meta = {} if meta is None else dict(meta) + + def absolute_magnitude( + self, + cosmo_obj: Cosmology, + z: FloatInput, + apparent_mag: FloatInput, + *, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Convert apparent magnitude to absolute magnitude. + + Args: + cosmo_obj: Cosmology object used for distance-modulus conversion. + z: Redshift values. + apparent_mag: Apparent magnitude values. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Absolute magnitudes using the LFKit convention + ``M = m - mu - K + E``. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return absolute_magnitude( + cosmo_obj, + z, + apparent_mag, + h=h, + k_correction=k_corr, + e_correction=e_corr, + ) + + def apparent_magnitude( + self, + cosmo_obj: Cosmology, + z: FloatInput, + absolute_mag: FloatInput, + *, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Convert absolute magnitude to apparent magnitude. + + Args: + cosmo_obj: Cosmology object used for distance-modulus conversion. + z: Redshift values. + absolute_mag: Absolute magnitude values. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Apparent magnitudes using the LFKit convention + ``m = M + mu + K - E``. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return apparent_magnitude( + cosmo_obj, + z, + absolute_mag, + h=h, + k_correction=k_corr, + e_correction=e_corr, + ) + + @classmethod + def schechter( + cls, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, + ) -> "LuminosityFunction": + """Create a standard Schechter luminosity function. + + Args: + phi_star: Normalization of the luminosity function. + m_star: Characteristic absolute magnitude. + alpha: Faint-end slope. + + Returns: + Luminosity-function API object using the standard Schechter model. + """ + return cls( + model="schechter", + parameters={ + "phi_star": phi_star, + "m_star": m_star, + "alpha": alpha, + }, + ) + + @classmethod + def evolving_schechter( + cls, + *, + phi_model: str = "linear_p", + phi_kwargs: Mapping[str, ParameterValue] | None = None, + m_star_model: str = "linear_q", + m_star_kwargs: Mapping[str, ParameterValue] | None = None, + alpha_model: str = "constant", + alpha_kwargs: Mapping[str, ParameterValue] | None = None, + ) -> "LuminosityFunction": + """Create a redshift-evolving Schechter luminosity function. + + Args: + phi_model: Parameter model used for the normalization evolution. + phi_kwargs: Keyword arguments for the normalization model. + m_star_model: Parameter model used for characteristic-magnitude evolution. + m_star_kwargs: Keyword arguments for the characteristic-magnitude model. + alpha_model: Parameter model used for faint-end-slope evolution. + alpha_kwargs: Keyword arguments for the faint-end-slope model. + + Returns: + Luminosity-function API object using an evolving Schechter model. + """ + return cls( + model="evolving_schechter", + parameters={ + "phi_model": phi_model, + "phi_kwargs": {} if phi_kwargs is None else dict(phi_kwargs), + "m_star_model": m_star_model, + "m_star_kwargs": {} if m_star_kwargs is None else dict(m_star_kwargs), + "alpha_model": alpha_model, + "alpha_kwargs": {} if alpha_kwargs is None else dict(alpha_kwargs), + }, + ) + + @classmethod + def double_schechter( + cls, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: float, + beta: float, + m_transition: ParameterValue, + ) -> "LuminosityFunction": + """Create a double-power-law Schechter luminosity function. + + Args: + phi_star: Normalization of the luminosity function. + m_star: Characteristic absolute magnitude. + alpha: Bright-end or main Schechter slope. + beta: Additional slope controlling the second power-law component. + m_transition: Transition magnitude for the second component. + + Returns: + Luminosity-function API object using the double Schechter model. + """ + return cls( + model="double_schechter", + parameters={ + "phi_star": phi_star, + "m_star": m_star, + "alpha": alpha, + "beta": beta, + "m_transition": m_transition, + }, + ) + + def phi( + self, + absolute_mag: FloatInput, + z: FloatInput | None = None, + ) -> FloatArray: + """Evaluate the luminosity function in absolute-magnitude space. + + Args: + absolute_mag: Absolute magnitude values where the LF is evaluated. + z: Redshift values. Required for evolving Schechter models. + + Returns: + Luminosity-function values evaluated at the input magnitudes. + """ + if self.model == "schechter": + return schechter( + np.asarray(absolute_mag, dtype=float), + **self.parameters_dict, + ) + + if self.model == "evolving_schechter": + if z is None: + raise ValueError("z is required for an evolving luminosity function.") + + return schechter_evolving( + np.asarray(absolute_mag, dtype=float), + np.asarray(z, dtype=float), + **self.parameters_dict, + ) + + if self.model == "double_schechter": + return schechter_double( + np.asarray(absolute_mag, dtype=float), + **self.parameters_dict, + ) + + raise ValueError(f"Unsupported luminosity-function model '{self.model}'.") + + def phi_from_m( + self, + cosmo_obj: Cosmology, + z: FloatInput, + apparent_mag: FloatInput, + *, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Evaluate the luminosity function from apparent magnitudes. + + Apparent magnitudes are converted to absolute magnitudes using the + supplied cosmology, optional reduced Hubble parameter, and optional + k- and e-correction model. + + Args: + cosmo_obj: Cosmology object used for distance-modulus conversion. + z: Redshift values. + apparent_mag: Apparent magnitude values. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Luminosity-function values evaluated from apparent magnitudes. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + if self.model == "schechter": + return schechter_from_m( + cosmo_obj, + np.asarray(z, dtype=float), + np.asarray(apparent_mag, dtype=float), + h=h, + k_correction=k_corr, + e_correction=e_corr, + **self.parameters_dict, + ) + + if self.model == "evolving_schechter": + return schechter_evolving_from_m( + cosmo_obj, + np.asarray(z, dtype=float), + np.asarray(apparent_mag, dtype=float), + h=h, + k_correction=k_corr, + e_correction=e_corr, + **self.parameters_dict, + ) + + if self.model == "double_schechter": + return schechter_double_from_m( + cosmo_obj, + np.asarray(z, dtype=float), + np.asarray(apparent_mag, dtype=float), + h=h, + k_correction=k_corr, + e_correction=e_corr, + **self.parameters_dict, + ) + + raise ValueError(f"Unsupported luminosity-function model '{self.model}'.") + + def parameters( + self, + z: FloatInput, + ) -> tuple[FloatArray, FloatArray, FloatArray]: + """Evaluate evolving Schechter parameters at redshift. + + Args: + z: Redshift values where the evolving LF parameters are evaluated. + + Returns: + Tuple containing ``phi_star(z)``, ``m_star(z)``, and ``alpha(z)``. + """ + if self.model != "evolving_schechter": + raise ValueError("parameters(z) is only defined for evolving_schechter.") + + return evaluate_lf_parameters( + np.asarray(z, dtype=float), + **self.parameters_dict, + ) + + def absolute_magnitude_limit( + self, + cosmo_obj: Cosmology, + z: FloatInput, + *, + m_lim: float, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Return the absolute-magnitude limit of a catalog apparent-magnitude cut. + + Args: + cosmo_obj: Cosmology object used for distance-modulus conversion. + z: Redshift values. + m_lim: Apparent-magnitude limit of the catalog. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Absolute-magnitude limits using the LFKit convention + ``M_lim = m_lim - mu - K + E``. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return absolute_magnitude_limit( + cosmo_obj, + z, + m_lim=m_lim, + h=h, + k_correction=k_corr, + e_correction=e_corr, + ) + + @staticmethod + def available_parameter_models() -> dict[str, list[str]]: + """Return available LF parameter evolution models.""" + return available_lf_parameter_models() + + @staticmethod + def register_phi_star_model( + name: str, + model: ParameterModel, + *, + overwrite: bool = False, + ) -> None: + """Register a phi_star evolution model.""" + register_phi_star_model(name, model, overwrite=overwrite) + + @staticmethod + def register_m_star_model( + name: str, + model: ParameterModel, + *, + overwrite: bool = False, + ) -> None: + """Register an M_star evolution model.""" + register_m_star_model(name, model, overwrite=overwrite) + + @staticmethod + def register_alpha_model( + name: str, + model: ParameterModel, + *, + overwrite: bool = False, + ) -> None: + """Register an alpha evolution model.""" + register_alpha_model(name, model, overwrite=overwrite) + + def integrated_number_density( + self, + z: FloatInput, + *, + m_bright: ParameterValue, + m_faint: ParameterValue, + n_m: int = 512, + ) -> FloatArray: + """Integrate the LF over an absolute-magnitude range. + + Args: + z: Redshift values. + m_bright: Bright absolute-magnitude integration limit. + m_faint: Faint absolute-magnitude integration limit. + n_m: Number of magnitude-grid points used in the integration. + + Returns: + Number density integrated over the requested magnitude range. + """ + return integrated_number_density( + z, + self._as_callable(), + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + ) + + def observed_number_density( + self, + cosmo_obj: Cosmology, + z: FloatInput, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Return the LF number density observable in a magnitude-limited catalog. + + Args: + cosmo_obj: Cosmology object used for apparent-to-absolute conversion. + z: Redshift values. + m_lim: Apparent-magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration limit. + m_faint: Faint absolute-magnitude integration limit. + n_m: Number of magnitude-grid points used in the integration. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Number density brighter than the catalog magnitude limit. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return observed_number_density( + cosmo_obj, + z, + self._as_callable(), + m_lim=m_lim, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + h=h, + k_correction=k_corr, + e_correction=e_corr, + ) + + def missing_number_density( + self, + cosmo_obj: Cosmology, + z: FloatInput, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Return the LF number density missing from a magnitude-limited catalog. + + Args: + cosmo_obj: Cosmology object used for apparent-to-absolute conversion. + z: Redshift values. + m_lim: Apparent-magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration limit. + m_faint: Faint absolute-magnitude integration limit. + n_m: Number of magnitude-grid points used in the integration. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Number density fainter than the catalog magnitude limit but inside + the requested absolute-magnitude range. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return missing_number_density( + cosmo_obj, + z, + self._as_callable(), + m_lim=m_lim, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + h=h, + k_correction=k_corr, + e_correction=e_corr, + ) + + def catalog_completeness( + self, + cosmo_obj: Cosmology, + z: FloatInput, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Return the observed LF fraction in a magnitude-limited catalog. + + Args: + cosmo_obj: Cosmology object used for apparent-to-absolute conversion. + z: Redshift values. + m_lim: Apparent-magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration limit. + m_faint: Faint absolute-magnitude integration limit. + n_m: Number of magnitude-grid points used in the integration. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Fraction of the LF number density observable in the catalog. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return catalog_completeness_fraction( + cosmo_obj, + z, + self._as_callable(), + m_lim=m_lim, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + h=h, + k_correction=k_corr, + e_correction=e_corr, + ) + + def out_of_catalog_fraction( + self, + cosmo_obj: Cosmology, + z: FloatInput, + *, + m_lim: float, + m_bright: float, + m_faint: float, + n_m: int = 512, + h: float | None = None, + corrections: Corrections | None = None, + ) -> FloatArray: + """Return the missing LF fraction for a magnitude-limited catalog. + + Args: + cosmo_obj: Cosmology object used for apparent-to-absolute conversion. + z: Redshift values. + m_lim: Apparent-magnitude limit of the catalog. + m_bright: Bright absolute-magnitude integration limit. + m_faint: Faint absolute-magnitude integration limit. + n_m: Number of magnitude-grid points used in the integration. + h: Optional reduced Hubble parameter used in the magnitude conversion. + corrections: Optional object providing k-correction and e-correction values. + + Returns: + Fraction of the LF number density missing from the catalog. + """ + k_corr, e_corr = self._correction_values(corrections, z) + + return out_of_catalog_fraction( + cosmo_obj, + z, + self._as_callable(), + m_lim=m_lim, + m_bright=m_bright, + m_faint=m_faint, + n_m=n_m, + h=h, + k_correction=k_corr, + e_correction=e_corr, + ) + + def _as_callable(self): + """Return this object as an ``lf(M, z)`` callable.""" + return lambda absolute_mag, z: self.phi(absolute_mag, z) + + @staticmethod + def _correction_values( + corrections: Corrections | None, + z: FloatInput, + ) -> tuple[FloatArray | None, FloatArray | None]: + """Evaluate optional correction values at redshift. + + Args: + corrections: Optional correction object with ``k(z)`` and ``e(z)`` methods. + z: Redshift values where corrections are evaluated. + + Returns: + Tuple of k-correction and e-correction arrays, or ``None`` values + when no correction object is supplied. + """ + if corrections is None: + return None, None + + return corrections.k(z), corrections.e(z) diff --git a/tests/test_api_lumfunc.py b/tests/test_api_lumfunc.py new file mode 100644 index 0000000..1f830ba --- /dev/null +++ b/tests/test_api_lumfunc.py @@ -0,0 +1,532 @@ +"""Unit tests for ``lfkit.api.lumfunc.py``.""" + +from __future__ import annotations + +from typing import cast + +import numpy as np +import pytest + +import lfkit.api.lumfunc as lf_api +from lfkit.api.corrections import Corrections +from lfkit.api.lumfunc import LuminosityFunction + + +class DummyCorrections: + """Small correction object used to test public API forwarding.""" + + def k(self, z): + """Tests that k-corrections can be evaluated at z.""" + return np.asarray(z, dtype=float) + 1.0 + + def e(self, z): + """Tests that e-corrections can be evaluated at z.""" + return np.asarray(z, dtype=float) - 1.0 + + +def test_schechter_constructor_stores_model_and_parameters(): + """Tests that the Schechter constructor stores the expected public state.""" + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + assert lf.model == "schechter" + assert lf.parameters_dict == { + "phi_star": 1.0e-3, + "m_star": -20.5, + "alpha": -1.2, + } + assert lf.meta == {} + + +def test_constructor_copies_parameter_and_metadata_mappings(): + """Tests that input mappings are copied rather than stored by reference.""" + parameters = {"phi_star": 1.0e-3, "m_star": -20.5, "alpha": -1.2} + meta = {"survey": "test"} + + lf = LuminosityFunction( + model="schechter", + parameters=parameters, + meta=meta, + ) + + parameters["phi_star"] = 9.0 + meta["survey"] = "changed" + + assert lf.parameters_dict["phi_star"] == 1.0e-3 + assert lf.meta["survey"] == "test" + + +def test_phi_dispatches_schechter_model(monkeypatch): + """Tests that phi dispatches to the standard Schechter implementation.""" + absolute_mag = np.array([-21.0, -20.0]) + + def fake_schechter(mag, *, phi_star, m_star, alpha): + assert np.allclose(mag, absolute_mag) + assert phi_star == 1.0e-3 + assert m_star == -20.5 + assert alpha == -1.2 + return np.array([1.0, 2.0]) + + monkeypatch.setattr(lf_api, "schechter", fake_schechter) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.phi(absolute_mag) + + assert np.allclose(result, [1.0, 2.0]) + + +def test_phi_requires_redshift_for_evolving_schechter(): + """Tests that evolving Schechter evaluation requires redshift input.""" + lf = LuminosityFunction.evolving_schechter() + + with pytest.raises(ValueError, match="z is required"): + lf.phi(np.array([-21.0, -20.0])) + + +def test_phi_dispatches_evolving_schechter_model(monkeypatch): + """Tests that phi dispatches to the evolving Schechter implementation.""" + absolute_mag = np.array([-21.0, -20.0]) + z = np.array([0.2, 0.8]) + + def fake_schechter_evolving( + mag, + redshift, + *, + phi_model, + phi_kwargs, + m_star_model, + m_star_kwargs, + alpha_model, + alpha_kwargs, + ): + assert np.allclose(mag, absolute_mag) + assert np.allclose(redshift, z) + 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([3.0, 4.0]) + + monkeypatch.setattr(lf_api, "schechter_evolving", fake_schechter_evolving) + + 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(absolute_mag, z) + + assert np.allclose(result, [3.0, 4.0]) + + +def test_phi_dispatches_double_schechter_model(monkeypatch): + """Tests that phi dispatches to the double Schechter implementation.""" + absolute_mag = np.array([-21.0, -20.0]) + + def fake_schechter_double( + mag, + *, + phi_star, + m_star, + alpha, + beta, + m_transition, + ): + assert np.allclose(mag, absolute_mag) + 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([5.0, 6.0]) + + monkeypatch.setattr(lf_api, "schechter_double", fake_schechter_double) + + 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(absolute_mag) + + assert np.allclose(result, [5.0, 6.0]) + + +def test_phi_raises_for_unsupported_model(): + """Tests that unsupported luminosity-function models fail clearly.""" + lf = LuminosityFunction(model="bad_model", parameters={}) + + with pytest.raises(ValueError, match="Unsupported luminosity-function model"): + lf.phi(np.array([-21.0, -20.0])) + + +def test_phi_from_m_forwards_corrections_to_schechter_from_m(monkeypatch): + """Tests that phi_from_m forwards correction arrays to magnitude evaluation.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + apparent_mag = np.array([22.0, 24.0]) + corrections = cast(Corrections, DummyCorrections()) + + def fake_schechter_from_m( + received_cosmo, + received_z, + received_apparent_mag, + *, + h, + k_correction, + e_correction, + phi_star, + m_star, + alpha, + ): + 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]) + assert phi_star == 1.0e-3 + assert m_star == -20.5 + assert alpha == -1.2 + return np.array([7.0, 8.0]) + + monkeypatch.setattr(lf_api, "schechter_from_m", fake_schechter_from_m) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.phi_from_m( + cosmo_obj, + z, + apparent_mag, + h=0.7, + corrections=corrections, + ) + + assert np.allclose(result, [7.0, 8.0]) + + +def test_phi_from_m_passes_none_corrections_when_not_supplied(monkeypatch): + """Tests that phi_from_m uses None corrections when no correction object is given.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + apparent_mag = np.array([22.0, 24.0]) + + def fake_schechter_from_m( + received_cosmo, + received_z, + received_apparent_mag, + *, + h, + k_correction, + e_correction, + phi_star, + m_star, + alpha, + ): + 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 + return np.array([1.0, 1.0]) + + monkeypatch.setattr(lf_api, "schechter_from_m", fake_schechter_from_m) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + result = lf.phi_from_m(cosmo_obj, z, apparent_mag) + + assert np.allclose(result, [1.0, 1.0]) + + +def test_parameters_only_work_for_evolving_schechter(monkeypatch): + """Tests that parameters delegates only for evolving Schechter models.""" + z = np.array([0.1, 0.5]) + + def fake_evaluate_lf_parameters(redshift, **kwargs): + assert np.allclose(redshift, z) + assert kwargs["phi_model"] == "constant" + return ( + np.array([1.0e-3, 1.0e-3]), + np.array([-20.5, -20.5]), + np.array([-1.2, -1.2]), + ) + + monkeypatch.setattr(lf_api, "evaluate_lf_parameters", fake_evaluate_lf_parameters) + + 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}, + ) + + phi_star, m_star, alpha = lf.parameters(z) + + assert np.allclose(phi_star, [1.0e-3, 1.0e-3]) + assert np.allclose(m_star, [-20.5, -20.5]) + assert np.allclose(alpha, [-1.2, -1.2]) + + +def test_parameters_raise_for_non_evolving_schechter(): + """Tests that parameters raises for non-evolving luminosity functions.""" + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + with pytest.raises(ValueError, match="only defined for evolving_schechter"): + lf.parameters(np.array([0.1, 0.5])) + + +def test_integrated_number_density_uses_api_callable(monkeypatch): + """Tests that integrated number density evaluates the public LF callable.""" + z = np.array([0.1, 0.5]) + + def fake_integrated_number_density(redshift, lf_callable, *, m_bright, m_faint, n_m): + mag = np.array([-21.0, -20.0]) + assert np.allclose(redshift, z) + assert m_bright == -24.0 + assert m_faint == -18.0 + assert n_m == 32 + assert np.allclose(lf_callable(mag, z), [2.0, 2.0]) + return np.array([10.0, 20.0]) + + monkeypatch.setattr( + lf_api, + "integrated_number_density", + fake_integrated_number_density, + ) + + lf = LuminosityFunction( + model="schechter", + parameters={"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_number_density( + z, + m_bright=-24.0, + m_faint=-18.0, + n_m=32, + ) + + assert np.allclose(result, [10.0, 20.0]) + + +def test_catalog_completeness_forwards_corrections(monkeypatch): + """Tests that catalog completeness forwards correction arrays.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + corrections = cast(Corrections, DummyCorrections()) + + def fake_catalog_completeness_fraction( + received_cosmo, + received_z, + lf_callable, + *, + m_lim, + m_bright, + m_faint, + n_m, + h, + k_correction, + e_correction, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert m_lim == 24.5 + assert m_bright == -24.0 + assert m_faint == -18.0 + assert n_m == 64 + assert h == 0.7 + assert np.allclose(k_correction, [1.1, 1.5]) + assert np.allclose(e_correction, [-0.9, -0.5]) + assert np.allclose(lf_callable(np.array([-21.0]), np.array([0.1])), [3.0]) + return np.array([0.8, 0.6]) + + monkeypatch.setattr( + lf_api, + "catalog_completeness_fraction", + fake_catalog_completeness_fraction, + ) + + 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, 3.0), + ) + + result = lf.catalog_completeness( + cosmo_obj, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=64, + h=0.7, + corrections=corrections, + ) + + assert np.allclose(result, [0.8, 0.6]) + + +def test_observed_and_missing_number_density_are_consistent(monkeypatch): + """Tests that observed and missing wrappers forward matching inputs.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + + def fake_observed_number_density( + received_cosmo, + received_z, + lf_callable, + *, + m_lim, + m_bright, + m_faint, + n_m, + h, + k_correction, + e_correction, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert m_lim == 24.5 + assert m_bright == -24.0 + assert m_faint == -18.0 + assert n_m == 128 + assert h is None + assert k_correction is None + assert e_correction is None + return np.array([4.0, 6.0]) + + def fake_missing_number_density( + received_cosmo, + received_z, + lf_callable, + *, + m_lim, + m_bright, + m_faint, + n_m, + h, + k_correction, + e_correction, + ): + assert received_cosmo is cosmo_obj + assert np.allclose(received_z, z) + assert m_lim == 24.5 + assert m_bright == -24.0 + assert m_faint == -18.0 + assert n_m == 128 + assert h is None + assert k_correction is None + assert e_correction is None + return np.array([1.0, 2.0]) + + monkeypatch.setattr(lf_api, "observed_number_density", fake_observed_number_density) + monkeypatch.setattr(lf_api, "missing_number_density", fake_missing_number_density) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + observed = lf.observed_number_density( + cosmo_obj, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=128, + ) + missing = lf.missing_number_density( + cosmo_obj, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + n_m=128, + ) + + assert np.allclose(observed, [4.0, 6.0]) + assert np.allclose(missing, [1.0, 2.0]) + + +def test_catalog_and_out_of_catalog_fractions_sum_to_one(monkeypatch): + """Tests that catalog and out-of-catalog wrappers preserve fraction semantics.""" + cosmo_obj = object() + z = np.array([0.1, 0.5]) + + monkeypatch.setattr( + lf_api, + "catalog_completeness_fraction", + lambda *args, **kwargs: np.array([0.75, 0.25]), + ) + monkeypatch.setattr( + lf_api, + "out_of_catalog_fraction", + lambda *args, **kwargs: np.array([0.25, 0.75]), + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.2, + ) + + completeness = lf.catalog_completeness( + cosmo_obj, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + ) + missing = lf.out_of_catalog_fraction( + cosmo_obj, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-18.0, + ) + + assert np.allclose(completeness + missing, 1.0) From d8c34f571f4a315d44f5915164eb4423c8d5c334 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:26:01 -0400 Subject: [PATCH 06/10] Add fake magnitude-limited demo catalog --- .../make_fake_magnitude_limited_catalog.py | 40 ++++ .../fake_magnitude_limited_catalog.csv | 201 ++++++++++++++++++ 2 files changed, 241 insertions(+) create mode 100644 scripts/make_fake_magnitude_limited_catalog.py create mode 100644 src/lfkit/data/demo_catalogs/fake_magnitude_limited_catalog.csv diff --git a/scripts/make_fake_magnitude_limited_catalog.py b/scripts/make_fake_magnitude_limited_catalog.py new file mode 100644 index 0000000..0b3695a --- /dev/null +++ b/scripts/make_fake_magnitude_limited_catalog.py @@ -0,0 +1,40 @@ +"""Create a small fake magnitude-limited galaxy catalog for examples.""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pandas as pd + + +def main() -> None: + """Create a deterministic fake catalog with redshifts and apparent magnitudes.""" + rng = np.random.default_rng(42) + + n_gal = 200 + z = np.sort(rng.uniform(0.01, 1.2, n_gal)) + + # Fake magnitude-redshift trend plus scatter. + m_app = 15.0 + 6.0 * np.log10(1.0 + 8.0 * z) + rng.normal(0.0, 0.35, n_gal) + + catalog = pd.DataFrame( + { + "galaxy_id": [f"fake-{i:05d}" for i in range(n_gal)], + "ra_deg": rng.uniform(0.0, 360.0, n_gal), + "dec_deg": rng.uniform(-30.0, 30.0, n_gal), + "z": z, + "m_app": m_app, + "band": "r", + } + ) + + repo_root = Path(__file__).resolve().parents[1] + output = repo_root / "src/lfkit/data/demo_catalogs/fake_magnitude_limited_catalog.csv" + + output.parent.mkdir(parents=True, exist_ok=True) + catalog.to_csv(output, index=False) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/lfkit/data/demo_catalogs/fake_magnitude_limited_catalog.csv b/src/lfkit/data/demo_catalogs/fake_magnitude_limited_catalog.csv new file mode 100644 index 0000000..308f9b1 --- /dev/null +++ b/src/lfkit/data/demo_catalogs/fake_magnitude_limited_catalog.csv @@ -0,0 +1,201 @@ +galaxy_id,ra_deg,dec_deg,z,m_app,band +fake-00000,32.55866171428779,28.188907002334165,0.01876110100369656,15.589552641209139,r +fake-00001,322.9320236027851,29.074684885897476,0.035718375057593205,15.51682132593073,r +fake-00002,10.439821153608992,-12.706305270725684,0.03702736702929397,15.674277516818826,r +fake-00003,86.69810089481294,14.02520984627497,0.037136606525340024,15.62062138214073,r +fake-00004,51.4878750599886,14.999012276116117,0.046673223135847895,15.944906380178335,r +fake-00005,279.6364586491987,-9.210428321217073,0.054520941849653934,16.435849281958593,r +fake-00006,71.35352179664324,-22.567813508146738,0.059520866895652044,16.046517557140955,r +fake-00007,327.8297617671187,-27.543182382856813,0.06212648128680224,16.27673000694812,r +fake-00008,236.25685415789644,16.640587651750714,0.06344363707687149,15.35206855336005,r +fake-00009,13.018575778090637,-0.6180154909856461,0.07938026260998855,16.264122198346048,r +fake-00010,1.9547402863321661,29.13241016969557,0.08594253476396863,16.068394289448463,r +fake-00011,18.596850125183614,-2.1015926283808173,0.1084486999077019,16.201104108744886,r +fake-00012,218.13306395318196,28.675018744430414,0.110488742156468,16.343007493241185,r +fake-00013,288.5334519093895,-5.305439868356407,0.11430340398616018,16.575289975543914,r +fake-00014,85.87901542585985,17.620929033598355,0.11715695432301368,17.04368748620109,r +fake-00015,305.78718345230624,-24.910843661290077,0.12207104398630293,16.311233572685676,r +fake-00016,20.60349845517374,3.3277026082658345,0.12241533499320316,16.789820628564502,r +fake-00017,288.34698731100343,18.123587219143452,0.12470524496266415,16.63364660809952,r +fake-00018,334.00635486202094,25.48210003319332,0.13304953158983646,16.774072789390623,r +fake-00019,277.959023677577,19.354985437045897,0.1382117256875169,17.291336556051437,r +fake-00020,251.3234822323094,-27.781756367929546,0.1392051319511778,17.1385278197845,r +fake-00021,301.6728787159379,-7.637859521595949,0.14629078750780838,17.48723953700004,r +fake-00022,14.454467836943913,-27.078091661168923,0.15042702352341322,17.004502227110738,r +fake-00023,72.64155985885994,-23.44306253312652,0.15608193956966837,16.867957657878367,r +fake-00024,44.97252455003151,10.518337720049047,0.16245522288389957,17.091610559540957,r +fake-00025,181.63115646947514,12.79549177870189,0.16460659134921124,17.274264451105562,r +fake-00026,268.267726198384,16.423240965097577,0.17630545549507867,17.35438256839901,r +fake-00027,226.80426404750352,21.92739288333977,0.17635842777191235,16.913504145234853,r +fake-00028,306.4071959142966,14.365888108590593,0.17689641590834682,17.329359088374137,r +fake-00029,55.87667727879272,18.052295526397923,0.1819837847505958,17.4211124542854,r +fake-00030,264.46359309665814,-27.062177410603244,0.1912514022286704,18.299850593940214,r +fake-00031,69.49493671280547,-15.927890973522086,0.1924083367831479,18.085145512216048,r +fake-00032,97.47315047440809,7.313866402514746,0.1936044955603819,17.13941614609369,r +fake-00033,255.5656910212694,21.48751827748939,0.2019134170499861,17.40455221593107,r +fake-00034,352.87372256352546,-29.72999250362045,0.2058833029226016,17.024384335028724,r +fake-00035,220.15569815174504,0.8777605269809885,0.20869777469016765,17.351912330133686,r +fake-00036,19.620113396866415,10.637243961863831,0.21300563255889446,17.702545297990156,r +fake-00037,221.87122917391181,-28.223562652283206,0.22035961169768573,18.070217973380316,r +fake-00038,15.246198568860155,-5.918666546796725,0.22332928670770952,17.41530052145709,r +fake-00039,318.2924560920542,23.73809288302181,0.23547091731029995,17.530810058219245,r +fake-00040,255.44818263622346,10.296768531351523,0.2416200623438414,17.05228612312066,r +fake-00041,62.32602471210116,-15.740498198364161,0.24375725220034214,17.76205027075003,r +fake-00042,33.019562096711304,21.16686777791476,0.2478907609453789,17.476184616099825,r +fake-00043,66.07196239691483,-9.118114595429141,0.25081240410632405,17.683062784614,r +fake-00044,352.8097846395093,21.200680268643566,0.26535576065523975,17.660402206553066,r +fake-00045,165.08183128505632,-12.063380932865668,0.2800221253705521,18.030420517827277,r +fake-00046,282.2691414029069,5.419215052988939,0.2804140789238845,17.450727943858567,r +fake-00047,229.10700315879967,-6.183595935129955,0.28395464916490815,17.575124984825862,r +fake-00048,206.06873398510945,-13.510496954695979,0.2883879881797555,18.861923530206145,r +fake-00049,52.24689167941789,23.193453730312154,0.29675498425676483,17.718303360938197,r +fake-00050,340.5688032809142,-18.74437908363052,0.31459845777544726,17.893001930316252,r +fake-00051,108.48334770529846,-24.911304525854213,0.322328013948209,18.965216237344166,r +fake-00052,208.08619769071518,-9.484383687260378,0.3263852541368283,19.362597635777924,r +fake-00053,251.9193400520638,13.058348858494053,0.32770103907935333,17.94336054982717,r +fake-00054,233.7239358883112,18.445896358576903,0.33396745859573923,18.260344154504644,r +fake-00055,338.61398749038665,29.92460220400589,0.34466820276841986,18.56883131758148,r +fake-00056,53.438036371220875,-12.218276579330105,0.34484683150617584,19.055321912258847,r +fake-00057,183.00698584089915,-5.52348293212307,0.35311044367699046,18.150309326757043,r +fake-00058,145.45238066667574,-21.790723318222422,0.35619222738674117,18.426604297725262,r +fake-00059,170.70074261374253,4.492315777382014,0.3593765717423535,18.801706848082123,r +fake-00060,42.918309431188064,29.854802235775033,0.36879938608597307,18.732011321983656,r +fake-00061,48.27405955419552,12.052860607239511,0.3717006166945085,18.463453492267554,r +fake-00062,100.10719648125266,5.712770424802304,0.37408838273927647,18.560766649939897,r +fake-00063,109.69365735389927,-6.45785442469236,0.3792851707511128,18.153383572686614,r +fake-00064,154.045156925286,24.91792561713329,0.38171630324462885,18.563768377390694,r +fake-00065,219.95551693318617,-0.1850042321790335,0.38595557167864364,18.575603135930262,r +fake-00066,228.4664822705747,-21.93798520121631,0.38739528245903126,18.757430118990438,r +fake-00067,148.2519227708053,-8.077292208051992,0.3977321761844008,18.533850795904637,r +fake-00068,147.16191939433097,-25.969999869133986,0.4025348426760145,18.91708530030669,r +fake-00069,78.34626961090632,-17.881257738673106,0.4045671068376371,19.116516469015416,r +fake-00070,211.79024942724297,-28.93987314000527,0.4227750673964561,18.904553752406574,r +fake-00071,114.13472807420564,-2.803205009990195,0.43188590207454336,19.016252006396957,r +fake-00072,12.981540339850932,8.07241589710783,0.4405570065585346,18.95200261639835,r +fake-00073,150.62401589646973,-9.402452257265534,0.45084705018149396,18.980412338375444,r +fake-00074,170.6877630324613,-4.777093709683373,0.45124964883677166,18.72965865265899,r +fake-00075,81.21343255351596,27.552556370674555,0.4513975177963021,19.093645524914116,r +fake-00076,206.08485592396164,15.117787294494931,0.454652697359715,18.963496779432095,r +fake-00077,203.67788415549154,2.451398216697074,0.4634152590548141,19.76925238923063,r +fake-00078,252.7207851997337,-12.927547493290486,0.4710992710459076,19.621125682273135,r +fake-00079,233.2614535924356,23.819807926379713,0.4936003651136839,19.30204724733981,r +fake-00080,234.87590034757136,-15.894172987160431,0.496149086032317,18.910645001746857,r +fake-00081,113.83709465665552,-10.479436168999376,0.5048108085830254,18.824456564163846,r +fake-00082,283.4755999149434,24.54388886537691,0.5241907173323734,19.709673763306824,r +fake-00083,197.6919780958846,1.772523317104342,0.5277655014036151,19.39904531432772,r +fake-00084,155.31055026476045,14.539076965078602,0.5296936931865112,19.48281905684554,r +fake-00085,225.36449313878055,5.444687649825967,0.5302107834580736,18.706221045085297,r +fake-00086,129.83664040209155,9.206352540569709,0.5322653433049422,19.649587789307187,r +fake-00087,184.58612805843296,-12.037002509110465,0.5323046377630055,19.484187219657237,r +fake-00088,265.2140478535993,-15.517676302400226,0.537662896604524,18.95764510961073,r +fake-00089,319.10503922053476,-10.650459169880085,0.5409259679330966,19.19406194314749,r +fake-00090,331.580591026269,-20.67350615499465,0.5447401875511085,19.466279100014585,r +fake-00091,181.3078530659157,22.458861917448445,0.5459592660957249,19.397080231211426,r +fake-00092,187.29904130660486,-13.005184004993156,0.5523737849019972,19.30125031672683,r +fake-00093,287.95334787498115,3.6893636648630768,0.5561097728906246,19.38162169888916,r +fake-00094,113.20224902791246,17.518465508482443,0.5596339039991509,19.343098252695007,r +fake-00095,301.45765044422467,17.029446562542812,0.5653979944351707,19.506533927844668,r +fake-00096,177.89099274627245,-3.6968244849484266,0.5687714154182114,19.980860525646122,r +fake-00097,41.70842076069662,-1.4245614847964605,0.5706044853103277,18.574383709845907,r +fake-00098,25.94129294625973,29.68210493610149,0.573452387257775,19.40046316582396,r +fake-00099,303.1175559742067,10.475848618741443,0.5760888622088611,19.554957676128566,r +fake-00100,20.004450086476044,18.878306587072963,0.5839410899620513,19.625798710111734,r +fake-00101,101.02011700811613,24.153238153153914,0.5891234136974108,19.411009452643384,r +fake-00102,120.2868145805565,17.255390173760354,0.5939413232818799,18.943847978227605,r +fake-00103,62.27800026897376,-18.889239302785946,0.5978490749009775,19.687624118546225,r +fake-00104,113.00161312871612,3.730244024254332,0.6062432823730026,20.207483833467066,r +fake-00105,267.3693240210308,-23.88635056112731,0.6229179837857403,19.12481228633064,r +fake-00106,5.285823681625477,9.175327593839903,0.6274414488598744,19.979715992795303,r +fake-00107,297.7824328274471,27.32096569126866,0.62960015894113,19.569857447900205,r +fake-00108,308.35728848684096,0.7639238580066787,0.668759486783015,19.795209974267486,r +fake-00109,134.01416634895733,-4.021650427933969,0.6693030108114687,19.449942695699125,r +fake-00110,55.30064366182853,-27.84943491491021,0.6699558965488434,19.703538471929846,r +fake-00111,216.3025468602291,27.58647133975322,0.6702744385758652,20.276657672539315,r +fake-00112,43.08212011136196,-23.819838952190473,0.6728682612861212,20.03405639332633,r +fake-00113,131.3709699881709,-27.535254401709185,0.6754565212870421,20.444875466751686,r +fake-00114,345.03450512435006,-15.236003309703692,0.6826309667126147,20.2739118586505,r +fake-00115,358.36721013010936,-26.068173964897056,0.6868020231839436,20.028921035600074,r +fake-00116,277.9577608876508,-2.6929294552942595,0.6903563823275297,20.497001856451543,r +fake-00117,111.94614357147583,0.965268075426744,0.6992712870408562,20.06860855794264,r +fake-00118,247.5594176997933,-11.245803058284965,0.70146275624347,20.211675507519416,r +fake-00119,253.94629157424458,-26.942371045506835,0.7050765830061553,19.829449062056497,r +fake-00120,139.62301025891392,-23.3039779050464,0.7092970408914773,19.96975517681343,r +fake-00121,230.71990845101737,-6.929737320604314,0.7131894348628859,19.714494678047178,r +fake-00122,3.8619521912434474,-26.36829295040783,0.7232655219098038,20.336085190788438,r +fake-00123,75.26075709680693,11.891383536594805,0.7564975204359073,19.67746022392296,r +fake-00124,189.031789068214,-17.57876750851128,0.7600362858114772,20.37412900262645,r +fake-00125,58.95046953298652,-11.87570587298169,0.7616806349552572,20.038415393463104,r +fake-00126,59.72647243555211,-6.3524533693141505,0.765314800800703,20.525737964158132,r +fake-00127,301.06954459808566,-5.003422015851127,0.7761994928959908,20.410269700056478,r +fake-00128,356.087880966098,-29.90035753951197,0.7833174558502357,20.805192039809295,r +fake-00129,200.14899409024142,-23.275777740078148,0.7947129352833631,20.456224831770975,r +fake-00130,302.0651031174168,21.765873824452683,0.7973769031696356,19.65777634915842,r +fake-00131,356.5157991574581,-29.92601624928551,0.7976806525250051,20.185414795960295,r +fake-00132,50.97451988475889,0.49004134877250527,0.8005331129895912,19.806689171805726,r +fake-00133,161.36842077156638,-0.6298796095675141,0.8011725193445182,20.03729332197187,r +fake-00134,141.32617770309116,-10.01435033936173,0.8053995245306612,20.75948675227579,r +fake-00135,28.817742112626608,-4.120378791213022,0.8070786536721872,20.458392127408118,r +fake-00136,271.91886220610803,16.834868251595058,0.8099657107173759,19.998688374065708,r +fake-00137,156.16044983613347,20.472222372897384,0.8221696497302208,19.922305484409776,r +fake-00138,168.95769630854446,-14.379087395084998,0.8228282543585209,20.29039155860299,r +fake-00139,54.242270572766465,-10.650573541936211,0.8386212463425058,19.89618207134707,r +fake-00140,65.13359480311253,-15.451028818300037,0.839867954580643,20.090448121376316,r +fake-00141,326.5573039762374,-1.208195966936092,0.8426514892347919,20.442057053022886,r +fake-00142,16.073672030740397,10.995501459528775,0.8433154713826764,20.73900036001062,r +fake-00143,83.82682257973426,-16.304827472194766,0.8491468005653388,20.56335288818212,r +fake-00144,105.1413589132302,-10.155855476708,0.8630993250992047,19.58539220989541,r +fake-00145,176.4711152645711,25.823077055349028,0.8661609175816147,20.501933291477222,r +fake-00146,211.12026228532653,-27.085842609182514,0.8696076272097758,20.42966238452063,r +fake-00147,177.58439130623998,-2.3538236856260264,0.8751235910013903,20.563723551272044,r +fake-00148,30.28152044047753,12.69348249580667,0.8962669655303024,21.03904901575279,r +fake-00149,87.72028346005688,-20.972759306551556,0.8986796375053409,19.757390912767498,r +fake-00150,303.6918185145148,-27.15755879967109,0.9008435431537696,20.278140498746918,r +fake-00151,229.53193217224526,-21.707677701688155,0.9121244107015948,20.720342272918522,r +fake-00152,233.69365804833876,25.129391611840973,0.9126382405039,19.96125769375506,r +fake-00153,241.27317192648778,-29.444413090732596,0.91575624536852,21.03921624285226,r +fake-00154,274.64508685423044,-18.70068151131889,0.9203486403250704,20.663031688556618,r +fake-00155,20.919053419721042,-28.12298922590133,0.9246905153820119,20.841211955565996,r +fake-00156,131.97901857574047,-23.362231946267805,0.9286947942120757,20.354998646943944,r +fake-00157,194.22987669420434,7.208957197697776,0.9310076977815963,20.845359717466856,r +fake-00158,121.84433398703737,-15.501665475654853,0.9362763615177766,20.94747591267198,r +fake-00159,304.0123943833871,4.152723025840608,0.9367135757192956,20.656091452285338,r +fake-00160,173.7261030914258,5.411722965967569,0.9390675469161417,20.66239542125845,r +fake-00161,276.7059322053765,20.96611998538674,0.9414549204278405,20.68081506969727,r +fake-00162,306.72558607722084,-29.715516425287564,0.942838868836632,20.28740342884963,r +fake-00163,181.7249338452593,21.202167560542005,0.945416523279575,20.544221710141844,r +fake-00164,327.43880779368595,7.14714518561844,0.9464400092275447,20.544964192998204,r +fake-00165,211.364618593313,-20.235390833702546,0.9676695854211942,20.783660902135146,r +fake-00166,306.09874758013984,16.376243288117042,0.9718192517565191,21.009288796795577,r +fake-00167,122.61268641294873,21.329492507435145,0.9786842577525812,20.30512138183417,r +fake-00168,179.57410507133469,-14.744190974750563,0.9890863197922877,20.65629980990538,r +fake-00169,191.30797476164798,25.127620854160853,0.9932873897912857,21.22836992962909,r +fake-00170,37.79269772258541,-2.732554182080531,0.9948810946711726,20.45331552072999,r +fake-00171,143.47890241255462,6.205467259999288,0.9974497943275716,20.431753964103905,r +fake-00172,330.2415621186034,29.078972017058618,0.9990240798286383,20.794000780110146,r +fake-00173,227.09960653130727,-8.307835281648039,1.0003891636602893,21.021891112796126,r +fake-00174,63.90236967730829,18.78956074920677,1.0008870533088265,20.731507992384806,r +fake-00175,121.98802881945778,-10.896052930292422,1.0255496571891176,21.249108332749024,r +fake-00176,68.97708347857194,17.95281133448946,1.0293709980699268,21.09249025899659,r +fake-00177,8.936327447486114,6.04403513175855,1.031731524694545,21.09257237544193,r +fake-00178,333.8857650614496,-17.018658286475915,1.0507288285965826,21.034323700181275,r +fake-00179,161.35463817247373,-5.158434372645817,1.0700430641454446,21.697519840955565,r +fake-00180,110.71262607148302,-10.9418619381863,1.0710959468433476,20.813329443110536,r +fake-00181,215.45178896078517,-25.313496415178168,1.0728141343734152,20.18764372081075,r +fake-00182,2.633204265485629,-28.20991379401815,1.0732015012404572,21.451208104806124,r +fake-00183,100.08795837629297,-9.211280675569144,1.0764313184893224,20.736538813885574,r +fake-00184,253.09204765160132,-28.857951130950806,1.0840968518871716,20.95106461988433,r +fake-00185,228.15711830493467,-20.070734545093497,1.0912110219420523,21.38693627104578,r +fake-00186,353.45014111531276,13.510999251052127,1.1000540957137535,20.38668068230072,r +fake-00187,223.32877547689586,12.485474141015274,1.1108929408683885,20.532348506949454,r +fake-00188,171.9021145059838,14.31482318744296,1.1128503367298361,20.414101478810068,r +fake-00189,274.11572274624274,-10.970809732201275,1.1238918483582663,20.71974120603635,r +fake-00190,325.1980339070591,23.40116742706173,1.1248076491351973,21.153471804039427,r +fake-00191,259.45054081860656,5.6298328817068555,1.1260134814202296,21.185577606031263,r +fake-00192,346.75604047777665,-22.43841949520241,1.1439502801329964,21.13590352042009,r +fake-00193,281.52186146574377,-21.37688031266795,1.1506854637573198,20.55853350360133,r +fake-00194,312.04851777429354,11.58572561578209,1.1546582208139222,20.25256766529105,r +fake-00195,41.07746562610939,-19.623363813475315,1.1613365815967098,21.093691400052204,r +fake-00196,263.66886107776355,0.4375620949225989,1.1633198890437915,20.91356057899252,r +fake-00197,158.43193182985272,29.50646855572264,1.165130649029935,21.243127098172096,r +fake-00198,199.1173688680591,-29.75954371599963,1.1709905984477396,21.339834797294525,r +fake-00199,235.47686741364842,-29.005211886146583,1.190926921226446,21.18231505402518,r From 1230c195bf07f9c96cc4f13177f67d6083eb87f4 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:26:05 -0400 Subject: [PATCH 07/10] Add catalog completeness API docs --- .../lfkit.photometry.catalog_completeness.rst | 7 ++++ docs/api/lfkit.photometry.rst | 1 + docs/api/lfkit.utils.rst | 1 + docs/api/lfkit.utils.types.rst | 7 ++++ docs/conf.py | 16 +++++++++ pyproject.toml | 33 ++++++++++++++----- 6 files changed, 57 insertions(+), 8 deletions(-) create mode 100644 docs/api/lfkit.photometry.catalog_completeness.rst create mode 100644 docs/api/lfkit.utils.types.rst diff --git a/docs/api/lfkit.photometry.catalog_completeness.rst b/docs/api/lfkit.photometry.catalog_completeness.rst new file mode 100644 index 0000000..791e459 --- /dev/null +++ b/docs/api/lfkit.photometry.catalog_completeness.rst @@ -0,0 +1,7 @@ +lfkit.photometry.catalog\_completeness module +============================================= + +.. automodule:: lfkit.photometry.catalog_completeness + :members: + :show-inheritance: + :undoc-members: diff --git a/docs/api/lfkit.photometry.rst b/docs/api/lfkit.photometry.rst index 1587a8a..890f6a8 100644 --- a/docs/api/lfkit.photometry.rst +++ b/docs/api/lfkit.photometry.rst @@ -7,6 +7,7 @@ Submodules .. toctree:: :maxdepth: 2 + lfkit.photometry.catalog_completeness lfkit.photometry.lf_parameter_models lfkit.photometry.luminosities lfkit.photometry.luminosity_function diff --git a/docs/api/lfkit.utils.rst b/docs/api/lfkit.utils.rst index 12d22ef..80b4e19 100644 --- a/docs/api/lfkit.utils.rst +++ b/docs/api/lfkit.utils.rst @@ -10,6 +10,7 @@ Submodules lfkit.utils.download_poggianti97_data lfkit.utils.interpolation lfkit.utils.io + lfkit.utils.types lfkit.utils.units lfkit.utils.validators diff --git a/docs/api/lfkit.utils.types.rst b/docs/api/lfkit.utils.types.rst new file mode 100644 index 0000000..6fbcd98 --- /dev/null +++ b/docs/api/lfkit.utils.types.rst @@ -0,0 +1,7 @@ +lfkit.utils.types module +======================== + +.. automodule:: lfkit.utils.types + :members: + :show-inheritance: + :undoc-members: diff --git a/docs/conf.py b/docs/conf.py index e462c9d..23e2593 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -4,6 +4,7 @@ import os import sys +import warnings import matplotlib @@ -14,6 +15,12 @@ # ----------------------------------------------------------------------------- sys.path.insert(0, os.path.abspath("../src")) +warnings.filterwarnings( + "ignore", + category=SyntaxWarning, + module=r"colorspacious\.comparison", +) + # ----------------------------------------------------------------------------- # Project information # ----------------------------------------------------------------------------- @@ -38,6 +45,7 @@ "sphinx_design", "sphinx_multiversion", "sphinx.ext.mathjax", + "sphinx_copybutton", ] #templates_path = ["_templates"] # if uncomment this removes the sidebar logo @@ -106,3 +114,11 @@ # ----------------------------------------------------------------------------- smv_tag_whitelist = r"^v\d+\.\d+\.\d+$" smv_branch_whitelist = "main" + + +# ----------------------------------------------------------------------------- +# Copybutton configuration +# ----------------------------------------------------------------------------- +copybutton_prompt_text = r">>> |\.\.\. " +copybutton_prompt_is_regexp = True +copybutton_copy_empty_lines = False diff --git a/pyproject.toml b/pyproject.toml index 888570b..af58c2e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ viz = [ docs = [ "sphinx<9", "sphinx-design", + "sphinx-copybutton", "furo", "sphinx-multiversion", #"sphinxawesome-theme", @@ -46,9 +47,7 @@ include-package-data = true [tool.setuptools.package-data] "lfkit" = [ - "data/kcorrect/grids/*.npz", - "data/kcorrect/responses/*.dat", - "data/poggianti1997/*.csv", + "data/**/*", ] [tool.setuptools.packages.find] @@ -138,12 +137,18 @@ commands = [ ], ] - [tool.tox.env.do] -description = "Run doctests + html build (uses current conda env; does not install deps)" +description = "Run doctests + html build and open local docs" +recreate = true extras = ["docs", "viz"] -allowlist_externals = ["open"] +allowlist_externals = ["open", "rm"] commands = [ + [ + "rm", + "-rf", + "{tox_root}/docs/_build/doctest", + "{tox_root}/docs/_build/html", + ], [ "sphinx-apidoc", "-d", "2", @@ -153,8 +158,20 @@ commands = [ "--no-toc", "-f", ], - ["sphinx-build", "-b", "doctest", "{tox_root}/docs", "{tox_root}/docs/_build/doctest", "--fail-on-warning"], - ["sphinx-build", "-b", "html", "{tox_root}/docs", "{tox_root}/docs/_build/html", "--fail-on-warning"], + [ + "sphinx-build", + "-b", "doctest", + "{tox_root}/docs", + "{tox_root}/docs/_build/doctest", + "--fail-on-warning", + ], + [ + "sphinx-build", + "-b", "html", + "{tox_root}/docs", + "{tox_root}/docs/_build/html", + "--fail-on-warning", + ], ["open", "{tox_root}/docs/_build/html/index.html"], ] From ce7d7e36ac37b003b95371e2f9bc5ca1ed4ad902 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:26:10 -0400 Subject: [PATCH 08/10] Add luminosity-function completeness examples --- docs/about/index.rst | 2 +- docs/about/lf_overview.rst | 482 ++++++++- .../catalog_completeness_examples.rst | 868 ++++++++++++++++ docs/examples/index.rst | 32 +- .../examples/luminosity_function_examples.rst | 959 ++++++++++++++++++ docs/index.rst | 66 +- 6 files changed, 2393 insertions(+), 16 deletions(-) create mode 100644 docs/examples/catalog_completeness_examples.rst create mode 100644 docs/examples/luminosity_function_examples.rst diff --git a/docs/about/index.rst b/docs/about/index.rst index f7f6765..a2ace1e 100644 --- a/docs/about/index.rst +++ b/docs/about/index.rst @@ -37,7 +37,7 @@ Core components :link-type: doc :shadow: md - **Luminosity Function Tools** + **Photometric Tools** ^^^ Utilities for constructing and analyzing galaxy luminosity diff --git a/docs/about/lf_overview.rst b/docs/about/lf_overview.rst index 3154893..97d3d84 100644 --- a/docs/about/lf_overview.rst +++ b/docs/about/lf_overview.rst @@ -2,5 +2,483 @@ :alt: LFKit logo black :width: 50px -|lfkitlogo| Luminosity function -=============================== \ No newline at end of file +|lfkitlogo| Luminosity Functions +================================ + +LFKit provides a public luminosity-function interface through +:class:`~lfkit.api.lumfunc.LuminosityFunction`. + +For most users, the recommended import is: + +.. code-block:: python + + from lfkit import LuminosityFunction + +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 +be loaded elsewhere and passed in as scalars, arrays, or correction objects. + + +Magnitude-space luminosity functions +------------------------------------ + +LFKit luminosity functions are evaluated in rest-frame absolute magnitude +:math:`M`. + +The standard magnitude-space Schechter luminosity function is + +.. math:: + + \phi(M) = + 0.4 \ln(10) \, \phi_\star \, + x^{\alpha + 1} \exp(-x), + +where + +.. math:: + + x = 10^{-0.4(M - M_\star)}. + +Here: + +- :math:`\phi_\star` is the normalization, +- :math:`M_\star` is the characteristic magnitude, +- :math:`\alpha` is the faint-end slope. + +By convention, more negative magnitudes are brighter. + + +Standard Schechter model +------------------------ + +Use :meth:`~lfkit.api.lumfunc.LuminosityFunction.schechter` to construct a +Schechter luminosity function with fixed parameters. + +.. code-block:: python + + from lfkit import LuminosityFunction + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.1, + ) + + phi = lf.phi(absolute_mag, z) + +The returned ``phi`` values are number densities per magnitude evaluated at the +input absolute magnitudes. + + +Evolving Schechter model +------------------------ + +Use :meth:`~lfkit.api.lumfunc.LuminosityFunction.evolving_schechter` to build a +Schechter luminosity function with redshift-dependent parameters. + +The evolving model evaluates + +.. math:: + + \phi_\star(z), \quad M_\star(z), \quad \alpha(z), + +and then evaluates the Schechter function at each redshift. + +.. code-block:: python + + from lfkit import LuminosityFunction + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 1.0}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.2, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + phi = lf.phi(absolute_mag, z) + +You can also evaluate the evolving parameters directly: + +.. code-block:: python + + phi_star, m_star, alpha = lf.parameters(z) + + +Double Schechter model +---------------------- + +Use :meth:`~lfkit.api.lumfunc.LuminosityFunction.double_schechter` to build a +double-power-law Schechter-style luminosity function. + +.. code-block:: python + + from lfkit import LuminosityFunction + + lf = LuminosityFunction.double_schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.0, + beta=-1.5, + m_transition=-18.0, + ) + + phi = lf.phi(absolute_mag, z) + +This model is useful when an additional slope or transition is needed beyond +the standard Schechter form. + + +Built-in parameter evolution models +----------------------------------- + +Redshift-dependent luminosity-function parameters are handled by parameter +evolution models. + +The built-in options are: + +.. list-table:: + :header-rows: 1 + :widths: 25 25 50 + + * - Parameter + - Model name + - Form + * - ``phi_star`` + - ``constant`` + - :math:`\phi_\star(z) = \phi_\star` + * - ``phi_star`` + - ``linear_p`` + - :math:`\phi_\star(z) = \phi_{0,\star} 10^{0.4 p z}` + * - ``M_star`` + - ``constant`` + - :math:`M_\star(z) = M_\star` + * - ``M_star`` + - ``linear_q`` + - :math:`M_\star(z) = M_{0,\star} - q(z - z_{\rm ref})` + * - ``alpha`` + - ``constant`` + - :math:`\alpha(z) = \alpha` + * - ``alpha`` + - ``linear`` + - :math:`\alpha(z) = \alpha_0 + \alpha_1(z - z_{\rm ref})` + +To inspect the available models: + +.. code-block:: python + + models = LuminosityFunction.available_parameter_models() + +Custom parameter-evolution models can be registered through: + +.. code-block:: python + + LuminosityFunction.register_phi_star_model(name, model) + LuminosityFunction.register_m_star_model(name, model) + LuminosityFunction.register_alpha_model(name, model) + +Each registered model should accept redshift values as its first argument and +return NumPy-compatible parameter values. + + +Evaluating 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 + +.. math:: + + M = m - \mu(z) - K(z) + E(z), + +and then evaluates :math:`\phi(M, z)`. + +.. code-block:: python + + from lfkit import Corrections, LuminosityFunction + + corr = Corrections.poggianti( + band="r", + gal_type="E", + ) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.1, + ) + + phi = lf.phi_from_m( + cosmo, + z, + apparent_mag, + corrections=corr, + ) + +The ``corrections`` argument is optional. If omitted, LFKit evaluates the +conversion without k- or e-corrections. + + +Magnitude conversion helpers +---------------------------- + +The same object can convert between apparent and absolute magnitudes using the +LFKit magnitude convention. + +.. code-block:: python + + absolute_mag = lf.absolute_magnitude( + cosmo, + z, + apparent_mag, + corrections=corr, + ) + + apparent_mag = lf.apparent_magnitude( + cosmo, + z, + absolute_mag, + corrections=corr, + ) + +The absolute-magnitude limit corresponding to an apparent-magnitude catalog cut +can be computed with: + +.. code-block:: python + + m_limit = lf.absolute_magnitude_limit( + cosmo, + z, + m_lim=24.5, + corrections=corr, + ) + + +Integrated number density +------------------------- + +Use :meth:`~lfkit.api.lumfunc.LuminosityFunction.integrated_number_density` to +integrate the luminosity function over an absolute-magnitude range. + +.. code-block:: python + + n_total = lf.integrated_number_density( + z, + m_bright=-24.0, + m_faint=-14.0, + ) + +This computes the number density inside the requested absolute-magnitude +interval. + + +Magnitude-limited catalog completeness +-------------------------------------- + +A luminosity function can be split into observed and missing components for a +magnitude-limited catalog. + +The core idea is: + +1. convert the apparent catalog limit :math:`m_{\rm lim}` into an + absolute-magnitude limit :math:`M_{\rm lim}(z)`, +2. integrate the luminosity function over a finite absolute-magnitude range, +3. split the population into observed and missing pieces. + +The limiting absolute magnitude is + +.. math:: + + M_{\rm lim}(z) = m_{\rm lim} - \mu(z) - K(z) + E(z). + + +Observed number density +^^^^^^^^^^^^^^^^^^^^^^^ + +The observed, or in-catalog, number density is + +.. math:: + + n_{\rm obs}(z) = + \int_{M_{\rm bright}}^{\min[M_{\rm lim}(z), M_{\rm faint}]} + \phi(M, z) \, dM. + +Use: + +.. code-block:: python + + n_obs = lf.observed_number_density( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-14.0, + corrections=corr, + ) + + +Missing number density +^^^^^^^^^^^^^^^^^^^^^^ + +The missing, or out-of-catalog, number density is + +.. math:: + + n_{\rm miss}(z) = + \int_{\max[M_{\rm lim}(z), M_{\rm bright}]}^{M_{\rm faint}} + \phi(M, z) \, dM. + +Use: + +.. code-block:: python + + n_miss = lf.missing_number_density( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-14.0, + corrections=corr, + ) + + +Completeness fractions +^^^^^^^^^^^^^^^^^^^^^^ + +The catalog completeness fraction is + +.. math:: + + f_{\rm obs}(z) = + \frac{n_{\rm obs}(z)} + {n_{\rm obs}(z) + n_{\rm miss}(z)}. + +The out-of-catalog fraction is + +.. math:: + + f_{\rm miss}(z) = 1 - f_{\rm obs}(z). + +Use: + +.. code-block:: python + + f_obs = lf.catalog_completeness( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-14.0, + corrections=corr, + ) + + f_miss = lf.out_of_catalog_fraction( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-14.0, + corrections=corr, + ) + + +Avoiding double-counted evolution +--------------------------------- + +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)`. + +For this reason, be careful when using an evolving Schechter model with +non-zero ``linear_q`` evolution together with an explicit evolution correction. + +This is not always wrong, but the two definitions should be intentionally +separated. + + +Typical workflow +---------------- + +A typical luminosity-function workflow is: + +1. define a cosmology, +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. + +For example: + +.. code-block:: python + + from lfkit import Corrections, LuminosityFunction + + corr = Corrections.poggianti( + band="r", + gal_type="E", + ) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 1.0}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.2, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + phi = lf.phi_from_m( + cosmo, + z, + apparent_mag, + corrections=corr, + ) + + n_obs = lf.observed_number_density( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-14.0, + corrections=corr, + ) + + +Lower-level functions +--------------------- + +The high-level API is recommended for most scripts, examples, notebooks, and +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, +- debugging the magnitude convention, +- building new public API objects, +- integrating LFKit into specialized workflows. + + +What LFKit does not do here +--------------------------- + +LFKit does not read survey catalogs, apply angular masks, or model survey-area +incompleteness in this layer. + +For magnitude-limited catalog applications, LFKit models the part of the +selection function that comes from the apparent-magnitude limit. Other effects, +such as unobserved sky area, angular masks, spectroscopic targeting, blending, +or color cuts, should be handled by the calling analysis code. diff --git a/docs/examples/catalog_completeness_examples.rst b/docs/examples/catalog_completeness_examples.rst new file mode 100644 index 0000000..48e9e47 --- /dev/null +++ b/docs/examples/catalog_completeness_examples.rst @@ -0,0 +1,868 @@ +.. |lfkitlogo| image:: /_static/logos/lfkit_logo-icon.png + :alt: LFKit logo + :width: 50px + +|lfkitlogo| Catalog completeness examples +========================================= + +This page shows how to use :class:`lfkit.LuminosityFunction` to compute +observed and missing galaxy number densities for a magnitude-limited catalog. + +Catalog completeness connects an apparent magnitude limit to the part of the +luminosity function that is actually observable. At fixed survey depth, nearby +galaxies can be observed to fainter intrinsic magnitudes, while at higher +redshift only brighter galaxies remain above the catalog limit. + +All examples below are executable via ``.. plot::``. + +The examples intentionally use ``corrections=None``. This means that no +:math:`K`-correction or evolution correction is applied. Users can pass their +own correction model through the ``corrections`` argument, for example a +:class:`lfkit.Corrections` object or any compatible correction callable used by +the LFKit magnitude-conversion methods. + +The number-density units follow the normalization of the luminosity function. +For example, if :math:`\phi_*` is given in comoving :math:`{\rm Mpc}^{-3}`, then +the integrated number densities are also in :math:`{\rm Mpc}^{-3}`. Fractions +are dimensionless and are always defined relative to the chosen intrinsic +absolute magnitude range. + + +Setup +----- + +First define a cosmology and a luminosity function. + +This setup plot shows the absolute magnitude limit implied by a fixed apparent +magnitude cut. The curve answers a simple question: at each redshift, how bright +does a galaxy need to be in absolute magnitude in order to enter the catalog? + +The limit becomes brighter at higher redshift because distant galaxies must be +intrinsically more luminous to remain above the same apparent magnitude +threshold. The y-axis is inverted so that brighter absolute magnitudes appear +higher on the plot. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + m_limit = lf.absolute_magnitude_limit( + cosmo, + z, + m_lim=24.5, + corrections=None, + ) + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + ax.plot(z, m_limit, lw=3, color=colors_blue[1]) + ax.invert_yaxis() + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"Absolute-magnitude limit $M_{\rm lim}(z)$", fontsize=LABEL_SIZE) + ax.set_title(r"Catalog limit for $m_{\rm lim}=24.5$", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + plt.tight_layout() + + +Observed and missing number densities +------------------------------------- + +The observed number density is the part of the luminosity function brighter +than the catalog limit. The missing number density is the part fainter than the +catalog limit but still inside the requested absolute magnitude range. + +This plot separates the full luminosity function integral into the part that is +included in the catalog and the part that falls below the magnitude limit. At +low redshift, the catalog limit is faint enough that most of the requested +magnitude range can be observed. At higher redshift, the observed contribution +shrinks and the missing contribution becomes more important. + +This is useful when diagnosing whether a catalog is tracing most of the galaxy +population described by the luminosity function, or only the bright tail of it. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + c_mid = 0.5 * (np.array(colors_blue[1]) + np.array(colors_red[1])) + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + m_bright = -24.0 + m_faint = -14.0 + + n_total = lf.integrated_number_density( + z, + m_bright=m_bright, + m_faint=m_faint, + ) + + n_obs = lf.observed_number_density( + cosmo, + z, + m_lim=24.5, + m_bright=m_bright, + m_faint=m_faint, + corrections=None, + ) + + n_miss = lf.missing_number_density( + cosmo, + z, + m_lim=24.5, + m_bright=m_bright, + m_faint=m_faint, + corrections=None, + ) + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + ax.plot(z, n_total, lw=3, color=c_mid, label="Total LF integral") + ax.plot(z, n_obs, lw=3, color=colors_blue[1], label=r"Observed: $M < M_{\rm lim}(z)$") + ax.plot(z, n_miss, lw=3, color=colors_red[1], label=r"Missing: $M_{\rm lim}(z) < M < M_{\rm faint}$") + ax.set_yscale("log") + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"Comoving number density $n(z)$ [$\mathrm{Mpc}^{-3}$]", fontsize=LABEL_SIZE) + ax.set_title(r"LF number densities for $-24 \leq M \leq -14$", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Completeness fraction +--------------------- + +The catalog completeness fraction is the observed fraction of the integrated +luminosity function over the chosen intrinsic absolute magnitude range. + +This plot shows the same information as a fraction rather than as an absolute +number density. The observed and missing fractions add up to one, so the figure +directly shows how much of the requested luminosity function range is retained +by the catalog limit. + +This is often the most useful summary diagnostic. It tells you where the sample +is effectively complete, where it is partially incomplete, and where the +catalog contains only a small fraction of the intrinsic galaxy population. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + m_bright = -24.0 + m_faint = -14.0 + + f_obs = lf.catalog_completeness( + cosmo, + z, + m_lim=24.5, + m_bright=m_bright, + m_faint=m_faint, + corrections=None, + ) + + f_miss = lf.out_of_catalog_fraction( + cosmo, + z, + m_lim=24.5, + m_bright=m_bright, + m_faint=m_faint, + corrections=None, + ) + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + ax.plot(z, f_obs, lw=3, color=colors_blue[1], label="Observable LF fraction") + ax.plot(z, f_miss, lw=3, color=colors_red[1], label="Missing LF fraction") + ax.set_ylim(-0.05, 1.05) + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel("Fraction of chosen LF integral", fontsize=LABEL_SIZE) + ax.set_title(r"Magnitude-limited completeness over $-24 \leq M \leq -14$", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="center right") + plt.tight_layout() + + +Comparing magnitude limits +-------------------------- + +A fainter apparent magnitude limit observes a larger fraction of the luminosity +function. + +This plot compares catalog completeness for several survey depths. Shallower +limits lose completeness earlier with redshift, while deeper limits retain a +larger fraction of the luminosity function over a wider redshift range. + +This is a useful way to visualize the gain from survey depth. A small change in +the apparent magnitude limit can produce a large change in completeness, +especially where the catalog is close to the transition between complete and +incomplete. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 4, cmap_range=(0.66, 0.98)) + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + limits = [22.5, 23.5, 24.5, 25.5] + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + + for m_lim, color in zip(limits, colors_blue): + f_obs = lf.catalog_completeness( + cosmo, + z, + m_lim=m_lim, + m_bright=-24.0, + m_faint=-14.0, + corrections=None, + ) + ax.plot(z, f_obs, lw=3, color=color, label=rf"$m_{{\rm lim}}={m_lim}$") + + ax.set_ylim(-0.05, 1.05) + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel("Observable LF fraction", fontsize=LABEL_SIZE) + ax.set_title(r"Magnitude-limit completeness over $-24 \leq M \leq -14$", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Absolute-magnitude limits for different depths +---------------------------------------------- + +The apparent magnitude limit can also be shown directly as an absolute magnitude +limit for several survey depths. + +This plot shows how the catalog threshold moves in intrinsic luminosity as the +survey becomes shallower or deeper. Deeper catalogs reach fainter absolute +magnitudes at the same redshift, while shallower catalogs only include brighter +galaxies. + +This is a useful diagnostic before looking at completeness fractions, because +it shows the selection boundary itself. The completeness calculation then asks +how much of the luminosity function lies on the observable side of this +boundary. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 4, cmap_range=(0.66, 0.98)) + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + limits = [22.5, 23.5, 24.5, 25.5] + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + + for m_lim, color in zip(limits, colors_blue): + m_limit = lf.absolute_magnitude_limit( + cosmo, + z, + m_lim=m_lim, + corrections=None, + ) + ax.plot(z, m_limit, lw=3, color=color, label=rf"$m_{{\rm lim}}={m_lim}$") + + ax.invert_yaxis() + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"Absolute-magnitude limit $M_{\rm lim}(z)$", fontsize=LABEL_SIZE) + ax.set_title("Catalog selection boundary", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Dependence on the faint-end integration limit +--------------------------------------------- + +Completeness is always defined relative to a chosen intrinsic magnitude range. + +This plot compares completeness fractions when the luminosity function is +integrated down to different faint-end limits. Extending the intrinsic range to +fainter galaxies usually lowers the completeness fraction, because the catalog +must now account for a larger population of faint galaxies that may fall below +the apparent magnitude limit. + +This is an important sanity check. A catalog can appear highly complete if the +requested intrinsic range is bright, but less complete if the calculation asks +about a much fainter underlying population. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_red = cmr.take_cmap_colors("cmr.guppy", 4, cmap_range=(0.03, 0.30)) + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + faint_limits = [-17.0, -16.0, -15.0, -14.0] + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + + for m_faint, color in zip(faint_limits, colors_red): + f_obs = lf.catalog_completeness( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=m_faint, + corrections=None, + ) + ax.plot( + z, + f_obs, + lw=3, + color=color, + label=rf"$M_{{\rm faint}}={m_faint}$", + ) + + ax.set_ylim(-0.05, 1.05) + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel("Observable LF fraction", fontsize=LABEL_SIZE) + ax.set_title(r"Completeness for varying faint-end LF limits", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Completeness as a redshift-depth surface +---------------------------------------- + +Completeness can also be viewed as a two-dimensional function of redshift and +survey depth. + +This plot shows how much of the chosen intrinsic galaxy population is visible +for different combinations of redshift and apparent magnitude limit. The +horizontal axis changes redshift, while the vertical axis changes catalog depth: +larger :math:`m_{\rm lim}` values correspond to deeper catalogs. + +The colour scale gives the observable fraction of the luminosity function +integral over :math:`-24 \leq M \leq -14`, after applying the apparent +magnitude cut. A value near one means that almost +all galaxies in this intrinsic magnitude range are included in the catalog. +A value near zero means that only a small part of that population is bright +enough to pass the apparent magnitude cut. + +The white contours mark fixed observable LF fractions of 0.2, 0.4, 0.6, and +0.8. In plain terms, the 0.8 contour is the line where the catalog keeps about +80 per cent of the chosen luminosity function integral, while the 0.2 contour +is where it keeps only about 20 per cent. + +This kind of summary is helpful when choosing survey depths or checking which +redshift range is safe for a magnitude-limited sample. + +.. plot:: + :include-source: True + :width: 560 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + z = np.linspace(0.05, 1.2, 160) + limits = np.linspace(22.0, 26.0, 120) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + completeness = [] + + for m_lim in limits: + f_obs = lf.catalog_completeness( + cosmo, + z, + m_lim=m_lim, + m_bright=-24.0, + m_faint=-14.0, + corrections=None, + ) + completeness.append(f_obs) + + completeness = np.asarray(completeness) + + fig, ax = plt.subplots(figsize=(7.2, 5.0)) + + mesh = ax.pcolormesh( + z, + limits, + completeness, + shading="auto", + cmap="cmr.guppy", + vmin=0.0, + vmax=1.0, + ) + + contour_levels = [0.2, 0.4, 0.6, 0.8] + contours = ax.contour( + z, + limits, + completeness, + levels=contour_levels, + colors="white", + linewidths=1.2, + ) + ax.clabel(contours, inline=True, fontsize=TICK_SIZE, fmt="%.1f") + + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"Apparent-magnitude limit $m_{\rm lim}$", fontsize=LABEL_SIZE) + ax.set_title(r"Completeness over $-24 \leq M \leq -14$", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + cbar = fig.colorbar(mesh, ax=ax) + cbar.set_label(r"Observable LF fraction", fontsize=LABEL_SIZE) + cbar.ax.tick_params(labelsize=TICK_SIZE) + + plt.tight_layout() + + +Comparing cosmologies +--------------------- + +The absolute magnitude limit depends on cosmology through the luminosity +distance and distance modulus. This means that the same apparent magnitude +limit can correspond to slightly different intrinsic magnitude cuts in +different cosmologies. + +The example below compares three cosmologies. The luminosity function is kept +fixed, so the plot isolates the effect of the cosmology-dependent magnitude +conversion. Users can replace the dictionary entries with any +:class:`pyccl.Cosmology` objects relevant to their analysis. + +The lower panel shows the residual relative to the reference cosmology, +:math:`\Omega_{\rm m}=0.30,\ h=0.70`. + +.. plot:: + :include-source: True + :width: 620 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + cosmologies = { + r"$\Omega_{\rm m}=0.25,\ h=0.70$": ccl.Cosmology( + Omega_c=0.20, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.30,\ h=0.70$": ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.35,\ h=0.70$": ccl.Cosmology( + Omega_c=0.30, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + } + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + m_limits = {} + + for label, cosmo in cosmologies.items(): + m_limits[label] = lf.absolute_magnitude_limit( + cosmo, + z, + m_lim=24.5, + corrections=None, + ) + + reference_label = r"$\Omega_{\rm m}=0.30,\ h=0.70$" + reference = m_limits[reference_label] + + fig, (ax_top, ax_bottom) = plt.subplots( + 2, + 1, + figsize=(7.0, 6.2), + sharex=True, + gridspec_kw={"height_ratios": [3, 1]}, + constrained_layout=True, + ) + + for (label, m_limit), color in zip(m_limits.items(), colors_red): + ax_top.plot(z, m_limit, lw=3, color=color, label=label) + ax_bottom.plot(z, m_limit - reference, lw=2.5, color=color) + + ax_top.invert_yaxis() + ax_top.set_ylabel(r"$M_{\rm lim}(z)$", fontsize=LABEL_SIZE) + ax_top.set_title(r"Catalog boundary for different cosmologies", fontsize=TITLE_SIZE) + ax_top.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + + ax_bottom.axhline(0.0, lw=1.0, color="0.3") + ax_bottom.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax_bottom.set_ylabel(r"$\Delta M_{\rm lim}$", fontsize=LABEL_SIZE) + + for ax in (ax_top, ax_bottom): + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + plt.tight_layout() + + +Completeness for different cosmologies +-------------------------------------- + +Changing cosmology also changes the completeness fraction through the +cosmology-dependent apparent-to-absolute magnitude conversion. + +This comparison keeps the luminosity function, apparent magnitude limit, and +intrinsic magnitude range fixed. Only the cosmology is changed. This makes the +plot useful for checking how sensitive the catalog-completeness calculation is +to the assumed distance-redshift relation. + +The lower panel shows the residual relative to the reference cosmology, +:math:`\Omega_{\rm m}=0.30,\ h=0.70`. + +.. plot:: + :include-source: True + :width: 620 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + cosmologies = { + r"$\Omega_{\rm m}=0.25,\ h=0.70$": ccl.Cosmology( + Omega_c=0.20, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.30,\ h=0.70$": ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.35,\ h=0.70$": ccl.Cosmology( + Omega_c=0.30, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + } + + z = np.linspace(0.05, 1.2, 250) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + completeness = {} + + for label, cosmo in cosmologies.items(): + completeness[label] = lf.catalog_completeness( + cosmo, + z, + m_lim=24.5, + m_bright=-24.0, + m_faint=-14.0, + corrections=None, + ) + + reference_label = r"$\Omega_{\rm m}=0.30,\ h=0.70$" + reference = completeness[reference_label] + + fig, (ax_top, ax_bottom) = plt.subplots( + 2, + 1, + figsize=(7.0, 6.2), + sharex=True, + gridspec_kw={"height_ratios": [3, 1]}, + constrained_layout=True, + ) + + for (label, f_obs), color in zip(completeness.items(), colors_red): + ax_top.plot(z, f_obs, lw=3, color=color, label=label) + ax_bottom.plot(z, f_obs - reference, lw=2.5, color=color) + + ax_top.set_ylim(-0.05, 1.05) + ax_top.set_ylabel("Observable LF fraction", fontsize=LABEL_SIZE) + ax_top.set_title(r"Cosmology dependence of LF completeness", fontsize=TITLE_SIZE) + ax_top.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + + ax_bottom.axhline(0.0, lw=1.0, color="0.3") + ax_bottom.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax_bottom.set_ylabel(r"$\Delta$ LF fraction", fontsize=LABEL_SIZE) + + for ax in (ax_top, ax_bottom): + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + plt.tight_layout() diff --git a/docs/examples/index.rst b/docs/examples/index.rst index 4dfcc1c..6237843 100644 --- a/docs/examples/index.rst +++ b/docs/examples/index.rst @@ -30,15 +30,43 @@ Working, executable examples for using LFKit. **Poggianti (1997) examples** ^^^ - Evaluate :math:`k(z)`, :math:`e(z)`, and :math:`k(z)+e(z)` from the + Evaluate :math:`k(z)`, :math:`e(z)`, and :math:`k(z)-e(z)` from the Poggianti tabulations; compare galaxy types and bands. +++ *Backends:* Poggianti (1997) + .. grid-item-card:: + :link: luminosity_function_examples + :link-type: doc + :shadow: md + + **Luminosity-function examples** + ^^^ + Build standard, evolving, and double Schechter luminosity functions, + evaluate :math:`\phi(M, z)`, and compare model behaviour. + + +++ + *API:* LuminosityFunction + + .. grid-item-card:: + :link: catalog_completeness_examples + :link-type: doc + :shadow: md + + **Catalog-completeness examples** + ^^^ + Compute observed and missing number densities for a magnitude-limited + catalog and visualize completeness as a function of redshift. + + +++ + *API:* LuminosityFunction + .. toctree:: :maxdepth: 1 :hidden: kcorrect_examples - poggianti_examples \ No newline at end of file + poggianti_examples + luminosity_function_examples + catalog_completeness_examples \ No newline at end of file diff --git a/docs/examples/luminosity_function_examples.rst b/docs/examples/luminosity_function_examples.rst new file mode 100644 index 0000000..b8f6e3e --- /dev/null +++ b/docs/examples/luminosity_function_examples.rst @@ -0,0 +1,959 @@ +.. |lfkitlogo| image:: /_static/logos/lfkit_logo-icon.png + :alt: LFKit logo + :width: 50px + +|lfkitlogo| Luminosity function examples +======================================== + +This page provides executable examples showing how to use +:class:`lfkit.LuminosityFunction` to construct, evaluate, and compare +luminosity function models. + +Luminosity functions describe how common galaxies are as a function of their +absolute magnitude. In these examples, brighter galaxies appear on the left +because astronomical magnitudes decrease for brighter objects. + +The luminosity function normalization sets the units. For example, if +:math:`\phi_\star` is supplied in comoving :math:`{\rm Mpc}^{-3}\,{\rm mag}^{-1}`, +then :math:`\phi(M, z)` has units of +:math:`{\rm Mpc}^{-3}\,{\rm mag}^{-1}`, and magnitude-integrated number +densities have units of :math:`{\rm Mpc}^{-3}`. + +The examples that connect apparent and absolute magnitude use +``corrections=None``. This means that no :math:`K`-correction or evolution +correction is applied. Users can pass their own correction model through the +``corrections`` argument, for example a :class:`lfkit.Corrections` object or any +compatible correction callable used by the LFKit magnitude-conversion methods. + +All examples below are executable via ``.. plot::``. + + +Standard Schechter luminosity function +-------------------------------------- + +A standard Schechter luminosity function has fixed parameters +:math:`\phi_\star`, :math:`M_\star`, and :math:`\alpha`. + +This plot shows the basic shape of a Schechter luminosity function as a +function of absolute magnitude. The function decreases rapidly at the bright end +and rises toward the faint end, reflecting the usual picture that very luminous +galaxies are rare while faint galaxies are more common. + +The y-axis is shown on a logarithmic scale because luminosity functions often +span several orders of magnitude. This makes both the bright-end cutoff and the +faint-end behaviour visible in the same figure. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + + absolute_mag = np.linspace(-24.0, -14.0, 500) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.1, + ) + + phi = lf.phi(absolute_mag) + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + ax.plot(absolute_mag, phi, lw=3, color=colors_blue[1]) + ax.set_yscale("log") + ax.invert_xaxis() + ax.set_xlabel("Absolute magnitude $M$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"$\phi(M)$ [$\mathrm{Mpc}^{-3}\,\mathrm{mag}^{-1}$]", fontsize=LABEL_SIZE) + ax.set_title("Standard Schechter luminosity function", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + plt.tight_layout() + + +Comparing Schechter slopes +-------------------------- + +Changing :math:`\alpha` modifies the faint-end behaviour of the luminosity +function. + +This comparison shows how the faint-end slope changes the abundance of faint +galaxies while keeping the other Schechter parameters fixed. More negative +values of :math:`\alpha` produce a steeper rise toward faint magnitudes. + +This is useful because the faint-end slope often controls how strongly low +luminosity galaxies contribute to integrated quantities, such as number density +or luminosity density. Even if the bright end is almost unchanged, the total +abundance can change noticeably when the faint end is modified. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + absolute_mag = np.linspace(-24.0, -14.0, 500) + alphas = [-0.8, -1.1, -1.4] + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + + for alpha, color in zip(alphas, colors_red): + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=alpha, + ) + phi = lf.phi(absolute_mag) + ax.plot(absolute_mag, phi, lw=3, color=color, label=rf"$\alpha={alpha}$") + + ax.set_yscale("log") + ax.invert_xaxis() + ax.set_xlabel("Absolute magnitude $M$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"$\phi(M)$ [$\mathrm{Mpc}^{-3}\,\mathrm{mag}^{-1}$]", fontsize=LABEL_SIZE) + ax.set_title("Effect of the faint-end slope", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Evolving Schechter luminosity function +-------------------------------------- + +An evolving Schechter model allows the luminosity function parameters to vary +with redshift. + +This plot compares the luminosity function at several redshifts. Instead of +using one fixed curve for all epochs, the evolving model allows the amplitude +and characteristic magnitude to change with redshift. + +This kind of plot is useful for visualizing galaxy evolution in a compact way. +Changes in normalization alter the overall abundance, while shifts in +:math:`M_\star` move the turnover of the luminosity function toward brighter or +fainter magnitudes. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + + absolute_mag = np.linspace(-24.0, -14.0, 500) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + redshifts = [0.1, 0.5, 1.0] + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + + for z_value, color in zip(redshifts, colors_blue): + z = np.full_like(absolute_mag, z_value) + phi = lf.phi(absolute_mag, z) + ax.plot(absolute_mag, phi, lw=3, color=color, label=rf"$z={z_value}$") + + ax.set_yscale("log") + ax.invert_xaxis() + ax.set_xlabel("Absolute magnitude $M$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"$\phi(M, z)$ [$\mathrm{Mpc}^{-3}\,\mathrm{mag}^{-1}$]", fontsize=LABEL_SIZE) + ax.set_title("Evolving Schechter luminosity function", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Inspecting evolving parameters +------------------------------ + +The evolving luminosity function parameters can also be evaluated directly. + +This plot shows the redshift evolution of the parameters that define the +luminosity function. Instead of plotting :math:`\phi(M, z)` itself, it shows how +:math:`\phi_\star`, :math:`M_\star`, and :math:`\alpha` change with redshift. + +This is useful when checking whether an evolving model behaves as expected +before using it in a larger calculation. For example, it can help verify that +the normalization evolves smoothly, that :math:`M_\star` shifts in the intended +direction, and that the faint-end slope remains in a sensible range. + +.. plot:: + :include-source: True + :width: 620 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + z = np.linspace(0.0, 1.5, 300) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="linear", + alpha_kwargs={"alpha_0": -1.1, "alpha_1": -0.1, "z_ref": 0.1}, + ) + + phi_star, m_star, alpha = lf.parameters(z) + + fig, axes = plt.subplots( + 3, + 1, + figsize=(7.0, 7.2), + sharex=True, + constrained_layout=True, + ) + + axes[0].plot(z, phi_star, lw=3, color=colors_blue[1]) + axes[0].set_ylabel(r"$\phi_\star$", fontsize=LABEL_SIZE) + axes[0].ticklabel_format(axis="y", style="sci", scilimits=(0, 0)) + axes[0].set_title("Evolving LF parameters", fontsize=TITLE_SIZE) + + axes[1].plot(z, m_star, lw=3, color=colors_red[1]) + axes[1].set_ylabel(r"$M_\star$", fontsize=LABEL_SIZE) + + axes[2].plot(z, alpha, lw=3, color=colors_blue[2]) + axes[2].set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + axes[2].set_ylabel(r"$\alpha$", fontsize=LABEL_SIZE) + + for ax in axes: + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + plt.tight_layout() + + +Double Schechter luminosity function +------------------------------------ + +A double Schechter-style model can be used when the luminosity function needs +additional structure beyond the standard Schechter form. + +This plot compares a standard Schechter model with a double-Schechter-style +model. The double model adds extra flexibility around the faint end, where a +single power-law slope may not describe the full galaxy population well. + +This type of comparison is useful when testing whether a simple one-component +model is sufficient. If the two curves differ strongly at faint magnitudes, +integrated quantities that depend on faint galaxies may also differ. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + absolute_mag = np.linspace(-24.0, -14.0, 500) + + standard = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.1, + ) + + double = LuminosityFunction.double_schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.0, + beta=-1.5, + m_transition=-18.0, + ) + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + + ax.plot( + absolute_mag, + standard.phi(absolute_mag), + lw=3, + color=colors_blue[1], + label="Standard Schechter", + ) + ax.plot( + absolute_mag, + double.phi(absolute_mag), + lw=3, + color=colors_red[1], + label="Double Schechter", + ) + + ax.set_yscale("log") + ax.invert_xaxis() + ax.set_xlabel("Absolute magnitude $M$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"$\phi(M)$ [$\mathrm{Mpc}^{-3}\,\mathrm{mag}^{-1}$]", fontsize=LABEL_SIZE) + ax.set_title("Standard and double Schechter models", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Integrated number density +------------------------- + +A luminosity function can be integrated over magnitude to estimate the number +density of galaxies brighter than a chosen absolute-magnitude limit. + +This plot shows how the cumulative number density changes as progressively +fainter galaxies are included. At very bright limits, only rare luminous +galaxies contribute. As the magnitude limit becomes fainter, more galaxies are +included and the integrated number density increases. + +This is one of the most common ways a luminosity function enters survey +calculations. Instead of using the value of :math:`\phi(M)` at one magnitude, +the model is integrated over the part of the galaxy population that the sample +selects. + +.. plot:: + :include-source: True + :width: 520 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.1, + ) + + magnitude_limits = np.linspace(-23.0, -15.0, 120) + + number_density = lf.integrated_number_density( + z=0.0, + m_bright=-25.0, + m_faint=magnitude_limits, + n_m=800, + ) + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + ax.plot(magnitude_limits, number_density, lw=3, color=colors_blue[1]) + ax.set_yscale("log") + ax.invert_xaxis() + ax.set_xlabel(r"Faint absolute-magnitude limit $M_{\rm lim}$", fontsize=LABEL_SIZE) + ax.set_ylabel(r"$n(M < M_{\rm lim})$ [$\mathrm{Mpc}^{-3}$]", fontsize=LABEL_SIZE) + ax.set_title("Integrated number density", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + plt.tight_layout() + + +Magnitude-redshift luminosity function surface +---------------------------------------------- + +For an evolving luminosity function, the abundance depends on both absolute +magnitude and redshift. + +This plot shows the luminosity function amplitude across the +magnitude-redshift plane. The horizontal direction shows which galaxies are +bright or faint, while the vertical direction shows how the model changes with +redshift. + +The filled colour scale shows :math:`\log_{10}\phi(M, z)`. The white contours +mark constant :math:`\log_{10}\phi(M, z)` levels at -5, -4, -3, and -2. These +contours make it easier to see where equal-abundance regions sit in the +magnitude-redshift plane. + +This view is helpful for checking the full two-dimensional behaviour of an +evolving model. It makes it easier to see whether the bright end, faint end, and +redshift evolution combine smoothly across the range of interest. + +.. plot:: + :include-source: True + :width: 560 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + absolute_mag = np.linspace(-24.0, -18.0, 220) + redshift = np.linspace(0.0, 1.5, 180) + + mag_grid, z_grid = np.meshgrid(absolute_mag, redshift) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.8}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 1.0, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + phi = lf.phi(mag_grid, z_grid) + log_phi = np.log10(phi) + + fig, ax = plt.subplots(figsize=(7.2, 5.0)) + + mesh = ax.pcolormesh( + absolute_mag, + redshift, + log_phi, + shading="auto", + cmap="cmr.guppy", + ) + + contour_levels = [-5.0, -4.0, -3.0, -2.0] + contours = ax.contour( + absolute_mag, + redshift, + log_phi, + levels=contour_levels, + colors="white", + linewidths=1.2, + ) + ax.clabel(contours, inline=True, fontsize=TICK_SIZE, fmt=r"$10^{%.0f}$") + + ax.invert_xaxis() + ax.set_xlabel("Absolute magnitude $M$", fontsize=LABEL_SIZE) + ax.set_ylabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_title("Evolving luminosity function", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + cbar = fig.colorbar(mesh, ax=ax) + cbar.set_label( + r"$\log_{10}\phi(M, z)$ [$\log_{10}(\mathrm{Mpc}^{-3}\,\mathrm{mag}^{-1})$]", + fontsize=LABEL_SIZE, + ) + cbar.ax.tick_params(labelsize=TICK_SIZE) + + plt.tight_layout() + + +Magnitude-limited LF-weighted redshift trend +-------------------------------------------- + +A survey magnitude limit selects different parts of the luminosity function at +different redshifts. + +This example shows LF-weighted redshift trends for several fixed +apparent-magnitude limits. The absolute-magnitude limit is computed from the +cosmology-dependent distance modulus. At lower redshift, a flux-limited sample +can include relatively faint galaxies. At higher redshift, the same +apparent-magnitude limit corresponds to a brighter absolute-magnitude cut, so +only intrinsically brighter galaxies remain in the selected sample. + +The result is not intended to be a full survey :math:`n(z)`, because it does not +include the survey volume element. Instead, it shows the luminosity function +selection factor that later enters magnitude-limited redshift-distribution +calculations. + +.. plot:: + :include-source: True + :width: 620 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_blue = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.72, 0.96)) + + cosmo = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.6}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 0.8, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + redshift = np.linspace(0.05, 1.5, 160) + apparent_limits = [23.5, 24.5, 25.5] + + fig, ax = plt.subplots(figsize=(7.0, 5.0)) + + for m_lim, color in zip(apparent_limits, colors_blue): + m_limit = lf.absolute_magnitude_limit( + cosmo, + redshift, + m_lim=m_lim, + corrections=None, + ) + + lf_selection = lf.integrated_number_density( + z=redshift, + m_bright=-25.0, + m_faint=m_limit, + n_m=700, + ) + lf_selection /= np.trapezoid(lf_selection, redshift) + + ax.plot( + redshift, + lf_selection, + lw=3, + color=color, + label=rf"$m_{{\rm lim}}={m_lim}$", + ) + + ax.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax.set_ylabel("LF selection", fontsize=LABEL_SIZE) + ax.set_title("Magnitude-limited LF selection", fontsize=TITLE_SIZE) + ax.tick_params(axis="both", labelsize=TICK_SIZE) + ax.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + plt.tight_layout() + + +Cosmology dependence of the absolute-magnitude limit +---------------------------------------------------- + +The apparent-to-absolute magnitude conversion depends on cosmology through the +luminosity distance and distance modulus. + +This plot compares the absolute-magnitude limit implied by the same apparent +magnitude cut in several cosmologies. The luminosity function is not used in +this plot; the comparison isolates the selection boundary itself. + +Users can replace the entries in the ``cosmologies`` dictionary with any +:class:`pyccl.Cosmology` objects relevant to their own analysis. + +.. plot:: + :include-source: True + :width: 620 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + cosmologies = { + r"$\Omega_{\rm m}=0.25,\ h=0.70$": ccl.Cosmology( + Omega_c=0.20, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.30,\ h=0.70$": ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.35,\ h=0.70$": ccl.Cosmology( + Omega_c=0.30, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + } + + z = np.linspace(0.05, 1.5, 250) + + lf = LuminosityFunction.schechter( + phi_star=1.0e-3, + m_star=-20.5, + alpha=-1.1, + ) + + m_limits = {} + + for label, cosmo in cosmologies.items(): + m_limits[label] = lf.absolute_magnitude_limit( + cosmo, + z, + m_lim=24.5, + corrections=None, + ) + + reference_label = r"$\Omega_{\rm m}=0.30,\ h=0.70$" + reference = m_limits[reference_label] + + fig, (ax_top, ax_bottom) = plt.subplots( + 2, + 1, + figsize=(7.0, 6.2), + sharex=True, + gridspec_kw={"height_ratios": [3, 1]}, + constrained_layout=True, + ) + + for (label, m_limit), color in zip(m_limits.items(), colors_red): + ax_top.plot(z, m_limit, lw=3, color=color, label=label) + ax_bottom.plot(z, m_limit - reference, lw=2.5, color=color) + + ax_top.invert_yaxis() + ax_top.set_ylabel(r"$M_{\rm lim}(z)$", fontsize=LABEL_SIZE) + ax_top.set_title(r"Cosmology dependence of $M_{\rm lim}(z)$", fontsize=TITLE_SIZE) + ax_top.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + + ax_bottom.axhline(0.0, lw=1.0, color="0.3") + ax_bottom.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax_bottom.set_ylabel(r"$\Delta M_{\rm lim}$", fontsize=LABEL_SIZE) + + for ax in (ax_top, ax_bottom): + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + plt.tight_layout() + + +Cosmology dependence of LF selection +------------------------------------ + +The LF-weighted selection factor also depends on cosmology because the same +apparent-magnitude cut maps to a different absolute-magnitude limit. + +This example keeps the luminosity function and apparent-magnitude limit fixed, +then changes only the cosmology. The curves are normalized to unit integral over +redshift, so the comparison emphasizes changes in shape rather than absolute +normalization. + +The lower panel shows the residual relative to the reference cosmology, +:math:`\Omega_{\rm m}=0.30,\ h=0.70`. + +This is still not a full survey :math:`n(z)`, because it does not include the +cosmological volume element. It is the luminosity function selection factor +alone. + +.. plot:: + :include-source: True + :width: 620 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + cosmologies = { + r"$\Omega_{\rm m}=0.25,\ h=0.70$": ccl.Cosmology( + Omega_c=0.20, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.30,\ h=0.70$": ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.35,\ h=0.70$": ccl.Cosmology( + Omega_c=0.30, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + } + + redshift = np.linspace(0.05, 1.5, 160) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.6}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 0.8, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + selections = {} + + for label, cosmo in cosmologies.items(): + m_limit = lf.absolute_magnitude_limit( + cosmo, + redshift, + m_lim=24.5, + corrections=None, + ) + + lf_selection = lf.integrated_number_density( + z=redshift, + m_bright=-25.0, + m_faint=m_limit, + n_m=700, + ) + lf_selection /= np.trapezoid(lf_selection, redshift) + + selections[label] = lf_selection + + reference_label = r"$\Omega_{\rm m}=0.30,\ h=0.70$" + reference = selections[reference_label] + + fig, (ax_top, ax_bottom) = plt.subplots( + 2, + 1, + figsize=(7.0, 6.2), + sharex=True, + gridspec_kw={"height_ratios": [3, 1]}, + constrained_layout=True, + ) + + for (label, lf_selection), color in zip(selections.items(), colors_red): + ax_top.plot(redshift, lf_selection, lw=3, color=color, label=label) + ax_bottom.plot(redshift, lf_selection - reference, lw=2.5, color=color) + + ax_top.set_ylabel("LF selection", fontsize=LABEL_SIZE) + ax_top.set_title(r"Cosmology dependence of LF selection", fontsize=TITLE_SIZE) + ax_top.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + + ax_bottom.axhline(0.0, lw=1.0, color="0.3") + ax_bottom.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax_bottom.set_ylabel(r"$\Delta$ selection", fontsize=LABEL_SIZE) + + for ax in (ax_top, ax_bottom): + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + plt.tight_layout() + + +Cosmology and volume-weighted LF redshift trend +----------------------------------------------- + +A full LF-based redshift trend for a magnitude-limited sample should include +both the luminosity function selection and the cosmological volume element. + +This example multiplies the magnitude-integrated luminosity function by a +simple comoving-volume weight per steradian, :math:`\chi^2(z) / H(z)`, up to an +overall constant. The result is closer to the ingredient used in LF-dependent +:math:`n(z)` construction. + +The curves are normalized to unit integral over redshift, so the comparison +shows how cosmology changes the shape of the redshift trend. The lower panel +shows the residual relative to the reference cosmology, +:math:`\Omega_{\rm m}=0.30,\ h=0.70`. + +The absolute normalization depends on survey area, LF normalization, and the +exact volume convention used by the calling analysis. + +.. plot:: + :include-source: True + :width: 620 + + import numpy as np + import matplotlib.pyplot as plt + import cmasher as cmr + import pyccl as ccl + + from lfkit import LuminosityFunction + + LABEL_SIZE = 15 + TICK_SIZE = 13 + TITLE_SIZE = 17 + LEGEND_SIZE = 15 + + colors_red = cmr.take_cmap_colors("cmr.guppy", 3, cmap_range=(0.03, 0.26)) + + cosmologies = { + r"$\Omega_{\rm m}=0.25,\ h=0.70$": ccl.Cosmology( + Omega_c=0.20, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.30,\ h=0.70$": ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + r"$\Omega_{\rm m}=0.35,\ h=0.70$": ccl.Cosmology( + Omega_c=0.30, + Omega_b=0.05, + h=0.70, + sigma8=0.8, + n_s=0.96, + transfer_function="bbks", + matter_power_spectrum="linear", + ), + } + + redshift = np.linspace(0.05, 1.5, 160) + + lf = LuminosityFunction.evolving_schechter( + phi_model="linear_p", + phi_kwargs={"phi_0_star": 1.0e-3, "p": 0.6}, + m_star_model="linear_q", + m_star_kwargs={"m_0_star": -20.5, "q": 0.8, "z_ref": 0.1}, + alpha_model="constant", + alpha_kwargs={"alpha": -1.1}, + ) + + trends = {} + + for label, cosmo in cosmologies.items(): + scale_factor = 1.0 / (1.0 + redshift) + + chi = ccl.comoving_radial_distance(cosmo, scale_factor) + h_over_h0 = ccl.h_over_h0(cosmo, scale_factor) + + volume_weight = chi**2 / h_over_h0 + + m_limit = lf.absolute_magnitude_limit( + cosmo, + redshift, + m_lim=24.5, + corrections=None, + ) + + lf_selection = lf.integrated_number_density( + z=redshift, + m_bright=-25.0, + m_faint=m_limit, + n_m=700, + ) + + weighted_trend = volume_weight * lf_selection + weighted_trend /= np.trapezoid(weighted_trend, redshift) + + trends[label] = weighted_trend + + reference_label = r"$\Omega_{\rm m}=0.30,\ h=0.70$" + reference = trends[reference_label] + + fig, (ax_top, ax_bottom) = plt.subplots( + 2, + 1, + figsize=(7.0, 6.2), + sharex=True, + gridspec_kw={"height_ratios": [3, 1]}, + constrained_layout=True, + ) + + for (label, weighted_trend), color in zip(trends.items(), colors_red): + ax_top.plot(redshift, weighted_trend, lw=3, color=color, label=label) + ax_bottom.plot(redshift, weighted_trend - reference, lw=2.5, color=color) + + ax_top.set_ylabel("Volume-weighted LF trend", fontsize=LABEL_SIZE) + ax_top.set_title(r"Cosmology dependence with volume weighting", fontsize=TITLE_SIZE) + ax_top.legend(frameon=True, fontsize=LEGEND_SIZE, loc="best") + + ax_bottom.axhline(0.0, lw=1.0, color="0.3") + ax_bottom.set_xlabel("Redshift $z$", fontsize=LABEL_SIZE) + ax_bottom.set_ylabel(r"$\Delta$ trend", fontsize=LABEL_SIZE) + + for ax in (ax_top, ax_bottom): + ax.tick_params(axis="both", labelsize=TICK_SIZE) + + plt.tight_layout() diff --git a/docs/index.rst b/docs/index.rst index 9dda441..c361bee 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,20 +1,64 @@ -.. LFKit documentation master file, created by - sphinx-quickstart on Fri Mar 6 21:59:59 2026. - You can adapt this file completely to your liking, but it should at least - contain the root `toctree` directive. +.. |lfkitlogo| image:: /_static/logos/lfkit_logo-icon.png + :alt: LFKit logo + :width: 60px -LFKit documentation -=================== +|lfkitlogo| LFKit +================= -Add your content using ``reStructuredText`` syntax. See the -`reStructuredText `_ -documentation for details. +**LFKit** is a toolkit for modelling galaxy luminosity functions, +photometric corrections, and magnitude-limited catalog completeness. +It provides a clean interface for building theoretical luminosity functions, +including redshift-dependent parameter models, and for connecting apparent +magnitude limits to observable and missing galaxy number densities. + +LFKit is designed to be science-use-case agnostic: the same luminosity-function +machinery can be used in photometric-redshift modelling, intrinsic-alignment +modelling, cluster science, GW-cosmology catalog completeness, or any other +analysis that needs luminosity-function-based number densities. + +Getting started +--------------- + +Use the examples section for runnable workflows with plots, or the API +reference for detailed documentation of the public classes and functions. + +.. grid:: 3 + :gutter: 2 + + .. grid-item-card:: + :link: about/index + :link-type: doc + :shadow: md + + **About** + ^^^ + Overview of the package, scope, and design choices. + + .. grid-item-card:: + :link: examples/index + :link-type: doc + :shadow: md + + **Examples** + ^^^ + Runnable examples for luminosity functions, corrections, and catalog + completeness. + + .. grid-item-card:: + :link: api/index + :link-type: doc + :shadow: md + + **API reference** + ^^^ + Public classes, functions, and modules. .. toctree:: :maxdepth: 2 - :caption: Contents: + :caption: Documentation + :hidden: about/index examples/index - api/index + api/index \ No newline at end of file From b49c508e492e2f221412af27a220fc5640ec5837 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:26:15 -0400 Subject: [PATCH 09/10] Update correction backend tests --- tests/test_corrections_kcorrect_backend.py | 1 - tests/test_corrections_responses.py | 1 - 2 files changed, 2 deletions(-) diff --git a/tests/test_corrections_kcorrect_backend.py b/tests/test_corrections_kcorrect_backend.py index d541598..1cae3d2 100644 --- a/tests/test_corrections_kcorrect_backend.py +++ b/tests/test_corrections_kcorrect_backend.py @@ -2,7 +2,6 @@ from __future__ import annotations -from pathlib import Path from types import SimpleNamespace import pytest diff --git a/tests/test_corrections_responses.py b/tests/test_corrections_responses.py index 93b233b..bf30002 100644 --- a/tests/test_corrections_responses.py +++ b/tests/test_corrections_responses.py @@ -7,7 +7,6 @@ import numpy as np import pytest -# Adjust import path to your package layout from lfkit.corrections.responses import write_kcorrect_response From 74af3451f02c6cbf68507158338487aa6931c904 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Mon, 11 May 2026 18:33:57 -0400 Subject: [PATCH 10/10] Remove pandas dependency from fake catalog tests --- ...st_photometry_completeness_fake_catalog.py | 28 +++++++++++-------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/tests/test_photometry_completeness_fake_catalog.py b/tests/test_photometry_completeness_fake_catalog.py index bbac27f..fa068b7 100644 --- a/tests/test_photometry_completeness_fake_catalog.py +++ b/tests/test_photometry_completeness_fake_catalog.py @@ -5,7 +5,6 @@ from importlib.resources import files import numpy as np -import pandas as pd import pyccl as ccl from lfkit.photometry.catalog_completeness import ( @@ -41,7 +40,7 @@ def make_cosmology() -> ccl.Cosmology: ) -def load_fake_catalog() -> pd.DataFrame: +def load_fake_catalog() -> np.ndarray: """Return the packaged fake magnitude-limited catalog.""" path = ( files("lfkit") @@ -49,7 +48,7 @@ def load_fake_catalog() -> pd.DataFrame: / "demo_catalogs" / "fake_magnitude_limited_catalog.csv" ) - return pd.read_csv(path) + return np.genfromtxt(path, delimiter=",", names=True, dtype=None, encoding="utf-8") def test_fake_catalog_is_packaged() -> None: @@ -57,9 +56,14 @@ def test_fake_catalog_is_packaged() -> None: catalog = load_fake_catalog() assert len(catalog) == 200 - assert {"galaxy_id", "ra_deg", "dec_deg", "z", "m_app", "band"} <= set( - catalog.columns - ) + assert set(catalog.dtype.names) == { + "galaxy_id", + "ra_deg", + "dec_deg", + "z", + "m_app", + "band", + } def test_completeness_fraction_on_fake_catalog_redshifts() -> None: @@ -67,7 +71,7 @@ def test_completeness_fraction_on_fake_catalog_redshifts() -> None: catalog = load_fake_catalog() cosmo = make_cosmology() - z = catalog["z"].to_numpy() + z = np.asarray(catalog["z"], dtype=float) completeness = catalog_completeness_fraction( cosmo, @@ -90,7 +94,7 @@ def test_observed_and_missing_fractions_sum_to_one() -> None: catalog = load_fake_catalog() cosmo = make_cosmology() - z = catalog["z"].to_numpy() + z = np.asarray(catalog["z"], dtype=float) completeness = catalog_completeness_fraction( cosmo, @@ -124,7 +128,7 @@ def test_deeper_catalog_limit_increases_completeness() -> None: catalog = load_fake_catalog() cosmo = make_cosmology() - z = catalog["z"] + z = np.asarray(catalog["z"], dtype=float) shallow = catalog_completeness_fraction( cosmo, @@ -154,7 +158,7 @@ def test_deeper_catalog_limit_decreases_missing_density() -> None: catalog = load_fake_catalog() cosmo = make_cosmology() - z = catalog["z"] + z = np.asarray(catalog["z"], dtype=float) shallow = missing_number_density( cosmo, @@ -183,7 +187,7 @@ def test_deeper_catalog_limit_increases_observed_density() -> None: catalog = load_fake_catalog() cosmo = make_cosmology() - z = catalog["z"] + z = np.asarray(catalog["z"], dtype=float) shallow = observed_number_density( cosmo, @@ -213,7 +217,7 @@ def test_observed_and_missing_number_densities_sum_to_total_density() -> None: catalog = load_fake_catalog() cosmo = make_cosmology() - z = catalog["z"].to_numpy() + z = np.asarray(catalog["z"], dtype=float) total = observed_number_density( cosmo,