function [bestever,gbpos,conv_curve,time_cost] = TL_SSLPSO_origfPB(lb,ub,max_FEs,ndims,npop,fobj,initPos,init_f)
%%%*********************************************************************************************************%%%
%% Implementation of Truncation-learning-driven surrogate assisted social learning particle swarm optimization for computationally expensive problem (TL-SSLPSO)
%% Haibo Yu , Li Kang , Ying Tan, Chaoli Sun , Jianchao Zeng,
%% Published in:Applied Soft Computing（Volume 97, Part A,?December 2020, 106812)
%%%*********************************************************************************************%%%
%% 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 was 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

% %
warning('off');

CE=zeros(max_FEs,2);     % sampling setting according to FES

time_begin=tic;



%  Invariable C3
c3 =0; % c3 >-1 guarantee convergence
% % %     beta=0.01; c3 = ndims/M*beta;     % epsilon
%%
PL = 1;
%
%initialization
p = initPos;
v = zeros(npop, ndims);

lu = [lb.* ones(1, ndims); ub.* ones(1, ndims)];

gen = 0;
FES = npop;
% % %     rand('seed', sum(100 * clock));

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);

%----------------End initialization------------------
%%
% 通过设置参数group_num 来切换对比算法(RBF-SLPSO vs. TL-SSLPSO)
group_num = 3; % // group_num=1 ---> RBF-SLPSO(S-SLPSO) // *** group_num>1--->Population division-based RBF-SLPSO (TL-SSLPSO) ***// take the value of 3 is of the best according to the sensitivity analysis
%% 灵敏度参数分析：3，5，7，10，15
%% 根据参数灵敏度实验结果，分成3组，整体效果较好

