%KHC
global cl;
cl=5;
SPDW=PDW;
format long
b=0.01;
sum1=0;
sum2=0;
sum3=0;
mean1=0;
mean2=0;
mean3=0;
var1=0;
var2=0;
var3=0;
c=1;
ii=0;
e=0;
D=zeros(500,cl);
D1=zeros(500,cl);
K=zeros(500,500);
kxw=zeros(500,cl);
kxw1=zeros(500,cl);
kww=zeros(cl,1);
kww1=zeros(cl,1);
sumk=zeros(500,1);
sumd=1;
sumd1=1;
sumx=0;
dd=zeros(500,cl);
ee11=zeros(cl,1);
ee1=0;
ee2=0;
ee3=0;
Emax=2;
Cluster=zeros(500,6*cl);
%for iii=1:(N(1)+N(2)+N(3)+N(4)+N(5)+N(6))
   %if PDW(iii,2)>0
    %   ii=ii+1;
 %   end
%end
for j=1:500
    sum1=sum1+PDW(j,1);
    sum2=sum2+PDW(j,2);
    sum3=sum3+PDW(j,3);
end
mean1=sum1/500;
mean2=sum2/500;
mean3=sum3/500;
for k=1:500
   var1=var1+(PDW(k,1)-mean1)^2;
   var2=var2+(PDW(k,2)-mean2)^2;
   var3=var3+(PDW(k,3)-mean3)^2;
end
var1=var1/(500-1);
var2=var2/(500-1);
var3=var3/(500-1);
for i=1:500
   SPDW(i,1)=(PDW(i,1)-mean1)/sqrt(var1);
   SPDW(i,2)=(PDW(i,2)-mean2)/sqrt(var2);
   SPDW(i,3)=(PDW(i,3)-mean3)/sqrt(var3);
end
for i=1:500
    for j=1:500
        K(i,j)=exp((-((SPDW(i,3)-SPDW(j,3))^2))/(2*b*b));%+(SPDW(i,2)-SPDW(j,2))^2+(SPDW(i,3)-SPDW(j,3))^2
    end
end
for j=1:500
    e=cl*rand(1);
    D(j,ceil(e))=1;
end
for k=1:cl
    for j=1:500
        sumk=sumk+D(j,k)*K(:,j);
        sumd=sumd+D(j,k);
        for i=1:500
            sumx=sumx+D(j,k)*D(i,k)*K(j,i);
        end
    end
    kxw(:,k)=sumk/sumd;
    kww(k,1)=sumx/((sumd)^2);
    sumk=zeros(500,1);
    sumd=1;
    sumx=0;
end
e=15;
%while(e>Emax)
for c=1:40
 kxw1=kxw;
 kww1=kww;
 D1=D;
 for j=1:500
        for k=1:cl
        dd(j,k)=sqrt(K(j,j)-2*kxw(j,k)+kww(k,1));
        end
 end
 for j=1:500
        for k=1:cl
            if dd(j,k)==min(dd(j,:));
               D(j,k)=1;
            else D(j,k)=0;
            end
        end
       % if dd(j,1)==dd(j,2);
           %     D(j,2)=0;
           % elseif dd(j,1)==dd(j,3);
            %   D(j,3)=0;
           % elseif dd(j,2)==dd(j,3);
           %     D(j,3)=0;
        %end
 end
 for k=1:cl
    for i=1:500
        for j=1:500
            sumx=sumx+D(i,k)*D(j,k)*K(i,j);   
        end
        sumk=sumk+D(i,k)*K(:,i);
        sumd=sumd+D(i,k); 
    end
   kxw(:,k)=sumk/(sumd);
   kww(k,1)=sumx/((sumd)^2);
   sumk=zeros(500,1);
   sumd=1;
   sumx=0;
 end
 for k=1:cl
    ee3=ee3+kww(k,1);
    ee2=ee2+kww1(k,1);
    for j=1:500
        for i=1:500
          ee11(k,1)=ee11(k,1)+D(j,k)*D1(i,k)*K(j,k);
        end
        sumd1=sumd1+D1(j,k);
        sumd=sumd+D(j,k);
    end
    ee1=ee1+ee11(k,1)/(sumd*sumd1);
    ee11=zeros(cl,1);
    sumd1=1;
    sumd=1;  
 end
 e=ee3-2*ee1+ee2;
 e=abs(e);
 ee1=0;
 ee2=0;
 ee3=0;
 %c=c+1;    
end

%***完成聚類指派***
des=zeros(cl,1);
Cluster=zeros(500,6*cl); 
for i=1:cl
   for j=1:500
        if D(j,i)==1
           des(i,1)=des(i,1)+1;
           Cluster(des(i,1),((i-1)*6+1):i*6)=PDW(j,1:6);
       end
    end
end