function [lb,ub,dim,fobj] = Get_Functions_details3(F)
%% Function with the tailored set of mathematical benchmark test functions for the surrogate problem (including the common F1-F7 functions with other number attributed)
% The function handle, i.e., the required objective function and the respective problem bounds are returned by this function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Adapted by Pedro Bento November 2023 
% When using please cite and refer this code and its author!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch F
    case 'F1'
        fobj = @F1;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F2'
        fobj = @F2;
        lb=-10;
        ub=10;
        dim=30;
        
    case 'F3'
        fobj = @F3;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F4'
        fobj = @F4;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F5'
        fobj = @F5;
        lb=-30;
        ub=30;
        dim=30;
        
    case 'F6'
        fobj = @F6;
        lb=-100;
        ub=100;
        dim=30;
        
    case 'F7'
        fobj = @F7;
        lb=-1.28;
        ub=1.28;
        dim=30;
        
    case 'F8'
        fobj = @F8;
        lb=-500;
        ub=500;
        dim=30;
        
    case 'F9'
        fobj = @F9;
        lb=-5.12;
        ub=5.12;
        dim=30;
        
    case 'F10'
        fobj = @F10;
        lb=-32;
        ub=32;
        dim=30;
        
    case 'F11'
        fobj = @F11;
        lb=-600;
        ub=600;
        dim=30;
        
    case 'F12'
        fobj = @F12;
        lb=-50;
        ub=50;
        dim=30;
        
    case 'F13'
        fobj = @F13;
        lb=-50;
        ub=50;
        dim=30;
        
    case 'F14' % newly added Ellipsoid Function
        fobj = @F14;
        lb=-100;
        ub=100;
        dim=30;
    
    case 'F15' % newly added Shifted Rotated Rastrigin's (F10 in CEC 2005)
        fobj = @F15;
        lb=-5;
        ub=5;
        dim=30;
        
    case 'F16' % newly added Rotated Hybrid Composition (F16 in CEC 2005)
        fobj = @F16;
        lb=-5;
        ub=5;
        dim=30;
        
    case 'F17' % newly added Ellipsoid Function
        fobj = @F17;
        lb=-5;
        ub=5;
        dim=30;
end
end
%------------------------------------------------------------
%% Defined functions
% F1 -> Sphere

function o = F1(x)
    o=sum(x.^2);
end

% F2 -> Schwefel 2.22

function o = F2(x)
    o=sum(abs(x))+prod(abs(x));
end

% F3 -> Schwefel 1.2

function o = F3(x)
    dim=size(x,2);
    o=0;
    for i=1:dim
        o=o+sum(x(1:i))^2;
    end
end

% F4 -> Schwefel 2.21

function o = F4(x)
    o=max(abs(x));
end

% F5 -> Rosenbrock

function o = F5(x)
    dim=size(x,2);
    o=sum(100*(x(2:dim)-(x(1:dim-1).^2)).^2+(x(1:dim-1)-1).^2);
end

% F6 -> Step 

function o = F6(x)
    o=sum(abs((x+.5)).^2);
end

% F7 -> Quartic

function o = F7(x)
    dim=size(x,2);
    o=sum([1:dim].*(x.^4))+rand;
end

% F8 -> Schwefel 2.26

function o = F8(x)
    o=sum(-x.*sin(sqrt(abs(x))));
end

% F9 -> Rastrigin

function o = F9(x)
    dim=size(x,2);
    o=sum(x.^2-10*cos(2*pi.*x))+10*dim;
end

% F10 -> Ackley

function o = F10(x)
    dim=size(x,2);
    o=-20*exp(-.2*sqrt(sum(x.^2)/dim))-exp(sum(cos(2*pi.*x))/dim)+20+exp(1);
end

% F11 -> Griewank

function o = F11(x)
    dim=size(x,2);
    o=sum(x.^2)/4000-prod(cos(x./sqrt([1:dim])))+1;
end

% F12 -> Penalized (1)

function o = F12(x)
    dim=size(x,2);
    o=(pi/dim)*(10*((sin(pi*(1+(x(1)+1)/4)))^2)+sum((((x(1:dim-1)+1)./4).^2).*...
    (1+10.*((sin(pi.*(1+(x(2:dim)+1)./4)))).^2))+((x(dim)+1)/4)^2)+sum(Ufun(x,10,100,4));
end

% F13 -> Penalized (2)

function o = F13(x)
    dim=size(x,2);
    o=.1*((sin(3*pi*x(1)))^2+sum((x(1:dim-1)-1).^2.*(1+(sin(3.*pi.*x(2:dim))).^2))+...
    ((x(dim)-1)^2)*(1+(sin(2*pi*x(dim)))^2))+sum(Ufun(x,5,100,4));
end

