function [xo,Ot,nS] = anneal(f,x0,mxit,ntmax,nxmax,Tin,alfa,deltax,solviz)
% SIMULATED ANNEALING - rotina que acha o mínimo de funções utilizando o método aleatório do
% 	             RECOZIMENTO SIMULADO
%
% [xo,Ot,nS] = anneal(f,x0,mxit,ntmax,nxmax,Tin,alfa,deltax,solviz)
%
% Parâmetros de entrada:
%		f   - string do nome da função objetivo a ser minimizada
%		x0     - chute inicial
%		mxit  - 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(nvi,1)-0.5), onde nvi é o número de variáveis independentes
%			2 - xnovo = x * deltax * (rand(nvi,1)-0.5)
%
%		Parâmetros de saida:
%		xo  - mínimo da função	
% 		Ot  - Energia mínima da funcao (xo aplicado à função objetivo) 		
%		xhist - Histórico das variáveis independentes utilizadas pela função (4 dimensões, xhist(ne,nt,nx,nvi))  
% 		Ehist - Energia da função para cada iteração (3 dimensões, Ehist(ne,nt,nx))

% Marcelo Farenzena
% 10 - 2004
%
% Modified by Giovani Tonel(giovani.tonel@ufrgs.br) on September 2006


ne = 0;								% numero de execucoes do algoritmo
nvi = length(x0);			% numero de variaveis independentes
nS=0;
E = feval(f,x0);	
Ot = E;
xo = x0;

while ne<mxit
   ne = ne + 1;
   nt = 0;					% numero de ciclos de resfriamento inicial
   nx = 0;					% numero inicial de iteracoes em cada temperatura
   
   E = feval(f,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,nvi) - 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(f,xnovo); % calcula a nova energia
         
         delta = Enovo - E;			% calcula a variacao de energia
         
         if delta < 0                %original: if delta < 0 
            x = xnovo;
            E = Enovo;
            
            if E < Ot
               Ot = E;
               xo = x;
            end
            
         else
            p = exp(-delta/T);
            if p > rand(1)
               x = xnovo;
               E = Enovo;         
            end
         end
       nS=nS+1; 
        
      end
      T = alfa*T;					% atualizacao da temperatura
   end   
end