Page MenuHomec4science
Files F61979581

dmd.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 10, 04:33

dmd.py
View Options

def svht(X, sv=None):
from numpy import squeeze, median
from scipy.linalg import svdvals
# svht for sigma unknown
m, n = sorted(X.shape) # ensures m <= n
beta = m / n # ratio between 0 and 1
if sv is None:
sv = svdvals(X)
sv = squeeze(sv)
omega_approx = 0.56 * beta**3 - 0.95 * beta**2 + 1.82 * beta + 1.43
return median(sv) * omega_approx
def dmd(X, Y, truncate=None, do_svht=False):
from numpy import dot, diag, power, array, zeros
from scipy.linalg import svd, svdvals
from numpy.linalg import inv, eig
if truncate == 0:
# return empty vectors
mu = np.array([], dtype='complex')
Phi = np.zeros([X.shape[0], 0], dtype='complex')
else:
# 2.8-2.11 in Kutz et al. 2016
U2, Sig2, Vh2 = svd(X, False) # SVD of input matrix
# determine rank-reduction
if do_svht:
_sv = svdvals(_D)
tau = svht(_D, sv=_sv)
r = sum(_sv > tau)
else:
r = len(Sig2)
U = U2[:, :r]
Sig = diag(Sig2)[:r, :r]
V = Vh2.conj().T[:, :r]
iSig = inv(Sig)
# S tilde in Scmid 2010, A tilde in Kou and Zhang (eq.16)
Atil = dot(dot(dot(U.conj().T, Y), V), iSig)
mu, W = eig(Atil)
# See 2.11 Kutz et al. 2016
Phi = dot(dot(dot(Y, V), iSig), W)
return mu, Phi
# multi resolution Dynamic Mode Decomposition
def mrdmd(D, max_levels=3, max_cycles=1, do_svht=True,
window="boxcar", overlap=0,
level=0, bin_num=0, offset=0, ol=0):
"""Compute the multi-resolution DMD on the dataset `D`, returning a list of
nodes in the hierarchy. Each node represents a particular "time bin"
(window in time) at a particular "level" of the recursion (time scale).
The node is an object consistingof the various data structures generated by
the DMD at its corresponding level and time bin.
max_levels: parameter which controls the maximum number of levels.
max_cycles: parameter controlling the subsampling of data at each level,
a larger number implies that data is subsampled using a
bigger step. It also controls which modes are considered
slow and which ones fast (larger number=higher frequency).
do_svht: parameter controlling whether or not to perform optimal singular
value hard thresholding.
window: selects the window used for splitting the dataset. By default,
a "boxcar" window is used with no overlap. Alternatively, any
window from scipy.signal module can be used.
overlap: Overlap between adiacent windows, given in fraction of the total
length, thus between 0 and 1 (default: 0).
NOTES:
max_levels, max_cycles, do_svht and ol are used for recusion bookeeping
and should not be modified
The window function is not used at the lowest level of decomposition."""
import numpy as np
import scipy.signal as ssig
from numpy import dot, multiply, diag, power, sum, abs
from numpy import pi, exp, sin, cos
from numpy.linalg import inv, eig, solve, norm
from scipy.linalg import svd, svdvals
from math import floor, ceil # python 3.x
if overlap >= 1 or overlap < 0:
raise ValueError("Overlap must be within 0"
"(included) and 1 (excluded)")
# 4 times nyquist limit to capture cycles
nyq = 8 * max_cycles
# time bin size
bin_size = D.shape[1]
if bin_size < nyq:
return []
# extract subsamples
step = int(floor(bin_size / nyq))
_D = D[:, ::step]
# window data
if window != "boxcar" and level > 0:
winfunc = getattr(ssig, window)
w = winfunc(_D.shape[1])
w /= w.sum()
_D *= w
X = _D[:, :-1]
Y = _D[:, 1:]
# determine rank-reduction
if do_svht:
_sv = svdvals(_D)
tau = svht(_D, sv=_sv)
r = sum(_sv > tau)
else:
r = min(X.shape)
# compute dmd
mu, Phi = dmd(X, Y, r)
# frequency cutoff (oscillations per sampling time)
rho = max_cycles / bin_size
# identify slow eigenvalues (as boolean mask)
slow = (np.abs(np.log(mu) / (2 * pi * step))) <= rho
n = sum(slow) # number of slow modes
# extract slow modes (perhaps empty)
mu = mu[slow]
Phi = Phi[:, slow]
if n > 0:
# find optimal b solution, 2.13 in Kutz et al. 2016
alpha = dot(np.linalg.pinv(Phi), D[:, 0])
# eq 2.12 in Kutz et al and Jovanovic et al 2014
Vand = np.vander(power(mu, 1/step), bin_size, True)
# time evolution
Psi = dot(diag(alpha), Vand)
else:
# zero time evolution
alpha = np.array([], dtype='complex')
Psi = np.zeros([0, bin_size], dtype='complex')
# dmd reconstruction
D_dmd = dot(Phi, Psi)
Ij = sum(abs(Psi), axis=1) * norm(Phi, axis=0)**2 * step
# remove influence of slow modes
D = D - D_dmd
# record keeping
node = type('Node', (object,), {})()
node.level = level # level of recursion
node.bin_num = bin_num # time bin number
node.bin_size = bin_size # time bin size
node.start = offset # starting index
node.stop = offset + bin_size # stopping index
node.overlap = ol # number of overlapping points
node.step = step # step size
node.rho = rho # frequency cutoff
node.r = r # rank-reduction
node.n = n # number of extracted modes
node.mu = np.log(mu) / step # extracted eigenvalues
node.Phi = Phi # extracted DMD modes
node.Psi = Psi # extracted time evolution
node.Ij = Ij # extracted optimal alpha vector
node.window = window
nodes = [node]
# split data into two and do recursion
if level < max_levels:
split = ceil(bin_size / 2) # where to split
ol = floor(split * overlap)
nodes += mrdmd(
D[:, :split + ol],
level=level+1,
bin_num=2*bin_num,
ol=ol,
offset=offset,
max_levels=max_levels,
max_cycles=max_cycles,
do_svht=do_svht,
overlap=overlap,
window=window
)
nodes += mrdmd(
D[:, split - ol:],
level=level+1,
bin_num=2*bin_num+1,
ol=-ol,
offset=offset+split-ol,
max_levels=max_levels,
max_cycles=max_cycles,
do_svht=do_svht,
overlap=overlap,
window=window
)
return nodes
def stitch(nodes, level):
import numpy as np
# get length of time dimension
start = min([nd.start +
max(0, -nd.overlap) # overlap left
for nd in nodes])
stop = max([nd.stop +
min(0, -nd.overlap) # overlap right
for nd in nodes])
t = stop - start
# extract relevant nodes
nodes = [nd for nd in nodes if nd.level == level]
nodes = sorted(nodes, key=lambda nd: nd.bin_num)
# stack DMD modes
Phi = np.hstack([nd.Phi for nd in nodes])
# stack eigenvalues
mu = np.hstack([nd.mu for nd in nodes])
# allocate zero matrix for time evolution
nmodes = sum([nd.n for nd in nodes])
Psi = np.zeros([nmodes, t], dtype='complex')
# copy over time evolution for each time bin
i = 0
for nd in nodes:
_nmodes = nd.Psi.shape[0]
# deal with overlap
Psi[i:i+_nmodes,
nd.start+max(0, -nd.overlap):nd.stop+min(0, -nd.overlap)] = \
nd.Psi[:, max(0, -nd.overlap):nd.Psi.shape[1]+min(0, -nd.overlap)]
i += _nmodes
return Phi, Psi, mu

Event Timeline

c4science · Help