function [f,v,P,R] = gkvol(x,c,m,e,sim,s)
% GK clustering with volume prototypes and 
% similarity driven cluster merging
%
% [F,V,P,R]  = GKVOL(X,C,m,e,sim,s)
%     Input:
%       X  ... M by N data matrix,
%              M is the number of data points and N data dimension
%       C  ... the initial number of clusters or the initial
%              cluster centers ([CxM]). Given an upper estimate
%              the merging routine will reduce the number of clusters.
%       m  ... optional parameter m > 1, determines the fuzziness
%              of clustering, for m close to 1 clusters become crisp
%              default value is 2
%       e  ... optional termination tolerance, the algorithm stops
%              when |F(k-1) - F(k)| <= e, default tolerance 1e-3
%       sim ... optional merging threshold 0 <= sim <=1
%              clusters are merged when their similarity is > sim, 
%              Special values:
%              sim == -1  --> ADAPTIVE, 1/(NC-1) (default)
%              sim == -2  --> ADAPTIVE, 1/NC
%              sim == -3  --> ADAPTIVE, 1/(NC+1)
%
%       s  ... (optional) parameter for plotting intermediate results
%              (projection onto 2-dimensions), default 0 - no plot, set
%              to 1 to view the clustering process on-line.
%     Output:
%       F  ... fuzzy partition matrix
%       V  ... centroid prototype matrix
%       P  ... cluster covariance matrices
%       R  ... cluster prototype radius

% (c) Kaymak & Setnes, March 1998 
% Setnes, Aug 1999, oct 1999
% Kaymak  Oct 1999
% Setnes, Oct 18, 1999
% Setnes, Jan 2000
% Kaymak, Jul 2000

if (nargin < 3), m = 2; elseif isempty(m), m = 2; end;
if (nargin < 4), e = 1e-3; elseif isempty(e), e = 1e-3; end;
if (nargin < 5), sim = -1; elseif isempty(sim), sim = -1; end;
if (nargin < 6), s=0;elseif isempty(s), s = 0; end;

%flag for adaptive similarity
if (sim == -1), 
   ADAPT=1;
   disp('Cluster merging with adaptive similarity 1/(NC-1)');
elseif (sim == -2),
   ADAPT=2;
   disp('Cluster merging with adaptive similarity 1/NC');
elseif (sim == -3 ),
   ADAPT=3;
   disp('Cluster merging with adaptive similarity 1/(NC+1)');
else
   ADAPT=0;
end;

if(size(x,1)<size(x,2)),x=x';end;
if(size(c,1)<size(c,2)),c=c';end;

[mx,nx] = size(x);
x1 = ones(mx,1);

%flag=0;
%if(nx<2),%IF INPUT DIM = 1, MAKE FAKE 2'ND DIMENSION
%   nx=2;x=[x,x1*0];
%   if(size(c,1)>1),c=[c,c*0];end;
%   flag=1;%SET FLAG input dim = 1
%end;

