Data-driven models

Summary

File: src/models/data_driven/. All use DiscreteIntegrator. ESN_model mixes in EchoStateNetwork — the reservoir core from the external echostatenetwork package (re-exported by romda.models.data_driven for convenience; see its own documentation for the reservoir API) — and POD_ESN mixes in ESN_model and POD.

Class	Key parameters	Notes
ESN_model	N_units, rho, sigma_in, N_wash	Inherits EchoStateNetwork; requires training data at init
POD_ESN	N_modes, sensor_locations	Inherits ESN_model + POD; sensor placement via QR
LinearModel	F (transition matrix), Q_noise	\(\boldsymbol{\psi}_{t+1} = \mathbf{F}\boldsymbol{\psi}_t + \boldsymbol{\eta}_t\)

The Projector hierarchy (POD/SPOD) and the standalone decomposition functions (pod_utils.py) live in the autoencoders/ subpackage and are documented below. Import everything from romda.models.data_driven (or its autoencoders subpackage).

Open-loop (training/washout) vs. closed-loop (forecasting) reservoir configurations.

---

Forecast models

romda.models.data_driven.esn.ESN_model(dt, **kwargs)

Bases: EchoStateNetwork, Model

Echo state network as a data-driven forecast model.

Wraps the EchoStateNetwork reservoir with the Model interface (state history, discrete integrator, observation operator), so a trained ESN can be used as the forecast model of an Ensemble — or as the forecaster inside ESN_bias. The model state is \([\mathbf{u}; \mathbf{r}]\): the physical outputs and the reservoir state. Training data is mandatory at construction (the network trains itself unless a cached configuration is found).

Parameters:

Name	Type	Description	Default
dt	float	Output time step (the internal ESN step is dt * upsample).	required
**kwargs		Supported keys include: data : np.ndarray, training data, shape \((L, N_t, N_\mathrm{dim})\) (\(L\) segments/experiments, \(N_t\) time steps spanning train + validation + test, \(N_\mathrm{dim}\) state dimensions). Required unless y0 is given (with pre-trained matrices, e.g. Wout). y0 : np.ndarray, initial state, used only if data is not given. plot_training : bool, whether to plot the training data and convergence. Default True. ESN hyperparameters (N_units, N_wash, rho_range, ...) and Model options.	{}

Source code in src/models/data_driven/esn.py

def __init__(self,
             dt,
             **kwargs):
    # See the class docstring for the meaning of dt and **kwargs.
    data = kwargs.pop('data', None)
    y0 = kwargs.pop('y0', None)
    plot_training = kwargs.pop('plot_training', True)

    # =================== STEP 0: ESTIMABLE INPUT PARAMETERS ======================
    # A parametric ESN is conditioned on a physical parameter vector that is
    # normally *given* at forecast time (`EchoStateNetwork.input_parameters`).
    # Naming the entries in `param_names` mirrors each one as a `Model` parameter:
    # it can be listed in `est_alpha`, joins the augmented state alongside the
    # reservoir (and the `Wout` singular values), and is read back per member
    # before every forecast (see `time_step`). `param_values` are the nominal
    # values (the ensemble mean at initialization). Avoid the names of the ESN's
    # own hyperparameters (`rho`/`sigma_in`/`tikh`).
    self.param_names = tuple(kwargs.pop('param_names', ()))
    param_labels = kwargs.pop('param_labels', None)
    values = np.atleast_1d(np.asarray(
        kwargs.pop('param_values', np.zeros(len(self.param_names))), dtype=float))
    # kept as attributes so the config store can save/rebuild a parametric ESN
    self.param_values = tuple(float(v) for v in values)
    self.param_labels = tuple(param_labels) if param_labels else None
    for name, val in zip(self.param_names, values):
        setattr(self, name, float(val))
    # `Model.__init__` builds alpha0 from `params`, so register them before it runs
    self.params = list(self.params) + list(self.param_names)

    # =================== STEP 1: EchoStateNetwork INITIALIZATION ======================

    [setattr(self, key, kwargs.pop(key)) for key in list(kwargs.keys()) if key in vars(ESN_model)]

    if data is not None:
        assert isinstance(data, np.ndarray), f"Expected data to be a numpy array, got {type(data)}"
        data = self._process_initialization_data(data, dt, kwargs) # type: np.ndarray # with shape (L, Nt, Ndim)
        y0 = data[0, 0]
    elif y0 is None:
        raise ValueError('Either training data or initial state y0 must be provided to initialize the ESN_model.')

    initial_dict = {key: kwargs.pop(key) for key in list(kwargs.keys()) if key in vars(EchoStateNetwork)}
    EchoStateNetwork.__init__(self,
                            y=y0,
                            dt=dt,
                            **initial_dict)


    # =================== STEP 2: EchoStateNetwork TRAINING ======================
    # Train the network if not already trained
    if not self.trained:
        print('Training ESN model...')
        if plot_training:
            self.plot_training_data(case=self, train_data=data, dt=dt)

        self.train(train_data=data, plot_training=plot_training, **kwargs)

        # save validation data for initialization
        Y_wtv = self._split_and_format_data(data)[1]
        self.validation_data = Y_wtv[-(self.N_wash + self.N_val):]


    # ================== STEP 3: DEFINE INITIAL STATE & PARAMS ======================
    # `Model.Nq = 1` is a class attribute, so `hasattr` is always True -- check the
    # instance dict instead, otherwise every ESN_model silently ends up with Nq=1
    # (POD_ESN sets its own Nq as an instance attribute before calling this).
    if 'Nq' not in vars(self):
        self.Nq = len(self.observed_idx)  # Number of observed dimensions (for the physical state)

    psi0 = self.initialize_from_val_data()  # shape (Ndim + N_units + Na, m)

    # Initialise SVD Wout terms if required
    if self.Wout_svd:
        [self.Wout_U, self.Wout_Sigma0, self.Wout_Vh] = sla.svd(self.Wout, full_matrices=False)
        self.Wout_Sigma = self.Wout_Sigma0

    # =================== STEP 4: Model INITIALIZATION ======================
    Model.__init__(self,
                   dt=dt,
                   psi0=psi0,
                   integrator_class=DiscreteIntegrator,
                   **kwargs)

    if self.param_names:
        # `Model.alpha_labels` defaults to a label->name mapping, so plotting an
        # estimated parameter by name raises KeyError until real labels are set.
        self.alpha_labels = dict(zip(self.param_names,
                                     param_labels or [f'${n}$' for n in self.param_names]))

t_transient property

float: Total time spanning training + validation + test, t_train + t_val + t_test.

dt_step property

float: Integrator time step, dt_ESN (the ESN advances in closed loop at its own upsampled time step, not dt).

t_CR property

float: Characteristic response time used by the base Model, aliased to t_val.

Wout_U property writable

np.ndarray: Left singular vectors of Wout (from scipy.linalg.svd(Wout, full_matrices=False)), shape Wout.shape i.e. (N_units + 1, N_dim). Only used/set when Wout_svd is True.

Wout_Vh property writable

np.ndarray: Right singular vectors of Wout (transposed), shape (N_dim, N_dim). Only used/set when Wout_svd is True.

Wout_Sigma property writable

np.ndarray: Ensemble of diagonal singular-value matrices used to reconstruct Wout as \(\mathbf{W}_\mathrm{out} \approx \mathbf{U}\,\boldsymbol{\Sigma}\,\mathbf{V}^\mathrm{h}\) (see reservoir_to_physical), shape (m, N_dim, N_dim). If Wout_svd, recomputed from the current svd_i ensemble parameters on every access (via alpha_to_Sigma); otherwise held fixed at whatever was last set.

alpha_to_Sigma property

np.ndarray: Per-ensemble-member diagonal singular-value matrices built from the current svd_i parameter estimates (get_alpha_matrix), falling back to the corresponding Wout_Sigma0 singular value for any svd_i not in est_alpha. Shape (m, N_dim, N_dim).

get_alpha_matrix property

np.ndarray: Current ensemble parameter estimates (est_alpha), shape (len(est_alpha), m), read from get_alpha.

Wout_Sigma0 property writable

np.ndarray: Reference (initial) singular values of Wout, shape (N_dim,), as computed by scipy.linalg.svd when Wout_svd is enabled.

N_ens property

int: Ensemble size, read from ensemble['m'] if an ensemble configuration is set, otherwise the trailing dimension of current_state.

state_labels property

list of str: LaTeX labels for the state vector, \(u_1, \dots, u_{N_\mathrm{dim}}\) (physical outputs) followed by \(r_1, \dots, r_{N_\mathrm{units}}\) (reservoir units).

obs_labels property

list of str: LaTeX labels for the observed physical outputs (observed_idx).

reservoir_state property

np.ndarray: Reservoir-state block of current_state, shape (N_units, m).

modify_settings(**kwargs)

Update existing attributes in place, switching to the SVD parametrization of Wout (see Wout_svd) if 'Wout' is requested as an ensemble parameter to estimate (est_alpha).

Parameters:

Name	Type	Description	Default
**kwargs		Attribute name/value pairs to set; each name must already exist on the instance.	{}

Returns:

Type	Description
None	Updates attributes (and, if applicable, est_alpha/alpha_labels/ alpha_lims/M) in place.

Raises:

Type	Description
ValueError	If a key in kwargs is not an existing attribute.

Source code in src/models/data_driven/esn.py

