can you fix the spread_metric

% 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 = 200;     % Number of iterations
dim = 30;    % Number of decision variables
lb = zeros(1, dim);  % Lower bound
ub = ones(1, dim);   % Upper bound

num_runs = 10;          % Number of runs for statistical evaluation

% Storage for metrics
GD_values = zeros(num_runs, 1);
IGD_values = zeros(num_runs, 1);
Spacing_values = zeros(num_runs, 1);
Spread_values = zeros(num_runs, 1);
DU_values = zeros(num_runs, 1);
MS_values = zeros(num_runs, 1);  % Added MS storage

% store best pareto
paretoFronts = cell(num_runs, 1);
paretoSets = cell(num_runs, 1);
%% Run MOEMA Multiple Times
for run = 1:num_runs
    fprintf('Run #%d/%d\n', run, num_runs);
    
    % Run MOEMA
    [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = ...
        MOEMA(NPop, T, lb, ub, dim, @zdt1);
    
    
    % Store Pareto front and set
    paretoFronts{run} = paretoFront;
    paretoSets{run} = paretoSet;
    
    % Get True Pareto Front
    true_pareto = get_true_pareto();
    
    % Update the metrics calculation in the main loop:
    [GD_values(run), IGD_values(run), Spacing_values(run), Spread_values(run), ...
        DU_values(run), MS_values(run)] = calculate_metrics(paretoFront, true_pareto);

end

metrics = {GD_values, IGD_values, Spacing_values, Spread_values, DU_values, MS_values};
metric_names = {'GD', 'IGD', 'Spacing', 'Spread', 'DU', 'MS'};

% Find the best run based on the chosen metric (e.g., IGD)
[~, best_run_idx] = min(IGD_values);
best_paretoFront = paretoFronts{best_run_idx};
best_paretoSet = paretoSets{best_run_idx};


% Print results
fprintf('\nPerformance Metrics over %d runs:\n', num_runs);
for i = 1:length(metrics)
    fprintf('\n%s Statistics:\n', metric_names{i});
    fprintf('  Best    = %.4f\n', min(metrics{i}));
    fprintf('  Average = %.4f\n', mean(metrics{i}));
    fprintf('  Median  = %.4f\n', median(metrics{i}));
    fprintf('  Worst   = %.4f\n', max(metrics{i}));
    fprintf('  Std Dev = %.4f\n', std(metrics{i}));
end

%% Plot Results

% Plot the best Pareto Front
figure;
plot(best_paretoFront(:,1), best_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;

% Display metrics for the best run
fprintf('\nBest Run #%d Metrics:\n', best_run_idx);
fprintf('GD    = %.4f\n', GD_values(best_run_idx));
fprintf('IGD   = %.4f\n', IGD_values(best_run_idx));
fprintf('MS    = %.4f\n', MS_values(best_run_idx));
% 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;

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
% Multi-Objective Evolutionary Mating Algorithm (MOEMA) with improvements

% 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:NPop/2, :);  % Eqn (1) in the paper
pop_mum = pop(NPop/2+1:end, :); % Eqn (2) in the paper

% Evaluation of objective function for initialization
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);

