clmsead.m

最新的模式识别分类工具箱,希望对朋友们有用!

% ---------------------------------------------------------------- 
% FUNCTION   clmsead.m    CLustering algorithm, MSE metric, adaptive
% ---------------------------------------------------------------- 
% 		builds a block training set. 
% Usage: 
%   cluster_centers = clmsead(orig_cluster_centers, training set,
%                       precision, threshsplit, threshmerge);
% 
% orig_cluster: matrix, each ROW is a cluster center.
% training set: matrix, each ROW is a sample.
% Precision = stopping criterion: If the difference between two following
%		mse's is less than epsilon, or if the decay in 
%		the dispersion is less than the threshold, then
%		inquire if better results come from splitting or merging,
%		then, if no splitting or merging occurs, terminate.
% threshsplit:  if the max average squared distance from the centroid
%		is greater than threshsplit* total dispersion, then split
%		CAVEAT: the higher the dimension, the higher Threshsplit 
%			should be !!
% threshmerge:  if two centroids are at a distance which is smaller than 
%		the threshold, then merge them
% The number of cluster centers is revised before stopping
%
% The strategy for merging is quite stupid
% We merge centers that are closer than a prefixed threshold,
% thus, a better strategy should be devised.
% __________________________
%  Author: Vittorio Castelli
%  Copyright IBM T. J. Watson Research Center
%      January 30, 1995;	Last modified: January 30, 1995

function clCen =  clmsead(origCen, trSet,prec, threshsplit, threshmerge);


lear_rate_orig = .7;
lear_rate_fin = .09;
lear_rate_dec = .09;
learning_rate = lear_rate_orig;
n_dist = 8;
maxcount = 64;
count = 0;
contatore = 0;
dis = 1;
firsttime = 1;
perc = 4;	% percentage of the training set shown in the plot
its1 = 9;
its2 = 18;	% indices of the coordinates of the pattern shown in
its3 = 27;	% plot3


[m,n] = size(trSet);
[n_cent,n1] = size(origCen);
maxcen = 3*n_cent;
if(its2 > n1)
   its2 = n1/2;
end

if(its3 > n1)
   its3 = n1;
end

if n ~= n1,
   errormes(4);
   return
end

aaa = ones(1,n);
aa = ones(1,m);
clCen = origCen;

centroid = mean(trSet);
origdisp = mean(sum((trSet'-centroid'*aa).^2))
splitdisp = origdisp*threshsplit;

figure(1);
clg;
figure(2);
clg;
figure(3);
clg;
asdfg = max(size(trSet));
plot3(trSet(1:perc:asdfg,its1),trSet(1:perc:asdfg,its2), ...
         trSet(1:perc:asdfg,its3),'.');
grid;
hold;
oldplCen = clCen;
plot3(clCen(:,its1),clCen(:,its2),clCen(:,its3),'g*');

contin =1;
disper = ones(1,n_cent)*1024;
oldisper = disper;


