Python Documentation

Top-level modules

boxes module

Functions and classes to manipulate N-Dimensional boxes

scampy.boxes.shift_periodic_box(box, shift, Lbox)

Given the cartesian coordinates from a periodic box applies a shift

Parameters:

box (ndarray) – Input comoving cartesian coordinates of particles in the box. Expects an array with with shape (Nparticles, Ndim). It will assume the coordinates are expressed in units of [cMpc/h] and belong to the interval (0,Lbox)
shift (scalar or iterable) – shift to apply along all or each direction. if it is a scalar, the same shift will be applied to all directions. If it is an iterable it should have the same length of Ndim, specifying the shift to apply along each dimension.
Lbox (scalar) – box side lenght in [Mpc/h]

Returns:

Copy of the input array with the coordinate shift applied

Return type:

ndarray

scampy.boxes.centre_box(box, Lbox)

Changes the coordinates of the particles in a comoving cubic box moving the origin to the box centre.

Parameters:

box (ndarray) – Input comoving coordinates of particles in the box. Expects an array with with shape (Nparticles, Ndim). It will assume the coordinates are expressed in units of [Mpc/h] and belong to the interval (0,Lbox)
Lbox (scalar) – box side lenght in [Mpc/h]

Returns:

Copy of the input array with the coordinate centred

Return type:

ndarray

scampy.boxes.cartesian_to_polar(coords, centre=(0.0, 0.0, 0.0))

Converts cartesian coordinates (x,y,z) to polar coordinates (rho, theta, phi) with physics convention.

Parameters:

coords (iterable) – an array or list or tuple with shape (3,Npoints), where Npoints is the number of coordinates to convert.
centre (iterable) – (Optional, default = (0,0,0)) an array or list or tuple with shape (3,) containing the coordinates of the centre of the polar coordinates system in cartesian coordinates

Returns:

spherical coordinates of the input positions

Return type:

iterable

See also

polar_to_cartesian: inverse transform

scampy.boxes.polar_to_cartesian(coords, centre=(0.0, 0.0, 0.0))

Converts polar coordinates (rho, theta, phi) to cartesian coordinates (x,y,z) with physics convention.

Parameters:

coords (iterable) – an array or list or tuple with shape (3,Npoints), where Npoints is the number of coordinates to convert.
centre (iterable) – (Optional, default = (0,0,0)) an array or list or tuple with shape (3,) containing the coordinates of the centre of the polar coordinates system in cartesian coordinates

Returns:

cartesian coordinates of the input positions

Return type:

iterable

See also

cartesian_to_polar: inverse transform

scampy.boxes.cartesian_to_equatorial(coords, centre=(0.0, 0.0, 0.0))

Converts cartesian coordinates (x,y,z) to astronomical coordinates (rho, theta, phi) with equatorial convention.

Parameters:

coords (iterable) – an array or list or tuple with shape (3,Npoints), where Npoints is the number of coordinates to convert.
centre (iterable) – (Optional, default = (0,0,0)) an array or list or tuple with shape (3,) containing the coordinates of the centre of the equatorial coordinates system (the observer’s position) in cartesian coordinates

Returns:

equatorial coordinates of the input positions

Return type:

iterable

See also

equatorial_to_cartesian: inverse transform

scampy.boxes.equatorial_to_cartesian(coords, centre=(0.0, 0.0, 0.0))

Converts astronomical coordinates (d, ra, dec) to cartesian coordinates (x,y,z) with equatorial convention.

Parameters:

coords (iterable) – an array or list or tuple with shape (3,Npoints), where Npoints is the number of coordinates to convert.
centre (iterable) – (Optional, default = (0,0,0)) an array or list or tuple with shape (3,) containing the coordinates of the centre of the equatorial coordinates system (the observer’s position) in cartesian coordinates

Returns:

cartesian coordinates of the input positions

Return type:

iterable

See also

cartesian_to_equatorial: inverse transform

scampy.boxes.angular_to_euclidean_dist(d, r=1.0)

Converts an angular distance to an euclidean distance, projected on a sphere of given radius.

Parameters:

d (float or array-like of floats) – angular distances to convert
r (float or array-like of floats) – (Optional, default=1.0) radius of the sphere. If an array is given, it should be broadcastable to the shape of d

Returns:

the projected euclidean distances

Return type:

float or ndarray

scampy.boxes.random_projection(RA_lim, Dec_lim, size=1, rng=None)

Produce random equatorial coordinates, projected in the unit-sphere

Parameters:

RA_lim (tuple)
Dec_lim (tuple)
size (int)
rng (random number generator or int) – (Optional, default = None) a random number generator with a uniform function.

Returns:

RA (float or ndarray) – random values of Right-Ascension (in radians). if size=1 this will be a float
Dec (float or ndarray) – random values of Declination (in radians). if size=1 this will be a float

catalogue module

Handling catalogues of haloes and subhaloes

class scampy.catalogue.kernelCat(*args, **kwargs)

Base catalogue container enforcing mandatory coordinate fields.

Extends fixedSizeDict by requiring that the fields listed in the class attribute fields ('X', 'Y', 'Z') are always present. Subclasses extend this list with additional mandatory fields.

Parameters:

*args – Forwarded to fixedSizeDict. Must include at least the mandatory fields in fields.
**kwargs – Forwarded to fixedSizeDict. Must include at least the mandatory fields in fields.

coord(): Return the 3D Cartesian coordinates as a (Nobj, 3) array.

return_sample(X_fields=[], Y_fields=[])

Return two arrays built from selected fields.

Parameters:

X_fields (list of str, optional) – Field names to stack into the first output array.
Y_fields (list of str, optional) – Field names to stack into the second output array.

Returns:

X (ndarray, shape (Nobj, len(X_fields)))
Y (ndarray, shape (Nobj, len(Y_fields)))

class scampy.catalogue.haloCat(*args, **kwargs)

Catalogue of host haloes.

Extends kernelCat with the mandatory halo-specific fields:

Mhalo — halo mass \([M_\odot\,h^{-1}]\)
Rhalo — halo radius \([h^{-1}\,\mathrm{Mpc}]\)
firstSub — index of the first sub-halo in the sub-halo table
numSubs — total number of sub-haloes in the halo

centrals()

Return list of indices in the sub-halo table locating the central sub-halo of each halo. (takes no arguments)

See also

satellites: same for satellite sub-haloes

satellites()

Return list of indices in the sub-halo table locating the satellite sub-haloes of each halo. (takes no arguments)

Warning

Computes the slices only for un-masked haloes

See also

centrals: same for central sub-haloes

class scampy.catalogue.subhaloCat(*args, **kwargs)

Catalogue of sub-haloes.

Extends kernelCat with the mandatory sub-halo-specific fields:

Msubh — sub-halo mass \([M_\odot\,h^{-1}]\)
Parent — index of the parent halo in the haloCat table

Nsub(mask=None)

Function returning the number of un-masked sub-haloes in not-empty parent halo A halo is considered ‘not-empty’ when it has at least one sub-halo which is un-masked

Parameters:

mask (bool or iterable or None) – It accepts an iterable with lenght equal to the size of the sub-halo catalogue or ‘None’ to use all the objects in the sub-halo catalogue.

Returns:

Nsub (numpy array) – array with the number of sub-haloes in each not-empty halo
hidx (numpy array) – indices in the halo catalogue of the not-empty haloes

See also

catalogue.Ncen: function returning the valid central galaxy number
catalogue.Nsat: function returning the valid satellite galaxies number
catalogue.NcenNsat: function returning a tuple with the number of valid galaxies

class scampy.catalogue.catalogue(haloes, subhaloes, Lbox, **kwargs)

class catalogue for managing halo/subhalo hierarchies

Parameters:

haloes (haloCat object)
subhaloes (subhaloCat object)
**kwargs –
any argument passed as keyword argument is considered as a metadatum of the catalogue (as long as it does not already exist in the internal

dictionary of the class).

save(outPath, hard=False)

Stores current state of catalogue in an hdf5 file

Parameters:

outPath (string) – ‘/path/to/name/of/file’ the function will append the ‘.hdf5’ extension to this path then check whether the file already exists
hard (bool) – if True eventually overwrites an already existing file with the same name (default = False)

See also

load: function for building catalogue from hdf5 file

Nsub(mask_subhaloes=None)

Function returning the number of un-masked sub-haloes in not-empty parent halo A halo is considered ‘not-empty’ when it has at least one sub-halo which is un-masked

Parameters:

mask_subhaloes (bool or iterable or None) – It accepts an iterable with lenght equal to the size of the sub-halo catalogue or ‘None’ to use all the objects in the sub-halo catalogue.

Returns:

Nsub (numpy array) – array with the number of sub-haloes in each not-empty halo
hidx (numpy array) – indices in the halo catalogue of the not-empty haloes

See also

Ncen: function returning the valid central galaxy number
Nsat: function returning the valid satellite galaxies number
NcenNsat: function returning a tuple with the number of valid galaxies

Ncen(mask=None)

Function returning the number of un-masked central galaxies

Parameters:: mask (bool or iterable or None) – Accepts an iterable with lenght equal to the lenght of the sub-halo catalogue mask or ‘None’ to use all the objects in the sub-halo catalogue.
Returns:: array of shape (haloes.size(),) with the number of centrals for each halo
Return type:: numpy array

See also

Nsub: function returning the number of un-masked sub-haloes in not-empty haloes
Nsat: function returning the valid satellite galaxies number
NcenNsat: function returning a tuple with the number of valid galaxies

Nsat(mask=None)

Function returning the number of un-masked satellite galaxies

Parameters:: add_mask (bool or iterable or None) – It accepts either a boolean (True = mask all the entries, False = unmask all the entries), an iterable with lenght equal to the lenght of the sub-halo catalogue mask or ‘None’ to only use the default mask of the sub-halo catalogue.
Returns:: array of shape (haloes.size(),) with the number of satellites for each halo
Return type:: numpy array

See also

Nsub: function returning the number of un-masked sub-haloes in not-empty haloes
Ncen: function returning the valid central galaxy number
NcenNsat: function returning a tuple with the number of valid galaxies

NcenNsat(mask=None)

Function returning a tuple with the number of un-masked galaxies (centrals and satellites)

Parameters:

add_mask (bool or iterable or None) – It accepts either a boolean (True = mask all the entries, False = unmask all the entries), an iterable with lenght equal to the lenght of the sub-halo catalogue mask or ‘None’ to only use the default mask of the sub-halo catalogue.

Returns:

Nc (numpy array) – array of shape (haloes.size(),) with the number of centrals for each halo
Ns (numpy array) – array of shape (haloes.size(),) with the number of satellites for each halo

See also

Nsub: function returning the number of un-masked sub-haloes in not-empty haloes
Ncen: function returning the valid central galaxy number
Nsat: function returning the valid satellite galaxies number

centrals(hmask=None, smask=None)

Return indices in the sub-halo table locating the central sub-halo of each halo.

Parameters:

hmask (ndarray of bool or None, optional) – Boolean mask of length haloes.size; only centrals whose parent halo passes this mask are returned. If None (default) all haloes are included.
smask (ndarray of bool or None, optional) – Boolean mask of length subhaloes.size; only sub-haloes passing this mask are considered as centrals. If None (default) all sub-haloes are included.

Returns:

cen – Indices into the sub-halo catalogue of the selected central sub-haloes.

Return type:

ndarray of int

See also

satellites: same for satellite sub-haloes

satellites(hmask=None, smask=None)

Return indices in the sub-halo table locating the satellite sub-haloes of each halo.

Parameters:

hmask (ndarray of bool or None, optional) – Boolean mask of length haloes.size; only satellites whose parent halo passes this mask are returned. If None (default) all haloes are included.
smask (ndarray of bool or None, optional) – Boolean mask of length subhaloes.size; only sub-haloes passing this mask are considered as satellites. If None (default) all sub-haloes are included.

Returns:

sat – Indices into the sub-halo catalogue of the selected satellite sub-haloes.

Return type:

ndarray of int

Warning

Computes the slices only for un-masked haloes.

See also

centrals: same for central sub-haloes

hod module

Halo Occupation Distribution (HOD) models.

Provides classes for the 5-parameter HOD prescription of Zheng, Coil & Zehavi (2007) and Zheng et al. (2009), as described in Sec. 2.1 and Eqs. 19–20 of Ronconi et al. (2020), together with a redshift-dependent extension. All classes expose Pcen(M) and Psat(M) methods compatible with halo_model.

scampy.hod.get_hosts(cat, Pcen, Psat, method='binomial', kw_pcen={}, kw_psat={}, rng=None, kw_rng={})

Select host sub-haloes from a catalogue using arbitrary occupation functions.

Low-level helper that applies user-supplied central and satellite occupation probabilities to a catalogue without requiring a full HOD class instance. For the standard HOD workflow prefer the get_hosts() method of HOD.

Parameters:

cat (scampy.catalogue.catalogue) – Halo/sub-halo catalogue to populate.
Pcen (callable or array-like) – Central occupation probability. If callable it is called as Pcen(Mhalo, **kw_pcen) and must return an array of probabilities in [0, 1]. If array-like it must have length cat.haloes.size.
Psat (callable or array-like) – Desired mean number of satellites per halo. If callable it is called as Psat(Mhalo, **kw_psat, **kw_pcen) and must return an array of non-negative values. If array-like it must have length cat.haloes.size.
kw_pcen (dict, optional) – Keyword arguments forwarded to Pcen when callable.
kw_psat (dict, optional) – Keyword arguments forwarded to Psat when callable.
rng (numpy.random.Generator or None, optional) – Random number generator. If None (default) a new generator is created via numpy.random.default_rng(**kw_rng).
kw_rng (dict, optional) – Keyword arguments forwarded to numpy.random.default_rng when rng is None.

