-
Notifications
You must be signed in to change notification settings - Fork 6
/
dmd.py
105 lines (86 loc) · 2.97 KB
/
dmd.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
import numpy as np
from numpy import pi
from matplotlib import pyplot as plt
def make_VandermondeT(omega, time):
"""
Compute the transpose of the Vandermonde
matrix from the times and frequencies
Parameters
----------
omega: 1D numpy array of complex values
(r = truncation rank of the DMD)
DMD frequencies
time: 1D numpy array of floats
(M = number of time samples)
Time range for the DMD
Returns
----------
VandermondeT: 2D numpy array of complex values
(r, M)
Transpose of the Vandermonde matrix for DMD reconstructions
"""
VandermondeT = np.exp(np.outer(time, omega))
return VandermondeT
def dmd(data, r, time, M_train):
"""
Compute the DMD of a set of time series
Parameters
----------
data: 2D numpy array of spacetime data
(N = number of spatial locations, M = number of time samples)
Data to use the DMD on
r: integer
Truncation rank for the DMD
time: 1D numpy array of floats
(M = number of time samples)
Time range for the DMD
M_train: integer
The length of the time range for building the DMD model,
the remainder is used for testing the model.
Returns
----------
Bfield: 2D numpy array of floats
(N = number of spatial locations, M = number of time samples)
DMD reconstruction of the data variable
"""
tsize = len(time)
time = time / 1e6
Qsize = int(np.shape(data)[0] / 6)
Bfield = np.zeros((np.shape(data)[0], tsize - M_train), dtype='complex')
dt = 1e-6
X = data[:, 0:M_train]
Xprime = data[:, 1:M_train + 1]
Udmd, Sdmd, Vdmd = np.linalg.svd(X, full_matrices=False)
Vdmd = np.transpose(Vdmd)
Udmd = Udmd[:, 0:r]
Sdmd = Sdmd[0:r]
Vdmd = Vdmd[:, 0:r]
S = np.diag(Sdmd)
A = np.dot(np.dot(np.transpose(Udmd), Xprime), Vdmd / Sdmd)
eigvals, Y = np.linalg.eig(A)
Bt = np.dot(np.dot(Xprime, Vdmd / Sdmd), Y)
omega = np.log(eigvals) / dt
VandermondeT = make_VandermondeT(omega, time - time[0])
Vandermonde = np.transpose(VandermondeT)
q = np.conj(np.diag(np.dot(np.dot(np.dot(
Vandermonde[:, :M_train], Vdmd), np.conj(S)), Y)))
P = np.dot(np.conj(np.transpose(Y)), Y) * np.conj(
np.dot(Vandermonde[:, :M_train], np.conj(
VandermondeT[:M_train, :])))
b = np.dot(np.linalg.inv(P), q)
c = 0.5 * (b + np.conj(b)).real
c = 500 * c / np.max(c)
energy_sort = np.flip(np.argsort(c))
plt.figure()
omega_sizes = c
for i in range(len(c)):
omega_sizes[i] = max(omega_sizes[i], 50)
# plot the frequencies
plt.scatter(omega.imag/(2 * pi * 1e3),
omega.real/(2 * pi * 1e3), s=omega_sizes, c='k')
plt.savefig('Pictures/dmd_freqs.pdf')
for mode in range(r):
Bfield += 0.5 * b[mode] * np.outer(Bt[:, mode],
Vandermonde[mode, M_train:])
Bfield += np.conj(Bfield)
return Bfield.real