Swarm Optimisation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Swarm Optimisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In this chapter we will consider search algorithms that fit into the general framework of evolutionary algorithms, but aim to imitate optimisation processes based on the swarm behaviour of different types of animals, such as flocks of birds or ant colonies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import random\n",
    "import math\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "\n",
    "# For presenting as slides\n",
    "#plt.rcParams['figure.figsize'] = [12, 8]\n",
    "#plt.rcParams.update({'font.size': 22})\n",
    "#plt.rcParams['lines.linewidth'] = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will continue to consider the one max problem we are now very familiar with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "n = 300\n",
    "\n",
    "def get_random_solution():\n",
    "    return [random.choice([0,1]) for _ in range(n)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def get_fitness(solution):\n",
    "    return sum(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will compare some new algorithms against baseline algorithms we are now familiar with, and so we need to set up the usual configuration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "max_steps = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "configuration = {\n",
    "    \"P_xover\": 0.7,\n",
    "    \"population_size\": 20,\n",
    "    \"tournament_size\": 2\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "fitness_values = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The first baseline and sanity check as always is random search."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def randomsearch():\n",
    "    fitness_values.clear()\n",
    "\n",
    "    best = None\n",
    "    best_fitness = -sys.maxsize\n",
    "\n",
    "    for step in range(max_steps):\n",
    "        # To make plots comparable between population based algorithms\n",
    "        candidates = [get_random_solution() for _ in range(configuration[\"population_size\"])]\n",
    "        candidate = max(candidates, key=lambda k: get_fitness(k))\n",
    "        fitness = get_fitness(candidate)\n",
    "        if fitness > best_fitness:\n",
    "            best = candidate\n",
    "            best_fitness = fitness\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    return best"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "As second baseline, we will use the familiar standard genetic algorithm."
   ]
  },
  {
   "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 tournament_selection(population):\n",
    "    # Make sure the sample isn't larger than the population\n",
    "    candidates = random.sample(population, min(len(population), configuration[\"tournament_size\"]))\n",
    "    winner = max(candidates, key=lambda x: get_fitness(x))    \n",
    "                \n",
    "    return winner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def crossover(parent1, parent2):\n",
    "    pos = random.randint(0, len(parent1))\n",
    "    offspring1 = parent1[:pos] + parent2[pos:]\n",
    "    offspring2 = parent2[:pos] + parent1[pos:]\n",
    "    return offspring1, offspring2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def ga():\n",
    "    fitness_values.clear()\n",
    "    population = [get_random_solution() for _ in range(configuration[\"population_size\"])]\n",
    "    best_solution = max(population, key=lambda k: get_fitness(k))\n",
    "    best_fitness = get_fitness(best_solution)\n",
    "    iteration = 0\n",
    "    \n",
    "    while iteration < max_steps:\n",
    "        new_population = []\n",
    "        iteration += 1\n",
    "\n",
    "        while len(new_population) < len(population):\n",
    "            parent1 = tournament_selection(population)\n",
    "            parent2 = tournament_selection(population)\n",
    "\n",
    "            if random.random() < configuration[\"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",
    "            new_population.append(offspring1)\n",
    "            new_population.append(offspring2)\n",
    "\n",
    "        population.clear()\n",
    "        population.extend(new_population)\n",
    "\n",
    "        best_solution = max(population + [best_solution], key=lambda k: get_fitness(k))\n",
    "        best_fitness = get_fitness(best_solution)\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's recall how these algorithms compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "ga()\n",
    "ga_values = fitness_values[:]\n",
    "\n",
    "randomsearch()\n",
    "random_values = fitness_values[:]\n",
    "\n",
    "plt.ylabel('Fitness', fontsize=15)\n",
    "plt.plot(ga_values, label = \"GA\")\n",
    "plt.plot(random_values, label = \"Random\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Ant Colony Optimisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Ants walking to and from a food source deposit pheromone on the ground. Other ants perceive the presence of pheromone and tend to follow paths where the pheromone concentration is higher.  Through this mechanism, ants are able to transport food to their nest in a remarkably effective way. The pheromones represent a form of stigmergic communication, i.e., an indirect, non-symbolic form of communication mediated by the environment. This stigmergic information is local, and can only be accessed by those ants that visit the locus in which it was released."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Ant colony optimisation algorithms imitate this process of stigmergic communication by allocating the edges of a graph with pheromones. An ant then represents a randomised traversal of the graph, influenced by the pheromone information. The general scheme of ant algorithms is as follows:\n",
    "\n",
    "- Step 1: Initialise the pheromone information\n",
    "- Step 2: For each ant, do the following:\n",
    "   -  Find a solution (a path) based on the current pheromone trail\n",
    "   -  Reinforcement: Add pheromone\n",
    "   -  Evaporation: Reduce pheromone\n",
    "- Step 3: Stop if terminating condition satisfied; return to Step 2 otherwise "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "There are several different popular flavours of ant colony optimisation algorithms, which mainly differ on how pheromones are reposited by ants and how the pheromones evaporate.\n",
    "\n",
    "In an _Ant System_ (AS) pheromone values are updated by all the ants. In an _Ant Colony System_ (ACS) ants update edges using a local pheromone updating rule, and only the best ant of an iteration is allowed to update the trails by applying a modified global pheromone updating rule. Finally, _Max-Min Ant Systems_ (MMAS) explicitly limit the minimum and maximum values of the pheromone, and only one ant per iteration updates the pheromones (best or iteration-best)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A classical problem to solve with ant colony optimisation is the travelling salesman problem. To solve combinatorial problems rather than graph problems, we require a _construction graph_ on which the ants can deposit pheromones, and which is then used to construct a solution based on these pheromones."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will create a construction graph for bitvectors. There is a node for each bit, with two outgoing alternative edges which will determine whether the bit is set to 0 or 1 by the ant. These two alternative edges lead to intermediate nodes, which both are then connected to node representing the successive bit. The initial pheromones are set to 1/2, which represents a uniform distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def construct_initial_graph():\n",
    "    num_nodes = 3 * n\n",
    "    construction_graph = {}\n",
    "    pheromones = {}\n",
    "    initial_probability = 1/2\n",
    "    for i in range(0, num_nodes, 3):\n",
    "\n",
    "        construction_graph[i]     = [i + 1, i + 2]\n",
    "        construction_graph[i + 1] = [i + 3]\n",
    "        construction_graph[i + 2] = [i + 3]\n",
    "\n",
    "        pheromones[(i, i + 1)]     = initial_probability\n",
    "        pheromones[(i, i + 2)]     = initial_probability\n",
    "        pheromones[(i + 1, i + 3)] = initial_probability\n",
    "        pheromones[(i + 2, i + 3)] = initial_probability\n",
    "\n",
    "    return construction_graph, pheromones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "This can be visualised using GraphViz as usual."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from graphviz import Digraph\n",
    "def plot_graph(construction_graph, pheromones):\n",
    "    num_cluster = 1\n",
    "\n",
    "    dot = Digraph()\n",
    "    for current_node in range(n-1):\n",
    "        for successor in construction_graph[current_node]:\n",
    "            dot.edge(str(current_node), str(successor), label=str(pheromones[(current_node, successor)]))\n",
    "    return dot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "graph, pheromones  = construct_initial_graph()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plot_graph(graph, pheromones)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 300"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To create a solution based on a construction graph, an ant traverses the graph starting from the initial node. For each node it chooses an edge proportional to the pheromones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_ant(construction_graph, pheromones):\n",
    "    ant = []\n",
    "    current_node = 0\n",
    "    while current_node in construction_graph:\n",
    "        successors = construction_graph[current_node]\n",
    "        successor_node = random.choices(successors, [pheromones[(current_node, successor)] for successor in successors])[0]\n",
    "        ant.append((current_node, successor_node))\n",
    "        current_node = successor_node\n",
    "\n",
    "    return ant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "ant = get_ant(graph, pheromones)\n",
    "ant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "To convert a graph traversal to a bitstring for onemax, we iterate over the edges contained in the path traversed by the ant, and check whether the edges represent 1 or 0 in our construction graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def get_bitstring(ant):\n",
    "    bitstring = []\n",
    "\n",
    "    for (from_node, to_node) in ant:\n",
    "        if from_node % 3 == 0:\n",
    "            if to_node == from_node + 1:\n",
    "                bitstring.append(1)\n",
    "            else:\n",
    "                bitstring.append(0)\n",
    "\n",
    "    return bitstring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "get_bitstring(ant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Ants that traverse an edge lead to an increase in pheromones, but pheromones also evaporate (decrease) at a certain, parameterised, rate $\\rho$.  \n",
    "\n",
    "Pheromones are evaporated with the rate $(1 − \\rho)$, although we need to make sure that pheromones do not evaporate completely, so we ensure a minimum value of $1/n$.\n",
    "\n",
    "To deposit pheromones for an ant in our bitstring construction graph, we increase the pheromones for traversed edges by $\\rho$, but again ensuring that the pheromones do not exceed $1 - 1/n$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "configuration[\"rho\"] = 0.2\n",
    "def update_pheromones(pheromones, ant):\n",
    "    updated_pheromones = {}\n",
    "\n",
    "    # Evaporation\n",
    "    for transition in pheromones.keys():\n",
    "        updated_pheromones[transition] = max((1 - configuration[\"rho\"]) * pheromones[transition], 1 / n)\n",
    "\n",
    "    # Update based on best ant\n",
    "    for transition in ant:\n",
    "        updated_pheromones[transition] = min(updated_pheromones[transition] + configuration[\"rho\"], 1 - 1 / n)\n",
    "\n",
    "    return updated_pheromones"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We consider a Max-Min Ant System (MMAS), which starts from a construction graph with all pheromones set to 1/2. For each iteration it creates an ant that traverses the construction graph, derives the bitstring represented by the path, and calculates the fitness of that bitstring. The pheromones are then updated by the overall best ant observed so far."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def mmas():\n",
    "    fitness_values.clear()\n",
    "    construction_graph, pheromones = construct_initial_graph()\n",
    "\n",
    "    best_ant = get_ant(construction_graph, pheromones)\n",
    "    best_fitness = get_fitness(get_bitstring(best_ant))\n",
    "\n",
    "    pheromones = update_pheromones(pheromones, best_ant)\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < max_steps * configuration[\"population_size\"]:\n",
    "        iteration += 1\n",
    "        if iteration % configuration[\"population_size\"] == 0:\n",
    "            fitness_values.append(best_fitness)\n",
    "\n",
    "        current_ant = get_ant(construction_graph, pheromones)\n",
    "        current_fitness = get_fitness(get_bitstring(current_ant))\n",
    "\n",
    "        if current_fitness > best_fitness:\n",
    "            best_ant = current_ant\n",
    "            best_fitness = current_fitness\n",
    "\n",
    "        pheromones = update_pheromones(pheromones, best_ant)\n",
    "\n",
    "    fitness_values.append(best_fitness)\n",
    "\n",
    "    return get_bitstring(best_ant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Our `mmas` implementation uses a pheromone update strategy where pheromones are updated only by the globally best ant observed so far. An alternative is to generate λ ants, and then to update pheromones by best of the λ ants. To keep things simple, we will set λ to our population size parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def mmas_iteration():\n",
    "    fitness_values.clear()\n",
    "    construction_graph, pheromones = construct_initial_graph()\n",
    "\n",
    "    best_ant = get_ant(construction_graph, pheromones)\n",
    "    best_fitness = get_fitness(get_bitstring(best_ant))\n",
    "\n",
    "    pheromones = update_pheromones(pheromones, best_ant)\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < max_steps:\n",
    "        iteration += 1\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "        iteration_ant = None\n",
    "        iteration_fitness = -sys.maxsize\n",
    "\n",
    "        for i in range(configuration[\"population_size\"]):\n",
    "            current_ant = get_ant(construction_graph, pheromones)\n",
    "            current_fitness = get_fitness(get_bitstring(current_ant))\n",
    "\n",
    "            if current_fitness > iteration_fitness:\n",
    "                iteration_ant = current_ant\n",
    "                iteration_fitness = current_fitness\n",
    "\n",
    "            if current_fitness > best_fitness:\n",
    "                best_ant = current_ant\n",
    "                best_fitness = current_fitness\n",
    "\n",
    "        pheromones = update_pheromones(pheromones, iteration_ant)\n",
    "\n",
    "    fitness_values.append(best_fitness)\n",
    "\n",
    "    return get_bitstring(best_ant)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Let's have a look how our Max-Min Ant Systems compare to our baseline algorithms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "randomsearch()\n",
    "random_values = fitness_values[:]\n",
    "\n",
    "ga()\n",
    "ga_values = fitness_values[:]\n",
    "\n",
    "mmas()\n",
    "mmas_values = fitness_values[:]\n",
    "\n",
    "mmas_iteration()\n",
    "immas_values = fitness_values[:]\n",
    "\n",
    "plt.ylabel('Fitness', fontsize=15)\n",
    "plt.plot(random_values, label = \"Random\")\n",
    "plt.plot(ga_values, label = \"GA\")\n",
    "plt.plot(mmas_values, label = \"MMAS*\")\n",
    "plt.plot(immas_values, label = \"MMAS-Iteration\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Particle Swarm Optimisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Particle swarm optimisation (PSO) is loosely based on a flock of birds or school of fish searching for food. Within a flock, the birds communicate by crying out. The intensity of the cry of a bird is proportional to how much food (insects) it can find at its current location. We also assume that birds know which other birds are nearby, and they can recognise who is crying out the loudest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The flock has a good chance of converging to the location with the most food, if each bird follows the direction which is the combination of (1) keeping the current direction; (2) returning to the location where it found the most food; and (3) moving towards the neighbouring bird whose cry is the loudest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "In PSO, the population is called swarm, and each individual is called a particle. To represent the three types of information that influence the movement direction, for each particle we record the current position $x$, the best position $p_b$ found so far, and the velocity $v$. In addition, we record the best particle (the global leader $p_g$) of the whole swarm found so far. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The algorithm itself is reminiscent of an evolutionary algorithm. The swarm is initialised randomly, and all particles are evaluated. Instead of the evolutionary operators we are familiar with from genetic algorithms, particles update their velocity and position before their fitness is evaluated. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We will start by considering a numerical problem where we can visualise the position and velocity, before moving on the one max. Our example problem consists simply of picking a target point in two dimensional space, and then the fitness of a particles is simply the distance to that point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "target = [1000, 1000]\n",
    "def get_distance(particle):\n",
    "    return math.sqrt((particle[0] - target[0]) ** 2 + (particle[1] - target[1]) ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We need some helper functions to visualise the swarm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.animation as animation\n",
    "from IPython.display import HTML\n",
    "\n",
    "ims = []  # global variable to store images of the animation\n",
    "    \n",
    "def initialise_plot():\n",
    "    global ims\n",
    "    global fig\n",
    "    global ax\n",
    "    \n",
    "    ims = []\n",
    "\n",
    "    %matplotlib agg\n",
    "    fig, ax = plt.subplots()\n",
    "    ims = []\n",
    "    %matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A plot consists of a number of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def addPoint(scat, new_point, c='k'):\n",
    "    old_off = scat.get_offsets()\n",
    "    new_off = numpy.concatenate([old_off,numpy.array(new_point, ndmin=2)])\n",
    "    old_c = scat.get_facecolors()\n",
    "    new_c = numpy.concatenate([old_c, numpy.array(matplotlib.colors.to_rgba(c), ndmin=2)])\n",
    "\n",
    "    scat.set_offsets(new_off)\n",
    "    scat.set_facecolors(new_c)\n",
    "\n",
    "    scat.axes.figure.canvas.draw_idle()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def plot(population):\n",
    "    function1_values = [x[0] for (x, xb, v) in population]\n",
    "    function2_values = [x[1] for (x, xb, v) in population]\n",
    "\n",
    "    x = ax.scatter(function1_values, function2_values, color=\"blue\")\n",
    "    addPoint(x, target, c='red')\n",
    "    ims.append((x,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "The velocity of a particle in the k-th iteration is updated as follows: \n",
    "\n",
    "Where k is the iteration number, 𝑎 is the inertia weight, 𝑏 and 𝑐 are the learning factors called personal influence and social influence, respectively, and w1 and w2 are random numbers taken from [0,1]. \n",
    "\n",
    "Each particle tries to learn from the current leader as well as the best position found by itself. \n",
    "\n",
    "Also depends on the random factors w1 and w2. Thus, we will obtain a different solution in a different run. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "Based on the new velocity, the new position is obtained as follows:  x(k+1) = x(k) + v(k+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "configuration[\"c1\"] = 2\n",
    "configuration[\"c2\"] = 2\n",
    "configuration[\"vmax\"] = 200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "def demo_pso():\n",
    "    initialise_plot()\n",
    "    particles = []\n",
    "    for i in range(configuration[\"population_size\"]):\n",
    "        x  = [random.randint(-10000, 10000) for _ in range(2)]\n",
    "        xb = x[:]\n",
    "        v  = [random.randint(-200, 200) for _ in range(2)]\n",
    "        particles.append((x, xb, v))\n",
    "\n",
    "    best_solution = min(particles, key=lambda x: get_distance(x[0]))[0]\n",
    "    best_fitness = get_distance(best_solution)\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < max_steps:\n",
    "        #print(f\"Iteration {iteration}: Fitness {best_fitness}: {best_solution}\")\n",
    "        iteration += 1\n",
    "\n",
    "        r1 = [random.uniform(0, configuration[\"c1\"]) for _ in range(n)]\n",
    "        r2 = [random.uniform(0, configuration[\"c2\"]) for _ in range(n)]\n",
    "\n",
    "        for (x, xb, v) in particles:\n",
    "            for i in range(2):\n",
    "                v[i] = v[i] + r1[i] * numpy.subtract(xb[i], x[i]) + r2[i] * numpy.subtract(best_solution[i], x[i])\n",
    "                v[i] = max(v[i], -configuration[\"vmax\"])\n",
    "                v[i] = min(v[i], configuration[\"vmax\"])\n",
    "                x[i] = numpy.add(x[i], v[i])\n",
    "\n",
    "            fitness_current = get_distance(x)\n",
    "            fitness_best    = get_distance(xb)\n",
    "\n",
    "            if fitness_current < fitness_best:\n",
    "                for i in range(2):\n",
    "                    xb[i] = x[i]\n",
    "                fitness_best = fitness_current\n",
    "\n",
    "            if fitness_best < best_fitness:\n",
    "                best_fitness = fitness_best\n",
    "                best_solution = xb[:]\n",
    "        plot(particles)\n",
    "        \n",
    "        if iteration % 100 == 0:\n",
    "            # Change target every 100 iterations\n",
    "            target[0] = random.randint(-10000, 10000)\n",
    "            target[1] = random.randint(-10000, 10000)\n",
    "            best_solution = min(particles, key=lambda x: get_fitness(x[0]))[0]\n",
    "            best_fitness = get_distance(best_solution)\n",
    "            #print(f\"New target: {target}\")\n",
    "            \n",
    "\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "demo_pso()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "im_ani = animation.ArtistAnimation(fig, ims, interval=50, repeat_delay=3000, blit=True)\n",
    "HTML(im_ani.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "We now consider how to use PSO to solve our onemax problem. To represent a bitvector problem, each particle $i$ is a triple:\n",
    "- Current position $x_i ∈ {0, 1}^n$\n",
    "- Best position so far $x_i* ∈ {0, 1}^n$\n",
    "- Velocity $v_i ∈ ℝ^n$\n",
    "\n",
    "The velocity $v_i$ is probabilistically translated into a new position for the particle (i.e., a new solution) using the sigmoid function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    return 1/(1 + math.exp(-x))\n",
    "\n",
    "xmin, xmax = -configuration[\"vmax\"], configuration[\"vmax\"]\n",
    "points = 1000\n",
    "xlist = list(map(lambda x: xmin + float(xmax - xmin)*x/points, range(points+1)))\n",
    "ylist = list(map(lambda y: sigmoid(y), xlist))\n",
    "plt.plot(xlist, ylist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "When updating a particle, we set $a = b = c = 1$, and $r1, r2 ∈ [0,2]$. \n",
    "\n",
    "- vi’ = vi  + r1 × (xi* - xi) + r2 × (xi** - xi)\n",
    "- v limited to [-vmax, vmax] with vmax often set to 4 or ln(n - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "configuration[\"c1\"] = 2\n",
    "configuration[\"c2\"] = 2\n",
    "configuration[\"vmax\"] = math.log(n - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def pso():\n",
    "    fitness_values.clear()\n",
    "    particles = []\n",
    "    for i in range(configuration[\"population_size\"]):\n",
    "        x  = [0 for _ in range(n)]\n",
    "        xb = [0 for _ in range(n)]\n",
    "        v  = [0 for _ in range(n)]\n",
    "        particles.append((x, xb, v))\n",
    "\n",
    "    best_solution = [0 for _ in range(n)]\n",
    "    best_fitness = get_fitness(best_solution)\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < max_steps:\n",
    "        iteration += 1\n",
    "\n",
    "        r1 = [random.uniform(0, configuration[\"c1\"]) for _ in range(n)]\n",
    "        r2 = [random.uniform(0, configuration[\"c2\"]) for _ in range(n)]\n",
    "\n",
    "        for (x, xb, v) in particles:\n",
    "            for i in range(n):\n",
    "                if random.random() < sigmoid(v[i]):\n",
    "                    x[i] = 1\n",
    "                else:\n",
    "                    x[i] = 0\n",
    "\n",
    "            fitness_current = get_fitness(x)\n",
    "            fitness_best    = get_fitness(xb)\n",
    "\n",
    "            if fitness_current > fitness_best:\n",
    "                for i in range(n):\n",
    "                    xb[i] = x[i]\n",
    "                fitness_best = fitness_current\n",
    "\n",
    "            if fitness_best > best_fitness:\n",
    "                best_fitness = fitness_best\n",
    "                best_solution = xb[:]\n",
    "\n",
    "        for (x, xb, v) in particles:\n",
    "            for i in range(n):\n",
    "                v[i] = v[i] + r1[i] * (xb[i] - x[i]) + r2[i] * (best_solution[i] - x[i])\n",
    "                v[i] = max(v[i], -configuration[\"vmax\"])\n",
    "                v[i] = min(v[i], configuration[\"vmax\"])\n",
    "\n",
    "        fitness_values.append(best_fitness)\n",
    "\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "pso()\n",
    "pso_values = fitness_values[:]\n",
    "\n",
    "plt.ylabel('Fitness', fontsize=15)\n",
    "plt.plot(random_values, label = \"Random\")\n",
    "plt.plot(ga_values, label = \"GA\")\n",
    "plt.plot(mmas_values, label = \"MMAS*\")\n",
    "plt.plot(immas_values, label = \"MMAS-Iteration\")\n",
    "plt.plot(pso_values, label = \"PSO\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "A special variant of PSO is 1-PSO, which works with just 1 particle. Own and global best are the same, so c2=0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "def onepso():\n",
    "    fitness_values.clear()\n",
    "    \n",
    "    v = [0 for _ in range(n)]\n",
    "    best_solution = [0 for _ in range(n)]\n",
    "    best_fitness = get_fitness(best_solution)\n",
    "\n",
    "    iteration = 0\n",
    "    while iteration < max_steps * configuration[\"population_size\"]:\n",
    "        iteration += 1\n",
    "\n",
    "        r = [random.uniform(0, configuration[\"c1\"]) for _ in range(n)]\n",
    "\n",
    "        x = []\n",
    "        for i in range(n):\n",
    "            p = 1/(1 + math.exp(-v[i]))\n",
    "            if random.random() < p:\n",
    "                x.append(1)\n",
    "            else:\n",
    "                x.append(0)\n",
    "\n",
    "        fitness_current = get_fitness(x)\n",
    "\n",
    "        if fitness_current > best_fitness:\n",
    "            best_fitness = fitness_current\n",
    "            best_solution = x[:]\n",
    "\n",
    "        for i in range(n):\n",
    "            v[i] = v[i] + r[i] * (best_solution[i] - x[i])\n",
    "            v[i] = max(v[i], -configuration[\"vmax\"])\n",
    "            v[i] = min(v[i], configuration[\"vmax\"])\n",
    "\n",
    "        if iteration % configuration[\"population_size\"] == 0:\n",
    "            fitness_values.append(best_fitness)\n",
    "\n",
    "\n",
    "    return best_solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "onepso()\n",
    "onepso_values = fitness_values[:]\n",
    "\n",
    "plt.ylabel('Fitness', fontsize=15)\n",
    "plt.plot(random_values, label = \"Random\")\n",
    "plt.plot(ga_values, label = \"GA\")\n",
    "plt.plot(mmas_values, label = \"MMAS*\")\n",
    "plt.plot(immas_values, label = \"MMAS-Iteration\")\n",
    "plt.plot(pso_values, label = \"PSO\")\n",
    "plt.plot(pso_values, label = \"1-PSO\")\n",
    "plt.legend()"
   ]
  }
 ],
 "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
}