function [Centers,Variances,Partition]=gkclus(data,n,m,metric,U);
% 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

treshold=0.001;
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
	
crit=1;

% Beginning of the algorithm

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);
		F=(((ones(dimension,1)*Ubis(i,:)).*Delta)*Delta')/Sums(i);
		
		theLim=1e-5;
		while (det(F)<BadConditioning),
			[eigVect,eigVal]=eig(F);
			eigVal=diag(eigVal);
			neweigVal=(eigVal<theLim)*theLim+(eigVal>=theLim).*eigVal;
			F=eigVect*diag(neweigVal)*eigVect';
			fprintf('.');
			theLim=theLim*10;
		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.
		fprintf('Nul distance\n');
		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));
	fprintf('Norme is %12.8f and Mean is %12.8f\n',crit,myMean);
end

Partition=U;	
