% kohonen.m - simulates a self-organizing map

figure(2)
clf
figure(1)
clf

% image
imsz=10;
im=zeros(imsz);
[imx imy]=meshgrid(1:imsz,1:imsz); % coordinates in image
showim=0;

% neuron array
SZ=10;
[X Y]=meshgrid(1:SZ,1:SZ);
sigma=1.0;

% weights
W=rand(imsz^2,SZ^2);
W=W*diag(1./sqrt(sum(W.^2)));

% learning rate
eta=0.05;

subplot(121)
h=imagesc(im,'EraseMode','none',[0 1]);
axis image

t=0;
while (1)
  
  % paint a Gaussian blob at random position in image
  % (you will need to change the procedure for setting x and y
  % to create a scotoma or over-stimulate retina)
  x=ceil(imsz*rand);
  y=ceil(imsz*rand);
  im=exp(-0.5*((x-imx).^2+(y-imy).^2));
  set(h,'CData',im)
  drawnow
  
  % compute output and do wta
  y=W'*im(:);
  [v ind]=max(y(:));
  i=mod(ind-1,SZ)+1;
  j=floor((ind-1)/SZ)+1;

  % spread activation to neighbors
  R=exp(-0.5*((Y-i).^2+(X-j).^2)/sigma^2);

  % Hebbian weight update
  W=W+eta*im(:)*R(:)';

  % display network every 100 trials
  if mod(t,100)==0
    subplot(122)
    mux=reshape((abs(W)'*imx(:))./sum(abs(W))',SZ,SZ);
    muy=reshape((abs(W)'*imy(:))./sum(abs(W))',SZ,SZ);
    plot(mux,muy,'k'), hold on
    plot(mux',muy','k'), hold off
    axis ij, axis([1 SZ 1 SZ])
    drawnow
    W=W*diag(1./sqrt(sum(W.^2))); % normalize weight vectors
    figure(2)
    showrfs(W)
  end
  figure(1)

%   pause(0.1)
  t=t+1;
end
