CAN YOU DISPLAY THE  best_GD_values, average_GD_values, median_GD_values, worst_GD_values, best_IGD_values, average_IGD_values, median_IGD_values, worst_IGD_values,  Convergence metric and Diversity metric. GIVE ME FULL CODE


clc; 
clear; 
close all;

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

% Run MOEMA for ZDT3
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, @zdt3);

% Generate true Pareto front for ZDT3
x1 = linspace(0, 1, 500);
x2 = 1 - sqrt(x1) - x1.*sin(10*pi*x1);
true_pareto_front = [x1', x2'];

% Plot the results
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'b-', 'DisplayName', 'True Pareto Front');
xlabel('f1');
ylabel('f2');
title('MOEMA vs True Pareto Front for ZDT3');
legend('show');
grid on;

% Calculate metrics
moema_min_f1 = min(paretoFront(:,1));
moema_min_f2 = min(paretoFront(:,2));
true_min_f1 = min(true_pareto_front(:,1));
true_min_f2 = min(true_pareto_front(:,2));

moema_max_f1 = max(paretoFront(:,1));
moema_max_f2 = max(paretoFront(:,2));
true_max_f1 = max(true_pareto_front(:,1));
true_max_f2 = max(true_pareto_front(:,2));

moema_solutions = size(paretoFront, 1);

% Calculate spread
moema_spread = calculate_spread(paretoFront);
true_spread = calculate_spread(true_pareto_front);

% Display comparison table
disp('Comparison Table:');
disp('------------------------------------------------------------');
disp('| Metric          | MOEMA          | True Pareto Front     |');
disp('------------------------------------------------------------');
fprintf('| Min f1           | %-14.4f | %-22.4f |\n', moema_min_f1, true_min_f1);
fprintf('| Min f2           | %-14.4f | %-22.4f |\n', moema_min_f2, true_min_f2);
fprintf('| Max f1           | %-14.4f | %-22.4f |\n', moema_max_f1, true_max_f1);
fprintf('| Max f2           | %-14.4f | %-22.4f |\n', moema_max_f2, true_max_f2);
fprintf('| Number of points | %-14d | %-22s |\n', moema_solutions, 'Discontinuous');
fprintf('| Spread           | %-14.4f | %-22.4f |\n', moema_spread, true_spread);
disp('------------------------------------------------------------');

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
    % Multi-Objective Evolutionary Mating Algorithm (MOEMA)
    % MOEMA attempts to solve multi-objective problems of the following form:
    %  min F(X) = [f1(X), f2(X), ..., fm(X)]  subject to: lb <= X <= ub
    %                  
    %  FITNESSFCN   - function handle of the multi-objective fitness function
    %  lb           - lower bounds on X
    %  ub           - upper bounds on X
    %  T            - number of iterations
    %  NPop         - size of the population (class size)  
    %  paretoFront  - set of non-dominated solutions in objective space
    %  paretoSet    - set of non-dominated solutions in decision space
    %  Trajectories - best solution trajectory over iterations
    %  fitness_history - history of fitness values
    %  position_history - history of population positions

    % 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 = rand();
            
            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
            
            % 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
            
            % 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
            if r < 0.2
                if ~isempty(archive.pop)
                    temp1 = archive.pop(randi(size(archive.pop, 1)), :);
                    pop_child = rand(1, dim) .* temp1;  % Eqn (9) in the paper
                    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
        
        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;
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, :);
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 f = zdt3(x)
    f1 = x(1);
    g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
    h = 1 - sqrt(f1 / g) - (f1 / g) * sin(10 * pi * f1);
    f2 = g * h;
    f = [f1, f2];
end

% New function to calculate spread
function spread = calculate_spread(front)
    sorted_front = sortrows(front);
    distances = sqrt(sum(diff(sorted_front).^2, 2));
    spread = std(distances) / mean(distances);
end