% SOM demo. Written by R. Murray-Smith
% 28th Feb 2000
% note: as this program is for visualisation purposes it is restricted to 2D-input spaces.

Xsize = 13;
Ysize = 13;
if 0
    % initialise map as a grid
   [MapXCoords, MapYCoords] = meshgrid([1:Xsize]/Xsize,[1:Ysize]/Ysize);
   MapXCoords = MapXCoords';
   MapYCoords = MapYCoords';
else
   % initialise map randomly
   MapXCoords = randn(Xsize, Ysize);
   MapYCoords = randn(Xsize, Ysize);
end


% set up inputs
if 0
   N = 100;
   [x,y] = meshgrid(1:sqrt(N));
   inputs = [x(:) y(:)];%randn(N,2);
else
   load test2d;
   N = length(inputs(:,1));
end

% initialise Self-Organising map data structures
Niters = 2000;
Nnodes = Xsize*Ysize;
H = zeros(Nnodes,Xsize,Ysize);
% create adjacency matrix for plotting graph.
LeftRight = diag(ones(Xsize-1,1),1) ...
            + diag(ones(Xsize-1,1),-1);
UpDown    = diag(ones(Nnodes-Xsize,1),-Xsize)...
            + diag(ones(Nnodes-Xsize,1),Xsize);
MapConnections = kron(diag(ones(Ysize,1)),LeftRight)...
                 + UpDown;

% pre-calculate node distance function
for node=1:Nnodes
   [nodex,nodey] = ind2sub(size(MapXCoords),node);
   for i=1:Xsize
      for j = 1:Ysize
         H(node,i,j) =  max(abs([(i-nodex) (j-nodey)]'));
      end
   end
end
stepsize = 1;
scale = 1;

for t=1:Niters
   stepsize = stepsize.*0.999; % make smaller steps as iterations inc.
   scale = scale*0.999;		% shrink neighbourhoods as iterations increase
   h = exp(-H/scale);
   %calculate distance function for all nodes of Map
   Dx = zeros(N,Nnodes);
   Dy = zeros(N,Nnodes);
   for i=1:Nnodes
      Dx(:,i) = inputs(:,1) - ones(N,1)*MapXCoords(i);
      Dy(:,i) = inputs(:,2) - ones(N,1)*MapYCoords(i);
   end
   
   [dummy,nearest] = min(sqrt(Dx.^2+Dy.^2)');
   for pos=1:N
      i = nearest(pos); % winner node for input 'pos'
      %for j = 1:Nnodes
      %   [nodex,nodey] = ind2sub(size(MapXCoords),j);
      %   MapXCoords(j) = MapXCoords(j) + (h(i,nodex,nodey).*Dx(pos,j))/N;
      %   MapYCoords(j) = MapYCoords(j) + (h(i,nodex,nodey).*Dy(pos,j))/N;
      %end
      hi = h(i,:,:); % find weighting function
      MapXCoords = MapXCoords + reshape((hi(:).*Dx(pos,:)')/N,Xsize, Ysize);
      MapYCoords = MapYCoords + reshape((hi(:).*Dy(pos,:)')/N,Xsize, Ysize);
   end
   figure(1)
   plot(MapXCoords(:),MapYCoords(:),'o',inputs(:,1), inputs(:,2),'.');
   hold on
   [x,y] = gplot(MapConnections, [MapXCoords(:) MapYCoords(:)]);
   plot(x,y,'-o')
   hold off
   pause(0.01);
end