cluster_ser.m

動態聚類k-means演算 將輸入在程式中的數據資料 給予適當的分群

function [dbelong, neach, clusters, niter, ST, DB] = ...
    cluster_ser(indata, nc, np, incnp, eta, deta, maxEpoch, info, psec);
%
% The size of the dataset
[ndata,nvar] = size(indata);
class=ones(1,ndata);
col=class;
leaf=0;
if info >= 1
  disp(sprintf('Number of datapoints: %d', ndata));
  disp(sprintf('Number of variables:  %d', nvar));
end

% Normalize data to unit standard deviation and zero mean
% using prestd.
[indata, datamean, datastd] = kmeans_norm(indata); 

%%%%%%%%%% Parameters for k-means %%%%%%%%%%

    npoints = ndata;

    clust(1,1)=0;
    clust(1,2)=0;

  fig = kmeans_plotdata(indata,col);
CP=zeros(1,200);
Ncl=2^20-1;  
clust=zeros(1000,2);
dotter=1000*ones(1000,2);
mor=dotter;
N=777*ones(1,500);
N(1)=200;
ant=1;
cl=0;
go2=1;
while leaf<200 & cl<900
    cl=cl+1;
       
    % The starting clusters
    u=find(class==cl);
    
    if length(u)>=3
        rand('state',sum(100*clock));
        clusters(1,:) = 0.0001 * (rand(1,2) - 0.5)+clust(cl,:);   % Slumpar startpunkter.
        clusters(2,:) = 0.0001 * (rand(1,2) - 0.5)+clust(cl,:);   % Slumpar startpunkter.
        bclusters=clusters;
        niter = 0;
        go = 1;
        qerr_old = 0.0;
        cold = zeros(nc,nvar);
      
    
        while go == 1

            niter = niter + 1;

            data = indata;
  
            % Compute the distance between each datapoint to the clusters
            dists = kmeans_dist(indata, clusters);
    
            % Assign datapoints to clusters and compute updates
            [minval, minidx] = min(dists');
            for i=1:nc

                idxtmp = find(minidx == i & class==cl); 
                if length(idxtmp > 0)
                    if cl<8
                        if i==1
                            plot (indata(idxtmp,1),indata(idxtmp,2),'.g')
                        else
                            plot (indata(idxtmp,1),indata(idxtmp,2),'.m')
                        end
                    end
                    
                    cpoints = indata(idxtmp,:);
	
	                % Compute the mean of the indatapoints belonging to cluster i
	                if length(idxtmp)==1
                        centroid=cpoints;
                    else
                        centroid = mean(cpoints);
                    end
                    update(i,:) = centroid - clusters(i,:);
                    
                
                else
	                update(i,:) = zeros(1,2);
                    
                    
                    
                    
                end
            end
            
            % Update the cluster centers
            clusters = clusters + eta * update;    
            update_old=update;
            
            
            % Compute changes
            dcluster = sum( sqrt( sum( ( (clusters-cold).*(clusters-cold) )') ) ) / ...
            nc;
            cold = clusters;
            %dcluster = sum( sqrt( sum( (update.*update)') ) ) / nc;

            % Compute some clusters statistics
            [dbelong, neach, avedist, quanterr] = ...
            cluster_stat(indata, clusters, nc,cl,class);
  
            % Some plotting if 2 dimenstions
            if nvar == 2 & info > 0 & cl<8  
                plot(clusters(:,1), clusters(:,2),'+r');
                %cluster_plotclusters(fig, clusters,cl,clust);

            end
  
            % Stopping criteria
            if abs(quanterr - qerr_old) < 0.00001 | niter >= maxEpoch 
                go = 0;
            end
            qerr_old = quanterr;
  
            % Update the learning rate
            %eta = eta * deta;
  
            % This for visualization purposes
            %pause(psec);
            if cl<8
                pause(.5)
            end
        end

        
        % Final plotting if in 2 dimenstions
        
        for i=1:2
            idxtmp = find(minidx == i & class==cl);
            ant=ant+1;
            class(idxtmp)=ant;
            N(ant)=length(idxtmp);  
            dotter(cl,i)=ant;    
            mor(ant)=cl;
        end

            
        
        clust(ant-1:ant,:)=clusters;

        if nvar == 2 & info > 0
            
            if cl<8
                close(fig);
                if cl==7 
                    col=ones(1,ndata);
                end
                fig=kmeans_plotdata(indata,col);
                drawnow    
            end
            cluster_plotclusters(fig, clusters,cl,clust,dotter);
        end
        ratio=num2str(floor(100*cl/399));
        xlabel([ratio,'%'])
    elseif length(u)==2
        
        ant=ant+1;
        clust(ant,:)=indata(u(1),:);
        dotter(cl,1)=ant;
        mor(ant)=cl;
        class(u(1))=ant;
        N(ant)=1;
        ant=ant+1;
        clust(ant,:)=indata(u(2),:);
        dotter(cl,2)=ant;
        mor(ant)=cl;
        class(u(2))=ant;
        N(ant)=1;
        cluster_plotclusters(fig, clusters,cl,clust,dotter);    
        
    elseif length(u)==1
        CP(u)=2;
        leaf =leaf+1;
        
    end

end



                
figure
hold on

for k=1:cl
    
    plot (clust(k,1),clust(k,2),'+k')
    d=dotter(k,1);
    if clust(d,1)
        plot ([clust(k,1) clust(d,1)],[clust(k,2) clust(d,2)],'k')
        plot(clust(d,1), clust(d,2),'+k');
    end
    
    d=dotter(k,2);
    if clust(d,1)
        plot ([clust(k,1) clust(d,1)],[clust(k,2) clust(d,2)],'k')
        plot(clust(d,1), clust(d,2),'+k');
        
    end
    
    
end
plot(indata(1:50,1), indata(1:50,2),'.g');
plot(indata(51:100,1), indata(51:100,2),'.b');
plot(indata(101:150,1), indata(101:150,2),'.r');
plot(indata(151:200,1), indata(151:200,2),'.m');


drawnow;

g=ones(1,1000);
y=zeros(1,1000);
x=y;
for i=1:cl
    d=dotter(i,:);
    if d(1)<1000
        
        g(d(1))=g(i)+1;
        x(d(1))=x(i)+.5^g(i);
        g(d(2))=g(i)+1;
        x(d(2))=x(i)-.5^g(i);
    end
end
u=find(N(1:cl)==1);

[tmp, q]=sort(x(u));
z=u(q);

y(z)=1:length(z);

[tmp, j]=sort(N(1:cl));

ccol=zeros(1,cl);
for r=1:50
    i=class(r);
    ccol(i)=1;
end
for r=51:100
    i=class(r);
    ccol(i)=2;
end
for r=101:150
    i=class(r);
    ccol(i)=3;
end
for r=151:200
    i=class(r);
    ccol(i)=4;
end



figure
hold on
for i=1:length(j)
    s=j(i);
    d=dotter(s,:);
    if d(1)<1000
        y(s)=(y(d(1))+y(d(2)))/2;
        plot ([g(d(1)), g(s), g(s), g(d(2))], [y(d(1)), y(d(1)), y(d(2)), y(d(2))],'k')
        if ccol(d(1)) == ccol(d(2))
            ccol(s)=ccol(d(1));
        end
    end
end
plot ([g(1)-1,g(1)],[y(1),y(1)],'k')

for i=1:cl
    if ccol(i)==1
        plot ([g(i)-1,g(i)],[y(i),y(i)],'g')
    end
    if ccol(i)==2
        plot ([g(i)-1,g(i)],[y(i),y(i)],'b')
    end
    if ccol(i)==3
        plot ([g(i)-1,g(i)],[y(i),y(i)],'r')
    end
    if ccol(i)==4
        plot ([g(i)-1,g(i)],[y(i),y(i)],'m')
    end
end
axis ([0 max(g) -10 210])