#ISODATA聚类算法
#ISODATA算法加入试探步骤，结合人机交互的结构

from numpy import*
from math import sqrt
import sys
#分类函数
def classification(x,z):
    (mx,nx)=x.shape
    
    #N存储各个聚类域包含的点号数
    N=[]
    x_z_dis=[]
    
    #考虑初始时一个聚类中心的情况
    if len(z) == 1:
        for i in range(0,mx):
            N.append(i)
    #考虑两个或两个以上聚类中心的情况
    else:
        for j in range(0,len(z)):
            N.append([])
        
        for i in range(0,mx):
            x_z_dis=[]
            for j in range(0,len(z)):
                x_z_dis.append((x[i][0]-z[j][0])**2+(x[i][1]-z[j][1])**2)
            
            #根据各样本点到聚类中心的距离情况分类
            x_z_min=x_z_dis[0]
            position=0
            for k in range(0,len(x_z_dis)):
                
                if x_z_dis[k] <= x_z_min:
                    position=k
                    x_z_min=x_z_dis[k] 
            
            N[position].append(i)
    
    return N

#根据样本分类情况获得聚类中心点     返回格式 z[[]]
def get_center(x,N):
    (mx,nx)=x.shape
    tem_z=[0,0]
    tem_1=0
    tem_2=0
    
    z=[]
    if len(N) == len(x):
        for i in range(0,mx):
            tem_z=tem_z+x[i]/mx
        z.append(tem_z)
    else:
        for i in range(0,len(N)):
            for j in range(0,len(N[i])):
                tem_1=tem_1+x[N[i][j]][0]/len(N[i])
                tem_2=tem_2+x[N[i][j]][1]/len(N[i])
                tem_z=[tem_1,tem_2]
            z.append(tem_z)
    return z 

#获得样本到相应聚类中心的距离 返回格式 dis[ , , ]
def get_distance(x,z,N):
    d=0
    dis=[]
    if len(N)==len(x):
        for j in range(0,len(N)):
            d=d+sqrt((z[0][0]-x[N[j]][0])**2+(z[0][1]-x[N[j]][1])**2)/len(N)
        dis.append(d)
    else:
        for i in range(0,len(N)):
            for j in range(0,len(N[i])):
                d=d+sqrt((z[i][0]-x[N[i][j]][0])**2+(z[i][1]-x[N[i][j]][1])**2)/len(N[i])
            dis.append(d)
    return dis
    
#计算各聚类域的标准差向量
def get_sigma(x,z,N):
    if len(N)==len(x):
        sig1=0
        sig2=0
        for j in range(0,len(N)):
            sig1=sig1+(z[0][0]-x[N[j]][0])**2/len(N)
            sig2=sig2+(z[0][1]-x[N[j]][1])**2/len(N)
        sig1=sqrt(sig1)
        sig2=sqrt(sig2)
        sigma=array([[sig1,sig2]])
    else:
        sigma=arange(2*len(N)).reshape(len(N),2)
        for i in range(0,len(N)):
            sig1=0
            sig2=0
            for j in range(0,len(N[i])):
                sig1=sig1+(z[i][0]-x[N[i][j]][0])**2/len(N[i])
                sig2=sig2+(z[i][1]-x[N[i][j]][1])**2/len(N[i])
            sig1=sqrt(sig1)
            sig2=sqrt(sig2)
            sigma[i]=array([sig1,sig2])
    return sigma
    
#获得列表最大值的位置  a=[[],[],[],[]] 返回最大值所处位置
def get_Max_position(a):
    position=0
    (pos_1,max_1)=get_Max_value(a[0])
    max_s=max_1
    for i in range(1,len(a)):
        (tem_p,tem_m)=get_Max_value(a[i])
        if (tem_m > max_s):
            max_s=tem_m
            position=i
    return position
    
#获得单层列表中的较大值的位置与数值 a=[] 
def get_Max_value(a):
    max_s=a[0]
    position=0
    for i in range(1,len(a)):
        if a[i]>max_s:
            max_s=a[i]
            position=i
    return position,max_s
    
#获得每个聚类中心之间的距离以及记录
def get_z_dis(z):
    z_dis=[]
    the_pos=[]
    for i in range(0,len(z)-1):
        for j in range(i+1,len(z)):
            z_dis.append(sqrt((z[i][0]-z[j][0])**2+(z[i][1]-z[j][1])**2))
            the_pos.append([i,j])
    return z_dis,the_pos
    
#获得所有样本到其相应聚类中心的距离的平均值
def average_dis(x,dis,N):
    ave_dis=0
    if len(N)==len(x):
        ave_dis=dis[0]
    else:
        for i in range(0,len(dis)):
            ave_dis=ave_dis+dis[i]*len(N[i])
        ave_dis=ave_dis/len(x)
    return ave_dis
            
#更新六大参数
def update_parameter():
    K=int(input("K is updated to: "))
    the_N=int(input("the_N is updated to: "))
    the_S=int(input("the_S is updated to: "))
    the_C=int(input("the_C is updated to: "))
    L=int(input("L is updated to: "))
    I=int(input("I is updated to: "))
    return K,the_N,the_S,the_C,L,I