for gen = 1:T  % Start iteration
    k1 = randperm(NPop/2);
    k2 = randperm(NPop/2);
    k1x = pop_dad(k1, :);
    k2x = pop_mum(k2, :);
    
    for k = 1:NPop/2
        % Mating process
        Imates = 1 + (var(k1x(k, :)) - var(k2x(k, :))); % Eqn (4) in the paper
        p = randn(1, dim);  % Eqn (6) in the paper
        r = levy_flights(1, dim);  % Using Levy flights
        
        if Imates >= 0  % Upper part of Eqn (5) in the paper
            temp1 = p .* k1x(k, :);
            temp2 = (1-p) .* k2x(k, :);
            pop_child = temp1 + temp2;
        else  % Lower part of Eqn (5) in the paper
            temp1 = p .* k2x(k, :);
            temp2 = (1-p) .* k1x(k, :);
            pop_child = temp1 + temp2;
        end
        
        % Enhanced Impact of best environment (using random non-dominated solution)
        if ~isempty(archive.pop)
            best_pop = archive.pop(randi(size(archive.pop, 1)), :);
            K = rand(1, dim) < 0.8;  % Eqn (8) in the paper
            pop_child = best_pop .* (1-K) + pop_child .* K;  % Eqn (7) in the paper
        end
        
                % Mutation operator
            if rand() < 0.1
                % Special handling for x1
                if rand() < 0.5
                    pop_child(1) = pop_child(1) + randn() * 0.1;
                    pop_child(1) = max(0, min(1, pop_child(1)));
                end
                % Regular mutation for other variables
                mutation_idx = rand(1, dim-1) < 0.1;
                pop_child(2:end) = pop_child(2:end) + mutation_idx .* randn(1, dim-1) * 0.1;
            end
        
        % Bounding of the solution
        pop_child = max(min(ub, pop_child), lb);
        
        % Evaluate child
        obj_child = FITNESSFCN(pop_child);
        
        % Update population using non-dominated sorting
        [pop_dad, obj_dad] = update_population(pop_dad, obj_dad, pop_child, obj_child, k1(k));
        [pop_mum, obj_mum] = update_population(pop_mum, obj_mum, pop_child, obj_child, k2(k));
        
        % Exploration process with levy flight
        if rand() < 0.2
            if ~isempty(archive.pop)
                temp1 = archive.pop(randi(size(archive.pop, 1)), :);
                pop_child = levy_flights(1, dim) .* temp1;  % Exploring with Levy flights
                pop_child = max(min(ub, pop_child), lb);
                obj_child = FITNESSFCN(pop_child);
                
                % Update population using non-dominated sorting
                [pop_dad, obj_dad] = update_population(pop_dad, obj_dad, pop_child, obj_child, randi(NPop/2));
                [pop_mum, obj_mum] = update_population(pop_mum, obj_mum, pop_child, obj_child, randi(NPop/2));
            end
        end
    end
    
    % Combine updated pool dad and mum
    pop = [pop_dad; pop_mum];
    obj = [obj_dad; obj_mum];
    
    % Update archive
    archive = update_archive(archive.pop, archive.obj, pop, obj);
    
    % Record the history at each iteration
    fitness_history(:, gen, :) = obj;
    position_history(:, gen, :) = pop;
    
    % Update Trajectories (using a random non-dominated solution)
    if ~isempty(archive.pop)
        Trajectories(:, gen) = archive.pop(randi(size(archive.pop, 1)), :)';
    end
    
        % Display progress every 10 iterations
    if mod(gen, 10) == 0
        disp(['Iteration #', num2str(gen), ' completed. Archive size: ', num2str(size(archive.pop, 1))]);
    end
end

% Return Pareto front and Pareto set
paretoFront = archive.obj;
paretoSet = archive.pop;

end

% Helper Functions

function pop = initialization(NPop, dim, ub, lb)
    pop = rand(NPop, dim) .* (ub - lb) + lb;
end

function obj = evaluate_population(pop, FITNESSFCN)
    obj = zeros(size(pop, 1), 2);  % Assuming 2 objectives
    for i = 1:size(pop, 1)
        obj(i, :) = FITNESSFCN(pop(i, :));
    end
end

function [pop, obj] = update_population(pop, obj, child, child_obj, index)
    if dominates(child_obj, obj(index, :))
        pop(index, :) = child;
        obj(index, :) = child_obj;
    elseif ~dominates(obj(index, :), child_obj)
        if rand() < 0.5
            pop(index, :) = child;
            obj(index, :) = child_obj;
        end
    end
end

function archive = update_archive(archive_pop, archive_obj, pop, obj)
    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, :);
    archive_size_max = 300;
    % Enhanced crowding distance (using NSGA-II inspired technique)
    if size(archive.pop, 1) > archive_size_max
        crowding_distances = calculate_crowding_distances(archive.obj);
        [~, sorted_indices] = sort(crowding_distances, 'descend');
        archive.pop = archive.pop(sorted_indices(1:archive_size_max), :);
        archive.obj = archive.obj(sorted_indices(1:archive_size_max), :);
    end
end

function distances = calculate_crowding_distances(objectives)
    num_solutions = size(objectives, 1);
    distances = zeros(num_solutions, 1);

    for m = 1:size(objectives, 2)
        [values, idx] = sort(objectives(:, m));
        distances(idx(1)) = inf;
        distances(idx(end)) = inf;

        for i = 2:num_solutions-1
            distances(idx(i)) = distances(idx(i)) + ...
                (values(i+1) - values(i-1)) / (max(values) - min(values));
        end
    end
end

function [is_dominated, domination_count] = pareto_dominance(objectives)
    [n_solutions, n_objectives] = size(objectives);
    is_dominated = false(n_solutions, 1);
    domination_count = zeros(n_solutions, 1);

    for i = 1:n_solutions
        for j = i+1:n_solutions
            if dominates(objectives(i, :), objectives(j, :))
                is_dominated(j) = true;
                domination_count(i) = domination_count(i) + 1;
            elseif dominates(objectives(j, :), objectives(i, :))
                is_dominated(i) = true;
                domination_count(j) = domination_count(j) + 1;
            end
        end
    end
end

function result = dominates(obj1, obj2)
    result = all(obj1 <= obj2) && any(obj1 < obj2);
