function [bestever,bestP,conv_curv,total_time]=SAEA_RFS_origfPB(lb,ub,ndims,npop,maxFEs,fobj,initPos,initFit)
%% SAEA-RFS Adapted by Pedro Bento November 2023 (retrieved from the repository https://github.com/h-hg/DDEA/blob/master/README.md)
%_____________start function_______________________________%
warning off

tic;

flag = 'cubic';
cycle = 0;

ts_num = 200;  % the number of initial training samples


K = 5;  % old 20! the number of sub-problems
sub_size = 20;  %old 100! maximum dimension of all sub-problems 
nIter = 5; % manter! maximum iterations of sub-problem optimization


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

% Initial population
Parent = initPos;
% generate initial archive

FES = ts_num; %start the FEcounter with the number of started positions for the archive

% the initial population goes to the archive, meaning initialize the remainder

archive=initPos;
init_l=ts_num-npop;

XRRmin = repmat(lu(1, :), init_l, 1);
XRRmax = repmat(lu(2, :), init_l ,1);
archive = [archive; XRRmin + (XRRmax - XRRmin) .* lhsdesign(init_l ,ndims)];
% an extra column is used to store the fitness of each position stored in the archive
archive(1:npop, ndims+1) = initFit;
for i=npop+1:ts_num
    archive(i, ndims+1) = fobj(archive(i,1:ndims));
end

bestever = min(archive(:,ndims+1));
bestP = inf;

array=[linspace(1,ts_num,ts_num)', bestever.*ones(ts_num,1)];        
               
% main loop
while (FES < maxFEs)
    %% form sub-problem and training
    ts = archive;
    [lambda, gamma, dv , sp ] = RFS(sub_size,ts(:, 1:ndims),...
        ts(:, end), flag, K);
    %% Optimization
    [~, index] = min(archive(:, end));
    bestX = archive(index, 1:ndims);
    sub_eval_min = [];

    for i = 1:K
        sub_Parent = Parent(:,dv{i});   % the population of sub-problem
        fit_Parent = RBF_eval(sub_Parent,ts(sp{i}, dv{i}), ...
                lambda{i}, gamma{i}, flag);
        subBound = lu(:, dv{i});    
        % use DE as a basic optimizer for each sub-problem
        for iter = 1:nIter
            sub_child = DE(sub_Parent, bestX(1,dv{i}), subBound);
            % Model evaluation
            fit_child = RBF_eval(sub_child, ts(sp{i}, dv{i}), ...
                lambda{i}, gamma{i}, flag);
            fit_pop = [fit_Parent;fit_child];
            [~,y] = sort(fit_pop);
            com_pop = [sub_Parent;sub_child];
            % Environment selection
            sub_Parent = com_pop(y(1:npop),:);
            fit_Parent = fit_pop(y(1:npop));
        end
     % save the best solution and its estimated value of sub-problem
        sub_eval_min(i,1) = fit_Parent(1);
        sub_eval_min(i,2:length(dv{i})+1) = sub_Parent(1,:); 

        % update the problem population
        Parent(:,dv{i}) = sub_Parent;  
    end

    %% Determine the solution to be evaluated using real funcion
    solution = bestX;
    [min_model,sign] =min(sub_eval_min(:,1));
    temp = RBF_eval(Parent(1,dv{sign}),ts(sp{sign},dv{sign}),lambda{sign},gamma{sign},flag);
    if temp <= min_model
        solution(:,dv{sign}) = Parent(1,dv{sign});
    else
        solution(:,dv{sign}) = sub_eval_min(sign,2:length(dv{sign})+1);
    end
    %% search difference
    diff = setdiff(solution,archive(:,1:ndims),'rows');
    for i=1:size(diff,1)
        diff(i,ndims+1) = fobj(diff(i,1:ndims));
    end
    FES = FES + 1;
    archive = [archive;diff];
    cycle = cycle+1;
    
    % update the global solution
    
%      bestever = min(min(archive(:,ndims+1)),bestever);  %old

    [bestever_pot,index] = min(archive(:,ndims+1));
    if (bestever_pot<bestever)
        bestever=bestever_pot;
        bestP=archive(index,ndims+1);
    else
        bestP=bestX;
    end
    
    array = [array;FES,bestever];
end
total_time=toc;
conv_curv=array(:,2);
end
%_____________end main function_______________________________%
%% Secondary Functions
%--------------------------------------------------------------------------
 function [lambda,gamma,c,r]=RFS(sub_size,S,Y,flag,K)
%
%Input:
%sub-size - maximum dimension of all sub-problems  
%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
%K - the number of sub-problems
%
%Output:
%lambda, gamma - vectors with RBF model parameters
%c - the desicion variables of each sub-problem
%r - the samples used for training models of each sub-problem

[n,m]=size(S);

%% RF training
i = 1;
while i<=K
    M = randperm(m,randperm(sub_size,1));  % form sub-problem
    d = length(M);
    L = randperm(n,2*d);
    S1 = S(L,M);     % training set
    Y1 = Y(L,:);
    [lambda{i},gamma{i}] = RBF(S1,Y1,flag);  % train a model for each sub-problem
    c{i} = M;
    r{i} = L;
    if sum(isnan(lambda{i}))==0 
         i=i+1;
    end  
end

end %function
%--------------------------------------------------------------------------
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);
R = pdist2(S,S,'euclidean');
%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));
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
%--------------------------------------------------------------------------
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

