Functions for random fields.

Random fields

The following functions provide functionality for simulating random fields on the sphere. This is done in the form of HEALPix maps.

Angular power spectra

All functions that process sets of two-point functions expect them as a sequence using the following “Christmas tree” ordering:

In other words, the sequence begins as such:

  • Index 0 describes the auto-correlation of field 0,

  • index 1 describes the auto-correlation of field 1,

  • index 2 describes the cross-correlation of field 1 and field 0,

  • index 3 describes the auto-correlation of field 2,

  • index 4 describes the cross-correlation of field 2 and field 1,

  • index 5 describes the cross-correlation of field 2 and field 0,

  • etc.

In particular, two-point functions for the first \(n\) fields are contained in the first \(T_n = n \, (n + 1) / 2\) entries of the sequence.

To easily generate or iterate over sequences of two-point functions in standard order, see the glass.enumerate_spectra() and glass.spectra_indices() functions.

Preparing inputs

glass.discretized_cls(cls, *, lmax=None, ncorr=None, nside=None)[source]

Apply discretisation effects to angular power spectra.

Depending on the given arguments, this truncates the angular power spectra to lmax, removes all but ncorr correlations between fields, and applies the HEALPix pixel window function of the given nside. If no arguments are given, no action is performed.

Parameters:

  • cls (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – Angular matter power spectra in GLASS ordering.

  • lmax (int | None) – The maximum mode number to truncate the spectra.

  • ncorr (int | None) – The number of correlated fields to keep.

  • nside (int | None) – The resolution parameter for the HEALPix maps.

Returns:

The discretised angular power spectra.

Raises:

ValueError – If the length of the Cls array is not a triangle number.

Return type:

Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]

glass.compute_gaussian_spectra(fields, spectra)[source]

Compute a sequence of Gaussian angular power spectra.

After transformation by fields, the expected two-point statistics should recover spectra when using a band-limited transform [Tessore23].

Parameters:
Returns:

Gaussian angular power spectra for simulation.

Return type:

Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]

glass.solve_gaussian_spectra(fields, spectra)[source]

Solve a sequence of Gaussian angular power spectra.

After transformation by fields, the expected two-point statistics should recover spectra when using a non-band-limited transform [Tessore23].

Parameters:
Returns:

Gaussian angular power spectra for simulation.

Return type:

Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]

Expectations

glass.effective_cls(cls, weights1, weights2=None, *, lmax=None)[source]

Compute effective angular power spectra from weights.

Computes a linear combination of the angular power spectra cls using the factors provided by weights1 and weights2. Additional axes in weights1 and weights2 produce arrays of spectra.

