is all the evaluation metric in here is used correctly?


clc;
clear;
close all;

% Set problem parameters
NPop = 200;  % 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 = 2;  % Number of runs for statistical evaluation

% Storage for metrics
metrics_storage = struct('GD', zeros(num_runs, 1), ...
                        'IGD', zeros(num_runs, 1), ...
                        'Spacing', zeros(num_runs, 1), ...
                        'Spread', zeros(num_runs, 1), ...
                        'HV', zeros(num_runs, 1));

%% Run NSGA-II Multiple Times
for run = 1:num_runs
    fprintf('Run #%d/%d\n', run, num_runs);
    
   % Run NSGA-II
   [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim);
    
    % Get True Pareto Front
    true_pareto = get_true_pareto();
    
     % Calculate metrics
    [metrics_storage.GD(run), ...
     metrics_storage.IGD(run), ...
     metrics_storage.Spacing(run), ...
     metrics_storage.Spread(run), ...
     metrics_storage.HV(run)] = calculate_metrics(paretoFront, true_pareto);
end

% Calculate statistics
metric_names = fieldnames(metrics_storage);
for i = 1:length(metric_names)
    metric_values = metrics_storage.(metric_names{i});
    fprintf('\n%s Statistics:\n', metric_names{i});
    fprintf('  Best    = %.4f\n', min(metric_values));
    fprintf('  Worst   = %.4f\n', max(metric_values));
    fprintf('  Mean    = %.4f\n', mean(metric_values));
    fprintf('  Median  = %.4f\n', median(metric_values));
    fprintf('  Std Dev = %.4f\n', std(metric_values));
end

% Visualize results
figure;
boxplot([metrics_storage.GD, metrics_storage.IGD, ...
         metrics_storage.Spacing, metrics_storage.Spread, ...
         metrics_storage.HV], ...
        'Labels', {'GD', 'IGD', 'Spacing', 'Spread', 'HV'});
title('Performance Metrics Distribution');
ylabel('Metric Value');
%}

%% Plot Results
% 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 zdt2 using NSGA-II');
legend('Obtained Pareto Front', 'True Pareto Front');
grid on;
hold off;

% Plot the results
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro', 'DisplayName', 'NSGA-II');
xlabel('f1');
ylabel('f2');
title('NSGA-II Pareto Front for zdt2');
legend('show');
grid on;

function [paretoFront, paretoSet] = NSGA2(NPop, T, lb, ub, dim)
    % Initialize population
    pop = initialization(NPop, dim, ub, lb);
    fitness = evaluate_population(pop, @zdt2);

    for gen = 1:T
        % Non-dominated sorting
        [fronts, ranks] = non_dominated_sort(fitness);
        
        % Create mating pool using binary tournament selection
        mating_pool = tournament_selection(pop, ranks, NPop);
        
        % Generate offspring using SBX crossover and polynomial mutation
        offspring = simulated_binary_crossover(mating_pool, dim, 20, 0.9);
        offspring = polynomial_mutation(offspring, 1/dim, 20);
        
        % Evaluate offspring
        offspring_fitness = evaluate_population(offspring, @zdt2);
        
        % 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] = select_new_population(comb_pop, comb_fitness, fronts, NPop);
        
        % Display generation progress
        disp(['Generation #', num2str(gen), ' completed.']);
    end

    % Return last population as Pareto front
    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);  % Assuming 2 objectives
    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 = simulated_binary_crossover(pop, num_vars, eta_c, crossover_prob)
    offspring = zeros(size(pop));
    for i = 1:2:size(pop, 1)
        parent1 = pop(i, :);
        parent2 = pop(i+1, :);
        child1 = parent1;
        child2 = parent2;
        if rand < crossover_prob
            for j = 1:num_vars
                if rand <= 0.5
                    if abs(parent1(j) - parent2(j)) > 1e-14 % Avoid division by nearly zero
                        beta = 1.0 + (2*min(parent1(j), parent2(j)) / abs(parent1(j) - parent2(j)));
                        alpha = 2.0 - beta^(-(eta_c+1));
                        rand_num = rand;
                        if rand_num <= 1.0/alpha
                            beta_q = (rand_num*alpha)^(1.0/(eta_c+1));
                        else
                            beta_q = (1.0/(2.0 - rand_num*alpha))^(1.0/(eta_c+1));
                        end
                        child1(j) = 0.5*((1 + beta_q)*parent1(j) + (1 - beta_q)*parent2(j));
                        child2(j) = 0.5*((1 - beta_q)*parent1(j) + (1 + beta_q)*parent2(j));
                    else
                        child1(j) = parent1(j);
                        child2(j) = parent2(j);
                    end
                else
                    child1(j) = parent1(j);
                    child2(j) = parent2(j);
                end
            end
        end
        offspring(i, :) = child1;
        offspring(i+1, :) = child2;
    end
end

function mutated_pop = polynomial_mutation(pop, mutation_prob, eta_m)
    num_vars = size(pop, 2);
    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
            range = sorted_values(end) - sorted_values(1);
            if range > 0
                distances(sorted_idx(i)) = distances(sorted_idx(i)) + ...
                (sorted_values(i+1) - sorted_values(i-1)) / range;
            end
        end
    end
end
function f = zdt2(x)
    f1 = x(1);
    g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
    h = 1 - (f1 / g)^2;
    f2 = g * h;
    f = [f1, f2];
end

