4.1. Concepts in MOO#

In this notebook we shall see how to

  1. Compute non-dominated set of solutions from a set of pool of points

  2. How to compute hyper volume given a non-dominated set of solutions and a reference point

!pip install matplotlib
Requirement already satisfied: matplotlib in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (3.9.1)
Requirement already satisfied: contourpy>=1.0.1 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (1.2.1)
Requirement already satisfied: cycler>=0.10 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (4.53.1)
Requirement already satisfied: kiwisolver>=1.3.1 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (1.4.5)
Requirement already satisfied: numpy>=1.23 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (2.0.0)
Requirement already satisfied: packaging>=20.0 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (24.1)
Requirement already satisfied: pillow>=8 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (10.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (3.1.2)
Requirement already satisfied: python-dateutil>=2.7 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from matplotlib) (2.9.0.post0)
Requirement already satisfied: six>=1.5 in /home/ksuresh/miniconda3/envs/env_jupyter-book/lib/python3.12/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)

4.1.1. 2D Case#

Let us create some 2D points to start with. We will generate 100 points. We will generate a random points in a circle with radius \(r=(0, 1)\) and \(\theta=(0., 2\pi)\). We will simply use (-x, -y) so that we can find the minimum

# Generate some points. 

import random, math
import matplotlib.pyplot as plt

n_points = 100
r_points = [-1.*random.random() for _ in range(n_points)] #0 - 1
th_points = [random.uniform(0., 2*math.pi) for _ in range(n_points)] #circle

points = [(r_points[i]*abs(math.cos(th_points[i])), r_points[i]*abs(math.sin(th_points[i]))) for i in range(n_points)]


# Print the results.
plt.figure(1, figsize=(8, 6))
plt.plot([p[0] for p in points], [p[1] for p in points], 'bo', label='All points')
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.legend()
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.show()
../_images/f96da24259f915be1e316aa0e3329c7306a87b944765b383399098694a6ce72a.png

4.1.2. Non dominated sorting#

Non-dominated sorting is a technique used in multi-objective optimization to identify the Pareto-optimal solutions from a set of candidate solutions. In multi-objective optimization, we aim to optimize multiple conflicting objectives simultaneously.

Mathematically, let’s consider a set of solutions S, where each solution s ∈ S is represented by a vector of objective values f(s) = (f1(s), f2(s), …, fn(s)). The goal of non-dominated sorting is to partition the set S into different fronts, where each front contains solutions that are not dominated by any other solution in the same front.

A solution s1 dominates another solution s2 if and only if:

  1. For all objectives i, \(f_{i}(s_{1}) ≤ f_{i}(s_{2})\), and

  2. There exists at least one objective j such that fj(s1) < fj(s2).

To perform non-dominated sorting, we can use the following steps:

  1. Initialize an empty list of fronts, F.

  2. Initialize an empty set of solutions that are dominated by other solutions, D.

  3. For each solution s in S, do the following: a. Initialize the domination count, n_dom, as 0. b. Initialize an empty set of solutions that s dominates, S_dom. c. For each solution p in S, excluding s, do the following: i. If s dominates p, add p to S_dom. ii. If p dominates s, increment n_dom by 1. d. If n_dom is 0, add s to the first front (F[0]). e. If n_dom > 0, add s to D and store S_dom as the set of solutions dominated by s.

  4. Initialize the front index, i, as 0.

  5. While F[i] is not empty, do the following: a. Initialize an empty list for the next front, F_next. b. For each solution s in F[i], do the following: i. For each solution q in the set of solutions dominated by s, do the following: - Decrement the domination count of q by 1. - If the domination count of q becomes 0, add q to F_next. c. Increment i by 1. d. Append F_next to F.

  6. Return the list of fronts, F, where each front contains non-dominated solutions.

Now, let’s build a function that performs non-dominated sorting:

In the above code, the dominates() function is a helper function that checks if one solution dominates another based on the defined dominance criteria.

You can use this non_dominated_sorting() function to perform non-dominated sorting on a set of solutions and obtain the Pareto-optimal fronts.

def dominates(point1, point2, epsilon=1e-4):
    """ Check if point1 dominates point2 based on two criteria (minimization). """
    return all(p1 <= p2 + epsilon for p1, p2 in zip(point1, point2)) and any(p1 < p2 - epsilon for p1, p2 in zip(point1, point2))

