"""
Module for DMD plotting.
"""
import warnings
from os.path import splitext
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.patches import Patch
from pydmd import MrDMD
from .bopdmd import BOPDMD
from .hankeldmd import HankelDMD
from .havok import HAVOK
from .preprocessing import PrePostProcessingDMD
mpl.rcParams["figure.max_open_warning"] = 0
[docs]def _enforce_ratio(goal_ratio, supx, infx, supy, infy):
"""
Computes the right value of `supx,infx,supy,infy` to obtain the desired
ratio in :func:`plot_eigs`. Ratio is defined as
::
dx = supx - infx
dy = supy - infy
max(dx,dy) / min(dx,dy)
:param float goal_ratio: the desired ratio.
:param float supx: the old value of `supx`, to be adjusted.
:param float infx: the old value of `infx`, to be adjusted.
:param float supy: the old value of `supy`, to be adjusted.
:param float infy: the old value of `infy`, to be adjusted.
:return tuple: a tuple which contains the updated values of
`supx,infx,supy,infy` in this order.
"""
dx = supx - infx
if dx == 0:
dx = 1.0e-16
dy = supy - infy
if dy == 0:
dy = 1.0e-16
ratio = max(dx, dy) / min(dx, dy)
if ratio >= goal_ratio:
if dx < dy:
goal_size = dy / goal_ratio
supx += (goal_size - dx) / 2
infx -= (goal_size - dx) / 2
elif dy < dx:
goal_size = dx / goal_ratio
supy += (goal_size - dy) / 2
infy -= (goal_size - dy) / 2
return (supx, infx, supy, infy)
[docs]def _plot_limits(dmd, narrow_view):
if narrow_view:
supx = max(dmd.eigs.real) + 0.05
infx = min(dmd.eigs.real) - 0.05
supy = max(dmd.eigs.imag) + 0.05
infy = min(dmd.eigs.imag) - 0.05
return _enforce_ratio(8, supx, infx, supy, infy)
return np.max(np.ceil(np.absolute(dmd.eigs)))
[docs]def plot_eigs(
dmd,
show_axes=True,
show_unit_circle=True,
figsize=(8, 8),
title="",
narrow_view=False,
dpi=None,
filename=None,
):
"""
Plot the eigenvalues.
:param dmd: DMD instance.
:type dmd: pydmd.DMDBase
:param bool show_axes: if True, the axes will be showed in the plot.
Default is True.
:param bool show_unit_circle: if True, the circle with unitary radius
and center in the origin will be showed. Default is True.
:param tuple(int,int) figsize: tuple in inches defining the figure
size. Default is (8, 8).
:param str title: title of the plot.
:param narrow_view bool: if True, the plot will show only the smallest
rectangular area which contains all the eigenvalues, with a padding
of 0.05. Not compatible with `show_axes=True`. Default is False.
:param dpi int: If not None, the given value is passed to
``plt.figure``.
:param str filename: if specified, the plot is saved at `filename`.
"""
if isinstance(dmd, MrDMD):
raise ValueError("You should use plot_eigs_mrdmd instead")
if dmd.eigs is None:
raise ValueError(
"The eigenvalues have not been computed."
"You have to call the fit() method."
)
if dpi is not None:
plt.figure(figsize=figsize, dpi=dpi)
else:
plt.figure(figsize=figsize)
plt.title(title)
plt.gcf()
ax = plt.gca()
points = ax.plot(dmd.eigs.real, dmd.eigs.imag, "bo", label="Eigenvalues")
if narrow_view:
supx, infx, supy, infy = _plot_limits(dmd, narrow_view)
# set limits for axis
ax.set_xlim((infx, supx))
ax.set_ylim((infy, supy))
# x and y axes
if show_axes:
endx = np.min([supx, 1.0])
ax.annotate(
"",
xy=(endx, 0.0),
xytext=(np.max([infx, -1.0]), 0.0),
arrowprops=dict(arrowstyle=("->" if endx == 1.0 else "-")),
)
endy = np.min([supy, 1.0])
ax.annotate(
"",
xy=(0.0, endy),
xytext=(0.0, np.max([infy, -1.0])),
arrowprops=dict(arrowstyle=("->" if endy == 1.0 else "-")),
)
else:
# set limits for axis
limit = _plot_limits(dmd, narrow_view)
ax.set_xlim((-limit, limit))
ax.set_ylim((-limit, limit))
# x and y axes
if show_axes:
ax.annotate(
"",
xy=(np.max([limit * 0.8, 1.0]), 0.0),
xytext=(np.min([-limit * 0.8, -1.0]), 0.0),
arrowprops=dict(arrowstyle="->"),
)
ax.annotate(
"",
xy=(0.0, np.max([limit * 0.8, 1.0])),
xytext=(0.0, np.min([-limit * 0.8, -1.0])),
arrowprops=dict(arrowstyle="->"),
)
plt.ylabel("Imaginary part")
plt.xlabel("Real part")
if show_unit_circle:
unit_circle = plt.Circle(
(0.0, 0.0),
1.0,
color="green",
fill=False,
label="Unit circle",
linestyle="--",
)
ax.add_artist(unit_circle)
# Dashed grid
gridlines = ax.get_xgridlines() + ax.get_ygridlines()
for line in gridlines:
line.set_linestyle("-.")
ax.grid(True)
# legend
if show_unit_circle:
ax.add_artist(
plt.legend(
[points, unit_circle],
["Eigenvalues", "Unit circle"],
loc="best",
)
)
else:
ax.add_artist(plt.legend([points], ["Eigenvalues"], loc="best"))
ax.set_aspect("equal")
if filename:
plt.savefig(filename)
else:
plt.show()
[docs]def plot_eigs_mrdmd(
dmd,
show_axes=True,
show_unit_circle=True,
figsize=(8, 8),
title="",
level=None,
node=None,
filename=None,
):
"""
Plot the eigenvalues.
:param bool show_axes: if True, the axes will be showed in the plot.
Default is True.
:param bool show_unit_circle: if True, the circle with unitary radius
and center in the origin will be showed. Default is True.
:param tuple(int,int) figsize: tuple in inches of the figure.
:param str title: title of the plot.
:param int level: plot only the eigenvalues of specific level.
:param int node: plot only the eigenvalues of specific node.
:param str filename: if specified, the plot is saved at `filename`.
"""
if not isinstance(dmd, MrDMD):
raise ValueError(f"Expected MrDMD, found {type(dmd)}")
if dmd.eigs is None:
raise ValueError(
"The eigenvalues have not been computed."
"You have to perform the fit method."
)
if level:
peigs = dmd.partial_eigs(level=level, node=node)
else:
peigs = dmd.eigs
plt.figure(figsize=figsize)
plt.title(title)
plt.gcf()
ax = plt.gca()
if not level:
cmap = plt.get_cmap("viridis")
colors = [cmap(i) for i in np.linspace(0, 1, len(dmd.dmd_tree.levels))]
points = []
for l in dmd.dmd_tree.levels:
eigs = dmd.partial_eigs(l)
points.append(
ax.plot(eigs.real, eigs.imag, ".", color=colors[l])[0]
)
else:
points = []
points.append(
ax.plot(peigs.real, peigs.imag, "bo", label="Eigenvalues")[0]
)
# set limits for axis
limit = np.max(np.ceil(np.absolute(peigs)))
ax.set_xlim((-limit, limit))
ax.set_ylim((-limit, limit))
plt.ylabel("Imaginary part")
plt.xlabel("Real part")
if show_unit_circle:
unit_circle = plt.Circle(
(0.0, 0.0), 1.0, color="green", fill=False, linestyle="--"
)
ax.add_artist(unit_circle)
# Dashed grid
gridlines = ax.get_xgridlines() + ax.get_ygridlines()
for line in gridlines:
line.set_linestyle("-.")
ax.grid(True)
ax.set_aspect("equal")
# x and y axes
if show_axes:
ax.annotate(
"",
xy=(np.max([limit * 0.8, 1.0]), 0.0),
xytext=(np.min([-limit * 0.8, -1.0]), 0.0),
arrowprops=dict(arrowstyle="->"),
)
ax.annotate(
"",
xy=(0.0, np.max([limit * 0.8, 1.0])),
xytext=(0.0, np.min([-limit * 0.8, -1.0])),
arrowprops=dict(arrowstyle="->"),
)
# legend
if level:
labels = [f"Eigenvalues - level {level}"]
else:
labels = [f"Eigenvalues - level {i}" for i in range(len(points))]
if show_unit_circle:
points += [unit_circle]
labels += ["Unit circle"]
ax.add_artist(plt.legend(points, labels, loc="best"))
if filename:
plt.savefig(filename)
else:
plt.show()
[docs]def plot_modes_2D(
dmd,
snapshots_shape=None,
index_mode=None,
filename=None,
x=None,
y=None,
order="C",
figsize=(8, 8),
):
"""
Plot the DMD Modes.
:param dmd: DMD instance.
:type dmd: pydmd.DMDBase
:param snapshots_shape: Shape of the snapshots.
:type tuple: A tuple of ints containing the shape of a single snapshot.
:param index_mode: the index of the modes to plot. By default, all
the modes are plotted.
:type index_mode: int or sequence(int)
:param str filename: if specified, the plot is saved at `filename`.
:param numpy.ndarray x: domain abscissa.
:param numpy.ndarray y: domain ordinate
:param order: read the elements of snapshots using this index order,
and place the elements into the reshaped array using this index
order. It has to be the same used to store the snapshot. 'C' means
to read/ write the elements using C-like index order, with the last
axis index changing fastest, back to the first axis index changing
slowest. 'F' means to read / write the elements using Fortran-like
index order, with the first index changing fastest, and the last
index changing slowest. Note that the 'C' and 'F' options take no
account of the memory layout of the underlying array, and only
refer to the order of indexing. 'A' means to read / write the
elements in Fortran-like index order if a is Fortran contiguous in
memory, C-like order otherwise.
:type order: {'C', 'F', 'A'}, default 'C'.
:param tuple(int,int) figsize: tuple in inches defining the figure
size. Default is (8, 8).
"""
if dmd.modes is None:
raise ValueError(
"The modes have not been computed."
"You have to perform the fit method."
)
if snapshots_shape is None:
snapshots_shape = dmd.snapshots_shape
if x is None and y is None:
if snapshots_shape is None:
raise ValueError(
"No information about the original shape of the snapshots."
)
if len(snapshots_shape) != 2:
raise ValueError("The dimension of the input snapshots is not 2D.")
# If domain dimensions have not been passed as argument,
# use the snapshots dimensions
if x is None and y is None:
x = np.arange(snapshots_shape[0])
y = np.arange(snapshots_shape[1])
xgrid, ygrid = np.meshgrid(x, y)
if index_mode is None:
index_mode = list(range(dmd.modes.shape[1]))
elif isinstance(index_mode, int):
index_mode = [index_mode]
if filename:
basename, ext = splitext(filename)
for idx in index_mode:
fig = plt.figure(figsize=figsize)
fig.suptitle(f"DMD Mode {idx}")
real_ax = fig.add_subplot(1, 2, 1)
imag_ax = fig.add_subplot(1, 2, 2)
mode = dmd.modes.T[idx].reshape(xgrid.shape, order=order)
real = real_ax.pcolor(
xgrid,
ygrid,
mode.real,
cmap="jet",
vmin=mode.real.min(),
vmax=mode.real.max(),
)
imag = imag_ax.pcolor(
xgrid,
ygrid,
mode.imag,
vmin=mode.imag.min(),
vmax=mode.imag.max(),
)
fig.colorbar(real, ax=real_ax)
fig.colorbar(imag, ax=imag_ax)
real_ax.set_aspect("auto")
imag_ax.set_aspect("auto")
real_ax.set_title("Real")
imag_ax.set_title("Imag")
# padding between elements
plt.tight_layout(pad=2.0)
if filename:
plt.savefig(f"{basename}.{idx}{ext}")
plt.close(fig)
if not filename:
plt.show()
[docs]def plot_snapshots_2D(
dmd,
snapshots_shape=None,
index_snap=None,
filename=None,
x=None,
y=None,
order="C",
figsize=(8, 8),
):
"""
Plot the snapshots.
:param dmd: DMD instance.
:type dmd: pydmd.DMDBase
:param snapshots_shape: Shape of the snapshots.
:type tuple: A tuple of ints containing the shape of a single snapshot.
:param index_snap: the index of the snapshots to plot. By default, all
the snapshots are plotted.
:type index_snap: int or sequence(int)
:param str filename: if specified, the plot is saved at `filename`.
:param numpy.ndarray x: domain abscissa.
:param numpy.ndarray y: domain ordinate
:param order: read the elements of snapshots using this index order,
and place the elements into the reshaped array using this index
order. It has to be the same used to store the snapshot. 'C' means
to read/ write the elements using C-like index order, with the last
axis index changing fastest, back to the first axis index changing
slowest. 'F' means to read / write the elements using Fortran-like
index order, with the first index changing fastest, and the last
index changing slowest. Note that the 'C' and 'F' options take no
account of the memory layout of the underlying array, and only
refer to the order of indexing. 'A' means to read / write the
elements in Fortran-like index order if a is Fortran contiguous in
memory, C-like order otherwise.
:type order: {'C', 'F', 'A'}, default 'C'.
:param tuple(int,int) figsize: tuple in inches defining the figure
size. Default is (8, 8).
"""
if dmd.snapshots is None:
raise ValueError("Input snapshots not found.")
if snapshots_shape is None:
snapshots_shape = dmd.snapshots_shape
if x is None and y is None:
if snapshots_shape is None:
raise ValueError(
"No information about the original shape of the snapshots."
)
if len(snapshots_shape) != 2:
raise ValueError("The dimension of the snapshots is not 2D.")
# If domain dimensions have not been passed as argument,
# use the snapshots dimensions
if x is None and y is None:
x = np.arange(snapshots_shape[0])
y = np.arange(snapshots_shape[1])
xgrid, ygrid = np.meshgrid(x, y)
if index_snap is None:
index_snap = list(range(dmd.snapshots.shape[1]))
elif isinstance(index_snap, int):
index_snap = [index_snap]
if filename:
basename, ext = splitext(filename)
for idx in index_snap:
fig = plt.figure(figsize=figsize)
fig.suptitle(f"Snapshot {idx}")
snapshot = dmd.snapshots.T[idx].real.reshape(xgrid.shape, order=order)
contour = plt.pcolormesh(
xgrid,
ygrid,
snapshot,
vmin=snapshot.min(),
vmax=snapshot.max(),
)
fig.colorbar(contour)
if filename:
plt.savefig(f"{basename}.{idx}{ext}")
plt.close(fig)
if not filename:
plt.show()
[docs]def plot_summary(
dmd,
*,
x=None,
y=None,
t=None,
d=1,
continuous=False,
snapshots_shape=None,
index_modes=(0, 1, 2),
filename=None,
order="C",
figsize=(12, 8),
dpi=200,
tight_layout_kwargs=None,
main_colors=("r", "b", "g"),
mode_color="k",
mode_cmap="bwr",
dynamics_color="tab:blue",
rank_color="tab:orange",
circle_color="tab:blue",
sval_ms=8,
max_eig_ms=10,
max_sval_plot=50,
title_fontsize=14,
label_fontsize=12,
plot_semilogy=False,
flip_continuous_axes=False,
):
"""
Generate a 3 x 3 summarizing plot that contains the following components:
- the singular value spectrum of the data
- the discrete-time and continuous-time DMD eigenvalues
- the DMD modes specified by the `index_modes` parameter
- the time dynamics that correspond with each plotted mode
The number of singular values used for the DMD fit are highlighted.
All eigenvalues, modes, and dynamics are sorted according to the magnitude
of their corresponding amplitude value, i.e. their significance in the fit.
Correspondence between eigenvalues, modes, and dynamics is indicated via
color coordination.
:param dmd: fitted DMD instance.
:type dmd: pydmd.DMDBase
:param x: Points along the 1st spatial dimension where data has been
collected.
:type x: np.ndarray or iterable
:param y: Points along the 2nd spatial dimension where data has been
collected. Note that this parameter is only applicable when the data
snapshots are 2-D, which must be indicated with `snapshots_shape`.
:type y: np.ndarray or iterable
:param t: The times of data collection, or the time-step between snapshots.
Note that time information must be accurate in order to accurately
visualize eigenvalues and times of the dynamics. For non-`BOPDMD`
models, the entries of t are assumed to be uniformly-spaced, and if
not provided, TimeDict information is used. This parameter is ignored
if an instance of `BOPDMD` is provided.
:type t: {numpy.ndarray, iterable} or {int, float}
:param d: Number of delays applied to the data passed to the DMD instance.
If `d` is greater than 1, then each plotted mode will be the average
mode taken across all `d` delays.
:type d: int
:param continuous: Whether or not the eigenvalues of the given DMD instance
are continuous-time. If `False`, the eigenvalues are assumed to be the
discrete-time eigenvalues. If `True`, the eigenvalues are taken to be
the continuous-time eigenvalues. Note that `continuous` is
automatically assumed to be `True` if a `BOPDMD` model is given.
:type continuous: bool
:param snapshots_shape: Shape of the snapshots. If not provided, the shape
of the snapshots and modes is assumed to be the flattened space dim of
the snapshot data. Provide as width, height dimension.
:type snapshots_shape: iterable
:param index_modes: Indices of the modes to plot after they have been
sorted based on significance. At most three may be provided.
By default, the first three leading modes are plotted.
:type index_modes: iterable
:param filename: If specified, the plot is saved at `filename`.
:type filename: str
:param order: Read the elements of snapshots using this index order,
and place the elements into the reshaped array using this index order.
It has to be the same used to store the snapshot. "C" means to
read/write the elements using C-like index order, with the last axis
index changing fastest, back to the first axis index changing slowest.
"F" means to read/write the elements using Fortran-like index order,
with the first index changing fastest, and the last index changing
slowest. Note that the "C" and "F" options take no account of the
memory layout of the underlying array, and only refer to the order of
indexing. "A" means to read/write the elements in Fortran-like index
order if a is Fortran contiguous in memory, C-like order otherwise.
"C" is used by default.
:type order: {"C", "F", "A"}
:param figsize: Width, height in inches.
:type figsize: iterable
:param dpi: Figure resolution.
:type dpi: int
:param tight_layout_kwargs: Dictionary of `tight_layout` parameters.
:type tight_layout_kwargs: dict
:param main_colors: Colors used to denote eigenvalue, mode, dynamics
associations.
:type main_colors: iterable
:param mode_color: Color used to plot the modes, if modes are 1-D.
:type mode_color: str
:param mode_cmap: Colormap used to plot the modes, if modes are 2-D.
:type mode_cmap: str
:param dynamics_color: Color used to plot the dynamics.
:type dynamics_color: str
:param rank_color: Color used to highlight the rank of the DMD fit and
all DMD eigenvalues aside from those highlighted by `index_modes`.
:type rank_color: str
:param circle_color: Color used to plot the unit circle.
:type circle_color: str
:param sval_ms: Marker size of all singular values.
:type sval_ms: int
:param max_eig_ms: Marker size of the most prominent eigenvalue. The marker
sizes of all other eigenvalues are then scaled according to eigenvalue
significance.
:type max_eig_ms: int
:param max_sval_plot: Maximum number of singular values to plot.
:type max_sval_plot: int
:param title_fontsize: Fontsize used for subplot titles.
:type title_fontsize: int
:param label_fontsize: Fontsize used for axis labels.
:type label_fontsize: int
:param plot_semilogy: Whether or not to plot the singular values on a
semilogy plot. If `True`, a semilogy plot is used.
:type plot_semilogy: bool
:param flip_continuous_axes: Whether or not to swap the real and imaginary
axes on the continuous eigenvalues plot. If `True`, the real axis will
be vertical and the imaginary axis will be horizontal, and vice versa.
:type flip_continuous_axes: bool
"""
# This plotting method is inappropriate for plotting HAVOK results.
if isinstance(dmd, HAVOK):
raise ValueError("You should use HAVOK.plot_summary() instead.")
# Check that the DMD instance has been fitted.
if dmd.modes is None:
raise ValueError("You need to perform fit() first.")
# By default, snapshots_shape is the flattened space dimension.
if snapshots_shape is None:
snapshots_shape = (len(dmd.snapshots) // d,)
# If provided, snapshots_shape must contain 2 entires.
elif len(snapshots_shape) != 2:
raise ValueError("snapshots_shape must be None or 2-D.")
# Check the length of index_modes.
if len(index_modes) > 3:
raise ValueError("index_modes must have a length of at most 3.")
# Get the actual rank used for the DMD fit.
rank = len(dmd.eigs)
# Ensure that at least rank-many singular values will be plotted.
if rank > max_sval_plot:
raise ValueError(f"max_sval_plot must be at least {rank}.")
# Override index_modes if there are less than 3 modes available.
if rank < 3:
warnings.warn(
"Provided DMD model has less than 3 modes."
"Plotting all available modes..."
)
index_modes = np.arange(rank)
# Indices cannot go past the total number of available modes.
if np.any(np.array(index_modes) >= rank):
raise ValueError(f"Cannot view past mode {rank}.")
# Sort eigenvalues, modes, and dynamics according to amplitude magnitude.
mode_order = np.argsort(-np.abs(dmd.amplitudes))
lead_eigs = dmd.eigs[mode_order]
lead_modes = dmd.modes[:, mode_order]
lead_dynamics = dmd.dynamics[mode_order]
lead_amplitudes = np.abs(dmd.amplitudes[mode_order])
# Get time information for eigenvalue conversions.
# The decisions that we make here depend on if we're dealing
# with a BOPDMD model or any other type of DMD model...
if isinstance(dmd, BOPDMD) or (
isinstance(dmd, PrePostProcessingDMD)
and isinstance(dmd.pre_post_processed_dmd, BOPDMD)
):
# BOPDMD models store time in the time attribute.
# BOPDMD models also always compute continuous-time eigenvalues.
cont_eigs = lead_eigs
time = dmd.time
dt = dmd.time[1] - dmd.time[0]
if not np.allclose(dmd.time[1:] - dmd.time[:-1], dt):
warnings.warn(
"Time step is not uniform. "
"No discrete-time eigenvalues to plot..."
)
disc_eigs = None
else:
disc_eigs = np.exp(lead_eigs * dt)
else:
# For all other dmd models, go to the TimeDict for time information,
# that or use the user-provided time information in t if available.
num_samples = dmd.snapshots.shape[-1]
if isinstance(t, (int, float)):
time = np.arange(num_samples) * t
dt = t
elif t is not None:
time = np.squeeze(np.array(t))
dt = time[1] - time[0]
if not np.allclose(time[1:] - time[:-1], dt):
warnings.warn(
"Time step is not uniform. DMD might produce unexpected "
"results. Consider using BOP-DMD instead."
)
else:
try:
time = dmd.original_timesteps
dt = dmd.original_time["dt"]
except AttributeError:
warnings.warn(
"No time information available. Using dt = 1 and t0 = 0."
)
time = np.arange(num_samples)
dt = 1.0
if continuous:
cont_eigs = lead_eigs
disc_eigs = np.exp(lead_eigs * dt)
else:
disc_eigs = lead_eigs
cont_eigs = np.log(lead_eigs) / dt
# Get mode averages across delays if time-delay was used.
if d > 1:
lead_modes = np.average(
lead_modes.reshape(
d,
lead_modes.shape[0] // d,
lead_modes.shape[1],
),
axis=0,
)
# Compute the singular values of the data matrix.
if isinstance(dmd, HankelDMD):
# Use time-delay data matrix to compute singular values.
snp = dmd.ho_snapshots
else:
# Use input data matrix to compute singular values.
snp = dmd.snapshots
s = np.linalg.svd(snp, full_matrices=False, compute_uv=False)
# Compute the percent of data variance captured by each singular value.
s_var = s * (100 / np.sum(s))
s_var = s_var[:max_sval_plot]
# Build a list of indices of the complex conjugate pairs to highlight.
# Example: If index_modes = [idx1, idx2, idx3, idx4], such that...
# idx1 has no complex conjugate pair
# idx2 and idx3 are complex conjugates
# idx4 and idx5 are complex conjugates
# Then index_modes_cc = [(idx1, idx1), (idx2, idx3), (idx4, idx5)]
index_modes_cc = []
for idx in index_modes:
eig = cont_eigs[idx]
if eig.conj() not in cont_eigs:
index_modes_cc.append((idx, idx))
elif idx not in np.array(index_modes_cc):
index_modes_cc.append((idx, list(cont_eigs).index(eig.conj())))
other_eigs = np.setdiff1d(np.arange(rank), np.array(index_modes_cc))
# Generate the summarizing plot.
fig, (eig_axes, mode_axes, dynamics_axes) = plt.subplots(
3,
3,
figsize=figsize,
dpi=dpi,
)
# PLOT 1: Plot the singular value spectrum.
eig_axes[0].set_title("Singular Values", fontsize=title_fontsize)
eig_axes[0].set_ylabel("% variance", fontsize=label_fontsize)
s_t = np.arange(len(s_var)) + 1
eig_axes[0].plot(
s_t[:rank],
s_var[:rank],
"o",
c=rank_color,
ms=sval_ms,
mec="k",
)
eig_axes[0].plot(
s_t[rank:],
s_var[rank:],
"o",
c="gray",
ms=sval_ms,
mec="k",
)
eig_axes[0].legend(
handles=[Patch(facecolor=rank_color, label="Rank of fit")]
)
if plot_semilogy:
eig_axes[0].semilogy()
# PLOTS 2-3: Plot the eigenvalues (discrete-time and continuous-time).
# Scale marker sizes to reflect their associated amplitude.
ms_vals = max_eig_ms * np.sqrt(lead_amplitudes / lead_amplitudes[0])
# PLOT 2: Plot the discrete-time eigenvalues on the unit circle.
eig_axes[1].axvline(x=0, c="k", lw=1)
eig_axes[1].axhline(y=0, c="k", lw=1)
eig_axes[1].axis("equal")
eig_axes[1].set_title("Discrete-time Eigenvalues", fontsize=title_fontsize)
t = np.linspace(0, 2 * np.pi, 100)
eig_axes[1].plot(np.cos(t), np.sin(t), c=circle_color, ls="--")
eig_axes[1].set_xlabel(r"$Re(\lambda)$", fontsize=label_fontsize)
eig_axes[1].set_ylabel(r"$Im(\lambda)$", fontsize=label_fontsize)
# PLOT 3: Plot the continuous-time eigenvalues.
eig_axes[2].axvline(x=0, c="k", lw=1)
eig_axes[2].axhline(y=0, c="k", lw=1)
eig_axes[2].axis("equal")
eig_axes[2].set_title(
"Continuous-time Eigenvalues",
fontsize=title_fontsize,
)
if flip_continuous_axes:
eig_axes[2].set_xlabel(r"$Im(\omega)$", fontsize=label_fontsize)
eig_axes[2].set_ylabel(r"$Re(\omega)$", fontsize=label_fontsize)
eig_axes[2].invert_xaxis()
cont_eigs = 1j * cont_eigs.real + cont_eigs.imag
else:
eig_axes[2].set_xlabel(r"$Re(\omega)$", fontsize=label_fontsize)
eig_axes[2].set_ylabel(r"$Im(\omega)$", fontsize=label_fontsize)
# Now plot the eigenvalues and record the colors used for each main index.
mode_colors = {}
for ax, eigs in zip([eig_axes[1], eig_axes[2]], [disc_eigs, cont_eigs]):
if eigs is not None:
# Plot the main indices and their complex conjugate.
for i, indices in enumerate(index_modes_cc):
for idx in indices:
ax.plot(
eigs[idx].real,
eigs[idx].imag,
"o",
c=main_colors[i],
ms=ms_vals[idx],
mec="k",
)
mode_colors[idx] = main_colors[i]
# Plot all other DMD eigenvalues.
for idx in other_eigs:
ax.plot(
eigs[idx].real,
eigs[idx].imag,
"o",
c=rank_color,
ms=ms_vals[idx],
mec="k",
)
# Build the spatial grid for the mode plots.
if x is None:
x = np.arange(snapshots_shape[0])
if len(snapshots_shape) == 2:
if y is None:
y = np.arange(snapshots_shape[1])
xgrid, ygrid = np.meshgrid(x, y)
# PLOTS 4-6: Plot the DMD modes.
for i, (ax, idx) in enumerate(zip(mode_axes, index_modes)):
ax.set_title(
f"Mode {idx + 1}",
c=mode_colors[idx],
fontsize=title_fontsize,
)
if len(snapshots_shape) == 1:
# Plot modes in 1-D.
ax.plot(x, lead_modes[:, idx].real, c=mode_color)
else:
# Plot modes in 2-D.
mode = lead_modes[:, idx].reshape(xgrid.shape, order=order)
vmax = np.abs(mode.real).max()
im = ax.pcolormesh(
xgrid,
ygrid,
mode.real,
vmax=vmax,
vmin=-vmax,
cmap=mode_cmap,
)
# Align the colorbar with the plotted image.
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="3%", pad=0.05)
fig.colorbar(im, cax=cax)
# PLOTS 7-9: Plot the DMD mode dynamics.
for i, (ax, idx) in enumerate(zip(dynamics_axes, index_modes)):
dynamics_data = lead_dynamics[idx].real
ax.set_title(
"Mode Dynamics",
c=mode_colors[idx],
fontsize=title_fontsize,
)
ax.plot(time, dynamics_data, c=dynamics_color)
ax.set_xlabel("Time", fontsize=label_fontsize)
# Re-adjust ylim if dynamics oscillations are extremely small.
dynamics_range = dynamics_data.max() - dynamics_data.min()
if dynamics_range / np.abs(np.average(dynamics_data)) < 1e-4:
ax.set_ylim(np.sort([0.0, 2 * np.average(dynamics_data)]))
# Padding between elements.
if tight_layout_kwargs is None:
tight_layout_kwargs = {}
plt.tight_layout(**tight_layout_kwargs)
# Save plot if filename is provided.
if filename:
plt.savefig(filename)
plt.close(fig)
else:
plt.show()