%% INITIALIZE FUZZY PARTITION MATRIX %%
if( max(size(c)) > 1 ),	% initial cluster centers assumed given when c is not a scalar
   v=c(:,1:nx);
   
   % add small amount of noise to prevent centers falling on data points
   v = v + e^2 * randn(size(v));
   
   c=size(c,1); % c is now number of clusters
   
   % Calculate f0 from the given centers with Euclidian distance norm
   for j = c:-1:1,    % for all clusters
      xv = x - x1*v(j,:);
      d(:,j) = sum([0*x1 (xv.^2)]')';
   end;
   d = d.^(-1/(m-1));
   f0 = (d ./ (sum(d')'*ones(1,c)));
   
else;					% Only initial number of clusters given
%   if(size(c,1)>1),	% initial cluster centers wrongly given with dim < dim of cluster space
%      c=size(c,1);	% extract number of wrongly given centers
%   end;
   f0  = zeros(c,mx);
   for j=1:mx
      cc = 1 + fix((j-1)/round(mx/c));
      while (cc > c),
         cc = cc - 1;
      end;
      f0(cc,j) = 0.6+0.4*rand;
      for i=1:c,      
         sumf0 = sum(f0(:,j));
         if (i ~= cc), 
            f0(i,j) = 1-sumf0;
            if (i ~= c),
               f0(i,j) = f0(i,j)*rand;
            end;
         end;
         if (c~=1 & cc==c & i==c),
            f0(i-1,j) = f0(i-1,j) + 1 - sumf0;      
         end;
      end;
   end;
   f0 = (f0.*(f0 > 0))';
   fm = f0.^m; sumf = sum(fm);
   v = (fm'*x)./(sumf'*ones(1,nx));	%calculate initial centers
end;

f = zeros(mx,c);                % partition matrix

%%% INITIALIZE GRAPHICS %%%
if (s ~= 0),
   figure(1);clf;
   if ( nx == 1 ), % Make a fake dimension for 2-D display
      vplot = [v(:,1) zeros(c,1)];
      xplot = [x(:,1) zeros(mx,1)];
   else
      vplot = v(:,1:2);
      xplot = x(:,1:2);
   end;
   
	H = plot(xplot(:,1),xplot(:,2),'k.',vplot(:,1),vplot(:,2),'ro','EraseMode','xor');
   minx = min(xplot(:,1)); maxx = max(xplot(:,1));
   miny = min(xplot(:,2)); maxy = max(xplot(:,2));
   if(miny==maxy),
       miny=miny-.01;
       maxy=maxy+.01;
   end;
   if(minx==maxx),
       minx=minx-.01;
       maxx=maxx+.01;
    end;
   axis([minx-.1*abs(maxx) maxx+.1*abs(maxx) miny-.1*abs(maxy) maxy+.1*abs(maxy)])
   set(H(2),'LineWidth',3);
   set(gcf,'UserData',H);
   drawnow;
end;

%%% ITERATE %%%
iter=0; ms0 = 1;
volfact=1; % scaling factor for the cluster volume (1 <= volfact <= c)
while( (max(max( abs(f0-f) > e ))) | (iter==0)),	%iter=0 means merging has just happened
   iter=iter+1;
   f = f0;
   
   % Calculate centers
   fm = f.^m; sumf = sum(fm);
   v = (fm'*x)./(sumf'*ones(1,nx));
   
   % Calculate cluster covariances
   P=[]; 
   for j = 1 : c,
      xv = x - x1*v(j,:);
      %covariance matrixes
      p = ones(nx,1)*fm(:,j)'.*xv'*xv/sumf(j);
      P((j-1)*nx+(1:nx),:) = p;
   end;
   
%   % Determine intercluster distances
%   % R cannot be larger than the distance to the nearest cluster
%   cldist = sqrt(gkdist(v,v,P)) + diag( inf*(1:size(v,1))); % distance to self is infinite
%   nearest = min(cldist);
   
%	% Distill main directions in P
   % Determine cluster radii from cluster volumes
	[nMatrix,sMatrix] = size(v);
%	% Get the cluster variances for volume centers
   R=[];
	for dummy = 1 : nMatrix,
	    dummyP = P(((dummy-1)*sMatrix+1):(dummy*sMatrix),:);
%       %     dummyP = max(diag(dummyP)) * eye(size(dummyP));
%       % dp = dummyP;
%       %dp = det(dummyP).^(1/nx)*inv(dummyP);
%       %[Vd,Dd] = eig(dp);
%       %Vd = Vd*sqrt(Dd);    %scale vectors
%       %Vd = Vd*sqrt(sqrt(inv(Dd)));
%       %dD = sqrt(gkdist(Vd',zeros(1,size(Dd,2)),dummyP)); %GKDIST returns d^2
%        %R = [R;dD(1)];  % all elements of dD are the same
%       %radius = sqrt(det(dummyP)^(1/nx)/c);
%       val1 = (1/sqrt(c))*nearest(dummy); val2=sqrt((det(dummyP)^(1/nx))/(c/volt));
%       radius=val2;
%       %tmp=[val1, val2*sqrt(c), val2*sqrt(iter)]; radius=min(tmp);
%       %radius = min( val1, sqrt(c)*val2);
%      radius = min( (1/c)*nearest(dummy), sqrt(det(dummyP)^(1/nx)));
       radius = sqrt((det(dummyP)^(1/nx))/(c/volfact));
       R = [R;radius];
%       % R=[R;1];
	end;
  
%   %% Avoid overlapping centers
%   temp=diag(inf*R);
%   temp=sqrt(gkdist(v,v,P))+temp;	%DIST returns d
%   temp=(min(temp)')/c; %maximum allowed radius
%   R=min(R,temp);
%R = 0*R;

%% Calculate distances
   d=sqrt(gkdist(x,v,P));         % Mahalonibis distance to cluster centers
   d = d-(ones(size(d,1),1)*R');	% Adjust for volume prototypes
	tt=(1e-8)*ones(size(d));		% Numerical stability
	d=max(d,tt);
   
%   % Alternative calculation for distances
%   %d = sqrt(gkdist(x,v,P));
%   %dd = d - (ones(size(d,1),1)*R');
%   %[mind,minx] = min(dd');
%   %tochange = find(mind < 0);
%   %d = dd;
%   %if ( ~isempty( tochange ) ),
%   %   subd = sub2ind(size(dd),tochange,minx(tochange));
%   %   
%   %   % First set all negative distances to 1.0 (a dummy value)
%   %   d(find(dd<0)) = ones(size(find(dd<0)));
%   %   d(subd) = 1e-8*ones(size(subd));
%   %end;
   
   %% Update f0
   d = d.^(2/(m-1));
   sumd = sum((d.^(-1))')'*ones(1,c);
   f0 = (d .* sumd).^(-1);
  
	%% Plot intermediate results
   if (s ~= 0),
      if ( nx == 1 ), % Make a fake dimension for 2-D display
         vplot = [v(:,1) zeros(c,1)];
      else
         vplot = v(:,1:2);
      end;
      H = get(gcf,'UserData');
      set(H(2),'xdata',vplot(:,1),'ydata',vplot(:,2)); drawnow;
   end;

	%% SIMILARITY MEASURE
   S=[]; 
%   %hvar=sum(f0)
  	for ii=1:(size(f0,2))
		for jj=1:(size(f0,2))
			if(ii==jj),
				S(ii,jj)=0;
         else
            S(ii,jj)=sum(min(f0(:,ii),f0(:,jj)))/sum(f0(:,ii));	   %inclusion
%            
%            % Do not merge if volumes are not similar
%            %if( abs(hvar(ii)/hvar(jj) - 1) > 0.2 ), S(ii,jj)=0.0; end;
%               
			end;
		end;
	end;

	%%% MERGING
   ms=max(S(:));
   
   %Adaptive similarity measures
   if(ADAPT == 1),
      sim=1/(c-1);	   
   elseif (ADAPT == 2),
      sim=1/(c);
   elseif (ADAPT == 3),
      sim=1/(c+1);
%     %sim=(1/c)^(1/size(x,2));
%     %sim=1/sqrt(((c^2)-c)); %adapt Uzay
%     %a=1/c;b=2/(c^2-c);sim=(b/nx + a*(nx-1))/nx;	% Adapt setnes
   end
   
   % Consider processing only if the memberships have sufficiently stabilized
%   if ( min(min( abs(f0-f) < 100*e )) ),
   if ( abs(ms0-ms) < e ),
      % Merge clusters if similarity > threshold
      if ( (ms>sim) & (c>2) ),
        iter=0; 
%        %volfact=1; % Reset after merge
         
        % Find the kandidates for merging
		 [ii,jj]=find(S==ms);
        % There may still be more than one candidate. Use Jaccard to select
		 S=[];
        for i=1:size(ii,1),
           S(i,1)= sum(min(f0(:,ii(i)),f0(:,jj(i))))/sum(max(f0(:,ii(i)),f0(:,jj(i))));	%jaccard
        end;
        i = find( S==max(S) );
        ii=ii(1);jj=jj(1);
        
   	     % merge the most similar cluster pair
	  	 f(:,ii) = f(:,ii) + f(:,jj);
        f(:,jj) = [];
  		 f0(:,ii) = f0(:,ii) + f0(:,jj);
        f0(:,jj) = [];
        c=c-1;
        
        % Reset similarity measure
        ms = 1;
        disp(['Number of clusters: ' num2str(c)]);
     else
        % You cannot merge but perhaps you can increase cluster volume
        volfact = min( c, volfact+1 );
     end;
  end;
  
  % Store maximum similarity
  ms0 = ms;
end;

% Clustering finished. Return iteration outcome
f = f0;

%S

%%% PLOT volume centers %%%
if ( (s ~= 0) ),
   hold on;
   if ( nx == 1 ), % Make a fake dimension for 2-D display
      vplot = [v(:,1) zeros(c,1)];
   else
      vplot = v(:,1:2);
   end;
   
   % Plot volume prototypes
   elps( vplot(:,1:2), P, R.^2 );   
     
   % plot also a coloured output
   % figure(2); colord = ['ymcrgbk'];
   % ff = f==(max(f')'*ones(1,size(f,2)));
   % for ( k = 1:size(ff,2) ),
   %   indx = find(ff(:,k) == 1); 
   %   plot(x(indx,1),x(indx,2),[colord(mod(k,7)+1) 'o']); grid; hold on;
   % end;
   
end;

%if(flag); %input dim = 1
%   v=v(:,1);
%end;

%%% End of function FCMVOL %%%


function d = gkdist(x,v,P)
% (c) Robert Babuska, April 1994
%     Uzay Kaymak, 2000
[mx,nx] = size(x);
[mv,nv] = size(v);
if (nx == 1),
   xv = x*ones(1,mv) - ones(mx,1)*v';
   d = xv.^2;
else 
   for j = mv : -1 : 1,
      p = P((j-1)*nx+(1:nx),:);
      M = det(p).^(1/nx)*inv(p);
      xv = x - ones(mx,1)*v(j,:);
      d(:,j) = sum((xv*M.*xv)')';
   end;
end;

%function d = dist(x,v)
%% (c) Robert Babuska, June 1993
%[mx,nx] = size(x);
%[mv,nv] = size(v);
%if (nx == 1)
%  d = abs(x*ones(1,mv) - ones(mx,1)*v');
%else
%  d = zeros(mx,mv);
%  for i = 1 : mv
%      d(:,i) = sqrt(sum((x - ones(mx,1)*v(i,:))'.^2)');
%  end;
%end;

%function circ(xc,yc,rc)
%% (c) Magne, 1999
%x=[];y=[];
%for i =1:length(xc),
%	x=[x,[-rc(i):(2*rc(i))./100:rc(i)]'];
%end;
%y=[([ones(length(x),1)*rc'].^2-x.^2).^(1/2)];
%y1=-y;
%plot(x+[ones(length(x),1)*xc'],y+[ones(length(x),1)*yc'],'k',x+[ones(length(x),1)*xc'],y1+[ones(length(x),1)*yc'],'k')

function elps(v,P,R)
% (c) Uzay Kaymak, 2000
[nm,dummy] = size(v);
[dummy,sm] = size(P);
phi = linspace(0,2*pi);
for cnt = 1:nm,
   if ( sm == 1 ),
      a = sqrt(R(cnt));
      b = a;
      Vdd = eye(2,2);
   else
      dp = P(((cnt-1)*sm+1):(cnt*sm),:);
      dp = det(dp).^(1/sm)*inv(dp);
      [Vdd,Dd] = eig(dp);
      a = sqrt(R(cnt)/Dd(1,1)); 
      b = sqrt(R(cnt)/Dd(2,2));
   end;
   
   % Plot the elliptic volume prototypes
   X = [a*cos(phi);b*sin(phi)];
   ec = (Vdd(1:2,1:2)*X)';
   XC = v(cnt,1) + ec(:,1);
   YC = v(cnt,2) + ec(:,2);
   plot(XC,YC);
end;