%rho=2;  
[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,nX]=size(X);
end

R=zeros(mX,mS); %compute pairwise distances of points in X and S
R = pdist2(X,S,'euclidean');
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));
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 %function
%--------------------------------------------------------------------------
function U=DE(X,bestX,Bound)
%变异向量用函数mutation(X,bestX,F,mutationStrategy)
%交叉向量用函数crossover(X,V,CR,crossStrategy)
%mutation
%mutationStrategy=1：DE/rand/1,
%mutationStrategy=2：DE/best/1,
%mutationStrategy=3：DE/rand-to-best/1,
%mutationStrategy=4：DE/best/2,
%mutationStrategy=5：DE/rand/2.
%crossover
%crossStrategy=1:binomial crossover
%crossStrategy=2:Exponential crossover

F=0.8;%scaling factor 缩放因子
CR=1;%crossover rate 交叉概率
mutationStrategy=2;%变异策略
crossStrategy=1;%交叉策略

% while Generation<maxIteration
     V=mutation(X,bestX,F,mutationStrategy);
     U=crossover(X,V,CR,crossStrategy,Bound);
% end
end
%--------------------------------------------------------------------------
function V=mutation(X,bestX,F,mutationStrategy)
[NP,D]=size(X);

for i=1:NP
    %在[1 NP]中产生nrandI个互不相等的随机数，且与i皆不相等
    nrandI=5;
    r=randi([1,NP],1,nrandI);
    for j=1:nrandI
    equalr(j)=sum(r==r(j));   % 计数：产生的随机数中相同的个数
    end
    equali=sum(r==i);           % 行号与产生的随机数相同的次数
    equalval=sum(equalr)+equali;     % 相同的值
    while(equalval>nrandI)
        r=randi([1,NP],1,nrandI);
        for j=1:nrandI
        equalr(j)=sum(r==r(j));
        end
        equali=sum(r==i);
        equalval=sum(equalr)+equali;
    end
    
    switch mutationStrategy
        case 1
            %mutationStrategy=1:DE/rand/1;
            V(i,:)=X(r(1),:)+F*(X(r(2),:)-X(r(3),:));
        case 2
            %mutationStrategy=2:DE/best/1;
            V(i,:)=bestX+F*(X(r(1),:)-X(r(2),:));
        case 3
            %mutationStrategy=3:DE/rand-to-best/1;
            V(i,:)=X(i,:)+F*(bestX-X(i,:))+F*(X(r(1),:)-X(r(2),:));
        case 4
            %mutationStrategy=4:DE/best/2;
            V(i,:)=bestX+F*(X(r(1),:)-X(r(2),:))+F*(X(r(3),:)-X(r(4),:));
        case 5
            %mutationStrategy=5:DE/rand/2;
            V(i,:)=X(r(1),:)+F*(X(r(2),:)-X(r(3),:))+F*(X(r(4),:)-X(r(5),:));
        otherwise
            error('没有所指定的变异策略，请重新设定mutationStrategy的值');
    end
   
    
end
end
%--------------------------------------------------------------------------
function U=crossover(X,V,CR,crossStrategy,Bound)
[NP,Dim]=size(X);
switch crossStrategy
    %crossStrategy=1:binomial crossover
    case 1
        for i=1:NP
            jRand=floor(rand*Dim);%由于jRand要在[0,1)*Dim中取值，故而用floor
            for j=1:Dim
                if rand<CR||j==jRand
                    U(i,j)=V(i,j);
                else
                    U(i,j)=X(i,j);
                end     
            end    
        end
    %crossStrategy=2:Exponential crossover
    case 2
        for i=1:NP
            j=floor(rand*Dim);%由于j在[0,1)*Dim中取值，故而用floor
            L=0;
            U=X;
            U(i,j)=V(i,j);
            j=mod(j+1,D);
            L=L+1;
            while(rand<CR&&L<Dim)
                U(i,j)=V(i,j);
                j=mod(j+1,D);
                L=L+1;
            end
        end
        
    otherwise
        error('没有所指定的交叉策略，请重新设定crossStrategy的值'); 
end
%% 越界重置
for i=1:NP
    for j=1:Dim
        while (U(i,j)>Bound(2,j)) || (U(i,j)<Bound(1, j))
            U(i,j)=Bound(1,j)+rand*(Bound(2,j) - Bound(1,j));
        end
    end

end
end
%--------------------------------------------------------------------------