Diemer & Kravtsov 2014 profile

This module implements the Diemer & Kravtsov 2014 form of the density profile. Please see Halo Density Profiles for a general introduction to the Colossus density profile module.

Basics

The DK14 profile (Diemer & Kravtsov 2014) is defined by the following density form:

\[ \begin{align}\begin{aligned}\rho(r) &= \rho_{\rm inner} \times f_{\rm trans} + \rho_{\rm outer}\\\rho_{\rm inner} &= \rho_{\rm Einasto} = \rho_{\rm s} \exp \left( -\frac{2}{\alpha} \left[ \left( \frac{r}{r_{\rm s}} \right)^\alpha -1 \right] \right)\\f_{\rm trans} &= \left[ 1 + \left( \frac{r}{r_{\rm t}} \right)^\beta \right]^{-\frac{\gamma}{\beta}}\end{aligned}\end{align} \]

This profile corresponds to an Einasto profile at small radii, and steepens around the virial radius. The profile formula has 6 free parameters, but most of those can be fixed to particular values that depend on the mass and mass accretion rate of a halo. The parameters have the following meaning:

Param. Symbol Explanation
rhos \(\rho_s\) The central scale density, in physical \(M_{\odot} h^2 / {\rm kpc}^3\)
rs \(r_{\rm s}\) The scale radius in physical kpc/h
alpha \(\alpha\) Determines how quickly the slope of the inner Einasto profile steepens
rt \(r_{\rm t}\) The radius where the profile steepens beyond the Einasto profile, in physical kpc/h
beta \(\beta\) Sharpness of the steepening
gamma \(\gamma\) Asymptotic negative slope of the steepening term

There are two ways to initialize a DK14 profile. First, the user can pass the fundamental parameters of the profile listed above. Second, the user can pass a spherical overdensity mass and concentration, the conversion to the native parameters then relies on the calibrations in DK14. In this case, the user can give additional information about the profile that can be used in setting the fundamental parameters.

In particular, the fitting function was calibrated for the median and mean profiles of two types of halo samples, namely samples selected by mass, and samples selected by both mass and mass accretion rate. The user can choose between those by setting selected_by = 'M' or selected_by = 'Gamma'. The latter option results in a more accurate representation of the density profile, but the mass accretion rate must be known.

If the profile is chosen to model halo samples selected by mass, we set \((\beta, \gamma) = (4, 8)\). If the sample is selected by both mass and mass accretion rate, we set \((\beta, \gamma) = (6, 4)\). Those choices result in a different calibration of the turnover radius \(r_{\rm t}\). In the latter case, both z and Gamma must not be None. See the deriveParameters() function for more details.

Note

The DK14 profile makes sense only if some description of the outer profile is added.

Adding outer terms is easy using the wrapper function getDK14ProfileWithOuterTerms():

from colossus.cosmology import cosmology
from colossus.halo import profile_dk14

cosmology.setCosmology('planck15')
p = profile_dk14.getDK14ProfileWithOuterTerms(M = 1E12, c = 10.0, z = 0.0, mdef = 'vir')

This line will return a DK14 profile object with a power-law outer profile and the mean density of the universe added by default. Alternatively, the user can pass a list of OuterTerm objects (see documentation of the HaloDensityProfile parent class). The user can pass additional parameters to the outer profile terms:

p = profile_dk14.getDK14ProfileWithOuterTerms(M = 1E12, c = 10.0, z = 0.0, mdef = 'vir',
                power_law_slope = 1.2)

or change the outer terms altogether:

p = profile_dk14.getDK14ProfileWithOuterTerms(M = 1E12, c = 10.0, z = 0.0, mdef = 'vir', 
                outer_term_names = ['mean', 'cf'], derive_bias_from = None, bias = 1.2)

