implement legacy mode

.gitignore

@@ -11,3 +11,4 @@ dist
 coverage*
 *.html
 .ipynb_checkpoints
+.jupyter

docs/examples.rst

@@ -45,3 +45,4 @@ More advanced examples doing multiple things at the same time.
 
    examples/2-advanced/cosmic_shear.ipynb
    examples/2-advanced/stage_4_galaxies.ipynb
+   examples/2-advanced/legacy-mode.ipynb

docs/modules/algorithm.rst

@@ -0,0 +1,26 @@
+.. module:: glass.core.algorithm
+
+:mod:`glass.core.algorithm` --- General purpose algorithms
+==========================================================
+
+.. currentmodule:: glass.core.algorithm
+
+This module contains general implementations of algorithms which are used by
+GLASS, but are otherwise unrelated to GLASS functionality.
+
+This module must be imported manually::
+
+    import glass.core.algorithm as algo
+
+
+Non-negative least squares
+--------------------------
+
+.. autofunction:: nnls
+
+
+Covariance matrix regularisation
+--------------------------------
+
+.. autofunction:: cov_clip
+.. autofunction:: cov_nearest

docs/modules/index.rst

@@ -5,4 +5,5 @@ Modules
 .. toctree::
    :maxdepth: 2
 
+   algorithm
    grf

docs/reference/fields.rst

@@ -61,13 +61,50 @@ Lognormal fields
 
 .. autofunction:: lognormal_fields
 
+GLASS comes with the following functions for setting accurate lognormal shift
+values:
+
+.. autofunction:: lognormal_shift_hilbert2011
+
+
+Regularisation
+--------------
+
+When sets of angular power spectra are used to sample random fields, their
+matrix :math:`C_\ell^{ij}` for fixed :math:`\ell` must form a valid
+positive-definite covariance matrix.  This is not always the case, for example
+due to numerical inaccuracies, or transformations of the underlying fields
+[Xavier16]_.
+
+Regularisation takes sets of spectra which are ill-posed for sampling, and
+returns sets which are well-defined and, in some sense, "close" to the input.
+
+.. autofunction:: regularized_spectra
+
+.. function:: regularized_spectra(..., method="nearest", niter=20)
+   :no-index:
+
+   Compute the (possibly defective) correlation matrices of the given spectra,
+   then find the nearest valid correlation matrices, using the algorithm from
+   [Higham02]_ for *niter* iterations.  This keeps the diagonals (i.e.
+   auto-correlations) fixed, but requires all of them to be nonnegative.
+
+.. function:: regularized_spectra(..., method="clip")
+   :no-index:
+
+   Compute the eigendecomposition of the given spectra's matrices and set all
+   negative eigenvalues to zero.
+
 
 Indexing
 --------
 
 .. autofunction:: getcl
 .. autofunction:: enumerate_spectra
 .. autofunction:: spectra_indices
+.. autofunction:: glass_to_healpix_spectra
+.. autofunction:: healpix_to_glass_spectra
+.. autofunction:: cov_from_spectra
 
 
 Deprecated

docs/user/bibliography.rst

@@ -0,0 +1,13 @@
+Bibliography
+============
+
+.. [Bro97] Bro, R. and De Jong, S. (1997), A fast non-negativity-constrained
+   least squares algorithm. J.  Chemometrics, 11, 393-401.
+
+.. [Higham02] Higham N. J., 2002, IMA J. Numer. Anal., 22, 329.
+
+.. [Lawson95] Lawson, C. L. and Hanson, R. J. (1995), Solving Least Squares
+   Problems. doi: 10.1137/1.9781611971217
+
+.. [Xavier16] Xavier H. S., et al., 2016, MNRAS, 459, 3693.
+   `doi:10.1093/mnras/stw874 <https://dx.doi.org/10.1093/mnras/stw874>`_

docs/user/index.rst

@@ -7,4 +7,5 @@ User guide
 
    how-glass-works
    publications
+   bibliography
    definitions

glass/__init__.py

