function W = som(a,k,nrruns,eta,h)
%SOM Simple routine computing a Self-Organizing Map (SOM)
%
%           W =  SOM(X,K)
%
% Train a 2D SOM on dataset X. In K the size of the map is defined. The
% map can maximally be 2D. When K contains just a single value, it is
% assumed that a 1D map should be trained. The output of the mapping
% contains the (negative) distances to all neurons. To obtain the index
% of the closest neuron, do A*W*LABELD.
%
%           W =  SOM(X,K,NRRUNS,ETA,H)
%
% Train the SOM for NRRUNS iterations, using learning rate ETA and a
% Gaussian neighborhood function with width H*sqrt(MAXD), where MAXD
% is the maximum distance in the dataset X.
%
% There is the extra feature, that NRRUNS, ETA and H can be vectors,
% such that it can be run several iterations using larger ETA and H,
% and after that with smaller values.
%
% Default: K=[5 5], NRRUNS = [20 40 40], ETA = [0.5 0.3 0.1],
% H = [0.6 0.2 0.01];
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% PCAM, PRKMEANS, PRPLOTSOM, PREX_SOM

% Copyright: D. Tax, davidt@ph.tn.tudelft.nl
% Faculty of Applied Physics, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

if nargin <5 | isempty(h)
	h = [0.6 0.2 0.01];
end
if nargin <4 | isempty(eta)
	eta = [0.5 0.3 0.1];
end
if nargin <3 | isempty(nrruns)
	nrruns = [20 40 40];
end
if nargin < 2| isempty(k),
	k = [5 5];
end
if nargin < 1 | isempty(a) 
    W = prmapping(mfilename,{k,nrruns,eta,h});
    W = setname(W,'Self-organising Map');
    return
end

if ~ismapping(k)           %training

	a = testdatasize(a);
	x = +a;     % use all the data!
	[nrx,dim] = size(x);
	% Set up the map size:
	if length(k)<2 % add a dummy dimension
		k = [k; 1];
	end

	% Some magic parameters:
	% The training consist of several runs, each differ in the width of
	% the neighborhood width, the learning rate and the number of SOM
	% updates:
	Dmax = distm(x,x); Dmax = max(Dmax(:));
	h = h*sqrt(Dmax);
	% (this should probably be something which can be set by the user)

	% Initialize the map
	w = 0.2*rand(k(1)*k(2),dim)+repmat(+mean(x),k(1)*k(2),1);

	% Run the different types of updates:
	for r1=1:length(nrruns)
		
		% Set up the neighbourhood:
		W = zeros(k(1),k(2));
		for i=1:k(1)
			for j=1:k(2)
				W(i,j) = exp(-((i-1)^2+(j-1)^2)/h(r1));
			end
		end

		% Do the real SOM training:
		for r2 = 1:nrruns(r1)

			%fprintf('%d/%d\t',r2,nrruns(r1));
			% Go randomly along the objects:
			I = randperm(nrx);
			for i=1:nrx
				% Find the winner:
				[mD,mI] = min(distm(w,x(I(i),:)),[],1);
				wx = rem(mI-1,k(1))+1;
				wy = floor((mI-1)/k(1))+1;
				% Update all weights:
				for j=1:k(1)*k(2)
					jx = rem(j-1,k(1))+1;
					jy = floor((j-1)/k(1))+1;
					w(j,:) = w(j,:) + eta(r1)*W(abs(jx-wx)+1,abs(jy-wy)+1)*...
						(x(I(i),:)-w(j,:));
				end
			end
		end  
	end

	clear W;
	% And save all useful data:
	V.k = k;  %(only used for plotting later...)
	V.neurons = w;
	W = prmapping(mfilename,'trained',V,(1:prod(k))',dim,prod(k));
	W = setname(W,'Self-organising Map');
else
    
	W = getdata(k); %unpack
	m = size(a,1); 
	% compute the distance to the nearest neuron in the map:
	%[mD,out] = min(distm(+a,W.w),[],2);
	D = distm(+a,W.neurons); % only compute the distances
	
	% Output a dataset with the negative distances:
	W = setdat(a,-D,k);
end