% F14 -> Ellipsoid (from https://www.mdpi.com/2076-3417/11/11/5053)
 
function o = F14(x)
    dim=size(x,2);
    outer = 0;

    for ii = 1:dim
        xj = x(ii);
        outer = outer+ii*(xj^2);
    end

    o = outer;
end

% F15 -> Shifted Rotated Rastrigin's (from F10 in CEC2005)
 
function o = F15(x)
    [ps,Dims]=size(x);
    load cec2005_matlab/rastrigin_func_data % load from the cec2005_matlab subfolder (original matlab code for cec2005 funcs)
    if length(o)>=Dims
         o=o(1:Dims);
    else
         o=-5+10*rand(1,Dims);
    end
    c=2;
    if Dims==2
        load cec2005_matlab/rastrigin_M_D2.mat
    elseif Dims==10
        load cec2005_matlab/rastrigin_M_D10
    elseif Dims==30
        load cec2005_matlab/rastrigin_M_D30
    elseif Dims==50
        load cec2005_matlab/rastrigin_M_D50
    else
        M=rot_matrix(Dims,c);
    end
    x=x-repmat(o,ps,1);
    x=x*M;
    o_partial=sum(x.^2-10.*cos(2.*pi.*x)+10,2);
    load cec2005_matlab/fbias_data;
    o=o_partial+f_bias(10);
end

% F16 -> Rotated Hybrid Composition Function 1(from F16 in CEC2005)

function o=F16(x)
    [~,Dims]=size(x);
    fun_num=10;
    load cec2005_matlab/hybrid_func1_data % saved the predefined optima
    if length(o(1,:))>=Dims
         o=o(:,1:Dims);
    else
         o=-5+10*rand(fun_num,Dims);
    end
    func.f1=@(x) frastrigin(x);
    func.f2=@(x) frastrigin(x);
    func.f3=@(x) fweierstrass(x);
    func.f4=@(x) fweierstrass(x);
    func.f5=@(x) fgriewank(x);
    func.f6=@(x) fgriewank(x);
    func.f7=@(x) fackley(x);
    func.f8=@(x) fackley(x);
    func.f9=@(x) fsphere(x);
    func.f10=@(x) fsphere(x);
    bias=((1:fun_num)-1).*100;
    sigma=ones(1,fun_num);
    lamda=[1; 1; 10; 10; 5/60; 5/60; 5/32; 5/32; 5/100; 5/100];
    lamda=repmat(lamda,1,Dims);
    c=[2,2,2,2,2,2,2,2,2,2,2];
    if Dims==2
        load cec2005_matlab/hybrid_func1_M_D2
    elseif Dims==10
        load cec2005_matlab/hybrid_func1_M_D10
    elseif Dims==30
        load cec2005_matlab/hybrid_func1_M_D30
    elseif Dims==50
        load cec2005_matlab/hybrid_func1_M_D50
    else
        for i=1:fun_num
            eval(['M.M' int2str(i) '=rot_matrix(Dims,c(i));']);
        end
    end
    o_partial=hybrid_composition_func(x,fun_num,func,o,sigma,lamda,bias,M);
    load cec2005_matlab/fbias_data;
    o=o_partial+f_bias(16);
end

% F17 -> Rotated Hybrid Composition Function 2 with a Narrow Basin for the Global Optimum (from F19 in CEC2005)

function o=F17(x)
    [~,Dims]=size(x);
    fun_num=10;
    load cec2005_matlab/hybrid_func2_data % saved the predefined optima
    if length(o(1,:))>=Dims
         o=o(:,1:Dims);
    else
         o=-5+10*rand(fun_num,Dims);
    end
    o(10,:)=0;
    func.f1=@(x) fackley(x);
    func.f2=@(x) fackley(x);
    func.f3=@(x) frastrigin(x);
    func.f4=@(x) frastrigin(x);
    func.f5=@(x) fsphere(x);
    func.f6=@(x) fsphere(x);
    func.f7=@(x) fweierstrass(x);
    func.f8=@(x) fweierstrass(x);
    func.f9=@(x) fgriewank(x);
    func.f10=@(x) fgriewank(x);
    bias=((1:fun_num)-1).*100;
    sigma=[0.1 2 1.5 1.5 1 1 1.5 1.5 2 2];
    lamda=[0.1*5/32; 5/32; 2*1; 1; 2*5/100; 5/100; 2*10; 10; 2*5/60; 5/60];
    lamda=repmat(lamda,1,Dims);
    c=[2 3 2 3 2 3 20 30 200 300];
    if Dims==2
        load cec2005_matlab/hybrid_func2_M_D2
    elseif Dims==10
        load cec2005_matlab/hybrid_func2_M_D10
    elseif Dims==30
        load cec2005_matlab/hybrid_func2_M_D30
    elseif Dims==50
        load cec2005_matlab/hybrid_func2_M_D50
    else
        for i=1:fun_num
            eval(['M.M' int2str(i) '=rot_matrix(Dims,c(i));']);
        end
    end
    o_partial=hybrid_composition_func(x,fun_num,func,o,sigma,lamda,bias,M);
    load cec2005_matlab/fbias_data;
    o=o_partial+f_bias(19);
end
%------------------------------------------------------------
%% Auxiliary Functions (called by the "main" functions
% Ufun-> used in F12-Penalized (1) and F13-Penalized (2)
function f=Ufun(x,a,k,m)
    f=k.*((x-a).^m).*(x>a)+k.*((-x-a).^m).*(x<(-a));
end
% rot_matrix-> used in F15-Shifted Rotated Rastrign's Function 
function M=rot_matrix(D,c)
  A=normrnd(0,1,D,D);
  P=cGram_Schmidt(A);
  A=normrnd(0,1,D,D);
  Q=cGram_Schmidt(A);
  u=rand(1,D);
  D=c.^((u-min(u))./(max(u)-min(u)));
  D=diag(D);
  M=P*D*Q;
end
% single classical functions for use in the composition/compositefunctions: F16 and F17
function f=fsphere(x)
    [ps,D]=size(x);
    f=sum(x.^2,2);
end
function f=fgriewank(x)
    [ps,D]=size(x);
    f=1;
    for i=1:D
        f=f.*cos(x(:,i)./sqrt(i));
    end
    f=sum(x.^2,2)./4000-f+1;
end
function f=fackley(x)
    [ps,D]=size(x);
    f=sum(x.^2,2);
    f=20-20.*exp(-0.2.*sqrt(f./D))-exp(sum(cos(2.*pi.*x),2)./D)+exp(1);
end
function f=frastrigin(x)
    [ps,D]=size(x);
    f=sum(x.^2-10.*cos(2.*pi.*x)+10,2);
end
function [f]=fweierstrass(x)
    [ps,D]=size(x);
    x=x+0.5;
    a = 0.5;
    b = 3;
    kmax = 20;
    c1(1:kmax+1) = a.^(0:kmax);
    c2(1:kmax+1) = 2*pi*b.^(0:kmax);
    f=0;
    c=-w(0.5,c1,c2);
    for i=1:D
    f=f+w(x(:,i)',c1,c2);
    end
    f=f+c*D;
end
% hybrid_composition_func computes the objective value in the composition/composite functions: F16 and F17
function fit=hybrid_composition_func(x,fun_num,func,o,sigma,lamda,bias,M)
    [ps,D]=size(x);
    for i=1:fun_num
        oo=repmat(o(i,:),ps,1);
        weight(:,i)=exp(-sum((x-oo).^2,2)./2./(D*sigma(i)^2));
        %[-sum((x-oo).^2,2), 2*D*sigma(i)^2, weight(:,i) ]
    end
    [tmp,~]=sort(weight,2);
    for i=1:ps
        weight(i,:)=(weight(i,:)==tmp(i,fun_num)).*weight(i,:)+(weight(i,:)~=tmp(i,fun_num)).*(weight(i,:).*(1-tmp(i,fun_num).^10));
    end
    weight=weight./repmat(sum(weight,2),1,fun_num);
    fit=0;
    for i=1:fun_num
        oo=repmat(o(i,:),ps,1);
        eval(['f=feval(func.f' int2str(i) ',((x-oo)./repmat(lamda(i,:),ps,1))*M.M' int2str(i) ');']);
        x1=5*ones(1,D);
        eval(['f1=feval(func.f' int2str(i) ',(x1./lamda(i,:))*M.M' int2str(i) ');']);
        %eval(['(x1./lamda(i,:))*M.M' int2str(i)])
        fit1=2000.*f./f1;
        %[fit1 f1 weight(:,i)]
        %weight(:,i).*(fit1+bias(i))
        fit=fit+weight(:,i).*(fit1+bias(i));
    end
end
% classical Gram Schmid (auxiliary for rot_matrix function)
 function [q,r] = cGram_Schmidt (A)
% computes the QR factorization of $A$ via
% classical Gram Schmid 
    % 
     [n,m] = size(A); 
     q = A;    
     for j=1:m
         for i=1:j-1 
             r(i,j) = q(:,j)'*q(:,i);
         end
         for i=1:j-1   
           q(:,j) = q(:,j) -  r(i,j)*q(:,i);
         end
         t =  norm(q(:,j),2 ) ;
         q(:,j) = q(:,j) / t ;
         r(j,j) = t  ;
     end
 end
 %(auxiliary for hybrid_composition_func function)
 function y = w(x,c1,c2) 
    y = zeros(length(x),1);
    for k = 1:length(x)
        y(k) = sum(c1 .* cos(c2.*x(:,k)));
    end
end
%_______________________________________________________________________________________________________________________________%
% by Pedro Bento November 2023
% When using please cite and refer this code and its author!