谱聚类的代码 测试版

import time
# 导入库函数
import numpy as np
import pandas as pd
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn import cluster, datasets, mixture
from sklearn.neighbors import kneighbors_graph
from sklearn.preprocessing import StandardScaler, LabelEncoder
from itertools import cycle, islice

import networkx as nx
from scipy import sparse
from sklearn.cluster import KMeans
import torch
np.random.seed(0)


def getClusterCentroids(X, spectral_labels):
    """
    Funcao auxiliar para obter os centroids dos clusters a partir dos dados X e das marcacoes de spectral_labels
    """
    tmp = pd.DataFrame(X)
    cols = tmp.columns
    tmp['spectral_labels'] = spectral_labels
    return tmp.groupby("spectral_labels")[cols].mean().values


def laplacian_matrix_old(G, nodelist=None, weight="weight"):
#未修改的拉普拉斯矩阵计算
    import scipy as sp
    import scipy.sparse  # call as sp.sparse

    if nodelist is None:
        nodelist = list(G)
    A = nx.to_scipy_sparse_array(G, nodelist=nodelist, weight=weight, format="csr")#A就是距离的矩阵
    a=A.toarray()
    print(a)
    n, m = A.shape
    # TODO: rm csr_array wrapper when spdiags can produce arrays
    D = sp.sparse.csr_array(sp.sparse.spdiags(A.sum(axis=1), 0, m, n, format="csr"))#D是权重矩阵
    d=D.toarray()
    print(d)
    diags = A.sum(axis=1)
    with sp.errstate(divide="ignore"):
        diags_sqrt = 1.0 / np.sqrt(diags)
    diags_sqrt[np.isinf(diags_sqrt)] = 0
    # TODO: rm csr_array wrapper when spdiags can produce arrays
    DH = sp.sparse.csr_array(sp.sparse.spdiags(diags_sqrt, 0, m, n, format="csr"))
    I = np.identity(1500)
    L=D-A
    equal1 = I - DH @ (A @ DH)
    print('DH @ (A @ DH):     ')
    print(DH @ (A @ DH))
    #equal11=equal1.toarray()
    print('equal1:   ')
    print(equal1)
    equal2 = DH @ (L @ DH)
    equal2 = equal2.toarray()
    print('equal2:   ')
    print(equal2)
    l=L.toarray()
    Dni = np.linalg.inv(a)
    equal3=Dni @ L
    print('equal3:      ')
    print(equal3)
    #equal3 = equal3.toarray()
    equal4 = I - Dni @ A
    equal2_L2 = np.linalg.norm(equal2,ord=2)
    print('equal2_L2:      ')
    print(equal2_L2)
    equal3_L2 = np.linalg.norm(equal3, ord=2)
    print('equal3_L2:      ')
    print(equal3_L2)
    print('正则化的差：  ')
    print(equal3_L2-equal2_L2)
    #print(cha)
    print((equal3 == equal2).all())
    return D - A

def laplacian_matrix_new(G, nodelist=None, weight="weight"):
#把权重矩阵和距离矩阵都倒数的拉普拉斯矩阵计算
    import scipy as sp
    import scipy.sparse  # call as sp.sparse

    if nodelist is None:
        nodelist = list(G)
    A = nx.to_scipy_sparse_array(G, nodelist=nodelist, weight=weight, format="csr")
    n,m=A.shape
    print(n,m)

    n, m = A.shape
    # TODO: rm csr_array wrapper when spdiags can produce arrays
    D = sp.sparse.csr_array(sp.sparse.spdiags(A.sum(axis=1), 0, m, n, format="csr"))

    a= A.toarray()
    print(a)
    a = 1/(a+10**(-10))
    d = D.toarray()
    print(d)
    diags = a.sum(axis=1)
    with sp.errstate(divide="ignore"):
        diags_sqrt = 1.0 / np.sqrt(diags)
    diags_sqrt[np.isinf(diags_sqrt)] = 0
    # TODO: rm csr_array wrapper when spdiags can produce arrays
    DH = sp.sparse.csr_array(sp.sparse.spdiags(diags_sqrt, 0, m, n, format="csr"))
    l = d-a
    I=np.identity(1500)
    equal1 = I-DH @ (A @ DH)
    L = sp.sparse.csr_matrix(l)
    equal2 = DH @ (L @ DH)
    Dni = np.linalg.inv(a)
    equal3 = Dni @ L
    equal4 = I - Dni @ A

    print((equal3==equal2).all())
    return L,DH