def modify_settings(self, **kwargs):
    """Update existing attributes in place, switching to the SVD parametrization
    of `Wout` (see `Wout_svd`) if ``'Wout'`` is requested as an ensemble parameter
    to estimate (`est_alpha`).

    Parameters
    ----------
    **kwargs
        Attribute name/value pairs to set; each name must already exist on the
        instance.

    Returns
    -------
    None
        Updates attributes (and, if applicable, `est_alpha`/`alpha_labels`/
        `alpha_lims`/`M`) in place.

    Raises
    ------
    ValueError
        If a key in `kwargs` is not an existing attribute.
    """
    for key, val in kwargs.items():
        if hasattr(self, key):
            setattr(self, key, val)
        else:
            raise ValueError(f'Key {key} not in ESN_model class')


    if self.ensemble_cfg is not None:
        # If Wout is being estimated, we need to update the est_alpha list to include the SVD components
        # and remove Wout. We do not directly estimate Wout, but rather its singular values.
        if 'Wout' in self.est_alpha:
            self.est_alpha = [a for a in self.est_alpha if a != 'Wout'] + \
                             [f'svd_{qi}' for qi in range(self.N_dim)]

        # The switch must key on the svd names, not on 'Wout': callers that build
        # the initial ensemble themselves (esn_ensemble) pass est_alpha already
        # expanded, and without Wout_svd the read-out keeps using the fixed Wout —
        # the estimated singular values never act on the forecast.
        svd_keys = [a for a in self.est_alpha if a.startswith('svd_')]
        if svd_keys:
            if not self.Wout_svd:
                self.Wout_svd = True
                [self.Wout_U, self.Wout_Sigma0, self.Wout_Vh] = sla.svd(self.Wout, full_matrices=False)
                self.Wout_Sigma = self.Wout_Sigma0
            # the setters merge, so labels/lims of other estimated parameters survive
            self.alpha_labels = {key: f'$\\sigma_{{{key.split("_")[1]}}}$' for key in svd_keys}
            self.alpha_lims = {key: (None, None) for key in svd_keys}
    self.M = None

init_ensemble(m=10, est_alpha=[], std_alpha=0.001, distribution_alpha='uniform', regimes=None, seed=0, measured=None, ensemble_psi0=None, **kwargs)

ESN override of Model.init_ensemble: members start from reservoir states visited during training (initialize_from_val_data) — the generic transient + multiplicative std_phi path pushes r outside the tanh range and the closed loop blows up.

regimes ((N_param, L), parametric ESN): washes each member out on one training segment and starts it from that segment's parameters, keeping state and parameter consistent — paired independently, members leave the learned attractor and diverge. 'Wout' in est_alpha estimates the read-out singular values (one svd_i per output dimension); measured restricts the observation operator to those state components.

Parameter perturbations draw from this instance's rng; an estimator builds the ensemble on its own copy (EnsembleEstimator copies parent_model first), so repeated builds from one parent network draw identical perturbations — the parent's rng is deliberately untouched.

Source code in src/models/data_driven/esn.py

def init_ensemble(self, m=10, est_alpha=[], std_alpha=0.001,
                  distribution_alpha='uniform', regimes=None, seed=0,
                  measured=None, ensemble_psi0=None, **kwargs):
    """ESN override of `Model.init_ensemble`: members start from reservoir states
    visited during training (`initialize_from_val_data`) — the generic transient +
    multiplicative `std_phi` path pushes `r` outside the tanh range and the closed
    loop blows up.

    `regimes` (`(N_param, L)`, parametric ESN): washes each member out on one
    training segment *and* starts it from that segment's parameters, keeping state
    and parameter consistent — paired independently, members leave the learned
    attractor and diverge. ``'Wout'`` in `est_alpha` estimates the read-out
    singular values (one ``svd_i`` per output dimension); `measured` restricts
    the observation operator to those state components.

    Parameter perturbations draw from *this* instance's `rng`; an estimator
    builds the ensemble on its own copy (`EnsembleEstimator` copies
    `parent_model` first), so repeated builds from one parent network draw
    identical perturbations — the parent's `rng` is deliberately untouched."""
    if ensemble_psi0 is None:
        nominal = {}
        if isinstance(std_alpha, dict) and 'Wout' in std_alpha:
            std_alpha = dict(std_alpha)
            std_alpha.update({f'svd_{qi}': std_alpha.pop('Wout')
                              for qi in range(self.N_dim)})
        if isinstance(std_alpha, dict) and not est_alpha:
            # base-class convention: an empty est_alpha with a dict std_alpha
            # estimates every parameter in the dict; an explicit est_alpha is
            # honored (the dict may carry spreads for more parameters)
            est_alpha = sorted(std_alpha)
        if 'Wout' in est_alpha:
            svd_names = [f'svd_{qi}' for qi in range(self.N_dim)]
            est_alpha = [a for a in est_alpha if a != 'Wout'] + svd_names
            nominal = dict(zip(svd_names, sla.svd(self.Wout, full_matrices=False)[1]))

        if regimes is None:
            phi, alpha = self.initialize_from_val_data(N_ens=m), {}
        else:
            phi, values = self._regime_matched_init(m, np.atleast_2d(regimes), seed=seed)
            alpha = dict(zip(self.param_names, values))

        other = [a for a in est_alpha if a not in alpha]
        if other:
            means = [nominal.get(a, getattr(self, a, None)) for a in other]
            assert not any(v is None for v in means), \
                f'init_ensemble: no nominal value to seed the ensemble for {other}'
            sub_std = ({a: std_alpha[a] for a in other}
                       if isinstance(std_alpha, dict) else std_alpha)
            alpha.update(zip(other, mean_vector_to_ensemble(
                self.rng, np.array(means, dtype=float), sub_std, m,
                method=distribution_alpha)))

        ensemble_psi0 = phi if not est_alpha else np.concatenate(
            [phi, np.stack([alpha[name] for name in est_alpha])], axis=0)

    out = super().init_ensemble(m=m, est_alpha=est_alpha, std_alpha=std_alpha,
                                distribution_alpha=distribution_alpha,
                                ensemble_psi0=ensemble_psi0, **kwargs)

    if measured is not None:
        if np.isscalar(measured):   # int count or digit string: resolve to indices
            from romda.observations import measured_idx
            idx = measured_idx(measured, self.N_dim)
        else:
            idx = [int(k) for k in measured]
        if idx != list(self.observed_idx):
            # Rows of `M` for the measured components only. Written to the private
            # attribute: the public setter requires exactly Nq rows, and a
            # restricted operator deliberately has fewer.
            obs = list(self.observed_idx)
            self._M = self.M[[obs.index(k) for k in idx]]
    return out

initialize_from_val_data(N_ens=1, seed=0)

Initialize the ESN state (physical output and reservoir) from an open-loop washout over a random time window of validation_data, so a fresh ensemble starts from an on-attractor reservoir state rather than zeros.

Parameters:

Name	Type	Description	Default
N_ens	int	Number of ensemble members (random washout windows) to draw. Default 1.	1
seed	int	Random seed for selecting the washout windows. Overridden by self.seed if set.	0

Returns:

Type	Description
ndarray	Initial full state (built via build_psi), shape (N_dim + N_units [+ Na], N_ens).

Raises:

Type	Description
AssertionError	If validation_data has not been set (i.e. before training).

Source code in src/models/data_driven/esn.py

def initialize_from_val_data(self, N_ens=1, seed=0):
    """Initialize the ESN state (physical output and reservoir) from an
    open-loop washout over a random time window of `validation_data`, so a fresh
    ensemble starts from an on-attractor reservoir state rather than zeros.

    Parameters
    ----------
    N_ens : int
        Number of ensemble members (random washout windows) to draw. Default 1.
    seed : int
        Random seed for selecting the washout windows. Overridden by `self.seed`
        if set.

    Returns
    -------
    np.ndarray
        Initial full state (built via `build_psi`), shape
        ``(N_dim + N_units [+ Na], N_ens)``.

    Raises
    ------
    AssertionError
        If `validation_data` has not been set (i.e. before training).
    """
    assert self.validation_data is not None

    data = self.validation_data.copy()

    if hasattr(self, 'seed'):
        seed = self.seed
    rng0 = np.random.default_rng(seed)


    # initialise state with a random sample from test data
    u_init, r_init = np.empty((self.N_dim, N_ens)), np.empty((self.N_units, N_ens))

    # Random time windows and dimension
    if data.shape[0] == 1:
        dim_ids = [0] * N_ens
    else:
        # Choose a random dimension from the data
        replace = N_ens >= data.shape[0]
        dim_ids = rng0.choice(data.shape[0], size=N_ens, replace=replace)

    # Choose random time indices from the data (with replacement when there are
    # fewer washout windows than members, e.g. a short validation record)
    n_windows = data.shape[1] - self.N_wash
    t_ids = rng0.choice(n_windows, size=N_ens, replace=N_ens >= n_windows)

    # validation_data is indexed by segment (dim_i matches a column of
    # input_parameters); condition the washout of each segment on its own
    # parameter vector, then restore the full (N_param, L) array afterwards.
    original_input_parameters = self.input_parameters

    for ii, ti, dim_i in zip(range(N_ens), t_ids, dim_ids):
        u_wash = data[dim_i, ti:ti+self.N_wash]
        r_open = np.zeros((self.N_units, 1))
        u_open = np.zeros((self.N_dim, 1))
        if original_input_parameters is not None:
            self.input_parameters = original_input_parameters[:, dim_i]
        # Open-loop reservoir
        for u_in in u_wash:
            u_open, r_open = self._single_step(u_in, r_open)

        #store final state into the initialization arrays
        u_init[:, ii] = u_open.squeeze()
        r_init[:, ii] = r_open.squeeze()

    self.input_parameters = original_input_parameters

    # Set physical and reservoir states as ensembles
    return self.build_psi(u=u_init, r=r_init)

