tutorials/03_Naive_approach.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# TSP + QAOA - naive approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to run this line to have a proper [bra-ket](https://en.wikipedia.org/wiki/Bra–ket_notation) formatting in markdown.\n",
    "$$\n",
    "\\newcommand{\\ket}[1]{\\left| #1 \\right\\rangle}\n",
    "\\newcommand{\\bra}[1]{\\left\\langle #1 \\right|}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the previous tutorials we learned what is the Traveling Salesman Problem and how to use QAOA to solve optimization problems.\n",
    "<br>\n",
    "Now it's finally time to solve TSP using quantum computing!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial you will learn the following:\n",
    "\n",
    "- How to encode TSP for QAOA\n",
    "- How to solve constrained optimization problems using QAOA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I realize that it might be hard to grasp all the concepts here if you are new to the field. It sometimes gets pretty low-level and mathematical. But that's not because it's inherently difficult, rather:\n",
    "- we (or I at least) have not yet discovered how to explain it in the most efficient way\n",
    "- the QAOA was invented in 2014 - it's pretty new stuff. \n",
    "\n",
    "So don't get discouraged, asked others for help if you something seems too complicated or just get some rest and come back with a fresh mind. And remember - you're pioneer!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Encoding the problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the important decisions we need to make with any problem we encounter is how to encode our input and output data. What will be the good choice depends on many factors:\n",
    "- what algorithm we use?\n",
    "- do we have any constraints (e.g. amount of data, memory)?\n",
    "- what hardware do we use?\n",
    "- many others"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the most convenient way of encoding solution to TSP is the one we used in the first part:\n",
    "\n",
    "0 -> 1 -> 2 -> 3 , which for Python translates into `[0,1,2,3]` or `(0,1,2,3)`.\n",
    "\n",
    "This means that we start from the city 0, then go to 1 and so on.\n",
    "\n",
    "However, in case of quantum computers, one of the main constraints is that we need to use binary representation. Qubits can take only 0 or 1 values during the measurement, so this format won't work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know that at any given time we can be only in one city. Hence, for N cities we have N time slots.\n",
    "<br>\n",
    "Using binary representation, we could use N strings (representing time slots) of N bits (representing cities), where the active bit in each string is the city that we visit.\n",
    "\n",
    "So in the case of [0,1,2,3] we will get the following matrix (row represents time slots, columns - cities):\n",
    "\n",
    "$\\begin{bmatrix}\n",
    "1 & 0 & 0 & 0 \\\\\n",
    "0 & 1 & 0 & 0 \\\\\n",
    "0 & 0 & 1 & 0 \\\\\n",
    "0 & 0 & 0 & 1 \\\\\n",
    "\\end{bmatrix}$\n",
    "\n",
    "Which translates into the following string of bits:\n",
    "\n",
    "`[1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1]`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a little bit unwieldy, so we will use the following helper methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "def points_order_to_binary_state(points_order):\n",
    "    \"\"\"\n",
    "    Transforms the order of points from the standard representation: [0, 1, 2],\n",
    "    to the binary one: [1,0,0,0,1,0,0,0,1]\n",
    "    \"\"\"\n",
    "    number_of_points = len(points_order)\n",
    "    binary_state = np.zeros((len(points_order))**2)\n",
    "    for j in range(len(points_order)):\n",
    "        p = points_order[j]\n",
    "        binary_state[(number_of_points) * (j) + (p)] = 1\n",
    "    return binary_state\n",
    "\n",
    "def binary_state_to_points_order(binary_state):\n",
    "    \"\"\"\n",
    "    Transforms the the order of points from the binary representation: [1,0,0,0,1,0,0,0,1],\n",
    "    to the binary one: [0, 1, 2]\n",
    "    \"\"\"\n",
    "    points_order = []\n",
    "    number_of_points = int(np.sqrt(len(binary_state)))\n",
    "    for p in range(number_of_points):\n",
    "        for j in range(number_of_points):\n",
    "            if binary_state[(number_of_points) * p + j] == 1:\n",
    "                points_order.append(j)\n",
    "    return points_order"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see if it works on an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Binary representation of  [0, 1, 2] is [1. 0. 0. 0. 1. 0. 0. 0. 1.]\n",
      "[0, 1, 2] should be equal to [0, 1, 2]\n"
     ]
    }
   ],
   "source": [
    "points_order = [0, 1, 2]\n",
    "binary_state = points_order_to_binary_state(points_order)\n",
    "points_order_again = binary_state_to_points_order(binary_state)\n",
    "\n",
    "print(\"Binary representation of \", points_order, \"is\", binary_state)\n",
    "print(points_order_again, \"should be equal to\", points_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving TSP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know how to encode our problem, we will use `ForestTSPSolverNaive` class to solve it. We will start with 3 cities since the calcultions are faster this way. We will go to more cities when we are sure that our code is working properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scripts.forest_tsp_solver_naive import ForestTSPSolverNaive"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFnJJREFUeJzt3X+M3HWdx/Hnq7O7pbToenQJTbtlkR+XqFDATa+EKBs4klJJayImJR5SUKs9UIwmdx53wZOcdzExelEuJQ00Vo8IpHC6kPa0d1AUkxa2XEspCCwEbUmxa6HFltLt7L7vj/mqy3S2892Z7+7M9Pt6JJN+5/v9fL/z/uwXXvvd73y+368iAjMzO/lNa3QBZmY2NRz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58AKkbaT3SQaS3kB5Cmt/osszMsqTcj8OXTgV2AEeBfwIC+BfgVOBCIg43sDozs8y0NbqAJvA54P3AXxIxCID0DPAS8HngO40rzcwsOz7Cl/4XOIWIy8rmPw5AxOUNqMrMLHM+hw8fBJ6tMH8X8IEprsXMbNI48OEvgDcrzH8DeN8U12JmNmkc+GZmOdGwc/izZ8+Onp6ehnz2WD/bsYPHOzv517POetf8v//tb/nrN9/kqgULGlSZmdnxtm3b9vuI6Kpl3YYFfm9vbwwMDDTks9/liitgeBieeOLd8/v6IAIef7whZZmZVSJpW0T01rKuT+ksXQpbtsArr/x53quvwq9+VVpmZnaScOB/7nPQ0wPLlsFPfwr9/aXp7m74/OcbXZ2ZWWZaOvCLI6O89c4xRkbrOC01cyY8+iicfz5cfz186lNw9tmlebNmZVesmVmDtdyVtkeLI2zYuZfVm1/mpX2HaJsmiqPB+WfM4gt957DkgjlMbytMbKPz58ODD05OwWZmTaKlvrTdvvsAK9Y+ybGRUQ4Pjxy3fGZHgfbCNNbdtJAF3Z1ZlWpm1jRy8aXtjt0HuG7NFg4cOVYx7AEOD49w4Mgxlq/Zwo7dB6a4QjOz5pY68CUVJP2fpEcqLJsu6X5Jg5K2SurJssijxRFuWPskR45VDvpyR46V2h8tpmtvZpYHEznCvxV4fpxlnwHejIhzge8C36q3sLE27NzLsZHRCa1zbGSUjTtfz7IMM7OWlirwJc0DPgbcPU6TZcC6ZHo9cKUk1V9eyerNL497Gmc8h4dHWL15MKsSzMxaXtoj/H8H/g4Y7zB7LrAbICKKwEHg9LqrA0ZGg5f2Happ3Rf3HapvyKaZ2UmkauBLugbYFxHb6v0wSSslDUgaGBoaSrXO4eEibdNq+2OhbZo4PFysaV0zs5NNmiP8y4Clkl4F7gOukPSfZW1eA7oBJLUB7wX2l28oItZERG9E9HZ1pbv3z8yONoo1HqUXR4OZHS13qYGZ2aSoGvgR8Q8RMS8ieoDlwKMR8TdlzfqBG5Lpa5M2mZxLKUwT551R2xWv558xi0KNfx2YmZ1sah6HL+kOSX+8u9g9wOmSBoGvAF/Lorg/WtV3DjM7Jnb17MyOAqv6zs2yDDOzljahwI+IzRFxTTJ9e0T0J9PvRMQnI+LciFgYEa+ceEsTs+SCObQXJva7qb0wjasvODPLMszMWlpLXGk7va3AupsWMqM93VH+jPZS+wnfU8fM7CTWEoEPsKC7k/tWLqJzRvu4p3dmdhTonNHOfSsX+V46ZmZlWmoIy4LuTrb+45Vs3Pk6qzcP8uK77pZ5Gqv6zuHqC870kb2ZWQUtFfhQOr3z8Yvn8vGL5zIyGhweLjKzo82jcczMqmi5wB+rME2855T2RpdhZtYSWuYcvpmZ1ceBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznEjzEPNTJD0paYekXZK+UaHNCklDkrYnr89OTrlmZlarNDdPOwpcERGHJLUDT0jaGBFbytrdHxG3ZF+imZlloWrgJw8jP5S8bU9emTyg3MzMpk6qc/iSCpK2A/uATRGxtUKzT0h6RtJ6Sd3jbGelpAFJA0NDQ3WUbWZmE5Uq8CNiJCIuAuYBCyV9qKzJw0BPRFwIbALWjbOdNRHRGxG9XV1d9dRtZmYTNKFROhFxAHgMWFw2f39EHE3e3g18OJvyzMwsK2lG6XRJ6kymZwBXAb8uazNnzNulwPNZFmlmZvVLM0pnDrBOUoHSL4gHIuIRSXcAAxHRD3xJ0lKgCLwBrJisgs3MrDYqDcKZer29vTEwMNCQzzYza1WStkVEby3r+kpbM7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTjjwzcxyIs0Tr06R9KSkHZJ2SfpGhTbTJd0vaVDSVkk9k1GsmZnVLs0R/lHgiohYAFwELJa0qKzNZ4A3I+Jc4LvAt7It08zM6lU18KPkUPK2PXmVPyZrGbAumV4PXClJmVVpZmZ1S3UOX1JB0nZgH7ApIraWNZkL7AaIiCJwEDi9wnZWShqQNDA0NFRf5WZmNiGpAj8iRiLiImAesFDSh2r5sIhYExG9EdHb1dVVyybMzKxGExqlExEHgMeAxWWLXgO6ASS1Ae8F9mdRoJmZZSPNKJ0uSZ3J9AzgKuDXZc36gRuS6WuBRyOi/Dy/mZk1UFuKNnOAdZIKlH5BPBARj0i6AxiIiH7gHuBHkgaBN4Dlk1axmZnVpGrgR8QzwMUV5t8+Zvod4JPZlmZmZlnylbZmZjnhwDczywkHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeVEmidedUt6TNJzknZJurVCmz5JByVtT163V9qWmZk1TponXhWBr0bE05JOA7ZJ2hQRz5W1+2VEXJN9iWZmloWqR/gRsTcink6m/wA8D8yd7MLMzCxbEzqHL6mH0uMOt1ZYfKmkHZI2SvpgBrWZmVmG0pzSAUDSLOBB4MsR8VbZ4qeBsyLikKQlwE+A8ypsYyWwEmD+/Pk1F21mZhOX6ghfUjulsL83Ih4qXx4Rb0XEoWR6A9AuaXaFdmsiojcieru6uuos3czMJiLNKB0B9wDPR8R3xmlzZtIOSQuT7e7PslAzM6tPmlM6lwHXAzslbU/m3QbMB4iIu4BrgVWSisARYHlExCTUa2ZmNaoa+BHxBKAqbe4E7syqKDMzy56vtDUzywkHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOZHmEYfdkh6T9JykXZJurdBGkr4naVDSM5IumZxyzcysVmkecVgEvhoRT0s6DdgmaVNEPDemzdXAecnrr4DVyb9mZtYkqh7hR8TeiHg6mf4D8Dwwt6zZMuCHUbIF6JQ0J/NqzcysZhM6hy+pB7gY2Fq2aC6we8z7PRz/SwFJKyUNSBoYGhqaWKVmZlaX1IEvaRbwIPDliHirlg+LiDUR0RsRvV1dXbVswszMapQq8CW1Uwr7eyPioQpNXgO6x7yfl8wzM7MmkWaUjoB7gOcj4jvjNOsHPp2M1lkEHIyIvRnWaWZmdUozSucy4Hpgp6TtybzbgPkAEXEXsAFYAgwCbwM3Zl+qmZnVo2rgR8QTgKq0CeDmrIoyM7Ps+UpbM7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTjjwzcxyIs0Tr9ZK2ifp2XGW90k6KGl78ro9+zLNzKxeaZ549QPgTuCHJ2jzy4i4JpOKzMxsUlQ9wo+IXwBvTEEtZmY2ibI6h3+ppB2SNkr64HiNJK2UNCBpYGhoKKOPNjOzNLII/KeBsyJiAfB94CfjNYyINRHRGxG9XV1dGXy0mZmlVXfgR8RbEXEomd4AtEuaXXdlZmaWqboDX9KZkpRML0y2ub/e7ZqZWbaqjtKR9GOgD5gtaQ/wdaAdICLuAq4FVkkqAkeA5RERk1axmZnVpGrgR8R1VZbfSWnYppmZNTFfaWtmlhMOfDOznHDgm5nlhAPfzCwnHPhmZjnhwDczywkHvplZTjjwzcxywoFvZpYTDnwzs5xw4JuZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU5UDXxJayXtk/TsOMsl6XuSBiU9I+mS7Ms0M7N6pTnC/wGw+ATLrwbOS14rgdX1l2VmZlmrGvgR8QvgjRM0WQb8MEq2AJ2S5mRVoJmZZSOLc/hzgd1j3u9J5pmZWROZ0i9tJa2UNCBpYGhoaCo/2sws97II/NeA7jHv5yXzjhMRayKiNyJ6u7q6MvhoMzNLK4vA7wc+nYzWWQQcjIi9GWzXzMwy1FatgaQfA33AbEl7gK8D7QARcRewAVgCDAJvAzdOVrFmZla7qoEfEddVWR7AzZlVZGZmk8JX2pqZ5YQD38wsJxz4ZmY54cA3M8sJB76ZWU448M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLCQe+mVlOOPDNzHLCgW9mlhMOfDOznEgV+JIWS3pB0qCkr1VYvkLSkKTtyeuz2ZdqZmb1SPOIwwLwH8BVwB7gKUn9EfFcWdP7I+KWSajRzMwykOYIfyEwGBGvRMQwcB+wbHLLMjOzrKUJ/LnA7jHv9yTzyn1C0jOS1kvqrrQhSSslDUgaGBoaqqFcMzOrVVZf2j4M9ETEhcAmYF2lRhGxJiJ6I6K3q6sro482M7M00gT+a8DYI/Z5ybw/iYj9EXE0eXs38OFsyjMzs6ykCfyngPMknS2pA1gO9I9tIGnOmLdLgeezK9HMzLJQdZRORBQl3QL8DCgAayNil6Q7gIGI6Ae+JGkpUATeAFZMYs1mZlYDRURDPri3tzcGBgYa8tlmZq1K0raI6K1lXV9pu2cPfPGLcOmlcOqpIMGrrza6KjOz461fz3/DOUi/QTqC9ALSvyGdlmZ1B/7gIDzwALzvffCRjzS6GjOz8X3724xAALcBi4HVwCpgE1LVPK96Dv+k99GPwu9+V5q++274+c8bW4+Z2XgefpiPnXHGKxFxbzLncaQ3KA2F7wMePdHqPsKf5h+BmbWIytcvPZX8W+mC2Hdx2pmZtbbLk3+rDod34JuZtSppLnAH8D9EVB326MA3M2tF0izgp5Suf7oxzSr+0tbMrNVIMyjdw+z9wOVE7EmzmgPfzKyFdICA9UAvcBURO9Ou68A3M2sVo6P8F5wNfAC4hogtE1m9pQO/ODLK28dGmNnRRmGaat/Q+vWlf7dtK/27cWNp+FNXF1x++fjrmZmllEle3XwzS+B9wDeBw0iLxizdU+3UTsvdS+docYQNO/eyevPLvLTvEG3TRHE0OP+MWXyh7xyWXDCH6W2FiW1U4/zwL78cNm+ecI1mZjAJedXTA7/5zXhLv0HEP59o9ZYK/O27D7Bi7ZMcGxnl8PDIcctndhRoL0xj3U0LWdDdmVWpZmYTNll5lYubp+3YfYDr1mzhwJFjFX94AIeHRzhw5BjL12xhx+4DU1yhmVlJs+ZVSwT+0eIIN6x9kiPHKv/gyh05Vmp/tJiuvZlZVpo5r1IFvqTFkl6QNCjpaxWWT5d0f7J8q6SeLIvcsHMvx0ZGJ7TOsZFRNu58PcsyzMyqaua8qhr4kgrAfwBXUxoKdJ2kD5Q1+wzwZkScC3wX+FaWRa7e/PK4fxaN5/DwCKs3D2ZZhplZVc2cV2mO8BcCgxHxSkQMA/cBy8raLKN0e04oXRBwpTTe0JeJGRkNXtp3qKZ1X9x3iJHRxnwpbWb50+x5lSbw5wK7x7zfw/G34fxTm4goAgeB08s3JGmlpAFJA0NDQ6kKPDxcpK3GMatt08Th4WJN65qZTVSz59WUfmkbEWsiojciersq39f5ODM72ijW+FuvOBrM7Gjpa8vMrIU0e16lCfzXgO4x7+cl8yq2kdQGvBfYn0WBhWnivDNm1bTu+WfMqu8KXDOzCWj2vEoT+E8B50k6W1IHsBzoL2vTD9yQTF8LPBoZXtG1qu8cZnZM7OrZmR0FVvWdm1UJZmapNHNeVQ385Jz8LcDPKD1R5YGI2CXpDklLk2b3AKdLGgS+Ahw3dLMeSy6YQ3thYmef2gvTuPqCM7Msw8ysqmbOq1RVRcSGiDg/Is6JiG8m826PiP5k+p2I+GREnBsRCyPilSyLnN5WYN1NC5nRnu635oz2UvsJ31PHzKxOzZxXLXGlLcCC7k7uW7mIzhnt4/65NLOjQOeMdu5bucj30jGzhmnWvGqpm6dB6bLljTtfZ/XmQV58193nTmNV3zlcfcGZPrI3s6YwGXlVz83TWi7wxxoZDQ4PF+u/H76Z2STLKq/qCfyWHqRemCbec0p7o8swM6uqGfKqZc7hm5lZfRz4ZmY50bBz+JKGgHGf1TUBs4HfZ7CdRnM/mov70Vzcjz87KyLS3ZumTMMCPyuSBmr9AqOZuB/Nxf1oLu5HNnxKx8wsJxz4ZmY5cTIE/ppGF5AR96O5uB/Nxf3IQMufwzczs3ROhiN8MzNLoSUCX9JiSS9IGpR03K2XJU2XdH+yfKuknqmvsroU/VghaUjS9uT12UbUWY2ktZL2SXp2nOWS9L2kn89IumSqa0wjRT/6JB0csz9un+oaq5HULekxSc9J2iXp1gptmn5/pOxH0+8PAEmnSHpS0o6kL9+o0KYxmRURTf0CCsDLwPuBDmAH8IGyNn8L3JVMLwfub3TdNfZjBXBno2tN0ZePApcAz46zfAmwERCwCNja6Jpr7Ecf8Eij66zShznAJcn0acCLFf67avr9kbIfTb8/kjoFzEqm24GtwKKyNg3JrFY4wl8IDEbEKxExDNwHLCtrswxYl0yvB66U1Gx3U0vTj5YQEb8A3jhBk2XAD6NkC9Apac7UVJdein40vYjYGxFPJ9N/oPSQorllzZp+f6TsR0tIfs6Hkrftyav8y9KGZFYrBP5cYPeY93s4/j+EP7WJ0hO6DgKnT0l16aXpB8Ankj+710vqrrC8FaTtayu4NPnTfKOkDza6mBNJTgtcTOmIcqyW2h8n6Ae0yP6QVJC0HdgHbIqIcffJVGZWKwR+njwM9ETEhcAm/nwEYI3xNKXL2BcA3wd+0uB6xiVpFvAg8OWIeKvR9dSqSj9aZn9ExEhEXATMAxZK+lCja4LWCPzXgLFHuvOSeRXbSGoD3gvsn5Lq0qvaj4jYHxFHk7d3Ax+eotqylmafNb2IeOuPf5pHxAagXdLsBpd1HEntlELy3oh4qEKTltgf1frRKvtjrIg4ADwGLC5b1JDMaoXAfwo4T9LZkjoofcHRX9amH7ghmb4WeDSSb0OaSNV+lJ1XXUrpPGYr6gc+nYwOWQQcjIi9jS5qoiSd+cfzqpIWUvr/pakOJJL67gGej4jvjNOs6fdHmn60wv4AkNQlqTOZngFcBfy6rFlDMqvpH4ASEUVJtwA/ozTSZW1E7JJ0BzAQpQep3wP8SNIgpS/hljeu4spS9uNLkpYCRUr9WNGwgk9A0o8pjZiYLWkP8HVKX0wREXcBGyiNDBkE3gZubEylJ5aiH9cCqyQVgSPA8iY8kLgMuB7YmZwzBrgNmA8ttT/S9KMV9geURhytk1Sg9EvpgYh4pBkyy1fampnlRCuc0jEzsww48M3McsKBb2aWEw58M7OccOCbmeWEA9/MLCcc+GZmOeHANzPLif8HMaV9ILW+Ur4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scripts import utilities\n",
    "from scripts import plots\n",
    "import matplotlib.pyplot as plt\n",
    "# cities = utilities.create_cities(3)\n",
    "cities = np.array([[0, 4],[0, 0],[3, 0]])\n",
    "distance_matrix = utilities.get_distance_matrix(cities)\n",
    "plots.plot_cities(cities)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: Fast method for expectation will be used. Noise\n",
      "                     models will be ineffective\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2225.902046395591\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2959.7191426814898\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2418.5086486322734\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2791.1897207746892\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2817.2141714650375\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2506.2053350309484\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2746.412355359732\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2461.5320208204657\n",
      "\tParameters: [0.32120728 5.81929383] \n",
      "\tE => -2917.102934101301\n",
      "\tParameters: [0.32171701 5.81900969] \n",
      "\tE => -2989.5439816908374\n",
      "\tParameters: [0.32071324 5.81900969] \n",
      "\tE => -2964.4265003113405\n",
      "\tParameters: [0.3212112  5.81915176] \n",
      "\tE => -2991.0777956864486\n",
      "\tParameters: [0.32020743 5.81915176] \n",
      "\tE => -2983.0174931466145\n",
      "\tParameters: [0.32071128 5.81908073] \n",
      "\tE => -2993.3252511963615\n",
      "\tParameters: [0.31970751 5.81908073] \n",
      "\tE => -2992.248956828172\n",
      "\tParameters: [0.31970751 5.81908073] \n",
      "\tE => -2992.9691102365696\n",
      "\tParameters: [0.31970751 5.81908073] \n",
      "\tE => -2993.3415997381435\n",
      "\tParameters: [0.31920513 5.81909848] \n",
      "\tE => -2993.240116895994\n",
      "\tParameters: [0.31870373 5.81908073] \n",
      "\tE => -2993.246348321256\n",
      "\tParameters: [0.31820136 5.81909848] \n",
      "\tE => -2993.250733394864\n",
      "\tParameters: [0.31769996 5.81908073] \n",
      "\tE => -2993.253492215863\n",
      "\tParameters: [0.31719758 5.81909848] \n",
      "\tE => -2993.2531862802393\n",
      "The solution is: [0, 0, 2]\n"
     ]
    }
   ],
   "source": [
    "tsp_solver = ForestTSPSolverNaive(distance_matrix, use_constraints=True)\n",
    "solution, naive_distribution = tsp_solver.solve_tsp()\n",
    "print(\"The solution is:\", solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wait a minute... This doesn't make sense at all. We wanted to get `[0, 1, 2]`, and I got `[0, 1, 2, 2]`. \n",
    "(Since the result is probabilistic, you might have gotten something different. Also every time I update this notebook I get something different, so sorry for some inconsistencies here.)\n",
    "<br>\n",
    "Maybe something is wrong with the `binary_state_to_points_order` in our solver? We should investigate the binary string."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Your binary solution is: (1, 0, 0, 1, 0, 0, 0, 0, 1)\n",
      "My initial binary solution is: (1, 1, 1, 0, 0, 0, 0, 0, 1)\n"
     ]
    }
   ],
   "source": [
    "print(\"Your binary solution is:\", tsp_solver.most_frequent_string)\n",
    "print(\"My initial binary solution is:\", (1,1,1,0,0,0,0,0,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This solution doesn't make sense at all! This string means, that at t=0 we are in all the cities at the same time, then at t=1 we are not present anywhere and at t=2 we are in the city 2.\n",
    "<br>\n",
    "It's not a surprise, that our order of traveling looks weird."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point you may have noticed something suspicious with our class. It's called `ForestTSPSolverNaive`.\n",
    "<br>\n",
    "\"Naive\"... It doesn't sound very promising to me. Let's see why I named it naive, what's actually inside it and why we got such a crappy solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating cost operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you remember from the previous part, in order to make QAOA solve the problem we want, we need to translate it into its native language - language of Pauli operators.\n",
    "<br>\n",
    "What we need to do for the TSP is pretty similar to what we have done for MaxCut, so let's code it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0. 4. 5.]\n",
      " [4. 0. 3.]\n",
      " [5. 3. 0.]]\n"
     ]
    }
   ],
   "source": [
    "distance_matrix = utilities.get_distance_matrix(cities)\n",
    "print(distance_matrix)\n",
    "number_of_cities = len(cities)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last time we learned two useful equations for MaxCut which will guide us now. \n",
    "<br>\n",
    "One for the classical cost function:\n",
    "<br>\n",
    "$C_{total} = \\sum C_{ij} = \\sum \\frac{1}{2} w_{ij} (1 - z_i z_j) $,\n",
    "\n",
    "And one showing how to construct the quantum cost operator:\n",
    "<br>\n",
    "$H_{cost} = \\sum \\frac{1}{2} w_{ij} (\\mathbb{1} - \\sigma^z_i \\sigma^z_j$)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, for TSP it's a little bit more complicated, since $i$ and $j$ correspond to indices of qubits, not the indices of cities.\n",
    "<br>\n",
    "I tried to write a nice mathematical formula for you, but after several minutes I have not produced anything concise and elegant.\n",
    "I guess that I'm so spoiled with the programming that it's much more intuitive for me to just show you how it works in code ;)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "City 0 to 1 at t = 0 costs -2.0 Qubits: 0 4\n",
      "City 0 to 1 at t = 1 costs -2.0 Qubits: 3 7\n",
      "City 0 to 2 at t = 0 costs -2.5 Qubits: 0 5\n",
      "City 0 to 2 at t = 1 costs -2.5 Qubits: 3 8\n",
      "City 1 to 2 at t = 0 costs -1.5 Qubits: 1 5\n",
      "City 1 to 2 at t = 1 costs -1.5 Qubits: 4 8\n"
     ]
    }
   ],
   "source": [
    "from pyquil.paulis import PauliTerm\n",
    "\n",
    "cost_operators = []\n",
    "number_of_cities = len(distance_matrix)\n",
    "\n",
    "for i in range(number_of_cities):\n",
    "    for j in range(i, number_of_cities):\n",
    "        for t in range(number_of_cities - 1):\n",
    "            weight = -distance_matrix[i][j]/2\n",
    "            if distance_matrix[i][j] != 0:\n",
    "                qubit_1 = t * number_of_cities + i\n",
    "                qubit_2 = (t + 1) * number_of_cities + j\n",
    "                cost_operators.append(PauliTerm(\"I\", 0, weight) - PauliTerm(\"Z\", qubit_1, weight) * PauliTerm(\"Z\", qubit_2))\n",
    "                print(\"City\", i, \"to\", j, \"at t =\", t, \"costs\", weight, \"Qubits:\", qubit_1, qubit_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several things going on here:\n",
    "1. We have \"-\" sign before Z terms - that's because we want to get the shortest, not the longest path.\n",
    "2. In principle we could omit the $\\frac{1}{2}$ term. Since constants don't influence the optimization process, we can do that, though it's good to be consistent with the equations, to avoid confusion, so we finally included it.\n",
    "3. We need to include every cost (N-1) times - once for each step between the cities. 3 time slots -> 2 steps.\n",
    "4. It's hard to see right away why we use these qubits and not the others. Take your time to understand where are those two lines defining `qubit_1` and `qubit_2` coming from. Pen and paper is useful for that :)\n",
    "5. We ommited terms with zeros - well, it didn't cost us anything ;)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this is how the cost operators look like. As you can see, the `Z` terms correspond to the qubit indices above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-2+0j)*I + (2-0j)*Z0*Z4\n",
      "(-2+0j)*I + (2-0j)*Z3*Z7\n",
      "(-2.5+0j)*I + (2.5-0j)*Z0*Z5\n",
      "(-2.5+0j)*I + (2.5-0j)*Z3*Z8\n",
      "(-1.5+0j)*I + (1.5-0j)*Z1*Z5\n",
      "(-1.5+0j)*I + (1.5-0j)*Z4*Z8\n"
     ]
    }
   ],
   "source": [
    "for single_cost_operator in cost_operators:\n",
    "    print(single_cost_operator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so this is basically the code we used (you can check out `forest_tsp_solver_naive.py` in `scripts` directory). \n",
    "But it doesn't answer the question, why the results we got don't make sense at all.\n",
    "\n",
    "Take your time to solve this puzzle by yourself and when you are ready, go to the next section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I leave some blank space, so you can try to solve this mystery by yourself, without the next header spoiling it.\n",
    "<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>\n",
    "<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spoiler-free header (just in case)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the main reasons why it didn't work is that we forgot about the constraints.\n",
    "<br>\n",
    "Constraints in TSP are pretty important and there are two of them:\n",
    "- we need to visit each city once\n",
    "- in every moment we need to be in exactly one city\n",
    "\n",
    "It wasn't immediately obvious to me how to encode that into our problem, so I used a standard optimization approach - adding a penalty."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is a penalty?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say we want to find the minimum of function $x^2$, but we want only values higher than 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl0VPX9//HnmySQBUgIYUsCBARBViEBwaXFrW5Fca1UBBXF1lZttbVa+6vt6SJ+rdalVYsKBRdQqRZb0aqoRVsLhB1JIBiWBBKSEBJCQvbP749MPWhRIJPJneX1OMczMzczua+B+OLmM/d+PuacQ0REwlcHrwOIiEhgqehFRMKcil5EJMyp6EVEwpyKXkQkzKnoRUTCnIpeRCTMqehFRMKcil5EJMxFex0AICUlxWVkZHgdQ0QkpKxevbrMOdfjaM8LiqLPyMggOzvb6xgiIiHFzHYey/M0dCMiEuZU9CIiYU5FLyIS5lT0IiJhTkUvIhLmjlr0ZjbXzErMbNNh25LN7B0zy/PddvNtNzN7zMy2mdkGMxsbyPAiInJ0x3JE/2fg/C9suxtY5pwbDCzzPQa4ABjs+28W8GTbxBQRkdY6atE755YD5V/YfAkw33d/PjDlsO0LXIv/AElm1qetwoqIhIvmZsdv3thMQXlNwPfV2jH6Xs65It/9YqCX734aUHDY8wp92/6Hmc0ys2wzyy4tLW1ljMgzadIkJk2a5HUMEfHT4+9t4+kPt/NhXlnA9+X3h7GuZXXx415h3Dk3xzmX5ZzL6tHjqFfwioiEjX9uLeWRZVu5bEwaU8f3Dfj+Wlv0e/87JOO7LfFt3w0cnjrdt01ERIDC/TXcvmgtQ3p14TeXjsTMAr7P1hb968AM3/0ZwJLDtk/3nX0zAag8bIhHRCSi1TU2ccsLa2hqcjw5LZO4jlHtst+jTmpmZguBSUCKmRUC9wGzgZfNbCawE7jK9/SlwIXANqAGuD4AmUVEQtIvXt/MhsJK/nRtJgNSEtptv0cteufc1C/50tlHeK4DvudvKBGRcLNo5S4WrtzFdyedwHnDe7frvnVlrIhIgK0rqODnSz7hjMEp/OgbQ9p9/yp6EZEAKjtYx3efX02PLp147OoxRHUI/IevXxQUC4+IiISjxqZmvv/iGsqr6/nLd0+lW0JHT3Ko6EVEAuSBt3L5T345D181mhFpiZ7l0NCNiEgA/G39Hp7+cDszJvbnsrHpnmZR0YuItLHc4gPctXgDWf27ce9Fw7yOo6IXEWlLlYcauPm51XSJjeaJa8bSMdr7mtUYvYhIG2ludvzwpXXsqTjEolkT6Nk11utIgI7oRUTazGPv5fFebgk//+YwMvsnex3nMyp6EZE2sCxnL4+8m8cVmelMm9Df6zifo6IXEfHTtpIqbl+0jhFpXfn1lBHtMiPl8VDRi4j4obKmgRvnZxMbE8Wca7OIjWmfGSmPh4peRKSVGpua+f7CNeyuOMRT08aSmhTndaQj0lk3IiKt9MBbuXyYV8YDl48kKyN4Pnz9Ih3Ri4i0wl9WF3525eu3xvXzOs5XUtGLiByndQUV3PPaRiYO7M7Pvun9la9Ho6IXETkOew/UMmtBNr26duKJa8YSExX8NRr8CUVEgkRtQxOznlvNwbpGnp6e5dm0w8dLH8aKiBwD5xz3vraJ9QUVPDUtk6G9u3od6ZjpiF5E5Bg8+9F2/rKmkB+ecyLnj2jfNV/9paIXETmK5VtL+e3SHC4Y0ZtbzxrkdZzjpqIXEfkK20oO8r0X13Biry787srRdPBgzVd/qehFRL7E/up6Zs5fRafoDjwzI4uETqH5sWZophYRCbD6xmZufn41RZW1LJo1gfRu8V5HajUd0YuIfEHLGTYbWbm9nAevGMXYft28juQXFb2IyBfMWZ7PK6sLue3swVxycprXcfymohcROcw/Pilm9lu5XDSqDz84e7DXcdqEil5ExGfT7kp+sGgdo9KTeChEz7A5EhW9iAhQcqCWmxZkkxQfw9PXZgblAiKtpbNuRCTiHapv4qYF2VQeauCV70ykZ9dYryO1Kb+O6M3sh2b2iZltMrOFZhZrZgPMbIWZbTOzl8wsNGb9EZGI1NzsuPOVdWzYXcmjV49heGqi15HaXKuL3szSgNuALOfcCCAKuBp4APi9c24QsB+Y2RZBRUQC4eF3trJ0YzH3XDCUc4f18jpOQPg7Rh8NxJlZNBAPFAFnAYt9X58PTPFzHyIiAfHSql384f1tXD2uLzedMdDrOAHT6qJ3zu0GfgfsoqXgK4HVQIVzrtH3tEIg9E9CFZGws3xrKT99bRNnDE7hV1NGYBYeZ9gciT9DN92AS4ABQCqQAJx/HK+fZWbZZpZdWlra2hgiIsdt854D3PLCGgb37Bwyq0T5w593dw6w3TlX6pxrAF4FTgOSfEM5AOnA7iO92Dk3xzmX5ZzL6tGjhx8xRESOXXFlLTf8eRWdO0Uz7/pxdImN8TpSwPlT9LuACWYWby2/85wNbAbeB67wPWcGsMS/iCIibaOqtoHr/7yKqtoG5l43jj6JcV5Hahf+jNGvoOVD1zXARt/3mgP8BLjDzLYB3YFn2yCniIhfGpqa+d6La9m6t4onpmUyLDV0lgL0l18XTDnn7gPu+8LmfGC8P99XRKQtOef4+ZJNLN9ayv2XjeTrJ0bWcHF4fwIhIgI88cGnLFxZwC2TTmDq+H5ex2l3KnoRCWtL1u3mwX9s4eLRqfzoG0O8juMJFb2IhK0V+fv48SsbGJ+RzINXjgqb2SiPl4peRMJSbvEBblyQTXpyHHOmZ9IpOnxmozxeKnoRCTu7Kw4xY+5K4jtGseCG8STFR/bcipqmWETCyv7qeqY/u4Ka+iZe+c7EkF7Uu63oiF5Ewsah+iZmzl9FQfkhnp6exdDekXOu/FfREb2IhIXGpmZuXbiWtQUV/PHbY5kwsLvXkYKGjuhFJOQ55/h/Szbxbs5efjF5OBeO7ON1pKCioheRkPf7d/NYuLKA7515AjNOzfA6TtBR0YtISHv+Pzt5bFkeV2WlR+wFUUejoheRkPXWpmJ+vmQTZw3tyW8vHRnWi4f4Q0UvIiFpRf4+blu0llHpSfzh22OIDvPFQ/yhPxkRCTmbdldy4/xs+naLY+5144jvqBMIv4qKXkRCyraSg0yfu5KucTE8f+MpJCdE9lWvx0JFLyIhY3fFIaY/u4IOBs/NHB8xK0T5S7/viEhIKDtYx7XPrKCqrpFFsyYwsEdnryOFDB3Ri0jQO1DbwIy5K9lTeYi5141jeGqi15FCiopeRIJabUMTN87PZktxFU9Oy2RcRrLXkUKOhm5EJGg1NDVzywtrWLWjnMeuHsOZQ3p6HSkk6YheRIJSc7PjR6+s573cEn49ZQSTR6d6HSlkqehFJOg45/jF3z5hybo93HX+EK45pb/XkUKail5Egopzjtlv5bLg453c/LWBfPfrJ3gdKeSp6EUkqDy6LI8//TOfayf05+4Lhmr+mjagoheRoPHUPz/lkXfzuDIznV9ePFwl30ZU9CISFOb9azuz38zl4tGpzL58FB06qOTbiopeRDy3cOUufvm3zZw3vBcPXTWaKJV8m1LRi4inXltbyE9f28iZQ3rw+NSxxGi64TanP1ER8cwbG4q48+X1TBzYnSenZdIxWpUUCPpTFRFPvLN5L7cvWktm/248MyOL2JgoryOFLRW9iLS75VtL+d4Laxie2lULh7QDv4rezJLMbLGZ5ZpZjplNNLNkM3vHzPJ8t93aKqyIhL7lW0u5aUE2J/TszPwbxtMlNsbrSGHP3yP6R4G3nHNDgdFADnA3sMw5NxhY5nssIvJZyQ/s0ZkXbjyFpHitDtUeWl30ZpYIfA14FsA5V++cqwAuAeb7njYfmOJvSBEJfR/mfb7ktQRg+/HniH4AUArMM7O1ZvaMmSUAvZxzRb7nFAO9jvRiM5tlZtlmll1aWupHDBEJdh/mlXLjfJW8V/wp+mhgLPCkc24MUM0Xhmmccw5wR3qxc26Ocy7LOZfVo0cPP2KISDD7b8kPSElQyXvEn6IvBAqdcyt8jxfTUvx7zawPgO+2xL+IIhKqPsor+6zkX7xpgkreI60ueudcMVBgZkN8m84GNgOvAzN822YAS/xKKCIh6aO8MmbOX6WSDwL+nrx6K/CCmXUE8oHrafnH42UzmwnsBK7ycx8iEmJU8sHFr6J3zq0Dso7wpbP9+b4iErr+ubWUWQs0XBNMdGWsiLSZdzbv5ab52ZzQo7NKPojoumMRaRNvbCji9kVrGZGWyPzrx5MYryteg4WO6EXEb6+uKeTWhWsY0y+J52aq5IONjuhFxC8vrtjFvX/dyKkndOfp6VmaoCwI6W9ERFpt3r+288u/bebMIT14clqmphoOUip6EWmVJz/4lAfeyuW84b14fOpYLRoSxFT0InJcnHM88m4ejy7L4+LRqTx01Wgt/xfkVPQicsycc8x+M5c/Lc/nysx0Zl8+Sgt5hwAVvYgck6Zmx72vbWTRqgKundCfX148nA4q+ZCgoheRo6prbOKHL61j6cZibj1rEHeceyJmKvlQoaIXka9UXdfId55fzYd5ZfzsopO48YyBXkeS46SiF5EvVVFTz3XzVrGhsIIHrxjFlVl9vY4kraCiF5Ej2nuglunPrmR7WTVPTsvkvOG9vY4kraSiF5H/sXNfNdOeXUH5wXrmXT+O0waleB1J/KCiF5HPySk6wPS5K2lsaubFmyYwum+S15HET7rKQUQ+s2pHOd/608dEmfHyzRNV8mFCR/QiAsBbm4q4bdE60pPimH/DePomx3sdSdqIil5EWPDxDu57/RNO7pvEszPGacGQMKOiF4lgzjke/McWnvjgU845qSePTx1LXEfNQBluVPQiEaqhqZmf/GUDr67ZzdTxffnVJSOI1uRkYUlFLxKBDtY1cssLa1i+tZQ7zj2RW88apCkNwpiKXiTClFbVccOfV7G56AAPXD6Sb43r53UkCTAVvUgEyS89yHXzVlFaVcfT0zM5a2gvryNJO1DRi0SIFfn7uPn51XQwY+GsCZysc+QjhopeJAK8traQuxZvoG9yPPOuG0f/7gleR5J2pKIXCWOHL/s3cWB3npqWSWJ8jNexpJ2p6EXCVF1jEz9ZvIG/rtvDFZnp/PbSkVrAO0Kp6EXC0P7qem5+bjUrd5Tz4/OGcMukE3T6ZART0YuEme1l1Vw/byV7Kmt5fOoYJo9O9TqSeExFLxJGPndmzU2nkNk/2etIEgT8HrAzsygzW2tmf/c9HmBmK8xsm5m9ZGaaHUmkHSxcuYtrnllBckJHXrvlVJW8fKYtPpm5Hcg57PEDwO+dc4OA/cDMNtiHiHyJxqZm7luyiXte3chpg1J47ZbTdPqkfI5fRW9m6cBFwDO+xwacBSz2PWU+MMWffYjIl6uoqWfGvJXM/3gnN54+gLnXjSMxTqdPyuf5O0b/CHAX0MX3uDtQ4Zxr9D0uBNL83IeIHMG2kipunJ/NnopaHrxiFFdm9fU6kgSpVh/Rm9k3gRLn3OpWvn6WmWWbWXZpaWlrY4hEpPdzS5jyx39zsK6JhbNOUcnLV/LniP404GIzuxCIBboCjwJJZhbtO6pPB3Yf6cXOuTnAHICsrCznRw6RiOGcY87yfGa/lcvw1K7MuTaL1KQ4r2NJkGv1Eb1z7h7nXLpzLgO4GnjPOXcN8D5whe9pM4AlfqcUEWobmrjz5fXc/2YuF47swys3n6qSl2MSiPPofwIsMrNfA2uBZwOwD5GIUlBew83PrSan+AB3nnsi39dCIXIc2qTonXMfAB/47ucD49vi+4oIfLClhNsXrcM5x9zrxnHmkJ5eR5IQoytjRYJUc7Pjj+9v4+F3tzKkVxf+dG2mzo+XVlHRiwShA7UN3PHSet7N2cuUk1O5/7JRxHWM8jqWhCgVvUiQ2VJcxXeeX01BeQ2/mDyMGadmaDxe/KKiFwkif9+wh7sWbyChUzQLZ01gXIbmqxH/qehFgkBdYxP3L83lz//eQVb/bjxxzVh6do31OpaECRW9iMd27avh+wvXsKGwkhtPH8Bd5w/VSlDSplT0Ih56a1MxP168HgPmXJvJN4b39jqShCEVvYgH6hubuf/NHOb9awej0xP5w7fH0jc53utYEqZU9CLtrKC8hu8vXMv6ggquPy2Dey44SUM1ElAqepF29PYnxfzolfU44KlpYzl/RB+vI0kEUNGLtIPahiZ+uzSHBR/vZGRaIn/89lj6dddQjbQPFb1IgG0pruK2hWvZsreKmacP4K7zh9ApWle5SvtR0YsEiHOO5/6zk1+/kUPX2Bjm3zCer5/Yw+tYEoFU9CIBsO9gHXct3sCy3BImDenB764cTUrnTl7HkgilohdpYx/mlXLHy+uprGngvsnDuE5z1YjHVPQibaS2oYmH39nKnOX5DOrZmfnXj2dYalevY4mo6EXawid7KrnjpfVs2VvFNaf042cXDdO0whI0VPQifmhsauapf37KI+/mkZzQkXnXjePMoVoBSoKLil6klT4tPcidL69nXUEFk0en8qtLhpMU39HrWCL/Q0Uvcpyamx0LPt7B7LdyiY2J4vGpY5g8OtXrWCJfSkUvchx2VxzirsXr+de2fZw5pAcPXD5K88ZL0FPRixyD5mbHwlW7uH9pLs457r9sJFeP66vTJiUkqOhFjmJHWTV3v7qB/+SXc9qg7sy+bJSmFJaQoqIX+RJNzY5nP8rnobe30jGqA7MvG8m3dBQvIUhFL3IEW4qruGvxetYXVnLOSb349ZQR9E7UWLyEJhW9yGHqG5t54oNt/PH9bXSJjeGxqWOYPKqPjuIlpKnoRXxW7Sjn3tc2snXvQS4encp9k4fRXRORSRhQ0UvEK6+uZ/abObycXUhaUhzPTM/inGG9vI4l0mZU9BKxnHMsXl3Ib5fmcKC2kZu/NpDbzxlMfEf9byHhRT/REpHy9lZx7183sXJ7OZn9u/GbS0cwtLdmmpTwpKKXiHKovonH38tjzvJ8OsdG88DlI7kysy8dOujDVglfrS56M+sLLAB6AQ6Y45x71MySgZeADGAHcJVzbr//UUVazznHm5uK+c0bOeyuOMQVmencc8FQfdgqEcGfI/pG4E7n3Boz6wKsNrN3gOuAZc652WZ2N3A38BP/o4q0Tm7xAX75+mY+zt/H0N5dWDRrAhMGdvc6lki7aXXRO+eKgCLf/SozywHSgEuASb6nzQc+QEUvHqioqef372zluf/spGtcDL+aMoKp4/oSHdXB62gi7apNxujNLAMYA6wAevn+EQAopmVo50ivmQXMAujXr19bxBABWqYuWLhyFw+9vYXKQw1cc0p/7jj3RLolaK54iUx+F72ZdQb+AvzAOXfg8CsInXPOzNyRXuecmwPMAcjKyjric0SO17+3lfHrN3LYXHSAUwYk84uLh3NSH51NI5HNr6I3sxhaSv4F59yrvs17zayPc67IzPoAJf6GFDmavL1V3P9mLu/llpCWFMcfvj2Gi0Zq6gIR8O+sGwOeBXKccw8f9qXXgRnAbN/tEr8SinyFkgO1/P7drby0qoCETtHcfcFQrjs1g9gYLcwt8l/+HNGfBlwLbDSzdb5tP6Wl4F82s5nATuAq/yKK/K/qukae/jCfOcvzqW9sZvrEDG47ezDJGocX+R/+nHXzEfBlvxef3drvK/JVGpqaWby6kIff2UppVR0XjuzNXecNJSMlwetoIkFLV8ZKSGhudvxtwx5+/85WduyrYWy/JJ6aNpbM/sleRxMJeip6CWrOOd7evJeH397Klr1VDO3dhaenZ3HOST31QavIMVLRS1ByzvFhXhkPvb2F9YWVDExJ4LGpY/jmyD6al0bkOKnoJeisyN/HQ+9sZeX2ctKS4vi/y0dx2dg0XdEq0koqegkKzjk+2lbG48u2sXJHOT26dOKXFw/n6vF96RStUyVF/KGiF08553gvt4TH39vGuoIKeneN5ReTh3H1+H46F16kjajoxRPNzY5/fFLM4+9tY3PRAdK7xfHbS0dyeWaajuBF2piKXtpVXWMTr6/bw5zl+eSVHGRASgIPXjGKKWPSiNEYvEhAqOilXVQeauDFFbuY96/tlFTVMbR3Fx69+mS+OSqVKJ1FIxJQKnoJqML9Ncz9aAcvrdpFdX0Tpw9K4XdXjuaMwSk6D16knajoJSA2FFbwzIfbeWNjEQZMHp3KjWcMYHhqotfRRCKOil7aTF1jE29uLGb+xztYu6uCzp2iueG0DK4/bQCpSXFexxOJWCp68VtR5SFeXLGLhSt3UXawngEpCdw3eRiXZ6bTNTbG63giEU9FL63inGPF9nIWfLyDf3yyl2bnOHtoT6ZPzOD0QSmapkAkiKjo5biUHazjtTW7WbRqF5+WVpMYF8PM0wdw7YT+9E2O9zqeiByBil6OqqnZ8WFeKS+tKuDdnL00NDnG9kvi/y4fxeTRqcR11AVOIsFMRS9fqqC8hsWrC3klu4A9lbV0i49hxsQMvjWuL4N7dfE6nogcIxW9fE5FTT1vbCzir2t3s2rHfszg9EEp3HvRMM4Z1lPTE4iEIBW9UNvQxPu5Jby2djfvbymhockxqGdnfnzeEC45OZX0bhp7FwllKvoI1dDUzL8/3cfSDUUs3VREVW0jPbp0YsbEDKaMSWN4aldduSoSJlT0EaSusYmP8spYurGYd3P2Unmogc6dovnG8F5cOiaNU09I0bwzImFIRR/mDtU3sTyvlDc3FrEsp4Squka6xkZz7rDeXDCiN6cPTtG87yJhTkUfhgr31/B+bgnv5Zbw70/3UdfYTFJ8DBeO7MMFI3tz6gkpdIzWlMAikUJFHwaamh3rCvazLKel3HOLqwDo3z2ea07pz9kn9WT8gGTN9y4SoVT0IcgB+aUH+de2Mj7aVsbHn+7jQG0jUR2McRnduPfCkzjrpJ4MTEnQB6oioqIPFSVVtXz86T7KBp7PocT+nPXQPwFIS4rjghF9OH1wCl87sQeJcZpETEQ+T0UfhJxz5JdVk72jnFU79pO9o5wd+2oA6JA8iNjKXfz0W6dy+qAU+neP11G7iHwlFX0QOFjXyCe7K1lfWEH2jv1k79xPeXU9AMkJHcns342p4/sx8YTu3DrtUgzHtAl3eJxaREKFir6d1dQ3snnPATYUVrJxdyUbCivIL6vGuZavZ3SP56yhPRmX0Y2sjOT/GWc3nEfJRSRUqegDpLGpmZ3lNWwtrmLL3iq27q1iS3EV28uqafZ1da+unRiZlsQlJ6cxMi2RkemJpHTu5G1wEQk7ASl6MzsfeBSIAp5xzs0OxH685pyjvLqeHftq2Lmv+rPbvL0H2VZ6kPrGZgDMoH9yPCf26sJFI/swKj2JkemJ9Ooa6/E7EJFI0OZFb2ZRwB+Bc4FCYJWZve6c29zW+wq0xqZmSqrqKKqspbiylqLKQ77bWnaWV7OzrIaqusbPnt/BIDUpjkE9O3PG4BRO7NWFIb27cEKPzpqzXUQ8E4gj+vHANudcPoCZLQIuATwpeucctQ3NVNc3UlPX1HJb30h1XROVhxrYX1NPeXU9+6vrKa9pYH91Pftr6ik7WEdpVd1nwyz/FRvTgdTEOPp1jyerfzL9u8eT0T2B/t3jSe8WrytORSToBKLo04CCwx4XAqcEYD+8vKqAPy3/lKZmR2Oz+/xtUzONzY5DDU2ffdD5VRLjYkhO6Ei3+Bj6JMYyPLUrvRPj6JMYS+/EWPokxtKnaxxd46J1OqOIhBTPPow1s1nALIB+/fq16nskxccwtHdXojoY0R2M6CgjqkMHojvYZ9viO0YR3ymahI5RxHeMJqFTy218xygS42LoltCRpLgYokNkeoAPPvjA6wgiEmICUfS7gb6HPU73bfsc59wcYA5AVlZWq84Z/Mbw3nxjeO/WvFREJGIE4jB2FTDYzAaYWUfgauD1AOxHRESOQZsf0TvnGs3s+8A/aDm9cq5z7pO23o+IiBybgIzRO+eWAksD8b1FROT4hMYnkCIi0moqehGRMKeiFxEJcyp6EZEwp6IXEQlz5o5lfoBAhzArBXa28uUpQFkbxgkFes+RQe85Mvjznvs753oc7UlBUfT+MLNs51yW1znak95zZNB7jgzt8Z41dCMiEuZU9CIiYS4cin6O1wE8oPccGfSeI0PA33PIj9GLiMhXC4cjehER+QohXfRmdr6ZbTGzbWZ2t9d5As3M+prZ+2a22cw+MbPbvc7UHswsyszWmtnfvc7SHswsycwWm1mumeWY2USvMwWamf3Q9zO9ycwWmlms15nampnNNbMSM9t02LZkM3vHzPJ8t90Cse+QLfrDFiG/ABgGTDWzYd6mCrhG4E7n3DBgAvC9CHjPALcDOV6HaEePAm8554YCownz925macBtQJZzbgQt05tf7W2qgPgzcP4Xtt0NLHPODQaW+R63uZAteg5bhNw5Vw/8dxHysOWcK3LOrfHdr6KlANK8TRVYZpYOXAQ843WW9mBmicDXgGcBnHP1zrkKb1O1i2ggzsyigXhgj8d52pxzbjlQ/oXNlwDzfffnA1MCse9QLvojLUIe1qV3ODPLAMYAK7xNEnCPAHcBzV4HaScDgFJgnm+46hkzS/A6VCA553YDvwN2AUVApXPubW9TtZtezrki3/1ioFcgdhLKRR+xzKwz8BfgB865A17nCRQz+yZQ4pxb7XWWdhQNjAWedM6NAaoJ0K/zwcI3Ln0JLf/IpQIJZjbN21Ttz7WcAhmQ0yBDueiPaRHycGNmMbSU/AvOuVe9zhNgpwEXm9kOWobmzjKz572NFHCFQKFz7r+/qS2mpfjD2TnAdudcqXOuAXgVONXjTO1lr5n1AfDdlgRiJ6Fc9BG3CLmZGS1jtznOuYe9zhNozrl7nHPpzrkMWv5+33POhfWRnnOuGCgwsyG+TWcDmz2M1B52ARPMLN73M342Yf4B9GFeB2b47s8AlgRiJwFZM7Y9ROgi5KcB1wIbzWydb9tPfWv0Svi4FXjBdwCTD1zvcZ6Acs6tMLPFwBpazixbSxheIWtmC4FJQIqZFQL3AbOBl81sJi0z+F4VkH3rylgRkfAWykM3IiJyDFT0IiJhTkUvIhLeAMZOAAAAI0lEQVTmVPQiImFORS8iEuZU9CIiYU5FLyIS5lT0IiJh7v8DhCgcPgVvtqQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0, 10, 1000)\n",
    "y = x**2\n",
    "plt.plot(x,y)\n",
    "plt.vlines(2, 0, 100)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's hard to do that using some simple classical optimization methods which work only for unconstrained problems. However, we can artificially add a penalty, for values lower than 2 or higher than 8:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "penalty = np.piecewise(x, [x<2, x>=2], [lambda x: 100*(-x+2)**3, 0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4lNXB/vHvyb4nhIRACDvIJrJFgkvVim1xxddaRaVFRdC6azftovbt8mJb29rNn4haRECU0oK1irtglTXsOwRCAoEEQkLInsz5/TEDL69FgdmeWe7PdXFNMszkuQfCncOZ85zHWGsREZHIFeN0ABERCSwVvYhIhFPRi4hEOBW9iEiEU9GLiEQ4Fb2ISIQ7ZdEbY14wxlQaYzaccF+2MeYdY8x2z20Hz/3GGPMHY8wOY8w6Y8yIQIYXEZFTO50R/V+BsZ+57xHgPWttP+A9z+cAlwP9PL+mAM/4J6aIiHjrlEVvrV0MVH/m7nHADM/HM4BrT7j/Jeu2FMgyxnTxV1gRETlzcV4+L89aW+H5eD+Q5/m4K1B2wuPKPfdV8BnGmCm4R/2kpqaOHDBggJdRHLJvDaTnQXpwf45t3boVgP79+wf1uCLiXzUNLZQdbqRLZhI5aYlefY1Vq1YdtNbmnupx3hb9cdZaa4w5430UrLXTgGkAhYWFduXKlb5GCa5fFsCIb8HYXwb1sJdccgkAH374YVCPKyL+U1bdwOVPL2FclwzmTBlNbIzx6usYY0pP53Herro5cGxKxnNb6bl/L9DthMcVeO6LPAmp0HLU6RQiEmba2l08NHcNBvjtjUO9Lvkz4W3RLwQmej6eCCw44f5veVbfjAZqT5jiiSwJqdBS73QKEQkzz3y4k5Wlh/nZtWdT0CElKMc85dSNMWYOcAmQY4wpBx4HpgKvGmMmAaXADZ6H/wu4AtgBNAC3BSBzaNCIXkTO0JqyGn7/3nauGZrPtcO7Bu24pyx6a+1Nn/NbY07yWAvc42uosJCYrhG9iJy2+uY2HnxlNZ0zkvjZtWcH9dg6M9ZbGtGLyBn42T83UVrdwFM3DCUzOT6ox1bRe0tz9CJymt7asJ9XVpRx18V9GN27Y9CPr6L3VkIqNGtELyJf7MCRJh6dv44hXTN56LKzHMmgovdWguboReSLtbssD81dQ1Ori9+PH0ZCnDOVq6L31rE5el1zV0Q+x7OLd/LJzkP89JrB9MlNcyyHit5bCalg26Gt2ekkIhKCVu85zG/f3saV53ThG4UFjmZR0XsrwfPTWdM3IvIZR5pauf+V1eRlJPHL/xqCMYE/+/WLqOi9lZDqvm2pczaHiIQUay0/+ccG9tU08YebhgV9KeXJqOi9lagRvYj8p/nFe1mwZh8PjunHyB7ZTscBVPTeOz6iV9GLiNuug/X8ZMEGinplc/eX+zod5zgVvbeOz9FrLb2IQEubi/vnrCYhLobfjx8WlF0pT5fP+9FHrWMjep00JSLAb97eyvq9tTz7zZF0yUx2Os7/oRG9t7TqRkQ8Fm+rYtriEiaM7s7XBnd2Os5/UNF7S1M3IgIcPNrMw6+u5ay8NH585SCn45yUpm68pTdjRaKey2X57mtrqWtqZdYdRSTFxzod6aQ0ovdWfDKYGBW9SBSbtqSED7dW8eMrB9K/c7rTcT6Xit5bxrinbzR1IxKVVpVW8+tFW7liSGcmjO7hdJwvpKL3hS4+IhKVDte3cN/s1XTNSmbq189xfIuDU9EcvS908RGRqGOte17+4NEW/vbt88lIcn6Lg1PRiN4XKnqRqDN9yS7e21LJD68YwJCCTKfjnBYVvS8S0qFZm5qJRIviPYd58q0tjB3cmYnn93Q6zmlT0fsiUUUvEi1qGtzz8l2yknjy+tCflz+R5uh9oaIXiQruefl1VNY1Me+u80Ni6+EzoRG9L1T0IlHhhX/v5t3NB3j08oEM7ZbldJwzpqL3RWKail4kwq0pq2Hqm5v56qA8brugp9NxvKKi90ViOrQ3Q1uL00lEJABqG1u5d3YxndKT+PX1Q8NqXv5EKnpfJGa4b3XSlEjEsdby/Xlr2V/bxJ9uHk5mSnjNy59IRe+LRM/eFs1HnM0hIn43fckuFm08wCOXD2B49w5Ox/GJit4Xx4te8/QikWRZySGmvrWFK4Z0ZtKFvZyO4zMVvS9U9CIRp7KuiXvnrKZHdgpPhsE+NqdD6+h9oaIXiSht7S7unb2ao01tvDypiPQw2MfmdKjofZGgoheJJL9+eyvLd1XzuxuHhvT+8mdKUze+0IheJGIs2rifZz9yX/f1v4YXOB3Hr1T0vlDRi0SE3Qfr+e6raxlakMlPrgrN6776QkXvi4RUwKjoRcJYY0s7d728ithYw59vGUFiXGhe99UXPhW9MeYhY8xGY8wGY8wcY0ySMaaXMWaZMWaHMWauMSbBX2FDjjHuk6ZU9CJhyVrLTxZsYOuBOn5/4zAKOqQ4HSkgvC56Y0xX4H6g0Fp7NhALjAeeBH5nre0LHAYm+SNoyNLGZiJha+6KMuatKue+S/txSf9OTscJGF+nbuKAZGNMHJACVACXAvM8vz8DuNbHY4S2xDSdGSsShjbsreWxhRv5Ur8cHhjTz+k4AeV10Vtr9wK/AfbgLvhaYBVQY61t8zysHOh6sucbY6YYY1YaY1ZWVVV5G8N5iena60YkzFTXt3DnzFXkpCbw9PjhxMaE/0lRX8SXqZsOwDigF5APpAJjT/f51tpp1tpCa21hbm6utzGcp6kbkbDS1u7ivjnFVB1t5pkJI8lOjdy3EY/xZermMmCXtbbKWtsKzAcuALI8UzkABcBeHzOGNhW9SFj59aKt/HvHIX5+7dlheRERb/hS9HuA0caYFOPeDGIMsAn4ALje85iJwALfIoY4Fb1I2Hh97T6eXVzCN0f34IbCbk7HCRpf5uiX4X7TtRhY7/la04AfAA8bY3YAHYHn/ZAzdGl5pUhY2FxxhO/PW0dhjw4ReVLUF/Fprxtr7ePA45+5uwQY5cvXDSsJnssJulwQo/PPREJRTYP7zdeM5Dj+MmEECXHR9W81ul5tICSmAxZa651OIiIn0e6y3P/KGipqG3lmwkg6pSc5HSnoVPS+Or7fjZZYioSip97eyuJtVfz3uLMZEeZXivKWit5X2thMJGS9ub6Cv3y4k5tGdeemUd2djuMYFb2vjl0gXEUvElK2HajjO6+tZXj3LJ64JrrefP0sFb2vdIFwkZBT29jKnTNXkZIQxzO3jIzIHSnPhIreV4lp7luN6EVCgvvM19WUH27gmQkj6JwZfW++fpYuJegrzdGLhJSpb25h8bYqpl43hHN7ZjsdJyRoRO+rY3P02thMxHGvrSxj+se7uPX8noyP4jdfP0tF76uEY1M3mqMXcdKq0mp+9PcNXNg3hx9fOdDpOCFFRe+ruASIS4ImFb2IU/bVNHLnzGLys5L4083DiYtVtZ1Ic/T+oP1uRBzT2NLO5JdW0tzazitTishKifxth8+Uit4fkjKhqdbpFCJRx1rLd+etZVPFEV6YeC59O6U7HSkk6f83/pCUoaIXccAf39/BG+sqeGTsAL48IHKv+eorFb0/aEQvEnRvbdjPb9/ZxnXDuzLlot5OxwlpKnp/UNGLBNXmiiM8/OoahnXL4pfXDcF97SP5PCp6f0jK1PJKkSCprGvijhkryUiKZ9o3R5IUH93bG5wOvRnrDxrRiwRFU2s7k19aRXV9C6/ddR6dMrS9wenQiN4fEjOgrQlam5xOIhKxXC7Ld15dy7ryGp4eP4yzu2Y6HSlsqOj9IcnzDafpG5GAeeqdrbyxvoIfXj6Qrw7u7HScsKKi94ekLPetpm9EAmLeqnL+/IH7AiJ3fKmX03HCjoreH46N6LUNgojfLS05xKPz13Fh3xz+e9xgrbDxgoreH44XfY2zOUQiTEnVUe6cuYoeHVP58y0jiNceNl7Rn5o/JHm2KtbUjYjf1DS0MGnGSmJjDC9MPJfM5HinI4UtLa/0h+MjehW9iD+0tLm4c+Yq9tY0MmdyEd07pjgdKaxpRO8PKnoRv7HW8uj89SzbVc2vrz+HkT10lShfqej9IT4FYuJU9CJ+8Mf3d/C34nIevKwf44Z1dTpORFDR+4Mx2gZBxA9eW1nm3qhsRFceGNPP6TgRQ0XvL9oGQcQni7dV8ej89VzYN4ep152jZZR+pKL3l0TtSS/irY37arl7VjF9O6XxzIQRJMSpmvxJf5r+ohG9iFf21jRy24srSE+K46+3jSI9Scso/U1F7y8qepEzVtvQyq0vLKextZ2/3jaKzpnajTIQVPT+oqIXOSPNbe1MmbmS3YfqefabI+nfWdd7DRSdMOUvSZna60bkNLlclu++to5lu6p5evwwzu+T43SkiKYRvb8kZUJrPbS3Op1EJOQ9uWgLr6/dx/fH9tda+SDwqeiNMVnGmHnGmC3GmM3GmPOMMdnGmHeMMds9tx38FTakaQdLkdPy0qe7efajEiaM7s63L+7jdJyo4OuI/mngLWvtAGAosBl4BHjPWtsPeM/zeeTTDpYip/TGugoeX7iRywZ24omrteVwsHhd9MaYTOAi4HkAa22LtbYGGAfM8DxsBnCtryHDgva7EflCn+w4yENz1zCyewf+eNMI4rTlcND48ifdC6gCXjTGrDbGTDfGpAJ51toKz2P2A3kne7IxZooxZqUxZmVVVZUPMUKELico8rnWl9cy+aWV9MpJ5fmJ55KcEOt0pKjiS9HHASOAZ6y1w4F6PjNNY621gD3Zk62106y1hdbawtzcXB9ihIhE7UkvcjK7DtZz64vLyUpJYMbto8hM0QlRweZL0ZcD5dbaZZ7P5+Eu/gPGmC4AnttK3yKGiWMj+kbN0YscU3mkiW8+vwwLvDRJJ0Q5xeuit9buB8qMMf09d40BNgELgYme+yYCC3xKGC6SPYuLGg87m0MkRNQ2tvKtF5ZTXd/Ci7eeS5/cNKcjRS1fT5i6D5hljEkASoDbcP/weNUYMwkoBW7w8RjhISEVYhNU9CJAU2s7k2esZGfVUV649VyGdstyOlJU86norbVrgMKT/NYYX75uWDLGPapX0UuUa2t3cd+c1aworeYP44fzpX4R8B5cmNP6Jn9S0UuUs9byo79v4J1NB3ji6sFcPTTf6UiCit6/VPQSxay1/OKNzcxdWcb9l/Zl4vk9nY4kHip6f0rOVtFL1PrDezuY/vEubj2/Jw995Syn48gJVPT+pBG9RKnpS0r43bvbuH5kAY9dNUhbG4QYFb0/JWep6CXqzFm+h5+/sZkrhnRm6nVDiIlRyYcaFb0/JXeA1gZobXI6iUhQLFy7jx/+fT0Xn5XL728crv1rQpT+VvwpJdt9q1G9RIF3Nx3g4blrOLdnNv9vwkhd0DuE6W/Gn3R2rESJT3Yc5O7ZxQzKz+D5iYXapCzEqej9SUUvUWBV6WHueGklPTumMOO2UaQnaZOyUKei96fjRV/tbA6RANmwt5bbXlxObnoiL08qokNqgtOR5DSo6P0pWXP0Erk2VxxhwvPLSEuM4+VJRXTK0E6U4UJF70+aupEIte1AHbdMX0ZSXCxzpoymW3aK05HkDKjo/SkhFWLiVfQSUXZUHuXm55YRF2OYPbmIHh1TnY4kZ0hF70/HdrBs0By9RIZdB+u5+bmlAMyePJre2lM+LKno/S1F+91IZCg9VM9N05bS7rLMnlxE304q+XClovc37XcjEaCsuoGbn1tGc1s7syYXcVZeutORxAcqen9L7qDrxkpY21vTyE3PLeVocxsv31HEgM4ZTkcSH6no/S25g9bRS9iqqG3k5ueWUtvYysuTihicn+l0JPEDFb2/aepGwtTemkbGT1tK9dEWZk4qYkiBSj5SqOj9TTtYShgqq27gxmc/pbq+hZl3FDFMF/OOKCp6f9NJUxJmSg/Vc+Ozn1LX1MbsO0ar5COQit7fVPQSRkqqjnLjs0tpbG1n9mRN10SqOKcDRJxje9I3HHI2h8gpuM94da+TnzNltFbXRDAVvb+l5LhvGw46m0PkC2w7UOc549XwypTR9NM6+YimqRt/S/UUfb2KXkLT5oojjJ+2lBijko8WGtH7W0pH962mbiQEbdhby4Tnl5EcH8vsyaPplaMNyqKBRvT+FhsPSZka0UvIWVVazU3PLSU1IY65U85TyUcRFX0gpORojl5CysfbDzJh+nJy0hJ59a7z6N5R+8lHE03dBEJqjkb0EjIWbdzPfbNX0zs3lZmTishNT3Q6kgSZRvSBkJKjOXoJCX9fXc7ds4oZlJ/B3CnnqeSjlIo+EFI7akQvjpv56W4emruWol7ZvHxHEZkp8U5HEodo6iYQjo3oXS6I0c9SCb6/fLiDX721lcsGduJPN48gKT7W6UjiIBV9IKTmgG2Hppr/PVNWJAistfxq0Vae+XAn1wzN56kbhhIfq8FGtFPRB8Lxs2MPqeglaNpdlscWbGDWsj3cXNSdn407m9gY43QsCQH6UR8IqZ6TpjRPL0HS1NrOPbOKmbVsD3de3JtfXKuSl//lc9EbY2KNMauNMf/0fN7LGLPMGLPDGDPXGJPge8wwo/1uJIiONLUy8YXlvLVxPz+5ahCPXj4QY1Ty8r/8MaJ/ANh8wudPAr+z1vYFDgOT/HCM8KL9biRIKo80ceOzSynec5inxw9j0oW9nI4kIcinojfGFABXAtM9nxvgUmCe5yEzgGt9OUZY0ohegmDXwXque+YTSg/V8/zEcxk3rKvTkSRE+fpm7O+B7wPHtr/rCNRYa9s8n5cDJ/3uM8ZMAaYAdO/e3ccYISY+CRLSoF4nTUlgrCuv4bYXV2CBOZNHM1RXhZIv4PWI3hhzFVBprV3lzfOttdOstYXW2sLc3FxvY4SulI4a0UtALNlexfhpS0lOiGXeXeep5OWUfBnRXwBcY4y5AkgCMoCngSxjTJxnVF8A7PU9ZhhKzYX6KqdTSIT5x+q9fG/eWvp2SmfGbefSKSPJ6UgSBrwe0VtrH7XWFlhrewLjgfettbcAHwDXex42EVjgc8pwlNYJjqroxT+stfzp/e08OHcNhT2ymXvnaJW8nLZArKP/AfCwMWYH7jn75wNwjNCX1gmOHnA6hUSA1nYXP/jbOn7z9jauG96VGbePIiNJ+9bI6fPLmbHW2g+BDz0flwCj/PF1w1panvvM2PZW98VIRLxwpKmVe2YVs2T7QR4Y048HL+unNfJyxrQFQqCk5QHWvZY+o4vTaSQM7atp5Pa/rmBH5VF+ff05fKOwm9ORJEyp6AMlLc99e/SAil7O2Ia9tdz+1xU0trQz4/ZRXNA3x+lIEsZU9IFyvOgrnc0hYeeDLZXcM7uYDikJzPx2Ef07p5/6SSJfQEUfKGmd3Ld6Q1bOwMxPd/P4wo0Mys/ghYlaPin+oaIPFBW9nIHWdhc/fX0jLy/dw2UDO/H0+OGkJuqfp/iHvpMCJT4ZEjM1dSOnVNPQwt2zivlk5yHuurgP3/taf20xLH6log8kraWXU9hReZQ7ZqxgX00TT31jKF8fWeB0JIlAKvpASsvTiF4+10fbqrh3djGJcTHMmVLEyB66GpkEhq4wFUga0ctJWGt54eNd3Pbicgo6pPCPey5QyUtAaUQfSBrRy2e0tLl4fOEG5iwv46uD8vjdjcP0pqsEnL7DAik9D1rqoKUeElKdTiMOq6xr4p5ZxazYfZh7vtyH73ylPzF601WCQEUfSCeeNJWtS7xFs+I9h/n2y6s40tjGH24azjVD852OJFFEc/SBpLX0Asxetocbn/2UxLhY5t99vkpegk4j+kBK9+xxU1fhbA5xRHNbO48v2MgrK8q4+Kxcnh4/jKyUBKdjSRRS0QdShmfkdmSfszkk6CpqG/n2y8WsKavhni/34eGv6CQocY6KPpCSsiA+RUUfZZbvqubuWatobGnn/00YwdiztXupOEtFH0jGuEf1R6LzsrnRxlrLc0tKePKtrfTITmHO5NH0y9POk+I8FX2gZeRrRB8Fahta+c5ra3l38wEuP7szT15/ji73JyFDRR9oGV1h98dOp5AAWldew92zijlwpInHrx7Eref31OX+JKSo6AMtI9+96sbVDjGxTqcRP7LWMnNpKT//52Zy0xN59c7zGN69g9OxRP6Dij7QMvLB1Qb1VZDe2ek04id1Ta08Mn89b6yr4NIBnXjqG0PpkKqlkxKaVPSBltHVfXtkr4o+Qmzad4R7Zhezp7qBH4wdwJ0X9dZWBhLSVPSBduJa+q4jnc0iPrHWMuOT3fzyzS1kJccz+44iinp3dDqWyCmp6APt+IheK2/C2aGjzXx/3jre21LJpQM68evrz6FjWqLTsUROi4o+0FI6QmyC1tKHsY+3H+ThV9dQ09jKE1cPYqJW1UiYUdEH2vGTpjSiDzctbS6eemcr0xaX0Cc3jRm3j2JglwynY4mcMRV9MGR0hVqN6MPJ7oP13P/KataV13LTqO48dtUgkhO0PFbCk4o+GDK7Qem/nU4hp8Fay2sry/np6xuJi43RXjUSEVT0wZDVHda/Cu2tEKvT4kNVVV0zj85fx7ubKxndO5vf3jCM/Kxkp2OJ+ExFHwwdeoB1QW0ZZPd2Oo2cxJvrK/jRPzZQ39zGY1e5tzHQ2niJFCr6YMjq4b49XKqiDzG1ja08sXAjf1+9lyFdM/ndjUPp20k7TkpkUdEHQwdP0deUOptD/o+Ptx/ke/PWUlnXzANj+nHvpX2Jj9XVNSXyqOiDIT0fTCzU7HE6iQD1zW08+dYWXvq0lD65qcz/9vkM7ZbldCyRgFHRB0NsHGQWuKduxFGLt1Xx6Pz17Ktt5LYLevKDsQNIiteySYlsXhe9MaYb8BKQB1hgmrX2aWNMNjAX6AnsBm6w1h72PWqY69BDUzcOqm1o5WdvbGLeqnJ656by2p3nUdgz2+lYIkHhy4i+DfiOtbbYGJMOrDLGvAPcCrxnrZ1qjHkEeAT4ge9Rw1xWD9i2yOkUUemtDfv5yYINVNe3cPclfbh/TD+N4iWqeF301toKoMLzcZ0xZjPQFRgHXOJ52AzgQ1T07qKvr4SWBkhIcTpNVKiqa+aJhRt5Y30Fg7pk8OKt53J210ynY4kEnV/m6I0xPYHhwDIgz/NDAGA/7qmdkz1nCjAFoHv37v6IEdqOrbypLYPc/s5miXDWWv5WvJefv7GJhuZ2vve1/ky5qLdW1EjU8rnojTFpwN+AB621R07c1c9aa40x9mTPs9ZOA6YBFBYWnvQxEeX4WvrdKvoA2n6gjh/9YwPLd1UzskcHnvz6EK2Ll6jnU9EbY+Jxl/wsa+18z90HjDFdrLUVxpguQKWvISNCxz7u20M7gK85GiUSNba088f3tzNtcQmpiXFMvW4INxR209mtIvi26sYAzwObrbW/PeG3FgITgame2wU+JYwUKR0hKQsObnc6ScR5f8sBHluwkfLDjVw/soBHLx+gi4KInMCXEf0FwDeB9caYNZ77foi74F81xkwCSoEbfIsYIYyBnH6eEb34Q0VtIz9duIm3Nu6nb6c0XpkymtG6tJ/If/Bl1c3HwOf9v3iMt183onXsByUfOJ0i7DW1tjN9SQl//mAnFsv3x/bnjgt7kxCnN1tFTkZnxgZTTl9YOxua6yBRbxCeKWstizYe4Bf/2kRZdSNjB3fmR1cOpFu2lquKfBEVfTB17Oe+PbQD8oc7myXMbDtQx09f38i/dxzirLw0Zt1RxAV9c5yOJRIWVPTBlOMp+oMq+tNV29DK797dxsylpaQmxPLE1YOYMLoHcVoTL3LaVPTBlN0bTAwc0sqbU2lpczFrWSl/eG87tY2t3DSqO9/5an+yUxOcjiYSdlT0wRSX6L6soJZYfi5rLW+sr+BXb21lT3UD5/fpyI+uHMjgfG1dIOItFX2wdeynov8cy0oO8cs3t7C2rIYBndP5623ncvFZuZx4trWInDkVfbB1Ggi7PoL2Nvc+9cKOyjqmvrmVdzcfoHNGEr+6/hy+PqKAWJ3VKuIXappgyzsb2lvc8/SdBjqdxlFl1Q388f3t/K14L8nxsXzva/25/YJeJCdoC2ERf1LRB1veYPftgY1RW/T7a5v40wfbmbuiDGMM3zqvB/d+ua+2LRAJEBV9sOWcBTHxcGADDLne6TRBVVXXzDMf7uTlZaW4XJYbz+3GvZf2pUtmstPRRCKaij7Y4hLc2xQf2Oh0kqA5XN/Cs4tLmPHJblraXVw3vCv3j+mnM1pFgkRF74S8wbD7Y6dTBFzlkSaeW1LCrGV7aGxt55qh+Twwph+9c9OcjiYSVVT0TsgbDOvmQv0hSI283RbLqht4dvFOXl1ZTlu7i2uG5nP3l/tyVp729xFxgoreCfkj3Lf7iqHfV5zN4kc7Kuv4y4c7WbBmHzEGrh9ZwF0X96FHx1Sno4lENRW9E/KHu7dCKF8ZEUW/qvQw05eU8NbG/STGxTDxvJ5MvqiX3mQVCREqeickpkHuANi7yukkXmtrd7Fo4wGmf1zC6j01ZCTFcfclfbj9gl5aJikSYlT0Tuk6Arb8C6x1X30qTNQ1tTJ3RRkv/ns3e2sa6dExhZ9eM5jrRxaQmqhvJ5FQpH+ZTulaCKtfhsO7IbuX02lOqfRQPTM/LeWVFWUcbW5jVM9sHrt6EJcNzNNWBSIhTkXvlIJz3bd7loZs0be7LO9vqeTlpaV8tK2K2BjDVed0YdKFvTinIMvpeCJymlT0Tuk0CJKzYfcSGHaT02n+j8q6JuYuL2PO8j3sq20iLyORBy/rx/hzu9M5M8npeCJyhlT0TomJgZ4XwK4lTicB3KP3T3Ye5JUVZSzasJ82l+XCvjk8dvUgxgzMI15XdBIJWyp6J/W8CDa/7p6n79DTkQi7D9Yzb1U584vL2VfbRGZyPBPP78ktRd11BqtIhFDRO6nXl9y3uxYHtejrm9t4Y30F81aWs3x3NTEGvtQvlx9eOZDLBuaRFK9tgkUiiYreSbkDID0fti2CEd8K6KFa2lws2V7F62v38famAzS0tNMrJ5Xvfa0/Xx9RoLl3kQimoneSMdD/clg7B1obId6/Z5K2uyzLSg6xcO0+3tywn9rGVjKT4xk3LJ+vjyhgZI8OukyfSBRQ0TsiiXQ0AAAHPElEQVRtwBWw8nko+Qj6j/X5y7W7LCt3V/Pmhv28sb6CqrpmUhNi+ergzlw9tAsX9s0lIU5vrIpEExW903peBIkZsOWfXhd9U2s7S7Yf5O2N+3lvSyXV9S0kxMVwaf9OXD00n0sHdNLl+USimIreaXEJ0P8K2LQALv8VJJzexTgOHW3mw61VvL1pP4u3HaSxtZ30pDjGDOjEVwd35qKzcknTlgQigoo+NAyfAOtecS+1HHrjSR/S2u5i9Z4aDhdcSGNWLwp/8S7WQueMJL5RWMBXB3VmVK9sTcuIyH9Q0YeCnhdCh16weubxorfWUnqogU92HuKjbZV8suMQdc1t0LWIxLp9PHzZWVzcP5chXTP1hqqIfCEVfSgwBjvyVsy7j/Pm2//izUNdWLbrEAeONAPQNSuZq4Z24eKzcvnFA7cS297MfX+50+HQIhIuVPQOaWptZ+O+I6wpq6G49DDrS/qw0KYSt+Q3LE36IUW9O1LUK5vRvTvSJzf1+Kh9anuzw8lFJNyo6IPA5bKUHKxnbVkNazy/Nlccoc1lAcjPTKLorG6UtU3kK9v/wmW3JGN6Dnc2tIhEDBW9n9U2trJ1fx2bK46wZf8RNlfUsXV/HY2t7QCkJsRyTkEWky/qzbBuWQzvlkWnDM9Zqc194Zl/Yl6/H+76N8TrbFUR8Z2K3gvtLsu+mkZ2Haw//qvkYD07K4+yt6bx+OMyk+MZ2CWd8aO6MbBLBsO6ZdEnN+3zL9SRmAZXPw0z/wveeBjG/Tmsrj4lIqEpIEVvjBkLPA3EAtOttVMDcZxAaWlzsb+2iX21jeyraaSitom9NY1U1DRSfriR0kMNtLS7jj8+NSGWXrmpjOzRgVtGd2dg5wwGdskgLyPxzFfE9LkULn4EPpoKMXHutfUa2YuID/xe9MaYWODPwFeAcmCFMWahtXaTv4/1WS6Xpc1laXdZ2lwuWtpcNLS0c7S5jYaWNo42t1Pf3EZ9cxt1TW1U17dwqL6F6vrmEz5uoaah9T++dnZqAl0yk+iZk8qlAzrRKyfV/Ss3ldw0Lwr9i1zyCLhaYclT7guTnHcP9BkDWd39dwwRiRrGWuvfL2jMecAT1tqveT5/FMBa+z+f95zCwkK7cuXKMz7W8x/v4jeLth4vdtcZvpQY4y7wY786pia6b9MSyM9KJj8zmfysJLpkJjuzhcDO92HRj6Fy4/G7WlyGNpchJeX0zqAVkRA39n+83r3WGLPKWlt4qscFYuqmK1B2wuflQNFnH2SMmQJM8Xx61Biz1cvj5QAHvXwuu7x9orNyoMbr1xymfPp7DlN6zdHgRxNzYKK3r7nH6TzIsTdjrbXTgGm+fh1jzMrT+YkWSfSao4Nec3QIxmsOxMYoe4FuJ3xe4LlPREQcEIiiXwH0M8b0MsYkAOOBhQE4joiInAa/T91Ya9uMMfcCi3Avr3zBWrvxFE/zhc/TP2FIrzk66DVHh4C/Zr+vuhERkdCizctFRCKcil5EJMKFddEbY8YaY7YaY3YYYx5xOk+gGWO6GWM+MMZsMsZsNMY84HSmYDDGxBpjVhtj/ul0lmAwxmQZY+YZY7YYYzZ7TkKMaMaYhzzf0xuMMXOMMRG374cx5gVjTKUxZsMJ92UbY94xxmz33HYIxLHDtuhP2GrhcmAQcJMxZpCzqQKuDfiOtXYQMBq4JwpeM8ADwGanQwTR08Bb1toBwFAi/LUbY7oC9wOF1tqzcS/iGO9sqoD4KzD2M/c9Arxnre0HvOf53O/CtuiBUcAOa22JtbYFeAUY53CmgLLWVlhriz0f1+EugK7OpgosY0wBcCUw3ekswWCMyQQuAp4HsNa2WGtrnE0VFHFAsjEmDkgB9jmcx++stYuB6s/cPQ6Y4fl4BnBtII4dzkV/sq0WIrr0TmSM6QkMB5Y5myTgfg98H3Cd6oERohdQBbzoma6aboxJdTpUIFlr9wK/AfYAFUCttfZtZ1MFTZ61tsLz8X4gLxAHCeeij1rGmDTgb8CD1tojTucJFGPMVUCltXaV01mCKA4YATxjrR0O1BOg/86HCs+89DjcP+TygVRjzARnUwWfda91D8h693Au+qjcasEYE4+75GdZa+c7nSfALgCuMcbsxj01d6kx5mVnIwVcOVBurT32P7V5uIs/kl0G7LLWVllrW4H5wPkOZwqWA8aYLgCe28pAHCSciz7qtlow7k3vnwc2W2t/63SeQLPWPmqtLbDW9sT99/u+tTaiR3rW2v1AmTGmv+euMUDAr+XgsD3AaGNMiud7fAwR/gb0CRYCEz0fTwQWBOIgYXspQQe2WggFFwDfBNYbY9Z47vuhtfZfDmYS/7sPmOUZwJQAtzmcJ6CstcuMMfOAYtwry1YTgVshGGPmAJcAOcaYcuBxYCrwqjFmElAK3BCQY2sLBBGRyBbOUzciInIaVPQiIhFORS8iEuFU9CIiEU5FLyIS4VT0IiIRTkUvIhLh/j8r/ly0MkgIugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y)\n",
    "plt.plot(x,penalty)\n",
    "plt.vlines(2, 0, 100)\n",
    "plt.ylim([0,100])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So now we solve for a slightly modified function, which has additional cost for values lower than 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VdW5x/Hvm5zMcyCEWZDRQAUxMji0KnrrTGtbFUccwGu12tFiW2vb2956a2vHWytOoCKiYKvVOlKpQxUJ86xhDgSSQAKZx3X/yMGmXlTIGfYZfp/n8UnO5pzsNxJ+WWfttd9lzjlERCR2JXhdgIiIhJaCXkQkxinoRURinIJeRCTGKehFRGKcgl5EJMZ9atCb2cNmVmFma7scyzezV83sA//HPP9xM7PfmVmpma02s3GhLF5ERD7dkYzoZwPnfOTYTGCRc24YsMj/GOBcYJj/vxnAfcEpU0REuutTg9459waw/yOHpwBz/J/PAb7Q5fijrtO7QK6Z9QlWsSIicvR83XxdoXOu3P/5HqDQ/3k/YGeX55X5j5XzEWY2g85RPxkZGSeOHDmym6V4Y2N5LZmpPvrnpYX1vJs2bQJgxIgRYT2viARXTUMLO6sb6ZOTSs/MlG59jWXLllU55wo+7XndDfoPOeecmR11HwXn3CxgFkBxcbErKSkJtJSwOvNXiynqk80fLg/vZYjTTz8dgMWLF4f1vCISPDv3N3Dub99kSp9s5s2YSGKCdevrmNn2I3led1fd7D00JeP/WOE/vgsY0OV5/f3HYk5Gso/65javyxCRKNPW3sHX56/EgHsvHdPtkD8a3Q3654Br/J9fAzzb5fjV/tU3E4EDXaZ4YkpGSiL1Le1elyEiUeaPizezbHs1P/3iaPrnpYflnJ86dWNm84DTgZ5mVgbcBdwNPGVm1wPbgUv8T/8bcB5QCjQA14ag5oiQkexjb22T12WISBRZvqOa3y76gClj+zJlbL+wnfdTg945N/Vj/mjyYZ7rgJsDLSoapKf4qK/SiF5EjkxdcxvfmL+S3tmp/GTK6LCeO+CLsfEqMyVRc/QicsR+8td17NzfwJMzJpGTlhTWc6sFQjelJ/to0By9iByBF9eU81RJGTedPoTxg/PDfn4FfTdlJCdS39KGdugSkU+y50ATM59Zw/H9c/j6WcM9qUFB300ZKT6cg8ZWjepF5PA6OhzffnoVLW0d/ObSsSQlehO5CvpuSk/pvLxR36ygF5HDe+DNLbxVWsWdFxRxbEGmZ3Uo6LspIzkRQBdkReSwVu2s4Z6XN3Hu6N5MHT/g018QQgr6bso4NKJvUdCLyL+rbWrl1idXUJidyt0XH49Z6O9+/SRaXtlNGcmd/+u08kZEunLOcedf1rJzfwNP3TiJnPTwLqU8HI3ouyk9pXPqpk5TNyLSxTPLd/GXlbv5+lnDKR4U/qWUh6Og76ZM/9RNgy7Giojflso67nx2LeMH53PzGUO9LudDCvpuSj90MVZz9CICtLR1cOuTK0j2JfDby8aGpSvlkdIcfTcdmqPXqhsRAbjn5Y2s3XWQ+686kT454d2Q6NNoRN9Nh1bd6GKsiCzeVMEDb27lqonH8PlRvb0u5/9R0HdTsi+BpETTiF4kzlXUNvHtp1cxsncW3z//OK/LOSxN3QQgXbtMicS1jg7Ht55aRV1zG/OmTyQ1KdHrkg5LI/oAZKb4tMuUSBx74M0tvPlBFT+8YBTDCrO8LudjKegDkJ6cSINW3YjEpeU7qiOmxcGnUdAHID3FR53W0YvEnZqGFr72xAp656Ry95e8b3HwaTRHH4CsFM3Ri8Qb5zpbD1fUNrHgP08O+25R3aERfQAyU3zUNSnoReLJQ29t5bUNFXzvvOMYMyDX63KOiII+AJmpPvW6EYkjy3dUc/eLG/n8qEKmnTzI63KOmII+AJkpPmqbWr0uQ0TCoOu8/C++PCbi5+W7UtAHIMs/ote+sSKxrXNefjUVtU387+XjomJevisFfQAyU3x0aN9YkZjXOS+/N6rm5btS0AcgM7Vz0ZIuyIrErhVROi/flYI+AId60tfqgqxITKppaOGWKJ2X70rr6AOQpRG9SMzqOi8fLevlP45G9AHITOn8i9cSS5HYE+3z8l0p6APw4dSNRvQiMaVk2/6on5fvSkEfgA+nbjSiF4kZlbXN3PzEcvrnpXHPV6J3Xr4rzdEH4NCIvk43TYnEhLb2Dr42bzkHGluZfe14slOjd16+KwV9AA5tJ6gRvUhs+OUr7/Pulv3ce8kYjuuT7XU5QaOpmwAk+xJI8SVoeaVIDHh53R7+9I/NXDFhIBeP6+91OUGloA+QOliKRL9tVfV8+6lVHN8/hx9eWOR1OUGnoA+QOliKRLfGlnb+8/FlJCYaf7xiHCm+yNz3NRABBb2ZfcPM1pnZWjObZ2apZjbYzJaYWamZzTez5GAVG4k0oheJXs45vv+XNWzaW8tvLh1L/7x0r0sKiW4HvZn1A24Fip1zo4FE4DLgf4BfO+eGAtXA9cEoNFJlpvg0Ry8SpZ54bwfPLN/FbZOHcfqIXl6XEzKBTt34gDQz8wHpQDlwJrDA/+dzgC8EeI6IlpWqEb1INFq1s4YfP7eezw0v4NYzh3ldTkh1O+idc7uAXwI76Az4A8AyoMY5dyj5yoB+h3u9mc0wsxIzK6msrOxuGZ7LTPFR36KgF4km1fUtfHXucgqyUvjNpWNJSIj+m6I+SSBTN3nAFGAw0BfIAM450tc752Y554qdc8UFBQXdLcNzmRrRi0SV9g7HrU+uoLK2mfuuHEdeRkxfRgQCm7o5C9jqnKt0zrUCzwCnALn+qRyA/sCuAGuMaJkpSZqjF4ki97y8iTc/qOInU0ZxfP/oblZ2pAIJ+h3ARDNLt85mEJOB9cDrwJf9z7kGeDawEiNbVqqPlrYOmtu0y5RIpHt+9e4Pb4q6bPxAr8sJm0Dm6JfQedF1ObDG/7VmAd8FvmlmpUAP4KEg1BmxDvW7qW9W0ItEso17DvKdp1dz4jF53HXhKK/LCauAet045+4C7vrI4S3A+EC+bjT5V2OzNvLjYK5PJBrVNLQw49FlZKX6uO+KcST74uteUTU1C9ChfWNrm9XBUiQSdV58XUn5gUaenDGJXtmpXpcUdgr6AGWlaDtBkUj2y1c28cb7lfz84s9w4jF5Xpfjifh6/xICmdp8RCRivbC6nPsWb2bq+IFMjaOLrx+loA9QpnrSi0SkTXtq+c6CVYwbmMuPLoq9jpRHQ0EfoA/n6DV1IxIxDjS0MuOxEjJSfNx35Ykx2ZHyaCjoA5SV0rnVmEb0IpGhvcNx2/wV7K5p5E9XjqMwDi++fpQuxgYoNSmBxATTxViRCPGLlzayeFMlP/viaE48Jt/rciKCRvQBMrPOVsXaIFzEcwuXlXH/G1u4auIxXDHhGK/LiRgK+iDITvNpjl7EY8u2V3PHM2uYdGyPmNwOMBAK+iDITk3ioEb0Ip7ZXdPIjY8to09uKn+8YhxJiYq2rjRHHwTZqUkcbNSIXsQLDS1tTH+0hKbWduZNnxAXbYePln7tBUF2mk8jehEPOOf4ztOrWV9+kN9PPYFhhVlelxSRFPRB0DmiV9CLhNvvFpXywppy7jh3JGeMjN09XwOloA+C7LQkDupirEhYvbimnF+/9j4Xj+vH9NOO9bqciKagD4Ls1CTqmttoa+/wuhSRuLBu9wG++dQqThiYy39/8TN07n0kH0dBHwTZaep3IxIulbXNTJ9TQm56EvdfdSKpSfHd3uBIKOiDIDu1sw2CVt6IhFZTazs3PlbC/oYWHri6mF5Zam9wJLS8Mgiy0/xBr5U3IiHT0eH49tOrWL6jhvuuGMfofjlelxQ1NKIPgmx/B0utvBEJnV+/9j7Pry5n5rkjOfczfbwuJ6oo6INAI3qR0FqwrIzf/72US4sHcONntcLmaCnog+DDoNccvUjQvbtlH3c8s5pThvbgp18crRU23aCgD4IPp240ohcJqi2Vddz42DIG5qfzxytOVA+bbtL/tSDISPaRYJqjFwmm6voWrpu9FF+C8ci08eT43znL0dOqmyBISDCyUnV3rEiwNLe1c+Njy9h9oIl50ycwsEe61yVFNY3ogyQ7zacRvUgQOOe4Y+Ea3tu2n19+ZYx2iQoCBX2QqCe9SHD8/u+lPLNiF986ezgXjenrdTkxQUEfJOpJLxK4P68o495XOxuV3XLmUK/LiRkK+iBRT3qRwLxdWsXtC1Yz6dge/PxiNSoLJgV9kKgnvUj3rd99kBsfW8aQgkzuv/pEUnxqVBZMCvog0aobke7ZXdPItbPfIzPFxyPXnvRhk0AJHgV9kGSn+dSTXuQoHWhsZdoj79HQ3M7s606iT06a1yXFJAV9kBwahagnvciRaW5rZ8ajJWytquf+q09kZO9sr0uKWQr6IFG/G5Ej19lyeDVLtnaulT95SE+vS4ppCvogUb8bkSN390sb+euq3cw8dyRTxvbzupyYF1DQm1mumS0ws41mtsHMJplZvpm9amYf+D/mBavYSPavEb2CXuSTPPL2Vma9sYWrJx2jlsNhEuiI/rfAS865kcAYYAMwE1jknBsGLPI/jnkfbieoEb3Ix3ppbTk/eX49ZxcVcteFo7RWPky6HfRmlgN8FngIwDnX4pyrAaYAc/xPmwN8IdAio8GhDcIPaEQvcljvbN7HrU+uZOyAXH532QkkJijkwyWQEf1goBJ4xMxWmNmDZpYBFDrnyv3P2QMUHu7FZjbDzErMrKSysjKAMiJDbnoyADUNCnqRj1q76wDTHy1hYH46D19zEmnJuiEqnAIJeh8wDrjPOXcCUM9Hpmmccw5wh3uxc26Wc67YOVdcUFAQQBmRISM5kaREo0YjepF/s62qnmmPvEd2qo/Hrh9PXkay1yXFnUCCvgwoc84t8T9eQGfw7zWzPgD+jxWBlRgdzIyctGRqGlq8LkUkYlQcbOKqh5fQ3uF49PoJuiHKI90OeufcHmCnmY3wH5oMrAeeA67xH7sGeDagCqNIXnqSpm5E/A40tnL1w++xr66F2deOZ2ivTK9LiluB7jD1NWCumSUDW4Br6fzl8ZSZXQ9sBy4J8BxRIzc9iWqN6EVobGnnhjlL2VxZxyPTxjNmQK7XJcW1gILeObcSKD7MH00O5OtGq9z0ZHbub/C6DBFPtbZ3cMsTyynZXs3vp57AqcN016vXdGdsEOWmaepG4ptzjpkL17BoYwU/mTKaC47XDlGRQEEfRHkZydQ0aupG4pNzjp+/uJGFy8v4+lnDuGriMV6XJH4K+iDKSUuiqbWDptZ2r0sRCbs//L30w9YGt00e5nU50oWCPojydNOUxKmH3trKr159n4tP6MeP1Nog4ijogyg3vbPfjVbeSDyZv3QH//X8es4Z1ZtffPl4EtTaIOIo6IPoUNBrRC/x4rlVu5n5zBo+N7yA3009AV+iIiUS6W8liHLTOqduDuiCrMSB19bv5ZvzV3LSoHz+dOWJJPsUJ5FKfzNB9K+pG43oJba9XVrFV59Yzqi+2Tx0TbGalEU4BX0Q6WKsxINl2/dzw5wSju2ZwZzrxpPl34tBIpeCPohSkxJI9iWosZnErLW7DjDtkaX0zknlsesnfNieWyKbgj6IzEyNzSRmvb+3lqsffo/s1CQev2ECBVkpXpckR0hBH2S5aclaXikx54O9tVz+wLv4Eoy5N0ygX67aDUcTBX2Q5aYnafMRiSmlFXVMfWAJZsa8GRMZ1DPD65LkKCnogyw3PUlz9BIztlbVc/kD7wKOedMnMKRAPeWjkYI+yPLSkzVHLzFh+756ps56l7YOxxPTJzK0V5bXJUk3KeiDLMd/MbZzu1yR6LRjXwNTZ71Lc1s7T0yfwPBChXw0U9AHWV56Mi3tHTSqg6VEqZ37G5j6wLs0tLYz94aJjOyd7XVJEiAFfZDlpunuWIleu2oaufzBd6ltauXx6ydQ1FchHwsU9EF26AaS6npdkJXoUn6gkcsfeJeahlYev2ECo/vleF2SBEmgm4PLR+Rn+INeK28kipRVN3D5A0uorm/h0evHc3x/beYdSxT0QXYo6PdrRC9RYse+zjn52qZWHrthAmMHKORjjYI+yHr4g76qTkEvkW9rVecSyqa2dp6YPlHTNTFKQR9kOWlJJCYY++ubvS5F5BOVVtRy+QNLaOtwzJs+keP66MJrrFLQB1lCgpGXnqypG4lom/bUcsWD7wLGkzMmap18jNOqmxDokZGsqRuJWOt2H+CyWe+QmGDMv1EhHw80og+BHpka0UtkWl1Ww1UPvUdGciJPTFeDsnihEX0I5Gco6CXyLN9RzRUPLiEr1cf8Gycp5OOIgj4EOqdudDFWIsfbpVVc+eAS8jOSmX/jJAbkp3tdkoSRgj4EemSmUNvURktbh9eliPDKuj1c+8hSBuSl8/SNk7RpSBxS0IeA7o6VSPHnFWXcNHc5x/XNZv6NE+mVnep1SeIBBX0I/OumKU3fiHcee2cb35i/igmD85l7gzbyjmdadRMCaoMgXnLO8cfFm7nn5U2cdVwhf7j8BFKTEr0uSzykoA+BHpkpgIJews85x90vbeT+f2zhC2P7cs9XxpCUqDfu8U5BHwLqdyNeaO9w/OAva5n33g6unDiQn1w0moQE87osiQAK+hBQvxsJt+a2dr751CpeWF3OV08fwnc+PwIzhbx0Cvg9nZklmtkKM3ve/3iwmS0xs1Izm29mcXcFSP1uJJwONrUy7eGlvLC6nO+dN5LbzxmpkJd/E4zJu9uADV0e/w/wa+fcUKAauD4I54g66ncj4VBxsIlL73+Xpdv2c+8lY5jx2SFelyQRKKCgN7P+wPnAg/7HBpwJLPA/ZQ7whUDOEa16ZunuWAmtLZV1XHzfP9m+r56Hpp3ExeP6e12SRKhAR/S/AW4HDt0C2gOocc61+R+XAf0O90Izm2FmJWZWUllZGWAZkacgM4XKWgW9hMbKnTV8+U/v0NjSzrzpE/nc8AKvS5II1u2gN7MLgArn3LLuvN45N8s5V+ycKy4oiL0f0oKszqB3znldisSY1zdVMHXWu2SkJLLgppMZo63/5FMEsurmFOAiMzsPSAWygd8CuWbm84/q+wO7Ai8z+vTKSqW5rYODTW3kpCV5XY7EiIXLyvjuwtUML8xi9nUn0StLLQ3k03V7RO+cu8M51985Nwi4DPi7c+4K4HXgy/6nXQM8G3CVUahXdudNU5q+kWBwzvG/r5fyradXMeHY/M6+NQp5OUKhuGXuu8A3zayUzjn7h0JwjohX4L87tqK2yeNKJNq1tncwc+Ea7nl5ExeN6cvD004iK1XvEuXIBeWGKefcYmCx//MtwPhgfN1ophG9BMPBplZunrucNz+o4mtnDuWbZw/XGnk5arozNkQKMjvfVivopbt21TRy3SNL2VxZxy++dDyXnDTA65IkSinoQyQ7zUeyL0FBL92ypuwA181ZSlNLO7OvHc+pw3p6XZJEMQV9iJgZBZkpVCjo5Si9tn4vX5u3gvyMZOZ+dQLDC7O8LkminII+hHpl66YpOTpz/rmNH/91HaP75fDgNcVaWSNBoaAPoYLMFLbtq/e6DIkCbe0d/PSFDcz+5zbOOq6Q300dS3qy/nlKcOgnKYR6ZaewdNt+r8uQCHegoZWbn1jOW6VV3HDqYO447zgS1UdegkhBH0IFmalUN7TS0tZBsk+7/Mj/V1pRx/RHSyirbuAXXz6eS4q1skaCT0EfQofW0lfVNdM3N83jaiTSLN5UwdfmrSDFl8C86RMpHpTvdUkSozTMDKF/3R2rC7LyL845HnprK9fNXkr/vHT+cvMpCnkJKY3oQ6gwu3PFxN6DaoMgnZrb2rnzL2t5qqSMz48q5N5LxpKRon+GElr6CQuhPrmdQV9e0+hxJRIJquqauenxZSzdVs2tZw7l62cN1+bdEhYK+hDqkZFMcmIC5Qc0oo93K3ZUc9Pjy6luaOH3U0/gwjF9vS5J4oiCPoTMjN45qQr6ODfvvR3c9ew6emWn8MxXT2ZU3xyvS5I4o6APsT45qZQf0NRNPGpua+euZ9fx5NKdfHZ4Ab+7bCy56clelyVxSEEfYn1yUinZXu11GRJmu2sauWnuclbtrOGWM4byjbOH6yYo8YyCPsT65Kaxd005HR1OF97ixDub93HLE8tpbuvgT1eeyDmje3tdksQ5BX2I9c1JpbXdUVXfrAZVMe7Q+vifv7iRQT3Suf+qYob2yvS6LBEFfaj1zum8I7a8pklBH8MONLZy+4JVvLxuL+eM6s0vLxlDptbHS4TQT2KI9cnxr6U/0MiYAbkeVyOhsLqshpufWE55TRM/OP84rj91sLb7k4iioA+xfwW9lljGGuccj76znZ+9sIGCrBSe+s9JjBuY53VZIv+Pgj7E8jOSSfbppqlYc7CplTsWruGFNeVMHtmLX10yRksnJWIp6EPMzPxr6RX0sWLd7gPcPHc5O6sbuePckUw/7VitqJKIpqAPgz45qexWv5uo55zj8SU7+K/n15Ofnsz8GWotLNFBQR8G/fPSefODSq/LkADsr2/h9gWreW3DXk4fUcCvvjKGHv421CKRTkEfBgPz09l7sJmm1nZSkxK9LkeO0tulVXxj/kpqGlr54QVFTDt5kKZqJKoo6MNgQH7nWvqy6kbdQBNFWto6+NWrm5j1xhaGFGQy+9rxFPXN9roskaOmoA+DgfnpAOzc36CgjxJbKuu47cmVrNl1gCsmDOQH5xeRlqx3YxKdFPRhMOBQ0Fc3eFyJfBrnHE+XlHHXc+tISUrg/qtO5POj1KtGopuCPgwKMlNITUpgxz4FfSSrqmvme8+s4ZX1e5l0bA9+felYeueobYVEPwV9GJgZA/LSNaKPYC+tLed7f15LXXMb3z/vOK47dbDaCkvMUNCHyYD8dHbs11r6SHOgoZUf/XUdf16xi9H9srn3krEML8zyuiyRoFLQh8nA/HSWbt2Pc04NryLEP96v5LsLVlNZ18xtk4dxy5lDSUpM8LoskaBT0IfJgPx0apvbqGloJS9DPVG8VN/cxn//bQNzl+xgWK9MHri6mM/01z6uErsU9GEyIK9zLf2O/Q0Keg+9XVrFzGdWU1bdyPTTBvOt/xihm9gk5nX7faqZDTCz181svZmtM7Pb/MfzzexVM/vA/1F9W4FjemQAsG1fvceVxKcDDZ0bg1zx4BJ8CQnMnzGJ759fpJCXuBDIiL4N+JZzbrmZZQHLzOxVYBqwyDl3t5nNBGYC3w281Oh2TI90Egw2Vyrow+2lteXc+ew69te3cNPpQ7ht8jAFvMSVbge9c64cKPd/XmtmG4B+wBTgdP/T5gCLUdCTmpTIgPx0NlfWeV1K3KiobeKuZ9fx4to9FPXJ5pFpJzG6n+biJf4EZY7ezAYBJwBLgEL/LwGAPUDhx7xmBjADYODAgcEoI+INKchkc4WCPtSccyxYVsZPX9hAY2s7t58zgumnHasVNRK3Ag56M8sEFgJfd84d7Lp00DnnzMwd7nXOuVnALIDi4uLDPifWDCnI4O3SKto7nG7GCZHSijru/Mta3tmyj/GD8vn5lz7DkAL1F5L4FlDQm1kSnSE/1zn3jP/wXjPr45wrN7M+QEWgRcaKob0yaW7rYFd1IwN7pHtdTkxpbGnnD69/wKw3tpCWlMjPvjiaqScNVDthEQIIeuscuj8EbHDO3dvlj54DrgHu9n98NqAKY8ihkeXmyjoFfRAt2rCXu55bR1l1I18a1587zhtJT20KIvKhQEb0pwBXAWvMbKX/2PfoDPinzOx6YDtwSWAlxo6uQX/GyF4eVxP9dtU08uPn1vHK+r0M65XJ/BkTmXBsD6/LEok4gay6eQv4uPfFk7v7dWNZXkYy+RnJlOqCbECa29p56K2t/H5RKQAzzx3JdacMJtmni60ih6M7Y8NsaEGmgr6bnHO8un4vP/vbBrbva+A/igr54YVF9M/TNJjIJ1HQh9mwwkyeW7Vbzc2O0qY9tfzX8+t5q7SKYb0yeez68Zw2rMDrskSigoI+zIr6ZjN3yQ7Kqhs/3HlKPl5NQwu/fvV9Hl+yg8wUHz++aBRXTBiIT2viRY6Ygj7Mivp0bi69vvyggv4TtLZ38MSSHfz6tfc52NjKlROP4RtnDVdDOJFuUNCH2cje2SQYrN99UHuRHoZzjhfX7uGelzextaqek4f04IcXFjGyd7bXpYlELQV9mKUlJzK4Zwbryw96XUrEWbJlHz9/cSMrd9YwvDCTh6cVc8aIXrqWIRIgBb0HivrmsHx7tddlRIz399byPy9uZNHGCnpnp/KLLx3Pl07srzYRIkGioPdAUZ9s/rpqN9X1LXE957yrppHfvvY+C5aVkZHs4/ZzRnDtyYNJS1YLYZFgUtB7YOyAXABW7qyJyztk9xxo4o+LS3nyvZ0AXHvKYG45Y2hc/9ITCSUFvQfGDMghMcFYtr06roK+oraJ+xZvZu6SHXR0OL5SPIBbzhxKv9w0r0sTiWkKeg+kJ/s4rk8Wy3fExzz9vrpm7n9jC4++s43WdsfFJ/Tj1snDtLxUJEwU9B4ZNzCPhcvKYro3fWVtMw+9tZVH39lGU2s7U8Z2BvzgnhlelyYSVxT0Hhk3MI9H39nOxj0HGdU3tra327m/gVlvbOGpkp20tndw/vF9uW3yUIb2yvK6NJG4pKD3yEmD8wF4Z/O+mAn6D/bWct/izTy7ajcJBl8a158bPzdEI3gRjynoPdIvN41je2bwVmkVN5x2rNflBGTlzhruW1zKy+v2kpaUyLSTB3HDaYPpk6OLrCKRQEHvoVOH9eTpkjKa29pJ8UXX2vH2Dscr6/bw0FtbKdleTU5aErdOHsa0kweRr2WSIhFFQe+hU4f25NF3trNiRw0To2RnpINNrTy1dCez/7nN34EzjR9eUMQlJw0gM0U/TiKRSP8yPTRxSA98CcbrmyoiPui376vnkbe38XTJTupb2hk/OJ8fnF/E2UWFMbtqSCRWKOg9lJ2axMlDe/Limj3MPGdkxDXvau9wLN5UwdwlO3h9UwWJZlw4pi/XnTKYz/SPjQvIIvFAQe+x8z/Tm+8uXMPaXQcjJjwrDjYxf+lO5r0C86YxAAAGWklEQVS3g90HmuiVlcItZwzlyonHUJid6nV5InKUFPQe+4+i3nzvz2t5YU25p0Hf0eH45+Z9zF2ynVfX76Wtw3Hq0J7ceUERZxUVkqQdnUSiloLeY3kZyXx2WE8WLi/jm2cPJ9kX3kDdvq+ehcvKWLh8F7tqGslLT+K6UwczdfxArX8XiREK+ghw9aRBXDt7KS+v28OFY/qG/Hx1zW38bXU5C5aV8d62/Zh1rgC6/ZwRfH5Ub1KTomupp4h8MgV9BPjc8AIG5qcz+5/buOD4PiG5KNva3sFbpVU8t3I3L63dQ2NrO8cWZPCdz4/g4nH9dHOTSAxT0EeAhATjhtMG88Nn1/H6pgrOHFkYlK/b1t7Bkq37+euq3by0bg81Da1kpfr44rh+fPnE/pwwIDfiVvqISPAp6CPEZScN5JG3t/Hff9vIacMKun3xs73DsWx7Nc+v3s3f1pRTVddCRnIiZxUVcuHxfTlteM+ouwtXRAKjoI8Qyb4Evn/ecdzwaAm/euV9Zp478ohf29jSzlulVby6fg+LNlSwr76F1KQEJo8s5ILj+3DGyF6adxeJYwr6CHJWUSFTxw/kT//YTL/cVK6aNOhjn1tV18zfN1Twyvq9vFVaSVNrB1mpPs4Y0Yuziwo5c2QvMtSSQERQ0EecH11URMXBJu58dh0l26u5/tTBDCnIpLG1nQ/21lHd/xQacwdz0s9ewznom5PKpcUDOLuoN+MH54d9eaaIRD5zznldA8XFxa6kpMTrMiJGa3sHv1/0Afe/sYXmto5//0PXQUpdOTd98XTOOq6QUX2zdUFVJE6Z2TLnXPGnPk9BH7n217fw5geV7DnQRFJiAscWZHDnV68isb2ZxYsXe12eiHjsSINeUzcRLD8jmSlj+/3bscT2Zo+qEZFopQldEZEYp6AXEYlxCnoRkRgXkqA3s3PMbJOZlZrZzFCcQ0REjkzQg97MEoH/Bc4FioCpZlYU7POIiMiRCcWIfjxQ6pzb4pxrAZ4EpoTgPCIicgRCsbyyH7Czy+MyYMJHn2RmM4AZ/od1Zrapm+frCVR187XRqqeZxd33TBz+PaPvOR4E8j0fcyRP8mwdvXNuFjAr0K9jZiVHcsNALNH3HB/0PceHcHzPoZi62QUM6PK4v/+YiIh4IBRBvxQYZmaDzSwZuAx4LgTnERGRIxD0qRvnXJuZ3QK8DCQCDzvn1gX7PF0EPP0ThfQ9xwd9z/Eh5N9zRDQ1ExGR0NGdsSIiMU5BLyIS46I66OOt1YKZDTCz181svZmtM7PbvK4pHMws0cxWmNnzXtcSDmaWa2YLzGyjmW0ws0le1xRqZvYN/8/0WjObZ2apXtcUbGb2sJlVmNnaLsfyzexVM/vA/zEvFOeO2qCP01YLbcC3nHNFwETg5jj4ngFuAzZ4XUQY/RZ4yTk3EhhDjH/vZtYPuBUods6NpnMRx2XeVhUSs4FzPnJsJrDIOTcMWOR/HHRRG/TEYasF51y5c265//NaOgOg3ye/KrqZWX/gfOBBr2sJBzPLAT4LPATgnGtxztV4W1VY+IA0M/MB6cBuj+sJOufcG8D+jxyeAszxfz4H+EIozh3NQX+4VgsxHXpdmdkg4ARgibeVhNxvgNuBjk97YowYDFQCj/inqx40swyviwol59wu4JfADqAcOOCce8XbqsKm0DlX7v98D1AYipNEc9DHLTPLBBYCX3fOHfS6nlAxswuACufcMq9rCSMfMA64zzl3AlBPiN7ORwr/vPQUOn/J9QUyzOxKb6sKP9e51j0k692jOejjstWCmSXRGfJznXPPeF1PiJ0CXGRm2+icmjvTzB73tqSQKwPKnHOH3qktoDP4Y9lZwFbnXKVzrhV4BjjZ45rCZa+Z9QHwf6wIxUmiOejjrtWCmRmdc7cbnHP3el1PqDnn7nDO9XfODaLz7/fvzrmYHuk55/YAO81shP/QZGC9hyWFww5gopml+3/GJxPjF6C7eA64xv/5NcCzoTiJZ90rA+VBq4VIcApwFbDGzFb6j33POfc3D2uS4PsaMNc/gNkCXOtxPSHlnFtiZguA5XSuLFtBDLZCMLN5wOl0thkvA+4C7gaeMrPrge3AJSE5t1ogiIjEtmieuhERkSOgoBcRiXEKehGRGKegFxGJcQp6EZEYp6AXEYlxCnoRkRj3fyctfid7bk+HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(x,y + penalty)\n",
    "plt.vlines(2, 0, 100)\n",
    "plt.ylim([0,100])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function has minimum at roughly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.8918918918918919"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x[np.argmin(y+penalty)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's not ideal, but:\n",
    "- 1.89 is a better answer than 0.\n",
    "- you can make it work much better by choosing a different penalty function.\n",
    "- if you consider only the integers, the minimum is actually at 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can read more about this approach here: https://www.rose-hulman.edu/~bryan/lottamath/penalty.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Going back to our case we have some constraints and we want to create a penalty function which will increase cost function for all the solution which don't satisfy the constraints. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Matrix representation of gates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, we will find the matrix representation of Pauli operators quite useful, so let's see how they work before we go deep into encoding the constraints.\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far, we used the `Z` operators for assigning costs. Take a look at the following code from the MaxCut example.\n",
    "```\n",
    "for i, j in graph.edges():\n",
    "    cost_operators.append(PauliTerm(\"Z\", i, 0.5)*PauliTerm(\"Z\", j) + PauliTerm(\"I\", 0, -0.5))\n",
    "```\n",
    "\n",
    "What we wanted to achieve? We wanted to mimic the classical version of adding +1 if it's inside our solution or adding -1 if it's not.\n",
    "<br>\n",
    "So in order to do that we need to find the combination of pauli operators, which will add +1 for the correct solutions and add -1 for the incorrect ones. +1/-1 values will be of course parametrized.\n",
    "<br>\n",
    "It will be helpful to see how the matrix representation of Pauli operators looks like."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single 'Z' matrix is given by the following matrix:\n",
    "<br>\n",
    "$Z_1 = \\begin{bmatrix}\n",
    "1 & 0 \\\\\n",
    "0 & -1 \\\\\n",
    "\\end{bmatrix}$\n",
    "<br>\n",
    "\n",
    "First column/row represents $\\ket{0}$ and the second one $\\ket{1}$.\n",
    "So if we apply $Z_1$ to given quantum state, it will change $\\ket{0}$ into $\\ket{0}$ (so no change) and $\\ket{1}$ into $-\\ket{1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we do an operation $Z_1 Z_2$, we actually do $Z_1 \\otimes Z_2$, which is a tensor multiplication. It's pretty straightforward to calculate, so don't be scared of the \"tensor\" word here. We just create a new matrix (in our case 4x4) where every element of the first matrix is multiplied by the second matrix. It's easier to get it with an example:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given:\n",
    "    $\\begin{bmatrix}\n",
    "    a_{1,1} & a_{1,2} \\\\\n",
    "    a_{2,1} & a_{2,2} \\\\\n",
    "  \\end{bmatrix}, \\qquad \\begin{bmatrix}\n",
    "    b_{1,1} & b_{1,2} \\\\\n",
    "    b_{2,1} & b_{2,2} \\\\\n",
    "  \\end{bmatrix}, $\n",
    "\n",
    "The tensor product of these two is:\n",
    "\n",
    "$\n",
    "  \\begin{bmatrix}\n",
    "    a_{1,1} & a_{1,2} \\\\\n",
    "    a_{2,1} & a_{2,2} \\\\\n",
    "  \\end{bmatrix}\n",
    "\\otimes\n",
    "  \\begin{bmatrix}\n",
    "    b_{1,1} & b_{1,2} \\\\\n",
    "    b_{2,1} & b_{2,2} \\\\\n",
    "  \\end{bmatrix}\n",
    "=\n",
    "  \\begin{bmatrix}\n",
    "    a_{1,1}  \\begin{bmatrix}\n",
    "              b_{1,1} & b_{1,2} \\\\\n",
    "              b_{2,1} & b_{2,2} \\\\\n",
    "            \\end{bmatrix} & a_{1,2}  \\begin{bmatrix}\n",
    "                                      b_{1,1} & b_{1,2} \\\\\n",
    "                                      b_{2,1} & b_{2,2} \\\\\n",
    "                                    \\end{bmatrix} \\\\\n",
    "     & \\\\\n",
    "    a_{2,1}  \\begin{bmatrix}\n",
    "              b_{1,1} & b_{1,2} \\\\\n",
    "              b_{2,1} & b_{2,2} \\\\\n",
    "            \\end{bmatrix} & a_{2,2}  \\begin{bmatrix}\n",
    "                                      b_{1,1} & b_{1,2} \\\\\n",
    "                                      b_{2,1} & b_{2,2} \\\\\n",
    "                                    \\end{bmatrix} \\\\\n",
    "  \\end{bmatrix}\n",
    "=\n",
    "  \\begin{bmatrix}\n",
    "    a_{1,1} b_{1,1} & a_{1,1} b_{1,2} & a_{1,2} b_{1,1} & a_{1,2} b_{1,2} \\\\\n",
    "    a_{1,1} b_{2,1} & a_{1,1} b_{2,2} & a_{1,2} b_{2,1} & a_{1,2} b_{2,2} \\\\\n",
    "    a_{2,1} b_{1,1} & a_{2,1} b_{1,2} & a_{2,2} b_{1,1} & a_{2,2} b_{1,2} \\\\\n",
    "    a_{2,1} b_{2,1} & a_{2,1} b_{2,2} & a_{2,2} b_{2,1} & a_{2,2} b_{2,2} \\\\\n",
    "  \\end{bmatrix}\n",
    "$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It may seem tedious, but we deal mainly with zeros and ones, so it's not that bad :)\n",
    "<br>\n",
    "So back to our Pauli operators. If we do $Z_1 \\otimes Z_2$ we have:\n",
    "\n",
    "$\n",
    "Z_1 \\otimes Z_2 = \n",
    "  \\begin{bmatrix}\n",
    "    1 & 0 \\\\\n",
    "    0 & -1 \\\\\n",
    "  \\end{bmatrix}\n",
    "\\otimes\n",
    "  \\begin{bmatrix}\n",
    "    1 & 0 \\\\\n",
    "    0 & -1 \\\\\n",
    "  \\end{bmatrix}\n",
    "=\n",
    "  \\begin{bmatrix}\n",
    "    1  \\begin{bmatrix}\n",
    "              1 & 0 \\\\\n",
    "              0 & -1 \\\\\n",
    "            \\end{bmatrix} & 0  \\begin{bmatrix}\n",
    "                                      1 & 0 \\\\\n",
    "                                      0 & -1 \\\\\n",
    "                                    \\end{bmatrix} \\\\\n",
    "     & \\\\\n",
    "    0  \\begin{bmatrix}\n",
    "              1 & 0 \\\\\n",
    "              0 & -1 \\\\\n",
    "            \\end{bmatrix} & -1  \\begin{bmatrix}\n",
    "                                      1 & 0 \\\\\n",
    "                                      0 & -1 \\\\\n",
    "                                    \\end{bmatrix} \\\\\n",
    "  \\end{bmatrix}\n",
    "=\n",
    "  \\begin{bmatrix}\n",
    "    1 & 0 & 0 & 0 \\\\\n",
    "    0 & -1 & 0 & 0  \\\\\n",
    "    0 & 0 & -1 & 0 \\\\\n",
    "    0 & 0 & 0 & 1 \\\\\n",
    "  \\end{bmatrix}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we once again take a look at our cost operator definition from MaxCut:\n",
    "<br>\n",
    "$H_{cost} = \\sum \\frac{1}{2} w_{ij} (\\mathbb{1} - \\sigma^z_i \\sigma^z_j$),\n",
    "<br>\n",
    "we can now write it in term of matrices. Let's do that for one pair:\n",
    "\n",
    "$\n",
    "\\mathbb{1} - Z_1 \\otimes Z_2 \n",
    "= \n",
    "  \\begin{bmatrix}\n",
    "    1 & 0 & 0 & 0 \\\\\n",
    "    0 & 1 & 0 & 0  \\\\\n",
    "    0 & 0 & 1 & 0 \\\\\n",
    "    0 & 0 & 0 & 1 \\\\\n",
    "  \\end{bmatrix}\n",
    "-\n",
    "  \\begin{bmatrix}\n",
    "    1 & 0 & 0 & 0 \\\\\n",
    "    0 & -1 & 0 & 0  \\\\\n",
    "    0 & 0 & -1 & 0 \\\\\n",
    "    0 & 0 & 0 & 1 \\\\\n",
    "  \\end{bmatrix}\n",
    "=\n",
    "  \\begin{bmatrix}\n",
    "    0 & 0 & 0 & 0 \\\\\n",
    "    0 & 2 & 0 & 0  \\\\\n",
    "    0 & 0 & 2 & 0 \\\\\n",
    "    0 & 0 & 0 & 0 \\\\\n",
    "  \\end{bmatrix}\n",
    "$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means, that having states $\\ket{00}$ and $\\ket{11}$ is associated with 0, so it doesn't add any value to our cost function. And if we have state $\\ket{01}$ or $\\ket{10}$, it does."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the beginning it might be confusing which row/colum represents which state. For 2x2 matrix it's (in order) $\\ket{00}$, $\\ket{01}$, $\\ket{10}$, $\\ket{11}$ - if you don't know why, follow carefully which states are getting multiplied by which in the tensor product."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can experiment with it using python using the `np.kron` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 0., 0., 0.],\n",
       "       [0., 2., 0., 0.],\n",
       "       [0., 0., 2., 0.],\n",
       "       [0., 0., 0., 0.]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "Z = np.array([[1,0],[0,-1]])\n",
    "np.eye(4) - np.kron(Z,Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating penalty matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does it mean for us?\n",
    "<br>\n",
    "That if we need to find a way to construct such matrix using the Pauli gates that it will add huge penalty for the states we don't want."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned earlier, there are two types of constraints we have:\n",
    "- we need to visit each city once\n",
    "- in every moment we need to be in exactly one city\n",
    "\n",
    "Let's take an example of 3 cities to make it clear how to encode these."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our solution will be the following vector state: \n",
    "<br>\n",
    "$\\ket{\\psi_{sol}} = \\ket{q_1 q_2 q_3 q_4 q_5 q_6 q_7 q_8 q_9}$\n",
    "<br>\n",
    "Let's say, that at $t=0$ we want to visit only one city. This means, that the only possible states are $\\ket{100}$, $\\ket{010}$ and $\\ket{001}$, so we want to penalize all the other states.\n",
    "It turns out, that the following operator almost gives us what we want:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Product of Z1, Z2 and Z3\n",
      "[[ 1  0  0  0  0  0  0  0]\n",
      " [ 0 -1  0  0  0  0  0  0]\n",
      " [ 0  0 -1  0  0  0  0  0]\n",
      " [ 0  0  0  1  0  0  0  0]\n",
      " [ 0  0  0  0 -1  0  0  0]\n",
      " [ 0  0  0  0  0  1  0  0]\n",
      " [ 0  0  0  0  0  0  1  0]\n",
      " [ 0  0  0  0  0  0  0 -1]]\n",
      "Our operator\n",
      "[[1. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "Z1Z2Z3 = np.kron(Z,np.kron(Z,Z))\n",
    "print(\"Product of Z1, Z2 and Z3\")\n",
    "print(Z1Z2Z3)\n",
    "print(\"Our operator\")\n",
    "basic_penalty = 0.5 * (np.eye(8) + Z1Z2Z3)\n",
    "print(basic_penalty)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This operator ignore the states we want and  will add penalty to all the others, namely: $\\ket{000}$, $\\ket{011}$, $\\ket{101}$ and $\\ket{110}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only problem is, that it doesn't penalize state $\\ket{111}$, but we can deal with that with the following addition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our helper operator:\n",
      "[[0. 0.]\n",
      " [0. 1.]]\n",
      "Our operator for 111 state:\n",
      "[[0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "helper_op = 0.5 * (np.eye(2) - Z)\n",
    "print(\"Our helper operator:\")\n",
    "print(helper_op)\n",
    "print(\"Our operator for 111 state:\")\n",
    "penalty_111 = np.kron(helper_op, np.kron(helper_op,helper_op))\n",
    "print(penalty_111)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can combine those operators to get our final penalty operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final penalty operator\n",
      "[[1. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 0. 1.]]\n"
     ]
    }
   ],
   "source": [
    "print(\"Final penalty operator\")\n",
    "final_penalty = basic_penalty + penalty_111\n",
    "print(final_penalty)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, so now we know how to deal with our very first constraint - being in only one city at the time, for the first time slot. We need to add penalty operator for qubits $q_1$, $q_2$, $q_3$, $P_{123}$, to our cost operators. With a proper weight, of course.\n",
    "<br>\n",
    "To deal with all the other time slots, we just need to work with other set of qubits, so have $P_{456}$ and $P_{789}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we want to deal with the second constraint, that we can visit each city only once. Well, it's actually exactly the same situation, but this time instead of using qubits associated with the same time slot, we need to use qubits associated with the same cities. \n",
    "<br>\n",
    "So operators: $P_{147}$, $P_{258}$ and $P_{369}$ will do the job."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to actually code it in pyQuil and add it to our cost operator."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementing penalties in QAOA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have all the information to code the penalty operator by yourself and apply it to different qubits. \n",
    "<br>\n",
    "Below are the empty functions that you should fill in.\n",
    "<br>\n",
    "If you encounter problems you can check the code in `scripts/forest_tsp_solver_naive.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Please implement function which returns an operator corresponding\n",
    "# to the final_penalty matrix we have defined above.\n",
    "# Qubits should be a list of integers, e.g. [0, 1, 2]\n",
    "def create_penalty_operators_for_qubits(qubits):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function should return a list of operators which penalize visiting\n",
    "# a city more than once. \n",
    "# You should use create_penalty_operators_for_qubits function.\n",
    "def create_penalty_operators_for_repetition(number_of_cities):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function should return a list of operators which penalize visiting\n",
    "# more than one city in a given time slot.\n",
    "# You should use create_penalty_operators_for_qubits function.\n",
    "def create_penalty_operators_for_bilocation(number_of_cities):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function should return a list of operators which \n",
    "# encode the distance matrix.\n",
    "# You can reuse the code we wrote earlier.\n",
    "def create_weights_cost_operators(distance_matrix):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function should return all the cost operators required to solve TSP.\n",
    "# It should have a part corresponding to the distance matrix\n",
    "# and to the penalties.\n",
    "def create_cost_operators(distance_matrix):\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Utilising new cost operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the cost operators with the penalties ready, we try using it "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: Fast method for expectation will be used. Noise\n",
      "                     models will be ineffective\n",
      "\tParameters: [2.36630561 2.87194104] \n",
      "\tE => -2659.729526109549\n",
      "\tParameters: [2.36630561 2.87194104] \n",
      "\tE => -2500.928153946556\n",
      "\tParameters: [2.36630561 2.87194104] \n",
      "\tE => -2649.549671687539\n",
      "\tParameters: [2.36630561 2.87194104] \n",
      "\tE => -2656.8578001498586\n",
      "\tParameters: [2.36630561 2.87194104] \n",
      "\tE => -2593.2672286683924\n",
      "\tParameters: [2.36630561 2.87194104] \n",
      "\tE => -2602.6375162862882\n",
      "\tParameters: [2.36630561 2.87194104] \n",
      "\tE => -2489.3044561287666\n",
      "\tParameters: [2.36977187 2.87250197] \n",
      "\tE => -2717.6451379467208\n",
      "\tParameters: [2.36977187 2.87250197] \n",
      "\tE => -2717.1651827131263\n",
      "\tParameters: [2.36977187 2.87250197] \n",
      "\tE => -2700.167597643281\n",
      "\tParameters: [2.36780766 2.87278243] \n",
      "\tE => -2690.609142205864\n",
      "\tParameters: [2.37150501 2.87278243] \n",
      "\tE => -2721.3134501326062\n",
      "\tParameters: [2.37150501 2.87278243] \n",
      "\tE => -2728.741501800037\n",
      "\tParameters: [2.37150501 2.87278243] \n",
      "\tE => -2731.8585776181585\n",
      "\tParameters: [2.36967078 2.87274737] \n",
      "\tE => -2732.3899979107464\n",
      "\tParameters: [2.37614836 2.87272984] \n",
      "\tE => -2732.780241912009\n",
      "\tParameters: [2.37614836 2.87272984] \n",
      "\tE => -2732.622686881663\n",
      "\tParameters: [2.38319193 2.8727408 ] \n",
      "\tE => -2733.439787870112\n",
      "\tParameters: [2.38319193 2.8727408 ] \n",
      "\tE => -2733.0264291847257\n",
      "\tParameters: [2.40319067 2.8727167 ] \n",
      "\tE => -2733.7516172683777\n",
      "\tParameters: [2.40530789 2.872766  ] \n",
      "\tE => -2734.5837580210487\n",
      "\tParameters: [2.40530789 2.872766  ] \n",
      "\tE => -2734.3407060767177\n",
      "\tParameters: [2.42742385 2.8727912 ] \n",
      "\tE => -2733.7879045368604\n",
      "\tParameters: [2.42083625 2.87276025] \n",
      "\tE => -2734.7191386661075\n",
      "\tParameters: [2.41471897 2.87277086] \n",
      "\tE => -2734.7560369013227\n",
      "\tParameters: [2.42260073 2.87277837] \n",
      "\tE => -2734.7757297740714\n",
      "\tParameters: [2.42260073 2.87277837] \n",
      "\tE => -2734.772663201975\n",
      "\tParameters: [2.41794668 2.87277188] \n",
      "\tE => -2734.779095761112\n",
      "The solution is: [2, 2, 0, 1, 2]\n"
     ]
    }
   ],
   "source": [
    "tsp_solver = ForestTSPSolverNaive(distance_matrix, use_constraints=True)\n",
    "solution, naive_distribution = tsp_solver.solve_tsp()\n",
    "print(\"The solution is:\", solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might got a proper solution or not, but it's good to see what's the actual distribution of the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAD8CAYAAADezxtfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHfdJREFUeJzt3X+MHfW53/H3p+baldOwgH2vs7JR7TRbkCP6g2wIkW4jNSRgp20MKUhGVVgkbqyWWG6EkLLRTQRKGhVIGlwE4cYYClQNxqGtWOIE13W4IkGYeA2JjeO6XhwTrwW4wbBAUCCYp3/MszAczu45u97d77H9eUmrPfPMd555zmi9j2fOd2cUEZiZmZXwd0oXYGZmJy83ITMzK8ZNyMzMinETMjOzYtyEzMysGDchMzMrxk3IzMyKcRMyM7Ni3ITMzKyYU0oX0Onmz58fixcvLl2GmdlxY8eOHb+PiD9vZ6ybUAuLFy9mcHCwdBlmZscNSc+2O9aX48zMrBg3ITMzK8ZNyMzMinETMjOzYtyEzMysGDchMzMrxk3IzMyKcRMyM7Ni/MeqLew6NMLi/k2lyzAzmzEHbvgXM7YvnwmZmVkxbkJmZlaMm5CZmRXjJmRmZsW01YQkLZO0V9KQpP6MLZH0RMbulzQ745+S9KSktyRd2pCnT9K+/OqrxT8maVfmukWSMn6ZpN2S3pbUWxs/T9Ijkl6TdGvDPr4t6aCk1xriY9ZlZmZltGxCkmYBtwHLgaXA5ZKWAjcCN0fER4CXgKtyk98BVwI/bMhzBnAd8AngPOA6Safn6tuBLwE9+bUs408DXwAebSjrj8A3gGublPxQ5m/UtC4zMyunnTOh84ChiNgfEW8CG4AVwKeBB3LMPcDFABFxICJ2Am835LkI2BIRRyLiJWALsExSN3BqRGyLiADureXaExF7GwuKiD9ExC+omlHjum0R8VyT+Fh1mZlZIe00oYXAwdrycMZejoi3GmKTybMwXzfGi5G0StKgpMGjr4+ULMXM7ITmiQlNRMS6iOiNiN5Zc7tKl2NmdsJqpwkdAs6sLS/K2GmSTmmITSbPoXzdGDczsxNcO01oO9CTs+FmAyuBAeARYHSWWR/wYIs8m4ELJZ2eExIuBDbn5zevSDo/Z8Vd0UYuMzM7AbRsQvm5z2qqJrIH2BgRu4GvAtdIGgLmAXcCSPq4pGHgMuAHknZnniPAt6ia2nbgmxkDuBpYDwwBzwA/zVyXZK5PApskbR6tS9IB4HvAlZKGc8Yekm7KbeZm/Prx6jIzs3JUTUizsczp7onuvrWlyzAzmzHHegNTSTsiorf1SE9MMDOzgtyEzMysGD9PqIVzFnYxOIPP1jAzO5n4TMjMzIpxEzIzs2LchMzMrBh/JtTCrkMjLO7fVLoMM+twxzqt+WTlMyEzMyvGTcjMzIpxEzIzs2LchMzMrJi2m5CkZZL2ShqS1J+x1bkckubXxkrSLblup6Rza+v6JO3Lr75a/NuSDkp6rWG/cyTdn7mekLQ44/MkPSLpNUm3NmzzMUm7cptb8u7cSDpD0pbc95ba48XNzKyAtpqQpFnAbcByYClwed61+jHgM8CzDZssB3ryaxVwe+Y5A7gO+ATVY8OvqzWChzLW6CrgpYj4CHAzcGPG/wh8A7i2yTa3A1+q1bAs4/3A1ojoAbbmspmZFdLumdB5wFBE7I+IN4ENwIqIeCoiDjQZvwK4NyrbqB6A1w1cBGyJiCMR8RKwhWwQEbEtny3ULNc9+foB4AJJiog/RMQvqJrRO3I/p2a+AO4FLm6S655a3MzMCmi3CS0EDtaWhzM20fETzfOeXPlsoxGq5xeNN354jH0sqDW654EFzRJIWiVpUNLg0ddHWpRnZmaTddJOTMizpKYPU4qIdRHRGxG9s+Z2zXBlZmYnj3ab0CHgzNryooxNdPxE87wnl6RTgC7gxRbjF42xjxfyct3oZbvDLfZtZmbTqN0mtB3okbRE0mxgJTAwzvgB4IqcJXc+MJKXwTYDF0o6PSckXJix8QwAo7PoLgV+FuM8Djb384qk83NW3BXAg01y9dXiZmZWQFtNKD+LWU3VMPYAGyNit6Q1koapzjZ2Slqfm/wE2A8MAXcAV2eeI8C3qJraduCbGUPSTZlrrqRhSddnrjuBeZKGgGuozWiTdAD4HnBlbrM0V10NrM/9PwP8NOM3AJ+VtI9qVt8NbR0lMzObFhrnpMKAOd090d23tnQZZtbhfAPTd0naERG97Yw9aScmmJlZeW5CZmZWjJuQmZkV44fatXDOwi4Gfa3XzGxa+EzIzMyKcRMyM7Ni3ITMzKwYfybUwq5DIyzu31S6DLNj5r9jsU7kMyEzMyvGTcjMzIpxEzIzs2Lafbz3Mkl7JQ1J6s/Y6lwOSfNrYyXplly3U9K5tXUPS3pZ0o8b8k8mV5+kffnVV4t/TNKu3OaWvJM2ki6TtFvS25LauqeRmZlNr5ZNSNIs4DZgObAUuDzvVv0Y1Z2on23YZDnQk1+rgNtr674DfLHJbiaUS9IZwHXAJ6gePX5dPhqCHPOl2nbLMv408AXg0Vbv2czMZkY7Z0LnAUMRsT8i3gQ2ACsi4qmIONBk/Arg3qhsA04bfZBcRGwFXm3cYBK5LgK2RMSRiHgJ2AIsy3WnRsS2fObQvcDFuY89EbG3jfdrZmYzpJ0mtBA4WFsezthUjZ/MvseLD0/Rvs3MbJp5YkITklZJGpQ0ePT1kdLlmJmdsNppQoeAM2vLizI2VeMns+/x4ouOdd8RsS4ieiOid9bcrgkXbWZm7WmnCW0HeiQtkTQbWAkMjDN+ALgiZ7adD4xExHOTrG+sXJuBCyWdnhMSLgQ257pXJJ2fs+KuAB6c5L7NzGyatWxCEfEWsJrqF/8eYGNE7Ja0RtIw1dnGTknrc5OfAPuBIeAO4OrRXJJ+DvwIuEDSsKSLMj6hXBFxBPgWVYPcDnwzY+SY9bnNM8BPcx+X5D4+CWyStHlCR8rMzKacqklkNpY53T3R3be2dBlmx8z3jrOZImlHRLT195iemGBmZsW4CZmZWTFuQmZmVoyfJ9TCOQu7GPS1dDOzaeEzITMzK8ZNyMzMinETMjOzYvyZUAu7Do2wuH9T6TLMJsx/F2THA58JmZlZMW5CZmZWjJuQmZkV03YTkrRM0l5JQ5L6M7Y6l0PS/NrYsyU9LukNSde2ypPxJZKeyPj9ecduJM3J5aFcv7i2zdcyvnf0ZqgZv0vSYUlPN+z7DElbJO3L76djZmbFtNWEJM0CbgOWA0uByyUtBR4DPgM827DJEWAN8N028wDcCNwcER8BXgKuyvhVwEsZvznHkdutBD4KLAO+n/kB7s5Yo35ga0T0AFtz2czMCmn3TOg8YCgi9kfEm8AGYEVEPBURBxoHR8ThiNgO/KmdPPnsn08DD+S4e4CL8/WKXCbXX5DjVwAbIuKNiPgt1aMbzsv9P0rVCBvVc9X3YWZmBbTbhBYCB2vLwxmbqLHyzANezmcXNeZ/Z5tcP5LjJ1PTgtoD9p4HFkziPZiZ2RQ5aScmRPUgpaYPU5K0StKgpMGjr4/McGVmZiePdpvQIeDM2vKijE3UWHleBE6TdEpD/D3b5PquHD+Zml6Q1J25uoHDzQZFxLqI6I2I3llzu9p8a2ZmNlHtNqHtQE/OYJtNNSFgYBL7a5onz0oeAS7NcX3Ag/l6IJfJ9T/L8QPAypw9twToAX7ZYv/1XPV9mJlZAW01ofwsZjWwGdgDbIyI3ZLWSBqmOgvZKWk9gKQPZfwa4OuShiWdOlae3M1XgWskDVF95nNnxu8E5mX8GnJGW263EfgN8DDw5Yg4mvu/D3gcOCv3PTrT7gbgs5L2Uc3qu2GiB8zMzKaOqpMKG8uc7p7o7ltbugyzCfO946wUSTsioredsSftxAQzMyvPTcjMzIpxEzIzs2L8PKEWzlnYxaCvrZuZTQufCZmZWTFuQmZmVoybkJmZFePPhFrYdWiExf2bSpdhNi7/TZAdr3wmZGZmxbgJmZlZMW5CZmZWTLuP975L0mFJT9diZ0jaImlffj8942dLelzSG5KubZWnRS5JukXSkKSdks6tbfOwpJcl/bgh1+ocH5Lm1+Jj5jIzszLaPRO6G1jWEOsHtkZED7A1l6F6rPYa4Ltt5hkv13KqRzT0AKuA22vbfAf4YpNcj1HdIfvZhvh4uczMrIB2H+XwKFVzqVsB3JOv7wEuzrGHI2I78Kc284yZK+P3RmUb1YPvujPXVuDVJvt4KiIOjLGPprnMzKyMY/lMaEFEPJevnwcWTEOuhcDB2rjhjE3GVOYyM7MpMCUTE/JJp1PyYKKpzDVZklZJGpQ0ePT1kZKlmJmd0I6lCb0wejkrvx+ehlyHgDNr4xZlbDLazhUR6yKiNyJ6Z83tmuTuzMyslWNpQgNAX77uAx6chlwDwBU5s+18YKR22W4y+5iqXGZmNgXanaJ9H/A4cJakYUlXATcAn5W0j2o22g059kOShoFrgK/n+FPHycNYuYCfAPuBIeAO4OpaTT8HfgRckLkuyvia3P8iYKek9a1ymZlZGao+grGxzOnuie6+taXLMBuX7x1nnUTSjojobWes75hgZmbFuAmZmVkxbkJmZlaMnyfUwjkLuxj09XYzs2nhMyEzMyvGTcjMzIpxEzIzs2L8mVALuw6NsLh/U+kyzJry3wfZ8c5nQmZmVoybkJmZFeMmZGZmxbgJmZlZMe3eRfsuSYclPV2LXSZpt6S3JfU2jP+apCFJe0fvbj1WnmPItSxjQ5L6a/HVGQtJ82vxsyU9LukNSde2877NzGx6tXsmdDewrCH2NPAF4NF6UNJSYCXw0dzm+5JmjZNnwrky323AcmApcHmOBXiM6nEQzzbs4wiwBvhuy3drZmYzoq0mFBGPUv0Sr8f2RMTeJsNXABsi4o2I+C3V83vOGyvPJHOdBwxFxP6IeBPYkGOJiKci4kCTfRyOiO3An9p5z2ZmNv2m4zOhhcDB2vJwxqYy11Tu430krZI0KGnw6OsjU5XWzMwaeGJCExGxLiJ6I6J31tyu0uWYmZ2wpuOOCYeAM2vLizI21bmmah9mZlbIdJwJDQArJc2RtAToAX45xbm2Az2SlkiaTTV5YWAKajczsxnU7hTt+4DHgbMkDUu6StIlkoaBTwKbJG0GiIjdwEbgN8DDwJcj4uhYeTI+oVwR8RawGtgM7AE25lgkrclci4CdktZn/EMZvwb4eu7/1GM5eGZmdmwUEaVr6Ghzunuiu29t6TLMmvINTK0TSdoREb2tR3pigpmZFeQmZGZmxfh5Qi2cs7CLQV/yMDObFj4TMjOzYtyEzMysGDchMzMrxp8JtbDr0AiL+zeVLsOsKU/RtuOdz4TMzKwYNyEzMyvGTcjMzIpxEzIzs2LavYHpXZIOS3q6FrtM0m5Jb0vqrcXnSXpE0muSbm3I821JByW91hCfI+l+SUOSnpC0uLbuaxnfK+mi8WrK+BmStkjal99Pz7gk3ZK5dko6t71DZGZm06XdM6G7gWUNsaeBLwCPNsT/CHwDuLZJnofIR303uAp4KSI+AtwM3AggaSnVYxo+mvv/vqRZ49QE0A9sjYgeYGsuAyynehRED7AKuL3JtmZmNoPaakIR8ShwpCG2JyL2Nhn7h4j4BVUzaly3LSKea7KLFcA9+foB4AJJyviGiHgjIn4LDJFNrFlNTXLdA1xci98blW3AaZK6x3vfZmY2vTrlM6GFwEGAfFbQCDCvHk/DGRvPglqjex5Y0LiPVrkkrZI0KGnw6OsjE3kfZmY2AZ3ShKZFVA9LmvADkyJiXUT0RkTvrLld01CZmZlB5zShQ8CZAJJOAbqAF+vxtChj43lh9DJbfj/cuI8J5DIzs2nUKU1oAOjL15cCP8uzmAFgZc6eW0I1qeCXE8jVBzxYi1+Rs+TOB0bG+HzKzMxmSLtTtO8DHgfOkjQs6SpJl0gaBj4JbJK0uTb+APA94MocvzTjN+U2czN+fW5yJzBP0hBwDTmjLSJ2AxuB3wAPA1+OiKNj1ZS5bgA+K2kf8JlcBvgJsJ9qcsMdwNUTOVBmZjb1VJ1w2FjmdPdEd9/a0mWYNeUbmFonkrQjInpbj+ycy3FmZnYSchMyM7Ni/DyhFs5Z2MWgL3mYmU0LnwmZmVkxbkJmZlaMm5CZmRXjz4Ra2HVohMX9m0qXYdaUp2jb8c5nQmZmVoybkJmZFeMmZGZmxbgJmZlZMS2bkKS7JB2W9HQtdpmk3ZLeltRbi8+T9Iik1yTd2pDn25IOSnqtIf4pSU9KekvSpQ3rHpb0sqQfN8RXSxqSFJLm1+JnS3pc0huSrm3Y5n3vw8zMymrnTOhuYFlD7GngC8CjDfE/At8AruX9HiIfzd3gd8CVwA+brPsO8MUm8ceo7pD9bEP8CLAG+G6Tbe7m/e/DzMwKatmEIuJRql/u9dieiNjbZOwfIuIXVM2ocd22Zs/viYgDEbETeLvJuq3Aq03iT0XEgSbxwxGxHfhTO+/DzMzK8mdCZmZWjJtQE5JWSRqUNHj09ZHS5ZiZnbDchJqIiHUR0RsRvbPmdpUux8zshOUmZGZmxbQzRfs+4HHgLEnDkq6SdImkYeCTwCZJm2vjDwDfA67M8UszflNuMzfj12f84xm/DPiBpN21XD8HfgRckNtclPE1uc0iYKek9Rn/UMavAb6e25w61vs4huNmZmZTQBFRuoaONqe7J7r71pYuw6wp38DUOpGkHRHR23qkL8eZmVlBbkJmZlaMm5CZmRXjh9q1cM7CLgZ93d3MbFr4TMjMzIpxEzIzs2LchMzMrBh/JtTCrkMjLO7fVLoMs3f4b4PsROIzITMzK8ZNyMzMinETMjOzYtpuQpKWSdoraUhSf8aWSHoiY/dLmp3xObk8lOsX1/J8LeN7R29ImvG7JB2W9HTDfs+QtEXSvvx+esbPlvS4pDckXduwzftqzfjqjIWk+RM7VGZmNtXaakKSZgG3AcuBpcDleXfsG4GbI+IjwEvA6J2prwJeyvjNOY7cZiXwUWAZ8P3MDXB3xhr1A1sjogfYmstQPap7DfDdNmsFeAz4DPBsO+/bzMymV7tnQucBQxGxPyLeBDYAK4BPAw/kmHuAi/P1ilwm118gSRnfEBFvRMRvgaHMTUQ8StVYGtVzvbOPiDgcEduBP7VZKxHxVEQcaPM9m5nZNGu3CS0EDtaWhzP2ckS81RB7z/hcPwLMGyfPeBZExHP5+nlgwSRrNTOzDnNcTUyI6uFH0/4AJEmrJA1KGjz6+sh0787M7KTVbhM6BJxZW16UsdMkndIQe8/4XN8FvDhOnvG8IKk7c3UDhydZa9siYl1E9EZE76y5XRPZ1MzMJqDdJrQd6MnZcLOpJhcMAI8Al+aYPuDBfD2Qy+T6n+VZzACwMmfPLQF6gF+22Hc9V30fE63VzMw6TFtNKD/XWQ1sBvYAGyNiN/BV4BpJQ1Sf+dyZm9wJzMv4NeSMttxmI/Ab4GHgyxFxFEDSfcDjwFmShiWNzrS7AfispH1UM9tuyPEfkjSc+b+e25w6Tq1IWpPbLAJ2Slo/8UNmZmZTRdUJio1lTndPdPetLV2G2Tt87zjrdJJ2RERvO2OPq4kJZmZ2YnETMjOzYtyEzMysGD9PqIVzFnYx6GvwZmbTwmdCZmZWjJuQmZkV4yZkZmbF+DOhFnYdGmFx/6bSZdhJzH8XZCcynwmZmVkxbkJmZlaMm5CZmRXT7uO9l0naK2lIUn/GVudySJpfGytJt+S6nZLOra17WNLLkn7ckH8yufok7cuvvlr825IOSnqtYR+fkvSkpLckXYqZmRXXsglJmgXcBiwHlgKXS1oKPEZ1V+tnGzZZTvWIhh5gFXB7bd13gC822c2Eckk6A7gO+ATV47yvk3R6bvNQxhr9DrgS+OF479fMzGZOO2dC5wFDEbE/It4ENgArIuKpiDjQZPwK4N6obKN68F03QERsBV5t3GASuS4CtkTEkYh4CdgCLMtc22qPA6/v40BE7ATebuM9m5nZDGinCS0EDtaWhzM2VeMns++p3IeZmRXiiQlNSFolaVDS4NHXR0qXY2Z2wmqnCR0CzqwtL8rYVI2fzL6nch/vExHrIqI3Inpnze2aqrRmZtagnSa0HeiRtETSbGAlMDDO+AHgipzZdj4w0uwzmjaNlWszcKGk03NCwoUZMzOz40jLJhQRbwGrqX7J7wE2RsRuSWskDVOdheyUtD43+QmwHxgC7gCuHs0l6efAj4ALJA1LuijjE8oVEUeAb1E1yO3ANzOGpJsy19zcx/UZ/3jGLwN+IGn3xA+XmZlNJUVE6Ro62pzunujuW1u6DDuJ+d5xdryRtCMietsZ64kJZmZWjJuQmZkV4yZkZmbF+HlCLZyzsItBX5M3M5sWPhMyM7Ni3ITMzKwYNyEzMyvGTcjMzIpxEzIzs2LchMzMrBg3ITMzK8ZNyMzMinETMjOzYnwX7RYkvQrsLV3HBM0Hfl+6iElw3TPLdc+c47FmmHzdfz8i/rydgb5tT2t7270leaeQNHi81Qyue6a57plzPNYMM1O3L8eZmVkxbkJmZlaMm1Br60oXMAnHY83gumea6545x2PNMAN1e2KCmZkV4zMhMzMr5oRvQpKWSdoraUhSf5P1cyTdn+ufkLS4tu5rGd8r6aJWOSUtyRxDmXN2p9Qt6UxJj0j6jaTdkv59bfz1kg5J+lV+fa5T6s74AUm7srbBWvwMSVsk7cvvp3dCzZLOqh3LX0l6RdJXcl3xYy1pXv4svCbp1oZtPpbHekjSLZKU8Sk51tNRt6S5kjZJ+j/5s31Dbd2Vkv5f7Xj/VafUnev+NnOO1vcX4+UqXbOkDzb8bP9e0tpcN7ljHREn7BcwC3gG+DAwG/g1sLRhzNXA3+TrlcD9+Xppjp8DLMk8s8bLCWwEVubrvwH+XQfV3Q2cm2M+CPzfWt3XA9d24vHOdQeA+U32dxPQn6/7gRs7peaG/M9T/d1EpxzrDwB/Cfxb4NaGbX4JnA8I+CmwfKqO9XTVDcwF/nm+ng38vFb3lY3vsVPqznV/C/Q22V/TXJ1Qc8P2O4BPHcuxPtHPhM4DhiJif0S8CWwAVjSMWQHck68fAC7I//2tADZExBsR8VtgKPM1zZnbfDpzkDkv7pS6I+K5iHgSICJeBfYACydZ34zV3WJ/9VyTPd7TXfMFwDMR8ewkapuWuiPiDxHxC+CP9cGSuoFTI2JbVL9V7uXdYzoVx3pa6o6I1yPikXz9JvAksGiS9c1Y3S2M9TPXMTVL+ofAX1A1/Uk70ZvQQuBgbXmY9//ifWdMRLwFjADzxtl2rPg84OXMMda+Stb9jjzl/qfAE7Xwakk7Jd11DJdapqvuAP6XpB2SVtXGLIiI5/L188CCDqp51ErgvoZY6WM9Xs7hMXJOxbF+T01N9vG+MW3W/Q5JpwH/CthaC//rPN4PSDqzA+v+L3n56hu1RjPpYzBDNcO7Z0712W0TPtYnehOyBpL+HvDfga9ExCsZvh34B8A/AZ4D/lOh8sbylxFxLrAc+LKkTzUOyH8IHTXVU9Vngp8HflQLd/qxbqkTjzWApFOoGv4tEbE/ww8BiyPiHwFbePd//Z3i30TEOcA/y68vFq5nIhr/gzWpY32iN6FDQL0bL8pY0zH5Q9wFvDjOtmPFXwROyxxj7atk3Uj6M6oG9N8i4n+MDoiIFyLiaES8DdxB68tgM1p3RIx+Pwz8z1p9L+QlpNFLSYc7pea0HHgyIl4YDXTIsR4vZ/0yVj3nVBzr6ap71DpgX0SsHQ1ExIsR8UYurgc+1kl11362XwV+yLs/D5M9BtNec479x8ApEbGj9l4mdaxP9Ca0HehRNWttNlXnHmgYMwD05etLgZ/l//QGgJU5e2QJ0EP1oW3TnLnNI5mDzPlgp9Sdp/l3Ansi4nv1RKO/XNIlwNMdVPcHJH0w6/wAcGGtvnquyR7v6fgZGXU5DZfiOuRYN5WX216RdH7+vFzBu8d0Ko71tNQNIOk/UP0C/UpDvH68P0/1WWhH1C3pFEnz8/WfAf+S5j/bbR2Dmai5ptXPdvvHeqIzGY63L+BzVDPBngH+OmPfBD6fr/8u1eWSIapfIB+ubfvXud1ecrbNWDkz/uHMMZQ553RK3VQzXQLYCfwqvz6X6/4rsCvXDQDdHVT3h6lm9fwa2N1wvOdRXfvfB/xv4IxOqDnjH6D6H2VXw7465VgfAI4Ar1F9VjA6U7KX6hfhM8CtvPsH7VNyrKejbqr/4QfVL73Rn+2/yvH/MX9ufk31n8SzO6juD1DNLtuZNf5n3p0ROmau0j8juW5/47Gc7LH2HRPMzKyYE/1ynJmZdTA3ITMzK8ZNyMzMinETMjOzYtyEzMysGDchMzMrxk3IzMyKcRMyM7Ni/j9Mm8lRvyutbwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "states_with_probs = []\n",
    "for state in tsp_solver.sampling_results:\n",
    "    states_with_probs.append([''.join(map(str,state)), tsp_solver.sampling_results[state]/10000])\n",
    "states_with_probs = sorted(states_with_probs,key=lambda x: -x[1])\n",
    "plots.plot_state_histogram(states_with_probs[:10][::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since it might be hard for you to decode the bitstrings, let's take a look at the human-readable representation.\n",
    "The two might actually differ, since many bitstrings may have the same representation in a case of an invalid solution. That's why we call the distribution coming from tsp_solver a `naive_distribution`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbYAAAD8CAYAAAAbvYHOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+0XWV95/H3Z6JEQoebINAJCdPAIoUJBkOJ6S9WRUEIVIMKtBAkAbXgGlnQhY1EmEbRxqWolAlgJeVHEanxR4GGBooUkoWtRrghaUIkwUjB5GIDWudamgEm8TN/7H1dp4dz7jn358nd+bzWOuvevZ/v+e5n73XIl2c/+55HtomIiKiK/9LpDkRERAynFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiU13W6A/uigw8+2NOmTet0NyIixpR169b9xPYhreJS2Dpg2rRpdHd3d7obERFjiqTn2onLrciIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiUFLaIiKiU/IF2B2zq6WXa4lWd7kZExKh69jO/PyrHyYgtIiIqJYUtIiIqJYUtIiIqJYUtIiIqZVQLm6S5krZK2iZp8SBzLJW0XdJLQ+zL5yRtkbRR0j2SJg4yz/lljk2SviPpzUPpV0REDM2oFTZJ44CbgNOBGcB5kmbUxewn6YAWqe4D5rRxvEktQh4C3mT7OOBp4GODzPMvwFttzwQ+BSxv1beIiBg5ozlimwNss/2M7VeBFcCZdTGTgM2Sbpb0lkZJbK+1/eM2jnevpJWS5kl6zZ812P6W7d3l5lpgapM8iyQ9JukSSQc2yPMd2z9rI09ERIyC0SxsU4DtNds7yn2/ZHsncDSwGlgqab2kyyQdNIjjnQRcB5wNPCXp05KOahL7fuCBRg22rwIuAI4EnpB0u6QTm+T5QLM8ki6W1C2pe8+u3gGcRkREDMRe9/CI7Vdsr7B9KsWI7hTgeUmHDTCPba+xvQA4ATCwRdJZtXGSrgZ2A3f1k2ur7Sspiu7DwCpJy+ryvI2isF3ZJMdy27Ntzx43oWsgpxIREQMwmt880gMcXrM9FXhR0oZye4ntlQCSDqUYJS2gGNnNB3Y2S1zO360rN1faXlLu3x94D8WIbCJwOcXcWt/7LgTeCZxs2+W+24Hjgedtn1HuE/C2Ms8cYBlwS02e48rt023/dCAXJSIihtdoFrbHgemSjqAocucC821f0xcgqQu4AzgGuBM4w3ZPq8S29wCzavdJuhY4B1gFLLK9vq59LvBRigc/dtXkuqgu7nxgCfAkcCuwsDxeX/t/B+4GLrD9dKu+RkTEyBq1wmZ7t6RLgQeBccBttjc3CF0GrO4bQdUrC9Z8YIKkHcAttj/RIHQNxSjw5SZduhEYDzxUDMhYa/tDDeKeA060/WKTPEuANwJfLPPstj27SWxERIwwNakfMYLGT57uyQuv73Q3IiJG1VC/BFnSunYGDnvdwyMRERFDkcIWERGVkvXYOmDmlC66R2ldooiIfU1GbBERUSkpbBERUSkpbBERUSmZY+uATT29TFu8qtPdiIgYsKE+sj8aMmKLiIhKSWGLiIhKSWGLiIhKSWGLiIhKGZOFTdJcSVslbZO0eJA5lkraLumlIfblc5K2SNoo6R5JE4eSLyIihmbMFbZy7bWbgNOBGcB5kmbUxewn6YAWqe6jWFut1fEmtQh5CHiT7eOAp4GPtcoZEREjZ8wVNopitM32M7ZfBVZQrLRdaxKwWdLNkt7SKInttbZ/3Mbx7pW0UtI8Sa/58wjb37K9u9xcS7GAakREdMhYLGxTgO012zvKfb9keydwNLAaWCppvaTLJB00iOOdBFwHnA08JenTko5qEvt+4IFGDZIultQtqXvPrt5BdCMiItoxFgtbW2y/YnuF7VMpRnSnAM9LOmyAeWx7je0FwAmAgS2SzqqNk3Q1sBu4q0me5bZn2549bkLXYE4pIiLaMBYLWw9weM32VOBFSRvK17y+BkmHSvoIxXzaOIqVt3c2SyxpXE2eT9bs31/SfOBu4DTgcoq5tb72C4F3Auc3W/k7IiJGx1j8Sq3HgemSjqAocucC821f0xcgqQu4AzgGuBM4w3ZPq8S29wCzavdJuhY4B1gFLLK9vq59LvBR4K22dw3lxCIiYujGXGGzvVvSpcCDFKOw22xvbhC6DFjdbARVFqz5wARJO4BbbH+iQegaYIntl5t06UZgPPCQJIC1tj80gFOKiIhhNOYKG4Dt+4H7+2nvBR5pkeOjFCOtdo7VX3uzB0kiIqIDxuIcW0RERFMpbBERUSkpbBERUSljco5trJs5pYvuMbBYX0TEWJQRW0REVEoKW0REVEoKW0REVErm2DpgU08v0xav6nQ3ImIf9Ow+ML+fEVtERFRKCltERFRKCltERFTKmCxskuZK2ippm6TFg3j/BEmrJG2RtFnSZ4bQl3PKHL+QNHuweSIiYniMucImaRxwE3A6MAM4T9KMupj9JB3QItXnbR8DHA/8rqTTmxxvUos8TwLvBR5tp/8RETGyxlxhA+YA22w/Y/tVYAXFCtm1JgGbJd0s6S31CWzvsr26/P1V4AmKBUsbuVfSSknzJL3mKVLbT9neOpQTioiI4TMWC9sUYHvN9o5y3y/Z3gkcDawGlkpaL+kySQfVJ5M0EXgX8HCT450EXAecDTwl6dOSslRNRMReaiwWtrbYfsX2CtunUozoTgGel3RYX0w5AvsqsMz2M03y2PYa2wuAEwADWySdNZD+SLpYUrek7j27egd7WhER0cJYLGw9wOE121OBFyVtKF/z+hokHSrpI8B9FKttzwd21rx3OfAD29eX8eNq8nyyJs/+kuYDdwOnAZcDDw2k07aX255te/a4CV0DOuGIiGjfWPzmkceB6ZKOoChy5wLzbV/TFyCpC7gDOAa4EzjDdk9tEkl/BnQBH+zbZ3sPMKsu7lrgHGAVsMj2+pE4qYiIGB5jrrDZ3i3pUuBBilHYbbY3NwhdBqy27foGSVOBq4EtwBOSAG60fUuDPGuAJbZfbtQfSe8BbgAOAVZJ2mD7tIGfWUREDIcxV9gAbN8P3N9Pey/wSD/tOwAN4Fj9td8D3NNOroiIGHljcY4tIiKiqRS2iIiolBS2iIiolDE5xzbWzZzSRfc+sCZSREQnZMQWERGVksIWERGVksIWERGVkjm2DtjU08u0xas63Y2IqLhn99G5/IzYIiKiUlLYIiKiUlLYIiKiUtoqbJLmStoqaZukxYM5kKSlkrZLemkw76/J8zlJWyRtlHRPuVDoYPJI0rLynDZK+o1B5nmjpNWSXpJ042ByRETE8GlZ2CSNA24CTgdmAOdJmlEXs5+kA1qkug+Y08bxJrUIeQh4k+3jgKeBjw0yz+nA9PJ1MfAXg8zzMvCnwJ+0iIuIiFHQzohtDrDN9jO2XwVWUKxIXWsSsFnSzZLe0iiJ7bW2f9zG8e6VtFLSvHKF6/o837K9u9xcS7HQaCOLJD0m6RJJBzZoPxP4crlC9lpgoqTJDeL+UNKTkj4i6ZAG/fkP2/9IUeAiIqLD2ilsU4DtNds7yn2/ZHsncDSwGlgqab2kyyQdNIg+nQRcB5wNPCXp05KOahL7fuCBRg22rwIuAI6kWHPtdkkn1oS0PK8yz5coRncTgEclfbO8NZv5yYiIvdCw/eNs+xXbK2yfSjEaOgV4XtJhA8xj22tsLwBOAAxskXRWbZykq4HdwF395Npq+0qKovswxUKgywZ0YkWe7bY/RXEr9rbyde9Acki6WFK3pO49u3oH2oWIiGhTO3+g3QMcXrM9FXhR0oZye4ntlQCSDqUYJS2gGAHNB3Y2S1zO360rN1faXlLu3x94D8WIbCJwOcXcWt/7LgTeCZzct0K2pNuB44HnbZ9R7hPwtjLPHIpVtftWyW50Xj2SlgK/D2B7Vs0x5wAXAe8Avg78ZbPzasT2cmA5wPjJ01+zqndERAyPdgrb48B0SUdQFINzgfm2r+kLkNQF3AEcA9wJnGG7p1Vi23uAWbX7JF0LnAOsAhbZXl/XPhf4KPBW27tqcl1UF3c+sAR4ErgVWFger89K4FJJK4DfBHrLOcCry1dfnlOBzwP/SlEULy/nGiMiYi/UsrDZ3i3pUuBBYBxwm+3NDUKXAav7RlD1yoI1H5ggaQdwi+1PNAhdQzEKbPYwxo3AeOChYkDGWtsfahD3HHCi7Reb5LkfOAPYBuyiGI018lPgXbafa9KOpGeBA4H9JL0bONX295vFR0TEyFGTOhQjaPzk6Z688PpOdyMiKq5q3xUpaZ3t2a3i8mRfRERUSgpbRERUSgpbRERUStZj64CZU7rorti974iIvUVGbBERUSkpbBERUSkpbBERUSmZY+uATT29TFu8qtPdiIgxomp/jzbSMmKLiIhKSWGLiIhKSWGLiIhKGdXCVi7QuVXSNkmLB/H+CZJWSdoiabOkzwyhLwdJekjSD8qfkwaZ5wpJ35e0UdLDkn5tsH2KiIihG7XCVq69dhPFatQzgPMkzaiL2U/SAS1Sfd72MRRrr/2upNObHK9VoVoMPGx7OsUipA0LbRt51gOzbR8HfBO4tkV8RESMoNEcsc0Bttl+plzPbAXFStu1JgGbJd0s6S31CWzvsr26/P1V4AmKBUIbuVfSSknzJDV6+vNMijXkKH++u0meGyQ9Iul8SW9o0KfVNevCre2nPxERMQpGs7BNAbbXbO8o9/2S7Z3A0cBqYKmk9ZIuk3RQfTJJE4F3UYy2GjkJuA44G3hK0qclHVXT/qvlwqJQLCL6q42S2H4fsAj4HYqie4OkNzc55geAB5q0RUTEKNjrHh6x/YrtFbZPpRhVnQI8L+mwvphyBPZVYJntZ5rkse01thcAJwAGtkg6q1Fs2d6sT+tsfxg4lmJh0sckXVEbI+l9wGzgc41ySLpYUrek7j27evu7BBERMQSjWdh6gMNrtqcCL0raUL7m9TVIOlTSR4D7KFbtng/srHnvcuAHtq8v48fV5PlkTZ79Jc0H7gZOAy4HHiqbd0qaXMZNBl4of3+wzHNLTZ7Xlf1bAfwRsAT4Sk37KcDVwDzbrzQ6edvLbc+2PXvchK62L1pERAzMaH7zyOPAdElHUBS5c4H5tq/pC5DURTHfdQxwJ3CG7Z7aJJL+DOgCPti3z/YeYFZd3LXAOcAqYJHt9XX9WQksBD5T/vzbMtdpdXmuAC4Fvg18wfa369qPB24G5tp+od2LERERI2PUCpvt3ZIuBR6kGIXdZntzg9BlwOry9uB/ImkqxchoC/CEJIAbbd9SHwusAZbYfrlJlz4DfF3SB4DngD9oErcRmGX7503aPwf8CvCNsj8/sj2vSWxERIywUf2uSNv3A/f3094LPNJP+w5AAzhWf+0/BU5uI88/tGg/pZ3+RETE6NjrHh6JiIgYihS2iIiolBS2iIiolKzH1gEzp3TRnfWVIiJGREZsERFRKSlsERFRKSlsERFRKZlj64BNPb1MW7yq092IiA54NvPrIy4jtoiIqJQUtoiIqJQUtoiIqJQUtoiIqJRRLWySbpP0gqQnh5DjCEnfk7RN0tck7TfIPOdI2izpF5JmD6E/50vaKGmTpO/0s7p2RESMgtEesf0VMLe/AEmTWuT4LPDnto8CfgZ8YJB5ngTeCzw6xP78C/BW2zOBT1EsghoRER0yqoXN9qPAv7UI65Z0l6S3q1zgrE+5/Xbgm+WuO4B3N8lzr6SVkuZJes2fNdh+yvbWNrq9SNJjki6RdGCDPN+x/bNycy3FyuAREdEhe+Mc268DX6VYtfr7kq6SdFjZ9kbg/9jeXW7vAKY0yXMScB1wNvCUpE9LOmqgnbF9FXABcCTF4qa3SzqxSfgHgAcaNUi6WFK3pO49u3oH2o2IiGjTXlfYbO+x/Xe23wv8HkVB+ZGkOQPMY9trbC8ATgAMbJF01iD6tNX2lcDRwMPAKknLamMkvY2isF3ZJMdy27Ntzx43oWugXYiIiDZ19JtHJB0O3Fdufsn2l8r9XcC5wIXAq8D7gY3AK8BESa8rR21TgR5J44B1ZZ6VtpeUefYH3lO+fyJwOfBQiz7dDhwPPG/7jHKfgLeVeeYAy4Bbat5zXLl9erkyd0REdEhHC5vt7cCs2n2SvgL8NvANYIHtH9S1r6a4vbgCWAj8re09DfJcC5wDrAIW2V7fZp8uqstzPrCE4mGTW4GF5fH62v87cDdwge2n2zlGRESMnFEtbJK+SjH3dbCkHcDHbd9aF/Z14MKaebR6VwIrJP0ZsJ6i2DSyBlhi++UmfXkPcANwCMWtxQ22T2sQ+hxwou0XmxxnCcXc3xfLZ1122x70nw9ERMTQyHan+7DPGT95uicvvL7T3YiIDsiXIA+epHXtDBz2uodHIiIihiKFLSIiKiXrsXXAzClddOd2RETEiMiILSIiKiWFLSIiKiWFLSIiKiVzbB2wqaeXaYtXdbobETGC8lh/52TEFhERlZLCFhERlZLCFhERlZLCFhERldKysEk6XNJqSd+XtFnS5YM5ULkq9lZJT0q6TdLrB5lnvKSvSdom6XuSpg0yzzGSvivpFUl/MpgcZZ5ZZZ7NkjZK+sPB5oqIiKFrZ8S2G/iI7RnAbwEfljSjPkjSpBZ57gKOAWYC+wMfbBTURp4PAD+zfRTw58BnB5nn34DLgM/3FySpS1J/12kXxfI6xwJzgeslTWxx7IiIGCEtC5vtH9t+ovz934GngCkNQrvLUdnby4U56/PcX65qbeAxikVCG7lB0iOSzpf0hgbtZwJ3lL9/Ezi50fGARZIek3SJpAMb9OcF248D/69JP/qcCGyV9Ily7bX6PE/3rRln+3ngBYqlcCIiogMGNMdW3vY7Hvheg+ZfB74KXAp8X9JVkg5rkOP1wAXA3zc6hu33AYuA3wE2S7pB0ptrQqYA28vY3UAvxXpo9XmuKo9zJPCEpNslndjemf6nPKsoFj7tBVZK+ntJ50jar8G5zQH2A37YoO1iSd2Suvfs6h1oNyIiok1tFzZJvwL8DfDHtn9e3257j+2/s/1e4PcoCsqPyn/sa30ReNT2t5sdy/Y62x8GjgW2AY9JuqLdvtbk2Wr7SuBo4GGKBUWXDSLPT2z/ue1ZwDXAJ4Hu2hhJk4E7gYts/6JBjuW2Z9uePW5C10C7EBERbWqrsJWjrL8B7rJ9d/lAyYby9aGauC5JlwArgenA+4GNNe0fp7hNd0XNvgfLPLfU7HudpHnACuCPKFap/krZ3AMc3hcHdAE/LUdkGyTdX5NHkt5OcetyCbAM+EKLc/1wzbkdVrN/hqTPAV8G/qnsV1/bgcAq4Grba1tczoiIGEEtv1KrnL+6FXjK9nUAtrcDs+rivkJxy+4bFA9T/KCu/YPAacDJtSMa26fVxV1BcTvz28AXGozsVgILge8CZwOPlPN2F9XlOZ+imD1Z9n+h7T2tztf2TcBNNXl+g2KU+Ysyz/G2X6pp3w+4B/iy7W+2yh8RESOrne+K/F2KuapNkjaU+66yfX9d3NeBC8t5r0a+BDwHfLd81uNu259sELcRmNXodmfpVuBOSdsonmw8t0ncc8CJtl9s1Cjpv1HcTjwQ+IWkPwZmNDju/6W4vfhUk+P8AcWt1zdKurDcd6HtDU3iIyJiBKkY7MRoGj95uicvvL7T3YiIEZQvQR5+ktbZnt0qLt88EhERlZLCFhERlZL12Dpg5pQuunObIiJiRGTEFhERlZLCFhERlZLCFhERlZI5tg7Y1NPLtMWrOt2NiBiCPM6/98qILSIiKiWFLSIiKiWFLSIiKiWFLSIiKmVMFjZJt0l6QdKTQ8hxl6Stkp4s871+kHl+T9ITknZLOnuw/YmIiOExJgsb8FfA3P4CJE1qkeMu4BhgJrA/8MFB5vkRcCHw1y3iIiJiFIzJwmb7UYola/rTXY7K3l6uKVef436XgMeAqU3y3CDpEUnnS3pDgzzP2t5IsV5bRER02JgsbG36deCrFIuWfl/SVbUrYvcpb0FeAPx9oyS23wcsAn4H2CzpBklvHrluR0TEUFS2sNneY/vvbL+XYiHQI4EfSZpTF/pF4NEGK3XX5lpn+8PAscA24LFype+2SbpYUrek7j27egd2MhER0bZKFDZJh0vaUL4+VLO/S9IlwEpgOvB+ihW6+9o/DhwCXFGz78Eyzy01+14naR6wAvgjYAnwlYH00fZy27Ntzx43oWtwJxoRES1V4iu1bG8HZtXuk/QV4LeBbwALbP+grv2DwGnAybZ/UZPrtLq4KyhuZ34b+EJ/I7uIiOi8MVnYJH0VOAk4WNIO4OO2b60L+zpwoe3dTdJ8CXgO+G75bMndtj/ZIG4jMMv2z5v05S3APcAk4F2SrrF97EDPKSIihseYLGy2z2sjZmWL9rbO3fY/tGh/nOZPVEZExCirxBxbREREnxS2iIiolBS2iIiolDE5xzbWzZzSRXcWKYyIGBEZsUVERKWksEVERKWksEVERKVkjq0DNvX0Mm3xqk53IyIG6NnMjY8JGbFFRESlpLBFRESlpLBFRESl7POFTdJcSVslbZO0eJA5lkraLuml4e5fREQMzD5d2CSNA24CTgdmAOdJmlEXs5+kA1qkug+oX8A0IiI6YJ8ubBTFaJvtZ2y/SrGQ6Jl1MZOAzZJuLpeoeQ3ba23/eIT7GhERbdjXC9sUYHvN9o5y3y/Z3gkcDawGlkpaL+kySQeNXjcjIqJd+3pha4vtV2yvsH0qxYjuFOB5SYe1m0PSxZK6JXXv2dU7Yn2NiNjX7euFrQc4vGZ7KvCipA3la15fg6RDJX2EYj5tHDAf2NnugWwvtz3b9uxxE7qGqfsREVFvX//mkceB6ZKOoChy5wLzbV/TFyCpC7gDOAa4EzjDdk8nOhsREa3t0yM227uBS4EHgaeAr9ve3CB0GfA/bC9tVNQkXStpBzBB0g5JnxjJfkdERHP7+ogN2/cD9/fT3gs80iLHR4GPDnPXIiJiEPbpEVtERFRPCltERFRKCltERFTKPj/H1gkzp3TRnXWdIiJGREZsERFRKSlsERFRKSlsERFRKZlj64BNPb1MW7yq092IiNKzmfOulIzYIiKiUlLYIiKiUlLYIiKiUtoqbJLmStoqaZukxYM5kKSlkrZLemkw76/J8zlJWyRtlHSPpImDzCNJy8pz2ijpNwaZ5x2S1knaVP58+2DyRETE8GhZ2CSNA24CTgdmAOdJmlEXs5+kA1qkug+Y08bxJrUIeQh4k+3jgKeBjw0yz+nA9PJ1MfAXg8zzE+BdtmcCCymWtomIiA5pZ8Q2B9hm+xnbrwIrKFaRrjUJ2CzpZklvaZTE9lrbP27jePdKWilpnqTXPLVp+1vlcjMAaykWB21kkaTHJF0i6cAG7WcCX3ZhLTBR0uQGcX8o6UlJH5F0SIP+rLf9fLm5Gdhf0viWZxkRESOincI2Bdhes72j3PdLtncCRwOrgaWS1ku6TNJBg+jTScB1wNnAU5I+LemoJrHvBx5o1GD7KuAC4EjgCUm3SzqxJqTleZV5vkQxupsAPCrpm+Wt2UbX7izgCduv9HeCERExcobt4RHbr9heYftUitHQKcDzkg4bYB7bXmN7AXACYGCLpLNq4yRdDewG7uon11bbV1IU3YeBVZKWDejEijzbbX+K4lbsbeXr3rr+HAt8FrikUQ5JF0vqltS9Z1fvQLsQERFtaucPtHuAw2u2pwIvStpQbi+xvRJA0qEUo6QFFCOg+cDOZonL+bt15eZK20vK/fsD76EYkU0ELqeYW+t734XAO4GTbbvcdztwPPC87TPKfQLeVuaZQ7ES9i39nFePpKXA7wPYnlVzzDnARcA7gK8Df1nTNhW4B1hg+4eNztX2cmA5wPjJ093smkRExNC0U9geB6ZLOoKiGJwLzLd9TV+ApC7gDuAYiocnzrDd0yqx7T3ArNp9kq4FzgFWAYtsr69rn0uxWvVbbe+qyXVRXdz5wBLgSeBWYGF5vD4rgUslrQB+E+gt5wCvLl99eU4FPg/8K0VRvLyca+xrn1j2dbHtf2p1zhERMbJaFjbbuyVdCjwIjANus725QegyYHXfCKpeWbDmAxMk7QBusf2JBqFrKEaBLzfp0o3AeOChYkDGWtsfahD3HHCi7Reb5LkfOAPYBuyiGI018lOKpx6fa9J+KXAUsETSknLfqbZfaBIfEREjSE3qUIyg8ZOne/LC6zvdjYgo5bsixwZJ62zPbhWXbx6JiIhKSWGLiIhKSWGLiIhKyXpsHTBzShfduacfETEiMmKLiIhKSWGLiIhKSWGLiIhKSWGLiIhKSWGLiIhKSWGLiIhKSWGLiIhKSWGLiIhKSWGLiIhKybf7d4Ckfwe2drofe7GDgZ90uhN7sVyf/uX69G8sX59fs31Iq6B8pVZnbG1n6YV9laTuXJ/mcn36l+vTv33h+uRWZEREVEoKW0REVEoKW2cs73QH9nK5Pv3L9elfrk//Kn998vBIRERUSkZsERFRKSlsw0DSXElbJW2TtLhB+3hJXyvbvydpWk3bx8r9WyWd1m7OsWSErs+zkjZJ2iCpe3TOZGQM9vpIeqOk1ZJeknRj3XtOKK/PNknLJGl0zmZ4jdC1WVPm3FC+Dh2dsxl+Q7g+75C0rvyMrJP09pr3jP3Pju28hvACxgE/BI4E9gP+GZhRF/M/gS+Vv58LfK38fUYZPx44oswzrp2cY+U1EtenbHsWOLjT59fh63MAcCLwIeDGuvc8BvwWIOAB4PROn+tedG3WALM7fX4dvj7HA4eVv78J6KnSZycjtqGbA2yz/YztV4EVwJl1MWcCd5S/fxM4ufy/oDOBFbZfsf0vwLYyXzs5x4qRuD5VMujrY/s/bP8j8HJtsKTJwIG217r4l+rLwLtH9CxGxrBfm4oZyvVZb/v5cv9mYP9ydFeJz04K29BNAbbXbO8o9zWMsb0b6AXe2M9728k5VozE9QEw8K3yNsrFI9Dv0TKU69Nfzh0tco4FI3Ft+txe3ob80zF5q60wXNfnLOAJ269Qkc9OvnkkxqoTbfeU8yMPSdpi+9FOdyrGhPPLz867Hz7qAAABk0lEQVR/Bf4GuIBiZLLPkXQs8Fng1E73ZThlxDZ0PcDhNdtTy30NYyS9DugCftrPe9vJOVaMxPXBdt/PF4B7GLu3KIdyffrLObVFzrFgJK5N7Wfn34G/Zh/97EiaSvHfzgLbP6yJH/OfnRS2oXscmC7pCEn7UUzQrqyLWQksLH8/G3ikvH+9Eji3vLd9BDCdYuK2nZxjxbBfH0kHlP+3jaQDKP5v88lROJeRMJTr05DtHwM/l/Rb5W22BcDfDn/XR9ywXxtJr5N0cPn764F3sg9+diRNBFYBi23/U19wZT47nX56pQov4AzgaYonlK4u930SmFf+/gbgGxQPPzwGHFnz3qvL922l5umjRjnH6mu4rw/FU2D/XL427+PX51ng34CXKOZDZpT7Z1P8g/1D4EbKL2MYa6/hvjYUT0uuAzaWn53/Tfmk7Vh8Dfb6AP8L+A9gQ83r0Kp8dvLNIxERUSm5FRkREZWSwhYREZWSwhYREZWSwhYREZWSwhYREZWSwhYREZWSwhYREZWSwhYREZXy/wEQBOOo+i66ggAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "solutions_with_probs = []\n",
    "for state in naive_distribution:\n",
    "    solutions_with_probs.append(['->'.join(map(str,state)), naive_distribution[state]/10000])\n",
    "solutions_with_probs = sorted(solutions_with_probs,key=lambda x: -x[1])\n",
    "\n",
    "plots.plot_state_histogram(solutions_with_probs[:10][::-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, still many (if not most) of the solutions don't make sense. In the next part you will see how to further improve this algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary and next steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we learned?\n",
    "- how to solve problem with constraints using QAOA\n",
    "- how to encode constraints and costs for TSP\n",
    "\n",
    "What's next?\n",
    "- improving the algorithm!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qenv",
   "language": "python",
   "name": "qenv"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}