Affinity Propagation (AP) 聚類是2007年在Science雜誌上提出的一種新的聚類演算法。它根據N個數據點之間的相似度進行聚類,這些相似度可以是對稱的,即兩個資料點互相之間的相似度一樣(如歐氏距離);也可以是不對稱的,即兩個資料點互相之間的相似度不等。這些相似度組成N×N的相似度矩陣S(其中N為有N個數據點)。

AP演算法不需要事先指定聚類數目,相反它將所有的資料點都作為潛在的聚類中心,稱之為 exemplar。以S矩陣的對角線上的數值s (k, k)作為k點能否成為聚類中心的評判標準,這意味著該值越大,這個點成為聚類中心的可能性也就越大,這個值又稱作參考度p ( preference) 。聚類的數量受到參考度p的影響,如果認為每個資料點都有可能作為聚類中心,那麼p就應取相同的值。如果取輸入的相似度的均值作為p的值,得到聚類數量是中等的。如果取最小值,得到類數較少的聚類。

AP演算法中傳遞兩種型別的訊息,(responsiility)和(availability) 。r(i,k)表示從點i傳送到候選聚類中心k的數值訊息,反映k點是否適合作為i點的聚類中心。a(i,k)則從候選聚類中心k傳送到i的數值訊息,反映i點是否選擇k作為其聚類中心。r (i, k)與a (i, k)越強,則k點作為聚類中心的可能性就越大,並且i點隸屬於以k點為聚類中心的聚類的可能性也越大。AP演算法通過迭代過程不斷更新每一個點的吸引度和歸屬度值,直到產生m個高質量的exemplar,同時將其餘的資料點分配到相應的聚類中。

在這裡介紹幾個文中常出現的名詞：

exemplar：指的是聚類中心。

similarity：資料點i和點j的相似度記為S(i，j)。是指點j作為點i的聚類中心的相似度。

preference：資料點i的參考度稱為P(i)或S(i,i)。是指點i作為聚類中心的參考度。一般取S相似度值的中值。

Responsibility：R(i,k)用來描述點k適合作為資料點i的聚類中心的程度。

Availability：A(i,k)用來描述點i選擇點k作為其聚類中心的適合程度。

Damping factor：阻尼係數，主要是起收斂作用的。

AP聚類演算法是將每個資料看成圖中的一個節點，迭代的過程即是在圖中通過傳播信 息來找到聚類集合。本文計算兩個資料點的相似度採用距離的負數，也就是說距離越近，相似度越大。相似矩陣S中i到j的相似度就是剛剛所說的距離的負數。但是主對角線上的那些數表示的是某個點和自身的相似度，但是這裡我們不能直接用0來表示。根據演算法要求，主對角線上的值s(k,k)一般稱為偏向引數，一般 情況下對所有k，s(k,k)都相等，取非主對角線上的所有數的中位數。這個值很重要，他的大小與最後得到的類的數目有關，一般而言這個數越大，得到的類的數目就越多。

這裡為什麼要設定一個偏向引數而不直接用0來算呢，估計是因為AP聚類演算法是要 用圖論的一些東西來理解的，它把所有的點都看成一個圖中的節點，通過節點之間的資訊傳遞來達到聚類的效果。具體比較複雜，形象一點說就是我告訴你我和這些人是死黨，如果你認為你也是我死黨的話，那你就加入我們這一堆人裡面來吧！

有一些詳細的原理上的東西就不說了，直接說計算過程吧。。聚類就是個不斷迭代的過程，迭代的過程主要更新兩個矩陣，代表(Responsibility)矩陣R = [r(i,k)]_N×N和適選(Availabilities)矩陣A=[a(i,k)]_N×N。這兩個矩陣才初始化為0，N是所有樣本的數目。r(i,k)表示第k個樣本適合作為第i個樣本的類代表點的代表程度，a(i,k)表示第i個樣本選擇第k個樣本作為類代表樣本的適合程度。迭代更新公式如下：

每次更新後就可以確定當前樣本i的代表樣本(exemplar)點k，k就是使{a(i,k)+r(i,k)}取得最大值的那個k，如果i=k的話，那麼說明樣本i就是自己這個cluster的類代表點，如果不是，那麼說明i屬於k所屬的那個cluster。

當然，迭代停止的條件就是所有的樣本的所屬都不在變化為止，或者迭代了n次都還沒有變化（n的值可以自己取）。

說起來還有一種判斷點屬於屬於哪一類的方法，就是找出所有決策矩陣主對角線元素{a(k,k)+r(k,k)}大於0的所有點，這些點全部都是類代表點，之後在決定其餘的點屬於這裡面的一類。這兩種方法的結果我沒比較過誒，不知是不是一樣的。

