"""
Module for the parametric Dynamic Mode Decomposition.
"""
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`).
"""
def __init__(self, dmd, spatial_pod, approximation, light=False):
self._dmd = dmd
self._spatial_pod = spatial_pod
self._approximation = approximation
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)
else:
spacemu_time = np.vstack(training_modal_coefficients)
self._dmd.fit(spacemu_time)
# ------------------------------------------------------------
# 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)