notebook.community

Edit and run 





In [1]:



    


%pylab inline



    


















    






Populating the interactive namespace from numpy and matplotlib

















In [2]:



    


rc('figure',figsize=(10,10))

Собственно, создаем данные 2D:





In [3]:



    


N=1000
sx,sy = normal(2, size=(N,)), normal(5, size=(N,))
S = array([sx, sy])             # unmixed data
A = array([[-2, 3.5], [1, 5]])  # mixing matrix
M = A@S                   # mixed data
M = M - M.mean(axis=1).reshape(-1,1)
X = M.T

print( X.shape)



    


















    






(1000, 2)

















In [ ]:



    






    











In [4]:



    


scatter(M[0],M[1],marker='.', alpha=0.5, color='gray')
axis('square')



    


















    
Out[4]:








(-14.406609606973387, 22.6952976771035, -20.441586618886937, 16.66032066518995)

Eigendecomposition of covariance matrix (PCA)





In [5]:



    


cov(X.T)



    


















    
Out[5]:








array([[16.78764606, 15.34667695],
       [15.34667695, 25.51193037]])

















In [6]:



    


var(X,axis=0)



    


















    
Out[6]:








array([16.77085842, 25.48641844])

















In [7]:



    


X.T.dot(X)/(N-1)



    


















    
Out[7]:








array([[16.78764606, 15.34667695],
       [15.34667695, 25.51193037]])

















In [8]:



    


eig?
  • eigh returns eigen values in ascending order, and we need descending
  • eigh returns normalized eigenvectors as columns 




In [9]:



    


eigen_values, eigen_vectors = eig(cov(X.T))
eigen_values = eigen_values[::-1]           
eigen_vectors = eigen_vectors[:,::-1]
print('eigen values:', eigen_values)
print('\n')
print('eigen vectors:')
eigen_vectors



    


















    






eigen values: [37.10437306  5.19520337]


eigen vectors:










    
Out[9]:








array([[-0.6027396 , -0.79793795],
       [-0.79793795,  0.6027396 ]])

Eigenvalues are variances of projections on each eigenvector. Square roots of eigenvalues are standard deviations.





In [10]:



    


projections = X@eigen_vectors
print (projections.std(0))
print( eigen_values**0.5)



    


















    






[6.08828947 2.27815894]
[6.0913359  2.27929888]

PCA via SVD

\begin{equation} X = USV^T \end{equation}

svd returns eigenvectors of cov(X.T) as rows of vh. At least, it promises to:

The rows of v are the eigenvectors of a.H a.





In [11]:



    


u,s,vh = svd(X, full_matrices=0)

print('u shape: ', u.shape)
print('vh shape: ', vh.shape)



    


















    






u shape:  (1000, 2)
vh shape:  (2, 2)

















In [12]:



    


s



    


















    
Out[12]:








array([192.52861785,  72.0417113 ])

















In [13]:



    


projections2 = X.dot(vh.T)

print (s/(N-1)**0.5)
print (eigen_values**0.5)
print ("")
print (projections.std(0))
print (projections2.std(0))



    


















    






[6.0913359  2.27929888]
[6.0913359  2.27929888]

[6.08828947 2.27815894]
[6.08828947 2.27815894]

Note, s is equal to standard deviations along projections if divided by $\sqrt{(N-1)}$.

Matrix U can be regarded as the "mixing" matrix or coefficient loadings. If $X=USV^T$, then $XV = US$





In [14]:



    


eigen_vectors



    


















    
Out[14]:








array([[-0.6027396 , -0.79793795],
       [-0.79793795,  0.6027396 ]])

















In [15]:



    


vh



    


















    
Out[15]:








array([[-0.6027396 , -0.79793795],
       [ 0.79793795, -0.6027396 ]])

Vizualization

  1. Plot the cloud of points
  2. Plot principal axes as lines
  3. Plot principal vectors scaled by standard deviation of projections onto each axis




In [16]:



    


scatter(M[0],M[1],marker='.', alpha=0.5, color='gray')
sigmas = s/sqrt(N-1)


xfit = linspace(-10,10,100)

k1 = vh[0,1]/vh[0,0]
k2 = vh[1,1]/vh[1,0]

arrow(0,0, -sigmas[0]*vh[0,0], -sigmas[0]*vh[0,1], color='red', width=0.2,head_width=1.,alpha=0.75)
arrow(0,0, sigmas[1]*vh[1,0], sigmas[1]*vh[1,1], color='blue', width=0.2,head_width=1., alpha=0.75)

plot(xfit, xfit*k1, 'r-',lw=0.25)
plot(xfit, xfit*k2, 'b-', lw=0.25)
axis('equal')



    


















    
Out[16]:








(-14.406609606973387,
 15.733690647952889,
 -20.441586618886937,
 16.66032066518995)

Данные в виде серии изображений





In [17]:



    


ls media



    


















    






aquarium.h5  hexbug1.mp4  reena-data-best.h5  two-windmills.h5  V_ch7_1.tif

















In [18]:



    


from imfun import fseq, ui



    


















    






Can't load imreg package, affine and homography registrations won't work

















In [19]:



    


frames = fseq.from_tiff("media/V_ch7_1.tif")
data = frames.as3darray()
print (data.shape)
ui.group_maps(data[:300:5],10, imkw=dict(cmap='gray'), colorbar=False, figscale=4)
tight_layout();



    


















    






(670, 240, 320)

На рис. показана развертка первых 300 кадров (показан каждый 5-й). Данные представляют собой эволюцию спиральной волны. Всего в записи 670 кадров.





In [20]:



    


xc = data[:,20:180,110:280]
xcn = array([a - a.mean() for a in (xc - xc.mean(axis=0))])
sh = xc.shape[1:]
print( sh)



    


















    






(160, 170)

















In [21]:



    


xcn.std()



    


















    
Out[21]:








23.13591004108793

















In [22]:



    


%matplotlib qt



    











In [23]:



    


X = xcn + randn(*xcn.shape)*50



    











In [24]:



    


ui.Picker(fseq.from_array(X)).start()



    


















    
Out[24]:








(<matplotlib.axes._subplots.AxesSubplot at 0x7ff2d0dae470>,
 <matplotlib.image.AxesImage at 0x7ff30c9e9390>,
 <imfun.ui.Picker_.Picker at 0x7ff30eeaf160>)










    
























In [25]:



    


%matplotlib inline



    











In [26]:



    


ui.group_maps(X[:50:5],5, imkw=dict(cmap='gray'), colorbar=False,figscale=3.5)
tight_layout();



    


















    
























In [27]:



    


sh[0]*sh[1]



    


















    
Out[27]:








27200

















In [28]:



    


U, S, Vh = svd(X.reshape(670,prod(sh)), full_matrices=False);
print('U shape:', U.shape)
print('Vh shape:', Vh.shape)



    


















    






U shape: (670, 670)
Vh shape: (670, 27200)

















In [29]:



    


ui.group_maps(Vh[:25].reshape(-1,*sh),5,colorbar=False,figscale=5)
tight_layout()



    


















    
























In [30]:



    


plot(U[:,0])
plot(U[:,1])
plot(U[:,3])



    


















    
Out[30]:








[<matplotlib.lines.Line2D at 0x7ff30ecd39b0>]










    
























In [31]:



    


fh, axs = subplots(4,5,sharex=True,sharey=True, figsize=(14,9))
for j,ax in enumerate(axs.ravel()):
    ax.plot(U[:,j], 'k-')
    ax.set_title(j+1)
setp(axs[-1,-1], xticks = [0,300,600])



    


















    
Out[31]:








[<matplotlib.axis.XTick at 0x7ff30ee8f4a8>,
 <matplotlib.axis.XTick at 0x7ff30ee6ad68>,
 <matplotlib.axis.XTick at 0x7ff30ec985c0>]










    
























In [32]:



    


plot(S, 's-')



    


















    
Out[32]:








[<matplotlib.lines.Line2D at 0x7ff30e559208>]










    
























In [33]:



    


figure(figsize=(12,6))
subplot(121)
plot(range(1,31),(S[:30]**2)/sum(S**2), 'ko')
title(u"Доля объясненной дисперсии")
xlabel(u"Номер компоненты")
ylabel("")

subplot(122)
title(u"Кумулятивная доля объясненной дисперсии")
plot(range(1,31),cumsum(S[:30]**2)/sum(S**2), 'ks')
xlabel(u"Номер компоненты")



    


















    
Out[33]:








Text(0.5, 0, 'Номер компоненты')

Optimal singular value threshold

See [Gavish, Donoho 2014] and http://www.pyrunner.com/weblog/2016/08/01/optimal-svht/





In [34]:



    


