Source code for camb.mathutils

"""
This module contains some fast utility functions that are useful in the same contexts as camb. They are entirely
independent of the main camb code.

"""

from ctypes import c_int, c_double, c_bool, POINTER
from .baseconfig import camblib, numpy_1d, numpy_2d, numpy_3d
import numpy as np

_chi2 = camblib.__mathutils_MOD_getchisquared
_chi2.argtypes = [numpy_2d, numpy_1d, POINTER(c_int)]
_chi2.restype = c_double


[docs]def chi_squared(covinv, x): """ Utility function to efficiently calculate x^T covinv x :param covinv: symmetric inverse covariance matrix :param x: vector :return: covinv.dot(x).dot(x), but parallelized and using symmetry """ if len(x) != covinv.shape[0] or covinv.shape[0] != covinv.shape[1]: raise ValueError('Wrong shape in chi_squared') return _chi2(covinv, x, c_int(len(x)))
int_arg = POINTER(c_int) _3j = camblib.__mathutils_MOD_getthreejs _3j.argtypes = [numpy_1d, int_arg, int_arg, int_arg, int_arg]
[docs]def threej(l2, l3, m2, m3): """ Convenience wrapper around standard 3j function, returning array for all allowed l1 values :param l2: L_2 :param l3: L_3 :param m2: M_2 :param m3: M_3 :return: array of 3j from max(abs(l2-l3),abs(m2+m3)) .. l2+l3 """ l1min = max(np.abs(l2 - l3), np.abs(m2 + m3)) result = np.zeros(int(l3 + l2 - l1min + 1)) l2in, l3in, m2in, m3in = c_int(l2), c_int(l3), c_int(m2), c_int(m3) _3j(result, l2in, l3in, m2in, m3in) return result
# Utils_3j_integrate(W,lmax_w, n, dopol, M, lmax) _coupling_3j = camblib.__mathutils_MOD_integrate_3j _coupling_3j.argtypes = [numpy_2d, POINTER(c_int), POINTER(c_int), POINTER(c_bool), numpy_3d, POINTER(c_int)]
[docs]def threej_coupling(W, lmax, pol=False): r""" Calculate symmetric coupling matrix :math`\Xi` for given weights :math:`W_{\ell}`, where :math:`\langle\tilde{C}_\ell\rangle = \Xi_{\ell \ell'} (2\ell'+1) C_\ell`. The weights are related to the power spectrum of the mask P by :math:`W_\ell = (2 \ell + 1) P_\ell / 4 \pi`. See e.g. Eq D16 of `arxiv:0801.0554 <http://arxiv.org/abs/0801.0554>`_. If pol is False and W is an array of weights, produces array of temperature couplings, otherwise for pol is True produces set of TT, TE, EE, EB couplings (and weights must have one spectrum - for same masks - or three). Use :func:`scalar_coupling_matrix` or :func:`pcl_coupling_matrix` to get the coupling matrix directly from the mask power spectrum. :param W: 1d array of Weights for each L, or list of arrays of weights (zero based) :param lmax: lmax for the output matrix (assumed symmetric, though not in principle) :param pol: if pol, produce TT, TE, EE, EB couplings for three input mask weights (or one if assuming same mask) :return: symmetric coupling matrix or array of matrices """ if not isinstance(W, (list, tuple)): W = [W] if pol: n = 4 if len(W) == 1: W = W * 3 assert len(W) == 3 else: n = len(W) M = np.zeros((n, lmax + 1, lmax + 1)) nW = len(W) lmax_w = min(2 * lmax, len(W[0]) - 1) for m in W[1:]: assert lmax_w == min(2 * lmax, len(m) - 1) Wmat = np.empty((nW, lmax_w + 1)) for i, m in enumerate(W): Wmat[i, :] = m[:lmax_w + 1] _coupling_3j(Wmat, c_int(lmax_w), c_int(nW), c_bool(pol), M, c_int(lmax)) if n == 1: return M[0, :, :] else: return [M[i, :, :] for i in range(n)]
[docs]def scalar_coupling_matrix(P, lmax): """ Get scalar Pseudo-Cl coupling matrix from power spectrum of mask, or array of power masks. Uses multiple threads. See Eq A31 of `astro-ph/0105302 <https://arxiv.org/abs/astro-ph/0105302>`_ :param P: power spectrum of mask, or list of mask power spectra :param lmax: lmax for the matrix (assumed square) :return: coupling matrix (square but not symmetric), or list of couplings for different masks """ if not isinstance(P, (list, tuple)): P = [P] elif any(x.size != P[0].size for x in P[1:]): raise ValueError('Mask power spectra must have same lmax') lmax_power = min(P[0].size - 1, 2 * lmax) if lmax_power < 2 * lmax: print('Warning: power spectrum lmax is less than 2*lmax') fac = (2 * np.arange(lmax_power + 1) + 1) / 4 / np.pi M = threej_coupling([fac * power for power in P], lmax) factor = 2 * np.arange(lmax + 1) + 1 if len(P) == 1: return M * factor else: return [m * factor for m in M]
[docs]def pcl_coupling_matrix(P, lmax, pol=False): """ Get Pseudo-Cl coupling matrix from power spectrum of mask. Uses multiple threads. See Eq A31 of `astro-ph/0105302 <https://arxiv.org/abs/astro-ph/0105302>`_ :param P: power spectrum of mask :param lmax: lmax for the matrix :param pol: whether to calculate TE, EE, BB couplings :return: coupling matrix (square but not symmetric), or list of TT, TE, EE, BB if pol """ lmax_power = min(P.size - 1, 2 * lmax) if lmax_power < 2 * lmax: print('Warning: power spectrum lmax is less than 2*lmax') W = (2 * np.arange(lmax_power + 1) + 1) * P / (4 * np.pi) M = threej_coupling(W, lmax, pol=pol) factor = 2 * np.arange(lmax + 1) + 1 if pol: return [mat * factor for mat in M] else: return M * factor
_gauss_legendre = camblib.__mathutils_MOD_gauss_legendre _gauss_legendre.argtypes = [numpy_1d, numpy_1d, int_arg] def gauss_legendre(xvals, weights, npoints): _gauss_legendre(xvals, weights, c_int(npoints))