wolffd@0: function [adm,admu,tdmu] = som_distortion(sM, D, arg1, arg2)
wolffd@0: 
wolffd@0: %SOM_DISTORTION Calculate distortion measure for the map.
wolffd@0: %
wolffd@0: % [adm,admu,tdmu] = som_distortion(sMap, D, [radius], ['prob'])
wolffd@0: %
wolffd@0: %  adm = som_distortion(sMap,D);
wolffd@0: %  [adm,admu] = som_distortion(sMap,D);
wolffd@0: %  som_show(sMap,'color',admu);
wolffd@0: %
wolffd@0: %  Input and output arguments: 
wolffd@0: %   sMap     (struct) a map struct
wolffd@0: %   D        (struct) a data struct
wolffd@0: %            (matrix) size dlen x dim, a data matrix
wolffd@0: %   [radius] (scalar) neighborhood function radius to be used.
wolffd@0: %                     Defaults to the last radius_fin in the 
wolffd@0: %                     trainhist field of the map struct, or 1 if
wolffd@0: %                     that is missing.
wolffd@0: %   ['prob'] (string) If given, this argument forces the 
wolffd@0: %                     neigborhood function values for each map
wolffd@0: %                     unit to be normalized so that they sum to 1.
wolffd@0: %
wolffd@0: %   adm      (scalar) average distortion measure (sum(dm)/dlen)
wolffd@0: %   admu     (vector) size munits x 1, average distortion in each unit 
wolffd@0: %   tdmu     (vector) size munits x 1, total distortion for each unit
wolffd@0: %
wolffd@0: % The distortion measure is defined as: 
wolffd@0: %                                           2
wolffd@0: %    E = sum sum h(bmu(i),j) ||m(j) - x(i)|| 
wolffd@0: %         i   j    
wolffd@0: % 
wolffd@0: % where m(i) is the ith prototype vector of SOM, x(j) is the jth data
wolffd@0: % vector, and h(.,.) is the neighborhood function. In case of fixed
wolffd@0: % neighborhood and discreet data, the distortion measure can be
wolffd@0: % interpreted as the energy function of the SOM. Note, though, that
wolffd@0: % the learning rule that follows from the distortion measure is
wolffd@0: % different from the SOM training rule, so SOM only minimizes the
wolffd@0: % distortion measure approximately.
wolffd@0: % 
wolffd@0: % If the 'prob' argument is given, the distortion measure can be 
wolffd@0: % interpreted as an expected quantization error when the neighborhood 
wolffd@0: % function values give the likelyhoods of accidentally assigning 
wolffd@0: % vector j to unit i. The normal quantization error is a special case 
wolffd@0: % of this with zero incorrect assignement likelihood. 
wolffd@0: % 
wolffd@0: % NOTE: when calculating BMUs and distances, the mask of the given 
wolffd@0: %       map is used.
wolffd@0: %
wolffd@0: % See also SOM_QUALITY, SOM_BMUS, SOM_HITS.
wolffd@0: 
wolffd@0: % Reference: Kohonen, T., "Self-Organizing Map", 2nd ed., 
wolffd@0: %    Springer-Verlag, Berlin, 1995, pp. 120-121.
wolffd@0: %
wolffd@0: %    Graepel, T., Burger, M. and Obermayer, K., 
wolffd@0: %    "Phase Transitions in Stochastic Self-Organizing Maps",
wolffd@0: %    Physical Review E, Vol 56, No 4, pp. 3876-3890 (1997).
wolffd@0: 
wolffd@0: % Contributed to SOM Toolbox vs2, Feb 3rd, 2000 by Juha Vesanto
wolffd@0: % Copyright (c) by Juha Vesanto
wolffd@0: % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0: 
wolffd@0: % Version 2.0beta juuso 030200
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% check arguments
wolffd@0: 
wolffd@0: % input arguments
wolffd@0: if nargin < 2, error('Not enough input arguments.'); end
wolffd@0: 
wolffd@0: % map
wolffd@0: M = sM.codebook;
wolffd@0: munits = prod(sM.topol.msize);
wolffd@0: 
wolffd@0: % data
wolffd@0: if isstruct(D), D = D.data; end
wolffd@0: [dlen dim] = size(D);
wolffd@0: 
wolffd@0: % arg1, arg2
wolffd@0: rad = NaN;
wolffd@0: normalize = 0;
wolffd@0: if nargin>2, 
wolffd@0:   if isnumeric(arg1), rad = arg1;
wolffd@0:   elseif ischar(arg1) & strcmp(arg1,'prob'), normalize = 0;
wolffd@0:   end
wolffd@0: end
wolffd@0: if nargin>3, 
wolffd@0:   if isnumeric(arg2), rad = arg2;
wolffd@0:   elseif ischar(arg2) & strcmp(arg2,'prob'), normalize = 0;
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % neighborhood radius
wolffd@0: if isempty(rad) | isnan(rad), 
wolffd@0:   if ~isempty(sM.trainhist), rad = sM.trainhist(end).radius_fin;
wolffd@0:   else rad = 1; 
wolffd@0:   end
wolffd@0: end
wolffd@0: if rad<eps, rad = eps; end
wolffd@0: 
wolffd@0: % neighborhood  
wolffd@0: Ud = som_unit_dists(sM.topol); 
wolffd@0: switch sM.neigh, 
wolffd@0:  case 'bubble',   H = (Ud <= rad);
wolffd@0:  case 'gaussian', H = exp(-(Ud.^2)/(2*rad*rad)); 
wolffd@0:  case 'cutgauss', H = exp(-(Ud.^2)/(2*rad*rad)) .* (Ud <= rad);
wolffd@0:  case 'ep',       H = (1 - (Ud.^2)/rad) .* (Ud <= rad);
wolffd@0: end  
wolffd@0: if normalize, 
wolffd@0:   for i=1:munits, H(:,i) = H(:,i)/sum(H(:,i)); end
wolffd@0: end
wolffd@0: 
wolffd@0: % total distortion measure
wolffd@0: mu_x_1 = ones(munits,1);
wolffd@0: tdmu = zeros(munits,1);
wolffd@0: hits = zeros(munits,1);
wolffd@0: for i=1:dlen,
wolffd@0:   x = D(i,:);                        % data sample
wolffd@0:   known = ~isnan(x);                 % its known components
wolffd@0:   Dx = M(:,known) - x(mu_x_1,known); % each map unit minus the vector
wolffd@0:   dist2 = (Dx.^2)*sM.mask(known);    % squared distances  
wolffd@0:   [qerr bmu] = min(dist2);           % find BMU
wolffd@0:   tdmu = tdmu + dist2.*H(:,bmu);     % add to distortion measure
wolffd@0:   hits(bmu) = hits(bmu)+1;           % add to hits
wolffd@0: end 
wolffd@0: 
wolffd@0: % average distortion per unit
wolffd@0: admu = tdmu; 
wolffd@0: ind = find(hits>0);
wolffd@0: admu(ind) = admu(ind) ./ hits(ind);
wolffd@0:   
wolffd@0: % average distortion measure
wolffd@0: adm = sum(tdmu)/dlen;
wolffd@0: 
wolffd@0: return;
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: 
wolffd@0: