function [f,c,P,V,D,J,M,fnc,S,hist,DT,DAV,FHV,PD,APD,F,FN] = amerge(x,f0,m,e,s,thr,flags)
% Compatible Cluster Merging algorithm for determining the optimal
% number of clusters in Gustafson-Kessel clustering. This algorithm 
% can be used for automatic generation of Takagi-Sugeno 
% fuzzy models from data.
%
% [F,C,P,V,D,J,M,FNC,S,HIST] = AMERGE(X,F0,m,e,s,THR,FLAGS)
%     Input:
%       X    ... M by N data matrix for GK-clustering,
%                M is the number of data points and N data dimension
%       F0   ... either an initial fuzzy partition matrix, or the
%                maximum number of clusters. In the latter case, 
%                a default partition matrix is generated.
%       m    ... m > 1, determines the fuzziness of clustering, 
%                for m close to 1 clusters become crisp
%       e    ... termination tolerance, the algorithm stops
%                when max(max(F(k-1) - F(k))) <= e
%       s    ... parameter for plotting intermediate results
%                (only for 2D data), set to 0 for no plot, set
%                to 1 to show the clustering process on-line
%       THR  ... threshold for merging. IF AUTO = 0, this threshold is
%                used for all merges. IF AUTO = 1, this is the minimum
%                threshold (i.e. the stopping threshold).
%       FLAGS .. optional flags vector for switching on various attributes 
%                of the AMERGE function
%                FLAGS(1): optional flag for automatic determination of 
%                          the merging threshold. Default 0 : use the 
%                          threshold supplied by the user, set to 1 for 
%                          automatic determination of the merging threshold.
%			   or to 2 to merge one cluster at a time
%                FLAGS(2): optional flag for switching off the heuristic
%                          step in the AMERGE algorithm. Default 1: use
%                          the heuristic step. Set to 0 for switching
%                          off the heuristic step.
%
%     Output:
%       F    ... final fuzzy partition matrix
%       C    ... final cluster means matrix
%       P    ... cluster covariance matrices from the final clustering step,
%                concatenated in one matrix: 
%                P = [P1;P2;...Pk], where k is number of clusters
%       V    ... eigenvectors of the covariance matrices from the final 
%                clustering step, corresponding to the smallest eigenvalues
%       D    ... eigenvalues of the covariance matrices from the final
%                clustering step
%       J    ... history of the clustering criterion J
%       M    ... matrices inducing the distance norm, calculated as
%                M = det(P).^(1/N)*inv(P), M = [M1;M2;...Mk]
%       FNC  ... final number of clusters
%       S    ... (upper triangular) matrix of final cluster similarity
%       HIST ... optional history of the merging process in terms of the 
%                number of clusters and the threshold values

% (c) Uzay Kaymak, September 1995-6
%     modified by Robert Babuska Feb 1996
%	- standartization of the centres before using CCM criteria
%	- FLAG(1)=2 does merging of one pair at a time down to two clusters

% initialization
if (nargin < 7) | isempty(flags), flags = [0 1]; end;
auto = flags(1);
heuristic = flags(2);
if (length(f0) > 1), NCmax = size(f0,2); else NCmax = f0; end;
if (auto ~= 0), th_min = thr; else th_min = 0; end; % stopping criterion
NClogger = [];                                  % logger for number of clusters
THlogger = [];                                  % logger for merging threshold
NCmin = 2;                                      % minimum number of clusters
show = 0;

% Start the main merging loop
stopC = 0;                                      % stopping flag, 1 : stop
NC = NCmax;                                     % initialization of number of clusters
f = f0;
rho = ones(1,NC);
nx = size(x,2);

while (NC >= NCmin) & ~stopC,                   % merging loop
  NClogger = [NClogger NC];                         % log history
  if (s ~= 0), fprintf('Number of clusters = %d\n\n',NC); end;
  
  [f,c,P,V,D,J,M] = gkfast(x,f,m,e,s,rho);         % GK-clustering step
%  [tr,tr,tr,tr,detP] = fcov(x,f,c,m);

% validity measures

[trash,ind] = sort(c(:,1));
c = c(ind,:);
f = f(:,ind);
V = V(ind,:);
D = D(ind,:);

