function [xo,Ot,nS,bank]=swarm_stopcr(S,x0,ip,isb,Lb,Ub,R,nger,npas,problem,c1,c2,w0,stopcrit,ish,model_data)
%   Hard constrained global optimization using particle swarm optimization.
%
%   [xo,Ot,nS,bank]=swarm_confreg(S,x0,ip,isb,Lb,Ub,R,nger,npas,problem,c1,c2,w0,stopcrit,ish,ibank,model_data)
%
%   S: objective function
%   x0: inital point
%   ip: (0): no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip
%   isb: isb(i)=1 is i-th variable is bounded, otherwise isb(i)=0 (default = 0)
%   Lb, Ub: lower and upper bound vectors (default = x0*(1+/-2) for plot if isb(i)=0)
%   R: axis vector of the hyperellipse centered in x0 (default = max(0.1*abs(x0+~x0),1))
%   nger: number of generations or flights (default = 50)
%   npas: number of birds (default = 10)
%   problem: (-1): minimum (default), (1): maximum
%   c1: local search parameter (default = 1)
%   c2: global search parameter (default = 1)
%   w0: initial inertial factor (default = 0.9)
%   stopcrit: tolerance variation for average swarm position (default = 1e-6)
%   ish: (0) no text output (default), (1) show swarm status and new optima
%   model_data: any necessary information/constant for the function S
%
%   xo: optimal point
%   Ot: optimal values of S
%   nS: number of objective function evaluations
%   bank: history of all evaluated points (optional)

%   $Revision: 1.0 $  $Date: 2002/09/15 23:10:05 $
%   Author: Esly Ferreira da Costa Junior (versao FORTRAN, 10/2001)
%   MATLAB version: Argimiro R. Secchi (arge@enq.ufrgs.br)

%   $Revision: 2.0 $  $Date: 2017/07/21 14:09:16 $
%   Author: Lucas Ferreira Bernardino
%   - added support to confidence region analysis, text output showing
%   algorithm status, structure for passing extra arguments to the
%   objective function, new inertial factor profile, and stop criterion 
%   which measures the swarm's average position (Moraes et al., 2015)

% References:
%  - A.O.S. Moraes, J.F. Mitre, P.L.C. Lage, A.R. Secchi,
%   A robust parallel algorithm of the particle swarm optimization
%   method for large dimensional engineering problems, Applied Mathematical
%   Modelling, v. 39, p. 4223-4241, 2015.

if nargin < 2,
    error('swarm_stopcr requires two input arguments');
end

x0=x0(:);
n=size(x0,1);

if nargin < 3 || isempty(ip),
    ip=0;
end
if nargin < 4 || isempty(isb),
    isb=zeros(n,1);
end
if nargin < 5 || isempty(Lb),
    Lb=-x0-~x0;
end
if nargin < 6 || isempty(Ub),
    Ub=2*x0+~x0;
end
if nargin < 7 || isempty(R),
    R=max(0.1*abs(x0+~x0),1);
end
if nargin < 8 || isempty(nger),
    nger=50;
end
if nargin < 9 || isempty(npas),
    npas=10;
end
if nargin < 10 || isempty(problem),
    problem=-1;
end
if nargin < 11 || isempty(c1),
    c1=1;
end
if nargin < 12 || isempty(c2),
    c2=1;
end
if nargin < 13 || isempty(w0),
    w0=0.9;
end
if nargin < 14 || isempty(stopcrit),
    stopcrit=1e-6;
end
if nargin < 15 || isempty(ish),
    ish=0;
end
if nargin < 16 || isempty(model_data),
    Fobj=@(x_)(feval(S,x_)*problem);
else
    Fobj=@(x_)(feval(S,x_,model_data)*problem);
end
ibank=(nargout > 3);

% Local variables:
%   ipg  = index of best position
%   tw   = inertial factor reduction rate
%   p    = matrix of best position of each bird
%   y    = vector of best objective function (y=Fobj(p))
%   v    = bird velocity matrix
%   x    = bird location matrix
%   w    = inertial factor
%   bank = memory of every evaluation