Returns:

cen_gxy (ndarray of bool) – Boolean mask over cat.subhaloes selecting central galaxies.
sat_gxy (ndarray of bool) – Boolean mask over cat.subhaloes selecting satellite galaxies.

class scampy.hod.HOD(Mmin=10000000000.0, sigma=1.0, M0=10000000000.0, M1=10000000000.0, alpha=1.0)

5-parameter Halo Occupation Distribution (Zheng et al. 2007).

Computes the average number of central and satellite galaxies hosted by a halo of mass \(M_h\) (Eqs. 19–20 of Ronconi et al. 2020):

\[\langle N_\mathrm{cen}(M_h)\rangle = \frac{1}{2}\left[1 + \mathrm{erf}\! \left(\frac{\log M_h - \log M_\mathrm{min}}{\sigma} \right)\right],\]

\[\langle N_\mathrm{sat}(M_h)\rangle = \langle N_\mathrm{cen}(M_h)\rangle \left(\frac{M_h - M_0}{M_1}\right)^\alpha.\]

The satellite term is conditioned on the central (i.e. a halo must host a central to host satellites). For the unconditioned variant see HOD_unconditioned_sat.

Parameters:

Mmin (float, optional) – Characteristic minimum halo mass \(M_\mathrm{min}\) for hosting a central galaxy (default: 1e10).
sigma (float, optional) – Width \(\sigma\) of the central occupation transition in \(\log_{10} M\) (default: 1.0).
M0 (float, optional) – Satellite mass offset \(M_0\) (default: 1e10).
M1 (float, optional) – Satellite mass normalisation \(M_1\) (default: 1e10).
alpha (float, optional) – Satellite power-law slope \(\alpha\) (default: 1.0).

Pcen(Mh, **kwargs)

Average number of central galaxies per halo of mass \(M_h\).

\[\langle N_\mathrm{cen}(M_h)\rangle = \frac{1}{2}\left[1 + \mathrm{erf}\! \left(\frac{\log M_h - \log M_\mathrm{min}}{\sigma} \right)\right].\]

Parameters:: Mh (scalar or array-like) – Halo mass \(M_h\) in \([M_\odot\,h^{-1}]\).
Keyword Arguments:: alpha (Mmin, sigma, M0, M1,) – Override the instance parameters for this call.
Returns:: Values in \([0, 1]\), same shape as Mh.
Return type:: ndarray

Psat(Mh, **kwargs)

Average number of satellite galaxies per halo of mass \(M_h\).

\[\langle N_\mathrm{sat}(M_h)\rangle = \langle N_\mathrm{cen}(M_h)\rangle \left(\frac{M_h - M_0}{M_1}\right)^\alpha.\]

Returns zero when \(M_h \leq M_0\).

Parameters:: Mh (scalar or array-like) – Halo mass \(M_h\) in \([M_\odot\,h^{-1}]\).
Keyword Arguments:: alpha (Mmin, sigma, M0, M1,) – Override the instance parameters for this call.
Returns:: Non-negative values, same shape as Mh.
Return type:: ndarray

get_hosts(cat, smask=None, method='binomial', rng=None, kw_rng={}, **kwargs)

Apply the HOD to select host sub-haloes from a catalogue.

Parameters:

cat (scampy.catalogue.catalogue) – Halo/sub-halo catalogue to populate.
smask (ndarray of bool or None, optional) – Boolean mask of length cat.subhaloes.size pre-selecting sub-halo candidates. True = include. If None (default) all sub-haloes are considered.
method ({'binomial', 'rankorder'}, optional) –
Galaxy assignment method (default: 'binomial'):
- 'binomial' — centrals drawn from a Binomial(1, Pcen) distribution; satellites drawn from a Binomial(1, Psat/Nsh) distribution.
- 'rankorder' — centrals as above; satellites chosen by keeping the first \(N_\mathrm{sat}\) sub-haloes in the input order (assumes sub-haloes are sorted by some property, e.g. mass).
rng (numpy.random.Generator or None, optional) – Random number generator. If None (default) one is created via numpy.random.default_rng(**kw_rng).
kw_rng (dict, optional) – Keyword arguments forwarded to numpy.random.default_rng (default: {}).

Keyword Arguments:

alpha (Mmin, sigma, M0, M1,) – Override the instance HOD parameters for this call.

Returns:

cen_gxy (ndarray of bool) – Boolean mask of length cat.subhaloes.size; True where a central galaxy is assigned.
sat_gxy (ndarray of bool) – Boolean mask of length cat.subhaloes.size; True where a satellite galaxy is assigned.

class scampy.hod.HOD_unconditioned_sat(Mmin=10000000000.0, sigma=1.0, M0=10000000000.0, M1=10000000000.0, alpha=1.0)

5-parameter HOD with unconditioned satellite occupation.

Identical to HOD except that the satellite term is not conditioned on the central:

\[\langle N_\mathrm{cen}(M_h)\rangle = \frac{1}{2}\left[1 + \mathrm{erf}\! \left(\frac{\log M_h - \log M_\mathrm{min}}{\sigma} \right)\right],\]

\[\langle N_\mathrm{sat}(M_h)\rangle = \left(\frac{M_h - M_0}{M_1}\right)^\alpha.\]

Parameters:

Mmin (float, optional) – Characteristic minimum halo mass \(M_\mathrm{min}\) (default: 1e10).
sigma (float, optional) – Width of the central occupation transition in \(\log_{10} M\) (default: 1.0).
M0 (float, optional) – Satellite mass offset \(M_0\) (default: 1e10).
M1 (float, optional) – Satellite mass normalisation \(M_1\) (default: 1e10).
alpha (float, optional) – Satellite power-law slope \(\alpha\) (default: 1.0).

Pcen(Mh, **kwargs)

Average number of central galaxies per halo of mass \(M_h\).

Same formula as HOD.Pcen().

Parameters:: Mh (scalar or array-like) – Halo mass in \([M_\odot\,h^{-1}]\).
Keyword Arguments:: alpha (Mmin, sigma, M0, M1,) – Override the instance parameters for this call.
Returns:: Values in \([0, 1]\), same shape as Mh.
Return type:: ndarray

Psat(Mh, **kwargs)

Average number of satellite galaxies per halo of mass \(M_h\).

\[\langle N_\mathrm{sat}(M_h)\rangle = \left(\frac{M_h - M_0}{M_1}\right)^\alpha.\]

Unlike HOD.Psat(), this is not conditioned on the central occupation. Returns zero when \(M_h \leq M_0\).

Parameters:: Mh (scalar or array-like) – Halo mass in \([M_\odot\,h^{-1}]\).
Keyword Arguments:: alpha (Mmin, sigma, M0, M1,) – Override the instance parameters for this call.
Returns:: Non-negative values, same shape as Mh.
Return type:: ndarray

get_hosts(cat, method='binomial', rng=None, kw_rng={}, **kwargs)

Apply the HOD to select host sub-haloes from a catalogue.

Parameters:

cat (scampy.catalogue.catalogue) – Halo/sub-halo catalogue to populate.
method ({'binomial', 'rankorder'}, optional) –
Galaxy assignment method (default: 'binomial'):
- 'binomial' — centrals drawn from a Binomial(1, Pcen) distribution; satellites drawn from a Binomial(1, Psat/Nsh) distribution.
- 'rankorder' — centrals as above; satellites chosen by keeping the first \(N_\mathrm{sat}\) sub-haloes in the input order (assumes sub-haloes are sorted by some property, e.g. mass).
rng (numpy.random.Generator or None, optional) – Random number generator. If None (default) one is created via numpy.random.default_rng(**kw_rng).
kw_rng (dict, optional) – Keyword arguments forwarded to numpy.random.default_rng (default: {}).

Keyword Arguments:

alpha (Mmin, sigma, M0, M1,) – Override the instance HOD parameters for this call.

Returns:

cen_gxy (ndarray of bool) – Boolean mask of length cat.subhaloes.size; True where a central galaxy is assigned.
sat_gxy (ndarray of bool) – Boolean mask of length cat.subhaloes.size; True where a satellite galaxy is assigned.

class scampy.hod.HOD_zdep(redshift, lMmin=10, lMmin_b=0.0, lsigma=0.0, lsigma_b=0.0, lM0=10, lM0_b=0.0, lM1=10, lM1_b=0.0, lalpha=0.0, lalpha_b=0.0)

Redshift-dependent 5-parameter HOD.

Extends HOD by allowing each parameter to evolve linearly with redshift in log-space:

\[\log_{10} X(z) = \ell_X + \ell_{X,b}\,z,\]

so that \(X(z) = 10^{\ell_X + \ell_{X,b}\,z}\). Setting all \(\ell_{X,b} = 0\) recovers a redshift-independent HOD.

The occupation functions follow the same HOD formulae with the redshift-evaluated parameters.

Parameters:

redshift (float or array-like) – Redshift(s) at which to evaluate the HOD parameters.
lMmin (float, optional) – \(\log_{10} M_\mathrm{min}\) at \(z=0\) (default: 10).
lMmin_b (float, optional) – Redshift slope of \(\log_{10} M_\mathrm{min}\) (default: 0.0).
lsigma (float, optional) – \(\log_{10} \sigma\) at \(z=0\) (default: 0.0).
lsigma_b (float, optional) – Redshift slope of \(\log_{10} \sigma\) (default: 0.0).
lM0 (float, optional) – \(\log_{10} M_0\) at \(z=0\) (default: 10).
lM0_b (float, optional) – Redshift slope of \(\log_{10} M_0\) (default: 0.0).
lM1 (float, optional) – \(\log_{10} M_1\) at \(z=0\) (default: 10).
lM1_b (float, optional) – Redshift slope of \(\log_{10} M_1\) (default: 0.0).
lalpha (float, optional) – \(\log_{10} \alpha\) at \(z=0\) (default: 0.0).
lalpha_b (float, optional) – Redshift slope of \(\log_{10} \alpha\) (default: 0.0).

set(redshift=None, **kwargs)

Update redshift and/or free parameters, then recompute self.params.

Parameters:

redshift (float or array-like or None, optional) – New redshift(s). If None the current redshift is kept.
**kwargs – Any subset of the free parameters (lMmin, lMmin_b, lsigma, lsigma_b, lM0, lM0_b, lM1, lM1_b, lalpha, lalpha_b). Unknown keys are silently ignored.

lightcone module

Functions and classes to generate and manipulate lightcones

scampy.lightcone.sequential_lightcones_limits(zmax, Lbox, d2z=None, centre=True)

Computes the central redshift of boxes of given size necessary to build a lightcone from the resulting snapshots up to some given redshift.

Parameters:

zmax (float) – Maximum redshift of the lightcone
Lbox (float) – comoving side size of the simulation box
d2z (function or dictionary or None) – (Optional, default = 'None') either a function taking a comoving distance as argument and returning the corresponding redshift or a dictionary with fields 'cosmo' with an instance of object of class scampy.cosmology.model and 'zgrid' the redshift grid on witch to interpolate the necessary function d2z (the thinner the grid, the more accurate the interpolation)
centre (bool) – (Optional, default = True) if true associate the values computed to the centre of the box, otherwise, they will be associated to the closest side.

Returns:

zeta (1d-array) – array with the sequential redshifts of boxes necessary to build a lightcone
dbox (1d-array) – array with the sequential comoving distances of the boxes necessary to build a lightcone

scampy.lightcone.ang_size(lC, z, cosmo, rad=True)

Converts a comoving distance, perpendicular to the LoS, at a given redshift to an angle.

Parameters:

lC (scalar or 1d-array) – Comoving distance in [Mpc/h], assumed perpendicular to the LoS. Note that if this is a 1D-array, then argument `z` must be a scalar.
z (scalar or 1d-array) – redshift along the LoS. Note that if this is a 1D-array, then argument `lC` must be a scalar.
cosmo (scampy.cosmology.model) – Instance of an object of type scampy.cosmology.model
rad (bool) – (Optional, default = True) whether to output an angle in radians (rad = True), or degrees (rad = False)

Returns:

Angle corresponding to the input distances at the input redshifts

Return type:

scalar or nd-array

Warning

Does not support broadcasting on two dimensions. Only one among lC and z can be a 1d-array, not both.

scampy.lightcone.size_ang(theta, z, cosmo, rad=True)

Converts an angle at a given redshift to a comoving distance, perpendicular to the LoS.

Parameters:

theta (scalar or 1d-array) – Viewing angle. Note that if this is a 1D-array, then argument `z` must be a scalar.
z (scalar or 1d-array) – redshift along the LoS. Note that if this is a 1D-array, then argument `theta` must be a scalar.
cosmo (scampy.cosmology.model) – Instance of an object of type scampy.cosmology.model
rad (bool) – (Optional, default = True) whether the input angle is given in radians (rad = True) or degrees (rad = False)

Returns:

Comoving distances in [Mpc/h], assumed perpendicular to the LoS, computed at the input redshifts.

Return type:

scalar or nd-array

Warning

Does not support broadcasting on two dimensions. Only one among theta and z can be a 1d-array, not both.

scampy.lightcone.to_cone(box, zbox, Lbox, cosmo, funcs=None, theta=None, zmin=None, zmax=None, return_indices=False)

Converts the coordinates of particles in a comoving cubic box to coordinates in the lightcone.

Paramters

