can you change this ngsa ii code to be like my moema code. show full new ngsa ii code.
ngsa code:
clc 
clear 
close all


nsga2_zdt1;


function nsga2_zdt1()
    %% NSGA-II Parameters
    pop_size = 100;        % Population size
    max_gen = 200;         % Maximum number of generations
    num_vars = 30;         % Number of decision variables
    crossover_prob = 0.9;  % Crossover probability
    mutation_prob = 1/num_vars; % Mutation probability

    %% Initialize Population
    pop = rand(pop_size, num_vars);
    fitness = evaluate_population(pop);

    for gen = 1:max_gen
        %% Non-dominated Sorting
        [fronts, ranks] = non_dominated_sort(fitness);
        
        %% Create Mating Pool
        mating_pool = tournament_selection(pop, ranks, pop_size);
        
        %% Generate Offspring
        offspring = binary_tournament_crossover(mating_pool, crossover_prob, num_vars);
        offspring = polynomial_mutation(offspring, mutation_prob, num_vars);
        
        %% Evaluate Offspring
        offspring_fitness = evaluate_population(offspring);
        
        %% Combine Population and Offspring
        comb_pop = [pop; offspring];
        comb_fitness = [fitness; offspring_fitness];
        
        %% Non-dominated Sorting on Combined Population
        [fronts, ~] = non_dominated_sort(comb_fitness);
        
        %% Create New Population
        pop = [];
        fitness = [];
        for i = 1:length(fronts)
            if size(pop, 1) + size(fronts{i}, 1) > pop_size
                break;
            end
            pop = [pop; comb_pop(fronts{i}, :)];
            fitness = [fitness; comb_fitness(fronts{i}, :)];
        end
        
        %% Crowd Distances and Final Selection
        remaining = pop_size - size(pop, 1);
        if remaining > 0
            selected_front = fronts{i};
            distances = crowding_distance(comb_fitness(selected_front, :));
            [~, sorted_indices] = sort(distances, 'descend');
            selected = selected_front(sorted_indices(1:remaining));
            pop = [pop; comb_pop(selected, :)];
            fitness = [fitness; comb_fitness(selected, :)];
        end
        
        %% Plot Pareto Front
        plot_pareto_front(fitness, gen);
        pause(0.1);
    end
end

function fitness = evaluate_population(pop)
    % Evaluate fitness using ZDT1 function
    fitness = zeros(size(pop, 1), 2);
    f1 = pop(:, 1);
    g = 1 + 9 * mean(pop(:, 2:end), 2);
    h = 1 - sqrt(f1 ./ g);
    f2 = g .* h;
    fitness(:, 1) = f1;
    fitness(:, 2) = f2;
end

function [fronts, rank] = non_dominated_sort(fitness)
    % Perform non-dominated sorting
    num_individuals = size(fitness, 1);
    domination_count = zeros(num_individuals, 1);
    dominated_set = cell(num_individuals, 1);
    rank = zeros(num_individuals, 1);

    fronts = cell(1);
    for i = 1:num_individuals
        for j = 1:num_individuals
            if i ~= j
                if dominates(fitness(i, :), fitness(j, :))
                    dominated_set{i} = [dominated_set{i}, j];
                elseif dominates(fitness(j, :), fitness(i, :))
                    domination_count(i) = domination_count(i) + 1;
                end
            end
        end
        
        if domination_count(i) == 0
            rank(i) = 1;
            fronts{1} = [fronts{1}, i];
        end
    end
    
    current_front = 1;
    while ~isempty(fronts{current_front})
        next_front = [];
        for i = fronts{current_front}
            for j = dominated_set{i}
                domination_count(j) = domination_count(j) - 1;
                if domination_count(j) == 0
                    rank(j) = current_front + 1;
                    next_front = [next_front, j];
                end
            end
        end
        current_front = current_front + 1;
        fronts{current_front} = next_front;
    end
    fronts = fronts(1:current_front-1);
end

function is_dominated = dominates(x, y)
    % Check if solution x dominates solution y
    is_dominated = all(x <= y) && any(x < y);
end

function selected_pop = tournament_selection(pop, ranks, num_selected)
    % Perform binary tournament selection
    selected_pop = zeros(num_selected, size(pop, 2));
    for i = 1:num_selected
        candidates = randperm(size(pop, 1), 2);
        if ranks(candidates(1)) < ranks(candidates(2))
            selected_pop(i, :) = pop(candidates(1), :);
        else
            selected_pop(i, :) = pop(candidates(2), :);
        end
    end
end

function offspring = binary_tournament_crossover(pop, crossover_prob, num_vars)
    % Perform crossover using a binary tournament selection
    offspring = zeros(size(pop));
    for i = 1:2:size(pop, 1)
        if rand < crossover_prob
            parent1 = pop(i, :);
            parent2 = pop(i+1, :);
            cross_point = randi([1, num_vars-1], 1, 1);
            offspring(i, :) = [parent1(1:cross_point), parent2(cross_point+1:end)];
            offspring(i+1, :) = [parent2(1:cross_point), parent1(cross_point+1:end)];
        else
            offspring(i, :) = pop(i, :);
            offspring(i+1, :) = pop(i+1, :);
        end
    end
end

function mutated_pop = polynomial_mutation(pop, mutation_prob, num_vars)
    % Perform polynomial mutation
    eta_m = 20;  % Mutation distribution index
    mutated_pop = pop;
    for i = 1:size(pop, 1)
        for j = 1:num_vars
            if rand < mutation_prob
                r = rand;
                if r < 0.5
                    delta = (2 * r)^(1/(eta_m + 1)) - 1;
                else
                    delta = 1 - (2 * (1 - r))^(1/(eta_m + 1));
                end
                mutated_pop(i, j) = pop(i, j) + delta;
                mutated_pop(i, j) = min(max(mutated_pop(i, j), 0), 1);
            end
        end
    end
end

function distances = crowding_distance(fitness)
    % Compute crowding distance
    num_individuals = size(fitness, 1);
    num_objectives = size(fitness, 2);
    distances = zeros(num_individuals, 1);
    
    for m = 1:num_objectives
        [sorted_values, sorted_idx] = sort(fitness(:, m));
        distances(sorted_idx(1)) = inf;
        distances(sorted_idx(end)) = inf;
        for i = 2:num_individuals-1
            if sorted_values(end) ~= sorted_values(1)
                distances(sorted_idx(i)) = distances(sorted_idx(i)) + ...
                    (sorted_values(i+1) - sorted_values(i-1)) / (sorted_values(end) - sorted_values(1));
            end
        end
    end
end

function plot_pareto_front(fitness, generation)
    % Plot the current Pareto front
    clf;
    plot(fitness(:, 1), fitness(:, 2), 'ro');
    xlim([0, 1]);
    ylim([0, 1]);
    title(['Generation ', num2str(generation)]);
    xlabel('f1');
    ylabel('f2');
    drawnow;
end

moema code:
% MOEMA Implementation and Comparison with True ZDT1 Pareto Front
clc; 
clear; 
close all;

% Set problem parameters
NPop = 100;  % 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

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

% Generate true Pareto front for ZDT1
x1 = linspace(0, 1, 100);
true_pareto_front = [x1', 1 - sqrt(x1)'];

% 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 ZDT1');
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, 'Continuous');
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 = 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

% 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