Merge pull request anyoptimization#80 from Peng-YM/NDSort

Implement the efficient non-dominated sorting

pymoo/usage/usage_non_dominated_sorting.py

@@ -0,0 +1,22 @@
+from timeit import timeit
+# noinspection PyUnresolvedReferences
+from pymoo.util.nds.non_dominated_sorting import NonDominatedSorting
+
+import numpy as np
+
+# generate random data samples
+F = np.random.random((1000, 2))
+
+# use fast non-dominated sorting
+res = timeit("NonDominatedSorting(method=\"fast_non_dominated_sort\").do(F)", number=10, globals=globals())
+print(f"Fast ND sort takes {res} seconds")
+
+# # use efficient non-dominated sorting with sequential search, this is the default method
+# res = timeit("NonDominatedSorting(method=\"efficient_non_dominated_sort\").do(F, strategy=\"sequential\")", number=10,
+#              globals=globals())
+# print(f"Efficient ND sort with sequential search (ENS-SS) takes {res} seconds")
+#
+#
+# res = timeit("NonDominatedSorting(method=\"efficient_non_dominated_sort\").do(F, strategy=\"binary\")", number=10,
+#              globals=globals())
+# print(f"Efficient ND sort with binary search (ENS-BS) takes {res} seconds")

pymoo/util/function_loader.py

@@ -5,13 +5,21 @@
 
 def get_functions():
     from pymoo.util.nds.fast_non_dominated_sort import fast_non_dominated_sort
