demo_multi-phase_lockdown.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using ESOP to design multi-phase 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 I2, P2, L2 to denote respectively, the initiation point, period and the level of the (third) lockdown that ESOP is trying to optimize. These are referred ensemble as IPL. The previous two lock-downs (that have already happened) are denoted using (I0, P0, L0) and (I1, P1, L1) respectively. In this experiment, all three parameters I2, P2, L2 are being optimized (but (I0, P0, L0) and (I1, P1, L1) are fixed) under constraints that I2 must be no less than I1 + P1 + 10 (i.e. the third lock-down must start no earlier than 10 days after the second lock-down ended) 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 I2, P2, L2 values\n",
    "def getIPL2( x ):\n",
    "    global IPL2c, IPL2w, offset, ul, ll\n",
    "    IPLVal = np.floor( IPL2c + (np.array(x) - 0.5) / 0.5 * IPL2w ).astype(int) + offset.astype(int)\n",
    "    correctedIPLVal = np.minimum( np.maximum( IPLVal, ll ), ul )\n",
    "    return correctedIPLVal\n",
    "\n",
    "# Get a lockdown schedule corresponding to a certain IPL value for the 3rd lock-down.\n",
    "# Make sure to incorporate the previous two lockdowns as well.\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",
    "    LKP = np.zeros((T,))\n",
    "    LKP[I0:I0+P0] = L0\n",
    "    LKP[I1:I1+P1] = L1\n",
    "    I2 = IPL[0]\n",
    "    P2 = IPL[1]\n",
    "    L2 = IPL[2] / 2\n",
    "    LKP[I2:I2+P2] = L2\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 the 3rd lock-down\n",
    "# is initiated as specified in the normalized parameter x\n",
    "def askSimulator( x ):\n",
    "    global 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 = getIPL2( 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 I2, P2, L2 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( getIPL2( 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"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "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",
      "subepoch 2 of 2\n",
      "epoch =  4\n",
      "subepoch 1 of 2\n",
      "subepoch 2 of 2\n"
     ]
    }
   ],
   "source": [
    "# The first two lock-downs have already happened\n",
    "# Their parameters are fixed and not optimizable anymore\n",
    "\n",
    "I0 = 15\n",
    "P0 = 15\n",
    "L0 = 4\n",
    "\n",
    "I1 = 56\n",
    "P1 = 35\n",
    "L1 = 4\n",
    "\n",
    "# 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 I2, P2, L2\n",
    "# have their own centers, widths and offsets specified below\n",
    "IPL2c = np.array( [25, 20, 5] )\n",
    "IPL2w = np.array( [25, 20, 5] )\n",
    "offset = np.array( [100, 0, 0] )\n",
    "\n",
    "# Upper and lower bounds on legal values of lock-down initiation points\n",
    "ul = np.array( [150, 40, 10] )\n",
    "ll = np.array( [100, 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",
    "IPL2w = IPL2w / 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",
    "        IPL2c = 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",
    "    IPL2w = IPL2w / 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 114, last for 25 days, at level 3.500000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XuUXGWZ7/Hvj7QhNAkk5MItQIIGHHBhwObmBRGGm0cTx8NoQJcRPWZm5CKugw4cZkBRlCHeZzw4qAGZwwIRQXIYhImInmGUSxIDEi4S7i2QCwkJEkxI8pw/9ttJdVHdqW6yu/be9fusVatrv/XuXU96VfLkfd633q2IwMzMrGi2a3UAZmZmjThBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZITlBmZlZIXW0OoA8jBs3LiZNmtTqMMzMrIEFCxasiIjxW+tXyQQ1adIk5s+f3+owzMysAUlPNdPPJT4zMyskJygzMyuk3BKUpP0lLap5rJF0dnrtTEmPSFos6dKac86TtCS9dkJN+4mpbYmkc/OK2czMiiO3OaiIeASYCiBpGPBH4EZJ7wGmAwdFxDpJE1KfA4AZwIHAHsAvJO2XLvdd4DigG7hX0tyIeDCv2M3MrPWGapHEscBjEfGUpNnAJRGxDiAilqU+04FrU/sTkpYAh6XXlkTE4wCSrk19naDMzCpsqOagZgDXpOf7Ae+SdLekX0s6NLXvCTxTc053auurvRdJsyTNlzR/+fLl2/wPYGZmQyv3EZSk4cA04Lya9xwDHAEcClwnaV9ADU4PGifR19wGOCIuBy4H6OrqGtRtgpcuXcoll1wymFOHzLBhwzjqqKOYNm1aq0MxM8vVUJT4TgIWRsTSdNwN3BDZvebvkbQJGJfa96o5byLwbHreV/s2tWrVKr71rW/lcelt6utf/zpXX301p556aqtDMTPLzVCU+E5hS3kP4GfAMQBpEcRwYAUwF5ghaXtJk4EpwD3AvcAUSZPTaGxG6tvWbrzxxlaHYGaWq1xHUJI6yVbf/U1N8xxgjqQHgPXAzDSaWizpOrLFDxuA0yNiY7rOGcBtwDBgTkQszjPuMvA8m5lVXa4JKiLWAmPr2tYDH+2j/8XAxQ3abwFuySPGWhMmTOAb3/hG3m8zKM899xyzZ8/efPzCCy+0MBozs/wpG7xUS1dXV1RtL76nn36affbZZ/Px7rvvzrPP5jIVZ2aWK0kLIqJra/281VFJjB3bayDKypUrqeJ/LszMejhBlURnZyfDhw/ffLxu3TrWrl3bwojMzPLlBFUSkl4zivI8lJlVmRNUiTQq85mZVZUTVInssssuvY49gjKzKnOCKhGX+MysnThBlYhLfGbWTpygSsQlPjNrJ05QJeISn5m1EyeoEnGJz8zaiRNUibjEZ2btxAmqRFziM7N24gRVIi7xmVk7cYIqEZf4zKydOEGVSKMR1KZNm1oUjZlZvpygSmT48OGMHDly8/GmTZtYs2ZNCyMyM8uPE1TJuMxnZu3CCapkvJLPzNqFE1TJeCWfmbULJ6iScYnPzNqFE1TJuMRnZu3CCapkXOIzs3bhBFUyLvGZWbtwgioZl/jMrF04QZWMS3xm1i6coErGIygzaxe5JShJ+0taVPNYI+nsmtfPkRSSxqVjSfqOpCWS7pd0SE3fmZIeTY+ZecVcBp6DMrN20ZHXhSPiEWAqgKRhwB+BG9PxXsBxwNM1p5wETEmPw4HLgMMl7QJcCHQBASyQNDciVuUVe5G5xGdm7WKoSnzHAo9FxFPp+JvA58kSTo/pwFWRuQsYLWl34ARgXkSsTElpHnDiEMVdOKNHj0bS5uPVq1ezYcOGFkZkZpaPoUpQM4BrACRNA/4YEffV9dkTeKbmuDu19dXei6RZkuZLmr98+fJtGXuhDBs2jNGjR/dq8yjKzKoo9wQlaTgwDfiJpE7gfOCCRl0btEU/7b0bIi6PiK6I6Bo/fvzrCbnwXOYzs3YwFCOok4CFEbEUeCMwGbhP0pPARGChpN3IRkZ71Zw3EXi2n/a25ZV8ZtYOhiJBnUIq70XE7yNiQkRMiohJZMnnkIh4HpgLfCyt5jsCWB0RzwG3AcdLGiNpDHB8amtbXslnZu0gt1V8AKmkdxzwN010vwV4L7AEWAucBhARKyV9Cbg39bsoItq6puURlJm1g1wTVESsBcb28/qkmucBnN5HvznAnG0dX1l5DsrM2oF3kighl/jMrB04QZWQS3xm1g6coErIJT4zawdOUCXkEp+ZtQMnqBJyic/M2oETVAm5xGdm7cAJqoRc4jOzduAEVUKjRo2io2PLV9heeeUVXnnllRZGZGa27TlBlZAkl/nMrPKcoErKZT4zqzonqJLySj4zqzonqJJyic/Mqs4JqqRc4jOzqnOCKimX+Mys6pygSsolPjOrOieoknKJz8yqzgmqpFziM7Oqc4IqKZf4zKzqnKBKyiU+M6s6J6iSconPzKrOCaqkGpX4IqJF0ZiZbXtOUCU1YsQIRo4cufl4w4YNvPjiiy2MyMxs23KCKrHx48f3Ol6+fHmLIjEz2/acoEqsPkEtW7asRZGYmW17TlAlNmHChF7HHkGZWZU4QZWYS3xmVmW5JShJ+0taVPNYI+lsSbMlPSzpfkk3Shpdc855kpZIekTSCTXtJ6a2JZLOzSvmsqkfQbnEZ2ZVkluCiohHImJqREwF3gasBW4E5gFviYiDgD8A5wFIOgCYARwInAj8b0nDJA0DvgucBBwAnJL6tj2PoMysyjqG6H2OBR6LiKeAp2ra7wJOTs+nA9dGxDrgCUlLgMPSa0si4nEASdemvg8OSeQF5kUSZlZlQzUHNQO4pkH7J4Cfp+d7As/UvNad2vpqb3teJGFmVdZUgpK0g6T9B/MGkoYD04Cf1LWfD2wAru5panB69NNe/z6zJM2XNL9d/qF2ic/MqmyrCUrS+4FFwK3peKqkuQN4j5OAhRGxtOaaM4H3AR+JLfvzdAN71Zw3EXi2n/ZeIuLyiOiKiK76f7iryoskzKzKmhlBfYFsLuhFgIhYBEwawHucQk15T9KJwN8D0yJibU2/ucAMSdtLmgxMAe4B7gWmSJqcRmMzUt+2V5+IV6xYwaZNm1oUjZnZttVMgtoQEasHc3FJncBxwA01zf8CjALmpeXn3wOIiMXAdWSLH24FTo+IjRGxATgDuA14CLgu9W173o/PzKqsmVV8D0g6FRgmaQpwFvCbZi6eRkhj69re1E//i4GLG7TfAtzSzHu2mwkTJvCnP/1p8/Hy5ctfc68oM7MyamYEdSbZd5PWkZXq1gBn5xmUNc8LJcysqrY6gkqjoPPTwwrG34Uys6raaoKSdAcNlnVHxDG5RGQD4u9CmVlVNTMHdU7N8xHAfyf7/pIVgEt8ZlZVzZT4FtQ1/ZekX+cUjw2QvwtlZlXVTImvdknYdmQbv+6WW0Q2IB5BmVlVNVPiW8CWLYc2AE8An8wzKGueF0mYWVU1U+KbPBSB2OB4kYSZVVWfCUrSB/s7MSJu6O91Gxou8ZlZVfU3gnp/P68FvbcvshZplKA2bdrEdtsN1Z1UzMzy0WeCiojThjIQG5wRI0YwatQoXnrpJQA2btzIiy++6O2OzKz0mrqjrqT/Rrbd0Yietoi4KK+gbGDGjx+/OUFBtlDCCcrMyq6Z+0F9D/gw2Z58Av4a2CfnuGwAvFDCzKqomYmKt0fEx4BVEfFF4Eh630DQWswLJcysippJUK+kn2sl7QG8CnjpeYF4Nwkzq6Jm5qBuljQamA0sJFvB9/1co7IB8QjKzKqomS/qfik9/amkm4ERg73DruXDu0mYWRU1s0jiPkn/S9IbI2Kdk1PxeJGEmVVRM3NQ08j24LtO0r2SzpG0d85x2QC4xGdmVbTVBBURT0XEpRHxNuBU4CCyDWOtILxIwsyqqNkv6k4CPkT2faiNwOfzC8kGyiMoM6uiZu4HdTfwBuA64K8j4vHco7IBqU9QK1as8H58ZlZ6zYygZkbEw7lHYoO2/fbbs9NOO7FmzRog249v1apVjB07tsWRmZkNXjNzUE5OJeAyn5lVjWtAFeGFEmZWNU5QFeERlJlVTTNf1O2U9I+Svp+Op0h6X/6h2UA4QZlZ1TQzgroCWEe2izlAN/DlrZ0kaX9Ji2oeaySdLWkXSfMkPZp+jkn9Jek7kpZIul/SITXXmpn6Pypp5iD+nJXnEp+ZVU0zCeqNEXEp2S7mRMQrZPeF6ldEPBIRUyNiKvA2YC1wI3AucHtETAFuT8cAJwFT0mMWcBmApF2AC4HDgcOAC3uSmm3hEZSZVU0zCWq9pB3IdjFH0hvJRlQDcSzwWEQ8BUwHfpTafwR8ID2fDlwVmbuA0ZJ2B04A5kXEyohYBcwDThzg+1eeN4w1s6pp5ntQXwBuBfaSdDXwDuDjA3yfGcA16fmuEfEcQEQ8J6mnNrUn8EzNOd2pra/2XiTNIht5sffe7bdVoDeMNbOqaeZ2G/8haQFwBFlp7zMRsaLZN5A0nGzD2fO21rXR2/fTXh/n5cDlAF1dXa95vepc4jOzqmlmFd9c4HjgVxFx80CSU3ISsDAilqbjpal0R/rZU4vqpvet5CcCz/bTbjW8SMLMqqaZOaivA+8CHpT0E0knSxoxgPc4hS3lPYC5QM9KvJnATTXtH0ur+Y4AVqdS4G3A8ZLGpMURx6c2qzFu3Lhexz378ZmZlVUzWx39OiI+DexLVkL7EFtGPf2S1AkcB9xQ03wJcJykR9Nrl6T2W4DHgSVkt5T/dHr/lcCXgHvT46LUZjV69uPrsWnTJlau9K/JzMqr2dtt7AC8n+x2G4ewZRVevyJiLTC2ru0FslV99X0DOL2P68wB5jTznu1swoQJmzeMhWweqn5kZWZWFs3MQf0YeAg4Bvgu2feizsw7MBs4L5QwsyppZgR1BXBqRGzMOxh7feoXSixdurSPnmZmxddngpJ0TET8EugEpku9V3tHxA0NT7SW2W233XodO0GZWZn1N4J6N/BLsrmnekHvhQ9WALvuumuv4+eff75FkZiZvX59JqiIuDA9vSginqh9TdLkXKOyQfEIysyqpJnvQf20Qdv12zoQe/08gjKzKulvDurNwIHAzpI+WPPSTsBAvqhrQ6R+BOUEZWZl1t8c1P7A+4DR9J6Hegn4VJ5B2eDUj6Bc4jOzMutvDuom4CZJR0bEb4cwJhukRiW+iKB+BaaZWRk0Mwf1t5JG9xykPfG8q0MBjRw5kpEjR24+fvXVV1m1alULIzIzG7xmEtRBEfFiz0G6aeDB+YVkr4fLfGZWFc0kqO1qb7GebsHe1B5+NvS8UMLMqqKZRPN14DeSrif7gu6HgItzjcoGzSMoM6uKZu6oe5Wk+WSbxQr4YEQ8mHtkNigeQZlZVTRT4gPYBXg5Iv4ZWO6dJIrLCcrMqqKZ221cCPw9cF5qegPwf/IMygbPJT4zq4pmRlB/BUwDXgaIiGeBUXkGZYPnEZSZVUUzCWp9utttAEjaMd+Q7PXwCMrMqqKZBHWdpH8FRkv6FPAL4Pv5hmWD5RGUmVVFM6v4vibpOGAN2f58F0TEvNwjs0GpH0EtW7aMTZs2sd12za6HMTMrhqa+cJsSkpNSCYwYMYKdd96Z1atXA7Bx40ZeeOEFxo8f3+LIzMwGps//Vku6M/18SdKaBo8nJH166EK1ZrnMZ2ZV0GeCioh3pp+jImKn+gfQBXxmqAK15nmhhJlVQVMlPkmHAO8kW8l3Z0T8LiJekHR0nsHZ4HgEZWZV0MwXdS8AfgSMBcYBV0r6B4CIeC7f8GwwfOt3M6uCZkZQpwAHR8SfASRdAiwEvpxnYDZ49SMol/jMrIyaWXv8JDCi5nh74LFcorFtwiU+M6uC/lbx/bOk7wDrgMWSrpR0BfAA8KdmLi5ptKTrJT0s6SFJR0qaKukuSYskzZd0WOorSd+RtETS/Wneq+c6MyU9mh4zX98fufq8SMLMqqC/Et/89HMBcGNN+68GcP1vA7dGxMmShgOdwHXAFyPi55LeC1wKHA2cBExJj8OBy4DD0w0SLyRbNRjAAklz0519rQGPoMysCvpMUBHxIwBJI4A3kSWHx3rmorZG0k7AUcDH0/XWA+slBbBT6rYz8Gx6Ph24Ku37d1cafe1OlrzmRcTKdN15wInANc3/MduLF0mYWRX0maAkdQBfAT4BPEVWDpyYynznR8SrW7n2vsBy4ApJbyUbiX0GOBu4TdLX0jXfnvrvCTxTc353auurvT7eWcAsgL333nsroVXbhAkTeh2vWLGCDRs20NHR1LcKzMwKob9FErPJblQ4OSLeFhEHA28ERgNfa+LaHcAhwGXp3JeBc4G/Az4bEXsBnwV+mPqrwTWin/beDRGXR0RXRHS1+7Y+w4cPZ+zYsZuPI4Lly5e3MCIzs4HrL0G9D/hURLzU0xARa8gSzHubuHY30B0Rd6fj68kS1kzghtT2E+Cwmv571Zw/kaz811e79cMLJcys7PpLUJHmg+obN9JgBNOg3/PAM5L2T03HAg+SJZd3p7ZjgEfT87nAx9JqviOA1emLwLcBx0saI2kMcHxqs354oYSZlV1/kxIPSvpYRFxV2yjpo8DDTV7/TODqtILvceA04Cbg22mO68+keSPgFrKR2RJgbepLRKyU9CXg3tTvop4FE9Y3L5Qws7LrL0GdDtwg6RNkCxwCOBTYgew28FsVEYvIlofXuhN4W4O+kd6z0XXmAHOaeU/LeDcJMyu7/paZ/5Hse0jHAAeSLVb4eUTcPlTB2eC5xGdmZdfMHXV/CfxyCGKxbciLJMys7Hwf8IryCMrMys4JqqK8SMLMys4JqqK8SMLMys4JqqLGjRuHtGUTjpUrV7J+/foWRmRmNjBOUBXV0dFB/ZZPy5Yta1E0ZmYD5wRVYV4oYWZl5gRVYV4oYWZl5gRVYV4oYWZl5gRVYR5BmVmZOUFVmOegzKzMnKAqrP7Our5poZmViRNUhdUvM3eCMrMycYKqMCcoMyszJ6gKc4IyszJzgqqwcePG9TpesWIF2X0hzcyKzwmqwjo7O+ns7Nx8vGHDBlavXt3CiMzMmucEVXEu85lZWTlBVZwTlJmVlRNUxdUnqBUrVrQoEjOzgXGCqrj6hRIeQZlZWThBVZxLfGZWVk5QFecSn5mVlRNUxXkEZWZl5QRVcZ6DMrOycoKqOI+gzKysck1QkkZLul7Sw5IeknRkaj9T0iOSFku6tKb/eZKWpNdOqGk/MbUtkXRunjFXjeegzKysOnK+/reBWyPiZEnDgU5J7wGmAwdFxDpJEwAkHQDMAA4E9gB+IWm/dJ3vAscB3cC9kuZGxIM5x14JHkGZWVnlNoKStBNwFPBDgIhYHxEvAn8HXBIR61L7snTKdODaiFgXEU8AS4DD0mNJRDweEeuBa1Nfa8LOO+9MR8eW/4e8/PLLvPLKKy2MyMysOXmW+PYFlgNXSPqdpB9I2hHYD3iXpLsl/VrSoan/nsAzNed3p7a+2nuRNEvSfEnzPUrYQpIXSphZKeWZoDqAQ4DLIuJg4GXg3NQ+BjgC+BxwnSQBanCN6Ke9d0PE5RHRFRFd9WWtdud5KDMrozwTVDfQHRF3p+PryRJWN3BDZO4BNgHjUvteNedPBJ7tp92a5HkoMyuj3BJURDwPPCNp/9R0LPAg8DPgGIC0CGI4sAKYC8yQtL2kycAU4B7gXmCKpMlpocWM1Nea5BKfmZVR3qv4zgSuTonlceA0slLfHEkPAOuBmZHd5nWxpOvIktgG4PSI2Agg6QzgNmAYMCciFuccd6V4BGVmZZRrgoqIRUBXg5c+2kf/i4GLG7TfAtyybaNrH56DMrMy8k4SbcAjKDMrIyeoNuA5KDMrIyeoNuARlJmVkRNUG/AclJmVkRNUG/AIyszKyAmqDYwdO7bX8cqVK9mwYUOLojEza44TVBvo6OhgzJgxvdpWrlzZomjMzJrjBNUmXOYzs7JxgmoTTlBmVjZOUG3CCcrMysYJqk3Uf1nXS83NrOicoNqER1BmVjZOUG3CCcrMysYJqk04QZlZ2ThBtQnPQZlZ2ThBtQmPoMysbJyg2oQTlJmVjRNUm2hU4ouIFkVjZrZ1TlBtorOzk87Ozs3Hr776KmvWrGlhRGZm/XOCaiMu85lZmThBtREnKDMrEyeoNuI765pZmThBtZH6hRIeQZlZkTlBtRGX+MysTJyg2ogTlJmViRNUG/EclJmVSa4JStJoSddLeljSQ5KOrHntHEkhaVw6lqTvSFoi6X5Jh9T0nSnp0fSYmWfMVeY5KDMrk46cr/9t4NaIOFnScKATQNJewHHA0zV9TwKmpMfhwGXA4ZJ2AS4EuoAAFkiaGxGrco69cupHUHfeeSfHHHNMi6IxszI7/PDD+epXv5rre+SWoCTtBBwFfBwgItYD69PL3wQ+D9xUc8p04KrI9t+5K42+dgeOBuZFxMp03XnAicA1ecVeVfUJas2aNdxxxx0tisbMymyHHXbI/T3yLPHtCywHrpD0O0k/kLSjpGnAHyPivrr+ewLP1Bx3p7a+2m2A9thjD97whje0Ogwzs6bkmaA6gEOAyyLiYOBl4AvA+cAFDfqrQVv00977ZGmWpPmS5ntupbEdd9yRs846q9VhmJk1Jc85qG6gOyLuTsfXkyWoycB9kgAmAgslHZb671Vz/kTg2dR+dF37r+rfLCIuBy4H6Orq8jbdfZg9ezazZs2iu7u71aGYWYmNHTs29/fILUFFxPOSnpG0f0Q8AhwLLIyIY3v6SHoS6IqIFZLmAmdIupZskcTqiHhO0m3AVySNSacdD5yXV9xVJ4n99tuP/fbbr9WhmJn1K+9VfGcCV6cVfI8Dp/XT9xbgvcASYG1P34hYKelLwL2p30U9CybMzKy6VMWb1nV1dcX8+fNbHYaZmTUgaUFEdG2tn3eSMDOzQnKCMjOzQnKCMjOzQnKCMjOzQnKCMjOzQqrkKj5Jy4GnBnjaOKAs959wrPlwrPlwrPkpU7y1se4TEeP76wwVTVCDIWl+M8sei8Cx5sOx5sOx5qdM8Q4mVpf4zMyskJygzMyskJygtri81QEMgGPNh2PNh2PNT5niHXCsnoMyM7NC8gjKzMwKyQnKzMwKqe0TlKQTJT0iaYmkc1sdTy1JcyQtk/RATdsukuZJejT9HNPfNYaKpL0k3SHpIUmLJX0mtRcuXkkjJN0j6b4U6xdT+2RJd6dYf5xuE1MIkoZJ+p2km9NxkWN9UtLvJS2SND+1Fe5zACBptKTrJT2cPrtHFjFWSfun32fPY42ks4sYK4Ckz6a/Ww9Iuib9nRvwZ7atE5SkYcB3gZOAA4BTJB3Q2qh6uRI4sa7tXOD2iJgC3J6Oi2AD8D8j4i+AI4DT0++yiPGuA46JiLcCU4ETJR0B/BPwzRTrKuCTLYyx3meAh2qOixwrwHsiYmrN916K+DkA+DZwa0S8GXgr2e+4cLFGxCPp9zkVeBvZPfNupICxStoTOIvsZrRvAYYBMxjMZzYi2vYBHAncVnN8HnBeq+Oqi3ES8EDN8SPA7un57sAjrY6xj7hvAo4rerxAJ7CQ7C7OK4CORp+NFsc4kewfn2OAmwEVNdYUz5PAuLq2wn0OgJ2AJ0iLxYoca118xwP/VdRYgT2BZ4BdyG6KezNwwmA+s209gmLLL7JHd2orsl0j4jmA9HNCi+N5DUmTgIOBuylovKlktghYBswDHgNejIgNqUuRPgvfAj4PbErHYylurAAB/IekBZJmpbYifg72BZYDV6Ty6Q8k7UgxY601A7gmPS9crBHxR+BrwNPAc8BqYAGD+My2e4JSgzavu38dJI0EfgqcHRFrWh1PXyJiY2TlkonAYcBfNOo2tFG9lqT3AcsiYkFtc4OuLY+1xjsi4hCy0vnpko5qdUB96AAOAS6LiIOBlylAiaw/ad5mGvCTVsfSlzQPNh2YDOwB7Ej2Wai31c9suyeobmCvmuOJwLMtiqVZSyXtDpB+LmtxPJtJegNZcro6Im5IzYWNFyAiXgR+RTZvNlpSR3qpKJ+FdwDTJD0JXEtW5vsWxYwVgIh4Nv1cRjZPchjF/Bx0A90RcXc6vp4sYRUx1h4nAQsjYmk6LmKsfwk8ERHLI+JV4Abg7QziM9vuCepeYEpaXTKcbOg8t8Uxbc1cYGZ6PpNsrqflJAn4IfBQRHyj5qXCxStpvKTR6fkOZH+hHgLuAE5O3QoRa0ScFxETI2IS2efzlxHxEQoYK4CkHSWN6nlONl/yAAX8HETE88AzkvZPTccCD1LAWGucwpbyHhQz1qeBIyR1pn8Xen6vA//MtnpCrdUP4L3AH8jmIM5vdTx1sV1DVsN9lex/e58km3+4HXg0/dyl1XGmWN9JNmS/H1iUHu8tYrzAQcDvUqwPABek9n2Be4AlZCWU7Vsda13cRwM3FznWFNd96bG45+9UET8HKa6pwPz0WfgZMKbAsXYCLwA717QVNdYvAg+nv1//Bmw/mM+stzoyM7NCavcSn5mZFZQTlJmZFZITlJmZFZITlJmZFZITlJmZFZITlFkdSV+VdLSkD6ifHe4lfSzt1rxY0oOSztnKdb/Q00fSlZJO7q9/zXlHS3r7AOKfJKlb0nZ17YskHdYgjifSawslHdmgfZGk36T2j0tantoelvTZZuMyGygnKLPXOpxsH8F3A//ZqIOkk4CzgeMj4kCyHQhW5xTP0WTfxG9KRDxJtsfku3raJL0ZGBUR9zQ45XORbft0LvCv9e3pUfv+P0793wGcL6l2NxazbcYJyiyRNFvS/cChwG+B/wFcJumCBt3PA86JLdv6/Dkivp+u8ylJ9yq739RPJXVu5X0vSSOw+yV9re61ScDfAp9No5Z3SdpH0u2p/+2S9m5w2WvIdp7oUbvBaF/+H/CmrfTZLCJeIPvS5e7NnmM2EE5QZklEfI4sKV1JlqTuj4iDIuKiBt3fQrZDcyM3RMShkd1v6iH6ue+NpF2AvwIOjIiDgC/XxfQk8D2y++hMjYj/BP4FuCr1vxr4ToNLXwd8oGbvsw+T7eXXn/cDv685nl1T4ru6Qex7AyPIdmEw2+Y6tt7FrK0cTLZN05vJ9g85lfXIAAABvElEQVQbjLdI+jIwGhgJ3NZP3zXAn4EfSPp3snvnbM2RwAfT838DLq3vEBHPS1oMHCtpKfBqRDxQ3y+ZLekfyG49UZtMPxcR1zfo/2FJ7wH2Bz4VEX9uImazAXOCMgMkTSUbOU0ku7FaZ9asRcCREfFK3SmLye5s+ssGl7sS+EBE3Cfp42RzSA1FxAZJh5FtqDkDOINsx/KB6Gu/sp4y31L6L+/1lYj68uOIOCMtqPh3ST+PbONVs23KJT4zICIWpYn/PwAHkCWeE1JZrT45AXwVuFTSbgCStpd0VnptFPBcuv3IR/p733T/rJ0j4hayRRdTG3R7KV2zx2/YMr/0EeDOPi7/U7INe5sp7w1YRPyWbAT3mW19bTPwCMpsM0njgVURsUnSmyOizxJfRNwiaVfgF+mWAgHMSS//I9kqwKfI5nRGNb4KpNdukjSC7EaEjZZt/1/geknTgTOBs4A5kj5HVpY7rY8YX5R0F9ldV5/oJ4a+9JT+ehzWoM8/AQslfSUiXhrEe5j1ybuZm5lZIbnEZ2ZmheQEZWZmheQEZWZmheQEZWZmheQEZWZmheQEZWZmheQEZWZmhfT/AaXECOQQUpcpAAAAAElFTkSuQmCC\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": [
    "**Results of partial lock-downs**: to understand what would have happened had the third lockdown not taken place, we perform hypothetical experiments with just the first and then just the first two lock-downs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get hold of outcomes with just the first lockdown\n",
    "pop.reset()\n",
    "np.random.seed(0)\n",
    "\n",
    "LKP = np.zeros((T,))\n",
    "LKP[I0:I0+P0] = L0\n",
    "\n",
    "(statsLKP1, tInfect, tQuarantine, tRecovery, tExpiry) = pop.simulate( T = T, LKP = LKP )\n",
    "\n",
    "# Get hold of outcomes with just the first two lockdowns\n",
    "pop.reset()\n",
    "np.random.seed(0)\n",
    "\n",
    "LKP = np.zeros((T,))\n",
    "LKP[I0:I0+P0] = L0\n",
    "LKP[I1:I1+P1] = L1\n",
    "\n",
    "(statsLKP2, tInfect, tQuarantine, tRecovery, tExpiry) = pop.simulate( T = T, LKP = LKP )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13977 Infected\n",
      "1999 Expired\n",
      "Peak 6052\n",
      "Recovery time 82.791667\n",
      "Expiry time 22.046523\n",
      "Quarantine time 11.152221\n",
      "13303 Quarantined\n",
      "Objective acheived is 6052.000000 (fEpi) + 350.000000 (fEco) = 6402.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(d) in the accompanying paper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXlYVNX7wD+HRYZlQAU0S3M3V0RFzQ1xSUvLzDQzzYVcv5o/00pLS8ylzbKsrKxcMkvT1FIrzRR3U1TENfcFRUNRBNnh/P64wzAsAwMMMOD5PM99nHvW91yH+8455z3vK6SUKBQKhUJha9iVtAAKhUKhUOSEUlAKhUKhsEmUglIoFAqFTaIUlEKhUChsEqWgFAqFQmGTKAWlUCgUCpvEoaQFUCgUCoVtI4S4CMQAqUCKlNKvOPpVCkqhUCgUltBJSnmzODtUS3wKhUKhsEmE8iRhHjs7O+ns7FzSYliNtLQ07OzK1m+SsjYmNR7bp6yNKS4uTgKHTJIWSikXmpYRQlwAbgMS+DprflGhFFQu6HQ6mZCQUNJiWI3g4GACAgJKWgyrUtbGpMZj+5S1MQkh4qSUrnmUeVBKeU0IUQn4C3hZSrmjqGUrOz8DFAqFQlEkSCmvGf79D1gLtCqOfpWCUigUCoVZhBCuQgh9+megG3CsOPpWVnwKhUKhyI3KwFohBGg640cp5Z/F0bFSUAqFGZKTkwkPD6c49yE9PDw4efJksfVX1JS18UDpHZNOp6Nq1ao4Ojrmq56U8jzQtGikyh2loBQKM4SHh6PX66lRowaGX49FTkxMDHq9vlj6Kg7K2nigdI5JSsmtW7cIDw+nZs2aJS2Oxag9KIXCDAkJCXh6ehabclIoigohBJ6ensW6GmANlIJSKHJBKSdFWaE0fpeVglIUP5f3wU8D4JsuuZe7fQlSU4pHJoVCYXMoBaUoPlJTYNNUWNQd/v0droZAXFTOZVOS4Ic+8HF92DgJ/jtVvLLaCPb29vj6+hqv9957r6RFMhIQEEBISEih2ggODubJJ5/Ms9z8+fNp0KABAwcOzHcfn3zyCXFxccb7Hj16cOfOnXy3YymHDx9m+PDhACxZsgRvb+9M/4cnTpwoVPtBQUHMnTsXgFdffZWtW7cWWmZbRRlJKIqP6MuwP4uHlKuHoG7X7GX3fw23zmqfD3wLB5fCS5vgoRZFL6cN4ezsTGhoaEmLUWAaN27MpUuXCt3OggUL+OOPPwq0wf/JJ58waNAgXFxcAPj9998LLU9uzJkzh2nTphnv+/fvz+eff25R3eDgYJYsWcKSJUssKv/yyy8zYsQIOnfuXBBRbR6loBTFR8Va0OFVCJ6TkXbtcHYFFfsfbP8gc1rdblDFt+hlzIkgjyJuPzpfxaOjo2nVqhW//fYbjzzyCAMGDKBz586MGDECNzc3Ro0axbZt26hQoQIrVqzA29ub0NBQRo8eTVxcHLVr12bRokVUqFCB+fPn89VXX+Hg4EDDhg1ZsWIF9+7d4+WXX+bo0aOkpKQQFBTE008/TXx8PMOGDePEiRM0aNCA+Ph46z6GoCAuX77M+fPnuXz5MhMmTGD8+PGMHj2a8+fP06tXLwIDAxk5cmSO8qWmpjJ58mQ2bdqEEIIRI0YgpeTatWt06tQJLy8vtm3bRo0aNQgJCcHLy4uPP/6YRYsWATB8+HAmTJjAxYsXefLJJzl2TDuLOnfuXGJjYwkKCmL+/PksWLCAcuXKGZ+XKTExMYSFhdG0afFYZVevXp1bt25x/fp1HnjggWLpszgp9iU+IcQ4IUSIECJRCLHETJnpQggphOhqkuYkhFgkhLgrhLguhJiYpU4XIcQpIUScEGKbEKK6pXUVxUj7CeD1COgfhIG/QIdJORQS0LgPOLlrtw/4wLPfgJ19sYpqC8THx2daHlq5ciUeHh58/vnnDB06lBUrVnD79m1GjBgBwL1792jevDmHDh2iY8eOzJgxA4DBgwfz/vvvExYWRpMmTYzp7733HocPHyYsLIyvvvoKgNmzZ9O5c2cOHDjAtm3beO2117h37x5ffvklLi4uhIWFMXXqVA4ePGj18Z46dYpNmzaxf/9+ZsyYQXJyMl999RUPPvgg27Zt45VXXjEr38KFC7lw4YJxPAMHDmT8+PHGutu2bcvU18GDB1m8eDH//PMP+/bt45tvvuHw4cO5yvfee++xa9euTM/LlJCQEBo3bpwpbeXKlZn+D62t2Js3b87u3but2qatUBIzqGvALKA7kM1VuBCiNtAXiMiSFQTUBaoDDwDbhBAnpJR/CiG8gDXAcGA9MBNYCTyaV12rjkyRNw5O8Oy34FYJ9GZ+8bl5w1OfwhMfQOiPUK87lMvVl2WZxdwS32OPPcaqVasYO3YsR44cMabb2dnRv39/AAYNGkSfPn2Ijo7mzp07dOzYEYAhQ4bQr18/AHx8fBg4cCC9e/emd+/eAGzevJnffvvNuM+RkJDA5cuX2bFjB+PHjzfW8/HxyVHmsWPHGl+YERER+PpqM99+/foxderUXMfbs2dPnJyccHJyolKlSty4cYOqVatmKmNOvi1btjB69GgcHLTXWsWKFXPta9euXTzzzDO4umrfrT59+rBz50569eplto6Pjw/Dhw+nb9++xudlSkREBN7e3pnSLFnia926NYmJicTGxhIVFWV8Zu+//z7du3fPtW6lSpW4du1armVKK8WuoKSUawCEEH5A1RyKfA5MBhZkSR8MDJNS3gZuCyG+AYYCfwJ9gONSylWGtoOAm0KI+lLKU3nUVRQ3VXJ+sWXDwQn8huVaxC61dJ3rsBZpaWmcPHkSZ2dnoqKisr3E08nLtHjjxo3s2LGD3377jZkzZ3L8+HGklPzyyy888sgj+W4P4IsvvjB+rl69er720JycnIyf7e3tSUnJbsVpTj4pZb5Mqc1FcnBwcCAtLc14b3p2aOPGjfz5559s2bLF+LzSFSJoPygKctbon3/+AfK/B5UuX1kKC2SKTe1BCSH6AUlSyt9Nv2hCiArAg8ARk+JHgPSfMI1M86SU94QQ54BGQogbedTNKsNIYCRoX9Tg4OBCjsp2iI2NLbbxOCbdIdnRA4rw7IUu/gZ+h9/gZORgbjwQYPX2PTw8iImJgUnhVm87EzExxo+pqalan5myY7LW4LPPPqNOnTpMmzaNoUOHsmXLFhwdHUlLS2PZsmX07duXxYsX06pVK+zs7PDw8GDTpk20bduWb7/9ljZt2hAdHc2VK1fw8/OjadOmLF++nIiICDp16sRHH33E3LlzEUJw5MgRmjZtSuvWrVmyZAl+fn6cOHGCsLAw7t27l6N86Ugpc82Pi4sjJSWFmJgYEhMTcXR0NJZPS0sjNjaWmJgYpJTExsbi5ORkVj5/f38+//xzWrRogYODA1FRUVSsWBFXV1ciIiKMyi+9rRYtWjBmzBjGjh1rVHoLFy7ExcWFGzducPHiRdzc3Pj111/p2rWr8Xm1a9eONm3aGJ9X+fLljeN5+OGH+ffff41jSEhIICkpKddnkPV5JCcn51o+63M6ceIEPXr0sKiPhISEUvVOsxkFJYRwA+agecrNipvhX9Pd5GhAb5IfmaVOen5edTNhCMS1ELR4UGUp7kuxxbFJS4Uv24KuPDw2Ax5+NO86+SU6HBaPh6RbNDizgAbtekK1llbt4uTJk8Xu0iarG534+Hg6dOhgvH/88ccJDAxk2bJl7N+/H71ez6ZNm/j000+ZMWMGrq6unD9/noCAADw8PFi5ciV6vZ5ly5YZjSRq1arF4sWLcXFxYfTo0URHRyOlZOLEiVSrVo2ZM2cyYcIE2rVrh5SSGjVqsGHDBiZMmMCwYcNo164dvr6+tGrVCldX11yfkRAi13wXFxccHBzQ6/XGpb308nZ2dri5uaHX6xFCGD+bk2/cuHFcvnyZdu3a4ejoyIgRIxg3bhyjR4+mX79+VKlShW3bthnb6tChA4GBgXTpop3HGzlyJO3btwdg+vTpdO3alZo1a9KoUSOcnJyMz+v27dsIIYzPy5QWLVoQGxsLgF6vR6fTsXbtWvbv328ss2DBAtq2bWv2eTg6Oub6zEyfU3JyMhcvXqRjx46ZZnLm0Ol0NGvWLM9ytkKJBSwUQswCqkophxruPwKipZTvGO4vAsOllFsMM6gooLIhHglCiGeBICllEyHEp4CjlPJ/Ju0fRdt72ppb3dxkVAELC8ihZfDbuIz7+k/Cc99bz8ghNQW+7gD/mZwnca0EI7eBR85LXQXh5MmTNGjQwGrtWUJh/by5ubkZX5C2QGn0W5cXeY1p3rx56PV641moomTt2rUcOnSImTNnWlQ+p++0JQELSwpbOqjbBRhvsLK7DlQDfhZCTDbsHUWQ2aNuU+C44fNx0zxDzJLaaPtSedVVWJPkeNg2J3Oak7t1LfDsHaB9FkPMSvXB0cV6fSgUBWTMmDGZ9tKKkpSUFCZNyskStmxQEmbmDkIIHWAP2AshdEIIBzQF1RjwNVzXgFFA+o7r98A0IUQFIUR9YASwxJC3FmgshHjW0PbbQJjBQCKvugprcucKuHpl3Ns7Qac3rd+PT78MJdVyOAxaAy65W23dD9jS7Ol+RafT8eKLLxZLX/369cu0B1bWKIk9qGnAdJP7QcAMKWWQaSEhRCpwW0qZ/hc3HfgSuATEA++nm4lLKSMNy3afAz8A/wDPmzRntq7CynjXg1E7NH97+7+G8tWhfLW86xWEzm8RFu2KT89Xi6Z9hUJRopSEmXkQ2t5QXuVqZLlPBAINV07ltwD1zeTlWldhZYSA6m20qyixsyPK069o+1AoFCWGLe1BKRQKhUJhRCkoRdmlhCxUFQqFdVAKSlG2kBIu7dXiTW2bXdLSFBprhdt4++232bJlS6FkcXNzyzE9Pj6ejh07kpqaavRFl5SUBMD58+epVasWd+/e5ejRowwdOrRQMijuLyzagzL4unORUl42SRuFZnW3SUq5oYjkU5QWEu6Czr1kZbh1DtaM1OJMAVzcDY/+r1Rb91kr3MY777yTY3pqair29oU7ArBo0SL69OmDvb09fn5++Pv7M3fuXN58800mTZrE7NmzcXd3p0mTJoSHh3P58mUefvjhQvWpuD+w1EhiERAO/A9ACPEWMAO4DfxPCPGClHJl0YiosHnSUuHTppq1Xr3HNeeuDzYvUjdHOaJ/AKLOZdwnRsPuTzVvFoXgZP2iPazb4NTJfJUvSLiNoUOH8uSTT9K3b19q1KhBYGAgmzdvZty4cbRs2ZKxY8cSGRmJi4sL33zzDfXr1+fChQu88MILpKSk8Pjjj5uVZ/ny5fz444/G+zlz5tC8eXMcHBxISUlhwIABxrynnnqKFStW8Prrr+f/QSnuOyxd4vMD/ja5Hw3MkVJ6op1TUuEr7meu7If4KIg4Atvfh2V9NKVV3JRzhbbjM6dd2AEmjj9LG9YKt5EVnU7Hrl27eP755xk5ciSfffYZBw8eZO7cufzvf5pDlv/7v/9jzJgxHDhwwGysoaSkJM6fP0+NGjWMaeXLl2fy5Mm88cYbfPzxx5nK+/n5sXPnTis8GcX9gKUzqIrADQAhRGO0kBVLDXnr0LyFK+5XTmc5Ulb3Mc3bQ0nQerQWgffBZtB6FNToUPwzOStijXAbOZFeJjY2lj179hjDb4DmjBRg9+7d/PLLLwC8+OKLTJ48OVs7N2/ezPGg6B9//EHlypU5deoUzZs3N6aX5dAQCutj6VvkFhmhMToD16SUZwz3jihji/ubmOuAAAxWc/XMLwcVOeVc4OWD4Fg2ww+kU9hwG+kxkNLS0ihfvrzZfa68wlfkFF5iw4YNREdHs2nTJp5++ml69+5tDLdelkNDKKyPpQpqCxBkMJaYhDZrSqc+mocGxf1Kn6+h2yw4s1mbTdXpUrLyWFk55XePqDiYN28eDRo0YM6cOQQGBrJ3715juI3Vq1fz/PPP8+OPPxq9c5vD3d2dmjVrsmrVKvr164eU0hiyvF27dqxYsYJBgwaxfPnyHOtXqFCB1NRUEhIS0Ol0xMfHM2nSJNauXUvDhg3p2bMns2fPZvZszaLy9OnT2SLOKhTmsHTm8zpwBXgXOIdmIJHOQGCXleVSlDbcvKHZQOi/DJwrlLQ0ZYase1BTpkzh9OnTfPvtt3z00Ud06NABf39/Zs2aBWgzo+PHj9OiRQu2bt3K22+/nWcfy5cv57vvvqNp06Y0atSIX3/9FYBPP/2UL774gpYtWxIdHW22frdu3di1S3sFzJw5k969e9OwYUMA3njjDVasWMGZM9qCy7Zt2+jZs2ehnoni/qHQ4TaEEO5AgpQyyToi2Q4q3IbtU5RjUuE2LOPw4cN8/PHHLFu2LFue6XgSExPp2LEju3btsih2ka1SmkOI3HfhNqSUd8uiclKUMaSEGyrCSlHQrFkzOnXqRGpq7pably9f5r333ivVyklRvJj9pgghFuWjHSmlfMkK8igU1iX+NoT+BCGLIOo8vHIc3KuUtFRFRkmF2wgMzNsPc926dalbt24xSKMoK+T2U6YzRrOsPFFOzxS2ybI+cO1Qxv2hpRAwpeTkUSgUFmNWQWUNd6FQZGP3p1qAwmqtoVpLLfaTrZ05ajogs4I6uBQ6vFpy57QUCoXFqPNLioJzZCUc+AbWDNdcHZ39O+86xU3T/hmh4O0c4OFHIcG8RZpCobAd8v0zUghRCdBlTTd1JKu4D4i/Df+dMEkQULVFiYljFp0HtH0Z7B2h2WDQVy5piRQKhYVYNIMSQtgJIeYIIW4BEcCFHC7F/cSV/WTaeqzU0HbPP3V6E/xfK5XKSQjBpEmTjPdz584lKCio5AQqAYYOHcrq1atzLZOYmEjXrl2N/grzw8WLFzM5uw0JCWH8+PG51Cg8EyZMYMeOHQAEBAQQEhKSr/rh4eE8/fTT1K1bl1q1ajFu3Diji6qyFNbE0iW+CcBY4CM0nzZzgFloiukcMMLSDoUQ44QQIUKIRCHEEpP0R4UQfwkhooQQkUKIVUKIKib5QgjxvhDiluH6QJj4YRFC+AohDgoh4gz/+lpaV1EAqreFgauh/UR4uA3UCihpicokTk5OrFmzhps3b5a0KBaTkpJS7H0ePnyY5ORkQkNDjX4GLSWrgvLz82P+/PnWFtFIVFQU+/btw9/fv0D1pZT06dOH3r17c+bMGc6cOUN8fLzRQ7xpWJPSjqVLfMOAd4BP0BTTWinlISHELGAzkJ/gLtcMbXQHTH3SVAAWApuAFOBzYDGQ7thtJNAbaIr20/0v4DzwlRCiHPCrQb4FwCjgVyFEXcMZLbN18yG3whQnveYUtu5jJS1JkfPF6K1F2v7YrzqbzXNwcGDkyJHMmzfP6C4onUuXLhEYGEhkZCTe3t4sXryYhx9+mKFDh+Lu7k5ISAjXr1/ngw8+oG/fvtnazhpOY968ecTGxhIcHMzcuXPZsEEL8zZu3Dj8/PwYOnQo77zzDuvXryc+Pp62bdvy9ddfI4QgICCAtm3bsnv3bnr16kW9evWYNWsWCQkJeHt7s3z5cipXrkxQUBCXL1/m/PnzXL58mQkTJuRrtlKjRg2GDBnC+vXrSU5OZtWqVVSsWJFBgwYRGRmJr68vv/zyC3fu3GHixInExsbi5eXFkiVLqFKlCmfPnmX06NFERkZib2/PqlWrmDJlCidPnsTX15chQ4bQrFkz4/ijoqIIDAzk/PnzuLi4sHDhQmrWrElQUBBubm68+uqrADRu3JgNGzbg7e3Nc889R3h4OKmpqbz11lvZFObq1atzDV+SF1u3bkWn0zFs2DBAC2o5b948qlevzuzZs3FzcyuSsCZCCHsgBLgqpXzSag3ngqUzqFpAiJQyFU15OANIKZPRlELehyAMSCnXSCnXoTmgNU3/Q0q5ynDwNw5NQbUzKTIE+EhKGS6lvIo2mxtqyAtAU7afSCkTpZTz0WZ6nS2oq1DYNGPHjmX58uXZ3A2NGzeOwYMHExYWxsCBAzO96CMiIti1axcbNmxgypSczeotCaeRlXHjxnHgwAGOHTtGfHy8UYkB3Llzh+3btzNp0iTat2/Pvn37jCE9PvjgA2O5U6dOsWnTJvbv38+MGTNITk7Oz+PAy8uLQ4cOMWbMGObOnUulSpX49ttv6dChA6GhoTz88MO8/PLLrF69moMHDxIYGMjUqVMBGDhwoNED/J49e6hSpQrvvfeese4rr7ySqa/p06fTrFkzwsLCmDNnDoMH5x644c8//+TBBx/kyJEjHDt2LEdFtHv3blq0KPh+bborK1Pc3d2pUaMGZ8+eBYosrMn/AcXqmNLSGVQ0GYYR14BHgN0mbRRFyFJ/wPTofyPgiMn9EUNael6YzOy3KcyQ/mcedTMhhBiJNuPCwcGB4ODggo/Axkj/dVyWKOiY7FPi8Ig+QZSnn9kyHh4exMTEFEI6yzDtIzU1NVufQgj69+/Phx9+iLOzM4mJicTExLBnzx6WLl1KTEwMvXv35rXXXiMmJobk5GS6d+/OvXv3qFatGjdu3MhxHLt27WLJkiXG+pMnTyYmJoa4uDhSUlKMdZKSkkhISCAmJobff/+dTz75hPj4eG7fvk2dOnUICAggNTWVp556yljn33//5c033+T69eskJydTvXp1YmJijHtFSUlJODk54eXlxblz53jooYfMPp/k5GTi4+OJiYlBSkm3bt2IiYmhfv36rFq1KpvMJ06c4NixY3Tp0sX4TCtXrsy1a9cIDw+na9eumZ5H1vGa3u/YsYNly5YRExNDy5YtuXnzJlFRUSQmJuLo6Gisk5aWRmxsLDVr1uSvv/7ilVde4fHHH6dt27bZnv2VK1dwcXExpqempnLv3j2Lv2vx8fEkJydnK2/ajouLC+Hh4dnKJCQkFOjvRQhRFegJzKYY4/9ZqqAOAw3Rlt82ATOEEPFos6nZwKFc6uYbIYQP8DbwtEmyG5qiTCcacDPsJWXNS8/X51U3i1JDSrkQbakRnU4nC+rnLTYplhuxNwpUt7JbZdzKuRWobm6Y81t3N/Eukfcird5fQXnI/SF0DtkMRXMk3774LuyAw8vh5G+QHA8TwqB8zivUJ0+eLBafa6Z95OTnTa/XM3nyZJo3b86wYcNwcnJCr9cjhECv1+Po6EhycjJ2dnbG+/LlyxvbkVKi1+uZOnUqGzduBCA0NBQhBO7u7jg4OJD+Z6DX63F3dze2BdrLV6fT4ejoyKRJkwgJCaFatWoEBQUZ27a3t8fb29tYZ8qUKUycOJFOnTpx8OBBgoKC0Ov1ODk54ebmZizn6OiITqfL9Tk7Ojri7OxsHLOnp6dRzvT+XVxccHBwMH5u1KgRe/fuzdTO3bt3jc/MFNO6We+FEJnkFULg4OCAq6srjo6OxvSkpCTc3NyoUaMGhw4d4vfff2fmzJl069Ytm8Pe9OeVXtfe3h5XV9dMcv3zzz+MGjUKgHfeeYdevXoZ81q0aMHGjRszlb979y6RkZE0b94cZ2dno4xZx6rT6WjWrFnWR+wghDC10lhoeA+a8gma0/BidUJoqYL6BG2ZD2A60BxI979/CRhnLYGEEHWAP4D/k1KazlFjAXeTe3cgVkophRBZ89LzY/Kqay25syKlJFUWLKpsEYpltr+CyloUFOn4twTB1YMZ90dXQ4fcfxDmtkdUXFSsWJHnnnuO7777zuhWqG3btqxYsYIXX3yR5cuX5xlawzTsBWA2nEb16tU5ceIEiYmJJCQk8Pfff9O+fXtj3CcvLy9iY2NZvXp1jntboIWlT58VLV26NMcyRcUjjzxCZGQke/fupU2bNiQnJ3P69GkaNWpE1apVWbduHb179yYxMZHU1FT0er3Z2Yu/vz/Lly/nrbfeIjg4GC8vL+NyWvry5qFDh7hwQTNkvnbtmnFPzM3NjSVLlmRrs0GDBpw9ezbXH1atW7c2G6OrS5cuTJkyhe+//57BgweTmprKpEmTGDdunDHWVj7DmqRIKc0uJQghngT+k1IeFEKYF7oIsGgPSkr5l5Tya8Pn60AroB7gC9STUoZZQxghRHW02FMzpZRZXSMfRzNySKcpGUuAxwGfLJZ5PlnyzdVV5Acp4fRmuH2pdIZS93k+8/3xNSUjRwGYNGlSJmu++fPns3jxYnx8fFi2bBmffvppvtozF06jWrVqPPfcc/j4+DBw4EDjL+7y5cszYsQImjRpQu/evWnZsqXZtoOCgujXrx/du3fHy8vLInl69OhhlWi75cqVY/Xq1UyePJmmTZvi6+vLnj17AFi2bBnz58/Hx8eHtm3bcv36dXx8fHBwcKBp06bMmzcv2zhCQkLw8fFhypQpRmX77LPPEhUVha+vL19++SX16tUDNBPvVq1a4evry+zZs5k2bVo2+Xr27FmopXYhBGvXrmX16tXUrVsXT09P7OzsjPtsYPWwJu2AXkKIi8AKoLMQ4gdrNZ4bhQ63ke8OhXBAm7lNR4vSOwJtqbAysAP4Skr5YQ71RqNt0nUlwxLvMylluhXfGeBjNMu8EcBrQF0pZVJudXOTtTDhNu4m3uV67PUC1a3iVgW9k/Vn0uaWw+4k3OG/e/9Z1sjda7DUYMDj6AoPNode+Xsx5kVV96q4pHt/yIN8L/HF3ICP64NMA30VaPQMPPaOdpA3C6Ux3EZhKIpQHaU5NIU5rDGm9u3bs2HDBsqXL19oefbs2cOAAQNYs2YNLVq0yDWsSWHDbRhmUK8WlxWfRUt8Qog8DfallDss7HMamnJKZxBaAESJtow4XQhhzJdSpm/GfG3IP2q4/9aQhkEJ9TakvYdmadLbJAyI2bqKfPKfiRFP8j1IKhnv2QVGXxm6zYYqTbXzW3bK25ei+Pnoo4+4fPmyVRRU27ZtuXQpI6h5WQprYukIgsnbY7m9JQ1JKYOAIDPZM8ykY9gvet1w5ZR/GMjRdjOvuraGtGXn8JFZrEy965eMHIWhzf9KWgKbpKRCddyPtG7dusjaLsqwJlLKYDR9UCxYqqA65ZDmCTwJdMSKRhIKG8ejmua9/MZJSLoLlYp3Cay4kVKinI4oygJFuZ0jhMjPaXYppexiSUGLFJSUcruZrDVCiHnAU2iWdwoDxb23VxjyJWvDp7VLSrh7FZyyGk8WszxFiE6n49aWISxhAAAgAElEQVStW3h6eiolpSjVSCm5desWOp1lxzcKgB2WxwW0+I/JGouUG9EsO9S6yf2EEOBRtaSlKFKqVq1KeHg4kZHFd0YsISGhKF8ixU5ZGw+U3jHpdDqqVi2av1kpZUBRtGsNBfUIUArtjW0XW5lB3FfcPKPFjfLI8Gjg6OhIzZo1i1WM4ODgnA5SllrK2nigbI7JVrHUii8nB1TlgMbAS0DpOUyiUKQTHQ5hP2tnoa4fhXYT4DGzdjoKhcJChBAPAZPQXNZ5Ak9JKY8JISYAe6WU/1jSjqUzqCVm0hOBlWhnjBSK0sWlPfC3iUI6vga6Btle2HqFohQhhGgE7ARSgb1AM7QJDUB1NEcPL1jSlqUKKqd1jgQpZcGczd0H2LSpeBZsTdZik+eRJ8BBBymGw9h3LsPVQ7YZGVihKD18hHYWtTuQACSZ5O0B3re0IUut+C7lXUpRpkmMgR+fgwq1oGIt8KqjWfOVZtJjWp1cD8IOanQAG/JJqFCUUtoDA6SUsYYYUqbcACyL7YJ1jCQUVsbWZjQARF2A2BvadWWvpqhKu4ICaDUSanbUxuJWqaSlUSjKArkZzXkB8ZY2ZFZBCSHSsNyuHSmlRZ4kFKWU2xcy31csXuu2IqOmv3YpFAprsR8tCvv6HPKeIyOWYJ7kNoN6hwwFJdCi5jobOk2fpj2Jpg2/s7TD+4UyZyoedT7zfYVaOZezAmXu2SkU9xczgS1CiM3Aj2h6pKsQ4v+AZ9As+yzCrIIy+MwDQAgxDS3uU3dDOPb0dFe0AIYp+RyAwoawSCE8OhbqP6XNpKLOa+6OFAqFIgtSyu0G592fAIsMye8BF9GceFtkYg6W70GNAsaaKieDIPeEEHOBz9Ai6yqsgE3OIBzKaYYRXnVKWhKFQmHjSCk3AhsNAWgrAbeklP/mtx1LYw14kWHHnpVyaAexFIqyQVJc5qi7CoXCYoQQPumfpZRnpZR7CqKcwHIFFQLMMJwONhXkIbTQGQcK0rlCYTOkpcLZv2HNKJhbF5Y9AymJJS2VQlEaCRVCHBFCTBJCVClMQ5YqqPHAg8A5IUSwEGKlECIYOIdmLKE8SWTBJk3FzWBrspaIPDIN1oyAsBVaEMaEaDjzV/HLoVCUfgag2Sy8C1wWQmwSQgwUQlgWJtsESw/qHjasJU4EHgWaABHAXGCelPJWfjtWKEqSu3chIcE0xREeGgyhyzKSdv0JngWPbO3qql0Kxf2ElHIlsFII4YWmrAYBy4BYIcRaYJmUcoslbVl8UNeghKYWQF5FPrG1GQ1n/gJ7R8203OMhsCv9R97i4jQllYmqfWDfb9pn54pgVwvuFLwPOzuloBT3L1LKm2gGdJ8JIeoCLwLDgYFYqHssXeKzGkKIcUKIECFEohBiSZa8LkKIU0KIOCHENiFEdZM8JyHEIiHEXSHEdSHERGvVVeTB3s9h40T4oTd82Q5unStpiYqGyo2hxTB45hsI3AQd1NdEoSgshqW9VoarEvk4lpSbJ4mtwP+klKcsCOdrcQhf4BowC82RoLNJf15oYTuGox0GnonmKf1RQ5EgoC6aN9wHgG1CiBNSyj8LU9dCmfONTZqKmyFXWZMTtLAU6aQlgb5Q+56Fk8dqfeSQaGcHbV8u8r4VirKO0EJQd0WbNfUG3NA8m49FezdbRG7TLNOYA3mF87U4PoGUcg2AEMIPMA3v2Ac4LqVcZcgPAm4KIepLKU8Bg4FhUsrbwG0hxDfAUODPQtZV5EbUBTL917tXhXL53uu8LylFv1EUCqshhPgQLZxGFTRDuo/Q9p3O51oxB3LzJNHJ5HNA/sXMN42AIyZ93hNCnAMaCSFuoFkRHjEpfwRNMxe2biaEECOBkQAODg4EBwcXaDDRydHcTc66yWEZ7o7ueDh6FKhubsTGxuY4njtJd4hJicmxji7hBpUrB+Acdw2XuKvE2lXk9D8WHwQvEBXKVcDNwc2isubGlBe3bpUjLq5o99L0+hTKl0/OV52CjsdWKWvjgbI5JisTCPyMppT2FKYhSyPqehaDpZ4bEJklLRrQG/LS77PmFbZuJqSUC4GFADqdTgYEBFg8AFNuxd3iVnzBHpmnsyeeLtY/+xwcHExO44m8F8nthNu51Oyl/SMlFVISae2os7psplR2rYyHzjIFbW5MeRERATE562SrUaECeHvnr05Bx2OrlLXxQNkck5WpIqVMyrtY3lhqJBEhhFgnhHhWCGHOo0RhiQXcs6S5AzGGPLLkp+cVtm6RYHOWeNZACChi5VRc5Gv5LTocUqzy96ZQlHnSlZMQwsdgFDddCPGAIa2OECLHyUFOWKqg3gJqAauA60KIr4QQ7fIreB4cB5qm3xgc0dZG21u6jXbuqqlJ+aaGOoWte99ja8rUJuS5cwVCFsOKgfB9LwgPyXcTag9KcT9isJpeBRwG5gNvo22zAHxAPo4rWaSgpJTvSyl9gObAYuApYIcQ4pwQIshg426p8A5CCB1gD9gLIXRCCAdgLdDYMEvTGQYVZjByAPgemCaEqCCEqA+MAJYY8gpT1+awiRf0/c7BJbD3M4g8qd2fVV4lFAoLmU2GBV9lMhvR/YFmwW0R+ToHJaUMlVJOAqoBT6AFnpqEFn/eUqahxZCagnbCOB6YJqWMBJ5FG9xtoDXwvEm96WgWIZeA7cCH6WbihamruD/Jc3ZTt2vm+/PBkJo/gweF4j5lANo7/UcgKkveBaCGpQ0VKOS7lDJNCHEPTbkkAxbbHRviTAWZydsC1DeTl4hmHRJo7boKM6SmgH2BviKln4f8wMkDEqMBARVrQ1wU6Ctb3IRa4lPcp3hiftJiBzhZ2lC+3j4m7ioGomnBq8DXaH6WFCaUiYO6u+ZpS1uetaFiHajfEyrl+BugeOQpTuwdwe8lLQ5WrU7glk9zPIXi/uUC0AbIycFDK8Di0BuWmpmPRVNMLYF7wC9o+zjbpE28TRRFQtRZiLupXVf+gQd9i0VB2QzNB5W0BApFaeR74E0hxEU0Dz8AUgjRCXgFMytoOWHpDOoTYAuaklorpYy3WFRFvrEZnZ/V555n2YmmayuPWKEog3yAZim9DPjWkLYL0AErpJSfWdqQpQqqmpTyer5EVJRuEqIhwcSVt1058KhqvrxCoVAAUspU4HkhxBdoFnuVgFvAn1LK7flpy9J4UEo53W/oPGD0Hoi+ArcvagqrDITZKE7ULE1xPyOl3AnsLEwbJeHN/L6gNJ1lMiurQznNQMKztm3IYwukJEFyHDiXL2lJFIoyT7F7M1fkjU2/oMsI+ZrdpCRpRiJnt2jnoeo+Bp2nFZVoCkWpQwiRRu46whQppbRo9c6WvJkrFLZJ+AHY8H8Z9+e2Qccp9+8ZMYUiO+9guYKyGPUXplDkRdWWUE4PSQb/wgm34dpBqNY612pqD0pRFjC4j9uBdsDWAVgtpZxuWsbggMHq5LYH5Z+fhqSUOwovTtnBZkzFLcDWZLU1ebTDugFwar12r6sA94o6+oxCYTMkAp2llLFCCEdglxDiDynlvqLuOLcZVDAZUzZB3tM3ZeJlJUr8BR1zHf79HSrU0C73qtpLugyR70dcvyfYl9N89D3YQi3vKe4bDM4Y0sMWORquYnlJ5fZX1snkc3ngM+AYsAK4gealdgBaNNuxRSWgogS4cRz2fp5xX60N9P6i5OSxBaq10i6FouzhIIQwjSez0BC41YgQwh44CNQBvpBSFm1Y7XTBzGWYHqgSQiwBNksph2cp9r0Q4jugD7C+SCRUFD+3L2a+r1C9RMQo7ZT0RFihsJAUKaVfbgUMh299hRDlgbVCiMZSymNFLZil4TaeBlaayVtpyFeYUKpNxbMpqBrF2n2pfnaFxCkhEo6uhnNb4W6E5lFeobARpJR30LZ/Hi+O/ixdSLdDm9rlFLWtLmr/qVSTTSHU9IdyrpqiunOp2BVUcWCrs5tK/+2CfUsyEvothUa9S0wehUII4Q0kSynvCCGc0YIRvl8cfVuqoDYC7wohbgJrpJSphjXJZ4FZwIaiEvB+pMRnEHUf0y5F3tw6BxVrgch+Vj1HJRh+EFw9zSr9ilGHMic4Vyi8jApF4agCLDW88+2An6WUub7zhRDuQA/gYTQnsaZIKeVMSzq2VEGNR4uiuxJIEULcBioY6u8y5CsU9wf3bsHpP+DURrj5Lzy7SAtFkhvJCbDhFTjyI/j0hz4Ls5dJjMUj+kTmtNxcKt2+pPYHFUWOlDIMaGZpeSFEOzSbBHNfXglYT0FJKW8CHYQQjwGPomnUCGCvIZKtQnH/sOMDLZBjOqc25K6gYiNhxQsQvl+7P7oK/F8j2aMG4eHhJCQkaOlpqaQ8sQ4HOyAtRbui7OBODsFJU5Mg5gaUu6I59rWzTbN3Dw8PTp40F1y1dFJax6TT6ahatSqOjo5F3dUnwEW0mIFHpZRJBW0oX99qKeVf5LwPZTWEEDWABWgRGROB1cAEKWWKEMIX+A5ogBZS+CUpZaihngDeA9ItDb8DJqcHVMytblFQ4meZ8oGtyVoc8hSqi0d6ZlZQZ/6CDq+CY9aVDAOxN+A/k5mRTIPt7xPefCp6vZ4aNWogDEuEMTEx6PX6vIW/dRYqpkf5TQWPB8DV9qL+WjSeUkZpHJOUklu3bhEeHk7NmjWLursGwHNSyoOFbchSKz4jQohKQoiHs16FFcSEBcB/aLM0X6Aj8D8hRDngV+AHtOXFpcCvhnSAkUBvtEBZPsCTwCiDzHnVVSgsp3obzZsEAPZQ1Q8So7MVMyrBBxpD30UgDH9uXvWg2SASEhLw9PQ0KieLSYiGpFiTBJnRtkKRA0IIPD09M2brRctlNLdIhcaib7UQwl0IsVgIEYe2tHchh8ta1ETbhEswxKH6E+0wcADajO8TKWWilHI+moeLzoZ6Q4CPpJThUsqrwEfAUENeXnVtihKd0SQXyxe4dGPvCM1ehNZjYNgG6PkRuFXOvU697tB9juYy6aW/tH8h/8oJtOU8B+eMeyd3cK6Y/3YU9xUF+q4VjBnAFIOhRKGw9GfXF0B/tCWyMUBgDpe1+BQtGqOLEOIh4AkylFSYzPz2DjOkY/j3iEnekSx5udVVpPNDX/iuO6wZCdve1X6tK7LjNxRajchbMZnSejQMWlP4WFJObuD9CHg8rLlf8qiWoxVhcTB06FBWr14NwPDhwzlx4kQeNcwTFBTE3Llz810vODiYJ598ssjrKCzmSTRPQxeEEBuEEN9nuZZa2pCle1DdgdeklMXh72Y72ubaXbTzVUuBdcA0IOvbMhpIXwx2y5IfDbgZ9qay5mWta0QIMRJtuRAHBweCg4MLNIhbibeIS40rUF0Xexc8nTwLVDc3YmNjcxxPZGIkCakJ2KUl0TL2mpYYF4m8eogQXQfS7Ip8UzUTegc95ctZ9hI3N6a8iIjQkZJi/Ze6XVoSaXbayrFOl4a3d6LZsh4eHsTExBjvz9w+Q1pqGna387FcJ4FrOR/or1uhruXtFIDU1FSSk5OJj48nJiaGefPmAWQaU2pqaqb7rPXt7TOOUCYmJuLo6Gi2vDni4uJISUnJV72C1EkntzHZOgkJCQV+p+WD9mjfzLvkPAmweInIUgUlgH8tbbSgCCHsgE3A10BbNMWyCO1QWASQdcroDqR/U2Kz5LsDsVJKKYTImpe1rhGDD6qFADqdTgYEBBRoLNdirhGbaZ/AcvTl9FTRVylQ3dwIDg4mp/GE3w0nLjkO/jsFJh65hPtDtGzT3upy5EV5XXkquVayqKy5MeXFhQuQnJzvarkTEQYb3tQMJuo/gbMzVKtmvvjJkyczbba7JbkRey8WN1c3q4ij1xl+WJgz3gB++OEH5s+fT1JSEq1bt2bBggXY29szZswYDhw4QHx8PH379mXGjBkA1KhRg8DAQDZv3sy4ceNwdHTE2dkZvV5PQEAAc+fOxc/Pj82bNzN9+nTi4+OpW7cuixcvxs3NLVv9559/3iiLk5MTTk5O6PV6QkNDGT16NHFxcdSuXZtFixZRoUIFzp49y+jRo4mMjMTe3p5Vq1bh4uKCg4MDer2eAwcOMHLkSH755Rdq1aqVaax//vknEyZMwMvLi+bNmxvrREVFERgYyPnz53FxcWHhwoX4+PjQpEkTdu7ciYeHB15eXsybN4/Bgwfz/PPPM3z4cMLDw/ntt9+Ii4vj3LlzPPPMM3zwwQekpqby0ksvERISghCCwMBAXnnlFav8nxYWnU5Hs2YWW4wXCCml1awwLP2ptgJ4ylqd5kJFtPNWnxv2im4Bi9EOfB0HfETmhVQfQzqGf5ua5DXNkpdbXZuixA7qxv4HwuQ3Sxn0IFFkxEbCH69rsaL+mgqbppbs8mhyItw6A1HnIS01xyInT55k5cqV7N69m9DQUOzt7Vm+fDkAs2fPJiQkhLCwMLZv305YWJixnk6nY9euXZmUiyk3b95k1qxZbNmyhZ07d+Ln58fHH39scX2AwYMH8/777xMWFkaTJk2MCnLgwIGMHTuWI0eOsGfPHqpUyfght2fPHkaPHs2vv/6aTTklJCQwYsQI1q9fz86dO7l+/boxb/r06TRr1oywsDDmzJnD4MGDAWjXrh27d+/m+PHj1KpVi507dwJw4MABHn30UQBCQ0NZuXIlR48eZeXKlVy5coXQ0FCuXr3KsWPHOHr0KMOGDTM7TkXuWDqD2gx8IoTQA78DUVkLSCm3FlYYKeVNIcQFYIwQYi7aDGoI2n5SMJAKjBdCfIW2DAiQ3u/3wEQhxO9oU8hJaB7YsaCuAqCWP4zZDXevaodAy7mUtESlg7Q0+GMy3PsvI+30H+DTDep2LH55khM1F1XuDwEpEB2e44Hev//+m4MHD9KyZUsA4uPjqVRJm7n+/PPPLFy4kJSUFCIiIjhx4gQ+Pj4A9O/fP9fu9+3bx4kTJ2jXrh1paWmkpKTQpk0bY35e9aOjo7lz5w4dO2rPbsiQIfTr14+YmBiuXr3KM888A2iKLp2TJ08ycuRINm/ezIMPPpitzVOnTlGzZk3q1tWWPQcNGsTChdph6V27dvHLL78A0LlzZ27dukV0dDQdOnRgx44dVK9enTFjxrBw4UKuXr1KhQoVcHPTZrldunTBw8MDgIYNG3Lp0iUaNWrE+fPnefnll+nZsyfdunXLdbxlDSHEPrR3azCwU0oZX9C2LFVQvxr+rUmGZRxoiiA9VpS1/PH1QTvoNRlNqWwDXpFSJgkhegPfop13Ogn0NjkE9jVQCzhquP/WkIYFdRXp2DtmxIEqw1jVUDLqHDhmsaptNhhqlYByAki4ox3yTSc+CnTu2dwmSSkZMmQI7777bqb0CxcuMHfuXA4cOECFChUYOnRoJvNkV1fXXLuXUvLYY4/x008/5XhmKK/6ubVrjipVqpCQkMDhw4eNCqp79+7cuHEDPz8/xo0bZ9aKLad2hRD4+/vzxRdfcPnyZWbPns3atWtZvXo1bdu2NZZzcsr4f7e3tyclJYUKFSpw5MgRNm3axBdffMHPP//MokWLCjTmUso5tInFFCBJCHEA7T2+Fc25g/mN2SxYusTXyczV2eRfqyClDJVSBkgpK0gpvaSU/aSU/xnyDkspW0gpnaWUzaWUh03qSSnl61LKiobrdVOrvdzqFgW2dvg1N2xNVluTJ0/2fQlXTMLjuHhC25dLziGtqzc4mChMnYcWsj4LXbp0YfXq1fz3nzbzi4qK4tKlS9y9exdXV1c8PDy4ceMGf/zxR766f/TRR9m9ezdnz54FNIOE06dPW1zfw8ODChUqGJfUli1bRseOHXF3d6dq1aqsW7cO0Iwq4uI0Q6Ty5cuzceNG3nzzTaMRwKZNmwgNDeXbb7+lfv36XLhwgXPnzgHw008/Gfvz9/c3Lm0GBwfj5eWFu7s71apV4+bNm5w5c4ZatWrRvn175s6dm0lB5cTNmzdJS0vj2WefZebMmRw6dCjX8mUNKeVAKeVDQENgIpr9wGjgb+C2EOJvS9uy1NXR9rxLKaxFqXtB3+/49IMLwRn3cXcgLgpcSsizg52dtrwnAfcq4OKVoxl6w4YNmTVrFt26dSMtLQ1HR0e++OILHn30UZo1a0ajRo2oVasW7dq1y1f33t7eLFmyhAEDBhAfH4+dnR2zZs2iXr16FrexdOlSo5FErVq1WLx4MaApq1GjRvH222/j6OjIqlWrjHUqV67M+vXreeKJJ1i0aBGtW7c25ul0OhYuXEjPnj3x8vKiffv2HDumWT8GBQUxbNgwfHx8cHFxYenSDCvo1q1bk5qq7eF16NCBN954w7j/ZI6rV68ybNgw0tLSALLNUO8XpJSngFNCiO/RzqL+H9DF8NkihHoZmken08mCnry+evcq95LvFaiuq6MrD7k/VKC6uWHO4u1K9BXiUwq8TGx1PJw8qGzh+aKCWvGdPw8p1gq1lJYGy/vCnYvgWhmaPAtN+qEr78HDufhYOXnyJA0aNMiUZlU3OmmpYFeykXBKo1ugvCjNY8rpOyeEiJNSFmzdNQeEEDo0U/P01bUWwD00x+LbgK2WupkzO4MyaD1LkVLKIfkor1CUKFb9XWZnB21fBpkKNQPA3kYct5awclLct9xGm7/vRLNfmAAckFKm5beh3P6S/LH8QJWahmWhxGM65QOJhAu7NIemFWtq8Y08qoFDybgqtPlnF31VcwD7UPOMtNqdshWz6cUJmQaIEvNAoSjTxKIdGaoMVDJcruRw7jQvzCooKWWNAgqnKI2c3won1mXct3sFmr9YcvLYKlLC9g/g0k6o/xS0+z9wKWV+8JLjtWjJrt7g6lXS0ijKGFJKbyGEDxlLfMPQvPocImOJb7MlbSkXyDZIicwgorL4+61Yu/hlKA2cD9aUE8Cp9bCsj3ZmrDQgJdy7CZGnISVBOx+VbDt7j4qyg5QyTEr5qZTyacATzTjiLvA6YLFZqI0slitKFCk1jwOmVCzymDElSoGW31KTYde8zGmetbXl0NJAWrJ2CJv0rQCpzaS86qn9KoVVEUI4ogW3TT+K1BotBMd/aAd4LULNoBSatVfb8eDzAlR7FMpXB/0DJS2V7SHsoc1YKF/DkGAPAW9oRhI5IK8fh1/HQfztYhMxV+zLgUfVzGmOyluIwroIIf5CM5TYDrwM3AJeA5pIKR+QUpr3cZUFpaCKiNJkvi/t7DXT6I6vQu8F8OLaEt08t9lnZ2enxXV64Wfo+g60HgledbKXO7cdfh4CPz0Ph5fB4eXFL6s5XDy12FHCDso/jFu1RgWaPQUEBBASEpJ3wULWKSwBAQEkJCQwYcIE9u3bZ0x//PHHadq0KY0aNWL06NHGs05Z+eSTT/j++++NbRVU/iVLluDt7Y2vr6/xOnHiBGlpaYwfP57GjRvTpEkTWrZsyYUL2nJ7dHQ0gwcPpnbt2tSuXZvBgwcTHa35d7x48SLOzs74+vrSsGFDRo8eTVpaGpGRkTz++OMFktGKJABvAc0Bbynls1LKz6WU+fZ9qpb4FPclhdKB9g7QIJdYQrfPwY2jGc6/DnwDj/7P7EwrndOnITZW4GYdZ+aYPRfrUVWLY5WLl/OyQHx8PPb29uh0Og4cOMCHH35ozPv5559xd3dHSknfvn1ZtWpVNue1KSkpLFq0yGqeIPr378/nn3+eKe2nn37i2rVrhIWFYWdnR3h4uNEV1EsvvUTjxo2NCnL69OkMHz7ceDi5du3ahIaGkpKSQufOnVm3bh19+vShSpUq7N69O98HrK2FlNJqjsXN/sUIIdYIIeoYPg8WQlg/QJEiR2x2BqGwjIa9wTSG1u2LmpWkrWBnn005SSl57bXXjL/kV65cacz74IMPaNKkCU2bNmXKlCmZ6qWlpTFkyBCmTZuWrZv4+HiGDh2Kj48P/fv3Jz4+wyDjp59+okmTJjRu3JjJkycDmtKYOHEiAJ9++qnRI/m5c+do314L+1KjRg2mT59O8+bNadKkCadOncpxiJ06daJJkyYcO3aMJk2acPToUVq2bMnvv/8OgLu7Fn0nJSWFpKSkHP30bd261RiWI50ffviBrl270rhxY/bv309aWhp169YlMjLS+Dzq1KnDzZs3c5QrKxEREVSpUgU7w4+XqlWrGsOKHDx4kLfeestY9u233yYkJMTorikdBwcH2rZta3Qt1bt3b6PrppJECPGkEOJDIcR3QogPhBA98ttGbj/pnkazZQct5IUy61IoLMGlItTrjpQCaneBASuhVvZzUrbEmjVrCA0N5ciRI2zZsoXXXnuNiIgI/vjjD9atW8c///zDkSNHeP311411UlJSGDhwIPXq1WPWrFnZ2vzyyy9xcXEhLCyMqVOncvDgQQCuXbvG5MmT2bp1K6GhoRw4cIB169bh7+9v9L+3c+dOPD09uXr1Krt27aJDhw7Gdr28vDh06BBjxowxG4F327ZtjBw5kgULFvDZZ58xatQoQkND6dEj4x3ZvXt3KlWqhF6vp2/fvtna2L17Ny1atMiUdu/ePbZs2cKCBQsIDAzEzs6OQYMGGRXCli1baNq0KV5e2c33V65cmWmJLz4+nueee47169fj6+vLpEmTOHxYcxF64sQJfH19MwV0tLe3x9fXl+PHM6+UxcXF8ffff9OkSRMA/Pz8jM+xJBBC6IUQ24Hf0Nwb9UA7rLteCBEshLB4jSA3BXUDSPeRn+6xXGEhNn/YVGEZqckFWw9sOQIG/wYvroFHHrd5K7ldu3YxYMAA7NOSqVzelY4dO3LgwAG2bNnCsGHDcHHRjCkqVsw48zVq1CgaN27M1KlTc2xzx44dxtAaPj4+xnAdBw4cICAgAG9vbxwcHBg4cCA7duzggQceIDY2lpiYGK5cucILL7zAjh072LlzZyYF1adPHwBatGjBxYsXzY7p8OHD+Pr6cvToUXx9fbPlb9q0iYiICBITE9m6NfsMNyIiAm/vzP4UBwwYAGgOZu/evcudO3cIDAw0LsMtWrTIbPyn/qNm3NoAACAASURBVP37ExoaarycnZ2pWrUq//77L++++y52dnZ06dKFv//+GylljrM60/Rz587h6+tLu3bt6NmzJ0888QQAlSpV4tq1a2afSzEwB23/6UXAWUpZBXAGBhvS51jaUG4K6mdgnhAiFU057RNCpJq5rOXVTFHcpCQi47KF9ypRikO5W6xzDn0Pq4bC5b35U1Tlq0HFGgWQrGSQaWmQGAM3/4XbFwyeJjD7ogRo27Yt27ZtM4biWLt2rXF2kG5MYO4la442bdqwePFiHnnkETp06MDOnTvZu3dvpv2U9BAX6eEtsvLtt9/i6+vL+vXrefbZZ5k+fTqzZs1i4MCB2crqdDp69erFr7/+mi3P2dmZrL44s45HCEG1atWoXLkyW7du5Z9//jEqCktxcnLiiSee4MMPP+TNN99k3bp1NGrUiMOHDxsdzoK2fHjkyBGjL730PajDhw8TFBRkLJeQkICzs3O+ZLAyzwLTpJTLpZSpAFLKVCnlcjTjiWctbSg3BfUK0B94B20GtRRN8+V03Z/ueouIYp19hR+Ahf7wTRf4ZQQcXFx8fds6CXc1BXXjKPw6FtaMhDtXSloq6yMl/i0asHLlClJTkomMjGTH9mBatWxJt27dWLRokTGsRVRUxo+Zl156iR49etCvXz9SUlJ45plnjLMDPz8//P39+fnnnwE4duyYMSpv69at2b59Ozdv3iQ1NZWffvrJGJzQ39+fuXPn4u/vT7Nmzdi2bRtOTk7GoICWMHz4cDZv3kznzp0JDQ2lTp06nDx50rgMFxsbS0REBKAtU/7+++/Ur18/WzsNGjQw7uukk743t2vXLjw8PIxyDR8+nEGDBvHcc89lWpbLi0OHDhlnO2lpaYSFhVG9enXq1KlDs2bNMi2dzpo1i+bNm1OnTg6WoyacPn2axo0bWyxDEeAJnDCTd8KQbxG5uTqSwCoAIcRQ4FMp5RHLZVSUCiINm8wJt+HaQc26S6ER+iMkmbgPu3lai61kIaXG1kUInnm6F3v37qHpY88jhOCDN1/mAQ8nHn/8caPCKVeuHD169GDOnIwVmokTJxIdHc2LL77I8uXLjZv9AGPGjGHQoEH4+Pjg6+tLq1atAC244LvvvkunTp2QUtKjRw+efvppQAtpceXKFfz9/bG3t6datWo5Ko+82LFjB+3bt+fKlStUr545mvC9e/fo1asXiYmJpKam0rlz5/9n77zj26ru/v8+WpYsbzve2XsvAgkQEgiUUcqmP2aBllVGn0Ln85S2QPs8hW4oo0CZhUJo2WFD4oQkQAiBTGc5iWM73luyrHl+fxzZlmzJkmzJlhO9X6/7su4959x7rmzf7z3nfL+fLzfffHOfc5x99tlcfbW/3FdmZiann346VqvVLwnheeedx3XXXddveveVK1eyfv367v1HHnmEtrY2brjhBux2lcPv+OOP57bbbgPgySef5Pbbb2fSpElIKVmyZAlPPvlkyHtfs2YN3/zmN0PWiyEHgXOBDwOUneMtD4tEuo1+GEy6jcOth+l0DaytUWdkTHo/eRoGSMDUFO/8hEOb/o6ja9S2+BZYdH3Urx0JKYYUClP7pu0OxEDTbYSVP+/rF2H9X0E61X6E341OB15HtB5cDnDbISk19uk2IkFKaNwHDquKkUotVDp9g4yHG8mpKQAuvPBCfv/733eniofA97R582buuOOOYXVO6OKUU07hjTfeIDMzs0/ZEKXbuAP4E8q57gVUwsJ84DLgeuBOKeUD4Zwr7DgoIUQB8CNgGcq7rxElWfFnKWVNBP0/Jhgxht/VidQYwOPNwpzV+4k69MTNdzfvcqVY/sGvoKMB5oYdAN8XlwO+fh4++QtMOwfOvj96/YwGQiiFjLZKlexQlxSyybHAfffdR3V1tZ+BClTn0UcfjQvX7vr6eu68886AxmmokFL+RQgxCrVMdK33sADswH3hGicIcwQlhJiCSjaVAWwAalAW8USUpMVSKeW+CO5hRDCYEVR5Szl2t31AbZO0SYzNGBu6YoQEG20cbNyPs/mQWhzPnwPJw/fHDZElbBzICEpK2BfJX6vLoaZCC+ZEdJ3uEVRdKTx/iXr4A2h0cMvnlNY742cEFSOOtvuBkX1PQzGC8jlvJkqPLwtoAj6TUkak+xWu1NH9QCswRUp5qpTycinlqcAU7/Govg4KIS4TQpQKIaxCiDIhxFLv8RVCiN1CiA4hxBohxFifNklCiKeEEG1CiBohxJ29zhm07TGNowP2vQ+734ZtK+G9n8HTZ8M/ToenzoSXvwOHNgx3L4cXnSFi4wQ+a1BZE/xVJDwu+OjX0elbggRxipSyWUr5rteb791IjROEP8V3KnCzlPJQrw6UCyHuBh6J9MLBEEKcgTJ4/w/YBBR4j+cAr6LmMN8CfgOsRFlogLuBycBY1OhujRBil5TyvTDaHps4rPDwCdBaTtAwN2s9vHU7zLlMCcoe5fI4MUGXBCt+Da98T+3rkyFPeVn158adIEE0ieXUuRDilEjqSynXhVMvXANlIHg2xHZvebS4B7hXStml7FgFIIS4EdgppezyLLwbaBBCTJNS7kYFgV3ntdLNQognUPOf7wEXhWgbdUZEoO72/0BrRXg93faSyoW04BqY/i0wxE4Fe0R8d5Ey62L4/DGV2uK0X0BaIcaDB2lsbCQ7Ozv+jVRHEyCV2GyCEYeUksbGRozGmL1gltDzlhuOsENYvvjhGqivgduFEO/65pUX6r/qFm/5oBFCaIHjgDeFEPsBI/A6Sqp9JtDt5i6ltAohyoCZQohaoNC33Pv5Au/noG0BPwPlNYQ3gtK4KikpGdC91HTW4PQ4B9RWr9Fz0Bi2J2bYWCwWv/uZsudtwvOV6zpBDay7H+eGh6jNW0ZN/grcuugbKqPWyKikUaEr0veewsHjgaqqvoGMpo5KUiyHaMw5AY+vlt4A0WgkFRU9a5hiwk+RGj1s2QvsRQiB2WymoqIntiruRlRSone1o3V3IhE4kjKRInyN6bi7nygwUu/J7XZjtVopL49Jgk1fLa8M4G/ADuAllCpRHnA56pl7a7gnDfcv7V5gFVAqhFhJj9vgpahptWg53ecBeuASYCngBN4A7gJSgPpe9VuBVG9Z137vMkK09UNK+TjwOCgniYG4MAMcajmEw+0YUFuD1sC47pxD0aOPQ8Hee/wrTFwBsy+B2l3g7ICWw1CzXRkmH/RuK8VH3qG4dTOceDtMi27MRbI+meK04tAVGZiThMcDveIv1YLR609D5SYm1L0Pc69Q30XSwKXFNRoIEVPZh4G6zccEjxsePQnqS3uOFcyD6z8CbXgGPK7uJ0ocjfc0WKSUa7s+CyGeAT6QUvaOyXhOCPEkakbrrXDOG5aThJTyPVTgVTvwC+BhlNGwAOeGm18+DLrkjv8mpayWUjYAf0YFd1mAtF7107x9svjs9y4jRNtjE1srHOmVRmDZz2H0CXDcdSox39n3wxUrYd5VoAkwi2utgw9/CWvjzGV6IJSthspN6nNHA3z6IFQf43HpGi2c9AP/Y3W7oOrL4elPgmFBCDHa61hWKoTYKYT4rxBNzket8Qdipbc8LMJOWCilfE9KeRxq1DEaSJVSHi+lfD/cc4RxjWagksDzlzuBuV07QggzSmF9p7ddtW+59/POUG2j1fcRx5Zn+x4zB1hfSEqFpXfCtW/Bcd8FQ287j/L++3CEe6XteNV/v3gRjD1xePoST8y9HKafpz7nz4Yb18KYY9u36BjEBfxISjkd5Vh2qxBiRj/1NUCwuYPJhLn+1HWiiJBSdkgpq6SUHZG2DZOnUetduV4/+h+iphdfA2YJIS4WQhiBXwHbfJwcngPuEkJkCiGmATcAz3jLQrWNOnETbBqMve91f5QAIoSfi3kULLkNrnsbJp7Rt3z3W7Dt5ah0LdbfXcDTn/sXOP5mEHoQOjj5R4NWUQjrNlqr4PVblKxSPCIEnPtXOPUuuH415PX3XEpwNOKdzdri/dwOlAL9BSq+DfxOCHGp168AIYRWCPFt4Leo53lYxGNG3d8AOcBeVOrgl4H/lVJ2CiEuBh4Cngc+R0lndPFr4FGgHDVVeL93ahIpZX2Itscedb1ssznM4FyDGc65Hz42w67X/cvW3gfWBjjh5pDZY+MOnUGlcJ94GtTugFHB0tFGCXs7bHgANj4ELpuaYpxxQeh2w4E5G5b9ZLh7kSB26IQQvrnsH/euxfdBCDEOmI96hgbjB6hZtpWASwjRDGSi7M16b3l4HQu34lAhpXSiPANvCVD2ERBQOVJKaQe+690ClQdtG28MyegrNRdsjT37+XOD1w3Eil8pl+PNvcQrN/8DmsrUG/cwK1IMiJxJaos1HU2w4UGlywfQXg2fPQwsiv21EyTwx+VdvukXb6LBV4AfSinbgtXz+g4s9ca0LkbFslYDn3qfw2Ezwl5zE0QFj6dv2ojei+HhsORWOOfPfR0oDqyB5y9S6hQJApM5Fk64yf/Ynve68zCNKDzu4e5BghgjhNCjjNMLUspXQ9UHkFJ+KKX8jZTyFu/PiIwTJAxUzIjrYNOmA+Cw9OwnpUFqwcDONXE5XPQ4mLL8j9tblYffS1dCR2QKJ7H+7qQEmg/F9Brd1+mPpT8CUyaYc9U6z3ffV0riIwV7O7x5O7z+/eHuSYIY4o13fRIolVL+OcK2GiHEaiFEcLXdfgj53yCEMAghtgghvjGQCySIQ6r946rlqGmDcwgomAPffg5yAsyg1pfCU99QruieOBkd7P0Qnr8Ytr40vP0wZcAVL8MPtijXfm3czbgHp/xTFSO15TnlxbnzteHuUYLYcRIqfftpQoivvds5YbYVwHICxJyGQ0gDJaV0AONRroYJhoCYj756x/eMisLSXFohfPtZFbir7ZWqQbrVQ+zlq6F84/Bm8qv6ClZeAUhY93slPzScjD5eufKPJDweJSrc4qNIsOoOaE9k3TkakVKul1IKKeUcKeU87/bOUFw73PmED4HECOpooWab/35ulHxHtHpYeB2c8ycQAUId6kvhzdvgpSuUS3pna986scRhg2fPVQazi02Pwb5AiT8TBEWjgQseBa3P2mNyDtgiFqtOkKBfwjVQfwMuF0L8UQhxshBiohBigu8Wy04miCJuNxz+zP/YqOmB6w6UcSfC1a8HP2/DHuWS/uSZ8PG9arpoKFh7HzgsSNlrOvO9/4Y1v1Pq7gnCI28mnHaX+nzCzXDzJ5Ab5b+jBEcDHpQA+JGBNA530rtLZ+lOVJbEQIQdHXwsELeBuntWgV8qeoFMHx3966QXwWUvwKGNapRSu71vHY9DxVLteh2MmbDgapj/neh/d1LCZ4/Chr8GqeCBHf+Gw5/Cqf8TVaUEKQewvCcllL4FbVWwOI4dEJbcBmNOhNEJ1/gEPQghiqWUlQBS/TPf06t8ma92X3+Ea6Cui6yLCQZDTI3btn/77wuNN6g2Rtccd6LaanbAlmegbE3ga3U2w8YH4ctnYdlPYdFNypBq9aAzKl24gdDRBK/dBPvCkItsq4Q3boHJZ8LJd0JKeIrq0cRoq4YXLoX9H6optElnDE1c1kDQaBPGKUEg3hdCnCSlbOld4M0btYownSbCMlBSygDCbQlGJNVf+e8bM4bmuvmz4Jw/Qls17F4FO18HS3XfevZW+OAu+OBXPceEBjLGqTLzKOXUUbQQdEmMPrwLttYoN3dzjpp60nmdNOr3wr++rVLZ+yDRwUVPqkX9zx/tScXexb73oXwDLL4VZl86dKoYHjdzt/4aOmvVvtsB7/wIrnpt5ClzJDiWsQBvCyFWSCm7p2uEECejZJDeDPdEEfm1CiE0wAwgG9gspUxM2o80egdV5gwoPGHgpBXA8TcoZ4qvX4AvnlCpPfpDeqD5gPrc0Qj1u7tlliYCHHimp64+WRldW5N6wPcOfE1Kg3MeBdN8tT9hOWx8ALb/G7+RncMC6+6Hg+vgjN+AuVecVyzQaDk4/ipmlP6p51jdbmg/AunhpR+JGxwd6veU4FjkmyhJo38LIc6XUnqEECcC76AM1FXhnijs1zIhxK1ADbANWA1M9R5/XQgxABmCo5u4DdTV9XIBX3Lb8PRVq4OF18DN62H5L1SwKmqxZlC9cXaoB7qrs69xKj5eXW/6uT3HDCZY/nO49NnAcVwVn8JLl0PF5r5lYRDpbG1d7lIYf0pPf29aO/KM05Gv4PFl8M8LMdgbQ9dPcFThlTo6C1gAPCmEWAK8C7wPXCkjWMMIy0AJIW4AHkBlt/02XU8SxSfAxeFeMMHwoXXZ/BUUhAYmrRi2/nQz+2L47nvK82/OFWoUBGrtSdc38+2AmH8VXPu2khgKRP4sFWy89CegN/uXddTD6zermKlYy/oIAec/rKYXr10FqfmxvV602fAA/ON0aNgLnS1M2/1A/ARoJxgypJSHgLNRWc3XAh8Bl0kpI/oHCneK707gT1LKn3XJp/uwG5WSPUGUiNWIxmw97H8gawIYop+yfcBkjIZlP4blP1N5p5Kz1dpLey3sekOtFVVvU5l+c6ZAxmgqq6opTpVq5FS7q0/2XxAqn9F5D3W70wV9f9PqYN7lStH8/f/ptV7nUd6I9XvgzP8DvTEW34AiYwyc9X+xO38sERrw9MT0ZzVvha+eg4XXDl+fEsQcIURAkW7UetPZwAfANaL7f1A+Fc55wzVQ41HDs0BYUTnoE8Q5Zmu5/4HcOM3to9H6e9Cl5qlUGAHYX1JCcVf6bY8H9ryjtAZNmZA3y+s0ESLXVW9S8+DCx5RB2vwUfpOOB0vgP9fBWfdD5pjIznsssPhW2Ps+HPoEgCMFZ1A465Jh7lSCIeAfIcof9fksgagaqAZgXJCyqUBVmOc5ZojHOKgUywH/A3mzhqcjsUKj8V9fGgxanVJrL1oAH/xSOV100bAHXv4OnPN7JVXUD1H/M/C4wdYSOPtxPKDRwIV/h2fPg2/8hr01ZgqTUoa7Vwliz/hYnDRcJ4m3gF/1UoyQQogcVODu64GbJYgn0tr2+B8omDM8HQlBXBn3MUtUwHHmRP/jjjZ44zbYunLotAWtDSqNyfMXgcs+NNccCOnFcNsXMO2bw92TBEOElLI8ki3c84ZroO4C7MAO1GKXBB5Epf51A/dGeD8J+iEmD+iWir4jqMKF8WUM4pWUPLjkSZh4uv9x6VKu6KvuVAHBseTw5/D3pXCgRKnRf/KnkE2GlYEGVicY8QghFnvTuwcqu1QIcUK45wrLQEkpG4HjgN8BeqAMNT34ELBESjnEqp8JIubTh+mjuJOaOxw9GXYGZJONaXD2/UrepzeH1sKLl/XVOIwmnz2i3Oe7+OTPULszdteLJbE25gmGm98BM4OUTfeWh0XYcVBSynZvVsSTpZRTpJRLpJT39Jf6dzAIISYLITqFEM/7HLtCCFEuhLB646+yfMqyhBCvecvKhRBX9Dpf0LbHBGW9kln2dqVOEBoh4Ljvwtl/COCK3gBv3K6m/Hxc0aM2QP3WX9VIrosxi72xYyOMA2vhgXnw9YvD3ZMEsWMuEOxtbRMQ9tpCRPopQog0IcSJ3mHaEiFEWiTtI+Rh4Aufa88EHkMlzsoDOoBHetV3eMuuBB71tgmnbdSJu0Dd1l5+LNlxqu9GHH53vZm0Ai77F+TN7lXgVlN+L10FNQHEcQeDKVON4ABO+qGKGRsGrcBBsfUllSjS3qrSrhwISy80wcjDSHDbogXCfjuOREniV0AFKjB3JbABqBBC3BXuOSK41mVAC/Cxz+ErgbeklOuklBbgl8BFQohUIYQZFSz8SymlRUq5HuV/f3WottHue1zidvZNJT7vSmAEGIN4JWM0XPwPWHQj9J48bdwD/74W1vwfWKOopDDzQrhxLZxxz8jKvgvQdgTe+i/wONW+xwUrr1bHExxtlALnBSk7D9gTpKwPYf2VCyHuQT3U/wG8BNSiRiKXA/cIIXRSyrvDvWiIa6WhnC5WAN/zKZoJbOzakVKWCSEcwBRUzhG3lHKvT/2twLIw2n7Z6/o3AjcC6HQ6SkpKBnQfFR0VA2rXRVVy9Dz301pLWeCwdO879BlstE2BkhI80kOVLb6iBLRCy2HT4dAVAYvFEvHvyG7XUFeXFLpiOIj5ZE66gYllT6OVTp8CCTv+g+OBx9k9/SYack7o+5IQgLDuZ2+I8hiicdsx2WowdtZgth7GbD2MsbMOvbMNjceBS2fGqU+nI7kYj0ZHnkxmV+06Oo25pI29nEllT3efq6zoAiq+3ANibz9XjD8G8jd3jPF34DEhRBvwBFAJFKGeq98Dbgn3ROG+ht2AUpLwVYzYCawWQrR6L3x3uBcNwW+AJ6WUFcI/kU4K0NsZoxUl2+7upyxUWz+klI8DjwMYjUa5vCsINEL2Ng7un25K9pRBtfejxH862DDtDJafeioAHulhf9P+6F0rCug0OiZkhpcDs6SkhEh/RzYbVAzu/aEXJ0DbebDxIdj3nl+J3tnOrJ33KXWIOZcpcdr0IkgtDBhAPJD7iRl2i0pTcmANNOxXqvDtARTo/do0AJDZ4pO12e/9x/s/PWYxEydPZWJKozeoeuaIkXWKq99RHCKlfEIIMRUVgnSnbxHwF+8zNizCNVDpBFeSeA+ISlY1IcQ84HRgfoBiC9B7zSsNaEeNoIKVhWp79NN7rn/C8uHoxdFNWqGSJ5pxHpTcB629RoAth2Hd79UGgFBOD6YMlS49axwYUhhf3QCaL5QThiFZ/dSbIClF1UvJ9UpA9XLjttSpKcUz/xcMA3CA8XigtQKaypRLe91O2PcRuGwD+Tb6wTulfPhTtfkyZoly/ph5EeTPHkCmxwTxgpTyx0KIR4EzgCyU2MNHUsoD/bf0J1wD9TmwCBUD1ZtF3vJosBylWHHYO3pKAbRCiBkoQzi3q6I3aDgJ2IsyUDohxGQp5T5vlbmoUR7en8HaRp24ii2yW6Byk/+x8cu6P8ZVX73Euk8xPf2YxXDFy/DlM7CpvxdFqXQDu7QDy9cDMBYg1Oym0ChjVTBH5cayt8PO18DeBg374IJHlMHU6oOfw25RI6LmcqjarLwP24d5PajLaK3/C6Tkw4Rl6m913Mkq+DcRWzWikFKWoUKSBkxQA+XN/dTFD4DXhBAu4N/0rEF9G/gucP5gOuHD46g1ri5+jDJY3wdygU+FEEuBLah1qlellO3e/r4K3CuEuB6Y5+3Tid7zvNBf26Oa8o1+4p1kT1IL/Alih86gtAPHnoisehDK3lDpP6KF9IC1DvZ/pDZfytfDA3OUEnzxImUwTZnKGaHlcM9mG2AsktAoY5E5HrInKj3HnClqZKcz9hjJlnIQWipLN1Oc4lFCv7YWdTwcLDWwbaXauq6bkq/yl828UCWSTEgoxS1CiGSUbViGGkE1AiXAM1LKEAngeuhvBOXCPzWPAO7zbvQ6vi3EucLC2/HuzgshLECnlLIeqBdC3IwyNtmo0ZxvKvpbUAKEdagv4/tSyp3e8+4M0TbukFIiojHF8U4voXmf0VOCGJM/CxY/Dp4/QOlbypg0l0NbFVhqY3ttV6cSbPWKtg6Y9DFq2nL8MmWQ0keHFt8t6J6sYL/LR8wXoLNNueBXb1WBxg6LmlYM5ZYvPWqE134EDq6FVXcoNf5Jp6vYtNwAubwSDAtCiHyUMZoClKPyCE4ALgFuF0Isl1KG9Q/Qn1G5l0HmjhssvT0DpZT/Av4VpG4TKvdIsHMFbXtU03LIf9/Hmy/BEGFMV/mo5vskEnXZ1ajGYVU5uqx14LRxcM9Oxhfnqoy0Tqv3p02NTKz1aq2psyU2/dToVb6swgVKBLdwgRLLjeZakDENxp2kNl/W/QFW/y/hP3KkMmybypTqfNFxamS14Gr1fScYTn4PZAJLpZQbug56s+q+AtwPXBvOiYIaqGi5jScYRrau9NuVgDjph8PTlwiIdWxWXCy76ZIgyysAnd+jKl9uL2F8KA8xl0OlFDmyRY3ErA09XnZavfK06wgRf6XRqam6rInK6WL8UvWATxqG0MDKL2H1b/2PpeTBvKvg8EaVOt7W3P85qjarreQ+OO5aOOHmkZeJ+OjhbOBnvsYJQEq50Rs323sWLigjLNpvZBA3wa8b/+a3K4UGkeefAypu+nqUEhNjqDOoKa1g01pSqlFZ+UY1lebqVOuO6WOUu3vGGOXSHS9OB8UL4eQ7lHMEgNDCil/5jzidndC4D7Y8p1zffTND++JoV3/3nz6iDO6p/6OmJhMMJSlAMI+bSm95WIRtoIQQ01FziKNRUha+SCnlNeGeK0FoJBLRV941MkyZfrs2U1H4GiMJRi5CqNFZ1niYrxRDkDK+3bZX/Bos9bDjFbj0GZh6ln+53qhcz8/5A/AHb4zW+/DFk1C+oe/5pBt2/EdlYj7hJjj5zvjNoXX0sQel4vNegLKrUFnYwyJcJYnvoBwQJMoJwdGrSuI1PN6wW6DCP0B317Q7WDRM3UkwjNTuUtp35z0EeXGaRVkI+NYDsOQWFbQbiqQUmHWx2tqOwFcvwOd/V6K9vnic8OlDsPlpWPRdWHK7ypicIJb8EXhOCJGHWvevBvKBy1Bxrlf309aPcLX4fgm8AYySUhZJKcf32sIL+08wdBwoAbfPe0TmeKypI+PXNKLjoOKNQ+vhqbOg6kt4+myo+CJ0m+FCqwvPOPUmrRCW/QTu2KGMXE4AFRanVU39PTAH3v1ZQgMwhkgpnwduBmah5PHeBp5EqZjf7HVYC4twDVQ+8IiUMkbuQ0cXcRH8uq+X8MeUM+N7iucoZVj/FKq+hH9eqNTDQXn/vXAxNB0cxk4Ngq0rob0meLneBAuvhVs+h4v+AWlFfeu4OtVI64G5ylW9JTzNxwSR4ZUzKkTpoC71/iySUj4RyXnCNVAbUImmEgwRgzJyLgeUrvI/NuXMkNcp+7qOVQ9tZe2/dtNp7T2Lm2DEUTAfFvRaGp5xfuAHd7yz5Tl4xW1lRAAAIABJREFU7UY1CmwJIaSo0cCcS+G2zfCN3wbOm+V2wOan4MH58MatI9doxzFSSo+UslRKucH70yOEOF0IsS10a0W4Buo24EYhxOVCiGwhhKb3NrBbSBATyj72VwowZsDYk4LXB1rqOvjs9QO01HVQsbuZL98LM+I/Qfyi0SingqU/Vp5xZ/8evvVg6EDbeGPby/DmD9TnpgPw9DnhGRRDMpx4O/xwm7r31MK+dTwu+Op5eGgRvPvz6KZHSRCIdIJn2+1DuIalEvgKeB7lJOHstSVet+OJ13pp9868UMXd9MO+TTVIT89o6uDWBhx2Vz8tRi7xMAM7ZAgBp90FN61T3mwjbZpXSuWJ5+uHZakJ7mYeCL1J3ft/fQ3n/kW52/fG44TPH1VTf2v+Dzp7Jz9IMByE62b+BPD/gNdRLoIJgxSv1O2Bzl5BjaMX99vE7fJwaEdfbbaqPc2MnzP0WVsTsVlRRgi/YOA+eNxKpSIe010IAZc8Da9erwyVRgeXPgsTT438XLokJYs0/2ql8ffJn9SIzBdHO6y9Hz5/DI6/AeZenoijGkbCNVDnAz+RUj4Qy84cLQzrA3ZVAKWIWRcFrS6RVO1txt7h7FN2eGfTsBioo4kRMVr75M/KFfu8vyndvXhDZ4CLnwLDD5T23rRzBnc+rV4FAc+5DLb/G9b8r0o14ktni5Jf+uRPypX9lJ/AqKmDu26CiAnXQFmBXbHsSAJ/Bmzkqjb776cVhVxzKPuqPuDxmgOtuN0etNrEEuNRS8UmKPmdCmx9+WrlBXfm/w0sp1Qs0epUGpFon3Pe5WoK/IsnlKHurfIuPcqIbf+PqnfKTwKf6xjHm8IoHCIapof75HkauCKSEycYBmp3+cc+gYqg7wdri50jewPrnDntbhoqjj5x2RExqhkK7O3wyveUcepi99sqyHuk4XHD2z+GxgGkH9Ibe5wpTr0LkgKJzUrY+So8uoQZO++Hmh2D7vJIQQjxlBCiTgjR303vB/aFsf0t2AkCEe4Iqhy4XAjxIUq+os8TTUr5VCQXThADPn3Yfz9rAhx/fb9N9nxe3e8Du6aslbxxvZMRJwiXuDaGejMsugE+vlc5CQCc/8jIVFpY+3s1Ctr6ogrWnX1J5OdISlUBvyfcCDtfV/FSdX0njnLrN8LfT4Jp58Kyn/qlFzlKeQZ4CHiunzoxSV8UroF61PtzLLAiQLlESSElYJgCddtrYfvL/sdOv6ffJlJKSj+t9juWXWSmscravX9kfwtzVwx9gsOo5cNKEByNBk76gVIyf+V6tb4z5RvD3avIOVCiHBtApZN55XtKzf2EmwZ2PmM6LLxGOVPsfksZv9oAg4fdq9Q29RxlqArnD/gW4hkp5TohxLgQdZ6NxbXDNVDjY3HxBFHki3/0kjYaB9O+2W+T6rJWWutt4BW11ugEiy+YyNsP98TRNR6xYLc6STL3kz48wcimcL5yQxf9qJvbWtQII14U0H0pXYWfG7o5F2YETQ0XPhqNCmye9i3Y844ygjUBYkz3vKO2yWfCsp8pdfaRhU4I4bt4/bhXCWLYCctASSkTUZtDTESjsM5WNb3hy+JbQj5MDm71F9YsnppFZp6ZzHwzzTXeUZSE6gOtjJudE35/4py4nnYbLvpzivC4YeVV6ou76LH4y7N0zh9U6pH3/ke9pF38RHSnKTUamH6ueuHb+x5tq+4irX1/33r73lfbxBWw/Ocq6ePIwCWlPG64OxGIhHvW0cDrt/ondEtKh3lXhmxWvt3fQI2ZkQVAwUT/ReLqskTQ4kA5KozhxgdV6vjy9fDoSd7A2ThCCFh0PVz/kTJWE5bH7jpTz2bLgj/CFf+GoiAjpbKP4ckz4LnzYe8H4D46A96HgrAMlBDioBDiQH9bNDojhEgSQjwphCgXQrQLIb4SQpztU75CCLFbCNEhhFgjhBjbq+1TQog2IUSNEOLOXucO2jbaDGkc1IYH1Ty5LyfertIR9ENbg43mmo7ufaERFE5Shqlgkr+BqilriQ8B3ARDT1WvbLedLWo6OR7/HgrmqODaYNRsh3V/BHffmL+IEEKt1V3/MVz1ChQHGSkdKIF/XQqPLlFG3ZXQN4iUcNeg1tI351M2cCJgAVZHsT8VwDLgMHAO8LIQYrb3Oq8C1wNvAb8BVgJdMgl3A5NRjhz5wBohxC4p5XtCiJwQbUcmrUfgw1/5HxMaWPz9wPV96B37lDsmFYNJrTONGpOKTqfB5fIAYG110FpvIyM3OTr9DoOoJGxMMHjSimDcUjiwRu2bsuDCx0aeZJLbBW/cBtVfKw+98x6EogWDO6cQyrFk4ooeR43Dn/at17AXXv6O+u4WXgtLf9T9AindbqTLhXQ6uzc8HnC7kQF+SocTj9WKp6MDpEeVeaSK1/J4lFyZ72ePp0+9Pp9D3qZ4EVgO5AghKoFfSymfHNyXFx7hrkFdG+i4ECID5Xb+UTQ6I6W0ogxNF6uEEAeBhSiDuFNK+W/vte8GGoQQ06SUu4HvANdJKZuBZiHEE8C13v5dFKJt3BHWKOypM+jz3jDxtJCjJyklu3t57xVP68m+q9NryR2fxpF9PdlVaspah9RAxZJ4fPmPW1Lz4apX4bOH4aN74PyHVP6lkcbnjyrjBFC7HZ78BvKOnUhdGu72djzezd1uwWNpx22xIO0OpMO7dRkQl4vUgwepLinxMyrS6UQ6xkJ7MrLhILKjDekB6RHezfv5Xy8gPS8ipdZrQOL/j1FKeflwXVsMdupGCHEx8L9SymnR6ZLfufNQMVjzgO8DBinl933KdwC/Ro3gmoB8KWWtt+wSlKWfLYR4IFhbKeUrva55I3AjgE6nW/jhhx9G3G+nx0lNZz95a8Ig35iPXhPcc67gyAdM3esf9+TSGFl/ysqgbSwWCykpKdgaJQc+VL93Jw6adLWMP02gM/a8FbcclNTv6vnbMOdC4aKhW7IsMhWhCUMkv+ueIsFi0dHcPDReiTk5Dkwmd+iKXgZyP0NFUmc9dmNw6ausxi00Z85G+vzdDsn9eDwIhwNhtyM67WisFjRt7Wja2xGWdrRtrYw6tAZsDtwODR6nBqcrCY8DhDv8383Ryow9uzuklHEmHaIId4qvPzqBqLv1CCH0wAvAs1LK3UKIFKC3Jk8rkAqk+Oz3LsNbHqytH173yscBjEajXL58ecR9d7gdHGo5FHE7X8ZljMOgDSBRJKWKy+hlnECgu3kdy3OD64WVlJSwfPlyPnxqJ1Dbfbx4SiYnLfN/v2gZ38GqXVu79+1NGhYuOA6dfmjcjCdlTQrLQHXdUyS0tEBd3QA7FiGFhRDJ83kg9xMXVG+Dx+6FnMkqtfyYE4CB3Y90uXA3N+NqasLd2IirsQlXYwPuxiZcTY3en96y5mZkR0fIc1rRAEafI+7EBPIIYMAGSgihQ6X0vRvYGa0Oec+tAf6JUk2/zXvYAvSWNEgD2r1lXfudvcpCtY06MXMocDvVgvWGv/of1xnhmrehH+PURVujjX2b/Z/Okxb2dclNH2XCnG7A2qoWdl0uD/Xl7RRMyhh4/yNASkniCTKC+PgeQKr1lqfOVJp1p/2iTzWPzYajogJnRQWOigo8bW14rFZc9Q04q6tx1tTgqquDY2hkIzTSuwFagTAkgzEFoTOAVovQaHp+6nVozSmIZBNCqwONQAiNcoXv77NGo9an/T5r1P/YL/r+nuKFsAyUEMJDXyeJLtqA/iNCI0Ao+YAngTzgHClll8vNTuAan3pmYCJqbalZCFENzAW65uTm0mM4g7aNVr+jTR8j99G9sPcdqCv1P67Rw0VPwOjwwhi2flzhl/cpI99I0ZS+RkcIQcHEDPZv6TFmR/a3DJmBiiWJNagoc2At7FfL0G6HwGnV4jzowPXSSsxfbOLIBx/iPHwYx+HDyvjEEUKvR5OaijY1FU1qKhrZjtZejWbyyWgy8xB6A8Kg9/mpZ/+hcqbMmI7Q67s3fD4Lnc9ng96vntAIRO1XiM8eRFRsABHM30SoJKOzL1EitaYY/t+NdAMF3EtfA9WJWh96V0oZzUCZR1Hp5U+XUtp8jr8G/MG75vU28Ctgm4+Tw3PAXd6I6DzgBnr0oUK1jU/cbpXmeudr/oKeXSSlwxUvwdgTwzqdyy7Zs8HfOWLuitEITeC31YJJ/gaqen8iHipSjhZjKB0OXPX1OOvqcFVX4zxyBOeRI7jq63E11OGunISrxYLH6Z2WfV/JbqXgP+8eK4TJhMZsRmMyoc3IQJeVhTYrC112FtqsbLRZmegyM9GkpaFNS1NGKS0NkZTUI6nVchgeOVHlhNJUwZJb4LT/Vuk5fNheUkLmYKZhc8+CWWdC6Zuw/i9w5KsAlaSKOytfD+/9t1K0WPAd9b8+0jwoB0G4Xnx3x7gfAHhjk24C7ECNjxbbTVLKF7wG5iFUZt/Pgct8mv8aZdzKARtwv5TyPW//60O0jS1SQkcDtNeAtQEOf65GQ0474HXzzJ8HGcXKrdeYBi1VcHAttFUFPmdyNlz+UkTR6s37wWXvMUbJ6QYmL8zniC3wNfInpCFEz0O2pa4Da6sDc/oISxmeICjS7cbV2Iirrh5XXZ13q1WGqK6u+7i7qW9Cy75Ez4lGm56ONidHGZrsbHTZ2Wizs9BleX9mZ3cbIY3ZjNAOcm3U44E3blXGCZR47oG1IfUsB4wQyujMOB+qt8KXz6rU9o4Aqw4uG2x7SW1ZE2HB1TD3ipEp6hsh0XCSiBpeSaWgrwdSyo+AgN6CUko78F3vFlHbmHCgRKWlrtkOtTvBGmJqo/pLtXWjIehXMeNC+OYfwRy+/JDL4aZxr//r/NzTRqPVa5Q5D0BSsp7sohQaKnvSL9SUtTBxQW7Y100wvEiHw2uA6nDW1uKqrcNZWYH94EEcBw/hPHJkaNd7hESf7MaQqUc/aQa66SejSU1Hm5GBvrAAfX4+uvx8NEZj6HNFE1eneunz5YSbhma0UjAXzv0znHEv7PgPbHlOBUgHoqkMProbPv4NTD1bCdpOOl3ltjoKCXpXQohfBSsLhJTy3sF35yjA7Ub+YRLYA+dYGhQanXpzOj+ilCoA7P6sBre9Z19v1DLzlCLc9B/dXjgpw89AVe8fGgMVSzWO4Zh2c1usuGprcLe24m5rU3E3nZ0qMNPpxONwgMuFef9+6rZ8pYI3XU5wuZAut/++06X23S7//a6towO3xYKnvR3pGFr1AmEwoC8oQFdYgD4vn6rOTiYuOg7D6NEYrNvQb/xFjyatyQa3/BP0piHtY0AMyXDpM0pk9t2fgscFMwNnohYetxrxzLpICehGi6QUFci78FqV12r7v+GrF6D1cN+60t2jpp5aoKTNZl4IuTOU88NRQn9m9+4w2vv+qycMFMDKK8A++Fn37i9Wa1Ciryf/SOWqGQAej+TrD/3/yGctLSLJpKPD2f8DrGByBttKKrv3qw+04fFINJqjYx7c45F0Wp0Yk3VoBpk52NnQgL2sTK3R1NTgqq7B1rYXY8MBPK3h/U2kAI2D6kWMEAJtTjb6Ubno8vLQFxWhLyxEn5+HNisLbWbXek+W8hLzsrekhKyu9Zp//rVbOR9QayrxYJx8mXmBGpHU71aJDAOQ3bgJ1t0H7/8C5l2hUnPkzYxuP7InKsHZU34KB0vUqKp0VU/eLl/aq+GTP6otrVhNAS66AczZfeuOMPozUKEiGecBvwXORGVKTOB2w973A5dpk1QKDPMoSM5R4q5aHSSlKZfPMSeAtRFaK6GtErQmpQqx8Fr1djcIDnxVr9JqeNFoBXNOCy/HU3ahGYNJh8OmBC/tHU6ajljJKY7PYNJIsLbaWffSHhqrrGQWmFl+xbSw19eklDhra+ncth3b9u10bt+Oq6ZvcHayth6D1hrgDPGDNj0dXa4yPLrcXHS5o9Dl5qLPze05np2N0A1yGmnZTyElF3a8okYoi74XvO6OV2HiqWDKDF4nViSlQHFwr9iiqrfVB0c7bHpMGdkzYrRWpdGo58DE09TzYdtLyljVB/HvaquEkt+pLX+OegmYfWlsvQBjSNC/OCkDuY2BEGIyarR0KVCFUl14Oia9G2m8+1P6ODsWLoST71ABjNoI1AvSx6j4pkHitLvZ+Ip/aoApx+eRkpkUVnuNVkP+hDQO7+xZJK8uax7xBqq6rIXP3ziApUXNezZXW/lk5R7OuG6mWpfrhZQSZ2Ultm3b6dyxA9v27bgbGvrUizt0OmWA8nK7Rz+6/DySxo3DMH48hrFj0SQPkYTVmMVqO+M3Sh09Y0zgevV74T/XgT4Z5l4Gx9+k0mnEA3W7yWzZ7n9s7hApAZmzYcmtKpVO5WbY8qwy5M4gL0A12+CdH8MHv1TTkVPPgQnLojstGWPCfiUSQoxGecp9B5Xy/cfAI1LKhERvF963mm4TpTGo3DTDyKa3DtDe1Nm9LwTMP6NHyD2coOLCSRn+Bmp/K7OXxTbLbqyCnTstTkpeOMDXG/p6pTVUWvjqzV1MH9uJq7YWV00Nzppa5Vxw5AgeiyXAGcNEr1cOAFlZytU5NVUFW/aKnymvrGT85ElqtKLTIbQ6hF7Xs6/TI3Q+x7Q6/32dDo0pGU2KWV3DaIy/zMQpo9QDMxibHlM/nR2w+Smo2QHXRy45FhOSszk09tuMa1gD1nqV7DGY8Ww7AmWr1bpWCH3MiBACRi9S21m/U0ZqzztwcJ36znrjssHXL6hNo1eu6pPPgMnfgJwp0etXDAhpoIQQucAvUCOlTtTo6S9eYdcEXZR/qmIWfFlyW+C6Q0T94Xa2flzhd2zOaaPJKoxMdit/ov/0QH2FhU6rA6M5uu7mUkpwOpFuN25NG0IKpNMFLqe/E4DPMf3u3Vi0Oq+QZ5dTga/jQM+x+hYd6w8U0Gg1EmwGe/fXbWhefJZ0S4CF6VDodCRNmoR+7BgMBYXo8vMompFN9pQ8dDk5fmszwdhVUkLOSJQ6iha2Fvj6Rf9jA03dHgtSRnFo/JWMu/oh5aCg72f0ufkpWPcHWHUnjF+qVMzDjFkMm6RUtQa28BrobFPGfdvL0LCPgNoKHqcKXzm4Fj64K/goNk7oz4svHfgZcDvK3/kBVGxRDNzTRjhuF7zby4GheBEsuAq3y0Pl7mZMaQZGjU4J+212sCMIj9vDmud3+3ms6ZPh+G+Nj/hcKRlJZOQm01Kn3s6kR3JoWyPTlhSEfQ5XSwv2vXux79ungjsbG3E3NeFqbFKpA1wuFYvixWHRo5ehv6ssVH6WULSkT2TrnFtxh3LHFYKyiRcwZ/vf0XlCTA4Y9BinTcM0azbGObMxTpvWxz3alAf69CDtE/RFaOCUH8Gmf0D7EeWhNuP8wHUtdSq/08wLYfQJQ+u9pkuCWRcHL3c5lKcfgNuulDYW3xLbPhnTlMTUKT+B1io1YvryWbUuFYyWAbyIDSH9/bceBNKBD1DOENVAphAi4KqllDIqSQtHJJufUvFOviy6EZvFwUfPlNLqfbCPGp3KSZdMIiUz9jEe29ZUUn/YP+iv4DiBwTiwhe5xc7L5+qOe6YOyLfVMXZzfr8F1NTXRXlKCdd0n2HcPn2hHc8Zkts+8EbfWu+7mNXzmjhrGH3wbt8ZA6fSru+vbkzI4NPYsJh180/9ERiOmGTMwzZ6FcfZskqZMQWNIBC1HFWOaGmmc+AM1QnE5gq/dlr6pRgybHlOGbMltcOLwzlp0U/qmf+yjIQXGnRy4btMB2PcRTFqhvPeiQXqRcko5+U41Wtr3Iez7QMVRjSD6e1p1zeucCXwjjHMNjcx1vGGp8884CsiZF+HOm0/JUzu6jRNAfUU7Hz9XylnXzyLJHLt0D22NNj5/66DfsUkLc0kqHPii/oR5uUrHzzsia661Urm7mdHTs/zqSSnp3LGDtlWrsKzfMOAgUCkIrv4YJi5tEvsmXkx14Ul9ygqrNzC6cjUa7w0VVH9GdUFP/sr6UfMonp3P6Ekp6PLz0Od5p+kGq1iQIDy0ejUy6o+dr/d8bq/GL9BvuBl9PJz0XyqTbvMh5YWnC+KYVPpWT+LRzPFw0g/guIB6A5Gj1SnDN2kFnH2fiq/a/5EyVofWqwDlOKY/A3VdP2UJABw2eHA+OHwWzw2pcNov2bn6EI1VfZfp2hs7WfvSHk77zvSYpK5wuz2sfrbUT9LIYNJx8rcn88VXfQ1UuAGxyWkGiqZkUrmnZ4Z3y/vlFE7OQKvTIN1u2levpvXV13AcPNjPmUKg1YJOi4ZktFoD6H0cA7yb77FWi4XMnBy/Y+j01MgCdtonY3H3XW+blttA/uQxaMy3oU1NQZOaSkFmNh+/1eTnjr+zuYjRC2eTnD3wWJ2jRYsv7mivVQ9YX2ZcELiu26Xcs6ecFZH6yqDIGKOUIU6/Rym8e1zB6+7/uOdz80HV31iRPVFtJ9wEjg71Hd5zZuyuN0j6czN/dig7MiJ5eJG/cQJY/nOO1BjZsS6Ihh5QV97OZ6+VcdIlkxFBAl4HoqQgPZK1L+yham+L3/ETL5qIOT08t/L+mHv6aKr2Nnc/dNubOtm7qYZxqY00/P3vOMr6n+U1TJxA0rRpGMaORdels5aVpbzNurzUvFOGQfNh9aKkpIS5Pk4FLbUdrHl+t182YF8WnDmWSSdPoLm57/d+0iVZvPf4DjxudYMOu5vVz5Xyje/NIjnt2JnKs9tcNFS0o9Vp0Bk0pGabcDnc2DtcSCmRHjCl6DFnDP5vasCYMuD/Pa+ElPe82/PgDUT5eqWzh1Bed3MvhxNuHJp+CgGj+kmDY7f0TRM/aUXgui0V8Nz5KkaraKEapRXOH3jfDMkwJZzJseHj6BRwGgpevxVaey3Pm0dRX3Al7//tCzw+KS0MyTrMaUk01/SMqA7taMScZWT+6dHxovF4JGueK2X3Z/7BokVTM5hxUnRSdGfmmZl8XB57v+hJdrj1/QPYdz5LqiXwQqy+qIjUs88m5ZSl6HNjJ5EkPZLN7x5i89uH/L777r7nJ7P8qmkUTsqgvnfqSi9ZBSks+MZYNr97qPuYpdnO6udKOeO6GTGdlh1KnA43zdVWGiosNFS001Jvw251Yu9wYbe5sFudYY38zOkG0kaZSMs2kZ5ronhaFmk5RpLTDLF3bdclwfRz1ea0KZfuYJS+5f0g4ciWiASWY470wGm/hLKPoXyjEosOZmirvlRrSE1lsG0ljF4M3wsiDHCUkDBQA6GlAra+2Ofwvul/Z80D23B0uvxW5JZcMJHcMal88ORO/ymkdVWkZCQx+bjBqRI7Ol18/EwpB772f/KmZCVx5vWzgo7SBsKc04o5tL0BR6eaQnRJLaXTrmLqnhdJby9XlYQgecli0s89F9O8eWG5Vw8Uj1uy85Mqtq6upLk6cOTDxAWjWHHtDPSG0FOqUxfn03TEyoGtPd9lS10Hq58vZcU1MyJ2MonFFJ/0SNqbOmmts+F2e5BSHWtrsNFcbaWlzoZWr8Fo1mO3OtHoNOgNGizNdlrrbXS0RSd00drqwNrq6E7Dssm77mk06xk9I4uCiel01Escna4BO+eEhd4U/KHu8SiJIF8mnBq4rpSw8irInw1jlihP3EGquITEmKbWnE76ATis/XvV9RaQ7UftgtW/VfdTtFDVSxmZAs8JAxUpjWXw3AV++ZlcUs+GpPvYsUoDuJGi56k0d8VoRk9TjgSnXjWN9x7fQae1R09r06qDmFINFE8dmKRL5Z5m1vyzlLYG/8VOc7qB834wD1Nq/1NT4bqzS7cb246d2L7czJiaZvZn9MRzuDUGSqddRXHVOiZO0pH73e9gyM+P/GYi6JfL6ebwjibK3pOUtu8JWKd4WibHnzs+oiSLQghOuGACjk6X33pbY5WVD5/axbLLpwyJF2YXbpeHxioLtQfbqD/cTuMRK03VVr81xnij0+pk3xe17POOtJ/4eB2mVDX6TEpWSfs8Lg/JaQayi1LILk4hq9CM3qDF5fRgNOvIyEuOzijMbVcOB6VvQe12Jbg8rq/TDKDWiroEWAFMWfCTsqFzXzeYIXd68PKqLf77RQsC15NSeRZ3+Kg63vJ5/KhxREDCQEXC4c/g5WvAoqbRpIQy+xI22a6h2d53FDRzaSGzTinq3k/JNLL8yml89NROXC4V8yM9krUv7mHhWWOZekKP23aoB7Td5mLjq/vZ9UnfqY2UzCTOv2M+GbmDe/vzOBx0lu7G+sk6LBs24mlWD+xRgL3QTkVxz5uoFFoqik+lXqtn1iHJaFNs8kZ1WpxsfvcQuzYcwdkZ+CFtStVz8rcnM/m4vAE95LRaDSd/ezJrnt9D7cEekdfmGivvPb6DxedPoGhqZkymsdxOSdlXdVTvb1VGqaIdt9MTumGUSc02kpSsw9buxNpiR2gEadlGNDoNSElrgw2PK/zhoa3d6fcToLXeRnVZYBFdg1FLUrIeg0mLPkmHMUVP3vg0CiakYzDpSPb+bSWnGvqfIdCbYPnP1GapU/I/waR+yjf67xctCG6c6vconcChHJlcsVLljqr6Um3FQaYqW8r9jZPerKTWArH3g+j3M4okDFQ4ODpg9W/gs0cBSY1jKl93fIvD9vk4ZV8joNEKFp0zninH931A5hSncNKlk1n70p5uN2rpkWx+5xB15e0cd844kvsZ9UgpOfh1A5+8vBdLc1+32uziFL55yxxSs8J/y5cOB87aWhyVVd5MqVU4D1fQuXcP2ANPBxUf+QThcXF4zBl+xzutTja/c4jN7xwiu8hMVmEKJrMeY4oeU4pBPXCMWgxGHXqDFp1RizaEiriUktpDbexcW8W+L+uCPrCFRjBpwShOvLj/WLNwBo06vZZlV0xh9bOlfulGOq1OSv61h8y8ZCYuyGX83BySkge3NtVa38GhbY0c2t5A1V7Jbs+OQZ0vFEIjSM1KIrsohZzRqeQUp5CcZiApWYfBpCMpWdftYSqlxN7hQqP1j6FzOz20N3XS1mijraGTyt1NNFU1pzHVAAARjElEQVR30N7UGZXRnaPT3T2N3MWhbX29UPVGLcZkPRqtQKPTkJqZhFavQZ+kJTXLSEqWkdQsIxl5yaTljEJMOj34RXsbqP5UH1bdAeUbwJwLeTPgvIcgI7byXySlqNFfsBFgF72nAgvnq4wIgThQEpWuxQoRK82zowFTkl7a/nEezrKNtDpH0ebOZ1vHN6lyzA7aJi3HyMnXjseZGSAzpg97Pq/hi7f7umMbTDqmn1jA0qXzyB3lH2N0ZF8Ln762n5oDbQHOKJmSb2FOQQNaj0NJ/zidSIf3p9NBdUUFeVlZSKcTT6cdd2srrS21VLdURr5YIgQpy5djPfliNq9r6lY7HwjKW0yLPklDUrIeo1nPmLQxpKWaQQgqdjUGdNn3ZfqJBSw6d3xYhrmuDloCO/n1wWF38elrZVTsCpxRVqvTUDwtk6wCM+aMJPWwT0/qTkcyahRk9pq9dbs81B5qo7K0iQNf14e8t0AkJevIKlBK83hjxkxpBrIKzKSPMuG0u/G4PZhSDNg7nLicHtKyTaSNUg/tUC8FA8XjkdQcaKV6fwsNlRYO763D2R4f7vbmjCRMqXq0Og1p2UaSM5KUx6aEJLMOk6YdU8c+TG3bSWrYjPOk/8aRPg1jip70HBPGFD1CCNau/ohlG6/0jyH6eYVaT+rN/o/gi6cgrQDSCtW61vhTYnuj7TUqG3DVZmWsJpwKK34ZuO4LlyKu+k+HlDIy/bMhImGg+sFsEHLHnSv4oOVHuAjhUitg6vH5LL1sCi6tnar24G7mXVTsbmLjK/txBnjjzJSjmDClgOxsDZaaFlrrbTS0BR5ZJXfUMm3PC2S0Rh4lbtW5qTeF98YrkpMxL1lM8qJFmObPR5euNHzsHU52f1rDzg1VEU379EeWOw89/U8RmtIMmPIcnHbRceSND/BwCEIkBgrUCHf72kq2l1SG9aAVQq21mFIM5OZBXoEOU5oBIZQ+Yk1ZK64Ipu3M6QbyxqeTNz6NnOIUsgpTMGcMgadcFCgpKWHp0lOwtTkAgb1DTfEJjaC13kZjpYXGKgtN1VakR2Iw6Wg6Yg34PzHc6JO0pOUYcdnrGWXZgkF0kKSxok/Sol12B1qdBq1OoDfqyClOwZRqIOmrR9Gu8cn9uvA6+NZfA1/gw1+pGKi0AkjJU+rj0RSZDcQD8xA/3JowUMONECILeBKlitEA/LeU8l/9tTEbhPzb9f/B5unHgUHA2FnZnHjRJLIK1O/Y4rBwpL0ft1cfWus7WP+f/X080DI8OSTJEAGi0sOYio8Yf+gdtIESmYWBReemoR8Dpc3OxjRvLuaTTyZ54cJ+pX0szZ2U72ikorTJb1psIPRnoFKzjMxeXszsU4tYv+ETlkcorhqpgeqiqdrCttWVfrFgoUgxOjEbIxtdZuQlM3ZmNgWTlFEaSqeMaFNSUhLx78fjkdg7nDhsbpx2Fw6bi+aaDqr2ttBa14HL6cHaYsfjkUHXIeMJLQ50wo5e2DBnmUkdp0b6BqNOGTW9ijfTvf8jdI4G9KITrXAhLv8XwpyF0Ag0GoHQgEajQVf3JfqtT6NPNqFNNiPGLg6uC9jofWk1pCgnDIPZP4196SrEjG/FrYE6ltagHgYcQB4q2eLbQoitUsqdwRq4pT6occrLdDJtjI3cdAd6bQOs20OzlIDE4rbR4mxS8xoS7/yGRNpsuNvacVva8bS1q9TfbW1Ma7NwxDyNyuJTcWvCcyzIadjGhANvktJRHeHXEBhtTg76okL0hd5MqcVFGMaORV9QEPabekqmkZlLi5i5tIiOdgf15e3Y2h10Wp3YLE7sVieOThfOTjcOuwtnpwen3RX2w37MjCxmn1rMmJnZg8roO9B3sqyCFJZfOQ1rq4MDX9dRtqUu4DpgpAgBhZMzGDcnh9qOMs48b3HoRkcxGo3AlGLA5DN4KJycycylRX71pJR0tDlwOz143BKnw01rnQrjcNhctDd1YmnqpLXBRu2htqiN7iPFjQG3NGCXqVgaoLahLkjNG/x3Hz6IkkQNRJd2pAetRqI1rEWr1yCEQAg1QhVCICzV4O5EIBF40GSNQWdOQafXoNNr0GhjvG42SI4JAyWEMAMXA7OklBZgvRDiTdRv+efB2nnoO0+f1bSLwuqNjKr/CkHg9NxWnYdGU2RvzYXtnzOqfisNOXOoGzUfjIElWTJa9jHh4FtktB1Am5ONbsw0pcqQnY0mNRVhMKjcQgY9Qm9Q8kAGA8KgZ09ZGdNnz1FlBgPatHSsyVpSU+mjwj1YklMNjJ0VOuW0lBKX04PL4cZhc6lg0Q4XWZ5cpE2LzeJECBgzM3vArvjRxpxuYPayYmYtLaKuvI2maivWZjvNtR0013Xg6Aj9uzdnJFE0JYOiqZmMn5PTHQ5QUnLsai5HihCij0LKqNGBPfRcTjeWJjtOuxu7zUVbg41OixOtTv2Pd1qd2Nod2CxOOi1O7B1OdAYtBqOWjjYHrQ3Rcf6IPhrcHnB3uiHgaLLX/0y9B+oDrWHHJ8eEgQKmAG4p5V6fY1uBZZGcZMFXfx7QOk+46N2dFNRuoqB2E2ZZgDNlMq3Z07FkjiPb1MGMiS6yp43FMP7PGMaOidiodJaUkN5rusVlb0Nj6ZuqfKgQQqA3aNEbtJhSekaPY9ILMUYho3AsERrhXRvyz6fhdnnotKoHnVFnJ0ljx9buwOXwkJplpHBKBumjTCNiDeloQafXkpHn43Eb4cuOlBJbu5O2RhtffLqFqZOn4/Aqb7gcHtxOD26X2tqbOmmps3W/cMkAyiYJwuNYMVApQO+Ai1agz+uWEOJGVHJG9D7qw2mtB0gP2zgNQEdPgCctHU9WJq7MLIyFU0maOI2sggKyvAra5d6N6iNqixCLxUJJSYn/MZeFZkf8pfg6mHSQJG1orbdA9xSKpiY9VuvQ/OmnpTlJT3eBETCCzQV1u4LXH8j9xDNH2/3w/9u71xi5yjqO49/f7NRt6e5GWhARkAoWkcY0ARETBI03wBir1hsFwhswQnjFG/UFplFjMbzSeG3CTRBrTAAxampIQFFSY9PGC0mpVEWQtgSF0l3abbvz98XzLBxmZ6d0L3POmf19kpPOOWe2+/z3PzP/Oec853kAjYzxzFiePibntWj4lFc+WCIgJkTrCEwcThPeHh5L/8ZEGtewNZHu+29NQBxJ48pGK3+KxMtXCPKIIcDhI7QOB62JBq3o79H1F0qBGgXau3mNAFP6gkfERmAjwOCiwTi1+QxLGy9xxqlPsfjNa9IFg5cXQCp8E07bmzrEQY2m1fQkkGgsWUxjaIiBkREaw8MMDA+nx5NTgC965X6ak4dOZni6GwpnqNMF630H97F3bG/nHyjRaSOnsWTR0UcRn8lF+D174MUeneVYvjwtr9VM4qmyfosHqhVTqxXpyC1fh4tWpAF98/BXEfmm/0jPbU20OHIonVKfOJLW+WHZUUxvoRSonUBT0sqI+HvethqYtoMEgAbEmu9cecy/bP/4fmJ0bjovWL0tkE6yVpJGQzTyKfJ+1MM5kssTEWPAvcBXJS2VdCGwBrir3JaZmdl0FkSByq4HlgDPAj8BruvWxXw2ZjKXk72a/4ZmtlBO8RER/wOmmXJz4VqIhcCn3czqYSEdQdXKQiwc/cjF0GzmXKDMzKwrSZdKelzSE5KmHdxgrrlAmZnZtCQNkIaKuww4B7hc0jm9+N0L5hpULzXUeE03mXYzoN50Gx3QwKzbOh8amr/vTs0mDPYo5KbfYVZ/7wKeiIh/AEjaROoF3eWW87nht08X4+PjIelA2e2YQ01g5hM3VVO/xeR4qq/fYjpO0tbC+sY8YMGkU4CnCutPAxf0omEuUN1ti4h3lt2IuSJpaz/FA/0Xk+Opvn6L6TXE02nQyJ50//E1KDMz6+ZpoDgvx6nAsQ8GOgMuUGZm1s2fgJWS3iLpdcDngAd68Yt9iq+7jUd/Sq30WzzQfzE5nurrt5i6xhMRRyTdAGwGBoDb5msUnnYLZsp3MzOrF5/iMzOzSnKBMjOzSnKBMjOzSnKB6kDSMkn3SRqT9KSkdWW36VhJeljSQUmjeXm8sG9djmtM0v2SlpXZ1k4k3SBpq6RxSXe07fuApB2SXpL0kKTTC/sGJd0m6UVJeyTd2PPGdzBdPJJWSIpCnkYl3VTYX9V4BiXdml9H+yVtl3RZYX+tctQtnrrmCEDS3ZJ257btlHRNYV/1c5SmB/ZSXEjzRf0UGALeA+wDVpXdrmOM4WHgmg7bV5Gmur84x3cPsKns9nZo5ydJ06N8H7ijsP2EnI9PA4uBW4Athf0bgEeA44G3A3uASysczwrSTY/NaX6uqvEsBdbn9jeAj+bX1Yo65ugo8dQyR7ltq4DB/Pjs3Lbz6pKj0v+AVVvyC/UQcFZh213AzWW37RjjmK5AfQO4p7B+Zo53uOw2TxPH19s+0D8PPNqWrwPA2Xn9P8CHC/u/RoUKcId4jvbhV+l42tr6F2Bt3XPUIZ6+yBHwNmA38Jm65Min+KY6C5iIiJ2FbX8mfROpmw2SnpP0B0nvy9tWkeIBICJ2kQtyCe2bifb2jwG7gFWSjgfeVNxPfXL3pKSnJd0u6QSAOsUj6STSa+gx+iBHbfFMqmWOJH1P0kvADlKB+hU1yZEL1FRDpEPfon3AcAltmY0vAmeQBnrcCPxC0pnUP75u7R8qrLfvq6rngPOB00mnXoaBH+d9tYhH0iJSm++MiB3UPEcd4ql1jiLielJ7LgLuBcapSY48ksRUo8BI27YR0vno2oiIPxZW75R0OfAR6h9ft/aPFtYPtu2rpIgYBSZHkt6b79jfLWmEGsQjqUE6BX4IuCFvrm2OOsVT9xwBRMQE8HtJVwLXUZMc+Qhqqp1AU9LKwrbVvPpQv46CNCrxY6R4AJB0BjBIirsO2tu/lHQd7bGIeJ50CmN14fl1y93k0C6qejySBNwKnASsjYjDeVctc9Qlnna1yVEHTXIuqEOOyr5wV8UF2ETqybcUuJCa9eIDXg9cQuqd0wSuAMZIF0lXAS+SDveXAndTzQu6zdz+DaRvtJOxnJjzsTZv+yav7n10M/BbUu+js0lvtNJ7VHWJ54KclwawnNR79KGqx5Pb9gNgCzDUtr2uOZounlrmCHgDaWDXIdIYepfkz4E1dclR6S/yKi7AMuD+nMx/A+vKbtMxtv9E0gjE+4EX8pvuQ4X963JcY8DPgWVlt7lDDOtJ31SLy/q874OkC74HSL0VVxR+bhC4jVSE9wI3lh1Lt3iAy4F/5lzsBn4EvLEG8ZyeYzhIOiU0uVxRxxx1i6fGOToxF5kXctv+Clxb2F/5HHmwWDMzqyRfgzIzs0pygTIzs0pygTIzs0pygTIzs0pygTIzs0pygTIzs0pygTKbhTxP0NGWf+Xn3jH52MyOzvdBmc2CpHe3bbqPNPLz+sK28YjYngfrHYmI7b1qn1mdebBYs1mIiC3FdUnjwHPt2/Nzd/WsYWZ9wKf4zHqk/RRfYSrxL0jakKfW3p+n6T5O0lslbc5TjD8h6eoO/+dqSQ9Iel7SgTz310U9DcxsnrhAmZXvy6QJ4q4GvgJ8ljRw6X3AL4FPkGZ3vV3Sy5PGSToXeJQ0duS1pIE//ws8KOm8XgZgNh98is+sfLsiYvLoaHM+AroKuCoi7gaQtBX4GPApXpn24BbSoL/vj4hD+Xmbgb8BNwEf710IZnPPR1Bm5ft12/qO/O/myQ2R5uh5FjgNQNIS4L3Az4CWpKakJmnOrweBi+e70WbzzUdQZuV7vm39UJfti/PjZaQ5fm7KyxSSGhHRmqtGmvWaC5RZPb0AtIDvkuYnmsLFyerOBcqshiJiTNIjpKm4t7kYWT9ygTKrrxuB35E6VtxKmu31BOBcYCAivlRm48xmy50kzGoqIrYB55O6ln8b+A3wLeAdpMJlVmse6sjMzCrJR1BmZlZJLlBmZlZJLlBmZlZJLlBmZlZJLlBmZlZJLlBmZlZJLlBmZlZJLlBmZlZJ/weN6fiSRn6VSAAAAABJRU5ErkJggg==\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( statsLKP1[SIDX_EI,:], color = \"#ff7f0e\", linestyle = \":\", linewidth = 4 )\n",
    "plt.plot( statsLKP2[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",
    "plt.legend( loc = \"upper right\", fontsize = 10 )\n",
    "ax.set_xlim( 0, 300 )\n",
    "ax2 = ax.twinx()\n",
    "fillLimit0 = getLKP( [0,0,0] )\n",
    "plb.fill( np.arange( T ), fillLimit0, facecolor = \"green\", label = \"earlier lock-downs\", alpha = 0.15 )\n",
    "fillLimit1 = LKP\n",
    "fillLimit1[:100] = 0\n",
    "plb.fill( np.arange( T ), fillLimit1, facecolor = \"blue\", label = \"lock-down #3 (by ESOP)\", alpha = 0.15 )\n",
    "ax2.set_ylim([0,5])\n",
    "plt.ylabel( \"Lock-down level\", size = 16 )\n",
    "plt.legend( loc = \"center right\", fontsize = 10 )\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
}