flowtorch.analysis

flowtorch.analysis.svd

Classes and functions wrapping around torch.linalg.svd.

class flowtorch.analysis.svd.SVD(data_matrix: Tensor, rank: int | None = None, robust: bool | dict = False)[source]

Bases: object

Compute and analyze the SVD of a data matrix.

Parameters:

  • U (pt.Tensor) – left singular vectors

  • s (pt.Tensor) – singular values

  • s_rel (pt.Tensor) – singular values normalized with their sum in percent

  • s_cum (pt.Tensor) – cumulative normalized singular values in percent

  • V (pt.Tensor) – right singular values

  • L (pt.Tensor) – low rank contribution to data matrix

  • S (pt.Tensor) – sparse contribution to data matrix

  • robust (Union[bool,dict]) – data_matrix is split into low rank and sparse contributions if True or if dictionary with options for Inexact ALM algorithm; the SVD is computed only on the low rank matrix; defaults to False

  • rank (int) – rank used for truncation

  • opt_rank (int) – optimal rank according to SVHT

  • required_memory (int) – memory required to store the truncation U, s, and V

Examples

>>> from flowtorch import DATASETS
>>> from flowtorch.data import FOAMDataloader
>>> from flowtorch.analysis import SVD
>>> loader = FOAMDataloader(DATASETS["of_cavity_ascii"])
>>> data = loader.load_snapshot("p", loader.write_times[1:])
>>> svd = SVD(data, rank=100)
>>> print(svd)
SVD of a 400x5 data matrix
Selected/optimal rank: 5/2
data type: torch.float32 (4b)
truncated SVD size: 7.9297Kb
>>> svd.s_rel
tensor([9.9969e+01, 3.0860e-02, 3.0581e-04, 7.8097e-05, 3.2241e-05])
>>> svd.s_cum
tensor([ 99.9687,  99.9996,  99.9999, 100.0000, 100.0000])
>>> svd.U.shape
torch.Size([400, 5])
>>> svd = SVD(data, rank=100, robust=True)
>>> svd.L.shape
torch.Size([400, 5])
>>> svd = SVD(data, rank=100, robust={"sparsity" : 1.0, "verbose" : True, "max_iter" : 100})
>>> svd.S.shape
torch.Size([400, 5])
property L: Tensor
property S: Tensor
property U: Tensor
property V: Tensor
property opt_rank: int
property rank: int
reconstruct(rank: int | None = None) → Tensor[source]

Reconstruct the data matrix for a given rank.

Parameters:

rank (int, optional) – rank used to compute a truncated reconstruction

Returns:

reconstruction of the input data matrix

Return type:

pt.Tensor

property required_memory: int

Compute the memory size in bytes of the truncated SVD.

Returns:

cumulative size of truncated U, s, and V tensors in bytes

Return type:

int

property robust: bool | dict
property s: Tensor
property s_cum: Tensor
property s_rel: Tensor
flowtorch.analysis.svd.inexact_alm_matrix_complection(data_matrix: Tensor, sparsity: float = 1.0, tol: float = 1e-06, max_iter: int = 100, verbose: bool = False) → Tuple[Tensor, Tensor][source]

Split a data matrix in low rank and sparse contributions.

This function implements the inexact augmented Lagrange multiplier matrix completion algorithm to solve the principal component pursuit problem. The implementation is based on the Matlab code of Isabel Scherl (link), which is in turn based on the LRSLibrary.

Parameters:

  • data_matrix (pt.Tensor) – input data matrix; snapshots must be organized as column vectors

  • sparsity (float, optional) – factor to compute Lagrangian multiplyer for sparsity (typically named lambda); lower values lead to more agressive filtering

  • tol (float, optional) – tolerance for the normalized Frobenius norm of the difference between original data and the sum of low rank and sparse contributions; defaults to 1.0e-6

  • max_iter (int, optional) – maximum number of iteration before to give up, defaults to 100

  • verbose (bool, optional) – residual is printed for every iteration if True; defaults to False

Returns:

tuple holding the low rank and the sparse matrices

Return type:

Tuple[pt.Tensor, pt.Tensor]

flowtorch.analysis.mssa

Implementation of multivariate singular spectrum analysis (MSSA).

The MSSA may be used as a preprocessing tool to reduce noise in time series data. For large state vectors, the memory requirements due to time delay embedding can be reduced by projecting the data onto a truncated POD basis before performing MSSA. This variant is referred to as projected MSSA (PMSSA). The implementation is based on the work of Ohmichi et al.