box3d-array: Input comoving coordinates of particles in the box. Expects an array with with shape (Nparticles, 3). It will assume the coordinates are expressed in units of [Mpc/h] and belong to the interval (0,Lbox)
zboxscalar: Redshift associated to the centre of the box
Lboxscalar: box side lenght in [Mpc/h]
cosmoscampy.cosmology.model: Instance of an object of type scampy.cosmology.model
funcstuple: (Optional, default = None) tuple (function,inverse) computing comoving distance from redshift and vice-versa. The functions should be vectorizable and returning scalar of numpy.ndarray. If the tuple is not provided attribute `cosmo` will be used to allocate a couple of interpolators serving this purpose.
thetascalar: (Optional, default = None) viewing angle in radians. If None is passed (default behaviour), the largest possible viewing angle associated with the provided Lbox and zbox is assumed.
zminscalar or None: (Optional, default = None) if passed, only selects coordinates with redshift larger than this value (exclusive)
zmaxscalar or None: (Optional, default = None) if passed, only selects coordinates with redshift smaller than this value (inclusive)
return_indicesbool: (Optional, default = False) when True, also return the mask that, applied to the input box argument, selects only coordinates within the light-cone

returns:

cone_coords (3d-array) – Spherical coordinates of the generated lightcone with shape (Nview, 3), where Nview <= Nparticles is the number of particles within the field of view. The light-cone provides, for each particle within the field of view, a (radians, radians, redshift) array.
cone_indices (bool 1d-array) – the boolean mask that, applied to the input box argument, selects only coordinates within the light cone. Only provided if return_indices is True.

Warning

This will just compute the radial coordinates of objects WITHOUT the effect of gravitational lensing, no light rays will be computed.

scampy.lightcone.collate_lightcones(zlist, zmin, zmax, Lbox, initshift=0.0, cosmo=None)

Computes automathically the tuple of quantities necessary for building a lightcone by collating coordinates from different snapshots of a cosmological N-body simulation. By shifting the position of the observer it is possible to obtain different realizations.

Parameters:

zlist (iterable) – list of redshifts associated to each box
zmin (scalar float) – minimum redshift of final lightcone
zmax (scalar float) – maximum redshift of final lightcone
Lbox (scalar float) – size of the periodic box
initshift (scalar or iterable) – (Optional, default = 0.0) either a scalar float or an iterable that can be unpacked in 3 values. In case of a float, it will be replicated along the 3 axis, shifting the periodic box by initshift Mpc/h in each dimension. In case iterable, it should strictly have length = 3, the first two elements are the shifts associated with the perpendicular direction while the 3rd is the shift applied along the LoS direction
cosmo (scampy.cosmology.model) – (Optional, default = None) Instance of an object of type scampy.cosmology.model. if none is passed, it will be built with default cosmological parameters (not an optimal solution)

Returns:

zcen (1d-array) – list of the central redshifts for each box (this differs from the input zlist iterable only by the first and last value, that are computed in order to have zmin < zcen[0] and zcen[-1] < zmax )
dcen (1d-array) – comoving distance of the observer from the centre of each box
shifts (2d-array) – shifts that are to be applied to each direction of each box in order to preserve continuity of the simulation (shape=(Nbox,3))
zlim (2d-array) – redshift limits of each box composing the lightcone (shape=(Nbox,2)). For each box j, these should be zlim[j,0] < zcen[j] < zlim[j,1]

scampy.lightcone.box_to_cone(box, Lbox, cosmo, central_redshift, redshift_limits, observer=(0.0, 0.0, 0.0), funcs=None, FoV=None, return_indices=False)

Converts the coordinates of particles in a comoving cubic box to coordinates in the lightcone.

Paramters

box3d-array: Input comoving coordinates of particles in the box. Expects an array with with shape (Nparticles, 3). It will assume the coordinates are expressed in units of [Mpc/h] and belong to the interval (0,Lbox)
Lboxscalar: box side lenght in [Mpc/h]
cosmoscampy.cosmology.model: Instance of an object of type scampy.cosmology.model
central_redshiftscalar: Redshift associated to the centre of the box
redshift_limitsiterable: an iterable with two values for the minimum and maximum redshift to include in the light-cone
funcstuple: (Optional, default = None) tuple (function,inverse) computing comoving distance from redshift and vice-versa. The functions should be vectorizable and returning scalar of numpy.ndarray. If the tuple is not provided attribute `cosmo` will be used to allocate a couple of interpolators serving this purpose.
FoVscalar: (Optional, default = None) viewing angle in radians. If None is passed (default behaviour), the largest possible viewing angle associated with the provided Lbox and central redshift is assumed.
return_indicesbool: (Optional, default = False) when True, also return the mask that, applied to the input box argument, selects only coordinates within the light-cone

returns:

cone_coords (3d-array) – Spherical coordinates of the generated lightcone with shape (Nview, 3), where Nview <= Nparticles is the number of particles within the field of view. The light-cone provides, for each particle within the field of view, a (radians, radians, redshift) array.
cone_indices (bool 1d-array) – the boolean mask that, applied to the input box argument, selects only coordinates within the light cone. Only provided if return_indices is True.

Warning

This will just compute the radial coordinates of objects WITHOUT the effect of gravitational lensing, no light rays will be computed.

power_spectrum module

Linear matter power spectrum and related statistics.

Provides the power_spectrum class, which wraps a tabulated linear \(P(k)\) at \(z=0\), normalises it to the cosmological \(\sigma_8\), evolves it to arbitrary redshift via the linear growth factor, and exposes derived quantities used by the halo model: the variance \(\sigma^2(R,z)\), its mass-keyed counterpart \(\sigma^2(M,z)\), and the real-space two-point correlation function \(\xi(r,z)\).

class scampy.power_spectrum.power_spectrum(kh, pk0, kmin=0.0001, kmax=100, cosmo=None, thinness=1000)

Container for the linear matter power spectrum and its derived statistics.

On construction the input tabulated \(P(k)\) is interpolated onto a uniform log-spaced grid, normalised so that \(\sigma_8\) matches the value stored in the cosmological model, and the linear growth factor \(D(z)/D(0)\) is pre-computed on a redshift grid for fast evaluation.

Parameters:

kh (array-like) – Wavenumbers of the input power spectrum in \([h\,\mathrm{Mpc}^{-1}]\).
pk0 (array-like) – Linear power spectrum at \(z=0\) in \([h^{-3}\,\mathrm{Mpc}^3]\), evaluated at the wavenumbers kh.
kmin (float, optional) – Minimum wavenumber of the internal interpolation grid (default: 1e-4).
kmax (float, optional) – Maximum wavenumber of the internal interpolation grid (default: 100).
cosmo (scampy.cosmology.model or None, optional) – Cosmological model instance. If None a default model is constructed.
thinness (int, optional) – Number of redshift points used to pre-compute the growth factor (default: 1000).

compute_sigma8()

Compute and store \(\sigma_8\) from the input power spectrum.

Evaluates

\[\sigma_8^2 = \frac{1}{2\pi^2} \int_0^\infty k^2\,P(k)\, \tilde{W}^2(8\,h^{-1}\,\mathrm{Mpc}\cdot k)\, \mathrm{d}k,\]

where \(\tilde{W}\) is the Fourier transform of a top-hat window of radius \(8\,h^{-1}\,\mathrm{Mpc}\). The result is stored as self.sigma8 and a correction factor self.sigma8correction = cosmo.sigma8 / sigma8 is set so that subsequent calls to Pz() return correctly normalised spectra.

Returns:: The unnormalised \(\sigma_8\) computed from the raw input pk0.
Return type:: float

Pz(kk, zz=0.0, comoving=True)

Linear matter power spectrum evolved to redshift \(z\).

\[P(k, z) = \sigma_8^{\mathrm{cosmo},2}\, \frac{D^2(z)}{D^2(0)}\, \frac{P_0(k)}{\sigma_8^2},\]

where \(P_0(k)\) is the input spectrum at \(z=0\) and \(D(z)\) is the linear growth factor.

Parameters:

kk (scalar or array-like) – Wavenumbers in \([h\,\mathrm{Mpc}^{-1}]\) (comoving) or \([\mathrm{Mpc}^{-1}]\) (physical).
zz (scalar or array-like, optional) – Redshift(s) (default: 0.0).
comoving (bool, optional) – If True (default) wavenumbers and the output are in comoving units; if False a factor of \(h\) is applied.

Returns:

Power spectrum values with shape (kk.size,) or (kk.size, zz.size).

Return type:

ndarray

sigma2R(rr, zz, comoving=True)

Variance of the linear density field smoothed on scale \(R\).

\[\sigma^2(R, z) = \frac{1}{2\pi^2} \int_0^\infty k^2\,P(k,z)\, \tilde{W}^2(kR)\,\mathrm{d}k,\]

where \(\tilde{W}(x) = 3(\sin x - x\cos x)/x^3\) is the Fourier transform of a spherical top-hat of radius \(R\).

Parameters:

rr (scalar or array-like) – Smoothing radii in \([h^{-1}\,\mathrm{Mpc}]\).
zz (scalar or array-like) – Redshift(s). The output is broadcast over (rr, zz).
comoving (bool, optional) – Passed to Pz() (default: True).

Returns:

Variance with shape (rr.size,) or (rr.size, zz.size).

Return type:

ndarray

dsigma2RdR(rr, zz, comoving=True)

Derivative of \(\sigma^2(R,z)\) with respect to \(R\).

\[\frac{\mathrm{d}\sigma^2}{\mathrm{d}R}(R, z) = \frac{1}{\pi^2} \int_0^\infty k^3\,P(k,z)\, \tilde{W}(kR)\,\tilde{W}'(kR)\,\mathrm{d}k,\]

where \(\tilde{W}'\) denotes the derivative of the top-hat window with respect to its argument. Used by sigma2M() and dsigma2MdM() via the chain rule.

Parameters:

rr (scalar or array-like) – Smoothing radii in \([h^{-1}\,\mathrm{Mpc}]\).
zz (scalar or array-like) – Redshift(s). The output is broadcast over (rr, zz).
comoving (bool, optional) – Passed to Pz() (default: True).

Returns:

Derivative with shape (rr.size,) or (rr.size, zz.size).

Return type:

ndarray

sigma2M(mm, zz, comoving=True)

Variance of the density field on the mass scale \(M\).

Converts a halo mass to the radius of an equivalent sphere of mean matter density,

\[R(M) = \left(\frac{3M}{4\pi\rho_0}\right)^{1/3},\]

and returns \(\sigma^2(R(M), z)\) via sigma2R().

Parameters:

mm (scalar or array-like) – Halo masses in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s). The output is broadcast over (mm, zz).
comoving (bool, optional) – If True (default) the mean density \(\rho_0\) is the comoving matter density; otherwise the physical value is used.

Returns:

Variance with shape (mm.size,) or (mm.size, zz.size).

Return type:

ndarray

dsigma2MdM(mm, zz, comoving=True)

Derivative of \(\sigma^2(M,z)\) with respect to \(M\).

Applies the chain rule through the \(R(M)\) relation:

\[\frac{\mathrm{d}\sigma^2}{\mathrm{d}M} = \frac{R}{3M}\,\frac{\mathrm{d}\sigma^2}{\mathrm{d}R}.\]

Used internally by halo mass function routines.

Parameters:

mm (scalar or array-like) – Halo masses in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s). The output is broadcast over (mm, zz).
comoving (bool, optional) – If True (default) uses the comoving matter density.

Returns:

Derivative with shape (mm.size,) or (mm.size, zz.size).

Return type:

ndarray

Xi(rr, zz=0.0, comoving=True)

Linear two-point correlation function.

Obtained by Fourier-transforming the linear power spectrum (Eq. 6 of Ronconi et al. 2020):

\[\xi(r, z) = \frac{1}{2\pi^2} \int_0^\infty k^2\,P(k,z)\, \frac{\sin(kr)}{kr}\,\mathrm{d}k.\]

The integral is evaluated via a Fast Spherical Transform.

Parameters:

rr (scalar or array-like) – Comoving separations in \([h^{-1}\,\mathrm{Mpc}]\).
zz (float, optional) – Redshift (default: 0.0). Broadcasting over zz is not yet supported.
comoving (bool, optional) – Passed to Pz() (default: True).

Returns:

Correlation function \(\xi(r,z)\) evaluated at rr.

Return type:

ndarray

sham module

Sub-Halo Abundance Matching (SHAM) utilities.

Provides functions implementing the abundance matching algorithm described in Sec. 3.1.2 of Ronconi et al. (2020): given a population of sub-haloes and an observed distribution of some observable property (e.g. UV luminosity), each sub-halo is assigned a value of that property by matching the two cumulative distributions.

scampy.sham.scatter_array(arr, method='normal', log=False, rng=None, kw_rng={'seed': 555}, **kwargs)

Given an input array adds scattering to the input values. Different methods are available:

‘normal’ : Gaussian scatter around the input array

Parameters:

arr (iterable)
method (str) – (Optional, default = 'normal')
log (bool) – (Optional, default = False)
rng (random number generator) – (Optional, default = None) if None is passed, it will build one using function numpy.random.default_rng(**kw_rng)
kw_rng (dict) – (Optional, default = { ‘seed’ : 555 }) keyword arguments to pass to the default random number generator. (if rng is not None, this input is ignored)

Keyword Arguments:

selected (Depending on the method)
'normal' (* method =) –

scalefloat or iterable
(Optional, default = 1.0) scattering magnitude (in dex). The default value is 0.0

locfloat or iterable
if not zero, adds a bias to the relation. The default value is 0.0.

Returns:

ndarray
a copy of the original array with scatter added

scampy.sham.matching_sample(InSize, OutSize, InProb=None, rng=None, kw_rng={'seed': 555})

Returns a boolean 1D mask of given size, with given number of occurrencies unmasked.

Parameters:

InSize (int) – Required mask size
OutSize (int) – Required number of un-masked occurencies in output mask
InProb (iterable) – (Optional, default = None) An iterable of length = InSize with probability associated to each occurrency in the output mask. If not given, the sample assumes a uniform distribution over all entries in the output mask.
rng (random number generator) – (Optional, default = None) if None is passed, it will build one using function numpy.random.default_rng(**kw_rng)
kw_rng (dict) – (Optional, default = { ‘seed’ : 555 }) keyword arguments to pass to the default random number generator. (if rng is not None, this input is ignored)

