"""
A module containing various smoothing filter Component models,
including the popular top-hat in real space.
"""
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline as _spline
import scipy.integrate as intg
import collections
from .._internals import _framework, _utils
import warnings
[docs]@_framework.pluggable
class Filter(_framework.Component):
r"""
Base class for Filter components.
Filters handle the calculation of the mass variance from the power
spectrum, via a window function. Subclasses of :class:`~Filer` implement
specific window functions.
The general design is to specify all quantities in terms of length scales,
rather than equivalent masses, but conversion methods are provided.
Any explicit subclass need only specify `k_space`, `mass_to_radius`,
`radius_to_mass` and `dw_dlnkr`, with `dlnr_dlnm` optional if the
mass assignment is unconventional.
Parameters
----------
k : array_like
Wavenumbers at which the power spectrum is defined.
power : array_like
The power spectrum at `k`.
\*\*model_parameters : unpacked-dict
As for any :class:`hmf._framework.Component` subclass, any particular
parameters of the model may be passed to the constructor. Allowed
parameters are found in the :attr:`~._defaults` attribute.
Notes
-----
Besides the raw filter function itself, two quantities are of primary interest:
firstly the mass variance (see :meth:`sigma`), which appears in many cosmological
applications, and secondly its logarithmic derivative with mass, which appears in
the Press-Schechter formalism for the halo mass function.
To remain extensible and general, the methodology in this class is to calculate the
latter quantity as
.. math:: \frac{d\ln\sigma}{d\ln m} = \frac{1}{2} \frac{d \ln \sigma^2}{d\ln R}\frac{d\ln R}{d\ln m}.
Each of the quantities on the right can be separately calculated, improving
extensibility.
The factor :math:`\frac{d\ln R}{d\ln m}` is typically 1/3, but this is not
necessarily the case for window functions of arbitrary shape.
"""
[docs] def __init__(self, k, power, **model_parameters):
self.k = k
self.power = power
super(Filter, self).__init__(**model_parameters)
[docs] def real_space(self, R, r):
r"""
Filter definition in real space.
Parameters
----------
R : float
The smoothing scale
r : array_like
The radial co-ordinate
"""
pass
[docs] def k_space(self, kr):
r"""
Fourier-transform of the real-space filter.
Parameters
----------
kr : array_like
The scales at which to return the filter function
Returns
-------
w : array_like
The filter in fourier space, ``len(kr)``
"""
pass
[docs] def mass_to_radius(self, m, rho_mean):
r"""
Return radius of a region of space given its mass.
Parameters
----------
m : array_like
Masses
rho_mean : float
Mean density of the Universe.
Returns
-------
r : array_like
The corresponding radii to m
Notes
-----
The units of `m` don't matter as long as they are consistent with `rho_mean`.
"""
pass
[docs] def radius_to_mass(self, r, rho_mean):
r"""
Return mass of a region of space given its radius
Parameters
----------
r : array_like
Radii
rho_mean : float
Mean density of the Universe.
Returns
-------
m : float or array of floats
The corresponding masses to r
Notes
-----
The units of `r` don't matter as long as they are consistent with `rho_mean`.
"""
pass
[docs] def dw_dlnkr(self, kr):
r"""
The derivative of the (fourier-transformed) filter with :math:`\ln(kr)`.
Parameters
----------
kr : array_like
Scale(s) at which the derivative is evaluated.
Notes
-----
In terms of :math:`\frac{dw^2}{dm}`, which is a commonly used
quantity, this has the relationship
.. math:: w\frac{dw}{d\ln r} = \frac{2}{r}\frac{dw^2}{dm}\frac{dm}{dr}.
"""
pass
[docs] def dlnss_dlnr(self, r):
r"""
The derivative of the mass variance with radius.
Parameters
----------
r : array_like
Radii
Returns
-------
dlnss_dlnr : array_like
The derivative of the the mass variance with radius.
Notes
-----
Given a prescription for how radius grows with mass (typically with
a log-slope of 1/3, and set in :meth:`dlnr_dlnm`), this specifies
the quantity :math:`\frac{d \ln \sigma^2}{d\ln m}`.
The general formula is
.. math:: \frac{d\ln \sigma^2}{d\ln R} = \frac{1}{\pi^2\sigma^2} \int_0^\infty W(kR) \frac{dW(kR)}{d\ln(kR)} P(k)k^2 dk
"""
dlnk = np.log(self.k[1] / self.k[0])
s = self.sigma(r)
rk = np.outer(r, self.k)
rest = self.power * self.k ** 3
w = self.k_space(rk)
dw = self.dw_dlnkr(rk)
integ = w * dw * rest
return intg.simps(integ, dx=dlnk, axis=-1) / (np.pi ** 2 * s ** 2)
[docs] def dlnr_dlnm(self, r):
r"""
The derivative of log radius with log mass.
For the usual :math:`m\propto r^3` mass assignment, this is just 1/3.
Parameters
----------
r : array_like
Radii.
"""
return 1.0 / 3.0
[docs] def dlnss_dlnm(self, r):
r"""
The logarithmic slope of mass variance with mass.
This is an important quantity, and is used directly to calculate
:math:`\frac{dn}{dm}`.
Parameters
----------
r : array_like
Radii.
"""
return self.dlnss_dlnr(r) * self.dlnr_dlnm(r)
[docs] def sigma(self, r, order=0, rk=None):
r"""
Calculate the nth-moment of the smoothed density field, :math:`\sigma_n(r)`.
.. note :: This is not :math:`\sigma_n^2(r)`!
Parameters
----------
r : float or array_like
The radii of the spheres at which to calculate the nth moment.
order : int, optional
The order of the moment. Default 0 corresponds to common
mass variance.
Returns
-------
sigma : array_like
The square root of the moment at `r`.
Notes
-----
The general definition for the nth-moment of the smoothed density field is (see Bardeen et al. 1986, Eq 4.6c)
.. math:: \sigma^2_n(R) = \frac{1}{2\pi^2} \int_0^\infty dk\ k^{2(1+n)} P(k) W^2(kR)
"""
if rk is None:
rk = np.outer(r, self.k)
dlnk = np.log(self.k[1] / self.k[0])
# we multiply by k because our steps are in logk.
rest = self.power * self.k ** (3 + order * 2)
integ = rest * self.k_space(rk) ** 2
sigma = (0.5 / np.pi ** 2) * intg.simps(integ, dx=dlnk, axis=-1)
return np.sqrt(sigma)
[docs] def nu(self, r, delta_c=1.686):
r"""
Peak height, :math:`\frac{\delta_c^2}{\sigma^2(r)}`.
Parameters
----------
r : array_like
Radii
delta_c : float, optional
Critical overdensity for collapse.
"""
return (delta_c / self.sigma(r)) ** 2
[docs]@_utils.inherit_docstrings
class TopHat(Filter):
r"""
Real-space top-hat window function.
This class is based on :class:`~Filter`, which can be consulted for
details of how to instantiate it.
Notes
-----
This filter specifically implements the real-space filter:
.. math:: F(r) = H(R-r)
for a filter size of *R*, where *H* is the Heaviside step-function.
Its fourier-transform is
.. math:: W(x=kR) = 3\frac{\sin x - x\cos x}{x^3}.
Furthermore the mass-assignment is
.. math:: m(R) = \frac{4\pi}{3} R^3 \bar{\rho}
and the derivative of the window function is
.. math:: \frac{dW}{d\ln x}(x=kR) = \frac{1}{x^3}[9x\cos x + 3(x^2-3)\sin x].
"""
[docs] def real_space(self, R, r):
a = np.where(r < R, 1, 0)
return np.where(r == R, 0.5, a)
[docs] def k_space(self, kr):
return np.where(kr > 1.4e-6, (3 / kr ** 3) * (np.sin(kr) - kr * np.cos(kr)), 1)
[docs] def mass_to_radius(self, m, rho_mean):
return (3.0 * m / (4.0 * np.pi * rho_mean)) ** (1.0 / 3.0)
[docs] def radius_to_mass(self, r, rho_mean):
return 4 * np.pi * r ** 3 * rho_mean / 3
[docs] def dw_dlnkr(self, kr):
return np.where(
kr > 1e-3,
(9 * kr * np.cos(kr) + 3 * (kr ** 2 - 3) * np.sin(kr)) / kr ** 3,
0,
)
[docs]@_utils.inherit_docstrings
class Gaussian(Filter):
r"""
Gaussian window function.
This class is based on :class:`~Filter`, which can be consulted for
details of how to instantiate it.
Notes
-----
The real-space filter is
.. math:: F(r) = \frac{\exp(-r^2/2R^2)}{R^3 (2\pi)^{3/2}}
for a filter scale of *R*.
The fourier-transform of the filter is
.. math:: W(x=kR) = \exp(-x^2/2).
The mass-assignment is
.. math:: m(R) = R^3(2\pi)^{3/2}\bar{\rho},
and the derivative of the window function is
.. math:: \frac{dW}{d\ln x}(x=kR) = -xW(x).
"""
[docs] def real_space(self, R, r):
return np.exp(-(r ** 2) / 2 / R ** 2) / (2 * np.pi) ** 1.5 / R ** 3
[docs] def k_space(self, kr):
return np.exp(-(kr ** 2) / 2.0)
[docs] def mass_to_radius(self, m, rho_mean):
return (m / rho_mean) ** (1.0 / 3.0) / np.sqrt(2 * np.pi)
[docs] def radius_to_mass(self, r, rho_mean):
return (2 * np.pi) ** 1.5 * r ** 3 * rho_mean
[docs] def dw_dlnkr(self, kr):
return -(kr ** 2) * self.k_space(kr)
[docs]@_utils.inherit_docstrings
class SharpK(Filter):
r"""
Fourier-space top-hat window function
This class is based on :class:`~Filter`, which can be consulted for
details of how to instantiate it.
Notes
-----
The real-space filter is
.. math:: F(r) = \frac{\sin(r/R) -(r/R)\cos(r/R)}{2\pi^2r^3},
for a filter scale of *R*.
The fourier-transform of the filter is
.. math:: W(x=kR) = H(kR-1)
where *H* is the Heaviside step-function. The mass-assignment is
.. math:: m(R) = \frac{4\pi}{3}[cR]^3\bar{\rho},
where *c* is a free parameter, typically c~2.5. The derivative of the window function is
.. math:: \frac{dW}{d\ln x}(x=kR) = \delta_D(x-1),
where :math:`\delta_D` is the Dirac delta. Furthermore, the derivative of the mass variance
takes a particularly simple form in this filter:
.. math:: \frac{d\ln \sigma^2}{d \ln R} = -\frac{P(1/R)}{2\pi^2\sigma^2(R)R^3}.
"""
_defaults = {"c": 2.5}
[docs] def k_space(self, kr):
a = np.where(kr > 1, 0, 1)
return np.where(kr == 1, 0.5, a)
[docs] def real_space(self, R, r):
return (np.sin(r / R) - (r / R) * np.cos(r / R)) / (2 * np.pi ** 2 * r ** 3)
[docs] def dw_dlnkr(self, kr):
return np.where(kr == 1, 1.0, 0.0)
[docs] def dlnss_dlnr(self, r):
sigma = self.sigma(r)
power = _spline(self.k, self.power)(1 / r)
return -power / (2 * np.pi ** 2 * sigma ** 2 * r ** 3)
[docs] def mass_to_radius(self, m, rho_mean):
return (1.0 / self.params["c"]) * (3.0 * m / (4.0 * np.pi * rho_mean)) ** (
1.0 / 3.0
)
[docs] def radius_to_mass(self, r, rho_mean):
return 4 * np.pi * (self.params["c"] * r) ** 3 * rho_mean / 3
[docs] def sigma(self, r, order=0):
if not isinstance(r, collections.Iterable):
r = np.atleast_1d(r)
if self.k.max() < 1 / r.min():
warnings.warn("Warning: Maximum r*k less than 1!")
# # Need to re-define this because the integral needs to go exactly kr=1
# # or else the function 'jitters'
sigma = np.zeros(len(r))
power = _spline(self.k, self.power)
for i, rr in enumerate(r):
k = np.logspace(
np.log10(self.k[0]),
min(np.log10(self.k.max()), np.log10(1.0 / rr)),
max(100, len(self.k) - i),
)
p = power(k)
dlnk = np.log(k[1] / k[0])
integ = p * k ** (3 + 2 * order)
sigma[i] = (0.5 / (np.pi ** 2)) * intg.simps(integ, dx=dlnk)
return np.sqrt(sigma)
[docs]@_utils.inherit_docstrings
class SharpKEllipsoid(SharpK):
"""
Fourier-space top-hat window function with ellipsoidal correction
See Schneider, Smith, Reed 2013.
Refer to :class:`~Filter` for more details.
"""
_defaults = {"c": 2.0}
[docs] def xm(self, g, v):
"""
Peak of the distribution of x, where x is the sum of the eigenvalues
of the inertia tensor (?) of an ellipsoidal peak, divided by the second
spectral moment.
Equation A6. in Schneider et al. 2013
"""
top = 3 * (1 - g ** 2) + (1.1 - 0.9 * g ** 4) * np.exp(
-g * (1 - g ** 2) * (g * v / 2) ** 2
)
bot = (3 * (1 - g ** 2) + 0.45 + (g * v / 2) ** 2) ** 0.5 + g * v / 2
return g * v + top / bot
[docs] def em(self, xm):
"""
The average ellipticity of a patch as a function of peak tensor
"""
return 1 / np.sqrt(5 * xm ** 2 + 6)
[docs] def pm(self, xm):
"""
The average prolateness of a patch as a function of peak tensor
"""
return 30.0 / (5 * xm ** 2 + 6) ** 2
[docs] def a3a1(self, e, p):
"""
The short to long axis ratio of an ellipsoid given its ellipticity and
prolateness
"""
return np.sqrt((1 - 3 * e + p) / (1 + 3 * e + p))
[docs] def a3a2(self, e, p):
"""
The short to medium axis ratio of an ellipsoid given its ellipticity and
prolateness
"""
return np.sqrt((1 - 2 * p) / (1 + 3 * e + p))
[docs] def gamma(self, r):
"""
Bardeen et al. 1986 equation 6.17
"""
sig_0 = self.sigma(r)
sig_1 = self.sigma(r, order=1)
sig_2 = self.sigma(r, order=2)
return sig_1 ** 2 / (sig_0 * sig_2)
[docs] def xi(self, pm, em):
return ((1 + 4 * pm) ** 2 / (1 - 3 * em + pm) / (1 - 2 * pm)) ** (1.0 / 6.0)
[docs] def a3(self, r):
g = self.gamma(r)
xm = self.xm(g, self.nu(r))
em = self.em(xm)
pm = self.pm(xm)
return r / self.xi(pm, em)
[docs] def r_a3(self, rmin, rmax):
r = np.logspace(np.log(rmin), np.log(rmax), 200, base=np.e)
a3 = self.a3(r)
s = _spline(a3, r)
return s
[docs] def dlnss_dlnr(self, r):
a3 = self.a3(r)
sigma = self.sigma(a3)
power = _spline(self.k, self.power)(1 / a3)
return -power / (2 * np.pi ** 2 * sigma ** 2 * a3 ** 3)
[docs] def dlnr_dlnm(self, r):
a3 = self.a3(r)
xi = r / a3
drda = _spline(a3, r).derivative()(a3)
return xi / 3 / drda