Changed hv computation to Monte-Carlo HypE algorithm

Author: Tanguy Damart

bluepyopt/deapext/CMA_MO.py

@@ -32,34 +32,36 @@
 
 from .stoppingCriteria import MaxNGen
 from . import utils
+from . import hype
 
 logger = logging.getLogger('__main__')
 
-from deap.tools._hypervolume import hv as hv_c
 
+def get_hyped(pop):
+    # Cap the obj at 250
+    points = numpy.array([ind.fitness.values for ind in pop])
+    points[points > 250.] = 250.
+    lbounds = numpy.min(points, axis=0)
+    ubounds = numpy.max(points, axis=0)
 
-def get_hv(to_evaluate):
-    i = to_evaluate[0]
-    wobj = to_evaluate[1]
-    ref = to_evaluate[2]
-    return hv_c.hypervolume(numpy.concatenate((wobj[:i], wobj[i + 1:])), ref)
+    # Remove the dimensions that do not show any improvement
+    to_remove = []
+    for i, (lb, ub) in enumerate(zip(lbounds, ubounds)):
+        if lb >= 240:
+            to_remove.append(i)
+            points = numpy.delete(points, to_remove, axis=1)
+            lbounds = numpy.delete(lbounds, to_remove)
+            ubounds = numpy.delete(ubounds, to_remove)
 
+    # Rescale the objective space
+    points = (points - lbounds) / numpy.max(ubounds.flatten())
+    ubounds = numpy.max(points, axis=0) + 2.
 
-def contribution(to_evaluate):
-    def _reduce_method(meth):
-        """Overwrite reduce"""
-        return (getattr, (meth.__self__, meth.__func__.__name__))
-
-    import copyreg
-    import types
-    copyreg.pickle(types.MethodType, _reduce_method)
-    import pebble
-
-    with pebble.ProcessPool(max_tasks=1) as pool:
-        tasks = pool.schedule(get_hv, kwargs={'to_evaluate': to_evaluate})
-        response = tasks.result()
-
-    return response
+    hv = hype.hypeIndicatorSampled(points=points,
+                                   bounds=ubounds,
+                                   k=5,
+                                   nrOfSamples=200000)
+    return hv
 
 
 class CMA_MO(cma.StrategyMultiObjective):
@@ -72,6 +74,7 @@ def __init__(self,
                  max_ngen,
                  IndCreator,
                  RandIndCreator,
+                 weight_hv=0.5,
                  map_function=None,
                  use_scoop=False):
         """Constructor
@@ -105,6 +108,8 @@ def __init__(self,
 
         self.population = []
         self.problem_size = len(starters[0])
+        
+        self.weight_hv = weight_hv
 
         self.map_function = map_function
         self.use_scoop = use_scoop
@@ -130,31 +135,6 @@ def __init__(self,
 
         self.stopping_conditions = [MaxNGen(max_ngen)]
 
-    def hyper_volume(self, front):
-        """Compute the hypervolume contribution of each individual"""
-        wobj = numpy.array([ind.fitness.values for ind in front])
-        obj_ranges = (numpy.max(wobj, axis=0) - numpy.min(wobj, axis=0))
-        ref = numpy.max(wobj, axis=0) + 1
-
-        # Above 23 dimension, the hypervolume computation is too slow,
-        # we settle for the 23 dimension showing the largest range of values
-        max_ndim = 23
-        if len(ref) > max_ndim:
-            idxs = list(range(len(ref)))
-            idxs = [idxs[k] for k in numpy.argsort(obj_ranges)]
-            idxs = idxs[::-1]
-            idxs = idxs[:max_ndim]
-            wobj = wobj[:, idxs]
-            ref = ref[idxs]
-
-        # Prepare the data and send it to multiprocess
-        to_evaluate = []
-        for i in range(len(front)):
-            to_evaluate.append([i, numpy.copy(wobj), numpy.copy(ref)])
-        contrib_values = self.map_function(contribution, to_evaluate)
-
-        return list(contrib_values)
-
     def _select(self, candidates):
         """Select the best candidates of the population
 