@@ -3,9 +3,11 @@
 __all__ = [
     "MultiPlaneConvergence",
     "RadialWindow",
+    "check_posdef_spectra",
     "cls2cov",
     "combine",
     "compute_gaussian_spectra",
+    "cov_from_spectra",
     "cubic_windows",
     "deflect",
     "density_weight",
@@ -29,14 +31,17 @@
     "generate_gaussian",
     "generate_lognormal",
     "getcl",
+    "glass_to_healpix_spectra",
     "grf",
+    "healpix_to_glass_spectra",
     "iternorm",
     "linear_bias",
     "linear_windows",
     "load_cls",
     "loglinear_bias",
     "lognormal_fields",
     "lognormal_gls",
+    "lognormal_shift_hilbert2011",
     "multalm",
     "multi_plane_matrix",
     "multi_plane_weights",
@@ -46,6 +51,7 @@
     "redshift_grid",
     "redshifts",
     "redshifts_from_nz",
+    "regularized_spectra",
     "restrict",
     "save_cls",
     "shear_from_convergence",
@@ -72,8 +78,10 @@
 # modules
 from glass import grf
 from glass.fields import (
+    check_posdef_spectra,
     cls2cov,
     compute_gaussian_spectra,
+    cov_from_spectra,
     discretized_cls,
     effective_cls,
     enumerate_spectra,
@@ -82,10 +90,14 @@
     generate_gaussian,
     generate_lognormal,
     getcl,
+    glass_to_healpix_spectra,
+    healpix_to_glass_spectra,
     iternorm,
     lognormal_fields,
     lognormal_gls,
+    lognormal_shift_hilbert2011,
     multalm,
+    regularized_spectra,
     solve_gaussian_spectra,
     spectra_indices,
     transform_cls,

glass/core/algorithm.py

@@ -16,7 +16,8 @@ def nnls(
     """
     Compute a non-negative least squares solution.
 
-    Implementation of the algorithm due to [1] as described in [2].
+    Implementation of the algorithm due to [Lawson95]_ as described by
+    [Bro97]_.
 
     Parameters
     ----------
@@ -42,14 +43,6 @@ def nnls(
     ValueError
         If the shapes of ``a`` and ``b`` do not match.
 
-    References
-    ----------
-    * [1] Lawson, C. L. and Hanson, R. J. (1995), Solving Least Squares
-          Problems. doi: 10.1137/1.9781611971217
-    * [2] Bro, R. and De Jong, S. (1997), A fast
-          non-negativity-constrained least squares algorithm. J.
-          Chemometrics, 11, 393-401.
-
     """
     a = np.asanyarray(a)
     b = np.asanyarray(b)
@@ -93,3 +86,71 @@ def nnls(
         x[p] = sp
         x[~p] = 0
     return x
+
+
+def cov_clip(cov: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
+    """Covariance matrix from clipping negative eigenvalues."""
+    xp = cov.__array_namespace__()
+
+    # Hermitian eigendecomposition
+    w, v = xp.linalg.eigh(cov)
+
+    # clip negative diagonal values
+    w = xp.clip(w, 0, None)
+
+    # put matrix back together
+    # enforce symmetry
+    v = xp.sqrt(w[..., None, :]) * v
+    return xp.matmul(v, xp.matrix_transpose(v))  # type: ignore[no-any-return]
+
+
+def cov_nearest(
+    cov: npt.NDArray[np.float64],
+    niter: int = 100,
+) -> npt.NDArray[np.float64]:
+    """
+    Covariance matrix from nearest correlation matrix.  Uses the
+    algorithm of [Higham02]_.  The diagonal of the result is unchanged.
+    """
+    xp = cov.__array_namespace__()
+
+    # size of the covariance matrix
+    n = cov.shape[-1]
+
+    # indices of the diagonal of the covariancd matrix
+    diag = (..., xp.arange(n), xp.arange(n))
+
+    # cannot fix negative diagonal
+    if xp.any(cov[diag] < 0):
+        raise ValueError("negative values on the diagonal")
+
+    # store the normalisation of the matrix
+    norm = xp.sqrt(cov[diag])
+    norm = norm[..., None, :] * norm[..., :, None]
+
+    # compute the correlation matrix
+    corr = cov / xp.where(norm > 0, norm, 1.0)
+
+    # initial correction is zero
+    dyks = xp.zeros_like(corr)
+
+    # find the nearest correlation matrix
+    for _ in range(niter):
+        # apply Dykstra's correction to current result
+        delt = corr - dyks
+
+        # project onto positive semi-definite matrices
+        corr = cov_clip(delt)
+
+        # compute Dykstra's correction
+        dyks = corr - delt
+
+        # if correction is zero, the matrix is a correlation matrix
+        if xp.all(dyks == 0):
+            break
+
+        # project onto matrices with unit diagonal
+        corr[diag] = 1
+
+    # apply normalisation to nearest correlation matrix
+    return corr * norm  # type: ignore[no-any-return]