class flowtorch.analysis.mssa.MSSA(data_matrix: Tensor, window_size: int | None = None, rank: int | None = None)[source]

Bases: object

Multivariate singular spectrum analysis (MSSA).

Examples

>>> import torch as pt
>>> from flowtorch.analysis import MSSA
>>> dm = pt.rand((10, 100))
>>> mssa = MSSA(dm, rank=20)
>>> rec = mssa.reconstruction
>>> err = mssa.reconstruction_error.norm() / dm.norm()
property delay: int
property reconstruction: Tensor

Compute reconstruction of the original data matrix.

In the rank-r approximation of the Hankel matrix, multiple instances of individual snapshots are reconstructed. To reconstruct the original data matrix, averaging over the anti-diagonals is performed.

Returns:

rank-r reconstruction of the original data matrix

Return type:

pt.Tensor

property reconstruction_error: Tensor

Compute elementwise reconstruction error.

Returns:

reconstruction error with one value for each element in the data matrix

Return type:

pt.Tensor

property svd: SVD

SVD of the Hankel matrix.

Returns:

truncated SVD of the Hankel matrix

Return type:

SVD

property window_size: int
class flowtorch.analysis.mssa.PMSSA(data_matrix: Tensor, window_size: int | None = None, rank: int | None = None)[source]

Bases: MSSA

Projected multivariate singular spectrum analysis (PMSSA).

Examples

>>> import torch as pt
>>> from flowtorch.analysis import PMSSA
>>> dm = pt.rand((1000, 200))
>>> pmssa = PMSSA(dm, rank=20)
>>> rec = pmssa.reconstruction
>>> err = pmssa.reconstruction_error.norm() / dm.norm()
property reconstruction: Tensor

Compute reconstruction error in the original basis

Returns:

reconstruction of the original data matrix

Return type:

pt.Tensor

property reconstruction_error: Tensor

Compute elementwise reconstruction error.

property svd_dr: SVD

SVD used for dimensionality reduction.

Returns:

SVD used for the initial projection step in PMSSA

Return type:

SVD

flowtorch.analysis.dmd

Classes and functions to compute the dynamic mode decomposition (DMD) of a data matrix.

class flowtorch.analysis.dmd.DMD(data_matrix: Tensor, dt: float, rank: int | None = None, robust: bool | dict = False, unitary: bool = False, optimal: bool = False, tlsq: bool = False, usecols: Tensor | None = None)[source]

Bases: object

Class computing the exact DMD of a data matrix.

Examples

>>> from flowtorch import DATASETS
>>> from flowtorch.data import FOAMDataloader
>>> from flowtorch.analysis import DMD
>>> path = DATASETS["of_cavity_binary"]
>>> loader = FOAMDataloader(path)
>>> data_matrix = loader.load_snapshot("p", loader.write_times)
>>> dmd = DMD(data_matrix, dt=0.1, rank=3)
>>> dmd.frequency
tensor([0., 5., 0.])
>>> dmd.growth_rate
tensor([-2.3842e-06, -4.2345e+01, -1.8552e+01])
>>> dmd.amplitude
tensor([10.5635+0.j, -0.0616+0.j, -0.0537+0.j])
>>> dmd = DMD(data_matrix, dt=0.1, rank=3, robust=True)
>>> dmd = DMD(data_matrix, dt=0.1, rank=3, robust={"tol": 1.0e-5, "verbose" : True})
property amplitude: Tensor
property dft_properties: Tuple[float, float, float]
property dynamics: Tensor
property eigvals: Tensor
property eigvals_cont: Tensor
property eigvecs: Tensor
property frequency: Tensor
property growth_rate: Tensor
property integral_contribution: Tensor

Integral contribution of individual modes according to J. Kou et al. 2017.

DOI: https://doi.org/10.1016/j.euromechflu.2016.11.015

property modes: Tensor
property operator: Tensor
partial_reconstruction(mode_indices: Set[int]) → Tensor[source]

Reconstruct data matrix with limited number of modes.

Parameters:

mode_indices (Set[int]) – mode indices to keep

Returns:

reconstructed data matrix

Return type:

pt.Tensor

predict(initial_condition: Tensor, n_steps: int) → Tensor[source]

Predict evolution over N steps starting from used-defined initial conditions.

Parameters:

  • initial_condition (pt.Tensor) – initial state vector

  • n_steps (int) – number of steps to predict

Returns:

predicted evolution including the initial state (N+1 states are returned)

Return type:

pt.Tensor

property projection_error: Tensor

