"""Models defining various definitions of halo mass.
This is primarily inspired by Benedikt Diemer's COLOSSUS code:
https://bdiemer.bitbucket.io/colossus/halo_mass_defs.html.
Definitions include both spherical overdensity (with various overdensities) and
friends-of-friends (with various linking lengths). The main class is
:class:`BaseMassDefinition`, which defines the interface for all mass definitions.
"""
import warnings
from typing import ClassVar, override
import astropy.units as u
import numpy as np
import scipy as sp
from astropy.cosmology import Planck15
from .._internals import _framework
from ..cosmology import Cosmology
__all__ = [
"FOF",
"BaseMassDefinition",
"MassDefinition",
"OptimizationError",
"SOCritical",
"SOMean",
"SOVirial",
"SphericalOverdensity",
]
[docs]
@_framework.pluggable
class BaseMassDefinition(_framework.Component):
"""A base class for a Mass Definition."""
[docs]
@staticmethod
def critical_density(z=0, cosmo=Planck15):
"""Get the critical density of the Universe at redshift z, [h^2 Msun/Mpc^3]."""
return (cosmo.critical_density(z) / cosmo.h**2).to(u.Msun / u.Mpc**3).value
[docs]
@classmethod
def mean_density(cls, z=0, cosmo=Planck15):
"""Get the mean density of the Universe at redshift z, [h^2 Msun / Mpc^3]."""
return cosmo.Om(z) * cls.critical_density(z, cosmo)
[docs]
def halo_density(self, z=0, cosmo=Planck15):
r"""
The density of haloes under this definition.
May not exist in some definitions. Units are :math:`M_\odot h^2/{\rm Mpc}^3`.
"""
raise AttributeError("halo_density does not exist for this Mass Definition")
@property
def colossus_name(self):
"""The name of the mass definition in Colossus format, if applicable."""
return None
[docs]
def halo_overdensity_mean(self, z=0, cosmo=Planck15):
"""Compute the halo overdensity with respect to the mean density."""
return self.halo_density(z, cosmo) / self.mean_density(z, cosmo)
[docs]
def halo_overdensity_crit(self, z=0, cosmo=Planck15):
"""Compute the halo overdensity with respect to the critical density."""
return self.halo_density(z, cosmo) / self.critical_density(z, cosmo)
[docs]
def m_to_r(self, m, z=0, cosmo=Planck15):
r"""
Return the radius corresponding to m for this mass definition.
Parameters
----------
m : float or array_like
The mass to convert to radius. Should be in the same units (modulo volume)
as :meth:`halo_density`.
Notes
-----
Computed as :math:`\left(\frac{3m}{4\pi \rho_{\rm halo}\right)`.
"""
try:
return (3 * m / (4 * np.pi * self.halo_density(z, cosmo))) ** (1.0 / 3.0)
except AttributeError as e:
raise AttributeError(f"{self.__class__.__name__} cannot convert mass to radius.") from e
[docs]
def r_to_m(self, r, z=0, cosmo=Planck15):
r"""
Return the mass corresponding to r for this mass definition.
Parameters
----------
r : float or array_like
The radius to convert to mass. Units should be compatible with
:meth:`halo_density`.
Notes
-----
Computed as :math:`\frac{4\pi r^3}{3} \rho_{\rm halo}`.
"""
try:
return 4 * np.pi * r**3 * self.halo_density(z, cosmo) / 3
except AttributeError as e:
raise AttributeError(f"{self.__class__.__name__} cannot convert radius to mass.") from e
def _duffy_concentration(self, m, z=0):
a, b, c, ms = 6.71, -0.091, 0.44, 2e12
return a / (1 + z) ** c * (m / ms) ** b
[docs]
def change_definition(self, m: np.ndarray, mdef, profile=None, c=None, z=0, cosmo=Planck15):
r"""
Change the spherical overdensity mass definition.
This requires using a profile, for which the `halomod` package must be used.
Parameters
----------
m : float or array_like
The halo mass to be changed, in :math:`M_\odot/h`. Must be
broadcastable with `c`, if provided.
mdef : :class:`MassDefinition` subclass instance
The mass definition to which to change.
profile : :class:`halomod.profiles.Profile` instance, optional
An instantiated profile object from which to calculate the expected
definition change. If not provided, a mocked NFW profile is used.
c : float or array_like, optional
The concentration(s) of the halos given. If not given, the concentrations
will be automatically calculated using the profile object.
Returns
-------
m_f : float or array_like
The masses of the halos in the new definition.
r_f : float or array_like
The radii of the halos in the new definition.
c_f : float or array_like
The concentrations of the halos in the new definition.
"""
if c is not None and not np.isscalar(c) and not np.isscalar(m) and len(m) != len(c):
raise ValueError("If both m and c are arrays, they must be of the same length")
if c is not None and np.isscalar(c) and not np.isscalar(m):
c = np.ones_like(m) * c
if c is not None and np.isscalar(m) and not np.isscalar(c):
m = np.ones_like(m) * m
if c is not None:
c = np.atleast_1d(c)
m = np.atleast_1d(m)
if profile is None:
try:
from halomod.concentration import Duffy08
from halomod.profiles import NFW
profile = NFW(
cm_relation=Duffy08(cosmo=Cosmology(cosmo)),
mdef=self,
z=z,
cosmo=cosmo,
)
except ImportError as e:
raise ImportError(
"Cannot change mass definitions without halomod installed!"
) from e
if profile.z != z:
warnings.warn(
f"Redshift of given profile ({profile.z})does not match redshift "
f"passed to change_definition(). Using the redshift directly passed.",
stacklevel=2,
)
profile.z = z
if c is None:
c = profile.cm_relation(m)
rs = self.m_to_r(m, z, cosmo) / c
rhos = profile._rho_s(c)
if not hasattr(rhos, "__len__"):
rhos = [rhos]
if not hasattr(c, "__len__"):
c = [c]
c_new = np.array(
[
_find_new_concentration(rho, mdef.halo_density(z, cosmo), profile._h, cc)
for rho, cc in zip(rhos, c, strict=True)
]
)
if len(c_new) == 1:
c_new = c_new[0]
r_new = c_new * rs
if len(r_new) == 1:
r_new = r_new[0]
return mdef.r_to_m(r_new, z, cosmo), r_new, c_new
def __eq__(self, other):
"""Test equality with another object."""
return self.__class__.__name__ == other.__class__.__name__ and self.params == other.params
# For backwards compatibility, alias MassDefinition to BaseMassDefinition.
MassDefinition = BaseMassDefinition
[docs]
class SphericalOverdensity(BaseMassDefinition):
"""An abstract base class for all spherical overdensity mass definitions."""
def __str__(self):
"""Describe the overdensity in standard notation."""
return f"{self.__class__.__name__}({self.params['overdensity']})"
[docs]
class SOGeneric(SphericalOverdensity):
"""A generic SO definition which can claim equality with any SO."""
[docs]
def __init__(self, preferred: SphericalOverdensity | None = None, **kwargs):
super().__init__(**kwargs)
self.preferred = preferred
def __eq__(self, other):
"""Test equality with another object."""
return isinstance(other, SphericalOverdensity)
def __str__(self):
return "SOGeneric"
[docs]
class SOMean(SphericalOverdensity):
"""A mass definition based on spherical overdensity wrt mean background density."""
_defaults: ClassVar[dict[str, float]] = {"overdensity": 200}
[docs]
def halo_density(self, z=0, cosmo=Planck15):
"""The density of haloes under this definition."""
return self.params["overdensity"] * self.mean_density(z, cosmo)
@override
@property
def colossus_name(self):
return f"{int(self.params['overdensity'])}m"
[docs]
class SOCritical(SphericalOverdensity):
"""A mass definition based on spherical overdensity wrt critical density."""
_defaults: ClassVar[dict[str, float]] = {"overdensity": 200}
[docs]
def halo_density(self, z=0, cosmo=Planck15):
"""The density of haloes under this definition."""
return self.params["overdensity"] * self.critical_density(z, cosmo)
@override
@property
def colossus_name(self):
return f"{int(self.params['overdensity'])}c"
[docs]
class SOVirial(SphericalOverdensity):
"""A mass definition based on spherical overdensity.
Density threshold isgiven by Bryan and Norman (1998).
"""
[docs]
def halo_density(self, z=0, cosmo=Planck15):
"""The density of haloes under this definition."""
x = cosmo.Om(z) - 1
overdensity = 18 * np.pi**2 + 82 * x - 39 * x**2
return overdensity * self.mean_density(z, cosmo) / cosmo.Om(z)
@override
@property
def colossus_name(self):
return "vir"
def __str__(self):
"""Describe the halo definition in standard notation."""
return "SOVirial"
[docs]
class FOF(BaseMassDefinition):
"""A mass definition based on FroF networks with given linking length."""
_defaults: ClassVar[dict[str, float]] = {"linking_length": 0.2}
[docs]
def halo_density(self, z=0, cosmo=Planck15):
r"""
The density of halos under this mass definition.
Note that for FoF halos, this is very approximate. We follow [1]_ and define
:math:`rho_{FOF} = 9/(2\pi b^3) \rho_m`, with *b* the linking length. This
assumes all groups are spherical and singular isothermal spheres.
References
----------
.. [1] White, Martin, Lars Hernquist, and Volker Springel. “The Halo Model and
Numerical Simulations.” The Astrophysical Journal 550, no. 2 (April 2001):
L129-32. https://doi.org/10.1086/319644.
"""
overdensity = 9 / (2 * np.pi * self.params["linking_length"] ** 3)
return overdensity * self.mean_density(z, cosmo)
@override
@property
def colossus_name(self):
return "fof"
def __str__(self):
"""Describe the halo definition in standard notation."""
return f"FoF(l={self.params['linking_length']})"
[docs]
def from_colossus_name(name):
if name == "vir":
return SOVirial()
if name.endswith("c"):
return SOCritical(overdensity=int(name[:-1]))
if name.endswith("m"):
return SOMean(overdensity=int(name[:-1]))
if name == "fof":
return FOF()
raise ValueError(f"name '{name}' is an unknown mass definition to colossus.")
def _find_new_concentration(rho_s, halo_density, h=None, x_guess=5.0):
r"""
Find :math:`x=r/r_{\\rm s}` where the enclosed density has a particular value.
.. note :: This is almost exactly the same code as profileNFW.xDelta from COLOSSUS.
It may one day be changed to literally just call that function. For now
it just sits here to be called whenever halomod is not installed
Parameters
----------
rho_s: float
The central density in physical units :math:`M_{\odot} h^2 / {\\rm Mpc}^3`.
halo_density: float
The desired enclosed density threshold in physical units
:math:`M_{\odot} h^2 / {\\rm Mpc}^3`.
h : callable
Return the enclosed density as function of r/r_s.
Returns
-------
x: float
The radius in units of the scale radius, :math:`x=r/r_{\\rm s}`, where the
enclosed density reaches ``density_threshold``.
"""
# A priori, we have no idea at what radius the result will come out, but we need to
# provide lower and upper limits for the root finder. To balance stability and
# performance, we do so iteratively: if there is no result within relatively
# aggressive limits, we try again with more conservative limits.
target = halo_density / rho_s
x = None
i = 0
XDELTA_GUESS_FACTORS = [5.0, 10.0, 20.0, 100.0, 10000.0]
if h is None:
def h(x):
return np.log(1.0 + x) - x / (1.0 + x)
def fnc(x):
return h(x) * 3.0 / x**3 - target
while x is None and i < len(XDELTA_GUESS_FACTORS):
try:
xmin = x_guess / XDELTA_GUESS_FACTORS[i]
xmax = x_guess * XDELTA_GUESS_FACTORS[i]
x = sp.optimize.brentq(fnc, xmin, xmax)
except Exception as e:
warnings.warn(f"raised following error: {e}", stacklevel=2)
i += 1
if x is None:
raise OptimizationError(
"Could not determine x where the density threshold "
f"{halo_density / rho_s:.2f} is satisfied. Largest x-range tried was "
f"x={xmin} -- {xmax} which brackets {fnc(xmin)} -- {fnc(xmax)}"
)
return x
class OptimizationError(Exception):
"""Exception class related to failed optimization."""