diff --git a/docs/core.rst b/docs/core.rst index 3f8178cd..2767ef5d 100644 --- a/docs/core.rst +++ b/docs/core.rst @@ -128,7 +128,11 @@ Domains (:mod:`hmf.core.domain`) Every model declares a *valid* domain (outside which it always raises a :class:`~hmf.core.domain.DomainError`) and a *calibration* domain (outside which a user-chosen :data:`~hmf.core.domain.DomainPolicy` applies: ``"ignore"``, - ``"warn"``, ``"mask"`` or ``"raise"``), or ``None`` for no calibration domain. A + ``"warn"``, ``"mask"`` or ``"raise"``), or ``None`` for no calibration domain. + Ask a model instance for its calibration domain, ``model.calibration`` + (:class:`~hmf.core.domain.Calibrated`): it can depend on the model's parameters + (e.g. :class:`~hmf.core.concentration.Duffy08`'s redshift range), and defaults to + the class's ``calibration_domain``, that of its default parameters. A :class:`~hmf.core.domain.Domain` bounds each variable by an interval whose ends are closed (``>=``, ``<=``) or open (``>``, ``<``): ``(0, None, "(]")`` is ``> 0``. :meth:`~hmf.core.domain.Domain.describe` writes it as inequalities, and diff --git a/docs/templates/core-model-class.rst b/docs/templates/core-model-class.rst index 0544064c..66dc6f80 100644 --- a/docs/templates/core-model-class.rst +++ b/docs/templates/core-model-class.rst @@ -8,7 +8,7 @@ {#- A model's parameters are listed in the "Parameters" section generated from its fields. They get no pages of their own: names that differ only in case (A_200, a_200) would clash. Only the model's metadata is listed here. #} - {% set metadata = ["requires", "valid_domain", "calibration_domain", + {% set metadata = ["requires", "valid_domain", "calibration_domain", "calibration", "measured_mass_definition", "measured_mass_definition_note", "simulations", "references", "parameter_source", diff --git a/src/hmf/core/_cosmology_models.py b/src/hmf/core/_cosmology_models.py index 6f569b72..ee025f70 100644 --- a/src/hmf/core/_cosmology_models.py +++ b/src/hmf/core/_cosmology_models.py @@ -3,7 +3,8 @@ :class:`_CosmologyModel` holds the class-level scaffolding that such a kind (:class:`~hmf.core.transfer_models.TransferModel`, :class:`~hmf.core.growth_models.GrowthModel`) declares: the Boltzmann ``backend``, -the ``valid_domain`` and ``calibration_domain``, and the ``check_cosmology`` hook. +the ``valid_domain``, the calibration domain (from +:class:`~hmf.core.domain.Calibrated`), and the ``check_cosmology`` hook. :class:`_BoltzmannBacked` is the cosmology check of the models computed by a Boltzmann code (CAMB or CLASS), for transfer and growth alike. @@ -18,10 +19,10 @@ from astropy.cosmology import FLRW from ._boltzmann import check_boltzmann_cosmology -from .domain import Domain +from .domain import Calibrated, Domain -class _CosmologyModel: +class _CosmologyModel(Calibrated): """The class-level scaffolding of a kind of models computed from a cosmology. The kind sets :attr:`valid_domain`; the other class variables default to @@ -37,9 +38,6 @@ class _CosmologyModel: #: model's class's domain, so a model can narrow it. valid_domain: ClassVar[Domain] - #: Where the model was calibrated, if that is stated by its source. - calibration_domain: ClassVar[Domain | None] = None - def check_cosmology(self, cosmology: FLRW) -> None: """Raise if the model does not apply to ``cosmology`` (by default it does). diff --git a/src/hmf/core/_kernels/mass_conversion.py b/src/hmf/core/_kernels/mass_conversion.py index c7cdd993..55511bc9 100644 --- a/src/hmf/core/_kernels/mass_conversion.py +++ b/src/hmf/core/_kernels/mass_conversion.py @@ -26,8 +26,9 @@ point, each elementwise, so the result does not depend on the batch size, and :func:`profile_convert` gives the masses and the analytic Jacobian. For NFW (:func:`nfw_ln_mu_and_slope`) the starting point is the fit of Hu & Kravtsov (2003, -eq. C11; :func:`hu_kravtsov_inverse`), and :func:`nfw_concentration_ratio_solve` and -:func:`nfw_convert` are the two with it. +eq. C11; :func:`nfw_start`). :class:`~hmf.core.profiles.NFW` calls the two with +them, through :meth:`DensityProfile.convert_concentration_kernel +`. The functions follow the rules of :mod:`hmf.core._kernels` and do not check their inputs: the kernel-level entry points of :mod:`hmf.core.profiles` do. @@ -45,12 +46,8 @@ "LnMuAndSlope", "concentration_ratio_solve", "hu_kravtsov_inverse", - "nfw_concentration_ratio_solve", - "nfw_convert", "nfw_ln_mean_density", - "nfw_ln_mu", "nfw_ln_mu_and_slope", - "nfw_ln_mu_slope", "nfw_start", "power_law_start", "profile_convert", @@ -61,14 +58,15 @@ #: A profile's :math:`x \mapsto (\ln\mu(x), \psi(x))`, on plain arrays. LnMuAndSlope = Callable[[FloatArray], tuple[FloatArray, FloatArray]] -#: The number of Newton steps that polish the Hu & Kravtsov (2003) starting point. +#: The number of Newton steps that polish the Hu & Kravtsov (2003) starting point +#: (:attr:`hmf.core.profiles.NFW.newton_steps`). #: Newton converges quadratically from its error (below 1% for 0 < c < 400): for #: 1e-3 < c < 1e4 and 1e-4 < q < 1e4 two steps leave 1e-12 and three round-off #: (~1e-13 in c); the fourth is a margin. NEWTON_STEPS = 4 -#: Below this x, :func:`nfw_ln_mu` uses the series of mu(x), whose closed form loses -#: digits to cancellation (relative error ~ 1e-16 / x^2). +#: Below this x, :func:`nfw_ln_mu_and_slope` uses the series of mu(x), whose closed +#: form loses digits to cancellation (relative error ~ 1e-16 / x^2). _SERIES_X = 0.01 #: The parameters (a1, a2, a3, a4) of Hu & Kravtsov (2003), eq. C11. @@ -78,7 +76,11 @@ def nfw_ln_mu_and_slope(x: npt.ArrayLike) -> tuple[FloatArray, FloatArray]: r"""The NFW :math:`\ln\mu(x)` and :math:`\psi(x)`, from one evaluation of :math:`\mu`. - See :func:`nfw_ln_mu` and :func:`nfw_ln_mu_slope`. + :math:`\mu(x) = \ln(1+x) - x/(1+x)`; below x = 0.01 it uses the series + :math:`\mu(x) = \sum_{n\ge2} (-1)^n (n-1)x^n/n` (to :math:`x^{9}`), which avoids + the cancellation in the closed form. The slope + :math:`\psi(x) = d\ln\mu/d\ln x = x^2/((1+x)^2\mu(x))` falls from 2 (as + :math:`x \to 0`) to 0 (as :math:`x \to \infty`). Parameters ---------- @@ -106,44 +108,6 @@ def nfw_ln_mu_and_slope(x: npt.ArrayLike) -> tuple[FloatArray, FloatArray]: return ln_mu, psi -def nfw_ln_mu(x: npt.ArrayLike) -> FloatArray: - r"""The log of the NFW mass function :math:`\mu(x) = \ln(1+x) - x/(1+x)`. - - Below x = 0.01 it uses the series - :math:`\mu(x) = \sum_{n\ge2} (-1)^n (n-1)x^n/n` (to :math:`x^{9}`), which avoids - the cancellation in the closed form. - - Parameters - ---------- - x - :math:`r/r_s` (> 0), dimensionless. - - Returns - ------- - numpy.ndarray - :math:`\ln\mu(x)`. - """ - return nfw_ln_mu_and_slope(np.asarray(x, dtype=np.float64))[0] - - -def nfw_ln_mu_slope(x: npt.ArrayLike) -> FloatArray: - r"""The logarithmic slope :math:`\psi(x) = d\ln\mu/d\ln x = x^2/((1+x)^2\mu(x))`. - - It falls from 2 (as :math:`x \to 0`) to 0 (as :math:`x \to \infty`). - - Parameters - ---------- - x - :math:`r/r_s` (> 0), dimensionless. - - Returns - ------- - numpy.ndarray - :math:`\psi(x)`, dimensionless. - """ - return nfw_ln_mu_and_slope(np.asarray(x, dtype=np.float64))[1] - - def nfw_ln_mean_density(x: npt.ArrayLike) -> FloatArray: r""":math:`\ln F(x) = \ln[\mu(x)/x^3]`: the mean enclosed density over :math:`3\rho_s`. @@ -161,7 +125,7 @@ def nfw_ln_mean_density(x: npt.ArrayLike) -> FloatArray: :math:`\ln F(x)`. """ x = np.asarray(x, dtype=np.float64) - out: FloatArray = nfw_ln_mu(x) - 3.0 * np.log(x) + out: FloatArray = nfw_ln_mu_and_slope(x)[0] - 3.0 * np.log(x) return out @@ -328,50 +292,3 @@ def nfw_start(c: FloatArray, ln_q: FloatArray) -> FloatArray: r"""The NFW starting point: Hu & Kravtsov (2003) at :math:`\ln F(c) + \ln q`.""" c = np.asarray(c, dtype=np.float64) return hu_kravtsov_inverse(nfw_ln_mean_density(c) + ln_q) - - -def nfw_concentration_ratio_solve(c: npt.ArrayLike, ln_q: npt.ArrayLike) -> FloatArray: - r"""The NFW concentration :math:`c'` with :math:`F(c') = q F(c)`, elementwise. - - :func:`concentration_ratio_solve` for NFW, from :func:`hu_kravtsov_inverse`, with - :data:`NEWTON_STEPS` steps, so the error is at round-off. - - Parameters - ---------- - c - The concentrations (> 0), dimensionless. - ln_q - :math:`\ln q`, broadcast with ``c``. - - Returns - ------- - numpy.ndarray - :math:`c'`, with the broadcast shape. - """ - c = np.asarray(c, dtype=np.float64) - ln_q = np.asarray(ln_q, dtype=np.float64) - return concentration_ratio_solve(nfw_ln_mu_and_slope, c, ln_q, nfw_start(c, ln_q), NEWTON_STEPS) - - -def nfw_convert( - c: npt.ArrayLike, slope: npt.ArrayLike, ln_q: npt.ArrayLike -) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray]: - r""":func:`profile_convert` for NFW haloes (see :func:`nfw_concentration_ratio_solve`). - - Parameters - ---------- - c - The concentrations in the old definition (> 0). - slope - :math:`s = d\ln c/d\ln M`, in the old definition. - ln_q - :math:`\ln q = \ln(\rho_{\rm new}/\rho_{\rm old})`. - - Returns - ------- - ln_mass_ratio, c_new, dln_m_new_dln_m, dln_c_new_dln_m : numpy.ndarray - As for :func:`profile_convert`. - """ - c = np.asarray(c, dtype=np.float64) - ln_q = np.asarray(ln_q, dtype=np.float64) - return profile_convert(nfw_ln_mu_and_slope, c, slope, ln_q, nfw_start(c, ln_q), NEWTON_STEPS) diff --git a/src/hmf/core/concentration.py b/src/hmf/core/concentration.py index 3abf3cf1..8afc5320 100644 --- a/src/hmf/core/concentration.py +++ b/src/hmf/core/concentration.py @@ -10,8 +10,9 @@ Every relation is a :class:`ConcentrationModel`, a model kind with its own registry (``ConcentrationModel.get("Duffy08")``), so other packages (e.g. halomod) add theirs by subclassing it. Each declares a :attr:`~ConcentrationModel.valid_domain` (outside -which it raises) and a :attr:`~ConcentrationModel.calibration_domain`, in its own -mass definition. +which it raises) and a calibration domain, in its own mass definition: an instance's +is :attr:`~ConcentrationModel.calibration` (it can depend on the instance's +parameters; see :class:`~hmf.core.domain.Calibrated`). Implemented: @@ -34,7 +35,7 @@ from . import _references as refs from ._fields import field from ._validators import one_of -from .domain import Domain, check_extent +from .domain import Calibrated, Domain, check_extent from .mass_definitions import MassDefinition, SOCritical, SOMean, SOVirial from .model import Model from .units import Msun_h, unit_boundary @@ -46,7 +47,7 @@ @attrs.frozen(kw_only=True) -class ConcentrationModel(Model, kind=True): +class ConcentrationModel(Calibrated, Model, kind=True): r"""The kind of concentration-mass relations: :math:`c(M, z)`. Subclasses implement :meth:`_ln_concentration_and_slope` (on checked plain arrays) @@ -60,19 +61,11 @@ class ConcentrationModel(Model, kind=True): {"m": (0 * Msun_h, None, "(]"), "z": (0, None)}, source="m > 0, z >= 0." ) - #: Where the relation was calibrated, with its source (``None`` if it was not). - calibration_domain: ClassVar[Domain | None] = None - @property @abc.abstractmethod def mass_definition(self) -> MassDefinition: """The mass definition the relation's masses and concentrations are in.""" - @property - def calibration(self) -> Domain | None: - """The calibration domain of this instance (the class's, unless it depends on it).""" - return type(self).calibration_domain - @abc.abstractmethod def _ln_concentration_and_slope( self, m: FloatArray, z: FloatArray @@ -80,7 +73,7 @@ def _ln_concentration_and_slope( """Ln c and dln c / dln m, on checked plain arrays (m in Msun/h).""" def ln_concentration_and_slope_kernel( - self, m: npt.ArrayLike, *, z: npt.ArrayLike + self, *, m: npt.ArrayLike, z: npt.ArrayLike ) -> tuple[FloatArray, FloatArray]: r"""Ln c(M, z) and its slope :math:`d\ln c/d\ln M`, at kernel level. @@ -112,7 +105,7 @@ def ln_concentration_and_slope_kernel( np.array(np.broadcast_to(slope, shape)), ) - def in_calibration_domain_kernel(self, m: npt.ArrayLike, *, z: npt.ArrayLike) -> BoolArray: + def in_calibration_domain_kernel(self, *, m: npt.ArrayLike, z: npt.ArrayLike) -> BoolArray: """Whether (m, z) is inside :attr:`calibration`, at kernel level. Parameters @@ -145,7 +138,7 @@ def in_calibration_domain_kernel(self, m: npt.ArrayLike, *, z: npt.ArrayLike) -> ) @unit_boundary(m=Msun_h) - def concentration(self, m: Any, *, z: npt.ArrayLike) -> FloatArray: + def concentration(self, *, m: Any, z: npt.ArrayLike) -> FloatArray: """The median concentration c(M, z), in :attr:`mass_definition`. Parameters @@ -165,7 +158,7 @@ def concentration(self, m: Any, *, z: npt.ArrayLike) -> FloatArray: DomainError Outside :attr:`valid_domain`. """ - return np.exp(self.ln_concentration_and_slope_kernel(m, z=z)[0]) + return np.exp(self.ln_concentration_and_slope_kernel(m=m, z=z)[0]) # --------------------------------------------------------------------------------- @@ -259,8 +252,8 @@ class Duffy08(ConcentrationModel, alias="Duffy08"): doc="The redshift range the set was fitted to: '0-2', or '0' (with C = 0).", ) - #: The calibration domain of the default set (z = 0-2); :attr:`calibration` is that - #: of the instance's set. + #: The calibration domain of the default set (z = 0-2): ask an instance for + #: :attr:`calibration`, that of its set. calibration_domain: ClassVar[Domain | None] = _DUFFY08_DOMAINS["0-2"] #: The pivot mass, in Msun/h. diff --git a/src/hmf/core/domain.py b/src/hmf/core/domain.py index b5801507..e2cf5676 100644 --- a/src/hmf/core/domain.py +++ b/src/hmf/core/domain.py @@ -1,14 +1,21 @@ """Model domains, the domain policy, and the errors of :mod:`hmf.core`. -Each model carries two domains, as class variables: +Each model carries two domains: ``valid_domain`` - Where the model's formula is mathematically defined or sensible (e.g. sigma > 0). - Outside it, evaluation **always raises** a :class:`DomainError`. -``calibration_domain`` - Where the model was calibrated against simulations, or ``None`` if it has none. - Outside it, results are extrapolations, handled according to the user's - :data:`DomainPolicy`. + Where the model's formula is mathematically defined or sensible (e.g. sigma > 0), + a class variable. Outside it, evaluation **always raises** a :class:`DomainError`. +``calibration`` + Where the model was calibrated against simulations, or ``None`` if it has none: + an instance property (:attr:`Calibrated.calibration`), since it can depend on the + model's parameters (e.g. :class:`~hmf.core.concentration.Duffy08`'s redshift + range). Outside it, results are extrapolations, handled according to the user's + :data:`DomainPolicy`. By default it is the class variable ``calibration_domain``, + the calibration of the model's default parameters. + +Every model kind with a calibration domain (fits, concentration-mass relations, +density profiles, filters, transfer and growth models) is a :class:`Calibrated`: +ask an instance for ``model.calibration``. This module provides the :class:`Domain` type describing either (with :meth:`Domain.describe` and :meth:`Domain.check`), :func:`apply_domain_policy`, @@ -52,7 +59,7 @@ from collections.abc import Mapping from functools import cached_property from types import MappingProxyType -from typing import Any, Literal, get_args +from typing import Any, ClassVar, Literal, get_args import astropy.units as u import attrs @@ -64,6 +71,7 @@ __all__ = [ "DOMAIN_POLICIES", + "Calibrated", "Domain", "DomainError", "DomainPolicy", @@ -402,6 +410,30 @@ def check(self, values: Mapping[str, Any], *, where: str) -> None: ) +class Calibrated: + """The calibration domain of a model kind (a mixin, listed before the kind's base). + + A kind sets :attr:`calibration_domain` (``None`` by default); a model whose + calibration depends on its parameters overrides :attr:`calibration`. + """ + + __slots__ = () + + #: The calibration domain of the class, with its source (``None`` if it was not + #: calibrated): that of its default parameters. Ask an instance for + #: :attr:`calibration`. + calibration_domain: ClassVar[Domain | None] = None + + @property + def calibration(self) -> Domain | None: + """Where this model was calibrated, with its source (``None`` if it was not). + + The class's :attr:`calibration_domain`, unless it depends on the instance's + parameters. + """ + return type(self).calibration_domain + + def check_extent( name: str, x: npt.ArrayLike, diff --git a/src/hmf/core/filters.py b/src/hmf/core/filters.py index 5225f86e..9f5c8ed5 100644 --- a/src/hmf/core/filters.py +++ b/src/hmf/core/filters.py @@ -27,7 +27,7 @@ from ._fields import field from ._kernels import filters as kernels from ._validators import positive -from .domain import Domain, check_extent +from .domain import Calibrated, Domain, check_extent from .model import Model from .units import unit_boundary @@ -65,7 +65,7 @@ class TailBounds(NamedTuple): @attrs.frozen(kw_only=True) -class Filter(Model, kind=True): +class Filter(Calibrated, Model, kind=True): r"""The kind of smoothing filters. A filter defines: diff --git a/src/hmf/core/fits.py b/src/hmf/core/fits.py index 1dee1562..6f8eeb3c 100644 --- a/src/hmf/core/fits.py +++ b/src/hmf/core/fits.py @@ -84,10 +84,12 @@ :attr:`FittingFunction.valid_domain` where its formula is defined, or sensible. :meth:`FittingFunction.fsigma` **always raises** :class:`~hmf.core.domain.DomainError` outside it. -:attr:`FittingFunction.calibration_domain` +:attr:`FittingFunction.calibration` where it was calibrated against simulations, with its source (``None`` for a fit that was not, such as :class:`PS`). Outside it, the result is an extrapolation, handled by a :data:`~hmf.core.domain.DomainPolicy` with :func:`evaluate_fsigma`. + It is an instance property (see :class:`~hmf.core.domain.Calibrated`); each fit + declares it as the class variable :attr:`FittingFunction.calibration_domain`. The bounds of both are in these variables: the inputs above, and the derived ``ln_sigma_inv`` (:math:`\ln\sigma^{-1}`), ``log10_sigma_inv`` @@ -135,6 +137,7 @@ from ._kernels import fits as _k from ._validators import check_finite_positive, less_than, positive from .domain import ( + Calibrated, Domain, DomainError, DomainPolicy, @@ -508,7 +511,7 @@ def _contains(domain: Domain | None, x: FitInputs) -> bool | BoolArray: # The kind # --------------------------------------------------------------------------------- @attrs.frozen(kw_only=True) -class FittingFunction(Model, kind=True): +class FittingFunction(Calibrated, Model, kind=True): r"""A halo mass function fit: the multiplicity function :math:`f(\sigma)`. Subclasses implement :meth:`_fsigma`, and declare, as class variables, the @@ -523,7 +526,8 @@ class FittingFunction(Model, kind=True): valid_domain: ClassVar[Domain] #: Where the fit was calibrated against simulations, with its source; ``None`` for - #: a fit that was not calibrated (e.g. an analytic one). + #: a fit that was not calibrated (e.g. an analytic one). Ask an instance for + #: :attr:`calibration`. calibration_domain: ClassVar[Domain | None] #: The halo mass definition the fit was measured in (see "Mass definitions" in the @@ -704,8 +708,10 @@ def check_valid(self, x: FitInputs) -> None: where=f"{type(self).__name__}'s valid domain", ) - def in_calibration_domain(self, x: FitInputs) -> bool | BoolArray: - """Whether each input is inside :attr:`calibration_domain`. + def in_calibration_domain_kernel(self, x: FitInputs) -> bool | BoolArray: + """Whether each input is inside :attr:`calibration`, at kernel level. + + It checks ``x`` against :attr:`valid_domain`, as :meth:`fsigma_kernel` does. Parameters ---------- @@ -722,8 +728,11 @@ def in_calibration_domain(self, x: FitInputs) -> bool | BoolArray: ------ ValueError If an input needed to check the domain was not given. + DomainError + If any input is outside :attr:`valid_domain`. """ - return _contains(self.calibration_domain, x) + self.check_valid(x) + return _contains(self.calibration, x) @unit_boundary(m=Msun_h) def fsigma( @@ -1038,14 +1047,15 @@ def evaluate_fsigma( """ x = inputs model.check_required(x) - missing = sorted(name for name in _domain_inputs(model.calibration_domain) if not x.given(name)) + missing = sorted(name for name in _domain_inputs(model.calibration) if not x.given(name)) if missing: raise ValueError( f"Checking {type(model).__name__}'s calibration domain needs the input(s) {missing}." ) f = model.fsigma_kernel(x) - inside = np.broadcast_to(model.in_calibration_domain(x), f.shape) - calibration = model.calibration_domain + # fsigma_kernel has checked the valid domain. + inside = np.broadcast_to(_contains(model.calibration, x), f.shape) + calibration = model.calibration # Without a calibration domain everything is inside: this only checks the policy. out = apply_domain_policy( f, diff --git a/src/hmf/core/mass_conversion.py b/src/hmf/core/mass_conversion.py index 94a5cb79..f607202d 100644 --- a/src/hmf/core/mass_conversion.py +++ b/src/hmf/core/mass_conversion.py @@ -149,8 +149,8 @@ def check_definitions(self, source: MassDefinition, target: MassDefinition) -> N @abc.abstractmethod def convert_kernel( self, - m: npt.ArrayLike, *, + m: npt.ArrayLike, z: npt.ArrayLike, omega_m_z: npt.ArrayLike, source: MassDefinition, @@ -185,8 +185,8 @@ def convert_kernel( @unit_boundary(m=Msun_h, returns=(Msun_h, None)) def convert( self, - m: Any, *, + m: Any, z: npt.ArrayLike, omega_m_z: npt.ArrayLike, source: MassDefinition, @@ -213,7 +213,7 @@ def convert( dm_new_dm : numpy.ndarray The Jacobian dM'/dM. """ - out = self.convert_kernel(m, z=z, omega_m_z=omega_m_z, source=source, target=target) + out = self.convert_kernel(m=m, z=z, omega_m_z=omega_m_z, source=source, target=target) return out.m, out.dm_dm @@ -281,30 +281,30 @@ def _concentration_in( model, profile = self.concentration, self.profile own = model.mass_definition if own == definition: - ln_c, slope = model.ln_concentration_and_slope_kernel(m, z=z) + ln_c, slope = model.ln_concentration_and_slope_kernel(m=m, z=z) return np.exp(ln_c), slope, np.array(np.broadcast_to(m, ln_c.shape)) - ln_q = np.log(definition.delta_halo_mean_kernel(omega_m_z)) - np.log( - own.delta_halo_mean_kernel(omega_m_z) + ln_q = np.log(definition.delta_halo_mean_kernel(omega_m_z=omega_m_z)) - np.log( + own.delta_halo_mean_kernel(omega_m_z=omega_m_z) ) ln_m = np.log(m) ln_m_own = np.array(np.broadcast_to(ln_m, np.broadcast_shapes(ln_m.shape, ln_q.shape))) for _ in range(OUTER_NEWTON_STEPS): - ln_c, slope = model.ln_concentration_and_slope_kernel(np.exp(ln_m_own), z=z) + ln_c, slope = model.ln_concentration_and_slope_kernel(m=np.exp(ln_m_own), z=z) ln_ratio, _, dln_m, _ = profile.convert_concentration_kernel( - np.exp(ln_c), slope=slope, ln_q=ln_q + c=np.exp(ln_c), slope=slope, ln_q=ln_q ) ln_m_own = ln_m_own - (ln_m_own + ln_ratio - ln_m) / dln_m m_own = np.exp(ln_m_own) - ln_c, slope = model.ln_concentration_and_slope_kernel(m_own, z=z) + ln_c, slope = model.ln_concentration_and_slope_kernel(m=m_own, z=z) _, c_new, dln_m, dln_c = profile.convert_concentration_kernel( - np.exp(ln_c), slope=slope, ln_q=ln_q + c=np.exp(ln_c), slope=slope, ln_q=ln_q ) return c_new, dln_c / dln_m, m_own def convert_kernel( self, - m: npt.ArrayLike, *, + m: npt.ArrayLike, z: npt.ArrayLike, omega_m_z: npt.ArrayLike, source: MassDefinition, @@ -352,11 +352,11 @@ def convert_kernel( ) self.check_definitions(source, target) c_source, slope, m_concentration = self._concentration_in(m, z, omega_m_z, source) - ln_q = np.log(target.delta_halo_mean_kernel(omega_m_z)) - np.log( - source.delta_halo_mean_kernel(omega_m_z) + ln_q = np.log(target.delta_halo_mean_kernel(omega_m_z=omega_m_z)) - np.log( + source.delta_halo_mean_kernel(omega_m_z=omega_m_z) ) ln_ratio, _, dln_m, _ = self.profile.convert_concentration_kernel( - c_source, slope=slope, ln_q=ln_q + c=c_source, slope=slope, ln_q=ln_q ) ratio = np.exp(ln_ratio) model = self.concentration @@ -365,7 +365,7 @@ def convert_kernel( m=np.array(np.broadcast_to(m * ratio, shape)), dm_dm=np.array(np.broadcast_to(ratio * dln_m, shape)), in_calibration_domain=np.array( - np.broadcast_to(model.in_calibration_domain_kernel(m_concentration, z=z), shape) + np.broadcast_to(model.in_calibration_domain_kernel(m=m_concentration, z=z), shape) ), calibration_description=( "" diff --git a/src/hmf/core/mass_definitions.py b/src/hmf/core/mass_definitions.py index 2949c6e2..41b20a91 100644 --- a/src/hmf/core/mass_definitions.py +++ b/src/hmf/core/mass_definitions.py @@ -28,7 +28,7 @@ :class:`~hmf.core.mass_variance.MassVariance`: a halo of mass m has the radius :math:`r = [3m/(4\pi\Delta_m\bar\rho_{\rm cb,0})]^{1/3}`, with :math:`\bar\rho_{\rm cb,0}` the comoving (z = 0) mean density -(:meth:`~MassDefinition.m_to_r_kernel`); the physical radius is +(:meth:`~MassDefinition.radius_from_m_kernel`); the physical radius is :math:`r/(1+z)`. The comoving halo density is :math:`\Delta_m\bar\rho_{\rm cb,0}` (:meth:`~MassDefinition.halo_density_kernel`). @@ -114,7 +114,8 @@ class MassDefinition(abc.ABC): def name(self) -> str: """A short name, e.g. ``"200m"``, ``"500c"``, ``"vir"``, ``"fof(b=0.2)"``. - :meth:`from_name` reads it back. + :meth:`from_name` reads it back to an equal definition: a number that ``:g`` + would round (e.g. 18 pi^2) is written in full. """ def __str__(self) -> str: @@ -151,8 +152,10 @@ def from_name(cls, name: str) -> MassDefinition: ---------- name ``"200m"`` (:class:`SOMean`), ``"500c"`` (:class:`SOCritical`), ``"vir"`` - (:class:`SOVirial`), ``"fof"`` or ``"fof(b=0.3)"`` (:class:`FOF`) or - ``"self_bound"`` (:class:`SelfBound`); case and spaces are ignored. + (:class:`SOVirial`), ``"so(preferred 200m)"`` (:class:`SOGeneric`, with + the name of its ``preferred`` definition), ``"fof"`` or ``"fof(b=0.3)"`` + (:class:`FOF`) or ``"self_bound"`` (:class:`SelfBound`); case and spaces + are ignored. Returns ------- @@ -170,6 +173,15 @@ def from_name(cls, name: str) -> MassDefinition: return SelfBound() if key == "fof": return FOF() + generic = re.fullmatch(r"so\(preferred(.+)\)", key) + if generic: + preferred = MassDefinition.from_name(generic.group(1)) + if not isinstance(preferred, SOMean | SOCritical | SOVirial): + raise ValueError( + f"Unknown mass definition name {name!r}: the preferred definition of " + "an SOGeneric must be an SOMean, SOCritical or SOVirial." + ) + return SOGeneric(preferred=preferred) fof = re.fullmatch(r"fof\(b=([0-9.eE+-]+)\)", key) if fof: return FOF(linking_length=float(fof.group(1))) @@ -181,7 +193,7 @@ def from_name(cls, name: str) -> MassDefinition: ) raise ValueError( f"Unknown mass definition name {name!r}: use e.g. '200m', '500c', 'vir', " - "'fof', 'fof(b=0.3)' or 'self_bound'." + "'so(preferred 200m)', 'fof', 'fof(b=0.3)' or 'self_bound'." ) # ------------------------------------------------------------------------------ @@ -192,7 +204,7 @@ def _delta_halo_mean(self, omega_m_z: FloatArray) -> FloatArray: """The overdensity relative to the mean, on checked Omega_m(z).""" raise ValueError(f"The mass definition {self} has no overdensity.") - def delta_halo_mean_kernel(self, omega_m_z: npt.ArrayLike) -> FloatArray: + def delta_halo_mean_kernel(self, *, omega_m_z: npt.ArrayLike) -> FloatArray: r"""The halo overdensity relative to the mean density, :math:`\Delta_m`. A pure, elementwise kernel of :math:`\Omega_m(z)` (of CDM + baryons), on plain @@ -251,11 +263,11 @@ def halo_density_kernel( rho = check_extent( "rho_mean0", rho_mean0, _POSITIVE, where=f"{type(self).__name__}.halo_density_kernel" ) - out: FloatArray = self.delta_halo_mean_kernel(omega_m_z) * rho + out: FloatArray = self.delta_halo_mean_kernel(omega_m_z=omega_m_z) * rho return out - def m_to_r_kernel( - self, m: npt.ArrayLike, *, omega_m_z: npt.ArrayLike, rho_mean0: npt.ArrayLike + def radius_from_m_kernel( + self, *, m: npt.ArrayLike, omega_m_z: npt.ArrayLike, rho_mean0: npt.ArrayLike ) -> FloatArray: r"""The comoving radius of haloes of mass m: :math:`[3m/(4\pi\rho_{\rm halo})]^{1/3}`. @@ -280,13 +292,13 @@ def m_to_r_kernel( ValueError For a definition without an overdensity. """ - m = check_extent("m", m, _POSITIVE, where=f"{type(self).__name__}.m_to_r_kernel") + m = check_extent("m", m, _POSITIVE, where=f"{type(self).__name__}.radius_from_m_kernel") rho = self.halo_density_kernel(omega_m_z=omega_m_z, rho_mean0=rho_mean0) out: FloatArray = np.cbrt(3.0 * m / (4.0 * np.pi * rho)) return out - def r_to_m_kernel( - self, r: npt.ArrayLike, *, omega_m_z: npt.ArrayLike, rho_mean0: npt.ArrayLike + def m_from_radius_kernel( + self, *, r: npt.ArrayLike, omega_m_z: npt.ArrayLike, rho_mean0: npt.ArrayLike ) -> FloatArray: r"""The mass of haloes of comoving radius r: :math:`4\pi r^3\rho_{\rm halo}/3`. @@ -311,7 +323,7 @@ def r_to_m_kernel( ValueError For a definition without an overdensity. """ - r = check_extent("r", r, _POSITIVE, where=f"{type(self).__name__}.r_to_m_kernel") + r = check_extent("r", r, _POSITIVE, where=f"{type(self).__name__}.m_from_radius_kernel") rho = self.halo_density_kernel(omega_m_z=omega_m_z, rho_mean0=rho_mean0) out: FloatArray = 4.0 * np.pi * r**3 * rho / 3.0 return out @@ -321,7 +333,7 @@ def r_to_m_kernel( # ------------------------------------------------------------------------------ @unit_boundary() - def delta_halo_mean(self, omega_m_z: npt.ArrayLike) -> FloatArray: + def delta_halo_mean(self, *, omega_m_z: npt.ArrayLike) -> FloatArray: r"""The halo overdensity relative to the mean density of CDM + baryons. See :meth:`delta_halo_mean_kernel`. @@ -336,7 +348,7 @@ def delta_halo_mean(self, omega_m_z: npt.ArrayLike) -> FloatArray: numpy.float64 or numpy.ndarray :math:`\Delta_m`. """ - return self.delta_halo_mean_kernel(omega_m_z) + return self.delta_halo_mean_kernel(omega_m_z=omega_m_z) @unit_boundary(rho_mean0=rho_unit, returns=rho_unit) def halo_density(self, *, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatArray: @@ -357,8 +369,8 @@ def halo_density(self, *, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatArra return self.halo_density_kernel(omega_m_z=omega_m_z, rho_mean0=rho_mean0) @unit_boundary(m=Msun_h, rho_mean0=rho_unit, returns=Mpc_h) - def m_to_r(self, m: Any, *, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatArray: - r"""The comoving radius of haloes of mass m (see :meth:`m_to_r_kernel`). + def radius_from_m(self, *, m: Any, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatArray: + r"""The comoving radius of haloes of mass m (see :meth:`radius_from_m_kernel`). Parameters ---------- @@ -374,11 +386,11 @@ def m_to_r(self, m: Any, *, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatAr Quantity Comoving radii in Mpc/h. """ - return self.m_to_r_kernel(m, omega_m_z=omega_m_z, rho_mean0=rho_mean0) + return self.radius_from_m_kernel(m=m, omega_m_z=omega_m_z, rho_mean0=rho_mean0) @unit_boundary(r=Mpc_h, rho_mean0=rho_unit, returns=Msun_h) - def r_to_m(self, r: Any, *, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatArray: - r"""The mass of haloes of comoving radius r (see :meth:`r_to_m_kernel`). + def m_from_radius(self, *, r: Any, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatArray: + r"""The mass of haloes of comoving radius r (see :meth:`m_from_radius_kernel`). Parameters ---------- @@ -394,7 +406,7 @@ def r_to_m(self, r: Any, *, omega_m_z: npt.ArrayLike, rho_mean0: Any) -> FloatAr Quantity Masses in Msun/h, in this definition. """ - return self.r_to_m_kernel(r, omega_m_z=omega_m_z, rho_mean0=rho_mean0) + return self.m_from_radius_kernel(r=r, omega_m_z=omega_m_z, rho_mean0=rho_mean0) @attrs.frozen(kw_only=True) @@ -402,6 +414,12 @@ class SphericalOverdensity(MassDefinition): """Base class of the spherical-overdensity definitions.""" +def _number(value: float) -> str: + """``value`` in a name: as ``:g`` (``200``) where that is exact, else in full.""" + short = f"{value:g}" + return short if float(short) == value else repr(value) + + def _overdensity(default: float) -> Any: return attrs.field(default=default, converter=float, validator=positive) @@ -423,7 +441,7 @@ class SOMean(SphericalOverdensity): @property def name(self) -> str: """``"m"``, e.g. ``"200m"``.""" - return f"{self.overdensity:g}m" + return f"{_number(self.overdensity)}m" def _delta_halo_mean(self, omega_m_z: FloatArray) -> FloatArray: return np.full_like(omega_m_z, self.overdensity) @@ -446,7 +464,7 @@ class SOCritical(SphericalOverdensity): @property def name(self) -> str: """``"c"``, e.g. ``"500c"``.""" - return f"{self.overdensity:g}c" + return f"{_number(self.overdensity)}c" def _delta_halo_mean(self, omega_m_z: FloatArray) -> FloatArray: out: FloatArray = self.overdensity / omega_m_z @@ -553,7 +571,7 @@ class FOF(MassDefinition): @property def name(self) -> str: """``"fof(b=)"``.""" - return f"fof(b={self.linking_length:g})" + return f"fof(b={_number(self.linking_length)})" def _delta_halo_mean(self, omega_m_z: FloatArray) -> FloatArray: return np.full_like(omega_m_z, 9 / (2 * math.pi * self.linking_length**3)) diff --git a/src/hmf/core/mass_function.py b/src/hmf/core/mass_function.py index b8afc0c6..0fc98a00 100644 --- a/src/hmf/core/mass_function.py +++ b/src/hmf/core/mass_function.py @@ -522,7 +522,7 @@ def _conversion(self, m: npt.ArrayLike, z: npt.ArrayLike, *, to_fit: bool) -> Co assert self.mass_conversion is not None out, fit = self.resolved_mass_definition, self.fit_mass_definition return self.mass_conversion.convert_kernel( - m, + m=m, z=z, omega_m_z=omega_m(self.cosmology, z, "cb"), source=out if to_fit else fit, @@ -646,7 +646,7 @@ def _fit_inputs( """The inputs of the fit at (m, z), from sigma and its slope there.""" omega_m_z = omega_m(self.cosmology, z, "cb") delta_halo = ( - self.fit_mass_definition.delta_halo_mean_kernel(omega_m_z) + self.fit_mass_definition.delta_halo_mean_kernel(omega_m_z=omega_m_z) if self._needs_delta_halo else None ) @@ -969,11 +969,11 @@ def in_calibration_domain_kernel(self, *, m: npt.ArrayLike, z: npt.ArrayLike) -> if not self.converts_masses: x, _ = self._inputs_at(m, z) shape = np.broadcast_shapes(np.shape(x.sigma), np.shape(m), np.shape(z)) - return np.array(np.broadcast_to(self.fit.in_calibration_domain(x), shape)) + return np.array(np.broadcast_to(self.fit.in_calibration_domain_kernel(x), shape)) conversion = self._conversion(m, z, to_fit=True) x, _ = self._inputs_at(conversion.m, z) shape = np.broadcast_shapes(np.shape(x.sigma), np.shape(conversion.m)) - inside = self.fit.in_calibration_domain(x) & conversion.in_calibration_domain + inside = self.fit.in_calibration_domain_kernel(x) & conversion.in_calibration_domain return np.array(np.broadcast_to(inside, shape)) def m_from_sigma_kernel(self, *, sigma: npt.ArrayLike, z: npt.ArrayLike) -> FloatArray: diff --git a/src/hmf/core/profiles.py b/src/hmf/core/profiles.py index a67fb72e..d4a2d328 100644 --- a/src/hmf/core/profiles.py +++ b/src/hmf/core/profiles.py @@ -38,7 +38,7 @@ from . import _references as refs from ._kernels import mass_conversion as _k -from .domain import Domain, Interval, check_extent +from .domain import Calibrated, Domain, Interval, check_extent from .model import Model from .units import unit_boundary @@ -50,7 +50,7 @@ @attrs.frozen(kw_only=True) -class DensityProfile(Model, kind=True): +class DensityProfile(Calibrated, Model, kind=True): r"""The kind of halo density profiles: their enclosed mass :math:`\mu(x)`. Subclasses implement :meth:`_ln_mu_and_slope`. See the module documentation. @@ -76,7 +76,7 @@ def _start(self, c: FloatArray, ln_q: FloatArray) -> FloatArray: _, psi = self._ln_mu_and_slope(c) return _k.power_law_start(c, ln_q, psi) - def ln_mu_and_slope_kernel(self, x: npt.ArrayLike) -> tuple[FloatArray, FloatArray]: + def ln_mu_and_slope_kernel(self, *, x: npt.ArrayLike) -> tuple[FloatArray, FloatArray]: r""":math:`\ln\mu(x)` and :math:`\psi(x) = d\ln\mu/d\ln x`, at kernel level. Parameters @@ -100,7 +100,7 @@ def ln_mu_and_slope_kernel(self, x: npt.ArrayLike) -> tuple[FloatArray, FloatArr return self._ln_mu_and_slope(x) def convert_concentration_kernel( - self, c: npt.ArrayLike, *, slope: npt.ArrayLike, ln_q: npt.ArrayLike + self, *, c: npt.ArrayLike, slope: npt.ArrayLike, ln_q: npt.ArrayLike ) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray]: r"""Convert haloes of concentration c to a definition of q times the halo density. @@ -139,7 +139,7 @@ def convert_concentration_kernel( ) @unit_boundary() - def mu(self, x: npt.ArrayLike) -> FloatArray: + def mu(self, *, x: npt.ArrayLike) -> FloatArray: r"""The enclosed mass in units of :math:`4\pi\rho_s r_s^3`, :math:`\mu(x)`. Parameters @@ -157,7 +157,7 @@ def mu(self, x: npt.ArrayLike) -> FloatArray: DomainError If an x is not finite and > 0. """ - return np.exp(self.ln_mu_and_slope_kernel(x)[0]) + return np.exp(self.ln_mu_and_slope_kernel(x=x)[0]) @attrs.frozen(kw_only=True) diff --git a/tests/core/test_core_concentration.py b/tests/core/test_core_concentration.py index 39d7137d..863861e3 100644 --- a/tests/core/test_core_concentration.py +++ b/tests/core/test_core_concentration.py @@ -7,10 +7,10 @@ import pytest from hmf.core import _references as refs -from hmf.core._kernels import mass_conversion as kern from hmf.core.concentration import DUFFY08_PARAMETERS, ConcentrationModel, Duffy08 from hmf.core.domain import DomainError from hmf.core.mass_definitions import SOCritical, SOMean, SOVirial +from hmf.core.profiles import NFW from hmf.core.units import Msun_h, UnitBoundaryError M = np.geomspace(1e11, 1e15, 9) @@ -46,9 +46,9 @@ def test_power_law_at_the_pivot(): """C = A (M / 2e12)^B (1 + z)^C: A at the pivot today, and slope B everywhere.""" for model in (Duffy08(), Duffy08(definition="vir", sample="relaxed")): a, b, c = model.parameters - assert model.concentration(2e12 * Msun_h, z=0.0) == pytest.approx(a, rel=1e-14) - assert model.concentration(2e12 * Msun_h, z=1.0) == pytest.approx(a * 2**c, rel=1e-14) - ln_c, slope = model.ln_concentration_and_slope_kernel(M, z=0.5) + assert model.concentration(m=2e12 * Msun_h, z=0.0) == pytest.approx(a, rel=1e-14) + assert model.concentration(m=2e12 * Msun_h, z=1.0) == pytest.approx(a * 2**c, rel=1e-14) + ln_c, slope = model.ln_concentration_and_slope_kernel(m=M, z=0.5) np.testing.assert_allclose(np.diff(ln_c) / np.diff(np.log(M)), b, rtol=1e-12) np.testing.assert_array_equal(slope, b) @@ -59,14 +59,14 @@ def test_concentrations_are_physical(): c_200 ~ 4-8 at 1e11-1e15 today, c_mean nearly twice c_200, relaxed haloes about 10% more concentrated, all falling with z. """ - c200 = Duffy08(redshifts="0").concentration(M * Msun_h, z=0.0) + c200 = Duffy08(redshifts="0").concentration(m=M * Msun_h, z=0.0) assert np.all((c200 > 3) & (c200 < 9)) - cmean = Duffy08(definition="200m", redshifts="0").concentration(M * Msun_h, z=0.0) + cmean = Duffy08(definition="200m", redshifts="0").concentration(m=M * Msun_h, z=0.0) assert np.all((cmean / c200 > 1.6) & (cmean / c200 < 2.0)) - relaxed = Duffy08(sample="relaxed", redshifts="0").concentration(M * Msun_h, z=0.0) + relaxed = Duffy08(sample="relaxed", redshifts="0").concentration(m=M * Msun_h, z=0.0) assert np.all((relaxed / c200 > 1.05) & (relaxed / c200 < 1.25)) z = np.array([0.0, 1.0, 2.0]) - assert np.all(np.diff(Duffy08().concentration(1e12 * Msun_h, z=z)) < 0) + assert np.all(np.diff(Duffy08().concentration(m=1e12 * Msun_h, z=z)) < 0) @pytest.mark.parametrize("sample", ["full", "relaxed"]) @@ -82,12 +82,14 @@ def test_duffy_relations_are_nfw_conversions_of_each_other(sample, definition): omega_m = 0.258 c200 = Duffy08(sample=sample, redshifts="0") other = Duffy08(definition=definition, sample=sample, redshifts="0") - ln_c, slope = c200.ln_concentration_and_slope_kernel(M, z=0.0) - ln_q = np.log(other.mass_definition.delta_halo_mean_kernel(omega_m)) - np.log( - SOCritical().delta_halo_mean_kernel(omega_m) + ln_c, slope = c200.ln_concentration_and_slope_kernel(m=M, z=0.0) + ln_q = np.log(other.mass_definition.delta_halo_mean_kernel(omega_m_z=omega_m)) - np.log( + SOCritical().delta_halo_mean_kernel(omega_m_z=omega_m) ) - ln_ratio, c_new, _, _ = kern.nfw_convert(np.exp(ln_c), slope, ln_q) - expected = other.concentration(M * np.exp(ln_ratio) * Msun_h, z=0.0) + ln_ratio, c_new, _, _ = NFW().convert_concentration_kernel( + c=np.exp(ln_c), slope=slope, ln_q=ln_q + ) + expected = other.concentration(m=M * np.exp(ln_ratio) * Msun_h, z=0.0) np.testing.assert_allclose(c_new, expected, rtol=0.015) @@ -97,16 +99,31 @@ def test_domains(): assert not Duffy08().calibration.contains(m=1e16 * Msun_h, z=0.0) assert "Duffy et al. 2008" in Duffy08().calibration.source assert Duffy08.calibration_domain is Duffy08().calibration + # The class's is the default set's; an instance's is its own set's. + assert Duffy08(redshifts="0").calibration is not Duffy08.calibration_domain + np.testing.assert_array_equal( + Duffy08(redshifts="0").in_calibration_domain_kernel(m=1e12, z=np.array([0.0, 1.5])), + [True, False], + ) np.testing.assert_array_equal( - Duffy08().in_calibration_domain_kernel(np.array([1e10, 1e12, 1e16]), z=0.0), + Duffy08().in_calibration_domain_kernel(m=np.array([1e10, 1e12, 1e16]), z=0.0), [False, True, False], ) with pytest.raises(DomainError, match="m in"): - Duffy08().ln_concentration_and_slope_kernel(np.array([1e12, -1.0]), z=0.0) + Duffy08().ln_concentration_and_slope_kernel(m=np.array([1e12, -1.0]), z=0.0) with pytest.raises(DomainError, match="z in"): - Duffy08().concentration(1e12 * Msun_h, z=-0.5) + Duffy08().concentration(m=1e12 * Msun_h, z=-0.5) with pytest.raises(UnitBoundaryError): - Duffy08().concentration(1e12, z=0.0) + Duffy08().concentration(m=1e12, z=0.0) + + +def test_arguments_are_keyword_only(): + with pytest.raises(TypeError, match="positional"): + Duffy08().concentration(1e12 * Msun_h, z=0.0) + with pytest.raises(TypeError, match="positional"): + Duffy08().ln_concentration_and_slope_kernel(1e12, z=0.0) + with pytest.raises(TypeError, match="positional"): + Duffy08().in_calibration_domain_kernel(1e12, z=0.0) def test_invalid_options(): @@ -125,7 +142,7 @@ def test_values_hash_and_pickle(): def test_scalar_rule_and_broadcasting(): - c = Duffy08().concentration(1e12 * Msun_h, z=0.0) + c = Duffy08().concentration(m=1e12 * Msun_h, z=0.0) assert isinstance(c, float) - shape = Duffy08().concentration(M[None, :] * Msun_h, z=np.array([0.0, 1.0])[:, None]).shape + shape = Duffy08().concentration(m=M[None, :] * Msun_h, z=np.array([0.0, 1.0])[:, None]).shape assert shape == (2, M.size) diff --git a/tests/core/test_core_domain.py b/tests/core/test_core_domain.py index 90fec28e..51e64359 100644 --- a/tests/core/test_core_domain.py +++ b/tests/core/test_core_domain.py @@ -420,3 +420,70 @@ def test_warn_once_forgets_collected_owners(): def test_warn_once_without_weak_references_warns_every_time(): caught = _record(lambda: [warn_once(1.5, "k", "m") for _ in range(2)]) assert len(caught) == 2 + + +# --------------------------------------------------------------------------------- +# The calibration domain of a model instance +# --------------------------------------------------------------------------------- + + +def _calibrated_kinds(): + from hmf.core.concentration import ConcentrationModel + from hmf.core.filters import Filter + from hmf.core.fits import FittingFunction + from hmf.core.growth_models import GrowthModel + from hmf.core.profiles import DensityProfile + from hmf.core.transfer_models import TransferModel + + return [FittingFunction, ConcentrationModel, DensityProfile, Filter, TransferModel, GrowthModel] + + +def test_every_kind_with_a_calibration_domain_is_calibrated(): + """Every such kind answers ``calibration`` on its instances: by default the class's.""" + from hmf.core.domain import Calibrated + + for kind in _calibrated_kinds(): + assert issubclass(kind, Calibrated), kind + for cls in kind.get_models().values(): + if cls.calibration is not Calibrated.calibration: + continue # it depends on the instance (Duffy08) + try: + model = cls() + except (TypeError, ValueError): + continue # a model that needs arguments (e.g. a table) + assert model.calibration is cls.calibration_domain, cls + + +def test_calibration_can_depend_on_the_instance(): + """Duffy08 at z = 0 alone is not calibrated at z = 1.5; the class default (z = 0-2) is.""" + from hmf.core.concentration import Duffy08 + + z0 = Duffy08(redshifts="0") + assert Duffy08.calibration_domain.contains(m=1e12 * Msun_h, z=1.5) + assert not z0.calibration.contains(m=1e12 * Msun_h, z=1.5) + assert z0.calibration.contains(m=1e12 * Msun_h, z=0.0) + + +def test_evaluate_fsigma_uses_the_instance_calibration(): + """A fit whose calibration depends on a parameter is masked by its own.""" + import attrs + + from hmf.core import fits + + @attrs.frozen(kw_only=True) + class _Narrow(fits.PS, abstract=True): + sigma_min: float = 0.5 + + @property + def calibration(self): + return Domain({"sigma": (self.sigma_min, None)}, source="A test.") + + x = fits.PS().inputs_from_arrays(np.array([0.3, 1.0]), delta_c=1.686) + # PS has no calibration domain: everything is inside. + assert fits.PS().calibration is None + np.testing.assert_array_equal(fits.PS().in_calibration_domain_kernel(x), [True, True]) + for sigma_min, inside in ((0.5, [False, True]), (0.1, [True, True])): + model = _Narrow(sigma_min=sigma_min) + np.testing.assert_array_equal(model.in_calibration_domain_kernel(x), inside) + out = fits.evaluate_fsigma(model, x, policy="mask") + np.testing.assert_array_equal(np.isnan(out.fsigma), np.logical_not(inside)) diff --git a/tests/core/test_core_fits.py b/tests/core/test_core_fits.py index e3918034..9eeff875 100644 --- a/tests/core/test_core_fits.py +++ b/tests/core/test_core_fits.py @@ -869,7 +869,7 @@ def test_fsigma_is_dimensionless_and_accepts_equivalent_mass_units(): other = [1e12, 1e13] * (u.Msun / cu.littleh) np.testing.assert_array_equal(model.fsigma([1.0, 2.0], m=other), f) x = model.inputs([1.0, 2.0], z=0.0, m=[1e15, 1e16] * Msun_h) - np.testing.assert_array_equal(model.in_calibration_domain(x), [True, False]) + np.testing.assert_array_equal(model.in_calibration_domain_kernel(x), [True, False]) def test_masses_without_h_need_a_hubble_constant(): @@ -1018,7 +1018,7 @@ def test_every_fit_has_a_delta_halo_if_it_needs_one(cls): """A fit that needs delta_halo can get it from its own measured mass definition.""" if "delta_halo" not in cls.domain_inputs(): return - delta = cls.measured_mass_definition.delta_halo_mean_kernel(_OMEGA_M_Z) + delta = cls.measured_mass_definition.delta_halo_mean_kernel(omega_m_z=_OMEGA_M_Z) assert np.all(np.isfinite(delta) & (delta > 0)) diff --git a/tests/core/test_core_kernel_safety.py b/tests/core/test_core_kernel_safety.py index 0e41c654..04b265bd 100644 --- a/tests/core/test_core_kernel_safety.py +++ b/tests/core/test_core_kernel_safety.py @@ -394,56 +394,64 @@ def _ups(model="EH"): ("UnnormalisedPower.ln_power_kernel", lambda: _ups().ln_power_kernel(np.array([np.inf]))), ("UnnormalisedPower.power_kernel", lambda: _ups().power_kernel(np.array([1.0, 0.0]))), ("FittingFunction.fsigma_kernel", _fsigma_kernel), + ( + "FittingFunction.in_calibration_domain_kernel", + lambda: fits.Tinker08().in_calibration_domain_kernel( + FitInputs(sigma=[1.0, -1.0], z=0.0, delta_halo=200.0) + ), + ), ( "Filter.window_derivatives_kernel", lambda: TopHat().window_derivatives_kernel(np.array([1.0, -0.5])), ), ( "MassDefinition.delta_halo_mean_kernel", - lambda: fits.Tinker08.measured_mass_definition.delta_halo_mean_kernel([0.3, -0.1]), + lambda: fits.Tinker08.measured_mass_definition.delta_halo_mean_kernel( + omega_m_z=[0.3, -0.1] + ), ), ( "MassDefinition.halo_density_kernel", lambda: SOMean().halo_density_kernel(omega_m_z=0.3, rho_mean0=[1.0, 0.0]), ), ( - "MassDefinition.m_to_r_kernel", - lambda: SOMean().m_to_r_kernel([1e12, -1.0], omega_m_z=0.3, rho_mean0=1e11), + "MassDefinition.radius_from_m_kernel", + lambda: SOMean().radius_from_m_kernel(m=[1e12, -1.0], omega_m_z=0.3, rho_mean0=1e11), ), ( - "MassDefinition.r_to_m_kernel", - lambda: SOMean().r_to_m_kernel([1.0, np.nan], omega_m_z=0.3, rho_mean0=1e11), + "MassDefinition.m_from_radius_kernel", + lambda: SOMean().m_from_radius_kernel(r=[1.0, np.nan], omega_m_z=0.3, rho_mean0=1e11), ), ( "ConcentrationModel.ln_concentration_and_slope_kernel", - lambda: Duffy08().ln_concentration_and_slope_kernel([1e12, 0.0], z=0.0), + lambda: Duffy08().ln_concentration_and_slope_kernel(m=[1e12, 0.0], z=0.0), ), ( "ConcentrationModel.in_calibration_domain_kernel", - lambda: Duffy08().in_calibration_domain_kernel(1e12, z=[0.0, -1.0]), + lambda: Duffy08().in_calibration_domain_kernel(m=1e12, z=[0.0, -1.0]), ), ( "DensityProfile.ln_mu_and_slope_kernel", - lambda: NFW().ln_mu_and_slope_kernel([1.0, 0.0]), + lambda: NFW().ln_mu_and_slope_kernel(x=[1.0, 0.0]), ), ( "DensityProfile.convert_concentration_kernel", - lambda: NFW().convert_concentration_kernel([5.0, -1.0], slope=-0.1, ln_q=0.5), + lambda: NFW().convert_concentration_kernel(c=[5.0, -1.0], slope=-0.1, ln_q=0.5), ), ( "DensityProfile.convert_concentration_kernel ln_q", - lambda: NFW().convert_concentration_kernel(5.0, slope=-0.1, ln_q=np.inf), + lambda: NFW().convert_concentration_kernel(c=5.0, slope=-0.1, ln_q=np.inf), ), ( "ProfileConversion.convert_kernel", lambda: ProfileConversion().convert_kernel( - [1e12, np.inf], z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical() + m=[1e12, np.inf], z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical() ), ), ( "ProfileConversion.convert_kernel omega_m_z", lambda: ProfileConversion().convert_kernel( - 1e12, z=0.0, omega_m_z=0.0, source=SOMean(), target=SOCritical() + m=1e12, z=0.0, omega_m_z=0.0, source=SOMean(), target=SOCritical() ), ), ("FittingFunction.modify_dndm_kernel", lambda: _modify_dndm(fits.PS(), m=-1.0)), @@ -526,7 +534,7 @@ def test_kernel_checks_do_not_change_values(): assert np.all(np.isfinite(out)) mdef = fits.Tinker08.measured_mass_definition om = np.array([0.3, 0.5]) - np.testing.assert_array_equal(mdef.delta_halo_mean_kernel(om), np.full(2, 200.0)) + np.testing.assert_array_equal(mdef.delta_halo_mean_kernel(omega_m_z=om), np.full(2, 200.0)) np.testing.assert_array_equal(n_eff_kernel([-0.5, 0.0]), [0.0, -3.0]) diff --git a/tests/core/test_core_mass_conversion.py b/tests/core/test_core_mass_conversion.py index f96109b9..977cfce1 100644 --- a/tests/core/test_core_mass_conversion.py +++ b/tests/core/test_core_mass_conversion.py @@ -46,31 +46,52 @@ def _ln_q(target, source, om=OMEGA_M0): - return np.log(target.delta_halo_mean_kernel(om)) - np.log(source.delta_halo_mean_kernel(om)) + return np.log(target.delta_halo_mean_kernel(omega_m_z=om)) - np.log( + source.delta_halo_mean_kernel(omega_m_z=om) + ) + + +def _ln_mu(x): + """The NFW ln mu(x), through the profile.""" + return NFW().ln_mu_and_slope_kernel(x=x)[0] + + +def _ln_mean_density(x): + """Ln F(x) = ln[mu(x) / x^3], the NFW mean enclosed density over 3 rho_s.""" + return _ln_mu(x) - 3 * np.log(x) + + +def _solve(c, ln_q): + """The NFW concentration c' with F(c') = q F(c), through the profile.""" + return NFW().convert_concentration_kernel(c=c, slope=0.0, ln_q=ln_q)[1] + + +def _nfw_convert(c, slope, ln_q): + return NFW().convert_concentration_kernel(c=c, slope=slope, ln_q=ln_q) # --------------------------------------------------------------------------------- -# Kernels +# The NFW profile's kernels # --------------------------------------------------------------------------------- def test_ln_mu_closed_form_and_series(): """mu(1) = ln 2 - 1/2; mu ~ x^2 / 2 as x -> 0; the series joins the closed form.""" - assert kern.nfw_ln_mu(1.0) == pytest.approx(math.log(math.log(2) - 0.5), rel=1e-15) + assert _ln_mu(1.0) == pytest.approx(math.log(math.log(2) - 0.5), rel=1e-15) x = np.array([1e-8, 1e-6, 1e-4]) np.testing.assert_allclose( - np.exp(kern.nfw_ln_mu(x)), x**2 / 2 * (1 - 4 * x / 3 + 1.5 * x**2), rtol=1e-11 + np.exp(_ln_mu(x)), x**2 / 2 * (1 - 4 * x / 3 + 1.5 * x**2), rtol=1e-11 ) below, above = np.nextafter(0.01, 0), 0.01 - assert kern.nfw_ln_mu(below) == pytest.approx(float(kern.nfw_ln_mu(above)), abs=1e-13) + assert _ln_mu(below) == pytest.approx(float(_ln_mu(above)), abs=1e-13) # psi = dln mu / dln x, from 2 at small x to 0 at large x. x = np.geomspace(1e-6, 1e6, 50) - psi = kern.nfw_ln_mu_slope(x) + psi = NFW().ln_mu_and_slope_kernel(x=x)[1] assert psi[0] == pytest.approx(2.0, rel=1e-5) assert np.all(np.diff(psi) < 0) assert psi[-1] < 0.1 h = 1e-6 - numeric = (kern.nfw_ln_mu(x * math.exp(h)) - kern.nfw_ln_mu(x * math.exp(-h))) / (2 * h) + numeric = (_ln_mu(x * math.exp(h)) - _ln_mu(x * math.exp(-h))) / (2 * h) np.testing.assert_allclose(psi, numeric, rtol=1e-7, atol=1e-9) @@ -79,8 +100,8 @@ def test_ln_mu_closed_form_and_series(): def test_root_solve_reaches_round_off_over_a_wide_range(): - c_new = kern.nfw_concentration_ratio_solve(C_GRID, LN_Q_GRID) - residual = kern.nfw_ln_mean_density(c_new) - kern.nfw_ln_mean_density(C_GRID) - LN_Q_GRID + c_new = _solve(C_GRID, LN_Q_GRID) + residual = _ln_mean_density(c_new) - _ln_mean_density(C_GRID) - LN_Q_GRID assert np.all(np.isfinite(c_new)) assert np.abs(residual).max() < 2e-13 # A denser definition has a smaller concentration (F decreases monotonically). @@ -90,13 +111,13 @@ def test_root_solve_reaches_round_off_over_a_wide_range(): def test_root_solve_is_elementwise_and_exact_for_the_same_density(): c = np.geomspace(0.5, 50, 101) ln_q = np.linspace(-2, 2, 101) - full = kern.nfw_concentration_ratio_solve(c, ln_q) + full = _nfw_convert(c, -0.1, ln_q) for i in (0, 17, 100): - np.testing.assert_array_equal( - kern.nfw_concentration_ratio_solve(c[i : i + 1], ln_q[i : i + 1]), full[i : i + 1] - ) - np.testing.assert_array_equal(kern.nfw_concentration_ratio_solve(c, 0.0), c) - ln_ratio, c_new, dln_m, _ = kern.nfw_convert(c, -0.1, 0.0) + part = _nfw_convert(c[i : i + 1], -0.1, ln_q[i : i + 1]) + for got, expected in zip(part, full, strict=True): + np.testing.assert_array_equal(got, expected[i : i + 1]) + np.testing.assert_array_equal(_solve(c, 0.0), c) + ln_ratio, c_new, dln_m, _ = _nfw_convert(c, -0.1, 0.0) np.testing.assert_array_equal(ln_ratio, 0.0) np.testing.assert_array_equal(c_new, c) np.testing.assert_array_equal(dln_m, 1.0) @@ -108,7 +129,7 @@ def test_root_solve_agrees_with_v3_brentq(): c = np.geomspace(0.5, 60, 25) for ln_q in (math.log(0.3), math.log(200 / 333), math.log(1600 / 200), math.log(5.0)): - rho_s = c**3 / (3 * np.exp(kern.nfw_ln_mu(c))) # for a halo density of 1 + rho_s = c**3 / (3 * np.exp(_ln_mu(c))) # for a halo density of 1 v3 = np.array( [ # A guess within v3's first bracket (x5) of the root: q^(-1/3) c. @@ -116,14 +137,14 @@ def test_root_solve_agrees_with_v3_brentq(): for r, cc in zip(rho_s, c, strict=True) ] ) - np.testing.assert_allclose(kern.nfw_concentration_ratio_solve(c, ln_q), v3, rtol=1e-10) + np.testing.assert_allclose(_solve(c, ln_q), v3, rtol=1e-10) def _hu_kravtsov_mass_error(c, ratio): """|M_HK / M - 1| for the mass ratio M_h / M_v = (D_h / D_v) (c_h / c_v)^3 (eq. C10).""" ln_q = np.log(ratio) c_hk = kern.hu_kravtsov_inverse(kern.nfw_ln_mean_density(c) + ln_q) - exact = np.exp(kern.nfw_ln_mu(kern.nfw_concentration_ratio_solve(c, ln_q)) - kern.nfw_ln_mu(c)) + exact = np.exp(_ln_mu(_solve(c, ln_q)) - _ln_mu(c)) return np.abs(ratio * (c_hk / c) ** 3 / exact - 1) @@ -150,11 +171,11 @@ def test_jacobian_is_the_derivative_of_the_mapping(): c = np.geomspace(1, 30, 20) slope = -0.1 ln_q = math.log(2.5) - _, _, dln_m, dln_c = kern.nfw_convert(c, slope, ln_q) + _, _, dln_m, dln_c = _nfw_convert(c, slope, ln_q) h = 1e-5 # c(M) = c0 (M / M0)^slope: a step h in ln M is a step slope h in ln c. - lo = kern.nfw_convert(c * math.exp(-slope * h), slope, ln_q) - hi = kern.nfw_convert(c * math.exp(slope * h), slope, ln_q) + lo = _nfw_convert(c * math.exp(-slope * h), slope, ln_q) + hi = _nfw_convert(c * math.exp(slope * h), slope, ln_q) np.testing.assert_allclose(dln_m, 1 + (hi[0] - lo[0]) / (2 * h), rtol=1e-8) np.testing.assert_allclose(dln_c, (np.log(hi[1]) - np.log(lo[1])) / (2 * h), rtol=1e-7) @@ -181,7 +202,7 @@ def test_registered_and_cited(nfw): @pytest.mark.parametrize("definition", [*SO_DEFINITIONS, FOF()], ids=str) def test_converting_to_the_same_definition_is_the_identity(nfw, definition): - out = nfw.convert_kernel(M, z=0.5, omega_m_z=0.6, source=definition, target=definition) + out = nfw.convert_kernel(m=M, z=0.5, omega_m_z=0.6, source=definition, target=definition) np.testing.assert_array_equal(out.m, M) np.testing.assert_array_equal(out.dm_dm, 1.0) @@ -190,7 +211,7 @@ def test_virial_is_so_mean_18pi2_in_einstein_de_sitter(nfw): """Omega_m(z) = 1: the same halo density, so the conversion is exactly the identity.""" eds = SOMean(overdensity=18 * math.pi**2) for source, target in ((SOVirial(), eds), (eds, SOVirial())): - out = nfw.convert_kernel(M, z=1.0, omega_m_z=1.0, source=source, target=target) + out = nfw.convert_kernel(m=M, z=1.0, omega_m_z=1.0, source=source, target=target) np.testing.assert_array_equal(out.m, M) np.testing.assert_array_equal(out.dm_dm, 1.0) @@ -205,8 +226,8 @@ def test_round_trip(a, b, concentration): """A -> B -> A returns the masses to 1e-13, and the Jacobians are reciprocal.""" conv = ProfileConversion(concentration=concentration) z, om = 0.7, float(omega_m(Planck18, 0.7)) - there = conv.convert_kernel(M, z=z, omega_m_z=om, source=a, target=b) - back = conv.convert_kernel(there.m, z=z, omega_m_z=om, source=b, target=a) + there = conv.convert_kernel(m=M, z=z, omega_m_z=om, source=a, target=b) + back = conv.convert_kernel(m=there.m, z=z, omega_m_z=om, source=b, target=a) np.testing.assert_allclose(back.m, M, rtol=1e-13) np.testing.assert_allclose(there.dm_dm * back.dm_dm, 1.0, rtol=1e-12) @@ -223,10 +244,10 @@ def test_round_trip(a, b, concentration): def test_jacobian_matches_a_finite_difference(source, concentration): conv = ProfileConversion(concentration=concentration) target = SOVirial() - out = conv.convert_kernel(M, z=0.3, omega_m_z=0.4, source=source, target=target) + out = conv.convert_kernel(m=M, z=0.3, omega_m_z=0.4, source=source, target=target) h = 1e-5 - hi = conv.convert_kernel(M * math.exp(h), z=0.3, omega_m_z=0.4, source=source, target=target) - lo = conv.convert_kernel(M * math.exp(-h), z=0.3, omega_m_z=0.4, source=source, target=target) + hi = conv.convert_kernel(m=M * math.exp(h), z=0.3, omega_m_z=0.4, source=source, target=target) + lo = conv.convert_kernel(m=M * math.exp(-h), z=0.3, omega_m_z=0.4, source=source, target=target) numeric = (np.log(hi.m) - np.log(lo.m)) / (2 * h) * out.m / M np.testing.assert_allclose(out.dm_dm, numeric, rtol=1e-8) @@ -235,10 +256,10 @@ def test_concentration_is_converted_to_the_source_definition(): """A halo converted to the relation's definition has the relation's concentration.""" conv = ProfileConversion(concentration=Duffy08()) # in 200c z, om = 0.5, 0.55 - out = conv.convert_kernel(M, z=z, omega_m_z=om, source=SOMean(), target=SOCritical()) + out = conv.convert_kernel(m=M, z=z, omega_m_z=om, source=SOMean(), target=SOCritical()) # Through the kernel: c in 200m from c in 200c, at the 200c masses found. - ln_c, slope = Duffy08().ln_concentration_and_slope_kernel(out.m, z=z) - ln_ratio, _, _, _ = kern.nfw_convert(np.exp(ln_c), slope, _ln_q(SOMean(), SOCritical(), om)) + ln_c, slope = Duffy08().ln_concentration_and_slope_kernel(m=out.m, z=z) + ln_ratio, _, _, _ = _nfw_convert(np.exp(ln_c), slope, _ln_q(SOMean(), SOCritical(), om)) np.testing.assert_allclose(out.m * np.exp(ln_ratio), M, rtol=1e-13) @@ -250,15 +271,15 @@ def test_concentration_and_mass_ordering(nfw): M_500c < M_200c < M_vir < M_200m. """ om = OMEGA_M0 - ln_c, slope = Duffy08().ln_concentration_and_slope_kernel(M, z=0.0) + ln_c, slope = Duffy08().ln_concentration_and_slope_kernel(m=M, z=0.0) order = [SOCritical(overdensity=500), SOCritical(), SOVirial(), SOMean()] concentrations, masses = [], [] for target in order: - ln_ratio, c_new, _, _ = kern.nfw_convert(np.exp(ln_c), slope, _ln_q(target, SOCritical())) + ln_ratio, c_new, _, _ = _nfw_convert(np.exp(ln_c), slope, _ln_q(target, SOCritical())) concentrations.append(c_new) masses.append(M * np.exp(ln_ratio)) np.testing.assert_allclose( - nfw.convert_kernel(M, z=0.0, omega_m_z=om, source=SOCritical(), target=target).m, + nfw.convert_kernel(m=M, z=0.0, omega_m_z=om, source=SOCritical(), target=target).m, masses[-1], rtol=1e-14, ) @@ -285,9 +306,9 @@ def _ln_concentration_and_slope(self, m, z): def test_a_user_concentration_model(): """With a constant c the mapping is a constant ratio mu(c') / mu(c), and dM'/dM it.""" conv = ProfileConversion(concentration=_ConstantConcentration()) - out = conv.convert_kernel(M, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical()) - c_new = kern.nfw_concentration_ratio_solve(5.0, _ln_q(SOCritical(), SOMean(), 0.3)) - ratio = math.exp(kern.nfw_ln_mu(c_new) - kern.nfw_ln_mu(5.0)) + out = conv.convert_kernel(m=M, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical()) + c_new = _solve(5.0, _ln_q(SOCritical(), SOMean(), 0.3)) + ratio = math.exp(_ln_mu(c_new) - _ln_mu(5.0)) np.testing.assert_allclose(out.m / M, ratio, rtol=1e-14) np.testing.assert_allclose(out.dm_dm, ratio, rtol=1e-14) # No calibration domain: everything is inside, and no policy applies to it. @@ -305,14 +326,14 @@ def test_a_user_concentration_model(): def test_fof_converts_through_its_overdensity(nfw): """FoF(b=0.2) is ~179 times the mean, so a little bigger than 200m: a few per cent.""" - out = nfw.convert_kernel(M, z=0.0, omega_m_z=OMEGA_M0, source=FOF(), target=SOMean()) + out = nfw.convert_kernel(m=M, z=0.0, omega_m_z=OMEGA_M0, source=FOF(), target=SOMean()) assert np.all((out.m / M < 1.0) & (out.m / M > 0.95)) def test_definitions_without_an_overdensity_raise(nfw): for bad in (SelfBound(), SOGeneric(preferred=SOMean())): with pytest.raises(ValueError, match="Can't convert masses"): - nfw.convert_kernel(M, z=0.0, omega_m_z=0.3, source=bad, target=SOMean()) + nfw.convert_kernel(m=M, z=0.0, omega_m_z=0.3, source=bad, target=SOMean()) with pytest.raises(ValueError, match="Can't convert masses"): nfw.check_definitions(SOMean(), bad) @@ -320,30 +341,36 @@ def test_definitions_without_an_overdensity_raise(nfw): def test_inputs_are_checked(nfw): args = {"source": SOMean(), "target": SOCritical()} with pytest.raises(DomainError, match="m in"): - nfw.convert_kernel([1e12, 0.0], z=0.0, omega_m_z=0.3, **args) + nfw.convert_kernel(m=[1e12, 0.0], z=0.0, omega_m_z=0.3, **args) with pytest.raises(DomainError, match="z in"): - nfw.convert_kernel(1e12, z=-1.0, omega_m_z=0.3, **args) + nfw.convert_kernel(m=1e12, z=-1.0, omega_m_z=0.3, **args) with pytest.raises(DomainError, match="omega_m_z"): - nfw.convert_kernel(1e12, z=0.0, omega_m_z=np.nan, **args) + nfw.convert_kernel(m=1e12, z=0.0, omega_m_z=np.nan, **args) def test_calibration_domain_is_reported(nfw): out = nfw.convert_kernel( - np.array([1e9, 1e12, 1e16]), z=0.0, omega_m_z=0.3, source=SOCritical(), target=SOMean() + m=np.array([1e9, 1e12, 1e16]), z=0.0, omega_m_z=0.3, source=SOCritical(), target=SOMean() ) np.testing.assert_array_equal(out.in_calibration_domain, [False, True, False]) assert "Duffy08's calibration domain" in out.calibration_description - high_z = nfw.convert_kernel(1e12, z=3.0, omega_m_z=0.9, source=SOMean(), target=SOVirial()) + high_z = nfw.convert_kernel(m=1e12, z=3.0, omega_m_z=0.9, source=SOMean(), target=SOVirial()) assert not high_z.in_calibration_domain def test_public_convert(nfw): - m, jac = nfw.convert(1e12 * Msun_h, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical()) + m, jac = nfw.convert( + m=1e12 * Msun_h, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical() + ) assert m.unit == Msun_h assert m.shape == () assert isinstance(jac, float) with pytest.raises(UnitBoundaryError): - nfw.convert(1e12, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical()) + nfw.convert(m=1e12, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical()) + with pytest.raises(TypeError, match="positional"): + nfw.convert(1e12 * Msun_h, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical()) + with pytest.raises(TypeError, match="positional"): + nfw.convert_kernel(1e12, z=0.0, omega_m_z=0.3, source=SOMean(), target=SOCritical()) # --------------------------------------------------------------------------------- @@ -421,14 +448,14 @@ def test_halo_number_is_conserved(smt, z, target): m_top = 1e16 om = omega_m(Planck18, z, "cb") m_meas_top = mf.mass_conversion.convert_kernel( - m_top, z=z, omega_m_z=om, source=target, target=SOVirial() + m=m_top, z=z, omega_m_z=om, source=target, target=SOVirial() ).m for m0 in m_lo: grid = np.geomspace(m0, m_top, 801) dndlnm = grid * mf.dndm(m=grid * Msun_h, z=z).value counted = simpson(dndlnm, x=np.log(grid)) m_meas = mf.mass_conversion.convert_kernel( - m0, z=z, omega_m_z=om, source=target, target=SOVirial() + m=m0, z=z, omega_m_z=om, source=target, target=SOVirial() ).m expected = ( smt.ngtm(m=m_meas * Msun_h, z=z).value - smt.ngtm(m=m_meas_top * Msun_h, z=z).value @@ -453,7 +480,7 @@ def test_rho_gtm_is_the_mass_in_the_new_definition(smt200): def test_converting_to_an_equal_overdensity_is_the_identity(smt): """SO-mean of the virial overdensity at z is the virial definition there (#362).""" z = 0.5 - delta = float(SOVirial().delta_halo_mean_kernel(omega_m(Planck18, z, "cb"))) + delta = float(SOVirial().delta_halo_mean_kernel(omega_m_z=omega_m(Planck18, z, "cb"))) mf = smt.evolve(mass_definition=SOMean(overdensity=delta), mass_conversion="profile") assert mf.converts_masses m = M * Msun_h @@ -508,7 +535,7 @@ def test_behroozi_uses_the_measured_ngtm(smt): z = 3.0 om = omega_m(Planck18, z, "cb") out = converted.mass_conversion.convert_kernel( - M, z=z, omega_m_z=om, source=SOMean(), target=SOVirial() + m=M, z=z, omega_m_z=om, source=SOMean(), target=SOVirial() ) np.testing.assert_allclose( converted.dndm(m=M * Msun_h, z=z).value, @@ -530,7 +557,7 @@ def test_bocquet_mass_ratio_is_applied_once(smt): z = 0.5 m = np.geomspace(1e13, 1e15, 9) out = converted.mass_conversion.convert_kernel( - m, z=z, omega_m_z=omega_m(Planck18, z, "cb"), source=SOMean(), target=SOCritical() + m=m, z=z, omega_m_z=omega_m(Planck18, z, "cb"), source=SOMean(), target=SOCritical() ) np.testing.assert_allclose( converted.dndm(m=m * Msun_h, z=z).value, @@ -655,7 +682,7 @@ def _rho_s(self, c): original = v3_mdefs.BaseMassDefinition.change_definition def profile_free(self, m, mdef, profile=None, c=None, z=0, cosmo=None): - c = duffy.concentration(np.asarray(m) * Msun_h, z=z) + c = duffy.concentration(m=np.asarray(m) * Msun_h, z=z) return original(self, m, mdef, profile=_NFW(self, z, cosmo), c=c, z=z, cosmo=cosmo) monkeypatch.setattr(v3_mdefs.BaseMassDefinition, "change_definition", profile_free) diff --git a/tests/core/test_core_mass_definitions.py b/tests/core/test_core_mass_definitions.py index 49266ae5..6bd2a9f5 100644 --- a/tests/core/test_core_mass_definitions.py +++ b/tests/core/test_core_mass_definitions.py @@ -89,9 +89,22 @@ def test_resolved(): assert d.resolved is d -@pytest.mark.parametrize("definition", [d for d in ALL if not isinstance(d, SOGeneric)], ids=str) +@pytest.mark.parametrize( + "definition", + [ + *ALL, + SOGeneric(preferred=SOCritical(overdensity=500)), + # Numbers that ``:g`` would round are written in full. + SOMean(overdensity=18 * math.pi**2), + SOGeneric(preferred=SOCritical(overdensity=18 * math.pi**2)), + FOF(linking_length=0.123456789), + ], + ids=str, +) def test_names_round_trip(definition): + """Every definition type reads its own name back, to an equal value.""" assert MassDefinition.from_name(definition.name) == definition + assert MassDefinition.from_name(definition.name.upper()) == definition assert coerce(definition.name) == definition assert coerce(definition) is definition @@ -103,8 +116,18 @@ def test_names(): MassDefinition() assert str(SOCritical(overdensity=500)) == "500c" assert str(SOGeneric(preferred=SOMean(overdensity=178))) == "so(preferred 178m)" + assert str(SOMean(overdensity=18 * math.pi**2)) == f"{18 * math.pi**2!r}m" with pytest.raises(ValueError, match="Unknown mass definition name"): MassDefinition.from_name("200x") + for bad in ( + "so(preferred fof)", + "so(preferred so(preferred 200m))", + "so(preferred self_bound)", + ): + with pytest.raises(ValueError, match="must be an SOMean, SOCritical or SOVirial"): + MassDefinition.from_name(bad) + with pytest.raises(ValueError, match="Unknown mass definition name"): + MassDefinition.from_name("so(preferred)") with pytest.raises(TypeError, match="Expected a MassDefinition"): coerce(200) @@ -136,8 +159,8 @@ def test_pickle_and_copy(definition): def test_so_mean_overdensity_is_unchanged(): mdef = SOMean(overdensity=200) - np.testing.assert_array_equal(mdef.delta_halo_mean_kernel(_OMEGA_M_Z), 200.0) - assert mdef.delta_halo_mean_kernel(0.3).shape == () + np.testing.assert_array_equal(mdef.delta_halo_mean_kernel(omega_m_z=_OMEGA_M_Z), 200.0) + assert mdef.delta_halo_mean_kernel(omega_m_z=0.3).shape == () def test_so_critical_overdensity_is_the_same_halo_density(): @@ -149,64 +172,66 @@ def test_so_critical_overdensity_is_the_same_halo_density(): """ mdef = SOCritical(overdensity=500) z = np.array([0.0, 0.5, 2.0, 6.0]) - delta_m = mdef.delta_halo_mean_kernel(omega_m(Planck18, z)) + delta_m = mdef.delta_halo_mean_kernel(omega_m_z=omega_m(Planck18, z)) rho_m = Planck18.Om0 * Planck18.critical_density0 * (1 + z) ** 3 np.testing.assert_allclose( (delta_m * rho_m).to_value(u.g / u.cm**3), (500 * Planck18.critical_density(z)).to_value(u.g / u.cm**3), rtol=1e-13, ) - assert mdef.delta_halo_mean_kernel(1.0) == 500.0 - assert np.all(mdef.delta_halo_mean_kernel(_OMEGA_M_Z[:-1]) > 500.0) + assert mdef.delta_halo_mean_kernel(omega_m_z=1.0) == 500.0 + assert np.all(mdef.delta_halo_mean_kernel(omega_m_z=_OMEGA_M_Z[:-1]) > 500.0) def test_bryan_norman_virial_overdensity(): """18 pi^2 for Omega_m = 1; ~330 times the mean density for Planck18 today.""" mdef = SOVirial() - assert mdef.delta_halo_mean_kernel(1.0) == 18 * np.pi**2 + assert mdef.delta_halo_mean_kernel(omega_m_z=1.0) == 18 * np.pi**2 # Bryan & Norman (1998): Delta_c ~ 100 for Omega_m ~ 0.3, so Delta_m ~ 330. - delta_m0 = mdef.delta_halo_mean_kernel(omega_m(Planck18, 0.0)) + delta_m0 = mdef.delta_halo_mean_kernel(omega_m_z=omega_m(Planck18, 0.0)) assert 300 < delta_m0 < 360 # It decreases towards 18 pi^2 as Omega_m(z) -> 1 at high z. - delta = mdef.delta_halo_mean_kernel(omega_m(Planck18, np.array([0.0, 1.0, 3.0, 10.0]))) + delta = mdef.delta_halo_mean_kernel( + omega_m_z=omega_m(Planck18, np.array([0.0, 1.0, 3.0, 10.0])) + ) assert np.all(np.diff(delta) < 0) assert delta[-1] == pytest.approx(18 * np.pi**2, rel=1e-2) def test_virial_is_so_mean_18pi2_in_einstein_de_sitter(): """With Omega_m(z) = 1 the two are the same overdensity, to the last bit.""" - assert SOVirial().delta_halo_mean_kernel(1.0) == SOMean( + assert SOVirial().delta_halo_mean_kernel(omega_m_z=1.0) == SOMean( overdensity=18 * math.pi**2 - ).delta_halo_mean_kernel(1.0) + ).delta_halo_mean_kernel(omega_m_z=1.0) def test_fof_overdensity(): """9 / (2 pi b^3): about 179 for b = 0.2, falling as b^-3, whatever Omega_m(z).""" - delta = FOF(linking_length=0.2).delta_halo_mean_kernel(_OMEGA_M_Z) + delta = FOF(linking_length=0.2).delta_halo_mean_kernel(omega_m_z=_OMEGA_M_Z) np.testing.assert_allclose(delta, 9 / (2 * np.pi * 0.2**3), rtol=1e-15) assert 178 < delta[0] < 180 wide = FOF(linking_length=0.4) - assert wide.delta_halo_mean_kernel(0.3) == pytest.approx(delta[0] / 8, rel=1e-14) + assert wide.delta_halo_mean_kernel(omega_m_z=0.3) == pytest.approx(delta[0] / 8, rel=1e-14) def test_so_generic_uses_its_preferred_definition_and_self_bound_has_none(): watson = fits.Watson.measured_mass_definition assert isinstance(watson, SOGeneric) np.testing.assert_array_equal( - watson.delta_halo_mean_kernel(_OMEGA_M_Z), - watson.preferred.delta_halo_mean_kernel(_OMEGA_M_Z), + watson.delta_halo_mean_kernel(omega_m_z=_OMEGA_M_Z), + watson.preferred.delta_halo_mean_kernel(omega_m_z=_OMEGA_M_Z), ) with pytest.raises(ValueError, match="no overdensity"): - SelfBound().delta_halo_mean_kernel(0.3) + SelfBound().delta_halo_mean_kernel(omega_m_z=0.3) with pytest.raises(ValueError, match="no overdensity"): - SelfBound().m_to_r_kernel(1e12, omega_m_z=0.3, rho_mean0=1e11) + SelfBound().radius_from_m_kernel(m=1e12, omega_m_z=0.3, rho_mean0=1e11) def test_delta_halo_mean_checks_omega_m(): with pytest.raises(DomainError, match="omega_m_z"): - SOVirial().delta_halo_mean_kernel([0.3, 0.0]) + SOVirial().delta_halo_mean_kernel(omega_m_z=[0.3, 0.0]) with pytest.raises(DomainError, match="omega_m_z"): - SOGeneric(preferred=SOMean()).delta_halo_mean_kernel(np.nan) + SOGeneric(preferred=SOMean()).delta_halo_mean_kernel(omega_m_z=np.nan) # --------------------------------------------------------------------------------- @@ -228,13 +253,13 @@ def test_radius_of_unit_mean_overdensity_is_the_lagrangian_radius(): """SO-mean(1) encloses the mean density: r = (3m / (4 pi rho_cb))^(1/3), the top hat's.""" m = np.logspace(8, 16, 9) rho = rho_mean0(Planck18, "cb") - r = SOMean(overdensity=1).m_to_r_kernel(m, omega_m_z=0.3, rho_mean0=rho) + r = SOMean(overdensity=1).radius_from_m_kernel(m=m, omega_m_z=0.3, rho_mean0=rho) np.testing.assert_allclose(4 * np.pi * r**3 * rho / 3, m, rtol=1e-14) # Haloes are denser, so smaller: r scales as overdensity^(-1/3). - r200 = SOMean().m_to_r_kernel(m, omega_m_z=0.3, rho_mean0=rho) + r200 = SOMean().radius_from_m_kernel(m=m, omega_m_z=0.3, rho_mean0=rho) np.testing.assert_allclose(r200 / r, 200 ** (-1 / 3), rtol=1e-14) # 1e12 Msun/h at 200m: about 0.2 Mpc/h (comoving). - assert 0.15 < SOMean().m_to_r_kernel(1e12, omega_m_z=0.3, rho_mean0=rho) < 0.25 + assert 0.15 < SOMean().radius_from_m_kernel(m=1e12, omega_m_z=0.3, rho_mean0=rho) < 0.25 @pytest.mark.parametrize("definition", [SOMean(), SOCritical(), SOVirial(), FOF()], ids=str) @@ -242,19 +267,35 @@ def test_mass_radius_round_trip(definition): m = np.logspace(6, 16, 11) om = omega_m(Planck18, 1.0) rho = rho_mean0(Planck18) - r = definition.m_to_r_kernel(m, omega_m_z=om, rho_mean0=rho) + r = definition.radius_from_m_kernel(m=m, omega_m_z=om, rho_mean0=rho) np.testing.assert_allclose( - definition.r_to_m_kernel(r, omega_m_z=om, rho_mean0=rho), m, rtol=1e-14 + definition.m_from_radius_kernel(r=r, omega_m_z=om, rho_mean0=rho), m, rtol=1e-14 ) def test_public_methods_take_and_give_quantities(): rho = rho_mean0(Planck18) * rho_unit - r = SOMean().m_to_r(1e12 * Msun_h, omega_m_z=0.3, rho_mean0=rho) + r = SOMean().radius_from_m(m=1e12 * Msun_h, omega_m_z=0.3, rho_mean0=rho) assert r.unit == Mpc_h assert r.shape == () - m = SOMean().r_to_m(r, omega_m_z=0.3, rho_mean0=rho) + m = SOMean().m_from_radius(r=r, omega_m_z=0.3, rho_mean0=rho) assert m.unit == Msun_h assert m.value == pytest.approx(1e12, rel=1e-14) assert SOMean().halo_density(omega_m_z=0.3, rho_mean0=rho).unit == rho_unit - assert isinstance(SOVirial().delta_halo_mean(1.0), float) + assert isinstance(SOVirial().delta_halo_mean(omega_m_z=1.0), float) + + +def test_arguments_are_keyword_only(): + rho = rho_mean0(Planck18) * rho_unit + with pytest.raises(TypeError, match="positional"): + SOMean().radius_from_m(1e12 * Msun_h, omega_m_z=0.3, rho_mean0=rho) + with pytest.raises(TypeError, match="positional"): + SOMean().m_from_radius(1.0 * Mpc_h, omega_m_z=0.3, rho_mean0=rho) + with pytest.raises(TypeError, match="positional"): + SOMean().radius_from_m_kernel(1e12, omega_m_z=0.3, rho_mean0=1e11) + with pytest.raises(TypeError, match="positional"): + SOMean().m_from_radius_kernel(1.0, omega_m_z=0.3, rho_mean0=1e11) + with pytest.raises(TypeError, match="positional"): + SOMean().delta_halo_mean_kernel(0.3) + with pytest.raises(TypeError, match="positional"): + SOMean().delta_halo_mean(0.3) diff --git a/tests/core/test_core_profiles.py b/tests/core/test_core_profiles.py index cd8dfd2f..fe7e8e34 100644 --- a/tests/core/test_core_profiles.py +++ b/tests/core/test_core_profiles.py @@ -8,7 +8,6 @@ import pytest from hmf.core import _references as refs -from hmf.core._kernels import mass_conversion as kern from hmf.core.domain import DomainError from hmf.core.mass_conversion import ProfileConversion from hmf.core.mass_definitions import SOCritical, SOMean, SOVirial @@ -47,26 +46,43 @@ def test_registered_and_cited(): def test_nfw_enclosed_mass(): """mu(1) = ln 2 - 1/2, and mu -> x^2 / 2 at small x (the density ~ 1/x).""" - assert NFW().mu(1.0) == pytest.approx(math.log(2) - 0.5, rel=1e-15) - assert isinstance(NFW().mu(1.0), float) - assert NFW().mu(1e-5) == pytest.approx(0.5e-10, rel=1e-4) - ln_mu, psi = NFW().ln_mu_and_slope_kernel(np.array([0.5, 2.0])) - np.testing.assert_array_equal(ln_mu, kern.nfw_ln_mu([0.5, 2.0])) - np.testing.assert_array_equal(psi, kern.nfw_ln_mu_slope([0.5, 2.0])) + assert NFW().mu(x=1.0) == pytest.approx(math.log(2) - 0.5, rel=1e-15) + assert isinstance(NFW().mu(x=1.0), float) + assert NFW().mu(x=1e-5) == pytest.approx(0.5e-10, rel=1e-4) + # The closed form, mu = ln(1 + x) - x / (1 + x), and its slope + # psi = x mu'(x) / mu = x^2 / [(1 + x)^2 mu]. + x = np.array([0.5, 2.0]) + mu = np.log1p(x) - x / (1 + x) + ln_mu, psi = NFW().ln_mu_and_slope_kernel(x=x) + np.testing.assert_allclose(ln_mu, np.log(mu), rtol=1e-14) + np.testing.assert_allclose(psi, x**2 / ((1 + x) ** 2 * mu), rtol=1e-14) with pytest.raises(DomainError, match="x in"): - NFW().mu(-1.0) + NFW().mu(x=-1.0) -def test_nfw_conversion_is_the_kernel_with_hu_kravtsov(): - out = NFW().convert_concentration_kernel(C, slope=-0.1, ln_q=LN_Q) - for got, expected in zip(out, kern.nfw_convert(C, -0.1, LN_Q), strict=True): - np.testing.assert_array_equal(got, expected) +def test_nfw_solve_reaches_round_off_from_hu_kravtsov(): + """From Hu & Kravtsov (2003, eq. C11), NFW needs fewer Newton steps than the default.""" + assert NFW.newton_steps < DensityProfile.newton_steps + _, c_new, _, _ = NFW().convert_concentration_kernel(c=C, slope=0.0, ln_q=LN_Q) + ln_mu_new, _ = NFW().ln_mu_and_slope_kernel(x=c_new) + ln_mu_old, _ = NFW().ln_mu_and_slope_kernel(x=np.broadcast_to(C, c_new.shape)) + residual = (ln_mu_new - 3 * np.log(c_new)) - (ln_mu_old - 3 * np.log(C)) - LN_Q + assert np.abs(residual).max() < 2e-13 + + +def test_kernels_take_keyword_arguments_only(): + with pytest.raises(TypeError, match="positional"): + NFW().mu(1.0) + with pytest.raises(TypeError, match="positional"): + NFW().ln_mu_and_slope_kernel(np.array([1.0])) + with pytest.raises(TypeError, match="positional"): + NFW().convert_concentration_kernel(5.0, slope=0.0, ln_q=0.5) def test_isothermal_sphere_is_analytic(): """F = x^-2, so c' = c q^(-1/2) and M'/M = q^(-1/2), for any c (dM'/dM = that).""" ln_ratio, c_new, dln_m, dln_c = _Isothermal().convert_concentration_kernel( - C, slope=-0.1, ln_q=LN_Q + c=C, slope=-0.1, ln_q=LN_Q ) np.testing.assert_allclose(c_new, C * np.exp(-LN_Q / 2), rtol=1e-14) np.testing.assert_allclose(ln_ratio, np.broadcast_to(-LN_Q / 2, ln_ratio.shape), atol=1e-13) @@ -77,7 +93,7 @@ def test_isothermal_sphere_is_analytic(): def test_hernquist_solve_reaches_round_off_from_the_default_start(): ln_mu_and_slope = _Hernquist()._ln_mu_and_slope - _, c_new, _, _ = _Hernquist().convert_concentration_kernel(C, slope=0.0, ln_q=LN_Q) + _, c_new, _, _ = _Hernquist().convert_concentration_kernel(c=C, slope=0.0, ln_q=LN_Q) ln_mu_new, _ = ln_mu_and_slope(c_new) ln_mu_old, _ = ln_mu_and_slope(np.broadcast_to(C, c_new.shape)) residual = (ln_mu_new - 3 * np.log(c_new)) - (ln_mu_old - 3 * np.log(C)) - LN_Q @@ -88,15 +104,15 @@ def test_hernquist_conversion_round_trips_and_its_jacobian_is_the_derivative(): conv = ProfileConversion(profile=_Hernquist()) m = np.geomspace(1e10, 1e15, 11) om = 0.3 - there = conv.convert_kernel(m, z=0.5, omega_m_z=om, source=SOMean(), target=SOVirial()) - back = conv.convert_kernel(there.m, z=0.5, omega_m_z=om, source=SOVirial(), target=SOMean()) + there = conv.convert_kernel(m=m, z=0.5, omega_m_z=om, source=SOMean(), target=SOVirial()) + back = conv.convert_kernel(m=there.m, z=0.5, omega_m_z=om, source=SOVirial(), target=SOMean()) np.testing.assert_allclose(back.m, m, rtol=1e-12) h = 1e-5 hi = conv.convert_kernel( - m * math.exp(h), z=0.5, omega_m_z=om, source=SOMean(), target=SOVirial() + m=m * math.exp(h), z=0.5, omega_m_z=om, source=SOMean(), target=SOVirial() ) lo = conv.convert_kernel( - m * math.exp(-h), z=0.5, omega_m_z=om, source=SOMean(), target=SOVirial() + m=m * math.exp(-h), z=0.5, omega_m_z=om, source=SOMean(), target=SOVirial() ) numeric = (np.log(hi.m) - np.log(lo.m)) / (2 * h) * there.m / m np.testing.assert_allclose(there.dm_dm, numeric, rtol=1e-8) @@ -106,13 +122,13 @@ def test_a_more_concentrated_profile_shifts_masses_less(): """Hernquist falls as x^-4 outside, NFW as x^-3: less mass between the two radii.""" m = np.geomspace(1e11, 1e15, 5) args = {"z": 0.0, "omega_m_z": 0.3, "source": SOCritical(), "target": SOMean()} - nfw = ProfileConversion().convert_kernel(m, **args).m / m - hernquist = ProfileConversion(profile=_Hernquist()).convert_kernel(m, **args).m / m + nfw = ProfileConversion().convert_kernel(m=m, **args).m / m + hernquist = ProfileConversion(profile=_Hernquist()).convert_kernel(m=m, **args).m / m assert np.all((hernquist > 1) & (hernquist < nfw)) def test_inputs_are_checked(): with pytest.raises(DomainError, match="c in"): - NFW().convert_concentration_kernel([5.0, 0.0], slope=0.0, ln_q=0.5) + NFW().convert_concentration_kernel(c=[5.0, 0.0], slope=0.0, ln_q=0.5) with pytest.raises(DomainError, match="slope must be finite"): - NFW().convert_concentration_kernel(5.0, slope=np.nan, ln_q=0.5) + NFW().convert_concentration_kernel(c=5.0, slope=np.nan, ln_q=0.5) diff --git a/tests/regression/v4_providers.py b/tests/regression/v4_providers.py index a153be45..dd00d173 100644 --- a/tests/regression/v4_providers.py +++ b/tests/regression/v4_providers.py @@ -163,7 +163,7 @@ def v4_fsigma( sigma, z=z, omega_m_z=om, - delta_halo=mdef.delta_halo_mean_kernel(om) if needs_delta else None, + delta_halo=mdef.delta_halo_mean_kernel(omega_m_z=om) if needs_delta else None, delta_c=delta_c, n_eff=n_eff, m=m * Msun_h,