function Ne = som_neighborhood(Ne1,n)
%SOM_NEIGHBORHOOD Calculate neighborhood matrix.
%
% Ne = som_neighborhood(Ne1,n)
%
% Ne = som_neighborhood(Ne1);
% Ne = som_neighborhood(som_unit_neighs(topol),2);
%
% Input and output arguments ([]'s are optional):
% Ne1 (matrix, size [munits m]) a sparse matrix indicating
% the units in 1-neighborhood for each map unit
% [n] (scalar) maximum neighborhood which is calculated, default=Inf
%
% Ne (matrix, size [munits munits]) neighborhood matrix,
% each row (and column) contains neighborhood
% values from the specific map unit to all other
% map units, or Inf if the value is unknown.
%
% For more help, try 'type som_neighborhood' or check out online documentation.
% See also SOM_UNIT_NEIGHS, SOM_UNIT_DISTS, SOM_UNIT_COORDS, SOM_CONNECTION.
%%%%%%%%%%%%% DETAILED DESCRIPTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% som_neighborhood
%
% PURPOSE
%
% Calculate to which neighborhood each map unit belongs to relative to
% each other map unit, given the units in 1-neighborhood of each unit.
%
% SYNTAX
%
% Ne = som_neighborhood(Ne1);
% Ne = som_neighborhood(Ne1,n);
%
% DESCRIPTION
%
% For each map unit, finds the minimum neighborhood to which it belongs
% to relative to each other map unit. Or, equivalently, for each map
% unit, finds which units form its k-neighborhood, where k goes from
% 0 to n.
%
% The neighborhood is calculated iteratively using the reflexivity of
% neighborhood.
% let N1i be the 1-neighborhood set a unit i
% and let N11i be the set of units in the 1-neighborhood of any unit j in N1i
% then N2i (the 2-neighborhood set of unit i) is N11i \ N1i
%
% Consider, for example, the case of a 5x5 map. The neighborhood in case of
% 'rect' and 'hexa' lattices (and 'sheet' shape) for the unit at the
% center of the map are depicted below:
%
% 'rect' lattice 'hexa' lattice
% -------------- --------------
% 4 3 2 3 4 3 2 2 2 3
% 3 2 1 2 3 2 1 1 2 3
% 2 1 0 1 2 2 1 0 1 2
% 3 2 1 2 3 2 1 1 2 3
% 4 3 2 3 4 3 2 2 2 3
%
% Because the iterative procedure is rather slow, the neighborhoods
% are calculated upto given maximal value. The uncalculated values
% in the returned matrix are Inf:s.
%
% REQUIRED INPUT ARGUMENTS
%
% Ne1 (matrix) Each row contains 1, if the corresponding unit is adjacent
% for that map unit, 0 otherwise. This can be calculated
% using SOM_UNIT_NEIGHS. The matrix can be sparse.
% Size munits x munits.
%
% OPTIONAL INPUT ARGUMENTS
%
% n (scalar) Maximal neighborhood value which is calculated,
% Inf by default (all neighborhoods).
%
% OUTPUT ARGUMENTS
%
% Ne (matrix) neighborhood values for each map unit, size is
% [munits, munits]. The matrix contains the minimum
% neighborhood of unit i, to which unit j belongs,
% or Inf, if the neighborhood was bigger than n.
%
% EXAMPLES
%
% Ne = som_neighborhood(Ne1,1); % upto 1-neighborhood
% Ne = som_neighborhood(Ne1,Inf); % all neighborhoods
% Ne = som_neighborhood(som_unit_neighs(topol),4);
%
% SEE ALSO
%
% som_unit_neighs Calculate units in 1-neighborhood for each map unit.
% som_unit_coords Calculate grid coordinates.
% som_unit_dists Calculate interunit distances.
% som_connection Connection matrix.
% Copyright (c) 1999-2000 by the SOM toolbox programming team.
% http://www.cis.hut.fi/projects/somtoolbox/
% Version 1.0beta juuso 141097
% Version 2.0beta juuso 101199
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Check arguments
error(nargchk(1, 2, nargin));
if nargin<2, n=Inf; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Action
% initialize
if issparse(Ne1), Ne = full(Ne1); else Ne = Ne1; end
clear Ne1
[munits dummy] = size(Ne);
Ne(find(Ne==0)) = NaN;
for i=1:munits, Ne(i,i)=0; end
% Calculate neighborhood distance for each unit using reflexsivity
% of neighborhood:
% let N1i be the 1-neighborhood set a unit i
% then N2i is the union of all map units, belonging to the
% 1-neighborhood of any unit j in N1i, not already in N1i
k=1;
if n>1,
fprintf(1,'Calculating neighborhood: 1 ');
N1 = Ne;
N1(find(N1~=1)) = 0;
end
while k1, fprintf(1,'\n'); end
% finally replace all uncalculated distance values with Inf
Ne(find(isnan(Ne))) = Inf;
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% faster version?
l = size(Ne1,1); Ne1([0:l-1]*(l+1)+1) = 1; Ne = full(Ne1); M0 = Ne1; k = 2;
while any(Ne(:)==0), M1=(M0*Ne1>0); Ne(find(M1-M0))=k; M0=M1; k=k+1; end
Ne([0:l-1]*(l+1)+1) = 0;