Returns:

1d-array
boolean mask with size = InSize and summing to OutSize.

scampy.sham.match_distribution(var, func, minF=-20, maxF=-10, minP=1e-20, maxP=1.0, nbinF=200, factF=1.0, minV=None, maxV=None, nbinV=100, factV=1.0, scatter_type=None, kw_scatter={}, rng=None, kw_rng={'seed': 555})

Implements the Abundance Matching algorithm. A value in the distribution func is associated to each element of the array var. Note that, the underline assumption is that all the elements of var can be linked to an element in the distribution func. This will not perform an HOD nor it will consider the possibility of not having counter parts in the func distribution.

Parameters:

var (iterable) – Input variables. Their cumulative distribution is computed in nbinV bins between minV and maxV. This cumulative abundance is then matched to the cumulative abundance of func.
func (function) – Differential distribution (i.e. density function) of the properties that have to be associated to variables in the input iterable var. The corresponding cumulative distribution is computed in nbinF bins between minF and maxF using function cumulative_from_differential of module scampy.measure.abundances
minF (scalar, optional) – Minimum value of the observable property grid (default: -20).
maxF (scalar, optional) – Maximum value of the observable property grid (default: -10).
minP (scalar, optional) – Minimum cumulative probability considered valid for inversion; bins with \(P < \text{minP}\) are excluded (default: 1e-20).
maxP (scalar, optional) – Maximum cumulative probability considered valid for inversion (default: 1.0).
nbinF (int, optional) – Number of bins for the observable property grid (default: 200).
factF (scalar, optional) – Normalisation factor applied to func when computing its cumulative distribution (default: 1.0).
minV (scalar or None, optional) – Minimum value of the input variable range. If None (default) it is set to var.min() * 0.99.
maxV (scalar or None, optional) – Maximum value of the input variable range. If None (default) it is set to var.max() * 1.01.
nbinV (int, optional) – Number of bins for the input variable grid (default: 100).
factV (scalar, optional) – Normalisation factor applied to var when computing its cumulative distribution (default: 1.0).
scatter_type (str) –
(Optional, default = None) Type of scatter to add to the output array of properties. Currently implemented:
- 'normal': Gaussian scatter around the predicted value
- None : none of the above, no scatter is added
kw_scatter (dict) –
(Optional, default = {}) Dictionary with the keyword arguments to pass to the function computing the scatter.
- method = 'normal': ‘scale’ = scattering magnitude (in dex). The default value is 1.0 ‘loc’ = if not zero, adds a bias to the relation. The default value is 0.0
rng (random number generator) – (Optional, default = None) if None is passed, it will build one using function numpy.random.default_rng(**kw_rng)
kw_rng (dict) – (Optional, default = { ‘seed’ : 555 }) keyword arguments to pass to the default random number generator. (if rng is not None, this input is ignored)

Returns:

array of properties sampled from func and matching in abundance the values in var

Return type:

ndarray

Note

The output reproducibility (i.e. the possibility to get the same result running the function twice on the same inputs) depends on whether a scattering has been added and, in that case, is guaranteed only if the random number generator passed is in the same state.

scampy.sham.match_cross(X, Y, Xlim=None, Xnbin=100, factX=1.0, logX=False, Ylim=None, Ynbin=100, factY=1.0, logY=False, Xscatter=None, kw_Xscatter={}, Yscatter=None, kw_Yscatter={}, rng=None, kw_rng={'seed': 555})

Find indices of the elements in the X array that match in abundance with elements of the Y array. Returns the two list of indices in bijective relation.

Parameters:

X (array-like) – First input property array (e.g. sub-halo masses).
Y (array-like) – Second input property array (e.g. observed luminosities).
Xlim (tuple or None, optional) – (Xmin, Xmax) range for X. If None (default) the range is set automatically from the data.
Xnbin (int, optional) – Number of bins for the X cumulative distribution (default: 100).
factX (float, optional) – Normalisation factor for the X cumulative distribution (default: 1.0).
logX (bool, optional) – If True bins X on a log scale and scatters in log-space (default: False).
Ylim (tuple or None, optional) – (Ymin, Ymax) range for Y. If None (default) the range is set automatically from the data.
Ynbin (int, optional) – Number of bins for the Y cumulative distribution (default: 100).
factY (float, optional) – Normalisation factor for the Y cumulative distribution (default: 1.0).
logY (bool, optional) – If True bins Y on a log scale and scatters in log-space (default: False).
Xscatter (str) –
(Optional, default = None) Type of scatter to add to the output array of properties. Currently implemented:
- 'normal': Gaussian scatter around the predicted value
- None : none of the above, no scatter is added to the X-array
kw_Xscatter (dict) –
(Optional, default = {}) Dictionary with the keyword arguments to pass to the function computing the scatter.
- method = 'normal': ‘scale’ = scattering magnitude (in dex). The default value is 1.0 ‘loc’ = if not zero, adds a bias to the relation. The default value is 0.0
Yscatter (str) – (Optional, default = None) same as Xscatter for the Y array
kw_Yscatter (dict) – (Optional, default = {}) same as kw_Xscatter for the Y array
rng (random number generator) – (Optional, default = None) if None is passed, it will build one using function numpy.random.default_rng(**kw_rng)
kw_rng (dict) – (Optional, default = { ‘seed’ : 555 }) keyword arguments to pass to the default random number generator. (if rng is not None, this input is ignored)

Returns:

idx_sort (int ndarray)
idy_sort (int ndarray)

Note

The output reproducibility (i.e. the possibility to get the same result running the function twice on the same inputs) depends on whether a scattering has been added and, in that case, is guaranteed only if the random number generator passed is in the same state on all runs.

cosmology extension module

Flat or curved FLRW cosmological model.

Provides the model class, which pre-tabulates cosmographic quantities on a log-spaced redshift grid at construction time for fast evaluation.

Available methods of `model`

Cosmographic distances

`Hz()`	Hubble parameter H(z) [km/s/Mpc].
`dH()`	Hubble distance c/H(z) [Mpc].
`dC()`	Line-of-sight comoving distance [Mpc].
`ddC()`	Derivative of comoving distance dd_C/dz [Mpc].
`dA()`	Angular diameter distance [Mpc].

Volume and time

`comoving_volume_unit()`	Comoving volume element \(dV/(dz\,d\Omega)\) \([(\mathrm{Mpc}/h)^3\,\mathrm{sr}^{-1}]\).
`comoving_volume()`	Total comoving volume V(z) \([(\mathrm{Mpc}/h)^3]\).
`cosmic_time()`	Age of the Universe at redshift z [Gyr].

Density and matter content

`critical_density_comoving()`	Critical density in comoving units \([M_\odot\,\mathrm{Mpc}^{-3}\,h^2]\).
`critical_density()`	Critical density in physical units \([M_\odot\,\mathrm{Mpc}^{-3}]\).
`OmegaM()`	Total matter density parameter \(\Omega_M(z)\).
`Omegab()`	Baryonic matter density parameter \(\Omega_b(z)\).

Structure growth

`deltac()`	Linear critical overdensity \(\delta_c(z)\).
`D()`	Linear growth factor D(z).
`gz()`	Unnormalised growth factor (1+z) D(z).
`Delta_c()`	Virial overdensity \(\Delta_c(z)\) (Nakamura & Suto 1998).

class scampy.cosmology.model(Om_M=0.3, Om_b=0.045, Om_L=0.7, Om_n=0.0, Om_r=0.0, Om_K=0.0, hh=0.7, sigma8=0.8, w_0=-1.0, w_a=0.0, zmin=0.0, zmax=100.0, thin=1000)

Flat or curved FLRW cosmological model.

Pre-tabulates E(z), 1/E(z), the comoving-distance integrand, the linear growth factor, and the cosmic-time integrand on a log-spaced redshift grid at construction time. All public methods are then fast interpolations or closed-form expressions.

Parameters:

Om_M (float, optional) – Total matter density parameter \(\Omega_M\) (default: 0.3).
Om_b (float, optional) – Baryonic matter density parameter \(\Omega_b\) (default: 0.045).
Om_L (float, optional) – Dark-energy density parameter \(\Omega_\Lambda\) (default: 0.7).
Om_n (float, optional) – Massive-neutrino density parameter (default: 0.0).
Om_r (float, optional) – Radiation density parameter (default: 0.0).
Om_K (float, optional) – Curvature density parameter (default: 0.0).
hh (float, optional) – Dimensionless Hubble constant \(h = H_0/100\) (default: 0.7).
sigma8 (float, optional) – RMS matter-fluctuation amplitude at \(8\,h^{-1}\,\mathrm{Mpc}\) (default: 0.8).
w_0 (float, optional) – Dark-energy equation-of-state constant term (default: -1.0).
w_a (float, optional) – Dark-energy equation-of-state linear term \(w(a) = w_0 + w_a(1-a)\) (default: 0.0).
zmin (float, optional) – Minimum redshift of the internal tabulation grid (default: 0.0).
zmax (float, optional) – Maximum redshift of the internal tabulation grid (default: 100.0).
thin (int, optional) – Number of points in the internal tabulation grid (default: 1000).

Hz(zz)

Hubble parameter H(z) in km/s/Mpc.

Parameters:: zz (float or array-like) – Redshift.
Returns:: H(z) [km/s/Mpc].
Return type:: float or ndarray

dH(zz)

Hubble distance c/H(z) in Mpc.

Parameters:: zz (float or array-like) – Redshift.
Returns:: c/H(z) [Mpc].
Return type:: float or ndarray

dC(zz)

Line-of-sight comoving distance in Mpc.

Parameters:: zz (float or array-like) – Redshift.
Returns:: d_C(z) [Mpc].
Return type:: float or ndarray

ddC(zz)

Derivative of the comoving distance with respect to redshift, dd_C/dz.

Parameters:: zz (float or array-like) – Redshift.
Returns:: dd_C/dz [Mpc].
Return type:: float or ndarray

dA(zz)

Angular diameter distance in Mpc.

Parameters:: zz (float or array-like) – Redshift.
Returns:: d_A(z) [Mpc].
Return type:: float or ndarray

comoving_volume_unit(zz)

Comoving volume element per unit redshift per unit steradian, dV/(dz dOmega).

Parameters:: zz (float or array-like) – Redshift.
Returns:: dV/(dz dOmega) [(Mpc/h)^3 / sr].
Return type:: float or ndarray

comoving_volume(zz)

Total comoving volume enclosed within redshift z in (Mpc/h)^3.

Parameters:: zz (float or array-like) – Redshift.
Returns:: V_C(z) [(Mpc/h)^3].
Return type:: float or ndarray

cosmic_time(zz)

Age of the Universe at redshift z (lookback time from z to infinity) in Gyr.

Parameters:: zz (float or array-like) – Redshift.
Returns:: t(z) [Gyr].
Return type:: float or ndarray

critical_density_comoving(zz)

Critical density in comoving units [Msol Mpc^-3 h^2].

Parameters:: zz (float or array-like) – Redshift.
Returns:: rho_crit(z) [Msol Mpc^-3 h^2].
Return type:: float or ndarray

critical_density(zz)

Critical density in physical units [Msol Mpc^-3].

Parameters:: zz (float or array-like) – Redshift.
Returns:: rho_crit(z) [Msol Mpc^-3].
Return type:: float or ndarray

OmegaM(zz)

Total matter density parameter Omega_M(z) = rho_M(z)/rho_crit(z).

Parameters:: zz (float or array-like) – Redshift.
Returns:: Omega_M(z).
Return type:: float or ndarray

Omegab(zz)

Baryonic matter density parameter Omega_b(z).

Parameters:: zz (float or array-like) – Redshift.
Returns:: Omega_b(z).
Return type:: float or ndarray

deltac(zz)

Linear critical overdensity for spherical collapse delta_c(z).

Parameters:: zz (float or array-like) – Redshift.
Returns:: delta_c(z).
Return type:: float or ndarray

D(zz)

Linear growth factor D(z), unnormalised (divide by D(0) for D(0)=1 convention).

Parameters:: zz (float or array-like) – Redshift.
Returns:: D(z).
Return type:: float or ndarray

gz(zz)

Unnormalised growth factor g(z) = (1+z)*D(z).

Parameters:: zz (float or array-like) – Redshift.
Returns:: (1+z)*D(z).
Return type:: float or ndarray

Delta_c(zz)

Virial overdensity Delta_c(z) from Nakamura & Suto (1998).

Parameters:: zz (float or array-like) – Redshift.
Returns:: Delta_c(z).
Return type:: float or ndarray

halo subpackage

halo.mass_function module

Halo mass functions.

Provides fitting functions for the halo mass function \(n(M_h, z)\), i.e. the comoving number density of dark matter haloes per unit mass interval. All functions share the same call signature and return \(\mathrm{d}n/\mathrm{d}M_h\) in units of \([M_\odot^{-1}\,h^3\,\mathrm{Mpc}^{-3}]\) (comoving) or \([M_\odot^{-1}\,\mathrm{Mpc}^{-3}]\) (physical).

scampy.halo.mass_function.ShethTormen01(mm, zz, pk, comoving=True)

Sheth & Tormen (2001) halo mass function.

Implements the fitting formula

\[\frac{\mathrm{d}n}{\mathrm{d}M} = \rho_0 \left|\frac{\mathrm{d}\ln\sigma}{\mathrm{d}M}\right| f(\nu)\,/\,M\]

where

\[f(\nu) = A\sqrt{\frac{2}{\pi}}\,0.3222 \left(1 + (a\nu^2)^{-0.3}\right)\nu\, e^{-a\nu^2/2}, \quad a = 0.707,\quad A \approx 0.8409,\]

