% somdemon.m: demonstrate SOM using a linear array of neurons
% with annimation
% copyright (c) 2001 by Yu Hen Hu
% created: 3/18/01 
%
% mfiles needed:  datagen1.m, randomize.m
%
clear
clf,figure(1)
% N=input('# of samples in each cluster = ');
disp('In the feature space, two Gaussian clusters are formed. The first cluster');
disp('locates at (.7 .7) and the second at (-.8, -.8). The standard deviation is 0.2.');
disp('There are 100 samples at each cluster');
N=100;
N2=N+N; B1=ceil(N/2); B2=N+B1;
eta=0.2;
means=[.7  -.8
       .7  -.8];
var= [0.2  0.2];
x=datagen1([N N],[means;var]); % x is 2N by 2
x=randomize(x);  % randomize the ordering of clusters
% disp('Choose initial conditions:');
% ncenter=input('# of clusters (neurons) to use: ');
disp('We use 11 neurons, forming a linear array, labeled 1 to 15');
disp('Initially, these neurons are randomly placed in the feature space.');
ncenter=11;
w=rand(ncenter,2)-0.5*ones(ncenter,2);  % neurons initially at random place
% but implicitly has labels as 1:ncenter

subplot(121),plot(x(:,1),x(:,2),'r.',w(:,1),w(:,2),'*-')
axis([-2 2 -2 2])
title('initial')
disp('Enter 1 to pause .1 second after each iteration, ');
pmode=input('Otherwise hit return key for continuous annimation: ');
if isempty(pmode),  pmode=0; end
i=1; iter=1; converge=0;
while converge==0,
   dn=ones(ncenter,1)*x(i,:)-w;
   ddn=sum((dn.*dn)')'; % ddn: ncenter by 1
   [tmp,istar]=min(ddn);
   if istar==1,
      w([istar:istar+1],:)=w([istar:istar+1],:)+eta*(ones(2,1)*x(i,:)-w([istar:istar+1],:));
   elseif istar==ncenter,
      w([istar-1:istar],:)=w([istar-1:istar],:)+eta*(ones(2,1)*x(i,:)-w([istar-1:istar],:));
   else
      w([istar-1:istar+1],:)=w([istar-1:istar+1],:)+eta*(ones(3,1)*x(i,:)-w([istar-1:istar+1],:));
   end
   subplot(122),
   plot(x(:,1),x(:,2),'r.',x(i,1),x(i,2),'o',w(:,1),w(:,2),'*-')
   drawnow
   if pmode==1,
      pause(.1)
   end
   i=rem(i+1,N2); iter=iter+1;
   if i==0,
      x=randomize(x);  % randomize the ordering of clusters
      i=1;             % start with a reordered input 
   end
   if rem(iter,50)==0,
      eta=eta*0.9;
      %converge=input('Enter return key to continue, 1 to quit! ');
      %if isempty(converge), converge=0; end
      if iter >=500, converge=1; end
   end
end
disp('after 500 iterations, you can see the chain of neurons are gradually');
disp('straightened out,(sometimes tangle remains) ');
disp('and gives a 1D approximation of the 2D distribution');
subplot(122),plot(x(:,1),x(:,2),'r.',w(:,1),w(:,2),'*-')
text(w(1,1),w(1,2)+0.2,'1'), text(w(ncenter,1),w(ncenter,2)+0.2,int2str(ncenter))
axis([-2 2 -2 2])
title(['at end of ' num2str(iter) ' iterations'])