function [ybest,xbest] = swarm(fun,nd,valmin,valmax,n,m,wfinal,winicial,cs1,cs2)
% SWARM acha o mínimo global de uma função utilizando o método do enxame 
%	de partículas
%
% Chamada mínima:
% 	[ybest] = swarm(fun,nd,valmin,valmax)
%
% Chamada completa
% 	[ybest,xbest] = swarm(fun,nd,valmin,valmax,n,m,wfinal,winicial,cs1,cs2)
%
% 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)
%	wfinal - ponderação da velocidade anterior da partícula na iteração 
%		final (default = 0)
%	winicial - ponderação da velocidade anterior da partícula na iteração 
%		inicial (default = 0.9)
%	cs1 - ponderação da crença da partícula em sua melhor posição
%		(default = .5)
%	cs2 - ponderação da crença da partícula na melhor posição das outras
%		(default = .5)
%
% 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;
   wfinal = 0;
   winicial = 0.9;
   cs1 = 0.5;
   cs2 = 0.5;  
end

if nargin == 5
   m = 20;
   wfinal = 0;
   winicial = 0.9;
   cs1 = 0.5;
   cs2 = 0.5;  
end

if nargin == 6
   wfinal = 0;
   winicial = 0.9;
   cs1 = 0.5;
   cs2 = 0.5;  
end

if nargin == 7
   winicial = 0.9;
   cs1 = 0.5;
   cs2 = 0.5;  
end

if nargin == 8
   cs1 = 0.5;
   cs2 = 0.5;  
end

if nargin == 9
   cs2 = 0.5;  
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(wfinal)
	wfinal = 0;   
end

if isempty(winicial)
	winicial = 0.9;   
end

if isempty(cs1)
	cs1 = 0.5;   
end

if isempty(cs2)
	cs2 = 0.5;   
end  
   
% etapa 0 - Inicializacao
% distribuido randomicamente
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 = winicial + (wfinal - winicial)*((i-1)/(m-1));	   
	for k =1:n
      Vel(k,:) = w*Vel(k,:) + cs1*rand(1,nd).*(vbest(k,:) - vin(k,:)) + ...
         cs2*rand(1,nd).*(vglobal - vin(k,:));           
   end
   vin = vin + Vel;   
   
   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      
   vglobal = vglobal(1,:);
end

ybest = fglobal;
xbest = vglobal;

