"""
Derived module from hankeldmd.py for higher order dmd.
Reference:
- S. L Clainche, J. M. Vega, Higher Order Dynamic Mode Decomposition.
Journal on Applied Dynamical Systems, 16(2), 882-925, 2017.
"""
import warnings
import numpy as np
from .hankeldmd import HankelDMD
from .utils import compute_svd
[docs]class HODMD(HankelDMD):
"""
Higher Order Dynamic Mode Decomposition
:param svd_rank: the rank for the truncation; If 0, the method computes the
optimal rank and uses it for truncation; if positive interger, 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 truncation.
:type svd_rank: int or float
:param int tlsq_rank: rank truncation computing Total Least Square. Default
is 0, that means no truncation.
:param bool exact: flag to compute either exact DMD or projected DMD.
Default is False.
:param opt: argument to control the computation of DMD modes amplitudes.
See :class:`DMDBase`. Default is False.
:type opt: bool or int
:param rescale_mode: Scale Atilde as shown in
10.1016/j.jneumeth.2015.10.010 (section 2.4) before computing its
eigendecomposition. None means no rescaling, 'auto' means automatic
rescaling using singular values, otherwise the scaling factors.
:type rescale_mode: {'auto'} or None or numpy.ndarray
:param bool forward_backward: If True, the low-rank operator is computed
like in fbDMD (reference: https://arxiv.org/abs/1507.02264). Default is
False.
:param int d: the new order for spatial dimension of the input snapshots.
Default is 1.
:param sorted_eigs: Sort eigenvalues (and modes/dynamics accordingly) by
magnitude if `sorted_eigs='abs'`, by real part (and then by imaginary
part to break ties) if `sorted_eigs='real'`. Default: False.
:type sorted_eigs: {'real', 'abs'} or False
:param reconstruction_method: Method used to reconstruct the snapshots of
the dynamical system from the multiple versions available due to how
HODMD is conceived. If `'first'` (default) the first version
available is selected (i.e. the nearest to the 0-th row in the
augmented matrix). If `'mean'` we compute the element-wise mean. If
`reconstruction_method` is an array of float values we compute the
weighted average (for each snapshots) using the given values as weights
(the number of weights must be equal to `d`).
:param svd_rank_extra: the rank for the initial reduction of the input
data, performed before the rearrangement of the input data to the
(pseudo) Hankel matrix format; If 0, the method computes the optimal
rank and uses it for truncation; if positive interger, 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
truncation.
:type svd_rank: int or float
"""
def __init__(self, svd_rank=0, tlsq_rank=0, exact=False, opt=False,
rescale_mode=None, forward_backward=False, d=1,
sorted_eigs=False, reconstruction_method="first",
svd_rank_extra=0,):
super().__init__(
svd_rank=svd_rank,
tlsq_rank=tlsq_rank,
exact=exact,
opt=opt,
rescale_mode=rescale_mode,
forward_backward=forward_backward,
d=d,
sorted_eigs=sorted_eigs,
reconstruction_method=reconstruction_method,
)
self.svd_rank_extra = svd_rank_extra # TODO improve names
self.U_extra = None
[docs] def reconstructions_of_timeindex(self, timeindex=None):
"""
Build a collection of all the available versions of the given
`timeindex`. The indexing of time instants is the same used for
:func:`reconstructed_data`. For each time instant there are at least
one and at most `d` versions. If `timeindex` is `None` the function
returns the whole collection, for all the time instants.
:param int timeindex: The index of the time snapshot.
:return: A collection of all the available versions for the requested
time instants, represented by a matrix (or tensor).
Axes:
0. Number of time instants;
1. Copies of the snapshot;
2. Space dimension of the snapshot.
The first axis is omitted if only one single time instant is
selected, in this case the output becomes a 2D matrix.
:rtype: numpy.ndarray
"""
snapshots = super().reconstructions_of_timeindex(timeindex)
if snapshots.ndim == 2: # single time instant
snapshots = self.U_extra.dot(snapshots.T).T
elif snapshots.ndim == 3: # all time instants
snapshots = np.array(
[self.U_extra.dot(snapshot.T).T for snapshot in snapshots]
)
else:
raise RuntimeError
return snapshots
[docs] def fit(self, X):
"""
Compute the Dynamic Modes Decomposition to the input data.
:param X: the input snapshots.
:type X: numpy.ndarray or iterable
"""
org_snp, snapshots_shape = self._col_major_2darray(X)
if org_snp.shape[0] == 1:
self.U_extra, _, _ = compute_svd(org_snp, -1)
warnings.warn((
"The parameter 'svd_rank_extra={}' has been ignored because "
"the given system is a scalar function").format(
self.svd_rank_extra))
else:
self.U_extra, _, _ = compute_svd(org_snp, self.svd_rank_extra)
snp = self.U_extra.T.dot(org_snp)
super().fit(snp)
self._snapshots_shape = snapshots_shape
self._snapshots = org_snp
return self