wolffd@0: function H = som_neighf(sMap,radius,neigh,ntype)
wolffd@0: 
wolffd@0: %SOM_NEIGHF Return neighborhood function values.
wolffd@0: %
wolffd@0: % H = som_neighf(sMap,[radius],[neigh],[ntype]);
wolffd@0: % 
wolffd@0: %  Input and output arguments ([]'s are optional): 
wolffd@0: %   sMap     (struct) map or topology struct
wolffd@0: %   [radius] (scalar) neighborhood radius (by default, the last used value
wolffd@0: %                     in sMap.trainhist is used, or 1 if that is unavailable)
wolffd@0: %   [neigh]  (string) neighborhood function type (by default, ..., or 
wolffd@0: %                     'gaussian' if that is unavailable)
wolffd@0: %   [ntype]  (string) 'normal' (default), 'probability' or 'mirror'
wolffd@0: %
wolffd@0: %   H        (matrix) [munits x munits] neighborhood function values from 
wolffd@0: %                     each map unit to each other map unit
wolffd@0: %
wolffd@0: % For more help, try 'type som_batchtrain' or check out online documentation.
wolffd@0: % See also  SOM_MAKE, SOM_SEQTRAIN, SOM_TRAIN_STRUCT.
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: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Check arguments 
wolffd@0: 
wolffd@0: % defaults
wolffd@0: rdefault = 1; 
wolffd@0: ndefault = 'gaussian';
wolffd@0: tdefault = 'normal';
wolffd@0: 
wolffd@0: % map 
wolffd@0: switch sMap.type, 
wolffd@0:  case 'som_map',   
wolffd@0:    sTopol = sMap.topol; 
wolffd@0:    sTrain = sMap.trainhist(end); 
wolffd@0:    if isempty(sTrain.radius_fin) | isnan(sTrain.radius_fin), 
wolffd@0:      rdefault = 1; 
wolffd@0:    else
wolffd@0:      rdefault = sTrain.radius_fin;
wolffd@0:    end
wolffd@0:    if ~isempty(sTrain.neigh) & ~isnan(sTrain.neigh), 
wolffd@0:      ndefault = sTrain.neigh;
wolffd@0:    end
wolffd@0:  case 'som_topol', sTopol = sMap; 
wolffd@0: end
wolffd@0: munits = prod(sTopol.msize); 
wolffd@0: 
wolffd@0: % other parameters
wolffd@0: if nargin<2 | isempty(radius), radius = rdefault; end
wolffd@0: if nargin<3 | isempty(neigh), neigh = ndefault; end
wolffd@0: if nargin<4 | isempty(ntype), ntype = tdefault; end
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% initialize
wolffd@0:                                    
wolffd@0: % basic neighborhood 
wolffd@0: Ud = som_unit_dists(sTopol);
wolffd@0: Ud = Ud.^2;
wolffd@0: radius = radius.^2;
wolffd@0: if radius==0, radius = eps; end % zero neighborhood radius may cause div-by-zero error
wolffd@0: 
wolffd@0: switch ntype, 
wolffd@0: case 'normal', 
wolffd@0:     H = neighf(neigh,Ud,radius); 
wolffd@0: case 'probability', 
wolffd@0:     H = neighf(neigh,Ud,radius); 
wolffd@0:     for i=1:munits, H(i,:) = H(i,:)/sum(H(i,:)); end
wolffd@0: case 'mirror', % only works for 2-dim grid!!!
wolffd@0:     H  = zeros(munits,munits);
wolffd@0:     Co = som_unit_coords(sTopol); 
wolffd@0:     for i=-1:1,
wolffd@0:         for j=-1:1,
wolffd@0:            Ud = gridmirrordist(Co,i,j);
wolffd@0:            H  = H + neighf(neigh,Ud,radius); 
wolffd@0:         end
wolffd@0:     end    
wolffd@0: end     
wolffd@0: 
wolffd@0: return;
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% subfunctions
wolffd@0: 
wolffd@0: function H = neighf(neigh,Ud,radius)
wolffd@0: 
wolffd@0:  switch neigh, 
wolffd@0:  case 'bubble',   H = (Ud<=radius); 
wolffd@0:  case 'gaussian', H = exp(-Ud/(2*radius)); 
wolffd@0:  case 'cutgauss', H = exp(-Ud/(2*radius)) .* (Ud<=radius);
wolffd@0:  case 'ep',       H = (1-Ud/radius) .* (Ud<=radius);
wolffd@0:  end  
wolffd@0:  return; 
wolffd@0: 
wolffd@0: function Ud = gridmirrordist(Co,mirrorx,mirrory)
wolffd@0: 
wolffd@0:  [munits,mdim] = size(Co); 
wolffd@0:  if mdim>2, error('Mirrored neighborhood only works for 2-dim map grids.'); end
wolffd@0:    
wolffd@0:  % width and height of the grid
wolffd@0:  dx = max(Co(:,1))-min(Co(:,1));
wolffd@0:  dy = max(Co(:,2))-min(Co(:,2));
wolffd@0: 
wolffd@0:  % calculate distance from each location to each other location
wolffd@0:  Ud = zeros(munits,munits);
wolffd@0:  for i=1:munits, 
wolffd@0:    inds = [i:munits]; 
wolffd@0:    coi = Co(i,:); % take hexagonal shift into account
wolffd@0:    coi(1) = coi(1)*(1-2*(mirrorx~=0)) + 2*dx*(mirrorx==1); % +mirrorx * step
wolffd@0:    coi(2) = coi(2)*(1-2*(mirrory~=0)) + 2*dy*(mirrory==1); % +mirrory * step
wolffd@0:    Dco = (Co(inds,:) - coi(ones(munits-i+1,1),:))'; 
wolffd@0:    Ud(i,inds) = sqrt(sum(Dco.^2));
wolffd@0:    Ud(inds,i) = Ud(i,inds)'; 
wolffd@0:  end
wolffd@0:  return; 
wolffd@0: 
wolffd@0: