A Showcase of Components in halomod

In this demo, we will showcase each and every one of the different models and components in halomod, to give a taste of its capabilities. Note that the plots generated in this notebook appear in THIS PAPER.

import hmf
import numpy as np

import halomod
from halomod.bias import make_colossus_bias
from halomod.concentration import make_colossus_cm

print(f"Using halomod v{halomod.__version__} and hmf v{hmf.__version__}")

Using halomod v2.0.2.dev51+geee0902 and hmf v3.5.0

Setup Matplotlib

%matplotlib inline

import matplotlib.pyplot as plt
from matplotlib import rc, rcParams

Since these plots will be included in a paper, we define some preset configurations for different kinds of plots, to achieve some consistency.

paper_font_size = 10  # elsarticle


def single_column(height_ratio=1 / 1.5):
    linewidth = 3.49  # Get this in your tex using \printinunitsof{in}\prntlen{\linewidth}
    paper_settings(linewidth, height_ratio)


def two_column(height_ratio=1 / 1.5):
    textwidth = 7.22  # Get this in your tex using \printinunitsof{in}\prntlen{\textwidth}
    paper_settings(textwidth, height_ratio)


def paper_settings(width, height_ratio=1 / 1.5):
    plt.style.use("dark_background")

    rc("figure", figsize=(width, width * height_ratio), dpi=200)
    rc(("xtick", "ytick"), labelsize=paper_font_size - 2, direction="in")
    rc("legend", fontsize=paper_font_size - 1, frameon=False, columnspacing=1, handletextpad=0.4)
    rc("axes", labelsize=paper_font_size - 1)

    rc("savefig", dpi=600, bbox="tight")

Setup default model

Some of the components presented here can be calculated using lower-level frameworks, but all of them can be calculated from the TracerHaloModel, so we use this.

from halomod import TracerHaloModel

base_model = TracerHaloModel(
    dr_table=0.03,
    dlnk=0.05,
    dlog10m=0.05,
    transfer_params={"kmax": 1e2, "extrapolate_with_eh": True},
)

# Pre-compute many of the relevant quantities
base_model.power_auto_tracer;

Since most of the plots we make will be comparing all possible models of a particular component, we define a simple function to find all the models for a particular component:

For bias models:

import importlib


def get_models(module, component, ignore=()):
    cmp = getattr(importlib.import_module(module), component)
    models = cmp.get_models()
    return [m for name, m in models.items() if name not in ignore]

get_models("halomod.bias", "Bias")

[halomod.bias.UnityBias,
 halomod.bias.Mo96,
 halomod.bias.Jing98,
 halomod.bias.ST99,
 halomod.bias.SMT01,
 halomod.bias.Seljak04,
 halomod.bias.Seljak04Cosmo,
 halomod.bias.Tinker05,
 halomod.bias.Mandelbaum05,
 halomod.bias.Pillepich10,
 halomod.bias.Manera10,
 halomod.bias.Tinker10,
 halomod.bias.Tinker10PBSplit]

Now, define a function that will do most of the plotting for us. It will grab all the models defined for a particular component, and plot a resulting quantity for all of them, with a lower panel showing the ratio to the fiducial model.

import random


