Comparison to other Libraries

Comparison to other Libraries#

This notebook compares halox implementations against the colossus and gala libraries to validate the accuracy of our calculations. Each comparison illustrate how to use halox to compute properties that are also present in colossus, and shows the agreement between the two libraries and their relative differences.

Note: These comparison are implemented as unit tests in the GitHub repository; see Tests workflow.

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib_inline.backend_inline import set_matplotlib_formats

# halox imports
from halox import cosmology, nfw, hmf, lss, bias
from halox.halo import einasto

# colossus imports
from colossus.halo import profile_nfw, profile_einasto
from colossus.lss import mass_function, peaks, bias as colossus_bias
import colossus.cosmology.cosmology as cc
from colossus.utils.constants import G as Gcol

# gala for potentials
import gala.potential.potential as gp
from gala.units import galactic

print(Gcol)

jax.config.update("jax_enable_x64", True)  # double precision is, th

print(plt.style.available)
plt.style.use(["seaborn-v0_8-darkgrid", "petroff10"])
plt.rcParams.update({"xtick.direction": "in", "ytick.direction": "in"})
set_matplotlib_formats("svg")

/Users/fkeruzore/Software/halox/.venv/lib/python3.12/site-packages/jax_cosmo/__init__.py:2: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  from pkg_resources import DistributionNotFound

4.300917270038e-06
['Solarize_Light2', '_classic_test_patch', '_mpl-gallery', '_mpl-gallery-nogrid', 'bmh', 'classic', 'cmap-cividis', 'cmap-inferno', 'cmap-magma', 'cmap-plasma', 'cmap-viridis', 'colors10', 'colors10-ls', 'colors10-markers', 'colors5', 'colors5-light', 'colorsblind10', 'colorsblind34', 'dark_background', 'fast', 'fivethirtyeight', 'ggplot', 'grayscale', 'ipynb', 'ls5', 'marker7', 'nature', 'petroff10', 'seaborn-v0_8', 'seaborn-v0_8-bright', 'seaborn-v0_8-colorblind', 'seaborn-v0_8-dark', 'seaborn-v0_8-dark-palette', 'seaborn-v0_8-darkgrid', 'seaborn-v0_8-deep', 'seaborn-v0_8-muted', 'seaborn-v0_8-notebook', 'seaborn-v0_8-paper', 'seaborn-v0_8-pastel', 'seaborn-v0_8-poster', 'seaborn-v0_8-talk', 'seaborn-v0_8-ticks', 'seaborn-v0_8-white', 'seaborn-v0_8-whitegrid', 'svg_no_fonttype', 'tableau-colorblind10', 'use_mathtext', 'use_tex']

We’ll use a consistent cosmology for both libraries:

# Setup cosmologies for both libraries
# Primary cosmology: Planck18
cosmo_halox_planck = cosmology.Planck18()
cc.addCosmology(  # not taking base planck18 because it differs sightly
    "my_planck18",
    params={
        "flat": True,
        "H0": 100 * cosmo_halox_planck.h,
        "Om0": cosmo_halox_planck.Omega_m,
        "Ob0": cosmo_halox_planck.Omega_b,
        "sigma8": cosmo_halox_planck.sigma8,
        "ns": cosmo_halox_planck.n_s,
    },
)
cosmo_colossus_planck = "my_planck18"

# Alternative cosmology: Flamingo's LS8
cosmo_halox_alt = cosmology.Planck18(
    h=0.682,
    Omega_c=0.305 - 0.0473,
    Omega_b=0.0473,
    sigma8=0.76,
    n_s=0.965,
)

cc.setCosmology(cosmo_colossus_planck)
cc.addCosmology(
    "low_s8",
    params={
        "flat": True,
        "H0": 100 * cosmo_halox_alt.h,
        "Om0": cosmo_halox_alt.Omega_m,
        "Ob0": cosmo_halox_alt.Omega_b,
        "sigma8": cosmo_halox_alt.sigma8,
        "ns": cosmo_halox_alt.n_s,
    },
)
cosmo_colossus_alt = "low_s8"

def create_comparison_plot():
    fig = plt.figure(figsize=(7.7, 5.5))
    gs = gridspec.GridSpec(2, 1, height_ratios=[2, 1], hspace=0.05)
    ax_main = fig.add_subplot(gs[0])
    ax_ratio = fig.add_subplot(gs[1], sharex=ax_main)

    # Remove x-axis labels from main panel
    ax_main.tick_params(labelbottom=False)

    # Style both panels
    for ax in [ax_main, ax_ratio]:
        ax.xaxis.set_ticks_position("both")
        ax.yaxis.set_ticks_position("both")
        ax.grid(True, alpha=1.0)

    return fig, ax_main, ax_ratio


