somtoolbox2/som_umat.m
4dbef185
 function U = som_umat(sMap, varargin)
 
 %SOM_UMAT Compute unified distance matrix of self-organizing map.
 %
 % U = som_umat(sMap, [argID, value, ...])
 %
 %  U = som_umat(sMap);  
 %  U = som_umat(M,sTopol,'median','mask',[1 1 0 1]);
 %
 %  Input and output arguments ([]'s are optional): 
 %   sMap     (struct) map struct or
 %            (matrix) the codebook matrix of the map
 %   [argID,  (string) See below. The values which are unambiguous can 
 %    value]  (varies) be given without the preceeding argID.
 %
 %   U        (matrix) u-matrix of the self-organizing map 
 %
 % Here are the valid argument IDs and corresponding values. The values which
 % are unambiguous (marked with '*') can be given without the preceeding argID.
 %   'mask'       (vector) size dim x 1, weighting factors for different 
 %                         components (same as BMU search mask)
 %   'msize'      (vector) map grid size
 %   'topol'     *(struct) topology struct
 %   'som_topol','sTopol' = 'topol'
 %   'lattice'   *(string) map lattice, 'hexa' or 'rect'
 %   'mode'      *(string) 'min','mean','median','max', default is 'median'
 %
 % NOTE! the U-matrix is always calculated for 'sheet'-shaped map and
 % the map grid must be at most 2-dimensional.
 % 
 % For more help, try 'type som_umat' or check out online documentation.
 % See also SOM_SHOW, SOM_CPLANE.
 
 %%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
 % som_umat
 %
 % PURPOSE
 %
 % Computes the unified distance matrix of a SOM.
 %
 % SYNTAX
 %
 %  U = som_umat(sM)  
 %  U = som_umat(...,'argID',value,...)
 %  U = som_umat(...,value,...)
 %
 % DESCRIPTION
 %
 % Compute and return the unified distance matrix of a SOM. 
 % For example a case of 5x1 -sized map:
 %            m(1) m(2) m(3) m(4) m(5)
 % where m(i) denotes one map unit. The u-matrix is a 9x1 vector:
 %    u(1) u(1,2) u(2) u(2,3) u(3) u(3,4) u(4) u(4,5) u(5) 
 % where u(i,j) is the distance between map units m(i) and m(j)
 % and u(k) is the mean (or minimum, maximum or median) of the 
 % surrounding values, e.g. u(3) = (u(2,3) + u(3,4))/2. 
 %
 % Note that the u-matrix is always calculated for 'sheet'-shaped map and
 % the map grid must be at most 2-dimensional.
 %
 % REFERENCES
 %
 % Ultsch, A., Siemon, H.P., "Kohonen's Self-Organizing Feature Maps
 %   for Exploratory Data Analysis", in Proc. of INNC'90,
 %   International Neural Network Conference, Dordrecht,
 %   Netherlands, 1990, pp. 305-308.
 % Kohonen, T., "Self-Organizing Map", 2nd ed., Springer-Verlag, 
 %    Berlin, 1995, pp. 117-119. 
 % Iivarinen, J., Kohonen, T., Kangas, J., Kaski, S., "Visualizing 
 %   the Clusters on the Self-Organizing Map", in proceedings of
 %   Conference on Artificial Intelligence Research in Finland,
 %   Helsinki, Finland, 1994, pp. 122-126.
 % Kraaijveld, M.A., Mao, J., Jain, A.K., "A Nonlinear Projection
 %   Method Based on Kohonen's Topology Preserving Maps", IEEE
 %   Transactions on Neural Networks, vol. 6, no. 3, 1995, pp. 548-559.
 % 
 % REQUIRED INPUT ARGUMENTS
 %
 %  sM (struct) SOM Toolbox struct or the codebook matrix of the map.
 %     (matrix) The matrix may be 3-dimensional in which case the first 
 %              two dimensions are taken for the map grid dimensions (msize).
 %
 % OPTIONAL INPUT ARGUMENTS
 %
 %  argID (string) Argument identifier string (see below).
 %  value (varies) Value for the argument (see below).
 %
 %  The optional arguments are given as 'argID',value -pairs. If the 
 %  value is unambiguous, it can be given without the preceeding argID.
 %  If an argument is given value multiple times, the last one is used. 
 %
 %  Below is the list of valid arguments: 
 %   'mask'      (vector) mask to be used in calculating
 %                        the interunit distances, size [dim  1]. Default is 
 %                        the one in sM (field sM.mask) or a vector of
 %                        ones if only a codebook matrix was given.
 %   'topol'     (struct) topology of the map. Default is the one
 %                        in sM (field sM.topol).
 %   'sTopol','som_topol' (struct) = 'topol'
 %   'msize'     (vector) map grid dimensions
 %   'lattice'   (string) map lattice 'rect' or 'hexa'
 %   'mode'      (string) 'min', 'mean', 'median' or 'max'
 %                        Map unit value computation method. In fact, 
 %                        eval-function is used to evaluate this, so 
 %                        you can give other computation methods as well.
 %                        Default is 'median'. 
 %
 % OUTPUT ARGUMENTS
 %
 %  U   (matrix) the unified distance matrix of the SOM 
 %               size 2*n1-1 x 2*n2-1, where n1 = msize(1) and n2 = msize(2)
 %
 % EXAMPLES
 %
 %  U = som_umat(sM);  
 %  U = som_umat(sM.codebook,sM.topol,'median','mask',[1 1 0 1]);
 %  U = som_umat(rand(10,10,4),'hexa','rect'); 
 % 
 % SEE ALSO
 %
 %  som_show    show the selected component planes and the u-matrix
 %  som_cplane  draw a 2D unified distance matrix
 
 % Copyright (c) 1997-2000 by the SOM toolbox programming team.
 % http://www.cis.hut.fi/projects/somtoolbox/
 
 % Version 1.0beta juuso 260997
 % Version 2.0beta juuso 151199, 151299, 200900
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% check arguments
 
 error(nargchk(1, Inf, nargin));  % check no. of input arguments is correct
 
 % sMap
 if isstruct(sMap), 
   M = sMap.codebook;
   sTopol = sMap.topol; 
   mask = sMap.mask;
 elseif isnumeric(sMap),
   M = sMap; 
   si = size(M);
   dim = si(end);
   if length(si)>2, msize = si(1:end-1);
   else msize = [si(1) 1];
   end
   munits = prod(msize);
   sTopol = som_set('som_topol','msize',msize,'lattice','rect','shape','sheet'); 
   mask = ones(dim,1);
   M = reshape(M,[munits,dim]);
 end
 mode = 'median';
 
 % varargin
 i=1; 
 while i<=length(varargin), 
   argok = 1; 
   if ischar(varargin{i}), 
     switch varargin{i}, 
       % argument IDs
      case 'mask',       i=i+1; mask = varargin{i}; 
      case 'msize',      i=i+1; sTopol.msize = varargin{i}; 
      case 'lattice',    i=i+1; sTopol.lattice = varargin{i};
      case {'topol','som_topol','sTopol'}, i=i+1; sTopol = varargin{i};
      case 'mode',       i=i+1; mode = varargin{i};
       % unambiguous values
      case {'hexa','rect'}, sTopol.lattice = varargin{i};
      case {'min','mean','median','max'}, mode = 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};
      case 'som_map',   sTopol = varargin{i}.topol;
      otherwise argok=0; 
     end
   else
     argok = 0; 
   end
   if ~argok, 
     disp(['(som_umat) Ignoring invalid argument #' num2str(i+1)]); 
   end
   i = i+1; 
 end
 
 % check
 [munits dim] = size(M);
 if prod(sTopol.msize)~=munits, 
   error('Map grid size does not match the number of map units.')
 end
 if length(sTopol.msize)>2, 
   error('Can only handle 1- and 2-dimensional map grids.')
 end
 if prod(sTopol.msize)==1,
   warning('Only one codebook vector.'); U = []; return;
 end
 if ~strcmp(sTopol.shape,'sheet'), 
   disp(['The ' sTopol.shape ' shape of the map ignored. Using sheet instead.']);
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% initialize variables
 
 y = sTopol.msize(1);
 x = sTopol.msize(2);
 lattice = sTopol.lattice;
 shape = sTopol.shape;
 M = reshape(M,[y x dim]);
 
 ux = 2 * x - 1; 
 uy = 2 * y - 1;
 U  = zeros(uy, ux);
 
 calc = sprintf('%s(a)',mode);
 
 if size(mask,2)>1, mask = mask'; end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% u-matrix computation
 
 % distances between map units
 
 if strcmp(lattice, 'rect'), % rectangular lattice
   
   for j=1:y, for i=1:x,
       if i<x, 
 	dx = (M(j,i,:) - M(j,i+1,:)).^2; % horizontal
 	U(2*j-1,2*i) = sqrt(mask'*dx(:));
       end 
       if j<y, 
 	dy = (M(j,i,:) - M(j+1,i,:)).^2; % vertical
 	U(2*j,2*i-1) = sqrt(mask'*dy(:));
       end
       if j<y & i<x,	
 	dz1 = (M(j,i,:) - M(j+1,i+1,:)).^2; % diagonals
 	dz2 = (M(j+1,i,:) - M(j,i+1,:)).^2;
 	U(2*j,2*i) = (sqrt(mask'*dz1(:))+sqrt(mask'*dz2(:)))/(2 * sqrt(2));
       end
     end
   end
 
 elseif strcmp(lattice, 'hexa') % hexagonal lattice
 
   for j=1:y, 
     for i=1:x,
       if i<x,
 	dx = (M(j,i,:) - M(j,i+1,:)).^2; % horizontal
 	U(2*j-1,2*i) = sqrt(mask'*dx(:));
       end
       
       if j<y, % diagonals
 	dy = (M(j,i,:) - M(j+1,i,:)).^2;
 	U(2*j,2*i-1) = sqrt(mask'*dy(:));	
 	
 	if rem(j,2)==0 & i<x,
 	  dz= (M(j,i,:) - M(j+1,i+1,:)).^2; 
 	  U(2*j,2*i) = sqrt(mask'*dz(:));
 	elseif rem(j,2)==1 & i>1,
 	  dz = (M(j,i,:) - M(j+1,i-1,:)).^2; 
 	  U(2*j,2*i-2) = sqrt(mask'*dz(:));
 	end
       end
     end
   end
   
 end
 
 % values on the units
 
 if (uy == 1 | ux == 1),
   % in 1-D case, mean is equal to median 
 
   ma = max([ux uy]);
   for i = 1:2:ma,
     if i>1 & i<ma, 
       a = [U(i-1) U(i+1)]; 
       U(i) = eval(calc);
     elseif i==1, U(i) = U(i+1); 
     else U(i) = U(i-1); % i==ma
     end
   end    
 
 elseif strcmp(lattice, 'rect')
 
   for j=1:2:uy, 
     for i=1:2:ux,
       if i>1 & j>1 & i<ux & j<uy,    % middle part of the map
 	a = [U(j,i-1) U(j,i+1) U(j-1,i) U(j+1,i)];        
       elseif j==1 & i>1 & i<ux,        % upper edge
 	a = [U(j,i-1) U(j,i+1) U(j+1,i)];
       elseif j==uy & i>1 & i<ux,       % lower edge
 	a = [U(j,i-1) U(j,i+1) U(j-1,i)];
       elseif i==1 & j>1 & j<uy,        % left edge
 	a = [U(j,i+1) U(j-1,i) U(j+1,i)];
       elseif i==ux & j>1 & j<uy,       % right edge
 	a = [U(j,i-1) U(j-1,i) U(j+1,i)];
       elseif i==1 & j==1,              % top left corner
 	a = [U(j,i+1) U(j+1,i)];
       elseif i==ux & j==1,             % top right corner
 	a = [U(j,i-1) U(j+1,i)];
       elseif i==1 & j==uy,             % bottom left corner
 	a = [U(j,i+1) U(j-1,i)];
       elseif i==ux & j==uy,            % bottom right corner
 	a = [U(j,i-1) U(j-1,i)];
       else
 	a = 0;
       end
       U(j,i) = eval(calc);
     end
   end
 
 elseif strcmp(lattice, 'hexa')
   
   for j=1:2:uy, 
     for i=1:2:ux,
       if i>1 & j>1 & i<ux & j<uy,      % middle part of the map
 	a = [U(j,i-1) U(j,i+1)];
 	if rem(j-1,4)==0, a = [a, U(j-1,i-1) U(j-1,i) U(j+1,i-1) U(j+1,i)];
 	else a = [a, U(j-1,i) U(j-1,i+1) U(j+1,i) U(j+1,i+1)]; end       
       elseif j==1 & i>1 & i<ux,        % upper edge
 	a = [U(j,i-1) U(j,i+1) U(j+1,i-1) U(j+1,i)];
       elseif j==uy & i>1 & i<ux,       % lower edge
 	a = [U(j,i-1) U(j,i+1)];
 	if rem(j-1,4)==0, a = [a, U(j-1,i-1) U(j-1,i)];
 	else a = [a, U(j-1,i) U(j-1,i+1)]; end
       elseif i==1 & j>1 & j<uy,        % left edge
 	a = U(j,i+1);
 	if rem(j-1,4)==0, a = [a, U(j-1,i) U(j+1,i)];
 	else a = [a, U(j-1,i) U(j-1,i+1) U(j+1,i) U(j+1,i+1)]; end
       elseif i==ux & j>1 & j<uy,       % right edge
 	a = U(j,i-1);
 	if rem(j-1,4)==0, a=[a, U(j-1,i) U(j-1,i-1) U(j+1,i) U(j+1,i-1)];
 	else a = [a, U(j-1,i) U(j+1,i)]; end
       elseif i==1 & j==1,              % top left corner
 	a = [U(j,i+1) U(j+1,i)];
       elseif i==ux & j==1,             % top right corner
 	a = [U(j,i-1) U(j+1,i-1) U(j+1,i)];
       elseif i==1 & j==uy,             % bottom left corner
 	if rem(j-1,4)==0, a = [U(j,i+1) U(j-1,i)];
 	else a = [U(j,i+1) U(j-1,i) U(j-1,i+1)]; end
       elseif i==ux & j==uy,            % bottom right corner
 	if rem(j-1,4)==0, a = [U(j,i-1) U(j-1,i) U(j-1,i-1)];
 	else a = [U(j,i-1) U(j-1,i)]; end
       else
 	a=0;
       end
       U(j,i) = eval(calc);
     end
   end
 end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% normalization between [0,1]
 
 % U = U - min(min(U)); 
 % ma = max(max(U)); if ma > 0, U = U / ma; end
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%