"""
A module defining many several halo mass function fits.
The :class:`~BaseFittingFunction` defines the interface for all fitting functions, and
the included fitting functions are subclasses of this. See the documentation of
:class:`~BaseFittingFunction` for details on how to define a new fitting function.
"""
import inspect
import re
import warnings
from copy import copy
from typing import Any, ClassVar, override
import numpy as np
import scipy.special as sp
from scipy.interpolate import CubicSpline
from .._internals import _framework
from .._internals import _references as refs
from ..cosmology import cosmo as csm
from ..exceptions import HMFExtrapolationWarning
from ..halos import mass_definitions as md
[docs]
class SimDetails:
"""
A description of a suite of simulations used to define a mass function.
The parameters given should describe the simulations used to *define*
the mass function in the given study, not all simulations run in the study
against which the fit was compared.
Several parametes take either scalar or list values. These should be provided
consistently, so that a single value refers to all simulations, and all lists
are the same length.
Parameters
----------
L : list of floats
The boxsizes of the simulations [Mpc/h]
N : list of ints
The number of particles in the simulations
halo_finder_type : str
Either "FoF" or "SO"
omegam : float or list of floats
Matter density used.
sigma_8 : float or list of floats
Normalisation used.
halo_overdensity : float
Halo overdensity used (linking length in case of FoF definition)
halo_finder : str, optional
Name of halo finding code.
softening : list of floats, optional
Softening length [kpc/h]
transfer : str or list of str, optional
An identifier for the transfer function calculator used.
z_start : float or list of floats, optional
The starting redshift of the simulation
z_meas : float or 2-tuple, optional
Either the redshift of HMF measurement, or (min,max).
ICS : str or list of str, optional
How the ICS were generated, either "1LPT" or "2LPT"
nmin : int, optional
The minimum number of particles per halo for haloes used in the fit.
hmf_analysis_notes : str, optional
A description of any pertinent details about how the HMF was
analysed in the study.
other_cosmo : dict, optional
Other cosmological parameters of interest.
"""
def __init__(
self,
L,
N,
halo_finder_type,
omegam,
sigma_8,
halo_overdensity,
halo_finder=None,
softening=None,
transfer=None,
z_start=None,
z_meas=None,
ICS=None,
nmin=None,
hmf_analysis_notes="",
other_cosmo=None,
):
# Possible multi-sims
self.L = np.atleast_1d(L)
self.N = np.atleast_1d(N)
self.omegam = np.atleast_1d(omegam)
self.sigma_8 = np.atleast_1d(sigma_8)
self.transfer = np.atleast_1d(transfer)
self.z_start = np.atleast_1d(z_start)
self.softening = np.atleast_1d(softening)
self.ICS = np.atleast_1d(ICS)
self.z_meas = z_meas
self.halo_finder_type = halo_finder_type
self.halo_overdensity = halo_overdensity
self.halo_finder = halo_finder
self.hmf_analysis_notes = hmf_analysis_notes
self.nmin = nmin
self.other_cosmo = other_cosmo or {}
# Derived
self.V = self.L**3
try:
self.mp = self.omegam * 2.7755e11 * self.V / self.N
self.mmin = self.mp * self.nmin
except TypeError:
self.mp = None
self.mmin = None
def _makedoc(
pdocs: str, lname: str, sname: str, eq: str, references: tuple[str, ...], note: str = ""
) -> str:
r"""Build the standard docstring of a fitting function.
A table of the model's default parameters is added automatically when the
class is defined (see :func:`_add_defaults_table`).
Parameters
----------
pdocs
The parameters section, usually ``BaseFittingFunction._pdocs``.
lname
Long name of the fit, used in the summary line.
sname
Short name of the fit, used as a subscript in the equation.
eq
LaTeX form of :math:`f(\sigma)`.
references
The references for the fit. The first is cited in the notes.
note
An optional extra paragraph placed after the summary line.
Returns
-------
str
The docstring.
"""
note = f"\n {note}\n" if note else ""
ref_lines = "\n".join(f" .. [{i}] {ref}" for i, ref in enumerate(references, start=1))
return (
rf"""
{lname} mass function fit.
{note}
For details on attributes, see documentation for :class:`BaseFittingFunction`.
"""
+ pdocs
+ rf"""
Notes
-----
The {lname} [1]_ form is:
.. math:: f_{{\rm {sname}}}(\sigma) = {eq}
References
----------
{ref_lines}
"""
)
[docs]
@_framework.pluggable
class BaseFittingFunction(_framework.Component):
r"""
Base-class for a halo mass function fit.
This class should not be called directly, rather use a subclass which is
specific to a certain fitting formula. The only property necessary to define
for any subclass is :attr:`fsigma`, as well as a dictionary of default
parameters as a class variable ``_defaults``. Model parameters defined here are
accessed through the ``params`` instance attribute (and may be overridden at
instantiation by the user). A subclass may optionally define a :attr:`cutmask`
property, to override the default behaviour of returning True for the whole
range.
In addition, the class attributes :attr:`req_mass` and :attr:`req_neff` say
whether a subclass needs ``m`` and ``n_eff`` to be passed. These must be set
accordingly.
Every subclass is registered as a plugin (unless defined with
``abstract=True``), so it can be used by name, e.g.
``MassFunction(hmf_model="MySMT")``. When a subclass is defined, a table of its
``_defaults`` is added to the Notes section of its docstring, and undocumented
:attr:`fsigma` and :attr:`cutmask` overrides get the docstrings defined here.
Examples
--------
The following defines the Sheth-Mo-Tormen mass function (which is already
included as :class:`SMT`), showing the basic idea of subclassing this class:
>>> from typing import ClassVar
>>> import numpy as np
>>> from hmf.mass_function.fitting_functions import BaseFittingFunction
>>> class MySMT(BaseFittingFunction):
... # Default parameters
... _defaults: ClassVar[dict[str, float]] = {"a": 0.707, "p": 0.3, "A": 0.3222}
...
... @property
... def fsigma(self):
... A = self.params["A"]
... a = self.params["a"]
... p = self.params["p"]
... return (
... A
... * np.sqrt(2.0 * a / np.pi)
... * self.nu
... * np.exp(-(a * self.nu2) / 2.0)
... * (1 + (1.0 / (a * self.nu2)) ** p)
... )
>>> fit = MySMT(nu2=np.array([1.0, 4.0, 9.0]), a=0.75)
>>> fit.params["a"]
0.75
>>> fit.fsigma.shape
(3,)
In that example, we did not specify :attr:`cutmask`, so all elements are
considered to be within the fitted range.
"""
_pdocs = r"""
Parameters
----------
nu2 : array_like
A vector of peak-heights, :math:`\delta_c^2/\sigma^2` corresponding to `m`
m : array_like, optional
A vector of halo masses [units M_sun/h]. Only necessary if :attr:`req_mass`
is True. Typically provides limits of applicability. Must correspond to
`nu2`.
z : float, optional
The redshift. Only used by models that depend on redshift. Default is 0.
n_eff : array_like, optional
The effective spectral index at `m`. Only required if :attr:`req_neff` is True.
mass_definition : :class:`hmf.halos.mass_definitions.MassDefinition` instance
A halo mass definition. Only required for fits which explicitly include a
parameterization for halo definition.
cosmo : :class:`astropy.cosmology.FLRW` instance, optional
A cosmology. Default is :data:`hmf.cosmology.cosmo.DEFAULT_COSMOLOGY`
(Planck18), as for the frameworks. Either `omegam_z` or `cosmo` is required if
:attr:`req_omz` is True. If both are passed, omegam_z takes precedence.
\*\*model_parameters : unpacked-dictionary
These parameters are model-specific. The available parameters and their
defaults are listed in the Notes section of each model (and are also
available as ``<model>._defaults``).
"""
__doc__ += _pdocs
_defaults: ClassVar[dict[str, Any]] = {}
# Subclass requirements
req_neff = False #: Whether `n_eff` is required for this subclass
req_mass = False #: Whether `m` is required for this subclass
sim_definition = None #: Details of the defining simulation, instance of :class:`SimDetails`
normalized = False #: Whether this model is normalized so that all mass is in halos
def __init__(
self,
nu2: np.ndarray,
m: np.ndarray | None = None,
z: float = 0.0,
n_eff: np.ndarray | None = None,
mass_definition: md.BaseMassDefinition | None = None,
cosmo: csm.FLRW = csm.DEFAULT_COSMOLOGY,
delta_c: float = 1.68647,
**model_parameters,
):
super().__init__(**model_parameters)
self.nu2 = nu2
self.z = z
self.n_eff = n_eff
self.mass_definition = mass_definition
self.m = m
self.delta_c = delta_c
self.cosmo = cosmo
# Simple Argument validation
if self.req_mass and m is None:
raise ValueError("This fitting function requires m as well as nu")
if self.req_neff and n_eff is None:
raise ValueError("This fitting function requires n_eff")
self.measured_mass_definition = self.get_measured_mdef()
# Set default mass definition. A generic SO definition has no density of its
# own, so use its preferred definition, as MassFunction does.
if self.mass_definition is None and self.measured_mass_definition is not None:
self.mass_definition = self.measured_mass_definition
if isinstance(self.mass_definition, md.SOGeneric):
self.mass_definition = self.mass_definition.preferred
[docs]
@classmethod
def get_measured_mdef(cls):
"""Get the mass definition used in the defining simulation."""
# Try to set the measured mass definition
measured = None
if cls.sim_definition is not None:
kind = cls.sim_definition.halo_finder_type
delta_h = cls.sim_definition.halo_overdensity
if kind.lower() == "fof":
measured = md.FOF(linking_length=float(delta_h))
elif kind.upper() == "SO":
if delta_h == "vir":
measured = md.SOVirial()
elif delta_h.endswith("c"):
measured = md.SOCritical(
overdensity=float(delta_h[:-1]),
)
elif delta_h.endswith("m"):
measured = md.SOMean(overdensity=float(delta_h[:-1]))
elif delta_h.startswith("*"):
# A Generic SO that will accept any SO definition, but has a
# preferred one.
measured = md.SOGeneric(
preferred=md.from_colossus_name(delta_h.split("(")[-1].split(")")[0])
)
else:
warnings.warn(
"Unrecognized overdensity criterion format. "
"Changing mass definitions will be impossible.",
stacklevel=2,
)
else:
warnings.warn(
"Unknown halo finder type in the sim_definition. "
"Changing mass definitions will be impossible.",
stacklevel=2,
)
return measured
@property
def omegam_z(self):
"""Normalised matter density at current redshift."""
return self.cosmo.Om(self.z)
@property
def nu(self):
r"""The peak height, :math:`\nu = \delta_c/\sigma`."""
return np.sqrt(self.nu2)
@property
def sigma(self):
"""The sqrt of mass variance as a function of mass."""
return self.delta_c / self.nu
@property
def lnsigma(self):
"""Negative log of sigma."""
return -np.log(self.sigma)
@property
def cutmask(self):
r"""
Boolean mask of the elements within the fit's calibrated range.
Has the same length as ``nu2``. An element is True if the corresponding
mass (or :math:`\sigma`, depending on the model) lies inside the range
over which the fit was calibrated against simulations, as given in the
reference paper. Elements outside that range are False: :attr:`fsigma`
is still computed there, but is an extrapolation. The default is True
everywhere. The mask is informational only; it is not applied anywhere
in ``hmf``.
"""
return np.ones(len(self.nu2), dtype=bool)
@property
def fsigma(self):
r"""The function :math:`f(\sigma)\equiv\nu f(\nu)`."""
#: Alias of :class:`BaseFittingFunction`, kept for backwards compatibility.
FittingFunction = BaseFittingFunction
_SECTION_HEADER = re.compile(r"^([A-Z][A-Za-z ]*)\n-+[ \t]*$", flags=re.MULTILINE)
def _defaults_table(defaults: dict[str, Any]) -> str:
"""Render a model's default parameters as an RST table.
Parameters
----------
defaults
The ``_defaults`` dictionary of a model.
Returns
-------
str
An unindented RST paragraph and simple table listing each key and default.
"""
if not defaults:
return "Default model parameters: none (this model has no free parameters)."
keys = [f"``{k}``" for k in defaults]
vals = [repr(v) for v in defaults.values()]
kw = max(len("Parameter"), *map(len, keys))
vw = max(len("Default"), *map(len, vals))
rule = f"{'=' * kw} {'=' * vw}"
rows = [f"{k:<{kw}} {v}" for k, v in zip(keys, vals, strict=True)]
return "\n".join(
[
(
"Default model parameters (override them with keyword arguments, or with "
"``hmf_params`` in :class:`~hmf.mass_function.hmf.MassFunction`):"
),
"",
rule,
f"{'Parameter':<{kw}} Default",
rule,
*rows,
rule,
]
)
def _add_defaults_table(doc: str | None, defaults: dict[str, Any]) -> str:
"""Add a table of default parameters to the Notes section of a docstring.
Parameters
----------
doc
The (numpydoc-style) docstring of a fitting function, possibly None.
defaults
The ``_defaults`` dictionary of the fitting function.
Returns
-------
str
The dedented docstring with the table at the end of its Notes section. If
there is no Notes section, one is created before the References or
Examples section (or at the end).
"""
doc = inspect.cleandoc(doc or "")
table = _defaults_table(defaults)
headers = list(_SECTION_HEADER.finditer(doc))
names = [h.group(1) for h in headers]
if "Notes" in names:
idx = names.index("Notes") + 1
insert = f"{table}\n\n"
else:
idx = next(
(i for i, name in enumerate(names) if name in ("References", "Examples")),
len(headers),
)
insert = f"Notes\n-----\n{table}\n\n"
pos = headers[idx].start() if idx < len(headers) else len(doc)
return f"{doc[:pos].rstrip()}\n\n{insert}{doc[pos:]}".rstrip() + "\n"
_register_plugin = BaseFittingFunction.__init_subclass__.__func__
def _init_fitting_function_subclass(cls: type[BaseFittingFunction], abstract: bool = False) -> None:
"""Register a new fitting function and complete its documentation.
Parameters
----------
cls
The newly-defined subclass.
abstract
If True, the subclass is not registered as a plugin.
"""
_register_plugin(cls, abstract=abstract)
cls.__doc__ = _add_defaults_table(cls.__doc__, cls._defaults)
# Let undocumented overrides inherit the base-class docstrings.
for name in ("fsigma", "cutmask"):
prop = cls.__dict__.get(name)
if isinstance(prop, property) and not prop.__doc__:
prop.__doc__ = getattr(BaseFittingFunction, name).__doc__
# The pluggable decorator sets __init_subclass__ on the class, so it has to be
# extended after the class is created.
BaseFittingFunction.__init_subclass__ = classmethod(_init_fitting_function_subclass)
[docs]
class PS(BaseFittingFunction):
"""Press-Schechter mass function fit."""
# Subclass requirements
req_sigma = False #: Whether sigma is required to compute this model.
req_z = False #: Whether redshift is required for this model.
_eq = r"\sqrt{\frac{2}{\pi}}\nu\exp(-0.5\nu^2)"
references: ClassVar[tuple[str, ...]] = (
(
"Press, W. H., Schechter, P., 1974. ApJ 187, 425-438. "
"http://adsabs.harvard.edu/full/1974ApJ...187..425P"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Press-Schechter", "PS", _eq, references)
normalized = True
@override
@property
def fsigma(self):
return np.sqrt(2.0 / np.pi) * self.nu * np.exp(-0.5 * self.nu2)
[docs]
class SMT(BaseFittingFunction):
"""Sheth-Mo-Tormen mass function fit."""
# Subclass requirements
req_sigma = False
req_z = False
_eq = r"A\sqrt{2a/\pi}\nu\exp(-a\nu^2/2)(1+(a\nu^2)^{-p})"
references: ClassVar[tuple[str, ...]] = (
(
"Sheth, R. K., Mo, H. J., Tormen, G., May 2001. MNRAS 323 (1), 1-12. "
"http://doi.wiley.com/10.1046/j.1365-8711.2001.04006.x"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Sheth-Mo-Tormen", "SMT", _eq, references)
_defaults: ClassVar[dict[str, Any]] = {"a": 0.707, "p": 0.3, "A": None}
normalized = True
sim_definition = SimDetails(
L=[84.5, 141.3],
N=[256**3, 256**3],
halo_finder_type="SO",
omegam=0.3,
sigma_8=0.9,
halo_overdensity="vir",
halo_finder=None,
softening=30.0,
transfer="BondEfs",
z_start=30.0,
z_meas=0.0,
ICS=None,
nmin=None,
hmf_analysis_notes="No details are given about measurement of HMF. ",
other_cosmo={"omegav": 0.7, "h": 0.7, "n": 1},
)
def __init__(self, *args, validate=True, **kwargs):
super().__init__(*args, **kwargs)
if validate:
if self.params["p"] >= 0.5:
raise ValueError(f"p in SMT must be < 0.5. Got {self.params['p']}")
if self.params["a"] <= 0:
raise ValueError(f"a in SMT must be > 0. Got {self.params['a']}.")
@override
@property
def fsigma(self):
A = self._norm()
a = self.params["a"]
p = self.params["p"]
return (
A
* np.sqrt(2.0 * a / np.pi)
* self.nu
* np.exp(-(a * self.nu2) / 2.0)
* (1 + (1.0 / (a * self.nu2)) ** p)
)
def _norm(self):
if self.params["A"] is not None:
return self.params["A"]
p = self.params["p"]
return 1.0 / (1 + 2**-p * sp.gamma(0.5 - p) / sp.gamma(0.5))
[docs]
class ST(SMT):
"""Alias of :class:`SMT`."""
__doc__ = _makedoc(
BaseFittingFunction._pdocs,
"Sheth-Mo-Tormen",
"ST",
SMT._eq,
SMT.references,
note="This is an alias of :class:`SMT` (Sheth-Tormen), kept for backwards compatibility.",
)
[docs]
class Jenkins(BaseFittingFunction):
"""Jenkins mass function fit."""
# Subclass requirements
req_z = False
_eq = r"A\exp\left(-\left|\ln\sigma^{-1}+b\right|^c\right)"
references: ClassVar[tuple[str, ...]] = (
(
"Jenkins, A. R., et al., Feb. 2001. MNRAS 321 (2), 372-384. "
"http://doi.wiley.com/10.1046/j.1365-8711.2001.04029.x"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Jenkins", "Jenkins", _eq, references)
_defaults: ClassVar[dict[str, float]] = {"A": 0.315, "b": 0.61, "c": 3.8}
normalized = False
sim_definition = SimDetails(
L=[84.5, 141.3, 479, 3000],
N=[256**3, 256**3, 134217728, 1000**3],
halo_finder_type="FoF",
omegam=0.3,
sigma_8=0.9,
halo_overdensity=0.2,
halo_finder=None,
softening=30.0,
# Jenkins et al. (2001) Sect. 2.1: the 512^3 and Hubble-volume LCDM runs
# used a CMBFAST transfer function, the others Bond & Efstathiou.
transfer=["BondEfs", "BondEfs", "CMBFAST", "CMBFAST"],
z_start=30.0,
z_meas=(0.0, 5.0),
ICS=None,
nmin=20,
hmf_analysis_notes="""
Many cosmologies used. Preferentially listed LCDM here.
Fit involves "smoothing" and deconvolving HMF.""",
other_cosmo={"omegav": 0.7, "h": 0.7, "n": 1},
)
@override
@property
def cutmask(self):
return np.logical_and(self.lnsigma > -1.2, self.lnsigma < 1.05)
@override
@property
def fsigma(self):
A = self.params["A"]
b = self.params["b"]
c = self.params["c"]
return A * np.exp(-(np.abs(self.lnsigma + b) ** c))
[docs]
class Warren(BaseFittingFunction):
"""Warren mass function fit."""
# Subclass requirements
req_z = False
req_mass = True
_eq = (
r"A\left[\left(\frac{e}{\sigma}\right)^b + c\right]\exp"
r"\left(\frac{d}{\sigma^2}\right)"
)
references: ClassVar[tuple[str, ...]] = (
(
"Warren, M. S., et al., Aug. 2006. ApJ 646 (2), 881-885. "
"http://adsabs.harvard.edu/abs/2006ApJ...646..881W"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Warren", "Warren", _eq, references)
_defaults: ClassVar[dict[str, float]] = {
"A": 0.7234,
"b": 1.625,
"c": 0.2538,
"d": 1.1982,
"e": 1,
}
normalized = False
uncertainties: ClassVar[dict[str, float]] = {
"A": 0.0073,
"a": 0.028,
"b": 0.0051,
"c": 0.0075,
} #: Quoted uncertainties of the model parameters.
sim_definition = SimDetails(
L=[96, 135, 192, 272, 384, 543, 768, 1086, 1536, 2172, 2583, 3072],
N=1024**3,
halo_finder_type="FoF",
omegam=0.3,
sigma_8=0.9,
halo_overdensity=0.2,
halo_finder=None,
softening=[
2.1,
134.0 / 31.0,
192 / 31.0,
272 / 31.0,
384 / 31.0,
543 / 31.0,
768 / 31.0,
1086 / 31.0,
1536 / 31.0,
2172 / 31.0,
2583 / 31.0,
98,
],
transfer="CMBFAST",
z_start=None,
z_meas=(0.0, 5.0),
ICS="1LPT",
nmin=400,
hmf_analysis_notes="FOF N-Correction applied. Fit uses ML of Poisson counts.",
other_cosmo={"omegav": 0.7, "omegab": 0.04, "h": 0.7, "n": 1},
)
@override
@property
def fsigma(self):
A = self.params["A"]
b = self.params["b"]
c = self.params["c"]
d = self.params["d"]
e = self.params["e"]
return A * ((e / self.sigma) ** b + c) * np.exp(-d / self.sigma**2)
@override
@property
def cutmask(self):
return np.logical_and(self.m > 1e10, self.m < 1e15)
[docs]
class Reed03(SMT):
"""Reed 2003 mass function fit."""
# Subclass requirements
req_sigma = True
_eq = r"f_{\rm SMT}(\sigma)\exp\left(-\frac{c}{\sigma \cosh^5(2\sigma)}\right)"
references: ClassVar[tuple[str, ...]] = (
(
"Reed, D., et al., Dec. 2003. MNRAS 346 (2), 565-572. "
"http://adsabs.harvard.edu/abs/2003MNRAS.346..565R"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Reed03", "R03", _eq, references)
_defaults: ClassVar[dict[str, float]] = {"a": 0.707, "p": 0.3, "A": 0.3222, "c": 0.7}
normalized = False
sim_definition = SimDetails(
L=50.0,
N=432**3,
halo_finder_type="FoF",
omegam=0.3,
sigma_8=1.0,
halo_overdensity=0.2,
halo_finder=None,
softening=5.0,
transfer="BBKS",
z_start=[69, 139],
z_meas=(0.0, 15.0),
ICS="1LPT",
nmin=64,
hmf_analysis_notes="HMF seems to be purely binned.",
other_cosmo={"omegav": 0.7, "omegab": 0.04, "h": None, "n": None},
)
@override
@property
def fsigma(self):
vfv = super().fsigma
return vfv * np.exp(-self.params["c"] / (self.sigma * np.cosh(2.0 * self.sigma) ** 5))
@override
@property
def cutmask(self):
return np.logical_and(self.lnsigma > -1.7, self.lnsigma < 0.9)
[docs]
class Reed07(BaseFittingFunction):
"""Reed 2007 mass function fit."""
req_neff = True
req_z = False
_eq = (
r"A\sqrt{2a/\pi}\left[1+(\frac{1}{a\nu^2})^p+0.6G_1+0.4G_2\right]\nu"
r"\exp\left(-ca\nu^2/2-\frac{0.03\nu^{0.6}}{(n_{\rm eff}+3)^2}\right)"
)
references: ClassVar[tuple[str, ...]] = (refs.REED07,)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Reed07", "R07", _eq, references)
_defaults: ClassVar[dict[str, float]] = {"A": 0.3222, "p": 0.3, "c": 1.08, "a": 0.764}
sim_definition = SimDetails(
L=[1.0, 2.5, 2.5, 2.5, 2.5, 4.64, 11.6, 20, 50, 100, 500, 1340, 3000],
N=[
400**3,
1000**3,
1000**3,
500**3,
200**3,
400**3,
1000**3,
400**3,
1000**3,
900**3,
2160**3,
1448**3,
1000**3,
],
halo_finder_type="FoF",
omegam=0.3,
sigma_8=0.9,
halo_overdensity=0.2,
halo_finder=None,
softening=[
0.125,
0.125,
0.125,
0.25,
0.625,
0.58,
0.58,
2.5,
2.4,
2.4,
5.0,
20,
100,
],
transfer="CMBFAST",
z_start=[299, 299, 299, 299, 299, 249, 249, 249, 299, 149, 127, 63, 35],
z_meas=[10, 10, 30, 10, 10, 10, 10, 10, 10, 10, 0, 0, 0],
ICS="1LPT",
nmin=100,
hmf_analysis_notes="Finite volume corrections applied.",
other_cosmo={"omegav": 0.7, "omegab": None, "h": 0.7, "n": 1.0},
)
@override
@property
def fsigma(self):
G_1 = np.exp(-((self.lnsigma - 0.4) ** 2) / (2 * 0.6**2))
G_2 = np.exp(-((self.lnsigma - 0.75) ** 2) / (2 * 0.2**2))
c = self.params["c"]
a = self.params["a"] / self.params["c"]
A = self.params["A"]
p = self.params["p"]
return (
A
* np.sqrt(2.0 * a / np.pi)
* (1.0 + (1.0 / (a * self.nu**2)) ** p + 0.6 * G_1 + 0.4 * G_2)
* self.nu
* np.exp(-c * a * self.nu**2 / 2.0 - 0.03 * self.nu**0.6 / (self.n_eff + 3) ** 2)
)
@override
@property
def cutmask(self):
return np.logical_and(self.lnsigma > -0.5, self.lnsigma < 1.2)
[docs]
class Peacock(BaseFittingFunction):
"""Peacock mass function fit."""
req_z = False
req_mass = True
_eq = r"\nu\exp(-c\nu^2)(2cd\nu+ba\nu^{b-1})/d^2"
references: ClassVar[tuple[str, ...]] = (
(
"Peacock, J. A., Aug. 2007. MNRAS 379 (3), 1067-1074. "
"http://adsabs.harvard.edu/abs/2007MNRAS.379.1067P"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Peacock", "Pck", _eq, references)
_defaults: ClassVar[dict[str, float]] = {"a": 1.529, "b": 0.704, "c": 0.412}
sim_definition = copy(Warren.sim_definition)
sim_definition.hmf_analysis_notes = "Fit directly to Warren+2006 fit."
normalized = True
@override
@property
def fsigma(self):
a = self.params["a"]
b = self.params["b"]
c = self.params["c"]
d = 1 + a * self.nu**b
return (
self.nu
* np.exp(-c * self.nu2)
* (2 * c * d * self.nu + b * a * self.nu ** (b - 1))
/ d**2
)
@override
@property
def cutmask(self):
# Peacock (2007) fits Warren+06, so is calibrated over the same mass range.
return np.logical_and(self.m > 1e10, self.m < 1e15)
[docs]
class Angulo(BaseFittingFunction):
"""Angulo mass function fit."""
req_mass = True
references: ClassVar[tuple[str, ...]] = (
(
"Angulo, R. E., et al., 2012. MNRAS 426, 2046. "
"https://doi.org/10.1111/j.1365-2966.2012.21830.x"
),
)
_eq = r"A \left[\left(\frac{d}{\sigma}\right)^b + 1 \right] \exp(-c/\sigma^2)"
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Angulo", "Ang", _eq, references)
_defaults: ClassVar[dict[str, float]] = {"A": 0.201, "b": 1.7, "c": 1.172, "d": 2.08}
sim_definition = SimDetails(
L=3000.0,
N=6720**3,
halo_finder_type="FoF",
omegam=0.25,
sigma_8=0.9,
halo_overdensity=0.2,
halo_finder=None,
softening=13.79,
transfer="CAMB",
z_start=63,
z_meas=0,
ICS="2LPT",
nmin=20,
hmf_analysis_notes="No corrections seem to be applied; no special techniques.",
other_cosmo={"omegav": 0.75, "omegab": 0.045, "h": 0.73, "n": 1.0},
)
@override
@property
def fsigma(self):
A = self.params["A"]
b = self.params["b"]
c = self.params["c"]
d = self.params["d"]
return A * ((d / self.sigma) ** b + 1) * np.exp(-c / self.sigma**2)
@override
@property
def cutmask(self):
return np.logical_and(self.m > 1e8, self.m < 1e16)
[docs]
class AnguloBound(Angulo):
"""Bounded version of Angulo mass function fit."""
__doc__ = _makedoc(
BaseFittingFunction._pdocs,
"Bounded Angulo",
"Ang",
Angulo._eq,
Angulo.references,
note="Same form as :class:`Angulo`, with the alternative (bound) parameter set.",
)
_defaults: ClassVar[dict[str, float]] = {"A": 0.265, "b": 1.9, "c": 1.4, "d": 1.675}
[docs]
class Watson_FoF(Warren):
"""Watson friend-of-friend mass function fit."""
req_mass = False
references: ClassVar[tuple[str, ...]] = (
("Watson, W. A., et al., 2013. MNRAS 433, 1230. https://doi.org/10.1093/mnras/stt791"),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Watson FoF", "WatF", Warren._eq, references)
_defaults: ClassVar[dict[str, float]] = {
"A": 0.282,
"b": 2.163,
"c": 1,
"d": 1.21,
"e": 1.406,
}
sim_definition = SimDetails(
L=[11.4, 20, 114, 425, 1000, 3200, 6000],
N=[3072**3, 5488**3, 3072**3, 5488**3, 3456**3, 4000**3, 6000**3],
halo_finder_type="FoF",
omegam=0.27,
sigma_8=0.8,
halo_overdensity=0.2,
halo_finder="GADGET3",
softening=[0.18, 0.18, 1.86, 3.87, 14.47, 40.0, 50.0],
transfer="CAMB",
z_start=[300, 300, 300, 300, 150, 120, 100],
z_meas=(0, 30),
ICS="1LPT",
nmin=1000,
hmf_analysis_notes="Warren FOF correction applied. Finite-box correction applied.",
other_cosmo={"omegav": 0.73, "omegab": 0.044, "h": 0.7, "n": 0.96},
)
@override
@property
def cutmask(self):
return np.logical_and(self.lnsigma > -0.55, self.lnsigma < 1.31)
[docs]
class Watson(BaseFittingFunction):
"""Watson mass function fit."""
req_cosmo = True
req_dhalo = True
req_omz = True
references: ClassVar[tuple[str, ...]] = Watson_FoF.references
_eq = r"\Gamma A \left(\left(\frac{\beta}{\sigma}\right)^\alpha+1\right)\exp(-\gamma/\sigma^2)"
_note = r"""
Three AHF fits from Watson+13 (v4/published) are used, by redshift. At
:math:`z=0` it uses the precise present-day fit (Table 2: ``A_0``, ``alpha_0``,
``beta_0``, ``gamma_0``). At :math:`z \geq` ``z_hi`` it uses the EoR fit (Table 2:
``*_hi``). In between it uses the redshift-dependent fit of eqs 14-16, with
:math:`\alpha(z) = \Omega_m(z)\{5.907(1+z)^{-3.599} + 2.344\}`,
:math:`\beta(z) = \Omega_m(z)\{3.136(1+z)^{-3.058} + 2.349\}`,
:math:`A(z) = \Omega_m(z)\{1.097(1+z)^{-3.216} + 0.074\}` and :math:`\gamma` =
``gamma_z``. (arXiv v1 swaps the labels of the :math:`\alpha` and :math:`\beta`
fits; v4 is the published version.)
As in the paper, the fits are not joined smoothly, so :math:`f(\sigma)` jumps
between :math:`z=0` and :math:`z \to 0^+`. For the simulation's
:math:`\Omega_m = 0.27` the jump is 0-9% over
:math:`-0.55 \leq \ln\sigma^{-1} \leq 1.05`, within the paper's ~10% accuracy.
Because :math:`\Omega_m(z)` multiplies :math:`A`, :math:`\alpha` and
:math:`\beta`, and the fit was calibrated at a single cosmology, it extrapolates
poorly to other :math:`\Omega_m`: at :math:`\Omega_m = 0.31` the jump is a factor
1.4-2.7. Bocquet+16 call this a "heightened and likely artificial cosmological
sensitivity".
"""
__doc__ = _makedoc(
BaseFittingFunction._pdocs, "Watson", "WatS", _eq, Watson_FoF.references, note=_note.strip()
)
sim_definition = copy(Watson_FoF.sim_definition)
sim_definition.halo_finder_type = "SO"
sim_definition.halo_finder = "AHF"
sim_definition.halo_overdensity = "*(vir)"
_defaults: ClassVar[dict[str, Any]] = {
"C_a": 0.023,
"d_a": 0.456,
"d_b": 0.139,
"p": 0.072,
"q": 2.13,
"A_0": 0.194,
"alpha_0": 1.805,
"beta_0": 2.267,
"gamma_0": 1.287,
"z_hi": 6,
"A_hi": 0.563,
"alpha_hi": 3.810,
"beta_hi": 0.874,
"gamma_hi": 1.453,
"A_a": 1.097,
"A_b": 3.216,
"A_c": 0.074,
"alpha_a": 5.907,
"alpha_b": 3.599,
"alpha_c": 2.344,
"beta_a": 3.136,
"beta_b": 3.058,
"beta_c": 2.349,
"gamma_z": 1.318,
}
[docs]
def gamma(self):
r"""Calculate :math:`\Gamma` for the Watson fit."""
if self.mass_definition is None:
delta_halo = 178.0
elif not isinstance(self.mass_definition, md.SphericalOverdensity):
raise ValueError("The Watson fitting function is a spherical-overdensity function.")
else:
delta_halo = self.mass_definition.halo_overdensity_mean(self.z, self.cosmo)
C = np.exp(self.params["C_a"] * (delta_halo / 178 - 1))
d = -self.params["d_a"] * self.omegam_z - self.params["d_b"]
p = self.params["p"]
q = self.params["q"]
return C * (delta_halo / 178) ** d * np.exp(p * (1 - delta_halo / 178) / self.sigma**q)
[docs]
def get_params(self):
r"""Get :math:`A`, :math:`\alpha`, :math:`\beta` and :math:`\gamma` at this redshift."""
if self.z == 0:
A = self.params["A_0"]
alpha = self.params["alpha_0"]
beta = self.params["beta_0"]
gamma = self.params["gamma_0"]
elif self.z >= self.params["z_hi"]:
A = self.params["A_hi"]
alpha = self.params["alpha_hi"]
beta = self.params["beta_hi"]
gamma = self.params["gamma_hi"]
else:
omz = self.omegam_z
A = omz * (
self.params["A_a"] * (1 + self.z) ** (-self.params["A_b"]) + self.params["A_c"]
)
alpha = omz * (
self.params["alpha_a"] * (1 + self.z) ** (-self.params["alpha_b"])
+ self.params["alpha_c"]
)
beta = omz * (
self.params["beta_a"] * (1 + self.z) ** (-self.params["beta_b"])
+ self.params["beta_c"]
)
gamma = self.params["gamma_z"]
return A, alpha, beta, gamma
@override
@property
def fsigma(self):
A, alpha, beta, gamma = self.get_params()
return (
self.gamma() * A * ((beta / self.sigma) ** alpha + 1) * np.exp(-gamma / self.sigma**2)
)
@override
@property
def cutmask(self):
return np.logical_and(self.lnsigma > -0.55, self.lnsigma < 1.05)
[docs]
class Crocce(Warren):
"""Crocce mass function fit."""
req_z = True
references: ClassVar[tuple[str, ...]] = (
(
"Crocce, M., et al., 2010. MNRAS 403, 1353. "
"https://doi.org/10.1111/j.1365-2966.2009.16194.x"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Crocce", "Cro", Warren._eq, references)
_defaults: ClassVar[dict[str, Any]] = {
"A_a": 0.58,
"A_b": 0.13,
"b_a": 1.37,
"b_b": 0.15,
"c_a": 0.3,
"c_b": 0.084,
"d_a": 1.036,
"d_b": 0.024,
"e": 1,
}
sim_definition = SimDetails(
L=[7680, 3072, 4500, 768, 384, 179],
N=[2048**3, 2048**3, 1200**3, 1024**3, 1024**3, 1024**3],
halo_finder_type="FoF",
omegam=0.25,
sigma_8=0.8,
halo_overdensity=0.2,
halo_finder=None,
softening=[50, 50, 100, 50, 50, 50],
transfer="CAMB",
z_start=[150, 50, 50, 50, 50, 50],
z_meas=(0, 1),
ICS=["1LPT", "1LPT", "2LPT", "2LPT", "2LPT", "2LPT"],
nmin=200,
hmf_analysis_notes="Warren FOF correction applied.",
other_cosmo={"omegav": 0.75, "omegab": 0.044, "h": 0.7, "n": 0.95},
)
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
self.params["A"] = self.params["A_a"] * (1 + self.z) ** (-self.params["A_b"])
self.params["b"] = self.params["b_a"] * (1 + self.z) ** (-self.params["b_b"])
self.params["c"] = self.params["c_a"] * (1 + self.z) ** (-self.params["c_b"])
self.params["d"] = self.params["d_a"] * (1 + self.z) ** (-self.params["d_b"])
@override
@property
def cutmask(self):
return np.logical_and(self.m > 10**10.5, self.m < 10**15.5)
[docs]
class Courtin(SMT):
"""Courtin mass function fit."""
req_sigma = True
references: ClassVar[tuple[str, ...]] = (
(
"Courtin, J., et al., 2011. MNRAS 410, 1911. "
"https://doi.org/10.1111/j.1365-2966.2010.17573.x"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Courtin", "Ctn", SMT._eq, references)
_defaults: ClassVar[dict[str, float]] = {"A": 0.348, "a": 0.695, "p": 0.1}
normalized = False
sim_definition = SimDetails(
L=[162, 648, 1296],
N=[512**3, 512**3, 512**3],
halo_finder_type="FoF",
omegam=0.26,
sigma_8=0.79,
halo_overdensity=0.2,
halo_finder=None,
softening=[2.47, 19.78, 39.55],
transfer="CAMB",
z_start=[93, 56, 41],
z_meas=0,
ICS="1LPT",
nmin=200,
hmf_analysis_notes="Many systematic effects tested but not applied.",
other_cosmo={"omegav": 0.74, "omegab": 0.044, "h": 0.72, "n": 0.963},
)
@override
@property
def cutmask(self):
return np.logical_and(self.lnsigma > -0.8, self.lnsigma < 0.7)
[docs]
class Bhattacharya(SMT):
"""Bhattacharya mass function fit."""
req_z = True
req_mass = True
_eq = r"f_{\rm SMT}(\sigma) (\nu\sqrt{a})^{q-1}"
references: ClassVar[tuple[str, ...]] = (
(
"Bhattacharya, S., et al., May 2011. ApJ 732 (2), 122. "
"http://labs.adsabs.harvard.edu/ui/abs/2011ApJ...732..122B"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Bhattacharya", "Btc", _eq, references)
_defaults: ClassVar[dict[str, Any]] = {
"A_a": 0.333,
"A_b": 0.11,
"a_a": 0.788,
"a_b": 0.01,
"p": 0.807,
"q": 1.795,
"normed": False,
}
normalized = False
sim_definition = SimDetails(
L=[1000 * 0.72, 1736 * 0.72, 2778 * 0.72, 178 * 0.72, 1300 * 0.72],
N=[1500**3, 1200**3, 1024**3, 512**3, 1024**3],
halo_finder_type="FoF",
omegam=None, # what is lower case omega??
sigma_8=0.8,
halo_overdensity=0.2,
halo_finder=None,
softening=[24, 51, 97, 14, 50],
transfer="CAMB",
z_start=[75, 100, 100, 211, 211],
z_meas=(0, 2),
ICS=["2LPT", "2LPT", "2LPT", "1LPT", "1LPT"],
nmin=400,
hmf_analysis_notes="Finite force correction. FOF Correction. Finite volume correction.",
other_cosmo={
"omegav": 0.74,
"omegab": None, # uses lower case omega without definition
"h": 0.72,
"n": 0.97,
},
)
def __init__(self, **kwargs):
super().__init__(validate=False, **kwargs)
if not self.params["normed"]:
self.params["A"] = self.params["A_a"] * (1 + self.z) ** -self.params["A_b"]
else:
self.params["A"] = self._norm()
self.params["a"] = self.params["a_a"] * (1 + self.z) ** -self.params["a_b"]
# To enable satisfying normalization to unity
if self.params["q"] <= 0:
raise ValueError("q in Bhattacharya must be > 0")
if self.params["p"] * 2 >= self.params["q"]:
raise ValueError("2p in Bhattacharya must be < q")
@property
@override
def fsigma(self):
r"""
Calculate :math:`f(\sigma)` for Bhattacharya form.
Bhattacharya, S., et al., May 2011. ApJ 732 (2), 122.
http://labs.adsabs.harvard.edu/ui/abs/2011ApJ...732..122B
.. note:: valid for :math:`10^{11.8}M_\odot < M <10^{15.5}M_\odot`
Returns
-------
vfv : array_like, len=len(pert.M)
The function :math:`f(\sigma)\equiv\nu f(\nu)`.
"""
vfv = super().fsigma
return vfv * (np.sqrt(self.params["a"]) * self.nu) ** (self.params["q"] - 1)
@override
@property
def cutmask(self):
return np.logical_and(self.m > 6 * 10**11, self.m < 3 * 10**15)
@override
def _norm(self):
if self.params.get("A") is not None:
return self.params["A"]
# A is the reciprocal of the integral of f/A over ln(1/sigma).
p, q = self.params["p"], self.params["q"]
return np.sqrt(np.pi) / (
2 ** (-1 / 2 - p + q / 2) * (2**p * sp.gamma(q / 2) + sp.gamma(-p + q / 2))
)
def _warn_if_delta_extrapolated(fit: BaseFittingFunction, delta_halo: float) -> None:
r"""Warn if ``delta_halo`` is outside the calibrated overdensity range of a Tinker fit.
The Tinker08/Tinker10 parameters are tabulated for :math:`200 \le \Delta_m \le 3200`
(Table 2 of Tinker et al. 2008; Table 4 of Tinker et al. 2010) and spline-interpolated
in between. Outside that range the spline extrapolates.
"""
lo, hi = fit.delta_virs[0], fit.delta_virs[-1]
if not lo * (1 - 1e-8) <= delta_halo <= hi * (1 + 1e-8):
warnings.warn(
f"{fit.__class__.__name__} is calibrated for {lo} <= Delta_mean <= {hi}, but "
f"Delta_mean={delta_halo:.4g} at z={fit.z}. The fit parameters are being "
"extrapolated.",
HMFExtrapolationWarning,
stacklevel=3,
)
def _log_delta_spline(delta_virs: np.ndarray, values: np.ndarray):
r"""Interpolate Tinker fit parameters between the tabulated overdensities.
Uses a natural cubic spline in :math:`\log_{10}\Delta`, as recommended in App. B
of Tinker et al. (2008); with the tabulated parameters it reproduces their
Table B3 of second derivatives.
Returns a callable of :math:`\Delta`.
"""
spl = CubicSpline(np.log10(delta_virs), values, bc_type="natural")
return lambda delta: spl(np.log10(delta))
[docs]
class Tinker08(BaseFittingFunction):
"""Tinker 2008 mass function fit."""
req_z = True
req_dhalo = True
_eq = r"A\left(\frac{\sigma}{b}^{-a}+1\right)\exp(-c/\sigma^2)"
references: ClassVar[tuple[str, ...]] = (
(
"Tinker, J., et al., 2008. ApJ 688, 709-728. "
"http://iopscience.iop.org/0004-637X/688/2/709"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Tinker08", "Tkr", _eq, references)
sim_definition = SimDetails(
L=[
768,
384,
271,
192,
96,
1280,
500,
250,
120,
80,
1000,
500,
500,
500,
384,
384,
120,
80,
],
N=[
1024**3,
1024**3,
1024**3,
1024**3,
1024**3,
640**3,
1024**3,
512**3,
512**3,
512**3,
1024**3,
512**3,
512**3,
512**3,
1024**3,
1024**3,
1024**3,
512**3,
],
halo_finder_type="SO",
omegam=[
0.3,
0.3,
0.3,
0.3,
0.3,
0.27,
0.3,
0.3,
0.3,
0.3,
0.27,
0.24,
0.24,
0.24,
0.26,
0.2,
0.27,
0.23,
],
sigma_8=[
0.9,
0.9,
0.9,
0.9,
0.9,
0.9,
0.9,
0.9,
0.9,
0.9,
0.79,
0.75,
0.75,
0.8,
0.75,
0.9,
0.79,
0.75,
],
halo_overdensity="*(200m)",
halo_finder=None,
softening=[
25,
14,
10,
4.9,
1.4,
120,
15,
7.6,
1.8,
1.2,
30,
15,
15,
15,
14,
14,
0.9,
1.2,
],
transfer=None,
z_start=[
40,
48,
51,
54,
65,
49,
40,
49,
49,
49,
60,
40,
40,
40,
35,
42,
100,
49,
],
z_meas=(0, 2.5),
ICS="1LPT",
nmin=None,
hmf_analysis_notes="No corrections applied.",
other_cosmo={
"omegav": [
0.7,
0.7,
0.7,
0.7,
0.7,
0.73,
0.7,
0.7,
0.7,
0.7,
0.73,
0.76,
0.76,
0.76,
0.74,
0.8,
0.73,
0.77,
],
"omegab": [
0.04,
0.04,
0.04,
0.04,
0.04,
0.04,
0.045,
0.04,
0.04,
0.04,
0.044,
0.042,
0.042,
0.042,
0.042,
0.044,
0.04,
0.044,
0.04,
],
"h": [
0.7,
0.7,
0.7,
0.7,
0.7,
0.7,
0.7,
0.7,
0.7,
0.7,
0.7,
0.73,
0.73,
0.73,
0.71,
0.7,
0.7,
0.73,
],
"n": [
1,
1,
1,
1,
1,
1,
1,
1,
1,
1,
0.95,
0.95,
0.95,
0.95,
0.94,
1,
0.95,
0.95,
],
},
)
_defaults: ClassVar[dict[str, Any]] = { # -- A
"A_200": 1.858659e-01,
"A_300": 1.995973e-01,
"A_400": 2.115659e-01,
"A_600": 2.184113e-01,
"A_800": 2.480968e-01,
"A_1200": 2.546053e-01,
"A_1600": 2.600000e-01,
"A_2400": 2.600000e-01,
"A_3200": 2.600000e-01,
# -- a
"a_200": 1.466904,
"a_300": 1.521782,
"a_400": 1.559186,
"a_600": 1.614585,
"a_800": 1.869936,
"a_1200": 2.128056,
"a_1600": 2.301275,
"a_2400": 2.529241,
"a_3200": 2.661983,
# --- b
"b_200": 2.571104,
"b_300": 2.254217,
"b_400": 2.048674,
"b_600": 1.869559,
"b_800": 1.588649,
"b_1200": 1.507134,
"b_1600": 1.464374,
"b_2400": 1.436827,
"b_3200": 1.405210,
# --- c
"c_200": 1.193958,
"c_300": 1.270316,
"c_400": 1.335191,
"c_600": 1.446266,
"c_800": 1.581345,
"c_1200": 1.795050,
"c_1600": 1.965613,
"c_2400": 2.237466,
"c_3200": 2.439729,
# -- others
"A_exp": 0.14,
"a_exp": 0.06,
}
delta_virs = np.array([200, 300, 400, 600, 800, 1200, 1600, 2400, 3200])
#: Whether to warn when the halo overdensity is outside the calibrated range.
warn_delta_extrapolation: ClassVar[bool] = True
def __init__(self, **model_parameters):
super().__init__(**model_parameters)
if not isinstance(self.mass_definition, md.SphericalOverdensity):
# Kept as ValueError (not TypeError): part of the public API contract, asserted
# verbatim (via Tinker08/Tinker10, which share this __init__) by
# tests/test_fitting_functions_extra.py::test_tinker08_non_so_raises and
# ::test_tinker10_non_so_raises.
raise ValueError( # noqa: TRY004
"The Tinker fitting function is a spherical-overdensity function."
)
delta_halo = self.mass_definition.halo_overdensity_mean(self.z, self.cosmo)
if self.warn_delta_extrapolation:
_warn_if_delta_extrapolated(self, delta_halo)
if delta_halo not in self.delta_virs:
A_array = np.array([self.params[f"A_{d}"] for d in self.delta_virs])
a_array = np.array([self.params[f"a_{d}"] for d in self.delta_virs])
b_array = np.array([self.params[f"b_{d}"] for d in self.delta_virs])
c_array = np.array([self.params[f"c_{d}"] for d in self.delta_virs])
A_func = _log_delta_spline(self.delta_virs, A_array)
a_func = _log_delta_spline(self.delta_virs, a_array)
b_func = _log_delta_spline(self.delta_virs, b_array)
c_func = _log_delta_spline(self.delta_virs, c_array)
A_0 = A_func(delta_halo)
a_0 = a_func(delta_halo)
b_0 = b_func(delta_halo)
c_0 = c_func(delta_halo)
else:
A_0 = self.params[f"A_{int(delta_halo)}"]
a_0 = self.params[f"a_{int(delta_halo)}"]
b_0 = self.params[f"b_{int(delta_halo)}"]
c_0 = self.params[f"c_{int(delta_halo)}"]
self.A = A_0 * (1 + self.z) ** (-self.params["A_exp"])
self.a = a_0 * (1 + self.z) ** (-self.params["a_exp"])
if self.z == 0:
self.b = b_0
elif delta_halo > 75.0:
# Tinker08 eq. 8, which is only defined for Delta > 75.
alpha = 10 ** (-((0.75 / np.log10(delta_halo / 75.0)) ** 1.2))
self.b = b_0 * (1 + self.z) ** (-alpha)
else:
raise ValueError(
f"The {self.__class__.__name__} redshift evolution of b is undefined for "
f"Delta_mean <= 75, got Delta_mean={delta_halo:.4g} at z={self.z}."
)
self.c = c_0
# The spline in Delta extrapolates beyond the calibrated range, and far enough
# out the parameters become unphysical (e.g. A < 0 for Delta of a few 10^4,
# which makes f(sigma) negative; b < 0 would make it NaN).
for name in ("A", "a", "b", "c"):
val = getattr(self, name)
if not np.isfinite(val) or val <= 0:
raise ValueError(
f"{self.__class__.__name__} parameter {name}={val} is unphysical (must be "
f"> 0) for Delta_mean={delta_halo:.4g} at z={self.z}. The fit is "
f"calibrated for {self.delta_virs[0]} <= Delta_mean <= "
f"{self.delta_virs[-1]}."
)
@override
@property
def fsigma(self):
return self.A * ((self.sigma / self.b) ** (-self.a) + 1) * np.exp(-self.c / self.sigma**2)
@override
@property
def cutmask(self):
if self.z == 0.0:
return np.logical_and(self.lnsigma / np.log(10) > -0.6, self.lnsigma / np.log(10) < 0.4)
return np.logical_and(self.lnsigma / np.log(10) > -0.2, self.lnsigma / np.log(10) < 0.4)
[docs]
class Tinker10(BaseFittingFunction):
"""Tinker 2010 mass function fit."""
req_z = True
req_dhalo = True
_eq = r"(1+(\beta\nu)^{-2\phi})\nu^{2\eta+1}\exp(-\gamma\nu^2/2)"
references: ClassVar[tuple[str, ...]] = (
(
"Tinker, J., et al., 2010. ApJ 724, 878. "
"http://iopscience.iop.org/0004-637X/724/2/878/pdf/apj_724_2_878.pdf"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Tinker10", "Tkr", _eq, references)
sim_definition = copy(Tinker08.sim_definition)
_defaults: ClassVar[dict[str, Any]] = { # --- alpha
"alpha_200": 0.368,
"alpha_300": 0.363,
"alpha_400": 0.385,
"alpha_600": 0.389,
"alpha_800": 0.393,
"alpha_1200": 0.365,
"alpha_1600": 0.379,
"alpha_2400": 0.355,
"alpha_3200": 0.327,
# --- beta
"beta_200": 0.589,
"beta_300": 0.585,
"beta_400": 0.544,
"beta_600": 0.543,
"beta_800": 0.564,
"beta_1200": 0.623,
"beta_1600": 0.637,
"beta_2400": 0.673,
"beta_3200": 0.702,
# --- gamma
"gamma_200": 0.864,
"gamma_300": 0.922,
"gamma_400": 0.987,
"gamma_600": 1.09,
"gamma_800": 1.2,
"gamma_1200": 1.34,
"gamma_1600": 1.5,
"gamma_2400": 1.68,
"gamma_3200": 1.81,
# --- phi
"phi_200": -0.729,
"phi_300": -0.789,
"phi_400": -0.910,
"phi_600": -1.05,
"phi_800": -1.2,
"phi_1200": -1.26,
"phi_1600": -1.45,
"phi_2400": -1.5,
"phi_3200": -1.49,
# -- eta
"eta_200": -0.243,
"eta_300": -0.261,
"eta_400": -0.261,
"eta_600": -0.273,
"eta_800": -0.278,
"eta_1200": -0.301,
"eta_1600": -0.301,
"eta_2400": -0.319,
"eta_3200": -0.336,
# --others
"beta_exp": 0.2,
"phi_exp": -0.08,
"eta_exp": 0.27,
"gamma_exp": -0.01,
"max_z": 3,
}
delta_virs = np.array([200, 300, 400, 600, 800, 1200, 1600, 2400, 3200])
terminate = True
normalized = True
def __init__(self, **model_parameters):
super().__init__(**model_parameters)
if self.mass_definition is None:
delta_halo = 200
elif not isinstance(self.mass_definition, md.SphericalOverdensity):
raise ValueError("The Tinker10 fitting function is a spherical-overdensity function.")
else:
delta_halo = self.mass_definition.halo_overdensity_mean(self.z, self.cosmo)
self.delta_halo = delta_halo
_warn_if_delta_extrapolated(self, delta_halo)
if int(delta_halo) not in self.delta_virs:
beta_array = np.array([self.params[f"beta_{d}"] for d in self.delta_virs])
gamma_array = np.array([self.params[f"gamma_{d}"] for d in self.delta_virs])
phi_array = np.array([self.params[f"phi_{d}"] for d in self.delta_virs])
eta_array = np.array([self.params[f"eta_{d}"] for d in self.delta_virs])
beta_func = _log_delta_spline(self.delta_virs, beta_array)
gamma_func = _log_delta_spline(self.delta_virs, gamma_array)
phi_func = _log_delta_spline(self.delta_virs, phi_array)
eta_func = _log_delta_spline(self.delta_virs, eta_array)
beta_0 = beta_func(delta_halo)
gamma_0 = gamma_func(delta_halo)
phi_0 = phi_func(delta_halo)
eta_0 = eta_func(delta_halo)
else:
beta_0 = self.params[f"beta_{int(delta_halo)}"]
gamma_0 = self.params[f"gamma_{int(delta_halo)}"]
phi_0 = self.params[f"phi_{int(delta_halo)}"]
eta_0 = self.params[f"eta_{int(delta_halo)}"]
self.beta = beta_0 * (1 + min(self.z, self.params["max_z"])) ** self.params["beta_exp"]
self.phi = phi_0 * (1 + min(self.z, self.params["max_z"])) ** self.params["phi_exp"]
self.eta = eta_0 * (1 + min(self.z, self.params["max_z"])) ** self.params["eta_exp"]
self.gamma = gamma_0 * (1 + min(self.z, self.params["max_z"])) ** self.params["gamma_exp"]
# The normalisation only works with specific conditions
# gamma > 0
if self.gamma <= 0:
if self.terminate:
raise ValueError("gamma must be > 0, got " + str(self.gamma))
self.gamma = 1e-3
# eta >-0.5
if self.eta <= -0.5:
if self.terminate:
raise ValueError("eta must be > -0.5, got " + str(self.eta))
self.eta = -0.499
# eta-phi >-0.5
if self.eta - self.phi <= -0.5:
if self.terminate:
raise ValueError("eta-phi must be > -0.5, got " + str(self.eta - self.phi))
self.phi = self.eta + 0.499
if self.beta <= 0:
if self.terminate:
raise ValueError("beta must be > 0, got " + str(self.beta))
self.beta = 1e-3
@property
def _normalise(self):
if int(self.delta_halo) in self.delta_virs and self.z == 0:
return self.params[f"alpha_{int(self.delta_halo)}"]
return 1 / (
2 ** (self.eta - self.phi - 0.5)
* self.beta ** (-2 * self.phi)
* self.gamma ** (-0.5 - self.eta)
* (
2**self.phi * self.beta ** (2 * self.phi) * sp.gamma(self.eta + 0.5)
+ self.gamma**self.phi * sp.gamma(0.5 + self.eta - self.phi)
)
)
@override
@property
def fsigma(self):
fv = (
(1 + (self.beta * self.nu) ** (-2 * self.phi))
* self.nu ** (2 * self.eta)
* np.exp(-self.gamma * (self.nu**2) / 2)
)
return fv * self._normalise * self.nu
@override
@property
def cutmask(self):
if self.z == 0.0:
return np.logical_and(self.lnsigma / np.log(10) > -0.6, self.lnsigma / np.log(10) < 0.4)
return np.logical_and(self.lnsigma / np.log(10) > -0.2, self.lnsigma / np.log(10) < 0.4)
[docs]
class Behroozi(Tinker08):
"""Behroozi mass function fit."""
references: ClassVar[tuple[str, ...]] = (
(
"Behroozi, P. S., Wechsler, R. H., Conroy, C., 2013. ApJ 770, 57. "
"https://doi.org/10.1088/0004-637X/770/1/57"
),
)
__doc__ = rf"""
Behroozi mass function fit [1]_.
This is an empirical modification to the :class:`Tinker08` fit, to improve
accuracy at high redshift.
{BaseFittingFunction._pdocs}
Notes
-----
:math:`f(\sigma)` is that of :class:`Tinker08`. The mass function is then
modified (Appendix G of [1]_; masses in :math:`M_\odot`) as
.. math::
\frac{{dn}}{{dM}} = \theta \frac{{dn_{{\rm T08}}}}{{dM}}
- n_{{\rm T08}}(>M) \frac{{d\theta}}{{dM}},
where
.. math::
\theta &= 10^{{\alpha (M/M_\star)^\gamma}}, \\
\alpha &= \frac{{0.144}}{{1 + \exp[14.79 (a - 0.213)]}}, \\
\gamma &= \frac{{0.5}}{{1 + \exp(6.5 a)}},
with :math:`a = 1/(1+z)` and :math:`M_\star = 10^{{11.5}} M_\odot`.
References
----------
.. [1] {references[0]}
"""
normalized = False
# Behroozi+13 deliberately applies Tinker08 at the virial overdensity at all z (which
# drops below Delta_mean=200 at z >~ 1), and calibrates its correction on top of that.
warn_delta_extrapolation = False
sim_definition = SimDetails(
L=[250, 1000, 420],
N=[2048**3, 2048**3, 1400**3],
halo_finder_type="SO",
omegam=0.27,
sigma_8=0.82,
halo_overdensity="vir",
halo_finder="Rockstar",
softening=[1, 7, 8],
transfer="CAMB",
z_start=None,
z_meas=(0, 8),
ICS=["1LPT", "1LPT", "2LPT"],
nmin=None,
hmf_analysis_notes="No corrections applied.",
other_cosmo={
"omegav": 0.73,
"omegab": None, # uses lower case omega without definition
"h": 0.7,
"n": 0.95,
},
)
def _modify_dndm(
self,
m: np.ndarray,
dndm: np.ndarray,
z: float,
ngtm_tinker: np.ndarray,
h: float = 1.0,
):
"""
Apply modifications to dndm in Appendix G of Behroozi+13.
Parameters
----------
m
Halo masses in Msun, NOT Msun/h (the fit's pivot mass is in Msun).
dndm
The Tinker08 mass function, per unit mass in Msun/h.
z
Redshift.
ngtm_tinker
The cumulative Tinker08 mass function, n(>m).
h
Hubble parameter, used to take dtheta/dM per unit Msun/h, consistent
with ``dndm``.
"""
a = 1 / (1 + z)
alpha = 0.144 / (1 + np.exp(14.79 * (a - 0.213))) # Behroozi+13 Eq G2.
gamma = 0.5 / (1 + np.exp(6.5 * a))
mstar = 10**11.5
mscale = (m / mstar) ** gamma
theta = 10 ** (alpha * mscale) # 10^(Eq G2), the multiplicative factor in Eq G3
# Derivative per unit mass in Msun/h, the units of dndm.
dthetadM = theta * np.log(10) * alpha * gamma * mscale / (m * h)
# if ngtm_tinker is very small (ie. 0), dthetadM will be nan.
res = dndm * theta - ngtm_tinker * dthetadM
res[np.isnan(res)] = 0
return res
[docs]
class Pillepich(Warren):
"""Pillepich mass function fit."""
references: ClassVar[tuple[str, ...]] = (
(
"Pillepich, A., Porciani, C., Hahn, O., 2010. MNRAS 402, 191. "
"https://doi.org/10.1111/j.1365-2966.2009.15914.x"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Pillepich", "Pillepich", Warren._eq, references)
_defaults: ClassVar[dict[str, float]] = {
"A": 0.6853,
"b": 1.868,
"c": 0.3324,
"d": 1.2266,
"e": 1,
}
normalized = False
sim_definition = SimDetails(
L=[1200, 1200, 150],
N=[1024**3, 1024**3, 1024**3],
halo_finder_type="FoF",
omegam=[0.279, 0.24, 0.279],
sigma_8=[0.817, 0.76, 0.817],
halo_overdensity=0.2,
halo_finder=None,
softening=[20, 20, 3],
transfer="LINGER",
z_start=[50, 50, 70],
z_meas=0,
ICS="1LPT",
nmin=100,
hmf_analysis_notes="No corrections applied.",
other_cosmo={
"omegav": [0.721, 0.76, 0.721],
"omegab": [0.0462, 0.042, 0.0462],
# uses lower case omega without definition
"h": [0.701, 0.73, 0.701],
},
)
[docs]
class Manera(SMT):
"""Manera mass function fit."""
references: ClassVar[tuple[str, ...]] = (
(
"Manera, M., Sheth, R. K., Scoccimarro, R., 2010. MNRAS 402, 589. "
"https://doi.org/10.1111/j.1365-2966.2009.15921.x"
),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Manera", "Man", SMT._eq, references)
# These are for z=0, new ML method, l_linnk = 0.2
_defaults: ClassVar[dict[str, Any]] = {"A": None, "a": 0.709, "p": 0.289}
sim_definition = SimDetails(
L=1280.0,
N=640**3,
halo_finder_type="FoF",
omegam=0.27,
sigma_8=0.9,
halo_overdensity=0.2,
halo_finder=None,
softening=20,
transfer="CMBFAST",
z_start=50,
z_meas=(0, 0.5),
ICS="2LPT",
nmin=105,
hmf_analysis_notes="FOF Correction applied.",
other_cosmo={
"omegav": 0.73,
"omegab": 0.046, # uses lower case omega without definition
"h": 0.72,
"n": 1.0,
},
)
[docs]
class Ishiyama(Warren):
"""Ishiyama mass function fit."""
_eq = r"A\left[\left(\frac{e}{\sigma}\right)^b + 1\right]\exp(\frac{d}{\sigma^2})"
references: ClassVar[tuple[str, ...]] = (
("Ishiyama, T., et al., 2015. PASJ 67, 61. https://doi.org/10.1093/pasj/psv021"),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Ishiyama", "Ishiyama", _eq, references)
_defaults: ClassVar[dict[str, float]] = {
"A": 0.193,
"b": 1.550,
"c": 1,
"d": 1.186,
"e": 2.184,
}
sim_definition = SimDetails(
L=[1120, 560, 280, 140, 70],
N=[8192**3, 4096**3, 2048**3, 2048**3, 2048**3],
halo_finder_type="FoF",
omegam=0.31,
sigma_8=0.83,
halo_overdensity=0.2,
halo_finder=None,
softening=[4.27, 4.27, 4.27, 2.14, 1.07],
transfer="CAMB",
z_start=None,
z_meas=0,
ICS=None,
nmin=40,
hmf_analysis_notes="No corrections applied.",
other_cosmo={
"omegav": 0.69,
"omegab": 0.048, # uses lower case omega without definition
"h": 0.68,
"n": 0.96,
},
)
@override
@property
def cutmask(self):
return np.logical_and(self.m > 1e8, self.m < 1e16)
[docs]
class Bocquet200mDMOnly(Warren):
"""Bocquet mass function fit for 200m definition with dark matter only."""
_eq = r"A\left[\left(\frac{e}{\sigma}\right)^b + 1\right]\exp(-\frac{d}{\sigma^2})"
references: ClassVar[tuple[str, ...]] = (
("Bocquet, S., et al., 2016. MNRAS 456, 2361. https://doi.org/10.1093/mnras/stv2657"),
)
__doc__ = _makedoc(BaseFittingFunction._pdocs, "Bocquet", "Bocquet", _eq, references)
_defaults: ClassVar[dict[str, Any]] = {
# Bocquet+16 v3 (published) Table 2. Paper (a, b, c) -> code (b, e, d).
"A": 0.175,
"b": 1.53,
"c": 1,
"d": 1.19,
"e": 2.55,
"A_z": -0.012,
"b_z": -0.04,
"d_z": -0.021,
"e_z": -0.194,
}
sim_definition = SimDetails(
L=[68.1, 181.8, 1274],
N=None,
halo_finder_type="SO",
omegam=0.272,
sigma_8=0.809,
halo_overdensity="200m",
halo_finder="Subfind",
softening=None,
transfer=None,
z_start=None,
z_meas=(0, 2),
ICS=None,
nmin=10000,
hmf_analysis_notes="Poisson bayesian likelihood and finite volume correction.",
other_cosmo={
"omegav": 0.69,
"omegab": 0.0456, # uses lower case omega without definition
"h": 0.704,
"n": 0.96,
},
)
[docs]
def get_params(self):
"""Get the redshift-dependent parameters."""
return (
self.params["A"] * (1 + self.z) ** self.params["A_z"],
self.params["b"] * (1 + self.z) ** self.params["b_z"],
self.params["d"] * (1 + self.z) ** self.params["d_z"],
self.params["e"] * (1 + self.z) ** self.params["e_z"],
)
[docs]
def convert_mass(self):
"""Function to compute mass in this definition compared to 200m.
This is an analytic approximation, not a full mass translation, and is calibrated
to the NFW profile with Duffy+08 concentration-mass relation. This ratio is
applied in :meth:`fsigma`.
"""
return 1
@override
@property
def fsigma(self):
A, b, d, e = self.get_params()
mass_conversion = self.convert_mass()
return A * ((e / self.sigma) ** b + 1) * np.exp(-d / self.sigma**2) * mass_conversion
[docs]
class Bocquet200mHydro(Bocquet200mDMOnly):
"""Bocquet mass function fit for 200m definition with hydrodynamics."""
__doc__ = _makedoc(
BaseFittingFunction._pdocs,
"Bocquet",
"Bocquet",
Bocquet200mDMOnly._eq,
Bocquet200mDMOnly.references,
)
_defaults: ClassVar[dict[str, Any]] = {
# Bocquet+16 v3 (published) Table 2. Paper (a, b, c) -> code (b, e, d).
"A": 0.228,
"b": 2.15,
"c": 1,
"d": 1.3,
"e": 1.69,
"A_z": 0.285,
"b_z": -0.058,
"d_z": -0.045,
"e_z": -0.366,
}
[docs]
class Bocquet200cDMOnly(Bocquet200mDMOnly):
"""Bocquet mass function fit for 200c definition with dark matter only."""
__doc__ = _makedoc(
BaseFittingFunction._pdocs,
"Bocquet",
"Bocquet",
Bocquet200mDMOnly._eq,
Bocquet200mDMOnly.references,
)
_defaults: ClassVar[dict[str, Any]] = {
# Bocquet+16 v3 (published) Table 2. Paper (a, b, c) -> code (b, e, d).
"A": 0.222,
"b": 1.71,
"c": 1,
"d": 1.46,
"e": 2.24,
"A_z": 0.269,
"b_z": 0.321,
"d_z": -0.153,
"e_z": -0.621,
}
sim_definition = copy(Bocquet200mDMOnly.sim_definition)
sim_definition.halo_overdensity = "200c"
[docs]
@override
def convert_mass(self):
g0 = 3.54e-2 + self.cosmo.Om0**0.09
g1 = 4.56e-2 + 2.68e-2 / self.cosmo.Om0
g2 = 0.721 + 3.5e-2 / self.cosmo.Om0
g3 = 0.628 + 0.164 / self.cosmo.Om0
d0 = -1.67e-2 + 2.18e-2 * self.cosmo.Om0
d1 = 6.52e-3 - 6.86e-3 * self.cosmo.Om0
g = g0 + g1 * np.exp(-(((g2 - self.z) / g3) ** 2))
d = d0 + d1 * self.z
# Eq. A2 takes ln(M/Msun); m is in Msun/h.
return g + d * np.log(self.m / self.cosmo.h)
[docs]
class Bocquet200cHydro(Bocquet200cDMOnly):
"""Bocquet mass function fit for 200c definition with hydrodynamics."""
__doc__ = _makedoc(
BaseFittingFunction._pdocs,
"Bocquet",
"Bocquet",
Bocquet200mDMOnly._eq,
Bocquet200mDMOnly.references,
)
_defaults: ClassVar[dict[str, Any]] = {
# Bocquet+16 v3 (published) Table 2. Paper (a, b, c) -> code (b, e, d).
"A": 0.202,
"b": 2.21,
"c": 1,
"d": 1.57,
"e": 2.0,
"A_z": 1.147,
"b_z": 0.375,
"d_z": -0.196,
"e_z": -1.074,
}
[docs]
class Bocquet500cDMOnly(Bocquet200cDMOnly):
"""Bocquet mass function fit for 500c definition with dark matter only."""
__doc__ = _makedoc(
BaseFittingFunction._pdocs,
"Bocquet",
"Bocquet",
Bocquet200mDMOnly._eq,
Bocquet200mDMOnly.references,
)
_defaults: ClassVar[dict[str, Any]] = {
# Bocquet+16 v3 (published) Table 2. Paper (a, b, c) -> code (b, e, d).
"A": 0.241,
"b": 2.18,
"c": 1,
"d": 2.02,
"e": 2.35,
"A_z": 0.37,
"b_z": 0.251,
"d_z": -0.31,
"e_z": -0.698,
}
sim_definition = copy(Bocquet200mDMOnly.sim_definition)
sim_definition.halo_overdensity = "500c"
[docs]
@override
def convert_mass(self):
alpha_0 = 0.880 + 0.329 * self.cosmo.Om0
alpha_1 = 1.0 + 4.31 * 1e-2 / self.cosmo.Om0
alpha_2 = -0.365 + 0.254 / self.cosmo.Om0
alpha = alpha_0 * (alpha_1 * self.z + alpha_2) / (self.z + alpha_2)
beta = -1.7e-2 + self.cosmo.Om0 * 3.74e-3
# Eq. 6 takes ln(M/Msun); m is in Msun/h.
return alpha + beta * np.log(self.m / self.cosmo.h)
[docs]
class Bocquet500cHydro(Bocquet500cDMOnly):
"""Bocquet mass function fit for 500c definition with hydrodynamics."""
__doc__ = _makedoc(
BaseFittingFunction._pdocs,
"Bocquet",
"Bocquet",
Bocquet200mDMOnly._eq,
Bocquet200mDMOnly.references,
)
_defaults: ClassVar[dict[str, Any]] = {
# Bocquet+16 v3 (published) Table 2. Paper (a, b, c) -> code (b, e, d).
"A": 0.18,
"b": 2.29,
"c": 1,
"d": 1.97,
"e": 2.44,
"A_z": 1.088,
"b_z": 0.15,
"d_z": -0.322,
"e_z": -1.008,
}
[docs]
class Yung24(BaseFittingFunction):
r"""Yung+24 GUREFT halo mass function fit, calibrated for :math:`z=6{-}19`.
The functional form follows Tinker08 / Rodriguez-Puebla16, with the four
parameters :math:`A, a, b, c` each expressed as a quadratic polynomial in
redshift :math:`\chi(z) = \chi_0 + \chi_1 z + \chi_2 z^2`. Best-fit
coefficients come from the combined gureft + MultiDark fit in Yung+24
Appendix A.
Parameters
----------
units : {"h", "physical"}, optional
Selects the parameter table. ``"h"`` (default) uses Table A1, calibrated
for masses in :math:`M_\odot/h` and volumes in
:math:`{\rm Mpc}^{-3} h^3 {\rm dex}^{-1}` -- matching ``hmf``'s
standard internal convention. ``"physical"`` uses Table A2, calibrated
for masses in :math:`M_\odot` and volumes in
:math:`{\rm Mpc}^{-3} {\rm dex}^{-1}`.
Raises
------
ValueError
If ``z`` is outside the calibration range :math:`[6, 19]`.
Notes
-----
The Yung+24 [1]_ form is:
.. math::
f_{\rm Yung24}(\sigma) = A(z)\left[\left(\frac{\sigma}{b(z)}\right)^{-a(z)}
+ 1\right]\exp\left(-\frac{c(z)}{\sigma^2}\right)
The ``*_0``, ``*_1`` and ``*_2`` parameters below are :math:`\chi_0, \chi_1, \chi_2`
for :math:`\chi \in \{A, a, b, c\}`. The defaults are those for
``units="h"``; with ``units="physical"`` any coefficients not passed explicitly
are taken from Table A2 instead.
References
----------
.. [1] Yung, L.Y.A., Somerville, R.S., Nguyen, T., Behroozi, P., Modi, C.,
Gardner, J.P., 2024. MNRAS 530, 4868. arXiv:2309.14408
"""
req_z = True
_eq = r"A(z)\left[(\sigma/b(z))^{-a(z)} + 1\right]\exp(-c(z)/\sigma^2)"
references: ClassVar[tuple[str, ...]] = (
(
"Yung, L.Y.A., Somerville, R.S., Nguyen, T., Behroozi, P., Modi, C., Gardner, "
"J.P., 2024. MNRAS 530, 4868. "
"arXiv:2309.14408"
),
)
sim_definition = SimDetails(
L=[5, 15, 35, 90, 250, 400],
N=[1024**3] * 6,
halo_finder_type="SO",
halo_overdensity="vir",
halo_finder="Rockstar",
omegam=0.307,
sigma_8=0.829,
z_meas=(6.0, 19.0),
ICS="2LPT",
hmf_analysis_notes="Combined gureft + MultiDark fit (Yung+24 Appendix A).",
)
_z_min = 6.0
_z_max = 19.0
_params_h: ClassVar[dict[str, float]] = {
"A_0": 0.11416632,
"A_1": -0.01486746,
"A_2": 0.00137191,
"a_0": 1.05274399,
"a_1": 0.02803087,
"a_2": -0.00306126,
"b_0": 8.62813020,
"b_1": 0.00384969,
"b_2": -0.02349983,
"c_0": 1.13138924,
"c_1": 0.01713172,
"c_2": -0.00113630,
}
_params_phys: ClassVar[dict[str, float]] = {
"A_0": 0.13765772,
"A_1": -0.01003821,
"A_2": 0.00102964,
"a_0": 1.06641384,
"a_1": 0.02475576,
"a_2": -0.00283342,
"b_0": 4.86693806,
"b_1": 0.09212356,
"b_2": -0.01426283,
"c_0": 1.19837952,
"c_1": -0.00142967,
"c_2": -0.00033074,
}
_defaults: ClassVar[dict[str, Any]] = {**_params_h, "units": "h"}
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
if self.z < self._z_min or self.z > self._z_max:
raise ValueError(
f"Yung24 fit is only valid for z in [{self._z_min}, {self._z_max}]; got z={self.z}."
)
units = self.params["units"]
if units == "physical":
user_overrides = set(kwargs) - {"units"}
for k, v in self._params_phys.items():
if k not in user_overrides:
self.params[k] = v
elif units != "h":
raise ValueError(f"units must be 'h' or 'physical', got {units!r}")
@override
@property
def fsigma(self):
z = self.z
A = self.params["A_0"] + self.params["A_1"] * z + self.params["A_2"] * z**2
a = self.params["a_0"] + self.params["a_1"] * z + self.params["a_2"] * z**2
b = self.params["b_0"] + self.params["b_1"] * z + self.params["b_2"] * z**2
c = self.params["c_0"] + self.params["c_1"] * z + self.params["c_2"] * z**2
return A * ((self.sigma / b) ** (-a) + 1) * np.exp(-c / self.sigma**2)
@override
@property
def cutmask(self):
if self.m is None:
return np.ones(len(self.nu2), dtype=bool)
log10m = np.log10(self.m)
lo = 6.0 if self.params["units"] == "h" else 5.0
return (log10m >= lo) & (log10m <= 13.0)