def model_comparison_plot(
    module,
    component,
    attr_name,
    x,
    y,
    ignore=(),
    height_ratio=1 / 1.5,
    ylab=None,
    y_unit=None,
    yratio_lab=None,
    ylog=True,
    comparison_limits=(0.8, 1.2),
    xlim=None,
    ylim=None,
    must_include=(),
    ncol_legend=1,
    legend_size=None,
    max_models=10,
    extra_models=None,
    extra_params=None,
):
    extra_params = extra_params or {}
    extra_models = extra_models or []
    # Get all the models for this component
    models = get_models(module, component, ignore)
    if extra_models:
        models += extra_models

    # Limit them to seven models (randomly)
    if len(models) > max_models:
        tmp_models = [m for m in models if m.__name__ in must_include or m in extra_models]
        for m in tmp_models:
            models.remove(m)

        models = random.sample(models, max_models - len(tmp_models))
        models.extend(tmp_models)

    # Most model comparison plots will be single column in the paper.
    single_column(height_ratio)

    subplot_kw = {"xscale": "log"}
    if xlim:
        subplot_kw["xlim"] = xlim

    # Create a two-panel plot with the lower panel slightly smaller.
    _fig, ax = plt.subplots(
        2, 1, sharex=True, gridspec_kw={"hspace": 0, "height_ratios": [3, 2]}, subplot_kw=subplot_kw
    )

    # Plot the different models.
    model = base_model.clone()  # Make a new updateable model.
    for i, cmp_model in enumerate(models):
        if cmp_model.__name__ != getattr(base_model, attr_name).__name__:
            model.update(
                **{attr_name.replace("model", "params"): extra_params.get(cmp_model.__name__, {})}
            )

        model.update(**{attr_name: cmp_model})
        ls = ":" if cmp_model in extra_models else ["-", "--", ":"][i // 10]

        xx = getattr(model, x)
        if x == "nu":
            xx = np.sqrt(xx)
        ax[0].plot(xx, getattr(model, y), label=cmp_model.__name__, lw=1, ls=ls)
        ax[1].plot(xx, getattr(model, y) / getattr(base_model, y), ls=ls)

    if ylog:
        ax[0].set_yscale("log")

    ax[0].legend(ncol=ncol_legend, fontsize=legend_size or rcParams["legend.fontsize"])

    ylab = ylab or y

    # Set axis labels and limits.
    ylabel = f"{ylab} [{y_unit}]" if y_unit else ylab

    ax[0].set_ylabel(ylabel)
    if x == "k":
        ax[1].set_xlabel("k [$h$ / Mpc]")
    elif x == "m":
        ax[1].set_xlabel(r"m [$h^{-1} M_\odot$]")
    elif x == "nu":
        ax[1].set_xlabel(r"Peak-height, $\nu$")

    if yratio_lab is None:
        yratio_lab = ylab

    ax[1].set_ylabel(
        r"{} / ${}_{{\rm {}}}$".format(
            yratio_lab, yratio_lab.replace("$", ""), getattr(base_model, attr_name).__name__
        )
    )
    ax[1].set_ylim(comparison_limits)

    if ylim:
        ax[0].set_ylim(ylim)

    # Save the figure
    plt.savefig(f"plots/{attr_name}s.pdf")

Transfer Function

First, transfer functions. We omit the models that require user input (eg. FromFile).

The transfer function determines the linear power spectrum, which in turn affects the halo mass function (more power on larger scales means more high-mass halos).

model_comparison_plot(
    module="hmf.density_field.transfer_models",
    component="TransferComponent",
    ignore=("FromFile", "FromArray", "EH"),
    attr_name="transfer_model",
    x="k",
    y="power",
    ylab=r"$P_{\rm lin}(k)$",
    y_unit=r"$h^{-3}{\rm Mpc}^3$",
    xlim=(1e-4, 1e2),
    yratio_lab=r"$P^{\rm lin}$",
    height_ratio=1 / 1.2,
    ncol_legend=2,
    ylim=(1e-4, 1e5),
    comparison_limits=(0.8, 1.25),
)

/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/hmf/density_field/transfer_models.py:233: UserWarning: 'extrapolate_with_eh' was not set. Defaulting to True, which is different behaviour than versions <=3.4.4. This warning may be removed in v4.0. Silence it by setting extrapolate_with_eh explicitly.
  warnings.warn(

First notice that the fiducial model is CAMB (this is only true if camb is installed, otherwise it is EH_BAO). Since CAMB has BAOs, all the models without BAOs show up in the ratio as having wiggles. Interestingly, CAMB has much higher power on small scales than the other models. CAMB is much slower to run than the other analytic models.

The normalization difference on large scales is due to the difference in the integral of the power spectrum, which is normalized to be the same value via \(\sigma_8\).

Mass Definitions

Here we plot the change in mass and concentration when moving from one mass definition to another.

from hmf.halos.mass_definitions import FOF, SOCritical, SOMean, SOVirial

single_column(height_ratio=1 / 1.5)


fig, ax = plt.subplots(
    2,
    1,
    sharex=True,
    gridspec_kw={"hspace": 0, "height_ratios": [3, 2]},
    subplot_kw={"xscale": "log"},
)

for mdef in [SOMean(overdensity=300), SOCritical(overdensity=300), SOVirial(), FOF()]:
    m, r, c = base_model.mdef.change_definition(
        base_model.m, mdef=mdef, profile=base_model.halo_profile, cosmo=base_model.cosmo
    )

    ax[0].plot(base_model.m, m / base_model.m, label=mdef.colossus_name)
    ax[1].plot(base_model.m, c / base_model.cmz_relation)

ax[0].grid(True)
ax[1].grid(True)

ax[0].legend(ncol=2, frameon=True)

ax[0].set_ylabel(r"$m_{\rm new} / m_{\rm old}$")
ax[1].set_ylabel(r"$c_{\rm new} / c_{\rm old}$")
ax[1].set_xlabel(r"m, [$h^{-1}M_{\odot}$]")

plt.savefig("plots/mass_conversion.pdf")

Notice that halos with a higher overdensity threshold move to lower masses and lower concentration.

Window Functions / Filters

The filter model, \(W(x)\), is convolved with the power spectrum to produce the mass variance (which in turn is used to define the halo mass function). Optimally, the shape of the filter model is matched to the shape of primordial patches that become halos. However, typically a simplified model such as a top-hat is used.

model_comparison_plot(
    module="hmf.density_field.filters",
    component="Filter",
    attr_name="filter_model",
    x="m",
    y="sigma",
    ylab=r"$\sigma(m)$",
    yratio_lab=r"$\sigma$",
    height_ratio=1 / 1.2,
    ncol_legend=1,
    ylim=(0, 10),
    ylog=False,
    xlim=(1e8, 1e16),
    comparison_limits=(0.7, 1.3),
)

Note that while any of these can be used to generate \(\sigma(m)\), the normalisation of the power spectrum, \(\sigma(R=8)\), is always calculated using the TopHat.

Mass Function

There are many fitting functions for the HMF included in hmf. Some are fitted against FOF-haloes, while others are fitted to SO halos of various definitions. Here we’ll only plot a few, so that the plot is not too confusing.

model_comparison_plot(
    module="hmf.mass_function.fitting_functions",
    component="FittingFunction",
    attr_name="hmf_model",
    x="nu",
    y="fsigma",
    ylab=r"$\nu f(\nu)$",
    height_ratio=1 / 1.2,
    ncol_legend=2,
    ylim=(5e-8, 1),
    # ylog=False,
    xlim=(6e-2, 1e1),
    comparison_limits=(0.2, 2.1),
    must_include=("PS", "SMT"),
    ignore=("ST", "Jenkins"),  # Jenkins gets in the way of the legend
    legend_size=7,
    max_models=14,
    yratio_lab="f",
)

/home/sgm/work/halos/halomod/src/halomod/halo_model.py:384: UserWarning: Requested mass definition 'FoF(l=0.2)' is not in native definitions for the 'Duffy08' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'FoF(l=0.2)'.
  return self.halo_concentration_model(
/home/sgm/work/halos/halomod/src/halomod/halo_model.py:384: UserWarning: Requested mass definition 'SOCritical(500.0)' is not in native definitions for the 'Duffy08' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOCritical(500.0)'.
  return self.halo_concentration_model(

Bias

The halo bias describes the clustering strength of halos of a given mass compared to the background matter. Many bias models are included in halomod, but even more can be accessed through the interface to colossus:

comparat = make_colossus_bias("comparat17")
bhat11 = make_colossus_bias("bhattacharya11")

model_comparison_plot(
    module="halomod.bias",
    component="Bias",
    attr_name="bias_model",
    x="nu",
    y="halo_bias",
    ylab=r"b($\nu$)",
    height_ratio=1 / 1.2,
    ncol_legend=2,
    comparison_limits=(0.5, 1.6),
    must_include=("Tinker05", "Mandelbaum05", "ST99", "SMT01", "Mo96", "Jing98"),
    ignore=("UnityBias", "Seljak04Cosmo"),
    legend_size=7,
    max_models=11,
    xlim=(6e-2, 1e1),
    ylim=(3e-1, 2e3),
    extra_models=[comparat, bhat11],
)

/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/cosmology/cosmology.py:3338: UserWarning: Astropy cosmology class contains massive neutrinos, which are not taken into account in Colossus.
  warnings.warn('Astropy cosmology class contains massive neutrinos, which are not taken into account in Colossus.')

Concentration

The concentration-mass relation affects the shape of the halo profile at different masses. halomod includes physically-motivated forms (eg. Bullock and Ludlow) and power-law forms. Since the power-law forms are really only applicable to a single cosmology, we focus on the analytic forms in this plot:

model_comparison_plot(
    module="halomod.concentration",
    component="CMRelation",
    attr_name="halo_concentration_model",
    x="m",
    y="cmz_relation",
    ylab=r"c(m)",
    yratio_lab="c",
    height_ratio=1 / 1.2,
    ncol_legend=2,
    comparison_limits=(0.3, 2.1),
    ignore=("Ludlow2016", "Ludlow2016Empirical", "Zehavi11", "Bullock01Power"),
    legend_size=7,
    ylim=(2, 3e2),
    xlim=(7e7, 5e15),
    extra_models=[
        make_colossus_cm("prada12"),
        make_colossus_cm("diemer15"),
        make_colossus_cm("diemer19"),
        #        make_colossus_cm('ishiyama20'),
    ],
)

/home/sgm/work/halos/halomod/src/halomod/halo_model.py:384: UserWarning: Requested mass definition 'SOMean(200)' is not in native definitions for the 'Bullock01' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOMean(200)'.
  return self.halo_concentration_model(
/home/sgm/work/halos/halomod/src/halomod/halo_model.py:384: UserWarning: Requested mass definition 'SOMean(200)' is not in native definitions for the 'Ludlow16' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOMean(200)'.
  return self.halo_concentration_model(
/home/sgm/work/halos/halomod/src/halomod/halo_model.py:384: UserWarning: Requested mass definition 'SOMean(200)' is not in native definitions for the 'Ludlow16Empirical' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOMean(200)'.
  return self.halo_concentration_model(
/home/sgm/work/halos/halomod/src/halomod/concentration.py:198: UserWarning: Requested mass definition 'SOMean(200)' is not in native definitions for the 'Diemer19' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOMean(200)'.
  super().__init__(*args, **kwargs)
/home/sgm/work/halos/halomod/src/halomod/concentration.py:198: UserWarning: Requested mass definition 'SOMean(200)' is not in native definitions for the 'Prada12' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOMean(200)'.
  super().__init__(*args, **kwargs)

ERROR
Requested overdensity 8.39e-12 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 7.56e-13 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 5.11e-14 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 2.51e-15 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 8.54e-17 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 1.93e-18 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 2.76e-20 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 2.34e-22 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 7.63e-12 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 6.80e-13 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 4.54e-14 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 2.20e-15 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 7.37e-17 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 1.64e-18 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 2.29e-20 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 1.90e-22 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 8.66e-25 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 2.02e-27 cannot be evaluated for scale density 2.46e-11, out of range.
ERROR
Requested overdensity 2.20e-30 cannot be evaluated for scale density 2.46e-11, out of range.

/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 1.12e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 1.26e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 1.41e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 1.58e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 1.78e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 2.00e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 2.24e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 2.51e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 2.82e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 3.16e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 3.55e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 3.98e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 4.47e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 5.01e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 5.62e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 6.31e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 7.08e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 7.94e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:421: UserWarning: Could not find concentration for model prada12, mass 8.91e+17, mdef 200m.
  warnings.warn('Could not find concentration for model %s, mass %.2e, mdef %s.' \
/home/sgm/work/halos/halomod/src/halomod/concentration.py:198: UserWarning: Requested mass definition 'SOMean(200)' is not in native definitions for the 'Diemer15' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOMean(200)'.
  super().__init__(*args, **kwargs)
/home/sgm/work/halos/halomod/src/halomod/concentration.py:198: UserWarning: Requested mass definition 'SOMean(200)' is not in native definitions for the 'Diemer19' CMRelation. No mass conversion will be performed, so results will be wrong. Using 'SOMean(200)'.
  super().__init__(*args, **kwargs)

Notice that both of the physical forms have discontinuities, which arise when the predicted redshift of formation for a given halo mass is less than zero. This typically does not present a problem for the halo model, for which intermediate masses are the most important.

Halo Exclusion

model = base_model.clone()

single_column(height_ratio=1 / 1.5)

fig, ax = plt.subplots(
    2,
    1,
    sharex=True,
    gridspec_kw={"hspace": 0, "height_ratios": [3, 2]},
    subplot_kw={"xscale": "log"},
)

for i, excl in enumerate(["NoExclusion", "Sphere", "DblSphere_", "DblEllipsoid_", "NgMatched_"]):
    model.exclusion_model = excl

    ax[0].plot(model.r, model.corr_auto_tracer, label=excl.replace("_", ""), color=f"C{i}")
    ax[0].plot(model.r, model.corr_2h_auto_tracer, ls="--", color=f"C{i}", lw=1)
    ax[1].plot(model.r, model.corr_auto_tracer / base_model.corr_auto_tracer)

model.exclusion_model = "NoExclusion"

model.hc_spectrum = "filtered-lin"
ax[0].plot(model.r, model.corr_auto_tracer, label=r"Filt. $P_{\rm lin}(k)$", color=f"C{i + 1}")
ax[0].plot(model.r, model.corr_2h_auto_tracer, ls="--", color=f"C{i + 1}", lw=1)
ax[1].plot(
    model.r,
    model.corr_auto_tracer / base_model.corr_auto_tracer,
)

model.hc_spectrum = "filtered-nl"
ax[0].plot(model.r, model.corr_auto_tracer, label=r"Filt. $P_{\rm nl}(k)$", color=f"C{i + 2}")
ax[0].plot(model.r, model.corr_2h_auto_tracer, ls="--", color=f"C{i + 2}", lw=1)
ax[1].plot(
    model.r,
    model.corr_auto_tracer / base_model.corr_auto_tracer,
)

ax[0].legend(ncol=3, frameon=False, fontsize=7, bbox_to_anchor=(0.0, 1.02, 1.0, 0.102))
ax[0].set_yscale("log")

ax[0].set_xlim(1e-1, 1e1)
ax[0].set_ylim(1e-1, 1e4)
ax[1].set_ylabel(r"$\xi / \xi_{\rm No Excl}$")
ax[0].set_ylabel(r"$\xi_{gg}(r)$")
ax[1].set_xlabel("r, [$h^{-1}$Mpc]")

/home/sgm/work/halos/halomod/src/halomod/tools.py:507: UserWarning: to use a power-law, y must be all positive or negative. Switching to zero extrapolation.
  self.lfunc = self._get_extension_func(
/home/sgm/work/halos/halomod/src/halomod/tools.py:507: UserWarning: to use a power-law, y must be all positive or negative. Switching to zero extrapolation.
  self.lfunc = self._get_extension_func(
/home/sgm/work/halos/halomod/src/halomod/tools.py:507: UserWarning: to use a power-law, y must be all positive or negative. Switching to zero extrapolation.
  self.lfunc = self._get_extension_func(
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/hmf/_internals/_cache.py:162: UserWarning: power_2h_auto_tracer_fnc not found in recalc cache.
  warnings.warn(f"{name} not found in recalc cache.")
/home/sgm/work/halos/halomod/.venv/lib/python3.12/site-packages/hmf/_internals/_cache.py:167: UserWarning: power_2h_auto_tracer_fnc not found in recalc_prop_par cache
  warnings.warn(f"{name} not found in recalc_prop_par cache")

HOD

For the HOD, we’ll only be showing models of point-source-like tracers (i.e. galaxies). A whole set of HOD models for continuous tracers (eg. HI gas) are also available. Also note that these models look very similar to each other. This is by design. Most of the models are essentially nested (i.e. you can get back Zehavi05 from setting some of the Contreras13 parameters to zero or one). By default, halomod sets the parameters so that the model returns the same or a similar occupation function. However, the HOD is not fundamental in nature – it depends on the actual galaxy sample being observed (and its depth/wavelength etc), so these models can vary widely for different parameters.

single_column(height_ratio=1 / 1.5)

model = base_model.clone()
for hod in ["Zehavi05", "Zheng05", "Tinker05", "Geach12", "Contreras13"]:
    model.hod_model = hod
    plt.plot(model.m, model.total_occupation, label=hod)

plt.xscale("log")
plt.yscale("log")
plt.ylim(1e-4, 1e4)
plt.xlim(1e9, 1e16)
plt.legend(fontsize=8)

plt.xlabel(r"m, [$h^{-1}M_\odot$]")
plt.ylabel("N(m)")

Text(0, 0.5, 'N(m)')

The key things to recognize are a low-mass cutoff, essentially corresponding to the depth of the survey, and a high-mass powerlaw with slope close to one (i.e. doubling halo mass doubles the occupancy). The differentiation between central and satellite halos causes the bump seen above at \(m\sim10^{12}\).

Profiles

two_column(height_ratio=1 / 2.3)

fig, ax = plt.subplots(
    2,
    2,
    sharex="col",
    gridspec_kw={"hspace": 0, "height_ratios": [3, 2]},
    subplot_kw={"xscale": "log"},
)

profiles = get_models(
    "halomod.profiles",
    "Profile",
    ignore=("ProfileInf", "Constant", "NFWInf", "HernquistInf", "MooreInf", "GeneralizedNFWInf"),
)

model = base_model.clone()
for profile in profiles:
    model.halo_profile_model = profile
    model.halo_profile_params = {"alpha": 1.5} if profile.__name__ == "GeneralizedNFW" else {}
    ax[0, 0].plot(
        model.r, model.halo_profile.rho(model.r, np.array([1e12])), label=profile.__name__
    )
    ax[1, 0].plot(
        model.r,
        model.halo_profile.rho(model.r, np.array([1e12]))
        / base_model.halo_profile.rho(model.r, np.array([1e12])),
    )

    ax[0, 1].plot(model.m, model.halo_profile.rho(0.03, model.m), label=profile.__name__)
    ax[1, 1].plot(
        model.m, model.halo_profile.rho(0.03, model.m) / base_model.halo_profile.rho(0.03, model.m)
    )

ax[0, 0].set_yscale("log")
ax[0, 1].set_yscale("log")

ax[0, 0].set_ylim(1e9, 1e17)
ax[0, 0].set_xlim(1e-2, 0.3)
ax[0, 0].legend(ncol=2, fontsize=8)
ax[1, 0].set_ylim(0.1, 2)

ax[0, 1].set_yscale("log")
ax[0, 1].set_ylim(1e9, 1e17)
ax[0, 1].set_xlim(1e9, 1e17)

ax[1, 0].set_xlabel(r"r, [$h^{-1}{\rm Mpc}$]")
ax[0, 0].set_ylabel(r"$\rho(r|m)$")

ax[0, 0].set_title(r"m = $10^{12} h^{-1}M_\odot$")
ax[0, 1].set_title("r = $0.03 h^{-1}$Mpc")

ax[1, 1].set_xlabel(r"m, [$h^{-1}M_\odot$]")
ax[1, 0].set_ylabel(r"$\rho/\rho_{\rm NFW}$")

/tmp/ipykernel_35768/2928401422.py:26: RuntimeWarning: invalid value encountered in divide
  model.halo_profile.rho(model.r, np.array([1e12]))
/tmp/ipykernel_35768/2928401422.py:32: RuntimeWarning: invalid value encountered in divide
  model.m, model.halo_profile.rho(0.03, model.m) / base_model.halo_profile.rho(0.03, model.m)
/tmp/ipykernel_35768/2928401422.py:26: RuntimeWarning: divide by zero encountered in divide
  model.halo_profile.rho(model.r, np.array([1e12]))
/tmp/ipykernel_35768/2928401422.py:32: RuntimeWarning: divide by zero encountered in divide
  model.m, model.halo_profile.rho(0.03, model.m) / base_model.halo_profile.rho(0.03, model.m)

Text(0, 0.5, '$\\rho/\\rho_{\\rm NFW}$')

We show here the profile for six models, plotted against mass and radius. We omit the “Inf” versions of each profile, since they are the same but extend to infinity (and are less useful for the halo model itself since some of them contain infinite mass).