som_batchtrain.m

酒吧老板价格表vbkkljl个ipfgik皮带 东风康帕斯飞 帕斯克哦佛i 日减肥假话骗人 欧冠欧赔健康的飞机

  if ischar(varargin{i}), 
    switch varargin{i}, 
     % argument IDs
     case 'msize', i=i+1; sTopol.msize = varargin{i}; 
     case 'lattice', i=i+1; sTopol.lattice = varargin{i};
     case 'shape', i=i+1; sTopol.shape = varargin{i};
     case 'mask', i=i+1; sTrain.mask = varargin{i};
     case 'neigh', i=i+1; sTrain.neigh = varargin{i};
     case 'trainlen', i=i+1; sTrain.trainlen = varargin{i};
     case 'tracking', i=i+1; tracking = varargin{i};
     case 'weights', i=i+1; weights = varargin{i}; 
     case 'radius_ini', i=i+1; sTrain.radius_ini = varargin{i};
     case 'radius_fin', i=i+1; sTrain.radius_fin = varargin{i};
     case 'radius', 
      i=i+1; 
      l = length(varargin{i}); 
      if l==1, 
        sTrain.radius_ini = varargin{i}; 
      else 
        sTrain.radius_ini = varargin{i}(1); 
        sTrain.radius_fin = varargin{i}(end);
        if l>2, radius = varargin{i}; end
      end 
     case {'sTrain','train','som_train'}, i=i+1; sTrain = varargin{i};
     case {'topol','sTopol','som_topol'}, 
      i=i+1; 
      sTopol = varargin{i};
      if prod(sTopol.msize) ~= munits, 
        error('Given map grid size does not match the codebook size.');
      end
      % unambiguous values
     case {'hexa','rect'}, sTopol.lattice = varargin{i};
     case {'sheet','cyl','toroid'}, sTopol.shape = varargin{i}; 
     case {'gaussian','cutgauss','ep','bubble'}, sTrain.neigh = varargin{i};
     otherwise argok=0; 
    end
  elseif isstruct(varargin{i}) & isfield(varargin{i},'type'), 
    switch varargin{i}(1).type, 
     case 'som_topol', 
      sTopol = varargin{i}; 
      if prod(sTopol.msize) ~= munits, 
        error('Given map grid size does not match the codebook size.');
      end
     case 'som_train', sTrain = varargin{i};
     otherwise argok=0; 
    end
  else
    argok = 0; 
  end
  if ~argok, 
    disp(['(som_batchtrain) Ignoring invalid argument #' num2str(i+2)]); 
  end
  i = i+1; 
end

% take only weights of non-empty vectors
if length(weights)>dlen, weights = weights(nonempty); end

% trainlen
if ~isempty(radius), sTrain.trainlen = length(radius); end

% check topology
if struct_mode, 
  if ~strcmp(sTopol.lattice,sMap.topol.lattice) | ...
	~strcmp(sTopol.shape,sMap.topol.shape) | ...
	any(sTopol.msize ~= sMap.topol.msize), 
    warning('Changing the original map topology.');
  end
end
sMap.topol = sTopol; 

% complement the training struct
sTrain = som_train_struct(sTrain,sMap,'dlen',dlen);
if isempty(sTrain.mask), sTrain.mask = ones(dim,1); end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% initialize

M        = sMap.codebook;
mask     = sTrain.mask;
trainlen = sTrain.trainlen;

% neighborhood radius
if trainlen==1, 
  radius = sTrain.radius_ini; 
elseif length(radius)<=2,  
  r0 = sTrain.radius_ini; r1 = sTrain.radius_fin;
  radius = r1 + fliplr((0:(trainlen-1))/(trainlen-1)) * (r0 - r1);
else
  % nil
end
                                   
% distance between map units in the output space
%  Since in the case of gaussian and ep neighborhood functions, the 
%  equations utilize squares of the unit distances and in bubble case
%  it doesn't matter which is used, the unitdistances and neighborhood
%  radiuses are squared.
Ud = som_unit_dists(sTopol);
Ud = Ud.^2;
radius = radius.^2;
% zero neighborhood radius may cause div-by-zero error
radius(find(radius==0)) = eps; 

% The training algorithm involves calculating weighted Euclidian distances 
% to all map units for each data vector. Basically this is done as
%   for i=1:dlen, 
%     for j=1:munits, 
%       for k=1:dim
%         Dist(j,i) = Dist(j,i) + mask(k) * (D(i,k) - M(j,k))^2;
%       end
%     end
%   end
% where mask is the weighting vector for distance calculation. However, taking 
% into account that distance between vectors m and v can be expressed as
%   |m - v|^2 = sum_i ((m_i - v_i)^2) = sum_i (m_i^2 + v_i^2 - 2*m_i*v_i)
% this can be made much faster by transforming it to a matrix operation:
%   Dist = (M.^2)*mask*ones(1,d) + ones(m,1)*mask'*(D'.^2) - 2*M*diag(mask)*D'
% Of the involved matrices, several are constant, as the mask and data do 
% not change during training. Therefore they are calculated beforehand.

% For the case where there are unknown components in the data, each data
% vector will have an individual mask vector so that for that unit, the 
% unknown components are not taken into account in distance calculation.
% In addition all NaN's are changed to zeros so that they don't screw up 
% the matrix multiplications and behave correctly in updating step.
Known = ~isnan(D);
W1 = (mask*ones(1,dlen)) .* Known'; 
D(find(~Known)) = 0;  

