function rho(caso,dim,vetk,simulacoes,salvar,op_morf)

%
% Raio da Bacia de Atração
%
% Esta rotina estima o raio da bacia de atração 
% dos modelos de memória associativa auto e hetero-associativas.
%
% Comando: rho(caso,dim,vetk,simulacoes,salvar,op_morf);
%
% Entradas: 
%       Caso - Caso auto-associativo (caso=1) ou caso heteroassociativo (caso=2);
%              Padrão: Caso = 1;
%        dim - Dimensão dos padrões de entrada e saída ;
%              Padrão: dim = [12,10];
%       vetk - Vetor contendo o número de pares de padrões que serão 
%                  armazenados nas memórias associativas (Padrão: vetk = [1:n]);
% simulacoes - Número de simulações (Padrão: simulacoes = 100)
%     salvar - Salva os resultados no arquivo fdi_auto.mat 
%                ou fdi_hetero.mat quando salvar=1 (padrão: salvar=1);
%    op_morf - Impõe X(:,1)=zeros(n,1) ou X(:,1)=ones(n,1) nas memórias 
%                associativas morfológicas quando op_morf=1 (Padrão: op_morf=1)
%

if nargin<6
   op_morf=1;
end

if nargin<5
   salvar=1;
end

if nargin<4
   simulacoes=100;
end

if nargin<2
   dim = [12,10];
end

if nargin<3
   vetk=[1:dim(1)];
end

if nargin<1
   caso=1;
end

n=dim(1);
if length(dim)<2
   m=dim(1);
else
   m=dim(2);
end

if length(vetk)==1
   vetk(2)=vetk(1);
end