Some of the outer term parameterizations (namely the 2-halo term) rely, in turn, on properties of the total profile such as the mass. In those cases, the constructor determines the mass iteratively, taking the changing contribution of the outer term into account. This procedure can make the constructor slow. Thus, it is generally preferred to initialize the outer terms with fixed parameters (e.g., pivot radius or bias). Please see the Tutorials for more code examples.

Module reference

class halo.profile_dk14.DK14Profile(rhos=None, rs=None, rt=None, alpha=None, beta=None, gamma=None, M=None, c=None, mdef=None, z=None, selected_by='M', Gamma=None, outer_terms=[], acc_warn=0.01, acc_err=0.05)

The Diemer & Kravtsov 2014 density profile.

The redshift must always be passed to this constructor, regardless of whether the fundamental parameters or a mass and concentration are given.

Parameters:
rhos: float

The central scale density, in physical \(M_{\odot} h^2 / {\rm kpc}^3\).

rs: float

The scale radius in physical kpc/h.

rt: float

The radius where the profile steepens, in physical kpc/h.

alpha: float

Determines how quickly the slope of the inner Einasto profile steepens.

beta: float

Sharpness of the steepening.

gamma: float

Asymptotic negative slope of the steepening term.

M: float

Halo mass in \(M_{\odot}/h\).

c: float

Concentration in the same mass definition as M.

mdef: str

The mass definition to which M corresponds. See Halo Mass Definitions for details.

z: float

Redshift

selected_by: str

The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma. This parameter influences how some of the fixed parameters in the profile are set, in particular those that describe the steepening term.

Gamma: float

The mass accretion rate as defined in DK14. This parameter only needs to be passed if selected_by == 'Gamma'.

acc_warn: float

If the function achieves a relative accuracy in matching M less than this value, a warning is printed.

acc_err: float

If the function achieves a relative accuracy in matching M less than this value, an exception is raised.

Methods

M4rs() The mass within 4 scale radii, \(M_{<4rs}\).
MDelta(z, mdef) The spherical overdensity mass of a given mass definition.
Msp([search_range]) The mass enclosed within \(R_{\rm sp}\), \(M_{\rm sp}\).
RDelta(z, mdef) The spherical overdensity radius of a given mass definition.
RMDelta(z, mdef) The spherical overdensity radius and mass of a given mass definition.
RMsp([search_range]) The splashback radius and mass within, \(R_{\rm sp}\) and \(M_{\rm sp}\).
Rsp([search_range]) The splashback radius, \(R_{\rm sp}\).
Vmax() The maximum circular velocity, and the radius where it occurs.
circularVelocity(r) The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).
cumulativePdf(r[, Rmax, z, mdef]) The cumulative distribution function of the profile.
deltaSigma(r[, interpolate, …]) The excess surface density at radius r.
deltaSigmaInner(r[, interpolate, …]) The excess surface density at radius r due to the inner profile.
deltaSigmaOuter(r[, interpolate, …]) The excess surface density at radius r due to the outer profile.
density(r) Density as a function of radius.
densityDerivativeLin(r) The linear derivative of density, \(d \rho / dr\).
densityDerivativeLinInner(r) The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).
densityDerivativeLinOuter(r) The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).
densityDerivativeLog(r) The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).
densityDerivativeLogInner(r) The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).
densityDerivativeLogOuter(r) The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).
densityInner(r) Density of the inner profile as a function of radius.
densityOuter(r) Density of the outer profile as a function of radius.
deriveParameters(selected_by[, nu200m, z, Gamma]) Calibration of the parameters \(\alpha\), \(\beta\), \(\gamma\), and \(r_{\rm t}\).
enclosedMass(r[, accuracy]) The mass enclosed within radius r.
enclosedMassInner(r[, accuracy]) The mass enclosed within radius r due to the inner profile term.
enclosedMassOuter(r[, accuracy]) The mass enclosed within radius r due to the outer profile term.
fit(r, q, quantity[, q_err, q_cov, method, …]) Fit the density, mass, or surface density profile to a given set of data points.
getParameterArray([mask]) Returns an array of the profile parameters.
setParameterArray(pars[, mask]) Set the profile parameters from an array.
surfaceDensity(r[, interpolate, accuracy, …]) The projected surface density at radius r.
surfaceDensityInner(r[, interpolate, …]) The projected surface density at radius r due to the inner profile.
surfaceDensityOuter(r[, interpolate, …]) The projected surface density at radius r due to the outer profile.
update() Update the profile options after a parameter change.
static deriveParameters(selected_by, nu200m=None, z=None, Gamma=None)

