Mercurial > hg > camir-aes2014

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [Err,sPropTotal,sPropMunits,sPropComps] = som_distortion3(sM,D,rad)
wolffd@0 2
wolffd@0 3 %SOM_DISTORTION3 Map distortion measures.
wolffd@0 4 %
wolffd@0 5 % [sE,Err] = som_distortion3(sM,[D],[rad]);
wolffd@0 6 %
wolffd@0 7 % sE = som_distortion3(sM);
wolffd@0 8 %
wolffd@0 9 % Input and output arguments ([]'s are optional):
wolffd@0 10 % sM (struct) map struct
wolffd@0 11 % [D] (matrix) a matrix, size dlen x dim
wolffd@0 12 % (struct) data or map struct
wolffd@0 13 % by default the map struct is used
wolffd@0 14 % [rad] (scalar) neighborhood radius, looked from sM.trainhist
wolffd@0 15 % by default, or = 1 if that has no valid values
wolffd@0 16 %
wolffd@0 17 % Err (matrix) size munits x dim x 3
wolffd@0 18 % distortion error elements (quantization error,
wolffd@0 19 % neighborhood bias, and neighborhood variance)
wolffd@0 20 % for each map unit and component
wolffd@0 21 % sPropTotal (struct) .n = length of data
wolffd@0 22 % .h = mean neighborhood function value
wolffd@0 23 % .err = errors
wolffd@0 24 % sPropMunits (struct) .Ni = hits per map unit
wolffd@0 25 % .Hi = sum of neighborhood values for each map unit
wolffd@0 26 % .Err = errors per map unit
wolffd@0 27 % sPropComps (struct) .e1 = total squared distance to centroid
wolffd@0 28 % .eq = total squared distance to BMU
wolffd@0 29 % .Err = errors per component
wolffd@0 30 %
wolffd@0 31 % See also SOM_QUALITY.
wolffd@0 32
wolffd@0 33 % Contributed to SOM Toolbox 2.0, January 3rd, 2002 by Juha Vesanto
wolffd@0 34 % Copyright (c) by Juha Vesanto
wolffd@0 35 % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0 36
wolffd@0 37 % Version 2.0beta juuso 030102
wolffd@0 38
wolffd@0 39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 40 %% arguments
wolffd@0 41
wolffd@0 42 % map
wolffd@0 43 [munits dim] = size(sM.codebook);
wolffd@0 44
wolffd@0 45 % neighborhood radius
wolffd@0 46 if nargin<3,
wolffd@0 47 if ~isempty(sM.trainhist),
wolffd@0 48 rad = sM.trainhist(end).radius_fin;
wolffd@0 49 else
wolffd@0 50 rad = 1;
wolffd@0 51 end
wolffd@0 52 end
wolffd@0 53 if rad<eps, rad = eps; end
wolffd@0 54 if isempty(rad) | isnan(rad), rad = 1; end
wolffd@0 55
wolffd@0 56 % neighborhood function
wolffd@0 57 Ud = som_unit_dists(sM.topol);
wolffd@0 58 switch sM.neigh,
wolffd@0 59 case 'bubble', H = (Ud <= rad);
wolffd@0 60 case 'gaussian', H = exp(-(Ud.^2)/(2*rad*rad));
wolffd@0 61 case 'cutgauss', H = exp(-(Ud.^2)/(2*rad*rad)) .* (Ud <= rad);
wolffd@0 62 case 'ep', H = (1 - (Ud.^2)/rad) .* (Ud <= rad);
wolffd@0 63 end
wolffd@0 64 Hi = sum(H,2);
wolffd@0 65
wolffd@0 66 % data
wolffd@0 67 if nargin<2, D = sM.codebook; end
wolffd@0 68 if isstruct(D), D = D.data; end
wolffd@0 69 [dlen dim] = size(D);
wolffd@0 70
wolffd@0 71 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 72 %% quality measures
wolffd@0 73
wolffd@0 74 % find Voronoi sets, and calculate their properties
wolffd@0 75
wolffd@0 76 [bmus,qerr] = som_bmus(sM,D);
wolffd@0 77 M = sM.codebook;
wolffd@0 78 Vn = M;
wolffd@0 79 Vm = M;
wolffd@0 80 Ni = zeros(munits,dim);
wolffd@0 81 for i=1:munits,
wolffd@0 82 inds = find(bmus==i);
wolffd@0 83 Ni(i,:) = sum(isfinite(D(inds,:)),1); % size of Voronoi set
wolffd@0 84 if any(Ni(i,:)), Vn(i,:) = centroid(D(inds,:),M(i,:)); end % centroid of Voronoi set
wolffd@0 85 Vm(i,:) = centroid(M,M(i,:),H(i,:)'); % centroid of neighborhood
wolffd@0 86 end
wolffd@0 87
wolffd@0 88 HN = repmat(Hi,1,dim).*Ni;
wolffd@0 89
wolffd@0 90 %% distortion
wolffd@0 91
wolffd@0 92 % quantization error (in each Voronoi set and for each component)
wolffd@0 93
wolffd@0 94 Eqx = zeros(munits,dim);
wolffd@0 95 Dx = (Vn(bmus,:) - D).^2;
wolffd@0 96 Dx(isnan(Dx)) = 0;
wolffd@0 97 for i = 1:dim,
wolffd@0 98 Eqx(:,i) = full(sum(sparse(bmus,1:dlen,Dx(:,i),munits,dlen),2));
wolffd@0 99 end
wolffd@0 100 Eqx = repmat(Hi,1,dim).*Eqx;
wolffd@0 101
wolffd@0 102 % bias in neighborhood (in each Voronoi set / component)
wolffd@0 103
wolffd@0 104 Enb = (Vn-Vm).^2;
wolffd@0 105 Enb = HN.*Enb;
wolffd@0 106
wolffd@0 107 % variance in neighborhood (in each Voronoi set / component)
wolffd@0 108
wolffd@0 109 Env = zeros(munits,dim);
wolffd@0 110 for i=1:munits, Env(i,:) = H(i,:)*(M-Vm(i*ones(munits,1),:)).^2; end
wolffd@0 111 Env = Ni.*Env;
wolffd@0 112
wolffd@0 113 % total distortion (in each Voronoi set / component)
wolffd@0 114
wolffd@0 115 Ed = Eqx + Enb + Env;
wolffd@0 116
wolffd@0 117 %% other error measures
wolffd@0 118
wolffd@0 119 % squared quantization error (to data centroid)
wolffd@0 120
wolffd@0 121 me = centroid(D,mean(M));
wolffd@0 122 Dx = D - me(ones(dlen,1),:);
wolffd@0 123 Dx(isnan(Dx)) = 0;
wolffd@0 124 e1 = sum(Dx.^2,1);
wolffd@0 125
wolffd@0 126 % squared quantization error (to map units)
wolffd@0 127
wolffd@0 128 Dx = D - M(bmus,:);
wolffd@0 129 Dx(isnan(Dx)) = 0;
wolffd@0 130 eq = sum(Dx.^2,1);
wolffd@0 131
wolffd@0 132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 133 %% output
wolffd@0 134
wolffd@0 135 % distortion error matrix
wolffd@0 136
wolffd@0 137 Err = zeros(munits,dim,5);
wolffd@0 138 Err(:,:,1) = Eqx;
wolffd@0 139 Err(:,:,2) = Enb;
wolffd@0 140 Err(:,:,3) = Env;
wolffd@0 141
wolffd@0 142 % total errors
wolffd@0 143
wolffd@0 144 sPropTotal = struct('n',sum(Ni),'h',mean(Hi),'e1',sum(e1),'err',sum(sum(Err,2),1));
wolffd@0 145
wolffd@0 146 % properties of map units
wolffd@0 147
wolffd@0 148 sPropMunits = struct('Ni',[],'Hi',[],'Err',[]);
wolffd@0 149 sPropMunits.Ni = Ni;
wolffd@0 150 sPropMunits.Hi = Hi;
wolffd@0 151 sPropMunits.Err = squeeze(sum(Err,2));
wolffd@0 152
wolffd@0 153 % properties of components
wolffd@0 154
wolffd@0 155 sPropComps = struct('Err',[],'e1',[],'eq',[]);
wolffd@0 156 sPropComps.Err = squeeze(sum(Err,1));
wolffd@0 157 sPropComps.e1 = e1;
wolffd@0 158 sPropComps.eq = eq;
wolffd@0 159
wolffd@0 160
wolffd@0 161 return;
wolffd@0 162
wolffd@0 163
wolffd@0 164 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
wolffd@0 165 %% subfunctions
wolffd@0 166
wolffd@0 167 function v = centroid(D,default,weights)
wolffd@0 168
wolffd@0 169 [n dim] = size(D);
wolffd@0 170 I = sparse(isnan(D));
wolffd@0 171 D(I) = 0;
wolffd@0 172
wolffd@0 173 if nargin==3,
wolffd@0 174 W = weights(:,ones(1,dim));
wolffd@0 175 W(I) = 0;
wolffd@0 176 D = D.*W;
wolffd@0 177 nn = sum(W,1);
wolffd@0 178 else
wolffd@0 179 nn = n-sum(I,1);
wolffd@0 180 end
wolffd@0 181
wolffd@0 182 c = sum(D,1);
wolffd@0 183 v = default;
wolffd@0 184 i = find(nn>0);
wolffd@0 185 v(i) = c(i)./nn(i);
wolffd@0 186
wolffd@0 187 return;
wolffd@0 188
wolffd@0 189
wolffd@0 190 function vis
wolffd@0 191
wolffd@0 192 figure
wolffd@0 193 som_show(sM,'color',{Hi,'Hi'},'color',{Ni,'hits'},...
wolffd@0 194 'color',{Ed,'distortion'},'color',{Eqx,'qxerror'},...
wolffd@0 195 'color',{Enb,'N-bias'},'color',{Env,'N-Var'});
wolffd@0 196
wolffd@0 197 ed = Eqx + Enb + Env;
wolffd@0 198 i = find(ed>0);
wolffd@0 199 eqx = 0*ed; eqx(i) = Eqx(i)./ed(i);
wolffd@0 200 enb = 0*ed; enb(i) = Enb(i)./ed(i);
wolffd@0 201 env = 0*ed; env(i) = Env(i)./ed(i);
wolffd@0 202
wolffd@0 203 figure
wolffd@0 204 som_show(sM,'color',Hi,'color',Ni,'color',Ed,...
wolffd@0 205 'color',eqx,'color',enb,'color',env);
wolffd@0 206
wolffd@0 207