Projectors and modal decompositions

The projector interface provides a common fit / encode / decode workflow for modal reduction. The package includes POD and spectral POD, together with standalone decomposition functions in pod_utils.

romda.models.data_driven.autoencoders.Projector

Bases: ABC

Abstract base for all dimensionality-reduction building blocks, linear (POD, SPOD) or nonlinear (autoencoders).

Every projector shares the same sklearn-style interface: fit learns the representation from data, encode maps state space to the latent space, and decode maps back. reconstruct and score are provided as concrete methods built on top of encode/decode.

Attributes:

Name	Type	Description
N_latent	int	Size of the latent (bottleneck) space. For POD/SPOD this is the number of modes retained.
fitted	bool	Whether fit has been called.
Q_mean	ndarray	Temporal mean removed from the data during preprocessing, shape \((N_x, 1)\). Raises AttributeError if accessed before fit.

fit(X) abstractmethod

Learn the projection from data.

Parameters:

Name	Type	Description	Default
X	ndarray	Snapshot data, shape \((N_x, N_t)\).	required

Returns:

Type	Description
Projector	The fitted instance (self).

Source code in src/models/data_driven/autoencoders/__init__.py

@abstractmethod
def fit(self, X: np.ndarray) -> Projector:
    r"""
    Learn the projection from data.

    Parameters
    ----------
    X : np.ndarray
        Snapshot data, shape $(N_x, N_t)$.

    Returns
    -------
    Projector
        The fitted instance (``self``).
    """

encode(X) abstractmethod

Map snapshot data to the latent representation.

Parameters:

Name	Type	Description	Default
X	ndarray	Snapshot data, shape \((N_x, N_t)\).	required

Returns:

Type	Description
ndarray	Latent coefficients \(\mathbf{Z}\), shape \((N_\mathrm{latent}, N_t)\).

Source code in src/models/data_driven/autoencoders/__init__.py

@abstractmethod
def encode(self, X: np.ndarray) -> np.ndarray:
    r"""
    Map snapshot data to the latent representation.

    Parameters
    ----------
    X : np.ndarray
        Snapshot data, shape $(N_x, N_t)$.

    Returns
    -------
    np.ndarray
        Latent coefficients $\mathbf{Z}$, shape $(N_\mathrm{latent}, N_t)$.
    """

decode(Z) abstractmethod

Map latent coefficients back to state space.

Parameters:

Name	Type	Description	Default
Z	ndarray	Latent coefficients, shape \((N_\mathrm{latent}, N_t)\).	required

Returns:

Type	Description
ndarray	Reconstructed state \(\hat{\mathbf{X}}\), shape \((N_x, N_t)\).

Source code in src/models/data_driven/autoencoders/__init__.py

@abstractmethod
def decode(self, Z: np.ndarray) -> np.ndarray:
    r"""
    Map latent coefficients back to state space.

    Parameters
    ----------
    Z : np.ndarray
        Latent coefficients, shape $(N_\mathrm{latent}, N_t)$.

    Returns
    -------
    np.ndarray
        Reconstructed state $\hat{\mathbf{X}}$, shape $(N_x, N_t)$.
    """

reconstruct(X)

Full round-trip: encode then decode.

Source code in src/models/data_driven/autoencoders/__init__.py

def reconstruct(self, X: np.ndarray) -> np.ndarray:
    """Full round-trip: `encode` then `decode`."""
    return self.decode(self.encode(X))

score(X)

Mean squared reconstruction error in the flat, zero-mean space,

\[ \mathrm{MSE} = \frac{1}{N_x N_t}\, \lVert \mathbf{Q} - \hat{\mathbf{Q}} \rVert_F^2, \]

where \(\mathbf{Q} = \mathrm{preprocess}(X)\) and \(\hat{\mathbf{Q}} = \mathrm{decode}(\mathrm{encode}(X)) - \bar{\mathbf{Q}}\) (with \(\bar{\mathbf{Q}}\) the stored Q_mean).

Parameters:

Name	Type	Description	Default
X	ndarray	Snapshot data, shape \((N_x, N_t)\).	required

Returns:

Type	Description
float	Mean squared reconstruction error.

Source code in src/models/data_driven/autoencoders/__init__.py

def score(self, X: np.ndarray) -> float:
    r"""
    Mean squared reconstruction error in the flat, zero-mean space,

    $$
    \mathrm{MSE} = \frac{1}{N_x N_t}\,
    \lVert \mathbf{Q} - \hat{\mathbf{Q}} \rVert_F^2,
    $$

    where $\mathbf{Q} = \mathrm{preprocess}(X)$ and
    $\hat{\mathbf{Q}} = \mathrm{decode}(\mathrm{encode}(X)) - \bar{\mathbf{Q}}$
    (with $\bar{\mathbf{Q}}$ the stored `Q_mean`).

    Parameters
    ----------
    X : np.ndarray
        Snapshot data, shape $(N_x, N_t)$.

    Returns
    -------
    float
        Mean squared reconstruction error.
    """
    Q = self.preprocess_snapshot(X)
    Q_hat = self.decode(self.encode(X)) - self.Q_mean
    return float(np.mean((Q - Q_hat) ** 2))

grid_index_to_flat_rows(grid_idx)

Map raw-grid indices to rows of the masked flat representation (Psi / Q_mean rows).

Both the raw grid and the flat representation use variable-block ordering: raw index = var * Nx * Ny + g (g = flattened (x, y) position), flat row = var * N_fluid + fluid_pos (position of g among the fluid points).

Parameters:

Name	Type	Description	Default
grid_idx	array-like of int	Raw-grid indices (e.g., sensor locations). Must correspond to fluid points.	required

Returns:

Type	Description
np.ndarray of int	Row indices into Psi / Q_mean corresponding to the requested grid points.

Source code in src/models/data_driven/autoencoders/__init__.py

def grid_index_to_flat_rows(self, grid_idx) -> np.ndarray:
    """
    Map raw-grid indices to rows of the masked flat representation (Psi / Q_mean rows).

    Both the raw grid and the flat representation use variable-block ordering:
    raw index = var * Nx * Ny + g (g = flattened (x, y) position),
    flat row  = var * N_fluid + fluid_pos (position of g among the fluid points).

    Parameters
    ----------
    grid_idx : array-like of int
        Raw-grid indices (e.g., sensor locations). Must correspond to fluid points.

    Returns
    -------
    np.ndarray of int
        Row indices into Psi / Q_mean corresponding to the requested grid points.
    """
    assert self.grid_shape is not None, 'grid_shape must be set to map grid indices.'
    Nu, Nx, Ny = self.grid_shape
    grid_idx = np.asarray(grid_idx).ravel()

    var = grid_idx // (Nx * Ny)
    g = grid_idx % (Nx * Ny)

    if not np.all(self.fluid_mask_flat[g]):
        raise ValueError('Some requested grid points are not fluid points.')

    fluid_idx = np.flatnonzero(self.fluid_mask_flat)
    N_fluid = fluid_idx.size
    fluid_pos = np.searchsorted(fluid_idx, g)
    return var * N_fluid + fluid_pos

preprocess_snapshot(X, subtract_mean=True)

Build the zero-mean data matrix from raw snapshot fields, automatically detecting and removing NaN-masked solid-body points.

Parameters:

Name	Type	Description	Default
X	ndarray	Raw snapshot data, either a single field \((N_t, N_x, N_y)\) or a stack of fields \((N_u, N_t, N_x, N_y)\).	required
subtract_mean	bool	If True (default), subtract the temporal mean row-wise.	True

Returns:

