Skip to content

Commit dcd10cd

Browse files
mtezzelemtezzele
authored
Merge pull request PyDMD#329 from sichinaga/master
Added PiDMD
2 parents 63e5c32 + dc02932 commit dcd10cd

File tree

6 files changed

+895
-2
lines changed

6 files changed

+895
-2
lines changed

Diff for: ‎docs/source/pidmd.rst

+14
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
Physics-informed DMD
2+
====================
3+
4+
.. currentmodule:: pydmd.pidmd
5+
6+
.. automodule:: pydmd.pidmd
7+
8+
.. autoclass:: PiDMD
9+
:members:
10+
:private-members:
11+
:undoc-members:
12+
:show-inheritance:
13+
:noindex:
14+

Diff for: ‎pydmd/__init__.py

+2-1
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
PyDMD init
33
"""
44
__all__ = ['dmdbase', 'dmd', 'fbdmd', 'mrdmd', 'cdmd', 'hodmd', 'dmdc',
5-
'optdmd', 'hankeldmd', 'rdmd', 'havok', 'bopdmd']
5+
'optdmd', 'hankeldmd', 'rdmd', 'havok', 'bopdmd', 'pidmd']
66

77

88
from .meta import *
@@ -22,3 +22,4 @@
2222
from .rdmd import RDMD
2323
from .havok import HAVOK
2424
from .bopdmd import BOPDMD
25+
from .pidmd import PiDMD

Diff for: ‎pydmd/pidmd.py

+322
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,322 @@
1+
"""
2+
Derived module from dmd.py for Physics-informed DMD.
3+
4+
References:
5+
- Peter J. Baddoo, Benjamin Herrmann, Beverley J. McKeon, J. Nathan Kutz, and
6+
Steven L. Brunton. Physics-informed dynamic mode decomposition (pidmd). 2021.
7+
arXiv:2112.04307.
8+
"""
9+
import numpy as np
10+
11+
from .dmd import DMD
12+
from .dmdoperator import DMDOperator
13+
from .utils import compute_svd
14+
from .pidmd_utils import (
15+
compute_unitary,
16+
compute_uppertriangular,
17+
compute_diagonal,
18+
compute_symmetric,
19+
compute_toeplitz,
20+
compute_circulant,
21+
compute_symtridiagonal,
22+
compute_BCCB,
23+
compute_BC,
24+
)
25+
26+
class PiDMDOperator(DMDOperator):
27+
"""
28+
DMD operator for Physics-informed DMD.
29+
30+
:param manifold: the matrix manifold for the full DMD operator A.
31+
:type manifold: str
32+
:param manifold_opt: option used to specify certain manifolds. If manifold
33+
is "diagonal", manifold_opt may be used to specify the width of the
34+
diagonal of A. If manifold starts with "BC", manifold_opt must be a
35+
2D tuple that specifies the desired dimensions of the blocks of A.
36+
:type manifold_opt: int, tuple(int,int), or numpy.ndarray
37+
:param compute_A: Flag that determines whether or not to compute the full
38+
Koopman operator A.
39+
:type compute_A: bool
40+
:param svd_rank: the rank for the truncation; If 0, the method computes the
41+
optimal rank and uses it for truncation; if positive integer, the
42+
method uses the argument for the truncation; if float between 0 and 1,
43+
the rank is the number of the biggest singular values that are needed
44+
to reach the 'energy' specified by `svd_rank`; if -1, the method does
45+
not compute truncation.
46+
:type svd_rank: int or float
47+
"""
48+
def __init__(
49+
self,
50+
manifold,
51+
manifold_opt,
52+
compute_A,
53+
svd_rank,
54+
):
55+
self._manifold = manifold
56+
self._manifold_opt = manifold_opt
57+
self._compute_A = compute_A
58+
self._svd_rank = svd_rank
59+
60+
self._A = None
61+
self._Atilde = None
62+
self._eigenvalues = None
63+
self._eigenvectors = None
64+
self._modes = None
65+
66+
67+
@property
68+
def A(self):
69+
"""
70+
Get the full Koopman operator A.
71+
72+
:return: the full Koopman operator A.
73+
:rtype: numpy.ndarray
74+
"""
75+
if not self._compute_A:
76+
msg = "A not computed during fit. " \
77+
"Set parameter compute_A = True to compute A."
78+
raise ValueError(msg)
79+
if self._A is None:
80+
raise ValueError("You need to call fit before")
81+
return self._A
82+
83+
84+
def _check_compute_A(self):
85+
"""
86+
Helper method that checks that compute_A is True.
87+
Throws an error if compute_A is False.
88+
"""
89+
if not self._compute_A:
90+
raise ValueError(
91+
"A must be computed for the chosen manifold."
92+
"Set compute_A = True to compute A."
93+
)
94+
95+
96+
def _compute_procrustes(self, X, Y):
97+
"""
98+
Private method that computes the best-fit linear operator A in the
99+
relationship Y = AX such that A is restricted to the family of matrices
100+
defined by the given manifold (and manifold option if applicable).
101+
Computes and returns a dictionary that contains either...
102+
(1) the reduced operator "atilde",
103+
(2) the full operator "A", or
104+
(3) the "eigenvalues" and eigenvectors of A, referred to as "modes",
105+
depending on the chosen manifold and the compute_A parameter.
106+
"""
107+
if self._manifold == "unitary":
108+
result_dict = compute_unitary(X, Y, self._svd_rank)
109+
elif self._manifold == "uppertriangular":
110+
self._check_compute_A()
111+
result_dict = compute_uppertriangular(X, Y)
112+
elif self._manifold == "lowertriangular":
113+
self._check_compute_A()
114+
A_rot = compute_uppertriangular(np.flipud(X), np.flipud(Y))["A"]
115+
result_dict = {"A": np.rot90(A_rot, 2)}
116+
elif self._manifold == "diagonal":
117+
result_dict = compute_diagonal(X, Y,
118+
self._svd_rank,
119+
self._manifold_opt,
120+
self._compute_A)
121+
elif self._manifold == "symmetric":
122+
result_dict = compute_symmetric(X, Y, self._svd_rank)
123+
elif self._manifold == "skewsymmetric":
124+
result_dict = compute_symmetric(X, Y, self._svd_rank,
125+
skew_symmetric=True)
126+
elif self._manifold == "toeplitz":
127+
self._check_compute_A()
128+
result_dict = compute_toeplitz(X, Y)
129+
elif self._manifold == "hankel":
130+
self._check_compute_A()
131+
result_dict = compute_toeplitz(X, Y, flipped=True)
132+
elif self._manifold in [
133+
"circulant",
134+
"circulant_unitary",
135+
"circulant_symmetric",
136+
"circulant_skewsymmetric"
137+
]:
138+
circulant_opt = self._manifold.replace("circulant_", "")
139+
result_dict = compute_circulant(X, Y,
140+
circulant_opt,
141+
self._svd_rank,
142+
self._compute_A)
143+
elif self._manifold == "symmetric_tridiagonal":
144+
result_dict = compute_symtridiagonal(X, Y,
145+
self._svd_rank,
146+
self._compute_A)
147+
elif self._manifold in [
148+
"BC",
149+
"BCTB",
150+
"BCCB",
151+
"BCCBunitary",
152+
"BCCBsymmetric",
153+
"BCCBskewsymmetric",
154+
]:
155+
# Specify the shape of the blocks in the output matrix A.
156+
self._check_compute_A()
157+
if self._manifold_opt is None:
158+
raise ValueError("manifold_opt must be specified.")
159+
if (not isinstance(self._manifold_opt, tuple)
160+
or len(self._manifold_opt) != 2):
161+
raise ValueError("manifold_opt is not in an allowable format.")
162+
block_shape = np.array(self._manifold_opt)
163+
164+
if self._manifold.startswith("BCCB"):
165+
bccb_opt = self._manifold.replace("BCCB", "")
166+
result_dict = compute_BCCB(X, Y,
167+
block_shape,
168+
bccb_opt,
169+
self._svd_rank)
170+
elif self._manifold == "BC":
171+
result_dict = compute_BC(X, Y, block_shape)
172+
else:
173+
result_dict = compute_BC(X, Y, block_shape,
174+
tridiagonal_blocks=True)
175+
else:
176+
raise ValueError("The selected manifold is not implemented.")
177+
178+
return result_dict
179+
180+
181+
def compute_operator(self, X, Y):
182+
"""
183+
Compute and store the low-rank operator and the full operator A.
184+
185+
:param X: matrix containing the snapshots x0,..x{n-1} by column.
186+
:type X: numpy.ndarray
187+
:param Y: matrix containing the snapshots x1,..x{n} by column.
188+
:type Y: numpy.ndarray
189+
:return: the (truncated) left-singular vectors matrix, the (truncated)
190+
singular values array, and the (truncated) right-singular vectors
191+
matrix of X.
192+
:rtype: numpy.ndarray, numpy.ndarray, numpy.ndarray
193+
"""
194+
U, s, V = compute_svd(X, self._svd_rank)
195+
196+
# Compute the corresponding Procrustes problem.
197+
result_dict = self._compute_procrustes(X, Y)
198+
199+
# Case 1: atilde was computed.
200+
if "atilde" in result_dict.keys():
201+
self._Atilde = result_dict["atilde"]
202+
self._eigenvalues, self._eigenvectors = np.linalg.eig(self._Atilde)
203+
self._modes = U.dot(self._eigenvectors)
204+
if self._compute_A:
205+
self._A = np.linalg.multi_dot([self._modes,
206+
np.diag(self._eigenvalues),
207+
np.linalg.pinv(self._modes)])
208+
else: # Cases 2 and 3.
209+
# Case 2: A was computed.
210+
if "A" in result_dict.keys():
211+
self._A = result_dict["A"]
212+
self._eigenvalues, self._modes = np.linalg.eig(self._A)
213+
else:
214+
# Case 3: eigenvalues and modes were computed.
215+
self._eigenvalues = result_dict["eigenvalues"]
216+
self._modes = result_dict["modes"]
217+
self._eigenvectors = U.conj().T.dot(self._modes)
218+
self._Atilde = np.linalg.multi_dot(
219+
[self._eigenvectors,
220+
np.diag(self._eigenvalues),
221+
np.linalg.pinv(self._eigenvectors)]
222+
)
223+
224+
return U, s, V
225+
226+
227+
class PiDMD(DMD):
228+
"""
229+
Physics-informed Dynamic Mode Decomposition.
230+
231+
:param manifold: the matrix manifold to restrict the full operator A to.
232+
The following matrix manifolds are permissible:
233+
- "unitary"
234+
- "uppertriangular"
235+
- "lowertriangular"
236+
- "diagonal"
237+
- "symmetric",
238+
- "skewsymmetric"
239+
- "toeplitz"
240+
- "hankel"
241+
- "circulant"
242+
- "circulant_unitary"
243+
- "circulant_symmetric"
244+
- "circulant_skewsymmetric"
245+
- "symmetric_tridiagonal"
246+
- "BC" (block circulant)
247+
- "BCTB" (BC with tridiagonal blocks)
248+
- "BCCB" (BC with circulant blocks)
249+
- "BCCBunitary" (BCCB and unitary)
250+
- "BCCBsymmetric" (BCCB and symmetric)
251+
- "BCCBskewsymmetric" (BCCB and skewsymmetric)
252+
:type manifold: str
253+
:param manifold_opt: option used to specify certain manifolds.
254+
If manifold is "diagonal", manifold_opt may be used to specify the
255+
width of the diagonal of A. If manifold_opt is an integer k, A is
256+
banded, with a lower and upper bandwidth of k-1. If manifold_opt is
257+
a tuple containing two integers k1 and k2, A is banded with a lower
258+
bandwidth of k1-1 and an upper bandwidth of k2-1. Finally, if
259+
manifold_opt is a numpy array of size (len(X), 2), the entries of
260+
manifold_opt are used to explicitly define the the upper and lower
261+
bounds of the indices of the non-zero elements of A.
262+
If manifold starts with "BC", manifold_opt must be a 2D tuple that
263+
specifies the desired dimensions of the blocks of A.
264+
Note that all other manifolds do not use manifold_opt.
265+
:type manifold_opt: int, tuple(int,int), or numpy.ndarray
266+
:param compute_A: Flag that determines whether or not to compute the full
267+
Koopman operator A.
268+
:type compute_A: bool
269+
:param svd_rank: the rank for the truncation; If 0, the method computes the
270+
optimal rank and uses it for truncation; if positive integer, the
271+
method uses the argument for the truncation; if float between 0 and 1,
272+
the rank is the number of the biggest singular values that are needed
273+
to reach the 'energy' specified by `svd_rank`; if -1, the method does
274+
not compute truncation. Default is -1, meaning no truncation.
275+
:type svd_rank: int or float
276+
:param tlsq_rank: rank truncation computing Total Least Square. Default is
277+
0, meaning no truncation.
278+
:type tlsq_rank: int
279+
:param opt: If True, amplitudes are computed like in optimized DMD (see
280+
:func:`~dmdbase.DMDBase._compute_amplitudes` for reference). If
281+
False, amplitudes are computed following the standard algorithm. If
282+
`opt` is an integer, it is used as the (temporal) index of the snapshot
283+
used to compute DMD modes amplitudes (following the standard
284+
algorithm).
285+
The reconstruction will generally be better in time instants near the
286+
chosen snapshot; however increasing `opt` may lead to wrong results
287+
when the system presents small eigenvalues. For this reason a manual
288+
selection of the number of eigenvalues considered for the analyisis may
289+
be needed (check `svd_rank`). Also setting `svd_rank` to a value
290+
between 0 and 1 may give better results. Default is False.
291+
:type opt: bool or int
292+
"""
293+
def __init__(
294+
self,
295+
manifold,
296+
manifold_opt=None,
297+
compute_A=False,
298+
svd_rank=-1,
299+
tlsq_rank=0,
300+
opt=False,
301+
):
302+
super().__init__(
303+
svd_rank=svd_rank,
304+
tlsq_rank=tlsq_rank,
305+
opt=opt,
306+
)
307+
self._Atilde = PiDMDOperator(
308+
manifold=manifold,
309+
manifold_opt=manifold_opt,
310+
compute_A=compute_A,
311+
svd_rank=svd_rank,
312+
)
313+
314+
@property
315+
def A(self):
316+
"""
317+
Get the full Koopman operator A.
318+
319+
:return: the full Koopman operator A.
320+
:rtype: numpy.ndarray
321+
"""
322+
return self.operator.A

0 commit comments

Comments
 (0)