#对类别进行分裂处理
def split_z(z,N,n,sigma,dis,the_S,C,K):
    ave_d=average_dis(x,dis,N)
    
    if sigma.max() > the_S and (C <= K/2 or ((dis[get_Max_position(sigma)]>ave_d) and len(N[get_Max_position(sigma)])>2*(the_N+1))):
        z_update=[]
        (pos,val)=get_Max_value(sigma[get_Max_position(sigma)])
        r=0.5*val
        
        if pos ==0:
            z.append([z[get_Max_position(sigma)][0]+r,z[get_Max_position(sigma)][1]])
            z.append([z[get_Max_position(sigma)][0]-r,z[get_Max_position(sigma)][1]])
        else:
            z.append([z[get_Max_position(sigma)][0],z[get_Max_position(sigma)][1]+r])
            z.append([z[get_Max_position(sigma)][0],z[get_Max_position(sigma)][1]-r])
        #更新聚类中心    
        
        del z[get_Max_position(sigma)]
        
        C=C+1
        #迭代次数加一
        n=n+1
        
        ISODATA(x,C,z,n,K,the_N,the_S,the_C,L,I)

#合并聚类中心
def gather_z(z,N,the_C,C):
    #获得每个聚类中心之间的距离以及记录
    (z_dis,the_pos)=get_z_dis(z)
    #z_pos用于存储距离小于阈值的z_dis位置
    z_pos=[]
    for i in range(0,len(z_dis)):
        if z_dis[i]<the_C:
            z_pos.append(i)    
            
    #防止某类别被多次归并 定义下变量
    used_num=[]
    if len(z_pos)==0:
        return z,C
        
    for i in range(0,len(z_pos)):
            if the_pos[z_pos[i]][0] in used_num or the_pos[z_pos[i]][1]in used_num:
                continue
            else:
                used_num.append(the_pos[z_pos[i]][0])
                used_num.append(the_pos[z_pos[i]][1])
                
                z_new=get_new_z(z[the_pos[z_pos[i]][0]],z[the_pos[z_pos[i]][1]],len(N[the_pos[z_pos[i]][0]]),len(N[the_pos[z_pos[i]][1]]))
                z.append(z_new)
                if the_pos[z_pos[i]][0]<the_pos[z_pos[i]][1]:
                    del z[the_pos[z_pos[i]][1]]
                    del z[the_pos[z_pos[i]][0]]
                else:
                    del z[the_pos[z_pos[i]][0]]
                    del z[the_pos[z_pos[i]][1]]
                
                C=C-1
    return z,C
        
#合并聚类中心中间操作  num_1为z_1类别样点数 num_2为z_2类别样点数
def get_new_z(z_1,z_2,num_1,num_2):
    z_new=[]
    z_new_1=(num_1*z_1[0]+num_2*z_2[0])/(num_1+num_2)
    z_new_2=(num_1*z_1[1]+num_2*z_2[1])/(num_1+num_2)
    z_new=[z_new_1,z_new_2]
    
    return z_new

#选择是否修改参数
def whether_modify(x,C,z,n,K,the_N,the_S,the_C,L,I):
    b=input("Need you change the parameter?(Y/N)\n")
    if b=="Y":
        (K,the_N,the_S,the_C,L,I)=update_parameter()
        n=n+1
        ISODATA(x,C,z,n,K,the_N,the_S,the_C,L,I)
    elif b=="N":
        n=n+1
        ISODATA(x,C,z,n,K,the_N,the_S,the_C,L,I)
    else:
        print("Enter a wrong word!")
        whether_modify(x,C,z,n,K,the_N,the_S,the_C,L,I)

#展示成果
def show_result(N,z):
    print("Clustering center: "+str(z))
    print("Classification: ")
    if len(N) ==len(x):
        print("Center include "+str(N))
    else:
        for i in range(0,len(N)):
            print("Center "+str(i)+" include "+str(N[i]))
    
#ISODATA聚类函数
def ISODATA(x,C,z,n,K,the_N,the_S,the_C,L,I):
    #N为各聚类域分类情况
    
    N=classification(x,z)
    #取消样本数低于阈值的样本区间
    
    if len(N) != len(x):
        for i in range(0,len(N)):
            if len(N[i]) <= the_N:
                del N[i]
                del z[i]
                C=C-1
                break
        
        
    #对各聚类域再确定聚类中心
    z=get_center(x,N)
    dis=get_distance(x,z,N)
    ave_d=average_dis(x,dis,N)
    
    if n >=I:
        the_C=0
        (z,C)=gather_z(z,N,the_C,C)
        N=classification(x,z)
        show_result(N,z)
        sys.exit(0)
        
    elif n%2==0 or C>=2*K:
        (z,C)=gather_z(z,N,the_C,C)
        
    else:
        
        sigma=get_sigma(x,z,N)
        split_z(z,N,n,sigma,dis,the_S,C,K)
        (z,C)=gather_z(z,N,the_C,C)
    
    if n >= I:
        N=classification(x,z)
        show_result(N,z)
        sys.exit(0)
    else:
        whether_modify(x,C,z,n,K,the_N,the_S,the_C,L,I)
            
    
#初始化算法需要确定并在计算中可以调整的参数：
#所要求聚类中心数
K=2
#每个类别至少应具有的样本数目
the_N=1
#每个类别标准差阈值
the_S=1
#聚类中心之间距离的阈值——归并系数
the_C=2
#一次迭代能归并的最多对数
L=0
#允许迭代的最多次数
I=4

x=array([[0,0],[3,8],[2,2],[1,1],[5,3],
         [4,8],[6,3],[5,4],[6,4],[7,5]])

#取初始聚类中心数C 聚类中心为z 函数迭代次数为n
C=1
z=[x[0]]
n=1

ISODATA(x,C,z,n,K,the_N,the_S,the_C,L,I)