Type	Description
ndarray	Zero-mean data matrix \(\mathbf{Q}\), shape \((N_\mathrm{fluid} \cdot n_\mathrm{fields}, N_t)\), ready for decomposition.

Source code in src/models/data_driven/autoencoders/__init__.py

def preprocess_snapshot(self, X: np.ndarray, subtract_mean=True):
    r"""
    Build the zero-mean data matrix from raw snapshot fields,
    automatically detecting and removing NaN-masked solid-body points.

    Parameters
    ----------
    X : np.ndarray
        Raw snapshot data, either a single field $(N_t, N_x, N_y)$ or a
        stack of fields $(N_u, N_t, N_x, N_y)$.
    subtract_mean : bool
        If True (default), subtract the temporal mean row-wise.

    Returns
    -------
    np.ndarray
        Zero-mean data matrix $\mathbf{Q}$, shape
        $(N_\mathrm{fluid} \cdot n_\mathrm{fields}, N_t)$, ready for decomposition.
    """

    if not self.fitted:
        if X.ndim == 2:
            # Already-flat data matrix (N_x, N_t): no grid/mask handling
            self.fluid_mask_flat = np.ones(X.shape[0], dtype=bool)
            if subtract_mean:
                self.Q_mean = X.mean(axis=1, keepdims=True)
            else:
                self.Q_mean = np.zeros_like(X[:, :1])
            Q = X - self.Q_mean
            self._TKE = 0.5 * float(np.sum(np.mean(Q**2, axis=1)))
            return Q

        # if the input is raw grid data, we need to detect the fluid points and flatten the data
        assert X.ndim == 4, f'Expected raw grid input with 4 dimensions, got {X.ndim}.'
        Nu, Nt, Nx, Ny = X.shape
        self.grid_shape = (Nu, Nx, Ny)

        ref = X[0]
        fluid_mask = ~np.isnan(ref[0])
        self.fluid_mask_flat  = fluid_mask.ravel()



        X_masked_flat = self._to_flat(X)                    # (N_fluid * n_fields, N_t)

        if subtract_mean:
            self.Q_mean = X_masked_flat.mean(axis=1, keepdims=True)
        else:
            self.Q_mean = np.zeros_like(X_masked_flat[:, :1])


        Q = X_masked_flat - self.Q_mean #shape (N_fluid * n_fields, N_t)
        # store the total kinetic energy for later use in relative error metrics
        self._TKE = 0.5 * float(np.sum(np.mean(Q**2, axis=1)))
        return Q

    elif X.shape[0] != self.Q_mean.shape[0]:

        # if the decomosition is already fitted, can expect 1 snapshot only
        assert X.ndim in (3, 4), f'Expected flat input with 2, 3 or 4 dimensions, got {X.ndim}.'
        if X.ndim == 3:
            X = X[:, np.newaxis]    # (n_fields, 1, Nx, Ny)

        #check grid
        Nu, _, Nx, Ny = X.shape
        grid_shape = (Nu, Nx, Ny)
        assert grid_shape == self.grid_shape, f'Expected grid shape {self.grid_shape}, got {grid_shape}.'

        X_masked_flat = self._to_flat(X)                    # (N_fluid * n_fields, N_t)
        return X_masked_flat - self.Q_mean #shape (N_fluid * n_fields, N_t)
    else:
        # already flat input, just check dimensions and remove mean
        assert X.ndim == 2, f'Expected flat input with 2 dimensions, got {X.ndim}.'

        return X - self.Q_mean #shape (N_fluid * n_fields, N_t)

romda.models.data_driven.autoencoders.POD(n_modes=20, method='randomized', n_iter=4, random_state=None, grid_shape=None, domain=None, **kwargs)

Bases: Projector

Snapshot POD.

Inherits the shared interface from Projector (fit/encode/decode/ reconstruct/score/N_latent) and adds linear-specific attributes, geometry helpers, and plotting utilities.

Two solvers are available via method:

• 'exact' — full eigendecomposition of the temporal correlation matrix \(\mathbf{C} = \mathbf{Q}^\mathrm{T}\mathbf{Q} / N_t\) (Sirovich 1987, snapshot method: see snapshot_pod). Exact but \(\mathcal{O}(N_t^3)\).
• 'randomized' (default) — randomized SVD of \(\mathbf{Q}\) (Halko, Martinsson & Tropp 2011: see snapshot_pod_randomized). Returns only the leading n_modes modes; fast and memory-efficient.

Either way, writing \(\mathbf{Q} = \mathbf{X} - \bar{\mathbf{Q}}\) for the zero-mean data matrix, fit(X) stores the orthonormal spatial modes \(\boldsymbol{\Psi}\) and singular values \(\boldsymbol{\Sigma}\) of \(\mathbf{Q}\), together with the temporal coefficients \(\boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q}\), so that

\[ \mathbf{X} \approx \boldsymbol{\Psi}\boldsymbol{\Phi} + \bar{\mathbf{Q}}, \]

with equality when all \(N_t\) modes are retained (n_modes >= N_t, method='exact').

Parameters:

Name	Type	Description	Default
n_modes	int	Number of modes retained (the latent-space size). Default 20.	20
method	str	'exact' or 'randomized'. Default 'randomized'.	'randomized'
n_iter	int	Power-iteration steps for the randomized solver. Default 4.	4
random_state	int	Seed for reproducibility of the randomized solver.	None
grid_shape	tuple	Grid shape (Nu, Nx, Ny) used to map flat vectors back to the grid.	None
domain	list	Physical domain [x0, x1, y0, y1] used by the plotting utilities.	None
**kwargs		Pre-set any instance attribute (e.g. pre-computed Sigma, Psi, Phi, Q_mean), or pass X=... to fit directly at construction.	{}

Attributes:

Name	Type	Description
Sigma	ndarray	Singular values (descending), shape \((N_\mathrm{latent},)\).
Psi	ndarray	Spatial modes with orthonormal columns, shape \((N_x, N_\mathrm{latent})\).
Phi	ndarray	Temporal coefficients, shape \((N_\mathrm{latent}, N_t)\).
Q_mean	ndarray	Temporal mean \(\bar{\mathbf{Q}}\), shape \((N_x, 1)\).

References

Sirovich (1987). Turbulence and the dynamics of coherent structures. Quart. Appl. Math., XLV(3), 561-590.

Halko, Martinsson & Tropp (2011). Finding structure with randomness. SIAM Review, 53(2), 217-288.

Source code in src/models/data_driven/autoencoders/pod.py

def __init__(self,
             n_modes:      int            = 20,
             method:       str            = 'randomized',
             n_iter:       int            = 4,
             random_state: int | None  = None,
             grid_shape:   tuple | None = None,
             domain:       list | None  = None,
             **kwargs):
    self.N_latent     = n_modes
    self.grid_shape   = grid_shape
    self.domain       = domain
    self.method       = method
    self.n_iter       = n_iter
    self.random_state = random_state
    for key, val in kwargs.items():
        if hasattr(type(self), key) or key in ('Sigma', 'Psi', 'Phi',
                                                'Q_mean', '_TKE',
                                                'indices_to_original_grid'):
            setattr(self, key, val)
    # infer latent size from pre-loaded Phi if provided
    if self._Phi is not None and self.N_latent == 20:
        self.N_latent = self._Phi.shape[0]
    # backward compat: auto-fit if raw data provided as kwarg 'X'
    if 'X' in kwargs and kwargs['X'] is not None:
        self.fit(kwargs['X'])

N_modes property writable

Backward-compatible alias for N_latent.

domain_mesh property

Meshgrid for the spatial domain.

Returns:

Name	Type	Description
X1	ndarray	First coordinate array, shape (Nx, Ny).
X2	ndarray	Second coordinate array, shape (Nx, Ny).

