Start of nsga-3 integration

wujiren · May 27, 2019 · 01130f0 · 01130f0
1 parent f37f955
commit 01130f0
Show file tree

Hide file tree

Showing 5 changed files with 372 additions and 36 deletions.
diff --git a/.gitignore b/.gitignore
@@ -3,6 +3,7 @@ build/
 dist/
 doc/_build/
 .vscode/
+.pytest_cache/
 __pycache__/
 env/
 *.egg-info/

diff --git a/deap/benchmarks/tools.py b/deap/benchmarks/tools.py
@@ -8,6 +8,8 @@
 except ImportError:
     numpy = False
 
+import scipy
+
 try:
     # try importing the C version
     from ..tools._hypervolume import hv
@@ -23,32 +25,32 @@ class translate(object):
     applied to the individual and the resulting list is given to the
     evaluation function. Thus, the evaluation function shall not be expecting
     an individual as it will receive a plain list.
-    
+
     This decorator adds a :func:`translate` method to the decorated function.
     """
     def __init__(self, vector):
         self.vector = vector
-    
+
     def __call__(self, func):
         # wraps is used to combine stacked decorators that would add functions
         @wraps(func)
         def wrapper(individual, *args, **kargs):
             # A subtraction is applied since the translation is applied to the
             # individual and not the function
-            return func([v - t for v, t in zip(individual, self.vector)], 
+            return func([v - t for v, t in zip(individual, self.vector)],
                 *args, **kargs)
         wrapper.translate = self.translate
         return wrapper
-    
+
     def translate(self, vector):
         """Set the current translation to *vector*. After decorating the
         evaluation function, this function will be available directly from
         the function object. ::
-            
+
             @translate([0.25, 0.5, ..., 0.1])
             def evaluate(individual):
                 return sum(individual),
-            
+
             # This will cancel the translation
             evaluate.translate([0.0, 0.0, ..., 0.0])
         """
@@ -63,13 +65,13 @@ class rotate(object):
     list is given to the evaluation function. Thus, the evaluation function
     shall not be expecting an individual as it will receive a plain list
     (numpy.array). The multiplication is done using numpy.
-    
+
     This decorator adds a :func:`rotate` method to the decorated function.
-    
+
     .. note::
-    
+
        A random orthogonal matrix Q can be created via QR decomposition. ::
-           
+
            A = numpy.random.random((n,n))
            Q, _ = numpy.linalg.qr(A)
     """
@@ -93,15 +95,15 @@ def rotate(self, matrix):
         """Set the current rotation to *matrix*. After decorating the
         evaluation function, this function will be available directly from
         the function object. ::
-            
+
             # Create a random orthogonal matrix
             A = numpy.random.random((n,n))
             Q, _ = numpy.linalg.qr(A)
-            
+
             @rotate(Q)
             def evaluate(individual):
                 return sum(individual),
-            
+
             # This will reset rotation to identity
             evaluate.rotate(numpy.identity(n))
         """
@@ -141,26 +143,26 @@ def wrapper(individual, *args, **kargs):
             return tuple(noisy)
         wrapper.noise = self.noise
         return wrapper
-    
+
     def noise(self, noise):
         """Set the current noise to *noise*. After decorating the
         evaluation function, this function will be available directly from
         the function object. ::
-        
+
             prand = functools.partial(random.gauss, mu=0.0, sigma=1.0)
-        
+
             @noise(prand)
             def evaluate(individual):
                 return sum(individual),
-        
+
             # This will remove noise from the evaluation function
             evaluate.noise(None)
         """
         try:
             self.rand_funcs = tuple(noise)
         except TypeError:
             self.rand_funcs = repeat(noise)
-            
+
 class scale(object):
     """Decorator for evaluation functions, it scales the objective function by
     *factor* which should be the same length as the individual size. When
@@ -169,7 +171,7 @@ class scale(object):
     individual and the resulting list is given to the evaluation function.
     Thus, the evaluation function shall not be expecting an individual as it
     will receive a plain list.
-    
+
     This decorator adds a :func:`scale` method to the decorated function.
     """
     def __init__(self, factor):
@@ -181,7 +183,7 @@ def __call__(self, func):
         # wraps is used to combine stacked decorators that would add functions
         @wraps(func)
         def wrapper(individual, *args, **kargs):
-            return func([v * f for v, f in zip(individual, self.factor)], 
+            return func([v * f for v, f in zip(individual, self.factor)],
                 *args, **kargs)
         wrapper.scale = self.scale
         return wrapper
@@ -190,11 +192,11 @@ def scale(self, factor):
         """Set the current scale to *factor*. After decorating the
         evaluation function, this function will be available directly from
         the function object. ::
-            
+
             @scale([0.25, 2.0, ..., 0.1])
             def evaluate(individual):
                 return sum(individual),
-            
+
             # This will cancel the scaling
             evaluate.scale([1.0, 1.0, ..., 1.0])
         """
@@ -213,7 +215,7 @@ class bound(object):
 
     The *type* determines how the attributes are brought back into the valid
     range
-    
+
     This decorator adds a :func:`bound` method to the decorated function.
     """
     def _clip(self, individual):
@@ -245,10 +247,10 @@ def __init__(self, bounds, type):
             self.bound = self._wrap
         elif type == "clip":
             self.bound = self._clip
-        
+
 def diversity(first_front, first, last):
-    """Given a Pareto front `first_front` and the two extreme points of the 
-    optimal Pareto front, this function returns a metric of the diversity 
+    """Given a Pareto front `first_front` and the two extreme points of the
+    optimal Pareto front, this function returns a metric of the diversity
     of the front as explained in the original NSGA-II article by K. Deb.
     The smaller the value is, the better the front is.
     """
@@ -269,13 +271,13 @@ def diversity(first_front, first, last):
     return delta
 
 def convergence(first_front, optimal_front):
-    """Given a Pareto front `first_front` and the optimal Pareto front, 
+    """Given a Pareto front `first_front` and the optimal Pareto front,
     this function returns a metric of convergence
     of the front as explained in the original NSGA-II article by K. Deb.
     The smaller the value is, the closer the front is to the optimal one.
     """
     distances = []
-    
+
     for ind in first_front:
         distances.append(float("inf"))
         for opt_ind in optimal_front:
@@ -285,7 +287,7 @@ def convergence(first_front, optimal_front):
             if dist < distances[-1]:
                 distances[-1] = dist
         distances[-1] = sqrt(distances[-1])
-        
+
     return sum(distances) / len(distances)
 
 
@@ -301,4 +303,7 @@ def hypervolume(front, ref=None):
     wobj = numpy.array([ind.fitness.wvalues for ind in front]) * -1
     if ref is None:
         ref = numpy.max(wobj, axis=0) + 1
-    return hv.hypervolume(wobj, ref)
+    return hv.hypervolume(wobj, ref)
+
+def igd(A, Z):
+    pass
diff --git a/deap/tools/_nsga_3_support.py b/deap/tools/_nsga_3_support.py
@@ -0,0 +1,152 @@
+# NSGA-3 is kindly provided by Luis Marti (IC/UFF) http://lmarti.com
+
+from __future__ import division  # making it work with Python 2.x
+import copy
+import random
+
+import numpy
+from deap import tools
+
+
+class ReferencePoint(list):
+    """A reference point exists in objective space an has a set of individuals
+    associated to it."""
+    def __init__(self, *args):
+        list.__init__(self, *args)
+        self.associations_count = 0
+        self.associations = []
+
+
+def generate_reference_points(num_objs, num_divisions_per_obj=4):
+    """Generates reference points for NSGA-III selection. This code is based on
+    `jMetal NSGA-III implementation <https://github.com/jMetal/jMetal>`_.
+    """
+    def gen_refs_recursive(work_point, num_objs, left, total, depth):
+        if depth == num_objs - 1:
+            work_point[depth] = left/total
+            ref = ReferencePoint(copy.deepcopy(work_point))
+            return [ref]
+        else:
+            res = []
+            for i in range(left):
+                work_point[depth] = i/total
+                res = res + gen_refs_recursive(work_point, num_objs, left-i, total, depth+1)
+            return res
+    return gen_refs_recursive([0]*num_objs, num_objs, num_objs*num_divisions_per_obj,
+                              num_objs*num_divisions_per_obj, 0)
+
+def find_ideal_point(individuals):
+    """Finds the ideal point from a set individuals."""
+    current_ideal = [numpy.infty] * len(individuals[0].fitness.values)
+    for ind in individuals:
+        # Use wvalues to accomodate for maximization and minimization problems.
+        current_ideal = numpy.minimum(current_ideal,
+                                   numpy.multiply(ind.fitness.wvalues, -1))
+    return current_ideal
+
+def find_extreme_points(individuals):
+    """Finds the individuals with extreme values for each objective function."""
+    return [sorted(individuals, key=lambda ind:ind.fitness.wvalues[o] * -1)[-1]
+            for o in range(len(individuals[0].fitness.values))]
+
+def construct_hyperplane(individuals, extreme_points):
+    """Calculates the axis intersects for a set of individuals and its extremes."""
+    def has_duplicate_individuals(individuals):
+        for i in range(len(individuals)):
+            for j in range(i+1, len(individuals)):
+                if individuals[i].fitness.values == individuals[j].fitness.values:
+                    return True
+        return False
+
+    num_objs = len(individuals[0].fitness.values)
+
+    if has_duplicate_individuals(extreme_points):
+        intercepts = [extreme_points[m].fitness.values[m] for m in range(num_objs)]
+    else:
+        b = numpy.ones(num_objs)
+        A = [point.fitness.values for point in extreme_points]
+        x = numpy.linalg.solve(A,b)
+        intercepts = 1/x
+    return intercepts
+
+def normalize_objective(individual, m, intercepts, ideal_point, epsilon=1e-20):
+    """Normalizes an objective."""
+    # Numeric trick present in JMetal implementation.
+    if numpy.abs(intercepts[m]-ideal_point[m] > epsilon):
+        return individual.fitness.values[m] / (intercepts[m]-ideal_point[m])
+    else:
+        return individual.fitness.values[m] / epsilon
+
+def normalize_objectives(individuals, intercepts, ideal_point):
+    """Normalizes individuals using the hyperplane defined by the intercepts as
+    reference. Corresponds to Algorithm 2 of Deb & Jain (2014)."""
+    num_objs = len(individuals[0].fitness.values)
+
+    for ind in individuals:
+        ind.fitness.normalized_values = list([normalize_objective(ind, m,
+                                                                  intercepts, ideal_point)
+                                                                  for m in range(num_objs)])
+    return individuals
+
+def perpendicular_distance(direction, point):
+    k = numpy.dot(direction, point) / numpy.sum(numpy.power(direction, 2))
+    d = numpy.sum(numpy.power(numpy.subtract(numpy.multiply(direction, [k] * len(direction)), point) , 2))
+    return numpy.sqrt(d)
+
+def associate(individuals, reference_points):
+    """Associates individuals to reference points and calculates niche number.
+    Corresponds to Algorithm 3 of Deb & Jain (2014)."""
+    pareto_fronts = tools.sortLogNondominated(individuals, len(individuals))
+    num_objs = len(individuals[0].fitness.values)
+
+    for ind in individuals:
+        rp_dists = [(rp, perpendicular_distance(ind.fitness.normalized_values, rp))
+                    for rp in reference_points]
+        best_rp, best_dist = sorted(rp_dists, key=lambda rpd:rpd[1])[0]
+        ind.reference_point = best_rp
+        ind.ref_point_distance = best_dist
+        best_rp.associations_count +=1 # update de niche number
+        best_rp.associations += [ind]
+
+def selNiching(individuals, k):
+    """Secondary niched selection based on reference points. Corresponds to
+    steps 13-17 of Algorithm 1 and to Algorithm 4."""
+    if len(individuals) == k:
+        return individuals
+
+    #individuals = copy.deepcopy(individuals)
+
+    ideal_point = find_ideal_point(individuals)
+    extremes = find_extreme_points(individuals)
+    intercepts = construct_hyperplane(individuals, extremes)
+    normalize_objectives(individuals, intercepts, ideal_point)
+
+    reference_points = generate_reference_points(len(individuals[0].fitness.values))
+
+    associate(individuals, reference_points)
+
+    res = []
+    while len(res) < k:
+        min_assoc_rp = min(reference_points, key=lambda rp: rp.associations_count)
+        min_assoc_rps = [rp for rp in reference_points if rp.associations_count == min_assoc_rp.associations_count]
+        chosen_rp = min_assoc_rps[random.randint(0, len(min_assoc_rps)-1)]
+
+        #print('Rps',min_assoc_rp.associations_count, chosen_rp.associations_count, len(min_assoc_rps))
+
+        associated_inds = chosen_rp.associations
+
+        if chosen_rp.associations:
+            if chosen_rp.associations_count == 0:
+                sel = min(chosen_rp.associations, key=lambda ind: ind.ref_point_distance)
+            else:
+                sel = chosen_rp.associations[random.randint(0, len(chosen_rp.associations)-1)]
+            res += [sel]
+            chosen_rp.associations.remove(sel)
+            chosen_rp.associations_count += 1
+            individuals.remove(sel)
+        else:
+            reference_points.remove(chosen_rp)
+    return res
+
+
+__all__ = ["selNiching"]