Mercurial > hg > camir-aes2014

annotate toolboxes/MIRtoolbox1.3.2/somtoolbox/som_batchtrain.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 [sMap,sTrain] = som_batchtrain(sMap, D, varargin)
wolffd@0 2
wolffd@0 3 %SOM_BATCHTRAIN Use batch algorithm to train the Self-Organizing Map.
wolffd@0 4 %
wolffd@0 5 % [sM,sT] = som_batchtrain(sM, D, [argID, value, ...])
wolffd@0 6 %
wolffd@0 7 % sM = som_batchtrain(sM,D);
wolffd@0 8 % sM = som_batchtrain(sM,sD,'radius',[10 3 2 1 0.1],'tracking',3);
wolffd@0 9 % [M,sT] = som_batchtrain(M,D,'ep','msize',[10 3],'hexa');
wolffd@0 10 %
wolffd@0 11 % Input and output arguments ([]'s are optional):
wolffd@0 12 % sM (struct) map struct, the trained and updated map is returned
wolffd@0 13 % (matrix) codebook matrix of a self-organizing map
wolffd@0 14 % size munits x dim or msize(1) x ... x msize(k) x dim
wolffd@0 15 % The trained map codebook is returned.
wolffd@0 16 % D (struct) training data; data struct
wolffd@0 17 % (matrix) training data, size dlen x dim
wolffd@0 18 % [argID, (string) See below. The values which are unambiguous can
wolffd@0 19 % value] (varies) be given without the preceeding argID.
wolffd@0 20 %
wolffd@0 21 % sT (struct) learning parameters used during the training
wolffd@0 22 %
wolffd@0 23 % Here are the valid argument IDs and corresponding values. The values which
wolffd@0 24 % are unambiguous (marked with '*') can be given without the preceeding argID.
wolffd@0 25 % 'mask' (vector) BMU search mask, size dim x 1
wolffd@0 26 % 'msize' (vector) map size
wolffd@0 27 % 'radius' (vector) neighborhood radiuses, length 1, 2 or trainlen
wolffd@0 28 % 'radius_ini' (scalar) initial training radius
wolffd@0 29 % 'radius_fin' (scalar) final training radius
wolffd@0 30 % 'tracking' (scalar) tracking level, 0-3
wolffd@0 31 % 'trainlen' (scalar) training length in epochs
wolffd@0 32 % 'train' *(struct) train struct, parameters for training
wolffd@0 33 % 'sTrain','som_train' = 'train'
wolffd@0 34 % 'neigh' *(string) neighborhood function, 'gaussian', 'cutgauss',
wolffd@0 35 % 'ep' or 'bubble'
wolffd@0 36 % 'topol' *(struct) topology struct
wolffd@0 37 % 'som_topol','sTopol' = 'topol'
wolffd@0 38 % 'lattice' *(string) map lattice, 'hexa' or 'rect'
wolffd@0 39 % 'shape' *(string) map shape, 'sheet', 'cyl' or 'toroid'
wolffd@0 40 % 'weights' (vector) sample weights: each sample is weighted
wolffd@0 41 %
wolffd@0 42 % For more help, try 'type som_batchtrain' or check out online documentation.
wolffd@0 43 % See also SOM_MAKE, SOM_SEQTRAIN, SOM_TRAIN_STRUCT.
wolffd@0 44
wolffd@0 45 %%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 46 %
wolffd@0 47 % som_batchtrain
wolffd@0 48 %
wolffd@0 49 % PURPOSE
wolffd@0 50 %
wolffd@0 51 % Trains a Self-Organizing Map using the batch algorithm.
wolffd@0 52 %
wolffd@0 53 % SYNTAX
wolffd@0 54 %
wolffd@0 55 % sM = som_batchtrain(sM,D);
wolffd@0 56 % sM = som_batchtrain(sM,sD);
wolffd@0 57 % sM = som_batchtrain(...,'argID',value,...);
wolffd@0 58 % sM = som_batchtrain(...,value,...);
wolffd@0 59 % [sM,sT] = som_batchtrain(M,D,...);
wolffd@0 60 %
wolffd@0 61 % DESCRIPTION
wolffd@0 62 %
wolffd@0 63 % Trains the given SOM (sM or M above) with the given training data
wolffd@0 64 % (sD or D) using batch training algorithm. If no optional arguments
wolffd@0 65 % (argID, value) are given, a default training is done. Using optional
wolffd@0 66 % arguments the training parameters can be specified. Returns the
wolffd@0 67 % trained and updated SOM and a train struct which contains
wolffd@0 68 % information on the training.
wolffd@0 69 %
wolffd@0 70 % REFERENCES
wolffd@0 71 %
wolffd@0 72 % Kohonen, T., "Self-Organizing Map", 2nd ed., Springer-Verlag,
wolffd@0 73 % Berlin, 1995, pp. 127-128.
wolffd@0 74 % Kohonen, T., "Things you haven't heard about the Self-Organizing
wolffd@0 75 % Map", In proceedings of International Conference
wolffd@0 76 % on Neural Networks (ICNN), San Francisco, 1993, pp. 1147-1156.
wolffd@0 77 %
wolffd@0 78 % KNOWN BUGS
wolffd@0 79 %
wolffd@0 80 % Batchtrain does not work correctly for a map with a single unit.
wolffd@0 81 % This is because of the way 'min'-function works.
wolffd@0 82 %
wolffd@0 83 % REQUIRED INPUT ARGUMENTS
wolffd@0 84 %
wolffd@0 85 % sM The map to be trained.
wolffd@0 86 % (struct) map struct
wolffd@0 87 % (matrix) codebook matrix (field .data of map struct)
wolffd@0 88 % Size is either [munits dim], in which case the map grid
wolffd@0 89 % dimensions (msize) should be specified with optional arguments,
wolffd@0 90 % or [msize(1) ... msize(k) dim] in which case the map
wolffd@0 91 % grid dimensions are taken from the size of the matrix.
wolffd@0 92 % Lattice, by default, is 'rect' and shape 'sheet'.
wolffd@0 93 % D Training data.
wolffd@0 94 % (struct) data struct
wolffd@0 95 % (matrix) data matrix, size [dlen dim]
wolffd@0 96 %
wolffd@0 97 % OPTIONAL INPUT ARGUMENTS
wolffd@0 98 %
wolffd@0 99 % argID (string) Argument identifier string (see below).
wolffd@0 100 % value (varies) Value for the argument (see below).
wolffd@0 101 %
wolffd@0 102 % The optional arguments can be given as 'argID',value -pairs. If an
wolffd@0 103 % argument is given value multiple times, the last one is
wolffd@0 104 % used. The valid IDs and corresponding values are listed below. The values
wolffd@0 105 % which are unambiguous (marked with '*') can be given without the
wolffd@0 106 % preceeding argID.
wolffd@0 107 %
wolffd@0 108 % Below is the list of valid arguments:
wolffd@0 109 % 'mask' (vector) BMU search mask, size dim x 1. Default is
wolffd@0 110 % the one in sM (field '.mask') or a vector of
wolffd@0 111 % ones if only a codebook matrix was given.
wolffd@0 112 % 'msize' (vector) map grid dimensions. Default is the one
wolffd@0 113 % in sM (field sM.topol.msize) or
wolffd@0 114 % 'si = size(sM); msize = si(1:end-1);'
wolffd@0 115 % if only a codebook matrix was given.
wolffd@0 116 % 'radius' (vector) neighborhood radius
wolffd@0 117 % length = 1: radius_ini = radius
wolffd@0 118 % length = 2: [radius_ini radius_fin] = radius
wolffd@0 119 % length > 2: the vector given neighborhood
wolffd@0 120 % radius for each step separately
wolffd@0 121 % trainlen = length(radius)
wolffd@0 122 % 'radius_ini' (scalar) initial training radius
wolffd@0 123 % 'radius_fin' (scalar) final training radius
wolffd@0 124 % 'tracking' (scalar) tracking level: 0, 1 (default), 2 or 3
wolffd@0 125 % 0 - estimate time
wolffd@0 126 % 1 - track time and quantization error
wolffd@0 127 % 2 - plot quantization error
wolffd@0 128 % 3 - plot quantization error and two first
wolffd@0 129 % components
wolffd@0 130 % 'trainlen' (scalar) training length in epochs
wolffd@0 131 % 'train' *(struct) train struct, parameters for training.
wolffd@0 132 % Default parameters, unless specified,
wolffd@0 133 % are acquired using SOM_TRAIN_STRUCT (this
wolffd@0 134 % also applies for 'trainlen', 'radius_ini'
wolffd@0 135 % and 'radius_fin').
wolffd@0 136 % 'sTrain', 'som_topol' (struct) = 'train'
wolffd@0 137 % 'neigh' *(string) The used neighborhood function. Default is
wolffd@0 138 % the one in sM (field '.neigh') or 'gaussian'
wolffd@0 139 % if only a codebook matrix was given. Other
wolffd@0 140 % possible values is 'cutgauss', 'ep' and 'bubble'.
wolffd@0 141 % 'topol' *(struct) topology of the map. Default is the one
wolffd@0 142 % in sM (field '.topol').
wolffd@0 143 % 'sTopol', 'som_topol' (struct) = 'topol'
wolffd@0 144 % 'lattice' *(string) map lattice. Default is the one in sM
wolffd@0 145 % (field sM.topol.lattice) or 'rect'
wolffd@0 146 % if only a codebook matrix was given.
wolffd@0 147 % 'shape' *(string) map shape. Default is the one in sM
wolffd@0 148 % (field sM.topol.shape) or 'sheet'
wolffd@0 149 % if only a codebook matrix was given.
wolffd@0 150 % 'weights' (vector) weight for each data vector: during training,
wolffd@0 151 % each data sample is weighted with the corresponding
wolffd@0 152 % value, for example giving weights = [1 1 2 1]
wolffd@0 153 % would have the same result as having third sample
wolffd@0 154 % appear 2 times in the data
wolffd@0 155 %
wolffd@0 156 % OUTPUT ARGUMENTS
wolffd@0 157 %
wolffd@0 158 % sM the trained map
wolffd@0 159 % (struct) if a map struct was given as input argument, a
wolffd@0 160 % map struct is also returned. The current training
wolffd@0 161 % is added to the training history (sM.trainhist).
wolffd@0 162 % The 'neigh' and 'mask' fields of the map struct
wolffd@0 163 % are updated to match those of the training.
wolffd@0 164 % (matrix) if a matrix was given as input argument, a matrix
wolffd@0 165 % is also returned with the same size as the input
wolffd@0 166 % argument.
wolffd@0 167 % sT (struct) train struct; information of the accomplished training
wolffd@0 168 %
wolffd@0 169 % EXAMPLES
wolffd@0 170 %
wolffd@0 171 % Simplest case:
wolffd@0 172 % sM = som_batchtrain(sM,D);
wolffd@0 173 % sM = som_batchtrain(sM,sD);
wolffd@0 174 %
wolffd@0 175 % To change the tracking level, 'tracking' argument is specified:
wolffd@0 176 % sM = som_batchtrain(sM,D,'tracking',3);
wolffd@0 177 %
wolffd@0 178 % The change training parameters, the optional arguments 'train','neigh',
wolffd@0 179 % 'mask','trainlen','radius','radius_ini' and 'radius_fin' are used.
wolffd@0 180 % sM = som_batchtrain(sM,D,'neigh','cutgauss','trainlen',10,'radius_fin',0);
wolffd@0 181 %
wolffd@0 182 % Another way to specify training parameters is to create a train struct:
wolffd@0 183 % sTrain = som_train_struct(sM,'dlen',size(D,1));
wolffd@0 184 % sTrain = som_set(sTrain,'neigh','cutgauss');
wolffd@0 185 % sM = som_batchtrain(sM,D,sTrain);
wolffd@0 186 %
wolffd@0 187 % By default the neighborhood radius goes linearly from radius_ini to
wolffd@0 188 % radius_fin. If you want to change this, you can use the 'radius' argument
wolffd@0 189 % to specify the neighborhood radius for each step separately:
wolffd@0 190 % sM = som_batchtrain(sM,D,'radius',[5 3 1 1 1 1 0.5 0.5 0.5]);
wolffd@0 191 %
wolffd@0 192 % You don't necessarily have to use the map struct, but you can operate
wolffd@0 193 % directly with codebook matrices. However, in this case you have to
wolffd@0 194 % specify the topology of the map in the optional arguments. The
wolffd@0 195 % following commads are identical (M is originally a 200 x dim sized matrix):
wolffd@0 196 % M = som_batchtrain(M,D,'msize',[20 10],'lattice','hexa','shape','cyl');
wolffd@0 197 % or
wolffd@0 198 % M = som_batchtrain(M,D,'msize',[20 10],'hexa','cyl');
wolffd@0 199 % or
wolffd@0 200 % sT= som_set('som_topol','msize',[20 10],'lattice','hexa','shape','cyl');
wolffd@0 201 % M = som_batchtrain(M,D,sT);
wolffd@0 202 % or
wolffd@0 203 % M = reshape(M,[20 10 dim]);
wolffd@0 204 % M = som_batchtrain(M,D,'hexa','cyl');
wolffd@0 205 %
wolffd@0 206 % The som_batchtrain also returns a train struct with information on the
wolffd@0 207 % accomplished training. This struct is also added to the end of the
wolffd@0 208 % trainhist field of map struct, in case a map struct was given.
wolffd@0 209 % [M,sTrain] = som_batchtrain(M,D,'msize',[20 10]);
wolffd@0 210 % [sM,sTrain] = som_batchtrain(sM,D); % sM.trainhist{end}==sTrain
wolffd@0 211 %
wolffd@0 212 % SEE ALSO
wolffd@0 213 %
wolffd@0 214 % som_make Initialize and train a SOM using default parameters.
wolffd@0 215 % som_seqtrain Train SOM with sequential algorithm.
wolffd@0 216 % som_train_struct Determine default training parameters.
wolffd@0 217
wolffd@0 218 % Copyright (c) 1997-2000 by the SOM toolbox programming team.
wolffd@0 219 % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0 220
wolffd@0 221 % Version 1.0beta juuso 071197 041297
wolffd@0 222 % Version 2.0beta juuso 101199
wolffd@0 223
wolffd@0 224 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 225 %% Check arguments
wolffd@0 226
wolffd@0 227 error(nargchk(2, Inf, nargin)); % check the number of input arguments
wolffd@0 228
wolffd@0 229 % map
wolffd@0 230 struct_mode = isstruct(sMap);
wolffd@0 231 if struct_mode,
wolffd@0 232 sTopol = sMap.topol;
wolffd@0 233 else
wolffd@0 234 orig_size = size(sMap);
wolffd@0 235 if ndims(sMap) > 2,
wolffd@0 236 si = size(sMap); dim = si(end); msize = si(1:end-1);
wolffd@0 237 M = reshape(sMap,[prod(msize) dim]);
wolffd@0 238 else
wolffd@0 239 msize = [orig_size(1) 1];
wolffd@0 240 dim = orig_size(2);
wolffd@0 241 end
wolffd@0 242 sMap = som_map_struct(dim,'msize',msize);
wolffd@0 243 sTopol = sMap.topol;
wolffd@0 244 end
wolffd@0 245 [munits dim] = size(sMap.codebook);
wolffd@0 246
wolffd@0 247 % data
wolffd@0 248 if isstruct(D),
wolffd@0 249 data_name = D.name;
wolffd@0 250 D = D.data;
wolffd@0 251 else
wolffd@0 252 data_name = inputname(2);
wolffd@0 253 end
wolffd@0 254 nonempty = find(sum(isnan(D),2) < dim);
wolffd@0 255 D = D(nonempty,:); % remove empty vectors from the data
wolffd@0 256 [dlen ddim] = size(D); % check input dimension
wolffd@0 257 if dim ~= ddim,
wolffd@0 258 error('Map and data input space dimensions disagree.');
wolffd@0 259 end
wolffd@0 260
wolffd@0 261 % varargin
wolffd@0 262 sTrain = som_set('som_train','algorithm','batch','neigh', ...
wolffd@0 263 sMap.neigh,'mask',sMap.mask,'data_name',data_name);
wolffd@0 264 radius = [];
wolffd@0 265 tracking = 1;
wolffd@0 266 weights = 1;
wolffd@0 267
wolffd@0 268 i=1;
wolffd@0 269 while i<=length(varargin),
wolffd@0 270 argok = 1;
wolffd@0 271 if ischar(varargin{i}),
wolffd@0 272 switch varargin{i},
wolffd@0 273 % argument IDs
wolffd@0 274 case 'msize', i=i+1; sTopol.msize = varargin{i};
wolffd@0 275 case 'lattice', i=i+1; sTopol.lattice = varargin{i};
wolffd@0 276 case 'shape', i=i+1; sTopol.shape = varargin{i};
wolffd@0 277 case 'mask', i=i+1; sTrain.mask = varargin{i};
wolffd@0 278 case 'neigh', i=i+1; sTrain.neigh = varargin{i};
wolffd@0 279 case 'trainlen', i=i+1; sTrain.trainlen = varargin{i};
wolffd@0 280 case 'tracking', i=i+1; tracking = varargin{i};
wolffd@0 281 case 'weights', i=i+1; weights = varargin{i};
wolffd@0 282 case 'radius_ini', i=i+1; sTrain.radius_ini = varargin{i};
wolffd@0 283 case 'radius_fin', i=i+1; sTrain.radius_fin = varargin{i};
wolffd@0 284 case 'radius',
wolffd@0 285 i=i+1;
wolffd@0 286 l = length(varargin{i});
wolffd@0 287 if l==1,
wolffd@0 288 sTrain.radius_ini = varargin{i};
wolffd@0 289 else
wolffd@0 290 sTrain.radius_ini = varargin{i}(1);
wolffd@0 291 sTrain.radius_fin = varargin{i}(end);
wolffd@0 292 if l>2, radius = varargin{i}; end
wolffd@0 293 end
wolffd@0 294 case {'sTrain','train','som_train'}, i=i+1; sTrain = varargin{i};
wolffd@0 295 case {'topol','sTopol','som_topol'},
wolffd@0 296 i=i+1;
wolffd@0 297 sTopol = varargin{i};
wolffd@0 298 if prod(sTopol.msize) ~= munits,
wolffd@0 299 error('Given map grid size does not match the codebook size.');
wolffd@0 300 end
wolffd@0 301 % unambiguous values
wolffd@0 302 case {'hexa','rect'}, sTopol.lattice = varargin{i};
wolffd@0 303 case {'sheet','cyl','toroid'}, sTopol.shape = varargin{i};
wolffd@0 304 case {'gaussian','cutgauss','ep','bubble'}, sTrain.neigh = varargin{i};
wolffd@0 305 otherwise argok=0;
wolffd@0 306 end
wolffd@0 307 elseif isstruct(varargin{i}) & isfield(varargin{i},'type'),
wolffd@0 308 switch varargin{i}(1).type,
wolffd@0 309 case 'som_topol',
wolffd@0 310 sTopol = varargin{i};
wolffd@0 311 if prod(sTopol.msize) ~= munits,
wolffd@0 312 error('Given map grid size does not match the codebook size.');
wolffd@0 313 end
wolffd@0 314 case 'som_train', sTrain = varargin{i};
wolffd@0 315 otherwise argok=0;
wolffd@0 316 end
wolffd@0 317 else
wolffd@0 318 argok = 0;
wolffd@0 319 end
wolffd@0 320 if ~argok,
wolffd@0 321 disp(['(som_batchtrain) Ignoring invalid argument #' num2str(i+2)]);
wolffd@0 322 end
wolffd@0 323 i = i+1;
wolffd@0 324 end
wolffd@0 325
wolffd@0 326 % take only weights of non-empty vectors
wolffd@0 327 if length(weights)>dlen, weights = weights(nonempty); end
wolffd@0 328
wolffd@0 329 % trainlen
wolffd@0 330 if ~isempty(radius), sTrain.trainlen = length(radius); end
wolffd@0 331
wolffd@0 332 % check topology
wolffd@0 333 if struct_mode,
wolffd@0 334 if ~strcmp(sTopol.lattice,sMap.topol.lattice) | ...
wolffd@0 335 ~strcmp(sTopol.shape,sMap.topol.shape) | ...
wolffd@0 336 any(sTopol.msize ~= sMap.topol.msize),
wolffd@0 337 warning('Changing the original map topology.');
wolffd@0 338 end
wolffd@0 339 end
wolffd@0 340 sMap.topol = sTopol;
wolffd@0 341
wolffd@0 342 % complement the training struct
wolffd@0 343 sTrain = som_train_struct(sTrain,sMap,'dlen',dlen);
wolffd@0 344 if isempty(sTrain.mask), sTrain.mask = ones(dim,1); end
wolffd@0 345
wolffd@0 346 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 347 %% initialize
wolffd@0 348
wolffd@0 349 M = sMap.codebook;
wolffd@0 350 mask = sTrain.mask;
wolffd@0 351 trainlen = sTrain.trainlen;
wolffd@0 352
wolffd@0 353 % neighborhood radius
wolffd@0 354 if trainlen==1,
wolffd@0 355 radius = sTrain.radius_ini;
wolffd@0 356 elseif length(radius)<=2,
wolffd@0 357 r0 = sTrain.radius_ini; r1 = sTrain.radius_fin;
wolffd@0 358 radius = r1 + fliplr((0:(trainlen-1))/(trainlen-1)) * (r0 - r1);
wolffd@0 359 else
wolffd@0 360 % nil
wolffd@0 361 end
wolffd@0 362
wolffd@0 363 % distance between map units in the output space
wolffd@0 364 % Since in the case of gaussian and ep neighborhood functions, the
wolffd@0 365 % equations utilize squares of the unit distances and in bubble case
wolffd@0 366 % it doesn't matter which is used, the unitdistances and neighborhood
wolffd@0 367 % radiuses are squared.
wolffd@0 368 Ud = som_unit_dists(sTopol);
wolffd@0 369 Ud = Ud.^2;
wolffd@0 370 radius = radius.^2;
wolffd@0 371 % zero neighborhood radius may cause div-by-zero error
wolffd@0 372 radius(find(radius==0)) = eps;
wolffd@0 373
wolffd@0 374 % The training algorithm involves calculating weighted Euclidian distances
wolffd@0 375 % to all map units for each data vector. Basically this is done as
wolffd@0 376 % for i=1:dlen,
wolffd@0 377 % for j=1:munits,
wolffd@0 378 % for k=1:dim
wolffd@0 379 % Dist(j,i) = Dist(j,i) + mask(k) * (D(i,k) - M(j,k))^2;
wolffd@0 380 % end
wolffd@0 381 % end
wolffd@0 382 % end
wolffd@0 383 % where mask is the weighting vector for distance calculation. However, taking
wolffd@0 384 % into account that distance between vectors m and v can be expressed as
wolffd@0 385 % |m - v|^2 = sum_i ((m_i - v_i)^2) = sum_i (m_i^2 + v_i^2 - 2*m_i*v_i)
wolffd@0 386 % this can be made much faster by transforming it to a matrix operation:
wolffd@0 387 % Dist = (M.^2)*mask*ones(1,d) + ones(m,1)*mask'*(D'.^2) - 2*M*diag(mask)*D'
wolffd@0 388 % Of the involved matrices, several are constant, as the mask and data do
wolffd@0 389 % not change during training. Therefore they are calculated beforehand.
wolffd@0 390
wolffd@0 391 % For the case where there are unknown components in the data, each data
wolffd@0 392 % vector will have an individual mask vector so that for that unit, the
wolffd@0 393 % unknown components are not taken into account in distance calculation.
wolffd@0 394 % In addition all NaN's are changed to zeros so that they don't screw up
wolffd@0 395 % the matrix multiplications and behave correctly in updating step.
wolffd@0 396 Known = ~isnan(D);
wolffd@0 397 W1 = (mask*ones(1,dlen)) .* Known';
wolffd@0 398 D(find(~Known)) = 0;
wolffd@0 399
wolffd@0 400 % constant matrices
wolffd@0 401 WD = 2*diag(mask)*D'; % constant matrix
wolffd@0 402 dconst = ((D.^2)*mask)'; % constant in distance calculation for each data sample
wolffd@0 403 % W2 = ones(munits,1)*mask'; D2 = (D'.^2);
wolffd@0 404
wolffd@0 405 % initialize tracking
wolffd@0 406 start = clock;
wolffd@0 407 qe = zeros(trainlen,1);
wolffd@0 408
wolffd@0 409 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 410 %% Action
wolffd@0 411
wolffd@0 412 % With the 'blen' parameter you can control the memory consumption
wolffd@0 413 % of the algorithm, which is in practive directly proportional
wolffd@0 414 % to munits*blen. If you're having problems with memory, try to
wolffd@0 415 % set the value of blen lower.
wolffd@0 416 blen = min(munits,dlen);
wolffd@0 417
wolffd@0 418 % reserve some space
wolffd@0 419 bmus = zeros(1,dlen);
wolffd@0 420 ddists = zeros(1,dlen);
wolffd@0 421
wolffd@0 422 for t = 1:trainlen,
wolffd@0 423
wolffd@0 424 % batchy train - this is done a block of data (inds) at a time
wolffd@0 425 % rather than in a single sweep to save memory consumption.
wolffd@0 426 % The 'Dist' and 'Hw' matrices have size munits*blen
wolffd@0 427 % which - if you have a lot of data - would be HUGE if you
wolffd@0 428 % calculated it all at once. A single-sweep version would
wolffd@0 429 % look like this:
wolffd@0 430 % Dist = (M.^2)*W1 - M*WD; %+ W2*D2
wolffd@0 431 % [ddists, bmus] = min(Dist);
wolffd@0 432 % (notice that the W2*D2 term can be ignored since it is constant)
wolffd@0 433 % This "batchy" version is the same as single-sweep if blen=dlen.
wolffd@0 434 i0 = 0;
wolffd@0 435 while i0+1<=dlen,
wolffd@0 436 inds = [(i0+1):min(dlen,i0+blen)]; i0 = i0+blen;
wolffd@0 437 Dist = (M.^2)*W1(:,inds) - M*WD(:,inds);
wolffd@0 438 [ddists(inds), bmus(inds)] = min(Dist);
wolffd@0 439 end
wolffd@0 440
wolffd@0 441 % tracking
wolffd@0 442 if tracking > 0,
wolffd@0 443 ddists = ddists+dconst; % add the constant term
wolffd@0 444 ddists(ddists<0) = 0; % rounding errors...
wolffd@0 445 qe(t) = mean(sqrt(ddists));
wolffd@0 446 trackplot(M,D,tracking,start,t,qe);
wolffd@0 447 end
wolffd@0 448
wolffd@0 449 % neighborhood
wolffd@0 450 % notice that the elements Ud and radius have been squared!
wolffd@0 451 % note: 'bubble' matches the original "Batch Map" algorithm
wolffd@0 452 switch sTrain.neigh,
wolffd@0 453 case 'bubble', H = (Ud<=radius(t));
wolffd@0 454 case 'gaussian', H = exp(-Ud/(2*radius(t)));
wolffd@0 455 case 'cutgauss', H = exp(-Ud/(2*radius(t))) .* (Ud<=radius(t));
wolffd@0 456 case 'ep', H = (1-Ud/radius(t)) .* (Ud<=radius(t));
wolffd@0 457 end
wolffd@0 458
wolffd@0 459 % update
wolffd@0 460
wolffd@0 461 % In principle the updating step goes like this: replace each map unit
wolffd@0 462 % by the average of the data vectors that were in its neighborhood.
wolffd@0 463 % The contribution, or activation, of data vectors in the mean can
wolffd@0 464 % be varied with the neighborhood function. This activation is given
wolffd@0 465 % by matrix H. So, for each map unit the new weight vector is
wolffd@0 466 %
wolffd@0 467 % m = sum_i (h_i * d_i) / sum_i (h_i),
wolffd@0 468 %
wolffd@0 469 % where i denotes the index of data vector. Since the values of
wolffd@0 470 % neighborhood function h_i are the same for all data vectors belonging to
wolffd@0 471 % the Voronoi set of the same map unit, the calculation is actually done
wolffd@0 472 % by first calculating a partition matrix P with elements p_ij=1 if the
wolffd@0 473 % BMU of data vector j is i.
wolffd@0 474
wolffd@0 475 P = sparse(bmus,[1:dlen],weights,munits,dlen);
wolffd@0 476
wolffd@0 477 % Then the sum of vectors in each Voronoi set are calculated (P*D) and the
wolffd@0 478 % neighborhood is taken into account by calculating a weighted sum of the
wolffd@0 479 % Voronoi sum (H*). The "activation" matrix A is the denominator of the
wolffd@0 480 % equation above.
wolffd@0 481
wolffd@0 482 S = H*(P*D);
wolffd@0 483 A = H*(P*Known);
wolffd@0 484
wolffd@0 485 % If you'd rather make this without using the Voronoi sets try the following:
wolffd@0 486 % Hi = H(:,bmus);
wolffd@0 487 % S = Hi * D; % "sum_i (h_i * d_i)"
wolffd@0 488 % A = Hi * Known; % "sum_i (h_i)"
wolffd@0 489 % The bad news is that the matrix Hi has size [munits x dlen]...
wolffd@0 490
wolffd@0 491 % only update units for which the "activation" is nonzero
wolffd@0 492 nonzero = find(A > 0);
wolffd@0 493 M(nonzero) = S(nonzero) ./ A(nonzero);
wolffd@0 494
wolffd@0 495 end; % for t = 1:trainlen
wolffd@0 496
wolffd@0 497 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 498 %% Build / clean up the return arguments
wolffd@0 499
wolffd@0 500 % tracking
wolffd@0 501 if tracking > 0, fprintf(1,'\n'); end
wolffd@0 502
wolffd@0 503 % update structures
wolffd@0 504 sTrain = som_set(sTrain,'time',datestr(now,0));
wolffd@0 505 if struct_mode,
wolffd@0 506 sMap = som_set(sMap,'codebook',M,'mask',sTrain.mask,'neigh',sTrain.neigh);
wolffd@0 507 tl = length(sMap.trainhist);
wolffd@0 508 sMap.trainhist(tl+1) = sTrain;
wolffd@0 509 else
wolffd@0 510 sMap = reshape(M,orig_size);
wolffd@0 511 end
wolffd@0 512
wolffd@0 513 return;
wolffd@0 514
wolffd@0 515
wolffd@0 516 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0 517 %% subfunctions
wolffd@0 518
wolffd@0 519 %%%%%%%%
wolffd@0 520 function [] = trackplot(M,D,tracking,start,n,qe)
wolffd@0 521
wolffd@0 522 l = length(qe);
wolffd@0 523 elap_t = etime(clock,start);
wolffd@0 524 tot_t = elap_t*l/n;
wolffd@0 525 fprintf(1,'\rTraining: %3.0f/ %3.0f s',elap_t,tot_t)
wolffd@0 526 switch tracking
wolffd@0 527 case 1,
wolffd@0 528 case 2,
wolffd@0 529 plot(1:n,qe(1:n),(n+1):l,qe((n+1):l))
wolffd@0 530 title('Quantization error after each epoch');
wolffd@0 531 drawnow
wolffd@0 532 otherwise,
wolffd@0 533 subplot(2,1,1), plot(1:n,qe(1:n),(n+1):l,qe((n+1):l))
wolffd@0 534 title('Quantization error after each epoch');
wolffd@0 535 subplot(2,1,2), plot(M(:,1),M(:,2),'ro',D(:,1),D(:,2),'b+');
wolffd@0 536 title('First two components of map units (o) and data vectors (+)');
wolffd@0 537 drawnow
wolffd@0 538 end
wolffd@0 539 % end of trackplot