def laplacian_matrix_newest(G, nodelist=None, weight="weight"):
#错误的矩阵，只倒数权重的对角线 运行不出来
    import scipy as sp
    import scipy.sparse  # call as sp.sparse

    if nodelist is None:
        nodelist = list(G)
    A = nx.to_scipy_sparse_array(G, nodelist=nodelist, weight=weight, format="csr")
    n, m = A.shape
    # TODO: rm csr_array wrapper when spdiags can produce arrays
    D = sp.sparse.csr_array(sp.sparse.spdiags(A.sum(axis=1), 0, m, n, format="csr"))
    a= A.toarray()
    a = 1/(a+10**(-10))
    d = D.toarray()
    print(d)
    for i in range(n):
        d[i][i]=1/d[i][i]
    print(d)
    l = d-a
    L = sp.sparse.csr_matrix(l)
    return L

def unnormalizedSpectralClustering(X, k, params):
    n = len(X)
    ## Construct a similarity graph by one of the ways described in Section 2. Let W be its weighted adjacency matrix
    # Computes the (weighted) graph of k-Neighbors for points in X. The default distance is 'euclidean'
    A = kneighbors_graph(X, params['k_neighbors'], mode='distance', metric='euclidean', include_self=True)
    G = nx.from_scipy_sparse_array(A)
    L ,DH= laplacian_matrix_old(G)
#只需要更换这个函数
    print(DH)

    #equal=DH@(L@DH)
    #print(equal!=L)

    print(L)

    print('hello')

    eigenvalues, eigenvectors = sparse.linalg.eigs(L, k=k, which='SM')
    eigenvectors = np.real_if_close(eigenvectors)
    eigenvalues = np.real_if_close(eigenvalues)
    kfirst_indices = np.argsort(eigenvalues)[:k]

    ## Let U be the matrix containing the vectors u_1,...,u_k as columns
    ## For i=1,...,n, let y_i be the vector corresponding to the i-th row of U
    Y = eigenvectors[:, kfirst_indices]

    ## Cluster the points (y_i)i=1,...,n with the k-means algorithm into clusters C_1,...,C_k
    kmeans = KMeans(n_clusters=k, random_state=0).fit(Y)

    ## Output: Clusters A_1,...,A_k with A_i = {j|y_j in C_i}
    return {
        'labels': kmeans.labels_,
        'centroids': getClusterCentroids(X, kmeans.labels_)
    }

n_samples = 1500
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5,
                                      noise=.05)
noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
no_structure = np.random.rand(n_samples, 2), None

simple_datasets = [
    (noisy_circles, {'name': 'Noisy Circles','n_clusters': 2}),
    (noisy_moons, {'name': 'Noisy Moons', 'n_clusters': 2}),
    (no_structure, {'name': 'No structure', 'n_clusters': 3})]

# ============
# Set up cluster parameters
# ============
plt.figure(figsize=(4 * 2 + 3, 12.5))
plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05,
                    hspace=.01)
plot_num = 1
for i_dataset, (dataset, params) in enumerate(simple_datasets):
    X, y = dataset
    # normalize dataset for easier parameter selection
    X = StandardScaler().fit_transform(X)
    kmeans = cluster.KMeans(n_clusters=params['n_clusters'])
    clustering_algorithms = (
        ('KMeans', kmeans),
        ('Unnormalized Spectral Clustering', unnormalizedSpectralClustering),
    )
    for name, algorithm in clustering_algorithms:
        t0 = time.time()
        if name == 'KMeans':
            algorithm.fit(X)
        else:
            k = params['n_clusters']
            spectral_params = {
                'k_neighbors': 12
            }
            result = algorithm(X, k, spectral_params)

        t1 = time.time()

        if name == 'KMeans':
            y_pred = algorithm.labels_.astype(int)
        else:
            y_pred = result['labels']

        plt.subplot(len(simple_datasets), len(clustering_algorithms), plot_num)
        if i_dataset == 0:
            plt.title(name)

        colors = np.array(list(islice(cycle(['#377eb8', '#ff7f00', '#4daf4a',
                                             '#f781bf', '#a65628', '#984ea3',
                                             '#999999', '#e41a1c', '#dede00']),
                                      int(max(y_pred) + 1))))

        plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])

        plt.xlim(-2.5, 2.5)
        plt.ylim(-2.5, 2.5)
        plt.xticks(())
        plt.yticks(())
        plt.text(.99, .01, ('%.2fs' % (t1 - t0)).lstrip('0'),
                 transform=plt.gca().transAxes, size=15,
                 horizontalalignment='right')
        plot_num += 1

plt.show()