Evolutionary Search - Part 2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Evolutionary Algorithms (Part 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In this chapter we take a closer look at some of the search operators and algorithmic choices that we took for granted so far. In particular, we reconsider selection, crossover, mutation, and the topology of the population itself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import math\n",
    "import sys\n",
    "from IPython.utils import io\n",
    "from statistics import mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will start by considering the one max problem again, where a solution is a vector of length _n_, consisting of binary numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "n = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our representation is a simple list of bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_solution():\n",
    "    individual = [random.choice([0,1]) for _ in range(n)]\n",
    "    return individual"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The fitness function as well as variation operators are still the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(solution):\n",
    "    return sum(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We start with the operators that we have used previously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mutate(individual):\n",
    "    P_mutate = 1/len(individual)\n",
    "    copy = individual[:]\n",
    "    for position in range(len(individual)):\n",
    "        if random.random() < P_mutate:\n",
    "            copy[position] = 1 - copy[position]\n",
    "    return copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def singlepoint_crossover(parent1, parent2):\n",
    "    pos = random.randint(0, len(parent1))\n",
    "    offspring1 = parent1[:pos] + parent2[pos:]\n",
    "    offspring2 = parent2[:pos] + parent1[pos:]\n",
    "    return offspring1, offspring2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The first selection operator we considered was tournament selection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "tournament_size = 3\n",
    "def tournament_selection(population, replacement=False):\n",
    "    if replacement:\n",
    "        candidates = random.choices(population, k = tournament_size)\n",
    "    else:\n",
    "        candidates = random.sample(population, tournament_size)\n",
    "        \n",
    "    winner = max(candidates, key = lambda x: get_fitness(x))\n",
    "    return winner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A general problem in evolutionary search is finding the right balance between exploration and exploitation. Using the wrong operators may lead to premature convergence, or the search may never converge at all at an optimum. As a first step towards understanding what is happening inside the population of a genetic algorithm, we will consider the _average fitness_ of the population as well as the _diversity_ within the population in addition to the best fitness value we already tracked in the past. For this, we first define the difference between two individuals as the hamming distance between the vector representations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def hamming_distance(individual1, individual2):\n",
    "    return sum(c1 != c2 for c1, c2 in zip(individual1, individual2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can now calculate an overall population diversity as the sum of pairwise hamming distances for all pairs of individuals in the population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def pairwise_distance(population):\n",
    "    distances = 0\n",
    "    for i in range(len(population)-1):\n",
    "        for j in range(i, len(population)):\n",
    "           distances += hamming_distance(population[i], population[j]) \n",
    "    return distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We need to set the parameters of our genetic algorithm first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "population_size = 20\n",
    "P_xover   = 0.7\n",
    "max_steps = 1000\n",
    "selection = tournament_selection\n",
    "crossover = singlepoint_crossover\n",
    "\n",
    "# These lists track values throughout the evolution\n",
    "# so we can compare the behaviour of different operators\n",
    "fitness_values = []\n",
    "diversity_values = []\n",
    "mean_fitness_values = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the lecture we discussed that there are many different variations of all the search operators involved in a genetic algorithm. We will now look at a couple of relevant search operators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Survivor Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Most evolutionary algorithms use a fixed population size, so we need a way of going from (parents + offspring) to the next generation. \n",
    "- In age-based selection one usually produces as many offspring as there are parents, and then replaces the parents with the offspring. \n",
    "- In fitness-based selection, we rank the parents and the offspring, and take the best of all.\n",
    "\n",
    "We saw both versions last time in the context of evolution strategies, but can also create versions of a genetic algorithm. In a generational genetic algorithm the offspring replaces the parents, while in a steady state genetic algorithm we apply fitness-based survivor selection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Elitism is a special case, where the best individuals of the population always survive, while other means are used for the rest of the population. We implement a simple version of elitism by simply ranking the population by diversity and taking the top `elite_size` elements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "elite_size = int(population_size * 0.05)\n",
    "\n",
    "def elitism(population):\n",
    "    population.sort(key=lambda k: get_fitness(k), reverse=True)\n",
    "    return population[:elite_size]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's revisit the standard genetic algorithm with a generational selection model, integrated with elitism:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def ga():\n",
    "    population = [get_random_solution() for _ in range(population_size)]\n",
    "    fitness_values.clear()\n",
    "\n",
    "    # This could probably be written in a single line, but let's keep it explicit\n",
    "    best_fitness = -1\n",
    "    for p in population:\n",
    "        fitness = get_fitness(p)\n",
    "        if fitness > best_fitness:\n",
    "            best_fitness = fitness\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    diversity_values.append(pairwise_distance(population))\n",
    "    mean_fitness_values.append(mean([get_fitness(x) for x in population]))\n",
    "\n",
    "    while len(fitness_values) < max_steps:\n",
    "\n",
    "        new_population = elitism(population)\n",
    "        while len(new_population) < len(population):\n",
    "            parent1 = selection(population)\n",
    "            parent2 = selection(population)\n",
    "\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",
    "            fitness1, fitness2 = get_fitness(offspring1), get_fitness(offspring2)\n",
    "\n",
    "            if fitness1 > best_fitness:\n",
    "                best_fitness = fitness1\n",
    "            fitness_values.append(best_fitness)\n",
    "            if fitness2 > best_fitness:\n",
    "                best_fitness = fitness2\n",
    "            fitness_values.append(best_fitness)\n",
    "            \n",
    "            new_population += [offspring1, offspring2]\n",
    "\n",
    "        population = new_population\n",
    "        diversity_values.append(pairwise_distance(population))\n",
    "        mean_fitness_values.append(mean([get_fitness(x) for x in population]))\n",
    "\n",
    "    return max(population, key=lambda k: get_fitness(k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The alternative was the steady state genetic algorithm, where we select two parents, derive their offspring, and then do fitness-based survivor selection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def steadystatega():\n",
    "    population = [get_random_solution() for _ in range(population_size)]\n",
    "    best_fitness = -1\n",
    "    fitness_values.clear()\n",
    "    \n",
    "    for p in population:\n",
    "        fitness = get_fitness(p)\n",
    "        if fitness > best_fitness:\n",
    "            best_fitness = fitness\n",
    "            best_solution = p\n",
    "        fitness_values.append(best_fitness)\n",
    "    diversity_values.append(pairwise_distance(population))\n",
    "    mean_fitness_values.append(mean([get_fitness(x) for x in population]))\n",
    "    \n",
    "    while len(fitness_values) < max_steps:                \n",
    "        parent1 = selection(population)\n",
    "        parent2 = selection(population)\n",
    "        \n",
    "        p1 = population.index(parent1)\n",
    "        p2 = population.index(parent2)\n",
    "\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",
    "        best1, best2 = sorted([parent1, parent2, offspring1, offspring2], key=lambda x: get_fitness(x), reverse=True)[:2]\n",
    "        population[p1] = best1\n",
    "        population[p2] = best2\n",
    "\n",
    "        fitness1, fitness2 = get_fitness(best1), get_fitness(best2)\n",
    "\n",
    "        if fitness1 > best_fitness:\n",
    "            best_fitness = fitness1\n",
    "        fitness_values.append(best_fitness)\n",
    "        if fitness2 > best_fitness:\n",
    "            best_fitness = fitness2\n",
    "        fitness_values.append(best_fitness)\n",
    "        \n",
    "        # To make plots comparable with the generational GA\n",
    "        if len(fitness_values) % population_size == 0:\n",
    "            diversity_values.append(pairwise_distance(population))\n",
    "            mean_fitness_values.append(mean([get_fitness(x) for x in population]))\n",
    "            \n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Note that there are further possibilities for variation here: Instead of replacing the selected parents in the population, we could select other individuals (e.g., the worst individuals)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To consider the effects on diversity, let's compare diversity and fitness throughout one run each. We'll increase `n` to make the problem slightly more challenging."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "n = 100\n",
    "max_steps = 1000\n",
    "population_size = 20\n",
    "\n",
    "with io.capture_output() as captured: \n",
    "    elite_size = 0\n",
    "    fitness_values = []\n",
    "    diversity_values = []\n",
    "    mean_fitness_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"No elitism\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"No elitism\")\n",
    "    axes[2].plot(fitness_values, label=f\"No elitism\")\n",
    "\n",
    "    elite_size = 1\n",
    "    fitness_values = []\n",
    "    diversity_values = []\n",
    "    mean_fitness_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"E1\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"E1\")\n",
    "    axes[2].plot(fitness_values, label=f\"E1\")\n",
    "    \n",
    "    elite_size = 10\n",
    "    fitness_values = []\n",
    "    diversity_values = []\n",
    "    mean_fitness_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"E10\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"E10\")\n",
    "    axes[2].plot(fitness_values, label=f\"E10\")\n",
    "    \n",
    "    fitness_values = []\n",
    "    diversity_values = []\n",
    "    mean_fitness_values = []\n",
    "    steadystatega()\n",
    "    axes[0].plot(diversity_values, label=f\"SS-GA\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"SS-GA\")\n",
    "    axes[2].plot(fitness_values, label=f\"SS-GA\")\n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As usual results may vary between runs, since these are randomised algorithms. However, a general trend we should see in the above plots is that the population diversity of the steady state GA reduces much slower than in a generational GA, and also the average fitness value in the population remains lower. The large elitism size of `10` means that the algorithm can run more generations with the same number of fitness evaluations, which is why it continues longer than the others in the first two plots (which shows generations rather than fitness evaluations). A small elitism size tends to generally lead to a better average fitness in this configuration. The best performing version will differ between runs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Parent Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A major difference between the evolution strategies we considered initially and the canonical genetic algorithm we looked at afterwards is the parent selection strategy. In classical evolution strategies all μ (mu) parents are involved in recombination, and the survivor selection is what drives the selective pressure. In genetic algorithms, instead, the parent selection applies selective pressure. \n",
    "\n",
    "We started off with tournament selection because it is the quickest to implement. A traditionally more common variant is fitness proportionate selection, where the probability of an individual to be selected is proportional to its fitness value. The selection thus first calculates the total fitness sum, and then probabilistically chooses an individual by sampling a number in the range between 0 and the total fitness sum. An important requirement is that the population is sorted by fitness values, starting with the best individual (largest fitness).  This selection operator is also known as _roulette wheel selection_.\n",
    "\n",
    "In our simple implementation, we create a list of tuples `fitness_list` that stores individuals with their fitness. Obviously, there's some redundant fitness calculations here; in practice one would cache fitness values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def roulette_selection(population):\n",
    "    fitness_sum = sum([get_fitness(x) for x in population])\n",
    "    population.sort(key=lambda x: get_fitness(x), reverse=True)\n",
    "    pick = random.uniform(0, fitness_sum)\n",
    "    current = 0\n",
    "    for x in population:\n",
    "        current += get_fitness(x)\n",
    "        if current > pick:\n",
    "            return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To evaluate this, let's create a simple example population for `n=5` with individuals with fitness 5, 4, 3, 2, 1, and 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "example_population = [ [1,1,1,1,1], [0,1,1,1,1], [0,0,1,1,1], [0,0,0,1,1], [0,0,0,0,1], [0,0,0,0,0]  ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Applying this selection operator will more likely select the best individual(s) (but may select worse individuals as well). We can do a simple experiment by sampling repeatedly from our `example_population` and looking at the resulting histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hideOutput": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = roulette_selection(example_population)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A problem with fitness proportionate selection is that individuals that have a much better fitness value will dominate the selection. For example, let's skew our example population:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "example_population = [ [1,1,1,1,1], [0,1,0,0,0], [1,0,0,0,0], [0,0,1,0,0], [0,0,0,1,0],[0,0,0,0,0] ]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In this population, the first individual has fitness value 5, while the other individuals have fitness 1 and 0. The sum of fitness values is 9, and so the first individual has a probability of 56% of being selected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = roulette_selection(example_population)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Tournament selection, which we implemented earlier, suffers less from this problem. In tournament selection, we can adjust the _selective pressure_ by adjusting the tournament size. With our example population of size 6, a tournament size of 2 without replacement would imply a 33% chance of the best of the six individuals being selected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "tournament_size = 2\n",
    "\n",
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = tournament_selection(example_population)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's also compare this to the population with a more equal spread:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "example_population = [ [1,1,1,1,1], [0,1,1,1,1], [0,0,1,1,1], [0,0,0,1,1], [0,0,0,0,1], [0,0,0,0,0]  ]\n",
    "\n",
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = tournament_selection(example_population)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To see the effects of the tournament size on the selection, let's repeat this with a larger tournament size, which means higher selective pressure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "tournament_size = 4\n",
    "\n",
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = tournament_selection(example_population, replacement=True)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Even with small tournament size the worst individual was not chosen in any of these cases. The reason is that we are using tournament selection _without_ replacement. If we pick any two individuals out of `example_population`, the individual with fitness 0 will _always_ be worse. If we use replacement, then there is a chance that the worst individual gets selected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "tournament_size = 2\n",
    "\n",
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = tournament_selection(example_population, replacement=True)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can also observe the effects of the selective pressure throughout evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n = 100\n",
    "from IPython.utils import io\n",
    "\n",
    "elite_size = 1\n",
    "selection = tournament_selection\n",
    "tournament_sizes = [1, 2, 5, 10]\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "for tournament_size in tournament_sizes:\n",
    "    fitness_values, diversity_values, mean_fitness_values = [], [], []\n",
    "    with io.capture_output() as captured: \n",
    "        ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Tournament size {tournament_size}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Tournament size {tournament_size}\")\n",
    "    axes[2].plot(fitness_values, label=f\"Tournament size {tournament_size}\")\n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Too little selective pressure (e.g., tournament size of 1) tends to be bad. Very large tournaments might be too eager (which is not so much of a problem in one max though)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Should we use fitness proportionate selection or tournament selection?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "n = 100\n",
    "with io.capture_output() as captured: \n",
    "    fitness_values = []\n",
    "    diversity_values = []\n",
    "    mean_fitness_values = []\n",
    "    selection = tournament_selection\n",
    "    tournament_size = 2\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Tournament size {tournament_size}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Tournament size {tournament_size}\")\n",
    "    axes[2].plot(fitness_values, label=f\"Tournament size {tournament_size}\")\n",
    "\n",
    "    selection = roulette_selection\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Roulette wheel selection\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Roulette wheel selection\")\n",
    "    axes[2].plot(fitness_values, label=f\"Roulette wheel selection\")\n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "An alternative selection operator is _rank selection_, which is similar to fitness proportionate selection, except that the probability is calculated based on the _rank_ in the population sorted by fitness, rather than the actual fitness value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "rank_bias = 2\n",
    "def rank_selection(population):\n",
    "    population.sort(key=lambda c: get_fitness(c), reverse=True)\n",
    "    \n",
    "    individuals = []\n",
    "    N = len(population)\n",
    "    for i in range(N):\n",
    "        f2 = rank_bias - (2 * i * (rank_bias - 1))/(N - 1)\n",
    "        individuals.append((population[i], f2))\n",
    "\n",
    "    # Now implement fitness proportionate selection using the f2 values\n",
    "    fitness_sum = sum([f for (c, f) in individuals])\n",
    "    pick = random.uniform(0, fitness_sum)\n",
    "    current = 0\n",
    "    for (chromosome, fitness) in individuals:\n",
    "        current += fitness\n",
    "        if current > pick:\n",
    "            return chromosome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The rank bias, which is in the range `[1,2]` allows us to adjust the selective pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "N = 100\n",
    "for rank_bias in [1, 1.5, 2]:\n",
    "    plt.plot([rank_bias - (2 * i * (rank_bias - 1)) / (N - 1) for i in range(100)], label = f\"Rank bias {rank_bias}\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "With a bias of `2`, the worst individual has a 0% chance of being selected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "rank_bias = 2\n",
    "\n",
    "example_population = [ [1,1,1,1,1], [0,1,0,0,0], [1,0,0,0,0], [0,0,1,0,0], [0,0,0,1,0],[0,0,0,0,0] ]\n",
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = rank_selection(example_population)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "With a bias of `1`, all individuals have the same probability of being selected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "rank_bias = 1\n",
    "\n",
    "example_population = [ [1,1,1,1,1], [0,1,0,0,0], [1,0,0,0,0], [0,0,1,0,0], [0,0,0,1,0],[0,0,0,0,0] ]\n",
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = rank_selection(example_population)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will select a reasonable default for the selective pressure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "rank_bias = 1.4\n",
    "\n",
    "example_population = [ [1,1,1,1,1], [0,1,0,0,0], [1,0,0,0,0], [0,0,1,0,0], [0,0,0,1,0],[0,0,0,0,0] ]\n",
    "counts = {x: 0 for x in range(len(example_population))}\n",
    "\n",
    "for i in range(1000):\n",
    "    selected = rank_selection(example_population)\n",
    "    index = example_population.index(selected)\n",
    "    counts[index] = counts[index] + 1\n",
    "    \n",
    "plt.bar(counts.keys(), counts.values())\n",
    "plt.xticks(list(counts.keys()));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can again observe the effects of the rank bias on the evolution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.utils import io\n",
    "selection = rank_selection\n",
    "rank_biases = [1, 1.3, 1.7, 2.0]\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "for rank_bias in rank_biases:\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    with io.capture_output() as captured: \n",
    "        ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Rank bias {rank_bias}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Rank bias {rank_bias}\")\n",
    "    axes[2].plot(fitness_values, label=f\"Rank bias {rank_bias}\")\n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Some variants of genetic algorithms use selection where each individual has the same chance of being selected. Although this removes selection pressure, this is usually compensated with a strong fitness-based survivor selection mechanism."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def uniform_selection(population):\n",
    "    return random.choice(population)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We would not usually use uniform selection without some other survival selection, but to see what the effects on the search are, let's put all the options together in one experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "n = 100\n",
    "with io.capture_output() as captured: \n",
    "    rank_bias = 1.4\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    selection = tournament_selection\n",
    "    tournament_size = 2\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Rank {rank_bias}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Rank {rank_bias}\")\n",
    "    axes[2].plot(fitness_values, label=f\"Rank {rank_bias}\")\n",
    "\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    selection = tournament_selection\n",
    "    tournament_size = 2\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Tournament size {tournament_size}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Tournament size {tournament_size}\")\n",
    "    axes[2].plot(fitness_values, label=f\"Tournament size {tournament_size}\")\n",
    "\n",
    "    selection = roulette_selection\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Roulette wheel selection\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Roulette wheel selection\")\n",
    "    axes[2].plot(fitness_values, label=f\"Roulette wheel selection\")\n",
    "\n",
    "    selection = uniform_selection\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Uniform selection\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Uniform selection\")\n",
    "    axes[2].plot(fitness_values, label=f\"Uniform selection\")\n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Survivor selection and parent selection are not independent: As an example, let's consider the effects of selective pressure on the generational GA and the steady state GA. We'll run tournament selection with two different tournament sizes to represent reasonable and high selective pressure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hideOutput": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "n = 100\n",
    "with io.capture_output() as captured: \n",
    "    tournament_size = 6\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    selection = tournament_selection\n",
    "    steadystatega()\n",
    "    axes[0].plot(diversity_values, label=f\"SS-GA Tournament {tournament_size}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"SS-GA Tournament {tournament_size}\")\n",
    "    axes[2].plot(fitness_values, label=f\"SS-GA Tournament {tournament_size}\")\n",
    "\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    selection = tournament_selection\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"GA Tournament {tournament_size}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"GA Tournament {tournament_size}\")\n",
    "    axes[2].plot(fitness_values, label=f\"GA Tournament {tournament_size}\")\n",
    "\n",
    "    tournament_size = 2\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    selection = tournament_selection\n",
    "    steadystatega()\n",
    "    axes[0].plot(diversity_values, label=f\"SS-GA Tournament {tournament_size}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"SS-GA Tournament {tournament_size}\")\n",
    "    axes[2].plot(fitness_values, label=f\"SS-GA Tournament {tournament_size}\")\n",
    "\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    selection = tournament_selection\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"GA Tournament {tournament_size}\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"GA Tournament {tournament_size}\")\n",
    "    axes[2].plot(fitness_values, label=f\"GA Tournament {tournament_size}\")\n",
    "\n",
    "        \n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "With high selective pressure, the steady state GA maintains very high diversity, even though the best fitness improves reasonably over time. The reason is that with a large tournament size of 6 (with a population size of 20) the parent selection will repeatedly select the same few individuals with a high probability. If the result leads to an improvement, these are replaced but will be picked again with high probability afterwards, while large parts of the population (worse individuals that keep losing tournaments) will simply remain unchanged. This is also reflected by the lower average fitness in the population. For the generational GA the high selective pressure has the opposite effect: The population loses diversity very quickly, because most offspring will be produced from the same few parents. For one max, this doesn't appear to be a bad thing though."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Alternative Crossover Operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our crossover operator so far only considers a single point for crossing two individuals. While this is the most common variant in practice, there is one potential downside: Only locally neighbouring genetic material is preserved; if a parent has relevant genes at the beginning and the end of the chromosome, these will not be inherited to the offspring directly. One way to circumvent this is by defining more than one crossover point. For example, we can define a two-point crossover operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def twopoint_crossover(parent1, parent2):\n",
    "    pos1 = random.randint(1, len(parent1))\n",
    "    pos2 = random.randint(pos1, len(parent1))\n",
    "    offspring1 = parent1[:pos1] + parent2[pos1:pos2] + parent1[pos2:]\n",
    "    offspring2 = parent2[:pos1] + parent1[pos1:pos2] + parent2[pos2:]\n",
    "    return offspring1, offspring2\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "parent1 = [0,0,0,0,0,0,0,0,0,0]\n",
    "parent2 = [1,1,1,1,1,1,1,1,1,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the single point crossover, the offspring of `parent1` and `parent2` will _always_ be either a sequence of `0` followed by a sequence of `1`, or vice versa:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "singlepoint_crossover(parent1, parent2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In the two point crossover, there will be some variation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "twopoint_crossover(parent1, parent2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "It is not common to increase the number of crossover points beyond two, but instead if more variation is required, it is simply possible to _uniformly_ select genes from either of the parents, resulting in _uniform crossover_:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def uniform_crossover(parent1, parent2):\n",
    "    offspring1 = []\n",
    "    offspring2 = []\n",
    "    for pos in range(len(parent1)):\n",
    "        if random.choice([True, False]):\n",
    "            offspring1.append(parent1[pos])\n",
    "            offspring2.append(parent2[pos])\n",
    "        else:\n",
    "            offspring1.append(parent2[pos])\n",
    "            offspring2.append(parent1[pos])\n",
    "    return offspring1, offspring2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Applying this to `parent1` and `parent2` from above, we will see offspring consisting of more variation in `1`s and `0`s, but these will always be chosen from parents and not random:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "uniform_crossover(parent1, parent2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The importance and influence of the crossover operator on the search is part of active research, and also depends on the problem we are trying to solve. To see whether there are any benefits to using crossover on our one max example, we can conduct a _headless chicken_ test: During crossover, we use a randomly generated individual as one of the parents. If the search still performs as well, then the crossover operator actually just serves as a kind of macro-mutation. If the search no longer performs as well, then it is the actual combination of parent genetic material that leads to an improvement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def chicken_crossover(parent1, parent2):\n",
    "    offspring1 = []\n",
    "    offspring2 = []\n",
    "    parent2 = get_random_solution()\n",
    "    for pos in range(len(parent1)):\n",
    "        if random.choice([True, False]):\n",
    "            offspring1.append(parent1[pos])\n",
    "            offspring2.append(parent2[pos])\n",
    "        else:\n",
    "            offspring1.append(parent2[pos])\n",
    "            offspring2.append(parent1[pos])\n",
    "    return offspring1, offspring2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will run the usual combination of experiments and analyses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "n = 100\n",
    "with io.capture_output() as captured: \n",
    "    crossover = singlepoint_crossover\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Singlepoint\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Singlepoint\")\n",
    "    axes[2].plot(fitness_values, label=f\"Singlepoint\")\n",
    "\n",
    "    crossover = uniform_crossover\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Uniform\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Uniform\")\n",
    "    axes[2].plot(fitness_values, label=f\"Uniform\")\n",
    "    \n",
    "    crossover = twopoint_crossover\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Twopoint\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Twopoint\")\n",
    "    axes[2].plot(fitness_values, label=f\"Twopoint\")\n",
    "    \n",
    "    crossover = chicken_crossover\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Chicken\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Chicken\")\n",
    "    axes[2].plot(fitness_values, label=f\"Chicken\")\n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This experiment just conducts a single run, and as usual we need to conduct an experiment with repetitions in order to draw conclusions. However, what most likely shows is that the headless chicken test leads to substantially worse results, while the uniform crossover tends to produce the best results. Consequently, it seems that crossover actually is useful for our problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Mutation Operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The mutation operator we considered so far flips each bit in a sequence of length $n$ with a probability of $1/n$. To be able to configure our genetic algorithm with alternative mutation operators, let's redefine this operator in a differently named function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def avg_mutate(individual):\n",
    "    P_mutate = 1/len(individual)\n",
    "    copy = individual[:]\n",
    "    for position in range(len(individual)):\n",
    "        if random.random() < P_mutate:\n",
    "            copy[position] = 1 - copy[position]\n",
    "    return copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "One of the reasons we used this operator so far was that it implicitly defines a mutation probability and we had one less parameter to worry about. However, when considering alternative operators, we will require a probability for applying mutation. This probability is usually fairly small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "P_mutate = 0.05"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A basic alternative operator would be to flip a single bit, with probability `P_mutate`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def mutate1(individual):\n",
    "    copy = individual[:]\n",
    "    if random.random() < P_mutate:\n",
    "        position = random.randint(0, len(copy) - 1)\n",
    "        copy[position] = 1 - copy[position]\n",
    "    return copy "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The scope for different operators on our bitvector representation is limited. Let's consider an alternative where we flip multiple bits at the same time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def mutate10(individual):\n",
    "    copy = individual[:]\n",
    "    if random.random() < P_mutate:\n",
    "        for _ in range(10):\n",
    "            position = random.randint(0, len(copy) - 1)\n",
    "            copy[position] = 1 - copy[position]\n",
    "    return copy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Given these three mutation operators, we can now run some comparative experiments again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "crossover = singlepoint_crossover\n",
    "selection = tournament_selection\n",
    "tournament_size = 2\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "n = 100\n",
    "max_steps = 5000\n",
    "with io.capture_output() as captured: \n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    mutate = avg_mutate\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Average Mutation\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Average Mutation\")\n",
    "    axes[2].plot(fitness_values, label=f\"Average Mutation\")\n",
    "\n",
    "    mutate = mutate1\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"1 Bit\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"1 Bit\")\n",
    "    axes[2].plot(fitness_values, label=f\"1 Bit\")\n",
    "    \n",
    "    mutate = mutate10\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"10 Bit\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"10 Bit\")\n",
    "    axes[2].plot(fitness_values, label=f\"10 Bit\")\n",
    "    \n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Mutating each bit with a probability of $1/n$ is beneficial for the diversity. Comparing the operators that flip 1 and 10 bits, we typically can observe that the 10 bit flip makes larger jumps in fitness improvements, but towards the end of the search it will struggle to home in on a solution as it flips _too much_. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Cellular Genetic Algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The last aspect of the genetic algorithm we will consider is the population: The population is usually a multiset, and we mostly implemented it as a list to impose an order on individuals. A cellular evolutionary algorithm imposes a topology on the population. During reproduction, an individual is only allowed to mate with its neighbours, as defined by the topology (commonly rings, grids, or two-dimensional torus graphs)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's put our population into a grid, i.e. a two-dimensional list. Let's assume we consider only the 4-neighbourhood (von Neumann neighbourhood) during selection, then we randomly select one of these 4 neighbours and the selected position itself:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def grid_selection(population, row, col):\n",
    "    neighbours = []\n",
    "    for dx in [-1, 1]:\n",
    "        if row + dx >= 0 and row + dx < len(population):\n",
    "            neighbours.append(population[row + dx][col])\n",
    "    for dy in [-1, 0, 1]:\n",
    "        if col + dy >= 0 and col + dy < len(population):\n",
    "            neighbours.append(population[row][col + dy])\n",
    "\n",
    "    return random.choice(neighbours)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The genetic algorithm needs to be modified such that the population is a properly initialised grid, and then uses the `grid_selection` we just defined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "grid_size = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def cellular_ga():\n",
    "    fitness_values.clear()\n",
    "    population = [[get_random_solution() for _ in range(grid_size)] for _ in range(grid_size)]\n",
    "\n",
    "    diversity_values.append(pairwise_distance([y for x in population for y in x]))\n",
    "    mean_fitness_values.append(mean([get_fitness(x) for x in [y for x in population for y in x]]))\n",
    "\n",
    "    best_fitness = -1\n",
    "    for p in [y for x in population for y in x]:\n",
    "        fitness = get_fitness(p)\n",
    "        if fitness > best_fitness:\n",
    "            best_fitness = fitness\n",
    "        fitness_values.append(best_fitness)\n",
    "        \n",
    "    while len(fitness_values) < max_steps:\n",
    "        new_population = []\n",
    "        for row in range(grid_size):\n",
    "            new_population.append([])\n",
    "            for col in range(grid_size):\n",
    "                parent1 = grid_selection(population, row, col)\n",
    "                parent2 = grid_selection(population, row, col)\n",
    "\n",
    "                if random.random() < P_xover:\n",
    "                    offspring1, offspring2 = crossover(parent1, parent2)\n",
    "                else:\n",
    "                    offspring1, offspring2 = parent1[:], parent2[:]\n",
    "\n",
    "                offspring = mutate(random.choice([offspring1, offspring2]))\n",
    "\n",
    "                if get_fitness(offspring) >= get_fitness(population[row][col]):\n",
    "                    new_population[row].append(offspring)\n",
    "                else:\n",
    "                    new_population[row].append(population[row][col])  \n",
    "                \n",
    "                if get_fitness(offspring) > best_fitness:\n",
    "                    best_fitness = get_fitness(offspring)\n",
    "                fitness_values.append(best_fitness)\n",
    "\n",
    "        population = new_population\n",
    "        \n",
    "        diversity_values.append(pairwise_distance([y for x in population for y in x]))\n",
    "        mean_fitness_values.append(mean([get_fitness(x) for x in [y for x in population for y in x]]))\n",
    "\n",
    "    return max([y for x in population for y in x], key=lambda k: get_fitness(k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When producing the next generation the cellular GA iterates over the grid, and produces an offspring at each grid location using only the neighbourhood for reproduction. We have implemented an elitist approach where the grid location is only replaced with the offspring if the offspring has the same or better fitness."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "mutate = avg_mutate\n",
    "crossover = singlepoint_crossover\n",
    "selection = tournament_selection\n",
    "tournament_size = 3\n",
    "\n",
    "fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 4))\n",
    "n = 100\n",
    "max_steps = 5000\n",
    "with io.capture_output() as captured: \n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    ga()\n",
    "    axes[0].plot(diversity_values, label=f\"GA\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"GA\")\n",
    "    axes[2].plot(fitness_values, label=f\"GA\")\n",
    "\n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    steadystatega()\n",
    "    axes[0].plot(diversity_values, label=f\"SS-GA\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"SS-GA\")\n",
    "    axes[2].plot(fitness_values, label=f\"SS-GA\")\n",
    "    \n",
    "    fitness_values = []\n",
    "    mean_fitness_values = []\n",
    "    diversity_values = []\n",
    "    cellular_ga()\n",
    "    axes[0].plot(diversity_values, label=f\"Cellular GA\")\n",
    "    axes[1].plot(mean_fitness_values, label=f\"Cellular GA\")\n",
    "    axes[2].plot(fitness_values, label=f\"Cellular GA\")\n",
    "        \n",
    "\n",
    "axes[0].set_title('Population Diversity')\n",
    "axes[0].legend()\n",
    "axes[1].set_title('Mean Population Fitness')\n",
    "axes[1].legend()\n",
    "axes[2].set_title('Best Fitness Evolution')\n",
    "axes[2].legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For larger grids the cellular GA may be slower initially -- it takes longer for genetic material to spread across the population of a large grid. This, however, is also its benefit: Premature convergence is less likely to occur."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Memetic Algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Since an improvement to a solution does not give us any guarantees that we are in fact moving towards the optimum, evolutionary algorithms allow non-local moves. However, evolution may take longer to fine-tune solutions, while local search does exactly this very efficiently. The term _Memetic Algorithms_ is commonly used to denote the combination of global and local search. The term sometimes also refers to the use of instance-specific knowledge in search operators, but we will focus on the combination of global and local search in this notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As an example problem we consider a variation of the one-max problem, such that the search operators are simplistic, and we understand the search landscape well. We are still trying to optimise the bitstring to contain all `1`s, but assume we have a somewhat more complex fitness landscape. Given a bitstring encoding, let the fitness function be the following _hurdle_ function:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "$f(x) = - \\Big\\lceil \\frac{z(x)}{w} \\Big\\rceil - \\frac{rem(z(x), w)}{w}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Here $z(x)$ is the number of zeros in the bitstring $x$; $w \\in \\{ 2, 3, \\ldots, n\\}$ is the hurdle width; $rem(z(x), w)$ is the remainder of $z(x)$ divided by $w$, and $\\lceil \\cdot \\rceil$ is the ceiling function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "It will be easiest to understand what this does by looking at the resulting fitness landscape over a limited input range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "w = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "values = []\n",
    "for z in range(0, 100):\n",
    "    values.append(-math.ceil(z/w) - (z % w)/w)\n",
    "plt.plot(values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The optimal value has the fitness value $0$ (i.e. $z(0) = 0$ implies all bits are $1$). The representation is our usual bitstring, and we will use a moderate length for our bitstring:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "n = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As usual, an individual is a random sequence of bits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_solution():\n",
    "    domain = [0,1]\n",
    "    return [random.choice(domain) for _ in range(n)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The fitness function calculates the hurdle function. In order to avoid redundantly calculating fitness values, we first check if the individual already has a cached fitness value, and if so, we simply return that. If we do have to calculate a new fitness value then we append an item to our `fitness_values` list; we will use this list to store the best fitness value seen to date, so that we can plot the result afterwards. We will also use it as a stopping condition, since it reliably keeps track of how many times the fitness was actually _calculated_."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "\n",
    "def get_fitness(solution):\n",
    "    z = len(solution) - sum(solution)\n",
    "    fitness = -math.ceil(z/w) - (z % w)/w\n",
    "    \n",
    "    if not fitness_values:\n",
    "        fitness_values.append(fitness)\n",
    "    else: \n",
    "        fitness_values.append(max(fitness, fitness_values[-1]))\n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_fitness(get_random_solution())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's first consider how local search fares with our hurdle problem. We define a neighbourhood and recall our basic hillclimber."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_neighbours(candidate):\n",
    "    neighbours = []\n",
    "    for pos in range(len(candidate)):\n",
    "        copy = candidate[:]\n",
    "        copy[pos] = 1 - copy[pos]\n",
    "        neighbours.append(copy)\n",
    "    return neighbours"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As a small change to prior implementations, we will use the `fitness_values` list as a stopping condition, and make sure that we don't exceed `max_steps` fitness evaluations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "max_steps = 50000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def hillclimbing_restart():\n",
    "    fitness_values.clear()\n",
    "    current = get_random_solution()\n",
    "    best = current[:]\n",
    "    best_fitness = get_fitness(current)\n",
    "\n",
    "    while len(fitness_values) < max_steps:\n",
    "\n",
    "        best_neighbour = None\n",
    "        neighbour_fitness = -sys.maxsize\n",
    "        for neighbour in get_neighbours(current):\n",
    "            fitness = get_fitness(neighbour)\n",
    "\n",
    "            if fitness > neighbour_fitness:\n",
    "                best_neighbour = neighbour\n",
    "                neighbour_fitness = fitness\n",
    "\n",
    "        # Random restart if no neighbour is better\n",
    "        if neighbour_fitness <= get_fitness(current):\n",
    "            current = get_random_solution()\n",
    "            neighbour_fitness = get_fitness(current)\n",
    "        else:\n",
    "            current = best_neighbour                        \n",
    "\n",
    "        if neighbour_fitness > best_fitness:\n",
    "            best = current[:]\n",
    "            best_fitness = neighbour_fitness\n",
    "\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "hillclimbing_restart()\n",
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The hillclimber is quite unlikely to reach the optimal solution: The fitness landscape contains many local optima, and a random (re-)start needs to luckily jump beyond the last hurdle for the gradient to point to the global optimum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's compare this to a basic evolutionary search algorithm; to keep things simple initially we will just use a (1+1)EA, so we just need to add a mutation function that implements the usual probabilistic bitflips:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def mutate(solution):\n",
    "    P_mutate = 1/len(solution)\n",
    "    mutated = solution[:]\n",
    "    for position in range(len(solution)):\n",
    "        if random.random() < P_mutate:\n",
    "            mutated[position] = 1 - mutated[position]\n",
    "    return mutated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def oneplusoneea():\n",
    "    fitness_values.clear()\n",
    "    current = get_random_solution()\n",
    "    fitness = get_fitness(current)\n",
    "\n",
    "    while len(fitness_values) < max_steps:\n",
    "        candidate = mutate(current)\n",
    "        candidate_fitness = get_fitness(candidate)\n",
    "        if candidate_fitness >= fitness:\n",
    "            fitness = candidate_fitness\n",
    "            current = candidate\n",
    "\n",
    "    return candidate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "oneplusoneea()\n",
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The hurdles are problematic also for this global search algorithm: The search not only needs to jump over all hurdles, but it also needs to jump higher than the preceding hurdles. Again the result tends to be that the algorithm does not find an optimal solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A memetic algorithm combines these two algorithms, and can overcome these problems: The global search can effectively jump over hurdles, and the local search can climb the next hurdle. As an initial memetic algorithm we create a simple (1+1)MA. For this we will adapt the hillclimber so that it takes a starting point of the search as a parameter. To avoid spending too much time exploring the neighbourhood, we will also make this a first-ascent hillclimber:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def hillclimbing(starting_point):\n",
    "\n",
    "    current =  starting_point\n",
    "    fitness = get_fitness(current)\n",
    "\n",
    "    improved = True\n",
    "    while improved and len(fitness_values) < max_steps:\n",
    "        improved = False\n",
    "        \n",
    "        for neighbour in get_neighbours(current):\n",
    "            neighbour_fitness = get_fitness(neighbour)\n",
    "\n",
    "            if neighbour_fitness > fitness:\n",
    "                current = neighbour\n",
    "                fitness = neighbour_fitness\n",
    "                improved = True\n",
    "                break\n",
    "        \n",
    "    return current, fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The (1+1)MA now applies local search after each mutation step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def oneplusonema():\n",
    "    fitness_values.clear()\n",
    "    current = get_random_solution()\n",
    "    fitness = get_fitness(current)\n",
    "\n",
    "    while len(fitness_values) < max_steps:\n",
    "        candidate = mutate(current)\n",
    "        candidate, candidate_fitness = hillclimbing(candidate)\n",
    "        if candidate_fitness >= fitness:\n",
    "            fitness = candidate_fitness\n",
    "            current = candidate\n",
    "\n",
    "    return candidate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now let's see how the three algorithms compare on our hurdle problem:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "hillclimbing_restart()\n",
    "hc_values = fitness_values[:]\n",
    "\n",
    "oneplusoneea()\n",
    "opo_values = fitness_values[:]\n",
    "\n",
    "oneplusonema()\n",
    "opoma_values = fitness_values[:]\n",
    "\n",
    "plt.ylabel('Fitness', fontsize=15)\n",
    "plt.plot(hc_values, label = \"Hillclimbing\")\n",
    "plt.plot(opo_values, label = \"(1+1)EA\")\n",
    "plt.plot(opoma_values, label = \"(1+1)MA\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The local improvement actually reflects Lamarckian evolution theory, as individuals can change their own genotype. An alternative way to implement a memetic algorithm would be to exploit the _Baldwin effect_: The genotype is not improved by the local search, but the fitness is evaluated on an improved genotype. The Baldwin effect describes how the preference of the locally improvable individuals leads to this establishing in the genotype eventually through evolution. However, when implementing algorithms we are not bound by biological realisties, so we can just implement Lamarckism."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When generalising our (1+1)MA to population-based memetic algorithms we are faced with a new parameter: Applying local improvement on _all_ offspring is intuitively very computationally expensive, as we need to explore the neighbourhoods of (potentially) many individuals. The local improvement is therefor usually only done probabilistically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "P_localsearch = 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def ma():\n",
    "    fitness_values.clear()\n",
    "    population = [get_random_solution() for _ in range(population_size)]\n",
    "    best_solution = max(population, key=lambda k: get_fitness(k))\n",
    "    best_fitness = get_fitness(best_solution)\n",
    "\n",
    "    while len(fitness_values) < max_steps:\n",
    "        new_population = []\n",
    "\n",
    "        while len(new_population) < len(population):\n",
    "            parent1 = tournament_selection(population)\n",
    "            parent2 = tournament_selection(population)\n",
    "\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",
    "            if random.random() < P_localsearch:\n",
    "                offspring1, offspring1_fitness = hillclimbing(offspring1)\n",
    "                offspring2, offspring2_fitness = hillclimbing(offspring2)\n",
    "\n",
    "            new_population.append(offspring1)\n",
    "            new_population.append(offspring2)\n",
    "\n",
    "        population = new_population\n",
    "\n",
    "        best_solution = max(population, key=lambda k: get_fitness(k))\n",
    "        best_fitness = get_fitness(best_solution)\n",
    "\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As a baseline for our experiments, we will also use a random search as usual:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def randomsearch():\n",
    "    fitness_values.clear()\n",
    "    best = get_random_solution()\n",
    "    best_fitness = get_fitness(best)\n",
    "\n",
    "    while len(fitness_values) < max_steps:\n",
    "        candidate = get_random_solution()\n",
    "        fitness = get_fitness(candidate)\n",
    "        if fitness > best_fitness:\n",
    "            best = candidate\n",
    "            best_fitness = fitness\n",
    "\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ga()\n",
    "ga_values = fitness_values[:]\n",
    "\n",
    "ma()\n",
    "ma_values = fitness_values[:]\n",
    "\n",
    "randomsearch()\n",
    "random_values = fitness_values[:]\n",
    "\n",
    "plt.ylabel('Fitness', fontsize=15)\n",
    "plt.plot(hc_values, label = \"Hillclimbing\")\n",
    "plt.plot(opo_values, label = \"(1+1)EA\")\n",
    "plt.plot(opoma_values, label = \"(1+1)MA\")\n",
    "plt.plot(ga_values, label = \"GA\")\n",
    "plt.plot(ma_values, label = \"MA\")\n",
    "plt.plot(random_values, label = \"Random\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Individual runs may vary (and we will do a more systematic comparison later on), but generally the MA-variants tend to find optimal solutions, while the others do not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "While the memetic versions seem to have an edge over the non-memetic versions, the performance is dependent on the size of the neighbourhood: For example, if we use a large `n` in a bitstring representation, then the neighbourhood may become very large, and the exploration used in the local search may be very expensive. An upper bound on the number of fitness evaluations could prevent the local search from wasting the search budget in those cases where the local search does not actually contribute to improving the solution. We can investigate the effect of the neighbourhood size by re-running the search for different values of `n`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.utils import io\n",
    "\n",
    "def run_times(algorithm, repetitions):\n",
    "    global fitness_values\n",
    "    result = []\n",
    "    for i in range(repetitions):\n",
    "        with io.capture_output() as captured: \n",
    "            algorithm()\n",
    "        result.append(fitness_values[-1])\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(18, 4))\n",
    "\n",
    "num_plot = 0\n",
    "for n in [100, 200, 300, 400, 500]:\n",
    "    results = {\n",
    "        \"Standard GA\"  : run_times(ga, 10),\n",
    "        \"MA\"           : run_times(ma, 10)\n",
    "    }\n",
    "    axes[num_plot].boxplot(results.values())\n",
    "    axes[num_plot].set_title(f\"n = {n}\")\n",
    "    axes[num_plot].set_xticklabels(results.keys())\n",
    "    num_plot += 1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As a point for comparison, let's also consider the standard one-max problem for the same sizes, where we should see a bit more spread in terms of fitness values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "hurdle_fitness = get_fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(solution):\n",
    "    fitness = sum(solution)\n",
    "    \n",
    "    if not fitness_values:\n",
    "        fitness_values.append(fitness)\n",
    "    else:\n",
    "        fitness_values.append(max(fitness, fitness_values[-1]))\n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "onemax_fitness = get_fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(18, 4))\n",
    "\n",
    "num_plot = 0\n",
    "for n in [100, 200, 300, 400, 500]:\n",
    "    results = {\n",
    "        \"Standard GA\"  : run_times(ga, 10),\n",
    "        \"MA\"           : run_times(ma, 10)\n",
    "    }\n",
    "    axes[num_plot].boxplot(results.values())\n",
    "    axes[num_plot].set_title(f\"n = {n}\")\n",
    "    axes[num_plot].set_xticklabels(results.keys())\n",
    "    num_plot += 1\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Restore \n",
    "get_fitness = hurdle_fitness\n",
    "n = 100"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Hide code",
  "hide_code_all_hidden": false,
  "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
}