def non_dominated_solutions(points, epsilon=1e-4):
    """ Find non-dominated solutions (Pareto front) from a list of points. """
    dominated = set()
    pareto_front = []

    for i, point1 in enumerate(points):
        if i in dominated:
            continue
        pareto_front.append(point1)
        for j, point2 in enumerate(points[i + 1:], start=i + 1):
            if dominates(point1, point2, epsilon):
                dominated.add(j)
            elif dominates(point2, point1, epsilon):
                pareto_front.pop()
                break

    return pareto_front

4.1.2.1. Computational complexity of this#

Outer Loop (First For Loop): The outer loop iterates over each point in the list points. If there are n points, this loop runs n times.

  • Complexity: \(O(n)\)

Inner Loop (Second For Loop):

The inner loop iterates over each subsequent point in the list points starting from i + 1 (where i is the current index of the outer loop).

For each point point1, it compares against all subsequent points point2 to determine dominance.

In the worst case, if all points are non-dominated, each point point1 might be compared against all n−i−1 points in the inner loop.

  • Complexity: \(O(n^{2})\) in the worst case, where each point is compared against all other points.

Domination Check (dominates Function):

The dominates function compares two points in three dimensions, which is O(3)O(3) or O(1)O(1) constant time operation (since the dimensions are fixed).

Space Complexity:

The space complexity is primarily determined by the dominated set and the pareto_front list. In the worst case, both can grow up to nn in size, resulting in \(O(n)\) space complexity.

Total Computational Complexity

Given the analysis above, the overall computational complexity of finding the Pareto front using the provided algorithm is dominated by the inner loop, which is \(O(n^{2})\).

Let us now plot and check if it works

pareto_front = non_dominated_solutions(points)
print("Non-dominated solutions (Pareto front):", pareto_front)

plt.figure(2, figsize=(8, 6))
plt.plot([p[0] for p in points], [p[1] for p in points], 'bo', label='All points')
plt.plot([p[0] for p in pareto_front], [p[1] for p in pareto_front], 'ro', label='Nondominated points')
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.legend()
plt.show()
Non-dominated solutions (Pareto front): [(-0.09676379681866248, -0.9553806067701416), (-0.5117974398510517, -0.7031760273411101), (-0.8158186848182191, -0.5165419084313226), (-0.5904123883075125, -0.5865480895269028), (-0.573098607513865, -0.6925197023102118), (-0.8355484588095629, -0.16109942723013831), (-0.4085392301659479, -0.882716216855052), (-0.3567279409792068, -0.909257444054572), (-0.9860305671983675, -0.10139295706907155)]
../_images/705b2fa45d79d14678c87c58ad923fa6abf2720bccd585bef330909009417d09.png

4.1.3. Excersise:#

4.1.3.1. 1. Can you implement the same function but for maximization?#

All we have to do is reverse the condition for checking in the dominates function which we will rename as dominates_max

### 1. Can you implement the same function but for maximization?


def dominates_max(point1, point2, epsilon=1e-4):
    """ Check if point1 dominates point2 based on two criteria (maximization). """
    return all(p1 >= p2 - epsilon for p1, p2 in zip(point1, point2)) and any(p1 > p2 + epsilon for p1, p2 in zip(point1, point2))

def non_dominated_solutions_max(points, epsilon=1e-4):
    """ Find non-dominated solutions (Pareto front) from a list of points. """
    dominated = set()
    pareto_front = []

    for i, point1 in enumerate(points):
        if i in dominated:
            continue
        pareto_front.append(point1)
        for j, point2 in enumerate(points[i + 1:], start=i + 1):
            if dominates_max(point1, point2, epsilon):
                dominated.add(j)
            elif dominates_max(point2, point1, epsilon):
                pareto_front.pop()
                break

    return pareto_front

pareto_front_max = non_dominated_solutions_max(points)
print("Non-dominated solutions (Pareto front):", pareto_front_max)

plt.figure(3, figsize=(8, 6))
plt.plot([p[0] for p in points], [p[1] for p in points], 'bo', label='All points')
plt.plot([p[0] for p in pareto_front_max], [p[1] for p in pareto_front_max], 'ro', label='Nondominated points')
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.legend()
plt.show()
Non-dominated solutions (Pareto front): [(-1.105930824760222e-05, -1.8573612026052515e-05, -0.0013778235202420222), (-0.00034035104544824333, -0.0002536689241251127, -5.548986617222313e-06)]
../_images/6cbb50a88bbfa1ea947dbd8d393c6bac41884ecbc715ddb84252513fb293deb6.png

