wolffd@0: function [sMap, sTrain] = som_seqtrain(sMap, D, varargin)
wolffd@0: 
wolffd@0: %SOM_SEQTRAIN  Use sequential algorithm to train the Self-Organizing Map.
wolffd@0: %
wolffd@0: % [sM,sT] = som_seqtrain(sM, D, [[argID,] value, ...])
wolffd@0: % 
wolffd@0: %  sM     = som_seqtrain(sM,D);
wolffd@0: %  sM     = som_seqtrain(sM,sD,'alpha_type','power','tracking',3);
wolffd@0: %  [M,sT] = som_seqtrain(M,D,'ep','trainlen',10,'inv','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: %   'alpha'       (vector) learning rates, length trainlen
wolffd@0: %   'alpha_ini'   (scalar) initial learning rate
wolffd@0: %   'tracking'    (scalar) tracking level, 0-3 
wolffd@0: %   'trainlen'    (scalar) training length
wolffd@0: %   'trainlen_type' *(string) is the given trainlen 'samples' or 'epochs'
wolffd@0: %   'train'      *(struct) train struct, parameters for training
wolffd@0: %   'sTrain','som_train '  = 'train'
wolffd@0: %   'alpha_type' *(string) learning rate function, 'inv', 'linear' or 'power'
wolffd@0: %   'sample_order'*(string) order of samples: 'random' or 'ordered'
wolffd@0: %   'neigh'      *(string) neighborhood function, 'gaussian', 'cutgauss',
wolffd@0: %                          'ep' or 'bubble'
wolffd@0: %   'topol'      *(struct) topology struct
wolffd@0: %   'som_topol','sTopo l'  = 'topol'
wolffd@0: %   'lattice'    *(string) map lattice, 'hexa' or 'rect'
wolffd@0: %   'shape'      *(string) map shape, 'sheet', 'cyl' or 'toroid'
wolffd@0: %
wolffd@0: % For more help, try 'type som_seqtrain' or check out online documentation.
wolffd@0: % See also  SOM_MAKE, SOM_BATCHTRAIN, SOM_TRAIN_STRUCT.
wolffd@0: 
wolffd@0: %%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %
wolffd@0: % som_seqtrain
wolffd@0: %
wolffd@0: % PURPOSE
wolffd@0: %
wolffd@0: % Trains a Self-Organizing Map using the sequential algorithm. 
wolffd@0: %
wolffd@0: % SYNTAX
wolffd@0: %
wolffd@0: %  sM = som_seqtrain(sM,D);
wolffd@0: %  sM = som_seqtrain(sM,sD);
wolffd@0: %  sM = som_seqtrain(...,'argID',value,...);
wolffd@0: %  sM = som_seqtrain(...,value,...);
wolffd@0: %  [sM,sT] = som_seqtrain(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 sequential SOM training algorithm. If no optional
wolffd@0: % arguments (argID, value) are given, a default training is done, the
wolffd@0: % parameters are obtained from SOM_TRAIN_STRUCT function. Using
wolffd@0: % optional arguments the training parameters can be specified. Returns
wolffd@0: % the 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. 78-82.
wolffd@0: % Kohonen, T., "Clustering, Taxonomy, and Topological Maps of
wolffd@0: %    Patterns", International Conference on Pattern Recognition
wolffd@0: %    (ICPR), 1982, pp. 114-128.
wolffd@0: % Kohonen, T., "Self-Organized formation of topologically correct
wolffd@0: %    feature maps", Biological Cybernetics 43, 1982, pp. 59-69.
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: %   '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: %   'alpha'      (vector) learning rate
wolffd@0: %                         length = 1: alpha_ini = alpha
wolffd@0: %                         length > 1: the vector gives learning rate
wolffd@0: %                                     for each step separately
wolffd@0: %                                     trainlen is set to length(alpha)
wolffd@0: %                                     alpha_type is set to 'user defined'
wolffd@0: %   'alpha_ini'  (scalar) initial learning rate
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 (see also 'tlen_type')
wolffd@0: %   'trainlen_type' *(string) is the trainlen argument given in 'epochs'
wolffd@0: %                         or in 'samples'. Default is 'epochs'.
wolffd@0: %   'sample_order'*(string) is the sample order 'random' (which is the 
wolffd@0: %                         the default) or 'ordered' in which case
wolffd@0: %                         samples are taken in the order in which they 
wolffd@0: %                         appear in the data set
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', 'alpha_type',
wolffd@0: %                         'alpha_ini', 'radius_ini' and 'radius_fin').
wolffd@0: %   'sTrain', 'som_train' (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: %   'alpha_type'*(string) learning rate function, 'inv', 'linear' or 'power'
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: %   
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_seqtrain(sM,D);  
wolffd@0: %  sM = som_seqtrain(sM,sD);  
wolffd@0: %
wolffd@0: % To change the tracking level, 'tracking' argument is specified:
wolffd@0: %  sM = som_seqtrain(sM,D,'tracking',3);
wolffd@0: %
wolffd@0: % The change training parameters, the optional arguments 'train', 
wolffd@0: % 'neigh','mask','trainlen','radius','radius_ini', 'radius_fin', 
wolffd@0: % 'alpha', 'alpha_type' and 'alpha_ini' are used. 
wolffd@0: %  sM = som_seqtrain(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),'algorithm','seq');
wolffd@0: %  sTrain = som_set(sTrain,'neigh','cutgauss');
wolffd@0: %  sM = som_seqtrain(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_seqtrain(sM,D,'radius',[5 3 1 1 1 1 0.5 0.5 0.5]);
wolffd@0: %
wolffd@0: % By default the learning rate (alpha) goes from the alpha_ini to 0
wolffd@0: % along the function defined by alpha_type. If you want to change this, 
wolffd@0: % you can use the 'alpha' argument to specify the learning rate
wolffd@0: % for each step separately: 
wolffd@0: %  alpha = 0.2*(1 - log([1:100]));
wolffd@0: %  sM = som_seqtrain(sM,D,'alpha',alpha);
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_seqtrain(M,D,'msize',[20 10],'lattice','hexa','shape','cyl');
wolffd@0: %
wolffd@0: %  M = som_seqtrain(M,D,'msize',[20 10],'hexa','cyl');
wolffd@0: %
wolffd@0: %  sT= som_set('som_topol','msize',[20 10],'lattice','hexa','shape','cyl');
wolffd@0: %  M = som_seqtrain(M,D,sT);
wolffd@0: %
wolffd@0: %  M = reshape(M,[20 10 dim]);
wolffd@0: %  M = som_seqtrain(M,D,'hexa','cyl');
wolffd@0: %
wolffd@0: % The som_seqtrain also returns a train struct with information on the 
wolffd@0: % accomplished training. This is the same one as is added to the end of the 
wolffd@0: % trainhist field of map struct, in case a map struct is given.
wolffd@0: %  [M,sTrain] = som_seqtrain(M,D,'msize',[20 10]);
wolffd@0: %
wolffd@0: %  [sM,sTrain] = som_seqtrain(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_batchtrain   Train SOM with batch 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 220997
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: D = D(find(sum(isnan(D),2) < dim),:); % remove empty vectors from the data
wolffd@0: [dlen ddim] = size(D);                % check input dimension
wolffd@0: if dim ~= ddim, error('Map and data input space dimensions disagree.'); end
wolffd@0: 
wolffd@0: % varargin
wolffd@0: sTrain = som_set('som_train','algorithm','seq','neigh', ...
wolffd@0: 		 sMap.neigh,'mask',sMap.mask,'data_name',data_name);
wolffd@0: radius     = [];
wolffd@0: alpha      = [];
wolffd@0: tracking   = 1;
wolffd@0: sample_order_type = 'random';
wolffd@0: tlen_type  = 'epochs';
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 'trainlen_type', i=i+1; tlen_type = varargin{i}; 
wolffd@0:      case 'tracking', i=i+1; tracking = varargin{i};
wolffd@0:      case 'sample_order', i=i+1; sample_order_type = 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}; tlen_type = 'samples'; end
wolffd@0:       end 
wolffd@0:      case 'alpha_type', i=i+1; sTrain.alpha_type = varargin{i};
wolffd@0:      case 'alpha_ini', i=i+1; sTrain.alpha_ini = varargin{i};
wolffd@0:      case 'alpha',     
wolffd@0:       i=i+1; 
wolffd@0:       sTrain.alpha_ini = varargin{i}(1);
wolffd@0:       if length(varargin{i})>1, 
wolffd@0:         alpha = varargin{i}; tlen_type = 'samples'; 
wolffd@0:         sTrain.alpha_type = 'user defined'; 
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 {'inv','linear','power'}, sTrain.alpha_type = varargin{i}; 
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:      case {'epochs','samples'}, tlen_type = varargin{i};
wolffd@0:      case {'random', 'ordered'}, sample_order_type = 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_seqtrain) Ignoring invalid argument #' num2str(i+2)]); 
wolffd@0:   end
wolffd@0:   i = i+1; 
wolffd@0: end
wolffd@0: 
wolffd@0: % training length
wolffd@0: if ~isempty(radius) | ~isempty(alpha), 
wolffd@0:   lr = length(radius);
wolffd@0:   la = length(alpha);
wolffd@0:   if lr>2 | la>1,
wolffd@0:     tlen_type = 'samples';
wolffd@0:     if     lr> 2 & la<=1, sTrain.trainlen = lr;
wolffd@0:     elseif lr<=2 & la> 1, sTrain.trainlen = la;
wolffd@0:     elseif lr==la,        sTrain.trainlen = la;
wolffd@0:     else
wolffd@0:       error('Mismatch between radius and learning rate vector lengths.')
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: if strcmp(tlen_type,'samples'), sTrain.trainlen = sTrain.trainlen/dlen; 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*dlen;
wolffd@0: 
wolffd@0: % neighborhood radius
wolffd@0: if length(radius)>2,
wolffd@0:   radius_type = 'user defined';
wolffd@0: else
wolffd@0:   radius = [sTrain.radius_ini sTrain.radius_fin];    
wolffd@0:   rini = radius(1); 
wolffd@0:   rstep = (radius(end)-radius(1))/(trainlen-1);
wolffd@0:   radius_type = 'linear';
wolffd@0: end    
wolffd@0: 
wolffd@0: % learning rate
wolffd@0: if length(alpha)>1, 
wolffd@0:   sTrain.alpha_type ='user defined';
wolffd@0:   if length(alpha) ~= trainlen, 
wolffd@0:     error('Trainlen and length of neighborhood radius vector do not match.')
wolffd@0:   end
wolffd@0:   if any(isnan(alpha)), 
wolffd@0:     error('NaN is an illegal learning rate.')
wolffd@0:   end
wolffd@0: else
wolffd@0:   if isempty(alpha), alpha = sTrain.alpha_ini; end
wolffd@0:   if strcmp(sTrain.alpha_type,'inv'), 
wolffd@0:     % alpha(t) = a / (t+b), where a and b are chosen suitably
wolffd@0:     % below, they are chosen so that alpha_fin = alpha_ini/100
wolffd@0:     b = (trainlen - 1) / (100 - 1);
wolffd@0:     a = b * alpha;
wolffd@0:   end
wolffd@0: end
wolffd@0:                                    
wolffd@0: % initialize random number generator
wolffd@0: rand('state',sum(100*clock));
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).^2;
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Action
wolffd@0: 
wolffd@0: update_step = 100; 
wolffd@0: mu_x_1 = ones(munits,1);
wolffd@0: samples = ones(update_step,1);
wolffd@0: r = samples; 
wolffd@0: alfa = samples;
wolffd@0: 
wolffd@0: qe = 0;
wolffd@0: start = clock;
wolffd@0: if tracking >  0, % initialize tracking
wolffd@0:   track_table = zeros(update_step,1);
wolffd@0:   qe = zeros(floor(trainlen/update_step),1);  
wolffd@0: end
wolffd@0: 
wolffd@0: for t = 1:trainlen, 
wolffd@0: 
wolffd@0:   % Every update_step, new values for sample indeces, neighborhood
wolffd@0:   % radius and learning rate are calculated. This could be done
wolffd@0:   % every step, but this way it is more efficient. Or this could 
wolffd@0:   % be done all at once outside the loop, but it would require much
wolffd@0:   % more memory.
wolffd@0:   ind = rem(t,update_step); if ind==0, ind = update_step; end
wolffd@0:   if ind==1, 
wolffd@0:     steps = [t:min(trainlen,t+update_step-1)];
wolffd@0:     % sample order    
wolffd@0:     switch sample_order_type, 
wolffd@0:      case 'ordered', samples = rem(steps,dlen)+1;
wolffd@0:      case 'random',  samples = ceil(dlen*rand(update_step,1)+eps);
wolffd@0:     end
wolffd@0: 
wolffd@0:     % neighborhood radius
wolffd@0:     switch radius_type, 
wolffd@0:      case 'linear',       r = rini+(steps-1)*rstep;
wolffd@0:      case 'user defined', r = radius(steps); 
wolffd@0:     end    
wolffd@0:     r=r.^2;        % squared radius (see notes about Ud above)
wolffd@0:     r(r==0) = eps; % zero radius might cause div-by-zero error
wolffd@0:     
wolffd@0:     % learning rate
wolffd@0:     switch sTrain.alpha_type,
wolffd@0:      case 'linear',       alfa = (1-steps/trainlen)*alpha;
wolffd@0:      case 'inv',          alfa = a ./ (b + steps-1);
wolffd@0:      case 'power',        alfa = alpha * (0.005/alpha).^((steps-1)/trainlen); 
wolffd@0:      case 'user defined', alfa = alpha(steps);
wolffd@0:     end    
wolffd@0:   end
wolffd@0:   
wolffd@0:   % find BMU
wolffd@0:   x = D(samples(ind),:);                 % pick one sample vector
wolffd@0:   known = ~isnan(x);                     % its known components
wolffd@0:   Dx = M(:,known) - x(mu_x_1,known);     % each map unit minus the vector
wolffd@0:   [qerr bmu] = min((Dx.^2)*mask(known)); % minimum distance(^2) and the BMU
wolffd@0: 
wolffd@0:   % tracking
wolffd@0:   if tracking>0, 
wolffd@0:     track_table(ind) = sqrt(qerr);
wolffd@0:     if ind==update_step, 
wolffd@0:       n = ceil(t/update_step); 
wolffd@0:       qe(n) = mean(track_table);
wolffd@0:       trackplot(M,D,tracking,start,n,qe);
wolffd@0:     end
wolffd@0:   end
wolffd@0:   
wolffd@0:   % neighborhood & learning rate
wolffd@0:   % notice that the elements Ud and radius have been squared!
wolffd@0:   % (see notes about Ud above)
wolffd@0:   switch sTrain.neigh, 
wolffd@0:   case 'bubble',   h = (Ud(:,bmu)<=r(ind));
wolffd@0:   case 'gaussian', h = exp(-Ud(:,bmu)/(2*r(ind))); 
wolffd@0:   case 'cutgauss', h = exp(-Ud(:,bmu)/(2*r(ind))) .* (Ud(:,bmu)<=r(ind));
wolffd@0:   case 'ep',       h = (1-Ud(:,bmu)/r(ind)) .* (Ud(:,bmu)<=r(ind));
wolffd@0:   end  
wolffd@0:   h = h*alfa(ind);  
wolffd@0:   
wolffd@0:   % update M
wolffd@0:   M(:,known) = M(:,known) - h(:,ones(sum(known),1)).*Dx;
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: if tracking, 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: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 errors for latest samples')    
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 for latest samples');
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
wolffd@0: