SSA奇异谱分析py代码

适合于初学者的SSA代码。

# -*- coding: utf-8 -*-
"""
Created on Fri Oct 16 15:10:29 2020

@author: ZHOU
"""

'''
                                  Python Setup
'''

# Load the usual suspects:
import numpy as np
from numpy import pi
import matplotlib.pyplot as plt
import pandas as pd
import os

# Fiddle with figure settings here:
plt.rcParams['figure.figsize'] = (10,8)
plt.rcParams['font.size'] = 14
plt.rcParams['image.cmap'] = 'plasma'
plt.rcParams['axes.linewidth'] = 2
# Set the default colour cycle (in case someone changes it...)
from cycler import cycler
cols = plt.get_cmap('tab10').colors
plt.rcParams['axes.prop_cycle'] = cycler(color=cols)

# A simple little 2D matrix plotter, excluding x and y labels.
def plot_2d(m, title=""):
    plt.imshow(m)
    plt.xticks([])
    plt.yticks([])
    plt.title(title)

'''
                                input tomeseries
'''
import csv
csv_path1 = "C:\\Users\\ZHOU\\Desktop\\成果\\rainata.csv"
csv_path2 = "C:\\Users\\ZHOU\\Desktop\\成果\\time.csv"
with open(csv_path1,'r',encoding='utf-8-sig')as fp:
    # 使用列表推导式，将读取到的数据装进列表
    data_rows1 = [i[0] for i in csv.reader(fp)]  # csv.reader 读取到的数据是list类型
with open(csv_path2,'r',encoding='utf8')as p:
    t_rows1 = [j[0] for j in csv.reader(p)]
#data_rows for data
list1 = data_rows1 
list2 = t_rows1
#F IS DATA

F = []
t = []

for i in list1:
     i =float(i)
     F.append(i)
for j in list2:
     j =float(j)
     t.append(j)
N = len(t) # The number of time 'moments' in our  series

# Plot everything
plt.plot(t, F, lw=2.5)
plt.legend(["GWS Series", "Trend", "Periodic #1", "Periodic #2", "Noise"])
plt.xlabel("Time")
plt.ylabel("Aonmaly Value(cm)")
plt.title("GWS Time Series");
plt.show()

'''
From Time Series to Trajectory Matrix
'''

L = 80 # The window length.
K = N - L + 1 # The number of columns in the trajectory matrix.
# Create the trajectory matrix by pulling the relevant subseries of F, and stacking them as columns.
X = np.column_stack([F[i:i+L] for i in range(0,K)])
# Note: the i+L above gives us up to i+L-1, as numpy array upper bounds are exclusive. 

#plot Trajectory Matrix
ax = plt.matshow(X)
plt.xlabel("$L$-Lagged Vectors")
plt.ylabel("$K$-Lagged Vectors")
plt.colorbar(ax.colorbar, fraction=0.03,pad=0.04)
ax.colorbar.set_label("$F(t)$")
plt.title("The Trajectory Matrix for the Time Series");


'''
Decomposing the Trajectory Matrix

'''


#Let's take a peak at the first 12 elementary matrices.

d = np.linalg.matrix_rank(X) # The intrinsic dimensionality of the trajectory space.

# For those interested in how to code up an SVD calculation, Numerical Recipes in Fortran 77
# has you covered: http://www.aip.de/groups/soe/local/numres/bookfpdf/f2-6.pdf
# Thankfully, we'll leave the actual SVD calculation to NumPy.
U, Sigma, V = np.linalg.svd(X)
V = V.T # Note: the SVD routine returns V^T, not V, so I'll tranpose it back here. This may seem pointless, 
# but I'll treat the Python representation of V consistently with the mathematical notation in this notebook.

# Calculate the elementary matrices of X, storing them in a multidimensional NumPy array.
# This requires calculating sigma_i * U_i * (V_i)^T for each i, or sigma_i * outer_product(U_i, V_i). 
# Note that Sigma is a 1D array of singular values, instead of the full L x K diagonal matrix.
X_elem = np.array( [Sigma[i] * np.outer(U[:,i], V[:,i]) for i in range(0,d)] )

# Quick sanity check: the sum of all elementary matrices in X_elm should be equal to X, to within a 
# *very small* tolerance:
if not np.allclose(X, X_elem.sum(axis=0), atol=1e-10):
    print("WARNING: The sum of X's elementary matrices is not equal to X!")