% ------------------- True Pareto Front Function -------------------
function true_pareto = get_true_pareto()
    % Generate the true Pareto front for ZDT2
    num_points = 1000;
    f1 = linspace(0, 1, num_points)';
    f2 = 1 - (f1).^2;
    true_pareto = [f1, f2];
end

%% Performance Metrics Functions
function [GD, IGD, Spacing, Spread, HV] = calculate_metrics(obtained_pf, true_pf)
    % Input validation
    if isempty(obtained_pf) || isempty(true_pf)
        GD = Inf; IGD = Inf; Spacing = Inf; Spread = Inf; HV = 0;
        return;
    end
    
    % Normalize both fronts using the same bounds
    [obtained_pf_norm, true_pf_norm] = normalize_both_fronts(obtained_pf, true_pf);
    
    % Calculate metrics
    GD = generational_distance(obtained_pf_norm, true_pf_norm);
    IGD = inverted_generational_distance(obtained_pf_norm, true_pf_norm);
    Spacing = spacing_metric(obtained_pf_norm);
    Spread = improved_spread_metric(obtained_pf_norm, true_pf_norm);
    HV = hypervolume(obtained_pf_norm, [1.1, 1.1]); % Reference point slightly above normalized maximum
end

function [normalized_obtained, normalized_true] = normalize_both_fronts(obtained_pf, true_pf)
    % Combine both fronts to find global min and max
    all_points = [obtained_pf; true_pf];
    min_vals = min(all_points);
    max_vals = max(all_points);
    
    % Avoid division by zero
    range = max_vals - min_vals;
    range(range == 0) = 1;
    
    % Normalize both fronts using the same bounds
    normalized_obtained = (obtained_pf - min_vals) ./ range;
    normalized_true = (true_pf - min_vals) ./ range;
end



function gd = generational_distance(obtained_pf, true_pf)
    % Calculates the average distance from each obtained point to the nearest true point
    distances = zeros(size(obtained_pf, 1), 1);
    
    for i = 1:size(obtained_pf, 1)
        % Calculate Euclidean distances to all true points
        dist_to_true = sqrt(sum((true_pf - obtained_pf(i,:)).^2, 2));
        distances(i) = min(dist_to_true);
    end
    
    % Use p-norm with p=2 (Euclidean norm)
    gd = sqrt(sum(distances.^2) / size(obtained_pf, 1));
end

function igd = inverted_generational_distance(obtained_pf, true_pf)
    % Calculates the average distance from each true point to the nearest obtained point
    distances = zeros(size(true_pf, 1), 1);
    
    for i = 1:size(true_pf, 1)
        % Calculate Euclidean distances to all obtained points
        dist_to_obtained = sqrt(sum((obtained_pf - true_pf(i,:)).^2, 2));
        distances(i) = min(dist_to_obtained);
    end
    
    % Use p-norm with p=2 (Euclidean norm)
    igd = sqrt(sum(distances.^2) / size(true_pf, 1));
end

function spacing = spacing_metric(obtained_pf)
    if size(obtained_pf, 1) < 2
        spacing = Inf;
        return;
    end
    n = size(obtained_pf, 1);
    distances = zeros(n, 1);
    for i = 1:n
        dist_to_others = sqrt(sum((obtained_pf - obtained_pf(i,:)).^2, 2));
        dist_to_others(i) = Inf;  % Exclude self-distance
        distances(i) = min(dist_to_others);
    end
    spacing = sqrt(sum((distances - mean(distances)).^2) / (n - 1));
end

function spread = improved_spread_metric(obtained_pf, true_pf)
    if size(obtained_pf, 1) < 2
        spread = Inf;
        return;
    end
    
    % Find extreme points
    df = zeros(1, size(obtained_pf, 2));
    for i = 1:size(obtained_pf, 2)
        [~, obtained_idx] = max(obtained_pf(:,i));
        [~, true_idx] = max(true_pf(:,i));
        df(i) = norm(obtained_pf(obtained_idx,:) - true_pf(true_idx,:));
    end
    
    % Sort points and calculate consecutive distances
    [~, idx] = sortrows(obtained_pf);
    sorted_pf = obtained_pf(idx,:);
    distances = sqrt(sum((sorted_pf(1:end-1,:) - sorted_pf(2:end,:)).^2, 2));
    
    d_bar = mean(distances);
    spread = (sum(df) + sum(abs(distances - d_bar))) / ...
            (sum(df) + length(distances) * d_bar);
end


function hv = hypervolume_nd(points, reference)
    if size(points, 2) == 2
        return hypervolume_2d(points, reference);
    end
    
    points = sortrows(points, size(points, 2));
    exclusive = zeros(size(points, 1), 1);
    
    for i = 1:size(points, 1)
        slice = points(1:i-1, 1:end-1);
        exclusive(i) = (reference(end) - points(i,end)) * ...
                      hypervolume_nd(slice, reference(1:end-1));
    end
    
    hv = sum(exclusive);
end

function min_vals = minimum_up_to(arr)
    min_vals = zeros(size(arr));
    min_vals(1) = arr(1);
    for i = 2:length(arr)
        min_vals(i) = min(min_vals(i-1), arr(i));
    end
end

function non_dominated = find_non_dominated(points)
    n = size(points, 1);
    is_dominated = false(n, 1);
    
    for i = 1:n
        for j = 1:n
            if i ~= j && all(points(j,:) <= points(i,:)) && any(points(j,:) < points(i,:))
                is_dominated(i) = true;
                break;
            end
        end
    end
    
    non_dominated = ~is_dominated;
end