From bbb2f40b7b8b42fde6de6f89bb002570a5798cf8 Mon Sep 17 00:00:00 2001 From: Steven Murray Date: Fri, 22 May 2026 21:40:46 +0200 Subject: [PATCH 1/2] feat: update growth factor threshold and improve documentation - Introduced a low radiation threshold for growth factor calculations. - Updated growth factor logic to use the new threshold. - Enhanced tests to validate growth factor behavior at the new threshold. - Added a utility script to determine safe radiation thresholds. - Created documentation for the comparison between hmf and Colossus. --- development/find_radiation_threshold.py | 263 +++++++++++++++++++++ docs/technical/colossus_comparison.rst | 124 ++++++++++ docs/technical/index.rst | 1 + src/hmf/cosmology/growth_factor.py | 12 +- src/hmf/density_field/filters.py | 2 +- src/hmf/mass_function/fitting_functions.py | 2 +- src/hmf/mass_function/hmf.py | 2 +- tests/test_fits.py | 121 +++++++++- tests/test_growth.py | 36 +++ 9 files changed, 552 insertions(+), 11 deletions(-) create mode 100644 development/find_radiation_threshold.py create mode 100644 docs/technical/colossus_comparison.rst 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_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 From 4d1c500902ae836d63b9cc837da55768c0902542 Mon Sep 17 00:00:00 2001 From: Steven Murray Date: Fri, 22 May 2026 21:46:21 +0200 Subject: [PATCH 2/2] test: revert to delta_c=1.686 for strict regression tests --- tests/test_genmf.py | 1 + 1 file changed, 1 insertion(+) 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