%% FOX: A Fox-inspired Optimization Algorithm
% Authors: Hardi M. Mohammed, Tarik A. Rashid.
% Implementation by: Hardi M. Mohammed
% 
% Cite as:
% Mohammed, H., Rashid, T. FOX: a FOX-inspired optimization algorithm.
% Appl Intell (2022). https://doi.org/10.1007/s10489-022-03533-0

function [sFeat, Sf, Nf, Convergence_curve] = F(feat, label, SearchAgents_no, Max_iter, HO)

%% Parameter Initialization
lb  = 0;  % Lower bound of the search space
ub  = 1;  % Upper bound of the search space

dim = size(feat, 2); % Number of features (dimension of search space)
Best_pos = zeros(1, dim); % Best position found
Best_score = inf; % Best fitness value (set to -inf for maximization problems)
MinT = inf; % Minimum time variable for motion calculations

% Initialize the positions of search agents
X = initialization(SearchAgents_no, dim, ub, lb);

% Initialize arrays
Distance_Fox_Rat = zeros(SearchAgents_no, dim); % Distance between fox and prey (rat)
Convergence_curve = zeros(1, Max_iter); % Stores best score in each iteration

% Loop counter
l = 0;

% FOX algorithm specific parameters
c1 = 0.18; % Range: [0, 0.18]
c2 = 0.82; % Range: [0.19, 1]

%% Main Optimization Loop
while l < Max_iter
    for i = 1:size(X, 1)
        % Boundary constraint handling
        Flag4ub = X(i, :) > ub;
        Flag4lb = X(i, :) < lb;
        X(i, :) = (X(i, :) .* ~(Flag4ub + Flag4lb)) + ub .* Flag4ub + lb .* Flag4lb;
        
        % Evaluate fitness function
        fun = @jFitnessFunction;
        fitness = fun(feat, label, (X(i, :) > 0.5), HO);
        
        % Update the best solution found so far
        if fitness < Best_score 
            Best_score = fitness; 
            Best_pos = X(i, :);
        end
    end
    
    % Adaptive parameter 'a' decreasing over iterations
    a = 2 * (1 - (l / Max_iter));
    Jump = 0; % Jump length initialization
    
    % Movement phase of FOX algorithm
    for i = 1:size(X, 1)
        r = rand;    % Random number for movement decision
        p = rand;    % Another random number for motion scaling
        
        if r >= 0.5  % Fox movement behavior
            Time(i, :) = rand(1, dim);
            sps = Best_pos ./ Time(i, :); % Speed calculation
            Distance_S_Travel(i, :) = sps .* Time(i, :);
            Distance_Fox_Rat(i, :) = 0.5 * Distance_S_Travel(i, :);
            
            tt = sum(Time(i, :)) / dim;
            t = tt / 2;
            Jump = 0.5 * 9.81 * t^2; % Jump calculation based on physics
            
            if p > 0.18
                X(i, :) = Distance_Fox_Rat(i, :) .* Jump * c1;
            else
                X(i, :) = Distance_Fox_Rat(i, :) .* Jump * c2;
            end
            
            % Update minimum time
            if MinT > tt
                MinT = tt;
            end 
        else
            % Random walk exploration step
            X(i, :) = Best_pos + randn(1, dim) .* (MinT * a);
        end
    end
    
    % Increment iteration counter
    l = l + 1;
    
    % Store best fitness value for convergence analysis
    Convergence_curve(l) = Best_score;
end

%% Feature Selection Output
Pos = 1:dim; % Feature indices
Sf = Pos((Best_pos > 0.5) == 1); % Selected features
Nf = length(Sf); % Number of selected features
sFeat = feat(:, Sf); % Final selected feature subset

end