%% GWO |--//-- Grey Wolf Optimizer --//--| %%
function [Alpha_score,Alpha_pos,convergence,time_GWO] = GWO_orig_corr(npop,Maxiter,lb,ub,dim,fobj,Pos,initFit)
%Adapted by Pedro Bento November 2023
tic
% Initialization
wolves=Pos;
% fitness=zeros(1,npop);
iter=1;
convergence=zeros(1,Maxiter);
%--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
% Social hierarchy (ranking):
% The fittest solution as the alpha, the second and third best solutions are named beta and delta respectively.
% The rest of the candidate solutions are assumed to be omega.
% The hunting (optimization) is guided by alpha, beta, and delta. The omega wolves follow these three wolves.
%--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------%
% Determine the fitness of each wolve
% Fitness Evaluation: iter=0
% for i=1:npop
%     fitness(i)=fobj(wolves(i,:));
% end
fitness=initFit;
% Update Alpha, Beta, and Delta
[sortedfit,indx]=sort(fitness,'ascend'); %minimizar->ascend / minimizar->ascend
% Update alpha
Alpha_score=sortedfit(1);
Alpha_pos=wolves(indx(1),:);
convergence(iter)=Alpha_score;
% Update beta
Beta_score=sortedfit(2);
Beta_pos=wolves(indx(2),:);
% Update delta
Delta_score=sortedfit(3);
Delta_pos=wolves(indx(3),:);
% Main Loop
while (iter<Maxiter)
    % Wolves Position Update
    a=2-iter*((2)/Maxiter); % linearly decreasing  from 2 to 0
    for i=1:npop %varrer a população de grey wolves
        % for Alpha
        % random vectors in [0,1]
        r1=rand(1,dim);
        r2=rand(1,dim);
        % coefficient vectors for Alpha
        A1=-a+2.*a.*r1; % Equation (3.3)
        C1=2.*r2; % Equation (3.4)
        D_Alpha=abs(C1.*Alpha_pos-wolves(i,:)); % Equation (3.5)
        X1=Alpha_pos-A1.*D_Alpha; % Equation (3.6)
        % for Beta
        % random vectors in [0,1]
        r1=rand(1,dim);
        r2=rand(1,dim);
        % coefficient vectors for Alpha
        A2=-a+2.*a.*r1; % Equation (3.3)
        C2=2.*r2; % Equation (3.4)
        D_beta=abs(C2.*Beta_pos-wolves(i,:)); % Equation (3.5)
        X2=Beta_pos-A2.*D_beta; % Equation (3.6)
        % for Alpha
        % random vectors in [0,1]
        r1=rand(1,dim);
        r2=rand(1,dim);
        % coefficient vectors for Alpha
        A3=-a+2.*a.*r1; % Equation (3.3)
        C3=2.*r2; % Equation (3.4)
        D_Delta=abs(C3.*Delta_pos-wolves(i,:)); % Equation (3.5)
        X3=Delta_pos-A3.*D_Delta; % Equation (3.6)
        % Update i_th wolf position
        wolves(i,:)=(X1+X2+X3)/3; % Equation (3.7)
        % Return back the search agents that go beyond the boundaries of the search space
        Flag4ub=wolves(i,:)>ub;
        Flag4lb=wolves(i,:)<lb;
        wolves(i,:)=(wolves(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb;
        % Wolf fitness
        fitness(i)=fobj(wolves(i,:));
    end
    % Update Alpha, Beta, and Delta
    [sortedfit,indx]=sort(fitness,'ascend'); %minimizar->ascend / minimizar->ascend
    % Update alpha
    Alpha_score=sortedfit(1);
    Alpha_pos=wolves(indx(1),:);
    % Update beta
    Beta_score=sortedfit(2);
    Beta_pos=wolves(indx(2),:);
    % Update delta
    Delta_score=sortedfit(3);
    Delta_pos=wolves(indx(3),:);
    % store & update values
    iter=iter+1;
    convergence(iter)=Alpha_score;
end %iterations
time_GWO=toc;
end %function