+    from pymoo.util.nds.efficient_non_dominated_sort import efficient_non_dominated_sort
+    from pymoo.util.nds.tree_based_non_dominated_sort import tree_based_non_dominated_sort
     from pymoo.decomposition.util import calc_distance_to_weights
     from pymoo.util.misc import calc_perpendicular_distance
 
     FUNCTIONS = {
         "fast_non_dominated_sort": {
             "python": fast_non_dominated_sort, "cython": "pymoo.cython.non_dominated_sorting"
         },
+        "efficient_non_dominated_sort": {
+            "python": efficient_non_dominated_sort, "cython": "pymoo.cython.non_dominated_sorting"
+        },
+        "tree_based_non_dominated_sort": {
+            "python": tree_based_non_dominated_sort, "cython": "pymoo.cython.non_dominated_sorting"
+        },
         "calc_distance_to_weights": {
             "python": calc_distance_to_weights, "cython": "pymoo.cython.decomposition"
         },

pymoo/util/nds/efficient_non_dominated_sort.py

@@ -0,0 +1,141 @@
+from math import floor
+
+import numpy as np
+
+from pymoo.util.dominator import Dominator
+
+
+def efficient_non_dominated_sort(F, strategy="sequential"):
+    """
+    Efficient Non-dominated Sorting (ENS)
+    Parameters
+    ----------
+    F: numpy.ndarray
+        objective values for each individual.
+    strategy: str
+        search strategy, can be "sequential" or "binary".
+    Returns
+    -------
+        indices of the individuals in each front.
+
+    References
+    ----------
+    X. Zhang, Y. Tian, R. Cheng, and Y. Jin,
+    An efficient approach to nondominated sorting for evolutionary multiobjective optimization,
+    IEEE Transactions on Evolutionary Computation, 2015, 19(2): 201-213.
+    """
+    assert (strategy in ["sequential", 'binary']), "Invalid search strategy"
+    N, M = F.shape
+    # sort the rows in F
+    indices = sort_rows(F)
+    F = F[indices]
+    # front ranks for each individual
+    fronts = []  # front with sorted indices
+    _fronts = []  # real fronts
+    for i in range(N):
+        if strategy == 'sequential':
+            k = sequential_search(F, i, fronts)
+        else:
+            k = binary_search(F, i, fronts)
+        if k >= len(fronts):
+            fronts.append([])
+            _fronts.append([])
+        fronts[k].append(i)
+        _fronts[k].append(indices[i])
+    return _fronts
+
+
+def sequential_search(F, i, fronts) -> int:
+    """
+    Find the front rank for the i-th individual through sequential search
+    Parameters
+    ----------
+    F: the objective values
+    i: the index of the individual
+    fronts: individuals in each front
+    """
+    num_found_fronts = len(fronts)
+    k = 0  # the front now checked
+    current = F[i]
+    while True:
+        if num_found_fronts == 0:
+            return 0
+        # solutions in the k-th front, examine in reverse order
+        fk_indices = fronts[k]
+        solutions = F[fk_indices[::-1]]
+        non_dominated = True
+        for f in solutions:
+            relation = Dominator.get_relation(current, f)
+            if relation == -1:
+                non_dominated = False
+                break
+        if non_dominated:
+            return k
+        else:
+            k += 1
+            if k >= num_found_fronts:
+                # move the individual to a new front
+                return num_found_fronts
+
+
+def binary_search(F, i, fronts):
+    """
+    Find the front rank for the i-th individual through binary search
+    Parameters
+    ----------
+    F: the objective values
+    i: the index of the individual
+    fronts: individuals in each front
+    """
+    num_found_fronts = len(fronts)
+    k_min = 0  # the lower bound for checking
+    k_max = num_found_fronts  # the upper bound for checking
+    k = floor((k_max + k_min) / 2 + 0.5)  # the front now checked
+    current = F[i]
+    while True:
+        if num_found_fronts == 0:
+            return 0
+        # solutions in the k-th front, examine in reverse order
+        fk_indices = fronts[k - 1]
+        solutions = F[fk_indices[::-1]]
+        non_dominated = True
+        for f in solutions:
+            relation = Dominator.get_relation(current, f)
+            if relation == -1:
+                non_dominated = False
+                break
+        # binary search
+        if non_dominated:
+            if k == k_min + 1:
+                return k - 1
+            else:
+                k_max = k
+                k = floor((k_max + k_min) / 2 + 0.5)
+        else:
+            k_min = k
+            if k_max == k_min + 1 and k_max < num_found_fronts:
+                return k_max - 1
+            elif k_min == num_found_fronts:
+                return num_found_fronts
+            else:
+                k = floor((k_max + k_min) / 2 + 0.5)
+
+
+def sort_rows(array, order='asc'):
+    """
+    Sort the rows of an array in ascending order.
+    The algorithm will try to use the first column to sort the rows of the given array. If ties occur, it will use the
+    second column, and so on.
+    Parameters
+    ----------
+    array: numpy.ndarray
+        array to be sorted
+    order: str
+        sort order, can be 'asc' or 'desc'
+    Returns
+    -------
+    the indices of the rows in the sorted array.
+    """
+    assert (order in ['asc', 'desc']), "Invalid sort order!"
+    ix = np.lexsort(array.T[::-1])
+    return ix if order == 'asc' else ix[::-1]

pymoo/util/nds/non_dominated_sorting.py

@@ -1,7 +1,7 @@
 import numpy as np
 
-from pymoo.util.function_loader import load_function
 from pymoo.util.dominator import Dominator
+from pymoo.util.function_loader import load_function
 
 
 class NonDominatedSorting:
@@ -11,19 +11,14 @@ def __init__(self, epsilon=0.0, method="fast_non_dominated_sort") -> None:
         self.epsilon = float(epsilon)
         self.method = method
 
-    def do(self, F, return_rank=False, only_non_dominated_front=False, n_stop_if_ranked=None):
+    def do(self, F, return_rank=False, only_non_dominated_front=False, n_stop_if_ranked=None, **kwargs):
         F = F.astype(np.float)
 
         # if not set just set it to a very large values because the cython algorithms do not take None
         if n_stop_if_ranked is None:
             n_stop_if_ranked = int(1e8)
-
-        if self.method == 'fast_non_dominated_sort':
-            func = load_function("fast_non_dominated_sort")
-        else:
-            raise Exception("Unknown non-dominated sorting method: %s" % self.method)
-
-        fronts = func(F, epsilon=self.epsilon)
+        func = load_function(self.method)
+        fronts = func(F, epsilon=self.epsilon, **kwargs)
 
         # convert to numpy array for each front and filter by n_stop_if_ranked if desired
         _fronts = []
@@ -65,5 +60,3 @@ def find_non_dominated(F, _F=None):
     M = Dominator.calc_domination_matrix(F, _F)
     I = np.where(np.all(M >= 0, axis=1))[0]
     return I
-
-

pymoo/util/nds/tree_based_non_dominated_sort.py

@@ -0,0 +1,133 @@
+import weakref
+
+import numpy as np
+
+from pymoo.util.nds.efficient_non_dominated_sort import sort_rows
+
+
+class Tree:
+    '''
+    Implementation of Nary-tree.
+    The source code is modified based on https://github.com/lianemeth/forest/blob/master/forest/NaryTree.py
+
+    Parameters
+    ----------
+    key: object
+        key of the node
+    num_branch: int
+        how many branches in each node
+    children: Iterable[Tree]
+        reference of the children
+    parent: Tree
+        reference of the parent node
+    Returns
+    -------
+        an N-ary tree.
+    '''
+
+    def __init__(self, key, num_branch, children=None, parent=None):
+        self.key = key
+        self.children = children or [None for _ in range(num_branch)]
+
+        self._parent = weakref.ref(parent) if parent else None
+
+    @property
+    def parent(self):
+        if self._parent:
+            return self._parent()
+
+    def __getstate__(self):
+        self._parent = None
+
+    def __setstate__(self, state):
+        self.__dict__ = state
+        for child in self.children:
+            child._parent = weakref.ref(self)
+
+    def traversal(self, visit=None, *args, **kwargs):
+        if visit is not None:
+            visit(self, *args, **kwargs)
+        l = [self]
+        for child in self.children:
+            if child is not None:
+                l += child.traversal(visit, *args, **kwargs)
+        return l
+
+
+def tree_based_non_dominated_sort(F):
+    """
+    Tree-based efficient non-dominated sorting (T-ENS).
+    This algorithm is very efficient in many-objective optimization problems (MaOPs).
+    Parameters
+    ----------
+    F: np.array
+        objective values for each individual.
+    Returns
+    -------
+        indices of the individuals in each front.
+    References
+    ----------
+    X. Zhang, Y. Tian, R. Cheng, and Y. Jin,
+    A decision variable clustering based evolutionary algorithm for large-scale many-objective optimization,
+    IEEE Transactions on Evolutionary Computation, 2018, 22(1): 97-112.
+    """
+    N, M = F.shape
+    # sort the rows in F
+    indices = sort_rows(F)
+    F = F[indices]
+
+    obj_seq = np.argsort(F[:, :0:-1], axis=1) + 1
+
+    k = 0
+
+    forest = []
+
+    left = np.full(N, True)
+    while np.any(left):
+        forest.append(None)
+        for p, flag in enumerate(left):
+            if flag:
+                update_tree(F, p, forest, k, left, obj_seq)
+        k += 1
+
+    # convert forest to fronts
+    fronts = [[] for _ in range(k)]
+    for k, tree in enumerate(forest):
+        fronts[k].extend([indices[node.key] for node in tree.traversal()])
+    return fronts
+
+
+def update_tree(F, p, forest, k, left, obj_seq):
+    _, M = F.shape
+    if forest[k] is None:
+        forest[k] = Tree(key=p, num_branch=M - 1)
+        left[p] = False
+    elif check_tree(F, p, forest[k], obj_seq, True):
+        left[p] = False
+
+
+def check_tree(F, p, tree, obj_seq, add_pos):
+    if tree is None:
+        return True
+
+    N, M = F.shape
+
+    # find the minimal index m satisfying that p[obj_seq[tree.root][m]] < tree.root[obj_seq[tree.root][m]]
+    m = 0
+    while m < M - 1 and F[p, obj_seq[tree.key, m]] >= F[tree.key, obj_seq[tree.key, m]]:
+        m += 1
+
+    # if m not found
+    if m == M - 1:
+        # p is dominated by the solution at the root
+        return False
+    else:
+        for i in range(m + 1):
+            # p is dominated by a solution in the branch of the tree
+            if not check_tree(F, p, tree.children[i], obj_seq, i == m and add_pos):
+                return False
+
+        if tree.children[m] is None and add_pos:
+            # add p to the branch of the tree
+            tree.children[m] = Tree(key=p, num_branch=M - 1)
+        return True