fit(X)

Fit the POD to data X.

Accepts a flat matrix \((N_x, N_t)\) or a raw grid array \((N_u, N_t, N_x, N_y)\) / \((N_t, N_x, N_y)\). For raw grid input the NaN solid-body mask is detected and stored automatically.

Parameters:

Name	Type	Description	Default
X	ndarray	Snapshot data, flat or raw grid (see above).	required

Returns:

Type	Description
POD	The fitted instance (self).

Source code in src/models/data_driven/autoencoders/pod.py

def fit(self, X: np.ndarray) -> POD:
    r"""
    Fit the POD to data X.

    Accepts a flat matrix $(N_x, N_t)$ or a raw grid array
    $(N_u, N_t, N_x, N_y)$ / $(N_t, N_x, N_y)$. For raw grid input the
    NaN solid-body mask is detected and stored automatically.

    Parameters
    ----------
    X : np.ndarray
        Snapshot data, flat or raw grid (see above).

    Returns
    -------
    POD
        The fitted instance (``self``).
    """

    Q = self.preprocess_snapshot(X)

    result = self._decompose(Q)
    self.Sigma = result[0]
    self.Psi   = result[1]
    self.Phi   = result[2]
    assert self._Sigma is not None and self._Psi is not None and self._Phi is not None, \
        "Decomposition must return Sigma, Psi, Phi."
    self.N_latent = self.Sigma.shape[0]
    self.fitted = True
    return self

encode(X)

Project X onto the spatial modes,

\[ \mathbf{Z} = \boldsymbol{\Psi}^\mathrm{T}(\mathbf{X} - \bar{\mathbf{Q}}). \]

Parameters:

Name	Type	Description	Default
X	ndarray	Snapshot data, shape \((N_x, N_t)\).	required

Returns:

Type	Description
ndarray	Latent (POD) coefficients \(\mathbf{Z}\), shape \((N_\mathrm{latent}, N_t)\).

Source code in src/models/data_driven/autoencoders/pod.py

def encode(self, X: np.ndarray) -> np.ndarray:
    r"""
    Project X onto the spatial modes,

    $$
    \mathbf{Z} = \boldsymbol{\Psi}^\mathrm{T}(\mathbf{X} - \bar{\mathbf{Q}}).
    $$

    Parameters
    ----------
    X : np.ndarray
        Snapshot data, shape $(N_x, N_t)$.

    Returns
    -------
    np.ndarray
        Latent (POD) coefficients $\mathbf{Z}$, shape $(N_\mathrm{latent}, N_t)$.
    """
    Q = self.preprocess_snapshot(X)
    return self.Psi.T @ Q

decode(Z, idx=None)

Reconstruct the state in the original space from latent coefficients,

\[ \hat{\mathbf{Q}} = \boldsymbol{\Psi}\mathbf{Z} + \bar{\mathbf{Q}}. \]

Parameters:

Name	Type	Description	Default
Z	ndarray	Latent coefficients, shape \((N_\mathrm{latent}, N_t)\).	required
idx	ndarray	Indices selecting a subset of rows of Psi and Q_mean (e.g. sensor locations) to reconstruct only those entries.	None

Returns:

Type	Description
ndarray	Reconstructed state, shape \((N_x, N_t)\), or (len(idx), N_t) if idx is given.

Source code in src/models/data_driven/autoencoders/pod.py

def decode(self, Z: np.ndarray, idx: np.ndarray | None = None) -> np.ndarray:
    r"""
    Reconstruct the state in the original space from latent coefficients,

    $$
    \hat{\mathbf{Q}} = \boldsymbol{\Psi}\mathbf{Z} + \bar{\mathbf{Q}}.
    $$

    Parameters
    ----------
    Z : np.ndarray
        Latent coefficients, shape $(N_\mathrm{latent}, N_t)$.
    idx : np.ndarray, optional
        Indices selecting a subset of rows of `Psi` and `Q_mean`
        (e.g. sensor locations) to reconstruct only those entries.

    Returns
    -------
    np.ndarray
        Reconstructed state, shape $(N_x, N_t)$, or ``(len(idx), N_t)``
        if `idx` is given.
    """

    if idx is not None:
        return self.Psi[idx, :] @ Z + self.Q_mean[idx, :]
    else:
        return self.Psi @ Z + self.Q_mean

reconstruct(X=None, n_modes=None, Phi=None)

Full round-trip: encode, decode, then map back to the physical grid (when a grid mask is available).

Parameters:

Name	Type	Description	Default
X	ndarray	Raw input, grid or flat. If None, the stored Phi is used.	None
n_modes	int	Retain only the first n_modes modes. Default: all fitted modes.	None
Phi	ndarray	Pre-computed latent coefficients, shape \((N_\mathrm{latent}, N_t)\); skips the encode step.	None

Returns:

Type	Description
ndarray	Reconstructed state, on the physical grid if grid_shape and a to_grid mapping are available, otherwise flat.

Source code in src/models/data_driven/autoencoders/pod.py

def reconstruct(self, X: np.ndarray | None = None,
                n_modes: int | None = None,
                Phi: np.ndarray | None = None) -> np.ndarray:
    r"""
    Full round-trip: encode, decode, then map back to the physical grid
    (when a grid mask is available).

    Parameters
    ----------
    X : np.ndarray, optional
        Raw input, grid or flat. If None, the stored `Phi` is used.
    n_modes : int, optional
        Retain only the first ``n_modes`` modes. Default: all fitted modes.
    Phi : np.ndarray, optional
        Pre-computed latent coefficients, shape $(N_\mathrm{latent}, N_t)$;
        skips the `encode` step.

    Returns
    -------
    np.ndarray
        Reconstructed state, on the physical grid if `grid_shape` and a
        ``to_grid`` mapping are available, otherwise flat.
    """
    if Phi is not None:
        Z = Phi
    elif X is not None:
        Z = self.encode(X)
    else:
        Z = self.Phi

    nm = n_modes if n_modes is not None else self.N_latent
    if nm < self.N_latent:
        X_hat = self.Psi[:, :nm] @ Z[:nm] + self.Q_mean
    else:
        X_hat = self.decode(Z)
    if getattr(self, 'to_grid', None) is not None and self.grid_shape is not None:
        return self._to_physical_grid(X_hat)
    return X_hat

energy_fraction()

Relative and cumulative energy per mode, from the eigenvalues \(\lambda_j = \Sigma_j^2\) of the temporal correlation matrix \(\mathbf{C}\):

\[ \mathrm{rel}_j = \frac{\lambda_j}{\sum_k \lambda_k}, \qquad \mathrm{cum}_j = \sum_{k \le j} \mathrm{rel}_k. \]

Returns:

Name	Type	Description
rel	ndarray	Relative energy per mode, shape \((N_\mathrm{latent},)\).
cum	ndarray	Cumulative energy fraction captured by the first \(j\) modes, shape \((N_\mathrm{latent},)\).

Source code in src/models/data_driven/autoencoders/pod.py

def energy_fraction(self):
    r"""
    Relative and cumulative energy per mode, from the eigenvalues
    $\lambda_j = \Sigma_j^2$ of the temporal correlation matrix
    $\mathbf{C}$:

    $$
    \mathrm{rel}_j = \frac{\lambda_j}{\sum_k \lambda_k}, \qquad
    \mathrm{cum}_j = \sum_{k \le j} \mathrm{rel}_k.
    $$

    Returns
    -------
    rel : np.ndarray
        Relative energy per mode, shape $(N_\mathrm{latent},)$.
    cum : np.ndarray
        Cumulative energy fraction captured by the first $j$ modes,
        shape $(N_\mathrm{latent},)$.
    """
    lam = self.Sigma ** 2  # eigenvalues of C = Q^T Q / N_t

    return lam / lam.sum(), np.cumsum(lam) / lam.sum()