% initialization
R=abs(R(:));
Lb=Lb(:);
Ub=Ub(:);
Ot=feval(Fobj,x0);
MaxF=Ot;
y=zeros(1,npas);
y(1)=Ot;
p=zeros(n,npas);
v=zeros(n,npas);
w=w0;
x=zeros(n,npas);
nS=1;
if ibank
    bank=zeros((nger+1)*npas,length(x0)+1);
    bank(nS,:)=[x0' Ot*problem];
else
    bank=[];
end
acompanha=[];

if ip && n == 2, %gero curvas de nível
    figure(abs(ip));
    [X1,X2]=meshgrid(Lb(1):(Ub(1)-Lb(1))/20:Ub(1),Lb(2):(Ub(2)-Lb(2))/20:Ub(2));
    [n1,n2]=size(X1);
    f=zeros(n1,n2);
    for i=1:n1,
        for j=1:n2,
            f(i,j)=feval(Fobj,[X1(i,j);X2(i,j)])*problem;
        end
    end
    mxf=max(max(f));
    mnf=min(min(f));
    df=mnf+(mxf-mnf)*(2.^(((0:10)/10).^2)-1);
    [vl,h]=contour(X1,X2,f,df); hold on;
    h1=plot(x0(1),x0(2),'ro');
    legend(h1,'start point');
    
    if ip > 0,
        disp('Pause: hit any key to continue...'); pause;
    end
end

del=Ub-Lb;

if ish, fprintf('Starting particles with random velocities...\n'), end
for j=1:npas, %gero posições iniciais para as partículas
    x(:,j)=Lb+rand(n,1).*del;
    v(:,j)=0.1*rand(n,1).*del; %0-10% of domain size
    p(:,j)=x(:,j);
end

x(:,1)=x0;
p(:,1)=x0;
ipg=1;

for j=2:npas, %calculo a função objetivo inicial para as partículas (conto tempo)
    y(j)=feval(Fobj,x(:,j));
    nS=nS+1;
    if ibank, bank(nS,:)=[x(:,j)' y(j)*problem]; end
    
    if y(j) > Ot,
        Ot=y(j);
        ipg=j;
        if ish, fprintf('**** New optimum found by particle %.0f: %.3e\n',j,Ot*problem), end
    elseif y(j) < MaxF,
        MaxF=y(j);
    end
end
val=y;

if ip && n == 2, %ploto o primeiro movimento do ótimo
    plotim=animatedline([x(1,ipg) x0(1)],[x(2,ipg) x0(2)],'Color','r');
    for iplot=1:numel(acompanha)
        plpath{iplot}=animatedline(x(1,acompanha(iplot)),x(2,acompanha(iplot)),'Color','k');
    end
    drawnow
    if ip > 0,
        disp('Pause: hit any key to continue...'); pause;
    end
end

xo=p(:,ipg);
xcomb=[(x-Lb*ones(1,npas))./(del*ones(1,npas)+eps); (val-Ot)/(MaxF-Ot)];
xcombg=[((xo-Lb)./(del+eps)); 0];
inorm=1./sqrt(sum((xcomb-xcombg*ones(1,npas)).^2,1));
inorm(isinf(inorm))=0;
wp=inorm/sum(inorm);
xcombmed_old=sum((ones(n+1,1)*wp).*xcomb,2);

for ig=1:nger,
    if ish, fprintf('-------------- Generation %.0f --------------\n',ig), end
    % gero velocidades e movo partículas
    rnd1=ones(n,1)*rand(1,npas);
    rnd2=ones(n,1)*rand(1,npas);
    v=w*v+c1*rnd1.*(p-x)+c2*rnd2.*(xo*ones(1,npas)-x);
    x=x+v;
    
    for i=1:n, %checo se as partículas bateram na parede i e reseto a velocidade daquela direção
        if isb(i),
            j=find(x(i,:) > Ub(i));
            if ~isempty(j),
                x(i,j)=Ub(i)-0.05*stopcrit*del(i);
                v(i,j)=0;
            end
            
            j=find(x(i,:) < Lb(i));
            if ~isempty(j),
                x(i,j)=Lb(i)+0.05*stopcrit*del(i);
                v(i,j)=0;
            end
        end
    end
    
    for j=1:npas, %calculo e guardo os valores da função objetivo (e conto tempo)
        val(j)=feval(Fobj,x(:,j));
        nS=nS+1;
        if ibank, bank(nS,:)=[x(:,j)' val(j)*problem]; end
        
        if val(j) > y(j),
            y(j)=val(j);
            p(:,j)=x(:,j);
            
            if val(j) > Ot,
                Ot=val(j);
                ipg=j;
                xo=p(:,ipg);
                if ish, fprintf('**** New optimum found by particle %.0f: %.3e\n',j,Ot*problem), end
            end
        end
    end
    
    if ip && n == 2, %adiciono pontos ao plot
        addpoints(plotim,p(1,ipg),p(2,ipg))
        for iplot=1:numel(acompanha)
            addpoints(plpath{iplot},x(1,acompanha(iplot)),x(2,acompanha(iplot)));
        end
        drawnow
        if ip > 0,
            disp('Pause: hit any key to continue...'); pause;
        end
    end
    
    xcomb=[(x-Lb*ones(1,npas))./(del*ones(1,npas)+eps); (val-Ot)/(MaxF-Ot)];
    xcombg=[((xo-Lb)./(del+eps)); 0];
    inorm=1./sqrt(sum((xcomb-xcombg*ones(1,npas)).^2,1));
    inorm(isinf(inorm))=0;
    wp=inorm/sum(inorm);
    xcombmed=sum((ones(n+1,1)*wp).*xcomb,2);
    if ish
        fprintf('delta_y_av: %.3e\n',norm(xcombmed-xcombmed_old))
        fprintf('delta_F_av: %.3e\n',abs(sum(Ot-val)/npas))
    end
    if (norm(xcombmed-xcombmed_old)<stopcrit), break, end
    xcombmed_old=xcombmed;
    w=(w0/0.9841)*(1/2+atan((-ig/nger+2/3)*30)/pi);
end

if ip && n == 2, %ploto ponto ótimo
    h2=plot(xo(1),xo(2),'r*');
    legend([h1,h2],'start point','optimum');
end

Ot=Ot*problem;
if ibank, bank=bank(1:nS,:); end

end
