ortools.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Capacitated Vehicle Routing Problem\n",
    "\n",
    "A short example of using Gurobi to model and solve the CRVP\n",
    "\n",
    "> Vehicle routing problems are inherently intractable: the length of time it takes to solve them grows exponentially with the size of the problem. For sufficiently large problems, it could take OR-tools (or any other routing software) years to find the optimal solution.. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import RandomState\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set()\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#or tools imports\n",
    "from ortools.constraint_solver import routing_enums_pb2\n",
    "from ortools.constraint_solver import pywrapcp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_PATIENTS = 15\n",
    "SEED = 101"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Patient Geographic Location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#random number generator for synthetic data.\n",
    "rand_gen = RandomState(seed=SEED)\n",
    "\n",
    "#n_patients = 1 so we include the depot\n",
    "x_coords = rand_gen.random(size=N_PATIENTS+1) * 200\n",
    "y_coords = rand_gen.random(size=N_PATIENTS+1) * 100\n",
    "\n",
    "coords = [[x, y] for x, y in zip(x_coords, y_coords)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_patient_coordinates(x_coords, y_coords):\n",
    "    '''\n",
    "    Plot a simple scatter of coordinates. \n",
    "    Depot is red, customers are blue\n",
    "    \n",
    "    Paramters:\n",
    "    x_coords: array-like\n",
    "        x coordinates for each patient. element 0 = depot.\n",
    "        \n",
    "    y_coords: array-like\n",
    "        y coordinates for each patient. element 0 = depot\n",
    "    '''\n",
    "    plt.plot(x_coords[0], y_coords[0], c='r', marker='^')\n",
    "    #plt.scatter(x_coords[1:], y_coords[1:], c='b')\n",
    "    \n",
    "    for i in range(1, len(x_coords[1:])):\n",
    "        x = x_coords[i]\n",
    "        y = y_coords[i]\n",
    "        plt.plot(x, y, 'b*')\n",
    "        plt.text(x * (1 + 0.01), y * (1 + 0.01) , i, fontsize=12)\n",
    "    \n",
    "    plt.xlim((0, 220))\n",
    "    plt.ylim((0, 120))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD/CAYAAAAQaHZxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de1hUdeIG8Hfm4HBJcAAB8bJimDpaSeE+ZDcVbadsBMt2p0hrzaxsybYyo7UAxcqxi/akrlot3cxtzcoYtTEj3c3H1Lz05M56B0MZuQwgF3VwLr8/+DmJXIe5nGHO+3menoc5B855/XKel9N35pwjczgcDhARUUCTix2AiIi8j2VPRCQBLHsiIglg2RMRSQDLnohIAlj2REQS0GHZ63Q6pKamYujQoThy5AgAoLq6GjNnzoRarcakSZOQmZmJqqoq588UFRVBq9VCrVZDq9WiuLjYa/8AIiLqWIdlP378eKxZswb9+vVzLpPJZHj00UdhMBhQUFCAAQMG4I033nCuz8nJQUZGBgwGAzIyMpCdne2d9ERE1Ckdlv2oUaMQHx/fbJlSqURKSorzdVJSEkpLSwEAZrMZRqMRGo0GAKDRaGA0Gpud+RMRkW+5PWdvt9uxdu1apKamAgBMJhPi4uIgCAIAQBAExMbGwmQyubsrIiLqIrfLPi8vD2FhYZg6daon8hARkRcEufPDOp0OJ0+exMqVKyGXN/3diI+PR1lZGWw2GwRBgM1mQ3l5eYupoM6orm6A3e5ft+6Jju4Js7le7BgtMJdrmMs1zOUasXLJ5TJERl7V6roul/2SJUtw8OBBrF69GgqFwrk8OjoaKpUKer0e6enp0Ov1UKlUiIqKcnkfdrvD78oegF9mApjLVczlGuZyjb/lknV018uFCxdiy5YtqKysRGRkJJRKJZYuXQqNRoOEhASEhIQAAPr374/ly5cDAI4fP46srCzU1tYiIiICOp0OV199tcvhzOZ6vxuwmJhwVFTUiR2jBeZyDXO5hrlcI1YuuVyG6Oiera7rsOzFxLLvPOZyDXO5hrlc449lzytoiYgkgGVPRCQBLHsiIglg2RMRSQDLnohIAlj2REQSwLInIpIAlj0RkQSw7ImIJIBlT0QkASx7IiIJYNkTEUkAy56ISAJY9kREEsCyJyKSAJY9EZEEsOyJiCSAZU9EJAEseyIiCWDZExFJAMueiEgCWPZERBLAsvcTxcVFmD37CajVY6DVTsb27d+LHYmIAgjL3g9YrVZkZT2Hm2++FZs2FWLu3HnIy3sZv/56UuxoRBQgWPZ+4Ndfi2E2V0CrfRCCICA5+fe47rqRMBg2iR2NiAIEy94POBytLXPgxInjvg9DRAGJZe8HBg5MgFIZhU8//QhWqxW7d/+IAwf2wWK5IHY0IgoQHZa9TqdDamoqhg4diiNHjjiXFxUVQavVQq1WQ6vVori4uFPrqKWgoCC89tob2LlzB9LS1PjnPz9BauodiImJFTsaEQWIDst+/PjxWLNmDfr169dseU5ODjIyMmAwGJCRkYHs7OxOraPWDR58DZYtW41Nm77DW28tQ2npaahUI8SORUQBosOyHzVqFOLj45stM5vNMBqN0Gg0AACNRgOj0Yiqqqp211FzNfUWLFqzF2frLTh27CgsFgsuXLiATz/9GGZzJSZOnCR2RCIKEEFd+SGTyYS4uDgIggAAEAQBsbGxMJlMcDgcba6LioryXPIAULCjCEdLzuLrHUWoPWZAQcFXsNmsuP76G7BkyXIoFAqxIxJRgOhS2ftKdHRPsSO0KiYm3K2fv/eFAly02p2vv99fCuA6DJwwEl/oun42724ub2Eu1zCXa5irc7pU9vHx8SgrK4PNZoMgCLDZbCgvL0d8fDwcDkeb61xlNtfDbm/lc4kiiokJR0VFnVvb0D0xGp8VHsP+IxVotNqhCJLjxiEx0KYO7vK2PZHLG5jLNczlGuZqTi6XtXmS3KWPXkZHR0OlUkGv1wMA9Ho9VCoVoqKi2l1HTZQ9gxGqEHDRakcPQY6LVjtCgwX06hksdjQiClAdntkvXLgQW7ZsQWVlJaZPnw6lUomNGzciNzcXWVlZWLFiBSIiIqDT6Zw/0946alJ7rhFjb+iHMUl9sf1AKc42WMSOREQBTOZwtHb9pn8I1Gkcb2Au1zCXa5jLNQEzjUNERN0Ly56IqIvWr/8MM2ZMw7hxo/HKK7mtfs8//rEat946Cnv27PJtuCv49UcviYj8We/eMXj44RnYvXsnLJaW77udPn0K27Z9h+jo3iKka45n9kREXTRmTCpuv30sIiJ6tbr+rbcWY9asp9CjRw8fJ2uJZU9E5AWFhVvRo0cQRo++VewoAFj2REQe19DQgNWrl2P27OfEjuLEsici8rB33nkHavVE9O3br+Nv9hG+QUtE5GE7d+6EyWTCl19+DgCoqalGdvaLePDBhzB16p9FycSyJyLqhJp6C1ZuOIhZ6dc6b21itVphs9lgt9tht9tgsVggCAI++OADlJXVOH925syHkZn5DG666Wax4rPsiYg64/Jbkk9TDwMAfPjh+8jPf9f5PQbDZkyfPhNZWXNgtf5Wr3K5HOHh4QgLC/N57ktY9kRE7Xj89W24aGt+S/Lv95eihyDHqucfx4wZj3e4jc8/L/BmxE7hG7RERO3QzRqNlOFxUAQ11aUiSI6bhsdh8azRIidzDcueiKgdgXJLck7jEBF1IBBuSc6yJyLqQOa91zu/nqYeKmKSruM0DhGRBLDsiYgkgGVPRCQBLHsiIglg2RMRSQDLnohIAlj2REQSwLInIpIAlj0RkQSw7ImIJIBlT0QkASx7IiIJYNkTEUmA22X//fffY/LkyUhPT8ekSZOwZcsWAEBRURG0Wi3UajW0Wi2Ki4vd3RUREXWRW7c4djgcmDt3LtasWYMhQ4bg0KFDeOCBBzBhwgTk5OQgIyMD6enp2LBhA7Kzs/HRRx95KjcREbnA7TN7uVyOuro6AEBdXR1iY2NRXV0No9EIjUYDANBoNDAajaiqqnJ3d0RE1AVundnLZDIsXboUTz75JMLCwtDQ0IBVq1bBZDIhLi4OgiAAAARBQGxsLEwmE6Kiojq9/ejonu7E85qYmHCxI7SKuVzDXK5hLtf4Wy63yt5qtWLVqlVYsWIFkpOTsXfvXjzzzDNYvHixR8KZzfWw2x0e2ZanxMSEo6KiTuwYLTCXa5jLNczlGrFyyeWyNk+S3Sr7//3vfygvL0dycjIAIDk5GaGhoQgODkZZWRlsNhsEQYDNZkN5eTni4+Pd2V3AWb/+M2zapMeJE8cwYYIa8+blAgAOHvwF7733dxw+fAiCIEdSUjL++tfn0bt3b3EDE1G35dacfZ8+fXDmzBmcOHECAHD8+HFUVlZi4MCBUKlU0Ov1AAC9Xg+VSuXSFE5nNDY24rXXFmDKFA3uuON2TJ+egZ07d3h0H97Uu3cMHn54Bu6+O63Z8rq6WqSl3YvPP/8an3+uR1hYGF59db5IKYkoELh1Zh8TE4Pc3Fw8/fTTkMlkAIDXXnsNSqUSubm5yMrKwooVKxAREQGdTueRwJez2WyIjY3DsmWrERfXBzt37kB29ov46KN/Ij6+r8f352ljxqQCAA4dMqKioty5fPToW5p935QpWmRmPubTbEQUWNwqewBIS0tDWlpai+WJiYlYt26du5tvV2hoKGbMeNz5+pZbbkPfvn1x+PD/ukXZd9bPP+/DoEFXix2DiLqxgLqCtqrKjJKSXzFoUKLYUTzm2LGjyM9/D3/5y9NiRyGibixgyt5qtWL+/Jdx5513Y+DABLHjeMSpUyWYM2c2nn76OYwceYPYcYioGwuIsrfb7cjLexk9egTh2WdfEDuOR5w5Y8Jf//ok/vznGbjzzrvFjkNE3Vy3L3uHw4FFi/JQVVWFV15ZjKAgt9+G8KqaegsWrdmLs/UWWK1WWCwW2O122O02WCxNyyoqyjF79hO4994/YvLk+8SOTEQBwL+bsQ019Ras3HAQs9Kvxeq/v4ni4iIsXboCwcEhYkfrUMGOIhwtOYuvdxSh8dR25Oe/61xnMGzG9OkzIZPJUFp6Gvn57zZb/+23/xEjMhEFgG5Z9pcKc83GPdiw4QsoFAqkp6ud659//m/4wx/uEjFhS4+/vg0XbXbn6+/3lwK4BiPS38Cq58e2+P5HHuFHLcnzSkp+xcMP34+xY8cjOztP7DjkQ92q7K8szJ+KGjFEsxg9BHmrhelPdLNG47PCY9h/pAKNVjsUQXLcOCQG2tTBYkcjCXnrLR2GDRsudgwSQbeas9fNGo2U4XFQBDXFVgTJcdPwOCyeNVrkZB1T9gxGqELARasdPQQ5LlrtCA0W0KtnsNjRSCK2bjWgZ89wJCf/XuwoJIJuVfbdvTBrzzVi7A39MO+hZIy9oR/ONjSKHYkkoqGhHu+9twqZmX8VOwqJpFtN4wC/FeaYpL7YfqAUZxssYkfqtMx7r3d+PU09VMQkJDXvvrsSGk0a4uL6iB2FRNLtyp6FSeSao0cP46efdiM/f43YUUhE3a7sicg1+/fvxZkzpZgypenJcefPn4PNZkdx8Qn84x/8AyAVLHuiAHP5dSi9egYjLe1ejB//B+f6tWs/wZkzpXjuuRdFTEm+1q3eoCWijl1+4R4AhISEIDq6t/O/0NBQKBTBiIyMFDkp+RLP7IkCRGsX7n2/v7TFdSiX3xacpINn9kQBojtfh0Lex7InChDd/ToU8i5O4xAFkO58HQp5F8ueKIDwOhRqC6dxiIgkgGVP5AXr13+GGTOmYdy40XjllVyx4xBxGofIG3r3jsHDD8/A7t07YbFw3pzEx7In8oIxY1IBAIcOGVFRUS5yGiJO4xARSQLLnohIAlj25DFtvSl58eJFvPTSXNx33yTceuso7Nv3k3ghiSSKZU8ec+lNybvvTmux7vrrk/Dyy3mIjo4WIRkRuf0GrcViwauvvoqdO3ciODgYSUlJyMvLQ1FREbKyslBTUwOlUgmdToeEhAQPRCZ/1dabkj169MCf/pQBAJDLBVGyeduVtxW2Wq2w2Wyw2+2w222wWCwQBAFBQfxMBInD7SPv9ddfR3BwMAwGA2QyGSorKwEAOTk5yMjIQHp6OjZs2IDs7Gx89NFHbgcm8keX31Z4mnoYPvzwfeTnv+tcbzBsxvTpM3nHSRKNW2Xf0NCAr776Ctu3b4dMJgMA9O7dG2azGUajEfn5+QAAjUaDvLw8VFVVISoqyv3URH6i7dsKD8UPP/C9CfIfbpV9SUkJlEolli1bhl27duGqq67C008/jZCQEMTFxUEQmv6XXRAExMbGwmQysewpoOhmjcZnhcew/0gFGq12KILkuHFIDLSpg8WORtSMW2VvtVpRUlKC4cOH44UXXsDPP/+MJ554Am+//bZHwkVH9/TIdjwtJiZc7AitEiNXVe0FLP74J7wwbRQiI0IAAGFhCoSE9HDmuTyXXC6DUhnmF2PoiQwxMeGI7HUaF2129AiS46LNjkhlKAYP6i1qLm/oaq5PPvkEX3zxBY4cOQKNRoNFixYBAI4dO4a5c+eipKQEADBixAi89NJLGDzYtT+UgTZe3uJW2fft2xdBQUHQaJoeZDxy5EhERkYiJCQEZWVlsNlsEAQBNpsN5eXliI+Pd2n7ZnM97HaHOxE9LiYmHBUVdWLHaEGsXB8bDsF4woz8r3/BA+MHw2azob7+As6du4BTpyrRp48S1dXn0djYCIfDAbvdgcrKWpw6VQmFQuGc/vM1T45XubkeY5N+u61wWWV9l7cdiMdXSEg4HnxwOnbv3okLFyzO7cjlYcjNfQ19+sTDbrfjiy/WYfbsp/Hhh//0SS5vEiuXXC5r8yTZrbKPiopCSkoKduzYgVtvvRVFRUUwm81ISEiASqWCXq9Heno69Ho9VCoVp3ACSGtz1ev++QGqjm51LjMYNiMzMxP33/9nZGRMwZkzJgDAs89mAgDWrfsa8fF9fRvcC3hb4fa19Smt8PBwhIc3nf06HA7I5XKcOlUiSkYpcPvTOPPnz8ff/vY36HQ6BAUFYfHixYiIiEBubi6ysrKwYsUKREREQKfTeSIv+YnW5qo190yDNnV+sycjXTrD+fzzAhHTkj+7886xOH/+POx2Oz+t5EVul/2AAQPw8ccft1iemJiIdevWubt58lN8BN5v1q//DJs26XHixDFMmKDGvHm5znU//bQbb72lQ1nZGQwffi3mzctFnz6uTWcGum++2Ybz589j82Y9x8aLeAUtddmlR+DNeygZY2/oh7MNjWJHEkVbVw7X1NRg3rzn8eijs7BpUyGGDRuO7OwXRUrp30JDQzF58hQsXJiD6uoq5/LMzMeQmnoz7rjjNtxxx2144IF7RUzZvfFyPuoyzlU3aWtOevv2QgwalIjU1AkAgEceeQx33z0BJ08WY+DABDGi+jW73Y4LFy6goqIckZG/vb/3zDNzMWnSZBGTBQae2RN5SVHRCQwefI3zdWhoKPr164eiouMipvKumnoLFq3Zi7P1vz2wxWq1wmKxNLt1hNVqxZ49P+LIkUOw2WxoaKjHsmVLEB4ejoEDB4n4LwhcPLMn8pLz589BqYxstqxnz544d+6cSIm878rbRgBo89YRgwYlYsmS11FRUY7g4GAMGzYcb775DoKDm7/vs2rVMqxc+Q5+97uBmDnzSdx44yif/psCBcueyEtCQ8PQ0NDQbFlDQwPCwsJESuQ9bd82Qo5Vzz/e5qdsLk1xtWXWrNkYNGgQgoJ64LvvtuCFF57FBx98in79+ns0vxRwGofISwYNuhrHjx9xvj5//jxOnz6FQYMSRUzlHbpZo5EyPA6KoKZKUQTJcdPwOCyeNdqt7Y4YcS3Cwq6CQqHAXXdpcN11I7Fz5w+eiCw5LHsiF105L93WnPTtt4/DiRPHsW3bd7BYLMjPfxeJidcE5Juznvgobmvz/VeSyWRw+NdF9d0Gy57IRZfPSwNNc9Ljx9+CTz75AAbDZowffws+/PB9REZGYuHCxVi9egXuuisVRuNBzJ//qsjpvcfdj+JeOa51dXXYtWun84/nli2b8fPP+5CScpM34gc8mcPhv38neW+czmMu13Ql15Xz0pc0zUuPFS2XL3gzV1vjKrOeg6z4Xzh5shiCIMfvfpeAmTOfwO9//1vZS3G82uO1e+MQSQlvZ+wd7Y1rr54aseMFDE7jEHUSbxHhHRxX3+CZPZELLs1LX7qd8dmGtt9MpM7juHofy57IBbxFhHdwXL2P0zhERBLAsicikgCWPRGRBLDsiYgkgGVPRCQBLHsiIglg2RMRSQDLnohIAlj2REQSwLInIpIAlj0RkQSw7ImIJIBlT0QkASx7IiIJ8FjZL1u2DEOHDsWRI0cAAEVFRdBqtVCr1dBqtSguLvbUroiIyEUeKfv//ve/OHDgAPr27etclpOTg4yMDBgMBmRkZCA7O9sTuyIioi5wu+wbGxuxYMEC5OTkQCaTAQDMZjOMRiM0mqbnR2o0GhiNRlRVVbm7OyIi6gK3y/7tt99GWloaBgwY4FxmMpkQFxcHQRAAAIIgIDY2FiaTyd3dERFRF7j1WML9+/fjl19+wZw5czyVp5no6J5e2a67YmLCxY7QKuZyDXO5hrlc42+53Cr7PXv24MSJExg/fjwA4MyZM5gxYwZefPFFlJWVwWazQRAE2Gw2lJeXIz4+3qXtm831sNsd7kT0uJiYcFRU1IkdowXmcg1zuYa5XCNWLrlc1uZJslvTOI899hh++OEHFBYWorCwEH369MH777+PiRMnQqVSQa/XAwD0ej1UKhWioqLc2R0REXWRW2f27cnNzUVWVhZWrFiBiIgI6HQ6b+2KiIg64NGyLywsdH6dmJiIdevWeXLzRETURbyClohIAlj2REQSwLInIpIAlj0RkQSw7ImIJIBlT0QkASx7IiIJYNkTEUmA166gpdYtWPAy9u7djfPnLyAqKhoPPvgQJk2aLHYsIgpwLHsfmzr1z8jKehkKhQInTxbjqacexzXXDMWwYSqxoxFRAOM0jo9dfXUiFAoFAEAma/rv9OlTIqciokDHM3sRvPHGImzeXACLxYIhQ4Zi9OhbxI5ERAGOZS+COXOy8Mwzz+PgwV+wf/9PzjN9IiJv4TSOSARBwMiRSaioKMeXX34udhwiCnAsex+oqbdg0Zq9OFtvabHOZrNxzp6IvI5l7wMFO4pwtOQsPttyAFu3GnDu3DnYbDbs2rUTW7cakJw8SuyIRBTgOGfvRY+/vg0XbXbn6x9+KcO6vR/DsjAPIQoBffr0wezZz+G228aKF5KIJIFl70W6WaPxWeEx7D9SgUarHWFXReCPj+ZCmzoYvXoGix2PiCSE0zhepOwZjFCFgItWO3oIcly02hEaLLDoicjneGbvZbXnGjH2hn4Yk9QX2w+U4mxDyzdpiYi8jWXvZZn3Xu/8epp6qIhJiEjKOI1DRCQBLHsiIglg2RMRSQDn7ImoGZOpFG++uQgHD/4ChUKBsWNTkZeXK3YschPP7ImomTffXITIyChs2PAN8vPX4MCBffj000/FjkVuYtkTUTMmUylSUycgODgY0dG9kZJyM44dOyZ2LHKTW2VfXV2NmTNnQq1WY9KkScjMzERVVRUAoKioCFqtFmq1GlqtFsXFxZ7IS0Re9sc/3o+tW7fgwoULqKgox48/7sBtt90mdixyk1tlL5PJ8Oijj8JgMKCgoAADBgzAG2+8AQDIyclBRkYGDAYDMjIykJ2d7ZHARORdSUnJKCo6AbV6DO65ZyKGDRuOCRMmiB2L3ORW2SuVSqSkpDhfJyUlobS0FGazGUajERqNBgCg0WhgNBqdZ/1E5J/sdjuefTYTY8aMw7ff/gcbN25FXV0tXn/9dbGjkZs8Nmdvt9uxdu1apKamwmQyIS4uDoIgAGh6UEdsbCxMJpOndkdEHnDlsxZqa2tRXl6GKVO0UCgU6NVLiYkT0/Dvf/9b5KTkLo999DIvLw9hYWGYOnUqjEajR7YZHd3TI9vxtJiYcLEjtIq5XMNcwLp/n8DRU2exZd9pPDllJGJiwtG/f398+20BHnnkEZw7dw6Fhd9g6NChHC8X+VsumcPhcLi7EZ1Oh8OHD2PlypVQKBQwm81Qq9XYtWsXBEGAzWZDSkoKtmzZgqioqE5v12yuh93udjyPiokJR0VFndgxWmAu10g915XPWrikhyDHnMnxePvtN3Hs2FEIghw33DAKr76aB7vd/56VLPXf45XkclmbJ8lun9kvWbIEBw8exOrVq50Pzo6OjoZKpYJer0d6ejr0ej1UKpVLRU9E3nPlsxYUQXLcOCTG+ayFZctWN/v+6Gj/LNXOuOOO5p8kslgsuOee+/DMM3NFSiQOt8r+6NGjWLlyJRISEnD//fcDAPr374/ly5cjNzcXWVlZWLFiBSIiIqDT6TwSmIjcJ6VnLXz77X+cX58/fx5paX/AuHHS+3SRW2V/zTXX4PDhw62uS0xMxLp169zZPBF5kRSftbBt23dQKqMwcuQNYkfxOd4bh0iipPishc2b9bjzzomQyWRiR/E53i6BiCThzJkzOHBgH+66SyN2FFGw7IlIEr75Ro/rr09C3779xI4iCpY9EQWUKy8Uu+SbbzbhzjvvFimV+Fj2RBRQCnYU4WjJWXy9o8i57JdffkZlZTlSU6X3KZxL+AYtEQWEKy8U+35/Kb7fX4oeghyJjh0YM2YcwsKuEjGhuFj2RBQQ2r9QbKzY8UTHaRwiCghSulCsK3hmT0QBQ4oXinUWy56IAoYULxTrLE7jEBFJAMueiEgCWPZERBLAsicCYK2pQcniV2E9WyN2FCKvYNlTm7ZuNeDBB+/DhAm34k9/SsfPP+8XO5LXmPUbcP7oUZgLvhY7CpFX8NM41Ko9e37E3//+DubPfw3Dh4+A2VwpdiSvsdbUoHbHD4DDgdod/0H0pDQE9VKKHYvIo3hmT616//3VmD79UVx77XWQy+WIiYlFTEys2LG8wqzfAIe96TJ7h93Os3sKSCx7asFms+HQISOqq2ug1U7GPfdMxFtv6WCxXBA7msc5z+pttqYFNhtqd/yHc/cUcFj21EJ1dRWsViu2bfsOy5e/h/z8T3H06GF88MH7YkfzuMvP6i/h2T0FIpY9taBQNN1L5L77tOjduzeUSiW02gfx4487RE7meReOH//trP4Smw0Xjh8TJxCRl/ANWgLQ9MCHlRsOYlb6tegVEYHY2DixI/nEwJwFYkcg8gme2ROAlg98mDhxEtav/xeqq6tQW1uLf/1rLW6++TaRUxJRV/HMXuLaeuBDkOwaDBtWgwceuBcKRTBSUyfgoYceETEpEbmDZS9x7T/wYTzmzMkSOyIReQCncSSOD3wgkgae2RMf+EAkASx74gMfiCSA0zhERBLg1bIvKiqCVquFWq2GVqtFcXGxN3dHRERt8GrZ5+TkICMjAwaDARkZGcjOzvbm7oiIqA1em7M3m80wGo3Iz88HAGg0GuTl5aGqqgpRUVGd2oZcLvNWPLcwl2uYyzXM5Rrm6tw+vVb2JpMJcXFxEAQBACAIAmJjY2EymTpd9pGRV3krnluio3uKHaFVzOUa5nINc7nG33LxDVoiIgnwWtnHx8ejrKwMtv+/o6DNZkN5eTni4+O9tUsiImqD18o+OjoaKpUKer0eAKDX66FSqTo9hUNERJ4jczgcDm9t/Pjx48jKykJtbS0iIiKg0+lw9dVXe2t3RETUBq+WPRER+Qe+QUtEJAEseyIiCWDZExFJAMueiEgC/O4Wx0VFRcjKykJNTQ2USiV0Oh0SEhJ8nqO6uhpz587Fr7/+CoVCgYEDB2LBggWIiopCamoqFAoFgoObHvAxZ84c3Hab757P2tb+xRy7U6dO4S9/+YvzdV1dHerr67F7926fj5dOp4PBYMDp06dRUFCAIUOGAGj/2PLF2LWWq73jDGj7d+3tXB3tW6zxau846yizJ7T3+xL7+OqQw89MmzbN8dVXXzkcDofjq6++ckybNk2UHNXV1Y4ff/zR+XrRokWOF1980eFwOBzjxo1zHD58WJRc7e3fX8bO4XA4Fi5c6Jg/f77D4fD9eO3Zs8dRWlraYr/tjaXNASgAAAK3SURBVI8vxq61XO0dZw6Hb8aurfFqb99ijdeVLj/OOsrsCe39vsQ+vjriV9M4l26eptFoADTdPM1oNKKqqsrnWZRKJVJSUpyvk5KSUFpa6vMcneVPY9fY2IiCggJMmTLF5/sGgFGjRrW4Uru98fHV2LWWyx+Os9ZytUfM8bqcGMdZW78vfzi+OuJX0zieuHmaN9jtdqxduxapqanOZXPmzIHD4UBycjKeffZZRERE+DTTlfv3p7ErLCxEXFwcRowY0WZeX49Xe+PjcDj8YuxaO84AcceutX37y7HW2nHWVmZvuPz31R2OL786s/dXeXl5CAsLw9SpUwEAa9aswddff43169fD4XBgwYIFPs0j9v47sn79+mZnW/6e119ceZwB4o6dv//erjzOAN9mbu335c/8quz98eZpOp0OJ0+exNKlSyGXy505AUChUCAjIwP79u3zaabW9u8vY1dWVoY9e/Zg0qRJ7eb1tfbGxx/GrrXj7FJuQJyxa2vf/jBerR1n7WX2tCt/X/5+fAF+Vvb+dvO0JUuW4ODBg1i+fDkUCgUA4Ny5c6irqwMAOBwObNq0CSqVymeZ2tq/v4zdl19+iTFjxiAyMrLdvL7W3viIPXatHWeAuGPX3r7FHi+g5XHWUWZPau335c/H1yV+d28cf7l52tGjR6HRaJCQkICQkBAAQP/+/ZGVlYWnnnoKNpsNdrsdiYmJeOmllxAbG+uTXCUlJW3u3x/GTq1WY968ebj99ts7zOstCxcuxJYtW1BZWYnIyEgolUps3Lix3fHxxdi1lmvp0qWtHmfLly/32di1lmvlypXt7lus8dq4cSOAlscZ4Jtjra1eWL58uejHV0f8ruyJiMjz/Goah4iIvINlT0QkASx7IiIJYNkTEUkAy56ISAJY9kREEsCyJyKSAJY9EZEE/B/1mPX7v/5F0QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#visualise\n",
    "plot_patient_coordinates(x_coords, y_coords)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Approaches to 'solving' the CVRP\n",
    "\n",
    "* Formulation as a Mixed Integer Programme (MIP) and use of a optimisation solver\n",
    "* **Heuristic solution**.  E.g. using Iterated Local Search."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution via meta-heuristic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Iterated Local Search - multiple runs of Ordinary Decent 2Opt\n",
    "objective = SimpleTSPObjective(matrix)\n",
    "local_search = OrdinaryDecent2OptNew(objective, tour)\n",
    "runner = IteratedLocalSearch(objective, local_search, maximisation=False)\n",
    "n = 20\n",
    "print(\"\\nRunning Iterated Local Search using Ordinary Decent 2-Opt ({} runs)...\"\\\n",
    "      .format(n))\n",
    "runner.run(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_data_model(coords):\n",
    "    data = {}\n",
    "    data['distance_matrix'] = gen_matrix(np.array(coords))\n",
    "    data['num_vehicles'] = round(N_PATIENTS / 4) \n",
    "    data['depot'] = 0\n",
    "    data['vehicle_capacities'] = [4 for i in range(data['num_vehicles'])]\n",
    "    data['demands'] = [1 for i in range(N_PATIENTS)]\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Vehicles Routing Problem (VRP).\"\"\"\n",
    "\n",
    "from __future__ import print_function\n",
    "from ortools.constraint_solver import routing_enums_pb2\n",
    "from ortools.constraint_solver import pywrapcp\n",
    "\n",
    "def print_solution(data, manager, routing, solution):\n",
    "    \"\"\"Prints solution on console.\"\"\"\n",
    "    max_route_distance = 0\n",
    "    total_distance = 0\n",
    "    routes = []\n",
    "    for vehicle_id in range(data['num_vehicles']):\n",
    "        index = routing.Start(vehicle_id)\n",
    "        plan_output = 'Route for vehicle {}:\\n'.format(vehicle_id)\n",
    "        route_distance = 0\n",
    "        route = []\n",
    "        while not routing.IsEnd(index):\n",
    "            plan_output += ' {} -> '.format(manager.IndexToNode(index))\n",
    "            route.append(manager.IndexToNode(index))\n",
    "            previous_index = index\n",
    "            index = solution.Value(routing.NextVar(index))\n",
    "            route_distance += routing.GetArcCostForVehicle(\n",
    "                previous_index, index, vehicle_id)\n",
    "            total_distance += route_distance\n",
    "        plan_output += '{}\\n'.format(manager.IndexToNode(index))\n",
    "        plan_output += 'Distance of the route: {}m\\n'.format(route_distance)\n",
    "        print(plan_output)\n",
    "        max_route_distance = max(route_distance, max_route_distance)\n",
    "        routes.append(route)\n",
    "    print('Maximum of the route distances: {}m'.format(max_route_distance))\n",
    "    print('total distance: {}m'.format(total_distance))\n",
    "    return routes\n",
    "\n",
    "\n",
    "def main():\n",
    "    \"\"\"Solve the CVRP problem.\"\"\"\n",
    "    # Instantiate the data problem.\n",
    "    data = create_data_model(coords)\n",
    "\n",
    "    # Create the routing index manager.\n",
    "    manager = pywrapcp.RoutingIndexManager(len(data['distance_matrix']),\n",
    "                                           data['num_vehicles'], data['depot'])\n",
    "\n",
    "    # Create Routing Model.\n",
    "    routing = pywrapcp.RoutingModel(manager)\n",
    "\n",
    "\n",
    "    # Create and register a transit callback.\n",
    "    def distance_callback(from_index, to_index):\n",
    "        \"\"\"Returns the distance between the two nodes.\"\"\"\n",
    "        # Convert from routing variable Index to distance matrix NodeIndex.\n",
    "        from_node = manager.IndexToNode(from_index)\n",
    "        to_node = manager.IndexToNode(to_index)\n",
    "        return data['distance_matrix'][from_node][to_node]\n",
    "\n",
    "    transit_callback_index = routing.RegisterTransitCallback(distance_callback)\n",
    "\n",
    "    # Define cost of each arc.\n",
    "    routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)\n",
    "\n",
    "\n",
    "    # Add Capacity constraint.\n",
    "    def demand_callback(from_index):\n",
    "        \"\"\"Returns the demand of the node.\"\"\"\n",
    "        # Convert from routing variable Index to demands NodeIndex.\n",
    "        from_node = manager.IndexToNode(from_index)\n",
    "        return data['demands'][from_node]\n",
    "\n",
    "    demand_callback_index = routing.RegisterUnaryTransitCallback(\n",
    "        demand_callback)\n",
    "    routing.AddDimensionWithVehicleCapacity(\n",
    "        demand_callback_index,\n",
    "        0,  # null capacity slack\n",
    "        data['vehicle_capacities'],  # vehicle maximum capacities\n",
    "        True,  # start cumul to zero\n",
    "        'Capacity')\n",
    "\n",
    "    # Setting first solution heuristic.\n",
    "    search_parameters = pywrapcp.DefaultRoutingSearchParameters()\n",
    "    search_parameters.first_solution_strategy = (\n",
    "        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)\n",
    "\n",
    "    #time limit\n",
    "    search_parameters.time_limit.seconds = 30\n",
    "    \n",
    "    # Solve the problem.\n",
    "    solution = routing.SolveWithParameters(search_parameters)\n",
    "\n",
    "\n",
    "    # Print solution on console.\n",
    "    if solution:\n",
    "        routes = print_solution(data, manager, routing, solution)\n",
    "        \n",
    "    print(routes)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Route for vehicle 0:\n",
      " 0 ->  2 ->  1 -> 0\n",
      "Distance of the route: 0m\n",
      "\n",
      "Route for vehicle 1:\n",
      " 0 ->  6 ->  5 ->  4 ->  3 -> 0\n",
      "Distance of the route: 0m\n",
      "\n",
      "Route for vehicle 2:\n",
      " 0 ->  10 ->  9 ->  8 ->  7 -> 0\n",
      "Distance of the route: 0m\n",
      "\n",
      "Route for vehicle 3:\n",
      " 0 ->  15 ->  14 ->  13 ->  12 ->  11 -> 0\n",
      "Distance of the route: 0m\n",
      "\n",
      "Maximum of the route distances: 0m\n",
      "total distance: 0m\n",
      "[[0, 2, 1], [0, 6, 5, 4, 3], [0, 10, 9, 8, 7], [0, 15, 14, 13, 12, 11]]\n"
     ]
    }
   ],
   "source": [
    "main()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_solution(x, instance):\n",
    "    active_arcs = [a for a in instance['A'] if x[a].x == 1.00]\n",
    "    xc = instance['x_coords']\n",
    "    yc = instance['y_coords']\n",
    "    \n",
    "    for i, j in active_arcs:\n",
    "        plt.plot([xc[i], xc[j]], [yc[i], yc[j]], c='green')\n",
    "        plt.plot(xc[0], yc[0], c='r', marker='s')\n",
    "        plt.scatter(xc[1:], yc[1:], c='b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixed Integer Programme (MIP) formulation\n",
    "\n",
    "* The advantage of a mathematical forumulation of the problem is that we can use a modern solver like **Gurobi**.\n",
    "* Gurobi is extremely powerful.  The downside is that a commerical tool and the source code is closed (but free to use in research)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$n$ is the number of patients\n",
    "\n",
    "$N$ is the set of patients $\\{1, 2, ..., n\\}$\n",
    "\n",
    "$V$ is the set of verticies (nodes) to visit. $V=\\{0\\}\\cup N$\n",
    "\n",
    "$A$ is the set of arcs (links) $A = \\left\\{ \\left(i,j\\right)  \\in V^2 : i\\neq j\\right\\}$\n",
    "\n",
    "$c_{ij}$ is cost of travel over arc $(i,j)\\in A$\n",
    "\n",
    "$Q$ is the vehicle capacity (for ambulances this is an integer)\n",
    "\n",
    "$q_i$ is the load represented by a patient (in our problem this always 1)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\begin{aligned}\n",
    "       \\text{min} &  \\sum_{i=1}^{n} \\sum_{j=1}^{n} c_{ij} x_{ij} \\\\\n",
    "       \\text{s.t.} & \\sum_{i=1}^{n} x_{ij} = 1, & & j \\in N, j \\neq i \\\\\n",
    "       &   \\sum_{j=1}^{n} x_{ij} = 1, & & i \\in N, i \\neq j \\\\\n",
    "       &  \\text{if} \\ x_{ij}=1 \\ \\Rightarrow \\ u_i + q_j = u_j && i,j \\in A : j\\neq 0, i\\neq 0 \\\\\n",
    "       &  q_i \\leq u_i \\leq Q && i \\in N \\\\\n",
    "       &  x_{ij} \\in \\{0,1\\} && i,j \\in A\n",
    "\\end{aligned}\n",
    "\\end{equation}\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def problem_instance(n_patients, vehicle_capacity, random_state=None):\n",
    "    '''\n",
    "    Create an instance of the CVRP in the formulation\n",
    "    of a MIP.\n",
    "    \n",
    "    Parameters:\n",
    "    ------\n",
    "    n_patient: int\n",
    "    vehicle_capacity: int\n",
    "    random_state: int, optional (default=None)\n",
    "    \n",
    "    Returns:\n",
    "    -----\n",
    "        dict, keys = N, V, A, c, Q, q\n",
    "    '''\n",
    "    rand_gen = RandomState(seed=random_state)\n",
    "    \n",
    "    #n_patients = 1 so we include the depot\n",
    "    x_coords = rand_gen.random(size=n_patients+1) * 200\n",
    "    y_coords = rand_gen.random(size=n_patients+1) * 100\n",
    "    \n",
    "    N = np.arange(1, n_patients+1)\n",
    "    V = np.arange(0, n_patients+1)\n",
    "    A = [(i, j) for i in V for j in V if i != j]\n",
    "    c = {(i, j): np.hypot(x_coords[i]-x_coords[j], y_coords[i]-y_coords[j]) \n",
    "     for i, j in A}\n",
    "    Q = vehicle_capacity\n",
    "    q = {i: 1 for i in N}\n",
    "    return {'x_coords':x_coords, 'y_coords':y_coords, \n",
    "            'N':N, 'V':V, 'A':A, 'c':c, 'Q':Q, 'q':q}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimisation_model(instance):\n",
    "    '''\n",
    "    Create and return an istance of a gurobi optimisation\n",
    "    model for solving the CVRP.\n",
    "    \n",
    "    Parameters:\n",
    "    -------\n",
    "    instance: dict\n",
    "        problem instance for the CVRP\n",
    "    '''\n",
    "    mdl = Model('dialysis_transport')\n",
    "    \n",
    "    N = instance['N']\n",
    "    V = instance['V']\n",
    "    A = instance['A']\n",
    "    c = instance['c']\n",
    "    Q = instance['Q']\n",
    "    q = instance['q']\n",
    "        \n",
    "    x = mdl.addVars(A, vtype=GRB.BINARY)\n",
    "    u = mdl.addVars(N, vtype=GRB.CONTINUOUS)\n",
    "    mdl.modelSense = GRB.MINIMIZE\n",
    "    mdl.setObjective(quicksum(x[a]*instance['c'][a] for a in A))\n",
    "    mdl.addConstrs(quicksum(x[i, j] for j in V if j != i) == 1 for i in N)\n",
    "    mdl.addConstrs(quicksum(x[i, j] for i in V if i != j) == 1 for j in N)\n",
    "    mdl.addConstrs((x[i, j] == 1) >> (u[i]+q[j] == u[j])\n",
    "                   for i, j in A if i != 0 and j != 0)\n",
    "    mdl.addConstrs(u[i] >= q[i] for i in N)\n",
    "    mdl.addConstrs(u[i] <= Q for i in N)\n",
    "    return mdl, x, u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_solution(x, instance):\n",
    "    active_arcs = [a for a in instance['A'] if x[a].x == 1.00]\n",
    "    xc = instance['x_coords']\n",
    "    yc = instance['y_coords']\n",
    "    \n",
    "    for i, j in active_arcs:\n",
    "        plt.plot([xc[i], xc[j]], [yc[i], yc[j]], c='green')\n",
    "        plt.plot(xc[0], yc[0], c='r', marker='s')\n",
    "        plt.scatter(xc[1:], yc[1:], c='b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Changed value of parameter TimeLimit to 10.0\n",
      "   Prev: inf  Min: 0.0  Max: inf  Default: inf\n",
      "Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (linux64)\n",
      "Optimize a model with 60 rows, 255 columns and 480 nonzeros\n",
      "Model fingerprint: 0x3fec67ff\n",
      "Model has 210 general constraints\n",
      "Variable types: 15 continuous, 240 integer (240 binary)\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  Objective range  [9e+00, 2e+02]\n",
      "  Bounds range     [1e+00, 1e+00]\n",
      "  RHS range        [1e+00, 4e+00]\n",
      "Presolve added 495 rows and 105 columns\n",
      "Presolve time: 0.01s\n",
      "Presolved: 555 rows, 360 columns, 3214 nonzeros\n",
      "Variable types: 120 continuous, 240 integer (240 binary)\n",
      "\n",
      "Root relaxation: objective 5.126874e+02, 66 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0  512.68736    0   33          -  512.68736      -     -    0s\n",
      "H    0     0                    2283.8105615  512.68736  77.6%     -    0s\n",
      "H    0     0                    2057.5596703  512.68736  75.1%     -    0s\n",
      "H    0     0                    1158.8027970  512.68736  55.8%     -    0s\n",
      "H    0     0                    1068.0495026  515.22171  51.8%     -    0s\n",
      "H    0     0                    1013.0455384  515.98030  49.1%     -    0s\n",
      "     0     0  515.98030    0   34 1013.04554  515.98030  49.1%     -    0s\n",
      "H    0     0                    1002.8454755  515.98030  48.5%     -    0s\n",
      "     0     0  516.00220    0   34 1002.84548  516.00220  48.5%     -    0s\n",
      "H    0     0                     963.9684188  516.00220  46.5%     -    0s\n",
      "     0     0  519.71221    0   40  963.96842  519.71221  46.1%     -    0s\n",
      "     0     0  521.68589    0   39  963.96842  521.68589  45.9%     -    0s\n",
      "H    0     0                     954.9533692  521.95062  45.3%     -    0s\n",
      "     0     0  521.95062    0   39  954.95337  521.95062  45.3%     -    0s\n",
      "     0     0  527.31085    0   33  954.95337  527.31085  44.8%     -    0s\n",
      "     0     0  527.39467    0   30  954.95337  527.39467  44.8%     -    0s\n",
      "H    0     0                     922.3994737  527.39467  42.8%     -    0s\n",
      "     0     2  527.40031    0   29  922.39947  527.40031  42.8%     -    0s\n",
      "H  390   370                     915.6841643  533.53734  41.7%   8.1    0s\n",
      "H42176 22741                     915.2833974  656.55819  28.3%  12.3    4s\n",
      " 46058 24685  770.80127   62   16  915.28340  656.55819  28.3%  12.3    5s\n",
      " 173350 78024  761.29099   53   20  915.28340  656.55819  28.3%  11.5   10s\n",
      "\n",
      "Cutting planes:\n",
      "  Gomory: 13\n",
      "  Implied bound: 1\n",
      "  MIR: 4\n",
      "  Flow cover: 10\n",
      "  Zero half: 14\n",
      "\n",
      "Explored 174541 nodes (2009223 simplex iterations) in 10.02 seconds\n",
      "Thread count was 16 (of 16 available processors)\n",
      "\n",
      "Solution count 10: 915.283 915.684 922.399 ... 2057.56\n",
      "\n",
      "Time limit reached\n",
      "Best objective 9.152833974417e+02, best bound 6.565581885100e+02, gap 28.2672%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD7CAYAAABnoJM0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd1iTVxuH7yTsJTJlqIh7L1y17oEL96qj7rZa96RqHXWC1tGqdY/WUWedVbS1zrZu0bplOBGQvUeS7w8KlU9UIAlJ4NzX5XW173ve5/wSIE/OOc+QKJVKJQKBQCAoski1LUAgEAgE2kU4AoFAICjiCEcgEAgERRzhCAQCgaCIIxyBQCAQFHGEIxAIBIIijnAEAoFAUMQx0LaA9xEVlYBCkT3NwdbWgoiIeC0p+jC6rg90X6Ou6wOhUR3ouj7QfY3/r08qlVC8uHme7ei0I1AolG85gszruoyu6wPd16jr+kBoVAe6rg90X6M69ImtIYFAICjiCEcgEAgERZwPOgIfHx9atmxJxYoVefjwYdb1oKAg+vTpg6enJ3369CE4ODhX9wQCgUCgW3zQEbRq1YodO3bg4uKS7frs2bPp168ffn5+9OvXj1mzZuXqnkAgEAh0iw8eFnt4eLx1LSIigrt377JlyxYAOnXqxLx584iMjESpVL7zno2NjZrlF26GLj791rXN3i21oEQgEBRm8nVGEBISgqOjIzKZDACZTIaDgwMhISHvvSfIPTk5gfddFwgEgvyi0+GjtrYWOV63t7csYCV5Q5367hhvJtjoV1rGr8NUaas2+0XpPdQUQqPq6Lo+0H2N6tCXL0fg5OREaGgocrkcmUyGXC4nLCwMJycnlErlO+/llYiI+LdiZO3tLQkPj8uP7AJBnfqeGJ4kyPgwAK8M/6JMaicAle0XpfdQUwiNqqPr+kD3Nf6/PqlU8s4v0O8jX1tDtra2VK5cmaNHjwJw9OhRKleujI2NzXvvCXLPC4Pz3Db5Afv0OljJ3XhheFbbkgQCQSHlg45g/vz5NG3alFevXjFkyBA6duwIwJw5c9i+fTuenp5s376duXPnZj3zvnuCD+MXfJybpiuwkVfBI3EaLmnNiJY9IkEizlkEAoH6kehyz+KiuDV0/vlZ+h3rSRXbqjg8nowhZiRJXvO7xQgqpPTlwsR1WteoaXRdHwiN6kDX9YHua1TX1pBOHxYXNa6+uszAX/viXqwsuzrtx8bENute14NbCUu8ilKpRCKRaFGlQCAobIgSEzrCP69v88mxnjiaO7LH62A2JwDQo0JvHkc/4lb4TS0pFAgEhRXhCHSAx1GP6H2kC+YG5uzrfBhH8xJvjenk3hlDqSH7H+3VgkKBQFCYEY5AyzyNfULPw50BCfu7HKakZakcxxU3saFV6bb88mgfcoW8YEUKBIJCjXAEWiQ04RU9D3cmIT2BPV4HKWtd/r3je5TvRWjiK/58eaGAFAoEgqKAcARaIjI5gl5HuhCWGMaujvuoZlf9g8+0dWuPhaEl+x/uKQCFAoGgqCAcgRaIS42l75HuBMUE8lOHn/EoUT9Xz5kamNLR3YujgYdJTk/WsEqBQFBUEI6ggElMS2TAr334J+I2Gz1/pIlrszw93718L2JTY/jtyUkNKRQIBEUN4QgKkFR5KkP9BvD3yz9Z3Wo9nm7t82yjiWsz7E0dOCCihwQCgZoQjqCASFek88WpYZx++hvLmn9Pt/I982XHQGpA13LdOfXkBLEpMWpWKRAIiiIis7gAUCgVTDwzhqOBh5jXeBH9q3yqkr0eFXqz4fZajgYepl/lgWpSKRAItIk2G1GJFYGGUSqVzLgwlZ/v72Bqvel8XvNLlW3WdqiLm1UZkVwmEBQS3nQCSpQ5XtckwhFomEWX5rHp9npG1hzDJI9parEpkUjoUaE3F56f5VWCqEgqEBQWbhmv4bhlr2zOoCAQjkCDfHd9GSuuL2VglcHM+Wi+WovF9SjfGyVKDj7erzabAoFAe6SRyHOjP1CiRELBFpYUjkBDbP5nA/P/nkP38j3xbbpc7RVDyxUvT0372ux/KLaHBILCQJDRERSSNCwULgU+t3AEGmD3/Z14n5uEp1t7vm+5DplUppF5upfvhX/4DR5HPdKIfYFAUDCkEkeg8SGkSiOs5GUKfH7hCNTMgXsHGPfHKJq4NGND220Yygw1Nle38j2QIGH/I1FyQiDQZyo1u0E6iSgkqZgrnLOui6ghPeT009/ou68vdRw82NZhFyYGJhqdr4S5Ex+7NOXAo73ocKM5gUDwHsISw9h4ey2tSrcFYKJXazZ7tywwJwDCEaiNv0P+YsiJ/lR1qMquTvuwMMx7u7j80KNCb4JiArkRdq1A5hMIBOrlu+vfkiJPoXXpNgC4Fytb4BqEI1AD/mE36H+sFy4WrvgN8KOYsXWBzd3JvTPGMmNRckIg0ENexD1n6z+b6FOxH4npSYBwBHrJ/ch79DnaDWtja/Z1PoyDuUOBzm9lXIzWpT355dF+0hXpBTq3QCBQjeXXlqJEyaR60wiOCcTO1A4r42IFrkM4AhUIigmk1+EuGEqN2Nf5MM4WBR/2BRk5BeFJYVx4cU4r8wsEgrwTHBPEzvs/MrDKYEpaliIwOgA3K3etaBGOIJ+8jH9Br8NdSJWnsLfzIcoU084PEKB16bZYGRUTDWsEAj1i6dXFGEgMmFB3CgCBMQG4Wxf8thAIR5AvwhPD6XW4C5HJkez2+oVKNpW1qsfEwIRO7p05FniEpH/3GQUCge7yMPIB+x7uZmj1z3A0L0FiWiIhCS+1cj4AwhHkmZiUaPoc7cbz+Gfs7LiXWg51tC0JgO4VehGfFsep4BPaliIQCD6A75WFmBqYMab2BCBjmxm0c1AMwhHkifi0eD452pMHkffY0m4HDZ0/0rakLBo7N8HRrAT7RHKZQKDT3A7353DAL3xecxS2prZAxrYQILaGdJ3k9GQGHe/H9bCrrGuzhZalWmtbUjZkUhldy/fg9ycniU6O0rYcgUDwDnwuL8Da2JqRNUdnXQvKdARiRaC7pMnT+OzUEM4/P8PKFmvoVLaztiXlSM/yvUlTpHEk8JC2pQgEghy4+uoyJ5+c4Mta47LlGwVGB2Bv6oCFkaVWdAlH8AHkCjljTn/BiaBjLGqylD6V+mlb0jupYV+LstblOCAqkgoEOsmiy/OxM7VneI0vsl0Pig3U2rYQCEfwXpRKJVPPTeTAo73MbDiHYdU/07ak9yKRSOhRvjd/vrzAy/gX2pYjEAje4MKLc5x/foZxdSZibmie7V5gdIDWtoVAOIJ3olQqmfPnTH66u4VxdSYxts5EbUvKFd0rZHQ3+uWRaFgjEOgKSqWSRZfm4WTuzKCqw7Ldi0+LJzTxlVZzkYQjeAfLrvnyg//3DK02gukNZmlbTq5xL1aWOg51RWlqgUCHOP30FFdeXWKix9S3qhJrO3QUwEBrM+sw6/xX43N5Ab0rfsLCJkvU3l1M0yQ/rcVdk0308t2GpaIkUHB1zQUCQXaUSiWLLs+nlJUbn1Qa8Nb9oOiMiKEy+nxG8Mcff9C1a1e6dOmCl5cXJ0+eBCAoKIg+ffrg6elJnz59CA4OVnWqAmHH3R/5+uJXdHTvzIoWq5FK9GvRNHTxaZzTPgallBeGZ7NdFwgEBc+xwCPcCr/JFA9vjGRGb93PzCHQ260hpVLJ1KlT8fX15dChQyxZsoRp06ahUCiYPXs2/fr1w8/Pj379+jFrlu5vrxx8tJ+JZ8bQslRr1rbZhIFUPxdMJsriWCicCTQ6jJwUbcsRCIoscoUcn8vzKW9dgZ4V+uQ4JjAmAEezEgXWwyQnVP66K5VKiYuLAyAuLg4HBweioqK4e/cunTp1AqBTp07cvXuXyMhIVafTGKeCTzDq9xE0cGrEZs/tGMuMtS1JJcqn9EUhSeWxsTg0Fgi0xS+P9/Eg6j7T6s94Z+/ywGjtFZvLRKWvvBKJhBUrVjBq1CjMzMxISEhg3bp1hISE4OjoiEyW8cJlMhkODg6EhIRgY2OTa/u2tjl7SHt79SZd/BH0B0P9BlKrRC38Pj2OlbGVSvbUrS8/uKR/TFjqFQKMDuCS1gwLhUs2Xbqg8X3ouj4QGtWBruuD/GtMk6fx7bXF1HSsyZCGA965zRwcF4hXBa98z6OO91AlR5Cens66detYs2YNdevW5dq1a0yYMAFfX1+VhQFERMSjUGTvxWtvb0l4eJxa7ANcC71Cz8NdKFPMne3t9pISKyGc/NtXtz5VqJwymFDDK/xjsp4GiXOydOmSxpzQdX0gNKoDXdcHqmncfncbAVEBbO+wm4jXCTmOiUuNJSwhDCfjUvma5//1SaWSd36Bfh8qbQ3du3ePsLAw6tatC0DdunUxNTXF2NiY0NBQ5HI5AHK5nLCwMJycnFSZTu3cef0PfY/2wN7Unr1eh7AxsdW2JJV5MzrIRFmcSskDeG3gT9eeMVpUJRAULVLkKXx71Ye6jh60Kd3uneMyQ0e1eVAMKq4ISpQowatXrwgMDMTd3Z2AgABev35N6dKlqVy5MkePHqVLly4cPXqUypUr52lbSNMERD+i15EumBmYsa/zYRzNS2hbktp40xnIFc3w3H+ZmRe8aVmqNZZGqm17CQSCD/PTnS28iH/Oihar3xt+Hhit3WJzmajkCOzt7ZkzZw7jxo3LerGLFi3C2tqaOXPm4O3tzZo1a7CyssLHx0ctgtXBs7in9DzcBVCyr/NhSlmV1rYkjSGTyvBtuoz2+1vhe2UR8xov0rYkgaBQk5CWwPJrS2ns3ISmrs3fO1YXQkdBDQllnTt3pnPnt6txli1blr17da/4WWhiKD0PdyYuNY5fuh6jfPEK2pakceo4evBp1aFsvLWWPhX70cJed/ooCASFjc3/bCA8KYzN7bZ/MBk1MCYAJ3NnzAzNCkhdzuhXtpSKRCVH0vtwF0ITQtnVaR/V7WpoW1KBMaPBLIqbFGfq2QkolAptyxEICiVxqbGsur6cVqXa0MCp4QfHa7vYXCZFxhHEp8bR92h3AmMC+LHDLuqVaKBtSQWKtUlxZjWax9XQy2y5sUXbcgSCQsla/9VEpUThXX9mrsYHabFh/ZsUCUeQlJ7EgF/7cCvcn42eP35w366w0qdiPxo6fcTU36YSkRShbTkCQaEiMjmCtf6r6ejemZoOtT84PiYlmojkCMqIFYHmSZWnMuzEQP56eZHVrdfj6dZe25K0hkQiwafpMmKSY1jw9xxtyxEIChWrb3xHfGoc0+rPyNV4Xag6mkmhdgTpinRG/TaC356eZGnzlXQv30vbkrROZdsqTGg4ge33tnHl1SVtyxEICgWhiaFsvL2W7uV7Ucmmcq6e0XbD+jfRz6pquUChVDDxzBgOB/zC3I8WMrDKYG1LyjU5VQpVZxnp2c1ns+PWTqaencipXmf1trieQKArfHftW1LlqUyp553rZzJzCEpbuWlIVe4pVJ8AmR+gSpTcMd5EsPFRJnt4M7LWaC0ryz3vKhc9dPFptTkDCyML5n/swzC/gWy+vZ7Pao5Si12BoCjyIu452+5spm+l/rhbl8v1c4ExAbhYuGJqYKpBdbmj0DiCNz9AHxrvJNj4KGVSOnPvdAOonz+bcoWcVEUqqfIUUuVppClSSZX/+0+RSpo8lVRFGqnylH/vpWEaLiM8Kpo0edr7n/13/JvPpilSuW4WgoJ04qUvKZvSFfe0t3M01EEn9860LNWaxZcX0LlcN0qY61b5D4FAX1h2LaO22iSPaXl6LihGN0JHoRA5gkxumfzAUyM/zOUupEriuGG6nGF+2/79QE4lTZH27wdyCqmKNNLkqaT8+6H+5odzqiJVI/H2hlJDDKVGGMuMMJQZYSQ1wlBmiLHMGEOpEWkkIcUABWncN9mOnbwGVgo3teuQSCQsbLKEZj83ZPbF6axrK0JKBYK8EhgTwM57PzG42jBcLUvm7dnoALzKdtOQsrxR6BxBlOw+UqUBaZJ4IgxuI1Ua8jCy2L8fuoYYyYwxNjDB0sgSI5lx1gexkdQII1nGP0OpEUayjLGZ/53x4W2ModQwa0xOH+aOdsWJj0nNwXbGsx/KNMxc2SRLIjlvPomrZov4OH4pRqi/XK97sbKMrTORJVcW0a/ypzQr2ULtcwgEhZmlVxZjJDNifJ3JeXouKjmSqJQonTgohkLoCJolrHzr2uYxBdev197eknCJ6qV1TZQ2eCRN40+zmdwwW0b9xNwlqOSVMbUnsPfBz3ifn8SZPn/pfUMegaCgeBB5n/0P9zCq1tg8F63MihjSka2hQh0+qo+8eSBcXF6JasmfEW5wg7LNLmpkPhMDExY3XUpA9GPW3PhOI3MIBIUR3ysLMTe0YEyd8Xl+VpdyCKAQrQg2e7fUeNhlQZFdc0smnUniuxvLqGFfk87l1L+n2LJUG7zKdmX5tSV0r9BLJ8LZBAJd5na4P0cCDjLJY1q++pgERgcgQULpYm7qF5cPCo0jAP380M8NC5v4ci/iDmNPj6J88YpUtq2i9jnmNV7E709OMf38FLZ32PPBswyBoCiz+PJ8rI2tGVkzf6HpgTEBuFqW1JmtWLE1pAcYy4zZ3O4nLIwsGHT8E6KTo9Q+h7OFC1PrT+fUEz9OBP+qdvsCQWHhyqtLnHrix+ja47EyLpYvG0ExATpRYyiTQrUiKMyUMHdik+dPdD/UkZG/DWd7hz3IpDK1zjGi+hfsvr+DGeen0tS1OeaG5mq1L9A/Cst2qzpZfGk+dqb2DKv+eb6eVyqVBMYE0q1cDzUryz9iRaBHNHBqyIKPffn96Sl8ryxQu31DmSE+zZbzPP4Zy68uUbt9gX7xviz3osr552c5/+Is4+tMyvcXpcjkSGJSonUmdBSEI9A7BlUdyoDKg1h+bSlHAw6r3X5Dp0b0rdSfNf7f8SDyvtrtC/SPKOkDggyPaVuG1lEqlSy6NA9ncxc+rTo033YCYx4DuhMxBMIR6B0SiYRFTZdS19GD0b9/zv3Ie2qfY1ajeVgYWuB9bhJKpVLt9gW6j1KpJFJ2lyumC7hoMY07JhtIJ1nbsrTKb0/8uBp6mYkeUzExMMm3nf8a1ue+LpGmEY5ADzGWGbPZczvmhuYMPt6PmJRotdq3M7VjRsM5XHx5nv2P9qjVtkC3kSvkHAs8QocDrfnTfDqRsvs4pjUACSRLi24zI4VSweLLCyht5cYnlQaoZCsoJgCpREopq9JqUqc6whHoKU4Wzmxq9xNP454w6rcRaq+LNKDyIOo41GX2xRlqdzQC3SM5PZkf72yh8S4PhpzoT3hSONWSPqN1/EbcU70yxkhea1ml9jgWeJjbr/2ZUu8rDGWGKtkKigmkpGUpjGRGalKnOsIR6DENnRox/2MfTj3xw/fKQrXalkll+DRdRkTyaxZfnq9W2wLdISo5kuVXl1Dnp6pMPjsOSyMrNrTdyt/9rnN60lJkGGOiyEiYSpJmOIKiFjUkV8jxubyACsUr0qN8b5XtBcYE6tT5AIjwUb1nSNXh3Aq7ybKrvlS3q0lHdy+12a7pUJsh1Yaz5Z+N9K3YP1d9WAX6wbO4p6zzX832uz+SmJ5Ay1KtGV17PI2dm2RLJtzs3ZIUeQol142kXTMrJnkULScAcODRXh5GPWCT548qh2xnhI4GUNfRQ03q1INwBHqORCJhcdNvuRd5h9G/f0556wpUsKmoNvve9Wdy+PFBpp6bwK/df1d77oKgYLn9+harb6zk0OMDSCQSupXryahaY6lqV+2dzxjLjLEztedl/IsCVKobpMnT8L2ykGp2NejornpvkNdJr4lLjdW5FYHYGioEmBiYsKXdDkwNTBl04hNiU2LUZruYsTVzGy/gRth1tt/bpja7goJDqVRy5tlpeh3uQqs9H+MXfJwRNUZypf8tVrde/14nkImLhSsv4p8XgFrdYuvNrTyJDear+jORSlT/uNSlPsVvIhxBIcHZwoVNnj/yJDZY7YfHPcr3prFzExb8PYfwxHC12RVolnRFOvsf7qHV3ib0PtKVe5F3mdlwDjc/vcs3jRfiYumaa1vOFi5FbkWQnJ7MN+e+oa5jPVqX9lSLzSAdKz+diXAEhYhGzo2Z13gRJ5+cYOmVxWqzK5FI8Gm6jPi0eOb9PUttdgWaIT4tnvX+ayj3XTlG/jaclPRkVrRYzbWB/zC2zkSKGVvn2aaLhQsvipgj+OnuFp7HPuerBl+rrQhjYHQAMomMkpa6EzoK4oyg0DG02mfcDLvB0quLqWFfi3ZlOqjFbgWbioyqOZbvbiyjX6WBNHT+SC12BeojLDGMTbfXsuWfjUSnRPNxqY+Z39iXNqU9Vd7WcLZwJS41lrjUWCyNrNSkWHdJSEtg+bWltHBrQVPX5mqzGxgTQCmr0iqHoKobsSIoZEgkEpY0W0Et+9qM+m0Ej6Ieqs32BI8puFqUZNq5iaTJ09RmV6AaAdGPmHRmHHV/qsqKa9/ykXMTjnU/xfkh5/F0a6+WvW0XCxeAIrMq2HR7Pa+TwpnfUr2h00E6GDoKwhEUSjIPj00MjBl0/BPiUmPVYtfc0JwFTXy5F3mXDbfXqsWmIP9cfXWZwcf789FOD/Y82Enviv34s99VtrbfQb0SDdQ6l7NFxnnCyyJwYBybEsOqG8tpXaotH5VU38o3M3RUOAJBgeFi6cqGttsIignky98/V9vhcTu3DrQt3Q7fywuL3OGhLqBQKvALPo7XL54ZZSBenmd83UlcG3iHb5uvpKx1eY3MW5RWBGv9VxOdEo13A/X2CQ9LCiMhLZ4yxdzValcdCEdQiGns0oRvGi/kRNAxll31VYtNiUTCgia+KJRyvr74lVpsCj5MijyFnfd+osmu+gz8tQ8v4p4zv/Firn96l68azMLBzEGj85cwd0IqkRb6ENLI5AjW+q+mk3sXatjXUqvtoGjdDB0FNRwWp6SksHDhQv766y+MjY2pVasW8+bNIygoCG9vb6Kjo7G2tsbHxwc3Nzc1SBbkheHVv+Bm2A18ryykhn1N2rq1V9lmaSs3JtSdwqLL8zj99BQtS7VRg1JBTnX+l0+ow7Y7W9hw6wdCE19Rza4GP7TeSOey3Qr0wNFAaoCjWYlCvwpcdWMlCWnxTK0/Xe22M3MIdKkzWSYqO4IlS5ZgbGyMn58fEomE168z6pHMnj2bfv360aVLFw4dOsSsWbP48ccfVRYsyBsSiYSlzVfyIOo+I38bgV+PP7C3r6Oy3VG1x7Ln4S68z03mXN9LKpXlFbztBJIkrwkyOkLljZ+QLkmimWsLvm+1lmauLbTWT9q5kIeQhiaGsun2OnpU6E0lm8pqtx8YHYCB1ICSlqXUbltVVNoaSkhI4ODBg4wbNy7rl9POzo6IiAju3r1Lp06dAOjUqRN3794lMjJSdcWCPGNqYMqWdtsxkhoy6PgnxKaofnhsLDPGp+kygmOD+P7GcjWoFABESR9y1dSH0xafE2R0BIf0evze6zx7Ox+iecmWWnMCkJFdXJgPi1deW0qqPJXJ9bw1Yj8wJoDSVm4YSHUval8lR/Ds2TOsra1ZtWoV3bt3Z+DAgVy9epWQkBAcHR2RyTLq0shkMhwcHAgJCVGLaEHeKWlZig2e2wiMCWDQwUFqOTxu6tqc7uV78t31ZVnLXoFqXDVbxCuDS7iltqdF/A/USZpIdfua2pYF/JddXBibFT2Pe8aPd7bwSaUBGovq0dWIIVBxayg9PZ1nz55RpUoVpk2bhr+/P1988QUrV65UizhbW4scr9vbW6rFvqbQVX3d7DuyNHkpE/wmUNfpe2Y2VT0qYpXXd5xa5cfsv7053v+42r6x6up7+Cbq1pgsiSJFEkW5lF5USu2vlnnUqbFiibIk+Schs0jD1sxWLTZ15ec8/a/lIIEFnt9gXyy7JnVoVCqVBMcG0rZca7W/ZnXYU8kRODs7Y2BgkLUFVLNmTYoXL46JiQmhoaHI5XJkMhlyuZywsDCcnJzyZD8iIh6FIvu3D3t7S8LD41SRrVF0XV8/96FcrX6VWX/Mwt20Im3c2qlkzwALvOvPZMaFaWy5tAOvsl1U1qjr7yFoRmOowRWQgHN642zX8zuPujVaYQeA/5P7VLerobI9Xfk5B8YEsOXmFoZUG45pavFsmtSl8VVCCIlpiZQwLqnW1/z/+qRSyTu/QL8PlbaGbGxsaNCgARcvXgQgKCiIiIgI3NzcqFy5MkePHgXg6NGjVK5cGRsbG1WmE6gBiUTCeq/1VLWrzsjfRhAY/Vhlm0OqjaCaXQ1mXphGfOr7f8mHLj791j9BRt3/UIPLmCocsFSUznZdV8jMJShskUNLLi/CSGbEuLqTNTZHZp/iMla6l0MAasgjmDt3LuvWrcPLy4uJEyfi6+uLlZUVc+bMYfv27Xh6erJ9+3bmzp2rDr0CNWBmaMbWdjswkMoYdLzfBz+8P4SB1ADfpssISXjJ0qs+7xyX+aEvJ5VU4t66XpSJT4snxuQ2A2v1YIt3KzZ7t9QpJwAZh8VAocoluB95jwOP9jKs+uc4mjlqbB5dLT+dicrH1yVLluSnn35663rZsmXZu3evquYFGqKUVWnWt91K7yNdGXN6JJs9f1Jpf9+jRH0GVhnMOv/V9KnYj8q2VXIcl0ocf1iMQilR0DbuR6SIRjcAZ56eJkWeQrsyHbUt5Z3YmzlgKDXkZVzhWRH4XF6AuaEFo2uP0+g8gTEBGEmNcLUoqdF58ovILC7CNHVtzqxG8zgWeJjvri9T2d6MhrMpZlyMqecm5BhZkk4Sl83mkS5JIF2SwEvDcyrPWVg4EXwMa2NrGjrpblVXqUSKk7kzLxMKhyPwD7vBscDDfFHzS2xM1HP4/S4CozNCR3W1w59wBEWckTVH0718TxZe+obfn5xUyZaNiS2zGs3jUshf7H6wM9u9FHkKV80WEy17TJ2kKVjJ3Xlo/DMK0lWaszCQrkjnVPAJWpf21MkY8zcpTA1qFl+eT3Hj4nxR80uNzxUUE6Cz20IgHEGRRyKRsKz5KqrYVuOL34arnA/Qt1J/6pVowDd/fU1UckYCoVwhZ+Sp4bw28Kdm8pc4pTeiYko/EqWhPDMU5wOXQ/4mKiWK9jq8LZRJRnax/p8RXAr5my3X6AEAACAASURBVN+fnuLL2uOxMi6m0bkUSgVBMYE6WVoiE+EIBBmHx+13IEXC4OP9iE+Lz7ctqUSKb9PlRCVHsfDSPJRKJZPPjuNo4CHmfrSQkmmtAHBIr4t1ekUeGe9mzWTd3Q4pCI4HH8NIakSLkq20LeWDuFi4EhL/Uq2tULWBz+X52Js6MKz6ZxqfKyT+JcnyZJ1NJgPRoUzwL6Wt3FjXdgt9j3Zn3OlRbGy7Ld+Hx1XtqjG8xhes919DYloCex/+zIS6kxlZazQj3yjoeO75Enoe7sz2u1sZXuMLNb0S/UKpVHIi6BhNXJthYaQbyVXvw9nChVRFKq+TXmu84mkmOUWVqRJRde75GS68OMeCj30wNzRXRVquCIoNBHQ3YgjEikDwBs1LtmRmw7kcCTjI9zdWqGRrar2vMDe0YO/Dn/m0yhC863/91pgmLs1o7NyE5deWkpiWqNJ8+sqDqPs8iQ3W6WihN3HOyiUomO2hd4UW5zfkWKlUsujSPJzNXRhYZYgq0nJNZg6BLq8IhCMQZOPLWmPpWq47C/6ew+mnv+XbzsHHB4hPy8gVqGhTOcfVhUQiYVqDmYQnhbH5nw35nkufORF0DABPNZQHLwi01aAm0PAwzwx+Q06qSnZOPTnBtdArTKo3rcAq5gbGBGAsM85yorqIcASCbEgkEpa3WE0lmyp8fmooQTGBebZx+PEvTD4zjlal2tDEpRk+lxcQmhia49iGTo1oWao1q24sV1tLTX3iRNAx6jjUpYR53sqvaAtttKx8Zniau6ab8TdbxUnLT7luupSXBhfyfJalUCpYfHkBblZl6Fux/4cfUBOBMQG4WZVRS+9oTaG7ygRaw9zQnG3tdyIBBh/vT0JaQq6f/ePp74z8bTj1SjRgk+dP+DZbRnJ6EnP/fHeBO+/6M4lMjmT9rR/UoF5/eJUQwvWwa3qzLQRga2KLicykwFYEoQZXuGWyCtv0angkzMAlrSmvZbe5braUypvL8Omvffn5/o6sCLX3cTTgEP+8vsWUel8VaFOfoOgAyujw+QAIRyB4B27FyrCuzRYeRN1j/Okvc1V6+MqrSww50Z8KxSuxo+MezAzNKGtdntF1xrPv4W4uvjif43O1HOrQvkwn1tz8Pld/0IWBoYtP039txjnMhd9t9abMhkQiwcnCuUBWBH+H/MU10yVYKcpQL3EGJeT1qJE8ijbxm2mUsIBBVYdy+/Utxp4eSZUtZel5uAtb/tlIaMKrbHaGLj7NkMWnGH98Bhbykvz6i2aTx95EoVQQHBuk0+cDIByB4D20KNWK6Q1mcyjgAKtvfvfesXcj7tDvWC8czBz52esAxYyts+6NqzOJUlZuTDs3kVR5znu80+rPID41jh9urlLra9BFMj/0XxlcwkzhiKWiVLbruo6LhavGVwR3I+4w8Nc+lCleivqJszDANOueBBmHpoxh/sc+XB94h5M9zzC69nhexD9j2rmJ1NhWkU4H2vLDzVX09dkFwAvDc8TLnlMx5RMkyArsvX4R/5wUeYrOOwIRPip4L2Nqj+dW+E3m/z2banbVaV7y7bC94Jggeh/piqmBKXs7H3qreJepgSmLmyyh37FerPVfzdg6E96yUcW2Kl3LdWf9rR/4quUUJG/84RdG0kkh3OAGjmkNkKC9rmP5wdnC5Z2rO3XwNPYJfY92x9TAlN2dfqGUVel3jpVIJNRyqEMthzpMbzCLB1H3ORZ4mGOBR5j953SwBCt5GZIk4VjKy1AivaHGdOdEoA43rH8TsSIQvBeJRMKKlqupWLwSn50cTHBMULb7oQmv6HWkC2nyVPZ4HaS0lVuOdlqX9qRDGS+WXfXhWdzTHMdMqTedZHkSPhfeXcG0sJAkDQWJEsd01ftHFzQuFi68SghBrpCr3fbrpNf0OdqNpPQkfu504L1O4P+RSCRUsqnMJI9pnO59gcv9/amcPBgJMtIk8ZRJ6YSkgD/ysqqO6viKQDgCwQexMLRgS/sdKIHBJ/47PI5KjqT3ka6EJ4azs9O+Dzb8nv/xYgBmXsi5J2y54uXpXfET1lxdQ0j8S7W+Bl0jUZqxj22pcNOukHzgbOGKXCknNPHVhwfngfjUOPod7cGLuOds77CHKrZVVbLnVqwMZVO70iRhKR3i9lMqveAzt4NiAjE1MNX5qDDhCAS5wr1YWda12cS9iDtM/GM08anx9D/Wm4Dox2xrv5O6jvU+aMPVsiST6nlzPOgoJ4OP5zhmksc00hXprLi+VN0vQadIkGaE05orS2hZSd75L5dAfQfGKfIUBp3oz+3Xt9jouY0GTurdwtFWufMgPQgdBeEIBHmgZak2TG8wi18e76fN3qZcD7vK2jabaVayRa5tfF5jFBWLV2L6+ak5ZhOXtnJjWO1hbL+7jaexT9QpX2fY7N2SRMkrDJRmGCots13XB/7LJVDPgbFcIWf0b59z/vkZlrdYRVs1Jte96z0tqPc6MDpAp4vNZSIOiwV54sta49j8zwYCYh7zRY3RdCrbOU/PG8mM8Gm6jK6HOrDy+lK+ajDrrTEzm85k682tLLvqy4qWq9UlXacoXyEdi4RybPlS9wvN/T/qzC5WKpXMuDCVQwEHmN1oPn0rqT/RS1sOVq6Q8yQ2GM8yHbQyf14QKwJBrlEqlUw7P4mQhJfYmzqw5+HOfH1r/8jlY3pV6MuqGyt5HPXorfuuVq4MqjqU3Q92EhD99v3CQHBM0DsP1nUdK6NimBtaqCWX4NurPmz+ZwNf1hrHl7XHqkGd7vA8/hmpilSdPygG4QgEeWDB33P56e4WxtaeyJFuJ0hXyBl8on++CsbN/mg+pgZmTDs/KcdktTF1JmIsM2bJlcXqkK5TKJQKnsY9wc2qjLal5AuJRIKLhYvKK4K1V9fie2UhfSr2Y1ajb9SkTnfQh2JzmQhHIMgVq26s5LsbyxhYZQgzGs7G3bocP7TewJ3Xt5l0ZmyuMo/fxMHMgekNZ3H++RkOPt7/1n1HM0eGVf+cXx7t417EXXW9DJ3gVUIIKfIUvV0RQGansvyvCI4EHGTUsVG0Ld2OZc2/V6lftq6i6w3r30Q4AsEH2XH3R77562u6lO2Ob9NlWX+0bdzaMa3+DPY/2sO6W3nfyx9UZSg17Wvz9cWviE2Jeev+l7XHYm5oge+VhSq/Bl3iSWwwkBHeqK+okl187vkZRp4azkclP2J9260FWvenIAmKCcDMwBxHM92PDBOOQPBejgQcYtLZsTQv2ZLVrde/1Xx7fN3JdCjjxdw/v+b887N5si2TyvBtuozwxLAcP+xtTGz5ouaXHAs8jH/YDZVehy6RmZSn7yuC8MSwd5YMeRf+YTcYdLwfZa3LceSTI5gZmmlIofbJaE/prherHeEIBO/k7LM/GHlqGHUcPNjSbgdGMqO3xkglUla1WktZ63J8dnLwO7OG30Vtx7oMqjqUjbfXcfv1rbfuf1HzS6yNrfG5vCDfr0PXeBIbhEwiw9WipLal5BsXC1eUKAlJyH3iX2D0Yz451gMbExt+7nSA4qbFNahQ+wTqeMP6NxGOQJAj10Kv/PvNrTw7O+59b0s/CyNLtrXfSaoijcHH+5OUnpSnuaY3mIWNiQ1Tz054qxeulXExRtcez29PT3I55FK+XouuERwbjItlSb3eEvmvU1nutodeJYTQ+0g3lEole7x+wcnCWZPytE66Ip0nscGUsXLXtpRcIRyB4C3uR97jk6M9sDezZ4/XL1ibfPibW1nr8vzQegO3X/vn+fDY2qQ4sxvN51roFXbe++mt+8Oqf46dqT0+l+fn6XXoKk9i9Td0NBOXf5PKcpNdHJMSTZ8j3YlIjmBXp/2UtS6vaXla51ncU9IV6WJFINBPnsQG0/tIV4xlJuz1OoSjee4Putq6tWdqvense7ibDXlsMtO74ic0dPqIeX/N4nXi62z3zA3NGVdnIudfnM3zOYQu8iQ2WG9DRzP5b0Xw/q2hpPQkBvzah8fRD9nabge1HPSvyF5+CNKTYnOZCEcgyCI0MZReh7uQnJ7EHq+D+YpqmegxlXZlOjL7zxl5KlUskUjwabqMuLQ4vH97uyjdoKrDcDJ3ZvHl+XkOVdUl4lPjeJ30Wu9XBOaG5lgbW783hDRdkc5nJwdzOeRv1rTekKdSJPpOZg6Brncmy0Q4AgGQuXzvRlhiKDs67qWybZV82ZFKpKxutY4yxdwZcXIQz+Oe5frZyrZV+LzGl2y6semt8wATAxMm1J3ClVeXOP30VL606QLB/4aOltHj0NFMnC1c33lGoFQqmXRmLH7Bx1nUdCldynUvYHXaJTAmAHNDCxxMHbQtJVcIRyAgMS2R/sd68yjqAVva7aBeiQYq2bM0smJb+10kp6cw5MSAPB0eT6o3DVcrV6aem0C6Ij3bvX6VB1LKsjSL9HhVUBhCRzN5X3bx/L/nsOv+diZ7eDO02ogCVqZ9AmMCcC9WVi9CR0E4giJPqjyVYX4DufLqEj+03kiLUuopgla+eAXWtN6Af/gNppwdn+sPbgtDC1a2W8ndiH/YdHtdtntGMiMm1/PmVvhNfg06qhadBU1WMpmenxEAOJnnnF285ub3fH9jOYOrDmNKva+0oEz7BEYH6M35AAhHUKSRK+SM+f1zfn96iqXNV9K5XDe12m9XpgOTPbzZ82AXm/9Zn+vnulXqRqtSbVh8ecFbDWp6VuhDWety+Fyer5EOWZrmSWwQxY2LY2VcTNtSVMbFwoWI5IhsK77d93cy588ZeJXtyqImS/XmG7E6SZOn8SzuKe7W+hE6CsIRFFmUSiVfnZ/ML4/3M7PhXAZWGayReSbX88bTrT1fX/yKv15ezNUzEomEhU2WkK5Iy+g7+wYGUgOm1pvO/ch7HAo4oAnJGiW4EISOZpIZORTy7/bQqeATjP/jS5q4NGNN6w1vZaEXFZ7FPUGulOtFH4JMhCMooiy+PI+tdzYxuvb4HJvJq4uMw+P1GQ1n/AbyIi53hcrKFHNnXJ1JHHx8gDPPTme716VcdyrbVMH38sK3zhF0nSexwXpdY+hNXCwzcwlecDnkEsNPDqKaXQ22td+JscxYy+q0R2axuSLpCFatWkXFihV5+PAhADdv3qRz5854enoydOhQIiIi1DWVQEWW/bWM5deWMqDyIL5uOFfj81kZF2Nbu10kpScz1G8AyenJuXpudO3xlCnmjve5SaTIU7KuSyVSptWfSWBMAHsf/Kwp2WonXZHOs7inlC4E5wPw34rgaugVBvzaCydzZ3Z23IeFkeUHnizc6FP56UzU4gju3LnDzZs3cXbOSBtXKpVMmTKFWbNm4efnh4eHB0uXFu4etPrCrnvbmXRyEl5lu7Kk2YoC28OtYFOR1a3WcyPsOlPPTcjV4bGJgQmLm3xLYEwAq2+szHavfZmO1LSvzdKri/Nc+ExbvIx/QboivVAcFAM4m2c4gu+vL8dYZsJur1+wN7PXsirtExgTgKWRFXamdtqWkmtUdgSpqal88803zJ49O+tD5fbt2xgbG+Ph4QFA3759OXHihKpTCVTk18CjTDgzmjbubbSyh9vBvRMTPaby8/0dbLmzMVfPtCjVis5lu7Hi2tKs0EvIOEf4qsFMnsU9zbEshS4SHPtv6GgxN+0KURPxafFIJTJS5Sns9vql0Jx9qIq+hY6CGhzBypUr6dy5MyVL/ldJMSQkJGt1AGBjY4NCoSA6OlrV6QT55Pzzs3x2cjC1HepwoM8Bre3hTq03nTalPZl5YRp/v/wzV8/Ma7wImdSA6eenZFtJtCjZmvolGrLsmm+eC91pg8IUOhqfFk//Yz1RKhVUt69JFduq2pakMwTGBOJeTH8ihkDF5vU3btzg9u3bTJ48WV16smFra5HjdXt73d6D1DV9V15cYdCJTyhvW56Tg/ywMLLAQosr+D19f6b+hvoMP/Up1z67hquV61tj3nwP7e0rMq/FN0w8OZE/I0/TtVLXrHs+notosa0FB57sZHzD8QWiPyeNuSEs7QWGUkNquFUssNWYJn4XU+Wp9N/Zg1vh/tR3qU9CWkK+59G1v5WcyIvGVHkqz+OeMqjWwAJ7beqYRyVHcOXKFQIDA2nVKiMJ6dWrVwwbNoyBAwfy8uV/8d+RkZFIJBKsra3zZD8iIh6FIvtesr29JeHhcarI1ii6pu9h5AM6H2yHjbEtuzocQB5vCKZoWaOMzW130G5/S7rs6MqhbieyrVByeg/7ug9mo+1mRh8bQy2rhlllsaua16WJa3MWnFtIl1J9sDDM+cuDusnPz/neq4eUtCxFZETeezznB038LiqUCr44NZRTgaf4ruUP+IffYO+D3fmaR9f+VnIirxofRz1CoVTgaOhaIK/t//VJpZJ3foF+HyptDX322WdcuHCB06dPc/r0aUqUKMGmTZsYPnw4ycnJXL16FYCff/6Z9u3bqzKVIB88i3tKryNdkEkM2NP5ICXMnbQtKYuKNpX4vuVaroddw/tczg3s38RAaoBP02W8iH/Osqu+2e5515/B66RwNt/OfdKaNtD30FGlUsmMC1M5+PgAXzf6hr6V+uNs4Upsagzxqbr9gV5QBMY8BvQrYgg0lEcglUrx9fVl7ty5tG3blitXrjBp0iRNTCV4B2GJYfQ63IXE9ET2eB3UyV/MTmU7M6HuZHbc+5FtdzZ/cHwDp4b0qzSQH/y/50Hk/azr9Uo0oE1pT1bdWJFj72NdQd+TyZZfW8Km2+sZWXMMo2uNAzKyi4F89y8ubOhjDgGo2RGcPn2aChUqAFCnTh2OHDnCyZMn2bJlC3Z2+hNKpe/EpsTQ92h3QhJesr3DHqraVdO2pHcytd4MWpVqw4wLU7kU8vcHx3/d6BssDS2Zdm5itlXEtPoziE6JZq3/ak3KzTfRyVHEpETjpicdq/6fbXc2s/jyfHpX/ITZH83LiohxzkODmqJAYHQAxYytsTGx0baUPCEyiwsZiWmJDPi1Dw8i77Gl3XYaODXUtqT3IpPK+KH1RlwsXBnmN5BXCSHvHW9rasvMRnP58+UF9j3cnXW9hn0tOrl3Ya3/aiKTdS95MSt0VA9XBEcCDjH17ATalPZkefNVSCX/fWy45LFlZWEnM2JIn0JHQTiCQkWaPI0RJwdxKeQvVrdaT8tSbbQtKVdYmxRnW/tdxKfGM+TEAFLSU947vn/lT6nr6MHsP2cQk/JfSPLU+tNJSItn9Y3vNC05z2SFjurZGcGFF+cYeWoYHiXqs6Httrf6LJcwc0KCRKwI/iUoJkDvtoVAOIJCg0KpYMzpLzj1xA/fZsvpWr6HtiXlicq2Vfi+1Q9cC73CmONj3jtWKpHi23Q5kckRLLo0L+t6JZvKdC/fi4231xKaGKppyXkiMxmulFVpLSvJPbfCb/Lpr5/gbl2W7R12Y2Zo9tYYQ5khjuYlxIoASE5P5nncM508j/sQwhEUAjKjOQ482suMBrMZVHWotiXlC6+yXRlbeyIbrm/gxztb3ju2un1NhlX7jC3/bMQ/7EbW9Sn1vEmVp/L99WWalpsnnsQGY2dqX2DhraoSGBNA36M9sDa2ZnenXyj+nj3v9zWoKUo8iQ1GiVJvGta/iXAEhQDfKwuzojnG1pmobTkq8VWDr/Es68lX5ydz5dWl946dVn8G9mYOTD03Ias3gbt1OfpW6s/WfzblutJpQaBPDetDE17R+0g3lEoFe7wO4mTh/N7xGS0rdee91haBetaw/k2EI9Bz1vuv4durPvSrNJA5H83Xu0Oq/0cmlbGzx06cLVwYeuL9h8dWxsX4pvFCboRd56e7W7OuT/SYihIly6/pTqFDfQkdjUmJps/R7rxODGdnx32UK17+g884W7jwMv6F3rYPVRdBMYGAcASCAmb3/Z3MvOhNhzJeLG2+Uu+dQCY2pjZsa7+LuNRYhp4Y+N7qot3K9aSJSzMWXJpLeGI4ACUtSzGwymB23v8xW6E6bZEqT+VF/HOdPyhOSk9i4K99eRT1gK3td1DbsW6unnOxcCExPZHolCgNK9RtAqMDKG5cHGuT4tqWkmeEI9BTTgT9mtENyrU5a9tswkCqUrUQnaOKbVVWtlzD1dDLzLgw7Z3jJBIJi5t+S2JaAt/89XXW9fF1J2MgMeDbqz4FIfe9PI97ikKp0OkVQboinc9PDc2KOGtesmWun3X5N5fg5f+1FS1qBMUE6OX5AAhHoJdcfHGeEScHUcO+Jtva7cDEwETbkjRCl3LdGV17PNvubGL73W3vHFe+eAW+rDWO3Q92ZrXDLGHuxJBqI9j78GceRT0sKMk5EpwVOqqbyWRKpZLJZ8ZxIugYC5v45jnizDkrl6BonxME6mnoKAhHoHf4h91g4K99cbMqw65O+wt9N6gZDWbTvGRLvM9N4uqry+8cN77uZEpalmLauYmkydMAGFNnAiYyU5ZcWVhQcnMkM5nMTUdXBAv+nsvO+z8x0WMqw6p/nufnXSz+a1lZVElKT+JF/HO9PB8A4Qj0ikdRD+l7tDvFTYqzx+sgNia22pakcWRSGevabKaEhTND/Qa+Mz/AzNCMhU2WcD/yHutv/QCAnakdn9ccycHHB7jz+p+ClJ2NJzHBmMhMcDBz1JqGd7HWfxXf3VjGp1WGMq3ejHzZsDd1wEBqUKRzCTLPosTWkECjPI97Rq/DXZBIpOzNRUhfYaK4iQ1b2+0gJiWaYe85PPZ0a087tw4subIoK3R0ZM0xWBkVw+fKgoKUnI3MiKE3SzPoAnsf/Mysi9Pp5N4Fn6bf5jvYQCaVUcLMqUhnF+tz6CgIR6AXhCeG0+tIF+LS4tjt9Qvu1uW0LanAqWZXnRUtVnP51d98fdH7nePmf+yDEgVfX/wKyChfMarWGE4EHeNG6LWCkpuNJ7HBOndQ/PuTk4z7YxQfuzTlhzYbVW6UkxlCWlQRjkCgUeJSY/nkWA9exr9gR8e9VLeroW1JWqNb+Z6MqjWWLf9sfGef4lJWpZlYdypHAw/x+5OTAHxWYyQ2JjYsvjy/IOUCGQexutaH4MqrSwz1G0hlm6psa79TLW1LM7KLi+6KIDgmEDtTO6yMi2lbSr4QjkCHyYzrvhvxD5s8f6ShUyNtS9I6MxvOoYlrc6aencD10Ks5jhlZawzlrSvgfX4ySelJWBhZMqb2RP549nuu+ySri9dJr0lIi9eZFcH9yHv0P9aLEuZO7Oq0H0sjK7XYdbZwJSThZZFNKguM1t+IIRCOQGdJk6fx2cnB/PXyIqtaraN1aU9tS9IJDKQGbGi7JSM89MQAwhLD3hpjJDNicdNveRIbzPfXlwMwpNpwHMwcWXR5XoF+WD3JihjS/orgedwz+hzphpHMmD1eB3Ewc1CbbRcLF1LkKbxOeq02m/pERuioboYH5wbhCHQQhVLB+D++xC/4OIuaLqV7+V7alqRT2JjYsqX9DqJTohju92mOh8dNXJvRvXwvvr+xnMDox5gZmjGh7mT+enmRc8/PFJjW//oQaNcRRCRF0OdINxLSEtjd6Re1r1Ccs5LKit72UGJaIiEJL/X2fACEI9A5lEolX1/wZu/Dn/GuP5Oh1UZoW5JOUt2uBsuaf8/fIX8y+8/pOY6Z+9ECjGTGfHV+CkqlkgFVBuNi4criAlwVZPYhKGlVqkDmy4n4tHj6H+vJs7inbO+wWyMd64pyy0p9rjGUSeGqS6CnDF18Ouu/Hxrt5qHJLj6v+SUT6k7Roirdp0eF3viH32St/ypq2temb6X+2e47mpfgq/ozmX4h4/DYq2xXJnlMY+KZMZx6coK2bu01rjE4Jggnc2dMDUw1PldOpMpTGXpiADfDb7C13U4aOn+kkXmK8oogK2JIT3MIQKwItM6bTiDI8CgPTXbhmtqC5+dbFZoicppkVqNvaOLSjClnx+cYHjq42nCq29Vk5gVv4lPj6FOxH25WZVh8eQEKpULj+rQZOqpQKhh7+gvOPDvNsubf065MB43NZWdqh7HMuIiuCPQ7dBSEI9AZXhr8yR3TjdikVaNG8mgk4keTKwykBqxruwV7UweGnBiQVYH0zfu+zZbxKiGEJVcWYygzZEq9r/jn9S2OBR7WuL7g2CCthI4qlUrGnxjPgUf7mNlwLv0qD9TofBKJBCdz56K5IogOwN7UQa/LvYhPGx3BUl4aA6UZCbLnJEnejoQRvBs7Uzu2tt9BZHIEI04Oyqo1lEldx3oMqDKY9bfWcDfiDt3L96JC8Yr4XF6Q1dBGEySlJ/EqIUQrK4IV15by/eXv+bzml4ypPb5A5nSxcC2aK4LYQL3eFgLhCHQGS6ULjRN8UCDnkvkckiWR2pakV9Swr8W3zb/jz5cXmPPn2zVzZjSchbWxNVPPTkAikTCt/gweRj3gwKO9GtP0LPYpUPChoz/d3cqiy/MYUGMAcz9aUGBbjEU1uzgwOkCvt4VAOAKdwlJRkvqJX5MiieGS2Vyik4t2o4+80qtiXz6vMYoNt9ey+/7ObPdsTGyZ1Wgel1/9zZ4Hu+jo3plqdjVYcmXRWysIdREcmxFNUpArgqMBh5lydjytSrVhc+fNBVrfyOXfpDJNrrJ0jfi0eEITX+l1DgEIR6B1NntnbwBSXFEBj8SvSDZ8Sf9fe5OQlqAlZfrJrEbzaOzchClnx2drag/Qp1I/6pdoyNw/ZxKTEs20+jMIjg1i94Od77CmGpmhowWVQ3DxxXm+ODWU2g512ej5I4YywwKZNxNnCxfkSjlh76gQWxgpDKGjIByBTrDZu2W2f0emTGBtm01cC73CML/3t2oUZMdQZsj6tluxNbVj8In+2TJdpRIpPk2XEZ0SzYK/v6Ft6XbUcajLt1d9SJGnqF1LcEwQ5oYW2Jnaqd32/3M73J+Bv/alTDF3dnbci7mhucbn/H/+yyUoOgfGQdEZEUNlxBmBQBN4le3KkmYrOP30N8ae/qJAQh0LC/Zm9mxtt4PXSeGM8BtEuiI9615Vu2qMqDGSn+5u4XrYOGrhAwAAHAtJREFUVbwbfM2L+Odsv7tV7ToyQ0c1vUcfGBNAn6PdKWZcjN1ev1DcxEaj872L/3IJis45QWYOgdgaEmiMgVUGM7PhXA482sdX5ycX2YJe+aGmQ22WNlvJxZfnmfvnzGz3ptb7CkfzEkw9O5GPnZvSyLkxy68tJTEtUa0agmODNH5QHJoYSp8j3VAo5ezxOpjVNlIbFMXs4sCYABzNSmBhaKFtKSohHIGOM6b2+KzSy75abrmob/Sp1I/h1T9n3a017H3wc9Z1CyNL5jdezO3X/my7uwnv+jMJSwxl651NaptboVTwNPaJRg+KY1Ni6HukO+GJ4ezsuI/yxStobK7cUMzYGjMD8yKVSxAYrb8N699EOAIdRyKRMLvRPPpVGsi3V33Y8G8bRkHumPvRQho5N2bSmbHcDvfPuu5VtivNS7Zk4aV5uFuXo3nJlnx/fRnxqXFqmTcsMZRkebLGksmS05MZeLwvD6Pus6Xdduo4emhknrwgkUj+7UtQtFYE+n5QDMIR6AUSiYSlzVfSoYwXMy5My/btVvB+DGWGbGi7DRsTWwaf6E9EUgSQ8Z4ubrKUlPRkZl+cjnf9mUQkR7Dh1lq1zJvZw1YTK4J0RTqfnxrKXy8v8n2rtbQo1Urtc+SXjFyCorEiiEuN5XVSuF73IchEOAI9wUBqwNo2m2ji0oyxp0dyMvi4tiXpDQ5mDmxpt52wxFA+Ozk46/DY3bocY+pM4MCjvSSkJdDOrQOrb36nlvyN4Kw+BG4q23oTpVLJ1LMTOB50lAUf++hciXIXC1deJrzUtowCobCEjoJwBHqFiYEJ29rvpLpdDYb7DeKvlxe1LUlvqO1YlyXNVnD+xVm++WtW1vWxdSZS2sqNaecmMtFjKrGpMaz1X6XyfMGxQUglUlwt1Vt+etGleWy/t42JdacwosZItdpWB84WLoQmvNJYkp4uERhdOCKGQEVHEBUVxYgRI/D09MTLy4vRo0cTGZlRGuHmzZt07twZT09Phg4dSkREhFoEF3UsjCzZ2Wk/rpYlGfBrH26/vqVtSXpD30r9GVptBGv9V7H/4R4ATA1MWdxkKY+jH3H22R90Kduddbd+ULnT1pOYYFwsXDGSGalDOgDr/dew4vpSBlYZwrT6Mz/8gBZwsXBFiZJXiSHalqJxCkvoKKjoCCQSCcOHD8fPz48jR45QsmRJli5dilKpZMqUKcyaNQs/Pz88PDxYunSpujQXeexM7djrdQgrIyv6HOlGYPRjbUvSG75pvIgGTo2YeGZMlhNtVbotHd07s+yaLwOqDCIpPZFVN1aoNI+6Q0f3PdzNzIvedHTvjG/TZTpboty5CIWQBsYE4GTujJmhmbalqIxKjsDa2poGDRpk/X+tWrV4+fIlt2/fxtjYGA+PjEiGvn37cuLECdWUCrLhYunKHq+DKJUKeh3pSkh80diXVRUjmRGbPH/C2rg4Q473JzI5Y6U6v/FiJEjZdHsdPSv0YfPt9YQmvMr3POrsQ3D66SnGnh5JY+cm/NB6IzKpTC12NUGmIygKB8aFodhcJmo7I1AoFOzatYuWLVsSEhKCs7Nz1j0bGxsUCgXR0dHqmk4AlC9egV2d9hOZHEnvI12JShYVS3ODg5kDm9v9xKuEEEacHEK6Ih0XS1em1PuKE8G/Us+xAenKdFZcz98qNj41jtdJ4WpxBNdCrzD0xEAq21Tlxw67MDEwUdmmJlmwPmO7ZMWRswxdfDpb46XCRlBM4cghADW2qpw3bx5mZmYMGDCAU6dOqcWmrW3O2Xr29rrdAKIg9bWxb8YRs8O039GeT/368Nunv2Fh9OEsx6L+Hrazb8kPaT8w7PAwlt9ahG8bX2a0msq+x7tY5b+cT2t8yo+3tvB1y+mUti6dJ40hoRkRQzVKVlHpddwLv0f/X3vhZOnEqUF+OFo45tlGQf6cvSYdwgBTDJRmJEv/O2MZuvg0R77tkuMzuv57CDlrjE6OJiI5gurOqv2M1YE65leLI/Dx8eHJkyesXbsWqVSK0//au/e4qMr8geOfYWC4DCLXgRk0EW9hpHkpzRTTVFRAvKxmJmmtu20Xt8uaa+lP10utaKtuRVtbWa26mqYi4L017eKWmlq65hXBZEaG+/06c35/ACPEMAoCMzDP+/Xy5TByznznnK/P95znec45ajVa7c2uiuzsbGQyGZ6eno1ab1ZWIUZj3dsq+Pl1ICOjeS76aQnWiC/UfSDvj/6Yp/bPJHLDBDZGbMVZ7tzg74ttWCWq01Rm33OU1UdX010ZwqQev+GvD/2NCfFjkVU6IkPGogNLWDPi7UbFeCrlfwB44d/k75FWcJ2IHaNxlDmxefwOHErcyChp3LqstZ+djZ4Uyep2q5mLw9bzEBqO8bS+6uJElWMnq36HX8fn4CBr8ADakjvuGlq7di1nz54lLi4OhaJqhkRoaCilpaWcOHECgC1btjBuXMs/KNyejQ+OZO3D73Dk+pc898Xv7eqe8HdixdBYHggYzItfPsfZzDMM1gzh0V4z+Ne5j4nuPpnN5zeaZofcrpvXEDRtsDi7NItpiRMpqChgS+QOqzzqsjEkSeJ89s9cVHzGEeWLFMm15MtTrB1Wi2oPD6yv7Y7OCC5dusR7771HUFAQ06dPB6BTp07ExcWxatUqlixZQllZGYGBgaxevbpZAhYa9ljITHLKcvjL0YV0dPbkzeHrbHZ2ia1QyBV8NHYDo7eFMXvf4xz8zWEWP7icfSl7uJafikKu4M3jK3l31Ae3vc7U/Kt0dPbE08Wr0fEUVRTx+O6pXCtIZWtUPKG+9zZ6Ha1BkiTOZP5I0pUEkpJ3cTn3EjjL8Dbczd2lMagrHrJ2iC2q5hoCazyGtCXcUSHo0aMHFy5cMPtv/fv3JzEx8U5WLzTBs/fNJbski7dOrcHHxYfXBi++9UJ2zt/Nn/XhG4iOH8fTB59ic8R2Fg5awvyvXmJM0Fi2X9zKC/3/RC/vu29rfSl5TZs6WmGo4Lf7YzilP8n68I08qLGtxtQoGfkh/ThJVxLYfTWRa/kpyGVyhgQO43d9nuHgbk9cJOvcAru1JeddIdC9E66OrtYOpVk022CxYDsWDl5CTlkO606+iZeLN8/c97y1Q7J5AwMeYGXY3/jT4T/yxvfLeG3QYjaf38DJ9B9wc1Ky6vgbfBT+r9taV2p+CqG+fRr1+UbJyB8PPcOha1+w9uF3GB8c2ZSv0ewMRgPf6Y6SlLyLPclJ6Iq0ODk4EdbpYV7qP4+xXSPwcfUB4MlQzM4S+vVT+NqDq+3kZnM1RCFoh2QyGavC1pBblsOSo6/h6ezJYyEzrR2WzYvpPZsf9ad5+9Ra+vj1ZVXYWsZ8/jD9VP1JvBLPmYwfudevr8V1GIwGfim4RkTwhNv+XEmSWPztq2y/tJWFg5bweO8n7vSr3JEKQwVfpx1hd3ICe68mkVmSiYvchRF3jWJR8F8YEzSWjs7mJ360x0bfnOTcK0R1m2TtMJqNKATtlNxBzrujPiC/LI+XDj9PR2dPmznKtGWvD4vlXNZZXjj0LHum/IcnQ+fw8dkP6eDUgdhjr7MxYqvF5bVFaVQYKxo1wPvWyTX886d/8HSfZ/lj/5fv9Cs0SWllKYd/OURS8i72p+wlrywXpZM7o7uMITI4mpFdRrf5h680l5zSbHLKctrNQDGIQtCuOcud+XjcJqYmTODpg0+yJXIHDwUOs3ZYNs1Z7sz6sRsYtS2MWXsf4/MJCSRciUchV3AgdR8nbhxjYMADDS5/84H1Qbf1eRvPfcrr3y9lSo9pLH3ojVYd3C+sKORQ6kGSkndxMPUARRWFdHT2JDxoHJHB0TzceaTNX8BmDaYZQ+2oa0jcfbSdc3dyZ1PENoI8uhKzZzo/6k9ZOySbF6BUsz58I9rCNOYfeYklDy5HW5iGu5M7K4+9bnHZmucQ3M5g8Z7kJOYdeYGRd43irZH/wEHW8v8d88vy2HZhC7P2zqD3+mDmHJjFN2lfMbnHb9gSuYP/zb7MO4+8z9iu40URaEDNjKH2VAjEGYEd8HbxYWtUPJE7xzA9aTLf+n+LN5pbL2jHHlAP4o1hq3nlyIvc69uXIZqhnEw/wVfXv+TbtK+Z6Dfe7HKp+Sk4Ojje8tnBR9O+4emDT9JP1Z+PwjfgJHdqia8BQFZJFvuu7iYpeRdfXT9MhbGCAKWax3s/QWRwNIPUD+LoIJqC2/HUykNccP4PKGS8/s8U5KS1i3ERcUZgJ9TuGrZFxSOTyRi9YTRpBe3/pmB36oneTzIzZBZvnVrD2K4RVBgrqm5bfWwFkiSZXSYl7yqdO9xlsWE9m3mGmL3T6eIRxKaIbSidlM0ee3rRDd49/i5TdkUR+kl3Xjr8PJdyLjLn3j+we/JBTj/xM38d9iYPBQ4TReA21cyIKnLQ4Sr5IsepzvttmcgAOxLs2Z3PIncyKSGCaYkTSZi03zT1T6hPJpPx17A3+Tn7f6z8fgWP9prBv89v4Hvdfzlw5QD9Ow6pt0xq/lWL4wMpeVd5NHESHZw68FnkTrxdmm/7Xy/4haTkXSRdSeD4je+RkOju2YO5/V4istsE7vXtKy4wbIAkSRRWFKAvTiejOAN9cTr64nSKZHmkZP6CvjidY8rLlMlyKZVl4WZUWzvkZiUKgZ25168viY8lEr4hnMeSJrMjOgl3he3f+MtanOXOrA/fyOjPh/Ot9msC3NRklWay8NBC9kw8VK9hTcm/SrSqv9l16Yv1TEucSKWxgp3Ruwns0OmO40vOvUxScgJJV3ZxOqNq/Ke3Tyiv3P8qTwycgR+d7brxL6ksqW7c9eiL9VWvS26+rvo5g4zidEoqS+otL5fJ8XNT4eeqQmHsiIcUhBEj6orBVvg2LUcUAjsU1iWMD8I/ZfbeGczaO4NNEdvEwKAFancNH4VvYNKu8fT2vocbxTp+0P3AvpQ9jOsaYfq93NIccstyzQ4U55flMT1pMvridLZHJ9LTu1eTYpEkiQs550m8Ek/SlQR+zq66wV0/VX8WDV5KZHAUwZ7dgbZxU7emKDeUk1mSUd24p5tt4GteF5Tnm12Hr6svfq4q/Nz86doxGJWbP36uKlRuqqrX1X/36tyFrMwioH10ATVEFAI7FR40jr+PfJfn//M0fzj4Wz4M/1T0FVswSD2Y14eu4s9fvUyQR1dS81NY8d8lhAeNM832aWjqaGllKbP2zuB89jk2jt/KAP/7G/XZkiTxU8ZpkpIT2J2cwOXcS8iQcX/AIJY99AYRwRPo3MzPRm5tBqOBrNKsWkfv6ehL9KbXGbUa+OwGnrvR0dkTlasKPzcV9/r2qdW4+9dp4H1cfG97cL41ZnLZAvE/345N6/UYuaU5LPp2Aa8ceZE1D79t190ItzL7nt/yo/4U/z6/AblMzqXciyRc3snEHlOAm4Wg9sVkBqOBZ76Yw7far/nHqA8Zedeo2/qsOvf1SU7gWkFq1X19NEOZ0+cPRHSNwl8Z0OzfsTlJkkRuWQ76Yn31EXq66Qg+o1hPbmUW13O1ZJToySzJwCgZ663DzdHN1IB39+rJkMCh9Rr4mq6blj6rXb9gZLu9hYYoBHbu932fJbssmzUnVuHp7MWSIcutHZLNkslkrAz7G+ezz3E26wwGg4G/HF1EZLdoHB0cSakpBNVnBJIkMf+rl9mdnMDrQ2OZ0nOaxfXXvq/P7uREbhTpbt7XZ8Arde7rYy2SJFFUUVh9xF7Vt27qa681yFrT+FcYK+qtQ+GgQOXmj9ojgE4dOtHff4CpMa9q3G828LZ2NXN7aPTNEYVA4M/3LySnNJu403/H29WHuf1etHZINsvF0YX1YzcyZvtwsoqz0Balsfn8RmJ6zyY1/yq+rr6mwffYYyvYcO5jXuw/j9/1ecbs+soN5XyT9pXZ+/pEBk+weF+f5lRSWVKr+yWj3tF7TVdNZrGe4sriesvLZXJ8Xf2qj979uNunNyrXqtc3G3d//Fz96OjsiUwma7djGG2RKARC1TTJYW+SW5rD8v8uxsvZi5m9Z1k7LJulcQ/k86mfM+LTEQAsPfp/PNprBim1Hlj/4U/vseaH1cwMmcWrg/6vzvKtdV+fCkMFukKt+ZkyxRnoS2429g0Nqvq4+KBy88fXTcX9Hg/cbNBrNfB+riq8XbyRO8jvOGbBOkQhEICqQbG3H3mfvPI85h15gY7OnkR1M/+cWQGGdRnGiqGxvPr1PPLL83j71DpS864yMOB+dlzaxmvfzGd81yhWDV+LTCZrtvv6GIwGskuza019rDVTplhfPcBa9W8NDap6KDqaBk9Dffqg6lx7pozK1Lj7uvq16BXPgu2QSQ1dImkDxDOLW4alGIsqipiWOJEf9afYFLGN4Z1HtHJ0bWcb6vX5XPzdbGTUH2CXkPD++yoOX/+SpOQEvrz2BaWGUnxdfRnXNZKI4AkMDQxDIVcgSRJ5ZblVDXpJusVpkVklmRik+o8hdXN0M015rJkG2dXvLpSSZ70G3lamCreV/WzLMTbXM4vFGYFQh9JJyabxW4mOH8esvTPYEZ1If/+B1g7LJslkMrNFAECGjL7/CqFSqsTbxYchgUPp1rEHSic3Mksy+fjsB8QeW2Hqgy83ltdbh5ODk2ngVKPUcJ9fv+pBVP86f6tcVSid3OvN+LL1RkywHaIQCPV4unixNSqeiJ1jeCxpComTDjT5Aqj2qKSyhMvZ6ZxNu4ifhd+rlCqBqofRH7r2BYf4AgeZg2lQVeWmoqfX3XVmyajc/FG5Vv1cM6gqCC1NFALBLH9lANui4onaGc7UxGiSJh9o8xct3UpNF42uSIeuKA1doQ5tURq6Qi26Ii3aQi03irTklOWYlvmGyQ2ub/GDy6sa+FrTIsWgqmCLRCEQGtS1YzCfRe4kOn4cUxOiSZx0AD83S8fAtssoGckoyUBXmIauSIe28GYDb/pTqDU7NdLPVYXaXUMXjy4M1jyIWqmhl7obSqMXnPykwc98vt8LLfiNBKH5iEIgWHSPbyibIrYxLTGa6UmTiZ+4mw4KD2uHVUe5oZwbRTq0RVpuFGrR1mrYtYVp3CjScaNYR6Wxss5yjg6OqJUaApRqQn36MKpLOBplIBp3DQFKDRp3Df5uASjkinqfWdP/fpFPWulbCkLLEYVAuKVB6sGsD99AzN7pxOyZzubI7bg6urbKZxdWFNbqmqlq1LWFadVH8VWvM0sy6i3n5uiGxj0QtVLDkMChqJUa1O4a1EoNGqUGtXsgvq6+dnMvGUGwRBQC4bY80mUM7zzyPs8cnMMDcVEMKPkzDtzs627spfeSJJFdml195F7dXVPdH59VoSc1+xq6Ih355Xn1lvV28a46Yldq6Ot3X1Xj7h5IgFJd3fir8VB0bJWB1p4fftLinyEILU0UAuG2Te4xlXcSvues6z/5SXqHvqVzkVU/5O6plYdMxaDSWElGsR5tUZppgFVbqz++5si+zFBWZ/0OMgf83QLo7NmJ7l49GdZpOOrqhl2jDETtXtWN01pnI4JgL0QhEBolqGI85bICLrpsplxWgLehN6WybEodMhn7+XJ0RTrSi2/Uu5Oks9zZdMQ+wH8g6l/1xauVGlRu/jg6OIr574LQykQhEBqtR/k08uWp3HA6it7pBI6SGy5GHzyce9DLO6ReX7xaqcHbxVvMiRcEGyUKgdBoMmQMLJlPYWkaLpI3jlR11ayPap+36BWE9k4UAqHJ3KVAa4cgCEIzEHPnhEZpaHZQe31ghyDYA3FGIDSaaPQFoX0RZwSCIAh2ThQCQRAEOycKgSAIgp0ThUAQBMHO2fRgsYOD+QuQGnrfVth6fGD7Mdp6fCBibA62Hh/Yfoy142tqrDb9zGJBEASh5YmuIUEQBDsnCoEgCIKdE4VAEATBzolCIAiCYOdEIRAEQbBzohAIgiDYOVEIBEEQ7JwoBIIgCHZOFAJBEAQ7Z9O3mKjt6tWrLFiwgNzcXDw9PYmNjSUoKMhq8eTk5DB//nyuXbuGQqGgS5cuLFu2DG9vb3r16kXPnj1xcKiqs6tWraJXr15WiXPkyJEoFAqcnZ0BmDdvHsOGDeP06dMsXryYsrIyAgMDWb16NT4+Pq0a2/Xr13nuuedMPxcUFFBYWMixY8cajLs1xMbGsn//ftLS0khMTKRnz56A5Rxs7fw0F6OlnARaNS8b2oaW9mtr56S5GC3l5K3ib26W9qelbdWk7Si1ETExMVJ8fLwkSZIUHx8vxcTEWDWenJwc6bvvvjP9vHLlSunVV1+VJEmSevbsKRUWFlortDpGjBghXbhwoc57RqNRGjVqlHT8+HFJkiQpLi5OWrBggTXCq2PFihXS0qVLJUkyH3drOX78uKTVauvFYCkHWzs/zcVoKSclqXXzsqFt2NB+tUZONhRjbbVzUpJaNy8b2p+WtlVTt2Ob6BrKysri3LlzREZGAhAZGcm5c+fIzs62Wkyenp4MGjTI9PN9992HVqu1WjyNcebMGZydnRk4cCAA06dPZ9++fVaNqby8nMTERKZMmWLVOAAGDhyIWq2u856lHLRGfpqL0ZZy0lx8llgjJ28Vo7VzsqH9aWlbNXU7tomuIZ1Oh7+/P3K5HAC5XI5KpUKn05lOe63JaDSyefNmRo68+QjHmJgYDAYDYWFhzJ07F4VCYbX45s2bhyRJDBgwgJdffhmdTodGozH9u7e3N0aj0dStYQ2HDh3C39+fe+65x/Ter+P28PCwSmxgOQclSbK5/DSXk2AbeWluv7aVnATr5GXt/WlpWzV1O7aJMwJbt3z5ctzc3Jg5cyYAhw8fZseOHWzatInLly8TFxdntdg2bdpEQkIC27dvR5Ikli1bZrVYLNm+fXudI6+2Eret+nVOgm3kZVvar7/OSbBe/Ob2Z3NqE4VArVaTnp6OwWAAwGAwoNfrG3Xq2VJiY2NJTU1l3bp1pkG4mrjc3d2ZOnUqJ0+etFp8NbEoFApmzJjByZMnUavVdboMsrOzkclkVjvySk9P5/jx40RFRZneMxe3NVnKQVvLT3M5WfMdwLp52dB+bQs5CdbJy1/vT0vbqqnbsU0UAh8fH0JCQkhKSgIgKSmJkJAQq3cLrV27lrNnzxIXF2c6xc7Ly6O0tBSAyspK9u/fT0hIiFXiKy4upqCgAABJktizZw8hISGEhoZSWlrKiRMnANiyZQvjxo2zSowAO3fuZPjw4Xh5eVmM25os5aAt5ae5nATbyEtL+9XWcxKsk5fm9qelbdXU7dhmHkxz5coVFixYQH5+Ph4eHsTGxhIcHGy1eC5dukRkZCRBQUG4uLgA0KlTJ+bMmcPixYuRyWRUVlbSr18/XnvtNZRKZavH+MsvvzB37lwMBgNGo5Fu3bqxaNEiVCoVJ0+eZMmSJXWmmPn6+rZ6jADh4eEsXLiQsLCwW8bdGlasWMGBAwfIzMzEy8sLT09Pdu/ebTEHWzs/zcW4bt06szkZFxfHqVOnWjUvzcX33nvvWdyvrZ2TDe1nqJ+T0Pp52VAbExcXZ3FbNWU7tplCIAiCILSMNtE1JAiCILQcUQgEQRDsnCgEgiAIdk4UAkEQBDsnCoEgCIKdE4VAEATBzolCIAiCYOdEIRAEQbBz/w8pYrFLlnd7XAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#model parameters\n",
    "N_PATIENTS = 15\n",
    "VEHICLE_CAPACITY = 4\n",
    "TIME_LIMIT = 10\n",
    "SEED = 101\n",
    "\n",
    "#create a problem instance\n",
    "instance = problem_instance(N_PATIENTS, VEHICLE_CAPACITY, random_state=SEED)\n",
    "\n",
    "#create an instance of a gurobi optimisation model\n",
    "model, x, u = optimisation_model(instance)\n",
    "\n",
    "#limit run time?\n",
    "if TIME_LIMIT > -1:\n",
    "    model.Params.TimeLimit = TIME_LIMIT  \n",
    "\n",
    "#execute the optimisation procedure\n",
    "model.optimize()\n",
    "\n",
    "#show the routes\n",
    "plot_solution(x, instance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}