Random and Local Search.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Random and Local Search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Random Search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Before we can code, some imports: We are considering randomised algorithms, and we want to plot results. To run this notebook, you will also need matplotlib, numpy, and scipy installed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will start with the One-Max problem: The aim of this optimisation problem is to find a binary vector of length `n` such that all elements are `1`. The solution is obvious to us, but not to a search algorithm. We first need to decide on a specific length of `n` for our experiments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "n = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to do any kind of random search, we need to be able to randomly sample solutions from the search space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_solution():\n",
    "    return [random.choice([0,1]) for _ in range(n)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now we can look at an example solution from our search space:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_random_solution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Every invocation will produce a different individual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_random_solution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To solve an optimisation problem, we also need an objective function, or _fitness function_, that quantifies how good a particular solution is with respect to the optimisation problem we are trying to solve. For our example one max problem, we can estimate how close we are to solving the problem by counting the number of `1`s in our solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(solution):\n",
    "    return sum(solution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_fitness([0,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_fitness([1,1,1,0,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will keep track of what the search does; for this I'm simply using a plain list that contains the best fitness value observed at each iteration of the search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In some cases, we know when we have found an optimal solution; sometimes we don't. Sometimes it may also take very long to find the optimal solution, therefore we limit the search with some _stopping criterion_. We will start by running the search for a fixed number of steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_steps = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The following implements our first metaheuristic: _Random Search_. It's not a particularly intelligent one, but it will be the baseline we are going to compare against from now on. If the fitness landscape is bad, this algorithm may actually be better than many other metaheuristics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def randomsearch():\n",
    "\n",
    "    best = None\n",
    "    best_fitness = -1\n",
    "\n",
    "    for step in range(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",
    "            print(f\"Iteration {step}, fitness {fitness}: {candidate}\")\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    print(f\"Solution fitness {best_fitness}: {best}\")\n",
    "    return best"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Before invoking it, we also have to reset our list of fitness values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "randomsearch()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We kept track of the best solution found throughout the search, and we will use this information to learn about the progress, and compare search algorithms. Therefore, we are going to visualise it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "If you re-run the search multiple times you will get different results, but even for small values of `n` the random search will rarely find an optimal solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Random search is a _global_ search algorithm -- it considers the entire search space. An alternative strategy is to consider the local neighbourhood of candidate solutions, resulting in _local_ search algorithms. As a naive baseline, we can change our random search to a local random search, also known as a _random walk_. A random walk consists of repeatedly picking a random neighbour for a given candidate solution, therefore we need to define a function that picks a random neighbour; that is, a solution which differs only in a single step (in our case, a bit-flip)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_neighbour(candidate):\n",
    "    pos = random.randint(0, len(candidate) - 1)\n",
    "    copy = candidate[:]\n",
    "    copy[pos] = 1 - copy[pos]\n",
    "    return copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_neighbour([0,0,0,0,0,1,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_neighbour([0,0,0,0,0,1,0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_neighbour([0,0,0,0,0,1,0,0,0,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To do a random walk, we start at a random solution, and repeatedly pick a random neighbour. Along the search, we keep track of the best solution found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def randomwalk():\n",
    "    current = get_random_solution()\n",
    "    fitness = get_fitness(current)\n",
    "\n",
    "    best = current[:]\n",
    "    best_fitness = fitness\n",
    "    print(f\"Starting at fitness {fitness}: {current}\")\n",
    "\n",
    "    for step in range(max_steps):\n",
    "        current = get_neighbour(current)\n",
    "        fitness = get_fitness(current)\n",
    "\n",
    "        if fitness > best_fitness:\n",
    "            best = current\n",
    "            best_fitness = fitness\n",
    "            print(f\"Iteration {step}, fitness {fitness}: {current}\")\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    print(f\"Solution fitness {best_fitness}: {best}\")\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "randomwalk()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A random walk usually isn't well suited as an actual search algorithm, it tends to produce even worse solutions than the random search (at least on one max). However, we can now generalize the idea of exploring the neighbourhood to construct some actual local search algoritms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Local Search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A local search algorithm consists of (1) a definition of the neighbourhood for a candidate solution, and (2) an exploration operator that picks the next step in the neighbourhood. The neighbourhood is the set of solutions that one can reach from a candidate solution in one step -- in our case, with a single bit flip:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "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": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_neighbours([0,0,0,0,0,1,0,0,0,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Hillclimbing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The first local search algorithm we consider is very simple: The exploration operator always picks the best neighbour. The result is called _hillclimbing_ because it climbs a hill in the fitness landscape (if we assume a maximisation problem, like in the case of one max). Starting from a random individual, we consider all neighbours, and pick the best one. Then, we repeat this process for this neighbour. The new point in the search space has to be strictly better."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def hillclimbing():\n",
    "\n",
    "    current = get_random_solution()\n",
    "    best = current[:]\n",
    "    best_fitness = get_fitness(current)\n",
    "    print(f\"Starting at fitness {best_fitness}: {current}\")\n",
    "\n",
    "    for step in range(max_steps):\n",
    "\n",
    "        best_neighbour = None\n",
    "        neighbour_fitness = -1\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",
    "            \n",
    "        if neighbour_fitness > best_fitness:\n",
    "            print(f\"Iteration {step}, fitness {best_fitness}: {best}\")\n",
    "            current = best_neighbour\n",
    "            best = current[:]\n",
    "            best_fitness = neighbour_fitness\n",
    "\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "\n",
    "    print(f\"Solution fitness {best_fitness}: {best}\")\n",
    "    return best"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When invoking the hillclimber, we need to reset the `fitness_values` list again so we can observe what is happening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "hillclimbing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The hillclimber will always find the solution very quickly, no matter how often you repeat the run, within at most `n` steps it will find an optimal solution. Plotting the `fitness_values` shows why the algorithm is called hillclimbing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### The n-queens problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The hillclimber seems perfect on the one max problem. However, not all problems are this easy. Let's consider an alternative problem: In the _n-queens_ problem the aim is to place _n_ queens on an _n_ x _n_ chessboard in\n",
    "such a way that they do not check each other."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As first step, we need to find a way to encode this problem such that we can apply the search on this encoding. A simple representation is a list of length _n_, where each element of the list is an integer up to _n_. That is, each element of the list maps to a column of the chess board, and the value of the element denotes the row in which a queen is located in this column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_random_solution():\n",
    "    return [random.randint(0, n - 1) for _ in range(n)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_random_solution()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to help us interpret such a vector, we can use a simple helper function that visualises the list as a matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def print_solution(solution):\n",
    "    n = len(solution)\n",
    "    for y in range(n):\n",
    "        for x in range(n):\n",
    "            if solution[x] == y:\n",
    "                print(\" Q \", end = '')\n",
    "            else:\n",
    "                print(\" . \", end = '')\n",
    "        print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "n = 8\n",
    "x = get_random_solution()\n",
    "print(x)\n",
    "print_solution(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Two queens check each other if they are in the same column, same row, or same diagonal. Since the goal is that no queens check each other, we can formulate an objective function based on how many queens check each other. To be consistent with our previous example, we will make this a _maximisation_ problem, such that the objective function calculates how many queens _don't_ check each other. Using our vector-encoding, there can only be one queen per column by construction. We thus only need to check if queens are in the same row or diagonal. For each pair of queens that does not check each other, we increase our fitness value by 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(solution):\n",
    "    fitness = 0\n",
    "    for i in range(len(solution) - 1):\n",
    "        for j in range(i + 1, len(solution)):\n",
    "            if solution[i] != solution[j] \\\n",
    "                and solution[i] != solution[j] + (j - i) \\\n",
    "                and solution[j] != solution[i] + (j - i):\n",
    "                fitness += 1\n",
    "\n",
    "    return fitness"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Here's a simple example of all queens in a row, such that they all check each other -- the fitness is thus 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x = [0,0,0,0,0,0,0,0]\n",
    "print_solution(x)\n",
    "get_fitness(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Here is a solution where no queens check each other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x = [4, 0, 3, 5, 7, 1, 6, 2]\n",
    "print_solution(x)\n",
    "get_fitness(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The optimal fitness value is n*(n-1)/2 = 28 for n = 8."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Before we can apply search on this problem, we need to define the neighbourhood of our problem. In the simplest case, a neighbour moves a queen up or down by one, taking into account the boundaries of the chess board."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def get_neighbours(solution):\n",
    "    neighbours = []\n",
    "    for i in range(len(solution)):\n",
    "        if solution[i] > 0:\n",
    "            copy = solution[:]\n",
    "            copy[i] = copy[i] - 1\n",
    "            neighbours.append(copy)\n",
    "        if solution[i] < n - 1:\n",
    "            copy = solution[:]\n",
    "            copy[i] = copy[i] + 1\n",
    "            neighbours.append(copy)\n",
    "\n",
    "    return neighbours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_neighbours([4,4,4,4,4,4,4,4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's first try how well random search copes with this problem. Since the problem is quite a bit more difficult than one max, let's also increase the maximum number of steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_steps = 1000\n",
    "fitness_values = []\n",
    "randomsearch()\n",
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Now let's compare this to our local search algorithm, hillclimbing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "hillclimbing()\n",
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Is hillclimbing still better than the random search? In many cases, it isn't. The reason is that, unlike for one max, there are _local optima_ in the search space. Our hillclimber is forced to consider only strictly better neighbours, hence it is stuck once it reaches a local optimum. (If we would allow it to move across fitness-plateaux, i.e., neighbours of equal fitness, it might get stuck in an infinite loop.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As an example, consider the following solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "x = [7, 1, 2, 0, 5, 3, 0, 4]\n",
    "print_solution(x)\n",
    "get_fitness(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's consider the fitness values of all the neighbours of this solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "for neighbour in get_neighbours([7, 1, 2, 0, 5, 3, 0, 4]):\n",
    "    print(f\"{neighbour} -> {get_fitness(neighbour)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As you can see, none of the neighbours of `x` has a better fitness value than `x` itself. The hillclimber has nowhere to go, and is stuck."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A simple solution to this problem is to use _random restarts_: If the hillclimber is stuck in an optimum, we restart the search from a new random point. With some luck, the next attempt will climb a different optimum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def hillclimbing_restart():\n",
    "\n",
    "    current = get_random_solution()\n",
    "    best = current[:]\n",
    "    best_fitness = get_fitness(current)\n",
    "\n",
    "    for step in range(max_steps):\n",
    "\n",
    "        best_neighbour = None\n",
    "        neighbour_fitness = -1\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",
    "        fitness_values.append(best_fitness)        \n",
    "\n",
    "\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "hillclimbing_restart()\n",
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Using the random restarts, the hillclimber tends to be better than before, and better than random search -- even if it still doesn't always find the best solution, depending on what value you choose for `max_steps`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The version we considered is called _steepest ascent_ hillclimbing, because it will always explore all neighbours, and then pick the best one. An alternative approach is _first ascent_, which will pick the first neighbour that has a better fitness value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def hillclimbing_first_ascent():\n",
    "\n",
    "    current = get_random_solution()\n",
    "    fitness = get_fitness(current)\n",
    "    best = current[:]\n",
    "    best_fitness = fitness\n",
    "\n",
    "    for step in range(max_steps):\n",
    "\n",
    "        replaced = False\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",
    "                replaced = True\n",
    "                break\n",
    "\n",
    "        # Random restart if no neighbour is better\n",
    "        if not replaced:\n",
    "            current = get_random_solution()\n",
    "            fitness = get_fitness(current)\n",
    "            \n",
    "        if fitness > best_fitness:\n",
    "            best = current[:]\n",
    "            best_fitness = fitness\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    return best"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For some problems, a first-ascent strategy may be less likely to run into local optima, because it is not so \"eager\" to climb the hill. However, this depends on the problem at hand, and luck. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "hillclimbing_first_ascent()\n",
    "plt.plot(fitness_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For our n-queens problem, is steepest ascent or first ascent better? And is either of them better than random search? To find out, we can run all three algorithms and compare the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "randomsearch()\n",
    "f_random = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "hillclimbing_restart()\n",
    "f_hill = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "hillclimbing_first_ascent()\n",
    "f_first = fitness_values[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(f_random, label=\"Random Search\")\n",
    "plt.plot(f_hill, label=\"Hill Climbing (Steepest Ascent)\")\n",
    "plt.plot(f_first, label=\"Hill Climbing (First Ascent)\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There are two problems with this comparison:\n",
    "\n",
    "1. Every time you run the algorithm, you will get a different result. Sometimes one algorithm performs better than another, sometimes it doesn't. There may be trends, though.\n",
    "2. The comparison actually isn't fair. In each iteration, random search considers exactly one candidate solution, steepest ascent considers up to _2*n_ solutions, and first ascent is somewhere in between. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To address the first problem, we need to compare multiple runs, not just individual runs. We will get back to this problem in a bit.\n",
    "\n",
    "To make the comparison fair, we need to use a stopping condition that amounts to the same computation time for each algorithm. In practice, one often uses _time_ (because users may not understand the algorithms applied but do have an opinion on how long they are prepared to wait). Another metric for fair comparisons, which we will use, is the number of _fitness evaluations_. This is already what our random search does, but we have to adapt our hillclimbers slightly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def hillclimbing_restart():\n",
    "    current = get_random_solution()\n",
    "    best, best_fitness = current[:], get_fitness(current)\n",
    "    step = 1\n",
    "    while step < max_steps:\n",
    "        best_neighbour, neighbour_fitness = None, -1\n",
    "        for neighbour in get_neighbours(current):\n",
    "            fitness = get_fitness(neighbour)\n",
    "            step += 1\n",
    "\n",
    "            if fitness > neighbour_fitness:\n",
    "                best_neighbour = neighbour\n",
    "                neighbour_fitness = fitness\n",
    "                if fitness > best_fitness:\n",
    "                    best = current[:]\n",
    "                    best_fitness = fitness\n",
    "            fitness_values.append(best_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",
    "            step += 1\n",
    "            if fitness > best_fitness:\n",
    "                best = current[:]\n",
    "                best_fitness = fitness\n",
    "                fitness_values.append(best_fitness)\n",
    "        else:\n",
    "            current = best_neighbour\n",
    "\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def hillclimbing_first_ascent():\n",
    "    current = get_random_solution()\n",
    "    fitness = get_fitness(current)\n",
    "    best, best_fitness = current[:], fitness\n",
    "    step = 1\n",
    "    while step < max_steps:\n",
    "        replaced = False\n",
    "        for neighbour in get_neighbours(current):\n",
    "            neighbour_fitness = get_fitness(neighbour)\n",
    "            step += 1\n",
    "            if neighbour_fitness > fitness:\n",
    "                current, fitness = neighbour, neighbour_fitness\n",
    "                replaced = True\n",
    "                if fitness > best_fitness:\n",
    "                    best = current[:]\n",
    "                    best_fitness = fitness\n",
    "                fitness_values.append(best_fitness)\n",
    "                break\n",
    "            else:\n",
    "                fitness_values.append(best_fitness)\n",
    "\n",
    "        # Random restart if no neighbour is better\n",
    "        if not replaced:\n",
    "            current = get_random_solution()\n",
    "            fitness = get_fitness(current)\n",
    "            step += 1\n",
    "            if fitness > best_fitness:\n",
    "                best, best_fitness = current[:], fitness\n",
    "                fitness_values.append(best_fitness)\n",
    "\n",
    "    return best"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []\n",
    "randomsearch()\n",
    "f_random = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "hillclimbing_restart()\n",
    "f_hill = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "hillclimbing_first_ascent()\n",
    "f_first = fitness_values[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(f_random, label=\"Random Search\")\n",
    "plt.plot(f_hill, label=\"Hill Climbing (Steepest Ascent)\")\n",
    "plt.plot(f_first, label=\"Hill Climbing (First Ascent)\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There probably isn't a clear winner yet. We will first introduce an algorithm that can cope better with local optima before turning to the question how to compare algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tabu Search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In tabu search, a tabu list prevents re-exploration of already visited parts of search space, but therefore the search is allowed to make moves in the search space that lead to worse fitness (whereas hillclimbing only allows moves that improve fitness)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The tabu list is often implemented as a circular buffer, where solutions are added as their fitness is evaluated, and they are removed once they have been in the buffer for a certain number of steps. Let's introduce a parameter to denote the size of our buffer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "tabu_size = 500"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The algorithm itself is very similar to hillclimbing, except for the following changes:\n",
    "- Only neighbours that are not in the tabu list are considered\n",
    "- The search moves to the best neighbour, even if that neighbour has worse fitness\n",
    "- The tabu list is updated at each step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def tabusearch():\n",
    "    current = get_random_solution()\n",
    "    fitness = get_fitness(current)\n",
    "    best_solution, best_fitness = current, fitness\n",
    "\n",
    "    tabu = [current]\n",
    "    step = 1\n",
    "    while step < max_steps:\n",
    "        \n",
    "        best_neighbour = None\n",
    "        neighbour_fitness = -1\n",
    "        for neighbour in get_neighbours(current):\n",
    "            if neighbour not in tabu:\n",
    "                new_fitness = get_fitness(neighbour)\n",
    "                step += 1\n",
    "                if new_fitness > neighbour_fitness:\n",
    "                    best_neighbour, neighbour_fitness = neighbour, new_fitness\n",
    "                    if neighbour_fitness > best_fitness:\n",
    "                        best_solution, best_fitness = best_neighbour, neighbour_fitness\n",
    "                fitness_values.append(best_fitness)                \n",
    "\n",
    "        # Append at the end of the tabu list\n",
    "        tabu.append(best_neighbour)\n",
    "        current = best_neighbour\n",
    "        \n",
    "        # Remove elements from the front of the tabu list\n",
    "        if len(tabu) > tabu_size:\n",
    "            tabu.pop(0)\n",
    "\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_steps = 5000\n",
    "\n",
    "fitness_values = []\n",
    "randomsearch()\n",
    "f_random = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "tabusearch()\n",
    "f_tabu = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "hillclimbing_first_ascent()\n",
    "f_first = fitness_values[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(f_random, label=\"Random Search\")\n",
    "plt.plot(f_tabu, label=\"Tabu Search\")\n",
    "plt.plot(f_first, label=\"Hill Climbing (First Ascent)\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Comparing Search Algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to draw a conclusion whether one randomised algorithm is better than another one, we cannot use individual runs, but need to run each algorithm multiple times, and then statistically check whether the observations result from different distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As a simple example, we will compare the algorithms with respect to the fitness values achieved. Does tabu search lead to better fitness values than hillclimbing? For this, we need to run an experiment in which we run each algorithm multiple times, and collect the fitness values we want to compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.utils import io\n",
    "\n",
    "max_steps = 1000\n",
    "repetitions = 30\n",
    "\n",
    "values_random = []\n",
    "values_hill = []\n",
    "values_tabu = []\n",
    "\n",
    "for i in range(repetitions):\n",
    "\n",
    "    print(f\"Repetition {i}\")\n",
    "    with io.capture_output() as captured:    \n",
    "        fitness_values = []\n",
    "        randomsearch()\n",
    "        values_random.append(fitness_values[-1])\n",
    "\n",
    "        fitness_values = []\n",
    "        hillclimbing_first_ascent()\n",
    "        values_hill.append(fitness_values[-1])\n",
    "        \n",
    "        fitness_values = []\n",
    "        tabusearch()\n",
    "        values_tabu.append(fitness_values[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(values_random)\n",
    "print(values_hill)\n",
    "print(values_tabu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As a result of this experiment, we now have three lists of fitness values, one for the results achieved by random search, one for those of our hillclimber, and one for tabu search. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Is our hillclimber any better than random search? To answer this question we compare the two corresponding datasets. First, we'll just look at the means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from numpy import mean\n",
    "mean(values_random)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mean(values_hill)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There may be some differences in the means, and by chance either of the two algorithms can have a higher mean. To really find out if there is a difference, we apply a statistical test. We will not make any assumptions on normality of data, and therefore apply a non-parametric test."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.stats import mannwhitneyu\n",
    "result = mannwhitneyu(values_random, values_hill)\n",
    "result[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "By convention, a p-value lower than 0.05 is usually considered statistically significant (though this is debated). If the p-value is larger, then either there really is no difference, or the difference is so small that more observations (i.e., more repetitions in our experiment) are required in order to detect it. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As a further way to quantify the difference, one typically uses effect size measurements. Again there are different statistical measures, in particular differing whether they are standardised (i.e., independent of the measurement unit) or non-parametric (i.e., make no assumptions on the underlying distribution). An example non-parametric effect size is the Vargha-Delaney A12 effect size measurement: Given a performance measure `M` seen in `m` measures of `X` and `n` measures of `Y`, the A12 statistics measures the probability that running algorithm `X` yields higher `M` values than running another algorithm `Y`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# https://gist.github.com/timm/5630491\n",
    "def a12(lst1,lst2,rev=True):\n",
    "  \"how often is x in lst1 more than y in lst2?\"\n",
    "  more = same = 0.0\n",
    "  for x in lst1:\n",
    "    for y in lst2:\n",
    "      if   x==y : same += 1\n",
    "      elif rev     and x > y : more += 1\n",
    "      elif not rev and x < y : more += 1\n",
    "  return (more + 0.5*same)  / (len(lst1)*len(lst2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a12(values_random, values_hill)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "For the comparison (random search, hillclimbing) we most likely see an A12 value < 0.5. The interpretation is that, if you run random search on our n-queens problem, with a probability of _A12_ it will produce a higher fitness value. If A12 < 0.5, then that means the hill climber is better. If A12 = 0, then random would _always_ produce worse results, and if A12 = 1, then hillclimbing would _always_ produce worse (lower) results; an A12 of 0.5 would mean both algorithms perform equally well. The value of A12 gives us an intuition of the magnitude, and a common interpretation is that a big effect is with A12 over 0.71, medium over 0.64, and small over 0.56 (and the same mirrored at 0.5)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We can now also check if tabu search is better than hillclimbing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = mannwhitneyu(values_hill, values_tabu)\n",
    "result[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mean(values_hill)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "mean(values_tabu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "a12(values_hill, values_tabu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The difference between hillclimbing and tabu search on our n-queens problems tends to be much smaller than the difference between these algorithms and random search, which tabu search usually performing _slightly_ better."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Simulated Annealing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The last local search algorithm we will consider is simulated annealing, which is another variant of hillclimbing with a twist in order to escape local optima. The algorithm simulates physical processes from metallurgy, where the atoms in heated metal undergo disordered movements of large amplitude, while cooling the metal down progressively reduces movement and stabilises atoms around fixed positions in a regular crystal structure with minimal energy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to simulate these effects, we will require some helper maths functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The central aspect of simulated annealing is that the local search accepts movements to neighbours with worse fitness values with a probability that is proportional to the current temperature (just like atoms move more freely in heated metal). To calculate the acceptance probability, we need to compare the current individual of the search and a neighbour. If the neighbour is better, the probability of accepting it is always `1`. If it is worse, then the probability is calculated as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def acceptance_probability(fitness, new_fitness, temperature):\n",
    "    if new_fitness > fitness:\n",
    "        return 1\n",
    "    else:\n",
    "        p = math.exp( (new_fitness - fitness) / temperature)\n",
    "        return p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The temperature is varied along the search, starting at a high value, and gradually cooling down. It is possible to use different cooling schedules; we will use a simple schedule where the temperature is calculated based on the fraction of the overall search budget (i.e., number of steps in the search) used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def temperature(fraction):\n",
    "    return max(0.01, min(1, 1 - fraction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Unlike hillclimbing and tabu search, in each step we only consider one random neighbour. The acceptance probability then decides whether the search moves to this individual, or remains at the current position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def simulatedannealing():\n",
    "    current = get_random_solution()\n",
    "    fitness = get_fitness(current)\n",
    "\n",
    "    best_solution = current\n",
    "    best_fitness = fitness\n",
    "\n",
    "    step = 1\n",
    "\n",
    "    while step < max_steps:\n",
    "        fraction = step / float(max_steps)\n",
    "        T = temperature(fraction)\n",
    "\n",
    "        neighbour = random.choice(get_neighbours(current))\n",
    "        neighbour_fitness = get_fitness(neighbour)\n",
    "        step += 1\n",
    "\n",
    "        if acceptance_probability(fitness, neighbour_fitness, T) > random.random():\n",
    "            current, fitness = neighbour, neighbour_fitness\n",
    "\n",
    "            if fitness > best_fitness:\n",
    "                best_fitness = fitness\n",
    "                best_solution = current\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In order to evaluate simulated annealing, let's first have a look at an individual run of the algorithm compared to individual runs of other algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "max_steps = 5000\n",
    "\n",
    "fitness_values = []\n",
    "randomsearch()\n",
    "f_random = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "tabusearch()\n",
    "f_tabu = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "hillclimbing_first_ascent()\n",
    "f_first = fitness_values[:]\n",
    "\n",
    "fitness_values = []\n",
    "simulatedannealing()\n",
    "f_anneal = fitness_values[:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(f_random, label=\"Random Search\")\n",
    "plt.plot(f_tabu, label=\"Tabu Search\")\n",
    "plt.plot(f_first, label=\"Hill Climbing\")\n",
    "plt.plot(f_anneal, label=\"Simulated Annealing\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To statistically compare the algorithms with each other, we need to run an experiment again for fair comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.utils import io\n",
    "\n",
    "max_steps = 10000\n",
    "repetitions = 30\n",
    "\n",
    "values_random = []\n",
    "values_hill = []\n",
    "values_tabu = []\n",
    "values_anneal = []\n",
    "\n",
    "for i in range(repetitions):\n",
    "\n",
    "    print(f\"Repetition {i}\")\n",
    "    with io.capture_output() as captured:    \n",
    "        fitness_values = []\n",
    "        randomsearch()\n",
    "        values_random.append(fitness_values[-1])\n",
    "\n",
    "        fitness_values = []\n",
    "        hillclimbing_first_ascent()\n",
    "        values_hill.append(fitness_values[-1])\n",
    "        \n",
    "        fitness_values = []\n",
    "        tabusearch()\n",
    "        values_tabu.append(fitness_values[-1])\n",
    "        \n",
    "        fitness_values = []\n",
    "        simulatedannealing()\n",
    "        values_anneal.append(fitness_values[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "print(values_random)\n",
    "print(values_hill)\n",
    "print(values_tabu)\n",
    "print(values_anneal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Is simulated annealing better than random search?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "result = mannwhitneyu(values_random, values_anneal)\n",
    "p = result[1]\n",
    "a = a12(values_random, values_anneal)\n",
    "print(\"Random search vs. Simulated Annealing: p={:f}, A12={:f}\".format(p, a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Is simulated annealing better than hillclimbing?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "result = mannwhitneyu(values_hill, values_anneal)\n",
    "p = result[1]\n",
    "a = a12(values_hill, values_anneal)\n",
    "print(\"Hill climbing vs. Simulated Annealing: p={:f}, A12={:f}\".format(p, a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Is simulated annealing better than tabu search?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "result = mannwhitneyu(values_tabu, values_anneal)\n",
    "p = result[1]\n",
    "a = a12(values_tabu, values_anneal)\n",
    "print(\"Tabu search vs. Simulated Annealing: p={:f}, A12={:f}\".format(p, a))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Most likely, simulated annealing is better than random search and hillclimbing, but there probably is no significant difference between tabu search and simulated annealing, on the n-queens problem. However, it generally is not known for a new problem which algorithm will work best."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}