visualdynamics.core.partial_coherence

partial_coherence

Partial, reference and virtual coherence of a multiple-input record.

In a test with several drives, the multiple coherence says how much of a response the drives account for together. It does not say which drive does, and it does not say whether the drives were independent enough for the identification to separate them. Three readings of the same cross-spectral matrix answer those questions:

• The reference coherence is the ordinary coherence between every pair of references. Two drives coherent near one cannot be told apart at that line: the reference matrix G_rr is nearly singular there, and an FRF computed from it is whichever of many answers the pseudo-inverse picks.

• The partial coherence of a response with one reference is the ordinary coherence between the two after the linear effect of every other reference has been removed from both: the coherence of the conditioned spectra. With references x_i and the rest r,

  G_ab·r = G_ab − G_ar G_rr⁻¹ G_rb γ²_iy·r = |G_iy·r|² / (G_ii·r G_yy·r)

It says how much of what the others leave of the response this reference explains. With one reference it is the ordinary coherence. For two, 1 − γ²_y:x = (1 − γ²_1y)(1 − γ²_2y·1): the multiple coherence is built from an ordinary and a partial coherence. - The virtual coherence is the coherence of a response with each principal reference: the uncorrelated combinations of the references given by the eigenvectors of G_rr, G_rr = U Λ Uᴴ. The principal references have autospectra Λ and cross spectra Uᴴ G_ry with the response, and their coherences with it sum to the multiple coherence. A principal reference whose eigenvalue is small next to the first is a direction the drives hardly excited, so the count of large eigenvalues is how many independent drives the test really had.

The spectra follow the convention of the rest of the package, G_ab = E[A* B], so G_ba is the conjugate of G_ab and every matrix is Hermitian.

A conditioned autospectrum that is nothing next to the one it came from is a channel the other references explain completely, and its partial coherence is undefined, not a number; so are the coherences of a channel that recorded nothing and of a line holding a value that is not a finite number. A dead channel or a line nobody drove is not a coherence of zero.

References

The partial and multiple coherence functions and the conditioned spectra they are read from are standard; this is where they are set out.

1. Bendat, J. S., & Piersol, A. G. (2010). Random Data: Analysis and Measurement Procedures, 4th ed. Wiley. Multiple-input systems: conditioned spectral density functions, the partial coherence function of an output with one input when the other inputs are removed, and the multiple coherence function written as a product of ordinary and partial coherences.

Classes:

Name	Description
PartialCoherencePoint	The coherences of one response at one line.
PartialCoherence	The partial, reference and virtual coherence of a multiple-input

Functions:

Name	Description
multiple_coherence	The multiple coherence of each response with a set of references,
conditioned	The conditioned cross-spectral matrix of the channels keep with
read_partial_coherence	Read the partial, reference and virtual coherence from a

Classes

PartialCoherencePoint dataclass

PartialCoherencePoint(frequency: float, multiple: float, ordinary: ndarray, partial: ndarray, virtual: ndarray, principal: ndarray)

The coherences of one response at one line.

Attributes:

Name	Type	Description
frequency	float	The line, Hz.
multiple	float	The multiple coherence with every reference.
ordinary, partial, virtual	ndarray	One value per reference (per principal reference for virtual).
principal	ndarray	The principal references' autospectra, largest first, as shares of the first.

PartialCoherence dataclass

PartialCoherence(frequency: ndarray, references: list[str], responses: list[str], multiple: ndarray, ordinary: ndarray, partial: ndarray, virtual: ndarray, reference: ndarray, principal: ndarray, averages: int = 0, notes: list[str] = list())

The partial, reference and virtual coherence of a multiple-input record, at every line.

Attributes:

Name	Type	Description
frequency	ndarray	The lines, Hz.
references, responses	list of str	The channels, by DOF (with the quantity where a DOF carries two).
multiple	ndarray	(responses, lines).
ordinary, partial, virtual	ndarray	(responses, references, lines): each response with each reference alone, with each conditioned on the others, and with each principal reference.
reference	ndarray	(references, references, lines): the ordinary coherence between every pair of references; one on the diagonal of a live channel.
principal	ndarray	(references, lines): the eigenvalues of G_rr, largest first, in the references' squared unit per hertz when they share one.
averages	int	The frames the spectra were averaged over.
notes	list of str	What the numbers should be read with.

Methods:

Name	Description
index	A response's row, by position or by name.
at	The coherences of one response at the line nearest
coupled	The pairs of references whose coherence passes limit at any
effective_references	How many principal references carry at least floor of the
describe	One line about a response: which reference explains most of

Methods:

index

index(response: int | str) -> int

A response's row, by position or by name.