#plot
n = min(20, d) # In case d is less than 12 for the toy series. Say, if we were to exclude the noise component...
for i in range(n):
    plt.subplot(4,3,i+1)
    title = "$\mathbf{X}_{" + str(i) + "}$"
    plot_2d(X_elem[i], title)
plt.tight_layout()


#Relative Contribution of Xi to Trajectory Matrix

sigma_sumsq = (Sigma**2).sum()
fig, ax = plt.subplots(1, 2, figsize=(14,5))
ax[0].plot(Sigma**2 / sigma_sumsq * 100, lw=2.5)
ax[0].set_xlim(0,12)
ax[0].set_title("Relative Contribution of $\mathbf{X}_i$ to Trajectory Matrix")
ax[0].set_xlabel("$i$")
ax[0].set_ylabel("Contribution (%)")
ax[1].plot((Sigma**2).cumsum() / sigma_sumsq * 100, lw=2.5)
ax[1].set_xlim(0,12)
ax[1].set_title("Cumulative Contribution of $\mathbf{X}_i$ to Trajectory Matrix")
ax[1].set_xlabel("$i$")
ax[1].set_ylabel("Contribution (%)");

'''
Reconstructing the Time Series

'''

#def Hankelisation 
def Hankelise(X):
    """
    Hankelises the matrix X, returning H(X).
    """
    L, K = X.shape
    transpose = False
    if L > K:
        # The Hankelisation below only works for matrices where L < K.
        # To Hankelise a L > K matrix, first swap L and K and tranpose X.
        # Set flag for HX to be transposed before returning. 
        X = X.T
        L, K = K, L
        transpose = True

    HX = np.zeros((L,K))
    
    # I know this isn't very efficient...
    for m in range(L):
        for n in range(K):
            s = m+n
            if 0 <= s <= L-1:
                for l in range(0,s+1):
                    HX[m,n] += 1/(s+1)*X[l, s-l]    
            elif L <= s <= K-1:
                for l in range(0,L-1):
                    HX[m,n] += 1/(L-1)*X[l, s-l]
            elif K <= s <= K+L-2:
                for l in range(s-K+1,L):
                    HX[m,n] += 1/(K+L-s-1)*X[l, s-l]
    if transpose:
        return HX.T
    else:
        return HX

n = min(d, 12)
for j in range(0,n):
    plt.subplot(4,3,j+1)
    title = r"$\tilde{\mathbf{X}}_{" + str(j) + "}$"
    plot_2d(Hankelise(X_elem[j]), title)
plt.tight_layout() 


def X_to_TS(X_i):
    """Averages the anti-diagonals of the given elementary matrix, X_i, and returns a time series."""
    # Reverse the column ordering of X_i
    X_rev = X_i[::-1]
    # Full credit to Mark Tolonen at https://stackoverflow.com/a/6313414 for this one:
    return np.array([X_rev.diagonal(i).mean() for i in range(-X_i.shape[0]+1, X_i.shape[1])])

'''
plot Fi

'''

n = min(11,d) # In case of noiseless time series with d < 12.

# Fiddle with colour cycle - need more colours!
fig = plt.subplot()
color_cycle = cycler(color=plt.get_cmap('tab20').colors)
fig.axes.set_prop_cycle(color_cycle)

# Convert elementary matrices straight to a time series - no need to construct any Hankel matrices.
for i in range(n):
    F_i = X_to_TS(X_elem[i])
    fig.axes.plot(t, F_i, lw=2)

fig.axes.plot(t, F, alpha=1, lw=1)
fig.set_xlabel("$t$")
fig.set_ylabel(r"$\tilde{F}_i(t)$")
legend = [r"$\tilde{F}_{%s}$" %i for i in range(n)] + ["$F$"]
fig.set_title("The First 12 Components of the  Time Series")
fig.legend(legend, loc=(1.05,0.1));


'''
展示12个分量
'''

fig,axes = plt.subplots(12,1,sharex = 'col')
fig.set_figwidth(20)
fig.set_figheight(40)

n = 20
for i in range(n):
    F_i = X_to_TS(X_elem[i])
    axes[0].set_title("The First 12 Components of the  Time Series")
    axes[i].plot(t, F_i, lw=1.5)
    axes[11].set_xlabel("$t$")
    axes[i].set_ylabel(r"$\tilde{F}_{%s}$" %i )