@@ -163,36 +143,60 @@ def _select(self, candidates):
         we rely on the hypervolume for this front. The remaining fronts are
         explicitly not chosen"""
 
-        if len(candidates) <= self.mu:
-            return candidates, []
-
-        pareto_fronts = deap.tools.sortLogNondominated(candidates,
-                                                       len(candidates))
-
-        chosen = list()
-        mid_front = None
-        not_chosen = list()
-
-        full = False
-        for front in pareto_fronts:
-            if len(chosen) + len(front) <= self.mu and not full:
-                chosen += front
-            elif mid_front is None and len(chosen) < self.mu:
-                mid_front = front
-                # With this front, we selected enough individuals
-                full = True
-            else:
-                not_chosen += front
-
-        # Hypervolume contribution to get the best candidates on the remaining
-        # front
-        k = self.mu - len(chosen)
-        if k > 0:
-            hyperv = self.hyper_volume(mid_front)
-            _ = [mid_front[k] for k in numpy.argsort(hyperv)]
-            chosen += _[:k]
-            not_chosen += _[k:]
+        #if len(candidates) <= self.mu:
+        #    return candidates, []
+
+        #pareto_fronts = deap.tools.sortLogNondominated(candidates,
+        #                                               len(candidates))
+
+        #chosen = list()
+        #mid_front = None
+        #not_chosen = list()
+
+        #full = False
+        #for front in pareto_fronts:
+        #    if len(chosen) + len(front) <= self.mu and not full:
+        #        chosen += front
+        #    elif mid_front is None and len(chosen) < self.mu:
+        #        mid_front = front
+        #        # With this front, we selected enough individuals
+        #        full = True
+        #    else:
+        #        not_chosen += front
+
+        # HypE contribution to get the best candidates on the remaining front
+        #k = self.mu - len(chosen)
+        #if k > 0:
+        #    contribution = get_hyped(mid_front)
+        #    print(contribution)
+        #    idxs = numpy.argsort(contribution)
+        #    ordered_front = [mid_front[i] for i in idxs[::-1]]
+        #    chosen += ordered_front[:k]
+        #    not_chosen += ordered_front[k:]
+        
+        if self.weight_hv == 0.: 
+            fit = [numpy.sum(ind.fitness.values) for ind in candidates]
+            idx_fit = list(numpy.argsort(fit))
+            idx_scores = idx_fit[:]
+
+        elif self.weight_hv == 1.:
+            hv = get_hyped(candidates)
+            idx_hv = list(numpy.argsort(hv))[::-1] 
+            idx_scores = idx_hv[:] 
 
+        else:
+            hv = get_hyped(candidates)
+            idx_hv = list(numpy.argsort(hv))[::-1] 
+            fit = [numpy.sum(ind.fitness.values) for ind in candidates]
+            idx_fit = list(numpy.argsort(fit))
+            scores = []
+            for i in range(len(candidates)):
+                score = (self.weight_hv * idx_hv.index(i)) + ((1.-self.weight_hv) * idx_fit.index(i))
+                scores.append(score)
+            idx_scores = list(numpy.argsort(scores))
+
+        chosen = [candidates[i] for i in idx_scores[:self.mu]]
+        not_chosen = [candidates[i] for i in idx_scores[self.mu:]]
         return chosen, not_chosen
 
     def get_population(self, to_space):

bluepyopt/deapext/hype.py

@@ -0,0 +1,85 @@
+import numpy
+
+def hypesub(l, A, actDim, bounds, pvec, alpha, k):
+
+    h = numpy.zeros(l)
+    i = numpy.argsort(A[:, actDim - 1])
+    S = A[i]
+    pvec = pvec[i]
+
+    for i in range(1, S.shape[0] + 1):
+
+        if i < S.shape[0]:
+            extrusion = S[i, actDim - 1] - S[i - 1, actDim - 1]
+        else:
+            extrusion = bounds[actDim - 1] - S[i - 1, actDim - 1]
+
+        if actDim == 1:
+            if i > k:
+                break
+            if alpha[i - 1] >= 0:
+                h[pvec[0:i]] += extrusion * alpha[i - 1]
+        elif extrusion > 0.:
+            h += extrusion * hypesub(l, S[0:i, :], actDim - 1, bounds,
+                                     pvec[0:i], alpha, k)
+
+    return h
+
+
+def hypeIndicatorExact(points, bounds, k):
+    """
+    points: objectives (to be minimized),
+    bounds: reference point,
+    k: parameter of HypE
+
+    Example: scores = hypeIndicatorExact([[1., 3.], [3., 1.]], [4., 4.], 1)
+    """
+
+    Ps = points.shape[0]
+    if k < 0:
+        k = Ps
+    actDim = points.shape[1]
+    pvec = numpy.arange(points.shape[0])
+
+    alpha = []
+    for i in range(1, k + 1):
+        j = numpy.arange(1, i)
+        alpha.append(numpy.prod((k - j) / (Ps - j) / i))
+    alpha = numpy.asarray(alpha)
+
+    return hypesub(points.shape[0], points, actDim, bounds, pvec, alpha, k)
+
+
+def hypeIndicatorSampled(points, bounds, k, nrOfSamples):
+
+    nrP = points.shape[0]
+    dim = points.shape[1]
+    F = numpy.zeros(nrP)
+
+    BoxL = numpy.min(points, axis=0)
+
+    alpha = []
+    for i in range(1, k + 1):
+        j = numpy.arange(1, i)
+        alpha.append(numpy.prod((k - j) / (nrP - j) / i))
+    alpha = numpy.asarray(alpha + [0.]*nrP)
+
+    S = numpy.random.uniform(low=BoxL, high=bounds, size=(nrOfSamples, dim))
+
+
+    dominated = numpy.zeros(nrOfSamples, dtype="uint")
+    for j in range(1, nrP+1):
+        B = S - points[j-1]
+        ind = numpy.sum(B >= 0, axis=1) == dim
+        dominated[ind] += 1
+
+    for j in range(1, nrP+1):
+        B = S - points[j-1]
+        ind = numpy.sum(B >= 0, axis=1) == dim
+        x = dominated[ind]
+        F[j-1] = numpy.sum(alpha[x-1])
+
+    F = F * numpy.prod(bounds - BoxL) / nrOfSamples
+
+    return F
+

bluepyopt/deapext/optimisationsCMA.py

@@ -54,6 +54,7 @@ def __init__(self,
                  map_function=None,
                  hof=None,
                  selector_name="single_objective",
+                 weight_hv=0.5,
                  fitness_reduce=numpy.sum):
         """Constructor
 
@@ -89,6 +90,7 @@ def __init__(self,
         self.offspring_size = offspring_size
         self.centroids = centroids
         self.sigma = sigma
+        self.weight_hv = weight_hv
 
         self.selector_name = selector_name
         if self.selector_name == 'single_objective':
@@ -245,6 +247,7 @@ def run(self,
                                       use_scoop=self.use_scoop)
 
             if self.selector_name == 'multi_objective':
+                CMA_es.weight_hv = self.weight_hv
                 to_evaluate = CMA_es.get_parents(self.to_space)
                 fitness = self.toolbox.map(self.toolbox.evaluate, to_evaluate)
                 fitness = list(map(list, fitness))