experiments/single_population/logistic_growth.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Logistic Growth\n",
    "\n",
    "Populations cannot grow forever without restrictions, like exponential growth. A population would more realistically grow logistically: exponential at small numbers but slows down as the population nears a carrying capacity. This example builds several logistic growth models, outlining how more advanced models can be created."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differential Equation\n",
    "\n",
    "Logistic growth model can be written as a differential equation,\n",
    "\n",
    "\\begin{equation*}\n",
    "\\frac{dN}{dt} = rN(t) \\left( 1-\\frac{N(t)}{k} \\right)\n",
    "\\end{equation*}\n",
    "\n",
    "where a population, $N$, increases by a growth rate, $r$, that is reduced as the population reaches its carrying capacity ($k$). There is a simple solution to this differential equation,\n",
    "\n",
    "\\begin{equation*}\n",
    "\\ N(t) = \\frac{kN_0 e^{rt}}{k+N_0(e^{rt}-1)}\n",
    "\\end{equation*}\n",
    "\n",
    "based on an initial population size, $N_0$. \n",
    "\n",
    "I'm going to make one change. The growth rate $r$ is the difference between  birth and death rates. ($r = b - d$) and both are density-dependent in the basic logistic model. Techincally, as the population approaches their carrying capacity, death rates decrease &mdash; that seems strange. Instead, the death rate is modelled as density-independent intrinsic factor.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\frac{dN}{dt} = bN(t) \\left( 1-\\frac{N(t)}{k} \\right) - dN(t)\n",
    "\\end{equation*}\n",
    "\n",
    "Here's how to plot the solution in Python:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgIAAAEMCAYAAAC2p94hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd5hb5Zn38e89mtGMpnqaZ9wbtsc23eCYDjEJEJaFhUBgdwPZNWyWdZaU5aUkedMISUhY4E0hjUBgkywhCUkgYCD05phiqrtxt6fPWNOkUXveP6SxhJnxHM+oP/fnunRJOkc6un9+bOvWqWKMQSmllFJ2Ksh0AUoppZTKHG0ElFJKKYtpI6CUUkpZTBsBpZRSymLaCCillFIW00ZAKaWUslhhpgtIp2effdYUFxdnugyllMopAwMDHY2Njf9dWFi4AP0BmWsioVBovd/v/+7ixYsjw73AqkaguLiYpqamMb13x44dzJgxI8kVZTfNbAfNbIfxZN60adNgbW3tJcXFxcN+kajsNjg4eERnZyfAd4abr52dQ0VFRZkuIe00sx00sx3Gk9ntdpdoE5C7iouLI7G1OcPSRsChqqqqTJeQdprZDprZDuPJLCJJrERlyIjf99oIONTR0ZHpEtJOM9tBM9sh1zNPnjz5mNFe893vfndiX19fyr/XNm/e7L7nnntqDvV9l19++cxf//rX1QdOj0QifPWrX5105JFHHn7UUUcdfuaZZ8574403SobmL1y48Ijjjjtu4XHHHbfw2GOPXXT99ddP9vl8MlRLY2PjsUuWLFk4dPv5z39eeyh1aSPgkP6CsINmtoNmzk+//OUvG/r7+5PyvRYMBkect3Xr1uIHH3zwkBuBkdxxxx31r7/+etnLL7+87q233nr3s5/9bMsnP/nJwwYGBvavinnkkUc2vfbaa+ueffbZ9Tt37iz+9Kc/vX+HjylTpgy+8sor64ZuV111VeehfL42Ag4FAoFMl5B2mtkOmtkO+ZL58ccfr1i2bNn8j3/847OPPvroRf/4j/84KxKJcOutt07s6OgoOvfcc+edeeaZ8wAefvjhylNPPbVp6dKlCy6++OLZPT09BQB//OMfq44++uhFp59++vwVK1ZMO++88w4D+NKXvjT5yiuvnHH22WfPvfzyy2dt3rzZfcYZZ8xfunTpgqVLly549tlnywC++c1vTnnzzTfLlyxZsvCWW26ZGAqF+PznPz/1xBNPXHD88ccv/MEPflAH0V/6V1999fRjjz120XnnnXdYV1fXsDvo//SnP51066237iovL48AnHfeeT1HH310/y9/+csP/LKvrKyM3HnnnTueeeaZCe3t7a5k/JladdTAePh8vkyXkHaa2Q6a2Q7JyvzRu95YnJQFHeCJK4953elrN23a5HnhhRe2Tps2LXjGGWc0PfPMM+XXXntt2913393wyCOPbGpoaAi1trYW3n777ZMeeeSRTRUVFZGbbrqp8Xvf+17DjTfe2HLDDTfMePjhhzfMmzcvcOmll85KXPbatWtLn3zyyQ1lZWWmr6+v4JFHHtlUWlpq1q1bV3zllVfOfvnll9d/+ctf3vPDH/6w4eGHH94C8IMf/KCusrIy/PLLL6/3+XyybNmyprPOOqvn9ddfL922bVvxK6+8snbv3r1FJ5988qLLLrvsfb/Wu7u7C/x+f0FTU9Ng4vSjjjqqf8OGDSUMo7q6OjJp0qTAhg0bShobG4N79uwpXrJkycKh+d/61rd2nnnmmX1O/zy1EXCosbEx0yWknWa2g2bODcYYIgZCEfO+W3joZoYeQyj2OBIxhA1EjMFVmD+bBhYtWtQ/c+bMIEBTU9PA9u3b3Qe+5sUXXyzbtm1bybJly5oAQqGQHHXUUX3vvvtuydSpUwfnzZsXALjooou67rvvvvqh9y1btmxfWVmZAQgEAnLNNdfM2LBhg6egoIBdu3YNeyKa5557rnLz5s2lK1eurAbo7+93bdy4seTll1+uOP/887sKCwuZPn168Ljjjut1mtEY43gnzaFNA06XfSBtBBxqaWmx7rhjzWwHzTw2xhj8oQi+YARfMBy9D0Uf+4MR/KHobTDxFjYMhiIEwhECYUMgFL0PhuP3wYghGDYEIxGC4YQv/bDBjKPe2VUufnLxhHFlhkP75Z4qRUVF+/8oXC4XoVDoA9+Yxhg+9KEP9dx///3bEqevXr3ac7Bll5aW7j9M8r//+78b6urqgqtXr94WiUSYNGnSsGtDjDFy00037fz7v//7nsTpTzzxRNVoX+bV1dWRkpKSyMaNG93z58/fv/3mnXfeKT3hhBOG/VXv9XoLmpub3QsWLPB3d3ePe/OANgIOud0faDjznma2g+2Zg+EIPf4w+/xBevxhvP4QPYMhegfD9Mbu+wJh+gbD9AdC9AciDATD9AfCRMbzzTwGBQKFBUJhgeAaupf446H5rqGbRKe5CoQ6d/6fBqC0tDTs9XoLGhoaOPnkk/u//OUvT1+/fn3xggULBvv6+gp27NhRdPjhh/t3795dvHnzZvfcuXMDf/zjH0fc6a+np8c1efLkgMvl4uc//3ltJBL9M6ysrAwPDAzs/wI+/fTTvffcc0/92Wef3et2u83atWuLp0+fHjzxxBN777vvvvorr7yys7m5uei1116ruOiii7oO/Jyrrrqq5dprr51+//33v1dWVmYeeeSRijVr1lT85Cc/2TFMTQUrVqyYccYZZ+yrq6sLayOQRhUVFZkuIe00sx3yNXMoYujoD9DeH6SjP0BHf5DOgeito3cQb6Cbfb4QfYHwmD/D7RI8RS5KiwrwFLnwFBVQUliw/76k0IW7UCgpLKA4dnO7CiguFNyugthNKNp/H31c5BLcBQW4Cog+T/hyH6u+PsebjHPWZZdd1nHJJZfMraurCz755JObbr/99u3Lly+fHQwGBeC6667bs2jRosGbb755x8UXXzx3woQJoSOPPLJ/pOVdffXVbVdcccWcRx99tHrp0qW9JSUlEYDFixf7XC6XWbJkycKLLrqo49prr23btWtX8YknnrjAGCPV1dXBBx544L1PfOIT+55//vnKJUuWLJoxY4Z/8eLFw24a+MIXvtC2b98+19KlSxcVFBRQV1cXvO+++7YMbaIAOPfcc+cZY8QYw5lnnrnva1/72t6heQfuI3DJJZd0XHvttW1O/9zEmDS3tBm0atUqo6cYdk4z2yFXMxtj2OcLsadnkL09gzT3BmjpHaSlN0Brb4DOgaCjVekFAlUlhftvlSWFVBa7qCwupKLYRXlxIeXFLsrd0VtZsYuyIhelbheF4/hiTrfxjPOOHTs6Z8yYsT25FWVOT09PQWVlZWRor/7Zs2f7r7/+esdfnLmora1t/fz58z853DxdI+BQdfUHzgGR9zSzHbI9szGGtr4g27t97Oj2s32fn137/Oz2DtJ/kF/zBQK1niLqy4uoK3NTV1pEbVkRNZ4iSggwta6Kak8RFcUuCiw4c162j3M6/fjHP677wx/+UBcKhaSpqWlgxYoVuX22pXHSRsAhn89HZWVlpstIK81sh2zKHI4Ydnv9bOoYYHOHj/c6fWzt8o34hV/mdjGlspgpVcVMqnAzqbKYxnI3DRVu6srcI/5ib21tpaH6oPuM5Z1sGudMu/7669vyfQ3AodBGwCG/35/pEtJOM9shk5n7BkOsa+tnbWs/69v62dg+gC/4wZ3aqkoKmVVTwowJHmZUlzB9QgnTqoqZ4Ckc03nwdZyVitNGwKFcPO54vDSzHdKZeSAQ5u2WPt7Y28s7zX281+n7wHb8ieVFzKsrZW5dKYfVljK71kPNGL/wR6LjfGhs2pcsj4142Ig2Ag7psdZ20MzJZYxhW5ef1bu8vLqrh/Vt/YQTvlMKC4T59aUsaihjwcToraY09ZcI1nE+NIFAwD84OFiglyLOTYODgwWhUGj9SPO1EXCopGTYMz3mNc1sh2RnDkcM77b08eL2fby8w0t7f/ziLQUCCyaWcszkCo6aXMGCiWWUFKb/kic6zoemr6+vpbOz8/exa9rrNWpySyQUCq33+/3fHekF2gg45PHYtWMRaGZbJCNzxBjWtvbzzJZuXti+D68/tH9ejaeQ46dV8qFpVRwzpYIyd1KukzIuOs6HbtGiRd9JUikqy2gj4FB3d7d1e9xqZjuMJ/PenkEe39jJU+910dYX/+U/ubKYU2ZWcdLMCcyrL826w/N0nJWK00bAodraD1wNMu9pZjscauZgOMIL2/bx6IZO3m6Jn61uYnkRZ8yu5vQ51cyu8SR1575k03FWKk4bAYd6e3spLy/PdBlppZnt4DRzZ3+Qv2zo4NENHXT7oqv+iwsLOHXWBM6aV8PhjeVZ98t/JDrOSsVpI+BQIBAY/UV5RjPbYbTM27t9/P7tNp5+r5tQ7Co7s6pLOG9hPWfMqc6Kbf6HSsdZqThtBBzS447toJnjNnUM8Os1Laza6QWie/yfMmsC5y+s54jGsqxe9T8aHWel4rQRcEiPO7aDZo42AL9a08zfdkYvrV7sEs6aX8uFh09kcmVxpspMKh1npeK0EXBIDzeyg82Z9/YMcs9re3lu6z4guv3/7xfU8fEjJ1LtSf1JftLJ5nFW6kDaCDjkdrszXULaaWY7BHBx56rd/GV9B6GIocglnL+wnovzsAEYYuM425hZOaONgENer5cJEyZkuoy00sz5LRwxPL6pk7tW76YvaBDgI3NruGLxJCaW5/eXhk3jPMTGzMoZbQQcqqury3QJaaeZ89eWjgHueHEXmzoGADiisZz/OGEKc2pLM1xZetgyzolszKycyeg5o0Vkmog8IyLrRWStiHw2Nr1GRP4qIptj99UJ77lRRLaIyEYROSth+mIReSc27/uS5F2avV5vMheXEzRz/hkMRbjrlT185s8b2dQxQF1pEf9+dBW3nnuYNU0A5P84D8fGzMqZTF88IgT8lzFmAbAUWCEiC4EbgKeMMXOBp2LPic27FFgEnA3cKSJDBzH/GPg3YG7sdnYyCw0Gg6O/KM9o5vyytqWPTz+4gQfebsMYuGBRPXd9fAGL6105fSjgWOTzOI/ExszKmYxuGjDGNAPNsce9IrIemAKcD5wee9m9wLPA9bHp9xtjBoFtIrIFWCIi24FKY8wqABG5D7gAWJmsWm08Blcz54dAOML/rGnhd2+3EjEws7qEz58ynQUTy4D8zDwazaxUXKbXCOwnIjOBY4DVQEOsSRhqFibGXjYF2JXwtt2xaVNijw+cnjQtLS3JXFxO0My5b2e3n2v+vInfvtUKwGVHNfCjC+bvbwIg/zI7oZmVisuKnQVFpBz4A/A5Y0zPQVZTDjfDHGT6+7S1tbF8+XIKCwsJh8NceOGFrFixgpaWFsrKynC5XPT09FBfX09XVxfGGOrr62ltbSUSidDZ2UlfXx8NDQ20t7cjItTU1NDe3k5lZSXhcJj+/n4aGxtpaWmhqKiIqqoqOjo6qKqqIhAI4PP59s93u91UVFTQ2dlJdXU1Pp8Pv9+/f35JSQkej4fu7m5qa2vp7e0lEAjsn+/xeHC73Xi9Xurq6vB6vQSDwf3zR8s0dN7xkTIZY+jo6MirTKON09Dfk3zItKZL+MkrLQQiMLHUxb8sKuXkplr27t71vkwiQmtra05kStbfPZfLRXNzc15lGm2cXC4Xe/fuHVMmld/EmA98X6a3AJEi4C/A48aY22LTNgKnG2OaRWQS8KwxZr6I3AhgjPl27HWPA18DtgPPGGOaYtMvi73/04mftWrVKtPU1DSmOru7u6murh79hXlEM+cmfyjCD1/axRObuwD48JxqrjlpGqUjXBMgHzIfKs18aNasWfP6smXLjktySSpLZPqoAQF+AawfagJiHgKuiD2+AvhzwvRLRaRYRGYR3Snwldjmg14RWRpb5uUJ70mKnp6eZC4uJ2jm3NPSO8jnH97EE5u7KHYJXzhlOtefPmPEJgByP/NYaGal4jK9aeAk4JPAOyLyZmzaF4HvAA+IyHJgJ3AxgDFmrYg8AKwjesTBCmNMOPa+q4FfAh6iOwkmbUdBgPr6+mQuLido5tzyxp5ebn56Gz2DYSZXFvPVM2cxq2b008rmcuax0sxKxWX6qIEXGX77PsCyEd5zM3DzMNNfAw5PXnXv19XVRWmpPcdZg2bOJQ+ta+fOVbuJGDh+aiU3njGD8mJn/7xzNfN4aGal4jK9RiBnZHpfikzQzNkvHDH8/JU9PPhuOwCfOKqBTy2ehKvA+XkBci1zMmhmpeK0EXDIxtVqmjm7+UMRbnlmOy/t8FJYIHzu5Gl8dF7tIS8nlzIni2ZWKi5rziOQ7VpbWzNdQtpp5uzVOxjihke38NIOL+VuF986e86YmgDInczJpJmVitM1Ag7ZeDytZs5OnQNBbly5he3dfurLivj22YcxvbpkzMvLhczJppmVitNGQKkcsrdnkBtWbqGlN8C0qmK+fc5heX/JYKVUaummAYf6+voyXULaaebsssfr59q/bKalN8D8+lJuO29eUpqAbM6cKppZqThdI+BQQ0NDpktIO82cPXbt8/N/Ht1M10CIIxrLuemjsw96kqBDka2ZU0kzKxWnawQcam9vz3QJaaeZs8PObj/XPhJtAo6aVM43z0peEwDZmTnVNLNScbpGwCHbrtcOmjkb7O0Z5LqVm+n2hTh6cjnf+OgcSgqT279nW+Z00MxKxekaAYdqamoyXULaaebMausLcP2jW/avCUhFEwDZlTldNLNScdoIOGTjajXNnDndA0FuWLmF1r4ATfWlfP0js1PSBED2ZE4nzaxUnDYCDlVWVma6hLTTzJnRHwhz42Pvsds7yOwaDzefPSep+wQcKBsyp5tmVipOGwGHwuHw6C/KM5o5/QLhCF9/citbu3xMrSrm2+fMocLhxYPGKtOZM0EzKxWnjYBD/f39mS4h7TRzekWM4XvP7eDNvX3UeAr51tlzqPYUpfxzdZztYGNm5Yw2Ag41NjZmuoS008zp9bPVe3hu6z5Kiwq4+ew5NFYUp+VzdZztYGNm5Yw2Ag61tLRkuoS008zp8/C6dh58t53CAuGrZ85mTm36rhuv42wHGzMrZ7QRcKioKPWraLONZk6P13b38KNVuwH43MnTOGZKRVo/X8fZDjZmVs5oI+BQVVVVpktIO82ceju6fXzzqW1EDFx6VMOYLyU8HjrOdrAxs3JGGwGHOjo6Ml1C2mnm1Orxh/jKE1sZCEY4eeYEPnXcpLR9diIdZzvYmFk5o42AQzZ205o5dcIRw7ee2U5zb4DDaj1cd/oMCjJ0ClgdZzvYmFk5o42AQ4FAINMlpJ1mTp27X93Lmj29VJUU8rUUnjXQCR1nO9iYWTmjjYBDPp8v0yWknWZOjWfe6+J377ThEvi/y2Yxsdyd8s88GB1nO9iYWTmjjYBDNh6Dq5mTb1uXj9ue3wnA1SdM5chJ5Sn9PCd0nO1gY2bljDYCDtl4DK5mTq7+QJhvPLmNwbDho3NrOG9BXco+61DoONvBxszKGW0EHHK7M7v6NhM0c/IYY7j9hZ3s6RlkVnUJnzlpWtZcH17H2Q42ZlbOaCPgUEVFek/ykg00c/L8aW07z2+Lnj74/545K6M7Bx5Ix9kONmZWzjj630iirhKRp0Xk7di0U0XkktSWlz06OzszXULaaebk2Njez89W7wHgC6dOZ2pVSdI/Yzx0nO1gY2bljNOfJd8AlgM/A6bHpu0Grk9FUdmouro60yWknWYev/5AmG89vZ2wgfMX1nPqrOz7M9VxtoONmZUzThuBTwF/Z4y5HzCxaduA2akoKhvZeOiNZh4fYwzff2nX/pMGXfWhyUlbdjLpONvBxszKGaeNgAvoiz0eagTKE6blPb/fn+kS0k4zj88Tm7t45r1uSgoL+OKHZ+J2Zc9+AYl0nO1gY2bljNP/mR4FbhORYojuMwDcBDycqsKyjY3H4Grmsdu1z88PX45eUfAzJ07Nuv0CEuk428HGzMoZp43AF4DJgBeoIromYAZwQ4rqyjo2HoOrmccmFDHc8uwOBkMRzphTzUfm1iShstTRcbaDjZmVM4VOXmSM6QEuEJGJRBuAXcYYq/5WlZRk7y+6VNHMY/OrNc1s6hhgYnkR12TR+QJGouNsBxszK2ecHj54t4icZYxpM8a8OtQEiMidqS0ve3g8nkyXkHaa+dCtbenj/rdaEeC602ZS5nYlp7AU0nG2g42ZlTNONw38M/ALEbl2mOlW6O7uznQJaaeZD81AIMwtz+0gYuCSIydmxXUEnNBxtoONmZUzThsBP7AUuFRE/kdEhs5Vmd3rPJOotrY20yWknWY+ND9dvYeW3gBzaj1cvnhSEqtKLR1nO9iYWTnj+HgmY8xu4BSihxK+JCJTiB9KOCaxTQ5tIvJuwrSvicgeEXkzdvtYwrwbRWSLiGwUkbMSpi8WkXdi874vKdgo29vbm+xFZj3N7Nwru7ys3NhJUYFw3WkzKMrSQwWHo+NsBxszK2ec/m8lAMYYnzHmH4E/AK8AxeP8/F8CZw8z/XZjzNGx26MAIrIQuBRYFHvPnSIytAH2x8C/AXNjt+GWOS6BQCDZi8x6mtmZ3sEQt7+wC4ArFk9iVk1ubYvVcbaDjZmVM4dyiuH9jDHfAa4EfjOeDzfGPA90OXz5+cD9xphBY8w2YAuwREQmAZXGmFXGGAPcB1wwnrqGY+MxuJrZmTtX7aZzIMjCiWVcdMTEFFSVWjrOdrAxs3LG6eGD3xtm2kpgZdIrivqMiFwOvAb8lzGmG5gC/C3hNbtj04KxxwdO/4C2tjaWL19OYWEh4XCYCy+8kBUrVtDS0kJZWRkul4uenh7q6+vp6urCGEN9fT2tra34fD7q6+vp6+ujoaGB9vZ2RISamhra29uprKwkHA7T399PY2MjLS0tFBUVUVVVRUdHB1VVVQQCAXw+3/75brebiooKOjs7qa6uxufz4ff7988vKSnB4/HQ3d1NbW0tvb29BAKB/fM9Hg9utxuv10tdXR1er5dgMLh//miZysujO7ONlMnn81FXV5dXmUYbJ7/fv79uJ5kee3snT23pw+0SLjusEN9Af9ZlGm2cBgcHqaqqyqlxGu/fvUAgQHl5eV5lGm2cAoEAZWVlY8qk8ptEf0QPM0PkMWPM2bHHLzDC/gDGmFPHVYDITOAvxpjDY88bgI7Y590ETDLG/KuI/AhYZYz5Vex1vyB6xsOdwLeNMWfGpp8CXGeMOe/Az1q1apVpamoaU51tbW1MnJh7v/bGQzMfXO9giKt+v54uX4irl07hHw7PzT8rHWc7jCfzmjVrXl+2bNlxSS5JZYmDrRG4L+HxXakuZIgxpnXosYj8HPhL7OluYFrCS6cCe2PTpw4zPancbvfoL8ozmvngfvq3PXT5QixqKOP8RfUprCq1dJztYGNm5cyIjYAx5jcJj+9NTzkgIpOMMc2xp/8ADB1R8BDwGxG5jejpjucCrxhjwiLSKyJLgdXA5cAPkl2X1+tlwoQJyV5sVtPMI3t1Vw9PbO7C7RK+cMp0CrL87IEHo+NsBxszK2cOuo+AiCwGBo0x78ae1wN3AIcDq4BrjTFjvgKhiPwvcDpQJyK7ga8Cp4vI0UQ3DWwHPg1gjFkrIg8A64AQsMIYE44t6mqiRyB4iO63kPR9F+rq6pK9yKynmYfXHwhzx4s7Abh88SSmTcjtU7fqONvBxszKmdF2FrwD+DrxX+V3Ef01/jPgMuC7wH+M9cONMZcNM/kXB3n9zcDNw0x/jWhzkjJer5eysrJUfkTW0czDu/vVvbT3B5lfX8pFObpfQCIdZzvYmFk5M1ojsAB4AUBEJgDnAIcbYzaJyEPAy4yjEcglwWAw0yWknWb+oLUtfTy8vgOXwBdOmY6rIHc3CQzRcbaDjZmVM6OdR6AQGDoLxVKgxRizCcAYswuwZoOTjcfgaub3C4Qj3PZCdJPAJUc15NyJg0ai42wHGzMrZ0ZrBNYCF8ceXwo8OTQjdophb4rqyjo2XstbM7/f/W+2sss7yNSqYv7p6Pz5T1XH2Q42ZlbOjLZp4HrgYRH5CRAGTk6Y9wngpVQVlm1s3LammeO2d/u4/63oka2fO3k67sLcuZbAaHSc7WBjZuXMQRsBY8yLIjIdmAdsMsYkXrXiEeD+VBaXTVyu7L+ufLJp5qiIMfy/F3cRihg+1lSbM5cXdkrH2Q42ZlbOjPqzxhjTa4x5/YAmAGPMRmNM0k/ck616enoyXULaaeaoxzZ2sra1n2pPIVcePzkDVaWWjrMdbMysnMmf9ZspVl+fu2eOGyvNDN0DQe56JdrvXr10KuXFji7PkVN0nO1gY2bljDYCDnV1Ob1IYv7QzPCT1XvoC4Q5bmoFp83Oz4NkdJztYGNm5Yw2Ag6NdHGmfGZ75td29/DMe90Uu4T/PGkaksOnET4Y28fZFjZmVs4c0npOEZkIvG9PKWPM1qRWlKVsXK1mc+bBUIQfvrwLgH86tpFJFcWZLCulbB5nm9iYWTnjaI2AiJwtInuAZmBLwm1zCmvLKq2traO/KM/YnPm3b7WytyfAjOoSPn5EQ4arSi2bx9kmNmZWzjjdNPAj4Cag3BhTkHCz5niU8vL8OmTMCVsz7/b6+W3snAHXnDSNwjw4jfDB2DrOtrExs3LG6aaBauCnRjcyqTxnjOEHL+0mGDGcNa+GIxr1P0+lVH5zukbgF8C/pLKQbNfXN+arLecsGzM/vaWLN/b2UlHs4solUzJdTlrYOM6aWak4p2sElgLXiMgNwPtOWG2MOTXpVWWhhob83k48HNsy9wfC/GGLH4Arj59MVUn+nTNgOLaNM2hmpRI5/Z/urtjNWu3t7UybNi3TZaSVbZnve72Zbn+YhRPLOGt+babLSRvbxhk0s1KJHDUCxph7U11ItsvXY8gPxqbM73UO8Od17QjwnydNpcCi7DaN8xDNrFSc4xMKici/iMjTIrIxdm/VPgM1NTWZLiHtbMkcie0gGDFw7rwJzKktzXRJaWXLOCfSzErFOT2PwJeAG4hebfCa2P11selWaG9vz3QJaWdL5r9u7mJdWz81nkI+kn/XFBqVLeOcSDMrFed0H4ErgdONMTuGJojI48DzwBFPsGgAACAASURBVM2pKCzbVFZWZrqEtLMhc48/tP+iQld9aAqNtfatPrVhnA+kmZWKc7ppoAw4sJ3sBDzJLSd7hcPhTJeQdjZk/uXrzXj9IY5sLOfDc6qtyHwgzWwHGzMrZ5w2Ao8BvxaR+SLiEZEm4F7g8dSVll36+/szXULa5XvmTe0DPLK+A5fAZ06aiojkfebhaGY72JhZOeO0EfgM0Au8BfQBbwL9wH+mqK6s09jYmOkS0i6fM0eM4Qcv78IA/3D4RGZWR1du5XPmkWhmO9iYWTnjqBEwxvQYYy4HSoFJQKkx5nJjzL6UVpdFWlpaRn9RnsnnzCs3drKxfYDa0iL++Zj4f5D5nHkkmtkONmZWzoy4s6CIzDTGbI89nn3A7PKhY1JtuQxxUVFRpktIu3zN3OMPcfer0R0EP/2hKZS649fOytfMB6OZ7WBjZuXMwY4aeAeoiD3eAhjgwF2qDWDFFQirqqoyXULa5Wvmu1/bS+9gmGMml3Pa7Anvm5evmQ9GM9vBxszKmRE3DRhjKhIeFxhjXAdcgtiqyxB3dHRkuoS0y8fMG9r6WbmhE5fAihOmfeBsa/mYeTSa2Q42ZlbOOD2h0PdHmH5HcsvJXjZ20/mWORwx/PDl3RjgoiMmMr265AOvybfMTmhmO9iYWTnj9KiBT40w/ZNJqiPrBQKBTJeQdvmWeeXGTjZ1DFBXWsQ/HTP8HtT5ltkJzWwHGzMrZw56ZkER+deh1yU8HjIbsGZdk8/ny3QJaZdPmb3+EPe8Ft1B8N+XTsFTNPxWrXzK7JRmtoONmZUzo51ieOgXv5v3//o3QCtwRSqKykY2HoObT5nvfnVoB8EKTpk1YcTX5VNmpzSzHWzMrJw56KYBY8wZxpgzgO8MPY7dPmyMucwY87c01ZlxNh6Dmy+Z17f1s3JjJ4UFwooTpx70cqz5kvlQaGY72JhZOePookPGmC8PPZbo/6KSMC+SgrqyjtvtznQJaZcPmcMRww9e2gXEdhCc8MEdBBPlQ+ZDpZntYGNm5YzTowYmi8gfRaQTCAHBhJsVKioqRn9RnsmHzI9s6GBLp4/6siL+8eiGUV+fD5kPlWa2g42ZlTNOjxr4KRAAlhG91sCxwEPAv6eorqzT2dmZ6RLSLtczdw8E+eVrzQBcvXTqiDsIJsr1zGOhme1gY2bljNNG4ETgX40xbwLGGPMWsBz4r/F8uIjcLSJtIvJuwrQaEfmriGyO3VcnzLtRRLaIyEYROSth+mIReSc27/tysI3AY1RdXT36i/JMrmf++at76QuEOW5qBSfNdHYMda5nHgvNbAcbMytnnDYCYaKbBAD2iUg90asPThnn5/8SOPuAaTcATxlj5gJPxZ4jIguBS4FFsffcKSJDP/F+DPwbMDd2O3CZ42bjoTe5nPnt5j6e3NxFkUuGPYPgSHI581hpZjvYmFk547QRWA18LPb4ceC3wIPAa+P5cGPM80DXAZPPB+6NPb4XuCBh+v3GmEFjzDai1z9YIiKTgEpjzCpjjAHuS3hP0vj9/mQvMuvlauZQJHqJYYBPHNnAlKpix+/N1czjoZntYGNm5YyjowaInkNgqGn4HHAtUA6k4hTDDcaYZgBjTLOITIxNnwIkHq64OzYtGHt84PSksvEY3FzN/Md329jR7WdypZtPHDX6DoKJcjXzeGhmO9iYWTnj9PDBfQmPfcBNKatoZMOt2x3uiohD0z+gra2N5cuXU1hYSDgc5sILL2TFihW0tLRQVlaGy+Wip6eH+vp6urq6MMZQX19Pa2srPp+P+vp6+vr6aGhooL29HRGhpqaG9vZ2KisrCYfD9Pf309jYSEtLC0VFRVRVVdHR0UFVVRWBQACfz7d/vtvtpqKigs7OTqqrq/H5fPj9/v3zS0pK8Hg8dHd3U1tbS29vL4FAYP98j8eD2+3G6/VSV1eH1+slGAzunz9apvLycoARM/l8Purq6nIq02BBCfe9Ht1B8KrFDbQ17zmkcfL7/fvrzpZMo43TeP/uDQ4OUlVVlVeZRhunQCBAeXl5XmUabZwCgQBlZWVjyqTym0TXpg8zQ+QbThZgjPnKuAoQmQn8xRhzeOz5RuD02NqAScCzxpj5InJj7PO+HXvd48DXgO3AM8aYptj0y2Lv//SBn7Vq1SrT1NQ0pjpbW1tpaDi0X5e5Lhczf/2vW3lph5eTZ07gK2fOOuT352Lm8dLMdhhP5jVr1ry+bNmy45JcksoSB9tHYJrDW7I9RPzUxVcAf06YfqmIFIvILKI7Bb4S24zQKyJLY0cLXJ7wnqTxeDzJXmTWy7XMf9vp5aUdXjxFBVx9wti2DuVa5mTQzHawMbNyZsRNA8aYf0n1h4vI/wKnA3Uishv4KvAd4AERWQ7sBC6O1bNWRB4A1hE9gmGFMSYcW9TVRI9A8AArY7ek6u7uprKyMtmLzWq5lNkfivCjl6O7ilx+7CTqy8Z2FrVcypwsmtkONmZWzjjaR0BEZo80zxizdawfboy5bIRZy0Z4/c3AzcNMfw04fKx1OFFbW5vKxWelXMr8mzdaaO0LMLvGwwWL6se8nFzKnCya2Q42ZlbOOD18cAuwOXa/JeH55hTVlXV6e3szXULa5Urm7d0+fvd2KwJ89uRpuArGfj6pXMmcTJrZDjZmVs44PWrgfQ2DiDQSXY3/QiqKykaBQCDTJaRdLmSOGMMdL+wibODcploWTCwb1/JyIXOyaWY72JhZOeN0jcD7GGNaiJ5P4NvJLSd72XgMbi5kXrmxk3Vt/dR4Cll+/ORxLy8XMiebZraDjZmVM2NqBGLmA6XJKiTb2Xgt72zP3DUQ5K5X9gJw9QlTKS92en6skWV75lTQzHawMbNyxunOgi/w/pP0lBI957+jcw3kAxsPvcn2zD/+2276A2GOn1rJqbMmJGWZ2Z45FTSzHWzMrJxx+hPqrgOe9wNvGWOs2VnQ7R7b4Wi5LJszr97p5bmt+yh2Cf950lTHFxUaTTZnThXNbAcbMytnnO4seO/or8pvXq+XCROS86szV2Rr5oFAmO+/FL2o0OWLJ9FY4fyiQqPJ1syppJntYGNm5YyjfQRExC0i3xCRzSLSH7u/SURKUl1gtqirq8t0CWmXrZnvea2Z9v4gc+s8XHj4xNHfcAiyNXMqaWY72JhZOeN0Z8EfAx8GrgGOj92fBtyZorqyjtfrzXQJaZeNmde19vPQunYKBL5wyvRxnTNgONmYOdU0sx1szKyccbqPwAXAnISrEK4TkdVETyz0rympLMsEg8FMl5B22ZY5GI5w+4s7McAnjmxgTm3yD1rJtszpoJntYGNm5YzTNQItfPBQQQ/QnNxyspeNx+BmW+bfvNnKjm4/UyqL+adjUlNbtmVOB81sBxszK2ecNgL/AzwmIleJyDki8m/Ao8B9IvLhoVvqysw8G4/BzabM73UOcP+bLQjwhVOnU1w4nlNgjCybMqeLZraDjZmVM043DXw6dv/FA6b/e+wG0fMMjHhxolxXVja+U9fmomzJHIoYbn1+J2ED5y+s54jG8pR9VrZkTifNbAcbMytnnB4+OCvVhWQ7l8uV6RLSLlsy//atVt7r9NFY4eZfj5+U0s/KlszppJntYGNm5Yzj9asiUigip4rIZSJyioiM/3yuOaSnpyfTJaRdNmTe1uXj129EV2l+/pTpeIpS+59ZNmRON81sBxszK2ecnmK4CXiY6A6Cu4BpgF9EzjPGrE9hfVmjvn7s17jPVZnOHAxH+N5zOwhFDOc21XLM5IqUf2amM2eCZraDjZmVM07XCNwJ/AyYZow5wRgzFfgJFp1HoKurK9MlpF2mM//mzVa2xDYJXLVkSlo+M9OZM0Ez28HGzMoZp43A0cBtxpjECw/dEZtuhfdHt0MmM29s7+d/Y0cJXHvqDErd6dm+qeNsB82sVJzTRmAv0TMJJjolNt0KNq5Wy1TmwVCE7z23k4iBCw+v58hJqTtK4EA6znbQzErFOW0Evgg8JCL3i8gtInI/8BAfPJwwb7W2tma6hLTLVOa7X93Lzn1+pk8o4VPHTU7rZ+s420EzKxXnqBEwxjwEHAO8C1TE7hcbY/6cwtqySnl5+n6VZotMZH5tdw9/XNuOS+C602ak7MRBI9FxtoNmViruoEcNiEgp8GXgcGAN8G1jzGA6ClP28fpD3PrcDiB6eeF59cm/loBSSqn3G+3n1g+B84ANwMeBW1NeUZbq6+vLdAlpl87Mxhhuf2EnXb4QRzSWc8mRDWn77EQ6znbQzErFjdYInAN81BhzXezx36W+pOzU0JCZL6ZMSmfmlRs7eXmHlzK3i+tOm5H0yws7peNsB82sVNxojUCZMaYZwBizC6hKfUnZqb29PdMlpF26Mm/r8nHnqt0AXHPSVBoq3Gn53OHoONtBMysVN9qZBQtF5AxARniOMebpVBWXTUQy8ws1k9KR2R+K8K2ntxMIG86aV8MZc2pS/pkHo+NsB82sVNxojUAbcHfC884Dnuf1FQcT1dRk9gsqE9KR+cerdrMjdqjgf5wwNeWfNxodZztoZqXiDrppwBgz0xgz6yA3K5oAsHO1WqozP/NeFys3dlLkEr704Zkpv6CQEzrOdtDMSsWl9yDtHFZZWZnpEtIulZl3dPu47YVdAFy9dCqzajwp+6xDoeNsB82sVJw2Ag6Fw+FMl5B2qco8EAjzjSe3MRiK8OE51ZzbVJuSzxkLHWc7aGal4rQRcKi/vz/TJaRdKjIbY7j9xZ3s8g4yo7qEz548Lat2YtJxtoNmVipOGwGHGhsbM11C2qUi85/WtvPc1n14igr4yrJZWbFfQCIdZztoZqXitBFwqKWlJdMlpF2yM7+xp5efrt4DwH+dMp1pE0qSuvxk0HG2g2ZWKk4bAYeKiooyXULaJTNzc+8g33x6GxEDnziqgVNnVydt2cmk42wHzaxUnDYCDlVV2XdSxWRl9gXDfP2vW+kdDHP81Eo+tXhSUpabCjrOdtDMSsVpI+BQR0dHpktIu2RkjhjD957bwdYuP1OrirnxjMxdR8AJHWc7aGal4rK2ERCR7SLyjoi8KSKvxabViMhfRWRz7L464fU3isgWEdkoImclux4bu+lkZL7n1b28uD16MaGvfWQ25cWjncwys3Sc7aCZlYrL2kYg5gxjzNHGmONiz28AnjLGzAWeij1HRBYClwKLgLOBO0UkqbujBwKBZC4uJ4w388oNHfz27TZcAl9ZNovpWbhz4IF0nO2gmZWKy/ZG4EDnA/fGHt8LXJAw/X5jzKAxZhuwBViSzA/2+XzJXFxOGE/mN/b08v2XomcOvOakaRwzpSJZZaWUjrMdNLNScdm8ntYAT4iIAX5qjPkZ0JBwWeRmEZkYe+0U4G8J790dm/Y+bW1tLF++nMLCQsLhMBdeeCErVqygpaWFsrIyXC4XPT091NfX09XVhTGG+vp6WltbKSkpobOzk76+PhoaGmhvb0dEqKmpob29ncrKSsLhMP39/TQ2NtLS0kJRURFVVVV0dHRQVVVFIBDA5/Ptn+92u6moqKCzs5Pq6mp8Ph9+v3///JKSEjweD93d3dTW1tLb20sgENg/3+Px4Ha78Xq91NXV4fV6CQaD++ePlqm8vBxgxEwej4eOjo5DztQRcvPVJ3cSNnDunHIWevoZHKzIikyjjVNpaSltbW05NU7j/btXVlZGa2trXmUabZwqKipobm7Oq0yjjVNlZSV79+4dUyaV38QYk+kahiUik40xe2Nf9n8F/hN4yBgzIeE13caYahH5EbDKGPOr2PRfAI8aY/6QuMxVq1aZpqamMdWzY8cOZsyYMdY4OWksmZt7B/n8Q5vo8oU4bfYEbjxjJgVZdObA0eg420EzH5o1a9a8vmzZsuNGf6XKRVm7acAYszd23wb8keiq/lYRmQQQu2+LvXw3MC3h7VOBvcmsx+12J3NxOeFQM+/zBfniyvfo8oU4ZnI5/+e0GTnVBICOsy00s1JxWdkIiEiZiFQMPQY+CrwLPARcEXvZFcCfY48fAi4VkWIRmQXMBV5JZk0VFbmxjTuZDiVz32CILz72Hnt6Bjms1sNXzpyN25WVf70OSsfZDppZqbhs/Z+6AXhRRN4i+oX+iDHmMeA7wEdEZDPwkdhzjDFrgQeAdcBjwApjTFIvtdXZ2ZnMxeUEp5l9wTBffnwrWzp9TK50882z5lDmzq5rCDil42wHzaxUXFbuLGiM2QocNcz0TmDZCO+5Gbg5VTVVV2fnKXFTyUnmwVCErzyxlXVt/UwsL+KWc+ZSU5q7pzLVcbaDZlYqLlvXCGQdGw+9GS3zYCjC15/cylvNfdSUFnLLOXNpqMjt7ZA6znbQzErFZeUagWzk9/szXULaHSyzPxThq09s5Y29vVSVFHLLOYcxpao4jdWlho6zHTSzUnHaCDhk47W8R8rsC4b5yhPRNQHVnkJu+dhhzKj2pLm61NBxtoNmVipONw04ZOO1vIfL3DcY4kuPvbd/c8D3zp3LzDxpAkDH2RaaWak4XSPgUElJ9p8nP9kOzNw5EORLj21ha5efurIivvexw5hSlV9/LjrOdtDMSsVpI+CQx5M/v3qdSsy8xzvIjY9toaU3wNSqYr5zzmFMLM/tHQOHY/s420IzKxWnmwYc6u7uznQJaTeUeX1bP59/eBMtvQHm15dy+3nz8rIJALvH2SaaWak4XSPgUG1tbaZLSLva2lqe39rNd5/bQSBsWDylgv+7bBalOXqyICdsHWfbaGal4nSNgEO9vb2ZLiGtjDH85o1mvvn0dgJhw8eaarnprDl53QSAfeMMmtkWNmZWzugaAYcCgUCmS0gbXzDMbc/v5LltfQhw5ZLJfPyIiUiOXUBoLGwa5yGa2Q42ZlbOaCPgkC3H4O7tGeTrf93Ktm4/nqICrj99BifOmDD6G/OELeOcSDPbwcbMyhndNOCQDcfgvrh9H5/500a2dfuZWlXMjcdVWNUEgB3jfCDNbAcbMytndI2AQ/l86E0gFOFnr+zhoXUdAJwwo4rrTptB/z77rlaWz+M8Es1sBxszK2e0EXDI7c7Pw+W2d/u45dkdvNfpo7BAuGrJZC5YVI+IEMzTzAeTr+N8MJrZDjZmVs7opgGHvF5vpktIqnDE8Lu3W1nxp4281+ljcqWbO86bxz8cHt8pMN8yO6GZ7aCZlYrTNQIO1dXVZbqEpNnt9XPbCzt5t6UfgHPm1/LpD035wKGB+ZTZKc1sB82sVJw2Ag55vV7KysoyXca4BMMRfvd2G79+s4Vg2FDjKeTzp0znQ9Orhn19PmQ+VJrZDppZqThtBBwKBoOZLmFc3m7u44cv72J7d/Sa5B+dW8O/fWgKlSUj/xXI9cxjoZntoJmVitNGwKFcPQa3rS/Az1fv4blt+wCYXFnMZ0+axjFTKkZ9b65mHg/NbAfNrFSc7izoUK4dg9sfCHPPa3tZ/rt1PLdtH8Uu4fJjG/nphU2OmgDIvczJoJntoJmVitM1Ag7lyra1QCjCXzZ08Js3WugZDANw2qwJXPWhKYd8xcBcyZxMmtkOmlmpOG0EHHK5svtiO4OhCI9u6OCBt9voHIhuCzy8sYwrj5/Cwoax/QeQ7ZlTQTPbQTMrFaeNgEM9PT1UV1dnuowP6A+EeWR9Bw++20aXLwTA7BoP/3LcJJZMqxzXhYKyNXMqaWY7aGal4rQRcKi+vj7TJbxPS+8gf17bzsqNnQwEIwAcVuvhn49t5ITpVUm5UmC2ZU4HzWwHzaxUnDYCDnV1dVFaWprRGiLG8PruXh5e387qnT2Y2PSjJ5fz8SMmcvzU8a0BOFA2ZE43zWwHzaxUnDYCDhljRn9RirT2BnhicydPbOqitS96TfGiAuGUWRO48IiJzKtLzT/uTGbOFM1sB82sVJw2Ag6le7Vajz/Ei9v38cx73bzd3Lf/1//E8iLObarj7Pm1VHuKUlqDjasSNbMdNLNScdoIONTa2sqMGTNS+hn7fEFW7ezh5e37eH1PL6FI9Ou/yCWcPHMCZ82r4ejJFRQkcfX/waQjc7bRzHbQzErFaSPgUHl5edKXaYxhe7efV3b1sHqXl3Wt/cS++ykQOHZKBR+eU81JMydQ5k7/oT+pyJztNLMdNLNScdoIpFnXQJA39vbyxp5e1uztpaM/fv7vwgLh2CnlnDxzAidMr6K6NLWr/pVSSiltBBzq6+ujtrb2kN5jjGFPzyDrWvtZ29rPOy197PYOvu81E0oKWTKtkiXTKlk8tTIjv/xHMpbMuU4z20EzKxWnjYBDDQ0NB51vjKFjIMjmjgE2d/jY1D7AhvZ+emOn+R1SUljA4Y1lHDu5gmOmVDCrxpO2bf6HarTM+Ugz20EzKxWnjYBD7e3tTJs2DYDewRA7u/3s2OdnR7efrV0+tnb5PvClD1DjKWTBxDIWNZRxxKRyDqstxVWQnV/8B0rMbAvNbAfNrFScNgIOPL+1m2c39rHvzU3s8g7i9YeGfV1FsYvDaj3MrStlbl0pTfVlTCwvSupJftIpV+seD81sB82sVJw2Ag681dzHi3v8+58Xu4RpE0qYWV3C9OoSZtd4mFXjoa40d7/0h1NTU5PpEtJOM9tBMysVl1eNgIicDfw/wAXcZYz5TjKWe9rsCVSYAY6cNYmpVSXUlRVl7Xb9ZGpvb7fuuGPNbAfNrFRc3jQCIuICfgR8BNgNvCoiDxlj1o132UdOqmBayUSqqyvHu6icUllpV17QzLbQzErFFWS6gCRaAmwxxmw1xgSA+4Hzk7XwcPiDOwLmO81sB81sBxszK2fyqRGYAuxKeL47Ni0p+vv7k7WonKGZ7aCZ7WBjZuVM3mwaAIbbaP++y221tbWxfPlyCgsLCYfDXHjhhaxYsYKWlhbKyspwuVz09PRQX19PV1cXxhjq6+tpbW2lpKSEzs5O+vr6aGhooL29HRGhpqaG9vZ2KisrCYfD9Pf309jYSEtLC0VFRVRVVdHR0UFVVRWBQACfz7d/vtvtpqKigs7OTqqrq/H5fPj9/v3zS0pK8Hg8dHd3U1tbS29vL4FAYP98j8eD2+3G6/VSV1eH1+slGAzunz9apqFTjo6UyePx0NHRkVeZRhun0tJS2tra8irTaONUVlZGa2trXmUabZwqKipobm7Oq0yjjVNlZSV79+4dUyaV3yRfLk0pIicAXzPGnBV7fiOAMebbQ69ZtWqVaWpqGtPyd+zYYd2ONprZDprZDuPJvGbNmteXLVt2XJJLUlkinzYNvArMFZFZIuIGLgUeStbC//SnPyVrUTlDM9tBM9vBxszKmbxpBIwxIeAzwOPAeuABY8zaZC3/wQcfTNaicoZmtoNmtoONmZUz+bSPAMaYR4FHU7HsUGj4swnmM81sB81sBxszK2fyZh8BJ5566ql2YMdY3tvV1VVXU1PTkeSSsppmtoNmtsM4M89YtmxZfVILUlnDqkZAKaWUUu+XN/sIKKWUUurQaSOglFJKWUwbAQdE5GwR2SgiW0TkhkzXM1YiMk1EnhGR9SKyVkQ+G5teIyJ/FZHNsfvqhPfcGMu9UUTOSpi+WETeic37vmT5ZRdFxCUib4jIX2LP8zqziEwQkd+LyIbYeJ9gQebPx/5evysi/ysiJfmWWUTuFpE2EXk3YVrSMopIsYj8NjZ9tYjMTGc+lSHGGL0d5Eb0SobvAbMBN/AWsDDTdY0xyyTg2NjjCmATsBD4LnBDbPoNwC2xxwtjeYuBWbE/B1ds3ivACUTP6LgSOCfT+UbJ/gXgN8BfYs/zOjNwL3Bl7LEbmJDPmYmeTnwb4Ik9fwD4VL5lBk4FjgXeTZiWtIzAfwA/iT2+FPhtpjPrLfU3XSMwupRezCidjDHNxpg1sce9RM+3MIVonntjL7sXuCD2+HzgfmPMoDFmG7AFWCIik4BKY8wqY4wB7kt4T9YRkanAucBdCZPzNrOIVBL9wvgFgDEmYIzZRx5njikEPCJSCJQCe8mzzMaY54GuAyYnM2Pisn4PLMumNSIqNbQRGF1KL2aUKbFVfscAq4EGY0wzRJsFYGLsZSNlnxJ7fOD0bHUHcB0QSZiWz5lnA+3APbHNIXeJSBl5nNkYswe4FdgJNANeY8wT5HHmBMnMuP89JnqSNi9Qm7LKVVbQRmB0o17MKNeISDnwB+Bzxpieg710mGnmINOzjoj8HdBmjHnd6VuGmZZTmYn+Mj4W+LEx5hign+gq45HkfObYdvHzia4CnwyUicg/H+wtw0zLqcwOjCVjPuVXDmkjMLrdwLSE51OJrnLMSSJSRLQJ+LUxZuico62x1YXE7tti00fKvjv2+MDp2egk4O9FZDvRzTofFpFfkd+ZdwO7jTGrY89/T7QxyOfMZwLbjDHtxpgg8CBwIvmdeUgyM+5/T2wTSxUf3BSh8ow2AqNL6cWM0im2re8XwHpjzG0Jsx4Crog9vgL4c8L0S2N7Es8C5gKvxFY/9orI0tgyL094T1YxxtxojJlqjJlJdOyeNsb8M/mduQXYJSLzY5OWAevI48xENwksFZHSWK3LiO4Dk8+ZhyQzY+KyPk7034uuEch3md5bMRduwMeI7mH/HvClTNczjhwnE13N9zbwZuz2MaLbAJ8CNsfuaxLe86VY7o0k7D0NHAe8G5v3Q2JnqczmG3A68aMG8jozcDTwWmys/wRUW5D568CGWL3/Q3Rv+bzKDPwv0X0ggkR/vS9PZkagBPgd0R0LXwFmZzqz3lJ/01MMK6WUUhbTTQNKKaWUxbQRUEoppSymjYBSSillMW0ElFJKKYtpI6CUUkpZTBsBpZRSymLaCCillFIW00ZAKaWUspg2AkoppZTFtBFQSimlLKaNgFJKKWUxbQSUSiMRWSsip2e6DqWUGlKY6QKUyici0pfwtBQYBMKx5582xixKf1VKKTUyvfqgUikiItuBK40xT2a6FqWUGoluGlAqjURku4icGXv8NRH5nYj8SkR6ReQdEZknIjeKrIQbYwAAATRJREFUSJuI7BKRjya8t0pEfiEizSKyR0S+KSKuzKVRSuUDbQSUyqzzgP8BqoE3gMeJ/rucAnwD+GnCa+8FQsBhwDHAR4Er01msUir/aCOgVGa9YIx53BgTAn4H1APfMcYEgfuBmSIyQUQagHOAzxlj+o0xbcDtwKUZq1wplRd0Z0GlMqs14bEP6DDGhBOeA5QDk4EioFlEhl5fAOxKR5FKqfyljYBSuWEX0SMQ6mJrD5RSKil004BSOcAY0ww8Afy3iFSKSIGIzBGR0zJdm1Iqt2kjoFTuuBxwA+uAbuD3wKSMVqSUynl6HgGllFLKYrpGQCmllLKYNgJKKaWUxbQRUEoppSymjYBSSillMW0ElFJKKYtpI6CUUkpZTBsBpZRSymLaCCillFIW00ZAKaWUstj/B3YmJ9fqdZixAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load required packages\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('bmh')\n",
    "from scipy import integrate\n",
    "\n",
    "# set parameters\n",
    "birth_rate = 0.0011\n",
    "death_rate = 0.0001\n",
    "init_size = 50\n",
    "xdim = 50\n",
    "ydim = 50\n",
    "k = xdim * ydim\n",
    "time = 10000\n",
    "de_times = np.linspace(0, time, 100) # inference times\n",
    "\n",
    "# create exponential growth equation as a function of N\n",
    "def dN_dt(N, t=0):\n",
    "    return birth_rate * (1 - N[0]/k) * N[0] - death_rate * N[0]\n",
    "\n",
    "# integrate growth equation\n",
    "X, info = integrate.odeint(dN_dt, init_size, de_times, full_output=True)\n",
    "# format and store results\n",
    "de_pop, = X.T\n",
    "\n",
    "# plot\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(de_times, de_pop, label='Integrated ODE')\n",
    "ax.set_ylabel('Population Size', labelpad=10)\n",
    "ax.set_xlabel('Time', labelpad=10)\n",
    "plt.legend(bbox_to_anchor=(1,1), loc='upper left', borderaxespad=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implicit IEBM\n",
    "\n",
    "One way to model logistic growth is applying a very simple limit to a population size. For whatever phenomological reason, the population is not allowed to grow beyond this limit. Individual birth events are added to the event heap at the same birth rate, *b*. But when birth events are handled, there's a chance $(1-N/k)$ they do not occur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style type='text/css'>\n",
       ".datatable table.frame { margin-bottom: 0; }\n",
       ".datatable table.frame thead { border-bottom: none; }\n",
       ".datatable table.frame tr.coltypes td {  color: #FFFFFF;  line-height: 6px;  padding: 0 0.5em;}\n",
       ".datatable .boolean { background: #DDDD99; }\n",
       ".datatable .object  { background: #565656; }\n",
       ".datatable .integer { background: #5D9E5D; }\n",
       ".datatable .float   { background: #4040CC; }\n",
       ".datatable .string  { background: #CC4040; }\n",
       ".datatable .row_index {  background: var(--jp-border-color3);  border-right: 1px solid var(--jp-border-color0);  color: var(--jp-ui-font-color3);  font-size: 9px;}\n",
       ".datatable .frame tr.coltypes .row_index {  background: var(--jp-border-color0);}\n",
       ".datatable th:nth-child(2) { padding-left: 12px; }\n",
       ".datatable .hellipsis {  color: var(--jp-cell-editor-border-color);}\n",
       ".datatable .vellipsis {  background: var(--jp-layout-color0);  color: var(--jp-cell-editor-border-color);}\n",
       ".datatable .na {  color: var(--jp-cell-editor-border-color);  font-size: 80%;}\n",
       ".datatable .footer { font-size: 9px; }\n",
       ".datatable .frame_dimensions {  background: var(--jp-border-color3);  border-top: 1px solid var(--jp-border-color0);  color: var(--jp-ui-font-color3);  display: inline-block;  opacity: 0.6;  padding: 1px 10px 1px 5px;}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import sys\n",
    "sys.path.append('../')\n",
    "from iebm.populations.population2D import Population2D\n",
    "from iebm.traits.static_trait import StaticTrait\n",
    "from iebm.events.birth import BirthEvent, BirthDiffusionEvent\n",
    "from iebm.events.death import DeathEvent\n",
    "from iebm.simulation import Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can recreate the exponential growth IEBM and add an implicit_capacity to the population \n",
    "# First lets create an empty 2D population. Each individual contains only a unique ID and x-y positions\n",
    "log_pop = Population2D(name='log_pop',\n",
    "                       init_size=init_size, \n",
    "                       implicit_capacity=k, # limit the population growth\n",
    "                       xdim=xdim, ydim=ydim)\n",
    "\n",
    "# create the same traits as the exponential growth example\n",
    "trait_list = [(StaticTrait, {'name' : 'birth_rate', \n",
    "                             'value' : birth_rate}), \n",
    "              (StaticTrait, {'name' : 'death_rate', \n",
    "                             'value' : death_rate})]\n",
    "log_pop.add_traits(trait_list)\n",
    "# create the same events as the exponential growth example\n",
    "event_list = [(BirthEvent, {'name' : 'birth',\n",
    "                            'is_primary' : True,\n",
    "                            'current_time' : 0}),\n",
    "              (DeathEvent, {'name' : 'death', \n",
    "                            'is_primary' : True,\n",
    "                            'current_time' : 0})]\n",
    "log_pop.add_events(event_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 36.8 s, sys: 50.1 ms, total: 36.8 s\n",
      "Wall time: 4.65 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Add the exp_pop to a dictionary \n",
    "pop_dict = {str(log_pop) : log_pop}\n",
    "\n",
    "# Create the simulation with the population dictionary\n",
    "sim = Simulation(pop_dict)\n",
    "\n",
    "# Run the simulation\n",
    "sim.run(time)\n",
    "\n",
    "# Store results\n",
    "imp_pop = sim.population_history[str(log_pop)]\n",
    "imp_times = sim.time_history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAEMCAYAAADd1S/AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3hUVd6A35NJJr1XAiH0ElCai4gKKKLYFXXVda3ori6W/dRFLKu7dnF1dW2r4lrXdXVt2MBeEEQRC9JbgADpIXUmk5k53x8zGRKSkEsy9Z7zPs883HvOnTu/1zPmN/feU4SUEo1Go9FoNJFPVKgD0Gg0Go1G4x90UtdoNBqNxiTopK7RaDQajUnQSV2j0Wg0GpOgk7pGo9FoNCZBJ3WNRqPRaExCdKgDCCaff/65jI2NDXUYGo1GE1E0NTVV5uXlPRAdHT0SfTEYatxOp3Ot3W6fP2HCBPe+lUol9djYWEaMGNGj927bto3CwkI/RxTeaGc10M5q0BvnDRs2NGdmZv46Nja2QxLRBJ/m5uaDqqqqAO7dt07/4jJITExMqEMIOtpZDbSzGvTG2Wq1xumEHj7Exsa6vXdNOqCTukFSU1NDHULQ0c5qoJ3VoDfOQgg/RqLxE53mb53UDVJZWRnqEIKOdlYD7awGke5cXFwcc9pppw0++OCDR48ePXr0lVdeWdDc3CwWL16cXFBQMPbQQw8tGjNmzOijjz56+Ouvv+77BXPzzTfnDxs27OCJEycWtb6qqqosoXQJJEo9U+8N+pe9GmhnNdDOkYXb7ea8884bcv7555dffvnlm51OJ5dddlnhvHnz+s6cObN27NixDe+8884mgG+//Tb+oosuGhIfH198wgkn1ANcdNFFZTfddFNZaC2Cg75SN4jD4Qh1CEFHO6uBdlaDSHZetGhRcmxsrPvyyy+vAoiOjubBBx/c8eabb2Y1NTW1y2MTJ060XXXVVbsWLFiQE5poQ4u+UjeIzWYLdQhBRzurgXZWA385H7vghwl+OdE+fHjpuO+7qluzZk38qFGjmtqWpaenu3Nzcx2bN2/uME55woQJTU899VRe6/5zzz2X+9Zbb2UCpKSkOD/++OMN/ow9nNBJ3SB5eXndH2QytLMaaOfQ4Ha00LhlBzHpKcTlZnV7vJQStwSnW7Z7uVpfsnUbnN5tt1vikuCWkuiYyL39LqVECNFhnXBveafHt0Wl2+86qRuktLRUuXGt2lkNtHPPkFJid7qxtbixtbg8/zo92/YWN3an59Xc5uUqLiFjwb9oiU8gffnyDucsGzmK8oFDOOj9t9ky7lcsOfVsGuMTPQncJemQ1YA+27dQnZ1Hc3zCfuMdnBrNE2el9coZ9n9FHSiKiopsH3zwQXrbspqamqjy8nLroEGDmj///PN2x69cuTJh4MCB9mDGGC7opG4Qq9Ua6hCCjnZWA9WdW1xuau3ONi8X9Y126lyChmYn9c0uGppdNDhcNDqcNNqdJGzdyo6sPrgRzHr+Ufpv2cBLf7iB8vz+AFicLRz7xkuM/HkFVmDLxCMY8+2SbuPKXbua3LWrARj0w3cM+uE7mhKT2DziYKpy8pj2wRtU5+RRn5FFZcEARn39GXFNje3OsXPsBLZNPYqaolFYhCBKgCVKkGWN3GHmxx9/fP1dd90V9fTTT2dedtllVU6nk+uvv77g1FNPrUxISGgntmLFivhHH300//777y8OUbghRex7m8LMLFu2TPZ0RrmGhgaSkpL8HFF4o53VwKzOLS43FY0OKhudVDY6qN64Hferb1F8wknstiZR63DTUNtInUuQWF/Lmc8+QnpVeZfnc0ZHY0tIIrluT5fHlBx+JM7MTAYsfKvb+MTQQST/cz7Ol18nOj6OukefBSAqKRF3Q2M37zZG/pkz2fW/RZ7zxsVybPFnPTrPtm3bqgoLC4v9ElQP2bp1a8w111xTWFxcHOd2u5kyZUrtAw88UPLll18mXXrppYPz8/Mddrs9KiMjo2XOnDmlZ555Zi14hrS99tprWWlpac7Wc/373//eNHTo0MjtOQiUl5evHT58+Pn7luukbhA9raQaaOfQ07RtJ3F9ctj+wps0rNlM/q+PJ21cEVGx7e8oSCmpaWphZ3UjO+9+nObtu2iyxiIrq7ELC/kb1obIoHMs8XGMuP1qcmZOwZqZhoja/+AjR00dO//zLo1btlPxyTKad1d0OEZYYxj71B00rN+Ky2andOGnNG3Z0fnnZ6QyY80HPYo9HJK6pj1dJXV9+90g6enp3R9kMrSzGgTSuWrJCmzbd9OwfivOJhuj7vtTh2RW/NR/WXfrw12eo+Tld9rtuy0WolyuDsfFel+9JWn4QPLPPI6y976g9sf2Pwyypx9G9jGTyZw6kXW3/QNnfSMHPXwLboeDpcdchLt578XfwKvOZ/jNVwAg3e5uk/i+WNNTGPiH3xg6NnfmFACGzfs9AHu+/4UdL7xFbF4W0uVGOl1kXzLrgD5fE5nopG4Qm81GSkpKqMMIKtpZDfzpLKWk7uf1/Hzl7TRuLO5QX/Li273+jM4Suo/xB8PKnzutKrrvT/Q772TsO8vYsfR7kmPjiMvLpmH9FpIPGkbahNHtelIPuuqC/cYx4cX72+0fu+3zLo890ITeW9ImjCZtwuh2ZWVlSnT+Vh6d1A1it6vXkVI7q0FvnaXbzY4X32bNDfd3f7ABVo89lOqcPFaPn0RTYjLWZjvj1qykaM0PpGzfhqXJM1w5OjONzMnjiUlNYuAVvyFxcH9D508o7EscTvK9jxwyJo/zS9zhjorfbRXRSd0g4TCuNdhoZzXoqbOUkk9HnUhLddcdx8Y8eQepY4Zjzc5g6+MvU/XNTzSmpdPy3se+Y34Zfxgfzvqtbz83ycrQrATOzopnaFYCgzLiyUiY7Kt3O51ERffuT5duZ41Z0UndIHosrxpoZ2PsfO0DVl11R4fy/F+fwMA//Aa3zU7quCKklGyptrF8Qx3fjZjC2owJuCVw2OkAREcJhmcncHZuIiNzExmZnUh6wv6XCO1tQgfdzhrzopO6QeLi4kIdQtDRzmpwoM7lHy7pkNCP3f4FUVZPMna5JT+XNvD10h0s3VZLZWOL77goAUU5iYzrm8yYPkmMzEkkNjr4S1DodtaYFb2gi0Hi4+NDHULQ0c5qYMTZ3exg+7OvsyhvMisvmOsrH//i/cwsXQox0fy8u56Hlmzn7H+v4ob3N7FwTSWVjS1kJERz/PBMbj1mIK+ffzAPnTKMCyf0YWx+ckgSOuh2jkTy8/P90vlh8eLFySeffPIQgNdeey319ttv3+9zialTp44A2Lhxo/XZZ5/N6OyYjRs3WidMmDCq9fwFBQVj2y71+t577yUDZGVlTWgtmzRp0sjPP/88sfX9GRkZE+bNm5ffes6ysrLo7Ozs8VdccYWxziJe9JW6QWpqapTrFa2d1aA755a6Bj4ZdmyH8mO3fc5uu5tnvtvFJ5uq212R90uN5YgBaRw+IJWhWQlEdTI/dyjR7awBOOuss2qB2v0d88UXX6wD2LJlS+wbb7yRcfHFF1d3d962S8G2xWq1ur/99ts1AAsXLky54447+k2bNm09QH5+fvNnn32WBuwCeOWVV9J7MtWtTuoGyczMDHUIQUc7q0F3zvsm9Ng+Obj+u4A/fbSNVaUNvvLcJCvTBqdz1KB0BmbEdbrQRrig2zlyWbx4cfL8+fPzMzMzW9avX58wY8aMmqKiItszzzyT29zcLF566aXNI0aMaL7gggsGxMbGujdt2hRfXV0dc+utt+4444wz2iXwJ598MvPHH39MfOKJJ7bv3Lkz+qqrriosKSmJBZg/f/62adOmNebn54/btWvXD3feeWff4uLiuIkTJxadccYZlTfccEPX0w8aoK6uzpKSkuKb5S42NlYOGjTItnTp0oTJkyc3LVy4MOPEE0+sLi0tPaB5nHVSN0h9fb0pp9LcH9pZDfZ1lm43a//8ENuf+V+74/KuOI8fTjyd99dVUfulZ9ay2Ogopg1K49hhmYzKTQy7K/Ku0O3ccxblTQ7I0qszS5caXihm48aN8c8999yWrKws55gxYw6yWq2VX3/99dr7778/57HHHst55JFHdgDs3Lkz9uOPP16/fv362NNPP3348ccfv6qrc1577bX9J02aVD937tzNTqeTuro6S9v6W265Zeejjz6a29kV+L78+OOPSRMnTixq3X/hhRc2jxgxotnhcERNnDixyOFwiKqqqphXX3213RKws2bNqn711Vcz8vPzWywWi8zLy2vRST1AOBwRPU1wj9DOatDWuatb7RUnncjD/Sbj+tEzgcmgjHhOKcpi2qB0EqyWDseHO6q3c6RTVFTUWFBQ0ALQr1+/5mOOOaYWYPTo0balS5cmtx538sknV1ssFoqKipr79u3b/Msvv3TZW/Dbb79Nfu6557YCREdHk5GRsZ9ZjvaPkdvvX3zxReKcOXMGrlixYnVr/SmnnFJ3//33983Ozm45+eSTu73N3xk6qRtExTGe2lkNYtYUs+jQc7usXzTrt6wZfxhRwJSBaZw2KptRuYlhfXu9O1RsZ385H8gVdaCIiYnxLVoihCAuLk4CREVF4XQ6Rdu6toTTd3bq1KmNtbW10aWlpb48HBsbK4uKipqeffbZvKVLl/7y5ptvHvBaubr3u0FKS0tDHULQ0c7mxmVvZt1fH+Xni2/sULft9Fk8eOdjPHjnY2z+1WROLcri2bOKuGX6QEbnJYXVH8eeoFI7t6Ki87vvvpvucrlYt25d7M6dO2NHjRrVZcezQw89tP6RRx7JBnA6ndTU1LTLjykpKa6mpia/3ZZatWpVnNvtJicnx9m2/JprrimdO3duSU5OTo/uFOgrdYNE+nCQnqCdzYmzoZHqpT+0G5rWSsmZZ/JG0WSc1ljioqM4tSiLMw7KIS1+/xPCRBoqtPO+qOg8cODA5mOOOWZ4dXV1zN13370tISGhy2VJH3jgge1z5swpPOSQQ7KioqK4//77t02dOtW3Bu6ECRNsFotFGukot+8z9WuuuWb3eeedV9P6TL3NZxZH7zOZ0rhx4+zjxo3r8Zy+eulVg+zZs4e0tAO+ExLRaGdzsXru/ex44c0O5Vmzz2L5iWfy7tpKXBKsFsFpo7I504TJvBUzt3NX9MY5EpdeveCCCwYcd9xxteedd15NqGMJBHrp1V5SW1ur3B8B7Wwetj//ZqcJ3TX/r9wls2lcU0mUgOOGZXDBhD5kJx5Qh9uIw6ztvD9UdFYRndQNkpWVFeoQgo52Ngf167a0W0Ht4EdvpW7ir/jHinI2VdsByZg+SfzhsH4MzFDjFq0Z27k7VHN+4YUXikMdQygIaUc5IUSBEOIzIcRaIcRqIcQ13vIMIcRHQoiN3n/T27znRiHEJiHEeiHEcW3KJwghVnnr/iH83JOntna/kw6ZEu0cuUi3m4+GHMOivMl8Pc2zAlraIaOZsv1L3i04iD9+tJ1N1XayE2P4w7hU5p8wRJmEDuZp5wNBRWcVCXXvdydwnZRyJDAJmCOEKALmAZ9IKYcCn3j38dadA4wCZgKPCyFaeyM+AfwOGOp9zfRnoC0tLd0fZDK0c2TisjWzOP8IXA1N7cpj75jH5W+u53+rPP17zhidzYIzRzIuyxLxvdkPFDO084HSG2eV+l5FEO7OCkOa1KWUu6WUK73b9cBaoC9wKvC897DngdO826cCr0gpm6WUW4FNwEQhRB8gRUq5THq+fS+0eY9f0ONa1SCSnV22ZkpeeY/PxpzcrtySGE/J8wuY930tpfUOBmXE8fApw/j9pH7Ex1gi2rmnaOcDw+Fw2Jubm0N9Eajx0tzcHOV0Otd2Vhc2z9SFEAOAccByIFdKuRs8iV8IkeM9rC/wTZu3lXjLWrzb+5b7DRXXItbOkcPaWx9m21P/bVc2dsFd2A+fxN2fFlO8w06UgHPH5HLeuDxiLHv/Pkeqc2/QzgdGQ0NDaVVV1f+io6NHEvo7vKrjdjqda+12+/zOKsMiqQshkoDXgT9KKev2cyuwswq5n/J2lJeXM3v2bKKjo3G5XMyaNYs5c+ZQWlpKYmIiFouFuro6srOzqa6uRkpJdnY2ZWVluN1uqqqqaGhoIDc3l4qKCoQQZGRkUFFRQUpKCi6Xi8bGRvLy8igtLSUmJobU1FQqKytJTU3F4XBgs9l89VarleTkZKqqqkhPT8dms2G32331cXFxxMfHU1NTQ2ZmJvX19TgcDl99fHw8VquV2tpasrKyqK2tpaWlxVffnVPrPNBdOUkpqaysNJVTd+3U+j2JJKfSWx6h5qOl7b7rhX+fx/L+w3n6zXU43JCXaOGiokQOH5HJrpId7ZyEEJSVlYWVU6D/f7JYLOzevdtUTt21k8ViYdeuXT1yAhg1atS9Xf1h1oQPIR+nLoSIAd4FFkspH/SWrQemea/S+wCfSymHCyFuBJBS3uM9bjHwF6AY+ExKOcJbfq73/b9v+1m9GadeU1NDenp69weaCO0c3ribHay/83G2Pf2qr2zS+wuIPWg4j3y9g483eYbnzhiawZWTPbfaOyOSnP2Fdj4wVq5c+f306dMP8XNImgAQ6t7vAngGWNua0L0sBC70bl8IvN2m/BwhRKwQYiCeDnHfem/V1wshJnnPeUGb9/iFuro6f54uItDO4Ynb0UL9ui18WDitXUI/budX2IYM5o8LN/Dxphpio6O4fkp//jS1sMuEDpHh7G+0s8ashPr2++HA+cAqIcSP3rKbgHuBV4UQs4HtwFkAUsrVQohXgTV4es7PkVK2zo97BfAcEA984H35jezsbH+eLiLQzuGFdLn4cODRSEfHXsxTv3+TlbsbufuzYuqbXfRLjeW2YwZSmN79MLVwdg4U2lljVkKa1KWUS+j8eTjA9C7ecxdwVyflK4DR/ouuPdXV1SQkJATq9GGJdg4vfrriLx0S+ui/30S/c0/izV/KeXL5ZtwSDi1IYd5RA0g0uCRqODsHCu2sMSuhvlKPGELd9yAUaOfwQLpc/Hzl7ZQu/MRXVnTPdfQ54ziikhJ5bGkJb6+pAOC8cXmcPz6PqAMYdx6OzoFGO2vMik7qBlHx1pV2Dj2uJjsfDTrat58xeTwT33gUAFuLi3s+3sI32+uIiRJcN6U/Rw/JOODPCDfnYKCdNWZFjzc0SFlZWahDCDraObRIt5uPR/hmQkZEW5jw0t8AqLM7mfv+Jr7ZXkdyrIV7TxjSo4QO4eUcLLSzxqzopG6QtuM1VUE7h449K1ezOP8I3zP0YTdfznElX2FJiKOy0cF1725kfUUTuUlWHj5lGAfl9TzucHEOJtpZY1Z0Utdowgi308mOl97mmxMu85UNmXsZg666AICdtXb+752NbNtjpzA9jodOHka/1LhQhavRaMIM/UzdIA0NDWRmZoY6jKCinYNLU3EJX076dbuyX732DzKP9Mz5sX2Pnbnvb6S6ycnInATuOHYwKXG9/19Yt7MaqOisIjqpGyQ3NzfUIQQd7Rw8GjYWs+TI37Qrm/CfB30JfVuNjbnvb6LG5mRMnyRuP3bQfieUORB0O6uBis4qom+/G6SioiLUIQQd7Rx43I4Wfvzdn9sl9MM+WMCxJV+SfdQkALZW2/jTe56EPi4/mTuOG+y3hA66nVVBRWcV0VfqBlFtvWnQzsHgw/5T2+0ftvhfpI7Zuz5BSa2dG97fxB67kwl9k/nLjEHERvv3t7huZzVQ0VlFdFI3SEZGz4YLRTLaObCsveXvvu2Egf04csl/EJa9V+Cl9c3M9Sb0cfnJ/HXGIKx+Tuig21kVVHRWEX373SAq3rrSzoGjfu1mti14zbc/Zdmr7RJ6VWML8z7YRGVjC6NyE/nLjIEBSeig21kVVHRWEZ3UDZKSkhLqEIKOdg4MzkYbXx91PgD5vz6BmaXt10Kvb3Zy46JN7KpzMDQrnjv9/Ax9X3Q7q4GKziqik7pBXC5X9weZDO0cGL6afLZve/D/XdSuzuF085ePtlJcY6cgNZa7Zw4xvDBLT9HtrAYqOquITuoGaWxsDHUIQUc7+59db3xIc1klALknTCVxYD9fncstuffzbawqbSAzIYZ7jh9Cqh/GoXeHbmc1UNFZRXRSN0heXl6oQwg62tm/NG4t4ec//AWAnOOOYNy/7vHVSSn55zclLCneQ6LVwt0zB5OTZA1YLG3R7awGKjqriE7qBiktLQ11CEFHO/uXr48+37c9dsHd7ereWl3B22sqiYkS/HXGQAZmxAcsjn3R7awGKjqriE7qBomJiQl1CEFHO/uPpm07cduaAZj41uNExey9rb58ey1PLt8JwLVT+nNwn+SAxNAVup3VQEVnFdFJ3SCpqamhDiHoaGf/UL92M9+c9HsA0iYeTMaksb66rdU27v6sGLeE347LY3oPl0/tDbqd1UBFZxXRSd0glZWVoQ4h6Gjn3uNqsrP8lMtxVFQDMOre6311e2wt3PrhFmwtbqYNSuP88aF55qnbWQ1UdFYRPaOcQVT8laude0/pe5/hrPf0Oj78sxdJHjkYAKdbctenxZQ1OBiencD1UwpDNo2nbmc1UNFZRfSVukEcDkeoQwg62rl3SJeL4n++AsDwW6/0JXSAp5fv5KfdDWTER3PbMYGbLc4Iup3VQEVnFdFJ3SA2my3UIQQd7dw7dr3+IfWrNxKXn0P/i8/wlX+4oYo3V1cQHSX48zEDyUoMztC1rtDtrAYqOquITuoGUXGMp3buOS57M6uuvgOAwt+djSU+FoDNVU08/PUOAOZM7seo3CS/fF5v0O2sBio6q4hO6gZRcYyndu4ZLbX1fDTgKN9+/wtOB6Ch2cntH2+lxSU5fngmJ47I6vVn+QPdzmqgorOK6KRuEKs1tLdIQ4F2PnCc9Y18Mvw43/7Qmy7HkhCHlJL7v9zO7noHQzLjmXNYv/2cJbjodlYDFZ1VRPd+N0hycnAnBAkHtPOBs+K863zbYxfcRd5Jniv2/60qZ9m2WhKtFv48PbQd4/ZFt7MaqOisIob+sggPlwkhPhVC/OwtmyKE+HVgwwsfqqqqQh1C0NHOB8a2Z/7Hnm9/BiB1wihfQl9T1sgz3+0CYO7UQvqkxPY+UD+i21kNVHRWEaOXC7cDs4GngP7eshLghkAEFY6kp6eHOoSgo52NU774K9be/KBv/9C3nwA8a6Pf450x7syDcjisMPzGCut2VgMVnVXEaFK/CDhJSvkKIL1lW4FBgQgqHFFxOIh2Ns7G+Qt820f/8h5R0dFIKfn7Vzt8E8xcfEgff4XpV3Q7q4GKzipi9Jm6BWjwbrcm9aQ2ZabHbreHOoSgo52NUfHJMupXb0REWzh69fvEpHqeXb63roolxXtIiInixqMGEGMJn+fobdHtrAYqOquI0b8y7wMPCiFiwfOMHbgDeCdQgYUbKo7x1M7d47I38723c1zB+af5EnpxjY1/flMCwDVH9Cc/zJ6jt0W3sxqo6KwiRpP6tUA+UAuk4rlCLwTmBSiusEPFMZ7auXvajkfvd97JADhcbu77fBsOl+S4YRkcNTi8n2XqdlYDFZ1VxNDtdyllHXCaECIHTzLfIaVU6hsSFxcX6hCCjnbeP0uPvdi3nXP8FFJGDwPghe93s7nKRp9kK1dMCp/x6F2h21kNVHRWEaND2v4lhDhOSlkupfyuNaELIR4PbHjhQ3x8fKhDCDrauWtsJaXU/bwe8KyRPv7ZewH4aVc9r/1cTpSAG6YNIMFqCVis/kK3sxqo6KwiRm+//xZ4RghxfSflSlBTUxPqEIKOdu6a3W997Ns+9C3Pb9uGZifzv9iGBH4zNo+i3MRAhOh3dDurgYrOKmI0qduBScA5QogXhRCt8w2GZgHoEJCZmRnqEIKOdu6a8kVfAjDq/rmIKM//Rk98s5OKxhaGZyfwm3GR0ylJt7MaqOisIobH2EgpS4Aj8Qxv+1oI0Ze9w9t6hPe2frkQ4pc2ZX8RQuwUQvzofZ3Qpu5GIcQmIcR6IcRxbconCCFWeev+4e2d71fq6+v9fcqwRzt3TtXXK9mzwvOV7TPL8zVctq2WjzZWY7UI5k4tJDoqcn7v6nZWAxWdVcRoUhcAUkqblPI3wOvAt0Bvx+k8B8zspPzvUsqx3tf7AEKIIuAcYJT3PY8LIVofWD4B/A4Y6n11ds5e4XA4/H3KsEc7d853Z13t245OjKfW7uShJdsBuORX+RSkRVaHJN3OaqCis4ocyDSxPqSU9wKXAi/35sOllF8C1QYPPxV4RUrZLKXcCmwCJgoh+gApUsplUkoJvACc1pu4OkPFMZ7auSOL8iaD2w14FmwBePTrHdTYnBycl8Rpo7IDHqO/0e2sBio6q4jRIW33d1L2AfCB3yPycKUQ4gJgBXCdlLIG6At80+aYEm9Zi3d73/IOlJeXM3v2bKKjo3G5XMyaNYs5c+ZQWlpKYmIiFouFuro6srOzqa6uRkpJdnY2ZWVl2Gw2srOzaWhoIDc3l4qKCoQQZGRkUFFRQUpKCi6Xi8bGRvLy8igtLSUmJobU1FQqKytJTU3F4XBgs9l89VarleTkZKqqqkhPT8dms2G32331cXFxxMfHU1NTQ2ZmJvX19TgcDl99fHw8VquV2tpasrKyqK2tpaWlxVffnVNSUhJAl042m42srCxTOXXXTna73Rf3vk5rbn3Y912KHVxA+ozJvP7tRr7Y2kBctOCcIRaaGhvDzqm7dmpubiY1NTWi2qm33z2Hw0FSUpKpnLprJ4fDQWJiYo+cNJGD8FzcdlIhxCIp5Uzv9ld08fxcSjmlVwEIMQB4V0o52rufC1R6P+8OoI+U8hIhxGPAMinlS97jnsEz09124B4p5THe8iOBuVLKk/f9rGXLlskRI0b0KM7y8nJycnJ69N5IRTvvxV5WyedjTvHtH7vjS+pdcOn/1lJrd3LV5H6cXBR5V+mg21kVeuO8cuXK76dPn36In0PSBID9Xam/0GZ7QZdH+RkpZVnrthDiaeBd724JUNDm0H7ALm95v07K/YrVau3+IJOhnffSNqEft/trhBA8vqSYWruTMX2SOHFkVrBC9Du6ndVARWcV6TKpSylfbrP9fHDCASFEHynlbu/u6UBrz/iFwMtCiCa076QAACAASURBVAfxTFk7FPhWSukSQtQLISYBy4ELgEf8HVdtbS1paWn+Pm1Yo509NG3b6dseMvcyhBAs3baHzzbXEBsdxf8d2Z8o/w+4CBq6ndVARWcV2e8zdSHEBKBZSvmLdz8beAgYDSwDrpdS9nilNiHEf4BpQJYQogS4DZgmhBiL5/Z7MfB7ACnlaiHEq8AawAnMkVK6vKe6Ak9P+ng8z/n9/qw/Kytyr8R6inYGZ30jXx56lm9/8P9dRH2zk398vQOASw7pE9aLtRhBt7MaqOisIt11lHsI+Ct7r5YX4LlKfgo4F5gP/KGnHy6lPLeT4mf2c/xdwF2dlK/A80MjYNTW1pKYGBkzhPkL1Z2llHwycu/oyIMfuw0hBE8v30V1k5OinEROidDn6G1RvZ1VQUVnFekuqY8EvgIQQqQBxwOjpZQbhBALgaX0IqlHEi0tLaEOIeio7OxucfJhwd4+oHmnTCf/jOP4cVc9izZUERMluPbI/lgiaJKZrlC5nVVCRWcV6W6cejTQOmPBJKBUSrkBQEq5A1DmAY2KYzxVdXbZmlky5Te+srRDRjP2qTtodrp5aInntvu54/Lonx5Zk8x0hartrBoqOqtId0l9NdD6QPEcwLeKhXea2NoAxRV2qLgWsarOa29+kKatnqkP+p59ApPefQqAl34oZVddM4XpcZx9sHmGQ6nazqqhorOKdHf7/QbgHSHEPwEXcESburOBrwMVWLih4rMoFZ3ZuouSl98BYMifLmXIdZcAsLmqidd+LkMA1x7ZnxiL4WUTwh4V21k7a8zKfpO6lHKJEKI/MAzYIKVsuyLAe8ArgQwunLBYwn9dbH+jmrN0uVh7znW+/daE7nJLHlqyA7eEU4uyGZljrj+OqrUzaGeNeen2ckNKWS+l/H6fhI6Ucr2U0u+TvIQrdXV1oQ4h6Kjm/PHQY33bkz961rf97tpK1lc0kZUQw8WH9AlFaAFFtXYG7awxL+a5hxhgsrMjf+jSgaKSs7vFiavJBkDBBaeRctBwACobHTy7wvPbdc7kfiRYzXe1o1I7t6KdNWZFJ3WDVFcbXUzOPKjk/MPF83zbo+bP9W0/vmwnTS1uDitM5fAB5hzsoVI7t6KdNWZFJ3WDdLXwjZlRxXnHS29T8fFSAPrd80df+Tfba1lSvIe46CjmHNavq7dHPKq0c1u0s8asGFp6tRUhRA7Qbi0+KeUWv0YUpqh460oF56olK1h9/X0A9DnjWIac61m4xdbi4rGlnmFtF07oQ06SeRfDUKGd90U7a8yKoSt1IcRMIcROYDewqc1rYwBjCyvKysq6P8hkqOC84Z4nAYjtk82o+Tf4nF/+oZSyBgeDM+M5bZS5/xiq0M77op01ZsXo7ffH8KxtniSljGrzMl+voS5ISkrq/iCTYXbn5vIqar9fDcCo+/5EdGI8SUlJFNfY+N+qcgRw9eEFppgKdn+YvZ07QztrzIrR2+/pwJNSP5TRmAS3o4XPDj7Zt58943DA89zxka9LcEk4cUSm6cakazQac2P0Sv0Z4OJABhLuNDT0eIXZiMXMzmtu/BsAwmJhyvL/IbzroX+0sZpVpQ2kxkVzya/yQxli0DBzO3eFdtaYFaNX6pOAq4UQ84B2EwhLKad0/hZzkZubG+oQgo5ZnTfev4CSf3ungp17KQmFnuRd3+zkjU12AH5/aF+SYw+oH2nEYtZ23h/aWWNWjP7VWuB9KUtFRQUFBQWhDiOomNG5qbiEzQ/8C4ARf72aAb8/x1f37Ird1Da7ODgvielD0kMVYtAxYzt3h3bWmBVDSV1K+XygAwl3Wm/PqoTZnFtq6/ly0q99+4W/O9u3vb6ikffWVmIRcOXh/Uznvj9Ucm1FO2vMiuHJZ4QQFwshPhVCrPf+q9Qz9oyMjFCHEHTM5vzJ8ON824d/9qLvj5zL7ekcJ4GTR6QzID0+RBGGBrO1sxG0s8asGB2nfjMwD8+qbFd7/53rLVeCioqKUIcQdMzk3OhdHx2g/8VnkDxysG///XWVbKhsIisxhul5oYgutJipnY2inTVmxegz9UuBaVLKba0FQojFwJfAXYEILNxISUkJdQhBx0zOP//hLwD0PfsEiu7Zu7zqHlsLz67YDcDlk/qSY87p3feLmdrZKNpZY1aM3n5PBPb9mVcFKHOf0uVyhTqEoGMW55KX36X2hzUADPzDee3qnvluFw0OFxP6JnPkgDTTOB8I2lkNVHRWEaNJfRHwbyHEcCFEvBBiBPA8sDhwoYUXjY2NoQ4h6JjB2dVk55dr7wYgadhAkoYP9NWtLmtg8YZqYqIEV072dI4zg/OBop3VQEVnFTGa1K8E6oGfgAbgR6ARuCpAcYUdeXnqPWw1g/OPv7vFtz3xzcd82y635FHvgi1nHpxD39Q4wBzOB4p2VgMVnVXEUFKXUtZJKS8AEoA+QIKU8gIp5Z6ARhdGlJaWdn+QyYh059J3P/MtqTrizj9izdz7wPzdtZVsrrKRm2Tl3LF7/9hFunNP0M5qoKKzinTZUU4IMUBKWezdHrRPdVLrcCBVll6NiYkJdQhBJ5Kd3U4nP166d3BG4eyzfNs1TS089/3eznFx0Xt/20ayc0/RzmqgorOK7K/3+yog2bu9CZDAvrMXSECJldpSU1NDHULQiWTntTf93bc9fcOH7SbeWPDdLhodLiYWpDC5sL1jJDv3FO2sBio6q0iXt9+llMlttqOklJZ9ll1VaunVysrKUIcQdCLVec2ND7DjhTcBGPfcvcSk7F1y8pfSBj7aWE2MRfCHwzrOHBepzr1BO6uBis4qYnTymX90Uf6Qf8MJX1T8lRuJzo1bdrD92dcByDl+Crkz96435OkctwOAsw/OJT8ltsP7I9G5t2hnNVDRWUWM9n6/qIvy8/0UR9jjcDhCHULQiTRnZ0MjX03eO5/7uAXt50VauKaCLdV2cpOsnD2m8xWrIs3ZH2hnNVDRWUX2O6OcEOKS1uPabLcyCFDmfo7NZgt1CEEn0pw3P7R33aGxT9+JsOx9OlTV1MLz3s5xfzisH7HRnf+ejTRnf6Cd1UBFZxXpbprY1itxK+2vyiVQBlwYiKDCERXHeEaSs21nGcVP/AeAovv+RN7JR7erX/DtTppa3BxakMJhhV3fhowkZ3+hndVARWcV2e/tdynlUVLKo4B7W7e9r6OllOdKKb8JUpwhR8UxnpHk/MWE05EuF9nHHkH/C09vV/fz7no+2VSD1ds5bn9EkrO/0M5qoKKzihhdT903LZfwdBcWbercAYgr7LBaraEOIehEivP62/fOFDfgd79uV+d0Sx7xzhx3zphc+nTSOa4tkeLsT7SzGqjorCKGkroQIh94DJgC7LuOlRLD2pKTk7s/yGREgvP3v73eN2tc1tGHkXnEIe3q31pdwbYaO32Srfz64M47x7UlEpz9jXZWAxWdVcRo7/cnAQcwHc/c7+OBhcDlAYor7Kiqqgp1CEEn3J1b9tT5EjrA+Ofva1df2ejgxZWeznFzJvfD2kXnuLaEu3Mg0M5qoKKzihhN6pOBS6SUPwJSSvkTMBu4bv9v2z9CiH8JIcqFEL+0KcsQQnwkhNjo/Te9Td2NQohNQoj1Qojj2pRPEEKs8tb9Q+w7o4gfSE9P7/4gkxHuzhvuftK3PX39YqJi2t94evKbndha3EwuTGVigbExuuHuHAi0sxqo6KwiRpO6C3B6t/cIIbLxrNLWt5ef/xwwc5+yecAnUsqhwCfefYQQRcA5wCjvex4XQrTe+n8C+B0w1Pva95y9RsXhIOHsvGflGt+scYd//hIxqe1vLX5fUscXW/cQGx3Vbee4toSzc6DQzmqgorOKGE3qy4ETvNuLgf8CbwArevPhUsovgep9ik/Fs1Y73n9Pa1P+ipSyWUq5Fc989BOFEH2AFCnlMimlBF5o8x6/Ybfb/X3KsCdcnaWUfHPCpQAIi4XkEe3XG3K43L5lVX87Lo+cJOMdhMLVOZBoZzVQ0VlFDHWUwzNGvfUHwB+B64EkIBDTxOZKKXcDSCl3CyFyvOV9gbZD6Eq8ZS3e7X3L/YqKYzzD1XnHc2/4tg9b/EyH+td+LmdnXTP90+KYNTr7gM4drs6BRDurgYrOKmJ0SNueNts24I6ARdQ1nT0n72zluNbyDpSXlzN79myio6NxuVzMmjWLOXPmUFpaSmJiIhaLhbq6OrKzs6murkZKSXZ2NmVlZdhsNrKzs2loaCA3N5eKigqEEGRkZFBRUUFKSgoul4vGxkby8vIoLS0lJiaG1NRUKisrSU1NxeFwYLPZfPVWq5Xk5GSqqqpIT0/HZrNht9t99XFxccTHx1NTU0NmZib19fU4HA5ffXx8PFarldraWrKysqitraWlpcVX351TUpJnoZOunGw2G1lZWWHl5F6zlfU3PgBA3jknUJeWQO2OHT6nRhHHyz94xuP+/pAcSnftPKB2stvtvrgjpZ16+91rbm4mNTXVVE7dtZPD4SApKclUTt21k8PhIDExsUdOmshBeO5Yd1IhxO1GTiClvLVXAQgxAHhXSjnau78emOa9Su8DfC6lHC6EuNH7efd4j1sM/AUoBj6TUo7wlp/rff/v9/2sZcuWyREjRvQozrKyMnJzux8SZSbC0XlR3mTf9owtn2JJiPPtSym5efFmVpTUc/TgdOYdNeCAzx+OzoFGO6tBb5xXrlz5/fTp0w/p/khNqNnfM/UCgy9/s5C9089eCLzdpvwcIUSsEGIgng5x33pv1dcLISZ5e71f0OY9fiM+Pt7fpwx7ws25ZvlPvu1pK99ql9ABvirew4qSepKsFn5/aM+ewISbczDQzmqgorOKdHn7XUp5caA/XAjxH2AakCWEKAFuA+4FXhVCzAa2A2d541kthHgVWIOnJ/4cKaXLe6or8PSkjwc+8L78Sk1NDSkpKf4+bVgTTs5up5Plp14BQPqkscTl57Srb3S4eGLZTgAu+VU+6QkxPfqccHIOFtpZDVR0VhGjM8oN6qpOSrmlpx8upTy3i6rpXRx/F3BXJ+UrgNE9jcMImZmZgTx9WBJOzlv+8aJve/wL8zvUv/D9bqqaWhiRncAJI3oedzg5BwvtrAYqOquI0SFtm4CN3n83tdnfGKC4wo76+vpQhxB0wsW59L3P2TT/aQAGX3cJMSntO+9sqGzi7TUVRAm45ogConox91C4OAcT7awGKjqriNHe7+2SvxAiD8+t8q8CEVQ44nA4Qh1C0AkX5x9n3wRA5pGHMOT62e3qXG7JQ19txy1h1uhsBmcm9OqzwsU5mGhnNVDRWUWMXqm3Q0pZime8+j3+DSd8UXGMZzg416/b+3Rn5D3Xse8MwG+vqWBTlY2cpBgunNCn158XDs7BRjurgYrOKtKjpO5lONC7y6IIQsW1iEPt7Gqy88NFNwBgSUwgaUhhu/ryBgfPrfAs2HLl5ALiY3q/YGConUOBdlYDFZ1VxGhHua9oP6FLAp452A2NZTcDKg4HCbXzx8NmIJ2eAQ6Hvv14uzopJY8tK8HudHPEgDQm9Te2YEt3hNo5FGhnNVDRWUWMThO7YJ/9RuAnKaUyHeWsVuPzh5uFUDqXLvzUl9BH3HENKaOHtatfUlzLsm21JMREMecAFmzpDt3OaqCdNWbFaEe557s/ytzU1taSlpYW6jCCSqic19/+GFsf/zcAQ+ZexoDLzm5XX9/s5LGlOwDPmPTMxJ6NSe8M3c5qoJ01ZsXQM3UhhFUIcbt3jfNG7793CCHiun+3OcjKygp1CEEnFM7VS3/wJfTk0UMZfM0FHY55evkuqm1ORucmctJI/8ao21kNtLPGrBjtKPcEcDRwNfAr779Tgcf39yYzUVtbG+oQgk4onDd6x6MDHLboGYSlfee3H3bVs2hDFTFRgj8e2b9XY9I7Q7ezGmhnjVkx+kz9NGBwm9Xa1gghluOZhOaSgEQWZrS0tIQ6hKATbOfqb36k5psfsSQmMO37N4iKbv/1tDvdPLxkOwC/GZdH/zT/3yjS7awG2lljVoxeqZfScfhaPLDbv+GELyqO8Qyms8vWzA8XzwMg94SpxKR1nKP6+RW72FXnYEB6HL8+OKdDvT/Q7awG2lljVowm9ReBRUKIy4QQxwshfge8D7wghDi69RW4MEOPimM8g+UspeSjgUfRUlMHwPA//6HDMWvLG3njF89UsNdPKSTG0pspFrpGt7MaaGeNWTF6+711bfKb9im/3PsCzzj2Lhd+iXQSExNDHULQCZbzT1fc5tue+NbjxOa0X3jC4XTzty+2IYFfH5TDsOzAzXmk21kNtLPGrBgd0jYw0IGEOxZL72crizSC4WzfVU7pWx8DUHjpWWRMGtvhmBd/KGVHbTMFqbGcP773U8HuD93OaqCdNWbF8D1MIUS0EGKKEOJcIcSRQgijV/mmoK6uLtQhBJ1gOJe8/I5ve+Sd/9ehfn1FI6/9XIYArptSiDU6MLfdW9HtrAbaWWNWjE4TOwJ4B0/nuB1AAWAXQpwspVwbwPjChuzs7FCHEHQC7Vzyn3fZ9LdnAJjwnwc71Dc73cz/fBtuCWcelENRbuBvH+p2VgPtrDErRi97HgeeAgqklIdJKfsB/0ShcerV1dWhDiHoBNLZvruCNTc9AEDuidPImjqxwzHPrdjFjtpm+qfFcZEfVmAzgm5nNdDOGrNiNKmPBR6UUrZd1OUhb7kStFdXg0A5u1ucrLxwLm5bMznHT2HcM3cjotp/FX/eXe/r7f6nqf0Dftu9Fd3OaqCdNWbF6F/KXXhmkGvLkd5yJVDx1lWgnD8bcwp1P68HYNi833eob3K4+NuX25HAuWPzGJ4dvF67up3VQDtrzIrRpH4TsFAI8YoQ4j4hxCvAQjoOcTMtZWVloQ4h6PjbWbrdfH3MhbRUeyYmHHDFb0ga3nFgxRPflFBa72BIZjy/GZvr1xi6Q7ezGmhnjVkxlNSllAuBccAvQLL33wlSyrcDGFtYkZSUFOoQgo6/nbc8/Dz1v3hW680/cyYjbruywzFLtu5h8YZqrBbB3GmBm2SmK3Q7q4F21piV/fZ+F0IkALcAo4GVwD1SyuZgBKYxFy57Mxvv8yzWUnDB6Yya/6cOx1Q2Ovi7d273Syf2ZUB6fFBj1Gg0mkinu8ugR4GTgXXAmcDfAh5RmNLQ0BDqEIKOv5yl281HA47y7Y+8q+N4dLeU3P/FduqbXRzSL5lTi0KzTKRuZzXQzhqz0l1SPx44Vko517t9UuBDCk9yc4P7bDcc8Jdz6cJPfNtDb7qcqJiON4jeWFXOD7vqSY2L5vophQg/L6lqFN3OaqCdNWalu6SeKKXcDSCl3AGkBj6k8KSioiLUIQQdfzg7aupYf4dnOoMBl5/L4Ksv6HDM+opGnvnOM5Di2iP7k5EQ0+vP7Sm6ndVAO2vMSnczykULIY4CRBf7SCk/DVRw4USorhxDiT+cvxh/Gi6bnZQxIxh2yxUd6hsdLu76tBiXhNNGZXNYYWh/N+p2VgPtrDEr3SX1cuBfbfar9tk39cpsbcnIyAh1CEGnt87Vy37AZbMDMPqBeURFt/+6SSl56KvtvuFrl07M79Xn+QPdzmqgnTVmZb+336WUA6SUA/fzUiKhg5q3rnrjXPb+F3x7+hwAkouGkDJ6WIdj3ltXxRdb9xAfE8XNRw/AGuTha52h21kNtLPGrIT+r2iEkJKSEuoQgk5PnfesXMMPl9wIQMKgAg77YEGHYzZUNPHEshIArj68gL6pcT0P1I/odlYD7awxKzqpG8TlcoU6hKDTE2fpdvPNCZf69ie98yRRsdZ2x9TZndzxyVZa3JKTRmQxfUj43BbU7awG2lljVnRSN0hjY2OoQwg6B+osXS4W5x/h2z/885ewZqa1O8YtJfO/2EZZg4NhWQlcflhfv8TqL3Q7q4F21pgVndQNkpeXF+oQgs6BOm+4+5++7Qkv/Y3kER27XLz8Yxnf7qgjOdbCn6cPDIvn6G3R7awG2lljVsLrL2oYU1paGuoQgs6BONt27GbrY/8GoM8Zx5J9zOQOxyzbVssL3+9GADdMKyQ32drhmFCj21kNtLPGrOikbpCYmNBNiBIqDsR55UXzAIgv6MPBj97WoX57jZ37Pi8G4KJD+jCxIDznMdLtrAbaWWNWdFI3SGpqeCahQGLUuX7tZupXe1ZfO+gft3SY5KK+2cltH22hqcXN1IFpnDMmfKer1O2sBtpZY1Z0UjdIZWVlqEMIOkadi598BYD+F80i47Bx7eqcbsldnxazs66ZQRnxXDulf1jPbKXbWQ20s8ashG1SF0IUCyFWCSF+FEKs8JZlCCE+EkJs9P6b3ub4G4UQm4QQ64UQx/k7HhV/5RpxLn3nU3a+8h4AfU6f0a5OSsljS3ewcqdnoZa/zBhIfIwlILH6C93OaqCdNWYlbJO6l6OklGOllId49+cBn0gphwKfePcRQhQB5wCjgJnA40IIv2YPh8Phz9NFBN05163eyE+Xe56f5585k/RDx7Srf31VOe+tqyLGIvjrjEHkJccGLFZ/odtZDbSzxqyEe1Lfl1OB573bzwOntSl/RUrZLKXcCmwCJvrzg202mz9PFxHsz7lp2y6WTr8Q6XKRddQkRv/9pnb1S4r38PS3npXX5k4tpCg3MaCx+gvdzmqgnTVmpbsFXUKJBD4UQkjgSSnlU0Bum6VgdwshcrzH9gW+afPeEm9ZO8rLy5k9ezbR0dG4XC5mzZrFnDlzKC0tJTExEYvFQl1dHdnZ2VRXVyOlJDs7m7KyMuLi4qiqqqKhoYHc3FwqKioQQpCRkUFFRQUpKSm4XC4aGxvJy8ujtLSUmJgYUlNTqaysJDU1FYfDgc1m89VbrVaSk5OpqqoiPT0dm82G3W731cfFxREfH09NTQ2ZmZnU19fjcDh89fHx8VitVmpra8nKyqK2tpaWlhZffXdOSUlJAF06xcfHU1lZ2c4pOjqahlc/ZOt9T/v+u2b8+Xe0uF2UbtuJ1WqlxB7NPZ+VIIGzRiQzwFJHc3NCWDh1104JCQmUl5dHVDv19ruXmJhIWVmZqZy6a6fk5GR2795tKqfu2iklJYVdu3b1yEkTOQgpZahj6BQhRL6Ucpc3cX8EXAUslFKmtTmmRkqZLoR4DFgmpXzJW/4M8L6U8vW251y2bJkcMWJEj+LZtm0bhYWFPdWJSDpz/njoDJz1npmpLEkJjH/uPjKPmOCr31pt47p3N9LgcHHCiEyuObwgrDvG7YtuZzXQzgfGypUrv58+ffoh3R+pCTVhe/tdSrnL+2858Cae2+llQog+AN5/y72HlwAFbd7eD9jlz3is1vCbKCXQ7Ou85ZEXfQm9/yVncsyGD9sl9PIGBzcv2kyDw8XkwlSumhxZCR10O6uCdtaYlbBM6kKIRCFEcus2cCzwC7AQuNB72IXA297thcA5QohYIcRAYCjwrT9jSk5O9ufpIoK2zlsefYkNdz0BQNKIQRTdfS0iau/Xp6aphXkfbKKyqYXReYnceNQALFGRldBBt7MqaGeNWQnLpA7kAkuEED/hSc7vSSkXAfcCM4QQG4EZ3n2klKuBV4E1wCJgjpTSr0sSVVVV+fN0EUGrc8Wn37DhzscByD5mMkd8/lK74+rsTuZ9sImS2mYGZcTx1xmDiI0O16/W/lG5nVVCO2vMSlh2lJNSbgHGdFJeBUzv4j13AXcFKqb09PTuDzIZ6enpuJsdrJ4731c2/oX57Y5pdLi4adFmttbYKUiN5Z7jh5AcG5ZfK0Oo2s6qoZ01ZiUyL6dCgIrDQep3lfFh4TTsJZ6FII5a9W67W+6NDhc3L9rMhsom+iRbue+EIaTHR/b80iq2s3ZWAxWdVUQndYPY7fZQhxBUapb/xKppF/r2J73/NLHZGb79hmbPLfc15Y1kJ8Zw3wlDyEqM/I44qrUzaGdVUNFZRSL3PmmQUWkt4oZN21h+6hW+/dEP3kTa+FG+/Tq7kxsXbWJjpY3cJCvzTxwSEbPFGUGldm5FO6uBis4qoq/UDaLCWsRSSrY98z+WHHGur+yYzZ/Q7zcn+farm1qY+/5GNlbayE+x8sBJQ+ljkoQOarTzvmhnNVDRWUX0lbpB4uLiQh1CQHE2NvHT7/5MxSfLfGXDX3mA6MR43/6uumZu/GATu+sd9EuNZb5Jbrm3xezt3BnaWQ1UdFYRndQNEh8f3/1BEUr92s38cv291H6/mqg4K4OvvYRBV51PfX2975jNVU3ctGgzNTYnw7ISuPO4QaRFeKe4zjBzO3eFdlYDFZ1VRN9+N0hNTU2oQ/A7Ukq2PvZvvp5+IbXfr8aalc5hHzzD4KsvQAjhc15RUsd1726kxuZkXH4S808YYsqEDuZs5+7QzmqgorOK6Ct1g2RmZoY6BL8ipWTtLX9n+zP/A6DPGccy4rariM3Z65mZmcm7ayt5dOkO3BKmDkzjT9MKsVrM+1vQbO1sBO2sBio6q4hO6gapr683zYpFzkYba+b9jV2vfYCItnDQI3+mz2kz2s3T7nJLnlq+k0VbmwA4d0wuFx7Sh6gIm8v9QDFTOxtFO6uBis4qopO6QRwOR6hD8As1363iuzOvwt3swBIfx0H/uIW8k49ud0yd3cldnxbzw64moqMEfzyigGOHqfEr3yztfCBoZzVQ0VlFdFI3iBnGeJYu/JQff3cLANHJiUx863FSRg1td8ymyib++vFWyhocpMZa+PMxgzi4jzq/7s3QzgeKdlYDFZ1VxLwPR/1MpI/xdFTXsu6vj/j2p3zzWruELqXk/XWV/PGdDZQ1OBiencBNE5OUSugQ+e3cE7SzGqjorCL6St0gkTocREpJ5SfL+GnOX3HW1pM0fCCHf/oCwmLxHdPocPHQku18sWUPADOHZXLl5H7sqa4MVdghI1LbuTdoZzVQ0VlFdFI3iNUaWZOsOBtt/HLdPZS+9bGvzBIfx4SXH2yX0FeXNjD/i23srncQHxPFoa1+BgAADkRJREFU1YcXMH2IZ473SHP2B9pZDbSzxqzopG6Q2tpa0tLSQh2GIWw7dvPFoWeB2+0riy/M55D//J34vrkAOFxuXvx+N6+tKsctYXBmPDcfPYB+qXtnnYokZ3+hndVAO2vMik7qBsnKygp1CIYo/+hrfph9ky+hH/ruk6SOGUlUzN6mXl/RyINfbmdrjZ0oAWePyeX88Xkdxp9HirM/0c5qoJ01ZkUndYPU1taSmJgY6jC6pOKzb9jy8PPUfPOTr+yoVe+2Wy7V1uLiue938/bqCtwS8lNi+dPU/ozK7bwzXLg7BwLtrAbaWWNWdFI3SEtLS6hD6JLip/7Lulsf9u1nzzic0Q/e6EvoUkq+Kt7Dk9/spKKxhSgBZx2Uw2/H5xEfY+nqtGHtHCi0sxpoZ41Z0UndIOE4xtNla2bV1XdQ+s6ngGeq16J7ricmZe+V99ZqG48vK+Gn3Q0ADMtK4P+OLGBwZkK35w9H50CjndVAO2vMih6nbpBwGuMppaRs0Zcsm3mJL6EPv/VKxjz2F19Cr2ps4aEl27nizXX8tLuB5FgLVx9ewMOnDDOU0CG8nIOFdlYD7awxK/pK3SDh8iyq4uOlrL/rCRrWbgYgKj6WcQvuJnv6YQDUNzt57edy3vylnGaXJErAKUVZXDC+DylxB9bc4eIcTLSzGmhnjVnRSd0gFkvXz54DjdvRwq7/LeaXa+9uV5455VeMvOv/SBo6gDq7kzd+Keet1RU0tXh6vh8xII2LD+lDQVpcZ6ftllA6hwrtrAbaWWNWdFI3SF1dHenp6UH9zIaNxex89QNKXnqblpo6T6EQ9L/wdIbfdhWW+FgqGh28vHwn762rxOZN5uPyk7jokHxG5vTul3konEONdlYD7awxKzqpGyQ7Ozton9WwfivFT/+XkpcW+sqiU5LIP3MmQ669GGtWOhsrm3hzeTGfba7BJT3HTOibzG/H5TEqzz/ztQfTOVzQzmqgnTVmRSd1g1RXV5OQYKyDWU+QLhdNxTvZ/eZHbPrbM77y5KIhDLzqt+SddDROEcWXW/ew8Ov1rKvwrHMeJWDqoDTOPCiH4dn+fWYWaOdwRDurgXbWmBWd1A0ipQzIeZu27WLX64vZ8dwbNJdX+crj+uUx7KbLyT3pKDbXtfDYt7v5bHMNDQ4XAElWC8cOy+C0UdnkJccGJLZAOYcz2lkNtLPGrOikbhB/3rpq2VPH+jseo+Tf7yCsMUiHZ1IIa2Ya8QP6Mmze72kecxCfb67hs3c2sa3G7nvvkMx4Ti7K5qjB6cRFB3ZEooq367SzGmhnjVnRSd0gZWVlFBYW9vj9UkrKF3/FjhfepvLTZXvLHS1kHDGB/hecjmPyoXxdUs+/i2vZ8Ooa3zEpsRamD8ng2GEZhseY+4PeOkci2lkNtLPGrOikbpCkpJ53Ptv99ids+tszNG4sbleeNeNw3Jecxw/xWfxzRx0lb23w1cVFRzG5MJWjBqczvm8yMZbgzxPUG+dIRTurgXbWmBWd1ANI5efLKX7qVd+VeXRyIqmzz2H75CP5wWbh590N2De7gQrA85x8UmEqRwxIZXzflIDfXtdoNBqNudBJ3SANDQ1kZmbu9xiXrZnShZ9Qs/wnan9eR/0vGwGQUVHsOPVUPjrkKGqlBTY0+t5TmB7HoQUpTCxIpSg3kegoEVCPA8GIs9nQzmqgnTVmRSd1g+Tm5nZaLt1utj/3JhVffEvl4q/a1bksFtYdfAhLp59EfVoGSMhOjGFMnyTG901hXH4ymYkxwQi/R3TlbGa0sxpoZ41Z0UndAFJKysvL6d+/Pw5bMxs/WELp92tp/HENbC7GumeP71inJZqvZ5xMRV5fKvoV0ic/kyk5CYzKTeKgvCRyk60hNDkwKioqKCgoCHUYQUU7q4F21pgVndQNsOjZDxA33cnS/L7EV1UR2+wZYtaanutS09lVOJiSY48jc9xIJuQkMSwrnqFZCftdrzzcESJ8HgUEC+2sBtpZY1Z0UjeA7b8LSQDSdu0EYE9eHxrGjiV2TBHZo4cw5qBBnJKVaLqObRkZGaEOIehoZzXQzhqzYqqkLoSYCTwMWIAFUsp7/XHeoifvYP3rixk6fiQFIwpJystU4ldvRUWFcuNatbMaaGeNWTFNUhdCWIDHgBlACfCdEGKhlHLN/t/ZPSMGZJN7yYnKrXCUkpIS6hCCjnZWA+2sMStmul88EdgkpdwipXQArwCn+uvkLpfLX6eKGLSzGmhnNVDRWUXMlNT7Ajva7Jd4y/xCY2Nj9weZDO2sBtpZDVR0VhHT3H4HOnvI3W5ZovLycmbPnk10dDQul4tZs2YxZ84cSkv/v717i7HqquM4/v11YDqXzgwMjEChEdAqwRfBxhQ12jhNa+ulPvjAQ9OaSGJiTKw+GEhfqjHxEm8xJmpSNNhqb0iwaWqoook+ENBStVCKbQMIZZgrDgMZOjP492Ev8IQChw5n5pyzzu+T7Mzea1+yfmzO/HPW3rP3Cdrb22lqauLUqVP09PQwMjJCRNDT00N/fz8tLS0MDw9z+vRpFi1axODgIJLo7u5mcHCQzs5Ozp07x5kzZ1i8eDEnTpxg7ty5dHV1MTQ0RFdXFxMTE4yPj19Y39zcTEdHB8PDw8yfP5/x8XHOnj17YX1LSwutra2cPHmSBQsWMDY2xsTExIX1ra2tNDc3Mzo6ysKFCxkdHWVycvLC+nKZzj828nKZWltbGRoayipTufPU1tbGwMBAVpnKnaf29nb6+/uzylTuPHV0dNDX15dVpnLnqbOzk+PHj08rk9UP5fI6PknrgIci4s60vAkgIr55fptdu3bFqlWrpnX8I0eONNxNJs7cGJy5MVxL5r179z7f29t7S4W7ZDMgp+H3vwI3S1ohqRlYDzxdqYNv3769UoeqG87cGJy5MTRi5kaUTVGPiCngi8AO4ADwZETsr9Txt23bVqlD1Q1nbgzO3BgaMXMjyumaOhHxLPDsTBx7ampqJg5b05y5MThzY2jEzI0om2vqV2Pnzp2DwJHp7DsyMrKwu7t7qMJdqmnO3BicuTFcY+a39/b29lS0QzYjGqqom5mZ5Syba+pmZmaNzkXdzMwsEy7qV0HSxyQdlPSqpI3V7s90SbpJ0p8kHZC0X9KXUnu3pN9LeiX9nF+yz6aU+6CkO0va3yfpxbTuR6rxN9xIapL0gqRn0nLWmSXNk7RV0svpfK9rgMxfTv+v90l6TFJLbpkl/VzSgKR9JW0VyyjpeklPpPbdkpbPZj6rgIjwdIWJ4o1vrwErKV6h/g9gdbX7Nc0sS4C1ab4D+BewGvgOsDG1bwS+neZXp7zXAyvSv0NTWrcHWEfxJL/fAXdVO1+Z7F8Bfg08k5azzgxsATak+WZgXs6ZKR4JfQhoTctPAp/NLTPwYWAtsK+krWIZgS8AP03z64Enqp3Z01ub/E29vBl9Ucxsioi+iNib5sco/p5/KUWeLWmzLcCn0/w9wOMR8UZEHAJeBd4vaQnQGRG7ovj0/7Jkn5ojaRnwceDhkuZsM0vqpPjlvxkgIiYi4j9knDmZA7RKmgO0AcfJLHNE/BkYuai5khlLj7UV6K2lkQorz0W9vBl9UUy1pGG1NcBuYFFE9EFR+IG3pc0ul31pmr+4vVb9EPgq8N+StpwzrwQGgV+kSw4PS2on48wR8TrwXeDfQB8wGhHPkXHmEpXMeGGfKB7oNQosmLGeW8W5qJdX9kUx9UbSDcBvgAci4tSVNr1EW1yhveZI+gQwEBHPX+0ul2irq8wU31jXAj+JiDXAGYph2cup+8zpOvI9FMPMNwLtku690i6XaKurzFdhOhlzyt+QXNTLOwbcVLK8jGJYry5JmktR0H8VEeefG9mfhuRIPwdS++WyH0vzF7fXog8Cn5J0mOLSyUclPUremY8BxyJid1reSlHkc858O3AoIgYjYhLYBnyAvDOfV8mMF/ZJlzG6ePNwv9UwF/XyZvRFMbMpXRvbDByIiO+XrHoauD/N3w/8tqR9fbojdgVwM7AnDfGNSbo1HfO+kn1qSkRsiohlEbGc4tz9MSLuJe/MJ4Cjkt6dmnqBl8g4M8Ww+62S2lJfeynuGck583mVzFh6rM9QfF78Tb2eVPtOvXqYgLsp7hR/DXiw2v25hhwfohhK+yfw9zTdTXHNbCfwSvrZXbLPgyn3QUruAgZuAfaldT8mPZ2wlifgNv5/93vWmYH3An9L53o7ML8BMn8NeDn19xGKu76zygw8RnHPwCTFt+rPVTIj0AI8RXFT3R5gZbUze3prkx8Ta2ZmlgkPv5uZmWXCRd3MzCwTLupmZmaZcFE3MzPLhIu6mZlZJlzUzczMMuGibmZmlgkXdTMzs0y4qJuZmWXCRd3MzCwTLupmZmaZcFE3m0WS9ku6rdr9MLM8zal2B8xyIul0yWIb8AZwLi1/PiLeM/u9MrNG4be0mc0QSYeBDRHxh2r3xcwag4ffzWaRpMOSbk/zD0l6StKjksYkvSjpXZI2SRqQdFTSHSX7dknaLKlP0uuSviGpqXppzKzWuKibVdcngUeA+cALwA6Kz+VS4OvAz0q23QJMAe8E1gB3ABtms7NmVttc1M2q6y8RsSMipoCngB7gWxExCTwOLJc0T9Ii4C7ggYg4ExEDwA+A9VXruZnVHN8oZ1Zd/SXz48BQRJwrWQa4AbgRmAv0STq//XXA0dnopJnVBxd1s/pwlOJO+oXpW72Z2Zt4+N2sDkREH/Ac8D1JnZKuk/QOSR+pdt/MrHa4qJvVj/uAZuAl4CSwFVhS1R6ZWU3x36mbmZllwt/UzczMMuGibmZmlgkXdTMzs0y4qJuZmWXCRd3MzCwTLupmZmaZcFE3MzPLhIu6mZlZJlzUzczMMvE/9uNcr6d7MlYAAAAASUVORK5CYII=\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(de_times, de_pop, label='ODE')\n",
    "ax.plot(imp_times, imp_pop, label='Implicit IEBM')\n",
    "ax.set_ylabel('Population Size', labelpad=10)\n",
    "ax.set_xlabel('Time', labelpad=10)\n",
    "plt.legend(bbox_to_anchor=(1,1), loc='upper left', borderaxespad=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explicit IEBM\n",
    "The other way to make a logistically growing population is to specifically add conspecific competition for space and not allow individuals to overlap. Individuals are given a *radius* to explicitly simulate space and a maximum distance offspring can spawn from their parent (*offspring_dist_max*). The *radius* is modelled as a simple *StaticTrait* that doesn't change. We could model the *offspring_dist_max* also as a *StaticTrait*, but this provides a simple example to use another *Trait*, a *LinkedTrait*. The *offspring_dist_max* depends on the *radius* of the individual and a factor. The *LinkedTrait* take the *radius* value and multiplies it by the factor to give the *offpsring_dist_max* value. These *LinkedTraits* are real helpful in future experiments. \n",
    "\n",
    "The *BirthEvent* is now a *BirthDiffusionEvent*, which is coded to spawn offpsring near parents based on the *radius* and *offspring_dist_max* *Trait* values. One option is that offspring are born near their parents and can overlap in space; in the other option, offspring cannot overlap. When individuals cannot overlap, their positions are stored in a *r-tree*, which is an efficient spatial data structure that quickly finds near by neighbours. When producing an offspring that cannot overlap, the parent neighbours are found with the *r-tree* and available space is determined by the vertices of a Voronoi diagram of the parents and neighbours positions. \n",
    "\n",
    "For a given environment size, there are different values of *radius* and *offspring_dist_max* that recreate the dynamics of the logistic growth ODE and implicit IEBM. In this example, from trial and error, a *radius* of about 0.4 and an *offspring_dist_max* of 20 have very similar dynamics as the above models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from iebm.traits.linked_trait import LinkedTrait"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def explicit_run(radius, odm_factor):\n",
    "    # we can recreate the exponential growth IEBM and add an implicit_capacity to the population \n",
    "    # First lets create an empty 2D population. Each individual contains only a unique ID and x-y positions\n",
    "    log_pop = Population2D(name='log_pop',\n",
    "                           init_size=init_size, \n",
    "                           implicit_capacity=k, # limit the population growth\n",
    "                           xdim=xdim, ydim=ydim)\n",
    "\n",
    "    trait_list = [(StaticTrait, {'name' : 'radius', \n",
    "                                 'value' : radius}), \n",
    "                  (LinkedTrait, {'name' : 'offspring_dist_max', \n",
    "                                 'link_trait' : 'radius', \n",
    "                                 'link_func' : lambda x: x*odm_factor}),\n",
    "                  (StaticTrait, {'name' : 'birth_rate', \n",
    "                                 'value' : birth_rate}), \n",
    "                  (StaticTrait, {'name' : 'death_rate', \n",
    "                                 'value' : death_rate})]\n",
    "    log_pop.add_traits(trait_list)\n",
    "\n",
    "    event_list = [(BirthDiffusionEvent, {'name' : 'birth',\n",
    "                                         'is_primary' : True,\n",
    "                                         'allow_overlap' : False,\n",
    "                                         'current_time' : 0}),\n",
    "                  (DeathEvent, {'name' : 'death', \n",
    "                                'is_primary' : True,\n",
    "                                'current_time' : 0})]\n",
    "\n",
    "    # add events\n",
    "    log_pop.add_events(event_list)    \n",
    "    \n",
    "    # add a trigger\n",
    "    log_pop.event_dict['birth'].triggers = log_pop.event_dict['birth'].insert_rtree\n",
    "    log_pop.event_dict['death'].triggers = log_pop.event_dict['birth'].remove_rtree\n",
    "    \n",
    "    # Add the exp_pop to a dictionary \n",
    "    pop_dict = {str(log_pop) : log_pop}\n",
    "\n",
    "    # Create the simulation with the population dictionary\n",
    "    sim = Simulation(pop_dict)\n",
    "\n",
    "    # Run the simulation\n",
    "    sim.run(time)\n",
    "\n",
    "    # Store results\n",
    "    iebm_pop = sim.population_history[str(log_pop)]\n",
    "    iebm_times = sim.time_history\n",
    "    \n",
    "    return iebm_times, iebm_pop\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6min 49s, sys: 392 ms, total: 6min 50s\n",
      "Wall time: 1min 1s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "iebm_times, iebm_pop = explicit_run(0.4, 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAEMCAYAAADd1S/AAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOydd3hUVfrHPydl0hPSC72HgLQoICpFEJC1AdZ1QV3UBbHsz4JY1rpi23XtrIprZy1rwwaCAlJVmiIgPUAI6W2STMnMnN8fMxkypA3JZGYy93yeZ57ce865d95v3kneOeeec14hpUShUCgUCkXHJ8jXBigUCoVCofAMKqgrFAqFQhEgqKCuUCgUCkWAoIK6QqFQKBQBggrqCoVCoVAECCqoKxQKhUIRIIT42gBvsnr1ahkWFuZrMxQKhaJDUVNTU5yWlvbPkJCQAajOoK+xWSyW3Uaj8ans7GzbyZWaCuphYWFkZma26trDhw/TvXt3D1vk3yjN2kBp1gZt0bx3715TYmLi5WFhYQ2CiML7mEym00pKSgCeOLlOfeNyk9DQUF+b4HWUZm2gNGuDtmjW6XThKqD7D2FhYTbHqEkDVFB3k7i4OF+b4HWUZm2gNGuDtmgWQnjQEoWHaDR+q6DuJsXFxb42wesozdpAadYGHV1zTk5O6CWXXNJ78ODBgwYNGjTo5ptv7moymcTy5ctjunbtOnTkyJFZQ4YMGXTuuef2//jjj53fYO67776Mfv36DR4xYkRW3aukpCTYl1raE009U28L6pu9NlCatYHS3LGw2WxcffXVfWbOnFk4Z86cAxaLhRtuuKH7ggULOk+ZMqVi6NChVV988cV+gJ9++ini2muv7RMREZEzdepUPcC1115bcO+99xb4VoV3UD11NzGbzb42wesozdpAadYGHVnzsmXLYsLCwmxz5swpAQgJCeGZZ545+umnnybV1NS4xLERI0YYbrnllrzFixen+MZa36J66m5iMBh8bYLXUZq1gdKsDTyledLibdkeudFJfHv9sC1N1e3atSti4MCBNfXL4uPjbampqeYDBw40WKecnZ1d8+qrr6bVnb/55pupn332WSJAbGysZeXKlXs9abs/oYK6m6SlpbXcKMBQmrWB0uwbbOZaqg8eJTQ+lvDUpBbbSymxSbDYpMvLWveSdcdgcRzbbBKrBJuUhIR23OF3KSVCiAZ5wh3ljbavj5aG31VQd5P8/HzNrWtVmrWB0tw6pJQYLTYMtTaqTRaMtTaMNomh1mo/tthfVWUGKnNKMZUbCA23EvfFF5iiOxG1fz+xObsRgCk2gZqULpT16k91fCI6m47gYzvZefZEasIjqNCFYEDQIKoB6UcOUpqchikisll7e8eFsOiyTm3SDM33qNuLrKwswzfffBNfv6ysrCyosLBQ16tXL9Pq1atd2m/dujWyZ8+eRm/a6C+ooO4mOp3O1yZ4HaVZG2hdc63VRoXRUu9lRV9tpNIqqDJZ0JusVJmsVJmtVJstVBstRB46RHFcCmfmlZ/S+9YCNVnjASjtNqzRNsLRrjbzLHoXGwD7sHmtuYau29dD2TH0SWlU9D6NkE6dgSg4rifu4EaCTQYMsTGU9+xNea+uiCCBLVKHCAshSddxl5mff/75+sceeyzotddeS7zhhhtKLBYLd955Z9eLL764ODIy0kXY5s2bI1588cWMp59+OsdH5voUFdTdJCYmxtcmeB2lWRsEquZaq42iajPF1RaKq82U7juC7cPPyJl6Acd10VSYy6iqqKbSKojSV3DpGy8QX1IIgA5IcrzqsISEYIiMpjatJ3lnX0jf6lML6O4QUl2JJSq20bpQXST5I8470fak+opeg5zHYWZI/T3feR6du58e6z6Dgys9aq+3CAoK4r333tt/2223dX/hhRfSbTYbY8aMqXjiiSeO/fDDD9Hbt2+PHjlyZJbRaAxKSEioffjhh4/UzXwH12fqAO+9997+vn37dtyZg80gTn72EMhs3LhRqm1i3Udp1gb+prnm8DHC01M48vanVO06QMbl59NpWBZBYa4jClJKympqOVZazbGFL2M6kkeNLgxZXIpRBJOxd7dH7ZJCkDt2mkvwrCNly/cYkjJABGELDiHm2AEiio5hTEglsjCXiJLjBEWE0/ehW4kdM4q4zklUVZnZtvEwleVG+g1Kpc+AVEJ1J5ZPG0sqOPb+V9QcOsLmwnBMhFKd3h0Z3Lq+WFCtidufvrhV1x4+fLike/fuOa26WNEuFBYW7u7fv//Mk8tVT91N4uPjW24UYCjN2qA9NZes24zhyHGq9hzCUmNg4JN3IYJcV9LmvPoBvz/wXJP3yF3yhcu5LTiYIKu1Qbswx8tdpBCU9z6NmCN7sURGE1RrRlddSXT/nmRcOpmCr9ZQuusQFT0HYkhMoyzTddJ3t8ocuuX+yuDn7sdmzmTDxGuxmU50/nreMpP+9821v5fN1kB3p4QQxv+h0Z0+AQhPjKP3vD8CcFqdzVJyeH8JUkq6904kKLjxVcllm3dw5J2lhKUm8mVVD/vvJzbqFH47io6KCupuYjAYiI1tfFgsUFGatYEnNUspqfx1D7/e/AjV+3Ia1Oe+83mb36OxgO5k+GBqfz9EZfdMrKFhyOAQUrb/AED4gvlEDs6iJLeUXbtLm32PvvGp7BueAsMb1l38p2H0zUptUD7p8Oom73dyQG8tQgh69G15pnz86acRf7r9q0AmIG2SwqJCj9ig8G9UUHcTo1F7EymVZm3QVs3SZuPoO5+z6+6nPWLPzqEjKU1JY+fwUdRExaAzGRm2aytZu7YRe+QwwTX25coyPZ3g0WdREZ1M8ukDSOqdztYNhykcWulyv6JhY+wHBcCK/W7ZsG+X6+qnlIxYevRJZPAZXemU2Pwsc39EBAlNfra1iArqbuIP61q9jdKsDVqrWUrJ9wP/QG1p0xPGhrzyKHFD+qNLTuDQy0so2fQL1Z3iqf3qxISt34afybfT/+Q8T43W0TcpkiuSIuibFEmvhAgSIkc768uL9Sx+Zr3L+xzeVgzbTm1vc11YMFNmnEZNlZn0rnHowkN47+VNGA21Lu1ufWgiOl3H/1epxc+2Fun4n1QvodbyagOl2T2OffQNO255tEF5xuVT6XnTH7EZjMQNy0JKycFSAz/ureTnzDHsTsjGJoEzpwEQEiTonxzJFalRDEiNYkByFPGRDVOE5uwrZsN3+8k70vyM84jIUHplJjNpmn0yW2lRNds2HiYpNYaBwzMIC7ffu6nJgTf/bcIp/R46Elr8bGsRFdTdJDw83NcmeB2lWRucqubCb9c1COiTjqwhSGcPmFab5Nf8KtZvOMqGwxUUV5/o+QYJyEqJYljnGIakRzMgJYqwkIbPm0uLqnjz+fXYrM2vzrnhrjHExduHw6VNgnBNE5qcFuMM8PVRflYEKiqhi5tERET42gSvozRrA3c020xmjrzxMcvSRrN11nxn+fB3nmZK/gYIDeHX43qeXXeEK97bwd1f72fprmKKq2tJiAzh/P6JPDCxJx/PHMyzF/Xjmux0hmbENBrQ846U859/rWs2oN/ywATuXDjFGdDB/tzY3bzfys8dj4yMjMZ36zlFli9fHnPhhRf2Afjoo4/iHnnkkWafS4wdOzYTYN++fbo33ngjobE2+/bt02VnZw+su3/Xrl2H1k/1+tVXX8UAJCUlZdeVjRo1asDq1auj6q5PSEjIXrBgQUbdPQsKCkKSk5OHz507t9up6FM9dTcpKyvT3KxopVkbtKS5trKK7/pNalA+6fBqjhttvP5zHt/tL3XpkXeJC+PsHp04q0ccfZMiCWom2NbWWinIrcBisbH809/Ql5+Y0NWjbxLxSZFkDetMehfP7V2u/KwAuOyyyyqAiubarFmz5neAgwcPhn3yyScJ1113XfNLJ4D6qWDro9PpbD/99NMugKVLl8Y++uijXcaNG7cHICMjw7Rq1apOQB7A+++/H9+arW5VUHeTxMTElhsFGEqzNmhJ88kBPSw9BesHi7lrxWF25Fc5y1OjdYzrHc/4XvH0TAhv0GuWNskn72wFKanSm+g3MJX1K5uejX7FDSPo2rPRjlGbUX7uuCxfvjzmqaeeykhMTKzds2dP5HnnnVeWlZVleP3111NNJpN49913D2RmZppmzZrVIywszLZ///6I0tLS0AceeODojBkzXAL4K6+8krh9+/aoRYsWHTl27FjILbfc0j03NzcM4Kmnnjo8bty46oyMjGF5eXnb/v73v3fOyckJHzFiRNaMGTOK77777jatEaysrAyOjY211J2HhYXJXr16GTZs2BA5evTomqVLlyb84Q9/KM3Pzz+lfZxVUHcTvV5PdHS0r83wKkqzNjhZs7TZ2P23Zzny+v9c2qXNvZptf5jG17+XUPHDUQDCQoIY16sTk/olMjA1qtEeeXlpDYv/8UOD8qLj+gZlAD37JTH9mmy3h9Jbg/Jz61mWNrpdUq9Oyd/gdqKYffv2Rbz55psHk5KSLEOGDDlNp9MVr1+/fvfTTz+d8tJLL6W88MILRwGOHTsWtnLlyj179uwJmzZtWv/zzz9/R1P3vP3227uNGjVKP3/+/AMWi4XKysrg+vX333//sRdffDG1sR74yWzfvj16xIgRWXXnb7/99oHMzEyT2WwOGjFiRJbZbBYlJSWhH374oUsK2OnTp5d++OGHCRkZGbXBwcEyLS2tVgX1dsJsDshtgptFadYG9TU3NdRedMEfeK7LaKzb7eu3eyVEcFFWEuN6xROpC27QHqC4QM/Kz3eRm1PW7PsHBQn++vB52Gz2Z+ghoY3fz5No3c8dnaysrOquXbvWAnTp0sU0ceLECoBBgwYZNmzY4ExmcOGFF5YGBweTlZVl6ty5s+m3335rcrbgTz/9FPPmm28eAggJCSEhIaGZXY6ax53h9zVr1kTNmzev5+bNm3fW1V900UWVTz/9dOfk5OTaCy+8sMVh/sZQQd1NtLjGU2nWBqG7clg28qom65dN/xO7hp9JEDCmZycuGZjMwNSoZnvSm1YdYN2KfQ3Kb5w/lthO9glb0iYRQSfuEdT+sdyJFv3sKc2n0qNuL0JDQ52zKIUQhIeHS7AnfrFYLKJ+XX3ac/TnVBk7dmx1RUVFSH5+vjMOh4WFyaysrJo33ngjbcOGDb99+umnp5wrVwV1N9HiGk+lObCxGk3se/I1chYtaVB3eNp0Ps62r9kOCxZc3D+R6YNSSI9tend1aZP8+MNB9vyaT1H+iaH1qZcNJmtYRoP29QO6t9GSn+vQouYvv/wy/oYbbijZt29f2LFjx8IGDhxoXLt2baPPIEaOHKl/4YUXkufPn19osVjQ6/VB8fHxzrSusbGx1pqaGo999dyxY0e4zWYjJSXFcujQIecQ+2233ZY/evRofUpKSqtGClRQd5OOvhykNSjNgYmlqprSDdtclqbVkXvppXySNRqLLozwkCAuzkpixmkpdIpouCFMfWprrTz34AqXsthO4cy8eTQRkf6Xr10Lfj4ZLWru2bOnaeLEif1LS0tDFy5ceDgyMrLJdZL//Oc/j8ybN6/76aefnhQUFMTTTz99eOzYsdV19dnZ2Ybg4GDpzkS5k5+p33bbbcevvvrqsrpn6vXeMyckxDUMDxs2zDhs2LBW7+mrUq+6SXl5OZ06nfJISIdGaQ4sds5/mqNvf9qgPGn2Zfz4h0v5cncxVgm6YMElA5O51I1gLqVk5ee7+OWno86ykNBg+p+WxrkXZDp3cPM3AtnPTdEWzR0x9eqsWbN6TJ48ueLqq69uflJHB0WlXm0jFRUVmvsnoDQHDkfe+rTRgG596mEek8lU7yomSMDkfgnMyk4nOar53rXNJqmqNLLk35uoqjQ5y7PP7sH4qa374uxNAtXPzaFFzVpEBXU3SUpqOd1hoKE0Bwb63w+6ZFAb/OIDVI44g+c3F7K/1AhIhqRHc9OZXeiZ0PQQba3Zys9rD7FvZ4HLM3OAqJgwrrnlLCKj/W+ovTEC0c8toTXNb7/9do6vbfAFPg3qQoiuwNtAGmADXpVSPieESAA+AHoAOcDlUsoyxzX3ALMBK3CrlHK5ozwbeBOIAL4GbpMefLZQUVFBVFSUp27XIVCaOy7SZmNlv0lYq2qcZZ1OH8TgT17m3W0FfLLiCDYJyVGhXNYvkouH92xyZnBBXiXvvLih2febe894j9rf3gSKn08FLWrWIr7uqVuAO6SUW4UQMcAWIcQK4FrgOynlE0KIBcAC4G4hRBZwJTAQyABWCiH6SSmtwCLgRmAT9qA+BfjGU4bW1ta23CjAUJo7JlaDiRU9GwbZsEcXMOfTPeTrzQQJmDEomVnZ6RTm5Tbc/U1Kjh4q5asPfqVab2pwr4TkKP5005nUVJuJiml6Rry/Egh+PlW0qFmL+DSoSymPA8cdx3ohxG6gM3AxMM7R7C1gNXC3o/x9KaUJOCSE2A+MEELkALFSyo0AQoi3gUvwYFBX61q1QUfWbDWYOP75Sn5/4DmX8uCoCA6//AL/2lKBBHolhPN/53Sjf7K911Zfs7RJ1n+3n02rDjS4f0hoELc8MJHg4BNJWHRhvu4XtI6O7OfWokXNWsRv/iKFED2AYcCPQKoj4COlPC6ESHE064y9J15HrqOs1nF8crnH0OIaT6W547D7gec4/OoHLmVDFz+G8axRLPw+h5yjRoIEXDUklauHpRFaLzDXaf7fm5vJ2Vvc4N7TZg2nV/9kv9q4o610VD+3BS1q1iJ+EdSFENHAx8BfpZSVzfzzaKxCNlPuQmFhIbNnzyYkJASr1cr06dOZN28e+fn5REVFERwcTGVlJcnJyZSWliKlJDk5mYKCAmw2GyUlJVRVVZGamkpRURFCCBISEigqKiI2Nhar1Up1dTVpaWnk5+cTGhpKXFwcxcXFxMXFYTabMRgMznqdTkdMTAwlJSXEx8djMBgwGo3O+vDwcCIiIigrKyMxMRG9Xo/ZbHbWR0REoNPpqKioICkpiYqKCmpra531LWmq2we6KU1SSoqLiwNKU0t+qvucdCRN+fe/QNkK12fe3f+1gB+79ee1T3/HbIO0qGCuzYrirMxE8nKPOjXt2JrD+mW5wO4Gf0CTZmQSlyRITIwkPz/fr/zU1s9ecHAwx48fDyhNLX32goODycvLa5UmfyApKSm7V69ehrrzCy64oPSBBx7IP9X71F/q9uc//7n7rbfeWjB06NBG14V/9NFHcbt374544IEH8j/44INO/fv3NzbW9r777suIioqy3nvvvQWzZs3qsWXLlpioqCgrQEREhG3NmjW/v/LKK4lPPPFEl+Tk5FqLxSJ69eplePPNN3Oio6Nt9913X8aiRYvS169f/9uAAQNMAE8++WTKk08+2fXLL7/cPXr06JqT37MpfL5OXQgRCnwJLJdSPuMo2wOMc/TS04HVUsr+jklySCkfd7RbDjyEfTLdKillpqP8Ksf1f6n/Xm1Zp15WVkZ8fHyrru2oKM3+jc1kZs/fX+bwax86y0Z9vZiw0/rzwvqjrNxvX557Xt8Ebh7dhYiT9lRf/c3vbF6b41IWHBLEnAXj/HLDGE/SkfzsKdqi2R/WqddlS2vrfVq7fr25604O6o21q58RDuCqq67qOXbs2Mo5c+aU3HfffRnffvttp6lTp5Y9/PDDx8Gex12v1we/+OKLhxoL6k2tUw86ucCbCHuX/HVgd11Ad7AUuMZxfA3web3yK4UQYUKInkBf4CfHUL1eCDHKcc9Z9a7xCJWVlZ68XYdAafZPbOZa9L8f5Nvu41wC+uRjazH06c1fl+5l5f4ywkKCuHNMN+4a290loEsp2bUtr0FAnzZzOH99+LyAD+jQMfzsaQJRc2lpafCQIUMG/fbbb2EAV155Zc8XX3wxCexfAm677bYuo0aNGjBp0qR+9fdYr2PChAn9N2zYEAnw+eefx44aNWrAiBEjsiZNmtQP7IF47ty53VavXh21Zs2aTgsXLuwyYsSIrN9//71Ns0Nra2sxGAxB8fHxzq1gJ06cWL5ixYpOAHv27NHFxMRY4uPjLU3fpXF8Pfx+FjAT2CGE2O4ouxd4AvhQCDEbOAJcBiCl3CmE+BDYhX3m/DzHzHeAuZxY0vYNHpwkB5CcnOzJ23UIlGb/QlqtfNvzXKS54SzmsVs+ZevxahauykFvstIlLowHJ/ake7zruvMjB0r48PWfXcrm3neO5pY6+bOf2wtPaf7HvcvaJfXqnQunNJso5uQtVm+++ebjs2bNKnvssceO3HTTTT1nz55dUFlZGXLzzTcXAxiNxqAhQ4bUPPfcc7l/+9vf0h9++OGMul7yyeTn54fMnz+/xxdffPF7v379zEVFRS7DWuPGjaseO3Zsubs9/IULF3Z57rnn0gH69OljWLJkySGAZcuWxY8YMSK6pKQktFu3bsYZM2aU110TExNjTUtLM2/dujV86dKlnS666KKyDz744JQ3F/D17Pd1NP48HGBCE9c8BjzWSPlmYJDnrHOltLSUyMjI9rq9X6I0+xe/zH2oQUAf9K976XLVBXz6WyGv/HgAm4SRXWNZML4HUSelRN2/q4DP3nUdvbxx/ljKyko0F9T92c/tRUfXXD9taX0uuOCCys8++yz+wQcf7L569WpnGtOgoCBmzpxZCvCnP/2pZObMmX2auve6deuisrOz9f369TMDJCcntzrtKsC9996b21jwnzJlStmiRYuO2Gw25s6d223hwoVp9ecFTJs2rfSDDz5IWLt2bdxXX321p8MF9Y6Er+ce+AKl2T+QViu/3vwI+Uu/c5ZlPX4H6TMmExQdxUsbcvl8VxEAVw9LY+bwNIJOmmz6+y/H+fKDX5znF141lP6n2Zc4Vej9T3N7449+bm88pbmlHrW3sVqtHDhwIDwsLMxWXFwc0qNHj0YX5De3ekNK6dXVHUFBQZx//vnlixcvTqlfPn369PKFCxd2ycrKqqmfIe6U7u0ZEwMfNVynDfxNs7XGyPLO53D8U3sGtITRw5mSv4Fu183AEhHBwysP8vmuIkKDBAvGdeea7PQGAb1ab3IJ6Lc/OskZ0MH/NHsDpTlweOqpp1J79+5tfOmllw7eeuutPcxmswCw2WwsWbIkHmDJkiWJ2dnZ+qbucc4551Rv3rw5Zu/evTqAk4ffAaKjo616vd5jMXPDhg0x3bp1c9nZKSoqSi5YsCD3rrvuOt7a+6qeupsUFBRobo2n0uxbpM3GyszJznMREkz2u/8AoNJo4b7lB9hTVENMWDAPndeL09IaLj+qqTaz6PFVzvPZd5xDULDr/yV/0uwtlOaOx8nP1MeMGVNx3XXXFX/44YdJq1ev3h0XF2dbunSp/qGHHkpfuHBhXnh4uG337t0RZ555Zlp0dLT1nXfeOdjUvdPS0ixPPvlkzsyZM/tIKUlISKhdtmzZvvptLr300tI777yzxxtvvJH61ltvHcjMzGy41aKD+s/UAdauXbsbTjxTl1KSmppqfvXVV3NOvvaaa65pU1Y5ny9p8yZtWdJWUlJCYmKihy3yb5Rm31G+dSebpt7gPO933xx63TILgOJqM/d8c4DD5UZSo3U8fn5vusSFO9vu3HqMb/63o8E9p8wYxKDsLg3K/UWzN1GaTw1/WNJ2qnhqCZy/olKvKhQdAJvFwrH3v2LnnU86y/rMv8EZ0I9VGFnwzQEKqsx0jw/niSl9SIyy5yyXUrJp1QHWr9zf4L7TZg2nd2ZKg3KFQhFYqKDuJlVVVZr7Zq80e5eanFx+GHW5S9kZHz1P4jmnA3Ck3Mj8r/dRWmNhQEokj07qTWy4/U+4vLSGxf/4ocE9R43rxYhxvdDpmv5TV37WBlrTHMi99OZQQd1NUlNTfW2C11GavUfVvhzWnfNHl7Ls/z7jDOiHywzM/3o/ZQYLQ9KjeWRSL+eGMq89vYaKMufumfTOTOaiq4e5JF5pDuVnbaBFzVpEzX53k6KiIl+b4HWU5vbHZq5l+41/cwnoZ36zmEm5P5A8fhQAh0oN3PWVPaAPy4jh0cm9iQgNxmyysOjxVS4B/ayJfZk2K9vtgA7Kz1qhLZq1NPeqA9HokjfVU3eTQMpQ5S5Kc/vzbbexLudnLv8PcUNOTObMrTBy99f7KTdayO4cw0Pn9SIsJIiKshpee9p1uP3/HplEcMipf09XftYGbdFsNpuNJpMpKCwsrFVrpxWexWQyBVksloZZmFBB3W0SEhJ8bYLXUZrbl933/8t5HNmzC+es+y8i+MTy2Hy9ifmOgD4sI4aHz+uFLiSIouN63nphvbPd6Wf3YNzU1q3qAOVnrdAWzVVVVfklJSX/CwkJGYAa4fU1NovFsttoND7VWKUK6m5SVFTUodd4tgaluf3Q7z7A4cUfOc/HbPzQpb6kupYF3+ynuLqWgalRPHReT3QhQWz8/gDrV55YPjtnwTiiY8NpC8rP2qCtmgcOHPiEB81RtBMqqLtJbGysr03wOkpz+2CpNrB+vH15acblUxn8/P0u9XqThXuW7Sev0kzfpAj+Prk34SFBvLzwe2qqzM52V/1lZJsDOig/awUtatYiKqi7idXapv39OyRKc/uwdvQVzuPe/3etS53ZYuOhFYfIKTPSNS6MhVP6YKo0suiZtc42cfERXH/nGI89F1Z+1gZa1KxF1LMRN6murva1CV5HafY8eZ98i6mgGIDUqWOJ6nlihzerTfLE6sPsyK8iMTKUx8/vQ6jVxn/qBXTAowEdlJ+1ghY1axHVU3eTtLS0lhsFGEqzZ6k+lMuvNz0EQMrksxn2n8eddVJK/r0pl3U55UTpgrlzSDJvL/ze5fpxU/tz+tk9PW6X8rM20KJmLaJ66m6Sn5/fcqMAQ2n2LOvPPbFN89DFC13qPttZxOe7igkNEszpHs2q91w3w5o8Y1C7BHRQftYKWtSsRVRP3U1CQ0N9bYLXUZo9R83hY9gM9qROIz57maDQE396Px6p4JUfjwFwY794dnx9Yvlp5+6duOovo9rFpjqUn7WBFjVrERXU3SQuLs7XJngdpdkz6Hcf4OfLbwOg04jBJIwa6qw7VGpg4aocbBIuz4hiX72A7onlau6g/KwNtKhZi6jhdzcpLi72tQleR2luO1x/a+cAACAASURBVNYaIz9eNAdzUSkAA5+401lXbqjlgW8PYjBbmXSwgPJ1J9I9/+Vu7wR0UH7WClrUrEVUT91NtPgtV2luO/lfrcKit886PmvVO8QM6A2AxSZ57PscCvQmzi6rcrnm9kcnEXQKe7e3FeVnbaBFzVpEBXU3MZvNLTcKMJTmtiGtVnL+/T4A/R+42RnQAV778Ri/5Ok5TV9DZHkNACGhQcy9Z7xXAzooP2sFLWrWIiqou4nBYGi5UYChNLeNvI+/Rb9zH+EZKXS7boaz/Nu9JXy3JY9JeaXOsikzBjEou0tjt2l3lJ+1gRY1axH1TN1NtLjGU2luPVajiR23PgpA9xuvIDgiDIADJTU8v/YII+oF9IHDO/ssoIPys1bQomYtooK6m2hxjafS3DpqK/Ss6DHeed5t1jQAqkwWHll5iGG5Jc4/vD5ZKUyeNrDN79kWlJ+1gRY1axE1/O4mOp3O1yZ4HaX51LHoq/mu/2Tned975xAcGY6Ukqd/OILteCVxJgsAF/9pGH2zUtv0fp5A+VkbaFGzFlFB3U1iYmJ8bYLXUZpPnc1X3+E8Hrr4MdIusPfY/7ejkJ8PlnFuYYWz3h8COig/awUtatYibg2/Czs3CCG+F0L86igbI4S4vH3N8x9KSkp8bYLXUZpPjcOv/4/yn34FIC57oDOg7yqo5vWf88g+XuZse/ujk9pmqAdRftYGWtSsRdx9pv4IMBt4FejmKMsF7m4Po/yR+Ph4X5vgdZRm9ylcvpbd9z3jPB/5+SLAnhv98VU5TDxQQJzZPux+3sVZXl+21hzKz9pAi5q1iLv/Wa4FLpBSvg9IR9khoFd7GOWPaHE5iNLsPvueWuw8Pve3rwgKCUFKyb/WHsVWUOmsGzaqG0NGdmvsFj5D+VkbaFGzFnE3qAcDddte1QX16HplAY/RaPS1CV5HaXaPou82ot+5DxESzIQ9y9El2XtEX/1ewrpDZQwoOfFncu6FAzxmq6dQftYGWtSsRdwN6l8DzwghwsD+jB14FPiivQzzN7S4xlNpbhmr0cQWx+S4rjMvITTOPhkpp8zAvzfl0qu8mnCLlZCQIP768HnY/3T8C+VnbaBFzVrE3aB+O5ABVABx2Hvo3YEF7WSX36HFNZ5Kc8vUX4/e5eoLATBbbTy5+jCjDxTQp8y+7/uZ5/YmJDTYc4Z6EOVnbaBFzVrErSVtUspK4BIhRAr2YH5USqmpT0h4uHcyZvkTSnPzbJh0nfM45fwxxA7qB8DbW45zrKCK3jbprD9jjP9OP1F+1gZa1KxF3F3S9h8hxGQpZaGU8ue6gC6EeLl9zfMfIiIifG2C11Gam8aQm0/lr3sAe4704W88AcAveXo++rWQMUdOpLmce894goL8b9i9DuVnbaBFzVrE3eH3PwGvCyHubKRcE5SVlbXcKMBQmpvm+GcrnccjP7N/t60yWXhqzWHCay3OuguuHEJUTJhnjfQwys/aQIuatYi7Qd0IjAKuFEK8I4So22/Qf7sfHiYxMdHXJngdpblpCpf9AMDAp+cjgux/Ros2HaOoupZB4sSwe+bgdM8b6WGUn7WBFjVrEbd3wJBS5gLnYF/etl4I0ZkTy9tahWNYv1AI8Vu9soeEEMeEENsdr6n16u4RQuwXQuwRQkyuV54thNjhqHtetMMUY71e7+lb+j1Kc+OUrN9K+Wb7RzZ9uv1juPFwBSv2lRIWBJ3L7ZPjLrhySPsZ6kGUn7WBFjVrEXeDugCQUhqklH8EPgZ+Ato6rvgmMKWR8n9JKYc6Xl8DCCGygCuBgY5rXhZC1E0nXgTcCPR1vBq7Z5swm82evqXfozQ3zs+X3eo8DomKoMJo4dl1RwAYf6SYmkoTnRIj6TfQP/Z2bwnlZ22gRc1a5FS2iXUipXwCuB5Y0pY3l1L+AJS22NDOxcD7UkqTlPIQsB8YIYRIB2KllBullBJ4G7ikLXY1hhbXeCrNDVmWNhpsNsCesAXgxfVHKTNYyI4KwWa2ApB9Vg+/2gq2OZSftYEWNWsRd5e0Pd1I2TfANx63yM7NQohZwGbgDillGdAZ2FSvTa6jrNZxfHJ5AwoLC5k9ezYhISFYrVamT5/OvHnzyM/PJyoqiuDgYCorK0lOTqa0tBQpJcnJyRQUFGAwGEhOTqaqqorU1FSKiooQQpCQkEBRURGxsbFYrVaqq6tJS0sjPz+f0NBQ4uLiKC4uJi4uDrPZjMFgcNbrdDpiYmIoKSkhPj4eg8GA0Wh01oeHhxMREUFZWRmJiYno9XrMZrOzPiIiAp1OR0VFBUlJSVRUVFBbW+usb0lTdHQ0QJOaDAYDSUlJAaWpJT8ZjUan3Sdr2vXAc87PUljvrsSfN5qPf9rHmkNVRARD4o5jzvqEdInJZPILTS35yWQyERcX16H81NbPntlsJjo6OqA0teQns9lMVFRUqzQpOg7C3rltpEKIZVLKKY7jtTTx/FxKOaZNBgjRA/hSSjnIcZ4KFDve71EgXUr5ZyHES8BGKeW7jnavY9/p7gjwuJRyoqP8HGC+lPLCk99r48aNMjMzs1V2FhYWkpKS0qprOypK8wmMBcWsHnKR83zS0R/QW+H6/+2mwmjhUmqpPGgfdJp3/7lERHac3NXKz9qgLZq3bt26ZcKECad72CRFO9BcT/3teseLm2zlYaSUBXXHQojXgC8dp7lA13pNuwB5jvIujZR7FJ2u4/yT9hRK8wnqB/TJx9cjhODldTlUGC0MTY6g8kf7x/b0c3p0qIAOys9aQYuatUiTQV1KuaTe8VveMQeEEOlSyuOO02lA3cz4pcASIcQz2Les7Qv8JKW0CiH0QohRwI/ALOAFT9tVUVFBp06dPH1bv0ZptlNz+MSwep/5NyCEYMPhclYdKCMsJIiuO45hctSPndzfi9Z6BuVnbaBFzVqk2WfqQohswCSl/M1xngw8CwwCNgJ3SilbnalNCPFfYByQJITIBR4ExgkhhmIffs8B/gIgpdwphPgQ2AVYgHlSSqvjVnOxz6SPwP6c3+PP+pOSkjx9S79HaQaLvpofRl7mPO/9f9eiN1l4fv1RAC5PCadgby0AF/1xKMKPd45rCuVnbaBFzVqkpYlyzwIPc6K3vBh7L/lV4CrgKeCm1r65lPKqRopfb6b9Y8BjjZRvxv5Fo92oqKggKiqqPd/C79C6Zikl3w04sTpy8EsPIoTgtR/zKK2xkJUcScG6Q876foM65uxirftZK2hRsxZpKagPANYCCCE6AecDg6SUe4UQS4ENtCGodyRqa2t9bYLX0bJmW62Fb7uemAOadtEEMmZMZnuenmV7SwgNEowornRO3ph9+zk+sNYzaNnPWkKLmrVISwtpQ4C6HQtGAflSyr0AUsqjgGYe0GhxjadWNVsNJtaN+aOzrNPpgxj66qOYLDaeXWcfdr8wQpB3oASAP1wxmPikjtsD0qqftYYWNWuRloL6TqDugeKVgDOLhWOb2Ip2ssvv0GIuYq1q3n3fM9Qcsm990PmKqYz68lUA3t2WT16lie6dwqjaYZ/LGRGlY8CQDJ/Z6wm06metoUXNWqSl4fe7gS+EEP8GrMDZ9equANa3l2H+hhafRWlRM4fyyF3yBQB97rqePnf8GYADJTV89GsBArgiJYLtjuazbh7tGzs9iBb9rDQrApVmg7qUcp0QohvQD9grpayfEeAr4P32NM6fCA4ObrlRgKE1zdJqZfeVdzjP6wK61SZ5dt1RbBIu7B7H9mX2POpjz+9PTFy4T2z1JFrzMyjNisClxc2ppZR6KeWWkwI6Uso9UkqPb/Lir1RWVvraBK+jNc0r+05yHo9e8Ybz+MvdxewpqiEpMpShwr7ve0JyFNmju3vdxvZAa34GpVkRuHSMjBN+QHJysq9N8Dpa0myrtWCtMQDQddYlxJ5m30SmuNrMG5vt311vOrMz+3+zP5cce37/DpOwpSW05Oc6lGZFoBIY/5W8QGmpu8nkAgctad523QLn8cCn5juPX954jJpaG2d2jyPDbKGksIrwiFB69AmcjTy05Oc6lGZFoKKCups0lfgmkNGK5qPvfk7Ryg0AdHn8r87yTUcqWJdTTnhIEDeN6swnb20BYODwDIJDAudPRyt+ro/SrAhU3Eq9WocQIgVwycUnpTzoUYv8FC0OXWlBc8m6zey880kA0mdMos9V9sQthlorL22wL2u7Jjuddx5f5bxm9IQ+3je0HdGCn09GaVYEKm51N4QQU4QQx4DjwP56r33taJtfUVBQ0HKjAEMLmvc+/goAYenJDHzqbqfmJdvyKagy0zsxgjOiT3z3FUGCsPBQn9jaXmjBzyejNCsCFXfHEF/Cnts8WkoZVO+lmTUS0dHRLTcKMAJds6mwhIotOwEY+ORdhERFEB0dTU6Zgf/tKEQAt57Vle0bjzivuf3RSU3creMS6H5uDKVZEai4G9TjgVeklIb2NEah8BY2cy2rBl/oPE8+7yzA/tzxhfW5WCVMzUwkWdrY+1s+wcGCG+ePRYiOl4VNoVBoB3eD+uvAde1piL9TVdXqDLMdlkDWvOuefwAggoMZ8+P/nMF6xb5SduRXERcewnWnp/PWC/YJdFnDOhPbKcJn9rYngeznplCaFYGKuxPlRgG3CiEWAC4bCEspxzR+SWCRmprqaxO8TqBq3vf0YnLfc2wFO/96Irvb927Xmyx8st8IwF9GdiZ/fwnSZp8xPHJcL98Y6wUC1c/NoTQrAhV3g/pix0uzFBUV0bVrV1+b4VUCUXNNTi4H/vkfADIfvpUef7nSWffG5uNUmKwMTovmnK4xvPW8PbXBqPG96ZQQ6RN7vUEg+rkllGZFoOJWUJdSvtXehvg7WnyWGmiaayv0/DDqcud59xuvcB7vKarmq93FBAu4+awufPSfzVSWG0lKi+bMc3v7wlyvEWh+dgelWRGouL2DhhDiOiHE90KIPY6fmnrGnpCQ4GsTvE6gaf6u/2Tn8Vmr3nH+k7Pa7JPjJHBhZjyfP7eO/Fx7VuHzLh5IcIBsB9sUgeZnd1CaFYGKu+vU7wMWYM/Kdqvj53xHuSYoKirytQleJ5A0VzvyowN0u24GMQNO9L6//r2YvcU1JEWF0qe8nFqzFYCktGg6d4/3uq3eJpD87C5KsyJQcfeZ+vXAOCnl4boCIcRy4AfgsfYwzN+IjY31tQleJ5A0/3rTQwB0vmIqWY+fSK9abqjljc3HAbh+aCpb37FvBZvWJY4/zhnldTt9QSD52V2UZkWg4m5QjwJO/ppXAgTmGp9GsFqtvjbB6wSK5twlX1KxbRcAPW+62qXu9Z/zqDJbye4cQ3x5DQBJqdFcPXeUZp5BBoqfTwWlWRGouPuwcBnwnhCivxAiQgiRCbwFLG8/0/yL6upqX5vgdQJBs7XGyG+3LwQgul9Povv3dNbtLKhi+d5SQoME1wxI5LsvdgMwfHR3zQR0CAw/nypKsyJQcTeo3wzogV+AKmA7UA3c0k52+R1paWm+NsHrBILm7Tfe7zwe8elLzmOrTfKiI2HLpYNT+OnL3c66PlnaWs8bCH4+VZRmRaDiVlCXUlZKKWcBkUA6ECmlnCWlLG9X6/yI/Pz8lhsFGB1dc/6Xq5wpVTP//ld0iZ2cdV/uLuZAiYHUaB0TUiIpzKsEYMT4DCKjdD6x11d0dD+3BqVZEag0+UxdCNFDSpnjOD55O63ouuFJraReDQ0NrMxc7tCRNdssFrZff2JxRvfZlzmPy2pqeXOLfXLc7KEpfPjvTc66XgMSvWekn9CR/dxalGZFoNLcRLkdQIzjeD8ggZMfNEpAE5na4uLifG2C1+nImnff+y/n8YS937o8I1/8cx7VZisjusby++c7neWXXnc68fGamfvppCP7ubUozYpApcnhdyllTL3jICll8ElpVzWVerW4uNjXJnidjqp51z3/5OjbnwIw7M0nCI09kXLyt/wqVuwrJTRYMNpkoqrCvtf75defQY++SR1Wc1tQmrWBFjVrEXc3n3m+ifJnPWuO/6LFb7kdUXP1waMceeNjAFLOH0PqlBP5huyT444CMC09il2b7HnS+5+WRrde9mH3jqi5rSjN2kCLmrWIu7Pfr22ifKaH7PB7zGazr03wOh1Ns6WqmrWjT+znPmyx675IS3cVcbDUSNdQQfkPB+xtzuzGhVcNdbbpaJo9gdKsDbSoWYs0u/mMEOLPde3qHdfRC9DMeI7BYPC1CV6no2k+8OyJvENDX/s7IvjE06GSmlre2nIcISUD9hQ4y8+a2NflHh1NsydQmrWBFjVrkZZ2lKvrietw7ZVLoAC4pj2M8ke0uMazI2k2HCsgZ9F/Ach68i7SLjzXpX7xT8eoqbVxdr0JwDNvHk14hOuM4I6k2VMozdpAi5q1SLPD71LK8VLK8cATdceO17lSyquklJuauz6Q0OIaz46keU32NKTVSvKks+l2zTSXul+P6/lufxm6IOjs2Aq236BUUjMa7oXdkTR7CqVZG2hRsxZxd/MZ57Zcwk5Q3av9TPMvdDptbUgCHUfznkdO7BTX48bLXeosNskLjp3jpmVEU1ZYRWSUjqmXDW70Xh1FsydRmrWBFjVrEbcSugghMoCXgDFAp5OqNbGsLSYmpuVGAUZH0LzlT3c6d41LOvdMEs8+3aX+s51FHC4zkh6jo1NBJWXAkJFdCQlt/GPbETR7GqVZG2hRsxZxt6f9CmAGJmDf+304sBSY0052+R0lJSW+NsHr+Lvm2vJKZ0AHGP7Wky71xdVm3tlq3zluZq84Du2xJxocOrJbk/f0d83tgdKsDbSoWYu4G9RHA3+WUm4HpJTyF2A2cEfzlzWPEOI/QohCIcRv9coShBArhBD7HD/j69XdI4TYL4TYI4SYXK88Wwixw1H3vGiHFFvx8fEtNwow/F3z3oWvOI8n7FlOUKjrwNMrm45hqLUxulss2z+1f8SS02OIiglr8p7+rrk9UJq1gRY1axF3g7oVsDiOy4UQydiztHVu4/u/CUw5qWwB8J2Usi/wneMcIUQWcCUw0HHNy0KIujHURcCNQF/H6+R7thktLgfxZ83lW3c5d407a/W7hMa5Di1uya1kzaFywoIFnTacSE9w+ewzmr2vP2tuL5RmbaBFzVrE3aD+IzDVcbwc+AD4BNjcljeXUv4AlJ5UfDH2XO04fl5Sr/x9KaVJSnkI+370I4QQ6UCslHKjlFICb9e7xmMYjUZP39Lv8VfNUko2Tb0eABEcTEyma74hs9XmTKt6WdcYLGYrANGxYURENj9ZyF81tydKszbQomYt4tZEOexr1Ou+APwVuBOIBtpjm9hUKeVxACnlcSFEiqO8M1B/CV2uo6zWcXxyuUfR4hpPf9V89M1PnMdnLn+9Qf1HvxZyrNJE91gd5esPAfaA/pe7x7V4b3/V3J4ozdpAi5q1iFtBvX7edCmlAXi03SxqmsaekzeWOa6uvAGFhYXMnj2bkJAQrFYr06dPZ968eeTn5xMVFUVwcDCVlZUkJydTWlqKlJLk5GQKCgowGAwkJydTVVVFamoqRUVFCCFISEigqKiI2NhYrFYr1dXVpKWlkZ+fT2hoKHFxcRQXFxMXF4fZbMZgMDjrdTodMTExlJSUEB8fj8FgwGg0OuvDw8OJiIigrKyMxMRE9Ho9ZrPZWR8REYFOp6OiooKkpCQqKiqora111rekKTranuikKU0Gg4GkpCS/0mTbdYg99/wTgLQrp1LZKZKKo0edmqpFOEu22dfjnmsxccxkf2p08czBHDlypEU/GY1Gp90dxU9t/eyZTCbi4uICSlNLfjKbzURHRweUppb8ZDabiYqKapUmRcdB2EesG6kQ4hF3biClfKBNBgjRA/hSSjnIcb4HGOfopacDq6WU/YUQ9zje73FHu+XAQ0AOsEpKmekov8px/V9Ofq+NGzfKzMzMVtlZUFBAampqq67tqPij5mVpo53H5x38nuDIcOe5lJL7lh9gc66eibZagnLsT3auuH4EXXsluHV/f9Tc3ijN2qAtmrdu3bplwoQJp7fcUuFrmnum3tXNl6dZyontZ68BPq9XfqUQIkwI0RP7hLifHEP1eiHEKMes91n1rvEYERHay7Ptb5rLfvzFeTxu62cuAR1gbU45m3P19DCYnAF9wJB0twM6+J9mb6A0awMtatYiTQ6/Symva+83F0L8FxgHJAkhcoEHgSeAD4UQs4EjwGUOe3YKIT4EdmGfiT9PSml13Gou9pn0EcA3jpdHKSsrIza24baigYw/abZZLPx48VwA4kcNJTwjxaW+2mxl0cZjpFQb6VdQ4Swf/4cBp/Q+/qTZWyjN2kCLmrWIuzvK9WqqTkp5sKm6lpBSXtVE1YQm2j8GPNZI+WZgUGvtcIfExMT2vL1f4k+aDz7/jvN4+NtPNah/e8tx9JVGxjkC+sDhnZkyYxCnumWBP2n2FkqzNtCiZi3i7pK2/cA+x8/99c73tZNdfoder/e1CV7HXzTnf7Wa/U+9BkDvO/5MaKzr5J29xTV8vquIXhU1zrKJF2WdckAH/9HsTZRmbaBFzVrE3dnvLsFfCJGGfah8bXsY5Y+YzWZfm+B1/EXz9tn3ApB4zun0uXO2S53VJnl27RGyCirIqLKvw51xbTahutalJPAXzd5EadYGWtSsRVqVZU1KmY99vfrjnjXHf9HiGk9/0Kz//cTTnQGP39Gg9/35riJs+4ucAR2gR5+kVr+fP2j2NkqzNtCiZi3SltSp/YFITxni72gxF7GvNVtrjGy79m4AgqMiie7T3aW+sMrM2z8do3dZtbPs1ocmIoJav/W/rzX7AqVZG2hRsxZxd6LcWlw3dInEvge7W2vZAwEtLgfxteaV/c5DWuwLHEZ+/rJLnZSSlzbm0q2wEsGJHePamsvH15p9gdKsDbSoWYu4u03s4pPOq4FfpJSamSin0zW/Z3gg4kvN+Uu/dwb0zEdvI3ZQP5f6dTkVbD5YxthKe5KKcVMz2xzQQflZKyjNikDF3Ylyb7XcKrCpqKigU6dOvjbDq/hK855HXuLQy+8B0Gf+DfS44QqXer3JwksbjpKhNxAEdO2ZQObgdI+8t/KzNlCaFYGKW8/UhRA6IcQjjhzn1Y6fjwohwlu+OjBISmr95KuOii80l27Y5gzoMYP60vu2WQ3avPZjHqU1tfQvrQJg+FndG7RpLcrP2kBpVgQq7k6UWwScC9wKnOH4ORZ4ubmLAomKioqWGwUYvtC8z7EeHeDMZa8jgl2Xpm3L07Nsbwn9602Oa8ts95NRftYGSrMiUHH3mfolQO962dp2CSF+xL4JzZ/bxTI/o7a21tcmeB1vay7dtJ2yTdsJjopk3JZPCApx/XgaLTaeW3eEyFoL3cvtQf3sSX1bvSa9MZSftYHSrAhU3O2p59Nw+VoEcNyz5vgvWlzj6U3NVoOJbdctACB16lhCOzXco/qtzXkUlBs5+2iJs2zk2CZ3MG4Vys/aQGlWBCruBvV3gGVCiBuEEOcLIW4EvgbeFkKcW/dqPzN9jxbXeHpLs5SSFT3HU1tWCUD/v93UoM3uwmo++a2ICTlFzrI5C9q+hO1klJ+1gdKsCFTcHX6vy01+70nlcxwvsK9j92y3yY+IiorytQlex1uaf5n7oPN4xGcvE5bimnjCbLHxjzWH6V9c6SybNms40bGen6ep/KwNlGZFoOLukrae7W2IvxMc7Lnnth0Fb2g25hWS/9lKALpffxkJo4Y2aPPOtnx0h0ro5liT3m9QKr0zUxq08wTKz9pAaVYEKm5vEyuECBFCjBFCXCWEOEcI4W4vPyCorKxsuVGA4Q3NuUu+cB4P+Pv/NajfU1TNp9uO09sxMW7EmJ5c9Mdh7WaP8rM2UJoVgYq728RmAl9gnxx3FOgKGIUQF0opd7ejfX5DcnKyr03wOu2tOfe/X7L/H68DkP3fZxrUmyw2nlp9mB71lq+dPalfg3aeRPlZGyjNikDF3Z76y8CrQFcp5ZlSyi7Av9HQOvXS0lJfm+B12lOz8XgRu+79JwCpfxhH0tgRDdq8uTmP/FID3SvtedL/OGcUQW1I1uIOys/aQGlWBCruBvWhwDNSyvpJXZ51lGsCV+naoL0022otbL1mPjaDiZTzxzDs9YWIINeP4q/H9Xyyo5BzjhQjJHTrnUhGt/bf4lL5WRsozYpAxd2gnod9B7n6nOMo1wRaHLpqL82rhlxE5a97AOi34C8N6mvMVv7xwxFSq42EOP4RnTO5fYfd61B+1gZKsyJQcTeo3wssFUK8L4R4UgjxPrCUhkvcApaCggJfm+B1PK1Z2mysn3gNtaX2jQl7zP0j0f0bLqxYtCmXqpIaBhfaJ/acd8lA0rvEedSWplB+1gZKsyJQcSuoSymXAsOA34AYx89sKeXn7WibXxEdHe1rE7yOpzUffO4t9L/Zs/VmXDqFzAdvbtBm3aFyNuwoYHSufde4xJRoBg7L8KgdzaH8rA2UZkWg0uzsdyFEJHA/MAjYCjwupTR5wzBFYGE1mtj3pD1ZS9dZ0xj41F0N2hRXm3lp1SGG5Zc7yy6ffQYhoWp9rUKhULhDSz31F4ELgd+BS4F/tLtFfkpVVZWvTfA6ntIsbTZW9BjvPB/wWMP16DYp+ddXe8neX0C41UZKRiy3PjiRqJgwj9jgLsrP2kBpVgQqLQX184FJUsr5juML2t8k/yQ1NdXXJngdT2nOX/qd87jvvXMICm04QLRk1SHit+cCEBYRwqXXno4uzPv7Gyk/awOlWRGotBTUo6SUxwGklEcB78xW8kOKiopabhRgeEKzuaySPY/atzPoMecqet86q0Gb3fl68lfudZ5fdt0ZREbr2vzerUH5WRsozYpApaWuUIgQYjwgmjhHSvl9exnnT3g6G1hHwBOa1wy/BKvBSOyQTPrdP7dBfbXZyusf7qAuKeRf7h5HTJznE7W4i/KzNlCaFYFKS0G9EPhPvfOSk84DOjNb5YmbNgAAGHJJREFUfRISEnxtgtdpq+bSjduwGowADPrnAoJCXD9uUkqe+2w3Kfn2pWsTpw30aUAH5WetoDQrApVmh9+llD2klD2beWkioIM2h67aorng6zX8NG0eADFZfYgd1HDzmM+2Hif4l1yCgN5D0hl6RtdWv5+nUH7WBkqzIlDRVKa1thAbG+trE7xOazWXb93Ftj/fA0Bkr66c+c1il3qLxcaK7/ZzYM1BgoHI5CgumnFaW831CMrP2kBpVgQqKqi7idVq9bUJXqc1mqXNxqap1zvPR33xCkFhJya9SSl568X1lBWeyLz2x1nZBIe4nQW4XVF+1gZKsyJQ8Y//pB2A6urqlhsFGKeqWVqtLM8423l+1up30SW6JmH5dXOuM6CbosK49o5z6JQY2XZjPYTyszZQmhWBiuqpu0laWlrLjQKMU9W8d+G/ncfZ7/6DmEzXKRc/rj7A2m/t28QWxUZw+7wzSYrxzdK1plB+1gZKsyJQUT11N8nPz/e1CV7nVDQbjh7n0EvvAZA+YxLJE0e71OsrjM6ADnDR9IGk+llAB+VnraA0KwIVFdTdJDQ01NcmeJ1T0bz12gUARHRNZ/CLD7rUSSn5ZMl253nGjMGc1S/JM0Z6GOVnbaA0KwIVNfzuJnFx2ttMz13N+t0H0O+098JPe/7+BptcrFmxj6Kj5dQGCWyjenLV8HSP2+oplJ+1gdKsCFRUT91NiouLfW2C13FXc84r7wPQ7drpJJw5zKWuxlDL5tUHASjvEs/tk3v79c5Wys/aQGlWBCp+G9SFEDlCiB1CiO1CiM2OsgQhxAohxD7Hz/h67e8RQuwXQuwRQkz2tD1a/Jbrjub8L77n2PtfAZA+7TyXOikl//7vL87z264aTISfp1FVftYGSrMiUPHboO5gvJRyqJTydMf5AuA7KWVf4DvHOUKILOBKYCAwBXhZCOHR6GE2mz15uw5BS5ord+7jlzn25+cZl04hfuQQl/pnHlyBbb+9d5B9YRYZPt4C1h2Un7WB0qwIVPw9qJ/MxcBbjuO3gEvqlb8vpTRJKQ8B+4ERnnxjg8Hgydt1CJrTXHM4jw0TrkFarSSNH8Wgf93rUv/Zyv1Iiw2AjGGdGX9mt3a11VMoP2sDpVkRqPjzRDkJfCuEkMArUspXgdR6qWCPCyFSHG07A5vqXZvrKHOhsLCQ2bNnExISgtVqZfr06cybN4/8/HyioqIIDg6msrKS5ORkSktLkVKSnJxMQUEB4eHhlJSUUFVVRWpqKkVFRQghSEhIoKioiNjYWKxWK9XV1aSlpZGfn09oaChxcXEUFxcTFxeH2WzGYDA463U6HTExMZSUlBAfH4/BYMBoNDrrw8PDiYiIoKysjMTERPR6PWaz2VkfERGBTqejoqKCpKQkKioqqK2tdda3pCk6OhqgSU0REREUFxe7aAoJCaHqw2859ORrzt9rwt9upNZmJf/wMXQ6HT9tLWP/+iMAhHaO4awRsZhMJr/Q1JKfIiMjKSws7FB+autnLyoqioKCgoDS1JKfYmJiOH78eEBpaslPsbGx5OXltUqTouMgpJS+tqFRhBAZUso8R+BeAdwCLJVSdqrXpkxKGS+EeAnYKKV811H+OvC1lPLj+vfcuHGjzMzMbJU9hw8fpnv37q2V0yFpTPPKvudh0dt3pgqOjmT4m0+SeHa2s/733Eo++/dGQmwSonT8391jCQ7x7+fo9VF+1gZK86mxdevWLRMmTDi95ZYKX+O3w+9SyjzHz0LgU+zD6QVCiHQAx89CR/NcoH6Kry5Anift0en8b6OU9uZkzQdfeMcZ0Lv9+VIm7v3WJaAXVpl5a8l2e0AHbr2rYwV0UH7WCkqzIlDxy6AuhIgSQsTUHQOTgN+ApcA1jmbXAJ87jpcCVwohwoQQPYG+wE+etCkmJsaTt+sQ1Nd88MV32fvYIgCiM3uRtfB2RNCJj8+RvEr+tfhnEstrkEGCmbedhU7XsQI6KD9rBaVZEaj46zP1VOBTx3rmEGCJlHKZEOJn4EMhxGzgCHAZgJRypxDiQ2AXYAHmSSk9mpKopKREc8+X6jQXfb+JvX9/GYDkiaPJfvcfLu2++ngHu7ccI9FxfsaYnqSmdsx/IFr2s5ZQmhWBil8GdSnlQWBII+UlwIQmrnkMeKy9bIqPj2+5UYARHx+PzWRm5/ynnGXD337Kpc2v2/LYveUYAMbwUKZM6Uf2GV28aqcn0aqftYbSrAhU/HL43R/R4nIQfV4B33YfhzHXnghi/I4vXYbcDx0q5duPfgWgIi6COXecw+kjuvr1jnEtoUU/K83aQIuatYgK6m5i/P/27jw8ruq84/j31TKa0WhG1mZJNvKCwcY2XjEYCiWLgUAWoGEpoQTakKYBQkoIDw8kJISmNElLoJACSYjbsAQMdlxjeKAG7KQhQZjdtfESbCxHQttosSzJkkbL6R9zPcgEW4M80syc836eZx7de+beq/Pzlfxq7j333t7eVHdhXLVv3MTmj18Rnz/5mQfIKyt+//3OPn79QGzYQmcwj6uvWkppMPMH4ri2n0Ezu8LFzC7Sop4gl55F3LVzDxvPuyo+f/yd32LC4rnx+ZaOXu6956X4/KWXLaQynP53i0uES/v5AM3sBhczu0iLeoJceBaxMYY9y1fx+9O+EG87Y9d6jrr0s/H52sZOfvmj3+Lv7gNg2V8vYNZUe87VubCfP0gzu8HFzC5Ky4Fy6cjvt+OT6KEMdO9n01e+Q2R9dbxt1oofkxMMxOdrmrtZcV91/Idm4enTWbQgfR+jOhq27+cPo5nd4GJmF2lRT1AgEBh5oQzVuW0XW274IR2vv02W38eM67/E0dd+kc7Ozvgyr21t4rkVm/B593M/+9KFHH+8fYfzbN7Ph6KZ3eBiZhdpUU9Qe3s74XA41d1IKmMMNfc9yo7b74ehIXylRZy48h5Cs2cAsczBYAErHttEw9YmfEB/0Mdlly+mqmrC4TeeoWzczyPRzG5wMbOLtKgnqKSkZOSFMogxhm233MWflq8CoPKCszju1mvJm/h+ztaGQVb+7Ln4/GBpkK9ftZRgIPNHuR+Kbfs5EZrZDS5mdpEW9QR1dnZaczemge4ett50B/Urn0Vyspn3k+9Qef6Z8evLh4YMj/70ZRrrOuLr5B9fyVe/MJ+sDL4GPRE27edEaWY3uJjZRVrUExSNRlPdhaRof3Uzr154LUN9UbIDfubdcwsVn/vkQcs8fF81kfp98fnjLpjPZ0+YNN5dTQlb9vNHoZnd4GJmF2lRT5AN13g2rt3AW1+5BYCcUJCT1txHeO6x8ff3d0dZ9cib8YL+bkUhX71oHvMr3fnr3ob9/FFpZje4mNlFep16gjL9Gs9oWwfbb/tJfP70l1ceVNC3vvke992+geY97QDsL8rn7z9V7lRBh8zfz6Ohmd3gYmYX6Sf1BGXq5SDGGFrWV7PpmtsY6OikYNZ0Tt3wEJIdeyzq4OAQTz2+iZ1bmuLryNJp3PSZmexta0lVt1MmU/fzkdDMbnAxs4u0qCfI58usEd8D3T1s+eYPaFzzQrwtO+DnhEfvjBd0gFUrNlH7dqygR0J+zrxwHmccGxslm2mZk0Ezu0EzK1tpUU9QR0cHEyZkxrXZPbUN/O/Si2BoKN4WmDqJJY/dRWByOcYY6mr3surRTQzuiz3kob0izDVfXMSUovf/ms+kzMmimd2gmZWttKgnqLS0NNVdSEjz83/gzSu/FS/oS5/+GYULZpOVm0O0b4BXX9zNa9V76N4bK+aDAqHFVXz9/Nn4sg8eYpEpmZNJM7tBMytbaVFPUEdHB8FgMNXdOKTIb17m3bsfpP3lTfG2T2x+Ov641M2v1bFu9Zb4e9EsoXtCPhddPI+FUz78r/d0zzwWNLMbNLOylRb1BPX396e6C4dU8/PH2f7du+PzZWeeyvF33kxeWTE9+6OsW72FnVubAejKzebd4gJOO6mKbyypJJCbfajNpnXmsaKZ3aCZla20qCcoHa/xHOzpY/PXv0/jUxuA2K1e5/zgBnLDBUQaO3nqzhdpa+mOL18XCjAwu4Jb/7KKGSX5I24/HTOPNc3sBs2sbKXXqScona7xNMbQ9D+/o/rsL8UL+qzvfo0F936PnFCQN6v38PC9Lx1U0N+ZNIFl583h7nNnJlTQIb0yjxfN7AbNrGyln9QTlC7noiIvvMSO2++na9suALICeSx44Hay5s7jpfU7qdnVSn1N7AYydaEA7xYHOWteOf+2uJKw/6Pt7nTJPJ40sxs0s7KVFvUEZWcf+tzzWBuK9lO/ah1brv+Xg9pLTj+Ryhuv5pnfR2jd8If3lwe2lYaYPn8S31xSSdUE/6i+byozp4pmdoNmVrbSop6gffv2UVRUNK7fs+udGt574lnqHnmS/nbvASsiTLnir5h167W0dUR5fPmr9HZHGcgSWv0+enOyKZhZyo2nTWX2xCP7yzwVmVNNM7tBMytbaVFPUFlZ2bh9r64du6l54HHqHlkbb8sJFzDpwrM55vq/YzA/yOo1W6l5uwkZMuzNy+WNigksmFLIVxdVMLciOfdrH8/M6UIzu0EzK1tpUU9QW1sb+fmJDTAbDTM4yP6a92j47+fZecfyeHtw7rGErrgUmXkczXt7eeGBt4h29gEgQGOBn6IlVdy1qIJZZck9ZzbWmdORZnaDZla20qKeIGPMmGx3/5566n+9jtpfrqavuTXeLsfMoOeCv2FP1E/n9l7Yvu2g9aLZWRQuncoNH5tGRShvTPo2VpnTmWZ2g2ZWttKinqBkHrrq37uPHd+/l7pfPYX4cjHRfowI0ZlziB49i965i2hqH4AIQC/RLKE1P49oVhYTzRBzT5/OOScdddgbxySDi4frNLMbNLOylRb1BDU1NTF16tRRr2+MoXndi9Q+9CQtG6pjbUD7lNl0LjmdvnAxfQPewu2xiTZ/Lo0FfrpKgnxiZilnzSxO+BrzZDjSzJlIM7tBMytbaVFPUEHB6AefNTy5np13LKf7nRp6iyay94RP0ltSwf5J0xnK8j5tD0B3TjZt+T66cnPoLshj4cxSzptRxOLJIXKzx/8+QUeSOVNpZjdoZmUrLepjqOW3G6n5+RM0VG8mMv9Uui4+l2jBnz88pdOXw+aJYQjmcfK0CZw2rZDFk8P4c/SGf0oppRKnRT1BXV1dlJSUHHaZwZ4+Gteup2njFhp3N9OyP4v9E6fRffEnMNmxf+oBEepDftoCPjp9OUwsC7J0SiEXVBUypzxITpaMR5yEJJLZNprZDZpZ2UqLeoLKy8s/tH2gf4Bty9dS+9ZuIn9qY3/ZZLor58Ls+Qct1+r3sasoSG5JPgsmhTh3cphFk0KUBHPHo/ujcqjMNtPMbtDMylZa1BNgjKG5uZnCUBm7dzSx57U/EqnvYH+PIZoXhKx8KJ4Lxe+v02cGaQ36ySov5KhpRSybVsR1FQWUh3ypC/IRRSIRqqqqUt2NcaWZ3aCZla20qCdg9f3PUre9k/6CwmGtIQjEprK7OxiK9tJVWYZ/agVVx5Ry/NFFHFuaP+aXnY0lkfQ5FTBeNLMbNLOylRb1BETeqqW/NHYpiL+lnpx9EQaLQgSmV1IxfxozFp7K1NKgdQPbiouLR17IMprZDZpZ2cqqoi4iZwN3A9nAL4wxP0zGdk+46hy2P7ORmfOnMGfJEgoqSpz4qzcSiTh3XatmdoNmVraypqiLSDZwL3AmUAe8KiJrjTFbj3TbJy6awjHTQs494SgcDqe6C+NOM7tBMytb2XS8+CRgpzHmXWNMFFgBnJesjQ8ODiZrUxlDM7tBM7vBxcwusqmoTwZqh83XeW1J0d3dnaxNZQzN7AbN7AYXM7vImsPvxJ5E+kEHPZaoubmZK6+8kpycHAYHB/n85z/PNddcQ2NjI8FgkOzsbPbt20dZWRltbW0YYygrK6OpqQm/309raytdXV2Ul5cTiUQQEYqLi4lEIoTDYQYHB+nu7qaiooLGxkZyc3MpLCykpaWFwsJCotEoPT098fd9Ph+hUIjW1laKioro6emht7c3/r7f7ycQCNDe3k5JSQmdnZ1Eo9H4+4FAAJ/PR0dHB6WlpXR0dNDf3x9/f6RMB24beahMgUCAlpYWqzKNtJ/y8/Npbm62KtNI+ykYDNLU1GRVppH2UygUoqGhwapMI+2ncDhMfX39qDKpzCG2PI5PRE4BvmeM+ZQ3fzOAMeYHB5aprq42xx133Ki2v2fPHucGmWhmN2hmNxxJ5jfeeOP1ZcuWLUlyl9QYsOnw+6vAsSIyXUR8wCXA2mRtfM2aNcnaVMbQzG7QzG5wMbOLrCnqxpgB4GvAOmAb8IQx5u1kbX/16tXJ2lTG0Mxu0MxucDGzi2w6p44x5hngmbHY9sDAwMgLWUYzu0Ezu8HFzC6y5px6ItavXx8B9oxm3ba2ttLi4uKWJHcprWlmN2hmNxxh5qnLli0rS2qH1JhwqqgrpZRSNrPmnLpSSinlOi3qSimllCW0qCdARM4WkR0islNEbkp1f0ZLRKpE5Dcisk1E3haRf/Tai0XkeRF5x/taNGydm73cO0TkU8PaTxCRzd5790iaP+FGRLJF5E0RedqbtzqziEwQkVUist3b36c4kPkb3s/1FhF5TET8tmUWkf8UkWYR2TKsLWkZRSRPRB732jeKyLTxzKeSwBijr8O8iD3xbRdwNOADNgFzUt2vUWapBBZ70yHgj8Ac4F+Bm7z2m4AfedNzvLx5wHTv3yHbe+8V4BRid/J7Fjgn1flGyH498CjwtDdvdWbgQeDL3rQPmGBzZmK3hN4NBLz5J4C/tS0zcDqwGNgyrC1pGYGrgZ9605cAj6c6s74+2ks/qY9sTB8UM56MMQ3GmDe86U5i1/NPJpbnQW+xB4HzvenzgBXGmD5jzG5gJ3CSiFQCYWNMtYn99j80bJ20IyJHAZ8BfjGs2drMIhIm9p//cgBjTNQYsxeLM3tygICI5AD5QD2WZTbG/A5o+0BzMjMO39YqYFk6HalQI9OiPrIxfVBMqniH1RYBG4FyY0wDxAo/MNFb7FDZJ3vTH2xPV/8O3AgMDWuzOfPRQAT4L++Uwy9EJIjFmY0x7wF3AH8CGoAOY8xzWJx5mGRmjK9jYjf06gBKxqznKum0qI9sxAfFZBoRKQB+DVxnjNl3uEU/pM0cpj3tiMhngWZjzOuJrvIhbRmVmdgn1sXA/caYRUA3scOyh5Lxmb3zyOcRO8w8CQiKyGWHW+VD2jIqcwJGk9Gm/E7Soj6yOqBq2PxRxA7rZSQRySVW0H9ljDlw38gm75Ac3tdmr/1Q2eu86Q+2p6NTgXNFpIbYqZNPisgj2J25Dqgzxmz05lcRK/I2Zz4D2G2MiRhj+oHVwF9gd+YDkpkxvo53GqOQPz/cr9KYFvWRjemDYsaTd25sObDNGHPnsLfWAld401cATw5rv8QbETsdOBZ4xTvE1ykiJ3vbvHzYOmnFGHOzMeYoY8w0YvtugzHmMuzO3AjUisgsr2kZsBWLMxM77H6yiOR7fV1GbMyIzZkPSGbG4du6kNjvi35SzySpHqmXCS/g08RGiu8Cvp3q/hxBjtOIHUr7P+At7/VpYufM1gPveF+Lh63zbS/3DoaNAgaWAFu89/4D7+6E6fwCPs77o9+tzgwsBF7z9vUaoMiBzLcB273+Pkxs1LdVmYHHiI0Z6Cf2qfrKZGYE/MBKYoPqXgGOTnVmfX20l94mVimllLKEHn5XSimlLKFFXSmllLKEFnWllFLKElrUlVJKKUtoUVdKKaUsoUVdKaWUsoQWdaWUUsoSWtSVUkopS2hRV0oppSyhRV0ppZSyhBZ1pZRSyhJa1JUaRyLytoh8PNX9UErZKSfVHVDKJiLSNWw2H+gDBr35fzDGzB3/XimlXKFPaVNqjIhIDfBlY8wLqe6LUsoNevhdqXEkIjUicoY3/T0RWSkij4hIp4hsFpGZInKziDSLSK2InDVs3UIRWS4iDSLynoj8s4hkpy6NUirdaFFXKrU+BzwMFAFvAuuI/V5OBv4J+NmwZR8EBoBjgEXAWcCXx7OzSqn0pkVdqdR60RizzhgzAKwEyoAfGmP6gRXANBGZICLlwDnAdcaYbmNMM3AXcEnKeq6USjs6UE6p1GoaNt0DtBhjBofNAxQAk4BcoEFEDiyfBdSORyeVUplBi7pSmaGW2Ej6Uu9TvVJK/Rk9/K5UBjDGNADPAT8WkbCIZInIDBH5WKr7ppRKH1rUlcoclwM+YCvQDqwCKlPaI6VUWtHr1JVSSilL6Cd1pZRSyhJa1JVSSilLaFFXSimlLKFFXSmllLKEFnWllFLKElrUlVJKKUtoUVdKKaUsoUVdKaWUsoQWdaWUUsoS/w/tvVMiZ5fzWAAAAABJRU5ErkJggg==\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(de_times, de_pop, label='ODE')\n",
    "ax.plot(imp_times, imp_pop, label='Implicit IEBM')\n",
    "ax.plot(iebm_times, iebm_pop, label='Explicit IEBM')\n",
    "ax.set_ylabel('Population Size', labelpad=10)\n",
    "ax.set_xlabel('Time', labelpad=10)\n",
    "plt.legend(bbox_to_anchor=(1,1), loc='upper left', borderaxespad=0)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}