i have this nsga ii code, can you change the calculate_hypervolume to be like the function that i give you? i try to replace it but error says that cant find the pareto_dominance in that code. my nsga ii doent use pareto_dominance. 

calculate_hypervolume function:
function hypervolume = calculate_hypervolume(obtained_pf, reference_point)
    if isempty(obtained_pf)
        hypervolume = 0;
        return;
    end
    
    % Remove dominated points first
    [is_dominated, ~] = pareto_dominance(obtained_pf);
    obtained_pf = obtained_pf(~is_dominated, :);
    
    % Sort by first objective (ascending)
    obtained_pf = sortrows(obtained_pf);
    
    % Calculate hypervolume
    hypervolume = 0;
    prev_x = obtained_pf(1,1);
    prev_y = reference_point(2);
    
    for i = 1:size(obtained_pf, 1)
        curr_x = obtained_pf(i,1);
        curr_y = obtained_pf(i,2);
        
        % Add rectangular area
        width = curr_x - prev_x;
        height = reference_point(2) - curr_y;
        hypervolume = hypervolume + (width * height);
        
        prev_x = curr_x;
    end
    
    % Add final rectangle
    width = reference_point(1) - prev_x;
    height = reference_point(2) - obtained_pf(end,2);
    hypervolume = hypervolume + (width * height);
end



my nsga code:
clc;
clear all;


nvars = 30;
dim = 30;              % Number of decision variables
num_runs = 10; % Number of independent runs
% Initialize metrics storage
% Update the metrics storage section at the beginning
GD_values = zeros(num_runs, 1);
IGD_values = zeros(num_runs, 1);
Spacing_values = zeros(num_runs, 1);
Spread_values = zeros(num_runs, 1);
Hypervolume_values = zeros(num_runs, 1);
SP_values = zeros(num_runs, 1);
DU_values = zeros(num_runs, 1);     % Distribution Uniformity

% Reference point for hypervolume calculation
reference_point = [1.1, 1.1];  % Adjust based on your problem

for run = 1:num_runs
% Set the number of variables and the objective functions
nvars = 2; % Number of variables
options = optimoptions('gamultiobj', 'PopulationSize', 300, 'Generations', 500, ...
                       'PlotFcn', {@gaplotpareto, @gaplotspread}, 'Display', 'iter');

% Define the bounds for the variables
dim = 30;              % Number of decision variables
lb = zeros(1, dim);    % Lower bound
ub = ones(1, dim);     % Upper bound
% Run the NSGA-II algorithm
[x, fval] = gamultiobj(@zdt1, nvars, [], [], [], [], lb, ub, options);


% Get True Pareto Front
    true_pareto = get_true_pareto();
    
    % Calculate all metrics
        [GD_values(run), IGD_values(run), Spacing_values(run), Spread_values(run),SP_values(run), DU_values(run)] = ...
        calculate_metrics(fval, true_pareto);
    % Calculate hypervolume
    Hypervolume_values(run) = calculate_hypervolume(fval, reference_point);
end
% Display results
figure;
scatter(fval(:,1), fval(:,2), 'filled');
xlabel('Objective 1');
ylabel('Objective 2');
title('Pareto Front');
grid on;


% Update the metrics array when displaying statistics
metrics = {GD_values, IGD_values, Spacing_values, Spread_values, Hypervolume_values, SP_values, DU_values};
metric_names = {'GD', 'IGD', 'Spacing', 'Spread', 'Hypervolume', 'SP', 'DU'};

% Print results
fprintf('\nPerformance Metrics over %d runs:\n', num_runs);
for i = 1:length(metrics)
    fprintf('\n%s Statistics:\n', metric_names{i});
    fprintf('  Best    = %.4f\n', min(metrics{i}));
    fprintf('  Average = %.4f\n', mean(metrics{i}));
    fprintf('  Median  = %.4f\n', median(metrics{i}));
    fprintf('  Worst   = %.4f\n', max(metrics{i}));
    fprintf('  Std Dev = %.4f\n', std(metrics{i}));
end