closed_loop(data, n_steps, input_parameters=None)

Open-loop washout on data[:N_wash] (sampled at dt_ESN), then a closed-loop forecast; returns (prediction, target), both (n_steps, N_dim).

The ESN maps u_t to u_{t+1}, so the last washout step already predicts data[N_wash] — an off-by-one here costs a full dt_ESN of drift.

Source code in src/models/data_driven/esn.py

def closed_loop(self, data, n_steps, input_parameters=None):
    """Open-loop washout on `data[:N_wash]` (sampled at `dt_ESN`), then a closed-loop
    forecast; returns ``(prediction, target)``, both `(n_steps, N_dim)`.

    The ESN maps `u_t` to `u_{t+1}`, so the *last* washout step already predicts
    `data[N_wash]` — an off-by-one here costs a full `dt_ESN` of drift."""
    original = self.input_parameters
    if input_parameters is not None:
        self.input_parameters = np.asarray(input_parameters, dtype=float).reshape(-1, 1)

    try:
        r = np.zeros((self.N_units, 1))
        u = np.zeros((self.N_dim, 1))
        for u_in in data[:self.N_wash]:
            u, r = self.step(self.outputs_to_inputs(np.asarray(u_in)[:, np.newaxis]), r)

        pred = np.empty((n_steps, self.N_dim))
        pred[0] = u[:, 0]
        for i in range(1, n_steps):
            u, r = self._single_step(u, r)
            pred[i] = u[:, 0]
    finally:
        # the training-time (N_param, L) array must survive a single-regime forecast
        self.input_parameters = original

    return pred, np.asarray(data[self.N_wash:self.N_wash + n_steps])

reset_ESN(data, u0=None, plot_training=False, **kwargs)

Reinitialize and retrain the underlying EchoStateNetwork from scratch on new data, then reset the Model state/history around the freshly trained network.

Parameters:

Name	Type	Description	Default
data	ndarray	New training data, shape (L, Nt, N_dim); forwarded to train.	required
u0	ndarray	Initial physical state for the new EchoStateNetwork. Defaults to the physical state corresponding to the current reservoir_state.	None
plot_training	bool	Whether to plot the (re-)training process. Default False.	False
**kwargs		ESN hyperparameters (forwarded to EchoStateNetwork.__init__), training options (forwarded to train) and Model reset options (forwarded to reset_model).	{}

Returns:

Type	Description
None	Reinitializes self in place.

Source code in src/models/data_driven/esn.py

def reset_ESN(self, data, u0=None, plot_training=False, **kwargs):
    """Reinitialize and retrain the underlying `EchoStateNetwork` from scratch
    on new `data`, then reset the `Model` state/history around the freshly
    trained network.

    Parameters
    ----------
    data : np.ndarray
        New training data, shape ``(L, Nt, N_dim)``; forwarded to `train`.
    u0 : np.ndarray, optional
        Initial physical state for the new `EchoStateNetwork`. Defaults to the
        physical state corresponding to the current `reservoir_state`.
    plot_training : bool
        Whether to plot the (re-)training process. Default False.
    **kwargs
        ESN hyperparameters (forwarded to `EchoStateNetwork.__init__`), training
        options (forwarded to `train`) and `Model` reset options (forwarded to
        `reset_model`).

    Returns
    -------
    None
        Reinitializes `self` in place.
    """
    if u0 is None:
        u0 = self.reservoir_to_physical(self.reservoir_state)

    EchoStateNetwork.__init__(self,
                              y=u0,
                              dt=self.dt,
                              figs_folder=self.results_folder,
                              **kwargs)
    # Train the network
    possible_args = inspect.getfullargspec(self.train)[0]
    train_args = {key: val for key, val in kwargs.items() if key in possible_args}

    # Train network
    self.train(train_data=data, plot_training=plot_training, **train_args)


    # Reset model class
    kwargs['psi0'] = self.build_psi()
    self.reset_model(**kwargs)

get_observables(Nt=1, **kwargs)

Observables are the observed physical outputs, which need not be the leading rows of psi (the base Model assumes psi[:Nq]).

Parameters:

Name	Type	Description	Default
Nt	int	Number of trailing history steps to return. Default 1.	1
**kwargs		Unused; accepted for interface compatibility.	{}

Returns:

Type	Description
ndarray	Observed outputs, shape (Nq, m) if Nt == 1 else (Nt, Nq, m).

Source code in src/models/data_driven/esn.py

def get_observables(self, Nt=1, **kwargs):
    """Observables are the *observed* physical outputs, which need not be the
    leading rows of psi (the base Model assumes `psi[:Nq]`).

    Parameters
    ----------
    Nt : int
        Number of trailing history steps to return. Default 1.
    **kwargs
        Unused; accepted for interface compatibility.

    Returns
    -------
    np.ndarray
        Observed outputs, shape ``(Nq, m)`` if ``Nt == 1`` else ``(Nt, Nq, m)``.
    """
    if Nt == 1:
        return self.hist[-1, self.observed_idx, :]
    return self.hist[-Nt:, self.observed_idx, :]

reservoir_to_physical(r)

Convert reservoir states to physical outputs via Wout (overrides EchoStateNetwork.reservoir_to_physical to also support the SVD parametrization of Wout, see Wout_svd).

When Wout_svd is False: \(\mathbf{u} = \mathbf{W}_\mathrm{out}^\mathrm{T}[\mathbf{r}; b_\mathrm{out}]\), as in the base class. When True, \(\mathbf{W}_\mathrm{out}\) is reconstructed (per ensemble member, if r has m members) from \(\mathbf{U}\,\boldsymbol{\Sigma}\,\mathbf{V}^\mathrm{h}\) (Wout_U, Wout_Sigma, Wout_Vh) before the same read-out is applied; if r does not have exactly m members (e.g. a single averaged state), Wout_Sigma is averaged over the ensemble first.

Parameters:

Name	Type	Description	Default
r	ndarray	Reservoir state, shape (N_units, N_ens).	required

Returns:

Type	Description
ndarray	Physical output, shape (N_dim, N_ens).

Source code in src/models/data_driven/esn.py

def reservoir_to_physical(self, r):
    r"""Convert reservoir states to physical outputs via `Wout` (overrides
    `EchoStateNetwork.reservoir_to_physical` to also support the SVD
    parametrization of `Wout`, see `Wout_svd`).

    When `Wout_svd` is False: $\mathbf{u} = \mathbf{W}_\mathrm{out}^\mathrm{T}[\mathbf{r}; b_\mathrm{out}]$,
    as in the base class. When True, $\mathbf{W}_\mathrm{out}$ is reconstructed
    (per ensemble member, if `r` has `m` members) from
    $\mathbf{U}\,\boldsymbol{\Sigma}\,\mathbf{V}^\mathrm{h}$ (`Wout_U`, `Wout_Sigma`,
    `Wout_Vh`) before the same read-out is applied; if `r` does not have exactly
    `m` members (e.g. a single averaged state), `Wout_Sigma` is averaged over the
    ensemble first.

    Parameters
    ----------
    r : np.ndarray
        Reservoir state, shape ``(N_units, N_ens)``.

    Returns
    -------
    np.ndarray
        Physical output, shape ``(N_dim, N_ens)``.
    """
    bias_out = self.bias_out * np.ones((1, r.shape[-1]))
    r_aug = np.concatenate((r, bias_out), axis=0)

    if not self.Wout_svd:
        return np.dot(r_aug.T, self.Wout).T
    else:

        if r.shape[-1] == self.m:
            Wout = np.einsum('ij,kjl,lm->imk', self.Wout_U, self.Wout_Sigma, self.Wout_Vh)
            return np.einsum('ij,ikj->kj', r_aug, Wout)
        else:
            # average the alpha values
            print('Averaging Wout_Sigma for reservoir_to_physical')
            Wout_Sigma_avg = np.mean(self.Wout_Sigma, axis=0)
            Wout = np.dot(self.Wout_U, np.dot(Wout_Sigma_avg, self.Wout_Vh))
            return np.dot(r_aug.T, Wout).T

time_step(Nt=10, averaged=False)

Advance the ESN in closed loop.

Parameters:

Name	Type	Description	Default
Nt	int	Number of forecast steps (in physical time steps, not dt_ESN).	10
averaged	bool	If True, the ensemble is forecast as its mean plus frozen deviations; otherwise each member is forecast individually.	False

Returns:

Type	Description
tuple	(psi, t) — forecasted state of shape (Nt+1, N, m) and the corresponding times.

Source code in src/models/data_driven/esn.py

def time_step(self, Nt=10, averaged=False):
    """Advance the ESN in closed loop.

    Parameters
    ----------
    Nt : int
        Number of forecast steps (in physical time steps, not ``dt_ESN``).
    averaged : bool
        If True, the ensemble is forecast as its mean plus frozen deviations;
        otherwise each member is forecast individually.

    Returns
    -------
    tuple
        ``(psi, t)`` — forecasted state of shape ``(Nt+1, N, m)`` and the
        corresponding times.
    """

    assert self.trained, 'ESN model not trained'

    if self.param_names:
        # Condition each member's forecast on its own current parameter value.
        # `get_alpha` returns one dict per member: the estimated value when the
        # parameter is in `est_alpha`, otherwise the nominal `alpha0` one (change
        # it in place with ``esn.alpha0[name] = ...`` to forecast at a fixed
        # parameter). Either way the `(N_param, m)` array replaces the
        # `(N_param, L)` one left over from training.
        alpha = self.get_alpha()
        self.input_parameters = np.array([[a[name] for a in alpha]
                                          for name in self.param_names])

    # 1. get initial condition


    t = np.round(self.current_time + np.arange(0, Nt + 1) * self.dt_ESN, self.precision_t)
    psi0 = self.current_state
    u, r_out = np.empty((Nt + 1, self.N_dim, self.m)), np.empty((Nt + 1, self.N_units, self.m))
    u[0], r_out[0] = self.unbuild_psi(psi0)

    if averaged:
        # Mean state
        u_m, r_m = (np.mean(yy[0], axis=-1, keepdims=True) for yy in [u, r_out])
        u_dev, r_dev = u[0] - u_m[0], r_out[0] - r_m[0]


        for i in range(Nt):
            u_m, r_m = self._single_step(u_m, r_m)
            u[i+1] = u_m + u_dev
            r_out[i+1] = r_m + r_dev

    else:

        for i in range(Nt):
            u[i+1], r_out[i+1] = self._single_step(u[i], r_out[i])


    psi = self.build_psi(u=u, r=r_out)

    return psi, t

