CAN YOU TRY ON ZDT2:
% NSGA-II optimization
% Main script to run NSGA-II optimization
clc; 
clear; 
close all;

% Set problem parameters
NPop = 100;  % Population size
T = 500;     % Number of iterations
[M, dim, lb, ub] = objective_description_function();  % Get problem parameters
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);

% 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 NSGA-II
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = NSGA2(NPop, T, lb, ub, dim, M);

    
    % 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)] = calculate_metrics(paretoFront, true_pareto);

end

metrics = {GD_values, IGD_values, Spacing_values};
metric_names = {'GD', 'IGD', 'Spacing'};

% 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};


% Plot the final Pareto front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro');
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT1 using NSGA-II');

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

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = NSGA2(NPop, T, lb, ub, dim, M)
% Non-dominated Sorting Genetic Algorithm II (NSGA-II)
% NSGA-II attempts to solve multi-objective problems of the following form:
%  min F(X) = [f1(X), f2(X), ..., fm(X)]  subject to: lb <= X <= ub
%                  
%  lb           - lower bounds on X
%  ub           - upper bounds on X
%  T            - number of iterations
%  NPop         - size of the population
%  M            - number of objectives
%  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, M);
position_history = zeros(NPop, T, dim);
Trajectories = zeros(dim, T);

% Initialize population
chromosome = initialize_variables(NPop, M, dim, lb, ub);

% Sort the initialized population
chromosome = non_domination_sort_mod(chromosome, M, dim);

% Start the evolution process
for i = 1:T
    % Selection
    pool = round(NPop/2);
    tour = 2;
    parent_chromosome = tournament_selection(chromosome, pool, tour);
    
    % Genetic operators (crossover and mutation)
    mu = 20;
    mum = 20;
    offspring_chromosome = genetic_operator(parent_chromosome, M, dim, mu, mum, lb, ub);
    
    % Merge populations
    [main_pop,~] = size(chromosome);
    [offspring_pop,~] = size(offspring_chromosome);
    intermediate_chromosome(1:main_pop,:) = chromosome;
    intermediate_chromosome(main_pop + 1:main_pop + offspring_pop,1:M+dim) = offspring_chromosome;
    
    % Non-dominated sorting
    intermediate_chromosome = non_domination_sort_mod(intermediate_chromosome, M, dim);
    
    % Selection for next generation
    chromosome = replace_chromosome(intermediate_chromosome, M, dim, NPop);
    
    % Record history
    fitness_history(:,i,:) = chromosome(:,dim+1:dim+M);
    position_history(:,i,:) = chromosome(:,1:dim);
    
    % Update Trajectories (using the best ranked solution)
    [~, best_idx] = min(chromosome(:,dim+M+1));
    Trajectories(:,i) = chromosome(best_idx,1:dim)';
    
    if ~mod(i,50)
        disp(['Iteration #', num2str(i), ' completed']);
    end
end

% Extract Pareto front and set
[~, unique_idx] = unique(chromosome(:,dim+1:dim+M), 'rows');
paretoFront = chromosome(unique_idx,dim+1:dim+M);
paretoSet = chromosome(unique_idx,1:dim);

end

function f = initialize_variables(N, M, V, min_range, max_range)

min = min_range;
max = max_range;

K = M + V;

for i = 1 : N
    for j = 1 : V
        f(i,j) = min(j) + (max(j) - min(j))*rand(1);
    end
    f(i,V + 1: K) = evaluate_objective(f(i,:), M, V);
end
end



function f = non_domination_sort_mod(x, M, V)

%% function f = non_domination_sort_mod(x, M, V)
[N, m] = size(x);
clear m

% Initialize the front number to 1.
front = 1;

% There is nothing to this assignment, used only to manipulate easily in
% MATLAB.
F(front).f = [];
individual = [];