%% Plot Results
% Plot Pareto Front
figure;
plot(fval(:,1), fval(:,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 ZDT1 using ngsaii');
legend('Obtained Pareto Front', 'True Pareto Front');
grid on;
hold off;



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

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

% Update the calculate_metrics function
function [GD, IGD, Spacing, Spread, SP, DU] = 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);
    Spread = spread_metric(obtained_pf, true_pf);
    SP = sigma_spacing(obtained_pf);
    DU = distribution_uniformity(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);
    distances = zeros(n, 1);
    
    for i = 1:n
        min_dist = Inf;
        for j = 1:n
            if i ~= j
                dist = sqrt(sum((obtained_pf(i,:) - obtained_pf(j,:)).^2));
                min_dist = min(min_dist, dist);
            end
        end
        distances(i) = min_dist;
    end
    
    d_bar = mean(distances);
    spacing = sqrt(sum((distances - d_bar).^2)/(n-1));
end

function spread = spread_metric(obtained_pf, true_pf)
    if size(obtained_pf, 1) < 2
        spread = Inf;
        return;
    end
    
    [~, idx] = sort(obtained_pf(:,1));
    obtained_pf = obtained_pf(idx,:);
    
    df = zeros(1, 2);
    for i = 1:2
        obtained_extreme = obtained_pf(i,:);
        true_extreme = true_pf(i,:);
        df(i) = sqrt(sum((obtained_extreme - true_extreme).^2));
    end
    
    distances = sqrt(sum(diff(obtained_pf).^2, 2));
    d_bar = mean(distances);
    
    spread = (sum(df) + sum(abs(distances - d_bar))) / (sum(df) + size(obtained_pf,1) * d_bar);
end

function hypervolume = calculate_hypervolume(obtained_pf, reference_point)
    if isempty(obtained_pf)
        hypervolume = 0;
        return;
    end
    
    [obtained_pf, ~] = sortrows(obtained_pf);
    
    hypervolume = 0;
    for i = 1:size(obtained_pf, 1)
        if i == 1
            width = reference_point(1) - obtained_pf(i,1);
        else
            width = obtained_pf(i,1) - obtained_pf(i-1,1);
        end
        height = reference_point(2) - obtained_pf(i,2);
        hypervolume = hypervolume + width * height;
    end
end


% Add the new metric functions
function sp = sigma_spacing(obtained_pf)
    if size(obtained_pf, 1) < 2
        sp = Inf;
        return;
    end
    
    n = size(obtained_pf, 1);
    d = zeros(n, 1);
    
    % Calculate distances between consecutive solutions
    for i = 1:n
        min_dist = Inf;
        for j = 1:n
            if i ~= j
                % Euclidean distance between solutions
                dist = sqrt(sum((obtained_pf(i,:) - obtained_pf(j,:)).^2));
                min_dist = min(min_dist, dist);
            end
        end
        d(i) = min_dist;
    end
    
    % Calculate mean distance
    d_bar = mean(d);
    
    % Calculate SP metric
    sp = sqrt(sum((d - d_bar).^2) / (n - 1));
end



% Add the Distribution Uniformity metric function
function du = distribution_uniformity(obtained_pf)
    if size(obtained_pf, 1) < 3
        du = 0;
        return;
    end
    
    % Sort solutions by first objective
    [~, idx] = sort(obtained_pf(:,1));
    sorted_pf = obtained_pf(idx,:);
    
    n = size(sorted_pf, 1);
    distances = zeros(n-1, 1);
    
    % Calculate Euclidean distances between consecutive solutions
    for i = 1:(n-1)
        distances(i) = sqrt(sum((sorted_pf(i+1,:) - sorted_pf(i,:)).^2));
    end
    
    % Calculate mean distance
    mean_dist = mean(distances);
    
    % Calculate standard deviation of distances
    std_dist = std(distances);
    
    % Calculate coefficient of variation (CV)
    cv = std_dist / mean_dist;
    
    % Calculate Distribution Uniformity
    % DU ranges from 0 (poor uniformity) to 1 (perfect uniformity)
    du = 1 / (1 + cv);
    
    % Alternative formulation (uncomment if you prefer):
    % du = exp(-cv);  % This also gives values between 0 and 1
end