wolffd@0: function [Neurons] = neural_gas(D,n,epochs,alpha0,lambda0)
wolffd@0: 
wolffd@0: %NEURAL_GAS Quantizes the data space using the neural gas algorithm.
wolffd@0: %
wolffd@0: % Neurons = neural_gas(D, n, epochs, [alpha0], [lambda0])
wolffd@0: %
wolffd@0: %   C = neural_gas(D,50,10);
wolffd@0: %   sM = som_map_struct(sD); 
wolffd@0: %   sM.codebook = neural_gas(sD,size(sM.codebook,1),10);
wolffd@0: %
wolffd@0: %  Input and output arguments ([]'s are optional):
wolffd@0: %   D          (matrix) the data matrix, size dlen x dim
wolffd@0: %              (struct) a data struct
wolffd@0: %   n          (scalar) the number of neurons
wolffd@0: %   epochs     (scalar) the number of training epochs (the number of
wolffd@0: %                       training steps is dlen*epochs)
wolffd@0: %   [alpha0]   (scalar) initial step size, 0.5 by default
wolffd@0: %   [lambda0]  (scalar) initial decay constant, n/2 by default
wolffd@0: %
wolffd@0: %   Neurons    (matrix) the neuron matrix, size n x dim
wolffd@0: %
wolffd@0: % See also SOM_MAKE, KMEANS.
wolffd@0: 
wolffd@0: % References: 
wolffd@0: %  T.M.Martinetz, S.G.Berkovich, and K.J.Schulten. "Neural-gas" network
wolffd@0: %  for vector quantization and its application to time-series prediction. 
wolffd@0: %  IEEE Transactions on Neural Networks, 4(4):558-569, 1993.
wolffd@0: 
wolffd@0: % Contributed to SOM Toolbox vs2, February 2nd, 2000 by Juha Vesanto
wolffd@0: % Copyright (c) by Juha Vesanto
wolffd@0: % http://www.cis.hut.fi/projects/somtoolbox/
wolffd@0: 
wolffd@0: % juuso 101297 020200
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Check arguments and initialize
wolffd@0: 
wolffd@0: error(nargchk(3, 5, nargin));  % check the number of input arguments
wolffd@0: 
wolffd@0: if isstruct(D), D = D.data; end
wolffd@0: [dlen,dim] = size(D);
wolffd@0: Neurons = (rand(n,dim)-0.5)*10e-5; % small initial values
wolffd@0: train_len = epochs*dlen;
wolffd@0: 
wolffd@0: if nargin<4 | isempty(alpha0) | isnan(alpha0), alpha0 = 0.5; end
wolffd@0: if nargin<5 | isempty(lambda0) | isnan(lambda0), lambda0 = n/2; end
wolffd@0: 
wolffd@0: % random sample order
wolffd@0: rand('state',sum(100*clock));
wolffd@0: sample_inds = ceil(dlen*rand(train_len,1));
wolffd@0: 
wolffd@0: % lambda
wolffd@0: lambda = lambda0 * (0.01/lambda0).^([0:(train_len-1)]/train_len);
wolffd@0: 
wolffd@0: % alpha
wolffd@0: alpha = alpha0 * (0.005/alpha0).^([0:(train_len-1)]/train_len);
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wolffd@0: %% Action
wolffd@0: 
wolffd@0: for i=1:train_len,
wolffd@0: 
wolffd@0:   % sample vector
wolffd@0:   x = D(sample_inds(i),:); % sample vector
wolffd@0:   known = ~isnan(x);       % its known components
wolffd@0:   X = x(ones(n,1),known);  % we'll need this 
wolffd@0: 
wolffd@0:   % neighborhood ranking
wolffd@0:   Dx = Neurons(:,known) - X;  % difference between vector and all map units
wolffd@0:   [qerrs, inds] = sort((Dx.^2)*known'); % 1-BMU, 2-BMU, etc.
wolffd@0:   ranking(inds) = [0:(n-1)];             
wolffd@0:   h = exp(-ranking/lambda(i));
wolffd@0:   H = h(ones(length(known),1),:)';
wolffd@0: 
wolffd@0:   % update 
wolffd@0:   Neurons = Neurons + alpha(i)*H.*(x(ones(n,1),known) - Neurons(:,known));
wolffd@0: 
wolffd@0:   % track
wolffd@0:   fprintf(1,'%d / %d \r',i,train_len);
wolffd@0:   if 0 & mod(i,50) == 0, 
wolffd@0:     hold off, plot3(D(:,1),D(:,2),D(:,3),'bo')
wolffd@0:     hold on, plot3(Neurons(:,1),Neurons(:,2),Neurons(:,3),'r+')
wolffd@0:     drawnow
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: fprintf(1,'\n');
wolffd@0: 
wolffd@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%