%% Non-Dominated sort. 
for i = 1 : N
    % Number of individuals that dominate this individual
    individual(i).n = 0;
    % Individuals which this individual dominate
    individual(i).p = [];
    for j = 1 : N
        dom_less = 0;
        dom_equal = 0;
        dom_more = 0;
        for k = 1 : M
            if (x(i,V + k) < x(j,V + k))
                dom_less = dom_less + 1;
            elseif (x(i,V + k) == x(j,V + k))
                dom_equal = dom_equal + 1;
            else
                dom_more = dom_more + 1;
            end
        end
        if dom_less == 0 && dom_equal ~= M
            individual(i).n = individual(i).n + 1;
        elseif dom_more == 0 && dom_equal ~= M
            individual(i).p = [individual(i).p j];
        end
    end   
    if individual(i).n == 0
        x(i,M + V + 1) = 1;
        F(front).f = [F(front).f i];
    end
end
% Find the subsequent fronts
while ~isempty(F(front).f)
   Q = [];
   for i = 1 : length(F(front).f)
       if ~isempty(individual(F(front).f(i)).p)
        	for j = 1 : length(individual(F(front).f(i)).p)
            	individual(individual(F(front).f(i)).p(j)).n = ...
                	individual(individual(F(front).f(i)).p(j)).n - 1;
        	   	if individual(individual(F(front).f(i)).p(j)).n == 0
               		x(individual(F(front).f(i)).p(j),M + V + 1) = ...
                        front + 1;
                    Q = [Q individual(F(front).f(i)).p(j)];
                end
            end
       end
   end
   front =  front + 1;
   F(front).f = Q;
end

[temp,index_of_fronts] = sort(x(:,M + V + 1));
for i = 1 : length(index_of_fronts)
    sorted_based_on_front(i,:) = x(index_of_fronts(i),:);
end
current_index = 0;

%% Crowding distance
% Find the crowding distance for each individual in each front
for front = 1 : (length(F) - 1)
%    objective = [];
    distance = 0;
    y = [];
    previous_index = current_index + 1;
    for i = 1 : length(F(front).f)
        y(i,:) = sorted_based_on_front(current_index + i,:);
    end
    current_index = current_index + i;
    % Sort each individual based on the objective
    sorted_based_on_objective = [];
    for i = 1 : M
        [sorted_based_on_objective, index_of_objectives] = ...
            sort(y(:,V + i));
        sorted_based_on_objective = [];
        for j = 1 : length(index_of_objectives)
            sorted_based_on_objective(j,:) = y(index_of_objectives(j),:);
        end
        f_max = ...
            sorted_based_on_objective(length(index_of_objectives), V + i);
        f_min = sorted_based_on_objective(1, V + i);
        y(index_of_objectives(length(index_of_objectives)),M + V + 1 + i)...
            = Inf;
        y(index_of_objectives(1),M + V + 1 + i) = Inf;
         for j = 2 : length(index_of_objectives) - 1
            next_obj  = sorted_based_on_objective(j + 1,V + i);
            previous_obj  = sorted_based_on_objective(j - 1,V + i);
            if (f_max - f_min == 0)
                y(index_of_objectives(j),M + V + 1 + i) = Inf;
            else
                y(index_of_objectives(j),M + V + 1 + i) = ...
                     (next_obj - previous_obj)/(f_max - f_min);
            end
         end
    end
    distance = [];
    distance(:,1) = zeros(length(F(front).f),1);
    for i = 1 : M
        distance(:,1) = distance(:,1) + y(:,M + V + 1 + i);
    end
    y(:,M + V + 2) = distance;
    y = y(:,1 : M + V + 2);
    z(previous_index:current_index,:) = y;
end
f = z();
end

function [number_of_objectives, number_of_decision_variables, min_range_of_decesion_variable, max_range_of_decesion_variable] = objective_description_function()
    number_of_objectives = 2;
    number_of_decision_variables = 30;
    min_range_of_decesion_variable = zeros(1, number_of_decision_variables);
    max_range_of_decesion_variable = ones(1, number_of_decision_variables);
end

function f  = replace_chromosome(intermediate_chromosome, M, V,pop)

%% function f  = replace_chromosome(intermediate_chromosome,pro,pop)

[N, m] = size(intermediate_chromosome);

% Get the index for the population sort based on the rank
[temp,index] = sort(intermediate_chromosome(:,M + V + 1));

clear temp m

% Now sort the individuals based on the index
for i = 1 : N
    sorted_chromosome(i,:) = intermediate_chromosome(index(i),:);