Compute the difference between Y and AX.

Returns:

projection error

Return type:

pt.Tensor

property reconstruction: Tensor

Reconstruct an approximation of the original data matrix.

Returns:

reconstructed data matrix

Return type:

pt.Tensor

property reconstruction_error: Tensor

Compute the reconstruction error.

Returns:

difference between reconstruction and data matrix

Return type:

pt.Tensor

property required_memory: int

Compute the memory size in bytes of the DMD.

Returns:

cumulative size of SVD, eigen values/vectors, and DMD modes in bytes

Return type:

int

property svd: SVD
property tlsq_error: Tuple[Tensor, Tensor]

Compute the noise in X and Y.

Returns:

noise in X and Y

Return type:

Tuple[pt.Tensor, pt.Tensor]

top_modes(n: int = 10, integral: bool = False, f_min: float = -inf, f_max: float = inf) → Tensor[source]

Get the indices of the first n most important modes.

Note that the conjugate complex modes for real data matrices are not filtered out by default. However, by setting the lower frequency threshold to a positive number, only modes with positive imaginary part are considered.

Parameters:

  • n (int) – number of indices to return; defaults to 10

  • integral (bool, optional) – if True, the modes are sorted according to their integral contribution; defaults to False

  • f_min (float, optional) – consider only modes with a frequency larger or equal to f_min; defaults to -inf

  • f_max (float, optional) – consider only modes with a frequency smaller than f_max; defaults to -inf

Returns:

indices of top n modes sorted by amplitude or integral contribution

Return type:

pt.Tensor

flowtorch.analysis.hodmd

Implementation of the higher-order DMD (HODMD).

class flowtorch.analysis.hodmd.HODMD(data_matrix: Tensor, dt: float, delay: int | None = None, rank_dr: int | None = None, **dmd_options: dict)[source]

Bases: DMD

Higher-order dynamic mode decomposition (HODMD).

For the theoretical background, refer to Clainche and Vega (link). The HODMD wraps around the standard DMD by adding an initial dimensionality reduction step and an enrichment of data matrix with delays. To reconstruct snapshots and modes in the original space, a few properties of the base class are overwritten. .. _link: https://doi.org/10.1137/15M1054924

Examples

>>> from flowtorch.analysis import HODMD
>>> dmd = HODMD(data_matrix, dt)
set time delay explicitly to 5 time levels
>>> dmd = HODMD(data_matrix, dt, delay=5)
set the rank for the initial dimensionality reduction
>>> dmd = HODMD(data_matrix, dt, delay=5, rank_dr=100)
use optimal mode coefficients
>>> dmd = HODMD(data_matrix, dt, delay=5, rank_dr=100, optimal=True)
property delay: int
property dynamics: Tensor

Get mode dynamics for the original data matrix.

Returns:

mode dynamics for the original snapshot sequence

Return type:

pt.Tensor

property modes: Tensor

Get DMD modes in the input space.

As suggested by Clainche and Vega, only the first set of modes corresponding to the first r rows of the reduced DMD modes are kept (r is the dimension after the initial dimensionality reduction).

Returns:

DMD modes in the input space

Return type:

pt.Tensor

predict(initial_condition: Tensor, n_steps: int) → Tensor[source]

Predict evolution over N steps starting from used-defined initial conditions.

The prediction is performed as follows: 1) the initial conditions are projected on the first r POD modes 2) the time delay embedding is computed in the reduced space 3) the prediction is computed in the reduced space 4) the first r rows of the prediction are reconstructed

Parameters:

  • initial_condition (pt.Tensor) – initial sequence of state vectors without time delay embedding; for d delays and a state vector of size M, the initial conditions should be given as a M x d matrix sorted such that the most recent state is contained in the last column

  • n_steps (int) – number of future steps to predict

Returns:

predicted states; due to the structure of the linear operator some states are predicted multiple times; only the first M rows of the reconstruction are returned

Return type:

pt.Tensor

property projection_error: Tensor

Compute the difference between Y and AX.

Returns:

projection error

Return type:

pt.Tensor

property reconstruction: Tensor

Compute reconstruction of original data matrix.

Due to the time delays, some states are contained multiple times in a reconstruction. Only the first occurrence of a state is kept when looping over the rows in steps of size r (rank).

Returns:

reconstruction of original data matrix

Return type:

pt.Tensor

property reconstruction_error: Tensor

Compute the point-wise reconstruction error.

Due to the time delay, not all snapshots from the input data matrix are reconstructed, so the error is only computed for the available reconstructed snapshots

