function C=NeuralGas01(X,C)
%FUNCTION C=NeuralGas(X,C)
%   This is the Neural Gas clustering algorithm.  There are
%   two ways to call this program:
%     [C,M]=NeuralGas(X);
%        X is dimension by (number of points)
%        With only one input argument, the program will
%        use the function initng.m for the initialization
%        procedures, and train the network.  
%
%     [C,M]=NeuralGas(X,C)
%        With two input parameters, we assume the centers
%        have been initialized in the structure C.  See the
%        initng file for structure specifications.

%We will use the following as the current update parameters:
%  lr = current learning rate
%  nbr= current neighborhood size (lambda)
%  t  = current time setting (for connection matrix)

%INITIALIZATION:

if nargin==1  %use the initng function if only 1 input
C=initng(X);
end

[m,n]=size(X)  %Dimension is m, number of points is n
[m,numcenters]=size(C.cen);

lr=C.lr(1);
nbr=C.nbr(1);
if C.tflag
   t=C.tcon(1);
   Age=zeros(numcenters,numcenters);
end

%MAIN LOOP:

for j=1:C.epochs
    for k=1:C.tmax
        
if rem(k,1000)==0
  disp('iterate =');
  disp(k)
  disp('out of')
  disp(C.tmax)
  
  %Graphical updates:
  clf
  plotng01(C);
  hold on
  plot(X(1,:),X(2,:),'k.');
  plot(C.cen(1,:),C.cen(2,:),'r*');
  hold off
  drawnow;
end


%****************************************************
%
%   Step 1:  Choose a data point at random and compute
%              distance vector.
%
%****************************************************

curidx=ceil(n*rand);
Nx=X(:,curidx);
D=(C.cen - repmat(Nx,1,numcenters)).^2;
dd=sum(D);
[Md,w]=min(dd);

D=(C.cen - repmat(C.cen(:,w),1,numcenters)).^2;
dd=sum(D);
[Md,Id]=sort(dd);

%********************************************************
%
%  Step 2:  Update all centers and training parameters
%
%*********************************************************

for s=1:numcenters
aa=lr*exp(-(s-1)/nbr);
C.cen(:,Id(s))=C.cen(:,Id(s))+aa*(X(:,curidx)-C.cen(:,Id(s)));
end

lr=paramUpdate(C.lr,k,C.tmax);
nbr=paramUpdate(C.nbr,k,C.tmax);

%************************************************************
%
%  Step 3:  If necessary, update Connection matrix and Time.
%
%************************************************************
if C.tflag
   ab=Id(2);
   C.M(w,ab)=1;
   Age(w,ab)=0;
   Ix=find(C.M>0);
   Age(Ix)=Age(Ix)+1;
   Iy=find(Age>=t);
   C.M(Iy)=0;

   t=paramUpdate(C.tcon,k,C.tmax);
   
end  %End of time and connection update

    end %End of C.tmax (k) loop
end %End of C.epochs (j) loop