function [xmin,Emin,xhist,Ehist] = recozimento(func,x0,nemax,ntmax,nxmax,Tin,alfa,deltax,solviz)
% RECOZIMENTO - rotina que acha o mínimo de funções utilizando o método aleatório do
% 	             RECOZIMENTO SIMULADO
%
% [xmin,Emin,xhist,Ehist] = recozimento(func,x0,nemax,ntmax,nxmax,Tin,alfa,deltax,solviz)
%
% Parâmetros de entrada:
%		func   - string do nome da função objetivo a ser minimizada
%		x0     - chute inicial
%		nemax  - número máximo de execução do algoritmo (default = 5)
%		ntmax  - numero máximo de ciclos de resfriamento	(default = 20)
%		nxmax  - numero máximo de iteracoes em cada temperatura (default = 100)  								
%		Tin    - temperatura incial (grau de aleatoriedade, quanto maior, mais aleatório) (default = 100)
% 		alfa   - taxa de decréscimo da temperatura (default = 0.85)
% 		deltax - máxima variação da variável independente (default = 0.1)
% 		solviz - política de atualização da variável independente	
%			1 - xnovo = x + deltax * (rand(n,1)-0.5), onde n é o número de variáveis independentes
%			2 - xnovo = x * deltax * (rand(n,1)-0.5)
%
% Parâmetros de saida:
%		xmin  - mínimo da função	
% 		Emin  - Energia mínima da funcao (xmin aplicado à função objetivo) 		
%		xhist - Histórico das variáveis independentes utilizadas pela função (4 dimensões, xhist(ne,nt,nx,n))  
% 		Ehist - Energia da função para cada iteração (3 dimensões, Ehist(ne,nt,nx))

% Marcelo Farenzena
% 10 - 2004

% checagem inicial - o q o usuário fez de bobagens e o q omitiu
if nargin < 1 | isempty(func)
   error('A função precisa de uma função objetivo um chute inicial');
   break
end   
if nargin < 2 | isempty(x0)
   error('A função precisa de um chute inicial');
   break
end

if nargin == 2 
   nemax = 5;				% numero máximo de execuções do algoritmo
   ntmax = 20;				% numero máximo de ciclos de resfriamento
   nxmax = 100;			% numero máximo de iteracoes em cada temperatura
   Tin = 100;				% temperatura incial
   alfa = 0.85;			% alfa (T(k+1) = alfa*T(k))
   solviz = 1;				
   deltax = 0.1;			 
   
elseif nargin == 3 
   ntmax = 20;				% numero máximo de ciclos de resfriamento
   nxmax = 100;			% numero máximo de iteracoes em cada temperatura
   Tin = 100;				% temperatura incial
   alfa = 0.85;			% alfa (T(k+1) = alfa*T(k))
   solviz = 1;				
   deltax = 0.1;			 
   
elseif nargin == 4 
   nxmax = 100;			% numero máximo de iteracoes em cada temperatura
   Tin = 100;				% temperatura incial
   alfa = 0.85;			% alfa (T(k+1) = alfa*T(k))
   solviz = 1;				
   deltax = 0.1;			 
   
elseif nargin == 5
   Tin = 100;				% temperatura incial
   alfa = 0.85;			% alfa (T(k+1) = alfa*T(k))
   solviz = 1;				
   deltax = 0.1;			 
   
elseif nargin == 6
   alfa = 0.85;			% alfa (T(k+1) = alfa*T(k))
   solviz = 1;				
   deltax = 0.1;			 
   
elseif nargin == 7
   solviz = 1;				
   deltax = 0.1;			 
   
elseif nargin == 8
   deltax = 0.1;			 
   
end

if isempty(nemax)
   nemax = 5;			 	% numero máximo de execuções do algoritmo
end   
if isempty(ntmax)   
   ntmax = 20;				% numero máximo de ciclos de resfriamento
end   
if isempty(nxmax)   
   nxmax = 100;			% numero máximo de iteracoes em cada temperatura
end   
if isempty(Tin)   
   Tin = 100;				% temperatura incial
end   
if isempty(alfa)   
   alfa = 0.85;			% alfa (T(k+1) = alfa*T(k))
end   
if isempty(solviz)   
   solviz = 1;				
end   
if isempty(deltax)   
   deltax = 0.1;			 
end

% ==============================================================
% O algortimo começa aqui
% ==============================================================
ne = 0;						% numero de execucoes do algoritmo
n = length(x0);			% numero de variaveis independentes

E = feval(func,x0);	
Emin = E;
xmin = x0;

while ne<nemax
   ne = ne + 1;
   nt = 0;					% numero de ciclos de resfriamento inicial
   nx = 0;					% numero inicial de iteracoes em cada temperatura
   
   E = feval(func,x0);	
   x = x0;
   T = Tin;
   
   while nt<ntmax						% criterio de parada - numero de ciclos de resfriamento
      nt = nt + 1;
      nx = 0;
      while nx<nxmax					% numero de iteracoes em cada temperatura
         nx = nx + 1;   
         
         deltax = 2 * (rand(1,n) - 0.5);
         
         switch solviz				% metodo de calculo da solucao vizinha
         case 1						% soma um delta ao ponto atual em todas as direcoes
            xnovo = x + deltax;
         case 2         			% multiplica um delta ao ponto atual em todas as direcoes
            xnovo = (x + 100 * eps) * (1 + deltax);
         end
         
         Enovo = feval(func,xnovo); % calcula a nova energia
         
         delta = Enovo - E;			% calcula a variacao de energia
         
         if delta < 0     
            x = xnovo;
            E = Enovo;
            
            if E < Emin
               Emin = E;
               xmin = x;
            end
            
         else
            p = exp(-delta/T);
            if p > rand(1)
               x = xnovo;
               E = Enovo;         
            end
         end
         
        Ehist(ne,nx,nt) = Enovo;      
        xhist(ne,nx,nt,:) = xnovo;
%      	Extab = [Extab; Enovo xnovo];   
      end
      T = alfa*T;					% atualizacao da temperatura
   end   
end

% Emin2 = min(min(min(Ehist)));
% xmin