plt.savefig('前12分量.png',dpi=300) 
print(os.getcwd())


    

'''
Assemble the grouped components of the time series.

'''

# Assemble the grouped components of the time series.
F_trend = X_to_TS(X_elem[[0]].sum(axis=0))
F_periodic1 = X_to_TS(X_elem[[1,4,2,3]].sum(axis=0))
F_periodic2 = X_to_TS(X_elem[[8,9,10,11]].sum(axis=0))
F_periodic3 = X_to_TS(X_elem[[6,7]].sum(axis=0))
F_noise = X_to_TS(X_elem[11:].sum(axis=0))

# Plot the toy time series and its separated components on a single plot.
plt.plot(t,F, lw=1)
plt.plot(t, F_trend)
plt.plot(t, F_periodic1)
#plt.plot(t, F_periodic2)
plt.plot(t, F_periodic3)
plt.plot(t, F_noise, alpha=0.5)
plt.xlabel("$t$")
plt.ylabel(r"$\tilde{F}^{(j)}$")
groups = ["trend", "periodic 1", "periodic 2","periodic 3", "noise"]
legend = ["$F$"] + [r"$\tilde{F}^{(\mathrm{%s})}$"%group for group in groups]
plt.legend(legend)
plt.title("Grouped Time Series Components")
plt.show()

# A list of tuples so we can create the next plot with a loop.
components = [("Original",  F),
              ("Trend",  F_trend), 
              ("Periodic1",  F_periodic1),
              ]

# Plot the separated components and original components together.
fig = plt.figure()
fig.set_figwidth(15)
fig.set_figheight(20)
n=1
for name, ssa_comp in components:
    ax = fig.add_subplot(4,1,n)
    # ax.plot(t, orig_comp, linestyle="--", lw=2.5, alpha=0.7) 
    ax.plot(t, ssa_comp,lw=2.5,alpha=1)
    ax.set_title(name, fontsize=16)
    ax.set_xlabel("$t$")
    ax.set_ylabel("cm")
    n += 1

fig.tight_layout()
plt.savefig('前12分量.png',dpi=300) 
print(os.getcwd())


'''
construct the w-correlation matrix for the time series
'''

# Get the weights w first, as they'll be reused a lot.
# Note: list(np.arange(L)+1) returns the sequence 1 to L (first line in definition of w), 
# [L]*(K-L-1) repeats L K-L-1 times (second line in w definition)
# list(np.arange(L)+1)[::-1] reverses the first list (equivalent to the third line)
# Add all the lists together and we have our array of weights.
w = np.array(list(np.arange(L)+1) + [L]*(K-L-1) + list(np.arange(L)+1)[::-1])

# Get all the components of the toy series, store them as columns in F_elem array.
F_elem = np.array([X_to_TS(X_elem[i]) for i in range(d)])

# Calculate the individual weighted norms, ||F_i||_w, first, then take inverse square-root so we don't have to later.
F_wnorms = np.array([w.dot(F_elem[i]**2) for i in range(d)])
F_wnorms = F_wnorms**-0.5

# Calculate the w-corr matrix. The diagonal elements are equal to 1, so we can start with an identity matrix
# and iterate over all pairs of i's and j's (i != j), noting that Wij = Wji.
Wcorr = np.identity(d)
for i in range(d):
    for j in range(i+1,d):
        Wcorr[i,j] = abs(w.dot(F_elem[i]*F_elem[j]) * F_wnorms[i] * F_wnorms[j])
        Wcorr[j,i] = Wcorr[i,j]
        
#Plot all the w-correlation matrix.

ax = plt.imshow(Wcorr)
plt.xlabel(r"$\tilde{F}_i$")
plt.ylabel(r"$\tilde{F}_j$")
plt.colorbar(ax.colorbar, fraction=0.045)
ax.colorbar.set_label("$W_{ij}$")
plt.clim(0,1)
plt.title("The W-Correlation Matrix for the Time Series");

#plot 0-6  wrr 
ax = plt.imshow(Wcorr)
plt.xlabel(r"$\tilde{F}_i$")
plt.ylabel(r"$\tilde{F}_j$")
plt.colorbar(ax.colorbar, fraction=0.045)
ax.colorbar.set_label("$W_{ij}$")
plt.xlim(-0.5,12.5)   #此处为plot窗口大小
plt.ylim(12.5,-0.5)
plt.clim(0,1)
plt.title(r"W-Correlation for Components 0–11");