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