`visualdynamics.core.matrix_checks`¶

matrix_checks ¶

Two checks on the matrices a multi-input test is built from, made at every frequency line before anything is inverted or played.

A multi-input test stands on two matrices per line. The plant matrix H, one row per control channel and one column per drive, is what the controller inverts to turn the control error into a drive correction. The cross-spectral density matrix S, one row and one column per channel, is what a specification asks the control channels to hold: their autospectra on the diagonal and, off it, how each pair moves together. Each can be wrong in a way no single record shows.

Conditioning

The singular values of H at a line are its gains in the combinations of drives it responds to most and least, largest first. Their ratio is the condition number,

kappa = sigma_max / sigma_min,

dimensionless. It bounds how much a relative error in the measured plant can grow in its inverse: an estimate good to one percent, at a condition number of 100, can give a drive wrong by up to a hundred percent. Two drives that reach the control channels in nearly the same proportions, such as two shakers under one stiff fixture below its first flexible mode, make two nearly equal columns, a small sigma_min and a large condition number, and the inverse that separates them is made of noise. CONDITION_LIMIT is the threshold a line is flagged over, a choice and not a law; a controller holds a setting of its own for it.

The effective rank at that threshold is the number of singular values at least sigma_max / limit, so a line is flagged exactly where the effective rank falls under the full rank. The same reading is made of the cross-spectral matrix of the drives, Gff, whose singular values are its eigenvalues: drives that are nearly copies of one another make it nearly singular, and an FRF estimate from them inverts it.

read_conditioning reads an FRF (several references) or a CSD and returns a Conditioning. The matrix is in one unit or the reading is refused: a condition number of a matrix whose rows are in different units changes when one row's unit does, so it would describe the units and not the plant.

Realizability

The cross-spectral matrix of any set of signals is Hermitian (the term for channel j against i is the conjugate of i against j) and positive semi-definite (no combination of the channels has negative power) at every line. A specification that is not cannot be produced by any drive, and a controller asked for one fails at those lines, however long it runs. Each entry can look right on its own: three channels with every coherence 1 and phases of 0, 0 and 180 degrees ask for 1 in step with 2, 2 in step with 3, and 1 opposite 3, which no motion does.

The check is made on the coherence matrix, the cross-spectral matrix scaled by the root of its diagonal,

C_ij = S_ij / sqrt(S_ii S_jj),

which has ones on its diagonal, is free of units, and is positive semi-definite exactly where S is (a scaling by a positive diagonal does not change the signs of the eigenvalues). Its eigenvalues sum to the number of channels. A line is realizable where the smallest is no lower than minus tolerance times the largest and the matrix is Hermitian to the same tolerance. The default REALIZABLE_TOLERANCE passes a matrix written to six significant figures and refuses any departure larger than that. The rank at a line counts the eigenvalues over tolerance times the largest: one for fully coherent channels, the number of channels for independent ones.

A matrix that fails is repaired to the nearest one that is realizable, in the Frobenius norm of the coherence matrix. Keeping the diagonal (the default), the autospectra stay exactly as written and only the cross terms move: the nearest unit-diagonal positive semi-definite matrix, found by alternating projections with Dykstra's correction. Not keeping it, the negative eigenvalues are set to zero and the matrix rebuilt, which also moves the autospectra. Lines that pass are left as they are, bit for bit. The reading reports the largest change in coherence and in phase the repair made, and where.

read_realizability reads a CSD and returns a Realizability; SpecificationDraft refuses a sheet whose coherence and phase entries cannot hold together, by the same check.

References

The singular value decomposition and the condition number are Golub and Van Loan's; the nearest positive semi-definite matrix is Higham's, in two papers; the properties of a cross-spectral density matrix are Bendat and Piersol's.

Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations, 4th ed. Johns Hopkins University Press. The singular value decomposition, the condition number, and the bound it puts on the error of a solution.
Higham, N. J. (1988). "Computing a nearest symmetric positive semidefinite matrix." Linear Algebra and its Applications, 103, 103-118. doi:10.1016/0024-3795(88)90223-6 The nearest positive semi-definite matrix in the Frobenius norm, by setting the negative eigenvalues to zero.
Higham, N. J. (2002). "Computing the nearest correlation matrix: a problem from finance." IMA Journal of Numerical Analysis, 22(3), 329-343. doi:10.1093/imanum/22.3.329 The alternating projections, with Dykstra's correction, that find the nearest unit-diagonal positive semi-definite matrix.
Bendat, J. S., & Piersol, A. G. (2010). Random Data: Analysis and Measurement Procedures, 4th ed. Wiley. The cross-spectral density matrix of several records, Hermitian and non-negative definite, and the ordinary and multiple coherence read from it.

Classes:

Name	Description
`ConditioningPoint`	The conditioning of the matrix at one line: what
`Conditioning`	The singular values and the condition number of a matrix at every
`RealizabilityPoint`	The realizability of the matrix at one line: what
`Realizability`	Whether a cross-spectral matrix can be produced by any set of

Functions:

