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