truncate(n_modes)

Truncate to the first n_modes modes in-place.

Source code in src/models/data_driven/autoencoders/pod.py

def truncate(self, n_modes: int) -> POD:
    """Truncate to the first n_modes modes in-place."""
    if n_modes >= self.N_latent:
        print(f"Requested n_modes={n_modes} >= N_latent={self.N_latent}. No truncation applied.")
        return self
    self.Psi      = self.Psi[:, :n_modes]
    self.Phi      = self.Phi[:n_modes, :]
    self.Sigma    = self.Sigma[:n_modes]
    self.N_latent = n_modes
    return self

original_data_to_domain_of_interest(original_data)

Crop original-grid data to the fitted domain of interest.

Source code in src/models/data_driven/autoencoders/pod.py

def original_data_to_domain_of_interest(self, original_data: np.ndarray):
    """Crop original-grid data to the fitted domain of interest."""
    if self.indices_to_original_grid is None:
        return original_data
    original_data = original_data.copy()
    try:
        if original_data.ndim == 2:
            return original_data[self.indices_to_original_grid]
        return original_data[:, self.indices_to_original_grid[0],
                                 self.indices_to_original_grid[1]]
    except Exception:
        raise ValueError(
            'Pass original_data in shape [(Nu) x Nx x Ny x (Nt)].')

compute_MSE(ROM_data, original_data, time_evolution=False) staticmethod

Mean Squared Error between ROM reconstruction and original data.

Source code in src/models/data_driven/autoencoders/pod.py

@staticmethod
def compute_MSE(ROM_data: np.ndarray, original_data: np.ndarray,
                time_evolution: bool = False):
    """Mean Squared Error between ROM reconstruction and original data."""
    ROM_data, original_data = POD.flatten(ROM_data, original_data)
    original_data[np.isnan(original_data)] = 0.
    if time_evolution:
        return np.mean((original_data - ROM_data) ** 2, axis=0)
    return float(np.mean((original_data - ROM_data) ** 2))

compute_RMS(ROM_data, original_data) staticmethod

Root Mean Square error (field).

Source code in src/models/data_driven/autoencoders/pod.py

@staticmethod
def compute_RMS(ROM_data: np.ndarray, original_data: np.ndarray):
    """Root Mean Square error (field)."""
    original_data[np.isnan(original_data)] = 0.
    return np.sqrt((original_data - ROM_data) ** 2)

flatten(*args) staticmethod

Flatten multi-dimensional arrays to 2-D (space × time).

Source code in src/models/data_driven/autoencoders/pod.py

@staticmethod
def flatten(*args):
    """Flatten multi-dimensional arrays to 2-D (space × time)."""
    return [a.reshape(-1, a.shape[-1]) if a.ndim > 2 else a.copy()
            for a in args]

romda.models.data_driven.autoencoders.SPOD(Nf=0, filter_kind='gaussian', n_modes=20, grid_shape=None, domain=None, **kwargs)

Bases: POD

Spectral POD via the filtered correlation matrix (Sieber, Paschereit & Oberleithner, 2016).

Inherits the full POD interface (fit/encode/decode/reconstruct); only _decompose is overridden, delegating to spod_sieber. The snapshot correlation matrix \(\mathbf{C} = \mathbf{Q}^\mathrm{T}\mathbf{Q}/N_t\) is replaced by a low-pass filtered version before the eigensolve,

\[ \tilde{\mathbf{C}} = \mathbf{G}^\mathrm{T} \mathbf{C}\, \mathbf{G}, \]

where \(\mathbf{G}\) is a banded symmetric Toeplitz filter matrix built from a normalised 1-D kernel of half-width Nf (filter_kind selects 'gaussian', 'box' or 'hann'). The eigendecomposition of \(\tilde{\mathbf{C}}\) then follows exactly as in snapshot_pod, giving orthonormal spatial modes \(\boldsymbol{\Psi}\), temporal coefficients \(\boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q}\) and singular values \(\boldsymbol{\Sigma}\). Setting Nf=0 skips the filtering step and exactly recovers snapshot POD; per Sieber et al. (2016), as Nf grows towards \(N_t/2\) the SPOD modes are reported to approach Fourier (DFT) modes.

Parameters:

Name	Type	Description	Default
Nf	int	Filter half-width (0 recovers POD; \(N_t/2\) approaches the DFT limit).	0
filter_kind	str	'gaussian', 'box' or 'hann'. Default 'gaussian'.	'gaussian'
n_modes	int	Number of modes to retain. Default 20.	20
grid_shape	tuple	Grid shape (Nu, Nx, Ny) for grid mapping.	None
domain	list	Physical domain [x0, x1, y0, y1] for the plotting utilities.	None
**kwargs		Forwarded to POD.__init__ (e.g. pre-set Sigma, Psi, Phi).	{}

Attributes:

Name	Type	Description
C_tilde	ndarray	Filtered correlation matrix \(\tilde{\mathbf{C}}\), shape \((N_t, N_t)\) (stored after fit).

References

Sieber, Paschereit & Oberleithner (2016). Spectral proper orthogonal decomposition. J. Fluid Mech., 792, 798–828.

Source code in src/models/data_driven/autoencoders/pod.py

def __init__(self,
             Nf:          int            = 0,
             filter_kind: str            = 'gaussian',
             n_modes:     int            = 20,
             grid_shape:  tuple | None = None,
             domain:      list | None  = None,
             **kwargs):
    super().__init__(n_modes=n_modes,
                     grid_shape=grid_shape,
                     domain=domain,
                     **kwargs)
    self.Nf                 = Nf
    self.filter_kind        = filter_kind
    self._n_modes_requested = n_modes   # None = keep all after fit

romda.models.data_driven.autoencoders.pod_utils.snapshot_pod(Q)

Snapshot POD — exact solver.

Solves the eigenvalue problem of the temporal correlation matrix,

\[ \mathbf{C}\mathbf{A} = \mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda}), \qquad \mathbf{C} = \mathbf{Q}^\mathrm{T} \mathbf{Q} / N_t, \qquad \lambda_1 \ge \lambda_2 \ge \cdots \ge 0, \]

and reconstructs the (large) spatial modes from the eigenvectors \(\mathbf{A}\) of the (small, \(N_t \times N_t\)) matrix \(\mathbf{C}\) — the "method of snapshots" of Sirovich (1987):

\[ \boldsymbol{\Psi} = \frac{1}{\sqrt{N_t}}\, \mathbf{Q}\mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda})^{-1/2}, \qquad \boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q}, \qquad \boldsymbol{\Sigma} = \sqrt{\boldsymbol{\lambda}}. \]

\(\boldsymbol{\Psi}\) has orthonormal columns and, for the modes with \(\lambda_j > 0\), \(\mathbf{Q} = \boldsymbol{\Psi}\boldsymbol{\Phi}\) exactly. Modes with \(\lambda_j \le 0\) (numerical noise) are set to zero.

Parameters:

Name	Type	Description	Default
Q	ndarray	Zero-mean data matrix, shape \((N_x, N_t)\).	required

Returns:

Name	Type	Description
Sigma	ndarray	Singular values (descending), shape \((N_t,)\).
Psi	ndarray	Spatial modes with orthonormal columns, shape \((N_x, N_t)\).
Phi	ndarray	Temporal coefficients, shape \((N_t, N_t)\).
C	ndarray	Temporal correlation matrix, shape \((N_t, N_t)\).

