from __future__ import print_function
from sklearn import datasets

import matplotlib.pyplot as plt
import matplotlib.cm as cmx
import matplotlib.colors as colors
import numpy as np

%matplotlib inline

def shuffle_data(X, y, seed=None):
   if seed:
     np.random.seed(seed)

   idx = np.arange(X.shape[0])
   np.random.shuffle(idx)
   return X[idx], y[idx]

# 正規(guī)化數(shù)據(jù)集 X

def normalize(X, axis=-1, p=2):
   lp_norm = np.atleast_1d(np.linalg.norm(X, p, axis))
   lp_norm[lp_norm == 0] = 1
   return X / np.expand_dims(lp_norm, axis)

# 標(biāo)準(zhǔn)化數(shù)據(jù)集 X

def standardize(X):
   X_std = np.zeros(X.shape)
   mean = X.mean(axis=0)
   std = X.std(axis=0)

   # 做除法運(yùn)算時(shí)請(qǐng)永遠(yuǎn)記住分母不能等于 0 的情形
   # X_std = (X - X.mean(axis=0)) / X.std(axis=0) 

   for col in range(np.shape(X)[1]):
     if std[col]:
       X_std[:, col] = (X_std[:, col] - mean[col]) / std[col]
   return X_std

# 劃分?jǐn)?shù)據(jù)集為訓(xùn)練集和測(cè)試集

def train_test_split(X, y, test_size=0.2, shuffle=True, seed=None):
   if shuffle:
     X, y = shuffle_data(X, y, seed)
   n_train_samples = int(X.shape[0] * (1-test_size))
   x_train, x_test = X[:n_train_samples], X[n_train_samples:]
   y_train, y_test = y[:n_train_samples], y[n_train_samples:]

   return x_train, x_test, y_train, y_test

# 計(jì)算矩陣 X 的協(xié)方差矩陣

def calculate_covariance_matrix(X, Y=np.empty((0,0))):
   if not Y.any():
      Y = X
   n_samples = np.shape(X)[0]
   covariance_matrix = (1 / (n_samples-1)) * (X - X.mean(axis=0)).T.dot(Y - Y.mean(axis=0))
   return np.array(covariance_matrix, dtype=float)

# 計(jì)算數(shù)據(jù)集 X 每列的方差

def calculate_variance(X):
   n_samples = np.shape(X)[0]
   variance = (1 / n_samples) * np.diag((X - X.mean(axis=0)).T.dot(X - X.mean(axis=0)))
   return variance

# 計(jì)算數(shù)據(jù)集 X 每列的標(biāo)準(zhǔn)差

def calculate_std_dev(X):
   std_dev = np.sqrt(calculate_variance(X))
   return std_dev

# 計(jì)算相關(guān)系數(shù)矩陣

def calculate_correlation_matrix(X, Y=np.empty([0])):
   # 先計(jì)算協(xié)方差矩陣
   covariance_matrix = calculate_covariance_matrix(X, Y)
   # 計(jì)算 X, Y 的標(biāo)準(zhǔn)差
   std_dev_X = np.expand_dims(calculate_std_dev(X), 1)

   std_dev_y = np.expand_dims(calculate_std_dev(Y), 1)

   correlation_matrix = np.divide(covariance_matrix, std_dev_X.dot(std_dev_y.T))
   return np.array(correlation_matrix, dtype=float)

class PCA():

   """

   主成份分析算法 PCA，非監(jiān)督學(xué)習(xí)算法.

   """
   def __init__(self):
     self.eigen_values = None
     self.eigen_vectors = None
     self.k = 2

   def transform(self, X):

     """ 

     將原始數(shù)據(jù)集 X 通過 PCA 進(jìn)行降維

     """

     covariance = calculate_covariance_matrix(X)
     # 求解特征值和特征向量
     self.eigen_values, self.eigen_vectors = np.linalg.eig(covariance)

     # 將特征值從大到小進(jìn)行排序，注意特征向量是按列排的，即 self.eigen_vectors 第 k 列是 self.eigen_values 中第 k 個(gè)特征值對(duì)應(yīng)的特征向量
     idx = self.eigen_values.argsort()[::-1]
     eigenvalues = self.eigen_values[idx][:self.k]
     eigenvectors = self.eigen_vectors[:, idx][:, :self.k]

     # 將原始數(shù)據(jù)集 X 映射到低維空間
     X_transformed = X.dot(eigenvectors)
     return X_transformed