def omega_approx(beta):
    return 0.56*beta**3 - 0.95*beta**2 + 1.82*beta + 1.43

def svht(sv, sh):
    m,n = sh
    if m>n: 
        m,n=n,m
    omg = omega_approx(m/n)
    return omg*median(sv)

def estimate_rank_gd(sv,sh):
    th = svht(sv,sh)
    return np.sum(sv>th)



    











In [35]:



    


r= estimate_rank_gd(S,Vh.shape)
r



    


















    
Out[35]:








15

Т.е. первых 15 компонент достаточно для описания динамики





In [36]:



    


close()
figure(figsize=(14,6))
plot(range(1,len(S)+1), cumsum(S**2)/sum(S**2), 'k-', lw=2)
title(u"Кумулятивная доля объясненной дисперии")
axvline(25, color='r')
axvline(r, color='m',ls='--')

text(27,0.9, "n=25", color='r')
text(r+3,0.95, "n=%d"%r, color='m')

xlabel(u"Номер компоненты")



    


















    
Out[36]:








Text(0.5, 0, 'Номер компоненты')

Приближенная реконструкция:





In [37]:



    


Xapprox = (U[:,:r]@diag(S[:r])@Vh[:r]).reshape(xcn.shape)



    











In [38]:



    


ui.group_maps(xcn[:25:5],5, imkw=dict(cmap='gray'), colorbar=False,figscale=3.5)
tight_layout();



    


















    
























In [39]:



    


ui.group_maps(Xapprox[:25:5],5, imkw=dict(cmap='gray'), colorbar=False,figscale=3.5)
tight_layout();



    


















    
























In [40]:



    


%matplotlib qt



    











In [41]:



    


p = ui.Picker(fseq.from_array(X));
p.start()
p.clims = [(-50,100)]



    


















    
























In [42]:



    


p = ui.Picker(fseq.from_array(Xapprox));
p.start()
p.clims = [(-50,100)]

Реконструкция с фильтрацией временных мод





In [43]:



    


%matplotlib inline



    











In [44]:



    


rc('figure', figsize=(14,10))



    











In [45]:



    


from imfun.filt import l2spline



    











In [46]:



    


plot(U[:,10],'gray')
plot(l2spline(U[:,10],5))



    


















    
Out[46]:








[<matplotlib.lines.Line2D at 0x7ff2f0392cc0>]










    
























In [47]:



    


plot(U[:,20],'gray')
plot(l2spline(U[:,20],5))



    


















    
Out[47]:








[<matplotlib.lines.Line2D at 0x7ff2f03005c0>]










    
























In [48]:



    


Usmoothed = array([l2spline(uv,5) for uv in U.T[:r]]).T



    











In [49]:



    


Xapprox2 = (Usmoothed@diag(S[:r])@Vh[:r]).reshape(xcn.shape)



    











In [50]:



    


%matplotlib qt



    











In [51]:



    


p = ui.Picker(fseq.FStackColl([fseq.from_array(X),
                               fseq.from_array(Xapprox),
                               fseq.from_array(Xapprox2),
                               fseq.from_array(Xapprox-Xapprox2)]));
p.start()
p.clims = [(-50,100)]*4

exact DMD (Fortran/Matlab data order)

if $X$ is a matrix of inputs or states at times $t=n\Delta t$, and $Y$ is a matrix of outputs or states at times $t+\Delta t$, with each data set (e.g. video frame) raveled as a column vector, then

Low-rank $r$ truncated SVD approximation \begin{equation} X = U\Sigma V^*\,, \end{equation}

\begin{equation} Y = AX \end{equation}\begin{equation} A = YX^\dagger = YV\Sigma^{-1}U^* \end{equation}\begin{equation} \tilde{A} = U^*AU = U^*YV\Sigma^{-1} \end{equation}\begin{equation} \tilde{A}W = W\lambda \end{equation}\begin{equation} \Phi = YV\Sigma^{-1}W \end{equation} 




In [344]:



    


