
% *********** sofm.m  **********  ver. 3.2 ***************
% Demonstration of Self-Organizing Feature Maps
% using  Kohonen's Algorithm in MATLAB  
%                                                       Enjoy.

clc;
disp('');
disp('');
disp('Formation of the Feature maps using Kohonen Algorithm ') 
disp('                                       26 June 2000 ')
disp('         Andrew P. Paplinski                        ')
disp('         Computer Science and Software Engineering  ')
disp('         Monash University, Melbourne, Australia    ')
disp('         Copyright (c) 1997  by A.P. Paplinski      ')
disp('         app@csse.monash.edu.au                     ')
disp('                                                    ')

czy = input('initialisation? Y/N [Y]: ','s'); 
      if isempty(czy), czy = 'y' ;  end

if (czy == 'y') | (czy == 'Y'),
   
  clear
  % Generation of the input training patterns.
  % A two-dimensional input pattern is stored as a complex number.
  % First, the form of the input domain is selected:
  
  indom = menu('Select the form of the input domain:',...
            'a rectangle', ...
            'a triangle', ...'
            'a circle', ...
            'a ring'  , ...
            'a cross' ,...
            'a letter A');
      if isempty(indom), indom = 2; end

  % Next, the dimensionality of the output space, l, is selected.
  % The output units ("neurons") can be arranged in a linear,
  % i.e. 1-D way, or in a rectangle, i.e., in a 2-D space.

  el = menu('Select the dimensionality of the output domain:',...
            '1-dimensional output domain', ...
            '2-dimensional output domain'); 
      if isempty(el), el = 1; end
  
  m1 = 12 ; m2 = 18;  % m1  by  m2  array of output units 
  if (el == 1), m1 = m1*m2 ; m2 = 1 ; end
  m = m1*m2 ;
  fprintf('The output lattice is %d by %d\n', m1, m2)
  mOK = input('would you like to change it? Y/N [N]: ','s'); 
  if isempty(mOK), mOK = 'n' ;  end
  if (mOK == 'y') | (mOK == 'Y')
    m = 1 ;  
    while ~((m1 > 1) & (m > 1) & (m < 4000)) 
      m1 = input('size of the output lattice: \n  m1 = ') ;
      if (el == 2)
       m2 = input('m2 = ') ;  
      end 
      m = m1*m2 ;
    end
  end
  fprintf('The output lattice is %d by %d\n', m1, m2) 
  
  % The position matrix V
  if el == 1
     V = (1:m1)' ;
  else
     [v1 v2] = meshgrid(1:m1, 1:m2); V = [v1(:) v2(:)];
  end
  
  % Creating input patterns  
  N = 20*m ;    % N is the number of input vectors
  X = rand(1, N)+j*rand(1, N) ;
  ix = 1:N;
  
  if (indom == 2),
    ix = find((imag(X)<=2*real(X))&(imag(X)<=2-2*real(X))) ;
    elseif (indom == 3),
      ix = find(abs(X-.5*(1+j))<= 0.5) ;
    elseif (indom == 4),
      ix = find((abs(X-.5*(1+j))<= 0.5) & (abs(X-.5*(1+j)) >= 0.3)) ;
    elseif (indom == 5),
      ix = find((imag(X)<(2/3)&imag(X)>(1/3))| ...
                   (real(X)<(2/3)&real(X)>(1/3))) ;
    elseif (indom == 6),
      ix = find((2.5*real(X)-imag(X)>0 & 2.5*real(X)-imag(X)<0.5) | ...
                (2.5*real(X)+imag(X)>2 & 2.5*real(X)+imag(X)<2.5) | ...
             (real(X)>0.2 & real(X)<0.8 & imag(X)>0.2 & imag(X)<0.4) );
  end
  X = X(ix); N = length(X);

  figure(1)
  clf reset,  hold off,  % resetting workspace

  plot(X, '.'), title('Input Distribution')
  drawnow

  % Initialisation of weights:
  % a pair of real weights (one complex weight) is associated with
  % each output unit. Therefore the weight matrix is an  m1 by m2
  % matrix of complex numbers. Initial values of the weights
  % are taken from the input distribution :
  
  W = X(1:m).' ;  X = X((m+1):N) ; N = N-m ;

  % as a check, the count of wins for each output unit is calculated
  % in the matrix  "hits". The expected value of  
  %                hits = N*#epochs/m  
  
  hits = zeros(m,1);
  
  % An Initial Feature Map: The map is represented by a grid spanned
  % by the weights of adjacent output units plotted in the input space.
  % In other words, the weights  W(j)  and  W(k) are connected
  % by a grid segment if  V(j)  and  V(k) are coordinates of the 
  % adjacent output units
  
  % Initial values of the training gain, eta, and the spread, sigma,
  % of the neighbourhood function  
  eta = 0.4 ;   % training gain
  sg2i = ((m1-1)^2+(m2-1)^2)/4 ; % sg2 = 2 sigma^2
  sg2 = sg2i ;

  figure(2)
  clf reset
  plot([0 1],[0 1],'.'), grid, hold on, 
  if el == 1
    plot(W, 'b'),  
  else
    FM = full(sparse(V(:,1), V(:,2), W)) ;
    plot(FM, 'b'),  plot(FM.', 'r') ;
  end
  title(['eta = ', num2str(eta,2), ...
        '  sigma^2 = ', num2str(sg2/2,3)])
  hold off, drawnow
  
  %  end of initialisation
else  % continuation
  eta = input('input the value of eta [0.4]: ') ;
     if isempty(eta), eta = 0.4; end
  sg2 = input(['input the value of 2sigma^2 [', ...
           num2str(sg2i), ']: ']) ;
     if isempty(sg2), sg2 = sg2i; end
end 

% rates of diminishing the values of eta and sigma during training
% We aim at final values:  eta = 0.1  and  sigma = 1 (sg2=2)
reta = (0.2)^(2/N); rsigma = (1/sg2)^(2/N) ;

% ********************** main loop ******************   

frm = 1;
for n = 1:N
      
  % for each input pattern, X(n),  and for each output unit 
  % which store the weight vector W(v1, v2), the distance  
  % between  X(n) and  W  is calculate
  
  WX = X(n) - W ; 
  
  % Coordinates of the winning neuron, V(kn, :), i.e.,
  % the neuron for which  abs(WX) attains minimum
  [mnm kn] = min(abs(WX));  vkn = V(kn, :) ;
 
  hits(kn) = hits(kn)+1; % utilization of neurons

  % The neighbourhood function, NB, of the "bell" shape,
  % is centered around the winning unit V(kn, :)
  rho2 = sum(((vkn(ones(m, 1), :) - V).^2), 2) ;
  NB = exp(-rho2/sg2) ;

  % Finally, the weights are updated according to the
  % Kohonen learning law:
  W = W + eta*NB.*WX ;

  % Values of  "eta", and "sigma" are reduced
  % if (n==floor(N/2)), eta = 0.8; end
  if (n<N/2), %ordering and convergence phase
    sg2 = sg2*rsigma; 
  else 
    eta = eta*reta; 
  end
 
  % Every  100  updates, the feature map is plotted
  if rem(n, 10) == 0
    plot([0 1],[0 1],'.'), grid, hold on, 
    if el == 1
      plot(W, 'b'), plot(W, '.r'),  
    else
      FM = full(sparse(V(:,1), V(:,2), W)) ;
      plot(FM, 'b'),  plot(FM.', 'r') ;
    end
    title(['eta = ',num2str(eta,2), ...
           '  sigma^2 = ', num2str(sg2/2,3), ...
           '   n = ', num2str(n)])
    hold off, drawnow
 end
  if sum(n==round([1, N/4, N/2, 3*N/4 N]))==1
    print('-depsc2', '-f2', ['Jsom2Dt', num2str(frm)])
    frm = frm+1;
  end
 % MV(frm) = getframe; frm = frm+1 ;
end

% Final presentation of the result
plot([0 1],[0 1],'.'), grid, hold on, 
if el == 1
  plot(W, 'b'), plot(W, '.r'),  
else
  FM = full(sparse(V(:,1), V(:,2), W)) ;
  plot(FM, 'b'),  plot(FM.', 'r') ;
end
title('A Feature Map'), hold off , drawnow
 
% ***************  end ************************

