function [Centers,Variances,Partition]=qkgkclus(data,n,m,metric,U)
%MAN_PAGE_BEGIN
%
%   @purpose	Gustafson-Kessel clustering algorithm
%
%   @synopsis	qkgkclus(data,n,m,metric,U)
%   @description  perform a Gustafson-Kessel clustering.
% 	If N is the number of training patterns, k is the dimension of the input space
% 	and p is the number of rules::
% 
% 	data is a k*N matrix containing the vectors [v1 v2 v3 ... vN]
% 	n is the number of cluster
% 	m is the fuzzyness (1<m)
% 	U is the initial partition matrix (optional)
% 
% 	Centers is a k*p matrix containing the centers of the clusters [c1 c2 c3 ... cp]
% 	Variances is a k*(p.k) matrix containing the variances of the clusters [s1 s2 s3 ... sp]
% 	Partition is a p*N matrix containing the membership degrees of each vector to each clusters
%
%MAN_PAGE_END

% By Antoine Duchateau. All Rights Reserved.
% Written: 31/05/99


% Gustafson-Kessel clustering algorithm
% [Centers,Variances,Partition]=gkclus(data,n,m);
% 
% where (if N is the number of training patterns, k is the dimension of the input space
% and p is the number of rules)
% 
% data is a k*N matrix containing the vectors [v1 v2 v3 ... vN]
% n is the number of cluster
% m is the fuzzyness (1<m)
% U is the initial partition matrix (optional)
% 
% Centers is a k*p matrix containing the centers of the clusters [c1 c2 c3 ... cp]
% Variances is a k*(p.k) matrix containing the variances of the clusters [s1 s2 s3 ... sp]
% Partition is a p*N matrix containing the membership degrees of each vector to each clusters

if (nargin < 3), m = 2;    elseif isempty(m),   m = 2; end;
if (nargin < 4), metric = 1;    elseif isempty(m),   metric = 1; end;
% Initialisation

warning off;
treshold=0.01;
BadConditioning=1e-40;
dimension=size(data,1);
N=size(data,2);
for i=1:n,
	Variances(:,:,i)=eye(dimension);
end
TotalSteps=0;

% U is an alias of Partition

if (nargin < 5),
	UnivCenter=mean(data');
	UnivStd=std(data');
	Centers=randn(dimension,n).*(UnivStd'*ones(1,n))+UnivCenter'*ones(1,n);
	F=ones(dimension);
	
	for i=1:n,
		Delta=data-Centers(:,i)*ones(1,N);
		Dsquare(i,:)=sum(Delta.*Delta);
	end
	
	if n>1,
		for i=1:n,	
			U(i,:)=sum(((ones(n,1)*Dsquare(i,:))./Dsquare).^(1/(m-1))).^(-1);
		end
	else
		U=ones(1,N);
	end;
	
end
	
% Beginning of the algorithm
for iteration = 1:3,
	crit=1;
	while ((crit>treshold) & (TotalSteps<500)),
		TotalSteps=TotalSteps+1;
		OldU=U;
		Ubis=U.^m;
		Sums=sum(Ubis');
		
		% Compute de centers
		
		Centers=(data*Ubis')./(ones(dimension,1)*Sums);
		
		% Compute the metrics
		
		OldVariances=Variances;
		
		for i=1:n,
			Delta=data-Centers(:,i)*ones(1,N);
			
			switch iteration
				case 1,	% Constant case %
				F = eye(dimension) * sum(sum(Delta .* Delta) .* Ubis(i,:),2)/Sums(i);
				case 2,	% Diagonal case %
				F = diag(sum(Delta .* Delta .* Ubis(i*ones(dimension,1),:),2)/Sums(i));
				case 3, % General case %
				F=(((Ubis(i*ones(dimension,1),:)).*Delta)*Delta')/Sums(i);
			end
			
			if (det(F)<BadConditioning),
				F=OldVariances(:,:,i);
			end
			
			Variances(:,:,i)=F;
			
		% Compute de distances
		
			IF=inv(F);
			if (metric ==1),
				Dsquare(i,:)=sum((Delta'*((det(F)^(1/dimension))*IF))'.*Delta);
			else
				Dsquare(i,:)=sum(Delta.^2);
			end
			
		end
		
		% Compute the new partition
		
		if (min(min(Dsquare))>BadConditioning),
			for i=1:n,
				U(i,:)=sum(((ones(n,1)*Dsquare(i,:))./Dsquare).^(1/(m-1))).^(-1);
			end
		else
			% We can't use the quick method. We have to perform a double loop.
			for j=1:N,
				if (min(Dsquare(:,j))>BadConditioning),
					for i=1:n,
						U(i,j)=1/sum(((ones(n,1)*Dsquare(i,j))./Dsquare(:,j)).^(1/(m-1)));
					end
				else
					for i=1:n,
						U(i,j)=(Dsquare(i,j)==0);
					end
					U(:,j)=U(:,j)/sum(U:j);
				end
			end
		end
		crit=abs(norm(OldU-U));
		myMean=sum(sum(abs(OldU-U)))/(size(U,1)*size(U,2));
	end
end

Partition=U;	
