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