4.1.3.2. 2. Can you implement the same function but for 3D points?#

### 2. Can you implement the same function but for 3D points?

n_points = 100
r_points = [-1.*random.random() for _ in range(n_points)] #0 - 1
th_points = [random.uniform(0., 2*math.pi) for _ in range(n_points)] #circle
phi_points = [random.uniform(0., 2*math.pi) for _ in range(n_points)] #sphere

points = [(r_points[i]*abs(math.cos(th_points[i]))*abs(math.sin(phi_points[i])), r_points[i]*abs(math.sin(th_points[i]))*abs(math.sin(phi_points[i])), r_points[i]*abs(math.cos(phi_points[i])) ) for i in range(n_points)]

# Print the results.
fig = plt.figure(4, figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter([p[0] for p in points], [p[1] for p in points], [p[2] for p in points], c='b', marker='o', label='All points')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.legend()
plt.show()
../_images/85f3473644348eb15c2b3c62d7ae8054d3610b520249de72de5dcd5a21841371.png
def dominates_3d(point1, point2, epsilon=1e-4):
    """ Check if point1 dominates point2 based on three criteria (minimization). """
    return all(p1 <= p2 + epsilon for p1, p2 in zip(point1, point2)) and any(p1 < p2 - epsilon for p1, p2 in zip(point1, point2))

def non_dominated_solutions_3d(points, epsilon=1e-4):
    """ Find non-dominated solutions (Pareto front) from a list of points. """
    dominated = set()
    pareto_front = []

    for i, point1 in enumerate(points):
        if i in dominated:
            continue
        pareto_front.append(point1)
        for j, point2 in enumerate(points[i + 1:], start=i + 1):
            if dominates_3d(point1, point2, epsilon):
                dominated.add(j)
            elif dominates_3d(point2, point1, epsilon):
                pareto_front.pop()
                break

    return pareto_front

pareto_front_3d = non_dominated_solutions_3d(points)
print("Non-dominated solutions (Pareto front):", pareto_front_3d)

fig = plt.figure(5, figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
ax.scatter([p[0] for p in points], [p[1] for p in points], [p[2] for p in points], c='b', marker='o', label='All points')
ax.scatter([p[0] for p in pareto_front_3d], [p[1] for p in pareto_front_3d], [p[2] for p in pareto_front_3d], c='r', marker='o', label='Nondominated points')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.legend()
plt.show()
Non-dominated solutions (Pareto front): [(-0.4711494960914555, -0.03280412382989566, -0.8087017144838647), (-0.10521106213075945, -0.05243209811668072, -0.9308992474820867), (-0.11341875697361008, -0.1714778796784776, -0.9217319117640743), (-0.15110270310817106, -0.8175149324635373, -0.4435783346636459), (-0.179552691921969, -0.6810295303255397, -0.36626627138967766), (-0.023369416669509747, -0.7030823159513367, -0.4803387611307129), (-0.8931082978803174, -0.1954611576723414, -0.24582123036924156), (-0.5464066126106423, -0.5340872760805806, -0.6405963992946936), (-0.4007825522672473, -0.21759152346234623, -0.8268683820135337), (-0.8437482457792418, -0.3298938991021348, -0.2831924788342274), (-0.7978669089329151, -0.3429411991217389, -0.015941090113372303), (-0.2818984498862414, -0.7082816164408218, -0.30637492531136923), (-0.768065610701401, -0.48462788888364583, -0.05685226394412791), (-0.13624168266641054, -0.9584437951517359, -0.044875187943089076), (-0.056903644806693514, -0.32769105892553607, -0.8248631389499478), (-0.41070692700609557, -0.4083400112482499, -0.7646712355933492), (-0.30235097591824145, -0.6136204661219998, -0.488188613837117), (-0.5063735070200795, -0.0007195311379880111, -0.7157694065003615)]
../_images/f2861599ca2112195977ec924f5fd822ab509f0157e195e3fc925220f2b7329b.png

4.1.3.3. 3. Can you have an isight on how to go about a N-D case?#

### 3. Can you have an isight on how to go about a `N-D` case?

# Even the 2D case the way it is written checks this by default. Hence, Domainated sorting is not a big issue

def dominates_nd(point1, point2, epsilon=1e-4):
    """ Check if point1 dominates point2 based on N criteria (minimization). """
    return all(p1 <= p2 + epsilon for p1, p2 in zip(point1, point2)) and any(p1 < p2 - epsilon for p1, p2 in zip(point1, point2))

def non_dominated_solutions_nd(points, epsilon=1e-4):
    """ Find non-dominated solutions (Pareto front) from a list of points. """
    dominated = set()
    pareto_front = []

    for i, point1 in enumerate(points):
        if i in dominated:
            continue
        pareto_front.append(point1)
        for j, point2 in enumerate(points[i + 1:], start=i + 1):
            if dominates_nd(point1, point2, epsilon):
                dominated.add(j)
            elif dominates_nd(point2, point1, epsilon):
                pareto_front.pop()
                break

    return pareto_front

4.1.4. Computing Hyper Volume#

Let us remind ourself about MOO and what hypervolume is

image.png

In multi-objective optimization (MOO), the goal is to find a set of solutions that represents the trade-off between multiple conflicting objectives. One way to evaluate the quality of these solutions is by using a metric called hypervolume.

4.1.4.1. Hypervolume#

Hypervolume is a measure of the dominated space in the objective space. It quantifies the extent of the Pareto-optimal front, which represents the best possible trade-off between the objectives. The hypervolume metric considers both the spread and density of solutions in the objective space.

The hypervolume is calculated by defining a reference point, which represents the ideal values for each objective. The hypervolume is then computed as the volume of the dominated space between the Pareto-optimal front and the reference point.

4.1.4.2. Convergence Tracking#

Convergence tracking is an important aspect of MOO algorithms. It helps us understand how well the algorithm is progressing towards finding the Pareto-optimal front. Hypervolume is a useful tool for tracking convergence because it provides a quantitative measure of the quality of the obtained solutions.

As the MOO algorithm iteratively improves the solutions, the hypervolume should increase. A higher hypervolume indicates a larger dominated space, which means that the algorithm is finding better solutions that are closer to the true Pareto-optimal front.

By monitoring the hypervolume over iterations, we can assess the convergence behavior of the algorithm. If the hypervolume increases rapidly and then plateaus, it may indicate that the algorithm has reached a local optimum. On the other hand, if the hypervolume continues to increase steadily, it suggests that the algorithm is exploring the solution space effectively and converging towards the true Pareto-optimal front.

4.1.4.3. Tracking Convergence with Hypervolume#

To track convergence using hypervolume, we can calculate the hypervolume at each iteration of the MOO algorithm and plot it over time. This allows us to visualize the progress of the algorithm and make informed decisions about when to stop the optimization process.

By analyzing the hypervolume plot, we can gain insights into the convergence behavior of the MOO algorithm. If the hypervolume increases steadily and reaches a plateau, it suggests that the algorithm has converged. On the other hand, if the hypervolume fluctuates or decreases, it indicates that the algorithm is not converging effectively and may require further adjustments.

Tracking convergence using hypervolume provides a quantitative measure of the optimization progress and helps in making informed decisions during the MOO process.

Tip

Can you think on how can you compute this hypervolume for a 2D case?

import numpy as np
import matplotlib.pyplot as plt

def compute_hypervolume_old(points, reference_point):
    """
    Compute the hypervolume from a given set of 2D points which are non-dominating solutions.
    
    Args:
        points (list): List of 2D points representing non-dominating solutions.
        reference_point (tuple): Reference point as a 2D point.
        
    Returns:
        float: The hypervolume of the non-dominating solutions.
    """
    # Sort the points based on the x-coordinate in ascending order
    sorted_points = sorted(points, key=lambda p: p[0])
    
    # Initialize the hypervolume to 0
    hypervolume = 0
    
    # Initialize the previous y-coordinate to the maximum y-coordinate
    prev_y = max(sorted_points, key=lambda p: p[1])[1]
    
    # Iterate over the sorted points
    for point in sorted_points:
        # Calculate the width of the rectangle
        width = point[0] - sorted_points[0][0]
        
        # Calculate the height of the rectangle
        height = point[1] - prev_y
        
        # Update the hypervolume
        hypervolume += width * height
        
        # Update the previous y-coordinate
        prev_y = point[1]
    
    # Calculate the width of the last rectangle
    width = reference_point[0] - sorted_points[0][0]
    
    # Calculate the height of the last rectangle
    height = reference_point[1] - prev_y
    
    # Update the hypervolume with the last rectangle
    hypervolume += width * height
    
    return hypervolume
import numpy as np
import matplotlib.pyplot as plt

def compute_hypervolume(points, reference_point):
    """
    Compute the hypervolume from a given set of 2D points which are non-dominating solutions.
    This method does 
    
    Args:
        points (list): List of 2D points representing non-dominating solutions.
        reference_point (tuple): Reference point as a 2D point.
        
    Returns:
        float: The hypervolume of the non-dominating solutions.
    """
    # Sort the points based on the x-coordinate in ascending order
    sorted_points = sorted(points, key=lambda p: p[0])
    
    # Initialize the hypervolume to 0
    hypervolume = 0
    
    
    # Iterate over the sorted points
    for i in range(len(sorted_points)-1):
        
        # Update the hypervolume
        # x1, y1 = reference_point
        # x2, y2 = point
        # x3, y3 = point[i+1]
        x1, y1 = reference_point
        x2, y2 = sorted_points[i]
        x3, y3 = sorted_points[i+1]
        if x1 < x2 or x1 < x3:
            continue
        else:
            hypervolume += (0.5*(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2)))
    
    return hypervolume
# draw the hypervolume

def draw_hypervolume(points, reference_point):
    """
    Draw the hypervolume from a given set of 2D points which are non-dominating solutions.
    
    Args:
        points (list): List of 2D points representing non-dominating solutions.
        reference_point (tuple): Reference point as a 2D point.
    """
    # Sort the points based on the x-coordinate in ascending order
    sorted_points = sorted(points, key=lambda p: p[0])
    
    # Initialize the hypervolume to 0
    hypervolume = 0
    
    # Initialize the plot
    fig, ax = plt.subplots()
    
    # Iterate over the sorted points
    for i in range(len(sorted_points)-1):
        
        # Update the hypervolume
        x1, y1 = reference_point
        x2, y2 = sorted_points[i]
        x3, y3 = sorted_points[i+1]
        if x1 < x2 or x1 < x3:
            continue
        else:
            hypervolume += (0.5*(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2)))
        
        # Plot the hypervolume
        polygon = plt.Polygon([reference_point, sorted_points[i], sorted_points[i+1]], edgecolor='r', facecolor='r', alpha=0.2)
        ax.add_patch(polygon)
    
    # Plot the points
    ax.scatter([p[0] for p in points], [p[1] for p in points], c='b', marker='o', label='All points')
    
    # Plot the reference point
    ax.scatter(reference_point[0], reference_point[1], c='g', marker='x', label='Reference point')
    
    # Set the labels
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    
    # Set the title
    ax.set_title(f'Hypervolume: {hypervolume}')
    
    # Show the plot
    plt.legend()
    plt.show()