and \(\nu = \delta_c(z)\,/\,[D(z)\,\sigma(M)]\) is the peak-height parameter.

Parameters:

mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s) at which to evaluate the mass function. The output is broadcast over (mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object carrying the cosmological model and the normalised growth factor.
comoving (bool, optional) – If True (default) masses and densities are in comoving units.

Returns:

Array of shape (mm.size,) or (mm.size, zz.size) containing \(\mathrm{d}n/\mathrm{d}M_h\).

Return type:

ndarray

scampy.halo.mass_function.Tinker08(mm, zz, pk, comoving=True)

Tinker et al. (2008) halo mass function.

Implements the fitting formula

\[\frac{\mathrm{d}n}{\mathrm{d}M} = \rho_0 \left|\frac{\mathrm{d}\ln\sigma}{\mathrm{d}M}\right| f(\sigma)\,/\,M,\]

where

\[f(\sigma) = A_z \left[\left(\frac{\sigma}{b_z}\right)^{-a_z} + 1\right] e^{-c_0/\sigma^2}\]

with redshift-dependent coefficients

\[A_z = A_0\,(1+z)^{-0.14},\quad a_z = a_0\,(1+z)^{-0.06},\quad b_z = b_0\,(1+z)^{-\alpha},\]

calibrated for a \(\Delta = 200\) overdensity threshold (\(A_0=0.186,\,a_0=1.47,\,b_0=2.57,\,c_0=1.19\)).

Parameters:

mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s) at which to evaluate the mass function. The output is broadcast over (mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object carrying the cosmological model and growth factor.
comoving (bool, optional) – If True (default) masses and densities are in comoving units.

Returns:

Array of shape (mm.size,) or (mm.size, zz.size) containing \(\mathrm{d}n/\mathrm{d}M_h\).

Return type:

ndarray

scampy.halo.mass_function.Behroozi13(mm, zz, pk, comoving=True)

Behroozi, Wechsler & Conroy (2013) halo mass function.

Applies a high-redshift correction on top of Tinker08():

\[n_\text{B13}(M,z) = 10^{N(a)}\cdot\left(3.16\times10^{-12}\,M\right)^{p(a)} \cdot n_\text{T08}(M,z),\]

where \(a = 1/(1+z)\) is the scale factor,

\[N(a) = \frac{0.144}{1 + e^{14.79\,(a - 0.213)}},\quad p(a) = \frac{0.5}{1 + e^{6.5\,a}}.\]

The correction is negligible at low redshift and boosts the abundance of low-mass haloes at high redshift.

Parameters:

mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s) at which to evaluate the mass function. The output is broadcast over (mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object carrying the cosmological model and growth factor.
comoving (bool, optional) – If True (default) masses and densities are in comoving units.

Returns:

Array of shape (mm.size,) or (mm.size, zz.size) containing \(\mathrm{d}n/\mathrm{d}M_h\).

Return type:

ndarray

scampy.halo.mass_function.Bhattacharya11(mm, zz, pk, comoving=True)

Generalised halo mass function of Bhattacharya et al. (2011) with scale- and redshift-dependent coefficients calibrated on the Euclid SUBFIND catalogues by Castro et al. (2023, Tab. 4).

The multiplicity function is (Eq. 3 of Castro et al. 2023):

\[\nu f(\nu) = A(p,q) \sqrt{\frac{2a\nu^2}{\pi}} e^{-a\nu^2/2} \left[1 + \frac{1}{(a\nu^2)^p}\right] (\nu\sqrt{a})^{q-1},\]

where the normalisation \(A(p,q)\) (Eq. 4) ensures \(\int f(\nu)\,\mathrm{d}\nu = 1\). The parameters \(a\), \(p\), \(q\) depend on both the local logarithmic slope of \(\sigma(R)\) and \(\Omega_M(z)\) via Eqs. 12–16.

Parameters:

mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s) at which to evaluate the mass function. The output is broadcast over (mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object carrying the cosmological model and growth factor.
comoving (bool, optional) – If True (default) masses and densities are in comoving units.

Returns:

Array of shape (mm.size,) or (mm.size, zz.size) containing \(\mathrm{d}n/\mathrm{d}M_h\).

Return type:

ndarray

halo.bias module

Halo bias functions.

Provides fitting functions for the large-scale linear halo bias \(b(M_h, z)\), defined such that the halo–matter power spectrum equals \(b^2(M_h)\,P_\mathrm{lin}(k,z)\) on large scales.

scampy.halo.bias.ShethMoTormen01(mm, zz, pk, comoving=True)

Sheth, Mo & Tormen (2001) halo bias.

Implements the fitting formula

\[b(M_h, z) = 1 + \frac{1}{\delta_c} \left[ a\nu^2 + \frac{a\nu^2}{(a\nu^2)^{0.6}} - \frac{(a\nu^2)^{0.6}}{(a\nu^2)^{0.6} + 0.14} \right],\quad a = 0.707,\]

where \(\nu = \delta_c(z)\,/\,[D(z)\,\sigma(M_h)]\) is the peak-height parameter and \(\delta_c(z)\) is the linear collapse threshold.

Parameters:

mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s) at which to evaluate the bias. The output is broadcast over (mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object carrying the cosmological model and growth factor.
comoving (bool, optional) – If True (default) masses are in comoving units (currently unused, reserved for future physical-unit support).

Returns:

Array of shape (mm.size,) or (mm.size, zz.size) containing \(b(M_h, z)\).

Return type:

ndarray

scampy.halo.bias.Tinker10(mm, zz, pk, Delta=200, comoving=True)

Tinker et al. (2010) halo bias.

Implements the fitting formula

\[b(M_h, z) = 1 - A\,\frac{\nu^a}{\nu^a + \delta_c^a} + B\,\nu^b + C\,\nu^c,\]

with coefficients depending on the overdensity threshold \(\Delta\) via \(y = \log_{10}\Delta\):

\[A = 1 + 0.24\,y\,e^{-(4/y)^4},\quad a = 0.44\,y - 0.88,\quad B = 0.183,\quad b = 1.5,\quad C = 0.019 + 0.107\,y + 0.19\,e^{-(4/y)^4},\quad c = 2.4.\]

Parameters:

mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s) at which to evaluate the bias. The output is broadcast over (mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object carrying the cosmological model and growth factor.
Delta (float, optional) – Overdensity threshold with respect to the mean matter density (default: 200).
comoving (bool, optional) – If True (default) masses are in comoving units (currently unused, reserved for future physical-unit support).

Returns:

Array of shape (mm.size,) or (mm.size, zz.size) containing \(b(M_h, z)\).

Return type:

ndarray

halo.density_profile module

NFW halo density profile and halo concentration models.

Provides the Fourier transform of the NFW density profile and several empirical fitting functions for the halo concentration–mass relation \(c(M_h, z)\).

scampy.halo.density_profile.concentration_giocoli12(mm, zz, pk, comoving=True): Giocoli et al. (2012) halo concentration–mass relation.

Note

Not yet implemented; returns None.

scampy.halo.density_profile.concentration_zhao09(mm, zz, pk, comoving=True): Zhao et al. (2009) halo concentration–mass relation.

Note

Not yet implemented; returns None.

scampy.halo.density_profile.concentration_shimizu03(mm, zz, pk, comoving=True)

Shimizu et al. (2003) halo concentration–mass relation.

Implements Eq. 4 of Shimizu et al. (2003):

\[c(M_h, z) = \frac{8}{1+z}\, \left(1.0204\times10^{-14}\,M_h\,h\right)^{-0.13}.\]

Parameters:

mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s). The output is broadcast over (mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object; used to retrieve \(h\) when comoving=False.
comoving (bool, optional) – If True (default) masses are already in \(M_\odot\,h^{-1}\) and no \(h\)-conversion is applied.

Returns:

Array of shape (mm.size,) or (mm.size, zz.size) containing the dimensionless concentration \(c = r_\mathrm{vir}/r_s\).

Return type:

ndarray

scampy.halo.density_profile.density_profile_FT(kk, mm, zz, pk, comoving=True)

Fourier transform of the NFW density profile.

Returns the normalised Fourier transform \(\tilde{u}(k,z|M_h)\) of the NFW profile for a halo of mass \(M_h\), as given by Eqs. 8–9 of Ronconi et al. (2020):

\[\tilde{u}(k, z|M_h) = \frac{ \cos(\mu)\,[\mathrm{Ci}(\mu(1+c)) - \mathrm{Ci}(\mu)] + \sin(\mu)\,[\mathrm{Si}(\mu(1+c)) - \mathrm{Si}(\mu)] - \dfrac{\sin(c\mu)}{\mu(1+c)} }{\ln(1+c) - c/(1+c)},\]

where \(\mu = k\,r_s\), \(c = c(M_h, z)\) is the halo concentration, \(r_s = r_\mathrm{vir}/c\) is the NFW scale radius, and \(\mathrm{Si}\), \(\mathrm{Ci}\) are the sine and cosine integrals. The concentration is computed via concentration_shimizu03().

Parameters:

kk (scalar or array-like) – Wavenumbers in \([h\,\mathrm{Mpc}^{-1}]\).
mm (scalar or array-like) – Halo masses \(M_h\) in \([M_\odot\,h^{-1}]\).
zz (scalar or array-like) – Redshift(s). The output is broadcast over (kk, mm, zz).
pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object; used to retrieve the cosmological model (critical density, \(\Delta_c\)).
comoving (bool, optional) – If True (default) the virial volume uses the comoving critical density; otherwise the physical critical density is used.

Returns:

Normalised profile transform \(\tilde{u}(k,z|M_h)\), dimensionless and in the range \([0, 1]\).

Return type:

ndarray

halo.model module

Halo model for the two-point statistics of a galaxy population.

Implements the halo-model framework (Cooray & Sheth 2002) as described in Sec. 2.1 of Ronconi et al. (2020). Given a power spectrum and an HOD, the halo_model class provides the galaxy number density, effective bias, and the 3D, projected, and angular two-point correlation functions decomposed into their 1-halo and 2-halo contributions.

class scampy.halo.model.halo_model(pk, comoving=True, Mlim=(5, 17), klim=(-3, 3), thin=256, hmf='Tinker08', hbias='Tinker10', density='NFW')

Halo model for galaxy two-point statistics.

Encapsulates the halo-model integrals over the halo mass function, bias, and NFW density profile (see Ronconi et al. 2020, Sec. 2.1). All power-spectrum and correlation-function methods accept an HOD instance that defines the galaxy occupation probabilities \(\langle N_\mathrm{cen}\rangle(M_h)\) and \(\langle N_\mathrm{sat}\rangle(M_h)\).

Parameters:

pk (scampy.power_spectrum.power_spectrum) – Power-spectrum object carrying the cosmological model and growth factor.
comoving (bool, optional) – Whether to work in comoving coordinates (default: True).
Mlim (tuple of float, optional) – Log10 of the minimum and maximum halo mass for the integration grid, in \(M_\odot\,h^{-1}\) (default: (5, 17)).
klim (tuple of float, optional) – Log10 of the minimum and maximum wavenumber for the internal \(k\)-grid in \(h\,\mathrm{Mpc}^{-1}\) (default: (-3, +3)).
thin (int, optional) – Number of points in both the mass and wavenumber grids (default: 256).
hmf (str, optional) – Halo mass function label (placeholder, currently always Tinker08()).
hbias (str, optional) – Halo bias label (placeholder, currently always Tinker10()).
density (str, optional) – Density profile label (placeholder, currently always NFW via density_profile_FT()).

ng(hod, zz=0.0)

Mean galaxy number density.

Integrates the HOD-weighted halo mass function (Eq. 5 of Ronconi et al. 2020):

\[n_g(z) = \int_{M_\mathrm{min}}^{M_\mathrm{max}} \langle N_g \rangle(M_h)\,n(M_h, z)\,\mathrm{d}M_h,\]

where \(\langle N_g\rangle = \langle N_\mathrm{cen}\rangle + \langle N_\mathrm{sat}\rangle\).

Parameters:

hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float or array-like, optional) – Redshift(s) (default: 0.0).

Returns:

Galaxy number density in \([h^3\,\mathrm{Mpc}^{-3}]\).

Return type:

float or ndarray

bias(hod, zz=0.0)

Effective large-scale halo bias of the galaxy population.

\[b_g(z) = \frac{1}{n_g(z)} \int \langle N_g\rangle(M_h)\, n(M_h, z)\,b(M_h, z)\,\mathrm{d}M_h.\]

Parameters:

hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float or array-like, optional) – Redshift(s) (default: 0.0).

Returns:

Dimensionless effective bias \(b_g(z)\).

Return type:

float or ndarray

Mhalo(hod, zz=0.0)

Mean host halo mass of the galaxy population.

\[\langle M_h \rangle(z) = \frac{1}{n_g(z)} \int \langle N_g\rangle(M_h)\, n(M_h, z)\,M_h\,\mathrm{d}M_h.\]

Parameters:

hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float or array-like, optional) – Redshift(s) (default: 0.0).

Returns:

Mean host halo mass in \([M_\odot\,h^{-1}]\).

Return type:

float or ndarray

dngdM(mm, hod, zz=0.0)

Differential galaxy number density per unit halo mass.

Returns the integrand of ng() evaluated at specific mass values:

\[\frac{\mathrm{d}n_g}{\mathrm{d}M}(M_h, z) = \langle N_g\rangle(M_h)\,n(M_h, z).\]

Parameters:

mm (scalar or array-like) – Halo masses in \([M_\odot\,h^{-1}]\).
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float or array-like, optional) – Redshift(s) (default: 0.0).

Returns:

