function [DF,Xb,Convergence_curve,elapsed_t,globDB]=SurrAssistedSMAv3_f2_4b_corr7(Max_FEs,pop_size,lb,ub,dim,fobj,Pos,initFitn,DB_maxsize)
%% Slime Mould Algorithm: Stochastic Optimization SurrAssistedSMAv3_f2_4b_corr7(Implementation by Pedro Bento)
% v3: Modified implementation of the SMA v1, with corrected new position generation
% f2: Implies that only positions evaluated with the real fobj are updated after the iteration, i.e., not use the approx fobj as a S(i) indication
% 4: DBmanagement based on a given set of rules (dynDB_managem_v4_A)
% corr: returns the convergence curve corrected and adjusted to the FEs counter during each iteration and not based on the iterations number
%       / also the function recieves the DB_maxsize parameter as an input from the main
% 7: globsurrogate_oper_vf7 (without newrb_mod, with
time_start=tic;
%_______________________________________________________________________________________________________________________%
% initialize variables working variables
[S,idx]=sort(initFitn,'ascend'); %sorted fitness to rank the best slime mould position, S(i)-> fitness
X=Pos(idx,:); % X(i)-> positions
Xtemp=zeros(pop_size,dim); %variable created to temporarily store the recently updated SMA positions
DF=S(1); % DF-> best fitness/destination fitness (in all iterations)
bF=S(1); % bF-> best fitness in the current iteration/swarm positioning
wF=S(end); % wF-> worst fitness in the current iteration/swarm positioning
Xb=X(1,:);  % Xb-> best position/search agent so far
W=zeros(pop_size,dim);
z=0.03; % parameter 3% of the times a new particle is generated
%% Other results related variables % % % % % % % % % % % % % %
% Convergence_curve=zeros(1,max_iters); %old
% Convergence_curve(1)=DF;
Convergence_curve=inf(1,Max_FEs); %correction
Convergence_curve(1:pop_size)=DF;
iter=1; % Iteration/Loop counter
FEs_counter=pop_size; % counter of the amount of times that the real fobjs was evaluated FEs
appr_fit=zeros(pop_size,1); %fitness de cada posicao do current swarm avaliado com a fobj aproximada (melhor delas)
%% -Init Global Database and Fmerit related variables
%-> Global database
% DB_maxsize=500;
[globDB,~]=dynDB_managem_v4_A('init',[],[],pop_size,[],X,S,[],[],iter); % Init DB with the sorted initial positioning from the SMA
%auxiliary variable for the weight given to the merit
wmer_min=0.35; % min weight factor
wmer_max=0.95; % max weight factor
%% - Estimate the initial global surrogate with a RBFNN
[globSurr,~,surr_trainset_size]=globsurrogate_oper_vf7(0,iter,FEs_counter,Max_FEs,[],globDB.Pos,globDB.Fit,[]);
%________________________________________________________%
%% Main metaheuristic loop
while (FEs_counter<Max_FEs) % stop crit
    %-----------------------------------------------------------------%
    if (iter>1) %for iter==1 it is unnecessary (already done)
        [S,idx]=sort(S,'ascend'); % (equation 2.6)
        X=X(idx,:); % X(i)-> positions ranked by fitness (best to worst)
        bF=S(1);
        wF=S(end);
    end
    if (S(1)<DF)  %Verify if the fitness of the current best position X(1) is better than the overall best fitness reached (seen) so far
        DF=S(1);
        Xb=X(1,:);
    end
    %-----------------------------------------------------------------%
    %- Weight of slime mould calculation (equation 2.6)
    for i=1:pop_size
        if (i<=round(pop_size/2)) % First half is the better ranked
            W(i,:)=1+rand*log10((bF-S(i))/(bF-wF+eps)+1);
        else %(i>round(pop_size/2)) %Lowest half
            W(i,:)=1-rand*log10((bF-S(i))/(bF-wF+eps)+1);
        end
    end
    %---------------------------------------------------------------------%
    %- Auxiliary variables
    a=atanh((-FEs_counter/Max_FEs)+1);  % (equation 2.4)
    b=(Max_FEs-FEs_counter)/Max_FEs;
    %---------------------------------------------------------------------%
    %- Search Agents Position Update
    for i=1:pop_size
        %** Approach food position update mechanism
        p=tanh(abs(S(i)-DF)); %treshold to select either position update formula (equation 2.2)
        %- Position Update (equations 2.7/2.1)
        if (rand<z) %generates a new slime mould position
            Xtemp(i,:)=lb+(ub-lb)*rand.*ones(1,dim); % biased SMA paper interpretation (all dims assume the same new random value)
            %             Xtemp(i,:)=lb+(ub-lb).*rand(1,dim); % my interpretation (worst results)
        else
            r=rand;
            if (r<p)
                vb=-a+2*a.*rand(1,dim); % vb is a vector where its a different parameter per dimension within the range [-a,a]
                %                 aux0=randi(pop_size,1,2); % SMA paper interpretation, i.e., two randomly selected positions, i.e., (by chance) can be the same ones!
                aux0=randperm(pop_size,2); % my interpretation, i.e., two different randomly selected positions
                XA=X(aux0(1),:);
                XB=X(aux0(2),:);
                Xtemp(i,:)=Xb+vb.*(W(i,:).*XA-XB);
            else %r>=p
                vc=-b+2*b.*rand(1,dim); %vc is a vector with random numbers (from continous uniform distribution) with it«s range linearly decreasing from 1 to 0
                Xtemp(i,:)=vc.*X(i,:);
            end
            if (~isreal(Xtemp(i,:)))
                disp('Some kind of error')
            end
        end
    end
    %---------------------------------------------------------------------%
    %- Verify if the new search agents positions are within the search space boundaries
    %- & Determine the search agents fitness and evaluate
    for i=1:pop_size
        Flag4ub=Xtemp(i,:)>ub;
        Flag4lb=Xtemp(i,:)<lb;
        Xtemp(i,:)=(Xtemp(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb; %corrected if needed
    end
    %-------------------------------------------------------------------------%
    count_realfit=zeros(pop_size,1); %control variable that checks if during a given (each) iteration each agent was updated or not using the real fobj
    %-------------------------------------------------------------------------%
    for i=1:pop_size % for all/each of the search agents
        %---------------------------------------------------------------------%
        %- Calculate the new objective function value for each updated search agent Xtemp(i,:) using the surrogates
        %- Approximated value using the global surrogate
        [~,appr_fit(i),~]=globsurrogate_oper_vf7(1,[],[],[],globSurr,[],[],Xtemp(i,:));
        %---------------------------------------------------------------------%
        %- Verify the updated positions approximated fitness vs the personal and global best recorded fitnesses
        if (appr_fit(i)<=S(i)) % test the personal improvement hypothesis of the updated search agent
            %- Given that the approximated fitness/surrogate value is lower than the prior (real) fitness this position/agent is considered promising
            X(i,:)=Xtemp(i,:); % We accept the new position as valid updated position, which is worthy of computing the new real fobj value
            % Always compute and assign  the fitness of all positions/agents using the real obj function to guarantee a correct convergence
            S(i)=feval(fobj,X(i,:)); % Compute the real objective function value using the real obj function, and this new S(i) can be better or worst than the previous
            count_realfit(i)=count_realfit(i)+1; % counts that at least one (this) ith particle real fitness was updated
            %else % Otherwise, i.e., in case the new temporary position don't see an improved approximated fitness/surrogate value
            % we do not accept the new updated position, so we mantain
        end
    end
    %-----------------------------------------------------------------%
    %- Check if during the course of this iteration (iter), the real obj function was used at least one time
    % i.e., the surrogates approximation found a value lower than the recorded Pbest(i), independent of the fact that the correspodent real obj function value was lower or not
    if (sum(count_realfit)==0) % meaning when none of the search agents/positions were evaluated using the real obj function
        %- Calculate the fitness of each particle i using real obj function
        for i=1:pop_size % for all/each of the search agents
            count_realfit(i)=1;
            %- Compute the real objective function value using the real obj function
            X(i,:)=Xtemp(i,:);
            S(i)=feval(fobj,X(i,:)); % Always computes the fitness of all Personal best and Global best positions/agents using the real obj function to guarantee a correct convergence
        end
    end
    %-------------------------------------------------------------------------%
    %- Dual based merit metric to evaluate the new (temporary) positions worthiness (surrogate value+distance to existing particles
    w_merit=wmer_min+((wmer_max-wmer_min)/Max_FEs)*FEs_counter; %weight to assign to each metric's component
    fmerit=merit_func(w_merit,Xtemp,appr_fit,globDB.Pos); % minimum fmerit (close to 0!) is better
    [~,idx_merit]=sort(fmerit,'ascend'); %order of the best updated positions
    %-------------------------------------------------------------------------%
    %- Update global database (DB) based on a dual based criteria
    [globDB,addit_coumtrealfit]=dynDB_managem_v4_A('update',globDB,DB_maxsize,pop_size,fobj,Xtemp,S,idx_merit,count_realfit,iter);
    %     [globDB,addit_coumtrealfit]=dynDB_managem_v4_A('update',globDB,DB_maxsize,pop_size,fobj,X,S,idx_merit,count_realfit,iter); % Alternative where the not accepted positions are not considered
    %-------------------------------------------------------------------------%
    %- Update global surrogate
    [globSurr,~,aux]=globsurrogate_oper_vf7(0,iter,FEs_counter,Max_FEs,[],globDB.Pos,globDB.Fit,Xtemp);
    %     [globSurr,~,aux]=globsurrogate_oper_vf7(0,iter,FEs_counter,Max_FEs,[],globDB.Pos,globDB.Fit,X); % Alternative where the not accepted positions are not considered
    %-------------------------------------------------------------------------%
    %- Update counter and FEs counter (stop criterion)
    iter=iter+1;
    %     Convergence_curve(iter)=DF;
    FEs_counter=FEs_counter+sum(count_realfit)+addit_coumtrealfit;
    Convergence_curve(FEs_counter)=DF;
    surr_trainset_size(iter+1)=aux;
    %-------------------------------------------------------------------------%
    %     figure(1),plot(X(:,1),X(:,2),'*')
    %     figure(2),plot(globPos(:,1),globPos(:,2),'*')
    %     hold on
    %     pause(0.1)
    %     clf
    %     if (FEs_counter>Max_FEs-50)
    %         iter
    %         disp('analyze here')
    %     end
    %-------------------------------------------------------------------------%
end %while
%________________________________________________________%
globDB.FEcounter=FEs_counter;
globDB.uncorrectedConvCurve=Convergence_curve;
%correct the convergenge curve based on the Max_FEs
excess=length(Convergence_curve)-Max_FEs;
repeat_values=find(Convergence_curve==inf);
if (excess) %~=0 %need to remove some repeated values to match the Max_FEs length for the convergence curve
    rem_idxs=randperm(length(repeat_values),excess); %randomly remove the excess values, always a small number
    Convergence_curve(repeat_values(rem_idxs))=[];
end
for i=1:Max_FEs
    if (Convergence_curve(i)==inf && i>1)
        Convergence_curve(i)=Convergence_curve(i-1);   % hold the same value until a new FE arrive
    elseif (Convergence_curve(i)==0 && Convergence_curve(Max_FEs)~=0) %significa que o Convergence_curve(i)=0 foi introduzido arbitrariamente pelo matlab
        Convergence_curve(i)=Convergence_curve(i-1);
    end
end
%________________________________________________________%
elapsed_t=toc(time_start); %store the computational cost in terms of time spent
%________________________________________________________%
end
%_______________________________________________________________________________________________________________________________%
% by Pedro Bento November 2023
% When using please cite and refer this code and its author!