Added nsga3 example

deap/tools/emo.py

@@ -532,7 +532,7 @@ def find_extreme_points(fitnesses, best_point, extreme_points=None):
 def find_intercepts(extreme_points, best_point, current_worst):
     """Find intercepts between the hyperplane and each axis with
     the ideal point as origin."""
-    # Construct hyperplane
+    # Construct hyperplane sum(f_i^n) = 1
     b = numpy.ones(extreme_points.shape[1])
     A = extreme_points - best_point
     try:

doc/code/examples/nsga3_ref_points.py

@@ -0,0 +1,28 @@
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+from deap import tools
+
+NOBJ = 3
+P = 12
+
+fig = plt.figure(figsize=(7, 7))
+ax = fig.add_subplot(111, projection="3d")
+
+# the coordinate origin (black + sign)
+ax.scatter(0, 0, 0, c="k", marker="+", s=100)
+
+# reference points (gray)
+ref_points = tools.uniform_reference_points(NOBJ, P)
+for rp in ref_points:
+    ax.scatter(rp[0], rp[1], rp[2], marker="o", color="darkblue", s=48)
+
+# final figure details
+ax.set_xlabel("$f_1(\mathbf{x})$", fontsize=15)
+ax.set_ylabel("$f_2(\mathbf{x})$", fontsize=15)
+ax.set_zlabel("$f_3(\mathbf{x})$", fontsize=15)
+ax.view_init(elev=11, azim=-25)
+ax.autoscale(tight=True)
+plt.tight_layout()
+
+plt.show()

doc/conf.py

@@ -45,7 +45,7 @@
 source_suffix = '.rst'
 
 # The encoding of source files.
-#source_encoding = 'utf-8'
+source_encoding = 'utf-8'
 
 # The master toctree document.
 master_doc = 'index'

doc/examples/index.rst

@@ -19,6 +19,7 @@ Genetic Algorithm (GA)
 	ga_onemax_numpy
 	ga_knapsack
 	coev_coop
+	nsga3
 
 .. _gpexamples:
 
@@ -60,6 +61,6 @@ Estimation of Distribution Algorithms (EDA)
 
 .. toctree::
    :maxdepth: 1
-   
+
    eda
 

doc/examples/nsga3.rst

@@ -0,0 +1,90 @@
+.. _nsga-3:
+
+======================================================
+Non-dominated Sorting Genetic Algorithm III (NSGA-III)
+======================================================
+The Non-dominated Sorting Genetic Algorithm III (NSGA-III) [Deb2014]_
+is implemented in the :func:`deap.tools.selNSGA3` function. This example
+shows how it can be used in DEAP for many objective optimization.
+
+Problem Definition
+------------------
+First we need to define the problem we want to work on. We will use
+the first problem tested in the paper, 3 objectives DTLZ2 with ``k = 10``
+and ``p = 12``. We will use `pymop <https://github.com/msu-coinlab/pymop>`_
+for problem implementation as it provides the exact Pareto front that we
+will use later for computing the performance of the algorithm.
+
+.. literalinclude:: /../examples/ga/nsga3.py
+   :start-after: # Problem definition
+   :end-before: ##
+
+Algorithm Parameters
+--------------------
+Then we define the various parameters for the algorithm, including the
+population size set to the first multiple of 4 greater than H, the
+number of generations and variation probabilities.
+
+.. literalinclude:: /../examples/ga/nsga3.py
+   :start-after: # Algorithm parameters
+   :end-before: ##
+
+Classes and Tools
+-----------------
+Next, NSGA-III selection requires a reference point set. The reference point
+set serves to guide the evolution into creating a uniform Pareto front in
+the objective space.
+
+.. literalinclude:: /../examples/ga/nsga3.py
+   :start-after: # Create uniform reference point
+   :lines: 1
+
+The next figure shows an example of reference point set with ``p = 12``.
+The cross represents the the utopian point (0, 0, 0).
+
+.. plot:: code/examples/nsga3_ref_points.py
+
+As in any DEAP program, we need to populate the creator with the type
+of individual we require for our optimization. In this case we will use
+a basic list genotype and minimization fitness.
+
+.. literalinclude:: /../examples/ga/nsga3.py
+   :start-after: # Create classes
+   :end-before: ##
+
+Moreover, we need to populate the evolutionary toolbox with initialization,
+variation and selection operators. Note how we provide the reference point
+set to the NSGA-III selection scheme.
+
+.. literalinclude:: /../examples/ga/nsga3.py
+   :start-after: # Toolbox initialization
+   :end-before: ##
+
+Evolution
+---------
+The main part of the evolution is mostly similar to any other
+DEAP example. The algorithm used is close to the
+:func:`~deap.algorithms.eaSimple` algorithm as crossover and mutation are
+applied to every individual (see variation probabilities above). However,
+the selection is made from the parent and offspring populations instead
+of completely replacing the parents with the offspring.
+
+.. literalinclude:: /../examples/ga/nsga3.py
+   :pyobject: main
+
+Finally, we can have a look at the final population
+
+.. image:: /_images/nsga3.png
+   :align: center
+
+Conclusion
+----------
+That's it for the NSGA-III algorithm using DEAP, now you can leverage
+the power of many-objective optimization with DEAP. If you're interrested,
+you can now change the evaluation function and try applying it to
+your own problem.
+
+.. [Deb2014] Deb, K., & Jain, H. (2014). An Evolutionary Many-Objective Optimization
+   Algorithm Using Reference-Point-Based Nondominated Sorting Approach,
+   Part I: Solving Problems With Box Constraints. IEEE Transactions on
+   Evolutionary Computation, 18(4), 577-601. doi:10.1109/TEVC.2013.2281535.

