5191066774

Committed 06 Jun 2023 04:42PM UTC coverage: 97.543% (+0.4%) from 97.136%

Build # 5191066774

Build Type

Pull #1087

github

Committed by

web-flow

Commit Message

Merge 987bfdbae into 17a0e5a2a

Pull Request Pull Request #1087: v3 preparation

Run Details

83 of 83 new or added lines in 36 files covered. (100.0%)

5122 of 5251 relevant lines covered (97.54%)

0.98 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

98.62

/pyccl/halos/profiles/profile_base.py

__all__ = ("HaloProfile", "HaloProfileMatter",
           "HaloProfilePressure", "HaloProfileCIB",)

import functools
from typing import Callable

import numpy as np

from ... import CCLAutoRepr, FFTLogParams, unlock_instance
from ... import physical_constants as const
from ...pyutils import resample_array, _fftlog_transform
from .. import MassDef


class HaloProfile(CCLAutoRepr):
    """ This class implements functionality associated to
    halo profiles. You should not use this class directly.
    Instead, use one of the subclasses implemented in CCL
    for specific halo profiles, or write your own subclass.
    ``HaloProfile`` classes contain methods to compute halo
    profiles in real (3D) and Fourier spaces, as well as other
    projected (2D) quantities.

    A minimal implementation of a ``HaloProfile`` subclass
    should contain a method ``_real`` that returns the
    real-space profile as a function of cosmology,
    comoving radius, mass and scale factor. The default
    functionality in the base ``HaloProfile`` class will then
    allow the automatic calculation of the Fourier-space
    and projected profiles, as well as the cumulative
    mass density, based on the real-space profile using
    FFTLog to carry out fast Hankel transforms. See the
    CCL note for details. Alternatively, you can implement
    a ``_fourier`` method for the Fourier-space profile, and
    all other quantities will be computed from it. It is
    also possible to implement specific versions of any
    of these quantities if one wants to avoid the FFTLog
    calculation.
    """

    def __init__(self, *, mass_def, concentration=None,
                 is_number_counts=False):
        # Verify that profile can be initialized.
        if not (hasattr(self, "_real") or hasattr(self, "_fourier")):
            name = type(self).__name__
            raise TypeError(f"Can't instantiate {name} with no "
                            "_real or _fourier implementation.")

        # Initialize FFTLog.
        self.precision_fftlog = FFTLogParams()

        self._is_number_counts = is_number_counts

        # Initialize mass_def and concentration.
        self.mass_def, *out = MassDef.from_specs(
            mass_def, concentration=concentration)
        if out:
            self.concentration = out[0]

    @property
    def is_number_counts(self):
        """If ``True``, this profile represents source-count overdensities,
        normalised by the mean number density within their survey window
        function. This must be propagated when estimated super-sample
        effects in the covariance matrix.
        """
        return self._is_number_counts

    @is_number_counts.setter
    @unlock_instance
    def is_number_counts(self, value):
        self._is_number_counts = value

    def get_normalization(self, cosmo=None, a=None, *, hmc=None):
        """Profiles may be normalized by an overall function of redshift
        (or scale factor). This function may be cosmology-dependent and
        often comes from integrating certain halo properties over mass.
        This method returns this normalizing factor. For example,
        to get the normalized profile in real space, one would call
        the :meth:`real` method, and then **divide** the result by the value
        returned by this method.

        Args:
            hmc (:class:`~pyccl.halos.halo_model.HMCalculator`): a halo
                model calculator object.
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            a (:obj:`float`): scale factor.

        Returns:
            float: normalization factor of this profile.
        """
        return 1.0

    @unlock_instance(mutate=True)
    @functools.wraps(FFTLogParams.update_parameters)
    def update_precision_fftlog(self, **kwargs):
        self.precision_fftlog.update_parameters(**kwargs)

    def _get_plaw_fourier(self, cosmo, a):
        """ This controls the value of `plaw_fourier` to be used
        as a function of cosmology and scale factor.

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            a (:obj:`float`): scale factor.

        Returns:
            float: power law index to be used with FFTLog.
        """
        return self.precision_fftlog['plaw_fourier']

    def _get_plaw_projected(self, cosmo, a):
        """ This controls the value of `plaw_projected` to be
        used as a function of cosmology and scale factor.

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            a (:obj:`float`): scale factor.

        Returns:
            float: power law index to be used with FFTLog.
        """
        return self.precision_fftlog['plaw_projected']

    _real: Callable       # implementation of the real profile

    _fourier: Callable    # implementation of the Fourier profile

    _projected: Callable  # implementation of the projected profile

    _cumul2d: Callable    # implementation of the cumulative surface density

    def real(self, cosmo, r, M, a):
        """
        real(cosmo, r, M, a)
        Returns the 3D real-space value of the profile as a
        function of cosmology, radius, halo mass and scale factor.

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            r (:obj:`float` or `array`): comoving radius in Mpc.
            M (:obj:`float` or `array`): halo mass in units of M_sun.
            a (:obj:`float`): scale factor.

        Returns:
            (:obj:`float` or `array`): halo profile. The shape of the
            output will be `(N_M, N_r)` where `N_r` and `N_m` are
            the sizes of `r` and `M` respectively. If `r` or `M`
            are scalars, the corresponding dimension will be
            squeezed out on output.
        """
        if getattr(self, "_real", None):
            return self._real(cosmo, r, M, a)
        return self._fftlog_wrap(cosmo, r, M, a, fourier_out=False)

    def fourier(self, cosmo, k, M, a):
        """
        fourier(cosmo, k, M, a)
        Returns the Fourier-space value of the profile as a
        function of cosmology, wavenumber, halo mass and
        scale factor.

        .. math::
           \\rho(k)=\\frac{1}{2\\pi^2} \\int dr\\, r^2\\,
           \\rho(r)\\, j_0(k r)

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            k (:obj:`float` or `array`): comoving wavenumber (in :math:`Mpc^{-1}`).
            M (:obj:`float` or `array`): halo mass.
            a (:obj:`float`): scale factor.

        Returns:
            (:obj:`float` or `array`): Fourier-space profile. The shape of the
            output will be ``(N_M, N_k)`` where ``N_k`` and ``N_m`` are
            the sizes of ``k`` and ``M`` respectively. If ``k`` or ``M``
            are scalars, the corresponding dimension will be
            squeezed out on output.
        """ # noqa
        if getattr(self, "_fourier", None):
            return self._fourier(cosmo, k, M, a)
        return self._fftlog_wrap(cosmo, k, M, a, fourier_out=True)

    def projected(self, cosmo, r_t, M, a):
        """
        projected(cosmo, r_t, M, a)
        Returns the 2D projected profile as a function of
        cosmology, radius, halo mass and scale factor.

        .. math::
           \\Sigma(R)= \\int dr_\\parallel\\,
           \\rho(\\sqrt{r_\\parallel^2 + R^2})

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            r_t (:obj:`float` or `array`): transverse comoving radius in Mpc.
            M (:obj:`float` or `array`): halo mass in units of M_sun.
            a (:obj:`float`): scale factor.

        Returns:
            (:obj:`float` or `array`): projected profile. The shape of the
            output will be `(N_M, N_r)` where `N_r` and `N_m` are
            the sizes of `r` and `M` respectively. If `r` or `M`
            are scalars, the corresponding dimension will be
            squeezed out on output.
        """
        if getattr(self, "_projected", None):
            return self._projected(cosmo, r_t, M, a)
        return self._projected_fftlog_wrap(cosmo, r_t, M, a, is_cumul2d=False)

    def cumul2d(self, cosmo, r_t, M, a):
        """
        cumul2d(cosmo, r_t, M, a)
        Returns the 2D cumulative surface density as a
        function of cosmology, radius, halo mass and scale
        factor.

        .. math::
           \\Sigma(<R)= \\frac{2}{R^2} \\int dR'\\, R'\\,
           \\Sigma(R')

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            r_t (:obj:`float` or `array`): transverse comoving radius in Mpc.
            M (:obj:`float` or `array`): halo mass in units of M_sun.
            a (:obj:`float`): scale factor.

        Returns:
            (:obj:`float` or `array`): cumulative surface density. The shape of the
            output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m`` are
            the sizes of ``r`` and ``M`` respectively. If ``r`` or ``M``
            are scalars, the corresponding dimension will be
            squeezed out on output.
        """ # noqa
        if getattr(self, "_cumul2d", None):
            return self._cumul2d(cosmo, r_t, M, a)
        return self._projected_fftlog_wrap(cosmo, r_t, M, a, is_cumul2d=True)

    def convergence(self, cosmo, r, M, *, a_lens, a_source):
        """
        convergence(cosmo, r, M, *, a_lens, a_source)
        Returns the convergence as a function of cosmology,
        radius, halo mass and the scale factors of the source
        and the lens.

        .. math::
           \\kappa(R) = \\frac{\\Sigma(R)}{\\Sigma_{\\rm crit}},

        where :math:`\\Sigma(R)` is the 2D projected surface mass density
        (see :meth:`projected`), and :math:`\\Sigma_{\\rm crit}` is
        the critical surface density (see
        :meth:`~pyccl.background.sigma_critical`).

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
            r (:obj:`float` or `array`): comoving radius.
            M (:obj:`float` or `array`): halo mass.
            a_lens (:obj:`float`): scale factor of lens.
            a_source (:obj:`float` or `array`): scale factor of source.
                If an array, it must have the same shape as ``r``.

        Returns:
            (:obj:`float` or `array`): convergence :math:`\\kappa`. The
            shape of the output will be ``(N_M, N_r)`` where ``N_r`` and
            ``N_m`` are the sizes of ``r`` and ``M`` respectively. If
            ``r`` or ``M`` are scalars, the corresponding dimension will
            be squeezed out on output.
        """
        Sigma = self.projected(cosmo, r, M, a_lens)
        Sigma /= a_lens**2
        Sigma_crit = cosmo.sigma_critical(a_lens=a_lens, a_source=a_source)
        return Sigma / Sigma_crit

    def shear(self, cosmo, r, M, *, a_lens, a_source):
        """
        shear(cosmo, r, M, *, a_lens, a_source)
        Returns the shear (tangential) as a function of cosmology,
        radius, halo mass and the scale factors of the
        source and the lens.

        .. math::
           \\gamma(R) = \\frac{\\Delta\\Sigma(R)}{\\Sigma_{\\mathrm{crit}}} =
           \\frac{\\overline{\\Sigma}(< R) -
           \\Sigma(R)}{\\Sigma_{\\mathrm{crit}}},

        where :math:`\\overline{\\Sigma}(< R)` is cumulative surface density
        (see :meth:`cumul2d`), :math:`\\Sigma(R)` is the projected 2D density
        (see :meth:`projected`), and :math:`\\Sigma_{\\rm crit}` is the
        critical surface density (see
        :meth:`~pyccl.background.sigma_critical`).

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
            r (:obj:`float` or `array`): comoving radius in Mpc.
            M (:obj:`float` or `array`): halo mass in units of M_sun.
            a_lens (:obj:`float`): scale factor of lens.
            a_source (:obj:`float` or `array`): source's scale factor.
                If array, it must have the same shape as `r`.

        Returns:
            (:obj:`float` or `array`): shear :math:`\\gamma`. The shape of the
            output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m`` are
            the sizes of ``r`` and ``M`` respectively. If ``r`` or ``M``
            are scalars, the corresponding dimension will be
            squeezed out on output.
        """
        Sigma = self.projected(cosmo, r, M, a_lens)
        Sigma_bar = self.cumul2d(cosmo, r, M, a_lens)
        Sigma_crit = cosmo.sigma_critical(a_lens=a_lens, a_source=a_source)
        return (Sigma_bar - Sigma) / (Sigma_crit * a_lens**2)

    def reduced_shear(self, cosmo, r, M, *, a_lens, a_source):
        """
        reduced_shear(cosmo, r, M, *, a_lens, a_source)
        Returns the reduced shear as a function of cosmology,
        radius, halo mass and the scale factors of the
        source and the lens.

        .. math::
           g_t (R) = \\frac{\\gamma(R)}{(1 - \\kappa(R))},

        where :math:`\\gamma(R)` is the shear and :math:`\\kappa(R)` is the
        convergence.

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
            r (:obj:`float` or `array`): comoving radius in Mpc.
            M (:obj:`float` or `array`): halo mass in units of M_sun.
            a_lens (:obj:`float`): scale factor of lens.
            a_source (:obj:`float` or `array`): source's scale factor.
                If array, it must have the same shape as `r`.

        Returns:
            (:obj:`float` or `array`): reduced shear :math:`g_t`. The shape
            of the output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m``
            are the sizes of ``r`` and ``M`` respectively. If ``r`` or
            ``M`` are scalars, the corresponding dimension will be
            squeezed out on output.
        """
        convergence = self.convergence(cosmo, r, M, a_lens=a_lens,
                                       a_source=a_source)
        shear = self.shear(cosmo, r, M, a_lens=a_lens, a_source=a_source)
        return shear / (1.0 - convergence)

    def magnification(self, cosmo, r, M, *, a_lens, a_source):
        """
        magnification(cosmo, r, M, a_lens, a_source)
        Returns the magnification for input parameters.

        .. math::
           \\mu(R) = \\frac{1}{\\left[(1 - \\kappa(R))^2 -
           \\vert \\gamma(R) \\vert^2 \\right]]},

        where :math:`\\gamma(R)` is the shear and :math:`\\kappa(R)` is the
        convergence (see :meth:`shear` and :meth:`convergence`).

        Args:
            cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
            r (:obj:`float` or `array`): comoving radius in Mpc.
            M (:obj:`float` or `array`): halo mass in units of M_sun.
            a_lens (:obj:`float`): scale factor of lens.
            a_source (:obj:`float` or `array`): source's scale factor.
                If array, it must have the same shape as `r`.

        Returns:
            (:obj:`float` or `array`): magnification :math:`\\mu`. The shape
            of the output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m``
            are the sizes of ``r`` and ``M`` respectively. If ``r`` or ``M``
            are scalars, the corresponding dimension will be
            squeezed out on output.
        """
        convergence = self.convergence(cosmo, r, M, a_lens=a_lens,
                                       a_source=a_source)
        shear = self.shear(cosmo, r, M, a_lens=a_lens, a_source=a_source)

        return 1.0 / ((1.0 - convergence)**2 - np.abs(shear)**2)

    def _fftlog_wrap(self, cosmo, k, M, a,
                     fourier_out=False,
                     large_padding=True, ell=0):
        # This computes the 3D Hankel transform
        #  \rho(k) = 4\pi \int dr r^2 \rho(r) j_ell(k r)
        # if fourier_out == True, and
        #  \rho(r) = \frac{1}{2\pi^2} \int dk k^2 \rho(k) j_ell(k r)
        # otherwise.

        # Select which profile should be the input
        if fourier_out:
            p_func = self._real
        else:
            p_func = self._fourier
        k_use = np.atleast_1d(k)
        M_use = np.atleast_1d(M)
        lk_use = np.log(k_use)
        nM = len(M_use)

        # k/r ranges to be used with FFTLog and its sampling.
        if large_padding:
            k_min = self.precision_fftlog['padding_lo_fftlog'] * np.amin(k_use)
            k_max = self.precision_fftlog['padding_hi_fftlog'] * np.amax(k_use)
        else:
            k_min = self.precision_fftlog['padding_lo_extra'] * np.amin(k_use)
            k_max = self.precision_fftlog['padding_hi_extra'] * np.amax(k_use)
        n_k = (int(np.log10(k_max / k_min)) *
               self.precision_fftlog['n_per_decade'])
        r_arr = np.geomspace(k_min, k_max, n_k)

        p_k_out = np.zeros([nM, k_use.size])
        # Compute real profile values
        p_real_M = p_func(cosmo, r_arr, M_use, a)
        # Power-law index to pass to FFTLog.
        plaw_index = self._get_plaw_fourier(cosmo, a)

        # Compute Fourier profile through fftlog
        k_arr, p_fourier_M = _fftlog_transform(r_arr, p_real_M,
                                               3, ell, plaw_index)
        lk_arr = np.log(k_arr)

        for im, p_k_arr in enumerate(p_fourier_M):
            # Resample into input k values
            p_fourier = resample_array(lk_arr, p_k_arr, lk_use,
                                       self.precision_fftlog['extrapol'],
                                       self.precision_fftlog['extrapol'],
                                       0, 0)
            p_k_out[im, :] = p_fourier
        if fourier_out:
            p_k_out *= (2 * np.pi)**3

        if np.ndim(k) == 0:
            p_k_out = np.squeeze(p_k_out, axis=-1)
        if np.ndim(M) == 0:
            p_k_out = np.squeeze(p_k_out, axis=0)
        return p_k_out

    def _projected_fftlog_wrap(self, cosmo, r_t, M, a, is_cumul2d=False):
        # This computes Sigma(R) from the Fourier-space profile as:
        # Sigma(R) = \frac{1}{2\pi} \int dk k J_0(k R) \rho(k)
        r_t_use = np.atleast_1d(r_t)
        M_use = np.atleast_1d(M)
        lr_t_use = np.log(r_t_use)
        nM = len(M_use)

        # k/r range to be used with FFTLog and its sampling.
        r_t_min = self.precision_fftlog['padding_lo_fftlog'] * np.amin(r_t_use)
        r_t_max = self.precision_fftlog['padding_hi_fftlog'] * np.amax(r_t_use)
        n_r_t = (int(np.log10(r_t_max / r_t_min)) *
                 self.precision_fftlog['n_per_decade'])
        k_arr = np.geomspace(r_t_min, r_t_max, n_r_t)

        sig_r_t_out = np.zeros([nM, r_t_use.size])
        # Compute Fourier-space profile
        if getattr(self, "_fourier", None):
            # Compute from `_fourier` if available.
            p_fourier = self._fourier(cosmo, k_arr, M_use, a)
        else:
            # Compute with FFTLog otherwise.
            lpad = self.precision_fftlog['large_padding_2D']
            p_fourier = self._fftlog_wrap(cosmo, k_arr, M_use, a,
                                          fourier_out=True, large_padding=lpad)
        if is_cumul2d:
            # The cumulative profile involves a factor 1/(k R) in
            # the integrand.
            p_fourier *= 2 / k_arr[None, :]

        # Power-law index to pass to FFTLog.
        if is_cumul2d:
            i_bessel = 1
            plaw_index = self._get_plaw_projected(cosmo, a) - 1
        else:
            i_bessel = 0
            plaw_index = self._get_plaw_projected(cosmo, a)

        # Compute projected profile through fftlog
        r_t_arr, sig_r_t_M = _fftlog_transform(k_arr, p_fourier,
                                               2, i_bessel,
                                               plaw_index)
        lr_t_arr = np.log(r_t_arr)

        if is_cumul2d:
            sig_r_t_M /= r_t_arr[None, :]
        for im, sig_r_t_arr in enumerate(sig_r_t_M):
            # Resample into input r_t values
            sig_r_t = resample_array(lr_t_arr, sig_r_t_arr,
                                     lr_t_use,
                                     self.precision_fftlog['extrapol'],
                                     self.precision_fftlog['extrapol'],
                                     0, 0)
            sig_r_t_out[im, :] = sig_r_t

        if np.ndim(r_t) == 0:
            sig_r_t_out = np.squeeze(sig_r_t_out, axis=-1)
        if np.ndim(M) == 0:
            sig_r_t_out = np.squeeze(sig_r_t_out, axis=0)
        return sig_r_t_out


