SIR model implementation tryout

refactor_stratified.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy.integrate import solve_ivp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import itertools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## With stratification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Population is divided into groups (e.g. age groups), optionally with group-dependent parameters (in the example `beta` is stratified) and with interaction among them, defined by `nc`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### As simple and scipy-close as possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def integrate(t, y, parameters, nc):\n",
    "    \"\"\"Basic SIR model with interaction nc\"\"\"\n",
    "    \n",
    "    # unpacking need to be done explicitly\n",
    "    S, I, R = y.reshape(3, 2)\n",
    "    beta, gamma = parameters\n",
    "    \n",
    "    # Model equations\n",
    "    N = S + I + R\n",
    "    dS = nc @ (-beta*S*I/N)\n",
    "    dI = nc @ (beta*S*I/N) - gamma*I\n",
    "    dR = gamma*I\n",
    "\n",
    "    return np.array([dS, dI, dR]).flatten()\n",
    "\n",
    "# ... time, parameters and initial conditions\n",
    "time = [0, 150]\n",
    "parameters = {\"beta\": np.array([0.5, 0.4]), \"gamma\": 0.3}  # same order as definition\n",
    "initial_states = {\"S\": [440000, 350000], \"I\": [20, 10], \"R\": [0, 0]} # same order as definition\n",
    "nc = np.array([[0.9, 0.2], [0.8, 0.1]])\n",
    "\n",
    "#  -> runs model with scipy directly.\n",
    "output_ = solve_ivp(integrate, time, \n",
    "                   list(itertools.chain(*initial_states.values())), \n",
    "                   args=(list(parameters.values()), nc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7feebc1577d0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD4CAYAAADy46FuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3wc1bnw8d/ZolXv1ZLce5ElF7kBNjEuYIMxzSYU0wPkBkhyuaS9CSGQkNybxg0J+FJcQjAlELDBNq7YxMZyk3vvsmRbVrX6lvP+MSNZMrIkyyutyvP9fIaZPXPO7LNG0rMz58wZpbVGCCGEuByLrwMQQgjRtkmiEEII0SBJFEIIIRokiUIIIUSDJFEIIYRokM3XAXhbdHS07t69u6/DEEKIdmXr1q3ntdYx9e3rcImie/fubNmyxddhCCFEu6KUOnG5fXLpSQghRIMkUQghhGiQJAohhBANkkQhhBCiQZIohBBCNEgShRBCiAZJohBCCNGgDncfRXOVVbl4be0R44VSKGo2USgsytxWCotSKAUWBRZl1LRaFDarBXv12qqwWSzYrOqSbQuBflbCAuyEBdgJdthQStUflBBCtAGSKEzlVW7+d81hWvvxHDaLIjTATniA3VgHGgkk3EwkYYF+RATaSU0Op0d0kCQVIUSrk0Rhigp2cOw30+qUaa3RGrS57dGgMco8+uLa4wG31rg8HlxujcutcZrbTrcHl0fjcntwuo06pZVuisqrKCp3UljmNNblTorLneSXVnE0t5SicifFFc46iSsxPIBrekdzbd9oxvWKJiLIr3X/kYQQnZIkigYo8xKT+arV39/t0ZRUuDh3oYJNx/JZfyiXz3fn8N6WUygFg7uEcW2faK7pE83wbhE4bNZWj1EI0fGpjvYo1BEjRuiOPNeTy+1h5+ki1h88z1eHc9l+shCXRxNgtzKqZ6RxxtEnhr5xwXKZSgjRZEqprVrrEfXuk0TRvl2ocPL10Xy+OpTL+kPnOXq+FIDYEAc3DUngJzcNwM8mg9uEEA1rKFHIpad2LsTfzqSBcUwaGAfA6cJyvjqUy5cHc5m34Tg5ReX85dvDsFslWQghmkf+enQwieEBzBrZlb/eM5yfTx/I8j1n+cH7O3B7OtaZoxCi9cgZRQf20DU9qHR5+O2y/ThsFn53ewoWi/RbCCGujCSKDu6JCb2ocLr586pD+Nst/GrGYOnkFkJcEUkUncAzN/ShwuXm9S+P4rBZ+dm0AZIshBBNJomiE1BK8aOp/al0enjzq2P42y08O6W/r8MSQrQTkig6CaUUv7h5IJUuN6+uOYK/zcr3JvbxdVhCiHZAEkUnopTipVuHUOn08PsVB/G3W3n0up6+DksI0cZJouhkLBbF7+5IodLl4aXP9+GwW7h/THdfhyWEaMMkUXRCNquFP81OpdLl4eef7MFhszBrZFdfhyWEaKOafMOdUsqqlNqulFpivo5USq1QSh0y1xG16v5YKXVYKXVAKTWlVvlwpdQuc98ryhx6o5RyKKXeM8s3KaW612ozx3yPQ0qpOd740ALsVguv3pPG+L4x/OijXfxr+2lfhySEaKOu5M7sp4F9tV7/CFilte4DrDJfo5QaCMwGBgFTgb8qpaqnNf0b8BjQx1ymmuUPAwVa697AH4HfmseKBH4BjALSgV/UTkji6jhsVl6/bzije0Txg/cz+XxXjq9DEkK0QU1KFEqpJGAa8Eat4hnAfHN7PnBrrfJFWutKrfUx4DCQrpRKAEK11hu1MRPhgkvaVB/rQ2CiebYxBVihtc7XWhcAK7iYXIQX+NutvDFnBMO6RvDUu9tZufesr0MSQrQxTT2j+BPwX4CnVlmc1joHwFzHmuWJwKla9bLMskRz+9LyOm201i6gCIhq4Fh1KKUeU0ptUUptyc3NbeJHEtWCHDbeenAkA7uE8uQ721h3UP4NhRAXNZoolFLTgXNa661NPGZ9t/zqBsqb2+ZigdZztdYjtNYjYmJimhimqC3U386Ch9LpFRvMYwu3sO1kga9DEkK0EU05oxgH3KKUOg4sAr6llPo7cNa8nIS5PmfWzwKSa7VPArLN8qR6yuu0UUrZgDAgv4FjiRYQHujH3x9OJ8Tfzl/XHPZ1OEKINqLRRKG1/rHWOklr3R2jk3q11vpe4FOgehTSHOATc/tTYLY5kqkHRqd1hnl56oJSarTZ/3D/JW2qj3WH+R4aWA5MVkpFmJ3Yk80y0UKigh3MTEtk7YFc8koqfR2OEKINuJrnUbwMTFJKHQImma/RWu8B3gf2AsuA72qt3WabJzA6xA8DR4ClZvmbQJRS6jDwA8wRVFrrfOBXwGZzecEsEy1oZloiLo9myU4ZBSWEkEehisuY+qd1+Nut/Ou743wdihCiFTT0KFR5wp2o18y0RDJPFXI0t8TXoQghfEwShajXjNRElELu2BZCSKIQ9YsP82dsryg+zjxNR7s8KYS4MpIoxGXNTEviVH45W0/IPRVCdGaSKMRlTR0cj7/dwkdy+UmITk0ShbisYIeNKYPi+WxnDpUud+MNhBAdkiQK0aBb0xIpKneyZr/M/yREZyWJQjTo2t7RRAc7+Hh7VuOVhRAdkiQK0SCb1cItQ7uwev85CsuqfB2OEMIHJFGIRs1MS8Tp1nwmDzYSolOSRCEaNTgxlN6xwXy8TUY/CdEZSaIQjVJKMTMtkS0nCjiZV+brcIQQrUwShWiSGaldAPhXppxVCNHZSKIQTZIUEcioHpF8vF2m9BCis5FEIZrstmGJHDtfSuapQl+HIoRoRZIoqnk8UJZfdykvgPJCY6kogsoLUFkCVWXgLAdXJbiqwO002ndwUwcn4GezyIyyQnQyNl8H0GaU58N/97q6YygLWGy1Fms9r+1g9wdHGDhCwD8UHKGXbIea2yHGdkA4hHQBi2/zeliAnUkD4li8M4efTR+I3SrfM4ToDCRRVPMLght/BzXX33U92xq0x9jWnm++9rhqLe6L225nrddOcFZAZTEUZcE580ylohh0A/MpBURAUjokp0PyKEgcZsTcymamJfLZrhy+PJDLDQPjWv39hRCtTxJFNXsAjPqO795fa3CWXUwaleZSUQxleZC9HU5lwKHlRn1lhfgh0HX0xeQRltTiYY7vF0NEoJ2PM09LohCik5BE0VYoZZwh+AVBSPzl65XlQ9YWOLXJWLYtgE2vGftCEy8mjeR0iE8Bq92rYdqtFm4e2oVFm09RXOEk1N+7xxdCtD2SKNqbwEjoO9lYwLisdXa3cbZxapOx3vOxsS84Dmb/A5LqfV56s81MS2TBxhMs3ZXDrJFdvXpsIUTbI72R7Z3VDl3SjMtmd7wF398N399rbNv8Yd502LfYq2+ZmhxOj+ggPpIpPYToFCRRdERhiTD4dnhkFcQNhPfug41/9drhlVLcmprIpmP5nC4s99pxhRBtkySKjiw4BuYsgf7TYPmPYelzxugrL5iZlggg91QI0QlIoujo/ALhrgUw+rtGp/d790JV6VUftmtUICO6RciUHkJ0ApIoOgOLFab+2rhP5OAyo9+i5NxVH/bWtEQOnythT3axF4IUQrRVkig6k1HfgVnvwLl98MZEyD1wVYebnpKAn9UindpCdHCSKDqb/jfBg58Zd4e/OQmOrW/2ocID/bi+fwyf7sjG5e74c10J0VlJouiMEofDIyshOB4WzoQd7zX7UDPTEjlfUslXh897MUAhRFvSKW64czqdZGVlUVFR4etQms3f35+kpCTsdi/dCR3RDR5ebgyd/fgxKDwB1z1r3CF+Ba7vH0tYgJ2Pt59mQr9Y78QmhGhTOkWiyMrKIiQkhO7du6Ou8A9hW6C1Ji8vj6ysLHr06OG9AwdEwL0fwaffgzUvGcli+p+uaNoPh83KtJQEPtqWRUmli2BHp/iREqJT6RSXnioqKoiKimqXSQKMG9yioqJa5ozI5gczX4Pxz8H2v8M7dxjP3rgCt6UlUuH0sHz3Ge/HJ4TwuU6RKIB2mySqtWj8SsH1P4EZf4XjX8FbU6E4u8nNh3eLIDkygI/l5jshOqROkyjagpdeeolBgwaRkpJCamoqmzZt4i9/+Qu9e/dGKcX58z7uEE67B+79JxQchxU/b3IzpRQzUxP595HznClqv/1AQoj6SaJoJRs3bmTJkiVs27aNnTt3snLlSpKTkxk3bhwrV66kW7duvg7R0HMCDJsDe/4FF5p+KenWtES0hk93yFmFEB2N9Dy2kpycHKKjo3E4HABER0cD0KVLF1+GVb/0R43pPra8ZVySaoKeMcEMTQ7no22neey6q3ykbCekPZrKcheVZS4qy5zm2oXb5cHt8uBx64vbLo3bba5dHjwuD263NtZmXY9Hmw9l1Mbj3LVG19k23lPX2TZjMTeqH+oIF/ddLK9V59LPUqdyPZvtYMqXdhBivaKTgrnpiRSvH7fRRKGU8gfWAQ6z/oda618opSKB94DuwHHgLq11gdnmx8DDgBt4Smu93CwfDswDAoDPgae11lop5QAWAMOBPGCW1vq42WYO8DMznBe11vOv5gP/cvEe9np5yomBXUL5xc2DGqwzefJkXnjhBfr27csNN9zArFmzGD9+vFfj8JqoXtB3ipEorv0h2BxNanZbWiK/+HQP+3KKGZAQ2sJBtn1VFS7ys0vJzyml/EJVzR//2omgZrvcVeePalNYbAqr1WKsbZaabYvVYjxeXSksFoVSoKrXSl3cthn1lLEDZV5fUFAzTLq6a8yoY+67WMnYX6v77OL++vvUaoob6HJTDe1sTW0kjCsRFhPQIsdtyhlFJfAtrXWJUsoOfKWUWgrcBqzSWr+slPoR8CPgOaXUQGA2MAjoAqxUSvXVWruBvwGPAV9jJIqpwFKMpFKgte6tlJoN/BaYZSajXwAjMH6NtiqlPq1OSO1JcHAwW7duZf369axZs4ZZs2bx8ssv88ADD/g6tPqNehwW3gq7P4LUu5vUZHpKAr9aspePt5/uVInC49EU55aTd7qE86dLyMsqIe90CcXn6/bXWGwKR6Ad/0AbjkAbgaF+RMQH4gi04zDLqrf9g2z4Bdix2WslApsFi1XVrNv7AA3RfjSaKLRxHllivrSbiwZmABPM8vnAWuA5s3yR1roSOKaUOgykK6WOA6Fa640ASqkFwK0YiWIG8Lx5rA+Bvyjjt2AKsEJrnW+2WYGRXN5t7gdu7Jt/S7JarUyYMIEJEyYwZMgQ5s+f33YTRc8JENMfNv0Nhs5u0o14UcEOxveN4ZPM0zw3tT9WS8f7Q1ZR4ryYEMykkJ9TiqvKmMJEKQiLDSSmaygDxiYQlRhMZJdgAsP8sNkt8sddtEtN6qNQSlmBrUBv4FWt9SalVJzWOgdAa52jlKq+LTcR44yhWpZZ5jS3Ly2vbnPKPJZLKVUERNUur6dN7fgewzhToWvXtvlozgMHDmCxWOjTpw8AmZmZbacDuz5KGZMILvm+8YjVrqOb1GzmsERW7T/HxiN5XNMnuoWDbB1ut4d9/85h+4qTFOdefFCTf7CdqMRgBl2TSFRSkJEUEoKw+Vl9GK0Q3tekRGFeNkpVSoUDHyulBjdQvb6vTLqB8ua2qR3fXGAuwIgRI9pkN1RJSQnf+973KCwsxGaz0bt3b+bOncsrr7zC7373O86cOUNKSgo33XQTb7zxhq/DNaTMgpXPGx3bTUwUNwyII8Rh46PtWe0+UXg8mkMZZ8hYcozi8xXE9wxl8LUXk0JgqJ+cIYhO4YpGPWmtC5VSazEu/5xVSiWYZxMJQPUDDrKA5FrNkoBsszypnvLabbKUUjYgDMg3yydc0mbtlcTcVgwfPpwNGzZ8o/ypp57iqaee8kFETeAXZAyV3fgqFGVBWFKjTfztVm4cEs9nO3Mou9VFoF/7G1intebo9lw2LT5GQU4p0cnBTPtuCt0Gt9+7+4W4Go3eR6GUijHPJFBKBQA3APuBT4E5ZrU5wCfm9qfAbKWUQynVA+gDZJiXqS4opUab/Q/3X9Km+lh3AKvNvpHlwGSlVIRSKgKYbJaJ1jLyEUDD5jeb3GRmWhKlVW5W7D3bcnG1AK01J3bn8cFvtrBs7m7QmimPDuauH4+k+5BoSRKi02rK170EYL7ZT2EB3tdaL1FKbQTeV0o9DJwE7gTQWu9RSr0P7AVcwHfNS1cAT3BxeOxScwF4E1hodnznY4yaQmudr5T6FbDZrPdCdce2aCUR3aDfTbB1Hoz/L7A3PvxuVI9IuoT589G208xI/UaXUpuUfaiArz85Ss7hIkKi/Jk4ZwB9R8Vj6YAd8kJcqaaMetoJpNVTngdMvEybl4CX6infAnyjf0NrXYGZaOrZ9xbwVmNxihY0+gnYvwR2fQDD7m+0usWimJGWyOtfHiH3QiUxIU27D8MXzh4vZtOnRzm1N5+gMD/G392XAeO6YLXJpAVCVJPfBtG4buMgbjB8/VqTb1m9LS0Rj4Z/tdGJAvNOl/D533by4ctbyD1xgbG39+beX41h8PgkSRJCXKL99TSK1qeUcQPep/9hzC7b49pGm/SJCyGtaziLNp/kkWt7tJnr+4XnyshYfIxDW87i57CSfnMPhk5Mxs9ffhWEuBz57RBNM+QOY0bZTa81KVEAzB6ZzHP/3MWWEwWM7B7ZwgE2rLSokozFx9i3IQerTTFscjfSJnfFP8hLTwwUogOTc+xWVN804/fccw/9+vVj8ODBPPTQQzidTl+HWT97AIx4EA58bkxD3gTTU7oQ5GdlUcapxiu3ILfLw+JXdrD/6xwGj0/k3l+NYczMXpIkhGgiSRSt5HLTjN9zzz3s37+fXbt2UV5e3nZutqvPiIcBBRn/16TqQQ4bt6Qm8tmubIrKfZcAt3x+nLzTJUx9bAjXzepLUFjb7VwXoi2SRNFK6ptmvEuXLtx0003GjJ5KkZ6eTlZWViNH8qGwRBg4A7YthMqSxusDd6cnU+H08Gmmbzq1c09eYOuyE/QfHU+PlPZ9p7gQvtL5+iiW/gjO7PLuMeOHwI0vN1ilsWnGnU4nCxcu5M9//rN3Y/O2UY/Dno9g5yLzZryGDUkMY0BCKIs2n+K+Md1bPr5a3E4PK+ftJTDEzjV39WnV9xaiI5EzilZSPc343LlziYmJYdasWcybN69m/5NPPsl1113Htdc2raPYZ5LToUsabHq9SUNllVLcnZ7MnuxidmUVtUKAF23+7Bj52aVMuLc/jkDpjxCiuTrfGUUj3/xb0uWmGf/lL39Jbm4ur7/+us9ia7LqobIffweOrIbe9d5zWceM1ERe+mwf724+yZCkIa0QpHEj3bYvTjJgbALdh8glJyGuhpxRtJIDBw5w6NChmtfV04y/8cYbLF++nHfffReLpZ387xg0E4JijbOKJggLsDNtSAKfZmZTVuVq4eDA5XSzav4+gsL8GHenXHIS4mq1k79M7V9JSQlz5sxh4MCBpKSksHfvXp5//nkef/xxzp49y5gxY0hNTeWFF17wdaiNszlgxENwaDnkHWlSk9npXSmpdLFkZ04LBweblxizvl5/b38cAZ3vpFkIb5PfolZyuWnGXa6W/4bdIkY8BOt/Dxlz4cbfNlp9ZPcIesYEsSjjJHeNSG60fnOdOVrE9i9OMvCaLnQdFNVi7yNEZyKJQjRPSBwMvg22vwPX/xT8G35GtlKK2SOT+fXn+zl49gJ940K8HpKrys3qBfsIinAw7vbeXj9+S/F43FSVl1NVVkZleZm5LqWqvBy304nH7cbjduF2ufC4XLjdbjwul1HmduNxOXG7LtbRbjdaa2PxeMxtD9qja7aps8/YBm2MT9AaXf18MG1Mv179vLDq/cYuXXdAwyVjG3Ttglr1mjhdWD3a5DPJ2pTo5G5MffL7Xj+uJArRfKO+Azvfg8x/wOjHG61++7Ak/nv5ARZlnOLnNw/0ejibFh+j4EwZtzydip8PLzlprSktyCc/+zQFOVkUnMmhouQCVeVlVJWXU1lWSlVZGVXlZVSWleGsrGjeGymF1WbDYrVhtVqx2GzGYrFisVhQFgXKUnOfjrJUbxv7qrexKBTKqG8c2JibSwEKLMoCStV6bLqqeYb6pXN4fWNOr1qv1WXKr+wjt405w9oq/2DvfwEDSRTiaiQOh6R0yHgd0h+DRjrjo4IdTB4Yz0fbs/ivqf3wt3vv2dI5R4rIXHmSQdclkjygdeaVqiovoyAnm/yc0xRkZxnb5tpZcfHZ2ja7H/6hoTgCAvELCMA/KJjQ6Fj8AgJxBAZeXAcGmnUullltdiw2q7G2Wo2kYLNhsVmxWOTZ3KJ1SKIQV2fUd+CfD8PhFdB3SqPVZ6cn89muHJbvOeO1hxo5q9ysmr+XkEh/xt7WyyvHvFRFSQn7N6wj98RRCnKyKcjOoqSg1jO0lCI0OpbILokk9h9IREIikQlJRHRJJCQyCtVeRrQJUQ9JFOLqDJwBX/zMmFW2CYliXK9okiICeG/zKa8lik2fHKXoXDkzvp/m9enC87NPs33Zp+xZuwpnZQX+QcFEdEmkW0oaEQmJRHRJJDIhkfD4Ltj8/Lz63kK0FZIoxNWx2mHkw7D6Rcg9ADH9GqxusShmjUjm9ysOciKvlG5RQVf19tmHCtmx+hRDxieS1C/iqo5VTWvNyV072Lb0E45u24zVZqP/uAkMu+kWYrv39Mp7CNGeyPlwKwkODq7Znjp1KuHh4UyfPt2HEXnR8AfB6mjyDXh3jkjGomDR5qubftxZ6WbVgn2ERvkzeubVX3JyVlWya/UXLHj2P/jwpZ9x5sghxtzxbR599W2mPvmMJAnRackZhQ88++yzlJWVtY8pO5oiKBqG3Ak73oWJ/w8CGv5mHx/mz/X9YvlwaxY/mNQXu7V531c2/usIxbnl3PqDq7vkVFKQz44vPmPHiqWUXygmplsPpj75ffqNvQ6bXeaIEkIShQ9MnDiRtWvX+joM7xr1Hcj8O2z/O4z9XqPVZ6d3ZdWCLazef44pg+Kv+O1OHyxg15osUq5PIrFv8y45nT16mG2ff8L+DevxeNz0Gj6K4TfdQtLAITIMU4haOl2i+G3Gb9mfv9+rx+wf2Z/n0p/z6jHbnYQU6DbOuFN79JPQyNDN6/vFEBfqYFHGyStOFFUVLlYv2EdYTACjb72yS04ej5sjmzex9fNPOL1/D3b/AFIn30Ta1JsJj0+4omMJ0Vl0ukQhWtCo78D798OBpTCg4f4Xm9XCncOT+evaw2QXltMlPKDJb7Px4yMU51Uw84fDsDuadi9BVXkZO1ctZ/uyJRTnniU0Jo4J9z/C4Osn4Qi8ug51ITq6TpcoOv03/5bUbxqEJRtDZRtJFAB3jUjmL2sO88GWLJ6+oWmzvGbtz2f3l6cZOjGZLr3Dm9TGVVXFB7/6KWeOHCJpwGAm3P8wvUaMkhvWhGgiGfUkvMdqM556d3w9nNndaPWuUYFc0zua97ecwu1pfB4f45LTfsLjAhk1o2kjkLTWrHzjr5w5coibv/8jZj3/Mn3Sx0qSEOIKSKLwgWuvvZY777yTVatWkZSUxPLly30dkvcMux9sAca0Hk0wOz2Z04XlfHX4fKN1N/zzMCUFFUycMwC7X9P+0O/44nP2fLmS0bffTd/R1zSpjRCirk536clXSkpKarbXr1/vw0haWGAkDJ0FOxbBxOchqOGpvicNjCMi0M6ijJOM7xtz2Xqn9uazZ302qZO6Et8zrEmhZO3fw5r5c+k5bCRj77j7Sj6FEKIWOaMQ3jfqcXBVwLb5jVZ12KzcPiyJFXvPknuhst46VeUuVi/cR0R8IKNu7tGkEC7kn2fxH35DWGwcN/7HD2WuJSGugvz2CO+LHQA9xsPmN8DtbLT67PRkXB7NP7dl1bv/3x8eorSwkm/NGYCtCZecXE4ni3//G5yVlcz4z5/hHxTcaBshxOVJohAtY/QTUHwa9n3aaNXesSGM6BbBe5tPmQ/JuejU/nz2/juHtMldie/RtEtOq99+jZzDB7jxye8TldS1WeELIS6SRCFaRp/JENUb/v3nJj3SbHZ6V46dL2XTsfw65bvWZBEY6sfI6U275LRjxVJ2rVrOqJl30WfU2GaFLoSoSxKFaBkWK4x7GnJ2wNE1jVafNiSBEH8bizJO1pRVlDo5sTuPPulx2JrwkKPTB/ax+u3X6Z46nLF33XNV4QshLpJEIVpOyiwISYCv/tho1QA/K7emJvL57jMUlRn9Goe3nsPj1vRLb3yKj5L8PBb/4deEREcz7XvPyn0SQniRJIpWUj3NeGZmJmPGjGHQoEGkpKTw3nvv+TiyFmRzGPM+HVsHp7c2Wn3WyGSqXB4+3m50ah/MOENEfCDRyQ13RrtdTj7942+oKi83Oq+DpfNaCG+SRNHKAgMDWbBgAXv27GHZsmU888wzFBYW+jqsljP8AfAPg6/+1GjVwYlhDEkMY9HmUxSfLyfncBF90+Mbncl1zby55Bzcz5QnniGma3fvxC2EqCGJopX17duXPn2MeY26dOlCbGwsubm5Po6qBfmHwshHYd9iOH+o0eqz05PZf+YC61aeAKBvelyD9XeuWs6OFUsZOeMO+o2RO6+FaAmN3pmtlEoGFgDxgAeYq7X+s1IqEngP6A4cB+7SWheYbX4MPAy4gae01svN8uHAPCAA+Bx4WmutlVIO8z2GA3nALK31cbPNHOBnZjgvaq0bv4urAWd+/Wsq93l3mnHHgP7E/+QnV9wuIyODqqoqevW6+qeztWmjHoeNfzFGQM34S4NVbxnahRcX7+Pw5rN07RVGaPTlZ5XNOXSA1W/9jW4paVwz+z5vRy2EMDXljMIF/FBrPQAYDXxXKTUQ+BGwSmvdB1hlvsbcNxsYBEwF/qqUqu5Z/BvwGNDHXKaa5Q8DBVrr3sAfgd+ax4oEfgGMAtKBXyilvPNgZB/Lycnhvvvu4+2338bS0e8aDo6BtPuMaT2KsxusGuJv57aeMdhL3XQddvkpPUoLC/j0D78mODKKaU//l3ReC9GCGj2j0FrnADnm9gWl1D4gEZgBTDCrzQfWAs+Z5Yu01pXAMaXUYSBdKXUcCNVabwRQSi0AbgWWmm2eN4/1IfAXZVyYngKs0Frnm21WYCSXd5v7gZvzzd/biouLmTZtGt++oisAACAASURBVC+++CKjR4/2dTitY+x/wJa34Ou/wuQXG6w6XDnIQnPQ7mZEPfvdLieL//gbKkpKuPtX/01AcEjLxCyEAK6wj0Ip1R1IAzYBcWYSqU4msWa1ROBUrWZZZlmiuX1peZ02WmsXUARENXCsS+N6TCm1RSm1pa1f76+qqmLmzJncf//93Hnnnb4Op/VEdIfBt8GWt6G84LLVPB5N8f5CzgUp3ttV/9nH2gVvcnr/XiY//hSx3Zs23bgQovmanCiUUsHAP4FntNbFDVWtp0w3UN7cNhcLtJ6rtR6htR4RE3P5yxVtwfvvv8+6deuYN28eqamppKamkpmZ6euwWse4Z6CqxJgD6jKyDxZQWlRF12Ex7DhVyL6cuj9qu9euJHP5EoZPn8mAceNbOmIhBE1MFEopO0aSeEdr/ZFZfFYplWDuTwDOmeVZQHKt5klAtlmeVE95nTZKKRsQBuQ3cKx2p3qa8XvvvRen00lmZmbNkpqa6uPoWkn8YGNqj69fg6qyeqscyDiL3d/KzGm98bNa6typfebIIVa+8SpdB6dw3bcfaKWghRCNJgqzr+BNYJ/W+g+1dn0KzDG35wCf1CqfrZRyKKV6YHRaZ5iXpy4opUabx7z/kjbVx7oDWK2N2eGWA5OVUhFmJ/Zks0y0V+OegbLzkPnON3a5nG6ObjtHr9QYYsIDmDo4no+3n6bC6aasqJBPf/9rgsIjmPb0c1is0nktRGtpyoOLxgH3AbuUUtXXSH4CvAy8r5R6GDgJ3Amgtd6jlHof2IsxYuq7Wmu32e4JLg6PXWouYCSihWbHdz7GqCm01vlKqV8Bm816L1R3bIt2qttYSEqHDa/A8AeNx6eaju/Mo6rCTd9RxpQds0cm8+mObD7fkYXzs79RXlzE7Bd+R2Bo02aRFUJ4R1NGPX1F/X0FABMv0+Yl4KV6yrcAg+spr8BMNPXsewt4q7E4RTuhFFzzfVh0N+z5GFIu/m8/mHGGwDA/EvsZI6BH94yiW1QgG96dR0LWbm787g+I69nbV5EL0Wl18AH8ok3qOxVi+huTBZpTkNfMFDsyDovF+F5isShmxRSRkLUZBl3LgGuv92XUQnRakihE67NYjL6Kc3vg0Arg8jPFhh5chys4mldLB/Cbpfu/8WAjIUTLk0QhfGPIHRCaVDMFeX0zxeYcPsDZIweZdMft3Du2J3PXHeXnn+zB45FkIURrkkTRSqxWK6mpqQwePJibb765zoyxU6dOJTw8nOnTp/swwlZmtRt3a5/cQPGujfXOFJu5bAl+AQEMnjCRX94yiMeu68nCr0/w3D934pZkIUSrkUTRSgICAsjMzGT37t1ERkby6quv1ux79tlnWbhwoQ+j85Fh90NABIc+/xKoO1NsWVEhBzauZ+B1E/ELCEQpxY9v7M/TE/vwwdYsnnkvE6fb46vIhehUJFH4wJgxYzh9+nTN64kTJxIS0gnnK/ILQqc/zsFT8SR0tdeZKXbnymW4XS7Spl48y1JK8f1JfXluan8W78jmu+9so9Llru/IQggvasp9FB3K+vcPcv5UiVePGZ0czLV39W1SXbfbzapVq3j44Ye9GkN7lZd0H/muA4wP+hq4FgC3y8WOFZ/TLSWNyC5J32jzxIReBNgtPL94L48t2Mpr9w4nwE9uwBOipcgZRSspLy8nNTWVqKgo8vPzmTRpkq9DahMO7CzHojz0yn8VCo35Hw9v3khJQT5pU2++bLsHxvXg5duGsO5QLg/Oy6C00tVaIQvR6XS6M4qmfvP3tuo+iqKiIqZPn86rr77KU0895ZNY2gqPR3No81m69g8hoOgCbHwVbnyZ7csWExYXT4+04Q22n53eFX+7lR9+sIP73tzE2w+mExZgb6Xoheg85IyilYWFhfHKK6/wP//zPzidTl+H41PZBwsoLayk77juMOQu2Dafc/u2cXr/XlInT2vSw4huTUvk1W+nset0Efe88TX5pVUtH7gQnYwkCh9IS0tj6NChLFq0CIBrr72WO++8k1WrVpGUlMTy5Z1j3sOD5kyx3VOiYdzT4Cxj+6K/YHM4GDyh6Zfmpg5OYO59Izh4toS7537NuQsVLRi1EJ1Pp7v05CvV04xXW7x4cc32+vXrWzscn3M53RwxZ4q1+1khtj/lPW5k/9IzDJwwCf/g4MYPUsv1/WN5+4GRPDJ/C7Nf/5p3Hh1FQtjln7cthGg6OaMQPlEzU2ytKTt2qVG4tIXUpOZdPhrXO5qFD6dz7kIld762kVP59T/zQghxZSRRCJ84mHGGwFA/EvsbM8V6PG52ZOwiOVITc3ABuJqXLEZ0j+SdR0ZxocLFna9t5Eiud4dCC9EZSaIQra6+mWKPbM2gOPccaVOmQ3EW7P6w2ccfmhzOosdG4/J4mPX6RvafaejJvUKIxkiiEK3uyDZzpthRFy87ZS5bQkhUDL2mPwpxg+GrP4Gn+VN0DEgIZdFjY7BaFLPnfs2urCJvhC5EpySJQrS6gxln68wUm5d1kpO7dzB00o1YbDZjCvLzB+Dg0kaO1LDescG8/50xBPnZ+Pb/fc36Q7kyTbkQzSCJQrSq4rxysg8V1pkpdvuyJVjtdoZMnGJUGjQTwrvWebBRc3WLCuKDx8cQFezHfW9mMOrXq3j2gx18tjOHovLOfR+LEE0lw2NbidVqZciQIbhcLnr06MHChQsJDw8nMzOTJ554guLiYqxWKz/96U+ZNWuWr8NtMYc2nwUuzhRbWVbK3nWr6T92/MVnYVttMPYp+Pw/4cQG6D7uqt6zS3gAn37vGr7Yc5a1B86xfM8ZPtiahdWiGNY1nAn9YhnfN4ZBXULrTHMuhDBIomgl1VN4AMyZM4dXX32Vn/70pwQGBrJgwQL69OlDdnY2w4cPZ8qUKYSHh/s4Yu/TWnMw4yzxPcNqZordvWYlzsqKOrPEApB2L6x92TiruMpEARDqb+eO4UncMTwJl9tD5qlC1h7IZe3Bc/z38gP89/IDxIY4GN83hgn9YrmmT7RMByKESRKFD4wZM4adO3cC0LfvxbmnunTpQmxsLLm5uR0yUeSdLiE/u5TxdxufWXs8ZH6xhC59BxDXs3fdyvYAGP04rH4RzuyC+CFei8NmtTCieyQjukfyn1P6ce5CBesOnmeNnG0IUa9OlyjWzJvLuRNHvXrM2G49uf6Bx5pUt6FpxjMyMqiqqqJXr15eja+tOLjpLBaLotfwWACO79hG4Zkcxt11b/0NRj5ijH7695/h9jdaLK7YEP9GzzZiQhxc1yeGXrFBxAQ7iAm5uEQFObBaJImIjqvTJQpfqZ5m/Pjx4wwfPvwb04zn5ORw3333MX/+fCyWjjfGwOPRHNx8lq6DowgI9gNg+7LFBIVH0GfU2PobBUTAiIdgwysQ3Q+u+09o4W/1lzvbWHvgHKv3n+Wf277ZAW5REBlUK3lckkiig/0I8rPhb7fib7fgb7fisF1cy5mKaOs6XaJo6jd/b2tomvHi4mKmTZvGiy++yOjRo30SX0vLPlRIaWEl4+4wLjEV5JzmWOZWxtzxbay2BvoCrv8JlJyFNS/CmZ1w69/AcWXzQF2N2mcbAOVVbs6XVHLuQgW5FyovLiUXtw+fvUBuSSVOd9NGbDlslprEUZ1MHDZjbbUorBaFRSlstbatFoXFYpYpY9uqFFarwqJAYa7NJKTqlBnlCuCScqPoYhuzCrULVN2XNfXrc7kc2NzUKDm1YbGh/tw1Itnrx+10icLXqqcZnzFjBk888QRaa2bOnMn999/PnXfe6evwWszBjDPYHeZMsUDm8s+wWG0MnXRjww3tATDzdYhPgRX/D948ArPfgcgerRD1NwX4WUmODCQ5MrDBelprisqdNUmkvMpNhdNDhdNNpctYV7iMssraZU6zzNzn9micbmNdvXi0ua1rldW8Bo/WaK3xaCMODaBBm3F5NGg02ixDG23MzZr46772/r+l8L7U5HBJFB1F7WnGlVKsW7eOvLw85s2bB8C8efNITU31bZBe5HK6ObL1HL3SjJliqyrK2b12JX1HjyMoPKLxAygFY/8D4gbCBw/C/10Pd7wNva5v+eCbSSlFeKAf4YF+9InreM9Dr0kkDSSQy+1q7k2Pkqsa11InXJIoWklD04zfe+9lOnM7iBO76s4Uu/fL1VSVlzX4qNN69foWPLYG3v02/P02mPwijH5Srkf4QO1LWs1o7dVYRMvreL2mos05sOniTLFaa7YvX0Jczz4k9Ol35QeL7AmPrIB+N8Hyn8C/ngBnufeDFkLUkEQhWlRFqZMTey7OFHty9w7yT58iber05o/2cYTAXQvh+p/Cjnfh7Ruh6LR3AxdC1JBEIVrUkW3n8Lh0zZQd25ctISAklH5jrr26A1ssMP6/YPY/4PwhmDsBTn599QELIb6h0ySK9j5raHuNv3qm2JiuIRSdO8vRrRmk3DAVm5+fd96g/zR4ZJUxZHbedNg6zzvHFULU6BSJwt/fn7y8vHb7x1ZrTV5eHv7+/r4O5YpcyK8wZ4qNQylF5hefgYKhk27y7hvF9odHV0OP62Dx07DkB81+Qp4Q4ps6xainpKQksrKyyM3N9XUozebv709SUpKvw7giBzPOANBnZDzOygp2r/6CPiPHEBIV7f03C4iAez6AVb80pvw4tw/uWgDBMd5/LyE6mU6RKOx2Oz16+OYGrc6seqbYsJgAdq3+gorSElIvnSXWmyxWmPSCcXPeJ981+i1mvwNdOs49KUL4Qqe49CRa35mjReRnl9JvVJwxJHbZYqK7didpwOCWf/Mhd8BDy41B/m9NgZ0ftPx7CtGBSaIQLWLL0uP4B9vpNzqB0/v3kHvi2NUNib1SXVLh0TWQOBw+egS++Bl43K3z3kJ0MI0mCqXUW0qpc0qp3bXKIpVSK5RSh8x1RK19P1ZKHVZKHVBKTalVPlwptcvc94oy/2IopRxKqffM8k1Kqe612swx3+OQUmqOtz60aFm5py5wYlceQ7+VjN1hZfuyJfgHBTPgmgmtG0hwDNz/CYx8FDb8L7w+HlY+D4dWQEVR68YiRDvWlDOKecDUS8p+BKzSWvcBVpmvUUoNBGYDg8w2f1VKWc02fwMeA/qYS/UxHwYKtNa9gT8CvzWPFQn8AhgFpAO/qJ2QRNu1bdkJ/PytDJmQyIW88xzK2MDgb03G7vDBqC2rHab9D9z6mjHB4Ib/hXfugN92h9euhaXPwd5PoKT9DnQQoqU12pmttV5X+1u+aQYwwdyeD6wFnjPLF2mtK4FjSqnDQLpS6jgQqrXeCKCUWgDcCiw12zxvHutD4C/m2cYUYIXWOt9sswIjubx75R9TtJaCM6Uc3naO4VO64Qi0s2XxUrTWpE728pDYK5V6t7FUlUHWZji5EU78G7bOh02vGXWi+kC3sdBtHHQbA+FdfRuzEG1Ec0c9xWmtcwC01jlKqVizPBGofXtsllnmNLcvLa9uc8o8lkspVQRE1S6vp00dSqnHMM5W6NpVfrl9advyE9hsFoZOTMbldLJz1XJ6DhtJWGy8r0Mz+AVCz/HGAsb9Fjk7jKRxYgPs+Rdsm2/sC0s2EkfXMUbyiO4jExCKTsnbw2Pr+y3SDZQ3t03dQq3nAnMBRowY0T7vqusAivPKObjpLIMnJBIQ4sfedaspKyq88lliW5PND5JHGss1zxgd3uf2GknjxAY4sgZ2vmfUDYwyzjICo8wlGgIjL74Oir64HRBhDNcVogNobqI4q5RKMM8mEoBzZnkWUPupGUlAtlmeVE957TZZSikbEAbkm+UTLmmztpnxilaQ+cVJUJA2yTir275sMZFdkug2pB3dx2CxQvwQYxn1HeOBC/lHjTOOU5vgwhkoy4PzB6E0D5yllzmQMpJFYJSRTOwBYPMHq5+xtplrqwNstRbrJdsWKyiLua69balbXlNmrX6EnRHDZdeWy+yr76Ooi5/pG2WXlF+2bT3/Ps0hZ3QNs/pBcGzj9a5QcxPFp8Ac4GVz/Umt8n8opf4AdMHotM7QWruVUheUUqOBTcD9wP9ecqyNwB3Aaq21VkotB35dqwN7MvDjZsYrWlhpUSV7/51D/9HxBEf4c3zHNs4cOcTEh55o38+EVgqiehnLsPu/ud9ZDmX5RvIoO19r21xKzxvrqlJjn7sKXBXGJS93JbjMxfPNZ3ELccUSR8Cjq7x+2EYThVLqXYxv9tFKqSyMkUgvA+8rpR4GTgJ3Amit9yil3gf2Ai7gu1rr6sHrT2CMoArA6MReapa/CSw0O77zMUZNobXOV0r9Cths1nuhumNbtD07Vp3C4/aQNqUbHrebtQveIDwugcHfmuzr0FqWPQDCEo3lang8dROHu9JIJtptXA6rvdb6kjKPuW2utcd89JxueK09xntXl1WrMyeavnzZN8ovdZl9zZ5zTa4qNyqwBabHoWmjnu6+zK6Jl6n/EvBSPeVbgG/clqu1rsBMNPXsewt4q7EYhW9VlDrZ/eVpeo+IIzw2kO3Ll5CXdZIZ//kzbHa7r8NrHywWsAQYiUeINkbuzBZXbeeaLJyVboZP7UZFSQkb3n+HroNT6DVilK9DE0J4gSQKcVWqKlzsXH2K7inRRCUGs/Gf71JZWsqE+x9t330TQogakijEVdmzLpvKMhcjbuxOfnYWmcuXMORbk4npJrP1CtFRSKIQzeZyuslceZKk/hHE9Qjly4VvYvPzY9yse30dmhDCiyRRiGbbvyGHsuIqht/YneM7tnF022ZG3zabwLBwX4cmhPAiSRSiWdxuD9uWnyS+ZygJvUJqhsOm3XiLr0MTQniZJArRLIc2n+VCfgXDb+zOrlXLycs6yXX3PSTDYYXogCRRiCumPZpty04QlRRMXHcH//7gHZIHpdB7xGhfhyaEaAGSKMQVO7I9l4IzZQyf2o1NHy+iouQCE+5/RIbDCtFBSaIQV0RrzdZlxwmPCyQywcn2ZcZw2NjuPX0dmhCihUiiEFfk5J58zp8qYdiUrqx75y1jOOxdMhxWiI5MEoVoMq01W5ceJzjSgSMwh6NbMxh922yCwuUJtUJ0ZJIoRJPlHC4k50gRQycmsW7hm4TFxctwWCE6AUkUosm2LD1BQIgdT+Uu8rJOMv5eGQ4rRGcgiUI0ydnjxZzam8+ga6PZ+NE/SB44hN4jx/g6LCFEK5BEIZpk27ITOAJtlJxfbwyHnSOzwwrRWUiiEI3Kyy7haGYuvYf7sXPlZwy5fpIMhxWiE5FEIRq1bdkJbA4reSeWm7PD3ufrkIQQrajRR6GKzq0ot5xDm8/SbXAZ+9dt5tpvP9Ck4bDny8/zzr53iAmIYVTCKHqG9ZRLVUK0U5IoRIO2fXECLJpzhz8nLC6eYTfNaLTNsuPLeOnrlyiqLEKjAYjyjyI9IZ1R8aNIT0gnOSS5pUMXQniJJApxWSUFlezfmENMl5Oc3HmSW37wkwaHwxZUFPDSppdYfnw5g6MGM2/qPBxWB5vPbGbTmU1sytnE0mNLAUgMTiQ9Pp30hHTS49OJDYxtrY8lhLhCkijEZWWuPInHVcG5o1+QNHAwvdMvPxx21clVvLDxBYqrinkq7Snm9L6b8pWrsUZGcUvKRGb2mYnWmmNFx9h0ZhMZORmsOrmKjw9/DECPsB6kx6czKmEUI+NGEu4vDz8Soq2QRCHqVV5SxZ71pwkJ20nuiRIm3F//cNiiyiJ+k/EbPjv6Gf0j+/P6Da+TsPEIJ565BVd2Tk09vx49CEhJIXxoCremDGX2NXegbVYO5B8g40wGm3I2sfjIYt478B4KRb/IfjWJY3jccILsQa358YUQtUiiEPXaseoUVRV5lOZuYMj1k4jr0esbddZlreP5Dc9TUFHAE0Of4D6dTt6Tz5O9YweO/v1J+OULKKuF8p07Kd+xk5KvvqLok08AUA4H/gMHEpmSwsyhKdyd8lO4Poa9eXvJOJNBRk4Gi/YvYsHeBfhb/bl34L08OPhBQv1CW/ufQohOT2mtfR2DV40YMUJv2bLF12G0a5XlLhb8ZAO6ajGVpcd56E9z64x0ulB1gd9t/h3/Ovwveof35qXezxD21mKKP/sMa3Q0sc88TdjMmSirtc5xtdY4T2dTsXMH5Tt2Ur5jBxV796KrqgCwRkcTkJJiLENTsAzsy66yI3x46EOWHltKqF8oDw95mLv7302ALaBV/02E6OiUUlu11iPq3SeJQlxq67Lj/PuDNThL/sm1336A9Bl31OzbcHoDP9/wc3LLc3m0933c/rWFovkLQGsiH3yQqEcfxRrc9MtEuqqKigMHKd+5gwrzzKPq+HFjp1L49epJ4IgRFE4bw/8WfcL60+uJDYjl8dTHubX3rdgtMteUEN4giUI02ZHt51i9YC8VhQsJCIYHfv83bH5+lDpL+f2W3/PBwQ/oGdKdl4pvwO/ND3Hnnif0ppuI/eEPsCcmeiUGd2Eh5bt2U75zB+U7d1K2KQNdUUHQ2LHkzRjHH2yryTy/g64hXfle2veY3H0yFiX3jgpxNSRRiEaVFlaybtFBjmbm4h94gMLTn3HzD35M31HjyMjJ4Ocbfk52STbfV5O47qMjVO0/gP/QFOJ+9CMC09JaNDZ3YSEF771Pwd//jis3F79evTh/y2j+J2oz+0uPMiByAE8Pe5qxXcbKTX1CNJMkCnFZ2qPZ++9sNvzzMC6Xk4TuJzi1exmx3Xsx/Sf/j1e2v8I/9v+DtKoE/isjDuv6LdgSEoj94Q8JnXZTq/5h1lVVFC9dSt68+VTu24c1MpLcqcP4Y7f9HOAMI+NH8vSwpxkaM7TVYhKio5BEIepVeLaMNX/fz+mD+YRGnaAsfx2lBXl0HzqMhJkTeGnf7zmfe4Kf7etH75UHUX5+RD/2KJEPPIDF399ncWutKduUQf68eZSsXYvy8yP3uoG82u8ke0KLuT75ep5Ke4reEb19FqMQ7Y0kClGH2+0hc8VJMhYfA88xLOprSvKyievZh/4zb2KF3sJ7e//B7XtDuf3LKtSFUsJum0nM009jj21bd1BXHj1G/oL5FP3rE3RFBflDu/H24PNkJFdyc+9beDL1SRKDvdN3IkRHJolC1Dh3opjVC/aTe+IgNuvXlBUeJywugZCJQ1nhv50tpzcx7KjiOxsCCc0uJjA9nbgfPYf/wIG+Dr1BroICCt97j/x33sGde54LSRG8O7SUDYMt3DZoNo8MeYSogChfhylEmyWJQuCsdJOx+Cjbv8jE49yAs+wQjtAQytPj+dx/C4mH8pl4OIDhB93YSyqwd+tK3LPPEjxxYrvqIPZUVVH82efkz5tH5YEDVIQ6+HSoi/UjA5iUegcj40eSFptGhH/jM+AK0ZlIoujkTu3NZ9X8DAqz1+Ku2oPys3FmgJVzpfsZfQhGHbbgX+bEEhRE8Le+ReiUyQRfdx3Kz8/XoTeb1pqyr78mb948Sr9ch8tuYXMfxcEEzbE40H27M7CrkTSGxQ4jKSSpXSVEIbxNEkUnVVHi5Mt3d7Bv/We4q7ajlZvC8BISc3MZdVgTUOFBBQcROnEiIVOmEjRuLBaHw9dhe13lkSPkz1/AhS/X4j57rqb8bKSFw3GaY/GKvK5hhKcMZ3CPUaTFpdEvoh82i8xwIzoPSRSdgNaasuIqCs6UkXe6mHPHcji46SsqijeidSU2Tykjjp0nqsSFDg4k7IbJhE6dQtDYsVja8ZnDlXKdP0/F3r1U7N1L+e49lOzeCWdqJY9wOBqvyEqwY+nfh/i0sQzuM5aU6BQC7YE+jFyIltXuE4VSairwZ8AKvKG1fvlydTt6ovB4NBfyKjh34jzZh06Qc/w4hWdyqCwqwOMsQnsK0Z4LgAeAsHI3Q05mE2izEDJxItHTbiFo9Oh2fVnJ21wFBUby2LOXwp1bKdu9G/uZvJr9uaFwLN5CSY8YrHGx2MMjCIiMJTAqluCoBMJjkogMjSXKP4oQvxC5S1y0S+06USilrMBBYBKQBWwG7tZa762vfntIFC6nm/LSci4UF1FaWERp8QXKiy9QfqGEitJSnOWVVJVV4KyswlVVhauiivKiC7jKy3C7ytCeItAVdY5p0VYcLiuBVU6Cy0sJKS8hQGkiR4+k2613Ezx6NKqBhw6JutxFRVTs20fRjm2cy9yEa98Bgs4UXbZ+uR+U+ENJgKIyyI4z2B9PSCCEhmALD8PiH4DFzw+rIwCrnwOrwx+bfwB2RwA2/wBsjgAcjiDs/gH4+Qdh9w/CERCE1eaHzWLHYrdhtdix2GzYrHYsyoJVWaVfRXhNQ4miPVyETQcOa62PAiilFgEzgHoTRXOd2Lubj174dQM1NKDMdd3/UmdbQ80DQDVoXVN2cXE3I0KFVfvj77bi8FhwYCfA4SY4zJ/QuDjCu3QjLKEbYQndscVEY4uJwRIUJH9ImskaFkbQ6NEEjR5NF54EwFNWhiu/AHdhIVWFeVw4n0NJbjblBblU5edhKSwguKiYkOISrGfL8Tt6Hv+yM1g9V/beTnMpa6COR5mLxVhr87W2KLQyfxpr1kYZUGutLnlt1KtW/dOsatW5WHJJWSv+iGn5eW5Qabdopi1c5fXjtodEkQicqvU6CxhVu4JS6jHgMYCuXbs2603sDn+sloAr+8Naq27tLaWqf8EUqnqfqrVPKSwWhcVmwWK1YLVZsdktWB027HY7dn8/7P52HIF+OAIDCAgOIr5nf6KTemOLiEBZ5NKGL1gCA/ELDISkRAKAsCa00R4PntJSdGUlzopyqipKqawooaqilKqKMpwVpTgrynFWlOGqrMBVWY67sgJXZQWeqgo8Hje4PWiPG22u8XjQHg+43cbaU12u0W638eWkzgI1X1pqv/aY6aC6vCboWhv6G4UX219S3GTNvYrRxq9+tAkJcS1y2PaQKOr7y13nJ+b/t3d3MXJWBRjH/w9UwGIJH5WvbmMLNgISFWIIIDFGNBZsWi9LIGmilxLQaICmCYnXGsALhTSAoV9jzAAABRNJREFUEkG4QNSG8NUAiVdWPgQslEqxBRarrTEokQRofLg4pzosMy8tlJ6z9Pklm533vLvbf97Zd87smemM7bXAWihLT+/lHznx5E9y2S9vei/fGjGRDjqIg+fNg3nzmAPkXTRiNpoNd02ngYUj21PAXxu1REQccGbDRPEIsETSYkmHACuBdY2bIiIOGN0vPdneJelS4H7K02Nvtv1046yIiANG9xMFgO17gHtad0REHIhmw9JTREQ0lIkiIiIGZaKIiIhBmSgiImJQ96/1tLck7QReeB8/Yj7wj32U80HovQ/6b+y9D9K4L/TeB301fsL2x8ft+NBNFO+XpEcnvTBWD3rvg/4be++DNO4LvffB7GiELD1FRMS7yEQRERGDMlG809rWAe+i9z7ov7H3PkjjvtB7H8yOxjxGERERw/IXRUREDMpEERERgzJRVJKWStosaYukq1r3AEhaKOlhSZskPS3p8jp+tKT1kp6rn49q3HmwpD9KurvTviMl3Snp2Xosz+mpUdJ36/W7UdLtkg5r3SfpZkk7JG0cGZvYJGl1PXc2S/paw8Yf1uv5KUm/lnRkq8ZxfSP7vi/Jkua36tsbmSgoN3TAT4ALgNOAiySd1rYKgF3A92yfCpwNfLt2XQU8aHsJ8GDdbulyYNPIdm99Pwbus30K8FlKaxeNkhYAlwGft3065aX0V3bQ93Ng6YyxsU31d3Il8On6PT+t51SLxvXA6bY/A/wZWN2wcVwfkhYCXwVeHBlrdQz3SCaK4ixgi+2/2H4DuANY0bgJ29ttP14vv0q5gVtAabulftktwDfaFIKkKeDrwI0jwz31HQF8EbgJwPYbtl+ho0bKy/1/VNIcYC7lHRyb9tn+HfDPGcOTmlYAd9h+3fZWYAvlnNrvjbYfsL2rbv6e8o6YTRonHEOAa4ErePtbOjc5hnsqE0WxAHhpZHu6jnVD0iLgDGADcJzt7VAmE+DYdmVcR/ml/+/IWE99JwE7gZ/V5bEbJR3eS6Ptl4EfUe5dbgf+ZfuBXvpmmNTU6/nzTeDeermLRknLgZdtPzljVxd9k2SiKDRmrJvnDUv6GPAr4Du2/926ZzdJy4Adth9r3TJgDnAmcL3tM4D/0H4p7H/qOv8KYDFwInC4pEvaVu217s4fSWsoS7e37R4a82X7tVHSXGANcPW43WPGurkNykRRTAMLR7anKH/+NyfpI5RJ4jbbd9Xhv0s6oe4/AdjRKO8LwHJJ2yjLdV+WdGtHfVCu22nbG+r2nZSJo5fGrwBbbe+0/SZwF3BuR32jJjV1df5IWgUsAy72//+jWA+NJ1PuEDxZz5kp4HFJx3fSN1EmiuIRYImkxZIOoTyotK5xE5JEWVvfZPuakV3rgFX18irgt/u7DcD2attTthdRjtlDti/ppQ/A9t+AlyR9qg6dDzxDP40vAmdLmluv7/Mpj0X10jdqUtM6YKWkQyUtBpYAf2jQh6SlwJXActuvjexq3mj7T7aPtb2onjPTwJn1d7R53yDb+Sh3Oi6kPEvieWBN657adB7lz8+ngCfqx4XAMZRnnTxXPx/dQeuXgLvr5a76gM8Bj9bj+BvgqJ4agR8AzwIbgV8Ah7buA26nPGbyJuUG7VtDTZQlleeBzcAFDRu3UNb6d58vN7RqHNc3Y/82YH7LY7inH3kJj4iIGJSlp4iIGJSJIiIiBmWiiIiIQZkoIiJiUCaKiIgYlIkiIiIGZaKIiIhBbwFQJsOv0ksqwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "S1, S2, I1, I2, R1, R2 = output_[\"y\"]\n",
    "plt.plot(output_[\"t\"], S1, label=\"S1\")\n",
    "plt.plot(output_[\"t\"], S2, label=\"S2\")\n",
    "plt.plot(output_[\"t\"], I1, label=\"I1\")\n",
    "plt.plot(output_[\"t\"], I2, label=\"I2\")\n",
    "plt.plot(output_[\"t\"], R1, label=\"R1\")\n",
    "plt.plot(output_[\"t\"], R2, label=\"R2\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Abstract away the boilerplate in Base class and define a new model as subclass "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BaseModel:\n",
    "    \n",
    "    state_names = None\n",
    "    parameter_names = None\n",
    "    parameters_stratified_names = None \n",
    "    stratification = None\n",
    "   \n",
    "    def __init__(self, states, parameters):\n",
    "        \"\"\"\"\"\"\n",
    "        self.parameters = parameters\n",
    "        self.initial_states = states\n",
    "        \n",
    "        if self.stratification:\n",
    "            self.stratification_size = parameters[self.stratification].shape[0]\n",
    "        else:\n",
    "            self.stratification_size = 1\n",
    "\n",
    "    def integrate(self):\n",
    "        \"\"\"to overwrite in subclasses\"\"\"\n",
    "        raise NotImplementedError        \n",
    "        \n",
    "    def create_fun(self):\n",
    "        \"\"\"Convert integrate statement to scipy-compatible function\"\"\"\n",
    "\n",
    "        def func(t, y, *pars):\n",
    "            \"\"\"As used by scipy -> flattend in, flattend out\"\"\"\n",
    "            \n",
    "            # for the moment assume sequence of parameters, vars,... is correct\n",
    "            y_reshaped = y.reshape((len(self.state_names), self.stratification_size)) \n",
    "            dstates = self.integrate(t, *y_reshaped, *pars) \n",
    "            return np.array(dstates).flatten()\n",
    "\n",
    "        return func\n",
    "    \n",
    "    def sim(self, time):\n",
    "        \"\"\"\"\"\"        \n",
    "        fun = self.create_fun()\n",
    "        output = solve_ivp(fun, time, \n",
    "                           list(itertools.chain(*self.initial_states.values())), \n",
    "                           args=list(self.parameters.values()))\n",
    "        return output[\"t\"], self.array_to_variables(output[\"y\"]) # map back to variable names\n",
    "            \n",
    "    def array_to_variables(self, y):\n",
    "        \"\"\"Convert array (used by scipy) to dictionary (used by model API)\"\"\"\n",
    "        return dict(zip(self.state_names, y.reshape(len(self.state_names), \n",
    "                                                    self.stratification_size, -1)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# User (model developer) specifies...\n",
    "\n",
    "# ... a model subclass\n",
    "class SIR_S(BaseModel):\n",
    "\n",
    "    # ...state variables and parameters\n",
    "    state_names = ['S', 'I', 'R']\n",
    "    parameter_names = ['gamma']\n",
    "    parameters_stratified_names = ['beta']  \n",
    "    stratification = 'nc'\n",
    "    \n",
    "    # ..transitions/equations\n",
    "    @staticmethod\n",
    "    def integrate(t, S, I, R, gamma, beta, nc):  \n",
    "        \"\"\"Basic SIR model\"\"\"\n",
    "\n",
    "        # Model equations\n",
    "        N = S + I + R\n",
    "        dS = nc @ (-beta*S*I/N)\n",
    "        dI = nc @ (beta*S*I/N) - gamma*I\n",
    "        dR = gamma*I\n",
    "        \n",
    "        return dS, dI, dR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model user defines...\n",
    "\n",
    "# ... parameters and initial conditions\n",
    "nc = np.array([[0.9, 0.2], [0.8, 0.1]])  \n",
    "parameters = {\"gamma\": 0.3, \"beta\": np.array([0.5, 0.4]), \"nc\": nc}  \n",
    "initial_states = {\"S\": [440000, 350000], \"I\": [20, 10], \"R\": [0, 0]}\n",
    "\n",
    "# -> user initiates the model\n",
    "sir_model = SIR_S(initial_states, parameters)\n",
    "\n",
    "# -> user runs a simulation for a defined time period\n",
    "time = [0, 150]\n",
    "t, output = sir_model.sim(time)\n",
    "\n",
    "# -> user can do fit, mc,... using the model class instance `sir_model`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7feebc02e950>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD4CAYAAADy46FuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3wc1bnw8d/ZolXv1ZLce5ElF7kBNjEuYIMxzSYU0wPkBkhyuaS9CSGQkNybxg0J+FJcQjAlELDBNq7YxMZyk3vvsmRbVrX6lvP+MSNZMrIkyyutyvP9fIaZPXPO7LNG0rMz58wZpbVGCCGEuByLrwMQQgjRtkmiEEII0SBJFEIIIRokiUIIIUSDJFEIIYRokM3XAXhbdHS07t69u6/DEEKIdmXr1q3ntdYx9e3rcImie/fubNmyxddhCCFEu6KUOnG5fXLpSQghRIMkUQghhGiQJAohhBANkkQhhBCiQZIohBBCNEgShRBCiAZJohBCCNGgDncfRXOVVbl4be0R44VSKGo2USgsytxWCotSKAUWBRZl1LRaFDarBXv12qqwWSzYrOqSbQuBflbCAuyEBdgJdthQStUflBBCtAGSKEzlVW7+d81hWvvxHDaLIjTATniA3VgHGgkk3EwkYYF+RATaSU0Op0d0kCQVIUSrk0Rhigp2cOw30+qUaa3RGrS57dGgMco8+uLa4wG31rg8HlxujcutcZrbTrcHl0fjcntwuo06pZVuisqrKCp3UljmNNblTorLneSXVnE0t5SicifFFc46iSsxPIBrekdzbd9oxvWKJiLIr3X/kYQQnZIkigYo8xKT+arV39/t0ZRUuDh3oYJNx/JZfyiXz3fn8N6WUygFg7uEcW2faK7pE83wbhE4bNZWj1EI0fGpjvYo1BEjRuiOPNeTy+1h5+ki1h88z1eHc9l+shCXRxNgtzKqZ6RxxtEnhr5xwXKZSgjRZEqprVrrEfXuk0TRvl2ocPL10Xy+OpTL+kPnOXq+FIDYEAc3DUngJzcNwM8mg9uEEA1rKFHIpad2LsTfzqSBcUwaGAfA6cJyvjqUy5cHc5m34Tg5ReX85dvDsFslWQghmkf+enQwieEBzBrZlb/eM5yfTx/I8j1n+cH7O3B7OtaZoxCi9cgZRQf20DU9qHR5+O2y/ThsFn53ewoWi/RbCCGujCSKDu6JCb2ocLr586pD+Nst/GrGYOnkFkJcEUkUncAzN/ShwuXm9S+P4rBZ+dm0AZIshBBNJomiE1BK8aOp/al0enjzq2P42y08O6W/r8MSQrQTkig6CaUUv7h5IJUuN6+uOYK/zcr3JvbxdVhCiHZAEkUnopTipVuHUOn08PsVB/G3W3n0up6+DksI0cZJouhkLBbF7+5IodLl4aXP9+GwW7h/THdfhyWEaMMkUXRCNquFP81OpdLl4eef7MFhszBrZFdfhyWEaKOafMOdUsqqlNqulFpivo5USq1QSh0y1xG16v5YKXVYKXVAKTWlVvlwpdQuc98ryhx6o5RyKKXeM8s3KaW612ozx3yPQ0qpOd740ALsVguv3pPG+L4x/OijXfxr+2lfhySEaKOu5M7sp4F9tV7/CFilte4DrDJfo5QaCMwGBgFTgb8qpaqnNf0b8BjQx1ymmuUPAwVa697AH4HfmseKBH4BjALSgV/UTkji6jhsVl6/bzije0Txg/cz+XxXjq9DEkK0QU1KFEqpJGAa8Eat4hnAfHN7PnBrrfJFWutKrfUx4DCQrpRKAEK11hu1MRPhgkvaVB/rQ2CiebYxBVihtc7XWhcAK7iYXIQX+NutvDFnBMO6RvDUu9tZufesr0MSQrQxTT2j+BPwX4CnVlmc1joHwFzHmuWJwKla9bLMskRz+9LyOm201i6gCIhq4Fh1KKUeU0ptUUptyc3NbeJHEtWCHDbeenAkA7uE8uQ721h3UP4NhRAXNZoolFLTgXNa661NPGZ9t/zqBsqb2+ZigdZztdYjtNYjYmJimhimqC3U386Ch9LpFRvMYwu3sO1kga9DEkK0EU05oxgH3KKUOg4sAr6llPo7cNa8nIS5PmfWzwKSa7VPArLN8qR6yuu0UUrZgDAgv4FjiRYQHujH3x9OJ8Tfzl/XHPZ1OEKINqLRRKG1/rHWOklr3R2jk3q11vpe4FOgehTSHOATc/tTYLY5kqkHRqd1hnl56oJSarTZ/3D/JW2qj3WH+R4aWA5MVkpFmJ3Yk80y0UKigh3MTEtk7YFc8koqfR2OEKINuJrnUbwMTFJKHQImma/RWu8B3gf2AsuA72qt3WabJzA6xA8DR4ClZvmbQJRS6jDwA8wRVFrrfOBXwGZzecEsEy1oZloiLo9myU4ZBSWEkEehisuY+qd1+Nut/Ou743wdihCiFTT0KFR5wp2o18y0RDJPFXI0t8TXoQghfEwShajXjNRElELu2BZCSKIQ9YsP82dsryg+zjxNR7s8KYS4MpIoxGXNTEviVH45W0/IPRVCdGaSKMRlTR0cj7/dwkdy+UmITk0ShbisYIeNKYPi+WxnDpUud+MNhBAdkiQK0aBb0xIpKneyZr/M/yREZyWJQjTo2t7RRAc7+Hh7VuOVhRAdkiQK0SCb1cItQ7uwev85CsuqfB2OEMIHJFGIRs1MS8Tp1nwmDzYSolOSRCEaNTgxlN6xwXy8TUY/CdEZSaIQjVJKMTMtkS0nCjiZV+brcIQQrUwShWiSGaldAPhXppxVCNHZSKIQTZIUEcioHpF8vF2m9BCis5FEIZrstmGJHDtfSuapQl+HIoRoRZIoqnk8UJZfdykvgPJCY6kogsoLUFkCVWXgLAdXJbiqwO002ndwUwcn4GezyIyyQnQyNl8H0GaU58N/97q6YygLWGy1Fms9r+1g9wdHGDhCwD8UHKGXbIea2yHGdkA4hHQBi2/zeliAnUkD4li8M4efTR+I3SrfM4ToDCRRVPMLght/BzXX33U92xq0x9jWnm++9rhqLe6L225nrddOcFZAZTEUZcE580ylohh0A/MpBURAUjokp0PyKEgcZsTcymamJfLZrhy+PJDLDQPjWv39hRCtTxJFNXsAjPqO795fa3CWXUwaleZSUQxleZC9HU5lwKHlRn1lhfgh0HX0xeQRltTiYY7vF0NEoJ2PM09LohCik5BE0VYoZZwh+AVBSPzl65XlQ9YWOLXJWLYtgE2vGftCEy8mjeR0iE8Bq92rYdqtFm4e2oVFm09RXOEk1N+7xxdCtD2SKNqbwEjoO9lYwLisdXa3cbZxapOx3vOxsS84Dmb/A5LqfV56s81MS2TBxhMs3ZXDrJFdvXpsIUTbI72R7Z3VDl3SjMtmd7wF398N399rbNv8Yd502LfYq2+ZmhxOj+ggPpIpPYToFCRRdERhiTD4dnhkFcQNhPfug41/9drhlVLcmprIpmP5nC4s99pxhRBtkySKjiw4BuYsgf7TYPmPYelzxugrL5iZlggg91QI0QlIoujo/ALhrgUw+rtGp/d790JV6VUftmtUICO6RciUHkJ0ApIoOgOLFab+2rhP5OAyo9+i5NxVH/bWtEQOnythT3axF4IUQrRVkig6k1HfgVnvwLl98MZEyD1wVYebnpKAn9UindpCdHCSKDqb/jfBg58Zd4e/OQmOrW/2ocID/bi+fwyf7sjG5e74c10J0VlJouiMEofDIyshOB4WzoQd7zX7UDPTEjlfUslXh897MUAhRFvSKW64czqdZGVlUVFR4etQms3f35+kpCTsdi/dCR3RDR5ebgyd/fgxKDwB1z1r3CF+Ba7vH0tYgJ2Pt59mQr9Y78QmhGhTOkWiyMrKIiQkhO7du6Ou8A9hW6C1Ji8vj6ysLHr06OG9AwdEwL0fwaffgzUvGcli+p+uaNoPh83KtJQEPtqWRUmli2BHp/iREqJT6RSXnioqKoiKimqXSQKMG9yioqJa5ozI5gczX4Pxz8H2v8M7dxjP3rgCt6UlUuH0sHz3Ge/HJ4TwuU6RKIB2mySqtWj8SsH1P4EZf4XjX8FbU6E4u8nNh3eLIDkygI/l5jshOqROkyjagpdeeolBgwaRkpJCamoqmzZt4i9/+Qu9e/dGKcX58z7uEE67B+79JxQchxU/b3IzpRQzUxP595HznClqv/1AQoj6SaJoJRs3bmTJkiVs27aNnTt3snLlSpKTkxk3bhwrV66kW7duvg7R0HMCDJsDe/4FF5p+KenWtES0hk93yFmFEB2N9Dy2kpycHKKjo3E4HABER0cD0KVLF1+GVb/0R43pPra8ZVySaoKeMcEMTQ7no22neey6q3ykbCekPZrKcheVZS4qy5zm2oXb5cHt8uBx64vbLo3bba5dHjwuD263NtZmXY9Hmw9l1Mbj3LVG19k23lPX2TZjMTeqH+oIF/ddLK9V59LPUqdyPZvtYMqXdhBivaKTgrnpiRSvH7fRRKGU8gfWAQ6z/oda618opSKB94DuwHHgLq11gdnmx8DDgBt4Smu93CwfDswDAoDPgae11lop5QAWAMOBPGCW1vq42WYO8DMznBe11vOv5gP/cvEe9np5yomBXUL5xc2DGqwzefJkXnjhBfr27csNN9zArFmzGD9+vFfj8JqoXtB3ipEorv0h2BxNanZbWiK/+HQP+3KKGZAQ2sJBtn1VFS7ys0vJzyml/EJVzR//2omgZrvcVeePalNYbAqr1WKsbZaabYvVYjxeXSksFoVSoKrXSl3cthn1lLEDZV5fUFAzTLq6a8yoY+67WMnYX6v77OL++vvUaoob6HJTDe1sTW0kjCsRFhPQIsdtyhlFJfAtrXWJUsoOfKWUWgrcBqzSWr+slPoR8CPgOaXUQGA2MAjoAqxUSvXVWruBvwGPAV9jJIqpwFKMpFKgte6tlJoN/BaYZSajXwAjMH6NtiqlPq1OSO1JcHAwW7duZf369axZs4ZZs2bx8ssv88ADD/g6tPqNehwW3gq7P4LUu5vUZHpKAr9aspePt5/uVInC49EU55aTd7qE86dLyMsqIe90CcXn6/bXWGwKR6Ad/0AbjkAbgaF+RMQH4gi04zDLqrf9g2z4Bdix2WslApsFi1XVrNv7AA3RfjSaKLRxHllivrSbiwZmABPM8vnAWuA5s3yR1roSOKaUOgykK6WOA6Fa640ASqkFwK0YiWIG8Lx5rA+Bvyjjt2AKsEJrnW+2WYGRXN5t7gdu7Jt/S7JarUyYMIEJEyYwZMgQ5s+f33YTRc8JENMfNv0Nhs5u0o14UcEOxveN4ZPM0zw3tT9WS8f7Q1ZR4ryYEMykkJ9TiqvKmMJEKQiLDSSmaygDxiYQlRhMZJdgAsP8sNkt8sddtEtN6qNQSlmBrUBv4FWt9SalVJzWOgdAa52jlKq+LTcR44yhWpZZ5jS3Ly2vbnPKPJZLKVUERNUur6dN7fgewzhToWvXtvlozgMHDmCxWOjTpw8AmZmZbacDuz5KGZMILvm+8YjVrqOb1GzmsERW7T/HxiN5XNMnuoWDbB1ut4d9/85h+4qTFOdefFCTf7CdqMRgBl2TSFRSkJEUEoKw+Vl9GK0Q3tekRGFeNkpVSoUDHyulBjdQvb6vTLqB8ua2qR3fXGAuwIgRI9pkN1RJSQnf+973KCwsxGaz0bt3b+bOncsrr7zC7373O86cOUNKSgo33XQTb7zxhq/DNaTMgpXPGx3bTUwUNwyII8Rh46PtWe0+UXg8mkMZZ8hYcozi8xXE9wxl8LUXk0JgqJ+cIYhO4YpGPWmtC5VSazEu/5xVSiWYZxMJQPUDDrKA5FrNkoBsszypnvLabbKUUjYgDMg3yydc0mbtlcTcVgwfPpwNGzZ8o/ypp57iqaee8kFETeAXZAyV3fgqFGVBWFKjTfztVm4cEs9nO3Mou9VFoF/7G1intebo9lw2LT5GQU4p0cnBTPtuCt0Gt9+7+4W4Go3eR6GUijHPJFBKBQA3APuBT4E5ZrU5wCfm9qfAbKWUQynVA+gDZJiXqS4opUab/Q/3X9Km+lh3AKvNvpHlwGSlVIRSKgKYbJaJ1jLyEUDD5jeb3GRmWhKlVW5W7D3bcnG1AK01J3bn8cFvtrBs7m7QmimPDuauH4+k+5BoSRKi02rK170EYL7ZT2EB3tdaL1FKbQTeV0o9DJwE7gTQWu9RSr0P7AVcwHfNS1cAT3BxeOxScwF4E1hodnznY4yaQmudr5T6FbDZrPdCdce2aCUR3aDfTbB1Hoz/L7A3PvxuVI9IuoT589G208xI/UaXUpuUfaiArz85Ss7hIkKi/Jk4ZwB9R8Vj6YAd8kJcqaaMetoJpNVTngdMvEybl4CX6infAnyjf0NrXYGZaOrZ9xbwVmNxihY0+gnYvwR2fQDD7m+0usWimJGWyOtfHiH3QiUxIU27D8MXzh4vZtOnRzm1N5+gMD/G392XAeO6YLXJpAVCVJPfBtG4buMgbjB8/VqTb1m9LS0Rj4Z/tdGJAvNOl/D533by4ctbyD1xgbG39+beX41h8PgkSRJCXKL99TSK1qeUcQPep/9hzC7b49pGm/SJCyGtaziLNp/kkWt7tJnr+4XnyshYfIxDW87i57CSfnMPhk5Mxs9ffhWEuBz57RBNM+QOY0bZTa81KVEAzB6ZzHP/3MWWEwWM7B7ZwgE2rLSokozFx9i3IQerTTFscjfSJnfFP8hLTwwUogOTc+xWVN804/fccw/9+vVj8ODBPPTQQzidTl+HWT97AIx4EA58bkxD3gTTU7oQ5GdlUcapxiu3ILfLw+JXdrD/6xwGj0/k3l+NYczMXpIkhGgiSRSt5HLTjN9zzz3s37+fXbt2UV5e3nZutqvPiIcBBRn/16TqQQ4bt6Qm8tmubIrKfZcAt3x+nLzTJUx9bAjXzepLUFjb7VwXoi2SRNFK6ptmvEuXLtx0003GjJ5KkZ6eTlZWViNH8qGwRBg4A7YthMqSxusDd6cnU+H08Gmmbzq1c09eYOuyE/QfHU+PlPZ9p7gQvtL5+iiW/gjO7PLuMeOHwI0vN1ilsWnGnU4nCxcu5M9//rN3Y/O2UY/Dno9g5yLzZryGDUkMY0BCKIs2n+K+Md1bPr5a3E4PK+ftJTDEzjV39WnV9xaiI5EzilZSPc343LlziYmJYdasWcybN69m/5NPPsl1113Htdc2raPYZ5LToUsabHq9SUNllVLcnZ7MnuxidmUVtUKAF23+7Bj52aVMuLc/jkDpjxCiuTrfGUUj3/xb0uWmGf/lL39Jbm4ur7/+us9ia7LqobIffweOrIbe9d5zWceM1ERe+mwf724+yZCkIa0QpHEj3bYvTjJgbALdh8glJyGuhpxRtJIDBw5w6NChmtfV04y/8cYbLF++nHfffReLpZ387xg0E4JijbOKJggLsDNtSAKfZmZTVuVq4eDA5XSzav4+gsL8GHenXHIS4mq1k79M7V9JSQlz5sxh4MCBpKSksHfvXp5//nkef/xxzp49y5gxY0hNTeWFF17wdaiNszlgxENwaDnkHWlSk9npXSmpdLFkZ04LBweblxizvl5/b38cAZ3vpFkIb5PfolZyuWnGXa6W/4bdIkY8BOt/Dxlz4cbfNlp9ZPcIesYEsSjjJHeNSG60fnOdOVrE9i9OMvCaLnQdFNVi7yNEZyKJQjRPSBwMvg22vwPX/xT8G35GtlKK2SOT+fXn+zl49gJ940K8HpKrys3qBfsIinAw7vbeXj9+S/F43FSVl1NVVkZleZm5LqWqvBy304nH7cbjduF2ufC4XLjdbjwul1HmduNxOXG7LtbRbjdaa2PxeMxtD9qja7aps8/YBm2MT9AaXf18MG1Mv179vLDq/cYuXXdAwyVjG3Ttglr1mjhdWD3a5DPJ2pTo5G5MffL7Xj+uJArRfKO+Azvfg8x/wOjHG61++7Ak/nv5ARZlnOLnNw/0ejibFh+j4EwZtzydip8PLzlprSktyCc/+zQFOVkUnMmhouQCVeVlVJWXU1lWSlVZGVXlZVSWleGsrGjeGymF1WbDYrVhtVqx2GzGYrFisVhQFgXKUnOfjrJUbxv7qrexKBTKqG8c2JibSwEKLMoCStV6bLqqeYb6pXN4fWNOr1qv1WXKr+wjt405w9oq/2DvfwEDSRTiaiQOh6R0yHgd0h+DRjrjo4IdTB4Yz0fbs/ivqf3wt3vv2dI5R4rIXHmSQdclkjygdeaVqiovoyAnm/yc0xRkZxnb5tpZcfHZ2ja7H/6hoTgCAvELCMA/KJjQ6Fj8AgJxBAZeXAcGmnUullltdiw2q7G2Wo2kYLNhsVmxWOTZ3KJ1SKIQV2fUd+CfD8PhFdB3SqPVZ6cn89muHJbvOeO1hxo5q9ysmr+XkEh/xt7WyyvHvFRFSQn7N6wj98RRCnKyKcjOoqSg1jO0lCI0OpbILokk9h9IREIikQlJRHRJJCQyCtVeRrQJUQ9JFOLqDJwBX/zMmFW2CYliXK9okiICeG/zKa8lik2fHKXoXDkzvp/m9enC87NPs33Zp+xZuwpnZQX+QcFEdEmkW0oaEQmJRHRJJDIhkfD4Ltj8/Lz63kK0FZIoxNWx2mHkw7D6Rcg9ADH9GqxusShmjUjm9ysOciKvlG5RQVf19tmHCtmx+hRDxieS1C/iqo5VTWvNyV072Lb0E45u24zVZqP/uAkMu+kWYrv39Mp7CNGeyPlwKwkODq7Znjp1KuHh4UyfPt2HEXnR8AfB6mjyDXh3jkjGomDR5qubftxZ6WbVgn2ERvkzeubVX3JyVlWya/UXLHj2P/jwpZ9x5sghxtzxbR599W2mPvmMJAnRackZhQ88++yzlJWVtY8pO5oiKBqG3Ak73oWJ/w8CGv5mHx/mz/X9YvlwaxY/mNQXu7V531c2/usIxbnl3PqDq7vkVFKQz44vPmPHiqWUXygmplsPpj75ffqNvQ6bXeaIEkIShQ9MnDiRtWvX+joM7xr1Hcj8O2z/O4z9XqPVZ6d3ZdWCLazef44pg+Kv+O1OHyxg15osUq5PIrFv8y45nT16mG2ff8L+DevxeNz0Gj6K4TfdQtLAITIMU4haOl2i+G3Gb9mfv9+rx+wf2Z/n0p/z6jHbnYQU6DbOuFN79JPQyNDN6/vFEBfqYFHGyStOFFUVLlYv2EdYTACjb72yS04ej5sjmzex9fNPOL1/D3b/AFIn30Ta1JsJj0+4omMJ0Vl0ukQhWtCo78D798OBpTCg4f4Xm9XCncOT+evaw2QXltMlPKDJb7Px4yMU51Uw84fDsDuadi9BVXkZO1ctZ/uyJRTnniU0Jo4J9z/C4Osn4Qi8ug51ITq6TpcoOv03/5bUbxqEJRtDZRtJFAB3jUjmL2sO88GWLJ6+oWmzvGbtz2f3l6cZOjGZLr3Dm9TGVVXFB7/6KWeOHCJpwGAm3P8wvUaMkhvWhGgiGfUkvMdqM556d3w9nNndaPWuUYFc0zua97ecwu1pfB4f45LTfsLjAhk1o2kjkLTWrHzjr5w5coibv/8jZj3/Mn3Sx0qSEOIKSKLwgWuvvZY777yTVatWkZSUxPLly30dkvcMux9sAca0Hk0wOz2Z04XlfHX4fKN1N/zzMCUFFUycMwC7X9P+0O/44nP2fLmS0bffTd/R1zSpjRCirk536clXSkpKarbXr1/vw0haWGAkDJ0FOxbBxOchqOGpvicNjCMi0M6ijJOM7xtz2Xqn9uazZ302qZO6Et8zrEmhZO3fw5r5c+k5bCRj77j7Sj6FEKIWOaMQ3jfqcXBVwLb5jVZ12KzcPiyJFXvPknuhst46VeUuVi/cR0R8IKNu7tGkEC7kn2fxH35DWGwcN/7HD2WuJSGugvz2CO+LHQA9xsPmN8DtbLT67PRkXB7NP7dl1bv/3x8eorSwkm/NGYCtCZecXE4ni3//G5yVlcz4z5/hHxTcaBshxOVJohAtY/QTUHwa9n3aaNXesSGM6BbBe5tPmQ/JuejU/nz2/juHtMldie/RtEtOq99+jZzDB7jxye8TldS1WeELIS6SRCFaRp/JENUb/v3nJj3SbHZ6V46dL2XTsfw65bvWZBEY6sfI6U275LRjxVJ2rVrOqJl30WfU2GaFLoSoSxKFaBkWK4x7GnJ2wNE1jVafNiSBEH8bizJO1pRVlDo5sTuPPulx2JrwkKPTB/ax+u3X6Z46nLF33XNV4QshLpJEIVpOyiwISYCv/tho1QA/K7emJvL57jMUlRn9Goe3nsPj1vRLb3yKj5L8PBb/4deEREcz7XvPyn0SQniRJIpWUj3NeGZmJmPGjGHQoEGkpKTw3nvv+TiyFmRzGPM+HVsHp7c2Wn3WyGSqXB4+3m50ah/MOENEfCDRyQ13RrtdTj7942+oKi83Oq+DpfNaCG+SRNHKAgMDWbBgAXv27GHZsmU888wzFBYW+jqsljP8AfAPg6/+1GjVwYlhDEkMY9HmUxSfLyfncBF90+Mbncl1zby55Bzcz5QnniGma3fvxC2EqCGJopX17duXPn2MeY26dOlCbGwsubm5Po6qBfmHwshHYd9iOH+o0eqz05PZf+YC61aeAKBvelyD9XeuWs6OFUsZOeMO+o2RO6+FaAmN3pmtlEoGFgDxgAeYq7X+s1IqEngP6A4cB+7SWheYbX4MPAy4gae01svN8uHAPCAA+Bx4WmutlVIO8z2GA3nALK31cbPNHOBnZjgvaq0bv4urAWd+/Wsq93l3mnHHgP7E/+QnV9wuIyODqqoqevW6+qeztWmjHoeNfzFGQM34S4NVbxnahRcX7+Pw5rN07RVGaPTlZ5XNOXSA1W/9jW4paVwz+z5vRy2EMDXljMIF/FBrPQAYDXxXKTUQ+BGwSmvdB1hlvsbcNxsYBEwF/qqUqu5Z/BvwGNDHXKaa5Q8DBVrr3sAfgd+ax4oEfgGMAtKBXyilvPNgZB/Lycnhvvvu4+2338bS0e8aDo6BtPuMaT2KsxusGuJv57aeMdhL3XQddvkpPUoLC/j0D78mODKKaU//l3ReC9GCGj2j0FrnADnm9gWl1D4gEZgBTDCrzQfWAs+Z5Yu01pXAMaXUYSBdKXUcCNVabwRQSi0AbgWWmm2eN4/1IfAXZVyYngKs0Frnm21WYCSXd5v7gZvzzd/biouLmTZtGt++oisAACAASURBVC+++CKjR4/2dTitY+x/wJa34Ou/wuQXG6w6XDnIQnPQ7mZEPfvdLieL//gbKkpKuPtX/01AcEjLxCyEAK6wj0Ip1R1IAzYBcWYSqU4msWa1ROBUrWZZZlmiuX1peZ02WmsXUARENXCsS+N6TCm1RSm1pa1f76+qqmLmzJncf//93Hnnnb4Op/VEdIfBt8GWt6G84LLVPB5N8f5CzgUp3ttV/9nH2gVvcnr/XiY//hSx3Zs23bgQovmanCiUUsHAP4FntNbFDVWtp0w3UN7cNhcLtJ6rtR6htR4RE3P5yxVtwfvvv8+6deuYN28eqamppKamkpmZ6euwWse4Z6CqxJgD6jKyDxZQWlRF12Ex7DhVyL6cuj9qu9euJHP5EoZPn8mAceNbOmIhBE1MFEopO0aSeEdr/ZFZfFYplWDuTwDOmeVZQHKt5klAtlmeVE95nTZKKRsQBuQ3cKx2p3qa8XvvvRen00lmZmbNkpqa6uPoWkn8YGNqj69fg6qyeqscyDiL3d/KzGm98bNa6typfebIIVa+8SpdB6dw3bcfaKWghRCNJgqzr+BNYJ/W+g+1dn0KzDG35wCf1CqfrZRyKKV6YHRaZ5iXpy4opUabx7z/kjbVx7oDWK2N2eGWA5OVUhFmJ/Zks0y0V+OegbLzkPnON3a5nG6ObjtHr9QYYsIDmDo4no+3n6bC6aasqJBPf/9rgsIjmPb0c1is0nktRGtpyoOLxgH3AbuUUtXXSH4CvAy8r5R6GDgJ3Amgtd6jlHof2IsxYuq7Wmu32e4JLg6PXWouYCSihWbHdz7GqCm01vlKqV8Bm816L1R3bIt2qttYSEqHDa/A8AeNx6eaju/Mo6rCTd9RxpQds0cm8+mObD7fkYXzs79RXlzE7Bd+R2Bo02aRFUJ4R1NGPX1F/X0FABMv0+Yl4KV6yrcAg+spr8BMNPXsewt4q7E4RTuhFFzzfVh0N+z5GFIu/m8/mHGGwDA/EvsZI6BH94yiW1QgG96dR0LWbm787g+I69nbV5EL0Wl18AH8ok3qOxVi+huTBZpTkNfMFDsyDovF+F5isShmxRSRkLUZBl3LgGuv92XUQnRakihE67NYjL6Kc3vg0Arg8jPFhh5chys4mldLB/Cbpfu/8WAjIUTLk0QhfGPIHRCaVDMFeX0zxeYcPsDZIweZdMft3Du2J3PXHeXnn+zB45FkIURrkkTRSqxWK6mpqQwePJibb765zoyxU6dOJTw8nOnTp/swwlZmtRt3a5/cQPGujfXOFJu5bAl+AQEMnjCRX94yiMeu68nCr0/w3D934pZkIUSrkUTRSgICAsjMzGT37t1ERkby6quv1ux79tlnWbhwoQ+j85Fh90NABIc+/xKoO1NsWVEhBzauZ+B1E/ELCEQpxY9v7M/TE/vwwdYsnnkvE6fb46vIhehUJFH4wJgxYzh9+nTN64kTJxIS0gnnK/ILQqc/zsFT8SR0tdeZKXbnymW4XS7Spl48y1JK8f1JfXluan8W78jmu+9so9Llru/IQggvasp9FB3K+vcPcv5UiVePGZ0czLV39W1SXbfbzapVq3j44Ye9GkN7lZd0H/muA4wP+hq4FgC3y8WOFZ/TLSWNyC5J32jzxIReBNgtPL94L48t2Mpr9w4nwE9uwBOipcgZRSspLy8nNTWVqKgo8vPzmTRpkq9DahMO7CzHojz0yn8VCo35Hw9v3khJQT5pU2++bLsHxvXg5duGsO5QLg/Oy6C00tVaIQvR6XS6M4qmfvP3tuo+iqKiIqZPn86rr77KU0895ZNY2gqPR3No81m69g8hoOgCbHwVbnyZ7csWExYXT4+04Q22n53eFX+7lR9+sIP73tzE2w+mExZgb6Xoheg85IyilYWFhfHKK6/wP//zPzidTl+H41PZBwsoLayk77juMOQu2Dafc/u2cXr/XlInT2vSw4huTUvk1W+nset0Efe88TX5pVUtH7gQnYwkCh9IS0tj6NChLFq0CIBrr72WO++8k1WrVpGUlMTy5Z1j3sOD5kyx3VOiYdzT4Cxj+6K/YHM4GDyh6Zfmpg5OYO59Izh4toS7537NuQsVLRi1EJ1Pp7v05CvV04xXW7x4cc32+vXrWzscn3M53RwxZ4q1+1khtj/lPW5k/9IzDJwwCf/g4MYPUsv1/WN5+4GRPDJ/C7Nf/5p3Hh1FQtjln7cthGg6OaMQPlEzU2ytKTt2qVG4tIXUpOZdPhrXO5qFD6dz7kIld762kVP59T/zQghxZSRRCJ84mHGGwFA/EvsbM8V6PG52ZOwiOVITc3ABuJqXLEZ0j+SdR0ZxocLFna9t5Eiud4dCC9EZSaIQra6+mWKPbM2gOPccaVOmQ3EW7P6w2ccfmhzOosdG4/J4mPX6RvafaejJvUKIxkiiEK3uyDZzpthRFy87ZS5bQkhUDL2mPwpxg+GrP4Gn+VN0DEgIZdFjY7BaFLPnfs2urCJvhC5EpySJQrS6gxln68wUm5d1kpO7dzB00o1YbDZjCvLzB+Dg0kaO1LDescG8/50xBPnZ+Pb/fc36Q7kyTbkQzSCJQrSq4rxysg8V1pkpdvuyJVjtdoZMnGJUGjQTwrvWebBRc3WLCuKDx8cQFezHfW9mMOrXq3j2gx18tjOHovLOfR+LEE0lw2NbidVqZciQIbhcLnr06MHChQsJDw8nMzOTJ554guLiYqxWKz/96U+ZNWuWr8NtMYc2nwUuzhRbWVbK3nWr6T92/MVnYVttMPYp+Pw/4cQG6D7uqt6zS3gAn37vGr7Yc5a1B86xfM8ZPtiahdWiGNY1nAn9YhnfN4ZBXULrTHMuhDBIomgl1VN4AMyZM4dXX32Vn/70pwQGBrJgwQL69OlDdnY2w4cPZ8qUKYSHh/s4Yu/TWnMw4yzxPcNqZordvWYlzsqKOrPEApB2L6x92TiruMpEARDqb+eO4UncMTwJl9tD5qlC1h7IZe3Bc/z38gP89/IDxIY4GN83hgn9YrmmT7RMByKESRKFD4wZM4adO3cC0LfvxbmnunTpQmxsLLm5uR0yUeSdLiE/u5TxdxufWXs8ZH6xhC59BxDXs3fdyvYAGP04rH4RzuyC+CFei8NmtTCieyQjukfyn1P6ce5CBesOnmeNnG0IUa9OlyjWzJvLuRNHvXrM2G49uf6Bx5pUt6FpxjMyMqiqqqJXr15eja+tOLjpLBaLotfwWACO79hG4Zkcxt11b/0NRj5ijH7695/h9jdaLK7YEP9GzzZiQhxc1yeGXrFBxAQ7iAm5uEQFObBaJImIjqvTJQpfqZ5m/Pjx4wwfPvwb04zn5ORw3333MX/+fCyWjjfGwOPRHNx8lq6DowgI9gNg+7LFBIVH0GfU2PobBUTAiIdgwysQ3Q+u+09o4W/1lzvbWHvgHKv3n+Wf277ZAW5REBlUK3lckkiig/0I8rPhb7fib7fgb7fisF1cy5mKaOs6XaJo6jd/b2tomvHi4mKmTZvGiy++yOjRo30SX0vLPlRIaWEl4+4wLjEV5JzmWOZWxtzxbay2BvoCrv8JlJyFNS/CmZ1w69/AcWXzQF2N2mcbAOVVbs6XVHLuQgW5FyovLiUXtw+fvUBuSSVOd9NGbDlslprEUZ1MHDZjbbUorBaFRSlstbatFoXFYpYpY9uqFFarwqJAYa7NJKTqlBnlCuCScqPoYhuzCrULVN2XNfXrc7kc2NzUKDm1YbGh/tw1Itnrx+10icLXqqcZnzFjBk888QRaa2bOnMn999/PnXfe6evwWszBjDPYHeZMsUDm8s+wWG0MnXRjww3tATDzdYhPgRX/D948ArPfgcgerRD1NwX4WUmODCQ5MrDBelprisqdNUmkvMpNhdNDhdNNpctYV7iMssraZU6zzNzn9micbmNdvXi0ua1rldW8Bo/WaK3xaCMODaBBm3F5NGg02ixDG23MzZr46772/r+l8L7U5HBJFB1F7WnGlVKsW7eOvLw85s2bB8C8efNITU31bZBe5HK6ObL1HL3SjJliqyrK2b12JX1HjyMoPKLxAygFY/8D4gbCBw/C/10Pd7wNva5v+eCbSSlFeKAf4YF+9InreM9Dr0kkDSSQy+1q7k2Pkqsa11InXJIoWklD04zfe+9lOnM7iBO76s4Uu/fL1VSVlzX4qNN69foWPLYG3v02/P02mPwijH5Srkf4QO1LWs1o7dVYRMvreL2mos05sOniTLFaa7YvX0Jczz4k9Ol35QeL7AmPrIB+N8Hyn8C/ngBnufeDFkLUkEQhWlRFqZMTey7OFHty9w7yT58iber05o/2cYTAXQvh+p/Cjnfh7Ruh6LR3AxdC1JBEIVrUkW3n8Lh0zZQd25ctISAklH5jrr26A1ssMP6/YPY/4PwhmDsBTn599QELIb6h0ySK9j5raHuNv3qm2JiuIRSdO8vRrRmk3DAVm5+fd96g/zR4ZJUxZHbedNg6zzvHFULU6BSJwt/fn7y8vHb7x1ZrTV5eHv7+/r4O5YpcyK8wZ4qNQylF5hefgYKhk27y7hvF9odHV0OP62Dx07DkB81+Qp4Q4ps6xainpKQksrKyyM3N9XUozebv709SUpKvw7giBzPOANBnZDzOygp2r/6CPiPHEBIV7f03C4iAez6AVb80pvw4tw/uWgDBMd5/LyE6mU6RKOx2Oz16+OYGrc6seqbYsJgAdq3+gorSElIvnSXWmyxWmPSCcXPeJ981+i1mvwNdOs49KUL4Qqe49CRa35mjReRnl9JvVJwxJHbZYqK7didpwOCWf/Mhd8BDy41B/m9NgZ0ftPx7CtGBSaIQLWLL0uP4B9vpNzqB0/v3kHvi2NUNib1SXVLh0TWQOBw+egS++Bl43K3z3kJ0MI0mCqXUW0qpc0qp3bXKIpVSK5RSh8x1RK19P1ZKHVZKHVBKTalVPlwptcvc94oy/2IopRxKqffM8k1Kqe612swx3+OQUmqOtz60aFm5py5wYlceQ7+VjN1hZfuyJfgHBTPgmgmtG0hwDNz/CYx8FDb8L7w+HlY+D4dWQEVR68YiRDvWlDOKecDUS8p+BKzSWvcBVpmvUUoNBGYDg8w2f1VKWc02fwMeA/qYS/UxHwYKtNa9gT8CvzWPFQn8AhgFpAO/qJ2QRNu1bdkJ/PytDJmQyIW88xzK2MDgb03G7vDBqC2rHab9D9z6mjHB4Ib/hXfugN92h9euhaXPwd5PoKT9DnQQoqU12pmttV5X+1u+aQYwwdyeD6wFnjPLF2mtK4FjSqnDQLpS6jgQqrXeCKCUWgDcCiw12zxvHutD4C/m2cYUYIXWOt9sswIjubx75R9TtJaCM6Uc3naO4VO64Qi0s2XxUrTWpE728pDYK5V6t7FUlUHWZji5EU78G7bOh02vGXWi+kC3sdBtHHQbA+FdfRuzEG1Ec0c9xWmtcwC01jlKqVizPBGofXtsllnmNLcvLa9uc8o8lkspVQRE1S6vp00dSqnHMM5W6NpVfrl9advyE9hsFoZOTMbldLJz1XJ6DhtJWGy8r0Mz+AVCz/HGAsb9Fjk7jKRxYgPs+Rdsm2/sC0s2EkfXMUbyiO4jExCKTsnbw2Pr+y3SDZQ3t03dQq3nAnMBRowY0T7vqusAivPKObjpLIMnJBIQ4sfedaspKyq88lliW5PND5JHGss1zxgd3uf2GknjxAY4sgZ2vmfUDYwyzjICo8wlGgIjL74Oir64HRBhDNcVogNobqI4q5RKMM8mEoBzZnkWUPupGUlAtlmeVE957TZZSikbEAbkm+UTLmmztpnxilaQ+cVJUJA2yTir275sMZFdkug2pB3dx2CxQvwQYxn1HeOBC/lHjTOOU5vgwhkoy4PzB6E0D5yllzmQMpJFYJSRTOwBYPMHq5+xtplrqwNstRbrJdsWKyiLua69balbXlNmrX6EnRHDZdeWy+yr76Ooi5/pG2WXlF+2bT3/Ps0hZ3QNs/pBcGzj9a5QcxPFp8Ac4GVz/Umt8n8opf4AdMHotM7QWruVUheUUqOBTcD9wP9ecqyNwB3Aaq21VkotB35dqwN7MvDjZsYrWlhpUSV7/51D/9HxBEf4c3zHNs4cOcTEh55o38+EVgqiehnLsPu/ud9ZDmX5RvIoO19r21xKzxvrqlJjn7sKXBXGJS93JbjMxfPNZ3ELccUSR8Cjq7x+2EYThVLqXYxv9tFKqSyMkUgvA+8rpR4GTgJ3Amit9yil3gf2Ai7gu1rr6sHrT2CMoArA6MReapa/CSw0O77zMUZNobXOV0r9Cths1nuhumNbtD07Vp3C4/aQNqUbHrebtQveIDwugcHfmuzr0FqWPQDCEo3lang8dROHu9JIJtptXA6rvdb6kjKPuW2utcd89JxueK09xntXl1WrMyeavnzZN8ovdZl9zZ5zTa4qNyqwBabHoWmjnu6+zK6Jl6n/EvBSPeVbgG/clqu1rsBMNPXsewt4q7EYhW9VlDrZ/eVpeo+IIzw2kO3Ll5CXdZIZ//kzbHa7r8NrHywWsAQYiUeINkbuzBZXbeeaLJyVboZP7UZFSQkb3n+HroNT6DVilK9DE0J4gSQKcVWqKlzsXH2K7inRRCUGs/Gf71JZWsqE+x9t330TQogakijEVdmzLpvKMhcjbuxOfnYWmcuXMORbk4npJrP1CtFRSKIQzeZyuslceZKk/hHE9Qjly4VvYvPzY9yse30dmhDCiyRRiGbbvyGHsuIqht/YneM7tnF022ZG3zabwLBwX4cmhPAiSRSiWdxuD9uWnyS+ZygJvUJqhsOm3XiLr0MTQniZJArRLIc2n+VCfgXDb+zOrlXLycs6yXX3PSTDYYXogCRRiCumPZpty04QlRRMXHcH//7gHZIHpdB7xGhfhyaEaAGSKMQVO7I9l4IzZQyf2o1NHy+iouQCE+5/RIbDCtFBSaIQV0RrzdZlxwmPCyQywcn2ZcZw2NjuPX0dmhCihUiiEFfk5J58zp8qYdiUrqx75y1jOOxdMhxWiI5MEoVoMq01W5ceJzjSgSMwh6NbMxh922yCwuUJtUJ0ZJIoRJPlHC4k50gRQycmsW7hm4TFxctwWCE6AUkUosm2LD1BQIgdT+Uu8rJOMv5eGQ4rRGcgiUI0ydnjxZzam8+ga6PZ+NE/SB44hN4jx/g6LCFEK5BEIZpk27ITOAJtlJxfbwyHnSOzwwrRWUiiEI3Kyy7haGYuvYf7sXPlZwy5fpIMhxWiE5FEIRq1bdkJbA4reSeWm7PD3ufrkIQQrajRR6GKzq0ot5xDm8/SbXAZ+9dt5tpvP9Ck4bDny8/zzr53iAmIYVTCKHqG9ZRLVUK0U5IoRIO2fXECLJpzhz8nLC6eYTfNaLTNsuPLeOnrlyiqLEKjAYjyjyI9IZ1R8aNIT0gnOSS5pUMXQniJJApxWSUFlezfmENMl5Oc3HmSW37wkwaHwxZUFPDSppdYfnw5g6MGM2/qPBxWB5vPbGbTmU1sytnE0mNLAUgMTiQ9Pp30hHTS49OJDYxtrY8lhLhCkijEZWWuPInHVcG5o1+QNHAwvdMvPxx21clVvLDxBYqrinkq7Snm9L6b8pWrsUZGcUvKRGb2mYnWmmNFx9h0ZhMZORmsOrmKjw9/DECPsB6kx6czKmEUI+NGEu4vDz8Soq2QRCHqVV5SxZ71pwkJ20nuiRIm3F//cNiiyiJ+k/EbPjv6Gf0j+/P6Da+TsPEIJ565BVd2Tk09vx49CEhJIXxoCremDGX2NXegbVYO5B8g40wGm3I2sfjIYt478B4KRb/IfjWJY3jccILsQa358YUQtUiiEPXaseoUVRV5lOZuYMj1k4jr0esbddZlreP5Dc9TUFHAE0Of4D6dTt6Tz5O9YweO/v1J+OULKKuF8p07Kd+xk5KvvqLok08AUA4H/gMHEpmSwsyhKdyd8lO4Poa9eXvJOJNBRk4Gi/YvYsHeBfhb/bl34L08OPhBQv1CW/ufQohOT2mtfR2DV40YMUJv2bLF12G0a5XlLhb8ZAO6ajGVpcd56E9z64x0ulB1gd9t/h3/Ovwveof35qXezxD21mKKP/sMa3Q0sc88TdjMmSirtc5xtdY4T2dTsXMH5Tt2Ur5jBxV796KrqgCwRkcTkJJiLENTsAzsy66yI3x46EOWHltKqF8oDw95mLv7302ALaBV/02E6OiUUlu11iPq3SeJQlxq67Lj/PuDNThL/sm1336A9Bl31OzbcHoDP9/wc3LLc3m0933c/rWFovkLQGsiH3yQqEcfxRrc9MtEuqqKigMHKd+5gwrzzKPq+HFjp1L49epJ4IgRFE4bw/8WfcL60+uJDYjl8dTHubX3rdgtMteUEN4giUI02ZHt51i9YC8VhQsJCIYHfv83bH5+lDpL+f2W3/PBwQ/oGdKdl4pvwO/ND3Hnnif0ppuI/eEPsCcmeiUGd2Eh5bt2U75zB+U7d1K2KQNdUUHQ2LHkzRjHH2yryTy/g64hXfle2veY3H0yFiX3jgpxNSRRiEaVFlaybtFBjmbm4h94gMLTn3HzD35M31HjyMjJ4Ocbfk52STbfV5O47qMjVO0/gP/QFOJ+9CMC09JaNDZ3YSEF771Pwd//jis3F79evTh/y2j+J2oz+0uPMiByAE8Pe5qxXcbKTX1CNJMkCnFZ2qPZ++9sNvzzMC6Xk4TuJzi1exmx3Xsx/Sf/j1e2v8I/9v+DtKoE/isjDuv6LdgSEoj94Q8JnXZTq/5h1lVVFC9dSt68+VTu24c1MpLcqcP4Y7f9HOAMI+NH8vSwpxkaM7TVYhKio5BEIepVeLaMNX/fz+mD+YRGnaAsfx2lBXl0HzqMhJkTeGnf7zmfe4Kf7etH75UHUX5+RD/2KJEPPIDF399ncWutKduUQf68eZSsXYvy8yP3uoG82u8ke0KLuT75ep5Ke4reEb19FqMQ7Y0kClGH2+0hc8VJMhYfA88xLOprSvKyievZh/4zb2KF3sJ7e//B7XtDuf3LKtSFUsJum0nM009jj21bd1BXHj1G/oL5FP3rE3RFBflDu/H24PNkJFdyc+9beDL1SRKDvdN3IkRHJolC1Dh3opjVC/aTe+IgNuvXlBUeJywugZCJQ1nhv50tpzcx7KjiOxsCCc0uJjA9nbgfPYf/wIG+Dr1BroICCt97j/x33sGde54LSRG8O7SUDYMt3DZoNo8MeYSogChfhylEmyWJQuCsdJOx+Cjbv8jE49yAs+wQjtAQytPj+dx/C4mH8pl4OIDhB93YSyqwd+tK3LPPEjxxYrvqIPZUVVH82efkz5tH5YEDVIQ6+HSoi/UjA5iUegcj40eSFptGhH/jM+AK0ZlIoujkTu3NZ9X8DAqz1+Ku2oPys3FmgJVzpfsZfQhGHbbgX+bEEhRE8Le+ReiUyQRfdx3Kz8/XoTeb1pqyr78mb948Sr9ch8tuYXMfxcEEzbE40H27M7CrkTSGxQ4jKSSpXSVEIbxNEkUnVVHi5Mt3d7Bv/We4q7ajlZvC8BISc3MZdVgTUOFBBQcROnEiIVOmEjRuLBaHw9dhe13lkSPkz1/AhS/X4j57rqb8bKSFw3GaY/GKvK5hhKcMZ3CPUaTFpdEvoh82i8xwIzoPSRSdgNaasuIqCs6UkXe6mHPHcji46SsqijeidSU2Tykjjp0nqsSFDg4k7IbJhE6dQtDYsVja8ZnDlXKdP0/F3r1U7N1L+e49lOzeCWdqJY9wOBqvyEqwY+nfh/i0sQzuM5aU6BQC7YE+jFyIltXuE4VSairwZ8AKvKG1fvlydTt6ovB4NBfyKjh34jzZh06Qc/w4hWdyqCwqwOMsQnsK0Z4LgAeAsHI3Q05mE2izEDJxItHTbiFo9Oh2fVnJ21wFBUby2LOXwp1bKdu9G/uZvJr9uaFwLN5CSY8YrHGx2MMjCIiMJTAqluCoBMJjkogMjSXKP4oQvxC5S1y0S+06USilrMBBYBKQBWwG7tZa762vfntIFC6nm/LSci4UF1FaWERp8QXKiy9QfqGEitJSnOWVVJVV4KyswlVVhauiivKiC7jKy3C7ytCeItAVdY5p0VYcLiuBVU6Cy0sJKS8hQGkiR4+k2613Ezx6NKqBhw6JutxFRVTs20fRjm2cy9yEa98Bgs4UXbZ+uR+U+ENJgKIyyI4z2B9PSCCEhmALD8PiH4DFzw+rIwCrnwOrwx+bfwB2RwA2/wBsjgAcjiDs/gH4+Qdh9w/CERCE1eaHzWLHYrdhtdix2GzYrHYsyoJVWaVfRXhNQ4miPVyETQcOa62PAiilFgEzgHoTRXOd2Lubj174dQM1NKDMdd3/UmdbQ80DQDVoXVN2cXE3I0KFVfvj77bi8FhwYCfA4SY4zJ/QuDjCu3QjLKEbYQndscVEY4uJwRIUJH9ImskaFkbQ6NEEjR5NF54EwFNWhiu/AHdhIVWFeVw4n0NJbjblBblU5edhKSwguKiYkOISrGfL8Tt6Hv+yM1g9V/beTnMpa6COR5mLxVhr87W2KLQyfxpr1kYZUGutLnlt1KtW/dOsatW5WHJJWSv+iGn5eW5Qabdopi1c5fXjtodEkQicqvU6CxhVu4JS6jHgMYCuXbs2603sDn+sloAr+8Naq27tLaWqf8EUqnqfqrVPKSwWhcVmwWK1YLVZsdktWB027HY7dn8/7P52HIF+OAIDCAgOIr5nf6KTemOLiEBZ5NKGL1gCA/ELDISkRAKAsCa00R4PntJSdGUlzopyqipKqawooaqilKqKMpwVpTgrynFWlOGqrMBVWY67sgJXZQWeqgo8Hje4PWiPG22u8XjQHg+43cbaU12u0W638eWkzgI1X1pqv/aY6aC6vCboWhv6G4UX219S3GTNvYrRxq9+tAkJcS1y2PaQKOr7y13nJ+b/t3d3MXJWBRjH/w9UwGIJH5WvbmMLNgISFWIIIDFGNBZsWi9LIGmilxLQaICmCYnXGsALhTSAoV9jzAAABRNJREFUEkG4QNSG8NUAiVdWPgQslEqxBRarrTEokQRofLg4pzosMy8tlJ6z9Pklm533vLvbf97Zd87smemM7bXAWihLT+/lHznx5E9y2S9vei/fGjGRDjqIg+fNg3nzmAPkXTRiNpoNd02ngYUj21PAXxu1REQccGbDRPEIsETSYkmHACuBdY2bIiIOGN0vPdneJelS4H7K02Nvtv1046yIiANG9xMFgO17gHtad0REHIhmw9JTREQ0lIkiIiIGZaKIiIhBmSgiImJQ96/1tLck7QReeB8/Yj7wj32U80HovQ/6b+y9D9K4L/TeB301fsL2x8ft+NBNFO+XpEcnvTBWD3rvg/4be++DNO4LvffB7GiELD1FRMS7yEQRERGDMlG809rWAe+i9z7ov7H3PkjjvtB7H8yOxjxGERERw/IXRUREDMpEERERgzJRVJKWStosaYukq1r3AEhaKOlhSZskPS3p8jp+tKT1kp6rn49q3HmwpD9KurvTviMl3Snp2Xosz+mpUdJ36/W7UdLtkg5r3SfpZkk7JG0cGZvYJGl1PXc2S/paw8Yf1uv5KUm/lnRkq8ZxfSP7vi/Jkua36tsbmSgoN3TAT4ALgNOAiySd1rYKgF3A92yfCpwNfLt2XQU8aHsJ8GDdbulyYNPIdm99Pwbus30K8FlKaxeNkhYAlwGft3065aX0V3bQ93Ng6YyxsU31d3Il8On6PT+t51SLxvXA6bY/A/wZWN2wcVwfkhYCXwVeHBlrdQz3SCaK4ixgi+2/2H4DuANY0bgJ29ttP14vv0q5gVtAabulftktwDfaFIKkKeDrwI0jwz31HQF8EbgJwPYbtl+ho0bKy/1/VNIcYC7lHRyb9tn+HfDPGcOTmlYAd9h+3fZWYAvlnNrvjbYfsL2rbv6e8o6YTRonHEOAa4ErePtbOjc5hnsqE0WxAHhpZHu6jnVD0iLgDGADcJzt7VAmE+DYdmVcR/ml/+/IWE99JwE7gZ/V5bEbJR3eS6Ptl4EfUe5dbgf+ZfuBXvpmmNTU6/nzTeDeermLRknLgZdtPzljVxd9k2SiKDRmrJvnDUv6GPAr4Du2/926ZzdJy4Adth9r3TJgDnAmcL3tM4D/0H4p7H/qOv8KYDFwInC4pEvaVu217s4fSWsoS7e37R4a82X7tVHSXGANcPW43WPGurkNykRRTAMLR7anKH/+NyfpI5RJ4jbbd9Xhv0s6oe4/AdjRKO8LwHJJ2yjLdV+WdGtHfVCu22nbG+r2nZSJo5fGrwBbbe+0/SZwF3BuR32jJjV1df5IWgUsAy72//+jWA+NJ1PuEDxZz5kp4HFJx3fSN1EmiuIRYImkxZIOoTyotK5xE5JEWVvfZPuakV3rgFX18irgt/u7DcD2attTthdRjtlDti/ppQ/A9t+AlyR9qg6dDzxDP40vAmdLmluv7/Mpj0X10jdqUtM6YKWkQyUtBpYAf2jQh6SlwJXActuvjexq3mj7T7aPtb2onjPTwJn1d7R53yDb+Sh3Oi6kPEvieWBN657adB7lz8+ngCfqx4XAMZRnnTxXPx/dQeuXgLvr5a76gM8Bj9bj+BvgqJ4agR8AzwIbgV8Ah7buA26nPGbyJuUG7VtDTZQlleeBzcAFDRu3UNb6d58vN7RqHNc3Y/82YH7LY7inH3kJj4iIGJSlp4iIGJSJIiIiBmWiiIiIQZkoIiJiUCaKiIgYlIkiIiIGZaKIiIhBbwFQJsOv0ksqwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t, output[\"S\"][0,:], label=\"S1\")\n",
    "plt.plot(t, output[\"S\"][1,:], label=\"S2\")\n",
    "plt.plot(t, output[\"I\"][0,:], label=\"I1\")\n",
    "plt.plot(t, output[\"I\"][1,:], label=\"I2\")\n",
    "plt.plot(t, output[\"R\"][0,:], label=\"R1\")\n",
    "plt.plot(t, output[\"R\"][1,:], label=\"R2\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Model developer not interested in supporting stratified...\n",
    "\n",
    "But using the same `BaseModel`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# User (model developer) specifies...\n",
    "\n",
    "# ... state transitions in a subclass\n",
    "class SIR(BaseModel):\n",
    "\n",
    "    # state variables and parameters\n",
    "    state_names = ['S', 'I', 'R']\n",
    "    parameter_names = ['beta', 'gamma']\n",
    "    \n",
    "    @staticmethod\n",
    "    def integrate(t, S, I, R, beta, gamma):  # All variables and parameters... will be long list or arguments(!)\n",
    "        \"\"\"Basic SIR model\"\"\"\n",
    "        N = S + I + R\n",
    "        dS = -beta*S*I/N\n",
    "        dI = beta*S*I/N - gamma*I\n",
    "        dR = gamma*I\n",
    "        \n",
    "        return dS, dI, dR\n",
    "\n",
    "# ... parameters and initial conditions\n",
    "parameters = {\"beta\": 0.5, \"gamma\": 0.3}\n",
    "initial_states = {\"S\": [7900000], \"I\": [10], \"R\": [0]} # states as array|list\n",
    "\n",
    "# -> user initiates the model\n",
    "sir_model = SIR(initial_states, parameters)\n",
    "\n",
    "# -> user runs a simulation for a defined time period\n",
    "time = [0, 150]\n",
    "t, output = sir_model.sim(time)\n",
    "\n",
    "# -> user can do fit, mc,... using the model class instance `sir_model`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7feebbfcce50>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEDCAYAAAAcI05xAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU9b3/8ddnZrISQoCELQn7vgsRpShVwRaXqtWKIrhV66+9va1ttb229vbetrd3ab3WeutaXMsiKtR9b13ABQmbLGFfA0jCkhDIOjOf3x9nEgIEGGAm50zyeT4c58w5Zw7vJJN3Tk7OfI+oKsYYY7zL53YAY4wxJ2ZFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHmdFbYwxHhe3ohaRJ0WkRERWRrn+ZBFZLSKrRGRWvHIZY0yikXidRy0i44GDwLOqOvQk6/YDngcuUtX9ItJJVUviEswYYxJM3PaoVfUjYF/jeSLSR0TeEpHFIjJfRAZGFn0HeEhV90eeayVtjDERzX2M+nHgB6o6GrgbeDgyvz/QX0Q+FpHPRGRSM+cyxhjPCjTXPyQiGcBXgBdEpH52SqMc/YALgDxgvogMVdWy5spnjDFe1WxFjbP3XqaqI5tYVgx8pqp1wGYRWYtT3IuaMZ8xxnhSsx36UNUDOCV8LYA4RkQWvwRcGJmfjXMoZFNzZTPGGC+L5+l5s4FPgQEiUiwitwFTgdtEZDmwCrgysvrbwF4RWQ28D/xUVffGK5sxxiSSuJ2eZ4wxJjbsnYnGGONxcfljYnZ2tvbs2TMemzbGmBZp8eLFe1Q1p6llcSnqnj17UlhYGI9NG2NMiyQiW4+3zA59GGOMx1lRG2OMx1lRG2OMx0VV1CLy48jwoytFZLaIpMY7mDHGGMdJi1pEcoEfAgWR4Ur9wPXxDmaMMcYR7aGPAJAmIgEgHdgZv0jGGGMaO2lRq+oO4D5gG7ALKFfVd45eT0TuEJFCESksLS2NfVJjjGmlTnoetYi0xxmToxdQhjNM6TRVndF4PVV9HGe8aQoKCk7rfekP/n09wVD4ZIFOZ9OHnx7F5gRBxFm3fnn90KzOfDk8v9E8n08I+AR/5L7xY3/DtA+/D/w+n7OOCEl+ISM1QGZqEm1TA2SkBJAz/DiNMS1HNG94mQhsVtVSABGZhzOu9IwTPus0PPrhRqrqQsdd3lqGJfEJtE1NIjPNKe/MxtNpxz4e2KUt+R3S3Y5tjImTaIp6G3CuiKQDVcAEIC5vO1z9m+a9sEtTA1KpgkaW6RHztOEHRcN9ZF79+uEwBMNhQqqEwkowpIRVCYadx/W34BHTYYIh5WBNkIrqOg5UBTlQXceBqjoOVAcj93Vs2VPZMP9Q7bE/zHp2TOe8ftmc3y+HsX06kpmaFJfPmTGm+Z20qFV1oYi8CCwBgsBSIoc4El1ThxcOz/LuoYdgKExFtVPo+yvrWLptP/PX72Hekh3M+Gwbfp8wIq8d5/fL4fx+2YzIzyLJb6fMG5Oo4jLMaUFBgdpYH82vNhhmybb9LFi/h/kb9rCiuIywQkZKgLF9OnJ+v2zO65tNr+w2dgzcGI8RkcWqWtDkMivqlqusspZPNu5l/vo9zF9fSvH+KgBys9Kc0o4Ud1Z6sstJjTFW1AZVZeveSuZv2MP8daV8unEvFTVB2qYE+Ovt5zAyP8vtiMa0albU5hjBUJhl28v4yfPLKa+q47k7zmVQ10y3YxnTap2oqO0vTK1UwO+joGcHZt5+DmlJfm584nM2lR50O5YxpglW1K1cfod0Ztx+DqrKtOkLKd5f6XYkY8xRrKgNfTtl8OxtY6ioCTJt+kJKKqrdjmSMacSK2gAwpFs7nr51DCUVNdw4/XP2H6p1O5IxJsKK2jQY3aM9028qYPPeQ9z81OdUVNe5HckYgxW1OcpX+mbz8A2jWL3zALc9U0hVE29XN8Y0Lytqc4yJgztz/3UjWbRlH9+dsZiaoJW1MW6yojZNumJEN/776mF8uK6UO2cvO/nws8aYuLGiNsd13dnd+dfLB/PWqi/52dwvCIdbyTizxnhMNMOcmlbstvN6cagmyP3vrqNNcoDfXDnEBnQypplZUZuT+sFFfTlUE+SxjzbRJiXAv0waYGVtTDOyojYnJSLcc8lADtYEefTDjbRNDfD9C/u6HcuYVsOK2kRFRPjtlUOprA3xh7fXkp7s59ZxvdyOZUyrEM3FbQcAcxrN6g38SlUfiFsq40k+n/CHbw2nsjbIr19dTZuUAJML8t2OZUyLd9KzPlR1raqOVNWRwGigEvhb3JMZTwr4fTw45SzO75fNPXO/4PUvdrkdyZgW71RPz5sAbFTVrfEIYxJDSsDP4zcWMLpHe+58bikbbXhUY+LqVIv6emB2UwtE5A4RKRSRwtLS0jNPZjwtLdnPQ1NHAfD8ou0upzGmZYu6qEUkGbgCeKGp5ar6uKoWqGpBTk5OrPIZD+vUNpWLBnZi7pId1Nk7F42Jm1PZo74EWKKqu+MVxiSeyQX57DlYw/trStyOYkyLdSpFPYXjHPYwrdcFA3LIaZvC84XFbkcxpsWKqqhFJB24GJgX3zgm0QT8Pq4elcv7a0vsyjDGxElURa2qlaraUVXL4x3IJJ7JBfmEwsrfluxwO4oxLZKNnmfOWJ+cDAp6tGdO4XZUbYQ9Y2LNitrExOSCfDaVHmLJtv1uRzGmxbGiNjFx6fCupCf7eX6R/VHRmFizojYxkZES4PLhXXnti50cqgm6HceYFsWK2sTM5IJ8DtWGeH2Fjf9hTCxZUZuYGd2jPb2z2/BCob2l3JhYsqI2MSMiXFuQz6It+9lkAzUZEzNW1CamrhmVi98nvLDY/qhoTKxYUZuY6pSZyoUDcpi7uJigDdRkTExYUZuYu7Ygn5KKGj5cZ8PdGhMLVtQm5i4a2InsjGSetz8qGhMTVtQm5pL8Pq4elcffi0rYc7DG7TjGJDwrahMX147OI2gDNRkTE1bUJi76dW7LWd2zeN4GajLmjFlRm7iZXJDP+pKDLNte5nYUYxKaFbWJm8uHdyUtyW9XfzHmDEV7hZcsEXlRRNaISJGIjI13MJP42qYmcemwrry6fCeVtTZQkzGnK9o96j8Bb6nqQGAEUBS/SKYlmVyQx8GaIG+u+NLtKMYkrJMWtYhkAuOBJwBUtVZV7aCjicqYXh3o2THdzqk25gxEs0fdGygFnhKRpSIyXUTaHL2SiNwhIoUiUlhaau9IM476gZoWbt7Hlj2H3I5jTEKKpqgDwCjgEVU9CzgE3HP0Sqr6uKoWqGpBTk5OjGOaRHbNqDx8Ai8str1qY05HNEVdDBSr6sLI4xdxituYqHRpl8pX++fw4uJiQmE7p9qYU3XSolbVL4HtIjIgMmsCsDquqUyLM7kgn90HavhovR0WM+ZURXvWxw+AmSLyBTAS+M/4RTIt0YRBnenQJtmu/mLMaQhEs5KqLgMK4pzFtGDJAR/fPCuXZz/dwt6DNXTMSHE7kjEJw96ZaJrN5IJ86kLKS8t2uh3FmIRiRW2azYAubRmR144XbKAmY06JFbVpVpPPzmfNlxWs2FHudhRjEoYVtWlW3xjRjZSAjzmL7I+KxkTLito0q8zIQE2vLNtJVW3I7TjGJAQratPsri3Io6ImyNurbKAmY6JhRW2a3bm9OtK9gw3UZEy0rKhNs/P5hGtH5/HJxr1s21vpdhxjPM+K2rjimtF5iMCLNlCTMScV1TsTjYm1bllpnN/PGajpzon98fvE7UjmKKpKSEOENdxwA5xpwqgqqnp4GiWsh6cb7htNOxumYbp+WcN0o3WOyHL0jEi+Y+Y1sV5UH+tpPu9oAQnQO6t3TLZ1xHZjvkVjonRdQT7fn7WEjzfsYXx/Gxr3ZA7WHuRA7QGqglVU1lVSGaw8fN9ouqqu6shldZVUBauoC9cR0hDBcJBgONgwHQqHCGrT88yp6ZjakQ+u+yDm27WiNq6ZOLgTWelJzCncbkUNVAWr2HlwJzsO7nBuFTsaposPFlNRWxHVdlL8KaQH0klPSictkNZwn+HLIOALEJAAfp+fgC+AX/wk+ZLwi79hXkACzjKf35kvfnziQ0Tw4dwL0jCvYRpx1mk0LTi/KdVPi0Qec/g3qCPWO846kRnHOGad48yLRv2/eyaSfclnvI2mWFEb16QE/Fw1MpdZC7ex/1At7dvE50XuFapK8cFiiiuKGwq5+GBxQynvrd57xPop/hS6ZXQjNyOX4TnDyc3IJSsli7SkNKeII2Xc+D41kErAZ9/WLY19RY2rJhfk8/QnW3h52Q5uGdfL7ThxsalsE69vfp03N7/J9orDfzz1i58ubbqQl5HHV/O/Sm5G7hG37LTsmOzlmcRnRW1cNbhbJsNy2/F8YXGLKupdB3fx5pY3eWPTG6zdvxaf+BjTZQw3D76Z3lm9yc3IpVN6J9v7NVGxV4lx3eSCPP715VWs3FHO0Nx2bsc5bfur9/POlnd4Y/MbLClZAsDw7OHcM+Yevt7z62SnZbuc0CSqqIpaRLYAFUAICKqqXUTAxMwVI3L57etFPF+4PeGK+lDdIf6x7R+8uflNPt35KUEN0rtdb35w1g+4pOcl5Gfmux3RtACnskd9oaruiVsS02q1S09i0pAuvLR0B7+4dBCpSX63I51QbaiWBTsW8ObmN/lg+wdUh6rp2qYrNw25iUt7XUr/9v3t2LKJKTv0YTzhurPzeWX5Tt5ZvZsrRnRzO06TVJWHlj3ErDWzqKitoH1Ke67seyWX9rqUkZ1G4hN7o6+Jj2iLWoF3RESBx1T18aNXEJE7gDsAunfvHruEplUY27sjuVlpPL9ouyeLWlX53cLfMWftHCZ2n8jV/a7m3G7nkuRLcjuaaQWi3QUYp6qjgEuA74vI+KNXUNXHVbVAVQtycuzNC+bU+HzCtQV5fLxxD9v3eWugJlXlPxf+J3PWzuHWIbdy/wX3c37e+VbSptlEVdSqujNyXwL8DRgTz1CmdfrW6DwA5i4pdjnJYarKf33+Xzy39jluGXILPx79Yzv+bJrdSYtaRNqISNv6aeBrwMp4BzOtT177dM7rm80LhcWEwu5f/FZV+Z9F/8PsNbO5afBN/GT0T6ykjSui2aPuDCwQkeXA58DrqvpWfGOZ1mrqOd3ZUVbFK8t3uJpDVfn9ot8zs2gmNw6+kbsL7raSNq456R8TVXUTMKIZshjD1wZ3YVDXTB54bz2XD+9Gkr/5z6SoL+kZRTOYNmgaPy34qZW0cZWdT2Q8xecT7rq4P1v3VjJ3cfMfq1ZV7iu8r6Gkf3b2z6ykjeusqI3nTBjUiZH5WTz49/XUBJvvSuWqyv8W/i/Prn6WGwbeYCVtPMOK2niOiHD31waws7ya5z5vnkt1qSr3L76fZ1Y/w5SBU7hnzD1W0sYzrKiNJ43r25FzenXgz+9voKo2vnvVqsofl/yRp1c9zXUDruPnY35uJW08xYraeJKIcNfXBlBaUcOzn26J27+jqjyw5AGeWvkU1w24jnvPuddK2niOFbXxrDG9OjC+fw6PfriRiuq6mG9fVXlw6YM8ufJJJvefzC/O+YWVtPEkK2rjaXdd3J/9lXU89fGWmG5XVfm/pf/H9BXTubb/tdx77r02qJLxLHtlGk8bkZ/FxYM785ePNlFWWRuTbdaX9F9W/IVr+l3DL8/9pZW08TR7dRrP+8nF/TlYG+Qv8zed8bbqhyqtL+lfjf2VlbTxPHuFGs8b1DWTy4d346mPt7DnYM0Zbevh5Q/z2BePcXW/q62kTcKwV6lJCD+a2I/quhCPfrDxtLfxyLJHeHT5o1zV9yr+bey/WUmbhGGvVJMQ+uRkcPWoPP762Va+LK8+5ec/svwRHl7+MFf2uZJff+XXVtImodir1SSMOyf0IxRW/vz++lN63mPLH+PhZQ9zRZ8rrKRNQrJXrEkY+R3Sue7sfOYs2h71VWDW7lvLn5f9mct7X85vvvIb/D5vXzjXmKZYUZuE8s8X9UVEePDv0e1VzyyaSVogjXvG3GMlbRKWFbVJKF3bpTHtnB7MXVLMptKDJ1x3X/U+Xt/0Ot/o/Q3apbRrpoTGxF7URS0ifhFZKiKvxTOQMSfzvQv6kBLw88B7J96rnrtuLrXhWm4YdEMzJTMmPk5lj/pOoCheQYyJVk7bFG4Z15NXv9jJmi8PNLlOXbiO59Y+x9iuY+mT1aeZExoTW1EVtYjkAZcB0+Mbx5jo/L/xvclIDvDHd9c1ufy9re9RUlnCtMHTmjmZMbEX7R71A8DPgPDxVhCRO0SkUEQKS0tLYxLOmOPJSk/m9vN78/aq3awoLj9m+cyimXRv253zcs9zIZ0xsXXSohaRy4ESVV18ovVU9XFVLVDVgpycnJgFNOZ4vn1eT7LSk7jvnbVHzF+5ZyXLS5dzw6Ab7Jxp0yJE8yoeB1whIluA54CLRGRGXFMZE4W2qUl896t9+HBdKYu27GuYP7NoJm2S2nBlnytdTGdM7Jy0qFX156qap6o9geuBf6iqHfgznnDT2B5kZ6Rw39trUVVKK0t5a8tbXNX3KjKSM9yOZ0xM2O+FJqGlJwf4/oV9WLh5H59s3MsL614gFA4xZeAUt6MZEzOnVNSq+oGqXh6vMMacjiljutO1XSp/eGclc9bO4fy88+mR2cPtWMbEjO1Rm4SXmuTnBxf1Y2XZfPZV72PqwKluRzImpqyoTYvwrdG5ZHT6lECoC+d0OdftOMbElBW1aRFW71tBKGk7B0vO5Z3Vu92OY0xMWVGbFmFG0QzaJrclP+k87n93HaGwuh3JmJixojYJ78tDX/Le1ve4pt813HXxcNaXHOSV5TvcjmVMzFhRm4Q3Z+0cFOX6gddzydAuDOqayQPvracudNwRD4xJKFbUJqFVB6t5cd2LXJh/IbkZufh8wl0X92fr3krmLi52O54xMWFFbRLaG5vfoKymjKmDDp+SN2FQJ0bkZ/Hg39dTEwy5mM6Y2LCiNglLVZlRNIP+7ftT0LmgYb6IcPfX+rOzvJrnPt/uYkJjYsOK2iSswt2FrN+/nqmDpiIiRyw7r282Y3p14M/vb6Cq1vaqTWKzojYJa8bqGWSlZHFpr0uPWebsVQ+gtKKGv362pfnDGRNDVtQmIRVXFPNB8Qd8q/+3SA2kNrnOmF4dGN8/h0c+2MjBmmAzJzQmdqyoTUJ6bs1zCMJ1A6474Xp3Xdyf/ZV1TP3LZ6ze2fT1FY3xOitqk3Aq6yqZt34eE3tMpEubLidcd0R+Fg/dMIri/VVc8ecF/M9ba6ius2PWJrFYUZuE8+rGV6moq2DaoOiuX3HZ8K6895OvctVZuTzywUYmPfARn2zYE+eUxsSOFbVJKGENM3PNTIZ0HMKInBFRP699m2Tuu3YEM28/BwVumL6Qn76wnLLK2viFNSZGorm4baqIfC4iy0VklYj8ujmCGdOUz3Z+xubyzU2ekheNcX2zeftH4/nuV/swb+kOJt7/Ia8s34mqDeJkvCuaPeoa4CJVHQGMBCaJiA34a1wxo2gGHVM78vWeXz/tbaQm+bnnkoG88s/j6JaVxg9nL+W2ZwrZUVYVw6TGxE40F7dVVT0YeZgUudnuh2l2Ww9sZf6O+UweMJlkf/IZb29It3b87Z/G8cvLBvHpxr1cfP+HPLlgsw2RajwnqmPUIuIXkWVACfCuqi5sYp07RKRQRApLS0tjndMYZhXNIuALMHnA5Jht0+8Tbj+/N+/8eDxn9+zAb15bzdWPfMKaL+1UPuMdURW1qoZUdSSQB4wRkaFNrPO4qhaoakFOTk6sc5pW7mDtQV7a8BKTek4iOy075tvP75DO07eezZ+uH0nxvkouf3ABv7dT+YxHnOpVyMuAD4BJcUljzHG8tOElKoOVR4ySF2siwpUjcxtO5Xu4/lS+jXYqn3FXNGd95IhIVmQ6DZgIrIl3MGPqhTXMrDWzGJEzgqHZx/wyF3PHnMr3l4X87EU7lc+4J5o96q7A+yLyBbAI5xj1a/GNZcxh84vns71ie9RvcImVcX2zeetO51S+uUucU/letVP5jAskHi+6goICLSwsjPl2Tev0nXe+w6byTbx1zVsk+ZJcybBqZzk/n7eCL4rL6dspg+F57RiW69wGd8skPTngSi7TcojIYlUtaGqZvbqMp20s28hnuz7jh2f90LWSBudUvnnf+wqzP9/GP9aU8NG6Pcxb4lxA1yfQJyeDYVbeJk7slWQ8bWbRTJJ9yVzT/xq3oxDw+7hxbE9uHNsTVWX3gRpW7ChnxY5yVu4oZ/76Jso7tx1Dc9sxLK8dg7tm0ibFvuXMqbNXjfGs8ppyXt34Kpf1vowOqR3cjnMEEaFLu1S6tEvl4sGdG+bvPlDNiuLD5b1gwx7mLd0ReQ706JBOhzbJZKUnk5WWRLv0JLLSkmnfJol2aUkN87Mi89umBvD5Tv2t8qZlsaI2njVv/TyqQ9VxPSUv1jpnptJ5cCoTG5V3yYHqhj3v9SUHKa+so6SimnW7KyivrKPiBBc18AlkpiXRPj2ZzLQk0pJ8JAf8JPuF5ICPJL+PZL+PpIBznxy5T4pMJ/mFlMh6Pp/gF8HnA59Iw83vc37w1C9rmK5/jCACgrPs8DQQWeZMRZZHlkmjZUdrPN95RhPzT/DzqfFzTrTt5lL/T/p8Qp+cjJhv34raeFIwHGT2mtkUdC5gQIcBbsc5I50yU5mQmcqEQZ2bXF4XClNeVUdZZR3lVbWUVTrTZVV1lFXWHjFdU+esWxcMUxsKUxcKUxt07muChx/bu+DdkZ2RQuEvJ8Z8u1bUxpM+2P4Buw7t4l/O/he3o8Rdkt9HdkYK2RkpMdtmKKzURsq8vshDYUUVQqqEVQmHlbBCWPWYZapKKOwsC6uC8x/OpEbuaThVUSP/a1imRw4I1PjssiPnN06tx5l/pBP9DHLjzEltlCgl4I/Lv2FFbTxpZtFMurXpxgX5F7gdJSH5fUJasp804lMcpnnZhQOM56zZt4bC3YVMGTgFv8+KxhgrauM5s4pmkRZI45v9vul2FGM8wYraeMq+6n28vul1vtH7G7RLaed2HGM8wYraeMrza5+nNlzLDYNucDuKMZ5hRW08o7KukplFMxmfN54+WX3cjmOMZ1hRG8+Yu34uZTVlfGfYd9yOYoynWFEbT6gL1fHMqmcY3Xk0IzuNdDuOMZ5iRW084bVNr7G7cje3D7vd7SjGeI4VtXFdKBziiZVPMKjDIMZ1G+d2HGM8J5pLceWLyPsiUiQiq0TkzuYIZlqP97a9x9YDW7lt2G2IGyPqGONx0byFPAjcpapLRKQtsFhE3lXV1XHOZloBVeWJFU/QM7MnE7vHfjAbY1qCk+5Rq+ouVV0Sma4AioDceAczrcPHOz+maF8R3x76bXu7uDHHcUrHqEWkJ3AWsLCJZXeISKGIFJaWlsYmnWnxpq+YTuf0zlze+3K3oxjjWVEXtYhkAHOBH6nqgaOXq+rjqlqgqgU5OTmxzGhaqKUlS1m8ezE3D7mZJL9710M0xuuiKmoRScIp6ZmqOi++kUxrMX3FdLJSsrimn/vXQzTGy6I560OAJ4AiVb0//pFMa7B231o+Kv6IqYOmkp6U7nYcYzwtmj3qccCNwEUisixyuzTOuUwL98TKJ0gPpDNl4BS3oxjjeSc9PU9VF8AJriRpzCnadmAbb295m5sH32xDmRoTBXtnoml2T616ioAEuHHwjW5HMSYhWFGbZlVSWcLLG17myr5XkpNuZwcZEw0ratOsnl31LCENcevQW92OYkzCsKI2zaa8ppzn1z3PpJ6TyG+b73YcYxKGFbVpNrOKZlEVrOK2Ybe5HcWYhGJFbZpFZV0lM9fM5IK8C+jfvr/bcYxJKFbUplm8uO5FymvKbW/amNNgRW3irjZUyzOrnqGgc4FdZsuY02BFbeLu1Y2vUlJVYhetNeY0WVGbuAqFQzy58kkGdRjE2G5j3Y5jTEKyojZx9e7Wd9lWsY3vDP+OXWbLmNNkRW3iRlWZvmI6PTN7MqH7BLfjGJOwrKhN3CzYsYC1+9fy7aHfxif2UjPmdNl3j4mb6Sum06VNF7vMljFnyIraxMWS3UtYUrKEW4bcYpfZMuYMWVGbuJi+YjrtU9pzdb+r3Y5iTMKL5lJcT4pIiYisbI5AJvGt2beG+TvmM23wNNICaW7HMSbhRbNH/TQwKc45TAvyxIonaJPUhusHXu92FGNahJMWtap+BOxrhiymBdh6YCvvbH2HyQMmk5mc6XYcY1qEmB2jFpE7RKRQRApLS0tjtVmTYJ5a6Vxm66bBN7kdxZgWI2ZFraqPq2qBqhbk5Nglllqj3Yd28/LGl/lmv2+SnZbtdhxjWgw768PEzLOrn0VVuWXILW5HMaZFsaI2MVFWXcYL617gkl6XkNc2z+04xrQo0ZyeNxv4FBggIsUiYiO/m2PMWhO5zNZQe3kYE2uBk62gqlOaI4hJXIfqDjGzaCYX5F9A3/Z9z2xjB0uhZBV0HwuBlNgENCbBnbSojTmZF9e9yIHaA9w+7PbT30jFbvjkQVj0BASrIDULhl0LZ02FriPBhkg1rZgVtTkj9ZfZGtNlDCNyRpz6Bip2w8d/gsInIVQDwybDgEug6BVY8iws+gt0GuIU9vDroI2dTWJaHytqc0Ze3vgypVWl/Md5/3FqTzywyynoxU9BqM4p4fF3Q8c+zvIhV0HVflg5F5bOhLd/Ae/+CvpPgpFTod/FYIM9mVbCitqctmA4yFMrn2JIxyGM7RrlZbYO7IQFD8DipyEchBFTYPxd0KH3seumtYezb3duu1fDspnwxRxY8xq0yXHK/axp0GlQTD8uY7zGitqclrLqMn772W/ZXrGdP17wx5NfZqt8Byz4Iyx5BjTsFPT5d0GHXtH9g50Hw9d/BxP/Hda/65T2wkfh0z9Dt1Ew8gYY9i2n3I1pYURVY77RgoICLSwsjPl2jTd8sP0D/v2Tf6e8tpzvj/w+tw297fhFXV4M8++HpX91CnrkVKeg2/c48yCH9sAXzzulvXsl+FNg4GXO8ezeF4LPf+b/hjHNREQWq2pBU8tsj9pEraK2gt8v+j0vbXiJ/u3789jFjzGgw4CmVy7bFinoGc7js6bB+T+BrO6xCxB2p3IAAAqtSURBVNQmG8b+E5z7Pdi13CnsFS/AqnmQmQsjrnd+MNQf9zYmQdketYnKZ7s+41cf/4rdlbu5behtfG/E95q+csv+rTD/f2HZLOeUurNuhPN+DFn5zRM0WANr33D+/Q3vOXvx+ec6h0aGXAWp7ZonhzGn6ER71FbU5oQq6yp5YMkDzF4zm56ZPfndeb9jeM7wY1fcv6VRQftg1M1w3o+gnYtvJz+wC754zjlrZO96CKQ6h0ZGTHEOjfjtF0rjHVbU5rQsK1nGvQvuZVvFNqYNmsYPR/3w2Cu27NsM8++D5c+B+GH0LU5BZ3ZzJXOTVGHHElg+G1a+6Jz2l9HZeUPNiCnQZajbCY2xojanpjZUy0PLHuLpVU/TJb0Lvx33W8Z0HXN4hWANrHvL2Xte/65zPvPoW2HcnZDZ1b3g0QjWwPp3nB8s695yThHsMswp7GHXQkYntxOaVsqK2kStaG8Rv1jwCzaUbeCaftdwd8HdZCRnOHulO5c65Vy/V9q2q/MHu3O+C227uB391B3a67yhZvks52MTP/Sd6HxMAy6FpFS3E5pWxM76MCcVDAeZvmI6jy1/jPap7XlowkOMzxsPFV86428smw2lRZHjvJfDyCmJfwpcm45wzh3OrWSNczx7+RxY/zaktIOh33T2tPPPsbFGjKtsj9qwqWwT9y64l5V7V3JJr0u4d9RdtNv6yVFnTpzjlNaQb0JaltuR4yccgs0fOYdGil6Bukpo38v52EdcB+17up3QtFB26MM0qTZUy+w1s3lwyYOkJ6Xzy3438PVd651DG9Xlh89FHnEDZJ/h8KWJqKYCil51fmBtme/M6zEO+k6Ajv0gu5/z1ncbjtXEgBV1K1cTqmFz+WY2lm1suG0q38S2im2ENcwF6fn8W0kJ2aXrIZAGg77hnHfca3xiH9qIpbJtzrsgv5gDe9Ydni8+5008HftGyrvv4enMbnbIxETtjItaRCYBfwL8wHRV/e8TrW9F7Y6qYBVbyrewoWwDm8o3sXH/OjbuW09x5ZeEcb7OfoTukkyfEPSpqWZ4+R7Or6pCuo91ynnwVZCa6fJH4nHVB2DvBti70Tk/e8/6w4/rDh1eLyndeVdkx35OeWf3cw6jpLZzPscpmZCUZmVugDP8Y6KI+IGHgIuBYmCRiLyiqqtjG7P1CmuYmlANNXXVVNcdpLr2INW1B6iurXQe1x2kuq6S6rpKauqqqA5WUh2spjpURVVdFTuqdrOhcjc7Qoeo/7EbUKVHXR0Da+u4rC5I77o6+tbW0UNSSGrb1TmNrlM3GNQXhl5tb7M+FamZkDvKuTWm6owOuHdDpMA3ONM7l8Dql5xj/UfzBSClrVPa9eV9xHTbI+cHUsCf7LxZx58MvqTDjxtPNyyL3HxJzm9H4rMfDAkomrM+xgAbVHUTgIg8B1wJxLyor3tqJNVNvZhPw4l+T4jmYI9G/q9HPG50L0duSAXCODdFI/cQapgXuZem5p3eN05AlVRVugSDDK0NcoUk0yc5i77pXcnP7EFSZq7z63fbrofvbW85fkSgXa5z6/3VI5cFa5w3B5VtdfbIa8qdY+DVB6DmwJHT5cVQEpmuPgAaikNW3wlucuRjJFLu0qjkG82r/9iRyMOj5h39OTr84ATLmlh+wnWjfF68pXeEb78Z881GU9S5wPZGj4uBc45eSUTuAO4A6N799Abe6ZXcntpwMGY/8OUEX7ATf5nlqMdHzq3PJw1LDj/2ieBD8IkPHyAi+PEdXiaCEJkn9esJKb4kUv0ppAZSnHt/KqmBVFICaaQF0kgNpJGS1IbUpDakJqWRkpRBUlK6s7eU3sF5p50NpO9dgRToNNC5nQpV58yT+iIP1UCoFkJBCNcdng7VRh7X32qdN/M0TNc529LwCW4nWq407JnUTzfsqGijeXrUvCM+mCM/ruMta3L5CdaN+nnNIE47QtEUdVOddsxnQ1UfBx4H5xj16YT576nvn87TjGm5RCC5jXNLxDcVmZjwRbFOMdB46LM8YGd84hhjjDlaNEW9COgnIr1EJBm4HnglvrGMMcbUO+mhD1UNisg/A2/jnJ73pKquinsyY4wxQJRjfajqG8Abcc5ijDGmCdEc+jDGGOMiK2pjjPE4K2pjjPE4K2pjjPG4uIyeJyKlwNbTfHo2sCeGcWLN6/nAMsaC1/OB9zN6PR94K2MPVc1pakFcivpMiEjh8UaQ8gKv5wPLGAtezwfez+j1fJAYGcEOfRhjjOdZURtjjMd5sagfdzvASXg9H1jGWPB6PvB+Rq/ng8TI6L1j1MYYY47kxT1qY4wxjVhRG2OMx3mmqEVkkoisFZENInKP23kARCRfRN4XkSIRWSUid0bmdxCRd0VkfeS+vcs5/SKyVERe82i+LBF5UUTWRD6XY72UUUR+HPn6rhSR2SKS6nY+EXlSREpEZGWjecfNJCI/j3zvrBWRr7uY8Q+Rr/MXIvI3EclyK2NT+Rotu1tEVESy3cp3KjxR1I0uoHsJMBiYIiKD3U0FQBC4S1UHAecC34/kugf4u6r2A/4eeeymO4GiRo+9lu9PwFuqOhAYgZPVExlFJBf4IVCgqkNxhvK93gP5ngYmHTWvyUyR1+T1wJDIcx6OfE+5kfFdYKiqDgfWAT93MWNT+RCRfJyLdW9rNM+tz2FUPFHUNLqArqrWAvUX0HWVqu5S1SWR6QqcgsnFyfZMZLVngKvcSQgikgdcBkxvNNtL+TKB8cATAKpaq6pleCgjznC/aSISANJxrmDkaj5V/QjYd9Ts42W6EnhOVWtUdTOwAed7qtkzquo7qhqMPPwM54pQrmQ8zucQ4I/AzzjykoKufA6j5ZWibuoCurkuZWmSiPQEzgIWAp1VdRc4ZQ50ci8ZD+C86Bpfvt1L+XoDpcBTkcMz00WkjVcyquoO4D6cvatdQLmqvuOVfEc5Xiavfv98G6i/JLcnMorIFcAOVV1+1CJP5DserxR1VBfQdYuIZABzgR+p6gG389QTkcuBElVd7HaWEwgAo4BHVPUs4BDuH4ppEDnOeyXQC+gGtBGRae6mOmWe+/4RkXtxDh3OrJ/VxGrNmlFE0oF7gV81tbiJeZ7pIK8UtWcvoCsiSTglPVNV50Vm7xaRrpHlXYESl+KNA64QkS04h4suEpEZHsoHzte2WFUXRh6/iFPcXsk4EdisqqWqWgfMA77ioXyNHS+Tp75/RORm4HJgqh5+o4YXMvbB+YG8PPI9kwcsEZEuHsl3XF4pak9eQFdEBOfYapGq3t9o0SvAzZHpm4GXmzsbgKr+XFXzVLUnzufsH6o6zSv5AFT1S2C7iAyIzJoArMY7GbcB54pIeuTrPQHnbxFeydfY8TK9AlwvIiki0gvoB3zuQj5EZBLwL8AVqlrZaJHrGVV1hap2UtWeke+ZYmBU5DXqer4TUlVP3IBLcf5KvBG41+08kUzn4fz68wWwLHK7FOiI81f39ZH7Dh7IegHwWmTaU/mAkUBh5PP4EtDeSxmBXwNrgJXAX4EUt/MBs3GOmdfhFMptJ8qE8yv9RmAtcImLGTfgHOut/3551K2MTeU7avkWINvNz2G0N3sLuTHGeJxXDn0YY4w5DitqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxOCtqY4zxuP8PvpAo84fDzDwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(t, output[\"S\"].flatten())\n",
    "ax.plot(t, output[\"I\"].flatten())\n",
    "ax.plot(t, output[\"R\"].flatten())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Age based SEIRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SEIRSAge(BaseModel):\n",
    "\n",
    "    # ...state variables and parameters\n",
    "    state_names = ['S', 'E', 'I', 'A', 'M', 'C', 'Cmirec', 'Cicurec', 'Mi', \n",
    "                   'ICU', 'R', 'D', 'SQ', 'EQ', 'IQ', 'AQ', 'MQ', 'RQ']\n",
    "    parameter_names = ['beta', 'sigma', 'omega', 'zeta', 'a', 'm', 'da', 'dm', 'dc', 'dmi', 'dICU', 'dICUrec', \n",
    "                  'dmirec', 'dhospital', 'maxICU', 'totalTests', 'psi_FP', 'psi_PP', 'dq']\n",
    "    parameters_stratified_names = ['h', 'c', 'm0','mi']  \n",
    "    stratification = 'nc'\n",
    "    \n",
    "    # ..transitions/equations\n",
    "    @staticmethod\n",
    "    def integrate(t, S, E, I, A, M, C, Cmirec, Cicurec, Mi, ICU, R, D, SQ, EQ, IQ, AQ, MQ, RQ, \n",
    "                  beta, sigma, omega, zeta, a, m, da, dm, dc, dmi, dICU, dICUrec, \n",
    "                  dmirec, dhospital, maxICU, totalTests, psi_FP, psi_PP, dq, h, c, m0, mi, Nc):  \n",
    "        \"\"\"Basic SIR model\"\"\"\n",
    "\n",
    "        # Model equations\n",
    "        Ctot = C + Cmirec + Cicurec\n",
    "        # calculate total population per age bin using 2D array\n",
    "        N = S + E + I + A + M + Ctot + Mi + ICU + R + SQ + EQ + IQ + AQ + MQ + RQ\n",
    "        # calculate the test rates for each pool using the total number of available tests\n",
    "        nT = S + E + I + A + M + R\n",
    "        theta_S = totalTests/nT\n",
    "        theta_S[theta_S > 1] = 1\n",
    "        theta_E = totalTests/nT\n",
    "        theta_E[theta_E > 1] = 1\n",
    "        theta_I = totalTests/nT\n",
    "        theta_I[theta_I > 1] = 1\n",
    "        theta_A = totalTests/nT\n",
    "        theta_A[theta_A > 1] = 1\n",
    "        theta_M = totalTests/nT\n",
    "        theta_M[theta_M > 1] = 1\n",
    "        theta_R = totalTests/nT\n",
    "        theta_R[theta_R > 1] = 1\n",
    "        # calculate rates of change using the 2D arrays\n",
    "        dS  = - beta*np.matmul(Nc,((I+A)/N)*S) - theta_S*psi_FP*S + SQ/dq + zeta*R\n",
    "        dE  = beta*np.matmul(Nc,((I+A)/N)*S) - E/sigma - theta_E*psi_PP*E\n",
    "        dI = (1/sigma)*E - (1/omega)*I - theta_I*psi_PP*I\n",
    "        dA = (a/omega)*I - A/da - theta_A*psi_PP*A\n",
    "        dM = (m/omega)*I - M*((1-h)/dm) - M*h/dhospital - theta_M*psi_PP*M\n",
    "        dC = c*(M+MQ)*(h/dhospital) - C*(1/dc)\n",
    "        dCmirec = Mi/dmi- Cmirec*(1/dmirec)\n",
    "        dCicurec = ((1-m0)/dICU)*ICU - Cicurec*(1/dICUrec)\n",
    "        dMi = mi*(M+MQ)*(h/dhospital) - Mi/dmi\n",
    "        dICUstar = (1-c-mi)*(M+MQ)*(h/dhospital) - ICU/dICU\n",
    "        dR  = A/da + ((1-h)/dm)*M + C*(1/dc) + Cmirec*(1/dmirec) + Cicurec*(1/dICUrec) + AQ/dq + MQ*((1-h)/dm) + RQ/dq - zeta*R\n",
    "        dD  = (m0/dICU)*ICU\n",
    "        dSQ = theta_S*psi_FP*S - SQ/dq\n",
    "        dEQ = theta_E*psi_PP*E - EQ/sigma\n",
    "        dIQ = theta_I*psi_PP*I + (1/sigma)*EQ - (1/omega)*IQ\n",
    "        dAQ = theta_A*psi_PP*A + (a/omega)*IQ - AQ/dq\n",
    "        dMQ = theta_M*psi_PP*M + (m/omega)*IQ - ((1-h)/dm)*MQ - (h/dhospital)*MQ\n",
    "        dRQ = theta_R*psi_FP*R - RQ/dq\n",
    "        \n",
    "        return (dS, dE, dI, dA, dM, dC, dCmirec, dCicurec, dMi, \n",
    "                dICUstar, dR, dD, dSQ, dEQ, dIQ, dAQ, dMQ, dRQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nc_total = np.loadtxt(\"../../data/raw/Interaction_matrices/Belgium/BELtotal.txt\", dtype='f', delimiter='\\t')\n",
    "initN = np.loadtxt(\"../../data/raw/Interaction_matrices/Belgium/BELagedist_10year.txt\", dtype='f', delimiter='\\t')\n",
    "\n",
    "h = np.array([0.0205,0.0205,0.1755,0.1755,0.2115,0.2503,0.3066,0.4033,0.4770])\n",
    "icu = np.array([0,0,0.0310,0.0310,0.055,0.077,0.107,0.1685,0.1895])\n",
    "r = icu/h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ... parameters and initial conditions\n",
    "levels = initN.size\n",
    "nc = Nc_total\n",
    "parameters = {'beta': 0.0622, 'sigma': 3.2, 'omega': 2.0, 'zeta': 0, 'a': 0.43, 'm': 1-0.43, 'da': 7, 'dm': 7, 'dc': 8, 'dmi': 8, 'dICU': 8, 'dICUrec': 7, \n",
    "              'dmirec': 7, 'dhospital': 4, 'maxICU': 2000, 'totalTests': 0, 'psi_FP': 0, 'psi_PP': 1, 'dq': 14, 'h': h, \n",
    "              'c': 1-r, 'm0': np.ones(levels)*0.50, 'mi': 0.5*r, 'nc': nc}  \n",
    "\n",
    "initial_states = {'S': initN, 'E': np.ones(levels), 'I': np.zeros(levels), 'A': np.zeros(levels), 'M': np.zeros(levels), \n",
    "                  'C': np.zeros(levels), 'Cmirec': np.zeros(levels), 'Cicurec': np.zeros(levels), 'Mi': np.zeros(levels), \n",
    "                  'ICU': np.zeros(levels), 'R': np.zeros(levels), 'D': np.zeros(levels), 'SQ': np.zeros(levels), \n",
    "                  'EQ': np.zeros(levels), 'IQ': np.zeros(levels), 'AQ': np.zeros(levels), 'MQ': np.zeros(levels), \n",
    "                  'RQ': np.zeros(levels)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# -> user initiates the model\n",
    "sir_model = SEIRSAge(initial_states, parameters)\n",
    "\n",
    "# -> user runs a simulation for a defined time period\n",
    "time = [0, 200]\n",
    "t, output = sir_model.sim(time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7feebbf87a50>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEDCAYAAAAlRP8qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXxU9b3w8c939iSTScgOBAmEsIkCEhesdalFwbrbVq11u/Txeh+17XP73Jfe26fW2tV67+1t3bi0VdvetlatrUtd2ipKKyoEBDRAICxCWLJCFrLOzO/5YyYYQpZJMjNnZvJ9+5rXnOV3znw5mXw9+Z3fIsYYlFJKJT+b1QEopZSKDk3oSimVIjShK6VUitCErpRSKUITulJKpQhN6EoplSIsTegi8riI1InIhxGU/ZGIbAy/tovIkXjEqJRSyUKsbIcuIucCbcAvjTHzRnDcXcBCY8w/xCw4pZRKMpbeoRtjVgNNfbeJSKmIvCoi60XkbyIye4BDrwd+G5cglVIqSTisDmAAK4HbjTE7RORM4FHgU707RWQqMA14w6L4lFIqISVUQhcRL3A28IyI9G529yt2HfCsMSYQz9iUUirRJVRCJ1QFdMQYs2CIMtcBd8QpHqWUShoJ1WzRGNMC7BaRzwFIyPze/SIyC5gAvGNRiEoplbCsbrb4W0LJeZaI1IjIcuAGYLmIbAIqgSv6HHI98JTRISKVUuoEljZbVEopFT3D3qEP1/lHRG4Qkc3h15q+VSRKKaXiZ9g79OE6/4jI2cBWY8xhEVkG3GeMOXO4D87LyzMlJSWji1oppcap9evXNxhj8gfaN2wrF2PMahEpGWL/mj6r7wLFkQRVUlJCRUVFJEWVUkqFichHg+2L9kPR5cArQwRym4hUiEhFfX19lD9aKaXGt6gldBG5gFBCv3uwMsaYlcaYcmNMeX7+gH8xKKWUGqWodCwSkVOBnwHLjDGN0TinUkqpkRnzHbqInAQ8B9xojNk+9pCUUkqNxrB36OHOP+cDeSJSA3wTcAIYY1YA9wK5wKPh8Vf8xpjyWAWslFJqYJG0crl+mP1fAr4UtYiUUkqNSkKN5aKUUmr0Em20xWFVHWrlT5sPjPzAj4fj/XhTv10S3mITsNkEu01CyyLYRHDYhTSnnQy3g3SXnXRX6D3D7WBCupPsdNco/1VKKTV2SZfQq+vaeGhV9YiOiddwNcUT0pg/JZsFxdnMn5LNvMk+0l1Jd4mVUknKssG5ysvLTSL0FO399/deBgMEjSEQNBgDAWMIGkMwaOgJGDp7Ahzt9tPeHaC9K7Tc0R2gtqWTzfub2bTvCDWHO4DQnf7MwkzOLs3jaxfNJMOtyV0pNTYisn6whifjPsP0zozUt0bGjuC0j/6cDW1dbK45wsZ9zWzcd4Qn1+xm7Z5GHr/5dAp8njFGrJRSAxv3d+jxsGpbHXf8ZgMT0l08cevpzCzMtDokpVSSGuoOXVu5xMEFswt4+h8X0xMIcs2ja3i7usHqkJRSKUgTepzMm5zFH+74BJOy07j58bU8U7HP6pCUUilGE3ocTc5O45l/WsxZ03P5l2c3859/rkJnjFJKRYsm9DjzeZw8cevpfG5RMT95o5rXKg9ZHZJSKkVoQreA027jB9ecSqHPzTMVNVaHo5RKEZrQLWK3CVefVsyb2+upa+20OhylVArQhG6ha04rJhA0PP/+KIYyUEqpfjShW2hGgZcFU7J5dn2NPhxVSo2ZJnSLfXZRMVW1rVQeaLE6FKVUkhv3Xf+tdtmpk7j/pS08u76GeZOzrA5HjYIxhq5AF12BLroD3XQHu+kJ9Az4HggG8Bs/gWCAgAm/wsv+oJ+ACRA0wWMvYwwBE8BgCJpgaNmYY+u9f9kFTRDT+58xx8oYzLEYgWP7e5f77xvo39Z/33HLcfrLcqDYktk5k89hydQlUT+vJnSLZaU7WTK3kOc37uffLpmDy6F/NMVbT7CHxo5GGjoaaO5qpqW75YT3lq4WjvYcpcPfQbu/nQ5/x7FXp7/TkoQjCCKCDRtIaN0mtmPb+5Y5thwqeGyo6L77jnsfcLjpj7f13d93+1CxjlkUTpEoir3FMTmvJvQE8NlFxfxp80He2FbH0nlFVoeTcroD3exp2cOuI7vY37af2vZa6trrqD1aS217LQ0dDYMm5DRHGj6XD5/bh9fpxevykp+eT5oj7YSXx+HBaXPisruOf7e5cNqdOG2hl91mxyY2HOLAbrNjFzsOmwO7hLbbxY6IHFvvXe5N3r3blOpPE3oC+OSMPAoy3Ty7vkYT+hgYY9jbupctjVuoPlLNriO7qD5Szb7WfQRM4Fi5TFcmhemFFKYXUjahjMKM0HJeWh5Z7iyyXFn43D58Lh8uu05aopKHJvQE4LDbuGrhZH7+9900tHWR53VbHVJSMMawr3Uf6w6tY+2htVQcqqCuow4Am9g4KfMkSrNLuajkIkqzSinNLmVK5hTSnekWR65UbGhCTxDXLCrmv1fv4vmNB1h+zjSrw0lY/qCfdw++yyu7X+Hdg+9S1x5K4HlpeZxeeDrlReXMz5/PtKxpenetxh1N6AliZmEm84uzeHZ9jSb0AWw/vJ0Xd77IS7teoqGjAZ/Lx9mTzub0olASn+abpvXKatzThJ5ArllUzL3PV7LlQAtzJ/msDsdynf5Ofr/j9zxf/Txbm7biEAefLP4kl5dezrnF5+oduFL9aEJPIBefXMS9z1fy3u7GcZ3Qe4I9/LH6j6zYuIK6jjrm5s7lnjPuYdm0ZeR4cqwOT6mEpQk9gRRkuslKc7Kjrs3qUCwRNEFe2/MaD7//MHtb97IgfwEPnPsA5UUDzrallOpn2IQuIo8DlwJ1xph5A+wX4MfAJUA7cIsxZkO0Ax0PRISZhV521LZaHUrcbarfxHfe/Q7bmrZRNqGMhz/1MOcWn6v14kqNQCTdEp8Elg6xfxlQFn7dBjw29rDGr7LCTLbXto2bwboCwQA/3fxTbn7lZpq7mvn+J7/PM5c+w3lTztNkrtQIDXuHboxZLSIlQxS5AvilCWWgd0UkW0QmGmMORinGcaWswEtzRw/1rV0U+DxWhxNTde11/Nvf/o33Dr3H0pKl3Lv4XjJdmVaHpVTSikYd+mSg74zHNeFtmtBHYWZhKKHtqGtL6YS+umY1/+/v/4/OQCf3n30/V864Uu/IlRqjaIwENdBv4YD1BSJym4hUiEhFfX19FD469ZQVegHYnsL16M9sf4Y7X7+T/PR8nvrMU1xVdpUmc6WiIBoJvQaY0me9GBhwCh5jzEpjTLkxpjw/Pz8KH5168r1ustOdbK9NzZYuv9ryK+5/537OmXwOv77k10zPnm51SEqljGgk9BeAmyTkLKBZ689HT0SYWZCZki1dVm5eyQ/X/ZAlU5fw4wt+jMeRulVKSlkhkmaLvwXOB/JEpAb4JuAEMMasAF4m1GSxmlCzxVtjFex4MaPQy0ubDmCMSZmqiIfef4iVm1dy6fRL+fYnvo3Dpl0glIq2SFq5XD/MfgPcEbWIFDMLvLR0+lOmpcsfdvyBlZtXcnXZ1Xxz8TexiU7ioVQs6G9WAupt6ZIK9egb6zZy/7v3s3jiYr5x1jc0mSsVQ/rblYDKjiX05K5HP3T0EF9d9VUmZUziwfMe1GoWpWJMf8MSUJ7XxYR0Jzvqkjehd/o7+cqqr9AZ6OTxix8ny60TYCsVa5rQE5CIUFaQyY4krnJ5cN2DbG3cykOfekibJioVJ1rlkqDKCr1sr21NyjFd1h1ax9Pbn+bGuTdy3pTzrA5HqXFDE3qCmlmYSUunn7rWLqtDGZFOfyf3rbmPYm8xdy680+pwlBpXNKEnqGQdAuCxTY+xt3Uv3zz7m6Q50qwOR6lxRRN6gjo2SFcS1aNXNlbyi8pfcHXZ1Zw18Syrw1Fq3NGEnqByM5KrpYs/6Oe+NfcxwTOBf170z1aHo9S4pAk9QYnIsckuksGre15lW9M27j79bm2iqJRFNKEnsJlJ0tIlaIL8bPPPmJE9g4tKLrI6HKXGLU3oCWxmYSatSdDSZdXeVexs3snyU5Zr136lLKS/fQmsrCDxhwAwxrDyg5UUe4tZWjLU1LNKqVjThJ7APm66mLj16GsOrGFL4xaWn7Jcx2pRymKa0BNYntdNToYroSe7WLl5JQXpBVxeernVoSg17mlCT3BlBV6q6xLzDn197Xo21G3g1pNvxWV3WR2OUuOeJvQEN3lCGgebO60OY0A/3fxTcjw5XDPzGqtDUUqhCT3hFfk81LZ0EgwmVtPFyoZK3j7wNjfOvVG7+CuVIDShJ7hCnwd/0NB4tNvqUI7z0w9+SqYrk+tmXWd1KEqpME3oCa4wPKdobUviVLtUH67m9b2v84XZX8Dr8lodjlIqTBN6givKCiX0QwlUj/6rrb8izZHGF+d80epQlFJ9aEJPcEW9d+itiZHQuwPd/GXPX1gydQnZnmyrw1FK9aEJPcHleV3YBGoT5A797f1v09rTyrJpy6wORSnVjyb0BOew28jPdHMoQerQX9nzCtnubM6ceKbVoSil+tGEngSKfB4OtVg/QFd7Tztv7nuTJVOX4LQ5rQ5HKdWPJvQkUOjzJESVy+r9q+nwd2h1i1IJKqKELiJLRaRKRKpF5J4B9meJyIsisklEKkXk1uiHOn4V+jwJUeXy2u7XyE/L57SC06wORSk1gGETuojYgUeAZcBc4HoRmduv2B3AFmPMfOB84D9ERAf3iJKiLA/NHT109gQsi6Gtu43VNau5qOQi7Da7ZXEopQYXyR36GUC1MWaXMaYbeAq4ol8ZA2SKiABeoAnwRzXScay3c5GVbdFX7VtFd7BbxzxXKoFFktAnA/v6rNeEt/X1MDAHOAB8AHzFGBPsfyIRuU1EKkSkor6+fpQhjz+9bdGtrHZ5ZfcrTMqYxPz8+ZbFoJQaWiQJXQbY1n+kqIuBjcAkYAHwsIj4TjjImJXGmHJjTHl+fv6Igx2virLcgHXd/490HuGdA+9w8bSLCf0RppRKRJEk9BpgSp/1YkJ34n3dCjxnQqqB3cDs6ISorB7P5a97/4rf+FlWoq1blEpkkST0dUCZiEwLP+i8DnihX5m9wIUAIlIIzAJ2RTPQ8SzT4yTDZedQszVt0V/d/SolvhJm5+j/o5VKZMMmdGOMH7gTeA3YCjxtjKkUkdtF5PZwsW8DZ4vIB8DrwN3GmIZYBT0eFWZ5LLlDb+hoYF3tOpZOW6rVLUoluIhm9TXGvAy83G/bij7LB4CLohua6qsw05q26H/e82eCJqitW5RKAtpTNEkUZXksabb46p5XKZtQRml2adw/Wyk1MprQk0Shz0Nda3ynomvsaOT9uve5eOrFcftMpdToaUJPEkU+Nz0BQ1N7/KaiW1e7DoDFkxbH7TOVUqOnCT1JWDFz0bqD68hwZjA3t/9ID0qpRKQJPUlY0RZ97aG1LCpchMMW0bNzpZTFNKEnid479No4jYte117HnpY9nFF0Rlw+Tyk1dprQk0S+141N4jeey9pDawE4vej0uHyeUmrsNKEnCYfdRp7XHbeJLtYdWofP5WPWhFlx+Tyl1Nhp5WgSiedEF+8dfI/ywnId+1ylnJ6eHmpqaujstH7SmKF4PB6Ki4txOiOf7lETehIp9HmoOdwe88850HaA/W37uXHujTH/LKXiraamhszMTEpKShJ2OAtjDI2NjdTU1DBt2rSIj9MqlyRSlOWOyx16b/25PhBVqaizs5Pc3NyETeYAIkJubu6I/4rQhJ5EinwejrTHfiq6tQfXkuPJYUb2jJh+jlJWSeRk3ms0MWpCTyLxaItujGHtobWUF5YnxZdeKfUxTehJJB69Rfe17qO2vVarW5SKse9+97ucfPLJnHrqqSxYsID33ntvzOfUh6JJJB5zix6rP5+oCV2pWHnnnXd46aWX2LBhA263m4aGBrq7xz5Okyb0JFIYvkOvi2Fv0bUH15Kflk+JryRmn6HUeHfw4EHy8vJwu0PzBefl5UXlvJrQk0im20Ga0x6zO/Te+vMzJ56p9edqXPjWi5VsOdAS1XPOneTjm5edPGSZiy66iPvvv5+ZM2fy6U9/mmuvvZbzzjtvzJ+tdehJRERCE13EKKHvbt5NY2cjZ048MybnV0qFeL1e1q9fz8qVK8nPz+faa6/lySefHPN59Q49yRT6Ytf9X8dvUePNcHfSsWS32zn//PM5//zzOeWUU/jFL37BLbfcMqZz6h16kimKYff/tYfWMjFjIsXe4picXykVUlVVxY4dO46tb9y4kalTp475vHqHnmQKszzUtXRhjIlqPXfQBFl3aB3nFp+r9edKxVhbWxt33XUXR44cweFwMGPGDFauXDnm82pCTzJFPg/dgSBNR7vJ9bqjdt4dh3dwpOuI1p8rFQeLFi1izZo1UT+vVrkkmVi1Ra+orQB0/Balkpkm9CQTq7boHzZ8SEF6AUUZRVE9r1IqfjShJ5nCGN2hVzZWcnKudU/8lVJjF1FCF5GlIlIlItUics8gZc4XkY0iUikib0U3TNWrIDNUbx7NAbrautvY07xHE7pSSW7Yh6IiYgceAZYANcA6EXnBGLOlT5ls4FFgqTFmr4gUxCrg8c5pt5HndUU1oW9t2orBcHKeJnSlklkkd+hnANXGmF3GmG7gKeCKfmW+ADxnjNkLYIypi26Yqq9Cn4faKNahb2kM/b95bu7cqJ1TKRV/kST0ycC+Pus14W19zQQmiMibIrJeRG4a6EQicpuIVIhIRX19/egiVqG5RaPYW7SyoZJJGZPI8eRE7ZxKqcF5vd6YnDeShD5QLxPTb90BLAI+A1wMfENEZp5wkDErjTHlxpjy/Pz8EQerQgp9Hupao5jQGyu1ukWpFBBJQq8BpvRZLwYODFDmVWPMUWNMA7AamB+dEFV/hT43DW3ddPuDYz5Xc1cze1v3anWLUikgkp6i64AyEZkG7AeuI1Rn3tfzwMMi4gBcwJnAj6IZqPpYb+ei+rYuJmenjelcW5u2AmgLFzU+vXIPHPoguucsOgWW/SC654zQsAndGOMXkTuB1wA78LgxplJEbg/vX2GM2SoirwKbgSDwM2PMh7EMfDw71ha9uXPMCb2yoRLQB6JKpYKIxnIxxrwMvNxv24p+6w8CD0YvNDWY3oReF4Wmi5WNlUzJnEKWO2vM51Iq6Vh0Jx0r2lM0CRX6Qp2LotFbdEvjFq1uUSpFaEJPQjkZLpx2GXNb9MOdh9nftl8TulIpQhN6EhIRCjI9Y+4t2tuhSJssKhVfbW1tMTmvJvQkVZQ19oRe2Rh6IDonZ040QlJKWUwTepIq9LnHXIde2VBJia8Erys2vdaUUvGlCT1JFfo8Yx4TXXuIKpVaNKEnqUKfh7YuP21d/lEd39DRQG17rT4QVSqFaEJPUr29RUdbj37sgagmdKVShib0JFUQboteO8pRFysbKrGJjdk5s6MZllLKQprQk9SxO/RRjrpY2VjJ9KzppDvToxmWUioCdrudBQsWMG/ePC677DKOHDkSlfNqQk9SH4/nMvIHo8YYKhsrdfwWpSySlpbGxo0b+fDDD8nJyeGRRx6Jynk1oSepDLeDTLdjVHXode11NHQ0aP25Uglg8eLF7N+/PyrnimhwLpWYCnzuUSX03g5F2mRRjXcPrH2AbU3bonrO2TmzufuMuyMqGwgEeP3111m+fHlUPlvv0JPYaHuLVjZWYhc7sybMikFUSqnhdHR0sGDBAnJzc2lqamLJkiVROa/eoSexwkwP7+1uGvFxlY2VlGaX4nF4YhCVUskj0jvpaOutQ29ububSSy/lkUce4ctf/vKYz6t36EmsMCs0t2gw2H+K18EZY9jSoEPmKpUIsrKy+MlPfsK///u/09PTM+bzaUJPYoWZbnoChqb27oiPOXj0IIe7DmsLF6USxMKFC5k/fz5PPfXUmM+lVS5JrCjr46no8rzuiI7pfQA0J1dHWFTKKv2Hz33xxRejcl69Q09iBb1T0Y2gc1FVUxWCUJZdFquwlFIW0YSexIpG0bmo6nAVU31TtYeoUilIE3oSy890IzKyAbq2NW1jVo42V1TjmzGRNySwymhi1ISexJx2G7kZkXcuau1uZX/bfm1/rsY1j8dDY2NjQid1YwyNjY14PCNrWqwPRZNc4Qh6i24/vB1A79DVuFZcXExNTQ319fVWhzIkj8dDcXHxiI7RhJ7kinweDkQ4hG5VUxWA3qGrcc3pdDJt2jSrw4gJrXJJcgU+D3UR3qFXHa5ignsCBekFMY5KKWWFiBK6iCwVkSoRqRaRe4Yod7qIBETks9ELUQ2lyOeh8Wg3Xf7AsGV7H4iKSBwiU0rF27AJXUTswCPAMmAucL2InNDNMFzuAeC1aAepBlcYnrmovnXopov+oJ/qw9Va3aJUCovkDv0MoNoYs8sY0w08BVwxQLm7gN8DdVGMTw2jMCuyuUX3NO+hO9itD0SVSmGRJPTJwL4+6zXhbceIyGTgKmDFUCcSkdtEpEJEKhL9CXOy6O1cdODI0Am96nD4gagmdKVSViQJfaAK1/4NOP8LuNsYM2RFrjFmpTGm3BhTnp+fH2mMaghTc0M9Pvc0HB2yXFVTFU6bk2lZqfl0XykVWbPFGmBKn/Vi4EC/MuXAU+GHbXnAJSLiN8b8MSpRqkGluxxMzk6jur5tyHJVh6uYkT0Dp80Zp8iUUvEWyR36OqBMRKaJiAu4DnihbwFjzDRjTIkxpgR4Fvjfmszjp7TAy85hErp2+Vcq9Q2b0I0xfuBOQq1XtgJPG2MqReR2Ebk91gGq4ZXmZ7Cz7uigE100dDTQ1NmkLVyUSnER9RQ1xrwMvNxv24APQI0xt4w9LDUSpfleOnoCHGrpZFJ22gn7e8dA1zt0pVKb9hRNATMKvABU1w1c7aIJXanxQRN6CijNDyX0werRtzdtZ1LGJHwuXzzDUkrFmSb0FJDndeHzOAZN6NsO6wNRpcYDTegpQERCLV3qTmyL3uHv4KOWjzShKzUOaEJPETPyvQO2Ra8+XE3QBJk9YbYFUSml4kkTeoooLfBS39pFc0fPcdt7u/zPzJlpRVhKqTjShJ4ieh+M7up3l761cSuZzkwmeycPdJhSKoVoQk8RpfkZAOysP74efUvjFubkzsEm+qNWKtXpb3mKOCknHaddjmuL3hPsYfvh7czJmWNhZEqpeNGEniIcdhsluRnHNV3cdWQX3cFu5uaeMB+JUioFaUJPIaX5xw/StaVxCwBzcvUOXanxQBN6CiktyGBvYzs9gSAQSujpjnSm+qZaHJlSKh40oaeQGQVe/EHDR42hB6Nbm7YyO2e2PhBVapzQ3/QU0tt0sbruKP6gn6qmKq0/V2ociWj4XJUcpvcZpGtGcxudgU5N6EqNI5rQU4jX7aDI52FnfRtbmvYCaEJXahzRhJ5iZhR42VnXRkHjFtIcaZT4SqwOSSkVJ1qHnmJK8zPYWX+ULY1bmDVhFnab3eqQlFJxogk9xZQWeGnr6mZr4zatblFqnNGEnmLmTPRhczXQGejQDkVKjTOa0FPMginZeH2HAH0gqtR4owk9xTjtNqYUNUHQwdTMaVaHo5SKI03oKciedoBA10Q217RaHYpSKo40oaeYoAlS17UL0zWZv26ptTocpVQcaUJPMfta93G0p41pmbP4y1ZN6EqNJxEldBFZKiJVIlItIvcMsP8GEdkcfq0RkfnRD1VFYnP9ZgAumLaIXfVHjxtOVymV2oZN6CJiBx4BlgFzgetFpH/zid3AecaYU4FvAyujHaiKzPra9WQ6M/n8/HIAXte7dKXGjUju0M8Aqo0xu4wx3cBTwBV9Cxhj1hhjDodX3wWKoxumilRFbQWnFZ7G1JxM5kz08dctdVaHpJSKk0gS+mRgX5/1mvC2wSwHXhlLUGp06tvr+ajlI8oLQ3fnS+YUUPFRE01Huy2OTCkVD5EkdBlgmxmwoMgFhBL63YPsv01EKkSkor6+PvIoVUTW164HoLwolNA/PbeQoIFV2/QuXanxIJKEXgNM6bNeDBzoX0hETgV+BlxhjGkc6ETGmJXGmHJjTHl+fv5o4lVDqKitIMOZweyc2QDMm5RFoc/NX7UeXalxIZKEvg4oE5FpIuICrgNe6FtARE4CngNuNMZsj36YKhIVhypYULAAhy00KrLNJlw4p5C3ttdzWKtdlEp5wyZ0Y4wfuBN4DdgKPG2MqRSR20Xk9nCxe4Fc4FER2SgiFTGLWA2oqbOJnc07j9Wf97pp8VS6/UG+9/JWiyJTSsVLRBNcGGNeBl7ut21Fn+UvAV+KbmhqJI7Vn/dL6LOLfNx27nQefXMnVy2czNkz8qwITykVB9pTNEVUHKogzZHGyXknn7DvyxeWMTU3nX/7wwd09gQsiE4pFQ+a0FNERW0F8/Pn47Q5T9jncdr53lWnsKexnYffqLYgOqVUPGhCTwHNXc3sOLzjhOqWvj4xI49rTitmxVs7qTqkozAqlYo0oaeADbUbMBgWFS4KbQgOXK3y9c/MwZfm5J7nNhMIDtiVQCmVxDShJ7vuo1S891+4EE559nb4XjH84CTY/ucTiuZkuPjGpXN4f+8Rvvb0Rtq7/RYErJSKFU3oye6N71BxeCunBh24C+bCwi/ChBJ45hY4uOmE4lcumMzXlszk+U0HuPKRt3U0RqVSiCb0ZLZvLa3vrWCb2035wuXw+V/Csh/AF38PaRPg15+HI/uOO0REuOvCMn75D2fQ0NbN5Q/9nT9tPmjRP0ApFU2a0JOVvwuev5NVuZMIAmdPOvvjfZlFcMMz0NMOv/k8dDafcPgny/J56a5zmFmUyR2/2cC3Xqyk2x+MX/xKqajThJ6s3vohNFTxQvEsir3FLMhfcPz+wrlw7a+gYTs8fdOAD0onZafxu9sWc+snSnji7T18bsUa/rqllqA+MFUqKWlCT0YHN8Pff8SBU65ibXM1l8+4HJEBBsWcfj585j9h15vw7mMDnsrlsPHNy07m4S8spLaliy/9soIL/uNNHv/7blo7e2L5r1BKRZkm9GQT6IHn74D0XF6cugCD4fLSywcvf9pNMHMZvPFtaNw5aLFLT53E3+6+gIeuX0ie1839L21h8fff4L4XKtm078GdLWEAAA3PSURBVIjetSuVBMQYa35Ry8vLTUWFjuE1Yu+thFf+BfPZJ7lsx8/JT8vniaVPDH1My0F49EzInwO3vgw2+7Afs2nfEZ5cs4eXNh+gJ2DIzXBx3qx8LphVwLll+WSln9gjVSkVeyKy3hgzYC/CiAbnUgmiqxXeegCmnsOm/BI+Wv8Ry+ctH/4430RY+gD88XZYuxLO+qdhD5k/JZsfXbuAb1w6l9Xb61lVVccb2+p4bsN+7DbhlMlZLJiSzanFWZxanMX0PC8220BzoSil4kUTejJZ8xC0N8CS+3l+1wukOdK4qOSiyI6dfx1UPgd//RaUXQS5pREdlpPh4sqFk7ly4WQCQcPGfYdZta2etbubeLpiH0+u2QOA1+1g3mQfswozKS3wMiPfy/R8LwWZbk30SsWJJvRk0VoLax6GuVfSWXQyr62+i0+f9GkynBmRHS8Cl/0YHjkLXrgLbn4JbCN7hGK3CYum5rBoag4AgaChuq6NzTVH2FzTzOb9zfx+w37auj7ugepy2CjOTqM4J53J2R4KfR6KfB4Kszzke93kZLjIyXDhcQ5fDaSUGpom9GTx1gMQ6IIL7+XNfW/S2tPK5TOGeBg6EN8kWPq90EPV1T+E8+8ZU0h2mzCrKJNZRZl8rjw0S6ExhtqWLnbWt7Gr4Sg1Te3sO9zOvqYOthxopqFt4JmTvG4HEzKcZKU58XnCrzQHPo+TDLeDdJeddLeDDJeddFdo3eO043HacDvsuB02PM7Qu8thw2EXnDab/nWgxhVN6MmgoRrWPwnlt0JuKX98/0GKMoo4o+iMkZ9rwQ2w52148/uQNQUW3hDVUEWEoiwPRVkePjHAZBrd/iB1rZ3UtnTR2NZF49Fumo5209DWxeGj3bR0+mnp6GFXQxstHX6aO3roGMMY7g6bhJK73YbTbgut2wRb+N1uExw2G/bwsk1CU/fZJLws4WVbaBlC7xLeJ+F/swjhZRDC6+Flju2TYzOu97Yy7d3eu3z8Qvj4fsf0KzKogVqynlAmgjNFcp5oiNfnJILzZuazdN7EqJ9XE3oyeP1b4EyD8+6m9mgt7xx4h+XzlmOTUbQ67a16aT0AL3451Kt0xoXRj3kQLoeN4gnpFE9Ij/iYQNDQ0ROgvdtPe1eAo91+2rsDdPUE6fIH6Ay/d/mDdPYE6PYH8QdN+D1IT8DQEwjSEwgSCBr8ARN6DxoCxhAIhJeDQYIGgsZgTOhzgyZUtjtgjm03xmAA06dsMNxazBgwhMvxcVnC64S3hTfR28gsXIq+jc4Ga4AWScu0SNquRdLAzUR0prGzqLGdZSZnp8XkvJrQE92+dbD1BTj/XzEZ+Xxn1Zexi52rZlw1+nM6XPD5X8ETy0K9SG99BSaeGr2Yo8xuE7xuB163AzKtjkapxKUdixJZVxv86f9ARj4svpNndzzLm/ve5KunfZUpviljO7fHFxrvxZMFv/7cCYN4KaWSjyb0RBUMwHP/C2or4crH2N1Zz4PrHuSsiWfxxblfjM5n+CbBDc9CT0coqXccic55lVKW0ISeqP5yL1S9DEsfoGf6+dzzt3tw2V1895zvjq7ufDCFc+G6/4HGavjt9dByIHrnVkrFlSb0RFTxBLzzMJxxG5x5G49uepQtjVv41uJvUZBeEP3Pm3YuXLUC9q+Hh08PtXcP6MBcSiUbTeiJZucq+NPXYMYSuPj7rDu0jp9/8HOuLruaC6fGsDXKKZ+FO96FqWfDn78O/30efPRO7D5PKRV1mtAThb8rNPDW0zdB/iw6r3yMhzY/xj/+5R+ZkjmFu0+/O/Yx5EyHLzwN1/5PaFKMJ5bCH/4J6qvGX7sypZKQNlu0WsAPm38Hb/4AmvfC1HP429nL+d5rN1HTVsNl0y/ja+VfI90ZebvtMRGBOZdB6adCk2i88zBs+k1ontKyi2HmxVByDjjc8YlHKRWxiIbPFZGlwI8BO/AzY8wP+u2X8P5LgHbgFmPMhqHOOa6HzzUGjtbD7tXw1gOYhu3sn3Qq759yKX/p2M+qfaso8ZXwjbO+wRkTR9EbNJqa98P2V2H7a7D7LfB3gjMDTjoLcmeE7upzS0Pv2SeBXYfVVSqWhho+d9iELiJ2YDuwBKgB1gHXG2O29ClzCXAXoYR+JvBjY8yZQ503ZRN6MBiay7OnA3qOQk8HwaP1tNZupqWukubGHTQ3f0RzTxv1djubfXlsTM+grqcFgExnJrfMu4VbTr4Fl91l8T+mn+522PO3UIKvqYCm3dDd+vF+sUN6buiVkReaqNrtA3dm6OXKAGc6OD2hd7srdKd/7N0NdgfYnGBzHL9ss4fexXb8stjD27T2UI0PYx0P/Qyg2hizK3yyp4ArgC19ylwB/NKE/u/wrohki8hEY0zUp5N/e93DPPjByhEdM5LaXxPBmiHcvfvYeqirtx9DAAgI+BH8AgGEHgHTO1CFDZjgATwATMwooLxgIQvDrxnZM7BHMAGFJVzpoSqXmReH1o2Bow3QtDM0G9Lh3dBWB+2NoVd9FXS3hcZx72plZD+JUeib4MUWHkzFRmgwld51Ca9H8k6/bcd9WPhNBlnvV+7YakSjsERQJErnSSTjaTCXhTfC2XdG/bSRJPTJQN9uhDWE7sKHKzMZOC6hi8htwG0AJ5100khjBSAjLYfpzqwRHzeS78rAAxaFB1ASCQ+4JH0GVrJhE8Fuc+Kwu7DbnThsLux2Fw67G6cznazsErKyTiLLnU2WOwuf28cE9wQmeCaM+N+SMETAmx96nXTW0GWDQfB3QE/nx3/B+Dsh0B16IBzoDr96INgT6lh1bNkfWjfBj5eD/tC6MWB694XfTSC8PRj6bBPs8wr973jYd3rf+m2DPg+IzcDrJ5Rj4P0Diejhc7TOk0iSLd4x8sag+TGRJfSBslv/qx9JGYwxK4GVEKpyieCzT7Bg3hdYMO8LozlUWclmC1W5uDKAXKujUSolRVLxWAP0HTikGOjfnTCSMkoppWIokoS+DigTkWki4gKuA17oV+YF4CYJOQtojkX9uVJKqcENW+VijPGLyJ3Aa4SaLT5ujKkUkdvD+1cALxNq4VJNqNnirbELWSml1EAi6lhkjHmZUNLuu21Fn2UD3BHd0JRSSo2ENt5VSqkUoQldKaVShCZ0pZRKEZrQlVIqRUQ0OFdMPlikHvholIfnAQ1RDCdaEjUuSNzYNK6R0bhGJhXjmmqMyR9oh2UJfSxEpGKwwWmslKhxQeLGpnGNjMY1MuMtLq1yUUqpFKEJXSmlUkSyJvSRjZ8bP4kaFyRubBrXyGhcIzOu4krKOnSllFInStY7dKWUUv1oQldKqRSRdAldRJaKSJWIVIvIPRbGMUVEVonIVhGpFJGvhLffJyL7RWRj+HWJBbHtEZEPwp9fEd6WIyJ/EZEd4fe4TpUkIrP6XJONItIiIl+14nqJyOMiUiciH/bZNuj1EZF/DX/fqkTk4jjH9aCIbBORzSLyBxHJDm8vEZGOPtdtxeBnjklcg/7cLL5ev+sT0x4R2RjeHs/rNVhuiP13zBiTNC9Cw/fuBKYDLmATMNeiWCYCp4WXMwlNpD0XuA/4vxZfpz1AXr9tPwTuCS/fAzxg8c/xEDDViusFnAucBnw43PUJ/0w3AW5gWvj7Z49jXBcBjvDyA33iKulbzoLrNeDPzerr1W//fwD3WnC9BssNMf+OJdsd+rEJq40x3UDvhNVxZ4w5aIzZEF5uBbYSmkc1UV0B/CK8/AvgSgtjuRDYaYwZbU/hMTHGrAaa+m0e7PpcATxljOkyxuwmNOb/GfGKyxjzZ2OMP7z6LqHZwOJqkOs1GEuvVy8JTfb7eeC3sfjsoQyRG2L+HUu2hD7YZNSWEpESYCHwXnjTneE/kR+Pd9VGmAH+LCLrwxNzAxSa8CxS4ffYzFIbmes4/hfN6usFg1+fRPrO/QPwSp/1aSLyvoi8JSKftCCegX5uiXK9PgnUGmN29NkW9+vVLzfE/DuWbAk9osmo40lEvMDvga8aY1qAx4BSYAFwkNCfffH2CWPMacAy4A4ROdeCGAYkoWkMLweeCW9KhOs1lIT4zonI1wE/8OvwpoPAScaYhcA/A78REV8cQxrs55YQ1wu4nuNvGuJ+vQbIDYMWHWDbqK5ZsiX0hJqMWkSchH5gvzbGPAdgjKk1xgSMMUHgp8Toz82hGGMOhN/rgD+EY6gVkYnhuCcCdfGOK2wZsMEYUxuO0fLrFTbY9bH8OyciNwOXAjeYcKVr+M/zxvDyekL1rjPjFdMQP7dEuF4O4Grgd73b4n29BsoNxOE7lmwJPZIJq+MiXEf3c2CrMeY/+2yf2KfYVcCH/Y+NcVwZIpLZu0zoodqHhK7TzeFiNwPPxzOuPo67c7L6evUx2PV5AbhORNwiMg0oA9bGKygRWQrcDVxujGnvsz1fROzh5enhuHbFMa7Bfm6WXq+wTwPbjDE1vRvieb0Gyw3E4zsWj6e+UX6CfAmhp8Y7ga9bGMc5hP4s2gxsDL8uAX4FfBDe/gIwMc5xTSf0xHwTUNl7jYBc4HVgR/g9x4Jrlg40All9tsX9ehH6H8pBoIfQ3dHyoa4P8PXw960KWBbnuKoJ1a/2fsdWhMteE/75bgI2AJfFOa5Bf25WXq/w9ieB2/uVjef1Giw3xPw7pl3/lVIqRSRblYtSSqlBaEJXSqkUoQldKaVShCZ0pZRKEZrQlVIqRWhCV0qpFKEJXSmlUsT/B7R73BfDBNPEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(t, output[\"S\"].sum(axis=0), label=\"S\")\n",
    "plt.plot(t, output[\"I\"].sum(axis=0), label=\"I\")\n",
    "plt.plot(t, output[\"R\"].sum(axis=0), label=\"R\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:COVID_MODEL]",
   "language": "python",
   "name": "conda-env-COVID_MODEL-py"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}

refactor_tryout.ipynb

Beautiful!