def dmdf(X,Y=None, r=None,sort_explained=True,smooth_coefs=False,smooth_tsigma=1.5):
    if Y is None:
        Y = X[:,1:]
        X = X[:,:-1]
    U,sv,Vh = svd(X,False)
    if r is None:
        tr = svht(sv,X.shape)
        r = np.sum(sv>tr)
        print('Setting rank to ', r)
        #if not r%2:
        #    r += 1
    if smooth_coefs:
        Vh = array([l2spline(v,smooth_sigma) for v in Vh[:r]])
    sv = sv[:r]
    V = Vh[:r].conj().T
    Uh = U[:,:r].conj().T
    B = Y@V@(diag(1/sv))
    Atilde = Uh@B
    lam, W = eig(Atilde)
    Phi = B@W
    #print(Vh.shape)
    # approx to b
    def _bPOD(i):
        alpha1 =diag(sv[:r])@Vh[:r,i]
        return lstsq(Atilde@W,alpha1,rcond=None)[0]
    bPOD = _bPOD(0)
    stats = (None,None)
    if sort_explained:
        proj_dmd = array([_bPOD(i) for i in range(Vh.shape[1])])
        dmd_std = proj_dmd.std(0)
        dmd_mean = abs(proj_dmd).mean(0)
        stats = (dmd_mean,dmd_std)
        kind = argsort(dmd_std)[::-1]
    else:
        kind = arange(r)[::-1] # from slow to fast
    Phi = Phi[:,kind]
    lam = lam[kind]
    bPOD=bPOD[kind]
    
    return lam, Phi,bPOD,stats

def plot_eigenvals(lam,cont_time=True,ax=None):
    if ax is None:
        f,ax = subplots(1,1)
    if cont_time:
        omega = log(lam)/2/pi
    else:
        omega = lam
    ax.plot(omega.real, omega.imag, 'o')
    omeg_str = r'\omega' if cont_time else r'\lambda'
    ax.set_xlabel('$\\Re(%s)$'%omeg_str)
    ax.set_ylabel('$\\Im(%s)$'%omeg_str)
    ax.grid(True)
    if not cont_time:
        c = Circle((0,0),radius=1.,ls='--',edgecolor='k',facecolor='none')
        ax.add_patch(c)
        ax.axis('square')
        ax.set_title('Discrete-time eigenvalues')
    else:
        ax.set_title('Continuous-time eigenvalues')
    #gcf()



    











In [53]:



    


Xdmd = X.reshape(-1,prod(sh)).T



    











In [54]:



    


lam, Phi,bx,_ = dmdf(Xdmd, r=r*2,smooth_coefs=True)



    











In [55]:



    


Phi.shape



    


















    
Out[55]:








(27200, 30)

















In [56]:



    


Phi_img = Phi.T.reshape(-1,*sh)



    











In [57]:



    


%matplotlib inline



    











In [58]:



    


ui.group_maps(Phi_img.real,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [59]:



    


ui.group_maps(Phi_img.imag,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [60]:



    


figure(figsize=(10,10))
plot_eigenvals(lam,False,gca())
xlim(0.75,1.1); ylim(-0.25,0.25)



    


















    
Out[60]:








(-0.25, 0.25)










    






/opt/anaconda/lib/python3.7/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['cmtt10'] not found. Falling back to DejaVu Sans.
  (prop.get_family(), self.defaultFamily[fontext]))
/opt/anaconda/lib/python3.7/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['cmb10'] not found. Falling back to DejaVu Sans.
  (prop.get_family(), self.defaultFamily[fontext]))
/opt/anaconda/lib/python3.7/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['cmss10'] not found. Falling back to DejaVu Sans.
  (prop.get_family(), self.defaultFamily[fontext]))
/opt/anaconda/lib/python3.7/site-packages/matplotlib/font_manager.py:1241: UserWarning: findfont: Font family ['DejaVu Sans Display'] not found. Falling back to DejaVu Sans.
  (prop.get_family(), self.defaultFamily[fontext]))

Reconstcrution of dynamics

\begin{equation} X_0 = \Phi\mathbf{z} \end{equation}\begin{equation} X_k = \Phi\Lambda^k\mathbf{z}\,,\;\text{for}\; k=0,\dots \end{equation} 




In [61]:



    


%time zfirst = lstsq(Phi, Xdmd[:,0],rcond=None)[0]
%time zlast = lstsq(Phi, Xdmd[:,-1],rcond=None)[0]



    


















    






CPU times: user 61.5 ms, sys: 0 ns, total: 61.5 ms
Wall time: 17.9 ms
CPU times: user 48 ms, sys: 0 ns, total: 48 ms
Wall time: 11.6 ms

















In [62]:



    


zlast.shape



    


















    
Out[62]:








(30,)

















In [63]:



    


