function [bestever,gbpos,conv_curve,time_cost] = GORS_SSLPSO_origfPB(lb,ub,max_FEs,ndims,npop,fobj,initPos,init_f)
%%%*********************************************************************************************************%%%
%% Implementation of A Generation-based Optimal Restart Strategy for Surrogate-assisted SLPSO (GORS-SSLPSO) adapted to one function by PedroBento
%% H. Yu, Y. Tan, C. Sun, J. Zeng, A generation-based optimal restart strategy for surrogate-assisted
%% social learning particle swarm optimization, Knowledge-Based Systems, (2018).
%%%*********************************************************************************************%%%
%% This paper and this code should be referenced whenever they are used to
%% generate results for the user's own research.
%%%*********************************************************************************************%%%
%% This matlab code was written by Haibo Yu and Adapted by Pedro Bento November 2023 (retrieved from the repository https://github.com/h-hg/DDEA/blob/master/README.md)
%% Please refer with all questions, comments, bug reports, etc. to tyustyuhaibo@126.com
% %
% rng('default');
% rng('shuffle');


warning('off');

CE=zeros(max_FEs,2);

time_begin=tic;

%---------------Initialization-----------------
%parameter setting
c3 = 0;
PL = ones(npop,1);

%initialization
p = initPos;
v = zeros(npop, ndims);
lu = [lb* ones(1, ndims); ub* ones(1, ndims)];
  
gen = 0;
FES = npop;  

fitness=init_f';

hx=p;   
hf=fitness;

[bestever,id] = min(fitness);
gbpos=p(id,:);

CE(1:npop,1)=1:1:npop;
CE(1:npop,2)=bestever.*ones(npop,1);


num_gen=30;% 1(SSLPSO), 10, 20, 30（Standard.GORS-SSLPSO）, 40, 50, 60, 80

%main loop
while(FES < max_FEs)
    
    if mod(gen, num_gen) == 0 % only meeting this condition, the top-ranking npop data are used
        
%         fprintf('Iteration: %d Fitness evaluation: %d Best fitness: %e\n', gen, FES, bestever);
        
        [~,idx]=sort(fitness);
        p_app=p(idx,:); 
        f_app=fitness(idx);
        [~,~,ip]=intersect(hx,p_app,'rows');
        p_app(ip,:)=[];
        f_app(ip)=[];
        if ~isempty(p_app)==1
            sbest_pos=p_app(1,:);
            
            sbesty=fobj(sbest_pos);
            FES=FES+1;
            
            if (sbesty<CE(FES-1,2)) % melhoria na conv curve é aceite, em teoria equivale a sbesty<bestever
                CE(FES,:)=[FES,sbesty];
            else
                CE(FES,:)=[FES,CE(FES-1,2)]; % se nao repete o valor já encontrado
            end
            

            hx=[hx;sbest_pos];   hf=[hf,sbesty];
            [bestever,ib] = min([sbesty, bestever]);
            if ib==1
                gbpos=sbest_pos;
            end
        end
        
        % select the best npop sample points to form the population
        [~,idx]=sort(hf);   idx=idx(1:npop);
        p=hx(idx,:);        fitness=hf(idx);
        
    end
    
    %% ********************* Phase 2 ************************ behavioral learning
    %population sorting
    [~,rank] = sort(fitness, 'descend');
    p = p(rank,:);
    v = v(rank,:);
    %center position
    center = ones(npop,1)*mean(p);
    %random matrix
    %         rand('seed', sum(100 * clock));
    randco1 = rand(npop, ndims);
    %         rand('seed', sum(100 * clock));
    randco2 = rand(npop, ndims);
    %         rand('seed', sum(100 * clock));
    randco3 = rand(npop, ndims);
    winidxmask = repmat([1:npop]', [1 ndims]);
    winidx = winidxmask + ceil(rand(npop, ndims).*(npop - winidxmask));
    %     winidx = npop - floor(0.5*rand(npop, ndims).*(npop - winidxmask));
    pwin = p;
    for j = 1:ndims
        pwin(:,j) = p(winidx(:,j),j);
    end
    %learning
    %         rand('seed', sum(100 * clock));
    lpmask = repmat(rand(npop,1) < PL, [1 ndims]);
    lpmask(npop,:) = 0;
    v1 =  1*(randco1.*v + randco2.*(pwin - p) + c3*randco3.*(center - p));
    p1 =  p + v1;
    v = lpmask.*v1 + (~lpmask).*v;
    p = lpmask.*p1 + (~lpmask).*p;
    %boundary
    for i = 1:npop
        p(i,:) = max(p(i,:), lu(1,:));
        p(i,:) = min(p(i,:), lu(2,:));
    end
    gen = gen + 1;
    
    %% ********************* RBF modeling ************************
    % Local model for 50D / 100D
    %         NS=2*(ndims+1);
    %         phdis=real(sqrt(p.^2*ones(size(hx'))+ones(size(p))*(hx').^2-2*p*(hx')));
    %         [~,sidx]=sort(phdis,2);
    %         nidx=sidx; nidx(:,NS+1:end)=[];
    %         nid=unique(nidx);
    %         trainx=hx(nid,:);
    %         trainf=hf(nid);
    
    % Global model for 10D / 20D / 30D
    trainx=hx;  trainf=hf;
    
    % radial basis function interpolation----(RBF-interpolation)
    flag='cubic';
    [lambda, gamma]=RBF(trainx,trainf',flag);
    FUN=@(x) RBF_eval(x,trainx,lambda,gamma,flag);
    %%
    fitness=FUN(p);     fitness=fitness';
end

conv_curve=CE(:,2);
time_cost=toc(time_begin);
end
% end main function!
%% Auxiliary functions
%______________________________________________________________________________
function [lambda, gamma]=RBF(S,Y,flag)
% function [lambda, gamma, rho]=RBF(S,Y,flag)   % For Gaussian basis function

%computes the parameters of the radial basis function interpolant
%--------------------------------------------------------------------------
%Copyright (c) 2012 by Juliane Mueller
%
% This file is part of the surrogate model module toolbox.
%
%--------------------------------------------------------------------------
%Author information
%Juliane Mueller
%Tampere University of Technology, Finland
%juliane.mueller2901@gmail.com
%--------------------------------------------------------------------------
%
%Input:
%S - Sample site matrix (m x d), m=number of sample sites, d=dimension
%Y - Objective function values corresponding to points in S
%flag - string determining what type of RBF model should be used
%
%Output:
%lambda, gamma - vectors with RBF model parameters
%--------------------------------------------------------------------------


[m,n]=size(S);
P=[S,ones(m,1)];

R=zeros(m,m);
for ii = 1:m
    for jj = ii:m
        R(ii,jj)=sum((S(ii,:)-S(jj,:)).^2,2);
        R(jj,ii)=R(ii,jj);
    end
end
R=sqrt(R);
% rho=2;

% Sd=real(sqrt(S.^2*ones(size(S'))+ones(size(S))*(S').^2-2*S*(S')));    % For Gaussian basis function
% rho=max(max(Sd))/(n*m)^(1/n);     % For Gaussian basis function

if strcmp(flag,'cubic') %cubic RBF
    Phi= R.^3;
elseif strcmp(flag,'TPS') %thin plate spline RBF
    R(R==0)=1;
    Phi=R.^2.*log(R);
elseif strcmp(flag,'linear') %linear RBF
    Phi=R;
% elseif strcmp(flag,'Gaussian')
%     Phi=exp(-R.^2/rho^2);
% elseif strcmp(flag,'multiquad')
%     Phi=sqrt((R.^2 + rho^2).^3);
% elseif strcmp(flag,'invmultiquad')
%     Phi=1./sqrt(R.^2 + rho^2);
end

A=[Phi,P;P', zeros(n+1,n+1)];
RHS=[Y;zeros(n+1,1)];
params=A\RHS;
lambda=params(1:m);
gamma=params(m+1:end);

end
%______________________________________________________________________________
function Yest=RBF_eval(X,S,lambda,gamma,flag)
% function Yest=RBF_eval(X,S,lambda,gamma,rho,flag)   % For Gaussian basis function

%--------------------------------------------------------------------------
%Copyright (c) 2012 by Juliane Mueller
%
% This file is part of the surrogate model module toolbox.
%
%--------------------------------------------------------------------------
%Author information
%Juliane Mueller
%Tampere University of Technology, Finland
%juliane.mueller2901@gmail.com
%--------------------------------------------------------------------------
%
%input:
%X are points where function values should be calculated
%S are points where the function values are known
%lambda parameter vector
%gamma contains the parameters of the optional polynomial tail
%flag is a string indicating which RBF to be used
%output: 
%the estimated function values at the points in X
%--------------------------------------------------------------------------

[mX,nX]=size(X); %dimensions of the points where function value should be predicted
[mS,nS]=size(S); %dimesnions of sample sites taken so far 
if nX~=nS %check that both matrices are of the right shape
    X=X';
    [mX,~]=size(X);
end

R=zeros(mX,mS); %compute pairwise distances of points in X and S
for ii = 1:mX
    for jj = 1:mS
        R(ii,jj)=norm(X(ii,:)-S(jj,:));
    end
end

if strcmp(flag,'cubic')
    Phi= R.^3;
elseif strcmp(flag,'TPS')
    R(R==0)=1;
    Phi=R.^2.*log(R);
elseif strcmp(flag,'linear')
    Phi=R;
% elseif strcmp(flag,'Gaussian')
%     Phi=exp(-R.^2/rho^2);
% elseif strcmp(flag,'multiquad')
%     Phi=sqrt((R.^2 + rho^2).^3);
% elseif strcmp(flag,'invmultiquad')
%     Phi=1./sqrt(R.^2 + rho^2);
end
    
Yest1=Phi*lambda; %first part of response surface
Yest2=[X,ones(mX,1)]*gamma; %optional polynomial tail
Yest=Yest1+Yest2; %predicted function value

end
%______________________________________________________________________________