References

Sirovich (1987). Turbulence and the dynamics of coherent structures. Quart. Appl. Math., XLV(3), 561–590.

Source code in src/models/data_driven/autoencoders/pod_utils.py

def snapshot_pod(Q):
    r"""Snapshot POD — exact solver.

    Solves the eigenvalue problem of the temporal correlation matrix,

    $$
    \mathbf{C}\mathbf{A} = \mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda}),
    \qquad \mathbf{C} = \mathbf{Q}^\mathrm{T} \mathbf{Q} / N_t,
    \qquad \lambda_1 \ge \lambda_2 \ge \cdots \ge 0,
    $$

    and reconstructs the (large) spatial modes from the eigenvectors
    $\mathbf{A}$ of the (small, $N_t \times N_t$) matrix $\mathbf{C}$ — the
    "method of snapshots" of Sirovich (1987):

    $$
    \boldsymbol{\Psi} = \frac{1}{\sqrt{N_t}}\,
    \mathbf{Q}\mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda})^{-1/2},
    \qquad
    \boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q},
    \qquad
    \boldsymbol{\Sigma} = \sqrt{\boldsymbol{\lambda}}.
    $$

    $\boldsymbol{\Psi}$ has orthonormal columns and, for the modes with
    $\lambda_j > 0$, $\mathbf{Q} = \boldsymbol{\Psi}\boldsymbol{\Phi}$ exactly.
    Modes with $\lambda_j \le 0$ (numerical noise) are set to zero.

    Parameters
    ----------
    Q : np.ndarray
        Zero-mean data matrix, shape $(N_x, N_t)$.

    Returns
    -------
    Sigma : np.ndarray
        Singular values (descending), shape $(N_t,)$.
    Psi : np.ndarray
        Spatial modes with orthonormal columns, shape $(N_x, N_t)$.
    Phi : np.ndarray
        Temporal coefficients, shape $(N_t, N_t)$.
    C : np.ndarray
        Temporal correlation matrix, shape $(N_t, N_t)$.

    References
    ----------
    Sirovich (1987). Turbulence and the dynamics of coherent structures.
    *Quart. Appl. Math.*, XLV(3), 561–590.
    """
    _, N_t = Q.shape
    C      = (Q.T @ Q) / N_t
    lam, A = np.linalg.eigh(C)
    idx      = lam.argsort()[::-1]
    lam, A   = lam[idx], A[:, idx]
    safe_lam = np.where(lam > 0, lam, np.inf)
    Psi      = Q @ A / (np.sqrt(N_t) * np.sqrt(safe_lam))
    Psi[:, lam <= 0] = 0.0
    Phi      = Psi.T @ Q
    Sigma    = np.sqrt(np.where(lam > 0, lam, 0.0))
    return Sigma, Psi, Phi, C

romda.models.data_driven.autoencoders.pod_utils.snapshot_pod_randomized(Q, n_modes=20, n_iter=4, random_state=None)

Randomized snapshot POD.

Computes a truncated randomized SVD of \(\mathbf{Q}\) (via sklearn.utils.extmath.randomized_svd, falling back to a full numpy.linalg.svd if scikit-learn is unavailable),

\[ \mathbf{Q} \approx \mathbf{U}\,\mathbf{S}\,\mathbf{V}^\mathrm{T}, \]

and returns

\[ \boldsymbol{\Psi} = \mathbf{U}, \qquad \boldsymbol{\Sigma} = \mathbf{S} / \sqrt{N_t}, \qquad \boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q}. \]

The \(\boldsymbol{\Sigma} = \mathbf{S}/\sqrt{N_t}\) scaling matches the eigenvalue-based normalisation of snapshot_pod, since \(\mathbf{Q}^\mathrm{T}\mathbf{Q}/N_t = \mathbf{V}(\mathbf{S}^2/N_t)\mathbf{V}^\mathrm{T}\).

Parameters:

Name	Type	Description	Default
Q	ndarray	Zero-mean data matrix, shape \((N_x, N_t)\).	required
n_modes	int	Leading modes to compute. Default 20.	20
n_iter	int	Power-iteration steps. Default 4.	4
random_state	int	Seed for reproducibility.	None

Returns:

Name	Type	Description
Sigma	ndarray	Singular values (descending), shape \((N_\mathrm{modes},)\).
Psi	ndarray	Spatial modes (approximately orthonormal), shape \((N_x, N_\mathrm{modes})\).
Phi	ndarray	Temporal coefficients, shape \((N_\mathrm{modes}, N_t)\).

References

Halko, Martinsson & Tropp (2011). Finding structure with randomness. SIAM Review, 53(2), 217–288.

Source code in src/models/data_driven/autoencoders/pod_utils.py

def snapshot_pod_randomized(Q, n_modes=20, n_iter=4, random_state=None):
    r"""Randomized snapshot POD.

    Computes a truncated randomized SVD of $\mathbf{Q}$ (via
    ``sklearn.utils.extmath.randomized_svd``, falling back to a full
    ``numpy.linalg.svd`` if scikit-learn is unavailable),

    $$
    \mathbf{Q} \approx \mathbf{U}\,\mathbf{S}\,\mathbf{V}^\mathrm{T},
    $$

    and returns

    $$
    \boldsymbol{\Psi} = \mathbf{U}, \qquad
    \boldsymbol{\Sigma} = \mathbf{S} / \sqrt{N_t}, \qquad
    \boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q}.
    $$

    The $\boldsymbol{\Sigma} = \mathbf{S}/\sqrt{N_t}$ scaling matches the
    eigenvalue-based normalisation of `snapshot_pod`, since
    $\mathbf{Q}^\mathrm{T}\mathbf{Q}/N_t = \mathbf{V}(\mathbf{S}^2/N_t)\mathbf{V}^\mathrm{T}$.

    Parameters
    ----------
    Q : np.ndarray
        Zero-mean data matrix, shape $(N_x, N_t)$.
    n_modes : int
        Leading modes to compute. Default 20.
    n_iter : int
        Power-iteration steps. Default 4.
    random_state : int, optional
        Seed for reproducibility.

    Returns
    -------
    Sigma : np.ndarray
        Singular values (descending), shape $(N_\mathrm{modes},)$.
    Psi : np.ndarray
        Spatial modes (approximately orthonormal), shape $(N_x, N_\mathrm{modes})$.
    Phi : np.ndarray
        Temporal coefficients, shape $(N_\mathrm{modes}, N_t)$.

    References
    ----------
    Halko, Martinsson & Tropp (2011). Finding structure with randomness.
    *SIAM Review*, 53(2), 217–288.
    """
    N_x, N_t = Q.shape
    n_modes  = min(n_modes, N_x, N_t)
    try:
        from sklearn.utils.extmath import randomized_svd  # type: ignore[import-untyped]
        U, s, _ = randomized_svd(Q, n_components=n_modes,
                                 n_iter=n_iter, random_state=random_state)
    except ImportError:
        U_full, s_full, _ = np.linalg.svd(Q, full_matrices=False)
        U, s = U_full[:, :n_modes], s_full[:n_modes]
    Psi   = U
    Sigma = s / np.sqrt(N_t)
    Phi   = Psi.T @ Q
    return Sigma, Psi, Phi

romda.models.data_driven.autoencoders.pod_utils.spod_sieber(Q, Nf, kind='gaussian')

Sieber spectral POD — filtered correlation matrix.

Low-pass filters the temporal correlation matrix \(\mathbf{C} = \mathbf{Q}^\mathrm{T}\mathbf{Q}/N_t\) with the banded symmetric Toeplitz matrix \(\mathbf{G}\) (see _toeplitz_filter_matrix, built from a normalised kernel of half-width Nf, see _filter_kernel),

