lessons/6-Other/21-GeneticAlgorithms/Genetic.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "collapsed": true
      },
      "source": [
        "# Genetic Algorithms\n",
        "\n",
        "This notebook is part of [AI for Beginners Curriculum](http://github.com/microsoft/ai-for-beginners)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "import random\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import math\n",
        "import time"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Some Theory\n",
        "\n",
        "**Genetic Algorithms** (GA) are based on **evolutionary approach** to AI, in which methods of evolution of population is used to obtain an optimal solution for a given problem. They were proposed in 1975 by [John Henry Holland](https://en.wikipedia.org/wiki/John_Henry_Holland).\n",
        "\n",
        "Genetic Algorithms are based on the following ideas:\n",
        "* Valid solutions to the problem can be represented as **genes**\n",
        "* **Crossover** allows us to combine two solutions together to obtain new valid solution\n",
        "* **Selection** is used to select more optimal solutions using some **fitness function**\n",
        "* **Mutations** are introduced to destabilize optimization and get us out of the local minimum \n",
        "\n",
        "If you want to implement a Genetic Algorithm, you need the following:\n",
        "\n",
        " * To find a method of coding our problem solutions using **genes** $g\\in\\Gamma$\n",
        " * On the set of genes $\\Gamma$ we need to define **fitness function** $\\mathrm{fit}: \\Gamma\\to\\mathbb{R}$. Smaller function values would correspond to better solutions.\n",
        " * To define **crossover** mechanism to combine two genes together to get a new valid solution $\\mathrm{crossover}: \\Gamma^2\\to\\Gamma$.\n",
        " * To define **mutation** mechanism $\\mathrm{mutate}: \\Gamma\\to\\Gamma$.\n",
        "In many cases, crossover and mutation are quite simple algorithms to manipulate genes as numeric sequences or bit vectors.\n",
        "\n",
        "Specific implementation of a genetic algorithm can vary from case to case, but overall structure is the following:\n",
        "\n",
        "1. Select initial population $G\\subset\\Gamma$\n",
        "2. Randomly select one of the operations that will be performed at this step: crossover or mutation \n",
        "3. **Crossover**:\n",
        "  * Randomly select two genes $g_1, g_2 \\in G$\n",
        "  * Compute crossover $g=\\mathrm{crossover}(g_1,g_2)$\n",
        "  * If $\\mathrm{fit}(g)<\\mathrm{fit}(g_1)$ or $\\mathrm{fit}(g)<\\mathrm{fit}(g_2)$ - replace corresponding gene in the population by $g$.\n",
        "4. **Mutation** - select random gene $g\\in G$ and replace it by $\\mathrm{mutate}(g)$\n",
        "5. Repeat from step 2, until we get sufficiently small value of $\\mathrm{fit}$, or until the limit on the number of steps is reached.\n",
        "\n",
        "Tasks typically solved by GA:\n",
        "1. Schedule optimization\n",
        "1. Optimal packing\n",
        "1. Optimal cutting\n",
        "1. Speeding up exhaustive search\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problem 1: Fair Treasure Split\n",
        "\n",
        "**Task**: \n",
        "Two people found a treasure that contains diamonds of different sizes (and, correspondingly, different price). They need to split the treasure in two parts in such a way that the difference in the price is 0 (or minimal).\n",
        "\n",
        "**Formal definition**: \n",
        "We have a set of numbers $S$. We need to split it into two subsets $S_1$ and $S_2$, such that $$\\left|\\sum_{i\\in S_1}i - \\sum_{j\\in S_2}j\\right|\\to\\min$$ and $S_1\\cup S_2=S$, $S_1\\cap S_2=\\emptyset$.\n",
        "\n",
        "First of all, let's define the set $S$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[8344 2197 9335 3131 5863 9429 3818 9791   15 5455 1396 9538 4872 6549\n",
            " 8587 5986 6021 9764 8102 5083 5739 7684 8498 3007 6599  820 7490 2372\n",
            " 9370 5235 3525 3154  859 1906 8159 3950 2173 2988 2050  349 8713 2284\n",
            " 4177 6033 1651 9176 5049 8201  171 5081 1216 3756 4711 2757 7738 1272\n",
            " 5650 6584 5395 9004 7797  969 8104 1283 1392 4001 5768  445  274  256\n",
            " 8239 8015 4381 9021 1189 8879 1411 3539 6526 8011  136 7230 2332  451\n",
            " 5702 2989 4320 2446 9578 8486 4027 2410 9588 8981 2177 1493 3232 9151\n",
            " 4835 5594 6859 8394  369 3200  126 4259 2283 7755 2014 2458 8327 8082\n",
            " 7413 7622 1206 5533 8751 3495 5868 8472 6850 3958 3149 4672 4810 6274\n",
            " 4700 6134 4627 4616 6656 9949  884 2256 7419 1926 7973 5319 5967 9158\n",
            " 3823 7697 9466 5675 5412 9784 5426 8209 3421 1136 6047 4429 8001 4417\n",
            " 1381  722 7350 6018 6235 7860 5853 7660 5937 6242    1 9552 3971 8302\n",
            " 2633 9227 7283  154 8599 4269 9392 8539 1630  368 2409 9351 3838 9814\n",
            " 6186 5743 5083 1325 1610  779 3643 3262 5768 8725  961 4611 6310 4788\n",
            " 1648 5951 8118 7779]\n"
          ]
        }
      ],
      "source": [
        "N = 200\n",
        "S = np.array([random.randint(1,10000) for _ in range(N)])\n",
        "print(S)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's encode each possible solution of the problem by a binary vector $B\\in\\{0,1\\}^N$, where the number on $i$-th position shows to which of the sets ($S_1$ or $S_2$) the $i$-th number in the original set $S$ belongs. `generate` function will generate those random binary vectors."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[1 0 0 1 1 1 1 1 0 1 1 1 0 0 1 0 1 1 1 0 0 1 1 0 1 1 0 0 1 0 1 0 1 0 1 1 1\n",
            " 0 1 1 1 0 1 0 0 1 0 0 1 1 0 1 0 1 1 0 0 1 0 0 0 1 1 0 1 1 0 0 0 0 1 0 1 0\n",
            " 1 0 0 0 0 0 1 1 0 1 0 0 1 0 1 0 0 1 1 0 0 1 1 1 0 0 1 1 0 1 1 0 0 0 0 1 1\n",
            " 1 0 1 0 0 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 0 1 0 1 0 1 0 0 1 1 1 0 0 1 1\n",
            " 0 1 1 0 1 1 0 0 0 1 1 0 0 0 0 0 0 0 0 1 1 0 1 1 1 0 0 1 1 0 1 1 0 0 1 1 0\n",
            " 0 0 0 1 0 1 1 0 1 1 0 1 0 0 0]\n"
          ]
        }
      ],
      "source": [
        "def generate(S):\n",
        "    return np.array([random.randint(0,1) for _ in S])\n",
        "\n",
        "b = generate(S)\n",
        "print(b)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's now define `fit` function that calculates the \"cost\" of the solution. It will be the difference between sum or two sets, $S_1$ and $S_2$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "133784"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def fit(B,S=S):\n",
        "    c1 = (B*S).sum()\n",
        "    c2 = ((1-B)*S).sum()\n",
        "    return abs(c1-c2)\n",
        "\n",
        "fit(b)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we need to define functions for mutation and crossover:\n",
        "* For mutation, we will select one random bit and negate it (change from 0 to 1 and vice versa)\n",
        "* For crossover, we will take some bits from one vector, and some bits from another one. We will use the same `generate` function to randomly select, which bits to take from which of the input masks."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "def mutate(b):\n",
        "    x = b.copy()\n",
        "    i = random.randint(0,len(b)-1)\n",
        "    x[i] = 1-x[i]\n",
        "    return x\n",
        "\n",
        "def xover(b1,b2):\n",
        "    x = generate(b1)\n",
        "    return b1*x+b2*(1-x)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's create initial population of the solutions $P$ of the size `pop_size`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "pop_size = 30\n",
        "P = [generate(S) for _ in range(pop_size)]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now, the main function to perform the evolution. `n` is the number of steps of evolution to undergo. At each step:\n",
        "* With the probability of 30% we perform a mutation, and replace the element with the worst `fit` function by the mutated element\n",
        "* With the probability of 70% we perform crossover\n",
        "\n",
        "The function returns the best solution (gene corresponding to the best solution), and the history of minimal fit function in the population on each iteration."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[0 0 0 1 1 0 0 0 0 1 1 1 0 1 0 0 0 1 0 1 0 1 0 1 0 1 1 0 0 0 0 0 1 0 1 1 0\n",
            " 0 0 0 1 1 0 0 1 0 0 0 0 0 1 0 0 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 0 1 0 0 0\n",
            " 0 1 1 1 0 1 0 1 1 1 1 1 0 0 0 1 1 0 1 0 0 1 0 0 1 1 1 1 1 1 1 1 0 1 0 1 1\n",
            " 0 1 1 0 0 0 0 1 1 1 1 0 1 0 0 1 0 1 1 1 0 1 0 0 0 0 0 0 1 1 0 0 0 1 1 0 0\n",
            " 1 0 1 1 1 1 1 0 1 0 1 0 1 1 1 0 0 0 1 1 0 0 0 0 0 0 1 1 1 0 1 0 0 0 1 0 1\n",
            " 0 1 0 1 0 0 1 1 1 0 1 1 0 0 1] 4\n"
          ]
        }
      ],
      "source": [
        "def evolve(P,S=S,n=2000):\n",
        "    res = []\n",
        "    for _ in range(n):\n",
        "        f = min([fit(b) for b in P])\n",
        "        res.append(f)\n",
        "        if f==0:\n",
        "            break\n",
        "        if random.randint(1,10)<3:\n",
        "            i = random.randint(0,len(P)-1)\n",
        "            b = mutate(P[i])\n",
        "            i = np.argmax([fit(z) for z in P])\n",
        "            P[i] = b\n",
        "        else:\n",
        "            i = random.randint(0,len(P)-1)\n",
        "            j = random.randint(0,len(P)-1)\n",
        "            b = xover(P[i],P[j])\n",
        "            if fit(b)<fit(P[i]):\n",
        "                P[i]=b\n",
        "            elif fit(b)<fit(P[j]):\n",
        "                P[j]=b\n",
        "            else:\n",
        "                pass\n",
        "    i = np.argmin([fit(b) for b in P])\n",
        "    return (P[i],res)\n",
        "\n",
        "(s,hist) = evolve(P)\n",
        "print(s,fit(s))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "You can see that we have managed to minimize the `fit` function quite a bit! Here is the graph that shows how the `fit` function for the whole population behaves during the process."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYD0lEQVR4nO3df5Bd5V3H8fdnd0MafgSSZhPDJjQBU8cENTRrTKUqSjUpo4aq7aTjNBllJspQbf0xI1i1dZyMrVo6ooKmwpA4tDROW4kVtDFDxdqUsGAgCSGwNEC2WbMrlCEUJySbr3/c5+69bO/+uLv3nns55/OauXPPPvece797dvPJs895zjmKCMzMrBg6Wl2AmZllx6FvZlYgDn0zswJx6JuZFYhD38ysQLpaXcBkFixYEMuWLWt1GWZmbyqPPvro/0ZE99j2tg/9ZcuW0dfX1+oyzMzeVCQ9X6vdwztmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFUjbz9Ofrrv/6xgvfed1OjrE+3uXcuklc1pdkplZy+U29D+7/wWePvkqAJ0Sv3HtihZXZGbWerkd3vnKb/0Ex/70OgBGfKMYMzMgx6FfzZlvZlZSiNA3M7OSXIe+JADc0TczK8l16JuZ2Rs59M3MCqQYoe8juWZmQAFCPw3rm5kZUwh9SW+RtF/S45IOS/rj1D5f0h5Jz6TneVXb3CKpX9JRSeur2tdIOpheu03KJpLdzzczK5lKT/808FMR8UPAamCDpHXAzcDeiFgB7E1fI2klsAlYBWwAbpfUmd7rDmArsCI9NjTuW6nNHX0zs4pJQz9KXk1fzkqPADYCO1L7DuD6tLwRuDciTkfEMaAfWCtpMTA3IvZFRAA7q7ZpKg/pm5mVTGlMX1KnpAPAELAnIh4GFkXEIEB6XphW7wGOV20+kNp60vLY9lqft1VSn6S+4eHhOr4dMzObyJRCPyJGImI1sIRSr/3KCVavNaISE7TX+rztEdEbEb3d3d1TKXH8Ynwk18xsVF2zdyLiZeCrlMbiT6YhG9LzUFptAFhatdkS4ERqX1KjvenCh3LNzICpzd7plnRJWp4DvBt4CtgNbEmrbQHuS8u7gU2SZktaTumA7f40BHRK0ro0a2dz1TZN436+mVnFVK6nvxjYkWbgdAC7IuLLkvYBuyTdALwAvA8gIg5L2gU8CZwFboqIkfReNwJ3A3OAB9Kj6Xwg18ysZNLQj4gngKtqtL8IXDvONtuAbTXa+4CJjgc0nIf0zcwqcn9GLvjkLDOzstyHvjyqb2Y2KvehDx7TNzMrK0Tom5lZSf5D36M7Zmaj8h/6+OQsM7Oy3Ie+O/pmZhW5D33AczbNzJLch75PzjIzq8h96IM7+mZmZbkPfZ+cZWZWkfvQNzOzikKEfviUXDMzoACh7wO5ZmYVuQ998LV3zMzKch/67uibmVXkPvTBUzbNzMpyH/ryoL6Z2ajch76ZmVXkPvSFD+SamZXlPvTNzKyiEKHv6+mbmZVMGvqSlkp6UNIRSYclfTi1f1zStyQdSI/rqra5RVK/pKOS1le1r5F0ML12m7I4yurjuGZmo7qmsM5Z4Hci4jFJFwGPStqTXvt0RPxF9cqSVgKbgFXApcC/S3p7RIwAdwBbgW8A9wMbgAca862Mz2P6ZmYlk/b0I2IwIh5Ly6eAI0DPBJtsBO6NiNMRcQzoB9ZKWgzMjYh9UboYzk7g+pl+A5NxR9/MrKKuMX1Jy4CrgIdT04ckPSHpLknzUlsPcLxqs4HU1pOWx7bX+pytkvok9Q0PD9dTopmZTWDKoS/pQuALwEci4hVKQzVXAKuBQeBT5VVrbB4TtH93Y8T2iOiNiN7u7u6pljhe3TPa3swsT6YU+pJmUQr8eyLiiwARcTIiRiLiHPAZYG1afQBYWrX5EuBEal9So93MzDIyldk7Au4EjkTErVXti6tWey9wKC3vBjZJmi1pObAC2B8Rg8ApSevSe24G7mvQ9zFB/b6evplZ2VRm71wNfBA4KOlAavt94AOSVlMaonkO+DWAiDgsaRfwJKWZPzelmTsANwJ3A3Mozdpp+swdMzOrmDT0I+Jr1B6Pv3+CbbYB22q09wFX1lNgI7ifb2ZWkvszcn0Y18ysIvehDz45y8ysLPeh7ymbZmYVuQ99MzOryH3oC19l08ysLPehb2ZmFbkP/dLJWa2uwsysPeQ+9M3MrKIQoe+OvplZSQFC31M2zczKChD6ZmZWlvvQ94FcM7OK3Ie+mZlV5D70SyP67uqbmUEBQt/MzCpyH/oe0zczq8h96JuZWYVD38ysQHIf+kIe3jEzS3If+mZmVpH70Jd8PX0zs7Lch76ZmVVMGvqSlkp6UNIRSYclfTi1z5e0R9Iz6Xle1Ta3SOqXdFTS+qr2NZIOptduUwY3sBWesmlmVjaVnv5Z4Hci4vuBdcBNklYCNwN7I2IFsDd9TXptE7AK2ADcLqkzvdcdwFZgRXpsaOD3YmZmk5g09CNiMCIeS8ungCNAD7AR2JFW2wFcn5Y3AvdGxOmIOAb0A2slLQbmRsS+iAhgZ9U2TSPJI/pmZkldY/qSlgFXAQ8DiyJiEEr/MQAL02o9wPGqzQZSW09aHtte63O2SuqT1Dc8PFxPiWZmNoEph76kC4EvAB+JiFcmWrVGW0zQ/t2NEdsjojcieru7u6daopmZTWJKoS9pFqXAvycivpiaT6YhG9LzUGofAJZWbb4EOJHal9RobzofyDUzK5nK7B0BdwJHIuLWqpd2A1vS8hbgvqr2TZJmS1pO6YDt/jQEdErSuvSem6u2MTOzDHRNYZ2rgQ8CByUdSG2/D3wC2CXpBuAF4H0AEXFY0i7gSUozf26KiJG03Y3A3cAc4IH0aCqfnGVmVjFp6EfE1xj/7uLXjrPNNmBbjfY+4Mp6CjQzs8bJ/Rm5Er5xlplZkvvQNzOzityHvsYdmTIzK57chz54dMfMrKwQoW9mZiW5D/3SjdHd1zczgwKEvpmZVUzl5Kw3NQEvvXaGA8dfnnTdRXNns/jiOU2vycysVXIf+nPO6+Khp4d56OnJr9Z50ewuDv7x+knXMzN7s8p96G//4Br6h16ddL1/OvAt7jtwgogggxt6mZm1RO5Df+n881k6//xJ13t84OXmF2Nm1mI+kJuUT+LyRB8zyzOHflIe0XHmm1meOfTH8Jx+M8szh35SPnTryDezPHPoJ56wY2ZF4NBPytM0PbpjZnnm0B/Dt1Y0szxz6I/hnr6Z5ZlDP/GYvpkVgUM/8R22zKwIHPpjeHjHzPJs0tCXdJekIUmHqto+Lulbkg6kx3VVr90iqV/SUUnrq9rXSDqYXrtNbXZVs8oZuU59M8uvqfT07wY21Gj/dESsTo/7ASStBDYBq9I2t0vqTOvfAWwFVqRHrfdsmdGTs5z5ZpZjk4Z+RDwEvDTF99sI3BsRpyPiGNAPrJW0GJgbEfuidJ2DncD106y5Kdrr7w4zs+aYyZj+hyQ9kYZ/5qW2HuB41ToDqa0nLY9tr0nSVkl9kvqGhye/+UkjuaNvZnk23dC/A7gCWA0MAp9K7bX6yzFBe00RsT0ieiOit7u7e5ol1qdyaWXHvpnl17RCPyJORsRIRJwDPgOsTS8NAEurVl0CnEjtS2q0tw1fWtnMimBaoZ/G6MveC5Rn9uwGNkmaLWk5pQO2+yNiEDglaV2atbMZuG8GdZuZ2TRMertESZ8DrgEWSBoAPgZcI2k1pY7xc8CvAUTEYUm7gCeBs8BNETGS3upGSjOB5gAPpEfb8eiOmeXZpKEfER+o0XznBOtvA7bVaO8DrqyrugzJ4ztmVgA+Izep3ETFqW9m+eXQT0Y7+s58M8sxh37ic7PMrAgc+mO4o29meebQTyq3S3Tsm1l+OfQTT94xsyJw6Cce0zezInDoj+HRHTPLM4d+WXlM3wM8ZpZjDv1kdHjHmW9mOebQT3wg18yKwKFvZlYgDv2kchOVFhdiZtZEDv2kMrzj1Dez/HLoJ6NX2XTmm1mOOfTHcOabWZ459BP5lFwzKwCHflI5kOu+vpnll0O/zDdRMbMCcOibmRWIQz/xkL6ZFYFDP6ncRKXFhZiZNdGkoS/pLklDkg5Vtc2XtEfSM+l5XtVrt0jql3RU0vqq9jWSDqbXbpPaa77M6Dx9T9o0sxybSk//bmDDmLabgb0RsQLYm75G0kpgE7AqbXO7pM60zR3AVmBFeox9z7bgnr6Z5dmkoR8RDwEvjWneCOxIyzuA66va742I0xFxDOgH1kpaDMyNiH1RmhO5s2qbtuCrbJpZEUx3TH9RRAwCpOeFqb0HOF613kBq60nLY9trkrRVUp+kvuHh4WmWWJ/2GmwyM2uORh/IrRWdMUF7TRGxPSJ6I6K3u7u7YcVNxCdnmVkRdE1zu5OSFkfEYBq6GUrtA8DSqvWWACdS+5Ia7W2j3NPfue955p1/3pS3W3f5fH7k8rc2qSozs8aabujvBrYAn0jP91W1f1bSrcCllA7Y7o+IEUmnJK0DHgY2A381o8obbMm88zmvs4O7v/5cXdv9QM/F/PNvvKs5RZmZNdikoS/pc8A1wAJJA8DHKIX9Lkk3AC8A7wOIiMOSdgFPAmeBmyJiJL3VjZRmAs0BHkiPtrHmbfN46k/qm1D0W7sO8N8vvNycgszMmmDS0I+ID4zz0rXjrL8N2FajvQ+4sq7qMtbRUd/R3IvnzOL4t1/jnX+6t+brF8zuYuevruXSS+Y0ojwzsxmb7vCOAe/vXcrpM+dqntA1fOo0Dx4d5rkXv+PQN7O24dCfgSt7LuaTv/SDNV/bf+wlHjw67In/ZtZWfO2dJinPBjrn0DezNuLQb5IO32jdzNqQQ79pSqnvnr6ZtROHfpN0jA7vOPXNrH049JtEvoKbmbUhh36TuKdvZu3Iod8klQu4tbgQM7MqDv0mkXv6ZtaGHPpN4iF9M2tHDv0m6ZCvz29m7ceh3ySjPX1nvpm1EYd+k5R7+j45y8zaiUO/ScoXafaBXDNrJw79JimfnOXIN7N24tBvksqYvmPfzNqHQ79JKrN3WlyImVkVh36TeEzfzNqRQ79J3NM3s3bk0G8SX4bBzNqRQ79JfHKWmbWjGYW+pOckHZR0QFJfapsvaY+kZ9LzvKr1b5HUL+mopPUzLb6djQ7veNKmmbWRRvT0fzIiVkdEb/r6ZmBvRKwA9qavkbQS2ASsAjYAt0vqbMDntyXfGN3M2lFXE95zI3BNWt4BfBX4vdR+b0ScBo5J6gfWAvuaUEPLlXv6t+55mju/dmxa73HjT1zBL65Z0siyzKzgZhr6AXxFUgB/FxHbgUURMQgQEYOSFqZ1e4BvVG07kNq+i6StwFaAyy67bIYltkb3hbP5lauXMfTK6Wlt/+DRIf7zmWGHvpk11ExD/+qIOJGCfY+kpyZYVzXaag5+pP88tgP09va+KQdIOjrEx35u1bS3v/ZTX+XMyJvyWzezNjajMf2IOJGeh4AvURquOSlpMUB6HkqrDwBLqzZfApyYyefn2azODs6MnGt1GWaWM9MOfUkXSLqovAz8DHAI2A1sSattAe5Ly7uBTZJmS1oOrAD2T/fz886hb2bNMJPhnUXAl9LVJLuAz0bEv0p6BNgl6QbgBeB9ABFxWNIu4EngLHBTRIzMqPoc6+oUZz31x8wabNqhHxHfBH6oRvuLwLXjbLMN2DbdzyySWZ0dvH7WPX0zayyfkdumZnWKh4+9xLe/83qrSzGzHHHot6llb70AgKdPnmpxJWaWJw79NvWzP3gp4DN6zayxHPptqsN33jKzJnDot6mOlPojDn0zayCHfpsqX7vHwztm1kgO/TbV4ZuwmFkTOPTbVOV2iw59M2sch36bKoe+r8RgZo3k0G9THekn4+EdM2skh36b8vCOmTWDQ79NefaOmTWDQ79NefaOmTWDQ79NjZ6c5a6+mTWQQ79NVcb0W1yImeWKQ79NeXjHzJrBod+mfCDXzJphJrdLtCYqj+l/+YkT9A+92uJqJnfVZZewftX3tLoMM5uEQ79NzTt/Fj2XzOHrz77I1599sdXlTOjsyDkWzX2LQ9/sTcCh36bOP6+L/7r5p1pdxpRs+5cn2bnvefYeOVnz9dldnay7fD5dnR5NNGs1h77N2KWXzOH02XPcsKNv3HU+s7mXn165KMOqzKwWh77N2OZ3LuOHl82vOdPo5dfOsPmu/Xz16NCkM5GWL7iAty+6qFllmhktCH1JG4C/BDqBv4+IT2RdgzVWZ4e4sufimq+dGTnHBed1cs/DL3DPwy9M+D4LLpxN3x+8uxklmlmSaehL6gT+BvhpYAB4RNLuiHgyyzosO7M6O3jwd69h+NXTE673+UeOs3Pf8/zhPx0izVbNnR+94q1suHJxq8uwgsu6p78W6I+IbwJIuhfYCDj0c2zh3LewcO5bJlxnw6ozPHDof/jyEycyqipb337tDDv3Pc+cWZ11bdfVIRbOnT163oYVy5d/813M7qrvd2YyWYd+D3C86usB4EfGriRpK7AV4LLLLsumMmupH/3eBTzy0fwO7fQPneIf+wao91y7V0+f5eXXXm9KTdb+ROP/s8869Gt9B9/17yAitgPbAXp7e31Oqr3pfe/Ci7jluu9vdRlmmV+GYQBYWvX1EiCff8+bmbWhrEP/EWCFpOWSzgM2AbszrsHMrLAyHd6JiLOSPgT8G6Upm3dFxOEsazAzK7LM5+lHxP3A/Vl/rpmZ+dLKZmaF4tA3MysQh76ZWYE49M3MCkTR5vdglTQMPD/NzRcA/9vAchrFddXHddXHddUnr3W9LSK6xza2fejPhKS+iOhtdR1jua76uK76uK76FK0uD++YmRWIQ9/MrEDyHvrbW13AOFxXfVxXfVxXfQpVV67H9M3M7I3y3tM3M7MqDn0zswLJZehL2iDpqKR+STdn/NlLJT0o6Yikw5I+nNo/Lulbkg6kx3VV29ySaj0qaX0Ta3tO0sH0+X2pbb6kPZKeSc/zsqxL0vdV7ZMDkl6R9JFW7S9Jd0kaknSoqq3ufSRpTdrX/ZJuk2Z2v8Nx6vpzSU9JekLSlyRdktqXSfq/qn33txnXVffPLqO6Pl9V03OSDqT2TPbXBNmQ7e9XROTqQemSzc8ClwPnAY8DKzP8/MXAO9LyRcDTwErg48Dv1lh/ZapxNrA81d7ZpNqeAxaMafsz4Oa0fDPwyazrGvOz+x/gba3aX8CPA+8ADs1kHwH7gXdSulvcA8B7mlDXzwBdafmTVXUtq15vzPtkUVfdP7ss6hrz+qeAP8pyfzF+NmT6+5XHnv7ozdcj4nWgfPP1TETEYEQ8lpZPAUco3Rt4PBuBeyPidEQcA/opfQ9Z2QjsSMs7gOtbWNe1wLMRMdEZ2E2tKyIeAl6q8ZlT3keSFgNzI2JflP6F7qzapmF1RcRXIuJs+vIblO5EN66s6ppAS/dXWeoVvx/43ETv0ei6JsiGTH+/8hj6tW6+PlHoNo2kZcBVwMOp6UPpT/G7qv6Ey7LeAL4i6VGVbj4PsCgiBqH0SwksbEFdZZt44z/EVu+vsnr3UU9azrLGX6XU4ytbLum/Jf2HpB9LbVnWVc/PLuv99WPAyYh4pqot0/01Jhsy/f3KY+hP6ebrTS9CuhD4AvCRiHgFuAO4AlgNDFL68xKyrffqiHgH8B7gJkk/PsG6me5HlW6f+fPAP6amdthfkxmvlqz33UeBs8A9qWkQuCwirgJ+G/ispLkZ1lXvzy7rn+kHeGPnItP9VSMbxl11nM+fUV15DP2W33xd0ixKP9R7IuKLABFxMiJGIuIc8BkqQxKZ1RsRJ9LzEPClVMPJ9Odi+c/ZoazrSt4DPBYRJ1ONLd9fVerdRwO8cailaTVK2gL8LPDL6U990nDAi2n5UUpjwW/Pqq5p/Oyy3F9dwC8An6+qN7P9VSsbyPj3K4+h39Kbr6fxwjuBIxFxa1X74qrV3guUZxXsBjZJmi1pObCC0kGaRtd1gaSLysuUDgIeSp+/Ja22Bbgvy7qqvKH31er9NUZd+yj9iX5K0rr0+7C5apuGkbQB+D3g5yPitar2bkmdafnyVNc3M6yrrp9dVnUl7waeiojR4ZGs9td42UDWv1/TPRLdzg/gOkpHxp8FPprxZ7+L0p9aTwAH0uM64B+Ag6l9N7C4apuPplqPMsNZCxPUdTmlmQCPA4fL+wV4K7AXeCY9z8+yrvQ55wMvAhdXtbVkf1H6j2cQOEOpR3XDdPYR0Esp7J4F/pp09nuD6+qnNOZb/j3727TuL6af8ePAY8DPZVxX3T+7LOpK7XcDvz5m3Uz2F+NnQ6a/X74Mg5lZgeRxeMfMzMbh0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFcj/Awy4I/IZjuXDAAAAAElFTkSuQmCC",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          },
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.plot(hist)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Problem 2: N Queens Problem\n",
        "\n",
        "**Task**:\n",
        "You need to place $N$ queens on a chess board of the size $N\\times N$ in such a way that they do not attack each other.\n",
        "\n",
        "First of all, let's solve the problem without genetic algorithms, using full search. We can represent the state of the board by the list $L$, where $i$-th number in the list is the horizontal position of the queen in $i$-th row. It is quite obvious that each solution will have only one queen per row, and each row would have a queen.\n",
        "\n",
        "Our goal would be to find the first solution to the problem, after which we will stop searching. You can easily extend this function to generate all possible positions for queens."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[1, 5, 8, 6, 3, 7, 2, 4]\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "True"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "N = 8\n",
        "\n",
        "def checkbeats(i_new,j_new,l):\n",
        "    for i,j in enumerate(l,start=1):\n",
        "        if j==j_new:\n",
        "            return False\n",
        "        else:\n",
        "            if abs(j-j_new) == i_new-i:\n",
        "                return False\n",
        "    return True\n",
        "\n",
        "def nqueens(l,N=8,disp=True):\n",
        "    if len(l)==N:\n",
        "        if disp: print(l)\n",
        "        return True\n",
        "    else:\n",
        "        for j in range(1,N+1):\n",
        "            if checkbeats(len(l)+1,j,l):\n",
        "                l.append(j)\n",
        "                if nqueens(l,N,disp): return True\n",
        "                else: l.pop()\n",
        "        return False\n",
        "            \n",
        "nqueens([],8)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's measure how long does it take to get a solution for 20-queens problem:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "10.6 s ± 2.17 s per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
          ]
        }
      ],
      "source": [
        "%timeit nqueens([],20,False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's solve the same problem using genetic algorithm. This solution is inspired by [this blog post](https://kushalvyas.github.io/gen_8Q.html).\n",
        "\n",
        "We will represent each solution by the same list of length $N$, and as a `fit` function we will take the number of queens that attack each other:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "def fit(L):\n",
        "    x=0\n",
        "    for i1,j1 in enumerate(L,1):\n",
        "        for i2,j2 in enumerate(L,1):\n",
        "            if i2>i1:\n",
        "                if j2==j1 or (abs(j2-j1)==i2-i1): x+=1\n",
        "    return x"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Since calculating fitness function is time consuming, let's store each solution in the population together with the value of fitness function. Let's generate the initial population:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[(array([2, 3, 8, 7, 5, 4, 1, 6]), 4),\n",
              " (array([3, 4, 5, 1, 2, 8, 6, 7]), 8),\n",
              " (array([1, 3, 7, 4, 5, 8, 6, 2]), 6),\n",
              " (array([1, 5, 4, 6, 8, 3, 7, 2]), 4),\n",
              " (array([3, 5, 7, 1, 8, 6, 4, 2]), 3)]"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def generate_one(N):\n",
        "    x = np.arange(1,N+1)\n",
        "    np.random.shuffle(x)\n",
        "    return (x,fit(x))\n",
        "\n",
        "def generate(N,NP):\n",
        "    return [generate_one(N) for _ in range(NP)]\n",
        "\n",
        "generate(8,5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now we need to define mutation and crossover functions. Crossover would combine two genes together by breaking them at some random point and concatenating two parts from different genes together."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([1, 2, 7, 8])"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "def mutate(G):\n",
        "    x=random.randint(0,len(G)-1)\n",
        "    G[x]=random.randint(1,len(G))\n",
        "    return G\n",
        "    \n",
        "def xover(G1,G2):\n",
        "    x=random.randint(0,len(G1))\n",
        "    return np.concatenate((G1[:x],G2[x:]))\n",
        "\n",
        "xover([1,2,3,4],[5,6,7,8])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We will enhance gene selection process by selecting more genes with better fitness function. The probability of selection of a gene would depend on the fitness function: "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "def choose_rand(P):\n",
        "    N=len(P[0][0])\n",
        "    mf = N*(N-1)//2 # max fitness fn\n",
        "    z = [mf-x[1] for x in P]\n",
        "    tf = sum(z) # total fitness\n",
        "    w = [x/tf for x in z]\n",
        "    p = np.random.choice(len(P),2,False,p=w)\n",
        "    return p[0],p[1]\n",
        "\n",
        "def choose(P):\n",
        "    def ch(w):\n",
        "        p=[]\n",
        "        while p==[]:\n",
        "            r = random.random()\n",
        "            p = [i for i,x in enumerate(P) if x[1]>=r]\n",
        "        return random.choice(p)\n",
        "    N=len(P[0][0])\n",
        "    mf = N*(N-1)//2 # max fitness fn\n",
        "    z = [mf-x[1] for x in P]\n",
        "    tf = sum(z) # total fitness\n",
        "    w = [x/tf for x in z]\n",
        "    p1=p2=0\n",
        "    while p1==p2:\n",
        "        p1 = ch(w)\n",
        "        p2 = ch(w)\n",
        "    return p1,p2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let's define the main evolutionary loop. We will make the logic slightly different from previous example, to show that one can get creative. We will loop until we get the perfect solution (fitness function=0), and at each step we will take the current generation, and produce the new generation of the same size. This is done using `nxgeneration` function, using the following steps:\n",
        "\n",
        "1. Discard the most unfit solutions - there is `discard_unfit` function that does that\n",
        "1. Add some more random solutions to the generation\n",
        "1. Populate new generation of size `gen_size` using the following steps for each new gene:\n",
        "    - select two random genes, with probability proportional to fitness function\n",
        "    - calculate a crossover\n",
        "    - apply a mutation with the probability `mutation_prob`"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(array([4, 7, 5, 3, 1, 6, 8, 2]), 0)"
            ]
          },
          "execution_count": 20,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "mutation_prob = 0.1\n",
        "\n",
        "def discard_unfit(P):\n",
        "    P.sort(key=lambda x:x[1])\n",
        "    return P[:len(P)//3]\n",
        "\n",
        "def nxgeneration(P):\n",
        "    gen_size=len(P)\n",
        "    P = discard_unfit(P)\n",
        "    P.extend(generate(len(P[0][0]),3))\n",
        "    new_gen = []\n",
        "    for _ in range(gen_size):\n",
        "        p1,p2 = choose_rand(P)\n",
        "        n = xover(P[p1][0],P[p2][0])\n",
        "        if random.random()<mutation_prob:\n",
        "            n=mutate(n)\n",
        "        nf = fit(n)\n",
        "        new_gen.append((n,nf))\n",
        "        '''\n",
        "        if (nf<=P[p1][1]) or (nf<=P[p2][1]):\n",
        "            new_gen.append((n,nf))\n",
        "        elif (P[p1][1]<P[p2][1]):\n",
        "            new_gen.append(P[p1])\n",
        "        else:\n",
        "            new_gen.append(P[p2])\n",
        "        '''\n",
        "    return new_gen\n",
        "    \n",
        "def genetic(N,pop_size=100):\n",
        "    P = generate(N,pop_size)\n",
        "    mf = min([x[1] for x in P])\n",
        "    n=0\n",
        "    while mf>0:\n",
        "        #print(\"Generation {0}, fit={1}\".format(n,mf))\n",
        "        n+=1\n",
        "        mf = min([x[1] for x in P])\n",
        "        P = nxgeneration(P)\n",
        "    mi = np.argmin([x[1] for x in P])\n",
        "    return P[mi]\n",
        "\n",
        "genetic(8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It is interesting that in most of the times we are able to get a solution pretty quickly, but in some rare cases optimization reaches local minimum, and the process is stuck for a long time. It is important to take that into account when you are measuring average time: while in most of the cases genetic algorithm will be faster than full search, in some cases it can take longer. To overcome this problem, it often makes sense to limit the number of generations to consider, and if we are not able to find the solution - we can start from scratch.  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "trusted": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The slowest run took 18.71 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
            "26.4 s ± 28.7 s per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
          ]
        }
      ],
      "source": [
        "%timeit genetic(10)"
      ]
    }
  ],
  "metadata": {
    "interpreter": {
      "hash": "86193a1ab0ba47eac1c69c1756090baa3b420b3eea7d4aafab8b85f8b312f0c5"
    },
    "kernelspec": {
      "display_name": "Python 3.6",
      "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.9.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}