%  p = [V(:,1:NI) -sum((V(:,1:NI+1).*c(:,1:NI+1))')']./(-V(:,NI+1)*ones(1,NI+1));
%  [ym,yl,ylm] = sugval(p,[d(:,1:NI) ones(size(d,1),1)],f);
 
  %%%%%%%%%%%% Merge clusters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  oo = ones(size(c,1),1);			% standardize cluster centers
  c = (c-oo*mean(c))./(oo*std(c));
%  c = stand(c);
  [c1,c2,c3] = ccm2(V,D,c(:,1:nx));           % evaluate modified CCM criteria

if show,
  c1,pause
  c3,pause
end;
                                                % Here comes the ccmm part
  c1=c1.^10;                                    % re-scale first criterion
                                                % for better resolution
                                                
  v1=c1(find(triu(ones(size(c1)),1)));          % extract column of c1 values
  v2=c2(find(triu(ones(size(c2)),1)));          % extract column of c2 values
  v3=c3(find(triu(ones(size(c3)),1)));          % extract column of c3 values

  aa = genmean(v1,1,[]);                        % determine membership function limits
  cc = genmean(v3,1,[]);
    
  mf1=mgrade(v1',[2 aa 1 1 1]);                 % define and evaluate 
  mf3=mgrade(v3',[2 0 0 0 cc]);                 % membership functions

  matr = [mf1' mf3'];                           % form MC decision matrix

  res=zeros(size(c1));                          % form similarity matrix RES
  res(find(triu(ones(size(res)),1)))=genmean(matr',0.5,[])';
if show, res, pause;  end;

  if (auto == 1),                               % auto determination threshold
    x_thr = res(find(triu(res)));
    [thr1,v_thr] = fisofast(x_thr,2,2,0.01);
    th = max(v_thr);
th
  elseif auto == 0,                             % user defined threshold
    th = thr;
  else
    th = max(res(find(triu(res))));
th
  end; 
  THlogger = [THlogger th];                     % log merging threshold

  if th < th_min,                               % stop if th_min is reached
    stopC = 1;
  else                                          % continue with merging
    merge_mat=(eye(size(res)) + res' + res);    % form merge matrix
                                                % RES is only upper triangular
                                                
    if (merge_mat > 0) == eye(size(res)),       % stop if you cannot merge
      stopC = 1;
    else
                                                % start relational clustering
      help_mat_1 = merge_mat;                   % initialize
if show, merge_mat , pause; end;
      help_mat = sup_min(merge_mat,merge_mat);
    
      while any(any(help_mat ~= help_mat_1)),   % loop until convergence
        help_mat_1 = help_mat;
if show,         help_mat_1, pause; end;
        help_mat = sup_min(help_mat,merge_mat);
      end;

      help_mat = 1.0*(help_mat >= th);          % Result of relational clustering
                                                % after thresholding
      NCc=NC;                                   % determine which clusters to merge
      NC = rank(help_mat,1e-9);                 % new number of clusters
      help_mat = rref(help_mat,1e-9);           % reduced row echelon form
      help_mat = (help_mat(1:NC,:));            % extract non-zero rows


      ii = 0;                                   % Heuristic step to prevent over merging
      help_mat_1 = help_mat(1,:);
      for cnt = 1 : size(help_mat,1),
        c_h1 = max(max(c2(find(help_mat(cnt,:)),find(help_mat(cnt,:)))));
        c_h2 = min(max(c2(find(help_mat(cnt,:)),find(~help_mat(cnt,:)))));
      
        if (c_h2 <= c_h1) & (heuristic == 1),   % Heuristic criterion activated, so do not merge
          h_m = diag(help_mat(cnt,:));
          help_mat_1(ii+1:ii+sum(help_mat(cnt,:)),:) = h_m(find(help_mat(cnt,:)),:);
          ii = ii + sum(help_mat(cnt,:));  
        else
          ii = ii + 1;                          % No need for heuristics, so continue (i.e. merge)
          help_mat_1(ii,:) = help_mat(cnt,:);
        end;
      end;                                      % end of FOR
    
      NC = ii;                                  % number of clusters after merging
      help_mat = help_mat_1(1:NC,:)';
      if NC >= NCc,                             % stop if no clusters are merged
        stopC = 1;
      end;
      f = f * help_mat;                         % merge clusters
%      [f,c,P,V,D,J,M] = fmle(x,f,m,e,s);        % FMLE clustering
%      c = (f.^m'*x)./(sum(f.^m)'*ones(1,nx));	% calculate centers
%      [tr,tr,tr,tr,detP] = fcov(x,f,c,m);
%      rho = detP/min(detP)
%      rho = detP1/min(detP)
    end;                                        % merge_mat <> identity  
  end;                                          % (th < th_min)
disp('Press a key to continue'); pause
end;                                            % end of WHILE
%%%%%%%%%%%%%% end of merge %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%[f,c,P,V,D,J,M] = gkfast(x,f,m,e,s,rho);         % GK-clustering step

% f, c, P, V, D, J, M are now defined
fnc = NC;
S = res;
if (nargout > 9), hist = [NClogger;THlogger]; end;