for s=1:simulacoes
   tic
   for indice=1:length(vetk);
      k=vetk(indice);
      % Gerando padroes bipolares;
      clear X;
      X=sign(2*rand(n,k)-1);
            
      if caso==1
         % Caso Auto-Associativo
         
         %
         % Memória Associativa de Hopfield;
         %
      	        
      	% Fase de Armazenamento;
      	clear W;
         W=corr0(X);
         
      	clear ind;
      	ind=randperm(n);
      	R_hopfield(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear xr2; xr2=mad(W,x);
      	% Erro da Recordação:
      	erro=norm(X(:,1)-xr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_hopfield(s,indice)<n)
            R_hopfield(s,indice)=R_hopfield(s,indice)+1;
            % R_hopfield representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_hopfield(s,indice)))=(-1)*x(ind(R_hopfield(s,indice)));
         
            % Fase de Recordação:
            clear xr2; xr2=mad(W,x);
            % Erro da recordação:
            erro=norm(X(:,1)-xr2);
         end
         
         %
         % Memória Associativa de Personnaz;
         %
      	        
      	% Fase de Armazenamento;
      	clear W;
         W=proj(X,X);
         
      	clear ind;
      	ind=randperm(n);
      	R_personnaz(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear xr2; xr2=mad(W,x);
      	% Erro da Recordação:
      	erro=norm(X(:,1)-xr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_personnaz(s,indice)<n)
            R_personnaz(s,indice)=R_personnaz(s,indice)+1;
            % R_personnaz representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_personnaz(s,indice)))=(-1)*x(ind(R_personnaz(s,indice)));
         
            % Fase de Recordação:
            clear xr2; xr2=mad(W,x);
            % Erro da recordação:
            erro=norm(X(:,1)-xr2);
         end
         
         %
         % Memória Associativa de Kanter-Sompolinsky;
         %
      	        
      	% Fase de Armazenamento;
      	clear W;
         W=proj0(X);
         
      	clear ind;
      	ind=randperm(n);
      	R_kanter(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear xr2; xr2=mad(W,x);
      	% Erro da Recordação:
      	erro=norm(X(:,1)-xr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_kanter(s,indice)<n)
            R_kanter(s,indice)=R_kanter(s,indice)+1;
            % R_personnaz representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_kanter(s,indice)))=(-1)*x(ind(R_kanter(s,indice)));
         
            % Fase de Recordação:
            clear xr2; xr2=mad(W,x);
            % Erro da recordação:
            erro=norm(X(:,1)-xr2);
         end
         
         %
         % Memória Associativa com capacidade Exponencial;
         %
         
      	clear ind;
      	ind=randperm(n);
      	R_ecam(s,indice)=0;
      	clear x; x=X(:,1);
        
         % Fase de Recordação; 
         %  Cuidado com os parametros da ECAM!
      	clear xr2; xr2=ecam(X,x);
      	% Erro da Recordação:
      	erro=norm(X(:,1)-xr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_ecam(s,indice)<n)
            R_ecam(s,indice)=R_ecam(s,indice)+1;
            % R_ecam representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_ecam(s,indice)))=(-1)*x(ind(R_ecam(s,indice)));
         
            % Fase de Recordação:
            clear xr2; xr2=ecam(X,x);
            % Erro da recordação:
            erro=norm(X(:,1)-xr2);
         end
         
         %
         % Modelos Morfológicos;
         %
         % Convertendo os padrões bipolares em padrões binários:
         X=(X+1)/2;
         
         %
         % Memória Auto-associativa Morfológica Wxx;
         %
         
         if op_morf==1
            X(:,1)=ones(n,1);
         end
                  
      	% Fase de Armazenamento;
      	clear W; 
      	W=boxmin(X,-X');
      	
      	clear ind;
      	ind=randperm(n);
      	R_Wxx(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear xr2;
      	xr2=boxmax(W,x);
      	% Erro da Recordação:
      	erro=norm(X(:,1)-xr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_Wxx(s,indice)<n)
            R_Wxx(s,indice)=R_Wxx(s,indice)+1;
            % R_Wxx representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_Wxx(s,indice)))=1-x(ind(R_Wxx(s,indice)));
            
            % Fase de Recordação:
            xr2=boxmax(W,x);
            % Erro da recordação:
            erro=norm(X(:,1)-xr2);
         end
         
         %
         % Memória Associativa Morfológica de Duas Camadas;
         %
         
         if op_morf==1
            X(:,1)=zeros(n,1);
         end
         
         clear Z;
      	Z=eye(k);
        
      	% Fase de Armazenamento;
      	clear Wzy; clear Mxxz;
      	Wzy=boxmin(X,-Z');
      	Mxxz=boxmin(boxmax(Z,-X'),boxmax(X,-X'));
      
      	clear ind;
      	ind=randperm(n);
      	R_mam2c1(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear xr2;
      	xr2=mam2c(Wzy,Mxxz,x);
      	% Erro da Recordação:
      	erro=norm(X(:,1)-xr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_mam2c1(s,indice)<n)
            R_mam2c1(s,indice)=R_mam2c1(s,indice)+1;
            % R_mam2c1 representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_mam2c1(s,indice)))=1-x(ind(R_mam2c1(s,indice)));
         
            % Fase de Recordação:
            xr2=mam2c(Wzy,Mxxz,x);
            % Erro da recordação:
            erro=norm(X(:,1)-xr2);
         end
      end
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      if caso==2
         % Caso Hetero-Associativo
         
         % Gerando padroes bipolares;
         clear Y; Y=sign(2*rand(m,k)-1);
         
         %
         % Memória Associativa Bidirecional;
         %
      	        
      	% Fase de Armazenamento;
      	clear W;
         W=corr(X,Y);
         
      	clear ind;
      	ind=randperm(n);
      	R_bam(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear yr2; yr2=madb(W,W',x);
      	% Erro da Recordação:
      	erro=norm(Y(:,1)-yr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_bam(s,indice)<n)
            R_bam(s,indice)=R_bam(s,indice)+1;
            % R_bam representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_bam(s,indice)))=(-1)*x(ind(R_bam(s,indice)));
         
            % Fase de Recordação:
            clear yr2; yr2=madb(W,W',x);
            % Erro da recordação:
            erro=norm(Y(:,1)-yr2);
         end
         
         %
         % Memória Associativa Bidirecional Assimétrica;
         %
      	        
      	% Fase de Armazenamento;
         clear W1; W1=proj(X,Y);
         clear W2; W2=proj(Y,X);
         
      	clear ind;
      	ind=randperm(n);
      	R_abam(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear yr2; yr2=madb(W1,W2,x);
      	% Erro da Recordação:
      	erro=norm(Y(:,1)-yr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_abam(s,indice)<n)
            R_abam(s,indice)=R_abam(s,indice)+1;
            % R_abam representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_abam(s,indice)))=(-1)*x(ind(R_abam(s,indice)));
         
            % Fase de Recordação:
            clear yr2; yr2=madb(W1,W2,x);
            % Erro da recordação:
            erro=norm(Y(:,1)-yr2);
         end
                
         %
         % Memória Associativa Bidirecional com capacidade Exponencial;
         %
         
      	clear ind;
      	ind=randperm(n);
      	R_becam(s,indice)=0;
      	clear x; x=X(:,1);
        
         % Fase de Recordação; 
         %  Cuidado com os parametros da BECAM!
      	clear yr2; yr2=becam(X,Y,x);
      	% Erro da Recordação:
      	erro=norm(Y(:,1)-yr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_becam(s,indice)<n)
            R_becam(s,indice)=R_becam(s,indice)+1;
            % R_becam representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_becam(s,indice)))=(-1)*x(ind(R_becam(s,indice)));
         
            % Fase de Recordação:
            clear yr2; yr2=becam(X,Y,x);
            % Erro da recordação:
            erro=norm(Y(:,1)-yr2);
         end
         
         %
         % Modelos Morfológicos;
         %
         % Convertendo os padrões bipolares em padrões binários:
         X=(X+1)/2;
         Y=(Y+1)/2;
         
         %
         % Memória Auto-associativa Morfológica Wxy;
         %
         
         if op_morf==1
            X(:,1)=ones(n,1);
         end
                  
      	% Fase de Armazenamento;
      	clear W; 
      	W=boxmin(Y,-X');
      	
      	clear ind;
      	ind=randperm(n);
      	R_Wxy(s,indice)=0;
      	clear x; x=X(:,1);
        
      	% Fase de Recordação;      
      	clear yr2;
      	yr2=boxmax(W,x);
      	% Erro da Recordação:
      	erro=norm(Y(:,1)-yr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_Wxy(s,indice)<n)
            R_Wxy(s,indice)=R_Wxy(s,indice)+1;
            % R_Wxy representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_Wxy(s,indice)))=1-x(ind(R_Wxy(s,indice)));
            
            % Fase de Recordação:
            yr2=boxmax(W,x);
            % Erro da recordação:
            erro=norm(Y(:,1)-yr2);
         end
         
         %
         % Memória Associativa Morfológica de Duas Camadas;
         %
         
         if op_morf==1
            X(:,1)=zeros(n,1);
         end
         
         clear Z;
      	Z=eye(k);
        
      	% Fase de Armazenamento;
      	clear Wzy; clear Mxxz;
      	Wzy=boxmin(Y,-Z');
      	Mxxz=boxmin(boxmax(Z,-X'),boxmax(X,-X'));
      
      	clear ind;
      	ind=randperm(n);
      	R_mam2c2(s,indice)=0;
      	clear x; x=X(:,1);
        
	      % Fase de Recordação;      
      	clear yr2;
      	yr2=mam2c(Wzy,Mxxz,x);
      	% Erro da Recordação:
      	erro=norm(Y(:,1)-yr2);
      
      	% Calculando o max (ou sup) da equação 7.29;
      	while (erro<1.e-4)&(R_mam2c2(s,indice)<n)
            R_mam2c2(s,indice)=R_mam2c2(s,indice)+1;
            % R_mam2c2 representa o número de componentes alteradas do padrão de entrada,
            % Alterando as compentes do vetor de entrada!
            x(ind(R_mam2c2(s,indice)))=1-x(ind(R_mam2c2(s,indice)));
         
            % Fase de Recordação:
            yr2=mam2c(Wzy,Mxxz,x);
            % Erro da recordação:
            erro=norm(Y(:,1)-yr2);
         end
      end
   end   
   disp(sprintf('Faltam %d simulações. O tempo estimado para completar a execução é %d min.',simulacoes-s,ceil((simulacoes-s)*toc/60)));
end

figure(1);
clf;
if caso==1
   % O resultado final é a média R^k (Eq. 7.30);
   rho_hopfield=max(0,mean(R_hopfield)-1);
   rho_personnaz=max(0,mean(R_personnaz)-1);
   rho_kanter=max(0,mean(R_kanter)-1);
   rho_ecam=max(0,mean(R_ecam)-1);
   rho_Wxx=max(0,mean(R_Wxx)-1);
   rho_mam2c1=max(0,mean(R_mam2c1)-1);
   
   figure(1);
   plot(vetk,rho_hopfield,vetk,rho_personnaz,vetk,rho_kanter,vetk,rho_ecam,vetk,rho_Wxx,vetk,rho_mam2c1);
   legend('MA Hopfield','MA Personnaz','MA Kanter-Somp.','ECAM','MAM Wxx','MAM 2 Camadas');
   xlabel('Número de Memórias Fundamentais Armazenadas','FontSize',12);
   ylabel('Raio de Atração','FontSize',12);
   
   if salvar==1
      save rho_auto rho_hopfield rho_personnaz rho_kanter rho_ecam rho_Wxx rho_mam2c1 vetk dim simulacoes;
      disp('Confira os resultados no arquivo rho_auto.mat')
   end   
end

if caso==2
   % O resultado final é a média R^k (Eq. 7.30);
   rho_bam=max(0,mean(R_bam)-1);
   rho_abam=max(0,mean(R_abam)-1);
   rho_becam=max(0,mean(R_becam)-1);
   rho_Wxy=max(0,mean(R_Wxy)-1);
   rho_mam2c2=max(0,mean(R_mam2c2)-1);
   
   figure(1);
   plot(vetk,rho_bam,vetk,rho_abam,vetk,rho_becam,vetk,rho_Wxy,vetk,rho_mam2c2);
   legend('BAM','ABAM','BECAM','MAM Wxy','MAM 2 Camadas');
   xlabel('Número de Memórias Fundamentais Armazenadas','FontSize',12);
   ylabel('Raio de Atração','FontSize',12);
   
   if salvar==1
      save rho_hetero rho_bam rho_abam rho_becam rho_Wxy rho_mam2c2 vetk dim simulacoes;
      disp('Confira os resultados no arquivo rho_hetero.mat')
   end   
end