Calibration of the parameters \(\alpha\), \(\beta\), \(\gamma\), and \(r_{\rm t}\).

This function determines the values of those parameters in the DK14 profile that can be calibrated based on mass, and potentially mass accretion rate. If the profile is chosen to model halo samples selected by mass (selected_by = 'M'), we set \((\beta, \gamma) = (4, 8)\). If the sample is selected by both mass and mass accretion rate (selected_by = 'Gamma'), we set \((\beta, \gamma) = (6, 4)\).

Those choices result in a different calibration of the turnover radius \(r_{\rm t}\). If selected_by = 'M', we use Equation 6 in DK14. Though this relation was originally calibrated for \(\nu = \nu_{\rm vir}\), but the difference is small. If selected_by = 'Gamma', \(r_{\rm t}\) is calibrated from Gamma and z.

Finally, the parameter that determines how quickly the Einasto profile steepens with radius, \(\alpha\), is calibrated according to the Gao et al. 2008 relation. This function was also originally calibrated for \(\nu = \nu_{\rm vir}\), but the difference is small.

Parameters:
selected_by: str

The halo sample to which this profile refers can be selected mass M or by accretion rate Gamma.

nu200m: float

The peak height of the halo for which the parameters are to be calibrated. This parameter only needs to be passed if selected_by == 'M'.

z: float

Redshift

Gamma: float

The mass accretion rate as defined in DK14. This parameter only needs to be passed if selected_by == 'Gamma'.

update()

Update the profile options after a parameter change.

The DK14 profile has one internal option, opt['R200m'], that does not stay in sync with the other profile parameters if they are changed (either inside or outside the constructor). This function adjusts \(R_{\rm 200m}\), in addition to whatever action is taken in the update function of the super class. Note that this adjustment needs to be done iteratively if any outer profiles rely on \(R_{\rm 200m}\).

densityInner(r)

Density of the inner profile as a function of radius.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
density: array_like

Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityDerivativeLinInner(r)

The linear derivative of the inner density, \(d \rho_{\rm inner} / dr\).

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
derivative: array_like

The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

RDelta(z, mdef)

The spherical overdensity radius of a given mass definition.

Parameters:
z: float

Redshift

mdef: str

The mass definition for which the spherical overdensity radius is computed. See Halo Mass Definitions for details.

Returns:
R: float

Spherical overdensity radius in physical kpc/h.

See also

MDelta
The spherical overdensity mass of a given mass definition.
RMDelta
The spherical overdensity radius and mass of a given mass definition.
M4rs()

The mass within 4 scale radii, \(M_{<4rs}\).

This mass definition was suggested by More et al. 2015, see the Advanced mass definition functions section for details.

Returns:
M4rs: float

The mass within 4 scale radii, \(M_{<4rs}\), in \(M_{\odot} / h\).

Rsp(search_range=5.0)

The splashback radius, \(R_{\rm sp}\).

See the Splashback radius section for a detailed description of the splashback radius. Here, we define \(R_{\rm sp}\) as the radius where the profile reaches its steepest logarithmic slope.

Parameters:
search_range: float

When searching for the radius of steepest slope, search within this factor of \(R_{\rm 200m}\) (optional).

Returns:
Rsp: float

The splashback radius, \(R_{\rm sp}\), in physical kpc/h.

See also

RMsp
The splashback radius and mass within, \(R_{\rm sp}\) and \(M_{\rm sp}\).
Msp
The mass enclosed within \(R_{\rm sp}\), \(M_{\rm sp}\).
RMsp(search_range=5.0)

