1#    This file is part of DEAP.
2#
3#    DEAP is free software: you can redistribute it and/or modify
4#    it under the terms of the GNU Lesser General Public License as
5#    published by the Free Software Foundation, either version 3 of
6#    the License, or (at your option) any later version.
7#
8#    DEAP is distributed in the hope that it will be useful,
9#    but WITHOUT ANY WARRANTY; without even the implied warranty of
10#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11#    GNU Lesser General Public License for more details.
12#
13#    You should have received a copy of the GNU Lesser General Public
14#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
15
16import array
17import random
18import json
19
20import numpy
21
22from math import sqrt
23
24from deap import algorithms
25from deap import base
26from deap import benchmarks
27from deap.benchmarks.tools import diversity, convergence, hypervolume
28from deap import creator
29from deap import tools
30
31creator.create("FitnessMin", base.Fitness, weights=(-1.0, -1.0))
32creator.create("Individual", array.array, typecode='d', fitness=creator.FitnessMin)
33
34toolbox = base.Toolbox()
35
36# Problem definition
37# Functions zdt1, zdt2, zdt3, zdt6 have bounds [0, 1]
38BOUND_LOW, BOUND_UP = 0.0, 1.0
39
40# Functions zdt4 has bounds x1 = [0, 1], xn = [-5, 5], with n = 2, ..., 10
41# BOUND_LOW, BOUND_UP = [0.0] + [-5.0]*9, [1.0] + [5.0]*9
42
43# Functions zdt1, zdt2, zdt3 have 30 dimensions, zdt4 and zdt6 have 10
44NDIM = 30
45
46def uniform(low, up, size=None):
47    try:
48        return [random.uniform(a, b) for a, b in zip(low, up)]
49    except TypeError:
50        return [random.uniform(a, b) for a, b in zip([low] * size, [up] * size)]
51
52toolbox.register("attr_float", uniform, BOUND_LOW, BOUND_UP, NDIM)
53toolbox.register("individual", tools.initIterate, creator.Individual, toolbox.attr_float)
54toolbox.register("population", tools.initRepeat, list, toolbox.individual)
55
56toolbox.register("evaluate", benchmarks.zdt1)
57toolbox.register("mate", tools.cxSimulatedBinaryBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0)
58toolbox.register("mutate", tools.mutPolynomialBounded, low=BOUND_LOW, up=BOUND_UP, eta=20.0, indpb=1.0/NDIM)
59toolbox.register("select", tools.selNSGA2)
60
61def main(seed=None):
62    random.seed(seed)
63
64    NGEN = 250
65    MU = 100
66    CXPB = 0.9
67
68    stats = tools.Statistics(lambda ind: ind.fitness.values)
69    # stats.register("avg", numpy.mean, axis=0)
70    # stats.register("std", numpy.std, axis=0)
71    stats.register("min", numpy.min, axis=0)
72    stats.register("max", numpy.max, axis=0)
73
74    logbook = tools.Logbook()
75    logbook.header = "gen", "evals", "std", "min", "avg", "max"
76
77    pop = toolbox.population(n=MU)
78
79    # Evaluate the individuals with an invalid fitness
80    invalid_ind = [ind for ind in pop if not ind.fitness.valid]
81    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
82    for ind, fit in zip(invalid_ind, fitnesses):
83        ind.fitness.values = fit
84
85    # This is just to assign the crowding distance to the individuals
86    # no actual selection is done
87    pop = toolbox.select(pop, len(pop))
88
89    record = stats.compile(pop)
90    logbook.record(gen=0, evals=len(invalid_ind), **record)
91    print(logbook.stream)
92
93    # Begin the generational process
94    for gen in range(1, NGEN):
95        # Vary the population
96        offspring = tools.selTournamentDCD(pop, len(pop))
97        offspring = [toolbox.clone(ind) for ind in offspring]
98
99        for ind1, ind2 in zip(offspring[::2], offspring[1::2]):
100            if random.random() <= CXPB:
101                toolbox.mate(ind1, ind2)
102
103            toolbox.mutate(ind1)
104            toolbox.mutate(ind2)
105            del ind1.fitness.values, ind2.fitness.values
106
107        # Evaluate the individuals with an invalid fitness
108        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
109        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
110        for ind, fit in zip(invalid_ind, fitnesses):
111            ind.fitness.values = fit
112
113        # Select the next generation population
114        pop = toolbox.select(pop + offspring, MU)
115        record = stats.compile(pop)
116        logbook.record(gen=gen, evals=len(invalid_ind), **record)
117        print(logbook.stream)
118
119    print("Final population hypervolume is %f" % hypervolume(pop, [11.0, 11.0]))
120
121    return pop, logbook
122
123if __name__ == "__main__":
124    # with open("pareto_front/zdt1_front.json") as optimal_front_data:
125    #     optimal_front = json.load(optimal_front_data)
126    # Use 500 of the 1000 points in the json file
127    # optimal_front = sorted(optimal_front[i] for i in range(0, len(optimal_front), 2))
128
129    pop, stats = main()
130    # pop.sort(key=lambda x: x.fitness.values)
131
132    # print(stats)
133    # print("Convergence: ", convergence(pop, optimal_front))
134    # print("Diversity: ", diversity(pop, optimal_front[0], optimal_front[-1]))
135
136    # import matplotlib.pyplot as plt
137    # import numpy
138
139    # front = numpy.array([ind.fitness.values for ind in pop])
140    # optimal_front = numpy.array(optimal_front)
141    # plt.scatter(optimal_front[:,0], optimal_front[:,1], c="r")
142    # plt.scatter(front[:,0], front[:,1], c="b")
143    # plt.axis("tight")
144    # plt.show()
145