Returns:

reconstruction error

Return type:

pt.Tensor

property svd_dr: SVD
property tlsq_error: Tensor

Compute the noise in X and Y.

Returns:

noise in X and Y

Return type:

Tuple[pt.Tensor, pt.Tensor]

flowtorch.analysis.optdmd

Optimized DMD via gradient descent and backpropagation.

class flowtorch.analysis.optdmd.EarlyStopping(patience: int = 40, min_delta: float = 0.0, checkpoint: str | None = None, model: Module | None = None)[source]

Bases: object

Provide stopping control for iterative optimization tasks.

class flowtorch.analysis.optdmd.OptDMD(*dmd_args, **dmd_kwargs)[source]

Bases: Module

Optimized DMD based on backpropagation and gradient descent.

For a detailed description this DMD variant refer to Weiner and Semaan (2023).

Examples

>>> import torch as pt
>>> from flowtorch.analysis import OptDMD
>>> dm = pt.rand((200, 100))
>>> dmd = DMD(dm, dt=1.0)
>>> dmd.train(stopping_options={'patience' : 80, 'checkpoint' : '/tmp/best_model.pt'})
>>> dmd.load_state_dict(pt.load('/tmp/best_model.pt'))
>>> train_loss = dmd.log['train_loss']
>>> val_loss = dmd.log['val_loss']
property amplitude: Tensor
property dmd_init: DMD
property dynamics: Tensor
property eigvals: Tensor
property eigvals_cont: Tensor
property eigvecs: Tensor
forward(time_indices: Tensor) → Tensor[source]

Defines the computation performed at every call.

Should be overridden by all subclasses.

Note

Although the recipe for forward pass needs to be defined within this function, one should call the Module instance afterwards instead of this since the former takes care of running the registered hooks while the latter silently ignores them.

property frequency: Tensor
property growth_rate: Tensor
property integral_contribution: Tensor
property log: dict
property modes: Tensor
partial_reconstruction(mode_indices: Set[int]) → Tensor[source]
predict(initial_condition: Tensor, n_steps: int) → Tensor[source]

Predict evolution over N steps starting from used-defined initial conditions.

Parameters:

  • initial_condition (pt.Tensor) – initial state vector

  • n_steps (int) – number of steps to predict

Returns:

predicted evolution including the initial state (N+1 states are returned)

Return type:

pt.Tensor

property reconstruction: Tensor
property reconstruction_error: Tensor
top_modes(n: int = 10, integral: bool = False, f_min: float = -inf, f_max: float = inf) → Tensor[source]
train(epochs: int = 1000, lr: float = 0.001, batch_size: int | None = None, train_size: int | float = 0.75, val_size: int | float = 0.25, loss_function: ~typing.Callable = <function l2_loss>, scheduler_options: dict = {}, stopping_options: dict = {}, loss_key: str = 'val_loss')[source]

Optimize modes and dynamics based on gradient descent.

Parameters:

  • epochs (int, optional) – number of training iterations, defaults to 1000

  • lr (float, optional) – initial learning rate, defaults to 1.0e-3

  • batch_size (int, optional) – batch size for batch training, defaults to None

  • train_size (Union[int, float], optional) – fraction or number of snapshots to use for training; defaults to 0.75

  • val_size (Union[int, float], optional) – fraction or number of snapshots to use for validation; defaults to 0.25

  • loss_function (Callable, optional) – user-defined loss function, e.g., to add sparsity promotion, defaults to l2_loss

  • scheduler_options (dict, optional) – options passed to learning rate scheduler; refer to PyTorch’s ReduceLROnPlateau documentation; defaults to {}

  • stopping_options (dict, optional) – options to modify early stopping behavior; refer to the EarlyStopping class; defaults to {}

  • loss_key (str, optional) – key of loss value based on which the learning rate schedule and early stopping criteria are evaluated; can be ‘train_loss’, ‘val_loss’, or ‘full_loss’; defaults to “val_loss”

flowtorch.analysis.optdmd.l2_loss(label: Tensor, prediction: Tensor, eigvecs: Tensor, eigvals: Tensor) → Tensor[source]

Compute the L2 norm of the prediction error.

Note: in contrast to the default norm function in PyTorch, the norm is normalized by the number of elements in the prediction tensor.

Note: the argument list is more generic than necessary for this loss function; the eigenvectors and eigenvalues may be used for regularization or to enforce physical constraints in more advanced loss functions.

Parameters:

  • label (pt.Tensor) – ground truth

  • prediction (pt.Tensor) – predictions

