"""
This module contains occupation components used by the ZuMandelbaum15 composite model.
"""
import numpy as np
from scipy.special import erf
from .occupation_model_template import OccupationComponent
from .. import model_defaults
from ..smhm_models import ZuMandelbaum15SmHm
__all__ = ("ZuMandelbaum15Cens", "ZuMandelbaum15Sats")
[docs]
class ZuMandelbaum15Cens(OccupationComponent):
"""HOD-style model for any central galaxy occupation that derives from
a stellar-to-halo-mass relation.
.. note::
The `~halotools.empirical_models.ZuMandelbaum15Cens` model is part of
the ``zu_mandelbaum15`` prebuilt composite HOD-style model.
For a tutorial on the ``zu_mandelbaum15``
composite model, see :ref:`zu_mandelbaum15_composite_model`.
"""
def __init__(
self,
threshold=model_defaults.default_stellar_mass_threshold,
prim_haloprop_key="halo_m200m",
**kwargs
):
"""
Parameters
----------
threshold : float, optional
Stellar mass threshold of the mock galaxy sample in h=1 solar mass units.
Default value is specified in the `~halotools.empirical_models.model_defaults` module.
prim_haloprop_key : string, optional
String giving the column name of the primary halo property governing
the occupation statistics of gal_type galaxies.
Default value is ``halo_m200m``, as in Zu & Mandelbaum (2015)
redshift : float, optional
Redshift of the stellar-to-halo-mass relation. Default is z=0.
Examples
--------
>>> cen_model = ZuMandelbaum15Cens()
>>> cen_model = ZuMandelbaum15Cens(threshold=11.25)
>>> cen_model = ZuMandelbaum15Cens(prim_haloprop_key='halo_m200b')
Notes
-----
Note also that the best-fit parameters of this model are based on the
``halo_m200m`` halo mass definition.
Using alternative choices of mass definition will require altering the
model parameters in order to mock up the same model published in Zu & Mandelbaum 2015.
The `Colossus python package <https://bitbucket.org/bdiemer/colossus/>`_
written by Benedikt Diemer can be used to
convert between different halo mass definitions. This may be useful if you wish to use an
existing halo catalog for which the halo mass definition you need is unavailable.
"""
upper_occupation_bound = 1.0
# Call the super class constructor, which binds all the
# arguments to the instance.
OccupationComponent.__init__(
self,
gal_type="centrals",
threshold=threshold,
upper_occupation_bound=upper_occupation_bound,
prim_haloprop_key=prim_haloprop_key,
)
self.smhm_model = ZuMandelbaum15SmHm(prim_haloprop_key=prim_haloprop_key)
for key, value in self.smhm_model.param_dict.items():
self.param_dict[key] = value
self._methods_to_inherit = [
"mc_occupation",
"mean_occupation",
"mean_stellar_mass",
"mean_halo_mass",
]
self.publications = ["arXiv:1103.2077", "arXiv:1104.0928", "1505.02781"]
self.publications.extend(self.smhm_model.publications)
self.publications = list(set(self.publications))
[docs]
def get_published_parameters(self):
""""""
return ZuMandelbaum15SmHm.get_published_parameters(self)
[docs]
def mean_occupation(self, **kwargs):
"""Expected number of central galaxies in a halo.
Parameters
----------
prim_haloprop : array, optional
Array of mass-like variable upon which occupation statistics are based.
If ``prim_haloprop`` is not passed, then ``table`` keyword argument must be passed.
table : object, optional
Data table storing halo catalog.
If ``table`` is not passed, then ``prim_haloprop`` keyword argument must be passed.
Returns
-------
mean_ncen : array
Mean number of central galaxies in the halo of the input mass.
Examples
--------
>>> cen_model = ZuMandelbaum15Cens(threshold=10.75)
>>> halo_masses = np.logspace(11, 15, 25)
>>> mean_ncen = cen_model.mean_occupation(prim_haloprop=halo_masses)
"""
for key, value in self.param_dict.items():
if key in list(self.smhm_model.param_dict.keys()):
self.smhm_model.param_dict[key] = value
if "table" in list(kwargs.keys()):
halo_mass = np.atleast_1d(kwargs["table"][self.prim_haloprop_key])
elif "prim_haloprop" in list(kwargs.keys()):
halo_mass = np.atleast_1d(kwargs["prim_haloprop"])
else:
raise KeyError(
"Must pass one of the following keyword arguments "
"to mean_occupation:\n``table`` or ``prim_haloprop``"
)
sigma = self.smhm_model.scatter_ln_mstar(halo_mass)
mean = self.smhm_model.mean_stellar_mass(prim_haloprop=halo_mass)
erfarg = (np.log(10**self.threshold) - np.log(mean)) / (sigma * np.sqrt(2))
return 0.5 * (1 - erf(erfarg))
[docs]
def mean_stellar_mass(self, **kwargs):
"""Return the stellar mass of a central galaxy as a function
of the input table.
Parameters
----------
prim_haloprop : array, optional
Array of mass-like variable upon which occupation statistics are based.
If ``prim_haloprop`` is not passed, then ``table`` keyword argument must be passed.
table : object, optional
Data table storing halo catalog.
If ``table`` is not passed, then ``prim_haloprop`` keyword argument must be passed.
Returns
-------
mstar : array_like
Array containing stellar masses living in the input table.
Examples
--------
>>> cen_model = ZuMandelbaum15Cens(threshold=10.75)
>>> halo_masses = np.logspace(11, 15, 25)
>>> mstar = cen_model.mean_stellar_mass(prim_haloprop=halo_masses)
"""
for key, value in self.param_dict.items():
if key in self.smhm_model.param_dict:
self.smhm_model.param_dict[key] = value
return self.smhm_model.mean_stellar_mass(**kwargs)
[docs]
def mean_halo_mass(self, stellar_mass):
"""Return the halo mass of a central galaxy as a function
of the input stellar mass.
Parameters
----------
stellar_mass : array
Array of stellar masses in h=1 solar mass units.
Returns
-------
halo_mass : array_like
Array containing halo mass in h=1 solar mass units.
Examples
--------
>>> cen_model = ZuMandelbaum15Cens(threshold=10.75)
>>> stellar_mass = np.logspace(9, 12, 25)
>>> halo_mass = cen_model.mean_halo_mass(stellar_mass)
"""
for key, value in self.param_dict.items():
if key in self.smhm_model.param_dict:
self.smhm_model.param_dict[key] = value
return self.smhm_model.mean_halo_mass(stellar_mass)
[docs]
class ZuMandelbaum15Sats(OccupationComponent):
r"""HOD-style model for a satellite galaxy occupation
based on Zu & Mandelbaum 2015.
.. note::
The `~halotools.empirical_models.ZuMandelbaum15Sats` model is part of
the ``zu_mandelbaum15`` prebuilt composite HOD-style model.
For a tutorial on the ``zu_mandelbaum15``
composite model, see :ref:`zu_mandelbaum15_composite_model`.
"""
def __init__(
self,
threshold=model_defaults.default_stellar_mass_threshold,
prim_haloprop_key="halo_m200m",
**kwargs
):
r"""
Parameters
----------
threshold : float, optional
Stellar mass threshold of the mock galaxy sample in h=1 solar mass units.
Default value is specified in the `~halotools.empirical_models.model_defaults` module.
prim_haloprop_key : string, optional
String giving the column name of the primary halo property governing
the occupation statistics of gal_type galaxies.
Default value is specified in the `~halotools.empirical_models.model_defaults` module.
Examples
--------
>>> sat_model = ZuMandelbaum15Sats()
>>> sat_model = ZuMandelbaum15Sats(threshold=11)
>>> sat_model = ZuMandelbaum15Sats(prim_haloprop_key='halo_mvir')
Notes
-----
Note also that the best-fit parameters of this model are based on the
``halo_m200m`` halo mass definition.
Using alternative choices of mass definition will require altering the
model parameters in order to mock up the same model published in Zu & Mandelbaum 2015.
The `Colossus python package <https://bitbucket.org/bdiemer/colossus/>`_
written by Benedikt Diemer can be used to
convert between different halo mass definitions. This may be useful if you wish to use an
existing halo catalog for which the halo mass definition you need is unavailable.
"""
self.central_occupation_model = ZuMandelbaum15Cens(
prim_haloprop_key=prim_haloprop_key, threshold=threshold
)
OccupationComponent.__init__(
self,
gal_type="satellites",
threshold=threshold,
upper_occupation_bound=float("inf"),
prim_haloprop_key=prim_haloprop_key,
)
self._initialize_param_dict()
self.param_dict.update(self.central_occupation_model.param_dict)
self.publications = self.central_occupation_model.publications
[docs]
def mean_occupation(self, **kwargs):
"""Expected number of satellite galaxies in a halo of mass halo_mass.
Parameters
----------
prim_haloprop : array, optional
array of masses of table in the catalog
table : object, optional
Data table storing halo catalog.
Returns
-------
mean_nsat : array
Mean number of satellite galaxies in the halo of the input mass.
Examples
--------
>>> sat_model = ZuMandelbaum15Sats()
>>> halo_masses = np.logspace(11, 15, 25)
>>> mean_nsat = sat_model.mean_occupation(prim_haloprop=halo_masses)
"""
# Retrieve the array storing the mass-like variable
if "table" in list(kwargs.keys()):
halo_mass = kwargs["table"][self.prim_haloprop_key]
elif "prim_haloprop" in list(kwargs.keys()):
halo_mass = np.atleast_1d(kwargs["prim_haloprop"])
else:
raise KeyError(
"Must pass one of the following keyword arguments "
"to mean_occupation:\n``table`` or ``prim_haloprop``"
)
self._update_satellite_params()
mean_ncen = self.central_occupation_model.mean_occupation(**kwargs)
mean_nsat = (
mean_ncen
* np.exp(-self._mcut / halo_mass)
* (halo_mass / self._msat) ** self.param_dict["alphasat"]
)
return mean_nsat
def _initialize_param_dict(self):
"""Set the initial values of ``self.param_dict`` according to
the SIG_MOD1 values of Table 5 of arXiv:1104.0928 for the
lowest redshift bin.
"""
self.param_dict["alphasat"] = 1.0
self.param_dict["bsat"] = 8.98
self.param_dict["betasat"] = 0.9
self.param_dict["bcut"] = 0.86
self.param_dict["betacut"] = 0.41
self.param_dict.update(self.central_occupation_model.param_dict)
self._update_satellite_params()
def _update_satellite_params(self):
"""Private method to update the model parameters."""
for key, value in self.param_dict.items():
if key in self.central_occupation_model.param_dict:
self.central_occupation_model.param_dict[key] = value
knee_threshold = self.central_occupation_model.mean_halo_mass(
10**self.threshold
)
knee_mass = 1.0e12
self._msat = (
knee_mass
* self.param_dict["bsat"]
* (knee_threshold / knee_mass) ** self.param_dict["betasat"]
)
self._mcut = (
knee_mass
* self.param_dict["bcut"]
* (knee_threshold / knee_mass) ** self.param_dict["betacut"]
)