2_ga/genetic_algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Stochastic Optimization\n",
    "\n",
    "## 2.1 Genetic Algorithms\n",
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-sa/4.0/\"><img alt=\"Creative Commons License\" align=\"left\" src=\"https://i.creativecommons.org/l/by-sa/4.0/80x15.png\" /></a>&nbsp;| Dennis G. Wilson | <a href=\"https://d9w.github.io/evolution/\">https://d9w.github.io/evolution/</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Outline\n",
    "\n",
    "1. [Problem](#problem)\n",
    "2. [Population](#population)\n",
    "3. [Evaluation](#evaluation)\n",
    "4. [Selection](#selection)\n",
    "5. [Crossover](#crossover)\n",
    "6. [Mutation](#mutation)\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/d9w/evolution/blob/master/2_ga/genetic_algorithm.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Today we'll look at an algorithm that can be adapted to meet problem constraints and which is often used in binary or discrete optimization: the Genetic Algorithm. This algorithm uses random selection and genetic recombination in a large population of individuals. By keeping a diverse set of individuals, the Genetic Algorithm can search around multiple **different solutions** to a problem **at the same time**. It has a flexible definition which allows for adaptation to different problems, but the basic components remain the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/ga.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <a id=\"problem\"></a>Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "To understand the Genetic Algorithm, we will use it to solve the Travelling Salesman Problem. This is a well-studied problem that simulates a salesman that must visit every city in a list of cities exactly once. The goal is to minimize the total distance traveled. Plenty of applications to this problem exist, such as supply chain management, manufacturing, and navigation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "[![TSP on XKCD](https://imgs.xkcd.com/comics/travelling_salesman_problem.png)](https://xkcd.com/399/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The TSP is an NP-hard and well-studied problem. The maximum complexity of the problem can be found by assuming an exhaustive search of all possible solutions - for n cities, this is O(n!). Various methods have reduced this computational complexity, including recent algorithms for quantum computers. Genetic Algorithms offer multiple benefits for this problem type, notably scaling up to a large number of cities, being flexible to modifications to the problem such as constraints, and easy parallelization. However, there is not guarantee of an exact solution, only an approximation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "To set up our TSP scenario, we will randomly generate some cities and use their Euclidean distance as the travel cost between them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "import seaborn as sns\n",
    "\n",
    "np.random.seed(1234)\n",
    "rng = np.random.default_rng()\n",
    "sns.set_theme(style='ticks')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n_cities = 10\n",
    "cities = np.random.rand(n_cities, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "plt.scatter(cities[:, 0], cities[:, 1])\n",
    "plt.title(\"Cities\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For our problem of 10 cities, we will calculate the different computational complexities. Remember that we will compare this against the number of **evaluations** in our algorithm (number of generations times the population size)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "n = n_cities\n",
    "print(f\"O(n!): {math.factorial(n)}\")\n",
    "print(f\"O(n^2 2^n): {n**2 * 2**n}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <a id=\"population\"></a>Population"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The first decision when using a Genetic Algorithm is the choice of solution representation. To create the first population, we create random solutions (individuals) in this representation (genome), and we'll later modify individuals by combining and mutating the genome. The genome representation depends on the problem, which for this example is the TSP. To represent a solution, we can use the order of cities traveled. Each solution must have each city *exactly once*, ie a permutation of the indices of the list of cities. This is a concise way to represent the path traveled which is also easy to mutate, as we'll see later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "order0 = rng.permutation(n_cities)\n",
    "order1 = rng.permutation(n_cities)\n",
    "print(order0)\n",
    "print(order1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(14, 4))\n",
    "ax = plt.subplot(1, 2, 1)\n",
    "ind = order0\n",
    "plt.title(\"Order 0\")\n",
    "ax.scatter(cities[ind, 0], cities[ind, 1])\n",
    "ax.plot(cities[ind, 0], cities[ind, 1], 'k')\n",
    "for i in range(n_cities):\n",
    "    ax.annotate(i, (cities[ind[i], 0]+0.01, cities[ind[i], 1]+0.01), c='g')\n",
    "ax = plt.subplot(1, 2, 2)\n",
    "ind = order1\n",
    "plt.title(\"Order 1\")\n",
    "ax.scatter(cities[ind, 0], cities[ind, 1])\n",
    "ax.plot(cities[ind, 0], cities[ind, 1], 'k')\n",
    "for i in range(n_cities):\n",
    "    ax.annotate(i, (cities[ind[i], 0]+0.01, cities[ind[i], 1]+0.01), c='g')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For our initial GA population, we will create 100 random permutations. Note that there are $10!$ possible solutions to the (asymmetric) TSP with 10 cities, so 100 is very small compared to the total size. However, it is much larger than the population sizes used by CMA-ES or others. This is one distinguishing feature of Genetic Algorithms: they tend to use larger population sizes in order to represent solutions which are diverse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "n_population = 100\n",
    "population = np.array([rng.permutation(n_cities) for i in range(n_population)])\n",
    "population[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The evaluation step of a Genetic Algorithm is the same as in other stochastic algorithms: each solution must be evaluated by the objective function. In this case, that's the total distance traveled when using the solution order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "d = np.zeros((n_cities, n_cities))\n",
    "for i in range(n_cities):\n",
    "    for j in range(i):\n",
    "        d[i, j] = np.sqrt((cities[i, 0] - cities[j, 0])**2 + (cities[i,1] - cities[j, 1])**2)\n",
    "        d[j, i] = d[i, j]\n",
    "d[:3, :3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def total_distance(order, distances):\n",
    "    t = 0\n",
    "    for i in range(1, len(order)):\n",
    "        t += distances[order[i-1], order[i]]\n",
    "    return t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(\"order 0:\", total_distance(order0, d))\n",
    "print(\"order 1:\", total_distance(order1, d))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def evaluate(population, distances):\n",
    "    fitness = np.zeros(len(population))\n",
    "    for i in range(len(population)):\n",
    "        fitness[i] = total_distance(population[i], distances)\n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fitness = evaluate(population, d)\n",
    "for i in range(3):\n",
    "    print(population[i], fitness[i])\n",
    "print(\"Minimum: \")\n",
    "print(population[np.argmin(fitness)], np.min(fitness), np.max(fitness))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <a id=\"selection\"></a>Selection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Which individuals should pass on their genetic information to the next generation? We could imagine a simple schemes of taking the best individuals globally, say 20% of them. This is known as truncation selection. In simple problems with no local minima, such a method might work. However, we would lose important genetic diversity, one of the main advantages of the large population in a Genetic Algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def truncation_selection(population, fitness, p=0.2):\n",
    "    n_elites = int(np.floor(len(population) * p))\n",
    "    elites = np.argsort(fitness)[:n_elites]\n",
    "    return population[elites], fitness[elites]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "elites, efit = truncation_selection(population, fitness)\n",
    "plt.scatter(fitness, np.zeros(len(fitness)), color='k', label='all')\n",
    "plt.scatter(efit, np.ones(len(efit)), color='b', label='truncation')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "To preserve genetic diversity, we'll instead allow any member of the population to have a chance of being selected, based on their fitness. One method would be to draw from the population randomly using a probability based on the fitness of each individual. This is known as fitness proportionate selection. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def fp_selection(population, fitness):\n",
    "    p = (np.max(fitness) - fitness)\n",
    "    p /= np.sum(p)\n",
    "    rng = np.random.default_rng()\n",
    "    ind = rng.choice(len(population), p=p)\n",
    "    return population[ind], fitness[ind]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fp_fits = np.zeros(len(efit))\n",
    "for i in range(len(efit)):\n",
    "    p, f = fp_selection(population, fitness)\n",
    "    fp_fits[i] = f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.scatter(fitness, np.zeros(len(fitness)), color='k', label='all')\n",
    "plt.scatter(efit, np.ones(len(efit)), color='b', label='truncation')\n",
    "plt.scatter(fp_fits, 2*np.ones(len(fp_fits)), color='r', label='FP')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "One of the desirable properties certain evolutionary strategies, notably CMA-ES, is invariance to certain transformations in the search space. Because CMA-ES is based on the **order** of fitness values instead of their value, any transformation of the search space which maintains population order does not affect the search. This makes the search more generalizable to search spaces which are not convex or which are ill-conditioned (think bumpy).\n",
    "\n",
    "Hansen, Nikolaus, et al. \"Impacts of invariance in search: When CMA-ES and PSO face ill-conditioned and non-separable problems.\" Applied Soft Computing 11.8 (2011): 5755-5769."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In a genetic algorithm, we can acheive the same invariance by using a random ordering of individuals. To do so, we randomly select a subset of individuals, called a tournament, and then take the best individual from that random subset. This is known as tournament selection. Because tournament selection doesn't depend on the absolute fitness value but rather the ranking of individuals in a tournament, it is also invariant to order-preserving transformations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 1</h3>\n",
    "\n",
    "Complete the following tournament selection definition which uses a tournament size of `t_size=3`. Plot the selected individuals and compare it to truncation selection and fitness proportionate selection.\n",
    "</div>\n",
    "\n",
    "[solution](https://github.com/d9w/evolution/blob/master/2_ga/solutions/ex1.py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def tournament_selection(population, fitness, t_size=3):\n",
    "    ind = rng.choice(len(population))\n",
    "    return population[ind], fitness[ind]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "t_fits = np.zeros(len(efit))\n",
    "for i in range(len(efit)):\n",
    "    p, f = tournament_selection(population, fitness)\n",
    "    t_fits[i] = f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "It is difficult to represent diversity for the TSP, but below there are two different visualizations of the paths present in a population. The first maps all paths on the map, and the second plots a heatmap of the links between different city indices. A diverse population will present a variety of links, so the distribution in the heatmap should be closer to uniform for a diverse selection scheme."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 2</h3>\n",
    "\n",
    "Using the following visualizations, compare the diversity of the selected population from tournament selection to that of truncation selection. Do you notice a change in the distribution of links present?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def plot_paths(pop):\n",
    "    q=np.linspace(0.0, 1.01, len(pop))\n",
    "    cmap = plt.cm.get_cmap('jet')\n",
    "    cmaplist = [cmap(x) for x in q]\n",
    "    plt.scatter(cities[:, 0], cities[:, 1])\n",
    "    for i in range(len(pop)):\n",
    "        ind = pop[i]\n",
    "        plt.plot(cities[ind, 0], cities[ind, 1], color=cmaplist[i], alpha=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.scatter(fitness, np.zeros(len(fitness)), color='k', label='all')\n",
    "plt.scatter(efit, np.ones(len(efit)), color='b', label='truncation')\n",
    "plt.scatter(fp_fits, 2*np.ones(len(fp_fits)), color='r', label='FP')\n",
    "plt.scatter(t_fits, 3*np.ones(len(t_fits)), color='g', label='tournament')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plot_paths(population)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plot_paths(elites)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_links(population):\n",
    "    links = np.zeros((10, 10))\n",
    "    for i in population:\n",
    "        for j in range(9):\n",
    "            links[i[j], i[j+1]] += 1\n",
    "        links[i[-1], i[0]] += 1\n",
    "    links /= np.max(links)\n",
    "    return links"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.imshow(get_links(population));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.imshow(get_links(elites));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <a id=\"crossover\"></a>Crossover\n",
    "\n",
    "Considering we have such a large population, is there some way to combine individual solutions to lead to better solutions? For example, could we make an individual which inherits information from two parent individuals? This is the idea behind crossover, the other operator in genetic algorithms besides mutation. It is based on sexual reproduction where the genetic information of two parent individuals is mixed to create an offspring individual. The idea of combining the information from multiple individuals together to create the next generation is something we'll explore in more detail next class when discussing evolutionary strategies. For now, let's look at ways to combine two individuals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://github.com/d9w/evolution/raw/master/imgs/crossover.png\" width=\"80%\" height=\"auto\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def one_point(p1, p2):\n",
    "    rng = np.random.default_rng()\n",
    "    x = rng.choice(np.arange(1, np.minimum(len(p1)-1, len(p2)-1)))\n",
    "    return np.concatenate((p1[:x], p2[x:])), np.concatenate((p2[:x],p1[x:]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "parent1, _ = tournament_selection(population, fitness)\n",
    "parent2, _ = tournament_selection(population, fitness)\n",
    "print(parent1, parent2)\n",
    "child1, child2 = one_point(parent1, parent2)\n",
    "print(\"crossover: \")\n",
    "print(child1, child2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "The one point crossover method is a common one, inspired by biology and useful for combining solutions in a number of problems. When increasing the number of swapping points beyond 1, this is known as k-point crossover, and it is useful when the problem dimensions are independent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In the Travelling Salesman problem, however, the dimensions of the problem are not independent; there is the constraint that all cities must be visited and also that no city can be visted twice. Instead, we'll use a crossover operator which respects those constraints. This is known as the edge recombination operator and it can be applied to any directed graph recombination problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/e/e5/Genetic_ero_crossover.svg\" alt=\"Genetic ero crossover.svg\" width=\"50%\" height=\"auto\">\n",
    "\n",
    "https://en.wikipedia.org/wiki/Edge_recombination_operator\n",
    "    \n",
    "Whitley, L. Darrell, Timothy Starkweather, and D'Ann Fuquay. \"Scheduling problems and traveling salesmen: The genetic edge recombination operator.\" ICGA. Vol. 89. 1989. [pdf](TSP_crossover.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "This operator randomly selects from the neighborhood of each node. This operator is rather long to code, so we'll use the [definition](https://github.com/anyoptimization/pymoo/blob/master/pymoo/operators/crossover/erx.py) from the [pymoo](https://pymoo.org/index.html) library. We'll use `pymoo` in the next notebook, so install it now if you haven't yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "!pip install pymoo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from pymoo.operators.crossover import erx\n",
    "\n",
    "print(\"Parents\")\n",
    "print(population[0], population[1])\n",
    "print(\"Adjacency Matrix 1\")\n",
    "print(erx.calc_adjency_matrix(population[0]))\n",
    "print(\"Adjacency Matrix 2\")\n",
    "print(erx.calc_adjency_matrix(population[1]))\n",
    "print(\"Child\")\n",
    "erx.erx(population[0], population[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <a id=\"mutation\"></a>Mutation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Crossover combines genetic information already existing in the population. In order to search outside the possible combinations of the population, we need to directly modify individuals. This is done using mutation, which is applied either directly to individuals from selection or, more commonly, from the new individuals created by crossover. In general, mutation operators should not greatly modify an individual, changing only a small percentage of genes. An example is uniform mutation, which randomly resamples some genes. A common mutation rate for this is `1/n_genes`, meaning that, on average, only one gene is changed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "As with crossover, we need to create a mutation operator which respects the constraints of the Travelling Salesman Problem. One such mutation operator would be to switch the order of a single random pair of cities, the same \"neighbor\" generating operator we discussed for Simulated Annealing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 3</h3>\n",
    "\n",
    "Complete the following mutation operator definition which inverts a single random pair of cities.\n",
    "</div>\n",
    "\n",
    "\n",
    "[solution](https://github.com/d9w/evolution/blob/master/2_ga/solutions/ex2.py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def mutate(ind):\n",
    "    return ind"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "print(population[0], population[1])\n",
    "child = erx.erx(population[0], population[1])\n",
    "print(child)\n",
    "mutate(child)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": [
     "hide-cell"
    ]
   },
   "source": [
    "#### Solution to exercise 3\n",
    "\n",
    "```python\n",
    "def mutate(ind):\n",
    "    rng = np.random.default_rng()\n",
    "    x = rng.choice(len(ind)-1)\n",
    "    child = np.copy(ind)\n",
    "    child[x] = ind[x+1]\n",
    "    child[x+1] = ind[x]\n",
    "    return child\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## <a id=\"ga\"></a>The Genetic Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now that we have all the different parts, we can combine them in the full genetic algorithm. We'll use two different selection methods: first, a truncation selection of the best few individuals which will pass directly to the next population. This is known as **elitism** and is done to preserve the best solutions between generations. Then, we'll use tournament selection to select parents for crossover. Finally, we'll mutate the result from crossover and pass this new individual into the next population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def ga_step(population):\n",
    "    fitness = evaluate(population, d)\n",
    "    next_pop, _ = truncation_selection(population, fitness)\n",
    "    while len(next_pop) < len(population):\n",
    "        parent1, _ = tournament_selection(population, fitness)\n",
    "        parent2, _ = tournament_selection(population, fitness)\n",
    "        child = erx.erx(parent1, parent2)\n",
    "        child = mutate(child)\n",
    "        next_pop = np.concatenate((next_pop, [child]))\n",
    "    return next_pop, fitness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n_cities = 20\n",
    "cities = np.random.rand(n_cities, 2)\n",
    "d = np.zeros((n_cities, n_cities))\n",
    "for i in range(n_cities):\n",
    "    for j in range(i):\n",
    "        d[i, j] = np.sqrt((cities[i, 0] - cities[j, 0])**2 + (cities[i,1] - cities[j, 1])**2)\n",
    "        d[j, i] = d[i, j]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n_population = 100\n",
    "n_gen = 100\n",
    "population = np.array([rng.permutation(n_cities) for i in range(n_population)])\n",
    "minfit = np.zeros(n_gen)\n",
    "for i in range(n_gen):\n",
    "    population, fitness = ga_step(population)\n",
    "    minfit[i] = np.min(fitness)\n",
    "    if i > 2 and minfit[i] < minfit[i-1]:\n",
    "        print(i, minfit[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(minfit);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 3</h3>\n",
    "\n",
    "Study the impact of the crossover operator. Change the definition of the genetic algorithm to only use mutation from a single parent. How does this compare to using the ERX? What about when using a one-point crossover?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "    <h3>Exercise 4</h3>\n",
    "\n",
    "Study two hyperparameters of the genetic algorithm: the population size and the number of elites. How do these affect the search? Are elites necessary?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "<div class=\"alert alert-info\">\n",
    "    <h3>Discussion</h3>\n",
    "\n",
    "In the Genetic Algorithm, how many evaluations were performed? Is this the same as the total number of permutations explored? Roughly, how many permutations were calculated compared to the total number of possible city permutations? Has the GA converged on the best possible solution? When should the GA stop?\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}