AP

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



'''
第一步：生成测试数据
    1.生成实际中心为centers的测试样本300个，
    2.Xn是包含150个(x,y)点的二维数组
    3.labels_true为其对应的真是类别标签
'''


def init_sample():
    ## 生成的测试数据的中心点
    centers = [[1, 1], [-2, -2]]
    ##生成数据
    Xn, labels_true = make_blobs(n_samples=6, centers=centers, cluster_std=0.3,
                                 random_state=0)
    # 3数据的长度，即：数据点的个数
    dataLen = len(Xn)

    return Xn, dataLen


'''
第二步：计算相似度矩阵
'''


def cal_simi(Xn):
    ##这个数据集的相似度矩阵，最终是二维数组
    simi = []
    for m in Xn:
        ##每个数字与所有数字的相似度列表，即矩阵中的一行
        temp = []
        for n in Xn:
            ##采用负的欧式距离计算相似度
            s = -np.sqrt((m[0] - n[0]) ** 2 + (m[1] - n[1]) ** 2)
            temp.append(s)
        simi.append(temp)

    ##设置参考度，即对角线的值，一般为最小值或者中值
    # p = np.min(simi)   ##11个中心
    # p = np.max(simi)  ##14个中心
    p = np.median(simi)  ##5个中心
    for i in range(dataLen):
        simi[i][i] = p
    return simi


'''
第三步：计算吸引度矩阵，即R
       公式1：r(n+1) =s(n)-(s(n)+a(n))-->简化写法，具体参见上图公式
       公式2：r(n+1)=(1-λ)*r(n+1)+λ*r(n)
'''


##初始化R矩阵、A矩阵
def init_R(dataLen):
    R = [[0] * dataLen for j in range(dataLen)]
    return R


def init_A(dataLen):
    A = [[0] * dataLen for j in range(dataLen)]
    return A


##迭代更新R矩阵
def iter_update_R(dataLen, R, A, simi):
    old_r = 0  ##更新前的某个r值
    lam = 0.1  ##阻尼系数,用于算法收敛
    ##此循环更新R矩阵
    for i in range(dataLen):
        for k in range(dataLen):
            old_r = R[i][k]
            if i != k:
                max1 = A[i][0] + R[i][0]  ##注意初始值的设置
                for j in range(dataLen):
                    if j != k:
                        if A[i][j] + R[i][j] > max1:
                            max1 = A[i][j] + R[i][j]
                ##更新后的R[i][k]值
                R[i][k] = simi[i][k] - max1
                ##带入阻尼系数重新更新
                R[i][k] = (1 - lam) * R[i][k] + lam * old_r
            else:
                max2 = simi[i][0]  ##注意初始值的设置
                for j in range(dataLen):
                    if j != k:
                        if simi[i][j] > max2:
                            max2 = simi[i][j]
                ##更新后的R[i][k]值
                R[i][k] = simi[i][k] - max2
                ##带入阻尼系数重新更新
                R[i][k] = (1 - lam) * R[i][k] + lam * old_r
    print("max_r:" + str(np.max(R)))
    # print(np.min(R))
    return R


'''
    第四步：计算归属度矩阵，即A
'''


##迭代更新A矩阵
def iter_update_A(dataLen, R, A):
    old_a = 0  ##更新前的某个a值
    lam = 0.1  ##阻尼系数,用于算法收敛
    ##此循环更新A矩阵
    for i in range(dataLen):
        for k in range(dataLen):
            old_a = A[i][k]
            if i == k:
                max3 = R[0][k]  ##注意初始值的设置
                for j in range(dataLen):
                    if j != k:
                        if R[j][k] > 0:
                            max3 += R[j][k]
                        else:
                            max3 += 0
                A[i][k] = max3
                ##带入阻尼系数更新A值
                A[i][k] = (1 - lam) * A[i][k] + lam * old_a
            else:
                max4 = R[0][k]  ##注意初始值的设置
                for j in range(dataLen):
                    ##上图公式中的i!=k 的求和部分
                    if j != k and j != i:
                        if R[j][k] > 0:
                            max4 += R[j][k]
                        else:
                            max4 += 0

                ##上图公式中的min部分
                if R[k][k] + max4 > 0:
                    A[i][k] = 0
                else:
                    A[i][k] = R[k][k] + max4

                ##带入阻尼系数更新A值
                A[i][k] = (1 - lam) * A[i][k] + lam * old_a
    print("max_a:" + str(np.max(A)))
    # print(np.min(A))
    return A


'''
   第5步：计算聚类中心
'''


##计算聚类中心
def cal_cls_center(dataLen, simi, R, A):
    ##进行聚类，不断迭代直到预设的迭代次数或者判断comp_cnt次后聚类中心不再变化
    max_iter = 200  ##最大迭代次数
    curr_iter = 0  ##当前迭代次数
    max_comp = 200  ##最大比较次数
    curr_comp = 0  ##当前比较次数
    class_cen = []  ##聚类中心列表，存储的是数据点在Xn中的索引
    while True:
        R = iter_update_R(dataLen, R, A, simi)
        A = iter_update_A(dataLen, R, A)

        ##开始计算聚类中心
        for k in range(dataLen):
            if R[k][k] + A[k][k] > 0:
                if k not in class_cen:
                    class_cen.append(k)
                else:
                    curr_comp += 1
        curr_iter += 1
        print(curr_iter)
        if curr_iter >= max_iter or curr_comp > max_comp:
            break
    return class_cen


if __name__ == '__main__':
    ##初始化数据
    Xn, dataLen = init_sample()
    ##初始化R、A矩阵
    R = init_R(dataLen)
    A = init_A(dataLen)
    ##计算相似度
    simi = cal_simi(Xn)
    ##输出聚类中心
    class_cen = cal_cls_center(dataLen, simi, R, A)
    # for i in class_cen:
    #    print(str(i)+":"+str(Xn[i]))
    # print(class_cen)

    ##根据聚类中心划分数据
    c_list = []
    for m in Xn:
        temp = []
        for j in class_cen:
            n = Xn[j]
            d = -np.sqrt((m[0] - n[0]) ** 2 + (m[1] - n[1]) ** 2)
            temp.append(d)
        ##按照是第几个数字作为聚类中心进行分类标识
        c = class_cen[temp.index(np.max(temp))]
        c_list.append(c)
    ##画图
    colors = ['red', 'blue', 'black', 'green', 'yellow']
    plt.figure(figsize=(8, 6))
    plt.xlim([-3, 3])
    plt.ylim([-3, 3])
    for i in range(dataLen):
        d1 = Xn[i]
        d2 = Xn[c_list[i]]
        c = class_cen.index(c_list[i])
        plt.plot([d2[0], d1[0]], [d2[1], d1[1]], color=colors[c], linewidth=1)

    plt.show()