Skip to content

Halos

Add dark matter halos to a density field, using Poisson statistics.

HaloDistribution

Source code in fastbox/halos.py 
class HaloDistribution(object):

    def __init__(self, box, mass_range, mass_bins):
        """
        An object to manage the placing of halos on top of a realisation of a 
        density field in a box.

        Parameters:
            box (CosmoBox):
                Object containing a simulation box.

            mass_range (tuple):
                Range of halo masses, in Msun.

            mass_bins (int):
                Number of mass bins.
        """
        self.box = box
        self.Mmin, self.Mmax = mass_range
        self.mass_bins = mass_bins


    def construct_bins(self, z):
        """
        Construct a binned halo mass function and bias function. This function 
        sets the ``self.dndlog10M`` and ``self.bias`` arrays.

        Parameters:
            z (float): Redshift at which to construct the mass function and bias.
        """
        a = 1./(1.+z)

        # Define mass bins
        Mh_edges = np.logspace(np.log10(self.Mmin), 
                               np.log10(self.Mmax), 
                               int(self.mass_bins)+1)
        Mh_centres = 0.5 * (Mh_edges[1:] + Mh_edges[:-1])

        # Get mass function and bias at mass bin centres
        self.dndlog10M = ccl.massfunction.massfunc(self.box.cosmo, Mh_centres, 
                                              a, overdensity=200)
        self.bias = ccl.halo_bias(cosmo, Mh_centres, a, overdensity=200)


    def halo_count_field(self, delta_x, nbar, bias, lognormal=False):
        """
        Generate a cube of halo counts, with a Poisson-generated realisation of 
        the number of halos per voxel.

        The expected number of halos per voxel is given by:

            N_halo(x) = vol_voxel * nbar(x) * (1 + bias(x) * delta_x(x))

        Parameters:
            delta_x (array_like):
                3D box of density contrast values.

            nbar (float or array_like):
                Mean number density (per unit comoving volume, Mpc^-3) of a chosen 
                (sub-)population of halos. 

                This can be specified as a constant across the whole box, a 1D 
                array which will be interpreted as a function of the z/LOS 
                direction, or a 3D array with the value of nbar per voxel.

            bias (float or array_like):
                Mean bias of the population of halos.

                This can be specified as a constant across the whole box, a 1D 
                array which will be interpreted as a function of the z/LOS 
                direction, or a 3D array with the value of nbar per voxel.

            lognormal (bool, optional):
                If True, apply a log-normal transform to the biased density field 
                before realising the halo counts. This should normally only be 
                applied to the *linear* density contrast.

        Returns:
            N_halo (array_like):
                3D box of halo counts per voxel.
        """
        # Build bias and number density functions with the right shape
        nbar = np.atleast_1d(nbar)
        bias = np.atleast_1d(bias)

        # Check for 1D shape (assumed to be either constant or fn. of redshift)
        if len(nbar.shape) == 1:
            nbar = nbar[np.newaxis,np.newaxis,:]
        if len(bias.shape) == 1:
            bias = bias[np.newaxis,np.newaxis,:]

        # Comoving volume of each voxel
        voxel_vol = self.box.Lx * self.box.Ly * self.box.Lz / self.box.N**3.

        # Calculate halo overdensity field
        delta_h = bias * delta_x
        if lognormal:
            delta_h = np.exp(delta_h)
            delta_h /= np.mean(delta_h)
            delta_h -= 1.

        # Expected number of halos in each voxel
        Nbar = voxel_vol * nbar * (1. + delta_h)
        if not lognormal:
            Nbar[np.where(Nbar < 0.)] = 0.

        # Realise halo count in each voxel
        Nhalo = np.random.poisson(lam=np.nan_to_num(Nbar))
        return Nhalo


    def realise_halo_catalogue(self, Nhalo, scatter=False, scatter_type='uniform'):
        """
        Generate a catalogue of halo positions (comoving) within the 3D box, 
        using a previously-calculated field of halo counts per voxel.

        Parameters:
            Nhalo (array_like):
                3D array containing the number of halos of this type in each voxel, 
                typically generated by `self.halo_count_field()`.

            scatter (bool, optional):
                Whether to randomly shift the positions of the halos within a 
                voxel. If False, halos in the same voxel have the same position. 

            scatter_type (str, optional):
                If `scatter=True`, how the positions within the voxel should be 
                assigned. Default: 'uniform' (uniform random positions).

        Returns:
            cat (array_like):
                Catalogue of halo positions, of shape (Nhalos, 3) (comoving 
                Cartesian positions x,y,z).
        """
        cat_x, cat_y, cat_z = [], [], []

        # Loop over sets of voxels with the same count
        unique_counts = np.unique(Nhalo)
        for count in unique_counts:
            if count < 1.: continue

            # Find voxels containing this number of halos and add to catalogue 
            # (repeated 'count' times)
            idx_x, idx_y, idx_z = np.where(Nhalo == count)
            cat_x.append(np.repeat(idx_x, count))
            cat_y.append(np.repeat(idx_y, count))
            cat_z.append(np.repeat(idx_z, count))

        # Concatenate into a single array
        cat = np.column_stack( (np.concatenate(cat_x),
                                np.concatenate(cat_y),
                                np.concatenate(cat_z)) ).astype(np.float64)

        # Add scatter in positions
        if scatter:
            if scatter_type == 'uniform':
                # Uniform scatter within each voxel
                cat += np.random.uniform(0., 1.-1e-8, cat.size).reshape(cat.shape)
            else:
                raise ValueError("scatter_type='%s' not recognised" 
                                  % scatter_type)

        # Scale coordinates with box size
        cat[:,0] *= self.box.Lx / self.box.N
        cat[:,1] *= self.box.Ly / self.box.N
        cat[:,2] *= self.box.Lz / self.box.N

        return cat

