wolffd@0: function [sM,sTrain] = som_prototrain(sM, D)
wolffd@0: 
wolffd@0: %SOM_PROTOTRAIN  Use sequential algorithm to train the Self-Organizing Map.
wolffd@0: %
wolffd@0: % [sM,sT] = som_prototrain(sM, D)
wolffd@0: % 
wolffd@0: %  sM = som_prototrain(sM,D);
wolffd@0: %
wolffd@0: %  Input and output arguments: 
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: %
wolffd@0: % This function is otherwise just like SOM_SEQTRAIN except that
wolffd@0: % the implementation of the sequential training algorithm is very 
wolffd@0: % straightforward (and slower). This should make it easy for you 
wolffd@0: % to modify the algorithm, if you want to. 
wolffd@0: %
wolffd@0: % For help on input and output parameters, try 
wolffd@0: % 'type som_prototrain' or check out the help for SOM_SEQTRAIN.
wolffd@0: % See also SOM_SEQTRAIN, SOM_BATCHTRAIN.
wolffd@0: 
wolffd@0: % Contributed to SOM Toolbox vs2, February 2nd, 2000 by Juha Vesanto
wolffd@0: % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0: 
wolffd@0: % Version 2.0beta juuso 080200 130300
wolffd@0:  
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Check input arguments
wolffd@0: 
wolffd@0: % map 
wolffd@0: struct_mode = isstruct(sM);
wolffd@0: if struct_mode, 
wolffd@0:   M = sM.codebook; 
wolffd@0:   sTopol = sM.topol; 
wolffd@0:   mask = sM.mask; 
wolffd@0:   msize = sTopol.msize;
wolffd@0:   neigh = sM.neigh;
wolffd@0: else  
wolffd@0:   M = sM; orig_size = size(M);
wolffd@0:   if ndims(sM) > 2, 
wolffd@0:     si = size(sM); dim = si(end); msize = si(1:end-1);
wolffd@0:     M = reshape(sM,[prod(msize) dim]);
wolffd@0:   else
wolffd@0:     msize = [orig_size(1) 1]; dim = orig_size(2);
wolffd@0:   end
wolffd@0:   sM = som_map_struct(dim,'msize',msize); sTopol = sM.topol;
wolffd@0:   mask = ones(dim,1);
wolffd@0:   neigh = 'gaussian';
wolffd@0: end
wolffd@0: [munits dim] = size(M); 
wolffd@0: 
wolffd@0: % data
wolffd@0: if isstruct(D), data_name = D.name; D = D.data; 
wolffd@0: else 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: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% initialize (these are default values, change as you will)
wolffd@0: 
wolffd@0: % training length
wolffd@0: trainlen = 20*dlen; % 20 epochs by default
wolffd@0: 
wolffd@0: % neighborhood radius
wolffd@0: radius_type = 'linear';
wolffd@0: rini = max(msize)/2;
wolffd@0: rfin = 1;
wolffd@0: 
wolffd@0: % learning rate
wolffd@0: alpha_type = 'inv'; 
wolffd@0: alpha_ini = 0.2;
wolffd@0: 
wolffd@0: % initialize random number generator
wolffd@0: rand('state',sum(100*clock));
wolffd@0: 
wolffd@0: % tracking 
wolffd@0: start = clock; trackstep = 100;
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Action
wolffd@0: 
wolffd@0: Ud = som_unit_dists(sTopol); % distance between map units on the grid
wolffd@0: mu_x_1 = ones(munits,1);     % this is used pretty often
wolffd@0: 
wolffd@0: for t = 1:trainlen, 
wolffd@0: 
wolffd@0:   %% find BMU
wolffd@0:   ind = ceil(dlen*rand(1)+eps);       % select one vector
wolffd@0:   x = D(ind,:);                       % pick it up
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:   dist2 = (Dx.^2)*mask(known);        % squared distances  
wolffd@0:   [qerr bmu] = min(dist2);            % find BMU
wolffd@0: 
wolffd@0:   %% neighborhood  
wolffd@0:   switch radius_type, % radius
wolffd@0:    case 'linear', r = rini+(rfin-rini)*(t-1)/(trainlen-1);
wolffd@0:   end
wolffd@0:   if ~r, r=eps; end % zero neighborhood radius may cause div-by-zero error  
wolffd@0:   switch neigh, % neighborhood function 
wolffd@0:   case 'bubble',   h = (Ud(:,bmu) <= r);
wolffd@0:   case 'gaussian', h = exp(-(Ud(:,bmu).^2)/(2*r*r)); 
wolffd@0:   case 'cutgauss', h = exp(-(Ud(:,bmu).^2)/(2*r*r)) .* (Ud(:,bmu) <= r);
wolffd@0:   case 'ep',       h = (1 - (Ud(:,bmu).^2)/(r*r)) .* (Ud(:,bmu) <= r);
wolffd@0:   end  
wolffd@0: 
wolffd@0:   %% learning rate
wolffd@0:   switch alpha_type,
wolffd@0:    case 'linear', a = (1-t/trainlen)*alpha_ini;
wolffd@0:    case 'inv',    a = alpha_ini / (1 + 99*(t-1)/(trainlen-1));
wolffd@0:    case 'power',  a = alpha_ini * (0.005/alpha_ini)^((t-1)/trainlen); 
wolffd@0:   end
wolffd@0:   
wolffd@0:   %% update
wolffd@0:   M(:,known) = M(:,known) - a*h(:,ones(sum(known),1)).*Dx;
wolffd@0: 			 
wolffd@0:   %% tracking
wolffd@0:   if t==1 | ~rem(t,trackstep),
wolffd@0:     elap_t = etime(clock,start); tot_t = elap_t*trainlen/t; 
wolffd@0:     fprintf(1,'\rTraining: %3.0f/ %3.0f s',elap_t,tot_t)
wolffd@0:   end
wolffd@0:   
wolffd@0: end; % for t = 1:trainlen
wolffd@0: fprintf(1,'\n');
wolffd@0: 
wolffd@0: % outputs
wolffd@0: sTrain = som_set('som_train','algorithm','proto',...
wolffd@0: 		 'data_name',data_name,...
wolffd@0: 		 'neigh',neigh,...
wolffd@0: 		 'mask',mask,...
wolffd@0: 		 'radius_ini',rini,...
wolffd@0: 		 'radius_fin',rfin,...
wolffd@0: 		 'alpha_ini',alpha_ini,...
wolffd@0: 		 'alpha_type',alpha_type,...
wolffd@0: 		 'trainlen',trainlen,...
wolffd@0: 		 'time',datestr(now,0));
wolffd@0: 
wolffd@0: if struct_mode, 
wolffd@0:   sM = som_set(sM,'codebook',M,'mask',mask,'neigh',neigh);
wolffd@0:   sM.trainhist(end+1) = sTrain;
wolffd@0: else
wolffd@0:   sM = reshape(M,orig_size);
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: 
wolffd@0: 
wolffd@0: 
wolffd@0: