demo_constrained_lockdown.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using ESOP to discover the optimal lock-downs under constraints\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\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 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",
    "T = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "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",
      "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] ) # 12 + 50 = 62\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": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Start on day 30, last for 40 days, at level 3.500000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuUXGWZ7/HvL517BAMkMJpEEjULB1gMl5arOkwYIDBIkINH0LOIyIHlHBiUERUOgzCgM3IRPIyKgxC5HIbLAEpUEGIAGY+AdBBCuJlW0DSJ0phwWRATEp7zx34bK5Xqqp3QVXt31e+zVq+q/da7dz21V7qfPO9+692KCMzMzMpmRNEBmJmZ1eIEZWZmpeQEZWZmpeQEZWZmpeQEZWZmpeQEZWZmpeQEZWZmpeQEZWZmpeQEZWZmpTSy6ABabdKkSTF9+vSiwzAz61iLFi16ISImN+rXcQlq+vTp9PT0FB2GmVnHkvTbPP08xGdmZqXkBGVmZqXkBGVmZqXkBGVmZqXkBGVmZqXUtAQlaZ6k5yUtqWg7T9JiSY9IukvSO1O7JF0qqTe9vnvFPnMlLU0/cyva95D0WNrnUklq1mcxM7PWa+Y086uAbwDXVLRdGBFnAUg6BfgS8GngEGBm+tkLuAzYS9LWwNlANxDAIknzI2JV6nMi8ABwOzAbuKMZH6Svr4+vfe1rzTh0XSNHjuSAAw5g9uzZLX9vM7OiNS1BRcR9kqZXtb1csTmBLOkAzAGuiez+8w9ImijpHcD+wIKIWAkgaQEwW9K9wJYRcX9qvwY4giYlqP7+fr7+9a8349ANXXTRRXz/+99nzpw5hby/mVlRWn4NStJXJC0DPkFWQQFMAZZVdOtLbfXa+2q0D/aeJ0rqkdTT39//1j9Ei/3gBz8oOgQzs5ZreYKKiDMjYhpwHXByaq51/Sg2o32w97w8Irojonvy5Iara5TOyy+/3LiTmVmbKXKpo/8AfkR2jakPmFbx2lRgeWrfv6r93tQ+tUb/ppgyZQoXX3xxsw6/kSeeeIIrrrjize3XXnutZe9tZlYWLU1QkmZGxNK0eTjwVHo+HzhZ0g1kkyReiogVku4E/kXSVqnfQcAZEbFS0iuS9gYeBI4F/q1ZcW+77baceuqpzTr8Ru6+++4NEtSrr77asvc2MyuLpiUoSdeTVT+TJPWRVUqHStoBeAP4LdkMPshm4R0K9AKvAccBpER0HvBQ6nfuwIQJ4O/JZgqOI5sc0ZQJEkUYP378BtuuoMysEzVzFt8xNZqvHKRvACcN8to8YF6N9h5g57cSY1lNmDBhg20nKDPrRF5JooRcQZmZOUGVUnWC8jUoM+tETlAl5ArKzMwJqpRqJajsMp2ZWedwgiqhUaNGMWrUqDe3169fz+uvv15gRGZmrecEVVIe5jOzTucEVVKeKGFmnc4JqqRcQZlZp3OCKiknKDPrdE5QJeXVJMys0zlBlZQrKDPrdE5QJeVJEmbW6ZygSsoVlJl1OieoknKCMrNO5wRVUp4kYWadzgmqpHwNysw6nRNUSXmIz8w6nRNUSTlBmVmnc4IqKScoM+t0TlAl5UkSZtbpnKBKypMkzKzTOUGVlIf4zKzTOUGVlBOUmXU6J6iScoIys07nBFVSniRhZp3OCaqkPEnCzDqdE1RJeYjPzDqdE1RJOUGZWadzgiopJygz63ROUCU1atQourq63txet24da9euLTAiM7PWcoIqKUmeyWdmHc0JqsQ8zGdmnaxpCUrSPEnPS1pS0XahpKckLZb0PUkTK147Q1KvpKclHVzRPju19Uo6vaJ9hqQHJS2VdKOk0c36LEVxgjKzTtbMCuoqYHZV2wJg54jYBfgVcAaApB2Bo4Gd0j7fktQlqQv4JnAIsCNwTOoLcD5wSUTMBFYBxzfxsxTCCcrMOlnTElRE3AesrGq7KyLWpc0HgKnp+RzghohYExHPAL3AnumnNyJ+ExFrgRuAOZIEzAJuTvtfDRzRrM9SFH9Z18w6WZHXoD4F3JGeTwGWVbzWl9oGa98GeLEi2Q201yTpREk9knr6+/uHKPzm8yQJM+tkhSQoSWcC64DrBppqdIvNaK8pIi6PiO6I6J48efKmhlsYD/GZWScb2eo3lDQXOAw4ICIGkkofMK2i21RgeXpeq/0FYKKkkamKquzfNpygzKyTtbSCkjQb+CJweERU/rWdDxwtaYykGcBM4BfAQ8DMNGNvNNlEivkpsd0DHJX2nwvc1qrP0SpOUGbWyZo5zfx64H5gB0l9ko4HvgFsASyQ9IikbwNExOPATcATwI+BkyJifaqOTgbuBJ4Ebkp9IUt0/yipl+ya1JXN+ixF8SQJM+tkTRvii4hjajQPmkQi4ivAV2q03w7cXqP9N2Sz/NqWJ0mYWSfzShIl5iE+M+tkTlAl5gRlZp3MCarEfA3KzDqZE1SJuYIys07mBFViniRhZp3MCarEXEGZWSdzgioxJygz62ROUCXmSRJm1smcoErMFZSZdbKWLxZr+VVPklixYgUXX3zxRv123HFHDj74YLLbZJmZtQcnqBKrrqBeeOEFPve5z9Xse/bZZ3POOee0ICozs9bwEF+JbbHFFrn7XnvttU2MxMys9ZygSmzLLbdk1qxZufquWrWqydGYmbWWh/hK7tZbb2XevHksW7Zsg/b169dz6aWXvrntCRRm1m7055vadobu7u7o6ekpOoy3LCIYMWLDAnj9+vUbtZmZlY2kRRHR3aif/5oNU5IYN27cBm2rV68uKBozs6HnBDWMOUGZWTtzghrG/EVeM2tnTlDDmCsoM2tnTlDDmCsoM2tnTlDDmCsoM2tnTlDDmCsoM2tnuRKUpHGSdmh2MLZpXEGZWTtrmKAkfRh4BPhx2t5V0vxmB2aNuYIys3aWp4I6B9gTeBEgIh4BpjcvJMvLFZSZtbM8CWpdRLzU9Ehsk7mCMrN2lmex2CWSPg50SZoJnAL8vLlhWR6uoMysneWpoP4B2AlYA1wPvAx8tplBWT7VCcoVlJm1k4YVVES8BpyZfqxEqof4XEGZWTtpmKAk3QNsdE+OiMh3Jz1rGldQZtbO8lyDOq3i+VjgvwHrmhOObQpXUGbWzvIM8S2qavp/kn7apHhsE7iCMrN2lmeIb+uKzRHAHsBfNC0iy80VlJm1szyz+BYBPenxfuBzwPGNdpI0T9LzkpZUtH1U0uOS3pDUXdX/DEm9kp6WdHBF++zU1ivp9Ir2GZIelLRU0o2SRuf4LG3FFZSZtbOGCSoiZkTEu9PjzIg4KCJ+luPYVwGzq9qWAEcC91U2StoROJpsOvts4FuSuiR1Ad8EDgF2BI5JfQHOBy6JiJnAKnIkzXbjCsrM2tmgQ3ySjqy3Y0Tc2uD1+yRNr2p7Mh27uvsc4IaIWAM8I6mXbHklgN6I+E3a7wZgjqQngVnAx1Ofq8mWZLqsXkztxhWUmbWzetegPlzntQDqJqhNNAV4oGK7L7UBLKtq3wvYBngxItbV6L8RSScCJwK8613vGqKQi+cKysza2aAJKiKOa2EcG5VUZEmw1hBk1OlfU0RcDlwO0N3dPWi/4cZLHZlZO8vzPSgk/R3Z9aGxA20Rce4QxtEHTKvYngosT89rtb8ATJQ0MlVRlf07hheLNbN2lud+UN8GPka2Jp+AjwLbD3Ec84GjJY2RNAOYCfwCeAiYmWbsjSabSDE/IgK4Bzgq7T8XuG2IYyo9V1Bm1s7yTDPfNyKOBVZFxD8D+7BhVVOTpOvJpqXvIKlP0vGSPiKpLx3jR5LuBIiIx4GbgCfIbox4UkSsT9XRycCdwJPATakvwBeBf0wTKrYBrsz/sduDKygza2fKipE6HaQHI2IvSQ+QTRH/I7AkTe8edrq7u6Onp6foMIZERDBixIb/x1i/fv1GbWZmZSJpUUR0N+qX5y/ZDyVNBC4EHgaeJbvthhVMkof5zKxt5VmL77z09BZJPwTG+g675TFu3LgNktLq1auZMGFCgRGZmQ2NPJMkHpX0vyW9JyLWODmVi69DmVm7yjPEdzjZ7TVukvSQpNMktc+3XYc5D/GZWbvKsxbfbyPigojYg2xpoV2AZ5oemeXiCsrM2lXeL+pOB/472feh1gNfaF5ItilcQZlZu8pzP6gHgVFk31P66MDCrVYOrqDMrF3lqaDmRsRTTY/ENosrKDNrV3muQTk5lZgrKDNrV15yYJhzBWVm7coJaphzBWVm7SrPF3XHSzpL0nfS9kxJhzU/NMvDFZSZtas8FdR3gTVkK5BDdu+mLzctItskrqDMrF3lSVDviYgLgNcBImI1te9oawVwBWVm7SpPgloraRzpluqS3kNWUVkJVCcoV1Bm1i7yfA/qHLKbCE6TdB2wH/DJJsZkm6B6iM8VlJm1izy327hL0iJgb7Khvc9ExAtNj8xycQVlZu0qz1JH88luUDg/Il5tfki2KVxBmVm7ynMN6mvAB4EnJP2npKMkjW1yXJaTKygza1d5hvh+CvxUUhcwCzgBmAds2eTYLAdXUGbWrvLebmMc8GGy223sDlzdzKAsP1dQZtau8lyDuhHYi2wm3zeBeyPijWYHZvm4gjKzdpWngvou8PGIWN/sYGzTuYIys3Y1aIKSNCsi7gbGA3OkDRePiIhbmxyb5eAKyszaVb0K6q+Bu8muPVULwAmqBLzUkZm1q0ETVEScnZ6eGxHPVL4maUZTo7LcvFismbWrPN+DuqVG281DHYhtnrFjN/xK2po1a3jjDc9hMbPhr941qPcBOwFvl3RkxUtbAv6ibklIYty4cRsM7a1evZoJEyYUGJWZ2VtX7xrUDsBhwEQ2vA71CtmXda0knKDMrB3VuwZ1G3CbpH0i4v4WxmSbaPz48axcufLNbV+HMrN2kOd7UJ+W9GREvAggaSvgaxHxqeaGZnlVz+Tba6+9GDVq1JC+x9ve9jaOO+44TjvtNKq/cmBm1gx5EtQuA8kJICJWSdqtiTHZJqqeyff73/++Ke/zhS98gf3224999923Kcc3M6uUZxbfiFQ1ASBpa/ItkTRP0vOSllTuK2mBpKXpcavULkmXSuqVtFjS7hX7zE39l0qaW9G+h6TH0j6XqoP/W7/99tu37L0ee+yxlr2XmXW2vLfb+Lmk8ySdC/wcuCDHflcBs6vaTgcWRsRMYGHaBjgEmJl+TgQugzeT4dlkawHuCZxdkSwvS30H9qt+r45x1llnsd1227XkvdauXduS9zEzy3O7jWsk9ZDdakPAkRHxRI797pM0vap5DrB/en41cC/wxdR+TUQE8ICkiZLekfouiIiVAJIWALMl3QtsOTB5Q9I1wBHAHY3iakfd3d0sX76c5557juwUDp3zzjuPK6644s3tNWvWDOnxzcwGk+t2G8DWwKsR8V1JkyXNqF5dIqftImIFQESskLRtap8CLKvo15fa6rX31WjvWCNGjGDatGlDftxJkyZtsO0KysxapeEQn6SzyaqcM1LTKOD/DnEcta4fxWa01z64dKKkHkk9/f39mxliZxozZswG266gzKxV8lyD+ghwOPAqQEQsB7bYzPf7Qxq6Iz0+n9r7gMr//k8Fljdon1qjvaaIuDwiuiOie/LkyZsZemcaPXr0BttOUGbWKnkS1Np0bSgAJL2VJQrmAwMz8eYCt1W0H5tm8+0NvJSGAu8EDpK0VZoccRBwZ3rtFUl7p9l7x1Ycy4ZQdQXlIT4za5U816BukvTvwERJJwCfAr7TaCdJ15NNcpgkqY9sNt5X0/GOB34HfDR1vx04FOgFXgOOA4iIlZLOAx5K/c4dmDAB/D3ZTMFxZJMjOnKCRLN5iM/MipJnFt9Fkg4EXiZbn+9LEbEgx37HDPLSATX6BnDSIMeZB8yr0d4D7NwoDntrqof4XEGZWavkmsWXElLDpGTtxxWUmRVl0GtQkn6WHl+R9HKNn2ck/a/WhWpF8CQJMytKvdXMP5Aea87Yk7QN2aoS32pOaFYGniRhZkXJNcSX1sb7ANlMvp9FxC8j4o+S9m9mcFY8D/GZWVHyfFH3S2TLEm0DTAKukvRPkK0G0dzwrGieJGFmRclTQR0D7BYRfwKQ9FXgYeDLzQzMysEVlJkVJc8XdZ8FxlZsjwF+3ZRorHRcQZlZUQatoCT9G9k1pzXA42kl8QAOBH7WmvCsaK6gzKwo9Yb4etLjIuB7Fe33Ni0aKx0nKDMrSr1p5lcDSBoLvJesevr1wLUo6wwe4jOzotT7ou5ISReQrRx+NdktNpZJukDSqFYFaMVyBWVmRak3SeJCshsVzoiIPSJiN+A9wETgolYEZ8VzBWVmRamXoA4DToiIVwYaIuJlslXED212YFYOrqDMrCj1ElSkVcarG9dT5+611l68Fp+ZFaVegnpC0rHVjZL+B/BU80KyMqm1Fl+N/7eYmQ25etPMTwJulfQpsqnmAbyf7AaBH2lBbFYCXV1ddHV1sX79egAignXr1jFqlOfJmFlz1Ztm/hywl6RZwE6AgDsiYmGrgrNyGD16NKtXr35ze+3atU5QZtZ0ee6oezdwdwtisZIaM2bMBglqzZo1TJgwocCIzKwT5FmLzzqcJ0qYWRGcoKwh37TQzIrgBGUN+btQZlYEJyhryKtJmFkRnKCsIVdQZlYEJyhryBWUmRXBCcoacgVlZkVwgrKGnKDMrAhOUNaQh/jMrAhOUNaQKygzK4ITlDXkCsrMiuAEZQ25gjKzIjhBWUNei8/MiuAEZQ15LT4zK4ITlDXkIT4zK0IhCUrSZyQtkfS4pM+mtq0lLZC0ND1uldol6VJJvZIWS9q94jhzU/+lkuYW8Vk6gSdJmFkRWp6gJO0MnADsCfwVcJikmcDpwMKImAksTNsAhwAz08+JwGXpOFsDZwN7pWOdPZDUbGi5gjKzIhRRQf0l8EBEvBYR64CfAh8B5gBXpz5XA0ek53OAayLzADBR0juAg4EFEbEyIlYBC4DZrfwgncKTJMysCEUkqCXAhyRtI2k8cCgwDdguIlYApMdtU/8pwLKK/ftS22DtG5F0oqQeST39/f1D+mE6gSdJmFkRWp6gIuJJ4HyyiufHwKPAujq7qNZh6rTXes/LI6I7IronT568iRGbh/jMrAiFTJKIiCsjYveI+BCwElgK/CEN3ZEen0/d+8gqrAFTgeV12m2IeZKEmRWhqFl826bHdwFHAtcD84GBmXhzgdvS8/nAsWk2397AS2kI8E7gIElbpckRB6U2G2KuoMysCCMLet9bJG0DvA6cFBGrJH0VuEnS8cDvgI+mvreTXafqBV4DjgOIiJWSzgMeSv3OjYiVrfwQncIVlJkVoZAEFREfrNH2R+CAGu0BnDTIceYB84Y8QNuAKygzK4JXkrCGnKDMrAhOUNaQh/jMrAhOUNaQKygzK4ITlDXkCsrMiuAEZQ25gjKzIjhBWUNei8/MiuAEZQ15LT4zK4ITlDXkIT4zK4ITlDXkSRJmVgQnKGvIFZSZFcEJyhryJAkzK4ITlDVUa5JEtkSimVnzOEFZQyNGjGDkyA3XFX799dcLisbMOoUTlOXiiRJm1mpOUJaLJ0qYWas5QVkurqDMrNWcoCwXV1Bm1mpOUJaLE5SZtZoTlOXiIT4zazUnKMvFFZSZtZoTlOXiCsrMWs0JynJxBWVmreYEZbl4PT4zazUnKMvFNy00s1ZzgrJcPMRnZq3mBGW5eJKEmbWaE5Tl4grKzFrNCcpy8SQJM2s1JyjLxZMkzKzVnKAsFw/xmVmrOUFZLp4kYWat5gRlubiCMrNWKyRBSTpV0uOSlki6XtJYSTMkPShpqaQbJY1Ofcek7d70+vSK45yR2p+WdHARn6VTuIIys1ZreYKSNAU4BeiOiJ2BLuBo4HzgkoiYCawCjk+7HA+sioj3ApekfkjaMe23EzAb+JakrlZ+lk7iCsrMWq2oIb6RwDhJI4HxwApgFnBzev1q4Ij0fE7aJr1+gCSl9hsiYk1EPAP0Anu2KP6O4wRlZq3W8gQVEc8BFwG/I0tMLwGLgBcjYl3q1gdMSc+nAMvSvutS/20q22vsY0PMQ3xm1mpFDPFtRVb9zADeCUwADqnRNQZ2GeS1wdprveeJknok9fT392960OYKysxabmQB7/m3wDMR0Q8g6VZgX2CipJGpSpoKLE/9+4BpQF8aEnw7sLKifUDlPhuIiMuBywG6u7trJjGrr7qCuuuuu5g1a1ZB0ZhZ0SSxcOHCpr5HEQnqd8DeksYDq4EDgB7gHuAo4AZgLnBb6j8/bd+fXr87IkLSfOA/JF1MVonNBH7Ryg/SSaorqP7+fu65556CojGzomVTAZqr5QkqIh6UdDPwMLAO+CVZdfMj4AZJX05tV6ZdrgSuldRLVjkdnY7zuKSbgCfScU6KiPUt/TAdZLvttis6BDPrMIrorBGv7u7u6OnpKTqMYWfdunUceOCB3HvvvUWHYmYlIIk33nhjc/ddFBHdjfoVMcRnw9DIkSP5yU9+wuLFi1m1alXR4ZhZB3CCsty6urrYbbfdig7DzDqE1+IzM7NScoIyM7NScoIyM7NScoIyM7NScoIyM7NScoIyM7NScoIyM7NS6riVJCT1A7/dxN0mAS80IZxWGK6xD9e4YfjG7rhbb7jG/lbj3j4iJjfq1HEJanNI6smzLEcZDdfYh2vcMHxjd9ytN1xjb1XcHuIzM7NScoIyM7NScoLK5/KiA3gLhmvswzVuGL6xO+7WG66xtyRuX4MyM7NScgVlZmal5ARlZmal5ATVgKTZkp6W1Cvp9KLjGYykaZLukfSkpMclfSa1by1pgaSl6XGromOtRVKXpF9K+mHaniHpwRT3jZJGFx1jLZImSrpZ0lPp3O8zHM65pFPTv5Mlkq6XNLas51zSPEnPS1pS0VbzHCtzafp9XSxp95LFfWH6t7JY0vckTax47YwU99OSDi4m6jdj2Sj2itdOkxSSJqXtpp1zJ6g6JHUB3wQOAXYEjpG0Y7FRDWod8LmI+Etgb+CkFOvpwMKImAksTNtl9BngyYrt84FLUtyrgOMLiaqx/wP8OCLeB/wV2Wco9TmXNAU4BeiOiJ2BLuBoynvOrwJmV7UNdo4PAWamnxOBy1oUYy1XsXHcC4CdI2IX4FfAGQDpd/VoYKe0z7fS35+iXMXGsSNpGnAg8LuK5qadcyeo+vYEeiPiNxGxFrgBmFNwTDVFxIqIeDg9f4XsD+UUsnivTt2uBo4oJsLBSZoK/B1wRdoWMAu4OXUpa9xbAh8CrgSIiLUR8SLD4JyT3U17nKSRwHhgBSU95xFxH7CyqnmwczwHuCYyDwATJb2jNZFuqFbcEXFXRKxLmw8AU9PzOcANEbEmIp4Besn+/hRikHMOcAnwBaBydl3TzrkTVH1TgGUV232prdQkTQd2Ax4EtouIFZAlMWDb4iIb1NfJ/tG/kba3AV6s+EUu63l/N9APfDcNT14haQIlP+cR8RxwEdn/glcALwGLGB7nfMBg53g4/c5+CrgjPS993JIOB56LiEerXmpa7E5Q9alGW6nn5Ut6G3AL8NmIeLnoeBqRdBjwfEQsqmyu0bWM530ksDtwWUTsBrxKyYbzaknXa+YAM4B3AhPIhmmqlfGcNzIs/u1IOpNsWP66gaYa3UoTt6TxwJnAl2q9XKNtSGJ3gqqvD5hWsT0VWF5QLA1JGkWWnK6LiFtT8x8Gyu30+HxR8Q1iP+BwSc+SDaHOIquoJqbhJyjvee8D+iLiwbR9M1nCKvs5/1vgmYjoj4jXgVuBfRke53zAYOe49L+zkuYChwGfiD9/EbXscb+H7D80j6bf1anAw5L+gibG7gRV30PAzDS7aTTZRcz5BcdUU7pucyXwZERcXPHSfGBuej4XuK3VsdUTEWdExNSImE52fu+OiE8A9wBHpW6lixsgIn4PLJO0Q2o6AHiCkp9zsqG9vSWNT/9uBuIu/TmvMNg5ng8cm2aW7Q28NDAUWAaSZgNfBA6PiNcqXpoPHC1pjKQZZBMOflFEjLVExGMRsW1ETE+/q33A7ul3oHnnPCL8U+cHOJRsts2vgTOLjqdOnB8gK6sXA4+kn0PJrucsBJamx62LjrXOZ9gf+GF6/m6yX9Be4D+BMUXHN0jMuwI96bx/H9hqOJxz4J+Bp4AlwLXAmLKec+B6smtlr6c/jMcPdo7Jhpu+mX5fHyObqVimuHvJrtcM/I5+u6L/mSnup4FDynbOq15/FpjU7HPupY7MzKyUPMRnZmal5ARlZmal5ARlZmal5ARlZmal5ARlZmal5ARlVoekf5W0v6QjVGc1e0nHppXBH5f0hKTTGhz3nIE+kq6SdFS9/hX77S9p302If7qkPkkjqtofkbRnjTieSa89LGmfGu2PSPp5av+kpP7U9pSkU/PGZZaHE5RZfXuRrWn418B/1eog6RDgs8BBEbET2WoSLzUpnv3JVn3IJSKeJfvezQcH2iS9D9giImp9EfTzEbEr2ZJN/17dnn4q3//G1H8/4My02rXZkHCCMqsh3bdnMfB+4H7gfwKXSaq1FtkZwGkRsRwgIv4UEd9JxzlB0kOSHpV0S1rTrN77fjVVYIslXVT12nTg08CpqWr5oKTtJS1M/RdKeleNw15PtkrHgKNTWz33Ae9t0OdNEfFHsi+hFrJyuLUnJyizGiLi82RJ6SqyJLU4InaJiHNrdN+ZbDXwWm6NiPdHxMC9oga9x5KkrYGPADtFdr+gL1fF9CzwbbJ7Nu0aEf8FfIPsVge7kC08emmNQ98EHFGxzt7HyNY9rOfDZKsCDLiwYojvuurOKTGOJVtRw2xIjGzcxaxj7Ua2HM37yNaq2xw7S/oyMBF4G3Bnnb4vA38CrpD0I+CHOY6/D3Bken4tcEF1h4j4vaTHgQMk/QF4PSI2ulNqcqGkfyK7jUhlMv18RNxco//HJP0NsANwQkT8KUfMZrk4QZlVkbQrWeU0FXiB7IZ+kvQIsE9ErK7a5XFgD+DuGoe7CjgiIh6V9Emya0g1RcQ6SXuSLd56NHAy2erum2KwtcsGhvn+QP3hvcES0WBujIiT04SKH0m6I7IFRM3eMg/xmVWJiEfShf9fATuSJZ6D07BadXIC+FfggnTrAdKK1Kek17YAVqRboXyi3vume3m9PSIRxJ/wAAAA1UlEQVRuJ5t0sWuNbq+kYw74OX++vvQJ4GeDHP4WssWD8wzvbbKIuJ+sgvvMUB/bOpcrKLMaJE0GVkXEG5LeFxGDDvFFxO2StgN+km5fEcC89PJZZLMAf0t2TWeL2keB9NptksaSrRBda9r2D4CbJc0B/gE4BZgn6fNkw3LHDRLji5IeILsT7TN1YhjMwNDfgFq3Iz+f7B5B/xIRr2zGe5htwKuZm5lZKXmIz8zMSskJyszMSskJyszMSskJyszMSskJyszMSskJyszMSskJyszMSun/A1X4YLIqNM9vAAAAAElFTkSuQmCC\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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15630 Infected\n",
      "2266 Expired\n",
      "Peak 7353\n",
      "Recovery time 82.490871\n",
      "Expiry time 21.882171\n",
      "Quarantine time 11.201635\n",
      "14928 Quarantined\n",
      "Objective acheived is 7353.000000 (fEpi) + 560.000000 (fEco) = 7913.000000 (overall)\n"
     ]
    }
   ],
   "source": [
    "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 4(b) in the accompanying paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd4FFXXwH83vRNIaAkhoYOUqCCoCIkIqDT5RFEpKgqoKIiKLxaQKiqC8oqVjlIEFaTYUUPRF0SlKB2UAAmQBAikl835/pjNZnezSTZls0mY3/PMQ+bOLWd22Tlz7z1FiQg6Ojo6OjpVDRdnC6Cjo6Ojo2MLXUHp6Ojo6FRJdAWlo6Ojo1Ml0RWUjo6Ojk6VRFdQOjo6OjpVEl1B6ejo6OhUSdycLYCOjo6OTtVGKXUSSAEMQK6IdKqMcXUFpaOjo6NjD7eKSFJlDqgv8eno6OjoVEnU1RxJwsXFRby9vZ0tRoWTl5eHi0vNe/eoifdVE+8J9PuqTqSnpwvwp1nRAhFZYF5HKfUvcAkQ4CPr647iqlZQXl5ekpmZ6WwxKpyYmBiio6OdLUaFUxPvqybeE+j3VZ1QSqWLiG8JdUJEJF4pVQ/4ARgrItscLVvNehXQ0dHR0alwRCTe+G8CsB7oXBnj6gpKR0dHR6dIlFK+Sin//L+B3sDflTG2bsWno6Ojo1Mc9YH1SinQdMYqEfm2MgbWFZSOjo6OTpGIyD9ApDPG1o0kdCMJp5OXl8eZM2dIS0srtl5mZiZeXl6VJFXlUBPvCfT7qmr4+vrSqFEjmxaI9hhJOAt9BqXjdJKSklBK0apVq2JNeFNSUvD3969EyRxPTbwn0O+rKpGXl0dcXBxJSUnUq1fP2eKUCl1BVWVE4OxeOH8ADDnQMBJCrgNtLbjGkJycTERERI3zL9HRqQq4uLhQv359YmNjdQWlU0HE74UNT8H5vyzLG3WGuxdAnSbOkcsBGAwG3N3dnS2Gjk6Nxd3dndzcXGeLUWr0V9aqSNwfsLx/YeUEcOY3+CgKDm6sfLkciKphs0IdnapEdf196QqqqpFyDlY/AFlXiq6TdRnWDoevn4fs9MqTTUdHR6cS0RVUVSIvD74YCannLcvDu2p7T9b8tgA+6g7JpytHPp1qxbJly7jlllucKsP69esJCwvDz8+PPXv2VFi/27dvp1WrVhXWn70sWbKE8ePHV/q4+WzcuJH777/faeNXNrqCqkr8sRRObrcsu3ksPPwVjPoZbn8NXKz2ai4c02ZcmcXMuHTKTEREBN7e3vj5+ZmOp556ytliOY2YmBgaNWpkd/0JEybw7rvvkpqaynXX2XjJshOlFMePHzedd+vWjSNHjpS5v7KQnZ3Nm2++yfPPPw/AyZMnUUpZ/N/w8/NjzZo1FTpuREQEW7ZsAWDAgAH8/fff7N+/v0LHqKroCqqqkH4Rtky1LGsSBT2naVZ7SsFNY+DR7yAw3LLe+b9gUU+4fKbSxL2a2LRpE6mpqabj3XffdbZIFU50dDQxMTEV3m9sbCxt27at8H6dwYYNG2jZsiWhoaEW5cnJyRb/P+67774S+1q2bBkPP/xwmeR44IEHWLCgUoKJOx1dQVUVdrxlue/k4Qd3vQcurpb1QjvCE79Ai96W5UlHYO2Dmjl6dWVqrWIP/7mNSqxTqqOcPPHEE9xzzz2m84kTJ3LbbbchIqaZxqxZswgODiYiIoKVK1ea6l6+fJkHH3yQJk2aEB4ezsyZM8nLywPg+PHjREVFUatWLYKDgy0eeIcPH6ZXr17UqVOHVq1asXbtWtO1CxcuMGDAAAICAujcuTMnTpwo9z2WRHR0NJMnT6Zr1674+/vTu3dvkpKSyMrKws/PD4PBQGRkJM2aNQMgPj6eQYMGUbduXZo0acI777xj6stgMDBr1iyaNWuGv78/HTt25PTp03Tv3h2AyMhI0wzFeiZ36NAhoqOjCQwMpG3btmzcuNFCxkWLFpnOzZc+RYRnnnmGevXqUatWLTp06MDff9sOM/fNN9/QtWvXivvwykh0dDRfffWVs8WoFHQFVRW4Eg+/LbQs6/48BIbZru/pD/cs1ZSVOXF/wNbZjpGxJiN5kJUKl2Lh3N9wdh8kHoHUhGKbzZ07l/3797Ns2TK2b9/O4sWLWb58ucli6ty5cyQlJREXF8fy5csZPXq0aVlq7NixXL58mf3797N161Y+/vhjli5dCsDkyZPp3bs3ly5d4syZM4wdOxaAtLQ0evXqxZAhQ0hISGD16tWMGTOGAwcOAPDkk0/i5eXF2bNnWbJkCUuWLHHUJ2bBqlWrWLp0KQkJCWRnZzNnzhw8PT1JTU0FYN++fZw4cYK8vDz69+9PZGQkcXFx/Pjjj8ybN4/vvvsOgLfeeovVq1fz9ddfc+XKFZYsWYKPjw/btm0z9WNrhpKTk0P//v3p3bs3CQkJzJ8/n6FDh9q1BPj999+zbds2jh49SnJyMmvWrCEoKMhm3b/++osWLVqU56OqENq0acPJkye5cqXmL+tXqoJSSqVaHQal1Hyz67cppQ4rpdKVUj8rpcLNrnkqpZYopa4opc4ppZ616rvItlWebW9CrlnIJf+G0OWx4tt4+ml7U817WZb/Mg/v9PiKl7EmIqItrZ4/oO3lZVyEvBxNYeWkw5U4MGQzcOBdBAbWMh6BLFyovUz4+PiwYsUKnn32WYYNG8b8+fML7c/MmDEDT09PoqKi6Nu3L2vXrsVgMLBmzRpee+01/P39iYiI4LnnnuOTTz4BNJ+V2NhY4uPj8fLyMr3tb968mYiICEaMGIGbmxvXX389gwYN4vPPP8dgMPDFF18wffp0fH19adeuHQ899FClfIwjRoygZcuWeHt7M3jwYPbu3Wuz3u7du0lMTOSVV17Bw8ODpk2bMmrUKD799FMAFi1axMyZM2nVqhVKKSIjI4tUFubs3LmT1NRUXnjhBTw8POjRowf9+vVj9erVJbZ1d3cnJSWFw4cPIyK0adOGhg0b2qybnJyMn59fofLg4GACAwNNx6FDh0octzzkR7JITk526DhVgUpVUCLil3+gRcjNAD4DUEoFA+uAyUAd4HfAfLdxKtACCAduBf6jlLrDzrZVl8tn4M+PLcui/gPudmT6dffWnHb96heUGbJpdWR+9V7qqyzOH4DkWMgr3oHxy0VzST4YYzx+ZtS9t5vadO7cmaZNmyIiDB482KJd7dq18fUtCHEWHh5OfHw8SUlJZGdnEx4ebnEtLi4OgNmzZyMidO7cmbZt25pmQrGxsezatcviYbhy5UrOnTtHYmIiubm5hIWFWfRZHOb97Nixg379+pnOX3/9dTs+QI0GDRqY/vbx8THNnKzJV7rm486aNYvz5zWr1dOnT5uWAktDfHw8YWFhFpFIzD/P4ujRowdPPfUUTz75JPXr12f06NFFzkxq165t896SkpJITk42HW3atLHZfsyYMab7HjNmDKtWrTKdd+jQwc671cItgfb91XScGUniHiAByDdbuxs4ICL5CmsqkKSUai0ih4EHgREicgm4pJRaCDwMfGtH26rLvk8tH5C1I+C64fa396kDvV+FdSNNRYGXD2oGF7e/WmFiVgpTLxd7uULioGVcguRT2iwpr4xKPP2CNvPy9Oe95Z+TlZVFSEgIs2fP5sUXXzRVu3TpEmlpaSYlderUKdq1a0dwcLBplpSvUE6dOmXafG/QoIFplrZjxw569uxJ9+7dCQsLIyoqih9++KGQSAaDATc3N06fPk3r1q1NfRaH+Rt4dHQ0U6dOdWiQ4bCwMJo0acKxY8eKvH7ixAnatWtXqn5DQkI4ffq0RTr2U6dO0bJlS0ALlJqeXuAveO7cOYv248aNY9y4cSQkJDB48GDefPNNZsyYUWicDh06WFgSlpb333+f999/H9D2wWJiYli2bFmp+zl06BAREREEBASUWZbqgjMV1EPAx1IQTr0tsC//ooikKaVOAG2VUueBEPPrxr8HltQWsFBQSqnRwGgANzc3h1gu2Y0IN+xeinkY4WNBvYnb/ksp+wmmQ+1rqXPJbGnlf+9yNCmX+JA7qnzsvlq1apneCovDYDDYVa8oXHMz8M6Io6hPI9fVhxx3fwyunnhkJ+OeY3smoCEcPfgXk6a/yk+fL8atVgjde/eje/fudOjQwfRAfPHFF5kyZQq///47mzdv5j//+Q/p6en83//9HxMnTuT999/n1KlTzJkzh3HjxpGSksL69evp3LkzoaGheHh4oJQiIyODqKgoJk6cyIIFC0zGGfv378fPz49WrVoxYMAAXn75ZVOfy5Yto3HjxnZ/tunp6SXWTU9PR0RM9QwGA5mZmabzzMxMDAaDxXeVmppKSkoKbdq0wc/Pj+nTp/P444/j4eHBkSNHyMjIoGPHjgwfPpyXX36Zxo0b06xZMw4cOEDDhg0JCgqiXr16/P3339SvX7+QHNdccw3e3t7MmDGDsWPHsnPnTjZt2kRMTIzp+meffcZ9993H2bNnWbhwIfXq1SMlJYU//viDvLw8rr32WvLy8nBzcyM3N9fm53DrrbeyePFii/sC7cXJza10j9LMzExycnLs+m5ExOK7+f7777nttttK/VvIzMx07vOuLIhIpR9AY8AANDErWwy8blXvF7RZUhgggJfZtV7AyZLaFieHp6enOJW4PSJTAgqOaXVEUhPL1ldKgsjcNpb9TQkQWdZfJPFYxcpdwRw8eNCueleuXCn7IJlXROL3icT9aXnE7xW5cEIkK61wm7w8CW/cWLy8PMXX10d8fbzF18dbBt5xq+TE/iY3XNtWXntxrKmv9996Xdq1ayeZmZny888/S2hoqMycOVOCgoIkLCxMPv74Y1PXFy9elKFDh0pQUJA0atRIpk2bJgaDQUREnn/+eQkJCRFfX19p2rSpfPTRR6Z2hw8flj59+khwcLDUqVNHbr31VtmzZ4+IiCQkJEjfvn3F399fbrjhBpk0aZJ07drVro8nKipKfv755xLr5d+XebuFCxeazpcuXSpdu3Y1fVeAHDtW8P8vLi5O7r//fqlfv74EBgZKly5d5IcffhARkdzcXJkxY4ZERESIn5+fdOrUSU6fPi0iIh988IE0aNBAatWqJWvWrCkkx99//y3du3eXgIAAadOmjaxbt850LTExUXr16iV+fn5y8803y5QpU0yfy5YtW6R9+/bi6+srQUFBMmTIEElJSbF579nZ2dKoUSOJi4sTEZF///1XAPH19bU45s6dW+LnuHTpUnnooYdKrCciEh4ebvqMRETatWsne/futautOUX9zoA0cYIesOdwSj4opdQkoJeIRJmV/RdwF5ExZmV/oe09/QRcBOqLSILx2iBgqoi0L66tiHxRlBxOzwf17Yuw8/2C85Z3wJBybJ2d+QOW3gGGbMtyVw+45Vm45Rlwr3q5bA4dOlTkur05ZV7iS7+oLeth9X89IAR864Kycys2L0+L8pGWoC0R2sI3GAIaEbN1K8OGDePMmeJ906pj+gZ7qKn39c477/DPP/8wb948p4y/adMmPvnkEwv3Ansp6ndWlfNBOcvM/EFguVXZAcyyNiqlfIFmaHtLl4CzWGZ1jDS2KbZthUteURhy4a/PLcs6lOzgVyyNOsLAD8hTVtEmDNmw9XVYEAVJZV9Dr5akntcMIayVk38DzbjEXuUE4OICAQ2hfnsIbgV+9cB6wTAtSTt0aiQjRoxwmnIC6N+/f5mUU3Wl0hWUUupmIBSj9Z4Z64F2SqlBSikv4BVgvxQYOXwMTFJK1VZKtQZGAcvsbFv1+DdGexPPxzMAWt1Z/n7b38PuG96BptGFryUehoU94Mi35R+nOpB+UfMxsyYgBPwaFC63FxcX8PCBgFAIagbKypn6ShxkFG/woaOjUzLOmEE9BKwTEYsdPhFJBAYBrwKXgC6AeVTEKcAJIBbYCrwpIt/a2bbqcWC95fk1A+wzLbeDDJ8QGP4lDFoMvlYJyrIuw+r74JuJ2iyuppKbBZetg+gqLUyUX/2KMxzx9IfgFlZKSoiObMyZQ39UzBg6OlcpdpmeGP2MfETklFnZY0A74DsR2WzvgCJSpAeqiGwBWhdxLQt4xHiUqm2Vw5ADh6w+snaDKnYMpaD9PdC8J3w9Af6ymrDu+hCyU2HAu1Xeyq/UiBSYkptQ2mzH0wH7Iu7eWtSPSycty1PitVBVvsEVP6aOzlWAvTOoJcAL+SdKqcnAB8AQYINSqpybJ1cZ/26FTDMvcO86ENHdMWN5B8LdC7VI6NZLUXtWaPH70i44ZmxnkZqgKV9zaoU6Rjnl413bGMTX6id1+QxkpzluXB2dGoy9CqoT8KPZ+ePALBEJAt4DnrXZSsc2cX9anrfuC64OdEnLj4Q+4hsIsEqVcGgjvNsR/limzTyqO2mJ2szFHA9/8KmEWYxPHW2WZvGzErhwwjKUlY6Ojl3Yq6DqAOcBlFLtgAYUWOF9CVR+5rDqTPpFy/O6lfTxNe4CD35Z+GGdcQk2PQ0bnoQ8Q+XI4gjSLxROOaJcteW3ylrG9PQrHORXDJCsp0LR0Skt9iqoC0D+q3cPIF5E8uOVuJeiHx2wXN4D8KrEmFrBLWDkFs1U2pq9K6uvksrNtJEPS0GdJuDmWbmy+NSxjI8IkJ2iL/Xp6JQSexXLFmCqUuop4Dm0WVM+rdEs63TsJcNKQXlXctDHOk1g1E9w2xRw97G8tm81bByrOaVWFyRPS5VhbRRRp6lj952Kw78huFv5PiafMsVdfPzxx23GeysrDz/8MJMmTSryemJiIq1atcJex/R33nmHF154oeSKOjoOxF4F9R/gNPAamqn3NLNrQ4EdFSxXzSbjkuW5d+3Kl8HNA7o9C49t1/x5zNm7En6YXH32pFLOaukxzKnVCLzKH0yzzCnfldKces3JzTRGtIAPP/yQyZMnl1s+e3n99dcZMWIEXl5eGAwGOnfuzKxZs0zXDQYDnTp1Ys6cOQCMHj2aFStWkJBQfE4sHR1HYpeCEpHzItJLRPxFpIeImLvK9wSedox4NRRnLvFZE9wcHt5c2Hjif+/CL87zmLebrJTCiQW9aoFPyXmE7KXMKd89/DRZzMm8DJnFJ5rLza1Y/7SsrCyWL1/OsGHDAHB1dWXJkiW8/vrrHD6s+bLPmTMHpRTPPPMMAF5eXtx55518/PHHRfaro+Noym06JiI1P61jRePsJT5r6jSFhzbCkts1K7h8tkzVTOA7Vk7iu0OtS47HV5G0OVy+xHJPPPEEiYmJfP65FrJq4sSJ/P7772zZsoWtxlh8Y554nLfemoufjzevTnySoXf3gctnePixN2gUFsbEiROJiYlh2LBhjB07lrfffptevXrxySefsHnzZiZNmsTJkye55ppr+PDDD015g/bs2cOjjz7KsWPH6NOnjymLry3yc0iZJ1Ns164dzz77LI8++iiLFi1i1qxZ7NixA1fXAleE/FTpEyZMKNfnpKNTVopUUEqp0uSLFhF5tALkuTqoSjOofIKawbAvYFk/yDJ75/jqOWh8E9Rt6TzZqihz587l2muvZdmyZTRr1ozFixezd+9ey5TvFy4Sd/IEO79fR5/hY+nU4RpaNY8otCR57tw5Ll68SGxsLHl5efz555888sgjbNq0iU6dOrFixQoGDBjAkSNHUEoxcOBAxo8fz1NPPcWGDRt44IEHmDhxok05//rrL1q1Kmwp+tJLL/HFF1/QrVs3xo8fT/v2loYzbdq0Yd++fYXa6ehUFsUt8fVAy1xr76FjDzmZlj4xLm7gUUUCCTeMhAdWg6uZ1VteDnw7sfrsRzmAgQMHWmSBLXXKd/86RN12O31v68baTcZkg7lZFlmPXVxcmDZtGp6ennh7e7Nw4UIee+wxunTpgqurKw899BCenp7s3LmTnTt3kpOTw/jx43F3d+eee+7hhhtuKFL+5ORkm5HFPTw86NKlCxcuXGDo0KGFrvv7+3P5sh5TUMd5FDmDEpGISpTj6sHW7KkqhRqKuAUGzIf1owvKTvwEJ3dAk27Ok8uJfPnll/Ts2dPmtfyU7/nZWM2xSPkeEEJ4WCjx582WULNSTIq/bt26eHkVpEKJjY1l+fLlzJ8/31SWnZ1NfHw8SilCQ0MtlvWKS+9eu3Ztm8nttm/fzpdffsmDDz7I008/zTfffGNxPSUlhVq1ahVqp6NTWTgzo+7VSSELviqwvGdNh8Gw5xM4ub2gbOf7DldQJe0JpaSk4O/np/k7pVultAhqoTnJVjLvvfeefSnfXVw5lXCZdk3MoqgbsnE1aLNp6z2ksLAwXn75ZV5++eVCY27dupW4uDgtoZux3alTp2jWrJlNGTt06MDbb79tUZaZmcmjjz7KnDlzGDx4MO3atWPFihUmQwrQ8gdFRkZad6ejU2mU2sFWKVVPKdXY+nCEcDUSawOJqrD/ZI1SEP2iZdmRbyDpmO36lYWIFqHcWjl5BTpFOR09epRJkyaxYsUKPvnkE2bPns3evXst6kyZMoXs7Gy2b9/O5m++496777K47pFtFVXEyKhRo/jwww/ZtWsXIkJaWhpfffUVKSkp3HTTTbi5ufHOO++Qm5vLunXr+O2334qUs3PnziQnJxMXF2cqe+WVVwgPD+fhhx/Gx8eHBQsW8Mwzz5CYWDDD27p1K3feWQEpYHR0yohdCkop5aKUmqWUuoCWOPBfG4eOPVgv8VXFGRRA+M3anpQJ0QwmnLgX5ZabooUzMsfFHWqF2W5QQfTv39/CD+r//u//yM3NZdiwYUycOJHIyEhatGjBrFmzGD58OFlZWQA0aNCA2rVrExISwtChQ/nwww9pfX1Xy3syZEB2RqExO3XqxMKFC3nqqaeoXbs2zZs3Z9myZYC2d7Ru3TqWLVtG7dq1WbNmDXfffXeR8nt4ePDwww+zYsUKAH7//Xc++ugjFixYYKrTs2dP+vXrx/jx4wFthvX111/z0EOVY8Gpo2MLu1K+K6WeRcvH9AYwEy3vUh6ak24e8LqIlMbqr0rglJTve1fDl48XnLe7B+5ZXKFDxMTEEB0dXf6O/l4Hn4+wLAu7Ee5fWaEpJOxK+Z6bhSQcRmEWLcLFDeo005IHVjHyTcdtpnxPOmoZ9sjdB4JbOnQvMjExkW7durFnzx68vUvOOzZ//nxOnz7N7NmzyzReTU35Xp3vqyanfB8BTEdTUADrRWQK0AaIA/QlPnupLjMogGsGQsNrLctO74SlfSo3RYcIXIq1VE4oCGpeJZVTifiHWJ7npBfem6xg6taty+HDh+1STgBjx44ts3LSqZkopVyVUnuUUnbn/ysv9iqopsDvImIAcgFvABHJAeZRRBJBHRsUctJ1Qpgje3Fxgb5zC8eUSzqiJUGsLDIuQo5VoNVaoRWWgbjS8bQRYSLlbPWKf6hzNfI0UD7v9lJir4K6DOTbwMZjmV7DDS0dh449WL8pV0UjCXMadYKRP0CjzpblB9bB/s9st6lIJA9SzlmWeQZUTn6nchAdHW17eS+fgBDAbEnPkG0ZxUNHpwqhlGoE9AUWVea49iqoPcA1xr+/A6YppR5QSt2LFkD2zyJb2kApdb9S6pBSKk0pdUIp1c1YfptS6rBSKl0p9bNSKtysjadSaolS6opS6pxxX8y8zyLbVikyrKy2fKqBbq/fFh7aZGU0AXz5BPwT49ixUxO0h7cJpQWCrUq+Y2XBzavwPl7qeQvnXR2dSsJNKfW72THaRp15aEHDK3Wab6+Cmgfkx2aZApwDVgJr0PJB2RHeWUMp1QttL2sE4A90B/5RSgUD64DJaDOy34395zMVaAGEo0Wu+I9S6g5jnyW1rTpYW6FVYFBTh+LuBQPe1azm8snLgbUPwcV/HDNmblbh2ZNvUOXnd3IUfg3IM/8JiqHw/eroOJ5cEelkdiwwv6iU6gckiMgflS2YvdHMfxCRj4x/nwM6Ay2Ba4GWIrK/FGNOA6aLyE4RyROROBGJA+4GDojIZyKSiaaQIpVSrY3tHgRmiMglETkELAQeNl4rqW3VwTqbrnc1mEHl07ADDPwAi6WpzGT48smKNz8XMSYgLOg3T7mCX8Oi21Q3XN3I9rT6/tOTtHBYOjpVh67AAKXUSeBToIdSakVlDFymSBKi2aYfL207pZQr0AnYqJQ6jrav9SXwPNAWMEWmFJE0pdQJoK1S6jwQYn7d+PdA499FtgUOW8kwGhgN4ObmRkxMTGlvo1zceDEeL7PzXX8dJ+NExWZaTU1NdeB91aVxk2E0/feTgqJTv3Jw7QwS6ncvU4+1atUqFIrHLScV7yzLQPmZ7rUxpBf2GarO5Ln64aEu4yIFS3s5F2PJ9AkpplXVx2Aw2AyvVN2pzveVmZlZpueCiLwIvAiglIoGJojIsGIbVRB2KSilVIlPHhHZZkdX9dGWBO8BugE5wAZgEuAHWO8SX0ZbBvQzO7e+RgltreVcACwAzQ+qQvyFSsOvllGsu9x6R4Vb8lWYH1RRSBSsPAvHt5iKrjn7GdcMekFLhFhKDh06ZOlbInmQcNqykocfBo9a1dYHpShSUlJwCQyFSydNZe6GdNw9cF424AqgOvsLFUd1vi8vLy+uu+46Z4tRKuzdg4oBfi7hsIf819/5InLWmPjwLaAPkApYp0ANAFKM17C6nn+NEtpWHXKzIDu14Fy5gmc1DMapFNzxhuV+VPIp2LeqYvpPOQeGLMuymmAYURRegYVN+a/EXVUR5NevX09YWBh+fn7s2bOnwvrdvn27zVQjjuajjz4yReWoDM6fP0+bNm1MUUwchYjEiEg/hw5ihr0K6la09Bvmx73AcuAkYJfAInIJsNxYKOAAYDITU0r5As3Q9pYuoYVYMjcjizS2KbatPXJVGoX2n2prvkbVkeDmhRMZbn8L8gzl6zfjsmbNZo5PkNN8niIiIqhfvz5paQXLsIsWLarYGapSml+XOTkZDnfedSTbt28vlHqkOCZMmMC7775Lampqud7ylVIcP16w+9CtWzeOHDlS5v7KQnZ2NjNnzuT5558H4OTJkyilypUp+eDBgwwYMIBatbRVhB49erBz507T9fr163PrrbdahK+qTJRSP5Xi+NHefu01kthq41gnIo8AG4H+pbiXpcBYY9DZ2sB4YDOwHmh8rXboAAAgAElEQVSnlBqklPICXgH2i0j+HtLHwCSlVG2j8cMoYJnxWkltqwbV1YKvKG55FlzNlvSSY+HYD2Xvz5Ct9WGOixv4O9cwIjc3l//+97+OHcTDt7BPnAP9oio6rXx5iY2NpW3bts4Wo0LYsGEDrVu3JjQ0tOTKdnDixAm6du1K+/bt+ffff4mPj2fgwIH06tXLIkjw0KFD+eijjypkzDLggmY9Zc9h91t5RaTb+ArNsmOMnfVnAMHAUSATWAu8KiKZSqlBwLvACmAXcL9ZuynAB0As2lLhGyLyLYCIJJbQtmpQHX2giqNWqBZL0Hxpb/dCaHVH6fsS4b0ndxRx0TEP6ic/7GFXveeff57Zs2czZswYAgMLO1b/+uuvPP300xw9epSWLVvy3//+l5tvvhnQHHa7devGTz/9xP79+7nppptYtWoVwcE2HI0DQvj36AEefmYKf/51mBuvb0+rtteSnJLKihUrbMb3i4iIYNGiRfTs2ZPffvuNp59+mkOHDuHt7c2gQYN466238PDQXiKUUrz77rvMmzeP3Nxc/v33X55++mnWrVvH5cuXadGiBfPmzaNbNy2tytSpUzl48CBeXl6sX7+exo0bs3z5cjp16lTaj7rIz8Hf35+goCAMBgORkZE0aNCAEydOEB8fz9ixY9m2bRt+fn4888wzjBs3DtAMFd544w0WL15MQkICLVu25MsvvzQlXYyMjEQpxeLFi6lfv77FZ3bo0CGeeOIJ9u7dS2hoKK+99hoDBgwwyThs2DBGjhwJwLJly1i0aBE7duxARHj22WdZsWIF2dnZhIeHs2rVKtq1a1foXr/55huioqJK/RkVxdSpU7npppt49dVXTWXjxo3j0KFDTJw4kZ9/1nZYunTpwj///ENsbGyx+cEcgYhEO6LfilhfakUpnLdEJEdExohIoIg0EJFxRtNwRGSLiLQWEW8RiRaRk2btskTkEREJEJH6IvKWVb9Ftq0y1LQZFMANIy3Pj2+Bc3+Vvp+UsxUjjwPo1KkT0dHRzJkzp9C1ixcv0rdvX8aNG8eFCxd49tln6du3LxcuFHzXq1atYunSpSQkJJCdnW2zHwDcPBkydjId27ch6a+fmDx+JMtX2G/N6+rqyttvv01SUhL/+9//+PHHH3n//fct6nz55Zfs2rWLgwcPAnDDDTewd+9eLl68yJAhQ7j33nsxD6C8ceNG7r//fpKTkxkwYABPPWW3y2MhbH0Onp6epKZq+7L79u3jxIkT5OXl0b9/fyIjI4mLi+PHH39k3rx5fPfddwC89dZbrF69mq+//porV66wZMkSfHx82LZtm6mf1NRU7rvvPovxc3Jy6N+/P7179yYhIYH58+czdOhQu5YAv//+e7Zt28aff/5JcnIya9asISjI9u/3r7/+qtB9rx9++IF77723UPngwYPZvn276ftyc3OjefPm7Nu3r1Dd6oq96TYetHGMVErNA14HvimpDx1s70FVdxp1hFCrN+ptb5auj6yUwvtOVYzp06czf/58i3xJAF999RUtWrRg+PDhuLm58cADD9C6dWs2bdpkqjNixAhatmyJt7c3gwcPLpQzKp9Tp06xe89fzPjPGDw9Peh+Y0f69+xm975ex44dufHGG3FzcyMiIoLHHnuMrVu3WtR58cUXqVOnjilo7LBhwwgKCsLNzY3nnnuOrKwsiwf2LbfcQp8+fXB1dWX48OHlevjZ+zns3r2bxMREXnnlFTw8PGjatCmjRo3i008/BbQ9wJkzZ9KqVSuUUkRGRhapLMzZuXMnqampvPDCC3h4eNCjRw/69evH6tWrS2zr7u5OSkoKR48eRURo06YNDRvaXnpOTk6uUEu/pKQkm2M1bNgQg8HAxYsFzxV/f3+Sk5ML1a1slFKhSqm3jJEp/lVKtTOWj1dKdbG3H3tnUMtsHAuAx4AvgHH2DnhVY62gasIMCqC7VeDYgxu0wx5O/lI4gG4VpF27dvTr14/XX3/dojw+Pr7Qckp4eLhFcsAGDQqy6Pr4+JhmDI8//jgNGzbEz8+PWbNmER8fr6WJDyhYRgxv1FAzmLCDo0eP0q9fPxo0aEBAQAAvvfQSSUmWyR3DwixzZ82dO5c2bdpQq1YtAgMDuXz5skUba9kzMzPLvH9V1OdgTWxsLPHx8QQGBpqOWbNmcf689hJz+vTpIrMHF0d8fDxhYWG4mBkmWX9XRdGjRw+eeuopJkyYQP369Rk9ejRXrlyxWbd27dp2+0pt377dlGesqD244OBgzp4tvMJw9uxZlFIWy8UpKSk2l6ErE6VUW+AvYDha7NbGQP5mdTha0Fm7sHcPqomNskwRqdqvvVWNmrYHlU/LO6BBe8ulvQ1jtdh9tSOKbnd2H3z6AEQtBODJV2qDcjGm0SicnsbZPijTpk3j+uuv57nnnjOVhYSEEBtradhx6tQp7rij5H24Dz/8kDfffNN0T7GxsVqaeJdAfNF+WqfizmnW9YZsfH19SU8v8KMzGAwWM7onnniC6667jtWrV+Pv78+8efP4/PPPLcY0Ty2/fft23njjDX788Ufatm2Li4sLtWvXxp4ccY4kLCyMJk2acOyY7QzOYWFhnDhxwub+T3GEhIRw+vRp8vLyTErq1KlTtGzZEqDQ53vunGXYqXHjxjFixAgyMjIYPHgwb775JjNmzCg0TocOHTh69KhdMnXr1q1IRZ1Pz549+eyzzxgxwjI329q1a7nxxhtNe4y5ubkcP36cyMhIW91UJnPRop7fjmZnYB5M81cK0jaViL1WfLE2Dl05lRbrPajqFOaoOJSCAfMt/aKyLsOKeyDe9jIOKedg5WDIvGxZXjvCpnKqCjRv3pz77ruPd955x1TWp08fjh49yqpVq8jNzWXNmjUcPHiQfv1K7yoSHh5Op06dmPLaXLLzXNjx2x42/WD0f0+/SMuWLcnMzOSrr74iJyeHmTNnWvi9pKSkEBAQgJ+fH4cPH+aDDz4odryUlBTc3NyoW7cuubm5TJ8+vchZQWXSuXNnAgICeOONN8jIyMBgMPD333+ze/duAEaOHMnkyZM5duwYIsL+/ftNe37169fnn39sx4bs0qULvr6+zJ49m5ycHGJiYti0aRP336/ZU1177bWsW7eO9PR0jh8/zuLFBYlEd+/eza5du8jJycHX1xcvLy9cXV1tjtOnT59CS6vlYcqUKfz666+8/PLLXLx4kZSUFObPn8/SpUuZNm2aqd5vv/1GREREpRtI2OAWtCS2qRR2KToPNCjcxDbV1AmnmpJmudyCb13nyOEIQq6DXtMtyy4cgwVRML8j7PygIN9R5hX4dAikWgVG9W9QOE9SFeOVV16x8IkKCgpi8+bNzJ07l6CgIGbPns3mzZttW+nZwapVq9j122/UaX0z095awIP39NUupF+gVkAA77//PiNHjiQ0NBRfX18LX6M5c+aYLONGjRpVyEjAmttvv50777yTli1bEh4ejpeXV6ElwOK48847mTVrVpnuszhcXV3ZtGkTe/fupUmTJgQHBzNy5EguX9ZeZp599lkGDx5M7969CQgI4NFHHyUjQ1sGnTp1Kg899BCBgYGsXbvWol8PDw82btzIN998Q3BwMGPGjOHjjz+mdWstZOczzzyDh4cH9evX56GHHjJZBQJcuXKFUaNG0bhxY8LDwwkKCmLCBNs50fr378/hw4eJj4+vkM+jRYsW7Nixg3379hEREUFgYCCTJ09m/fr19OrVy1Rv5cqVPP7448X0VGkUZzQXTEHAhhIpMuW7UioP2w61NhER268TVZhKT/n+YTc4ZxZXd+RPmpFBBePwUEdFIaIpniNf274e2klbDvz7c0gscFE7dPta2lzTFgJCi40W4ewlPkdQ5D3l5cK5v5k69wOOnzzNivmvauntvawDplRNauJ3Bfbf14IFCzh48CDz5s2rcBnOnDnDjTfeyLRp03j00UcBSEhIICoqij179uDl5WWzXWWlfFdKbQGuiMjdxvirOUAnEflTKfUp4CMiA+zpq7g9qOkUKCiFljXXG9hEwTStH5o2XGyrAx0rCs2gaoiRRD5Kwd0LYcMY20YScb9rhzVuXiUqp6sOF7fCVp5pSdVGQV3tjB5tK6VSxdCoUSO++eYbNmzYQGpqKn5+ftSrV49Dhyo12W1xzAC2KKW+B1ah6ZGeSqmngf9DS7FkF0UqKBGZmv+3UmoSmoPs7SKSblbui5bAsGq5pVdFRLRUCubUpCW+fDz94N7lcPQ72PUh/FNCmMZGN2iJ+3TlVBjrhIZZl7VoG66lD8irU7No37497du3d7YYNhGRrUqpgWh5BJcYi19HC4s3UER22duXvVZ8jwFPmisnoyBpSqk5wHzgVZstdTSyrlhmhnX3qbLGAOVGKS2aRKs7IOk4fPEonLVhLNG8J9yzFP4t2cz3qsTdh6kTx0Ou2TJ0+gWnh37S0SkJEfkK+Eop1RyoB1wQkVIHRbRXQQVTYMdujQdQw9aqHID18p5P2TbRqx3BzWHUzxC/B/75SbPq8/DVlFPbu8HVDdAVlE2U0v6fXCkIbUTaBfCtBy7VbstX5ypBKdUhP4mtiBynDLkD87FXQf0OTFNK/c+Y/TZfkFC07LW7yyrAVUOh/aerREGBFrG9UcdiDUJExMJHR8eIT21IiddyZAHk5WhBZP3tttTV0als37a9Sqm/0AJ8rxKRMscxs9fMfBxaRtsTSqkYpdQapVQMcALNWMJuz+CrFuvI1DVx/6mMeHl5ceHCBac7iFZJXNwKz7ZTz5c/rYnOVYOIcOHChSKt+xzAA2g2C68Bp5RS3ymlhiqlfErbkV0zKBHZY1xLfBa4EWiPlp9pDvC2iFworr0ONgwkgklLg7SKzfYOwKVL7iQkVFx/AQHgyP/bjRo14syZM4Xi3FmTmZlZmT+ySsGue5I8uJJUMIsCOJ9VpTPu1sTvCqrvfXl5eZUqP1d5EJE1wBqlVDCashoGfAKkKqXWA5+IyJbi+sjH7nQbRiX0chnk1QEbM6hgMjPBEXEdU1PdKrRfLy/HKih3d3eaNLEVTcuSmJiYapeyuiTsvqeYTRDzWsF57QgY+2eV3Yuqid8V1Nz7cgTGjOnzgflKqRZosflGAkOxU/fokSQqi6vVSEKnYuj0iKV5+aWTcPRbp4mjo2MvxqW9zsajHqVwSypSiymlfgLGiMhh49/FISJym72DXpVU4zBH+tZQFcCvHrS/F/auLCjb+QG07us8mXR0ikBpFk890WZNAwE/4H/Ak8Aae/spbgZlblJVUjpffSZWEjaMJPQHv06p6GIVZ+3kdi0ivI5OFUIp9SZwBi2Iw01o0c2bi0hXEflIROzegChSsYjIrSJy2Ph3tPG8yKMUwscopTKVUqnG44jZtSFKqVilVJpS6kulVB2za3WUUuuN12KVUkOs+i2ybZXAOpJ5TQtzpON4GnaAiG6WZb/81zmy6OgUzSPARuAWEWkhItNExHaI+RKwN6NuRT9NnxIRP+PRyjhGW+AjtClhfSAdMM9X/R5aXpH6aJtsHxjb2NPW+ehm5joVwc1jLc8PrNf2o3R0qg4NReQJEfm1vB3ZuzR31jgrGaSUclQgsKHAJhHZZswjMhm4Wynlb4z5NwiYLCKpIrIDTUMPL6mtg2QtHXl51dpIQl+KrEI07wV1zSJSSx78+q7z5NHRsUJEskGLKKGUekopNUUp1cBY1rw0z2V7FdRkoCnwGXBOKfWhUqpraQU34zWlVJJS6helVLSxrC1gWlAXkRNoM6aWxsMgIuZpKvcZ25TU1vlkJoOYOVZ6+IO7l/7g1yk9Li7Q1covfs+Kwi9AOjpOQinlqZT6DNgDvAO8ghboAWA2pXBXstdR9w3gDaXUtWizlvuBUUqpk2gOWCtFxHZ+5sJMBA6iKZD7gU3Gfv0Aq/SqXAb8AUMx1yihrQVKqdHAaAA3NzdiYmLsFLvs+KSdobPZeYaLL7tiYkhOdiclxW5XNLtJT09n1y67AwaXSO3aOfj5OT9gfWpqaqV8X5VJWe5J5QXTxTMIryzjvmZuBifXvszJJg9UvIBlpCZ+V1Bz76uCeZUCC74f0NIz5fMNMAZ4wZ6OSvV0FJG9aHGWnjcKMAx4Dphkb19WodaXK6UeAPoAqYB1spsAIAUtQ2NR1yihrfX4C4AFoCUsrJTEfrG/WkQr9A5uTHR0NImJcOlSxQ+3a9cuunTpUmH91asHgYEV1l2ZcVoiRgdS5nvyfg6+e8l0GpHwHRFD36oyEfJr4ncFNfe+KpgHgEkissqYsNCcf4EIezsqk3m4iOQBaWjJCnOwNEkvdXfG9geAyPxCpVRTwBM4ajzcjN7I+UQa21BCW+djI4pEdUJfiqyCXP8geNUqOM+4pC316eg4nyCgqOyJLmjPZrsolYJSSrVQSk1XSp0AtgF3olnP2ZU5SykVqJS6XSnlpZRyU0oNRcuu+B2wEuivlOpmNIqYDqwTkRQRSQPWAdOVUr7G/a+70JYXKa5tae7PYRShoPQHv06Z8fSHG0ZZlv36LhhynCOPjk4B/6L5P9miM2B3Xih7zcyfVErtBA4D49GUUy8gXEReFJGDdo7nDswEEoEkYCxahsUjInIAeBxN2SSg7R+NMWs7Bi3lfAKwGnjC2AY72jqXNCsfqGpkwadThenyGLiavYxePgUHNzhPHh0djY+BF4wTkHyrb1FK3Qo8Q0GW3RKxdw9qHrAFbdNrvYhklEJYEyKSCNxQzPVVaDnsbV27iBYyo9RtnY7uA6XjCPzqwbVD4I+lBWW/LYT29zhPJh0dzVIvEm2Fa5GxbAfgBXwqIvPt7cheBRUmIudKJaJOAUUoqOqyxFdd5LwquelJSwV1eiec+xsatHOeTDpXNSJiAO5XSr0H3I4x5TvwrYhsLU1f9lre6cqpPOhhjnQcRXALaBIF/5r97n9fDP3edp5MOjqAiGwHtpenDz2aeWWgL/HpOJIbHrVUUPvWQM9p4GXteaGjU73Qo5lXBtU41YZONaBVH/BvWHCekwb77c5ooKNTbpRSeUopg51H+fNBmUcoF5Hocsp/9ZJnKLzE56Mt8VWXvZ3qIudVi6s7dHzYMuPu7kVww0hQ5XFR1NGxm+loPq0VSsXH2dGxJP0iFt+bV6D2QNHRqUiufwi2zi6I+Zh4GGJ/gYhbnCuXTrVHKeWF5lrkiaYzPheRKeZ1RGSqI8Yubg+qe2k6EpFt5RenBpJuvbyn+0DpOICAhlp23UMbC8p2L9IVlE5FkAX0EJFUpZQ7sEMp9Y2I7HT0wMXNoGIoePVXlDx9s465pAM1wkBCX+KrJtww0lJBHdoEKefAv4HzZNKp9oiIoMU7BS3YgjsOWM6zRXHGDbcCPYzH/wFxaCGJRqAFdx0BfG8sL9KB9qqnkIFEwQxKf/DrVChNukOQWbjKvFz4Y7nz5NGpLrgppX43O0ZbV1BKuSql9qJF6vnBKui34wQr6oK5Q5VSahnwvYiMtKr2sVJqMXA3sMkhElZ3qnGiQp1qhlLaLOrbiQVluxdp+aPcvZwnl05VJ1dEOhVXweh8e61SKhBYr5RqJyJ/O1owe83D7wKKsltdY7yuY4sasMSnU42IvB88/ArO0xJgX9WMAKZT/RCRZLTtnzsqYzx7FZQL0LyIay3Q95+KppCRRPVTUPpSZDXCO1AzOTdn2xzITneKODrVH6VUXePMCaWUN1ouwMOVMba9ZuZfYUzTjpbGwmBMRDUILTr5ZkcJWO0pNIMqCHOkP/h1HMKNY2DXR5BnTL1xJQ5+mQe3vlR8Ox0d2zRESy7rijZZWSsixT7zlVIBaLYKjdGCxJojIjLDnoHtVVDjgDC05bxcpdQloLax/Q7jdR1bWKfaqIYzKJ1qRq1QLRXH/94tKPv1XejyOPjUcZ5cOtUSEdkPXGdvfWO+vk1AUXm4BbBLQdm1xCciSSLSDS0y7QzgC+O/vUWku4hcKLaDq5m0BMtz3UhCpzKI+o/ly1BOGuz8wHny6FxNzANOoqVW8hIRF6vD7i2hUkWSEJEfgB9KJerVjAhcOWtZZuaTUl2W+KqLnDpmeNWCm8fBD5MLyv73npYqPjDMeXLpXA20AQaLyB/l7ajUQV6VUvWUUo2tj/IKUiPJuqK9uebj5gXetZ0nj87VRadHLP+/5aTB1887Tx6dq4VTaGGRyo29Kd8DlFJLlVLpwFm0nPPWh441hWZPDfXgnTqVh6cf9JhkWXb0Gzj5i3Pk0blamIaW8r3c+V7snUG9B9wHLAaeAB6xcZQKpVQLpVSmUmqFWdkQpVSsUipNKfWlUqqO2bU6Sqn1xmuxSqkhVv0V2dZpXImzPA8IdY4c5URf4qvGdHwEwrpYlv38qv6l6jiSfkB94F+l1Gal1MdWh93hTezdg7odeF5E3iuLtEXwHrA7/0Qp1Rb4COgL/AksAN4H7jern41249cCXyml9onIATvaOocUqxlUQEOLU/0ZoeNIcnJyOHPmDJmdZ8M15y0v/rXXodElatWqxaFDhxzWv7Oorvfl5eVFo0aNcHevlEwKt6BZ6l0B2tq4bveTz14FpYAj9nZaYmdK3Q8kA79S4AA8FNiUHxVdKTUZOKSU8gfy0Hyu2olIKlo03Y3AcOCF4tqKSEpFyV1qbC3x6ehUEmfOnMHf35+IiAjUxROQZfZTcPeB4JYOW3JOSUnB39/fIX07k+p4XyLChQsXOHPmDE2aNKmM8SpsEHsV1KdAf2BLeQc0rktOB24DHjW71BZNYQEgIieUUtlASzQFZRCRo2b19wFRdrS1sCQxBkIcDeDm5kZMTEx5b6lIWh7ZTYjZ+bGEdOLMxktM9CQzs+KTEaenp7NrV8XFcvT1zaVOnZwK66+spKamOvT7cgaOvKdatWoRFBREamoqLq618MVMQeWkk558HoObr0PGNhgMpKQ4793QUVTX+/Lw8CA5Obna/X7sVVDfA/OMs5mvgYvWFUTkJzv7mgEsFpHTyvLtzQ+4bFX3MuAPGIq5VlJbazkXoC0B4uXlJdHR0XaKXQbiP9BMSoy0uL47La4pGO/MGUh3QASaXbt20aVLl5Ir2klAADSoAhkbYmJicOj35QQceU+HDh0iICB/n9ofDCmQVfAz8clJhsD6DplFVceZhj1U5/vy8vLiuuvs9rctM0qpncBPaDH7totIRln7sldBbTD+2wR42KxcKMgVVaLzlVLqWrQ4TrY+pVTA2uojAEhBm0EVda2kts7jSrzluX+Ixam+B6VTqQQ0gESz97jcDMi8rMXv09GpOE4AD6Ftv2QrpXYDP6Mprf+JSJa9Hdm7vnRrEUcPs3/tIRqIAE4ppc4BE4BBSqk/gQNAZH5FpVRTNFv6o8bDTSllluyGSGMbSmjrPKwVVECI7Xo6OpWBuw94WSmjlLNOeVOKiIhgy5by7RhER0ezaNEiu+pOnTqVYcOGlWu88lAZ48fExNCoUSOHjmEPIjJUREKBa4Bn0daRHgd+BC4ppX60ty+7ZlDmuaHKyQK0/ax8JqAprCeAesD/lFLd0CzxpqMFpk0BUEqtA6YrpUaiWfHdBdxs7GdlcW2dQm6WZSRz5QJ+9W1WzUjNJvbvCwQEexHSXHfk1XEcR8+GwEXzFReBsxU3i2rZskK60akBiMhh4LBS6mO0ycnTaLYH0fb2UapQR+VFRNIB066LUioVyBSRRCBRKfU4mrIJQjPIGGHWfAywBC2j4wXgCRE5YOz3QAltK5+Uc5bnvvXAtfDHbcjJ49uP/iLtcjYANw1sRrPr61WGhHajL0XWINw9wSsAMq8UlKUlaaGRdCdynQpAKeWFZmqev7rWEUhDCyw+AW2pzy6KXOKz4VxV3FGmvNIiMlVEhpmdrxKRxiLiKyJ3ichFs2sXRWSg8VpjEVll1VeRbZ2CHct7IhB/LNmknAD++C7W0ZLpXO1YR9Q3ZGt7UU4iKyuL8ePHExISQkhICOPHjycrq2CbYsOGDVx77bUEBATQrFkzvv3220J9nD17lg4dOjBnzhwA/v33X6KiovD396dXr14kJVnmZdu4cSNt27YlMDCQ6Ohok2/T0qVL6d+/v6le8+bNGTx4sOm8TZs27N27FwClFB9++CEtWrSgdu3aPPnkk4idb3M7d+7k5ptvJjAwkMjISJN13aeffkqnTpbJbd9++20GDBhg+qwmTJhA48aNqV+/Po8//jgZGWW2QXAUl4CNQCc0+4VuQJCI9BeRt0Rkr70dFbcH1d3Ysb2Hjjkp9u0/JZy2XIXMzshF8vQpi44DcfMsvKSXcck5sgCvvvoqO3fuZO/evezbt4/ffvuNmTNnAvDbb7/x4IMP8uabb5KcnMy2bduIiIiwaH/y5EmioqJ46qmnmDBhAgBDhgyhY8eOJCUlMXnyZJYvL3iHPnr0KA888ADz5s0jMTGRPn360L9/f7Kzs4mKimL79u3k5eVx9uxZcnJy+OUXLTTUP//8Q1paGh06dDD1tXnzZnbv3s2+fftYu3Yt3333XYn3GxcXR9++fZk0aRIXL15kzpw5DBo0iMTERAYMGMCRI0c4duyYqf6qVasYMkQLnDNx4kSOHj3K3r17OX78OHFxcUyfPr1sH7zjSEWzAaiPtnVTDyiTP0ORCkpEIkSkiZ1H0zLdRk3GTiddMeQVKktNznSERDo6BXgHWZ7nZEC2c97EV65cySuvvEK9evWoW7cuU6ZM4ZNPPgFg8eLFPPLII/Tq1QsXFxdCQ0Np3bq1qe3BgweJjo5m2rRpjB49GoBTp06xe/duZsyYgaenJ927d7eYFa1Zs4a+ffvSq1cv3N3dmTBhAhkZGfz66680bdoUf39/9u7dy9atW7n99tsJDQ3l8OHDbN26lZtuugkXl4LH5gsvvEBgYCCNGzfm1ltvNc2uimPFihX06dOHPn364OLiQq9evejUqRNff/01Pj4+3HXXXaxevRqAY8eOcfjwYQYMGICIsMCgG0kAACAASURBVHDhQt5++23q1KmDv78/L730Ep9++mkJI1YuIlIXzVJ7KdAM+Bi4qJTapZR6XSnV296+Kt5LVEfDTgu+1EuFLS4vna1a6bn1PagaiLsneFi91KYn2a7rYOLj4wkPDzedh4eHEx+v/X5Onz5Ns2bNimy7cuVKQkNDueeeeyz6q127Nr6+Bfdn3r/1eC4uLoSFhREXp8XOjIqKIiYmhm3bthEVFUV0dDRbt25l69at3HLLLRbjNzBzEPTx8SE1NbXE+42NjeWzzz4jMDDQdOzYsYOzZ7WX2iFDhpgU1KpVqxg4cCA+Pj4kJiaSnp5Ox44dTe3uuOMOEhMTixvOKYjIfhH5r4jchWYXcBta6KP/AN/Y24+uoByFHUt8IkUoqHNVS0Hp1FCsU79kpUB2mu26DiQkJITY2IK911OnThESov1ewsLCOHHiRJFtp06dSnBwMEOGDMFgMADQsGFDLl26RFpawb2cOnWqyPFEhNOnTxMaqgVzzldQ27dvJyoqiqioKJOC6tq1a7nvNywsjOHDh5OcnGw60tLSeOGFFwDo3bs3SUlJ7N27l9WrV5uW94KDg/H29ubAgQOmdpcvX7ZLKVY2Sil3pVQ3pdQUNKOI79CUVALwmb396ArKUdixxCcipF4qvJx36VzlPyR0rkI8/QsHjE09X+lT5gceeICZM2eSmJhIUlIS06dPN/kMPfrooyxdupQff/yRvLw84uLiOHz4sKmtu7s7n332GWlpaQwfPpy8vDzCw8Pp1KkTU6ZMITs7mx07drBp0yZTm8GDB/PVV1/x448/kpOTw9y5c/H09OTmmzWvlaioKH7++WcyMjJo1KgR3bp149tvv+XChQtERkZSXoYNG8amTZv47rvvMBgMZGZmEhMTw5kzZwAtBNs999zD888/z8WLF+nVqxegzfRGjRrFM888Q0KClqk7Li7Orn2vykQp9QOaocRWYCya1fXzQHsRaSAidgfx1hWUoyi0xFc41UZGSja5OYX3oC4nVTmrHJ2aiFKFffNyMivdom/SpEl06tSJDh060L59e66//nomTdLyWHXu3JmlS5fyzDPPUKtWLaKioixmP6DFmVu3bh0JCQk88sgj5OXlsWrVKnbt2kWdOnWYNm0aDz74oKl+q1atWLFiBWPHjiU4OJhNmzaxadMmPDw8AGjZsiV+fn5066bZfgUEBNC0aVO6du2Kq6vd2cqLJCwsjA0bNjBr1izq1q1LWFgYb775Jnl5Bc+CIUOGsGXLFu69917c3ArcU9544w2aN2/OjTfeSEBAAD179uTIkQqL411RZAKTgeuBuiIySETezXcLKg3KXrPImoiXl5dkZjrAICEvD2bWgzyzAKsvxRda89+97Qqb3i/8nbm4Ku6f3AUXl7L5pVR0LD5fXwitAqms9Fh8pePQoUO0adOm5IoX/7FUSsoV6rbSrP3KSHWOWVcc1fm+ivr/oJRKFxHHRA0uJ8X5Qa1TSjU3/v2gUiqoqLo6VqRfsFROnrUKb0gDlxNtK8c8g5B+2e5wVTo65cM/BC2kphExFE62qaNTSpRS/ZRSbyqlFiulZiul+pS2j+KW+O4C8rPS5psL6thDoUy6ti34LicUbQyRclE3NdepJNy9Cv8fzbysLffp6JQSpZS/UmormrPu0/D/7Z13fFzFtfi/Z3clrXqvli13XDHGODbd9BJaQgkmoSVAyo+8hIS8F5IQSHtppJD3AgkJoRMIjxI6mGJjwBhcMLZxL7JVLFm970q78/tjVtLe1a60kiXtrjTfz+d+dO+0O6Mr3XPnzJlzOB/4NvCCiKwUkZRw2+pPQFUBx3ffk0FEQRz3DBBJt5uGfgRUMOs+g2HESM6FuERrWktV8LIGQ//8N3r96WogUSlVCCQC1/jS/zvchvoTUP8C/iAiHrRw+kBEPCGOriEPZSwSOINKDT6DaqgObQwRTTOocbxMOX4IZjDRXhcRs3NDzHMp8COl1GNKKQ+AUsqjlHoMbTxxabgN9ecs9hbgPbTL9DuABwGjmA6HQBPzICo+d0eXb5YU3BCiJYoElCF2UUoh4TqBdWZow4guv9l7w0EdGt5mDH5jmVE2hssGPg2R96kvPyxCCiilR/QUgIhcB9ytlNoUfh/HMWGo+Bqq+t+M22xUfIYjxOl0UltbS3Z2dnhCSgTSiqHOb2NsV7vedJ4e+ThDhqGhlKK2than0zlw4eFhH3ABsCJI3vm+/LAINx7UlHAbNBDESKKvjXZdhVV1klOcQk1Z747w1rqOwX39GgwBFBcXU1ZWNnhXOG2tAaq9w5BaB/a4sJvo6OgYzRfiqBGr43I6naMZzPCvwO98xhCPoQMWFgBXAjeggxiGRdjxoESkEPgucCrauq8WHXP+90qpQ/1UHX+E4UUiUEDlT0mnoaqtZ+Ou2+XB1daJMzl+xLoZLmYNKjaJi4tjypQhfFt2NMK9J0Fjr3sg5l0Kl/0j7CZWrlzJwoULB3/vKGesjms4UUr9QURy0ctE1/mSBXABv1JK3R1uW2EplkVkJrAJ+A+0K/UP0QGovgV8HBCK3dBHxdd3DaouwJ1Rel4iKZnWL7PmOqPmM0QAZzqceYc1bcvTsPvIQrQbxg9KqR8ARWhV3zXAZ4EipdQPB9NOuCufvwYagZlKqdOUUsuVUqcBM33pvx7MTcc0rmZw+UUrtcdDUt81wcAZVEZeIimZ1p37wfz0GQyjwtzPQ/58a9oLt5i9UYawUUrVK6Ve8VnzvaKUGnTQsXBVfKcBX1NK7Q/oQKmI3AncM9gbj1kay6zXaUV9Qml3ujw013YAesYkAmnZiaRkBc6gouNlYFR84xCbDS78I9x/Fiifj7jGA/DR3+CEb0a2b4aoQ0ROGUx5pdQ74ZQLdwYVDzSHyGv25YeFiDwqIpUi0iQiO0XkBr+8M0Rku4i0icjbIlLil5cgIv/w1TskIt8JaDdk3VGl4YD1On1inyLd3sqV0oIrJdOJI95OaoCAajUqPkMkKT4OPnOTNW3176C1NjL9MUQzK4G3fYf/eagjLMIVUB8D3xQRS3nRJmbf8OWHyy+ByUqpNOAi4OciskhEcoBn0Bu5soB1wJN+9e4EZgAl6Bndf4rIub5+DFR39AgUUBmT+hSpq+yr3gNIybKq+JqNis8QaU75HsT7OUdtr4fXB7WMYBgfnAac7js+h94z+xpwPdq0/HrgdV/6JeE2Gq6K76fAi8A2EXmSXrPBy9FC47Ph3jDA5bryHdOARcBWpVT33qs7gRoRmaWU2o5eaLvep8esF5G/oS1EXgU+P0Dd0aPxoPU6yAwqcP0pLTcJgNQ+RhJGQBkiTHIOnHwLvPnT3rRN/4Rjr4GSEyLXL0NUoZRa1X0uIg8Cryulbggo9rCI3I9+X79AGIS7D+pVEbkA+DnwQ3p9860HLlBKvR5OO92IyD1o4ZIIbAReBn6BthTsvmeriOwB5opIFdoixH+j8CZ6JfHcUHUBi4ASkZuAm0AHBlu5cuVguj4gc3auJ8/vevuhVg4F3KN0izUGVE1zBWvXVuL1Whd72ps7WbPmg0GH3Whra2Pt2rWDqtMfcXFe9u6NvLqxpaVl2J9XpImFMYn3aBYll5DS2huHqfmpm1m/6C6Q4EqYWBjXUBir4+oPEZkIPIyelHiB+wYwFb8YuCJE3pPAE+HeO+x9UEqpV4FXRSQJyATqlVJDik2ulPqGiHwT7Yx2Gdo+PgUI3FHYCKT68rqvA/MYoG7gve8D7gMdD2rYY/Hs/rnlctaSs5g19VRL2iNvvI+O6aVZtHQeWUV6iIfeX09ro7snb870o0n3zbDCZbjjQSUkQElkVvQsmHhQEWTaX+GBc3suU1v2sCyzCo5ZHrR4zIxrkIzVcQ1AF/BdpdQGEUkF1ovICqVUKHdGNmA6wT1JzADCjvo4aAdbSqk2pVT5UIWTXzsepdS7QDHwdfT+qrSAYmloI4wWv+vAPAaoO7o0BKj4AtagOl0emmp6hZMIpOX0epFOzrCuQ7U2RH7mYjBQcjzMudia9uZPjDPZcYBSqlIptcF33gxsA/oLYfoS8EsRuVxE7AAiYheRK9BauBfDvXc0eIB0oNegtgILuhNFJLk73bfuVOmf7zvvXs8KWXdEex5Ilwta/J1qSB83R/UBG3S7Lfh6rjOs61AtDW4MhqjgzJ/ofX3dNFfCe3+KXH8Mw4VDRNb5HTeFKigik4GFQH9rCP8BbEar89p9SzTtaNXeZl9+WIyqgBKRPBG5UkRSfBL1HGA58BbwLDBPRC4VESfwY+ATPyOHh4EfiUimiMwCbkR7WCeMuqND4B6o1EJwWC3wAy340nOtMXiicQZl9kEZAMiaAku/bk17725oNEEOYpwupdRxfsd9wQr5fOs9DXxbKdUUrAyAUqpGKXUycA7wM1+dnwFnK6VOUUqFvU9htGdQCq3OKwPqgbvQg/23UuowOk7IL3x5S9DOBbu5A9gDlAKrgN/61sUIo+7oEGjBlxFkD5SfgFJAep51fSkl0yrQjDcJQ1Rx8q2QlNN73dVutfAzjElEJA4taB5TSj0TTh2l1Aql1M+UUt/w/Ry0r6ywjSSGA58gObWf/DeAWSHyXMCXfceg6o4agetPYZiYp+cFzqACNutGwQzKYOjBmQan/xBevKU37ZMnYMlNMGFR5PplGDF8+13vB7YppX4/yLo24A3gq0qpXYO994AzKBGJF5ENInL2YBsfd4Qxg+qzSTdAxdfHH58RUIZoY+E1kDfHmvbaD40ueOxyIjp8++ki8rHvOD/MuoK21O5jUR0OAwoopZQbmII2NTT0xwAzqE6Xh6Zaq8rO34IPICktHvHb99TR0kmn2zO8/Rwk5r1jsGB3wDm/sKYdWANbw9L8GGIMpdS7SilRSh2tlDrGd7w8GvcOdw1qBWBmUAPRZwZlNTGvP9SqF558JAdY8AHY7DaSU63rUEbNZ4g6pp0OMwJeCa/epl0hGQzDRLgC6n+A5SJyl4icJCLTRGSq/zGSnYwZBnAUWz+Aeq+b5Mzos+QzGPpw9i/A5hdlt6UKVtwRurxhPOIFfgJUDKVyuAJqFTARHap3FbAT2BVwjG+8nr6h3gPWoPqYmOcF9xCREoWm5gZDH3Jnwkm3WNM2PAT7VkemP4aoQER6YssrzU/8o66LSEhDuUDCteK7fhD9G580V4LXb5kuMQviky1FDh+wOrYItODrJnAG1VIfWQFl1qAMITn5u7D1Waj1+0Z97hs45v0ycn0yRJrXROREpVRDYIYvbtSLhGk0Ea6z2IcG179xSJ8wG9bZk/Iqqvb17m1TCrKLrAKsm8DNusaSzxC1xDnhoj/BA+f1pjUeYNb2P8Lp52qDCsN4owV4SUTOUEr1WIWJyEloN0jPh9vQoDbqiohNROaJyKk+d0KGbur2Wa8zJ1su66vacHf0WuMlJDlIyw4+g+qr4jObdQ1RTMkJsMTqYSKn9iN4/UcR6pAhwnwWyAae6o4hKCInoKNWvAR8KdyGwhZQIvL/gEPAJ2jXREf50p8TkbB9K41Z6vZar7OsdiNV+xot1/lT0izm5P4kB8SFaqk3/vgMUc5ZP4GCo61pa++FLcb0fLyhlKoBzgWOBe4XkeOBV9ABDL+oVPiLBmEJKBG5EbgbeA4d58P/zboa7WZofFMfMIMKEFCH9lldV+VPDnS+3ktSapwlBpSrrRO3K3Lb0MwalGFAHAlw1ZOQVmxNf+k7xvR8HKKU2g+ch47ZtwrtTeJKpdSgNnWGqyD+DvA7pdR/dbtP92M78L3B3HRMMuAMyiqg8qak0xmiKZvdRlJ6vMU4oq3BTXy+0ecbopi0IvjCI/CPc8Djm/W318M7d/Xd2GsYU4hIUBd06PWm89Dh3q/VXpNAKfWPcNoN9403BT09C0YrkBFmO2OXQAGVOaXn1N3RRV15iyU7f3IaZZWhm0tOT7AIqNYmFxn5gwtcaDAMFuXx4G3vQHW0423Xh+rosKZ1dEBXF6qzE9Xp+9nVherqRHV2wuHTUPvXoJSgPKA+fBz1Uh3K5uwt09ml63TX7exEdfnSvd4QnQsxlR9kuiJU+VC/lN6MHJeLXQkJ/d93lPt0ROnDx98HyL/X71wBwyqgaoDJIfKOAsa3v/22OujwW2NyOHWoDR+HS5stfx8Z+Uk4k/02OAYhOd1qKNHWaNahDOGhvF68zc14Ghv10dCIp0mfe3uum3rzGxt8eU0o93D9naVYL/euHKZ2I4sd4/MtBFMGLjJ4whVQLwA/FpGV6HAXAEpEcoBb0GtT45c+FnxTwNa7vFe1P2D9aUro9aduoikulFmDii6UUnjq6+ksr6Czwu8oL8e9bx/ugwehy7xGDaOHUqp04FKDJ1wB9SPgdGALOpKiAv6EDm9RDYzvgDA1O6zXgQYSe60WfAVT0gZ86SelB/jjazR7ocYj3o4OXNu30751Kx1bP6Xj009x79+P6jBbDwzRiYgsBSYppf4VJO9y4IBSqr+IvD2Eu1G3VkSOA76NjpK4x1f3f4E/9BddcVxQ/an1Os8alqrmoHX9Ka8fC75ukgMElFHxjW28Lheubdtw7d2Ha9cuOjZvxn3gAF3V1aPbEREkMRGb04ktMRFJdGJz6mtJSsTmTEScCUhcHOKI8/106J9xDvCd7y0tZXpxFvLBnxBvB2IDsSlIykQ++yskJaenTncbOBy6TUegHZa1f4NKJ3h66OL9t79mzRqOP/74MBoafJ9CJg92zINNz84Onj50fgm8EyJvNjpo7enhNBS2WZhSqhkdtvdn4dYZN1Rvs177xcpxtXXSXNf7tSs2IbsoQD8fhKR0o+Iby3RWVdG6Zg2pr71G6UMP0b7pE1Rb27C1b0tJwZ6eji09DXt6Ovb0DP0zLQ17RrrOS0vT6Rm+9LQ0JCkp9AtxEGxZuZKsZcvglGnwr2v8cirh0H1wzb/B3v86bDTizcoirrBw4ILjmwXAb0LkfQiEvW92UHbLIpIGzAMmoMO2bx3M7ElEEoB7gDOBLGA38AOl1Cu+/DOAPwOT0KrE67p1m7669wKXAW3Ab/yjO/ZXd8TpR0DVBljvZRYkYY+zDbhEELgG1dbkxutVlv1RhtjB09JK29oPaH1/Da1r1uDeq60+k9B/zIPFlpRE3IQiHEVFxBUWElc0gbiiIuJLSoifMhl7ysAfQaPCnIvhxG/De3/sTSt9D16/Hc77VeT6ZRhJnITeY2sHwvZCFLaAEpEfA99Fm+d0vyWbReS3SqmfD+J+B9Fh3w8A5wP/EpH5aP9NzwA3oI0yfgY8CSz11b0TmAGUAAXA2yLyqVLqVZ+xRn91R472BqsXc5sDsqf3XNaUWQVU9oTwXhxx8XYSkuJwtendUl6vor3Z3ce6zxCdKKVw7dxJ6+rVtLyzmrYNG4ZmuCBCXHExzrlzcc6dQ+LcuSTMno09I2NYZjqjwum3Q+XHVku+tfdCZgks/XrIaoaYZRtwEdqtUSAXATuCpAclLAElIj8Bbkfbuj8BVAH5wHLgJyLiUErdOVA7SqlWtKDp5kUR2QcsQvtu2qqUesp3zzuBGhGZpZTaDlwDXK+UqgfqReRvwHXAq8DnB6g7chwOaD57Bjh6148CBVROcfhftskZ8T0CCrRXcyOgopuOHTtoePppml99bdDrR468PBIXLCB+6lScs2fhnDULR1ERtvj4gStHM3YHXPoPuO9Ua1DPV78PaRNgzkWR65thJPgL8FcRaQL+hta2TQBuAr4CfCPchiQct0giUgE8ppTq4zFCRO4CrlJKFYV7U7+6+Wiz9WPQC2fxSqmv++VvAe5A+/6rAwqUUlW+vMuAO5RS80Xk7lB1lVJPB9zzJvQvCofDsWjFihWD7baFovJXmLnrLz3X1bkn8enc3l/Tnte9dNT1li85VUgpFDweoaLC6nMvkMoNXlr8NvPmHy2kTRz4q7mtrY2kpOHd1DtxYvuwtjcUWlpaSIkW1ZUf9vJyErZuxbl+A3Gl4WmVlc1G5+TJtEydgm3aNLoKC/Hk5Vm2J8QywZ5VSvNuFm78AXZv73pql93JhmN/R1tycWATUUm0/g0eCaeddlqbUmpYnX/75MK3sZp+KLRRXdieh8JV8aUT2pPEq2jhMihEJA54DHhIKbVdRFKAwwHFGtFxQ1L8rgPz8OWHqmtBKXUfcB+A0+lUy5YtG2zXrbz0ouUyb/4y8k7RbXo8Xrb/3zvooJKa084/geT0BDo7YV/A9qlAPm46wJbKXvVhbnoRxyyZ1E8Nzdq1a1myZEnYQwiHmTOHtbkhsXLlSo74eQ0TXrebllWrqH/kUdo+/DCsOvGTJ5N88skkH388SZ9ZjD0lJarGNJwEH9cymDMVHrsCvFoz4PB08Jl9d8ONb0FCWCGCIspYfV7DjVLqVhG5FzgLbW9QA7yhlNrbf00r4QqotcBitMO/QBb78sPG54L9EcAN3OxLbgEC7a/TgGZfXvd1R0DeQHVHln4MJOor2/B09QqnpPT4QanoUrOtZZvqzN6XSNO+ZSv1jzxC85tv4m1p6besJCaSvHQpySefRMrJJxM/cWK/5ccF006Hc38JL9/am1azE/51rXY2G4OWfYbgKKX2oLckDZmQAqo7joeP/wCeFZEu4Cl616CuAL4MXBzuDUWv7N7vq3++Uqp7kWUrcK1fuWRgGnptqV5EKtHmi906uQW+Ov3WDbdfQ0KpIHugZvecBkbQzZ00uC/ElCxrvKhmI6Aigtftpunll2l87t+0ffBB/4VFSD7+eDIuu5SUM87AlmDWDPuw+AYoWwefPNGbtudNeP6bcMm9/e8tMsQEIpKElg2nomdQtcBK4EGlVNiGq/3NoLqwuioU4Fe+g4D0TwZoy5970Zu1zlRK+S9sPAv8VkQuRVt//Bj4xM/I4WHgRyKyDi3cbqQ3FP1AdUeGlmpo91tgciRCxuSey5qDAQJq4uAEVGp2QFyoug6UUrFjvTUG6Nixk/JbbukxCw9F0uLFpJ59Nqmnn0bchAmj1LsYRQQu+IP2wFKxsTd90z8htQDOvDNSPTMMAyJSgBZGM9E2BoeAqegtQt8UkWXdtgQD0Z9Q+SmhfekOCREpAb4KuIBDfi/aryqlHvMJmP8FHkWrDa/0q34HWriVAu3Ar5VSrwIopQ4PUHdkCOZBwm+R+3A/Aiqcza+JKXE44u10uXUIlU6Xh47WThJTRt+qS6nx9WHbVVtLzZ/vof7JJ8ETPISNPSeHtLPPIuPyy3HOnh20jCEE8Ulw1VNw/1nWWGrv/gHSi/UsyxCr/AbIBE5WSr3XneiLqvs08Gu0BfaAhBRQ4ZiNDxbfxtmQrzml1Bto/37B8lzoKWPQuCP91R0xqrZYr/3Wn5RXcTjAxVHOpMFZ/4gIqVlO6g+19qQ11XRERECNF7xtbVT96tc0PPssdAaP2JW0eDHZN3yF5JNOQuz9uOUx9E9KLnzpabj/bGir6U1/+XuQlA1zPxe5vhmOhPOA//IXTgBKqfdF5Ef01cKFxETAOxL81RMABfN7Thuq2+hy9X55JyQ7SM3q36w8GOm5iRYBVX+otd9ovIah4XW7aX7lFWru+xvuPcHXddPOP4/sG280s6XhJHsaXPUveOgC6PQtTSgvPHW9DmOz+CuR7Z9hKKQAFSHyyugTiyU0g/EkMRutQ5yIdmXhj1JKXdu31hgnUEAVHdtzGuggNndi6pDWjrIKk9m/uffrsr6ytZ/ShsGivF4a//08h//4R7qqgqvF4yZOJO97t5J61llm/W8kKF4Elz0ATyzXwgkApWdSOTNhyskR7Z5h0OwArkZvQQrkS+go7GERrieJa9AREBU6vEaga+3x5060vd4aRVfslhlUHwu+AAOJcB2wZhVaN9zWVQyfQ9HBMBbXoNo++oiqX/6Kjk8/DZrvKCwk+7pryVi+PPa9OUQ7R50LF98D//5/oHyaB+WBf16phdfMsyPbP8NguAt42OeI4XGgEu2e7kq0H9arw20o3BnU7cC/ga8opRoG19cxSsXH1uu82Xrh10cfA4lBmph3k1lo3eDdcLiNrk4Pjjiz9jFU3AcOUP3bu2gO5UVEhLQLL6Dgx3dgTxnWDfaG/jhmOTjT4Ymr6PnmdbfAP78A5//WGE7ECEqpR31m5j/FGgq+CviaUurxcNsKV0AV+Bo2wqmbPuq9Y3pOlVLDJqASkuJIyUygpV67h1FeRUNV+6B8+hk0XfX11P7979Q//AgqiAGExMeTccUVZF39JeJLSiLQQwOzzofTfwhv+fmfVl546bvQVKEdz461qfwYRCl1n4j8HTgKvQ+qDtihlPL2X9NKuALqPfTepTcH1cuxTMUG63XRwp7Txup2XK29nqvjnHbSc62bbgdDVmFyj4ACqKtsHXUBFcsxodq3bqX+scdpevFFlDt44MfU884l77u3El9s9jBFnFO+B/Gp8NptfmtSwOrfQfMhuPBu43EiBvAJI4urHRE5E/i9UurocNoIV0DdDDwjIrXA60B9iM6MHwJVfH4GEpV7+oZ4l4A4ToN54WcWJnPg094NwXWVLei9yoZQqK4uml9/nbpHHqV948aQ5Zzz55N/2/dJOvbYkGUMEWDp13Q4jv/7cq91H8DHj0FLFVz+ECQYLUIMkg7MDbdwuAKqDNiI3gQbDDWItmKflsPWsAG2OMjv/Z0f2hsgoKamH9Htsoqs6yDGki80npZWGp95mroHH6KzIpSlKzgKCsj7zi2kXXABMkY8iI85jjoPrn0BHr8C2mp703e/oc3Sr3pK76UyjFnCFSp/A74APIc2EQyuJxkvBK4/5c8FR6/PtT4CatoRCqgAQ4n6qjY8Hi92u3mxdtNZVUX9o49S/8STeJtD+wiOmzCBzKuuIvOq5dgSh652NYwSxcfBl1+HRz8PDX6hTCo2ai8UX3pa76UyjEnCFVAXA99TSt09kp2JGQLXnyb0qofam93U+c9wBPKng5+FiAAAIABJREFUHJmASkyJJyk1jrZmX3TdLkXT4TYyC0ZPxRGta1AdO3ZS98ADNL70UkjPDwDJJ5xA5pe+SMqppxrvD7FGznS44Q147DKo3NSbXr9Pe6H44lOW/0HD2CFcAdUKBN8sMh45sMZ67WcgsX9zjWVXWE5xCgmJfX/Ng33hZxUm09bca0RZV9E6qgIqmlAeD42r17D30Zdp2FVOXGcrmR6FHXDHJaPETry7CVt8POkXX0zWddeSMM18Zcc0KXlw3Uvwr2tgz1u96W018OAF8KX/g5ITIte/cY6ITA2zaMFg2g1XQD0AXEVvqIvxi6cTDgYEqCs5sed036YaS9bko3OG5baZRSmU7fQTUJVtjKdXrlIKx4EDVP7sF2xdW8PugjPxOM7UtqWAvasDu6cDd0IGAPG2LvKmZJA7LYuEpAJM0IsxQEIqLH9Sh+XwD9XR2QpPfglueBOypkSuf2MUEfkHcAFQrZSaF6LYbsJz2CBhlgPCF1ClwHIRWYF2XxHMiu8f4d40pqnYaLUqSimALP3x0On2cHBbnaX41AXDs4jbZx1qnBhKuPbupf6FV9j27kGqVT47k6fTXnx8n3IehxOPo9cDl9vroGxPC2V7Wvj4jYOcdPl05p9a3Mea0hBjOOLhc3+BtELt+bybtlr4+5nwhUfMTGr4eRAdKeLhfspc30/ekBlMDCeAEuCMIPkK7Qpp7LN/tfV68ok9GwdLN9fS5e61tk/JTCBnYnA13GBVfIEeJeoOteH1Kmyj9MIdjTUob2srHTt30rF5C+2bN9O2aRMHXQXsnXIhruyhO2hVXsXqJ3fxyVtlzF9WzKzjC0hIMvtoYhYRHTPKHg+rft2b3lYDD10I598Fx43I+3JcopR6R0QmD1DmoZG4d7gCysybu9n/nvV68kk9p7vXW52NTluYN2zORZPT40lIisPVpg0ButwemmvbSc9NGqDm6KGUQrndqI4OvB0ulKsDb0cHyuXC296u01tb8ba24mlpoetQFZ0VFXSWl9NZXo6noVeF6bE52DnjSiqn9p0tdSN4ySuIp65e0ekKHrPJn8bD7bz71C7WvrCX0744ixmLzV6ymObU70PdPtj8r940bxe8+G2o2qpDy5sNveHg8AWC7eY+pdR9EeuNH2EJKF8cJ4OnCw6utaaVaAHl7uiidHOtJWv6cXnDdmsRIasomcrdvS/x2orWYRdQXrebrsOH9VF9mK7D1XTV1hHn3Y3d3aaFjculhY2rA9XRfe5CuVxHPNVyx6WwZ8pFVBad2G+5/MmpnHbNbLKLUvB0eWmu7cDrUaTnJdLR0kn1gWZqy5pZ90opnk7rHvLODg+v37+VtiY3C86YeET9NUQQmw0+91etYl8VEGLoo79B2Yfw2T9ob+mG/uhSSh0X6U4EY/xsrh0OKjdp55XdJOdBzgwA9n9SQ5ffizA1y0n+lOGN25RVaBVQVfuahrzGpTweXLt349q9B9fePbj37qOz6hDeuj7LiwA0xZUTL6HNuIeDtsRcNi74Fi5nZp88m02RXuxh0RlHkzMxheyiXtWp3WEjI79XUCdnJDAlI4EpR+cw5ZhcNrxayu711Xg9VuH57lO72PfJYY45cxIlc7PN+lQsYrPBabdB/hx49mvW9eHKTfD3M2DRdXDGjyEpK2LdNAyNcMNt7GMAywulVLhmhrFL6bvW65ITetafdq2rtmRNX9S/em8oE42CqWlsXV3ec125uwGlVNhqROX10rFtGy0rV9Hy7rt464MLo0jgFRtb53w5qHDKyE/ivK/O55OdH3HUkkFZqZJdlMJZX57LiZfNYPPKMja8VmoRVOU7Gijf0UBmQRInXT6DSXOzj3gshggw52I9k/rnVdB4wC9DwfoHYOerOmxHSWiVsSH6CHcGtYq+AiobOAFoAd7qUyMEInIzOh79fOCfSqnr/PLOAP4MTALWAtd1qxdFJAFtrHEZ0Ab8Rin1+3DqDht7V1mvfetPrrZODnw6cuq9bvImpeGIs/XM1Nqa3DRUt5GZ339ICOX10rp6NXWPPkrnwbJh75c/EheHOJ2IMwFbgt/PRCe2BCe25OSew5GTTdyECaicQt79QNG8o68HiFnHF3DSFTP1XrKdQ+9XUlo8Sy6aStGMDF6655M+ar/6Q228+OdPOPP62cxcPDghaIgSCubDTW/Dy7fC1metec2V8OBn4ayfwPE3G4/og0BE/gksA3JEpAy4Qyl1/2jcO9w1qOuCpYtIBtrs/I1B3LMC+DlwDtDja0ZEcoBngBuAF4CfAU8CS31F7gRmoC0JC4C3ReRTpdSrYdQ9cjoaYd871rQppwB675O3q1d+p+cmDjm8Rn/Y42wUTEmz7Icq31EfUkAppWj9YA31Dz+Ke9++8G4igj07G0duLnH5eThyc7FnZzNhlhNneiKS4MTmTNBCKCEBm9OJOJ3YEhKQhIRBe2noaOnk+T993CfA48TZmZz3taOJSxherw8TZ2dx4TcXsOrxHdQfsgZ/VF7Fivs/ZcuqchacMZFpC4f/I8MwwiTnwOUPwrHXwEu3Qt2e3jzlgdd/BHvehmXfh4mfiVg3Ywml1PJI3fuI1qCUUg0i8lvgF+jIieHUeQZARI4Div2yPg9sVUo95cu/E6gRkVlKqe3ANcD1Sql6oF5E/oaeib0aRt0jZ9cK8PqtwWRP1+GogV3rrNZ7048bPuu9QApnZFoEVOnmOuadUtynnPvAAap+93vcO3aEbEsSE0laeAwJ06cTN3Uq8ZMnE5edjcT1tXxKLYGEYd7t6vV4efkvn/QRTs6UOM68fu6wC6duJszMZPkdSzh8oJmPVxzoo56t3N1I5e5Gpi3M5dhzS8grGd61RMMoMO10+MYaePsX8F6Ah7Y9b+rjlP+EU261+NE0RBeijtDqSkQ+CzyplBqU3x0R+TlQ3D07E5G7gXil1Nf9ymwB7kCrEOuAAqVUlS/vMvRUc35/dZVSTwfc9ybgJgCHw7FoRaioqgHM2fob8g73mpgfmPh59k67li6XYsdzyqIAnXau4MzoX0C1t9upqRl8GPGuDsW+t6z3m3SykJDWez/18cdkv/QSNndfowYl4J49B9f8+bhnzgBHeN8o+fkdxMcP72aoqk1earZZ0xLSofgEwZne9/fX0tJCSsrwu3dq2Kco/1CFXGXNmQN582REjChGakyRJprGlV3zIbO2/5G4rr6b213xmZSWfIHKwrNQtoH/F6JpXMPFaaed1qaUisrQ0UOeQYmIA5iHVr1tHYa+pACHA9IagVRfXvd1YN5AdS347PvvA3A6nWrZsmUD96yzA97fZEmadPY3mDRxMVtWlbFD9S6OZBYkcc7FSwacQbW0QD/RIPrFVbqNCj9rvmRPAYuWTAag4bnnqH32ub6VBFJOPZXML36R+ImDN60uGeYZVOnWWrZus/5OC6enc+E3jwk5c1q5ciVhPa/BsgyqT29i9ZO7+niiB6j5FJzeDE66fAbZE4b35TRiY4ow0TWuZVB/BTx9ozY99yPBXc/MXX9hZu0Kbek35+J+16eia1xjn7DiNYiIV0Q8/gfgAtYD04FbhqEvLUCgLiUNaPblEZDfnTdQ3SNn3ztW8/KUApig91Zse7/SUnTG4vwRU+91M3mB1b/fno2H6XR7aF23jtq//rVP+aSlSym+5x7yv//9IQmn4UQpxe711bx63xZLenJ6POfeNH/E1HoDkVeSxqX/uYjP3XosqVnOPvll2+t54ucf8sG/99B4uB13R1eQVgxRS+Zk+MrrcOU/IbGvpSh1e+Cpa7VZ+tZnoWt8RxSKFsKdQf2UvgqQDrSPvleUUn0/OwfPVuDa7gsRSQamodeW6kWkElhAr8PaBfTO3ELWHYZ+wfYXrNezzgebjdqKFqpL/WSgwFFLR94CbNLsLNb7eZVwt3exe9UuEv76K+tTstvJ/tpXSb/ggiMWmsPh6qir08MbD3zKng3Wya4InH3DXJLSBq/yHG6Kpmdw1U+WsGd9Ne89s4f2Jr8XlYL1r5Sy/pVSEMidmMqck4qYtbQAR7wJ4RH1iOj/3aI18P6fYN0D0NVuLVO+Hp66Tu9xXPhFOPZa44A2goRrxXfncN3Qpxp0AHbALiJOoAt4FvitiFwKvAT8GPjEz8jhYeBHPpcc+cCN9DooHKju0Olsh0+ft6bN+iwA29ccsiQXH5VJWnZ4QfCO5IXviLczfVEuW1f36gg3r6rk6PYueswbRCi48w6SFy8e+o2GEXdHF6/8ZTNl2/vuvVr6uWkUzQjyVRshHHF2jlpaSPHsLF7/+1YqdjX0LaTg8IFmVj2+gw9f2Mv8ZcXMP7UYZ4pxrRP1pBVqN0gn/If2QLHhEW3h509rtXZG++4ftMHF4hth5rmR6e84JhIhWX8EtAPfB77kO/+RUuowcCnaIrAeWAJc6VfvDmAPeta2CvitUupVgDDqDp1Pn4cOvxdUYhZMPgWPx8uOtVYBNev4wmG5ZTjMWFyA3dH7+Nw2JweLe/34Zl1/XdQIpy63hxf/Z1Mf4SQCJ10xg2PPLolQz/onOT2BS25ZyLlfndfv7K69uZMPX9jHQz94j7XP78Xr8YYsa4gi0grhwrvhGx/ArAtCl9vzFjyxHP44n6l7HoLmqtBlDcNKyBmUiPx4MA0ppX4aZrk70YYVwfLeAGaFyHMBX/Ydg6p7RKx/0Hp9zFXgiOfAJzUW9U+8087UhcMTWiMcUjISmHtKEZ+81bvxtipvEXk1HxM3KY2Myy8ftb70h1KKlY/toHKPVQucWZjM2V+ZS05xdFtEiU2YtjCPohkZfPzGQfZuPEx7sxtXW981qC63l3Uv72fPhmrmLytm9omFOOKM6i/qyZ0JVz4Gh7ZorxOf/AtcTX3LNZUxqakM/viS3md14rcgw/hyHEn6U/HdGUZ9f0VVWAIqpji8Ew68b007Vi91bX77oCV5+nH5xI3yOsTsEwrZ/fZO2pTPD50IO6dfRvHJ9mE31BiqSnLzyrI+M83cSalc+B8LSEyJ/JpTuCSmxHP8JdM4/hIdJrKlvoNP3ipj6+py3B1W9VD9oTbeeWIn61/Zz7HnljDnpCIjqGKBgnnw2d/BWT+FLc9oYVW+vm85j0s7o13/ACxYrkN/JA9PYFKDlf5UfHEDHIuB19EREnePbDcjxIaAECclJ0LuTGrLWzi4zaqumnty0aCaHg6jg/a33qBk2zOWNFdCBmU7UnsMKCLJlnfKee8p659GZkESF3/7mJgSTsFIyXRywqXTueaXJ3LC56eTGEQF2NroZvWTu3j09g/Ys6E6SCuGqCQ+GY69Gm58C25aBfOvAEeQtWVvF2x8BO49Ad65CxpH1o3YeCSkgFJKeYIdwFTgUbS/uznoTa9zRqe7o0hnB3wc4Bxj0XUArHtlvyW5cHr6qHsbcFdWcvjev5DRtJf8qnWWPFcjrHxsO13ugWMkjRTrXtnPqsd34PX2SuJ4p53zvjZ/TAULTEh0sPDsSXzhB4spnhXc0KO1wcWr923hg3/vCZpviGKKjoFL/wb/uRfOv4uOhCAzpZYqeOtn8Ie58NdTfdaBrtHv6xgkbCMJEZkoIn9Hm26fDtwKzFBK/d0nuMYW21+Edr/w7c4MmH0R1aVN7A5wjTPaMYWUUtT86U/g0v8Ekw+8QnrzfkuZwwdbeOfJnRHZr7NlVRlr/73XmihwxnVzyCyIyg3rR0xyRgIXf3shV//8eI45cyKOuL7/WutfKeXVv26mvcXssYk54pPgMzeydslf4KL/0Z7Tg1H5sQ6Y+Mej4b0/QUeQtSxD2AwooEQkz+dKaCfaUu6nwFSl1B+VUmP3Py2IcYSHeN56xGq5njMxZcgxmYZKw5P/on3jxz3XNqU48cwM8idbZ3EVuxp47W9baKppD2xi0ISrktz54SFWPWF1Ox6XYOecG+Yx9ZjR/T1FgrScRE68bAZX/+IEFpw5EVuAe6Q9Gw/z3O83mo2+MYqyxWkDif/3EZz6XzrsfDBaDsGK2+GumdqDxbYXwTU8fgPGEyEFlIiki8h/o027vwLcjRZMP1dK9XVqNZYofR/2r7amHXstH764l9qyFkvy8ZdMG5KPtqGuQbWuW0/dQw9a0pwLFpB5/rmcunwmGXnWCLuNh9t56Z5P2LjiAIcPNltUbsPNng3VvPngNovpjCPOxgXfXMD0RePLM3hSWjwnXTaDz916bB/vGHUVrTx2xwdsXlmGu1lxpP4wDRHA7oDTfgDf3QHn/bbHs0wfutp1SPonvwi/mwUvfge2v2yEVZj0Z8W3D0hHG0L8HKgEMkUkqKJdKbU3WHrMoRS8frs1reREKhrz2fD6RkvyjOPyRjXAXVd9PdW//51FAEhKCrnf/hYiQnxiHKdfO5uX7luPy8+q29PlZevqcrauLiclM4FZSwvJK0kloyC5zxf+UNn01kHefWqXpW82m3DuV+dTND1jWO4RixRMTeeS7yxkxT8+paGqN7xHW6Obd3wzzdqNH/aYpdvtkdiaaBgySVmw5CZ9NJbBmnu0dV9nW9+y7hZYd78+HE446jw4+kqYfgbYx8667HDSn4DqfqucA5wdRltjw45267NQHmB0cOLtvPGAdWaQnJHAKcuPGrVuqa4uqn71K2tIdhEKfvgD4gt7NwgnpcZTfLzgKc9m/+aaPu201Lt6jDziE+zkTEolrySV3JI0copSsAdZO+m3X17F+8/s5uM3rGb3CJz55TmUzDMRavNK0rj8tuP4v1+vp76yr/KhrqKVVY/vYMuqcs69aZ4lfL0hhkgvhnP/W4fw+PBv2sKv8WDwsl0d+l2z9VlIyoH5l8HRX4CihSaYoh/9Cajr+8kbm7ha4I07LUnu6Zfw2itOmuusZuVnXDcbZ/LoffXUPfwwHZs+saRlLl9O0sKFfcra7MLSy6ZTOCOdja8doKM1uMm52+WhYldDjysfR5yNgqnpFM3MJHdSCum5ST0zrGBaqK5OD28+uI3d661GI2ITTr96FjOOyx/KUMck8U4Hl9yykHef2sXu9dWoIKrW2vIWHrvjAybMzGD2iUUUz8okOd3EKoo5krJg2X/BKd/TSwU7XoYtT0NrYMAFH201sPYv+siapuPMTVgEk0/UP8dxvKqQAkop9VCovDGJUvD8zdDQGyW+yZPP89u+QmOtVTgtOH0iE2dlHfHtwqX5rbdp+NdTljTn/PlkXhU60KWIMO2YPCbNyWbfphoqdtRRsbux3zWork4vZTvqKduhxxsXbydrQjITZ2eTsDSFxOm9KsHGw+28fv9WqvdbrZTiEuycc9M8SkZR9RkrJKXFc/ZX5nLC56ez7f0KynfUU76roY8b5vKdDZTvbMBmkx7VX1ZR8oh7yTcMMzYbTD1VH2f+RLtM2vs27HwVGg4Er1O3Rx87X9HXDicUL9Z7MCefqM/jwvP3ORY4ooi6Y4o1f9bTbcCrbLzb/BU2t52P9mPbS/aEZJZeEsLEdARoWf0u1XfdZUmzZ2WSf9v3kTCCDcbF25m5OJ+Zi/Npqe/g4LY6qvc3UX2gZcDNvJ1uD1X7mqja18TuVS4ycxwc9ZkC0nKcvP/sHjoDPCgkpcVzwc0LRiTc/VgiJTOBxZ+dwuLPTmHFy2/j3pfF/s21fcp5vYpNbx1k01sHySxIYurCXCbOyqJoRsaIBE80jCBxTu1Jfdb5cN5v4MAHsOmfsPU5LAvGgXR16FnY/tXaA6k9XqsBi47VP4uP0ybvY/TjxQgogH2rYYV2Pdip4nmt4XuUuo7rUyx3UioX3Lxg1EIrNL2+gsN33w1eP+ejdjv5P/whjqzBz+BSMp3MPqGI2ScUobyKxsPtVB9o1gJrfyNtzf0LrLZGNxtXBP/yyyxI4oKbF5CWM36+7oaDuCThjK/N5+M3DvLpexU0VgffElB/qK0n1Ed2cQpLLpzCpHnZxqgiFhGBkuP1cd5v9Ixq0xOwe4X2TtEfHjccXKuPbrKmwszzYObZWmg500e2/6OIEVCN5fB/14Py4FEOXqj7MZWdc/sUS89NHFb/cf2p+DytbdTc82da3nzLmiFC3vduJXFu3/4NFrEJGflJZOQnMXNxPkop6g+1Ur6zgZqDzdSUtQbMsEJ/oU1ZkMPp14zumtxYwma3cew5JRx7TgnVpU3s+qiKnR9V0dYYfJthbVkLL9+7GZtDmLogl8lH55CRn0RmQRLxTvMvHVPEOWHuJfpoq4OanXB4B5S+B/vfg6Yw3CfV7YUP/qwPgLRiyJvtO+bon7lHxaRqcFz/Ndu9btTv5/Jhy5Xs7jiRBs+EoOXmnFTEZy6YMuL+41RXF61rPqDuwQfpLC+3ZgrkfucWUkco3LSIkFWYQlah9i6ulKKlvoPSLbVUlzbjqqnpEzPH7rBx4mXTmXfqBLM+MkzklaSRV5LGZy6ayu51Vexef5iybXVB1w69XTo6sb+RSlpuIkUzMpg8P5vM/GRSc5yj7sTYMESSsmDSUn0sulZ/xTaUakFV+h7sf9eyRh6SpjJ97F7RmyY2yJzSK6wyp+gow5mTR2o0w8K4FlCg2O86jnWtVwTNzS5O4cJvLhhxSyrl8dCyahV1jz5KV0Vl3wIJ8eR961uknn76iPbDHxEhNSuReacUA5Cb7aF692G2r6nk0N5GcopTOfWqmeQUm/WmkSAu3t6jju1o7WT/5hr2baph78YQlmA+mg6303S4ne3v9/4dJabFM2FmBjnFKZTMyyF7gjG4iAlEeoXIwi/qtMZyqNgA5Rv0dpjSNeANwzG08vYaYGx/cSR7PayMcwEFB91HB01Py3FyyS0LR1Rt5XW5aH7zTRqffqbvjMlHfMkk8m+7jfjJk0esH+HgiLdz1JICjlpSgFLKvOBGEWdyHLOWFjJraSHVpU2sf7WU8h31QWNSBaO9yc3uddXsXlfNB8/tJT7RQbzTTnpeEjkTUsguTiZ7QgpZhckmdH20kz5BH7Mv1NeuZm0duPM1HRqkZlff6MAxzLgXUDWdfS3yUrOcnP/1o0dEOHVWVNC8ciVlr31E+do94A799ZN67rnkfO2r2JzOYe/HkWCEU+TIK0njvK/OR3kVpVtrKdtRT2N1O/WHWmk83N7HZD0Y7vYu3O1dtNS7KN/Ru4VCBDLyk8gqSiZnYio5E1JwJNhJSosnKTWehCSHsR6MNhJSYc7F+gDtRb12N1Rvg+pPe3/W749oN4fKmBJQIpIF3I/2fFED3KaUejxUeY8tnsMyF///6s/feiwFU9OP+B9ReTx0VlTgLj2A+0ApnaUHaP3wQ1zbtgHQ6kkHTxCvUTYbKWedScZFF5EwbdoR9cEwdhGbMHl+DpPn94Z/6HR7qNzVwP4ttdRXttJQ1UZbkxuvJ7xNd0ppa8H6Q23s2dBXlSgCCclxOJPjSEyJ8507etK6j4RkB+31iua6DpzJcTjibHg8XpRXbyK32cV85IwUjgTIn6sPf1wtULNDC6y6fVpg1e+H+n1A9HpcH1MCCvgz4AbygWOAl0Rkk1Jqa9DSHuhy9/7zJsR5cW5eSeMWX4L//3W32Z3F/K73XHV20llWhmvfftz799N54ACqcxBBA+12Uk47jcyrlhNfNLjghwYD6HWrSXOzLf4hPR4v1fuaqClroXRrLQe31eHtGppzWqWgo6WTjpZOGqoGLr/3tfdD5tlsWlDZHDZsdsFuF2x2GzaH/ml3iK+M77y7rE2w++r01vX9dPgEn/jN8n3nPZeiE8XWfe1XLlS6X/263Yqt9nLEptN68mJC4E7QRwr66IkSVBiyRqQZMwJKRJLR4UDmKaVagHdF5HngauD7wep4A0ynk6p2cOi2/x3prlqw5+aScvLJpF9yMXF50evx2zjcjk3sdhuF0zMonJ7B/GXFdLo8tDW58HoUteWt1Ja3UFveQk1ZC821HaPWL69XacvETu/AhaOMynU7It2FccOYEVDATMCjlPIPRrQJODXcBlJbQjh2HE5ESFy4EPvic0iZczLxkycbdYdh1IhLsJOeq53RZhYkW8KguNu7qK3QQqtiVwOutk46Ozy0Nblpb+nE3W5iWBlGl7EkoFKAQJ8hjYDFDlpEbkKHqScuwAljbo3VGeuR4k1NoSs3D09eLp7cXLry8nAfdRQqLY3GRgdN1dVQXT1wQ4Okra2NtWvXDlxwEGRluUlOjqx1UEtLCytXroxoH4abaB1T/DTo3vXX7bPE6xE8bvC4tEMD//Mul7KkdXZ4UV02PC5t4Sw2QPR5OIYcBgOMLQHVAqQFpKUBlshgSqn7gPsAEuISVLGjApt4KXIcomDZPGBewIzG71yCnPvpqh15eSRMnkz8lCnEl5RgTw/tcqS2Vh8jwdq1a1myZMmwtllQAGmBv91RZuXKlSwboY3KkWIsjgl6x9UdjNH/f0r51HueLi9ej/Id+tw/zePx4u3yy+su50vzdJfrqeNFeZV1uVgpLQ+V3nyOIsS1r5wXFNZyKN2mAirKyyksKNTX3elKacEbq/w10h0IzVgSUDsBh4jMUErt8qUtAIIbSABiFy7+3y+NSucCibU1nVjrryE6CKa+Fptg9xk7xBorV1aybNnsSHdjeLkh0h0ITez9hYTAF4b+GeCnIpIsIicCFwOPRLZnBoPBYBgKY0ZA+fgGkAhUA/8Evh7SxNxgMBgMUc1YUvGhlKoDLol0PwwGg8Fw5Iy1GVTMEGtrOrHWX4PBEPsYAWUwGAyGfhGRc0Vkh4jsFpGgjg9GAiOgDAaDwRASEbGj3cidB8wBlovInNG495hag4olHA5IGKEwU3Fx3mFv226iMBgM45XPALuVUnsBROQJtIX0pyN943EtoFwulxKR9kj3YwRwAGPRL81YHNdYHBOYccUSSSKyzu/6Pp9Dg24mAP5+4MqA4fUEEIJxLaCADUqp4yLdieFGRNaZccUGY3FMYMYVS4QxpmDOQkfFbMqsQRkMBoOhP8rwC84BFAMVo3FjI6AMBoOFwHSIAAAGjElEQVTB0B8fATNEZIqIxANXAs+Pxo3Hu4rvvoGLxCRmXLHDWBwTmHHFEv2OSSnVJSI3A68BduAfo+WhR5TZgWkwGAyGKMSo+AwGg8EQlRgBZTAYDIaoxAgog8FgMEQl41JAiUiWiDwrIq0iUioiV0W6T0NBRFaKSIeItPiOHX55V/nG1ioiz4lIVn9tRRIRuVlE1omIS0QeDMg7Q0S2i0ibiLwtIiV+eQki8g8RaRKRQyLynVHvfAhCjUlEJouI8ntmLSJyu19+NI8pQUTu9/1dNYvIRhE5zy8/Vp9VyHHF8vMCEJFHRaTS17+dInKDX170Py8dtnh8HehYUU8CKcBJQCMwN9L9GsI4VgI3BEmfiw51f4pvjI8DT0S6v/2M4/PoMCn3Ag/6pef4ns3lgBP4LfCBX/4vgdVAJjAbOAScG+nxDDCmyehNjo4Q9aJ5TMnAnb4x2IALfH9nk2P8WfU3rph9Xr7+zQUSfOezfP1bFCvPK+K/wAg8sGTADcz0S3sE+FWk+zaEsYQSUP8NPO53Pc035tRI93mA8fw84GV+E/B+wLNrB2b5rsuBs/3yf0aUCeIgYxrohRf1Ywro7yfApWPhWYUY15h5XsBRQCVwRaw8r/Go4psJeJRSO/3SNqG/NGKRX4pIjYi8JyLLfGlz0WMCQCm1B59QjkD/joTAcbQCe4C5IpIJFPnnE1vPsVREykTkARHJAYi1MYlIPvpvaitj6FkFjKubmH1eInKPiLQB29EC6mVi5HmNRwGVgp7a+tMIpEagL0fKfwFT0c4c7wNeEJFpjJ0x9jeOFL/rwLxopgZYDJSgVS2pwGO+vJgZk4jEofv9kFJqO2PkWQUZV8w/L6XUN9B9Ohl4BnARI89rPHqSaAHSAtLS0DrnmEIptdbv8iERWQ6cz9gZY3/jaPG77gjIi1qUUi1At+foKt8O/UoRSSNGxiQiNrRa3A3c7EuO+WcVbFxj4XkBKKU8wLsi8iXg68TI8xqPM6idgENEZvilLcA6nY9VFNrz8Fb0mAAQkalAAnrssUTgOJLR62lblVL1aHXFAr/ysfgcu125SCyMSUQEuB/IBy5VSnX6smL6WfUzrkBi6nkFwYHvuRALzyvSC3cRWix8Am3JlwycSAxa8QEZwDloCxwH8EWgFb0QOhdoQk/pk4FHidKFW99YHL5x/BL9Bds9plzfs7nUl/ZrrJZGvwJWoS2NZqH/qaLCgqqfMS3xPSMbkI22Jn07Fsbk699fgA+AlID0mH1WA4wrZp8XkId27JqC9qF3ju8dcXGsPK+I/xIj9OCygOd8D+sAcFWk+zSEMeSivQw3Aw2+f66z/PKv8o2tFfg3kBXpPvczljvRX6b+x52+vDPRi7vtaKvFyX71EoB/oIVxFfCdSI9loDEBy4F9vudSCTwMFMTImEp84+hAq4G6jy/G+LMKOa4Yf165PiHT4OvfZuBGv/yof17GWazBYDAYopLxuAZlMBgMhhjACCiDwWAwRCVGQBkMBoMhKjECymAwGAxRiRFQBoPBYIhKjIAyGAwGQ1RiBJTBMAh8sYEGOvb7yj7YfW4wGAaP2QdlMAwCEVkakPQs2tPznX5pLqXURp/j3jSl1MbR6p/BMJYYj85iDYYho5T6wP9aRFxATWC6r+yeUeuYwTAGMSo+g2GECFTx+YUP/5qI/NIXSrvZF5Y7SUSmi8hrvrDiu0Xk2iBtLhCR50WkXkTafXHATh7VgRkMo4QRUAbD6HMbOiDctcCPgS+gnZU+C7wEfA4d0fUBEekJEicixwLvo31J3oh29FkLvCEii0ZzAAbDaGBUfAbD6LNHKdU9O3rNNwO6GrhaKfUogIisAy4CLqM3zMFv0Q6AT1dKuX3lXgO2ALcDl4zeEAyGkcfMoAyG0eeVgOvtvp+vdScoHZOnGpgIICKJwKnAU4BXRBwi4kDH/3oDOGWkO20wjDZmBmUwjD71AdfuftKdvvMsdEyf231HH0TEppTyDlcnDYZIYwSUwRAbNABe4M/omER9MMLJMNYwAspgiAGUUq0ishodenuDEUaG8YARUAZD7PAd4B20YcX96AivOcCxgF0p9f1Ids5gGG6MkYTBECMopTYAi9Gm5X8CXgfuBuajBZfBMKYwro4MBoPBEJWYGZTBYDAYohIjoAwGg8EQlRgBZTAYDIaoxAgog8FgMEQlRkAZDAaDISoxAspgMBgMUYkRUAaDwWCISoyAMhgMBkNU8v8BczL/KQtU2sIAAAAASUVORK5CYII=\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",
    "fig = plt.figure( )\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
}