function [Err,sPropTotal,sPropMunits,sPropComps] = som_distortion3(sM,D,rad)%SOM_DISTORTION3 Map distortion measures.%  % [sE,Err] = som_distortion3(sM,[D],[rad]);% %  sE = som_distortion3(sM); %%  Input and output arguments ([]'s are optional): %   sM          (struct) map struct%   [D]         (matrix) a matrix, size dlen x dim%               (struct) data or map struct%                        by default the map struct is used%   [rad]       (scalar) neighborhood radius, looked from sM.trainhist%                        by default, or = 1 if that has no valid values%                           %   Err         (matrix) size munits x dim x 3%                        distortion error elements (quantization error, %                        neighborhood bias, and neighborhood variance)%                        for each map unit and component%   sPropTotal  (struct) .n   = length of data%                        .h   = mean neighborhood function value%                        .err = errors%   sPropMunits (struct) .Ni  = hits per map unit%                        .Hi  = sum of neighborhood values for each map unit%                        .Err = errors per map unit%   sPropComps  (struct) .e1  = total squared distance to centroid%                        .eq  = total squared distance to BMU%                        .Err = errors per component%% See also  SOM_QUALITY.% Contributed to SOM Toolbox 2.0, January 3rd, 2002 by Juha Vesanto% Copyright (c) by Juha Vesanto% http://www.cis.hut.fi/projects/somtoolbox/% Version 2.0beta juuso 030102%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% arguments% map[munits dim] = size(sM.codebook);% neighborhood radiusif nargin<3,   if ~isempty(sM.trainhist),     rad = sM.trainhist(end).radius_fin;   else     rad = 1;   endendif rad<eps, rad = eps; endif isempty(rad) | isnan(rad), rad = 1; end% neighborhood functionUd = 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  Hi = sum(H,2); % dataif nargin<2, D = sM.codebook; endif isstruct(D), D = D.data; end[dlen dim] = size(D);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% quality measures% find Voronoi sets, and calculate their properties[bmus,qerr] = som_bmus(sM,D); M  = sM.codebook; Vn = M; Vm = M; Ni = zeros(munits,dim);for i=1:munits,   inds    = find(bmus==i);     Ni(i,:) = sum(isfinite(D(inds,:)),1);                      % size of Voronoi set  if any(Ni(i,:)), Vn(i,:) = centroid(D(inds,:),M(i,:)); end % centroid of Voronoi set    Vm(i,:) = centroid(M,M(i,:),H(i,:)');                      % centroid of neighborhoodendHN = repmat(Hi,1,dim).*Ni; %% distortion% quantization error (in each Voronoi set and for each component)Eqx           = zeros(munits,dim); Dx            = (Vn(bmus,:) - D).^2; Dx(isnan(Dx)) = 0; for i = 1:dim,   Eqx(:,i)    = full(sum(sparse(bmus,1:dlen,Dx(:,i),munits,dlen),2)); endEqx           = repmat(Hi,1,dim).*Eqx;   % bias in neighborhood (in each Voronoi set / component)Enb = (Vn-Vm).^2;Enb = HN.*Enb; % variance in neighborhood (in each Voronoi set / component)Env = zeros(munits,dim);for i=1:munits, Env(i,:) = H(i,:)*(M-Vm(i*ones(munits,1),:)).^2; endEnv = Ni.*Env; % total distortion (in each Voronoi set / component)Ed = Eqx + Enb + Env;%% other error measures% squared quantization error (to data centroid)me            = centroid(D,mean(M));Dx            = D - me(ones(dlen,1),:); Dx(isnan(Dx)) = 0; e1            = sum(Dx.^2,1); % squared quantization error (to map units)Dx            = D - M(bmus,:);Dx(isnan(Dx)) = 0; eq            = sum(Dx.^2,1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output% distortion error matrixErr  = zeros(munits,dim,5); Err(:,:,1) = Eqx; Err(:,:,2) = Enb; Err(:,:,3) = Env; % total errorssPropTotal = struct('n',sum(Ni),'h',mean(Hi),'e1',sum(e1),'err',sum(sum(Err,2),1));% properties of map unitssPropMunits = struct('Ni',[],'Hi',[],'Err',[]); sPropMunits.Ni  = Ni; sPropMunits.Hi  = Hi; sPropMunits.Err = squeeze(sum(Err,2));% properties of componentssPropComps = struct('Err',[],'e1',[],'eq',[]);sPropComps.Err = squeeze(sum(Err,1));sPropComps.e1  = e1; sPropComps.eq  = eq;return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55%% subfunctionsfunction v = centroid(D,default,weights)    [n dim] = size(D);  I       = sparse(isnan(D));  D(I)    = 0;    if nargin==3,     W    = weights(:,ones(1,dim));     W(I) = 0;     D    = D.*W;    nn   = sum(W,1);  else    nn   = n-sum(I,1);  end   c    = sum(D,1);  v    = default;   i    = find(nn>0);   v(i) = c(i)./nn(i);        return; function vis  figure  som_show(sM,'color',{Hi,'Hi'},'color',{Ni,'hits'},...           'color',{Ed,'distortion'},'color',{Eqx,'qxerror'},...           'color',{Enb,'N-bias'},'color',{Env,'N-Var'});  ed = Eqx + Enb + Env;  i = find(ed>0);   eqx = 0*ed; eqx(i) = Eqx(i)./ed(i);  enb = 0*ed; enb(i) = Enb(i)./ed(i);  env = 0*ed; env(i) = Env(i)./ed(i);  figure  som_show(sM,'color',Hi,'color',Ni,'color',Ed,...           'color',eqx,'color',enb,'color',env);