function [ybest,xbest] = swarm_mod(fun,nd,valmin,valmax,n,m,Dt,K)
% SWARM_MOD acha o mínimo global de uma função utilizando o método do enxame 
%	de partículas modificado
%
% Chamada mínima:
% 	[ybest] = swarm(fun,nd,valmin,valmax)
%
% Chamada completa
% 	[ybest,xbest] = swarm(fun,nd,valmin,valmax,n,m,Dt,K)
%
% Entradas
%	fun - função a ser otimizada (string)
% 	nd - número de dimensões
%	valmin - valor mínimo para as variáveis de decisão 
%	valmax - valor máximo para as variáveis de decisão	
%	n - número de partículas (default = 10)
%	m - número de iterações ou instantesd e tempo (default = 20)
%	Dt - Intervalo de tempo (default = 1)
%	K - Constante de amortecimento (K >= 1) (default = 2)
%
% Saídas
%	ybest - valor mínimo da função encontrado
%	xbest - valor das variáveis de decisão que correspondem à ybest

% Marcelo Farenzena - 01/2005

% Teste anti - usuário
if nargin < 4
   error('Número de entradas insuficiente')
end

if nargin == 4
   n = 10;
   m = 20;
   K = 2;
   Dt = 1;
end

if nargin == 5
   m = 20;
   K = 2;
   Dt = 1;
end

if nargin == 6
   K = 2;
   Dt = 1;
end

if nargin == 7
   K = 2;
end

if isempty(fun)
	error('Insira corretamente o nome da função!')   
end

if isempty(nd)
	error('Insira corretamente o número de dimensões!')   
end

if isempty(valmin)
	error('Insira corretamente o vetor de valores mínimos!')   
end

if isempty(valmax)
	error('Insira corretamente o vetor de valores máximos!')   
end

if isempty(n)
   n = 10;
end

if isempty(m)	
   m = 20;
end
       
if isempty(Dt)
	Dt = 1;   
end

if isempty(K)
	K = 2;   
end

% etapa 0 - Inicializacao
for k = 1:n  
   lambda = rand(1,nd);   
   vin(k,:) = valmin + lambda.*(valmax-valmin);
end

% distribuido uniformemente
% nn = floor(sqrt(n));

% for k = 1:nn  
%    delta = (valmax(1) - valmin(1))/(nn-1);
%    vin((k-1)*nn+1:k*nn,1) = [valmin(1):delta:valmax(1)]';
%    vin((k-1)*nn+1:k*nn,2) = valmin(2) + (k-1)/(nn-1)*(valmax(2)-valmin(2));   
% end

% for k = nn^2:n  
%     lambda = rand(1,nd);   
%     vin(k,:) = valmin + lambda.*(valmax-valmin);
% end
% 


for i = 1:n
   fbest(i) = feval(fun,vin(i,:));
end
vbest = vin;

Vel = zeros(n,nd);

[temp,pxmin] = find(fbest == min(fbest));
fglobal = min(fbest);
vglobal = vin(pxmin(1),:);

% etapa 1 - iterando
for i=1:m
   w = sqrt(K - rand(n,nd).^2);	   
   psi = rand(n,nd);
   X = psi.*vbest + vecmmat(1-psi,vglobal);                 
   A = vin - X;
   B = (rand(n,nd).*A + Vel)./w;      
   f1 = exp(-rand(n,nd).*Dt).*cos(w*Dt);
   f2 = exp(-rand(n,nd).*Dt).*sin(w*Dt);
   vin = X + f1.*A + f2.*B;
   Vel = Vel.*f1 - (w.*A + rand(n,nd).*B).*f2;
   
   for k=1:nd
      [temp,vgt] = find(vin(:,k)>valmax(k));
      [temp,vlw] = find(vin(:,k)<valmin(k));
      vin(vgt,k) = valmax(k);
      vin(vlw,k) = valmin(k);      
      Vel(vgt,k) = 0;   
      Vel(vlw,k) = 0;
   end   
   
	for k =1:n
   	fatual(k) = feval(fun,vin(k,:));      
   	if fatual(k) < fbest(k)   
         fbest(k) = fatual(k);
         vbest(k,:) = vin(k,:);
      end      
   end
   
   [temp,pxmin] = find(fbest == min(fbest));
   if fglobal > min(fbest) 
      fglobal = min(fbest);
      vglobal = vin(pxmin(1),:);
   end           
end

ybest = fglobal;
xbest = vglobal;

function mmat = vecmmat(mat,vec)
% funcao que multiplica cada elemento do vetor vec por uma coluna
% da matriz mat

for i=1:length(vec)
	mmat(:,i) = mat(:,i)*vec(i);	      
end