style_kw = [
    dict(color="C1", ls="-", lw=1.5, zorder=1),
    dict(color="C0", ls="--", lw=2.5, zorder=2),
    dict(color="C3", ls="-", lw=1.5, zorder=1),
    dict(color="C2", ls="--", lw=2.5, zorder=2),
    dict(color="C5", ls="-", lw=1.5, zorder=1),
    dict(color="C4", ls=":", lw=2.5, zorder=2),
    dict(color="C0", ls="-", lw=2.5, zorder=2),
    dict(color="C2", ls="-", lw=2.5, zorder=2),
    dict(color="C4", ls=":", lw=2.5, zorder=2),
]

NFW Profile Comparisons#

NFW Density Profile#

# Halo parameters
m_delta = 1e14  # h-1 Msun
c_delta = 5.0
deltas = [200.0]
redshifts = [0.0, 2.0]

# Radial range
r = jnp.logspace(-2, 1, 32)  # h-1 Mpc

# Compute density profiles for both redshifts
rho_halox = {}
rho_colossus = {}

for z in redshifts:
    nfw_halox = nfw.NFWHalo(
        m_delta, c_delta, z, cosmo_halox_planck, delta=deltas[0]
    )
    nfw_colossus = profile_nfw.NFWProfile(
        M=m_delta, c=c_delta, z=z, mdef=f"{deltas[0]:.0f}c"
    )

    rho_halox[z] = nfw_halox.density(r)
    rho_colossus[z] = nfw_colossus.density(r * 1000) * 1e9

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.loglog(r, rho_colossus[0.0], label="colossus, $z=0$", **style_kw[0])
ax_main.loglog(r, rho_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=2.0
ax_main.loglog(r, rho_colossus[2.0], label="colossus, $z=2$", **style_kw[2])
ax_main.loglog(r, rho_halox[2.0], label="halox, $z=2$", **style_kw[3])

ax_main.set_ylabel(r"Density $\rho(r)$ [$h^2 \, M_\odot \, {\rm Mpc}^{-3}$]")
ax_main.set_title("NFW Density Profile")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (rho_halox[0.0] / rho_colossus[0.0] - 1.0) * 100  # Convert to %
ratio_z1 = (rho_halox[2.0] / rho_colossus[2.0] - 1.0) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=2$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.1, 0.1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/641ce6054b303baca0444c69511c4004470ab76bdcc266171f9ae0036af172fc.svg

NFW Enclosed Mass#

# Compute enclosed mass profiles for both redshifts
mass_halox = {}
mass_colossus = {}

for z in redshifts:
    nfw_halox = nfw.NFWHalo(
        m_delta, c_delta, z, cosmo_halox_planck, delta=deltas[0]
    )
    nfw_colossus = profile_nfw.NFWProfile(
        M=m_delta, c=c_delta, z=z, mdef=f"{deltas[0]:.0f}c"
    )

    mass_halox[z] = nfw_halox.enclosed_mass(r)
    mass_colossus[z] = nfw_colossus.enclosedMass(r * 1000)

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.loglog(r, mass_colossus[0.0], label="colossus, $z=0$", **style_kw[0])
ax_main.loglog(r, mass_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=2.0
ax_main.loglog(r, mass_colossus[2.0], label="colossus, $z=2$", **style_kw[2])
ax_main.loglog(r, mass_halox[2.0], label="halox, $z=2$", **style_kw[3])

ax_main.set_ylabel(r"Enclosed Mass $M(<r)$ [$h^{-1} \, M_\odot$]")
ax_main.set_title("NFW Enclosed Mass")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (mass_halox[0.0] / mass_colossus[0.0] - 1.0) * 100  # Convert to %
ratio_z1 = (mass_halox[2.0] / mass_colossus[2.0] - 1.0) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=2$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.1, 0.1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/578072acf66cea279c692b98ca0fea0097f17fd6ce48b9319b8d5ed29ae453a4.svg

NFW Surface Density#

# Compute surface density profiles for both redshifts
sigma_surf_halox = {}
sigma_surf_colossus = {}

for z in redshifts:
    nfw_halox = nfw.NFWHalo(
        m_delta, c_delta, z, cosmo_halox_planck, delta=deltas[0]
    )
    nfw_colossus = profile_nfw.NFWProfile(
        M=m_delta, c=c_delta, z=z, mdef=f"{deltas[0]:.0f}c"
    )

    sigma_surf_halox[z] = nfw_halox.surface_density(r)
    sigma_surf_colossus[z] = nfw_colossus.surfaceDensity(r * 1000) * 1e6

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.loglog(
    r, sigma_surf_colossus[0.0], label="colossus, $z=0$", **style_kw[0]
)
ax_main.loglog(r, sigma_surf_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=2.0
ax_main.loglog(
    r, sigma_surf_colossus[2.0], label="colossus, $z=2$", **style_kw[2]
)
ax_main.loglog(r, sigma_surf_halox[2.0], label="halox, $z=2$", **style_kw[3])

ax_main.set_ylabel(
    r"Surface Density $\Sigma(r)$ [$h \, M_\odot \, {\rm Mpc}^{-2}$]"
)
ax_main.set_title("NFW Surface Density")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (
    sigma_surf_halox[0.0] / sigma_surf_colossus[0.0] - 1.0
) * 100  # Convert to %
ratio_z1 = (
    sigma_surf_halox[2.0] / sigma_surf_colossus[2.0] - 1.0
) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=1$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Projected Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.1, 0.1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/291efb10c762127640d61eaa52f77ae083d92b83f5d08975219f04f02f4d5159.svg

NFW Circular Velocity#

# Compute circular velocity profiles for both redshifts
v_c_halox = {}
v_c_colossus = {}

for z in redshifts:
    nfw_halox = nfw.NFWHalo(
        m_delta, c_delta, z, cosmo_halox_planck, delta=deltas[0]
    )
    nfw_colossus = profile_nfw.NFWProfile(
        M=m_delta, c=c_delta, z=z, mdef=f"{deltas[0]:.0f}c"
    )

    v_c_halox[z] = nfw_halox.circular_velocity(r)
    v_c_colossus[z] = nfw_colossus.circularVelocity(r * 1000)

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.semilogx(r, v_c_colossus[0.0], label="colossus, $z=0$", **style_kw[0])
ax_main.semilogx(r, v_c_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=2.0
ax_main.semilogx(r, v_c_colossus[2.0], label="colossus, $z=2$", **style_kw[2])
ax_main.semilogx(r, v_c_halox[2.0], label="halox, $z=2$", **style_kw[3])

ax_main.set_ylabel(r"Circular Velocity $v_c(r)$ [${\rm km} \, {\rm s}^{-1}$]")
ax_main.set_title("NFW Circular Velocity")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (v_c_halox[0.0] / v_c_colossus[0.0] - 1.0) * 100  # Convert to %
ratio_z1 = (v_c_halox[2.0] / v_c_colossus[2.0] - 1.0) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=1$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.1, 0.1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/a5b11e4a17670ef30044a6880ae8d152d6fef023af35cba6e8152495ae61d8da.svg

NFW Potential#

import astropy.units as u
import numpy as np

phi_gala = {}
phi_halox = {}

kmperkpc = 30856775999999956
secpermyr = 3.15576e13
confactor = kmperkpc**2 / secpermyr**2
# print(confactor) #conversion from kpc^2 per Myr^2 to km^2 per s^2

for z in redshifts:
    haloxhalo = nfw.NFWHalo(m_delta, c_delta, z, cosmo_halox_planck)
    m = (m_delta / (np.log(1 + c_delta) - c_delta / (1 + c_delta))) * u.Msun
    r_s = (
        haloxhalo.r_delta * 1000 / c_delta * u.kpc
    )  # converting r_delta to kpc
    galahalo = gp.NFWPotential(m=m, r_s=r_s, units=galactic)

    r_kpc = r * 1000 * u.kpc  # if r was in Mpc
    xyz = np.zeros((3, len(r_kpc))) * u.kpc
    xyz[0] = r_kpc
    phi_halox[z] = haloxhalo.potential(r)
    phi_gala[z] = galahalo.energy(xyz).value * confactor

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.semilogx(r, phi_gala[0.0], label="gala, $z=0$", **style_kw[0])
ax_main.semilogx(r, phi_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=2.0
ax_main.semilogx(r, phi_gala[2.0], label="gala, $z=2$", **style_kw[2])
ax_main.semilogx(r, phi_halox[2.0], label="halox, $z=2$", **style_kw[3])

ax_main.set_ylabel(r"Potential $\phi(r)$ [${\rm km^2} \, {\rm s}^{-2}$]")
ax_main.set_title("NFW Potential")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (phi_halox[0.0] / phi_gala[0.0] - 1.0) * 100  # Convert to %
ratio_z1 = (phi_halox[2.0] / phi_gala[2.0] - 1.0) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=2$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - gala [%]")
ax_ratio.set_ylim(-0.1, 0.1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/7e4fe6012c69b880517710bcbf0b42519ff1c239ce6e681ff6bb7e9b51ecad14.svg

Einasto Profile Comparisons#

Einasto Density Profile#

# Halo parameters
m_delta = 1e14  # h-1 Msun
c_delta = 5.0
alpha = 0.5
deltas = [200.0]
redshifts = [0.0, 2.0]

# Radial range
r = jnp.logspace(-3, 1, 128)  # h-1 Mpc

# Compute density profiles for both redshifts
rho_halox = {}
rho_colossus = {}
for z in redshifts:
    einasto_halox = einasto.EinastoHalo(
        m_delta, c_delta, z, alpha, cosmo_halox_planck, delta=deltas[0]
    )
    einasto_colossus = profile_einasto.EinastoProfile(
        M=m_delta, c=c_delta, z=z, mdef=f"{deltas[0]:.0f}c", alpha=alpha
    )

    rho_halox[z] = einasto_halox.density(r)
    rho_colossus[z] = einasto_colossus.density(r * 1000) * 1e9

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.loglog(r, rho_colossus[0.0], label="colossus, $z=0$", **style_kw[0])
ax_main.loglog(r, rho_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=1.0
ax_main.loglog(r, rho_colossus[2.0], label="colossus, $z=2$", **style_kw[2])
ax_main.loglog(r, rho_halox[2.0], label="halox, $z=2$", **style_kw[3])

ax_main.set_ylabel(r"Density $\rho(r)$ [$h^2 \, M_\odot \, {\rm Mpc}^{-3}$]")
ax_main.set_title("Einasto Density Profile")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (rho_halox[0.0] / rho_colossus[0.0] - 1.0) * 100  # Convert to %
ratio_z1 = (rho_halox[2.0] / rho_colossus[2.0] - 1.0) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=2$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-1, 1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/7374e35d5ceffc0932c727f51219dbe49691de163bf3c2546ee2daab2688da5b.svg

Einasto Enclosed Mass#

# Compute enclosed mass profiles for both redshifts
mass_halox = {}
mass_colossus = {}

for z in redshifts:
    einasto_halox = einasto.EinastoHalo(
        m_delta, c_delta, z, alpha, cosmo_halox_planck, delta=deltas[0]
    )
    einasto_colossus = profile_einasto.EinastoProfile(
        M=m_delta, c=c_delta, z=z, mdef=f"{deltas[0]:.0f}c", alpha=alpha
    )

    mass_halox[z] = einasto_halox.enclosed_mass(r)
    mass_colossus[z] = einasto_colossus.enclosedMass(r * 1000)

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.loglog(r, mass_colossus[0.0], label="colossus, $z=0$", **style_kw[0])
ax_main.loglog(r, mass_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=2.0
ax_main.loglog(r, mass_colossus[2.0], label="colossus, $z=1$", **style_kw[2])
ax_main.loglog(r, mass_halox[2.0], label="halox, $z=1$", **style_kw[3])

ax_main.set_ylabel(r"Enclosed Mass $M(<r)$ [$h^{-1} \, M_\odot$]")
ax_main.set_title("Einasto Enclosed Mass")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (mass_halox[0.0] / mass_colossus[0.0] - 1.0) * 100  # Convert to %
ratio_z1 = (mass_halox[2.0] / mass_colossus[2.0] - 1.0) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=1$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.1, 0.1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/2261d185cd9c80cc326fb6455b7e6e500eb9284a222f08cefac432f449384b83.svg

Einasto Circular Velocity#

# Compute circular velocity profiles for both redshifts
v_c_halox = {}
v_c_colossus = {}

for z in redshifts:
    einasto_halox = einasto.EinastoHalo(
        m_delta, c_delta, z, alpha, cosmo_halox_planck, delta=deltas[0]
    )
    einasto_colossus = profile_einasto.EinastoProfile(
        M=m_delta, c=c_delta, z=z, mdef=f"{deltas[0]:.0f}c", alpha=alpha
    )

    v_c_halox[z] = einasto_halox.circular_velocity(r)
    v_c_colossus[z] = einasto_colossus.circularVelocity(r * 1000)

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.semilogx(r, v_c_colossus[0.0], label="colossus, $z=0$", **style_kw[0])
ax_main.semilogx(r, v_c_halox[0.0], label="halox, $z=0$", **style_kw[1])

# z=2.0
ax_main.semilogx(r, v_c_colossus[2.0], label="colossus, $z=2$", **style_kw[2])
ax_main.semilogx(r, v_c_halox[2.0], label="halox, $z=2$", **style_kw[3])

ax_main.set_ylabel(r"Circular Velocity $v_c(r)$ [${\rm km} \, {\rm s}^{-1}$]")
ax_main.set_title("Einasto Circular Velocity")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (v_c_halox[0.0] / v_c_colossus[0.0] - 1.0) * 100  # Convert to %
ratio_z1 = (v_c_halox[2.0] / v_c_colossus[2.0] - 1.0) * 100  # Convert to %

ax_ratio.semilogx(r, ratio_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(r, ratio_z1, label="$z=2$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Radius $r$ [$h^{-1} \, {\rm Mpc}$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.1, 0.1)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/6ef255dc7ff9be09ca2806c15d354a83ff1ec7120bc1762a0700cb701dcc888e.svg

LSS and Halo Mass Function Comparisons#

RMS Variance σ(R,z)#

# Mass range and redshifts
masses = jnp.logspace(10, 15.5, 32)  # h-1 Msun
redshifts = [0.0, 1.0]

# Compute sigma values for both redshifts and cosmologies
sigma_halox_planck = {}
sigma_halox_alt = {}
sigma_colossus_planck = {}
sigma_colossus_alt = {}


@jax.jit
def compute_sigma_halox(masses, z, cosmo):
    """JIT-compiled batch computation of sigma values"""
    R = lss.mass_to_lagrangian_radius(masses, cosmo)
    return lss.sigma_R(R, z, cosmo)


for z in redshifts:
    # Planck cosmology
    sigma_halox_planck[z] = compute_sigma_halox(masses, z, cosmo_halox_planck)
    cosmo_col_obj_planck = cc.setCosmology(cosmo_colossus_planck)
    sigma_colossus_planck[z] = cosmo_col_obj_planck.sigma(
        peaks.lagrangianR(masses), z=z
    )

    # Alternative cosmology
    sigma_halox_alt[z] = compute_sigma_halox(masses, z, cosmo_halox_alt)
    cosmo_col_obj_alt = cc.setCosmology(cosmo_colossus_alt)
    sigma_colossus_alt[z] = cosmo_col_obj_alt.sigma(
        peaks.lagrangianR(masses), z=z
    )

# Reset to Planck cosmology for subsequent cells
cc.setCosmology(cosmo_colossus_planck)

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# Original Planck cosmology comparisons - z=0.0
ax_main.loglog(
    masses, sigma_colossus_planck[0.0], label="colossus, $z=0$", **style_kw[0]
)
ax_main.loglog(
    masses, sigma_halox_planck[0.0], label="halox, $z=0$", **style_kw[1]
)

# Original Planck cosmology comparisons - z=1.0
ax_main.loglog(
    masses, sigma_colossus_planck[1.0], label="colossus, $z=1$", **style_kw[2]
)
ax_main.loglog(
    masses, sigma_halox_planck[1.0], label="halox, $z=1$", **style_kw[3]
)

# Alternative cosmology - z=0 only
ax_main.loglog(
    masses,
    sigma_colossus_alt[0.0],
    label="colossus Low-S8, $z=0$",
    **style_kw[4],
)
ax_main.loglog(
    masses, sigma_halox_alt[0.0], label="halox Low-S8, $z=0$", **style_kw[5]
)

ax_main.set_ylabel(r"RMS Variance $\sigma(M,z)$")
ax_main.set_title("RMS Variance of Density Fluctuations")
ax_main.legend()

# Ratio panel: relative differences
# Planck cosmology ratios
ratio_planck_z0 = (
    sigma_halox_planck[0.0] / sigma_colossus_planck[0.0] - 1.0
) * 100
ratio_planck_z1 = (
    sigma_halox_planck[1.0] / sigma_colossus_planck[1.0] - 1.0
) * 100

# Alternative cosmology ratio (z=0 only)
ratio_alt_z0 = (sigma_halox_alt[0.0] / sigma_colossus_alt[0.0] - 1.0) * 100

ax_ratio.semilogx(masses, ratio_planck_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(masses, ratio_planck_z1, label="$z=1$", **style_kw[-2])
ax_ratio.semilogx(masses, ratio_alt_z0, label="Low-S8 $z=0$", **style_kw[-1])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Mass $M_{200c}$ [$h^{-1} \, M_\odot$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.3, 0.3)
ax_ratio.legend(ncol=3)

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/2c54d38f6ebac8784dcc3e1875172b750d09194f96136d68fed56bedd3262278.svg

Tinker08 Mass Function dn/d ln M#

# Mass range and redshifts
masses = jnp.logspace(11, 15.5, 32)  # h-1 Msun
redshifts = [0.0, 1.0]
delta_c = 200.0

# JIT-compiled function for mass function computation
halox_compute_hmf = jax.jit(hmf.tinker08_mass_function)

# Compute mass functions for both redshifts and cosmologies
dn_dlnM_halox_planck = {}
dn_dlnM_halox_alt = {}
dn_dlnM_colossus_planck = {}
dn_dlnM_colossus_alt = {}

for z in redshifts:
    # Planck cosmology
    dn_dlnM_halox_planck[z] = halox_compute_hmf(
        masses, z, cosmo_halox_planck, delta_c
    )
    cc.setCosmology(cosmo_colossus_planck)
    dn_dlnM_colossus_planck[z] = mass_function.massFunction(
        masses,
        z,
        mdef=f"{delta_c:.0f}c",
        model="tinker08",
        q_in="M",
        q_out="dndlnM",
    )

    # Alternative cosmology
    dn_dlnM_halox_alt[z] = halox_compute_hmf(
        masses, z, cosmo_halox_alt, delta_c
    )
    cc.setCosmology(cosmo_colossus_alt)
    dn_dlnM_colossus_alt[z] = mass_function.massFunction(
        masses,
        z,
        mdef=f"{delta_c:.0f}c",
        model="tinker08",
        q_in="M",
        q_out="dndlnM",
    )

# Reset to Planck cosmology
cc.setCosmology(cosmo_colossus_planck)

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# Original Planck cosmology comparisons - z=0.0
ax_main.loglog(
    masses,
    dn_dlnM_colossus_planck[0.0],
    label="colossus, $z=0$",
    **style_kw[0],
)
ax_main.loglog(
    masses, dn_dlnM_halox_planck[0.0], label="halox, $z=0$", **style_kw[1]
)

# Original Planck cosmology comparisons - z=1.0
ax_main.loglog(
    masses,
    dn_dlnM_colossus_planck[1.0],
    label="colossus, $z=1$",
    **style_kw[2],
)
ax_main.loglog(
    masses, dn_dlnM_halox_planck[1.0], label="halox, $z=1$", **style_kw[3]
)

# Alternative cosmology - z=0 only
ax_main.loglog(
    masses,
    dn_dlnM_colossus_alt[0.0],
    label="colossus Low-S8, $z=0$",
    **style_kw[4],
)
ax_main.loglog(
    masses, dn_dlnM_halox_alt[0.0], label="halox Low-S8, $z=0$", **style_kw[5]
)

ax_main.set_ylabel(
    r"Mass Function $\frac{{\rm d}n}{{\rm d}\ln M}$ [$h^3 \, {\rm Mpc}^{-3}$]"
)
ax_main.set_title("Tinker08 Mass Function")
ax_main.legend()

# Ratio panel: relative differences
# Planck cosmology ratios
ratio_planck_z0 = (
    dn_dlnM_halox_planck[0.0] / dn_dlnM_colossus_planck[0.0] - 1.0
) * 100
ratio_planck_z1 = (
    dn_dlnM_halox_planck[1.0] / dn_dlnM_colossus_planck[1.0] - 1.0
) * 100

# Alternative cosmology ratio (z=0 only)
ratio_alt_z0 = (dn_dlnM_halox_alt[0.0] / dn_dlnM_colossus_alt[0.0] - 1.0) * 100

ax_ratio.semilogx(masses, ratio_planck_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(masses, ratio_planck_z1, label="$z=1$", **style_kw[-2])
ax_ratio.semilogx(masses, ratio_alt_z0, label="Low-S8 $z=0$", **style_kw[-1])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Mass $M_{200c}$ [$h^{-1} \, M_\odot$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-3, 3)
ax_ratio.legend(ncol=3)

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/0d579e30a756aeba5596514d8ad020470049d8650e72cf1e568dce4f8319ec31.svg

Tinker10 Halo Bias#

# Mass range and redshifts
masses = jnp.logspace(10, 15.5, 32)  # h-1 Msun
redshifts = [0.0, 1.0]
delta_c = 200.0

# JIT-compiled function for halo bias computation
halox_compute_bias = jax.jit(bias.tinker10_bias)

# Compute halo bias for both redshifts and cosmologies
bias_halox_planck = {}
bias_halox_alt = {}
bias_colossus_planck = {}
bias_colossus_alt = {}

for z in redshifts:
    # Planck cosmology
    bias_halox_planck[z] = halox_compute_bias(
        masses, z, cosmo_halox_planck, delta_c
    )
    cc.setCosmology(cosmo_colossus_planck)
    bias_colossus_planck[z] = colossus_bias.haloBias(
        masses, z, mdef=f"{delta_c:.0f}c", model="tinker10"
    )

    # Alternative cosmology
    bias_halox_alt[z] = halox_compute_bias(masses, z, cosmo_halox_alt, delta_c)
    cc.setCosmology(cosmo_colossus_alt)
    bias_colossus_alt[z] = colossus_bias.haloBias(
        masses, z, mdef=f"{delta_c:.0f}c", model="tinker10"
    )

# Reset to Planck cosmology
cc.setCosmology(cosmo_colossus_planck)

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# Original Planck cosmology comparisons - z=0.0
ax_main.loglog(
    masses, bias_colossus_planck[0.0], label="colossus, $z=0$", **style_kw[0]
)
ax_main.loglog(
    masses, bias_halox_planck[0.0], label="halox, $z=0$", **style_kw[1]
)

# Original Planck cosmology comparisons - z=1.0
ax_main.loglog(
    masses, bias_colossus_planck[1.0], label="colossus, $z=1$", **style_kw[2]
)
ax_main.loglog(
    masses, bias_halox_planck[1.0], label="halox, $z=1$", **style_kw[3]
)

# Alternative cosmology - z=0 only
ax_main.loglog(
    masses,
    bias_colossus_alt[0.0],
    label="colossus Low-S8, $z=0$",
    **style_kw[4],
)
ax_main.loglog(
    masses, bias_halox_alt[0.0], label="halox Low-S8, $z=0$", **style_kw[5]
)

ax_main.set_ylabel(r"Halo Bias $b(M,z)$")
ax_main.set_title("Tinker10 Halo Bias")
ax_main.legend()

# Ratio panel: relative differences
# Planck cosmology ratios
ratio_planck_z0 = (
    bias_halox_planck[0.0] / bias_colossus_planck[0.0] - 1.0
) * 100
ratio_planck_z1 = (
    bias_halox_planck[1.0] / bias_colossus_planck[1.0] - 1.0
) * 100

# Alternative cosmology ratio (z=0 only)
ratio_alt_z0 = (bias_halox_alt[0.0] / bias_colossus_alt[0.0] - 1.0) * 100

ax_ratio.semilogx(masses, ratio_planck_z0, label="$z=0$", **style_kw[-3])
ax_ratio.semilogx(masses, ratio_planck_z1, label="$z=1$", **style_kw[-2])
ax_ratio.semilogx(masses, ratio_alt_z0, label="Low-S8 $z=0$", **style_kw[-1])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Mass $M_{200c}$ [$h^{-1} \, M_\odot$]")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.5, 0.5)
ax_ratio.legend(ncol=3)

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/2d01a519e831cdbb480bf2cfb452e31b87863926ed7613d5496ab7d23c470377.svg

Mass-concentration relations#

from colossus.cosmology import cosmology
from colossus.halo import concentration
from colossus.halo import mass_defs

from halox.cm import duffy08, prada12, klypin11, child18all, child18relaxed

mray = jnp.logspace(11, 15, 32)

zray = jnp.array([0.0, 1.0])


cduffy_col = np.zeros((len(zray), len(mray)))
cprada_col = np.zeros_like(cduffy_col)
ccall_col = np.zeros_like(cduffy_col)
ccrel_col = np.zeros_like(cduffy_col)


def vectorize_model(model, mray, zray):
    return jax.vmap(lambda z: model(mray, z))(zray)


cduffy = vectorize_model(duffy08(), mray, zray)
cprada = vectorize_model(prada12(cosmo_halox_planck), mray, zray)
ccall = vectorize_model(child18all(cosmo_halox_planck), mray, zray)
ccrel = vectorize_model(child18relaxed(cosmo_halox_planck), mray, zray)

for i, z in enumerate(zray):
    cduffy_col[i] = concentration.concentration(
        mray, "200c", z, model="duffy08"
    )

    cprada_col[i] = concentration.concentration(
        mray, "200c", z, model="prada12"
    )

    ccall_col[i] = concentration.concentration(
        mray, "200c", z, model="child18", halo_sample="individual_all"
    )

    ccrel_col[i] = concentration.concentration(
        mray, "200c", z, model="child18", halo_sample="individual_relaxed"
    )

# halox klypin calculation
x = cosmo_halox_planck.Omega_m - 1
# Below: Brian and Norman virial overdensity calculation
deltaBN = 18 * jnp.pi**2 + 82 * x - 39 * x**2
mklyp, rklyp, cklyp = nfw.delta_delta(
    mray,
    cduffy[0],
    z=zray[0] * jnp.ones(len(mray)),
    cosmo=cosmo_halox_planck,
    delta_old=200,
    delta_new=deltaBN,
)
cklypin = klypin11()(M=mklyp)

# colossus klypin calculation
# Brian & Norman overdensity
x = cc.getCurrent().Om(0) - 1
deltaBN = 18 * np.pi**2 + 82 * x - 39 * x**2

# Convert Duffy (z=0 slice like yours)
mklyp_col, rklyp, cklyp_conv = mass_defs.changeMassDefinition(
    mray, cprada[0], zray[0], "200c", "vir"
)

cklypin_col = concentration.concentration(
    mklyp_col, "vir", zray[0], model="klypin11"
)

halox_models = {
    "duffy08": cduffy,
    "prada12": cprada,
    "child18all": ccall,
    "child18relaxed": ccrel,
}

colossus_models = {
    "duffy08": cduffy_col,
    "prada12": cprada_col,
    "child18all": ccall_col,
    "child18relaxed": ccrel_col,
}

/Users/fkeruzore/Software/halox/.venv/lib/python3.12/site-packages/colossus/halo/concentration.py:442: UserWarning: Some masses or redshifts are outside the validity of the concentration model.
  warnings.warn('Some masses or redshifts are outside the validity of the concentration model.')

# Create comparison plot
fig, ax_main, ax_ratio = create_comparison_plot()

# Main panel: plot colossus first (behind), then halox (on top)
# z=0.0
ax_main.loglog(mklyp_col, cklypin_col, label="colossus, $z=0$", **style_kw[0])
ax_main.loglog(mklyp, cklypin, label="halox, $z=0$", **style_kw[1])

ax_main.set_ylabel(r"$c_{200}$")
ax_main.set_title("Klypin 2011 et al. cM relation")
ax_main.legend()

# Ratio panel: relative differences for both redshifts
ratio_z0 = (cklypin / cklypin_col - 1.0) * 100  # Convert to %

ax_ratio.semilogx(mray, ratio_z0, label="$z=0$", **style_kw[-3])
# ax_ratio.semilogx(mray, ratio_z1, label="$z=0.5$", **style_kw[-2])

ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
ax_ratio.set_xlabel(r"Mass $[M_\odot]$")
ax_ratio.set_ylabel("halox - colossus [%]")
ax_ratio.set_ylim(-0.5, 0.5)
ax_ratio.legend()

# Align y-axis labels
fig.align_labels([ax_main, ax_ratio])
../_images/4887b98bc33e0f1c1580b6d85568ec1f6274617a835bd6dff73648068e3da28a.svg 
for name in halox_models.keys():
    fig, ax_main, ax_ratio = create_comparison_plot()

    halox_vals = halox_models[name]
    col_vals = colossus_models[name]

    for i, z in enumerate(zray):
        # Main panel
        ax_main.loglog(
            mray, col_vals[i], label=f"colossus, z={z}", **style_kw[2 * i]
        )
        ax_main.loglog(
            mray, halox_vals[i], label=f"halox, z={z}", **style_kw[2 * i + 1]
        )

        # Ratio panel
        ratio = (halox_vals[i] / col_vals[i] - 1.0) * 100
        ax_ratio.semilogx(mray, ratio, label=f"z={z}", **style_kw[-3 + i])

    ax_main.set_ylabel(r"$c_{200}$")
    ax_main.set_title(f"{name} c-M relation")
    ax_main.legend()

    ax_ratio.axhline(0.0, color="k", linestyle="--", alpha=0.5)
    ax_ratio.set_xlabel(r"Mass $[M_\odot]$")
    ax_ratio.set_ylabel("halox - colossus [%]")
    ax_ratio.set_ylim(-0.1, 0.1)
    ax_ratio.legend()

    fig.align_labels([ax_main, ax_ratio])
../_images/2def4b883f6c67794db0efbaf7965992ef9d293d7214805e25fd52146a87efc1.svg ../_images/559d84fa1efbd14aaa42742db81c8556c57528f13de75e53ee7e528f6bd151ab.svg ../_images/916142d57007cdac9f92f7e51322119dd83b8fa028015dc4608f1ac5ccef9972.svg ../_images/d89631a4b1045f4a32e45103f1bd98d292b85be2963f109a07731bdbc8c355f9.svg