"""
Module for the parametric Dynamic Mode Decomposition.
References:
- A Dynamic Mode Decomposition Extension for the Forecasting of Parametric
Dynamical Systems, F. Andreuzzi, N. Demo, G. Rozza, 2023, SIAM Journal on
Applied Dynamical Systems
"""
import pickle
import numpy as np
# roll by one position the shape of X. if X.shape == (a,b,c), the returned
# NumPy array's shape is (b,c,a)
def back_roll_shape(X):
"""
Roll by one position the shape of `X`. if `X.shape == (a,b,c)`, the returned
NumPy array's shape is `(b,c,a)`.
"""
return np.swapaxes(np.swapaxes(X, 0, 1), 1, 2)
def roll_shape(X):
"""
Roll by one position the shape of `X`. if `X.shape == (a,b,c)`, the returned
NumPy array's shape is `(c,a,b)`.
"""
return np.swapaxes(np.swapaxes(X, 0, 2), 1, 2)
[docs]class ParametricDMD:
"""
Implementation of the parametric Dynamic Mode Decomposition proposed in
arXiv:2110.09155v1. Both the *monolithic* and *partitioned* methods are
available, see the documentation of the parameter `dmd` for more details.
:param dmd: Instance(s) of :class:`dmdbase.DMDBase`, used by the
paramtric DMD for the prediction of future spatial modal coefficients.
If `dmd` is a `list` the *partitioned* approach is selected, in this
case the number of parameters in the training set should be equal to
the number of DMD instances provided. If `dmd` is not a list, we employ
the monolithic approach.
:type dmd: DMDBase or list
:param spatial_pod: Instance of an object usable for the generation of a
ROM of the dataset (see for instance the class
`POD <https://mathlab.github.io/EZyRB/pod.html>`_ from the Python
library `EZyRB <https://github.com/mathLab/EZyRB>`_).
:param approximation: An interpolator following the standard
learning-prediction pattern (`fit()` -> `predict()`). For some
convenient wrappers see those implemented in
`EZyRB <https://github.com/mathLab/EZyRB>`_).
:param bool light: Whether this instance should be light or not. A light
instance uses less memory since it caches a smaller number of resources.
Setting `light=True` might invalidate several properties (see also
:meth:`training_modal_coefficients`).
:param dmd_fit_args: Positional arguments to be passed to the `fit` method
of the given DMD instance.
:param dmd_fit_kwargs: Keyword arguments to be passed to the `fit` method
of the given DMD instance.
"""
def __init__(
self,
dmd,
spatial_pod,
approximation,
light=False,
dmd_fit_args=None,
dmd_fit_kwargs=None,
):
self._dmd = dmd
self._spatial_pod = spatial_pod
self._approximation = approximation
if dmd_fit_args is None:
dmd_fit_args = tuple()
if not isinstance(dmd_fit_args, (list, tuple)):
raise TypeError("Expected list, tuple or None for dmd_fit_args")
self._dmd_fit_args = dmd_fit_args
if dmd_fit_kwargs is None:
dmd_fit_kwargs = {}
if not isinstance(dmd_fit_kwargs, dict):
raise TypeError("Expected dict or None for dmd_fit_kwargs")
self._dmd_fit_kwargs = dmd_fit_kwargs
self._training_parameters = None
self._parameters = None
self._ntrain = None
self._time_instants = None
self._space_dim = None
self._light = light
self._training_modal_coefficients = None
@property
def is_partitioned(self):
"""
Return `True` if this instance is partitioned, `False` if it is
monolithic.
:type: bool
"""
return self._dmd is not None and isinstance(self._dmd, list)
@property
def _reference_dmd(self):
"""
An object used as a reference for several properties like
:func:`dmd_time` and :func:`dmd_timesteps`. If this instance is
monolithic the returned value is `self._dmd`, otherwise it is the first
item of the list `self._dmd`.
:return: The object used as a reference.
:rtype: pydmd.DMDBase
"""
if self.is_partitioned:
return self._dmd[0]
return self._dmd
@property
def dmd_time(self):
"""
The time dictionary used by the reference DMD instance (see also
:func:`_reference_dmd`). Note that when you set this attribute the
value is set only for the reference DMD (see :func:`_reference_dmd`),
however when :func:`_predict_modal_coefficients` is called the values
of all DMDs become consistent.
:getter: Return the time dictionary used by the reference DMD instance.
:setter: Set the given time dictionary in the field `dmd_time` for all
DMD instances.
:type: pydmd.dmdbase.DMDTimeDict
"""
return self._reference_dmd.dmd_time
@dmd_time.setter
def dmd_time(self, value):
self._reference_dmd.dmd_time = value
@property
def dmd_timesteps(self):
"""
The timesteps in the output of this instance, which coincides with the
timesteps in the output of the reference of this instance (see
:func:`_reference_dmd`).
:return: The timesteps in the output of this instance.
:rtype: list
"""
return self._reference_dmd.dmd_timesteps
@property
def original_time(self):
"""
The original time dictionary used by this instance, which coincides
with the original dictionary used by the reference of this instance
(see :func:`_reference_dmd`).
:return: The original time dictionary used by this instance.
:rtype: dict
"""
return self._reference_dmd.original_time
@property
def original_timesteps(self):
"""
The original timesteps in the input fed to this instance, which
coincides with the original timesteps in the input fed to the reference
of this instance (see :func:`_reference_dmd`).
:return: The original timesteps in the input fed to this instance.
:rtype: list
"""
return self._reference_dmd.original_timesteps
@property
def training_parameters(self):
"""
The original parameters passed when `self.fit` was called, represented
as a 2D array (the index of the parameter vary along the first
dimension).
:type: numpy.ndarray
"""
return self._training_parameters
def _set_training_parameters(self, params):
"""
Set the value of `self._original_parameters`, while checking that the
value provided is a 2D array.
:param numpy.ndarray: A 2D array which contains the original
parameters.
"""
if isinstance(params, list):
params = np.array(params)
if params.ndim == 1:
params = params[:, None]
if params.ndim > 2:
raise ValueError("Parameters must be stored in 2D arrays.")
self._training_parameters = params
@property
def parameters(self):
"""
The new parameters to be used in `reconstructed_data`, represented
as a 2D array (the index of the parameter vary along the first
dimension). For, instance, the following feeds a set of four 3D
parameters to `ParametricDMD`:
>>> from pydmd import ParametricDMD
>>> pdmd = ParametricDMD(...)
>>> pdmd.fit(...)
>>> p0 = [0.1, 0.2, 0.1]
>>> p1 = [0.1, 0.2, 0.3],
>>> p2 = [0.2, 0.2, 0.2],
>>> p3 = [0.1, 0.2, 0.2]
>>> pdmd.parameters = np.array([p0,p1,p2,p3])
Therefore, when we collect the results from `reconstructed_data`:
>>> result = pdmd.reconstructed_data
>>> # reconstruction corresponding to p0
>>> rec_p0 = result[0]
>>> # reconstruction corresponding to p1
>>> rec_p1 = result[1]
>>> ...
:getter: Return the current parameters.
:setter: Change the current parameters.
:type: numpy.ndarray
"""
return self._parameters if hasattr(self, "_parameters") else None
@parameters.setter
def parameters(self, value):
if isinstance(value, list):
value = np.array(value)
if value.ndim == 1:
value = value[:, None]
elif value.ndim > 2:
raise ValueError("Parameters must be stored in 2D arrays.")
self._parameters = value
def _arrange_parametric_snapshots(self, X):
"""
Arrange the given parametric snapshots (see :func:`fit` for an overview
of the shape of `X`) into a 2D matrix such that the shape is distributed
as follows:
- 0: Space;
- 1: Time/Parameter.
Time varies faster than the parameter along the columns of the matrix.
An overview of the shape of the resulting matrix:
.. math::
M = \\begin{bmatrix}
x_1(t_1,\\mu_1) & \dots & x_1(t_n,\\mu_1) & x_1(t_1,\\mu_1)
& \dots & x_1(t_{n-1},\\mu_k) & x_1(t_n,\\mu_k)\\\\
\\vdots & \\dots & \\vdots & \\vdots & \\dots & \\vdots
& \\dots\\\\
x_m(t_1,\\mu_1) & \dots & x_m(t_n,\\mu_1) & x_m(t_1,\\mu_1)
& \dots & x_m(t_{n-1},\\mu_k) & x_m(t_n,\\mu_k)
\\end{bmatrix}
:math:`x(t, \mu) \in \mathbb{R}^m` is the functon which represents the
parametric system at time :math:`t` with the parameter :math:`\\mu`.
:param X: Parametric snapshots (distribition of axes like in
:func:`fit`).
:type X: numpy.ndarray
:return: Parametric snapshots arranged in a 2D matrix like explained
above.
:rtype: numpy.ndarray
"""
# swap parameters dimension and space dimension
X = np.swapaxes(X, 0, 1)
return X.reshape((X.shape[0], -1), order="C")
def _compute_training_modal_coefficients(self, space_timemu):
"""
Compute the POD modal coefficient from the given matrix, and put
the resulting coefficients (along with their time evolution in matrix
form) into a list.
In symbols, from the given matrix :math:`X^x_{t,\mu} \in
\mathbb{R}^{m \\times nk}` we compute the modal
coefficients corresponding to its columns. At this point we have
something like this:
.. math::
\\widetilde{X}^s_{t,\mu} = \\begin{bmatrix}
\\widetilde{x}_1(t_1,\\mu_1), & \dots &
\\widetilde{x}_1(t_n,\\mu_1), &
\\widetilde{x}_1(t_1,\\mu_1), & \dots &
\\widetilde{x}_1(t_{n-1},\\mu_k), &
\\widetilde{x}_1(t_n,\\mu_k)\\\\
\\vdots & \\dots & \\vdots & \\vdots & \\dots & \\vdots &
\\dots\\\\
\\widetilde{x}_p(t_1,\\mu_1), & \dots & x_p(t_n,\\mu_1) &
\\widetilde{x}_p(t_1,\\mu_1), & \dots &
\\widetilde{x}_p(t_{n-1},\\mu_k), &
\\widetilde{x}_p(t_n,\\mu_k)
\\end{bmatrix} \in \mathbb{R}^{p \\times nk}
Detecting the sub-matrices corresponding to the time evolution of the
POD modal coefficients corresponding to a particular realization of the
system for some parameter :math:`\\mu_i`, we may rewrite this matrix as
follows:
.. math::
\\widetilde{X}^s_{t,\mu} = \\begin{bmatrix}
\\widetilde{X}_{\\mu_1}, & \dots & \\widetilde{X}_{\\mu_1}
\\end{bmatrix}
The returned list contains the matrices
:math:`\\widetilde{X}_{\\mu_i} \in \\mathbb{p \\times n}`.
:param space_timemu: A matrix containing parametric/time snapshots like
the matrix returned by :func:`_arrange_parametric_snapshots`. The
input size should be `p x nk` where `p` is the dimensionality of
the full-dimensional space, `k` is the number of training parameters
and `n` is the number of time instants used for the training.
:type space_timemu: numpy.ndarray
:return: A list of `k` matrices. Each matrix has shape `r x n` where `r`
is the dimensionality of the reduced POD space, and `n`, `k` are the
same of the parameter `space_timemu`.
:rtype: list
"""
spatial_modal_coefficients = self._spatial_pod.fit(space_timemu).reduce(
space_timemu
)
return np.split(spatial_modal_coefficients, self._ntrain, axis=1)
[docs] def fit(self, X, training_parameters):
"""
Compute the parametric Dynamic Modes Decomposition from the input data
stored in the array `X`. The shape of the parameter `X` must be
used as follows:
- 0: Training parameters;
- 1: Space;
- 2: Training time instants.
The parameter `training_parameters` contains the list of training
parameters corresponding to the training datasets in `X`. For instance,
`training_parameters[0]` is the parameter which generated the dataset
in `X[0]`. For this reason `len(training_parameters)` should be equal
to `X.shape[0]`.
:param numpy.ndarray X: Training snapshots of the parametric system,
observed for two or more parameters and in multiple time instants.
:param numpy.ndarray training_parameters: Training parameters
corresponding to the snapshots in `X`.
"""
if X.shape[0] != len(training_parameters):
raise ValueError(
"Unexpected number of snapshots for the given"
"parameters. Received {} parameters, and {} snapshots".format(
len(training_parameters), X.shape[0]
)
)
# we store these values for faster access
self._ntrain, self._space_dim, self._time_instants = X.shape
if self.is_partitioned and self._ntrain != len(self._dmd):
raise ValueError(
"Invalid number of DMD instances provided: "
"expected n_train={}, got {}".format(
self._ntrain, len(self._dmd)
)
)
# store the training parameters: they will be used in
# `reconstructed_data`
self._set_training_parameters(training_parameters)
# arrange the parametric snapshots in a convenient way to perform POD
space_timemu = self._arrange_parametric_snapshots(X)
# obtain POD modal coefficients from the training set
training_modal_coefficients = self._compute_training_modal_coefficients(
space_timemu
)
if not self._light:
self._training_modal_coefficients = np.array(
training_modal_coefficients
)
# fit DMD(s) with POD modal coefficients
if self.is_partitioned:
# partitioned parametric DMD
for dmd, data in zip(self._dmd, training_modal_coefficients):
dmd.fit(data, *self._dmd_fit_args, **self._dmd_fit_kwargs)
else:
spacemu_time = np.vstack(training_modal_coefficients)
self._dmd.fit(
spacemu_time, *self._dmd_fit_args, **self._dmd_fit_kwargs
)
# ------------------------------------------------------------
# getter properties for intermediate values of the computation
@property
def training_modal_coefficients(self):
"""
Modal coefficients of the input dataset. Since this is cached after
calls to :func:`fit` this property needs to be called after :func:`fit`,
and `light` should be set to `False` in the constructor of the class.
The tensor returned has the following shape:
- 0: Training parameters;
- 1: Dimensionality of the POD sub-space;
- 2: Time.
"""
if self._light:
raise RuntimeError(
"""Light instances do not cache the property
`training_modal_coefficients`."""
)
if self._training_modal_coefficients is None:
raise RuntimeError(
"""
Property not available now, did you call `fit()`?"""
)
return self._training_modal_coefficients
@property
def forecasted_modal_coefficients(self):
"""
Modal coefficients forecasted for the input parameters.
The tensor returned has the following shape:
- 0: Training parameters;
- 1: Dimensionality of the POD sub-space;
- 2: Time.
"""
forecasted = self._predict_modal_coefficients()
return forecasted.reshape((self._ntrain, -1, forecasted.shape[1]))
@property
def interpolated_modal_coefficients(self):
"""
Modal coefficients forecasted and then interpolated for the untested
parameters.
The tensor returned has the following shape:
- 0: Parameters;
- 1: Dimensionality of the POD sub-space;
- 2: Time.
"""
forecasted = self._predict_modal_coefficients()
return self._interpolate_missing_modal_coefficients(forecasted)
# ------------------------------------------------------------
def _predict_modal_coefficients(self):
"""
Predict future spatial modal coefficients in the time instants in
`dmd_time`.
:return: Predicted spatial modal coefficients. Shape: `rk x n` (`r`:
dimensionality of POD subspace, `k`: number of training parameters,
`n`: number of snapshots).
:rtype: numpy.ndarray
"""
if self.is_partitioned:
for dmd in self._dmd:
# we want to "bound" this DMD objects' dmd_time
dmd.dmd_time = self._reference_dmd.dmd_time
return np.vstack(
list(map(lambda dmd: dmd.reconstructed_data, self._dmd))
)
return self._dmd.reconstructed_data
def _interpolate_missing_modal_coefficients(
self, forecasted_modal_coefficients
):
"""
Interpolate spatial modal coefficients for the (untested) parameters
stored in `parameters`. The interpolation uses the interpolator
provided in the constructor of this instance.
The returned value is a 3D tensor, its shape is used as follows:
- 0: Parameters;
- 1: Reduced POD space;
- 2: Time.
:param numpy.ndarray forecasted_modal_coefficients: An array of spatial
modal coefficients for tested parameters. The shape is used like in
the matrix returned by :func:`_predict_modal_coefficients`.
:return: An array of (interpolated) spatial modal coefficients for
untested parameters.
:rtype: numpy.ndarray
"""
if self.parameters is None or len(self.parameters) == 0:
raise ValueError(
"""
Unknown parameters not found. Did you set `ParametricDMD.parameters`?"""
)
approx = self._approximation
forecasted_modal_coefficients = forecasted_modal_coefficients.reshape(
(self._ntrain, -1, forecasted_modal_coefficients.shape[1]),
order="C",
)
def interpolate_future_pod_coefficients(time_slice):
approx.fit(self.training_parameters, time_slice)
return approx.predict(self.parameters)
return np.dstack(
[
interpolate_future_pod_coefficients(time_slice)[..., None]
for time_slice in roll_shape(forecasted_modal_coefficients)
]
)
@property
def reconstructed_data(self):
"""
Get the reconstructed data, for the time instants specified in
`dmd_time`, and the parameters stored in `parameters`.
The shape of the returned data is distributed as follows:
- 0: Parameters;
- 1: Space;
- 2: Time.
:return: Snapshots predicted/interpolated using parametric DMD and the
given method for ROM.
:rtype: numpy.ndarray
"""
forecasted_modal_coefficients = self._predict_modal_coefficients()
interpolated_modal_coefficients = (
self._interpolate_missing_modal_coefficients(
forecasted_modal_coefficients
)
)
return np.apply_along_axis(
self._spatial_pod.expand, 1, interpolated_modal_coefficients
)
[docs] def save(self, fname):
"""
Save the object to `fname` using the pickle module.
:param str fname: the name of file where the reduced order model will
be saved.
Example:
>>> from pydmd import ParametricDMD
>>> pdmd = ParametricDMD(...) # Construct here the rom
>>> pdmd.fit(...)
>>> pdmd.save('pydmd.pdmd')
"""
with open(fname, "wb") as output:
pickle.dump(self, output, pickle.HIGHEST_PROTOCOL)
[docs] @staticmethod
def load(fname):
"""
Load the object from `fname` using the pickle module.
:return: The `ReducedOrderModel` loaded
Example:
>>> from pydmd import ParametricDMD
>>> pdmd = ParametricDMD.load('pydmd.pdmd')
>>> print(pdmd.reconstructed_data)
"""
with open(fname, "rb") as output:
return pickle.load(output)