The splashback radius and mass within, \(R_{\rm sp}\) and \(M_{\rm sp}\).

See the Splashback radius section for a detailed description of the splashback radius. Here, we define \(R_{\rm sp}\) as the radius where the profile reaches its steepest logarithmic slope.

Parameters:
search_range: float

When searching for the radius of steepest slope, search within this factor of \(R_{\rm 200m}\) (optional).

Returns:
Rsp: float

The splashback radius, \(R_{\rm sp}\), in physical kpc/h.

Msp: float

The mass enclosed within the splashback radius, \(M_{\rm sp}\), in \(M_{\odot} / h\).

See also

Rsp
The splashback radius, \(R_{\rm sp}\).
Msp
The mass enclosed within \(R_{\rm sp}\), \(M_{\rm sp}\).
Msp(search_range=5.0)

The mass enclosed within \(R_{\rm sp}\), \(M_{\rm sp}\).

See the Splashback radius section for a detailed description of the splashback radius. Here, we define \(R_{\rm sp}\) as the radius where the profile reaches its steepest logarithmic slope.

Parameters:
search_range: float

When searching for the radius of steepest slope, search within this factor of \(R_{\rm 200m}\) (optional).

Returns:
Msp: float

The mass enclosed within the splashback radius, \(M_{\rm sp}\), in \(M_{\odot} / h\).

See also

Rsp
The splashback radius, \(R_{\rm sp}\).
RMsp
The splashback radius and mass within, \(R_{\rm sp}\) and \(M_{\rm sp}\).
MDelta(z, mdef)

The spherical overdensity mass of a given mass definition.

Parameters:
z: float

Redshift

mdef: str

The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns:
M: float

Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta
The spherical overdensity radius of a given mass definition.
RMDelta
The spherical overdensity radius and mass of a given mass definition.
RMDelta(z, mdef)

The spherical overdensity radius and mass of a given mass definition.

This is a wrapper for the RDelta() and MDelta() functions which returns both radius and mass.

Parameters:
z: float

Redshift

mdef: str

The mass definition for which the spherical overdensity mass is computed. See Halo Mass Definitions for details.

Returns:
R: float

Spherical overdensity radius in physical kpc/h.

M: float

Spherical overdensity mass in \(M_{\odot} /h\).

See also

RDelta
The spherical overdensity radius of a given mass definition.
MDelta
The spherical overdensity mass of a given mass definition.
Vmax()

The maximum circular velocity, and the radius where it occurs.

Returns:
vmax: float

The maximum circular velocity in km / s.

rmax: float

The radius where fmax occurs, in physical kpc/h.

See also

circularVelocity
The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).
circularVelocity(r)

The circular velocity, \(v_c \equiv \sqrt{GM(<r)/r}\).

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
vc: float

The circular velocity in km / s; has the same dimensions as r.

See also

Vmax
The maximum circular velocity, and the radius where it occurs.
cumulativePdf(r, Rmax=None, z=None, mdef=None)

The cumulative distribution function of the profile.

Some density profiles do not converge to a finite mass at large radius, and the distribution thus needs to be cut off. The user can specify either a radius (in physical kpc/h) where the profile is cut off, or a mass definition and redshift to compute this radius (e.g., the virial radius \(R_{vir}\) at z = 0).

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Rmax: float

The radius where to cut off the profile in physical kpc/h.

z: float

Redshift

mdef: str

The radius definition for the cut-off radius. See Halo Mass Definitions for details.

Returns:
pdf: array_like

The probability for mass to lie within radius r; has the same dimensions as r.

deltaSigma(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r.

This quantity is useful in weak lensing studies, and is defined as \(\Delta\Sigma(R) = \Sigma(<R)-\Sigma(R)\) where \(\Sigma(<R)\) is the averaged surface density within R weighted by area,

\[\Delta\Sigma(R) = \frac{1}{\pi R^2} \int_0^{R} 2 \pi r \Sigma(r) dr - \Sigma(R)\]
Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

interpolate: bool

Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.

interpolate_surface_density: bool

Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.

accuracy: float

The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).