end

function steps = levy_flights(n, dim, beta)
    % Generates n Levy flight steps in dim dimensions based on Mantegna's algorithm
    if nargin < 3
        beta = 1.5;
    end
    sigma_u = (gamma(1 + beta) * sin(pi * beta / 2) / ...
        (gamma((1 + beta) / 2) * beta * 2^((beta - 1) / 2)))^(1/beta);
    u = random('Normal', 0, sigma_u, n, dim);
    v = random('Normal', 0, 1, n, dim);
    steps = 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

%% Performance Metrics Functions

function [GD, IGD, Spacing, Spread, DU, MS] = calculate_metrics(obtained_pf, true_pf)
    GD = generational_distance(obtained_pf, true_pf);
    IGD = inverted_generational_distance(obtained_pf, true_pf);
    Spacing = spacing_metric(obtained_pf);
    Spread = spread_metric(obtained_pf, true_pf);
    DU = distribution_uniformity(obtained_pf);
    MS = maximum_spread(obtained_pf);
end

function gd = generational_distance(obtained_pf, true_pf)
    if isempty(obtained_pf)
        gd = Inf;
        return;
    end
    distances = zeros(size(obtained_pf, 1), 1);
    for i = 1:size(obtained_pf, 1)
        dist_to_true = sqrt(sum((true_pf - obtained_pf(i,:)).^2, 2));
        distances(i) = min(dist_to_true);
    end
    p = 2;
    gd = (sum(distances.^p) / size(obtained_pf, 1))^(1/p);
end

function igd = inverted_generational_distance(obtained_pf, true_pf)
    if isempty(obtained_pf)
        igd = Inf;
        return;
    end
    distances = zeros(size(true_pf, 1), 1);
    for i = 1:size(true_pf, 1)
        dist_to_obtained = sqrt(sum((obtained_pf - true_pf(i,:)).^2, 2));
        distances(i) = min(dist_to_obtained);
    end
    p = 2;
    igd = (sum(distances.^p) / size(true_pf, 1))^(1/p);
end

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

function spread = spread_metric(obtained_pf, true_pf)
    if size(obtained_pf, 1) < 2
        spread = Inf;
        return;
    end
    
    % Identify extreme points in the true Pareto front
    f1_min = min(true_pf(:,1));
    f2_min = min(true_pf(:,2));
    extreme_points = [f1_min, max(true_pf(:,2)); max(true_pf(:,1)), f2_min];
    
    % Calculate distances to the extreme points
    df = zeros(1, 2);
    for i = 1:2
        dist_to_obtained = sqrt(sum((obtained_pf - extreme_points(i,:)).^2, 2));
        df(i) = min(dist_to_obtained);
    end
    
    % Calculate distances between consecutive points in the obtained Pareto front
    sorted_pf = sortrows(obtained_pf);
    distances = sqrt(sum(diff(sorted_pf).^2, 2));
    d_bar = mean(distances);
    
    spread = (sum(df) + sum(abs(distances - d_bar))) / (sum(df) + length(distances) * d_bar);
end

function du = distribution_uniformity(obtained_pf)
    if size(obtained_pf, 1) < 3
        du = 0;
        return;
    end
    
    % Sort based on the first objective for sequential distance calculation
    sorted_pf = sortrows(obtained_pf);
    
    % Calculate pairwise distances between consecutive solutions
    distances = sqrt(sum(diff(sorted_pf).^2, 2));
    
    % Calculate mean and standard deviation
    mean_dist = mean(distances);
    std_dist = std(distances);
    
    % Original implementation:
    % cv = std_dist / mean_dist;
    % du = 1 / (1 + cv);
    
    % Corrected implementation:
    if mean_dist == 0
        du = 0;
    else
        % Calculate coefficient of variation
        cv = std_dist / mean_dist;
        % Normalize to [0,1] where higher values indicate better uniformity
        du = exp(-cv);  % Alternative: du = 1 / (1 + cv)
    end
end

% New Maximum Spread metric function
function MS = maximum_spread(obtained_pf)
    if isempty(obtained_pf)
        MS = 0;
        return;
    end
    
    % Number of objectives
    o = size(obtained_pf, 2);
    
    % Initialize sum
    sum_spread = 0;
    
    % For each objective
    for i = 1:o
        % Find maximum and minimum values for current objective
        ai = max(obtained_pf(:,i));
        bi = min(obtained_pf(:,i));
        
        % Calculate Euclidean distance between max and min
        spread_i = (ai - bi)^2;
        
        % Add to sum
        sum_spread = sum_spread + spread_i;
    end
    
    % Calculate final Maximum Spread
    MS = sqrt(sum_spread);
end