Xrec = Phi@diag(zfirst)@vander(lam,len(X),True)
Xrec = Xrec.T.real.reshape(X.shape)



    











In [64]:



    


Xrec2 = Phi@diag(zlast)@vander(lam,300,True)
Xrec2 = Xrec2.T.real.reshape(300,*sh)



    











In [65]:



    


%matplotlib qt



    











In [66]:



    


p = ui.Picker(fseq.FStackColl([fseq.from_array(X),
                               fseq.from_array(Xapprox),
                               fseq.from_array(Xrec)]));
p.start()
p.clims = [(-20,100)]*10



    


















    
























In [67]:



    


p=ui.Picker(fseq.from_array(Xrec2));
p.start()
p.clims = [(-20,100)]*10



    


















    
























In [68]:



    


def one_step_predictions(X,lam,Phi):
    out = []
    for i in range(X.shape[1]-1):
        Y = X[:,i+1]
        b = lstsq(Phi, X[:,i],rcond=None)[0]
        Ypred = Phi@diag(lam)@b.reshape(-1,1)
        out.append(norm(Y-Ypred.ravel())/len(Y))
    return array(out)



    











In [69]:



    


v = one_step_predictions(Xdmd,lam,Phi)



    











In [70]:



    


plot(v)



    


















    
Out[70]:








[<matplotlib.lines.Line2D at 0x7ff2f02b26d8>]










    
























In [71]:



    


norm?

Windmills





In [72]:



    


frames = fseq.from_hdf5('media/two-windmills.h5')



    











In [73]:



    


frames.frame_shape



    


















    
Out[73]:








(90, 160)

















In [74]:



    


ui.group_maps(frames[:100:10],5,colorbar=False,imkw=dict(cmap='gray'))
gcf()



    


















    
Out[74]:

















    
























In [75]:



    


frame_crop = (slice(10,None),slice(0,80))



    











In [76]:



    


cropped_frames = array([f[frame_crop] for f in frames])



    











In [77]:



    


ui.group_maps(cropped_frames[:100:10],5,colorbar=False,imkw=dict(cmap='gray'))
gcf()



    


















    
Out[77]:

















    
























In [78]:



    


sh = cropped_frames[0].shape
sh



    


















    
Out[78]:








(80, 80)

















In [79]:



    


Xdmd = cropped_frames.reshape(-1,prod(sh)).T



    











In [80]:



    


lam, Phi,bx,_ = dmdf(Xdmd, r=r*2,smooth_coefs=False)



    











In [81]:



    


%matplotlib inline



    











In [82]:



    


Phi_img = Phi.T.reshape(-1,*sh)



    











In [83]:



    


