clear all
x(1,1)=1.1;x(1,2)=2.1;
x(2,1)=1.5;x(2,2)=2.8;
x(3,1)=2.1;x(3,2)=2.5;
x(4,1)=0.9;x(4,2)=3.1;
x(5,1)=1.8;x(5,2)=3.4;
x(6,1)=4.2;x(6,2)=4.1;
x(7,1)=4.3;x(7,2)=5.5;
x(8,1)=5.1;x(8,2)=4.6;
x(9,1)=5.3;x(9,2)=3.1;
x(10,1)=5.6;x(10,2)=5.8;
r1=0;r2=0;m=0;d=0;
c=5;                                  %種群數
for i=1:10
    r1(i)=x(i,1);
    r2(i)=x(i,2);
end
figure(1)
plot(r1,r2,'o')
axis([0,10,0,10])
clus=0;clus_flag=0;
index=randperm(10);                    %選擇代表點
clus1=r1(index);
clus2=r2(index);
for i=1:c   
    clus_flag(i)=i;
end
clus_flag(c+1:10)=0;
m1=r1(index(1:c));
m2=r2(index(1:c));
min=0;

for i=1:c
    index_1(i)=1;
end

for i=c+1:10                         %初始分類           
    for j=1:c
        d(j)=((r1(index(i))-m1(j))^2+(r2(index(i))-m2(j))^2)^(1/2);
    end
    k=1;
    min=d(k);p=k;
    for k=2:j
        if(min>d(k))
            min=d(k);
            p=k;
        end
    end
    clus_flag(i)=p;                 %clus中第i個屬于第p類
    index_1(p)=index_1(p)+1;        %第p類的個數加1
    sum1=0;sum2=0;
    for m=1:10                      %修改均值
        if(clus_flag(m)==p)
             sum1=sum1+clus1(m);
             sum2=sum2+clus2(m);
        end
    end
    m1(p)=sum1/index_1(p);
    m2(p)=sum2/index_1(p);
end
clus_flag1=0;
clus_flag1=clus_flag;
Je1=0;                              %求代價函數
for n=1:10
    Je1=Je1+(clus1(n)-m1(clus_flag(n)))^2+(clus2(n)-m2(clus_flag(n)))^2;
end


g=1;Je_min=0;
Je2=0;
h=1;
Je=0;
Je(1)=Je1;
while(flag<50)
for l=1:10                         %重新分類
   if(index_1(clus_flag(l))~=0)
    for q=1:c
        if(clus_flag(l)~=q)
          
           if(index_1(clus_flag(l))~=1)
               
               Je2=Je1+(index_1(q)/(index_1(q)+1))*...
                   ((clus1(l)-m1(q))^2+(clus2(l)-m2(q))^2)-...
                   (index_1(clus_flag(l))/(index_1(clus_flag(l))-1))*...
                   ((clus1(l)-m1(clus_flag(l)))^2+(clus2(l)-m2(clus_flag(l)))^2);               
           
           if(Je2<Je1)               
               u=q;
               Je_min=Je2;
           else
               u=clus_flag(l);
           end 
           end
        end
    end
    Je1=Je_min;
    index_1(clus_flag(l))=index_1(clus_flag(l))-1;  %歸類
    clus_flag(l)=u;                 
    index_1(u)=index_1(u)+1;
    
    sum3=0;sum4=0;sum5=0;sum6=0;
    for m=1:10                      %修改均值
        if(clus_flag(m)==u)
             sum3=sum3+clus1(m);
             sum4=sum4+clus2(m);
        end
    end
    m1(u)=sum3/index_1(u);
    m2(u)=sum4/index_1(u);
    for n=1:10                      %修改均值
        if(clus_flag(n)==clus_flag(l))
             sum5=sum5+clus1(n);
             sum6=sum6+clus2(n);
        end
    end
    m1(clus_flag(l))=sum5/index_1(clus_flag(l));
    m2(clus_flag(l))=sum6/index_1(clus_flag(l));
    
    Je(h+1)=Je1;
    h=h+1;
   end  
end
flag=flag+1;
end
 clus=[clus1;clus2;clus_flag1;clus_flag]';