min_r_interpolate: float

The minimum radius in physical kpc/h from which the surface density profile is averaged.

max_r_interpolate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.

max_r_integrate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns:
DeltaSigma: array_like

The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaInner(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the inner profile.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

interpolate: bool

Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.

interpolate_surface_density: bool

Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.

accuracy: float

The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).

min_r_interpolate: float

The minimum radius in physical kpc/h from which the surface density profile is averaged.

max_r_interpolate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.

max_r_integrate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns:
DeltaSigma: array_like

The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

deltaSigmaOuter(r, interpolate=True, interpolate_surface_density=True, accuracy=0.0001, min_r_interpolate=1e-06, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The excess surface density at radius r due to the outer profile.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

interpolate: bool

Use an interpolation table for the surface density during the integration. This can speed up the evaluation significantly, as the surface density can be expensive to evaluate.

interpolate_surface_density: bool

Use an interpolation table for density during the computation of the surface density. This should make the evaluation somewhat faster, but can fail for some density terms which are negative at particular radii.

accuracy: float

The minimum accuracy of the integration (used both to compute the surface density and average it to get DeltaSigma).

min_r_interpolate: float

The minimum radius in physical kpc/h from which the surface density profile is averaged.

max_r_interpolate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.

max_r_integrate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns:
DeltaSigma: array_like

The excess surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

density(r)

Density as a function of radius.

Abstract function which must be overwritten by child classes.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
density: array_like

Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

densityDerivativeLin(r)

The linear derivative of density, \(d \rho / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer terms.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
derivative: array_like

The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLinOuter(r)

The linear derivative of the outer density, \(d \rho_{\rm outer} / dr\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
derivative: array_like

The linear derivative in physical \(M_{\odot} h / {\rm kpc}^2\); has the same dimensions as r.

densityDerivativeLog(r)

The logarithmic derivative of density, \(d \log(\rho) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from the inner and outer profile terms.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
derivative: array_like

The dimensionless logarithmic derivative; has the same dimensions as r.

densityDerivativeLogInner(r)

The logarithmic derivative of the inner density, \(d \log(\rho_{\rm inner}) / d \log(r)\).

This function evaluates the logarithmic derivative based on the linear derivative. If there is an analytic expression for the logarithmic derivative, child classes should overwrite this function.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
derivative: array_like

The dimensionless logarithmic derivative; has the same dimensions as r.

densityDerivativeLogOuter(r)

The logarithmic derivative of the outer density, \(d \log(\rho_{\rm outer}) / d \log(r)\).

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from outer profile terms.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
derivative: array_like

The dimensionless logarithmic derivative; has the same dimensions as r.

densityOuter(r)

Density of the outer profile as a function of radius.

This function should generally not be overwritten by child classes since it handles the general case of adding up the contributions from all outer profile terms.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

Returns:
density: array_like

Density in physical \(M_{\odot} h^2 / {\rm kpc}^3\); has the same dimensions as r.

enclosedMass(r, accuracy=1e-06)

The mass enclosed within radius r.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

accuracy: float

The minimum accuracy of the integration.

Returns:
M: array_like

The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassInner(r, accuracy=1e-06)

The mass enclosed within radius r due to the inner profile term.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

accuracy: float

The minimum accuracy of the integration.

Returns:
M: array_like

The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

enclosedMassOuter(r, accuracy=1e-06)

The mass enclosed within radius r due to the outer profile term.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

accuracy: float

The minimum accuracy of the integration.

Returns:
M: array_like

The mass enclosed within radius r, in \(M_{\odot}/h\); has the same dimensions as r.

fit(r, q, quantity, q_err=None, q_cov=None, method='leastsq', mask=None, verbose=True, tolerance=1e-05, maxfev=0, initial_step=0.1, nwalkers=100, random_seed=None, convergence_step=100, converged_GR=0.01, best_fit='median', output_every_n=100)

Fit the density, mass, or surface density profile to a given set of data points.

This function represents a general interface for finding the best-fit parameters of a halo density profile given a set of data points. These points can represent a number of different physical quantities: quantity can either be density, enclosed mass, surface density, or Delta Sigma (rho, M, Sigma, or DeltaSigma).

The data points q at radii r can optionally have error bars, and the user can pass a full covariance matrix. Please note that not passing any estimate of the uncertainty, i.e. q_err = None and q_cov = None, can lead to very poor fit results: the fitter will minimize the absolute difference between points, strongly favoring the high densities at the center.

There are two fundamental methods for performing the fit, a least-squares minimization (method = 'leastsq') and a Markov-Chain Monte Carlo (method = 'mcmc'). The MCMC method has some specific options (see below). In either case, the current parameters of the profile instance serve as an initial guess. Finally, the user can choose to vary only a sub-set of the profile parameters through the mask parameter.

The function returns a dictionary with outputs that depend on which method is chosen. After this function has completed, the profile instance represents the best-fit profile to the data points (i.e., its parameters are the best-fit parameters). Note that all output parameters are bundled into one dictionary. The explanations below refer to the entries in this dictionary.

Parameters:
r: array_like

The radii of the data points, in physical kpc/h.

q: array_like

The data to fit; can either be density in physical \(M_{\odot} h^2 / {\rm kpc}^3\), enclosed mass in \(M_{\odot} /h\), or surface density in physical \(M_{\odot} h/{\rm kpc}^2\). Must have the same dimensions as r.

quantity: str

Indicates which quantity is given in as input in q, can be rho, M, Sigma, or DeltaSigma.

q_err: array_like

Optional; the uncertainty on the values in q in the same units. If method == 'mcmc', either q_err or q_cov must be passed. If method == 'leastsq' and neither q_err nor q_cov are passed, the absolute different between data points and fit is minimized. In this case, the returned chi2 is in units of absolute difference, meaning its value will depend on the units of q.

q_cov: array_like

Optional; the covariance matrix of the elements in q, as a 2-dimensional numpy array. This array must have dimensions of q**2 and be in units of the square of the units of q. If q_cov is passed, q_err is ignored since the diagonal elements of q_cov correspond to q_err**2.

method: str

The fitting method; can be leastsq for a least-squares minimization of mcmc for a Markov-Chain Monte Carlo.

mask: array_like

Optional; a numpy array of booleans that has the same length as the variables vector of the density profile class. Only variables where mask == True are varied in the fit, all others are kept constant. Important: this has to be a numpy array rather than a list.

verbose: bool / int

If true, output information about the fitting process. The flag can also be set as a number, where 1 has the same effect as True, and 2 outputs large amounts of information such as the fit parameters at each iteration.

tolerance: float

Only active when method == 'leastsq'. The accuracy to which the best-fit parameters are found.

maxfev: int

Only active when method == 'leastsq'. The maximum number of function evaluations before the fit is aborted. If zero, the default value of the scipy leastsq function is used.

initial_step: array_like

Only active when method == 'mcmc'. The MCMC samples (“walkers”) are initially distributed in a Gaussian around the initial guess. The width of the Gaussian is given by initial_step, either as an array of length N_par (giving the width of each Gaussian) or as a float number, in which case the width is set to initial_step times the initial value of the parameter.

nwalkers: int

Only active when method == 'mcmc'. The number of MCMC samplers that are run in parallel.

random_seed: int

Only active when method == 'mcmc'. If random_seed is not None, it is used to initialize the random number generator. This can be useful for reproducing results.

convergence_step: int

Only active when method == 'mcmc'. The convergence criteria are computed every convergence_step steps (and output is printed if verbose == True).

converged_GR: float

Only active when method == 'mcmc'. The maximum difference between different chains, according to the Gelman-Rubin criterion. Once the GR indicator is lower than this number in all parameters, the chain is ended. Setting this number too low leads to very long runtimes, but setting it too high can lead to inaccurate results.

best_fit: str

Only active when method == 'mcmc'. This parameter determines whether the mean or median value of the likelihood distribution is used as the output parameter set.

output_every_n: int

Only active when method == 'mcmc'. This parameter determines how frequently the MCMC chain outputs information. Only effective if verbose == True.

Returns:
results: dict

A dictionary bundling the various fit results. Regardless of the fitting method, the dictionary always contains the following entries:

x: array_like

The best-fit result vector. If mask is passed, this vector only contains those variables that were varied in the fit.

q_fit: array_like

The fitted profile at the radii of the data points; has the same units as q and the same dimensions as r.

chi2: float

The chi^2 of the best-fit profile. If a covariance matrix was passed, the covariances are taken into account. If no uncertainty was passed at all, chi2 is in units of absolute difference, meaning its value will depend on the units of q.

ndof: int

The number of degrees of freedom, i.e. the number of fitted data points minus the number of free parameters.

chi2_ndof: float

The chi^2 per degree of freedom.

If method == 'leastsq', the dictionary additionally contains the entries returned by scipy.optimize.leastsq as well as the following:

nfev: int

The number of function calls used in the fit.

x_err: array_like

An array of dimensions [2, nparams] which contains an estimate of the lower and upper uncertainties on the fitted parameters. These uncertainties are computed from the covariance matrix estimated by the fitter. Please note that this estimate does not exactly correspond to a 68% likelihood. In order to get more statistically meaningful uncertainties, please use the MCMC samples instead of least-squares. In some cases, the fitter fails to return a covariance matrix, in which case x_err is None.

If method == 'mcmc', the dictionary contains the following entries:

x_initial: array_like

The initial positions of the walkers, in an array of dimensions [nwalkers, nparams].

chain_full: array_like

A numpy array of dimensions [n_independent_samples, nparams] with the parameters at each step in the chain. In this thin chain, only every nth step is output, where n is the auto-correlation time, meaning that the samples in this chain are truly independent.

chain_thin: array_like

Like the thin chain, but including all steps. Thus, the samples in this chain are not indepedent from each other. However, the full chain often gives better plotting results.

R: array_like

A numpy array containing the GR indicator at each step when it was saved.

x_mean: array_like

The mean of the chain for each parameter; has length nparams.

x_median: array_like

The median of the chain for each parameter; has length nparams.

x_stddev: array_like

The standard deviation of the chain for each parameter; has length nparams.

x_percentiles: array_like

The lower and upper values of each parameter that contain a certain percentile of the probability; has dimensions [n_percentages, 2, nparams] where the second dimension contains the lower/upper values.

getParameterArray(mask=None)

Returns an array of the profile parameters.

The profile parameters are internally stored in an ordered dictionary. For some applications (e.g., fitting), a simply array is more appropriate.

Parameters:
mask: array_like

Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). Only those parameters that correspond to True values are returned.

Returns:
par: array_like

A numpy array with the profile’s parameter values.

setParameterArray(pars, mask=None)

Set the profile parameters from an array.

The profile parameters are internally stored in an ordered dictionary. For some applications (e.g., fitting), setting them directly from an array might be necessary. If the profile contains values that depend on the parameters, the profile class must overwrite this function and update according to the new parameters.

Parameters:
pars: array_like

The new parameter array.

mask: array_like

Optional; must be a numpy array (not a list) of booleans, with the same length as the parameter vector of the profile class (profile.N_par). If passed, only those parameters that correspond to True values are set (meaning the pars parameter must be shorter than profile.N_par).

surfaceDensity(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r.

The surface density is computed by projecting the 3D density along the line of sight,

\[\Sigma(R) = 2 \int_R^{\infty} \frac{r \rho(r)}{\sqrt{r^2-R^2}} dr\]
Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

interpolate: bool

Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.

accuracy: float

The minimum accuracy of the integration.

max_r_interpolate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.

max_r_integrate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns:
Sigma: array_like

The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityInner(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the inner profile.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

interpolate: bool

Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.

accuracy: float

The minimum accuracy of the integration.

max_r_interpolate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.

max_r_integrate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns:
Sigma: array_like

The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

surfaceDensityOuter(r, interpolate=True, accuracy=0.0001, max_r_interpolate=100000000.0, max_r_integrate=1e+20)

The projected surface density at radius r due to the outer profile.

This function checks whether there are explicit expressions for the surface density of the outer profile terms available, and uses them if possible. Note that there are some outer terms whose surface density integrates to infinity, such as the mean density of the universe which is constant to infinitely large radii.

Parameters:
r: array_like

Radius in physical kpc/h; can be a number or a numpy array.

interpolate: bool

Use an interpolation table for density during the integration. This should make the evaluation somewhat faster, depending on how large the radius array is.

accuracy: float

The minimum accuracy of the integration.

max_r_interpolate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using interpolating density.

max_r_integrate: float

The maximum radius in physical kpc/h to which the density profile is integrated when using exact densities.

Returns:
Sigma: array_like

The surface density at radius r, in physical \(M_{\odot} h/{\rm kpc}^2\); has the same dimensions as r.

halo.profile_dk14.getDK14ProfileWithOuterTerms(outer_term_names=['mean', 'pl'], power_law_norm=1.0, power_law_slope=1.5, power_law_max=1000.0, derive_bias_from='R200m', bias=1.0, **kwargs)

A wrapper function to create a DK14 profile with one or many outer profile terms.

The DK14 profile only makes sense if some description of the outer profile is added. This function provides a convenient way to construct such profiles without having to set the properties of the outer terms manually. Valid keys for outer terms include the following.

  • mean: The mean density of the universe at redshift z (see the documentation of OuterTermMeanDensity).
  • pl: A power-law profile in radius (see the documentation of OuterTermPowerLaw). For the DK14 profile, the chosen pivot radius is \(5 R_{\rm 200m}\). Note that \(R_{\rm 200m}\) is set as a profile option in the constructor once, but not adjusted thereafter unless the update() function is called. Thus, in a fit, the fitted norm and slope refer to a pivot of the original \(R_{\rm 200m}\) until update() is called which adjusts these parameters. Furthermore, the parameters for the power-law outer profile (norm and slope, called \(b_{\rm e}\) and \(s_{\rm e}\) in the DK14 paper) exhibit a complicated dependence on halo mass, redshift and cosmology. At low redshift, and for the cosmology considered in DK14, power_law_norm = 1.0 and power_law_slope = 1.5 are reasonable values over a wide range of masses (see Figure 18 in DK14), but these values are by no means universal or accurate.
  • cf: The matter-matter correlation function times halo bias (see the documentation of OuterTermCorrelationFunction). Here, the user has a choice regarding halo bias: it can enter the profile as a parameter (if derive_bias_from == None or it can be derived according to the default model of halo bias based on \(M_{\rm 200m}\) (in which case derive_bias_from = 'R200m' and the bias parameter is ignored). The latter option can make the constructor slow because of the iterative evaluation of bias and \(M_{\rm 200m}\).
Parameters:
outer_term_names: array_like

A list of outer profile term identifiers, can be mean, pl, or cf.

power_law_norm: float

The normalization of a power-law term (called \(b_{\rm e}\) in DK14).

power_law_slope: float

The negative slope of a power-law term (called \(s_{\rm e}\) in DK14).

power_law_max: float

The maximum density contributed by a power-law term.

derive_bias_from: str

See cf section above.

bias: float

See cf section above.

kwargs: kwargs

The arguments passed to the DK14 profile constructor (i.e., the fundamental parameters or M, c etc).