\[ \tilde{\mathbf{C}} = \mathbf{G}^\mathrm{T} \mathbf{C} \mathbf{G}, \]

then solves the same "method of snapshots" eigenvalue problem as snapshot_pod, with \(\tilde{\mathbf{C}}\) in place of \(\mathbf{C}\):

\[ \tilde{\mathbf{C}}\mathbf{A} = \mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda}), \qquad \boldsymbol{\Psi} = \frac{1}{\sqrt{N_t}}\, \mathbf{Q}\mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda})^{-1/2}, \qquad \boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q}, \qquad \boldsymbol{\Sigma} = \sqrt{\boldsymbol{\lambda}}. \]

Nf=0 skips the filtering step (\(\tilde{\mathbf{C}} = \mathbf{C}\)) and recovers standard snapshot POD exactly.

Parameters:

Name	Type	Description	Default
Q	ndarray	Zero-mean data matrix, shape \((N_x, N_t)\).	required
Nf	int	Filter half-width (0 recovers POD; \(N_t/2\) approaches the DFT).	required
kind	str	'gaussian', 'box' or 'hann'.	'gaussian'

Returns:

Name	Type	Description
Sigma	ndarray	Singular values (descending), shape \((N_t,)\).
Psi	ndarray	Spatial modes with orthonormal columns, shape \((N_x, N_t)\).
Phi	ndarray	Temporal SPOD coefficients, shape \((N_t, N_t)\).
C_tilde	ndarray	Filtered correlation matrix, shape \((N_t, N_t)\).

References

Sieber, Paschereit & Oberleithner (2016). Spectral proper orthogonal decomposition. J. Fluid Mech., 792, 798–828.

Source code in src/models/data_driven/autoencoders/pod_utils.py

def spod_sieber(Q, Nf, kind='gaussian'):
    r"""Sieber spectral POD — filtered correlation matrix.

    Low-pass filters the temporal correlation matrix
    $\mathbf{C} = \mathbf{Q}^\mathrm{T}\mathbf{Q}/N_t$ with the banded
    symmetric Toeplitz matrix $\mathbf{G}$ (see `_toeplitz_filter_matrix`,
    built from a normalised kernel of half-width ``Nf``, see
    `_filter_kernel`),

    $$
    \tilde{\mathbf{C}} = \mathbf{G}^\mathrm{T} \mathbf{C} \mathbf{G},
    $$

    then solves the same "method of snapshots" eigenvalue problem as
    `snapshot_pod`, with $\tilde{\mathbf{C}}$ in place of $\mathbf{C}$:

    $$
    \tilde{\mathbf{C}}\mathbf{A} = \mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda}),
    \qquad
    \boldsymbol{\Psi} = \frac{1}{\sqrt{N_t}}\,
    \mathbf{Q}\mathbf{A}\,\mathrm{diag}(\boldsymbol{\lambda})^{-1/2},
    \qquad
    \boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q},
    \qquad
    \boldsymbol{\Sigma} = \sqrt{\boldsymbol{\lambda}}.
    $$

    ``Nf=0`` skips the filtering step ($\tilde{\mathbf{C}} = \mathbf{C}$) and
    recovers standard snapshot POD exactly.

    Parameters
    ----------
    Q : np.ndarray
        Zero-mean data matrix, shape $(N_x, N_t)$.
    Nf : int
        Filter half-width (0 recovers POD; $N_t/2$ approaches the DFT).
    kind : str
        ``'gaussian'``, ``'box'`` or ``'hann'``.

    Returns
    -------
    Sigma : np.ndarray
        Singular values (descending), shape $(N_t,)$.
    Psi : np.ndarray
        Spatial modes with orthonormal columns, shape $(N_x, N_t)$.
    Phi : np.ndarray
        Temporal SPOD coefficients, shape $(N_t, N_t)$.
    C_tilde : np.ndarray
        Filtered correlation matrix, shape $(N_t, N_t)$.

    References
    ----------
    Sieber, Paschereit & Oberleithner (2016). Spectral proper orthogonal
    decomposition. *J. Fluid Mech.*, 792, 798–828.
    """
    N_x, N_t = Q.shape
    C        = (Q.T @ Q) / N_t
    if Nf == 0:
        C_tilde = C
    else:
        G       = _toeplitz_filter_matrix(N_t, Nf, kind)
        C_tilde = G.T @ C @ G
    lam, A    = np.linalg.eigh(C_tilde)
    idx       = lam.argsort()[::-1]
    lam, A    = lam[idx], A[:, idx]
    safe_lam  = np.where(lam > 0, lam, np.inf)
    Psi       = Q @ A / (np.sqrt(N_t) * np.sqrt(safe_lam))
    Psi[:, lam <= 0] = 0.0
    Phi       = Psi.T @ Q
    Sigma     = np.sqrt(np.where(lam > 0, lam, 0.0))
    return Sigma, Psi, Phi, C_tilde

romda.models.data_driven.autoencoders.pod_utils.spod_towne(Q, dt=1.0, n_fft=None, n_ovlp=None, window='hamming', weight=None, conf_level=0.95)

Spectral POD via Welch-averaged cross-spectral density (Towne et al., 2018).