n_points = 100
r_points = [-1.*random.random() for _ in range(n_points)] #0 - 1
th_points = [random.uniform(0., 2*math.pi) for _ in range(n_points)] #circle

points = [(r_points[i]*abs(math.cos(th_points[i])), r_points[i]*abs(math.sin(th_points[i]))) for i in range(n_points)]

pareto_front = non_dominated_solutions(points)

reference_point = (0., 0.)

hypervolume = compute_hypervolume(pareto_front, reference_point)
print("Hypervolume:", hypervolume)
draw_hypervolume(pareto_front, reference_point)
Hypervolume: 0.647175261766771
../_images/ef8c16ba4a4cf7cfd663a09814e768aeb6e293589beff8d8621550036da09af8.png

4.1.5. Excersise#

4.1.5.1. 0. Could you perhaps check if HV is monotonic as you go through different fronts? Should it increase or decrease?#

4.1.5.2. 1. Can you identify the computational complexity of this function?#

4.1.5.3. 2. Can you extend this to 3D and how would you generalize this to n-D objectives?#

4.1.5.4. 3. This way of computing hypervolume can result in having smaller non monotonic functions. How can one circumvent that (Thanks and Credits to Shannon Zelitch)#

4.1.5.5. 0. Could you perhaps check if HV is monotonic as you go through different fronts? Should it increase or decrease?#

### 0. Could you perhaps check if HV is monotonic as you go through different fronts? Should it increase or decrease?

hv_list = []
points_cpy = points.copy()
reference_point = (0.5, 0.5)
while(len(points_cpy) > 1):
    pareto_front = non_dominated_solutions(points_cpy)
    hv = compute_hypervolume(pareto_front, reference_point)
    hv_list.append(hv)
    for p in pareto_front:
        points_cpy.remove(p)
plt.plot(hv_list)
plt.xlabel("Iterations")
plt.xticks(range(len(hv_list)), range(len(hv_list)-1, -1, -1))
plt.ylabel("Hypervolume")
plt.show()
../_images/a7cdc3e397e9780ce8902925e3d280b4823741dbce9f232795116c254caa50f7.png

4.1.5.6. 2. Can you extend this to 3D and how would you generalize this to n-D objectives?#

### 2. 3D

def compute_hypervolume_3d(points, reference_point):
    """
    Compute the hypervolume from a given set of 3D points which are non-dominating solutions.
    
    Args:
        points (list): List of 3D points representing non-dominating solutions.
        reference_point (tuple): Reference point as a 3D point.
        
    Returns:
        float: The hypervolume of the non-dominating solutions.
    """
    # Sort the points based on the x-coordinate in ascending order
    sorted_points = sorted(points, key=lambda p: p[0])
    
    # Initialize the hypervolume to 0
    hypervolume = 0
    
    # Initialize the previous y-coordinate to the maximum y-coordinate
    prev_y = max(sorted_points, key=lambda p: p[1])[1]
    
    # Initialize the previous z-coordinate to the maximum z-coordinate
    prev_z = max(sorted_points, key=lambda p: p[2])[2]
    
    # Iterate over the sorted points
    for point in sorted_points:
        # Calculate the width of the rectangle
        width = point[0] - sorted_points[0][0]
        
        # Calculate the height of the rectangle
        height = point[1] - prev_y
        
        # Calculate the depth of the rectangle
        depth = point[2] - prev_z
        
        # Update the hypervolume
        hypervolume += width * height * depth
        
        # Update the previous y-coordinate
        prev_y = point[1]
        
        # Update the previous z-coordinate
        prev_z = point[2]
    
    # Calculate the width of the last rectangle
    width = reference_point[0] - sorted_points[0][0]
    
    # Calculate the height of the last rectangle
    height = reference_point[1] - prev_y
    
    # Calculate the depth of the last rectangle
    depth = reference_point[2] - prev_z
    
    # Update the hypervolume with the last rectangle
    hypervolume += width * height * depth
    
    return hypervolume