def main():
   # Load the dataset
   data = datasets.load_iris()
   X = data.data
   y = data.target

   # 將數(shù)據(jù)集 X 映射到低維空間
   X_trans = PCA().transform(X)

   x1 = X_trans[:, 0]

   x2 = X_trans[:, 1]
   cmap = plt.get_cmap('viridis')
   colors = [cmap(i) for i in np.linspace(0, 1, len(np.unique(y)))]
   class_distr = []

   # Plot the different class distributions
   for i, l in enumerate(np.unique(y)):
       _x1 = x1[y == l]
       _x2 = x2[y == l]
       _y = y[y == l]
       class_distr.append(plt.scatter(_x1, _x2, color=colors[i]))
   # Add a legend
   plt.legend(class_distr, y, loc=1)

   # Axis labels
   plt.xlabel('Principal Component 1')
   plt.ylabel('Principal Component 2')
   plt.show()

if __name__ == "__main__":

   main()

# coding:utf-8
# 實(shí)現(xiàn)KPCA

from sklearn.datasets import load_iris
from sklearn.decomposition import KernelPCA
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import pdist, squareform

'''
author: heucoder
email: 812860165@qq.com
date: 2019.6.13
'''


def sigmoid(x, coef = 0.25):
    x = np.dot(x, x.T)
    return np.tanh(coef*x+1)

def linear(x):
    x = np.dot(x, x.T)
    return x

def rbf(x, gamma = 15):
    sq_dists = pdist(x, 'sqeuclidean')
    mat_sq_dists = squareform(sq_dists)
    return np.exp(-gamma*mat_sq_dists)

def kpca(data, n_dims=2, kernel = rbf):
    '''
    :param data: (n_samples, n_features)
    :param n_dims: target n_dims
    :param kernel: kernel functions
    :return: (n_samples, n_dims)
    '''

    K = kernel(data)
    #
    N = K.shape[0]
    one_n = np.ones((N, N)) / N
    K = K - one_n.dot(K) - K.dot(one_n) + one_n.dot(K).dot(one_n)
    #
    eig_values, eig_vector = np.linalg.eig(K)
    idx = eig_values.argsort()[::-1]
    eigval = eig_values[idx][:n_dims]
    eigvector = eig_vector[:, idx][:, :n_dims]
    print(eigval)
    eigval = eigval**(1/2)
    vi = eigvector/eigval.reshape(-1,n_dims)
    data_n = np.dot(K, vi)
    return data_n


if __name__ == "__main__":
    data = load_iris().data
    Y = load_iris().target
    data_1 = kpca(data, kernel=rbf)


    sklearn_kpca = KernelPCA(n_components=2, kernel="rbf", gamma=15)
    data_2 = sklearn_kpca.fit_transform(data)

    plt.figure(figsize=(8,4))
    plt.subplot(121)
    plt.title("my_KPCA")
    plt.scatter(data_1[:, 0], data_1[:, 1], c = Y)

    plt.subplot(122)
    plt.title("sklearn_KPCA")
    plt.scatter(data_2[:, 0], data_2[:, 1], c = Y)
    plt.show()

#coding:utf-8
import numpy as np
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt

'''
author: heucoder
email: 812860165@qq.com
date: 2019.6.13
'''


