diff --git a/development/find_radiation_threshold.py b/development/find_radiation_threshold.py
new file mode 100644
index 0000000..3abb517
--- /dev/null
+++ b/development/find_radiation_threshold.py
@@ -0,0 +1,263 @@
+"""Determine a safe radiation threshold for the growth-factor selector.
+
+This utility compares `Tinker08` halo mass functions computed with
+`Eisenstein97GrowthFactor` and `ODEGrowthFactor` over a redshift and mass grid, then
+reports the largest radiation-density threshold for which `dndlnm` stays within a
+chosen relative tolerance for all masses below a configurable upper limit.
+"""
+
+from __future__ import annotations
+
+import argparse
+import sys
+from dataclasses import dataclass
+
+import numpy as np
+from astropy.cosmology import FlatLambdaCDM
+
+from hmf import MassFunction
+
+
+@dataclass(frozen=True)
+class ThresholdConfig:
+ """Configuration for the radiation-threshold scan.
+
+ Parameters
+ ----------
+ z_min
+ Minimum redshift included in the scan.
+ z_max
+ Maximum redshift included in the scan.
+ dz
+ Redshift step size.
+ mass_limit
+ Maximum halo mass, in `Msun / h`, included in the error check.
+ mmin
+ Minimum log10 halo mass passed to `MassFunction`.
+ mmax
+ Maximum log10 halo mass passed to `MassFunction`.
+ dlog10m
+ Log10 halo-mass spacing passed to `MassFunction`.
+ tolerance
+ Maximum allowed relative difference in `dndlnm`.
+ hmf_model
+ Halo mass function model to evaluate.
+ H0
+ Hubble constant for the comparison cosmology.
+ Om0
+ Matter density parameter at `z=0`.
+ Ob0
+ Baryon density parameter at `z=0`.
+ Tcmb0
+ CMB temperature for the comparison cosmology.
+ sigma_8
+ Present-day amplitude of matter fluctuations.
+ n
+ Primordial scalar spectral index.
+ """
+
+ z_min: float = 0.0
+ z_max: float = 12.0
+ dz: float = 0.1
+ mass_limit: float = 1.0e14
+ mmin: float = 10.0
+ mmax: float = 14.2
+ dlog10m: float = 0.02
+ tolerance: float = 0.01
+ hmf_model: str = "Tinker08"
+ H0: float = 67.74
+ Om0: float = 0.3089
+ Ob0: float = 0.0486
+ Tcmb0: float = 2.7255
+ sigma_8: float = 0.8159
+ n: float = 0.9667
+
+ @property
+ def cosmo(self) -> FlatLambdaCDM:
+ """Comparison cosmology used for the scan."""
+ return FlatLambdaCDM(
+ H0=self.H0,
+ Om0=self.Om0,
+ Ob0=self.Ob0,
+ Tcmb0=self.Tcmb0,
+ )
+
+ @property
+ def redshifts(self) -> np.ndarray:
+ """Redshift grid used for the scan."""
+ return np.arange(self.z_min, self.z_max + self.dz / 2, self.dz)
+
+
+def radiation_density(cosmo: FlatLambdaCDM, z: float) -> float:
+ """Compute the fractional radiation density at redshift `z`.
+
+ Parameters
+ ----------
+ cosmo
+ Cosmology used for the comparison.
+ z
+ Redshift at which to evaluate the radiation fraction.
+
+ Returns
+ -------
+ float
+ Radiation density fraction, matching the selector logic used by
+ `GrowthFactor.radiation_density`.
+ """
+ return float((cosmo.Ogamma0 + cosmo.Onu0) * (1 + z) ** 4 * cosmo.inv_efunc(z) ** 2)
+
+
+def build_mass_function(z: float, growth_model: str, config: ThresholdConfig) -> MassFunction:
+ """Construct a comparison `MassFunction` instance.
+
+ Parameters
+ ----------
+ z
+ Redshift of the calculation.
+ growth_model
+ Growth-factor model name to pass to `MassFunction`.
+ config
+ Configuration for the threshold scan.
+
+ Returns
+ -------
+ MassFunction
+ Configured mass function instance.
+ """
+ return MassFunction(
+ z=z,
+ hmf_model=config.hmf_model,
+ growth_model=growth_model,
+ transfer_model="EH",
+ mdef_model="SOMean",
+ mdef_params={"overdensity": 200},
+ Mmin=config.mmin,
+ Mmax=config.mmax,
+ dlog10m=config.dlog10m,
+ cosmo_params={
+ "H0": config.H0,
+ "Om0": config.Om0,
+ "Ob0": config.Ob0,
+ "Tcmb0": config.Tcmb0,
+ },
+ sigma_8=config.sigma_8,
+ n=config.n,
+ )
+
+
+def max_relative_dndlnm_error(z: float, config: ThresholdConfig) -> float:
+ """Return the maximum relative `dndlnm` error at one redshift.
+
+ Parameters
+ ----------
+ z
+ Redshift at which to compare the growth models.
+ config
+ Configuration for the threshold scan.
+
+ Returns
+ -------
+ float
+ Maximum relative difference between the Eisenstein and ODE `dndlnm`
+ predictions for masses below `config.mass_limit`.
+ """
+ eisenstein = build_mass_function(z, "Eisenstein97GrowthFactor", config)
+ ode = build_mass_function(z, "ODEGrowthFactor", config)
+
+ mask = eisenstein.m < config.mass_limit
+ rel = np.abs(eisenstein.dndlnm[mask] - ode.dndlnm[mask]) / ode.dndlnm[mask]
+ return float(np.max(rel))
+
+
+def find_safe_threshold(config: ThresholdConfig) -> list[tuple[float, float, float]]:
+ """Scan the grid and collect redshift, radiation density, and HMF error.
+
+ Parameters
+ ----------
+ config
+ Configuration for the threshold scan.
+
+ Returns
+ -------
+ list of tuple
+ Tuples of `(z, radiation_density, max_relative_error)`.
+ """
+ return [
+ (
+ float(z),
+ radiation_density(config.cosmo, float(z)),
+ max_relative_dndlnm_error(float(z), config),
+ )
+ for z in config.redshifts
+ ]
+
+
+def write_line(text: str = "") -> None:
+ """Write one line to standard output."""
+ sys.stdout.write(f"{text}\n")
+
+
+def build_parser() -> argparse.ArgumentParser:
+ """Create the command-line parser for the utility."""
+ parser = argparse.ArgumentParser(description=__doc__)
+ parser.add_argument("--z-min", type=float, default=ThresholdConfig.z_min)
+ parser.add_argument("--z-max", type=float, default=ThresholdConfig.z_max)
+ parser.add_argument("--dz", type=float, default=ThresholdConfig.dz)
+ parser.add_argument("--mass-limit", type=float, default=ThresholdConfig.mass_limit)
+ parser.add_argument("--mmin", type=float, default=ThresholdConfig.mmin)
+ parser.add_argument("--mmax", type=float, default=ThresholdConfig.mmax)
+ parser.add_argument("--dlog10m", type=float, default=ThresholdConfig.dlog10m)
+ parser.add_argument("--tolerance", type=float, default=ThresholdConfig.tolerance)
+ parser.add_argument("--hmf-model", type=str, default=ThresholdConfig.hmf_model)
+ return parser
+
+
+def main() -> int:
+ """Run the threshold scan and print a compact report."""
+ args = build_parser().parse_args()
+ config = ThresholdConfig(
+ z_min=args.z_min,
+ z_max=args.z_max,
+ dz=args.dz,
+ mass_limit=args.mass_limit,
+ mmin=args.mmin,
+ mmax=args.mmax,
+ dlog10m=args.dlog10m,
+ tolerance=args.tolerance,
+ hmf_model=args.hmf_model,
+ )
+
+ rows = find_safe_threshold(config)
+ safe_rows = [row for row in rows if row[2] <= config.tolerance]
+ failing_rows = [row for row in rows if row[2] > config.tolerance]
+
+ write_line("z radiation_density max_rel_dndlnm_error")
+ for z, rad, err in rows:
+ write_line(f"{z:5.2f} {rad:17.10f} {err:20.10f}")
+
+ if not safe_rows:
+ write_line("\nNo safe redshifts found on the requested grid.")
+ return 1
+
+ z_safe, threshold, err_safe = safe_rows[-1]
+ write_line(
+ "\nLargest safe threshold on scanned grid: "
+ f"{threshold:.10f} at z={z_safe:.2f} with max_rel_error={err_safe:.6f}"
+ )
+
+ if failing_rows:
+ z_fail, rad_fail, err_fail = failing_rows[0]
+ write_line(
+ "First failing point: "
+ f"z={z_fail:.2f}, radiation_density={rad_fail:.10f}, max_rel_error={err_fail:.6f}"
+ )
+ else:
+ write_line(
+ "No failing points found on the scanned grid; increase --z-max to probe further."
+ )
+
+ return 0
+
+
+if __name__ == "__main__":
+ raise SystemExit(main())
diff --git a/docs/technical/colossus_comparison.rst b/docs/technical/colossus_comparison.rst
new file mode 100644
index 0000000..c23b02a
--- /dev/null
+++ b/docs/technical/colossus_comparison.rst
@@ -0,0 +1,124 @@
+Colossus Comparison Notes
+=========================
+
+This note documents the known causes of differences between ``hmf`` and
+`Colossus `_ for the native
+``Tinker08`` halo mass function.
+In versions pre-3.6.1, a major difference was the growth factor computation, which was
+definitively less accurate in ``hmf`` than in Colossus. However, after
+fixing the growth factor implementation in
+`this PR `_ for v3.6.0, and then tightening
+the growth-selector threshold in v3.6.1, some small residual differences remain,
+particularly at high redshift.
+
+The goal of this note is not to argue that either code is definitively "more
+correct". Instead, it records the main implementation choices that explain the
+observed residual differences so that users and developers understand where
+agreement is expected and where small systematic offsets are normal.
+
+Setup of the comparison
+-----------------------
+
+The comparisons discussed here were made with:
+
+- the native ``200m`` form of ``Tinker08``,
+- matched flat cosmologies with ``H0 = 67.74``, ``Om0 = 0.3089``,
+ ``Ob0 = 0.0486``, ``sigma8 = 0.8159``, and ``ns = 0.9667``,
+- the ``EH`` transfer model in ``hmf``, and
+- direct comparisons of ``dndlnm`` at representative masses
+ :math:`10^{11}`, :math:`10^{12}`, and :math:`10^{13}\,M_\odot/h`.
+
+With this setup, the mismatch is small at low redshift and grows
+toward high redshift, particularly in the rare-halo tail.
+
+What is *not* driving the mismatch
+----------------------------------
+
+Several obvious suspects were checked and found not to be the dominant cause:
+
+- **Mass-definition conversion:** this comparison uses native ``200m``
+ ``Tinker08`` predictions, so it does not rely on the mass-definition
+ conversion path.
+- **Transfer-function normalization at** :math:`z=0`: ``hmf`` and Colossus
+ agree on :math:`\sigma(M,0)` at about the :math:`10^{-4}` level for the
+ matched setup.
+- **Slope term:** the logarithmic slope entering ``dndlnm``,
+ :math:`d \ln \sigma / d \ln R`, agrees at the sub-:math:`0.1\%` level.
+
+The main causes of the residual difference
+------------------------------------------
+
+Two effects dominate the remaining mismatch.
+
+High-redshift growth implementation
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+After the selector update, ``hmf`` uses the full ODE growth solution whenever
+the radiation fraction exceeds the calibrated threshold (essentially for z>1.5).
+Colossus uses a different hybrid approach for LCDM cosmologies:
+
+- an analytic matter-radiation approximation at high redshift, and
+- an integral solution at low redshift,
+
+with a transition regime between them.
+
+These are both reasonable algorithmic choices, but they do not produce exactly
+the same high-redshift growth history. In the matched comparison used here,
+``hmf`` ends up with a slightly larger :math:`\sigma(M, z)` than Colossus at
+high redshift, typically by about :math:`0.26`--:math:`0.34\%` over
+``z = 6``--``10`` for the masses tested.
+
+That difference is tiny in :math:`\sigma` itself, but it is evaluated in the
+exponential tail of the halo mass function, where very small shifts in
+:math:`\sigma` can produce multi-percent shifts in abundance.
+
+Tinker08 coefficient precision
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The two codes also do not use numerically identical ``Tinker08`` coefficients.
+
+``hmf`` stores the more precise coefficient values, for example:
+
+- ``A_200 = 0.1858659``
+- ``a_200 = 1.466904``
+- ``b_200 = 2.571104``
+- ``c_200 = 1.193958``
+
+Colossus uses the rounded table values:
+
+- ``A_200 = 0.186``
+- ``a_200 = 1.47``
+- ``b_200 = 2.57``
+- ``c_200 = 1.19``
+
+At fixed :math:`\sigma`, this coefficient rounding changes :math:`f(\sigma)` by
+only a little at low redshift, but by several percent in the high-redshift
+tail. In the matched tests performed here, the coefficient choice alone
+accounts for approximately:
+
+- :math:`2`--:math:`6\%` at ``z = 6``,
+- :math:`3`--:math:`10\%` at ``z = 8``, and
+- :math:`4`--:math:`15\%` at ``z = 10``,
+
+depending on mass.
+
+How the effects combine
+-----------------------
+
+The two dominant effects push in opposite directions:
+
+- the slightly larger high-redshift :math:`\sigma(M, z)` in ``hmf`` tends to
+ increase the abundance relative to Colossus, while
+- the more precise ``Tinker08`` coefficients in ``hmf`` tend to decrease the
+ abundance relative to Colossus' rounded table.
+
+Because these effects partially cancel, the final mismatch is significantly
+smaller than either ingredient by itself. For the matched setup used in the
+external regression test, the resulting ``hmf`` versus Colossus difference is
+typically of order:
+
+- ``z = 6``: :math:`\sim 0.1\%` to :math:`2.7\%`,
+- ``z = 8``: :math:`\sim 1\%` to :math:`7\%`,
+- ``z = 10``: :math:`\sim 2\%` to :math:`14\%`,
+
+over the range :math:`10^{11}`--:math:`10^{13}\,M_\odot/h`.
diff --git a/docs/technical/index.rst b/docs/technical/index.rst
index 6cfd38a..d0cfa36 100644
--- a/docs/technical/index.rst
+++ b/docs/technical/index.rst
@@ -6,3 +6,4 @@ Technical Documentation and Derivations
:maxdepth: 1
growth_factors
+ colossus_comparison
diff --git a/src/hmf/cosmology/growth_factor.py b/src/hmf/cosmology/growth_factor.py
index 2ccd783..00df338 100644
--- a/src/hmf/cosmology/growth_factor.py
+++ b/src/hmf/cosmology/growth_factor.py
@@ -49,6 +49,12 @@
HAVE_CAMB = False
+# The low-radiation analytic/integral branches are only used while they keep the
+# Tinker08 HMF within 1% of the full ODE solution below 1e14 Msun/h, based on the
+# calibration in development/find_radiation_threshold.py.
+LOW_RADIATION_THRESHOLD = 5.0e-4
+
+
@pluggable
class BaseGrowthFactor(Cmpt):
r"""General class for a growth factor calculation.
@@ -233,7 +239,7 @@ def _choose_solution(self, z):
return self._ode_gf
# Low radiation component solutions.
- if np.all(self.radiation_density(z) < 0.02):
+ if np.all(self.radiation_density(z) < LOW_RADIATION_THRESHOLD):
if self.cosmo.is_flat:
return self._eisenstein_gf
if self.cosmo.Ode0 == 0:
@@ -412,7 +418,7 @@ class IntegralGrowthFactor(BaseGrowthFactor):
def _validate_assumptions(self, z: float | np.ndarray):
- if np.any(self.radiation_density(z) > 0.02):
+ if np.any(self.radiation_density(z) > LOW_RADIATION_THRESHOLD):
warnings.warn(
f"The {self.__class__.__name__} is not accurate when the radiation "
f"density is significant. Consider using the "
@@ -495,7 +501,7 @@ def _growth_rate(self, z):
# This is because the coefficients in front of the second term are chosen
# based on the normalization specified in _d_plus_unnormalized (i.e. 5/2 Om0).
# This is a convenient choice so that D(a) -> a for a -> 0.
- if np.any(self.radiation_density(z) > 0.02):
+ if np.any(self.radiation_density(z) > LOW_RADIATION_THRESHOLD):
# Need to use the basic spline because that is always consistent
# with the growth factor. The analytic equation is only consistent with
# the growth factor under the assumption of negligible radiation.
diff --git a/src/hmf/density_field/filters.py b/src/hmf/density_field/filters.py
index 5cc58b4..394d11f 100644
--- a/src/hmf/density_field/filters.py
+++ b/src/hmf/density_field/filters.py
@@ -274,7 +274,7 @@ def sigma(self, r, order=0, rk=None):
sigma = (0.5 / np.pi**2) * intg.simpson(integ, dx=dlnk, axis=-1)
return np.sqrt(sigma)
- def nu(self, r, delta_c=1.686):
+ def nu(self, r, delta_c=1.68647):
r"""
Peak height, :math:`\frac{\delta_c^2}{\sigma^2(r)}`.
diff --git a/src/hmf/mass_function/fitting_functions.py b/src/hmf/mass_function/fitting_functions.py
index e7be975..4b127a5 100644
--- a/src/hmf/mass_function/fitting_functions.py
+++ b/src/hmf/mass_function/fitting_functions.py
@@ -224,7 +224,7 @@ def __init__(
n_eff: None | np.ndarray = None,
mass_definition: None | md.BaseMassDefinition = None,
cosmo: csm.FLRW = csm.Planck15,
- delta_c: float = 1.686,
+ delta_c: float = 1.68647,
**model_parameters,
):
super().__init__(**model_parameters)
diff --git a/src/hmf/mass_function/hmf.py b/src/hmf/mass_function/hmf.py
index 1bc196e..927d736 100644
--- a/src/hmf/mass_function/hmf.py
+++ b/src/hmf/mass_function/hmf.py
@@ -104,7 +104,7 @@ def __init__(
hmf_params: dict[str, Any] | None = None,
mdef_model: None | str | MassDef = None,
mdef_params: dict | None = None,
- delta_c: float = 1.686,
+ delta_c: float = 1.68647,
filter_model: str | BaseFilter = TopHat,
filter_params: dict | None = None,
disable_mass_conversion: bool = True,
diff --git a/tests/test_fits.py b/tests/test_fits.py
index 8f26c86..9384f23 100644
--- a/tests/test_fits.py
+++ b/tests/test_fits.py
@@ -3,6 +3,8 @@
import numpy as np
import pytest
+from colossus.cosmology.cosmology import setCosmology
+from colossus.lss.mass_function import massFunction
from hmf import MassFunction
from hmf.mass_function import fitting_functions as ff
@@ -21,6 +23,15 @@
]
+def _conversion_is_active(hmf: MassFunction) -> bool:
+ """Whether `MassFunction.dndm` will apply mass-definition conversion."""
+ return (
+ hmf.hmf.measured_mass_definition is not None
+ and hmf.hmf.measured_mass_definition != hmf.mdef
+ and not hmf.disable_mass_conversion
+ )
+
+
@pytest.fixture(scope="module")
def hmf():
return MassFunction(
@@ -124,16 +135,16 @@ def test_tinker10_neg_beta():
@pytest.mark.filterwarnings("ignore:.*does not match the mass definition.*:UserWarning")
@pytest.mark.filterwarnings("ignore:.*Halo-Exclusion models.*:UserWarning")
-@pytest.mark.parametrize("fit", ["Behroozi", "SMT", "Tinker08"])
+@pytest.mark.parametrize("fit", ["Behroozi", "SMT"])
@pytest.mark.parametrize("z", [0.0, 1.0, 2.0])
def test_hmf_mass_definition_consistency(z, fit):
"""The dndm must be consistent between equivalent SO mass definitions.
SOCritical(200) and SOMean(200/Omega_m(z)) define the same physical density
- threshold, so any HMF with mass-definition conversion must give the same dndm
- for both. This was broken before because the mass-definition conversion in dndm
- always used z=0 and the default Planck15 cosmology instead of the actual
- redshift and cosmology of the computation.
+ threshold, so any HMF with mass-definition conversion must give the same dndm for
+ both. This was broken before because the mass-definition conversion in dndm always
+ used z=0 and the default Planck15 cosmology instead of the actual redshift and
+ cosmology of the computation.
"""
from astropy.cosmology import FlatLambdaCDM
@@ -167,7 +178,107 @@ def test_hmf_mass_definition_consistency(z, fit):
**common,
)
+ assert _conversion_is_active(h_crit)
+ assert _conversion_is_active(h_mean)
+
# Allow ~2% tolerance: a small residual comes from the floating-point
# imprecision of equiv_overdensity = 200/Om(z) and the NFW c-M approximation.
# The original bug produced 20-50% errors, so this tolerance is well above that.
np.testing.assert_allclose(h_crit.dndm, h_mean.dndm, rtol=2e-2)
+
+
+@pytest.mark.filterwarnings("ignore:.*does not match the mass definition.*:UserWarning")
+@pytest.mark.parametrize("z", [0.0, 1.0, 2.0])
+def test_tinker08_native_so_definition_consistency(z):
+ """Tinker08 should agree between equivalent native SO definitions.
+
+ Tinker08 uses `SOGeneric` as its measured mass definition, so it natively supports
+ spherical-overdensity definitions without entering the mass-definition conversion
+ branch used by the PR #281 regression above.
+ """
+ from astropy.cosmology import FlatLambdaCDM
+
+ cosmo_params = {"Om0": 0.3, "H0": 70.0, "Ob0": 0.05}
+ cosmo = FlatLambdaCDM(**cosmo_params)
+ equiv_overdensity = 200.0 / cosmo.Om(z)
+
+ common = {
+ "hmf_model": "Tinker08",
+ "transfer_model": "EH",
+ "Mmin": 10,
+ "Mmax": 15,
+ "dlog10m": 0.1,
+ "z": z,
+ "disable_mass_conversion": False,
+ "cosmo_params": cosmo_params,
+ }
+
+ h_crit = MassFunction(
+ mdef_model="SOCritical",
+ mdef_params={"overdensity": 200},
+ **common,
+ )
+ h_mean = MassFunction(
+ mdef_model="SOMean",
+ mdef_params={"overdensity": equiv_overdensity},
+ **common,
+ )
+
+ assert not _conversion_is_active(h_crit)
+ assert not _conversion_is_active(h_mean)
+ np.testing.assert_allclose(h_crit.dndm, h_mean.dndm, rtol=2e-2)
+
+
+@pytest.mark.parametrize(
+ ("z", "rtol"),
+ [(0.0, 5e-3), (2.0, 1e-2), (4.0, 5e-3), (6.0, 3e-2), (8.0, 8e-2), (10.0, 1.5e-1)],
+)
+def test_tinker08_matches_colossus(z, rtol):
+ """Tinker08 should remain reasonably close to the Colossus implementation.
+
+ This compares the native `200m` Tinker08 prediction against Colossus at several
+ redshifts using a matched cosmology. The tolerance widens with redshift because the
+ remaining mismatch is still driven mostly by different high-z growth treatments and
+ by the fact that `hmf` uses more precise Tinker08 coefficients while Colossus uses
+ rounded table values. After tightening the selector threshold, the agreement is good
+ enough to use a meaningfully stricter external regression than before. See
+ `docs/technical/colossus_comparison.rst` for a longer explanation.
+ """
+ cosmo_params = {"H0": 67.74, "Om0": 0.3089, "Ob0": 0.0486}
+ setCosmology(
+ "hmf-tinker08-test",
+ {
+ "flat": True,
+ "H0": cosmo_params["H0"],
+ "Om0": cosmo_params["Om0"],
+ "Ob0": cosmo_params["Ob0"],
+ "sigma8": 0.8159,
+ "ns": 0.9667,
+ },
+ )
+
+ hmf = MassFunction(
+ hmf_model="Tinker08",
+ transfer_model="EH",
+ mdef_model="SOMean",
+ mdef_params={"overdensity": 200},
+ Mmin=10,
+ Mmax=14.2,
+ dlog10m=0.02,
+ z=z,
+ cosmo_params=cosmo_params,
+ sigma_8=0.8159,
+ n=0.9667,
+ )
+
+ for mass in (1e11, 1e12, 1e13):
+ hmf_dndlnm = np.exp(np.interp(np.log(mass), np.log(hmf.m), np.log(hmf.dndlnm)))
+ colossus_dndlnm = massFunction(
+ mass,
+ z,
+ q_in="M",
+ q_out="dndlnM",
+ mdef="200m",
+ model="tinker08",
+ )
+ assert hmf_dndlnm == pytest.approx(colossus_dndlnm, rel=rtol)
diff --git a/tests/test_genmf.py b/tests/test_genmf.py
index 993894b..458929c 100644
--- a/tests/test_genmf.py
+++ b/tests/test_genmf.py
@@ -99,6 +99,7 @@ def hmf(self):
z=0.0,
transfer_model="FromFile",
growth_model="GenMFGrowth",
+ delta_c=1.686,
)
@staticmethod
diff --git a/tests/test_growth.py b/tests/test_growth.py
index 1e26871..feec5e1 100644
--- a/tests/test_growth.py
+++ b/tests/test_growth.py
@@ -5,6 +5,7 @@
from astropy import cosmology
from astropy.cosmology import Planck13, w0waCDM
+from hmf import MassFunction
from hmf.cosmology import growth_factor
@@ -240,6 +241,41 @@ def test_growth_rate_uses_ode_at_highz():
gf._growth_rate(1000)
+def test_growth_selector_switches_at_calibrated_radiation_threshold():
+ """The default selector should swap to ODE once the calibrated threshold is exceeded."""
+ cosmo = cosmology.FlatLambdaCDM(H0=67.74, Om0=0.3089, Ob0=0.0486, Tcmb0=2.7255)
+ gf = growth_factor.GrowthFactor(cosmo)
+
+ assert gf.radiation_density(1.0) < growth_factor.LOW_RADIATION_THRESHOLD
+ assert gf._choose_solution(1.0) is gf._eisenstein_gf
+
+ assert gf.radiation_density(2.0) > growth_factor.LOW_RADIATION_THRESHOLD
+ assert gf._choose_solution(2.0) is gf._ode_gf
+
+
+def test_growth_selector_keeps_tinker08_within_one_percent_below_threshold():
+ """Below the selector threshold, default Tinker08 should stay within 1% of ODE."""
+ common = {
+ "hmf_model": "Tinker08",
+ "transfer_model": "EH",
+ "mdef_model": "SOMean",
+ "mdef_params": {"overdensity": 200},
+ "Mmin": 10,
+ "Mmax": 14.2,
+ "dlog10m": 0.02,
+ "z": 1.0,
+ "cosmo_params": {"H0": 67.74, "Om0": 0.3089, "Ob0": 0.0486},
+ "sigma_8": 0.8159,
+ "n": 0.9667,
+ }
+ default = MassFunction(**common)
+ ode = MassFunction(growth_model="ODEGrowthFactor", **common)
+
+ mask = default.m < 1e14
+ relative_error = np.abs(default.dndlnm[mask] - ode.dndlnm[mask]) / ode.dndlnm[mask]
+ assert np.max(relative_error) < 1e-2
+
+
def test_expected_warnings():
"""Test that we get expected warnings for unsupported cosmologies."""
cosmo = Planck13