%--------------------------------------------------------------------------
% MRFO code v1.0.
% Developed in MATLAB R2011b
% The code is based on the following papers.
% W. Zhao, Z. Zhang and L. Wang, Manta ray foraging optimization: An
% effective bio-inspired optimizer for engineering applications.
% Engineering Applications of Artifical Intelligence (2019),
% https://dio.org/10.1016/j.engappai.2019.103300.
% -------------------------------------------------------------------------
% FunIndex: Index of function.
% MaxIt: The maximum number of iterations.
% PopSize: The size of population.
% PopPos: The position of population.
% PopFit: The fitness of population.
% Dim: The dimensionality of prloblem.
% Alpha: The weight coefficient in chain foraging.
% Beta: The weight coefficient in cyclone foraging.
% S: The somersault factor.
% BestF: The best fitness corresponding to BestX. 
% HisBestFit: History best fitness over iterations. 
% Low: The low bound of search space.
% Up: The up bound of search space.

function [optimalChart]=HMRFO_func(fhd,D,max_nfes,Xmin,Xmax,func,FW,PR)

%          [Xmin,Xmax,D]=FunRange(F_index); 
    nfes = 0;
    optimalChart = [];
    optimalChart_interval = 100;
    optimum = func * 100.0;
    val_2_reach = 10^(-8);
    rand('seed', sum(100 * clock));
    nPop = 100;
    pbest_number = PR*nPop;
    pbest_number2 = 0.6*nPop;