Differential density with shape (mm.size,) or (mm.size, zz.size) in \([M_\odot^{-1}\,h^4\,\mathrm{Mpc}^{-3}]\).

Return type:

ndarray

Pk_1halo(kk, hod, zz=0.0, fact=-1.0)

1-halo term of the galaxy power spectrum.

Pairs within the same halo (central–satellite and satellite–satellite), Eq. 10 of Ronconi et al. (2020):

\[P_\mathrm{1h}(k,z) = \frac{1}{n_g^2(z)} \int \left(2\langle N_\mathrm{cen}\rangle + \langle N_\mathrm{sat}\rangle\right) \langle N_\mathrm{sat}\rangle\, n(M_h,z)\,|\tilde{u}(k,z|M_h)|^2\,\mathrm{d}M_h.\]

Parameters:

kk (scalar or array-like) – Wavenumbers in \([h\,\mathrm{Mpc}^{-1}]\).
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float or array-like, optional) – Redshift(s) (default: 0.0).
fact (float or array-like, optional) – Pre-computed normalisation factor \(1/n_g^2\). If negative (default) it is computed internally.

Returns:

1-halo power spectrum with shape (kk.size,) or (kk.size, zz.size) in \([h^{-3}\,\mathrm{Mpc}^3]\).

Return type:

ndarray

Pk_2halo(kk, hod, zz=0.0, fact=-1.0)

2-halo term of the galaxy power spectrum.

Pairs from different haloes, Eq. 11 of Ronconi et al. (2020):

\[P_\mathrm{2h}(k,z) = \frac{P_\mathrm{lin}(k,z)}{n_g^2(z)} \left[\int \langle N_g\rangle(M_h)\, n(M_h,z)\,b(M_h,z)\, \tilde{u}(k,z|M_h)\,\mathrm{d}M_h\right]^2.\]

Parameters:

kk (scalar or array-like) – Wavenumbers in \([h\,\mathrm{Mpc}^{-1}]\).
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float or array-like, optional) – Redshift(s) (default: 0.0).
fact (float or array-like, optional) – Pre-computed normalisation factor \(1/n_g^2\). If negative (default) it is computed internally.

Returns:

2-halo power spectrum with shape (kk.size,) or (kk.size, zz.size) in \([h^{-3}\,\mathrm{Mpc}^3]\).

Return type:

ndarray

Pk(kk, hod, zz=0.0, fact=-1.0)

Total halo-model galaxy power spectrum.

Sum of the 1-halo and 2-halo terms (Eq. 7 of Ronconi et al. 2020):

\[P(k,z) = P_\mathrm{1h}(k,z) + P_\mathrm{2h}(k,z).\]

Parameters:

kk (scalar or array-like) – Wavenumbers in \([h\,\mathrm{Mpc}^{-1}]\).
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float or array-like, optional) – Redshift(s) (default: 0.0).
fact (float or array-like, optional) – Pre-computed normalisation factor \(1/n_g^2\). If negative (default) it is computed internally.

Returns:

Total power spectrum with shape (kk.size,) or (kk.size, zz.size) in \([h^{-3}\,\mathrm{Mpc}^3]\).

Return type:

ndarray

Xi_1halo(rr, hod, zz=0.0, fact=-1.0)

1-halo term of the 3D two-point correlation function.

Obtained by Fourier-transforming Pk_1halo() via a Fast Spherical Transform (Eq. 6 of Ronconi et al. 2020).

Parameters:

rr (scalar or array-like) – Comoving separations in \([h^{-1}\,\mathrm{Mpc}]\).
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float, optional) – Redshift (default: 0.0). Broadcasting over zz is not yet supported.
fact (float, optional) – Pre-computed normalisation factor. If negative (default) it is computed internally.

Returns:

1-halo correlation function \(\xi_\mathrm{1h}(r)\) evaluated at rr.

Return type:

ndarray

Xi_2halo(rr, hod, zz=0.0, fact=-1.0)

2-halo term of the 3D two-point correlation function.

Fourier transform of Pk_2halo(). Same signature as Xi_1halo().

Xi(rr, hod, zz=0.0, fact=-1.0)

Total 3D two-point correlation function.

Sum of Xi_1halo() and Xi_2halo(). Same signature as Xi_1halo().

Wr_1halo(rp, hod, zz=0.0, fact=-1.0)

1-halo term of the projected two-point correlation function.

Obtained by Abel-transforming Pk_1halo() via a Fast Projected Spherical Transform (Eq. 15 of Ronconi et al. 2020, zeroth-order Hankel transform).

Parameters:

rp (scalar or array-like) – Projected comoving separations in \([h^{-1}\,\mathrm{Mpc}]\).
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float, optional) – Redshift (default: 0.0). Broadcasting over zz is not yet supported.
fact (float, optional) – Pre-computed normalisation factor. If negative (default) it is computed internally.

Returns:

1-halo projected correlation function \(w_\mathrm{1h}(r_p)\) evaluated at rp.

Return type:

ndarray

Wr_2halo(rp, hod, zz=0.0, fact=-1.0)

2-halo term of the projected two-point correlation function.

Abel transform of Pk_2halo(). Same signature as Wr_1halo().

Wr(rp, hod, zz=0.0, fact=-1.0)

Total projected two-point correlation function.

Sum of Wr_1halo() and Wr_2halo(). Same signature as Wr_1halo().

Wt_1halo(th, hod, zz=0.0, fact=-1.0)

1-halo term of the angular two-point correlation function.

Converts the angular separation \(\theta\) to a projected comoving separation \(r_p = \theta\,d_C(z)\) and returns Wr_1halo() (flat-sky / Limber approximation, Eqs. 15–16 of Ronconi et al. 2020).

Parameters:

th (scalar or array-like) – Angular separations in radians.
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zz (float, optional) – Redshift (default: 0.0). Broadcasting over zz is not yet supported.
fact (float, optional) – Pre-computed normalisation factor. If negative (default) it is computed internally.

Returns:

1-halo angular correlation function \(\omega_\mathrm{1h}(\theta)\) evaluated at th.

Return type:

ndarray

Wt_2halo(th, hod, zz=0.0, fact=-1.0)

2-halo term of the angular two-point correlation function.

Converts \(\theta\) to \(r_p\) and returns Wr_2halo(). Same signature as Wt_1halo().

Wt(th, hod, zz=0.0, fact=-1.0)

Total angular two-point correlation function.

Sum of Wt_1halo() and Wt_2halo(). Same signature as Wt_1halo().

Nz(hod, zbins, solid_angle=12.566370614359172)

Normalised redshift distribution of the galaxy population.

Computes the fraction of sources per unit redshift interval:

\[N(z) = \frac{n_g(z)\,\frac{\mathrm{d}V_C}{\mathrm{d}z}\, \Omega\,\Delta z}{\Delta z\,\sum_i n_g(z_i)\, \frac{\mathrm{d}V_C}{\mathrm{d}z_i}\,\Omega\,\Delta z_i},\]

where \(\mathrm{d}V_C/\mathrm{d}z\) is the comoving volume element and \(\Omega\) is the solid angle.

Parameters:

hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zbins (array-like) – Bin edges in redshift. The bin centres are taken as the midpoints.
solid_angle (float, optional) – Solid angle of the survey in steradians (default: \(4\pi\), full sky).

Returns:

Normalised \(N(z)\) with zbins.size - 1 elements, such that \(\sum_i N(z_i)\,\Delta z_i = 1\).

Return type:

ndarray

Wt_zdist(th, hod, zbins, Nzdist=None, fact=-1, valid_rp=(0.01, 80.0), solid_angle=12.566370614359172, component='total')

Angular two-point correlation function averaged over a redshift distribution.

Implements Eq. 17 of Ronconi et al. (2020) in the Limber approximation:

\[\omega(\theta) = \int \frac{\mathrm{d}V}{\mathrm{d}z}\, N^2(z)\, \omega\!\left[r_p(\theta, z),\, z\right]\mathrm{d}z,\]

where \(N(z)\) is the normalised source redshift distribution, \(r_p(\theta, z) = \theta\,d_C(z)\) is the projected comoving separation at the comoving distance \(d_C(z)\), and \(\mathrm{d}V/\mathrm{d}z\) is the comoving volume element. The integral is evaluated by computing \(\omega(r_p, z)\) via a Fast Projected Spherical Transform at each redshift bin centre and summing with trapezoidal weights.

See also

Nz: Computes \(N(z)\) from the HOD and cosmology; called internally when Nzdist is not provided.

Parameters:

th (array-like) – Angular separations in radians.
hod (scampy.hod.HOD or compatible) – Object providing Pcen(M) and Psat(M) methods.
zbins (array-like) – Redshift bin edges. Bin centres are taken as midpoints.
Nzdist (array-like or None, optional) – Normalised redshift distribution \(N(z_i)\) as returned by \(Nz\), with Nzdist.size == zbins.size - 1. If None (default) it is computed internally via \(Nz\).
fact (float or array-like, optional) – Pre-computed normalisation factor \(1/n_g^2\). If negative (default) it is computed internally.
valid_rp (tuple of float, optional) – Range of projected separations \((r_{p,\min}, r_{p,\max})\) in \([h^{-1}\,\mathrm{Mpc}]\) over which the interpolation of the fixed-\(z\) angular correlation is trusted (default: (1e-2, 8e1)).
solid_angle (float, optional) – Survey solid angle in steradians, used when Nzdist is computed internally (default: \(4\pi\), full sky).
component (str, optional) – Which halo-model term to include: 'total' (default), '1halo', or '2halo'.

Returns:

Angular correlation function \(\omega(\theta)\) with the same shape as th.

Return type:

ndarray

io subpackage

io.hdf5 module

scampy.io.hdf5.write_to_hdf5(outfile, metadata, hard=False, **groups)

Automatic dump of several dictionaries into hdf5 file as separate groups in file. Also support generic metadata

Parameters:

outfile (string) – /path/to/outputfile.hdf5 destination file. The path should exist on the filesystem
metadata (dict) – dictionary with global metadata
hard (bool) – (Optional, default=False) if True eventually overwrites an already existing file with the same name
**groups (Keyword arguments) – (Optional) Dictionaries to save to the file. Each dictionary will be saved to the output file as a group with its corresponding keyword name as keyword.

Return type:

None

See also

load_from_hdf5: loads all the groups in a .hdf5 file

scampy.io.hdf5.load_from_hdf5(infile)

Loads the groups in a hdf5 file into a dictionary.

Parameters:: infile (string) – /path/to/inputfile.hdf5 input file
Returns:: nested dictionary
Return type:: dict

See also

write_to_hdf5: dump of several dictionaries into hdf5 file as separate groups in file

Warning

This is a simplified reader intended to work with hdf5 files generate using the scampy.io.hdf5.write_to_hdf5 method. It is not guaranteed to work with any .hdf5. It assumes that at the top level all the fields can be cast to dictionaries. If this is not possible it raises a warning.

io.subfind_table module

Reader for GADGET SUBFIND binary subgroup table files.

Provides subfind_table, which parses the binary format written by the GADGET SUBFIND halo-finder and exposes halo and sub-halo properties as plain NumPy arrays.

class scampy.io.subfind_table.subfind_table(filebase, byteorder='little', ids_lenght=8, masstab=True)

A class for reading GaDGET SUBFIND subgroup tables

Parameters:

filebase (string) – common name of the files to be read
byteorder (string) – ‘little’ or ‘big’ for little or big endian, respectively
ids_lenght (int) – number of bytes of the integer storing the IDs (8 or 16)
masstab (bool) – whether the file contains mass tables

generate_content_dict()

Allocate an empty content dictionary for all catalogue fields.

Returns:: A dictionary mapping each field name in fields to an empty numpy.ndarray with the shape and dtype specified in values.
Return type:: dict

read_header(num=0)

Reads only the Header of the num`^th file in `base

Parameters:

num (the file number to read)

Returns:

loc (dictionary with meta-data local to file num)
glob (dictionary with meta-data global to all files in base)

read_file(num, scale_mass=10000000000.0, scale_lenght=0.001, add_to_content=False, verbose=False)

Reads the num`^th file in `base

Parameters:

num (int) – file to read
scale_mass (float) – mass unit (in terms of solar masses, default = 1.e+10)
scale_lenght (float) – lenght unit (in terms of Mpc/h, default = 1.e-3)
add_to_content (bool) – whether to add the data to the internal storage array of the class (default = False)
verbose (bool) – whether to print on screen additional info (default = False)

Return type:

None

read_all_files(reset=True, **kwargs)

Read all files in the set and accumulate their data.

Resets content to an empty dictionary, then iterates over all Nfiles files calling read_file() with add_to_content=True so that each file’s data is appended to the shared arrays.

Parameters:

reset (bool, optional) – Reserved for future use; currently ignored (default: True).
**kwargs – Additional keyword arguments forwarded to read_file() (e.g. scale_mass, scale_lenght, verbose).

Return type:

None

measure subpackage

measure.abundances module

Functions computing abundance statistics of different quantities that can be extracted from catalogues.

scampy.measure.abundances.get_abundances(Ax, Nc, Ns, bins, statistic='mean')

Compute the mean number of central and satellite objects within bins along a given axis.

Parameters:

Axarray-like: the axis along which we want to compute the average values of Ncen and Nsat
Nc, Nsarray-like: the value of Ncen and Nsat for all the values in Ax
binsarray-like: edges of the bins
statisticstring or callable: (Optional, default = 'mean') argument statistic to be passed to function scipy.stats.binned_statistic used to compute the abundances

Returns:

: array-like: bin-centre (computed linearly)
<Nc>, <Ns>array-like: tuple with the chosen statistic of Ncen and Nsat (respectively, default is mean), within the provided bins. It has the same shape of bins

scampy.measure.abundances.differential_counts(var, bins, fact=1.0)