ui.group_maps(Phi_img.real,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [84]:



    


plot_eigenvals(lam,False)



    


















    
























In [85]:



    


%matplotlib qt



    











In [86]:



    


ui.Picker(fseq.from_array(cropped_frames)).start()



    


















    
Out[86]:








(<matplotlib.axes._subplots.AxesSubplot at 0x7ff3323e5780>,
 <matplotlib.image.AxesImage at 0x7ff30bcc6080>,
 <imfun.ui.Picker_.Picker at 0x7ff3323e55f8>)










    
























In [87]:



    


%time zfirst = lstsq(Phi, Xdmd[:,0],rcond=None)[0]



    


















    






CPU times: user 16.5 ms, sys: 3.55 ms, total: 20 ms
Wall time: 5.25 ms

















In [88]:



    


Xrec = Phi@diag(zfirst)@vander(lam,len(cropped_frames),True)
Xrec = Xrec.T.real.reshape(cropped_frames.shape)



    











In [89]:



    


%matplotlib tk



    


















    






Warning: Cannot change to a different GUI toolkit: tk. Using qt instead.

















In [90]:



    


p = ui.Picker(fseq.from_array(Xrec))
p.start()



    


















    
Out[90]:








(<matplotlib.axes._subplots.AxesSubplot at 0x7ff30bcfe4e0>,
 <matplotlib.image.AxesImage at 0x7ff30bb7f0b8>,
 <imfun.ui.Picker_.Picker at 0x7ff30eb39d68>)

Full frame





In [91]:



    


sh = frames.frame_shape



    











In [92]:



    


Xdmd = reshape(frames.data, (-1,prod(sh))).T



    











In [93]:



    


lam, Phi,bx,_ = dmdf(Xdmd, r=r*2,smooth_coefs=False)



    











In [94]:



    


%matplotlib inline



    











In [95]:



    


Phi_img = Phi.T.reshape(-1,*sh)



    











In [96]:



    


ui.group_maps(Phi_img.real,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [97]:



    


plot_eigenvals(lam,False)



    


















    
























In [98]:



    


%time zfirst = lstsq(Phi, Xdmd[:,0],rcond=None)[0]



    


















    






CPU times: user 34.6 ms, sys: 0 ns, total: 34.6 ms
Wall time: 8.78 ms

















In [99]:



    


Xrec = Phi@diag(zfirst)@vander(lam,len(frames),True)
Xrec = Xrec.T.real.reshape(frames.data.shape)



    











In [100]:



    


%matplotlib qt



    











In [101]:



    


p = ui.Picker(fseq.FStackColl([frames,fseq.from_array(Xrec)]))
p.start()



    


















    
Out[101]:








(<matplotlib.axes._subplots.AxesSubplot at 0x7ff30b549978>,
 <matplotlib.image.AxesImage at 0x7ff30b5875f8>,
 <imfun.ui.Picker_.Picker at 0x7ff30bcc6f60>)

















In [ ]:

Aquarium





In [285]:



    


%matplotlib inline



    











In [286]:



    


frameseq = fseq.FStackColl_hdf5('/home/data/consumer/videos/aquarium.h5')



    


















    






The file /home/data/consumer/videos/aquarium.h5 has the following data sets: 0, 1, 2

















In [287]:



    


frames = frameseq[:].astype(float32)



    











In [288]:



    


frames = frames[...,0] + frames[...,1]



    











In [289]:



    


frames = frames[:,::2,::2]



    











In [290]:



    


#frames  = frames - frames.mean(0)



    











In [291]:



    


ui.Picker(fseq.from_array(frames)).start()



    


















    
Out[291]:








(<matplotlib.axes._subplots.AxesSubplot at 0x7ff2a0703a58>,
 <matplotlib.image.AxesImage at 0x7ff2c7d5dd68>,
 <imfun.ui.Picker_.Picker at 0x7ff2a0703dd8>)










    
























In [292]:



    


frames.shape



    


















    
Out[292]:








(657, 270, 480)

















In [293]:



    


sh = frames[0].shape
sh



    


















    
Out[293]:








(270, 480)

















In [294]:



    


Xdmd = reshape(frames, (-1,prod(sh))).T



    











In [295]:



    


Xdmd.shape



    


















    
Out[295]:








(129600, 657)

















In [296]:



    


u,s,vh = svd(Xdmd,full_matrices=False)



    











In [297]:



    


ui.group_plots(vh[:15],ncols=1,figsize=(20,20))



    


















    
























In [345]:



    


ui.group_plots([l2spline(v_,1.5) for v_ in vh[:15]],ncols=1,figsize=(20,20))



    











In [299]:



    


u.shape



    


















    
Out[299]:








(129600, 657)

















In [300]:



    


ui.group_maps([f.reshape(sh) for f in u.T[:15]],colorbar=False)
tight_layout()



    


















    
























In [222]:



    


estimate_rank_gd(s,vh.shape)



    


















    
Out[222]:








44

















In [223]:



    


plot(s[1:],'s-')



    


















    
Out[223]:








[<matplotlib.lines.Line2D at 0x7ff2efea5390>]










    
























In [316]:



    


%time lam, Phi,bx,_ = dmdf(Xdmd, r=60, smooth_coefs=True, smooth_sigma=1.5)



    


















    






CPU times: user 20.7 s, sys: 2.12 s, total: 22.8 s
Wall time: 6.42 s

















In [368]:



    


%matplotlib inline



    











In [361]:



    


Phi.shape



    


















    
Out[361]:








(129600, 60)

















In [318]:



    


Phi_img = Phi.T.reshape(-1,*sh)



    











In [319]:



    


ui.group_maps(Phi_img.real,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [348]:



    


ui.group_maps(Phi_img.imag,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [349]:



    


plot_eigenvals(lam,False)



    


















    
























In [350]:



    


from scipy import ndimage as ndi



    











In [369]:



    


Phi_img_x = array([l2spline(phi,3) for phi in Phi_img])



    











In [370]:



    


ui.group_maps(Phi_img_x.real,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [353]:



    


ui.group_maps(Phi_img_x.imag,colorbar=False, figscale=5,imkw=dict(cmap='bwr'))
tight_layout()



    


















    
























In [363]:



    


Phi_x = array([ravel(phi) for phi in Phi_img_x]).T



    











In [321]:



    


lam.imag



    


















    
Out[321]:








array([-0.06920962,  0.06920962, -0.06194206,  0.06194206,  0.07388886,
       -0.07388886,  0.02544111, -0.02544111, -0.01736248,  0.01736248,
       -0.0887362 ,  0.0887362 ,  0.00937419, -0.00937419,  0.01097757,
       -0.01097757,  0.        ,  0.02371824, -0.02371824, -0.0047406 ,
        0.0047406 , -0.0483461 ,  0.0483461 ,  0.06899043, -0.06899043,
       -0.01284443,  0.01284443,  0.        ,  0.01950454, -0.01950454,
       -0.02001715,  0.02001715, -0.06651811,  0.06651811, -0.02166664,
        0.02166664,  0.        ,  0.        ,  0.04459919, -0.04459919,
       -0.02835171,  0.02835171, -0.03856583,  0.03856583,  0.        ,
        0.        ,  0.        , -0.01286546,  0.01286546,  0.        ,
       -0.06999981,  0.06999981,  0.        , -0.1422803 ,  0.1422803 ,
       -0.02233337,  0.02233337,  0.        ,  0.        ,  0.        ])

















In [328]:



    


def project_to_unitary_circle(x,th=0.05):
    if abs(x.imag) > th:
        r = np.sqrt(x.real**2 + x.imag**2)
        phi = arctan2(x.imag, x.real)
        return cos(phi) + 1j*sin(phi)
    else:
        return x



    











In [329]:



    


lamx = array([project_to_unitary_circle(l) for l in lam])



    











In [330]:



    


plot_eigenvals(lamx,False)



    


















    
























In [364]:



    


%time zfirst = lstsq(Phi, Xdmd[:,0],rcond=None)[0]

%time zfirstx = lstsq(Phi_x, Xdmd[:,0],rcond=None)[0]

%time zlast = lstsq(Phi, Xdmd[:,-1],rcond=None)[0]



    


















    






CPU times: user 1.04 s, sys: 104 ms, total: 1.15 s
Wall time: 288 ms
CPU times: user 873 ms, sys: 71.8 ms, total: 944 ms
Wall time: 234 ms
CPU times: user 1.02 s, sys: 76 ms, total: 1.09 s
Wall time: 273 ms

















In [332]:



    


Xrec = Phi@diag(zfirst)@vander(lam,len(frames),True)
Xrec = Xrec.T.real.reshape(frames.data.shape)



    











In [371]:



    


Xrecx = Phi_x@diag(zfirstx)@vander(lamx,len(frames),True)
Xrecx = Xrecx.T.real.reshape(frames.data.shape)



    











In [373]:



    


%matplotlib qt



    











In [342]:



    


imshow(frames[0])



    


















    
Out[342]:








<matplotlib.image.AxesImage at 0x7ff2a0437ba8>

















In [374]:



    


p = ui.Picker(fseq.FStackColl([fseq.from_array(frames),
                               fseq.from_array(Xrec),
                               fseq.from_array(Xrecx)
                              ]))
p.start()
p.clims = [(-450, 550)]*4



    


















    






180.80
area <itertools._grouper object at 0x7ff2a0996c18>
signal shape: (3, 657)
area <itertools._grouper object at 0x7ff2a494e518>
signal shape: (3, 657)
area <itertools._grouper object at 0x7ff2a4815358>
signal shape: (3, 657)

















In [ ]:



    






    











In [141]:



    


Xrec2 = Phi@diag(zlast)@vander(lam,300,True)
Xrec2 = Xrec2.T.real.reshape(300,*sh)



    











In [142]:



    


p2 = ui.Picker(fseq.FStackColl([fseq.from_array(frames),fseq.from_array(Xrec2)]))
p2.start()



    


















    
Out[142]:








(<matplotlib.axes._subplots.AxesSubplot at 0x7ff3097b9eb8>,
 <matplotlib.image.AxesImage at 0x7ff30ad49710>,
 <imfun.ui.Picker_.Picker at 0x7ff3096e1748>)

















In [160]:



    


close('all')



    











In [203]:



    


#import pydmd



    











In [ ]:

Content source: abrazhe/sipc

Similar notebooks:

notebook.community | gallery | about