% SOM demo. Written by R. Murray-Smith
% 28th Feb 2000

Xsize = 10;
Ysize = 10;
if 1
   MapXCoords = randn(Xsize, Ysize);
   MapYCoords = randn(Xsize, Ysize);
   MapZCoords = randn(Xsize, Ysize);
else
   [MapXCoords, MapYCoords] = meshgrid([1:Xsize]/Xsize,[1:Ysize]/Ysize);
   MapXCoords = MapXCoords';
   MapYCoords = MapYCoords';
   MapZCoords = 0.01*randn(Xsize, Ysize);
end
N = 125;
%[x,y] = meshgrid(1:sqrt(N));
[X] = 5*rand(N,2); x = X(:,1);y =X(:,2);
%inputs = [x(:) y(:) sqrt(x(:).^2+y(:).^2)];%randn(N,2);
%[b,inputs] = bucky;
%[x,y,z] = cylinder(1+cos([0:0.5:2*pi]),50); inputs = [x(:) y(:) z(:)];
%load test1d3d;
load test3d;

N = length(inputs(:,1));
% normalise the inputs
inputs  = inputs - ones(N,1)*mean(inputs);
inputs  = inputs./(ones(N,1)*std(inputs));
Niters = 1000;
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;

% 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;
h = exp(-h);

for t=1:Niters
   stepsize = stepsize.*0.999;
   %calculate distance function for all nodes of Map
   Dx = zeros(N,Nnodes);
   Dy = zeros(N,Nnodes);
   Dz = 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);
      Dz(:,i) = inputs(:,3) - ones(N,1)*MapZCoords(i);
   end
   
   [dummy,nearest] = min(((Dx.^2+Dy.^2+Dz.^2)').^0.33333);
   for pos=1:N
      i = nearest(pos); % winner node for input 'pos'
      hi = h(i,:,:); % find weighting function
      MapXCoords = MapXCoords + reshape((hi(:).*Dx(pos,:)')/N,Xsize, Ysize);
      MapYCoords = MapYCoords + reshape((hi(:).*Dy(pos,:)')/N,Xsize, Ysize);
      MapZCoords = MapZCoords + reshape((hi(:).*Dz(pos,:)')/N,Xsize, Ysize);
   end
   if rem(t,1) == 0
      figure(1)
      %surf(MapXCoords,MapYCoords,MapZCoords)
      [x,y,z] = gplot3(MapConnections, [MapXCoords(:) MapYCoords(:) MapZCoords(:)]);
      plot3(x,y,z,'-o', inputs(:,1), inputs(:,2),inputs(:,3),'.');
      hold on
      plot3(inputs(:,1), inputs(:,2),inputs(:,3),'.');
      hold off
      pause(0.01);
   end
end