% ----------------------------------------------------------------------- % 
% Flower pollenation algorithm (FPA), or flower algorithm                 %
% Programmed by Xin-She Yang @ Feb to May 2012                            % 
% Updated and modified by XS Yang in Sept 2012                            %
% - Adapted by Pedro Bento November 2023 
% ----------------------------------------------------------------------- %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Notes: This demo program contains the very basic components of          %
% the flower pollination algorithm (FPA), or flower algorithm (FA),       % 
% for single objective optimization.    It usually works well for         %
% unconstrained functions only. For functions/problems with               % 
% limits/bounds and constraints, constraint-handling techniques           %
% should be implemented to deal with constrained problems properly.       %
%                                                                         %    
% References -- citation details:                                         % 
% 1) Xin-She Yang, Flower pollination algorithm for global optimization,  %
%    Unconventional Computation and Natural Computation,                  %
%    Lecture Notes in Computer Science, Vol. 7445, pp. 240-249 (2012).    %
% 2) X. S. Yang, M. Karamanoglu, X. S. He, Multi-objective flower         %
%    algorithm for optimization, Procedia in Computer Science,            %
%    vol. 18, pp. 861-868 (2013).                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fmin,best,convergence,time_FPA]=FPA_orig_corr(n,N_iter,lb,ub,dim,fobj,Pos,InitFit)
tic
p=0.8;           % The probabibility of switch
% Dimension of the search variables
Lb=lb*ones(1,dim);     % Lower bounds
Ub=ub*ones(1,dim);      % Upper bounds
% Initialize the population/solutions
Sol=Pos;
Fitness=InitFit;
% Find the current best among the initial population
[fmin,I]=min(Fitness);
best=Sol(I,:);
S=Sol;
convergence=zeros(1,N_iter);
t=1;
convergence(t)=fmin;
% Start the main iterations -- Flower Pollination Algorithm
for t=2:N_iter,
    % Loop over all bats/solutions
    for i=1:n,
        % Pollens are carried by insects and thus can move on large scale
        % This L should be drawn from the Levy distribution
        % Formula: x_i^{t+1}=x_i^t+ L (x_i^t-gbest)
        if rand>p,     % Probability (p) is checked after drawing a rand
            L=Levy(dim);     % Draw random numbers from a Levy distribution
            % The main search mechanism via flower pollination
            dS=L.*(Sol(i,:)-best);      % Caclulate the step increments
            S(i,:)=Sol(i,:)+dS;         % Update the new solutions
            % Check new solutions to satisfy the simple limits/bounds
            S(i,:)=simplebounds(S(i,:),Lb,Ub);
            % Another search move via local pollenation of neighbor flowers
        else
            epsilon=rand;
            % Find random pollen/flowers in the neighbourhood
            JK=randperm(n);
            % As they are random, the first two entries also random
            % If the flower are the same or similar species, then
            % they can be pollenated, otherwise, no action is needed.
            % Formula: x_i^{t+1}+epsilon*(x_j^t-x_k^t)
            S(i,:)=S(i,:)+epsilon*(Sol(JK(1),:)-Sol(JK(2),:));
            % Check if the simple limits/bounds are met
            S(i,:)=simplebounds(S(i,:),Lb,Ub);
        end
        
        % Evaluate the objective values of the new solutions
        Fnew=fobj(S(i,:));
        % If fitness improves (better solutions found), update then
        if (Fnew<=Fitness(i)),
            Sol(i,:)=S(i,:);
            Fitness(i)=Fnew;
        end
        
        % Update the current global best among the population
        if Fnew<=fmin,
            best=S(i,:);
            fmin=Fnew;
        end
    end
    convergence(t)=fmin;
end  % End of the main iterations
time_FPA=toc;
%________________________________________________________________%
%% Auxiliary nested functions
%________________________________________________________________%
% Application of simple lower bounds and upper bounds
    function s=simplebounds(s,Lb,Ub)
        % Apply the lower bound
        ns_tmp=s;
        Is=ns_tmp<Lb;
        ns_tmp(Is)=Lb(Is);
        % Apply the upper bounds
        J=ns_tmp>Ub;
        ns_tmp(J)=Ub(J);
        % Update this new move
        s=ns_tmp;
    end
%________________________________________________________________%
% Draw n samples for the Levy flight from the Levy distribution
    function L=Levy(d)
        % For details of the Levy flights, see Chapter 11 of the following book:
        % Xin-She Yang, Nature-Inspired Optimization Algorithms, Elsevier, (2014).
        beta=3/2;
        sigma=(gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
        % Mantegna's algorithm for Levy random numbers
        u=randn(1,d)*sigma;
        v=randn(1,d);
        step=u./abs(v).^(1/beta);
        L=0.01*step;             % Final Levy steps
    end
%________________________________________________________________%
%________________________________________________________________%
end %function