Source code in src/visualdynamics/core/partial_coherence.py

def index(self, response: int | str) -> int:
    """A response's row, by position or by name."""
    if isinstance(response, str):
        if response not in self.responses:
            raise ValueError(f'{response!r} is not a response of this '
                             'reading: ' + ', '.join(self.responses))
        return self.responses.index(response)
    response = int(response)
    if not 0 <= response < len(self.responses):
        raise ValueError(f'no response {response}: the reading has '
                         f'{len(self.responses)}')
    return response

at

at(frequency: float, response: int | str = 0) -> PartialCoherencePoint

The coherences of one response at the line nearest frequency.

Source code in src/visualdynamics/core/partial_coherence.py

def at(self, frequency: float, response: int | str = 0
       ) -> PartialCoherencePoint:
    """The coherences of one response at the line nearest
    `frequency`."""
    k = self.index(response)
    line = nearest_line(self.frequency, frequency)
    first = self.principal[0, line]
    with np.errstate(divide='ignore', invalid='ignore'):
        shares = (self.principal[:, line] / first if first > 0.0
                  else np.full(len(self.references), np.nan))
    return PartialCoherencePoint(
        float(self.frequency[line]), float(self.multiple[k, line]),
        self.ordinary[k, :, line].copy(), self.partial[k, :, line].copy(),
        self.virtual[k, :, line].copy(), shares)

coupled

coupled(limit: float = COUPLED) -> list[tuple[str, str, float, float]]

The pairs of references whose coherence passes limit at any line, as (one, other, share of the lines past it, the largest), the most coupled first.

Source code in src/visualdynamics/core/partial_coherence.py

def coupled(self, limit: float = COUPLED
            ) -> list[tuple[str, str, float, float]]:
    """The pairs of references whose coherence passes `limit` at any
    line, as (one, other, share of the lines past it, the largest),
    the most coupled first."""
    found = []
    q = len(self.references)
    for i in range(q):
        for j in range(i + 1, q):
            values = self.reference[i, j]
            defined = np.isfinite(values)
            if not defined.any():
                continue
            past = values[defined] > limit
            if past.any():
                found.append((self.references[i], self.references[j],
                              float(np.mean(past)),
                              float(np.max(values[defined]))))
    return sorted(found, key=lambda entry: -entry[2])

effective_references

effective_references(floor: float = 0.01) -> ndarray

How many principal references carry at least floor of the first one's power, at every line: the number of drives the test had in effect.

Source code in src/visualdynamics/core/partial_coherence.py

def effective_references(self, floor: float = 0.01) -> np.ndarray:
    """How many principal references carry at least `floor` of the
    first one's power, at every line: the number of drives the test
    had in effect."""
    first = self.principal[0]
    with np.errstate(divide='ignore', invalid='ignore'):
        shares = self.principal / first
    return np.where(first > 0.0, np.sum(shares >= floor, axis=0), 0)

describe

describe(response: int | str = 0) -> str

One line about a response: which reference explains most of it once the others are removed, and the drives' worst coupling.

Source code in src/visualdynamics/core/partial_coherence.py

def describe(self, response: int | str = 0) -> str:
    """One line about a response: which reference explains most of
    it once the others are removed, and the drives' worst
    coupling."""
    k = self.index(response)
    with np.errstate(all='ignore'):
        medians = [np.nanmedian(row) if np.isfinite(row).any()
                   else np.nan for row in self.partial[k]]
    parts = [(f'{self.responses[k]} against {len(self.references)} '
              f'reference{"s" * (len(self.references) != 1)}')]
    if np.isfinite(medians).any():
        best = int(np.nanargmax(medians))
        parts.append(f'most explained by {self.references[best]} '
                     f'(median partial coherence {medians[best]:.2f})')
        # the median is over the lines where it is defined; where it
        # is not is part of the answer, so it is said beside it
        undefined = int(np.sum(np.isnan(self.partial[k]).any(axis=0)))
        if undefined:
            parts.append(f'undefined at {undefined} of '
                         f'{self.partial.shape[-1]} lines')
    else:
        parts.append('no partial coherence defined')
    coupled = self.coupled()
    if coupled:
        one, other, share, _top = coupled[0]
        parts.append(f'{one} and {other} coherent over {COUPLED:g} on '
                     f'{share:.0%} of the lines')
    return ', '.join(parts)

Functions:

multiple_coherence

multiple_coherence(grr: ndarray, gry: ndarray, gyy: ndarray) -> ndarray

The multiple coherence of each response with a set of references, γ² = g_ryᴴ G_rr⁺ g_ry / G_yy, at every line.