Name	Description
`channel_labels`	A label per (DOF, quantity) channel: the DOF, with the quantity
`runs`	(first, last) index of each run of consecutive true flags.
`nearest_line`	The index of the line of `frequency` nearest `value`.
`frf_matrix`	The records of an FRF assembled into the response-by-reference
`entry_channels`	(row channel, column channel) of one record of a spectral matrix.
`matrix_entries`	Each record of a spectral matrix with the entry it holds.
`mirrors_disagree`	Whether a cross term and its mirror, read at the same lines, are
`spectral_matrix`	The records of a CSD assembled into its square matrix, at every
`flagged_runs`	(first, last) frequency of each run of consecutive flagged
`condition_matrix`	The singular values of a matrix at every line, as a
`read_conditioning`	The singular values and condition number of an FRF matrix, or of
`repair_coherence`	The nearest realizable matrix to each Hermitian coherence matrix,
`judge_lines`	Whether a cross-spectral matrix is realizable at every line, with
`check_spectral_matrix`	Whether a cross-spectral matrix is realizable at every line, and
`read_realizability`	Whether a CSD (a specification, or a measured cross-spectral
`made_realizable`	A copy of a CSD with the matrix the picked records make repaired
`entries_refusal`	Why a specification sheet's coherence and phase entries cannot

Classes¶

ConditioningPoint `dataclass` ¶

ConditioningPoint(frequency: float, singular: tuple[float, ...], condition: float, rank: int, size: int, limit: float)

The conditioning of the matrix at one line: what Conditioning.at returns.

singular holds the singular values, largest first, in the reading's units (SI); condition is the largest over the smallest, infinite where the smallest is zero and not a number where the line was not judged; rank is the effective rank at limit, of a full rank of size.

Methods:

Name	Description
`describe`	One line for a status bar: the condition number against the

Attributes:

Name	Type	Description
`ill`	`bool`	Whether the condition number is over the limit, or the

Attributes¶

ill `property` ¶

ill: bool

Whether the condition number is over the limit, or the smallest singular value is zero.

Methods:¶

describe ¶

describe(unit: str = '', scale: float = 1.0) -> str

One line for a status bar: the condition number against the limit, the effective rank, and the singular values in unit (scale converts the reading's SI into it).

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

def describe(self, unit: str = '', scale: float = 1.0) -> str:
    """One line for a status bar: the condition number against the
    limit, the effective rank, and the singular values in `unit`
    (`scale` converts the reading's SI into it)."""
    if np.isnan(self.condition):
        return f'{self.frequency:.4g} Hz: not judged'
    unit = f' {unit}' if unit else ''
    shown = ', '.join(f'{value * scale:.4g}' for value in self.singular)
    verdict = 'over' if self.ill else 'within'
    return (f'{self.frequency:.4g} Hz: condition number '
            f'{self.condition:.4g}, {verdict} {self.limit:g}; rank '
            f'{self.rank} of {self.size}; singular values {shown}{unit}')

Conditioning `dataclass` ¶

Conditioning(frequency: ndarray, singular: ndarray, rows: tuple[str, ...], columns: tuple[str, ...], limit: float = CONDITION_LIMIT, dimension: str = UNKNOWN, kind: str = 'FRF', notes: tuple[str, ...] = ())

The singular values and the condition number of a matrix at every line: what read_conditioning returns.

frequency holds the lines in Hz and singular the singular values (one row per value, largest first; not a number on a line that was not judged), in the matrix's units (SI, dimension). rows and columns name the matrix's channels, responses and references for an FRF, the channels twice for a CSD; kind is 'FRF' or 'CSD'. limit is the condition number lines are flagged over. notes says what the reading should be read with.

Methods:

Name	Description
`bands`	The runs of flagged lines, as (first, last) frequencies.
`line`	The index of the line nearest `frequency`.
`at`	The reading at the line nearest `frequency`.
`describe`	The verdict in a sentence: where the condition number is over

Attributes:

Name	Type	Description
`size`	`int`	The full rank, the smaller of the rows and the columns.
`judged`	`ndarray`	The lines that were judged: every entry a finite number and
`condition`	`ndarray`	The condition number at every line, the largest singular value
`rank`	`ndarray`	The effective rank at every line: how many singular values are
`ill`	`ndarray`	The lines flagged: judged, with a condition number over the
`title`	`str`	What the matrix is, in words.

Attributes¶

size `property` ¶

size: int

The full rank, the smaller of the rows and the columns.

judged `property` ¶

judged: ndarray

The lines that were judged: every entry a finite number and the matrix not all zero.

condition `property` ¶

condition: ndarray

The condition number at every line, the largest singular value over the smallest: infinite where the smallest is zero, not a number where the line was not judged.

rank `property` ¶

rank: ndarray

The effective rank at every line: how many singular values are at least the largest over limit. Zero on a line not judged.

ill `property` ¶

ill: ndarray

The lines flagged: judged, with a condition number over the limit (where the effective rank is under the full rank).

title `property` ¶

title: str

What the matrix is, in words.

Methods:¶

bands ¶

bands() -> list[tuple[float, float]]

The runs of flagged lines, as (first, last) frequencies.

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

def bands(self) -> list[tuple[float, float]]:
    """The runs of flagged lines, as (first, last) frequencies."""
    return flagged_runs(self.frequency, self.ill)

line ¶

line(frequency: float) -> int

The index of the line nearest frequency.

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

def line(self, frequency: float) -> int:
    """The index of the line nearest `frequency`."""
    return nearest_line(self.frequency, frequency)

at ¶

at(frequency: float) -> ConditioningPoint

The reading at the line nearest frequency.

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

def at(self, frequency: float) -> ConditioningPoint:
    """The reading at the line nearest `frequency`."""
    k = self.line(frequency)
    return ConditioningPoint(
        frequency=float(self.frequency[k]),
        singular=tuple(float(v) for v in self.singular[:, k]),
        condition=float(self.condition[k]), rank=int(self.rank[k]),
        size=self.size, limit=float(self.limit))

describe ¶

describe() -> str

The verdict in a sentence: where the condition number is over the limit, or that it is nowhere, or that nothing was judged.

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

def describe(self) -> str:
    """The verdict in a sentence: where the condition number is over
    the limit, or that it is nowhere, or that nothing was judged."""
    judged = int(self.judged.sum())
    if not judged:
        return f'{self.title}: no line could be judged'
    runs = self.bands()
    worst = float(np.nanmax(np.where(self.judged, self.condition,
                                     np.nan)))
    if not runs:
        return (f'{self.title}: condition number within '
                f'{self.limit:g} at all {judged} lines judged '
                f'(largest {worst:.4g}){_unjudged(self.judged)}')
    flagged = int(self.ill.sum())
    return (f'{self.title}: condition number over {self.limit:g} at '
            f'{flagged} of {judged} lines, {_span_text(runs)}'
            f'{_unjudged(self.judged)}')

RealizabilityPoint `dataclass` ¶

RealizabilityPoint(frequency: float, eigenvalues: tuple[float, ...], asymmetry: float, rank: int, realizable: bool, tolerance: float)

The realizability of the matrix at one line: what Realizability.at returns. eigenvalues are those of the coherence matrix, smallest first; asymmetry is its largest departure from Hermitian; realizable says whether the line passes.

Methods:

Name	Description
`describe`	One line for a status bar.

Methods:¶

describe ¶

describe() -> str

One line for a status bar.

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

def describe(self) -> str:
    """One line for a status bar."""
    if not self.eigenvalues or np.isnan(self.eigenvalues[0]):
        return f'{self.frequency:.4g} Hz: not judged'
    verdict = 'realizable' if self.realizable else 'not realizable'
    return (f'{self.frequency:.4g} Hz: {verdict}; smallest eigenvalue '
            f'{self.eigenvalues[0]:.4g} against largest '
            f'{self.eigenvalues[-1]:.4g}, rank {self.rank} of '
            f'{len(self.eigenvalues)}')

Realizability `dataclass` ¶

Realizability(frequency: ndarray, coherence: ndarray, scale: ndarray, eigenvalues: ndarray, asymmetry: ndarray, negative: ndarray, repaired: ndarray, channels: tuple[str, ...], tolerance: float = REALIZABLE_TOLERANCE, keep_diagonal: bool = True, notes: tuple[str, ...] = ())

Whether a cross-spectral matrix can be produced by any set of signals, line by line, and the nearest one that can: what read_realizability returns.

frequency holds the lines in Hz. coherence is the coherence matrix at every line (lines x n x n, the cross-spectral matrix scaled by the root of its diagonal) and scale that root (lines x n), so the cross-spectral matrix is recovered as scale_i scale_j coherence_ij. eigenvalues (n x lines, smallest first) are those of the coherence matrix's Hermitian part; asymmetry is its largest departure from Hermitian at each line; negative marks lines where an autospectrum is below zero. repaired is the coherence matrix after the repair, the same as coherence on every line that passes. channels names the rows; tolerance is the check's.

Methods:

Name	Description
`bands`	The runs of lines that are not realizable, as (first, last)
`matrix`	The cross-spectral matrix at every line (lines x n x n), as
`pair`	(coherence before, after, phase before, after) of channels `i`
`worst`	(i, j, line) of the cross term the repair moved most, as a
`changes`	(coherence, phase in degrees): the largest change the repair
`line`	The index of the line nearest `frequency`.
`at`	The reading at the line nearest `frequency`.
`describe`	The verdict in a sentence: realizable, or the bands that are

Attributes:

Name	Type	Description
`judged`	`ndarray`	The lines judged: every entry a finite number.
`smallest`	`ndarray`	The smallest eigenvalue at every line.
`largest`	`ndarray`	The largest eigenvalue at every line.
`failing`	`ndarray`	The lines judged and not realizable: an eigenvalue below minus
`realizable`	`bool`	Whether every line was judged realizable. False when no line
`rank`	`ndarray`	The rank at every line: the eigenvalues over the tolerance

Attributes¶

judged `property` ¶

judged: ndarray

The lines judged: every entry a finite number.

smallest `property` ¶

smallest: ndarray

The smallest eigenvalue at every line.

largest `property` ¶

largest: ndarray

The largest eigenvalue at every line.

failing `property` ¶

failing: ndarray

The lines judged and not realizable: an eigenvalue below minus the tolerance times the largest, a departure from Hermitian over the tolerance, or a negative autospectrum.

realizable `property` ¶

realizable: bool

Whether every line was judged realizable. False when no line could be judged: a check that judged nothing has not passed.

rank `property` ¶

rank: ndarray

The rank at every line: the eigenvalues over the tolerance times the largest. Zero on a line not judged.

Methods:¶

bands ¶

bands() -> list[tuple[float, float]]

The runs of lines that are not realizable, as (first, last) frequencies.

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

def bands(self) -> list[tuple[float, float]]:
    """The runs of lines that are not realizable, as (first, last)
    frequencies."""
    return flagged_runs(self.frequency, self.failing)

matrix ¶

matrix(repaired: bool = False) -> ndarray

The cross-spectral matrix at every line (lines x n x n), as read or, with repaired, after the repair.

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

def matrix(self, repaired: bool = False) -> np.ndarray:
    """The cross-spectral matrix at every line (lines x n x n), as
    read or, with `repaired`, after the repair."""
    c = self.repaired if repaired else self.coherence
    return c * self.scale[:, :, None] * self.scale[:, None, :]

pair ¶

pair(i: int, j: int) -> tuple[ndarray, ndarray, ndarray, ndarray]

(coherence before, after, phase before, after) of channels i and j at every line, the phase in degrees: what the repair did to that pair. Coherence is data.ordinary_coherence of the matrix before and after.

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

def pair(self, i: int, j: int
         ) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """(coherence before, after, phase before, after) of channels `i`
    and `j` at every line, the phase in degrees: what the repair did
    to that pair. Coherence is `data.ordinary_coherence` of the
    matrix before and after."""
    from .data import ordinary_coherence

    out = []
    for c in (self.coherence, self.repaired):
        out.append(ordinary_coherence(c[:, i, i], c[:, j, j], c[:, i, j]))
    before, after = self.coherence[:, i, j], self.repaired[:, i, j]
    return (out[0], out[1], np.degrees(np.angle(before)),
            np.degrees(np.angle(after)))

worst ¶

worst() -> tuple[int, int, int] | None

(i, j, line) of the cross term the repair moved most, as a complex coherence entry, or None when nothing moved.

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

def worst(self) -> tuple[int, int, int] | None:
    """(i, j, line) of the cross term the repair moved most, as a
    complex coherence entry, or None when nothing moved."""
    return self._worst

changes ¶

changes() -> tuple[float, float]

(coherence, phase in degrees): the largest change the repair made to any pair's coherence and to any phase, over lines where the pair's coherence is not zero after (a phase is undefined without it).

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

def changes(self) -> tuple[float, float]:
    """(coherence, phase in degrees): the largest change the repair
    made to any pair's coherence and to any phase, over lines where
    the pair's coherence is not zero after (a phase is undefined
    without it)."""
    return self._changes

line ¶

line(frequency: float) -> int

The index of the line nearest frequency.

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

def line(self, frequency: float) -> int:
    """The index of the line nearest `frequency`."""
    return nearest_line(self.frequency, frequency)

at ¶

at(frequency: float) -> RealizabilityPoint

The reading at the line nearest frequency.

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

def at(self, frequency: float) -> RealizabilityPoint:
    """The reading at the line nearest `frequency`."""
    k = self.line(frequency)
    return RealizabilityPoint(
        frequency=float(self.frequency[k]),
        eigenvalues=tuple(float(v) for v in self.eigenvalues[:, k]),
        asymmetry=float(self.asymmetry[k]), rank=int(self.rank[k]),
        realizable=bool(self.judged[k] and not self.failing[k]),
        tolerance=float(self.tolerance))

describe ¶

describe() -> str

The verdict in a sentence: realizable, or the bands that are not and the pair the repair moves most, or that nothing was judged.

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

def describe(self) -> str:
    """The verdict in a sentence: realizable, or the bands that are
    not and the pair the repair moves most, or that nothing was
    judged."""
    judged = int(self.judged.sum())
    n = len(self.channels)
    if not judged:
        return f'{n} channels: no line could be judged'
    if self.realizable:
        return (f'{n} channels: realizable at all {judged} lines'
                f'{_unjudged(self.judged)}')
    flagged = int(self.failing.sum())
    text = (f'{n} channels: not realizable at {flagged} of {judged} '
            f'lines, {_span_text(self.bands())}'
            f'{_unjudged(self.judged)}')
    worst = self.worst()
    if worst is not None:
        i, j, line = worst
        text += (f'; the repair moves {self.channels[i]} and '
                 f'{self.channels[j]} most, at '
                 f'{self.frequency[line]:.4g} Hz')
    return text

Functions:¶

channel_labels ¶

channel_labels(channels: Sequence[tuple[str, str]]) -> list[str]

A label per (DOF, quantity) channel: the DOF, with the quantity added only where one DOF carries two channels (an accelerometer and a force gauge at one drive point). The matrix checks, the matrix comparison and the drive prediction all name their channels here.

Parameters:

Name	Type	Description	Default
`channels`	`sequence of (str, str)`	The channels, as (DOF, quantity).	required

Returns:

Type	Description
`list of str`	One label per channel.

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

def channel_labels(channels: Sequence[tuple[str, str]]) -> list[str]:
    """A label per (DOF, quantity) channel: the DOF, with the quantity
    added only where one DOF carries two channels (an accelerometer and
    a force gauge at one drive point). The matrix checks, the matrix
    comparison and the drive prediction all name their channels here.

    Parameters
    ----------
    channels : sequence of (str, str)
        The channels, as (DOF, quantity).

    Returns
    -------
    list of str
        One label per channel.
    """
    dofs = [dof for dof, _quantity in channels]
    return [dof if dofs.count(dof) == 1 else f'{dof} {quantity}'
            for dof, quantity in channels]

runs ¶

runs(flags: Any) -> list[tuple[int, int]]

(first, last) index of each run of consecutive true flags.

Parameters:

Name	Type	Description	Default
`flags`	`array_like`	Booleans, one per line.	required

Returns:

Type	Description
`list of (int, int)`	Each run's first and last index, inclusive.

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

def runs(flags: Any) -> list[tuple[int, int]]:
    """(first, last) index of each run of consecutive true flags.

    Parameters
    ----------
    flags : array_like
        Booleans, one per line.

    Returns
    -------
    list of (int, int)
        Each run's first and last index, inclusive.
    """
    flags = np.asarray(flags, dtype=bool).ravel()
    if not flags.any():
        return []
    edges = np.diff(np.concatenate([[0], flags.astype(np.int8), [0]]))
    return [(int(a), int(b) - 1) for a, b in
            zip(np.flatnonzero(edges == 1), np.flatnonzero(edges == -1))]

nearest_line ¶

nearest_line(frequency: Any, value: float) -> int

The index of the line of frequency nearest value.

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

def nearest_line(frequency: Any, value: float) -> int:
    """The index of the line of `frequency` nearest `value`."""
    return int(np.argmin(np.abs(np.asarray(frequency, dtype=float)
                                - float(value))))

frf_matrix ¶

frf_matrix(frf: Any, records: Any = None, values: Any = None) -> tuple[list[tuple[str, str]], list[tuple[str, str]], ndarray, ndarray]

The records of an FRF assembled into the response-by-reference matrix they belong to, at every line.

A row is a response channel and a column a reference channel, each named (DOF, quantity) from the record's dimension, so a drive point's accelerometer and force gauge are not one channel. records picks the records (all of them when omitted) and values, when given, are their values in that order (the display units, for a plot); otherwise the stored SI values are used. Returns (rows, columns, matrix of lines x rows x columns, filled), where filled says which entries a record supplied; an entry no record supplied is zero in the matrix. The CMIF and the conditioning reading both assemble their matrix here.

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

def frf_matrix(frf: Any, records: Any = None, values: Any = None
               ) -> tuple[list[tuple[str, str]], list[tuple[str, str]],
                          np.ndarray, np.ndarray]:
    """The records of an FRF assembled into the response-by-reference
    matrix they belong to, at every line.

    A row is a response channel and a column a reference channel, each
    named (DOF, quantity) from the record's dimension, so a drive
    point's accelerometer and force gauge are not one channel. `records`
    picks the records (all of them when omitted) and `values`, when
    given, are their values in that order (the display units, for a
    plot); otherwise the stored SI values are used. Returns (rows,
    columns, matrix of lines x rows x columns, filled), where `filled`
    says which entries a record supplied; an entry no record supplied
    is zero in the matrix. The CMIF and the conditioning reading both
    assemble their matrix here.
    """
    from .data import channel_quantities

    wanted = (list(range(frf.num_records)) if records is None
              else [int(i) for i in records])
    values = (np.asarray(frf.ordinate)[wanted] if values is None
              else np.asarray(values))
    rows: list[tuple[str, str]] = []
    columns: list[tuple[str, str]] = []
    keys = []
    for i in wanted:
        out, into = channel_quantities(frf.known_dim(i))
        row = (str(frf.response_dof[i]), out)
        column = (str(frf.reference_dof[i]), into)
        if row not in rows:
            rows.append(row)
        if column not in columns:
            columns.append(column)
        keys.append((rows.index(row), columns.index(column)))
    matrix = np.zeros((len(frf.abscissa), len(rows), len(columns)),
                      dtype=np.complex128)
    filled = np.zeros((len(rows), len(columns)), dtype=bool)
    for position, (r, c) in enumerate(keys):
        matrix[:, r, c] = values[position]
        filled[r, c] = True
    return rows, columns, matrix, filled

entry_channels ¶

entry_channels(data: Any, record: int) -> tuple[tuple[str, str], tuple[str, str]]

(row channel, column channel) of one record of a spectral matrix.

A channel is (DOF, quantity), read from the record's dimension (channel_quantities): a cross term 'acceleration*force/frequency' runs from an acceleration channel to a force channel, so a drive point's accelerometer and its force gauge are two channels. A record with no reference is its own channel's autospectrum.

Parameters:

Name	Type	Description	Default
`data`	`Psd or Specification`	The matrix's records.	required
`record`	`int`	Which record.	required

Returns:

Type	Description
`tuple of ((str, str), (str, str))`	The row channel and the column channel.

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

def entry_channels(data: Any, record: int
                   ) -> tuple[tuple[str, str], tuple[str, str]]:
    """(row channel, column channel) of one record of a spectral matrix.

    A channel is (DOF, quantity), read from the record's dimension
    (`channel_quantities`): a cross term 'acceleration*force/frequency'
    runs from an acceleration channel to a force channel, so a drive
    point's accelerometer and its force gauge are two channels. A
    record with no reference is its own channel's autospectrum.

    Parameters
    ----------
    data : Psd or Specification
        The matrix's records.
    record : int
        Which record.

    Returns
    -------
    tuple of ((str, str), (str, str))
        The row channel and the column channel.
    """
    from .data import channel_quantities

    response, reference = data.record_pair(record)
    out, into = channel_quantities(data.known_dim(record))
    return (str(response), out), (str(reference), into)

matrix_entries ¶

matrix_entries(data: Any, records: Any = None, quantity: str | None = None) -> list[tuple[int, tuple[str, str], tuple[str, str]]]

Each record of a spectral matrix with the entry it holds.

Parameters:

Name	Type	Description	Default
`data`	`Psd or Specification`	The matrix's records.	required
`records`	`sequence of int`	Which records; all of them when omitted.	`None`
`quantity`	`str`	Keep only the records whose two channels are both of this quantity.	`None`

Returns:

Type	Description
`list of (int, (str, str), (str, str))`	The record, its row channel and its column channel (`entry_channels`), in record order.

Raises:

Type	Description
`ValueError`	If two records hold one entry: which of them is the matrix is not in the data.

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

def matrix_entries(data: Any, records: Any = None,
                   quantity: str | None = None
                   ) -> list[tuple[int, tuple[str, str], tuple[str, str]]]:
    """Each record of a spectral matrix with the entry it holds.

    Parameters
    ----------
    data : Psd or Specification
        The matrix's records.
    records : sequence of int, optional
        Which records; all of them when omitted.
    quantity : str, optional
        Keep only the records whose two channels are both of this
        quantity.

    Returns
    -------
    list of (int, (str, str), (str, str))
        The record, its row channel and its column channel
        (`entry_channels`), in record order.

    Raises
    ------
    ValueError
        If two records hold one entry: which of them is the matrix is
        not in the data.
    """
    wanted = (range(data.num_records) if records is None
              else [int(i) for i in records])
    entries = []
    held: dict[tuple[tuple[str, str], tuple[str, str]], int] = {}
    for i in wanted:
        row, column = entry_channels(data, i)
        if quantity is not None and (row[1], column[1]) != (quantity,
                                                             quantity):
            continue
        if (row, column) in held:
            raise ValueError(
                f'{data.record_label(held[(row, column)])} and '
                f'{data.record_label(i)} are both the entry of '
                f'{row[0]} against {column[0]}; keep one')
        held[(row, column)] = i
        entries.append((i, row, column))
    return entries

mirrors_disagree ¶

mirrors_disagree(term: Any, mirror: Any) -> bool

Whether a cross term and its mirror, read at the same lines, are not each other's conjugate, to a part in a million of the larger, on the lines where both are numbers. A line that is not a number is left to the note that says so.

Parameters:

Name	Type	Description	Default
`term`	`array_like`	The two halves of one pair, complex, one value per line.	required
`mirror`	`array_like`	The two halves of one pair, complex, one value per line.	required

Returns:

Type	Description
`bool`	True when they disagree.

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

def mirrors_disagree(term: Any, mirror: Any) -> bool:
    """Whether a cross term and its mirror, read at the same lines, are
    not each other's conjugate, to a part in a million of the larger, on
    the lines where both are numbers. A line that is not a number is
    left to the note that says so.

    Parameters
    ----------
    term, mirror : array_like
        The two halves of one pair, complex, one value per line.

    Returns
    -------
    bool
        True when they disagree.
    """
    term = np.asarray(term, dtype=complex)
    other = np.conj(np.asarray(mirror, dtype=complex))
    finite = np.isfinite(term) & np.isfinite(other)
    if not finite.any():
        return False
    term, other = term[finite], other[finite]
    scale = float(np.max(np.abs(np.concatenate([term, other]))))
    return not np.allclose(term, other, rtol=1e-6, atol=1e-9 * scale)

spectral_matrix ¶

spectral_matrix(psd: Any, records: Any = None, quantity: str | None = None) -> tuple[list[tuple[str, str]], ndarray, ndarray]

The records of a CSD assembled into its square matrix, at every line.

A channel is (DOF, quantity), read from each record's dimension: a cross term 'acceleration*force/frequency' runs from an acceleration channel to a force channel. records picks the records (all when omitted) and quantity keeps only the records whose two channels are both of it. Returns (channels, matrix of lines x n x n, filled), with an entry no record supplied zero and not filled; two records for one entry are refused (matrix_entries). The principal shapes and the realizability and conditioning readings all assemble their matrix here.

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

def spectral_matrix(psd: Any, records: Any = None,
                    quantity: str | None = None
                    ) -> tuple[list[tuple[str, str]], np.ndarray, np.ndarray]:
    """The records of a CSD assembled into its square matrix, at every
    line.

    A channel is (DOF, quantity), read from each record's dimension: a
    cross term 'acceleration*force/frequency' runs from an acceleration
    channel to a force channel. `records` picks the records (all when
    omitted) and `quantity` keeps only the records whose two channels
    are both of it. Returns (channels, matrix of lines x n x n, filled),
    with an entry no record supplied zero and not filled; two records
    for one entry are refused (`matrix_entries`). The principal shapes
    and the realizability and conditioning readings all assemble their
    matrix here.
    """
    entries = matrix_entries(psd, records, quantity)
    channels: list[tuple[str, str]] = []
    for _i, row, _column in entries:
        if row not in channels:
            channels.append(row)
    for _i, _row, column in entries:
        if column not in channels:
            channels.append(column)
    index = {key: k for k, key in enumerate(channels)}
    n = len(channels)
    ordinate = np.asarray(psd.ordinate)
    matrix = np.zeros((len(psd.abscissa), n, n), dtype=np.complex128)
    filled = np.zeros((n, n), dtype=bool)
    for i, row, column in entries:
        matrix[:, index[row], index[column]] = ordinate[i]
        filled[index[row], index[column]] = True
    return channels, matrix, filled

flagged_runs ¶

flagged_runs(frequency: ndarray, flags: ndarray) -> list[tuple[float, float]]

(first, last) frequency of each run of consecutive flagged lines.

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

def flagged_runs(frequency: np.ndarray, flags: np.ndarray
                 ) -> list[tuple[float, float]]:
    """(first, last) frequency of each run of consecutive flagged
    lines."""
    return [(float(frequency[a]), float(frequency[b])) for a, b in runs(flags)]

condition_matrix ¶

condition_matrix(frequency: Any, matrix: Any, rows: Any, columns: Any, *, limit: float = CONDITION_LIMIT, dimension: str = UNKNOWN, kind: str = 'FRF', notes: Any = (), progress: Callable[[int, int], None] | None = None) -> Conditioning

The singular values of a matrix at every line, as a Conditioning.

matrix is lines x rows x columns, complex, one unit throughout. A line holding any value that is not a finite number is not judged (its singular values are not a number) and is noted with the first entry that holds one; so is a line where the matrix is all zero, and a column (or, for an FRF, a row) that is zero at every line is named. limit is the condition number lines are flagged over, more than one. progress, when given, is told (0, core.progress.BUSY) before the decomposition, which cannot count itself.

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

def condition_matrix(frequency: Any, matrix: Any, rows: Any, columns: Any,
                     *, limit: float = CONDITION_LIMIT,
                     dimension: str = UNKNOWN, kind: str = 'FRF',
                     notes: Any = (),
                     progress: Callable[[int, int], None] | None = None
                     ) -> Conditioning:
    """The singular values of a matrix at every line, as a
    `Conditioning`.

    `matrix` is lines x rows x columns, complex, one unit throughout.
    A line holding any value that is not a finite number is not judged
    (its singular values are not a number) and is noted with the first
    entry that holds one; so is a line where the matrix is all zero, and
    a column (or, for an FRF, a row) that is zero at every line is named.
    `limit` is the condition number lines are flagged over, more than
    one. `progress`, when given, is told ``(0, core.progress.BUSY)``
    before the decomposition, which cannot count itself.
    """
    limit = float(limit)
    if not limit > 1.0:
        raise ValueError('a condition-number limit is more than 1; every '
                         'matrix has a condition number of at least 1')
    f = np.asarray(frequency, dtype=np.float64)
    matrix = np.asarray(matrix, dtype=np.complex128)
    if matrix.ndim != 3 or matrix.shape[0] != len(f):
        raise ValueError('the matrix is one rows-by-columns matrix per '
                         'line')
    size = min(matrix.shape[1], matrix.shape[2])
    if size == 0:
        raise ValueError('the matrix is empty; there is nothing to read')
    finite = np.isfinite(matrix).all(axis=(1, 2))
    singular = np.full((size, len(f)), np.nan)
    if finite.any():
        Ticker(progress).busy()
        singular[:, finite] = np.linalg.svd(matrix[finite],
                                            compute_uv=False).T
    notes = list(notes)
    if not finite.all():
        notes.append(_not_finite_note(matrix, finite, rows, columns))
    zero = finite & ~(singular[0] > 0.0)
    if zero.any():
        notes.append(f'{int(zero.sum())} lines are zero throughout and '
                     'are not judged')
    dead = [str(columns[c]) for c in range(matrix.shape[2])
            if finite.any() and not np.any(matrix[finite][:, :, c])]
    if dead and finite.any():
        notes.append(f'{", ".join(dead)} reach{"es" * (len(dead) == 1)} '
                     'no row at any line: the matrix is singular wherever '
                     'it is judged')
    # a response that nothing reaches is a dead channel too; in a CSD it
    # is its own column, named above
    silent = [str(rows[r]) for r in range(matrix.shape[1])
              if kind == 'FRF' and finite.any()
              and not np.any(matrix[finite][:, r, :])]
    if silent:
        notes.append(f'{", ".join(silent)} '
                     f'{"is" if len(silent) == 1 else "are"} reached by no '
                     'reference at any line: a dead channel gives the '
                     'matrix a zero row')
    if kind == 'FRF' and matrix.shape[2] > matrix.shape[1]:
        notes.append(f'{matrix.shape[2]} references drive '
                     f'{matrix.shape[1]} responses: the references can '
                     'make one motion in more than one way, and the '
                     f'condition number is of the {size} gains there are')
    return Conditioning(frequency=f, singular=singular,
                        rows=tuple(str(r) for r in rows),
                        columns=tuple(str(c) for c in columns),
                        limit=limit, dimension=str(dimension), kind=kind,
                        notes=tuple(notes))

read_conditioning ¶

read_conditioning(data: Any, records: Any = None, *, limit: float = CONDITION_LIMIT, progress: Callable[[int, int], None] | None = None) -> Conditioning

The singular values and condition number of an FRF matrix, or of a cross-spectral matrix, at every line.

An Frf gives the plant matrix, responses by references; it needs two references or more, since a condition number compares drives. A Psd with cross terms gives the square cross-spectral matrix of its channels, such as the drives' Gff. records picks the records (all of them when omitted): the matrix they make must be complete (a CSD term held only by its mirror is taken as the mirror's conjugate) and in one unit. limit is the condition number lines are flagged over. progress is passed to condition_matrix.

Raises TypeError for an object of neither kind, and ValueError for an empty selection, records in different units, an incomplete matrix, an FRF with one reference, or a CSD of one channel.

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

def read_conditioning(data: Any, records: Any = None, *,
                      limit: float = CONDITION_LIMIT,
                      progress: Callable[[int, int], None] | None = None
                      ) -> Conditioning:
    """The singular values and condition number of an FRF matrix, or of
    a cross-spectral matrix, at every line.

    An `Frf` gives the plant matrix, responses by references; it needs
    two references or more, since a condition number compares drives.
    A `Psd` with cross terms gives the square cross-spectral matrix of
    its channels, such as the drives' `Gff`. `records` picks the records
    (all of them when omitted): the matrix they make must be complete
    (a CSD term held only by its mirror is taken as the mirror's
    conjugate) and in one unit. `limit` is the condition number lines
    are flagged over. `progress` is passed to `condition_matrix`.

    Raises TypeError for an object of neither kind, and ValueError for
    an empty selection, records in different units, an incomplete
    matrix, an FRF with one reference, or a CSD of one channel.
    """
    from .data import Frf, Psd

    wanted = (list(range(data.num_records)) if records is None
              else [int(i) for i in records])
    if not wanted:
        raise ValueError('no records are selected; there is no matrix to '
                         'read')
    if isinstance(data, Frf):
        if data.reference_dof is None:
            raise ValueError('the FRF names no references')
        dimension, notes = _one_unit(data, wanted)
        rows, columns, matrix, filled = frf_matrix(data, wanted)
        if len(columns) < 2:
            raise ValueError(
                'conditioning compares references, and these records have '
                f'one ({channel_labels(columns)[0]}); select an FRF from two '
                'references or more')
        if not filled.all():
            r, c = np.argwhere(~filled)[0]
            raise ValueError(
                f'the matrix is incomplete: {channel_labels(rows)[r]} has no FRF '
                f'to {channel_labels(columns)[c]}. Select every response against '
                'every reference')
        return condition_matrix(data.abscissa, matrix, channel_labels(rows),
                                channel_labels(columns), limit=limit,
                                dimension=dimension, kind='FRF',
                                notes=notes, progress=progress)
    if isinstance(data, Psd):
        if data.reference_dof is None:
            raise ValueError('an autospectrum set has no cross terms; the '
                             'matrix needs them')
        channels, matrix, filled = spectral_matrix(data, wanted)
        quantities = list(dict.fromkeys(q for _dof, q in channels))
        if len(quantities) > 1:
            raise ValueError(
                f'the channels measure different quantities '
                f'({" and ".join(quantities[:2])}): the singular values of '
                'a matrix in different units change when one unit does. '
                'Select the records of one quantity')
        if len(channels) < 2:
            raise ValueError('conditioning compares channels, and these '
                             'records have one')
        labels = channel_labels(channels)
        matrix, notes = _completed(matrix, filled, labels)
        dimension, unit_notes = _one_unit(data, wanted)
        return condition_matrix(data.abscissa, matrix, labels, labels,
                                limit=limit, dimension=dimension,
                                kind='CSD', notes=unit_notes + notes,
                                progress=progress)
    raise TypeError(f'a {type(data).__name__} is neither an FRF nor a CSD; '
                    'there is no matrix to condition')

repair_coherence ¶

repair_coherence(matrix: Any, keep_diagonal: bool = True, progress: Callable[[int, int], None] | None = None) -> ndarray

The nearest realizable matrix to each Hermitian coherence matrix, in the Frobenius norm.

matrix is lines x n x n. With keep_diagonal, the diagonal is held where it is (zero stays zero, a negative entry goes to zero) and the nearest positive semi-definite matrix with that diagonal is found by alternating projections with Dykstra's correction (Higham, 2002), then scaled onto the diagonal exactly, which keeps it positive semi-definite. A line still moving after REPAIR_ITERATIONS passes is realizable all the same, if not quite the nearest. Without it, the negative eigenvalues are set to zero (Higham, 1988) and the diagonal moves too. progress, when given, is told (done, total) pass by pass, out of REPAIR_ITERATIONS; whatever it raises stops the repair.

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

def repair_coherence(matrix: Any, keep_diagonal: bool = True,
                     progress: Callable[[int, int], None] | None = None
                     ) -> np.ndarray:
    """The nearest realizable matrix to each Hermitian coherence matrix,
    in the Frobenius norm.

    `matrix` is lines x n x n. With `keep_diagonal`, the diagonal is
    held where it is (zero stays zero, a negative entry goes to zero)
    and the nearest positive semi-definite matrix with that diagonal is
    found by alternating projections with Dykstra's correction (Higham,
    2002), then scaled onto the diagonal exactly, which keeps it positive
    semi-definite. A line still moving after `REPAIR_ITERATIONS` passes
    is realizable all the same, if not quite the nearest. Without it,
    the negative eigenvalues are set to zero (Higham, 1988) and the
    diagonal moves too. `progress`, when given, is told ``(done,
    total)`` pass by pass, out of `REPAIR_ITERATIONS`; whatever it
    raises stops the repair.
    """
    matrix = np.asarray(matrix, dtype=np.complex128)
    hermitian = 0.5 * (matrix + np.conj(np.swapaxes(matrix, -2, -1)))
    if not keep_diagonal:
        return _psd_part(hermitian)
    target = np.maximum(np.real(np.diagonal(hermitian, axis1=1, axis2=2)),
                        0.0)
    n = hermitian.shape[-1]
    at = np.arange(n)
    y = hermitian.copy()
    y[:, at, at] = target
    correction = np.zeros_like(y)
    x = y.copy()
    # each line stops when it has converged, so a target of thousands
    # of lines pays for its slowest lines alone
    active = np.arange(len(y))
    ticker = Ticker(progress)
    ticker.add(REPAIR_ITERATIONS)
    for _ in range(REPAIR_ITERATIONS):
        if not len(active):
            break
        ticker.tick()
        r = y[active] - correction[active]
        projected = _psd_part(r)
        correction[active] = projected - r
        x[active] = projected
        moved = projected.copy()
        moved[:, at, at] = target[active]
        change = np.max(np.abs(moved - y[active]), axis=(1, 2))
        y[active] = moved
        active = active[change >= REPAIR_CONVERGED]
    ticker.tick(ticker.total - ticker.done)
    # x is positive semi-definite and its diagonal is within the
    # convergence of the target: a scaling by a positive diagonal puts
    # it on the target exactly and keeps it positive semi-definite. A
    # channel held at zero has its row and column zero
    diagonal = np.real(np.diagonal(x, axis1=1, axis2=2))
    with np.errstate(divide='ignore', invalid='ignore'):
        factor = np.where((target > 0.0) & (diagonal > 0.0),
                          np.sqrt(target / diagonal), 0.0)
    out = x * factor[:, :, None] * factor[:, None, :]
    out[:, at, at] = target
    return out

judge_lines ¶

judge_lines(frequency: Any, matrix: Any, channels: Any, *, tolerance: float = REALIZABLE_TOLERANCE, keep_diagonal: bool = True) -> Realizability

Whether a cross-spectral matrix is realizable at every line, with nothing repaired and nothing noted: check_spectral_matrix without its repair. repaired is the coherence matrix as it stands.

Parameters:

Name	Type	Description	Default
`frequency`	`array_like`	The lines, in Hz.	required
`matrix`	`array_like`	Lines x n x n, complex, any units.	required
`channels`	`sequence`	A label per channel.	required
`tolerance`	`float`	As for `check_spectral_matrix`.	`REALIZABLE_TOLERANCE`
`keep_diagonal`	`bool`	Carried into the reading.	`True`

Returns:

Type	Description
`Realizability`	Its `failing` says which lines fail.

Raises:

Type	Description
`ValueError`	If the tolerance is not between 0 and 1, or the matrix is not one square matrix per line, or is empty.

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

def judge_lines(frequency: Any, matrix: Any, channels: Any, *,
                tolerance: float = REALIZABLE_TOLERANCE,
                keep_diagonal: bool = True) -> Realizability:
    """Whether a cross-spectral matrix is realizable at every line, with
    nothing repaired and nothing noted: `check_spectral_matrix` without
    its repair. `repaired` is the coherence matrix as it stands.

    Parameters
    ----------
    frequency : array_like
        The lines, in Hz.
    matrix : array_like
        Lines x n x n, complex, any units.
    channels : sequence
        A label per channel.
    tolerance : float, optional
        As for `check_spectral_matrix`.
    keep_diagonal : bool, optional
        Carried into the reading.

    Returns
    -------
    Realizability
        Its `failing` says which lines fail.

    Raises
    ------
    ValueError
        If the tolerance is not between 0 and 1, or the matrix is not
        one square matrix per line, or is empty.
    """
    tolerance = float(tolerance)
    if not 0.0 < tolerance < 1.0:
        raise ValueError('a realizability tolerance lies between 0 and 1')
    f = np.asarray(frequency, dtype=np.float64)
    matrix = np.asarray(matrix, dtype=np.complex128)
    if (matrix.ndim != 3 or matrix.shape[1] != matrix.shape[2]
            or matrix.shape[0] != len(f)):
        raise ValueError('the cross-spectral matrix is one square matrix '
                         'per line')
    n = matrix.shape[1]
    if n == 0:
        raise ValueError('the matrix is empty; there is nothing to judge')
    finite = np.isfinite(matrix).all(axis=(1, 2))
    safe = np.where(finite[:, None, None], matrix, 0.0)
    # a line of zeros is judged by nothing: its eigenvalues are zero and
    # would read as realizable, a pass on no evidence
    zero = finite & ~np.any(safe, axis=(1, 2))
    judged = finite & ~zero
    coherence, scale = _normalized(safe)
    conjugate = np.conj(np.swapaxes(coherence, -2, -1))
    asymmetry = 0.5 * np.max(np.abs(coherence - conjugate), axis=(1, 2))
    diagonal = np.real(np.diagonal(safe, axis1=1, axis2=2))
    negative = finite & (diagonal < 0.0).any(axis=1)
    eigenvalues = np.full((n, len(f)), np.nan)
    if judged.any():
        eigenvalues[:, judged] = np.linalg.eigvalsh(
            0.5 * (coherence[judged] + conjugate[judged])).T
    coherence = np.where(finite[:, None, None], coherence, np.nan + 0j)
    return Realizability(
        frequency=f, coherence=coherence, scale=scale,
        eigenvalues=eigenvalues,
        asymmetry=np.where(finite, asymmetry, np.nan), negative=negative,
        repaired=coherence, channels=tuple(str(c) for c in channels),
        tolerance=tolerance, keep_diagonal=bool(keep_diagonal), notes=())

check_spectral_matrix ¶

check_spectral_matrix(frequency: Any, matrix: Any, channels: Any, *, tolerance: float = REALIZABLE_TOLERANCE, keep_diagonal: bool = True, notes: Any = (), progress: Callable[[int, int], None] | None = None) -> Realizability

Whether a cross-spectral matrix is realizable at every line, and the nearest that is, as a Realizability.

matrix is lines x n x n, complex, any units (the check is on the coherence matrix). A line holding a value that is not a finite number is not judged, and is noted with the first entry that holds one; so is a line where the matrix is zero throughout, since a matrix of nothing says nothing about how channels move together. A channel with no power and no cross term at any line judged is named. Lines that fail are repaired (repair_coherence, which progress is passed to); the rest are left exactly as they are.

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

def check_spectral_matrix(frequency: Any, matrix: Any, channels: Any, *,
                          tolerance: float = REALIZABLE_TOLERANCE,
                          keep_diagonal: bool = True,
                          notes: Any = (),
                          progress: Callable[[int, int], None] | None = None
                          ) -> Realizability:
    """Whether a cross-spectral matrix is realizable at every line, and
    the nearest that is, as a `Realizability`.

    `matrix` is lines x n x n, complex, any units (the check is on the
    coherence matrix). A line holding a value that is not a finite
    number is not judged, and is noted with the first entry that holds
    one; so is a line where the matrix is zero throughout, since a
    matrix of nothing says nothing about how channels move together. A
    channel with no power and no cross term at any line judged is named.
    Lines that fail are repaired (`repair_coherence`, which `progress`
    is passed to); the rest are left exactly as they are.
    """
    reading = judge_lines(frequency, matrix, channels, tolerance=tolerance,
                          keep_diagonal=keep_diagonal)
    f, coherence, scale = reading.frequency, reading.coherence, reading.scale
    eigenvalues, negative = reading.eigenvalues, reading.negative
    matrix = np.asarray(matrix, dtype=np.complex128)
    n = matrix.shape[1]
    finite = np.isfinite(matrix).all(axis=(1, 2))
    safe = np.where(finite[:, None, None], matrix, 0.0)
    zero = finite & ~np.any(safe, axis=(1, 2))
    judged = finite & ~zero
    asymmetry = reading.asymmetry
    failing = reading.failing
    repaired = coherence.copy()
    if failing.any():
        repaired[failing] = repair_coherence(coherence[failing],
                                             keep_diagonal, progress)
    notes = list(notes)
    labels = [str(c) for c in channels]
    if not finite.all():
        notes.append(_not_finite_note(matrix, finite, labels, labels))
    if zero.any():
        count = int(zero.sum())
        notes.append(f'{count} line{"s" * (count != 1)} '
                     f'{"are" if count != 1 else "is"} zero throughout and '
                     f'{"are" if count != 1 else "is"} not judged')
    if judged.any():
        live = np.any(safe[judged] != 0.0, axis=(0, 2))
        dead = [labels[k] for k in range(n) if not live[k]]
        if dead:
            notes.append(f'{", ".join(dead)} '
                         f'{"has" if len(dead) == 1 else "have"} no power '
                         'and no cross term at any line judged: a dead '
                         'channel is realizable, and says nothing about how '
                         'the others move with it')
    if negative.any():
        notes.append(f'{int(negative.sum())} lines have an autospectrum '
                     'below zero, which no signal has; the repair takes '
                     'it to zero')
    unsymmetric = finite & (asymmetry > tolerance)
    if unsymmetric.any():
        notes.append(f'{int(unsymmetric.sum())} lines are not Hermitian: '
                     'a cross term is not the conjugate of its mirror; '
                     'the repair takes their mean')
    if failing.any() and not keep_diagonal:
        notes.append('the repair was not held to the diagonal, so it '
                     'moves the autospectra too')
    return Realizability(
        frequency=f, coherence=coherence, scale=scale,
        eigenvalues=eigenvalues, asymmetry=reading.asymmetry,
        negative=negative, repaired=repaired, channels=reading.channels,
        tolerance=tolerance, keep_diagonal=bool(keep_diagonal),
        notes=tuple(notes))

read_realizability ¶

read_realizability(psd: Any, records: Any = None, *, tolerance: float = REALIZABLE_TOLERANCE, keep_diagonal: bool = True, progress: Callable[[int, int], None] | None = None) -> Realizability

Whether a CSD (a specification, or a measured cross-spectral matrix) is realizable at every line, and the nearest that is.

records picks the records (all of them when omitted); the matrix they make must be complete, a term held only by its mirror being taken as the mirror's conjugate. The channels may be of different quantities, since the check is on the coherence matrix, which has no units; but the records are declared in units throughout or not at all, since a coherence of a declared autospectrum against a cross term in whatever units a file held is not a coherence. A specification written at breakpoints is judged at its breakpoints. progress is passed to the repair (repair_coherence).

Raises TypeError for an object that is not a Psd, and ValueError for an empty selection, an object with no cross terms, records some declared in units and some not, fewer than two channels, or an incomplete matrix.

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

def read_realizability(psd: Any, records: Any = None, *,
                       tolerance: float = REALIZABLE_TOLERANCE,
                       keep_diagonal: bool = True,
                       progress: Callable[[int, int], None] | None = None
                       ) -> Realizability:
    """Whether a CSD (a specification, or a measured cross-spectral
    matrix) is realizable at every line, and the nearest that is.

    `records` picks the records (all of them when omitted); the matrix
    they make must be complete, a term held only by its mirror being
    taken as the mirror's conjugate. The channels may be of different
    quantities, since the check is on the coherence matrix, which has
    no units; but the records are declared in units throughout or not
    at all, since a coherence of a declared autospectrum against a cross
    term in whatever units a file held is not a coherence. A
    specification written at breakpoints is judged at its breakpoints.
    `progress` is passed to the repair (`repair_coherence`).

    Raises TypeError for an object that is not a `Psd`, and ValueError
    for an empty selection, an object with no cross terms, records some
    declared in units and some not, fewer than two channels, or an
    incomplete matrix.
    """
    from .data import Psd

    if not isinstance(psd, Psd):
        raise TypeError(f'a {type(psd).__name__} is not a cross-spectral '
                        'density; there is no matrix to judge')
    wanted = (list(range(psd.num_records)) if records is None
              else [int(i) for i in records])
    if not wanted:
        raise ValueError('no records are selected; there is no matrix to '
                         'judge')
    if psd.reference_dof is None:
        raise ValueError('an autospectrum set has no cross terms, and '
                         'realizability is a statement about them')
    undeclared = [i for i in wanted if str(psd.ordinate_dim[i]) == UNKNOWN]
    if undeclared and len(undeclared) < len(wanted):
        i = undeclared[0]
        raise ValueError(
            f'{psd.response_dof[i]} against {psd.reference_dof[i]} has no '
            'declared units and other records do: a coherence of terms in '
            'SI against terms in whatever units a file held is not a '
            'coherence. Declare the units of every record, or of none')
    channels, matrix, filled = spectral_matrix(psd, wanted)
    if len(channels) < 2:
        raise ValueError('realizability is about how channels move '
                         'together, and these records have one channel')
    labels = channel_labels(channels)
    matrix, notes = _completed(matrix, filled, labels)
    return check_spectral_matrix(psd.abscissa, matrix, labels,
                                 tolerance=tolerance,
                                 keep_diagonal=keep_diagonal, notes=notes,
                                 progress=progress)

made_realizable ¶

made_realizable(psd: Any, records: Any = None, *, tolerance: float = REALIZABLE_TOLERANCE, keep_diagonal: bool = True, progress: Callable[[int, int], None] | None = None) -> tuple[Any, Realizability]

A copy of a CSD with the matrix the picked records make repaired at every line that is not realizable, and the reading that says what moved.

Every record keeps its place, its DOFs and its units; only the values at failing lines change, the cross terms and, without keep_diagonal, the autospectra. A specification's limits follow an autospectrum that moved, so its bands stay where they were in decibels about it. progress is passed to the repair.

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

def made_realizable(psd: Any, records: Any = None, *,
                    tolerance: float = REALIZABLE_TOLERANCE,
                    keep_diagonal: bool = True,
                    progress: Callable[[int, int], None] | None = None
                    ) -> tuple[Any, Realizability]:
    """A copy of a CSD with the matrix the picked records make repaired
    at every line that is not realizable, and the reading that says what
    moved.

    Every record keeps its place, its DOFs and its units; only the
    values at failing lines change, the cross terms and, without
    `keep_diagonal`, the autospectra. A specification's limits follow an
    autospectrum that moved, so its bands stay where they were in
    decibels about it. `progress` is passed to the repair.
    """
    from .data import _recast, channel_quantities

    reading = read_realizability(psd, records, tolerance=tolerance,
                                 keep_diagonal=keep_diagonal,
                                 progress=progress)
    # every field copied, as reading a PSD as a specification does, so
    # nothing passes through a constructor that might drop one
    out = _recast(psd, psd.__class__)
    ordinate = np.array(psd.ordinate, dtype=np.complex128)
    failing = reading.failing
    if failing.any():
        before = reading.matrix(False)
        after = reading.matrix(True)
        wanted = (list(range(psd.num_records)) if records is None
                  else [int(i) for i in records])
        references = (psd.reference_dof if psd.reference_dof is not None
                      else psd.response_dof)
        # the reading's rows are the channels in `spectral_matrix`'s order
        channels, _matrix, _filled = spectral_matrix(psd, wanted)
        index = {key: k for k, key in enumerate(channels)}
        limits = getattr(out, 'limits', None)
        moved_limits = ({name: np.array(values, dtype=float)
                         for name, values in limits.items()}
                        if limits else None)
        for i in wanted:
            out_q, into_q = channel_quantities(psd.known_dim(i))
            r = index[(str(psd.response_dof[i]), out_q)]
            c = index[(str(references[i]), into_q)]
            values = after[failing, r, c]
            # an autospectrum is real; the repair's rounding is not
            ordinate[i, failing] = np.real(values) if r == c else values
            if r == c and moved_limits is not None:
                with np.errstate(divide='ignore', invalid='ignore'):
                    ratio = np.real(after[failing, r, r]) / np.real(
                        before[failing, r, r])
                ratio = np.where(np.isfinite(ratio), ratio, 1.0)
                for values in moved_limits.values():
                    values[i, failing] = values[i, failing] * ratio
        if moved_limits is not None:
            out.limits = moved_limits
    if not np.iscomplexobj(psd.ordinate):
        ordinate = np.real(ordinate)
    out.ordinate = ordinate
    return out, reading

entries_refusal ¶

entries_refusal(channels: Any, pairs: Any, frequencies: Any, tolerance: float = REALIZABLE_TOLERANCE) -> str | None

Why a specification sheet's coherence and phase entries cannot hold together, or None when they can.

pairs is {(i, j): (coherence, phase_degrees) or None} for i < j, as SpecificationDraft holds it. The entries do not change with frequency and a scaling by the autospectra does not change the verdict, so the coherence matrix they make is judged once, for every line from the first breakpoint to the last. A pair not entered is not judged: each largest group of channels whose pairs are all entered is judged on its own. The reason names the channels the failure is in, the eigenvalue, the band, and the nearest entries that would hold.

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

def entries_refusal(channels: Any, pairs: Any, frequencies: Any,
                    tolerance: float = REALIZABLE_TOLERANCE) -> str | None:
    """Why a specification sheet's coherence and phase entries cannot
    hold together, or None when they can.

    `pairs` is `{(i, j): (coherence, phase_degrees) or None}` for `i <
    j`, as `SpecificationDraft` holds it. The entries do not change with
    frequency and a scaling by the autospectra does not change the
    verdict, so the coherence matrix they make is judged once, for every
    line from the first breakpoint to the last. A pair not entered is
    not judged: each largest group of channels whose pairs are all
    entered is judged on its own. The reason names the channels the
    failure is in, the eigenvalue, the band, and the nearest entries
    that would hold.
    """
    channels = [str(c) for c in channels]
    stated = {(int(i), int(j)): value for (i, j), value in pairs.items()
              if value is not None}
    for group in _complete_groups(len(channels), set(stated)):
        if len(group) < 3:
            # two channels with a coherence in [0, 1] always hold
            continue
        m = len(group)
        c = np.eye(m, dtype=np.complex128)
        for a in range(m):
            for b in range(a + 1, m):
                coherence, phase = stated[(group[a], group[b])]
                c[a, b] = np.sqrt(coherence) * np.exp(1j * np.radians(phase))
                c[b, a] = np.conj(c[a, b])
        values, vectors = np.linalg.eigh(c)
        if not values[0] < -tolerance * values[-1]:
            continue
        weight = np.abs(vectors[:, 0]) ** 2
        named = [channels[group[k]] for k in range(m)
                 if weight[k] >= 0.5 / m]
        repaired = repair_coherence(c[None])[0]
        moves = []
        for a in range(m):
            for b in range(a + 1, m):
                moves.append((abs(repaired[a, b] - c[a, b]), a, b))
        # sizes equal but for rounding keep the channels' order
        moves.sort(key=lambda move: -round(float(move[0]), 9))
        shown = []
        for _size, a, b in moves[:3]:
            entry = repaired[a, b]
            # to a hundredth of a degree, so a phase the repair kept at
            # zero does not read as 1.7e-14
            phase = (round(float(np.degrees(np.angle(entry))), 2) + 0.0
                     if abs(entry) > 0.0 else 0.0)
            shown.append(f'{channels[group[a]]}–{channels[group[b]]} '
                         f'coherence {abs(entry) ** 2:.3g} at {phase:g}°')
        more = (f' and {len(moves) - 3} more' if len(moves) > 3 else '')
        who = (', '.join(named[:-1]) + ' and ' + named[-1]
               if len(named) > 1 else named[0])
        return (f'the cross terms of {who} cannot hold together: no set of '
                'signals has those coherences and phases (the matrix they '
                f'make has an eigenvalue of {values[0]:.3g} against a '
                f'largest of {values[-1]:.3g}), at every line from '
                f'{float(frequencies[0]):g} to {float(frequencies[-1]):g} '
                'Hz, since the entries do not change with frequency. The '
                f'nearest entries that hold are {"; ".join(shown)}{more}')
    return None

visualdynamics.core.matrix_checks¶

matrix_checks ¶

Classes¶

ConditioningPoint dataclass ¶

Attributes¶

ill property ¶

Methods:¶

describe ¶

Conditioning dataclass ¶

Attributes¶

size property ¶

judged property ¶

condition property ¶

rank property ¶

ill property ¶

title property ¶

Methods:¶

bands ¶

line ¶

at ¶

describe ¶

RealizabilityPoint dataclass ¶

Methods:¶

describe ¶

Realizability dataclass ¶

Attributes¶

judged property ¶

smallest property ¶

largest property ¶

failing property ¶

realizable property ¶

rank property ¶

Methods:¶

bands ¶

matrix ¶

pair ¶

worst ¶

changes ¶

line ¶

at ¶

describe ¶

Functions:¶

channel_labels ¶

runs ¶

nearest_line ¶

frf_matrix ¶

entry_channels ¶

matrix_entries ¶

mirrors_disagree ¶

spectral_matrix ¶

flagged_runs ¶

condition_matrix ¶

read_conditioning ¶

repair_coherence ¶

judge_lines ¶

check_spectral_matrix ¶

read_realizability ¶

made_realizable ¶

entries_refusal ¶

`visualdynamics.core.matrix_checks`¶

ConditioningPoint `dataclass` ¶

ill `property` ¶

Conditioning `dataclass` ¶

size `property` ¶

judged `property` ¶

condition `property` ¶

rank `property` ¶

ill `property` ¶

title `property` ¶

RealizabilityPoint `dataclass` ¶

Realizability `dataclass` ¶

judged `property` ¶

smallest `property` ¶

largest `property` ¶

failing `property` ¶

realizable `property` ¶

rank `property` ¶