```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:
%  '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
%             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.

%%%%%%%%%%%%% 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
%             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');
%
%
%  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 arguments

error(nargchk(2, 3, nargin));  % check no. of input args is correct

if 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.')
end
else
munits = size(sMap,1);
end
hits = 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),
else
end
switch sTrain.neigh,
end

% weight hits with neighborhood kernel
hits = sum(H(bmus,:),1)';

case 'fuzzy',

% extract the two matrices (M, D) and the mask
if isstruct(sMap),
if strcmp(sMap.type,'som_data'), M = sMap.data;
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;
case 'som_data', D = sData.data;
otherwise, error('Illegal 2nd argument.');
end
else D = sData;
end
[dlen dim] = size(D);

% 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) ...