Havok¶

Module for the Hankel alternative view of Koopman (HAVOK) analysis.

References: - S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N. Kutz, Chaos as an intermittently forced linear system, Nature Communications, 8 (2017), pp. 1-9. - S. M. Hirsh, S. M. Ichinaga, S. L. Brunton, J. N. Kutz, and B. W. Brunton, Structured time-delay models for dynamical systems with connections to frenet-serret frame, Proceedings of the Royal Society A, 477 (2021). art. 20210097.

class HAVOK(svd_rank=0, delays=10, lag=1, num_chaos=1, structured=False, lstsq=True, dmd=None)[source]

Bases: object

Hankel alternative view of Koopman (HAVOK) analysis.

Parameters:

svd_rank (int or float) – the rank for the truncation; if 0, the method computes the optimal rank and uses it for the truncation; if positive integer, the method uses the argument for the truncation; if float between 0 and 1, the rank is the number of the biggest singular values that are needed to reach the ‘energy’ specified by svd_rank; if -1, the method does not compute a truncation.
delays (int) – the number of consecutive time-shifted copies of the data to use when building Hankel matrices. Note that if examining an n-dimensional data set, this means that the resulting Hankel matrix will contain n * delays rows.
lag (int) – the number of time steps between each time-shifted copy of data in the Hankel matrix. This means that each row of the Hankel matrix will be separated by a time-step of dt * lag.
num_chaos (int) – the number of forcing terms to use in the HAVOK model.
structured (bool) – whether to perform standard HAVOK or structured HAVOK (sHAVOK). If True, sHAVOK is performed, otherwise HAVOK is performed. Note that sHAVOK cannot be performed with a BOPDMD model.
lstsq (bool) – method used for computing the HAVOK operator if a DMD method is not provided. If True, least-squares is used, otherwise the pseudo- inverse is used. This parameter is ignored if dmd is provided.
dmd (DMDBase) – DMD instance used to compute the HAVOK operator. If None, least-squares or the pseudo-inverse is used depending on lstsq.

property A

Get the matrix A in the HAVOK relationship dv/dt = Av + Bu, where v denotes the linear HAVOK embeddings and u denotes the forcing terms.

Returns:: linear dynamics matrix A.
Return type:: numpy.ndarray

property B

Get the matrix B in the HAVOK relationship dv/dt = Av + Bu, where v denotes the linear HAVOK embeddings and u denotes the forcing terms.

Returns:: forcing dynamics matrix B.
Return type:: numpy.ndarray

_compute_embeddings(forcing, time, V0)[source]: Helper function that uses the fitted HAVOK model to reconstruct the time-delay embeddings for a generic forcing term, set of times, and initial condition for the time-delay embeddings.

_embeddings_to_original(V)[source]: Helper function that uses SVD and Hankel parameter information stored in the HAVOK model to convert data in time-delay embedding space back to the space of the original input data.

static _get_index_slices(x, min_jump_dist)[source]

Helper function that, given an array x of indices at which to plot, computes and returns the beginning and ending index for each consecutive set of indices.

Example:

>>> a = np.array([2, 3, 4, 5, 10, 11, 12, 25, 26, 28])
>>> _get_index_slices(a, min_jump_dist=2)
[(2, 5), (10, 12), (25, 28)]

compute_threshold(forcing=0, p=0.01, bins=50, plot=False, plot_kwargs=None)[source]

Use the distribution of forcing terms to determine a threshold at which the absolute value of the forcing is large enough to be considered “active”. This method uses a histogram of the forcing signal values and a forcing event probability in order to estimate this threshold.

Parameters:

forcing ({numpy.ndarray, iterable} or int) – (m,) array of forcing inputs to be thresholded. Alternatively, users may provide an integer, which will be used to index the stored forcing terms. By default, the first forcing term stored will be used.
p (int or float) – desired approximate probability that a forcing event occurs. Note that p must be a float between 0.0 and 1.0, and that smaller values of p will result in larger threshold values. If p is an integer instead, p will be used to index candidate thresholds that are located at the intersection of the forcing term histogram and a fitted Gaussian distribution.
bins (int or sequence of scalars or str) – bins input to the numpy.histogram function.
plot (bool) – whether or not to plot the computed histogram of forcing values and the computed threshold. A Gaussian distribution fitted to the computed histogram is also plotted if plot=True.
plot_kwargs (dict) – optional dictionary of plot parameters. Currently, one may set the figure size, the y-axis limits, and whether or not to use a semilogy scale.

Returns:

active threshold for the absolute value of the forcing terms.

Return type:

float

dehankel(H)[source]

Given a Hankel matrix H as a 2-D numpy.ndarray, uses the delays and lag attributes to unravel the data in the Hankel matrix.

Parameters:: H (numpy.ndarray) – 2-D Hankel matrix of data.
Returns:: de-Hankeled (m,) or (n, m) array of data.
Return type:: numpy.ndarray

property delay_embeddings

Get all of the HAVOK embeddings (linear dynamics and forcing). Coordinates are stored as columns of the returned matrix. Note that this is the V matrix from the SVD of the Hankel matrix.

Returns:: matrix containing all of the HAVOK embeddings.
Return type:: numpy.ndarray

property delays

Get the number of delays used when building Hankel matrices.

