ADD TRUE PARETO

clc;
clear;
close all;

% Problem parameters
dim = 30; 
lb = zeros(1, dim);
ub = ones(1, dim);

% Define the fitness function
fitnessFunction = @zdt6;

% Run MOEMA
NPop_MOEMA = 50;
T_MOEMA = 200;
[paretoFront_MOEMA, ~] = MOEMA(NPop_MOEMA, T_MOEMA, lb, ub, dim, fitnessFunction);

% Run NSGA-II
NPop_NSGA2 = 50;
T_NSGA2 = 200;
[paretoFront_NSGA2, ~] = NSGA2(NPop_NSGA2, T_NSGA2, lb, ub, dim, fitnessFunction);

% True Pareto front for ZDT6


% Reference point for hypervolume calculation
reference_point = [1.1, 1.1];

% Calculate metrics for MOEMA
gd_MOEMA = calculate_gd(paretoFront_MOEMA, true_pareto_front);
spread_MOEMA = calculate_spread(paretoFront_MOEMA);
hv_MOEMA = calculate_hypervolume(paretoFront_MOEMA, reference_point);
convergence_MOEMA = calculate_convergence(paretoFront_MOEMA, true_pareto_front);

% Calculate metrics for NSGA-II
gd_NSGA2 = calculate_gd(paretoFront_NSGA2, true_pareto_front);
spread_NSGA2 = calculate_spread(paretoFront_NSGA2);
hv_NSGA2 = calculate_hypervolume(paretoFront_NSGA2, reference_point);
convergence_NSGA2 = calculate_convergence(paretoFront_NSGA2, true_pareto_front);

% Display Comparison
disp('Comparison Table:');
disp('-----------------------------------------------------------');
disp('| Metric       | MOEMA          | NSGA-II       |');
disp('-----------------------------------------------------------');
fprintf('| GD           | %-14.4f | %-14.4f |\n', gd_MOEMA, gd_NSGA2);
fprintf('| Spread       | %-14.4f | %-14.4f |\n', spread_MOEMA, spread_NSGA2);
fprintf('| Hypervolume  | %-14.4f | %-14.4f |\n', hv_MOEMA, hv_NSGA2);
fprintf('| Convergence  | %-14.4f | %-14.4f |\n', convergence_MOEMA, convergence_NSGA2);
disp('-----------------------------------------------------------');

% Plotting Results
figure;
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(paretoFront_NSGA2(:,1), paretoFront_NSGA2(:,2), 'bo', 'DisplayName', 'NSGA-II');
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'k-', 'DisplayName', 'Estimated True Pareto Front');
xlabel('f1');
ylabel('f2');
title('Comparison of MOEMA and NSGA-II on ZDT6');
legend('show');
grid on;

% Plotting Results
figure;

% MOEMA vs True Pareto Front
subplot(1, 2, 1); % 1 row, 2 columns, first plot
plot(paretoFront_MOEMA(:,1), paretoFront_MOEMA(:,2), 'ro', 'DisplayName', 'MOEMA');
hold on;
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'k-', 'DisplayName', 'Estimated True Pareto Front');
xlabel('f1');
ylabel('f2');
title('MOEMA vs Estimated True Pareto Front for ZDT6');
legend('show');
grid on;

% NSGA-II vs True Pareto Front
subplot(1, 2, 2); % 1 row, 2 columns, second plot
plot(paretoFront_NSGA2(:,1), paretoFront_NSGA2(:,2), 'bo', 'DisplayName', 'NSGA-II');
hold on;
plot(true_pareto_front(:,1), true_pareto_front(:,2), 'k-', 'DisplayName', 'Estimated True Pareto Front');
xlabel('f1');
ylabel('f2');
title('NSGA-II vs Estimated True Pareto Front for ZDT6');
legend('show');
grid on;