while contin==1,

  learning_rate =   (learning_rate-lear_rate_fin)* ...
			exp(-lear_rate_dec *contatore) + lear_rate_fin;

  oldCen=clCen;
  contatore = contatore +1;
  olddisper = disper;
  for ii = 1:n_cent,
      distan(ii,:) = sum((trSet'-((clCen(ii,:))'*aa)).^2);
  end;

  

  [distan,ord] = sort(distan);
  ord = ord';

  ord = ord(:,1);
  chancent=0;
  
  for i=1:n_cent,			% Calculates the number of pts. 
        nelcloud(i) = sum(ord==i);      % in the cloud. If it is > 0
        if nelcloud(i) ==0,		% it moves the center adaptively
             fprintf(' %g',i)		% with a greedy algorithm
  
	   % Calculate the new center, with a greedy algorithm
 	   % calculate the dispersion of the points around the 
	   % centers of the clouds
        else
           aver = sum((((ord==i)*aaa).*trSet))/nelcloud(i);
           clCen(i,:) = clCen(i,:).*(1-learning_rate)+aver.*learning_rate;
           disper(i) =  sum((ord == i).* ...
                   (sum(((clCen(i,:)'*aa)-trSet').^2))')/nelcloud(i);
        end
        d_i = sum((clCen(i,:)-oldCen(i,:)).^2);
        chancent=max(chancent, d_i);
        disi(i) = d_i;         
   end
 
   figure(1);
   subplot(211)
   plot(1:length(disper), sqrt(disper),'o',[1,length(disper)], ...
          [sqrt(splitdisp), sqrt(splitdisp)],'-');
   title('Dispersion');
   subplot(212)
   plot(disi,'x')
   titlestr = sprintf('Changes, iter. %g',contatore);
   title(titlestr);
   
   count = count+1;
   distance(count) = sum(nelcloud.*disper)/m;


   if rem(count,n_dist)==0, % Plot the maximum change in the centroids
      figure(2);
      plot(distance,'+');
      grid
      figure(3)
      if firsttime == 0,
		plot3(oldplCen(:,its1),oldplCen(:,its2),oldplCen(:,its3),'b+');
      end

      plot3(clCen(:,its1),clCen(:,its2),clCen(:,its3),'rx');
      oldplCen = clCen;
   end
   if count==maxcount,      % Plots at most maxcount changes
       distance(1:count-n_dist+1) = distance(n_dist:count);
       count = count-n_dist+1;
   end


   % termination conditions: small change in the centroids, 
   % or small change in MSE 

   if firsttime == 0,
      if chancent < prec,
         contin = 0;
         fprintf('\n Should stop, centroids');        
      end
      if max(abs(disper-oldisper))/max(abs(disper))<prec/100,
         contin =0;
         fprintf('\n Should stop: dispersion!');        
      end
   else
      firsttime = 0;
      oldisper = disper;
   end

   % Split and or merge

   if contin == 0,

      % split first      
        numcen = n_cent;
        [auxdisper ,inddisper] = sort(disper);
        i = inddisper(length(inddisper));

        if ( disper(i) > splitdisp) & numcen<maxcen, %then split cluster i !!
             fprintf('\n Aha, we are splitting cluster %g',i);        
             figure(3);		 % Mark the  center we split
             plot3(clCen(i,its1),clCen(i,its2),clCen(i,its3),'ro');
                                            
             diffmat = trSet - aa'*clCen(i,:);  % Calclulate cov. matrix
             diffmat = diffmat .*((ord==1)*aaa); % Identify direction of max
	     covmat = diffmat'*diffmat;		 % variation, split the 
	     [eivec, eival] = eig(covmat);	 % center
             eival = diag(eival);
	     [maxeig,maxpos] = max(eival);
             clCen(i,:) = clCen(i,:) + (eivec(:,maxpos))'.*sqrt(threshsplit);
             numcen = numcen+1;
             clCen(numcen,:) = clCen(i,:)-(eivec(:,maxpos))'.* ...
                                  sqrt(threshsplit); 
             contin = 1;
             firsttime = 1;
        end,

        distcen = ones(n_cent-1,n_cent-1)*(1+threshmerge);

        for i = 1:n_cent -1,
	  distcen(i,i:n_cent-1) = sum(((clCen(1:n_cent-i,:) - ...
                                       clCen(i+1:n_cent,:)).^2)');
        end

     % The strategy for merging is quite stupid
     % We merge centers that are closer than a prefixed threshold,
     % thus, a better strategy should be devised.

        [mindiscen,colu] = min(distcen);

        if min(mindiscen) < threshmerge,
            [cen1,cen2] = min(colu);      %find the closest centers
             fprintf('\n Aha, we are merging clusters %g, %g',cen1,cen2);
            figure(3);
            plot(clCen(cen1,1),clCen(cen1,2),clCen(cen1,3),'ro');
            plot(clCen(cen2,1),clCen(cen2,2),clCen(cen2,3),'ro');
            clCen(cen1,:) = (clCen(cen1,:)+clCen(cen2,:))/2;
            clCen(cen2,:) = clCen(numcen,:);
            numcen = numcen -1;
            contin =1 ;
            firsttime = 1;
        end;

        if n_cent < numcen
	   oldisper = ones(1,numcen)*threshmerge;
        end

        n_cent = numcen;     
        clCen = clCen(1:n_cent,:);
 	learning_rate = min(learning_rate*5,lear_rate_orig*.7);
   end; % if contin .. 
  
end
fprintf('\n Stop! \n');
figure(3)
hold;

figure(1)
clg;
asdfg = max(size(trSet));
plot3(trSet(1:perc:asdfg,its1),trSet(1:perc:asdfg,its2), ...
               trSet(1:perc:asdfg,its3),'.');
grid;
hold;
plot3(clCen(:,its1),clCen(:,its2),clCen(:,its3),'mx');
hold;

figure(3)