Returns:: the number of Hankel matrix delays.
Return type:: int

property eigs

Get the eigenvalues of the linear HAVOK operator A.

Returns:: the eigenvalues of the operator A.
Return type:: numpy.ndarray

fit(X, t)[source]

Perform the HAVOK analysis.

Parameters:

X (numpy.ndarray or iterable) – the input snapshots.
t ({numpy.ndarray, iterable} or {int, float}) – the input time vector or uniform time-step between snapshots.

property forcing

Get the HAVOK embeddings that force the linear dynamics. Coordinates are stored as columns of the returned matrix.

Returns:: matrix containing the chaotic forcing terms.
Return type:: numpy.ndarray

hankel(X)[source]

Given a data matrix X as a 1-D or 2-D numpy.ndarray, uses the delays and lag attributes to return the data as a 2-D Hankel matrix.

Parameters:: X (numpy.ndarray) – (m,) or (n, m) array of data.
Returns:: Hankel matrix of data.
Return type:: numpy.ndarray

property ho_snapshots

Get the time-delay data matrix (i.e. the Hankel matrix).

Returns:: the matrix that contains the time-delayed data.
Return type:: numpy.ndarray

property lag

Get the lag used when building Hankel matrices.

Returns:: the number of time-steps used for Hankel matrix lag.
Return type:: int

property linear_dynamics

Get the HAVOK embeddings that are governed by linear dynamics. Coordinates are stored as columns of the returned matrix.

Returns:: matrix containing the linear HAVOK embeddings.
Return type:: numpy.ndarray

property modes

Get the U matrix from the SVD of the Hankel matrix. Note that the columns of this matrix are referred to as the eigen-time-delay modes.

Returns:: matrix containing the eigen-time-delay modes.
Return type:: numpy.ndarray

property operator

Get the full HAVOK regression model, which contains A, B, and the bad fit.

Returns:: the full HAVOK regression model.
Return type:: numpy.ndarray

plot_summary(num_plot=None, index_linear=(0, 1, 2), index_forcing=0, forcing_threshold=None, min_jump_dist=None, true_switch_indices=None, figsize=(20, 4), dpi=200, filename=None)[source]

Generate a 5-element summarizing plot that contains the following: - the time-series used to apply HAVOK - the full linear operator, which contains A, B, and the bad fit - the first linear embedding term and the first forcing term - the HAVOK embeddings, along with active forcing times - the HAVOK reconstruction of the embeddings.

Parameters:

num_plot (int) – The number of time points to plot across all subplots. By default, all available data points are plotted.
index_linear (iterable) – Tuple of indices of the linear embeddings to be plotted. May contain either 2 or 3 valid indices. The final two subplots will be plotted in 2-D or 3-D depending on the number of indices provided.
index_forcing (int) – Index of the forcing term to be plotted. Note that this index refers to indices of the forcing term itself rather than the full matrix of time-delay embeddings. Hence if 0, the first forcing term will be plotted, and so on.
forcing_threshold (float) – Threshold value at which the absolute value of the forcing signal is considered large enough to be “active”.
min_jump_dist (int) – The minimum number of indices used to separate distinct forcing events. Decreasing this parameter will lead to many short forcing events, while increasing this parameter will lead to fewer longer forcing events.
true_switch_indices (numpy.ndarray or iterable) – Optional vector that contains the indices at which true chaotic bursting occurs. If provided, true bursting times are plotted on top of the forcing term.
figsize (tuple(int, int)) – Tuple in inches defining the figure size.
dpi (int) – Figure resolution.
filename (str) – If specified, the plot is saved at filename.

predict(forcing, time, V0=0)[source]

Use a custom forcing input to make system predictions.

Parameters:

forcing (numpy.ndarray) – (m, num_chaos) array of forcing inputs.
time (numpy.ndarray) – (m,) array that contains the times that correspond with the provided forcing inputs. These will also be the times at which system predictions are computed.
V0 ({numpy.ndarray, iterable} or int) – (r - num_chaos,) array that contains the initial condition of the linear dynamics. This array should contain the linear dynamics evaluated at the first time in the time array. If not provided, this initial condition is assumed to be the initial condition stored in this HAVOK model instance. If V0 is an int, V0 is assumed to be the index of the stored linear dynamics to use as an initial condition.

Returns:

system predictions evaluated at the times in time.

Return type:

numpy.ndarray

property r

Get the number of HAVOK embeddings utilized by the HAVOK model. Note that this is essentially the integer rank truncation used.

Returns:: rank of the HAVOK model.
Return type:: int

property reconstructed_data

Get the reconstructed data.

Returns:: the matrix that contains the reconstructed snapshots.
Return type:: numpy.ndarray

property reconstructed_embeddings

Get the reconstructed time-delay embeddings.

Returns:: the matrix that contains the reconstructed embeddings.
Return type:: numpy.ndarray

property singular_vals

Get the singular value spectrum of the Hankel matrix.

Returns:: the singular values of the Hankel matrix.
Return type:: numpy.ndarray

property snapshots

Get the input data (time-series or space-flattened).

Returns:: the matrix that contains the original input data.
Return type:: numpy.ndarray

property time

Get the times of the input data.

Returns:: the vector that contains the times of the input data.
Return type:: numpy.ndarray