def lda(data, target, n_dim):
    '''
    :param data: (n_samples, n_features)
    :param target: data class
    :param n_dim: target dimension
    :return: (n_samples, n_dims)
    '''

    clusters = np.unique(target)

    if n_dim > len(clusters)-1:
        print("K is too much")
        print("please input again")
        exit(0)

    #within_class scatter matrix
    Sw = np.zeros((data.shape[1],data.shape[1]))
    for i in clusters:
        datai = data[target == i]
        datai = datai-datai.mean(0)
        Swi = np.mat(datai).T*np.mat(datai)
        Sw += Swi

    #between_class scatter matrix
    SB = np.zeros((data.shape[1],data.shape[1]))
    u = data.mean(0)  #所有樣本的平均值
    for i in clusters:
        Ni = data[target == i].shape[0]
        ui = data[target == i].mean(0)  #某個(gè)類別的平均值
        SBi = Ni*np.mat(ui - u).T*np.mat(ui - u)
        SB += SBi
    S = np.linalg.inv(Sw)*SB
    eigVals,eigVects = np.linalg.eig(S)  #求特征值，特征向量
    eigValInd = np.argsort(eigVals)
    eigValInd = eigValInd[:(-n_dim-1):-1]
    w = eigVects[:,eigValInd]
    data_ndim = np.dot(data, w)

    return data_ndim

if __name__ == '__main__':
    iris = load_iris()
    X = iris.data
    Y = iris.target
    data_1 = lda(X, Y, 2)

    data_2 = LinearDiscriminantAnalysis(n_components=2).fit_transform(X, Y)


    plt.figure(figsize=(8,4))
    plt.subplot(121)
    plt.title("my_LDA")
    plt.scatter(data_1[:, 0], data_1[:, 1], c = Y)

    plt.subplot(122)
    plt.title("sklearn_LDA")
    plt.scatter(data_2[:, 0], data_2[:, 1], c = Y)
    plt.savefig("LDA.png")
    plt.show()

# coding:utf-8
import numpy as np
from sklearn.datasets import make_s_curve
import matplotlib.pyplot as plt
from sklearn.manifold import Isomap
from mpl_toolkits.mplot3d import Axes3D

def floyd(D,n_neighbors=15):
    Max = np.max(D)*1000
    n1,n2 = D.shape
    k = n_neighbors
    D1 = np.ones((n1,n1))*Max
    D_arg = np.argsort(D,axis=1)
    for i in range(n1):
        D1[i,D_arg[i,0:k+1]] = D[i,D_arg[i,0:k+1]]
    for k in range(n1):
        for i in range(n1):
            for j in range(n1):
                if D1[i,k]+D1[k,j]<D1[i,j]:
                    D1[i,j] = D1[i,k]+D1[k,j]
    return D1

def cal_pairwise_dist(x):
    '''計(jì)算pairwise 距離, x是matrix
    (a-b)^2 = a^2 + b^2 - 2*a*b
    '''
    sum_x = np.sum(np.square(x), 1)
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    #返回任意兩個(gè)點(diǎn)之間距離的平方
    return dist

def my_mds(dist, n_dims):
    # dist (n_samples, n_samples)
    dist = dist**2
    n = dist.shape[0]
    T1 = np.ones((n,n))*np.sum(dist)/n**2
    T2 = np.sum(dist, axis = 1)/n
    T3 = np.sum(dist, axis = 0)/n

    B = -(T1 - T2 - T3 + dist)/2

    eig_val, eig_vector = np.linalg.eig(B)
    index_ = np.argsort(-eig_val)[:n_dims]
    picked_eig_val = eig_val[index_].real
    picked_eig_vector = eig_vector[:, index_]

    return picked_eig_vector*picked_eig_val**(0.5)

def my_Isomap(data,n=2,n_neighbors=30):
    D = cal_pairwise_dist(data)
    D[D < 0] = 0
    D = D**0.5
    D_floyd=floyd(D, n_neighbors)
    data_n = my_mds(D_floyd, n_dims=n)
    return data_n

def scatter_3d(X, y):
    fig = plt.figure(figsize=(6, 5))
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y, cmap=plt.cm.hot)
    ax.view_init(10, -70)
    ax.set_xlabel("$x_1$", fontsize=18)
    ax.set_ylabel("$x_2$", fontsize=18)
    ax.set_zlabel("$x_3$", fontsize=18)
    plt.show()