The one statement of the multiple coherence: TimeHistory.compute_multiple_coherence and read_partial_coherence both read it here. A pseudo-inverse is used rather than a solve, because two drives can be nearly the same and G_rr is then close to singular. Limited to 0 to 1; not a number where the response's autospectrum is not above zero or a value is not finite.

Parameters:

Name	Type	Description	Default
grr	ndarray	(lines, q, q), the references' cross-spectral matrix.	required
gry	ndarray	(lines, q, p), the cross spectra from each reference to each response, E[R* Y].	required
gyy	ndarray	(p, lines), the responses' autospectra.	required

Returns:

Type	Description
ndarray	(p, lines).

Source code in src/visualdynamics/core/partial_coherence.py

def multiple_coherence(grr: np.ndarray, gry: np.ndarray,
                       gyy: np.ndarray) -> np.ndarray:
    """The multiple coherence of each response with a set of references,
    `γ² = g_ryᴴ G_rr⁺ g_ry / G_yy`, at every line.

    The one statement of the multiple coherence:
    `TimeHistory.compute_multiple_coherence` and `read_partial_coherence`
    both read it here. A pseudo-inverse is used rather than a solve,
    because two drives can be nearly the same and `G_rr` is then close
    to singular. Limited to 0 to 1; not a number where the response's
    autospectrum is not above zero or a value is not finite.

    Parameters
    ----------
    grr : numpy.ndarray
        `(lines, q, q)`, the references' cross-spectral matrix.
    gry : numpy.ndarray
        `(lines, q, p)`, the cross spectra from each reference to each
        response, `E[R* Y]`.
    gyy : numpy.ndarray
        `(p, lines)`, the responses' autospectra.

    Returns
    -------
    numpy.ndarray
        `(p, lines)`.
    """
    grr = np.asarray(grr, dtype=np.complex128)
    gry = np.asarray(gry, dtype=np.complex128)
    gyy = np.real(np.asarray(gyy)).astype(np.float64)
    finite = (np.all(np.isfinite(grr), axis=(1, 2))[np.newaxis, :]
              & np.all(np.isfinite(gry), axis=1).T & np.isfinite(gyy))
    safe_rr = np.where(np.isfinite(grr), grr, 0.0)
    safe_ry = np.where(np.isfinite(gry), gry, 0.0)
    inverse = np.linalg.pinv(safe_rr)
    explained = np.real(np.einsum('lip,lij,ljp->pl', np.conj(safe_ry),
                                  inverse, safe_ry))
    with np.errstate(divide='ignore', invalid='ignore'):
        value = explained / gyy
    return np.where(finite & (gyy > 0.0), np.clip(value, 0.0, 1.0), np.nan)

conditioned

conditioned(matrix: ndarray, keep: Sequence[int], given: Sequence[int]) -> ndarray

The conditioned cross-spectral matrix of the channels keep with the linear effect of the channels given removed from them, G_kk·g = G_kk − G_kg G_gg⁺ G_gk, at every line.

Parameters:

Name	Type	Description	Default
matrix	ndarray	(lines, n, n), Hermitian.	required
keep	sequence of int	Channel indices; given may be empty, and then the block of keep comes back as it is.	required
given	sequence of int	Channel indices; given may be empty, and then the block of keep comes back as it is.	required

Returns:

Type	Description
ndarray	(lines, len(keep), len(keep)).

Source code in src/visualdynamics/core/partial_coherence.py

def conditioned(matrix: np.ndarray, keep: Sequence[int],
                given: Sequence[int]) -> np.ndarray:
    """The conditioned cross-spectral matrix of the channels `keep` with
    the linear effect of the channels `given` removed from them,
    `G_kk·g = G_kk − G_kg G_gg⁺ G_gk`, at every line.

    Parameters
    ----------
    matrix : numpy.ndarray
        `(lines, n, n)`, Hermitian.
    keep, given : sequence of int
        Channel indices; `given` may be empty, and then the block of
        `keep` comes back as it is.

    Returns
    -------
    numpy.ndarray
        `(lines, len(keep), len(keep))`.
    """
    keep, given = list(keep), list(given)
    block = matrix[:, keep][:, :, keep]
    if not given:
        return block.copy()
    gkg = matrix[:, keep][:, :, given]
    ggg = matrix[:, given][:, :, given]
    return block - gkg @ np.linalg.pinv(ggg, hermitian=True) @ \
        np.conj(np.swapaxes(gkg, 1, 2))

read_partial_coherence

read_partial_coherence(frequency: ndarray, matrix: ndarray, references: Sequence[str], responses: Sequence[str], *, averages: int = 0) -> PartialCoherence

