demo_India_demographics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Effect of Demographics and Population Clusters on Lock-down Outcomes\n",
    "### Accompanying paper available at https://arxiv.org/abs/2005.11257"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from skopt import gp_minimize as gpm\n",
    "import numpy as np\n",
    "import time\n",
    "from utils import (Population, Demographics)\n",
    "from matplotlib import pyplot as plt\n",
    "import pylab as plb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Getting started**: we load the generic population file, replace the demographics (RST and SUS values) based on distributions aligned to Indian statistics, reorganize portion of the population into urban clusters, and set the simulations to run for 500 time steps.\n",
    "\n",
    "**Notation**: We use I0, P0, L0 to denote respectively, the initiation point, period and the level of the lockdown that ESOP is trying to optimize. These are referred ensemble as IPL. In this experiment, all three parameters are being optimized under constraints that I0 must be no less than 12 (i.e. no lock-down before t = 12) and P0 must be no larger than 40 (i.e. no lock-down longer than 40 days)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pop = Population.fileInit( \"pop_generic\" )\n",
    "N = pop.N\n",
    "\n",
    "# Load the patch file\n",
    "demo = Demographics.fileInit( \"patch_India_fit\" )\n",
    "\n",
    "# Use city clusters as initial locations of the individuals\n",
    "pop.X = demo.loc[:, 0]\n",
    "pop.Y = demo.loc[:, 1]\n",
    "\n",
    "# Use RST and SUS values derived from Indian demographics (see accompaying paper)\n",
    "newSUS = np.sum( demo.demoData[:,1:], axis = 1 ) * 0.2 + 0.1\n",
    "newRST = 1 - newSUS\n",
    "pop.SUS = newSUS\n",
    "pop.RST = newRST\n",
    "\n",
    "pop.reset()\n",
    "T = 500"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Age structure in Indian population**: The following code is meant to reproduce Fig 6(a) in the accompanying paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX9//HXRzaVBNk0iihUZbHYak2sWqUNarV1bdXWHaiV4FJtS+3mRkRbtbX6tbUqoVJ3UL9uX/CHtVbSutWvoOhXlFA3VLayGQmbBT+/P86dMBkzk5thMjMh7+fjcR+Ze88995w7k+Qz99xzzzF3R0REpNhsU+gKiIiINEcBSkREipIClIiIFCUFKBERKUoKUCIiUpQUoEREpCgpQImISFFSgBIRkaKkACUiIkWpc6ErUEh9+/b1gQMHZp1/zZo1dO/ePXcVUpkqsw3Nnj27yXp5eXleyoWt/71Vma0ze/bs5e6+Y4s7unuHXcrLy31LzJw5c4vyq0yVmc8yJ06c6BMnTvRx48b5xIkT81au+9b/3qrM1gFmeYz/0R36CkqkI6mqqgKgtraWysrKwlZGJAbdgxIRkaKkACUiIkVJAUpERIqSApSIiBQlBSgRESlK6sUn0kGYWZN112zaUuR0BSUiIkVJAUpERIqSmvi2wOKGxVTXVue1zEoq81qeiEih6ApKpINIDB8zc+ZM3X+SdkEBSkREipIClIiIFCUFKBERKUoKUCIiUpQUoEREpCipm7lIB1FTUwNAXV0d8+fPb5wfSqRY5TVAmVlDyqbtgFvc/cIo/XDgj8DuwIvAaHdfEKV1A24FTgbWAr9x9xuSjp02r4jA2LFjm6wrQEmxy2sTn7uXJBagDFgHPAhgZn2Bh4HLgd7ALOD+pOzVwCBgADAC+JmZfSNmXhERaWcKeQ/qZODfwDPR+onAXHd/0N3XEwLSvmY2NEofCVzl7qvc/U1gEjA6Zl4REWlnrFBPlJvZ08A/3L06Wr8J6Oru5yXt8zowHngaWAns7O5Lo7STgfHu/oVMed39oZRyq4AqgLKysvKpU6dmfQ4r61eyodOGrPNno5RSSkpK8lpmQ0ODytwKyrz++usB2LhxI507d+biiy/OS7mw9b+3KrN1RowYMdvdK1raryCdJMxsd+BrwPeTNpcAy1J2rQdKo7TEempaS3mbcPcaoAagoqLCKysrW38CkSnTp1BXUpd1/mxUUsmW1DkbtbW1KnMrKDNRztZ+nipz6ymzUE18I4Fn3f3dpG0NQI+U/XoAq6M0UtITaS3lFRGRdihWgDKzvtFVT/K2sWb2BzM7NotyRwJ3pmybC+ybdPzuwJ6Ee0urgMXJ6dHruS3lzaJuIiJSBOJeQU0GfpFYMbPLCV2+TwceM7NT4hZoZl8BdiXqvZfkEWAfMzvJzLYFrgBec/d5UfpdwGVm1ivq/DAGuCNmXhERaWfiBqgK4G9J6+cCv3b3PoRnj8a1osxRwMPu3qT5zd2XAScBvwJWAQcCpybtMh54G1gA/B34rbs/ETOviIi0M3E7SfQGEr3n9gF2ZnMT3aOEJrtY3H1shrSngGa7hrv7BuDsaGlVXhERaX/iBqgVQP/o9WHAInf/V7TeBY3pJ1L0ysvLAVi9ejWlpaXMnj27wDUSySxugHoKqI5GbPgJ4aopYSih2U1EitjLL79c6CqItErcK5+fAR8A1xDuA12ZlHYG8GyO6yUiIh1crCuoaPSGr6dJPgJYn7MaiYiIkIORJNz941xURETa1qxZsxp/VlS0OMqMSMGlDVBmNrkVx3F3/37Lu4lIoSR3kki8Filmma6gDgPijiRbmBFnRURkq5U2QLn7wDzWQ0REpAk9vyQiIkWp1Z0kzGwnYNvU7e7+fk5qJCIiQswAZWbbAFcDY4GeaXbrlKtKiUjuJUaOqKuro7S0VB0lpOjFvYL6EXABcB0hUP0K+JTwkO6nwLVtUjsRyZnUruWFmk1bJK6496C+B0wgBCiAR9x9PLA3sBDYPV1GERGRbMQNUHsAs9x9E7AR2A7A3f8D/BdpRhgXERHJVtwmvno2d4xYBAwBnks6Ru8c10tEcmz//fcHNo9mLlLs4gaoV4DPA3+JlivNbB3haupXgIZJFilyiU4StbW1VFZWFrYyIjHEDVD/RWjmgzCz7f7AvdH6AuAHOa6XiIh0cHFHM/9r0uslZvZlYE9ge+DN6F6UiIhIzmQ1mrmH/qlv5bguIiIijWL14jOzr7a0tKZQMzvVzN40szVm9raZDY+2H25m88xsrZnNNLMBSXm6mdlkM/vYzJaY2biUY6bNKyIi7U/cK6haWh6xPNZIEmb2dcLzVKcA/wvsEm3vCzwMnANMA64C7gcOirJWA4OAAcDOwEwze8Pdn4iRV0RE2pm4AWpEM9v6AMcCX6N1nSSuBCa4+z+j9YUAZlYFzHX3B6P1amC5mQ1193nASOB77r4KWGVmk4DRwBPAiS3kFenwqqqqAFi0aBH33XcfNTU1Ba6RSGa2pcOdmNmNQDd3Pz/Gvp2AdcAVhKudbYFHgZ8Shkvq6u7nJe3/OqHX4NPASmDnaPp5zOxkYLy7f8HMbkqX190fSqlDFVAFUFZWVj516tSsz31l/Uo2dNqQdf5slFJKSUlJXstsaGhQmVtBmSNGNP2eOXPmzLyUC1v/e6syW2fEiBGz3b3FaZ23eMp34HFgKtBigALKgC7AycBw4D/AY8BlQAmwLGX/eqA0Skusp6bRQt4m3L0GqAGoqKjwLXkeZMr0KdSV1GWdPxuVVOb9GZZCPDejMttePsvtKO+tysytXMwHNYQwYGwc66Kff3D3xe6+HLgBOBpoAHqk7N8DWB2lkZKeSKOFvCIi0g7FnW5jZDObuwL7AN8ndFBokbuvMrMPab7DxVxgVFKZ3QnPWs2N8i0G9gUSz2TtG+XJmDdOvUQ6gokTJwJhuo0hQ4YUuDYiLYvbxHdHmu0bCL3lftiKMv8MXGhmTxCa+H4ETAceAX5rZicRmg2vAF5L6uRwF3CZmc0iNBWOIYyyToy8Ih1eopOEhjqS9iJugPpcM9vWJzostNJVQF9gPrAeeAD4lbuvjwLMzcA9wIvAqUn5xgO3EoZWWgdc5+5PALj7shbyiohIOxN3qKMFuSowGhbpfJrpVOHuTwFD0+TbQJjWo9mpPTLlFRGR9icXnSRERERyLu0VlJl9SsujRzRy91gjSYiIiMSRqYlvApsDlBGa1rYjDCW0lDDc0LGE+0G3t2EdRUSkA0oboNy9OvHazC4jdE44yt3XJm3vTpjAcGMb1lFEcsDMmqxv6SgyIm0t7j2oscBvk4MTgLuvAa4Hzs11xUREpGOLG6D6Eh7MbU5XwsCxIiIiORM3QM0CrjSzXZM3RuvVwEs5rpeIiHRwcQPURUA/4G0zqzWz+82sFnib0FmiNSNJiEgBuDvuzsyZM3X/SdqFWAHK3V8B9gJ+B2wCvhD9vB4Y5O5z2qyGIiLSIcWebsPdVwCXtmFdREREGmkkCRERKUqZRpJ4Gjjf3edFrzNxdz88t1UTEZGOLFMTX/JTfduQedgjy5AmIiLSaplGkhiR9LoyL7URkTZTU1MDhAkL58+f3zg/lEixijujbp+ok4SItFNjx45tsq4AJcUubieJxWb2qJmdZGbpRpQQERHJmbgB6nJgD+BBYImZ3WZmh7RdtUREpKOLO6PudcB1ZrYfcBZhOvUxZvYecDdwr7v/q81qKSJbbMyYMQAsWrSIfv36Fbg2Ii1r1XNQ7j7H3X8C7AZ8E3gO+AnwZtxjREMlrTezhmipS0o73cwWmNmaqEmxd1JabzN7JEpbYGanpxw3bV4RCZ0kampquPjiixs7TIgUs6we1HX3T4E1hMkK/0Pru5n/wN1LomUIgJkNAyYSrtDKgLXALUl5/gh8EqWdAdwa5YmTV0RE2pnYQx0BmNkgQhA4AxgILCQEhrtzUJczgGnu/o+orMuBN82sFPgUOAnYx90bgGfN7H+iuvwiU153X52DuomISJ7FuoIyswvM7J/APOBHwD+ArwMD3P2X7v5GK8u9xsyWm9lzZlYZbRsGvJrYwd3fJlwxDY6WTe4+P+kYr0Z5WsorIiLtkMUZdt/M/gM8RbhSesTd12VdoNmBwBuEAHIqcDOwH1ADPOjutyXtu5BwdbQpSts5KW0McIa7V5rZ39LldffalPKrgCqAsrKy8qlTp2Z7KqysX8mGThuyzp+NUkopKSnJa5kNDQ0qU2WqTJWZszJHjBgx290rWtovbhPfbu6+JOvaJHH3F5NW7zSz04CjgQagR8ruPYDVhCa+dGm0kDe1/BpCMKSiosIrKytbfxKRKdOnUFdS1/KOOVRJJVtS52zU1taqTJWpMlVm3suM2808J8Ep3eEJnSzmAvsmNprZHkA3YD4hQHU2s0FJ3dn3jfLQQl4RAcrLywFYvXo1paWlzJ49u8A1Esksr6OZm1lP4EDg78BG4BTgq4T7Wp2BF8xsOPAyMAF4ONHJwcweBiaY2TmEJsETgK9Eh743U14RgZdffrnQVRBplXyPZt4FuBoYSrivNA/4lrvXAZjZuYRg04dwz+t7SXnPByYD/wZWAOe5+1wAd5/bQl4REWln8jqaubsvAw7IkH4fcF+atJXAt7LJKyIi7U+rnoMSkfZr1qxZjT8rKlrsQCVScJnuQX21NQdKPCQrIsUpuZNE4rVIMct0BVXL5vtORuZ7UACdclEhERERyBygRiS97gn8AXgdmAosJYx5dxphFIcL2qqCIiLSMWXqJPH3xGszuwN40t3PSdntLjO7HTgRmNYmNRQRkQ4p7mjmJwD3p0m7P0oXERHJmbi9+LYB9gL+2kzaIHT/SaToJUaOqKuro7S0VB0lpOjFDVCPE41AThihYZOZdSJMgXE1ML2tKigiuZHatTzOQNEihRQ3QF1EmEX3fmCjma0CekX5n43SRUREcibuYLHLgeFm9nXgIGAXYDHwgrs/1Yb1ExGRDqpVI0m4+19p/j6UiBS5/fffH9g8mrlIsWv1UEdmthOwbep2d38/JzUSkTaR6CRRiPmDRLIRK0CZWQ/gJsL0GN3S7KaefCIikjNxn4P6IyE43Q6cB5zdzCIikrV7772XI488Muv83/zmN7nzzjtzWKPm1dbW0r9//zYvB2DgwIE89VR2t/nNjLfeeqvZtDvuuINDDz10S6qWF3ED1FHAT939Qnevcfc7U5e2rKSIFI9nn32Wr3zlK+ywww707t2bQw45hJdeemmLj3vGGWfw5JNPZp1/xowZjBo1Cgj/gC+88MKsjvPNb36TkpISSkpK6NKlC127dm1cP/fcc7Oun7Re3HtQBtS1ZUVEpPh9/PHHHHvssdx6661897vf5ZNPPuGZZ56hW7d0Lf9tz91xd7bZJu737cxmzJjR+Hr06NH079+fq6++unFbbW1tq463adMmOnXSHZBsxP1EpwLHtWVFRKT4zZ8/H4DTTjuNTp06sd1223HkkUfyxS9+sXGfyZMns/fee9OrVy+OOuooFixY0JhmZtx2220MGjSIXr16ccEFFzQ+MJza7PT8889zwAEHsMMOO3DAAQfw/PPPN6ZVVlZy6aWXcsghh7D99tvzzjvvUFlZyZ/+9CfefPNNzj33XN544w1KSkro2bMnL730EmVlZWzcuLHxGA899BD77bdf1u/F7373O3baaSd22WUX/vznPzduHz16NOeddx5HH3003bt3Z+bMmWzYsIGLL76Y3XffnbKyMs4991zWrVsHwPLlyzn22GPp2bMnvXv3Zvjw4Xz66aeNx5szZw5f/OIX2WGHHTjllFNYv359Y9qkSZM444wz6N27N8cffzyLFi1qtq4rVqzg+OOPp0ePHnz5y1/m7bffzvq88ylugHoSONbMJpvZyWZ2WOrSlpUUkS1XVVVFVVUV119/PVVVVVkdY/DgwXTq1IlRo0YxY8YMVq1a1ST90Ucf5de//jUPP/wwy5YtY/jw4Zx22mlN9pk+fTovvfQSr776Kg888AB/+ctfPlPOypUrOeaYY7joootYsWIF48aN45hjjmHFihWN+9x9993U1NSwevVqBgwY0Lh977335rbbbuPzn/88DQ0NfPTRRxxwwAH06dOHv/5181My99xzD2eddVZW78OSJUuor69n4cKF3H777VxwwQVN3ov77ruPSy+9lNWrV3PooYfy85//nPnz5zNnzhzeeustFi5cyIQJE4AQ6Pr378+yZctYunQpv/71rzGzxmM98MADPPHEE7z77ru89tpr3HHHHQA8/fTT/PKXv2T8+PEsXryYAQMGcOqppzZb3wsuuIBtt92WxYsXM3nyZCZPnpzVeedb3AD1GPA5YDTwAPBUtPw16aeIFLFJkyYxadIkHn/8cSZNmpTVMXr06MGzzz6LmTFmzBh23HFHjj/+eJYuXQrAxIkT+eUvf8nee+9N586dueSSS5gzZw5LlixpPMYvfvELevbsye67786IESOYM2fOZ8p5/PHHGTRoEGeddRadO3fmtNNOY+jQoUybtnnShNGjRzNs2DA6d+5Mly5dWqz7qFGjuOeee4AQAP/yl79w+umnZ/U+dOnShSuuuIIuXbpw9NFHU1JSQl3d5rsgJ5xwAocccgjbbLMN3bp1Y9KkSdx444307t2b0tJSLrnkEqZOndp4rMWLF7NgwQK6dOnC8OHDmwSoiy66iH79+tG7d2+OO+64xvfr3nvv5eyzz2bw4MF069aNa665hhdeeIH33nuvSV03bdrEQw89xIQJE+jevTv77LNP4726Yhc3QI1IsxyW9LNVzGyQma03s3uStp1uZgvMbI2ZPWpmvZPSepvZI1HaAjM7PeV4afOKSO7svffe3HHHHXz44Ye8/vrrLFq0iB/96EcALFiwgB/+8If07NmzscnK3Vm+fHlj/p133rnx9fbbb09DQ8Nnyli0aFGTqyKAAQMGsHDhwsb13XbbrVX1PvPMM5k2bRoNDQ088MADDB8+nF122aVVx0jo06cPnTtvvoWfeh7JdVu2bBlr166lvLy88X35xje+wbJlywD46U9/yl577cWRRx7JHnvswbXXXtukrHTvV+p7VFJSQp8+fZq8R4nyN27c2KROqe9tsYoVoNz97y0tWZT9R6Cx64+ZDQMmAmcRJkNcC9ySsv8nUdoZwK1Rnjh5RaQNDB06lNGjR/P6668D4R/zxIkT+eijjxqXdevWsc8++7TquP369Wty7wrg/fffZ9ddd21cT77KSNVc2q677srBBx/MI488wt133511814cyeX37duX7bbbjrlz5za+J/X19Y2BprS0lN/97ne88847TJs2jRtuuIG//e1vLZaR+h6tWbOGFStWNHmPAHbccUc6d+7MBx980Ljt/ffbx7gKrR5JIhfM7FTgI+B5wjQeEILONHf/R7TP5cCbZlYKfEoYOX0fd28AnjWz/yEEpF9kyuvuq/N4am1uccNiqmur81pmJZV5LU/axsSJE4Ew3caQIUOyOsa8efN4/PHHOeWUU+jfvz8ffPABU6ZM4aCDDgLg3HPP5fLLL2e//fZj2LBh1NfX8+STT7Ljjju2qpyjjz6aCy+8kPvuu4/vfve7PPTQQ7zxxhsce+yxsfKXlZWxbNkyPvnkE7p27dq4feTIkVx77bUsWLCAb3/7262qU7a22WYbxowZw49//GNuvvlmdtppJxYuXMjrr7/OUUcdxfTp0xk6dCh77rknPXr0oFOnTrF6/Z1++umceuqpDB48mIMPPphLLrmEAw88kIEDBzbZr1OnTpx44olUV1czefJk3nvvPe68887P7FeM0gYoM7urFcdxd4/VqBmNSjEBOBz4flLSMELAShzwbTP7BBhMCFCb3H1+0v6vAl+LkXd2SvlVQBWEX+LWdhlN1m1TN4Y0ZPeH3p7KbKBhi96nrMpsUJm5NnjwYCB88y4pKcmq3GXLljFt2jSuvfZaGhoaKCkp4aCDDuJb3/oWtbW19OrVi+OPP57jjjuOpUuX0r17d8rLy7ngggsay3vxxRf58MMPgdDZYNOmTdTW1jJv3jzq6+sb97vqqquorq6mqqqKfv36cdVVVzVeqX300UfU1dU1OYfkbZ06dWK33Xajb9++mBmPPfYYEJrm3nnnHQ499NBYz24l1y9hzpw5bNiwocm29evX8+qrrzJkyJBm8xxzzDHceeed7LffftTX19O3b1+OP/54unXrxowZMzjnnHOor6+npKSE444LHaZra2sbj5toTnzvvfdYunRp4zmeddZZXH755YwbN45hw4Yxbty4JuUm3uvTTjuN6667jh133JHdd9+dww47jFdeeSXr3718/d5aujlhzOw9IO6EMe7ue8Qq0OwmYJG7X2dm1cBe7n6mmf0NeNDdb0vadyHh6mhTlLZzUtoY4Ax3r8yU191r09WloqLCZ82aFfMUP2vK9CnUleT38bAhDUPyXmYllXkfu60Q48WpzI5R5p577snEiRM54ogj8lZmW2qPZZrZbHevaGm/tFdQ7j4w69LTMLP9gCOALzWT3AD0SNnWA1hNuIJKl9ZSXhERIDz7ZGYcdpiejGkP8n0PqhIYCLwf3UQsATqZ2eeBJ4B9Ezua2R6EgWnnEwJUZzMb5O7/inbZF5gbvZ6bIa+ICJWVlbzxxhvcfffdORt1QtpWvgNUDWFUioSLCQHrPGAn4AUzGw68TLhP9XCik4OZPQxMMLNzgP2AE4CvRMe5N1NeEZF831+ULZfXrxHuvtbdlyQWQtPcendf5u5zgXMJwebfQClwflL284HtorQpwHlRHmLkFRGRdqYg3cwT3L06Zf0+4L40+64EvpXhWGnzishnnw1K10FKpFioIVZERIpS2gBlZg+b2V7R65Fm1id/1RIRkY4u0xXUCUBiPLs/A3u2fXVERESCTAFqKXBw9NqI/9CuiBShxMR+M2fO1P0naRcyBagHgBvNbBMhOP3TzDalWTZmOI6IiEirZerF92PgOeDzwHjgDmBhhv1FRERyJtNQRw48CGBmo4Gb3P3VPNVLREQ6uFjPQbn759q6IiIiIsliPwdlZruY2fVm9pKZvW1m/2tmvzGznVvOLSIi0jqxrqDMbDDwLNCTcF/qLWBn4IfASDMbnjSIq4gUoZqaGiBMWDh//nyqqqoKXCORzOIOdXQdUA982d3fS2w0swHAk1H6iTmvnYjkzNixY5usK0BJsYvbxDcCuDw5OAG4+wKgOkoXERHJmbgBqivpJ/9bHaWLiIjkTNwmvjnAhWY2w90/TWy0MDzy+VG6iBSxMWPGALBo0SL69etX4NqItCxugJoATAfeNLP7gcWEThLfAQYBx7RN9UQkVxKdJGpra6msrCxsZURiiPsc1BNmdixwNXApm8fmmw0c6+5Ptl0VRUSkI4o9YaG7PwE8YWbbA72AVe6+ts1qJiIiHVqrZ9SNgpICk4iItCnNqCsiIkUp7wHKzO4xs8Vm9rGZzTezc5LSDjezeWa21sxmRg8CJ9K6mdnkKN8SMxuXcty0eUVEpP1pdRNfDlwDfN/dN5jZUKDWzF4BFgAPA+cA04CrgPuBg6J81YQegwMIPQhnmtkbUQeOvi3kFenwysvLAVi9ejWlpaXMnj27wDUSySzvAcrd5yavRsueQDkw190TU3xUA8vNbKi7zwNGAt9z91XAKjObBIwGniAMs5Qpr2yBxQ2Lqa6tzmuZlVTmtbyO4OWXXy50FURaxVqa+tnMugL/BH6Rq+7kZnYLIbhsB7wCfBX4FdDV3c9L2u91wmSJTwMrgZ3dfWmUdjIw3t2/YGY3pcvr7g+llF0FVAGUlZWVT506NevzWFm/kg2dNmSdPxvdNnXrEGWWUkpJSUley2xoaNiqyxwxoumIZDNnzsxLubD1v7cqs3VGjBgx290rWtqvxSsod//EzD4H5Gxad3c/38wuBA4GKoENQAmwLGXXeqA0Skusp6bRQt7UsmuAGoCKigrfkgcWp0yfQl1JXdb5szGkYUiHKLOSyrw/TFqIB1gL+dBsPsvtKO+tysytuE18fwWOJFzJ5IS7bwKeNbMzgfOABqBHym49CGP9NSStr09Jo4W8IgLMmjWr8WdFRYtfXkUKLm6A+gNwj5l1Bh4lDHXUpG3Q3d/ZgjrsCcwFRiU2mln3xHZ3X2Vmi4F9CcGS6HXiflbavFnWSWSrk9xJIvFapJjF7Wb+d2A3YFz0ej7wr5SlRWa2k5mdamYlZtbJzI4CTiNcmT0C7GNmJ5nZtsAVwGtJnRzuAi4zs15R778xwB1RWkt5RUSknYl7BfW9HJXnhOa82wjBcQHwI3d/DMDMTgJuBu4BXgROTco7Hrg1yrMOuC4afgl3X9ZCXhERaWfiDhZ7Zy4Kc/dlwNcypD8FDE2TtgE4O1palVdERNqfVo0kYWbbmNk+Zva16D6PiIhIm4gdoMzsAmAJ8BrhntGQaPujZnZR21RPRHJl9uzZzJ49m7q6Oo0iIe1CrCY+MxsD3ARMBp4EHkhKfgY4Cfh9zmsnIjmT2rW8pYf0RQot7hXUOOB37l5F6DGXbB7R1ZSIiEiuxA1QnwP+kiZtDdAzN9UREREJ4nYzXw4MTJM2BFiYk9qISJvZf//9gc2jmYsUu7gBahpwhZnVEp5DAvBomosfE0aXEJEilugYUcjx/0RaI24T32WEAV1fB54iPHD7e+BNYBMwoU1qJyIiHVasAOXuK4AKwmSDXYC3CVdfNwMHu3t9huwiIiKtFnvCQndfTZip9qq2q46IiEjQqhl1zawHsA+wK/AhYaTxj9uiYiIi0rHFDlBmdgXwE8LkgBZtXm1mv3X3q9uiciIi0nHFHUniSuBy4E/AVGApUEaYKuNKM+vs7tVtVUkR2XJVVVUALFq0iPvuu4+ampoC10gks7hXUGMII0n8NGnbXOBpM6sHqoDqHNdNRHJo0qRJTdYVoKTYxe1mvgPpR5J4IkoXERHJmbgB6kXggDRpB0TpIiIiOZO2ic/MkoPXRcAjZrYReJDN96C+S5hA8IS2rKSIbLmJEycCUFdXx5AhGt9Zil+me1AbCSNGJBhwbbSQsv21Fo4lIgWW6CShoY6kvcgUVCbQNECJiIjkTdqNrBXMAAAVxElEQVQA1Rbdxs2sG3ALcATQG3gLuMTdZ0TphwN/BHYn3Nca7e4LkvLeCpwMrAV+4+43JB07bV4REWl/Yk/5niOdgQ+ArxF6/l0OPGBmA6OR0R+OtvUGZgH3J+WtBgYBA4ARwM/M7BsAMfKKiEg705qRJPYmXL3sBmybkuzuPqqlY7j7Gpo+LzXdzN4FyoE+hKGTHozKqwaWm9lQd58HjAS+5+6rgFVmNgkYTejmfmILeUVEpJ0x95ZvM5nZSGAy4Z7Uv4FPUnZxd9+j1YWblRHml9oPOA/o6u7nJaW/DowHngZWAju7+9Io7WRgvLt/wcxuSpfX3R9KKbOK8GAxZWVl5VOnTm1ttRutrF/Jhk4bss6fjW6bunWIMksppaSkJK9lNjQ0qEyVqTLzUOaIESNmu3tFS/vFvYK6HHgM+L67f5R1rZKYWRfgXuBOd59nZiXAspTd6oFSwvh/ifXUNKL0dHmbcPcaoAagoqLCt6Q305TpU6grqcs6fzaGNAzpEGVWUpn3nmaF6N2WzzLNrMl6nC+nubK1v7cqs23EvQe1M3BLDoPTNsDdhCuxH0SbG4AeKbv2AFZHaaSkJ9JayisiIu1Q3AD1HLB3Lgq08DXudsKDvie5+3+ipLnAvkn7dQf2JNxbWgUsTk6PXs9tKW8u6iwiIvkXN0D9AKgys9PMrI+ZbZO6tKLMWwnB7jh3X5e0/RFgHzM7ycy2Ba4AXkvq5HAXcJmZ9TKzoYQBbO+ImVdERNqZuIHlQ+AV4B5CJ4n/pCypnSaaZWYDgLGEThFLzKwhWs5w92XAScCvgFXAgcCpSdnHE6aaXwD8Hfituz8BECOvSIfn7rg7M2fOzOv9J5Fsxe0kMQk4BXgUmEfMgJQqenDWMqQ/BQxNk7aBMO7f2a3NKyIi7U/cAHUC8FN3v6ktKyMiIpIQt4lvDfBGW1ZEREQkWdwA9Wfg9LasiIiISLK4TXwLgNPM7K+EoYVWpe7g7pNzWTEREenY4gaoW6OfA4DDm0l3wlBIIlKkampqgDBh4fz58xvnhxIpVnED1OfatBYi0ubGjh3bZF0BSopdrACleZVERCTf8j0flIiISCyxrqCiOZsyPnqezXQbIpI/Y8aMAWDRokX069evwLURaVnce1B/57MBqg/wFcJI4k/nslIiknuJThKFmJ5BJBtx70GNbm67mfUkdDt/Kod1EhER2bJ7UNH8UL8ljB4uIiKSM7noJLEe6J+D44iIiDSKew/qM8ysM7APUI0mBhQRkRyL24vvU9L34vsYOCZnNRIRESH+FdQEPhug1hPG6Jvh7vU5rZWI5Fx5eTkAq1evprS0lNmzZxe4RiKZxe3FV93G9RBpYnHDYqprq/NaZiWVeS0v315++eVCV0GkVTSShIiIFKW0V1Bm1qqu4+4+YcurIyIiEmRq4quOkT/5vlSsAGVmPwBGA18ApiQ/BGxmhwN/BHYHXgRGJwaqNbNuhGk/TgbWAr9x9xvi5BWJY2tvVpw1a1bjz4qKiryVK5KtTE18XVpYDgCeBAx4qxVlLgKuJmX+KDPrCzwMXA70BmYB9yftUg0MIsxJNQL4mZl9I2ZekQ6vvLyc8vJyhgwZ0thhQqSYpQ1Q7r6puQXYA7iHcJXyeaAq+hmLuz/s7o8CK1KSTgTmuvuD7r6eEJD2NbOhUfpI4Cp3X+XubwKTCFdicfKKiEg7Y+4ZBynfvKPZbsB4QqBYBVwD3OLun2RVsNnVQP9EE5+Z3QR0dffzkvZ5PSrzaWAlsLO7L43STgbGu/sXMuV194dSyq0iBFXKysrKp06dmk31AVhZv5INnTZknT8b3TZ1U5lbUZmllFJSUpLXMhsaGlSmyixomSNGjJjt7i22M7fYzdzMdgIuJfxTX0+413Sju6/JunbNKwGWpWyrB0qjtMR6alpLeZtw9xqgBqCiosK3ZFTnKdOnUFdSl3X+bAxpGKIyt6IyK6nM+8jihRjNXGWqzGxk6sW3A/Bz4ELCfaabgOvcfVUb1aUB6JGyrQewOkpLrK9PSWspr4iItEOZrqDeBXYgdIS4GlgM9DKzXs3t7O7vbGFd5gKjEitm1h3Yk3BvaZWZLQb2Bf4a7bIvm8cATJt3C+sk0qby2XNwUd0iAMopp7S0VB0lpOhlClA9o59HAUfGOFanOAVGg8x2jvbvZGbbAhuBR4DfmtlJwOOEKTxec/d5Uda7gMvMbBZQBowBvheltZRXpMObdO6k8JPwM+79Z5FCyRSgvpchbUtcRuj4kHAmcKW7V0cB5mY29xI8NWm/8YTnoBYA6wjNjU8AuPuyFvKKiEg7kzZAufudbVFgNK5fdZq0p4Bmu4a7+wbg7GhpVV4REWl/sp4PSkTal10G7QJACSWUln6mg6tI0VGAEukgqmqqgM3d6fM5rNPWPlK8tA2NZi4iIkVJAUpERIqSApSIiBQl3YMSkTa3tU9lIm1DV1AiIlKUdAUl0kFMu34aAP/Y+A/qO9dz3MXHFbhGIpkpQIl0EC8//nKTdQUoKXYKUCKyVdJ9r/ZP96BERKQo6QpKpIM4dtyxAJRtKGNpt6UFro1IyxSgRDqI8uPC/E+FmDlY2s7W3JSpJj4RESlKuoISEcmRQlzNDGFIXsvLJ11BiYhIUVKAEhGRoqQAJSIiRUn3oEQ6iCtHXNlkffzM8QWqiUg8W9UVlJn1NrNHzGyNmS0ws9MLXScREcnO1nYF9UfgE6AM2A943Mxedfe5ha2WiIi01lZzBWVm3YGTgMvdvcHdnwX+BzirsDUTEZFsmLsXug45YWZfAp539+2Stl0MfM3dj0vaVgVURatDgC15pL4vsHwL8qtMlakyVWZHLHOAu+/Y0k5bUxNfCVCfsq0eKE3e4O41QE0uCjSzWe5ekYtjqUyVqTJVpspsaqtp4gMagB4p23oAqwtQFxER2UJbU4CaD3Q2s0FJ2/YF1EFCRKQd2moClLuvAR4GJphZdzM7BDgBuLsNi81JU6HKVJkqU2WqzM/aajpJQHgOCpgMfB1YAfzC3e8rbK1ERCQbW1WAEhGRrcdW08QnIiJbFwUoEREpSgpQKeKO52dmO5nZFDNbZGb1ZvacmR2Yss/p0THWmNmj0T2ygjCze8xssZl9bGbzzeycpLTDzWyema01s5lmNiDmMf9sZm5meyVtK4rxEM1skJmtN7N7ovVKM/vUzBqSllEZ8o82s00p+1cmpQ+M3qu10Xt3RB5OK7l+tdH5JepWF21v1XlGefYws+lmttrMlpvZb5LSCv55mtmpZvZmVIe3zWx49P57ynlenuEYl6Tsuy56n/pG6d3MbHL097HEzMbl8fwaUpZNZvaHKK1V5xnludDM3o3OZZaZHZqUZmZ2nZmtiJbfmJm19Tlmzd21JC3AFOB+woO/hxIe9h3WzH57AOOAXYBOhNEplgMlUfowwjNYX42OdR8wtYDnNQzoFr0eCiwByglPhNcD3wG2BX4L/DPG8Q4F/gE4sFdr3788nO+TwDPAPdF6JfBhK/KPBp7NkP4CcAOwHWGIrY+AHfN4frXAOc1sb+15dgXejn6Xu0e/A18sls+T0OFpAXAQ4Qv1rtEyMPrd65zlcauBp5PWr4l+X3oBe0d/H98owO9td8IznV+N1lt1nsCBwJrob9uA84BlQKcofSxh9Jz+0fv4BnBuvs8z9vtR6AoU0xL9cnwCDE7adjdwbcz8HwPl0etfA/clpe0ZHbu0CM5zCLAY+C4hsD6f8h6sA4ZmyN8ZeAX4YnKA2tL3L4fndyrwQPRPKOcBChgMbEj+LKN/bnn7Qyd3AaoKeCZNWsE/T+B54PvNbM86QEX/uN8GRiVtWwgcmbR+FQX4QgmMAt5hcwe21gaoU4D/TfkMHdgl6f2sSkr/PjG+kBZqURNfU4OBTe4+P2nbq4Srj4zMbD/Ct9G3ok3DorwAuPvbRH/sOattK5nZLWa2FphHCFD/j8/Wcw3hjzfTOf8Y+Ie7v5ayPev3L1fMrAcwAfhJM8k7mdnSqPnjRgsDDGfypajJa76ZXW5miaHBhgHvuHvyKCV5Pc/INVH9nktufqR153kQ8J6ZzYiOVWtmX4jSCvp5mlknoALY0czeMrMPzexmM9suabcF0fY/J5rrYhhOmPHgoaicXkA/kv4OKMznCSFA3eVR9EgS9zxnAJ3M7MDo/TsbmEO4IoSUv3cKd56xKEA1FWs8v1TRP8W7gSvdPZE/q2O1JXc/Pyp/OOGh5g20sp5mthuhmeCKZpKL4ZyvAm539w9Sts8jTMGyC3AYoQnkhgzH+QewD7AToQnvNOCnUVoxnOfPCc3MuxIempxmZnvS+vPsT7ji/D3hn/TjwGNm1pXCn2cZ0AU4mfA7ux/wJeAyQnP6AcAAwjmWAvfGPO4o4L/dvSFaL4l+Jp9r3v9WzWx34GvAnUmbW3ueqwmB91nC3/d4whVTIuClfqb1QEmx3odSgGoq7Xh+ZjY36Sbl8ERi9G1uGuEy+Zo4x2qDesfm7ps8TEXSn9A+nemchyedc2LIqP8CJiQF4mQFPefoKvYI4MbUNHdf4u5vuPun7v4u8DPCPz7M7Iyk85wR7f+Ou78b7f9/hKuyk6PDFfyzdfcX3X21u29w9zuB54CjW3uehObcZ919hrt/AlwP9CHchyn0ea6Lfv7B3Re7+3JCsD3aw5Q6s9x9o7svBX4AHGlmPdL83gKNf6/foWkQSASq5HMtxN/qSMJn8W5iQxbneQ7hqmkYoUXnTGC6mfWL0lM/0x5AQzNXbEVBAaqptOP5ufswdy+Jlmcg9PwBHiW0X49NOdbcKC/RvnsA3aIyikFnwn2x1Hp2T2x392eSzjnRDHA48Nuop1Oi2eAFC727Cj0eYiWhzf79qG4XAyeZ2cvN7OuEexG4+71J5/nNNMdu3J9wPnuYWfI37EKP+5hcv2a3pznP16J9mlPQz9PdVwEfkr5+TXaPflqa39uEE4GVhHt4yeUsJunvgMJ8niNpGjib09J57gtMc/f50ZeUJwjn9pUovcnfO4X/vc2s0DfBim0BphJ6LnUHDiF9L74uhCunR2nmBibhG8zHhKaJ7sA9FKgXH6GZ6lTC5X0n4ChCT58TgB2jczyJ0IPrOjLcNI2OtXPS4oT7GNu15v1ro/PcPqVu1wP/HZ1jJbA74Z/1bsBM4M8ZjvVNoCx6PRR4HRiflP7P6PjbAt8mj734gJ7RZ7gt4YvGGdHnOSSL8xwCrCVceXYi3F98G+ha6M8zKn8C8FL0e9eL0BnlKkJvtSGEL9l9CD0NZ8Y43pOEFoDU7dcCf4/KGEr4p563XnyEALKGlE5UrT1PQvPlfELzrxF6Qa4l6vQEnAu8SWga7kcITurF114WoDch6KwB3gdOT7Pf16J/zmsJl82JZXjSPqdHx1gDPAb0LtA57Rj98X1ECJr/B4xJSj+CcO9iHeGb5cBWHDu1m3ms9y9P513N5l584whXumuBD4A/pP4zSMl7PbA0Oo93on+UXZLSB0bv1TpCt90j8vx5vkRogvqIECy/ns15RnlOJHTu+Tg6p2FJaQX9PAlfBG+JznMJ4V7ZtoR7gu9G9VoM3AXs3MKxdgU2Jv++JqV1I4zj+XH0uY/L83lOBO5uZnurzpMQlCZEn9VqQjA6KyX9N4SryJXRa8vnubZm0Vh8IiJSlHQPSkREipIClIiIFCUFKBERKUoKUCIiUpQUoEREpCgpQImISFFSgBLJwMz+FM3Hk2k8u3zVZXcz+72Z1VmYz2ithXmSbjOz/QtdP5Fc03NQImlE47YtIYxX9m9gV3ffWKC6VBIe9l5KeHD1NcJDl18kTA0y2N23S5dfpD1SgBJJIxpf8F7CtCRHA8e5+/QC1KMPYbSKOsKIEWtT0g04z91vyXAMI4yE8UmbVlYkh9TEJ5LeKGAV4QplHWEwz88ws9MsTPu+3sz+z8yOj+ZVqk3Zr6+Z3WpmC81sQ5SnKkY9xhDGYftBanAC8KBJcDKz98zsHjM728zmEeYiOyZK28XM7ormf9pgZq+Z2Zkp+avN7DPfXs3sDjN7L2k9MSX5+WZ2g5n9O2p6nG5mA2Ocm0hanVveRaTjiaYnOAKocfdlZvYocKKZ9fIw+nViv68TrrL+hzBJYl/ClCTbkjRyfTRn2HOEKeKrCeOrHQXcambd3P0PGapzOLDY3V9p5WmMIMyhdCWhifK9aLT6xKColxDG6zsTuNvMtnf3mlaWkfBLwsR43yMM7Ppr4EkzG+bu/8nymNLBKUCJNO8sQgvDXdH6nYSBO08Bbkva70rgDeDbnhiN0+z/gNk0nVrlh4RJ577g7v+Ktj1lZj2B8WZ2a4b7W/2BBakboxlTk6fY2ORN2+x7AeXuviQpzw+AQcAId6+NNs8wszLgajO73d03palHJquBE9z906ic+YRJ80YCt2dxPBE18YmkMRL4l7u/EK0/BSwiqZnPNk9J/lByYHD3lwlXSMm+AbwIvGtmnRML8BdC893nM9Ql3Wync4H/JC2Hp6T/Mzk4Rb4KLEwKTgn3EEZJz1SPTP47EZwA3P05wlxOB2d5PBFdQYmkMrMDCP+or4uucBIeBn5gZoPdfT6hOa8Lofks1dKU9Z2AvQiBpDl9MlTpA5oPHCcRmgzLaXpVl7C4mW2902xfkpSejdTzTWzbNcvjiShAiTRjVPTz59GSaiRwGbCcEHB2amafMsKcPAkrCIHsh2nKrMtQn6eBI8zsS8n3odx9LoCZlaTJ11wX3ZWECfBS7ZxUT4D10bG7pvT8SxdIy9Jsm5Nmf5EWqYlPJImZdSXMPvwioZNB6jIHOMvMLLpXM4swrbwlHaMc+FzKoZ8gzNT6vrvPamZZnaFakwi9Cf9gZttv4Sn+HehvZoekbD+dEEDfjNYT97z2SewQXU1+headbGbbJO17COHe2Qtp9hdpka6gRJo6lnCV8JNm7tNgZhOBWwlTq88ExhOmEX/EzGoIzX7VhCazT5Oy3kjoYPGMmd1IuGLqTghaw939hHQVcvflZvYd4BFgjpklHtT9lDCt+0jC1dKaGOd3B+Eq7mEzu5Rwn+gMwtTgY5M6SMwgTO8+yczGE2ac/Rlh1ujmlAKPRu/PjsA1wL/Y3MlEpPUKPaWvFi3FtBBGa/gY2D5N+g6E6dTvSNp2OiHgbCB0XPg28ArwSEreXoRA9S7huaR/A88AP4pZt4HAzYR//OsJz2a9SQiY+6Xs+x7RdPfNHGcX4G5CE+UGQrA7s5n9DiVMLb+W0CPxTEKAey+lTg6cD9wALIv2fxz4XKE/Ty3te9FIEiI5Zmb9gbeAX7n7VYWuT1uKHsZ9Fxjj7n8qbG1ka6MmPpEtEI3XdwOhG/pyYA9CU9haQP+wRbaAApTIltlE6AF3M+He1RpCs9133L257twiEpOa+EREpCipm7mIiBQlBSgRESlKClAiIlKUFKBERKQoKUCJiEhR+v/Ar8tYwQz+lQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.grid()\n",
    "plt.hist( demo.demoData[:,0], 11, alpha = 0.5, color = 'g' )\n",
    "bars = ( \"0-20\", \"\", \"30-40\", \"\", \"45-50\", \"\", \"55-60\", \"\", \"65-70\", \"\", \"75-80\", \"\" )\n",
    "plt.xlabel( \"Age Group\", size = 16 )\n",
    "plt.ylabel( \"Number of Individuals\", size = 16 )\n",
    "y_pos = np.arange( len( bars ) )\n",
    "plt.xticks( y_pos, bars )\n",
    "plt.axvline( x = 6, linestyle = \":\", linewidth = 4, color = 'k' )\n",
    "plt.text( 6.5, 4000, \"Seniority Threshold\", fontsize = 12 )\n",
    "ax = plt.gca()\n",
    "# ax.set_xlim( -0.5, 11 )\n",
    "ax.tick_params( axis='both', which='major', labelsize = 12 )\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Distribution of RST values in population**: The following code is meant to reproduce Fig 6(b) in the accompanying paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuYXFWd7vHvayIJphO5DY1XEqANCmM4Q6ujHqAjKt4YEBwNoOGiREBERVHmyKUJeGEcRxGvcbiLBhkDiJcoDDSKjiiIIJFAxBgFEwgQQjrQaQy/88faBZWiq2pXV1X37vT7eZ56umqvvWq/3Unnl733qrUUEZiZmRXNs0Y7gJmZ2VBcoMzMrJBcoMzMrJBcoMzMrJBcoMzMrJBcoMzMrJBcoMzMrJBcoMzMrJBcoMzMrJAmjnaA0bTddtvF9OnTh91//fr1TJkypXWBWsz5muN8zSlyviJng80/3y233PJgRPxD3R0jYtw+9txzz2jG9ddf31T/dnO+5jhfc4qcr8jZIjb/fMDNkePfaF/iMzOzQnKBMjOzQnKBMjOzQnKBMjOzQnKBMjOzQnKBMjOzQnKBMjOzQnKBMjOzQnKBMjOzQhrXUx2Z2fAtX95bs31wcGbdfUZLs9lmzBh+X8vPZ1BmZlZILlBmZlZILlBmZlZILlBmZlZILlBmZlZILlBmZlZILlBmZlZILlBmZlZII16gJE2X9CNJayStkvRlSROztj0k3SLpsezrHmX9JOlsSQ9lj3+XpLL2qn3NzGzsGY0zqK8CDwDPA/YA9gGOk7QFcBXwLWBr4CLgqmw7wDzgQGAW8HLgbcD7AXL0NTOzMWY0CtQM4LsRMRARq4DFwG5AD2nqpS9GxIaI+BIg4HVZv8OBz0fEvRFxH/B54IisrV5fMzMbYxQRI3tA6RjgNcAxpLOdnwCnAtOBN0bEm8v2/QFwfUR8XtLarP2mrK07a5sq6SO1+lYcfx7pbIzOzs49Fy5cOOzvpb+/n46OjmH3bzfna47z1TY4uLJm+8DAJCZP3jBCaRrTbLYttnheC9M802j/2dbTbL7Zs2ffEhHd9fYbjclibwCOBh4FJpAux10JnAKsrdh3LTA1e95R0b4W6MjuQ1W2VfZ9SkQsABYAdHd3R09Pz7C/kb6+Pprp327O1xznq63eZKvLls2kq+uukQnToGazzZhxSAvTPNNo/9nWM1L5cl3ik7SdpBdXbHu/pHMlvS3vwSQ9i3TGtAiYAmxHOos6G+gHplV0mQasy55Xtk8D+iOdAtbra2ZmY0zee1DnAyeXXkg6FfgacChpMMK7cr7PNsCLgC9n94oeAi4A3gIsAV5ePjKPNBhiSfZ8CWmARMmsirZafc3MbIzJW6C6gf8pe30M8OmI2Bb4CnBinjeJiAeB5cCxkiZK2oo0+OE2oA/YCJwgaZKk47Nu12VfLwZOlPQCSc8HPgpcmLXV62tmZmNM3gK1DXA/gKTdgR1I944g3T+a2cAxDwLeBKwG/gj8HfhIRAyShpHPBR4BjgIOzLYDfAO4Gvg9cAfww2wbOfqamdkYk3eQxEPAC7PnrwP+FhHLstfPpoHh6hHxO9Kw8KHabgX2rNIWwMezR0N9zcxs7MlboK4FeiVtR7q0dmVZ267AilYHMzOz8S3vmc/Hgb8CnwHuAc4oazsMuLHFuczMbJzLdQYVEfcDb6jS/HpgoGWJzMzMaMEHdSPi0VYEMTMzK1e1QEk6v4H3iYh4bwvymJmZAbXPoF4H5J2ob2Qn9DMzs81e1QIVEdNHMIeZmdkmvKKumZkVUsODJCRtD0yu3B4Rf2lJIjMzM3IWqGwW8rNIK9huVWW3Ca0KZWZmlvcS34eBD5BWsRXwaVLBWk764O7RbUlnZmbjVt4CdSQwn7RuE8AVEXE68FLgPuDF1TqamZkNR94CtRNwc0RsJM0+viVARDwBfJE0e7iZmVnL5C1Qa3l6YMTf2HR5jYmk5TjMzMxaJu8ovluBl5GWa/8JcIakx0lnU58CftueeGZmNl7lPYP6IvBY9vx0YBVwKXAZaT2o46v024Sk/orHRknnlrXvK2mppMckXS9px7K2SZLOl/SopFWSTqx476p9zcxs7MlVoCLimogorV67Cngl8BJgD+AlEXF7zvfpKD2ATuBx4HKAbK2pRcCppEuGN5MKYEkv0AXsCMwGPi7pTTn7mpnZGDOs2cyz1W3/2OSx3wE8APw8e30QsCQiSgWrF3hQ0q4RsZS0nPuREbEGWCPpm8ARwOIcfc3MbIxRqjV1dpL2rrdPRPysoQNL1wE/i4je7PU5wBYRcWzZPneQLileBzwM7JCtTYWkdwCnR8Q/1uobEd+rOO48YB5AZ2fnngsXLmwk9ib6+/vp6OgYdv92c77mOF9tg4Mra7YPDExi8uQNI5SmMc1m22KL57UwzTON9p9tPc3mmz179i0R0V1vv7xnUH3Un7E890wSkl4M7AOUL9HRAayu2HUtMDVrK72ubKvXdxMRsQBYANDd3R09PT15Yz9DX18fzfRvN+drjvPVtnx5b832Zctm0tV118iEaVCz2WbMOKSFaZ5ptP9s6xmpfHkL1Owhtm0LvI1UaHINkigzF7gxIpaXbesHplXsNw1Yl7WVXg9UtNXra2bWUvWKc7MGB2e2/RjN6RmRo+Rd8v2GKk2LJH0B2B/4cQPHnQt8tmLbEuDw0gtJU4CdSfeW1khaCcwCrsl2mZX1qdm3gUxmZlYgrVhu44fAO/PuLOk1wAvIRu+VuQLYXdLBkiYDpwG3lw1yuBg4RdLWknYlzf93Yc6+ZmY2xrSiQM0Enmxg/8OBRRGxyeW3iFgNHEz64O8a4FXAnLJdTidNTLsCuAH4XEQsztnXzMzGmLzLbcwdYvMWwO6kgQ6L8h4wIt5fo+1aYNcqbRtIc/4NOe9frb5mZjb25B0kcWGV7RtIH4j9UEvSmJmZZfIWqBlDbBsofSbJzMys1fKO4lvR7iBmZmblWjFIwszMrOWqnkFJepL6s0c8JSJyzyRhZmZWT61LfPN5ukCJNHpuS+Bq4H5gB9JMEo8D57Uxo5mZjUNVC1RpElcASaeQPn+0X0Q8VrZ9CmkBw7+3MaOZmY1Dee9BvZ/0wdjHyjdGxHrgP4BjWh3MzMzGt7wFajvSB3OHsgVp4lgzM7OWyVugbgbOkPSC8o3Z617gNy3OZWZm41zeD+qeQFo08B5JvyINkugE/hl4DDi0PfHMzGy8ynUGFRG3ArsAnwc2Av+Yff0PoCsifte2hGZmNi7lPYMiIh4CPtnGLGZmZk/xTBJmZlZItWaSuA44LiKWZs9riYjYt7XRzMxsPKt1BqWK/VTj0dCZmKQ5ku6UtF7SPZL2yrbvK2mppMckXS9px7I+kySdL+lRSasknVjxnlX7mpnZ2FNrJonZZc97WnVASW8AzgbeBfwaeF62fTvSwofvI02ndCZpral/zrr2Al3AjqRplq6X9IeIWJyjr5mZjTG5znwktfKDuGcA8yPiVxHxZETcFxH3AQcBSyLi8ogYIBWkWZJKq+TOBc6MiDURcSfwTeCIrK1eXzMzG2MUUX/CckmDwI+AS4CrI2JwWAeTJpAmlz2NdLYzGbgSOAn4LLBFRBxbtv8dwOmkz2A9DOxQWiRR0juA0yPiHyWdU61vRHyvIsM8YB5AZ2fnngsXLhzOtwJAf38/HR0dw+7fbs7XHOerbXBwZc32gYFJTJ68YYTSNKbI2aD4+QYHpzb1d2/27Nm3RER3vf3yDjM/FTgMuBx4RNJ3gUsi4hcN5uoEng28A9gLeAK4CjgF6ABWV+y/FpiatZVeV7ZRp+8mImIBsACgu7s7enp6GvwWntbX10cz/dvN+ZrjfLUtX95bs33Zspl0dd01MmEaVORsUPx8K1b0jMjfvbwf1D07Il4O/BNwAbA/8LNsgEOvpK6cx3s8+3puRKyMiAeB/wTeAvQD0yr2nwasy9qoaC+1UaevmZmNQQ2NvouI30XER4EXAW8GfgF8FLgzZ/81wL0MvRDiEmBW6UW2lMfOpHtLa4CV5e3Z8yX1+ub6xszMrHCG9UHdiHgSWE86I3qCTYek13MB8EFJ20vaGvgw8APgCmB3SQdLmky6T3V7RCzN+l0MnCJp62zww9HAhVlbvb5mZjbGNPr5pS5J8yXdA/yMdBb1DdLcfHmdSZr9/G7SmdetwKciYjVwMPApYA3wKmBOWb/TgXtICyfeQFqfajFAjr5mZjbG5BokIekDwHuAV5DOnL5HOoO5PvIMAywTEU8Ax2WPyrZrgSGHhkfEBtKy80dVaa/a18zMxp68o/i+CFxLKlJXRMTjdfY3MzNrSt4C9aKIWNXWJGZmZmXyDjN3cTIzsxHl2czNzKyQal3iq5zNvNZgiEaGmZuZmdU14rOZm5mZ5eEVdc3MrJBq3YPau5E3ioifNR/HzMwsqXUPqo+n7zuJ2vegACa0IpCZmRnULlCzy55vBZwL3AEsBO4nLZ1xCLAb8IF2BTQzs/Gp1iCJG0rPJV0I/DQi3lex28WSziOtaHt1WxKamdm4lHeQxAHAZVXaLsvazczMWiZvgXoWsEuVti58/8nMzFosb4H6IfAZSf8qaQKApAmS3gmcRVrPyczMrGXyThZ7AmkV3cuAv0taA2yd9b8xazczM2uZvJPFPhgRewH7kRYc/F729Y0RsXdEPJT3gJL6JA1I6s8ed5W1HSpphaT1kq6UtE1Z2zaSrsjaVkg6tOJ9q/Y1M7Oxp6GZJCLimog4MyKOy75eO8zjHh8RHdljJoCk3Uir876HNIT9MeCrZX2+AgxmbYcBX8v65OlrZmZjTN5LfE+RtD0wuXJ7RPylySyHAVeXZqSQdCpwp6SpwJOkJd13j4h+4EZJ3ycVpJNr9Y2IdU3mMjOzUaA8K7ZLmgacA7wLmDTUPhGRaySfpD7Sh3sF3AV8MiL6JF0F/DIizi7btx/Yh1SgfhkRW5a1fQzYJyL2r9U3Im6pOP48YB5AZ2fnngsXLswTe0j9/f10dHQMu3+7OV9znK+2wcGVNdsHBiYxefKGEUrTmCJng+LnGxyc2tTfvdmzZ98SEd319st7BvUV0hnMecDvgWZ+cp8A/kC6XDcHuFrSHkAHsLZi37XAVGBjjTbq9N1ERCwAFgB0d3dHT0/PcL8P+vr6aKZ/uzlfc5yvtuXLe2u2L1s2k66uu2ruM1qKnA2Kn2/Fip4R+buXt0DtB5wUEV9p9oARcVPZy4skHQK8BegHplXsPg1YRzqDqtZGnb5mZjYG5R0kUboc1w6Rvf8SYNZTB5R2Il1OvDt7TJTUVdZvVtaHOn3NzGwMylugFgL7N3swSVtJ2k/SZEkTJR0G7A38BLgU2F/SXpKmAPOBRRGxLiLWA4uA+ZKmSHotaXqlS7K3rtq32cxmZjY68l7i+ynwxWxE3Y+Ahyt3iIjrcrzPs0kzT+xKuq+0FDgwIu4CkHQMqdhsC1wLHFnW9zjgfOAB4CHg2IhYkh17SZ2+ZmY2xuQtUFdlX2cAR5RtL12eC3LMxxcRq4FX1Gj/NvDtKm0PAwcOp6+ZmY09eQvU7Pq7mJmZtU6uAlW+NpSZmdlIaGiqIzMzs5FS9QxK0sUNvE9ExOEtyGNmZgbUvsS3N2nwQx559zMzM8ulaoGKiOkjmMPMzGwTvgdlZmaF5AJlZmaF5AJlZmaF5AJlZmaF5AJlZmaFVLVASVokaZfs+VxJ245cLDMzG+9qnUEdAGyTPb8A2Ln9cczMzJJaBep+4NXZ89KM5WZmZiOiVoH6LvAFSRtJxelXkjZWefx9ZOKamdl4UWuqo48AvwBeBpwOXAjc16oDZ8u3/x7474h4d7btUOAzwHbANcBR2TpQSNoGOA94I/Ag8G/ZGlDU62tmZmNPramOArgcQNIRwDkRcVsLj/0V4DelF5J2A74BvBX4LbAA+Cowp2z/QaAT2AP4oaTbstV06/U1M7MxJu96UDNaeVBJc4BHgF8Cu2SbDwOujoifZfucCtyZLTP/JHAwsHtE9AM3Svo+8B7g5Fp9I2JdK7ObmdnIUDpRyrGj9Dzgo8A+pNF9DwF9wH9GxKrcB5SmATcD+wLvBXaJiHdLugr4ZUScXbZvf3a8J7O2LcvaPgbsExH71+obEbdUHH8eMA+gs7Nzz4ULF+aN/gz9/f10dHQMu3+7OV9znK+2wcGVNdsHBiYxefKGEUrTmCJng+LnGxyc2tTfvdmzZ98SEd319st1BiXpJcCNwFak+1J/BHYAPgTMlbRXRCzLme1M4LyI+Kuk8u0dwNqKfdcCU4GNNdrq9d1ERCwgXQKku7s7enp6csZ+pr6+Pprp327O1xznq2358t6a7cuWzaSr666RCdOgImeD4udbsaJnRP7u5SpQwNmkf/BfGRF/Lm2UtCPw06z9oHpvImkP4PXA/xmiuR+YVrFtGrCOdAZVra1eXzMzG4PyFqjZwDHlxQkgIlZI6iUNSMijB5gO/CU7e+oAJkh6GbAYmFXaUdJOwCTgblKBmiipq+xMbRawJHu+pEZfMzMbg/IWqC2ofjayLmvPYwFQftPnY6SCdSywPfC/kvYijcSbDywqDXKQtAiYL+l9pFF8BwCvyd7n0lp9zcxs7Mk7WezvgA9K2mR/pdOg47L2uiLisYhYVXqQLs0NRMTqiFgCHEMqNg+Q7h8dV9b9OGDLrO07wLFZH3L0NTOzMSbvGdR84AekoduXAStJgyT+Fegiff6oYRHRW/H628C3q+z7MHBgjfeq2tfMzMaevJ+DWizpbcBZwCd5em6+W4C3RcRP2xfRzMzGo7xnUETEYmCxpOcAWwNrIuKxtiUzM7NxLXeBKsmKkguTmZm1lVfUNTOzQnKBMjOzQnKBMjOzQmr4HpSZJfXmomu3wcGZo57BrJ3qnkFJ2kLSbyW9cSQCmZmZQY4CFRGDwAzAy7qbmdmIyXsP6hrSUutmZmYjIu89qHOBb0maCFxJmupok5UOI+JPLc5mZmbjWN4CdUP29UTgI1X2mdB8HDMzsyRvgTqyrSnMzMwq5J0s9qJ2BzEzMyvX0Ad1JT1L0u6S9pE0pV2hzMzMchcoSR8AVgG3A9cBM7PtV0o6oYH3+ZaklZIelXR3tkJuqW1fSUslPSbpekk7lrVNknR+1m+VpBMr3rdqXzMzG3tyFShJRwPnkEbwvZO0HlTJz4GDGzjmZ4DpETEN+BfgLEl7StoOWAScCmwD3AxcVtavl7Q44o7AbODjkt6U5avX18zMxpi8Z1AnAp+PiHnAFRVtS8nOpvKIiCURsaH0MnvsDBwELImIyyNigFSQZknaNdt3LnBmRKyJiDuBbwJHZG31+pqZ2RiTdxTfDOAnVdrWA1s1clBJXyUVly2BW4EfAZ8CbivtExHrJd0D7CbpfuD55e3Z89IS8LtV60sqoOXHngfMA+js7KSvr6+R6Jvo7+9vqn+7OV9z6uUbHMz9/7K2GBiYxLJlo5uhliLnK3I2KH6+wcGR+d3NW6AeBKZXaZsJ3NfIQSPiOEkfBF4N9AAbgA5gdcWua4GpWVvpdWUbdfpWHnsBsACgu7s7enp6Gom+ib6+Pprp327O15x6+UZ7otZly2bS1XXXqGaopcj5ipwNip9vxYqeEfndzXuJ72rgNEk7lW2L7N7PR0j3phoSERsj4kbghcCxQD8wrWK3acC6rI2K9lIbdfqamdkYlLdAnUI6y7kDuJZ03+hLwJ3ARmB+Exkmku5BLQFmlTZmw9h3Jt1bWkOaXmlWWb9ZWR9q9W0il5mZjaJcBSoiHgK6SSPwng3cQyosXwZeHRFra3R/iqTtJc2R1CFpgqT9gENIw9avAHaXdLCkycBpwO0RUbqHdDFwiqSts8EPRwMXZm31+pqZ2RiTe8HCiFgHnJk9hitIl/O+TiqOK4APR8RVAJIOJhW9bwE3AXPK+p4OfC3r8zhwdkQszrKtrtPXzMzGmIZW1JU0DdgdeAFwL+ny26N5+0fEamCfGu3XAkMODc+Gph+VPRrqa2ZmY0/uAiXpNOCjpBFzpQ/qrpP0uYg4qx3hzMxs/MpVoCSdQZql4b+AhcD9QCfp/tEZkiZGRG+7QpqZ2fiT9wzqaNJMEieVbVsCXCdpLemDr70tzmZmZuNY3mHmz6X6TBKLs3YzM7OWyVugbgJeUaXtFVm7mZlZy1S9xCepvHidAFwh6e/A5Tx9D+qdpFF1B7QzpJmZjT+17kH9nfS5pRIBn80eVGy/vc57mZmZNaRWUZnPpgXKzMxsxFQtUB42bmZmoyn3ku9mZmYjqZGZJF4KvAN4ETC5ojki4vBWBjMzs/Et70wSc4HzSfekHgAGK3bxvSozM2upvGdQpwJXAe+NiEfamMfMzAzIX6B2AI5xcTIzs5GSd5DEL4CXtjOImZlZubwF6nhgnqRDJG0r6VmVjzxvImmSpPMkrZC0TtKtkt5c1r6vpKWSHpN0vaQdK/qeL+lRSasknVjx3lX7mpnZ2JO3QN0L3EparfYB4ImKR+WgiWomAn8lLVr4XNK9re9Kmi5pO2BRtm0b4GbgsrK+vUAXsCMwG/i4pDcB5OhrZmZjTN57UN8E3gVcCSwlf0HaRESsZ9NlOX4gaTmwJ7AtaYXeywEk9QIPSto1IpYCc4EjI2INsEbSN4EjSLOpH1Snr5mZjTF5C9QBwEkRcU4rDy6pE3gJaW2pY4HbSm0RsV7SPcBuku4Hnl/enj0/MHu+W7W+pIJqZmZjTN4CtR74QysPLOnZwKXARRGxVFIHsLpit7XAVNIy86XXlW1k7dX6Vh53HmmBRTo7O+nr6xv299Df399U/3ZzvubUyzc4OHPkwgxhYGASy5aNboZaipyvyNmg+PkGB0fmdzdvgboAOBS4phUHzQZVXEK6VHh8trkfmFax6zRgXdZWej1Q0Vav7yYiYgGwAKC7uzt6enqG+23Q19dHM/3bzfmaUy/f8uW9I5ZlKMuWzaSr665RzVBLkfMVORsUP9+KFT0j8rubt0CtAA6RdA3pns+ayh0i4vw8byRJwHmk9aTeEhFPZE1LgMPL9psC7Ey6t7RG0kpgFk8XyVlZn5p9c35/ZmZWMHkL1NeyrzsC+w7RHqSpkPK+10uB10fE42XbrwA+J+lg4IfAacDtZYMcLgZOkXQzqbgdDRyZs6+ZmY0xeQvUjFYcLPts0vuBDcCqdDIFwPsj4tKswHyZNJz9JmBOWffTScVtBfA4cHZELAaIiNV1+pqZ2RiTq0BFxIpWHCx7H9VovxbYtUrbBtLy8kc12tfMzMYerwdlZmaFlHe5jeXUWVIjInZqSSIzMzPy34O6gWcWqG2B15CGeF/XylBmZmZ570EdMdR2SVuRhp1f28JMZmZmzd2DytaH+hxpWLeZmVnLtGKQxADwwha8j5mZ2VPy3oN6BkkTgd1Js5N7xgYzM2upvKP4nqT6KL5Hgbe2LJGZmRn5z6Dm88wCNUCa1eHHEbH2mV3MzMyGL+8ovt425zAzM9uEZ5IwM7NCqnoGJamhoeMRMb/5OGZmZkmtS3y9OfqX35dygTIzs5apdYnv2XUerwB+Spqd/I/tjWlmZuNN1QIVERuHegA78fSaSy8D5mVfzczMWib3IAlJL5L0X6QP5b4O+BjQFRH/lRWuvO9zvKSbJW2QdGFF276Slkp6TNL12QKHpbZJks6X9KikVZJOzNvXzMzGnroFStL2ks4B7gYOJt1r2ikivhgRg8M45t+As6hYIl7SdsAi4FRgG+Bm4LKyXXqBLtKy87OBj0t6U86+ZmY2xtQaxfdc4BPAB0n3mc4hLbO+ppkDRsSi7P272XQOv4OAJRFxedbeCzwoadeIWArMBY7Mjr9G0jeBI0izqdfra2ZmY0ytUXzLgeeSBkKcBawEtpa09VA7R8SfmsyyG3Bb2futl3QPsJuk+4Hnl7dnzw+s1xdwgTIzG4NqFaitsq/7AW/M8V4TmszSAayu2LYWmJq1lV5XttXruwlJ80gDO+js7KSvr2/Ygfv7+5vq327O15x6+QYHZ45cmCEMDExi2bLRzVBLkfMVORsUP9/g4Mj87tYqUEe2/eib6gemVWybBqzL2kqvByra6vXdREQsABYAdHd3R09Pz7AD9/X10Uz/dnO+5tTLt3x574hlGcqyZTPp6rprVDPUUuR8Rc4Gxc+3YkXPiPzuVi1QEXFR24++qSXA4aUXkqYAO5PuLa2RtBKYBVyT7TKLp5f5qNp3BHKPW83+Az04OHPU/5Gvpej5zDZ3Iz4Xn6SJkiaTLglOkDQ5W1vqCmB3SQdn7acBt5cNcrgYOEXS1pJ2BY4GLsza6vU1M7MxZjQmiz0FeBw4GXh39vyUiFhNGsb+KWAN8CpgTlm/04F7SEt83AB8LiIWA+Toa2ZmY8ywV9Qdrmzpjt4qbdcCu1Zp2wAclT0a6mtmZmOPl9swM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCcoEyM7NCGvHlNjYng4MrC73iqleENbOxzGdQZmZWSJtVgZK0jaQrJK2XtELSoaOdyczMhmdzu8T3FWAQ6AT2AH4o6baIWDK6sczMrFGbzRmUpCnAwcCpEdEfETcC3wfeM7rJzMxsODanM6iXABsj4u6ybbcB+5TvJGkeMC972S/priaOuR3wYBP92835muN8zSlyviJng8LnO6PZfDvm2WlzKlAdwNqKbWuBqeUbImIBsKAVB5R0c0R0t+K92sH5muN8zSlyviJnA+cr2Wwu8QH9wLSKbdOAdaOQxczMmrQ5Fai7gYmSusq2zQI8QMLMbAzabApURKwHFgHzJU2R9FrgAOCSNh62JZcK28j5muN8zSlyviJnA+cDQBExEscZEZK2Ac4H3gA8BJwcEd8e3VRmZjYcm1WBMjOzzcdmc4nPzMw2Ly5QZmZWSC5QdeSd30/SbEnXS1or6c8FzHeSpDskrZO0XNJJBcr2YUl/kvSopL9J+oKktn9Gr9G5GyVtIWmppHvbna2RfJJ6JT0hqb/ssVNR8mX7/pOkn2XZ7pf0oaLkk/Tjip/doKTfFyjfJElfz35uD0v4Iw6eAAAIhklEQVS6WtILCpJtK0kXSXoge/S2MocLVH3l8/sdBnxN0m5D7LeeNECj7f/wV8ibT8BcYGvgTcDxkuYUJNvVwD9FxDRgd9LHA05oc7ZG8pWcBDwwArlKGsl3WUR0lD3+VJR8krYDFgPfALYFdgF+WpR8EfHm8p8d8Evg8qLkAz4EvBp4OfB84BHg3IJk+wLwHGA68ErgPZKObFmKiPCjygOYkv0hvaRs2yXAZ2v0eT3w56LmK9vvS8C5RctG+gfsWuCrRfrZATOAO4E3A/cW6c8W6AW+NRJ/54aZ79PAJUXNV9FvOrARmFGUfMDXgH8ve/1W4K6CZHsQeEXZ6/8H/LxVWXwGVVu1+f1q/S97JA0rnyQBe9HeDzE3lE3SoZIeJf2Fn0X633Y7NfqzO5f0y/d4m3OVNJpv/+zyzxJJx7Y/XkP5/hl4WNIvs8tAV0t6cYHylZtL+gd2eduSJY3kOw94raTnS3oO6YzmxwXJBunqTPnz3VsVxAWqtlzz+42i4ebrJf3ZX9CGTCUNZYuIb0e6xPcS4OvA/W3MBg3kk/R2YGJEXNHmTOUa+fl9F3gp8A/A0cBpkg5pb7yG8r0QOJx0qerFwHLgO21NN/zfjbnAhe0IVKGRfHcDfwHuAx4l/VnPL0i2xcDJkqZK2gU4inTJryVcoGor+vx+DeeTdDzpl/CtEbGhSNkAImIZ6czuq23KVZIrn9IyLv8OfLDNeSrl/vlFxB8i4m8RsTEifgmcA7yjKPlIZ51XRMRvImIAOAN4jaTnFiQfAJL+L7AD8N9tzFXSSL6vAZNJl7+nkGbMaecZVCPZTiD9+S4DriL9x6Nlg4hcoGor+vx+DeWTdBRwMrBvRLR7JFozP7uJwM5tSfW0vPm6SPclfi5pFekfh+dJWiVpegHyDSXY9LJLOzSS7/YsU0npeTszDufndziwKCL625irpJF8s4ALI+Lh7D+V5wKvzAafjGq2LNNhEbFDROxGqim/blmSdt4I3BwewELS/wqmAK8lneruNsR+zyL9L+fNwIrs+RYFyncYsAp4aQF/du8Dts+evyz7RfjPIuQjFcsdyh4HAX/Lnk8Y7XzZfgeQRmeKNJLqPuDwIvz8sv1eB6whrXL9bNLIr5bdSG82X7bvlqTRca9rd65h/PwuAL4HPDf7+f0/4L6CZNuZdGY3Ifu378FqP+Nh5RipP4yx+gC2Aa4kDSP/C3Botn0voL9svx7S/wzLH30FyrcceIJ0+l56fL0g2S4g3XNaD/wZ+BwwuSg/u4o+PYzAKL4Gf37fIc092Q8sBU4oUr5s27GkwrmG9LGCFxUs3yGk/1hqJH52Df75bgtcSvqIwyPAjcArC5LtnaT/sD0G/A7Yr5U5PBefmZkVku9BmZlZIblAmZlZIblAmZlZIblAmZlZIblAmZlZIblAmZlZIblAmVUh6QhJUfYYlHSPpE9LmjzE/vtIWpytaTUg6d7s9WFZe1/F+1V7TB/ivT+Stb2yRt6fK631lXuGBknvy973hXn7mI2Uti8KZ7YZ+FfS/GJTgbcD/5Y9f2p+PkkHkqZB+j5wPPAwsCPwBuAtpA9aHsemc5ydCrwC+JeK460cIsOlpDkB38MQU8lImkH6xP+Z4Q832mbCBcqsvt9FxB+z59dkc5S9V9KHIuLJbPuJwK3A2ysKxEWSngVpUtfyN5W0GhiMiF/VCxARD0haDMyRdGJEPFGxy1zSVEcXN/zdmRWUL/GZNe63pLnbyifr3AZ4YKizl7Ii1qyLsmO+eYi2dwO/iIh7ACRtKemcbH2o9ZJWSvq+pJm1DiBpYnbJ75SK7btk299dsX22pOv09HLpP5b0sua+TbPEBcqscdNJk2c+VLbt18AbJZ0l6eWN3AdqwNWkuezeU75R0mtIy6hfVLZ5y+wxn3SJ8QOkiT//V9L2rQgj6QDgGtL8cIeSiuTWpJnfX9CKY9j45kt8ZvVNkDSRp+9BHQx8OCI2lu1zMqlIfDJ7PCrpeuDbEfHdVoSIiA2SFgJHSdoqIh7JmuYCA6SFC0v7PgzMK72WNAH4CbAaeBdpyYZhywrwOcD/RMRBZdv7gD8BHwE+1swxzHwGZVbfUtJM8A+Tlt/+RkR8uXyHiHggIvYmLXdxGvBz4PXAZZK+2cIsFwGTSLNII6n0/KqI2GQVVElzJP1a0lrg76TZzrcEal7my2lX0iCQS7PLghOzIt4P3ATs3YJj2DjnAmVW39tJo+3eAlwLHCdp7lA7Rlo19syIeBtpqfP/Ad4nafdWBImIm0gFs3T8/UmX1cov75WWqf8OcAdpKYlXZd/Dw6S1yppVukx4Eal4lz/eRFoiwqwpvsRnVt8dpVF8kq4jrRD7OUnfi4j11TpFxCOSvgTsS1qI8Y4W5bkY+LSknUj3o1YBP63YZw6wNCKOKm3IPru1VZ333kgqMltUbK8sOKX7bx8Hrh/ifTbUOY5ZXT6DMmtApCW3TyKdQRxX2i7pRVW67Jp9HeqzTcN1CfAk8GHSiL5LK+6HATyHdFmv3Fzq/M5noxD/ClSe8b214vUfsv1eFhE3D/H4ff5vx2xoPoMya1BEfF/Sb4CPSfpyRDwO/EjSA8BlwF2kez17Ax8F/hf4RQuPf292Jnc86bNPFw2x22Lgy5L+A/gx6fLeB4BHcxxiIfAJSf8G/Ib0fcypyPCkpOOBRdmZ2eWks6odgNcAf4qIc4bz/ZmV+AzKbHhOIZ1FHZO9Pov0D/QnSAXhCuBA4POkZbBb9VmokotIxel3Vc5Wvg58hjT8+2pgP+BtwLoc731W1v8E0uwYXTx9z+spEfF9YB/S7BjnkUYJfpb0c7mpsW/H7Jm85LuZmRWSz6DMzKyQXKDMzKyQXKDMzKyQXKDMzKyQXKDMzKyQXKDMzKyQXKDMzKyQXKDMzKyQ/j+kSxNhNmYJoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.grid()\n",
    "plt.hist( 1 - newSUS, 5, alpha = 0.5, color = 'y' )\n",
    "plt.xlabel( \"RST Value\", size = 16 )\n",
    "plt.ylabel( \"Number of Individuals\", size = 16 )\n",
    "ax = plt.gca()\n",
    "ax.tick_params( axis='both', which='major', labelsize = 12 )\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIDX_S = 0    # Number of susceptible but non-recovered individuals\n",
    "SIDX_E = 1    # Number of exposed individuals\n",
    "SIDX_I = 2    # Number of infectious individuals\n",
    "SIDX_Q = 3    # Number of quanrantined individuals\n",
    "SIDX_R = 4    # Number of recovered individuals\n",
    "SIDX_X = 5    # Number of expired individuals\n",
    "SIDX_V = 6    # Average virulence of the viral strains in E and I populations\n",
    "SIDX_EI = 7    # Number of infected individuals\n",
    "SIDX_D = 8    # Number of individuals infected each day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From the normalized parameters being optimized by ESOP, recover the I0, P0, L0 values\n",
    "def getIPL0( x ):\n",
    "    global IPL0c, IPL0w, offset, ul, ll\n",
    "    IPLVal = np.floor( IPL0c + (np.array(x) - 0.5) / 0.5 * IPL0w ).astype(int) + offset.astype(int)\n",
    "    # Make sure that the recovered values are legal\n",
    "    correctedIPLVal = np.minimum( np.maximum( IPLVal, ll ), ul )\n",
    "    return correctedIPLVal\n",
    "\n",
    "# Get a lockdown schedule corresponding to a certain IPL value\n",
    "# The IPL parameter stores lockdown levels from 0 to 10 instead of 0 to 5\n",
    "# Thus, before using them, divide by 2\n",
    "def getLKP( IPL ):\n",
    "    I0 = IPL[0]\n",
    "    P0 = IPL[1]\n",
    "    L0 = IPL[2] / 2\n",
    "    LKP = np.zeros((T,))\n",
    "    LKP[I0:I0+P0] = L0\n",
    "    return LKP\n",
    "\n",
    "# Find out the objective value corresponding to the stats sent as input\n",
    "def obj( stats, IPL ):\n",
    "    global SIDX_EI\n",
    "    fEpi = np.max( stats[ SIDX_EI, : ] )\n",
    "    # One day of lockdown at level l causes (l/50)% of population to lose their jobs\n",
    "    daysLockdown = IPL[1]\n",
    "    # Remember, the IPL stores levels from 0 to 10 instead of 0 to 5\n",
    "    level = IPL[2] / 2\n",
    "    fEco = level / 5 * daysLockdown * N / 1000\n",
    "    return [fEpi, fEco, fEpi + fEco]\n",
    "\n",
    "# Ask the VIPER simulator what does it think will happen if lock-down is\n",
    "# initiated as specified in the normalized parameter x\n",
    "def askSimulator( x ):\n",
    "    global IPL0c, IPL0w, evalCount, globalDict, freshMask\n",
    "    # Before starting a simulation, turn back time to reset everything\n",
    "    pop.reset()\n",
    "    # Also seed the RNG so that simulations are replicatable\n",
    "    np.random.seed(0)\n",
    "    IPL = getIPL0( x )\n",
    "    # Avoid simulating again for previously queried points\n",
    "    if tuple(IPL) in globalDict:\n",
    "        freshMask.append( False )\n",
    "        return globalDict[tuple(IPL)][-1]\n",
    "    LKP = getLKP( IPL )\n",
    "    stats = pop.simulate( T = T, LKP = LKP, minimal = True )\n",
    "    globalDict[tuple(IPL)] = obj( stats, IPL )\n",
    "    freshMask.append( True )\n",
    "    evalCount += 1\n",
    "    return obj( stats, IPL )[-1]\n",
    "\n",
    "# Once ESOP is done, find out which parameter, i.e. which values of I0, P0, L0 won!\n",
    "def getWinner():\n",
    "    global globalDict\n",
    "    keys = list( globalDict.keys() )\n",
    "    vals = np.array( list( globalDict.values() ) )[:,-1]\n",
    "    winnerIdx = np.argmin( vals )\n",
    "    winner = np.array( keys[ winnerIdx ] )\n",
    "    return winner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Multi-scale Optimization**: ESOP performs Bayesian optimization in a multi-scale manner in several epochs. Initial epochs perform coarse optimization to roughly identify the region with promising objective values. Later epochs zoom into those regions to more and more finely search to obtain points offering highly optimal objective values.\n",
    "\n",
    "In multi-dimensional optimization cases such as this, where 3 parameters are simultaneously getting optimized, this procedure is made slightly more robust by retaining a diverse _active set_ of good-performing candidates from each epoch and searching around them in the next epoch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getAllPairsDistances( A, B ):\n",
    "    weights = np.array( [ 0.01, 0.01, 0.1 ] )\n",
    "    Aw = weights * A\n",
    "    Bw = weights * B\n",
    "    squaredNormsA = np.square( np.linalg.norm( Aw, axis = 1 ) )\n",
    "    squaredNormsB = np.square( np.linalg.norm( Bw, axis = 1 ) )\n",
    "    dists = squaredNormsA[:, np.newaxis] + squaredNormsB - 2 * Aw.dot( Bw.T )\n",
    "    dists[ dists < 0 ] = 0\n",
    "    return dists\n",
    "\n",
    "def addCandidates( res ):\n",
    "    global branchBound, candidates\n",
    "    perm = np.argsort( res.func_vals )\n",
    "    for i in range( 2 * branchBound ):\n",
    "        key = tuple( getIPL0( res.x_iters[ perm[i] ] ) )\n",
    "        val = res.func_vals[ perm[i] ]\n",
    "        if key not in candidates:\n",
    "            candidates[key] = val\n",
    "\n",
    "def updateActive():\n",
    "    global branchBound, candidates\n",
    "    keys = list( candidates.keys() )\n",
    "    vals = list( candidates.values() )\n",
    "    perm = np.argsort( vals )\n",
    "    # Screen the top few best performing keys\n",
    "    n = min( 2 * branchBound, len( candidates ) )\n",
    "    keysScreened = np.array( [ keys[ perm[j] ] for j in range(n) ] )\n",
    "    activeSet = np.zeros( (branchBound, keysScreened.shape[1]) )\n",
    "    # Definitely include the best performing candidate in the active set\n",
    "    i = 0\n",
    "    activeSet[0,:] = keysScreened[i,:]\n",
    "    # To select the rest of the active set, perform the k-means++ algorithm\n",
    "    for t in range( 1 , min( branchBound, len( candidates ) ) ):\n",
    "        dist = np.min( getAllPairsDistances( keysScreened, activeSet[0:t,:] ), axis = 1 )\n",
    "        probs = dist/np.sum(dist)\n",
    "        i = np.argmax( probs ) # Alternatively, we could have used np.random.choice( np.arange(n), p = probs )\n",
    "        activeSet[t,:] = keysScreened[i,:]\n",
    "    return activeSet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch =  0\n",
      "subepoch 1 of 2\n",
      "subepoch 2 of 2\n",
      "epoch =  1\n",
      "subepoch 1 of 2\n",
      "subepoch 2 of 2\n",
      "epoch =  2\n",
      "subepoch 1 of 2\n",
      "subepoch 2 of 2\n",
      "epoch =  3\n",
      "subepoch 1 of 2\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "subepoch 2 of 2\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch =  4\n",
      "subepoch 1 of 2\n",
      "subepoch 2 of 2\n"
     ]
    }
   ],
   "source": [
    "# All parameters are optimized in various ranges that are described using\n",
    "# the center point of the range, the half-width of the range and the initial\n",
    "# offset of the range e.g. the range with center c, width w and offset o is\n",
    "# the range [o + c - w, o + c + w]. Each of the three parameters I0, P0, L0\n",
    "# have their own centers, widths and offsets specified below\n",
    "IPL0c = np.array( [25, 20, 5] )\n",
    "IPL0w = np.array( [25, 20, 5] )\n",
    "offset = np.array( [12, 0, 0] )\n",
    "\n",
    "# Upper and lower bounds on legal values of lock-down initiation points\n",
    "ul = np.array( [62, 40, 10] )\n",
    "ll = np.array( [12, 0, 0] )\n",
    "\n",
    "branchBound = 2\n",
    "activeSet = []\n",
    "evalCount = 0\n",
    "globalDict = {}\n",
    "freshMask = []\n",
    "\n",
    "candidates = {}\n",
    "res = gpm( askSimulator, [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)],  acq_func = \"EI\", n_calls = 20, n_random_starts = 5, random_state = 0 )\n",
    "addCandidates( res )\n",
    "activeSet = updateActive()\n",
    "# Halve the width of the region in which we are searching\n",
    "IPL0w = IPL0w / 2\n",
    "\n",
    "for epoch in range(5):\n",
    "    print( \"epoch = \", epoch )\n",
    "    candidates = {}\n",
    "    for subepoch in range( len( activeSet ) ):\n",
    "        print( \"subepoch %d of %d\" % ( subepoch + 1, len( activeSet ) ) )\n",
    "        IPL0c = activeSet[subepoch]\n",
    "        res = gpm( askSimulator, [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)],  acq_func = \"EI\", n_calls = 20 - epoch, n_random_starts = 5, random_state = 0 )\n",
    "        addCandidates( res )\n",
    "    # Refine the search in the next epoch\n",
    "    # Choose new candidates to populate the active set\n",
    "    activeSet = updateActive()\n",
    "    # Halve the width of the region in which we are searching\n",
    "    IPL0w = IPL0w / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Convergence Plot**: The following plot shows how fast ESOP converges to a near-optimal point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Start on day 27, last for 40 days, at level 3.000000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X2cVnWd//HXmxu5CRLFsRRwQcXwJkOcUMu2olJ0U9yixO2nbGvyWNc2bbM2c8skzbuy1v2Vu5asuj83MjQl1zJKarMSBEUUER1Tg6BEAcWAQfDz++N8B68Zrpk5c3Ou65qZ9/PxuB6c8z3fc85nzoPhw/fm+h5FBGZmZrWmX7UDMDMzK8cJyszMapITlJmZ1SQnKDMzq0lOUGZmVpOcoMzMrCY5QZmZWU1ygjIzs5rkBGVmZjVpQLUDKMI+++wTY8eOrXYYZmZWxtKlS1+IiLr26vXKBDV27FiWLFlS7TDMzKwMSc/lqecuPjMzq0lOUGZmVpOcoMzMrCY5QZmZWU1ygjIzs5rkBGVmZjWpV04z76w//elPXHnllR06RxL19fWcccYZSCooMjOzvscJqsTGjRv55je/2alz165dy4UXXtjNEZmZ9V3u4usmd911V7VDMDPrVZygusmmTZuqHYKZWa/iLr4S++67L9dee22uuhs2bOCyyy7btf/KK68UFZaZWZ+kiKh2DN2uvr4+il6L7/nnn+dNb3rTrv199tmH9evXF3pPM7PeQNLSiKhvr567+Dpp+PDhzfY3b95cpUjMzHonJ6hOGjx4MP36vf74GhsbefXVV6sYkZlZ7+IE1UmSGDZsWLOyP//5z1WKxsys93GC6oKWCcoTJczMuo8TVBd4HMrMrDhOUF3gFpSZWXGcoLrACcrMrDhOUF3gBGVmVhwnqC7wGJSZWXEKS1CSBktaLOkRSSskXZrKJelySU9KWinpUyXl10lqkLRc0qSSa82U9FT6zCwq5o5yC8rMrDhFrsXXCEyJiFckDQTul/Rj4FBgDDAhIl6TtG+qfxIwPn2OAa4HjpG0N3AJUA8EsFTS/IjYWGDsuThBmZkVp7AWVGSa/sUemD4BnAvMjojXUr3nU51pwC3pvAeAEZL2A04EFkTEhpSUFgBTi4q7I5ygzMyKU+gYlKT+kpYBz5MlmUXAQcDpkpZI+rGk8an6KGB1yelrUllr5S3vNStdc0mlFm31GJSZWXEKTVARsTMiJgKjgcmSjgAGAdvSSrbfAeak6uXelx5tlLe81w0RUR8R9XV1dd3zA7TDLSgzs+JUZBZfRGwCfkHWNbcGuD0d+iFwZNpeQzY21WQ0sLaN8qpzgjIzK06Rs/jqJI1I20OA9wNPAHcCU1K1dwNPpu35wFlpNt+xwEsRsQ64FzhB0l6S9gJOSGVV5wRlZlacImfx7QfcLKk/WSK8LSLulnQ/cKukTwOvAJ9I9e8BTgYagC3AxwEiYoOkrwAPpnqzI2JDgXHn5jEoM7PiFJagImI5cFSZ8k3AX5UpD+C8Vq41h9fHqmqGW1BmZsXxShJd4ARlZlYcJ6guaJmg3MVnZtZ9nKC6wC0oM7PiOEF1QctJEk5QZmbdxwmqCwYNGkT//v137W/fvp3t27dXMSIzs97DCaoLJLmbz8ysIE5QXeQEZWZWDCeoLvI4lJlZMZygusgtKDOzYjhBdZG/C2VmVgwnqC5yC8rMrBhOUF3kMSgzs2I4QXWRW1BmZsVwguoij0GZmRXDCaqL3IIyMyuGE1QXeQzKzKwYTlBd5BaUmVkxCktQkgZLWizpEUkrJF3a4vi/SXqlZH+QpO9LapC0SNLYkmMXpfJVkk4sKubO8BiUmVkximxBNQJTIuJtwERgqqRjASTVAyNa1D8b2BgRBwPfAK5KdQ8DZgCHA1OBb0vqT41wC8rMrBiFJajINP1rPTB9IiWXa4DPtThlGnBz2p4HvE+SUvnciGiMiGeABmByUXF3lMegzMyKUegYlKT+kpYBzwMLImIR8ElgfkSsa1F9FLAaICJ2AC8BI0vLkzWprCa4BWVmVowBRV48InYCEyWNAH4o6S+BjwDvKVNd5S7RRnnzk6VZwCyAAw44oLMhd5jHoMzMilGRWXwRsQn4BfBe4GCgQdKzwFBJDanaGmAMgKQBwJ7AhtLyZDSwtsw9boiI+oior6urK+gn2V3LBLVy5UqOP/54FixYULEYzMx6oyJn8dWllhOShgDvB5ZGxJsjYmxEjAW2pEkRAPOBmWl7OnBfREQqn5Fm+Y0DxgOLi4q7o1qOQQH8+te/5qMf/SgbN26sQkRmZr1DkS2o/YCFkpYDD5KNQd3dRv0bgZGpRfVPwOcBImIFcBvwOPAT4LzUdVgT9txzT0aMaDkhETZt2sT8+fOrEJGZWe9Q5Cy+5RFxVEQcGRFHRMTsMnWGlWxvi4iPRMTBETE5In5XcuzyiDgoIt4SET8uKubOGDBgAFdddRWDBg3a7di8efOqEJGZWe/glSS6waxZs9i4cSN33XVXs/Kf/vSnvPTSS1WKysysZ3OC6iZDhgzhlFNO4eCDD95Vtn37du6+u61eTTMza02h08z7GklMnz6dK6+8clfZTTfdxIEHHrhb3UMOOYSRI0dWMjwzsx5F2US53qW+vj6WLFlSlXsvXbqU+vr6duv169eP//7v/+b000+vQFRmZrVD0tKIaPcfSnfxdbNJkyYxbty4duu99tprXHHFFRWIyMysZ3KC6maSmDVrVq6669a1XO3JzMyaeAyqABdeeCGNjY387Gc/Y8eOHbvKd+7cyYMPPrhrf9u2bdUIz8ysR/AYVAVt27aNIUOG7NrfY489aGxsrGJEZmaV5zGoGtTyy7zbt29n586aWRTDzKymOEFVkCQGDx7crMwtKDOz8pygKqxlgvI4lJlZeU5QFeYEZWaWjxNUhTlBmZnl4wRVYaWz+AC2bt1apUjMzGqbE1SFuQVlZpaPE1SFOUGZmeXjBFVhTlBmZvk4QVWYE5SZWT6FJShJgyUtlvSIpBWSLk3lt0paJekxSXMkDUzlknSdpAZJyyVNKrnWTElPpc/MomKuBE+SMDPLp8gWVCMwJSLeBkwEpko6FrgVmAC8FRgCfCLVPwkYnz6zgOsBJO0NXAIcA0wGLpG0V4FxF8otKDOzfApLUJF5Je0OTJ+IiHvSsQAWA6NTnWnALenQA8AISfsBJwILImJDRGwEFgBTi4q7aE5QZmb55EpQkoZIektHLy6pv6RlwPNkSWZRybGBwJnAT1LRKGB1yelrUllr5S3vNUvSEklL1q9f39FQK8YJyswsn3YTlKRTgGWkRCJpoqT5eS4eETsjYiJZK2mypCNKDn8b+N+I+FXTrcpdoo3ylve6ISLqI6K+rq4uT3hV4QRlZpZPnhbUl8nGfjYBRMQyYGxHbhIRm4BfkLrmJF0C1AH/VFJtDTCmZH80sLaN8h7JkyTMzPLJk6B2RMRLHb2wpDpJI9L2EOD9wBOSPkE2rnRGRLxWcsp84Kw0m+9Y4KWIWAfcC5wgaa80OeKEVNYjuQVlZpZPnle+Pybpb4D+ksYDnwJ+k+O8/YCbJfUnS4S3RcTdknYAzwG/lQRwR0TMBu4BTgYagC3AxwEiYoOkrwBN70qfHREbcv+ENcYJyswsnzwJ6h+Bi8mmjX+PrPXylfZOiojlwFFlysveM83qO6+VY3OAOTlirXlOUGZm+bSboCJiC1mCurj4cHq/lmNQTlBmZuW1m6AkLaT8rLkphUTUy7VsQXmShJlZeXm6+C4s2R4MfBjYUUw4vZ+7+MzM8snTxbe0RdGvJf2yoHh6PScoM7N88nTx7V2y2w84GnhzYRH1ck5QZmb55OniW8rrKzrsAJ4Bzi4yqN7MkyTMzPLJ08U3rhKB9BWeJGFmlk+rCUrSh9o6MSLu6P5wej938ZmZ5dNWC+qUNo4F4ATVCU5QZmb5tJqgIuLjlQykr3CCMjPLJ88kCST9FXA42fegAEjr51kHeZKEmVk+ed4H9e/A6WRr8gn4CPAXBcfVa3mShJlZPnlet/GOiDgL2BgRlwLH0fz9TNYBgwYNarbf2NhItk6umZmVypOgmv6Lv0XS/sCrgKeed1K/fv3YY489mpU1NjZWKRozs9qVJ0HdnV48eA3wEPAs2Ws3rJM8UcLMrH15vqjb9O6n2yXdDQzuzBt27XVDhgzh5Zdf3rXvBGVmtrs8kyQekfQFSQdFRKOTU9e5BWVm1r48XXynkq3Bd5ukByVdKOmAguPq1TyTz8ysfe0mqIh4LiKujoijgb8BjiRbMLZNkgZLWpxaYCskXZrKx0laJOkpSd+XtEcqH5T2G9LxsSXXuiiVr5J0Yid/1prhFpSZWfvytKCQNFbS54C5wATgczlOawSmRMTbgInAVEnHAlcB34iI8cBGXl8Z/WyyqewHA99I9ZB0GDCD7IvCU4FvS+qf8+erSf6yrplZ+/KMQS0iW3evH/CRiJgcEV9v77zIvJJ2B6ZPAFOAean8ZuC0tD0t7ZOOv0+SUvncNP71DNAATM7zw9Uqt6DMzNqXZ6mjmRHxRGcunlo6S4GDgW8BTwObIqLplfFrgFFpexSwGiAidkh6CRiZyh8ouWzpOaX3mgXMAjjggNoeInOCMjNrX54xqE4lp3TuzoiYCIwma/UcWq5a+lOtHGutvOW9boiI+oior6ur62zIFeFJEmZm7cs1BtVVEbEJ+AVwLDBCUlPLbTSwNm2vIS2hlI7vCWwoLS9zTo/kFpSZWfsKS1CS6tIKFEgaArwfWAksBKanajOBu9L2/LRPOn5fZIvUzQdmpFl+44DxwOKi4q4ET5IwM2tfu2NQkoYCnwEOiIhzJI0H3hIRd7dz6n7AzWkcqh9wW0TcLelxYK6ky4CHgRtT/RuB/5LUQNZymgEQESsk3QY8TvZ9rPMiYmeHf9Ia4haUmVn78kyS+E+yiQ7Hpf01wA+ANhNURCwHjipT/jvKzMKLiG1kr/Iod63LgctzxNojOEGZmbUvTxffQRFxNdkq5kTEVspPXLCcPEnCzKx9eRLU9jSGFACSDiL7Eq51kltQZmbty9PF92XgJ8AYSbcC7wT+tsCYej1PkjAza1+e1238VNJSsiniAs6PiBcKj6wXcwvKzKx9eWbxzSd7QeH8iPhz8SH1fk5QZmbtyzMG9XXgXcDjkn4gabqkwe2dZK3zJAkzs/bl6eL7JfDL9H2mKcA5wBzgjQXH1mt5DMrMrH15Jkk0rQRxCnA6MInXVx23TnAXn5lZ+/KMQX0fOIZsJt+3gF9ExGtFB9abOUGZmbUv70oSf9PTlxeqJU5QZmbtazVBSZoSEfcBQ4Fp2bsDXxcRdxQcW6/lSRJmZu1rqwX1buA+srGnloLsLbvWCZ4kYWbWvlYTVERckjZnp1et75Jee2Gd5C4+M7P25fke1O1lyuZ1dyB9iROUmVn72hqDmgAcDuwp6UMlh94I+Iu6XeAEZWbWvrbGoN4CfBAYQfNxqM1kX9a1TvIkCTOz9rU1BnUXcJek4yLitxWMqdcr14LasmULLWdKAkjarb6ZWV+QZwzq7yWNaNqRtJekOe2dJGmMpIWSVkpaIen8VD5R0gOSlklaImlyKpek6yQ1SFouaVLJtWZKeip9Znbi56wp/fv3Z+DAgbv2I4I3vOENDB06dLfPkCFDmDBhAo8++mgVIzYzq7w8CerIiNjUtBMRGynzKvcydgCfiYhDyV7VcZ6kw4CrgUsjYiLwpbQPcBIwPn1mAdcDSNobuIRsNYvJwCWS9spx/5rWkVbRqlWr+MIXvlBgNGZmtSdPgupXmhBSwsizyOy6iHgobW8GVgKjyL5D1bTQ7J7A2rQ9DbglMg8AIyTtB5wILIiIDSk5LgCm5vrpatiECRM6VH/lypUFRWJmVpvyLHX0deA3kuaRJZePApd35CaSxpK1uhYBFwD3SvoaWYJ8R6o2ClhdctqaVNZaect7zCJreXHAAQd0JLyquP766zn77LNZtWpV2eMRQWNj4679F198sVKhmZnVhDwtoVskLSF71YaAD0XE43lvIGkY2XepLoiIlyVdBnw6Im6X9FHgRuD96dq73b6N8pZx3gDcAFBfX7/b8Vpz9NFHs2zZslaP79ixo9k41aZNm9i5cyf9+/evRHhmZlWXp4sPYG/gzxHxb8D6vCtJSBpIlpxuLVm7byavL5P0A7JxJchaRmNKTh9N1v3XWnmvNmDAAEaMGNGsbOPGjVWKxsys8tpNUJIuAf4ZuCgVDQT+X47zRNY6WhkR15YcWku2zh9krbKn0vZ84Kw0m+9Y4KWIWAfcC5yQZg/uBZyQynq9kSNHNtt3N5+Z9SV5xqD+mmz8qGnCw1pJw3Oc907gTOBRSU19WV8g+5Lvv0oaAGwjjRsB9wAnAw3AFuDj6X4bJH0FeDDVmx0RG3Lcv8cbOXIkTz/99K59Jygz60vyJKjtERGSAkDSG/JcOCLup/z4EcDRZeoHcF4r15pD9pr5PsUtKDPry/KMQd0m6T/Ipn2fA/wM+E6xYRnA3nvv3WzfCcrM+pI8s/i+JukDwMtk6/N9KSIWFB6Z7daC2rChT/RsmpkB+br4SAnJSanC3MVnZn1Zq118ku5Pf26W9HKZzzOS/qFyofY9TlBm1pe1tZr58enPsjP2JI0EfgN8u5jQzAnKzPqyXF18aWXx48lWcLg/Ih6OiBclvafI4Po6Jygz68vyfFH3S8DNwEhgH+AmSf8C2YKwxYbXtzlBmVlflqcFdQZwVERsA5B0JdmXdi8rMjDzNHMz69vyfA/qWaD05UWDgKfLV7Xu5GnmZtaXtdqCkvRvZGNOjcAKSQvS/geA+ysTXt82bNgwBg4cyKuvvgrA1q1b2bp1K0OGDKlyZGZmxWuri29J+nMp8MOS8l8UFo01I4mRI0fyxz/+cVfZiy++yOjRo6sYlZlZZbQ1zfxmAEmDgYPJWk9PN41FWWU4QZlZX9XWF3UHSLqa7H1MN5O9YmO1pKvTe56sAjyTz8z6qrYmSVxD9qLCcRFxdEQcBRwEjAC+VongzAnKzPquthLUB4FzImJzU0FEvAycS/beJquAllPNPZPPzPqKthJUpHc0tSzcSTYeZRXgFpSZ9VVtJajHJZ3VslDS/wGeKC4kK+UEZWZ9VVvTzM8D7pD0d2RTzQN4OzCE7DXwVgFOUGbWV7XagoqIP0TEMcBsstUkfg/MjojJEfGH9i4saYykhZJWSloh6fySY/8oaVUqv7qk/CJJDenYiSXlU1NZg6TPd/Jn7ZGcoMysr8rzRt37gPs6ce0dwGci4iFJw4GlaTWKNwHTgCMjolHSvgCSDgNmAIcD+wM/k3RIuta3yFawWAM8KGl+RDzeiZh6HCcoM+urcr1uozPSSufr0vZmSSuBUcA5wJUR0ZiOPZ9OmQbMTeXPSGoAJqdjDRHxOwBJc1PdPpGgvGCsmfVVeRaL7TJJY4GjgEXAIcC7JC2S9EtJb0/VRgGrS05bk8paK295j1mSlkhasn79+u7/IarEC8aaWV9VeIKSNAy4HbggfY9qALAXcCzwWeA2SQJU5vRoo7x5QcQNEVEfEfV1dXXdFn+1levicyvKzPqCQhNUWhLpduDWiLgjFa8B7ojMYuA1shchrgHGlJw+GljbRnmfMHDgQA4//PBd+xHB/PnzqxiRmVllFJagUqvoRmBlRFxbcuhOYEqqcwiwB/ACMB+YIWmQpHHAeGAx8CAwXtI4SXuQTaToU/9Cf/jDH262P2/evCpFYmZWOUW2oN4JnAlMkbQsfU4G5gAHSnoMmAvMTK2pFcBtZJMffgKcFxE7I2IH8EngXmAlcFuq22dMnz692f6CBQvYtGlTlaIxM6sMlVnNqMerr6+PJUuWtF+xh4gIJkyYwJNPPrmr7JZbbuHMM8+sYlRmZp0jaWlE1LdXr7Bp5tZ9JDF9+nS++tWv7iq75JJLuPPOOysex/HHH8+5557LoEGDKnpvM+t73ILqIR5++GEmTZpU7TAA+OIXv8js2bOrHYaZ9VB5W1AV+R6Udd3EiRM58MADqx0GAN/97nfpjf+xMbPa4gTVQ0ji61//erXDAGDdunWsWrWq2mGYWS/nMage5LTTTmP16tUsXry44i2Yq6++msWLF+/aX7hwIRMmTKhoDGbWtzhB9TCjR49m9OjRFb9vQ0PDbgnq3HPPrXgcZtZ3uIvPcnnve9/bbH/hwoW89tprVYrGzPoCJyjLZdKkSbzxjW/ctf/CCy+wYkWf+r60mVWYE5TlMmDAAN797nc3K5s3bx6rVq1i1apVrFu3rkqRmVlv5QRlubXs5ps9ezYTJkxgwoQJ7L///px66qm8+uqrVYrOzHobJyjLbcqUKW0e/9GPfsSvfvWrCkVjZr2dE5Tl9ta3vpVDDz20zTrPPfdchaIxs97OCcpy69evH3PnzmXatGkccsghHHLIIYwYMaJZnVdeeaVK0ZlZb+MEZR1y5JFHcuedd+6aHDFr1qxmxzdv3lylyMyst3GCsi4ZPnx4s323oMysuzhBWZcMGzas2b4TlJl1Fyco6xInKDMrihOUdUnLLj6PQZlZdyksQUkaI2mhpJWSVkg6v8XxCyWFpH3SviRdJ6lB0nJJk0rqzpT0VPrMLCpm6zi3oMysKEWuZr4D+ExEPCRpOLBU0oKIeFzSGOADwO9L6p8EjE+fY4DrgWMk7Q1cAtQDka4zPyI2Fhi75eQEZWZFKawFFRHrIuKhtL0ZWAmMSoe/AXyOLOE0mQbcEpkHgBGS9gNOBBZExIaUlBYAU4uK2zrGXXxmVpSKjEFJGgscBSySdCrwh4h4pEW1UcDqkv01qay18pb3mCVpiaQl69ev78borS1uQZlZUQpPUJKGAbcDF5B1+10MfKlc1TJl0UZ584KIGyKiPiLq6+rquhCxdYQTlJkVpdAEJWkgWXK6NSLuAA4CxgGPSHoWGA08JOnNZC2jMSWnjwbWtlFuNcBdfGZWlCJn8Qm4EVgZEdcCRMSjEbFvRIyNiLFkyWdSRPwRmA+clWbzHQu8FBHrgHuBEyTtJWkv4IRUZjVg6NChzfa3bNnCzp07qxSNmfUmRbag3gmcCUyRtCx9Tm6j/j3A74AG4DvAPwBExAbgK8CD6TM7lVkN6N+/f9kkZWbWVYVNM4+I+yk/flRaZ2zJdgDntVJvDjCnO+Oz7jN8+PBmSWnz5s27df2ZmXWUV5KwLvNECTMrghOUdZkTlJkVwQnKuqxlgvJMPjPrDk5Q1mV+J5SZFcEJyrrMXXxmVgQnKOsyd/GZWRGcoKzL3MVnZkVwgrIucxefmRXBCcq6zF18ZlYEJyjrMnfxmVkRnKCsy9zFZ2ZFcIKyLnMXn5kVwQnKusxdfGZWBCco6zJ38ZlZEZygrMucoMysCE5Q1mV+7buZFcEJyrrMLSgzK0JhCUrSGEkLJa2UtELS+an8GklPSFou6YeSRpScc5GkBkmrJJ1YUj41lTVI+nxRMVvnOEGZWRGKbEHtAD4TEYcCxwLnSToMWAAcERFHAk8CFwGkYzOAw4GpwLcl9ZfUH/gWcBJwGHBGqms1YujQoUjatb9t2zZ27NhRxYjMrDcoLEFFxLqIeChtbwZWAqMi4qcR0fSv1wPA6LQ9DZgbEY0R8QzQAExOn4aI+F1EbAfmprpWIyS5FWVm3a4iY1CSxgJHAYtaHPo74MdpexSwuuTYmlTWWnnLe8yStETSkvXr13dP4JabE5SZdbfCE5SkYcDtwAUR8XJJ+cVk3YC3NhWVOT3aKG9eEHFDRNRHRH1dXV3XA7cO8Uw+M+tuA4q8uKSBZMnp1oi4o6R8JvBB4H0R0ZRs1gBjSk4fDaxN262VW41wC8rMuluRs/gE3AisjIhrS8qnAv8MnBoRW0pOmQ/MkDRI0jhgPLAYeBAYL2mcpD3IJlLMLypu6xwnKDPrbkW2oN4JnAk8KmlZKvsCcB0wCFiQZn49EBF/HxErJN0GPE7W9XdeROwEkPRJ4F6gPzAnIlYUGLd1grv4zKy7FZagIuJ+yo8f3dPGOZcDl5cpv6et86z63IIys+7mlSSsWzhBmVl3K3SShPUdLRPUNddcw9y5c6sUjZkV7ZhjjuGKK64o9B5OUNYtWo5BNTQ00NDQUKVozKxoQ4YMKfwe7uKzbrH//vtXOwQz62WcoKxbTJ8+nYMOOqjaYZhZL+IuPusWdXV1PPbYYzz88MNs3bq12uGYWcFGjhxZ+D2coKzbDB48mOOOO67aYZhZL+EuPjMzq0lOUGZmVpOcoMzMrCY5QZmZWU1ygjIzs5rkBGVmZjXJCcrMzGqSXn+hbe8haT3wXCdP3wd4oRvDKVJPihUcb9F6Urw9KVZwvN3tLyKirr1KvTJBdYWkJRFRX+048uhJsYLjLVpPircnxQqOt1rcxWdmZjXJCcrMzGqSE9Tubqh2AB3Qk2IFx1u0nhRvT4oVHG9VeAzKzMxqkltQZmZWk5ygzMysJjlBJZKmSlolqUHS56sdT0uSxkhaKGmlpBWSzk/le0taIOmp9Ode1Y61iaT+kh6WdHfaHydpUYr1+5L2qHaMTSSNkDRP0hPpGR9X48/20+nvwWOSvidpcC09X0lzJD0v6bGSsrLPU5nr0u/eckmTaiTea9Lfh+WSfihpRMmxi1K8qySdWO1YS45dKCkk7ZP2q/5su8IJiuwfUuBbwEnAYcAZkg6rblS72QF8JiIOBY4Fzksxfh74eUSMB36e9mvF+cDKkv2rgG+kWDcCZ1clqvL+FfhJREwA3kYWd00+W0mjgE8B9RFxBNAfmEFtPd+bgKktylp7nicB49NnFnB9hWIsdRO7x7sAOCIijgSeBC4CSL93M4DD0znfTv+GVMpN7B4rksYAHwB+X1JcC8+205ygMpOBhoj4XURsB+YC06ocUzMRsS4iHkrbm8n+AR1FFufNqdrNwGnVibA5SaOBvwK+m/YFTAHmpSq1FOsbgb8EbgSIiO0RsYkafbbJAGCIpAHAUGAdNfR8I+J/gQ0tilt7ntOAWyLzADBC0n6ViTRTLt6I+GlE7Ei7DwCj0/Y0YG5ENEbEM0AD2b8hVYs1+QbwOaB05lvVn21XOEFlRgGrS/bXpLKaJGkscBSTAluUAAAF70lEQVSwCHhTRKyDLIkB+1Yvsma+SfbL8lraHwlsKvmFr6VnfCCwHvjP1CX5XUlvoEafbUT8Afga2f+U1wEvAUup3efbpLXn2RN+//4O+HHarrl4JZ0K/CEiHmlxqOZi7QgnqIzKlNXk/HtJw4DbgQsi4uVqx1OOpA8Cz0fE0tLiMlVr5RkPACYB10fEUcCfqZHuvHLS2M00YBywP/AGsq6clmrl+banlv9uIOlisi72W5uKylSrWryShgIXA18qd7hMWc082/Y4QWXWAGNK9kcDa6sUS6skDSRLTrdGxB2p+E9NTfb05/PViq/EO4FTJT1L1l06haxFNSJ1SUFtPeM1wJqIWJT255ElrFp8tgDvB56JiPUR8SpwB/AOavf5Nmntedbs75+kmcAHgY/F618arbV4DyL7z8oj6XduNPCQpDdTe7F2iBNU5kFgfJoFtQfZAOj8KsfUTBrDuRFYGRHXlhyaD8xM2zOBuyodW0sRcVFEjI6IsWTP8r6I+BiwEJieqtVErAAR8UdgtaS3pKL3AY9Tg882+T1wrKSh6e9FU7w1+XxLtPY85wNnpRlnxwIvNXUFVpOkqcA/A6dGxJaSQ/OBGZIGSRpHNgFhcTViBIiIRyNi34gYm37n1gCT0t/rmny2uUWEP9l/jE4mm6nzNHBxteMpE9/xZE3z5cCy9DmZbGzn58BT6c+9qx1ri7jfA9ydtg8k+0VuAH4ADKp2fCVxTgSWpOd7J7BXLT9b4FLgCeAx4L+AQbX0fIHvkY2PvUr2D+bZrT1Psm6ob6XfvUfJZifWQrwNZOM3Tb9v/15S/+IU7yrgpGrH2uL4s8A+tfJsu/LxUkdmZlaT3MVnZmY1yQnKzMxqkhOUmZnVJCcoMzOrSU5QZmZWk5ygzFqQdIWk90g6TW2sbC/prLSa+ApJj0u6sJ3rfrmpjqSbJE1vq37Jee+R9I4OxD9W0hpJ/VqUL5M0uUwcz6RjD0k6rkz5Mkm/SeV/K2l9KntC0qfzxmXWUU5QZrs7hmydw3cDvypXQdJJwAXACRFxONnKEy8VFM97yFaKyCUiniX7/s67msokTQCGR0S5L5R+NiImki3v9B8ty9On9P7fT/XfCVycVtE263ZOUGZJev/PcuDtwG+BTwDXSyq3xtlFwIURsRYgIrZFxHfSdc6R9KCkRyTdntZKa+u+V6YW2HJJX2txbCzw98CnU6vlXZL+QtLPU/2fSzqgzGW/R7aKR5MZqawt/wsc3E6dXSLiRbIvs/aY1bGtZ3GCMksi4rNkSekmsiS1PCKOjIjZZaofQbaCeDl3RMTbI6LpvVKtvpdJ0t7AXwOHR/beoctaxPQs8O9k73maGBG/Av4v2SsUjiRbwPS6Mpe+DTitZG2+08nWRWzLKWSrDTS5pqSL79aWlVNiHEy2+oZZtxvQfhWzPuUosmVtJpCtb9cZR0i6DBgBDAPubaPuy8A24LuS/ge4O8f1jwM+lLb/C7i6ZYWI+KOkFcD7JP0JeDUidnsDa3KNpH8he+VIaTL9bETMK1P/dEnvBd4CnBMR23LEbNZhTlBmgKSJZC2n0cALZC8BlKRlwHERsbXFKSuAo4H7ylzuJuC0iHhE0t+SjSGVFRE7JE0mW/B1BvBJstXfO6K19cqauvn+RNvde60lotZ8PyI+mSZU/I+kH0e2MKlZt3IXnxkQEcvSwP+TwGFkiefE1K3WMjkBXAFcnV5pQFrZ+lPp2HBgXXo9ysfaum96v9eeEXEP2aSLiWWqbU7XbPIbXh9f+hhwfyuXv51sQeE83XsdFhG/JWvBnd/d1zYDt6DMdpFUB2yMiNckTYiIVrv4IuIeSW8CfpZeeRHAnHT4i2SzAJ8jG9MZXv4qkI7dJWkw2crT5aZt/wiYJ2ka8I/Ap4A5kj5L1i338VZi3CTpAbI32T7TRgytaer6a1LuteZXkb176KsRsbkT9zBrlVczNzOzmuQuPjMzq0lOUGZmVpOcoMzMrCY5QZmZWU1ygjIzs5rkBGVmZjXJCcrMzGrS/wdePnfLlHvwiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Find out the suggestion made by ESOP\n",
    "IPLBest = getWinner()\n",
    "print( \"Start on day %d, last for %d days, at level %f\" % ( IPLBest[0], IPLBest[1], IPLBest[2]/2 ) )\n",
    "\n",
    "validVals = np.array( list( globalDict.values() ) )\n",
    "mins = np.minimum.accumulate( validVals[:,2] )\n",
    "plt.figure()\n",
    "plt.plot( np.arange(1, len( validVals ) + 1), mins, color = 'k', linewidth = 4 )\n",
    "plt.xlabel( \"# Calls to VIPER\" )\n",
    "plt.ylabel( \"Objective value\" )\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10849 Infected\n",
      "0 Expired\n",
      "Peak 1653\n",
      "Recovery time 17.019725\n",
      "Expiry time nan\n",
      "Quarantine time 12.008547\n",
      "117 Quarantined\n",
      "Objective acheived is 1653.000000 (fEpi) + 480.000000 (fEco) = 2133.000000 (overall)\n"
     ]
    }
   ],
   "source": [
    "# Get hold of outcomes had no lock-down been applied at all\n",
    "pop.reset()\n",
    "np.random.seed(0)\n",
    "\n",
    "LKP = np.zeros((T,))\n",
    "(statsNoLockDown, tInfect, tQuarantine, tRecovery, tExpiry) = pop.simulate( T = T, LKP = LKP )\n",
    "\n",
    "# Now get hold of outcomes with the lock-down suggested by ESOP\n",
    "pop.reset()\n",
    "np.random.seed(0)\n",
    "\n",
    "LKP = getLKP( IPLBest )\n",
    "\n",
    "(stats, tInfect, tQuarantine, tRecovery, tExpiry) = pop.simulate( T = T, LKP = LKP )\n",
    "\n",
    "outcome = obj( stats, IPLBest )\n",
    "print( \"%d Infected\" % len(tInfect) )\n",
    "print( \"%d Expired\" % len(tExpiry) )\n",
    "print( \"Peak %d\" % max(stats[SIDX_EI,:]) )\n",
    "print( \"Recovery time %f\" % np.mean( tRecovery ) )\n",
    "print( \"Expiry time %f\" % np.mean( tExpiry ) )\n",
    "print( \"Quarantine time %f\" % np.mean( tQuarantine ) )\n",
    "print( \"%d Quarantined\" % len( tQuarantine ) )\n",
    "print( \"Objective acheived is %f (fEpi) + %f (fEco) = %f (overall)\" % (outcome[0], outcome[1], outcome[2]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Outcome of Lock-down suggested by ESOP**: The following code is meant to reproduce Fig 7(b) in the accompanying paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXlYldXWwH+bSZRBUBRFCRzRNLE0bTLRq3Yz9XazrNSctTJzyr66paVlVmZm2e2WqWnODVqi2aCFQ6XZgJZizoigAioKKDKt74/3cOAgw8t0DuD+Pc/7eN5377322kfOWWfvvfZaSkTQaDQajaay4eRoBTQajUajKQhtoDQajUZTKdEGSqPRaDSVEm2gNBqNRlMp0QZKo9FoNJUSbaA0Go1GUylxcbQCGo1Go6ncKKWOA8lAFpApIh3t0a82UBqNRqMxQzcRSbRnh3qJT6PRaDSVEnUtR5JwcnKSmjVrOlqNcic7Oxsnp+r326M6jqs6jgn0uKoSly5dEuD3PI8WiMiCvHWUUseA84AAH+QvryiuaQPl7u4uaWlpjlaj3ImIiCAsLMzRapQ71XFc1XFMoMdVlVBKXRIRj2LqBIhInFKqPvAd8KSIbKto3arXTwGNRqPRlDsiEmf5Nx5YB3SyR7/aQGk0Go2mUJRSHkopr5zXQC/gL3v0rb34NBqNRlMU/sA6pRQYNmOliHxtj461gdJoNBpNoYjIUSDUEX1rJwntJOFwsrOzOXnyJKmpqUXWS0tLw93d3U5a2YfqOCbQ46pseHh40Lhx4wI9EM04STgKPYPSOJzExESUUoSEhBTpwpucnIyXl5cdNat4quOYQI+rMpGdnU1sbCyJiYnUr1/f0eqUCO0koXE4SUlJ+Pv7V7vzJRpNZcDJyQl/f38uXLjgaFVKjP5GuNaIj4INk+HsEUdrYiUrKwtXV1dHq6HRVFtcXV3JzMx0tBolRhuoa4mUeFgzGH5dBP+7HeL+cLRGViweQhqNpgKoqp8vbaCuFS6fh2X3wdnDxn3mZdjysmN10mg0miLQBupa4dBmOPOn7bP0VMi47Bh9NBXOkiVLuOOOOxyqw7p16wgMDMTT05M//ii/Gfv27dsJCQkpN3lmWbx4MRMnTrR7vzmsX7+ehx56yGH92xttoK4V2j0AQ9YbrwM7w8jNMPIbcK1+wXLLk+DgYGrWrImnp6f1GjdunKPVchgRERE0btzYdP0pU6bw7rvvkpKSwo033ljqfpVSHD582HrfpUsX/v7771LLKw3p6em88cYbPP300wAcP34cpZTN34anpydr1qwp136Dg4PZvHkzAP369eOvv/5i79695dpHZUUbqGpIg1PfwVf/BxfjbAuadoVpiTDyWwi82THKVUHCw8NJSUmxXu+++66jVSp3wsLCiIiIKHe50dHRtGnTptzlOoIvv/ySli1b0qhRI5vnSUlJNn8fDz74YLGylixZwrBhw0qlx8MPP8yCBXYJJu5wtIGqbmReIfj4KvjlA3i7PXz1NKQk5JY7V2Jvuem1i7y83mxcbJ0SXWXk8ccf5/7777feP/PMM/zjH/9ARKwzjVmzZuHn50dwcDArVqyw1r1w4QJDhgyhSZMmBAUFMXPmTLKzswE4fPgwXbt2pXbt2vj5+dl84R04cICePXtSp04dQkJC+OSTT6xlZ8+epV+/fnh7e9OpUyeOHKl4T82wsDCmTZvG7bffjpeXF7169SIxMZErV67g6elJVlYWoaGhNGvWDIC4uDj69+9PvXr1aNKkCe+8845VVlZWFrNmzaJZs2Z4eXnRoUMHYmJiuPPOOwEIDQ21zlDyz+SioqIICwvDx8eHNm3asH79ehsdFy5caL3Pu/QpIkyaNIn69etTu3Zt2rVrx19/FRxmbtOmTdx+++3l9+aVkrCwMDZu3OhoNeyCNlDVjT+W4X7lrPE66wr8/jFIlmN1qqa8+eab7N27lyVLlrB9+3YWLVrE0qVLrR5Tp0+fJjExkdjYWJYuXcqYMWOsy1JPPvkkFy5cYO/evWzdupWPP/6Yjz76CIBp06bRq1cvzp8/z8mTJ3nyyScBSE1NpWfPngwcOJD4+HhWrVrF2LFj2bdvHwBPPPEE7u7unDp1isWLF7N48WK7vA8rV67ko48+Ij4+nvT0dObMmUONGjVISUkBYM+ePRw5coTs7Gz69u1LaGgosbGxbNmyhXnz5vHNN98AMHfuXFatWsVXX33FxYsXWbx4MbVq1WLbtm1WOQXNUDIyMujbty+9evUiPj6e+fPnM2jQIFNLgN9++y3btm3j4MGDJCUlsWbNGurWrVtg3T///JMWLVqU5a0qF1q3bs3x48e5ePGio1WpcLSBqk6IwE/5lp86DAevBo7Rp5pw77334uPjY70+/PBDAGrVqsXy5cuZPHkygwcPZv78+Vftz7z88svUqFGDrl27cs899/DJJ5+QlZXFmjVrePXVV/Hy8iI4OJinnnqKZcuWAcaZlejoaOLi4nB3d7f+2t+wYQPBwcEMHz4cFxcXbrrpJvr3789nn31GVlYWn3/+OS+99BIeHh60bduWoUOH2uX9GT58OC1btqRmzZoMGDCAyMjIAuvt3r2bhIQEXnjhBdzc3GjatCmjR49m9erVACxcuJCZM2cSEhKCUorQ0NBCjUVedu7cSUpKCs8++yxubm50796dPn36sGrVqmLburq6kpyczIEDBxARWrduTcOGDQusm5SUhKen51XP/fz8bP4+oqKiiu23LOREskhKSqrQfioD2kBVNx5YwtEmg+G628DVA+4w4XEkAgkHISuj4vWrgnzxxRckJSVZr9GjR1vLOnXqRNOmTRERBgwYYNPO19cXD4/cEGdBQUHExcWRmJhIeno6QUFBNmWxsbEAzJ49GxGhU6dOtGnTxjoTio6OZteuXTZfhitWrOD06dMkJCSQmZlJYGCgjcyiyCtnx44d9OnTx3r/2muvmX5/GjTI/QFUq1Yt68wpPzlGN2+/s2bN4syZMwDExMRYlwJLQlxcHIGBgTaRSPK+n0XRvXt3xo0bxxNPPIG/vz9jxowpdGbi6+tb4NgSExNt/j5at25dYPuxY8daxz127FhWrlxpvW/Xrp3J0RrhlsD4/6vu6Fh81QmlIKA9J4IeoGlYmOFCXpSX3tEIiAqHQ99C0gkYthGCHeiWPL3oUCyVMQ7af//7X65cuUJAQACzZ8/mP//5j7Xs/PnzpKamWo3UiRMnaNu2LX5+ftZZUo5BOXHihHXzvUGDBtZZ2o4dO+jRowd33nkngYGBdO3ale++++4qPbKysnBxcSEmJoZWrVpZZRZF3l/gYWFhTJ8+vUKDDAcGBtKkSRMOHTpUaPmRI0do27ZtieQGBAQQExNjk479xIkTtGzZEjACpV66dMla//Tp0zbtx48fz/jx44mPj2fAgAG88cYbvPzy1WcE27VrZ+NJWFLee+893nvvPcDYB4uIiGDJkiUllhMVFUVwcDDe3t6l1qWqoGdQ1ZniXMgjV8HuhYZxAjh09RefpnAOHjzI1KlTWb58OcuWLWP27NlXLW+9+OKLpKens337djZs2MADDzyAs7MzAwYM4Pnnnyc5OZno6Gjmzp3L4MGDAfj00085efIkYPxqV0rh7OxMnz59OHjwIMuWLSMjI4OMjAx2795NVFQUzs7O3HfffUyfPp1Lly6xf/9+li5davf3pCg6deqEt7c3r7/+OpcvXyYrK4u//vqL3bt3AzBq1CimTZvGoUOHEBH27t3L2bPGfqq/vz9Hjx4tUG7nzp3x8PBg9uzZZGRkEBERQXh4uPW8UPv27Vm7di2XLl3i8OHDLFq0yNp29+7d7Nq1i4yMDDw8PHB3d8fZ2bnAfnr37s2OHTvK8y0pFVu3buXuu+92tBp2QRuoa5kWPW3vD292jB6VnL59+9qcc/n3v/9NZmYmgwcP5plnniE0NJQWLVowa9YsHnnkEa5cuQIYMyFfX18CAgIYNGgQ77//vnV2M3/+fDw8PGjXrh133HEHAwcOZMSIEYDxpdm5c2c8PT3p168fb7/9Nk2aNMHLy4tvv/2W1atXExAQQIMGDXjmmWes/eWcN2rQoAHDhg1j+PDhjnnDCsHZ2Znw8HAiIyNp0qQJfn5+jBo1yhrEdPLkyQwYMIBevXrh7e3NyJEjuXzZOEg+ffp0hg4dio+Pj43nIoCbmxvr169n06ZN+Pn5MXbsWD7++GPrez1p0iTc3Nzw9/dn6NChDBo0yNr24sWLjB49Gl9fX4KCgqhbty5TpkwpUP++ffty6NAh4uJsj2/4+PjY/H3MnTu33N6zgli1ahWPPvpohfZRWdD5oK7lfFCXzsEbzUCywbmGsbz38CpwqVHhOuYlKiqq0HX7vFTGJb7CiIiIYPDgwdaZUGFUpTGVhOo6rnfeeYejR48yb948h/QfHh7OsmXLrjLSZijsc6bzQWkqJ7XqQI/p4BcCTbqAW6X8G9VoKg3Dhw93qOHt27cvffv2dVj/9kYbqOrCj+9AYCdo1LFk7W6fUDH6aDQaTRnRBqo6kHQCvptmvK5Rm9a120PXroZXn8YhhIWFFbu8p9FoisaUk4RSyk8pdV2+Z48qpeYrpfqUpmOlVAulVJpSanmeZwOVUtFKqVSl1BdKqTp5yuoopdZZyqKVUgPzySu0bbUn+qfc11cuUONKgjZOGo2mymPWi28x8GzOjVJqGvA/YCDwpVKq+OiIV/NfYHcemW2AD4BHAH/gEvBevvrplrJBwP8sbcy0rd7kNVDAhdrVIzinRqO5tjFroDoCW/LcPwbMEpG6GIZjckk6VUo9BCTlkzkICBeRbSKSAkwD7lNKeSmlPID+wDQRSRGRHcB6DINUZNuS6FVlCe4C1/8LPOoBkOSjDZRGo6n6mN2DqgOcAVBKtQUaADmnAL8AhpjtUCnlDbwE/AMYmaeoDWCdCojIEaVUOtASyAayRORgnvp7gK4m2v6Wr/8xwBgAFxeXCkkxYH/qQf0RUG84NS/HcjazFudLOK7aSX/hffEgXsmH8Uo+wm8d3iTT9eq4YxVB7dq1reFbiiIrK8tUvapEdRwT6HFVRtLS0qrc951ZA3UWyImC2R2IE5GceCWulOzA78vAIhGJUbb7JJ5A/lg3FwAvIKuIsuLa2iAiC4AFYJyDqsjQLo7C9DmovLz3HMTvs97e0cwTmpZQRimJiooy5bpbHc/WVMcxgR5XZcTd3b1MSSMdgVnDshmYrpQaBzyFMWvKoRUQbUaIUqo90AN4q4DiFCB/cClvILmYsuLaaswQkO8PN6780nNrruaxxx4rMN5baRk2bBhTp04ttDwhIYGQkBDMHkx/5513ePbZZ4uvqNFUIGYN1P8BMcCrwBFgRp6yQYDZAFVhQDBwQil1GpgC9FdK/Q7sA0JzKiqlmgI1gIOWy0UplTcZS6ilDcW01ZghoL3tvTZQQMWlfH///feZNm1aOWhojtdee43hw4fj7u5OVlYWnTp1YtasWdbyrKwsOnbsyJw5cwAYM2YMy5cvJz4+3m46ajT5MbXEJyJngJ6FFPcAzMYLWgCsznM/BcNgPQ7UB35WSnUBfsfYp1orIskASqm1wEtKqVFAe+BfwG0WOSuKaqsxQdBtEDrQmEkF3AgNShZRujoTHh5Ojx497NZfZmYmLi7ld0TxypUrLF261BrI1tnZmcWLF3Pbbbdx33330apVK+bMmYNSikmTJgHGctDdd9/Nxx9/XGhsOo2moinzp0BETKd1FJFLGC7gACilUoA0EUkAEpRSj2EYm7oYy4p5o12OxXB3j8fYE3tcRPZZ5O4rpm315MJJI4aeZ72yy/JvA//+X9nllIGoVsXH4ytPWh8oW2K5xx9/nISEBD777DPASPn+66+/snnzZrZu3crgwYMZO3Ysc+fOxdPTk1deecUaqHTYsGE0btyYZ555xhq378knn+Stt96iZ8+eLFu2jA0bNjB16lSOHz/O9ddfz/vvv2/NG/THH38wcuRIDh06RO/evVFFnHvLySGVN5li27ZtmTx5MiNHjmThwoXMmjWLHTt22ETyzkmVrg2UxlEUaqCUUiXJFy0iMrL4alc1mp7vfiWwspC654B7i5BVaNtqy7Y58NsSY8bT8i5o9yDUaeJora4Z3nzzTdq3b8+SJUto1qwZixYtIjIyssCU7zt37qR379507NiRkJCQq2SdPn2ac+fOER0dTXZ2Nr///jsjRowgPDycjh07snz5cvr168fff/+NUop7772XiRMnMm7cOL788ksefvhhnnnmmQL1/PPPPwvs87nnnuPzzz+nS5cuTJw4kRtuuMGmvHXr1uzZs6cc3imNpnQUtQfVHehWgktjT0SMRIMIxP0OEa/C2SOO1qpaUt4p3wvCycmJGTNmUKNGDWrWrMmHH37Io48+SufOnXF2dmbo0KHUqFGDnTt3snPnTjIyMpg4cSKurq7cf//93HzzzYXqn5SUVKDnmZubG507d+bs2bM2KShy8PLysqbC0GgcQaEzKBEJtqMempISvx8u5klp7VITgm93nD7VmC+++KLQPaiclO852VjzUljK94KoV68e7u7u1vvo6GiWLl3K/Pnzrc/S09OJi4tDKUWjRo1slvWKSu/u6+tb4Nmd7du388UXXzBkyBAmTJjApk2bbMqTk5OpXbt2oXI1mopGB4utqqScAd8mcP6Ycd/kzuIz6FZyitsTqoxnUEqT8r0g8u8hBQYG8vzzz/P8889fVXfr1q3ExsYiItZ2J06coFmzZgXKbteuHW+9ZXuyIy0tjZEjRzJnzhwGDBhA27ZtWb58uTWrLxjn00JDQ/OL02jsRokz6iql6iulrst/VYRymiJo1h0mRMKUQ/DgCrh9fPnJTrsIx7YbKTxSz5af3GpGaVO+m2H06NG8//777Nq1CxEhNTWVjRs3kpyczK233oqLiwvvvPMOmZmZrF27ll9++aVQWZ06dSIpKYnY2NwZ9wsvvEBQUBDDhg2jVq1aLFiwgEmTJpGQkGCtcy2lFtdUTsxGM3dSSs1SSp0FTgHHCrg0jsCzPrTuY2TDLQ8+GQqvBcLSPkYKj9hfy0duFaYiUr4XR8eOHfnwww8ZN24cvr6+NG/enCVLlgDG3tHatWtZsmQJvr6+rFmzhvvuu69QWW5ubgwbNozly43EAb/++isffPABCxYssNbp0aMHffr0YeLEiYAxw/rqq68YOnRoad4yjaZcMJXyXSk1GXgReB2YCbyCER9vkOXf10SkJF5/lYJrPuV7QWx8CnYvzL0P+w+EVWxEAZ3yveLHlJCQQJcuXfjjjz+oWbP4peD58+cTExPD7NmzS9VfVfq/KglVeVzVOeX7cIzDr/MwDNQ6EfldKTUT+BbQS3zVBR3yqFpSr149Dhw4YLr+k08+WYHaaKoiSiln4FcgVkRKlQewpJjdg2oK/CoiWUAmUBNARDIwjNaIilFPY3cCbgQU1GsFoQ9D636O1kij0VQOJgBlO91eQszOoC4AOT6wcUAI8GMeGddO9trqTr1W8J+TUMM+qTaqKzrlu6Y6oZRqDNyDsb1Tovx/ZcGsgfoDuB74xnLNUEpdxphNvYIR/05jDzIuw8Ke0OgmCOxsXH7Ny0++k7M2ThrNtYWLUiqvN9QCS1qivMzDCBpu1w04swZqHsYyHxjOEjdhxL0DI9VG2cM7a8wR+zuc+dO4fl8KvsEwQYej0Wg0pSZTRDoWVqiU6gPEi8hvSqkw+6llPpr5d3len1ZKdQKaAbWAKMtelMYexP5mex/Y2TF6aDSaa4XbgX5Kqd4YWz3eSqnlIjK4mHZlpsQHdcGIDCsih0VkrzZOdibxb9v7hu0LrqfRaDTlgIj8R0QaW8LfPQR8bw/jBCZnUEqpO4urIyLbyq6Opli6TIEWvSDxICQchMaFzszLj6xMY2+qiJQOGo1GU96Y3YOKAIo70etcTLmmPKjTxD4pNbbPhVN7DEN49jBM2mdErdBoKoB169Yxfvx4zp8/z/bt27nxxhuLb2SC7du3M2rUKP7+++/iK5cjH3zwAVFRUcybN88u/Z05c4awsDAiIyOpUaNGhfUjIhEY9sAumF3i64aRfiPv9QCwFDgO2OXQlsaO7FsH+78woqZnpUOCfT/glYXg4GD8/f1JTU21Plu4cGHpI3VcI2zfvv2q1CNFMWXKFN59911SUlLKZJyUUhw+fNh636VLF7sbp/T0dGbOnMnTTz8NwPHjx1FKkZmZWWqZ+/fvp1+/ftSuXRsvLy+6d+/Ozp07reX+/v5069bNJnyVPVFKfV+Ca4tZuaYMlIhsLeBaKyIjgPVA39IOTFNJqZcvwV3+va9riMzMTN5++21Hq1GulOXLsiKIjo6mTZs2jlajXPjyyy9p1aoVjRo1Khd5R44c4fbbb+eGG27g2LFjxMXFce+999KzZ0+bIMGDBg3igw8+KJc+S4EToExepn0fyiPdxkZgNUZKdk11wS+fgUqKqfAu//vY9xXeR16eeL+7qXpPP/00s2fPZuzYsfj4+FxV/tNPPzFhwgQOHjxIy5Ytefvtt7ntttsA48Buly5d+P7779m7dy+33norK1euxM/Pr8C+jh07xrBhw/j999+55ZZbCAkJISkpieXLlxcY3y84OJiFCxfSo0cPfvnlFyZMmEBUVBQ1a9akf//+zJ07Fzc3N8CYXbz77rvMmzePzMxMjh07xoQJE1i7di0XLlygRYsWzJs3jy5dugAwffp09u/fj7u7O+vWreO6665j6dKldOxY8n3Pwt4HLy8v6tatS1ZWFqGhoTRo0IAjR44QFxfHk08+ybZt2/D09GTSpEmMH29E7M/KyuL1119n0aJFxMfH07JlS7744gtr0sXQ0FCUUixatAh/f3+b9ywqKorHH3+cyMhIGjVqxKuvvkq/fv2sOg4ePJhRo0YBsGTJEhYuXMiOHTsQESZPnszy5ctJT08nKCiIlStXFpg+ZdOmTXTt2rXE71FhTJ8+nVtvvZVXXnnF+mz8+PFERUXxzDPP8MMPPwDQuXNnjh49SnR0dJH5wSoCEQmrCLml8uLLRwhGwFhNdSLkbug3H0Z8C/93DHrOcLRGDqNjx46EhYUxZ86cq8rOnTvHPffcw/jx4zl79iyTJ0/mnnvu4ezZ3DQlK1eu5KOPPiI+Pp709PQC5eQwcOBAOnToQGJiItOmTWPp0qWm9XR2duatt94iMTGRn3/+mS1btvDee+/Z1Pniiy/YtWsX+/fvB+Dmm28mMjKSc+fOMXDgQB544AHyBlBev349Dz30EElJSfTr149x40p/5LGg96FGjRqkpKQAsGfPHo4cOUJ2djZ9+/YlNDSU2NhYtmzZwrx58/jmm28AmDt3LqtWreKrr77i4sWLLF68mFq1arFt2zarnJSUFB588EGb/jMyMujbty+9evUiPj6e+fPnM2jQIFNLgN9++y3btm3j999/JykpiTVr1lC3bt0C6/7555+EhIQUWFYavvvuuwLTtAwYMIDt27db/79cXFxo3rw5e/ZUn3ORZtNtDCngGqWUmge8BmwqToamHDi6FeKjIMsOnv0N2sJNQ+C6zlBLR7J66aWXmD9/vk2+JICNGzfSokULHnnkEVxcXHj44Ydp1aoV4eHh1jrDhw+nZcuW1KxZkwEDBlyVMyqHEydOsHv3bmua+DvvvJO+fc2vnnfo0IFbbrkFFxcXgoODefTRR9m6datNnf/85z/UqVPHGtF88ODB1K1bFxcXF5566imuXLli84V9xx130Lt3b5ydnXnkkUfK9OVn9n3YvXs3CQkJvPDCC7i5udG0aVNGjx7N6tWrAWMPcObMmYSEhKCUIjQ0tFBjkZedO3eSkpLCs88+i5ubG927d6dPnz6sWrWq2Laurq4kJydz8OBBRITWrVvTsGHDAusmJSWVa8TzxMTEAvtq2LAhWVlZnDt3zvrMy8uLpKSkcuu7tCilGiml5iqlflVKHVNKtbU8n6iUMn140+wMakkB1wLgUeBzoByz5WkK5bMR8N4tMCsA/nc7pMQ7WqNrhrZt29KnTx9ee+01m+dxcXFXLacEBQXZJAds0KCB9XWtWrWsM4bHHnuMhg0b4unpyaxZs4iLiyswTbxZDh48SJ8+fWjQoAHe3t4899xzJCYm2tQJDAy0uX/zzTdp3bo1tWvXxsfHhwsXLti0ya97WlpaqfevCnsf8hMdHU1cXBw+Pj7Wa9asWZw5cwaAmJiYQrMHF0VcXByBgYE4OeV+7eX/vyqM7t27M27cOKZMmYK/vz9jxozh4sWLBdb19fUlOTnZlE7bt2+35hkrbA/Oz8+PU6dOXfX81KlTKKVslouTk5MLXIa2J0qpNsCfwCMYsVuvA9wsxUEYQWdNYXYPqiC/5jQROWO2I00ZSUmAS5Yvjqx0w/W7VvG/GqsSxe0JOToXz4wZM7jpppt46qmnrM8CAgKIjo62qXfixAn++c9/Fivv/fff54033rCOKTo6usA08Tlp3T08PLh06ZK1fVZWls2M7vHHH+fGG29k1apVeHl5MW/ePD777DObPvOmlt++fTuvv/46W7ZsoU2bNjg5OeHr64uZHHEVSWBgIE2aNOHQoUOFlh85cqTA/Z+iCAgIICYmhuzsbKuROnHiBC1btgSufn9Pnz5t0378+PEMHz6cy5cvM2DAAN544w1efvnlq/pp164dBw8eNKVTly5dCjXUOfTo0YNPP/2U4cOH2zz/5JNPuOWWW6x7jJmZmRw+fJjQ0FBTfVcgb2JEPb8LSAPS85T9hJFX0BRmvfiiC7i0cbIn8fts7+uFGIdnNXajefPmPPjgg7zzzjvWZ7179+bgwYOsXLmSzMxM1qxZw/79++nTp+QnL4KCgujYsaM1TfyOHTtslgpbtmxJWloaGzduJCMjg5kzZ1qz94JhwL29vfH09OTAgQP873//K7K/5ORkXFxcqFevHpmZmbz00kuFzgrsSadOnfD29ub111/n8uXLZGVl8ddff7F7924ARo0axbRp0zh06BAiwt69e617fv7+/hw9erRAuZ07d8bDw4PZs2eTkZFBREQE4eHhPPTQQwC0b9+etWvXcunSJQ4fPsyiRYusbXfv3s2uXbvIyMh5Tq64AAAgAElEQVTAw8MDd3d3nJ0L/vz17t37qqXVsvDiiy/y008/8fzzz3Pu3DmSk5OZP38+H330ETNm5O4N//LLLwQHB9vdQaIA7sBIYpvC1ednzwANrm5SMOXhJKGxB85u0Kw7eFr+b+tf71h9rlFeeOEFmzNRdevWZcOGDbz55pvUrVuX2bNns2HDhkK99Ipj5cqV7Nq1izp16jBjxgyGDBliLatduzbvvfceo0aNolGjRnh4eNicNZozZ47VM2706NFXOQnk56677uLuu++mZcuWBAUF4e7uftUSYFHcfffdzJo1q+SDLAZnZ2fCw8OJjIykSZMm+Pn5MWrUKC5cuADA5MmTGTBgAL169cLb25uRI0dy+fJlwPB4Gzp0KD4+PnzyySc2ct3c3Fi/fj2bNm3Cz8+PsWPH8vHHH9OqVSsAJk2ahJubG/7+/gwdOtTqFQhw8eJFRo8ezXXXXUdQUBB169ZlypQpBerft29fDhw4QFxcXLm8Hy1atGDHjh3s2bOH4OBgfHx8mDZtGuvWraNnz57WeitWrOCxxx4rlz7LSFFOc37AZdOSRKTAy9JJltmrMDmV+apRo4ZUSVLPilyILbT4hx9+KL++Lp0TObZdZOf7IjvmlZ/cPOzfv99UvYsXL1ZI/46kuDG9+OKLMmjQIDtpU35Ux/8rEfPj+uCDD2TChAkVokNMTIw0atRIFi5caH125swZadWqlVy+fLnQdoV9zoBUKcfvVWAzsNby2tliS26y3K8G1puVVdQe1EvkTs8URtbcmkA4udO0PhZruKggAZoKwl5edeePw9t51rNr+sJt43VMPo2mGMaMGVNhshs3bsymTZv48ssvSUlJwdPTk/r16xMVZddkt0XxMrBZKfUtsBLDjvRQSk0A/g0UG9s1h0INlIhMz3mtlJqKkffpLhG5lOe5B0YCw8p1LF1TPtS+Dlw9IMOypHX5PFyMg9rlc0Jeo9GUjhtuuIEbbrjB0WoUiIhsVUrdi5FHcLHl8WsYYfHuFZFdZmWZ9eJ7FHgir3GyKJKqlJoDzMfIrKupTjg5gf/1cHJ37rMzf2kDZUemT5/uaBU0mhIjIhuBjUqp5kB94KyIlDhemlkD5UeuH3t+3IDq5e+syaXBDXDpHDTqYFzaOUOj0RSBUqqdiOwFEJHDwOFimhSKWQP1KzBDKfWziFhPtSmlGgHTgd2FNdRUcXrPsYs7u4jYnNHRaDTlh9j3bFukUupP4GNgpYhcfcrYJGYN1Hjge+CIUmonhpOEP3ALcAkYWFoFNCbYMMlId1G/tTGDCekN3gWHWSl37GCc3N3dOXv2LHXr1tVGSqMpZ0SEs2fP4u7ubq8uHwYGAa8CrymlvscwVuvybxMVhykDJSJ/WNYSJ2MYpRuAU8Ac4C0ROVtUe00Zif4JEg5A9I/Gff3W9jNQdqBx48acPHnyqjh3+UlLS7Pnh8wuVMcxgR5XZcPd3b1E+bnKgoisAdYopfwwjNVgYBmQopRaBywTkc1mZJlOt2ExQs+XQl9NWci8Aon5Qr5Us30gV1dXmjQpPktwREREuWVarSxUxzGBHpcGRCQRw4FuvlKqBUZsvlEYsytTtkdHkqjsJB4Eycq9rx0INR0bDFKj0WjMopSqBXSyXPUpwbGkQq2YZd1wrIgcsLwuChGRf5jtVFMC6reBCXvgzH44sw+cyyPHZBnIuGychapb8mjSGo3m2kAZm8k9MGZN9wKewM/AE8Aas3KK+rbLu1vtxNVB/wqrqylPnJzAN9i4WvV2jA6Xz8PGKXD6Tzh7CDzqwRRz0Zo1Gs21hVLqDQzHuYbAEYzo5stEpOAovkVQVCSJbnleh5VczYJRSi0H/gF4AKeB2SKy0FL2D+C/GPlDdgHDRCTaUlYD+B9wP4bn4GwRmZtHbqFtNWXEzQuiwiHLEjk75YyRi8qzvmP10mg0lZERwCcYRumnsggym1G3PA/ivgoEi4g30A+YqZTqYPH4WAtMA+pgnL3KOxWcDrTASHjVDfg/pdQ/LfoV11ZTFpxdDM/BvJz+0zG6aDSayk5DEXm8rMYJzDtJnFJKfaGU6q+UKiyihClEZJ+I5CSxEcvVDLgP2Ccin4pIGoZBClVKtbLUHQK8LCLnRSQK+BAYZikrrq2mrDTIE/fLNxgySnScQaPRXCOISDoYESWUUuOUUi8qpRpYnjVXSpnOOmp2x30ahmvgp0CSUipn+vZjCXUHQCn1HoZxqQn8AXyFEctvT04dS5y/I0AbpdQZICBvueX1vZbXbQprCxwojY6afHR+DNoPAv824O7taG00Gk0lxbIdsxxj4qAwJiHhWLZ0gIPAs2ZkmT2o+zrwulKqPYZXxkPAaKXUcYwDWCtEpOD8zAXLG6uUehK4FQgDrmB4eeQ/qXkB8LKU5dznL6OYtjYopcYAYwBcXFyIiIgwq7bd8b4QhXtaAqkeQVyq1QhxMvd7IiUlpeLGdez3ipFrggodl4OojmMCPa5rnFfI9eD7DiPyUA6bgLGUp4HKQUQiMeIsPW1RYDDwFDC1FLKygB1KqcHA40AKkP+nuTeQbCnLuU/LV0YxbfP3uwBYAODu7i5hYWElUdu+fPkZRC0zXju7wT1z4aZHim0WERFBpR5XKamO46qOYwI9rmuch4GpIrJSKZU/VtoxINisoFId1BWRbCAVI1lhBmVzM3fB2IPaB1iz41lyTTXD2Fs6jxFaKU/2PEItbSiqbRn0cjxn8qiflV6twhtpNJpqS12gsOyJTkANs4JKZKCUUi2UUi9Z9ne2AXcDH2DE5jPTvr5S6iGllKdSylkpdReGtf0eWAe0tThiuAMvAHtFJGcP6WNgqlLK1+L8MBpYYikrrm3VIzPdyL2UF/+2jtFFo9FozHMMY/umIDoBpvNCmVqWU0o9gbGeeDPGzOlzDAPxg5QsjrtgLOe9j2Eco4GJIvKlpZ/+wLsYG2y7MPa6cngR4xxUNMbM7XUR+RpARBKKaVv1yLoCd0yGuN8h9ndjic+rgaO10mg0muL4GHjO4qOw1vJMlFLdgEkYXtamMLtvNA/YjGGk1onIZdOq5kFEEoCuRZRvBgp0Dbe4po+wXCVqWyWp4QXd/mO8FjGiOTiay0kQ+5sRvDbxb6jbHG59wtFaaTSaysVsjC2XZcBCy7MdgDuwWkTmmxVk1kAFisjpEqmoKT+Uglp1HK0FxOyClQNy7xuGagOl0WhssDjAPaSU+i9wF5aU78DXIrK1JLLMuplr46SBRh1t70//BVdSoIZnwfU1Gs01i4hsB7aXRYaOZq4xj0ddY1nv7GHjXrKMPbImdzpWL41GUy3R0cw1JSOkN5w/BoGdjathaPFtNBpNtUYplU3RNiIvIiKmVu/sHs1cY4LsLHDKf76tktDrZUdroNFoKh8vYd5AmcbB2e80BfLjPPjpXagXAn4toW1/aFqo86NGo9FUGJazpdswDti6AJ+JyIt564jI9Irou6g9qBJtLIjItrKrowEg4SBcPgcnfjYu/7baQGk0GkdxBeguIilKKVeMEHWbRGRnRXdc1Awqgtwpm6L46VslXZOqXFy5AtnZxVSKOwYZ7rn3Xq2Mo8mm+3DicqlOquXi7m54t2s0mmsbSzCGnHiorpar3JfzCqIoA9Utz2sfYD7wF7AaIzqtP0aYojYYeeY1JkhMhNTUIiqIQHQWpDXKfZbWCmLM9xEfX4OYEtQviOBgcCtT5i+NRlNFcFFK/ZrnfoElqLYVS9DX34DmwH9FZJddFCusIO+BKqXUEuBbERmVr9rHSqlFGHk/witEw2sNpWDkt3AxDs4fh6Ro8KjnaK0KJuMynDsGru5Qp6mjtdFoNKUjU0Q6FlXBcvi2vVLKB1inlGorIn8V1aY8MOsk8S9gQCFlazBmVRoTmIpc6OQMPoHGRZeKVqnk/PU5fPciXDgJCNw8Cu5509FaaTSaCkZEkpRSEcA/MVbUKhSz0cydMKZ2BdECvf9U7SjSkDq7wYUYrMvQOQd3NRpNtUMpVc8yc0IpVRMjF6BdMkWYnUFtBF5VSiUCa0Uky7Im2R+YCWyoKAU1lZC6LWzv4wtL/aLRaKoBDYGllu98J+ATESnyO18p5Q30Bq7DCBKbFxERUwcqzRqo8UAgxnJeplLqPOBrab/DUq4xQYmSk1RW6jYzZlFZ6eDbBPzbQOYVcDGdh0yj0VQRRGQvcKPZ+kqp2zF8EnwKEwmUn4ESkUSgi1KqJ3ALhkU9BfxsSXOhuZZwdoXR34NvsJEWRKPRaHKZBxzHyBn4p4ikl1ZQiSJJiMh3wHel7UxTDGkX4I+VUKeJcfkEGR5yDqDYmV4DU0mUNRrNtUdrYICI/FZWQSUOdaSUqs/Va4qIyImyKnMtUOQXf+Jh+PXD3Pu6ITBwVYXrpNFoNOXICYywSGXGbMp3b+Bt4MEiOtaefGXl/DHbe98gx+ih0Wg0pWcG8KxSaouIXCyLILMzqP9ieOwtAv7EiM2kKW/OHbe9d+Dh12rhzKHRaBxBH4xIQ8eUUj8D5/KVi4gMNSPIrIG6C3haRP5rXkdNQRT5xR98Gzg5GREkzh2Fujo6g6b0ZGRkcPLkSdLS0uzed+3atYmKqn7HD6rquNzd3WncuDGurq726O4ODE+9ixih8PJj+uevWQOlgL/NCtWUkqDbjKuqcGInxOyCxIPG/lm353TU9UrEyZMn8fLyIjg4GGXnyL/Jycl4eVU/D8+qOC4R4ezZs5w8eZImTZrYo79y68SsgVoN9AW0S/k1gqklvj+Wwx/Lcu/P/KUNVCUiLS3NIcZJU7lQSlG3bl0SEhIcrUqJMWugvgXmKaW8gK+4ek0REfm+PBXTVAH8WtreJx50jB6aQtHGSQP2/TtQSu0EvsdI2bRdREqd/MesgfrS8m8TYFie50JurijtxWeCauV8cJWBOuQYPTQaTWXiCDAUeBZIV0rtBn7AMFo/i4hpJzuzwWK7FXJ1z/OvphphypA2aAsdhsNdr8Kgz+HfH1S4XprqQXBwMJs3l23HICwsjIULF5qqO336dAYPHlym/sqCPfqPiIigcePGFdqHGURkkIg0Aq4HJmNEHXoM2AKcV0ptMSvLbKijrcXX0pSJi3Hg2cDw4qsq1G4Mfec5WgtNCThYwauwLVsWX0dzbSAiB4ADSqmPgTBgAvAPy2tTlDiShKZsFDgzyUyHj+81gq36hUC9EOgy2cgLpdFoNFUIpZQ7hqt5zupaByAVI7D4FIylPlMU+nNdKfVxCa6lZRvSNc7ZwyCZkJEKp36HoxHaOGmuGa5cucLEiRMJCAggICCAiRMncuVK7jbFl19+Sfv27fH29qZZs2Z8/fXXV8k4deoU7dq1Y86cOQAcO3aMrl274uXlRc+ePUlMTLSpv379etq0aYOPjw9hYWHWs00fffQRffv2tdZr3rw5Awbk5mpt3bo1kZGRgOF48P7779OiRQt8fX154oknEJObzDt37uS2227Dx8eH0NBQIiIiAFi9ejUdO9omt33rrbfo16+f9b2aMmUK1113Hf7+/jz22GNcvlxqH4SK4jywHuiI4b/QBagrIn1FZK6IRJoVVNR60p0WwWYvTWlJyJf7q34rx+iRh2rlzKGp1Lzyyivs3LmTyMhI9uzZwy+//MLMmTMB+OWXXxgyZAhvvPEGSUlJbNu2jeDgYJv2x48fp2vXrowbN44pU6YAMHDgQDp06EBiYiLTpk1j6dLc39AHDx7k4YcfZt68eSQkJNC7d2/69u1Leno6Xbt2Zfv27WRnZ3Pq1CkyMjL48ccfATh69Cipqam0a9fOKmvDhg3s3r2bPXv28Mknn/DNN98UO97Y2Fjuuecepk6dyrlz55gzZw79+/cnISGBfv368ffff3PoUK7D0cqVKxk4cCAAzzzzDAcPHiQyMpLDhw8TGxvLSy+9VLo3vuJIwQiJ5w/Ut1wepRFUqIESkWARaWLy0iEPTFLwEl8auHnn3tdzvIHSaOzFihUreOGFF6hfvz716tXjxRdfZNky43zdokWLGDFiBD179sTJyYlGjRrRqlXu52P//v2EhYUxY8YMxowZA8CJEyfYvXs3L7/8MjVq1ODOO++0mRWtWbOGe+65h549e+Lq6sqUKVO4fPkyP/30E02bNsXLy4vIyEi2bt3KXXfdRaNGjThw4ABbt27l1ltvxSnPPvGzzz6Lj48P1113Hd26dbPOropi+fLl9O7dm969e+Pk5ETPnj3p2LEjX331FbVq1eJf//oXq1YZQaIPHTrEgQMH6NevHyLChx9+yFtvvUWdOnXw8vLiueeeY/Xq1eXy/1BeiEg9jPxRHwHNgI+Bc0qpXUqp15RSvczK0ntQlYH2AyH0YUg+ZcymfCv+tHe5cyXFCHZ7+Tw0udPR2miqEHFxcQQF5QZGDgoKIi4uDoCYmBh69+5daNsVK1bQvHlz7r//fht5vr6+eHjk/mgPCgoiJiamwP6cnJwIDAwkNjYWgK5duxIREcHhw4fp2rUrPj4+bN26lZ9//pk77rjDpv8GDRpYX9eqVYuUlJRixxsdHc2nn35KeHi49VlGRgbdunUDjNnfU089xQsvvMDKlSu59957qVWrFvHx8Vy6dIkOHTpY24kIWVlZxfZpbyxJDvcCbyulnDD2pKYB/wc8jcljSdpAVRaUAu8A46oEmF7iuxgHC8Ig5Yxx7+kPU/SBXY15AgICiI6Opk0bI2zbiRMnCAgwPgeBgYEcOXKk0LbTp0/n66+/ZuDAgaxevRpnZ2caNmzI+fPnSU1NtRqpEydOWA+rBgQE8Oeff1pliAgxMTE0atQIMAxUeHg4x44d47nnnsPHx4cVK1bw888/M3z48DKPNzAwkEceeYQPP/ywwPJevXqRmJhIZGQkq1at4q233gLAz8+PmjVrsm/fPquulRWllCtGctuco0idMZb94jEO8JqiCvk0ayolHvXg0tnc+5QzRuJFjcYkDz/8MDNnziQhIYHExEReeukl65mhkSNH8tFHH7Flyxays7OJjY3lwIHcPVtXV1c+/fRTUlNTeeSRR8jOziYoKIiOHTvy4osvkp6ezo4dO2xmKwMGDGDjxo1s2bKFjIwM3nzzTWrUqMFttxlxMLt27coPP/zA5cuXady4MV26dOHrr7/m7NmzhIaGlnm8gwcPJjw8nG+++YasrCzS0tKIiIjg5MmTALi4uHD//ffz9NNPc+7cOXr27AkYM73Ro0czadIk4uPjAWM/y8y+lz1RSn2H4SixFXgSOIsxa7pBRBqIyENmZWkDZWeqnfOBsyvUaWb7LEHHFdaYZ+rUqXTs2JF27dpxww03cNNNNzF16lQAOnXqxEcffcSkSZOoXbs2Xbt2JTo62qa9m5sba9euJT4+nhEjRpCdnc3KlSvZtWsXderUYcaMGQwZMsRaPyQkhOXLl/Pkk0/i5+dHeHg44eHhuLm5AdCyZUs8PT3p0sXw/fL29qZp06bcfvvtODuX3bs2MDCQL7/8klmzZlGvXj0CAwN54403yM7OttYZOHAgmzdv5oEHHsDFJXeh6/XXX6d58+bccssteHt706NHD/7+u9J93tIwlvNuAuqJSH8ReVdE9pVUkDLrFlkdcXd3F3unIjhyBCp6yXjXrl107ty5TDIaNABv7+LrAfDJUNj/pZGmvv710OUpaHRTmfoviIiICMLCwspdriOpyDFFRUXRunXrCpFdHFUx6rcZqvK4Cvt7UEpdEpFSedlVNIXuQSml1gL/JyKHlVJDgI0icraw+pprmLtnw73/A7dajtZEo9FUEpRSfYCuQB2MZb4IEfmqJDKKWuL7l0Uw5LoLlgmlVA2l1CKlVLRSKlkp9YdS6u485f9QSh1QSl1SSv2glArK13axUuqiUuq0UmpyPtmFtq1M2ExYs7Nh/QTY8Rb8/TWcO1Y11wC9/LVx0mg0ACilvJRSWzEO604AegMTgXClVIRSytOsrKIM1Bng1pw+KUEWxCJwAWIwrGptjHXKT5RSwUopP2Ct5Vkd4FdgTZ6204EWQBCGV8j/KaX+CWCibeXkwkmI3m7kVPr2Ofh0mKM10mg0mrIyC2P/6RGgpog0BGoCQyzPZ5kVVJSB+gR4SymVhWGcdiqlsgq5Ms10JiKpIjJdRI6LSLaIbACOYcRqug/YJyKfikgahkEKVUrlnMobArwsIudFJAr4kNzUH8W1rZwk5tvcrBdiuJtXAqriRE6j0VQK+gNTRWSFiGQBiEiWiKzAmET0NyuoqHNQk4AfMUKmvwgsAWJLq3FBKKX8gZbAPuBxYE9OmYikKqWOAG2UUmeAgLzlltf3Wl63KawtkC+OkGOx+eLPH+JIR5DQlDMiopMWakzHCCwn6gL7Cynbbyk3RaEGSowRfQqglBoGvC0iewqrX1IsB7lWAEtF5IBlXTJ/TuILgBfgmec+fxmW8sLa5u93DDAGjPMGOUEa7cXJkzWtRsotszmezUbikRpDrUsxxKd4cH7XrjL3cenSJXaVUY6vbzqenqV3N1TZGYiTa5l0yE9KSord/78qmoock6enJydPnqR27dp2N1JZWVkkJyfbtU97UBXHJSJcuHCB1NRUe31+jgF9gO8KKOttKTeF2XxQ5Rp7xxL6YhmQDoyzPE4B8js2ewPJlrKc+7R8ZcW1tUFEFgALwHAzt7fb8qFDhS+f+ZRTH+XhZl6/PviURKG0i7D/CzixC2J2Qg1vGPNDmXTIj3YzLxkZGRmcPHnSGsLHnqSlpeHu7m73fiuaqjoud3d3QkNDcXUt3x+NhfAB8KZl0rECI2FhA+AhYBRGEkNTmA51pJRqCDxFPrdBYK6InC6BHAUswoh021tEMixF+zDSBOfU88DwHNwnIueVUqeAUHKtcqilTZFtzeplL6rt3k5mGqx/MvdeOUN6KrhVyuMV1wSurq40aeKYuI4RERHceOONDum7Iqmu4ypPROQtpVQ9jG2iYZbHCrgCvCYib5uVZSqShFKqJcYez3iM2covGAmoJgCRSqkWprWH/wGtgb4ikjeRyTqgrVKqvyXh1QvAXktWRjAi4k5VSvlanB9GY+yLmWmrqWg860OdPEHtJQvOVLrfBxqNxg6IyHMYfgN9MBzc7gECROT5ksgxG+rodYw9nZYi0k1EHhaRbhgODhcs5cViOZv0KNAeOK2USrFcg0QkAcO74xWMOE6dMaaEObwIHAGiMWI8vSEiXwOYaKspIaWa6TW4wfY+8VDB9TQaTbXH4nG9yeLNt0lEzpdUhtklvm7AYyJyPJ8C0Uqp6cB7ZoSISDTGVK+w8s1Aga5sInIFGGG5StRWYydCB8J1t0H91uDfBjz8HK2RRqOxA0qpEuXYEZFtZuqZNVBuFOBwYCHZUq4pBptZyZGtULcp1G5cac4+lZmQfzpaA41G4xgiyA3mYCawQ7nmg4oEnlRKbRIRa8hdi8PDWEu5xiyp5+CrScZrNy/wbwv/erdSGapq68yh0Wgqgm55XvsA84G/gNUYUYn8gYcxzqY+YVaoWQP1ErABiFJKrSHXbfABjPBD95jtUAMk5vHdSE+GtPOVyjhpNBpNSRCRrTmvlVJLgG9FZFS+ah8rpRZhRP4JxwRmz0F9bYlMOxN4ntwp3G9AHxH51owcjQUdQUKj0VQRlFKBGF7UDYBsYEExruL/AgYUUrYGY1ZlCtMJC0XkaxHpiBGdIRDwEpFOIlK50jlWYqzLZh71IaADuFkCZGgDpdFoKi+ZwFMi0hojjfsTSqnri6jvBDQvpKwFJvefoAQHdXMQkUvApZK20+ShdR/jys6Gi7G5hqoSUaY9qIw0OPMXxP4O2Rlwq+klZ41GU8kQkVMY2zqISLJSKgpoROHx9jYCryqlEoG1IpKllHLGOAo0E2O7yBQlNlCacsTJCXwCHa1F+XL2CPy3s2GYwJgt3jJW77FpNJUXF6XUr3nuF1hCwl2FUioYuBEoKtjneIxVtjVAplLqPOCLYW92WMrNKWa2oqbsXBOecT5B4OSSa6BS4428V9XBEGdegX3rIDUR/K+HJl3ByfRqhUZTWcm0bN8UiSW23ufARBG5WFg9EUkEuiilemIsCTbEmIH9bDmvahptoDQFUmpj6uwCAe3hxM+5z+J+r9oGKu0ipKfA6oEQ90fu84ah0LY/uNaCNveBh+ksAhpNlcKSfeJzYIWIrDXTRkS+o+CI5qbRBkpT/jTqABdiodFNxmv/to7WqHRkpMGGibBnVcHlp/YYF8DW2TBsI9RraT/9NBo7kCfAd5SIzC1hWydgM/CoiJQ49lmxBkop5QbsBJ7V7uRlQ7KyIf0KuNV0tCoVS48ZcNcrjtaidIjA+WNwbBv88qHh7GGG1Hj4uB+M/h68AypWR43GvtyOkb79T6VUTlCG50TkKxNtFRBGAbn5zFCsgRKRdKVUEwxXQ01ZOPMXfDgaGrSD+m2MA7pJMeDVAPyag3djaHwz1PJ1tKZl2y9zrqIT8+QzsG4MHI0oZftTML8D9JoJHUdoxxBNtUBEdlBEDNWKxOw3yXdAL+D7CtSl+nNsu+E8EPebceVwGsiZ/LrUgq5Pw/X/coSG1x4p8RD7G6QmwJaXjZlQYTi7QZcpcMdEOPiNcUUut62TcQk2TjaSN7Z7CILvAN+gih2DRlNNMWug5gPLlVIuwBcYHhk2v7FF5Gg561btkJhfiq+UeQm2zIBTe40vwhqlmhlriiM7y3iff3rXyF1VFCH3GD8a6jQF99rGs+v7GVe/dwzniYNf27Y5ts24nGvAwNXQrHvFjEOjqdxkAzOAuNI0NmugcuIsTcbIklgQ2t+2OM6VwIbvXwfxUXDfB9pIVQSbnoHdHxZe7t0IOg43jJN/EYfmnZyh/yL4+lmIXAG5sZQNsq7A+gkw7hdwreZ7jxoNoJRqLCInAUREMAxU3vKueWP3FYVZAzW8ZCpqCmR0BLw2CM4etH3erIfhxhyz0/Z54gFY9xj0mQee9eymJj92liQAACAASURBVJTzma2My8YB3gaVxJvv6NaijVNIb7j/I3B1NyevhqcRjf72ifDpMDjzp235hROw6324o7DfdhpNteIbpdTtIpKUv8CSN2oDJp0mzAaLXVoy/TQFIb8tudo4/ftDaNzBeH3yN2PZ6eLJ3PKEKNgwAQYsq1qHQrMy4aspxhmoM/sABc/FgksNx+p1+k/DiOTH/wZj9nPjIOj8uBHlo6T4NYcxP8BvS4yx52X3YrhtQunkajRVixRgo1LqHyKSlvNQKXUHRhik9WYFlejTopRyUkq1VUp1VUp5lKTtNc+x7cYyUF6adss1TmC8fngVNAi1rZdwAI5UMf8UZxc4vNk4J5SdaTiHnDbpsl1RnI+GZf+Gy+fyPFQw4lt4fAeM/cmIG1gWI+LsCp1Gw/8dA5c8M7ALJ+DET6WXq9FUHe4B6gKfWs5BoZS6DfgKw0ANNivI9CdRKfUEhr/ZXgxvvhDL8y+UUqZjK12z/DgPG78S5xrGL+r8uHnAv+Ybv+jzsnuxsbFvJ8plia/RTbb3sb8WXK+iObHTmDW9e7PhrZeXLpPhus7l32etOtAqX5q0yJXl349GU8mwhDr6J3ATsEgpdSuwCfgGGGTZlzKFKQOllBoNvI3hwTcAW5/47RhRajWFkXkFjv9o+6z7C+B7XcH13Tyhx3Rs3uazf8OviytIwQqi8c25r+u1sp1R2IsTO2FJHyOGXtYV27LOj0H3aRXXd7uHbO8jV1z9d6DRVENE5DhwN3AvhpPdZuAhkeJcZm0x6yQxGXhTRJ6xhE3PywHg6ZJ0es1xcjdkXkbEzbj38IeQfxbdpk4TaN4DDucJZbXrA+NLP6B9xelanrTua0RVCO4CnvXt23d2Fvz4NvwwKzdwbV7a/BvuerViD9M26w6+wXD+eO6zJb3hgaXQ5t6K61ejsTNKqRGFFK3HMFTfAkOV5fMmIqZ+bZs1UE0wpmcFkYqRg15TGNvzha8KvNncF2OXp41DpNY9k2z45nkY9Cm41Sp3Nf+/vfMOj6s49//n3aJVb5Ysucmy3KtcMG4YG2NjwMEU0zElATsUX3LJTW7IDwgEcgMJIYWQQEzopoRiDAaMaTbFGIN7wd1Ytlzkpt5XO78/zsq7Z3clrWyV3dV8nuc80pmZc3ZGoz3vmZl3vq83LTLFl5JtHG1NxQl4/QbY+2Xg/In/A5N/fdoOCxUlNVSW1pDcORarPcC9rDZjJOzrlPH2T6HrCL2BVxNJ/LuJ/Ce9fldAixqoY0B2A3n9gQNB3qdj4us+HpVkOj1xsIxda49wZG8psUlRnHFhNomdYiA+DaY+CIvneQqXHTI2hQ65rA0qHoaUHzc08QJp6GWNg6tePmXVceVSHNxZxMFdRWz+/AAVJTUA2KIsjLogm1Hn90R8XzwGXQIDZ8JWL8clZ5VhQC8LGHJHowlHerXGTYM1UIuB34jIciDPnaZEJA1j4+6iVqhbZFCcDzXlACj3mpKr/wWsfm8PB7YXUl5cYypedKSC9/6+gYETujD4rK5EZY831jI2vuYptHWxNlCBqC6FV6/yN04xKTDpbhjz02ZN6dXW1LFn3VEO7iyirLCaooJySo5V+ZVz1rhY9c4eqiucjL+st9lIiRh7qt79L9jg5SRxaD08dRbds2eDmqR1+zRhjVIqr+lSzSdYA3UvMAXYjBFJUQGPAwOAI8CDrVG5iGD1c+Zzi529h9PY8e3uBi9xuRRbvjzIkb2lTP3JIKy515gN1OENhipFak4rVToMgysW58OCy419Y95kjYcrnoeEjGbdbsd3h1nx5i4qfF4gGmP9x/vY8e1hxl/am/5ju3gyrDa4+B+GSrp3nKy6GvrsfhY+ijUEZrWR0kQAIjIWyFJKvR4g7wpgn1KqsYi8JwlqEl4pdRw4A3gYsAO7MYzbE8A4pVRxkHXveFT5bKZO6UX+tsKgLj26v5SNn7mj0XYdZc5c+Y8WqmAEUH4MXrzE3zhlT4TZbzXLOCmlWL1kLx8/832zjFM9FcU1fPL8Vr57/wdM3rQWC1z9CvS7wP+ilU/Ahtf80zWa8ORhYHADeQPd+UER9CqxUqpUKfWQUuospVQ/pdQ4pdRvGwv92+FRym+DrRo9l8N7/f9kGb2SSEzz12rb8tUBDu0qgqFXmDP2LIMDa1u0uq3Gka2Gi/y7d8JTEyGvBTesulzw1s1w3CcWWteRcPXLQTmTFB+tYP+2E1SV17LirV2seqcRzUSBXrlpDBzfhSGTunHeLYOJivGfiPh28Q/887ZlvPmH1Rw/UGYkxqYaG7EvfxYc5nVI3v8fYyOxRhP+5GLEEAzEt8CwYG/UrMA9IpIIDAG6AfnAFm2gGiH/O5NArLI4KLQNpqbC/AD80bxckjsbD9LCgnKWzt+Ms9YtOqpgxcJdzLhtMjEZQ806b6uf9d8MG4p89VfzFOWBNdBzfMvc+9v5/vGb+s+AWf9u0jjVOV2seGsXm5blN1qunsS0aCZc3pec4WZdxITUaN7/x0aqyv3d2Qt+KOHNR1Yz445hdB+QakzjDZkFnfrAM9PBWWkUrC2Hpf/PMKoaTXgTTcODHysQtApRc5QkfgPsx9iY+x9gBbBfRO4N9h4dDl/lgL7TOHTAHPexx8DUk8YJICUjjnOuG2Dao1tVVsvXb+9GjfcR7Nj3NRzd3tK1Blp4Dcp339a+hl6umklVMSz/vTktazxc8VyTxkkpxacvbG3UOEVFW7n0FyO546kp/PTvk5j94Dg/4wSQmZPE7N+NY9ylvbFY/NeRnLUuPnx6MyXHKj2JXXL9ow5vew926KDVmrBnKzCzgbyZQNAPrWCVJH4LPIBhmKYBQ4GpwOvAb0XkgWA/sMPgcpldiwE19CoK9pgHnJk5PlM9QEZOEkMndTelHdpVxMbtnVEZPjp968Ng7SJrnPn8wJqWsYDfPGkYqXociTDr6UYFaetqXWz9+hAv3vM1O78raLBcTIKdi+8aQdc+xhY/m92KBDA+Jz86xsbI6T256r4zSc/yF2quLney8NE1HD9Y5kkc9WOKE/ubCy75JdT6ewpqNGHEU8AcEXlURPqJSKyI9BWRR4GbgX8Ge6NgR1BzMJQk5iqlPlNKbXH/nAP8BZjb7CZEOke+h4rjnnNHEnU9J3Mkz9dAJQa8fOikbnTuaX7Qbfo8n1U1N5uf7Ts/hEo/VfvQInOYoaow8Rdw43vwsw2n77FWmGcoRXgz/k5I6h64POCqc/Hu4+v57MWtlJ2oDljGEWtj+LQsrvnNGDr3DNw3jZHaJY4rfn0G1z04ltwpPUx55cU1LHlqE84at9qLxcLOvreCeH0NC/f6t0ujCSOUUk8Df8bYgrQVKMVQHLoL+ItSKugNgMGuQSXRsJLEh8BtwX5gh+Fbn3hDPcdzZH+FZ20JiE2wB3SMALBYLYyf1ZcPntpITYVnWnDXDivdMybR3eWO91VXDVvfg5FBCwQHRYtO8VkscP3bLXe/8mPw+vVGePV6YlKMfU6NsPqDvRzc6W/Mo6KtzPzvEVisQmqXOKy201OYEBGSO8cy4fI+lJdUs2u1J4x88ZFKvnt/L+Mu7Q1AWUIOjL7FWEur58vHjGi9nQeeVj00mvZCKfULEXkSY8YtFUPs4ZPmRl4P9pu4ChjdQN5od77Gmx1LzOdx6eRvNz8cM3OS/ZUHvIhPdnDW5X380nfJheaE7e+fcjXDDlcdvHyFEcbDm6m/hWj/EU9ZYTWblufz0b838937e/3yHbE2ZtwxjIzsRNJ7JJy2cfJGLMK0nwym72izm/uGT/dTXuQ1gjvnHojzWtuqq4YXZhpBHjWaMEUptVsp9ZRS6vdKqfnNNU7QiIFyx36yuON53An8RER+KSLZIhLj/vm/wE+AeQ3dp0NSVQJl5vUN1Wsyu9eZQz1kNDC9503XPilMvLKfKS3/cDxHant7Eo5t7zgPs/UvG0EQvcmeCCOuNyXVOV1s+fIAr/z2G754bQc7vUYxJy8blsbld59B174prVZdi0WYfG1/4pI962J1Thffvf+Dp1BMMkx7yHxh+RF49Rrjf0mjCTPc607zROQNEflURF4XkdtFpFkioo29LjqBWvexAegNPIKxSbfM/fNhd/rGU2lExPLln3wShBPJUykq8ExJWSxC937BPRizBps9/QCWlc6j2On1Zr65BafQCE0lCVttGXz2O3Ni99Fw1QKT8GvJ8UpeeeAblr+8ndoqf3V/q93CzJ8NZ8btw/z+rq1BVIyNMy7MNqVt+fIgi/68lqpi9x8692o442bzhce2+we51GhCHBHJBNZiqA2dAcRizLQ9AawRkaB3zje2BvUgpgh7mqA56hPWvfNAdq49ZkrK7J2EI84e1O1EhIHjM1m5yDNCrq2L4tuyq5ma9DfD32Djf2DAhZAx6HRrH3o4q2HNC5y1wieqi9Vh6NzFeMT0a6vrWPLUpoCaeWB4502fM4RuQb4ctBQDx3dhzYd7Tc4ZB3YUIbuhILeEjOxEmPGYoSe4yUshZv3LMPRyw8lEowkP/gikABOVUicDoLmj6r4F/AG4KZgbNWiglFIPnFYVG0BE5mFUbijwqlLqJq+8c4F/AFkY61o31YsQiogDQ7L9cqAC+KNS6s/BXNumlB83Qp17oabcz86Xj+C9uann0OYpaueM6Ezh4Qq2fXP4ZFpBbT/yqkeRHb0GqINPH4KrXjyd2rc+rjo4vAnyVkDJQf+9QL5Ul8JLl0H+t/55E+40ZKDc1FQ6Wfz39RzbX+ZX1BFrY8zMHPqPzSQquln701sEq83C1BsHsfjvG6hzehxlVB0s/vt6pt8yhB4DU2Hm48bfx1u2afHP4LaV4Ihv83prNKfABcCvvI0TgFLqa/e+2UeCvVHLrQgHz0Hgd/jEA3Eroy8E7sPw+liNse+qngeAvkBP4Bzgf0Xk/CCvbTu2LDQHyOvUl6PR4yg56tmkabEJPQakNuu2IsKoC7Lp0sccemtTxYW4lLsbj2+Hb/91ylVvdapL4Y+9YP4kQzXhm3+a9zEF4otHAxun9AFwtnlE9fmr2znss8/MFmVh5p3DueXPZzN0cvd2MU71dOufwow7hhGbGGVKry538u7f1rNswTbqLA64+AlMO7WL9sE7d0BdgMCLGk3oEY/xnA9Evjs/KJqjJDFQRO4Tkfki8qLP8UKw91FKLVRKLQKO+2RdhiGd9IZSqgrDIOWKyAB3/g3AQ0qpQqXUVuBpPMPEpq5tOza8aj7Pvdpvgb5bn+RTelCKCKMvzDZtISquy2RftZdSw+pnSTt2+koNrbIG5UiA+EyvD3HBvkYcQI/vhpUB9vSl9YMrXzRtyN24LJ8d35odU9KzEvjxH8+ix6DmvQy0Jj0GpnLjw+MZNKGLX973Xx3ksxe3orqOgrG3+2QuMvT6NJrQZztwfQN5szH2RAVFUE9JEbkBY8SjMMJr+Mo8t8TjbDCGM4ZxQ6XKRWQ3MFhECoCu3vnu3y9p6lqa8cc4bQrzDJUEL+oGXcn2xXnuehlpPYemnfJHJKbF0Cs3nT3rPR6Be2onkh3t8WzLynsTqm8OzSmh7LOMxf968r6CfucFLvvRvabRaK0tDvvcT43Rk9tKV5XXsvTpzX4K8Z26xXHxXSPadcTUEBarhcmzB5DQKYZV75o9b3esKqBzz0Ryp9xjbFXw0nJk7Qsw/FrIGtvGNdZomsWfgBfdzhCvAIeATOBqDAWihoyXH8F+e+8D3gFuVkq1lmxBPHDUJ60YSMAzJCwOkNfUtSZEZC5u5Qubzcby5ctPq9LedDn4Ed7CNUVJg/j8/d1Ulnrst8UGh0p3UbDq1N3CnQnm94FD1X2pSEgi1v3nsdeVs/+DP3Gwa4DQDkHicLjYvTuw2sLpkF6RQj9bAkXJgylKHkJhTT8qAvRByon15G7/wJS2qcdNlHxfABgjJaUUez9TVPj0vFggeUgFK1d91eL1b1FiIWN8FSc2RFNb7kle8dZODpUKyX3vJnfDfcRUeUaGpa/fzppRj4V87KiysrIW/W6FCpHarpZEKbXA7U7+IOZQ8AXArUqpVwJf6U+wBirTfePW1NQpA3w3BiViyGSUeZ1X+eQ1da0Jt8zGfIDo6Gg1efLk0623h6fNcRuTR16Ga3US4Pmz9R2VyejxpxcdWSnFe9s2UOxe11JKOJh5G30Oe9Yeexz9nB4zfgVRgZUqmiI6GrKyTquagambAHIP6RYL/rKr9WWc8NSvzGk9xlKSNQ3v/tq35TjfHzVv2LXaLUyfM4Rew059lNqWLF++nGkXnMnr//cdNW6XeFUH+z6z0O3SCURf9Sy8MONk+YSy3UzOtkKvs9urykGxfPlyWvS7FSJEartaGqXUfBH5N9Afwy/gBLBdKeVq/Eozwa5BrcAINNWabMGIIwKAiMRh7LHaopQqxBgmeiul5rqvafTaVq6zhzqnX3ymE7ahftI6fc7ofNofJSJkDTJ7AW7el0ON3euhXF0MPiOQkMBqN+1Z8kMp+Pg+OOozM3v+w6ZRg6vOxTc+cZtsdguzfjkqbIxTPUnpsUy4oq8prba6ji9e28Gyrzqh+s0wX/D1E21YO43m1FBKuZRSW5VSK9w/XSIyVUSC3jcbrIGaB8wVkWtEpJO3yoSX2kRQiIhNRKIx4oJYRSRaRGzA28AQEZnlzv8NsFEpVf+kehG4V0RS3M4Pc4Dn3XlNXdv6rH8Z36W4zfvMQ5DM3kmkZAQdCqVRsod1Mjl6lRXXssHqs7C+4RXDrTtcKNoHL1xkePd5M3z2ybhXRQUVfPrC9zx5x3KO7jMPkC/9xciASuLhwMDxXRg8satf+tYVh/j0+Fycysvzb+dS/1AuGk14kETD0Xb9CNaw5APrgAUYThK1PkdzYmPfC1QCd2N4dFQC9yqljgKzgP8DCoExGItq9dyPoV6RB3wOPKqU+hAgiGtbn6//bjp1xXdlx2rz5txB4/0fQKdKUnosA8aZPcG2788gv2aoJ6HwB9juowkYJG2qJOGqg1Xz4V9nw94vzXlxnWHq/QA4qxVvP7aWbSsP+92izxmdT0l9PFQQEc6+pj/Dp/nPq27f4uLDyoc82wkA3vs5FLTdBIFG0x4Euwb1NHAVsAjDK645BsmEewPwAw3kfQIEdA1XSlVj6P79pLnXtjqFe6ForynpYLfbqN7lUSF3xNrIHpbGUV/n+tMgd0oP8rcVUnbCo5rwRclPmZjwND0c7rWZVU9B3/PAFtXAXdqZulp488ewdbF/XnQSXL8Q4o1p0UNrFBUl/v96UTE2zvKZIgtHLBZhwqw+9Dszg3f/ut4UoTevpB8rXLcwMd6teu6shNdvhNtWNBr/SqMJZ4IdQV0M/FIpNUspdY9S6re+R2tWMqSpKTeUDrw3UcaksIeppmK9hqVhaUGlbAB7lJVxF+eY0lxKWFk6m8o691RX6UFj83Co8t5dgY1TzjkwdzlkGiPC3WuPULLPv1hiWjTnzx1CXFLkPKTTeyRw6S9Gkpxh1gncWHYBWyqmeRKO79RTfZqIJtgnZjnwfWtWJGxZtwBOmF3Ga0ffye51J0xpvQKECm8JMnolMeI887RQjYplXfmlnoRV86Fof7Pu2yZTfMsehnUvmdNsMTDtQZi9EFIN43ssv4zlr5ijRCemRTP38Ulc/7vxhkRQhJHaJY5L7hphUkEHWFE+h2qXl+H66s9aYULT7ohITjAHhkd40ARroJ4Drm12rTsCvg/YPtNYW3qRaSrKFmVpVTWDwWd1I/dccyTZPdVjKKhxT3tVFxl6bqEUStxZA9sCjJyuXwQTfnbS0+/Q7mJD9bvM8xC2WIULbh2GPcraVrVtF+KSHcy4fRi2KM/XtLbOzuaqH3kKFe2DTW+2Q+00HQUReVZEjojI5kaK7QJ2BnH8vaEbBCLYNag84BoR+Rgjgm6hbwGl1LN+V0U6hzYawp71iIWKcx5j/SM/mIoNndQde5QVVUGrMeisbuRtOkHREc+HfFt+FTPsv8ciLmONbNObLR5595RwuYypPd9FfrFASf7J031bjrPkX5tw1pi3Toye0Yu07iGoktEKpGclMOqCbFZ5udSvLr+CLPu3pNv3GglfPgbDrgRLZBtsTbvxPEaojMaUqH/cGh8crIF60v2zJ3BugHyFj/hrh2Dbe+bz3uey8TunKax7TILdLxZQa2C1Wjjzol589IznoV/s7MLOyon0j3WHh1/zHAy+pH0lkOqchvDpxtfM6SnZcNMHkNQNgIM7i3j/nxtx1ZnnGoed051R5/dso8qGBkPO7sbaD/OorTa2DDjrbHxQeA/Xps3Dbqk21qK+eBQm69hRmpZHKfWFiGQ3USZoPdbmEOwUX68mjpyGL41gVpmVw2v6z2LT8nxT2qjzs4mKaRs9uM49E+k9wrzWtbFyhmfNoqoQVgY3wm6VNShnDbx1s79xSs6CWz47aZyO7itl6dOb/YxT+hDhrCv7IpbQlvlpaaLj7IyZaf6KlbnS2FBxkSdh+SOQt7KNa6aJEGwistrrmNveFaonqCdnu8RVCnUK86DKrBKx62hvaqo8sYii4+wMOsuz96ktHA9yp2bxw6ajuNwe7tWuODZVnM8Z8W5Pvk1vwMCZkBH0XrmWobYK3rjJEED1Jj4TrnsL4jpRcrySZS9t8xN+BTjryr4UWnYjIa5B11oMm9KdooIKNn9x4GTa2vLLGBDzGfHWE4CCJb+EuZ/rqT5Nc3Eqpc5o70oEoj3iQUUGq58xn1sdbN1sfjAMmtgVu6NtHxaxCVGk5Jgf4tsrz+FQjdcWsXUL2rRO1FTAq1f7G6ekHvCTJZDej7LCKt756/qAxmnUBT3JndLDL70jISKMvbQ3jjjPO2WtiuHLkjmeQoc3wes3GFsfNJoIICgDJSI/iMiexo7WrmhIUVPu95AvSpjA4T3m4HsDx/vH/GkLknMgLsmzMVdh4cuSn3j2Ru38CAoac8hpwdFeaQEsuAz2LDOnp+bAj5dAag6H9xTz+u+/MwV1rGfg+C6ceVHHnEH2xRFj48wfmf8We6rHsrvKK/zGtvfg+RlQcQKNJtwJdgT1eYBjM4ZiuAVY3hqVC0mObIOHe0CFWRJiV/RVpvOufZNJ7mzeaNlWWKxiPNS9BlI1Ko7Nlee7zxS8fRtsfL115x2P7YR/TYR9Pmsj6QMM45TcgyN5JSx+fD2Vpea9PAmpDqbfMogpNwzE0sHWnBpjyKRuZPQySzp9UXorxU6v7SUH18Gi29tYr0oTqYjIq8BKoL+I5IvIzW312UEZKKXUTUqpH/scM4E+wGHgk1atZShQsA3+OhT+OcaIh+CFM64nq9abVcr7js7wu0VbPi+69Uth2DnmvVE7KydyrCbbOKkth88fgbUv+V/cEpQfh5evgDJzlFsyh8JN70NCJvu3nWDx4xtOhpmop1dmAdel3EKfKB9dPg0Wi3DO7AEmo11Rl8TCoj9SVuelcL9jCbz33+Bs+Zhemo6FUuoapVQXpZRdKdVdKfVM01e1DKe1BuWOD/Uohnp4ZPPKFcamSB+qXbG8eewhU5oI5LSSckRzGDyhG/FeSgQurHxeOtcz1Qew8gnY87nftadlTJ3V8J/ZhlitN9kT4cbFlDsTWPLUJj+9OYCRcW9xgboNa/lBw3Xa1azwMR2CTt3iGenjal/hTGBZ9d3mflvzPDxzXsD/W40mHGgJJ4kqoHuTpcKZov1QHPhLvrpqNseLzSEeuvZLITax/cVZrXaL34Os0pXE+opLPAnKCe/fBV/9pWUkc6pKjJHTvq/N6bnXwg3vUFTq4M1HVptC1tczYmIiYxNeQcT9lD22Hba+c/p1ikBGz8im3xjzKH1fSR++qriVOuXlmHNoPTw9xZj202jCjFM2UO64TsMxlMkjW/d/zXMBk4vjRrGx3BxW3Vj/CRwxtz2WBLIGdTK5ugPsrhrLtopzzPVZ9xJ88tvTq2TZUWOB/gefEVn2RLjobxTsK+ftP62lrNB/2il3Sg/GXTsKGXalJzG1N9iiT70+EYzFamHqjYPo2jfZlL6xdDr/KXzCvCZVftQQND66HY0mnAjWi88lInXeB1ANrMFYh7qrNSvZ7vhEygUgoSubuj1uigdotVu45v4xdO2T7F++Hck9twdJ6WaHjdXll/NN2WxzjKEdH8C6l0/tQ0oOwkuXwGGfYJlp/eDKF9m1oYi3H1vrFy6jc3YiF981wtiAKwJn/wLSB8Jl/4Z530F/8wuAxoNYhCk3DMARa97OWFiTybtlf6DQ6fViUnnCMFLF+Wg04UKwEgcP4hsu1pjaywOWKKWK/S+JEH74wt9FeugVOM++j20Pm/cvT7qmX7t57jVGvQzSJ89/j3J5unF31TjioqsZZn/DU3jFXyGtDyp7bIA7BaC6FD7/A3zzJCd3B9fTdSRF0xewbuERvv/qoN+lA8Zmco6vl15aX7h9pSm8u6ZhktJjuejO4Sx+fD3VFZ6/f0lVIq9V/53Rca8yMm6hocdYkm/sR5uzDKz2dqy1RhMcwSpJPNDK9QhNXC746F5zWtY4XJfMZ/3SfaYHQnScPaDnnjft6fWbkZ3IpKv7seKtXSc13QC2lJ1DfHwBWbaV2KQWcMGHv4arXwD8o7uaKD4AC2bB0a3+eTmT2dLjL3zxxx1+kkUAI87LYtwlvQPLFmnj1CwyshO56t4z+eS57zm406Nu4lIWVpVdR37NMM5LeoxYa7Gxmffbp2Hc7e1YY40mOLSSRGNseh0OGZFp65SNlaXX8dbBX/PkHctZ9a55b/LA8V2w2UNbYqb7gFSmzxlCdJzn7bnOqfi66CoWnXiQwzX9jcTqYnjjx7BnecM3Kz9mrDcFMk5dhrMh408s/89eP+NksRrTUuMv69PhNPVak4TUaH40L5fMnCS/vAM1Q3m38H5qXW6Pzo9/0/ZqIhrNKdDgCEpEmuU6rpR68PSrE0LUVsKnniZ9WjyPnVWToLzOr6jNVWeUjgAAFQtJREFUbmHYlPBwZEzuHMuIaT1YuchsYKtciSwv+SnnJj5OetReVMUJePFiI7P7mZAxyPBmPLTBUNJw+qs+kNAVNXoO68pnsvJtf6/HmAQ70+cMoVu/lFOqu9VZbgxD9QgrIHaHlUt+PoLtqw6z8u3dphhax529eOPEHzkv6c+k2fMMRfkti+DsX0LWmHastUbTMI1N8T0QxPXer8eRZaDeuhlKDGHOwzX9DOPUACPOyyI+JXy8zXKGd6awoJJtKw+Z0p3KwacldzI58V90idrmycj/1jgaY8CPqJ3xFJ+9lseuNWbjZIuyMGZmDoMmdD01ZXel4PtFjFn1M8j4Mwy9vPn36CBYbRYGTehKzyGd+OjfW0xTfoXOLBadeIhLU++hk30/7PoYdn0CUx8wgkRqw68JMRqb4rM3cYwGPsIQ1NnVutVsYw6sg23vA6CUsKL0pgaL9h6R7rfXqCFCRXlGLMIZF2Qz445h9MpNM+U5lYMvSuZQWteMjca9zqbmR0+z+F872LXmiCnLYhEuvHUYw6dmnZpxKj1sbPp94yaiaouNUa2zpunrOjhxSQ5m3DGM1K5xpvRqlcDrxx/ju7Ir3B6cCj65HxbO0X9XTcjRoIFSStUFOjBiPy0AVgGDgLnun5HDfzzR7bdWTuFw7UC/IharcNW9ozn/p0NDfu2pIVIy4pgwqy9DJ5mnJ2tUDJ8U38mR2t5N36TveRyb+DQL/7KZQ7vMzpyOWBsX3jHs9MLd15TDjqWe86I8Q/1C0yRR0TYuuWsE2cPMLyEu7Hxbdi3vFj5AsdPt2LPpDXj9eq2Ergkpgn6lFZEewP3ADRgh338B/FMpFVmvXVveNvb0ADWuaL4pu96vSN/RGQw5uxtp3RP88sKR3HN7YLEJGz7dfzKtrK4Tbxz/EwNHRDNkYAnptj2I1Wpo6dXV4vrha3YXDWTTnu4c+uP3fvdM7RrHhbcN9dt/1Ww69YZRN8F3TxvnjkTo7P/CoAlMTEIUF942lE+f38r2VYdNeQdqhvLasb8xJuFlcmPfQ3Z8CE+fC1e+AOn926nGGo2HJg2UiHQG7sEYKVVhrDX9RSkVea9aShkuuG42VFxEpcvjFWVzb8RNTIs55duHKkPO7kZRQQV5m80q7VvXVbF1XRSpXUcx7JzuFK+ppKywmoIfJlFyrAoo8btX5+xEZt6ZiyO2hfbanP1L2PQ6FRJP7M3vGHulNEEjIky5cSDpPRNYs2SvSTneiYMVpT9hf/Vwzk16nNijW2H+OfCjv8CwK/W6lKZdacyLLwn4FfBfGOtMfwP+oJTyjygXKaxbAHkrACh0dmNt2SxT9vBpWadsnEIdEWH8ZX2IirGx87sCv/wTB8tZ/nLTUjm9R6Zz7o2DWjZQY0IGzF7Imm0FTNTG6ZSwWITcKT3oPaIznzz/PQe2m7/G+2pG8uqxvzEm4VUGxXyM5e25xlTqxU9Al9x2qrWmo9PYCOoHIAnDEeJ3wCEgRUQC+ggrpcI7aGHB9/Dh3QDUuGJYWvQ/OPEogTvibAyf1sTG1TDHarMw5qIcuvdPYc2HeajayqBfoDNzEhk5vSfZw9JaJyx79zOo27W85e/bwYhPcXDxfw9n69eH+PqtXabN5lUqic9LbmVTxQVMT/4TqYc3wrMXwBXPQ7/z2q/Smg5LYwaqXlBuOhDMf2d4egqAIXL6ypVQU4ZLWfiw6H857jQLvk6Y1RfHqXiheRHKU3zedOuXQmZOEnXHDvHNot2mh5g3FqvQf2wmI6ZlkZIZF7BMm1BTDrYYsOh958EgIgya0JUeA1P5+Nktfs4tJ5w9WXj894xPeIEBMcuwvHo1nPsbGDcPrKf3HdBomkNj/20/brNatCe1VYb3UrHhIPBN2Wz21ww3Fek9sjMDxmUGujpisdosDJjYjT6jOvP5K9s5sLOIpLRoMnsnk9ollrhkB52zEomOb2dNt9oqWHA5JPeAmU+Arf3DnIQLCanRXHLXCNYu3ceaJXtx1npib1WrBJaVzGNn1UQmJz5J0if3w9oXYeL/wJDLwB6ZU92a0KJBA6WUeqEtK9IeWFy18H8ZFNT0ZU/1bPZUjaGozuxyndErkXNvHNg601YhjlKGxuD0OUPauyqBcdXB23ON2FP7MKL3XvkSRCc2eanGwGK1cMaF2QwY14UVb+1k12rzPrb8mlxePvYPukVtZmjVB/RadAfy3l0w/FqYej/EnJoqiEYTDB16vC64OF6bxcIT/4cL/5FAXLKDC28b1mIL/uEyxRc2LPs9fO8V0HDPcvh2vhGyQ9Ms4lMcnHfzYFIyYlmzJA+Xl+q9wkp+TS75NbmkWPfT07GGvt98QeftY2DaQ8aISquja1qBDj9pv7f6jIDGyWITzv/pkJCIjKtpgBGzITXHc545FMb/V/vVJ8wREc68KIfrHhzrFwixnsK6HqyvuIQ3jv+Jd/fdyqaXF1L653Nh62L9BqZpcTq8gSqt6+yfKDDpmv5k9vJXhu5IhPzzJrUX3PwxdBsF9ji4/DmwOZq+TtMoiWkxXHzXCKbeNJC0HvENlttfM5wvSn7Kgt338PW/P6LmiSmw8Q1MUTw1mtOgQ0/xKbFSGjMIvIS5u/VLZtxlfcjI1usYYUFcGty4GA6ub3gDr7MalvwKuo6AwZdAdMd+8QgGi0XoP7YL/cd24Vh+KWuX7mPX6oKALy0u7Kwrv5RtW89h7P4FDOz2CL0ShkNOjPHyoD3/NKdIh/7PcYmNsugBgEcUY/ysPnTu2TrGKeRHJOFKVBxkT2g4f+2LsOY541jyK5j1NAy8qO3qF+akdU/gvJsHM/6yPhzYfoLNXxzk8B7/INqVrmSWlcxjU8Ue+sZ8idrx/0iKKSExqyvWhE6Q1t+QqUofALGpEJNqOLZYoyAmGSzhu1NF0zp0aAMFUHq8ynSekBo+YTM0QVBTAV886jl3VkJsp/arTxgTn+LwGlWVsXfjUbYsz6OsxGUqd8yZw7FS99pgEVgOOUm17SPdvpt0+8vEWopQWIixFFPjiqHMlYbgwuFQRMXH4kjLJCq9O7EZmUQnJ4IjwT2Nq7+bHY0Ob6C8w59b7Zb239cTQkTEiG/tC8Zbej0WuzHVF4gdS+GDX0LX4Ya8T69J0P2MtqlnmJHWPZ607vHknpvF2o/yWLc0jzpn4H8YFzbDaDlz2BogzqWJoxgaNm5iLSeIsRykV/ybdO0VA44kElLtRMc7sLtKsNjtSHQidfZEVHQK1sROSFy64f7uSNBagmFORBkoEUkFnsFQvjgG/Fop9Uqw1yekRrfqfqeIeOCHGyNvgNTekP8dbHgV4jMa3mR6cL0RzqMoz3BfHzcvsIFSCo7vMqaoomLBFt1hH4R2h5UxF+UwcHwXVi7c7RcP7HSpcKVS4UrleFE2rPPPF+qwSTW1ylDNt3ACuxwgSiqxW6qIslYTZanCbqkmymEhKjaK2MQo4lLjSMhIIbFHF6JS0rDFxmFxxGhFkhAjogwU8A+gBsgAhgPvi8gGpdSWQIWVeWaC+BTtARZxRMUZOnL9zoPJd5+MkhyQQ+vN5w2F9agpgye8DFd0MtydF7hs/hrYt9JznjPJcIcPxMp/YgpSPea2wA/Lgu+NPV8n6zkAek8JfM9Nb0L5Uc/54MsM8V1fyo8ZZeuJS2s4cvEPX8KRrZ7z7AkkZgxm+pwhjLu0kh82HGPLmp3E2pIp2l9AeWXrfa8U1pPGCYzRWrVKoFolgAvwVukqB0743qEYKMZCLTGWEqxSi1Wc2CxOoqzV2K21YLHhwgZiQ7kqWfbxamKjynHYnca6WWJ3sNgMSW3cLyoCUlMBlV6ivI44Y+0Nr5cZEeO8rACpq+bkxQmZbo9U8RQDcNW6wwEJWASs0ZDcI/CLdflRqHJHGxAFsekQk+T+hPB4WxYVIa/1IhKHEadqiFJqhzvtJeCAUuruQNdE2aPVX27+4OR5l0MrGbh9QavV8YgznQpXO2rWNZNu9nzsEliHT6MJllpbHCUJPShNyKIiNgOnLQaFhWpHMnZnOdFVRogXpzUGp63+iKUquhNKO060OvP+dW6FUiokH0yRNILqB9TVGyc3G4BJwd7AUe33etWiKNUxp4E0HRu7s5xOhdvoVLitWdfVWexURXfiWKehFCf1wmmNQYmFqphOOK3RuKxRKHEbMOVClAtliaRHmiaSejMeY7zuTTFgCnsrInMxgi9i99nU2em4f2RYjUbTPlhdtcRVHCau4jDsD1zGJVbqrFFY62qwqDpcYsNpi6bO6sBpjabOFm38tDpw2qJx2mKpiUqiKjqFquhUqqI7UWex47Lq6f1QJJIMVBngu4EpESj1TlBKzQfmAzhsUSqj4Dssrlo6H11HUuneNqlouKDQIz5NaGNRdViclV7nTqJqy6C2rFn3cVodOG0xuMSGstjcBi0Gp9WBuEdngsJpjabWHke1IwWXabRm/q74T5b4JjT83VJNOtw0fq+mP9urbIh/xyNxDWqwUmqnO+1F4GBDa1DR0dGqqqoqUFarkJ8PFRWt/zmrVq1izJgxp32f7GyICiEpwuXLlzN58uT2rkaLEoltAt2ucEJEQnYNKmL8KZVS5cBC4EERiRORCcDFwEvtWzONRqPRnAoRY6Dc3A7EAEeAV4HbGnIx1zRNhAyuNRpNmBJJa1AopU4Al7R3PTQajUZz+kTaCCqk0SMSjUajCR5toDQajUbTKCJyvohsF5FdIhLQ6aw10AZK0yB6xKfRaETEiiEjdwEwCLhGRAa1xWdH1BpUqGO3g8vVdLnT/xwXjhbYd9hB9U81Go2ZM4FdSqk9ACLyGoaHdKsrG3RoA1VdXa1EpKkAAOGIDbNMZqQQie2KxDaBblc4ESsiq73O57sFDerphlnLIx84/Y2WQdChDRSwVikVcQF/RGS1bld4EIltAt2ucCKINgWaS2mTBQC9BqXRaDSaxsgHeniddwcOtsUHawOl0Wg0msb4DugrIr1EJAq4Gni3LT64o0/xzW+6SFii2xU+RGKbQLcrnGi0TUopp4jMA5YCVuDZtlLoiRixWI1Go9FEFnqKT6PRaDQhiTZQGo1GowlJtIHSaDQaTUjSIQ2UiKSKyNsiUi4ieSJybXvX6VQQkeUiUiUiZe5ju1fete62lYvIIhFJbc+6NoaIzBOR1SJSLSLP++SdKyLbRKRCRJaJSE+vPIeIPCsiJSJyWER+3uaVb4CG2iQi2SKivPqsTETu88oP5TY5ROQZ9/9VqYisE5ELvPLDta8abFc49xeAiCwQkUPu+u0QkVu88kK/v5RSHe7AiBX1HyAeOAsoxojE2+51a2Y7lgO3BEgfjBHq/mx3G18BXmvv+jbSjsswwqQ8CTzvlZ7m7psrgGjgUeAbr/yHgS+BFGAgcBg4v73b00SbsjE2OdoauC6U2xQHPOBugwX4kfv/LDvM+6qxdoVtf7nrNxhwuH8f4K7fqHDpr3b/A7ZDh8UBNUA/r7SXgEfau26n0JaGDNTvgVe8znu725zQ3nVuoj2/83mYzwW+9um7SmCA+/wAcJ5X/kOEmCEO0KamHngh3yaf+m4EZkVCXzXQrojpL6A/cAi4Mlz6qyNO8fUD6pRSO7zSNmC8aYQjD4vIMRFZISKT3WmDMdoEgFJqN26j3A71Ox1821EO7AYGi0gK0NU7n/DqxzwRyReR50QkDSDc2iQiGRj/U1uIoL7yaVc9YdtfIvJPEakAtmEYqA8Ik/7qiAYqHmNo600xkNAOdTldfgXkYIg5zgcWi0hvIqeNjbUj3uvcNy+UOQaMBnpiTLUkAC+788KmTSJix6j3C0qpbURIXwVoV9j3l1Lqdow6TQQWAtWESX91RCWJMiDRJy0RY845rFBKrfI6fUFErgEuJHLa2Fg7yrzOq3zyQhalVBlQrxxd4N6hf0hEEgmTNomIBWNavAaY504O+74K1K5I6C8ApVQd8JWIzAZuI0z6qyOOoHYANhHp65WWi3k4H64oDOXhLRhtAkBEcgAHRtvDCd92xGGsp21RShViTFfkepUPx36sl3KRcGiTiAjwDJABzFJK1bqzwrqvGmmXL2HVXwGw4e4XwqG/2nvhrp0WC1/D8OSLAyYQhl58QDIwHcMDxwZcB5RjLIQOBkowhvRxwAJCdOHW3Rabux0PY7zB1rcp3d03s9xpf8DsafQI8DmGp9EAjC9VSHhQNdKmMe4+sgCdMLxJl4VDm9z1ewr4Boj3SQ/bvmqiXWHbX0BnDGHXeAwNvenuZ8TF4dJf7f5HbKeOSwUWuTtrH3Bte9fpFNqQjqEyXAoUub9c07zyr3W3rRx4B0ht7zo30pYHMN5MvY8H3HlTMRZ3KzG8FrO9rnMAz2IY4wLg5+3dlqbaBFwD/ODul0PAi0BmmLSpp7sdVRjTQPXHdWHeVw22K8z7K91tZIrc9dsEzPHKD/n+0mKxGo1GowlJOuIalEaj0WjCAG2gNBqNRhOSaAOl0Wg0mpBEGyiNRqPRhCTaQGk0Go0mJNEGSqPRaDQhiTZQGk0zcMcGaurY6y77fP3vGo2m+eh9UBpNMxCRsT5Jb2MoPT/glVatlFrnFu5NVEqta6v6aTSRREcUi9VoThml1Dfe5yJSDRzzTXeX3d1mFdNoIhA9xafRtBK+U3xe4cNvFZGH3aG0S91huWNFpI+ILHWHFd8lIjcGuGeuiLwrIoUiUumOAzaxTRum0bQR2kBpNG3PrzECwt0I/Aa4CkOs9G3gfeBSjIiuz4nIySBxIjIS+BpDS3IOhtDnceATERnVlg3QaNoCPcWn0bQ9u5VS9aOjpe4R0PXA9UqpBQAishqYCVyOJ8zBoxgCwFOUUjXuckuBzcB9wCVt1wSNpvXRIyiNpu1Z4nO+zf1zaX2CMmLyHAF6AIhIDDAJeANwiYhNRGwY8b8+Ac5u7UprNG2NHkFpNG1Poc95TSPp0e7fUzFi+tznPvwQEYtSytVSldRo2httoDSa8KAIcAH/wIhJ5Ic2TppIQxsojSYMUEqVi8iXGKG312pjpOkIaAOl0YQPPwe+wHCseAYjwmsaMBKwKqXubs/KaTQtjXaS0GjCBKXUWmA0hmv548BHwN+AoRiGS6OJKLTUkUaj0WhCEj2C0mg0Gk1Iog2URqPRaEISbaA0Go1GE5JoA6XRaDSakEQbKI1Go9GEJNpAaTQajSYk0QZKo9FoNCGJNlAajUajCUn+P+7b7TImZ/s6AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "names = [ \"S\", \"E\", \"I\", \"Q\", \"R\", \"X\", \"V\", \"EI\", \"D\" ]\n",
    "fullnames = [ \"Susceptible (S)\", \"Exposed (E)\", \"Infectious (I)\", \"Quarantined (Q)\", \"Recovered (R)\", \"Expired (X)\", \"V\", \"Exposed + Infectious (E + I)\", \"Daily\" ]\n",
    "\n",
    "fig = plt.figure()\n",
    "plt.plot( statsNoLockDown[SIDX_EI,:], color = \"#ff7f0e\", linestyle = \":\", linewidth = 4 )\n",
    "plt.plot( stats[SIDX_EI,:], label = fullnames[SIDX_EI], color = \"#ff7f0e\", linewidth = 4 )\n",
    "plt.plot( stats[SIDX_X,:], label = fullnames[SIDX_X], color = \"#d62728\", linewidth = 4 )\n",
    "plt.plot( stats[SIDX_I,:] - stats[SIDX_Q], label = \"Non-quaran. Infectious (I - Q)\", color = \"#9467bd\", linewidth = 4 )\n",
    "ax = plt.gca()\n",
    "plt.grid()\n",
    "plt.xlabel(\"Time\", size = 16)\n",
    "plt.ylabel(\"Number of Individuals\", size = 16)\n",
    "ax.tick_params( axis='both', which='major', labelsize = 12 )\n",
    "ax.set_xlim( 0, 300 )\n",
    "plt.legend( loc = \"upper right\", fontsize = 12 )\n",
    "ax2 = ax.twinx()\n",
    "fillLimit = LKP\n",
    "fillLimit[0] = fillLimit[-1] = 0\n",
    "plb.fill( np.arange( T ), fillLimit, facecolor = \"blue\", label = \"lockdown level\", alpha = 0.15 )\n",
    "ax2.set_ylim([0,5])\n",
    "plt.ylabel( \"Lock-down level\", size = 16 )\n",
    "plt.legend( loc = \"center right\", fontsize = 12 )\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "author": {
   "email": "purushot@cse.iitk.ac.in",
   "institution": "IIT Kanpur",
   "name": "Puru"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}