build_psi(u=None, r=None)

Assemble the full model state from physical output u and reservoir state r: concatenate([u, r, alpha]) along the state axis, keeping only u or only r if update_state/update_reservoir is False, and appending the ensemble parameter block (get_alpha_matrix) if Na > 0.

Parameters:

Name	Type	Description	Default
u	ndarray	Physical output, shape (N_dim, m) or (Nt, N_dim, m). Defaults to reservoir_to_physical(r).	None
r	ndarray	Reservoir state, shape (N_units, m) or (Nt, N_units, m). Defaults to reservoir_state.	None

Returns:

Type	Description
ndarray	Full state psi, shape (N [+ Na], m) or (Nt, N [+ Na], m).

Raises:

Type	Description
ValueError	If u and r have incompatible number of dimensions, or (for 3D inputs) a mismatched number of time steps.

Source code in src/models/data_driven/esn.py

def build_psi(self, u=None, r=None):
    """Assemble the full model state from physical output `u` and reservoir
    state `r`: ``concatenate([u, r, alpha])`` along the state axis, keeping only
    `u` or only `r` if `update_state`/`update_reservoir` is False, and appending
    the ensemble parameter block (`get_alpha_matrix`) if `Na > 0`.

    Parameters
    ----------
    u : np.ndarray, optional
        Physical output, shape ``(N_dim, m)`` or ``(Nt, N_dim, m)``. Defaults to
        `reservoir_to_physical(r)`.
    r : np.ndarray, optional
        Reservoir state, shape ``(N_units, m)`` or ``(Nt, N_units, m)``. Defaults
        to `reservoir_state`.

    Returns
    -------
    np.ndarray
        Full state ``psi``, shape ``(N [+ Na], m)`` or ``(Nt, N [+ Na], m)``.

    Raises
    ------
    ValueError
        If `u` and `r` have incompatible number of dimensions, or (for 3D
        inputs) a mismatched number of time steps.
    """
    if r is None:
        r = self.reservoir_state
    if u is None:
        u = self.reservoir_to_physical(r)


    if u.ndim == 2 and r.ndim == 2:
        ax_dim = 0
    elif u.ndim == 3 and r.ndim == 3:
        ax_dim = 1
        if u.shape[0] != r.shape[0]:
            raise ValueError(f'Incompatible time steps for u ({u.shape[0]}) and r ({r.shape[0]})')
    else:
        raise ValueError(f'Incompatible dimensions for u ({u.ndim}) and r ({r.ndim})')

    if self.update_state and self.update_reservoir:
        phi = np.concatenate((u, r), axis=ax_dim)
    elif self.update_state:
        phi = u
    else:
        phi = r


    if self.Na > 0:
        alph = self.get_alpha_matrix
        if u.ndim == 3:
            alph = np.tile(alph, reps=(u.shape[0], 1, 1)) # repeat for all time steps (alpha is constant in time)
        return np.concatenate((phi, alph), axis=ax_dim)
    else:
        return phi

unbuild_psi(psi=None)

Inverse of build_psi: split a full state vector into its physical (u) and reservoir (r) blocks, assuming they occupy the leading N_dim + N_units rows of psi (as build_psi lays them out when both update_state and update_reservoir are True).

Parameters:

Name	Type	Description	Default
psi	ndarray	Full state, shape (N, m) or (Nt, N, m). Defaults to current_state.	None

Returns:

Name	Type	Description
u	ndarray	Physical state, shape (N_dim, m) (or (Nt, N_dim, m)).
r	ndarray	Reservoir state, shape (N_units, m) (or (Nt, N_units, m)).

Raises:

Type	Description
AssertionError	If psi has fewer than N_units + 1 rows.

Source code in src/models/data_driven/esn.py

def unbuild_psi(self, psi=None):
    """Inverse of `build_psi`: split a full state vector into its physical
    (`u`) and reservoir (`r`) blocks, assuming they occupy the leading
    ``N_dim + N_units`` rows of `psi` (as `build_psi` lays them out when both
    `update_state` and `update_reservoir` are True).

    Parameters
    ----------
    psi : np.ndarray, optional
        Full state, shape ``(N, m)`` or ``(Nt, N, m)``. Defaults to
        `current_state`.

    Returns
    -------
    u : np.ndarray
        Physical state, shape ``(N_dim, m)`` (or ``(Nt, N_dim, m)``).
    r : np.ndarray
        Reservoir state, shape ``(N_units, m)`` (or ``(Nt, N_units, m)``).

    Raises
    ------
    AssertionError
        If `psi` has fewer than `N_units` + 1 rows.
    """
    if psi is None:
        psi = self.current_state
    if psi.ndim == 2:
        psi = np.expand_dims(psi, axis=0)
        squeeze = True
    else:
        squeeze = False

    assert psi.shape[1] > self.N_units, f"Expected psi shape (N x m) with N > {self.N_units}, got {psi.shape}"
    u = psi[:, :self.N_dim]
    r = psi[:, self.N_dim:self.N_dim+self.N_units]

    if squeeze:
        u = u.squeeze(axis=0)
        r = r.squeeze(axis=0)

    return u, r

plot_training_data(case, train_data, dt=None) staticmethod

Plot each dimension of train_data, shading the training/validation/test windows (case.t_train/t_val/t_test).

Parameters:

Name	Type	Description	Default
case	ESN_model	Instance providing t_train, t_val, t_test (and dt if not given).	required
train_data	ndarray	Data to plot, shape (L, Nt, N_dim) (or 1D/2D, reshaped accordingly).	required
dt	float	Time step for the x-axis. Defaults to case.dt.	None

Returns:

Type	Description
None	Displays the figure with matplotlib.pyplot.show.

Source code in src/models/data_driven/esn.py

@staticmethod
def plot_training_data(case, train_data, dt=None):
    """Plot each dimension of `train_data`, shading the training/validation/test
    windows (`case.t_train`/`t_val`/`t_test`).

    Parameters
    ----------
    case : ESN_model
        Instance providing `t_train`, `t_val`, `t_test` (and `dt` if not given).
    train_data : np.ndarray
        Data to plot, shape ``(L, Nt, N_dim)`` (or 1D/2D, reshaped accordingly).
    dt : float, optional
        Time step for the x-axis. Defaults to ``case.dt``.

    Returns
    -------
    None
        Displays the figure with `matplotlib.pyplot.show`.
    """
    if train_data.ndim == 1:
        train_data = train_data[np.newaxis, :, np.newaxis]
    elif train_data.ndim == 2:
        train_data = train_data[np.newaxis, :]

    L, Nt, Ndim = train_data.shape
    if dt is None:
        dt = case.dt
    t_data = np.arange(0, Nt) * dt
    nrows = min(Ndim*L, 10)


    _, axs = plt.subplots(nrows=nrows, ncols=1,
                            figsize=(8, nrows), sharex=True,
                            layout='constrained')
    if nrows * L > 1 and isinstance(axs, np.ndarray):
        axs = axs.T.flatten()
    else:
        axs = [axs]


    for l, data_l in enumerate(train_data):
        axs_dim = axs[l*Ndim:(l+1)*Ndim]

        for kk, ax in enumerate(axs_dim):

            ax.plot(t_data, data_l[:, kk], lw=1., color='k')
            ax.axvspan(0, case.t_train, facecolor='orange',
                       alpha=0.3, zorder=-100, label='Train')
            ax.axvspan(case.t_train, case.t_train + case.t_val,
                       facecolor='red', alpha=0.3, zorder=-100, label='Validation')
            ax.axvspan(case.t_train + case.t_val,
                       case.t_train + case.t_val + case.t_test, facecolor='navy',
                       alpha=0.2, zorder=-100, label='Test')

            ax.legend(ncols=1, loc='upper left', bbox_to_anchor=(1., 1.), frameon=False, title=f'L={l}, dim={kk}', fontsize='x-small', title_fontsize='small')
    axs[-1].set(xlabel='time')
    plt.show()

visualize_config()

Plot the trained read-out matrix (plot_Wout).

Returns:

Type	Description
None

Source code in src/models/data_driven/esn.py

def visualize_config(self):
    """Plot the trained read-out matrix (`plot_Wout`).

    Returns
    -------
    None
    """
    self.plot_Wout()

visualize_spatiotemporal_hist(y_hist=None, t=None, averaged=False, reference_y=1.0, reference_t=1.0, **kwargs)

Plot the physical and reservoir state history as space-time heat-maps (one row per state component vs. time), either per ensemble member or as the ensemble mean and standard deviation.

Parameters:

