From da6f37689582692bfbafab5875ce9495a426b5c3 Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Fri, 24 Jun 2022 15:55:09 +0200 Subject: [PATCH 01/10] Adding function to compute projected correlation functions w_ab. - Compute the 3D correlation function xi_ab(r,z) integrating using FFTLog. - Calculate the window function and use it to obtain w_ab. --- pyccl/correlations.py | 154 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 154 insertions(+) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index acd1de15f..9f7533119 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -10,6 +10,7 @@ from . import constants as const from .core import check from .pk2d import parse_pk2d +from .pyutils import resample_array, _fftlog_transform import numpy as np import warnings @@ -330,3 +331,156 @@ def correlation_pi_sigma(cosmo, a, beta, pi, sig, if scalar: return xis[0] return xis + + +def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, + dndz: np.ndarray, dndz2=None, + p_of_k_a=None, + type='gg', + precision_fftlog: dict = None): + """Computes :math:`w_{ab}`. + .. math:: + \\w_{ab}(r_p)=\\int dz\\,(W)(z) + \\int\\frac{dk k}{2\\pi^2}\\,P(k,a)\\,J_n(k r_p) + + Args: + cosmo (:class:`~pyccl.core.Cosmology`): A Cosmology object. + r_p (float or array-like): Projected radial separation where the + correlation function will be evaluated. + z (array-like): Redshift values where the redshift distribution + has been evaluated. + dndz (array-like): Redshift distribution to be used when computing + the window function (see e.g. Eq. (3.11) in (1811.09598)). + dndz2 (array-like or None): Redshift distribution corresponding to + sample `b` of tracers to be correlated. If tracers `a` and `b` + are the same, this should be left as None. + Default value: None. + p_of_k_a (:class:`~pyccl.pk2d.Pk2D`): 3D Power spectrum to + integrate. + type (string): Type of `ab` correlation. This changes the Bessel + function in the integral above. Choices: 'gg' (`J_0`), + 'g+' (`J_2`), '++' (`J_0+J_4`). + precision_fftlog (dict): Dictionary containing the precision + parameters used by FFTLog to compute Hankel transforms. The + dictionary should contain the keys: `padding_lo_fftlog`, + `padding_lo_fftlog`, `padding_lo_extra`, `padding_hi_fftlog`, + `padding_hi_extra`, `n_per_decade`, `extrapol`, + `plaw_fourier`. Default values are 0.01, 0.1, 10., 10., + 100, 'linx_liny' and -1.5. For more information look at + `pyccl.halos.profiles.HaloProfile.update_precision_fftlog`. + + Returns: + array-like: Value(s) of the correlation function at the input + projected radial separation(s). + """ + if dndz2 is None: + dndz2 = dndz + a = 1/(1+z) + wz = dndz*dndz2*(1+z)**2/( + cosmo.comoving_radial_distance(a)**2. * + _compute_dchi_da_num(a, cosmo)) + wz /= np.trapz(wz, z) + if type == 'gg': + xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0) + for a_ in a]) + elif type == 'g+': + xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=2) + for a_ in a]) + elif type == '++': + xi = np.array([(_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0) + + _fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=4)) + for a_ in a]) + else: + raise ValueError + if np.ndim(xi) == 2: + wz = wz.reshape((len(z), 1)) + w = np.trapz(wz*xi, z, axis=0) + return w + + +def _fftlog_wrap(r_p: np.ndarray, a: float, cosmo, p_of_k_a=None, + precision_fftlog: dict = None, + n=0, large_padding=True): + '''Computes the integral + .. math:: + \\xi(r_p,a)=\\int\\frac{dk k}{2\\pi^2}\\,P(k,a)\\,J_n(k r_p) + Adapted from pyccl.halos.profiles.HaloProfile._fftlog_wrap. + ''' + # Note: The _fftlog_transform solves for + # f(r) = \int 4 pi k^2 f(k) j_l(kr) dk + # To compute the integral shown above we use _fftlog_transform + # and set fk = rp^1/2 P(k,a:float)/(2^5/2 pi^5/2 k^1/2). + if precision_fftlog is None: + precision_fftlog = {'padding_lo_fftlog': 0.01, + 'padding_lo_extra': 0.1, + 'padding_hi_fftlog': 10., + 'padding_hi_extra': 10., + 'n_per_decade': 100, + 'extrapol': 'linx_liny', + 'plaw_fourier': -1.5} + else: + for key in ['padding_lo_fftlog', 'padding_lo_extra', + 'padding_hi_fftlog', 'padding_hi_extra', 'n_per_decade', + 'extrapol', 'plaw_fourier']: + if key not in precision_fftlog.keys(): + raise ValueError('Error.') + l = n - 1 / 2 + + r_use = np.atleast_1d(r_p) + lr_use = np.log(r_use) + + # k-range to be used with FFTLog and its sampling. + if large_padding: + k_min = precision_fftlog['padding_lo_fftlog'] * np.amin(r_use) + k_max = precision_fftlog['padding_hi_fftlog'] * np.amax(r_use) + else: + k_min = precision_fftlog['padding_lo_extra'] * np.amin(r_use) + k_max = precision_fftlog['padding_hi_extra'] * np.amax(r_use) + n_k = (int(np.log10(k_max / k_min)) * + precision_fftlog['n_per_decade']) + k_arr = np.geomspace(k_min, k_max, n_k) # Array to be used by FFTLog. + + # What needs to go in _fftlog_transform to get out \xi_gg. + fk = p_of_k_a.eval(k_arr, a, cosmo)/((2*np.pi)**(5./2)*k_arr**(1/2)) + r_arr, xi_fourier = _fftlog_transform(k_arr, fk, 3, l, + precision_fftlog['plaw_fourier']) + lr_arr = np.log(r_arr) + + # Resample into input k values. + xi_out = resample_array(lr_arr, xi_fourier, lr_use, + precision_fftlog['extrapol'], + precision_fftlog['extrapol'], + 0, 0) + # Multiply to get the correct output. + xi_out *= (2*np.pi)**3*np.sqrt(r_use) + + if np.ndim(r_p) == 0: + xi_out = np.squeeze(xi_out, axis=-1) + return xi_out + + +def _compute_dchi_da_num(a, cosmo, a_step=0.01): + a_use = np.atleast_1d(a) + # If the a-array has values higher than 1, throw an error. + if len(a_use[a_use > 1]) > 0: + raise ValueError('Cannot accept scale factor larger than 1.') + # If the a-array has values that go beyond 1 when the step size is added + # reduce the step size until that doesn't happen anymore and + # throw a warning. + if len(a_use[(a_use + a_step > 1) & (a_use < 1)] > 0): + while len(a_use[(a_use + a_step > 1) & (a_use < 1)] > 0): + a_step /= 2. + warnings.warn('Step size causes scale factor > 1. Using step=%f' + % a_step) + dxda = np.empty(a_use.shape) + # Compute for all values of a that do not go above 1. + dxda[a_use < 1] = np.array([ + (cosmo.comoving_radial_distance(a + a_step) - + cosmo.comoving_radial_distance(a - a_step)) / (2 * a_step) + for a in a_use[a_use < 1]]) + # For the values that go above 1, + # use a constant value that is very close to 1. + dxda[a_use == 1] = (cosmo.comoving_radial_distance(0.9999 + 0.00001) - + cosmo.comoving_radial_distance( + 0.9999 - 0.00001)) / (2 * 0.00001) + return dxda From 0700e2a0a126453ffb2b13075e767fb953943284 Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Fri, 24 Jun 2022 16:45:11 +0200 Subject: [PATCH 02/10] Set window function to zero when input redshift is 0. - Not the best way to do this. --- pyccl/correlations.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index 9f7533119..5e9cf1dd5 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -379,6 +379,7 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, wz = dndz*dndz2*(1+z)**2/( cosmo.comoving_radial_distance(a)**2. * _compute_dchi_da_num(a, cosmo)) + wz[np.isnan(wz)] = 0 wz /= np.trapz(wz, z) if type == 'gg': xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0) From 691103759bdf3ee5123ff05d121ef89a304e512d Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Thu, 8 Sep 2022 15:31:15 +0200 Subject: [PATCH 03/10] Added functionality for when input redshift is a float. --- pyccl/correlations.py | 52 +++++++++++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index 5e9cf1dd5..57cf34997 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -334,9 +334,9 @@ def correlation_pi_sigma(cosmo, a, beta, pi, sig, def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, - dndz: np.ndarray, dndz2=None, + dndz: np.ndarray = None, dndz2: np.ndarray = None, p_of_k_a=None, - type='gg', + type: str = 'gg', precision_fftlog: dict = None): """Computes :math:`w_{ab}`. .. math:: @@ -348,13 +348,13 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, r_p (float or array-like): Projected radial separation where the correlation function will be evaluated. z (array-like): Redshift values where the redshift distribution - has been evaluated. + has been evaluated. If float, the window function is set to + a delta function. dndz (array-like): Redshift distribution to be used when computing the window function (see e.g. Eq. (3.11) in (1811.09598)). dndz2 (array-like or None): Redshift distribution corresponding to - sample `b` of tracers to be correlated. If tracers `a` and `b` - are the same, this should be left as None. - Default value: None. + sample `b` of tracers to be correlated. If None, tracers `a` + and `b` are the same. Default value: None. p_of_k_a (:class:`~pyccl.pk2d.Pk2D`): 3D Power spectrum to integrate. type (string): Type of `ab` correlation. This changes the Bessel @@ -373,29 +373,43 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, array-like: Value(s) of the correlation function at the input projected radial separation(s). """ + # TODO: add functionality for p_of_k_a=None. + # TODO: add some unit tests for Value Errors. if dndz2 is None: dndz2 = dndz + if dndz is None and np.iterable(z): + raise ValueError('Cannot have iterable redshift array while ' + 'dndz is None.') a = 1/(1+z) - wz = dndz*dndz2*(1+z)**2/( - cosmo.comoving_radial_distance(a)**2. * - _compute_dchi_da_num(a, cosmo)) - wz[np.isnan(wz)] = 0 - wz /= np.trapz(wz, z) + a = np.atleast_1d(a) if type == 'gg': - xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0) + xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0, + precision_fftlog=precision_fftlog) for a_ in a]) elif type == 'g+': - xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=2) + xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=2, + precision_fftlog=precision_fftlog) for a_ in a]) elif type == '++': - xi = np.array([(_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0) + - _fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=4)) + xi = np.array([(_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0, + precision_fftlog=precision_fftlog) + + _fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=4, + precision_fftlog=precision_fftlog)) for a_ in a]) else: - raise ValueError - if np.ndim(xi) == 2: - wz = wz.reshape((len(z), 1)) - w = np.trapz(wz*xi, z, axis=0) + raise ValueError('Correlation type not recognised. Accepted' + 'values: "gg", "g+", "++".') + if np.iterable(z): + wz = dndz*dndz2*(1+z)**2/( + cosmo.comoving_radial_distance(a)**2. * + _compute_dchi_da_num(a, cosmo)) + wz[np.isnan(wz)] = 0 # For values at z=0 + wz /= np.trapz(wz, z) + if np.ndim(xi) == 2: + wz = wz.reshape((len(z), 1)) + w = np.trapz(wz*xi, z, axis=0) + else: + w = xi[0] return w From b4d54be170dddf8a2baa1bc98e0b56d47ec8fa9c Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Mon, 27 Feb 2023 14:45:04 +0100 Subject: [PATCH 04/10] Analytical derivative for dchi/dz. Made FFTLog parameters easier for the user. --- pyccl/correlations.py | 57 +++++++++---------------------------------- 1 file changed, 12 insertions(+), 45 deletions(-) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index fd5f9f77f..a9379e4d5 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -11,6 +11,7 @@ from .core import check from .pk2d import parse_pk2d from .pyutils import resample_array, _fftlog_transform +from .parameters import physical_constants import numpy as np import warnings @@ -372,10 +373,10 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, 'g+' (`J_2`), '++' (`J_0+J_4`). precision_fftlog (dict): Dictionary containing the precision parameters used by FFTLog to compute Hankel transforms. The - dictionary should contain the keys: `padding_lo_fftlog`, - `padding_lo_fftlog`, `padding_lo_extra`, `padding_hi_fftlog`, - `padding_hi_extra`, `n_per_decade`, `extrapol`, - `plaw_fourier`. Default values are 0.01, 0.1, 10., 10., + dictionary should contain the keys: + `padding_lo_fftlog`, `padding_hi_fftlog`, + `n_per_decade`, `extrapol`, `plaw_fourier`. + Default values are 0.01, 0.1, 10., 10., 100, 'linx_liny' and -1.5. For more information look at `pyccl.halos.profiles.HaloProfile.update_precision_fftlog`. @@ -410,10 +411,10 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, raise ValueError('Correlation type not recognised. Accepted' 'values: "gg", "g+", "++".') if np.iterable(z): - wz = dndz*dndz2*(1+z)**2/( - cosmo.comoving_radial_distance(a)**2. * - _compute_dchi_da_num(a, cosmo)) - wz[np.isnan(wz)] = 0 # For values at z=0 + speed_of_light_kmps = physical_constants.CLIGHT * 1e-3 + dchi_dz = speed_of_light_kmps/(cosmo.h_over_h0(a)*cosmo['H0']) + wz = dndz * dndz2 / ( + cosmo.comoving_radial_distance(a) ** 2. * dchi_dz) wz /= np.trapz(wz, z) if np.ndim(xi) == 2: wz = wz.reshape((len(z), 1)) @@ -437,15 +438,12 @@ def _fftlog_wrap(r_p: np.ndarray, a: float, cosmo, p_of_k_a=None, # and set fk = rp^1/2 P(k,a:float)/(2^5/2 pi^5/2 k^1/2). if precision_fftlog is None: precision_fftlog = {'padding_lo_fftlog': 0.01, - 'padding_lo_extra': 0.1, 'padding_hi_fftlog': 10., - 'padding_hi_extra': 10., 'n_per_decade': 100, 'extrapol': 'linx_liny', 'plaw_fourier': -1.5} else: - for key in ['padding_lo_fftlog', 'padding_lo_extra', - 'padding_hi_fftlog', 'padding_hi_extra', 'n_per_decade', + for key in ['padding_lo_fftlog', 'padding_hi_fftlog', 'n_per_decade', 'extrapol', 'plaw_fourier']: if key not in precision_fftlog.keys(): raise ValueError('Error.') @@ -455,12 +453,8 @@ def _fftlog_wrap(r_p: np.ndarray, a: float, cosmo, p_of_k_a=None, lr_use = np.log(r_use) # k-range to be used with FFTLog and its sampling. - if large_padding: - k_min = precision_fftlog['padding_lo_fftlog'] * np.amin(r_use) - k_max = precision_fftlog['padding_hi_fftlog'] * np.amax(r_use) - else: - k_min = precision_fftlog['padding_lo_extra'] * np.amin(r_use) - k_max = precision_fftlog['padding_hi_extra'] * np.amax(r_use) + k_min = precision_fftlog['padding_lo_fftlog'] * np.amin(r_use) + k_max = precision_fftlog['padding_hi_fftlog'] * np.amax(r_use) n_k = (int(np.log10(k_max / k_min)) * precision_fftlog['n_per_decade']) k_arr = np.geomspace(k_min, k_max, n_k) # Array to be used by FFTLog. @@ -482,30 +476,3 @@ def _fftlog_wrap(r_p: np.ndarray, a: float, cosmo, p_of_k_a=None, if np.ndim(r_p) == 0: xi_out = np.squeeze(xi_out, axis=-1) return xi_out - - -def _compute_dchi_da_num(a, cosmo, a_step=0.01): - a_use = np.atleast_1d(a) - # If the a-array has values higher than 1, throw an error. - if len(a_use[a_use > 1]) > 0: - raise ValueError('Cannot accept scale factor larger than 1.') - # If the a-array has values that go beyond 1 when the step size is added - # reduce the step size until that doesn't happen anymore and - # throw a warning. - if len(a_use[(a_use + a_step > 1) & (a_use < 1)] > 0): - while len(a_use[(a_use + a_step > 1) & (a_use < 1)] > 0): - a_step /= 2. - warnings.warn('Step size causes scale factor > 1. Using step=%f' - % a_step) - dxda = np.empty(a_use.shape) - # Compute for all values of a that do not go above 1. - dxda[a_use < 1] = np.array([ - (cosmo.comoving_radial_distance(a + a_step) - - cosmo.comoving_radial_distance(a - a_step)) / (2 * a_step) - for a in a_use[a_use < 1]]) - # For the values that go above 1, - # use a constant value that is very close to 1. - dxda[a_use == 1] = (cosmo.comoving_radial_distance(0.9999 + 0.00001) - - cosmo.comoving_radial_distance( - 0.9999 - 0.00001)) / (2 * 0.00001) - return dxda From d9c57f16807d355f1981cb9a665440448853f3ea Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Tue, 13 Jun 2023 17:18:57 +0200 Subject: [PATCH 05/10] Updated correlation_ab to v3. --- pyccl/correlations.py | 53 +++++++++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 22 deletions(-) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index c69f9fb5d..6ad42af0c 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -1,12 +1,12 @@ __all__ = ("CorrelationMethods", "CorrelationTypes", "correlation", "correlation_3d", "correlation_multipole", "correlation_3dRsd", - "correlation_3dRsd_avgmu", "correlation_pi_sigma",) + "correlation_3dRsd_avgmu", "correlation_pi_sigma", + "correlation_ab") from enum import Enum import numpy as np from .pyutils import resample_array, _fftlog_transform -from .parameters import physical_constants -from . import DEFAULT_POWER_SPECTRUM, check, lib +from . import DEFAULT_POWER_SPECTRUM, check, lib, physical_constants class CorrelationMethods(Enum): @@ -354,12 +354,13 @@ def correlation_pi_sigma(cosmo, *, pi, sigma, a, beta, return xis -def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, +def correlation_ab(cosmo, *, r_p: np.ndarray, z: np.ndarray, dndz: np.ndarray = None, dndz2: np.ndarray = None, p_of_k_a=None, type: str = 'gg', precision_fftlog: dict = None): - """Computes :math:`w_{ab}`. + """ + Computes :math:`w_{ab}`. .. math:: \\w_{ab}(r_p)=\\int dz\\,(W)(z) \\int\\frac{dk k}{2\\pi^2}\\,P(k,a)\\,J_n(k r_p) @@ -368,16 +369,18 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, cosmo (:class:`~pyccl.core.Cosmology`): A Cosmology object. r_p (float or array-like): Projected radial separation where the correlation function will be evaluated. - z (array-like): Redshift values where the redshift distribution - has been evaluated. If float, the window function is set to - a delta function. + z (float or array-like): Redshift values where the redshift + distribution has been evaluated. If float, the window function + is set to a delta function. dndz (array-like): Redshift distribution to be used when computing the window function (see e.g. Eq. (3.11) in (1811.09598)). dndz2 (array-like or None): Redshift distribution corresponding to sample `b` of tracers to be correlated. If None, tracers `a` and `b` are the same. Default value: None. - p_of_k_a (:class:`~pyccl.pk2d.Pk2D`): 3D Power spectrum to - integrate. + p_of_k_a (:class:`~pyccl.pk2d.Pk2D`, :obj:`str` or :obj:`None`): 3D Power spectrum + to integrate. If a string, it must correspond to one of the + non-linear power spectra stored in `cosmo` (e.g. + `'delta_matter:delta_matter'`). type (string): Type of `ab` correlation. This changes the Bessel function in the integral above. Choices: 'gg' (`J_0`), 'g+' (`J_2`), '++' (`J_0+J_4`). @@ -393,9 +396,13 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, Returns: array-like: Value(s) of the correlation function at the input projected radial separation(s). - """ - # TODO: add functionality for p_of_k_a=None. + """ # noqa + cosmo.compute_nonlin_power() + cosmo_in = cosmo + psp = cosmo_in.parse_pk(p_of_k_a) # TODO: add some unit tests for Value Errors. + # TODO: add the benchmark test. + if dndz2 is None: dndz2 = dndz if dndz is None and np.iterable(z): @@ -404,17 +411,17 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, a = 1/(1+z) a = np.atleast_1d(a) if type == 'gg': - xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0, + xi = np.array([_fftlog_wrap(r_p, a_, psp, n=0, precision_fftlog=precision_fftlog) for a_ in a]) elif type == 'g+': - xi = np.array([_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=2, + xi = np.array([_fftlog_wrap(r_p, a_, psp, n=2, precision_fftlog=precision_fftlog) for a_ in a]) elif type == '++': - xi = np.array([(_fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=0, + xi = np.array([(_fftlog_wrap(r_p, a_, psp, n=0, precision_fftlog=precision_fftlog) + - _fftlog_wrap(r_p, a_, cosmo, p_of_k_a, n=4, + _fftlog_wrap(r_p, a_, psp, n=4, precision_fftlog=precision_fftlog)) for a_ in a]) else: @@ -422,9 +429,11 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, 'values: "gg", "g+", "++".') if np.iterable(z): speed_of_light_kmps = physical_constants.CLIGHT * 1e-3 - dchi_dz = speed_of_light_kmps/(cosmo.h_over_h0(a)*cosmo['H0']) - wz = dndz * dndz2 / ( - cosmo.comoving_radial_distance(a) ** 2. * dchi_dz) + dchi_dz = speed_of_light_kmps/(cosmo_in.h_over_h0(a)*cosmo_in['H0']) + denominator = cosmo_in.comoving_radial_distance(a) ** 2. * dchi_dz + wz = np.divide(dndz * dndz2, + denominator, + where=denominator > 0) wz /= np.trapz(wz, z) if np.ndim(xi) == 2: wz = wz.reshape((len(z), 1)) @@ -434,9 +443,9 @@ def correlation_ab(cosmo, r_p: np.ndarray, z: np.ndarray, return w -def _fftlog_wrap(r_p: np.ndarray, a: float, cosmo, p_of_k_a=None, +def _fftlog_wrap(r_p: np.ndarray, a: float, p_of_k_a, precision_fftlog: dict = None, - n=0, large_padding=True): + n=0): '''Computes the integral .. math:: \\xi(r_p,a)=\\int\\frac{dk k}{2\\pi^2}\\,P(k,a)\\,J_n(k r_p) @@ -470,7 +479,7 @@ def _fftlog_wrap(r_p: np.ndarray, a: float, cosmo, p_of_k_a=None, k_arr = np.geomspace(k_min, k_max, n_k) # Array to be used by FFTLog. # What needs to go in _fftlog_transform to get out \xi_gg. - fk = p_of_k_a.eval(k_arr, a, cosmo)/((2*np.pi)**(5./2)*k_arr**(1/2)) + fk = p_of_k_a(k_arr, a)/((2*np.pi)**(5./2)*k_arr**(1/2)) r_arr, xi_fourier = _fftlog_transform(k_arr, fk, 3, l, precision_fftlog['plaw_fourier']) lr_arr = np.log(r_arr) From 6bfb92f4349be87ac0f2010738b4f35fb013f365 Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Wed, 18 Oct 2023 11:26:25 +0200 Subject: [PATCH 06/10] Added unit test. --- pyccl/correlations.py | 2 -- pyccl/tests/test_correlations.py | 9 +++++++++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index cb8f1c31d..21eb5271c 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -411,8 +411,6 @@ def correlation_ab(cosmo, *, r_p: np.ndarray, z: np.ndarray, cosmo.compute_nonlin_power() cosmo_in = cosmo psp = cosmo_in.parse_pk(p_of_k_a) - # TODO: add some unit tests for Value Errors. - # TODO: add the benchmark test. if dndz2 is None: dndz2 = dndz diff --git a/pyccl/tests/test_correlations.py b/pyccl/tests/test_correlations.py index c421ec5b7..01d441881 100644 --- a/pyccl/tests/test_correlations.py +++ b/pyccl/tests/test_correlations.py @@ -133,4 +133,13 @@ def test_correlation_zero_ends(): ccl.correlation(COSMO, ell=ell, C_ell=C_ell, theta=theta) +def test_correlation_ab(): + r_p = np.geomspace(0.1, 100, 128) + z = np.linspace(0, 0.5, 128) + dndz = z ** 2 * np.exp(-(z / 0.11) ** 0.68) + with pytest.raises(ccl.CCLError): + COSMO.correlation_ab(r_p=r_p, z=z, dndz=None) + COSMO.correlation_ab(r_p=r_p, z=z, dndz=dndz, type='somethingelse') + + ccl.gsl_params.reload() # reset to the default parameters From ff94fb116e5907402d0fd2f4618d8ffcf8ed6775 Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Wed, 18 Oct 2023 11:48:37 +0200 Subject: [PATCH 07/10] Small fixes. --- pyccl/correlations.py | 4 ++-- pyccl/tests/test_correlations.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index 21eb5271c..9647355d1 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -402,7 +402,7 @@ def correlation_ab(cosmo, *, r_p: np.ndarray, z: np.ndarray, `n_per_decade`, `extrapol`, `plaw_fourier`. Default values are 0.01, 0.1, 10., 10., 100, 'linx_liny' and -1.5. For more information look at - `pyccl.halos.profiles.HaloProfile.update_precision_fftlog`. + `pyccl.halos.profiles.profile_base.update_precision_fftlog`. Returns: array-like: Value(s) of the correlation function at the input @@ -458,7 +458,7 @@ def _fftlog_wrap(r_p: np.ndarray, a: float, p_of_k_a, '''Computes the integral .. math:: \\xi(r_p,a)=\\int\\frac{dk k}{2\\pi^2}\\,P(k,a)\\,J_n(k r_p) - Adapted from pyccl.halos.profiles.HaloProfile._fftlog_wrap. + Adapted from pyccl.halos.profiles.profile_base._fftlog_wrap. ''' # Note: The _fftlog_transform solves for # f(r) = \int 4 pi k^2 f(k) j_l(kr) dk diff --git a/pyccl/tests/test_correlations.py b/pyccl/tests/test_correlations.py index 01d441881..26db69b1d 100644 --- a/pyccl/tests/test_correlations.py +++ b/pyccl/tests/test_correlations.py @@ -137,7 +137,7 @@ def test_correlation_ab(): r_p = np.geomspace(0.1, 100, 128) z = np.linspace(0, 0.5, 128) dndz = z ** 2 * np.exp(-(z / 0.11) ** 0.68) - with pytest.raises(ccl.CCLError): + with pytest.raises(ValueError): COSMO.correlation_ab(r_p=r_p, z=z, dndz=None) COSMO.correlation_ab(r_p=r_p, z=z, dndz=dndz, type='somethingelse') From 296a99e4a909aa4bea7b0a3727bf298a21091617 Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Wed, 13 Mar 2024 17:24:28 +0100 Subject: [PATCH 08/10] Added benchmark test and data. --- benchmarks/data/wgplus.out | 16 ++++++++++++++++ benchmarks/test_correlation_ab.py | 21 +++++++++++++++++++++ pyccl/correlations.py | 2 +- 3 files changed, 38 insertions(+), 1 deletion(-) create mode 100644 benchmarks/data/wgplus.out create mode 100644 benchmarks/test_correlation_ab.py diff --git a/benchmarks/data/wgplus.out b/benchmarks/data/wgplus.out new file mode 100644 index 000000000..7a8ec0fcd --- /dev/null +++ b/benchmarks/data/wgplus.out @@ -0,0 +1,16 @@ +1.4925e+00 3.6363e-01 +2.0289e+00 3.0676e-01 +2.7580e+00 2.3818e-01 +3.7491e+00 1.8313e-01 +5.0963e+00 1.3000e-01 +6.9277e+00 9.4586e-02 +9.4173e+00 7.0554e-02 +1.2801e+01 5.3728e-02 +1.7402e+01 4.3597e-02 +2.3655e+01 3.5296e-02 +3.2156e+01 2.8372e-02 +4.3711e+01 2.2266e-02 +5.9419e+01 1.6761e-02 +8.0772e+01 1.1925e-02 +1.0980e+02 7.9172e-03 +1.4925e+02 4.5258e-03 diff --git a/benchmarks/test_correlation_ab.py b/benchmarks/test_correlation_ab.py new file mode 100644 index 000000000..309782826 --- /dev/null +++ b/benchmarks/test_correlation_ab.py @@ -0,0 +1,21 @@ +import numpy as np +import pyccl as ccl +import os + + +def test_correlation_ab(): + cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, A_s=2.1e-9, n_s=0.96) + z_eval = 0. + + dirdat = os.path.dirname(__file__) + '/data/' + benchmark_file = os.path.join(dirdat, 'wgplus.out') + rp_benchmark, wgplus_benchmark = np.loadtxt(benchmark_file, unpack=True, delimiter=' ') + + C1rhocrit=0.0134 # standard IA normalisation + Pk_GI = ccl.Pk2D.from_function(pkfunc=lambda k, a: C1rhocrit * cosmo['Omega_m'] + / cosmo.growth_factor(a) * + cosmo.nonlin_matter_power(k, a), + is_logp=False) + wgplus = cosmo.correlation_ab(r_p=rp_benchmark, z=z_eval, p_of_k_a=Pk_GI, type='g+') + + assert np.all(np.abs((wgplus-wgplus_benchmark)/wgplus_benchmark)<5e-2) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index 5b1774c0b..fc146dd6b 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -465,7 +465,7 @@ def _fftlog_wrap(r_p: np.ndarray, a: float, p_of_k_a, # To compute the integral shown above we use _fftlog_transform # and set fk = rp^1/2 P(k,a:float)/(2^5/2 pi^5/2 k^1/2). if precision_fftlog is None: - precision_fftlog = {'padding_lo_fftlog': 0.01, + precision_fftlog = {'padding_lo_fftlog': 0.001, 'padding_hi_fftlog': 10., 'n_per_decade': 100, 'extrapol': 'linx_liny', From 5b611e1729b41e4892f04ccab979a5fedabad243 Mon Sep 17 00:00:00 2001 From: Christos Georgiou Date: Wed, 13 Mar 2024 17:32:44 +0100 Subject: [PATCH 09/10] Flake8 fixing. --- benchmarks/test_correlation_ab.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/benchmarks/test_correlation_ab.py b/benchmarks/test_correlation_ab.py index 309782826..33f5c1da9 100644 --- a/benchmarks/test_correlation_ab.py +++ b/benchmarks/test_correlation_ab.py @@ -4,18 +4,22 @@ def test_correlation_ab(): - cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, A_s=2.1e-9, n_s=0.96) + cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67, + A_s=2.1e-9, n_s=0.96) z_eval = 0. dirdat = os.path.dirname(__file__) + '/data/' benchmark_file = os.path.join(dirdat, 'wgplus.out') - rp_benchmark, wgplus_benchmark = np.loadtxt(benchmark_file, unpack=True, delimiter=' ') + rp_benchmark, wgplus_benchmark = np.loadtxt(benchmark_file, + unpack=True, delimiter=' ') - C1rhocrit=0.0134 # standard IA normalisation - Pk_GI = ccl.Pk2D.from_function(pkfunc=lambda k, a: C1rhocrit * cosmo['Omega_m'] - / cosmo.growth_factor(a) * - cosmo.nonlin_matter_power(k, a), - is_logp=False) - wgplus = cosmo.correlation_ab(r_p=rp_benchmark, z=z_eval, p_of_k_a=Pk_GI, type='g+') + C1rhocrit = 0.0134 # standard IA normalisation + Pk_GI = ccl.Pk2D.from_function( + pkfunc=lambda k, a: C1rhocrit * cosmo['Omega_m'] / + cosmo.growth_factor(a) * + cosmo.nonlin_matter_power(k, a), + is_logp=False) + wgplus = cosmo.correlation_ab(r_p=rp_benchmark, z=z_eval, + p_of_k_a=Pk_GI, type='g+') - assert np.all(np.abs((wgplus-wgplus_benchmark)/wgplus_benchmark)<5e-2) + assert np.all(np.abs((wgplus-wgplus_benchmark)/wgplus_benchmark) < 5e-2) From a53338893d86e61ce1dc14a6d707711aeebb6052 Mon Sep 17 00:00:00 2001 From: chrgeorgiou Date: Thu, 24 Apr 2025 15:34:02 +0200 Subject: [PATCH 10/10] Correcting the w++ correlation. --- pyccl/correlations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyccl/correlations.py b/pyccl/correlations.py index fc146dd6b..8ec6dd886 100644 --- a/pyccl/correlations.py +++ b/pyccl/correlations.py @@ -432,7 +432,7 @@ def correlation_ab(cosmo, *, r_p: np.ndarray, z: np.ndarray, precision_fftlog=precision_fftlog) + _fftlog_wrap(r_p, a_, psp, n=4, precision_fftlog=precision_fftlog)) - for a_ in a]) + for a_ in a])/2 else: raise ValueError('Correlation type not recognised. Accepted' 'values: "gg", "g+", "++".')