Read the partial, reference and virtual coherence from a cross-spectral matrix whose first channels are the references.

Parameters:

Name	Type	Description	Default
frequency	ndarray	The lines, Hz.	required
matrix	ndarray	(lines, q + p, q + p), Hermitian: the q references, then the p responses, G_ab = E[A* B].	required
references	sequence of str	The channels' names, in the matrix's order.	required
responses	sequence of str	The channels' names, in the matrix's order.	required
averages	int	The frames averaged, carried for the notes.	0

Returns:

Type	Description
PartialCoherence

Raises:

Type	Description
ValueError	When there is no reference or no response, or the matrix's shape does not match the names.

Source code in src/visualdynamics/core/partial_coherence.py

def read_partial_coherence(frequency: np.ndarray, matrix: np.ndarray,
                           references: Sequence[str],
                           responses: Sequence[str], *,
                           averages: int = 0) -> PartialCoherence:
    """Read the partial, reference and virtual coherence from a
    cross-spectral matrix whose first channels are the references.

    Parameters
    ----------
    frequency : numpy.ndarray
        The lines, Hz.
    matrix : numpy.ndarray
        `(lines, q + p, q + p)`, Hermitian: the `q` references, then
        the `p` responses, `G_ab = E[A* B]`.
    references, responses : sequence of str
        The channels' names, in the matrix's order.
    averages : int, optional
        The frames averaged, carried for the notes.

    Returns
    -------
    PartialCoherence

    Raises
    ------
    ValueError
        When there is no reference or no response, or the matrix's
        shape does not match the names.
    """
    references, responses = list(references), list(responses)
    q, p = len(references), len(responses)
    matrix = np.asarray(matrix, dtype=np.complex128)
    frequency = np.asarray(frequency, dtype=np.float64)
    if not q or not p:
        raise ValueError('partial coherence needs at least one reference '
                         'and one response')
    if matrix.shape != (len(frequency), q + p, q + p):
        raise ValueError(f'the matrix is {matrix.shape}, not '
                         f'({len(frequency)}, {q + p}, {q + p})')
    lines = len(frequency)
    finite = np.all(np.isfinite(matrix), axis=(1, 2))
    safe = np.where(finite[:, np.newaxis, np.newaxis], matrix, 0.0)
    auto = np.real(np.diagonal(safe, axis1=1, axis2=2)).T      # (n, lines)
    auto = np.where(finite[np.newaxis, :], auto, np.nan)

    reference = np.full((q, q, lines), np.nan)
    for i in range(q):
        for j in range(q):
            reference[i, j] = ordinary_coherence(auto[i], auto[j],
                                                 safe[:, i, j])
    ordinary = np.full((p, q, lines), np.nan)
    for y in range(p):
        for i in range(q):
            ordinary[y, i] = ordinary_coherence(auto[i], auto[q + y],
                                                safe[:, i, q + y])

    multiple = multiple_coherence(safe[:, :q, :q], safe[:, :q, q:],
                                  auto[q:])

    # each reference with the responses, the other references removed
    # from all of them; the coherence of what is left
    partial = np.full((p, q, lines), np.nan)
    for i in range(q):
        others = [j for j in range(q) if j != i]
        keep = [i, *range(q, q + p)]
        block = conditioned(safe, keep, others)
        drive = _floored(block[:, 0, 0], auto[i])
        for y in range(p):
            response = _floored(block[:, 1 + y, 1 + y], auto[q + y])
            partial[y, i] = ordinary_coherence(
                np.where(finite, drive, np.nan), response, block[:, 0, 1 + y])

    # the principal references: G_rr = U Λ Uᴴ, largest first; their
    # cross spectra with a response are Uᴴ g_ry
    values, vectors = np.linalg.eigh(safe[:, :q, :q])
    values, vectors = values[:, ::-1], vectors[:, :, ::-1]
    first = values[:, :1]
    principal_auto = np.where(values > PRINCIPAL_FLOOR * first, values, 0.0)
    cross = np.conj(np.swapaxes(vectors, 1, 2)) @ safe[:, :q, q:]
    virtual = np.full((p, q, lines), np.nan)
    for y in range(p):
        for k in range(q):
            virtual[y, k] = ordinary_coherence(
                np.where(finite, principal_auto[:, k], np.nan), auto[q + y],
                cross[:, k, y])
    principal = np.where(finite[np.newaxis, :], np.maximum(values.T, 0.0),
                         np.nan)

    notes = _notes(references, responses, auto, finite, partial, q)
    return PartialCoherence(frequency, references, responses, multiple,
                            ordinary, partial, virtual, reference, principal,
                            int(averages), notes)