Merge pull request DEAP#203 from ericjster/master

Fix ClosestValidPenalty to use distance

deap/tests/test_algorithms.py

@@ -31,7 +31,7 @@
 FITCLSNAME = "FIT_TYPE"
 INDCLSNAME = "IND_TYPE"
 
-HV_THRESHOLD = 119.0
+HV_THRESHOLD = 116.0        # 120.777 is Optimal value
 
 
 def setup_func_single_obj():
@@ -115,6 +115,10 @@ def test_nsga2():
 
     assert hv > HV_THRESHOLD, "Hypervolume is lower than expected %f < %f" % (hv, HV_THRESHOLD)
 
+    for ind in pop:
+        assert not (any(numpy.asarray(ind) < BOUND_LOW) or any(numpy.asarray(ind) > BOUND_UP))
+
+
 @unittest.skipIf(platform.python_implementation() == "PyPy", "PyPy has no support for eigen decomposition.")
 @with_setup(setup_func_multi_obj_numpy, teardown_func)
 def test_mo_cma_es():
@@ -141,12 +145,14 @@ def valid(individual):
     MU, LAMBDA = 10, 10
     NGEN = 500
 
+    numpy.random.seed(128)
+
     # The MO-CMA-ES algorithm takes a full population as argument
     population = [creator.__dict__[INDCLSNAME](x) for x in numpy.random.uniform(BOUND_LOW, BOUND_UP, (MU, NDIM))]
 
     toolbox = base.Toolbox()
     toolbox.register("evaluate", benchmarks.zdt1)
-    toolbox.decorate("evaluate", tools.ClosestValidPenalty(valid, closest_feasible, 1.0e-6, distance))
+    toolbox.decorate("evaluate", tools.ClosestValidPenalty(valid, closest_feasible, 1.0e+6, distance))
 
     for ind in population:
         ind.fitness.values = toolbox.evaluate(ind)
@@ -168,5 +174,17 @@ def valid(individual):
         # Update the strategy with the evaluated individuals
         toolbox.update(population)
 
+    # Note that we use a penalty to guide the search to feasible solutions,
+    # but there is no guarantee that individuals are valid.
+    # We expect the best individuals will be within bounds or very close.
+    num_valid = 0
+    for ind in strategy.parents:
+        dist = distance(closest_feasible(ind), ind)
+        if numpy.isclose(dist, 0.0, rtol=1.e-5, atol=1.e-5):
+            num_valid += 1
+    assert num_valid >= len(strategy.parents)
+
+    # Note that NGEN=500 is enough to get consistent hypervolume > 116,
+    # but not 119. More generations would help but would slow down testing.
     hv = hypervolume(strategy.parents, [11.0, 11.0])
     assert hv > HV_THRESHOLD, "Hypervolume is lower than expected %f < %f" % (hv, HV_THRESHOLD)

deap/tools/constraint.py

@@ -117,11 +117,12 @@ def wrapper(individual, *args, **kwargs):
 
             dists = tuple(0 for w in individual.fitness.weights)
             if self.dist_fct is not None:
-                dist = self.dist_fct(f_ind, individual)
+                dists = self.dist_fct(f_ind, individual)
                 if not isinstance(dists, Sequence):
                     dists = repeat(dists)
 
-            # print("returned", tuple(f - w * self.alpha * dist for f, w in zip(f_fbl, weights)))
+            # print("penalty ", tuple(  - w * self.alpha * d for f, w, d in zip(f_fbl, weights, dists)))
+            # print("returned", tuple(f - w * self.alpha * d for f, w, d in zip(f_fbl, weights, dists)))
             return tuple(f - w * self.alpha * d for f, w, d in zip(f_fbl, weights, dists))
 
         return wrapper

examples/es/cma_mo.py

@@ -29,6 +29,7 @@
 # ZDT1, ZDT2, DTLZ2
 MIN_BOUND = numpy.zeros(N)
 MAX_BOUND = numpy.ones(N)
+EPS_BOUND = 2.e-5
 
 # Kursawe
 # MIN_BOUND = numpy.zeros(N) - 5
@@ -54,9 +55,15 @@ def valid(individual):
         return False
     return True
 