n_points = 1000
r_points = [-1.*random.random() for _ in range(n_points)] #0 - 1
th_points = [random.uniform(0., 2*math.pi) for _ in range(n_points)] #circle
phi_points = [random.uniform(0., 2*math.pi) for _ in range(n_points)] #sphere

points = [(r_points[i]*abs(math.cos(th_points[i]))*abs(math.sin(phi_points[i])), r_points[i]*abs(math.sin(th_points[i]))*abs(math.sin(phi_points[i])), r_points[i]*abs(math.cos(phi_points[i])) ) for i in range(n_points)]

pareto_front = non_dominated_solutions_3d(points)

reference_point = (0., 0., 0.)
# Print the results.
fig = plt.figure(11, figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')

plt.plot([p[0] for p in points], [p[1] for p in points], [p[-1] for p in points], 'bo', label='All points')
plt.plot(pareto_front[0], pareto_front[1], pareto_front[2], 'ro', label='Nondominated points')
plt.xlabel('X', fontsize=14)
plt.ylabel('Y', fontsize=14)
plt.legend()
plt.show()
../_images/0196e7c5015b74d00203cfa160e27034f9b62f278ba1fb106751ad6a49433498.png
## compute hv trend

hv_list = []
points_cpy = points.copy()
reference_point = (10., 10., 10.)
while(len(points_cpy) > 1):
    pareto_front = non_dominated_solutions_3d(points_cpy)
    hv = compute_hypervolume_3d(pareto_front, reference_point)
    hv_list.append(hv)
    for p in pareto_front:
        points_cpy.remove(p)
plt.plot(hv_list)
plt.xlabel("Iterations")
plt.xticks(range(len(hv_list)), range(len(hv_list)-1, -1, -1))
plt.ylabel("Hypervolume")
Text(0, 0.5, 'Hypervolume')
../_images/d6d7587dbec6f62a9b10eda128bfbbc173c04c779fb9a96a88316089e9c3c029.png