examples/ga/nsga3.py

@@ -1,8 +1,9 @@
 from math import factorial
 import random
 
+import matplotlib.pyplot as plt
 import numpy
-from pymop.factory import get_problem
+import pymop.factory
 
 from deap import algorithms
 from deap import base
@@ -11,35 +12,40 @@
 from deap import tools
 
 # Problem definition
-PROBLEM = "dtlz1"
+PROBLEM = "dtlz2"
 NOBJ = 3
-K = 5
+K = 10
 NDIM = NOBJ + K - 1
 P = 12
 H = factorial(NOBJ + P - 1) / (factorial(P) * factorial(NOBJ - 1))
+BOUND_LOW, BOUND_UP = 0.0, 1.0
+problem = pymop.factory.get_problem(PROBLEM, n_var=NDIM, n_obj=NOBJ)
+##
+
+# Algorithm parameters
 MU = int(H + (4 - H % 4))
 NGEN = 400
 CXPB = 1.0
 MUTPB = 1.0
+##
 
-# Problems dtlz have bounds [0, 1]
-BOUND_LOW, BOUND_UP = 0.0, 1.0
-
+# Create uniform reference point
 ref_points = tools.uniform_reference_points(NOBJ, P)
 
+# Create classes
 creator.create("FitnessMin", base.Fitness, weights=(-1.0,) * NOBJ)
 creator.create("Individual", list, fitness=creator.FitnessMin)
+##
 
-problem = get_problem(PROBLEM, n_var=NDIM, n_obj=NOBJ)
-
-toolbox = base.Toolbox()
 
+# Toolbox initialization
 def uniform(low, up, size=None):
     try:
         return [random.uniform(a, b) for a, b in zip(low, up)]
     except TypeError:
         return [random.uniform(a, b) for a, b in zip([low] * size, [up] * size)]
 
+toolbox = base.Toolbox()
 toolbox.register("attr_float", uniform, BOUND_LOW, BOUND_UP, NDIM)
 toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.attr_float)
 toolbox.register("population", tools.initRepeat, list, toolbox.individual)
@@ -48,13 +54,16 @@ def uniform(low, up, size=None):
 toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=30.0)
 toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM)
 toolbox.register("select", tools.selNSGA3, ref_points=ref_points)
+##
+
 
 def main(seed=None):
     random.seed(seed)
 
+    # Initialize statistics object
     stats = tools.Statistics(lambda ind: ind.fitness.values)
-    # stats.register("avg", numpy.mean, axis=0)
-    # stats.register("std", numpy.std, axis=0)
+    stats.register("avg", numpy.mean, axis=0)
+    stats.register("std", numpy.std, axis=0)
     stats.register("min", numpy.min, axis=0)
     stats.register("max", numpy.max, axis=0)
 
@@ -69,10 +78,7 @@ def main(seed=None):
     for ind, fit in zip(invalid_ind, fitnesses):
         ind.fitness.values = fit
 
-    # This is just to assign the crowding distance to the individuals
-    # no actual selection is done
-    pop = toolbox.select(pop, len(pop))
-
+    # Compile statistics about the population
     record = stats.compile(pop)
     logbook.record(gen=0, evals=len(invalid_ind), **record)
     print(logbook.stream)
@@ -87,20 +93,20 @@ def main(seed=None):
         for ind, fit in zip(invalid_ind, fitnesses):
             ind.fitness.values = fit
 
-        # Select the next generation population
+        # Select the next generation population from parents and offspring
         pop = toolbox.select(pop + offspring, MU)
 
+        # Compile statistics about the new population
         record = stats.compile(pop)
         logbook.record(gen=gen, evals=len(invalid_ind), **record)
         print(logbook.stream)
 
     return pop, logbook
 
-if __name__ == "__main__":
-    pf = problem.pareto_front(ref_points)
 
+if __name__ == "__main__":
     pop, stats = main()
-
     pop_fit = numpy.array([ind.fitness.values for ind in pop])
 
-    print(igd(pop_fit, pf)) # 0.033909687899703514
+    pf = problem.pareto_front(ref_points)
+    print(igd(pop_fit, pf))