// Kohonen network - topological mapping
function [nn]=SOM(a,b,n,m)
[x,n]=generare(a,b,n); // data generation
nn=tlist(["retea Kohonen","m", "W"],m,zeros(m,m,2));
nn.W=init(a,b,m); // initialization
disp(nn.W);
eta=1.25; s=m;
for t=1:100
  for l=1:n
  // find the winning unit
  [imin,jmin]=winner(nn,x(:,l));
  // weights adjustment
  for i=1:nn.m   
    for j=1:nn.m
     nn.W(i,j,1)=nn.W(i,j,1)+eta*(x(1,l)-nn.W(i,j,1))*lambda(i,j,imin,jmin,s);
     nn.W(i,j,2)=nn.W(i,j,2)+eta*(x(2,l)-nn.W(i,j,2))*lambda(i,j,imin,jmin,s);
    end 
  end  
  end
//  eta=eta/log(t+1);
  eta=0.995*eta;
  if modulo(t,10) ==0 & s>1
      s=s/1.25;
//      pause
  end    
//  disp(t,"t="); disp(eta,"eta=");
end

plot(x(1,:),x(2,:),'ko');
for i=1:nn.m do
 for k=1:nn.m-1 do
    plot([nn.W(i,k,1),nn.W(i,k+1,1)],[nn.W(i,k,2),nn.W(i,k+1,2)],'r*-');
end
end
for k=1:nn.m do
 for i=1:nn.m-1 do
    plot([nn.W(i,k,1),nn.W(i+1,k,1)],[nn.W(i,k,2),nn.W(i+1,k,2)],'r*-');
end
end
endfunction

function [imin,jmin]=winner(nn,x)
        imin=1; jmin=1;
        for i=1:nn.m
          for j=1:nn.m
            w1=[nn.W(i,j,1);nn.W(i,j,2)];
            w2=[nn.W(imin,jmin,1);nn.W(imin,jmin,2)];
            if dist(x,w1)<dist(x,w2)
               imin=i; jmin=j;
            end
          end   
        end 
endfunction    

function W=init(a,b,m)
for i=1:m
 for j=1:m
  for k=1:2
     W(i,j,k)=grand(1,1,'unf',a,b); 
  end
 end
end 
endfunction

function d=dist(x,y)
n=length(x);
d=0;
for i=1:n
    d=d+(x(i,1)-y(i,1))^2;
end    
endfunction

// neighborhood function
function f=lambda(i,j,imin,jmin,s)
   f=exp(-((i-imin)^2+(j-jmin)^2)/(2*s^2));
endfunction

function [x,n]=generare(a,b,n)
x=[]; 
for k=1:n
    x=[x [grand(1,1,"unf",a,b);grand(1,1,"unf",a,b)]];
  
end
//plot(x(1,:),x(2,:),'ko');
endfunction

SOM(0,1,20,5)