class HaloProfileMatter(HaloProfile):
    """Base for matter halo profiles."""

    def get_normalization(self, cosmo, a, *, hmc=None):
        """Returns the normalization of all matter overdensity
        profiles, which is simply the comoving matter density.

        Args:
            hmc (:class:`~pyccl.halos.halo_model.HMCalculator`): a halo
                model calculator object.
            cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
            a (:obj:`float`): scale factor.

        Returns:
            :obj:`float`: normalization factor of this profile.
        """
        return const.RHO_CRITICAL * cosmo["Omega_m"] * cosmo["h"]**2


class HaloProfilePressure(HaloProfile):
    """Base for pressure halo profiles."""


class HaloProfileCIB(HaloProfile):
    """Base for CIB halo profiles."""

1	__all__ = ("HaloProfile", "HaloProfileMatter",	1✔
2	"HaloProfilePressure", "HaloProfileCIB",)
3
4	import functools	1✔
5	from typing import Callable	1✔
6
7	import numpy as np	1✔
8
9	from ... import CCLAutoRepr, FFTLogParams, unlock_instance	1✔
10	from ... import physical_constants as const	1✔
11	from ...pyutils import resample_array, _fftlog_transform	1✔
12	from .. import MassDef	1✔
13
14
15	class HaloProfile(CCLAutoRepr):	1✔
16	""" This class implements functionality associated to	1✔
17	halo profiles. You should not use this class directly.
18	Instead, use one of the subclasses implemented in CCL
19	for specific halo profiles, or write your own subclass.
20	``HaloProfile`` classes contain methods to compute halo
21	profiles in real (3D) and Fourier spaces, as well as other
22	projected (2D) quantities.
23
24	A minimal implementation of a ``HaloProfile`` subclass
25	should contain a method ``_real`` that returns the
26	real-space profile as a function of cosmology,
27	comoving radius, mass and scale factor. The default
28	functionality in the base ``HaloProfile`` class will then
29	allow the automatic calculation of the Fourier-space
30	and projected profiles, as well as the cumulative
31	mass density, based on the real-space profile using
32	FFTLog to carry out fast Hankel transforms. See the
33	CCL note for details. Alternatively, you can implement
34	a ``_fourier`` method for the Fourier-space profile, and
35	all other quantities will be computed from it. It is
36	also possible to implement specific versions of any
37	of these quantities if one wants to avoid the FFTLog
38	calculation.
39	"""
40
41	def __init__(self, *, mass_def, concentration=None,	1✔
42	is_number_counts=False):
43	# Verify that profile can be initialized.
44	if not (hasattr(self, "_real") or hasattr(self, "_fourier")):	1✔
45	name = type(self).__name__	×
46	raise TypeError(f"Can't instantiate {name} with no "	×
47	"_real or _fourier implementation.")
48
49	# Initialize FFTLog.
50	self.precision_fftlog = FFTLogParams()	1✔
51
52	self._is_number_counts = is_number_counts	1✔
53
54	# Initialize mass_def and concentration.
55	self.mass_def, *out = MassDef.from_specs(	1✔
56	mass_def, concentration=concentration)
57	if out:	1✔
58	self.concentration = out[0]	1✔
59
60	@property	1✔
61	def is_number_counts(self):	1✔
62	"""If ``True``, this profile represents source-count overdensities,
63	normalised by the mean number density within their survey window
64	function. This must be propagated when estimated super-sample
65	effects in the covariance matrix.
66	"""
67	return self._is_number_counts	1✔
68
69	@is_number_counts.setter	1✔
70	@unlock_instance	1✔
71	def is_number_counts(self, value):	1✔
72	self._is_number_counts = value	1✔
73
74	def get_normalization(self, cosmo=None, a=None, *, hmc=None):	1✔
75	"""Profiles may be normalized by an overall function of redshift
76	(or scale factor). This function may be cosmology-dependent and
77	often comes from integrating certain halo properties over mass.
78	This method returns this normalizing factor. For example,
79	to get the normalized profile in real space, one would call
80	the :meth:`real` method, and then divide the result by the value
81	returned by this method.
82
83	Args:
84	hmc (:class:`~pyccl.halos.halo_model.HMCalculator`): a halo
85	model calculator object.
86	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
87	a (:obj:`float`): scale factor.
88
89	Returns:
90	float: normalization factor of this profile.
91	"""
92	return 1.0	1✔
93
94	@unlock_instance(mutate=True)	1✔
95	@functools.wraps(FFTLogParams.update_parameters)	1✔
96	def update_precision_fftlog(self, **kwargs):	1✔
97	self.precision_fftlog.update_parameters(**kwargs)	1✔
98
99	def _get_plaw_fourier(self, cosmo, a):	1✔
100	""" This controls the value of `plaw_fourier` to be used
101	as a function of cosmology and scale factor.
102
103	Args:
104	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
105	a (:obj:`float`): scale factor.
106
107	Returns:
108	float: power law index to be used with FFTLog.
109	"""
110	return self.precision_fftlog['plaw_fourier']	1✔
111
112	def _get_plaw_projected(self, cosmo, a):	1✔
113	""" This controls the value of `plaw_projected` to be
114	used as a function of cosmology and scale factor.
115
116	Args:
117	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
118	a (:obj:`float`): scale factor.
119
120	Returns:
121	float: power law index to be used with FFTLog.
122	"""
123	return self.precision_fftlog['plaw_projected']	1✔
124
125	_real: Callable # implementation of the real profile	1✔
126
127	_fourier: Callable # implementation of the Fourier profile	1✔
128
129	_projected: Callable # implementation of the projected profile	1✔
130
131	_cumul2d: Callable # implementation of the cumulative surface density	1✔
132
133	def real(self, cosmo, r, M, a):	1✔
134	"""
135	real(cosmo, r, M, a)
136	Returns the 3D real-space value of the profile as a
137	function of cosmology, radius, halo mass and scale factor.
138
139	Args:
140	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
141	r (:obj:`float` or `array`): comoving radius in Mpc.
142	M (:obj:`float` or `array`): halo mass in units of M_sun.
143	a (:obj:`float`): scale factor.
144
145	Returns:
146	(:obj:`float` or `array`): halo profile. The shape of the
147	output will be `(N_M, N_r)` where `N_r` and `N_m` are
148	the sizes of `r` and `M` respectively. If `r` or `M`
149	are scalars, the corresponding dimension will be
150	squeezed out on output.
151	"""
152	if getattr(self, "_real", None):	1✔
153	return self._real(cosmo, r, M, a)	1✔
154	return self._fftlog_wrap(cosmo, r, M, a, fourier_out=False)	1✔
155
156	def fourier(self, cosmo, k, M, a):	1✔
157	"""
158	fourier(cosmo, k, M, a)
159	Returns the Fourier-space value of the profile as a
160	function of cosmology, wavenumber, halo mass and
161	scale factor.
162
163	.. math::
164	\\rho(k)=\\frac{1}{2\\pi^2} \\int dr\\, r^2\\,
165	\\rho(r)\\, j_0(k r)
166
167	Args:
168	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
169	k (:obj:`float` or `array`): comoving wavenumber (in :math:`Mpc^{-1}`).
170	M (:obj:`float` or `array`): halo mass.
171	a (:obj:`float`): scale factor.
172
173	Returns:
174	(:obj:`float` or `array`): Fourier-space profile. The shape of the
175	output will be ``(N_M, N_k)`` where ``N_k`` and ``N_m`` are
176	the sizes of ``k`` and ``M`` respectively. If ``k`` or ``M``
177	are scalars, the corresponding dimension will be
178	squeezed out on output.
179	""" # noqa
180	if getattr(self, "_fourier", None):	1✔
181	return self._fourier(cosmo, k, M, a)	1✔
182	return self._fftlog_wrap(cosmo, k, M, a, fourier_out=True)	1✔
183
184	def projected(self, cosmo, r_t, M, a):	1✔
185	"""
186	projected(cosmo, r_t, M, a)
187	Returns the 2D projected profile as a function of
188	cosmology, radius, halo mass and scale factor.
189
190	.. math::
191	\\Sigma(R)= \\int dr_\\parallel\\,
192	\\rho(\\sqrt{r_\\parallel^2 + R^2})
193
194	Args:
195	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
196	r_t (:obj:`float` or `array`): transverse comoving radius in Mpc.
197	M (:obj:`float` or `array`): halo mass in units of M_sun.
198	a (:obj:`float`): scale factor.
199
200	Returns:
201	(:obj:`float` or `array`): projected profile. The shape of the
202	output will be `(N_M, N_r)` where `N_r` and `N_m` are
203	the sizes of `r` and `M` respectively. If `r` or `M`
204	are scalars, the corresponding dimension will be
205	squeezed out on output.
206	"""
207	if getattr(self, "_projected", None):	1✔
208	return self._projected(cosmo, r_t, M, a)	1✔
209	return self._projected_fftlog_wrap(cosmo, r_t, M, a, is_cumul2d=False)	1✔
210
211	def cumul2d(self, cosmo, r_t, M, a):	1✔
212	"""
213	cumul2d(cosmo, r_t, M, a)
214	Returns the 2D cumulative surface density as a
215	function of cosmology, radius, halo mass and scale
216	factor.
217
218	.. math::
219	\\Sigma(<R)= \\frac{2}{R^2} \\int dR'\\, R'\\,
220	\\Sigma(R')
221
222	Args:
223	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
224	r_t (:obj:`float` or `array`): transverse comoving radius in Mpc.
225	M (:obj:`float` or `array`): halo mass in units of M_sun.
226	a (:obj:`float`): scale factor.
227
228	Returns:
229	(:obj:`float` or `array`): cumulative surface density. The shape of the
230	output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m`` are
231	the sizes of ``r`` and ``M`` respectively. If ``r`` or ``M``
232	are scalars, the corresponding dimension will be
233	squeezed out on output.
234	""" # noqa
235	if getattr(self, "_cumul2d", None):	1✔
236	return self._cumul2d(cosmo, r_t, M, a)	1✔
237	return self._projected_fftlog_wrap(cosmo, r_t, M, a, is_cumul2d=True)	1✔
238
239	def convergence(self, cosmo, r, M, *, a_lens, a_source):	1✔
240	"""
241	convergence(cosmo, r, M, *, a_lens, a_source)
242	Returns the convergence as a function of cosmology,
243	radius, halo mass and the scale factors of the source
244	and the lens.
245
246	.. math::
247	\\kappa(R) = \\frac{\\Sigma(R)}{\\Sigma_{\\rm crit}},
248
249	where :math:`\\Sigma(R)` is the 2D projected surface mass density
250	(see :meth:`projected`), and :math:`\\Sigma_{\\rm crit}` is
251	the critical surface density (see
252	:meth:`~pyccl.background.sigma_critical`).
253
254	Args:
255	cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
256	r (:obj:`float` or `array`): comoving radius.
257	M (:obj:`float` or `array`): halo mass.
258	a_lens (:obj:`float`): scale factor of lens.
259	a_source (:obj:`float` or `array`): scale factor of source.
260	If an array, it must have the same shape as ``r``.
261
262	Returns:
263	(:obj:`float` or `array`): convergence :math:`\\kappa`. The
264	shape of the output will be ``(N_M, N_r)`` where ``N_r`` and
265	``N_m`` are the sizes of ``r`` and ``M`` respectively. If
266	``r`` or ``M`` are scalars, the corresponding dimension will
267	be squeezed out on output.
268	"""
269	Sigma = self.projected(cosmo, r, M, a_lens)	1✔
270	Sigma /= a_lens**2	1✔
271	Sigma_crit = cosmo.sigma_critical(a_lens=a_lens, a_source=a_source)	1✔
272	return Sigma / Sigma_crit	1✔
273
274	def shear(self, cosmo, r, M, *, a_lens, a_source):	1✔
275	"""
276	shear(cosmo, r, M, *, a_lens, a_source)
277	Returns the shear (tangential) as a function of cosmology,
278	radius, halo mass and the scale factors of the
279	source and the lens.
280
281	.. math::
282	\\gamma(R) = \\frac{\\Delta\\Sigma(R)}{\\Sigma_{\\mathrm{crit}}} =
283	\\frac{\\overline{\\Sigma}(< R) -
284	\\Sigma(R)}{\\Sigma_{\\mathrm{crit}}},
285
286	where :math:`\\overline{\\Sigma}(< R)` is cumulative surface density
287	(see :meth:`cumul2d`), :math:`\\Sigma(R)` is the projected 2D density
288	(see :meth:`projected`), and :math:`\\Sigma_{\\rm crit}` is the
289	critical surface density (see
290	:meth:`~pyccl.background.sigma_critical`).
291
292	Args:
293	cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
294	r (:obj:`float` or `array`): comoving radius in Mpc.
295	M (:obj:`float` or `array`): halo mass in units of M_sun.
296	a_lens (:obj:`float`): scale factor of lens.
297	a_source (:obj:`float` or `array`): source's scale factor.
298	If array, it must have the same shape as `r`.
299
300	Returns:
301	(:obj:`float` or `array`): shear :math:`\\gamma`. The shape of the
302	output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m`` are
303	the sizes of ``r`` and ``M`` respectively. If ``r`` or ``M``
304	are scalars, the corresponding dimension will be
305	squeezed out on output.
306	"""
307	Sigma = self.projected(cosmo, r, M, a_lens)	1✔
308	Sigma_bar = self.cumul2d(cosmo, r, M, a_lens)	1✔
309	Sigma_crit = cosmo.sigma_critical(a_lens=a_lens, a_source=a_source)	1✔
310	return (Sigma_bar - Sigma) / (Sigma_crit * a_lens**2)	1✔
311
312	def reduced_shear(self, cosmo, r, M, *, a_lens, a_source):	1✔
313	"""
314	reduced_shear(cosmo, r, M, *, a_lens, a_source)
315	Returns the reduced shear as a function of cosmology,
316	radius, halo mass and the scale factors of the
317	source and the lens.
318
319	.. math::
320	g_t (R) = \\frac{\\gamma(R)}{(1 - \\kappa(R))},
321
322	where :math:`\\gamma(R)` is the shear and :math:`\\kappa(R)` is the
323	convergence.
324
325	Args:
326	cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
327	r (:obj:`float` or `array`): comoving radius in Mpc.
328	M (:obj:`float` or `array`): halo mass in units of M_sun.
329	a_lens (:obj:`float`): scale factor of lens.
330	a_source (:obj:`float` or `array`): source's scale factor.
331	If array, it must have the same shape as `r`.
332
333	Returns:
334	(:obj:`float` or `array`): reduced shear :math:`g_t`. The shape
335	of the output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m``
336	are the sizes of ``r`` and ``M`` respectively. If ``r`` or
337	``M`` are scalars, the corresponding dimension will be
338	squeezed out on output.
339	"""
340	convergence = self.convergence(cosmo, r, M, a_lens=a_lens,	1✔
341	a_source=a_source)
342	shear = self.shear(cosmo, r, M, a_lens=a_lens, a_source=a_source)	1✔
343	return shear / (1.0 - convergence)	1✔
344
345	def magnification(self, cosmo, r, M, *, a_lens, a_source):	1✔
346	"""
347	magnification(cosmo, r, M, a_lens, a_source)
348	Returns the magnification for input parameters.
349
350	.. math::
351	\\mu(R) = \\frac{1}{\\left[(1 - \\kappa(R))^2 -
352	\\vert \\gamma(R) \\vert^2 \\right]]},
353
354	where :math:`\\gamma(R)` is the shear and :math:`\\kappa(R)` is the
355	convergence (see :meth:`shear` and :meth:`convergence`).
356
357	Args:
358	cosmo (:class:`~pyccl.cosmology.Cosmology`): A Cosmology object.
359	r (:obj:`float` or `array`): comoving radius in Mpc.
360	M (:obj:`float` or `array`): halo mass in units of M_sun.
361	a_lens (:obj:`float`): scale factor of lens.
362	a_source (:obj:`float` or `array`): source's scale factor.
363	If array, it must have the same shape as `r`.
364
365	Returns:
366	(:obj:`float` or `array`): magnification :math:`\\mu`. The shape
367	of the output will be ``(N_M, N_r)`` where ``N_r`` and ``N_m``
368	are the sizes of ``r`` and ``M`` respectively. If ``r`` or ``M``
369	are scalars, the corresponding dimension will be
370	squeezed out on output.
371	"""
372	convergence = self.convergence(cosmo, r, M, a_lens=a_lens,	1✔
373	a_source=a_source)
374	shear = self.shear(cosmo, r, M, a_lens=a_lens, a_source=a_source)	1✔
375
376	return 1.0 / ((1.0 - convergence)2 - np.abs(shear)2)	1✔
377
378	def _fftlog_wrap(self, cosmo, k, M, a,	1✔
379	fourier_out=False,
380	large_padding=True, ell=0):
381	# This computes the 3D Hankel transform
382	# \rho(k) = 4\pi \int dr r^2 \rho(r) j_ell(k r)
383	# if fourier_out == True, and
384	# \rho(r) = \frac{1}{2\pi^2} \int dk k^2 \rho(k) j_ell(k r)
385	# otherwise.
386
387	# Select which profile should be the input
388	if fourier_out:	1✔
389	p_func = self._real	1✔
390	else:
391	p_func = self._fourier	1✔
392	k_use = np.atleast_1d(k)	1✔
393	M_use = np.atleast_1d(M)	1✔
394	lk_use = np.log(k_use)	1✔
395	nM = len(M_use)	1✔
396
397	# k/r ranges to be used with FFTLog and its sampling.
398	if large_padding:	1✔
399	k_min = self.precision_fftlog['padding_lo_fftlog'] * np.amin(k_use)	1✔
400	k_max = self.precision_fftlog['padding_hi_fftlog'] * np.amax(k_use)	1✔
401	else:
402	k_min = self.precision_fftlog['padding_lo_extra'] * np.amin(k_use)	1✔
403	k_max = self.precision_fftlog['padding_hi_extra'] * np.amax(k_use)	1✔
404	n_k = (int(np.log10(k_max / k_min)) *	1✔
405	self.precision_fftlog['n_per_decade'])
406	r_arr = np.geomspace(k_min, k_max, n_k)	1✔
407
408	p_k_out = np.zeros([nM, k_use.size])	1✔
409	# Compute real profile values
410	p_real_M = p_func(cosmo, r_arr, M_use, a)	1✔
411	# Power-law index to pass to FFTLog.
412	plaw_index = self._get_plaw_fourier(cosmo, a)	1✔
413
414	# Compute Fourier profile through fftlog
415	k_arr, p_fourier_M = _fftlog_transform(r_arr, p_real_M,	1✔
416	3, ell, plaw_index)
417	lk_arr = np.log(k_arr)	1✔
418
419	for im, p_k_arr in enumerate(p_fourier_M):	1✔
420	# Resample into input k values
421	p_fourier = resample_array(lk_arr, p_k_arr, lk_use,	1✔
422	self.precision_fftlog['extrapol'],
423	self.precision_fftlog['extrapol'],
424	0, 0)
425	p_k_out[im, :] = p_fourier	1✔
426	if fourier_out:	1✔
427	p_k_out = (2 np.pi)**3	1✔
428
429	if np.ndim(k) == 0:	1✔
430	p_k_out = np.squeeze(p_k_out, axis=-1)	1✔
431	if np.ndim(M) == 0:	1✔
432	p_k_out = np.squeeze(p_k_out, axis=0)	1✔
433	return p_k_out	1✔
434
435	def _projected_fftlog_wrap(self, cosmo, r_t, M, a, is_cumul2d=False):	1✔
436	# This computes Sigma(R) from the Fourier-space profile as:
437	# Sigma(R) = \frac{1}{2\pi} \int dk k J_0(k R) \rho(k)
438	r_t_use = np.atleast_1d(r_t)	1✔
439	M_use = np.atleast_1d(M)	1✔
440	lr_t_use = np.log(r_t_use)	1✔
441	nM = len(M_use)	1✔
442
443	# k/r range to be used with FFTLog and its sampling.
444	r_t_min = self.precision_fftlog['padding_lo_fftlog'] * np.amin(r_t_use)	1✔
445	r_t_max = self.precision_fftlog['padding_hi_fftlog'] * np.amax(r_t_use)	1✔
446	n_r_t = (int(np.log10(r_t_max / r_t_min)) *	1✔
447	self.precision_fftlog['n_per_decade'])
448	k_arr = np.geomspace(r_t_min, r_t_max, n_r_t)	1✔
449
450	sig_r_t_out = np.zeros([nM, r_t_use.size])	1✔
451	# Compute Fourier-space profile
452	if getattr(self, "_fourier", None):	1✔
453	# Compute from `_fourier` if available.
454	p_fourier = self._fourier(cosmo, k_arr, M_use, a)	1✔
455	else:
456	# Compute with FFTLog otherwise.
457	lpad = self.precision_fftlog['large_padding_2D']	1✔
458	p_fourier = self._fftlog_wrap(cosmo, k_arr, M_use, a,	1✔
459	fourier_out=True, large_padding=lpad)
460	if is_cumul2d:	1✔
461	# The cumulative profile involves a factor 1/(k R) in
462	# the integrand.
463	p_fourier *= 2 / k_arr[None, :]	1✔
464
465	# Power-law index to pass to FFTLog.
466	if is_cumul2d:	1✔
467	i_bessel = 1	1✔
468	plaw_index = self._get_plaw_projected(cosmo, a) - 1	1✔
469	else:
470	i_bessel = 0	1✔
471	plaw_index = self._get_plaw_projected(cosmo, a)	1✔
472
473	# Compute projected profile through fftlog
474	r_t_arr, sig_r_t_M = _fftlog_transform(k_arr, p_fourier,	1✔
475	2, i_bessel,
476	plaw_index)
477	lr_t_arr = np.log(r_t_arr)	1✔
478
479	if is_cumul2d:	1✔
480	sig_r_t_M /= r_t_arr[None, :]	1✔
481	for im, sig_r_t_arr in enumerate(sig_r_t_M):	1✔
482	# Resample into input r_t values
483	sig_r_t = resample_array(lr_t_arr, sig_r_t_arr,	1✔
484	lr_t_use,
485	self.precision_fftlog['extrapol'],
486	self.precision_fftlog['extrapol'],
487	0, 0)
488	sig_r_t_out[im, :] = sig_r_t	1✔
489
490	if np.ndim(r_t) == 0:	1✔
491	sig_r_t_out = np.squeeze(sig_r_t_out, axis=-1)	1✔
492	if np.ndim(M) == 0:	1✔
493	sig_r_t_out = np.squeeze(sig_r_t_out, axis=0)	1✔
494	return sig_r_t_out	1✔
495
496
497	class HaloProfileMatter(HaloProfile):	1✔
498	"""Base for matter halo profiles."""	1✔
499
500	def get_normalization(self, cosmo, a, *, hmc=None):	1✔
501	"""Returns the normalization of all matter overdensity
502	profiles, which is simply the comoving matter density.
503
504	Args:
505	hmc (:class:`~pyccl.halos.halo_model.HMCalculator`): a halo
506	model calculator object.
507	cosmo (:class:`~pyccl.cosmology.Cosmology`): a Cosmology object.
508	a (:obj:`float`): scale factor.
509
510	Returns:
511	:obj:`float`: normalization factor of this profile.
512	"""
513	return const.RHO_CRITICAL * cosmo["Omega_m"] * cosmo["h"]**2	1✔
514
515
516	class HaloProfilePressure(HaloProfile):	1✔
517	"""Base for pressure halo profiles."""	1✔
518
519
520	class HaloProfileCIB(HaloProfile):	1✔
521	"""Base for CIB halo profiles."""	1✔

LSSTDESC / CCL / 5191066774

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous