solprop.processing.coherence#

First-order coherence across an ensemble of noise realisations.

The quantity that says whether a supercontinuum is reproducible from shot to shot: the modulus of the normalised cross-correlation of the spectra of independent realisations, at each frequency (J. M. Dudley and S. Coen, Optics Letters 27, 1180 (2002), doi:10.1364/OL.27.001180).

The definition here is the unbiased one over distinct pairs, which is one for identical realisations and zero, in the mean, for independent ones. A sum over every ordered pair, the diagonal i == j included, would carry the self-terms in its numerator and tend to 1/N rather than to zero for a fully incoherent ensemble.

Functions

`band_coherence`(spectra, reference, *[, ...])	Return an ensemble's coherence averaged over the bands of a reference spectrum.
`g12`(spectra)	Return the first-order coherence of an ensemble of complex spectra.
`g12_of`(outputs, *[, plane, z, modes, band])	Return the first-order coherence of an ensemble of saved runs.
`spectral_bands`(weights[, n_bands])	Return the bin indices of equal bands across a spectrum's occupied span.

solprop.processing.coherence.band_coherence(spectra, reference, *, n_bands=32, db_min=-30.0)[source]#

Return an ensemble’s coherence averaged over the bands of a reference spectrum.

The coherence g12() of spectra per bin, averaged over each band’s bins weighted by the reference ensemble’s mean spectrum, \(\sum_b \bar S_b g_{12,b} / \sum_b \bar S_b\), in the n_bands equal bands of spectral_bands() across the reference mean spectrum’s occupied span. A band whose mean spectral energy is below db_min of the largest band’s is not scored. The reference selects the bands and weighs them, so two ensembles judged against one reference (for example two precisions or two solvers with the same seeds) are averaged over the same bins with the same weights. The reference’s own band coherence is band_coherence(reference, reference).

Parameters:

spectranumpy.ndarray: Complex spectra of the ensemble judged, shape (n_realisations, ..., n_x), at least two realisations.
referencenumpy.ndarray: Complex spectra of the reference ensemble, shape (n_reference, ..., n_x) with the same trailing shape. Its mean of \(|\tilde E|^2\) over the realisations, summed over the axes between, is the weight.
n_bandsint, optional: How many bands, 32 by default.
db_minfloat, optional: Bands below this level of the largest band’s energy are not scored, dB; -30 by default.

Returns:

coherencenumpy.ndarray: The weighted mean coherence of each band, shape (n_bands,); nan for a band not scored.
energynumpy.ndarray: The reference mean spectrum’s energy in each band, shape (n_bands,).

Raises:

solprop.exceptions.OutputError: If the trailing shapes differ, or spectra has fewer than two realisations.

Examples

>>> import numpy as np
>>> from solprop.processing.coherence import band_coherence
>>> one = np.ones(64, dtype=complex)
>>> coherence, energy = band_coherence(np.stack([one, one]), np.stack([one, one]))
>>> bool(np.allclose(coherence, 1.0)), energy.shape
(True, (32,))

solprop.processing.coherence.g12(spectra)[source]#

Return the first-order coherence of an ensemble of complex spectra.

\[g_{12}(\omega) = \frac{\left| \sum_{i \neq j} \tilde E_i^{*}(\omega) \tilde E_j(\omega) \right|} {(N - 1) \sum_i \left| \tilde E_i(\omega) \right|^2}\]

evaluated with the identity \(\sum_{i \neq j} \tilde E_i^{*} \tilde E_j = \left| \sum_i \tilde E_i \right|^2 - \sum_i \left| \tilde E_i \right|^2\), which costs one pass over the ensemble rather than \(N^2\) slices. Identical realisations give exactly one; independent ones give zero in the mean, because the self-terms the identity removes are the only ones that survive the average.

Parameters:

spectranumpy.ndarray: Complex spectra of the ensemble, shape (n_realisations, ..., n_x), with at least two realisations along the leading axis.

Returns:

numpy.ndarray: The coherence, shape (...), between 0 and 1. A frequency at which every realisation is empty gives nan.

Raises:

solprop.exceptions.OutputError: If fewer than two realisations are given.

Examples

>>> import numpy as np
>>> from solprop.processing.coherence import g12
>>> one = np.array([1.0 + 0.0j, 2.0 + 0.0j])
>>> g12(np.stack([one, one])).round(12).tolist()
[1.0, 1.0]

solprop.processing.coherence.g12_of(outputs, *, plane=None, z=None, modes=None, band='active')[source]#

Return the first-order coherence of an ensemble of saved runs.

Every run must have been computed on the same grid and saved at the same planes; the spectra are read through solprop.processing.fields.get_spectrum(), so the ensemble may mix memory and HDF5 handlers.

Parameters:

outputssequence of solprop.output.schema.Output: The ensemble, at least two opened handlers.
planeint or None: Index of the plane wanted, negative from the end; with neither plane nor z, every plane is used and the result carries a leading n_z axis.
zfloat or None: Position of the plane wanted, m.
modesint, sequence of int, ‘sum’ or None: Which modes to report; 'sum' is refused, as it is for any complex amplitude.
band{‘active’, ‘all’}: Which part of the stored axis to use.

Returns:

omeganumpy.ndarray: Angular frequency, shape (n_x,), rad/s, ascending.
coherencenumpy.ndarray: The coherence, between 0 and 1, and nan at a frequency where every run of the ensemble is empty.

Raises:

solprop.exceptions.OutputError: If fewer than two runs are given, or their spectra do not have the same shape.

Examples

>>> import numpy as np
>>> from solprop.processing._examples import example_output
>>> from solprop.processing.coherence import g12_of
>>> ensemble = [example_output(), example_output()]
>>> omega, coherence = g12_of(ensemble, plane=-1)
>>> bool(np.all(coherence[np.isfinite(coherence)] == 1.0))
True

solprop.processing.coherence.spectral_bands(weights, n_bands=32)[source]#

Return the bin indices of equal bands across a spectrum’s occupied span.

SEP 33’s spectral metric: the span from the first to the last bin with non-zero weight, cut into n_bands bands of as nearly equal width as the bins allow.

Parameters:

weightsnumpy.ndarray: Non-negative spectral weight per bin, shape (n_x,), zero outside the band of interest (a mean spectrum restricted to the active band).
n_bandsint, optional: How many bands, 32 by default.

Returns:

list of numpy.ndarray: The bins of each band, in order; empty when every weight is zero.

Examples

>>> import numpy as np
>>> from solprop.processing.coherence import spectral_bands
>>> [band.tolist() for band in spectral_bands(np.array([0.0, 1, 1, 1, 1, 0]), 2)]
[[1, 2], [3, 4]]