Splits \(\mathbf{Q}\) into n_blks overlapping blocks of length n_fft (Welch's method, step n_fft - n_ovlp), windows and Fourier-transforms each block, then — at every frequency — solves a "method of snapshots" eigenvalue problem across blocks (analogous to snapshot_pod, but with blocks in place of time snapshots) to avoid ever forming the full \(N_x \times N_x\) cross-spectral density (CSD) matrix.

For block \(b = 0, \dots, n_\mathrm{blk}-1\) starting at snapshot \(b\,(n_\mathrm{fft}-n_\mathrm{ovlp})\), the windowed block DFT at frequency \(f_k\) is

\[ \hat{\mathbf{q}}^{(b)}_k = \frac{1}{n_\mathrm{fft}\,\bar{w}} \sum_{n=0}^{n_\mathrm{fft}-1} w[n]\,\mathbf{q}^{(b)}[n]\, e^{-\mathrm{i}2\pi kn/n_\mathrm{fft}}, \qquad \bar{w} = \mathrm{mean}(w), \]

with \(w\) the window. Stacking the blocks, \(\hat{\mathbf{Q}}_k = [\hat{\mathbf{q}}^{(1)}_k, \dots, \hat{\mathbf{q}}^{(n_\mathrm{blk})}_k] \in \mathbb{C}^{N_x \times n_\mathrm{blk}}\), the (weighted) cross-block Gram matrix and its eigendecomposition give the SPOD modes and modal energies at frequency \(f_k\):

\[ \mathbf{M}_k = \frac{1}{n_\mathrm{blk}}\, \hat{\mathbf{Q}}_k^\mathrm{H}\, \mathbf{W}\, \hat{\mathbf{Q}}_k, \qquad \mathbf{M}_k \boldsymbol{\Theta}_k = \boldsymbol{\Theta}_k\, \mathrm{diag}(\boldsymbol{\lambda}_k), \]

\[ \boldsymbol{\Psi}_k = \frac{1}{\sqrt{n_\mathrm{blk}}}\, \hat{\mathbf{Q}}_k\, \boldsymbol{\Theta}_k\, \mathrm{diag}(\boldsymbol{\lambda}_k)^{-1/2}, \]

with \(\mathbf{W} = \mathrm{diag}(\mathrm{weight})\) the spatial weight matrix (\(\mathbb{I}\) by default). The modal energy spectrum is \(L_k = \boldsymbol{\lambda}_k\), doubled (\(L_k = 2\boldsymbol{\lambda}_k\)) at interior frequency bins of a real, one-sided spectrum to account for the folded negative-frequency energy.

Parameters:

Name	Type	Description	Default
Q	ndarray	Zero-mean data matrix, shape \((N_x, N_t)\).	required
dt	float	Time step.	1.0
n_fft	int	Block/FFT length. Default \(2^{\lfloor \log_2 (N_t / 10) \rfloor}\).	None
n_ovlp	int	Block overlap. Default n_fft // 2.	None
window	str or ndarray	Window name (passed to scipy.signal.get_window) or an array of length n_fft.	'hamming'
weight	ndarray	Spatial integration weights \(\mathrm{diag}(\mathbf{W})\), shape \((N_x,)\). Uniform weights (no integration) by default.	None
conf_level	float	Target confidence level for the chi-squared-based interval Lc (Welch's method, nominally \(2\,n_\mathrm{blk}\) degrees of freedom).	0.95

Returns:

Name	Type	Description
L	ndarray	Modal energy spectrum, shape (n_freq, n_blks).
Psi	ndarray	Complex SPOD spatial modes, shape (n_freq, N_x, n_blks).
f	ndarray	Frequency vector, shape (n_freq,).
Lc	ndarray	Nominal confidence bounds [lower, upper] for L, shape (n_freq, n_blks, 2). See Notes.
info	dict	Effective n_fft, n_ovlp, n_blks, n_freq and window used.

Notes

Lc is computed from scipy.special.gammaincinv(1 - conf_level, n_blks) and scipy.special.gammaincinv(conf_level, n_blks), intended as a chi-squared confidence interval in the spirit of Welch's method. However, scipy.special.gammaincinv(a, y) requires its second argument \(y\) (a probability) in \([0, 1]\), whereas here it is called with \(y =\) n_blks (an integer \(\ge 2\)); numerically this returns nan for the n_blks values produced by this function. Treat Lc as unverified until this is checked against the intended formula.

References

Towne, Schmidt & Colonius (2018). Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. J. Fluid Mech., 847, 821–867.

Source code in src/models/data_driven/autoencoders/pod_utils.py

def spod_towne(Q, dt=1.0, n_fft=None, n_ovlp=None, window='hamming',
               weight=None, conf_level=0.95):
    r"""Spectral POD via Welch-averaged cross-spectral density (Towne et al., 2018).

    Splits $\mathbf{Q}$ into ``n_blks`` overlapping blocks of length ``n_fft``
    (Welch's method, step ``n_fft - n_ovlp``), windows and Fourier-transforms
    each block, then — at every frequency — solves a "method of snapshots"
    eigenvalue problem across blocks (analogous to `snapshot_pod`, but with
    blocks in place of time snapshots) to avoid ever forming the full
    $N_x \times N_x$ cross-spectral density (CSD) matrix.

    For block $b = 0, \dots, n_\mathrm{blk}-1$ starting at snapshot
    $b\,(n_\mathrm{fft}-n_\mathrm{ovlp})$, the windowed block DFT at
    frequency $f_k$ is

    $$
    \hat{\mathbf{q}}^{(b)}_k = \frac{1}{n_\mathrm{fft}\,\bar{w}}
    \sum_{n=0}^{n_\mathrm{fft}-1} w[n]\,\mathbf{q}^{(b)}[n]\,
    e^{-\mathrm{i}2\pi kn/n_\mathrm{fft}}, \qquad \bar{w} = \mathrm{mean}(w),
    $$

    with $w$ the ``window``. Stacking the blocks,
    $\hat{\mathbf{Q}}_k = [\hat{\mathbf{q}}^{(1)}_k, \dots,
    \hat{\mathbf{q}}^{(n_\mathrm{blk})}_k] \in \mathbb{C}^{N_x \times n_\mathrm{blk}}$,
    the (weighted) cross-block Gram matrix and its eigendecomposition give the
    SPOD modes and modal energies at frequency $f_k$:

    $$
    \mathbf{M}_k = \frac{1}{n_\mathrm{blk}}\,
    \hat{\mathbf{Q}}_k^\mathrm{H}\, \mathbf{W}\, \hat{\mathbf{Q}}_k,
    \qquad
    \mathbf{M}_k \boldsymbol{\Theta}_k
    = \boldsymbol{\Theta}_k\, \mathrm{diag}(\boldsymbol{\lambda}_k),
    $$

    $$
    \boldsymbol{\Psi}_k = \frac{1}{\sqrt{n_\mathrm{blk}}}\,
    \hat{\mathbf{Q}}_k\, \boldsymbol{\Theta}_k\,
    \mathrm{diag}(\boldsymbol{\lambda}_k)^{-1/2},
    $$

    with $\mathbf{W} = \mathrm{diag}(\mathrm{weight})$ the spatial weight matrix
    ($\mathbb{I}$ by default). The modal energy spectrum is
    $L_k = \boldsymbol{\lambda}_k$, doubled ($L_k = 2\boldsymbol{\lambda}_k$) at
    interior frequency bins of a real, one-sided spectrum to account for the
    folded negative-frequency energy.

    Parameters
    ----------
    Q : np.ndarray
        Zero-mean data matrix, shape $(N_x, N_t)$.
    dt : float
        Time step.
    n_fft : int, optional
        Block/FFT length. Default $2^{\lfloor \log_2 (N_t / 10) \rfloor}$.
    n_ovlp : int, optional
        Block overlap. Default ``n_fft // 2``.
    window : str or np.ndarray
        Window name (passed to ``scipy.signal.get_window``) or an array of
        length ``n_fft``.
    weight : np.ndarray, optional
        Spatial integration weights $\mathrm{diag}(\mathbf{W})$, shape
        $(N_x,)$. Uniform weights (no integration) by default.
    conf_level : float
        Target confidence level for the chi-squared-based interval `Lc`
        (Welch's method, nominally $2\,n_\mathrm{blk}$ degrees of freedom).

    Returns
    -------
    L : np.ndarray
        Modal energy spectrum, shape ``(n_freq, n_blks)``.
    Psi : np.ndarray
        Complex SPOD spatial modes, shape ``(n_freq, N_x, n_blks)``.
    f : np.ndarray
        Frequency vector, shape ``(n_freq,)``.
    Lc : np.ndarray
        Nominal confidence bounds ``[lower, upper]`` for `L`, shape
        ``(n_freq, n_blks, 2)``. See Notes.
    info : dict
        Effective ``n_fft``, ``n_ovlp``, ``n_blks``, ``n_freq`` and window used.

    Notes
    -----
    `Lc` is computed from ``scipy.special.gammaincinv(1 - conf_level, n_blks)``
    and ``scipy.special.gammaincinv(conf_level, n_blks)``, intended as a
    chi-squared confidence interval in the spirit of Welch's method. However,
    `scipy.special.gammaincinv(a, y)` requires its second argument $y$
    (a probability) in $[0, 1]$, whereas here it is called with $y =$
    ``n_blks`` (an integer $\ge 2$); numerically this returns ``nan`` for the
    ``n_blks`` values produced by this function. Treat `Lc` as unverified
    until this is checked against the intended formula.

    References
    ----------
    Towne, Schmidt & Colonius (2018). Spectral proper orthogonal decomposition and
    its relationship to dynamic mode decomposition and resolvent analysis.
    *J. Fluid Mech.*, 847, 821–867.
    """
    N_x, N_t = Q.shape
    is_real  = np.isrealobj(Q)

    if n_fft is None:
        n_fft = int(2 ** np.floor(np.log2(N_t / 10)))
    if n_ovlp is None:
        n_ovlp = n_fft // 2
    if n_ovlp >= n_fft:
        raise ValueError('n_ovlp must be < n_fft.')

    if isinstance(window, str):
        win = get_window(window, n_fft)
    else:
        win = np.asarray(window, dtype=float)
        if win.size != n_fft:
            raise ValueError(f'window length ({win.size}) must equal n_fft ({n_fft}).')
    win_norm = 1.0 / win.mean()
    win_col  = win[:, None]

    W = np.ones(N_x) if weight is None else np.asarray(weight, dtype=float).ravel()
    if W.size != N_x:
        raise ValueError('weight must have length N_x.')

    n_step = n_fft - n_ovlp
    n_blks = int(np.floor((N_t - n_ovlp) / n_step))
    if n_blks < 2:
        raise ValueError(
            f'Too few blocks ({n_blks}). Reduce n_fft/n_ovlp or use more snapshots.')

    if is_real:
        n_freq = n_fft // 2 + 1
        f      = np.arange(n_freq) / (n_fft * dt)
    else:
        n_freq = n_fft
        f      = np.fft.fftfreq(n_fft, d=dt)

    Q_hat = np.zeros((n_freq, N_x, n_blks), dtype=complex)
    for b in range(n_blks):
        i0     = b * n_step
        Q_blk  = Q[:, i0:i0 + n_fft]
        Q_fft  = np.fft.fft(Q_blk * win_col.T, axis=1) * win_norm / n_fft
        if is_real:
            Q_hat[:, :, b] = Q_fft[:, :n_freq].T
        else:
            Q_hat[:, :, b] = Q_fft.T

    L   = np.zeros((n_freq, n_blks))
    Psi = np.zeros((n_freq, N_x, n_blks), dtype=complex)
    for k in range(n_freq):
        Qf         = Q_hat[k]
        M          = (Qf * W[:, None]).conj().T @ Qf / n_blks
        lam, Theta = np.linalg.eigh(M)
        idx        = lam.argsort()[::-1]
        lam        = np.abs(lam[idx]); Theta = Theta[:, idx]
        Psi[k]     = Qf @ Theta / (np.sqrt(lam) * np.sqrt(n_blks))
        if is_real and 0 < k < n_freq - 1:
            L[k] = 2.0 * lam
        else:
            L[k] = lam

    xi2_up = 2 * gammaincinv(1 - conf_level, n_blks)
    xi2_lo = 2 * gammaincinv(    conf_level, n_blks)
    Lc     = np.stack([L * 2 * n_blks / xi2_lo,
                       L * 2 * n_blks / xi2_up], axis=-1)
    info = dict(n_fft=n_fft, n_ovlp=n_ovlp, n_blks=n_blks,
                window=win, n_freq=n_freq)
    return L, Psi, f, Lc, info

romda.models.data_driven.autoencoders.pod_utils.spod_towne_reconstruct(Psi, A_blk, n_fft, n_ovlp, N_t)

Reconstruct snapshots from Towne SPOD modes [stub — not yet implemented].

Parameters:

Name	Type	Description	Default
Psi	ndarray	SPOD spatial modes, as returned by spod_towne.	required
A_blk	ndarray	Block expansion coefficients.	required
n_fft	int	Block/FFT length used by spod_towne.	required
n_ovlp	int	Block overlap used by spod_towne.	required
N_t	int	Number of snapshots in the reconstructed series.	required

Raises:

Type	Description
NotImplementedError	Always — inverse SPOD is not yet implemented.

References

Nekkanti & Schmidt (2021). Frequency-time analysis, low-rank reconstruction and denoising of turbulent flows using SPOD. J. Fluid Mech., 926, A26.

Source code in src/models/data_driven/autoencoders/pod_utils.py

def spod_towne_reconstruct(Psi, A_blk, n_fft, n_ovlp, N_t):  # noqa: ARG001
    """
    Reconstruct snapshots from Towne SPOD modes  [stub — not yet implemented].

    Parameters
    ----------
    Psi : np.ndarray
        SPOD spatial modes, as returned by `spod_towne`.
    A_blk : np.ndarray
        Block expansion coefficients.
    n_fft : int
        Block/FFT length used by `spod_towne`.
    n_ovlp : int
        Block overlap used by `spod_towne`.
    N_t : int
        Number of snapshots in the reconstructed series.

    Raises
    ------
    NotImplementedError
        Always — inverse SPOD is not yet implemented.

    References
    ----------
    Nekkanti & Schmidt (2021). Frequency-time analysis, low-rank reconstruction
    and denoising of turbulent flows using SPOD. *J. Fluid Mech.*, 926, A26.
    """
    raise NotImplementedError(
        "Full reconstruction (inverse SPOD) is not yet implemented. "
        "See Nekkanti & Schmidt (JFM 2021) for details.")

romda.models.data_driven.autoencoders.pod_utils.print_spod_towne_summary(info)

Pretty-print the block/frequency parameters of a spod_towne run.

Parameters:

Name	Type	Description	Default
info	dict	The info dictionary returned by spod_towne (must contain n_fft, n_ovlp, n_blks and n_freq).	required

Source code in src/models/data_driven/autoencoders/pod_utils.py

def print_spod_towne_summary(info):
    """
    Pretty-print the block/frequency parameters of a `spod_towne` run.

    Parameters
    ----------
    info : dict
        The ``info`` dictionary returned by `spod_towne` (must contain
        ``n_fft``, ``n_ovlp``, ``n_blks`` and ``n_freq``).
    """
    print('SPOD (Towne / Welch) parameters')
    print('────────────────────────────────')
    print(f'  Snapshots per block (n_fft)  : {info["n_fft"]}')
    print(f'  Block overlap (n_ovlp)       : {info["n_ovlp"]}')
    print(f'  Number of blocks             : {info["n_blks"]}')
    print(f'  Resolved frequencies         : {info["n_freq"]}')

romda.models.data_driven.autoencoders.pod_utils.energy_fraction(Sigma)

Relative energy fraction and cumulative energy per mode.

From the eigenvalues \(\lambda_j = \Sigma_j^2\) of the (possibly filtered) temporal correlation matrix,

\[ \mathrm{rel}_j = \frac{\lambda_j}{\sum_k \lambda_k}, \qquad \mathrm{cum}_j = \sum_{k \le j} \mathrm{rel}_k. \]

Parameters:

Name	Type	Description	Default
Sigma	ndarray	Singular values, as returned by snapshot_pod, snapshot_pod_randomized or spod_sieber.	required

Returns:

Name	Type	Description
rel	ndarray	Relative energy per mode (sums to 1).
cum	ndarray	Cumulative relative energy.

Source code in src/models/data_driven/autoencoders/pod_utils.py

def energy_fraction(Sigma):
    r"""
    Relative energy fraction and cumulative energy per mode.

    From the eigenvalues $\lambda_j = \Sigma_j^2$ of the (possibly filtered)
    temporal correlation matrix,

    $$
    \mathrm{rel}_j = \frac{\lambda_j}{\sum_k \lambda_k}, \qquad
    \mathrm{cum}_j = \sum_{k \le j} \mathrm{rel}_k.
    $$

    Parameters
    ----------
    Sigma : np.ndarray
        Singular values, as returned by `snapshot_pod`,
        `snapshot_pod_randomized` or `spod_sieber`.

    Returns
    -------
    rel : np.ndarray
        Relative energy per mode (sums to 1).
    cum : np.ndarray
        Cumulative relative energy.
    """
    lam = Sigma ** 2
    rel = lam / lam.sum()
    return rel, np.cumsum(rel)