__init__(self, box, mass_range, mass_bins) special

An object to manage the placing of halos on top of a realisation of a density field in a box.

Parameters:

Name Type Description Default
box CosmoBox

Object containing a simulation box.

 required
mass_range tuple

Range of halo masses, in Msun.

 required
mass_bins int

Number of mass bins.

 required
Source code in fastbox/halos.py 
def __init__(self, box, mass_range, mass_bins):
    """
    An object to manage the placing of halos on top of a realisation of a 
    density field in a box.

    Parameters:
        box (CosmoBox):
            Object containing a simulation box.

        mass_range (tuple):
            Range of halo masses, in Msun.

        mass_bins (int):
            Number of mass bins.
    """
    self.box = box
    self.Mmin, self.Mmax = mass_range
    self.mass_bins = mass_bins

construct_bins(self, z)

Construct a binned halo mass function and bias function. This function sets the self.dndlog10M and self.bias arrays.

Parameters:

Name Type Description Default
z float

Redshift at which to construct the mass function and bias.

 required
Source code in fastbox/halos.py 
def construct_bins(self, z):
    """
    Construct a binned halo mass function and bias function. This function 
    sets the ``self.dndlog10M`` and ``self.bias`` arrays.

    Parameters:
        z (float): Redshift at which to construct the mass function and bias.
    """
    a = 1./(1.+z)

    # Define mass bins
    Mh_edges = np.logspace(np.log10(self.Mmin), 
                           np.log10(self.Mmax), 
                           int(self.mass_bins)+1)
    Mh_centres = 0.5 * (Mh_edges[1:] + Mh_edges[:-1])

    # Get mass function and bias at mass bin centres
    self.dndlog10M = ccl.massfunction.massfunc(self.box.cosmo, Mh_centres, 
                                          a, overdensity=200)
    self.bias = ccl.halo_bias(cosmo, Mh_centres, a, overdensity=200)

halo_count_field(self, delta_x, nbar, bias, lognormal=False)

Generate a cube of halo counts, with a Poisson-generated realisation of the number of halos per voxel.

The expected number of halos per voxel is given by:

N_halo(x) = vol_voxel * nbar(x) * (1 + bias(x) * delta_x(x))

Parameters:

Name Type Description Default
delta_x array_like

3D box of density contrast values.

 required
nbar float or array_like

Mean number density (per unit comoving volume, Mpc^-3) of a chosen (sub-)population of halos.

This can be specified as a constant across the whole box, a 1D array which will be interpreted as a function of the z/LOS direction, or a 3D array with the value of nbar per voxel.

 required
bias float or array_like

Mean bias of the population of halos.

This can be specified as a constant across the whole box, a 1D array which will be interpreted as a function of the z/LOS direction, or a 3D array with the value of nbar per voxel.

 required
lognormal bool

If True, apply a log-normal transform to the biased density field before realising the halo counts. This should normally only be applied to the linear density contrast.

 False

Returns:

Type Description
N_halo (array_like)

