EnonceTP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modèles à base d'équations différentielles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The classical prey predator model\n",
    "\n",
    "def model(C,t):\n",
    "    prey=C[0]\n",
    "    predator=C[1]\n",
    "    dprey = prey*((2/3)-(4/3)*predator)\n",
    "    dpredator = predator*(1*prey-1)\n",
    "    return [dprey,dpredator]\n",
    "\n",
    "t = np.linspace(0,100,100000)\n",
    "ds = odeint(model,[1,2],t)\n",
    "\n",
    "X=ds[:,0]\n",
    "Y=ds[:,1]\n",
    "\n",
    "%matplotlib inline\n",
    "plt.figure(1)\n",
    "plt.plot(t,X,'b-',t,Y,'r')\n",
    "plt.xlabel('Temps')\n",
    "plt.ylabel('Abondance')\n",
    "plt.legend(['proies','prédateurs'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question 1 : tracer le plan de phase courbe (X(t),Y(t)) et déterminer l'équilibre du système"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question 2 : tracer les courbes des espèces du modèle BIOMD0000000006 - Tyson1991 - Cell Cycle 2 var\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exemple tiré de Daun, Silvia et al. “Equation-based models of dynamic biological systems.” Journal of critical care vol. 23,4 (2008): 585-94. doi:10.1016/j.jcrc.2008.02.003\n",
    "\n",
    "\n",
    "\n",
    "The viral kinetics of influenza A during infection within an individual can be modeled using a simple 3-dimensional ODE System, which incorporates viral growth, target cell limitation and the interferon response:\n",
    "\n",
    "$$\n",
    "\\frac{dT(t)}{dt}=−βT(t)V(t),\n",
    "$$\n",
    "$$\n",
    "\\frac{dI(t)}{dt}=βT(t)V(t)−δI(t),\n",
    "$$\n",
    "$$\n",
    "\\frac{dV(t)}{dt}=pI(t)−cV(t),\n",
    "$$\n",
    "\n",
    "where T(t) is the number of uninfected target cells, I(t) is the number of infected cells, and V(t) is the infectious viral titer (or total number of free viral particles or concentration of viral particles per unit volume of tissue) at time t. It is assumed in this model that infection is initiated by introduction of virus into the upper respiratory tract at a concentration equivalent to V0 viral particles per unit volume of tissue. Susceptible cells become infected by virus at rate βT(t)V(t), where β is the rate constant quantifying the probability that the interaction between a healthy cell and a viral particle results in infection of the cell. Infected cells, I(t), increase viral titers at an average rate of p per cell and die at a rate of δ per unit time, where 1/δ is the average life span of an infected cell. Free virus is cleared at a rate of c per unit time. The effects of immune responses are not explicitly described in this simple model, but they are implicitly included in the death rate of infected cells (δ) and the clearance rate of virus (c). Such simple systems may still admit complex solutions and several equilibria, depending on the choice of parameters.\n",
    "\n",
    "\n",
    "*Question* : Modéliser ce système.\n",
    "\n",
    "*Question* : On suppose que β=0.1 , δ=0.2 , p=0.5 ,c=1 , V0=1. Simuler ce système, puis tracer l'espace des phases pour différentes valeurs de T0 et I0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "attachments": {
    "image-2.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyAAAADhCAYAAAAwLNHkAAAAAXNSR0IArs4c6QAAAGhlWElmTU0AKgAAAAgABAEGAAMAAAABAAIAAAESAAMAAAABAAEAAAEoAAMAAAABAAIAAIdpAAQAAAABAAAAPgAAAAAAA6ABAAMAAAABAAEAAKACAAQAAAABAAADIKADAAQAAAABAAAA4QAAAAAQivvBAAACC2lUWHRYTUw6Y29tLmFkb2JlLnhtcAAAAAAAPHg6eG1wbWV0YSB4bWxuczp4PSJhZG9iZTpuczptZXRhLyIgeDp4bXB0az0iWE1QIENvcmUgNS40LjAiPgogICA8cmRmOlJERiB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMiPgogICAgICA8cmRmOkRlc2NyaXB0aW9uIHJkZjphYm91dD0iIgogICAgICAgICAgICB4bWxuczp0aWZmPSJodHRwOi8vbnMuYWRvYmUuY29tL3RpZmYvMS4wLyI+CiAgICAgICAgIDx0aWZmOk9yaWVudGF0aW9uPjE8L3RpZmY6T3JpZW50YXRpb24+CiAgICAgICAgIDx0aWZmOlBob3RvbWV0cmljSW50ZXJwcmV0YXRpb24+MjwvdGlmZjpQaG90b21ldHJpY0ludGVycHJldGF0aW9uPgogICAgICAgICA8dGlmZjpSZXNvbHV0aW9uVW5pdD4yPC90aWZmOlJlc29sdXRpb25Vbml0PgogICAgICAgICA8dGlmZjpDb21wcmVzc2lvbj4xPC90aWZmOkNvbXByZXNzaW9uPgogICAgICA8L3JkZjpEZXNjcmlwdGlvbj4KICAgPC9yZGY6UkRGPgo8L3g6eG1wbWV0YT4Kpl32MAAAQABJREFUeAHsnQe8XUW5t0du6EUh9BIIQgKREnoNBghFcxGkKiJdQcVGFRGFK/ciqCAIeulIEWkC0iGUECQUwdB7QsilE3oN935853nxvzNZ7L3PPufsnbPLf36/vWfW9HnWrFnzrmmf+bhLJSsTMAETMAETMAETqEJA3YVXX3013XHHHWnKlCnpySefTI899liaOHFiGj9+fFpwwQWrxGAnEzABE/iEwACDMAETMAETMAETMIHuCHzmM58JLx9++GEaOHBgeuONN9KJJ54Ydttss42Fj+4A2t0ETKBE4DMeASmxsMEETMAETMAETCAjwKiHBA9Z53YrrLBCjIKcdtppaY899kj/7//9v/Rv//Zv8mrdBEzABMoS8AhIWSy2NAETMAETMAETKAofEEHIQAi5++670xNPPJFmmWWWtMEGG4Sggv9cQDFBEzABEyhHYJZylrYzARMwARMwARMwgXIEEDAGDBiQbr/99hA21lprrTR48OCS4FFOaCkXj+1MwAQ6l4AFkM699y65CZiACZiACdREAKEDJeGDUZAbbrgh7EaNGpVmm222cAsL/5mACZhANwQ8BasbQHY2ARMwARMwgU4nkI9qIHywA9bYsWNj+tXIkSND+Pi///u/uMYvflCY81+nc3T5TcAEPiFgAcQ1wQRMwARMwARMoCYCCBYsMr/zzjvT//7v/6aFF144rbbaaiFkzDrrrKU4WBdiZQImYAKVCFgAqUTG9iZgAiZgAiZgAmUJ3HzzzWG/6aabpgUWWCBNnTo11oQw2jHXXHPFSMg777yTpk2bluaff/608cYbh13ZyGxpAibQcQT8iaLjbrkLbAImYAImYAK9I8DoB4LFjTfeGBEggLAuhNEQBI9bbrklbbbZZgn7E044IT3//PNp7rnnDuEjX0fSu9QdygRMoF0IeASkXe6ky2ECJmACJmACDSKA8MD0K6ZWPfTQQ2nSpElhHjFiRKS40EILheDBCemchn7ppZemjTbaqJQbwjM6gpJecrTBBEyg4whYAOm4W+4Cm4AJmIAJmEDPCCA0SIhg8Tlq3XXXTcsss0xpkfmhhx6azj///HTXXXelZZddNvwjtBDWa0J6xtu+TaDdCVgAafc77PKZgAmYgAmYQB0IMP0KgULb77L7FdvvMv3qgAMOSBMmTEgPP/xwmnfeeRM7YuEfwcMjHnWA7yhMoM0IWABpsxvq4piACZiACZhAvQlo+tXkyZNj+12Eiq233jq9/fbbac8994zRDtaFIJBoO16NmNQ7L47PBEyg9Ql4EXrr30OXwARMwARMwAQaSgABBDV+/PgQMIYMGRLCBovNL7nkkrTccsvNIHwgoPBDCLEyARMwgSIBCyBFIr42ARMwARMwARMoEUCI0BoOdrlCMe1q9913jy14mWr1l7/8Jb333nufmnJVbvqVhZISWhtMoGMJfKarIfDniY69/S64CZiACZiACVQnoG7Cu+++m4YOHRpb6zICcvHFF8doyOqrrx6jHY888khaYYUVSrtllYuVuMoJJeX82s4ETKB9CXgEpH3vrUtmAiZgAiZgAn0mIKHhwQcfDOGDEQ+Ej1VWWSUtssgiab755oupVixAR+GfKVusBdE1doonLGv4w39RlbMr+vG1CZhA8xOwANL898g5NAETMAETMIF+I4AgQcef6VeMXqy55ppp5ZVXDiEDAUQjIJwBIgGBKVsIKlwTptqoh8KogLouFwY7ucu/dRMwgdYjYAGk9e6Zc2wCJmACJmACM4UAnf1ZZ501Ov1jxowJfcsttwyB4qOPPoo1HxtuuGHY33777WnatGlhx+5Yv/71r0t5JJ6iQCFBIhcq5E9uiiC/LsYjP9ZNwARah4AFkNa5V86pCZiACZiACcxUAur4T5kyJSFgoDj/A6UzPtZZZ524vvfee9NLL70Ugsbee++dhg0bFvZMxyonNGCHm4QOPMsfOvb66Toi9J8JmEDLE/A5IC1/C10AEzABEzABE6g/AQkf6EyvYsRj4MCBafjw4ZEYU6xQq622WggOCAkcSIi/zTbbLI0ePfpTC9KJC3/5QYXEIfvcjD+U/OoaoUW7coUH/5mACbQcAY+AtNwtc4ZNwARMwARMoPEE1OFHZ/RjrrnmStttt1363Oc+FwIDQgCCw+KLL5522GGH2JqXM0G22GKL9L3vfS+ED8Wh3HJNGIQXdtVi4fpjjz0WQgn2KPl555130htvvBF+X3311TRx4sRPCTSK17oJmEBrEfA2vK11v5xbEzABEzABE5jpBKZOnRqCwxxzzJHmmWeeUvoauXj55Zdjkfraa6+dlllmmZKAUvLYZZBf7M4888z05ptvJqZvPf300zHCctJJJ8XIxnXXXZeuvPLKdNVVV6VTTz01hJArrrgiXXjhhWnHHXdM55xzTqxLKQo3eVo2m4AJNDcBj4A09/1x7pqUAC9SbTOJmZ9V7wjk7GSW3rsYHcoEZh6BSnW1kv3My1l9U2Lq1YILLjiD8EEKEgIWXnjhtNNOO6XBgweHXbkpUjCh3dxjjz3SpEmT0g9+8IO0/vrrp5tuuimdcsop6cMPP4ywrDFh/QjrTh566KG0xBJLpPPOOy9+CCHXX399+OMwxKJqN+7F8tVybQa1ULKf/ibgNSD9fQecfssQ0GJJpg7w0tWLl8ZeZgqja+ktU8B+yqjYiZf0fsqOkzWBHhGg/lJnUarLRXM4tvkfDHIOYqHnWes4Dj/88PToo4+ms846K4QR1otwnggjHkzxwt+cc86ZWNDO+SIINUsttVTQW2uttWKE5LnnnotrpZGjLWeXu7eiWQzzvOt9lNuV85e7t4s5r2e1lqkd2VQqUyX7WlnVyx/Pou6V4sSOussHCgsgomLdBCoQ0AOUf9H74IMPElMSnn322fg6N2jQoNJDpYeuHV+EFRDVxTrnpgaqLhE7kqYn0CwvzN6AauW896a8lcLw/ObPsPzJjg83zzzzTDrqqKPSjTfeGO0lzzlb/O6///7yHus9WPtx/vnnp6OPPjqED0ZGZpttthgRIcyKK64Y/om7nVVet3IzZdb6GzHI3XNzO/Nx2ZqbAPUw7zcpt6q7FkBExLoJVCCgBv7OO+9Md999d8xzfuKJJ9KTTz4Ziy75moeSPzf+FUB2Y51zK9dodRPczi1MQM9OKxahlfPeCN7leCA0IIBwejrP9qqrrho6fvPnXqMk999/f+yktf3220cWNerMmpCFFloosc4kD9eIcjRDnGIpXWX+n//5n/T3v/89RotghrtY4kcdvGYog/PQuQQq1Unq54ABAzwC0rlVwyWvlQANOr9FF100FkCy4wsLIbHjQK6hQ4fOEFXxZTGDoy8+RUB80WGn0Q+urRpPQNxJKTc3PuXpKfRHuuXSxA6lZ3h6Dsub5J8XKmFqDVc+tva0zTk/8MADiUXs8847b3y8gRfCxYsvvpjmn3/+6JTg/9prr03LLbdcYmSZDjZ833rrrXT88ceHEMMULdZ/0InpJCVBbsKECelrX/ta6YBI1cO8Dsquk/i4rI0jQN1SncrNtaZIGJ5Znnet3eqsp7dWUvZnAhkBNerLdO3sgmIeMg8RL0a2m0QVX4Y8qISz6p6A+OITbnQq0Ol0WJmACbQuAbWDdJxRSy65ZHrvvffSPffck0aMGBF2nC/y1FNPxVoP2lX8sj7kwAMPjDYUO9S+++6bvvOd78Q2wPjpBOFD/ALAv/6wYzoa5X/++edjnQxtpfyi85MdwbjOVa3vJsWpsHk8tcSR+ycOwhTtFHct8cmvdOVPcSoO2ctfLXoxDsIoHrlhl6eRX2NGya/K2p3/T0JN/1ea022m5wO7cu7yW81NfsrpyrPclGeu8zhzs9zQc/9cS+Ff9RAzzy3Xf/vb30KAtgAiUtZNoAoBHh6EDOYr33XXXSF8YLfRRhtFKB4qFHY8jJUeyPDUwX/iAwI1RrxE2dUGvtjhZ2Yq7pXSlFn6zMqH6ovyMbPSrSWd7ljk7nk5cvvu0snDyW8xfDU/ud/crLh6q1eLS260CewOtfXWW5ee/7ye9zbtdggHIymYfPOb30wnnnhi+uIXv5i++tWvxkgIhxqyG5YU60RoE1hnxxkgnBXCDlmMQP/2t78Nb2pvFaZd9ZyfyogdH79oLxFE+JVTroPlqNhuZhNQPZRO+gjP1GELIDP7bji9liRAo8+XOB4iFlCiDxkyJLaKpEB6UUhvyULO5EzDEMUhZHvuuWd0RuhYiKEaLPnL7Qmna8zyg7mningUXmbp3cVVzIOulffuwufpym9uJ7PcGqErz3ncSlduXBfNuZ/cTDzyL3vFrTh0LXfZ5+GwK7qXC5eHwV3hpCuM0tA1uuLHnLvLXnFIlz+lSZvAZhQrrLBCCCCKQzr+rVK0nXxcYGvdRx55JHHOB9xGjRoV6zlgRIea53/MmDExTYt1HmPHjg18W221VYw8i7uZfkIApih9zPnE9pP/etbB/HlQGrLjuphW7iZ33buim+IrxiH7ol4pfNEf18Qp/3n8yovCyI+u5Te3L9op7lxXeIWTG/blwpfzn9spvMIW3XQtvViu3B5zHo/yKD+VdOVB4fNweXwKX86d+olfnnGNcuK/rQQQFbwcFMFplE7a1dLtzr1R+XK89SHA/ePlyMFZV199dUTK+g/mM6vx1z2WjqfcnOekkn3upx3N+TMiM1/wYMs0DBaYiqfca+EAT5TC5NdiLR1/cq8UpuiHaymloes8Ltmh46+YpsLWap/HIzN6Hj6/ln2uKxx25ZTyVHRTHEX7/LpSnLkfxa/4cj331xtzkbHiyNOQH9xys/xWs5cflUHXxB9f8Lq+5LFT07nnnhtOSlf+rE8nwDOOWnrppdM+++xTcig+75ykzmjSpptuWvIDV3ViSpYdbFB9lC62jUKidPL4y9nJvZyb7KTLb0/1noYv579oV7xWnsrZ53YySy8XrppbOf9Fu2J43MvZ9dS+UhxKP9dzv7k59yNz0Z1r2cms67YSQFQogUAv90LADpX7L+cvPNXwV0tYpVWL3xqStJeZTICXH5I7iyhfeOGFSH3jjTcOXS9G3WPpOObmvN7l9hFJB//BJWcD556qIs/8WmbpxJ2blVbRrngtf0W9mr/crV5m0s/jyq9lX9RzP5hrUYqjmt9a/Ci8/BZ1ufdWV3x5eNkVdfzILvdfzb7oT34VD1/1+LltL0fq03a0l3rm6TjzU/vKtCtGPZgjjuKMEPmBt5h/OlbbmIAJtBqBthJAysEv12AV7fTikK54iteyz3X89EQV0+5JWPvtPwK6z7feemtk4rOf/WzaYIMNwqzFkPhRndELNc9xfu/lL3fvVLO4SYeDWIkTuuxl1nU4/OtPbgqfu2HO45OfYpjidbk4ZEccRf/5NeZKfsqlXwxLOtXCK3785WF1nYfFTmnKHb2Rqlz6eT5lVh7y/Mmumq7yS8evzEVdbujKl3TsUHn6hJfK7bHDTe646Sf/1isTKPe1Hpavv/56OvPMM9O0adNidImzQDiYsJz/yrHbxQRMoFUItL0AQsNWfHno5shN7tLlXryWfa7nfhSf3PNrmdFReTj5t968BBAy+MrJ+g/UJptsEgtPdV+x457qvha/4uf+5BfdajoB8dOIEi7iKT23mx5yuin3N912uknu0nHJzeWup4f+xNSd/9xdZumKK7/uqTnPY6Wwlfwo/dw9t2uEOc9jMd2iW0/TV3jpefyyk5675ebcPU+/kr3C5u5q1/PwNtdOAJYPPvhg7Oh09tlnx+JzzrpgbY1UsQ2VvXUTMIHWJND2Akj+kshvUbXGLH+ZKHxul8cjM/6KfvNrmdG7i0txWm8OAqorkydPTmwZyT1kfjI6L8nzzjsvvfrqq7FzC1/r+JKHwLLEEkukDTfcMO28887hV/EU9eYoZfPkQs9K8+TIOTGB7gnoue7ep30UCdBusqOgdhUsx9LtQpGar02gtQm0vQDC7ckbMxYO0tih06Bpqgx+NNRbrqErZ6dbT9hcKV7pRXelX/xKnsdhc/MQ0GLT22+/PeoL95NtJFELL7xwbCHJQVocnIUbO7awreSKK64YoyQqieoQOv50LfdO1cVBz4nZdGpNaL1yq66qDrdeCZonx7DkfYziHYlZ7+TmyaVzYgImUC8CbSeA0IjpByTMeUdfZs3b5xo/xYbuww8/DHuBxo8UZr1wpOdu7IxEfHm8eRj85unhph9uxKkf11b9S0B15Zprron7tNpqq6Xll18+XpC4sYsT20pyz04//fTYUlY55r4WVbEuFN076TrnU3yWOomDy9qaBKizarvz+utnvHf3k/eimObvyN7F5lAmYALNTKDlBRB9MQEyDZd+OXR1ct5///309ttvp5dffjk9/vjjsdiN6TJPP/10vETY4Qh3djmaMmVKKQo1iCWLf6VV7iWDHVNvFl988cjLKqusEvNa55577rTSSitFOksttVQaOHBgWmyxxcKtXJ6VFusO5C5dbtZnHgH2rL/lllsiQaZfzT777KXEjzzyyHTsscfG/vXsjKXRtUr3C/tydacUYQcZYFFU5eyKfvpybfZ9oeew5Qiozqpu6bqcX9uVJyBmYljel21NwATahcBMFUB60rBU8os9QgeNFV9I9JVE/nFDgHjiiSfSc889lx5++OHQH3rooXT//feXvlZxA4mDEZA11lgjzIMHD06cysrJtggLdDKJtzuFH/KhfCHQ0GFlN4+777473Fg/oO1bc6Fp0KBBCYGEdJmyQx6WXXbZ2Ct9zjnnjLUEefqEVZ40mpO7yyweuu6J3pewPUmnFfxqyt59992XXnnllcjy5ptvHjpbRu6///6xBmTChAkxKoLAyH3Ry7RSGbtzrxSuHe1Vn3O9kXwaGXc73h+XqXsCqruuW92z6s6HGXZH6JOZHd37sg8TaE4Cai9nqgCihoXEZa6EB3f5oxOIopOPvTre2DP3ntGMe++9N40fPz79/e9/jxEO/OOPDj0HH/Flet99941Rh89//vNpkUUWSfPMMw/eZviaHRa9+FNeCZqb82s6rBKQWLT87LPPxs4fLGSmg3vqqafGvueEoZzrr79+/DgR9gtf+EJJKIEDaaAkkOBfP+wxoyrlJRwr/ClsBeeOs4YH2++izzfffLH9LqNou+66a+zW8swzz6QFFlggdsnSdK2Og+QCm4AJmIAJzFQC6gcU3/MzNRNOzAR6SWCmCiB6SIodXNlLpyz68owZQYKONuE4LZlRhTFjxqQrr7wyRhoIxxQndtD40Y9+lJijz2gCU5yYn19O5WlhliKd3iryp3zmcWiUhrUhKAQgfuuss07aYYcdwo48cOgSIycSTNhxiQOZfv3rX4cfODD9Z9SoUbG7EqMmjJLkCm7kQ2mKdV7e3L/NlQnADI7clxtuuCGEuW984xux29WXvvSlqIunnHLKDMKHOVfmaRcTMAETMAETMAETgEDDBZBaOmTquKvDx7VGOd544400bty4+AKNwPHkk0/GnWPa1He+8534Go2wwZoLvj4rPel4llCBnRQdS+zRSU9uSlf+eqITR56HPCxu/HKBAHfZkQ+EJYSmRRddNDHqsffee0cUr732WpSbfdIRvA4//PDEInn8jx49OjEliHMphgwZUuJGQAkjmIkfpfTjwn9VCVA/qA8IvYywwZBpWCuvvHJiNAuWLEznPqnemG9VpHY0ARMwARPoIwG9Z4p6H6N1cBOYqQQaLoDoAVFHu1g6CQfqIOOPDjdbnv7lL39JF1xwQXT0EDLU2WbkgCkvuSIcCp2OtwQB7Igbe6WBHYq8YY+ufGIvO8w9UYojj5fwxIdSPorphWPXn1jgX3mgYzv//PPHaAnlprOL8MEaF4QRdl/63ve+F2FXXXXVtNtuuwUndmkirOKpdW2C8mL9EwLwgzOK+3PJJZek66+/PuxYeM60vzfffDPuEX7EG7NVzwjo+elZKPs2ARMwgc4iwHvGygRanUDDBRABKnbKNcVKQgEjHWPHjg2hA8EDte666yamuHz5y1+OnaWKHRTiQOUdesyahy//xU5hfo2f/Do3R+Q1/imcdAXLr3Oz3HNdLHI7zITjRwcYPyyO5ys8vx//+MfprbfeSnfddVe68MIL04EHHpgOOOCAhDDCGgWENkZGJJAprnJpdZe/Yr7a6bpc2VV/brrppuDPKNtVV10VU/xYA4JCRxhEOKQ+SujDTeExW9VOwNxqZ2WfJmACnUdAbSTvLVS591fnUXGJW43AJ/NyGpxrPSQkQyeNazpqPET//Oc/49A2vvJvs802sf3tf//3f6dJkybF1+Vvf/vbackllwy/hOWnkQI60frpgZSepyk7FbPaddFNYbrTCZc3AopHOuFzc3fxyZ04KSNhESIkOMCAUQ3cWRjNuhDOoOBrPKMia665ZjrooIPS0KFD0+qrr57+9Kc/xQndio/4dS+wQ/UmfxGwDf6KZZewx5QrBBDcv//978duZfBiUwCFYWqc/OcsZW4DPC6CCZiACZhAkxDwu6VJboSz0ScCDRdAJCyg02FD8KCTfO6556a11lorOsd8VT7mmGPSxIkTY+rVPvvsEzs+UTLCKQ4634RXJ5z4+MldJJQW1/mDmpvlt5xeq79iWHVIexu+GB/XxEl8ilt+YIBAInf8IFCwsxdrQk477bQQOFg8zZbCu+++e5w98sMf/jDWMxCPOBJWeZaudDpRz3ncc8896Z133gk+I0aMKHFiitsyyywTeG677bZgqXrHPWGaXPGedSJLl9kETMAETKAxBPSOkd6YVByrCTSGQF0FkGLnlQ4Zdjwc/FhATgeY9RtMD+LrPIe7MYXl4IMPToO7tsxF0ZHGP2HpJKujXOkhk7sQ5dd5mNwsv+X0Wv2VC4tdX8MX4+0uPtwps6b/aJQIYWSzzTZL55xzTpyFcuKJJ5ZGR1jEzw5bYkwchOtEVay3sBBzHT7IWS1MeYMz9Zpd1zbccMPAxdbP7777bgkdI09sq4yqlWkxD4QtZ4d9u6pOK2+73sdOKhd1Nm8v8rK7Puc0bK5GIK8rublamNytN2Hy8DabwMwkoP5VTQJItcqduxEp13S6+KlTzDkd3/3ud2MtwuWXX57oCHP2xXnnnZdGjhwZX/LxT1h+6kgrkzMTTKunBTO4o/Ojswxb1i8wheixxx6LtTYrrLBCTHmjU40ggj+4o3Qfiiywl3vRrVWvKVOxnsEChoxisOAchSA377zzBhuuCceICP4YuWMqIfzOOuus2J2M9UtwF1PC5KrIWGzxI3MxX3n4djar/O1cRpetPQjwjBafZZWsU59fld96zwgwpRpVrDfl2kPeOyi5FcOEo/9MoMkJ1CSAqJGlLKrw0vOKT8cNezpd/BA8mE5FZ/faa6+NBeaMguy3336xqBy/dNJQ+CcuxSc9HP3XIwI5RwLCVoIIZs5LOf/88+NkeNaHbL311vF1X4KIwnN/dJ/RsZdOvHLD3KpKZSrmn7KxtuPhhx+Ocm/Stc0xflVfMbPOBq6zzjpr+u1vf5uOOuqoOAjzl7/8ZfjFj5RY6RnBLXfXCwX/slcYxdEpusrfKeV1OVubgOtra9+//sx93sZr8xy9Y+SW1y/Z8R6xMoFWJ1CTAEIh84eg3DUPBJ0ofpMnT44zOhA8brzxxnTRRRfF9Kuddtopvg5r8TRx0iEuKj1k0ovuvq6dAPcDjugSRNQJXmWVVWItzkMPPRSdaQQRDkhkyhZ+uD/Fe6B6gL3MteemOX3m5RAb7Nh+l3Ly22CDDSLzcMQNu2HDhsWBlxxUeMUVV4SwwtobuYs9AYt2nCsyYcKEEpAXX3wxmBMvCp0wnaBUTumdUGaXsX0IUG/56dltn5K5JI0moHpz9NFHx0wE0uM9TV2qVp/UVkqv5rfRZXD8JtBbAjUJIKrc6MUKj52EDyR3dmIa3LWWgy1hETyY8sNp30j3uONfi6fJtOLOC6A0pOduNvecQM4Rc96Jhj87OrFLFoIIJ3xzlggdbhZgqxPNPc5VHmdu36pm1UPKyw+14447pgceeCA9+uijadCgQWEnd/xzEOSf//zn9LOf/SwOJOTMGinxkS6BDoEcvkzZuvnmm4M70xP32GMPBQ1d4WawbNMLsVfxiteyt24CzUhAX6x5Zqm7rr/NeJeaM096L7DWkA+AvH817VfvAE3NogSqY3JrzlI5VyZQG4GazgFRZc8rP2YaXqR1zHfccUdMraJjdcQRR8T5FGwPi5K/cqMdiptGW+basm5fvSGQM5aZ+4MZQeSPf/xj2nfffeP+cRo754z85Cc/SQsvvHBJ0GyneyXhWSxgKvOyyy4biMuVV0LK8OHDS1vz6n4QPg+j+s9Ob6eeemqMChI3fjhY8swzz4xRJ+KUX8XVCbp4d0JZXcb2I6DnXSXrSX3O2wnMKIWv5oa/ov9ydvKDW57PPA3cKin5q+TenX2ePn4rxSd/uTt2XFfSiU/hMFdTeby5v1rCK6z85tcy53HmZoXJ7cqZWV/Ih1m20OfHTpa8e1l7qKlZvKsUH2aUrrvLR7k0bWcC/U2gphGQPJOq6EjlCBRsUcquP+wIxELnRx55JP385z+PsynoTKEkeOhhkY6bzMSLWdfS8WPVWAISImFOw8Yhhnyd50DIP/zhD7FehylG6nTrPrXDPVKZRFhlUhmpw2rs5SfXcZMfwuj5kB8JFGeffXYIcgjqCB8ffPBB1HXOvyEMC9pRxfCKp5N0M+iku936Zc3bT+pu3oZ0V7q8rmMuXit80Q37WuzkR/Hquhiv7Iu6/PVWrzU++cvTwQ5VSZebwlbTI6Iyf9XCyE3Byl3LrZKuMJX0/P1DnwphAzu2z2c0hB8CCXUKe8WjOlYpXdubQCsQqGkEpFgQOl08KEyvYl0Hi3WZbrX99tvHg6KOGA2zzMTBw5Prlcx5mAjgvz4TKMc0t9O9ISHuL40d95YvMEceeWTsmHXYYYeln/70p2muueZq2y/14iBdwnOlG5C/QOApRXgJH0zhYooV07UWWWSROECS6Vv4YQetZZZZJs69ye+H4ulE3Rw68a63bpn5CPf000+nddZZJ02bNm2Gd169SqX2qNjGEH81O4VTPop+82v5mZk6+cvzoPzmduRH/op6d3mtFF934eSu9LguxpW7yX9vdTbsQSGEEK+EDabs8hs5cmTMTPjqV79acsN/nieZsbcygVYg8CkBhAe/UkVWx5SHg2kje+21VzwY6kTR4cKNn+KpFFc1OL0JUy0+u01vqHIW5ThjJ3vuJ2e2nHDCCemLX/xi2m677WKxNet7WICt+y3/edydZlZ9z8st4eWQQw5Jiy22WPDjGZL922+/na6++uo4AweGer7yODrBrM6GdNenTrjr7VNGzgTimV5yySWjA9mpz3H73NGZVxLVlSlTpsSoOCnTDvLTe4L2cNy4cdHXQkBh18WicptZJOLrViDwKQFEFZkHQGYKogfl/fffj7UeCCC/+tWv0gEHHBCjIfraq0LnYWVnvbUIcA+57/y23XbbNGnSpFizwFqRiy++OEa8LIR8ck+L9V3PC7tbXXnllemIrnVRjHrACzdGEE866aTEDlqMMqH04inG9UkKnfNfbHs6p+QuaSsSoL4O7tp45a9//WsrZt95bgICfNyj/iB0UJ/03l100UWjj8XBzVqH2envh/66XX4v1Z98CCACK51kVMmxo9NEh+mNN96IE8yvueaadOutt8ZXcfzysEha5zqPh2ur1iTASBYKnTqwTNdUIXbo+M///M/Y2Yx1DezopLmrKqXv//RpEewsBj/ODIEhimeJZwiOc845Z+JUepjlz5BYdpKuNqeTyuyytjYBnlt+6jCq7dMsgEaXjvSkeH50LXN3z5T8K46e6n2Nv7/Dd1feRuePekMa6LwX9I6Q4PHNb34zpu2ST9y6y0935bF77wmYfe/ZVQo5wwgIgNWAKgCVngeDk8s33XTT9Morr8R5B0OHDo0HgoZWja1ukHTFYb31CXCPqQvobDKw4IILpt133z099dRTiYP3VG+kt36J+14CWKhjMn78+DR69OgQ1s4444zQWYy+xRZbxIGd88wzTxoyZEj4h3EnKbcXnXS326usqrs868XnVm6NLHExjfw6N1fKQy1+KoWtxb6v8fd3+O7KWI/8UW/4iMf7lam6+++/f9pll11mEDzwwwcq/Fg1BwH3dfp+H0IA4SESzPyBklTOAimmibCAduzYsQnpnAch/2KrcIpHet+z6BiagQD3Vy9Y6gVnV7DrGYviWHjJdDzs8aN7L70Z8j+z8yBWjG4grDPacd9994Xg9rWvfS2NGjUqFhXyLDH8vtJKK5X4kddOZjez75XTM4G+EuBd6Ge2rxQ7Lzx1hvcm2/DyIY8t8Fl3yftWMwvyfhb+UepvdR6x5imx70Hf70VpBKQIUwIGi6OYPsK8f7aD+9znPhcPTN7g5g8F8ehaHdK+Z9MxNAMB1RF06sc222wTO3QgnGKHEII9Zv2aId/9kQfKz3MwcODAmK547bXXBhNGPBA6eDYuu+yymIK13nrrzZBFwhG+ExRlVXshvRPK7TK2F4Fyz2w5u/YqtUvTVwL0oxA0OB+KUXDqjPpNueChdDrlvaDyWm9vAiUBJG8sMVP5X3vttTSya/u3VVZZJV111VUl4UMPQVEHVR6P3NsbYXuXLr+feUmpHwgbfMkfM2ZM6GzPe/jhh0cDit9KYfN42tms+o/AkZ90zgsGfghwUjkrwuXX8tOOuhhRNpmlt2N5Xab2I6DntVgy1+MiEV+XI8AUdwkfuJcTPIrheIdY9R+BTnk/N5pwCCCCmevsdrXzzjun119/Pd11112JA9M0KqJMyb+u0d3o5jRa35zfz+L9lhCyySabxK5YO+ywQxo0aFDavWttiIaPW59AfUrAswM/pmZpepZGi7jOOZNi8bo+uWjuWOBjZQKtRqDYLhavW608zu/MJ5DXmdxcKSfl3hmV/Nq+7wS4Jyjey5jp35TbDrnvKXVWDLHaNe/s0Cmicn/nO9+J07DvvPPOmLcu4UM3QjejE3DlZaa8xetOYEAZ83qiMiOE8DWGdQycms6XfqYb8VWn0leaTuQHJ5hI+IAfdvm1mHaaXqxXnVg/GnnPi88hfIt2jUy/3eNWp0TlLNZn2Vs3gUoE8jqTmyv55xnWr5If29ePQP6M84F1woQJEbnb0b4xniV/2UvIOOuss9I555wTC87ZmQdpj84SKr8RfUu6NULDJ28QdJ1za42SNCaXcKBuoCO07rfffmmrrbaK3dLgVu4BzXk2JleOtZUI+Flq7N0qCrk8f0W7xuag/WN3m9b+97gZS6h65za093cnZ1fOjJ1mK/z+97+PvvHgrnN/rPpOYBYqMID50ZF85pln4oTzI488MrE4Np9Ko5ujMH1PvvljUFlVdg7j+8lPfjKDUNL8pWhcDsUHnYf02GOPjd2x2CULO1QuhIhj43LkmE3ABKoR4Dwntvpk9zqrvhOgTXO71neOjsEE+oOA+imkXTTr2aZvzJleP/jBD2LHMn10zf33R95bNU21l6UpWFgAk6/YK664YjrkkEOi48i0EZTcFbDTwNOJ5rd71/AbZ6LAIe9Yt2pF6Gu+VS+IhzrBwXoXXHBBuuSSS9JFF10Udqoz8oOe23FtZQIm0FgCfExC0cYztXa22Wbzc1gH5LR7nfY+rAM2R1EHAnqPuv71HaZY5jGJ65tvvhlns+DGjrBsmyy33L/NtREQuxBA6EgzJH/SSSfFVrvnn39+mn322WeIiQDcIAWcwbHNLyg3Eu9xxx2Xxo0bl5ZbbrlPdazbHEHF4qk+oFOHGAXZYIMN0kEHHRSbGLz00ksleyLRQ65wFSO2Q8cScN2o/63XSPa5556b/vKXv6Rhw4bFs6jnsf4pOkYTMAETaH4CagOL7x3s6c/QrznwwAPTAw88EIWZY445Yj0nF8UwzV/a5sphrAEB8FtvvRWQDzvssDR8+PCYeoW9FDejE2EjnFHuW2+9NR188MGBY4kllhAW6wUCGppkCh9nxhx99NHBT3UHXQ98IagvTcAEGkCAlyjP5YMPPhibRGBecskl47n0s9h34DA0x75zdAw9J6D3qutfz9kphBjqGl3PNG0lH2xOP/30hOCB4mBhFO2qVd8IzCKIZ5xxRswJPuCAAwJ+OeEjr+S5uW9ZaM7QlE/CB+eh7LnnnvESp7LON998zZnpJsiVuDEV6ze/+U068cQT0/PPPx9fETRlzUJIE9woZ6FjCPBMfvDBB2mfffaJlybPIduqW9WHAO1ZuU5MfWJ3LCZgAo0iUK0fyzP98MMPp2984xvp0EMPLa2Zm3vuuSM71cI2Kr/tFu8srPFgUSKAf/GLX8SLiRdULoCocZUOhNzcblBUHioY5dxrr70Si88xY8dZF3KTX+ufEIARXw1QnCPDIXyMguTK7HIaNptA4wjwgYk2no0zxo8fX1r3QRtWSfF8WvWMgJn1jJd914eA6l0n9MfqQ2zGWMpxgyk/2s4vfelL6Y9//GNae+21S2t+NQJSLuyMsfuqEgHV25hjdfnll6cPP/wwtlCV8CEPlSJod3tx+NWvfpUef/zx9NWvfjV99NFHUWwvQKp+9/XwssiVqVisLZo6dWppFEQPLv46vZ5VJ9kZrqoP0nn2UK4bfbv/cORjAFuqsyHErrvuWmrDqo3i6j70LXWHNgETaBSB4jPqtrJ+pNX3Yx3rYostlr797W+nZ599tvTRXVOx6pdi58Wk+hsCyFVXXZU4zXqBBRYojXzIQ6egyR9gvbjvvvvu9NOf/jRddtlliZPhYcKLe8EFF+wULL0qJ5w0gvaVr3wl4rj99ttDzznjr1jP5C69VxlwoH4noPvHVyRN86yUKfmVXqwTlcLZ/hMC4saVzHqJPv3002m33XZLV199dWwsgjtCyeKLL258JmACLUpAz3mLZr9ps827ivbxb3/7W/rd734XfT9Ys35Oaqmllgqj31Mi0nt9FjrWN9xwQxo9enR0BtktpVPBUtF4caNYlD9q1Kh03nnnpRVWWCE9+uij8XKfddZZY6tZ/HQqJ8peSalhhA1mpmCtvPLK6Yorroggsi+Gz8PhVslfMZyvm5sAjTk/GnbdY3Kcm/UcSVeJiteytz6dABzFSWZ0frTlW265ZWIUd/XVV0+PPPJIBORezDXXXGFW2Okx2mQCJmACnUVAbSbvqhdeeCFtvfXW6frrr5/hQw1+UPPMM0/ouo4L//WKwCxPPPFEevvtt9OIESMigk58IVGRKDc/zHy9373rvI/NN9881jFMmTIlFlIDaOmll47dEBBUOpFVd7UsZwJLfl//+tdjCgjT/DQyonhwRxXDFe3Ck/9ajgCjhxMnTgwhhHssQUT3W/c/L1g1t9yfzZ88N/poIm5c8yJlCsFnP/vZONOJDyrPPfdcIOOjgKZgKYxZmoAJmECnEqAd5N3Eb9ttt00/+tGPov+n9xX9ZBTTr2g/3f+rT02ZhYWJQOXwQVTeQSzXOahPss0Vi17CVDZe3L///e/jvI8//elPkVHWfuCGGjhwYCzkzDmFg//KEoDtZpttlt599934AlusU7iXsysbmS1bjsAJJ5yQhgwZko444oj0zDPPzCCI8Ezp2VPBqAv6Fd3kx/qMBGiL9AxJ+GBdH1MIrrzyyvBMG8ZmIzBlqq0FkBkZ+soEWpGA2kjprViGZsgz7SObdbBhDu3kr3/96xmOouADjhQ7fLr/Jxp90weMHTs2rbvuujGsxMsrB9tJlVrCx0MPPZR++MMfpvvuu6801KaTz0HNCAhciqz6dhvaM7TqD8Itgh3T2FZbbbUQ5riWwh8dKHQJerjlHSv5tV4fAuqw6h7VJ9bpseh+LrTQQnFP/+M//iMa9f333z+2tB48eHDUCfwV80KeGpWv6TlsLxO84Mgzw7bXbJrB1FoWUdJWvfzyyzHSjR/tgIXZnNurHrg0JmACPSNA+8jU+ttuuy39/Oc/j34K7Sj2tI+0neyCivr85z/v6as9w1vV94Cbbropfetb3wpP6lR32otJ5WU9zMYbb5z++7//Ow5jnDZtWox2sIe+lPbPJ4xVdQLq3LAbFuZ33nmnFEDMZSG/uWCCm+zlz3prEMjvL/eQr0s8X0cddVScD8PJsnvssUdadtllS2uqNCLiZ6t391gvTDZ+YNtdRh75ssfLlbZMwj2baPCcyX/vUnMoEzABE2htAurzImTQXnLoINvs8g7inYVS3wWdw5W9fq5+93zAK6+8El+lBTrX65dMc8dEJeSFvPfee8d+zwhk+WJ8FiWpU5R/PWzuUjVH7uAG21VXXTUxuiSOuWAhO9baMFULN9k1RynaKxfl2Ne7hNw/0uEQTxSdX+z4ssRaIASR4447Lh1++OHxlZ46ghtKYePCfzURgC8MOc+JDyZsf00bpueIUREpt2EiYd0EWpNA3oa3ZgmaI9dqH+nzMWq81VZbpddffz0+2tB+0qY+9dRT0Y7ilxkwKAkuzVGK1s1FiHhAzlU7dgAqlYkXNx2fM888M7EdMcIGD7ekX8JxhoWUtuAtNgCV4le4ot5T/8XwrXQNX85OgTXcKLuUzNjvu+++6brrrpvBXf6stzYB1XcJ+7Q57733XnSYF1544agXNPgo6ov8t3apZ07uJXxce+21seMVU0YZdcwVG41IsQZEz53srPeegOtq79k5ZO8IFJ9f18HaOOac1G6eeuqpsVaO6cK0jWo7eVfRL2EUWYo1wKgif7lbr42A+IUAAmip/AbJrtV1lUmFlvCgCsguPZx2zl752L366qsh+cpdOyDAQVOwciaKP7erZpZ/6dX8tosbLPOOpcque8E1i/8fe+yxkpCCHUp+MJezwz5X8oNdHjb3Y3N9CIi1OOua2LFj7Qfn6XDvZUdd4LfBBhukX/ziF/GsseW1hH49dxHAf1UJwBu2fCRh95YzzjgjhH1GthHyxJL1V1K8RHW/ZGfdBEygdQn4ea7t3okTH7t43zzwwANpn332STfffHNaZJFFZvjwpVEONvRgpB6FkGLVdwK6D58aAZFD35NonhgoEy/qvGy8mLlmusIOO+wQmeXES+YCMmdaYXiJ8/WQlzwVki3YUIpL8UoPx27+FFZ6N95b2llcKCvMpYpl55oFXvys2ocAO8n94x//iM6wviQheBxwwAHpy1/+chyOxyGV1BP9inWjfWj0viR6jhSDrvWSZD0NU9sQ6Pbbb7/wlnNk/Y3asCWWWCLcicPKBEyg9QmoPWj9ktS/BDkbtXkagWfNL7sFoldSd955Z/T9cKftVByV/Nu+dgID2IKXA6q22Wabtgabv4yFhxfyj3/84zT33HPHLgcIIwgfqmCEodOEgML6BYbm2FcfpfhyXZ0BxY9OXPpR6TErTO6vHc2UFcZ8bZgwYUJwpOzYF5XsYCjVSaxU5lbWdb+4x9xH6Zj5bbjhhiF4/Pu//3t8fZK9PgZQ9kr1o5W59DXv4ipdnDS68V//9V/RPjFSix1f9vDLT36/+93vphtvvDH4spASxbNpZQImYALtSiBvMykj7xe1m4x8rLPOOrHrKX0U2kPc5Y++H322fPSYKfjyU4w7AvqvRwQGsAXvHXfcEYHa+YWUVxZVwAsvvDCdcsop6cUXX0zMQ6+kxIXdD/TyViXMw8hfnhb+cr+5OQ/bzmZ4sANWNXbiokYgZ9jObNqpbPk9U71nT/U111wzhrAZ8VDnGOED/zTw+MVsVZ6A+OgZwRf8YMfWkYcddlhC+FhuueUigvw+YME13BWPphHk8UVA/5mACbQkAT/L5W9bkYumXv35z39O5513Xml9r/oduX/aV+x5h0kxTQuFv9yv3K33jMCATTfdNJ188smxTaMW3xRfYD2Lsjl9q7LoxT158uT0ta99Ld16660hfCDtqnOEX4QUKh+7+LAwHbXUUkvFCEl3JVRaL730UoRHwCFuvgDDlnjbkXGRi8rICcyYNXqke5D7l1+xk577sbm5CeT3TGaGt5dffvkYPcROzxXPAOZcKQy66kPu3olmcciZ8PxwTds0cuTIxIcUGOsrHpzkBzNTs9gFi7hY/zHPPPNgbWUCJtAmBNROtElx6laMnAttIv0wdrX6xje+kTiCgkXn+iBdTJQ2ljCsCUYhkHAIoVQet+ys10YAdqhZNtlkkxgBePbZZ8OiHaGqsFQmzOijR49OP/3pT9NGG20UFZCpVwGkq2MknU4SX+5VAflyqD2gw1OVPyr1P//5z3TssccmGDOFSxWaYJg7RSHkwZ25/yix1n3BDh75NXZWrU9gpZVWinUeevbyEQ89A9Lz0pazy907xZxzkJnnBPPuu+8eL1KmiNLe8HKlzdIPP5g5A4SdsVBMIWBHOj9rnVKDXE4T6FwCajNpHzGzFo728ktf+lL0y/howzspbw8x4x+d34MPPhgAZ5999jjEFTfeZ4q7c+n2vuRiN8vKK68cL67x48dHbIBtN6XCUpmobL/85S+jYqHnFUn+0FX53nzzzTDDRDtg4VaL2nLLLRO7+5Am+0sr/VrD15JGK/i55JJLYsExXxvUEc3zLR7in7vZ3NoEdG8leKg01e61wsiv9ekE9MI87bTTYvoVh6aiEDRypTYMu7feeivO18FMGya3avcAv1YmYAKtQcDPcvX7xPuHj8m77rprrEfVcQr5rJc8BvzTpo4ZMybWB7NWmpFkNkzBjXC54p1V7NsgqFhVJzALX8OYhsUCxXZVVAwqAxWHrXaPOOKIdNJJJ0VxcaOiUYH0k38e6kmTJpWwMM0BO6ZrqXJV6iwpTs4XWXHFFdOSSy4ZYQnfCY0FXGCAAMf5BNppDHuVv6iXQNvQNgR0j9umQP1UENokzV++5557Ejv2/eEPf4gNNHArcuY5o41CZwokG2yg2GWOqbZyKxaHuFCEk8rNsrM+nQDse8qop/6np2aTCXyagOvTp5nIhq39mYlCP4yPoahzzz03jl649957o+3M+TF1Xlvv8hGZNpH2kzaTD8n0H3HPz4ejDaC/o3ZYH4pIC8HFqjyBEONYHMpCRiRE5gdzMwSyfLDWsaUsEgY45G7rrbcOQYT1CJSxuOsVJcNe62FYrMQ1v2uuuSYW0zIPkGv5JQ1VYNJCUWk5UfOWW25JxxxzTLjneVH48NyGf5QfFmPHjg0WbHNH+du93G14KxtWJD0zDUugjSLmWeLHmSrsIsbHFI3Iqs2huPkzpqmOF1xwQYTFH1/02Gq8uOkG4XhmiVdx8MLlWs+s7NsIa92KIka1RthT/7XGa3+dRYBnEuX6VPm+M3WKjVAuvvji6NfRzjGCweJydkDN+cGTNo8P8/RZEDhoN+GLGz/6yeh535E1dvQP2e2TDz5nnXVWGjduXAg6bPLELoSHHnpopJW315Vz3SEuXTfj465F0tTij7u+qHVx/fjjLuntY+xR0uOiBf+6XqIfX3rppR93LTr6uKsSRTkpa9eC8ihvlxRbKhVl5de1s8zHXSMXH2+xxRYl/4Th17Vr2Md//OMfP+5a1zBDOF10nRnycVdljMsugSfCdK0FiWu4doKCOQqWw4YN+3innXYKs+w7gYHL2D0B1QeeN54t2iGU7LuPof19qE3qerl9vO+++wYn2jF+XR9RPj7uuOOivcFfrrpGST4+55xzPt5xxx1LYWBMuBVWWOHjE0888ePrr78+WBfbJe5D18hlHl3LvwdmKEydLsTtqKOO+rhrdLxOsToaE6hMQHWO9oDnuesjZ3guPv+VY7BLkYDYFXX5k325a9z4dU1z/bhrrUjcE/qNxx9/fLS/XWvvPu6a6h/2XaMuEYXuYTFexd/OuspOn5z6O6ALQpwAecghh6SDDz445sgx302q64UlY0vq5J/pT9/5znfSQQcdFNJs1w2OaVTvvvtuSMRc409lZXSEqQqc4vyb3/wmht6QWvHHUBxDappD2NVZijiZatRV6UpzA9lyloXniy++eAz9ERbJGqX0WhJoDZmmTsGL6W6cMdNV2cqGancOZQtty08RoB6g9Px9ykMHW8AEPosttlgsOO8SQkptWFdjHms7aK+L7PiCt/TSS6cf/ehHMbotxjyXtF/vvfdebIuNvUY9mC5Km0U4phcMHz48Tlcn7mL8HXxLXHQTaBoCfi67vxW0cfRJxEptISFpD2UvHXv6dbrGP/4UB9e45WFpbxmRZpOip59+Oo0YMSKtvvrq4Y82mwNimQ2yyy67hJ3iIK1OVgOACAxOJmaqECcXM1zEDVCHuZUBUVHWXnvtUqe/0o3HHoX/VVZZpaYiUyFhxG5X2223Xfr5z3+edu/amQb1+9//Ps4YOfDAA2cQcnBTxcbcbgqO6tAg0HaNfqSuL67x8FLXctXOHPJy2lyegJ658q62FQGeE4QBVLH9Kl4rzNChQxM/uUsvxoE9U0Vpv9gmnPnN8803Xzr77LPT9ttvHx9cNB2VsPkzm8epdDtVz7l0x8DcuiNk91oJUJdQPal/tcbdLv5g05O+bCX/xTj0HEtYufnmm+PDDpsbrbHGGvGRm2lafOhG6RgCwqnf3en3LdaAAIcv+qwD+eEPfxhb1PIVTF+ycW9VpQeU/FMeKdlTAfSTG/5UuVRZcMNeFQZ7KuT9998flY21IpwrotMzqYAotnsjjEZKwvJff0ojt2tlM+WBEVy6pqnFCaIXXXRRsCyWq93KXiyfr7snoGdJPqkTqKK93DtZZ7RDL0DaEilY5bz0XKkNw5/cxTf3w4gum5CwRTYjvijOLdpzzz3jpzQVR3j41185u9y9k8xiW0uZxU33oZYw9mMCJtA7AsXnrHjd01gVXs8x4REoWJjOTp98zJEf3B544IHoF40cOZLLksrDlyw7zDAACMBCHX744QGRHVZYsM1LDKUv1znUcGiBP91k8q5ykG3ZF4tQ9Cd37PUyFgckWxaE7r333jMIH6Rz5ZVXRhoSRBRW8aFXykPupxXM4kFeKfujjz6avve976UTTjghcQ5EudG0dil7K9wf57H1CWjbR7VD+TOXl07PFXre3uV+MCPQEGfX+oV4Qf79738vfSRhgfrEiRPTEkssEX5Iq1J6xXh9XZ4A/FC6P+V92dYEekZA9cnPZ3luOReZxawYQu7Y5+aiv2J4+nZMwUcAYSYRox4aFSEePsJyhtwmXefBKV7pxbg75VoMY04MFwgb7Bbwl7/8Jd1www3p9NNPjw63hJBWB6MCVyuHKgW6lMyElxkmXJ9yyilxwNd//ud/hhsVEXumNHAC9DbbbBPDbopXcbabLrY8dJR1jz32SOutt14IIeWED3FsNw4ujwk0moDaIT1z5dKr1t7gxg/hgwNWmXb7k5/8pHTCL/EiuAwePLi0EyB21dIrl4dOs6vGR8yLTKqFKfr1tQmUI0DdQrkulaMzIxcxEjNC5Ga5Y5+buUblfj+x+eRDDmaNcrBtL/7U75k8eXI6//zzoz/N+hDcysWt+DpNL03Kp/MMNL5Y/9d//VfaZ599EqejYy/wAqfrdoOl8kmnfEUzZecFzUJOBA+mWCHdYi/BhIXrTG1gy1/s2kWIq3bf9cCx9uWuu+6KbeioO+W+wuZM260OuTw9J+D60DNm3fGq5o4bzyqKQ7VQI/81NUDhqj3nEcB/NRGg3Rdr2kGm57bLu6AmAPbUMAJ6Rsu9XxuWaJtErHaO4uTm3hQP/tyLyy67LJYxsPicOFk3x7PetSNWbATyla98Ja7xjzs/3cPepNsuYWbJIQCHBpNhpHXWWSdOr+YrGbAYspfiulNUzgczlYryI5y99tprsVc0bHi58FXxb3/7W0jD8GFuNUz55SqPM7dvVrPyW+mhoc5QRna9Ytcd9sAeMmRI1KVOqivNev+aPV+qX82ez3bIH6zVHr3wwgvRlvcsBYgAADJFSURBVLH+D3vaMbVv7GX//vvvt0ORG16GYv2FIT848xGGHcfYYYy54bwnrEygpwSKdUzv1aJ9T+O1/+oE4CvGYp6H4Bmn3WTKPTOImH7etSVvbEy07bbbxi6CXVull9qDPGy5+HL3TjDPUoQAUKS3v/71r2natGkxjYgOJp3r/GZ0AhzKmPPBzA8OvGBQCCKwofIhBVMZWcC/0UYbxVxATuEsKsWpil10b7Zr5Zd85WauKS8vWU5n5tAe7QRGncHeygS6I1CsU935t3vPCOTtDKzVjrO1L2baLeyZu0z7z3aRHKDq57c2zuKJro8xhGSN4Nldu4mtv/76cepy1zkBpfYTv1YmUCsBPbe5LnOtcdhfzwnAmJ/6e3kMsnvqqafSpEmT0je/+c1YB9J1nl6soWYH1KOPPrrsh1g//5+QLO2CxSWgUTSivJy6DqqKrRx32223aEj9Qpq+IH9w1xxpTtc86aST4sWNsMb2vSzg32uvvdJyyy2XbrrpptLCI7iKL5VPFRv7VlUSTDnrg5csi/GPPPLIsg9cq5bR+W48ATfGjWWsdkepIGRgN7Jr6hVtGLtf8eWO9mvKlClxbhEjmVbdE6DuwhMFU96RnJR8ySWXpN/+9rdxrgpu2M8zzzzRkaHdRBGur3Vf75KI8F9/5exwUlrF+pCHrWSuFLZSWkqvN2lVykOj7FU24s/zm9sr7dxddn3Ri/yUptLRteqKBFw6v8WwfcmHw1YmAGfx576IOzo/Ptigf//7349z3/KYdL/ycLjr/uZ+O9EcAoiACgCwAbfsssumu+++O7aZ5brrZN3SCn/8oDoRJLw0SnTqqafGS5vON/voo1hDwxoQRgQ4oCbnm5vDc4v+aeQD4YOOC3McWf9BPaFOqH60aPGc7ZlIQG2I9JmYdNsnpfZGOgWGM9cclspuh2whTjvPTi4cmrXaaqsFlzxM24PqZQFhSWeQ0SPWBbKJy7HHHhsHsOLG6Djuahd1yK/Yus73EnwHBqOu6L3KgmauqUdWjSWgZzTX4a5nmBFkPjizeyD9Ij3zuOcf7RVe4Rqb69aIPQQQVWQBQpcQwsuIA1Y23njjgHnGGWfEdCOBbo1i9j6XqixFnRiXWmqpxKEzUrxo4Map61IKp2sx1nUr6rr39957bxzy+NWvfjWEU6ahWZlATwnwjFg1lkCx3eGa9mrgwIHx5U6pcy/UZhXDyI/16QRgOOecc8bpxxy4ypRc7OiEIHTwU/1+4oknYlEqnRIxlk6M8ifuup6eWv1NSqtSzMpD0V/RXteKp5L/onsxnNybVe+uXL3Nt+LNeWDHtXTFTR+Ds3p43+b+5W69sQTEnGecdXTsGks/UM81qUtQLJeT/F7LXM5fJ9iVpmAVQXANUDqbI7uG6u+7774Yrme+21VXXRUr/mlcBboYvl3gUS5VOMpULCcMZC8W2Mkfdnkjgr0ElQjYQn+Ugx8PHovt2WaYwxeZ58wXQCkxEwPZWzeBcgRcT8pRqY9dtWdR3Gnj1XZhJ/v65KC9Y4Eb/JhexXS2Z555JgpMGw9HdCmumdos1rJvBb3aO7Ba/vXuK+enmls5/51opzpEnVFdwsyU7/yedCKb/ioz94TNJMaNGxebSzCSzPquO+64I6277rozZCuv47kZT8TT6aokgOQgclD6ksNIyJNPPhk7Y3Fq7o033pgGDRoUDwUg8zB5XO1gVkUpV8Z8iE1llSSscEW9FV9ACFWUi7Iwv3mHHXaI9R4cXoldUajCrhwvMbJuAiJAPbFqHAG1P8UUZE8bL+V7IRK163QGF1lkkXTxxRcnDnQ8/vjj06WXXhoRwFgf8pZffvn4WFN7zPZpAtMJ5O9TOr9sdGPVPwT44MrzzHTLE088MQRChBK1qcoV17pvRTf56WR9+psno1AEpQZ0cNfC69tuuy22E2SnJ054pCMKYEHOomk7Y5FLtQL2xG+1ePrLLb+fEj7YXo5Dy/74xz/GCcqHHXZY6atMufKWs+uv8jjd5iWgekKds+pfAroX/ZuL1ko9/yLNxzl+CCInnHBCCCUaEYct60T0ISdvY1urxM5tfxGgLvHB4O23344s6GOm69LMvyNMwa9FuU2tTKmsAFLOO5Wer9wcusfuTkcccUTaaaedYtiJXVQYglZHtdzDUM6uXDq2qz8B2KN68iDgV6MaNHL/+Mc/0je+8Y00ceLE2B1t8803j/utBrAncde/hI6xFQmoXkpXGdxWiIT1ViCg+qo2E51dASWIsDkHH+teeumlED7yEadWKJ/z2FwEeOdSx/zO7f/7Unx3+Z707J58spVVjWHU2WT4iZ2eOHzl5JNPTquuumocvod78YYo6vzG0LGt5E/+rdePgNjnzDHn18XUJHzg55RTTklrrbVW7HbFHGeED7kTTvEX4/C1CVQjUKw3xetqYe1mAs1AgPaReqt3I2Z+eschhFxwwQWxVeeOO+4Y60XId7W2txnK5Tw0PwHXof6/R3repfd/jlorBz0SQPKiMdoxevToOICF3T+GDx8eoyJs5cjNKKd4YPhJgi/nx3aNIVDuASl3n1hQieIe/fOf/0wjuzYgYFcvTSdYYokl4iWqe1gujsaUwLF2CgHXqU65061fTrWrCBwo1V2mWaFkP2LEiNimXJt14M8dyEDkPxMwgQ4l0GMBRA0sHVAaVzqk7Ip11llnxboA1olwiroUgooaYTXWcrM+cwjkLzqZ8xcgdtwjdKYHvPbaa+mggw5Kq6++emwxyQ4PP/jBDyKz+PMUgplz39o9lbwutntZXb72JaB6nJcQO9pYfajhmp8EE/zqXZqHs9kETMAEOoVATQKIBAigqGFV48o1jSynpT///PNp1113Tdtvv33acsst4xRYGlwJK2qEOwVus5RTLzrdO91HdARE3Plxn88999w4zfO0006LnVw4qGyllVYKf/jP4+LaygR6S0B1SeGpn1Ym0A4E8rqdt7sqm+u6SFg3ARPoVALdCiA0lAgQUmpY1YByjZnOK1sR/u53v0t33nlnmjp1alp55ZXTd7/73cQBTPoSRDz4VXjFa72xBHLe4s+9Q0Dk+pprrokTkBEgDznkkDR58uS07bbblu6tvtzp/ktvbK4deycSyOtqJ5bfZW4tAnoHVsq12krp1G+ZK4WxvQmYgAm0O4HpkkWFkuYNZd4xyO0x00HFnd/aa6+dbr/99jgd++qrr06sEeG07Lvvvrv0tZ0wdHytGk9A9w0d5giD/N5///04SOfzn/982mqrrdKKK64Y6z5+8YtfpM9+9rOl0ZFcAG18bp1CpxFQ/czblE5j4PK2LgG1qZRAdTkvTdGOel60y/3bbAK1EHB7WQsl+2lmAt0KIHnmu6vwuPOjQZ599tnTLrvskjg5nQOaXnnllTgl8stf/nJs40rnN+/YavFznl5uVoMtvZxbbtfp5qJwp3sD86effjodc8wxcSrvXnvtFZsJcNL9hRdeGDuaiZ1GPXRt3QQaQaDYrhSvG5Gm4zSBehHI62tuVvy12sm/dRPojgB1Ku8Llatj3cVhdxPobwI9EkBqySwPgh4OOsHs+rHddtvFiAinp3Nq7Je+9KXEIS5s5UtnGKWFzYQpdp550PSASc/zovRyu1Yy5w1JPfJd5AcfBL7rr78+bb311nGC589+9rOYHse2uieddFLsYkY48pILhvXIj+MwgWoEVP+lV/NrNxNoNgLl3knNlkfnpz0JqM2U3p6ldKnalUDdBJD8AaBB5kdHFnsWOqM23XTThBDyyCOPxNauRx11VHSG//3f/z0xVeudd96JMGrQCadF0uVugNJEV5hy/prdTnmnHCpTpTxXcsdeggfcJUQ8/vjjIegh8CH4cZo5I1Kvv/562HOiPeHgnIerlL7tTaDeBFT/pdc7fsdnAo0kUKlNbmSajtsETMAEWpWA2sy6CSDlOg/Y8dNUHnWQWRPyy1/+Mr366qux09IHH3yQEELmn3/+tM8++8RJ6wgjhFNYdbAlzHCtNNEVd6UbogJXcu9ve5VHZSqXH/nBDTMsVG7CSeh47LHH0m9/+9uYTgXrI488Mn37299ODz/8cLrlllticTkn1xOWeAgnzoo71zFbmUCjCFAHrUzABEzABHpGQP0F6T0Lbd8m0L8EBvQ1+bxT3F1cPCT86Pjym2uuuWJx+jbbbJPoNHOeyNlnn51OPfXUmJK1xx57pJ122ilxmuwcc8xREjiUjkZHiFOdb9zK5anZH1Dlr1zeVaaisEGZFY6RDkaROI/loYceCn677757CCLrrLNOmnfeeQMb8RMPAoeYVUpTcUdA/5mACZiACXyKgNvJTyGxhQmYgAlUJKA2s88CiCKqmFIZh7zjq84wOzDxO/DAA0MYoTONMMJ5FHPPPXcslN5kk03SRhttlIYOHRod7+JX+7yD3pt8lclqv1hJSMgFDMpDeSUssKifXcVuvvnmxFkdTGtjHc2ee+6Zjj/++FjwzygHSvFhJs48HuykFHcrs1NZrLcOAde31rlXzumnCdBuWplAfxBQ3dO7uz/y4DRNoLcE+iyA1JpwuQck71TjjgBBJ1rCyAEHHJAeeOCBWDfC2pF99903klt00UVjPcPo0aPjrJEll1wyRlOKAgnxSyjRg0oEece+1vzXy5/yIV35I37yr85YsSwvvPBCYpSD7Y2vvPLKdM8990SW2EKXncVOOOGEtN566wUH4iB+/YrlxV7pqFy5XW6Wu3UTaBQB6htKeqPScbwmYAImYAImYALNQaDuAkilzmuxw5sXHzd+dJTVCaFjTid81VVXjR8jIyygZrvYW2+9tTTdiHgIu/7668dULc4gGTZsWBo8eHBM2yLOaor09Mv9KU/kQ2Z0KeWz3HXuhjkPrzikK3/y9+GHH8baGA5v/Oc//5nGjx+fxo0bl1566aVIaqGFFgqB4wc/+EEaMWJEYhF5roinUp7lT2nXei1/1k2gEQRUH6U3Ig3HaQKNIuB62yiyjrc7Aqp70rvzb3cTaAYC6iPXXQDp64Og8BoB0KJzOurzzTdfGjlyZPw4LO/5559PU6ZMSffff3+666670l//+td07LHHBl/8cxI7IwTsALXKKqukZZddNi244IKJERMO2lNaEgLK3RjlQ24SFBRW9rrGPY9P/uWP63fffTe9+eab6cUXX0xPPvlkeu655+K0+HvvvTc9+uij6b333ou8MfUMIWO//fZLa6yxRowMIXCQVh6v1sKQLm7Ki9K0bgKtREAfH1opz85r5xKgLbYygZlJQHUOXeb8vV/Obmbmz2mZQDkCqpdyq7sAoojrpRcFADon/Jiqtfjii6clllgi1juwexZq6tSpcfjhhAkT0rPPPhvCydixY9PJJ5+cOOxQnfdFFlkkDRw4MA5MpHOPPYIKwslss80WnX090FxjDzyFJy2u6fRzxgbCkK6ZLqURC0YxSJc1Gw8++GCEZ8F9rlirsfrqq0c5dt555xCclltuuciPzkfJ/SNwkK7yI8Ej92OzCbQaAT1b1GcrE2h2Ano/5PlUm5zbdbrZTBpXA6iDajfzVMrVzdzdZhPoDwKql9KbXgDJIamDr843bhohoUDYL7DAAolpWOz8lCumNiEUvPbaayGYsGgbwQGBhelOqMsuuyy9/fbbcVgi16RXTpFW0U1AZT/nnHOGn+HDhycEDH5f//rXIzqmhy222GIxTYwpVbjlghZxKD4CaIQDM2WUX6WV+8WPlQm0GgE9U+iaQthqZXB+O4OA2udy7a/b4up1QOyq+7JrdwSoZ6prtJfFjzbm3B1Bu88sAuXqotrOlhJA9MABTmY643kBMeuHHx5M9Nlnnz0NGjQofggFX/nKV8ryZ7SCKVJvvPFGuBOW0QvOLEGpc0S8SgehhxEVrrFny1umi3FN/pRXwmOXX2OHIl7ccqW8S+CQm+IgHpnlZt0EWo2A6r30Vsu/89tZBPL2W3WW9tuqPIGcV24u79u2RQLl3vHY8VN/JK9/2Of9k2J8M/O6u/tNXqupRoevlnbupnwqP7rGT9Euv5Y5jys35/Hk9jITHj/SsZdZccsdN8UnP9hVU3kcuT/Z53ZKRzpuSi/3l5vL5SOvqy0jgOSFzguIOYeFmYevnNJIQg5NYdH5MeUJ4YE1IkoTwaWo5Fa0L17jD+BKMw+XCxZ5nnM/xfi4Vp6781curO1MoBkJ5HVZ9bsZ8+k8dS4B6igqr5+Yse9J+925BF3ynhLI20WFzesb56OhVP9yN+zzusp1s6m+5q+v4WvlUUyneE08Rbvidbm0euIn95ub87Rz+9xcLu3crha/8iM9TzePq2jGv+ox9ZTrWWedNby1jABSLFSla8FRgeWPa3X4EQgAUc4P4RUHYfGjH/Yy4yZ/eRjc87iLbgpDeCnC5Pa5WX6KejFM0d3XJtAqBFTf8y8jrZJ357OzCKiuUmraYOm5vcxuoz9hJB4BK/sznwxGBSPsipxkR5BJkyaVZltwrb5HMQxuVr0jkLPMzb2LrXNDwU78WDNNPW4JAUSZLnfr5CZdfoqNnh5adD2k5fwoPHruTphKKk9b6ZTzq/hy//iTfbkw2Elgwqyw3YXBr5UJtBoB6rdGKvM6npspE/5qUcVwtYSppx/ls5iPSvb1TLvWuMiL8pfnK7evNa6iP8VXtNc16eZ+lA/cZ0b6ykclXfnJ8ygzuoRm2VWKp1Ps4aAfZa70rsVNbDFbVSagZ0Q6PvmYyldkzv6aNm1aBM7duQfmW5mpXWYOAd7l1FXVTenUT+xbQgCp9iCpQNX8CHXuJzfLvZJeLQ096NKJo1rcuT+l152dhJ9y/hSHdRNoVQKq12zcwLOjnd9kX65c1Z6xcv77y65SPivZ90c+87xUMvc2X3l8leKo5KeSfaV4ytnXIw7iJZ68PnKtdpmOIG5FVa+0i/E28zVlzssNFzoheqaV95yl7KxXJiCm4sj2/JMnT56BdeXQdmlXAvlzVMlca9kJj1JdqzVcd/7yfOV+SaclBJA80+XM9QbWkzSUtvRyYXO7cv76YpfHbbMJtBIBGiY1TjwDnInDOT35iF+zladRjXRPyyluhFOe9LW5p3HZf88IUFepo/z4iscmJWzVnqv8/uT27WqmvGzgwpQgtqFn05ann346DtP9whe+kA477LAouthJeGtXHo0uF2s/2EnTygRamUBbCCCtfAOcdxPoVAJ0RtQhoTO31lprdSqKquWGEUqCBtdFs+yKetWI7VgTATGV5+I19muuuaacS0J1yaLNDfpgwLzuO+64Iw7YPfTQQ0ulvvzyy+M513QMCx8lNH02qB2oFFF37sW6zLUUYfNr2ed6LfHn/ntq7mv8xfDdlaen+eur/2L+ivEV89ud/2L47q6L8Xfnv+jeXX66i/8zXRF8euy4mIqvTcAETKDOBGh6aKDeeuutODy02Fjl1/KbZyG3q2TO/ffFTPyoPE99iS8Pq7hzO8ykVSyX7NHl9vLLL6ctttgi/elPf0qrrLJKfIlmqobC4zdXldLL/cicl1fp4VbOrHiVrnTFVYuep1eL/2bwQ54ZAVlhhRUaUj+aoYzd5UH1Qfq3vvWtdNZZZ8UZV5yztfDCCzf1yGZ35bO7CZhA/Ql4BKT+TB2jCZhADQTU2WTb69VXX/1TIdSZ+ZSDLYKA+HC4KgrhgzOOpOSua+v1I1Bkm19jRql+1y/V5owpLzvTsLjmoF9GRkaOHBnCBzkvjnzk4ZqzZM6VCZhAIwlYAGkkXcdtAibQLQE6InRWUMV1DLjVovLOTG6uJWw1P/WMq1I6KqM6rMXrauEI884770Rn9/333w+OmhJDOMVVNCtOpanr3uhipLSIs2hXKV75lXueH8Uht2bQlV+VVfVV9uh5GZohz43MQ36PqHcsyJ84cWK6++67ow5svvnmkTyCiRZQK0wncWrkPXDcJtCqBCyAtOqdc75NoE0I0BGhI4dO56T4pbRNiln3Yqgjx0JomdUhNsO64y5FqHqqDrT0kocug+5HbteO5pwFZeZ36623xjRA3DbccMMothh1Cpd2vNcukwnUm0Dlwy3qnZLjMwETMIEKBNRBkU5Hxao2AurUiZkY1hbavnpDoMhY7IlL96M38bZiGLFA5zdmzJgoxrBhw9LQoUPDnAvJeRlhlbPL3Ww2ARNobwIWQNr7/rp0JtBSBNQhoSPjjkltt06sNPohbtJri8W++kKAeyCVm2XXCTr1b+rUqenGG28MQWSrrbZKs88+e/roo4/iWaY+shMWU7XyX5GN622RiK9NoD0JeApWe95Xl8oEWpJA3nnLzS1ZmAZnmo5azqh4nbs1OCuOvsMJaI3H/fffH2eAgGPTTTcNKqwLQVEfi3U0HAp/tforBPOlCZhAixGwANJiN8zZNQETMIFqBCx4VKNjt3oTyIWKm266KaIfOHBgaWe7e+65J0ZG5pprrliIzkgJGyZ8+OGHIZSst956ad55551BmHYdrvddcnwm0HwELIA03z1xjkzABEyg1wToEFqZwMwgIOEDgYFREKZfoTbbbLM0//zzx1QrBI333nsv/exnP0v33ntvuA8ZMiRtv/32pW2jPeoRWPxnAh1FwAJIR91uF9YETKDdCEjg8FfjdruzzV8e6hzrOhjVeOyxx9I//vGPyPSoUaNiRAOhZIMNNgjzeeedl+6777502mmnpd122620La9KKSGEa9dlUbFuAu1LwAJI+95bl8wETKCDCEgQ6aAiu6hNQgCBYdy4caU1Huuvv36Mfmj3K0Y/rrvuujR27NjYmleCCzrCC0qjKU1SJGfDBEygwQQsgDQYsKM3ARMwgUYS8NfiRtJ13N0RQIBAeND0K4SPZZddNgQLdsDaYYcd0oQJE2KEZNCgQSGYECfCSa7yemxhJCdjswm0JwFvw9ue99WlMgETMAETMIGGEpCg8PLLL8f5HwgRG2+8cWy/+/zzz6ctt9wyFpszNWuppZaKdSIILLmwQQaJJ9eL7uHoPxMwgbYi4BGQtrqdLowJmIAJmIAJzBwCnOfBSAaLy99+++1IdOedd47Rji9/+cshcEycODHWe+B3wIDyXQ4JHNJnTu6digmYQH8S8AhIf9J32iZgAibQSwKVOmuV7HuZjIOZQLcEbr311vCz2mqrpfHjx6cVV1wxTZo0KX3wwQfp9ddfj2lXWuvRbWT2YAIm0BEELIB0xG12IU3ABEzABEygvgQY/WCL3csuuyymVT388MPpoIMOSt/73vdiZOSVV15JHE6IUMxuWbnStKvczmYTMIHOIWABpHPutUtqAibQhgTckWvDm9oCRWJKFYrtd5966qlYx7H88svHblhHHXVUqQR33nlnCCDFeopQUrQrBbLBBEyg7QlYAGn7W+wCmoAJmIAJmEB9CUh4uOWWW2LHq9lnnz1GQoYNG5bmmGOOxHQsFGd/oPIpWBJePF0w0PjPBDqSgAWQjrztLrQJmEC7EHAnrl3uZPOXQ0IHOdX2uzfffHOs8dhwww0TIyAcPogAsvrqq0eB7rjjjvTmm2+GfwkeuTDS/KV2Dk3ABBpBwAJII6g6ThMwARMwARNoUwIIIgi+bLV70003hXmzzTaLKVUSMkaMGBGlZ4veRx55pESC3bKuuOKK0rUNJmACnUnAAkhn3neX2gRMoE0I6Ku0R0La5IY2cTFUxyRk3HPPPWnatGkheIwcOTJyrtGNNdZYo7T247bbbiu5/eIXv0ivvvpqXBOP6m9Y+M8ETKBjCFgA6Zhb7YKagAm0KwF1DCmfO3Ttepf7r1zUqXL1asyYMTG1arnllksrrbRSCBwSQDgNne14qZvnnntueuONN9Lxxx+fPvvZz6a99tqrtDVvXndVf9kxS7tmIaToulweRAU3+ZOAJDfrJmACzUfAAkjz3RPnyARMwARqJkAHLu8gFjt0NUdkjyZQgQB1SvWKzj1CBgLFKaecEp3+L37xi2nuuecOM274YVH6RhttFHWT7XkXWmihxLa8jIAgKEhQKSZJOmzvyw+FP8zouBWFEAkqeTj8Fv0V0/G1CZhA/xIofyxp/+bJqZuACZhAXQjQCVHHqS4ROhIT6EACxeeIZ4ppVN/61rfi+dphhx2Cip41dMLsu+++6fLLL0+zzTZbOvjgg9N3v/vdsJfwUYz3/fffT5MnT45T1dnad5tttolF7RdddFEcaLjlllvGSIuEIIQPhBPiuf3229Pjjz+eiGNk13SwlVdeOfJUTKMDb5+LbAJNSeAzXQ/nx02ZM2fKBEzABEygIgF1rP7nf/4nLbXUUonzFtZZZ53S1JaKAe1gAnUkoC6EhA4JIUrirbfeCiFlnnnmKfsxgPD8EEqmTJmSbrjhhrT33nvH9K2TTz45XX/99WnRRRdNJ510Urjfddddafjw4SGYDBgwID3zzDPpgAMOSKuuumraaqutEumg//nPf46duCSsKD/WTcAEmoOAp2A1x31wLkzABEygLgTUIaxLZI7EBP5FoFivuGYEgg6+VC58yP98882X5p133pI/2edhNCKCIL3rrruGEEE4pm4xZeuHP/xhCCYseD/jjDMiKMIHhyCyzmTnnXdOP//5z+PsEQQWRkIQfKxMwASal4CnYDXvvXHOTMAETMAETKApCBSFC661TqNcBov+JWTIHkFEZsJzfghCBVOp3n333cTC9v322y9GRxByPve5z4V/3FDorD359re/nbbbbrv00UcfRX44g+R3v/tdTMPCX54G11YmYALNQcACSHPcB+fCBEzABOpCwB2uumB0JBkBCQvSVcd0jY6SPWZNfZIf7KSKdvk1J6tzfeyxx4Z3rfNgbQiKEQ/Uqaeemjhj5NBDD420EF5In+lbKMWZ5ykc/GcCJtAUBCyANMVtcCZMwAQaQYBOEB0RFB0R/RqRVn/FqQ6Wytlf+XC67UuAOqYOfV5K1T3scjN+NeKR+5c59ys7RlMQNi644II0atSotPjii4dggV9+nLhOvNtuu23oxx13XOL0ddaHKG/oPPP4r5a+0rRuAibQfwS8BqT/2DtlEzCBBhOgE0LHhh/mch2fBmfB0ZtAyxKgQy9V7dkpunGtsEU3xZfr8suCcna/2mWXXcIZgURx/eY3v0l77rln+vznP5/YeIHfuuuuG1O39KFBoy7VhA+lladvswmYwMwn4BGQmc/cKZqACdSZAJ2KvKPDeQNsE8pZBZMmTUr3339/WmCBBdIhhxxSmhpS5yw4OhNoOwL5M9XTwvUkrISCsWPHRjLrr79+jIZwgTBxzDHHxInr6CjF/fzzz8e6kbDs+uNDA4cjLrbYYukLX/iCrGfQFXYGS1+YgAnMdAIWQGY6cidoAiZQbwJ0KujE8HvhhRfSFVdcEV9Ijz766FJSp512Wpjlzx2REhobTKDfCPA8ImSgM/1qySWXTMsvv3zkB4GC7XRZD8I20wsuuGAIJggYI7vO+sBtvfXWSzvttFOaOnVquvTSS0PwYJ2IRkP6rWBO2ARMoCoBnwNSFY8dTcAEWo0AHRkJFxyEdvrpp8dX0ieffDLOy8jdW61seX5Vjueeey46bePHj48pKe545ZRsbnYCqseMWDK9asiQIWnjjTeOxeZMx0KwOOKII2KtB1OyNJUSN57vm266KU5Z59R1DkbcfPPNS89/s5fd+TOBTibgEZBOvvsuuwm0AQF1YKQjfGhO+LPPPhtfTEeMGBGLWimuhJM2KLqLYAItT4BnlZGO++67L87uOPHEE+NZfe2110KYWGKJJUqjm/iTYptepltxeCFxDBo0KJ5t2gErEzCB5idgAaT575FzaAImUIWABAp0CSF8JX3ppZfSbbfdFiG32GKL6ORw1gCdGIWpEq2dTMAEZgIBCQzXXXddPJsrr7xyHES49NJLx/OsLBSfWY30cXihlM4S0bV1EzCB5iXgXbCa9944ZyZgAj0kQCeFaRp0asaNGxeHlREFIyAoBBMrEzCB5iHA+R08rxdffHHafffdQ/hAuOAnoQMdPxJWyL2eZezwi664mqd0zokJmEAlAh4BqUTG9iZgAi1HgE4InRV+zA1H8YV0+PDhYVZHBt3KBEygfwnwsYCpVhdeeGFsGoEgwfW8886bZp111lLm9FyXLLoMuV0ujOjZzt3zcDabgAk0BwELIM1xH5wLEzCBOhCg88Hv3XffTddff33EOHr06DTXXHPFyAgdFdzpnPCTUqeF69wsd+smYAL1J8DzOGHChBi5OPPMMyMB1oKM7NrhKlflnknZoUvYkJmwcs/jsdkETKB5CFgAaZ574ZyYgAn0gYCECjoejzzySJz/QXSs/8AuX8CKvTor7qhAw8oEZj4Bnr3NNtssfn1JPX+Gc3Nf4nRYEzCBxhKwANJYvo7dBExgJhGg48GUDuaB33LLLZHqHHPMkdZaa60wP/HEE7HLztxzzx1+8M+Uj3feeSe+oHL2wHzzzVf6mjqTsl33ZPKRnbpH7ghNoI4E9NFAOs9k/qtjUo7KBEygyQhYAGmyG+LsmIAJ9JwAHRg6LggfqBtvvDF0pnKwjSeKE9E5M+OUU05Jd999d/jndPTtttsurbrqqmmZZZaZQfhQnBG4hf7gYGUCrUBAwkYr5NV5NAETqC8BCyD15enYTMAE+oEAHRkJDE8//XS69dZbIxdM70CxPSdTsVCcqIwActhhh6Uf/ehHaeDAgWHPH3FId0c+UPjPBPqFgJ7nfknciZqACTScgAWQhiN2AiZgAjODgDosCBgIHKgvfvGLIVTgxu+0006L3wUXXJC+9rWvhR9t4ZmfD5ILNOGpBf4kMFFOKxNodQKqz61eDuffBEygPAELIOW52NYETKAFCEjoIKus/2BXHW2/O2zYsLTCCitEKdjS88c//nE6+eST0z333JPWXHPN8I9jvjhd8UmPwC3yZ8GjRW6Us2kCJmACJpAsgLgSmIAJtCwBfSWl842Q8eabb8b6D+xHjRqVWHDOlrzf//730x133JHuv//+tOKKK6aPPvqodM5ALmwQTtfSWw2OmLRavp1fEzABEzCBziFgAaRz7rVLagJtR0BCgvQHH3wwDjSjoDvuuGN65ZVX0g477BBrPliAPv/888fIR37IWbHDngshrQgMFlYmYAImYAIm0MwELIA0891x3kzABKoSKAoP2n6XqVdTp05Nq6yySnrxxRfToosuGlvust6DKVcSWCpFXoy3kj/bm4AJmIAJmIAJ9JzALD0P4hAmYAIm0DwEJExMmzYt/e1vf4uMvfzyy2mnnXZKa6+9dmzN+9JLLyVOWGaNCNOveitgkFZ3IwzdudebXDG93pat3vlyfCZgAiZgAiZQiYAFkEpkbG8CJtAyBOh0c87HP/7xjxAu6JRfffXV6eyzzy4JG2y9i0IIyVWxA48bduXsSYcfboymsPAdXQp7CwCiYd0ETMAETMAEyhOY8U1c3o9tTcAETKBpCWgb3b///e+RRwSAK6+8Mm2yySZpzjnnTGussUYIDHfddVfJXYUhrAQK7CR4YFdNkMANQYbpXOgIIqhqYcJDA/6KaZYTnBqQrKM0ARMwARMwgV4TsADSa3QOaAIm0AwEEADohI8ZMyayg8CxwQYbxFkgc8wxRxo+fHjYs/3ua6+9NsMIiEZDCK+Oey0d+g8//DBdfvnl6T/+4z9ikfu1115bGhVpBibOgwmYgAmYgAk0MwELIM18d5w3EzCBqgQ05en111+P7XfxvPnmm0cYjUogjKBYjP7YY4+FmYMK+Y0bNy6uFU9cFP4QSPJpVjhjxwL3Z555Jl1yySVpyJAhnxr9kEBTiK7ul8V0igJU3RN0hCZgAiZgAibQRwIWQPoI0MFNwAT6j4CmX02YMCF2vSInTL1CMbpB55xDBzXSwVkg2A0YMCAdf/zx6fbbbw+/2PGr1HlX+PDc9TfbbLOlwYMHxxkjQ4cODQEEN/mrFpfisG4CJmACJmACnUrAAkin3nmX2wTahABCww033BDCwyKLLFKacoWQgRuCgk5EP++882I049JLL00TJ05Mhx56aAgeCA74RXCQQrhhFKXcQnP8cOjhRRddlHbeeecIohEXLohrZqliWnkZZlYenI4JmIAJmIAJ9ISABZCe0LJfEzCBpiLAIvD3338/nXvuuSFYbLjhhmmBBRYIMx1zhIjZZ589jRgxIvL9wAMPpOWXXz4WqTMCUuysKwyeEUqInx/28iudHbdQW265Zega/YiLfviTICK9H7LgJE3ABEzABEygJgIWQGrCZE8mYALNSABh4IUXXoidrkaPHp2+/vWvz5BNCQt77rlnSZjYeuut0+mnn55YoE5nPe+wI7BIkEDAuPDCC9Nf//rX9MEHH5TiVZzYEwdrQaxMwARMwARMwARqJ/CZrpfp9DkHtYezTxMwARPodwI0X7kAUcwQ7vwQKiZPnhx+Bw0aVPRWGt0gLs4T+elPf5o23XTT9IUvfCHtt99+ad5550033XRTCDEERiBZdtll06hRo9Kf/vSniK9aPj6VYB0sVPbnnnsuLbnkkmn8+PFp3XXXjVEfCVF1SMZRmIAJmIAJmEDdCQyoe4yO0ARMwARmEoG8068OeZ507r700kuHk/xJl3/8Pv7442mjjTZK559/fggX+GFXrWOOOabUsccfu2kx8rL99tuHUMP6Dzr9eXqK17oJmIAJmIAJmMCMBDwFa0YevjIBE2ghAggI3SmEAvzpJyFBusKzLe9WW22VDjrooBA+sGc05Fe/+lU688wz06yzzpo++uij8M6idwQORhyI18KHKFo3ARMwARMwge4JWADpnpF9mIAJNCmBXIjIzXl2ERBw0y93w8zoBW7nnHNOnOux//77p3feeSfWf2y22Wbp97//fdptt91iBAQhBP+sIRk5cmRaeOGFQwAhPOlYmYAJmIAJmIAJdE/AU7C6Z2QfJmACLUBAgoayqmuEg2qK0QuECgQNDhQ85ZRTEusq5pxzznTdddfFrlkSLohrypQp6cknn0w//vGPS6MqSqtaOo12666cjU7f8ZuACZiACZhArQQsgNRKyv5MwASalkAuAMhcS4dcu169+uqr6aGHHkq77LJLYjctzhNh+14p4tI6jzFjxoQ1W/tir5/8WjcBEzABEzABE6hOwFOwqvOxqwmYQAsQyIWN3Fxr1ln/wY+w7JKF8KFRj/vuuy+98cYbpR2wLrvssrTccsullVZaKUZKlIb867rRejG94nWj03f8JmACJmACJtBbAhZAekvO4UzABFqegIQV1nJwiOFZZ52VLrnkkjRt2rT07LPPpuOOOy7de++9MR0Lv4yUXHPNNTFSwra3Dz74YCxAZyRFcbU8FBfABEzABEzABBpMwAJIgwE7ehMwgeYlgNCA8MDi8j/84Q+x3mPHHXeM8z8OPvjgNGzYsPStb30rDRjwyWzVqVOnRmGuvvrq9NZbb8Up6Iw8lBM+PCLRvPfdOTMBEzABE+hfAl4D0r/8nboJmEA/E2AROsLCyiuvnJhuxRkfTMHiEEItUP+3f/u3yOXQoUPTo48+GiegL7PMMhGukvBRzr6fi+rkTcAETMAETKApCFgAaYrb4EyYgAn0JwGEBYSQeeaZJ62xxhqlEQ1GRxA+8lEOhBAJF9KLea9kX/TnaxMwARMwARPoRAKegtWJd91lNgETmIFALmCwGB3BA8UICEoCBf74yT0csz/crEzABEzABEzABKoT8AhIdT52NQETaGMCEjwkYFBU1oNUEiTwl/stoqnmVvTraxMwARMwARPoVAIeAenUO+9ym4AJVBQmygkSEkqkG58JmIAJmIAJmEDvCFgA6R03hzIBE+gwAhJKpDdb8ZUvC0jNdmecHxMwARMwgSIBCyBFIr42ARMwARMwARMwARMwARNoGAELIA1D64hNwARMwARMwARMwARMwASKBCyAFIn42gRMwARakICmXmkqVgsW4f+3d8c4jQNRGIBnIQKJijMgEBIoPRegBYmGA4DoQYIC6OAaVJyIjkvQUVAAu/tm8yIvSyQiJZtR+CwFO7bHnnyDdvXzbEeXCRAgQOCbCAgg32SgfUwCBOZTIIPHfH46n4oAAQIE5lFAAJnHUfWZCBAgQIAAAQIECDQqIIA0OjC6RYAAAQIECBAgQGAeBQSQeRxVn4kAAQIECBAgQIBAowICSKMDo1sECBD4yv0dCwt//hmPm89j/2yTc4oECBAgQKA1AQGktRHRHwIECAwEvhIqImjEfu/v7399s7unYfk1IkCAAIFWBXqtdky/CBAg8N0FMlykQ4aKWJ8VjggeMeW2nGcbcwIECBAg0JqACkhrI6I/BAgQGAhEmMigkfPYFOszaMQlWLEtL8Xq7geSAAECBAi0KCCAtDgq+kSAAIGBQIaQDBxR8YiQkUGjO89qSDTN9SAJECBAgEBrAi7Bam1E9IcAAQIfBDJ8RKjISsfi4mLdq9fr1WpIrM9tsaHbJpc/HNZbAgQIECAwEwEBZCbsTkqAAIGvCUToyAAR8+fn5/Ly8jJs/PT0VKsdMY9XVEEiiERAWV1drduy/bCRBQIECBAgMEOBH7//c/s5w/M7NQECBAiMEOiGjwwWj4+PZWdnp7Z4e3srr6+vNZAsLy+XqIYsLS3VkHJ/f18ODw9L7JPVku7xRpzSagIECBAgMHUBAWTqxE5AgACByQhkCDk+Pi53d3e1MvLZ35DW19fLw8NDWVlZGVkBEUYmMyaOQoAAAQLjC7gJfXwzLQgQIDATgQwbFxcXtaoRl1Z1L6/KSsfNzU0NHxFYok2263a626673jIBAgQIEJi2gAAybWHHJ0CAwIQEImBEqNjY2ChHR0fD+z3i8Llta2ur7O3t1dCRN6YLGxMaAIchQIAAgYkICCATYXQQAgQITF8gKxkRKM7Ozmr1IwJJvM9Kx/X19fDSq+n3yBkIECBAgMD4AgLI+GZaECBA4L8JZOiIE0bQiFeEjs3NzXJyclKDR9x8HtP29nbZ39+vy7GfiQABAgQItCgggLQ4KvpEgACBgcBnQSLXnZ6e1kfuZvXj6uqqVj/iyVi5D0gCBAgQINCagADS2ojoDwECBDoC3QpIrI5gEa94vG5WQWK53++Xg4OD2jIrIp3DWCRAgAABAs0ICCDNDIWOECBA4F+BUZWMWB/h5Pz8vDa6vLws8V0gUf0wESBAgACBlgV8E3rLo6NvBAgQGCEQT7iKaW1trdze3pbd3d36XvWjMvhBgAABAg0L/AJXsYBWD4frLgAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A representative pharmacokinetic model is depicted in the following figure where drug is administered in a depot compartment and diffuses following concentration gradients to a central and a peripheral compartment. $A_d(t)$ represents the amount of drug in the depot compartment, $A_c(t)$ the amount of drug in the central compartment and $A_p(t)$ the amount of drug in the peripheral compartment. The rate constant $k_{cp}$ describes the diffusion of the drug from the central compartment into the peripheral compartment and $k_{pc}$ the diffusion from the peripheral into the central compartment, while $k_{el}$ is the elimination constant and $k_a$ is the infusion rate of the drug. \n",
    "\n",
    "![image-2.png](attachment:image-2.png)\n",
    "\n",
    "fig. A pharmacokinetic model. Drug is administered in a depot compartment and diffuses following concentration gradients to a central and peripheral compartment.\n",
    "\n",
    "Question : On suppose que le système suit la loi d'action de masse. Ecrire le système d'ODE et programmer le modèle correspondant.\n",
    "\n",
    "Question : simuler le système en fixant les différentes constantes.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Exercice : Reprendre et simuler les modèles SIR et SEIR de l'article [\"Why is it difficult to accurately predict the COVID-19 epidemic?\"](https://www.sciencedirect.com/science/article/pii/S2468042720300075?via%3Dihub)  de Roda et al. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation Stochastique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Installer StochPy, une bibliothèque de simulation stochastique\n",
    "\n",
    "pip install stochpy\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# installation de stochpy\n",
    "!pip install stochpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import stochpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = stochpy.SSA()\n",
    "model.Model(\"Burstmodel.psc\")\n",
    "model.DoStochSim()\n",
    "model.PlotSpeciesTimeSeries()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Créons un petit modèle et sauvons le dans Stochpy/pscmodels/test.psc\n",
    "\n",
    "model = \"\"\"\n",
    "# Reactions\n",
    "R1:\n",
    "    A > B + byp\n",
    "    k1*A\n",
    "\n",
    "R2:\n",
    "    B > A + byp2\n",
    "    k2*B\n",
    "\n",
    "#InitPar\n",
    "k1 = 0.5\n",
    "k2 = 0.5\n",
    "\n",
    "#InitVar\n",
    "A = 10000\n",
    "B = 10000\n",
    "byp = 0\n",
    "byp2 = 0\n",
    "\"\"\"\n",
    "\n",
    "fd=open(\"Stochpy/pscmodels/test.psc\",\"w\")\n",
    "fd.write(model)\n",
    "fd.close()\n",
    "\n",
    "# Question : Dessiner ce modèle et le comprendre. Puis le simuler en utilisant StochPy\n",
    "\n",
    "# Question : Coder tous les modèles la première partie et les simuler.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modèles booléens"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reprendre l'exemple de PBN donné en cours et le simuler dans les modes synchrone, puis asynchrone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Question : Construire le graphe de transition (utiliser networkx)\n",
    "\n",
    "# Question : trouver les cycles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question : Etude de l'article https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4567168/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Réseaux bayésiens"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Afin d'étudier les réseaux bayésiens, nous allons nous baser sur une bibliothèque Python puissante, pgmpy\n",
    "\n",
    "pip install pgmpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install pgmpy\n",
    "!pip install pandas\n",
    "!pip install wrapt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "L'objectif est de créer un modèle bayésien simple issu de (http://www.bnlearn.com/bnrepository/#cancer) en utilisant pgmpy.\n",
    "\n",
    "Le principe général consiste à définir la structure du réseau, puis de lui ajouter des paramètres probabilistes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Définition de la structure du réseau\n",
    "from pgmpy.models import BayesianModel\n",
    "\n",
    "cancer_model = BayesianModel([('Pollution', 'Cancer'), \n",
    "                              ('Smoker', 'Cancer'),\n",
    "                              ('Cancer', 'Xray'),\n",
    "                              ('Cancer', 'Dyspnoea')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Définition des paramètres\n",
    "from pgmpy.factors.discrete import TabularCPD\n",
    "\n",
    "cpd_poll = TabularCPD(variable='Pollution', variable_card=2,\n",
    "                      values=[[0.9], [0.1]])\n",
    "cpd_smoke = TabularCPD(variable='Smoker', variable_card=2,\n",
    "                       values=[[0.3], [0.7]])\n",
    "cpd_cancer = TabularCPD(variable='Cancer', variable_card=2,\n",
    "                        values=[[0.03, 0.05, 0.001, 0.02],\n",
    "                                [0.97, 0.95, 0.999, 0.98]],\n",
    "                        evidence=['Smoker', 'Pollution'],\n",
    "                        evidence_card=[2, 2])\n",
    "cpd_xray = TabularCPD(variable='Xray', variable_card=2,\n",
    "                      values=[[0.9, 0.2], [0.1, 0.8]],\n",
    "                      evidence=['Cancer'], evidence_card=[2])\n",
    "cpd_dysp = TabularCPD(variable='Dyspnoea', variable_card=2,\n",
    "                      values=[[0.65, 0.3], [0.35, 0.7]],\n",
    "                      evidence=['Cancer'], evidence_card=[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Association des parameters à la structure.\n",
    "cancer_model.add_cpds(cpd_poll, cpd_smoke, cpd_cancer, cpd_xray, cpd_dysp)\n",
    "\n",
    "# Test du model.\n",
    "cancer_model.check_model()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Inference des lois de probabilités (algorithme Variable Elimination)\n",
    "from pgmpy.inference import VariableElimination\n",
    "infer = VariableElimination(cancer_model)\n",
    "\n",
    "# Calcul de la probabilité de faire un test RayonX sachant qu'il y a de la pollution\n",
    "q = infer.query(variables=['Xray'], evidence={'Pollution': 1})\n",
    "print(q['Xray'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question : reprendre le modèle du cours et calculer les probabilités demandées"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tar czfv ArchiveTP.tgz *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "TestBravo",
   "language": "python",
   "name": "testbravo"
  },
  "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}