+def close_valid(individual):
+    """Determines if the individual is close to valid."""
+    if any(individual < MIN_BOUND-EPS_BOUND) or any(individual > MAX_BOUND+EPS_BOUND):
+        return False
+    return True
+
 toolbox = base.Toolbox()
 toolbox.register("evaluate", benchmarks.zdt1)
-toolbox.decorate("evaluate", tools.ClosestValidPenalty(valid, closest_feasible, 1.0e-6, distance))
+toolbox.decorate("evaluate", tools.ClosestValidPenalty(valid, closest_feasible, 1.0e+6, distance))
 
 def main():
     # The cma module uses the numpy random number generator
@@ -65,6 +72,7 @@ def main():
     MU, LAMBDA = 10, 10
     NGEN = 500
     verbose = True
+    create_plot = False
 
     # The MO-CMA-ES algorithm takes a full population as argument
     population = [creator.Individual(x) for x in (numpy.random.uniform(0, 1, (MU, N)))]
@@ -83,6 +91,8 @@ def main():
     logbook = tools.Logbook()
     logbook.header = ["gen", "nevals"] + (stats.fields if stats else [])
 
+    fitness_history = []
+
     for gen in range(NGEN):
         # Generate a new population
         population = toolbox.generate()
@@ -91,6 +101,7 @@ def main():
         fitnesses = toolbox.map(toolbox.evaluate, population)
         for ind, fit in zip(population, fitnesses):
             ind.fitness.values = fit
+            fitness_history.append(fit)
 
         # Update the strategy with the evaluated individuals
         toolbox.update(population)
@@ -102,22 +113,55 @@ def main():
 
     if verbose:
         print("Final population hypervolume is %f" % hypervolume(strategy.parents, [11.0, 11.0]))
-
-    # import matplotlib.pyplot as plt
-
-    # valid_front = numpy.array([ind.fitness.values for ind in strategy.parents if valid(ind)])
-    # invalid_front = numpy.array([ind.fitness.values for ind in strategy.parents if not valid(ind)])
-
-    # fig = plt.figure()
-
-    # if len(valid_front) > 0:
-    #     plt.scatter(valid_front[:,0], valid_front[:,1], c="g")
 
-    # if len(invalid_front) > 0:
-    #     plt.scatter(invalid_front[:,0], invalid_front[:,1], c="r")
+        # Note that we use a penalty to guide the search to feasible solutions,
+        # but there is no guarantee that individuals are valid.
+        # We expect the best individuals will be within bounds or very close.
+        num_valid = 0
+        for ind in strategy.parents:
+            dist = distance(closest_feasible(ind), ind)
+            if numpy.isclose(dist, 0.0, rtol=1.e-5, atol=1.e-5):
+                num_valid += 1
+        print("Number of valid individuals is %d/%d" % (num_valid, len(strategy.parents)))
+
+        print("Final population:")
+        print(numpy.asarray(strategy.parents))
+
+    if create_plot:
+        interactive = 0
+        if not interactive:
+            import matplotlib as mpl_tmp
+            mpl_tmp.use('Agg')   # Force matplotlib to not use any Xwindows backend.
+        import matplotlib.pyplot as plt
+
+        fig = plt.figure()
+        plt.title("Multi-objective minimization via MO-CMA-ES")
+        plt.xlabel("First objective (function) to minimize")
+        plt.ylabel("Second objective (function) to minimize")
+
+        # Limit the scale because our history values include the penalty.
+        plt.xlim((-0.1, 1.20))
+        plt.ylim((-0.1, 1.20))
+
+        # Plot all history. Note the values include the penalty.
+        fitness_history = numpy.asarray(fitness_history)
+        plt.scatter(fitness_history[:,0], fitness_history[:,1],
+            facecolors='none', edgecolors="lightblue")
+
+        valid_front = numpy.array([ind.fitness.values for ind in strategy.parents if close_valid(ind)])
+        invalid_front = numpy.array([ind.fitness.values for ind in strategy.parents if not close_valid(ind)])
+
+        if len(valid_front) > 0:
+            plt.scatter(valid_front[:,0], valid_front[:,1], c="g")
+        if len(invalid_front) > 0:
+            plt.scatter(invalid_front[:,0], invalid_front[:,1], c="r")
+
+        if interactive:
+            plt.show()
+        else:
+            print("Writing cma_mo.png")
+            plt.savefig("cma_mo.png")
 
-    # plt.show()
-
     return strategy.parents
 
 if __name__ == "__main__":