1from __future__ import division 2import bisect 3from collections import defaultdict, namedtuple 4from itertools import chain 5import math 6from operator import attrgetter, itemgetter 7import random 8 9import numpy 10 11###################################### 12# Non-Dominated Sorting (NSGA-II) # 13###################################### 14 15def selNSGA2(individuals, k, nd='standard'): 16 """Apply NSGA-II selection operator on the *individuals*. Usually, the 17 size of *individuals* will be larger than *k* because any individual 18 present in *individuals* will appear in the returned list at most once. 19 Having the size of *individuals* equals to *k* will have no effect other 20 than sorting the population according to their front rank. The 21 list returned contains references to the input *individuals*. For more 22 details on the NSGA-II operator see [Deb2002]_. 23 24 :param individuals: A list of individuals to select from. 25 :param k: The number of individuals to select. 26 :param nd: Specify the non-dominated algorithm to use: 'standard' or 'log'. 27 :returns: A list of selected individuals. 28 29 .. [Deb2002] Deb, Pratab, Agarwal, and Meyarivan, "A fast elitist 30 non-dominated sorting genetic algorithm for multi-objective 31 optimization: NSGA-II", 2002. 32 """ 33 if nd == 'standard': 34 pareto_fronts = sortNondominated(individuals, k) 35 elif nd == 'log': 36 pareto_fronts = sortLogNondominated(individuals, k) 37 else: 38 raise Exception('selNSGA2: The choice of non-dominated sorting ' 39 'method "{0}" is invalid.'.format(nd)) 40 41 for front in pareto_fronts: 42 assignCrowdingDist(front) 43 44 chosen = list(chain(*pareto_fronts[:-1])) 45 k = k - len(chosen) 46 if k > 0: 47 sorted_front = sorted(pareto_fronts[-1], key=attrgetter("fitness.crowding_dist"), reverse=True) 48 chosen.extend(sorted_front[:k]) 49 50 return chosen 51 52 53def sortNondominated(individuals, k, first_front_only=False): 54 """Sort the first *k* *individuals* into different nondomination levels 55 using the "Fast Nondominated Sorting Approach" proposed by Deb et al., 56 see [Deb2002]_. This algorithm has a time complexity of :math:`O(MN^2)`, 57 where :math:`M` is the number of objectives and :math:`N` the number of 58 individuals. 59 60 :param individuals: A list of individuals to select from. 61 :param k: The number of individuals to select. 62 :param first_front_only: If :obj:`True` sort only the first front and 63 exit. 64 :returns: A list of Pareto fronts (lists), the first list includes 65 nondominated individuals. 66 67 .. [Deb2002] Deb, Pratab, Agarwal, and Meyarivan, "A fast elitist 68 non-dominated sorting genetic algorithm for multi-objective 69 optimization: NSGA-II", 2002. 70 """ 71 if k == 0: 72 return [] 73 74 map_fit_ind = defaultdict(list) 75 for ind in individuals: 76 map_fit_ind[ind.fitness].append(ind) 77 fits = map_fit_ind.keys() 78 79 current_front = [] 80 next_front = [] 81 dominating_fits = defaultdict(int) 82 dominated_fits = defaultdict(list) 83 84 # Rank first Pareto front 85 for i, fit_i in enumerate(fits): 86 for fit_j in fits[i+1:]: 87 if fit_i.dominates(fit_j): 88 dominating_fits[fit_j] += 1 89 dominated_fits[fit_i].append(fit_j) 90 elif fit_j.dominates(fit_i): 91 dominating_fits[fit_i] += 1 92 dominated_fits[fit_j].append(fit_i) 93 if dominating_fits[fit_i] == 0: 94 current_front.append(fit_i) 95 96 fronts = [[]] 97 for fit in current_front: 98 fronts[-1].extend(map_fit_ind[fit]) 99 pareto_sorted = len(fronts[-1]) 100 101 # Rank the next front until all individuals are sorted or 102 # the given number of individual are sorted. 103 if not first_front_only: 104 N = min(len(individuals), k) 105 while pareto_sorted < N: 106 fronts.append([]) 107 for fit_p in current_front: 108 for fit_d in dominated_fits[fit_p]: 109 dominating_fits[fit_d] -= 1 110 if dominating_fits[fit_d] == 0: 111 next_front.append(fit_d) 112 pareto_sorted += len(map_fit_ind[fit_d]) 113 fronts[-1].extend(map_fit_ind[fit_d]) 114 current_front = next_front 115 next_front = [] 116 117 return fronts 118 119def assignCrowdingDist(individuals): 120 """Assign a crowding distance to each individual's fitness. The 121 crowding distance can be retrieve via the :attr:`crowding_dist` 122 attribute of each individual's fitness. 123 """ 124 if len(individuals) == 0: 125 return 126 127 distances = [0.0] * len(individuals) 128 crowd = [(ind.fitness.values, i) for i, ind in enumerate(individuals)] 129 130 nobj = len(individuals[0].fitness.values) 131 132 for i in xrange(nobj): 133 crowd.sort(key=lambda element: element[0][i]) 134 distances[crowd[0][1]] = float("inf") 135 distances[crowd[-1][1]] = float("inf") 136 if crowd[-1][0][i] == crowd[0][0][i]: 137 continue 138 norm = nobj * float(crowd[-1][0][i] - crowd[0][0][i]) 139 for prev, cur, next in zip(crowd[:-2], crowd[1:-1], crowd[2:]): 140 distances[cur[1]] += (next[0][i] - prev[0][i]) / norm 141 142 for i, dist in enumerate(distances): 143 individuals[i].fitness.crowding_dist = dist 144 145def selTournamentDCD(individuals, k): 146 """Tournament selection based on dominance (D) between two individuals, if 147 the two individuals do not interdominate the selection is made 148 based on crowding distance (CD). The *individuals* sequence length has to 149 be a multiple of 4. Starting from the beginning of the selected 150 individuals, two consecutive individuals will be different (assuming all 151 individuals in the input list are unique). Each individual from the input 152 list won't be selected more than twice. 153 154 This selection requires the individuals to have a :attr:`crowding_dist` 155 attribute, which can be set by the :func:`assignCrowdingDist` function. 156 157 :param individuals: A list of individuals to select from. 158 :param k: The number of individuals to select. 159 :returns: A list of selected individuals. 160 """ 161 162 if len(individuals) % 4 != 0: 163 raise ValueError("selTournamentDCD: individuals length must be a multiple of 4") 164 165 if k % 4 != 0: 166 raise ValueError("selTournamentDCD: number of individuals to select must be a multiple of 4") 167 168 def tourn(ind1, ind2): 169 if ind1.fitness.dominates(ind2.fitness): 170 return ind1 171 elif ind2.fitness.dominates(ind1.fitness): 172 return ind2 173 174 if ind1.fitness.crowding_dist < ind2.fitness.crowding_dist: 175 return ind2 176 elif ind1.fitness.crowding_dist > ind2.fitness.crowding_dist: 177 return ind1 178 179 if random.random() <= 0.5: 180 return ind1 181 return ind2 182 183 individuals_1 = random.sample(individuals, len(individuals)) 184 individuals_2 = random.sample(individuals, len(individuals)) 185 186 chosen = [] 187 for i in xrange(0, k, 4): 188 chosen.append(tourn(individuals_1[i], individuals_1[i+1])) 189 chosen.append(tourn(individuals_1[i+2], individuals_1[i+3])) 190 chosen.append(tourn(individuals_2[i], individuals_2[i+1])) 191 chosen.append(tourn(individuals_2[i+2], individuals_2[i+3])) 192 193 return chosen 194 195####################################### 196# Generalized Reduced runtime ND sort # 197####################################### 198 199def identity(obj): 200 """Returns directly the argument *obj*. 201 """ 202 return obj 203 204def isDominated(wvalues1, wvalues2): 205 """Returns whether or not *wvalues1* dominates *wvalues2*. 206 207 :param wvalues1: The weighted fitness values that would be dominated. 208 :param wvalues2: The weighted fitness values of the dominant. 209 :returns: :obj:`True` if wvalues2 dominates wvalues1, :obj:`False` 210 otherwise. 211 """ 212 not_equal = False 213 for self_wvalue, other_wvalue in zip(wvalues1, wvalues2): 214 if self_wvalue > other_wvalue: 215 return False 216 elif self_wvalue < other_wvalue: 217 not_equal = True 218 return not_equal 219 220def median(seq, key=identity): 221 """Returns the median of *seq* - the numeric value separating the higher 222 half of a sample from the lower half. If there is an even number of 223 elements in *seq*, it returns the mean of the two middle values. 224 """ 225 sseq = sorted(seq, key=key) 226 length = len(seq) 227 if length % 2 == 1: 228 return key(sseq[(length - 1) // 2]) 229 else: 230 return (key(sseq[(length - 1) // 2]) + key(sseq[length // 2])) / 2.0 231 232def sortLogNondominated(individuals, k, first_front_only=False): 233 """Sort *individuals* in pareto non-dominated fronts using the Generalized 234 Reduced Run-Time Complexity Non-Dominated Sorting Algorithm presented by 235 Fortin et al. (2013). 236 237 :param individuals: A list of individuals to select from. 238 :returns: A list of Pareto fronts (lists), with the first list being the 239 true Pareto front. 240 """ 241 if k == 0: 242 return [] 243 244 #Separate individuals according to unique fitnesses 245 unique_fits = defaultdict(list) 246 for i, ind in enumerate(individuals): 247 unique_fits[ind.fitness.wvalues].append(ind) 248 249 #Launch the sorting algorithm 250 obj = len(individuals[0].fitness.wvalues)-1 251 fitnesses = unique_fits.keys() 252 front = dict.fromkeys(fitnesses, 0) 253 254 # Sort the fitnesses lexicographically. 255 fitnesses.sort(reverse=True) 256 sortNDHelperA(fitnesses, obj, front) 257 258 #Extract individuals from front list here 259 nbfronts = max(front.values())+1 260 pareto_fronts = [[] for i in range(nbfronts)] 261 for fit in fitnesses: 262 index = front[fit] 263 pareto_fronts[index].extend(unique_fits[fit]) 264 265 # Keep only the fronts required to have k individuals. 266 if not first_front_only: 267 count = 0 268 for i, front in enumerate(pareto_fronts): 269 count += len(front) 270 if count >= k: 271 return pareto_fronts[:i+1] 272 return pareto_fronts 273 else: 274 return pareto_fronts[0] 275 276def sortNDHelperA(fitnesses, obj, front): 277 """Create a non-dominated sorting of S on the first M objectives""" 278 if len(fitnesses) < 2: 279 return 280 elif len(fitnesses) == 2: 281 # Only two individuals, compare them and adjust front number 282 s1, s2 = fitnesses[0], fitnesses[1] 283 if isDominated(s2[:obj+1], s1[:obj+1]): 284 front[s2] = max(front[s2], front[s1] + 1) 285 elif obj == 1: 286 sweepA(fitnesses, front) 287 elif len(frozenset(map(itemgetter(obj), fitnesses))) == 1: 288 #All individuals for objective M are equal: go to objective M-1 289 sortNDHelperA(fitnesses, obj-1, front) 290 else: 291 # More than two individuals, split list and then apply recursion 292 best, worst = splitA(fitnesses, obj) 293 sortNDHelperA(best, obj, front) 294 sortNDHelperB(best, worst, obj-1, front) 295 sortNDHelperA(worst, obj, front) 296 297def splitA(fitnesses, obj): 298 """Partition the set of fitnesses in two according to the median of 299 the objective index *obj*. The values equal to the median are put in 300 the set containing the least elements. 301 """ 302 median_ = median(fitnesses, itemgetter(obj)) 303 best_a, worst_a = [], [] 304 best_b, worst_b = [], [] 305 306 for fit in fitnesses: 307 if fit[obj] > median_: 308 best_a.append(fit) 309 best_b.append(fit) 310 elif fit[obj] < median_: 311 worst_a.append(fit) 312 worst_b.append(fit) 313 else: 314 best_a.append(fit) 315 worst_b.append(fit) 316 317 balance_a = abs(len(best_a) - len(worst_a)) 318 balance_b = abs(len(best_b) - len(worst_b)) 319 320 if balance_a <= balance_b: 321 return best_a, worst_a 322 else: 323 return best_b, worst_b 324 325def sweepA(fitnesses, front): 326 """Update rank number associated to the fitnesses according 327 to the first two objectives using a geometric sweep procedure. 328 """ 329 stairs = [-fitnesses[0][1]] 330 fstairs = [fitnesses[0]] 331 for fit in fitnesses[1:]: 332 idx = bisect.bisect_right(stairs, -fit[1]) 333 if 0 < idx <= len(stairs): 334 fstair = max(fstairs[:idx], key=front.__getitem__) 335 front[fit] = max(front[fit], front[fstair]+1) 336 for i, fstair in enumerate(fstairs[idx:], idx): 337 if front[fstair] == front[fit]: 338 del stairs[i] 339 del fstairs[i] 340 break 341 stairs.insert(idx, -fit[1]) 342 fstairs.insert(idx, fit) 343 344def sortNDHelperB(best, worst, obj, front): 345 """Assign front numbers to the solutions in H according to the solutions 346 in L. The solutions in L are assumed to have correct front numbers and the 347 solutions in H are not compared with each other, as this is supposed to 348 happen after sortNDHelperB is called.""" 349 key = itemgetter(obj) 350 if len(worst) == 0 or len(best) == 0: 351 #One of the lists is empty: nothing to do 352 return 353 elif len(best) == 1 or len(worst) == 1: 354 #One of the lists has one individual: compare directly 355 for hi in worst: 356 for li in best: 357 if isDominated(hi[:obj+1], li[:obj+1]) or hi[:obj+1] == li[:obj+1]: 358 front[hi] = max(front[hi], front[li] + 1) 359 elif obj == 1: 360 sweepB(best, worst, front) 361 elif key(min(best, key=key)) >= key(max(worst, key=key)): 362 #All individuals from L dominate H for objective M: 363 #Also supports the case where every individuals in L and H 364 #has the same value for the current objective 365 #Skip to objective M-1 366 sortNDHelperB(best, worst, obj-1, front) 367 elif key(max(best, key=key)) >= key(min(worst, key=key)): 368 best1, best2, worst1, worst2 = splitB(best, worst, obj) 369 sortNDHelperB(best1, worst1, obj, front) 370 sortNDHelperB(best1, worst2, obj-1, front) 371 sortNDHelperB(best2, worst2, obj, front) 372 373def splitB(best, worst, obj): 374 """Split both best individual and worst sets of fitnesses according 375 to the median of objective *obj* computed on the set containing the 376 most elements. The values equal to the median are attributed so as 377 to balance the four resulting sets as much as possible. 378 """ 379 median_ = median(best if len(best) > len(worst) else worst, itemgetter(obj)) 380 best1_a, best2_a, best1_b, best2_b = [], [], [], [] 381 for fit in best: 382 if fit[obj] > median_: 383 best1_a.append(fit) 384 best1_b.append(fit) 385 elif fit[obj] < median_: 386 best2_a.append(fit) 387 best2_b.append(fit) 388 else: 389 best1_a.append(fit) 390 best2_b.append(fit) 391 392 worst1_a, worst2_a, worst1_b, worst2_b = [], [], [], [] 393 for fit in worst: 394 if fit[obj] > median_: 395 worst1_a.append(fit) 396 worst1_b.append(fit) 397 elif fit[obj] < median_: 398 worst2_a.append(fit) 399 worst2_b.append(fit) 400 else: 401 worst1_a.append(fit) 402 worst2_b.append(fit) 403 404 balance_a = abs(len(best1_a) - len(best2_a) + len(worst1_a) - len(worst2_a)) 405 balance_b = abs(len(best1_b) - len(best2_b) + len(worst1_b) - len(worst2_b)) 406 407 if balance_a <= balance_b: 408 return best1_a, best2_a, worst1_a, worst2_a 409 else: 410 return best1_b, best2_b, worst1_b, worst2_b 411 412def sweepB(best, worst, front): 413 """Adjust the rank number of the worst fitnesses according to 414 the best fitnesses on the first two objectives using a sweep 415 procedure. 416 """ 417 stairs, fstairs = [], [] 418 iter_best = iter(best) 419 next_best = next(iter_best, False) 420 for h in worst: 421 while next_best and h[:2] <= next_best[:2]: 422 insert = True 423 for i, fstair in enumerate(fstairs): 424 if front[fstair] == front[next_best]: 425 if fstair[1] > next_best[1]: 426 insert = False 427 else: 428 del stairs[i], fstairs[i] 429 break 430 if insert: 431 idx = bisect.bisect_right(stairs, -next_best[1]) 432 stairs.insert(idx, -next_best[1]) 433 fstairs.insert(idx, next_best) 434 next_best = next(iter_best, False) 435 436 idx = bisect.bisect_right(stairs, -h[1]) 437 if 0 < idx <= len(stairs): 438 fstair = max(fstairs[:idx], key=front.__getitem__) 439 front[h] = max(front[h], front[fstair]+1) 440 441###################################### 442# Non-Dominated Sorting (NSGA-III) # 443###################################### 444 445NSGA3Memory = namedtuple("NSGA3Memory", ["best_point", "worst_point", "extreme_points"]) 446 447 448class selNSGA3WithMemory(object): 449 """Class version of NSGA-III selection including memory for best, worst and 450 extreme points. Registering this operator in a toolbox is a bit different 451 than classical operators, it requires to instanciate the class instead 452 of just registering the function:: 453 454 >>> from deap import base 455 >>> ref_points = uniform_reference_points(nobj=3, p=12) 456 >>> toolbox = base.Toolbox() 457 >>> toolbox.register("select", selNSGA3WithMemory(ref_points)) 458 459 """ 460 def __init__(self, ref_points, nd="log"): 461 self.ref_points = ref_points 462 self.nd = nd 463 self.best_point = numpy.full((1, ref_points.shape[1]), numpy.inf) 464 self.worst_point = numpy.full((1, ref_points.shape[1]), -numpy.inf) 465 self.extreme_points = None 466 467 def __call__(self, individuals, k): 468 chosen, memory = selNSGA3(individuals, k, self.ref_points, self.nd, 469 self.best_point, self.worst_point, 470 self.extreme_points, True) 471 self.best_point = memory.best_point.reshape((1, -1)) 472 self.worst_point = memory.worst_point.reshape((1, -1)) 473 self.extreme_points = memory.extreme_points 474 return chosen 475 476 477def selNSGA3(individuals, k, ref_points, nd="log", best_point=None, 478 worst_point=None, extreme_points=None, return_memory=False): 479 """Implementation of NSGA-III selection as presented in [Deb2014]_. 480 481 This implementation is partly based on `lmarti/nsgaiii 482 <https://github.com/lmarti/nsgaiii>`_. It departs slightly from the 483 original implementation in that it does not use memory to keep track 484 of ideal and extreme points. This choice has been made to fit the 485 functional api of DEAP. For a version of NSGA-III see 486 :class:`~deap.tools.selNSGA3WithMemory`. 487 488 :param individuals: A list of individuals to select from. 489 :param k: The number of individuals to select. 490 :param ref_points: Reference points to use for niching. 491 :param nd: Specify the non-dominated algorithm to use: 'standard' or 'log'. 492 :param best_point: Best point found at previous generation. If not provided 493 find the best point only from current individuals. 494 :param worst_point: Worst point found at previous generation. If not provided 495 find the worst point only from current individuals. 496 :param extreme_points: Extreme points found at previous generation. If not provided 497 find the extreme points only from current individuals. 498 :param return_memory: If :data:`True`, return the best, worst and extreme points 499 in addition to the chosen individuals. 500 :returns: A list of selected individuals. 501 :returns: If `return_memory` is :data:`True`, a namedtuple with the 502 `best_point`, `worst_point`, and `extreme_points`. 503 504 505 You can generate the reference points using the :func:`uniform_reference_points` 506 function:: 507 508 >>> ref_points = tools.uniform_reference_points(nobj=3, p=12) # doctest: +SKIP 509 >>> selected = selNSGA3(population, k, ref_points) # doctest: +SKIP 510 511 .. [Deb2014] Deb, K., & Jain, H. (2014). An Evolutionary Many-Objective Optimization 512 Algorithm Using Reference-Point-Based Nondominated Sorting Approach, 513 Part I: Solving Problems With Box Constraints. IEEE Transactions on 514 Evolutionary Computation, 18(4), 577-601. doi:10.1109/TEVC.2013.2281535. 515 """ 516 if nd == "standard": 517 pareto_fronts = sortNondominated(individuals, k) 518 elif nd == "log": 519 pareto_fronts = sortLogNondominated(individuals, k) 520 else: 521 raise Exception("selNSGA3: The choice of non-dominated sorting " 522 "method '{0}' is invalid.".format(nd)) 523 524 # Extract fitnesses as a numpy array in the nd-sort order 525 # Use wvalues * -1 to tackle always as a minimization problem 526 fitnesses = numpy.array([ind.fitness.wvalues for f in pareto_fronts for ind in f]) 527 fitnesses *= -1 528 529 # Get best and worst point of population, contrary to pymoo 530 # we don't use memory 531 if best_point is not None and worst_point is not None: 532 best_point = numpy.min(numpy.concatenate((fitnesses, best_point), axis=0), axis=0) 533 worst_point = numpy.max(numpy.concatenate((fitnesses, worst_point), axis=0), axis=0) 534 else: 535 best_point = numpy.min(fitnesses, axis=0) 536 worst_point = numpy.max(fitnesses, axis=0) 537 538 extreme_points = find_extreme_points(fitnesses, best_point, extreme_points) 539 front_worst = numpy.max(fitnesses[:sum(len(f) for f in pareto_fronts), :], axis=0) 540 intercepts = find_intercepts(extreme_points, best_point, worst_point, front_worst) 541 niches, dist = associate_to_niche(fitnesses, ref_points, best_point, intercepts) 542 543 # Get counts per niche for individuals in all front but the last 544 niche_counts = numpy.zeros(len(ref_points), dtype=numpy.int64) 545 index, counts = numpy.unique(niches[:-len(pareto_fronts[-1])], return_counts=True) 546 niche_counts[index] = counts 547 548 # Choose individuals from all fronts but the last 549 chosen = list(chain(*pareto_fronts[:-1])) 550 551 # Use niching to select the remaining individuals 552 sel_count = len(chosen) 553 n = k - sel_count 554 selected = niching(pareto_fronts[-1], n, niches[sel_count:], dist[sel_count:], niche_counts) 555 chosen.extend(selected) 556 557 if return_memory: 558 return chosen, NSGA3Memory(best_point, worst_point, extreme_points) 559 return chosen 560 561 562def find_extreme_points(fitnesses, best_point, extreme_points=None): 563 'Finds the individuals with extreme values for each objective function.' 564 # Keep track of last generation extreme points 565 if extreme_points is not None: 566 fitnesses = numpy.concatenate((fitnesses, extreme_points), axis=0) 567 568 # Translate objectives 569 ft = fitnesses - best_point 570 571 # Find achievement scalarizing function (asf) 572 asf = numpy.eye(best_point.shape[0]) 573 asf[asf == 0] = 1e6 574 asf = numpy.max(ft * asf[:, numpy.newaxis, :], axis=2) 575 576 # Extreme point are the fitnesses with minimal asf 577 min_asf_idx = numpy.argmin(asf, axis=1) 578 return fitnesses[min_asf_idx, :] 579 580 581def find_intercepts(extreme_points, best_point, current_worst, front_worst): 582 """Find intercepts between the hyperplane and each axis with 583 the ideal point as origin.""" 584 # Construct hyperplane sum(f_i^n) = 1 585 b = numpy.ones(extreme_points.shape[1]) 586 A = extreme_points - best_point 587 try: 588 x = numpy.linalg.solve(A, b) 589 except numpy.linalg.LinAlgError: 590 intercepts = current_worst 591 else: 592 intercepts = 1 / x 593 594 if (not numpy.allclose(numpy.dot(A, x), b) or 595 numpy.any(intercepts <= 1e-6) or 596 numpy.any((intercepts + best_point) > current_worst)): 597 intercepts = front_worst 598 599 return intercepts 600 601 602def associate_to_niche(fitnesses, reference_points, best_point, intercepts): 603 """Associates individuals to reference points and calculates niche number. 604 Corresponds to Algorithm 3 of Deb & Jain (2014).""" 605 # Normalize by ideal point and intercepts 606 fn = (fitnesses - best_point) / (intercepts - best_point) 607 608 # Create distance matrix 609 fn = numpy.repeat(numpy.expand_dims(fn, axis=1), len(reference_points), axis=1) 610 norm = numpy.linalg.norm(reference_points, axis=1) 611 612 distances = numpy.sum(fn * reference_points, axis=2) / norm.reshape(1, -1) 613 distances = distances[:, :, numpy.newaxis] * reference_points[numpy.newaxis, :, :] / norm[numpy.newaxis, :, numpy.newaxis] 614 distances = numpy.linalg.norm(distances - fn, axis=2) 615 616 # Retrieve min distance niche index 617 niches = numpy.argmin(distances, axis=1) 618 distances = distances[range(niches.shape[0]), niches] 619 return niches, distances 620 621 622def niching(individuals, k, niches, distances, niche_counts): 623 selected = [] 624 available = numpy.ones(len(individuals), dtype=numpy.bool) 625 while len(selected) < k: 626 # Maximum number of individuals (niches) to select in that round 627 n = k - len(selected) 628 629 # Find the available niches and the minimum niche count in them 630 available_niches = numpy.zeros(len(niche_counts), dtype=numpy.bool) 631 available_niches[numpy.unique(niches[available])] = True 632 min_count = numpy.min(niche_counts[available_niches]) 633 634 # Select at most n niches with the minimum count 635 selected_niches = numpy.flatnonzero(numpy.logical_and(available_niches, niche_counts == min_count)) 636 numpy.random.shuffle(selected_niches) 637 selected_niches = selected_niches[:n] 638 639 for niche in selected_niches: 640 # Find the individuals associated with this niche 641 niche_individuals = numpy.flatnonzero(niches == niche) 642 numpy.random.shuffle(niche_individuals) 643 644 # If no individual in that niche, select the closest to reference 645 # Else select randomly 646 if niche_counts[niche] == 0: 647 sel_index = niche_individuals[numpy.argmin(distances[niche_individuals])] 648 else: 649 sel_index = niche_individuals[0] 650 651 # Update availability, counts and selection 652 available[sel_index] = False 653 niche_counts[niche] += 1 654 selected.append(individuals[sel_index]) 655 656 return selected 657 658 659def uniform_reference_points(nobj, p=4, scaling=None): 660 """Generate reference points uniformly on the hyperplane intersecting 661 each axis at 1. The scaling factor is used to combine multiple layers of 662 reference points. 663 """ 664 def gen_refs_recursive(ref, nobj, left, total, depth): 665 points = [] 666 if depth == nobj - 1: 667 ref[depth] = left / total 668 points.append(ref) 669 else: 670 for i in range(left + 1): 671 ref[depth] = i / total 672 points.extend(gen_refs_recursive(ref.copy(), nobj, left - i, total, depth + 1)) 673 return points 674 675 ref_points = numpy.array(gen_refs_recursive(numpy.zeros(nobj), nobj, p, p, 0)) 676 if scaling is not None: 677 ref_points *= scaling 678 ref_points += (1 - scaling) / nobj 679 680 return ref_points 681 682 683###################################### 684# Strength Pareto (SPEA-II) # 685###################################### 686 687def selSPEA2(individuals, k): 688 """Apply SPEA-II selection operator on the *individuals*. Usually, the 689 size of *individuals* will be larger than *n* because any individual 690 present in *individuals* will appear in the returned list at most once. 691 Having the size of *individuals* equals to *n* will have no effect other 692 than sorting the population according to a strength Pareto scheme. The 693 list returned contains references to the input *individuals*. For more 694 details on the SPEA-II operator see [Zitzler2001]_. 695 696 :param individuals: A list of individuals to select from. 697 :param k: The number of individuals to select. 698 :returns: A list of selected individuals. 699 700 .. [Zitzler2001] Zitzler, Laumanns and Thiele, "SPEA 2: Improving the 701 strength Pareto evolutionary algorithm", 2001. 702 """ 703 N = len(individuals) 704 L = len(individuals[0].fitness.values) 705 K = math.sqrt(N) 706 strength_fits = [0] * N 707 fits = [0] * N 708 dominating_inds = [list() for i in xrange(N)] 709 710 for i, ind_i in enumerate(individuals): 711 for j, ind_j in enumerate(individuals[i+1:], i+1): 712 if ind_i.fitness.dominates(ind_j.fitness): 713 strength_fits[i] += 1 714 dominating_inds[j].append(i) 715 elif ind_j.fitness.dominates(ind_i.fitness): 716 strength_fits[j] += 1 717 dominating_inds[i].append(j) 718 719 for i in xrange(N): 720 for j in dominating_inds[i]: 721 fits[i] += strength_fits[j] 722 723 # Choose all non-dominated individuals 724 chosen_indices = [i for i in xrange(N) if fits[i] < 1] 725 726 if len(chosen_indices) < k: # The archive is too small 727 for i in xrange(N): 728 distances = [0.0] * N 729 for j in xrange(i + 1, N): 730 dist = 0.0 731 for l in xrange(L): 732 val = individuals[i].fitness.values[l] - \ 733 individuals[j].fitness.values[l] 734 dist += val * val 735 distances[j] = dist 736 kth_dist = _randomizedSelect(distances, 0, N - 1, K) 737 density = 1.0 / (kth_dist + 2.0) 738 fits[i] += density 739 740 next_indices = [(fits[i], i) for i in xrange(N) 741 if not i in chosen_indices] 742 next_indices.sort() 743 #print next_indices 744 chosen_indices += [i for _, i in next_indices[:k - len(chosen_indices)]] 745 746 elif len(chosen_indices) > k: # The archive is too large 747 N = len(chosen_indices) 748 distances = [[0.0] * N for i in xrange(N)] 749 sorted_indices = [[0] * N for i in xrange(N)] 750 for i in xrange(N): 751 for j in xrange(i + 1, N): 752 dist = 0.0 753 for l in xrange(L): 754 val = individuals[chosen_indices[i]].fitness.values[l] - \ 755 individuals[chosen_indices[j]].fitness.values[l] 756 dist += val * val 757 distances[i][j] = dist 758 distances[j][i] = dist 759 distances[i][i] = -1 760 761 # Insert sort is faster than quick sort for short arrays 762 for i in xrange(N): 763 for j in xrange(1, N): 764 l = j 765 while l > 0 and distances[i][j] < distances[i][sorted_indices[i][l - 1]]: 766 sorted_indices[i][l] = sorted_indices[i][l - 1] 767 l -= 1 768 sorted_indices[i][l] = j 769 770 size = N 771 to_remove = [] 772 while size > k: 773 # Search for minimal distance 774 min_pos = 0 775 for i in xrange(1, N): 776 for j in xrange(1, size): 777 dist_i_sorted_j = distances[i][sorted_indices[i][j]] 778 dist_min_sorted_j = distances[min_pos][sorted_indices[min_pos][j]] 779 780 if dist_i_sorted_j < dist_min_sorted_j: 781 min_pos = i 782 break 783 elif dist_i_sorted_j > dist_min_sorted_j: 784 break 785 786 # Remove minimal distance from sorted_indices 787 for i in xrange(N): 788 distances[i][min_pos] = float("inf") 789 distances[min_pos][i] = float("inf") 790 791 for j in xrange(1, size - 1): 792 if sorted_indices[i][j] == min_pos: 793 sorted_indices[i][j] = sorted_indices[i][j + 1] 794 sorted_indices[i][j + 1] = min_pos 795 796 # Remove corresponding individual from chosen_indices 797 to_remove.append(min_pos) 798 size -= 1 799 800 for index in reversed(sorted(to_remove)): 801 del chosen_indices[index] 802 803 return [individuals[i] for i in chosen_indices] 804 805def _randomizedSelect(array, begin, end, i): 806 """Allows to select the ith smallest element from array without sorting it. 807 Runtime is expected to be O(n). 808 """ 809 if begin == end: 810 return array[begin] 811 q = _randomizedPartition(array, begin, end) 812 k = q - begin + 1 813 if i < k: 814 return _randomizedSelect(array, begin, q, i) 815 else: 816 return _randomizedSelect(array, q + 1, end, i - k) 817 818def _randomizedPartition(array, begin, end): 819 i = random.randint(begin, end) 820 array[begin], array[i] = array[i], array[begin] 821 return _partition(array, begin, end) 822 823def _partition(array, begin, end): 824 x = array[begin] 825 i = begin - 1 826 j = end + 1 827 while True: 828 j -= 1 829 while array[j] > x: 830 j -= 1 831 i += 1 832 while array[i] < x: 833 i += 1 834 if i < j: 835 array[i], array[j] = array[j], array[i] 836 else: 837 return j 838 839 840__all__ = ['selNSGA2', 'selNSGA3', 'selNSGA3WithMemory', 'selSPEA2', 'sortNondominated', 'sortLogNondominated', 841 'selTournamentDCD', 'uniform_reference_points'] 842