Obtain the differential distribution of a variable \(v\). Computes \(\frac{dN}{dx} \cdot const\), where \(dN\) is the number of \(v\) values in the bin, \(dv\) it the size of the bin and \(const\) is a constant value.

Parameters:

var (array-like) – the list of \(X\) values
bins (array-like) – it defines the bin edges, including the rightmost edge, allowing for non-uniform bin widths.
fact (scalar) – the value of \(const\), constant factor by which to multiply the result

Returns:

array-like – bin center (shape = len(bins) - 1)
array-like – sequence of values with the values of \(\frac{dN}{dv} \cdot const\) (shape = len(bins) - 1)
array-like – Poisson errors on the differential counts, defined as \(\frac{\sqrt{dN}}{dv} \cdot const\) (shape = len(bins) - 1)

scampy.measure.abundances.cumulative_counts(X, bins, fact=1.0)

Obtain the cumulative distribution of a variable \(X\). Computes \(N(X>x_i) \cdot const\), where \(N(X>x_i)\) is the number of \(X\) values greater than \(x_i\) and \(const\) is a constant value.

Parameters:

var (array-like) – the list of \(X\) values
binning (array-like) – defines the binned values \(x_i\), allowing for non-uniform bin widths.
fact (scalar) – the value of \(const\), constant factor by which to multiply the result

Returns:

array-like – sequence of values with the result of \(N(X>x_i) \cdot const\) (shape = binning.shape)
array-like – Poisson errors on the cumulative counts, defined as \(\sqrt{N(X>x_i)} \cdot const\) (shape = binning.shape)

scampy.measure.abundances.cumulative_from_differential(func, bins, fact=1.0)

Obtain the cumulative distribution from a differential distribution \(f(X)\). Computes

\[F(X > x) \equiv const \cdot \int_{-\infty}^x f(X) dX\]

where \(F(X>x)\) is the cumulative distribution and \(const\) is a constant value.

Parameters:

func (function) – the differential function \(f(X)\), it should take only one argument (a.k.a. the value of \(X\)). It allows lambda expressions.
bins (array-like) – defines the list \(x_i\) values, allowing for non-uniform bin widths.
fact (scalar) – the value of \(const\), constant factor by which to multiply the result

Returns:

sequence of values with the result of \(F(X>x_i)\) (shape = bins.shape)

Return type:

array-like

measure.clustering module

Two-point correlation function estimators.

Provides standard and Landy–Szalay estimators for the two-point correlation function \(\xi(r)\) (Eq. 23 of Ronconi et al. 2020), together with a bootstrap error estimator. Pair counting is delegated to the compiled C++ extension scampy.measure.clustering_core.

scampy.measure.clustering.two_point_standard(data, rand, rbins, omp=True, angular=False)

Two-point correlation function with the standard estimator.

Computes \(\xi(r) = DD/RR - 1\), where \(DD\) and \(RR\) are the normalised data–data and random–random pair counts in each separation bin.

Parameters:

data (ndarray, shape (Ndim, Nobj)) – Coordinates of the data catalogue. Ndim must be 2 or 3.
rand (ndarray, shape (Ndim, Nrand)) – Coordinates of the random catalogue, same dimensionality as data.
rbins (array-like) – Bin edges for the separation \(r\).
omp (bool, optional) – Use the OpenMP-parallel pair counter (default: True).
angular (bool, optional) – Reserved for future angular-space counting (currently unused).

Returns:

Correlation function \(\xi(r)\) in each bin, shape (rbins.size - 1,).

Return type:

ndarray

scampy.measure.clustering.two_point_landyszalay(data, rand, rbins, omp=True, return_error=False, angular=False)

Two-point correlation function with the Landy–Szalay estimator.

Implements Eq. 23 of Ronconi et al. (2020):

\[\xi(r) = \frac{DD(r) - 2\,DR(r) + RR(r)}{RR(r)},\]

where \(DD\), \(DR\), and \(RR\) are the normalised data–data, data–random, and random–random pair counts. This is the preferred estimator used in the paper’s validation tests.

Parameters:

data (ndarray, shape (Ndim, Nobj)) – Coordinates of the data catalogue. Ndim must be 2 or 3.
rand (ndarray, shape (Ndim, Nrand)) – Coordinates of the random catalogue.
rbins (array-like) – Bin edges for the separation \(r\).
omp (bool, optional) – Use the OpenMP-parallel pair counter (default: True).
return_error (bool, optional) – If True, also return a per-bin error estimate derived from the Poisson fluctuations of the pair counts (default: False).
angular (bool, optional) – Reserved for future angular-space counting (currently unused).

Returns:

xi (ndarray) – Correlation function \(\xi(r)\), shape (rbins.size - 1,).
exi (ndarray) – Per-bin error estimate, only returned when return_error=True.

scampy.measure.clustering.bootstrap_two_point(data, rand, rbins, standard=True, Nboots=10, return_boots=False, omp=True, verbose=True, angular=False, rng=None, kw_rng={'seed': 555})

Two-point correlation function with bootstrap error estimate.

Computes a baseline \(\xi(r)\) and estimates its uncertainty by resampling the data catalogue with replacement Nboots times. The random catalogue is fixed across all resamples (only \(RR\) is computed once). Either the standard or the Landy–Szalay estimator can be used for each resample.

Parameters:

data (ndarray, shape (Ndim, Nobj)) – Coordinates of the data catalogue. Ndim must be 2 or 3, or 2 when angular=True.
rand (ndarray, shape (Ndim, Nrand)) – Coordinates of the random catalogue.
rbins (array-like) – Bin edges for the separation \(r\) (or angular separation in radians when angular=True).
standard (bool, optional) – If True (default) use the standard estimator; otherwise use the Landy–Szalay estimator.
Nboots (int, optional) – Number of bootstrap resamples (default: 10).
return_boots (bool, optional) – If True, also return the full array of bootstrap realisations (default: False).
omp (bool, optional) – Use the OpenMP-parallel pair counter (default: True).
verbose (bool, optional) – Print progress messages (default: True).
angular (bool, optional) – If True, interpret data and rand as equatorial coordinates (ra, dec) in radians, convert to 3D Cartesian unit vectors, and convert rbins from angular to chord distances (default: False).
rng (numpy.random.Generator or None, optional) – Random number generator. If None (default) one is created via numpy.random.default_rng(**kw_rng).
kw_rng (dict, optional) – Keyword arguments forwarded to numpy.random.default_rng (default: {'seed': 555}).

Returns:

xi (ndarray) – Baseline correlation function, shape (rbins.size - 1,).
sigma (ndarray) – Bootstrap standard deviation per bin, same shape as xi.
boots (ndarray, shape (Nboots, rbins.size - 1)) – Individual bootstrap realisations; only returned when return_boots=True.

measure.clustering_core extension module

Pair-counting kernels for two-point statistics.

Provides brute-force \(O(N^2)\) pair counters used by the two-point estimators in scampy.measure.clustering. Pairs are binned in log-spaced separation bins; the _omp variants are provided for API compatibility and are identical to their serial counterparts.

Available functions

2D Cartesian pair counts

`d2D_DD()`, `d2D_DD_omp()`	Data–data pair counts on a 2D Cartesian catalogue.
`d2D_DR()`, `d2D_DR_omp()`	Data–random pair counts on two 2D Cartesian catalogues.

2D angular pair counts

`dA2D_DD()`, `dA2D_DD_omp()`	Data–data pair counts using angular separations.
`dA2D_DR()`, `dA2D_DR_omp()`	Data–random pair counts using angular separations.

3D Cartesian pair counts

`d3D_DD()`, `d3D_DD_omp()`	Data–data pair counts on a 3D Cartesian catalogue.
`d3D_DR()`, `d3D_DR_omp()`	Data–random pair counts on two 3D Cartesian catalogues.

scampy.measure.clustering_core.d2D_DD(X, Y, rbin)

Count data-data pairs in 2D separation bins.

Parameters:

X (array-like of float) – X-coordinates of the catalogue.
Y (array-like of float) – Y-coordinates of the catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.d2D_DR(X1, Y1, X2, Y2, rbin)

Count data-random cross-pairs in 2D separation bins.

Parameters:

X1 (array-like of float) – X-coordinates of the first (data) catalogue.
Y1 (array-like of float) – Y-coordinates of the first catalogue.
X2 (array-like of float) – X-coordinates of the second (random) catalogue.
Y2 (array-like of float) – Y-coordinates of the second catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.dA2D_DD(RA, Dec, thetabin)

Count data-data pairs in 2D angular separation bins.

Parameters:

