Multi-Objective Optimisation - Part 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Multi Objective Evolutionary Algorithms (Part 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In this chapter we continue with multi-objective search algorithms, and consider some alternatives to NSGA-II. First we need to import the usual stuff, set up our wrapper class and plotting infrastructure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "## Setting the Scene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from IPython.display import HTML\n",
    "from math import sqrt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As an example problem instance, we consider the Next Release Problem (NRP) again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "profit_map = {}        # customer id => weight\n",
    "requirements_map = {}  # customer id => [ requirement_id * ]\n",
    "cost_map = {}          # requirement id => cost\n",
    "dependency_map = {}    # requirement id => [ requirement_id * ]\n",
    "num_requirements = 140"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We'll need our parser for the standard format again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def parse_nrp(filename):\n",
    "\n",
    "    with open(filename) as fp:\n",
    "        levels   = int(fp.readline())\n",
    "        req_id = 1\n",
    "        for level in range(levels):\n",
    "            num_reqs = int(fp.readline())\n",
    "            for cost in fp.readline().split():\n",
    "                cost_map[req_id] = int(cost)\n",
    "                req_id += 1\n",
    "        total_deps = int(fp.readline())\n",
    "        for num_dep in range(total_deps):\n",
    "            r1, r2 = fp.readline().split()\n",
    "            r1_id = int(r1)\n",
    "            if r1_id in dependency_map:\n",
    "                dependency_map[int(r1)].append(int(r2))\n",
    "            else:\n",
    "                dependency_map[int(r1)] = [int(r2)]\n",
    "\n",
    "        total_customers = int(fp.readline())\n",
    "        for num_cust in range(total_customers):\n",
    "            customer = fp.readline().split()\n",
    "            profit_map[num_cust + 1] = int(customer[0])\n",
    "            num_reqs = int(customer[1])\n",
    "            requirements = []\n",
    "            for num_req in range(num_reqs):\n",
    "                requirements.append(int(customer[2+num_req]))\n",
    "            requirements_map[num_cust + 1] = requirements\n",
    "\n",
    "        num_requirements = len(cost_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As running example, we will use the small instance `nrp1.txt`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "parse_nrp(\"data/nrp/nrp1.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our first objective is to maximise the profit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Profit is the sum of the weights of customers whose requirements are satisfied\n",
    "def function1(solution):\n",
    "    fitness = 0\n",
    "\n",
    "    requirements = set([x+1 for x in range(len(solution)) if solution[x] == 1])\n",
    "    for customer_id in profit_map.keys():\n",
    "        reqs = set(requirements_map[customer_id])\n",
    "        if reqs.issubset(requirements):\n",
    "            fitness += profit_map[customer_id]\n",
    "\n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The second objective is to minimise the costs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Cost is the sum of costs of the implemented requirements\n",
    "def function2(solution):\n",
    "    cost = 0\n",
    "\n",
    "    for i in range(num_requirements):\n",
    "        if solution[i]:\n",
    "            cost += cost_map[i+1]\n",
    "\n",
    "    return cost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "total_costs = sum([cost_map[i+1] for i in range(num_requirements)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "total_profit = sum(profit_map.values())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To avoid redundantly calculating fitness values, we will use the wrapper class for the list to cache values in attributes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "class L(list):\n",
    "    \"\"\"\n",
    "    A subclass of list that can accept additional attributes.\n",
    "    Should be able to be used just like a regular list.\n",
    "    \"\"\"\n",
    "    def __new__(self, *args, **kwargs):\n",
    "        return super(L, self).__new__(self, args, kwargs)\n",
    "\n",
    "    def __init__(self, *args, **kwargs):\n",
    "        if len(args) == 1 and hasattr(args[0], '__iter__'):\n",
    "            list.__init__(self, args[0])\n",
    "        else:\n",
    "            list.__init__(self, args)\n",
    "        self.__dict__.update(kwargs)\n",
    "\n",
    "    def __call__(self, **kwargs):\n",
    "        self.__dict__.update(kwargs)\n",
    "        return self"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can now store the fitness values of individuals as attributes of the objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def evaluate(individual):\n",
    "    individual.fitness1 = function1(individual)\n",
    "    individual.fitness2 = function2(individual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Individuals of NRP solutions are instances of `L` rather than lists, and we can define our usual search operators for bitvector representation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_individual():\n",
    "    individual = L(random.choice([0,1]) for _ in range(num_requirements))\n",
    "    evaluate(individual)\n",
    "    return individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def mutate(solution):\n",
    "    P_mutate = 1/len(solution)\n",
    "    mutated = L(solution[:])\n",
    "    for position in range(len(solution)):\n",
    "        if random.random() < P_mutate:\n",
    "            mutated[position] = 1 - mutated[position]\n",
    "    evaluate(mutated)\n",
    "    return mutated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def crossover(parent1, parent2):\n",
    "    pos = random.randint(0, len(parent1))\n",
    "    offspring1 = L(parent1[:pos] + parent2[pos:])\n",
    "    offspring2 = L(parent2[:pos] + parent1[pos:])\n",
    "    return (offspring1, offspring2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In addition to the representation, we also want to keep using animations to study how the algorithms explore the search space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "ims = []  # global variable to store images of the animation\n",
    "    \n",
    "def initialise_plot():\n",
    "    global ims\n",
    "    global fig\n",
    "    global ax\n",
    "    \n",
    "    ims = []\n",
    "\n",
    "    %matplotlib agg\n",
    "    fig, ax = plt.subplots()\n",
    "    plt.xlabel('Profit', fontsize=15)\n",
    "    plt.ylabel('Cost', fontsize=15)\n",
    "    ims = []\n",
    "    %matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For each iteration of the algorithm, we will update this animation with a snapshot of the current population and their fitness values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def plot(population):\n",
    "    function1_values = [x.fitness1 for x in population]\n",
    "    function2_values = [x.fitness2 for x in population]\n",
    "    \n",
    "    ims.append((ax.scatter(function1_values, function2_values, color=\"blue\"),))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Baseline 1: Random Search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When studying different search algorithms, random search always servers as a sanity check to compare against. To generalise random search to multi-objective random search, we'll need to define the dominance relation again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def dominates(solution1, solution2):\n",
    "    \"\"\"\n",
    "     A solution x(1) is said to dominate the other solution x(2) if both condition 1 and 2 below are true:\n",
    "\n",
    "     Condition 1: x(1) is no worse than x(2) for all objectives\n",
    "     Condition 2: x(1) is strictly better than x(2) in at least one objective\n",
    "     \n",
    "     We are maximising fitness 1, but minimising fitness 2\n",
    "    \"\"\"\n",
    "    if solution1.fitness1 < solution2.fitness1 or solution1.fitness2 > solution2.fitness2:\n",
    "        return False\n",
    "\n",
    "    if solution1.fitness1 > solution2.fitness1 or solution1.fitness2 < solution2.fitness2:\n",
    "        return True\n",
    "\n",
    "    return False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "An operation we will frequently need in this chapter, and also for random search, is to extract the non-dominated solutions for a given collection of solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_nondominated(population):\n",
    "    nondominated = []\n",
    "\n",
    "    for x in population:\n",
    "        dominated = False\n",
    "        for y in population:\n",
    "            if dominates(y, x):\n",
    "                dominated = True\n",
    "                break\n",
    "        if not dominated:\n",
    "            nondominated.append(x)\n",
    "\n",
    "    return nondominated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given these ingredients, we can generalise random search to multi-objective random search by repeatedly sampling random individuals, and keeping those that are not dominated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_gen = 100\n",
    "population_size = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def random_moo():\n",
    "    initialise_plot()\n",
    "    result = []\n",
    "    for iteration in range(max_gen * population_size):\n",
    "        candidate = get_random_individual()\n",
    "        result = get_nondominated(result + [candidate])\n",
    "        if iteration % population_size == 0:\n",
    "            plot(result)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Similar to the last chapter, we can look at the evolution for each algorithm by plotting the solutions found for each iteration (or in this case, after each `population_size` individuals have been evaluated, to speed up the animation)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = random_moo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To systematically compare sarch algorithms, we need to run experiments with multiple repetitions again. For this we use our usual helper function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.utils import io\n",
    "\n",
    "def run_times(algorithm, repetitions):\n",
    "    global ims\n",
    "    result = []\n",
    "    for i in range(repetitions):\n",
    "        ims = []\n",
    "        with io.capture_output() as captured: \n",
    "            front = algorithm()\n",
    "        result.append(front)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We also need a metric to compare algorithms with. In the last chapter we considered several different metrics, in particular the hypervolume."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def hypervolume(front, r):\n",
    "    front.sort(key=lambda i: i.fitness1)\n",
    "\n",
    "    hv = (abs(r[0] - front[0].fitness1) / total_profit) * (abs(r[1] - front[0].fitness2) / total_costs)\n",
    "    \n",
    "    for i in range(1, len(front)):\n",
    "        hv += (abs(front[i-1].fitness1 - front[i].fitness1) / total_profit) * (abs(r[1] - front[i].fitness2) / total_costs)\n",
    "    \n",
    "    return hv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will keep all experiment results in a dictionary `results` such that we can produce comparative boxplots throughout this chapter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fronts_random = run_times(random_moo, 10)\n",
    "results = {}\n",
    "results[\"Random\"] = [hypervolume(front, (0,total_costs)) for front in fronts_random]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Baseline 2: NSGA-II"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the last chapter we introduced NSGA-II as our first multi-objective search algorithm, and indeed it is one of the most popular multi-objective algorithms and has been shown to be effective and efficient for many different problems. Hence we definitely need to include it in our comparison in this chapter. In the following, the individual bits and pieces of NSGA-II from the previous chapter are re-introduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def fast_non_dominated_sort(solutions):\n",
    "    front = [[]]\n",
    "\n",
    "    S = [[] for _ in range(len(solutions))]\n",
    "    n = [0 for _ in range(len(solutions))]\n",
    "\n",
    "    for p in range(len(solutions)):\n",
    "        S[p] = []\n",
    "        n[p] = 0\n",
    "        for q in range(len(solutions)):\n",
    "            if dominates(solutions[p], solutions[q]):\n",
    "                S[p].append(q)\n",
    "            elif dominates(solutions[q], solutions[p]):\n",
    "                n[p] = n[p] + 1\n",
    "\n",
    "        if n[p] == 0:\n",
    "            front[0].append(p)\n",
    "            solutions[p].rank = 0\n",
    "\n",
    "    i = 0\n",
    "    while front[i]:\n",
    "        Q = []\n",
    "        for p in front[i]:\n",
    "            for q in S[p]:\n",
    "                n[q] = n[q] - 1\n",
    "                if n[q] == 0:\n",
    "                    Q.append(q)\n",
    "                    solutions[q].rank = i + 1\n",
    "        i = i + 1\n",
    "        front.append(Q)\n",
    "\n",
    "    del front[-1]\n",
    "    return front"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def calculate_crowding_distance_and_sort(front):\n",
    "\n",
    "    data = [(x, front[x].fitness1, front[x].fitness2) for x in range(len(front))]\n",
    "    sorted1 = [(x, y) for (x, y, z) in sorted(data, key=lambda tup: tup[1])]\n",
    "    sorted2 = [(x, z) for (x, y, z) in sorted(data, key=lambda tup: tup[2])]\n",
    "\n",
    "    distance = [0 for _ in range(0,len(front))]\n",
    "    range_fitness1 = max(x.fitness1 for x in front) - min(x.fitness1 for x in front)\n",
    "    range_fitness2 = max(x.fitness2 for x in front) - min(x.fitness2 for x in front)\n",
    "\n",
    "    distance[sorted1[0][0]] = sys.maxsize\n",
    "    distance[sorted1[-1][0]] = sys.maxsize\n",
    "\n",
    "    distance[sorted2[0][0]] = sys.maxsize\n",
    "    distance[sorted2[-1][0]] = sys.maxsize\n",
    "\n",
    "    for k in range(1,len(front)-1):\n",
    "        index = sorted1[k][0]\n",
    "        distance[index] = distance[index] + (sorted1[k+1][1] - sorted1[k-1][1]) / range_fitness1\n",
    "\n",
    "    for k in range(1,len(front)-1):\n",
    "        index = sorted2[k][0]\n",
    "        distance[index] = distance[index] + (sorted2[k+1][1] - sorted2[k-1][1]) / range_fitness2\n",
    "\n",
    "    for k in range(0, len(front)):\n",
    "        front[k].distance = distance[k]\n",
    "\n",
    "    front.sort(key = lambda i: i.distance, reverse=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def binary_rank_tournament(population):\n",
    "    individual1 = random.choice(population)\n",
    "    individual2 = random.choice(population)\n",
    "\n",
    "    if individual1.rank < individual2.rank:\n",
    "        return individual1\n",
    "    elif individual1.rank > individual2.rank:\n",
    "        return individual2\n",
    "    else:\n",
    "        return max([individual1, individual2], key = lambda i: i.distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def generate_offspring(population):\n",
    "    offspring_population = []\n",
    "    while len(offspring_population) < len(population):\n",
    "        parent1 = binary_rank_tournament(population)\n",
    "        parent2 = binary_rank_tournament(population)\n",
    "        if random.random() < P_xover:\n",
    "            offspring1, offspring2 = crossover(parent1, parent2)\n",
    "        else:\n",
    "            offspring1, offspring2 = parent1, parent2\n",
    "            \n",
    "        offspring1 = mutate(offspring1)\n",
    "        offspring2 = mutate(offspring2)\n",
    "\n",
    "        offspring_population.append(offspring1)\n",
    "        offspring_population.append(offspring2)\n",
    "\n",
    "    return offspring_population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_initial_population(population_size):\n",
    "    population = [get_random_individual() for _ in range(population_size)]\n",
    "    fronts = fast_non_dominated_sort(population)\n",
    "    \n",
    "    for front_indices in fronts:\n",
    "        front = [population[index] for index in front_indices]\n",
    "        calculate_crowding_distance_and_sort(front)\n",
    "\n",
    "    return population"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def nsgaii():\n",
    "    initialise_plot()\n",
    "    population = get_initial_population(population_size)\n",
    "    offspring_population = generate_offspring(population)\n",
    "\n",
    "    for iteration in range(max_gen):\n",
    "        combined = population + offspring_population\n",
    "        #plot(combined)\n",
    "        fronts = fast_non_dominated_sort(combined)\n",
    "        population = []\n",
    "\n",
    "        for front_indices in fronts:\n",
    "            front = [combined[index] for index in front_indices]\n",
    "            calculate_crowding_distance_and_sort(front)\n",
    "\n",
    "            for i in front:\n",
    "                population.append(i)\n",
    "                if len(population) == population_size:\n",
    "                    break\n",
    "            if len(population) == population_size:\n",
    "                break\n",
    "\n",
    "        plot(population)\n",
    "        offspring_population = generate_offspring(population)\n",
    "\n",
    "    non_dominated_sorted_solution = fast_non_dominated_sort(population)\n",
    "    result = [population[x] for x in non_dominated_sorted_solution[0]]\n",
    "    plot(result)\n",
    "\n",
    "    return result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "First let's consider the evolution of an individual run again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "P_xover = 0.7\n",
    "result = nsgaii()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We also need to collect some more data for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fronts_nsgaii = run_times(nsgaii, 10)\n",
    "results[\"NSGA-II\"] = [hypervolume(front, (0,total_costs)) for front in fronts_nsgaii]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our first sanity check is whether NSGA-II is indeed better than random search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## PAES"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As first alternative example of a multi-objective search algorithm, we consider the Pareto Archived Evolution Strategy (PAES):\n",
    "\n",
    "J. Knowles, D. Corne. “The pareto archived evolution strategy: A new baseline algorithm for pareto multiobjective optimisation.” In Congress on Evolutionary Computation (CEC99) (Vol. 1, pp. 98-105). 1999\n",
    "\n",
    "PAES is intended to be the simplest possible non-trivial algorithm capable of generating diverse solutions in the Pareto optimal set. It essentiall is a (1+1)ES adapted for multi-objective search, and is intended as baseline approach for evaluation of other algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Like a (1+1)ES the algorithm has a population of size 1 and produces 1 offspring. In addition, there is an archive which stores non-dominated solutions (subject to diversity criteria). At the end of the run, the archive is the final result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "1. Generate random current solution\n",
    "2. Evaluate and add to archive\n",
    "3. While not done:\n",
    "4. Create candidate by mutating current solution\n",
    "5. If candidate is dominated by current, reject it\n",
    "6. Else if current is dominated by the candidate, accept candidate and add it to the archive\n",
    "7. Else Compare candidate solutions with archive members, update archive and current"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The archive stores only non-dominated solutions, and has a maximum size. The archive is also used in order to decide whether or not to accept an offspring:\n",
    "- Candidates which dominate the archive are always accepted and archived.\n",
    "- Candidates which are dominated by the archive are always rejected.\n",
    "- Non-dominated are accepted/archived based on how many similar individuals already exist\n",
    "\n",
    "When an individual is added to the archive, we need to remove all individuals from the archive that are dominated by this new member. However, when the maximum archive size is reached, we may have to decide which individuals to keep in the archive. In this case, PAES checks if the new candidate would increase the diversity in the archive. If so, it replaces the individual risiding in the most crowded part of the archive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We start by defining some helper functions since PAES works with an archive. First, we add a function that tells us if a candidate individual is dominated by any of the individuals in the archive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def is_dominated(candidate, archive):\n",
    "    for i in archive:\n",
    "        if dominates(i, candidate):\n",
    "            return True\n",
    "    return False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We also need an update operation that removes all dominated individuals from an archive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def remove_dominated(candidate, archive):\n",
    "    copy = [x for x in archive if not dominates(candidate, x)]\n",
    "    archive.clear()\n",
    "    archive.extend(copy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Third,we need to select all solutions in the archive that are dominated by an individual, as we want to remove all of them when adding a new non-dominated solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_dominated(candidate, archive):\n",
    "    return [x for x in archive if dominates(candidate, x)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Finally, we need an operation that updates the archive. As we are using lists but technically the archive should be an set for PAES, we'll have to add a check here (inefficient, but shorter than redefining previous operations to work on sets rather than lists):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_archive(individual, archive):\n",
    "    if not individual in archive:\n",
    "        archive.append(individual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When there is dominance, the choice of what to do is easy. The more tricky case is when the new individual is not dominated by the archive. As the archive has a maximum size, we need to decide which individuals to keep once we hit the maximum. In PAES, this is done using a grid from which we can infer which areas of the objective space are more crowded. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "grid_size = 10\n",
    "\n",
    "def get_grid(archive):\n",
    "    grid = {}\n",
    "    max_profit = max([i.fitness1 for i in archive])+1\n",
    "    max_effort = max([i.fitness2 for i in archive])+1\n",
    "    f1step = max_profit / grid_size\n",
    "    f2step = max_effort / grid_size\n",
    "\n",
    "    for x in range(grid_size):\n",
    "        grid[x] = {}\n",
    "        for y in range(grid_size):\n",
    "            grid[x][y] = []\n",
    "\n",
    "    for individual in archive:\n",
    "        x = individual.fitness1 // f1step\n",
    "        y = individual.fitness2 // f2step\n",
    "        grid[x][y].append(individual)\n",
    "\n",
    "    return grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We want to add individuals to the grid if they help us preserve diversity. Thus, we check if we would add the individual to the most crowded grid cell. If we are, then we are not increasing diversity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def increases_diversity(candidate, grid, archive):\n",
    "    max_crowd = 0\n",
    "    target_crowd = 0\n",
    "\n",
    "    # We have to redundantly calculate this to figure out the step size...\n",
    "    max_profit = max([i.fitness1 for i in archive]) + 1\n",
    "    max_cost = max([i.fitness2 for i in archive]) + 1\n",
    "    f1step = max_profit / grid_size\n",
    "    f2step = max_cost / grid_size\n",
    "\n",
    "    target_x = candidate.fitness1 // f1step\n",
    "    target_y = candidate.fitness2 // f2step\n",
    "    \n",
    "    for x in range(grid_size):\n",
    "        for y in range(grid_size):\n",
    "            num = len(grid[x][y])\n",
    "            if num > max_crowd:\n",
    "                max_crowd = num\n",
    "            if target_x == x and target_y == y:\n",
    "                target_crowd = num\n",
    "\n",
    "    return target_crowd < max_crowd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can also compare individuals in terms of how crowded their grid cells are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def less_crowded_than(individual1, individual2, grid):\n",
    "    individual1_crowd = 0\n",
    "    individual2_crowd = 0\n",
    "\n",
    "    for x in range(grid_size):\n",
    "        for y in range(grid_size):\n",
    "            if individual1 in grid[x][y]:\n",
    "                individual1_crowd = len(grid[x][y])\n",
    "            if individual2 in grid[x][y]:\n",
    "                individual2_crowd = len(grid[x][y])\n",
    "\n",
    "    return individual1_crowd < individual2_crowd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If we need to reduce the size of the archive, we randomly pick one individual from the most crowded grid cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def remove_from_grid(grid, archive):\n",
    "    pos_x, pos_y = 0, 0\n",
    "    max_crowd = 0\n",
    "\n",
    "    for x in range(grid_size):\n",
    "        for y in range(grid_size):\n",
    "            num = len(grid[x][y])\n",
    "            if num > max_crowd:\n",
    "                max_crowd = num\n",
    "                pos_x, pos_y = x, y\n",
    "\n",
    "    selected = random.choice(grid[pos_x][pos_y])\n",
    "    archive.remove(selected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The PAES algorithm itself mainly consists of the logic to decide what to do in case of non-domination:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def paes():\n",
    "    initialise_plot()\n",
    "    archive = []\n",
    "    current = get_random_individual()\n",
    "    add_to_archive(current, archive)\n",
    "\n",
    "    for step in range(population_size * max_gen):\n",
    "        candidate = mutate(current)\n",
    "        while candidate == current or candidate in archive:\n",
    "            candidate = mutate(current)\n",
    "        if step % population_size == 0:\n",
    "            plot(archive + [candidate])\n",
    "\n",
    "        if not dominates(current, candidate) and not is_dominated(candidate, archive):\n",
    "            dominated_archive = get_dominated(candidate, archive)\n",
    "            if dominated_archive:\n",
    "                # If the candidate dominates something in the archive\n",
    "                # we keep it in the archive and make it the new current\n",
    "                remove_dominated(candidate, archive)\n",
    "                add_to_archive(candidate, archive)\n",
    "                current = candidate\n",
    "            else:\n",
    "                if len(archive) == population_size:\n",
    "                    # Maximum archive size reached\n",
    "                    grid = get_grid(archive)\n",
    "                    if increases_diversity(candidate, grid, archive):\n",
    "                        remove_from_grid(grid, archive)\n",
    "                        add_to_archive(candidate, archive)\n",
    "                        grid = get_grid(archive)\n",
    "                        if less_crowded_than(candidate, current, grid):\n",
    "                            current = candidate\n",
    "                else:\n",
    "                    # Enough space in archive\n",
    "                    add_to_archive(candidate, archive)\n",
    "                    grid = get_grid(archive)\n",
    "                    if less_crowded_than(candidate, current, grid):\n",
    "                        current = candidate\n",
    "\n",
    "    return get_nondominated(archive)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = paes()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As always, we add some data and compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fronts_paes = run_times(paes, 10)\n",
    "results[\"PAES\"] = [hypervolume(front, (0,total_costs)) for front in fronts_paes]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## SPEA-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The Strength Pareto Evolutionary Algorithm (SPEA) version 2 is often treated as the main competitor of NSGA-II. It was defined in the following paper:\n",
    "\n",
    "E. Zitzler, M. Laumanns and L. Thiele, \"SPEA2: Improving the strength pareto evolutionary algorithm for multiobjective optimization\", in Evolutionary Methods for Design, Optimisation and Control with Application to Industrial Problems, 2002.\n",
    "\n",
    "This algorithm assigns a fitness to an individual based on the strength of its dominators. This fitness value is then minimised to produce solutions that are not dominated.\n",
    "\n",
    "Another central distinguishing feature is the use of an archive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The overall workflow of the algorithm is as follows:\n",
    "1. Create initial population $P_0$ and create an empty archive $\\bar{P}_0$.\n",
    "2. Evaluate fitness values of the individuals in $P_t$ and $\\bar{P}_t$\n",
    "3. Copy all non-dominated individuals in $P_t$ and $\\bar{P}_t$ to $\\bar{P}_{t+1}$. If $\\bar{P}_{t+1}$ exceeds the maximum archive size then truncate; if it is smaller then fill up with dominated individuals from $P_t$ and $\\bar{P}_t$.\n",
    "4. Perform binary tournament selection with replacement on $\\bar{P}_{t+1}$ to fill the mating pool\n",
    "5. Apply recombination and mutation operators to the mating pool and set $P_{t+1}$ to the resulting population.\n",
    "6. Repeat from step 2 until done\n",
    "\n",
    "The solution is represented by the non-dominated individuals in $\\bar{P}$ at the end."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's first look at how the fitness values are determined. The score of an individual is the number of solutions that it dominates:\n",
    "\n",
    "$Score(i) = |j \\; | \\; j \\in P_t \\cup \\bar{P}_t \\wedge i \t\\succ j   |$\n",
    "\n",
    "Here, $i \\succ j$ denotes that $i$ dominates $j$. (Note that Zitzler et al. use the operator $\\succ$ the other way round as Deb et al.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def score(individual, combined):\n",
    "    score = 0\n",
    "    for other in combined:\n",
    "        if dominates(individual, other):\n",
    "            score += 1\n",
    "\n",
    "    individual.score = score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The raw fitness value of an individual is calculated as the sum of its dominators' strengths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def raw_fitness(individual, combined):\n",
    "    raw_fitness = 0\n",
    "    for other in combined:\n",
    "        if dominates(other, individual):\n",
    "            raw_fitness += other.score\n",
    "\n",
    "    individual.raw_fitness = raw_fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The fitness value considers not only the raw fitness, but also the density, such that individuals in less populated areas of the search space are preferred:\n",
    "\n",
    "$Density(i) = \\frac{1}{\\sigma_i^k + 2}$\n",
    "\n",
    "Here, $\\sigma_i^k$ denotes the distance of $i$ to its $k$th neighbour, and $k = \\sqrt{N + \\bar{N}}$, where $N$ is the size of the population and $\\bar{N}$ is the size of the archive."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The distance $\\sigma$ is defined in terms of the objective values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_distance(individual, other):\n",
    "    sum = 0\n",
    "\n",
    "    sum += (individual.fitness1 - other.fitness1) ** 2\n",
    "    sum += (individual.fitness2 - other.fitness2) ** 2\n",
    "\n",
    "    return sqrt(sum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we just need to calculate these distances, sort, and assign to the individuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def density(individual, combined):\n",
    "    k = int(sqrt(len(combined)))\n",
    "    \n",
    "    distances = []\n",
    "    for j in range(len(combined)):\n",
    "        other = combined[j]\n",
    "        if individual == other:\n",
    "            continue\n",
    "        distances.append(1/(2.0 + get_distance(individual, other)))\n",
    "    distances.sort()\n",
    "    individual.distance = distances[k]        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To put it all together, the following function calculates the fitness (strength) of all individuals given a population and archive:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def calculate_strength(population, archive):\n",
    "    combined = population + archive\n",
    "\n",
    "    for x in combined:\n",
    "        score(x, combined)\n",
    "    for x in combined:\n",
    "        raw_fitness(x, combined)\n",
    "    for x in combined:\n",
    "        density(x, combined)\n",
    "    \n",
    "    for individual in combined:\n",
    "        individual.fitness = individual.raw_fitness + individual.distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The other unique aspect of this algorithm is the handling of the archive. The size of the archive is constant over time, which means that sometimes individuals need to be removed or added in order to adjust the size. \n",
    "\n",
    "The _truncate_ function removes individuals from the archive; at each step, the individual which has the minimum distance to another individual is chosen for removal until the archive is no longer too large. When multiple individuals have the same distance to their closest neighbour, they are compared against the next closest neighbour until a difference is found. For this, we need a custom comparison operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def compare_distances(distances1, distances2):\n",
    "    for i in range(len(distances1)):\n",
    "        if distances1[i] != distances2[i]:\n",
    "            return distances1[i] - distances2[i]\n",
    "    return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The truncate function now needs to determine all distances and then sort the individuals in order to decide which ones to drop:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from functools import cmp_to_key\n",
    "def truncate(archive):\n",
    "    # Remove individual with minimum distance\n",
    "    num_remove = abs(population_size - len(archive))\n",
    "\n",
    "    for i in range(len(archive)):\n",
    "        individual = archive[i]\n",
    "        distances = []\n",
    "        for j in range(len(archive)):\n",
    "            if i == j:\n",
    "                continue\n",
    "            other = archive[j]\n",
    "            distances.append(get_distance(individual, other))\n",
    "        distances.sort()\n",
    "        individual.distances = distances\n",
    "        \n",
    "        # This ignores the case that multiple individuals have the same distance\n",
    "        # in which case we need to look at the next distance\n",
    "\n",
    "    def compare(x, y):\n",
    "        return compare_distances(x.distances, y.distances)\n",
    "    archive.sort(key=cmp_to_key(compare))\n",
    "    for i in range(num_remove):\n",
    "        del archive[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Since Python >= 3 only supports the use of keys to sort lists, this function uses some extra code to wrap out `compare_distances` function as a key function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If the archive is too small, then we fill it with the best dominated individuals from the population. Thus, we sort the population by fitness (i.e., strength), and then pick from the sorted list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def pad_archive(next_archive, population, archive):\n",
    "    num_missing = population_size - len(next_archive)\n",
    "\n",
    "    # fill with dominated individuals in archive and population\n",
    "    candidates = [i for i in population+archive if i.fitness >= 1]\n",
    "    candidates.sort(key = lambda r : r.fitness)\n",
    "    next_archive.extend(candidates[:num_missing])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The mating is the same we already know from other evolutionary algorithms. The selection operator used is a binary tournament similar to what NSGA-II uses. However, unlike the binary tournament in NSGA-II, the tournament is now decided by the strength (fitness) function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def binary_tournament_fitness(population):\n",
    "    individual1 = random.choice(population)\n",
    "    individual2 = random.choice(population)\n",
    "\n",
    "    if individual1.fitness < individual2.fitness:\n",
    "        return individual1\n",
    "    elif individual1.fitness > individual2.fitness:\n",
    "        return individual2\n",
    "    else:\n",
    "        return random.choice([individual1, individual2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Selected individuals are subjected to crossover and mutation as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def generate_offspring_spea2(population):\n",
    "    offspring_population = []\n",
    "    while len(offspring_population) < population_size:\n",
    "        parent1, parent2 = binary_tournament_fitness(population), binary_tournament_fitness(population)\n",
    "        if random.random() < P_xover:\n",
    "            offspring1, offspring2 = crossover(parent1, parent2)\n",
    "        else:\n",
    "            offspring1, offspring2 = parent1, parent2\n",
    "        offspring1 = mutate(offspring1)\n",
    "        offspring2 = mutate(offspring2)\n",
    "\n",
    "        offspring_population.append(offspring1)\n",
    "        offspring_population.append(offspring2)\n",
    "\n",
    "    return offspring_population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we have all the components in place and can implement the overall algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def spea2():\n",
    "    initialise_plot()\n",
    "    population = [get_random_individual() for _ in range(population_size)]\n",
    "    archive = []\n",
    "\n",
    "    for iteration in range(max_gen):\n",
    "        plot(population+archive)\n",
    "\n",
    "        calculate_strength(population, archive)\n",
    "\n",
    "        # Copy all non-dominated individuals in population and archive to next archive\n",
    "        next_archive = []\n",
    "        for i in population + archive:\n",
    "            if i.fitness < 1:\n",
    "                next_archive.append(i)\n",
    "\n",
    "        if len(next_archive) < population_size:\n",
    "            pad_archive(next_archive, population, archive)\n",
    "\n",
    "        if len(next_archive) > population_size:\n",
    "            truncate(next_archive)\n",
    "\n",
    "\n",
    "        # Mating selection:\n",
    "        population = generate_offspring_spea2(next_archive)\n",
    "        archive = next_archive\n",
    "\n",
    "        # Termination: return non-dominated individuals in next_archive\n",
    "        result = [p for p in next_archive if p.fitness < 1]\n",
    "\n",
    "    plot(result)\n",
    "\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = spea2()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As usual we add some datapoints for this algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fronts_spea2 = run_times(spea2, 10)\n",
    "results[\"SPEA2\"] = [hypervolume(front, (0,total_costs)) for front in fronts_spea2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Two Archives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The Two Archives algorithm generalises the idea of using archives and uses one for convergence, and another one for diversity:\n",
    "\n",
    "- If a new solution is not dominated by both archives and dominates at least one solution in either archive, it goes into the convergence archive. \n",
    "- If a new solution is not dominated by both archives but fails to dominate any solution in either archive, it goes into the diversity archive. \n",
    "- When the archives get full, the convergence archive is preserved while the diversity archive is pruned based on crowding distance. \n",
    "\n",
    "Praditwong, K., & Yao, X. (2006, November). A new multi-objective evolutionary optimisation algorithm: The two-archive algorithm. In 2006 International Conference on Computational Intelligence and Security (Vol. 1, pp. 286-291). IEEE."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The general workflow is as follows:\n",
    "1. Initialise the population (as usual)\n",
    "2. Initialise both archives as empty sets\n",
    "3. Evaluate the initial population\n",
    "4. Repeat until done:\n",
    "5. Collect non-dominated individuals to the archives\n",
    "6. Generate new population with parents from archives\n",
    "7. Evaluate new population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "During selection, one of the archives is probabilistically chosen. The probability is a pre-defined parameter that is a ratio to choose members from the convergence archive to the diversity archive. A member in the chosen archive is selected uniformly at random."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def twoarchive_selection(convergence_archive, diversity_archive):\n",
    "    if not convergence_archive:\n",
    "        return random.choice(diversity_archive)\n",
    "    if not diversity_archive:\n",
    "        return random.choice(convergence_archive)\n",
    "    if random.random() < 0.5:\n",
    "        return random.choice(convergence_archive)\n",
    "    else:\n",
    "        return random.choice(diversity_archive)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The selection is the only difference in terms of reproduction compared to other evolutionary algorithms we have seen previously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def twoarchive_generate_offspring(convergence_archive, diversity_archive):\n",
    "    offspring_population = []\n",
    "    while len(offspring_population) < population_size:\n",
    "        parent1 = twoarchive_selection(convergence_archive, diversity_archive)\n",
    "        parent2 = twoarchive_selection(convergence_archive, diversity_archive)\n",
    "        if random.random() < P_xover:\n",
    "            offspring1, offspring2 = crossover(parent1, parent2)\n",
    "        else:\n",
    "            offspring1, offspring2 = parent1, parent2\n",
    "        offspring1 = mutate(offspring1)\n",
    "        offspring2 = mutate(offspring2)\n",
    "\n",
    "        offspring_population.append(offspring1)\n",
    "        offspring_population.append(offspring2)\n",
    "\n",
    "    return offspring_population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We need to select all individuals from the resulting population that are not dominated by the archives. For this, we need some helper functions. First, we define a helper function that tells us whether a candidate solution is dominated by an archive:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A new individual is first compared with all members in the current archives. If it is dominated by a member of the archives, it is discarded, otherwise it becomes a new member of the archives.\n",
    "\n",
    "The remainder of the archives are compared with the new member and two cases are possible:\n",
    "- If the new member dominates a member of the archives then the dominated member is removed and the new member is received by the _convergence_ archive \n",
    "- If the new member does not dominate any members and is not dominated by any archive members then it becomes a member of the diversity archive, and the size of the diversity archive is increased. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def collect_nondominated(individual, ca, da):\n",
    "    # TODO: During this stage, any duplicated member is deleted.\n",
    "\n",
    "    if is_dominated(individual, ca) or is_dominated(individual, da):\n",
    "        return\n",
    "\n",
    "    dominated_ca = get_dominated(individual, ca)\n",
    "    dominated_da = get_dominated(individual, da)\n",
    "    if len(dominated_ca) + len(dominated_da) > 0:\n",
    "        for i in dominated_ca:\n",
    "            ca.remove(i)\n",
    "        for i in dominated_da:\n",
    "            da.remove(i)\n",
    "        ca.append(individual)\n",
    "    else:\n",
    "        # No individuals are dominated\n",
    "        da.append(individual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If the total size of the archives overflows, then we need to remove individuals. The removal operator deletes only members in the diversity archive and has no impact on the convergence archive. \n",
    "\n",
    "To select which members to remove from the diversity archive, we calculate the shortest distance to the nearest member in the convergence archive for all members in the diversity archive. The member with the shortest distance among the diversity members is deleted until the total size equals the capacity. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "archive_size = 40\n",
    "\n",
    "def remove(ca, da):\n",
    "    if len(ca) + len(da) > archive_size:\n",
    "        for individual in da:\n",
    "            individual.length = sys.maxsize\n",
    "            for other in ca:\n",
    "                d = get_distance(individual, other)\n",
    "                if individual.length > d:\n",
    "                    individual.length = d\n",
    "        da.sort(key=lambda r: r.length)\n",
    "        while len(ca) + len(da) > archive_size:\n",
    "            del da[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we just need to put everything together. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def twoarchives():\n",
    "    initialise_plot()\n",
    "    population = [get_random_individual() for _ in range(population_size)]\n",
    "    convergence_archive = []\n",
    "    diversity_archive   = []\n",
    "\n",
    "    for iteration in range(max_gen):\n",
    "        plot(population+convergence_archive+diversity_archive)\n",
    "\n",
    "        # Collect non-dominated individuals\n",
    "        for i in get_nondominated(population):\n",
    "            collect_nondominated(i, convergence_archive, diversity_archive)\n",
    "\n",
    "        # Truncate archive sizes\n",
    "        remove(convergence_archive, diversity_archive)\n",
    "\n",
    "        next_generation = twoarchive_generate_offspring(convergence_archive, diversity_archive)\n",
    "        population = next_generation\n",
    "\n",
    "    result = get_nondominated(convergence_archive)\n",
    "    plot(result)\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = twoarchives()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fronts_ta = run_times(twoarchives, 10)\n",
    "results[\"TwoArchives\"] = [hypervolume(front, (0,total_costs)) for front in fronts_ta]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## SMS-EMOA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We've been comparing the multi-objective algorithms in terms of the hypervolume, although they used domination and diversity measurements to guide the search. An alternative is to use the hypervolume _directly_ to drive the search. The hypervolume is often referred to as S metric, and the idea to optimise the hypervolume is implemented in the S metric selection evolutionary multi objective algorithm (SMS-EMOA):\n",
    "\n",
    "Beume, N., Naujoks, B., & Emmerich, M. (2007). SMS-EMOA: Multiobjective selection based on dominated hypervolume. European Journal of Operational Research, 181(3), 1653-1669."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The overall principle of SMS-EMOA is the following:\n",
    "1. Generate a random population of individuals\n",
    "2. Repeat until done:\n",
    "3. Generate an offspring by variation\n",
    "4. Replace existing member of population if it improves the hypervolume"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The _reduce_ function sorts the population using fast non-dominated sort (as introduced with NSGA-II), and then discards one individual from the worst front. The choice of which individual to discard is determined by the contribution the hypervolume: The individual contributing least to the hypervolume is the one discarded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def reduce(population):\n",
    "    fronts = fast_non_dominated_sort(population)\n",
    "    last_front = [population[index] for index in fronts[-1]]\n",
    "    if len(last_front) > 1:\n",
    "        hypervolumes = []\n",
    "        for i in range(len(last_front)):\n",
    "            front = last_front[:i] + last_front[i+1:]\n",
    "            hypervolumes.append(hypervolume(front, (0,total_costs)))\n",
    "        chosen = last_front[hypervolumes.index(max(hypervolumes))]\n",
    "        population.remove(chosen)\n",
    "    else:\n",
    "        population.remove(last_front[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Since the calculation of the hypervolume can be computationally expensive,the algorithm is implemented as a steady state algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def sms():\n",
    "    initialise_plot()\n",
    "    population = [get_random_individual() for _ in range(population_size)]\n",
    "\n",
    "    for step in range(max_gen * population_size):\n",
    "        if step % population_size == 0:\n",
    "            plot(population)\n",
    "\n",
    "        parent1, parent2 = random.choice(population), random.choice(population)\n",
    "        if random.random() < P_xover:\n",
    "            offspring1, offspring2 = crossover(parent1, parent2)\n",
    "        else:\n",
    "            offspring1, offspring2 = parent1, parent2\n",
    "        \n",
    "        population.append(mutate(random.choice([offspring1, offspring2])))\n",
    "        reduce(population)\n",
    "\n",
    "    result = get_nondominated(population)\n",
    "    plot(result)\n",
    "    return result "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = sms()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fronts_sms = run_times(sms, 10)\n",
    "results[\"SMS-EMOA\"] = [hypervolume(front, (0,total_costs)) for front in fronts_sms]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This tells us which algorithm is best overall on our first NRP problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Alternative Next Release Problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "profit_map = {}        # customer id => weight\n",
    "requirements_map = {}  # customer id => [ requirement_id * ]\n",
    "cost_map = {}          # requirement id => cost\n",
    "dependency_map = {}    # requirement id => [ requirement_id * ]\n",
    "num_requirements = 620\n",
    "parse_nrp(\"data/nrp/nrp2.txt\")\n",
    "total_costs = sum([ cost_map[i+1] for i in range(num_requirements)])\n",
    "total_profit = sum(profit_map.values())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "results = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fronts_random = run_times(random_moo, 10)\n",
    "fronts_nsgaii = run_times(nsgaii, 10)\n",
    "fronts_paes = run_times(paes, 10)\n",
    "fronts_spea2 = run_times(spea2, 10)\n",
    "fronts_ta = run_times(twoarchives, 10)\n",
    "fronts_sms = run_times(sms, 10)\n",
    "\n",
    "results[\"Random\"]   = [hypervolume(front, (0,total_costs)) for front in fronts_random]\n",
    "results[\"NSGA-II\"]  = [hypervolume(front, (0,total_costs)) for front in fronts_nsgaii]\n",
    "results[\"PAES\"] = [hypervolume(front, (0,total_costs)) for front in fronts_paes]\n",
    "results[\"SPEA2\"] = [hypervolume(front, (0,total_costs)) for front in fronts_spea2]\n",
    "results[\"TwoArchives\"] = [hypervolume(front, (0,total_costs)) for front in fronts_ta]\n",
    "results[\"SMS-EMOA\"] = [hypervolume(front, (0,total_costs)) for front in fronts_sms]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There is a range of benchmark problems to compare multi-objective optimisation algorithms on. For example, we will consider the ZDT set of functions defined in the following paper:\n",
    "\n",
    "Chase, N., Rademacher, M., Goodman, E., Averill, R., & Sidhu, R. (2009). A benchmark study of multi-objective optimization methods. BMK-3021, Rev, 6, 1-24."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The ZDT1 function has a convex Pareto-optimal front. The objective functions are\n",
    "\n",
    "$f_1(x) = x_1$\n",
    "$f_2(x) = g(x) [ 1 - \\sqrt{x_1 / g(x)} ]$\n",
    "\n",
    "where $g(x)$ is defined as:\n",
    "\n",
    "$g(x) = 1 + 9 (\\Sigma^n_{i=2} x_i) / (n-1)$.\n",
    "\n",
    "Individuals are vectors of floating point numbers of size $n$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n = 30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def function1(individual):\n",
    "    return individual[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def function2(individual):\n",
    "    x0 = individual[0]\n",
    "\n",
    "    sum = 0.0;\n",
    "    for x in individual[1:]:\n",
    "        sum += x\n",
    "        \n",
    "    g = 1.0 + 9.0 * sum / (len(individual) - 1)\n",
    "\n",
    "    return g * (1 - sqrt(x0 / g))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our representation currently only gives us bitvectors, so we need to adapt the search operators to produce vectors of floating point numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_individual():\n",
    "    individual = L(random.random() for _ in range(n))\n",
    "    evaluate(individual)\n",
    "    return individual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To mutate a floating point number, we add some random noise using a Gaussian distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mutate(solution):\n",
    "    P_mutate = 1/len(solution)\n",
    "    mutated = L(solution[:])\n",
    "    for position in range(len(solution)):\n",
    "        if random.random() < P_mutate:\n",
    "            mutated[position] = max(0, min(1, mutated[position] + random.gauss(0, 0.05)))\n",
    "    evaluate(mutated)\n",
    "    return mutated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Both our objectives are minimisation problems, so we need to update the dominance relation accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def dominates(solution1, solution2):\n",
    "    if solution1.fitness1 > solution2.fitness1 or solution1.fitness2 > solution2.fitness2:\n",
    "        return False\n",
    "\n",
    "    if solution1.fitness1 < solution2.fitness1 or solution1.fitness2 < solution2.fitness2:\n",
    "        return True\n",
    "\n",
    "    return False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We also need to remove the normalisation based on total profit and costs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_fitness1 = function1([1 for _ in range(n)])\n",
    "max_fitness2 = function2([1 for _ in range(n)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def hypervolume(front, r):\n",
    "    front.sort(key=lambda i: i.fitness1)\n",
    "\n",
    "    hv = (abs(r[0] - front[0].fitness1) / max_fitness1) * (abs(r[1] - front[0].fitness2) / max_fitness2)\n",
    "    \n",
    "    for i in range(1, len(front)):\n",
    "        hv += (abs(front[i-1].fitness1 - front[i].fitness1) / max_fitness1) * (abs(r[1] - front[i].fitness2) / max_fitness2)\n",
    "    \n",
    "    return hv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Finally, we also need to fix the axis labels on our animation plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "ims = []  # global variable to store images of the animation\n",
    "    \n",
    "def initialise_plot():\n",
    "    global ims\n",
    "    global fig\n",
    "    global ax\n",
    "    \n",
    "    ims = []\n",
    "\n",
    "    %matplotlib agg\n",
    "    fig, ax = plt.subplots()\n",
    "    plt.xlabel('Function 1', fontsize=15)\n",
    "    plt.ylabel('Function 2', fontsize=15)\n",
    "    ims = []\n",
    "    %matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "random_moo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = nsgaii()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "results = {}\n",
    "\n",
    "fronts_random = run_times(random_moo, 10)\n",
    "fronts_nsgaii = run_times(nsgaii, 10)\n",
    "fronts_paes = run_times(paes, 10)\n",
    "fronts_spea2 = run_times(spea2, 10)\n",
    "fronts_ta = run_times(twoarchives, 10)\n",
    "fronts_sms = run_times(sms, 10)\n",
    "\n",
    "results[\"Random\"]   = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_random]\n",
    "results[\"NSGA-II\"]  = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_nsgaii]\n",
    "results[\"PAES\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_paes]\n",
    "results[\"SPEA2\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_spea2]\n",
    "results[\"TwoArchives\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_ta]\n",
    "results[\"SMS-EMOA\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_sms]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For ZDT4, $f_1$ remains the same, but $f_2 = g(x)[1 - \\sqrt{x_1/g(x)}]$, where $g(x) = 1 + 10(n - 1) + \\Sigma^n_{i=2} [x_i^2 - 10 cos (4\\pi{}x_i)]$. This results in a highly multi-modal fitness landscape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from math import cos, pi\n",
    "\n",
    "def function2(individual):\n",
    "    x0 = individual[0]\n",
    "\n",
    "    g = 0\n",
    "    for x in individual[1:]:\n",
    "        g += pow(x, 2) - 10.0 * cos(4.0 * pi * x / 180.0)\n",
    "        \n",
    "    g += 1.0 + 10.0 * (len(individual) - 1)\n",
    "\n",
    "    return g * (1.0 - sqrt(x0 / g))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "max_fitness1 = function1([1 for _ in range(n)])\n",
    "max_fitness2 = function2([1 for _ in range(n)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "results = {}\n",
    "\n",
    "fronts_random = run_times(random_moo, 10)\n",
    "fronts_nsgaii = run_times(nsgaii, 10)\n",
    "fronts_paes = run_times(paes, 10)\n",
    "fronts_spea2 = run_times(spea2, 10)\n",
    "fronts_ta = run_times(twoarchives, 10)\n",
    "fronts_sms = run_times(sms, 10)\n",
    "\n",
    "results[\"Random\"]   = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_random]\n",
    "results[\"NSGA-II\"]  = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_nsgaii]\n",
    "results[\"PAES\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_paes]\n",
    "results[\"SPEA2\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_spea2]\n",
    "results[\"TwoArchives\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_ta]\n",
    "results[\"SMS-EMOA\"] = [hypervolume(front, (max_fitness1,max_fitness2)) for front in fronts_sms]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.boxplot(results.values())\n",
    "ax.set_xticklabels(results.keys())"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}