support for Gaussian random fields

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2023 Nicolas Tessore
+Copyright (c) 2023-2024 Nicolas Tessore
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

docs/grf.rst

@@ -0,0 +1,34 @@
+:mod:`angst.grf` --- Gaussian random fields
+===========================================
+
+.. currentmodule:: angst.grf
+.. module:: angst.grf
+
+
+Gaussian angular power spectra
+------------------------------
+
+.. autofunction:: solve
+
+.. class:: Transformation(Protocol)
+
+   .. automethod:: __call__
+   .. automethod:: inv
+   .. automethod:: der
+
+
+Transformations
+---------------
+
+.. class:: Lognormal
+
+   Implements the :class:`Transformation` for lognormal fields.
+
+.. class:: LognormalXNormal
+
+   Implements the :class:`Transformation` for the cross-correlation between
+   :class:`Lognormal` and Gaussian fields.
+
+.. class:: SquaredNormal
+
+   Implements the :class:`Transformation` for squared normal fields.

docs/index.rst

@@ -7,13 +7,19 @@
 
 .. toctree::
    :maxdepth: 2
-   :caption: Contents:
+   :caption: Contents
 
    points
-   spectra
+   twopoint
    core
    glossary
 
+.. toctree::
+   :maxdepth: 2
+   :caption: Modules
+
+   grf
+
 
 Indices and tables
 ==================

docs/spectra.rst → docs/twopoint.rst

@@ -1,22 +1,31 @@
-Angular power spectra
-=====================
+Two-point functions
+===================
 
 .. currentmodule:: angst
 
 
-Sets of angular power spectra
------------------------------
+Spectra and correlation functions
+---------------------------------
 
-.. autofunction:: spectra_indices
-.. autofunction:: enumerate_spectra
+.. autofunction:: cl2corr
+.. autofunction:: corr2cl
 
+.. autofunction:: cl2var
 
-.. _spectra_order:
+
+Sets of two-point functions
+---------------------------
+
+.. autofunction:: indices2
+.. autofunction:: enumerate2
+
+
+.. _twopoint_order:
 
 Standard order
 --------------
 
-All functions that process sets of angular power spectra expect them as a
+All functions that process sets of two-point functions expect them as a
 sequence using the following "Christmas tree" ordering:
 
 .. raw:: html
@@ -32,9 +41,8 @@ In other words, the sequence begins as such:
 * index 5 describes the cross-correlation of field 2 and field 0,
 * etc.
 
-In particular, cross-correlations for the first :math:`n` fields are contained
+In particular, two-point functions for the first :math:`n` fields are contained
 in the first :math:`T_n = n \, (n + 1) / 2` entries of the sequence.
 
-To easily generate or iterate over sequences of angular power spectra in
-standard order, see the :func:`enumerate_spectra` and :func:`spectra_indices`
-functions.
+To easily generate or iterate over sequences of two-point functions in standard
+order, see the :func:`enumerate2` and :func:`indices2` functions.

pyproject.toml

@@ -18,6 +18,7 @@ classifiers = [
 requires-python = ">=3.8"
 dependencies = [
     "numpy>=1.20.0",
+    "flt>=2022.7.27",
 ]
 dynamic = ["version"]
 

src/angst/__init__.py

@@ -1,26 +1,38 @@
 """angst -- angular statistics"""
 
 __all__ = [
+    "cl2corr",
+    "cl2var",
+    "corr2cl",
     "displace",
     "displacement",
-    "enumerate_spectra",
+    "enumerate2",
+    "grf",
     "inv_triangle_number",
-    "spectra_indices",
+    "indices2",
     "__version__",
     "__version_tuple__",
 ]
 
+from . import grf
+
 from ._core import (
     inv_triangle_number,
 )
+
 from ._points import (
     displace,
     displacement,
 )
-from ._spectra import (
-    enumerate_spectra,
-    spectra_indices,
+
+from ._twopoint import (
+    cl2corr,
+    cl2var,
+    corr2cl,
+    enumerate2,
+    indices2,
 )
+
 from ._version import (
     __version__,
     __version_tuple__,

src/angst/_twopoint.py

@@ -0,0 +1,147 @@
+"""
+Module for two-point functions.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+# typing
+from typing import Any, Iterable, Iterator
+from numpy.typing import ArrayLike, NDArray
+
+
+def enumerate2(
+    entries: Iterable[ArrayLike | None],
+) -> Iterator[tuple[int, int, ArrayLike | None]]:
+    """
+    Iterate over a set of two-point functions in :ref:`standard order
+    <twopoint_order>`, returning a tuple of indices and their associated entry
+    from the input.
+
+    >>> spectra = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
+    >>> list(angst.enumerate2(spectra))
+    [(0, 0, [1, 2, 3]), (1, 1, [4, 5, 6]), (1, 0, [7, 8, 9])]
+
+    """
+
+    for k, cl in enumerate(entries):
+        i = int((2 * k + 0.25) ** 0.5 - 0.5)
+        j = i * (i + 3) // 2 - k
+        yield i, j, cl
+
+
+def indices2(n: int) -> Iterator[tuple[int, int]]:
+    """
+    Return an iterator over indices in :ref:`standard order <twopoint_order>`
+    for a set of two-point functions for *n* fields.  Each item is a tuple of
+    indices *i*, *j*.
+
+    >>> list(angst.indices2(3))
+    [(0, 0), (1, 1), (1, 0), (2, 2), (2, 1), (2, 0)]
+
+    """
+
+    for i in range(n):
+        for j in range(i, -1, -1):
+            yield i, j
+
+
+def cl2corr(cl: NDArray[Any], closed: bool = False) -> NDArray[Any]:
+    r"""transform angular power spectrum to 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.
+
+    """
+
+    from flt import idlt  # type: ignore [import-not-found]
+
+    # 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
+    c = np.arange(1, 2 * n + 1, 2, dtype=float)
+    c /= 4 * np.pi
+    c *= cl
+
+    # perform the inverse DLT
+    corr: NDArray[Any] = idlt(c, closed=closed)
+
+    # done
+    return corr
+
+
+def corr2cl(corr: NDArray[Any], closed: bool = False) -> NDArray[Any]:
+    r"""transform angular correlation function to 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.
+    closed : bool
+        Compute correlation function over open (``closed=False``) or closed
+        (``closed=True``) interval.
+
+    Returns
+    -------
+    cl : (n,) array_like
+        Angular power spectrum from :math:`0` to :math:`\mathtt{lmax}`.
+
+    """
+
+    from flt import dlt
+
+    # length n of the transform
+    if corr.ndim != 1:
+        raise TypeError("corr must be 1d array")
+    n = corr.shape[-1]
+
+    # compute the DLT coefficients
+    cl: NDArray[Any] = dlt(corr, closed=closed)
+
+    # DLT coefficients = (2l+1)/(4pi) * Cl
+    cl /= np.arange(1, 2 * n + 1, 2, dtype=float)
+    cl *= 4 * np.pi
+
+    # done
+    return cl
+
+
+def cl2var(cl: NDArray[Any]) -> float:
+    """
+    Compute the variance of the spherical random field in a point from the
+    given angular power spectrum.  The input can be multidimensional, with
+    the last axis representing the modes.
+    """
+    ell = np.arange(np.shape(cl)[-1])
+    return np.sum((2 * ell + 1) / (4 * np.pi) * cl)  # type: ignore