% Neural Gas Script File

%Load in the data set in X (number of points by dimension)

load Data6;
X=XX{1};     %2 circles data
clear XX group_num

%X=rand(300,2);

%
% Initialization:
%

MaxIters=10000;      %Maximum number of iterations

MaxAge=100;         %Maximum age for edges between centers
MaxCenters=50;     %Maximum number of centers allowed
NumCenters=2;
Lambda=300;         %For Center Insertion
ew=0.05;            %Parameters for updating centers
en=0.0006;
alpha=0.5;          %Parameters for dealing with Error
beta=0.0005;

[numX,dimX]=size(X);
temp=randperm(numX);

temp=[105,143];     %Only for the double circle data

fprintf('Using Centers %d %d \n',temp(1),temp(2));

C=X(temp(1:NumCenters),:);
M=ones(NumCenters,NumCenters)-eye(NumCenters);
Age=zeros(NumCenters,NumCenters);
E=zeros(1,NumCenters);

for j=1:MaxIters
    
    %
    %  Pick a data point at random, find distances to centers
    %
    k=ceil(rand*numX);
    P=edm(X(k,:),C);
    [vals,idx]=sort(P);
    
    %
    %  Update the Error
    %
    
    w=idx(1); s=idx(2);
    E(w)=E(w)+ew*(vals(1)^2);
    
    %
    %  Update the centers
    %
    
    C(w,:)=C(w,:)+ew*(X(k,:)-C(w,:));
    idx=find(M(w,:)==1);
    if length(idx)>0
        C(idx,:)=C(idx,:)+en*(repmat(X(k,:),length(idx),1)-C(idx,:));
    end

    %
    % Age the existing connections
    %
    
    Age(w,idx)=Age(w,idx)+1;
    Age(idx,w)=Age(idx,w)+1;
    Age=Age-diag(diag(Age));
    
    %
    % Be sure there is a connection between the current winner and the
    % second closest center.
    %
    
    M(w,s)=1; M(s,w)=1;
    Age(w,s)=0; Age(s,w)=0;
    
    %
    %  If an edge is too old, disconnect the centers:
    %
    AA=Age(:);
    MM=M(:);
    idx=find(AA>MaxAge);
    if length(idx)>0
        AA(idx)=zeros(length(idx),1);
        Age=reshape(AA,NumCenters,NumCenters);
        MM(idx)=zeros(length(idx),1);
        M=reshape(MM,NumCenters,NumCenters);
    end
    clear AA MM
    
    %
    %  If a center is not connected to anything by an edge, remove it.
    %   - Corresponds to a zero row (or column)
    
    temp=sum(M,2);
    idx=find(temp==0);
    if length(idx>0)  %Remove empty centers
        fprintf('Removing empty centers\n');
        C(idx,:)=[];
        M(idx,:)=[];
        M(:,idx)=[];
        Age(idx,:)=[];
        Age(:,idx)=[];
        E(idx)=[];
    end
    
    %
    %  If j is a multiple of Lambda, and we have not reached the max number
    %  of centers, add a new center
    %
    if mod(j,Lambda)==0
    fprintf('Adding a new center\n');
        %
        %  Find the center with largest error
        %
        [vals,idx]=max(E);
        %
        %  Insert a new center as the midpoint of the winner and its
        %  neighbor
        %
        w=idx(1);
        idx=find(M(w,:)==1);
        [vals,idx2]=max(E(idx));
        w2=idx(idx2(1));
        
        %fprintf('Adding a center between centers %d and %d\n',w,w2);
        
        NumCenters=NumCenters+1;
        C(NumCenters,:)=0.5*(C(w,:)+C(w2,:));
        
        E(w)=alpha*E(w);
        E(w2)=alpha*E(w2);
        E(NumCenters)=alpha*E(w);
        
        %
        % Connect the new center in, reset ages, delete old connection
        %
        
        M(w,w2)=0; M(w2,w)=0;
        M=[M,zeros(NumCenters-1,1)];
        M=[M;zeros(1,NumCenters)];
        
        M(w,NumCenters)=1; M(NumCenters,w)=1;
        M(w2,NumCenters)=1; M(NumCenters,w2)=1;
        
        Age(w,w2)=0; Age(w2,w)=0;
        Age=[Age,zeros(NumCenters-1,1)];
        Age=[Age;zeros(1,NumCenters)];
    
    end
    
    %
    %  Decrease all errors by beta
    % 
    
    E=(1-beta)*E;
    
    %
    % Stopping Criteria here
    %
    
    NumCenters=length(E);
    if NumCenters>=MaxCenters
        fprintf('Maximum Number of Centers Reached\n');
        break
    end
    if mod(j,20)==0
      plotNG
      pause(0.1)
    end
end