function [Z,order,Md] = som_linkage(sM,varargin)

%SOM_LINKAGE Make a hierarchical linkage of the SOM map units.
%
% [Z,order,Dist] = som_linkage(sM, [[argID,] value, ...])
%  
%  Z = som_linkage(sM);
%  Z = som_linkage(D,'complete');
%  Z = som_linkage(sM,'single','ignore',find(~som_hits(sM,D)));
%  Z = som_linkage(sM,pdist(sM.codebook,'mahal'));
%  som_dendrogram(Z); 
%
%  Input and output arguments ([]'s are optional):
%   sM       (struct) map or data struct to be clustered
%            (matrix) size dlen x dim, a data set: the matrix must not
%                     contain any NaN's!
%   [argID,  (string) See below. The values which are unambiguous can 
%    value]  (varies) be given without the preceeding argID.
%
%   Z        (matrix) size dlen-1 x 3, the linkage info
%                     Z(i,1) and Z(i,2) hold the indeces of clusters 
%                     combined on level i (starting from bottom). The new
%                     cluster has index dlen+i. The initial cluster 
%                     index of each unit is its linear index in the 
%                     original data matrix. Z(i,3) is the distance
%                     between the combined clusters. See LINKAGE
%                     function in the Statistics Toolbox.
%                     The ignored samples are listed at the 
%                     end of the Z-matrix and have Z(*,3) == Inf
%   Dist     (matrix) size dlen x dlen, pairwise distance matrix
%
% Here are the valid argument IDs and corresponding values. The values 
% which are unambiguous (marked with '*') can be given without the
% preceeding argID.
%   'linkage' *(string) the linkage criteria to use: 'single' (the
%                       default), 'average' or 'complete' 
%   'topol'   *(struct) topology struct
%   'connect' *(string) 'neighbors' or 'any' (default), whether the
%                       connections should be allowed only between 
%                       neighbors or between any vectors
%              (matrix) size dlen x dlen indicating the connections
%                       between vectors
%              (scalar) the N-neighborhood upto which the connections
%                       should be formed (implies 'neighbors')
%   'ignore'   (vector) the units/vectors which should be ignored 
%   'dist'     (matrix) size dlen x dlen, pairwise distance matrix to 
%                       be used instead of euclidian distances
%              (vector) as the output of PDIST function
%              (scalar) distance norm to use (euclidian = 2)
%   'mask'     (vector) size dim x 1, the search mask used to 
%                       weight distance calculation. By default 
%                       sM.mask or a vector of ones is used.
%
% Note that together 'connect'='neighbors' and 'ignore' may form
% areas on the map which will never be connected: connections
% across the ignored map units simply do not exist.
%
% See also KMEANS_CLUSTERS, LINKAGE, PDIST, DENDROGRAM. 

% Copyright (c) 2000 by Juha Vesanto
% Contributed to SOM Toolbox on June 16th, 2000 by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/
 
% Version 2.0beta juuso 160600

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% input arguments

% the data
if isstruct(sM), 
  if isfield(sM,'data'), D = sM.data; sT = []; mask = []; 
  else D = sM.codebook; sT = sM.topol; mask = sM.mask;
  end
else
  D = sM; sT = []; mask = []; 
end
[dlen dim] = size(D);
if isempty(mask), mask = ones(dim,1); end
if any(isnan(D(:))), error('Data matrix must not have any NaNs.'); end

% varargin
Md = 2; 
linkage = 'single';
ignore_units = []; 
constrained = 0;
i=1; 
while i<=length(varargin), 
  argok = 1; 
  if ischar(varargin{i}), 
    switch varargin{i}, 
      % argument IDs
     case {'topol','som_topol','sTopol'}, i=i+1; sT = varargin{i};
     case 'connect', i=i+1; 
       if ischar(varargin{i}), constrained = ~strcmp(varargin{i},'any');
       else constrained = varargin{i}; end
     case 'ignore',  i=i+1; ignore_units = varargin{i}; 
     case 'dist',    i=i+1; Md = varargin{i};
     case 'linkage', i=i+1; linkage = varargin{i};
     case 'mask',    i=i+1; mask = varargin{i};
     case 'tracking',i=i+1; tracking = varargin{i}; 
      % unambiguous values
     case 'neighbors', constrained = 1; 
     case 'any',       constrained = 0; 
     case {'single','average','complete'}, linkage = varargin{i};
     otherwise argok=0; 
    end
  elseif isstruct(varargin{i}) & isfield(varargin{i},'type'), 
    switch varargin{i}(1).type, 
     case 'som_topol', sTopol = varargin{i}; 
     otherwise argok=0; 
    end
  else
    argok = 0; 
  end
  if ~argok, 
    disp(['(som_linkage) Ignoring invalid argument #' num2str(i+1)]); 
  end
  i = i+1; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% distance matrix

% given distance matrix % jh corrected this place totally 27.3. 03
if (prod(size(Md))==1), % no explicit distance matrix, set flag
  q=2; % 17.2.03 kr added some brackets
else
  if (prod(size(Md))<dlen^2), % check pdist form
    Md = squareform(Md);   % transform to ordinary square diastance matrix
  end
  % jh: 27.3. 03 "calculate pairwise dist. matrix" see approx. 20 lines below
  % sets self-distance to Inf! This must be set here also,
  % otherwise clustering fails for user defined dist. matrix! 
  Md(eye(dlen)==1)=Inf; 
end

% neighborhood constraint
if length(constrained)==1 & constrained>0, 
  Ne1 = som_unit_neighs(sT); 
  Conn = som_neighborhood(Ne1,constrained); 
  Conn(~isfinite(Conn(:))) = 0; 
else Conn = constrained; end
if ~isempty(Conn), for i=1:dlen, C(i,i) = 1; end, end

% pairwise distance matrix across connected units
n = size(D,1);
if prod(size(Md))>1,   
  % remove distances between non-neighbors
  if constrained, for i = 1:n, Md(i,find(Conn(i,:)==0)) = Inf; end, end
else    
  % calculate pairwise distance matrix
  q = Md; 
  Md = zeros(n,n)+Inf;
  if ~constrained & q==2, % fast for the default case 
    for i = 1:n-1,
      x = D(i,:);
      inds = [(i+1):n]; 
      Diff = D(inds,:) - x(ones(n-i,1),:);
      Md(inds,i) = sqrt((Diff.^2)*mask);
      Md(i,inds) = Md(inds,i)';
    end  
  else
    for i = 1:n-1, 
      inds = find(Conn(i,:)==1); 
      inds = inds(find(inds>i)); 
      Diff = abs(D(inds,:) - D(i*ones(length(inds),1),:));
      switch q, 
       case 1,    dist = Diff*mask;
       case 2,    dist = sqrt((Diff.^2)*mask);
       case Inf,  dist = max(Diff,[],2);
       otherwise, dist = ((Diff.^q)*mask).^(1/q);
      end
      Md(inds,i) = dist; 
      Md(i,inds) = dist'; 
    end
  end
end

% set distances to ignored units to Inf
if ~isempty(ignore_units), 
  Md(ignore_units,:) = Inf;
  Md(:,ignore_units) = Inf;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% construct dendrogram

Z = zeros(n-1,3)+NaN;  % merged clusters and distance for each step
clusters = 1:dlen;     % each vector is at first in its own cluster
Cd = Md;               % distances between clusters

h = waitbar(0,'Constructing hierarchical clustering'); 
  
for i=1:n-1,   

  % tracking
  waitbar(i/(n-1),h); 

  %% combine two closest clusters  
  % find the clusters which are closest to each other (c1 and c2)
  [d,ind] = min(min(Cd));  
  if ~isfinite(d), break; end   % no more connected clusters
  [d,c1] = min(Cd(:,ind));      % cluster1
  c2 = clusters(ind);           % cluster2    
  % combine the two clusters
  c1_inds = find(clusters==c1); % vectors belonging to c1
  c2_inds = find(clusters==c2); % vectors belonging to c2
  c_inds = [c1_inds, c2_inds];  % members of the new cluster  
  % new cluster index = bigger cluster
  if length(c2_inds)>length(c1_inds), c=c2; k=c1; else c=c1; k=c2; end
  clusters(c_inds) = c;         % update cluster info
  Z(i,:) = [c, k, d];           % save info into Z
  
  %% update cluster distances  
  % remove the subclusters from the Cd table  
  Cd(c_inds,c_inds) = Inf;      % distance of the cluster to its members = Inf
  k_inds = c_inds(c_inds ~= c); % vectors of the smaller cluster
  Cd(k_inds,:) = Inf;           % distance of the subclusters to 
  Cd(:,k_inds) = Inf;           % other clusters = Inf
  % update the distance of this cluster to the other clusters
  cl = unique(clusters(clusters ~= c)); % indeces of all other clusters
  if ~isempty(cl), % added since v6.5 works differently than 6.1
    for l=cl,   
      o_inds = find(clusters==l); % indeces belonging to cluster k
      vd = Md(c_inds,o_inds);     % distances between vectors in c and k
      vd = vd(isfinite(vd(:)));   % remove infinite distances (no connection)
      len = length(vd);
      if ~len, % if the two clusters are not connected, their distance in Inf
	sd = Inf;
      else   % otherwise, calculate the distance between them
	switch linkage,
	 case 'single',   sd = min(vd);
	 case 'average',  sd = sum(vd)/len; 
	 case 'complete', sd = max(vd);
	 otherwise, error(['Unknown set distance: ' linkage]);
	end
      end
      Cd(c,l) = sd; Cd(l,c) = sd;
    end  
  end
end
close(h); 

last = Z(i,1); 
if isnan(last), 
  last = Z(i-1,1); 
  rest = setdiff(unique(clusters),last); 
  Z(i:n-1,1) = rest'; 
  Z(i:n-1,2) = last; 
  Z(i:n-1,3) = Inf; 
  i = i - 1; 
else 
  rest = []; 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% return values

% calculate the order of samples
order = last; 
% go through the combination matrix from top to down
for k=i:-1:1, 
  c = Z(k,1); k = Z(k,2);                     % what (k) change to what (c)
  j = find(order==c);                         % the occurance of c in order    
  if j == length(order), order = [order k];   % put k behind c
  else order = [order(1:j) k order(j+1:end)]; 
  end
end  
order = [rest, order]; 

% to maintain compatibility with Statistics Toolbox, the values in 
% Z must be yet transformed so that they are similar to the output
% of LINKAGE function

Zs = Z;
current_cluster = 1:dlen;
for i=1:size(Z,1),
  Zs(i,1) = current_cluster(Z(i,1));
  Zs(i,2) = current_cluster(Z(i,2));
  current_cluster(Z(i,[1 2])) = dlen+i;  
end

Z = Zs;

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%