
% *********** somko.m  **********  ver. 3.2 ***************
% Demonstração dos mapas característicos self-Organizing   usando
%o algoritmo de Kohonen em MATLAB


clc;
disp('');
disp('');
disp('Formação de mapas característicos usando algoritmo de Kohonen  ') 
disp('                                       Setembro 2000 ')
disp('         Maritza Egoavil Churampi                       ')
disp('         Redes Neurais  ')
disp('         Posgraduação em Ciencias da Computação  ')
disp('         UFSC    ')

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

if (czy == 'y') | (czy == 'Y'),
   
  clear
  % Geração dos testes padrões do treinamento da entrada
  % Um teste padrão bidimensional da entrada é armazenado como um número complexo.
  
  % Primeiro, a forma do dominio entrada é selecionado:
  
  indom = menu('Selecione a forma do dominio entrada:',...
            ' rectangulo', ...
            ' triangulo', ...'
            ' circulo', ...
            ' anel'  , ...
            ' cruz' ,...
            ' letra A');
         if isempty(indom), indom = 2; end
        % Em seguida, o dimensão do espaço da saída, l, é selecionado.
  %As unidades de saída ("neurônios") podem ser arranjadas em um linear. 
 % p.e. de uma 1-D, ou em um retângulo, p.e., em um 2-D espaço.

  
  el = menu('Seleciona a dimensão do dominio saída :',...
            'dominio de saída de 1-dimensão', ...
            'dominio de saída de 2-dimensão'); 
      if isempty(el), el = 1; end
  
  m1 = 13 ; m2 = 17;  % m1  pela  m2  unidades de saída 
  if (el == 1), m1 = m1*m2 ; m2 = 1 ; end
  m = m1*m2 ;
  fprintf('A  lattice de saída é %d por %d\n', m1, m2)
  mOK = input('você gosta de mudar? 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('tamanho de uma lattice de saída: \n  m1 = ') ;
      if (el == 2)
       m2 = input('m2 = ') ;  
      end 
      m = m1*m2 ;
    end
  end
  fprintf('A lattice de saída é  %d by %d\n', m1, m2) 
  
  % A matriz posição V
  if el == 1
     V = (1:m1)' ;
  else
     [v1 v2] = meshgrid(1:m1, 1:m2); V = [v1(:) v2(:)];
  end
  
  % Creando teste padrão da entrada   
  N = 20*m ;    % N é numero de vectores entrada.
  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, %  restaura o espaço de trabalho

  plot(X, '.'), title(' Distribuição da  Entrada ')
  drawnow

  % Inicio dos  pesos:
% um o par de pesos reais (um peso complexo) é associado com  
%cada unidade de saída  Conseqüentemente a matriz do peso é um m1
%por  da matriz do m2 de números complexos.  Os valores iniciais dos
% pesos  são feitos exame da distribuição da entrada:  
  
  W = X(1:m).' ;  X = X((m+1):N) ; N = N-m ;

% como uma verificação, a contagem das vitórias para cada unidade
% de saída é  calculado na matriz " bate "  O valor previsto de
 % bate = N*#epochs/m.
    
  hits = zeros(m,1);
  
 % Um Mapa  Característico inicial :  O mapa é representado por um 
%medido grade pelos pesos das unidades de saída adjacentes traçadas
% no espaço da entrada   ou seja os pesos W(j) e W(k) são 
%conectados por um segmento da grade se V(j) e V(k) forem coordenadas
%dos  de unidades de saída adjacentes

% de valores iniciais do ganho do treinamento, do eta, e da propagação, sigma,
 % da função da vizinhança.

  epsilon = 0.5 ;   %  ganho do treinamento
  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(['epsilon = ', num2str(epsilon,2), ...
        '  sigma^2 = ', num2str(sg2/2,3)])
  hold off, drawnow
  
  %  fim de inicio
else  % continua
  epsilon = input('valor de entrada de epsilon [0.5]: ') ;
     if isempty(epsilon), epsilon = 0.5; end
  sg2 = input([' valor de entrada de  2sigma^2 [', ...
           num2str(sg2i), ']: ']) ;
     if isempty(sg2), sg2 = sg2i; end
end 
%  taxas de diminuir os valores do eta e do sigma durante o treinamento
% nós visamos valores finais:  eta = 0,1 e sigma = 1 (sg2=2)

reta = (0.2)^(2/N); rsigma = (1/sg2)^(2/N) ;

% **********************  malha principal ******************   

frm = 1;
for n = 1:N
      
% para cada teste padrão da entrada, X(n), e para cada unidade de saída
 % que armazenam o vetor W(v1 do peso, v2), a distância 
% entre X(n) e W é calcula
    
  WX = X(n) - W ; 
  
 %  coordenadas do neurônio ganhador, V(kn,:),  isto é, 
% o neurônio para que abs(WX) alcança o mínimo
  
[mnm kn] = min(abs(WX));  vkn = V(kn, :) ;
 
  hits(kn) = hits(kn)+1; % utilização de  neurônios

  % a função da vizinhança, NB, da forma do " sino ", 
% é centrado em torno da unidade V(kn ganhador,:)

  rho2 = sum(((vkn(ones(m, 1), :) - V).^2), 2) ;
  NB = exp(-rho2/sg2) ;

% finalmente, os pesos são atualizados de acordo com   
 % a lei de  aprendizagem de Kohonen:
  
  W = W + epsilon*NB.*WX ;

  % Valores de  "epsilon", e "sigma" são reduzidos
  % if (n==floor(N/2)), epsilon = 0.10; end
  if (n<N/2), %ordenamento e fase de convergencia
    sg2 = sg2*rsigma; 
  else 
    epsilon = epsilon*reta; 
  end
 
  % Cada  100  atualizações, o mapa caracteristico é ilustrado 
  if rem(n, 10) == 0
    plot([0 1],[0 1],'.'), grid, hold on, 
    if el == 1
      plot(W, 'b'), plot(W, '.r'),  
    else
      MF = full(sparse(V(:,1), V(:,2), W)) ;
      plot(MF, 'b'),  plot(MF.', 'r') ;
    end
    title(['epsilon = ',num2str(epsilon,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
    fprintf('-depsc2', '-f2', ['Jsom2Dt', num2str(frm)])
    frm = frm+1;
  end
 % MV(frm) = getframe; frm = frm+1 ;
end

% Apresentação Final  dos resultados
plot([0 1],[0 1],'.'), grid, hold on, 
if el == 1
  plot(W, 'b'), plot(W, '.r'),  
else
  MF = full(sparse(V(:,1), V(:,2), W)) ;
  plot(MF, 'b'),  plot(MF.', 'r') ;
end
title(' Mapa Caracteristico '), hold off , drawnow
 
% ***************  fim ************************