3D box of halo counts per voxel.
Source code in fastbox/halos.py 
def halo_count_field(self, delta_x, nbar, bias, lognormal=False):
    """
    Generate a cube of halo counts, with a Poisson-generated realisation of 
    the number of halos per voxel.

    The expected number of halos per voxel is given by:

        N_halo(x) = vol_voxel * nbar(x) * (1 + bias(x) * delta_x(x))

    Parameters:
        delta_x (array_like):
            3D box of density contrast values.

        nbar (float or array_like):
            Mean number density (per unit comoving volume, Mpc^-3) of a chosen 
            (sub-)population of halos. 

            This can be specified as a constant across the whole box, a 1D 
            array which will be interpreted as a function of the z/LOS 
            direction, or a 3D array with the value of nbar per voxel.

        bias (float or array_like):
            Mean bias of the population of halos.

            This can be specified as a constant across the whole box, a 1D 
            array which will be interpreted as a function of the z/LOS 
            direction, or a 3D array with the value of nbar per voxel.

        lognormal (bool, optional):
            If True, apply a log-normal transform to the biased density field 
            before realising the halo counts. This should normally only be 
            applied to the *linear* density contrast.

    Returns:
        N_halo (array_like):
            3D box of halo counts per voxel.
    """
    # Build bias and number density functions with the right shape
    nbar = np.atleast_1d(nbar)
    bias = np.atleast_1d(bias)

    # Check for 1D shape (assumed to be either constant or fn. of redshift)
    if len(nbar.shape) == 1:
        nbar = nbar[np.newaxis,np.newaxis,:]
    if len(bias.shape) == 1:
        bias = bias[np.newaxis,np.newaxis,:]

    # Comoving volume of each voxel
    voxel_vol = self.box.Lx * self.box.Ly * self.box.Lz / self.box.N**3.

    # Calculate halo overdensity field
    delta_h = bias * delta_x
    if lognormal:
        delta_h = np.exp(delta_h)
        delta_h /= np.mean(delta_h)
        delta_h -= 1.

    # Expected number of halos in each voxel
    Nbar = voxel_vol * nbar * (1. + delta_h)
    if not lognormal:
        Nbar[np.where(Nbar < 0.)] = 0.

    # Realise halo count in each voxel
    Nhalo = np.random.poisson(lam=np.nan_to_num(Nbar))
    return Nhalo

realise_halo_catalogue(self, Nhalo, scatter=False, scatter_type='uniform')

Generate a catalogue of halo positions (comoving) within the 3D box, using a previously-calculated field of halo counts per voxel.

Parameters:

Name Type Description Default
Nhalo array_like

3D array containing the number of halos of this type in each voxel, typically generated by self.halo_count_field().

 required
scatter bool

Whether to randomly shift the positions of the halos within a voxel. If False, halos in the same voxel have the same position.

False
scatter_type str

If scatter=True, how the positions within the voxel should be assigned. Default: 'uniform' (uniform random positions).

 'uniform'

Returns:

Type Description
cat (array_like)

Catalogue of halo positions, of shape (Nhalos, 3) (comoving Cartesian positions x,y,z).
Source code in fastbox/halos.py 
def realise_halo_catalogue(self, Nhalo, scatter=False, scatter_type='uniform'):
    """
    Generate a catalogue of halo positions (comoving) within the 3D box, 
    using a previously-calculated field of halo counts per voxel.

    Parameters:
        Nhalo (array_like):
            3D array containing the number of halos of this type in each voxel, 
            typically generated by `self.halo_count_field()`.

        scatter (bool, optional):
            Whether to randomly shift the positions of the halos within a 
            voxel. If False, halos in the same voxel have the same position. 

        scatter_type (str, optional):
            If `scatter=True`, how the positions within the voxel should be 
            assigned. Default: 'uniform' (uniform random positions).

    Returns:
        cat (array_like):
            Catalogue of halo positions, of shape (Nhalos, 3) (comoving 
            Cartesian positions x,y,z).
    """
    cat_x, cat_y, cat_z = [], [], []

    # Loop over sets of voxels with the same count
    unique_counts = np.unique(Nhalo)
    for count in unique_counts:
        if count < 1.: continue

        # Find voxels containing this number of halos and add to catalogue 
        # (repeated 'count' times)
        idx_x, idx_y, idx_z = np.where(Nhalo == count)
        cat_x.append(np.repeat(idx_x, count))
        cat_y.append(np.repeat(idx_y, count))
        cat_z.append(np.repeat(idx_z, count))

    # Concatenate into a single array
    cat = np.column_stack( (np.concatenate(cat_x),
                            np.concatenate(cat_y),
                            np.concatenate(cat_z)) ).astype(np.float64)

    # Add scatter in positions
    if scatter:
        if scatter_type == 'uniform':
            # Uniform scatter within each voxel
            cat += np.random.uniform(0., 1.-1e-8, cat.size).reshape(cat.shape)
        else:
            raise ValueError("scatter_type='%s' not recognised" 
                              % scatter_type)

    # Scale coordinates with box size
    cat[:,0] *= self.box.Lx / self.box.N
    cat[:,1] *= self.box.Ly / self.box.N
    cat[:,2] *= self.box.Lz / self.box.N

    return cat