Parameters:

  • cls (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – Angular matter power spectra in GLASS ordering.

  • weights1 (ndarray[tuple[Any, ...], dtype[float64]] | Array | Array) – Weight factors for spectra. The first axis must be equal to the number of fields.

  • weights2 (ndarray[tuple[Any, ...], dtype[float64]] | Array | Array | None) – Second set of weights. If not given, weights1 is used.

  • lmax (int | None) – Truncate the angular power spectra at this mode number. If not given, the longest input in cls will be used.

Returns:

A dictionary of effective angular power spectra, where keys correspond to the leading axes of weights1 and weights2.

Raises:

  • ValueError – If the length of cls is not a triangle number.

  • ValueError – If the shapes of weights1 and weights2 are incompatible.

Return type:

ndarray[tuple[Any, ...], dtype[float64]] | Array | Array

Generating fields

glass.generate(fields, gls, nside, *, ncorr=None, rng=None)[source]

Sample random fields from Gaussian angular power spectra.

Iteratively sample HEALPix maps of transformed Gaussian random fields with the given angular power spectra gls and resolution parameter nside.

The random fields are sampled from Gaussian random fields using the transformations in fields.

The gls array must contain the angular power power spectra of the Gaussian random fields in standard order.

The optional number ncorr limits how many realised fields are correlated. This saves memory, as only ncorr previous fields are kept.

Parameters:

  • fields (Sequence[Transformation]) – Transformations for the random fields.

  • gls (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – Gaussian angular power spectra.

  • nside (int) – Resolution parameter for the HEALPix maps.

  • ncorr (int | None) – Number of correlated fields. If not given, all fields are correlated.

  • rng (Generator | None) – Random number generator. If not given, a default RNG is used.

Yields:

x – Sampled random fields.

Return type:

Iterator[ndarray[tuple[Any, ...], dtype[Any]]]

Lognormal fields

glass.lognormal_fields(shells, shift=None)[source]

Create lognormal fields for radial windows shells. If shifts is given, it must be a callable that returns a lognormal shift (i.e. the scale parameter) at the nominal redshift of each shell.

Parameters:
Returns:

A sequence describing the lognormal fields.

Return type:

Sequence[Lognormal]

GLASS comes with the following functions for setting accurate lognormal shift values:

glass.lognormal_shift_hilbert2011(z)[source]

Lognormal shift of Hilbert et al. (2011) for convergence fields.

Lognormal shift parameter for the weak lensing convergence from the fitting formula of [Hilbert11].

Parameters:

z (float) – Redshift.

Returns:

Lognormal shift.

Return type:

float

Regularisation

When sets of angular power spectra are used to sample random fields, their matrix \(C_\ell^{ij}\) for fixed \(\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.

glass.regularized_spectra(spectra, *, lmax=None, method='nearest', **method_kwargs)[source]

Regularise a set of angular power spectra.

Regularises a complete set spectra of angular power spectra in standard order such that at every angular mode number \(\ell\), the matrix \(C_\ell^{ij}\) is a valid positive semi-definite covariance matrix.

The length of the returned spectra is set by lmax, or the maximum length of the given spectra if lmax is not given. Shorter input spectra are padded with zeros as necessary. Missing spectra can be set to None or, preferably, an empty array.

The method parameter determines how the regularisation is carried out.

Parameters:

  • spectra (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – Spectra to be regularised.

  • lmax (int | None) – Maximum angular mode number. If not given, the maximum is taken from the provided spectra.

  • method (Literal['nearest', 'clip']) – Regularisation method.

  • method_kwargs (float | None)

Return type:

Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]

glass.regularized_spectra(..., method="nearest", tol=None, niter=100)[source]

Compute the (possibly defective) correlation matrices of the given spectra, then find the nearest valid correlation matrices, using the alternating projections algorithm of [Higham02] with tolerance tol for niter iterations. This keeps the diagonals (i.e. auto-correlations) fixed, but requires all of them to be nonnegative.

See also

glass.algorithm.cov_nearest()

Equivalent function for covariance matrices.

glass.algorithm.nearcorr()

Nearest correlation matrix.

Parameters:
Return type:

Sequence[ndarray[tuple[Any, …], dtype[Any]] | Array | Array]

glass.regularized_spectra(..., method="clip", rtol=None)[source]

Clip negative eigenvalues of the spectra’s covariance matrix to zero. This is a simple fix that guarantees positive semi-definite spectra, but can affect the spectra significantly.

See also

glass.algorithm.cov_clip()

Equivalent function for covariance matrices.

Parameters:
Return type:

Sequence[ndarray[tuple[Any, …], dtype[Any]] | Array | Array]

Indexing

glass.getcl(cls, i, j, lmax=None)[source]

Return a specific angular power spectrum from an array in standard order.

Parameters:

  • cls (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – Angular matter power spectra in GLASS ordering.

  • i (int) – Indices to return.

  • j (int) – Indices to return.

  • lmax (int | None) – Truncate the returned spectrum at this mode number.

Returns:

The angular power spectrum for indices i and j from an array in GLASS ordering.

Return type:

ndarray[tuple[Any, ...], dtype[float64]] | Array | Array

glass.enumerate_spectra(entries)[source]

Iterate over a set of two-point functions in standard order, yielding a tuple of indices and their associated entry from the input.

Return type:

Iterator[tuple[int, int, ndarray[tuple[Any, ...], dtype[Any]]]]

Parameters:

entries (Iterable[ndarray[tuple[Any, ...], dtype[Any]]])

Examples

>>> spectra = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
>>> list(enumerate_spectra(spectra))
[(0, 0, [1, 2, 3]), (1, 1, [4, 5, 6]), (1, 0, [7, 8, 9])]
glass.spectra_indices(n)[source]

Return an array of indices in standard order for a set of two-point functions for n fields. Each row is a pair of indices i, j.

Return type:

ndarray[tuple[Any, ...], dtype[integer]]

Parameters:

n (int)

Examples

>>> spectra_indices(3)
array([[0, 0],
       [1, 1],
       [1, 0],
       [2, 2],
       [2, 1],
       [2, 0]])
glass.glass_to_healpix_spectra(spectra)[source]

Reorder spectra from GLASS to HEALPix order.

Reorder spectra in GLASS order to conform to (new) HEALPix order.

Parameters:

spectra (Sequence[TypeVar(T)]) – Sequence of spectra in GLASS order.

Returns:

Sequence of spectra in HEALPix order.

Return type:

list[TypeVar(T)]

glass.healpix_to_glass_spectra(spectra)[source]

Reorder spectra from HEALPix to GLASS order.

Reorder HEALPix spectra (in new order) to conform to GLASS order.

Parameters:

spectra (Sequence[TypeVar(T)]) – Sequence of spectra in HEALPix order.

Returns:

Sequence of spectra in GLASS order.

Return type:

list[TypeVar(T)]

glass.cov_from_spectra(spectra, *, lmax=None)[source]

Construct covariance matrix from spectra.

Construct a covariance matrix from the angular power spectra in spectra.

Parameters:

  • spectra (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – Sequence of angular power spectra.

  • lmax (int | None) – Maximum angular mode number. If not given, the maximum is taken from the provided spectra.

Returns:

Covariance matrix from the given spectra.

Return type:

ndarray[tuple[Any, ...], dtype[Any]]

Deprecated

glass.lognormal_gls(cls, shift=1.0)[source]

Compute Gaussian Cls for a lognormal random field.

Deprecated since version 2025.1: Use glass.lognormal_fields() and glass.compute_gaussian_spectra() or glass.solve_gaussian_spectra() instead.

Parameters:
Returns:

The Gaussian angular power spectra for a lognormal random field.

Return type:

Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]

glass.generate_gaussian(gls, nside, *, ncorr=None, rng=None)[source]

Sample Gaussian random fields from Cls iteratively.

Deprecated since version 2025.1: Use glass.generate() instead.

A generator that iteratively samples HEALPix maps of Gaussian random fields with the given angular power spectra gls and resolution parameter nside.

The optional argument ncorr can be used to artificially limit now many realised fields are correlated. This saves memory, as only ncorr previous fields need to be kept.

The gls array must contain the angular power power spectra of the Gaussian random fields in standard order.

Parameters:

  • gls (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – The Gaussian angular power spectra for a random field.

  • nside (int) – The resolution parameter for the HEALPix maps.

  • ncorr (int | None) – The number of correlated fields. If not given, all fields are correlated.

  • rng (Generator | None) – Random number generator. If not given, a default RNG is used.

Yields:

fields – The Gaussian random fields.

Raises:

ValueError – If all gls are empty.

Return type:

Generator[ndarray[tuple[Any, ...], dtype[float64]]]

glass.generate_lognormal(gls, nside, shift=1.0, *, ncorr=None, rng=None)[source]

Sample lognormal random fields from Gaussian Cls iteratively.

Deprecated since version 2025.1: Use glass.generate() instead.

Parameters:

  • gls (Sequence[ndarray[tuple[Any, ...], dtype[Any]] | Array | Array]) – The Gaussian angular power spectra for a lognormal random field.

  • nside (int) – The resolution parameter for the HEALPix maps.

  • shift (float) – The shift parameter for the lognormal transformation.

  • ncorr (int | None) – The number of correlated fields. If not given, all fields are correlated.

  • rng (Generator | None) – Random number generator. If not given, a default RNG is used.

Yields:

fields – The lognormal random fields.

Return type:

Generator[ndarray[tuple[Any, ...], dtype[float64]]]