if __name__ == '__main__':
    X, Y = make_s_curve(n_samples = 500,
                           noise = 0.1,
                           random_state = 42)

    data_1 = my_Isomap(X, 2, 10)

    data_2 = Isomap(n_neighbors = 10, n_components = 2).fit_transform(X)

    plt.figure(figsize=(8,4))
    plt.subplot(121)
    plt.title("my_Isomap")
    plt.scatter(data_1[:, 0], data_1[:, 1], c = Y)

    plt.subplot(122)
    plt.title("sklearn_Isomap")
    plt.scatter(data_2[:, 0], data_2[:, 1], c = Y)
    plt.savefig("Isomap1.png")
    plt.show()

# coding:utf-8
import numpy as np
from sklearn.datasets import make_s_curve
import matplotlib.pyplot as plt
from sklearn.manifold import LocallyLinearEmbedding
from mpl_toolkits.mplot3d import Axes3D

'''
author: heucoder
email: 812860165@qq.com
date: 2019.6.13
'''

def make_swiss_roll(n_samples=100, noise=0.0, random_state=None):
    #Generate a swiss roll dataset.
    t = 1.5 * np.pi * (1 + 2 * np.random.rand(1, n_samples))
    x = t * np.cos(t)
    y = 83 * np.random.rand(1, n_samples)
    z = t * np.sin(t)
    X = np.concatenate((x, y, z))
    X += noise * np.random.randn(3, n_samples)
    X = X.T
    t = np.squeeze(t)
    return X, t

def cal_pairwise_dist(x):
    '''計(jì)算pairwise 距離, x是matrix
    (a-b)^2 = a^2 + b^2 - 2*a*b
    '''
    sum_x = np.sum(np.square(x), 1)
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    #返回任意兩個(gè)點(diǎn)之間距離的平方
    return dist

def get_n_neighbors(data, n_neighbors = 10):
    '''
    :param data: (n_samples, n_features)
    :param n_neighbors: n nearest neighbors
    :return: neighbors indexs
    '''

    dist = cal_pairwise_dist(data)
    dist[dist < 0] = 0
    dist = dist**0.5
    n = dist.shape[0]
    N = np.zeros((n, n_neighbors))

    for i in range(n):
        index_ = np.argsort(dist[i])[1:n_neighbors+1]
        N[i] = N[i] + index_

    return N.astype(np.int32)

def lle(data, n_dims = 2, n_neighbors = 10):
    '''
    :param data:(n_samples, n_features)
    :param n_dims: target n_dims
    :param n_neighbors: n nearest neighbors
    :return: (n_samples, n_dims)
    '''
    N = get_n_neighbors(data, n_neighbors)
    n, D = data.shape

    # prevent Si to small
    if n_neighbors > D:
        tol = 1e-3
    else:
        tol = 0

    # calculate W
    W = np.zeros((n_neighbors, n))
    I = np.ones((n_neighbors, 1))
    for i in range(n):
        Xi = np.tile(data[i], (n_neighbors, 1)).T
        Ni = data[N[i]].T

        Si = np.dot((Xi-Ni).T, (Xi-Ni))
        # magic and why????
        Si = Si+np.eye(n_neighbors)*tol*np.trace(Si)

        Si_inv = np.linalg.pinv(Si)
        wi = (np.dot(Si_inv, I))/(np.dot(np.dot(I.T, Si_inv), I)[0,0])
        W[:, i] = wi[:,0]

    print("Xi.shape", Xi.shape)
    print("Ni.shape", Ni.shape)
    print("Si.shape", Si.shape)

    W_y = np.zeros((n, n))
    for i in range(n):
        index = N[i]
        for j in range(n_neighbors):
            W_y[index[j],i] = W[j,i]

    I_y = np.eye(n)
    M = np.dot((I_y - W_y), (I_y - W_y).T)

    eig_val, eig_vector = np.linalg.eig(M)
    index_ = np.argsort(np.abs(eig_val))[1:n_dims+1]
    print("index_", index_)
    Y = eig_vector[:, index_]
    return Y

