problemSet3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem Set 3: Simulating the Spread of Disease and Virus Population Dynamics \n",
    "\n",
    "import random\n",
    "import pylab\n",
    "from numpy import mean\n",
    "\n",
    "#random.seed(0)\n",
    "\n",
    "''' \n",
    "Begin helper code\n",
    "'''\n",
    "\n",
    "class NoChildException(Exception):\n",
    "    \"\"\"\n",
    "    NoChildException is raised by the reproduce() method in the SimpleVirus\n",
    "    and ResistantVirus classes to indicate that a virus particle does not\n",
    "    reproduce. You can use NoChildException as is, you do not need to\n",
    "    modify/add any code.\n",
    "    \"\"\"\n",
    "\n",
    "'''\n",
    "End helper code\n",
    "'''\n",
    "\n",
    "#\n",
    "# PROBLEM 1\n",
    "#\n",
    "class SimpleVirus(object):\n",
    "    \"\"\"\n",
    "    Representation of a simple virus (does not model drug effects/resistance).\n",
    "    \"\"\"\n",
    "    def __init__(self, maxBirthProb, clearProb):\n",
    "        \"\"\"\n",
    "        Initialize a SimpleVirus instance, saves all parameters as attributes\n",
    "        of the instance.        \n",
    "        maxBirthProb: Maximum reproduction probability (a float between 0-1)        \n",
    "        clearProb: Maximum clearance probability (a float between 0-1).\n",
    "        \"\"\"\n",
    "        self.maxBirthProb = maxBirthProb\n",
    "        self.clearProb = clearProb\n",
    "\n",
    "    def getMaxBirthProb(self):\n",
    "        \"\"\"\n",
    "        Returns the max birth probability.\n",
    "        \"\"\"\n",
    "        return self.maxBirthProb\n",
    "\n",
    "    def getClearProb(self):\n",
    "        \"\"\"\n",
    "        Returns the clear probability.\n",
    "        \"\"\"\n",
    "        return self.clearProb\n",
    "\n",
    "    def doesClear(self):\n",
    "        \"\"\" Stochastically determines whether this virus particle is cleared from the\n",
    "        patient's body at a time step. \n",
    "        returns: True with probability self.getClearProb and otherwise returns\n",
    "        False.\n",
    "        \"\"\"\n",
    "        if self.getClearProb() >= random.random():\n",
    "            return True\n",
    "        else:\n",
    "            return False\n",
    "        \n",
    "    \n",
    "    def reproduce(self, popDensity):\n",
    "        \"\"\"\n",
    "        Stochastically determines whether this virus particle reproduces at a\n",
    "        time step. Called by the update() method in the Patient and\n",
    "        TreatedPatient classes. The virus particle reproduces with probability\n",
    "        self.maxBirthProb * (1 - popDensity).\n",
    "        \n",
    "        If this virus particle reproduces, then reproduce() creates and returns\n",
    "        the instance of the offspring SimpleVirus (which has the same\n",
    "        maxBirthProb and clearProb values as its parent).         \n",
    "\n",
    "        popDensity: the population density (a float), defined as the current\n",
    "        virus population divided by the maximum population.         \n",
    "        \n",
    "        returns: a new instance of the SimpleVirus class representing the\n",
    "        offspring of this virus particle. The child should have the same\n",
    "        maxBirthProb and clearProb values as this virus. Raises a\n",
    "        NoChildException if this virus particle does not reproduce.               \n",
    "        \"\"\"        \n",
    "        if (self.getMaxBirthProb() * (1 - popDensity)) >= random.random():\n",
    "            return SimpleVirus(self.getMaxBirthProb(), self.getClearProb())\n",
    "        else:\n",
    "            raise NoChildException\n",
    "\n",
    "\n",
    "class Patient(object):\n",
    "    \"\"\"\n",
    "    Representation of a simplified patient. The patient does not take any drugs\n",
    "    and his/her virus populations have no drug resistance.\n",
    "    \"\"\"    \n",
    "\n",
    "    def __init__(self, viruses, maxPop):\n",
    "        \"\"\"\n",
    "        Initialization function, saves the viruses and maxPop parameters as\n",
    "        attributes.\n",
    "\n",
    "        viruses: the list representing the virus population (a list of\n",
    "        SimpleVirus instances)\n",
    "\n",
    "        maxPop: the maximum virus population for this patient (an integer)\n",
    "        \"\"\"\n",
    "        self.viruses = viruses\n",
    "        self.maxPop = maxPop\n",
    "\n",
    "\n",
    "    def getViruses(self):\n",
    "        \"\"\"\n",
    "        Returns the viruses in this Patient.\n",
    "        \"\"\"\n",
    "        return self.viruses\n",
    "\n",
    "\n",
    "    def getMaxPop(self):\n",
    "        \"\"\"\n",
    "        Returns the max population.\n",
    "        \"\"\"\n",
    "        return self.maxPop\n",
    "\n",
    "\n",
    "    def getTotalPop(self):\n",
    "        \"\"\"\n",
    "        Gets the size of the current total virus population. \n",
    "        returns: The total virus population (an integer)\n",
    "        \"\"\"\n",
    "        return len(self.viruses)\n",
    "\n",
    "\n",
    "    def update(self):\n",
    "        \"\"\"\n",
    "        Update the state of the virus population in this patient for a single\n",
    "        time step. update() should execute the following steps in this order:\n",
    "        \n",
    "        - Determine whether each virus particle survives and updates the list\n",
    "        of virus particles accordingly.   \n",
    "        \n",
    "        - The current population density is calculated. This population density\n",
    "          value is used until the next call to update() \n",
    "        \n",
    "        - Based on this value of population density, determine whether each \n",
    "          virus particle should reproduce and add offspring virus particles to \n",
    "          the list of viruses in this patient.                    \n",
    "\n",
    "        returns: The total virus population at the end of the update (an\n",
    "        integer)\n",
    "        \"\"\"\n",
    "       \n",
    "        surviving_viruses = []\n",
    "        for virus in self.getViruses():\n",
    "            if not virus.doesClear():\n",
    "                surviving_viruses.append(virus)\n",
    "                \n",
    "        # print('printing surviving viruses list in .update():', surviving_viruses)\n",
    "        self.viruses = surviving_viruses\n",
    "        popDensity = self.getTotalPop() / self.getMaxPop()\n",
    "        # print('printint getTotalPop', self.getTotalPop())\n",
    "        # print('printint getTotalPop', self.getMaxPop())\n",
    "        # print('printing popDensity:', popDensity)\n",
    "        virus_offspring = []\n",
    "        \n",
    "        for virus in self.getViruses():\n",
    "            try:\n",
    "                virus_offspring.append((virus.reproduce(popDensity)))\n",
    "            except NoChildException:\n",
    "                continue\n",
    "                \n",
    "        self.viruses += virus_offspring\n",
    "        \n",
    "        return self.getTotalPop()\n",
    "\n",
    "#\n",
    "# PROBLEM 2\n",
    "#\n",
    "def simulationWithoutDrug(numViruses, maxPop, maxBirthProb, clearProb,\n",
    "                          numTrials):\n",
    "    \"\"\"\n",
    "    Run the simulation and plot the graph for problem 3 (no drugs are used,\n",
    "    viruses do not have any drug resistance).    \n",
    "    For each of numTrials trial, instantiates a patient, runs a simulation\n",
    "    for 300 timesteps, and plots the average virus population size as a\n",
    "    function of time.\n",
    "\n",
    "    numViruses: number of SimpleVirus to create for patient (an integer)\n",
    "    maxPop: maximum virus population for patient (an integer)\n",
    "    maxBirthProb: Maximum reproduction probability (a float between 0-1)        \n",
    "    clearProb: Maximum clearance probability (a float between 0-1)\n",
    "    numTrials: number of simulation runs to execute (an integer)\n",
    "    \"\"\"\n",
    "     #Ensuring numTrials stays under 100 (as directed in problem outline)\n",
    "    if numTrials > 100:\n",
    "        raise Exception('numTrials must be less than 100')\n",
    "        \n",
    "    #creating initial viruses\n",
    "    virus_list = []\n",
    "    for i in range(numViruses):\n",
    "        virus_list.append(SimpleVirus(maxBirthProb, clearProb))   \n",
    "    \n",
    "    #creating patient list\n",
    "    patient_list = []\n",
    "    for patient in range(numTrials):\n",
    "        patient_list.append(Patient(virus_list, maxPop))\n",
    "    \n",
    "    #creating x and y values to be plotted\n",
    "    timesteps = range(300)     # x-values\n",
    "    avg_pop_per_timestep = []  # y-values\n",
    "    \n",
    "    #creating list of avg population size per timestep\n",
    "    for timestep in timesteps:\n",
    "        virus_pops_per_timestep = []\n",
    "        for patient in patient_list:\n",
    "            patient.update()\n",
    "            virus_pops_per_timestep.append(patient.getTotalPop())\n",
    "        avg_pop_per_timestep.append(mean(virus_pops_per_timestep))\n",
    "    \n",
    "    # print(avg_pop_per_timestep)\n",
    "    # pylab.plot(timesteps, avg_pop_per_timestep)\n",
    "    pylab.plot(avg_pop_per_timestep, label = \"SimpleVirus\")\n",
    "    pylab.title(\"SimpleVirus simulation\")\n",
    "    pylab.xlabel(\"Time Steps\")\n",
    "    pylab.ylabel(\"Average Virus Population\")\n",
    "    pylab.legend(loc = \"best\")\n",
    "    pylab.show()\n",
    "    \n",
    "\n",
    "#simulationWithoutDrug(100, 1000, abs(random.random()-0.5), random.random(), 100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA980lEQVR4nO3deXxU1fnH8c83ISSEnRCQfUfZFDDiiuKK4oK74lK3Vtva1tZaq6222pbWtfVnW22xWmjVUtqqRRRkqbiBICD7vhMIIYQ1LFmf3x/3EgMkYQJMJsvzfr3mlZkzd3nuXJhnzjn3niMzwznnnAOIi3UAzjnnqg5PCs4554p5UnDOOVfMk4JzzrlinhScc84V86TgnHOumCcFVykk3SppYpS2PVLSr47TthZJGnQ8tnWU+x8oaVmUtv2EpNePYf2YfjaucnhScMeNpHMkTZO0U9I2SZ9JOg3AzN4ws0tiHN8Hkn5RSvlQSZsl1TGzXmY2NQbhAWBmn5jZibHa/wGlJdpYfzaucnhScMeFpEbAOOD3QDOgDfAkkBvLuA4xErhdkg4pvx14w8wKyltZUp1oBeZcVeFJwR0v3QHM7B9mVmhm+8xsopnNB5B0p6RPDywsySR9W9IKSbsl/VJSF0nTJe2SNEZS3XDZQZLSJf1E0lZJayXdWlYgkq6QNFfSjrDmcnL41jsECWtgiWWbAlcAfwtfr5V0Ufj8CUn/lvS6pF3AnYf+gj4QW4nXP5a0MTymZZIuLCPGIZIWh8ttlPRQGdtbK+lHkuZL2iPpVUktJY0P150cHsNh6x56PKXE8K+whrRT0seSeoXl9wK3Ag9LypH0bimfTaKkFyRtCh8vSEo85Hz9UNIWSRmS7irrfLmqxZOCO16WA4WSRkm67MAX1RFcCpwKnAE8DIwg+DJqB/QGhpVY9gSgOUEN5A5ghKTDmlkk9QdeA+4DUoA/A2MlJZrZPmAM8LUSq9wILDWzeWXEOBT4N9AEeKO8gwnj+Q5wmpk1BAYDa8tY/FXgvnC53sD/ytn0dcDFBIn3SmA88BOCzyMO+F55cZVjPNANaAHMITw+MxsRPn/GzBqY2ZWlrPtTgvPWFzgFGAA8VuL9E4DGBOfrHuCPEf6bcDHmScEdF2a2CzgHMOAVIEvSWEkty1ntaTPbZWaLgIXARDNbbWY7Cb6w+h2y/ONmlmtmHwHvEXyhH+obwJ/NbEZYYxlF0IR1Rvj+KOAGSfXC118Ly8oy3czeMbOiMKmUpxBIBHpKSjCztWa2qoxl88PlGpnZdjObU852f29mmWa2EfgEmGFmX5pZLvA2h39OETGz18xsd7idJ4BTJDWOcPVbgV+Y2RYzyyJoKry9xPv54fv5ZvY+kAPEvK/EHZknBXfcmNkSM7vTzNoS/PptDbxQziqZJZ7vK+V1gxKvt5vZnhKv14XbP1QH4Idh09EOSTsIah6twxg/BbKAoZI6A6cBb5YT44Zy3juIma0Evk/wBbtF0mhJpcUIwa//IcA6SR9JOrOcTVfkc4qIpHhJT0laFTaNrQ3fah7hJloTnIMDDj0f2Yf00ew9mjhd5fOk4KLCzJYSdOz2Pk6bbCqpfonX7YFNpSy3ARhuZk1KPJLN7B8llvkbQQ3hdoLaSWYp2zng0GGE9wDJJV6fcNDCZm+a2TkEycmAp0vdqNkXZjaUoOnmHYJmrWN1UGyS4oHUMpa9haBp7CKCZp6OB1Y7EOIR9rWJ4BgPKOt8uGrGk4I7LiSdFHYstg1ftyPoE/j8OO7mSUl1JQ0k6Bz+VynLvAJ8U9LpCtSXdLmkhiWW+RvBl+E3KL/pqDRzgSGSmkk6gaBmAAR9CpIuCDtc9xP8ii88dAPhMdwqqbGZ5QO7SlvuKCwHksLjTSBo408sY9mGBM1q2QSJ5NeHvJ8JdC5nX/8AHpOUKqk58DPgqO+BcFWHJwV3vOwGTgdmSNpDkAwWAj88TtvfDGwn+DX6BvDNsDZyEDObRfBl/4dw+ZXAnYcssxaYBtQHxlYwjr8D8wiaWyYC/yzxXiLwFLA1jLcFQYdwaW4H1oZNN98EbqtgHIcJ+2K+DfwF2EhQc0gvY/G/ETT5bAQWc3jyfpWgz2OHpHdKWf9XwCxgPrCAoKP6uNxA6GJLPsmOq+oU3EX7ethX4ZyLIq8pOOecK+ZJwTnnXDFvPnLOOVfMawrOOeeKVesBvpo3b24dO3aMdRjOOVetzJ49e6uZlXoPS7VOCh07dmTWrFmxDsM556oVSevKes+bj5xzzhXzpOCcc66YJwXnnHPFqnWfQmny8/NJT09n//79sQ6l1klKSqJt27YkJCTEOhTn3FGqcUkhPT2dhg0b0rFjR3TYrIsuWsyM7Oxs0tPT6dSpU6zDcc4dpRrXfLR//35SUlI8IVQySaSkpHgNzblqrsYlBcATQoz45+5c9Vcjk4Jz7ujs2JvHjr15sQ6jRissMv7++TqmrdpKVRxmyJNClAwfPpxevXpx8skn07dvX2bMmMHXv/51Fi9efFy236BB+TMbDho0iA8++OCgshdeeIFvf/vbjB07lqeeeuq4xOGqv0WbdvLmjPU89K95nPGbKVz+4qfs2p8f67BqrBEfr+bxdxZyyysz+L8pK/jf0swqlYhrXEdzVTB9+nTGjRvHnDlzSExMZOvWreTl5fGXv/yl0mIYNmwYo0ePZvDgwcVlo0eP5tlnn2XgwIFcddVVh61TUFBAnTr+T6KmyisoYmtOLic0SmJPXgEJ8XE8MXYRo78IpqFulFSHi3q0ZPzCzVzw3FQGndiC5244pdRtmRmvfbaW6au28vth/alXN74yD6Xa+nx1Nr+btJzBvVpSt048L0xeAUCn5vV55vqTmbV2O+d0bU6fto0B2L4nj/qJdahbp/J+v0f1G0DSWoIZuQqBAjNLk9SMYLaqjgSzV91oZtvD5R8F7gmX/56ZfVDKZqu8jIwMmjdvTmJiMBNi8+bBXOiDBg3iueeeIy0tjQYNGnD//fczefJkmjZtyq9//Wsefvhh1q9fzwsvvMBVV13FyJEjefvtt8nNzWXNmjXccsst/PznPz9sf88++yxjxowhNzeXa665hieffJLrr7+exx57jNzcXBITE1m7di2bNm3inHPOYeTIkcyaNYs//OEP3HnnnTRr1owvv/yS/v3707BhQxo0aMBDDz0EQO/evRk3bhypqanceOONpKenU1hYyOOPP85NN91UeR/qUdqyaz9Tl2Vx3omptGyUFOtwYmbnvnyueekzVmftYUDHZqRv38v+giK27cnjzrM68o1zO9OyYSJ14uN4Y8Y6Xv98Pf+enc4PLu5Omyb1WLt1D89PWs6KzN10SElm4459LNy4C4Bfv7+EX159vKbirrnWZ+/lnpFf0D4lmaeuPZn4eFFUZJx0QkNGTV/LDX+aDsDv4uO4qGcLzuicwtPjl5LaMJF2zZK56pTW3JDWDoB12XuoVzeeFg2P/7/pyvhZeL6ZbS3x+hFgipk9JemR8PWPJfUEbgZ6Aa2ByZK6m9lRz1375LuLWLxp17HEfpierRvx8yt7lbvMJZdcwi9+8Qu6d+/ORRddxE033cR555130DJ79uxh0KBBPP3001xzzTU89thjTJo0icWLF3PHHXcU/5KfOXMmCxcuJDk5mdNOO43LL7+ctLS04u1MnDiRFStWMHPmTMyMq666io8//phzzz2XAQMGMGHCBIYOHcro0aO56aabSu0MXr58OZMnTyY+Pp4nnnii1GOaMGECrVu35r333gNg586dFfnYYmLxpl1c9/I09uUX0jm1Pt8e1JUzOjfjv3M3sTevgOS6dUhKiGf7njzO7tqcM7ukRDWeLbv2Ex8nUhqUNW1y6fILi/hw6Ra6tGjAr8Yt5hsDO3NW1+YRr79zbz4PjpnLuuy93HNOJ179dA1JCXE0rpfAgI7NePyKnsTHffXv4tbTOzCwayrnPvshv3x3Me1Tkvn37HTyCoro264Jq7P20LxBIj8d0oONO/Yxctpa7jirA11bfDUNdmGRsTUnt9IS8c69+Yz4ZBWnd0qhTpxIqBPHaR2bVcq+I/Wb8UsoMhh19wCa1q8LwB9v7Q/A187qyNtz0unSogHvztvEx8u38v6CzbRunERy3Tos3rSLz1dn0zS5LnFx8MMx8zilXRNG3jXguMcZi7aCocCg8PkoYCrw47B8tJnlAmskrQQGANNjEOMxadCgAbNnz+aTTz7hww8/5KabbjqsDb9u3bpceumlAPTp04fExEQSEhLo06cPa9euLV7u4osvJiUl+LK69tpr+fTTTw9LChMnTqRfv34A5OTksGLFCs4999ziJqQDSeG1114rNd4bbriB+Pjyq/99+vThoYce4sc//jFXXHEFAwcOrPDncrxs2LaXlo2SmLwkk7O7NufDpVtIrhvPRT1a8v7CDOas28H6bXuYl76Thkl1+M21fXjkrfk89K95SGAG8XGisOirTr6XP1rFr67uzc2ntePTlVtJbZjISSc0Om4xPz1hKS9PXUV8nLjv3M4M6dOK3m0aH3G9Lbv3c+drX7A4Y1dxzDPWbOM7F3Tl8j6t6JBSv9z1i4qM6/80jVVZOfzsip7ceXYnTmnXhBYNE+nXvglCByWEA9qnJNO/fRMmLNpMnKB9s2RG3jWAjs0P3l92Ti7/mLmeVz5eQ1rHpvz1s7U8fkVPnnx3ESu25DDmvjM5tUPTiD6jgsIi6sQf3ExiZkhi5958/m/KCk5q1ZBr+7U5aLms3blc89JnpG/fxztfbiJ7Ty7784vo06Yxp3Zoyq2nt6dby4aH7q5SvTFjHeMXbubBsOZ1qMb1Erjz7OD+noHdUskvLGLc/E2c2r4Z7VOS2bE3j6F//Iyv/y0YALRz8/o8cYQfp0cr2knBgImSDPizmY0AWppZBoCZZUhqES7bhoMnD08Pyw4i6V7gXoD27duXu/Mj/aKPpvj4eAYNGsSgQYPo06cPo0aNOuj9hISE4l/tcXFxxU1NcXFxFBQUFC936C/7Q1+bGY8++ij33XffYTFcffXVPPjgg8yZM4d9+/bRv3//UmOtX/+r/+h16tShqKio+PWB+w66d+/O7Nmzef/993n00Ue55JJL+NnPfnbEz+F4+8fM9Tz61gLaN0tm/ba9JCXEsT8/iHfQialMXZZFvYR4OqQkc9IJDXnw4u70a9+U809swYbte3l56irOP6kFV57SCjPYl1eIBA+Mnsujby3gxSkryNi5n8Q6cfzhlv5c3LNlRHEVFRnLt+xmb14hW3blkr0nl1sGtEcS783P4OWpqxjatzVm8NLUVfzpo1W8972B9Gh1eOIpLDI+XLqFUzs05duvz2H11hx+NPhE/jt3I/ec04n/zN7IMxOW8cyEZfxiaC++dmbHg9Y/8EVaUFjErHXbWbElh2euP5kbw6aHq05pHdExPXfDKazbtpeBXZsTJxFXSvJIaZDIjWnteH3GOv4zJ52CImPYK5/TvEFdUurX5eujviBO4uKeLbnvvC6M/mI99w7sfFhtadGmndw84nO+eV4XWjZKYnCvluTkFnDZ/31CQnwcufmF7Nof/L9YnbWHRy47CTPjL5+s4fUZ68jOyePeczsz4uPVANxzTieWZOzizZnrGf3Fen53Y18u69MqouM+FhMWZvDS1FX0b9+Uhy89keS6dZi2ais/fXshF5zUgm+e1yWi7STEx3FNv6+mJG+SXJcJD5zLlKWZ1ImLY2C35tRPjM7Xd7STwtlmtin84p8kaWk5y5Z2kfth12uFiWUEQFpaWtW7ngtYtmwZcXFxdOvWDYC5c+fSoUMHFi5cWOFtTZo0iW3btlGvXj3eeeedw37tDx48mMcff5xbb72VBg0asHHjRhISEmjRogUNGjRg0KBB3H333QwbNiyi/XXs2JFx48YBMGfOHNasWQPApk2baNasGbfddhsNGjRg5MiRFT6WY7Uuew+/GreYNk3qsWH7XoYNaM+8DTuCL8o56UxdlsUJjZL4+OHzD+uYa5ycQOPkxsXV9QOSEoIa0l/uSOPFKStYm72XAZ2a8e/Z6XznzTn88Zb+tG1Wjz25BTRJrkv7Zsn85ZM17Mkt4Nvnd2Hzzv0UFhk/H7uIaauyD9r2lCVbSIgXn6/eRt92TXjuhlNIiI/jh5d05+LffcybM9aX2hb/xox1/Oy/i0iIFwVFxos39+PKU1pz//ldAbjptPakb9/LY+8s5FfjltC1RQM+X5XNhT1a0jS5LsNe+ZzWTZJYtnk3BiTXjeeKkyv+hdg5tQGdU8u/yg3gkctOIiE+jvnpOxh0Yiojp61j5F2nkZWTy8/+u5DerRvzz1kbePvLjeQWBE1h3Vs2JE7iw6VbQNC8QSK79xfw7AfLAJi9rh1FRbA3t5ArT2lJnThx04B2jJq2lpHT1vD1gZ2YtDiT4e8v4ZR2Tfjl0N6c0TmFsXM30btNYx6/oicQ1CLu+/ssvjf6S35bZHRISSalQWKpv9aP1Wcrt/KtN+bQrmkyI6etxcw4s0sKz09cTtum9Xjp1v7H1GFcr248V5wcWUI/FlFNCma2Kfy7RdLbBM1BmZJahbWEVsCWcPF0oF2J1dsCm6IZX7Tk5OTw3e9+lx07dlCnTh26du3KiBEjuP766yu8rXPOOYfbb7+dlStXcssttxzUdARB/8WSJUs488wzgaDp6vXXX6dFi6ACNmzYMK699lpGjx4d0f6uu+46/va3v9G3b19OO+00unfvDsCCBQv40Y9+RFxcHAkJCbz88ssVPpZjsT+/kG++Poc68XGM+eaZNK6XQIMSv5RaN6nHtFXZ3H1Ox6P6jxd8WZ9Y/Pqy3icw9A9fVddL8+ePV5FfGPwuiY8Tj13eg07N6xMfJ8bO28TYuZtoklyX/MIinr8xSAgAHVLqc0WfVoyZtYEN2/fykyE9+M+cdOZv2Enmrv1s35tH95YNSG2YyH3nduHc7ofPhdK2aTK/vbEvQ//4Kbe8MgOAv3y6hnoJ8RQUGRk799OtZUPmbtjB1X3bkFw3ev/V6yfW4WdX9ix+ff/5XYtrtJ88fAEAf/jfCp6buJyvndmBSYszmZ++k5zcAi7u1ZLMXfv5bGU2D17cnfqJdVi4cSf/mBlcEXX32Z0O2nbT5Lq8O28TV//xMzbt2MfAbs0ZedeA4iawsd89+6BjTW2YyF/vHMC1L3/Gd//xZXF5Sv26NE5O4LLeJ5DWoRk9Wzc65v6PF6es4IRGSUz8wbk8MXYRo6avY9T0ddSJE3+67dTiHyBVXdTmaJZUH4gzs93h80nAL4ALgewSHc3NzOxhSb2ANwkSR2tgCtCtvI7mtLQ0O3SSnSVLltCjR4+oHFNlK3mVUHURrc//5amreHrCUl67M40LTiq9SWfhxp30aNWo1Dbyo7Fzbz6z129jf34RyXXj2bxzP5t37ad368Zs3LGPRZt2ckbnFOLjRJfUBgf1EZgZ+/ILqRsfx57cQhonHzxI4IrM3fz6/SVMW5Vd3LfRs3UjkhLimbV2G298/YyIOr637NrPk+8uDvpWlm0hXuI7F3QtjmXt1j2kNkyMWlNDRWTtziW14eGd7EVFxhdrt5HWsRnxcWLX/nx+9K959GnTmK8P7HzYl+n/lmbyzIRl9O/QlJ8M6XHQj4Oy7M0rYM66HezNK2B55m427tjHhm37mL46+PybJicw+cHzWJu9h8IiOK1jU5Zn5vDhsi2c07U5iXXiaJycUHy1T9buXP41ewOtGidxWe9WzN2wg5tHfM5jl/fg6wM7s3NfPq9+spqzujanW4sGFb64INokzTaztFLfi2JS6Ay8Hb6sA7xpZsMlpQBjgPbAeuAGM9sWrvNT4G6gAPi+mY0vbx+eFKqeaHz+e3ILOOfp/3Fy2yaMuvv4X20RS//8Yj1PjF3Mi8P6Ffdf7MktqBJf4rXBtj15fLF2G996fTZJCfHszQt+g57StjFbc/LYuGNf8bLxceLeczvzo0tO5FtvzOaDRZlAUBuJl6gTLz74/rnV4tzFJClUhpqeFKqjY/38N2zby7RVWxk3P4NOzevzkyE9eH7iMl75ZA1vf/ss+rWP7EqW6iS/sKi4acnFxq/fX8LERZv53oXdyC8s4tkPlrFzXz4v3XoqeQVFGMbUZVn8e3Y6A7s155MVW3nw4u6kdWzKi1NWMGfdDt74xulV7jLYspSXFKp+SjsKB66+cJXr0B8YRUVG5u79JMTHBe3n2/bxyGUn0bhe6fMt/Gb8Ev78UXD1SMeUZD5ZsZXsnDzGL8zg1tPb18iEAHhCqAJ+MqQHPxny1Y+ZS3u1Iisnl64tvupov7xPK4rMeGvORgb3asl953UmsU48Z3Vpzr68whpzV3eNSwpJSUlkZ2f78NmV7MB8CklJQZtrXkERX3ttBp+v3la8TJxg2qqtjCrlevfxCzL480erubZ/G67r35azuqTw9VGzeG9BBie2bMgjl51UqcfjarfgarWDf7xI4tnrT+H7F3anfUryQe/VlIQANbD5yGdei50DM69l5uTz+DsL+XBZFt8a1IXUBon0bd8EM+Pro2ZRUGh0bdmAddl7ubT3CTx62Unc8KfpFBYZ4x8YWHxj0rrsPfzpo1X84OLuUbmd37naqlb1KbjKt3TzLv5v8goaJtXhjrM6cvfIL9i9v4AfDT6Ru84+eBa2tVv38OwHy1iVlUPXFg0Yv3AzTZMT2JqTd9ANVs656Kl1fQqu8uQXFnHPyFns2p/P7v0FjJmVTtPkBN7+9tmceMLhQwt0bF7/oBvIpq/K5tVPVxMnMbRv9G/Mcc6Vz5OCOyZj525i4459vHpHGvFxYuWWHC4/uRWtGkd2x+iZXVKiPhCdcy5ynhTcUSsqMl7+aBUnndCQC05qgSQGndjiyCs656osvxbOHbXJSzJZuSWHbw3q4ld6OVdDeFJwRyWvoIjfTlpOu2b1uLwSRp90zlUOTwruqPzxw5Us3bybxy/vedgY+M656sv/N7sKyy8sYuS0tQzu1ZJLep0Q63Ccc8eRJwVXYZ+t3MrOffnccKrfU+BcTeNJwVXYuPkZNEysw8Dukc8T7JyrHvySVBexWWu3sWDjTt6ak86wAe1JrFNzxntxzgU8KbiIZOfk8rXXZrI3r5BOzev7AHXO1VBRbz6SFC/pS0njwtdPSNooaW74GFJi2UclrZS0TNLgaMfmIpO1O5ffjF/K/vxC/nRbf/557xk0TCp9+GvnXPVWGTWFB4AlQKMSZb8zs+dKLiSpJ3Az0ItgOs7JkrqXNx2ni678wiK+9focJi8JZpi6/YwOXNrb70lwriaLalKQ1Ba4HBgOPHiExYcCo80sF1gjaSXBfM3ToxmjK9vvp6xg8pJMvnleF648pRU9WzU68krOuWot2s1HLwAPA0WHlH9H0nxJr0k6MJ1WG2BDiWXSwzIXAy9MXs6L/1vJdf3b8shlJ9GrdWMfysK5WiBqSUHSFcAWM5t9yFsvA12AvkAG8PyBVUrZzGGTPUi6V9IsSbOysrKOY8QOYNGmnYxfkMELk1dwTb82PHVdn1iH5JyrRNFsPjobuCrsSE4CGkl63cxuO7CApFeAceHLdKDk3VBtgU2HbtTMRgAjIJhkJ0qx10oLN+5k6B8/o7DIaNU4ieHX9Pb5g52rZaL2P97MHjWztmbWkaAD+X9mdpukkj2V1wALw+djgZslJUrqBHQDZkYrPnewwiLjR/+eT0r9unztzA48f+MpJNf1K5adq21i8b/+GUl9CZqG1gL3AZjZIkljgMVAAXC/X3lUed7+ciNLMnbxh1v6ccXJPgOac7VVpSQFM5sKTA2f317OcsMJrlRylcTM+NNHq3nlk9X0adPYh8F2rpbzBuNabs3WPTw9YSntmiXz9HUn+xVGztVy3mhcy322KhuAF27qS6fm9WMcjXMu1rymUMtNW7mV1o2T6JiSHOtQnHNVgCeFWqygsIhpq7I5u2tzbzZyzgERJAVJ10paIWmnpF2SdkvaVRnBueh6d/4mdu7LZ7DPnuacC0XSp/AMcKWZLYl2MK7y5BcW8dKHqzixZUMuOKlFrMNxzlURkTQfZXpCqHmGv7eEFVty+MHF3YmL86Yj51wgkprCLEn/BN4Bcg8Umtlb0QrKRdeLU1Ywctpa7j67E5f29qYj59xXIkkKjYC9wCUlygzwpFANrdySw28nLeeafm346eU9Yh2Oc66KOWJSMLO7KiMQVzkmLt4MwMOXnki8Nxs55w4RydVHbSW9LWmLpExJ/wknz3HV0KTFmfRp05hWjevFOhTnXBUUSUfzXwlGMG1NMOnNu2GZq0aKioxnJizly/U7uLhny1iH45yroiJJCqlm9lczKwgfI4HUKMfljrNPVm7lpamruOLkVtx1dsdYh+Ocq6IiSQpbJd0mKT583AZkRzswd3y9PSedRkl1eP7GU2iYlBDrcJxzVVQkSeFu4EZgM8H0mdeHZa6a2JqTyweLMrn85NYk1omPdTjOuSoskquP1gNXVUIs7jjbsTePH/xzLquy9lBkxtfO7BDrkJxzVVyZSUHSw2b2jKTfE9yXcBAz+14kO5AUD8wCNprZFZKaAf8EOhLMvHajmW0Pl30UuAcoBL5nZh9U7HBcSW/N2ciHy7Lo06Yxv7q6Nz1aNYp1SM65Kq68msKBoS1mHeM+Hgi3deAb6RFgipk9JemR8PWPJfUkmMu5F8GVTpMldfcpOY/e219upFfrRrz73XNiHYpzrpoos0/BzN4Nn+41s1ElHwR3OB9ReD/D5cBfShQPBUaFz0cBV5coH21muWa2BlgJDIj4SNxBVmTuZsHGnVzb328pcc5FLpKO5kcjLCvNC8DDQFGJspZmlgEQ/j0wRGcbYEOJ5dLDsoNIulfSLEmzsrKyIgyj9nnry43Ex4mrTmkd61Ccc9VIeX0KlwFDgDaSXizxViOg4EgblnQFsMXMZksaFEEspY25UFpfxghgBEBaWtph77vgRrV3vtzIud2ak9owMdbhOOeqkfL6FDYR9CdcBcwuUb4b+EEE2z4buErSECAJaCTpdSBTUiszy5DUCtgSLp8OtCuxftswBldBHy7bQsbO/fxkiA9455yrmDKTgpnNA+ZJetPM8iu6YTN7lLCZKawpPGRmt0l6FrgDeCr8+99wlbHAm5J+S9DR3A2YWdH91nZmxv9NWUG7ZvV8WGznXIVFMnR2R0m/AXoS/OIHwMw6H+U+nwLGSLoHWA/cEG5vkaQxwGKC5qn7/cqjisnctZ/b/jKDFVtyeOa6k0mI9ym4nXMVE0lS+Cvwc+B3wPnAXZTe/l8mM5sKTA2fZwMXlrHccGB4RbbtvvLv2ems2JLDs9efzPWn+lVHzrmKi+SnZD0zmwLIzNaZ2RPABdENyx2NCQs307ddE25Ia4fkcyU45youkqSwX1IcsELSdyRdw1eXkboqYsO2vSzYuJPLvB/BOXcMIkkK3weSge8BpwK3E3QQuyrkg0XBjGreueycOxaRDIj3Rfg0h6A/wVVB4xdupkerRnRIqR/rUJxz1Vh5N6+9Syk3jx1gZj5yahWRuWs/s9dt58GLu8c6FOdcNVdeTeG5SovCHZOJYdOR9yc4545VeTevfVSZgbijN37hZrqk1qdby4axDsU5V80dsU9B0hpKH4PoaG9ec8fRtj15zFizjW+d1yXWoTjnaoBIbl5LK/E8ieAO5GbRCcdV1Lj5mygsMi7r401Hzrljd8RLUs0su8Rjo5m9gN+8ViWYGW/OWE/vNo3o1bpxrMNxztUAkTQf9S/xMo6g5uCN11XA9FXZLN28m99c2yfWoTjnaohImo+eL/G8gHBe5ahE4yK2c18+P/r3fNo1q8fQvj6RjnPu+Ijk5rXzKyMQVzF//HAlm3bu4z/fOovkupHkduecO7Ij9ilISpH0oqQ5kmZL+j9JKZURnCvd5p37GTVtLdf0a0P/9k1jHY5zrgaJZOyj0UAWcB1wffj8n9EMypXvHzPXk1dYxA8u8juYnXPHVyRJoZmZ/dLM1oSPXwFNjrSSpCRJMyXNk7RI0pNh+ROSNkqaGz6GlFjnUUkrJS2TNPioj6oGMzPembuRs7qk0K5ZcqzDcc7VMJE0Rn8o6WZgTPj6euC9CNbLBS4wsxxJCcCnksaH7/3OzA4aRkNST+BmoBfBdJyTJXX32dcONi99J+uy93L/+V1jHYpzrgaKpKZwH/AmkBc+RgMPStotaVdZK1kgJ3yZED7KHGAPGAqMNrNcM1sDrAQGRBBfrTJlSSZxgsG9/GY159zxF8nNaw3NLM7M6oSPuLCsoZk1Km9dSfGS5gJbgElmNiN86zuS5kt6TdKBntI2wIYSq6eHZYdu815JsyTNysrKiuQYa5Rpq7I5uW0TGtdLiHUozrkaKKKZ3SVdJem58HFFpBs3s0Iz6wu0BQZI6g28DHQB+gIZfHUfRGnzR5Y25tIIM0szs7TU1NRIQ6kRcnILmLdhB2d39Yu/nHPREcklqU8BDwCLw8cDYVnEzGwHMBW41Mwyw2RRBLzCV01E6UC7Equ1BTZVZD813cw12RQUGWd1aR7rUJxzNVQkNYUhwMVm9pqZvQZcGpaVS1KqpCbh83rARcBSSa1KLHYNsDB8Pha4WVKipE5AN2BmxEdSC0xavIXkuvGc2sHvTXDORUekt8I2AbaFzyMdea0VMEpSPEHyGWNm4yT9XVJfgqahtQQd2ZjZIkljCGojBcD9fuXRVwoKi5i4aDMX9mhJUkJ8rMNxztVQkSSF3wBfSvqQoN3/XODRI61kZvOBfqWU317OOsOB4RHEVOvMXLuN7D15DPHZ1ZxzURRJUpgEfEQwOqqAH5vZ5qhG5Q7zzpcbqV83nkEntoh1KM65GqzMPgVJV0rKAhYQtO1nm9l/PSFUvr15Bbw3P4PLT25FvbredOSci57yOpqHAwPNrBXBuEe/rpyQ3KEmLspkT14h15/a7sgLO+fcMSgvKRSY2VKA8KYzn1gnRqYs3UJqw0RO6+hXHTnnoqu8PoUWkh4s67WZ/TZ6YbkDCouMT1dkcf5JLZBKu7/POeeOn/KSwiscXDs49LWrBAs27mT73nzO61677t52zsVGmUnBzJ6szEBc6T5enoUEA7t5UnDORV9EYx+52PloeRZ92jSmWf26sQ7FOVcLeFKownbuy2fuhh3edOScqzSeFKqwaSu3UlhknOtJwTlXSSIZJfUBSY0UeFXSHEmXVEZwtd2kJZk0SqpD33ZNYh2Kc66WiKSmcLeZ7QIuAVKBu4AKDZ3tKi6voIjJizO5uOcJJMR7hc45Vzki+bY5cHH8EOCvZjaP0ifEccfR9NXZ7NpfwGU+AJ5zrhJFkhRmS5pIkBQ+kNQQKIpuWLVbfmERz09cRtPkBM7p5hPqOOcqTySjpN5DMHXmajPbKymFoAnJRckbn69jfvpOXrq1v8+d4JyrVJHUFM4BGgAnSzoX6EUw6U65JCVJmilpnqRFkp4My5tJmiRpRfi3aYl1HpW0UtIySYOP7pCqv4+WZ9EltT5D+rQ68sLOOXccRVJT+FGJ50kEcyrPBi44wnq5wAVmliMpAfhU0njgWmCKmT0l6RHgEeDHknoCNxMkndbAZEnda9vsa0VFxux12z0hOOdi4ohJwcyuLPlaUjvgmQjWMyAnfJkQPgwYCgwKy0cBU4Efh+WjzSwXWCNpJUECmh7BcdQYK7Ny2LW/wOdhds7FxNFc65gO9I5kQUnxkuYCW4BJ4RDcLc0sAyD8e2AqsTbAhkP206aUbd4raZakWVlZWUcRftU2c00wFXZax2YxjsQ5VxsdsaYg6fcEv/AhSCJ9gXmRbDxs+ukrqQnwtqTykklpl7naYQVmI4ARAGlpaYe9X51l7trPC5NX0K1FAzqmJMc6HOdcLRRJn8KsEs8LgH+Y2WcV2YmZ7ZA0FbgUyJTUyswyJLUiqEVAUDMoObVYW2BTRfZT3T09YSk5ufm8+Y3Tfe4E51xMlNt8JCkeuNjMRoWPNyJNCJJSwxoCkuoBFwFLgbHAHeFidwD/DZ+PBW6WlCipE9CNYG7oWiFrdy7j5mVwU1o7urf0aSucc7FRbk3BzArDL/e6ZpZXwW23AkaFiSUOGGNm4yRNB8ZIugdYD9wQ7muRpDHAYoIayf216cqjN2esJ6+wiK+d1THWoTjnarFImo/WAp9JGgvsOVB4pOk4zWw+0K+U8mzgwjLWGQ4MjyCmGiWvoIjXZ6zjvO6pdEltEOtwnHO1WCRJYVP4iMOn44yK8QszyNqdy53Xd4x1KM65Wi6S+xR8Ws4oGzt3E22a1OM8n3LTORdjZSYFSS+Y2fclvUvpl4ZeFdXIaoncgkKmrcrmhrS2xMX5FUfOudgqr6bwt/Dvc5URSG31xZrt7Msv9Ck3nXNVQnlJ4ZeS3gTeMbM95SznjsFHy7dQNz6OM7ukxDoU55wr9z6FEcAVBOMQ/VPS1ZLqVlJctcbnq7fRr30TkutG0ufvnHPRVWZSMLP/mtkwoAPwFsGNZuslvSbp4soKsCbbtT+fRZt2cnpnryU456qGIw6IZ2b7zOyfZnYNwTzN/YAJUY+sFpi9djtFBmd08sHvnHNVwxGTgqSWkr4r6TPgHWAicGq0A6sNPl+dTUK86Nfeh8l2zlUN5V2S+g1gGHAiQfPRwxUdCM+Vzcz4YNFmBnRqRr26PuWmc65qKK938yzgKWCymRVVUjy1xsKNu1ibvZdvntcl1qE451yxMpOCmd1VmYHUNv+Zk06dOHFp7xNiHYpzzhU7mpnX3DFaunkXr3++jmv6taFJsl/l65yrOjwpVLKiIuMnby2gYVIdHh3SI9bhOOfcQSJKCpLOkXRX+Dw1nATHHYU3Z65nzvodPHZ5T5rV91qCc65qieSS1J8DPwYeDYsSgNcjWK+dpA8lLZG0SNIDYfkTkjZKmhs+hpRY51FJKyUtkzT46A6p6srJLeB3k5ZzeqdmXNu/TazDcc65w0QytsI1BDeszQEws02SIplXoQD4oZnNCZefLWlS+N7vzOyggfYk9QRuBnoBrYHJkrrXpNnXRn62huw9efzlspN8DmbnXJUUSfNRnpkZ4fDZkupHsmEzyzCzA4lkN7AEKO/n8VBgtJnlmtkaYCUwIJJ9VRfj5mdwRudmfrOac67KiiQpjJH0Z6BJeEPbZOCViuxEUkeC2saMsOg7kuaH4ygd+IZsA2wosVo6pSQRSfdKmiVpVlZWVkXCiKmd+/JZlrmbMzs3j3UozjlXpkjGPnoO+DfwH4K7m39mZr+PdAeSGoTrft/MdgEvA12AvkAG8PyBRUvbfSnxjDCzNDNLS02tPnMQzF63DTM4rZPXEpxzVVdE4zWb2SRg0hEXPISkBIKE8IaZvRVuK7PE+68A48KX6UC7Equ3JZgbukaYuWZ7MM5RO08KzrmqK5Krj3ZL2nXIY4OktyV1Lmc9Aa8CS8zstyXKW5VY7BpgYfh8LHCzpMTwktduwMyjOaiqpqjImLhoM/3aNfVxjpxzVVokNYXfEvxif5Ogiedm4ARgGfAaMKiM9c4GbgcWSJoblv0EGCapL0HT0FrgPgAzWyRpDLCY4Mql+2vKlUcfLc9i9dY9PHBRt1iH4pxz5VJwYVE5C0gzzOz0Q8o+N7MzJM0zs1OiGmE50tLSbNasWbHafUTMjBv+NJ0N2/fyycMXULeO30TunIstSbPNLK209yL5hiqSdKOkuPBxY4n3ys8ojre/3Misddv5/kXdPSE456q8SL6lbiVoBtoCZIbPb5NUD/hOFGOr9j5fnc1P315Iv/ZNuCmt3ZFXcM65GDtin4KZrQauLOPtT49vODXH/vxCHhj9Ja2bJDHi9jTi4vwOZudc1XfEpCApCbiHYPiJpAPlZnZ3FOOq9kZNW0vmrlxevLcfqQ0TYx2Oc85FJJLmo78TXG00GPiI4P6B3dEMqiYY/cUGzuycwumdU2IdinPORSySpNDVzB4H9pjZKOByoE90w6reMnbuY83WPVzYo0WsQ3HOuQqJJCnkh393SOoNNAY6Ri2iGmD6qmwAzvBagnOumonk5rUR4aB1jxHcddwAeDyqUVVzn6/OpnG9BHq2ahTrUJxzrkLKTQqS4oBdZrYd+Bgoc1gLFygqMj5evpUzO6f4FUfOuWqn3OYjMyvC70WokLnpO9i8az+De7eMdSjOOVdhkfQpTJL0UDi9ZrMDj6hHVk1NWLiZhHhxwUmeFJxz1U8kfQoH7ke4v0SZ4U1JpZq8OJMzuzSncb2EWIfinHMVFskdzZ0qI5CaYPPO/azeuodhA9rHOhTnnDsqkcynkCzpMUkjwtfdJF0R/dCqn89XB5eintnFL0V1zlVPkfQp/BXIA84KX6cDv4paRNXY9FXZNEqqQw+/FNU5V01FkhS6mNkzhDexmdk+Sp9P+SBhx/SHkpZIWiTpgbC8maRJklaEf5uWWOdRSSslLZM0+CiPKSbMjE9XbuX0zinE+6WozrlqKpKkkBcOk20AkroAuRGsVwD80Mx6AGcA90vqCTwCTDGzbsCU8DXhezcTDLx3KfCSpGozd+WCjTvZuGMfl/T0q46cc9VXJEnhCWAC0E7SGwRf5A8faSUzyzCzOeHz3cASoA0wFBgVLjYKuDp8PhQYbWa5ZrYGWAkMiPhIYuz9BZupEycu9qTgnKvGIrn6aKKk2QS/9gU8YGZbK7ITSR2BfsAMoKWZZYTbzpB0YNS4NsDnJVZLD8uqPDNj/MIMzuranCbJdWMdjnPOHbVIrj4aC1wCTDWzcUeREBoA/wG+b2a7ylu0lLLDpvuUdK+kWZJmZWVlVSSUqFmcsYt12XsZ0vuEWIfinHPHJJLmo+eBgcBiSf+SdH048c4RSUogSAhvmNlbYXGmpFbh+60IpvmEoGZQcs7KtsCmQ7dpZiPMLM3M0lJTUyMJI+rGL9hMfJy4pJcnBedc9XbEpGBmH5nZtwnuYB4B3MhXX+RlkiTgVWCJmf22xFtjgTvC53cA/y1RfrOkREmdgG7AzEgPJFbMjPcXZHBG52Y0q+9NR8656i2SYS4Irz66ErgJ6M9XHcXlORu4HVggaW5Y9hPgKWCMpHuA9cANAGa2SNIYYDHBlUv3m1lh5IcSG8szc1i9dQ93n+M3fjvnqr9I5mj+J3A6wRVIfyToWyg60npm9ill389wYRnrDAeGH2nbVcn7CzKQYLA3HTnnaoBIagp/BW458Ktd0tmSbjGz+4+wXq0wfmEGAzo2I7VhYqxDcc65YxZJn8IEoI+kpyWtJRjiYmm0A6sOVm7JYXlmDkP6tIp1KM45d1yUWVOQ1J3gDuNhQDbwT0Bmdn4lxVblTViYAcClfimqc66GKK/5aCnwCXClma0EkPSDSomqmnh/wWZO7dCUlo0iukLXOeeqvPKaj64DNgMfSnpF0oVEMBBebbF26x4WZ+ziMq8lOOdqkDKTgpm9bWY3AScBU4EfAC0lvSzpkkqKr8oav3AzAJd5f4JzrgaJpKN5j5m9YWZXENxlPJdwZNPabPzCDE5p14Q2TerFOhTnnDtuIhnmopiZbTOzP5vZBdEKqDrYsG0v89N3+lhHzrkap0JJwQUmHGg66u1NR865msWTwlF4b0EGvds0on1KcqxDcc6548qTQgV9uX47czfs4Oq+1WKqB+ecqxBPChX0+/+tpGlyAsMGtI91KM45d9x5UqiAhRt38r+lW7jnnE7UT4xogFnnnKtWPClUwItTVtAoqQ5fO6tjrENxzrmo8KQQoazduUxaksltZ3SgUVJCrMNxzrmo8KQQoSlLMjGDK05uHetQnHMuaqKWFCS9JmmLpIUlyp6QtFHS3PAxpMR7j0paKWmZpMHRiutoTVqcSZsm9ejRqmGsQ3HOuaiJZk1hJHBpKeW/M7O+4eN9AEk9CYbp7hWu85Kk+CjGViHb9+TxycqtXNyzJcHU0845VzNFLSmY2cfAtggXHwqMNrNcM1sDrAQGRCu2inpjxjryCoq45XS/DNU5V7PFok/hO5Lmh81LTcOyNsCGEsukh2WHkXSvpFmSZmVlZUU7VnILChk1fR3ndU+le0tvOnLO1WyVnRReBroAfYEM4PmwvLQ2GSttA2Y2wszSzCwtNTU1KkGWNHbuJrJ25/KNgZ2jvi/nnIu1Sk0KZpZpZoVmVgS8wldNROlAuxKLtgU2VWZspTEzXv10DSed0JCzu6bEOhznnIu6Sk0KkkoOK3oNcODKpLHAzZISJXUCugEzKzO20sxZv52lm3dz19kdvYPZOVcrRG2sBkn/AAYBzSWlAz8HBknqS9A0tBa4D8DMFkkaAywGCoD7zawwWrFFaswX6STXjfd7E5xztUbUkoKZDSul+NVylh8ODI9WPBWVnZPLuPmbuOLkVj7OkXOu1vA7mssw/L0l5BUWeQezc65W8aRQik9XbOWtLzfyzfO60M0vQ3XO1SKeFA6RW1DIY+8soFPz+tx/ftdYh+Occ5XKG8sP8eaM9azN3svf7h5AUkKVGWnDOecqhdcUStiXV8gfP1zJWV1SOLd79G+Mc865qsaTQglj521ka04e37uwW6xDcc65mPCkEDIz/v75Orq3bMDpnZrFOhznnIsJTwqh2eu2s3DjLm4/o4Pfveycq7U8KYRenrqKpskJXHdq21iH4pxzMeNJAVi4cSdTlm7hjrM6klzXL8hyztVetT4pmBnD31tC0+QE7j6nU6zDcc65mKr1SWHaqmymr87mexd2o1FSQqzDcc65mKr1SeEP/1tJi4aJDBvgU20651ytTgofLtvC9NXZ3HtuZ7972TnnqMVJYdf+fB75z3y6t2zAbWd0iHU4zjlXJUQtKUh6TdIWSQtLlDWTNEnSivBv0xLvPSpppaRlkgZHK64DXvpwFVt25/LcDad4LcE550LRrCmMBC49pOwRYIqZdQOmhK+R1BO4GegVrvOSpKh9U6dv38trn63hmn5tOLltk2jtxjnnqp2oJQUz+xjYdkjxUGBU+HwUcHWJ8tFmlmtma4CVwIBoxZZbUMQZnVN46JITo7UL55yrlir7Tq2WZpYBYGYZklqE5W2Az0sslx6WHUbSvcC9AO3bH90VQ11SG/C3u6OWc5xzrtqqKh3NpQ02ZKUtaGYjzCzNzNJSU314a+ecO54qOylkSmoFEP7dEpanA+1KLNcW2FTJsTnnXK1X2UlhLHBH+PwO4L8lym+WlCipE9ANmFnJsTnnXK0XtT4FSf8ABgHNJaUDPweeAsZIugdYD9wAYGaLJI0BFgMFwP1mVhit2JxzzpUuaknBzIaV8daFZSw/HBgerXicc84dWVXpaHbOOVcFeFJwzjlXzJOCc865YjIr9XaAakFSFrDuGDbRHNh6nMKJpZpyHODHUlX5sVRNR3ssHcys1Bu9qnVSOFaSZplZWqzjOFY15TjAj6Wq8mOpmqJxLN585JxzrpgnBeecc8Vqe1IYEesAjpOachzgx1JV+bFUTcf9WGp1n4JzzrmD1faagnPOuRI8KTjnnCtWK5OCpEvDuaBXSnok1vFUlKS1khZImitpVlhW5vzXVUlVn7u7Iso4lickbQzPzVxJQ0q8VyWPRVI7SR9KWiJpkaQHwvJqd17KOZbqeF6SJM2UNC88lifD8uieFzOrVQ8gHlgFdAbqAvOAnrGOq4LHsBZofkjZM8Aj4fNHgKdjHWcZsZ8L9AcWHil2oGd4fhKBTuF5i4/1MRzhWJ4AHipl2Sp7LEAroH/4vCGwPIy32p2Xco6lOp4XAQ3C5wnADOCMaJ+X2lhTGACsNLPVZpYHjCaYI7q6K2v+6yrFqvDc3RVVxrGUpcoei5llmNmc8PluYAnBdLjV7ryUcyxlqcrHYmaWE75MCB9GlM9LbUwKbYANJV6XOR90FWbAREmzwzmr4ZD5r4EWZa5d9ZQVe3U9V9+RND9sXjpQta8WxyKpI9CP4FdptT4vhxwLVMPzIile0lyCWSonmVnUz0ttTAoRzwddhZ1tZv2By4D7JZ0b64CipDqeq5eBLkBfIAN4Piyv8sciqQHwH+D7ZrarvEVLKavqx1Itz4uZFZpZX4IpigdI6l3O4sflWGpjUqj280Gb2abw7xbgbYIqYlnzX1cHNWbubjPLDP8jFwGv8FX1vUofi6QEgi/RN8zsrbC4Wp6X0o6lup6XA8xsBzAVuJQon5famBS+ALpJ6iSpLnAzwRzR1YKk+pIaHngOXAIspOz5r6uDGjN394H/rKFrCM4NVOFjkSTgVWCJmf22xFvV7ryUdSzV9LykSmoSPq8HXAQsJdrnJdY97DHq1R9CcFXCKuCnsY6ngrF3JrjCYB6w6ED8QAowBVgR/m0W61jLiP8fBNX3fIJfNveUFzvw0/A8LQMui3X8ERzL34EFwPzwP2mrqn4swDkEzQzzgbnhY0h1PC/lHEt1PC8nA1+GMS8EfhaWR/W8+DAXzjnnitXG5iPnnHNl8KTgnHOumCcF55xzxTwpOOecK+ZJwTnnXDFPCq7WkJRSYpTMzSVGzcyR9FIU9neipKnhPpZIGhGW9y05SqdzVUmdWAfgXGUxs2yCYQ6Q9ASQY2bPRXGXLwK/M7P/hvvsE5b3BdKA96O4b+eOitcUXK0naZCkceHzJySNkjRRwbwV10p6RsH8FRPCIRSQdKqkj8JBCT845I7ZA1oR3NQGgJktCO+i/wVwU1iDuCm8S/01SV9I+lLS0HAfd0r6b7jfZZJ+HpbXl/ReOM7+Qkk3RfszcrWH1xScO1wX4HyC8emnA9eZ2cOS3gYul/Qe8HtgqJllhV/Kw4G7D9nO74D/SZoGTAT+amY7JP0MSDOz7wBI+jXwPzO7OxzWYKakyeE2BgC9gb3AF+G+OwCbzOzycP3GUfocXC3kNQXnDjfezPIJhkWIByaE5QuAjsCJBF/Uk8JhjR8jGHzsIGb2V6AH8C9gEPC5pMRS9ncJ8Ei4ralAEtA+fG+SmWWb2T7gLYJhHBYAF0l6WtJAM9t5jMfrXDGvKTh3uFwAMyuSlG9fjQVTRPB/RsAiMzvzSBuyYETb14DXFEzbWdrQxyKojSw7qFA6ncOHPjYzWy7pVIIxfX4jaaKZ/aICx+dcmbym4FzFLQNSJZ0JwVDNknodupCCucAP9EGcQDCQ2UZgN8FUkQd8AHw3HOETSf1KvHexgjl56xHMsPWZpNbAXjN7HXiOYEpQ544LTwrOVZAF07heDzwtaR7BSJxnlbLoJcDCcJkPgB+Z2WbgQ6DngY5m4JcEUy3OD2sTvyyxjU8JRvicC/zHzGYBfQj6HeYSjIr5q+N+kK7W8lFSnauiJN1JiQ5p5yqD1xScc84V85qCc865Yl5TcM45V8yTgnPOuWKeFJxzzhXzpOCcc66YJwXnnHPF/h917ygVEnNFQgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "simulationWithoutDrug(100, 1000, 0.1, 0.05, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}