% Initialize the parameters
N = 50; % Number of bees
D = 10; % Number of dimensions
L = 100; % Limit for abandoning a food source
MaxIter = 1000; % Maximum number of iterations

% Load the F2.mat file that contains the data for F2 function
% load('F2.mat'); % This will create variables z, M, and o in the workspace

% Define the F2 function as a function handle
f = @(x) sum(100*(x(2:end)-x(1:end-1).^2).^2+(x(1:end-1)-1).^2);

% Set the lower and upper bounds of the search space
lb = -100; % Lower bound
ub = 100; % Upper bound

% Initialize the population
X = lb + (ub - lb) * rand(N, D); % Randomly generate N bees in D dimensions
Fit = f(M'*(X-z)+o); % Evaluate the fitness of each bee using the F2 function

% Find the best bee
[BestFit, BestIndex] = min(Fit); % Find the minimum fitness and its index
BestX = X(BestIndex, :); % Find the best bee position

% Initialize the counters
iter = 1; % Iteration counter
trial = zeros(N, 1); % Abandonment counter

% Initialize some performance indicators
BestFitHistory = zeros(MaxIter, 1); % Best fitness value in each iteration
MeanFitHistory = zeros(MaxIter, 1); % Mean fitness value in each iteration
StdFitHistory = zeros(MaxIter, 1); % Standard deviation of fitness values in each iteration

% Main loop
while iter <= MaxIter
    
    % Employed bees phase
    for i = 1:N
        % Choose a partner bee randomly
        j = randi(N);
        while j == i
            j = randi(N);
        end
        
        % Generate a new bee position by perturbing the current one
        phi = -1 + 2 * rand(1, D); % Random number in [-1, 1]
        v = X(i, :) + phi .* (X(i, :) - X(j, :)); % New position
        
        % Apply boundary constraints
        v(v < lb) = lb;
        v(v > ub) = ub;
        
        % Evaluate the new bee position using the F2 function
        NewFit = f(M*(v-z)+o); % Apply shift, rotation, and origin transformation
        
        % Apply greedy selection
        if NewFit < Fit(i)
            X(i, :) = v; % Replace the old bee with the new one
            Fit(i) = NewFit; % Update the fitness
            trial(i) = 0; % Reset the abandonment counter
            
            % Apply local search with a small probability 
            p_ls = 0.01; % Local search probability 
            if rand < p_ls 
                u = v + 0.001 * randn(1, D); % Add a small Gaussian noise to the new position 
                u(u < lb) = lb; 
                u(u > ub) = ub; 
                NewFit2 = f(M*(u-z)+o); % Apply shift, rotation, and origin transformation 
                if NewFit2 < NewFit 
                    X(i, :) = u; 
                    Fit(i) = NewFit2; 
                end 
            end 
            
        else
            trial(i) = trial(i) + 1; % Increment the abandonment counter
        end
        
    end
    
    % Calculate the probability of each bee based on its fitness
    P = Fit / sum(Fit); % Probability vector
    
    % Onlooker bees phase
    for i = 1:N
        % Choose a food source based on roulette wheel selection
        r = rand; % Random number in [0, 1]
        k = find(r <= cumsum(P), 1); % Selected food source index
        
        % Choose a partner bee randomly
        j = randi(N);
        while j == k
            j = randi(N);
        end
        
        % Generate a new bee position by perturbing the current one
        phi = -1 + 2 * rand(1, D); % Random number in [-1, 1]
        v = X(k, :) + phi .* (X(k, :) - X(j, :)); % New position
        
        % Apply boundary constraints
        v(v < lb) = lb;
        v(v > ub) = ub;
        
        % Evaluate the new bee position using the F2 function
        NewFit = f(M*(v-z)+o); % Apply shift, rotation, and origin transformation
        
        % Apply greedy selection
        if NewFit < Fit(k)
            X(k, :) = v; % Replace the old bee with the new one
            Fit(k) = NewFit; % Update the fitness
            trial(k) = 0; % Reset the abandonment counter
            
            % Apply local search with a small probability 
            p_ls = 0.01; % Local search probability 
            if rand < p_ls 
                u = v + 0.001 * randn(1, D); % Add a small Gaussian noise to the new position 
                u(u < lb) = lb; 
                u(u > ub) = ub; 
                NewFit2 = f(M*(u-z)+o); % Apply shift, rotation, and origin transformation 
                if NewFit2 < NewFit 
                    X(k, :) = u; 
                    Fit(k) = NewFit2; 
                end 
            end 
            
        else
            trial(k) = trial(k) + 1; % Increment the abandonment counter
        end
        
    end
    
    % Scout bees phase
    for i = 1:N
        if trial(i) > L 
            X(i, :) = lb + (ub - lb) * rand(1, D); % Generate a new random position for the scout bee 
            Fit(i) = f(M*(X(i, :)'-z)+o); % Evaluate its fitness using the F2 function
            trial(i) = 0; % Reset its abandonment counter 
        end 
    end 
    
    % Find the best bee 
    [NewBestFit, NewBestIndex] = min(Fit); % Find the minimum fitness and its index 
    if NewBestFit < BestFit 
       BestFit = NewBestFit; % Update the best fitness 
       BestX = X(NewBestIndex, :); % Update the best position 
    end 
    
    % Update the performance indicators
    BestFitHistory(iter) = BestFit; % Record the best fitness value in this iteration
    MeanFitHistory(iter) = mean(Fit); % Record the mean fitness value in this iteration
    StdFitHistory(iter) = std(Fit); % Record the standard deviation of fitness values in this iteration
    
    % Increment the iteration counter 
    iter = iter + 1; 
    
end 

% Display the results 
disp('Best fitness value:'); 
disp(BestFit); 
disp('Best position:'); 
disp(BestX);

% Plot the performance indicators
figure;
plot(1:MaxIter, BestFitHistory, 'r', 'LineWidth', 2); % Plot the best fitness value vs iteration
hold on;
plot(1:MaxIter, MeanFitHistory, 'b', 'LineWidth', 2); % Plot the mean fitness value vs iteration
hold on;
plot(1:MaxIter, StdFitHistory, 'g', 'LineWidth', 2); % Plot the standard deviation of fitness values vs iteration
hold off;
xlabel('Iteration');
ylabel('Fitness value');
legend('Best', 'Mean', 'Std');
title('Performance of ABC algorithm on F2 function');

% Compare with other ABC variants from 
% Define the other ABC variants as functions
ABC_GA = @(f, D, lb, ub, MaxIter) abc_ga(f, D, lb * ones(1, D), ub * ones(1, D), MaxIter);
ABC_LS = @(f, D, lb, ub, MaxIter) abc_ls(f, D, lb * ones(1, D), ub * ones(1, D), MaxIter);
ABC_DE = @(f, D, lb, ub, MaxIter) abc_de(f, D, lb * ones(1, D), ub * ones(1, D), MaxIter);

% Run the other ABC variants on the same function and parameters
[ABC_GA_BestX, ABC_GA_BestFit] = ABC_GA(f, D, lb, ub, MaxIter);
[ABC_LS_BestX, ABC_LS_BestFit] = ABC_LS(f, D, lb, ub, MaxIter);
[ABC_DE_BestX, ABC_DE_BestFit] = ABC_DE(f, D, lb, ub, MaxIter);

% Display the results of other ABC variants
disp('ABC-GA best fitness value:');
disp(ABC_GA_BestFit);
disp('ABC-GA best position:');
disp(ABC_GA_BestX);
disp('ABC-LS best fitness value:');
disp(ABC_LS_BestFit);
disp('ABC-LS best position:');
disp(ABC_LS_BestX);
disp('ABC-DE best fitness value:');
disp(ABC_DE_BestFit);
disp('ABC-DE best position:');
disp(ABC_DE_BestX);

% Plot the convergence curves of ABC and its variants
figure;
plot(1:MaxIter, BestFitHistory, 'r', 'LineWidth', 2); % Plot the best fitness value of ABC vs iteration
hold on;
plot(1:MaxIter, ABC_GA_BestFit, 'b', 'LineWidth', 2); % Plot the best fitness value of ABC-GA vs iteration
hold on;
plot(1:MaxIter, ABC_LS_BestFit, 'g', 'LineWidth', 2); % Plot the best fitness value of ABC-LS vs iteration
hold on;
plot(1:MaxIter, ABC_DE_BestFit, 'm', 'LineWidth', 2); % Plot the best fitness value of ABC-DE vs iteration
hold off;
xlabel('Iteration');
ylabel('Fitness value');
legend('ABC', 'ABC-GA', 'ABC-LS', 'ABC-DE');
title('Convergence curves of ABC and its variants on F2 function');