function som
%
% NSCI 850 - Computational Approaches to Neuroscience
%
% runs Kohonen's self-organizing feature map algorithm
% click on the plot to display results every 100 iterations
%
% G. Blohm (2010)
% adapted from Trappenberg (2010)
%-------------------------------------------------------------------------

clear;

%% initialize
nn = 10; % number of units along each dimension
eps = .1; % learning rate
sig = 2; % SD of Gaussian centroids
sig2 = 1/(2*sig^2); % this is just convenient for computations...
[X,Y] = meshgrid(1:nn,1:nn); 
ntrial = 0;

% Initial centres of prefered features:
c1 = 0.5-.3*(2*rand(nn)-1);
c2 = 0.5-.3*(2*rand(nn)-1);

%% training session
while(true)
    if(mod(ntrial,100) == 0) % Plot grid of feature centres
        clf; hold on; 
        axis square; 
        axis([0 1 0 1]);
        plot(c1,c2,'k'); 
        plot(c1',c2','k');
        xlabel('c_{ij1}')
        ylabel('c_{ij2}')
        tstring = ['After ' int2str(ntrial) ' training steps']; 
        title(tstring);
        waitforbuttonpress;
    end
    r_in = [rand;rand]; % new random input
    r = exp(-(c1-r_in(1)).^2-(c2-r_in(2)).^2); % activation of cortical sheet before winner-take-all
    [rmax,x_winner] = max(max(r)); % compute winner
    [rmax,y_winner] = max(max(r'));
    r = exp(-((X-x_winner).^2+(Y-y_winner).^2)*sig2); % final activiation of cortical sheet around winning node
    c1 = c1+eps*r.*(r_in(1)-c1); % input connections weight updates
    c2 = c2+eps*r.*(r_in(2)-c2);
    ntrial = ntrial+1;
end