function [yout,xout,it]=becam(X,Y,xin,it_max,a)

%
% Memória Associativa Bidirecional com Capacidade Exponencial (BECAM)
%
% Comando: [xout,it]=becam(X,Y,xin,it_max,a)
% 
% Entradas:
%      X,Y  - Matrizes contendo as memórias fundamentais como coluna,
%       xin - Padrão de entrada,
%    it_max - Número máximo de iteração (padrão: it_max=N, núm. de neurônios),
%        a  - Coeficiente da exponencial (padrão: a=exp(1/N))
%   
% Saídas:
% yout,xout - Padrão recordado pela memória associativa,
%        it - Número de iterações efetuadas,
%

% Número de neurônios da primeira camada,
M=size(Y,1);
% Número de neurônios da segunda camada,
N=size(X,1);

% Define o número máximo de iterações,
if nargin<4
    it_max=max(N,M);
end

if nargin<5
   a=1/N;
else
   a=log(a);
end

% Entrada yin,
yin=sign(Y*exp(a*(X'*xin)));

% it, conta o numero de iterações,
it=0;
% Usado como critério de parada:
Erro=1;
% Calcular até convergir ou até o número máximo de iterações,
while (Erro>1.e-4)&(it<it_max)
    % Conta o número da iterações,
    it=it+1;
    % xout,yout, serão o resultado no final da iteração,     
    xout=xin; yout=yin;
    % Escolhe uma seqüência de índices aleatoriamente,
    ind=randperm(N+M);
    % Para cada índice da seqüência gerada aleatoriamente, faça:
    for i=1:N+M
       % Calcula a ativação do neurônio com certo índice,
       if ind(i)<=M
          u=Y(ind(i),:)*exp(a*(X'*xin));
          yout(ind(i))=sign(u);
          % Calcula a saída do neurônio e atualiza xout.
          % Se u=0, então y(t+1)_i = y(t)_i,
          if abs(u)<1.e-9 % Evitar erros de arredondamento! 
             yout(ind(i))=yin(ind(i));
          end
       else
          ind(i)=ind(i)-M;
          u=X(ind(i),:)*exp(a*(Y'*yin));
          xout(ind(i))=sign(u);
          % Calcula a saída do neurônio e atualiza xout.
          % Se u=0, então x(t+1)_i = x(t)_i,
          if abs(u)<1.e-9 % Evitar erros de arredondamento! 
             xout(ind(i))=xin(ind(i));
          end
       end
    end
    % Calcula o erro,    
    Erro=max(norm(xout-xin),norm(yout-yin));
    % Atualiza x(t+1) e y(t+1),
    xin=xout; yin=yout;
 end