通过快速搜索和寻找密度峰值进行聚类

你好，欢迎你来看这篇小文章！
本篇的内容是对2014年发表在science上的论文Clustering by fast search and find of density peaks论文的概述和用Python代码的实现，希望能对你有所帮助。本篇文章会做到尽量简洁。

关于论文

本篇将不会对论文原理和具体内容进行具体描述，附上原文链接
link.
这篇文章的要点内容我将其整理为了一张思维导图，如下：

关于代码

本代码使用的数据共有三列，前两列分别是1号点和2号点的代号，最后一列是两点距离,数据链接: https://pan.baidu.com/s/1hF6Y0zf7c3mdmRdJzwTOjg 提取码: 7x45
接下来是代码部分，实现了其核心功能，最后输出了聚类中心的文件

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

#load data
print("load data file")
xx = pd.read_csv('D:\matlab program\example_distances.dat',sep=" ",header=None)
print(len(xx))

NL=max(xx[0])
ND=max(xx[1])
print(NL)
if NL>ND:
    ND=NL #确保 ND 取为第一二列最大值中的较大者，并将其作为数据点总数
N = len(xx)#距离个数  xx 第一个维度的长度，相当于文件的行数（即距离的总个数）

#初始化为dist数组为0
print(ND)
dist=np.zeros([ND,ND])
print(dist.shape[0])
#print(xx[2][2000])

#为dist数组赋值
# print(1999*1999)
for i in range(N):
    ii = xx[0][i]-1
    jj = xx[1][i]-1
    dist[ii][jj]=xx[2][i]
    dist[jj][ii]=xx[2][i]
# 测试赋值所需时间
# print(dist[1][30])

#确定dc
percent=2.0
print("average percentage of neighbours:",percent)

position=round(N*percent/100); #round 是一个四舍五入函数
sda=xx[2].sort_values().tolist()
#print(type(sda))
# print(len(sda))
# print(sda[0:10])
dc=sda[position]
print("截断距离是:",dc)

#计算局部密度 rho(用Gaussian高斯核）
from math import exp

#1.将每个数据点的rho值初始化为0
rho=np.zeros(ND)
# print(rho)

#高斯核计算局部密度
for i in range(0,ND-1):
  for j in range(i+1,ND):
     rho[i]=rho[i]+exp(-(dist[i,j]/dc)*(dist[i,j]/dc))
     rho[j]=rho[j]+exp(-(dist[i,j]/dc)*(dist[i,j]/dc))

 #"Cut off" kernel

#for i in range(1,ND-1):
#  for j in range(i+1,ND):
#    if dist[i,j]<dc:
#       rho[i]=rho[i]+1.;
#       rho[j]=rho[j]+1.;

# 先求矩阵列最大值，再求最大值，最后得到所有距离值中的最大值
maxd=dist.max()
print(maxd)

#将密度rho降序排列，ordho保持序
ordrho = np.argsort(-1*rho)
rho_sorted =rho[ordrho]

#初始化delta数组和nneigh数组
delta=[maxd]*ND
nneigh=[0]*ND
 #处理 rho 值最大的数据点
delta[ordrho[0]] = -1
nneigh[ordrho[0]] = 0

#delta 和 nneigh 数组赋值
for ii in range(1,ND):
    for jj in range(0,ii):
        if dist[ordrho[ii],ordrho[jj]]<delta[ordrho[ii]]:
            delta[ordrho[ii]]=dist[ordrho[ii],ordrho[jj]]
            nneigh[ordrho[ii]]=ordrho[jj]
# 记录 rho 值更大的数据点中与 ordrho(ii) 距离最近的点的编号 ordrho(jj)

# 生成 rho 值最大数据点的 delta 值
delta[ordrho[0]]=max(delta)

# 决策图
print('Generated file:DECISION GRAPH')
print('column 1:Density')
print('column 2:Delta')
#开画图文件写画图文件
with open('D:/matlab program/DECISION_GRAPH','w') as fid:
    for i in range(ND):
        fid.write('%6.2f %6.2f\n'%(rho[i],delta[i]))