Name	Type	Description	Default
y_hist	ndarray	State history to plot, shape (Nt, N, m). Defaults to the last t_CR of hist.	None
t	ndarray	Time points for y_hist. Defaults to the corresponding hist_t slice.	None
averaged	bool	If True, plot the ensemble mean and standard deviation (2 rows); otherwise plot up to 10 individual state components. Default False.	False
reference_y	float	Value to normalize y_hist by. Default 1.0 (no normalization).	1.0
reference_t	float	Reference time used to normalize/label the time axis (see romda.utils.normalized_time). Default 1.0.	1.0
**kwargs		nrows : int, optional, number of state-component rows to plot when not averaged (default min(10, m)).	{}

Returns:

Type	Description
None	Displays the figure(s) in place.

Source code in src/models/data_driven/esn.py

def visualize_spatiotemporal_hist(self, y_hist=None, t=None, averaged=False,
                                  reference_y=1.0, reference_t: float = 1.0, **kwargs):
    """Plot the physical and reservoir state history as space-time heat-maps
    (one row per state component vs. time), either per ensemble member or as
    the ensemble mean and standard deviation.

    Parameters
    ----------
    y_hist : np.ndarray, optional
        State history to plot, shape ``(Nt, N, m)``. Defaults to the last
        `t_CR` of `hist`.
    t : np.ndarray, optional
        Time points for `y_hist`. Defaults to the corresponding `hist_t` slice.
    averaged : bool
        If True, plot the ensemble mean and standard deviation (2 rows);
        otherwise plot up to 10 individual state components. Default False.
    reference_y : float
        Value to normalize `y_hist` by. Default 1.0 (no normalization).
    reference_t : float
        Reference time used to normalize/label the time axis (see
        `romda.utils.normalized_time`). Default 1.0.
    **kwargs
        ``nrows`` : int, optional, number of state-component rows to plot when
        not `averaged` (default ``min(10, m)``).

    Returns
    -------
    None
        Displays the figure(s) in place.
    """
    if y_hist is None:
        n_t = int(self.t_CR // self.dt)
        y_hist = self.hist[-n_t:, :self.Nphi]

    if t is None:
        t = self.hist_t[-len(y_hist):]

    (t,), t_lbl = normalized_time(reference_t, t)
    assert t is not None
    if reference_y != 1.0:
        y_hist = y_hist / reference_y

    if y_hist.shape[1] > self.N_dim:
        y_hist_list = [y_hist[:, :self.N_dim], y_hist[:, self.N_dim:self.N_dim + self.N_units]]
        titles = ['Physical state', 'Reservoir state']
        labels = [self.state_labels[:self.N_dim], self.state_labels[self.N_dim:self.N_dim + self.N_units]]
        cmaps = ['RdBu_r', 'PRGn']
    else:
        y_hist_list = [y_hist]
        titles = ['Physical state']
        labels = [self.state_labels[:self.N_dim]]
        cmaps = ['RdBu_r']

    if not averaged:
        nrows_kw = kwargs.get('nrows', None)
        if nrows_kw is None:
            nrows = min(10, y_hist.shape[-1])
        else:
            nrows = int(nrows_kw)

        for y_hist, ttl, lbl, cmap in zip(y_hist_list, titles, labels, cmaps):
            fig, axs = plt.subplots(nrows=nrows, figsize=(10, 1.5 * nrows), sharey=True, sharex=True)
            axs_arr = np.atleast_1d(axs).ravel()
            lim = np.max(abs(y_hist))
            im = None

            for mi, ax in zip(range(nrows), axs_arr):
                im = ax.imshow(y_hist[:, :, mi].T,
                            aspect='auto', origin='lower',
                            cmap=cmap, vmin=-lim, vmax=lim,
                            extent=[t[0], t[-1], 0, y_hist.shape[1]])


            axs_arr[0].set(title=rf"ESN_model {ttl} spatiotemporal evolution. $N_\text{{units}}={self.N_units}$")
            axs_arr[-1].set(xlabel=t_lbl)
            ytx = np.arange(len(lbl))+.5
            if len(lbl) > 6:
                lbl, ytx = [zz[::len(lbl)//5] for zz in (lbl, ytx)]

            [ax.set(yticks=ytx, yticklabels=lbl) for ax in axs_arr]
            assert im is not None
            fig.colorbar(im, ax=axs_arr.tolist(), orientation='vertical', shrink=1/nrows)
    else:

        for y_hist, ttl, lbl, cmap in zip(y_hist_list, titles, labels, cmaps):
            # Averaged ensemble visualization
            y_mean_hist = np.mean(y_hist, axis=-1)

            fig, axs = plt.subplots(nrows=2, figsize=(10, 6), sharex=True)

            # Mean evolution
            lim_mean = np.max(abs(y_mean_hist))
            im0 = axs[0].imshow(y_mean_hist.T,
                                aspect='auto', origin='lower',
                                cmap=cmap, vmin=-lim_mean, vmax=lim_mean,
                                extent=[t[0], t[-1], 0, y_hist.shape[1]]
                                )
            axs[0].set(title=rf"{ttl} spatiotemporal evolution (mean and std). $N_\text{{units}}={self.N_units}$")
            fig.colorbar(im0, ax=axs[0], orientation='vertical')

            # Deviation covariance evolution

            var_ensemble = np.var(y_hist, axis=-1, ddof=1).T            # (Nt, Nx)
            var_ensemble = np.sqrt(var_ensemble)                     # Standard deviation

            lim_dev = np.max(abs(var_ensemble))
            im1 = axs[1].imshow(var_ensemble,  # Plot covariance of deviations
                                aspect='auto', origin='lower',
                                cmap='magma', vmin=0, vmax=lim_dev,
                                extent=[t[0], t[-1], 0, y_hist.shape[1]]
                                )

            fig.colorbar(im1, ax=axs[1], orientation='vertical')
            # Add ticks and labels

            ytx = np.arange(len(lbl))+.5
            if len(lbl) > 6:
                lbl, ytx = [zz[::len(lbl)//5] for zz in (lbl, ytx)]

            axs[1].set(xlabel=t_lbl)
            [ax.set(yticks=ytx, yticklabels=lbl) for ax in axs]

plot_Wout()

Visualize the trained read-out matrix Wout (overrides EchoStateNetwork.plot_Wout): a single heat-map if Wout_svd is False, or the (ensemble-averaged) SVD factors Wout_U, Wout_Sigma, Wout_Vh and their reconstruction if True.

Returns:

Type	Description
Figure

Source code in src/models/data_driven/esn.py

def plot_Wout(self):
    """Visualize the trained read-out matrix `Wout` (overrides
    `EchoStateNetwork.plot_Wout`): a single heat-map if `Wout_svd` is False, or
    the (ensemble-averaged) SVD factors `Wout_U`, `Wout_Sigma`, `Wout_Vh` and
    their reconstruction if True.

    Returns
    -------
    matplotlib.figure.Figure
    """
    if not self.Wout_svd:
        # Visualize the output matrix
        fig, ax = plt.subplots()
        im = ax.matshow(self.Wout.T, cmap="PRGn", aspect=4., vmin=-np.max(self.Wout), vmax=np.max(self.Wout))
        ax.tick_params(axis="x", bottom=True, top=False, labelbottom=True, labeltop=False)
        plt.colorbar(im, orientation='horizontal', extend='both')
        ax.set(ylabel='$N_u$', xlabel='$N_r$', title='$\\mathbf{W}_\\mathrm{out}$')

    else:
        fig, axs = plt.subplots(1, 4, figsize=(15, 15), width_ratios=[1, 1, 1, 1])
        eigs = self.Wout_Sigma
        if eigs.ndim >2:
            eigs = np.mean(eigs, axis=0)

        Wout = np.dot(self.Wout_U, np.dot(eigs, self.Wout_Vh))

        for W, ax, title in zip([Wout, self.Wout_U, eigs, self.Wout_Vh], axs,
                                ['$\\bar{\\mathbf{W_{out}}} = $', '$\\mathbf{U}$', '$\\bar{\\Sigma}$', '$\\mathbf{V}^\\mathrm{T}$']):
            cmap = 'PuOr'
            im = ax.imshow(W, cmap=cmap, vmin=-np.max(W), vmax=np.max(W))
            ax.set(title=title)
            # set the same colorbar for all the matrices
            fig.colorbar(im, ax=ax, shrink=.9, orientation='horizontal')


    return fig

romda.models.data_driven.pod_esn.POD_ESN(data, dt, plot_case=False, pdf_file=None, skip_sensor_placement=False, train_ESN=True, domain_of_measurement=None, down_sample_measurement=None, **kwargs)

Bases: ESN_model, POD

POD-projected echo state network: a POD decomposition reduces the (spatial) field to a handful of temporal coefficients, and an ESN_model is trained to forecast those coefficients in time.

Following the POD convention (see its docstring), writing \(\mathbf{Q} = \mathbf{X} - \bar{\mathbf{Q}}\) for the zero-mean data, the field is approximated as

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

with \(\boldsymbol{\Psi}\) (Psi) the orthonormal spatial modes and \(\boldsymbol{\Phi}\) (Phi) the temporal coefficients (\(\boldsymbol{\Phi} = \boldsymbol{\Psi}^\mathrm{T}\mathbf{Q}\), already scaled by the singular values Sigma). POD_ESN.__init__ runs the POD decomposition first, then trains the ESN on \(\boldsymbol{\Phi}\) (transposed to the (L, Nt, N_modes) layout ESN_model expects) to forecast the temporal coefficients forward in time; get_observables maps the ESN's closed-loop Phi forecast back to physical sensor readings (decode) or returns the modal coefficients directly if measure_modes.

Run the POD decomposition, train the ESN on the resulting temporal coefficients, and select sensor locations for observation.

Parameters:

Name	Type	Description	Default
data	ndarray	Data for the POD decomposition and ESN training, shape (Nu, Nt, Nx, Ny) or (Nt, Ndim*Nx*Ny).	required
dt	float	Time step of data; forwarded to ESN_model.__init__.	required
plot_case	bool	Whether to plot the POD modes/spectrum/reconstruction (and the ESN training process, if train_ESN). Default False.	False
pdf_file	str	If given (and plot_case), save the resulting figures to f'{pdf_file}.pdf'.	None
skip_sensor_placement	bool	If True, skip sensor placement and observe the POD modes directly (equivalent to measure_modes). Default False.	False
train_ESN	bool	Whether to train the ESN on the POD temporal coefficients. Default True.	True
domain_of_measurement	list	Sub-domain [x0, x1, y0, y1] to restrict candidate sensor locations to, if sensor_locations is not already given. Defaults to domain.	None
down_sample_measurement	int or (int, int)	Grid down-sampling factor(s) applied to the measurement domain before sensor selection.	None
**kwargs		Additional keyword arguments to configure the parent Model/ESN/POD classes, e.g. domain (physical domain of the data), grid_shape, Nq (number of sensors), sensor_locations, N_modes, t_train, t_val, N_units, etc.	{}

Source code in src/models/data_driven/pod_esn.py

def __init__(self,
             data,
             dt,
             plot_case=False,
             pdf_file=None,
             skip_sensor_placement=False,
             train_ESN=True,
             domain_of_measurement=None,
             down_sample_measurement=None,
             **kwargs):
    """Run the POD decomposition, train the ESN on the resulting temporal
    coefficients, and select sensor locations for observation.

    Parameters
    ----------
    data : np.ndarray
        Data for the POD decomposition and ESN training, shape
        ``(Nu, Nt, Nx, Ny)`` or ``(Nt, Ndim*Nx*Ny)``.
    dt : float
        Time step of `data`; forwarded to `ESN_model.__init__`.
    plot_case : bool
        Whether to plot the POD modes/spectrum/reconstruction (and the ESN
        training process, if `train_ESN`). Default False.
    pdf_file : str, optional
        If given (and `plot_case`), save the resulting figures to
        ``f'{pdf_file}.pdf'``.
    skip_sensor_placement : bool
        If True, skip sensor placement and observe the POD modes directly
        (equivalent to `measure_modes`). Default False.
    train_ESN : bool
        Whether to train the ESN on the POD temporal coefficients. Default True.
    domain_of_measurement : list, optional
        Sub-domain ``[x0, x1, y0, y1]`` to restrict candidate sensor locations
        to, if `sensor_locations` is not already given. Defaults to `domain`.
    down_sample_measurement : int or (int, int), optional
        Grid down-sampling factor(s) applied to the measurement domain before
        sensor selection.
    **kwargs
        Additional keyword arguments to configure the parent Model/ESN/POD
        classes, e.g. ``domain`` (physical domain of the data), ``grid_shape``,
        ``Nq`` (number of sensors), ``sensor_locations``, ``N_modes``,
        ``t_train``, ``t_val``, ``N_units``, etc.
    """

    for key in list(kwargs.keys()):
        if key in vars(POD_ESN):
            setattr(self, key, kwargs.pop(key))

    # __________________________ Init POD ___________________________ #
    POD.__init__(self,
                 X=data,
                 **kwargs)  # Initialize POD class and run decomposition

    # __________________________ Init ESN ___________________________ #
    # Initialize ESN to forecast the POD coefficients
    if train_ESN:
        ESN_model.__init__(self,
                           data=phi_to_esn_layout(self.Phi.copy()),
                           dt = dt,
                           plot_training=plot_case,
                           **kwargs)

    # __________________________ Select sensors ___________________________ #
    if self.measure_modes or skip_sensor_placement:
        self.Nq = self.N_modes
    elif self.sensor_locations is None:
        self.domain_of_measurement = domain_of_measurement
        self.down_sample_measurement = down_sample_measurement
        self.sensor_locations = self.define_sensors(N_sensors=self.Nq)
        self.Nq = len(self.sensor_locations)
    else:
        # If the sensors are already defined, use them
        self.Nq = len(self.sensor_locations)


    if plot_case:
        rec = self.reconstruct(X=data[:, -1])
        rec = self._to_physical_grid(rec)
        err = np.sqrt((rec - data[:, -1])**2) / np.nanmax(data[:, -1]**2)
        datasets={'Input': data[:, -1],
                  f'POD {self.N_modes} modes': rec,
                  'Error': err}

        POD_ESN.plot_case(case=self, num_modes=self.N_modes, datasets=datasets)


        if pdf_file is not None:
            self.pdf_file = pdf_file
            if isinstance(self.pdf_file, str):
                self.pdf_file = plt_pdf.PdfPages(f'{self.pdf_file}.pdf')

            figs = [plt.figure(ii) for ii in plt.get_fignums()]
            for fig in figs:
                add_pdf_page(self.pdf_file, fig_to_add=fig, close_figs=True)


    print('========= POD-ESN model complete =========')

obs_labels property

list of str: LaTeX labels for the observed quantities -- POD coefficients \(\Phi_1, \dots, \Phi_{N_\mathrm{modes}}\) if measure_modes, otherwise the sensor readings' \(u_x\)/\(u_y\) components.

state_labels property

list of str: LaTeX labels for the state vector, \(\Phi_1, \dots, \Phi_{N_\mathrm{modes}}\) (POD temporal coefficients, if update_state) followed by \(r_1, \dots, r_{N_\mathrm{units}}\) (reservoir units, if update_reservoir).

N_sensors property

int: Number of sensor locations per velocity component (\(u_x\) or \(u_y\)), i.e. Nq // 2 (each location contributes 2 observables); 0 if measure_modes.

sensor_rows property

np.ndarray or None: Rows of Psi/Q_mean corresponding to sensor_locations (cached after first access; invalidated by select_sensors). sensor_locations are raw-grid indices (var * Nx * Ny + g), whereas the POD basis rows follow the masked flat ordering -- this property maps between the two, via grid_index_to_flat_rows. None if sensor_locations is None.

domain_of_measurement property writable

list: Sub-domain [x0, x1, y0, y1] sensors may be placed in. Defaults to the full domain if never set.

down_sample_measurement property writable

list of int or None: Grid down-sampling factors [step_x, step_y] applied to the measurement grid before sensor placement. None (no down-sampling) if never set.

grid_of_measurement property

np.ndarray: Flat-grid indices (fluid cells only, all velocity components) eligible for sensor placement, i.e. domain_of_measurement down-sampled by down_sample_measurement and intersected with fluid_mask_flat.

Raises:

Type	Description
ValueError	If domain_of_measurement does not overlap the model domain.

get_POD_coefficients(Nt=1)

Read the forecasted POD temporal coefficients off the state history.

Parameters:

Name	Type	Description	Default
Nt	int	Number of trailing history steps to return. Default 1.	1

Returns:

Type	Description
ndarray	\(\boldsymbol{\Phi}\)-block of hist, shape (N_modes, m) if Nt == 1 else (Nt, N_modes, m).

Source code in src/models/data_driven/pod_esn.py

def get_POD_coefficients(self, Nt=1):
    r"""Read the forecasted POD temporal coefficients off the state history.

    Parameters
    ----------
    Nt : int
        Number of trailing history steps to return. Default 1.

    Returns
    -------
    np.ndarray
        $\boldsymbol{\Phi}$-block of `hist`, shape ``(N_modes, m)`` if
        ``Nt == 1`` else ``(Nt, N_modes, m)``.
    """
    if Nt == 1:
        Phi = self.hist[-1, :self.N_modes]
    else:
        Phi = self.hist[-Nt:, :self.N_modes]
    return Phi

get_observables(Nt=1, Phi=None, **kwargs)

Map the (forecasted) POD coefficients to observables: the coefficients themselves if measure_modes, otherwise physical-space sensor readings obtained via decode at sensor_rows.

Parameters:

Name	Type	Description	Default
Nt	int	Number of trailing history steps to return. Default 1.	1
Phi	ndarray	POD coefficients to decode, shape (N_modes, Nt, m). Defaults to get_POD_coefficients(Nt).	None
**kwargs		Unused; accepted for interface compatibility.	{}

Returns:

Type	Description
ndarray	Observables, shape (Nq, m) (or (Nt, Nq, m) for the reshaped case).

Source code in src/models/data_driven/pod_esn.py

def get_observables(self, Nt=1, Phi=None, **kwargs):
    """Map the (forecasted) POD coefficients to observables: the coefficients
    themselves if `measure_modes`, otherwise physical-space sensor readings
    obtained via `decode` at `sensor_rows`.

    Parameters
    ----------
    Nt : int
        Number of trailing history steps to return. Default 1.
    Phi : np.ndarray, optional
        POD coefficients to decode, shape ``(N_modes, Nt, m)``. Defaults to
        `get_POD_coefficients(Nt)`.
    **kwargs
        Unused; accepted for interface compatibility.

    Returns
    -------
    np.ndarray
        Observables, shape ``(Nq, m)`` (or ``(Nt, Nq, m)`` for the reshaped case).
    """
    if self.measure_modes:
        obs = self.get_POD_coefficients(Nt=Nt)
    else:
        if Phi is None:
            Phi = self.get_POD_coefficients(Nt=Nt) # Nt x N_modes x Ndim

        og_shape = Phi.shape
        reshape = Phi.ndim == 3
        if reshape:
            Phi = Phi.transpose(1, 0, 2)  # N_modes x Nt x Ndim
            Phi = Phi.reshape(self.N_modes, -1)  # N_modes x Nt*Ndim

        obs = self.decode(Z=Phi, idx=self.sensor_rows) #shape (N, Nt*Ndim)

        if reshape:
            obs = obs.reshape(self.Nq, og_shape[0], og_shape[2])  # Nq x Nt x Ndim
            obs = obs.transpose(1, 0, 2)  # Nt x Nq x Ndim


    return obs # Nt x Nq x m

reset_case(reset_POD=False, reset_ESN=False, Phi0=None, **kwargs)

Optionally rerun the POD decomposition and/or reset (retrain) the ESN.

Parameters:

Name	Type	Description	Default
reset_POD	bool	If True, rerun the POD decomposition (rerun_POD_decomposition) and force reset_ESN to True (since the modes it forecasts change). Default False.	False
reset_ESN	bool	If True, reset the EchoStateNetwork/Model state via reset_ESN (ESN_model.reset_ESN). Default False.	False
Phi0	ndarray	Passed to reset_ESN as psi0 (defaults to Phi[0]); note reset_ESN itself recomputes psi0 from the freshly trained network before resetting the model, so this value is not the final initial state.	None
**kwargs		Forwarded to rerun_POD_decomposition and/or reset_ESN.	{}

Returns:

Type	Description
None

Source code in src/models/data_driven/pod_esn.py

def reset_case(self, reset_POD=False, reset_ESN=False, Phi0=None, **kwargs):
    """Optionally rerun the POD decomposition and/or reset (retrain) the ESN.

    Parameters
    ----------
    reset_POD : bool
        If True, rerun the POD decomposition (`rerun_POD_decomposition`) and
        force `reset_ESN` to True (since the modes it forecasts change).
        Default False.
    reset_ESN : bool
        If True, reset the `EchoStateNetwork`/`Model` state via `reset_ESN`
        (`ESN_model.reset_ESN`). Default False.
    Phi0 : np.ndarray, optional
        Passed to `reset_ESN` as ``psi0`` (defaults to `Phi[0]`); note
        `reset_ESN` itself recomputes ``psi0`` from the freshly trained network
        before resetting the model, so this value is not the final initial state.
    **kwargs
        Forwarded to `rerun_POD_decomposition` and/or `reset_ESN`.

    Returns
    -------
    None
    """
    if reset_POD:
        self.rerun_POD_decomposition(**kwargs)
        reset_ESN = True  # The ESN must be reset to account for the change in POD modes

    if reset_ESN:
        if Phi0 is None:
            Phi0 = self.Phi[0]
        self.reset_ESN(psi0=Phi0, **kwargs)

select_sensors(measure_modes=False, domain_of_measurement=None, down_sample_measurement=None, N_sensors=None, qr_selection=False)

(Re)configure how the model is observed: either the raw POD coefficients (measure_modes) or physical-space point sensors placed by define_sensors.

Parameters:

Name	Type	Description	Default
measure_modes	bool	If True, observe the POD coefficients directly (Nq = N_modes, no sensors). Default False.	False
domain_of_measurement	list	Sub-domain [x0, x1, y0, y1] to restrict candidate sensor locations to. Defaults to domain.	None
down_sample_measurement	int or (int, int)	Grid down-sampling factor(s) for the measurement domain.	None
N_sensors	int	Number of sensor locations to place. Defaults to N_sensors.	None
qr_selection	bool	Whether to use QR-pivoting sensor placement (define_sensors). Default False.	False

Returns:

Type	Description
None	Sets measure_modes, sensor_locations and Nq in place.

Source code in src/models/data_driven/pod_esn.py

def select_sensors(self, measure_modes=False,
                  domain_of_measurement=None,
                  down_sample_measurement=None,
                  N_sensors=None, qr_selection=False):
    """(Re)configure how the model is observed: either the raw POD coefficients
    (`measure_modes`) or physical-space point sensors placed by `define_sensors`.

    Parameters
    ----------
    measure_modes : bool
        If True, observe the POD coefficients directly (`Nq = N_modes`, no
        sensors). Default False.
    domain_of_measurement : list, optional
        Sub-domain ``[x0, x1, y0, y1]`` to restrict candidate sensor locations
        to. Defaults to `domain`.
    down_sample_measurement : int or (int, int), optional
        Grid down-sampling factor(s) for the measurement domain.
    N_sensors : int, optional
        Number of sensor locations to place. Defaults to `N_sensors`.
    qr_selection : bool
        Whether to use QR-pivoting sensor placement (`define_sensors`).
        Default False.

    Returns
    -------
    None
        Sets `measure_modes`, `sensor_locations` and `Nq` in place.
    """
    self.measure_modes = measure_modes
    self._sensor_rows = None  # invalidate the cached Psi-row mapping
    if measure_modes:
        self.Nq = self.N_modes
        self.sensor_locations = None
    else:
        self.domain_of_measurement = domain_of_measurement
        self.down_sample_measurement = down_sample_measurement
        self.qr_selection = qr_selection
        self.sensor_locations = self.define_sensors(N_sensors=N_sensors)
        self.Nq = len(self.sensor_locations)

define_sensors(N_sensors=None, plot=False)

Choose sensor grid locations within grid_of_measurement.

If qr_selection, uses column-pivoted QR on the (physical-grid) spatial modes Psi restricted to the candidate locations, so the chosen sensors best condition the mode-reconstruction problem (a greedy, deterministic sensor-placement heuristic); otherwise picks N_sensors random locations (rng).

Parameters:

Name	Type	Description	Default
N_sensors	int	Number of sensor locations to place. Defaults to N_sensors.	None
plot	bool	Show the debug scatter of grid/measurement domain/sensors (blocks in interactive backends). Default False.	False

Returns:

Type	Description
ndarray	Flat-grid sensor indices, one block per velocity component, shape (grid_shape[0] * N_sensors,).

Source code in src/models/data_driven/pod_esn.py

def define_sensors(self, N_sensors=None, plot=False):
    """Choose sensor grid locations within `grid_of_measurement`.

    If `qr_selection`, uses column-pivoted QR on the (physical-grid) spatial
    modes `Psi` restricted to the candidate locations, so the chosen sensors
    best condition the mode-reconstruction problem (a greedy, deterministic
    sensor-placement heuristic); otherwise picks `N_sensors` random locations
    (`rng`).

    Parameters
    ----------
    N_sensors : int, optional
        Number of sensor locations to place. Defaults to `N_sensors`.
    plot : bool, optional
        Show the debug scatter of grid/measurement domain/sensors (blocks in
        interactive backends). Default False.

    Returns
    -------
    np.ndarray
        Flat-grid sensor indices, one block per velocity component, shape
        ``(grid_shape[0] * N_sensors,)``.
    """
    # Define the measurement grid

    Nu, Nx, Ny = self.grid_shape
    measure_grid_idx = np.asarray(self.grid_of_measurement)

    one_dom = measure_grid_idx[measure_grid_idx < Nx * Ny]  # only the first variable (e.g., ux) for sensor placement



    if N_sensors is None:
        N_sensors = self.N_sensors


    if self.qr_selection:

        Psi = self._to_physical_grid(self.Psi).transpose(1, 0, 2, 3)  # (r, Nu, Nx, Ny)

        Psi = np.nan_to_num(Psi, nan=0.0)
        Psi = Psi.reshape(Psi.shape[0], Nu, Nx * Ny)

        # choose one variable block for placement, e.g. variable 0
        A = Psi.reshape(Psi.shape[0], -1)  # shape (n_candidates, r)
        A = A[:, measure_grid_idx].T # shape (r, n_candidates)


        if N_sensors > A.shape[1]:
            A = np.dot(A, A.T)  # shape (n_candidates, n_candidates)

        qr_idx = sla.qr(A.T, pivoting=True)[-1]

        sensor_idx = measure_grid_idx[qr_idx[:N_sensors]]
        sensor_idx = sensor_idx.ravel() % (Nx * Ny)  # only the first variable (e.g., ux) for sensor placement

        if np.unique(sensor_idx).size < N_sensors:
            print(f'Warning: QR selection returned {np.unique(sensor_idx).size} unique sensors, less than requested {N_sensors}.')
            sensor_idx = np.unique(sensor_idx)
            extra_needed = N_sensors - sensor_idx.size
            if extra_needed > 0:
                extra_sensor_idx = measure_grid_idx[qr_idx[N_sensors:N_sensors + extra_needed]]
                extra_sensor_idx = extra_sensor_idx.ravel() % (Nx * Ny)
                sensor_idx = np.concatenate([sensor_idx, extra_sensor_idx])
    else:
        if N_sensors < len(one_dom):
            sensor_idx = np.sort(self.rng.choice(one_dom, size=N_sensors, replace=False), axis=None)
        else:
            sensor_idx = one_dom.copy()

    if N_sensors > len(measure_grid_idx):
        print(f'Requested number of sensors {N_sensors} >= grid size in domain of measurement ({len(measure_grid_idx)})')


    if plot:
        #plot the sensors against the og grid and measurement grid for debugging
        plt.figure()
        # original grid
        x_idx, y_idx = np.unravel_index(np.arange(Nx*Ny), (Nx, Ny))
        plt.scatter(x_idx, y_idx, label='Original grid', alpha=0.01
                    )
        #measurement grid
        x_idx, y_idx = np.unravel_index(measure_grid_idx[:len(measure_grid_idx)//2], (Nx, Ny))
        plt.scatter(x_idx, y_idx, label='Measurement grid')

        #sensors
        x_idx, y_idx = np.unravel_index(sensor_idx, (Nx, Ny))
        plt.scatter(x_idx, y_idx, label='Sensors')
        plt.legend()
        plt.show()

    # sensors fro all u
    sensor_idx = [sensor_idx + Nx*Ny*i for i in range(self.grid_shape[0])]

    return np.array(sensor_idx).reshape((-1,))

plot_case(case, datasets=None, num_modes=None) staticmethod

Plot the POD modes, temporal coefficients, spectrum and (if datasets is given) flow/reconstruction/error fields with sensor locations overlaid.

Parameters:

Name	Type	Description	Default
case	POD_ESN	Instance to plot.	required
datasets	dict	Named fields to pass to romda.plotting.pod.plot_flows_rms (e.g. {'Input': ..., 'Reconstruction': ..., 'Error': ...}).	None
num_modes	int	Number of POD modes to show. Defaults to case.N_modes.	None

Returns:

Type	Description
None

Source code in src/models/data_driven/pod_esn.py

@staticmethod
def plot_case(case, datasets: Optional[dict]=None, num_modes=None):
    """Plot the POD modes, temporal coefficients, spectrum and (if `datasets`
    is given) flow/reconstruction/error fields with sensor locations overlaid.

    Parameters
    ----------
    case : POD_ESN
        Instance to plot.
    datasets : dict, optional
        Named fields to pass to `romda.plotting.pod.plot_flows_rms` (e.g.
        ``{'Input': ..., 'Reconstruction': ..., 'Error': ...}``).
    num_modes : int, optional
        Number of POD modes to show. Defaults to `case.N_modes`.

    Returns
    -------
    None
    """
    from romda.plotting.pod import plot_flows_rms, plot_modes, plot_spectrum, plot_time_coefficients

    if num_modes is None:
        num_modes = case.N_modes

    plot_modes(case=case, num_modes=num_modes, cmap='viridis', n_col=2)
    plot_time_coefficients(case=case, num_modes=num_modes)
    plot_spectrum(case=case, max_mode=num_modes)


    if datasets is not None:
        display_sensors = case.sensor_locations is not None
        plot_flows_rms(case=case, datasets=datasets, display_sensors=display_sensors)

State estimation via data assimilation on a POD-ESN reduced-order model.

Joint state and parameter estimation on the POD-ESN.

romda.models.data_driven.linear_model.LinearModel(F, M_obs=None, psi0=None, dt=1.0, Q=None, **model_dict)

Bases: Model

Simple linear state-space forecast model.

\[ \boldsymbol{\psi}_{t+1} = \mathbf{F}\boldsymbol{\psi}_t + \boldsymbol{\eta}_t, \]

where \(\mathbf{F}\) is the \((N_\phi \times N_\phi)\) state transition matrix and \(\boldsymbol{\eta}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{Q})\) is optional zero-mean Gaussian process noise with covariance \(\mathbf{Q}\) (Q_noise).

Uses DiscreteIntegrator (a fixed-step map), so it only requires a time_step method instead of a continuous time_derivative.

Build a LinearModel from a fixed state transition matrix F.

Parameters:

Name	Type	Description	Default
F	ndarray	State transition matrix, shape (Nphi, Nphi) (or scalar/1D, promoted to a square 2D array via numpy.atleast_2d).	required
M_obs	ndarray	Observation (measurement) operator, shape (Nq, Nphi). Defaults to the identity (full-state observation).	None
psi0	ndarray	Initial state, shape (Nphi,) or (Nphi, m) for an ensemble. Defaults to zeros.	None
dt	float	Time step. Default 1.0.	1.0
Q	ndarray	Process noise covariance Q_noise, shape (Nphi, Nphi). Defaults to the zero matrix (no noise).	None
**model_dict		Additional Model options, forwarded to Model.__init__.	{}

Raises:

Type	Description
AssertionError	If F is not square.

Source code in src/models/data_driven/linear_model.py

def __init__(self, F, M_obs=None, psi0=None, dt=1.0, Q=None, **model_dict):
    """Build a `LinearModel` from a fixed state transition matrix `F`.

    Parameters
    ----------
    F : np.ndarray
        State transition matrix, shape ``(Nphi, Nphi)`` (or scalar/1D,
        promoted to a square 2D array via `numpy.atleast_2d`).
    M_obs : np.ndarray, optional
        Observation (measurement) operator, shape ``(Nq, Nphi)``. Defaults to
        the identity (full-state observation).
    psi0 : np.ndarray, optional
        Initial state, shape ``(Nphi,)`` or ``(Nphi, m)`` for an ensemble.
        Defaults to zeros.
    dt : float
        Time step. Default 1.0.
    Q : np.ndarray, optional
        Process noise covariance `Q_noise`, shape ``(Nphi, Nphi)``. Defaults to
        the zero matrix (no noise).
    **model_dict
        Additional `Model` options, forwarded to `Model.__init__`.

    Raises
    ------
    AssertionError
        If `F` is not square.
    """
    F = np.atleast_2d(np.array(F, dtype=float))
    Nphi = F.shape[0]
    assert F.shape == (Nphi, Nphi), f"F must be square, got {F.shape}"

    self.F = F
    self.Q_noise = np.zeros((Nphi, Nphi)) if Q is None else np.array(Q, dtype=float)
    self._has_noise = not np.allclose(self.Q_noise, 0)

    if psi0 is None:
        psi0 = np.zeros(Nphi)
    psi0 = np.atleast_1d(np.array(psi0, dtype=float))

    # Measurement operator M  (Nq x Nphi)
    if M_obs is None:
        M_obs = np.eye(Nphi)
    self._M_obs = np.atleast_2d(np.array(M_obs, dtype=float))
    self.Nq = self._M_obs.shape[0]

    super().__init__(psi0=psi0, dt=dt,
                     integrator_class=DiscreteIntegrator,
                     **model_dict)

state_labels property

list of str: LaTeX labels for the state vector, \(x_0, \dots, x_{N_\phi-1}\).

obs_labels property

list of str: LaTeX labels for the observed outputs, \(y_0, \dots, y_{N_q-1}\).

time_step(Nt)

Propagate the state forward Nt steps, \(\boldsymbol{\psi}_{k+1} = \mathbf{F}\boldsymbol{\psi}_k + \boldsymbol{\eta}_k\) with \(\boldsymbol{\eta}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{Q})\) added if Q_noise is non-zero.

Parameters:

Name	Type	Description	Default
Nt	int	Number of steps to propagate.	required

Returns:

Name	Type	Description
psi_out	ndarray	State trajectory, shape (Nt + 1, Nphi, m) (psi_out[0] is the current state).
t_out	ndarray	Corresponding time points, shape (Nt + 1,).

Source code in src/models/data_driven/linear_model.py

def time_step(self, Nt):
    r"""Propagate the state forward `Nt` steps,
    $\boldsymbol{\psi}_{k+1} = \mathbf{F}\boldsymbol{\psi}_k + \boldsymbol{\eta}_k$
    with $\boldsymbol{\eta}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{Q})$ added if
    `Q_noise` is non-zero.

    Parameters
    ----------
    Nt : int
        Number of steps to propagate.

    Returns
    -------
    psi_out : np.ndarray
        State trajectory, shape ``(Nt + 1, Nphi, m)`` (``psi_out[0]`` is the
        current state).
    t_out : np.ndarray
        Corresponding time points, shape ``(Nt + 1,)``.
    """
    psi0 = self.current_state
    t0   = self.current_time
    dt   = self.dt

    t_out = np.round(t0 + np.arange(Nt + 1) * dt, self.precision_t)

    m = psi0.shape[1] if psi0.ndim == 2 else 1
    psi_out = np.empty((Nt + 1, self.Nphi, m))
    psi_out[0] = psi0

    for k in range(1, Nt + 1):
        psi_out[k] = self.F @ psi_out[k - 1]
        if self._has_noise:
            psi_out[k] += self.rng.multivariate_normal(
                np.zeros(self.Nphi), self.Q_noise, size=m
            ).T

    # shape expected by DiscreteIntegrator: (Nt+1, Nphi, m=1)
    return psi_out, t_out

get_observables(Nt=1, **kwargs)

Map the trailing states to observables, \(\mathbf{y} = \mathbf{M}_\mathrm{obs}\boldsymbol{\psi}\).

Parameters:

Name	Type	Description	Default
Nt	int	Number of trailing time steps to return (same convention as Model.get_observables). If 1 (default), the leading Nt axis is dropped; 0 returns the full history.	1
**kwargs		Unused; accepted for interface compatibility.	{}

Returns:

Type	Description
ndarray	Observed outputs, shape (Nq, m) if Nt == 1, else (Nt, Nq, m).

Source code in src/models/data_driven/linear_model.py

def get_observables(self, Nt=1, **kwargs):
    r"""Map the trailing states to observables, $\mathbf{y} = \mathbf{M}_\mathrm{obs}\boldsymbol{\psi}$.

    Parameters
    ----------
    Nt : int
        Number of trailing time steps to return (same convention as
        `Model.get_observables`). If 1 (default), the leading ``Nt`` axis
        is dropped; 0 returns the full history.
    **kwargs
        Unused; accepted for interface compatibility.

    Returns
    -------
    np.ndarray
        Observed outputs, shape ``(Nq, m)`` if ``Nt == 1``, else
        ``(Nt, Nq, m)``.
    """
    if Nt == 1:
        return self._M_obs @ self.hist[-1, :self.Nphi, :]
    psi_hist = self.hist[-Nt:, :self.Nphi, :]              # (Nt, Nphi, m)
    return np.einsum('qp,tpm->tqm', self._M_obs, psi_hist)

---

Modal decompositions (projectors)

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)