function [adm,dm] = som_distortion(sM, D)

%SOM_DISTORTION Calculate distortion measure for the map.
%
% [adm,dm] = som_distortion(sMap, D)
%
%  adm = som_distortion(sMap,D);
%  [adm,dm] = som_distortion(sMap,D);
%  som_show(sM,'color',dm);
%
%  Input and output arguments: 
%   sMap     (struct) a map struct
%   D        (struct) a data struct
%            (matrix) size dlen x dim, a data matrix
%
%   adm      (scalar) average distortion measure
%   dm       (vector) size munits x 1, total distortion for each unit
%
% The distortion measure is defined as: 
%                                           2
%    E = sum sum h(i,bmu(j)) ||m(i) - x(j)|| 
%         i   j    
% 
% where m(i) is the ith prototype vector of SOM, x(j) is the jth data
% vector, and h(.,.) is the neighborhood function. In case of fixed
% neighborhood and discreet data, the distortion measure can be
% interpreted as the energy function of the SOM. Note, though, that
% the learning rule that follows from the distortion measure is
% different from the SOM training rule, so SOM only minimizes the
% distortion measure approximately.
% 
% NOTE: when calculating BMUs and distances, the mask of the given 
%       map is used.
%
% See also SOM_QUALITY, SOM_BMUS, SOM_HITS.

% Reference: Kohonen, T., "Self-Organizing Map", 2nd ed., 
%    Springer-Verlag, Berlin, 1995, pp. 120-121.

% Contributed to SOM Toolbox vs2, Feb 3rd, 2000 by Juha Vesanto
% http://www.cis.hut.fi/projects/somtoolbox/

% Version 2.0beta juuso 030200

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% check arguments

% input arguments
if nargin < 2, error('Not enough input arguments.'); end

% map
M = sM.codebook;
munits = prod(sM.topol.msize);

% data
if isstruct(D), D = D.data; end
[dlen dim] = size(D);

% neighborhood radius
rad = NaN;
if ~isempty(sM.trainhist), 
  rad = sM.trainhist(end).radius_fin;
end
if isnan(rad) | isempty(rad), rad = 1; end
if rad<eps, rad = eps; end

% neighborhood  
Ud = som_unit_dists(sM.topol); 
switch sM.neigh, 
 case 'bubble',   H = (Ud <= rad);
 case 'gaussian', H = exp(-(Ud.^2)/(2*rad*rad)); 
 case 'cutgauss', H = exp(-(Ud.^2)/(2*rad*rad)) .* (Ud <= rad);
 case 'ep',       H = (1 - (Ud.^2)/rad) .* (Ud <= rad);
end  

% total distortion measure
mu_x_1 = ones(munits,1);
dm = zeros(munits,1);
ns = zeros(munits,1);
for i=1:dlen,
  x = D(i,:);                        % data sample
  known = ~isnan(x);                 % its known components
  Dx = M(:,known) - x(mu_x_1,known); % each map unit minus the vector
  dist2 = (Dx.^2)*sM.mask(known);    % squared distances  
  [qerr bmu] = min(dist2);           % find BMU
  dm = dm + dist2.*H(:,bmu);         % add to distortion measure
  ns = ns + H(:,bmu);
end 
  
% average distortion measure
adm = sum(dm)/dlen;

return;

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