end

% Find the maximum rank in the current population
max_rank = max(intermediate_chromosome(:,M + V + 1));

% Start adding each front based on rank and crowing distance until the
% whole population is filled.
previous_index = 0;
for i = 1 : max_rank
    % Get the index for current rank i.e the last the last element in the
    % sorted_chromosome with rank i. 
    current_index = max(find(sorted_chromosome(:,M + V + 1) == i));
    % Check to see if the population is filled if all the individuals with
    % rank i is added to the population. 
    if current_index > pop
        % If so then find the number of individuals with in with current
        % rank i.
        remaining = pop - previous_index;
        % Get information about the individuals in the current rank i.
        temp_pop = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        % Sort the individuals with rank i in the descending order based on
        % the crowding distance.
        [temp_sort,temp_sort_index] = ...
            sort(temp_pop(:, M + V + 2),'descend');
        % Start filling individuals into the population in descending order
        % until the population is filled.
        for j = 1 : remaining
            f(previous_index + j,:) = temp_pop(temp_sort_index(j),:);
        end
        return;
    elseif current_index < pop
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
    else
        % Add all the individuals with rank i into the population.
        f(previous_index + 1 : current_index, :) = ...
            sorted_chromosome(previous_index + 1 : current_index, :);
        return;
    end
    % Get the index for the last added individual.
    previous_index = current_index;
end
end

function f = tournament_selection(chromosome, pool_size, tour_size)

%% function tournament_selection(chromosome, pool_size, tour_size) 

%% Tournament selection process

% Get the size of chromosome. The number of chromosome is not important
% while the number of elements in chromosome are important.
[pop, variables] = size(chromosome);
% The peunltimate element contains the information about rank.
rank = variables - 1;
% The last element contains information about crowding distance.
distance = variables;

% Until the mating pool is filled, perform tournament selection
for i = 1 : pool_size
    % Select n individuals at random, where n = tour_size
    for j = 1 : tour_size
        % Select an individual at random
        candidate(j) = round(pop*rand(1));
        % Make sure that the array starts from one. 
        if candidate(j) == 0
            candidate(j) = 1;
        end
        if j > 1
            % Make sure that same candidate is not choosen.
            while ~isempty(find(candidate(1 : j - 1) == candidate(j)))
                candidate(j) = round(pop*rand(1));
                if candidate(j) == 0
                    candidate(j) = 1;
                end
            end
        end
    end
    % Collect information about the selected candidates.
    for j = 1 : tour_size
        c_obj_rank(j) = chromosome(candidate(j),rank);
        c_obj_distance(j) = chromosome(candidate(j),distance);
    end
    % Find the candidate with the least rank
    min_candidate = ...
        find(c_obj_rank == min(c_obj_rank));
    % If more than one candiate have the least rank then find the candidate
    % within that group having the maximum crowding distance.
    if length(min_candidate) ~= 1
        max_candidate = ...
        find(c_obj_distance(min_candidate) == max(c_obj_distance(min_candidate)));
        % If a few individuals have the least rank and have maximum crowding
        % distance, select only one individual (not at random). 
        if length(max_candidate) ~= 1
            max_candidate = max_candidate(1);
        end
        % Add the selected individual to the mating pool
        f(i,:) = chromosome(candidate(min_candidate(max_candidate)),:);
    else
        % Add the selected individual to the mating pool
        f(i,:) = chromosome(candidate(min_candidate(1)),:);
    end
end
end

function f  = genetic_operator(parent_chromosome, M, V, mu, mum, l_limit, u_limit)

[N,m] = size(parent_chromosome);

clear m
p = 1;
% Flags used to set if crossover and mutation were actually performed. 
was_crossover = 0;
was_mutation = 0;


for i = 1 : N
    % With 90 % probability perform crossover
    if rand(1) < 0.9
        % Initialize the children to be null vector.
        child_1 = [];
        child_2 = [];
        % Select the first parent
        parent_1 = round(N*rand(1));
        if parent_1 < 1
            parent_1 = 1;
        end
        % Select the second parent
        parent_2 = round(N*rand(1));
        if parent_2 < 1
            parent_2 = 1;
        end
        % Make sure both the parents are not the same. 
        while isequal(parent_chromosome(parent_1,:),parent_chromosome(parent_2,:))
            parent_2 = round(N*rand(1));
            if parent_2 < 1
                parent_2 = 1;
            end
        end
        % Get the chromosome information for each randomnly selected
        % parents
        parent_1 = parent_chromosome(parent_1,:);
        parent_2 = parent_chromosome(parent_2,:);
        % Perform corssover for each decision variable in the chromosome.
        for j = 1 : V
            % SBX (Simulated Binary Crossover).
            % For more information about SBX refer the enclosed pdf file.
            % Generate a random number
            u(j) = rand(1);
            if u(j) <= 0.5
                bq(j) = (2*u(j))^(1/(mu+1));
            else
                bq(j) = (1/(2*(1 - u(j))))^(1/(mu+1));
            end
            % Generate the jth element of first child
            child_1(j) = ...
                0.5*(((1 + bq(j))*parent_1(j)) + (1 - bq(j))*parent_2(j));
            % Generate the jth element of second child
            child_2(j) = ...
                0.5*(((1 - bq(j))*parent_1(j)) + (1 + bq(j))*parent_2(j));
            % Make sure that the generated element is within the specified
            % decision space else set it to the appropriate extrema.
            if child_1(j) > u_limit(j)
                child_1(j) = u_limit(j);
            elseif child_1(j) < l_limit(j)
                child_1(j) = l_limit(j);
            end
            if child_2(j) > u_limit(j)
                child_2(j) = u_limit(j);
            elseif child_2(j) < l_limit(j)
                child_2(j) = l_limit(j);
            end
        end
        % Evaluate the objective function for the offsprings and as before
        % concatenate the offspring chromosome with objective value.
        child_1(:,V + 1: M + V) = evaluate_objective(child_1, M, V);
        child_2(:,V + 1: M + V) = evaluate_objective(child_2, M, V);
        % Set the crossover flag. When crossover is performed two children
        % are generate, while when mutation is performed only only child is
        % generated.
        was_crossover = 1;
        was_mutation = 0;
    % With 10 % probability perform mutation. Mutation is based on
    % polynomial mutation. 
    else
        % Select at random the parent.
        parent_3 = round(N*rand(1));
        if parent_3 < 1
            parent_3 = 1;
        end
        % Get the chromosome information for the randomnly selected parent.
        child_3 = parent_chromosome(parent_3,:);
        % Perform mutation on eact element of the selected parent.
        for j = 1 : V
           r(j) = rand(1);
           if r(j) < 0.5
               delta(j) = (2*r(j))^(1/(mum+1)) - 1;
           else
               delta(j) = 1 - (2*(1 - r(j)))^(1/(mum+1));
           end
           % Generate the corresponding child element.
           child_3(j) = child_3(j) + delta(j);
           % Make sure that the generated element is within the decision
           % space.
           if child_3(j) > u_limit(j)
               child_3(j) = u_limit(j);
           elseif child_3(j) < l_limit(j)
               child_3(j) = l_limit(j);
           end
        end
        % Evaluate the objective function for the offspring and as before
        % concatenate the offspring chromosome with objective value.    
        child_3(:,V + 1: M + V) = evaluate_objective(child_3, M, V);
        % Set the mutation flag
        was_mutation = 1;
        was_crossover = 0;
    end
    % Keep proper count and appropriately fill the child variable with all
    % the generated children for the particular generation.
    if was_crossover
        child(p,:) = child_1;
        child(p+1,:) = child_2;
        was_cossover = 0;
        p = p + 2;
    elseif was_mutation
        child(p,:) = child_3(1,1 : M + V);
        was_mutation = 0;
        p = p + 1;
    end
end
f = child;
end


function f = evaluate_objective(x, M, V)
    f = [];
    f(1) = x(1);
    sum = 0;
    for i = 2:V
        sum = sum + x(i);
    end
    g = 1 + 9 * sum / (V - 1);
    h = 1 - sqrt(f(1) / g);
    f(2) = g * h;
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] = 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);
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