function [PotSetsHierarchy,PotSetsMetric,PotSets]=gkreduc(PotNodes,PotNodesMetric,PotNodesPartition,Treshold,m,data,aexp,bexp,cexp);
% Reduction of the number of clusters
% 
% Syntax:
%
% [ClusSetsHierarchy,ClusSetsMetric,ClusSets]=gkreduc(Centers,Variances,Partition,Threshold,m,data,aexp,bexp,cexp)
%
% Inputs:
% 
% Centers is a k*n matrix containing the centers of the clusters [c1 c2 c3 ... cn]
% Variances is a k*(n.k) matrix containing the variances of the clusters [s1 s2 s3 ... sn]
% Partition is a n*N matrix containing the membership degrees of each vector to each cluster
% Threshold is the minimal number of clusters
% m is scalar representing the fuzziness of the partition
% data is a k*N matrix
% aexp is a scalar representing the relative importance of the parallelism criterion
% bexp is a scalar representing the relative importance of the closeness criterion
% cexp is a scalar representing the relative importance of the colinearity criterion
%
% Outputs:
%
% ClusSetsHierarchy is a matrix containing the hierarchy relationships between the clusters
% ClusSetsMetric is a matrix containing the variances of the merged clusters
% ClusSets is a matrix containing the centers of the merged clusters

if (nargin < 7), aexp = 1;    elseif isempty(m),   aexp = 1; end;
if (nargin < 8), bexp = 3;    elseif isempty(m),   bexp = 3; end;
if (nargin < 9), cexp = 2;    elseif isempty(m),   cexp = 2; end;

dimension=size(PotNodes,1);
numrules=size(PotNodes,2);
PotSetsHierarchy=[];
PotSets=[];
PotSetsMetric=[];
Limit=0;

crit=1;
while ((numrules>Treshold)&(crit>Limit)),
	[SecondSetsMetrics,SecondSets,IDs]=get1lev(PotSetsHierarchy,PotSetsMetric,PotSets,PotNodesMetric,PotNodes);
	
	% S1 tient compte du parallelisme
	
	Phi=[];
	for i=1:size(SecondSets,2),
		[Vectors,values]=eig(SecondSetsMetrics(:,:,i));
		[Minvalue,MinIndex]=min(abs(diag(values)));
		Phi(:,i)=Vectors(:,MinIndex);
	end;
	S1=abs(Phi'*Phi);
	
	% S2 tient compte de la proximite
	S2=zeros(size(SecondSets,2));
	for i=1:size(SecondSets,2),
		Psi=SecondSets(:,i)*ones(1,size(SecondSets,2));
		S2(i,:)=sum((Psi-SecondSets).^2);
	end;
	S2=sqrt(S2);
	
	% S3 tient compte de la colinearite
	S3=zeros(size(SecondSets,2));
	for i=1:size(SecondSets,2),
		dPsi=SecondSets(:,i)*ones(1,size(SecondSets,2))-SecondSets;
		NormDist=sqrt(sum(dPsi.^2));
		NormDist=NormDist+(NormDist==0);
		S3(i,:)=(abs(sum((Phi(:,i)*ones(1,size(SecondSets,2))).*dPsi)./NormDist)+abs(sum(Phi.*dPsi)./NormDist))/2;
	end;
	
	a=1/(size(SecondSets,2)*(size(SecondSets,2)-1))*sum(sum(S1-diag(diag(S1))));
	b=1/(size(SecondSets,2)*(size(SecondSets,2)-1))*sum(sum(S2-diag(diag(S2))));
	c=1/(size(SecondSets,2)*(size(SecondSets,2)-1))*sum(sum(S3-diag(diag(S3))));
	S1=S1.*(S1>a)+(S1<=a)*a;
	S2=S2.*(S2<b)+(S2>=b)*b;
	S3=S3.*(S3<c)+(S3>=c)*c;
	S1=(S1-a)/(1-a)*0.9+0.1;
	S2=1-0.9*S2/b;
	S3=1-0.9*S3/c;
	
	S123=triu(((S1.^aexp).*(S2.^bexp).*(S3.^cexp)).^(1/(aexp+bexp+cexp)),1);
	
	[MaxColS123,iTable]=max(S123);
	[MaxLineS123,j]=max(MaxColS123);
	
	if (MaxLineS123>Limit),
		i=iTable(j);
		% If S123=0, on peut se retrouver avec l'element 1 1
		if ~(i==j),
			numrules=numrules-1;
			[PotSetsHierarchy,PotSetsMetric,PotSets]=gkmerge(PotNodes,PotNodesMetric,PotNodesPartition,...
			PotSetsHierarchy,PotSetsMetric,PotSets,i,j,IDs,m,data);
		end
	end
	
	crit=MaxLineS123;
end
