wolffd@0: function [hits] = som_hits(sMap, sData, mode)
wolffd@0: 
wolffd@0: %SOM_HITS Calculate the response of the given data on the map.
wolffd@0: %
wolffd@0: % hits = som_hits(sMap, sData, [mode])
wolffd@0: %
wolffd@0: %   h = som_hits(sMap,sData);
wolffd@0: %   h = som_hits(sMap,sData,'fuzzy');
wolffd@0: %
wolffd@0: %  Input and output arguments ([]'s are optional): 
wolffd@0: %   sMap     (struct) map struct
wolffd@0: %            (matrix) codebook matrix, size munits x dim
wolffd@0: %   sData    (struct) data struct
wolffd@0: %            (matrix) data matrix, size dlen x dim
wolffd@0: %   [mode]   (string) 'crisp' (default), 'kernel', 'fuzzy'
wolffd@0: %
wolffd@0: %   hits     (vector) the number of hits in each map unit, length = munits
wolffd@0: %
wolffd@0: % The response of the data on the map can be calculated e.g. in
wolffd@0: % three ways, selected with the mode argument: 
wolffd@0: %  'crisp'    traditional hit histogram
wolffd@0: %  'kernel'   a sum of dlen neighborhood kernels, where kernel
wolffd@0: %             is positioned on the BMU of each data sample. The 
wolffd@0: %             neighborhood function is sMap.neigh and the
wolffd@0: %             neighborhood width is sMap.trainhist(end).radius_fin
wolffd@0: %             or 1 if this is empty or NaN
wolffd@0: %  'fuzzy'    fuzzy response calculated by summing 1./(1+(q/a)^2)
wolffd@0: %             for each data sample, where q is a vector containing
wolffd@0: %             distance from the data sample to each map unit and 
wolffd@0: %             a is average quantization error
wolffd@0: % 
wolffd@0: % For more help, try 'type som_hits' or check out online documentation.
wolffd@0: % See also SOM_AUTOLABEL, SOM_BMUS.    
wolffd@0: 
wolffd@0: %%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %
wolffd@0: % som_hits
wolffd@0: %
wolffd@0: % PURPOSE
wolffd@0: %
wolffd@0: % Calculate the response of the given data on the map.
wolffd@0: %
wolffd@0: % SYNTAX
wolffd@0: %
wolffd@0: %  hits = som_hits(sMap, sData)
wolffd@0: %  hits = som_hits(M, D)
wolffd@0: %  hits = som_hits(..., mode)
wolffd@0: %
wolffd@0: % DESCRIPTION
wolffd@0: %
wolffd@0: % Returns a vector indicating the response of the map to the data.
wolffd@0: % The response of the data on the map can be calculated e.g. in
wolffd@0: % three ways, selected with the mode argument: 
wolffd@0: %  'crisp'    traditional hit histogram: how many times each map unit 
wolffd@0: %             was the BMU for the data set
wolffd@0: %  'kernel'   a sum of neighborhood kernels, where a kernel
wolffd@0: %             is positioned on the BMU of each data sample. The 
wolffd@0: %             neighborhood function is sMap.neigh and the
wolffd@0: %             neighborhood width is sMap.trainhist(end).radius_fin
wolffd@0: %             or 1 if this is not available 
wolffd@0: %  'fuzzy'    fuzzy response calculated by summing 
wolffd@0: %
wolffd@0: %                            1
wolffd@0: %                       ------------
wolffd@0: %                       1 +  (q/a)^2
wolffd@0: %
wolffd@0: %             for each data sample, where q is a vector containing
wolffd@0: %             distance from the data sample to each map unit and 
wolffd@0: %             a is average quantization error
wolffd@0: %
wolffd@0: % REQUIRED INPUT ARGUMENTS
wolffd@0: %
wolffd@0: %   sMap              The vectors from among which the BMUs are searched
wolffd@0: %                     for. These must not have any unknown components (NaNs).
wolffd@0: %            (struct) map struct
wolffd@0: %            (matrix) codebook matrix, size munits x dim
wolffd@0: %                     
wolffd@0: %   sData             The data vector(s) for which the BMUs are searched.
wolffd@0: %            (struct) data struct
wolffd@0: %            (matrix) data matrix, size dlen x dim
wolffd@0: %
wolffd@0: % OPTIONAL INPUT ARGUMENTS
wolffd@0: %
wolffd@0: %   mode     (string) The respond mode: 'crisp' (default), 'kernel'
wolffd@0: %                     or 'fuzzy'. 'kernel' can only be used if 
wolffd@0: %                     the first argument (sMap) is a map struct.                     
wolffd@0: % 
wolffd@0: % OUTPUT ARGUMENTS
wolffd@0: % 
wolffd@0: %   hits     (vector) The number of hits in each map unit.
wolffd@0: %
wolffd@0: % EXAMPLES
wolffd@0: %
wolffd@0: %  hits = som_hits(sM,D);
wolffd@0: %  hits = som_hits(sM,D,'kernel');
wolffd@0: %  hits = som_hits(sM,D,'fuzzy');
wolffd@0: %
wolffd@0: % SEE ALSO
wolffd@0: % 
wolffd@0: %  som_bmus      Find BMUs and quantization errors for a given data set.
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 161199
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% check arguments
wolffd@0: 
wolffd@0: error(nargchk(2, 3, nargin));  % check no. of input args is correct
wolffd@0: 
wolffd@0: if isstruct(sMap), 
wolffd@0:   switch sMap.type, 
wolffd@0:    case 'som_map', munits = prod(sMap.topol.msize);
wolffd@0:    case 'som_data', munits = size(sMap.data,1);
wolffd@0:    otherwise, 
wolffd@0:     error('Illegal struct for 1st argument.')
wolffd@0:   end
wolffd@0: else 
wolffd@0:   munits = size(sMap,1); 
wolffd@0: end
wolffd@0: hits = zeros(munits,1);
wolffd@0: 
wolffd@0: if nargin<3, mode = 'crisp'; end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% action
wolffd@0: 
wolffd@0: % calculate BMUs
wolffd@0: [bmus,qerrs] = som_bmus(sMap,sData,1);
wolffd@0: 
wolffd@0: switch mode, 
wolffd@0: case 'crisp',    
wolffd@0: 
wolffd@0:  % for each unit, check how many hits it got
wolffd@0:  for i=1:munits, hits(i) = sum(bmus == i); end
wolffd@0:    
wolffd@0: case 'kernel',
wolffd@0: 
wolffd@0:  % check that sMap really is a map 
wolffd@0:  if ~isstruct(sMap) & ~strcmp(sMap.type,'som_map'), 
wolffd@0:    error('Kernel mode can only be used for maps.');
wolffd@0:  end	   
wolffd@0: 
wolffd@0:  % calculate neighborhood kernel  
wolffd@0:  Ud = som_unit_dists(sMap.topol).^2;
wolffd@0:  sTrain = sMap.trainhist(end); 
wolffd@0:  if ~isempty(sTrain), 
wolffd@0:    rad = sTrain.radius_fin; 
wolffd@0:    if isempty(rad) | isnan(rad), rad = 1; end 
wolffd@0:  else 
wolffd@0:    rad = 1; 
wolffd@0:  end    
wolffd@0:  rad = rad^2;
wolffd@0:  if rad==0, rad = eps; end % to avoid divide-by-0 errors
wolffd@0:  switch sTrain.neigh, 
wolffd@0:   case 'bubble',   H = (Ud<=rad); 
wolffd@0:   case 'gaussian', H = exp(-Ud/(2*rad)); 
wolffd@0:   case 'cutgauss', H = exp(-Ud/(2*rad)) .* (Ud<=rad);
wolffd@0:   case 'ep',       H = (1-Ud/rad) .* (Ud<=rad);
wolffd@0:  end
wolffd@0:  
wolffd@0:  % weight hits with neighborhood kernel
wolffd@0:  hits = sum(H(bmus,:),1)';
wolffd@0:    
wolffd@0: case 'fuzzy',	
wolffd@0: 
wolffd@0:  % extract the two matrices (M, D) and the mask
wolffd@0:  mask = [];
wolffd@0:  if isstruct(sMap), 
wolffd@0:    if strcmp(sMap.type,'som_data'), M = sMap.data; 
wolffd@0:    else M = sMap.codebook; mask = sMap.mask;
wolffd@0:    end
wolffd@0:  else M = sMap;
wolffd@0:  end
wolffd@0:  if any(isnan(M(:))), 
wolffd@0:    error('Data in first argument must not have any NaNs.'); 
wolffd@0:  end
wolffd@0: 
wolffd@0:  if isstruct(sData), 
wolffd@0:    switch sData.type, 
wolffd@0:     case 'som_map', 
wolffd@0:      D = sData.codebook; 
wolffd@0:      if isempty(mask), mask = sData.mask; end
wolffd@0:     case 'som_data', D = sData.data;
wolffd@0:     otherwise, error('Illegal 2nd argument.');
wolffd@0:    end	
wolffd@0:  else D = sData;
wolffd@0:  end
wolffd@0:  [dlen dim] = size(D);   
wolffd@0:  if isempty(mask), mask = ones(dim,1); end
wolffd@0: 
wolffd@0:  % scaling factor   
wolffd@0:  a = mean(qerrs).^2;
wolffd@0:  
wolffd@0:  % calculate distances & bmus
wolffd@0:  % (this is better explained in som_batchtrain and som_bmus)
wolffd@0:  Known = ~isnan(D); D(find(~Known)) = 0; % unknown components  
wolffd@0:  blen = min(munits,dlen); % block size 
wolffd@0:  W1 = mask*ones(1,blen); W2 = ones(munits,1)*mask'; D = D'; Known = Known';
wolffd@0:  i0 = 0; 
wolffd@0:  while i0+1<=dlen,    
wolffd@0:    inds = [(i0+1):min(dlen,i0+blen)]; i0 = i0+blen; % indeces   
wolffd@0:    Dist = (M.^2)*(W1(:,1:length(inds)).*Known(:,inds)) ...
wolffd@0: 	  + W2*(D(:,inds).^2) ...
wolffd@0: 	  - 2*M*diag(mask)*D(:,inds); % squared distances
wolffd@0:    hits = hits + sum(1./(1+Dist/a),2);   
wolffd@0:  end  
wolffd@0:  
wolffd@0:  otherwise, 
wolffd@0:   error(['Unknown mode: ' mode]);
wolffd@0:   
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: