{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Individual based modelling\n",
    "\n",
    "Support for modelling alive or active particles, whose properties/state change in time, is being developed. At the current time, only support for modelling particle *mortality* is included. Examples describing how to include these processes in lagrangian experiments are given below. To ensure the workflow is clean and easy to follow, the examples are run in isolation using a sample particle set. Transport processes are not included. However, the same processes can be easily included when doing transport modelling.\n",
    "Configuration options are specified in the run config file under the secion `BIO_MODEL`.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "**Note:**\n",
    "    \n",
    "Where appropriate, the term *individual* will be used to describe particles that are alive. \n",
    "\n",
    "</div>\n",
    "\n",
    "## Mortality\n",
    "\n",
    "PyLag includes an abstract base class called `MortalityCalculator`. New mortality calculators can be added to the model by subclassing `MortalityCalculator`. At the current time, the model includes the following mortality calculators:\n",
    "\n",
    "1. **FixedTimeMortalityCalculator** - individuals are killed after a fixed, specified period of time.\n",
    "\n",
    "2. **ProbabilisticMortalityCalculator** - individuals are killed probabilistically, governed by a specific rate.\n",
    "\n",
    "Below we demonstrate how both of these can be used to control individual mortality. To begin, we do some initial setup work that is common to each example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "from configparser import ConfigParser\n",
    "\n",
    "import pylag.random as random\n",
    "from pylag.data_reader import DataReader\n",
    "from pylag.particle_cpp_wrapper import ParticleSmartPtr\n",
    "from pylag.mortality import get_mortality_calculator\n",
    "from pylag.processing.plot import create_figure\n",
    "\n",
    "# Ensure inline plotting\n",
    "%matplotlib inline\n",
    "\n",
    "# Parameters\n",
    "seconds_per_day = 86400.\n",
    "\n",
    "# Seed the random number generator\n",
    "random.seed(10)\n",
    "\n",
    "# Create the config\n",
    "config = ConfigParser()\n",
    "config.add_section('NUMERICS')\n",
    "config.add_section('BIO_MODEL')\n",
    "config.add_section('FIXED_TIME_MORTALITY_CALCULATOR')\n",
    "config.add_section('PROBABILISTIC_MORTALITY_CALCULATOR')\n",
    "\n",
    "# We need a data reader to pass to the mortality calculator. It\n",
    "# can be used to draw out environmental variables (e.g. temperature)\n",
    "# that affect mortality. In both cases below, it isn't used, so we\n",
    "# use the base class.\n",
    "data_reader = DataReader()\n",
    "\n",
    "# Set time stepping params\n",
    "simulation_duration_in_days = 20.0\n",
    "time_step = 100\n",
    "time_end = simulation_duration_in_days * seconds_per_day\n",
    "times = np.arange(0.0, time_end, time_step)\n",
    "\n",
    "# Helper function in which the model is run and mortality computed\n",
    "def run(config, n_particles=1000):\n",
    "    \"\"\" Run the model to compute mortality through time \"\"\"\n",
    "    \n",
    "    # Create the mortality calculator\n",
    "    mortality_calculator = get_mortality_calculator(config)\n",
    "\n",
    "    # Create the living particle seed\n",
    "    particle_set = []\n",
    "    for i in range(n_particles):\n",
    "        # Instantiate a new particle\n",
    "        particle = ParticleSmartPtr(age=0.0, is_alive=True)\n",
    "        \n",
    "        # Initialise particle mortality parameters\n",
    "        mortality_calculator.set_initial_particle_properties_wrapper(particle)\n",
    "        \n",
    "        # Append it to the particle set\n",
    "        particle_set.append(particle)\n",
    "\n",
    "    # Store the number of living particles in a list\n",
    "    n_alive_arr = []\n",
    "\n",
    "    # Run the model\n",
    "    n_alive = n_particles\n",
    "    for t in times:\n",
    "        n_alive_arr.append(n_alive)\n",
    "\n",
    "        n_deaths = 0\n",
    "        for particle in particle_set:\n",
    "            if particle.is_alive:\n",
    "                mortality_calculator.apply_wrapper(data_reader, t, particle)\n",
    "                if particle.is_alive == False:\n",
    "                    n_deaths += 1\n",
    "            particle.set_age(t)\n",
    "\n",
    "        n_alive -= n_deaths\n",
    "\n",
    "    return n_alive_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### FixedTimeMortalityCalculator\n",
    "\n",
    "The use of a `FIXED_TIME_MORTALITY_CALCULATOR` is speficied by setting the option `mortality_calculator` to `fixed_time` in the section `BIO_MODEL` in the run configuration file. The mortality calculator kills particles after they reach a given age. The age at which particles die can be set in one of three different ways. The choice is made through the option `initialisation_method` in the section `FIXED_TIME_MORTALITY_CALCUATOR` in the run configuration file. The choices are:\n",
    "\n",
    "1. **common_value**: It is set to a fixed, common value which is the same for all particles.\n",
    "2. **uniform_random**: It is set to a value drawn from a uniform random distribution with specified limits.\n",
    "3. **gaussian_random**: It is set to a value drawn from a normal distribution with specified mean and standard deviation.\n",
    "\n",
    "All values are given in days. The use of each option is demonstrated below. In the demonstration, we create a population of $1000$ individuals. In the **common_value** scenario, all particles are killed once they reach an age, $t_{\\mathrm{age}}$, of 10 days. In the **uniform_random** scenario, the age of death is drawn from a uniform random distribution with $t_{\\mathrm{age}} \\in [8, 12]$ days. In the **gaussian_random** scenario, $t_{\\mathrm{age}}$ is drawn from a normal distribution with a mean and standard deviation of 10 and 1 days respectively.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f4512016a30>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAHZCAYAAABw5MliAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABU5klEQVR4nO3de1yUZf7/8ddw9AQqCogiAuIBEUVFUTOPkVqestbcX7WapW22tW1tbbttW1vtaoctrfy22Ukz07QsLUu3PGSapng+pJKCCiKiguIB5TC/P24ZHQHxwMw9M7yfj4ePS65r7pnPOMKb+3RdFqvVakVERETcmpfZBYiIiMj1U6CLiIh4AAW6iIiIB1Cgi4iIeAAFuoiIiAdQoIuIiHgAH7MLcJSGDRsSGRlpdhkiIiJVJj09nSNHjpQ75rGBHhkZSUpKitlliIiIVJnExMQKx3TIXURExAMo0EVERDyAAl1ERMQDeOw5dBERT1NYWEhGRgYFBQVmlyIOVqNGDcLDw/H19b3ibRToIiJuIiMjg4CAACIjI7FYLGaXIw5itVo5evQoGRkZREVFXfF2OuQuIuImCgoKaNCggcLcw1ksFho0aHDVR2IU6CIibkRhXj1cy+esQBcREfEACnQREZFyREZGVjgrmytSoIuIiHgAXeUuIuKGHn0UNm2q2udMSIBJky7/mI8++ohXX30Vi8VCu3btePHFFxkzZgw5OTkEBwfz4YcfEhERwejRo6lZsyY7d+5k3759fPjhh0yfPp3Vq1eTlJTEtGnTAKhTpw4PPfQQ33//PfXr1+ff//43Tz75JPv372fSpEkMGTKEgoICHnzwQVJSUvDx8eG1116jT58+TJs2jQULFnD69Gn27NnDbbfdxssvv1xu3W+//TZpaWm28WnTprF+/XrefPNNhg0bxoEDBygoKOCPf/wj48aNs9s2PT2dQYMGsW3bNgBeffVVTp48yXPPPceePXt46KGHyMnJoVatWrz77ru0bt36ej6Ga6Y9dBERuSLbt2/nX//6F0uXLmXz5s1MnjyZP/zhD/zud79jy5Yt3HXXXTzyyCO2x+fm5rJ06VJef/11Bg8ezJ/+9Ce2b9/O1q1b2XT+t5FTp07Ru3dv1q9fT0BAAH//+9/57rvv+OKLL/jHP/4BwJQpUwDYunUrs2bNYtSoUbYrwDdt2sSnn37K1q1b+fTTTzlw4EC5td9xxx3MmzfP9vWnn37KnXfeCcAHH3zA+vXrSUlJ4Y033uDo0aNX/G8ybtw43nzzTdavX8+rr77K+PHjr/wftIppD11ExA1VtiftCEuXLuWOO+6gYcOGAAQFBbF69WpbUN5zzz08+eSTtscPHjwYi8VCfHw8oaGhxMfHAxAXF0d6ejoJCQn4+fkxYMAAAOLj4/H398fX15f4+HjS09MBWLlyJQ8//DAArVu3plmzZuzevRuAfv36UbduXQDatGnDvn37aNq0aZnag4ODiY6OZs2aNbRo0YJdu3Zxww03APDGG2/wxRdfAHDgwAFSU1Np0KBBpf8eJ0+e5KeffuI3v/mNre/s2bNX+K9Z9RToIiJyRaxWa6W3U1087u/vD4CXl5ft76VfFxUVAeDr62vb5uLHXfwYq9Va4etd/Lze3t62bcpz5513MmfOHFq3bs1tt92GxWJh+fLlfP/996xevZpatWrRu3fvMvd/+/j4UFJSYvu6dLykpIR69erZjjaYzWGH3MeMGUNISAht27a19R07dozk5GRatGhBcnIyubm5trEJEyYQExNDq1atWLx4sa1//fr1xMfHExMTwyOPPHLZD1ZERBynX79+zJkzx3ZI+tixY3Tv3p3Zs2cDMHPmTHr06FHlr9uzZ09mzpwJwO7du9m/fz+tWrW66ucZPnw4X375JbNmzbIdbj9+/Dj169enVq1a7Ny5kzVr1pTZLjQ0lMOHD3P06FHOnj3L119/DUBgYCBRUVHMnTsXMH7x2Lx587W+zevmsEAfPXo0ixYtsuubOHEi/fr1IzU1lX79+jFx4kQAduzYwezZs9m+fTuLFi1i/PjxFBcXA/Dggw8ydepUUlNTSU1NLfOcIiLiHHFxcTz99NP06tWL9u3b89hjj/HGG2/w4Ycf0q5dO2bMmMHkyZOr/HVLMyE+Pp4777yTadOm2e2ZX6n69evbDst36dIFgAEDBlBUVES7du145pln6Nq1a5ntfH19+cc//kFSUhKDBg2yu+ht5syZvP/++7Rv3564uDjmz59/7W/0OlmsDtzlvfTKwFatWrF8+XLCwsLIysqid+/e7Nq1iwkTJgDw17/+FYD+/fvz3HPPERkZSZ8+fdi5cycAs2bNYvny5bzzzjuVvnZiYiIpKSnX/R4OZx7mT38bet3PI54vNq8xPQ63rPLnPetXh5VdHqfIp0aVP3dVCAqCxx8HL11i63C//PILsbGxZpchTlLe5325bHPqOfTs7GzCwsIACAsL4/DhwwBkZmba/VYUHh5OZmYmvr6+hIeHl+mvyNSpU5k6dSoAOTk5VVJz3tFcPosoewhG5GLnzn8ndT1gwe+SU3hBZ2DITiPtLFbou9eL0FNXNq2jFyX4UsQbP3Xme6+bq7LkKlFcbPwZNAiUMyLmcomL4so7SGCxWCrsr8i4ceNs9w8mJiZWSW0t27XibDudt5fL+znjZ5747gkszez/f67cv5ISawlftim+qLeYns16Mrr9aH7X/nd4e3lX/MSpqdCyJQtf2ABPuV6gf/YZ/OY3cNH1QiKmS0pKKnO1+YwZM2xX2XsqpwZ6aGgoWVlZtkPuISEhgLHnffG9gxkZGTRu3Jjw8HAyMjLK9Iu4mqTwJFbcu6JMf3FJMel56VgxfilclraMvy75Kyv2rWDFvhW8tOolHkl6hK7hXekY1rHsE0dGGm1qqgOrF/EsP//8s9klmMKpZ72GDBnC9OnTAZg+fTpDhw619c+ePZuzZ8+SlpZGamoqXbp0ISwsjICAANasWYPVauWjjz6ybSPiDry9vGke1JyYoBhigmIY22ksR548wv5H9+Pv7c+uo7t46JuH6DS1E6/+9Co/pP/AibMnLjyBry/4+MCPP5r3JkTELThsD/23v/0ty5cv58iRI4SHh/PPf/6Tp556ihEjRvD+++8TERFhu9Q/Li6OESNG0KZNG3x8fJgyZQre3sZhyLfffpvRo0dz5swZBg4cyMCBAx1VsojTNK3blJN/O8mxM8eYuWUmj/3vMZ747gnb+J5H9hBdP9r4olcvWLLEpEpFxF049Cp3M1XVVe4izpCVn8XOIzt5d8O7zNo2C4ChrYyjUW9tCSd8whRIT4dmzUyssqzSc+jbtkFcnNnVeD5d5V69uPRV7iJSvrCAMMICwugT1Ye2IW2Zs30OW7K3kJaXRpeY0fwNYOVKlwt0EXEdunNUxMX87ca/sen3m9j7x734evmyyuegMfDNN+YWJnINUlJSbAu2nD17lptuuomEhAQ+/fRTkyu7NnXq1DG7hAppD13EhUXUjWBd9kZo0AA++QQ+/hgqmUtbxJUkJibabiPeuHEjhYWFVzX3eXFxse2aqutRVc/jyrSHLuLCOoR1IOd0Dsf7GqtCoamPpdSjj0Lv3lX759FHK33Z9PR0uzU6Xn31VZ577jl69+7NX/7yF7p06ULLli358fydGcuXL2fQoEEcPnyYu+++m02bNpGQkMCePXtYsmQJHTp0ID4+njFjxtjuHY+MjOT555+nR48ezJ07l8jISP72t7/RrVs3EhMT2bBhA/3796d58+b897//rbDW5cuX06dPH/7f//t/tnvQhw0bRqdOnYiLi7NNRAbGnvfTTz9N+/bt6dq1K9nZ2QCkpaXRrVs3OnfuzDPPPGN7vNVq5YknnqBt27bEx8fbjjgsX76cXr16MWLECFq2bMlTTz3FzJkz6dKlC/Hx8ezZs6fSf+NrpUAXcWHDWg0DoNcNxlKRvPuuecWIVKKoqIi1a9cyadIk/vnPf9qNhYSE8N5773HjjTeyadMmmjRpwujRo21rmRcVFfH222/bHl+jRg1WrlzJyJEjAWjatCmrV6/mxhtvZPTo0Xz22WesWbPGtmZ6RdauXcu//vUvduzYAVS89vmpU6fo2rUrmzdvpmfPnrx7/nvtj3/8Iw8++CDr1q2jUaNGtuedN28emzZtYvPmzXz//fc88cQTZGVlAdjWit+6dSszZsxg9+7drF27lvvvv58333zzOv+VK6ZD7iIu7K52d/HEd0+wOW8nD/42gFE/fkFXq1WH3cWcBdErMXz4cAA6depkW8u8Irt27SIqKoqWLY31D0aNGsWUKVN49PxRgtLV0EoNGTIEMNZMP3nyJAEBAQQEBFCjRg3y8vKoV69eua/TpUsXoqKibF9XtPa5n58fgwYNstX/3XffAbBq1So+//xzwFjv/S9/+QtgrNH+29/+Fm9vb0JDQ+nVqxfr1q0jMDCQzp0726Y5b968OTfffLOt9mXLll3+H/E6aA9dxMXN/Y0xX8N/W+XT7X5YuOJ9kyuS6qyitcHhwtrkla1LDpdf4xygdu3adl9fydrqlT3PxWufb968mQ4dOtjqv3hd9kvrL2/K8Stdo72iNd4dQYEu4uJuiLiBwmcK+V+kcWhx0PKx5J7JNbkqe545m4WUp6K1wa9W69atSU9P59dffwWMudZ79epVlaWWcSVrn1/qhhtusFvvvVTPnj359NNPKS4uJicnhxUrVtiWZDWLAl3EDfh4+ZDcczT3bTC+Dno5iDOFZ8wtSqqly60NfjVq1KjBhx9+yG9+8xvi4+Px8vLi97//fRVXa+9K1j6/1OTJk5kyZQqdO3fm+PHjtv7bbruNdu3a0b59e/r27cvLL79sd47dDJopTsRdnD4NtWvT/++R/M8nnd6RvVk2ynHn467E55/DHXfA1q1w0YXP4iCaKa56udqZ4rSHLuIuatWCBg346n8N8PP245ecX8yuSERciAJdxJ00a4bf2vX8NfFPZJ/K5n97/md2RSKm27p1KwkJCXZ/kpKSzC7L6XTbmog7uf9+GD+eEXtr8k+g/8f92fjARhIaJZhdmYhp4uPjr2r2OU+lPXQRd3L//QC0+WoNE/pNAGDUl6PMrEhEXIQCXcSd+PpC/fqwaBFP9XiKkNohbMneQuaJTLMrExGTKdBF3E3pDFpZWbw3+D0ABs4caGJBIuIKFOgi7qZ08o1t27ilxS3Ur1GfrYe3UmItufx2IuLRFOgi7qZFC6PNzMTby5s/d/8zAD8d+MnEokQc4x//+Afff/+92WVU6NLV58ykQBdxN23aGO3PPwMwPNZYEGPG5hlmVSTiMM8//zw33XSTQ57bkfOqm0G3rYm4m5o1jTYvD4DWDVvj4+XD1A1T+e+g/5a7kIR4nkcXPcqmQ5uq9DkTGiUwacCkyz7mhRdeYObMmTRt2pSGDRvSqVMn6taty9SpUzl37hwxMTHMmDGDWrVqMXr0aAYNGsQdd9wBGGuOnzx5kqysLO68805OnDhhWza1e/fu3HfffaSkpGCxWBgzZgx/+tOf7J7j+eef56uvvuLMmTN0796dd955B4vFQu/evUlKSmLZsmXk5eXx/vvvc+ONN5Zb/7Rp01i4cCEFBQWcOnWKBQsWMHToUHJzcyksLOTFF19k6NChpKenM3DgQHr06MFPP/1EkyZNmD9/PjVr1mT9+vWMGTOGWrVq0aNHD9tzFxQU8OCDD5KSkoKPjw+vvfYaffr0Ydq0aXz55ZcUFxezbds2Hn/8cc6dO8eMGTPw9/fnm2++ISgo6Lo/P+2hi7ij7t3h/IIRcGEv/Zcj5swe55kTSMulUlJS+Pzzz9m4cSPz5s2zTUE6fPhw1q1bx+bNm4mNjeX99y+/IuAnn3xC//79beuJJyQksGnTJjIzM9m2bRtbt27l3nvvLbPdH/7wB9atW8e2bds4c+aM3cIwl1uL/VKrV69m+vTpLF26lBo1avDFF1+wYcMGli1bxuOPP25bSS01NZWHHnqI7du3U69ePdsyqvfeey9vvPEGq1evtnveKVOmAMZEN7NmzWLUqFG21dy2bdvGJ598wtq1a3n66aepVasWGzdupFu3bnz00UeXrfdKaQ9dxB21bQs//QTHjkFQEL/v9HvmbJ/DrK2zeKHvC2ZXJ05Q2Z60I6xcuZKhQ4dS8/xRosGDBwNGWP39738nLy+PkydP0r9//8s+T+fOnRkzZgyFhYUMGzaMhIQEoqOj2bt3Lw8//DC33nqrbQ3xiy1btoyXX36Z06dPc+zYMeLi4mw1XM1a7MnJybY9YqvVyt/+9jdWrFiBl5cXmZmZZGdnAxAVFUVCQoLd8x4/fpy8vDzbynD33HMP3377re3f5+GHHwaM1eSaNWvG7t27AejTp49tDfe6deva6o6Pj2fLli2XrfdKaQ9dxB2VTmu5ahUAiY0TAXjlp1ecWoaO7lcvFa3lNXr0aN566y22bt3Ks88+a9srvXjtdKvVyrlz5wBj6dEVK1bQpEkT7rnnHj766CPq16/P5s2b6d27N1OmTOH+85MolSooKGD8+PF89tlnbN26lbFjx17zWuwXr5E+c+ZMcnJyWL9+PZs2bSI0NNT2vBeva176vFartcLTWmavka5AF3FHpeft9u0DIMA/gLvi7+Js8Vl2HdllYmHiyXr06MFXX31FQUEBJ0+eZOHChQDk5+cTFhZGYWGh3ZrhkZGRrF+/HoD58+dTWFgIwL59+wgJCWHs2LHcd999bNiwgSNHjlBSUsLtt9/OCy+8wIYNG+xeuzRkGzZsyMmTJ/nss8+q5D0dP36ckJAQfH19WbZsGfvOf09VpF69etStW5eVK1cCZddIL/169+7d7N+/n1atWlVJnVdCh9xF3FHz5uDlBYcO2brGdRrHzK0z+Xr317Rq6LwfIlJ9dO7cmSFDhtC+fXuaNWtGYmIidevW5YUXXiApKYlmzZoRHx9Pfn4+AGPHjmXo0KF06dKFfv362faMly9fziuvvIKvry916tTho48+IjMzk3vvvde2Rz9hwgS7165Xrx5jx44lPj6eyMhIOnfuXCXv6a677mLw4MEkJiaSkJBwReu7f/jhh7aL4i4+vTB+/Hh+//vfEx8fj4+PD9OmTbPbM3c0rYcu4q5CQ6F2bdi7F4CikiJ8X/DlmZ7P8Hyf551Swrx5cPvtsGULxMc75SWrNVdYD/3kyZPUqVOH06dP07NnT6ZOnUrHjh1NrclTaT10keoiOhrS0myXmPt4+RBaO5RZ22aZXJh4snHjxpGQkEDHjh25/fbbFeYuRIfcRdzVDTfAmjXGle4NGgAQ6B9I6rFUCosL8fX2NblA8USffPKJ2SVckcWLF/OXv/zFri8qKoovvvjCpIocT3voIu6qWzejffddW9e9Cca9u9tztptRkTiBh54lrXKl97lf/MedwvxaPmcFuoi7On8fKxdNrjEgZgAAk9ZMMqEgcbQaNWpw9OhRhbqHs1qtHD16lBo1alzVdjrkLuKu/PwgLs64Iu28hEYJAEzfPJ0Ph36oaWA9THh4OBkZGeTk5JhdijhYjRo1CA8Pv6ptFOgi7iwxEbZvh3PnwM8Pi8XC3e3u5uMtH/Ptr99yS4tbzK5QqpCvry9RUVFmlyEuSofcRdxZz55Gm5lp63rpppcA+Hr31+VtISIeSoEu4s6aNjXaixaJCKsTBkBaXprTytApXRHzKdBF3Fnp2uipqbYui8VCj4geLPp1ESfPnTSpMBFxNgW6iDtr3NhoL1mtaXBL4wr4eb/Mc3ZFImISBbqIO7NYIDbWbg8d4KHODwHwc8bPDn95EXENCnQRdxcdDVu3wunTtq7afrWp6VOTxXsWm1iYiDiTAl3E3ZXOGHfJso+tG7ZmT+4eCooKytlIRDyNAl3E3SUlGe2RI3bdI9uOBCD1aOqlW4iIB1Kgi7i7hg2N9pLZw5KaGEG/LH2ZsysSERMo0EXcXXCw0X7+uV33jc1uBOC55c85uSARMYMCXcTdNWpktJs22XV7WbyICYohtyCXvII8p5clIs6lQBdxd97ecM89sGMHnDplN/RMz2cAWJ6+3ITCRMSZFOginqBrV6P92n7+9r5RfQHYkLXB2RWJiJMp0EU8wd13G+2iRXbdTQKaAGj9bJFqQIEu4gkCA402L8+u22KxEOgfSP65fIe+vH5fEDGfAl3EUyQnw48/lumu41eHL3d+6fx6RMSpFOginqJmTTh6FI4ds+sOqxPGvuP72JK9pYINRcQTKNBFPMWwYUa7bZtd92v9XwPgr0v+WuUvqcVZRFyHAl3EU3TsaLSXzBjXs1lPgmsF803qN5RYS0woTEScQYEu4ilKZ4w7fLjM0L0J9wLwxS9fOLMiEXEiBbqIpyid033p0jJDT97wJADf7f3OmRWJiBMp0EU8hZ8f1KgBv/xSZqhBrQYE+AVoClgRD6ZAF/Ek3brB9u3l3hgeHxrPkdNHytlIRDyBAl3Ek3TvbrTflT203rBWQ3JO55TpFxHPoEAX8SS/+53RpqWVGQquFax70UU8mAJdxJNERhrt99+XGQr0N6aH/SWn7Dl2EXF/CnQRT+LnB0FBsHt3maHBLQcDMO+XeVX+sprLXcR8CnQRT9OyJezcWaa7Z7OeAKw6sMrZFYmIEyjQRTxNhw5w7hzk5tp1e3t5ExMUw6ZDm8ypS0QcSoEu4ml69zbacu5HT2ycSNbJLI6ePurcmkTE4RToIp4mJsZoV6woM1R6Hv3r3V9XyUtpcRYR16FAF/E0CQlGO2NGmaEhrYYAMHfHXCcWJCLOoEAX8TReXsZe+o4dsHy53VAdvzoE+gey6+guc2oTEYdRoIt4ovnzjfa118oMDWs9jF+P/cr+4/udXJSIOJICXcQTtWkDzZrBV19Bfr7d0NBWQwF4/ofnzahMRBxEgS7iqe66y2iXLLHrLj2P/uXOL51ckIg4kimB/vrrrxMXF0fbtm357W9/S0FBAceOHSM5OZkWLVqQnJxM7kX30E6YMIGYmBhatWrF4sWLzShZxP08/LDRrl9v1+3j5UP/5v05euYoVk3xJuIxnB7omZmZvPHGG6SkpLBt2zaKi4uZPXs2EydOpF+/fqSmptKvXz8mTpwIwI4dO5g9ezbbt29n0aJFjB8/nuLiYmeXLeJ+QkKMds6cMkM3Rd8EQG5BbpkxEXFPpuyhFxUVcebMGYqKijh9+jSNGzdm/vz5jBo1CoBRo0bx5ZdfAjB//nxGjhyJv78/UVFRxMTEsHbtWjPKFnEvXl4QF2fM6378uN1Qbd/aAOzN3WtGZSLiAE4P9CZNmvDnP/+ZiIgIwsLCqFu3LjfffDPZ2dmEhYUBEBYWxuHDhwFjj75p06a27cPDw8nMzHR22SLu6bHHjHbdOrvuhEYJAKzNrJpfjnXkXsR8Tg/03Nxc5s+fT1paGgcPHuTUqVN8/PHHFT6+vHN8lgqmp5o6dSqJiYkkJiaSk5NTZTWLuK2kJKPNyLDrbhvSFoBdR3Q/uoincHqgf//990RFRREcHIyvry/Dhw/np59+IjQ0lKysLACysrIIOX/+Lzw8nAMHDti2z8jIoHHjxuU+97hx40hJSSElJYXg4GDHvxkRV1e6Pvpc+5nhAvwDCK0dyre/fuv8mkTEIZwe6BEREaxZs4bTp09jtVpZsmQJsbGxDBkyhOnTpwMwffp0hg417pUdMmQIs2fP5uzZs6SlpZGamkqXLl2cXbaIe6pdG5o0gZSUMkNNApuQeixVV7qLeAgfZ79gUlISd9xxBx07dsTHx4cOHTowbtw4Tp48yYgRI3j//feJiIhg7vk9iri4OEaMGEGbNm3w8fFhypQpeHt7O7tsEfcVHw+LFkFJiXGh3Hk3Rd3EhqwNHMw/SJPAJiYWKCJVwWL10F/PExMTSSlnr0Sk2vm//4OHHoKtW6FtW1v3V7u+YsjsIXx/z/f0i+53TU+9YAEMHQobNhjLsIuIY10u2zRTnIini4012ktu92xa17h7pKqudBcRcynQRTxdx45Gm2s/iUy70HYAbM7e7OyKRMQBFOgini4wEHx9Ydkyu24vixdhdcL46cBPJhUmIlVJgS7i6SwWI9DPT9Z0seTmyRw4cYCCogITChORqqRAF6kORowwZou7ZB2ETmGdAPi/df9nRlUiUoUU6CLVwflplS+9MO7ehHsB+PyXz51dkYhUMQW6SHVw551G+7//2XUH+AfQLrQdGScyytlIRNyJAl2kOoiLM9r09DJD/aL6cTD/4HXNGOeZs1mIuBcFukh14OMDnTrBjBllhsLqhFFUUsS2w9tMKExEqooCXaS6CAoyLoo7d86uu2ezngC8uvpVM6oSkSqiQBepLgYNMtpLJphJCjeWWD1y+oizKxKRKqRAF6kuSpcUTk0tM9S/eX9NASvi5hToItVFRITRZmWVGfL19sWC5aqf0nL1m4iIgyjQRaqL0kVaXnyxzFBiWCI5p3M4ee6kk4sSkaqiQBepLoKCjPaXX8rcZxZSOwSAzYe0UIuIu1Kgi1Qn//wnFBbCoUN23Tc2uxHQUqoi7kyBLlKdxMQY7fLl9t1BRv/0zdOdXJCIVBUFukh1kpxstHv32nXX8KlBSO0QNmdv1nl0ETelQBepToKDoVYt2Fz2XPmLfYyL5X7c96OzqxKRKqBAF6lu/Pxgz54y3be0uAXQeXQRd6VAF6lu2raFDRvKXOkeFmAssbrjyI6rfkotziJiPgW6SHXTp4/Rbt9u1+1l8aJ+jfqcLjxtQlEicr0U6CLVTa9eRrtrV5mhLk26kH0y28kFiUhVUKCLVDfx8Ub7889lhkLrhLLu4DonFyQiVUGBLlLdNGxotNvKrn/u7+0PQFZ+2fneRcS1KdBFqhsvL0hIgG+/LXM1W3K0cZ/60TNHTShMRK6HAl2kOurUyWh32F/RHlTTmO/92JljV/Q0Wm1NxHUo0EWqo7Fjjfbbb+26SwP96GntoYu4GwW6SHWUmGi0S5bYdTeo1QC48j10EXEdCnSR6sjbGwID4fvv7bpte+g6hy7idhToItXV/fdDUREcPmzrqu1bG2+LN4dOHrrMhiLiihToItVVz55G+9FHti6LxUKgfyCLfl1kUlEicq0U6CLVVffuRvvOO3bd0fWjsejydRG3o0AXqa6Cg+E3v4Fff7W7H71reFd25OzAehUrrmhxFhHzKdBFqrPwcKM9ftzWFVwrGIDl6ctNKEhErpUCXaQ6K51gJvvCgix3t7sbgDUZa8yoSESukQJdpDoLCTHazZttXZH1IgE4cvqICQWJyLVSoItUZ7GxRvvrr7Yuby9vIupGcPDkQZOKEpFroUAXqc5Kz6Efs58ZrlGdRny580vn1yMi10yBLlLdNW4Mubl2Xf7e/hQUFXDy3MnLbqq720RchwJdpLoLCyszBezohNGA1kUXcScKdJHqrmFDyMmxu5m8SUATAH7Y94NZVYnIVVKgi1R3N90EZ85AZqatq1dkLwDWZa4zqyoRuUoKdJHqLjLSaJcutXXV8KkBQEFxgQkFici1UKCLVHc33WS0FwU6QKewTny/9/tyNhARV6RAF6nu6tWDoCCYPt3usLuftx8WdBm7iLtQoIsI/PWvRvvDhYvgekT04MjpI1e0SIsWZxExnwJdRGDsWKPduNHWFVYnjLPFZ8ktyK1gIxFxJQp0EYHAQKP96itbV5NA49a1bYe3mVGRiFwlBbqIGFO+deoEe/ZASQkAbYLbAPDxlo/NrExErpACXUQMw4dDURHs3QtA25C2AOSfyzezKhG5Qgp0ETEkJRntO+/Yum5oeoMmlxFxEwp0ETH07Wu0aWm2Li+L12UXaNHiLCKuQ4EuIgaLBfr1g0WLbF3dm3Yn+1T2Fd26JiLmUqCLyAX16sGpU3D8OGDcugaw8dDGy2wkIq5AgS4iFwwZYrTp6QDcEHEDAGsy1phUkIhcKQW6iFwQHW20W7cCkNAoAYBVB1aZVJCIXCkFuohc0Na4VY3NmwHw8fKhrn9d7aGLuAEFuohcUK+e0Z45Y+uKC4njyOkj5tQjIldMgS4i9tq0gTUX9sh7N+vNibMnOF14usJNdBG8iPkU6CJiLyAA1q+HwkIAgmoGAXAw/6CZVYlIJRToImKvdIKZffsAiA+NB+DQyUNmVSQiV0CBLiL2brrJaM/PGBdaOxTQqmsiru6ygb569Woeeugh2rVrR3BwMBEREdxyyy1MmTKF4+cnnhARDxMVZbRr1wIQUTcC4LLn0EXEfBUG+sCBA3nvvffo378/ixYtIisrix07dvDiiy9SUFDA0KFDWbBggTNrFRFniIw02u++A6BejXr4efuRejTVvJpEpFI+FQ3MmDGDhg0b2vXVqVOHjh070rFjRx5//HGOHNGtLCIex2KBhg0hNfX8lxZq+tRk/4n9JhcmIpdT4R76pWEO8PXXX1f6GBHxALfdBgcP2u5Hi6wXSUFRQZmHabU1EddxVRfF/eMf/3BUHSLiSiKM8+alc7pH1Y/i8KnD5tUjIpW6qkDXEooi1US3bkZ7finVkFohbDu8TT8DRFzYVQX6O++8UyUvmpeXxx133EHr1q2JjY1l9erVHDt2jOTkZFq0aEFycjK5ubm2x0+YMIGYmBhatWrF4sWLq6QGEbmMzp2N9tQpAAL9AwHIPpVtVkUiUokKA33lypVl+rp06WL7+4kTJ9i27druS/3jH//IgAED2LlzJ5s3byY2NpaJEyfSr18/UlNT6devHxMnTgRgx44dzJ49m+3bt7No0SLGjx9PcXHxNb2uiFyhgADjBPn57/HExokA/Jzxs5lVichlVBjon3/+Od27d+f5559n4cKFrF27lhUrVvDBBx9wzz33MGjQIM5ctIDDlTpx4gQrVqzgvvvuA8DPz4969eoxf/58Ro0aBcCoUaP48ssvAZg/fz4jR47E39+fqKgoYmJiWHv+/lgRcRCLBWrVgh9+AKBHRA8AtmRvMbMqEbmMCm9be/3118nNzeWzzz5j7ty5ZGVlUbNmTWJjY3nggQfo0aPHNb3g3r17CQ4O5t5772Xz5s106tSJyZMnk52dTVhYGABhYWEcPmxcgJOZmUnXrl1t24eHh5OZmXlNry0iV6F7d1i3DoAmgU2o6VOTY2eOlftQnVoXMV+FgQ5Qv359xo4dy9ixY6vsBYuKitiwYQNvvvkmSUlJ/PGPf7QdXi9PeRfhWCq4V2bq1KlMnToVgJycnKopWKS66tHDmFzmxAkIDKRJYBOyTmaZXZWIVMDpc7mHh4cTHh5OUlISAHfccQcbNmwgNDSUrCzjh0VWVhYhISG2xx84cMC2fUZGBo0bNy73uceNG0dKSgopKSkEBwc7+J2IeLjoaKOdMQOApoFNWZa+zMSCRORynB7ojRo1omnTpuzatQuAJUuW0KZNG4YMGcL06dMBmD59OkOHDgVgyJAhzJ49m7Nnz5KWlkZqaqrdxXki4iB33mm0X3wBGFe66150Edd12UPujvLmm29y1113ce7cOaKjo/nwww8pKSlhxIgRvP/++0RERDB37lwA4uLiGDFiBG3atMHHx4cpU6bg7e1tRtki1YuvrzGv+5IlYLUSHxLP/F3z2XlkJ60btja7OhG5RKV76HPnziU/Px+AF198keHDh7Nhw4bretGEhARSUlLYsmULX375JfXr16dBgwYsWbKE1NRUlixZQlBQkO3xTz/9NHv27GHXrl0MHDjwul5bRK7CzTcb7d699I/pD8DmQ5tNLEhEKlJpoL/wwgsEBASwcuVKFi9ezKhRo3jwwQedUZuImO3224124UJaBLUA4OvUry+zgYiYpdJALz28vXDhQh588EGGDh3KuXPnHF6YiLiA0ltGFywgtE4oAEdOX1hlUYuziLiOSgO9SZMmPPDAA8yZM4dbbrmFs2fPUlJS4ozaRMRsgYGQkADHjwNwc/ObyT2Te/ltRMQUlQb6nDlz6N+/P4sWLaJevXocO3aMV155xRm1iYgraN8ezt86GlQzqMLJZUTEXBVe5X7s2IVv2t69e9v6/P39SUxMdHhhIuIi6tWD7Gw4e5agGgp0EVdVYaB36tQJi8WC1Wq1taUsFgt79+51SoEiYrLStdEPHiSoZhC5BbmUWEvwsjh9GgsRuYwKAz0tLc2ZdYiIq2rZ0mg3bSKoURAl1hJOnD1BvRr1TC1LROxd0cQyubm5pKamUlBQYOvr2bOnw4oSERfSqZPRvv8+9V+8A4DcM7l2ga7FWUTMV2mgv/fee0yePJmMjAwSEhJYs2YN3bp1Y+nSpc6oT0TMFhYGtWvDwoUEvWos1HTszDGi6keZXJiIXKzSk2CTJ09m3bp1NGvWjGXLlrFx40YtfCJS3QwfDkBQobEPkH0q28xqRKQclQZ6jRo1qFGjBgBnz56ldevWtoVVRKSaOD8FbP2fjWlf9+XtM7MaESlHpYEeHh5OXl4ew4YNIzk5maFDh1a4fKmIeKjze+hNs43raHYe2WlmNSJSjkrPoX9xfunE5557jj59+nD8+HEGDBjg8MJExIXUqgUNGhDwv+XQF4qtxWZXJCKXqDTQ9+/fb/t7VJRxEcyhQ4eIKL03VUSqh8BALDt+ofntzckt0PSvIq6m0kC/9dZbbRPLFBQUkJaWRqtWrdi+fbsz6hMRV3HrrfDWW9T3j9BscSIuqNJA37p1q93XGzZs4J133nFYQSLiomJiAAiy1rQFulZbE3EdVz13Y8eOHVm3bp0jahERV9a5MwBBmUe14pqIC6p0D/21116z/b2kpIQNGzboPnSR6igpCYD6p0p0yF3EBVUa6Pn5+Rce7OPDrbfeyu233+7QokTEBXl7Q9euBB0/RG49Y4GWazjIJyIOUmmgP/vss86oQ0TcQWgo9fdsoCSihPyz+UBdsysSkfMqDPTBgwdjucwVLwsWLHBIQSLiwgIDCdp3DuD8YXcj0LU4i4j5Kgz0P//5zwDMmzePQ4cOcffddwMwa9YsIiMjnVKciLiYRx6h7u9mAMZscV5ogRYRV1FhoPfq1QuAZ555hhUrVtj6Bw8erKVTRaqrxERahrYBdrAjZwdtGWh2RSJyXqVXtOTk5LB3717b12lpaeTk5Di0KBFxXS1DYgF0pbuIi6n0orjXX3+d3r17Ex0dDUB6eromlhGpxvwC6hF8Co7mZ0NNs6sRkVKVBvqAAQNITU1l505jdaXWrVvj7+/v8MJExEXdeCOhP7/P4b1bIM7sYkSkVIWBvnTpUvr27cu8efPs+vfs2QPA8PPLKYpINTNwIA2Xwc/5u3jA7FpExKbCQP/hhx/o27cvX331VZkxi8WiQBeproKD8bVaOMiJ88uoeptdkYhwmUD/5z//CcB7772Ht7e+YUXkPIuFhGN+fBd9loyCnei4u4hrqPQq96ioKMaNG8eSJUuwavYIEQEGRd0MQPqJDSZXIiKlKg30Xbt2cdNNNzFlyhSioqL4wx/+wMqVK51Rm4i4qObh7QHI+WWWyZWISKlKA71mzZqMGDGCefPmsXHjRk6cOGGbdEZEqqfGv38CgML8gyZXIiKlrmippB9++IHx48fTsWNHCgoKmDNnjqPrEhEXZgkMJDbXh63eeyt/sIg4RaX3oUdFRZGQkMCIESN45ZVXqF27tjPqEhEX5+Plw37fk4AWZxFxBZUG+ubNmwkMDHRGLSLiRm7yjuGNGtsI5SDQ2OxyRKq9CgP95Zdf5sknn+Tvf/97ueNvvPGGw4oSEdfXLKoDxTnb6F17LvBHs8sRqfYqDPTYWGMBhk6dOjmtGBFxH7XadYIlM+gdMAsFuoj5Kgz0wYMHAzBq1CinFSMi7qN1044ABNfabnIlIgKVBLrFYqlwwwULFjikIBFxDyG1QwA4W/skvkeygDBzCxKp5ioM9D//+c8AzJs3j0OHDnH33XcDMGvWLCIjI51SnIi4ruDawQCsawy9P3kFhrxmckUi1VuFgV46ecwzzzzDihUrbP2DBw+mZ8+ejq9MRFxa/Rr1ASixQNCKLwAFuoiZKp1YJicnh717L0wekZaWRk5OjkOLEhHXZ7FYaODXmF/8mlIjKx2ysswuSaRaq/Q+9Ndff53evXsTHR0NQHp6Ou+8847DCxMR1+fvVYvdvsHAAfjiCxg/3uySRKqtSgN9wIABpKamsnPnTgBat26Nv7+/wwsTEddXw6s2+3yDjC+OHTO3GJFqrtJAB1i/fj3p6ekUFRWxefNmAH73u985tDARcX01vGtRHLKTooB6+GzaZHY5ItVapYF+zz33sGfPHhISEvD29gaMc2cKdBE5W3IGztXmXH0ffD7/3JjU/TK3u4qI41Qa6CkpKezYseOy96SLSPXUsnYie2ovJD/uZmrt3wVHjkBwsNlliVRLlV7l3rZtWw4dOuSMWkTEzdTzDYGALI4mDTA6/vc/cwsSqcYq3UM/cuQIbdq0oUuXLnYXw2mmOBGp4V0HgF3tWtIGYNUquOsuU2sSqa4qDfTnnnvOCWWIiDtqVqsNAFl+pyEoCN5+GyZNAj8/cwsTqYYqDfTSGeNERC7VyD8KgN0nNsDIkfB//weLF8P5xZ1ExHkqPIfeo0cPAAICAggMDLT9Kf1aRCTUPxKAHSd+hr/8xejcscO8gkSqsQr30FeuXAlAfn6+04oREfdS2ycQ8huRUy8TIiKgZk1Yt87sskSqpUqvchcRuawjrTldfP4X/8BAmDfP3HpEqikFuohcn7wodudvMP5+883G5DKaNU7E6RToInLNLBagxDhzV1BUAE88YQz8+9/mFSVSTSnQReT6ZHYG4MDxAxAfb/T99JOJBYlUT5UG+qVXuQcGBtK0aVNuu+02u3XSRaSaOt0QgC3ZW4yv77wTMjMhJcXEokSqn0oD/bHHHuOVV14hMzOTjIwMXn31VcaOHcvIkSMZM2aMM2oUEVeW0RWAnzN/Nr5+9FGjnT7dnHpEqqlKA33RokU88MADtj31cePG8c0333DnnXeSm5vrjBpFxJWdDAMgryDP+LprV2jYELZsMa8mkWqo0kD38vJizpw5lJSUUFJSwpw5c2xjWoFNRACa12lH9qnsCx2xsbBiBRQWmleUSDVTaaDPnDmTGTNmEBISQmhoKDNmzODjjz/mzJkzvPXWW86oUURcXH2/EL7a9dWFjr59jXbPHnMKEqmGKp3LPTo6mq+++qrcsdLpYUWkeqvtHYjFYsFqtRpH7rp0MQb27IHWrc0tTqSaqDTQc3JyePfdd0lPT6eoqMjW/8EHHzi0MBFxH7F1k/ghZx55BXnUr1kfOnQwBj7+GG691dziRKqJSgN96NCh3Hjjjdx00014e3s7oyYRcTMBPvUB2JGzgxsiboCwMGNe940bTa5MpPqoNNBPnz7NSy+95IxaRMRNta93IwDTNk0zAh3ghhvg+++NqWB1Aa2Iw1V6UdygQYP45ptvnFGLiLipiFqtsGBhTeaaC52xsUZ7+LA5RYlUM5UG+uTJkxk0aBA1a9bUeugiUi6LxUKPiB5sO7ztQmdXY8IZNF+FiFNUGuj5+fmUlJRw5swZTpw4QX5+PidOnHBGbSLi4i4+kh4XHAfA6cLTRkeDBkb74YdOrkqkeqow0Hfu3AnAhg0byv1zvYqLi+nQoQODBg0C4NixYyQnJ9OiRQuSk5PtZqGbMGECMTExtGrVisWLF1/3a4tI1esabuyR7z++3+i40TivjtZ8EHGKCi+Ke+2115g6dSqPP/54mTGLxcLSpUuv64UnT55MbGysbW9/4sSJ9OvXj6eeeoqJEycyceJEXnrpJXbs2MHs2bPZvn07Bw8e5KabbmL37t264l7ExYTUDgFga/ZWWjdsDbVqQZ8+sHu3yZWJVA8VBvrUqVMBWLZsWZW/aEZGBgsXLuTpp5/mtddeA2D+/PksX74cgFGjRtG7d29eeukl5s+fz8iRI/H39ycqKoqYmBjWrl1Lt27dqrwuEbl2HcKMe8935Oy40FmvHpz/vhYRx6r0HHr79u2ZMGECe6pwCsdHH32Ul19+GS+vCy+fnZ1NWJixyENYWBiHz18Zm5mZSdOmTW2PCw8PJzMzs8pqEZHrZ7Ve2EO3m9O9RQtjMC3NpMpEqo9KA33BggV4e3szYsQIOnfuzKuvvsr+/fuv+QW//vprQkJC6NSp0xU93mq1lumraFGYqVOnkpiYSGJiIjk5Oddco4hcPS+LF60btmbdwXUXOrt3N1rN6S7icJUGerNmzXjyySdZv349n3zyCVu2bCEqKuqaX3DVqlUsWLCAyMhIRo4cydKlS7n77rsJDQ0lKysLgKysLEJCjN/2w8PDOXDggG37jIwMGjduXO5zjxs3jpSUFFJSUggODr7mGkXk2tTxq8OmQ5sudLRpY7S//mpKPSLVSaWBDpCens7LL7/MyJEj2blzJy+//PI1v+CECRPIyMggPT2d2bNn07dvXz7++GOGDBnC9OnTAZg+fTpDhw4FYMiQIcyePZuzZ8+SlpZGamoqXUoXfhARl9KhUQeKSoooKCowOkp/+b7ol3IRcYxKp35NSkqisLCQ3/zmN8ydO5fo6GiHFPLUU08xYsQI3n//fSIiIpg7dy4AcXFxjBgxgjZt2uDj48OUKVN0hbuIi+oa3pV3N7xL9slsmtVrBrVrG1e7Hz9udmkiHq/SQJ8+fTqtHbT8Ye/evenduzcADRo0YMmSJeU+7umnn+bpp592SA0iUnVKL4zbnL3ZCHSAiAj44QcTqxKpHioM9I8//pi7776bb775pty53B977DGHFiYi7qd1Q+OX/+/2fMeQVkOMzsBA3Ysu4gQVBvqpU6cAY+rXS1V0lbmIVG/N6zcH4Nfciy6Ca9EC1q7VqmsiDlZhoD/wwAMAPPvss2XGJk2a5LCCRMR9WSwWuoZ35bs9313oTEqCmTPh0CFjnXQRcYgrusr9UqWzu4mIXKpZ3WYUW4vJPXN+PYbmxl47779vXlEi1cA1BXp5k72ISPVT3hH0Ya2HAfDd3vN76QMHGu1//uOcokSqqWsKdJ1DF5GK9I7sDcC0TdOMDovFmDEuLw/Or+IoIlWvwnPoAQEB5Qa31WrlzJkzDi1KRNxXozqN8LJ48euxiy6Me/FF6NvXWKjFQbfBilR3Fe6h5+fnc+LEiTJ/8vPzKSoqcmaNIuLiLj0L90iXR0g9lkpRyfmfFT17Gu2iRc4tTKQauaZD7iIil1O/Zn0AMk5kGB3e3tCwISxcWDb9RaRKKNBFpMq1DWkLwJqMNRc6BwyAoiLYu9ekqkQ8mwJdRKpcUpMkAPIK8i50Dh9utNu2Ob8gkWpAgS4iVS64trF8cc6pnAudHTsa7Zw5JlQk4vkU6CJS5fy8/Qj0D+TH/T9e6GzWzDiP/sknoDtlRKqcAl1EHMLb4s36rPX2nWPHGu1335XdQESuiwJdRBxibMexHDtz7MKV7gB/+pPRvvOOOUWJeDAFuog4RK/IXgBMWjPpQmdwsHEL28aN5hQl4sEU6CLiEANjjDncV2esth8YMAAKCkyoSMSzKdBF5JpdblkHi8VC64at7a90B0hMNOZ114yTIlVKgS4iDtM3si+px1LtO0NCjNnijh41pygRD6VAFxGH8fX2BeDkuZMXOkNDjXb9+nK2EJFrpUAXketW0fTsXZp0AWBj1kUXwcXHG+3q1eVsISLXSoEuIg7TuqGxVOoXO7+40BkTY7THjplQkYjnUqCLiMMkNEoA4Id9P1zo9PIy9tK/+cacokQ8lAJdRBzGy+JFr2a92JC1wX6gVi3toYtUMQW6iDhUXHAcAMcLjl/o7N0bTpyAI0fMKUrEAynQRcShSg+7292+1q2b0T7wgPMLEvFQCnQRcajY4FgAMk9kXugcOhR8fOCLLyrYSkSulgJdRBwqtLZx33mZCWaGDzfud8vPN6EqEc+jQBcRh4qsFwnAsvRl9gPDhhntokVOrUfEUynQRcShfL19iaoXRUHRJQuy9OljtC+95PyiRDyQAl1ErtnlFme5WIewDmVvXWvUCJo1gw0byt9IRK6KAl1EHM7Hy4dT505hvXSO2ORk4zz64cPmFCbiQRToIuJwXZt0pbCkkLyCPPuBW2812o8/dnpNIp5GgS4i162ixVlKhdYxrnQ/fOqSPfGBA412/34HVCVSvSjQRcThQmqHAOUEur8/tGgBhw6ZUJWIZ1Ggi4jDlQb65uzNZQcbNYLly51bkIgHUqCLiMPFBBlLpi5PX1520Nu78mP2IlIpBbqIOFwt31p4WbzYdnhb2cHu3eHoUSgpcX5hIh5EgS4iTjGk1RB2Hd1FUUmR/UBwMBQXw/Hj5W8oIldEgS4iTtGhUQcAsk9m2w80bGi0mZmIyLVToIuIU8Q2NFZd2310t/1AUJDR/vKLkysS8SwKdBFxiqj6UQAcO3PMfqBTJ6P9+WcnVyTiWRToIuIU4YHhACxJW2I/EBoKXl7w6acmVCXiORToIuIUpeuil3sv+k03QUaGkysS8SwKdBG5Zle62prxWAsxQTGcLjxddrBlS6PNz6+awkSqIQW6iDhN72a92XRoU9lb1zp3NtqcHOcXJeIhFOgict2udKK36PrRACzYtcB+IMSYGlbLqIpcOwW6iDjN7xN/D5QzBWxpoH//vXMLEvEgCnQRcZr6NesDsD1nu/1AfLzR7tzp5IpEPIcCXUScqkVQC06dO2Xf6esLzZvDtnLmeheRK6JAFxGn6ta0G78cKWdWuMhIKChwej0inkKBLiJOFeAXwImzJzh57qT9QOvWkJVlTlEiHkCBLiJO1SnMmOp1beZa+4GGDeHECeOPiFw1BbqIOFXX8K4ALNl7yRSwERFGu3q1kysS8QwKdBFxqthgY9W11RmXBHeXLkarZVRFrokCXUScrnFAY1bsW2Hf2aKF0X79tfMLEvEACnQRcbr7O9xPsbXY/vY1f38IDIS0NPMKE3FjCnQRuWZXszjLxULrGCuv/bj/R/uBzp1h06Yrn0tWRGwU6CJy3a42f2+MuBGA9Lx0+4EOHYx2//7rL0qkmlGgi4jTRdaLBCD/7CXLpfbrZ7Q/XrLnLiKVUqCLiNPV9qsNUHZymcREo9UUsCJXTYEuIk7nZfGitm9t8s9dsofesCHUrAmffWZOYSJuTIEuIqYI8A/gjZ/fKDvQuTPs2QMnT5YdE5EKKdBFxBTBtYIpthaXvTBuwACj1e1rIldFgS4ippjQbwIAvx771X7ghhuMNifHyRWJuDcFuoiYomWDlgCsy1xnPxASYrS7djm5IhH3pkAXEVNE148G4Jtfv7EfaNzYaPfudXJFIu5NgS4ipvD28iY8MJyV+1faDwQGQo0a2kMXuUoKdBExTb8oYyKZI6eP2A/ExcFXX5lQkYj7cnqgHzhwgD59+hAbG0tcXByTJ08G4NixYyQnJ9OiRQuSk5PJzc21bTNhwgRiYmJo1aoVixcvdnbJIuIg/Zv3ByDlYIr9QEyM0R444OSKRNyX0wPdx8eH//znP/zyyy+sWbOGKVOmsGPHDiZOnEi/fv1ITU2lX79+TJw4EYAdO3Ywe/Zstm/fzqJFixg/fjzFxcXOLltEynGti7OUSmxszAxXZg999GijTbkk6EWkQk4P9LCwMDp27AhAQEAAsbGxZGZmMn/+fEaNGgXAqFGj+PLLLwGYP38+I0eOxN/fn6ioKGJiYli7dq2zyxaRy7jWxdEa1WkEwKr9q+wH2rY1Ws0YJ3LFTD2Hnp6ezsaNG0lKSiI7O5uwsDDACP3Dhw8DkJmZSdOmTW3bhIeHk5mZaUq9IlK16vjVAWB7znb7gdIr3T/5xMkVibgv0wL95MmT3H777UyaNInAwMAKH2ct51d/SwXH+aZOnUpiYiKJiYnkaFIKEZdnsVjo3LgzPx34yX7Aywsee8z4+1tvOb8wETdkSqAXFhZy++23c9dddzF8+HAAQkNDycrKAiArK4uQ85NLhIeHc+CiC2MyMjJoXPrb+yXGjRtHSkoKKSkpBAcHO/hdiEhV6N60O8XWYvbmXnLf+XPPGa1OsYlcEacHutVq5b777iM2NpbHSn8DB4YMGcL06dMBmD59OkOHDrX1z549m7Nnz5KWlkZqaipdunRxdtki4iBDWg0BYMbmGfYDAQGQlARbt5pQlYj78XH2C65atYoZM2YQHx9PQkICAP/+97956qmnGDFiBO+//z4RERHMnTsXgLi4OEaMGEGbNm3w8fFhypQpeHt7O7tsEXGQ3pG9ATh86nDZwRo1jJXXRKRSTg/0Hj16lHteHGDJkiXl9j/99NM8/fTTjixLREziZfEiLjiOdQfXlR1MSoIffoDTp6FWLecXJ+JGNFOciLiEzPxy7l4pvX3tqaecW4yIG1Kgi4jp+kb15diZY2UH7rrLaD/80LkFibghBbqImC48MJyCogL2H99vP+DlBbffDidPguafELksBbqImC4uOA6AXUfKWWGtdBpYBbrIZSnQRcR0bUOMc+WLfl1UdjA01Gg3bHBiRSLuR4EuIqZrVq8ZAAVFBWUHW7Uy2vx8J1Yk4n4U6CJyzUpnYb7WxVku1ia4Deuz1pcdCAgw7kdfvfr6X0TEgynQRcQl1PSpyc+ZP5cdsFjA2xu0bLLIZSnQRcQldGliTOlcZm10gG7d4HA5M8mJiI0CXURcQt+ovgCsyyxnxriQEFhXTr+I2CjQRcQlJDRKAODH/T+WHQwIMA65b97s3KJE3IgCXURcQkxQDIH+gXy9++uyg/ffb7Q/lhP2IgIo0EXEhYQHhpN9KrvsQLt2Rrt8uVPrEXEnCnQRcRm3xNzC4VOH2X10t/2Anx8EB2tyGZHLUKCLiMsovTDux33lHFpv185YRlVEyqVAFxGX0S+6H0D596PfeCNkZxsLtYhIGQp0EXEZft5++Hr58u6GdykqKbIfDAszWs0YJ1IuBbqIuJRhrYcB8Nrq1+wHbrjBaPftc25BIm5CgS4iLmXGbTMA+HDTh/YDLVsa7cqVTq5IxD0o0EXkmlXl4iyl/H38uSn6JnYe2Yn14if29TUWafn886p9QREPoUAXEZeTEJoAQG5B7iUDCcZFcd995/SaRFydAl1EXE7nJp0B+OKXL+wHpk832q1bnVyRiOtToIuIy0mOTgZg1rZZ9gMtWoC/PyxYYEJVIq5NgS4iLqd+zfo0qtOIlfsvuQDOYoE6dWDFCsjNLX9jkWpKgS4iLql/8/6cLT5LXkGe/cCkSUY7e7azSxJxaQp0EXFJg1oOAuDRRY/aD4wYYbRz5ji3IBEXp0AXEZc0MGYgANM3T7e/fc3PD2rV0iF3kUso0EXEJdX2q80zPZ8B4K21b9kP9u8P+fkmVCXiuhToIuKynuj+BAA/7r9k9bWICGMKWE0wI2KjQBcRlxXgH0C70HZszt5sPxAcDMXFcOiQOYWJuCAFuoi4tFYNWrH76G5Onrto2dTWrY02Pd2UmkRckQJdRFxa78jeALy++vULnc2bG+0HHzi/IBEXpUAXkevmyFPZDyY+CMBra167cLV7fLzRbtjguBcWcTMKdBG5ZqWrrTn2NSz0iexDXkEee3P3Gp3e3jBqlBHoBQWOL0LEDSjQRcTl/anrnwDYdXTXhc7S9dE3by5nC5HqR4EuIi6vY1hHAF5fc9F59J49jTYz04SKRFyPAl1EXF6TwCYAHMw/eKGzdA/9k09MqEjE9SjQRcQtPNLlEXbk7GBL9hajIzjYaA8fNq8oEReiQBcRtzCs9TAAVuxbYXRYLPCb38CPP0JJiXmFibgIBbqIuIWezXri4+XD35f+/UJnRITR/vqrOUWJuBAFuoi4BW8vbzo06sDxs8f56cBPRmdystGuWmVeYSIuQoEuIm7jtf6vATDg4wFGR/fuRrtihUkVibgOBbqIuI0eET3o1awX+efyyT+bDwEBxp9p02DLFrPLEzGVAl1E3Mp9He4DYN/xfUbH9OlG+/77JlUk4hoU6CLiVhoHNAbg022fGh3DhhntG2+YU5CIi1Cgi8h1c+TiLJdKCk8C4MUfX+Tr3V8bt6/9v/9nDGp9dKnGFOgics2csTjLper41eGH0T8AcMecO4zOm2822r17nV+QiItQoIuI2+nZrCedG3fmbPFZPtr8kZZTFUGBLiJuavno5QBszd4KsbFG5w8/mFeQiMkU6CLilmr51iK6fjT/Xf9frDVqQJ06sG+f2WWJmEaBLiJuK6xOGCfPnWRp2lLo1w/WrTO7JBHTKNBFxG1NGzYNgHfWv3PhsPvu3eYVJGIiBbqIuK2YoBh8vHyYu2Mup3qfnwb2jjvMLUrEJAp0EXFrD3d5GID/q7MTGjeGrVvh1CmTqxJxPgW6iLi1/9z8HwCWpC2BV14xOl991cSKRMyhQBcRt2axWGhevzmL9yym6PbbjM6vvjK3KBETKNBFxO3d0cY4b377FyOhTx9ITze3IBETKNBFxO090/MZABbsWsD67lFw9CgUFJhclYhzKdBF5Lo5c3GW8tT2q827g98FINH3A074o4VapNpRoIvINTNjcZaK3N/xfp7r9RwAieOg6MhhcwsScTIFuoh4jH/0+gd+Fl9SG4Dfwq4cPqVQl+pDgS4iHsNisXD6D1l0PQBWrMz7ZZ7ZJYk4jQJdRDyKd1ADvpth/P3BhQ/SZ3of9uVp0RbxfAp0EfE4dfoO4N11jWjdsDXL05czbdM0s0sScTgFuoh4niZNuH/hIXY8uB0vixfP/fAcr/70KsUlxWZXJuIwCnQR8TzNmwNg2bOHyQMmA/DEd08QMCGADVkbzKxMxGEU6CLiebp2NdrXXuMPXf7A8aeO0zSwKWeKztBpaiceX/y4ufWJOICP2QWIiFS5Pn2M9uefAQj0DyT14VRW7FvBzR/fzGtrXuPX3F+xYMGKlVtb3MrYjmOxuNKN9SJXSXvoIuKZevSAjRttU8D6+/iT3DyZeSPm0T60Pfvy9rEndw8Ldi3gga8fYMaWGSYXLHJ93CbQFy1aRKtWrYiJiWHixIlmlyMirm7YMKMtXVL1vNtib2PT7zex6feb2PrgVn59+FcARn05iiGzhjB3+1zOFJ5xcrEi189itZo9C3PliouLadmyJd999x3h4eF07tyZWbNm0aZNmwq3SUxMJCUlxYlVilQ/a9ZAt27w7bcwYIDZ1VyiqAh8fSEqCvbuvexDn132LC//9DIFRRcWdPlD5z+U+1iLxcK4TuOIC47TIXpxustlm1ucQ1+7di0xMTFER0cDMHLkSObPn3/ZQBcR5/n6a1dcsdSH4c0SCUlLYeWINyjx8avwkf0Iow//4aDvEd4MmcsvNdL4aM2HZR53wvsUAG+ufZPA4tqEFTawjRVTwsDj3fC1usWPVXGSwFoNePY/bzrltdzif15mZiZNmza1fR0eHs7P5y92udjUqVOZOnUqADk5OU6rT6S6CgkBb2+YMsXsSsr3I48yk7vpMfePV7zN/7vMmBVYHANvJIGX9RRgBHxaPdgVAm/WmHs95YoHijrqy7Mo0G3KOytQ3qGucePGMW7cOMA4LCEijhUdDceOwenTZldSkbvIzh2IpfBclT1jR2BaOf155/I4V1J1ryOewcfX13mv5bRXug7h4eEcOHDA9nVGRgaNGzc2sSIRKRUYaPxxWY2CnPIyITRyyuuIVMQtrnLv3LkzqamppKWlce7cOWbPns2QIUPMLktERMRluMUeuo+PD2+99Rb9+/enuLiYMWPGEBcXZ3ZZIiIiLsMtAh3glltu4ZZbbjG7DBEREZfkFofcRURE5PIU6CIiIh5AgS4iIuIBFOgiIiIeQIEuIiLiARToIiIiHkCBLiIi4gEU6CIiIh5AgS4iIuIBFOgiIiIeQIEuIiLiARToIiIiHkCBLiIi4gEU6CIiIh7AYrVarWYX4QgNGzYkMjKyyp4vJyeH4ODgKns+s+h9uA5PeA+g9+Fq9D5cS1W/j/T0dI4cOVLumMcGelVLTEwkJSXF7DKum96H6/CE9wB6H65G78O1OPN96JC7iIiIB1Cgi4iIeAAF+hUaN26c2SVUCb0P1+EJ7wH0PlyN3odrceb70Dl0ERERD6A9dBEREQ+gQL/IokWLaNWqFTExMUycOLHMuNVq5ZFHHiEmJoZ27dqxYcMGE6q8vAMHDtCnTx9iY2OJi4tj8uTJZR6zfPly6tatS0JCAgkJCTz//PMmVFq5yMhI4uPjSUhIIDExscy4O3weu3btsv07JyQkEBgYyKRJk+we46qfx5gxYwgJCaFt27a2vmPHjpGcnEyLFi1ITk4mNze33G0r+15ypvLexxNPPEHr1q1p164dt912G3l5eeVuW9n/QWcq730899xzNGnSxPZ/55tvvil3W1f/PO68807be4iMjCQhIaHcbV3p86joZ62p3yNWsVqtVmtRUZE1OjraumfPHuvZs2et7dq1s27fvt3uMQsXLrQOGDDAWlJSYl29erW1S5cuJlVbsYMHD1rXr19vtVqt1hMnTlhbtGhR5n0sW7bMeuutt5pR3lVp1qyZNScnp8Jxd/g8LlZUVGQNDQ21pqen2/W76ufxww8/WNevX2+Ni4uz9T3xxBPWCRMmWK1Wq3XChAnWJ598ssx2V/K95EzlvY/FixdbCwsLrVar1frkk0+W+z6s1sr/DzpTee/j2Weftb7yyiuX3c4dPo+LPfbYY9Z//vOf5Y650udR0c9aM79HtId+3tq1a4mJiSE6Oho/Pz9GjhzJ/Pnz7R4zf/58fve732GxWOjatSt5eXlkZWWZVHH5wsLC6NixIwABAQHExsaSmZlpclWO4Q6fx8WWLFlC8+bNadasmdmlXJGePXsSFBRk1zd//nxGjRoFwKhRo/jyyy/LbHcl30vOVN77uPnmm/Hx8QGga9euZGRkmFHaVSnvfVwJd/g8SlmtVubMmcNvf/tbJ1d19Sr6WWvm94gC/bzMzEyaNm1q+zo8PLxMEF7JY1xJeno6GzduJCkpqczY6tWrad++PQMHDmT79u0mVFc5i8XCzTffTKdOnZg6dWqZcXf7PGbPnl3hDyp3+DwAsrOzCQsLA4wfaIcPHy7zGHf7XD744AMGDhxY7lhl/wddwVtvvUW7du0YM2ZMuYd33enz+PHHHwkNDaVFixbljrvq53Hxz1ozv0d8rvsZPIS1nIv9LRbLVT/GVZw8eZLbb7+dSZMmERgYaDfWsWNH9u3bR506dfjmm28YNmwYqampJlVasVWrVtG4cWMOHz5McnIyrVu3pmfPnrZxd/o8zp07x4IFC5gwYUKZMXf5PK6UO30u//rXv/Dx8eGuu+4qd7yy/4Nme/DBB3nmmWewWCw888wzPP7443zwwQd2j3Gnz2PWrFmX3Tt3xc/jcj9rK+Koz0R76OeFh4dz4MAB29cZGRk0btz4qh/jCgoLC7n99tu56667GD58eJnxwMBA6tSpA8Att9xCYWFhhXMDm6n03zYkJITbbruNtWvX2o27y+cB8O2339KxY0dCQ0PLjLnL5wEQGhpqO62RlZVFSEhImce4y+cyffp0vv76a2bOnFnhD9PK/g+aLTQ0FG9vb7y8vBg7dmy59bnL51FUVMS8efO48847K3yMq30e5f2sNfN7RIF+XufOnUlNTSUtLY1z584xe/ZshgwZYveYIUOG8NFHH2G1WlmzZg1169a1HVpxFVarlfvuu4/Y2Fgee+yxch9z6NAh22+Ia9eupaSkhAYNGjizzEqdOnWK/Px829//97//2V0VC+7xeZS63J6HO3wepYYMGcL06dMBIxCHDh1a5jFX8r1ktkWLFvHSSy+xYMECatWqVe5jruT/oNkuvmbkiy++KLc+d/g8AL7//ntat25NeHh4ueOu9nlU9LPW1O+R676szoMsXLjQ2qJFC2t0dLT1xRdftFqtVuvbb79tffvtt61Wq9VaUlJiHT9+vDU6Otratm1b67p168wst1w//vijFbDGx8db27dvb23fvr114cKFdu/jzTfftLZp08barl07a1JSknXVqlUmV13Wnj17rO3atbO2a9fO2qZNG7f9PKxWq/XUqVPWoKAga15enq3PHT6PkSNHWhs1amT18fGxNmnSxPree+9Zjxw5Yu3bt681JibG2rdvX+vRo0etVqvVmpmZaR04cKBt2/K+l8xS3vto3ry5NTw83PY98sADD1itVvv3UdH/QbOU9z7uvvtua9u2ba3x8fHWwYMHWw8ePGi1Wt3v87BardZRo0bZvidKufLnUdHPWjO/RzRTnIiIiAfQIXcREREPoEAXERHxAAp0ERERD6BAFxER8QAKdBEREQ+gQBcREfEACnSRauTo0aO2ZSobNWpkW3qzTp06jB8/3iGvOWnSJD766KMy/enp6baJQbZu3cro0aMd8voi1YXmchepRho0aMCmTZsAYy3tOnXq8Oc//9lhr1dUVMQHH3xQ6Vr18fHxZGRksH//fiIiIhxWj4gn0x66iLB8+XIGDRoEGEE/atQobr75ZiIjI5k3bx5PPvkk8fHxDBgwgMLCQgDWr19Pr1696NSpE/379y936dqlS5fSsWNH21Kl69evp3379nTr1o0pU6bYPXbw4MHMnj3bwe9UxHMp0EWkjD179rBw4ULmz5/P3XffTZ8+fdi6dSs1a9Zk4cKFFBYW8vDDD/PZZ5+xfv16xowZw9NPP13meVatWkWnTp1sX99777288cYbrF69usxjExMT+fHHHx36vkQ8mQ65i0gZAwcOxNfXl/j4eIqLixkwYABgHBpPT09n165dbNu2jeTkZACKi4vLXRgnKyuL2NhYAI4fP05eXh69evUC4J577uHbb7+1PTYkJISDBw86+q2JeCwFuoiU4e/vD4CXlxe+vr625UW9vLwoKirCarUSFxdX7p72xWrWrElBQQFgrE51uTWfCwoKqFmzZhW9A5HqR4fcReSqtWrVipycHFugFxYWsn379jKPi42N5ddffwWgXr161K1bl5UrVwIwc+ZMu8fu3r3b5ZYnFXEnCnQRuWp+fn589tln/OUvf6F9+/YkJCTw008/lXncwIEDWbFihe3rDz/8kIceeohu3bqV2RtftmwZt956q8NrF/FUWj5VRBzqtttu4+WXX6ZFixYVPubs2bP06tWLlStX2q6IF5Gro0AXEYfatWsX2dnZ9OzZs8LHpKamkpmZSe/evZ1XmIiHUaCLiIh4AJ1DFxER8QAKdBEREQ+gQBcREfEACnQREREPoEAXERHxAP8fagBgBcwHACMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 566.929x566.929 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Specify a fixed time mortality calculator\n",
    "config.set('BIO_MODEL', 'mortality_calculator', 'fixed_time')\n",
    "\n",
    "# 1) Fixed time scenario\n",
    "age_of_death_in_days = 10.\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'initialisation_method', 'common_value')\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'common_value', str(age_of_death_in_days))\n",
    "n_alive_common_value = run(config)\n",
    "           \n",
    "# 2) Uniform random\n",
    "minimum_bound = 8.\n",
    "maximum_bound = 12.\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'initialisation_method', 'uniform_random')\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'minimum_bound', str(minimum_bound))\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'maximum_bound', str(maximum_bound))\n",
    "n_alive_uniform_random = run(config)\n",
    "\n",
    "# 2) Gaussian random\n",
    "mean = 10.\n",
    "standard_deviation = 1.\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'initialisation_method', 'gaussian_random')\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'mean', str(mean))\n",
    "config.set('FIXED_TIME_MORTALITY_CALCULATOR', 'standard_deviation', str(standard_deviation))\n",
    "n_alive_gaussian_random = run(config)\n",
    "\n",
    "# Set the bio time step\n",
    "config.set('NUMERICS', 'time_step_bio', str(time_step))\n",
    "    \n",
    "    \n",
    "# Plot\n",
    "font_size = 10\n",
    "fig, ax = create_figure(figure_size=(20, 20), font_size=font_size)\n",
    "plt.plot(times/seconds_per_day, n_alive_common_value, 'b', label='common_value')\n",
    "plt.plot(times/seconds_per_day, n_alive_uniform_random, 'r', label='uniform_random')\n",
    "plt.plot(times/seconds_per_day, n_alive_gaussian_random, 'g', label='gaussian_random')\n",
    "plt.ylabel('Living individuals (-)', fontsize=font_size)\n",
    "plt.xlabel('Time (d)', fontsize=font_size)\n",
    "\n",
    "# Add legend\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ProabilisticMortalityCalculator\n",
    "\n",
    "The mortality calculator kills particles at a rate $\\mu_{D} \\Delta t$, where $\\mu_{D}$ is a fixed mortality rate which is set in the run configuraiton file and $\\Delta t$ is the model time step for biological processes. The model computes a uniform random deviate in the range (0, 1). If the number is less than the computed death rate, the particle is killed. Below, we create a population of $1000$ individuals. We apply a death rate of $0.1$ per day and use a time step of $100$ seconds. The model is run forward for $20$ days and the number of living individuals plotted as a function of time. The result is compared with a simple analytical solution of exponential decay."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f456c4712e0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAHZCAYAAABw5MliAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABThklEQVR4nO3de3yO9R/H8de9mTmOHBtzGFsMQyxUzppDTuHnFFoh/UIphH6lVORQP4ciJYdItRCmiMgxhzSkEM1hv2zmFHKI2eH6/fHVIps5bPe1+977+Xjsse26rvu+P5d729v3ur4Hh2VZFiIiIuLSPOwuQERERO6cAl1ERMQNKNBFRETcgAJdRETEDSjQRURE3IACXURExA3ksLuAzFKkSBHKli1rdxkiIiIZJjo6mpMnT6a6z20DvWzZskRGRtpdhoiISIYJCQlJc58uuYuIiLgBBbqIiIgbUKCLiIi4Abe9hy4iIllHQkICMTExXLp0ye5SXEKuXLnw8/PDy8vrph+jQBcRkUwXExND/vz5KVu2LA6Hw+5ysjTLsvj999+JiYnB39//ph+nS+4iIpLpLl26ROHChRXmN8HhcFC4cOFbvpqhQBcREadQmN+82/m3UqCLiIg4wZIlSxgzZkymPb/uoYuIiGSyxMRE2rRpQ5s2bTLtNdRCFxGRbCE6OpqgoCCefPJJKleuTNOmTbl48SINGzZMmVn05MmTKdOGf/TRRzzyyCO0bt0af39/Jk+ezPjx47n33nupU6cOp06dAuDAgQM0b96cmjVrUq9ePfbu3QvA448/zsCBA2nUqBFDhw7lo48+on///gAcO3aMdu3aUa1aNapVq8amTZvu+PzUQhcREad67jn48ceMfc7q1WHixPSPi4qK4rPPPuPDDz+kU6dOfPHFFzc8fteuXezYsYNLly4REBDA2LFj2bFjB88//zxz5szhueeeo0+fPrz//vsEBgby/fff07dvX1avXg3Ar7/+yqpVq/D09OSjjz5Ked5nn32WBg0asGjRIpKSkjh//vztn/wVCnQREck2/P39qV69OgA1a9YkOjr6hsc3atSI/Pnzkz9/fgoUKEDr1q0BCA4O5qeffuL8+fNs2rSJjh07pjwmPj4+5euOHTvi6el53fOuXr2aOXPmAODp6UmBAgXu8MwU6CIi4mQ305LOLN7e3ilfe3p6cvHiRXLkyEFycjLAdUPFrj7ew8Mj5XsPDw8SExNJTk6mYMGC/JjGJYe8efNm8BmkTffQRUQkWytbtizbtm0DYMGCBbf0WB8fH/z9/Zk/fz5gJoXZuXNnuo9r0qQJU6dOBSApKYmzZ8/eYtXXy7RA79mzJ8WKFaNKlSop206dOkVoaCiBgYGEhoZy+vTplH2jR48mICCAChUqsGLFipTt27ZtIzg4mICAAJ599lksy8qskkVEJBsaPHgwU6dO5YEHHkhzrfEb+eSTT5gxYwbVqlWjcuXKREREpPuYSZMmsWbNGoKDg6lZsya7d+++ndKvZWWSdevWWdu2bbMqV66csu2FF16wRo8ebVmWZY0ePdoaMmSIZVmWtXv3bqtq1arWpUuXrIMHD1rlypWzEhMTLcuyrPvuu8/atGmTlZycbDVv3txatmzZTb1+zZo1M/iMRETkdu3Zs8fuElxOav9mN8q2TGuh169fn0KFCl2zLSIigrCwMADCwsJYvHhxyvYuXbrg7e2Nv78/AQEBbN26lbi4OM6ePcv999+Pw+HgscceS3mMs1yIOc3aB//D4eUZ8L8nERGRTOLUe+jHjh3D19cXAF9fX44fPw5AbGwspUqVSjnOz8+P2NhYYmNj8fPzu257WqZNm0ZISAghISGcOHEiQ2q+9GcytTZNJO6F/2bI84mIiGSGLNEpzkrlvrjD4Uhze1r69OlDZGQkkZGRFC1aNENqK3xPYRb49KT6rrlsWRSXIc8pIiKS0Zwa6MWLFycuzoRiXFwcxYoVA0zL+/DhwynHxcTEUKJECfz8/IiJibluu7NdfOp5cpDIkRffdfpri4iI3AynBnqbNm2YPXs2ALNnz6Zt27Yp28PDw4mPj+fQoUNERUVRq1YtfH19yZ8/P1u2bMGyLObMmZPyGGd6alx5VuRpR8N975P0x53P5iMiIpLRMi3Qu3btyv3338++ffvw8/NjxowZDBs2jJUrVxIYGMjKlSsZNmwYAJUrV6ZTp05UqlSJ5s2bM2XKlJSZdaZOnUrv3r0JCAigfPnytGjRIrNKvqHtjQZTiNNMrD7LltcXERG5EYeV2o1qNxASEpIy2X5GSE6GzZ4P4kscM4ZGMWrM9VP5iYhI6n755ReCgoLsLiNDfPTRR0RGRjJ58uQbHtO0adOU28S9e/dm4MCBVKpU6aZfJ7V/sxtlW5boFOcKPDyg/JRBlOMQ+8YusrscERHJwj766COOHDmS8v306dNvKcxvhwL9Ftz9VFtOFijPYN4mNsYtL2yIiLi1Rx55hJo1a1K5cmWmTZsGQL58+XjppZeoVq0aderU4dixYwB8+eWX1K5dm3vvvZeHHnooZftfzp07h7+/PwkJCQCcPXuWsmXLMn/+fCIjI+nWrRvVq1e/bonW5cuXU6NGDapVq0aTJk0y7Ny0OMut8PQkptNA6nzYj83hmyg5+EG7KxIRcT02rp86c+ZMChUqxMWLF7nvvvvo0KEDFy5coE6dOowaNYohQ4bw4Ycf8vLLL1O3bl22bNmCw+Fg+vTpjBs3jv/+9+85SfLnz0/Dhg1ZunQpjzzyCOHh4XTo0IGOHTsyZcoU3n77bUJCQq55/RMnTvDkk0+yfv16/P39U9ZUzwhqod+ifP0f53cKUTJcE82IiLiad955J6UlfvjwYaKiosiZMyetWrUCrl1SNSYmhmbNmhEcHMxbb72V6nzrvXv3ZtYs01l61qxZPPHEEzd8/S1btlC/fn38/f0BrptR9U6ohX6LipbJwzv05aVtoyAqCgID7S5JRMS12LR+6tq1a1m1ahWbN28mT548NGzYkEuXLuHl5ZUyaZmnpyeJiYkAPPPMMwwcOJA2bdqwdu1aRowYcd1zPvjgg0RHR7Nu3TqSkpKuWZAsNZZl3XCCtDuhFvot8vGBKfQjAS+S3p5gdzkiInKT/vjjD+666y7y5MnD3r172bJlS7rHlyxZEiBlDpXUPPbYY3Tt2vWa1nn+/Pk5d+7cdcfef//9rFu3jkOHDgHokrudHA64/5G7+ZgeJM/6CG5jqT0REXG+5s2bk5iYSNWqVRk+fDh16tS54fEjRoygY8eO1KtXjyJFiqR5XLdu3Th9+jRdu3ZN2fb444/z73//O6VT3F+KFi3KtGnTaN++PdWqVaNz5853fmJXaBz6bUhOhiqee9hDZXj9dRg+PFNeR0TEXbjTOPR/WrBgAREREXz88ccZ+rwah+4EHh6QEFCJZY6HYfJkuHTJ7pJERMQGzzzzDMOGDWN4FmjYKdBvU+PG8JY1GI4fhwz+X5mIiLiGd999l/3793PPPffYXYoC/Xb17QtracgPhJA87i1ISrK7JBERycYU6LepWjUYMMDBWIbisT8KFmk6WBGRG3HTLluZ4nb+rRTod2DCBPgmTzt+JZDTw8aCflhFRFKVK1cufv/9d4X6TbAsi99//51cuXLd0uM0scwdcDjg3fc8Gff4EKYfeJL4r1fj/XDGzcsrIuIu/Pz8iImJ4cSJE3aX4hJy5cqFn5/fLT1GgX6HwsLg0N4eHBnzCrwwhhIKdBGR63h5eaVMdyqZQ5fcM8DAF72ZwPOU2LMKtm2zuxwREcmGFOgZwMcHlpd+ijMUYHuXsfzyi90ViYhIdqNAzyDPveLDzNz9qL5/AWN6RdldjoiIZDMK9AzSqxcMPPQslx3ePLj5LbvLERGRbEaBnpGKF+f7oCcIYzYxP8TZXY2IiGQjCvQMlmPYYHKQSESjiSQk2F2NiIhkFwr0DHZf53J8kaMzPS5MZf6HZ+wuR0REsgkFegbLmRNarhuKD+f4ud9UzpyxuyIREckOFOiZIO8D1fjx7uY8x0QeaXYx/QeIiIjcIQV6Jqn22TCKc5xKW2dpuXQREcl0CvRM4mhQn5gyDzCUsTR4UL3jREQkcynQM4vDQcn3XqYMv1Fp+1xOn7a7IBERcWcK9EzkaNGc38vU4D+8yZxZSXaXIyIibkyBnpkcDjxffZlA9rNl0DySk+0uSERE3JUCPZMVDGtL7F2VeYlRzJyuRBcRkcyhQM9sHh4UGPcSVdjNtlcisCy7CxIREXekQHeCfE90IjpnIL2PjaR3LyW6iIhkPAW6M3h6ctfYF6nJdo7MWk6S+seJiEgGU6A7SYF+3TlbsDTDeYOP56iVLiIiGUuB7ixeXuQaMYwH2Ez0R2vtrkZERNyMAt2Jcj71BEcdvtRbP5IRI9CUsCIikmEU6M6UKxdHe7xAE1bzzWubaNXK7oJERMRdKNCdrPp7fUguXISXGMW338Lnn9tdkYiIuAMFurPlzYvHoIG0ZBk12EaXLmhsuoiI3DEFuh369YOCBXnv7jcAqFdPoS4iIndGgW4HHx8YOJDaRyO4l+1s3AgHDthdlIiIuDIFul2efRbuuouv67wGQGAgHDlic00iIuKyFOh2KVAABg6k+JYl9K2zDYC2bW2uSUREXJYC3U5XWulTir5G8eIQGQnx8XYXJSIirkiBbicfHxg0CL78kondIwF49FGbaxIREZekQLfbM89AoUJ03mfupS9cCP/7n801iYiIy1Gg2+1KK93x1Vcse+0HAMqWhd9/t7csERFxLQr0rOBKK7359yNSpoPt08fekkRExLUo0LOC/Plh8GAcy5ax6MWtAHz5pc01iYiIS1GgZxX9+0PhwuQYOYKxYyEhASZNsrsoERFxFQr0rOJKK52vv6Zj6e8BeO45+OUXe8sSERHXoEDPSq600v1nj2DuXLNpyxZ7SxIREdegQM9K8uWDF16A5ctp52uSvF8/LdwiIiLpU6BnNf36QZEi5BnzCr6+cPEirFtnd1EiIpLVKdCzmnz5YNgwWLmSjW+aJF+61OaaREQky1OgZ0V9+0KJEpSd/hJg8c03dhckIiJZnQI9K8qdG4YPx7FxIyMf+JqffoJz5+wuSkREsjIFelbVsyf4+9Pnt5dxkMz06XYXJCIiWZkCPavKmRNGjKBozA7as5Bly+wuSEREsjIFelbWrRsEBTHK4xUORyfZXY2IiGRhCvSszNMTXn+dCsm/UGv/J6xZY3dBIiKSVSnQs7r27Tlb/l5GMIJmjS8TH293QSIikhUp0LM6Dw983hlJOQ7Rk5m8/LLdBYmISFakQHcFLVqQVOdBhvMGk9++yJ9/2l2QiIhkNQp0V+Bw4DlmFCU5wtNMpW1buwsSEZGsRoHuKho0IKlxKC8ymi2rzvHjj3YXJCIiWYkC3YV4jh5JUU7yHBO59144dMjuikREJKtQoLuSWrXgkUd4OedbFOYk5cppSlgRETEU6K7mzTfxTrzAtDKjANTrXUREAAW66wkKgp49aXdkCmU5xDvvoLHpIiKiQHdJI0bg8PRkWc1XAOjUyeZ6RETEdgp0V1SyJDz3HBW3f0I1fmTJEjh1yu6iRETETgp0VzV0KI6CBfmy8osA9OgBa9faW5KIiNhHge6qChaEl16i1O7lNGI1y5ZBo0bw0Ud2FyYiInZQoLuyfv2gVClWVB/K2jUWAJMm2VyTiIjYQoHuynLlgjfewOvHSBqcWMD998OPP8Lvv9tdmIiIOJsC3dV17w5VqsB//kPHRxIAqFPH5ppERMTpFOiuztMTxoyB/ft5Lt90ataE/fs1g5yISHajQHcHDz8M9evjeP01+j9+HoCQEEhIsLkuERFxGgW6O3A4YOxYOHaM7sfH4+UFv/4KlSvbXZiIiDiLAt1d1KkD7duTY8JbXDx0lFKlICpKl95FRLILWwJ9woQJVK5cmSpVqtC1a1cuXbrEqVOnCA0NJTAwkNDQUE6fPp1y/OjRowkICKBChQqsWLHCjpJdw5gxcOkSnq+/yuuvm00dO9pbkoiIOIfTAz02NpZ33nmHyMhIdu3aRVJSEuHh4YwZM4YmTZoQFRVFkyZNGDNmDAB79uwhPDyc3bt3s3z5cvr27UtSUpKzy3YNgYFmbPr06TxWYxceHrBihYaxiYhkB7a00BMTE7l48SKJiYn8+eeflChRgoiICMLCwgAICwtj8eLFAERERNClSxe8vb3x9/cnICCArVu32lG2a3jlFfDxwWPoC4wbZzbNn29vSSIikvmcHuglS5Zk8ODBlC5dGl9fXwoUKEDTpk05duwYvr6+APj6+nL8+HHAtOhLlSqV8ng/Pz9iY2NTfe5p06YREhJCSEgIJ06cyPyTyYoKFYLhw2H5cp6r9A1gvhUREffm9EA/ffo0ERERHDp0iCNHjnDhwgXmzp2b5vGWZV23zeFwpHpsnz59iIyMJDIykqJFi2ZYzS6nXz8oVw7PoYOpHpzEyZNw5YKHiIi4KacH+qpVq/D396do0aJ4eXnRvn17Nm3aRPHixYmLiwMgLi6OYsWKAaZFfvjw4ZTHx8TEUKJECWeX7Vq8vc0wtp9/ZlnHWQC0awcXL9pcl4iIZBqnB3rp0qXZsmULf/75J5Zl8e233xIUFESbNm2YPXs2ALNnz6Zt27YAtGnThvDwcOLj4zl06BBRUVHUqlXL2WW7ng4d4MEH8Z3yMk93N2PXRo+2uSYREck0Tg/02rVr869//YsaNWoQHBxMcnIyffr0YdiwYaxcuZLAwEBWrlzJsGHDAKhcuTKdOnWiUqVKNG/enClTpuDp6enssl2PwwH//S8cO8bk0qZ33Btv2FyTiIhkGoeV2k1qNxASEkJkZKTdZdiva1eIiKBXvV+Z+Y0f1arBt99C4cJ2FyYiIrfqRtmmmeLc3ejRkJzM1MIv4eMDO3fCl1/aXZSIiGQ0Bbq7K1sWBgwg52dzOPLVdsAMYzt40N6yREQkYynQs4P//AeKFCHvK4OoWMEiJsZsEhER96FAzw4KFIDXXoO1a/np9UU88IAuu4uIuBsFenbRpw9UqYLX0EH4FbnEn3/Clcn4RETEDSjQs4scOWDiRIiO5uU84wHTaBcREfegQM9OmjSBdu2o8uWblCCW996DKVPsLkpERDKCAj27efttHImJRD5kJu7p3x+mTrW5JhERuWMK9OymXDkYNAjfVXMJH7AZgL59YfJkm+sSEZE7okDPjl58EXx96bxpAEePJAPwzDNw6pTNdYmIyG1ToGdH+fKZ1dh++IHi33zMwIFm89NP21uWiIjcPs3lnl0lJ8MDD8D//ge//orDJz+enpCYaHdhIiKSFs3lLtfz8IBJk+DoURg1iv79ISkJrlp6XkREXIgCPTurXRvCwmDCBNpU2g/ARx/ZW5KIiNweBXp2N3o05MxJw68GA/DKK3DmjL0liYjIrVOgZ3e+vvDSS3gti2B6xxUALFlic00iInLLFOgCzz8PgYGEbXuGnMQzd67dBYmIyK1SoAt4e8O775LjYBRjivyXlSvBPcc+iIi4LwW6GM2aQfv29D09klL8Rt++dhckIiK3QoEuf5swgZw5YQLP8/77sHq13QWJiMjNUqDL30qXxjF8OB1YSFNW0KQJLFhgd1EiInIzFOhyrYEDITCQRSVNB7nu3e0uSEREboYCXa51pYNcntgopt3zX+LjYfFiu4sSEZH0KNDles2aQYcO9PjNdJBr1w4uX7a7KBERuREFuqRu/Hg8PBws8HsegOnTba5HRERuSIEuqStdGl5+mVoxpoPcmDF2FyQiIjeiQJe0DRwI99zDNO9nOHY4nqQkuwsSEZG0KNAlbVc6yJWJj2II47QSm4hIFqZAlxtr2pTEDp14iVGM6R3Fnj12FyQiIqlRoEu6crw7EXJ68x59afeIJnkXEcmKFOiSPl9fck0YTSirqBn1me6li4hkQQp0uTlPPUVc6VpM4Hk+GHPa7mpEROQfFOhyczw9KbLgAwrzO3lHDrO7GhER+QcFutw0r/uq83Gh5wi7NI2zyzfZXY6IiFxFgS63pPC7I/iNUpx79ClISLC7HBERuUKBLrekddd89GcyJU/vwvrveLvLERGRKxTockscDvDt04ZFPELC8Nfg0CG7SxIRERTochsmTIBneYfLSZ7Qrx9YGpsuImI3Bbrcsjx54P6OpXjZegO+/hrmz7e7JBGRbE+BLrfl2WdhMv3ZRg2sZ5+FU6fsLklEJFtToMttqVsX2nfMQS9mwMmTMHiw3SWJiGRrCnS5bf36wU6qs6v5CzBrFqxaZXdJIiLZlgJdbltAgPlca+krnC95Dzz5JFy4YG9RIiLZlAJdblvJkvDZZ3CJ3Dwc+yFER8Pw4XaXJSKSLSnQ5Y506QKDBsEG6jOVf2NNmgRbt9pdlohItqNAlzv29tvQqhUMZSxHHb5YPXvB5ct2lyUikq0o0CVDLFwIRfx96JM0FcfuXTBmjN0liYhkKwp0yRBeXrBrF3xFaz6jC9bIkbBnj91liYhkGwp0yTB58sDo0TCASSTkyg+9e0NSkt1liYhkCwp0yVDt2sEJivHvixNh82Z47z27SxIRyRZuGOibN2+mX79+VK1alaJFi1K6dGkefvhhpkyZwh9//OGsGsWFVKgAjRrBrMTunKvbHF58EQ4etLssERG3l2agt2jRgunTp9OsWTOWL19OXFwce/bsYeTIkVy6dIm2bduyZMkSZ9YqLuL11wEcLGo+DTw9oWdPSE62uywREbfmsKzU1748efIkRYoUueGDb+YYu4SEhBAZGWl3GdnS5cvg7Q1BQbBn0AxzL/3dd6F/f7tLExFxaTfKtjRb6KkF9VdffZXuMSI5c0KRIvDLL7CpYk9o3hyGDoUDB+wuTUTEbd1Sp7hXXnkls+oQN/Pdd+bzt6sd8OGHZlybLr2LiGSaWwr0NK7Oi1znnnvM56++Avz8YOJEWL8eJk+2sywREbd1S4H+wQcfZFYd4mYcDvD3N9O6JyUBYWHQsiUMGwZRUXaXJyLidtIM9O/+umZ6lVq1aqV8ffbsWXbt2pU5VYlbePRR83nsWEzCf/CBucH+xBOacEZEJIOlGehffPEFDzzwAK+//jpLly5l69atrF+/npkzZ9KjRw9atWrFxYsXnVmruJhhw8znjz++sqFkSXjnHdi40XwWEZEMk+awNYDTp0+zYMECNm7cSFxcHLlz5yYoKIiWLVtSt25dZ9Z5yzRsLWto1AjWroUFC6BDB8CyoE0bWLUKdu78+2a7iIik60bZdsNAd2UK9Kxh5Upo2tR8/dJLMHIkcOQIVKkCFSvChg1m8hkREUnXbY1DF8kIoaHw/ffm61Gj4M8/gRIlzCX3zZvNYuoiInLHFOiS6WrVutIxDtNiB6BbN3MNfvhwc+ldRETuiAJdnOKvHu/z5l3Z4HDA++9D4cLQvTtcumRbbSIi7iDdQJ8/fz7nzp0DYOTIkbRv357t27dnemHiXvz8oGhR2Lbtqo1FisDMmbBrF7z8sm21iYi4g3QD/Y033iB//vx89913rFixgrCwMJ5++mln1CZuJjgY9u0za7WkaNECnn4axo833eFFROS2pBvonld6IC9dupSnn36atm3bcvny5UwvTNzPjBl/f1669Kodb70FAQFmNrk//rClNhERV5duoJcsWZKnnnqKefPm8fDDDxMfH0+yFtiQ21C2LPz0k/l65MirduTNa2afiY2FZ5+1ozQREZeXbqDPmzePZs2asXz5cgoWLMipU6d46623nFGbuKHgYBPsW7aYcekpatc2G+bMMbPQiIjILUlzYplTp07d8IGFChXKlIIyiiaWybp27TLBDmbEWtWqV3YkJMADD8DBg+YgX1/bahQRyYpulG050npQzZo1cTgcWJaV8vkvDoeDgwcPZnylki1UqWJmfn3oIVi48KpA9/KCuXPh3nuhVy9zo93hsLVWERFXkWagHzp0yJl1SDZTv775PGcOjBhx1Y4KFUwnuf79zTh1jagQEbkpaQb61U6fPk1UVBSXrpr8o/5ff5FFboOXF9SoAdu3w+nTcNddV+3s2xe++goGDoQGDaBSJdvqFBFxFel2ips+fTr169enWbNmvPrqqzRr1owR1zSpRG7PoEHm89Ch/9jhcMCsWeDjA126aBY5EZGbkG6gT5o0iR9++IEyZcqwZs0aduzYQdGiRZ1Rm7i5rl3N5w8/hJMn/7Hz7rtNqP/8MwwZ4vTaRERcTbqBnitXLnLlygVAfHw8FStWZN++fZlemLg/h8OswAZ/93q/xsMPw4AB8O67/5iJRkRE/indQPfz8+PMmTM88sgjhIaG0rZtW0qUKOGM2iQbePFF06n96FHYvz+VA8aOhWrV4PHHIS7O2eWJiLiMNMehp2bdunX88ccfNG/enJw5c2ZmXXdM49Bdx6pVZt30fv1g8uRUDvjlF6hZEx58EFasAA8tEigi2dONsi3dv4y//fZbyoe/vz/Vq1fn6NGjGV6kZF8PPWRWYpsy5arlVa8WFAQTJ5rkHz/e2eWJiLiEdIettWzZMmVimUuXLnHo0CEqVKjA7t27nVGfZBMdOphh5507w7lzZl6Zazz5pGmd/+c/0KiRabGLiEiKdFvoP//8Mz/99BM///wzUVFRbN26lbp16zqjNslGpk6F6Gjz9dtvp3KAw2G6wxcvbrrHnz/vzPJERLK8W74ZWaNGDX744Yc7etEzZ87wr3/9i4oVKxIUFMTmzZs5deoUoaGhBAYGEhoayunTp1OOHz16NAEBAVSoUIEVK1bc0WtL1lWmjJlHZu9eeOONVA4oVMhMDbt/v1ZlExH5h3QvuY+/6p5lcnIy27dvv+Nx6AMGDKB58+YsWLCAy5cv8+eff/Lmm2/SpEkThg0bxpgxYxgzZgxjx45lz549hIeHs3v3bo4cOcJDDz3Er7/+mrJOu7iXsWOhTh145RUoVcp0br9GgwbmsvuoUdCkCXTrZkeZIiJZTrot9HPnzqV8xMfH07JlSyIiIm77Bc+ePcv69evpdeUmac6cOSlYsCARERGEhYUBEBYWxuLFiwGIiIigS5cueHt74+/vT0BAAFu3br3t15esrXZt+KsD5xNPwJ9/pnLQiBFQty489RRoTgQREeAmWuivvvpqhr7gwYMHKVq0KE888QQ7d+6kZs2aTJo0iWPHjuF7ZblMX19fjh8/DkBsbCx16tRJebyfnx+xsbGpPve0adOYNm0aACdOnMjQusV5atY0a7NMngyvvmrWarlGjhzw2WdmAHvHjvD995A7ty21iohkFWkGeuvWrXHcYOnKJUuW3NYLJiYmsn37dt59911q167NgAEDGDNmTJrHpzZMPq26+vTpQ58+fQAzVk9c19ixJtDfftv0fL/u7fTzg48/hhYt4Lnn4IMP7ChTRCTLSPOS++DBgxk0aBD+/v7kzp2bJ598kieffJJ8+fJRpUqV235BPz8//Pz8qF27NgD/+te/2L59O8WLFyfuykxgcXFxFCtWLOX4w4cPpzw+JiZGM9VlA3nywKJF5uv77oOYmFQOat4chg2DadNMi11EJBtLM9AbNGhAgwYN2LFjB59//jmtW7emdevWfPrpp3z33Xe3/YJ33303pUqVSpkP/ttvv6VSpUq0adOG2bNnAzB79mzatm0LQJs2bQgPDyc+Pp5Dhw4RFRVFrVq1bvv1xXW0bWuGnANUrZrGQW+8YWaQ69MHfv3VabWJiGQ16d5DP3HiBAcPHqRcuXIAHDp06I7vT7/77rt069aNy5cvU65cOWbNmkVycjKdOnVixowZlC5dmvnz5wNQuXJlOnXqRKVKlciRIwdTpkxRD/dswuGAb781t8p37jTrs7Rs+Y+DcuSA8HCoXh06dYItW+DKYkIiItlJunO5L1++nD59+qQEenR0NB988AHNmjVzSoG3S3O5u4/jx818Mu3awcKFaRy0bJlJ+3//28xSIyLihm6Ubem20Js3b05UVBR79+4FoGLFinh7e2dshSI3UKyY+fjxxxsc9PDDZt30ceOgYUPTk05EJBtJM9BXr15N48aNWfiPJtGBAwcAaN++feZWJnKVmjXh66/NVfVUF3ABGDkSvvvOzPteowYEBjq1RhERO6UZ6OvWraNx48Z8+eWX1+1zOBwKdHGq9983U8POnw8ffZTKDHIAXl5/30//179g82bTXV5EJBtI9x56UlKSS3ZC0z1093PqFBQuDGXLmvne07zzs2KFGZ/evTvMnm1614mIuIE7Wg/d39+fPn368O2336Y6yYuIsxQqZHq8R0dDmzY3OLBZMzM97Mcfm6a9iEg2kG6g79u3j4ceeogpU6bg7+9P//7972gcusidWLnSjFT75hu4cOEGB778sukoN2AAaO5/EckG0g303Llz06lTJxYuXMiOHTs4e/YsDRo0cEZtItcpXPjvud3z5YMZM9I40MPDtNBLljT300+edFqNIiJ2uKn10NetW0ffvn2pUaMGly5dYl6a3YxFMt+AAWb6doDeveG11+C33+Ds2X8cWKgQfPGFGcjetSskJTm7VBERp7mpe+gTJ06kXr167Nq1i3nz5tGhQwdn1CaSKocDJkwwWQ3mdnmZMlCgAPz00z8OrlEDpkyBVavM0m0iIm4q3Ylldu7ciY+PjzNqEbkl7dub9dIXLIAffoB33zXb9u//x4G9epkhbKNGmQXXW7e2pV4RkcyU5rC1cePGMWTIEJ599tlUH/jOO+9kamF3SsPWsp8yZcyl96+/NguxXePSJbOIy4EDsG0blC9vS40iInfitqZ+DQoKAqBmzZqZU5VIBpszx8z62qKFWczlmhXacuUyTfmaNaFDB9i0SZPOiIhbSXdiGVelFnr2NH26mfkVzCywDz74jwO+/tos4tK5M3z6qSadERGXclst9NatW+O4wR+7JUuW3HllIhmsd284eBBGj4a6dVO5/N6ihbmX/p//mFlqhgyxrVYRkYyUZqAPHjwYgIULF3L06FG6d+8OwGeffUbZsmWdUpzI7XjzTahSBbp1M/n98cdQubLJbwCGDTNLtw0bZq7LX3fDXUTE9aR7yb1+/fqsX78+3W1ZjS65S/v2sGjR39/HxkKJEle+uXABHnjA9KL74QcICLClRhGRW3FHc7mfOHGCgwcPpnx/6NAhTpw4kXHViWSSefMgKsrMAgtm0rivvoKYGEjImRcWLwZPT2jbFs6ds7VWEZE7le449AkTJtCwYUPKlSsHQHR0NB988EGmFyZyp3LkMA3vN96AYsXg2Wf/HoJ+770QGemPx+efm8VcHnvMzFTjcVOTJ4qIZDk31cs9Pj6evXv3AlCxYkW801y3MuvQJXf5p7VrzaQzQ4eapVgrVDC30nO9PxGef95MOafZ5EQkC7utXu5X27ZtG9HR0SQmJrJz504AHnvssYyrUMQJGjY0H+3bm0Ve9u0zHeaeHDAAduwwgV69urkELyLiYtIN9B49enDgwAGqV6+Op6cnAA6HQ4EuLqtQIYiPB29v6NMHgoMd1Hn/fdizB7p3h++/h0qV7C5TROSWpBvokZGR7Nmz54Zj0kVcTc6cZsK4L76A+++HpKTceCxaBCEh5kb71q2mGS8i4iLS7QFUpUoVjh496oxaRJxqwQL497/N15Urw0+n/Mw4t9hYc13+8mV7CxQRuQXpttBPnjxJpUqVqFWr1jWd4TRTnLiD4cNh+3bTIH/lFVi8+H6YMcNcen/6aTOXrK5OiYgLSDfQR4wY4YQyROxRooS5ZV61KkREmE5yPXp0g717YeRIcy990CC7yxQRSZcWZxEBVq2C0FDz9cWLkCtnslnA5YsvYMkSaNXK3gJFRLjNmeLq1q0LQP78+fHx8Un5+Ot7EXfy0EN/309v1gwzwczs2VCjBnTtCj/9ZGt9IiLpUQtd5IrkZPDzg7i4K630XJgOcrVqmWnntm6F4sXtLlNEsrE7mstdJLvw8Pi7lT5z5pWNJUuaS+4nTkC7dnDpkm31iYjciAJd5CovvWQ+DxgA/ftf2VizJsyZA5s3mwXX3fOiloi4OAW6yFU8PeHTTyExEaZMuSq7//Uv0+v9k0/g9ddtrVFEJDUKdJF/6NrVjE8HeOQRE+4A/Oc/8PjjZs732bPtKU5EJA3pBvo/e7n7+PhQqlQp2rVrd8066SLuZMgQ83nJEnjnnSsbHQ744ANo0sRcev/2W9vqExH5p3R7ub/66quUKFGCRx99FMuyCA8P5+jRo1SoUIGpU6eydu1aJ5V6a9TLXe7U0aPg62tyPCnpqgnj/vgD6taF336DjRuhShVb6xSR7OOOerkvX76cp556KqWl3qdPH5YtW0bnzp05ffp0hhcrklXcfbe5dW5ZUKrUVffTCxSAZcsgb154+GE4csTWOkVE4CYC3cPDg3nz5pGcnExycjLz5s1L2acV2MTdffABVKxohqNXrQqjR18J9lKlYOlSOH3azCJ3/rzdpYpINpduoH/yySd8/PHHFCtWjOLFi/Pxxx8zd+5cLl68yOTJk51Ro4htChWCdesgKAh27TL94saPv7Lz3nth3jwzi1znzlf1nhMRcT7NFCdyk/66p/7AA+bWeYpp0+Cpp8zH1KlanU1EMs2Nsi3d1dZOnDjBhx9+SHR0NIlXtUBmpkylJZI93H23uey+aZOZMC5Xris7+vSBQ4dgzBjw94ehQ22tU0Syp3QDvW3bttSrV4+HHnoIT09PZ9QkkmV16mSusL/4IkyYcNWOUaMgOhqGDTNrsvboYVeJIpJNpXvJvXr16vz4449OKifj6JK7ZAbLMnO+g+ngPncu3HXXlZ3x8Wbj+vVmAHuLFrbVKSLu6Y6GrbVq1Yply5ZleFEirsjhMLfMS5Y0I9cKFYKVK6/s9PaGRYsgONiMd/v+e1trFZHsJd0Wev78+blw4QLe3t54eXlhWRYOh4OzZ886q8bboha6ZLann4b33zdf79oFlStf2XHsmOk598cfpvdchQq21Sgi7uWOWujnzp0jOTmZixcvcvbsWc6dO5flw1zEGaZONTPAAgwefNWO4sVhxQqz0kvTpmYQu4hIJkuzU9zevXupWLEi27dvT3V/jRo1Mq0oEVfx4Ydw4AAsXw6LF5vFXAAICICvv4YGDaB5c9iwAQoWtK9QEXF7aQb6+PHjmTZtGoMGDbpun8PhYPXq1ZlamIir6N8f1qyBdu3MVXYfnys7atQw99QffhjatDGt9ty5ba1VRNyXJpYRyQAvvmiGoYOZBTZv3qt2fv65WZO1bVtYsMBcihcRuQ13dA+9WrVqjB49mgMHDmR4YSLu4s03zaqqAB06/GNn584waZK5Jv/vf1+1youISMZJN9CXLFmCp6cnnTp14r777uPtt9/mt99+c0ZtIi7D4YBVq8zXK1bAk0/+44BnnoGXXoLp081i6wp1Eclg6QZ6mTJlGDJkCNu2bePTTz/lp59+wt/f3xm1ibicmBjzefp0iIv7x8433jA33N9+28wsJyKSgdKd+hUgOjqaefPm8fnnn+Pp6cm4ceMyuy4Rl1SypAnz3r2hbl3TAz6Fw2EuvZ89C8OHm95zzz5rW60i4l7SDfTatWuTkJBAx44dmT9/PuXKlXNGXSIu64kn4JVX4OBBMwx9woSrJp3x8IAZM+DcORgwwIT644/bWa6IuIl0A3327NlUrFjRGbWIuAUPDzOde4UKZlrYKlWgWTMYOdJsy58/B3z2GbRuDb16Qf78qfSkExG5NWkG+ty5c+nevTvLli1LdS73gQMHZmphIq6sfHm4fNnMLdO+vekot2KF2demDRQq5M3kuYvI266pGdL25Zcm9UVEblOaneIuXLgAmKlf//lx/vx5pxUo4qo8PKBlS7MI26pVMGgQ3HOPWYjto48gX/G8zO60FCpVMrPSfPed3SWLiAu7rYllJk6cyHPPPZcJ5WQcTSwjWVVCgukL9/77ZjbY0/uOQ716cPSomXJO0yqLSBruaGKZ1IwfP/6OChLJzry8zMIuXbvCmTPw4oRipgl/110QGgo//WR3iSLigm4r0N10tlgRpxo71nweMwYmflEKVq+GPHnMlHO7d9tbnIi4nNsKdIfDkdF1iGQ7pUr93Rh//nn4Jb6cCXUvLxPqe/faW6CIuJQ0Az1//vz4+Phc95E/f36OHDnizBpF3FZwsFmCFUzfuFnfBXJqwZWVDBs3hqgo+4oTEZeSZqCfO3eOs2fPXvdx7tw5EhMTnVmjiFvr2RP+6mPasycUfrAiK4Z8i5WYaEL94EFb6xMR13Bbl9xFJON4eJjZ5A4fNtPFAjQfVJlaf6wi+cKf0KgRREfbWqOIZH0KdJEsws/PzDD3V7BHXq5KyOmVXDpxluRGjc0OEZE0KNBFshCHwwT76tVmWNsOalDv4jfEx/5uLr/HxtpdoohkUQp0kSzIyws+/dRcad+V6z4aJazg4v+OQYMGaqmLSKoU6CJZWJkyMH8+fE8dGiV8w8XDJzhXowHJB6PtLk1EshgFukgW16oVrF0LkZ51qH95FYknT3MksAGX96r3u4j8TYEu4gIaNIDERFh+8j663/0tuZPPc7JyAzbM2k9yst3ViUhWoEAXcSGFC8Mnv9SgT/nV5Ey+SLmeDfjizX12lyUiWYACXcTFFCwIX+yvxrHwNXiRQN3hDZnw5B7efRcuXrS7OhGxiwJdxEVV7hzM5/9eiwOLbtMb8sGzu7j/ftDaSSLZkwJdxIU9M7UShX5aRxFfL9bSEHb+SP/+dlclInZQoIu4uJzBFfBYv467SuRhLQ3Z8d4m1q+3uyoRcTYFuog7CAjAc9MGPIoXZSWhvN5gFYGB0K+f3YWJiLMo0EXcRZky+Py4gXi/8nzt0ZJK+yN47z2z1vqhQ3YXJyKZTYEu4k7uvptCO9fiFVKdRR4deJRPmDgRypXT0uoi7k6BLuJuChWCVavwqF+PTxw9WNHufQA6djSLvhw9anN9IpIpFOgi7ih/fli2DB5+mKaLnua1vOPYuROaNIHy5SEhwe4CRSSjKdBF3FXu3LBoEXTpwisXhvK/7i9Rp7bFn3/CqFGQnKwx6yLuRIEu4s68vGDuXOjdm9Jz32Rj9X54kMRrr4GnJ3h4wIYNdhcpIhkhh90FiEgm8/SEadOgUCE8xo0jrt5xPmwwl/8dy8WHH0L9+rB8OTRrZnehInInbGuhJyUlce+999KqVSsATp06RWhoKIGBgYSGhnL69OmUY0ePHk1AQAAVKlRgxYoVdpUs4rocDhg7FsaPp9iGL3hpQ3OmvfUHgwaZ3c2bQ40asHGjvWWKyO2zLdAnTZpEUFBQyvdjxoyhSZMmREVF0aRJE8aMGQPAnj17CA8PZ/fu3Sxfvpy+ffuSlJRkV9kiru35580l+I0boUED3h4Ux8qVpg/djh1Qty58843dRYrI7bAl0GNiYli6dCm9e/dO2RYREUFYWBgAYWFhLF68OGV7ly5d8Pb2xt/fn4CAALZu3WpH2SLuoVs3WLoU9u+HBx7goTJRnDljNoG59L5nj60VishtsCXQn3vuOcaNG4eHx98vf+zYMXx9fQHw9fXl+PHjAMTGxlKqVKmU4/z8/IiNjU31eadNm0ZISAghISGcOHEiE89AxMU1bQpr1sD58/Dgg3hsj+Thh2H0aLP7vffsLU9Ebp3TA/2rr76iWLFi1KxZ86aOt1IZV+NwOFI9tk+fPkRGRhIZGUnRokXvqE4Rt3fffebSe9680LAhfPMNQ4eanu9TpsDZs3YXKCK3wumBvnHjRpYsWULZsmXp0qULq1evpnv37hQvXpy4uDgA4uLiKFasGGBa5IcPH055fExMDCVKlHB22SLu6Z57TKiXLw8tW+L4ZC4DBphdaqWLuBanB/ro0aOJiYkhOjqa8PBwGjduzNy5c2nTpg2zZ88GYPbs2bRt2xaANm3aEB4eTnx8PIcOHSIqKopatWo5u2wR91WiBKxbZ3rE9ejBOJ+RgMWSJWZeGhFxDVlmHPqwYcPo1KkTM2bMoHTp0syfPx+AypUr06lTJypVqkSOHDmYMmUKnp6eNlcr4mYKFjSD0Xv1Isdrw1lZ+hAtNr9P+/Ze/PwzVKlid4Eikh6HldpNajcQEhJCZGSk3WWIuBbLgldegZEjORUSStnIBTz6lA/vv293YSICN842Tf0qIn9zOOCNN2DGDAr9uIbNHnVZ+sFhLlywuzARSY8CXUSu17MnLF1KOc9otlCHB/P9SHi4FnMRycoU6CKSuqZNyRW5EU8vDzZQj9ldv2bIELuLEpG0KNBFJE2OqsEUP7gFzwoBfElrzr79AQMH2l2ViKRGgS4iN+TwK0meH9ZzuWFTPuDflJ7wHF07JhIdbXdlInI1BbqIpC9/fvKsXMLJbgN4jkmELWhF7Yp/2F2ViFxFgS4iNydHDorMnUj85Gk85PiWNfH3c3zzAburEpErFOgicku8+z3J9jErKc4xvOrWgrVr7S5JRFCgi8htqP5cQxrn+Z6jycVIfigUPvzQ7pJEsj0Fuojcspw5YdaGAOqwhW+SmkCfPpx6/HlISrK7NJFsS4EuIrelRg2Y9UUBWvEVExlAodkT2VayNSeizthdmki2pEAXkdvWvj2cPJ2DAjMn0ocPqHpsJafvqUWrcnvYvBni4+2uUCT7UKCLyB0pWBCeeAKmXO7D1y+s4S6Ps3x2qDZvPbCQXLkgIgL27oXLl+2uVMS9KdBFJEN4eUGbcXUp+ts2HJUrs5AOvMHLtH8kiaAgKFfO7gpF3JsCXUQyVsmS5Nu2Dnr14mVGcbh6G/zvOkNsLMybZ3dxIu5LgS4iGc/b2wxle+89Suz6hn0F7qMSu+ncGcaOtbs4EfekQBeRzOFwwNNPw5o1eF08x4+56tCOhQwbBuPHa4SbSEZToItI5qpbF7Ztw6uaua/+Ji8yZFAinTrZXZiIe1Ggi0jmK1kS1q2DJ5/kRcawyhHKxoVHiYy0uzAR96FAFxHn8PaGadPgo4+o6/U926nBq43W212ViNtQoIuIc4WFkSPye5Jy5yPifGN2PPoWiQmW3VWJuDwFuog4X3AwUZ9Gsoh23PvZEL7K2Y5DO87YXZWIS1Ogi4gtGj/iQ53/zWNW9Ym0ZCke99Xkwnc77C5LxGUp0EXENqVKO3hixwB6B67HI+kynvXuJ2rIh2DpErzIrVKgi4jtJmy5n3alt7Oe+gS+1Yf53t2Y+PpZVq+2uzIR16FAFxHbFSoEm/cXJfear3mvxEjaJcyj1as1GNIkkg0bIDHR7gpFsj4FuohkCV5eUK+hJ31jXyLp23WUKHKZTTzAovrjyemVzOnTdlcokrUp0EUky/Fu/CB59v3IqTotGc8gvqQ13ZudsLsskSxNgS4iWVOhQty9aSHWu5Npwrd8+EM1nqu2hmnTIC7O7uJEsh4FuohkXQ4Hjv79iPr4e857+DD+pyYcfeoVSpVIZMQI2L/f7gJFsg4FuohkecHdqxFwZhuXOoXxCm+wlobMfu0QgYHo3rrIFQp0EXEJHvnzkufzWfDJJ9T1+Zk9XtV4jNkUKmTRogXExNhdoYi9FOgi4loefRR++gnvOvcym8dZ4t2Jrct/p1QpWLxY66xL9qVAFxHXU6YMHmtWw9ixtE6O4NdcVXmIlbRrB+3bw6lTdhco4nwKdBFxTZ6eMGQIfP89hfwLsJKmjOd5Viy5xEMP2V2ciPMp0EXEtd17L45t2+CZZ3ieiez0CiFpx05depdsR4EuIq4vd2545x1YvpyS3r+zlVoM9x7H6JFJXL5sd3EizqFAFxH30awZ3r/+TOTdrXkzaSgNh9cl2HsfW7faXZhI5lOgi4hb8fItwoNH5nPmvU+pnGMfP1Kdz2pPIOmyrsGLe1Ogi4j7cTgo+HRXfH7bza+lQ5nAQHYVbci5HZpaTtyXAl1E3JevL8W3RNDLaw6lz+7Cs0ZV9vV/F5KT7a5MJMMp0EXErd3t6+DDSz14ue0u1tKQClOe5XhwYzh40O7SRDKUAl1E3J6HB0xZXJKC3y2lJzPw3rOD5OCqMHmyWuviNhToIpJtPPCggzz9elKFXay8WBeeeQbq1oU9e+wuTeSOKdBFJFt55x0oVqMUza2v6cEckvfug3vvhddfR4PWxZUp0EUkW/HwgG3b4IUXHMylBx0r/wIdOsCrr0KNGrBli90litwWBbqIZEvjxkGePLDwu2L4rf+UBWFfknT6D3jgARgwAM6ft7tEkVuiQBeRbGvzZmjQAGJjoePsVhQ8soc/uvWFd9+FypVh+XK7SxS5aQp0Ecm2qlaFtWvh9GkT7OfJT8G5kxnbagNxZ/NAixZm/fWjR+0uVSRdCnQRyfYKFjTB3q8f+PrCpMgHKXvmR0bwKknzv4AKFcwQNy3hJlmYAl1E5IrJk+HIEfPx7QZvXmMEQYk/E12slhniVrs2REbaXaZIqhToIiKpqFsXZs+GKO7Bf/83TGv0mbnZXqsW9O8PZ87YXaLINRToIiJpeOwxOHYMwMFTa7qwa8FeE+ZTp0LFivDpp2BZdpcpAijQRURuqFgxWLzYfN24XQEzM83WrVC6NHTrBqGhsG+frTWKgAJdRCRdbdvCPffAiRNXOrzXrGnGvE2ZYu6pBwfDkCFw9qzdpUo2pkAXEbkJb71lPvv6mk5zeHpC376mdd6jhzmgQgWYM0cLvogtFOgiIjehVSuoU8d8XbIkzJt3ZUfx4jBjBnz/vbkMHxYGDz6o3vDidAp0EZGb4OFhrrLPmWO+f/zxf/SHq1XLHDBrFhw6ZL7v3RuOH7ejXMmGFOgiIregRw/TUr940QxN/9//rtrp4WGSft8+GDjQjHu75x6YOBESEmyqWLILBbqIyC365hvzecoUKFsWdu36xwEFCsDbb8PPP5vJaJ5/HqpVg2XLNMxNMo0CXUTkFuXPD7/+Cq+9Zr4PDjad5Xx9oXNnWL/+yoEVK5oFXhYvNi30li2haVPYudOu0sWNKdBFRG5DYCC88gpMmAB9+phOc0ePms5yDRrAzJlXDnQ4zLi33bvNwdu2wb33Qq9eV7rLi2QMh2W55/WfkJAQItXLVEScKDHRBPlTT5nvmzSBokWhTRvo2vXKQadOwciRZuJ4Ly8zfn3wYMib17a6xXXcKNvUQhcRySA5cpjW+u7dpgV/+DCEh5sVWKdNu3JQoUIwfjz88otZnnXECNNxbtYsreYmd0SBLiKSwSpVMvfY9+2DJUvMtqeegqAg+OAD+OEHoHx5WLAANmwAPz/o2RNCQmDVKltrF9elQBcRyUStW8PXX4O/P+zdC//+txmiPmQIjB4Nl2vVNePXP/sMTp82c8OHhmpiGrllCnQRkUzWvDkcPGjmmJk61Wx76y34z3/A2xt+2uUBXbqYJv2ECfDjj3DffdCpkxZ+kZumQBcRcZKiRU0L3bLMdO9Nmpjt1aqZ2+mRP3vDc8/BgQPw6qumaV+5srkxHxtra+2S9SnQRURs4HCY2+Xvvw933WWGq993H7RvD+HLfDg9YIQJ9n794KOPICAAhg41veRFUqFAFxGx0VNPwe+/m0nkABYtMkPcChWC8XOLcXDAJHPZvWNHc52+fHkYMwYuXLC3cMlyFOgiIjZzOMwl9/PnTXa3bGm2DxoEL78M54v6m1Vhdu6EunXhxRehXDkz/O3PP+0tXrIMBbqISBaRN68Zkv7VVybcH3zQdH7Pn9+05A/mDYYvv4SNG6FqVZP45cvDO+/ApUt2ly82U6CLiGRBefPChx+adV3ATExTvrzZxgMPwMqVsG4dVKgAAwaYnVOmQHy8rXWLfRToIiJZVFCQuap+4QJMn2629eljZqGrWBEWn6oPa9fC6tXmEnz//mbnBx/A5cu21i7Op0AXEcni8uQxa7msWmU6zN13n7nX3q4dPPkkfH68EfEr15t1XUuWNGPj7rkHZszQOuzZiAJdRMRFNGkCn35qPv5qsU+fbuakyZXbwROfhmJt3GS6zBctCr17m+Fu772ne+zZgAJdRMQF9eplJqg5csQMTwczXL1YcQet32vBlLCtJCz6yrTY+/Uzc8/+97+mt524JQW6iIgL8/U1w9IvXTJ95UqVMr3k+z/jwO+plsSv3mjusVeqZJZpLVvWLN965ozdpUsGU6CLiLgBb28zmm37djh0CEqXNnPHlyvvYOLORhz79FvYtAnq1IHhw6FMGXjpJThxwu7SJYMo0EVE3EzZsmZN9lKlzCX555+H4GD4o9L9pvm+fTs0bWqWeytbFgYOhJgYu8uWO6RAFxFxQ/nywf/+B2fPmhFtJ05AwYIwcSJsS76XszPmm9Tv0MFMTOPvD48/Drt22Vy53C6nB/rhw4dp1KgRQUFBVK5cmUmTJgFw6tQpQkNDCQwMJDQ0lNOnT6c8ZvTo0QQEBFChQgVWrFjh7JJFRFySw2FmmfvlF3OfHUxrPSQEatfGDHSfMwf274e+fWH+fNOUf/hhWLPG9LoTl+H0QM+RIwf//e9/+eWXX9iyZQtTpkxhz549jBkzhiZNmhAVFUWTJk0Yc+Wnb8+ePYSHh7N7926WL19O3759SUpKcnbZIiIuK2dO0xN+925YsgRKlIC9e2HUqCsHlC0LkybBb7/BG2/Atm3QuDHUqgXz5kFiop3ly01yeqD7+vpSo0YNAPLnz09QUBCxsbFEREQQFhYGQFhYGIsXLwYgIiKCLl264O3tjb+/PwEBAWzdutXZZYuIuLxKlaB1azNBDZiFXwoXhhEjzMfrUwrz22MvQ3S0Wdf1jz+gc2czSc2UKVoIJouz9R56dHQ0O3bsoHbt2hw7dgxfX1/AhP7x48cBiI2NpVSpUimP8fPzIzY2NtXnmzZtGiEhIYSEhHBCPTdFRFIVFGT6wFWpYpZXf+018/Hqq6bz+469uc1qML/8AgsXQvHiZlrZ0qVND/m4OLtPQVJhW6CfP3+eDh06MHHiRHx8fNI8zkrlHo7D4Uj12D59+hAZGUlkZCRFixbNsFpFRNxNyZLw88/mNvlfH2+/bfbVqGFa8o/38mTtXe1Y++YmdryzgZMVHsQaNQqrTBno3h1++MHek5Br5LDjRRMSEujQoQPdunWjffv2ABQvXpy4uDh8fX2Ji4ujWLFigGmRHz58OOWxMTExlChRwo6yRUTc2qBBZiK5JUvM6DaA2bMBHEBdoC7l2c/4ku/SJmImfPKJmc1mwABo3x5y2BIpcoXTW+iWZdGrVy+CgoIYOHBgyvY2bdow2/zkMHv2bNq2bZuyPTw8nPj4eA4dOkRUVBS1atVydtkiItnCq6+aPnEXL8L69aaz+9UfBwigbfQkZo+KgQkT4OhRc5+9XDkYO9ZcwxdbOKzUrmlnou+++4569eoRHByMh4f5/8Sbb75J7dq16dSpE7/99hulS5dm/vz5FCpUCIBRo0Yxc+ZMcuTIwcSJE2nRokW6rxMSEkJkZGSmnouISHazf79ZoRXM8LcGdZOodGgp9bZPJPDwGi7nyE1kpce49OSzNO5fyd5i3dCNss3pge4sCnQRkcyxZYsZqn75shnr/pfKST/xdMI7dE6cSy7iORzYiMu9+lLgsbYU8fWyr2A3cqNs00xxIiJyS+rUMVfWz5+Hc+f+/tjyZ1XCEqbz3aeHeZE3SYo6SPlhHblcogyz/V/l9M+aXjYzKdBFRCRDPdS1KE8eeJGfFx3gow5fsi93dXpEv0H+qmU5Wb+9GQifnGx3mW5HXRJFRCTDlSsH5cp5wiOtsKxWdAo5SMj2D+i5YSaELuJ8yXvI0e/fJHR7HO66K+Vx+fJdexlfbp5a6CIikqkcDpi/rRxFpo+lFIfpzsfsjC1Crv8MxLNMSeb79KSpz2Z8fCyKFTOT0v3+u91Vux51ihMREac6cAAWL4aisT9SbfNUgnZ8Ss748xzIVYnJl3rzMT34nSI0awZeXmbK2lGjNMwd1MtdRESysnPnzCIw06fDli0kenix2qcdi4v05v39TbDwwNMTBg+GatXMQypUMDPaZTcKdBERcQ27dsGMGWZZ11OnSC5Vhs/z9WLIL48TQ6lrDu3TBx56CDp2tKlWG2jYmoiIuIYqVcwMdEeOQHg4HhUC6frLK/zmUZZz9R8mZtIXfD4nHoBp06BTJ3OP/q+Pzz/Pvsu4q4UuIiJZ26FDMHOm+ThyBAoVIqljF062DGPajvtISnZw8SKMG2cOL1kS6tc3wd6+Pdx3n1ny3R3okruIiLi+xEQzhn3OHFi0CC5dgooV4bHHoHt3vv21FM8/bzYnJcHBg38/NCQE8uY1Xz//PHj84/r0/fdDkSLOO5XbpUAXERH38scfMH++CfcNG8z19saNISzMNMvz5mXnTvjuOwgPB09Pc3v+RsPh6teHbt2gaFFo1855p3IrFOgiIuK+Dh6Ejz824X7woGmK/+tfJtzr1zdpjrkEv3s3xMdf+/ClS83l+gsX/t726KNmEZoCBczqsP9s0dtFgS4iIu7PsmDjRrOI+7x5cPYslChhlnd99FGoWTPNaeiSk81KsPv3m4Z+UtLf+zp1Mp3tsgIFuoiIZC8XL8KXX8Jnn8GyZWZpuIAA6NrVfAQFpfsUiYlmYhswDX0fH7j3XvN/g4oVM7n+NGjYmoiIZC+5c5um9aJFcOyYGdtepoyZcq5SJZPM48bBb7+l+RQ5csD69WYym/Pn4auv4I03oGpVMxdOVqNAFxER91awIPTsaXrIx8TAxIng7Q1Dh5qQr1sX3nsPjh+/7qH16sGPP8K2beb/Ba+8AgkJ5mEffGBa8VmFLrmLiEj2dOCA6QL/2Wemt5yHBzRoYDrUtW8Pd9993UNOnTJD3H799drts2ZBs2am/12xYplXsu6hi4iI3MjPP5uOdPPnw759pvNcvXp/h3vJktccfvasmanu9Gl4881rn6pMGXjgAfO1ry/8978ZV6YCXURE5GZYFuzZY4J9wQLTcgd48EET7h06QKlr55Rfs+bvFvs775hL8n/x94cVKzKuPAW6iIjI7fjlF/jiCxPwP/1kttWu/Xe4+/s7tRz1chcREbkdQUHw8suwc6e5FP/mm6YJ/sILUK6c6QL/6quwY4ftq8Io0EVERG7GPffAiy+aLu8HDpib4wUKmLFsNWqYFWAGDDDX4G3o/q5AFxERuVXlysHAgWag+tGjZpx7tWpmLFvjxlC8uJl6duFCp7XcFegiIiJ3olgxM859yRI4edLcc2/Z0sxU9/LLaU43m9FyOOVVREREsoN8+cwwt/btzb32mBinvbRa6CIiIpnBy8upveAV6CIiIm5AgS4iIuIGFOgiIiJuQIEuIiLiBhToIiIibkCBLiIi4gYU6CIiIm5AgS4iIuIGFOgiIiJuQIEuIiLiBhToIiIibkCBLiIi4gYU6CIiIm5AgS4iIuIGFOgiIiJuQIEuIiLiBhToIiIibsBhWZZldxGZoUiRIpQtWzbDnu/EiRMULVo0w57PLjqPrMMdzgF0HlmNziNryejziI6O5uTJk6nuc9tAz2ghISFERkbaXcYd03lkHe5wDqDzyGp0HlmLM89Dl9xFRETcgAJdRETEDSjQb1KfPn3sLiFD6DyyDnc4B9B5ZDU6j6zFmeehe+giIiJuQC10ERERN6BAv8ry5cupUKECAQEBjBkz5rr9lmXx7LPPEhAQQNWqVdm+fbsNVd7Y4cOHadSoEUFBQVSuXJlJkyZdd8zatWspUKAA1atXp3r16rz++us2VJq+smXLEhwcTPXq1QkJCbluvyu8H/v27Uv5d65evTo+Pj5MnDjxmmOy6vvRs2dPihUrRpUqVVK2nTp1itDQUAIDAwkNDeX06dOpPja93yVnSu08XnjhBSpWrEjVqlVp164dZ86cSfWx6f0MOlNq5zFixAhKliyZ8rOzbNmyVB+b1d+Pzp07p5xD2bJlqV69eqqPzUrvR1p/a239HbHEsizLSkxMtMqVK2cdOHDAio+Pt6pWrWrt3r37mmOWLl1qNW/e3EpOTrY2b95s1apVy6Zq03bkyBFr27ZtlmVZ1tmzZ63AwMDrzmPNmjVWy5Yt7SjvlpQpU8Y6ceJEmvtd4f24WmJiolW8eHErOjr6mu1Z9f1Yt26dtW3bNqty5cop21544QVr9OjRlmVZ1ujRo60hQ4Zc97ib+V1yptTOY8WKFVZCQoJlWZY1ZMiQVM/DstL/GXSm1M7j1Vdftd56660bPs4V3o+rDRw40HrttddS3ZeV3o+0/tba+TuiFvoVW7duJSAggHLlypEzZ066dOlCRETENcdERETw2GOP4XA4qFOnDmfOnCEuLs6milPn6+tLjRo1AMifPz9BQUHExsbaXFXmcIX342rffvst5cuXp0yZMnaXclPq169PoUKFrtkWERFBWFgYAGFhYSxevPi6x93M75IzpXYeTZs2JUeOHADUqVOHmJgYO0q7Jamdx81whffjL5ZlMW/ePLp27erkqm5dWn9r7fwdUaBfERsbS6lSpVK+9/Pzuy4Ib+aYrCQ6OpodO3ZQu3bt6/Zt3ryZatWq0aJFC3bv3m1DdelzOBw0bdqUmjVrMm3atOv2u9r7ER4enuYfKld4PwCOHTuGr68vYP6gHT9+/LpjXO19mTlzJi1atEh1X3o/g1nB5MmTqVq1Kj179kz18q4rvR8bNmygePHiBAYGpro/q74fV/+ttfN3JMcdP4ObsFLp7O9wOG75mKzi/PnzdOjQgYkTJ+Lj43PNvho1avC///2PfPnysWzZMh555BGioqJsqjRtGzdupESJEhw/fpzQ0FAqVqxI/fr1U/a70vtx+fJllixZwujRo6/b5yrvx81ypfdl1KhR5MiRg27duqW6P72fQbs9/fTTDB8+HIfDwfDhwxk0aBAzZ8685hhXej8+++yzG7bOs+L7caO/tWnJrPdELfQr/Pz8OHz4cMr3MTExlChR4paPyQoSEhLo0KED3bp1o3379tft9/HxIV++fAA8/PDDJCQkpDk3sJ3++rctVqwY7dq1Y+vWrdfsd5X3A+Drr7+mRo0aFC9e/Lp9rvJ+ABQvXjzltkZcXBzFihW77hhXeV9mz57NV199xSeffJLmH9P0fgbtVrx4cTw9PfHw8ODJJ59MtT5XeT8SExNZuHAhnTt3TvOYrPZ+pPa31s7fEQX6Fffddx9RUVEcOnSIy5cvEx4eTps2ba45pk2bNsyZMwfLstiyZQsFChRIubSSVViWRa9evQgKCmLgwIGpHnP06NGU/yFu3bqV5ORkChcu7Mwy03XhwgXOnTuX8vU333xzTa9YcI334y83anm4wvvxlzZt2jB79mzABGLbtm2vO+Zmfpfstnz5csaOHcuSJUvIkydPqsfczM+g3a7uM7Jo0aJU63OF9wNg1apVVKxYET8/v1T3Z7X3I62/tbb+jtxxtzo3snTpUiswMNAqV66cNXLkSMuyLGvq1KnW1KlTLcuyrOTkZKtv375WuXLlrCpVqlg//PCDneWmasOGDRZgBQcHW9WqVbOqVatmLV269JrzePfdd61KlSpZVatWtWrXrm1t3LjR5qqvd+DAAatq1apW1apVrUqVKrns+2FZlnXhwgWrUKFC1pkzZ1K2ucL70aVLF+vuu++2cuTIYZUsWdKaPn26dfLkSatx48ZWQECA1bhxY+v333+3LMuyYmNjrRYtWqQ8NrXfJbukdh7ly5e3/Pz8Un5HnnrqKcuyrj2PtH4G7ZLaeXTv3t2qUqWKFRwcbLVu3do6cuSIZVmu935YlmWFhYWl/E78JSu/H2n9rbXzd0QzxYmIiLgBXXIXERFxAwp0ERERN6BAFxERcQMKdBERETegQBcREXEDCnQRERE3oEAXyUZ+//33lGUq77777pSlN/Ply0ffvn0z5TUnTpzInDlzrtseHR2dMjHIzz//zOOPP54pry+SXWgud5FspHDhwvz444+AWUs7X758DB48ONNeLzExkZkzZ6a7Vn1wcDAxMTH89ttvlC5dOtPqEXFnaqGLCGvXrqVVq1aACfqwsDCaNm1K2bJlWbhwIUOGDCE4OJjmzZuTkJAAwLZt22jQoAE1a9akWbNmqS5du3r1amrUqJGyVOm2bduoVq0a999/P1OmTLnm2NatWxMeHp7JZyrivhToInKdAwcOsHTpUiIiIujevTuNGjXi559/Jnfu3CxdupSEhASeeeYZFixYwLZt2+jZsycvvfTSdc+zceNGatasmfL9E088wTvvvMPmzZuvOzYkJIQNGzZk6nmJuDNdcheR67Ro0QIvLy+Cg4NJSkqiefPmgLk0Hh0dzb59+9i1axehoaEAJCUlpbowTlxcHEFBQQD88ccfnDlzhgYNGgDQo0cPvv7665RjixUrxpEjRzL71ETclgJdRK7j7e0NgIeHB15eXinLi3p4eJCYmIhlWVSuXDnVlvbVcufOzaVLlwCzOtWN1ny+dOkSuXPnzqAzEMl+dMldRG5ZhQoVOHHiREqgJyQksHv37uuOCwoKYv/+/QAULFiQAgUK8N133wHwySefXHPsr7/+muWWJxVxJQp0EbllOXPmZMGCBQwdOpRq1apRvXp1Nm3adN1xLVq0YP369Snfz5o1i379+nH//fdf1xpfs2YNLVu2zPTaRdyVlk8VkUzVrl07xo0bR2BgYJrHxMfH06BBA7777ruUHvEicmsU6CKSqfbt28exY8eoX79+msdERUURGxtLw4YNnVeYiJtRoIuIiLgB3UMXERFxAwp0ERERN6BAFxERcQMKdBERETegQBcREXED/wfVnJttvyijjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 566.929x566.929 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Specify a probabilistic mortality calculator\n",
    "config.set('BIO_MODEL', 'mortality_calculator', 'probabilistic')\n",
    "\n",
    "# Set the death rate - currently the same for all particles.\n",
    "death_rate_per_day = 0.1\n",
    "config.set('PROBABILISTIC_MORTALITY_CALCULATOR', 'death_rate_per_day', str(death_rate_per_day))\n",
    "\n",
    "# Set the bio time step\n",
    "config.set('NUMERICS', 'bio_time_step', str(time_step))\n",
    "\n",
    "# Number of particles\n",
    "n_particles = 1000\n",
    "\n",
    "# Run the model\n",
    "n_alive_numeric = run(config, n_particles=n_particles)\n",
    "\n",
    "# Compute the equivalen analytical solution\n",
    "death_rate_per_second = death_rate_per_day / seconds_per_day\n",
    "n_alive_analytic = n_particles * np.exp(-death_rate_per_second * times)\n",
    "\n",
    "# Plot\n",
    "font_size = 10\n",
    "fig, ax = create_figure(figure_size=(20, 20), font_size=font_size)\n",
    "plt.plot(times/seconds_per_day, n_alive_numeric, 'b', label='numeric')\n",
    "plt.ylabel('Living individuals (-)', fontsize=font_size)\n",
    "plt.xlabel('Time (d)', fontsize=font_size)\n",
    "\n",
    "# Add equivalent analytical solution\n",
    "plt.plot(times/seconds_per_day, n_alive_analytic, 'r', label='analytic')\n",
    "\n",
    "# Add legend\n",
    "plt.legend()"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}