RA (array-like of float) – Right-ascension coordinates [rad].
Dec (array-like of float) – Declination coordinates [rad].
thetabin (array-like of float) – Angular separation bin edges [rad] (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.dA2D_DR(RA1, Dec1, RA2, Dec2, thetabin)

Count data-random cross-pairs in 2D angular separation bins.

Parameters:

RA1 (array-like of float) – Right-ascension of the first (data) catalogue [rad].
Dec1 (array-like of float) – Declination of the first catalogue [rad].
RA2 (array-like of float) – Right-ascension of the second (random) catalogue [rad].
Dec2 (array-like of float) – Declination of the second catalogue [rad].
thetabin (array-like of float) – Angular separation bin edges [rad] (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.d3D_DD(X, Y, Z, rbin)

Count data-data pairs in 3D separation bins.

Parameters:

X (array-like of float) – X-coordinates of the catalogue.
Y (array-like of float) – Y-coordinates of the catalogue.
Z (array-like of float) – Z-coordinates of the catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.d3D_DR(X1, Y1, Z1, X2, Y2, Z2, rbin)

Count data-random cross-pairs in 3D separation bins.

Parameters:

X1 (array-like of float) – X-coordinates of the first (data) catalogue.
Y1 (array-like of float) – Y-coordinates of the first catalogue.
Z1 (array-like of float) – Z-coordinates of the first catalogue.
X2 (array-like of float) – X-coordinates of the second (random) catalogue.
Y2 (array-like of float) – Y-coordinates of the second catalogue.
Z2 (array-like of float) – Z-coordinates of the second catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.d2D_DD_omp(X, Y, rbin)

Count data-data pairs in 2D separation bins.

Parameters:

X (array-like of float) – X-coordinates of the catalogue.
Y (array-like of float) – Y-coordinates of the catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.d2D_DR_omp(X1, Y1, X2, Y2, rbin)

Count data-random cross-pairs in 2D separation bins.

Parameters:

X1 (array-like of float) – X-coordinates of the first (data) catalogue.
Y1 (array-like of float) – Y-coordinates of the first catalogue.
X2 (array-like of float) – X-coordinates of the second (random) catalogue.
Y2 (array-like of float) – Y-coordinates of the second catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.dA2D_DD_omp(RA, Dec, thetabin)

Count data-data pairs in 2D angular separation bins.

Parameters:

RA (array-like of float) – Right-ascension coordinates [rad].
Dec (array-like of float) – Declination coordinates [rad].
thetabin (array-like of float) – Angular separation bin edges [rad] (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.dA2D_DR_omp(RA1, Dec1, RA2, Dec2, thetabin)

Count data-random cross-pairs in 2D angular separation bins.

Parameters:

RA1 (array-like of float) – Right-ascension of the first (data) catalogue [rad].
Dec1 (array-like of float) – Declination of the first catalogue [rad].
RA2 (array-like of float) – Right-ascension of the second (random) catalogue [rad].
Dec2 (array-like of float) – Declination of the second catalogue [rad].
thetabin (array-like of float) – Angular separation bin edges [rad] (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.d3D_DD_omp(X, Y, Z, rbin)

Count data-data pairs in 3D separation bins.

Parameters:

X (array-like of float) – X-coordinates of the catalogue.
Y (array-like of float) – Y-coordinates of the catalogue.
Z (array-like of float) – Z-coordinates of the catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

scampy.measure.clustering_core.d3D_DR_omp(X1, Y1, Z1, X2, Y2, Z2, rbin)

Count data-random cross-pairs in 3D separation bins.

Parameters:

X1 (array-like of float) – X-coordinates of the first (data) catalogue.
Y1 (array-like of float) – Y-coordinates of the first catalogue.
Z1 (array-like of float) – Z-coordinates of the first catalogue.
X2 (array-like of float) – X-coordinates of the second (random) catalogue.
Y2 (array-like of float) – Y-coordinates of the second catalogue.
Z2 (array-like of float) – Z-coordinates of the second catalogue.
rbin (array-like of float) – Separation bin edges (log-spaced).

Returns:

Pair counts per bin.

Return type:

ndarray of int

utilities subpackage

utilities.base_classes module

Utility dictionary containers for fixed-size array fields.

Provides fixedSizeDict and maskedDict, which enforce that all stored values are numpy arrays of the same fixed length. These classes underpin the catalogue data structures in scampy.catalogue.

class scampy.utilities.base_classes.fixedSizeDict(Nobj, *args, **kwargs)

Dictionary-like container where all values have the same fixed length.

Inherits from collections.abc.MutableMapping. All values are coerced to numpy.ndarray on assignment and must have exactly Nobj elements; attempts to store arrays of a different length raise TypeError. Deletion of keys is not allowed.

Parameters:

Nobj (int) – Fixed length that all stored arrays must have.
*args (MutableMapping, optional) – Initial key-value pairs copied from another mapping.
**kwargs – Additional key-value pairs to store on construction.

Examples

>>> d = fixedSizeDict(3, x=[1, 2, 3], y=[4, 5, 6])
>>> d['z'] = [7, 8, 9]   # OK
>>> d['w'] = [1, 2]      # raises TypeError

keys() → a set-like object providing a view on D's keys

values() → an object providing a view on D's values

items() → a set-like object providing a view on D's items

pop(k[, d]) → v, remove specified key and return the corresponding value.: If key is not found, d is returned if given, otherwise KeyError is raised.

class scampy.utilities.base_classes.maskedDict(Nobj, *args, **kwargs)

Dictionary-like container with a shared boolean mask over all fields.

Extends the concept of fixedSizeDict by storing every value as a numpy.ma.MaskedArray with a single shared hard mask of length Nobj. Masking an element hides it across all fields simultaneously, which is useful for applying selection cuts to catalogue columns without copying data.

The effective length reported by len() is the number of unmasked elements; use the size property for the full (unmasked + masked) length.

Parameters:

Nobj (int) – Total number of elements (masked and unmasked).
*args (MutableMapping, optional) – Initial key-value pairs copied from another mapping.
**kwargs – Additional key-value pairs to store on construction.

Examples

>>> d = maskedDict(4, x=[1, 2, 3, 4])
>>> d.add_mask(numpy.array([False, True, False, False]))
>>> len(d)       # 3 unmasked elements
3
>>> d.size()     # 4 total elements
4

keys() → a set-like object providing a view on D's keys

values() → an object providing a view on D's values

items() → a set-like object providing a view on D's items

pop(k[, d]) → v, remove specified key and return the corresponding value.: If key is not found, d is returned if given, otherwise KeyError is raised.

reset_mask(): Set all mask entries to False, making all elements visible.

add_mask(mask)

Apply an additional boolean mask (logical OR with the current mask).

Parameters:: mask (array-like of bool) – Boolean array of length Nobj; True hides the element.

iter_masked(*keys)

Iterate over unmasked rows for a subset of keys.

Parameters:: *keys (str) – Field names to include in each yielded tuple.
Yields:: tuple – One tuple per unmasked element containing the values of the requested fields in the order given.

utilities.functions module

Mathematical and utility helper functions.

Provides coordinate distance, angle formatting, Fourier-space window functions, trapezoidal integration, rejection sampling, and a simple linear interpolator used throughout the halo-model pipeline.

scampy.utilities.functions.dist_periodic(coord1, coord2, side=None)

Returns the distance in a periodic box

Parameters:

coord1 (array-like or scalar) – position of the 2 sets of coordinates
coord2 (array-like or scalar) – position of the 2 sets of coordinates
side (scalar) – side lenght of the periodic box, if it is not None (default), the periodic distance is computed

Returns:

distance between coord1 and coord2

Return type:

array-like or scalar

scampy.utilities.functions.str_format_deg(theta, rad=False)

String formatting of angles in (degrees-primes-seconds)

Parameters:

theta (scalar float) – angle to be formatted
rad (bool) – (Optional, default = False) whether theta has been provided in radians or not

Returns:

formatted string

Return type:

str

scampy.utilities.functions.FT_tophat(kR)

Fourier transform of a spherical top-hat window function.

\[\tilde{W}(kR) = \frac{3\,[\sin(kR) - kR\cos(kR)]}{(kR)^3}.\]

Parameters:: kR (scalar or array-like) – Dimensionless product of wavenumber and smoothing radius.
Returns:: Window function values, same shape as kR.
Return type:: ndarray

scampy.utilities.functions.FT_tophat_D1(kR)

Derivative of the spherical top-hat window with respect to \(kR\).

\[\tilde{W}'(kR) = \frac{3\,[(kR)^2 - 3]\sin(kR) + 9\,kR\cos(kR)}{(kR)^4}.\]

Used by dsigma2RdR() to evaluate \(\mathrm{d}\sigma^2/\mathrm{d}R\).

Parameters:: kR (scalar or array-like) – Dimensionless product of wavenumber and smoothing radius.
Returns:: Derivative values, same shape as kR.
Return type:: ndarray

scampy.utilities.functions.mod_erf(x, factor=1.0)

Normalised cumulative-like transform of an array via the error function.

Centres x at its midpoint and rescales by factor * std(x), then applies the standard CDF of the normal distribution:

\[\mathrm{mod\_erf}(x) = \frac{1}{2}\left[1 + \mathrm{erf}\! \left(\frac{x - x_M}{\mathrm{factor}\cdot\sigma_x} \right)\right],\]

where \(x_M = (x_\min + x_\max)/2\) and \(\sigma_x = \mathrm{std}(x)\).

Parameters:

x (array-like) – Input array.
factor (float, optional) – Scaling factor applied to the standard deviation (default: 1.0).

Returns:

Values in \([0, 1]\), same shape as x.

Return type:

ndarray

scampy.utilities.functions.truncated_gaussian(mean, std, lower, upper, size, rng=None, kw_rng={'seed': 555})

Draw samples from a Gaussian distribution truncated to [lower, upper).

Internally converts the bounds to the standard-normal parameterisation and delegates to scipy.stats.truncnorm().

Parameters:

mean (float) – Mean of the untruncated Gaussian.
std (float) – Standard deviation of the untruncated Gaussian.
lower (float) – Lower bound of the truncation interval.
upper (float) – Upper bound of the truncation interval (exclusive).
size (int) – Number of samples to draw.
rng (numpy.random.Generator or None, optional) – Random number generator. If None (default) one is created via numpy.random.default_rng(**kw_rng).
kw_rng (dict, optional) – Keyword arguments forwarded to numpy.random.default_rng (default: {'seed': 555}).

Returns:

1-D array of size samples drawn from \(\mathcal{N}(\mu,\sigma^2)\) restricted to [lower, upper).

Return type:

ndarray

scampy.utilities.functions.trap_int(xx, yy)

Trapezoid integration

Parameters:

xx (array-like) – x-domain grid
yy (array-like) – y-domain grid

Returns:

the integral along the whole x-domain

Return type:

float

scampy.utilities.functions.rejection_sampling(pdf, xmin=0.0, xmax=1.0, umin=0.0, umax=1.0, logu=False, size=1, rng=None, verbose=False)

Samples the X random variable from a custom target PDF using rejection sampling. Proposals for the random variable are drawn from a uniform distribution defined between (xmin, xmax). A proposal x_j is accepted if the corresponding value of the proposal PDF u_j is lower than the value of the target PDF computed at the proposal random variable value, i.e. if u_j < PDF(x_j). Values of the proposal PDF are drawn from a uniform distribution defined between (umin, umax). The higher the value of umax with respect to the maximum of the target PDF, the lower the acceptance fraction (but a higher value of umax might produce more reliable samples).

Parameters:

pdf (callable) – the callable function defining the custom PDF (it must have a __call__ method)
xmin (float) – (Optional, default = 0.0) minimum value of the random variable X
xmax (float) – (Optional, default = 1.0) maximum value of the random variable X
umin (float) – (Optional, default = 0.0) minimum value of the proposal PDF
umax (float) – (Optional, default = 1.0) maximum value of the proposal PDF
logu (bool) – (Optional, default = False) if True, umin and umax are interpreted as log10(umin) and log10(umax)
size (int) – (Optional, default = 1) size of the output sample
rng (int or random generator) – (Optional, default = None) either an integer or a random number generator. If an integer is passed, a generator will be instantiated with seed=rng. If None is passed, a random number generator will be instantiated with randomly generated seed.
verbose (bool) – (Optional, default = False) if True prints on stdout informations about the sampling

Returns:

A 1D array with size samples of the random variable extracted from the target PDF

Return type:

1d-array

Note

One can run the Kolmogorov-Smirnov Test to verify that the sampled random variable has a PDF that approximates the target PDF.

scampy.utilities.functions.linear_interpolation(x, xgrid, ygrid)

Interpolates values of some input array on a user defined grid, where the x-domain must be evenly spaced. Specifically, computes:

\[y = y_i + ( x - x_i ) \dfrac{y_{i+1}-y_i}{x_{i+1}-x_i}\]

Parameters:

x (1d-array) – x-values
xg (1d-array) – x-axis of the interpolation grid
yg (1d-array) – y-axis of the interpolation grid

Returns:

y – interpolated y-values

Return type:

1d-array

scampy.utilities.functions.repeated_mask(trues, falses)

Generate a boolean mask with repetitions (equivalent to numpy.repeat with alternating sequences of booleans)

Parameters:

trues (ndarray of int) – number of true values
false (ndarray of int) – number of false values

Returns:

mask

Return type:

ndarray of bool

Examples

>>> from scampy.utilities.functions import repeated_mask
>>> repeated_mask( [1,2], [2,3] )
[ True False False True True False False False ]

utilities.fft module

Fast Hankel transforms on log-spaced grids.

Provides wrappers around scipy.fft.fht() for computing the 3D Fourier transform pair \(P(k) \leftrightarrow \xi(r)\) and the projected (Abel) transform pair \(P(k) \leftrightarrow w(r_p)\), both on log-spaced grids.

scampy.utilities.fft.fstl(xn, yn, lk0=0.0, bias=0.0, mu=0.5)

Fast Spherical Transform on a log-spaced grid (direct).

Evaluates the 3D two-point correlation function from the power spectrum via a Fast Hankel Transform (FHT) of order \(\mu = 1/2\), implementing Eq. 6 of Ronconi et al. (2020):

\[\xi(r) = \frac{1}{2\pi^2} \int_0^\infty k^2\,P(k)\,\frac{\sin(kr)}{kr}\,\mathrm{d}k.\]

The input grid xn must be log-spaced; the output \(r\)-grid is the conjugate grid returned by scipy.fft.fht().

Parameters:

xn (array-like) – Log-spaced wavenumber grid \(k\) in \([h\,\mathrm{Mpc}^{-1}]\).
yn (array-like) – Power spectrum \(P(k)\) values, shape (xn.size,) or (xn.size, Nz) for multiple redshifts.
lk0 (float, optional) – Initial log-offset passed to scipy.fft.fhtoffset() (default: 0.0).
bias (float, optional) – Power-law bias for the FHT (default: 0.0).
mu (float, optional) – Order of the Hankel transform (default: 0.5, corresponding to the spherical Bessel function \(j_0\)).

Returns:

kn (ndarray) – Output \(r\)-grid in \([h^{-1}\,\mathrm{Mpc}]\).
xi (ndarray) – Correlation function \(\xi(r)\), shape (xn.size,) or (Nz, xn.size).

scampy.utilities.fft.fpstl(xn, yn, lk0=0.0, bias=0.0, mu=0.0)

Fast Projected Spherical Transform on a log-spaced grid (direct).

Evaluates the projected two-point correlation function from the power spectrum via a zeroth-order Fast Hankel Transform (FHT), implementing the Abel transform of Eq. 15 of Ronconi et al. (2020):

\[w(r_p) = \frac{1}{2\pi} \int_0^\infty k\,P(k)\,J_0(k\,r_p)\,\mathrm{d}k,\]

where \(J_0\) is the zeroth-order Bessel function of the first kind.

Parameters:

xn (array-like) – Log-spaced wavenumber grid \(k\) in \([h\,\mathrm{Mpc}^{-1}]\).
yn (array-like) – Power spectrum \(P(k)\) values, shape (xn.size,) or (xn.size, Nz) for multiple redshifts.
lk0 (float, optional) – Initial log-offset passed to scipy.fft.fhtoffset() (default: 0.0).
bias (float, optional) – Power-law bias for the FHT (default: 0.0).
mu (float, optional) – Order of the Hankel transform (default: 0.0, corresponding to \(J_0\)).

Returns:

kn (ndarray) – Output \(r_p\)-grid in \([h^{-1}\,\mathrm{Mpc}]\).
wp (ndarray) – Projected correlation function \(w(r_p)\), shape (xn.size,) or (Nz, xn.size).

utilities.constants module

utilities.interpolation extension module

Piecewise-linear interpolators on pre-tabulated grids.

Provides two interpolator classes sharing a common interface.

Available classes

`lin_interp`	Piecewise-linear interpolation in (x, y) space.
`log_interp`	Piecewise-linear interpolation in \((\log_{10} x,\,\log_{10} y)\) space.

Both classes expose the same methods: __call__, get_x, get_y, and integrate.

class scampy.utilities.interpolation.lin_interp(x, y)

Piecewise-linear interpolator built from two equal-length arrays.

Parameters:

x (array-like) – Strictly increasing x-axis values.
y (array-like) – Corresponding y-axis values (same length as x).

get_x(): Return the x-axis array.

get_y(): Return the y-axis array.

integrate(aa, bb)

Integrate the interpolated function over [aa, bb] using the trapezoidal rule.

Parameters:

aa (float) – Lower integration limit.
bb (float) – Upper integration limit.

Returns:

Approximate integral.

Return type:

float

class scampy.utilities.interpolation.log_interp(x, y)

Log-space piecewise-linear interpolator built from two equal-length arrays.

Interpolation is performed in log10(x) vs log10(y) space.

Parameters:

x (array-like) – Strictly increasing positive x-axis values.
y (array-like) – Corresponding positive y-axis values (same length as x).

get_x(): Return the x-axis array.

get_y(): Return the y-axis array.

integrate(aa, bb)

Integrate the interpolated function over [aa, bb] using the trapezoidal rule.

Parameters:

aa (float) – Lower integration limit.
bb (float) – Upper integration limit.

Returns:

Approximate integral.

Return type:

float

Python Documentation

Top-level modules

boxes module

catalogue module

hod module

lightcone module

Paramters

Paramters

power_spectrum module

sham module

cosmology extension module

Available methods of model

halo subpackage

halo.mass_function module

halo.bias module

halo.density_profile module

halo.model module

io subpackage

io.hdf5 module

io.subfind_table module

measure subpackage

measure.abundances module

Parameters:

Returns:

measure.clustering module

measure.clustering_core extension module

Available functions

utilities subpackage

utilities.base_classes module

utilities.functions module

utilities.fft module

utilities.constants module

utilities.interpolation extension module

Available classes

Available methods of `model`