% constant matrices
WD = 2*diag(mask)*D';    % constant matrix
dconst = ((D.^2)*mask)'; % constant in distance calculation for each data sample 
                         % W2 = ones(munits,1)*mask'; D2 = (D'.^2); 		      

% initialize tracking
start = clock;
qe = zeros(trainlen,1); 
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Action

% With the 'blen' parameter you can control the memory consumption 
% of the algorithm, which is in practive directly proportional
% to munits*blen. If you're having problems with memory, try to 
% set the value of blen lower. 
blen = min(munits,dlen);

% reserve some space
bmus = zeros(1,dlen); 
ddists = zeros(1,dlen); 

for t = 1:trainlen,  

  % batchy train - this is done a block of data (inds) at a time
  % rather than in a single sweep to save memory consumption. 
  % The 'Dist' and 'Hw' matrices have size munits*blen
  % which - if you have a lot of data - would be HUGE if you 
  % calculated it all at once. A single-sweep version would 
  % look like this: 
  %  Dist = (M.^2)*W1 - M*WD; %+ W2*D2 
  %  [ddists, bmus] = min(Dist);
  % (notice that the W2*D2 term can be ignored since it is constant)
  % This "batchy" version is the same as single-sweep if blen=dlen. 
  i0 = 0;     
  while i0+1<=dlen, 
    inds = [(i0+1):min(dlen,i0+blen)]; i0 = i0+blen;      
    Dist = (M.^2)*W1(:,inds) - M*WD(:,inds);
    [ddists(inds), bmus(inds)] = min(Dist);
  end  
  
  % tracking
  if tracking > 0,
    ddists = ddists+dconst; % add the constant term
    ddists(ddists<0) = 0;   % rounding errors...
    qe(t) = mean(sqrt(ddists));
    trackplot(M,D,tracking,start,t,qe);
  end
  
  % neighborhood 
  % notice that the elements Ud and radius have been squared!
  % note: 'bubble' matches the original "Batch Map" algorithm
  switch sTrain.neigh, 
   case 'bubble',   H = (Ud<=radius(t)); 
   case 'gaussian', H = exp(-Ud/(2*radius(t))); 
   case 'cutgauss', H = exp(-Ud/(2*radius(t))) .* (Ud<=radius(t));
   case 'ep',       H = (1-Ud/radius(t)) .* (Ud<=radius(t));
  end  
  
  % update 

  % In principle the updating step goes like this: replace each map unit 
  % by the average of the data vectors that were in its neighborhood.
  % The contribution, or activation, of data vectors in the mean can 
  % be varied with the neighborhood function. This activation is given 
  % by matrix H. So, for each map unit the new weight vector is
  %
  %      m = sum_i (h_i * d_i) / sum_i (h_i), 
  % 
  % where i denotes the index of data vector.  Since the values of
  % neighborhood function h_i are the same for all data vectors belonging to
  % the Voronoi set of the same map unit, the calculation is actually done
  % by first calculating a partition matrix P with elements p_ij=1 if the
  % BMU of data vector j is i.

  P = sparse(bmus,[1:dlen],weights,munits,dlen);
       
  % Then the sum of vectors in each Voronoi set are calculated (P*D) and the
  % neighborhood is taken into account by calculating a weighted sum of the
  % Voronoi sum (H*). The "activation" matrix A is the denominator of the 
  % equation above.
  
  S = H*(P*D); 
  A = H*(P*Known);
  
  % If you'd rather make this without using the Voronoi sets try the following: 
  %   Hi = H(:,bmus); 
  %   S = Hi * D;            % "sum_i (h_i * d_i)"
  %   A = Hi * Known;        % "sum_i (h_i)"
  % The bad news is that the matrix Hi has size [munits x dlen]... 
    
  % only update units for which the "activation" is nonzero
  nonzero = find(A > 0); 
  M(nonzero) = S(nonzero) ./ A(nonzero); 

end; % for t = 1:trainlen

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Build / clean up the return arguments

% tracking
if tracking > 0, fprintf(1,'\n'); end

% update structures
sTrain = som_set(sTrain,'time',datestr(now,0));
if struct_mode, 
  sMap = som_set(sMap,'codebook',M,'mask',sTrain.mask,'neigh',sTrain.neigh);
  tl = length(sMap.trainhist);
  sMap.trainhist(tl+1) = sTrain;
else
  sMap = reshape(M,orig_size);
end

return;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% subfunctions

%%%%%%%%
function [] = trackplot(M,D,tracking,start,n,qe)

  l = length(qe);
  elap_t = etime(clock,start); 
  tot_t = elap_t*l/n;
  fprintf(1,'\rTraining: %3.0f/ %3.0f s',elap_t,tot_t)  
  switch tracking
   case 1, 
   case 2,   
    plot(1:n,qe(1:n),(n+1):l,qe((n+1):l))
    title('Quantization error after each epoch');
    drawnow
   otherwise,
    subplot(2,1,1), plot(1:n,qe(1:n),(n+1):l,qe((n+1):l))
    title('Quantization error after each epoch');
    subplot(2,1,2), plot(M(:,1),M(:,2),'ro',D(:,1),D(:,2),'b+'); 
    title('First two components of map units (o) and data vectors (+)');
    drawnow
  end
  % end of trackplot