% Self-organizing feature map (SOM al la Kohonen) in two dimensions
clear; clf; 
n_out= 20; % number of nodes in the output layer 
k=.5;    % learning rate
sig=2;   % (initial) width of `learning window'
sig_inv=1/(2*sig^2); % this is actually a factor we need below
[X,Y] = meshgrid(1:n_out, 1:n_out); % X is a matrix that contains the 'x' index of a node 
                                    % Y is a matrix that contains the 'y' index of a node 
feature_minmax = [0 2]; % limits on of the feature domain for plot purposes
                                    
                                    
r=zeros(n_out); % initial activity (rate) of the nodes 

% 2 input nodes: w1 is the weight vector from the first input node to all SOM nodes 
%                w2 is the weight vector from the second input node to all SOM nodes 
% Initial random weights:
w1=0.5-.1*(2*rand(n_out)-1);
w2=0.5-.1*(2*rand(n_out)-1);

% Initial orderly weights
%for i=1:n_out; for j=1:n_out; w1(i,j)=.4+i/n_out/10; w2(i,j)=.5+j/n_out/10; end; end;   

% Plot grid of receptive fields (weight values for eacxh node)
hold on;
plot(w1,w2,'k')
plot(w1',w2','k')
%axis([feature_minmax feature_minmax]); xlabel('w1'); ylabel('w2');
tstring=['Training examples:' INT2Str(0)]; title(tstring);
waitforbuttonpress

for epochs=1:100;
    %sig_inv=1/(2*sig^2);
    %if sig<1.5; sig=1.5; end % or any other decreasing function
    for i_example=1:100
        %k=0.1*i_example.^(-0.1);  % dreacresing learning rate 
        % training example
        if epochs <= 100;
            r_in=[rand;rand];
        else
            r_in=[rand;rand]+floor(2*rand)*[1;1];
        end

        % calculate winner: Node that is most excited (receptive filed fits
        % best weight!)
        r=exp(-(w1-r_in(1)).^2-(w2-r_in(2)).^2);
        [rmax,x_winner]=max(max(r));
        [rmax,y_winner]=max(max(r'));

        % update weight vector (train winner the most and neigbours
        % according to Gaussian distance
        r=exp(-((X-x_winner).^2+(Y-y_winner).^2)*sig_inv);
        w1=w1+k*r.*(r_in(1)-w1);
        w2=w2+k*r.*(r_in(2)-w2);
    end

% plot organization after 100 additional training examples = 1 additional epoch
clf; hold on;
plot(w1,w2,'k')
plot(w1',w2','k')
axis([feature_minmax feature_minmax]); xlabel('w1'); ylabel('w2');
tstring=['Training examples:' INT2Str(epochs*100)]; title(tstring);
disp([epochs,sig])
waitforbuttonpress

end