Daniel@0: function [sMap,sTrain] = som_batchtrain(sMap, D, varargin) Daniel@0: Daniel@0: %SOM_BATCHTRAIN Use batch algorithm to train the Self-Organizing Map. Daniel@0: % Daniel@0: % [sM,sT] = som_batchtrain(sM, D, [argID, value, ...]) Daniel@0: % Daniel@0: % sM = som_batchtrain(sM,D); Daniel@0: % sM = som_batchtrain(sM,sD,'radius',[10 3 2 1 0.1],'tracking',3); Daniel@0: % [M,sT] = som_batchtrain(M,D,'ep','msize',[10 3],'hexa'); Daniel@0: % Daniel@0: % Input and output arguments ([]'s are optional): Daniel@0: % sM (struct) map struct, the trained and updated map is returned Daniel@0: % (matrix) codebook matrix of a self-organizing map Daniel@0: % size munits x dim or msize(1) x ... x msize(k) x dim Daniel@0: % The trained map codebook is returned. Daniel@0: % D (struct) training data; data struct Daniel@0: % (matrix) training data, size dlen x dim Daniel@0: % [argID, (string) See below. The values which are unambiguous can Daniel@0: % value] (varies) be given without the preceeding argID. Daniel@0: % Daniel@0: % sT (struct) learning parameters used during the training Daniel@0: % Daniel@0: % Here are the valid argument IDs and corresponding values. The values which Daniel@0: % are unambiguous (marked with '*') can be given without the preceeding argID. Daniel@0: % 'mask' (vector) BMU search mask, size dim x 1 Daniel@0: % 'msize' (vector) map size Daniel@0: % 'radius' (vector) neighborhood radiuses, length 1, 2 or trainlen Daniel@0: % 'radius_ini' (scalar) initial training radius Daniel@0: % 'radius_fin' (scalar) final training radius Daniel@0: % 'tracking' (scalar) tracking level, 0-3 Daniel@0: % 'trainlen' (scalar) training length in epochs Daniel@0: % 'train' *(struct) train struct, parameters for training Daniel@0: % 'sTrain','som_train' = 'train' Daniel@0: % 'neigh' *(string) neighborhood function, 'gaussian', 'cutgauss', Daniel@0: % 'ep' or 'bubble' Daniel@0: % 'topol' *(struct) topology struct Daniel@0: % 'som_topol','sTopol' = 'topol' Daniel@0: % 'lattice' *(string) map lattice, 'hexa' or 'rect' Daniel@0: % 'shape' *(string) map shape, 'sheet', 'cyl' or 'toroid' Daniel@0: % 'weights' (vector) sample weights: each sample is weighted Daniel@0: % Daniel@0: % For more help, try 'type som_batchtrain' or check out online documentation. Daniel@0: % See also SOM_MAKE, SOM_SEQTRAIN, SOM_TRAIN_STRUCT. Daniel@0: Daniel@0: %%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Daniel@0: % Daniel@0: % som_batchtrain Daniel@0: % Daniel@0: % PURPOSE Daniel@0: % Daniel@0: % Trains a Self-Organizing Map using the batch algorithm. Daniel@0: % Daniel@0: % SYNTAX Daniel@0: % Daniel@0: % sM = som_batchtrain(sM,D); Daniel@0: % sM = som_batchtrain(sM,sD); Daniel@0: % sM = som_batchtrain(...,'argID',value,...); Daniel@0: % sM = som_batchtrain(...,value,...); Daniel@0: % [sM,sT] = som_batchtrain(M,D,...); Daniel@0: % Daniel@0: % DESCRIPTION Daniel@0: % Daniel@0: % Trains the given SOM (sM or M above) with the given training data Daniel@0: % (sD or D) using batch training algorithm. If no optional arguments Daniel@0: % (argID, value) are given, a default training is done. Using optional Daniel@0: % arguments the training parameters can be specified. Returns the Daniel@0: % trained and updated SOM and a train struct which contains Daniel@0: % information on the training. Daniel@0: % Daniel@0: % REFERENCES Daniel@0: % Daniel@0: % Kohonen, T., "Self-Organizing Map", 2nd ed., Springer-Verlag, Daniel@0: % Berlin, 1995, pp. 127-128. Daniel@0: % Kohonen, T., "Things you haven't heard about the Self-Organizing Daniel@0: % Map", In proceedings of International Conference Daniel@0: % on Neural Networks (ICNN), San Francisco, 1993, pp. 1147-1156. Daniel@0: % Daniel@0: % KNOWN BUGS Daniel@0: % Daniel@0: % Batchtrain does not work correctly for a map with a single unit. Daniel@0: % This is because of the way 'min'-function works. Daniel@0: % Daniel@0: % REQUIRED INPUT ARGUMENTS Daniel@0: % Daniel@0: % sM The map to be trained. Daniel@0: % (struct) map struct Daniel@0: % (matrix) codebook matrix (field .data of map struct) Daniel@0: % Size is either [munits dim], in which case the map grid Daniel@0: % dimensions (msize) should be specified with optional arguments, Daniel@0: % or [msize(1) ... msize(k) dim] in which case the map Daniel@0: % grid dimensions are taken from the size of the matrix. Daniel@0: % Lattice, by default, is 'rect' and shape 'sheet'. Daniel@0: % D Training data. Daniel@0: % (struct) data struct Daniel@0: % (matrix) data matrix, size [dlen dim] Daniel@0: % Daniel@0: % OPTIONAL INPUT ARGUMENTS Daniel@0: % Daniel@0: % argID (string) Argument identifier string (see below). Daniel@0: % value (varies) Value for the argument (see below). Daniel@0: % Daniel@0: % The optional arguments can be given as 'argID',value -pairs. If an Daniel@0: % argument is given value multiple times, the last one is Daniel@0: % used. The valid IDs and corresponding values are listed below. The values Daniel@0: % which are unambiguous (marked with '*') can be given without the Daniel@0: % preceeding argID. Daniel@0: % Daniel@0: % Below is the list of valid arguments: Daniel@0: % 'mask' (vector) BMU search mask, size dim x 1. Default is Daniel@0: % the one in sM (field '.mask') or a vector of Daniel@0: % ones if only a codebook matrix was given. Daniel@0: % 'msize' (vector) map grid dimensions. Default is the one Daniel@0: % in sM (field sM.topol.msize) or Daniel@0: % 'si = size(sM); msize = si(1:end-1);' Daniel@0: % if only a codebook matrix was given. Daniel@0: % 'radius' (vector) neighborhood radius Daniel@0: % length = 1: radius_ini = radius Daniel@0: % length = 2: [radius_ini radius_fin] = radius Daniel@0: % length > 2: the vector given neighborhood Daniel@0: % radius for each step separately Daniel@0: % trainlen = length(radius) Daniel@0: % 'radius_ini' (scalar) initial training radius Daniel@0: % 'radius_fin' (scalar) final training radius Daniel@0: % 'tracking' (scalar) tracking level: 0, 1 (default), 2 or 3 Daniel@0: % 0 - estimate time Daniel@0: % 1 - track time and quantization error Daniel@0: % 2 - plot quantization error Daniel@0: % 3 - plot quantization error and two first Daniel@0: % components Daniel@0: % 'trainlen' (scalar) training length in epochs Daniel@0: % 'train' *(struct) train struct, parameters for training. Daniel@0: % Default parameters, unless specified, Daniel@0: % are acquired using SOM_TRAIN_STRUCT (this Daniel@0: % also applies for 'trainlen', 'radius_ini' Daniel@0: % and 'radius_fin'). Daniel@0: % 'sTrain', 'som_topol' (struct) = 'train' Daniel@0: % 'neigh' *(string) The used neighborhood function. Default is Daniel@0: % the one in sM (field '.neigh') or 'gaussian' Daniel@0: % if only a codebook matrix was given. Other Daniel@0: % possible values is 'cutgauss', 'ep' and 'bubble'. Daniel@0: % 'topol' *(struct) topology of the map. Default is the one Daniel@0: % in sM (field '.topol'). Daniel@0: % 'sTopol', 'som_topol' (struct) = 'topol' Daniel@0: % 'lattice' *(string) map lattice. Default is the one in sM Daniel@0: % (field sM.topol.lattice) or 'rect' Daniel@0: % if only a codebook matrix was given. Daniel@0: % 'shape' *(string) map shape. Default is the one in sM Daniel@0: % (field sM.topol.shape) or 'sheet' Daniel@0: % if only a codebook matrix was given. Daniel@0: % 'weights' (vector) weight for each data vector: during training, Daniel@0: % each data sample is weighted with the corresponding Daniel@0: % value, for example giving weights = [1 1 2 1] Daniel@0: % would have the same result as having third sample Daniel@0: % appear 2 times in the data Daniel@0: % Daniel@0: % OUTPUT ARGUMENTS Daniel@0: % Daniel@0: % sM the trained map Daniel@0: % (struct) if a map struct was given as input argument, a Daniel@0: % map struct is also returned. The current training Daniel@0: % is added to the training history (sM.trainhist). Daniel@0: % The 'neigh' and 'mask' fields of the map struct Daniel@0: % are updated to match those of the training. Daniel@0: % (matrix) if a matrix was given as input argument, a matrix Daniel@0: % is also returned with the same size as the input Daniel@0: % argument. Daniel@0: % sT (struct) train struct; information of the accomplished training Daniel@0: % Daniel@0: % EXAMPLES Daniel@0: % Daniel@0: % Simplest case: Daniel@0: % sM = som_batchtrain(sM,D); Daniel@0: % sM = som_batchtrain(sM,sD); Daniel@0: % Daniel@0: % To change the tracking level, 'tracking' argument is specified: Daniel@0: % sM = som_batchtrain(sM,D,'tracking',3); Daniel@0: % Daniel@0: % The change training parameters, the optional arguments 'train','neigh', Daniel@0: % 'mask','trainlen','radius','radius_ini' and 'radius_fin' are used. Daniel@0: % sM = som_batchtrain(sM,D,'neigh','cutgauss','trainlen',10,'radius_fin',0); Daniel@0: % Daniel@0: % Another way to specify training parameters is to create a train struct: Daniel@0: % sTrain = som_train_struct(sM,'dlen',size(D,1)); Daniel@0: % sTrain = som_set(sTrain,'neigh','cutgauss'); Daniel@0: % sM = som_batchtrain(sM,D,sTrain); Daniel@0: % Daniel@0: % By default the neighborhood radius goes linearly from radius_ini to Daniel@0: % radius_fin. If you want to change this, you can use the 'radius' argument Daniel@0: % to specify the neighborhood radius for each step separately: Daniel@0: % sM = som_batchtrain(sM,D,'radius',[5 3 1 1 1 1 0.5 0.5 0.5]); Daniel@0: % Daniel@0: % You don't necessarily have to use the map struct, but you can operate Daniel@0: % directly with codebook matrices. However, in this case you have to Daniel@0: % specify the topology of the map in the optional arguments. The Daniel@0: % following commads are identical (M is originally a 200 x dim sized matrix): Daniel@0: % M = som_batchtrain(M,D,'msize',[20 10],'lattice','hexa','shape','cyl'); Daniel@0: % or Daniel@0: % M = som_batchtrain(M,D,'msize',[20 10],'hexa','cyl'); Daniel@0: % or Daniel@0: % sT= som_set('som_topol','msize',[20 10],'lattice','hexa','shape','cyl'); Daniel@0: % M = som_batchtrain(M,D,sT); Daniel@0: % or Daniel@0: % M = reshape(M,[20 10 dim]); Daniel@0: % M = som_batchtrain(M,D,'hexa','cyl'); Daniel@0: % Daniel@0: % The som_batchtrain also returns a train struct with information on the Daniel@0: % accomplished training. This struct is also added to the end of the Daniel@0: % trainhist field of map struct, in case a map struct was given. Daniel@0: % [M,sTrain] = som_batchtrain(M,D,'msize',[20 10]); Daniel@0: % [sM,sTrain] = som_batchtrain(sM,D); % sM.trainhist{end}==sTrain Daniel@0: % Daniel@0: % SEE ALSO Daniel@0: % Daniel@0: % som_make Initialize and train a SOM using default parameters. Daniel@0: % som_seqtrain Train SOM with sequential algorithm. Daniel@0: % som_train_struct Determine default training parameters. Daniel@0: Daniel@0: % Copyright (c) 1997-2000 by the SOM toolbox programming team. Daniel@0: % http://www.cis.hut.fi/projects/somtoolbox/ Daniel@0: Daniel@0: % Version 1.0beta juuso 071197 041297 Daniel@0: % Version 2.0beta juuso 101199 Daniel@0: Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Daniel@0: %% Check arguments Daniel@0: Daniel@0: error(nargchk(2, Inf, nargin)); % check the number of input arguments Daniel@0: Daniel@0: % map Daniel@0: struct_mode = isstruct(sMap); Daniel@0: if struct_mode, Daniel@0: sTopol = sMap.topol; Daniel@0: else Daniel@0: orig_size = size(sMap); Daniel@0: if ndims(sMap) > 2, Daniel@0: si = size(sMap); dim = si(end); msize = si(1:end-1); Daniel@0: M = reshape(sMap,[prod(msize) dim]); Daniel@0: else Daniel@0: msize = [orig_size(1) 1]; Daniel@0: dim = orig_size(2); Daniel@0: end Daniel@0: sMap = som_map_struct(dim,'msize',msize); Daniel@0: sTopol = sMap.topol; Daniel@0: end Daniel@0: [munits dim] = size(sMap.codebook); Daniel@0: Daniel@0: % data Daniel@0: if isstruct(D), Daniel@0: data_name = D.name; Daniel@0: D = D.data; Daniel@0: else Daniel@0: data_name = inputname(2); Daniel@0: end Daniel@0: nonempty = find(sum(isnan(D),2) < dim); Daniel@0: D = D(nonempty,:); % remove empty vectors from the data Daniel@0: [dlen ddim] = size(D); % check input dimension Daniel@0: if dim ~= ddim, Daniel@0: error('Map and data input space dimensions disagree.'); Daniel@0: end Daniel@0: Daniel@0: % varargin Daniel@0: sTrain = som_set('som_train','algorithm','batch','neigh', ... Daniel@0: sMap.neigh,'mask',sMap.mask,'data_name',data_name); Daniel@0: radius = []; Daniel@0: tracking = 1; Daniel@0: weights = 1; Daniel@0: Daniel@0: i=1; Daniel@0: while i<=length(varargin), Daniel@0: argok = 1; Daniel@0: if ischar(varargin{i}), Daniel@0: switch varargin{i}, Daniel@0: % argument IDs Daniel@0: case 'msize', i=i+1; sTopol.msize = varargin{i}; Daniel@0: case 'lattice', i=i+1; sTopol.lattice = varargin{i}; Daniel@0: case 'shape', i=i+1; sTopol.shape = varargin{i}; Daniel@0: case 'mask', i=i+1; sTrain.mask = varargin{i}; Daniel@0: case 'neigh', i=i+1; sTrain.neigh = varargin{i}; Daniel@0: case 'trainlen', i=i+1; sTrain.trainlen = varargin{i}; Daniel@0: case 'tracking', i=i+1; tracking = varargin{i}; Daniel@0: case 'weights', i=i+1; weights = varargin{i}; Daniel@0: case 'radius_ini', i=i+1; sTrain.radius_ini = varargin{i}; Daniel@0: case 'radius_fin', i=i+1; sTrain.radius_fin = varargin{i}; Daniel@0: case 'radius', Daniel@0: i=i+1; Daniel@0: l = length(varargin{i}); Daniel@0: if l==1, Daniel@0: sTrain.radius_ini = varargin{i}; Daniel@0: else Daniel@0: sTrain.radius_ini = varargin{i}(1); Daniel@0: sTrain.radius_fin = varargin{i}(end); Daniel@0: if l>2, radius = varargin{i}; end Daniel@0: end Daniel@0: case {'sTrain','train','som_train'}, i=i+1; sTrain = varargin{i}; Daniel@0: case {'topol','sTopol','som_topol'}, Daniel@0: i=i+1; Daniel@0: sTopol = varargin{i}; Daniel@0: if prod(sTopol.msize) ~= munits, Daniel@0: error('Given map grid size does not match the codebook size.'); Daniel@0: end Daniel@0: % unambiguous values Daniel@0: case {'hexa','rect'}, sTopol.lattice = varargin{i}; Daniel@0: case {'sheet','cyl','toroid'}, sTopol.shape = varargin{i}; Daniel@0: case {'gaussian','cutgauss','ep','bubble'}, sTrain.neigh = varargin{i}; Daniel@0: otherwise argok=0; Daniel@0: end Daniel@0: elseif isstruct(varargin{i}) & isfield(varargin{i},'type'), Daniel@0: switch varargin{i}(1).type, Daniel@0: case 'som_topol', Daniel@0: sTopol = varargin{i}; Daniel@0: if prod(sTopol.msize) ~= munits, Daniel@0: error('Given map grid size does not match the codebook size.'); Daniel@0: end Daniel@0: case 'som_train', sTrain = varargin{i}; Daniel@0: otherwise argok=0; Daniel@0: end Daniel@0: else Daniel@0: argok = 0; Daniel@0: end Daniel@0: if ~argok, Daniel@0: disp(['(som_batchtrain) Ignoring invalid argument #' num2str(i+2)]); Daniel@0: end Daniel@0: i = i+1; Daniel@0: end Daniel@0: Daniel@0: % take only weights of non-empty vectors Daniel@0: if length(weights)>dlen, weights = weights(nonempty); end Daniel@0: Daniel@0: % trainlen Daniel@0: if ~isempty(radius), sTrain.trainlen = length(radius); end Daniel@0: Daniel@0: % check topology Daniel@0: if struct_mode, Daniel@0: if ~strcmp(sTopol.lattice,sMap.topol.lattice) | ... Daniel@0: ~strcmp(sTopol.shape,sMap.topol.shape) | ... Daniel@0: any(sTopol.msize ~= sMap.topol.msize), Daniel@0: warning('Changing the original map topology.'); Daniel@0: end Daniel@0: end Daniel@0: sMap.topol = sTopol; Daniel@0: Daniel@0: % complement the training struct Daniel@0: sTrain = som_train_struct(sTrain,sMap,'dlen',dlen); Daniel@0: if isempty(sTrain.mask), sTrain.mask = ones(dim,1); end Daniel@0: Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Daniel@0: %% initialize Daniel@0: Daniel@0: M = sMap.codebook; Daniel@0: mask = sTrain.mask; Daniel@0: trainlen = sTrain.trainlen; Daniel@0: Daniel@0: % neighborhood radius Daniel@0: if trainlen==1, Daniel@0: radius = sTrain.radius_ini; Daniel@0: elseif length(radius)<=2, Daniel@0: r0 = sTrain.radius_ini; r1 = sTrain.radius_fin; Daniel@0: radius = r1 + fliplr((0:(trainlen-1))/(trainlen-1)) * (r0 - r1); Daniel@0: else Daniel@0: % nil Daniel@0: end Daniel@0: Daniel@0: % distance between map units in the output space Daniel@0: % Since in the case of gaussian and ep neighborhood functions, the Daniel@0: % equations utilize squares of the unit distances and in bubble case Daniel@0: % it doesn't matter which is used, the unitdistances and neighborhood Daniel@0: % radiuses are squared. Daniel@0: Ud = som_unit_dists(sTopol); Daniel@0: Ud = Ud.^2; Daniel@0: radius = radius.^2; Daniel@0: % zero neighborhood radius may cause div-by-zero error Daniel@0: radius(find(radius==0)) = eps; Daniel@0: Daniel@0: % The training algorithm involves calculating weighted Euclidian distances Daniel@0: % to all map units for each data vector. Basically this is done as Daniel@0: % for i=1:dlen, Daniel@0: % for j=1:munits, Daniel@0: % for k=1:dim Daniel@0: % Dist(j,i) = Dist(j,i) + mask(k) * (D(i,k) - M(j,k))^2; Daniel@0: % end Daniel@0: % end Daniel@0: % end Daniel@0: % where mask is the weighting vector for distance calculation. However, taking Daniel@0: % into account that distance between vectors m and v can be expressed as Daniel@0: % |m - v|^2 = sum_i ((m_i - v_i)^2) = sum_i (m_i^2 + v_i^2 - 2*m_i*v_i) Daniel@0: % this can be made much faster by transforming it to a matrix operation: Daniel@0: % Dist = (M.^2)*mask*ones(1,d) + ones(m,1)*mask'*(D'.^2) - 2*M*diag(mask)*D' Daniel@0: % Of the involved matrices, several are constant, as the mask and data do Daniel@0: % not change during training. Therefore they are calculated beforehand. Daniel@0: Daniel@0: % For the case where there are unknown components in the data, each data Daniel@0: % vector will have an individual mask vector so that for that unit, the Daniel@0: % unknown components are not taken into account in distance calculation. Daniel@0: % In addition all NaN's are changed to zeros so that they don't screw up Daniel@0: % the matrix multiplications and behave correctly in updating step. Daniel@0: Known = ~isnan(D); Daniel@0: W1 = (mask*ones(1,dlen)) .* Known'; Daniel@0: D(find(~Known)) = 0; Daniel@0: Daniel@0: % constant matrices Daniel@0: WD = 2*diag(mask)*D'; % constant matrix Daniel@0: dconst = ((D.^2)*mask)'; % constant in distance calculation for each data sample Daniel@0: % W2 = ones(munits,1)*mask'; D2 = (D'.^2); Daniel@0: Daniel@0: % initialize tracking Daniel@0: start = clock; Daniel@0: qe = zeros(trainlen,1); Daniel@0: Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Daniel@0: %% Action Daniel@0: Daniel@0: % With the 'blen' parameter you can control the memory consumption Daniel@0: % of the algorithm, which is in practive directly proportional Daniel@0: % to munits*blen. If you're having problems with memory, try to Daniel@0: % set the value of blen lower. Daniel@0: blen = min(munits,dlen); Daniel@0: Daniel@0: % reserve some space Daniel@0: bmus = zeros(1,dlen); Daniel@0: ddists = zeros(1,dlen); Daniel@0: Daniel@0: for t = 1:trainlen, Daniel@0: Daniel@0: % batchy train - this is done a block of data (inds) at a time Daniel@0: % rather than in a single sweep to save memory consumption. Daniel@0: % The 'Dist' and 'Hw' matrices have size munits*blen Daniel@0: % which - if you have a lot of data - would be HUGE if you Daniel@0: % calculated it all at once. A single-sweep version would Daniel@0: % look like this: Daniel@0: % Dist = (M.^2)*W1 - M*WD; %+ W2*D2 Daniel@0: % [ddists, bmus] = min(Dist); Daniel@0: % (notice that the W2*D2 term can be ignored since it is constant) Daniel@0: % This "batchy" version is the same as single-sweep if blen=dlen. Daniel@0: i0 = 0; Daniel@0: while i0+1<=dlen, Daniel@0: inds = [(i0+1):min(dlen,i0+blen)]; i0 = i0+blen; Daniel@0: Dist = (M.^2)*W1(:,inds) - M*WD(:,inds); Daniel@0: [ddists(inds), bmus(inds)] = min(Dist); Daniel@0: end Daniel@0: Daniel@0: % tracking Daniel@0: if tracking > 0, Daniel@0: ddists = ddists+dconst; % add the constant term Daniel@0: ddists(ddists<0) = 0; % rounding errors... Daniel@0: qe(t) = mean(sqrt(ddists)); Daniel@0: trackplot(M,D,tracking,start,t,qe); Daniel@0: end Daniel@0: Daniel@0: % neighborhood Daniel@0: % notice that the elements Ud and radius have been squared! Daniel@0: % note: 'bubble' matches the original "Batch Map" algorithm Daniel@0: switch sTrain.neigh, Daniel@0: case 'bubble', H = (Ud<=radius(t)); Daniel@0: case 'gaussian', H = exp(-Ud/(2*radius(t))); Daniel@0: case 'cutgauss', H = exp(-Ud/(2*radius(t))) .* (Ud<=radius(t)); Daniel@0: case 'ep', H = (1-Ud/radius(t)) .* (Ud<=radius(t)); Daniel@0: end Daniel@0: Daniel@0: % update Daniel@0: Daniel@0: % In principle the updating step goes like this: replace each map unit Daniel@0: % by the average of the data vectors that were in its neighborhood. Daniel@0: % The contribution, or activation, of data vectors in the mean can Daniel@0: % be varied with the neighborhood function. This activation is given Daniel@0: % by matrix H. So, for each map unit the new weight vector is Daniel@0: % Daniel@0: % m = sum_i (h_i * d_i) / sum_i (h_i), Daniel@0: % Daniel@0: % where i denotes the index of data vector. Since the values of Daniel@0: % neighborhood function h_i are the same for all data vectors belonging to Daniel@0: % the Voronoi set of the same map unit, the calculation is actually done Daniel@0: % by first calculating a partition matrix P with elements p_ij=1 if the Daniel@0: % BMU of data vector j is i. Daniel@0: Daniel@0: P = sparse(bmus,[1:dlen],weights,munits,dlen); Daniel@0: Daniel@0: % Then the sum of vectors in each Voronoi set are calculated (P*D) and the Daniel@0: % neighborhood is taken into account by calculating a weighted sum of the Daniel@0: % Voronoi sum (H*). The "activation" matrix A is the denominator of the Daniel@0: % equation above. Daniel@0: Daniel@0: S = H*(P*D); Daniel@0: A = H*(P*Known); Daniel@0: Daniel@0: % If you'd rather make this without using the Voronoi sets try the following: Daniel@0: % Hi = H(:,bmus); Daniel@0: % S = Hi * D; % "sum_i (h_i * d_i)" Daniel@0: % A = Hi * Known; % "sum_i (h_i)" Daniel@0: % The bad news is that the matrix Hi has size [munits x dlen]... Daniel@0: Daniel@0: % only update units for which the "activation" is nonzero Daniel@0: nonzero = find(A > 0); Daniel@0: M(nonzero) = S(nonzero) ./ A(nonzero); Daniel@0: Daniel@0: end; % for t = 1:trainlen Daniel@0: Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Daniel@0: %% Build / clean up the return arguments Daniel@0: Daniel@0: % tracking Daniel@0: if tracking > 0, fprintf(1,'\n'); end Daniel@0: Daniel@0: % update structures Daniel@0: sTrain = som_set(sTrain,'time',datestr(now,0)); Daniel@0: if struct_mode, Daniel@0: sMap = som_set(sMap,'codebook',M,'mask',sTrain.mask,'neigh',sTrain.neigh); Daniel@0: tl = length(sMap.trainhist); Daniel@0: sMap.trainhist(tl+1) = sTrain; Daniel@0: else Daniel@0: sMap = reshape(M,orig_size); Daniel@0: end Daniel@0: Daniel@0: return; Daniel@0: Daniel@0: Daniel@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Daniel@0: %% subfunctions Daniel@0: Daniel@0: %%%%%%%% Daniel@0: function [] = trackplot(M,D,tracking,start,n,qe) Daniel@0: Daniel@0: l = length(qe); Daniel@0: elap_t = etime(clock,start); Daniel@0: tot_t = elap_t*l/n; Daniel@0: fprintf(1,'\rTraining: %3.0f/ %3.0f s',elap_t,tot_t) Daniel@0: switch tracking Daniel@0: case 1, Daniel@0: case 2, Daniel@0: plot(1:n,qe(1:n),(n+1):l,qe((n+1):l)) Daniel@0: title('Quantization error after each epoch'); Daniel@0: drawnow Daniel@0: otherwise, Daniel@0: subplot(2,1,1), plot(1:n,qe(1:n),(n+1):l,qe((n+1):l)) Daniel@0: title('Quantization error after each epoch'); Daniel@0: subplot(2,1,2), plot(M(:,1),M(:,2),'ro',D(:,1),D(:,2),'b+'); Daniel@0: title('First two components of map units (o) and data vectors (+)'); Daniel@0: drawnow Daniel@0: end Daniel@0: % end of trackplot