
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% A new firefly algorithm with improved global exploration and
%%%    convergence with application to engineering optimization
%%%                 Decision Analytics Journal 2022
%%%           Mojtaba Ghasemi,Soleiman kadkhoda Mohammadi,
%%%       Mohsen Zare, Seyedali Mirjalili, Milad Gil,Rasul Hemmati
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear all;

for trr=1:30
    func_num=trr;
    kkk=0;
    for iijj=1:25
        
        fhd=str2func('cec14_func');%if you want to use form  Cost Functions CEC 2014
        
        nVar=30;                 % Number of Decision Variables
        
        VarSize=[1 nVar];       % Decision Variables Matrix Size
        
        %% Firefly Algorithm Parameters
        nPop=30;            % Number of Fireflies (Swarm Size)
        MaxIt=300000;         % Maximum Number of Iterations
        VarMin=-100;             % Decision Variables Lower Bound
        VarMax= -VarMin;             % Decision Variables Upper Bound
        
        %% Firefly Algorithm Parameters
        
        gamma=1;            % Light Absorption Coefficient
        
        beta0=2;            % Attraction Coefficient Base Value
        
        alpha=0.2;          % Mutation Coefficient
        
        alpha_damp=0.98;    % Mutation Coefficient Damping Ratio
        
        delta=0.05*(VarMax-VarMin);     % Uniform Mutation Range
        
        m=2;
        
        if isscalar(VarMin) && isscalar(VarMax)
            dmax = (VarMax-VarMin)*sqrt(nVar);
        else
            dmax = norm(VarMax-VarMin);
        end
        
        %% Initialization
        
        % Empty Firefly Structure
        firefly.Position=[];
        firefly.Cost=[];
        
        % Initialize Population Array
        pop=repmat(firefly,nPop,1);
        
        % Initialize Best Solution Ever Found
        BestSol.Cost=inf;
        
        % Create Initial Fireflies
        for i=1:nPop
            pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
            pop(i).Cost=feval(fhd,pop(i).Position',func_num);%if you want to use form  Cost Functions CEC 2014
            if pop(i).Cost<=BestSol.Cost
                BestSol=pop(i);
            end
            CCost(i)=BestSol.Cost-func_num*100;
        end
        
        
        BestCost=zeros(MaxIt,1);
        
        %% Firefly Algorithm Main Loop
        
        it=nPop;
        % for it=1:MaxIt
        while it<=MaxIt
            
            newpop=repmat(firefly,nPop,1);
            for i=1:nPop
                newpop(i).Cost = inf;
                
                for j=1:3:nPop
                    
                    rij=norm(pop(i).Position-pop(j).Position)/dmax;
                    beta=beta0*exp(-gamma*rij^m);
                    e=delta*unifrnd(-1,+1,VarSize);
                    
                    rij2=norm(pop(i).Position-pop(j+1).Position)/dmax;
                    beta2=beta0*exp(-gamma*rij2^m);
                    
                    rij3=norm(pop(i).Position-pop(j+2).Position)/dmax;
                    beta3=beta0*exp(-gamma*rij3^m);
                    
                    ER1=1;
                    ER2=1;
                    ER3=1;
                    
                    A1=sign(pop(i).Cost-pop(j).Cost);
                    A2=sign(pop(i).Cost-pop(j+1).Cost);
                    A3=sign(pop(i).Cost-pop(j+2).Cost);
                    newsol.Position = pop(i).Position+ER1*1*A1*rand(VarSize).*(pop(j).Position-pop(i).Position)+ER2*1*A2*rand(VarSize).*(pop(j+1).Position-pop(i).Position)+ER3*1*A3*rand(VarSize).*(pop(j+2).Position-pop(i).Position) ...
                        + alpha*e;
                    newsol.Position=max(newsol.Position,VarMin);
                    newsol.Position=min(newsol.Position,VarMax);
                    
                    %                 newsol.Cost=CostFunction(newsol.Position);
                    newsol.Cost=feval(fhd,newsol.Position',func_num);%if you want to use form  Cost Functions CEC 2014
                    if newsol.Cost <= newpop(i).Cost
                        newpop(i) = newsol;
                        if newpop(i).Cost<=BestSol.Cost
                            BestSol=newpop(i);
                        end
                    end
                    it=it+1;
                    CCost(it)=BestSol.Cost-func_num*100;
                    
                end
            end
            
            
            pop=[pop
                newpop];
            
            % Sort
            [~, SortOrder]=sort([pop.Cost]);
            pop=pop(SortOrder);
            
            pop=pop(1:nPop);
            
            % Store Best Cost Ever Found
            BestCost(it)=BestSol.Cost-func_num*100;
            
            % Show Iteration Information
            disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
            
            % Damp Mutation Coefficient
            alpha = alpha*alpha_damp;
            
        end
        
        Cost_Rsult(1, iijj)=BestSol.Cost-func_num*100;
        Rsult(iijj,:)=CCost(1,1:300000);
        %  plot(log(BestCost),'r--');   hold on;
    end
    trr
    Mean(trr)=mean(Cost_Rsult);
    Best(trr)=min(Cost_Rsult);
    Std(trr)=std(Cost_Rsult);
    hold on
    plot(log(mean(Rsult)),'k--','LineWidth',4); hold on
end
