From a60b936cc096942aaf17ffc87c00b77216b55ad5 Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Sun, 14 Jun 2026 23:32:06 -0400 Subject: [PATCH 1/2] added additional lf models --- .../conditional_models.py | 3 + src/lfkit/luminosity_functions/integrals.py | 6 +- .../luminosity_functions/models/gamma.py | 145 ++++ .../models/non_parametric.py | 718 ++++++++++++++++++ .../luminosity_functions/models/saunders.py | 230 ++++++ .../luminosity_functions/models/schechter.py | 122 +++ src/lfkit/utils/validators.py | 178 +++++ tests/test_lumfuncs_models_gamma.py | 0 tests/test_lumfuncs_models_non_parametric.py | 0 tests/test_lumfuncs_models_saunders.py | 0 10 files changed, 1400 insertions(+), 2 deletions(-) create mode 100644 src/lfkit/luminosity_functions/models/gamma.py create mode 100644 src/lfkit/luminosity_functions/models/non_parametric.py create mode 100644 src/lfkit/luminosity_functions/models/saunders.py create mode 100644 tests/test_lumfuncs_models_gamma.py create mode 100644 tests/test_lumfuncs_models_non_parametric.py create mode 100644 tests/test_lumfuncs_models_saunders.py diff --git a/src/lfkit/luminosity_functions/conditional_models.py b/src/lfkit/luminosity_functions/conditional_models.py index 426d6cd5..093ddfbd 100644 --- a/src/lfkit/luminosity_functions/conditional_models.py +++ b/src/lfkit/luminosity_functions/conditional_models.py @@ -72,6 +72,9 @@ def _conditional_model_name(name: str) -> str: Returns: Conditional model name prefixed with ``"conditional_"``. """ + if name.endswith("_lf"): + name = name.removesuffix("_lf") + return f"conditional_{name}" diff --git a/src/lfkit/luminosity_functions/integrals.py b/src/lfkit/luminosity_functions/integrals.py index 9c3254bc..d3025703 100644 --- a/src/lfkit/luminosity_functions/integrals.py +++ b/src/lfkit/luminosity_functions/integrals.py @@ -320,9 +320,11 @@ def luminosity_weight( raise ValueError("m_reference must be finite.") absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") - weight = 10.0 ** (-0.4 * (absolute_mag_arr - m_reference)) - weight = np.clip(weight, 1e-300, 1e300) + exponent = -0.4 * (absolute_mag_arr - m_reference) + exponent = np.clip(exponent, -300.0, 300.0) + + weight = 10.0**exponent return np.asarray(weight, dtype=float) diff --git a/src/lfkit/luminosity_functions/models/gamma.py b/src/lfkit/luminosity_functions/models/gamma.py new file mode 100644 index 00000000..b899a4c8 --- /dev/null +++ b/src/lfkit/luminosity_functions/models/gamma.py @@ -0,0 +1,145 @@ +r"""Gamma-family luminosity function models. + +This module provides standalone functions for evaluating gamma-like luminosity +function models in rest frame absolute magnitude space. + +Implemented models: + - Gamma luminosity function + - Generalized gamma luminosity function +""" + +from __future__ import annotations + +import numpy as np + +from lfkit.utils.integrators import safe_power10 +from lfkit.utils.types import FloatArray, FloatInput, ParameterValue +from lfkit.utils.validators import validate_array + + +__all__ = [ + "gamma_lf", + "generalized_gamma_lf", +] + + +def gamma_lf( + absolute_mag: FloatInput, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, +) -> FloatArray: + r"""Return a gamma luminosity function in magnitude space. + + This computes + + .. math:: + + \phi(M) = + 0.4 \ln(10) \, \phi_\star \, + x^{\alpha + 1} \exp(-x), + + where + + .. math:: + + x = 10^{-0.4(M - M_\star)}. + + This is mathematically equivalent to the standard Schechter luminosity + function, but is provided as a gamma-family model for users who want to + describe the luminosity function as a gamma distribution in luminosity + space. + + Args: + absolute_mag: Absolute magnitude value(s). + phi_star: Gamma luminosity function normalization. + m_star: Characteristic magnitude. + alpha: Low-luminosity power-law slope. + + Returns: + NumPy array containing the gamma luminosity function evaluated at + ``absolute_mag``. + + Raises: + ValueError: If ``phi_star`` is negative. + """ + absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") + phi_star_arr = validate_array(phi_star, name="phi_star") + alpha_arr = validate_array(alpha, name="alpha") + + if np.any(phi_star_arr < 0.0): + raise ValueError("phi_star must be non-negative.") + + x = safe_power10(-0.4 * (absolute_mag_arr - m_star)) + + phi = 0.4 * np.log(10.0) * phi_star_arr * x ** (alpha_arr + 1.0) * np.exp(-x) + + return np.asarray(phi, dtype=float) + + +def generalized_gamma_lf( + absolute_mag: FloatInput, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, + beta: ParameterValue, +) -> FloatArray: + r"""Return a generalized gamma luminosity function in magnitude space. + + This computes + + .. math:: + + \phi(M) = + 0.4 \ln(10) \, \beta \, \phi_\star \, + x^{\alpha + 1} \exp(-x^\beta), + + where + + .. math:: + + x = 10^{-0.4(M - M_\star)}. + + The gamma or Schechter form is recovered when ``beta = 1``. Values + ``beta < 1`` produce a shallower bright-end cutoff, while values + ``beta > 1`` produce a sharper bright-end cutoff. + + Args: + absolute_mag: Absolute magnitude value(s). + phi_star: Generalized gamma luminosity function normalization. + m_star: Characteristic magnitude. + alpha: Low-luminosity power-law slope. + beta: Positive bright-end cutoff shape parameter. + + Returns: + NumPy array containing the generalized gamma luminosity function + evaluated at ``absolute_mag``. + + Raises: + ValueError: If ``phi_star`` is negative or if ``beta`` is not positive. + """ + absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") + phi_star_arr = validate_array(phi_star, name="phi_star") + alpha_arr = validate_array(alpha, name="alpha") + beta_arr = validate_array(beta, name="beta") + + if np.any(phi_star_arr < 0.0): + raise ValueError("phi_star must be non-negative.") + + if np.any(beta_arr <= 0.0): + raise ValueError("beta must be positive.") + + x = safe_power10(-0.4 * (absolute_mag_arr - m_star)) + + phi = ( + 0.4 + * np.log(10.0) + * beta_arr + * phi_star_arr + * x ** (alpha_arr + 1.0) + * np.exp(-(x**beta_arr)) + ) + + return np.asarray(phi, dtype=float) diff --git a/src/lfkit/luminosity_functions/models/non_parametric.py b/src/lfkit/luminosity_functions/models/non_parametric.py new file mode 100644 index 00000000..1a3d96a4 --- /dev/null +++ b/src/lfkit/luminosity_functions/models/non_parametric.py @@ -0,0 +1,718 @@ +"""Tabulated and binned luminosity function models. + +This module provides non-parametric luminosity function models. Unlike +Schechter, Saunders, Gaussian, or power-law forms, these models do not assume a +fixed analytic shape. Instead, the luminosity function is supplied directly as +values on a magnitude grid or inside magnitude bins. + +Two related representations are provided: + +``tabulated_lf`` + Interpolates luminosity function values sampled at magnitude grid points. + This is useful when an LF has already been measured, simulated, or fitted + elsewhere and should be evaluated smoothly between tabulated points. + +``binned_lf`` + Treats the luminosity function as piecewise constant inside magnitude bins. + This is useful for directly representing measured binned luminosity + functions without imposing interpolation structure inside each bin. + +Redshift-dependent versions are also provided. These allow the LF to vary over +a two-dimensional grid in absolute magnitude and redshift. + +These models are especially useful when the bright end, faint end, or redshift +evolution should be data-driven rather than forced into a Schechter-like or +Saunders-like form. +""" + +from __future__ import annotations + +import numpy as np + +from lfkit.utils.types import FloatArray, FloatInput, ParameterValue +from lfkit.utils.validators import ( + validate_array, + validate_tabulated_grid, + validate_binned_grid, + validate_2d_tabulated_grid, + validate_2d_binned_grid, +) + + +__all__ = [ + "tabulated_lf", + "binned_lf", + "redshift_tabulated_lf", + "redshift_binned_lf", + "distance_tabulated_lf", + "distance_binned_lf", +] + + +def tabulated_lf( + absolute_mag: FloatInput, + *, + magnitude_grid: ParameterValue, + phi_grid: ParameterValue, + fill_value: float = 0.0, + log_phi: bool = False, +) -> FloatArray: + r"""Return a luminosity function interpolated from tabulated values. + + This evaluates a non-parametric luminosity function by interpolating + tabulated values, + + .. math:: + + \phi(M_i) = \phi_i, + + where ``magnitude_grid`` contains the sampled absolute magnitudes + :math:`M_i` and ``phi_grid`` contains the corresponding luminosity function + values :math:`\phi_i`. + + For magnitudes between tabulated points, this function performs linear + interpolation, + + .. math:: + + \phi(M) = + \mathrm{interp}\left(M; \{M_i\}, \{\phi_i\}\right). + + If ``log_phi`` is True, interpolation is instead performed in + :math:`\log_{10}\phi`, + + .. math:: + + \log_{10}\phi(M) = + \mathrm{interp}\left( + M; \{M_i\}, \{\log_{10}\phi_i\} + \right), + + and the result is transformed back to linear space. Logarithmic + interpolation is often better for luminosity functions because LF values + can vary by many orders of magnitude across the magnitude range. + + This model is useful when the luminosity function is known from a table, + simulation, external fit, or observational measurement and no analytic + Schechter-like form should be imposed. It gives a smooth representation of + the tabulated curve, but it does not extrapolate beyond the supplied + magnitude range. Values outside the table are set to ``fill_value``. + + Args: + absolute_mag: Absolute magnitude value(s) where the luminosity function + should be evaluated. + magnitude_grid: One-dimensional strictly increasing grid of absolute + magnitudes where the LF is tabulated. + phi_grid: One-dimensional array of non-negative LF values evaluated at + ``magnitude_grid``. + fill_value: Value returned outside the tabulated magnitude range. + log_phi: If True, interpolate in ``log10(phi)`` instead of directly in + ``phi``. + + Returns: + NumPy array containing the interpolated luminosity function evaluated + at ``absolute_mag``. + + Raises: + ValueError: If the magnitude grid is not one-dimensional, if the grid + is not strictly increasing, if the grid and value arrays have + inconsistent lengths, if any LF value is negative, if + ``fill_value`` is not finite, or if ``log_phi`` is True and any + tabulated LF value is not positive. + """ + absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") + magnitude_grid_arr, phi_grid_arr = validate_tabulated_grid( + magnitude_grid, + phi_grid, + coordinate_name="magnitude_grid", + values_name="phi_grid", + positive_values=log_phi, + ) + + if not np.isfinite(fill_value): + raise ValueError("fill_value must be finite.") + + if log_phi: + if np.any(phi_grid_arr <= 0.0): + raise ValueError("phi_grid must be positive when log_phi is True.") + + log_phi_grid = np.log10(phi_grid_arr) + + if fill_value > 0.0: + log_fill_value = np.log10(fill_value) + elif fill_value == 0.0: + log_fill_value = -np.inf + else: + raise ValueError("fill_value must be non-negative.") + + interpolated = np.interp( + absolute_mag_arr, + magnitude_grid_arr, + log_phi_grid, + left=log_fill_value, + right=log_fill_value, + ) + phi = 10.0**interpolated + phi = np.where(np.isfinite(phi), phi, fill_value) + else: + if fill_value < 0.0: + raise ValueError("fill_value must be non-negative.") + + phi = np.interp( + absolute_mag_arr, + magnitude_grid_arr, + phi_grid_arr, + left=fill_value, + right=fill_value, + ) + + return np.asarray(phi, dtype=float) + + +def binned_lf( + absolute_mag: FloatInput, + *, + magnitude_bin_edges: ParameterValue, + phi_bin_values: ParameterValue, + fill_value: float = 0.0, +) -> FloatArray: + r"""Return a piecewise constant binned luminosity function. + + This evaluates a non-parametric luminosity function defined by constant + values inside absolute magnitude bins. If the bin edges are + + .. math:: + + M_0 < M_1 < \cdots < M_N, + + and the supplied bin values are :math:`\phi_j`, then + + .. math:: + + \phi(M) = \phi_j + \quad \mathrm{for} \quad + M_j \le M < M_{j+1}. + + The final edge is treated as the upper boundary of the last bin. Values + outside the supplied bin range are set to ``fill_value``. + + This representation is closest to how many observational luminosity + functions are reported: as measurements in finite magnitude bins. Unlike + ``tabulated_lf``, this function does not smooth or interpolate between + measurements. It preserves the step-like structure of the binned estimate. + + This is useful when the bin values themselves are the model parameters, + when one wants a deliberately non-parametric LF, or when interpolation + between bins would imply more information than the data actually provide. + + Args: + absolute_mag: Absolute magnitude value(s) where the luminosity function + should be evaluated. + magnitude_bin_edges: One-dimensional strictly increasing absolute + magnitude bin edges. This array must have one more element than + ``phi_bin_values``. + phi_bin_values: One-dimensional array of non-negative LF values inside + each magnitude bin. + fill_value: Value returned outside the supplied magnitude bin range. + + Returns: + NumPy array containing the binned luminosity function evaluated at + ``absolute_mag``. + + Raises: + ValueError: If the bin edges or bin values are not one-dimensional, if + the number of edges is not one larger than the number of values, if + the edges are not strictly increasing, if any LF value is negative, + or if ``fill_value`` is not finite and non-negative. + """ + absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") + edges, values = validate_binned_grid( + magnitude_bin_edges, + phi_bin_values, + edges_name="magnitude_bin_edges", + values_name="phi_bin_values", + ) + + if not np.isfinite(fill_value): + raise ValueError("fill_value must be finite.") + + if fill_value < 0.0: + raise ValueError("fill_value must be non-negative.") + + indices = np.searchsorted(edges, absolute_mag_arr, side="right") - 1 + valid = (indices >= 0) & (indices < values.size) + + phi = np.full(absolute_mag_arr.shape, fill_value, dtype=float) + phi[valid] = values[indices[valid]] + + return np.asarray(phi, dtype=float) + + +def redshift_tabulated_lf( + absolute_mag: FloatInput, + redshift: FloatInput, + *, + magnitude_grid: ParameterValue, + redshift_grid: ParameterValue, + phi_grid: ParameterValue, + fill_value: float = 0.0, + log_phi: bool = False, +) -> FloatArray: + r"""Return a luminosity function interpolated over magnitude and redshift. + + This evaluates a redshift-dependent tabulated luminosity function, + + .. math:: + + \phi(M_i, z_k) = \phi_{k i}, + + where ``magnitude_grid`` contains absolute magnitudes :math:`M_i`, + ``redshift_grid`` contains redshifts :math:`z_k`, and ``phi_grid`` has + shape + + .. math:: + + (N_z, N_M). + + For an input pair :math:`(M, z)`, the model performs bilinear interpolation + on the supplied grid. Operationally, it first interpolates in magnitude at + each tabulated redshift and then interpolates those values in redshift. + + If ``log_phi`` is True, the interpolation is performed in + :math:`\log_{10}\phi`, + + .. math:: + + \log_{10}\phi(M, z) + = + \mathrm{interp}_{M,z} + \left[ + \log_{10}\phi(M_i, z_k) + \right], + + and the result is converted back to linear space. + + This model is useful when the LF is measured or simulated in several + redshift slices and one wants a smooth LF surface :math:`\phi(M,z)` without + assuming an analytic redshift evolution law. It is more flexible than an + evolving Schechter or evolving Saunders form because each redshift slice can + have an arbitrary shape. + + Values outside the supplied magnitude or redshift grid are not extrapolated. + They are set to ``fill_value``. + + Args: + absolute_mag: Absolute magnitude value(s) where the luminosity function + should be evaluated. + redshift: Redshift value(s) where the luminosity function should be + evaluated. Must be non-negative. + magnitude_grid: One-dimensional strictly increasing absolute magnitude + grid. + redshift_grid: One-dimensional strictly increasing non-negative + redshift grid. + phi_grid: Two-dimensional array of non-negative LF values with shape + ``(redshift_grid.size, magnitude_grid.size)``. + fill_value: Value returned outside the tabulated magnitude or redshift + range. + log_phi: If True, interpolate in ``log10(phi)`` instead of directly in + ``phi``. + + Returns: + NumPy array containing the redshift-dependent tabulated luminosity + function evaluated at ``absolute_mag`` and ``redshift``. + + Raises: + ValueError: If the magnitude or redshift grid is invalid, if + ``phi_grid`` has the wrong shape, if any LF value is negative, if + any requested redshift is negative, if ``fill_value`` is not finite, + or if ``log_phi`` is True and any tabulated LF value is not + positive. + """ + absolute_mag_arr, redshift_arr = np.broadcast_arrays( + validate_array(absolute_mag, name="absolute_mag"), + validate_array(redshift, name="redshift"), + ) + magnitude_grid_arr, redshift_grid_arr, phi_grid_arr = validate_2d_tabulated_grid( + magnitude_grid, + redshift_grid, + phi_grid, + x_name="magnitude_grid", + y_name="redshift_grid", + values_name="phi_grid", + allow_negative_y=False, + positive_values=log_phi, + ) + + if np.any(redshift_arr < 0.0): + raise ValueError("redshift must be non-negative.") + + if not np.isfinite(fill_value): + raise ValueError("fill_value must be finite.") + + if fill_value < 0.0: + raise ValueError("fill_value must be non-negative.") + + if log_phi: + if np.any(phi_grid_arr <= 0.0): + raise ValueError("phi_grid must be positive when log_phi is True.") + + table = np.log10(phi_grid_arr) + outside_value = np.log10(fill_value) if fill_value > 0.0 else -np.inf + else: + table = phi_grid_arr + outside_value = fill_value + + out = np.full(absolute_mag_arr.shape, outside_value, dtype=float) + + for index in np.ndindex(absolute_mag_arr.shape): + mag = absolute_mag_arr[index] + z = redshift_arr[index] + + outside_mag = mag < magnitude_grid_arr[0] or mag > magnitude_grid_arr[-1] + outside_z = z < redshift_grid_arr[0] or z > redshift_grid_arr[-1] + + if outside_mag or outside_z: + continue + + phi_at_mag_by_redshift = np.array( + [ + np.interp(mag, magnitude_grid_arr, table[z_index]) + for z_index in range(redshift_grid_arr.size) + ], + dtype=float, + ) + out[index] = np.interp( + z, + redshift_grid_arr, + phi_at_mag_by_redshift, + ) + + if log_phi: + out = 10.0**out + out = np.where(np.isfinite(out), out, fill_value) + + return np.asarray(out, dtype=float) + + +def redshift_binned_lf( + absolute_mag: FloatInput, + redshift: FloatInput, + *, + magnitude_bin_edges: ParameterValue, + redshift_bin_edges: ParameterValue, + phi_bin_values: ParameterValue, + fill_value: float = 0.0, +) -> FloatArray: + r"""Return a piecewise constant binned luminosity function in magnitude and redshift. + + This evaluates a two-dimensional binned luminosity function. If the + magnitude bin edges are + + .. math:: + + M_0 < M_1 < \cdots < M_N, + + and the redshift bin edges are + + .. math:: + + z_0 < z_1 < \cdots < z_K, + + then the supplied values define + + .. math:: + + \phi(M,z) = \phi_{k j} + \quad \mathrm{for} \quad + M_j \le M < M_{j+1}, + \quad + z_k \le z < z_{k+1}. + + The ``phi_bin_values`` array must therefore have shape + + .. math:: + + (K, N) + = + (N_z - 1, N_M - 1). + + This model is the redshift-dependent analogue of ``binned_lf``. It is + useful for representing luminosity functions measured independently in + redshift and magnitude bins. It does not assume smoothness in either + direction. That makes it useful for conservative non-parametric analyses, + survey selection studies, or tests where each LF bin is treated as an + independent degree of freedom. + + Values outside the supplied magnitude or redshift bin ranges are set to + ``fill_value``. + + Args: + absolute_mag: Absolute magnitude value(s) where the luminosity function + should be evaluated. + redshift: Redshift value(s) where the luminosity function should be + evaluated. Must be non-negative. + magnitude_bin_edges: One-dimensional strictly increasing absolute + magnitude bin edges. + redshift_bin_edges: One-dimensional strictly increasing non-negative + redshift bin edges. + phi_bin_values: Two-dimensional array of non-negative LF values with + shape ``(redshift_bin_edges.size - 1, + magnitude_bin_edges.size - 1)``. + fill_value: Value returned outside the supplied magnitude or redshift + bin ranges. + + Returns: + NumPy array containing the binned luminosity function evaluated at + ``absolute_mag`` and ``redshift``. + + Raises: + ValueError: If the bin edges are invalid, if ``phi_bin_values`` has the + wrong shape, if any LF value is negative, if any requested redshift + is negative, or if ``fill_value`` is not finite and non-negative. + """ + absolute_mag_arr, redshift_arr = np.broadcast_arrays( + validate_array(absolute_mag, name="absolute_mag"), + validate_array(redshift, name="redshift"), + ) + mag_edges, z_edges, values = validate_2d_binned_grid( + magnitude_bin_edges, + redshift_bin_edges, + phi_bin_values, + x_edges_name="magnitude_bin_edges", + y_edges_name="redshift_bin_edges", + values_name="phi_bin_values", + allow_negative_y_edges=False, + ) + + if np.any(redshift_arr < 0.0): + raise ValueError("redshift must be non-negative.") + + if not np.isfinite(fill_value): + raise ValueError("fill_value must be finite.") + + if fill_value < 0.0: + raise ValueError("fill_value must be non-negative.") + + mag_indices = np.searchsorted(mag_edges, absolute_mag_arr, side="right") - 1 + z_indices = np.searchsorted(z_edges, redshift_arr, side="right") - 1 + + valid = ( + (mag_indices >= 0) + & (mag_indices < mag_edges.size - 1) + & (z_indices >= 0) + & (z_indices < z_edges.size - 1) + ) + + phi = np.full(absolute_mag_arr.shape, fill_value, dtype=float) + phi[valid] = values[z_indices[valid], mag_indices[valid]] + + return np.asarray(phi, dtype=float) + + +def distance_tabulated_lf( + absolute_mag: FloatInput, + comoving_distance: FloatInput, + *, + magnitude_grid: ParameterValue, + distance_grid: ParameterValue, + phi_grid: ParameterValue, + fill_value: float = 0.0, + log_phi: bool = False, +) -> FloatArray: + r"""Return a luminosity function interpolated over magnitude and comoving distance. + + This evaluates a non-parametric luminosity function tabulated on a + two-dimensional grid, + + .. math:: + + \phi(M_i, \chi_k) = \phi_{k i}, + + where :math:`M_i` are absolute magnitude grid points and :math:`\chi_k` + are comoving distance grid points. + + This is analogous to ``redshift_tabulated_lf``, but uses comoving distance + instead of redshift as the second coordinate. It can be useful for + projection calculations, tomographic kernels, or simulation analyses where + the natural radial coordinate is distance rather than redshift. + + Args: + absolute_mag: Absolute magnitude value(s). + comoving_distance: Comoving distance value(s). Must be non-negative. + magnitude_grid: One-dimensional strictly increasing absolute magnitude + grid. + distance_grid: One-dimensional strictly increasing non-negative + comoving distance grid. + phi_grid: Two-dimensional array of non-negative LF values with shape + ``(distance_grid.size, magnitude_grid.size)``. + fill_value: Value returned outside the tabulated magnitude or distance + range. + log_phi: If True, interpolate in ``log10(phi)``. + + Returns: + NumPy array containing the interpolated luminosity function evaluated at + ``absolute_mag`` and ``comoving_distance``. + + Raises: + ValueError: If the grids are invalid, if ``phi_grid`` has the wrong + shape, if LF values are negative, if requested distances are + negative, or if ``log_phi`` is True and any LF value is not + positive. + """ + absolute_mag_arr, distance_arr = np.broadcast_arrays( + validate_array(absolute_mag, name="absolute_mag"), + validate_array(comoving_distance, name="comoving_distance"), + ) + magnitude_grid_arr, distance_grid_arr, phi_grid_arr = validate_2d_tabulated_grid( + magnitude_grid, + distance_grid, + phi_grid, + x_name="magnitude_grid", + y_name="distance_grid", + values_name="phi_grid", + allow_negative_y=False, + positive_values=log_phi, + ) + + if np.any(distance_arr < 0.0): + raise ValueError("comoving_distance must be non-negative.") + + if not np.isfinite(fill_value): + raise ValueError("fill_value must be finite.") + + if fill_value < 0.0: + raise ValueError("fill_value must be non-negative.") + + if log_phi: + if np.any(phi_grid_arr <= 0.0): + raise ValueError("phi_grid must be positive when log_phi is True.") + + table = np.log10(phi_grid_arr) + outside_value = np.log10(fill_value) if fill_value > 0.0 else -np.inf + else: + table = phi_grid_arr + outside_value = fill_value + + out = np.full(absolute_mag_arr.shape, outside_value, dtype=float) + + for index in np.ndindex(absolute_mag_arr.shape): + mag = absolute_mag_arr[index] + distance = distance_arr[index] + + outside_mag = mag < magnitude_grid_arr[0] or mag > magnitude_grid_arr[-1] + outside_distance = ( + distance < distance_grid_arr[0] or distance > distance_grid_arr[-1] + ) + + if outside_mag or outside_distance: + continue + + phi_at_mag_by_distance = np.array( + [ + np.interp(mag, magnitude_grid_arr, table[distance_index]) + for distance_index in range(distance_grid_arr.size) + ], + dtype=float, + ) + out[index] = np.interp( + distance, + distance_grid_arr, + phi_at_mag_by_distance, + ) + + if log_phi: + out = 10.0**out + out = np.where(np.isfinite(out), out, fill_value) + + return np.asarray(out, dtype=float) + + +def distance_binned_lf( + absolute_mag: FloatInput, + comoving_distance: FloatInput, + *, + magnitude_bin_edges: ParameterValue, + distance_bin_edges: ParameterValue, + phi_bin_values: ParameterValue, + fill_value: float = 0.0, +) -> FloatArray: + r"""Return a piecewise constant binned LF in magnitude and comoving distance. + + This evaluates a luminosity function defined in bins of absolute magnitude + and comoving distance, + + .. math:: + + \phi(M,\chi) = \phi_{k j} + + for + + .. math:: + + M_j \le M < M_{j+1}, + \qquad + \chi_k \le \chi < \chi_{k+1}. + + This is useful when an LF is measured or assigned directly in radial + distance shells, for example in projection-space calculations or mock + light-cone analyses. + + Args: + absolute_mag: Absolute magnitude value(s). + comoving_distance: Comoving distance value(s). Must be non-negative. + magnitude_bin_edges: One-dimensional strictly increasing absolute + magnitude bin edges. + distance_bin_edges: One-dimensional strictly increasing non-negative + comoving distance bin edges. + phi_bin_values: Two-dimensional array of non-negative LF values with + shape ``(distance_bin_edges.size - 1, + magnitude_bin_edges.size - 1)``. + fill_value: Value returned outside the supplied magnitude or distance + bin ranges. + + Returns: + NumPy array containing the binned luminosity function evaluated at + ``absolute_mag`` and ``comoving_distance``. + + Raises: + ValueError: If bin edges are invalid, if ``phi_bin_values`` has the + wrong shape, if LF values are negative, or if requested distances + are negative. + """ + absolute_mag_arr, distance_arr = np.broadcast_arrays( + validate_array(absolute_mag, name="absolute_mag"), + validate_array(comoving_distance, name="comoving_distance"), + ) + mag_edges, distance_edges, values = validate_2d_binned_grid( + magnitude_bin_edges, + distance_bin_edges, + phi_bin_values, + x_edges_name="magnitude_bin_edges", + y_edges_name="distance_bin_edges", + values_name="phi_bin_values", + allow_negative_y_edges=False, + ) + + if np.any(distance_arr < 0.0): + raise ValueError("comoving_distance must be non-negative.") + + if not np.isfinite(fill_value): + raise ValueError("fill_value must be finite.") + + if fill_value < 0.0: + raise ValueError("fill_value must be non-negative.") + + mag_indices = np.searchsorted(mag_edges, absolute_mag_arr, side="right") - 1 + distance_indices = np.searchsorted(distance_edges, distance_arr, side="right") - 1 + + valid = ( + (mag_indices >= 0) + & (mag_indices < mag_edges.size - 1) + & (distance_indices >= 0) + & (distance_indices < distance_edges.size - 1) + ) + + phi = np.full(absolute_mag_arr.shape, fill_value, dtype=float) + phi[valid] = values[distance_indices[valid], mag_indices[valid]] + + return np.asarray(phi, dtype=float) diff --git a/src/lfkit/luminosity_functions/models/saunders.py b/src/lfkit/luminosity_functions/models/saunders.py new file mode 100644 index 00000000..cc6cffc0 --- /dev/null +++ b/src/lfkit/luminosity_functions/models/saunders.py @@ -0,0 +1,230 @@ +"""Saunders luminosity function models.""" + +from __future__ import annotations + +import numpy as np + +from lfkit.photometry.luminosities import luminosity_ratio +from lfkit.utils.types import FloatArray, FloatInput, ParameterValue +from lfkit.utils.validators import validate_array + + +__all__ = [ + "saunders_lf", + "evolving_saunders_lf", + "double_saunders_lf", + "generalized_saunders_lf", +] + + +def saunders_lf( + absolute_mag: FloatInput, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, + sigma: ParameterValue, +) -> FloatArray: + r"""Return a Saunders luminosity function in magnitude space.""" + return generalized_saunders_lf( + absolute_mag, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + sigma=sigma, + beta=2.0, + ) + + +def generalized_saunders_lf( + absolute_mag: FloatInput, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, + sigma: ParameterValue, + beta: ParameterValue, +) -> FloatArray: + r"""Return a generalized Saunders luminosity function in magnitude space. + + This computes + + .. math:: + + \phi(M) = + 0.4 \ln(10) \, \phi_\star \, + x^{\alpha} + \exp\left[ + -\left( + \frac{\log_{10}(1 + x)} + {\sqrt{2}\sigma} + \right)^\beta + \right], + + where + + .. math:: + + x = 10^{-0.4(M - M_\star)} = L/L_\star. + + For ``beta = 2``, this reduces to the standard Saunders form. + + Args: + absolute_mag: Absolute magnitude value(s). + phi_star: Non-negative normalization. + m_star: Characteristic magnitude used to define ``x``. + alpha: Faint-end luminosity slope parameter. + sigma: Positive width of the bright-end logarithmic cutoff. + beta: Positive exponent controlling the bright-end cutoff shape. + + Returns: + NumPy array containing the generalized Saunders luminosity function + evaluated at ``absolute_mag``. + + Raises: + ValueError: If ``phi_star`` is negative, ``sigma`` is not positive, + or ``beta`` is not positive. + """ + absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") + phi_star_arr = validate_array(phi_star, name="phi_star") + alpha_arr = validate_array(alpha, name="alpha") + sigma_arr = validate_array(sigma, name="sigma") + beta_arr = validate_array(beta, name="beta") + + if np.any(phi_star_arr < 0.0): + raise ValueError("phi_star must be non-negative.") + + if np.any(sigma_arr <= 0.0): + raise ValueError("sigma must be positive.") + + if np.any(beta_arr <= 0.0): + raise ValueError("beta must be positive.") + + x = luminosity_ratio(absolute_mag_arr, m_star) + + cutoff_argument = np.log10(1.0 + x) / (np.sqrt(2.0) * sigma_arr) + cutoff = np.exp(-(cutoff_argument**beta_arr)) + phi = 0.4 * np.log(10.0) * phi_star_arr * x**alpha_arr * cutoff + + return np.asarray(phi, dtype=float) + + +def evolving_saunders_lf( + absolute_mag: FloatInput, + redshift: FloatInput, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, + sigma: ParameterValue, + p: ParameterValue = 0.0, + q: ParameterValue = 0.0, +) -> FloatArray: + r"""Return an evolving Saunders luminosity function. + + This computes a Saunders luminosity function with redshift evolution in + normalization and characteristic magnitude, + + .. math:: + + \phi_\star(z) = \phi_\star (1 + z)^p, + + and + + .. math:: + + M_\star(z) = M_\star - qz. + + Args: + absolute_mag: Absolute magnitude value(s). + redshift: Redshift value(s). + phi_star: Non-negative normalization at ``z = 0``. + m_star: Characteristic magnitude at ``z = 0``. + alpha: Faint-end luminosity slope parameter. + sigma: Positive width of the bright-end logarithmic cutoff. + p: Density evolution exponent. + q: Linear magnitude evolution coefficient. + + Returns: + NumPy array containing the evolving Saunders luminosity function + evaluated at ``absolute_mag`` and ``redshift``. + + Raises: + ValueError: If ``redshift`` is negative. + """ + redshift_arr = validate_array(redshift, name="redshift") + p_arr = validate_array(p, name="p") + q_arr = validate_array(q, name="q") + + if np.any(redshift_arr < 0.0): + raise ValueError("redshift must be non-negative.") + + phi_star_z = ( + validate_array(phi_star, name="phi_star") * (1.0 + redshift_arr) ** p_arr + ) + m_star_z = validate_array(m_star, name="m_star") - q_arr * redshift_arr + + return saunders_lf( + absolute_mag, + phi_star=phi_star_z, + m_star=m_star_z, + alpha=alpha, + sigma=sigma, + ) + + +def double_saunders_lf( + absolute_mag: FloatInput, + *, + phi_star_1: ParameterValue, + m_star_1: ParameterValue, + alpha_1: ParameterValue, + sigma_1: ParameterValue, + phi_star_2: ParameterValue, + m_star_2: ParameterValue, + alpha_2: ParameterValue, + sigma_2: ParameterValue, +) -> FloatArray: + r"""Return a two-component Saunders luminosity function. + + This computes + + .. math:: + + \phi(M) = \phi_1(M) + \phi_2(M), + + where each component is a Saunders luminosity function with its own + normalization, characteristic magnitude, faint-end slope, and bright-end + cutoff width. + + Args: + absolute_mag: Absolute magnitude value(s). + phi_star_1: Non-negative normalization of the first component. + m_star_1: Characteristic magnitude of the first component. + alpha_1: Faint-end luminosity slope parameter of the first component. + sigma_1: Positive bright-end cutoff width of the first component. + phi_star_2: Non-negative normalization of the second component. + m_star_2: Characteristic magnitude of the second component. + alpha_2: Faint-end luminosity slope parameter of the second component. + sigma_2: Positive bright-end cutoff width of the second component. + + Returns: + NumPy array containing the summed Saunders luminosity function + evaluated at ``absolute_mag``. + """ + phi_1 = saunders_lf( + absolute_mag, + phi_star=phi_star_1, + m_star=m_star_1, + alpha=alpha_1, + sigma=sigma_1, + ) + phi_2 = saunders_lf( + absolute_mag, + phi_star=phi_star_2, + m_star=m_star_2, + alpha=alpha_2, + sigma=sigma_2, + ) + + return np.asarray(phi_1 + phi_2, dtype=float) diff --git a/src/lfkit/luminosity_functions/models/schechter.py b/src/lfkit/luminosity_functions/models/schechter.py index 809db221..8d7ebd47 100644 --- a/src/lfkit/luminosity_functions/models/schechter.py +++ b/src/lfkit/luminosity_functions/models/schechter.py @@ -38,6 +38,9 @@ "double_schechter_from_m", "schechter_cumulative", "schechter_cumulative_evolving", + "modified_schechter", + "truncated_schechter", + "multi_schechter", ] @@ -673,3 +676,122 @@ def schechter_cumulative_evolving( n_total = np.asarray(phi_star * total_gamma, dtype=float) return np.asarray(n_total - n_brighter, dtype=float) + + +def modified_schechter( + absolute_mag: FloatInput, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, + beta: ParameterValue, +) -> FloatArray: + r"""Return a modified Schechter luminosity function in magnitude space.""" + absolute_mag = validate_array(absolute_mag, name="absolute_mag") + phi_star_arr = validate_array(phi_star, name="phi_star") + alpha_arr = validate_array(alpha, name="alpha") + beta_arr = validate_array(beta, name="beta") + + if np.any(phi_star_arr < 0): + raise ValueError("phi_star must be non-negative.") + + if np.any(phi_star_arr == 0): + warnings.warn("phi_star is zero; LF will be identically zero.", stacklevel=2) + + if np.any(beta_arr <= 0.0): + raise ValueError("beta must be positive.") + + x = luminosity_ratio(absolute_mag, m_star) + x = np.clip(x, 1e-300, 1e300) + + prefactor = 0.4 * np.log(10.0) * beta_arr * phi_star_arr + + return np.asarray( + prefactor * x ** (alpha_arr + 1.0) * np.exp(-(x**beta_arr)), + dtype=float, + ) + + +def truncated_schechter( + absolute_mag: FloatInput, + *, + phi_star: ParameterValue, + m_star: ParameterValue, + alpha: ParameterValue, + m_bright: float | None = None, + m_faint: float | None = None, +) -> FloatArray: + r"""Return a Schechter luminosity function truncated outside magnitude limits.""" + absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") + + if (m_bright is not None) and (not np.isfinite(m_bright)): + raise ValueError("m_bright must be finite when provided.") + + if (m_faint is not None) and (not np.isfinite(m_faint)): + raise ValueError("m_faint must be finite when provided.") + + if (m_bright is not None) and (m_faint is not None) and (m_bright >= m_faint): + raise ValueError("m_bright must be less than m_faint.") + + phi = schechter( + absolute_mag_arr, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + ) + + mask = np.ones_like(absolute_mag_arr, dtype=bool) + + if m_bright is not None: + mask &= absolute_mag_arr >= m_bright + + if m_faint is not None: + mask &= absolute_mag_arr <= m_faint + + return np.asarray(np.where(mask, phi, 0.0), dtype=float) + + +def multi_schechter( + absolute_mag: FloatInput, + *, + phi_stars: FloatInput, + m_stars: FloatInput, + alphas: FloatInput, +) -> FloatArray: + r"""Return a sum of Schechter luminosity function components.""" + absolute_mag_arr = validate_array(absolute_mag, name="absolute_mag") + phi_stars_arr = validate_array(phi_stars, name="phi_stars") + m_stars_arr = validate_array(m_stars, name="m_stars") + alphas_arr = validate_array(alphas, name="alphas") + + if not (phi_stars_arr.shape == m_stars_arr.shape == alphas_arr.shape): + raise ValueError("phi_stars, m_stars, and alphas must have matching shapes.") + + if phi_stars_arr.ndim != 1: + raise ValueError("phi_stars, m_stars, and alphas must be 1D arrays.") + + if np.any(phi_stars_arr < 0): + raise ValueError("phi_stars must be non-negative.") + + if np.any(phi_stars_arr == 0): + warnings.warn( + "At least one phi_star is zero; that component is identically zero.", + stacklevel=2, + ) + + total = np.zeros_like(absolute_mag_arr, dtype=float) + + for phi_star, m_star, alpha in zip( + phi_stars_arr, + m_stars_arr, + alphas_arr, + strict=True, + ): + total += schechter( + absolute_mag_arr, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + ) + + return np.asarray(total, dtype=float) diff --git a/src/lfkit/utils/validators.py b/src/lfkit/utils/validators.py index cb4c85e4..b86bc1e6 100644 --- a/src/lfkit/utils/validators.py +++ b/src/lfkit/utils/validators.py @@ -10,6 +10,11 @@ "validate_array", "validate_luminosity_distance", "validate_magnitude_range", + "validate_strictly_increasing_1d", + "validate_tabulated_grid", + "validate_binned_grid", + "validate_2d_tabulated_grid", + "validate_2d_binned_grid", ] @@ -61,3 +66,176 @@ def validate_magnitude_range( if m_faint <= m_bright: raise ValueError("m_faint must be larger than m_bright.") + + +def validate_strictly_increasing_1d( + x: FloatInput, + *, + name: str, + min_size: int = 2, + allow_negative: bool = True, +) -> FloatArray: + """Return a finite strictly increasing one-dimensional float array.""" + arr = validate_array(x, name=name, allow_negative=allow_negative) + + if arr.ndim != 1: + raise ValueError(f"{name} must be one-dimensional.") + + if arr.size < min_size: + raise ValueError(f"{name} must contain at least {min_size} values.") + + if np.any(np.diff(arr) <= 0.0): + raise ValueError(f"{name} must be strictly increasing.") + + return arr + + +def validate_tabulated_grid( + coordinate_grid: FloatInput, + values: FloatInput, + *, + coordinate_name: str, + values_name: str, + allow_negative_coordinate: bool = True, + positive_values: bool = False, +) -> tuple[FloatArray, FloatArray]: + """Return validated one-dimensional tabulated grid coordinates and values.""" + coordinate_arr = validate_strictly_increasing_1d( + coordinate_grid, + name=coordinate_name, + allow_negative=allow_negative_coordinate, + ) + values_arr = validate_array(values, name=values_name) + + if values_arr.ndim != 1: + raise ValueError(f"{values_name} must be one-dimensional.") + + if coordinate_arr.size != values_arr.size: + raise ValueError( + f"{coordinate_name} and {values_name} must have the same length." + ) + + if positive_values: + if np.any(values_arr <= 0.0): + raise ValueError(f"{values_name} must be positive.") + elif np.any(values_arr < 0.0): + raise ValueError(f"{values_name} must be non-negative.") + + return coordinate_arr, values_arr + + +def validate_binned_grid( + bin_edges: FloatInput, + bin_values: FloatInput, + *, + edges_name: str, + values_name: str, + allow_negative_edges: bool = True, + positive_values: bool = False, +) -> tuple[FloatArray, FloatArray]: + """Return validated one-dimensional bin edges and bin values.""" + edges_arr = validate_strictly_increasing_1d( + bin_edges, + name=edges_name, + allow_negative=allow_negative_edges, + ) + values_arr = validate_array(bin_values, name=values_name) + + if values_arr.ndim != 1: + raise ValueError(f"{values_name} must be one-dimensional.") + + if edges_arr.size != values_arr.size + 1: + raise ValueError(f"{edges_name} must have one more value than {values_name}.") + + if positive_values: + if np.any(values_arr <= 0.0): + raise ValueError(f"{values_name} must be positive.") + elif np.any(values_arr < 0.0): + raise ValueError(f"{values_name} must be non-negative.") + + return edges_arr, values_arr + + +def validate_2d_tabulated_grid( + x_grid: FloatInput, + y_grid: FloatInput, + values: FloatInput, + *, + x_name: str, + y_name: str, + values_name: str, + allow_negative_x: bool = True, + allow_negative_y: bool = True, + positive_values: bool = False, +) -> tuple[FloatArray, FloatArray, FloatArray]: + """Return validated two-dimensional tabulated coordinates and values.""" + x_arr = validate_strictly_increasing_1d( + x_grid, + name=x_name, + allow_negative=allow_negative_x, + ) + y_arr = validate_strictly_increasing_1d( + y_grid, + name=y_name, + allow_negative=allow_negative_y, + ) + values_arr = validate_array(values, name=values_name) + + if values_arr.ndim != 2: + raise ValueError(f"{values_name} must be two-dimensional.") + + if values_arr.shape != (y_arr.size, x_arr.size): + raise ValueError( + f"{values_name} must have shape ({y_name}.size, {x_name}.size)." + ) + + if positive_values: + if np.any(values_arr <= 0.0): + raise ValueError(f"{values_name} must be positive.") + elif np.any(values_arr < 0.0): + raise ValueError(f"{values_name} must be non-negative.") + + return x_arr, y_arr, values_arr + + +def validate_2d_binned_grid( + x_bin_edges: FloatInput, + y_bin_edges: FloatInput, + values: FloatInput, + *, + x_edges_name: str, + y_edges_name: str, + values_name: str, + allow_negative_x_edges: bool = True, + allow_negative_y_edges: bool = True, + positive_values: bool = False, +) -> tuple[FloatArray, FloatArray, FloatArray]: + """Return validated two-dimensional bin edges and bin values.""" + x_edges_arr = validate_strictly_increasing_1d( + x_bin_edges, + name=x_edges_name, + allow_negative=allow_negative_x_edges, + ) + y_edges_arr = validate_strictly_increasing_1d( + y_bin_edges, + name=y_edges_name, + allow_negative=allow_negative_y_edges, + ) + values_arr = validate_array(values, name=values_name) + + if values_arr.ndim != 2: + raise ValueError(f"{values_name} must be two-dimensional.") + + if values_arr.shape != (y_edges_arr.size - 1, x_edges_arr.size - 1): + raise ValueError( + f"{values_name} must have shape " + f"({y_edges_name}.size - 1, {x_edges_name}.size - 1)." + ) + + if positive_values: + if np.any(values_arr <= 0.0): + raise ValueError(f"{values_name} must be positive.") + elif np.any(values_arr < 0.0): + raise ValueError(f"{values_name} must be non-negative.") + + return x_edges_arr, y_edges_arr, values_arr diff --git a/tests/test_lumfuncs_models_gamma.py b/tests/test_lumfuncs_models_gamma.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_lumfuncs_models_non_parametric.py b/tests/test_lumfuncs_models_non_parametric.py new file mode 100644 index 00000000..e69de29b diff --git a/tests/test_lumfuncs_models_saunders.py b/tests/test_lumfuncs_models_saunders.py new file mode 100644 index 00000000..e69de29b From ea83772d39ab4dc6a62dcf93e08ddccf766cd26e Mon Sep 17 00:00:00 2001 From: "niko.sarcevic" Date: Sun, 14 Jun 2026 23:32:28 -0400 Subject: [PATCH 2/2] updated test suite for new models --- ...est_api_conditional_luminosity_function.py | 29 +- tests/test_api_luminosity_function.py | 113 +++- tests/test_lumfuncs_conditional_models.py | 5 +- tests/test_lumfuncs_models_gamma.py | 387 +++++++++++++ tests/test_lumfuncs_models_non_parametric.py | 515 ++++++++++++++++++ tests/test_lumfuncs_models_saunders.py | 505 +++++++++++++++++ tests/test_lumfuncs_models_schechter.py | 340 ++++++++++++ tests/test_lumfuncs_registry.py | 26 + tests/test_utils_validators.py | 392 +++++++++++++ 9 files changed, 2302 insertions(+), 10 deletions(-) diff --git a/tests/test_api_conditional_luminosity_function.py b/tests/test_api_conditional_luminosity_function.py index e1282190..9311fd49 100644 --- a/tests/test_api_conditional_luminosity_function.py +++ b/tests/test_api_conditional_luminosity_function.py @@ -7,6 +7,7 @@ from lfkit.api.conditional_luminosity_function import ConditionalLuminosityFunction from lfkit.api.luminosity_function import LuminosityFunction +from lfkit.luminosity_functions.registry import CONDITIONAL_LF_MODELS def test_conditional_schechter_constructor_stores_parameters(): @@ -87,9 +88,9 @@ def test_lognormal_constructor_uses_default_amplitude(): assert lf.meta == {} -def test_modified_schechter_constructor_is_not_registered() -> None: - """Tests that modified Schechter is not a standalone conditional model.""" - assert not hasattr(ConditionalLuminosityFunction, "modified_schechter") +def test_luminosity_cutoff_modifier_is_not_registered() -> None: + """Tests that luminosity cutoff modifiers are not registered as conditional LF models.""" + assert "apply_luminosity_cutoff" not in CONDITIONAL_LF_MODELS def test_two_component_constructor_stores_parameters(): @@ -205,3 +206,25 @@ def test_constructor_rejects_unexpected_parameter(): alpha=-1.1, bad_parameter=123, ) + + +@pytest.mark.parametrize( + "name", + [ + "saunders", + "evolving_saunders", + "double_saunders", + "generalized_saunders", + "gamma", + "generalized_gamma", + "tabulated", + "binned", + "redshift_tabulated", + "redshift_binned", + "distance_tabulated", + "distance_binned", + ], +) +def test_available_models_includes_extended_model_families(name: str) -> None: + """Tests that extended model families are available as conditional models.""" + assert name in ConditionalLuminosityFunction.available_models() diff --git a/tests/test_api_luminosity_function.py b/tests/test_api_luminosity_function.py index a93eff10..494cf2d9 100644 --- a/tests/test_api_luminosity_function.py +++ b/tests/test_api_luminosity_function.py @@ -704,16 +704,82 @@ def _minimal_parameter_payload(model_name, function): "alpha_kwargs": {"alpha": -1.1}, } + if model_name == "tabulated": + return { + "magnitude_grid": np.array([-24.0, -22.0, -20.0, -18.0]), + "phi_grid": np.array([1.0e-5, 5.0e-4, 1.0e-3, 2.0e-4]), + } + + if model_name == "binned": + return { + "magnitude_bin_edges": np.array([-24.0, -22.0, -20.0, -18.0]), + "phi_bin_values": np.array([1.0e-5, 5.0e-4, 1.0e-3]), + } + + if model_name == "redshift_tabulated": + return { + "magnitude_grid": np.array([-24.0, -22.0, -20.0, -18.0]), + "redshift_grid": np.array([0.1, 0.5, 1.0]), + "phi_grid": np.array( + [ + [1.0e-5, 5.0e-4, 1.0e-3, 2.0e-4], + [2.0e-5, 6.0e-4, 1.1e-3, 3.0e-4], + [3.0e-5, 7.0e-4, 1.2e-3, 4.0e-4], + ] + ), + } + + if model_name == "redshift_binned": + return { + "magnitude_bin_edges": np.array([-24.0, -22.0, -20.0, -18.0]), + "redshift_bin_edges": np.array([0.1, 0.5, 1.0, 1.5]), + "phi_bin_values": np.array( + [ + [1.0e-5, 5.0e-4, 1.0e-3], + [2.0e-5, 6.0e-4, 1.1e-3], + [3.0e-5, 7.0e-4, 1.2e-3], + ] + ), + } + + if model_name == "distance_tabulated": + return { + "magnitude_grid": np.array([-24.0, -22.0, -20.0, -18.0]), + "distance_grid": np.array([100.0, 500.0, 1000.0]), + "comoving_distance": 500.0, + "phi_grid": np.array( + [ + [1.0e-5, 5.0e-4, 1.0e-3, 2.0e-4], + [2.0e-5, 6.0e-4, 1.1e-3, 3.0e-4], + [3.0e-5, 7.0e-4, 1.2e-3, 4.0e-4], + ] + ), + } + + if model_name == "distance_binned": + return { + "magnitude_bin_edges": np.array([-24.0, -22.0, -20.0, -18.0]), + "distance_bin_edges": np.array([100.0, 500.0, 1000.0, 1500.0]), + "comoving_distance": 500.0, + "phi_bin_values": np.array( + [ + [1.0e-5, 5.0e-4, 1.0e-3], + [2.0e-5, 6.0e-4, 1.1e-3], + [3.0e-5, 7.0e-4, 1.2e-3], + ] + ), + } + payload = {} for name, parameter in inspect.signature(function).parameters.items(): - if name in {"absolute_mag", "z"}: + if name in {"absolute_mag", "z", "redshift"}: continue if parameter.default is not inspect.Parameter.empty: continue - if name in {"phi_star", "modified_phi_star"}: + if name in {"phi_star", "modified_phi_star", "phi_star_1", "phi_star_2"}: payload[name] = 1.0e-3 elif name == "log_phi_star": payload[name] = -3.0 @@ -723,28 +789,42 @@ def _minimal_parameter_payload(model_name, function): "mean_absolute_mag", "lognormal_mean_absolute_mag", "m_transition", + "m_star_1", + "m_star_2", }: payload[name] = -20.5 elif name in { "alpha", - "beta", "alpha_faint", "alpha_bright", "modified_alpha", + "alpha_1", + "alpha_2", }: payload[name] = -1.1 + elif name == "beta": + payload[name] = 1.0 + elif name == "phi_stars": + payload[name] = np.array([1.0e-3, 5.0e-4]) + elif name == "m_stars": + payload[name] = np.array([-20.5, -19.5]) + elif name == "alphas": + payload[name] = np.array([-1.1, -0.5]) elif name in { "sigma_absolute_mag", "sigma_log_luminosity", "lognormal_sigma_log_luminosity", + "sigma", + "sigma_1", + "sigma_2", }: - payload[name] = 0.2 + payload[name] = 0.7 elif name in {"amplitude", "lognormal_amplitude"}: payload[name] = 1.0 elif name in {"fraction", "modified_luminosity_fraction"}: payload[name] = 0.5 else: - pytest.skip(f"No test default for required parameter {name!r}") + pytest.fail(f"No test default for required parameter {name!r}") return payload @@ -795,3 +875,26 @@ def test_registered_luminosity_function_models_expose_generic_integrals(model_na for name in expected_methods: assert callable(getattr(lf.integrals, name)) + + +def test_available_models_includes_extended_model_families() -> None: + """Tests that available models includes extended model families.""" + models = LuminosityFunction.available_models() + + expected = [ + "gamma", + "generalized_gamma", + "saunders", + "evolving_saunders", + "double_saunders", + "generalized_saunders", + "tabulated", + "binned", + "redshift_tabulated", + "redshift_binned", + "distance_tabulated", + "distance_binned", + ] + + for name in expected: + assert name in models diff --git a/tests/test_lumfuncs_conditional_models.py b/tests/test_lumfuncs_conditional_models.py index a442cf68..39150202 100644 --- a/tests/test_lumfuncs_conditional_models.py +++ b/tests/test_lumfuncs_conditional_models.py @@ -461,8 +461,9 @@ def test_conditional_model_registry_exports_generated_names() -> None: assert "conditional_schechter" in __all__ assert "conditional_double_schechter" in __all__ - assert "conditional_lognormal_lf" in __all__ - assert "conditional_two_component_lf" in __all__ + assert "conditional_lognormal" in __all__ + assert "conditional_two_component" in __all__ + assert "conditional_distance_binned" in __all__ def test_conditional_schechter_accepts_scalar_condition_with_callable_parameter() -> None: diff --git a/tests/test_lumfuncs_models_gamma.py b/tests/test_lumfuncs_models_gamma.py index e69de29b..ef3741df 100644 --- a/tests/test_lumfuncs_models_gamma.py +++ b/tests/test_lumfuncs_models_gamma.py @@ -0,0 +1,387 @@ +"""Unit tests for ``lfkit.luminosity_functions.models.gamma``.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from lfkit.luminosity_functions.models.gamma import ( + gamma_lf, + generalized_gamma_lf, +) +from lfkit.luminosity_functions.models.schechter import schechter +from lfkit.utils.integrators import safe_power10 + + +def test_gamma_lf_returns_positive_finite_values() -> None: + """Tests that gamma_lf returns positive finite values.""" + m = np.linspace(-22.0, -18.0, 10) + + result = gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + ) + + assert result.shape == m.shape + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) + + +def test_gamma_lf_matches_schechter() -> None: + """Tests that gamma_lf matches the standard Schechter form.""" + m = np.array([-22.0, -21.0, -20.0, -19.0]) + + result = gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + ) + expected = schechter( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + ) + + np.testing.assert_allclose(result, expected) + + +def test_gamma_lf_matches_manual_formula() -> None: + """Tests that gamma_lf matches the analytic magnitude-space formula.""" + m = np.array([-22.0, -21.0, -20.0]) + phi_star = 1.0e-3 + m_star = -20.0 + alpha = -1.2 + + result = gamma_lf( + m, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + ) + + x = safe_power10(-0.4 * (m - m_star)) + expected = ( + 0.4 + * np.log(10.0) + * phi_star + * x ** (alpha + 1.0) + * np.exp(-x) + ) + + np.testing.assert_allclose(result, expected) + + +def test_gamma_lf_rejects_negative_phi_star() -> None: + """Tests that gamma_lf rejects negative phi_star.""" + with pytest.raises(ValueError, match="phi_star must be non-negative"): + gamma_lf( + [-20.0], + phi_star=-1.0e-3, + m_star=-20.0, + alpha=-1.0, + ) + + +def test_gamma_lf_rejects_nonfinite_absolute_mag() -> None: + """Tests that gamma_lf rejects non-finite absolute magnitudes.""" + with pytest.raises(ValueError, match="absolute_mag contains NaN or infinite values"): + gamma_lf( + [np.nan], + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + ) + + +def test_gamma_lf_rejects_nonfinite_phi_star() -> None: + """Tests that gamma_lf rejects non-finite phi_star.""" + with pytest.raises(ValueError, match="phi_star contains NaN or infinite values"): + gamma_lf( + [-20.0], + phi_star=np.inf, + m_star=-20.0, + alpha=-1.0, + ) + + +def test_gamma_lf_rejects_nonfinite_alpha() -> None: + """Tests that gamma_lf rejects non-finite alpha.""" + with pytest.raises(ValueError, match="alpha contains NaN or infinite values"): + gamma_lf( + [-20.0], + phi_star=1.0e-3, + m_star=-20.0, + alpha=np.nan, + ) + + +def test_gamma_lf_accepts_array_parameters() -> None: + """Tests that gamma_lf supports array-valued parameters.""" + m = np.array([-22.0, -21.0, -20.0]) + + result = gamma_lf( + m, + phi_star=np.array([1.0e-3, 2.0e-3, 3.0e-3]), + m_star=-20.0, + alpha=np.array([-1.2, -1.0, -0.8]), + ) + + assert result.shape == m.shape + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) + + +def test_gamma_lf_broadcasts_scalar_parameters() -> None: + """Tests that gamma_lf broadcasts scalar parameters over magnitudes.""" + m = np.array([-22.0, -21.0, -20.0]) + + result = gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + ) + + assert result.shape == m.shape + + +def test_gamma_lf_extreme_magnitudes_remain_finite() -> None: + """Tests that gamma_lf remains finite for extreme magnitudes.""" + m = np.array([-1000.0, 1000.0]) + + result = gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + ) + + assert result.shape == m.shape + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) + + +def test_generalized_gamma_lf_returns_positive_finite_values() -> None: + """Tests that generalized_gamma_lf returns positive finite values.""" + m = np.linspace(-22.0, -18.0, 10) + + result = generalized_gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + beta=0.8, + ) + + assert result.shape == m.shape + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) + + +def test_generalized_gamma_lf_matches_gamma_when_beta_is_one() -> None: + """Tests that generalized_gamma_lf reduces to gamma_lf when beta=1.""" + m = np.array([-22.0, -21.0, -20.0, -19.0]) + + result = generalized_gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + beta=1.0, + ) + expected = gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + ) + + np.testing.assert_allclose(result, expected) + + +def test_generalized_gamma_lf_matches_manual_formula() -> None: + """Tests that generalized_gamma_lf matches the analytic formula.""" + m = np.array([-22.0, -21.0, -20.0]) + phi_star = 1.0e-3 + m_star = -20.0 + alpha = -1.2 + beta = 0.8 + + result = generalized_gamma_lf( + m, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + beta=beta, + ) + + x = safe_power10(-0.4 * (m - m_star)) + expected = ( + 0.4 + * np.log(10.0) + * beta + * phi_star + * x ** (alpha + 1.0) + * np.exp(-(x**beta)) + ) + + np.testing.assert_allclose(result, expected) + + +def test_generalized_gamma_lf_rejects_negative_phi_star() -> None: + """Tests that generalized_gamma_lf rejects negative phi_star.""" + with pytest.raises(ValueError, match="phi_star must be non-negative"): + generalized_gamma_lf( + [-20.0], + phi_star=-1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=1.0, + ) + + +def test_generalized_gamma_lf_rejects_zero_beta() -> None: + """Tests that generalized_gamma_lf rejects zero beta.""" + with pytest.raises(ValueError, match="beta must be positive"): + generalized_gamma_lf( + [-20.0], + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=0.0, + ) + + +def test_generalized_gamma_lf_rejects_negative_beta() -> None: + """Tests that generalized_gamma_lf rejects negative beta.""" + with pytest.raises(ValueError, match="beta must be positive"): + generalized_gamma_lf( + [-20.0], + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=-0.5, + ) + + +def test_generalized_gamma_lf_rejects_nonfinite_beta() -> None: + """Tests that generalized_gamma_lf rejects non-finite beta.""" + with pytest.raises(ValueError, match="beta contains NaN or infinite values"): + generalized_gamma_lf( + [-20.0], + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=np.inf, + ) + + +def test_generalized_gamma_lf_rejects_nonfinite_absolute_mag() -> None: + """Tests that generalized_gamma_lf rejects non-finite absolute magnitudes.""" + with pytest.raises(ValueError, match="absolute_mag contains NaN or infinite values"): + generalized_gamma_lf( + [np.inf], + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=1.0, + ) + + +def test_generalized_gamma_lf_rejects_nonfinite_phi_star() -> None: + """Tests that generalized_gamma_lf rejects non-finite phi_star.""" + with pytest.raises(ValueError, match="phi_star contains NaN or infinite values"): + generalized_gamma_lf( + [-20.0], + phi_star=np.nan, + m_star=-20.0, + alpha=-1.0, + beta=1.0, + ) + + +def test_generalized_gamma_lf_rejects_nonfinite_alpha() -> None: + """Tests that generalized_gamma_lf rejects non-finite alpha.""" + with pytest.raises(ValueError, match="alpha contains NaN or infinite values"): + generalized_gamma_lf( + [-20.0], + phi_star=1.0e-3, + m_star=-20.0, + alpha=np.inf, + beta=1.0, + ) + + +def test_generalized_gamma_lf_accepts_array_parameters() -> None: + """Tests that generalized_gamma_lf supports array-valued parameters.""" + m = np.array([-22.0, -21.0, -20.0]) + + result = generalized_gamma_lf( + m, + phi_star=np.array([1.0e-3, 2.0e-3, 3.0e-3]), + m_star=-20.0, + alpha=np.array([-1.2, -1.0, -0.8]), + beta=np.array([0.8, 1.0, 1.2]), + ) + + assert result.shape == m.shape + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) + + +def test_generalized_gamma_lf_broadcasts_scalar_parameters() -> None: + """Tests that generalized_gamma_lf broadcasts scalar parameters over magnitudes.""" + m = np.array([-22.0, -21.0, -20.0]) + + result = generalized_gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.2, + beta=0.8, + ) + + assert result.shape == m.shape + + +def test_generalized_gamma_lf_extreme_magnitudes_remain_finite() -> None: + """Tests that generalized_gamma_lf remains finite for extreme magnitudes.""" + m = np.array([-1000.0, 1000.0]) + + result = generalized_gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=1.0, + ) + + assert result.shape == m.shape + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) + + +def test_generalized_gamma_lf_beta_changes_bright_end_shape() -> None: + """Tests that beta changes the bright-end cutoff shape.""" + m = np.array([-24.0]) + + shallow = generalized_gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=0.5, + ) + steep = generalized_gamma_lf( + m, + phi_star=1.0e-3, + m_star=-20.0, + alpha=-1.0, + beta=2.0, + ) + + assert shallow[0] > steep[0] diff --git a/tests/test_lumfuncs_models_non_parametric.py b/tests/test_lumfuncs_models_non_parametric.py index e69de29b..c147807a 100644 --- a/tests/test_lumfuncs_models_non_parametric.py +++ b/tests/test_lumfuncs_models_non_parametric.py @@ -0,0 +1,515 @@ +"""Unit tests for ``lfkit.luminosity_functions.models.non_parametric``.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from lfkit.luminosity_functions.models.non_parametric import ( + binned_lf, + distance_binned_lf, + distance_tabulated_lf, + redshift_binned_lf, + redshift_tabulated_lf, + tabulated_lf, +) + + +MAG_GRID = np.array([-24.0, -22.0, -20.0, -18.0]) +PHI_GRID = np.array([1.0e-5, 5.0e-4, 1.0e-3, 2.0e-4]) + +MAG_EDGES = np.array([-24.0, -22.0, -20.0, -18.0]) +PHI_BINS = np.array([1.0e-5, 5.0e-4, 1.0e-3]) + +Z_GRID = np.array([0.1, 0.5, 1.0]) +Z_EDGES = np.array([0.1, 0.5, 1.0, 1.5]) + +DIST_GRID = np.array([100.0, 500.0, 1000.0]) +DIST_EDGES = np.array([100.0, 500.0, 1000.0, 1500.0]) + +PHI_2D_GRID = np.array( + [ + [1.0e-5, 5.0e-4, 1.0e-3, 2.0e-4], + [2.0e-5, 6.0e-4, 1.1e-3, 3.0e-4], + [3.0e-5, 7.0e-4, 1.2e-3, 4.0e-4], + ] +) + +PHI_2D_BINS = np.array( + [ + [1.0e-5, 5.0e-4, 1.0e-3], + [2.0e-5, 6.0e-4, 1.1e-3], + [3.0e-5, 7.0e-4, 1.2e-3], + ] +) + + +def test_tabulated_lf_interpolates_linearly() -> None: + """Tests that tabulated_lf linearly interpolates between grid values.""" + result = tabulated_lf( + -21.0, + magnitude_grid=MAG_GRID, + phi_grid=PHI_GRID, + ) + + assert np.asarray(result).shape == () + assert result == pytest.approx(7.5e-4) + + +def test_tabulated_lf_accepts_array_input() -> None: + """Tests that tabulated_lf preserves array input shape.""" + absolute_mag = np.array([-23.0, -21.0, -19.0]) + + result = tabulated_lf( + absolute_mag, + magnitude_grid=MAG_GRID, + phi_grid=PHI_GRID, + ) + + assert result.shape == absolute_mag.shape + np.testing.assert_allclose(result, np.array([2.55e-4, 7.5e-4, 6.0e-4])) + + +def test_tabulated_lf_uses_fill_value_outside_grid() -> None: + """Tests that tabulated_lf returns fill_value outside the magnitude grid.""" + result = tabulated_lf( + np.array([-25.0, -21.0, -17.0]), + magnitude_grid=MAG_GRID, + phi_grid=PHI_GRID, + fill_value=9.0, + ) + + np.testing.assert_allclose(result, np.array([9.0, 7.5e-4, 9.0])) + + +def test_tabulated_lf_interpolates_in_log_phi() -> None: + """Tests that tabulated_lf can interpolate in log10(phi).""" + result = tabulated_lf( + -21.0, + magnitude_grid=np.array([-22.0, -20.0]), + phi_grid=np.array([1.0e-4, 1.0e-2]), + log_phi=True, + ) + + assert result == pytest.approx(1.0e-3) + + +def test_tabulated_lf_log_phi_zero_fill_value_outside_grid() -> None: + """Tests that log interpolation supports zero fill_value.""" + result = tabulated_lf( + np.array([-25.0, -21.0]), + magnitude_grid=np.array([-22.0, -20.0]), + phi_grid=np.array([1.0e-4, 1.0e-2]), + fill_value=0.0, + log_phi=True, + ) + + np.testing.assert_allclose(result, np.array([0.0, 1.0e-3])) + + +def test_tabulated_lf_rejects_negative_fill_value() -> None: + """Tests that tabulated_lf rejects negative fill_value.""" + with pytest.raises(ValueError, match="fill_value must be non-negative"): + tabulated_lf( + -21.0, + magnitude_grid=MAG_GRID, + phi_grid=PHI_GRID, + fill_value=-1.0, + ) + + +def test_tabulated_lf_rejects_nonfinite_fill_value() -> None: + """Tests that tabulated_lf rejects non-finite fill_value.""" + with pytest.raises(ValueError, match="fill_value must be finite"): + tabulated_lf( + -21.0, + magnitude_grid=MAG_GRID, + phi_grid=PHI_GRID, + fill_value=np.inf, + ) + + +def test_tabulated_lf_rejects_nonpositive_values_for_log_phi() -> None: + """Tests that log interpolation requires positive LF values.""" + with pytest.raises(ValueError, match="phi_grid must be positive"): + tabulated_lf( + -21.0, + magnitude_grid=MAG_GRID, + phi_grid=np.array([1.0e-5, 0.0, 1.0e-3, 2.0e-4]), + log_phi=True, + ) + + +def test_tabulated_lf_rejects_mismatched_grid_lengths() -> None: + """Tests that tabulated_lf rejects mismatched grid and value lengths.""" + with pytest.raises(ValueError, match="magnitude_grid and phi_grid"): + tabulated_lf( + -21.0, + magnitude_grid=MAG_GRID, + phi_grid=np.array([1.0, 2.0]), + ) + + +def test_binned_lf_returns_piecewise_constant_values() -> None: + """Tests that binned_lf returns the value assigned to each bin.""" + result = binned_lf( + np.array([-23.0, -21.0, -19.0]), + magnitude_bin_edges=MAG_EDGES, + phi_bin_values=PHI_BINS, + ) + + np.testing.assert_allclose(result, PHI_BINS) + + +def test_binned_lf_uses_fill_value_outside_bins() -> None: + """Tests that binned_lf returns fill_value outside the bin range.""" + result = binned_lf( + np.array([-25.0, -23.0, -17.0]), + magnitude_bin_edges=MAG_EDGES, + phi_bin_values=PHI_BINS, + fill_value=7.0, + ) + + np.testing.assert_allclose(result, np.array([7.0, 1.0e-5, 7.0])) + + +def test_binned_lf_treats_upper_edge_as_outside() -> None: + """Tests that the final bin edge is treated as outside the binned range.""" + result = binned_lf( + np.array([-24.0, -18.0]), + magnitude_bin_edges=MAG_EDGES, + phi_bin_values=PHI_BINS, + fill_value=9.0, + ) + + np.testing.assert_allclose(result, np.array([1.0e-5, 9.0])) + + +def test_binned_lf_rejects_negative_fill_value() -> None: + """Tests that binned_lf rejects negative fill_value.""" + with pytest.raises(ValueError, match="fill_value must be non-negative"): + binned_lf( + -21.0, + magnitude_bin_edges=MAG_EDGES, + phi_bin_values=PHI_BINS, + fill_value=-1.0, + ) + + +def test_binned_lf_rejects_nonfinite_fill_value() -> None: + """Tests that binned_lf rejects non-finite fill_value.""" + with pytest.raises(ValueError, match="fill_value must be finite"): + binned_lf( + -21.0, + magnitude_bin_edges=MAG_EDGES, + phi_bin_values=PHI_BINS, + fill_value=np.nan, + ) + + +def test_binned_lf_rejects_wrong_number_of_bin_values() -> None: + """Tests that binned_lf requires one fewer value than bin edges.""" + with pytest.raises(ValueError, match="magnitude_bin_edges must have one more"): + binned_lf( + -21.0, + magnitude_bin_edges=MAG_EDGES, + phi_bin_values=np.array([1.0, 2.0]), + ) + + +def test_redshift_tabulated_lf_interpolates_in_magnitude_and_redshift() -> None: + """Tests bilinear interpolation for redshift_tabulated_lf.""" + result = redshift_tabulated_lf( + -21.0, + 0.3, + magnitude_grid=MAG_GRID, + redshift_grid=Z_GRID, + phi_grid=PHI_2D_GRID, + ) + + assert np.asarray(result).shape == () + assert result == pytest.approx(8.0e-4) + + +def test_redshift_tabulated_lf_accepts_broadcastable_inputs() -> None: + """Tests that redshift_tabulated_lf broadcasts magnitude and redshift.""" + result = redshift_tabulated_lf( + np.array([-21.0, -21.0]), + 0.3, + magnitude_grid=MAG_GRID, + redshift_grid=Z_GRID, + phi_grid=PHI_2D_GRID, + ) + + assert result.shape == (2,) + np.testing.assert_allclose(result, np.array([8.0e-4, 8.0e-4])) + + +def test_redshift_tabulated_lf_uses_fill_value_outside_grid() -> None: + """Tests that redshift_tabulated_lf fills outside magnitude or redshift grid.""" + result = redshift_tabulated_lf( + np.array([-25.0, -21.0, -21.0]), + np.array([0.3, 2.0, 0.3]), + magnitude_grid=MAG_GRID, + redshift_grid=Z_GRID, + phi_grid=PHI_2D_GRID, + fill_value=4.0, + ) + + np.testing.assert_allclose(result, np.array([4.0, 4.0, 8.0e-4])) + + +def test_redshift_tabulated_lf_interpolates_in_log_phi() -> None: + """Tests that redshift_tabulated_lf supports log interpolation.""" + result = redshift_tabulated_lf( + -21.0, + 0.3, + magnitude_grid=np.array([-22.0, -20.0]), + redshift_grid=np.array([0.1, 0.5]), + phi_grid=np.array([[1.0e-4, 1.0e-2], [1.0e-3, 1.0e-1]]), + log_phi=True, + ) + + assert result == pytest.approx(np.sqrt(1.0e-3 * 1.0e-2)) + + +def test_redshift_tabulated_lf_rejects_negative_redshift() -> None: + """Tests that redshift_tabulated_lf rejects negative requested redshift.""" + with pytest.raises(ValueError, match="redshift must be non-negative"): + redshift_tabulated_lf( + -21.0, + -0.1, + magnitude_grid=MAG_GRID, + redshift_grid=Z_GRID, + phi_grid=PHI_2D_GRID, + ) + + +def test_redshift_tabulated_lf_rejects_negative_redshift_grid() -> None: + """Tests that redshift_tabulated_lf rejects negative redshift grid values.""" + with pytest.raises(ValueError, match="redshift_grid contains negative values"): + redshift_tabulated_lf( + -21.0, + 0.3, + magnitude_grid=MAG_GRID, + redshift_grid=np.array([-0.1, 0.5, 1.0]), + phi_grid=PHI_2D_GRID, + ) + + +def test_redshift_tabulated_lf_rejects_wrong_phi_grid_shape() -> None: + """Tests that redshift_tabulated_lf validates phi_grid shape.""" + with pytest.raises(ValueError, match=r"phi_grid must have shape"): + redshift_tabulated_lf( + -21.0, + 0.3, + magnitude_grid=MAG_GRID, + redshift_grid=Z_GRID, + phi_grid=np.ones((2, 4)), + ) + + +def test_redshift_binned_lf_returns_piecewise_constant_values() -> None: + """Tests that redshift_binned_lf returns values from magnitude-redshift bins.""" + result = redshift_binned_lf( + np.array([-23.0, -21.0, -19.0]), + np.array([0.2, 0.7, 1.2]), + magnitude_bin_edges=MAG_EDGES, + redshift_bin_edges=Z_EDGES, + phi_bin_values=PHI_2D_BINS, + ) + + np.testing.assert_allclose(result, np.array([1.0e-5, 6.0e-4, 1.2e-3])) + + +def test_redshift_binned_lf_uses_fill_value_outside_bins() -> None: + """Tests that redshift_binned_lf fills outside magnitude or redshift bins.""" + result = redshift_binned_lf( + np.array([-25.0, -21.0, -21.0]), + np.array([0.2, 2.0, 0.7]), + magnitude_bin_edges=MAG_EDGES, + redshift_bin_edges=Z_EDGES, + phi_bin_values=PHI_2D_BINS, + fill_value=8.0, + ) + + np.testing.assert_allclose(result, np.array([8.0, 8.0, 6.0e-4])) + + +def test_redshift_binned_lf_rejects_negative_redshift() -> None: + """Tests that redshift_binned_lf rejects negative requested redshift.""" + with pytest.raises(ValueError, match="redshift must be non-negative"): + redshift_binned_lf( + -21.0, + -0.1, + magnitude_bin_edges=MAG_EDGES, + redshift_bin_edges=Z_EDGES, + phi_bin_values=PHI_2D_BINS, + ) + + +def test_redshift_binned_lf_rejects_wrong_phi_bin_shape() -> None: + """Tests that redshift_binned_lf validates phi_bin_values shape.""" + with pytest.raises(ValueError, match=r"phi_bin_values must have shape"): + redshift_binned_lf( + -21.0, + 0.7, + magnitude_bin_edges=MAG_EDGES, + redshift_bin_edges=Z_EDGES, + phi_bin_values=np.ones((2, 3)), + ) + + +def test_distance_tabulated_lf_interpolates_in_magnitude_and_distance() -> None: + """Tests bilinear interpolation for distance_tabulated_lf.""" + result = distance_tabulated_lf( + -21.0, + 300.0, + magnitude_grid=MAG_GRID, + distance_grid=DIST_GRID, + phi_grid=PHI_2D_GRID, + ) + + assert np.asarray(result).shape == () + assert result == pytest.approx(8.0e-4) + + +def test_distance_tabulated_lf_accepts_broadcastable_inputs() -> None: + """Tests that distance_tabulated_lf broadcasts magnitude and distance.""" + result = distance_tabulated_lf( + np.array([-21.0, -21.0]), + 300.0, + magnitude_grid=MAG_GRID, + distance_grid=DIST_GRID, + phi_grid=PHI_2D_GRID, + ) + + assert result.shape == (2,) + np.testing.assert_allclose(result, np.array([8.0e-4, 8.0e-4])) + + +def test_distance_tabulated_lf_uses_fill_value_outside_grid() -> None: + """Tests that distance_tabulated_lf fills outside magnitude or distance grid.""" + result = distance_tabulated_lf( + np.array([-25.0, -21.0, -21.0]), + np.array([300.0, 2000.0, 300.0]), + magnitude_grid=MAG_GRID, + distance_grid=DIST_GRID, + phi_grid=PHI_2D_GRID, + fill_value=4.0, + ) + + np.testing.assert_allclose(result, np.array([4.0, 4.0, 8.0e-4])) + + +def test_distance_tabulated_lf_interpolates_in_log_phi() -> None: + """Tests that distance_tabulated_lf supports log interpolation.""" + result = distance_tabulated_lf( + -21.0, + 300.0, + magnitude_grid=np.array([-22.0, -20.0]), + distance_grid=np.array([100.0, 500.0]), + phi_grid=np.array([[1.0e-4, 1.0e-2], [1.0e-3, 1.0e-1]]), + log_phi=True, + ) + + assert result == pytest.approx(np.sqrt(1.0e-3 * 1.0e-2)) + + +def test_distance_tabulated_lf_rejects_negative_distance() -> None: + """Tests that distance_tabulated_lf rejects negative requested distance.""" + with pytest.raises(ValueError, match="comoving_distance must be non-negative"): + distance_tabulated_lf( + -21.0, + -1.0, + magnitude_grid=MAG_GRID, + distance_grid=DIST_GRID, + phi_grid=PHI_2D_GRID, + ) + + +def test_distance_tabulated_lf_rejects_negative_distance_grid() -> None: + """Tests that distance_tabulated_lf rejects negative distance grid values.""" + with pytest.raises(ValueError, match="distance_grid contains negative values"): + distance_tabulated_lf( + -21.0, + 300.0, + magnitude_grid=MAG_GRID, + distance_grid=np.array([-100.0, 500.0, 1000.0]), + phi_grid=PHI_2D_GRID, + ) + + +def test_distance_tabulated_lf_rejects_wrong_phi_grid_shape() -> None: + """Tests that distance_tabulated_lf validates phi_grid shape.""" + with pytest.raises(ValueError, match=r"phi_grid must have shape"): + distance_tabulated_lf( + -21.0, + 300.0, + magnitude_grid=MAG_GRID, + distance_grid=DIST_GRID, + phi_grid=np.ones((2, 4)), + ) + + +def test_distance_binned_lf_returns_piecewise_constant_values() -> None: + """Tests that distance_binned_lf returns values from magnitude-distance bins.""" + result = distance_binned_lf( + np.array([-23.0, -21.0, -19.0]), + np.array([200.0, 700.0, 1200.0]), + magnitude_bin_edges=MAG_EDGES, + distance_bin_edges=DIST_EDGES, + phi_bin_values=PHI_2D_BINS, + ) + + np.testing.assert_allclose(result, np.array([1.0e-5, 6.0e-4, 1.2e-3])) + + +def test_distance_binned_lf_uses_fill_value_outside_bins() -> None: + """Tests that distance_binned_lf fills outside magnitude or distance bins.""" + result = distance_binned_lf( + np.array([-25.0, -21.0, -21.0]), + np.array([200.0, 2000.0, 700.0]), + magnitude_bin_edges=MAG_EDGES, + distance_bin_edges=DIST_EDGES, + phi_bin_values=PHI_2D_BINS, + fill_value=8.0, + ) + + np.testing.assert_allclose(result, np.array([8.0, 8.0, 6.0e-4])) + + +def test_distance_binned_lf_rejects_negative_distance() -> None: + """Tests that distance_binned_lf rejects negative requested distance.""" + with pytest.raises(ValueError, match="comoving_distance must be non-negative"): + distance_binned_lf( + -21.0, + -1.0, + magnitude_bin_edges=MAG_EDGES, + distance_bin_edges=DIST_EDGES, + phi_bin_values=PHI_2D_BINS, + ) + + +def test_distance_binned_lf_rejects_negative_distance_edges() -> None: + """Tests that distance_binned_lf rejects negative distance bin edges.""" + with pytest.raises(ValueError, match="distance_bin_edges contains negative values"): + distance_binned_lf( + -21.0, + 300.0, + magnitude_bin_edges=MAG_EDGES, + distance_bin_edges=np.array([-100.0, 500.0, 1000.0, 1500.0]), + phi_bin_values=PHI_2D_BINS, + ) + + +def test_distance_binned_lf_rejects_wrong_phi_bin_shape() -> None: + """Tests that distance_binned_lf validates phi_bin_values shape.""" + with pytest.raises(ValueError, match=r"phi_bin_values must have shape"): + distance_binned_lf( + -21.0, + 700.0, + magnitude_bin_edges=MAG_EDGES, + distance_bin_edges=DIST_EDGES, + phi_bin_values=np.ones((2, 3)), + ) diff --git a/tests/test_lumfuncs_models_saunders.py b/tests/test_lumfuncs_models_saunders.py index e69de29b..ac818256 100644 --- a/tests/test_lumfuncs_models_saunders.py +++ b/tests/test_lumfuncs_models_saunders.py @@ -0,0 +1,505 @@ +"""Unit tests for ``lfkit.luminosity_functions.models.saunders``.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from lfkit.luminosity_functions.models.saunders import ( + double_saunders_lf, + evolving_saunders_lf, + generalized_saunders_lf, + saunders_lf, +) +from lfkit.photometry.luminosities import luminosity_ratio + + +def test_saunders_lf_matches_manual_formula() -> None: + """Tests that saunders_lf matches the analytic formula.""" + absolute_mag = np.array([-23.0, -22.0, -21.0]) + phi_star = 0.01 + m_star = -21.0 + alpha = -0.3 + sigma = 0.7 + + result = saunders_lf( + absolute_mag, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + sigma=sigma, + ) + + x = luminosity_ratio(absolute_mag, m_star) + cutoff = np.exp(-(np.log10(1.0 + x) ** 2) / (2.0 * sigma**2)) + expected = 0.4 * np.log(10.0) * phi_star * x**alpha * cutoff + + np.testing.assert_allclose(result, expected) + + +def test_saunders_lf_accepts_scalar_input() -> None: + """Tests that saunders_lf accepts scalar magnitude input.""" + result = saunders_lf(-21.0, phi_star=0.01, m_star=-21.0, alpha=-0.3, sigma=0.7) + + assert np.shape(result) == () + assert np.isfinite(result) + + +def test_saunders_lf_preserves_array_shape() -> None: + """Tests that saunders_lf preserves input shape.""" + absolute_mag = np.array([[-23.0, -22.0], [-21.0, -20.0]]) + + result = saunders_lf( + absolute_mag, + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + ) + + assert result.shape == absolute_mag.shape + + +def test_saunders_lf_zero_phi_star_returns_zero() -> None: + """Tests that zero phi_star returns zero values.""" + result = saunders_lf( + np.array([-23.0, -22.0, -21.0]), + phi_star=0.0, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + ) + + np.testing.assert_allclose(result, np.zeros(3)) + + +def test_saunders_lf_rejects_negative_phi_star() -> None: + """Tests that negative phi_star is rejected.""" + with pytest.raises(ValueError, match="phi_star must be non-negative"): + saunders_lf( + np.array([-23.0, -22.0]), + phi_star=-0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + ) + + +def test_saunders_lf_rejects_nonpositive_sigma() -> None: + """Tests that non-positive sigma is rejected.""" + with pytest.raises(ValueError, match="sigma must be positive"): + saunders_lf( + np.array([-23.0, -22.0]), + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.0, + ) + + +def test_generalized_saunders_lf_matches_manual_formula() -> None: + """Tests that generalized_saunders_lf matches the analytic formula.""" + absolute_mag = np.array([-23.0, -22.0, -21.0]) + phi_star = 0.01 + m_star = -21.0 + alpha = -0.3 + sigma = 0.7 + beta = 1.5 + + result = generalized_saunders_lf( + absolute_mag, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + sigma=sigma, + beta=beta, + ) + + x = luminosity_ratio(absolute_mag, m_star) + cutoff_argument = np.log10(1.0 + x) / (np.sqrt(2.0) * sigma) + cutoff = np.exp(-(cutoff_argument**beta)) + expected = 0.4 * np.log(10.0) * phi_star * x**alpha * cutoff + + np.testing.assert_allclose(result, expected) + + +def test_generalized_saunders_lf_beta_two_matches_saunders_lf() -> None: + """Tests that beta=2 gives the standard Saunders model.""" + absolute_mag = np.array([-23.0, -22.0, -21.0]) + + result = generalized_saunders_lf( + absolute_mag, + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + beta=2.0, + ) + + expected = saunders_lf( + absolute_mag, + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + ) + + np.testing.assert_allclose(result, expected) + + +def test_generalized_saunders_lf_rejects_nonpositive_beta() -> None: + """Tests that non-positive beta is rejected.""" + with pytest.raises(ValueError, match="beta must be positive"): + generalized_saunders_lf( + np.array([-23.0, -22.0]), + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + beta=0.0, + ) + + +def test_evolving_saunders_lf_matches_saunders_lf_at_zero_evolution() -> None: + """Tests that zero evolution matches saunders_lf.""" + absolute_mag = np.array([-23.0, -22.0, -21.0]) + redshift = np.array([0.1, 0.2, 0.3]) + + result = evolving_saunders_lf( + absolute_mag, + redshift, + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + p=0.0, + q=0.0, + ) + + expected = saunders_lf( + absolute_mag, + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + ) + + np.testing.assert_allclose(result, expected) + + +def test_evolving_saunders_lf_applies_density_and_magnitude_evolution() -> None: + """Tests that evolving_saunders_lf applies p and q evolution.""" + absolute_mag = np.array([-23.0, -22.0, -21.0]) + redshift = np.array([0.1, 0.2, 0.3]) + phi_star = 0.01 + m_star = -21.0 + p = 2.0 + q = 1.0 + + result = evolving_saunders_lf( + absolute_mag, + redshift, + phi_star=phi_star, + m_star=m_star, + alpha=-0.3, + sigma=0.7, + p=p, + q=q, + ) + + expected = saunders_lf( + absolute_mag, + phi_star=phi_star * (1.0 + redshift) ** p, + m_star=m_star - q * redshift, + alpha=-0.3, + sigma=0.7, + ) + + np.testing.assert_allclose(result, expected) + + +def test_evolving_saunders_lf_rejects_negative_redshift() -> None: + """Tests that negative redshift is rejected.""" + with pytest.raises(ValueError, match="redshift must be non-negative"): + evolving_saunders_lf( + np.array([-23.0, -22.0]), + np.array([0.1, -0.2]), + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + ) + + +def test_double_saunders_lf_matches_sum_of_components() -> None: + """Tests that double_saunders_lf matches the sum of two Saunders models.""" + absolute_mag = np.array([-23.0, -22.0, -21.0]) + + result = double_saunders_lf( + absolute_mag, + phi_star_1=0.01, + m_star_1=-21.0, + alpha_1=-0.3, + sigma_1=0.7, + phi_star_2=0.005, + m_star_2=-20.5, + alpha_2=-0.8, + sigma_2=0.5, + ) + + expected = saunders_lf( + absolute_mag, + phi_star=0.01, + m_star=-21.0, + alpha=-0.3, + sigma=0.7, + ) + saunders_lf( + absolute_mag, + phi_star=0.005, + m_star=-20.5, + alpha=-0.8, + sigma=0.5, + ) + + np.testing.assert_allclose(result, expected) + + +def test_double_saunders_lf_preserves_array_shape() -> None: + """Tests that double_saunders_lf preserves input shape.""" + absolute_mag = np.array([[-23.0, -22.0], [-21.0, -20.0]]) + + result = double_saunders_lf( + absolute_mag, + phi_star_1=0.01, + m_star_1=-21.0, + alpha_1=-0.3, + sigma_1=0.7, + phi_star_2=0.005, + m_star_2=-20.5, + alpha_2=-0.8, + sigma_2=0.5, + ) + + assert result.shape == absolute_mag.shape + + +@pytest.mark.parametrize( + ("function", "kwargs"), + [ + ( + saunders_lf, + {"phi_star": 0.01, "m_star": -21.0, "alpha": -0.3, "sigma": 0.7}, + ), + ( + generalized_saunders_lf, + { + "phi_star": 0.01, + "m_star": -21.0, + "alpha": -0.3, + "sigma": 0.7, + "beta": 1.5, + }, + ), + ( + evolving_saunders_lf, + { + "redshift": 0.1, + "phi_star": 0.01, + "m_star": -21.0, + "alpha": -0.3, + "sigma": 0.7, + }, + ), + ( + double_saunders_lf, + { + "phi_star_1": 0.01, + "m_star_1": -21.0, + "alpha_1": -0.3, + "sigma_1": 0.7, + "phi_star_2": 0.005, + "m_star_2": -20.5, + "alpha_2": -0.8, + "sigma_2": 0.5, + }, + ), + ], +) +def test_saunders_models_accept_scalar_inputs(function, kwargs: dict[str, object]) -> None: + """Tests that all Saunders models accept scalar magnitude input.""" + result = function(-21.0, **kwargs) + + assert np.shape(result) == () + assert np.isfinite(result) + + +@pytest.mark.parametrize( + ("function", "kwargs"), + [ + ( + saunders_lf, + {"phi_star": 0.01, "m_star": -21.0, "alpha": -0.3, "sigma": 0.7}, + ), + ( + generalized_saunders_lf, + { + "phi_star": 0.01, + "m_star": -21.0, + "alpha": -0.3, + "sigma": 0.7, + "beta": 1.5, + }, + ), + ( + evolving_saunders_lf, + { + "redshift": 0.1, + "phi_star": 0.01, + "m_star": -21.0, + "alpha": -0.3, + "sigma": 0.7, + }, + ), + ( + double_saunders_lf, + { + "phi_star_1": 0.01, + "m_star_1": -21.0, + "alpha_1": -0.3, + "sigma_1": 0.7, + "phi_star_2": 0.005, + "m_star_2": -20.5, + "alpha_2": -0.8, + "sigma_2": 0.5, + }, + ), + ], +) +def test_saunders_models_reject_nonfinite_absolute_mag( + function, + kwargs: dict[str, object], +) -> None: + """Tests that all Saunders models reject non-finite magnitudes.""" + with pytest.raises(ValueError, match="absolute_mag contains NaN or infinite values"): + function(np.array([-22.0, np.nan]), **kwargs) + + +@pytest.mark.parametrize( + ("function", "kwargs", "bad_key", "match"), + [ + ( + saunders_lf, + {"phi_star": 0.01, "m_star": -21.0, "alpha": -0.3, "sigma": 0.7}, + "phi_star", + "phi_star contains NaN or infinite values", + ), + ( + saunders_lf, + {"phi_star": 0.01, "m_star": -21.0, "alpha": -0.3, "sigma": 0.7}, + "alpha", + "alpha contains NaN or infinite values", + ), + ( + saunders_lf, + {"phi_star": 0.01, "m_star": -21.0, "alpha": -0.3, "sigma": 0.7}, + "sigma", + "sigma contains NaN or infinite values", + ), + ( + generalized_saunders_lf, + { + "phi_star": 0.01, + "m_star": -21.0, + "alpha": -0.3, + "sigma": 0.7, + "beta": 1.5, + }, + "beta", + "beta contains NaN or infinite values", + ), + ( + evolving_saunders_lf, + { + "redshift": 0.1, + "phi_star": 0.01, + "m_star": -21.0, + "alpha": -0.3, + "sigma": 0.7, + }, + "redshift", + "redshift contains NaN or infinite values", + ), + ( + double_saunders_lf, + { + "phi_star_1": 0.01, + "m_star_1": -21.0, + "alpha_1": -0.3, + "sigma_1": 0.7, + "phi_star_2": 0.005, + "m_star_2": -20.5, + "alpha_2": -0.8, + "sigma_2": 0.5, + }, + "sigma_2", + "sigma contains NaN or infinite values", + ), + ], +) +def test_saunders_models_reject_nonfinite_validated_parameters( + function, + kwargs: dict[str, object], + bad_key: str, + match: str, +) -> None: + """Tests that validated Saunders parameters reject non-finite values.""" + kwargs = dict(kwargs) + kwargs[bad_key] = np.nan + + with pytest.raises(ValueError, match=match): + function(np.array([-22.0, -21.0]), **kwargs) + + +def test_saunders_models_return_float_arrays() -> None: + """Tests that all Saunders model outputs are floating point.""" + absolute_mag = np.array([-23, -22, -21]) + + results = [ + saunders_lf( + absolute_mag, + phi_star=1, + m_star=-21, + alpha=-0.3, + sigma=0.7, + ), + generalized_saunders_lf( + absolute_mag, + phi_star=1, + m_star=-21, + alpha=-0.3, + sigma=0.7, + beta=1.5, + ), + evolving_saunders_lf( + absolute_mag, + np.array([0.1, 0.2, 0.3]), + phi_star=1, + m_star=-21, + alpha=-0.3, + sigma=0.7, + ), + double_saunders_lf( + absolute_mag, + phi_star_1=1, + m_star_1=-21, + alpha_1=-0.3, + sigma_1=0.7, + phi_star_2=0.5, + m_star_2=-20, + alpha_2=-0.8, + sigma_2=0.5, + ), + ] + + for result in results: + assert result.dtype.kind == "f" diff --git a/tests/test_lumfuncs_models_schechter.py b/tests/test_lumfuncs_models_schechter.py index d133b1d3..f21deab1 100644 --- a/tests/test_lumfuncs_models_schechter.py +++ b/tests/test_lumfuncs_models_schechter.py @@ -8,6 +8,9 @@ schechter, evolving_schechter, double_schechter, + modified_schechter, + truncated_schechter, + multi_schechter, schechter_cumulative, schechter_cumulative_evolving, ) @@ -473,3 +476,340 @@ def test_double_schechter_accepts_array_phi_star() -> None: assert result.shape == (2,) assert np.all(result >= 0.0) + + +def test_modified_schechter_matches_schechter_when_beta_is_one() -> None: + """Tests that modified_schechter reduces to schechter for beta=1.""" + m = np.array([-22.0, -21.0, -20.0, -19.0]) + + result = modified_schechter( + m, + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + beta=1.0, + ) + expected = schechter( + m, + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + ) + + np.testing.assert_allclose(result, expected) + + +def test_modified_schechter_matches_manual_formula() -> None: + """Tests that modified_schechter matches the analytic formula.""" + m = np.array([-22.0, -21.0, -20.0]) + phi_star = 1e-3 + m_star = -20.0 + alpha = -1.2 + beta = 0.8 + + result = modified_schechter( + m, + phi_star=phi_star, + m_star=m_star, + alpha=alpha, + beta=beta, + ) + + x = luminosity_ratio(m, m_star) + expected = ( + 0.4 + * np.log(10.0) + * beta + * phi_star + * x ** (alpha + 1.0) + * np.exp(-(x**beta)) + ) + + np.testing.assert_allclose(result, expected) + + +def test_modified_schechter_rejects_negative_phi_star() -> None: + """Tests that modified_schechter rejects negative phi_star.""" + with pytest.raises(ValueError, match="phi_star must be non-negative"): + modified_schechter( + [-20.0], + phi_star=-1e-3, + m_star=-20.0, + alpha=-1.0, + beta=1.0, + ) + + +def test_modified_schechter_zero_phi_star_warning() -> None: + """Tests that modified_schechter warns when phi_star is zero.""" + with pytest.warns(UserWarning): + result = modified_schechter( + [-20.0], + phi_star=0.0, + m_star=-20.0, + alpha=-1.0, + beta=1.0, + ) + + np.testing.assert_allclose(result, np.array([0.0])) + + +def test_modified_schechter_rejects_nonpositive_beta() -> None: + """Tests that modified_schechter rejects beta <= 0.""" + with pytest.raises(ValueError, match="beta must be positive"): + modified_schechter( + [-20.0], + phi_star=1e-3, + m_star=-20.0, + alpha=-1.0, + beta=0.0, + ) + + +def test_modified_schechter_rejects_nonfinite_beta() -> None: + """Tests that modified_schechter rejects non-finite beta.""" + with pytest.raises(ValueError, match="beta contains NaN or infinite values"): + modified_schechter( + [-20.0], + phi_star=1e-3, + m_star=-20.0, + alpha=-1.0, + beta=np.inf, + ) + + +def test_modified_schechter_accepts_array_parameters() -> None: + """Tests that modified_schechter supports array-valued parameters.""" + m = np.array([-22.0, -21.0, -20.0]) + + result = modified_schechter( + m, + phi_star=np.array([1e-3, 2e-3, 3e-3]), + m_star=-20.0, + alpha=np.array([-1.2, -1.0, -0.8]), + beta=np.array([0.8, 1.0, 1.2]), + ) + + assert result.shape == m.shape + assert np.all(np.isfinite(result)) + assert np.all(result >= 0.0) + + +def test_truncated_schechter_matches_schechter_inside_limits() -> None: + """Tests that truncated_schechter matches schechter inside the limits.""" + m = np.array([-21.0, -20.0, -19.0]) + + result = truncated_schechter( + m, + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + m_bright=-22.0, + m_faint=-18.0, + ) + expected = schechter( + m, + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + ) + + np.testing.assert_allclose(result, expected) + + +def test_truncated_schechter_zeroes_values_outside_limits() -> None: + """Tests that truncated_schechter is zero outside the magnitude limits.""" + m = np.array([-23.0, -21.0, -20.0, -19.0, -17.0]) + + result = truncated_schechter( + m, + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + m_bright=-22.0, + m_faint=-18.0, + ) + + assert result[0] == 0.0 + assert result[-1] == 0.0 + assert np.all(result[1:-1] > 0.0) + + +def test_truncated_schechter_allows_only_bright_limit() -> None: + """Tests that truncated_schechter supports only a bright limit.""" + m = np.array([-23.0, -22.0, -21.0]) + + result = truncated_schechter( + m, + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + m_bright=-22.0, + ) + + assert result[0] == 0.0 + assert np.all(result[1:] > 0.0) + + +def test_truncated_schechter_allows_only_faint_limit() -> None: + """Tests that truncated_schechter supports only a faint limit.""" + m = np.array([-21.0, -20.0, -17.0]) + + result = truncated_schechter( + m, + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + m_faint=-18.0, + ) + + assert np.all(result[:2] > 0.0) + assert result[-1] == 0.0 + + +def test_truncated_schechter_rejects_nonfinite_bright_limit() -> None: + """Tests that truncated_schechter rejects non-finite m_bright.""" + with pytest.raises(ValueError, match="m_bright must be finite"): + truncated_schechter( + [-20.0], + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + m_bright=np.nan, + ) + + +def test_truncated_schechter_rejects_nonfinite_faint_limit() -> None: + """Tests that truncated_schechter rejects non-finite m_faint.""" + with pytest.raises(ValueError, match="m_faint must be finite"): + truncated_schechter( + [-20.0], + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + m_faint=np.inf, + ) + + +def test_truncated_schechter_rejects_reversed_limits() -> None: + """Tests that truncated_schechter rejects m_bright >= m_faint.""" + with pytest.raises(ValueError, match="m_bright must be less than m_faint"): + truncated_schechter( + [-20.0], + phi_star=1e-3, + m_star=-20.0, + alpha=-1.2, + m_bright=-18.0, + m_faint=-22.0, + ) + + +def test_truncated_schechter_propagates_schechter_validation() -> None: + """Tests that truncated_schechter propagates base Schechter validation.""" + with pytest.raises(ValueError, match="phi_star must be non-negative"): + truncated_schechter( + [-20.0], + phi_star=-1e-3, + m_star=-20.0, + alpha=-1.2, + m_bright=-22.0, + m_faint=-18.0, + ) + + +def test_multi_schechter_matches_sum_of_components() -> None: + """Tests that multi_schechter equals the sum of individual components.""" + m = np.array([-22.0, -21.0, -20.0]) + + result = multi_schechter( + m, + phi_stars=np.array([1e-3, 2e-3]), + m_stars=np.array([-20.0, -19.5]), + alphas=np.array([-1.2, -0.5]), + ) + + expected = ( + schechter(m, phi_star=1e-3, m_star=-20.0, alpha=-1.2) + + schechter(m, phi_star=2e-3, m_star=-19.5, alpha=-0.5) + ) + + np.testing.assert_allclose(result, expected) + + +def test_multi_schechter_single_component_matches_schechter() -> None: + """Tests that one-component multi_schechter matches schechter.""" + m = np.array([-22.0, -21.0, -20.0]) + + result = multi_schechter( + m, + phi_stars=np.array([1e-3]), + m_stars=np.array([-20.0]), + alphas=np.array([-1.2]), + ) + expected = schechter(m, phi_star=1e-3, m_star=-20.0, alpha=-1.2) + + np.testing.assert_allclose(result, expected) + + +def test_multi_schechter_rejects_mismatched_component_shapes() -> None: + """Tests that multi_schechter rejects mismatched component arrays.""" + with pytest.raises( + ValueError, + match="phi_stars, m_stars, and alphas must have matching shapes", + ): + multi_schechter( + [-20.0], + phi_stars=np.array([1e-3, 2e-3]), + m_stars=np.array([-20.0]), + alphas=np.array([-1.2, -0.5]), + ) + + +def test_multi_schechter_rejects_non_1d_component_arrays() -> None: + """Tests that multi_schechter rejects non-1D component arrays.""" + with pytest.raises( + ValueError, + match="phi_stars, m_stars, and alphas must be 1D arrays", + ): + multi_schechter( + [-20.0], + phi_stars=np.array([[1e-3, 2e-3]]), + m_stars=np.array([[-20.0, -19.5]]), + alphas=np.array([[-1.2, -0.5]]), + ) + + +def test_multi_schechter_rejects_negative_phi_stars() -> None: + """Tests that multi_schechter rejects negative component normalizations.""" + with pytest.raises(ValueError, match="phi_stars must be non-negative"): + multi_schechter( + [-20.0], + phi_stars=np.array([1e-3, -2e-3]), + m_stars=np.array([-20.0, -19.5]), + alphas=np.array([-1.2, -0.5]), + ) + + +def test_multi_schechter_zero_phi_star_component_warning() -> None: + """Tests that multi_schechter warns when a component has zero normalization.""" + with pytest.warns(UserWarning): + result = multi_schechter( + [-20.0], + phi_stars=np.array([1e-3, 0.0]), + m_stars=np.array([-20.0, -19.5]), + alphas=np.array([-1.2, -0.5]), + ) + + expected = schechter([-20.0], phi_star=1e-3, m_star=-20.0, alpha=-1.2) + np.testing.assert_allclose(result, expected) + + +def test_multi_schechter_rejects_nonfinite_component_values() -> None: + """Tests that multi_schechter rejects non-finite component arrays.""" + with pytest.raises(ValueError, match="m_stars contains NaN or infinite values"): + multi_schechter( + [-20.0], + phi_stars=np.array([1e-3]), + m_stars=np.array([np.nan]), + alphas=np.array([-1.2]), + ) diff --git a/tests/test_lumfuncs_registry.py b/tests/test_lumfuncs_registry.py index 20abb234..367b0078 100644 --- a/tests/test_lumfuncs_registry.py +++ b/tests/test_lumfuncs_registry.py @@ -435,3 +435,29 @@ def test_get_lf_from_m_model_raises_for_unknown_name( match="phi_from_m is not defined for luminosity function model 'bad'", ): registry.get_lf_from_m_model("bad") + + +@pytest.mark.parametrize( + "name", + [ + "saunders", + "evolving_saunders", + "double_saunders", + "generalized_saunders", + ], +) +def test_saunders_models_are_registered(name: str) -> None: + """Tests that Saunders models are discovered by the LF registry.""" + assert name in registry.LF_MODELS + + +@pytest.mark.parametrize( + "name", + [ + "gamma", + "generalized_gamma", + ], +) +def test_gamma_models_are_registered(name: str) -> None: + """Tests that gamma models are discovered by the LF registry.""" + assert name in registry.LF_MODELS diff --git a/tests/test_utils_validators.py b/tests/test_utils_validators.py index 674bceec..ca6abfd9 100644 --- a/tests/test_utils_validators.py +++ b/tests/test_utils_validators.py @@ -7,6 +7,11 @@ validate_array, validate_luminosity_distance, validate_magnitude_range, + validate_2d_binned_grid, + validate_2d_tabulated_grid, + validate_binned_grid, + validate_strictly_increasing_1d, + validate_tabulated_grid, ) @@ -269,3 +274,390 @@ def test_validate_magnitude_range_rejects_negative_infinite_faint_bound() -> Non """Tests that negative infinite faint 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_strictly_increasing_1d_accepts_valid_grid() -> None: + """Tests that a valid increasing one-dimensional grid is accepted.""" + result = validate_strictly_increasing_1d([0.0, 1.0, 2.0], name="z") + + np.testing.assert_allclose(result, np.array([0.0, 1.0, 2.0])) + assert result.dtype == np.float64 + + +def test_validate_strictly_increasing_1d_rejects_multidimensional_input() -> None: + """Tests that multidimensional inputs are rejected.""" + with pytest.raises(ValueError, match="z must be one-dimensional"): + validate_strictly_increasing_1d([[0.0, 1.0], [2.0, 3.0]], name="z") + + +def test_validate_strictly_increasing_1d_rejects_too_few_values() -> None: + """Tests that grids with too few values are rejected.""" + with pytest.raises(ValueError, match="z must contain at least 2 values"): + validate_strictly_increasing_1d([0.0], name="z") + + +def test_validate_strictly_increasing_1d_respects_custom_min_size() -> None: + """Tests that the custom minimum size is enforced.""" + with pytest.raises(ValueError, match="z must contain at least 3 values"): + validate_strictly_increasing_1d([0.0, 1.0], name="z", min_size=3) + + +def test_validate_strictly_increasing_1d_rejects_equal_values() -> None: + """Tests that repeated values are rejected.""" + with pytest.raises(ValueError, match="z must be strictly increasing"): + validate_strictly_increasing_1d([0.0, 1.0, 1.0], name="z") + + +def test_validate_strictly_increasing_1d_rejects_decreasing_values() -> None: + """Tests that decreasing grids are rejected.""" + with pytest.raises(ValueError, match="z must be strictly increasing"): + validate_strictly_increasing_1d([0.0, 2.0, 1.0], name="z") + + +def test_validate_strictly_increasing_1d_rejects_negative_values_when_disallowed() -> None: + """Tests that negative grid values are rejected when requested.""" + with pytest.raises(ValueError, match="z contains negative values"): + validate_strictly_increasing_1d([-1.0, 0.0, 1.0], name="z", allow_negative=False) + + +def test_validate_tabulated_grid_accepts_valid_nonnegative_values() -> None: + """Tests that a valid tabulated grid is accepted.""" + x, y = validate_tabulated_grid( + [0.0, 1.0, 2.0], + [0.0, 2.0, 4.0], + coordinate_name="x", + values_name="y", + ) + + np.testing.assert_allclose(x, np.array([0.0, 1.0, 2.0])) + np.testing.assert_allclose(y, np.array([0.0, 2.0, 4.0])) + + +def test_validate_tabulated_grid_rejects_multidimensional_values() -> None: + """Tests that tabulated values must be one-dimensional.""" + with pytest.raises(ValueError, match="y must be one-dimensional"): + validate_tabulated_grid( + [0.0, 1.0], + [[1.0, 2.0]], + coordinate_name="x", + values_name="y", + ) + + +def test_validate_tabulated_grid_rejects_length_mismatch() -> None: + """Tests that coordinate and value lengths must match.""" + with pytest.raises(ValueError, match="x and y must have the same length"): + validate_tabulated_grid( + [0.0, 1.0, 2.0], + [1.0, 2.0], + coordinate_name="x", + values_name="y", + ) + + +def test_validate_tabulated_grid_rejects_negative_values_by_default() -> None: + """Tests that negative tabulated values are rejected by default.""" + with pytest.raises(ValueError, match="y must be non-negative"): + validate_tabulated_grid( + [0.0, 1.0, 2.0], + [1.0, -1.0, 2.0], + coordinate_name="x", + values_name="y", + ) + + +def test_validate_tabulated_grid_rejects_zero_when_positive_values_required() -> None: + """Tests that zero values are rejected when positive values are required.""" + with pytest.raises(ValueError, match="y must be positive"): + validate_tabulated_grid( + [0.0, 1.0, 2.0], + [1.0, 0.0, 2.0], + coordinate_name="x", + values_name="y", + positive_values=True, + ) + + +def test_validate_tabulated_grid_rejects_negative_coordinate_when_disallowed() -> None: + """Tests that negative coordinates are rejected when requested.""" + with pytest.raises(ValueError, match="x contains negative values"): + validate_tabulated_grid( + [-1.0, 0.0, 1.0], + [1.0, 2.0, 3.0], + coordinate_name="x", + values_name="y", + allow_negative_coordinate=False, + ) + + +def test_validate_binned_grid_accepts_valid_nonnegative_values() -> None: + """Tests that a valid binned grid is accepted.""" + edges, values = validate_binned_grid( + [0.0, 1.0, 2.0], + [0.0, 2.0], + edges_name="edges", + values_name="counts", + ) + + np.testing.assert_allclose(edges, np.array([0.0, 1.0, 2.0])) + np.testing.assert_allclose(values, np.array([0.0, 2.0])) + + +def test_validate_binned_grid_rejects_multidimensional_values() -> None: + """Tests that binned values must be one-dimensional.""" + with pytest.raises(ValueError, match="counts must be one-dimensional"): + validate_binned_grid( + [0.0, 1.0, 2.0], + [[1.0, 2.0]], + edges_name="edges", + values_name="counts", + ) + + +def test_validate_binned_grid_rejects_wrong_edge_count() -> None: + """Tests that bin edges must have one more value than bin values.""" + with pytest.raises( + ValueError, + match="edges must have one more value than counts", + ): + validate_binned_grid( + [0.0, 1.0, 2.0], + [1.0, 2.0, 3.0], + edges_name="edges", + values_name="counts", + ) + + +def test_validate_binned_grid_rejects_negative_values_by_default() -> None: + """Tests that negative binned values are rejected by default.""" + with pytest.raises(ValueError, match="counts must be non-negative"): + validate_binned_grid( + [0.0, 1.0, 2.0], + [1.0, -1.0], + edges_name="edges", + values_name="counts", + ) + + +def test_validate_binned_grid_rejects_zero_when_positive_values_required() -> None: + """Tests that zero binned values are rejected when positive values are required.""" + with pytest.raises(ValueError, match="counts must be positive"): + validate_binned_grid( + [0.0, 1.0, 2.0], + [1.0, 0.0], + edges_name="edges", + values_name="counts", + positive_values=True, + ) + + +def test_validate_binned_grid_rejects_negative_edges_when_disallowed() -> None: + """Tests that negative bin edges are rejected when requested.""" + with pytest.raises(ValueError, match="edges contains negative values"): + validate_binned_grid( + [-1.0, 0.0, 1.0], + [1.0, 2.0], + edges_name="edges", + values_name="counts", + allow_negative_edges=False, + ) + + +def test_validate_2d_tabulated_grid_accepts_valid_values() -> None: + """Tests that a valid two-dimensional tabulated grid is accepted.""" + x, y, values = validate_2d_tabulated_grid( + [0.0, 1.0, 2.0], + [10.0, 20.0], + [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], + x_name="x", + y_name="y", + values_name="phi", + ) + + np.testing.assert_allclose(x, np.array([0.0, 1.0, 2.0])) + np.testing.assert_allclose(y, np.array([10.0, 20.0])) + np.testing.assert_allclose(values, np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])) + + +def test_validate_2d_tabulated_grid_rejects_non_2d_values() -> None: + """Tests that tabulated 2D values must be two-dimensional.""" + with pytest.raises(ValueError, match="phi must be two-dimensional"): + validate_2d_tabulated_grid( + [0.0, 1.0], + [10.0, 20.0], + [1.0, 2.0], + x_name="x", + y_name="y", + values_name="phi", + ) + + +def test_validate_2d_tabulated_grid_rejects_wrong_shape() -> None: + """Tests that tabulated 2D values must match y-by-x shape.""" + with pytest.raises( + ValueError, + match=r"phi must have shape \(y.size, x.size\)", + ): + validate_2d_tabulated_grid( + [0.0, 1.0, 2.0], + [10.0, 20.0], + [[1.0, 2.0], [3.0, 4.0]], + x_name="x", + y_name="y", + values_name="phi", + ) + + +def test_validate_2d_tabulated_grid_rejects_negative_values_by_default() -> None: + """Tests that negative tabulated 2D values are rejected by default.""" + with pytest.raises(ValueError, match="phi must be non-negative"): + validate_2d_tabulated_grid( + [0.0, 1.0], + [10.0, 20.0], + [[1.0, -1.0], [2.0, 3.0]], + x_name="x", + y_name="y", + values_name="phi", + ) + + +def test_validate_2d_tabulated_grid_rejects_zero_when_positive_values_required() -> None: + """Tests that zero tabulated 2D values are rejected when positivity is required.""" + with pytest.raises(ValueError, match="phi must be positive"): + validate_2d_tabulated_grid( + [0.0, 1.0], + [10.0, 20.0], + [[1.0, 0.0], [2.0, 3.0]], + x_name="x", + y_name="y", + values_name="phi", + positive_values=True, + ) + + +def test_validate_2d_tabulated_grid_rejects_negative_x_when_disallowed() -> None: + """Tests that negative x coordinates are rejected when requested.""" + with pytest.raises(ValueError, match="x contains negative values"): + validate_2d_tabulated_grid( + [-1.0, 0.0], + [10.0, 20.0], + [[1.0, 2.0], [3.0, 4.0]], + x_name="x", + y_name="y", + values_name="phi", + allow_negative_x=False, + ) + + +def test_validate_2d_tabulated_grid_rejects_negative_y_when_disallowed() -> None: + """Tests that negative y coordinates are rejected when requested.""" + with pytest.raises(ValueError, match="y contains negative values"): + validate_2d_tabulated_grid( + [0.0, 1.0], + [-10.0, 20.0], + [[1.0, 2.0], [3.0, 4.0]], + x_name="x", + y_name="y", + values_name="phi", + allow_negative_y=False, + ) + + +def test_validate_2d_binned_grid_accepts_valid_values() -> None: + """Tests that a valid two-dimensional binned grid is accepted.""" + x_edges, y_edges, values = validate_2d_binned_grid( + [0.0, 1.0, 2.0], + [10.0, 20.0, 30.0], + [[1.0, 2.0], [3.0, 4.0]], + x_edges_name="x_edges", + y_edges_name="y_edges", + values_name="counts", + ) + + np.testing.assert_allclose(x_edges, np.array([0.0, 1.0, 2.0])) + np.testing.assert_allclose(y_edges, np.array([10.0, 20.0, 30.0])) + np.testing.assert_allclose(values, np.array([[1.0, 2.0], [3.0, 4.0]])) + + +def test_validate_2d_binned_grid_rejects_non_2d_values() -> None: + """Tests that binned 2D values must be two-dimensional.""" + with pytest.raises(ValueError, match="counts must be two-dimensional"): + validate_2d_binned_grid( + [0.0, 1.0, 2.0], + [10.0, 20.0, 30.0], + [1.0, 2.0], + x_edges_name="x_edges", + y_edges_name="y_edges", + values_name="counts", + ) + + +def test_validate_2d_binned_grid_rejects_wrong_shape() -> None: + """Tests that binned 2D values must match y-bin by x-bin shape.""" + with pytest.raises( + ValueError, + match=r"counts must have shape \(y_edges.size - 1, x_edges.size - 1\)", + ): + validate_2d_binned_grid( + [0.0, 1.0, 2.0], + [10.0, 20.0, 30.0], + [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]], + x_edges_name="x_edges", + y_edges_name="y_edges", + values_name="counts", + ) + + +def test_validate_2d_binned_grid_rejects_negative_values_by_default() -> None: + """Tests that negative binned 2D values are rejected by default.""" + with pytest.raises(ValueError, match="counts must be non-negative"): + validate_2d_binned_grid( + [0.0, 1.0, 2.0], + [10.0, 20.0, 30.0], + [[1.0, -1.0], [2.0, 3.0]], + x_edges_name="x_edges", + y_edges_name="y_edges", + values_name="counts", + ) + + +def test_validate_2d_binned_grid_rejects_zero_when_positive_values_required() -> None: + """Tests that zero binned 2D values are rejected when positivity is required.""" + with pytest.raises(ValueError, match="counts must be positive"): + validate_2d_binned_grid( + [0.0, 1.0, 2.0], + [10.0, 20.0, 30.0], + [[1.0, 0.0], [2.0, 3.0]], + x_edges_name="x_edges", + y_edges_name="y_edges", + values_name="counts", + positive_values=True, + ) + + +def test_validate_2d_binned_grid_rejects_negative_x_edges_when_disallowed() -> None: + """Tests that negative x bin edges are rejected when requested.""" + with pytest.raises(ValueError, match="x_edges contains negative values"): + validate_2d_binned_grid( + [-1.0, 0.0, 1.0], + [10.0, 20.0, 30.0], + [[1.0, 2.0], [3.0, 4.0]], + x_edges_name="x_edges", + y_edges_name="y_edges", + values_name="counts", + allow_negative_x_edges=False, + ) + + +def test_validate_2d_binned_grid_rejects_negative_y_edges_when_disallowed() -> None: + """Tests that negative y bin edges are rejected when requested.""" + with pytest.raises(ValueError, match="y_edges contains negative values"): + validate_2d_binned_grid( + [0.0, 1.0, 2.0], + [-10.0, 20.0, 30.0], + [[1.0, 2.0], [3.0, 4.0]], + x_edges_name="x_edges", + y_edges_name="y_edges", + values_name="counts", + allow_negative_y_edges=False, + )