%     for i=1:nPop   
%         PopPos(i,:)=rand(1,D).*(Xmax-Xmin)+Xmin;
%         PopFit(i)=BenFunctions(PopPos(i,:),F_index,D);   
%     end
    if length(Xmin) ==1
        Rmin = repmat(Xmin,1,D);
        Rmax = repmat(Xmax,1,D);
    end
    VRmin = repmat(Rmin,nPop,1);
    VRmax = repmat(Rmax,nPop,1);
    PopPos = VRmin + (VRmax-VRmin).*rand(nPop,D);
    PopFit = feval(fhd,PopPos',func);
    PopFit = PopFit';
       BestF=inf;
       BestX=[];

    for i=1:nPop
        if PopFit(i)<=BestF
           BestF=PopFit(i);
           BestX=PopPos(i,:);
        end
    end
    
    bsf_fit_var = 1e+30;
    for i = 1 : nPop
        nfes = nfes + 1;
        if PopFit(i) < bsf_fit_var
            bsf_fit_var = PopFit(i);
        end
        if mod(nfes, optimalChart_interval) == 0
            optimalChart = [optimalChart;bsf_fit_var];
        else
            if nfes == max_nfes
                optimalChart = [optimalChart;bsf_fit_var];
            end
        end
        if nfes > max_nfes; break; end
    end
    
        switch FW
            case 1
                SortedIndex = fitnessDistanceBalanceWithWeightGaussian(PopPos, func);
            case 0 
                SortedIndex = fitnessDistanceBalanceWithWeightCauchy(PopPos, func);   
        end
        SortedIndex = SortedIndex';  

while nfes < max_nfes

     Coef=nfes/max_nfes; 
     
       if rand<0.5
          r1=rand;                         
          Beta=2*exp(r1*((max_nfes-nfes+1)/max_nfes))*(sin(2*pi*r1));    
          if  Coef>rand                                                      
              newPopPos(1,:)=BestX+rand(1,D).*(BestX-PopPos(1,:))+Beta*(BestX-PopPos(1,:)); %Equation (4)
          else
              IndivRand=rand(1,D).*(Xmax-Xmin)+Xmin;                                
              newPopPos(1,:)=IndivRand+rand(1,D).*(IndivRand-PopPos(1,:))+Beta*(IndivRand-PopPos(1,:)); %Equation (7)         
          end              
       else 
            Alpha=2*rand(1,D).*(-log(rand(1,D))).^0.5;           
            newPopPos(1,:)=PopPos(1,:)+rand(1,D).*(BestX-PopPos(1,:))+Alpha.*(BestX-PopPos(1,:)); %Equation (1)
       end
     
    for i=2:nPop
   
        if rand<0.5
           r1=rand;                         
           Beta=2*exp(r1*((max_nfes-nfes+1)/max_nfes))*(sin(2*pi*r1));    
             if  Coef>rand                                                      
                 newPopPos(i,:)=BestX+rand(1,D).*(PopPos(i-1,:)-PopPos(i,:))+Beta*(BestX-PopPos(i,:)); %Equation (4)change
             else
                 IndivRand=rand(1,D).*(Xmax-Xmin)+Xmin;                                
                 newPopPos(i,:)=IndivRand+rand(1,D).*(PopPos(i-1,:)-PopPos(i,:))+Beta*(IndivRand-PopPos(i,:));  %Equation (7)       
             end              
        else
            Alpha=2*rand(1,D).*(-log(rand(1,D))).^0.5;           
            newPopPos(i,:)=PopPos(i,:)+rand(1,D).*(PopPos(i-1,:)-PopPos(i,:))+Alpha.*(BestX-PopPos(i,:)); %Equation (1)change
       end         
    end
         
%            for i=1:nPop        
%                newPopPos(i,:)=SpaceBound(newPopPos(i,:),Xmax,Xmin);
%                newPopFit(i)=BenFunctions(newPopPos(i,:),F_index,D);    
%               if newPopFit(i)<PopFit(i)
%                  PopFit(i)=newPopFit(i);
%                  PopPos(i,:)=newPopPos(i,:);
%               end
%            end 
           newPopPos = BoundaryDetection1(newPopPos,Xmin,Xmax);
           newPopFit = feval(fhd,newPopPos',func);
           newPopFit = newPopFit';
           for i=1:nPop
               nfes = nfes + 1;
               if PopFit(i) < bsf_fit_var
                   bsf_fit_var = PopFit(i);
%                    bsf_solution = pos(i, :);
               end
               if newPopFit(i)<PopFit(i)
                   PopFit(i)=newPopFit(i);
                   PopPos(i,:)=newPopPos(i,:);
               end
               if mod(nfes, optimalChart_interval) == 0
                   optimalChart = [optimalChart;bsf_fit_var];
               else
                   if nfes == max_nfes
                       optimalChart = [optimalChart;bsf_fit_var];
                   end
               end
               if nfes > max_nfes; break; end
           end
     
            S=2;
        for i=1:60  
            for ii = 1:999
                pop_choose = SortedIndex(randi(pbest_number));
                if pop_choose ~= i
                    break;
                end
            end
            newPopPos(i,:)=PopPos(i,:)+S*(rand*PopPos(pop_choose,:)-rand*PopPos(i,:)); %Equation (8)
        end
        for i=61:90
            for ii = 1:999
                pop_choose = SortedIndex(randi(pbest_number2));
                if pop_choose ~= i
                    break;
                end
            end
            newPopPos(i,:)=PopPos(i,:)+S*(rand*PopPos(pop_choose,:)-rand*PopPos(i,:));
        end
        for i=91:100
            newPopPos(i,:)=PopPos(i,:)+S*(rand*BestX-rand*PopPos(i,:));
        end
%      for i=1:nPop        
%          newPopPos(i,:)=SpaceBound(newPopPos(i,:),Xmax,Xmin);
%          newPopFit(i)=BenFunctions(newPopPos(i,:),F_index,D);    
%          if newPopFit(i)<PopFit(i)
%             PopFit(i)=newPopFit(i);
%             PopPos(i,:)=newPopPos(i,:);
%          end
%      end
     
     newPopPos = BoundaryDetection1(newPopPos,Xmin,Xmax);
     newPopFit = feval(fhd,newPopPos',func);
     newPopFit = newPopFit';
     for i=1:nPop
         nfes = nfes + 1;
         if PopFit(i) < bsf_fit_var
             bsf_fit_var = PopFit(i);
             %                    bsf_solution = pos(i, :);
         end
         if newPopFit(i)<PopFit(i)
             PopFit(i)=newPopFit(i);
             PopPos(i,:)=newPopPos(i,:);
         end
         if mod(nfes, optimalChart_interval) == 0
             optimalChart = [optimalChart;bsf_fit_var];
         else
             if nfes == max_nfes
                 optimalChart = [optimalChart;bsf_fit_var];
             end
         end
         if nfes > max_nfes; break; end
     end
     
     for i=1:nPop
        if PopFit(i)<BestF
           BestF=PopFit(i);
           BestX=PopPos(i,:);            
        end
     end

       switch FW
            case 1
                SortedIndex = fitnessDistanceBalanceWithWeightGaussian(PopPos, func);
            case 0 
                SortedIndex = fitnessDistanceBalanceWithWeightCauchy(PopPos, func);   
        end
        SortedIndex = SortedIndex';  
  
%      HisBestFit(nfes)=BestF;
end
    optimalChart = optimalChart - optimum;
    optimalChart(optimalChart <= val_2_reach) = 0;
end

