add archive max to 100

% MOEMA optimization with Levy Flights and Enhanced Crowding Distance
% Main script to run MOEMA optimization with improvements
clc;
clear;
close all;

% Set problem parameters
NPop = 150;  % Population size
T = 300;     % Number of iterations
dim = 30;    % Number of decision variables
lb = zeros(1, dim);  % Lower bound
ub = ones(1, dim);   % Upper bound

% Run MOEMA with improvements
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, @zdt1);


% Get True Pareto Front
true_pareto = get_true_pareto();
% Plot the Pareto front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro');
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT1 using Enhanced MOEMA');

% Plot Pareto Front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'r');
hold on;
plot(true_pareto(:,1), true_pareto(:,2), 'b--', 'LineWidth', 1.5);
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT1 (MOEMA)');
legend('Obtained Pareto Front', 'True Pareto Front');
grid on;
hold off;

% Calculate and display metrics
disp(['Number of solutions in Pareto Front: ', num2str(size(paretoFront, 1))]);

% Plot the convergence of the Pareto front
figure;
for i = 1:T
    front = squeeze(fitness_history(:,i,:));
    [~, unique_idx] = unique(front, 'rows');
    front = front(unique_idx,:);
    plot(front(:,1), front(:,2), 'bo', 'MarkerSize', 2);
    xlabel('f1');
    ylabel('f2');
    title(['Pareto Front Evolution at Iteration ', num2str(i)]);
    xlim([0 1]);
    ylim([0 1]);
    drawnow;
    pause(0.05);
end

% Plot the trajectory of a random solution
figure;
plot(1:T, Trajectories(1,:));
xlabel('Iteration');
ylabel('x1 value');
title('Trajectory of x1 for a Random Non-dominated Solution');

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history, spacing] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
    % MOEMA: Multi-Objective Evolutionary Mating Algorithm with Lévy Flights
    
    % Preallocation
    fitness_history = zeros(NPop, T, 2); % Assuming 2 objectives
    position_history = zeros(NPop, T, dim);
    Trajectories = zeros(dim, T);
    
    % Generation of initial population
    pop = initialization(NPop, dim, ub, lb);
    pop_dad = pop(1:floor(NPop/2), :); % Parent subset 1
    pop_mum = pop(floor(NPop/2)+1:end, :); % Parent subset 2
    
    % Evaluation of objective function for initial population
    obj_dad = evaluate_population(pop_dad, FITNESSFCN);
    obj_mum = evaluate_population(pop_mum, FITNESSFCN);
    obj = [obj_dad; obj_mum];
    
    % Initialize archive with non-dominated solutions
    archive = update_archive([], [], pop, obj);
    
    % Initialize additional parameters for Lévy flights
    beta = 1.5;  % Lévy flights parameter
    levy_prob = 0.3;  % Probability of applying Lévy flight
    
    % Evolutionary Loop
    for gen = 1:T
        % Shuffle parents
        k1 = randperm(size(pop_dad,1));
        k2 = randperm(size(pop_mum,1));
        k1x = pop_dad(k1, :);
        k2x = pop_mum(k2, :);
        
        % Generate children
        for k = 1:size(k1x,1)
            % Mating process
            imates = randn(1, dim); % Mating factor
            p = rand(1, dim); % Crossover probability
            r = rand();
            
            if mean(imates) >= 0 % Decide mating based on imates
                child = p .* k1x(k, :) + (1 - p) .* k2x(k, :);
            else
                child = p .* k2x(k, :) + (1 - p) .* k1x(k, :);
            end
            
            % Impact of best environment
            if ~isempty(archive.pop)
                best_pop = archive.pop(randi(size(archive.pop,1)), :);
                K = rand(1, dim) < 0.8; % Modification probability
                child = best_pop .* (1 - K) + child .* K;
            end
            
            % Ensure child is within bounds
            child = max(min(child, ub), lb);
            
            % Apply Lévy flight with a certain probability
            if rand() < levy_prob
                levy_step = levy_flight(dim, beta);
                child = child + levy_step;
                % Ensure child is within bounds after Lévy flight
                child = max(min(child, ub), lb);
            end
            
            % Evaluate child
            obj_child = FITNESSFCN(child);
            
            % Update population using non-dominated sorting
            [pop_dad, obj_dad] = update_population(pop_dad, obj_dad, child, obj_child, k);
            [pop_mum, obj_mum] = update_population(pop_mum, obj_mum, child, obj_child, k);
            
            % Exploration process
            if r < 0.2
                if ~isempty(archive.pop)
                    temp1 = archive.pop(randi(size(archive.pop,1)), :);
                    temp_child = rand(1, dim) .* temp1; % Exploration child
                    temp_child = max(min(temp_child, ub), lb);
                    obj_temp_child = FITNESSFCN(temp_child);
                    
                    replace_idx_dad = randi(size(pop_dad,1));
                    replace_idx_mum = randi(size(pop_mum,1));
                    
                    pop_dad(replace_idx_dad, :) = temp_child;
                    obj_dad(replace_idx_dad, :) = obj_temp_child;
                    
                    pop_mum(replace_idx_mum, :) = temp_child;
                    obj_mum(replace_idx_mum, :) = obj_temp_child;
                end
            end
        end
        
        % Combine updated populations
        pop = [pop_dad; pop_mum];
        obj = [obj_dad; obj_mum];
        
        % Update archive with non-dominated solutions
        archive = update_archive(archive.pop, archive.obj, pop, obj);
        
        % Record fitness and positions
        fitness_history(:, gen, :) = obj;
        position_history(:, gen, :) = pop;
        
        % Update Trajectories with a random non-dominated solution
        if ~isempty(archive.pop)
            Trajectories(:, gen) = archive.pop(randi(size(archive.pop,1)), :)';
        end
        
        % Display iteration status
        disp(['Iteration #', num2str(gen), ' completed. Archive size: ', num2str(size(archive.pop,1))]);
    end
    
    % Return Pareto front and Pareto set
    paretoFront = archive.obj;
    paretoSet = archive.pop;
    
    % Calculate spacing metric
    spacing = spacing_metric(paretoFront);