另外還有一點就是AP聚類演算法迭代過程很容易產生震盪，所以一般每次迭代都加上一個阻尼係數λ：

r_new(i,k) = λ*r_old(i,k) + (1-λ)*r(i,k)

a_new(i,k) = λ*a_old(i,k) + (1-λ)*a(i,k)

# coding:utf-8

from sklearn.datasets.samples_generator import make_blobs
import matplotlib.pyplot as plt
import random
import numpy as np

##############################################################################
# Generate sample data
centers = [[1, 1], [-1, -1], [1, -1]]
X, labels_true = make_blobs(n_samples=1000, centers=centers, cluster_std=0.4,
                            random_state=0)

##############################################################################


def euclideanDistance(X, Y):
	"""計算每個點與其他所有點之間的歐幾里德距離"""
	X = np.array(X)
	Y = np.array(Y)
	# print X
	return np.sqrt(np.sum((X - Y) ** 2))



def computeSimilarity(datalist):

	num = len(datalist)

	Similarity = []
	for pointX in datalist:
		dists = []
		for pointY in datalist:
			dist = euclideanDistance(pointX, pointY)
			if dist == 0:
				dist = 1.5
			dists.append(dist * -1)
		Similarity.append(dists)

	return Similarity


def affinityPropagation(Similarity, lamda):

	#初始化 吸引矩陣 和 歸屬 矩陣
	Responsibility = np.zeros_like(Similarity, dtype=np.int)
	Availability = np.zeros_like(Similarity, dtype=np.int)

	num = len(Responsibility)

	count = 0
	while count < 10:
		count += 1
		# update 吸引矩陣

		for Index in range(num):
			# print len(Similarity[Index])
			kSum = [s + a  for s, a in zip(Similarity[Index], Availability[Index])]
			# print kSum
			for Kendex in range(num):
				kfit = delete(kSum, Kendex)
				# print fit
				ResponsibilityNew = Similarity[Index][Kendex] - max(kfit)
				Responsibility[Index][Kendex] = lamda * Responsibility[Index][Kendex] + (1 - lamda) * ResponsibilityNew

		# print "Responsibility", Responsibility


		# update 歸屬矩陣

		ResponsibilityT = Responsibility.T

		# print ResponsibilityT, Responsibility

		for Index in range(num):

			iSum = [r for r in ResponsibilityT[Index]]

			for Kendex in range(num):

				# print Kendex
				# print "ddddddddddddddddddddddddddd", ResponsibilityT[Kendex]
				#
				ifit = delete(iSum, Kendex)
				ifit = filter(isNonNegative, ifit)   #上面 iSum  已經全部大於0  會導致  delete 下標錯誤

				#   k == K 對角線的情況
				if Kendex == Index:
					AvailabilityNew  = sum(ifit)
				else:
					result = Responsibility[Kendex][Kendex] + sum(ifit)
					AvailabilityNew = result if result > 0 else 0
				Availability[Kendex][Index] = lamda * Availability[Kendex][Index] + (1 - lamda) * AvailabilityNew
		print "###############################################"
		print Responsibility
		print Availability
		print "###############################################"
	return Responsibility + Availability

def computeCluster(fitable, data):
	clusters = {}
	num = len(fitable)
	for node in range(num):
		fit = list(fitable[node])
		key = fit.index(max(fit))
		if not clusters.has_key(key):
			clusters[key] = []
		point = tuple(data[node])
		clusters[key].append(point)

	return clusters
##############################################################################

"""切片刪除 返回新陣列"""
def delete(lt, index):
    lt = lt[:index] + lt[index+1:]
    return lt

def isNonNegative(x):
	return x >= 0


##############################################################################

Similarity = computeSimilarity(X)

Similarity = np.array(Similarity)

print "Similarity", Similarity

fitable = affinityPropagation(Similarity, 0.34)

print fitable

clusters = computeCluster(fitable, X)

# print clusters

##############################################################################
clusters = clusters.values()

print len(clusters)

##############################################################################
def plotClusters(clusters, title):
	""" 畫圖 """
	plt.figure(figsize=(8, 5), dpi=80)
	axes = plt.subplot(111)
	col=[]
	r = lambda: random.randint(0,255)
	for index in range(len(clusters)):
		col.append(('#%02X%02X%02X' % (r(),r(),r())))
	color = 0
	for cluster in clusters:
		cluster = np.array(cluster).T
		axes.scatter(cluster[0],cluster[1], s=20, c = col[color])
		color += 1
	plt.title(title)
	# plt.show()
##############################################################################
plotClusters(clusters, "clusters by affinity propagation")
plt.show()

##############################################################################