function [fmin,bsf_solution,convergence_curve,time_GSK]=GSK_orig_corr(pop_size,Max_iter,lb,ub,dim,fobj,Pos_old,InitFit)
%Adapted by Pedro Bento November 2023
tic
fitness=InitFit;
lu=[lb.*ones(1,dim); ub.*ones(1,dim)];
iter=1;
fmin=1e+300;
convergence_curve=zeros(1,Max_iter);
convergence_curve(iter)=min(fitness);
%%%%%%%%%%%%%%%%%%%%%%%% GSK: Parameter settings%%%%%%%%%%
G_Max=fix(Max_iter/pop_size);
KF=0.5;% Knowledge Factor
KR=0.9;%Knowledge Ratio
K=10*ones(pop_size,1);%Knowledge Rate
g=0;
children_fitness=zeros(pop_size,1);
%% main loop
while (iter<Max_iter)
    g=g+1;
    D_Gained_Shared_Junior=ceil((dim)*(1-g/G_Max).^K);
    D_Gained_Shared_Senior=dim-D_Gained_Shared_Junior;
    Pos = Pos_old; % the old population becomes the current population
    
    [valBest, indBest] = sort(fitness, 'ascend');
    [Rg1, Rg2, Rg3] = Gained_Shared_Junior_R1R2R3(indBest,pop_size);
    
    [R1, R2, R3] = Gained_Shared_Senior_R1R2R3(indBest,pop_size);
    R01=1:pop_size;
    Gained_Shared_Junior=zeros(pop_size, dim);
    ind1=fitness(R01)>fitness(Rg3);
    
    if(sum(ind1)>0)
        Gained_Shared_Junior (ind1,:)= Pos(ind1,:) + KF*ones(sum(ind1), dim) .* (Pos(Rg1(ind1),:) - Pos(Rg2(ind1),:)+Pos(Rg3(ind1), :)-Pos(ind1,:)) ;
    end
    ind1=~ind1;
    if(sum(ind1)>0)
        Gained_Shared_Junior(ind1,:) = Pos(ind1,:) + KF*ones(sum(ind1), dim) .* (Pos(Rg1(ind1),:) - Pos(Rg2(ind1),:)+Pos(ind1,:)-Pos(Rg3(ind1), :)) ;
    end
    R0=1:pop_size;
    Gained_Shared_Senior=zeros(pop_size, dim);
    ind=fitness(R0)>fitness(R2);
    if(sum(ind)>0)
        Gained_Shared_Senior(ind,:) = Pos(ind,:) + KF*ones(sum(ind), dim) .* (Pos(R1(ind),:) - Pos(ind,:) + Pos(R2(ind),:) - Pos(R3(ind), :)) ;
    end
    ind=~ind;
    if(sum(ind)>0)
        Gained_Shared_Senior(ind,:) = Pos(ind,:) + KF*ones(sum(ind), dim) .* (Pos(R1(ind),:) - Pos(R2(ind),:) + Pos(ind,:) - Pos(R3(ind), :)) ;
    end
    Gained_Shared_Junior = boundConstraint(Gained_Shared_Junior,lu,pop_size);
    Gained_Shared_Senior = boundConstraint(Gained_Shared_Senior,lu,pop_size);
    
    
    D_Gained_Shared_Junior_mask=rand(pop_size, dim)<=(D_Gained_Shared_Junior(:, ones(1, dim))./dim);
    D_Gained_Shared_Senior_mask=~D_Gained_Shared_Junior_mask;
    
    D_Gained_Shared_Junior_rand_mask=rand(pop_size, dim)<=KR*ones(pop_size, dim);
    D_Gained_Shared_Junior_mask=and(D_Gained_Shared_Junior_mask,D_Gained_Shared_Junior_rand_mask);
    
    D_Gained_Shared_Senior_rand_mask=rand(pop_size, dim)<=KR*ones(pop_size, dim);
    D_Gained_Shared_Senior_mask=and(D_Gained_Shared_Senior_mask,D_Gained_Shared_Senior_rand_mask);
    ui=Pos;
    
    ui(D_Gained_Shared_Junior_mask) = Gained_Shared_Junior(D_Gained_Shared_Junior_mask);
    ui(D_Gained_Shared_Senior_mask) = Gained_Shared_Senior(D_Gained_Shared_Senior_mask);
    
    for i = 1 : pop_size
        
        children_fitness(i) = fobj(ui(i,:));
        
        if children_fitness(i) < fmin
            fmin = children_fitness(i);
            bsf_solution = ui(i, :);
        end
    end
    
    [fitness, Child_is_better_index] = min([fitness, children_fitness], [], 2);
    
    Pos_old = Pos;
    Pos_old(Child_is_better_index == 2, :) = ui(Child_is_better_index == 2, :);
    
    iter=iter+1;
    convergence_curve(iter)=fmin;
end % end while loop
time_GSK=toc;
%________________________________________________________________%
%This function ensures that the bounds are respected
    function Pos=boundConstraint(Pos,lu,pop_size)
        % Check if solutions go outside the search space and bring them back
        for j=1:pop_size
            Flag4ub=Pos(j,:)>lu(2,:); %lu(2,:)==ub
            Flag4lb=Pos(j,:)<lu(1,:); %lu(1,:)==lb
            Pos(j,:)=(Pos(j,:).*(~(Flag4ub+Flag4lb)))+lu(2,:).*Flag4ub+lu(1,:).*Flag4lb;
        end
    end
%________________________________________________________________%
    function [R1, R2, R3] = Gained_Shared_Senior_R1R2R3(indBest,pop_size)
        R1=indBest(1:round(pop_size*0.1));
        R1rand = ceil(length(R1) * rand(pop_size, 1));
        R1 = R1(R1rand);
        R2=indBest(round(pop_size*0.1)+1:round(pop_size*0.9));
        R2rand = ceil(length(R2) * rand(pop_size, 1));
        R2 = R2(R2rand);
        R3=indBest(round(pop_size*0.9)+1:end);
        R3rand = ceil(length(R3) * rand(pop_size, 1));
        R3 = R3(R3rand);
    end
%________________________________________________________________%
    function [R1, R2, R3] = Gained_Shared_Junior_R1R2R3(indBest,pop_size)
        R0n=1:pop_size;
        R1=[];
        R2=[];
        R3=[];
        for hj=1:pop_size
            INDX=find(indBest==hj);
            if(INDX==1)% best
                R1(hj)=indBest(2);
                R2(hj)=indBest(3);
            elseif(INDX==pop_size)% worst
                R1(hj)=indBest(pop_size-2);
                R2(hj)=indBest(pop_size-1);
            else
                R1(hj)=indBest(INDX-1);
                R2(hj)=indBest(INDX+1);
            end
        end
        R3 = floor(rand(1, pop_size) * pop_size) + 1;
        for hj = 1 : 99999999
            pos = ((R3 == R2) | (R3 == R1) | (R3 == R0n));
            if sum(pos)==0
                break;
            else % regenerate r2 if it is equal to r0 or r1
                R3(pos) = floor(rand(1, sum(pos)) * pop_size) + 1;
            end
            if hj > 1000, % this has never happened so far
                error('Can not genrate R3 in 1000 iterations');
            end
        end
    end
%________________________________________________________________%

end