if __name__ == '__main__':
    # X, Y = make_s_curve(n_samples = 500,
    #                        noise = 0.1,
    #                        random_state = 42)
    X, Y = make_swiss_roll(n_samples = 500, noise=0.1, random_state=42)
    
    data_1 =lle(X, n_neighbors = 30)
    print(data_1.shape)

    data_2 = LocallyLinearEmbedding(n_components=2, n_neighbors = 30).fit_transform(X)


    plt.figure(figsize=(8,4))
    plt.subplot(121)
    plt.title("my_LLE")
    plt.scatter(data_1[:, 0], data_1[:, 1], c = Y)

    plt.subplot(122)
    plt.title("sklearn_LLE")
    plt.scatter(data_2[:, 0], data_2[:, 1], c = Y)
    plt.savefig("LLE.png")
    plt.show()

# coding:utf-8

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from mpl_toolkits.mplot3d import Axes3D

'''
author: heucoder
email: 812860165@qq.com
date: 2019.6.13
'''

def make_swiss_roll(n_samples=100, noise=0.0, random_state=None):
    #Generate a swiss roll dataset.
    t = 1.5 * np.pi * (1 + 2 * np.random.rand(1, n_samples))
    x = t * np.cos(t)
    y = 83 * np.random.rand(1, n_samples)
    z = t * np.sin(t)
    X = np.concatenate((x, y, z))
    X += noise * np.random.randn(3, n_samples)
    X = X.T
    t = np.squeeze(t)
    return X, t

def rbf(dist, t = 1.0):
    '''
    rbf kernel function
    '''
    return np.exp(-(dist/t))

def cal_pairwise_dist(x):

    '''計(jì)算pairwise 距離, x是matrix
    (a-b)^2 = a^2 + b^2 - 2*a*b
    '''
    sum_x = np.sum(np.square(x), 1)
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    #返回任意兩個(gè)點(diǎn)之間距離的平方
    return dist

def cal_rbf_dist(data, n_neighbors = 10, t = 1):

    dist = cal_pairwise_dist(data)
    dist[dist < 0] = 0
    n = dist.shape[0]
    rbf_dist = rbf(dist, t)

    W = np.zeros((n, n))
    for i in range(n):
        index_ = np.argsort(dist[i])[1:1+n_neighbors]
        W[i, index_] = rbf_dist[i, index_]
        W[index_, i] = rbf_dist[index_, i]

    return W

def le(data,
          n_dims = 2,
          n_neighbors = 5, t = 1.0):
    '''
    :param data: (n_samples, n_features)
    :param n_dims: target dim
    :param n_neighbors: k nearest neighbors
    :param t: a param for rbf
    :return:
    '''
    N = data.shape[0]
    W = cal_rbf_dist(data, n_neighbors, t)
    D = np.zeros_like(W)
    for i in range(N):
        D[i,i] = np.sum(W[i])

    D_inv = np.linalg.inv(D)
    L = D - W
    eig_val, eig_vec = np.linalg.eig(np.dot(D_inv, L))

    sort_index_ = np.argsort(eig_val)

    eig_val = eig_val[sort_index_]
    print("eig_val[:10]: ", eig_val[:10])

    j = 0
    while eig_val[j] < 1e-6:
        j+=1

    print("j: ", j)

    sort_index_ = sort_index_[j:j+n_dims]
    eig_val_picked = eig_val[j:j+n_dims]
    print(eig_val_picked)
    eig_vec_picked = eig_vec[:, sort_index_]

    # print("L: ")
    # print(np.dot(np.dot(eig_vec_picked.T, L), eig_vec_picked))
    # print("D: ")
    # D not equal I ???
    print(np.dot(np.dot(eig_vec_picked.T, D), eig_vec_picked))

    X_ndim = eig_vec_picked
    return X_ndim

if __name__ == '__main__':
    # X, Y = make_swiss_roll(n_samples = 2000)
    # X_ndim = le(X, n_neighbors = 5, t = 20)
    #
    # fig = plt.figure(figsize=(12,6))
    # ax1 = fig.add_subplot(121, projection='3d')
    # ax1.scatter(X[:, 0], X[:, 1], X[:, 2], c = Y)
    #
    # ax2 = fig.add_subplot(122)
    # ax2.scatter(X_ndim[:, 0], X_ndim[:, 1], c = Y)
    # plt.show()

    X = load_digits().data
    y = load_digits().target

    dist = cal_pairwise_dist(X)
    max_dist = np.max(dist)
    print("max_dist", max_dist)
    X_ndim = le(X, n_neighbors = 20, t = max_dist*0.1)
    plt.scatter(X_ndim[:, 0], X_ndim[:, 1], c = y)
    plt.savefig("LE2.png")
    plt.show()