%main loop
while(FES < max_FEs)
    %         fprintf('Iteration: %ndims Fitness evaluation: %ndims Best fitness: %e\n', gen, FES, bestever);% 应该每评价一次，显示一次，而不是每迭代一次显示一次（应放在下面if后面）
    
    %% population segmentation
    if gen >= 1 && group_num~=1  % group_num=1 denotes the comparison algorithm: RBF-SLPSO algorithm
        % Equal division
        num_c=group_num;
        [~,idx]=sort(fitness); % 分割前，将数据按适应值从小到大进行排序
        idx=buffer(idx,ceil(npop/num_c)); % idx 的每一列对应每一分割所包含样本的指标集
        
        p_id=cell(num_c,1);% 各类中种群个体的指标集[初始化]
        pos_mean=cell(num_c,1);% 各类中种群个体集的平均位置
        centroid=cell(num_c,1);% 各类的质心位置
        nc=zeros(1,num_c);% 各类的类内样本数
        cf_mean=zeros(1,num_c);
        for i=1:num_c % divide database into npop levels according to the fitness value
            % Equal division
            k=idx(:,i); k(k==0)=[];
            pos_i=p(k,:);% positions of samples in ith cluster/group
            cf=fitness(k);   cf_mean(i)=mean(cf); % mean fitness value of ith cluster/group
            centroid{i}=mean(pos_i); % the centroid of ith cluster/group
            [~,~,ip]=intersect(pos_i,p,'rows');
            if ~isempty(ip) == 1
                p_id{i}=ip; % record the index of candidate individual belonging to the ith cluster/group
                pos_mean{i}=mean(p(ip,:));
            else
                p_id{i}=[];
                pos_mean{i}=[];
            end
            
            nc(i)=length(p_id{i});% number of samples in ith cluster/group
            
            % count the number of hx and p in ith cluster
            if i == 1 % the group of the highest-level sub-population
                [~,ihh,~]=intersect(hx,pos_i,'rows');
                if ~isempty(ihh)==1
                    num_hx=length(ihh);
                    num_p=nc(i)-num_hx;
                else
                    num_hx=0;
                    num_p=nc(i);
                end
                num_app(gen,:)=[num_hx,num_p];
            end
            
        end
        
        [cf_mean,icf]=sort(cf_mean,'descend'); % 降序排列：according to fitness mean value of each cluster
        
        nc=nc(icf);                  % 排序后对应各类中样本数量
        p_id=p_id(icf);              % 排序后对应各类中个体的指标集
        pos_mean=pos_mean(icf);      % 排序后对应各类中种群个体的平均位置
        centroid=centroid(icf);      % 排序后对应各类样本的平均位置（质心）
        
    end
    
    %% ********************* Phase 2 ************************ behavioral learning
    if gen < 1
        %population sorting
        [fitness,rank] = sort(fitness, 'descend');                 % 按降序（由大到小）进行排列
        p = p(rank,:);
        v = v(rank,:);
        %center position
        center = ones(npop,1)*mean(p);
        %random matrix
        randco1 = rand(npop, ndims);
        randco2 = rand(npop, ndims);
        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);              % 每一个个体的第j维元素贡献学习向量的第j维
        end
        %learning
        %             rand('seed', sum(100 * clock));
        lpmask = repmat(rand(npop,1) < PL, [1 ndims]);
        lpmask(npop,:) = 0;                              % 最优的第m个个体不学习
        v1 =  1*(randco1.*v + randco2.*(pwin - p) + c3*randco3.*(center - p));
        p1 =  p + v1;
        v = lpmask.*v1 + (~lpmask).*v;                % 除最优个体外剩余个体都学习
        p = lpmask.*p1 + (~lpmask).*p;                % 保留最优个体到下一代：第m个个体不学习
    else
        for i = 1:group_num
            p_index=p_id{i}(:,:);
            nn=length(p_index); % number of candidate inviduales in ith cluster/group
            if i == group_num % 排序后最优类种群个体的学习
                if nn == 1 % 只包含一个种群个体
                    continue;% 最优类内的个体不学习
                elseif group_num == 1 % *(comparison algorithm：RBF-SLPSO算法)* 自学习:实验发现，自学习容易破坏潜在最优解，误导种群的搜索
                    sam=[]; fit=[];
                    sam=p(p_index,:);
                    fit=fitness(p_index);
                    v_tmp=v(p_index,:);
                    [~,rank]=sort(fit,'descend');
                    sam=sam(rank,:);
                    fit=fit(rank);
                    %random matrix
                    randco1 = rand(nn, ndims);
                    randco2 = rand(nn, ndims);
                    randco3 = rand(nn, ndims);
                    winidxmask = repmat([1:nn]', [1 ndims]);
                    winidx = winidxmask + ceil(rand(nn, ndims).*(nn - winidxmask));
                    pwin = sam;
                    for j = 1:ndims
                        pwin(:,j) = sam(winidx(:,j),j);              % 每一个个体的第j维元素贡献学习向量的第j维
                    end
                    %learning
                    lpmask = repmat(rand(nn,1) < PL, [1 ndims]);
                    lpmask(nn,:) = 0;                              % 最优的第m个个体不学习
                    v1 =  1*(randco1.*v_tmp + randco2.*(pwin - sam) + c3*randco3.*(ones(nn,1)*pos_mean{i}(:,:) - sam)); % Main strategy:既向优秀个体学习保证算法的收敛性，同时又共享种群的平均信息保证算法的多样性,recall c3=0;
                    p1 =  sam + v1;
                    v_tmp = lpmask.*v1 + (~lpmask).*v_tmp;                % 除最优个体外剩余个体都学习
                    sam = lpmask.*p1 + (~lpmask).*sam;                % 保留最优个体到下一代：第m个个体不学习
                    
                    idd = p_index;
                    idd = idd(rank);
                    p(idd,:)=sam;
                    v(idd,:)=v_tmp;
                    
                    continue;
                else
                    continue
                end
                
            end
            
            %% 需要找到比当前类层次高的所有类 (Demonstrator selection)
            id_demo=i+1:group_num; % 包含DB样本和种群个体的所有高层类 (Main strategy)
            num_democluster=length(id_demo);
            
            for j = 1:ndims
                %% competitive selection ---> RBF-SLPSO-division (TL-SSLPSO-B)
                %                     if num_democluster == 1
                %                         dim_cluster = num_c;
                %                     else
                %                         id=randperm(num_democluster);
                %                         dim_cluster=id_demo(id(1:2)); % j 维度上 demonstrator *类指标*
                %                         % 不同于level-based pso中比较高层类中随机对（random pair）个体的适应值,这里我们选取比较高层类间的适应值均值，选取较优类中的个体作为当前个体的demonstrators
                %                         if cf_mean(dim_cluster(1)) < cf_mean(dim_cluster(2))
                %                             dim_cluster=dim_cluster(1);
                %                         else
                %                             dim_cluster=dim_cluster(2);
                %                         end
                %                     end
                
                %% random selection ---> RBF-SLPSO-division-random (TL-SSLPSO-R)
                id=randi(num_democluster); % randomly select one---> for comparison
                dim_cluster=id_demo(id);
                %%
                demo_p=p(p_id{dim_cluster}(:,:),:);
                demonstrator=demo_p; % (Main strategy)：demonstrator太多，考虑选择性的保留若干代表性的历史样本
                pwin(p_index,j) = demonstrator(randi(size(demonstrator,1),1,size(p_index,1)),j);
            end
            
            %learning
            lpmask = repmat(rand(nn,1) < PL, [1 ndims]);
            %random matrix
            randco1 = rand(nn, ndims);
            randco2 = rand(nn, ndims);
            randco3 = rand(nn, ndims);
            v1 =  1*(randco1.*v(p_index,:) + randco2.*(pwin(p_index,:)  - p(p_index,:)) + c3*randco3.*(ones(nn,1)*pos_mean{i}(:,:) - p(p_index,:))); % Main strategy
            
            p1 =  p(p_index,:) + v1;
            v(p_index,:) = lpmask.*v1 + (~lpmask).* v(p_index,:);                % 除最优个体外剩余个体都学习
            p(p_index,:) = lpmask.*p1 + (~lpmask).* p(p_index,:);                % 保留最优个体到下一代：第m个个体不学习
        end
    end
    %boundary
    for i = 1:npop
        p(i,:) = max(p(i,:), lu(1,:)); % 下界
        p(i,:) = min(p(i,:), lu(2,:)); % 上界
    end
    
    gen = gen + 1;
    
    % % %         % compute swarm diversity
    % % %         p_mean=mean(p);
    % % %         distance=real(sqrt(p_mean.^2*ones(size(p'))+ones(size(p_mean))*(p').^2-2*p_mean*(p')));
    % % %         diversity(gen)=max(distance);
    
    %% RBF-modeling and evaluation
    % select training samples% 包围当前搜索种群的 DB 样本 ------> 局部模型（local）
    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);
    
    % 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';% model values
    
    %% *** Expensive evaluation: individual-based evolution control **************** %%
    %% 计算适应度估值优于当前最优的个体的真实适应值，否则计算当前估值最优个体的真实适应值
    %% 实验发现，受限于代理模型精度及问题解空间特性的影响，算法的寻优过程（若干次迭代搜索）存在受困于某些真实的历史最优区域的风险。
    %% 尽管如此，算法结合分组学习机制，有能力跳出这些影响域,换言之，分组学习机制间接减小了算法陷入模型局部极值的概率，因为其保留了部分局部最优样本，这些样本不进行学习
    %% Two alternative infill sampling criteria
    pid = find(fitness < bestever);
    p_tmp = p(pid,:);
    f_tmp = fitness(pid);
    if ~isempty(p_tmp) == 1
        for i = 1 : size(p_tmp,1)
            [~,ih,~]=intersect(hx,p_tmp(i,:),'rows');
            if ~isempty(ih)==1
                f_tmp(i)=hf(ih);
            else
                
                f_tmp(i)=fobj(p_tmp(i,:));
                FES=FES+1;
                if FES <= max_FEs
                    if (f_tmp(i)<CE(FES-1,2)) % melhoria na conv curve é aceite, em teoria equivale a sbesty<bestever
                        CE(FES,:)=[FES,f_tmp(i)];
                    else
                        CE(FES,:)=[FES,CE(FES-1,2)]; % se nao repete o valor já encontrado
                    end
                end
                hx=[hx;p_tmp(i,:)];   hf=[hf,f_tmp(i)];                  % 更新历史数据库
            end
            [bestever,ib] = min([f_tmp(i), bestever]);             % 更新全局最优
            if ib==1
                gbpos=p_tmp(i,:);
            end
        end
        fitness(pid) = f_tmp;
    else
        %% 计算适应度估值最优的个体的真实适应值
        [~,idx]=sort(fitness);
        p_app=p(idx,:); f_app=fitness(idx);
        [~,~,ip]=intersect(hx,p_app,'rows');
        p_app(ip,:)=[]; 
        f_app(ip)=[]; % delete history samples
        
        sbest_pos=p_app(1,:);% performance ----> under limited computational cost, we prefer fast convergence performance
        
        sbesty=fobj(sbest_pos);
        FES=FES+1;
        
        if FES <= max_FEs
            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
        end
        
        hx=[hx;sbest_pos];   hf=[hf,sbesty];                  % 更新历史数据库
        [bestever,ib] = min([sbesty, bestever]);             % 更新全局最优
        if ib==1
            gbpos=sbest_pos;
        end
        [~,ip,~]=intersect(p,p_app(1,:),'rows');
        fitness(ip)=sbesty;
    end
    
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
%______________________________________________________________________________