end

% ---------- Spacing Metric Function ----------
function spacing = spacing_metric(obtained_pf)
    if size(obtained_pf, 1) < 2
        spacing = Inf;
        return;
    end
    
    n = size(obtained_pf, 1);
    d = zeros(n, 1);
    
    % Calculate di for each solution
    for i = 1:n
        % Calculate distances to all other solutions
        min_dist = Inf;
        for j = 1:n
            if i ~= j
                % Calculate |f1(x) - f1(x')| + |f2(x) - f2(x')|
                dist = abs(obtained_pf(i,1) - obtained_pf(j,1)) + ...
                       abs(obtained_pf(i,2) - obtained_pf(j,2));
                min_dist = min(min_dist, dist);
            end
        end
        d(i) = min_dist;
    end
    
    % Calculate mean of distances
    d_bar = mean(d);
    
    % Calculate spacing metric SP
    spacing = sqrt(sum((d_bar - d).^2) / (n - 1));
end

% ---------- Initialization Function ----------
function pop = initialization(NPop, dim, ub, lb)
    % INITIALIZATION Generates an initial population within bounds
    pop = rand(NPop, dim) .* (ub - lb) + lb;
end

% ---------- Evaluate Population Function ----------
function obj = evaluate_population(pop, FITNESSFCN)
    % EVALUATE_POPULATION Evaluates the objective function for each individual
    num_individuals = size(pop,1);
    obj = zeros(num_individuals, 2);  % Assuming 2 objectives
    parfor i = 1:num_individuals
        obj(i, :) = FITNESSFCN(pop(i, :));
    end
end

% ---------- Update Population Function ----------
function [pop, obj] = update_population(pop, obj, child, child_obj, index)
    % UPDATE_POPULATION Updates the population based on dominance
    if dominates(child_obj, obj(index, :))
        pop(index, :) = child;
        obj(index, :) = child_obj;
    elseif dominates(obj(index, :), child_obj)
        % Do nothing, parent dominates child
    else
        % Neither dominates; randomly decide to replace
        if rand() < 0.5
            pop(index, :) = child;
            obj(index, :) = child_obj;
        end
    end
end

% ---------- Update Archive Function ----------
function archive = update_archive(archive_pop, archive_obj, pop, obj)
    % UPDATE_ARCHIVE Updates the archive with non-dominated solutions
    combined_pop = [archive_pop; pop];
    combined_obj = [archive_obj; obj];
    
    [is_dominated, ~] = pareto_dominance(combined_obj);
    non_dominated = ~is_dominated;
    
    archive.pop = combined_pop(non_dominated, :);
    archive.obj = combined_obj(non_dominated, :);
end

% ---------- Pareto Dominance Function ----------
function [is_dominated, domination_count] = pareto_dominance(objectives)
    % PARETO_DOMINANCE Determines domination relationships among solutions
    % 
    % Outputs:
    %   is_dominated - Logical array indicating if a solution is dominated
    %   domination_count - Number of times a solution is dominated
    
    n_solutions = size(objectives,1);
    is_dominated = false(n_solutions,1);
    domination_count = zeros(n_solutions, 1);
    
    for i = 1:n_solutions
        for j = 1:n_solutions
            if i ~= j
                if dominates(objectives(j,:), objectives(i,:))
                    is_dominated(i) = true;
                    domination_count(j) = domination_count(j) + 1;
                end
            end
        end
    end
end

% ---------- Dominates Function ----------
function result = dominates(obj1, obj2)
    % DOMINATES Checks if obj1 dominates obj2
    result = all(obj1 <= obj2) && any(obj1 < obj2);
end

function step = levy_flight(dim, beta)
    % LEVY_FLIGHT Generates a Lévy flight step in a "dim"-dimensional space.
    % "beta" is the distribution parameter, typically between 1.5 and 2.
    
    sigma_u = (gamma(1 + beta) * sin(pi * beta / 2) / ...
              (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2)))^(1 / beta);
    u = randn(1, dim) * sigma_u;
    v = randn(1, dim);
    step = u ./ abs(v).^(1 / beta);
end
function f = zdt1(x)
    f1 = x(1);
    g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
    h = 1 - sqrt(f1 / g);
    f2 = g * h;
    f = [f1, f2];
end


function true_pareto = get_true_pareto()
    num_points = 100;
    f1 = linspace(0, 1, num_points)';
    f2 = 1 - sqrt(f1);
    true_pareto = [f1, f2];
end