transformcl/transformcl.py at main · glass-dev/transformcl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
"""Transform angular power spectra."""

__all__ = [
    "cl",
    "corr",
    "theta",
    "var",
]

from array_api_compat import array_namespace

import flt


def corr(cl, closed=False):
    r"""
    Transform angular power spectrum to angular correlation function.

    Takes an angular power spectrum with :math:`\mathtt{n} =
    \mathtt{lmax}+1` coefficients and returns the corresponding angular
    correlation function in :math:`\mathtt{n}` points.

    The correlation function values can be computed either over the
    closed interval :math:`[0, \pi]`, in which case :math:`\theta_0 = 0`
    and :math:`\theta_{n-1} = \pi`, or over the open interval :math:`(0,
    \pi)`.

    Parameters
    ----------
    cl : (n,) array_like
        Angular power spectrum from :math:`0` to :math:`\mathtt{lmax}`.
    closed : bool
        Compute correlation function over open (``closed=False``) or closed
        (``closed=True``) interval.

    Returns
    -------
    corr : (n,) array_like
        Angular correlation function.

    See Also
    --------
    transformcl.cl :
        the inverse operation
    transformcl.theta :
        angles at which the correlation function is evaluated

    """

    xp = array_namespace(cl)

    # length n of the transform
    if cl.ndim != 1:
        raise TypeError("cl must be 1d array")
    n = cl.shape[-1]

    # DLT coefficients = (2l+1)/(4pi) * Cl
    a = (2 * xp.arange(n) + 1) / (4 * xp.pi) * cl

    return flt.idlt(a, closed)


def cl(corr, closed=False):
    r"""
    Transform angular correlation function to angular power spectrum.

    Takes an angular function in :math:`\mathtt{n}` points and returns
    the corresponding angular power spectrum from :math:`0` to
    :math:`\mathtt{lmax} = \mathtt{n}-1`.

    The correlation function must be given at the angles returned by
    :func:`transformcl.theta`.  These can be distributed either over the
    closed interval :math:`[0, \pi]`, in which case :math:`\theta_0 = 0`
    and :math:`\theta_{n-1} = \pi`, or over the open interval :math:`(0,
    \pi)`.

    Parameters
    ----------
    corr : (n,) array_like
        Angular correlation function.

    Returns
    -------
    cl : (n,) array_like
        Angular power spectrum from :math:`0` to :math:`\mathtt{lmax}`.
    closed : bool
        Compute correlation function over open (``closed=False``) or
        closed (``closed=True``) interval.

    See Also
    --------
    transformcl.corr :
        the inverse operation
    transformcl.theta :
        angles at which the correlation function is evaluated

    """

    xp = array_namespace(corr)

    # length n of the transform
    if corr.ndim != 1:
        raise TypeError("corr must be 1d array")
    n = corr.shape[-1]

    # DLT coefficients = (2l+1)/(4pi) * Cl
    fl = (2 * xp.arange(n) + 1) / (4 * xp.pi)

    return flt.dlt(corr, closed) / fl


def var(cl):
    r"""
    Compute variance from angular power spectrum.

    Given the angular power spectrum, compute the variance of the
    spherical random field in a point.

    Parameters
    ----------
    cl : array_like
        Angular power spectrum.  Can be multidimensional, with the last
        axis representing the modes.

    Returns
    -------
    var: float
        The variance of the given power spectrum.

    Notes
    -----
    The variance :math:`\sigma^2` of the field with power spectrum
    :math:`C_l` is

    .. math::

        \sigma^2 = \sum_{l} \frac{2l + 1}{4\pi} \, C_l \;.

    """
    xp = array_namespace(cl)
    # ell cannot be an integer here as, within the array api
    # only floating-point dtypes are allowed in __truediv__
    ell = xp.arange(cl.shape[-1], dtype=xp.float64)
    return xp.sum((2 * ell + 1) / (4 * xp.pi) * cl, axis=-1)


def theta(n, closed=False):
    r"""
    Return the angles :math:`\theta_1, \ldots, \theta_n` of the
    correlation function with *n* points.
    """
    return flt.theta(n, closed)


cltocorr = corr
corrtocl = cl
cltovar = var