#选择一个围住类中心的矩形
#print('Select a rectangle enclosing cluster centers')
#scrsz = get(0,'ScreenSize');
#figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]);
plt.xlabel('密度rho')
plt.ylabel('距离delta')
plt.scatter(rho,delta)
plt.show()


#输入选取聚类中心的范围（最小值）
rhomin=int(input("please input rho密度最小值"))
deltamin=float(input("please input delta距离最小值"))

#初始化 cluster 个数
NCLUST=0;

# cl 为归属标志数组，cl[i]=j 表示第 i 号数据点归属于第 j 个 cluster
#先统一将 cl 初始化为 -1
cl=[-1]*ND
#初始化 点映射到类别列表
icl=[-1]*10

#在矩形区域内统计数据点（即聚类中心）的个数
for i in range(ND):
    if  (rho[i]>rhomin) and (delta[i]>deltamin):
     NCLUST=NCLUST+1;
     cl[i]=NCLUST # 第 i 号数据点属于第 NCLUST 个 cluster
     #icl[NCLUST]=i # 逆映射,第 NCLUST 个 cluster 的中心为第 i 号数据点

print('NUMBER OF CLUSTERS: %i \n'%NCLUST)
print('Performing assignation')

#将其他数据点归类 (assignation)
for i in range(ND):
    if (cl[ordrho[i]]==-1):
        cl[ordrho[i]]=cl[nneigh[ordrho[i]]]

#由于是按照 rho 值从大到小的顺序遍历,循环结束后, cl 应该都变成正的值了

halo=[-1]*ND
rho_aver=0
#处理光晕点：
if (NCLUST>1):
    for i in range(ND):
        halo[i]=cl[i]
    bord_rho=np.zeros(NCLUST+1)
#初始化数组 bord_rho 为 0,每个 cluster 定义一个 bord_rho 值
#获取每一个 cluster 中平均密度的一个界 bord_rho
    for i in range(ND - 1):
        for j in range(i+1,ND):
# 距离足够小但不属于同一个cluster的i和j
            if (cl[i]!=cl[j]) and (dist[i,j]<= dc):
                rho_aver = (rho[i] + rho[j]) / 2
# 取i, j两点的平均局部密度
                if rho_aver > bord_rho[cl[i]]:
                    bord_rho[cl[i]] = rho_aver

                if rho_aver > bord_rho[cl[j]]:
                    bord_rho[cl[j]] = rho_aver



#halo 值为 0 表示为 outlier
for i in range(ND):
    if rho[i]<bord_rho[cl[i]]:
        halo[i]=0

#逐一处理每个 cluster
for i in range(NCLUST):
  nc=0; #用于累计当前 cluster 中数据点的个数
  nh=0; # 用于累计当前 cluster 中核心数据点的个数
  for j in range(ND):
    if cl[j]==i:
      nc=nc+1

    if halo[j]==i:
      nh=nh+1

 #print('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i \n'% (i,icl[i],nc,nh,nc-nh))

#写聚类后的打标文件
with open('D:/matlab program/clustered data','w') as fp:
    for i in range(ND):
        fp.write('%d,%6.2f,%6.2f,%d\n'%(i,rho[i],delta[i],cl[i]))

# colors = ['b','g','r','orange','c','m','y','k','w','p','gray']

# plt.xlabel('密度rho')
# plt.ylabel('距离delta')
# for i in range(ND):
#     for j in range(len(colors)):
#         if int(cl[i]) == j:
#             plt.plot(rho[i], delta[i], colors[j])
#             continue
# plt.show()



print("file:D:/matlab program/clustered data")

print('column 1:element id')
print('column 2:cluster assignation without halo control')
print('column 3:cluster assignation with halo control')

本文地址：https://blog.csdn.net/qq_43305009/article/details/107553581