# coding:utf-8

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D
from sklearn.datasets import load_digits, load_iris

'''
author: heucoder
email: 812860165@qq.com
date: 2019.6.13
'''

def make_swiss_roll(n_samples=100, noise=0.0, random_state=None):
    #Generate a swiss roll dataset.
    t = 1.5 * np.pi * (1 + 2 * np.random.rand(1, n_samples))
    x = t * np.cos(t)
    y = 83 * np.random.rand(1, n_samples)
    z = t * np.sin(t)
    X = np.concatenate((x, y, z))
    X += noise * np.random.randn(3, n_samples)
    X = X.T
    t = np.squeeze(t)
    return X, t

def rbf(dist, t = 1.0):
    '''
    rbf kernel function
    '''
    return np.exp(-(dist/t))

def cal_pairwise_dist(x):

    '''計(jì)算pairwise 距離, x是matrix
    (a-b)^2 = a^2 + b^2 - 2*a*b
    '''
    sum_x = np.sum(np.square(x), 1)
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    #返回任意兩個(gè)點(diǎn)之間距離的平方
    return dist

def cal_rbf_dist(data, n_neighbors = 10, t = 1):

    dist = cal_pairwise_dist(data)
    dist[dist < 0] = 0
    n = dist.shape[0]
    rbf_dist = rbf(dist, t)

    W = np.zeros((n, n))
    for i in range(n):
        index_ = np.argsort(dist[i])[1:1 + n_neighbors]
        W[i, index_] = rbf_dist[i, index_]
        W[index_, i] = rbf_dist[index_, i]

    return W

def lpp(data,
        n_dims = 2,
        n_neighbors = 30, t = 1.0):
    '''
    :param data: (n_samples, n_features)
    :param n_dims: target dim
    :param n_neighbors: k nearest neighbors
    :param t: a param for rbf
    :return:
    '''
    N = data.shape[0]
    W = cal_rbf_dist(data, n_neighbors, t)
    D = np.zeros_like(W)

    for i in range(N):
        D[i,i] = np.sum(W[i])

    L = D - W
    XDXT = np.dot(np.dot(data.T, D), data)
    XLXT = np.dot(np.dot(data.T, L), data)

    eig_val, eig_vec = np.linalg.eig(np.dot(np.linalg.pinv(XDXT), XLXT))

    sort_index_ = np.argsort(np.abs(eig_val))
    eig_val = eig_val[sort_index_]
    print("eig_val[:10]", eig_val[:10])

    j = 0
    while eig_val[j] < 1e-6:
        j+=1

    print("j: ", j)

    sort_index_ = sort_index_[j:j+n_dims]
    # print(sort_index_)
    eig_val_picked = eig_val[j:j+n_dims]
    print(eig_val_picked)
    eig_vec_picked = eig_vec[:, sort_index_]

    data_ndim = np.dot(data, eig_vec_picked)

    return data_ndim

if __name__ == '__main__':
    X = load_digits().data
    y = load_digits().target
    # X, y = make_swiss_roll(n_samples = 1000)

    dist = cal_pairwise_dist(X)
    max_dist = np.max(dist)
    print("max_dist", max_dist)

    data_2d = lpp(X, n_neighbors = 5, t = 0.01*max_dist)
    data_2 = PCA(n_components=2).fit_transform(X)


    plt.figure(figsize=(12,6))
    plt.subplot(121)
    plt.title("LPP")
    plt.scatter(data_2d[:, 0], data_2d[:, 1], c = y)

    plt.subplot(122)
    plt.title("PCA")
    plt.scatter(data_2[:, 0], data_2[:, 1], c = y)
    plt.show()