Returns:

L2 norm normalized by number of elements

Return type:

pt.Tensor

flowtorch.analysis.hooptdmd

Higher-order optimized DMD version.

class flowtorch.analysis.hooptdmd.HOOptDMD(data_matrix: Tensor, dt: float, delay: int = 1, rank_dr: int | None = None, **dmd_options: dict)[source]

Bases: OptDMD

Optimized DMD extended with data projection and time delay embedding.

Parameters:

OptDMD (flowtorch.analysis.OptDMD) – Optimized DMD base class

property dynamics: Tensor
property modes
partial_reconstruction(mode_indices: Set[int]) → Tensor[source]
predict(initial_condition: Tensor, n_steps: int) → Tensor[source]

Predict evolution over N steps starting from used-defined initial conditions.

The prediction is performed as follows: 1) the initial conditions are projected on the first r POD modes 2) the time delay embedding is computed in the reduced space 3) the prediction is computed in the reduced space 4) the first r rows of the prediction are reconstructed

Parameters:

  • initial_condition (pt.Tensor) – initial sequence of state vectors without time delay embedding; for d delays and a state vector of size M, the initial conditions should be given as a M x d matrix sorted such that the most recent state is contained in the last column

  • n_steps (int) – number of future steps to predict

Returns:

predicted states; due to the structure of the linear operator some states are predicted multiple times; only the first M rows of the reconstruction are returned

Return type:

pt.Tensor

property reconstruction: Tensor
property reconstruction_error: Tensor
property svd_dr: SVD

flowtorch.analysis.psp_explorer

Module with classes and functions to explore and analyse iPSP data.

class flowtorch.analysis.psp_explorer.PSPExplorer(file_path: str)[source]

Bases: object

Explore iPSP data interactively.

interact(zone: str, field_name: str, times: list, mask: bool = True, width: int = 1024, every: int = 5, cmin: float = -2, cmax: float = 0.5) → Figure[source]

Show selected snapshots of an iPSP file as interactive surface plot.

Parameters:

  • zone (str) – name of the zone to display

  • field_name (str) – name of the field to display

  • times (list) – list of times instances to show

  • mask (bool, optional) – if True, the data is masked using the mask provided in the iPSP dataset, defaults to True

  • width (int, optional) – width of the plot, defaults to 1024

  • every (int, optional) – plot only every nth point, defaults to 5

  • cmin (float, optional) – lower bound for color range, defaults to -2

  • cmax (float, optional) – upper bound for color range, defaults to 0.5

Returns:

interactive iPSP surface plot

Return type:

go.Figure

property loader: PSPDataloader

Get the PSPdataloader instance used to access data.

This propertie allows accessing the loader’s metadata via the explorer without having to create a new loader instance.

Returns:

PSP dataloader instance

Return type:

PSPDataloader

mean(zone: str, field_name: str, times: list, mask: bool = True, width: int = 1024, every=5, cmin: float = -2, cmax: float = 0.5) → Figure[source]

Show temporal mean of an iPSP file as interactive surface plot.

Parameters:

  • zone (str) – name of the zone to display

  • field_name (str) – name of the field to display

  • times (list) – snapshot times over which to compute the mean

  • mask (bool, optional) – if True, the data is masked using the mask provided in the iPSP dataset, defaults to True

  • width (int, optional) – width of the plot, defaults to 1024

  • every (int, optional) – use only every nth point for plotting, defaults to 5

  • cmin (float, optional) – lower bound for color range, defaults to -2

  • cmax (float, optional) – upper bound for color range, defaults to 0.5

Returns:

interactive iPSP surface plot showing the temporal mean

Return type:

go.Figure

std(zone: str, field_name: str, times: list, mask: bool = True, width: int = 1024, every=5, cmin: float = 0, cmax: float = 0.5) → Figure[source]

Show temporal standard deviation as interactive surface plot.

Parameters:

  • zone (str) – name of the zone to display

  • field_name (str) – name of the field to display

  • times (list) – snapshot times over which to compute the standard deviation

  • mask (bool, optional) – if True, the data is masked using the mask provided in the iPSP dataset, defaults to True

  • width (int, optional) – width of the plot, defaults to 1024

  • every (int, optional) – use only every nth point for plotting, defaults to 5

  • cmin (float, optional) – lower bound for color range, defaults to -2

  • cmax (float, optional) – upper bound for color range, defaults to 0.5

Returns:

interactive iPSP surface plot

Return type:

go.Figure