comparison toolboxes/MIRtoolbox1.3.2/somtoolbox/som_distortion3.m @ 0:e9a9cd732c1e tip

first hg version after svn

+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 radius
+if nargin<3, 
+  if ~isempty(sM.trainhist), 
+    rad = sM.trainhist(end).radius_fin; 
+  else 
+    rad = 1; 
+  end
+end
+if rad<eps, rad = eps; end
+if isempty(rad) | isnan(rad), rad = 1; end
+
+% neighborhood function
+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  
+Hi = sum(H,2); 
+
+% data
+if nargin<2, D = sM.codebook; end
+if 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 neighborhood
+end
+
+HN = 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)); 
+end
+Eqx           = 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; end
+Env = 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 matrix
+
+Err  = zeros(munits,dim,5); 
+Err(:,:,1) = Eqx; 
+Err(:,:,2) = Enb; 
+Err(:,:,3) = Env; 
+
+% total errors
+
+sPropTotal = struct('n',sum(Ni),'h',mean(Hi),'e1',sum(e1),'err',sum(sum(Err,2),1));
+
+% properties of map units
+
+sPropMunits = struct('Ni',[],'Hi',[],'Err',[]); 
+sPropMunits.Ni  = Ni; 
+sPropMunits.Hi  = Hi; 
+sPropMunits.Err = squeeze(sum(Err,2));
+
+% properties of components
+
+sPropComps = struct('Err',[],'e1',[],'eq',[]);
+sPropComps.Err = squeeze(sum(Err,1));
+sPropComps.e1  = e1; 
+sPropComps.eq  = eq;
+
+
+return; 
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
+%% subfunctions
+
+function 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); 
+
+