% Function Definitions

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
    fitness_history = zeros(NPop, T, 2);
    position_history = zeros(NPop, T, dim);
    Trajectories = zeros(dim, T);
    pop = initialization(NPop, dim, ub, lb);
    pop_dad = pop(1:NPop/2, :);
    pop_mum = pop(NPop/2+1:end, :);
    obj_dad = evaluate_population(pop_dad, FITNESSFCN);
    obj_mum = evaluate_population(pop_mum, FITNESSFCN);
    obj = [obj_dad; obj_mum];
    archive = update_archive([], [], pop, obj);

    for gen = 1:T
        k1 = randperm(NPop/2);
        k2 = randperm(NPop/2);
        k1x = pop_dad(k1,:);
        k2x = pop_mum(k2,:);
        
        for k = 1:NPop/2
            Imates = 1 + (var(k1x(k,:)) - var(k2x(k,:)));
            p = randn(1, dim);
            r = rand();
            
            if Imates >= 0
                temp1 = p .* k1x(k,:);
                temp2 = (1-p) .* k2x(k,:);
                pop_child = temp1 + temp2;
            else
                temp1 = p .* k2x(k,:);
                temp2 = (1-p) .* k1x(k,:);
                pop_child = temp1 + temp2;
            end
            
            if ~isempty(archive.pop)
                best_pop = archive.pop(randi(size(archive.pop, 1)), :);
                K = rand(1, dim) < 0.8;
                pop_child = best_pop .* (1-K) + pop_child .* K;
            end
            
            pop_child = max(min(ub, pop_child), lb);
            obj_child = FITNESSFCN(pop_child);
            [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));
            
            if r < 0.2
                if ~isempty(archive.pop)
                    temp1 = archive.pop(randi(size(archive.pop, 1)), :);
                    pop_child = rand(1, dim) .* temp1;
                    pop_child = max(min(ub, pop_child), lb);
                    obj_child = FITNESSFCN(pop_child);
                    [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
        
        pop = [pop_dad; pop_mum];
        obj = [obj_dad; obj_mum];
        archive = update_archive(archive.pop, archive.obj, pop, obj);
        fitness_history(:, gen, :) = obj;
        position_history(:, gen, :) = pop;
        
        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

    paretoFront = archive.obj;
    paretoSet = archive.pop;
end

function [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim, FITNESSFCN)
    pop = initialization(NPop, dim, ub, lb);
    fitness = evaluate_population(pop, FITNESSFCN);

    for gen = 1:T
        [fronts, ranks] = non_dominated_sort(fitness);
        mating_pool = tournament_selection(pop, ranks, NPop);
        offspring = crossover(mating_pool, 0.9, dim);
        offspring = mutation(offspring, 1/dim, dim);
        offspring_fitness = evaluate_population(offspring, FITNESSFCN);
        comb_pop = [pop; offspring];
        comb_fitness = [fitness; offspring_fitness];
        [fronts, ~] = non_dominated_sort(comb_fitness);
        [pop, fitness] = select_new_population(comb_pop, comb_fitness, fronts, NPop);
        disp(['Generation #', num2str(gen), ' completed.']);
    end

    paretoFront = fitness;
    paretoSet = pop;
end

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

function fitness = evaluate_population(pop, FITNESSFCN)
    fitness = zeros(size(pop, 1), 2);
    for i = 1:size(pop, 1)
        fitness(i, :) = FITNESSFCN(pop(i, :));
    end
end

function [fronts, rank] = non_dominated_sort(fitness)
    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)
    is_dominated = all(x <= y) && any(x < y);
end

function selected_pop = tournament_selection(pop, ranks, num_selected)
    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 = crossover(pop, crossover_prob, num_vars)
    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 = mutation(pop, mutation_prob, num_vars)
    eta_m = 20;
    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 [pop, fitness] = select_new_population(comb_pop, comb_fitness, fronts, NPop)
    pop = [];
    fitness = [];
    for i = 1:length(fronts)
        if size(pop, 1) + length(fronts{i}) > NPop
            break;
        end
        pop = [pop; comb_pop(fronts{i}, :)];
        fitness = [fitness; comb_fitness(fronts{i}, :)];
    end
    
    remaining = NPop - 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
end

function distances = crowding_distance(fitness)
    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 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 [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 f = zdt6(x)
    f1 = 1 - exp(-4 * x(1)) * (sin(6 * pi * x(1)))^6;
    g = 1 + 9 * (sum(x(2:end)) / (length(x) - 1))^0.25;
    h = 1 - (f1 / g)^2;
    f2 = g * h;
    f = [f1, f2];
end

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

function gd = calculate_gd(approx_front, true_front)
    sum_of_distances = 0;
    for i = 1:size(approx_front, 1)
        distances = sqrt(sum((true_front - approx_front(i, :)).^2, 2));
        min_distance = min(distances);
        sum_of_distances = sum_of_distances + min_distance^2;
    end
    gd = sqrt(sum_of_distances) / size(approx_front, 1);
end

function hv = calculate_hypervolume(front, reference_point)
    sorted_front = sortrows(front);
    hypervolume = 0;
    for i = 1:size(sorted_front, 1)
        if i == 1
            width = reference_point(1) - sorted_front(i, 1);
        else
            width = sorted_front(i-1, 1) - sorted_front(i, 1);
        end
        height = reference_point(2) - sorted_front(i, 2);
        hypervolume = hypervolume + width * height;
    end
    hv = hypervolume;
end

function convergence = calculate_convergence(approx_front, true_front)
    distances = zeros(size(approx_front, 1), 1);
    for i = 1:size(approx_front, 1)
        dist = sqrt(sum((true_front - approx_front(i, :)).^2, 2));
        distances(i) = min(dist);
    end
    convergence = mean(distances);
end