r"""Modal Assurance Criterion and mode-to-mode matching.
The Modal Assurance Criterion (MAC) measures the correlation between two mode
shapes; a value of 1 indicates perfect correlation and 0 indicates orthogonality
[1]_. For complex (state-space) modes the pole-weighted MACXP criterion [2]_ is
also provided. ``match_modes`` performs a one-to-one Hungarian assignment over a MAC
(optionally frequency-penalized) matrix between two operating points; this is a
pairwise building block. Cross-operating-point *tracking* over a whole sweep is done
separately by :func:`vane.modal.identify_modes`, which extracts globally optimal
continuity paths rather than chaining pairwise matches.
References
----------
.. [1] Allemang, R. J. (2003). The modal assurance criterion — twenty years of
use and abuse. *Sound and Vibration*, 37(8), 14-23.
.. [2] Vacher, P., Jacquier, B., & Bucharles, A. (2010). Extensions of the MAC
criterion to complex modes. *Proceedings of ISMA 2010*, 2713-2726.
"""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
from scipy.optimize import linear_sum_assignment
if TYPE_CHECKING:
import numpy.typing as npt
__all__ = ["compute_mac", "compute_macxp", "match_modes"]
[docs]
def compute_mac(
phi_ref: npt.NDArray[np.complex128],
phi_test: npt.NDArray[np.complex128],
) -> npt.NDArray[np.float64]:
r"""Compute the Modal Assurance Criterion (MAC) matrix.
Parameters
----------
phi_ref : npt.NDArray[np.complex128]
Reference mode-shape matrix, shape ``(n_dof, n_modes_ref)``.
phi_test : npt.NDArray[np.complex128]
Test mode-shape matrix, shape ``(n_dof, n_modes_test)``.
Returns
-------
npt.NDArray[np.float64]
MAC matrix of shape ``(n_modes_ref, n_modes_test)`` with values in
``[0, 1]``.
Raises
------
ValueError
If the two mode-shape matrices have incompatible first dimensions.
Notes
-----
The MAC between modes :math:`i` and :math:`j` is
.. math::
\text{MAC}_{ij} = \frac{|\phi_i^H \phi_j|^2}
{(\phi_i^H \phi_i)(\phi_j^H \phi_j)}.
"""
if phi_ref.shape[0] != phi_test.shape[0]:
msg = (
f"Mode shapes have incompatible DOF dimensions: "
f"{phi_ref.shape[0]} vs {phi_test.shape[0]}"
)
raise ValueError(msg)
cross = np.abs(phi_ref.conj().T @ phi_test) ** 2
ref_norm = np.real(np.einsum("ij,ij->j", phi_ref.conj(), phi_ref))
test_norm = np.real(np.einsum("ij,ij->j", phi_test.conj(), phi_test))
denom = np.outer(ref_norm, test_norm)
with np.errstate(divide="ignore", invalid="ignore"):
mac = np.where(denom > 0.0, cross / denom, 0.0)
# Clip to [0, 1] to absorb floating-point overshoot on near-collinear modes.
return np.asarray(np.clip(mac, 0.0, 1.0), dtype=np.float64)
[docs]
def compute_macxp(
phi_ref: npt.NDArray[np.complex128],
lambda_ref: npt.NDArray[np.complex128],
phi_test: npt.NDArray[np.complex128],
lambda_test: npt.NDArray[np.complex128],
) -> npt.NDArray[np.float64]:
r"""Compute the pole-weighted MACXP matrix for complex modes.
MACXP augments the MAC with the transpose product and weights both by the
modal poles, making it more robust than MAC for complex state-space modes.
Parameters
----------
phi_ref, phi_test : npt.NDArray[np.complex128]
Reference and test mode-shape matrices, shape ``(n_dof, n_modes)``.
lambda_ref, lambda_test : npt.NDArray[np.complex128]
Corresponding modal poles (eigenvalues), shape ``(n_modes,)``.
Returns
-------
npt.NDArray[np.float64]
MACXP matrix of shape ``(n_modes_ref, n_modes_test)``; a mode correlated
with itself yields 1.
Raises
------
ValueError
If shapes are inconsistent between mode matrices and pole vectors.
Notes
-----
The criterion is defined for damped modes (poles with a non-zero real part);
an undamped pole makes the pole-weighting denominator zero, so such entries
are returned as NaN rather than raising.
"""
if phi_ref.shape[0] != phi_test.shape[0]:
msg = (
f"Mode shapes have incompatible DOF dimensions: "
f"{phi_ref.shape[0]} vs {phi_test.shape[0]}"
)
raise ValueError(msg)
if (
phi_ref.shape[1] != lambda_ref.shape[0]
or phi_test.shape[1] != lambda_test.shape[0]
):
msg = "Number of poles must match the number of mode shapes"
raise ValueError(msg)
hermitian = np.abs(phi_ref.conj().T @ phi_test)
transpose = np.abs(phi_ref.T @ phi_test)
l_ref = lambda_ref[:, np.newaxis]
l_test = lambda_test[np.newaxis, :]
with np.errstate(divide="ignore", invalid="ignore"):
numerator = hermitian / np.abs(l_ref.conj() + l_test) + transpose / np.abs(
l_ref + l_test
)
denom = np.outer(
_macxp_self(phi_ref, lambda_ref), _macxp_self(phi_test, lambda_test)
)
result = numerator**2 / denom
return np.asarray(result, dtype=np.float64)
[docs]
def match_modes(
mac: npt.NDArray[np.float64],
ref_frequencies: npt.NDArray[np.float64] | None = None,
test_frequencies: npt.NDArray[np.float64] | None = None,
frequency_weight: float = 0.0,
) -> list[tuple[int, int]]:
"""Match reference modes to test modes by maximum-weight assignment.
A Hungarian (optimal) assignment maximizes the total correlation between the
two mode sets. When frequencies are supplied, the MAC of each pair is scaled
by ``1 - frequency_weight * |df| / span`` so that frequency proximity breaks
ties between similar shapes (the strategy used for Campbell-diagram tracking).
Parameters
----------
mac : npt.NDArray[np.float64]
Correlation matrix, shape ``(n_ref, n_test)`` (e.g. from
:func:`compute_mac`).
ref_frequencies, test_frequencies : npt.NDArray[np.float64] or None, optional
Per-mode natural frequencies for the reference and test sets. Both must be
given to apply the frequency penalty.
frequency_weight : float, optional
Strength of the frequency penalty in ``[0, 1]``; ``0`` (default) uses MAC
alone.
Returns
-------
list[tuple[int, int]]
``(ref_index, test_index)`` pairs, one per matched mode (``min(n_ref,
n_test)`` pairs), ordered by reference index.
"""
weight = np.array(mac, dtype=np.float64, copy=True)
if (
ref_frequencies is not None
and test_frequencies is not None
and frequency_weight > 0.0
):
both = np.concatenate([ref_frequencies, test_frequencies])
span = float(both.max() - both.min())
if span > 0.0:
delta = np.abs(
ref_frequencies[:, np.newaxis] - test_frequencies[np.newaxis, :]
)
penalty = 1.0 - frequency_weight * delta / span
weight = weight * np.clip(penalty, 0.0, 1.0)
row_indices, col_indices = linear_sum_assignment(weight, maximize=True)
pairs = sorted(zip(row_indices.tolist(), col_indices.tolist(), strict=True))
return [(int(r), int(c)) for r, c in pairs]
def _macxp_self(
phi: npt.NDArray[np.complex128], lam: npt.NDArray[np.complex128]
) -> npt.NDArray[np.float64]:
"""Return the MACXP self-correlation denominator term for each mode."""
hermitian = np.real(np.einsum("ij,ij->j", phi.conj(), phi))
transpose = np.abs(np.einsum("ij,ij->j", phi, phi))
term = hermitian / (2.0 * np.abs(lam.real)) + transpose / (2.0 * np.abs(lam))
return np.asarray(term, dtype=np.float64)