Add an example of regular hypervolume-base algorithm

This algorithm illustrates the use of the “fit_attr” parameters on selection methods.

examples/ga/mo_rhv.py

@@ -0,0 +1,192 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+# Regular Hypervolume-based Algorithm (greedy version)
+
+import array
+import random
+import json
+
+import numpy
+
+from math import sqrt
+
+from deap import algorithms
+from deap import base
+from deap import benchmarks
+from deap.benchmarks.tools import diversity, convergence, hypervolume
+from deap.tools.indicator import hv
+from deap import creator
+from deap import tools
+
+
+creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0))
+# Hypervolume contribution
+creator.create("FitnessHV", base.Fitness, weights=(1.0,))
+creator.create("Individual", array.array, typecode='d',
+               fitness=creator.FitnessMin, fitness_hv=creator.FitnessHV)
+
+toolbox = base.Toolbox()
+
+# Problem definition
+# Functions zdt1, zdt2, zdt3, zdt6 have bounds [0, 1]
+BOUND_LOW, BOUND_UP = 0.0, 1.0
+
+# Functions zdt4 has bounds x1 = [0, 1], xn = [-5, 5], with n = 2, ..., 10
+# BOUND_LOW, BOUND_UP = [0.0] + [-5.0]*9, [1.0] + [5.0]*9
+
+# Functions zdt1, zdt2, zdt3 have 30 dimensions, zdt4 and zdt6 have 10
+NDIM = 30
+
+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)]
+
+
+def hypervolume_contrib(front, **kargs):
+    """Returns the hypervolume contribution of each individual. The provided
+    *front* should be a set of non-dominated individuals having each a
+    :attr:`fitness` attribute.
+    """
+    # Must use wvalues * -1 since hypervolume use implicit minimization
+    # And minimization in deap use max on -obj
+    wobj = numpy.array([ind.fitness.wvalues for ind in front]) * -1
+    ref = kargs.get("ref", None)
+    if ref is None:
+        ref = numpy.max(wobj, axis=0) + 1
+
+    total_hv = hv.hypervolume(wobj, ref)
+
+    def contribution(i):
+        # The contribution of point p_i in point set P
+        # is the hypervolume of P without p_i
+        return total_hv - hv.hypervolume(numpy.concatenate((wobj[:i], wobj[i+1:])), ref)
+
+    # Parallelization note: Cannot pickle local function
+    return map(contribution, range(len(front)))
+
+
+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)
+
+toolbox.register("evaluate", benchmarks.zdt1)
+toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0)
+toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM)
+toolbox.register("sort", tools.sortLogNondominated)
+# Selection is based on HV fitness
+toolbox.register("select", tools.selBest, fit_attr="fitness_hv")
+
+def main(seed=None):
+    random.seed(seed)
+
+    NGEN = 250
+    MU = 100
+    CXPB = 0.9
+
+    stats = tools.Statistics(lambda ind: ind.fitness.values)
+    # 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)
+
+    logbook = tools.Logbook()
+    logbook.header = "gen", "evals", "std", "min", "avg", "max"
+
+    pop = toolbox.population(n=MU)
+
+    # Evaluate the individuals with an invalid fitness
+    invalid_ind = [ind for ind in pop if not ind.fitness.valid]
+    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+    for ind, fit in zip(invalid_ind, fitnesses):
+        ind.fitness.values = fit
+
+    record = stats.compile(pop)
+    logbook.record(gen=0, evals=len(invalid_ind), **record)
+    print(logbook.stream)
+
+    # Begin the generational process
+    for gen in range(1, NGEN):
+        # Vary the population
+        offspring = tools.selRandom(pop, len(pop))
+        offspring = [toolbox.clone(ind) for ind in offspring]
+
+        for ind1, ind2 in zip(offspring[::2], offspring[1::2]):
+            if random.random() <= CXPB:
+                toolbox.mate(ind1, ind2)
+
+            toolbox.mutate(ind1)
+            toolbox.mutate(ind2)
+            del ind1.fitness.values, ind2.fitness.values
+
+        # Evaluate the individuals with an invalid fitness
+        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
+        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+        for ind, fit in zip(invalid_ind, fitnesses):
+            ind.fitness.values = fit
+
+        # Select the next generation population
+        pop = pop + offspring
+        fronts = toolbox.sort(pop, len(pop))
+        chosen = []
+        for i, front in enumerate(fronts):
+            # Move is front to chosen population til it is almost full
+            if len(chosen) + len(front) <= MU:
+                chosen.extend(front)
+            else:
+                # Assign hypervolume contribution to individuals of front that
+                # cannot be completely move over to chosen individuals
+                fitness_hv = hypervolume_contrib(front)
+                for ind, fit_hv in zip(front, fitness_hv):
+                    ind.fitness_hv.values = (fit_hv,)
+                # Fill chosen with best indiviuals from inspect front
+                # (based on hypervolume contribution)
+                chosen.extend(toolbox.select(front, MU - len(chosen)))
+                break
+
+        pop = chosen
+
+        record = stats.compile(pop)
+        logbook.record(gen=gen, evals=len(invalid_ind), **record)
+        print(logbook.stream)
+
+    print("Final population hypervolume is %f" % hypervolume(pop, [11.0, 11.0]))
+
+    return pop, logbook
+
+if __name__ == "__main__":
+    # with open("pareto_front/zdt1_front.json") as optimal_front_data:
+    #     optimal_front = json.load(optimal_front_data)
+    # Use 500 of the 1000 points in the json file
+    # optimal_front = sorted(optimal_front[i] for i in range(0, len(optimal_front), 2))
+
+    pop, stats = main()
+    # pop.sort(key=lambda x: x.fitness.values)
+
+    # print(stats)
+    # print("Convergence: ", convergence(pop, optimal_front))
+    # print("Diversity: ", diversity(pop, optimal_front[0], optimal_front[-1]))
+
+    # import matplotlib.pyplot as plt
+    # import numpy
+
+    # front = numpy.array([ind.fitness.values for ind in pop])
+    # optimal_front = numpy.array(optimal_front)
+    # plt.scatter(optimal_front[:,0], optimal_front[:,1], c="r")
+    # plt.scatter(front[:,0], front[:,1], c="b")
+    # plt.axis("tight")
+    # plt.show()