%       function RP = NGAS(DP, NCLU, TMAX, LI, LF, EI, EF)
%
% Neural-Gas clustering.
%
% DP 		Data Point matrix. Every row vector represents a data point.
% NCLU		Number of clusters.
% TMAX		Number of adaptation steps to be performed.
% LI, LF	Initial and final value of Lambda.
% EI, EF	Initial and final value of Epsilon.

% Neural-Gas is a vector quantization procedure developed by T. Martinetz.
% This version was implemented by
% Stefano Ferrari 	ferrari@inb.mi.cnr.it
% Istituto di Neuroscienze e Bioimmagini - CNR
% via Mario Bianco, 9 - 20131 Milano ITALY

function RP = ngas(DP, NCLU, Tmax, Li, Lf, Ei, Ef)

if nargin < 7,
	help ngas;
	break;
end;
NCLU = round(NCLU);
if NCLU<1, 
	disp('NGAS: NCLU must be greater than or equal to 1');
	break;
end;
ns = size(DP, 1);
DIM = size(DP, 2);
SeedTemp = rand('seed');
rand('seed', 0);

UNI = ones(NCLU, 1);
UNI2 = ones(1, DIM);
mincoord = min(DP);
maxcoord = max(DP);
RP = (rand(NCLU, DIM).*(UNI*(maxcoord-mincoord)))+(UNI*mincoord);
disp('***************************************************');
disp(' ');
disp('                Neural-Gas algorithm');
disp('           Implemented by Stefano Ferrari');
disp('      Istituto di Neuroscienze e Bioimmagini');
disp('                CNR - Milano - Italy');
disp(' ');

LAMBDA=Li;
EPSIL=Ei;
rr = ones(1, NCLU);
dtpt = ones(1, DIM);
Dist = ones(size(RP));
Distvet = ones(1, NCLU);
Hl = ones(1, NCLU);

time1 = cputime;
for T = 0:Tmax, 
  dtpt=DP(rand*(ns-1)+1, :);                   

  Dist = UNI*dtpt-RP;
  Distvet = sum((Dist.^2)'); 	
  [ORD,IND]=sort(Distvet);
  rr(IND) = 0:NCLU-1;
  Hl=exp(-rr/LAMBDA);
  RP = RP+(EPSIL*Hl'*UNI2).* Dist;

  LAMBDA=Li*(Lf/Li).^(T/Tmax); 
  EPSIL=Ei*(Ef/Ei).^(T/Tmax);
  
  if ~rem(T, round(Tmax/48)),
	fprintf('*');
  end;
  
end;
time2 = cputime-time1;
fprintf('\n\n');
disp(time2);
rand('seed', SeedTemp);


