{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# S\u00e9rialisation, pickle, COVID\n", "\n", "Ce notebook aborde les sujets \u00e9voqu\u00e9s dans le titre."]}, {"cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [{"data": {"text/html": ["<div id=\"my_id_menu_nb\">run previous cell, wait for 2 seconds</div>\n", "<script>\n", "function repeat_indent_string(n){\n", "    var a = \"\" ;\n", "    for ( ; n > 0 ; --n)\n", "        a += \"    \";\n", "    return a;\n", "}\n", "// look up into all sections and builds an automated menu //\n", "var update_menu_string = function(begin, lfirst, llast, sformat, send, keep_item, begin_format, end_format) {\n", "    var anchors = document.getElementsByClassName(\"section\");\n", "    if (anchors.length == 0) {\n", "        anchors = document.getElementsByClassName(\"text_cell_render rendered_html\");\n", "    }\n", "    var i,t;\n", "    var text_menu = begin;\n", "    var text_memo = \"<pre>\\nlength:\" + anchors.length + \"\\n\";\n", "    var ind = \"\";\n", "    var memo_level = 1;\n", "    var href;\n", "    var tags = [];\n", "    var main_item = 0;\n", "    var format_open = 0;\n", "    for (i = 0; i <= llast; i++)\n", "        tags.push(\"h\" + i);\n", "\n", "    for (i = 0; i < anchors.length; i++) {\n", "        text_memo += \"**\" + anchors[i].id + \"--\\n\";\n", "\n", "        var child = null;\n", "        for(t = 0; t < tags.length; t++) {\n", "            var r = anchors[i].getElementsByTagName(tags[t]);\n", "            if (r.length > 0) {\n", "child = r[0];\n", "break;\n", "            }\n", "        }\n", "        if (child == null) {\n", "            text_memo += \"null\\n\";\n", "            continue;\n", "        }\n", "        if (anchors[i].hasAttribute(\"id\")) {\n", "            // when converted in RST\n", "            href = anchors[i].id;\n", "            text_memo += \"#1-\" + href;\n", "            // passer \u00e0 child suivant (le chercher)\n", "        }\n", "        else if (child.hasAttribute(\"id\")) {\n", "            // in a notebook\n", "            href = child.id;\n", "            text_memo += \"#2-\" + href;\n", "        }\n", "        else {\n", "            text_memo += \"#3-\" + \"*\" + \"\\n\";\n", "            continue;\n", "        }\n", "        var title = child.textContent;\n", "        var level = parseInt(child.tagName.substring(1,2));\n", "\n", "        text_memo += \"--\" + level + \"?\" + lfirst + \"--\" + title + \"\\n\";\n", "\n", "        if ((level < lfirst) || (level > llast)) {\n", "            continue ;\n", "        }\n", "        if (title.endsWith('\u00b6')) {\n", "            title = title.substring(0,title.length-1).replace(\"<\", \"&lt;\")\n", "         .replace(\">\", \"&gt;\").replace(\"&\", \"&amp;\");\n", "        }\n", "        if (title.length == 0) {\n", "            continue;\n", "        }\n", "\n", "        while (level < memo_level) {\n", "            text_menu += end_format + \"</ul>\\n\";\n", "            format_open -= 1;\n", "            memo_level -= 1;\n", "        }\n", "        if (level == lfirst) {\n", "            main_item += 1;\n", "        }\n", "        if (keep_item != -1 && main_item != keep_item + 1) {\n", "            // alert(main_item + \" - \" + level + \" - \" + keep_item);\n", "            continue;\n", "        }\n", "        while (level > memo_level) {\n", "            text_menu += \"<ul>\\n\";\n", "            memo_level += 1;\n", "        }\n", "        text_menu += repeat_indent_string(level-2);\n", "        text_menu += begin_format + sformat.replace(\"__HREF__\", href).replace(\"__TITLE__\", title);\n", "        format_open += 1;\n", "    }\n", "    while (1 < memo_level) {\n", "        text_menu += end_format + \"</ul>\\n\";\n", "        memo_level -= 1;\n", "        format_open -= 1;\n", "    }\n", "    text_menu += send;\n", "    //text_menu += \"\\n\" + text_memo;\n", "\n", "    while (format_open > 0) {\n", "        text_menu += end_format;\n", "        format_open -= 1;\n", "    }\n", "    return text_menu;\n", "};\n", "var update_menu = function() {\n", "    var sbegin = \"\";\n", "    var sformat = '<a href=\"#__HREF__\">__TITLE__</a>';\n", "    var send = \"\";\n", "    var begin_format = '<li>';\n", "    var end_format = '</li>';\n", "    var keep_item = -1;\n", "    var text_menu = update_menu_string(sbegin, 2, 4, sformat, send, keep_item,\n", "       begin_format, end_format);\n", "    var menu = document.getElementById(\"my_id_menu_nb\");\n", "    menu.innerHTML=text_menu;\n", "};\n", "window.setTimeout(update_menu,2000);\n", "            </script>"], "text/plain": ["<IPython.core.display.HTML object>"]}, "execution_count": 2, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": ["%load_ext pyensae"]}, {"cell_type": "markdown", "metadata": {}, "source": ["##  S\u00e9rialisation and pickle\n", "\n", "La s\u00e9rialisation consiste \u00e0 rassembler tous les informations qu'on souhaite conserver dans un seul bloc contig\u00fc. Parmi les formats utilis\u00e9es, on distingue les formats qu'on peut lire (json, csv, ...) et ceux qui ne peuvent \u00eatre interpr\u00e9t\u00e9s qu'avec le m\u00eame logiciel ou langage qui les a produits.\n", "\n", "La s\u00e9rialisation est plus efficace quand le fichier produit n'est pas lisible car l'information est stock\u00e9e telle qu'elle est dans la m\u00e9moire. Elle est plus facile \u00e0 restaurer.\n", "\n", "On cr\u00e9e un dataframe al\u00e9atoire et on \u00e9crire et lit sous diff\u00e9rents formats."]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [{"data": {"text/plain": ["(10000, 50)"]}, "execution_count": 4, "metadata": {}, "output_type": "execute_result"}], "source": ["import pandas\n", "import numpy\n", "\n", "mat = numpy.random.randn(10000, 50)\n", "df = pandas.DataFrame(mat)\n", "df.shape"]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [{"data": {"text/plain": ["500000"]}, "execution_count": 5, "metadata": {}, "output_type": "execute_result"}], "source": ["df.size"]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": ["for c in range(1, 51):\n", "    df[\"t%d\" % c] = [\"a\" * c for i in range(df.shape[0])]"]}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [{"data": {"text/html": ["<div>\n", "<style scoped>\n", "    .dataframe tbody tr th:only-of-type {\n", "        vertical-align: middle;\n", "    }\n", "\n", "    .dataframe tbody tr th {\n", "        vertical-align: top;\n", "    }\n", "\n", "    .dataframe thead th {\n", "        text-align: right;\n", "    }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", "  <thead>\n", "    <tr style=\"text-align: right;\">\n", "      <th></th>\n", "      <th>0</th>\n", "      <th>1</th>\n", "      <th>2</th>\n", "      <th>3</th>\n", "      <th>4</th>\n", "      <th>5</th>\n", "      <th>6</th>\n", "      <th>7</th>\n", "      <th>8</th>\n", "      <th>9</th>\n", "      <th>...</th>\n", "      <th>t41</th>\n", "      <th>t42</th>\n", "      <th>t43</th>\n", "      <th>t44</th>\n", "      <th>t45</th>\n", "      <th>t46</th>\n", "      <th>t47</th>\n", "      <th>t48</th>\n", "      <th>t49</th>\n", "      <th>t50</th>\n", "    </tr>\n", "  </thead>\n", "  <tbody>\n", "    <tr>\n", "      <th>0</th>\n", "      <td>-0.687078</td>\n", "      <td>-1.317984</td>\n", "      <td>0.485271</td>\n", "      <td>-0.274160</td>\n", "      <td>0.298888</td>\n", "      <td>0.485624</td>\n", "      <td>0.190280</td>\n", "      <td>0.938187</td>\n", "      <td>-0.934170</td>\n", "      <td>-0.200593</td>\n", "      <td>...</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa...</td>\n", "    </tr>\n", "    <tr>\n", "      <th>1</th>\n", "      <td>-1.836182</td>\n", "      <td>0.003514</td>\n", "      <td>-0.591336</td>\n", "      <td>-0.404873</td>\n", "      <td>1.002684</td>\n", "      <td>-0.320763</td>\n", "      <td>0.846543</td>\n", "      <td>0.800089</td>\n", "      <td>-1.545944</td>\n", "      <td>1.314889</td>\n", "      <td>...</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa</td>\n", "      <td>aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa...</td>\n", "    </tr>\n", "  </tbody>\n", "</table>\n", "<p>2 rows \u00d7 100 columns</p>\n", "</div>"], "text/plain": ["          0         1         2         3         4         5         6  \\\n", "0 -0.687078 -1.317984  0.485271 -0.274160  0.298888  0.485624  0.190280   \n", "1 -1.836182  0.003514 -0.591336 -0.404873  1.002684 -0.320763  0.846543   \n", "\n", "          7         8         9  ...  \\\n", "0  0.938187 -0.934170 -0.200593  ...   \n", "1  0.800089 -1.545944  1.314889  ...   \n", "\n", "                                         t41  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                          t42  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                           t43  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                            t44  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                             t45  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                              t46  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                               t47  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                                t48  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                                 t49  \\\n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa   \n", "\n", "                                                 t50  \n", "0  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa...  \n", "1  aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa...  \n", "\n", "[2 rows x 100 columns]"]}, "execution_count": 7, "metadata": {}, "output_type": "execute_result"}], "source": ["df.head(n=2)"]}, {"cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [{"data": {"text/plain": ["(10000, 100)"]}, "execution_count": 8, "metadata": {}, "output_type": "execute_result"}], "source": ["df.shape"]}, {"cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [{"data": {"text/plain": ["1000000"]}, "execution_count": 9, "metadata": {}, "output_type": "execute_result"}], "source": ["df.size"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Format csv"]}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["1.24 s \u00b1 57.7 ms per loop (mean \u00b1 std. dev. of 7 runs, 1 loop each)\n"]}], "source": ["%timeit df.to_csv('df.csv')"]}, {"cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [{"data": {"text/html": ["<pre>\n", ",0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18,t19,t20,t21,t22,t23,t24,t25,t26,t27,t28,t29,t30,t31,t32,t33,t34,t35,t36,t37,t38,t39,t40,t41,t42,t43,t44,t45,t46,t47,t48,t49,t50\n", "0,-0.6870776957633313,-1.3179835340609687,0.48527087365816,-0.27416029146654963,0.29888766592509347,0.4856236853659831,0.1902802951340063,0.938187409922834,-0.9341699028583568,-0.20059263086238005,-3.0312237407101383,-0.14792803072340885,-0.17514585644531186,0.5576225316885989,0.872080172085751,-0.6336555533330944,0.43342340019405573,0.12480280413020967,0.07634055595551763,1.1870404644728039,-0.5311587375657408,1.3739042153908319,0.6798255965813582,0.41067936759788853,0.14085251582137018,0.6429927805056169,2.1060860575958418,-0.35819378122856993,0.11921391604054045,-0.8446502736222339,-1.4067923027345046,0.6705390864871307,0.33940633716443686,1.9139672970623614,-0.4321594542966359,-0.47425010019867475,-0.5011479466303276,1.0011488927144743,-0.2162670400341441,-0.9041989514092752,-0.04800111927101552,0.5699567652022954,-1.8262540292255849,-2.394157132140209,0.016152511443535616,1.197895466553288,0.22436193047863318,0.024432777931183668,-1.249625722518517,0.49733677581926555,a,aa,aaa,aaaa,aaaaa,aaaaaa,aaaaaaa,aaaaaaaa,aaaaaaaaa,aaaaaaaaaa,aaaaaaaaaaa,aaaaaaaaaaaa,aaaaaaaaaaaaa,aaaaaaaaaaaaaa,aaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa,aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa\n", "\n", "</pre>"], "text/plain": ["<IPython.core.display.HTML object>"]}, "execution_count": 11, "metadata": {}, "output_type": "execute_result"}], "source": ["%head df.csv -n 2"]}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["318 ms \u00b1 11.8 ms per loop (mean \u00b1 std. dev. of 7 runs, 1 loop each)\n"]}], "source": ["%timeit df2 = pandas.read_csv('df.csv')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["La lecture est plus rapide que l'\u00e9criture. C'est souvent le cas avec les disques SSD, la lecture est plus rapide que l'\u00e9criture."]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Format pickle\n", "\n", "Le format pickle n'est pas \u00e0 proprement parler un format, python impl\u00e9mente un m\u00e9canisme pour s\u00e9rialiser tout objet python, cela revient \u00e0 rassembler tous les fragments de m\u00e9moires que cet objet occupe en un seul bloc qu'on \u00e9crit ensuite sur le disque, on qu'on envoie \u00e0 une autre ordinateur connect\u00e9 \u00e0 internet."]}, {"cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["132 ms \u00b1 3.25 ms per loop (mean \u00b1 std. dev. of 7 runs, 10 loops each)\n"]}], "source": ["import pickle\n", "\n", "def pickle_obj(name, obj):\n", "    with open(name, 'wb') as f:\n", "        pickle.dump(obj, f)\n", "\n", "%timeit pickle_obj(\"df.pickle\", df)"]}, {"cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["67.3 ms \u00b1 1.86 ms per loop (mean \u00b1 std. dev. of 7 runs, 10 loops each)\n"]}], "source": ["def load_pickle(name):\n", "    with open(name, 'rb') as f:\n", "        return pickle.load(f)\n", "    \n", "%timeit df = load_pickle('df.pickle')"]}, {"cell_type": "markdown", "metadata": {}, "source": ["La lecture est l'\u00e9criture sont tr\u00e8s rapides car ce qui est \u00e9crit ou lu est directement ce qu'il y a en m\u00e9moire. En contrepartie, ces informations ne peuvent \u00eatre relus que depuis python et souvent depuis la m\u00eame version de python."]}, {"cell_type": "markdown", "metadata": {}, "source": ["## Mod\u00e8le \u00e9pid\u00e9miologique (SIRD)\n", "\n", "Etrangement, il n'\u00e9tait pas possible d'acc\u00e9der aux donn\u00e9es de l'\u00e9pid\u00e9mie le jour du TD car le site [data.gouv.fr](https://www.data.gouv.fr/fr/) \u00e9tait inacessible suite \u00e0 une [panne chez OVH](https://www.lemonde.fr/pixels/article/2021/10/13/le-fournisseur-de-services-cloud-et-hebergeur-ovhcloud-touche-par-une-panne_6098154_4408996.html).\n", "\n", "* $\\frac{dS}{dt} = - \\beta \\frac{S I}{N}$\n", "* $\\frac{dI}{dt} = \\frac{\\beta S I}{N} - \\mu I - \\nu I$\n", "* $\\frac{dD}{dt} = \\nu I$\n", "* $\\frac{dR}{dt} = \\mu I$\n", "\n", "Comme les vraies donn\u00e9es sont inaccessibles, on s'int\u00e9resse \u00e0 des donn\u00e9es simul\u00e9es."]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Simulation"]}, {"cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [{"data": {"text/plain": ["(9.909090909090908, 0.8909090909090909, 0.1, 0.1)"]}, "execution_count": 15, "metadata": {}, "output_type": "execute_result"}], "source": ["def iteration(S, I, R, D, beta, mu, nu, N):\n", "    dR = I * mu\n", "    dD = I * nu\n", "    dS = - S * I * beta / N\n", "    dI = -dS - dR - dD\n", "    return S + dS, I + dI, R + dR, D + dD\n", "\n", "iteration(10, 1, 0, 0, 0.1, 0.1, 0.1, 11)"]}, {"cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[10.        ,  1.        ,  0.        ,  0.        ],\n", "       [ 9.90909091,  0.89090909,  0.1       ,  0.1       ],\n", "       [ 9.82883546,  0.79298272,  0.18909091,  0.18909091]])"]}, "execution_count": 16, "metadata": {}, "output_type": "execute_result"}], "source": ["def simulation(n_iter, S, I, R, D, beta, mu, nu):\n", "    sim = numpy.zeros((n_iter + 1, 4))\n", "    sim[0, :] = numpy.array([S, I, R, D]).ravel()\n", "    N = S + I + R + D\n", "    for i in range(n_iter):\n", "        S, I, R, D = iteration(S, I, R, D, beta, mu, nu, N)\n", "        sim[i + 1, :] = numpy.array([S, I, R, D]).ravel()\n", "    return sim\n", "\n", "simulation(2, 10, 1, 0, 0, 0.1, 0.1, 0.1)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Visualisation"]}, {"cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": ["%matplotlib inline"]}, {"cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAACMCAYAAABPlLzyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA97klEQVR4nO3deXgcxZn48W9Nz6Xbum1Jlm3ZMr7xzR2MTQgEh8OQmGXDsSHrsIRNsrk2x17JJtkk+8u5kMM5NkA2CQlgrkAIGAgJEPCBsY0v2bIsW7IlS7Ksc86u3x/VkkbyyJYtjUbH+3mefrq6unq6WqPRvKqqrlZaa4QQQgghxNByJbsCQgghhBBjkQRZQgghhBAJIEGWEEIIIUQCSJAlhBBCCJEAEmQJIYQQQiSABFlCCCGEEAngTnYF+srLy9NTp05NdjWEEEIIIc5oy5YtDVrr/Hj7zhhkKaV+DqwG6rXW85y8HOBhYCpQBXxAa30izrF3AP/ibH5Fa/3Amc43depUNm/efKZiQgghhBBJp5Q61N++gXQX/gK4uk/e54CNWutyYKOz3fekOcC/AxcAy4F/V0plD7DOQgghhBCj2hmDLK31K0BTn+zrga5WqQeAG+Ic+h7gea11k9PK9TynBmtCCCGEEGPSuY7JKtRaH3XSx4DCOGWKgcMx20ecvKSKRG1u/ckbzCnKZFHpBBaXZlOSnYJSKtlVE0IIIcQYMuiB71prrZQa1AMQlVLrgHUApaWlg63SaZ3oCONywcObDvOL16oAyM/wsdgJuBZPyWZ+cRZ+j5XQegghhBBibDvXIKtOKTVJa31UKTUJqI9TpgZYEbNdArwc78W01uuB9QBLly5N6BOr8zN8/GbdRUSiNnuOtfJW9Qm2HDrB1upmnnunDgCPpZgzKZNFTtC1uHQCxROktUsIIYQQA6e0PnNMo5SaCjwdc3fhfwONWuuvK6U+B+RorT/b55gcYAuw2MnaCizRWvcd39XL0qVLdbLuLmxoC7LVCbi2Vp9g+5FmAmEbgIIMn9PSZVq85klrlxBCCDHuKaW2aK2Xxts3kCkcfo1pkcpTSh3B3DH4deC3Sqm7gEPAB5yyS4G7tdYf1lo3KaX+E9jkvNSXzxRgJVteuo+r5k7kqrkTAQhHbfYcbWVr9Ynu5Q/vHAOc1q6irF7djEVZfmntEkIIIQQwwJas4ZTMlqyBON4a7A643jrUzPaantauwkyntctp8ZpXnIXPLa1dQgghxFg1qJYs0Vt+ho/3zJ3Ie2Jau3YfbenVzfjsTtPalel3c8OiYj6wdDLzirOSWW0hhBBCDDNpyUqA+tYAWw+d4Jkdx/jDO8cIRWzmFmWydtlkrj+/mKxUT7KrKIQQQoghcLqWLAmyEqy5I8QT22p5eNNhdh1twed2cfW8iaxdOpkLy3JxuWQMlxBCCDFaSZA1QuysOcnDmw7z+LYaWgMRSnNSef+SEm5eWsKkrJRkV08IIYQQZ0mCrBEmEI7yh53HeHjTYV6vbMSl4PKZ+axdNpmVswrxugfySEkhhBBCJJsEWSPYocZ2frf5CI9sOcKxlgC5aV7WLC5m7bLJzCjISHb1hBBCCHEaEmSNAlFb88q+4zy86TAv7K4jYmsWl05g7bLJrF5QRJpPbgQVQgghRhoJskaZhrYgG7bW8PDmw+yvbyPVa7F6wSTWLpvM4tJsmfBUCCGEGCEkyBqltNZsrW7mt5sO89T2WjpCUWYUpPOBpSWsWVxCXrov2VUUQgghxjUJssaAtmCE3283U0FsrW7G7VJcObuQtcsns2JmvrRuCSGEEEkgQdYYs7++lYc3HeaxrTU0toe4ZEYu/3n9PMry05NdNSGEEGJckSBrjApFbB7eVM03n9tLMGxz94rp3LNiOn6PPC9RCCGEGA6nC7JkQqZRzOt2cdtFU9n4qcu5Zv5Evr+xgqu/+wqv7Due7KoJIYQQ496oaMkKh8McOXKEQCCQpFoNnN/vp6SkBI9n+J9P+Or+Bv718Z1UNrSzesEk/nX1HAoz/cNeDyGEEGK8GPXdhQcPHiQjI4Pc3NwRPcBba01jYyOtra1MmzYtKXUIRqL8+E+V3PfSfryWi09fNZPbLpqKJc9IFEIIIYbcqO8uDAQCIz7AAlBKkZubm9QWN5/b4mOryvnjJ97FotIJ/MdTu7jh/lfZfqQ5aXUSQgghxqNzDrKUUucppbbFLC1KqU/0KbNCKXUypsy/DeJ853rosBop9Zyal8aDH1rOfbcuoq4lwPX3v8q/PbGTk53hZFdNCCGEGBfOOcjSWu/VWi/UWi8ElgAdwIY4Rf/cVU5r/eVzPd9I8NWvfpW5c+eyYMECFi5cyBtvvJHsKp2WUorVC4rY+KnLueOiqfzyr4dY9a0/8cS2GkZaN7EQQggx1gzVA/FWAQe01oeG6PVGnNdff52nn36arVu34vP5aGhoIBQKJbtaA5Lh9/Af183l5iUlfHHDDj7+m238dvNhmVtLCCGESKChGpN1C/DrfvZdpJR6Wyn1rFJq7hCdb9gdPXqUvLw8fD7zKJu8vDyKioqSXKuzM684i8fuuYT/vH4u24+c5Orv/plvP7+PQDia7KoJIYQQY86g7y5USnmBWmCu1rquz75MwNZatyml3gt8T2tdHuc11gHrAEpLS5ccOtS7QWz37t3Mnj0bgC899Q67alsGVee+5hRl8u/vO33819bWxqWXXkpHRwdXXnkla9eu5fLLL49bNra+I1V9a4Cv/n43T2yrZWpuKl++fh7vmpmf7GoJIYQQo0qi7y68BtjaN8AC0Fq3aK3bnPQzgEcplRen3Hqt9VKt9dL8/JH5RZ+ens6WLVtYv349+fn5rF27ll/84hfJrtY5K8jw871bFvF/H74Al1Lc/vM3ufdXW6lrGflzkQkhhBCjwVCMyfob+ukqVEpNBOq01loptRwT1DUO5mRnanFKJMuyWLFiBStWrGD+/Pk88MAD3HnnnUmrz1C4ZEYez37isu65tV7ee1zm1hJCCCGGwKBaspRSacC7gcdi8u5WSt3tbN4M7FRKvQ18H7hFj9Lb2vbu3UtFRUX39rZt25gyZUoSazR04s2tdf39f+Htw83JrpoQQggxag2qJUtr3Q7k9sn7UUz6PuC+wZxjpGhra+Mf//EfaW5uxu12M2PGDNavX5/sag2prrm1fr/jKF9+ahc3/OBVPnjBFD79nvPIShn+xwQJIYRILq01traxsXvSXUufPE1MWmtsYtKnKa/RoOmVB/1va+0sMcf3t53rz2Vx4eKk/fyGagqHMW/JkiW89tprya5GwnXNrXX5zHy+9cd9PPh6Fc/uPMY3b57PylmFya6eEEIMC1vbROwIYTtMOBo2a2fpLz9qR4noCFE7SlRHidiRnnWffV35XWXiHdvfvqiOYmubqI6ite61bdt2720noImX198xsQHSaHdp8aX8sPCHSTu/BFkirti5tT77yHbuemAzn7t6FuveVTZiZrUXQowtUTtKMBokFA31WscuA8nrtR0JErJDcYOj7u2udEzQFNXDN7WNW7lxu9xYLgtLWSatrO5tj8uDpSxcLhdu5calXGZbuboXr8tr0q6eff2tu5bufFf8/S7lwoULpVT3tkKdUiY2TymFi5h0P/tj93XtV8rZcsqg6HU80G9ZpVR3XtcxAOne5M4FKUGWOK15xVk8+g8X8+nfvc1/PbuHivo2vnrjPHxuK9lVE0IkQdgO0xnppCPcYdaRjl7pzrCzjilzuvKxwVHEjgyqbm6XG5/lw2f58Fpe/JYfr+XF6/LisTx4XB78Hj8enwe3y43H5enO797uSjv5ffd5rD7bMfvdLnevAMmtegdOfbe7ghr5x3XskiBLnFGK1+K+WxdRvjGd775QQVVDOz+6bQl56b5kV00IMUBaazojnbSEWmgJtdAaaqU11Nqdbgm10BLsyW8Nt8YNjML2wJ9/6lIuUt2ppLhTSPU4a3cqWf4sJrknkeJO6Q6E/G6z7gqSYoOlM+V1bVsu+edPjCwSZIkBUUrxiStnMqMgnU/99m2uv+9VfnbnUmZNzEx21YQYN2xtczJ4kuZgc68AqTtI6gqYgqcGUa2hViL69C1Fqe5UMrwZZPoyyfBkMME/gSJ3UdxAKTad4nHyYtIp7hR8lk9aacS4JkGWOCurFxRRmpPK3z+4mZt+8BrfvWUR754jA+KFOFdRO8qJ4AkaOxtpDDTS2NlIU6Cp13bX+kTgxGkDJY/LQ6Y3sztQyvJnMTljsgmavBkm39uTzvJmdaczvBm4XfKVIMRQkk+UOGsLSibwxEcvZd1Dm1n30Gb++epZfEQGxAvRLWyHORE4ETdQ6rtuDjbHvYvL4/KQm5JLrj+XgtQCZufM7t6e4J9Apjeze+kKkvxufxKuVgjRHwmyBig9PZ22trZkV2PEmJjl5+F1F/HpR97m68/uYV9dK/+1Zr4MiBdjntaak8GT1LTXUNtWS21bLUdaj1DbbtINnQ00B5vjHpviTiHHn0NuSi4l6SUsyFvQHTj1WqfkkuHJkH9chBjlJMgS5yzFa3Hf3yxiZkEG33lhH4caO/ixDIgXo5zWmpZQC7VttdS01VDT1hNMdQVW7eH2XsdkeDIozihmcsZklhQuOTVoctapntQkXZUQIhkkyBKDopTi41eWmwHxv9vG9fe9yk/vWMrsSTIgXoxcraHWUwKoI21HutNt4d6t1mmeNIrTiylOL2b5xOUUpRVRnGG2i9KLyPTK77sQ4lSjL8h69nNwbMfQvubE+XDN14f2NceZaxdM6hkQ/8PX+O7ahVw1d2KyqyXGsXA0TFVLFRUnKqhorqCyuZLadtM61Rpq7VU2xZ3SHUQtLVxKUXpRdwBVnF5MpjdTuu6EEGdt9AVZYsSaX5LFE/dewroHN/ORX27hM+85j3+4fLp8OYmE0lpT217L/hP7qWiuYN+JfVScqKCqpap7cku3clOaWUpJRgkL8xf2BFAZxRSnFZPly5LfUyHEkBt9QZa0OI1ohZl+Hv7IRXzmke188w972V/XxtfWzMfvkQHxYvCaA83dgdT+5v1UnKhgf/P+XmOkJqVNojy7nMtLLqc8u5wZE2ZQllWGx5KHnAshhtfoC7LEiOf3WHz/loWUF6Tz7ef3cbCxnfW3LSU/QwbEi4EJRAIcOHnAdPU5gVTFiQqOdx7vLpPpzaQ8u5z3lb2P8uxyZmbPZPqE6WR4M5JYcyGE6CFBlkgIpRQfW1VOeUE6n/zt21x/31/4yR1LmVuUleyqiRFEa01NWw27m3b3CqiqW6u7547yWT7Kssq4qOgiyieUU55tlvyUfOniE0KMaBJkDZDMkXVurpk/ick5qXz4gc3c/MPX+e4tC3mPDIgft2xtU9lcydb6rWyu28zWuq3UddQBoFCUZpZSPqGca6Zdw4wJMyjPLqc0o1SeSSeEGJUGFWQppaqAViAKRLTWS/vsV8D3gPcCHcCdWuutgzmnGH3mFWfx5L2X8PcPbeEjD5kB8feskAHx40HEjrCnaQ9b6rawpW4Lb9W/1T1RZ35KPksKl7CkcAnz8+ZTNqGMFHdKcisshBBDaChasq7QWjf0s+8aoNxZLgB+6KzFOFOQ6efhdRfy2Ue289/P7aWirpWv37RABsSPMYFIgB0NO9hat5UtdVt4+/jbdEQ6ACjNKGXF5BUmsCpYQklGiQTaQogxLdHdhdcDD2qtNfBXpdQEpdQkrfXRBJ9XjEB+j8X3blnIzMJ0/t8f91HV2MH625dQkDFMz1trPgz7nwc7CpYHXB5nbcWkPWC5Y7bdPet+9znbygXjLGhoDbWyrX4bW+tNULWzYSdhO4xCUZ5dznXTr2PJxCUsLlhMQWpBsqsrhBDDarBBlgb+qJTSwI+11uv77C8GDsdsH3HyegVZSql1wDqA0tLSQVZJjGRKKe5daWaI/6eH3+aG+15l/e1LmVecoAHx4U7Y/TRs+yVU/gnzK5tALg9YXiicA3PXwNwbILMoseccRo2djWyt39rdUrX3xF5sbeNWbubkzeGDsz/IksIlLCxYSJZPbnIQQoxvgw2yLtVa1yilCoDnlVJ7tNavnO2LOMHZeoClS5cm+FtQjARXz5tESbaZIf79P3qd76w9n6vnTRqaF9cajmw2gdXOxyDYAhNK4fJ/hvnvB38W2GGIhsGOOOuu7WhMOgzRSPyydqT/fZEgVP0Znvs8PPcFKL0I5q2BOddD+uhqzTnadpTNdZu7x1RVtVQB4Lf8LMhfwEcWfITFhYtZkLdAnssnhBB9DCrI0lrXOOt6pdQGYDkQG2TVAJNjtkucPCGYV9w1Q/wW7v7lVj591Uw+esWMcx+n03IUtv8Gtv0KGvaBJ9UENgtvhSmXgss1tBdwJg0VJsh75zF45tPw7Gdh6mUm4Jp9HaTmDG99BkBrTUVzBRurN7Lx0Eb2ntgLmAcgLypcxI3lN7K4YDFzc+fK5J5CCHEGygyXOocDlUoDXFrrVif9PPBlrfUfYspcC9yLubvwAuD7Wuvlp3vdpUuX6s2bN/fK2717N7Nnzz6neg4Vy7KYP38+kUiEadOm8dBDDzFhwoS4ZUdCfUeTQDjK5x7dzuPbalmzuJj/WjMfn3uAA+IjQdj7DLz1f3BgI2jbtBwtvBXm3AD+EfLg3rpdJtja+Rg0HTDjucpWmC7FWddCyoSkVc3WNjsadnQHVtWt1SgUCwsWsqp0FRdOupAZE2bINApCCBGHUmpL39kVugymJasQ2OC0OriBX2mt/6CUuhtAa/0j4BlMgLUfM4XD3w3ifEmVkpLCtm3bALjjjju4//77+eIXv5jcSo0Rfo/Fd9YupCzfzBB/pKmTH922hJw0b/wDtIaj20xgteN3EGiGzGK49JMmuMqdPpzVH5jCOWa54otw9G0n4NoAT9wDT3thxpUm4DrvavAlfsbysB1mS90WXjj0Ai9Vv0R9Zz1u5eaCSRdwx9w7WFm6kryUvITXQwghxrJzDrK01pXA+XHyfxST1sBHz/UcI9VFF13E9u3bk12NMaVrhvipeWl8+ndvc+MPXuXndy5jen56T6G247D9YdMdWP8OWD6Y/T4TWJWtMHcJjnRKQdFCs1z5JajZ4nQpbjAtcm4/lF8F824ya+/QjXMKRAK8Xvs6L1S/wMuHX6Yl1ILf8nNp8aWsLF3J5ZMvJ9M7Qlr+hBBiDBh1M75/481vsKdpz5C+5qycWfzz8n8eUNloNMrGjRu56667hrQOwrju/CKKJ6Sw7sHN3Hj/q/z4bxdwUWSLCawqnjMDzIuXwrXfNoFIErvZBk0pKFlqlqu+Aof/agKuXY/D7ifBkwbnXWPGcM24Etxn/+zH1lArrxx5hY3VG/lLzV/ojHSS4c1gRckKVk1ZxcVFF8sEoEIIkSCjLshKls7OThYuXEhNTQ2zZ8/m3e9+d7KrNGYtmZLN79dm8+oj32PmL18G1QLphXDhPbDwb6FgVrKrOPRcLphysVmu+QZU/cV0Ke56AnY+Ar5MM3Zr3k2m1e40g84bOxt56fBLvFD9Am8cfYOIHSEvJY/rpl/HytKVLJu4DI9LBq0LIUSinfPA90QZqQPf09PTaWtro6Ojg/e85z28//3v52Mf+1jcsiOhvqNSR5MZY7Xt/+Do22iXh02+C/jhyYuYdekNfObqubhc42uyT6JhM7/XO4+Z+b6CJyEl23STzrvZ3K3oclHbVsvG6o28cOgF3qp/C42mJL2EK6dcyarSVSzIX4BLDfPdlUIIMQ4kauD7uJSamsr3v/99brjhBu655x7cbvkRDlrTQXjlv02AFQ3BxAVwzTdR89/PIt8Eip58hx++coiDjUG+s3YhKd5RMPZqqFgeKL/SLKu/AwdehJ2PoXc+RuWOX7Ext4gXsnLYHWoCYGb2TO4+/25Wla5iZvZMeWyNEEIkkUQI52DRokUsWLCAX//619x2223Jrs7odfKICa7e+qWZ0mDJnbD4dpg4v7uIB/jKDfMoy0/nK7/fxdr1r/PT25dSkDlMj+IZSdw+Tk69mKftJja4G9nbXAHA+S1H+WRHJ6sKllJ63odhxrvN43+EEEIklfwlHqC2trZe20899VSSajIGtB6DP38LtvzCTMew5O/gsk9BZvwZ35VS3HXpNKbkpPKx37zFDfe/yk/vWMacovFxJ5ytbd44+gYbKjawsXojITvEnNw5fH7557lyypUUdLbCWw+ZKS1+fQtkTDJj1xZ9EHKmJbv6QggxbsmYrAQYbfUdNu0N8Op34c2fmLFGi/4W3vUZ88ibAXqn9iR3/WIzrYEw/3PrIlbOKkxcfZPsWPsxNuzfwBP7n6CmrYZMbyary1azpnwN5+Wcd+oB0TDsew62PmgehK1tM0h+8e0wa/U53Z0ohBDi9GRMlkiuzhPw2v/AX38EkU6Y/wG4/LPnNGno3CLzKJ67HtjEhx/YzL+unsOdF08dM2OPQtEQLx1+iQ0VG3it9jU0mgsmXcDHFn2MVVNW4bNOEyhZHpi92iwnj5hpL7Y+BI98CFJy4Py/MQHXWLw7UwghRiAJskTiBFrgrz+E1+83d8XNvRFWfB7y47TCnIXCTD+//chFfOI32/jSU7uoPN7Ov79vDm5r9N49V3GigscqHuPpyqdpDjZTmFrIugXruGHGDZRklJz9C2aVmED2sk9B5cuw9QF4cz389X6YfAEsvgPm3gDetKG+FCHEOKe1NkNB+iza7Oy9OHk9nWqn7jP7dewJeu3r2jbndV7D2a98PtzZ2Qm93tORIEsMvVC76RJ89bumFeu8a+GKL8DEeUN2ilSvmx99cAnf+MMefvxKJdVNHdx36yIy/KNn/qe2UBvPVj3LhooN7GjYgdvl5orJV7CmfA0XTbpoaJ4V6LJgxiqztB2Ht39tuhOfuAf+8DmYf7Np3SpaNPhzCTGEtNYQjaKjUYhE0JEIOhpFRyK986NRdCQKUScdjvSkY/OjUYjaYEfRsWttm3223V1G26dZ27rn2GgUbXetnX3aNnm2Nq+pzy5tjneuvd8yuieQcbZBm+N0/2V6BT+2jUY7de6nzOmCo37KjjRp77qM0vXrk3Z+CbLE0AkHYPPP4S/fhvbj5i63K74AxYsTcjqXS/H5985mal4a//r4Tm7+4ev87M6llGQP3aNohprWmq31W3ms4jGeP/Q8nZFOZkyYwWeWfobV01eT489J3MnT8+GSj8HF/wjVr5tga9uvzHs2cYEJthZ8APxZiauDSDpt2+hQCB0MYgeDPelAAB0MoUPBnn3hcPdCJNJrW4fCJvCJzQuH0ZHe24TD6HCfcrHHRSI9wVJM4EQ0muwf1amUApcL5XKBZfVe9027FEqdKe0yQx36SSvLDZ6Y8kr1HN9VVpmbg1Cx+0+zPdAyXa+LislTp57PyUOpOOVj8uPs68mnZ1/Xz7n7/PQ+NmZ/Vzru6zibnknxb6gaLjLwPQFGW30HLRKCtx6EV74FrbUw7V1wxb9A6QXDVoVX9zdw9y+34HO7+MntS1lUmrzm4XgaOht4Yv8TPL7/capaqkjzpHH11KtZU76G+XnzkzemrLPZzE+29QE4tgPcKaYbcfEdUHphzB8ukUhaaxPYdHRgd3SiOztMurOzO8/uaEd3bXcGTPlQ0ARGwSB2MNCdNvtC6ECgJ+3k61Bo6Cru8aC6Fre7J32aPDwxeW6njNsCtxtlOWnLikm7UZbVk3bH5LstlNWTjj1Oud09+S4LZblOXVuW+exZlglseq1jyrhiygrRx+kGvkuQNUCWZTF//nzC4TBut5vbb7+df/qnf8LlOnUc0Eio77CIRkz30yvfhOZqmHwhrPyiCbKSYH99Gx/6xSbqWgJ86wPns3pBUVLq0SViR/jzkT/z2P7H+PORPxPVURYXLObG8hu5aspVpHpGUIub1nB0G2x5AHY8AqFWyJtpWrfm3QSZyf1ZjjRaa3RHB9HWVqItLdgxa7u93QRF3QFRB7o7UOrolW93dKCdstj2wCtgWbh8PpSzxKaVz4vL60P5/T3prnJ+H8rbdYzXyff3pL1OGZ8P5fWapZ/gCY9Hgg4hkLsLh0RKSgrbtm0DoL6+nltvvZWWlha+9KUvJbdiyWBHYeej8PLXoemAGc9z7XfMuJ8k/tGdUZDO4x+9hHUPbubeX73FocYO7lkxfdi/CA61HGJDxQaePPAkxzuPk+vP5fa5t3PjjBuZljVC561SyryPRYvgPV+FdzaY7sQ//otZSpbDnOtg9nWQPSXZtR00rTU6ECDa0ord2tJrHW1twW5pOc2+VqKtrRCJnPE8yu/HlZJilrRUVEoqrtRUPFlZuFJTcaWm4EpNRaWk4EpNM+Xi5TvbXa+lvN5h+CkJIQZLWrIGqOvZhV0qKytZtmwZDQ0Np3yJj4T6JoRtw56n4KWvwfE9UDjPjLk6770jqlspGInyz49s5/Fttdy0uISvrZmHz53YR/EEIgGeP/Q8j1U8xua6zVjK4rLiy7ix/EYuK7ls9D6QuaHCPKR61xNwbLvJm7QQ5lxvlnOYhiMRdCRC9MQJIk1NRBoaiDY1EWlsJNrYSKSxyaxPnHCCpxYTJIXDp31N5fdjZWTgysx01hlYGZndayszA1fXOjMTKzMTV3o6rrQ0J2Dym24nIcSYlpCWLKXUZOBBoBBzv+R6rfX3+pRZATwBHHSyHtNaf/lczwlw7GtfI7h7z2Be4hS+2bOY+IUvnNUxZWVlRKNR6uvrKSwcuxNiAqYrad8f4KWvmnE7eTPh5v+FOTeYgZcjjM9t8Z21C5mWl853XtjH4aYOfnzbErLThv6//z1Ne3h036P8vvL3tIZbmZwxmY8v/jjXTb+OgtSCIT/fsMsrh3d92ixNB2H3k7DrSdj4JbMUzjOtW3OuH/L5t+yOjphAySzRpiYiDY1Em0zwFGlsINrYRLS5Of6dTR4P7pwc3Lm5WNnZeEuKnaCpn2Cpe52ByyeTtwohBmcw3YUR4FNa661KqQxgi1Lqea31rj7l/qy1Xj2I84hkiQRNC8ZffwC1b0H2NLjxxzD//WZqgBFMKcXHryxnal4qn3lkOzf+4FV+ducypuenD/q1W0OtPHvwWR6teJRdjbvwury8e+q7uan8JpYULsGlRl7gOSRypsElHzfLySOw+ynz+/Hyf8HLXzPB95zrTdA1cX6/rZt2MEi4tpZwbS2Ro0dNur6eaGMTkaZGog2NRJqa0J2dcY93ZWTgzsnBysvDN60Ma+lS3Ll5WLk5uHNyceflYuXk4s7NwZWZKeOGhBBJc85Bltb6KHDUSbcqpXYDxUDfIGtInW2LU6JUVlZiWRYFBWOgtaKvllrY/L+w5X/NVAy5M+B934eFt5pZxUeR6xcWU5KdwroHt7DmB6/xow8u4aLpuWf9Olpr3qp/i0crHuWPVX8kEA0wM3smn1/+ea4tu5Ys3zib9iCrBC78B7O0HjMB1+4n0a98C/uF/0fYmkx4wnLC3hmEO1yEa51g6uhRog0NvV/L5TIBUm4e7pwcvFOmnBIsWbl5Zp2TIy1MQohRY0gGviulpgKLgDfi7L5IKfU2UAt8Wmv9zlCcM5mOHz/O3Xffzb333jt2/kvW2syd9OZ684VpR2Hm1bD876HsihHZLThQS6bk8PhHL+FDv9jEbT97g7suncZHV84gcwATlzZ2NvLUgad4tOLR7qkX3jf9fdxUfhNzcueMnfd/gLRtEzl+nHBNrRM01fa0StVCuLYMu70D09D9GvAaygJPTjqeyVPwr7gcT3Ex7kmT8BQV4SkqxlNYYG7tF0KIMWbQQZZSKh14FPiE1rqlz+6twBStdZtS6r3A40B5nNdYB6wDKC0d+MOCh1NnZycLFy7snsLhtttu45Of/GSyqzV4oQ4zT9KbP4G6HWYiygvuhmUfNt1DY8TknFQevedivvzULtb/uZJHthzhn949k1uWTT7lcTxRO8rrR1/nsYrHeKn6JSI6wqKCRdw1/66RN/VCAtiBAKGDBwlWVhKqPEi4pqa7FSp87NgpA8atrCzcRUV4SqeQeuFFJniaNAlPXjqe9l1YR15AVb4E0X2QvgNK3wdz5kLpYrDkBmchxNg1qLsLlVIe4GngOa31twdQvgpYqrVu6K/MSL278GyMivqeqIJNPzUPEA40Q8FcuGCdeXizd2wHETtrTvLlp3fx5sEmZham8y/XzuFdM/Opbavl8f2Ps2H/Bo61HyPbl811069jTfkayiaUJbvaQy7a3EywspLggQOEKg8SrDxA6EAl4ZqankHkSuEuLDRBU1GRWYqLurfdk4qw0gfw/MNAC1T80YzhqnjePCg8NQ9mXWumhpj6LnDLtARCiNEnIZORKtNP8gDQpLX+RD9lJgJ1WmutlFoOPIJp2er3pBJkJZDWUPkSvLHe3C2oXDD7fbB8HUy5eERNw5BoWmuee+cYX31mJ0fDWygsfptWZXqyLy66mDXla7hi8hV4RtkYtL601kSOHSN4oJJQ5QFnXUmwspJoY2N3OeXz4Z02DV9ZGd7pZWZdNh3v1ClDPwYq1A77XzAB177nINQGnlSYvBymXApTL4HiJeCWsVdCiJEvUZORXgLcBuxQSm1z8r4AlAJorX8E3Az8g1IqAnQCt5wuwBIJEmiBt39jxls1VkBavrklf8nfQVZxsmuXFAdPHmRn52NQ+iQpwRO0RCYQbl7F+8qu4wsXX5KQ6R4SSYfDhA4fNq1SByq7W6WCBw+iOzq6y7mysvCVlZF+xQp8ZdPxlk3DN306nqKi4ZvTyZvWM89WOAAHXoTKl+HQq/DSV0wZtx9KlsGUS0zQVbIMPCnDUz8hhBgiMhlpAoyY+h7fB5t+Yh4CHGozrQPLP2KeTTcOWwlaQi1sPLSRDfs38Fb9W7iVmytKr2BN+RrKMxfz/RcO8Os3q0n3ufn4lTO57cIpeN0ja8C/tm3C1dUEdu0isHdfT+tUdXWvGcjdEyc6rVLT8U0vw1tmWqes3NyRPVi/o8ncgFH1Khz6i5mXTdtgeaF4qQm4plxiWr28A+imFEKIBBsTzy6cNWvWyP5ycGit2bNnT/KCLDtqumDeXG+6Bi0vzF1jugRLliSnTklU31HPi9Uv8mL1i2w6tomIjjA1cyo3ld/E6umryUvJ61V+77FWvvL7Xfy5ooFpeWl84b2zuXJ2QVJ+93Q4TLCyksCu3Sao2r2L4O492O3tpoDbjbe01LRGlXUFU9PxTps2sHFSo0HgJFT/Far+Ypajb4OOgssNRYudoOtS8zByX0ayayuEGIdGfZB18OBBMjIyyB3h/4VrrWlsbKS1tZVp04b5zryOJnjrl2Ywe/MhyCiCZR+CxXdCev7w1iXJDp48yMbqjbxU/RLbG8yjYKZmTmVl6UpWlq5kQd6C0/4eaa15ee9xvvL7XRw43s7F03P5l2vnMKcoM2F1tgMBgnv3Eti9uzuoCu7bhw6FAFCpqfjPOw//nDn458zGP2cOvunTx98z7IKtUP2GaeWqehVqt4IdAWXBpPNjgq4LIWVCsmsrhBgHRn2QFQ6HOXLkCIFAIEm1Gji/309JSQme4Zr359gO02q1/Xfmjq0pl5q5rWatHje3x2uteafxHTZWb+TF6hepPFkJwNzcuawqXcWq0lVMy5p21gF6OGrzqzeq+c4L+zjZGWbt0sl86qrzyM8YXFdrtLXVCaZ2EexaVx6EaBQw46b8c2bjnz2nO6jyTpkiz8GLJ9QOh98047mqXoWazRANAcrMOj/1UtO9OOViSM1Jdm2FEGPQqA+yRIxgm+k22f8CHNgITZXgToHz18Kyv4eJ85Jdw2ERtsNsqdvCxkMbeenwS9R11GEpi6WFS7tbrCamTRySc53sCPP9Fyt44LUq/B6Le66YzocumYbfc+agJ9LYaLr6urv8dhOuru7e7y4owD97Nv65c/DNnk3KnDm4i4pGdIvtiBbuhCObnaDrL3BkE0Scf87yZ5vAa+I8KJxrnruYXjiu7qoVQgw9CbJGM62h7p2eoOrQ62CHzS3vUy+D8nfD/JshJTvZNU24zkgnr9W8xsbqjfzpyJ9oCbXgt/xcXHQxq6as4vKSyxP6eJuDDe187ZndPL+rjuIJKXzumlmsXjAJpVT3VAmBXbsIvLPLCax2Eamv7z7eM3myaZlygir/7Nm48/JOc0YxaJEg1Gw13YuHN5nPUsuRnv2puT0BV+Fcs+TPkjsZhRADJkHWaNPRZG5rP/Ai7N8IbcdMfuE8mL4SZlxpxpyMgzsEmwPN/OnIn9hYvZHXa18nEA2Q6c1kxeQVrCxdycVFF5PiHt4vxNcqjnP/r15B7d/LZbqRS3Ujnsp9RE+cMAVcLrxl05yuvjmm22/2LKzMxI3pEmeh8wTU7TIBV91Os67fBWFnqgvlMs/r7Aq6ugKwrMnS6iWEOIUEWSNdNAI1W0xL1f4XzH/eaNM6VXaFCaqmr4TMScmu6bA42naUFw+bOwK31G0hqqMUphaysnQlq0pXsbhwMR7X8Ix509Eooaqq3i1Uu3djt7YCEHFZVGVMJDJ9JotWXUDhkvPxnXcerhRpCRlV7Kh5CkJX0NUVgJ2o6injy4wJvJzgq2C23NUoxDiXqMlIxWCcrOkJqipfNreqK5eZC2jF52HGKihaBK6xP9g5FA2xt2kvrx99nY3VG9nVuAuA6VnT+dC8D7GqdNWwPIxZh8MEDxzo1d0X2LMH3dkJmFnRfbPOI3P1td2tVJHSaTz7ajU//ctBXFWwrjSVu+d4GNsPJhqDXBbkTjfLnOt78oOtUL+7d/C1/bcQjHlMa/bU3t2NOdMhe4oEX0IIackaNuEAVL9muv/2b4Tju01+RpEJqGasgrIVY35sla1tqluq2dGwgx0NO9jZsJM9TXsI2+ahwwvyF7Bysmmxmpo1NXH1CAYJ7tvXK6AK7t2Ldh5+7EpNxedMldC1+MrKUO74/5ccburgG3/Yw9Pbj1KY6eNvL5jC8mk5LJw8YUAD5MUoojU0V/du8ap7B5oOmIlTu6TmwoQpJuDqWmdPNemsyfKsRiHGCOkuTAatoaGip7Wq6lUzxYLlM7eTz1hlugHzZ43pcR4NnQ3sbNhpgqrjO9jZuJPWkOlqS3GnMDd3LvPz5jM/fz4L8xeSnzq0c3rpSITwkSMEKw8SOlhJsGI/gd27Ce7ff+qUCTEBlXfKFJTr7Gd733Koia8/u4dNVWZ8ltdysaAki2XTclg+LYclU7LJ9I/u5yGKfoQ64PgeOHEQThwy89V1rZsPmxtWuiiX+QcrNvCKDcbSJ8I5/P4JIYafBFmJpDW01Zk/rsf39l53OA/gzS3vCaqmXALesdmZ1BnpZFfjLnY27GT78e3sbNhJbXstAC7lonxCOfPy5rEgfwHz8uYxPWs61hB1h9rt7QQPVplAqrKS0IFKQgcrCVUd6m6dArDy88zdfd0B1Vw8xUM/ZUJzR4jNVSfYVNXEm1VN7DhykoitcSmYPSmTZVNN0LVsas6g590So4AdhZba3oFX97oKWo/2Lm/5YEJp/Faw7CljvsVbiNFEgqyhoDW01MQPpgIne8r5s8x8PPnnQdFCmL7K/FEcY6J2lAMnD/QKqPY37yeqTetQUVoR8/PnMz9vPvPy5jE7ZzapnsEFl1prog0N5ll9BytN69SBAwQPHiRyNOZLyrLwTp5sntc3vQzvtDJ8ZdPwlpUl7Q6/jlCEbdXNvHGwiU1VTWytPkEgbLqWyvLSuoOu5dNyKMlOkXmyxptwAE4edgKvqt4B2IlDEGjuXd6bDukFkFZg1umFztpJd+cXjIu7kIVIJgmyzoZtw8nqUwOp43vNQ5a7pOaZrr7883qv0wvGXPef1pq6jrruLr8dDTt4p/EdOiNmQHiGN6M7mFqQt4C5eXNPeSbgWZ0vEiF0+DChgwcJHjhAqPIgoUrTQtV1Vx+YcVPesrLuZ/eZdRne0tIR/7iZUMRmZ+1JNh1s4k0n8GoJmAc8T8z0m1auaTlcMC2HGfnpuFxj63dKnKXAyd6BV8tR04LeVgftx82680T8Y/1ZThBWCGn5MQFZn8AsNW/cPCVCiKEkQVY8XbdsH9/TJ5jaZ8ZOdUmfeGoglX8epI2tSSRD0RC1bbXUtNX0WmrbajnSeoQTQfMH3OPyMCtnFvPy5pmxVHnzKc0sxaUGPn5Eh0JEjh8nXFdHpK7OWdebsVMHKwkdqoaYLj53QYFplSqbhrdsenerlLuwcMy0+Ni2Zl99K286QdebB5uobw0CMCHVw9IpJuBaNi2HuUWZeCwZryP6iASdgKveWeqgPSbdnV8PodY4L6DMYP30QvO807QC8/xH/4T+1/4s8KaNuX8shTgbEmTFajsOD91gBqVHgz35mSVxgqmZY2bsQ8SOUNdRR01rzSmBVE1rDcc7j6Pp+V1wu9wUpRVRlF5EcXoxM7NnMj9vPuflnIfXit9KpLXGbm11Aqd6Z32MSFe63gRT0cbGU45Vfj+eoiLTGjWtDO/0MtMqNW0aVsb4uxVea011U0d3wLWpqomqRjNZZqrXYnFpNgtKssjP8JGb7iMvzUtuuo/cdC/ZqV4safkSpxPq6BOA1Zm/jV3BWNe+QDMEWoDTfE+4PCbYOl0g1t8+X4YEaGLUkyArlh2Fhz9oZnTOn2WWvHLwj+7ZuG1tc7zjeNyWqJq2Go61H+seLwVmIHphaiHF6cUUpRdRkl5CcUYxRWlFlGSUkJ+S32tQuo5EiDQ2Ejl2rLvlKVLf0woVOXaMcH1995xSsazsbNyFhXgKC3EXFuIuLIhJm3xXZuaYaZVKlPqWAG9WNbHpYBNvHGxiX10rdpyPr1KQk+olN91LjhN8xQZhuWk+8tLNdk6al0y/W372on921MwL1tlsgq6+68DJ0++LndaiL2WBL92MMfOmOUtsuu/2afZ5utapcmemGFYJC7KUUlcD3wMs4Kda66/32e8DHgSWAI3AWq111eleM+ljskYIW9u0hlppCbbQHGzmZOikWQdP9spr6myitr2W2rba7rmmuuSl5FGcXtyzpEykROUwUWeQE/bjamsn2tJKtLUFu6WFaEsrdmsL0ZMtTl4r0ZYWszQ1mfFqsTwePPn5JliaWIinICaImjjRpAsKcI3w8VGjlW1rmjvDNLYFaWwP0dgWorE9SENbyOQ5241tIRragt1jvvryWIrcNCcA6w7IvOSk+chN85Lqs/C5LfweF36Phd9Jd+X5PGbttVwSrInebNt0TXY2m4ArXiAWbINQuxnzGmqPk26HcPvpg7W+PH0DsVRw+83i8fekh2Lb7ZegbpxLyIzvSikLuB94N3AE2KSUelJrvSum2F3ACa31DKXULcA3gLXnes7RSGtNa7i1d3AUPNkdNPXNOxk0S0uoBVvbWFGNN4JZwnSnJ+hUJqgUJuk0Ful0CqILyQ17yQy7SQuCvyMCre1EW1uJtmzCbtmI3d4OQIeznMKysDIycGVmYmVmYmVm4C4sxMrMwMrLM61PBYV4JppgysrOPqe5pMTQcLkUOWmmtap8AOVDEZum9p7AqycAC8UEakEO1LfR2B7svvtxoJQCn7t3IOb3WPjcXYGYhd/dk+f39JTxWi4sS+F2KSyXy1mbbbfVe9tyKdyWKefps32641wuhaUULqVQLnrSClxK4VJguZQEikPJ5TLdhf5BPrhda4gEzhCM9QnM+m5HAiaoCwdMumsJB3oPHTkXlhfcKWB5TLp77TWTznal4+0/27TLbRbLY55U0LXt8jhrK2Z/zHZ/+5VLumwTaDC3kiwH9mutKwGUUr8Brgdig6zrgf9w0o8A9ymllB5BfZS2tglEAgSiATqD7QQC7QSCbQSC7QQD7QSDHQRDHYSCHYSCnYRDnYRCHYSDAcLhIJFwgHAoQCQUJBoOEg2FiEZC2OEQ0UCAaLATT0jjjehegZIvAulhyI9apNgW/ogLnxNAecIad1jjCtm4ov190bU6S2+u9HRcmRlEM0yg5CkpwZ+RgZWViSvDBE7dQVRGBq7MLJOXkYkrLVW+YMYwr9vFxCw/E7P8ZyyrtaYjFKWpPURnOEogHCUQtgmEowQjtrMdJRCxCfbNC9sEIz3lA86+k51h6uOViUQZOX8RDFd34KVwuXrSygnEuoIy1RWcKROcuVyg6AncFEBMuisfeo5VqucY1adcVzo2v6dsT7pL73Kn5jk53fu78vq+ppPZ6zWcLGebPtu999PndeMfq3pt903HvFrvMt15KUAKkN9T575l3ZhvudSeOsajtI2bMG47aBYdwmuHcNsh3NrkeZx8jx3EctZuu/e2pcPOEsGyw1jRMFYwgkuHsXQrlh3p2a/DpoyTdukwlh3BrUOo042BS4CosrCVG1tZ2Ji1Vi60cmFzatpWFpo4+2PznbStXGhlxeyzTJ5Trvdada/h1DyNZdbKhcas6VPGxoKY7ZT8aVxyza3D+vOMNZggqxg4HLN9BLigvzJa64hS6iSQCzQM4ryDcqKhhj3XXIUrqrFsjRUFywYrCl1tMi4g1VmGinYptNeL8nlRfj+W34/lT8VKT0H5/bh8PrP2+1D+FLP2+VF+H66utT/FWftRvp6yVmYGVmYmrvT0fh/7IsTZUEqR5nOT5kv875PWmqitidh91lG7T57ZjkRjtqOnHnu647TW2BpsrZ3FpLWGqN2T11WnrnRXftTu+xqm27b3a5lbSLqO1QAx59F0va7ZEVtOO+WISfc6xgaN7eTFvL5TXjuJ3nm6O4iNPY44ZXqOi7kNRvdadR8bew6zrXtvx8QJZzrmlPK9fj96/bacoazuJ59e+v6ff0pIowE8zpIWt8wZX+OU8mcoQM/Pw4WNh8gpi1eHMWFQBDc2FlEsonhU1EmbPLezmLJd6f7z4uW7lO28no2rn3VX2myHsQjiQp9S1k0Ud5zXA7q3nbDIWWK3u87hbKuzC0C31y6DURpkDRml1DpgHUBpaWlCz+VPSadx4RSU24PL48HyeHF5vFjO4vb6cXt9uD0+3F4/Xo8fjy8Fj9eP15uC15eK2+vH5fGg3G4T1LjdKLcH5XHH5HlQXo8JiPx+lMcjrURCxKGU6e5zyyMehRBnop3/ULQdZ4k66579C1zJDXMGc/YaYHLMdomTF6/MEaWUG8jCDIDvRWu9HlgPZuD7IOp0RilpWaz+yTOJPIUQQgghEqGrX5vRMR54MLXcBJQrpaYppbzALcCTfco8CdzhpG8GXhxJ47GEEEIIIRLlnFuynDFW9wLPYaZw+LnW+h2l1JeBzVrrJ4GfAQ8ppfYDTZhATAghhBBizBtUZ6XW+hngmT55/xaTDgDvH8w5hBBCCCFGoxE347tS6jhwaBhOlUcS73JMsvF87TC+r1+uffwaz9c/nq8dxvf1D8e1T9Fa58fbMeKCrOGilNrc3wytY914vnYY39cv1z4+rx3G9/WP52uH8X39yb720TE8XwghhBBilJEgSwghhBAiAcZzkLU+2RVIovF87TC+r1+uffwaz9c/nq8dxvf1J/Xax+2YLCGEEEKIRBrPLVlCCCGEEAkzpoMspdTVSqm9Sqn9SqnPxdnvU0o97Ox/Qyk1NQnVTAil1GSl1EtKqV1KqXeUUh+PU2aFUuqkUmqbs/xbvNcajZRSVUqpHc51bY6zXymlvu+899uVUouTUc9EUEqdF/OeblNKtSilPtGnzJh575VSP1dK1Suldsbk5SilnldKVTjr7H6OvcMpU6GUuiNemZGun+v/b6XUHud3e4NSakI/x572czLS9XPt/6GUqon53X5vP8ee9vthpOvn2h+Oue4qpdS2fo4d1e879P8dN+I++7rrqfFjbMHMQn8AKAO8wNvAnD5l7gF+5KRvAR5Odr2H8PonAYuddAawL871rwCeTnZdE3T9VUDeafa/F3gWUMCFwBvJrnOCfg4WcAwzj8uYfO+BdwGLgZ0xed8EPuekPwd8I85xOUCls8520tnJvp4huv6rALeT/ka863f2nfZzMtKXfq79P4BPn+G4M34/jPQl3rX32f8t4N/G4vvuXEPc77iR9tkfyy1Zy4H9WutKrXUI+A1wfZ8y1wMPOOlHgFVKKTWMdUwYrfVRrfVWJ90K7AaKk1urEeV64EFt/BWYoJSalOxKJcAq4IDWejgm+E0KrfUrmMd2xYr9bD8A3BDn0PcAz2utm7TWJ4DngasTVc9EiXf9Wus/aq0jzuZfgZJhr9gw6Oe9H4iBfD+MaKe7dud77APAr4e1UsPoNN9xI+qzP5aDrGLgcMz2EU4NMrrLOH+QTgK5w1K7YeR0gy4C3oiz+yKl1NtKqWeVUnOHt2YJpYE/KqW2KKXWxdk/kN+PseAW+v9DO1bfe4BCrfVRJ30MKIxTZrz8DnwI02obz5k+J6PVvU5X6c/76S4a6+/9ZUCd1rqin/1j6n3v8x03oj77YznIEoBSKh14FPiE1rqlz+6tmG6k84H/AR4f5uol0qVa68XANcBHlVLvSnaFhptSygtcB/wuzu6x/N73ok3/wLi8jVop9UUgAvxfP0XG4ufkh8B0YCFwFNNtNt78DadvxRoz7/vpvuNGwmd/LAdZNcDkmO0SJy9uGaWUG8gCGoeldsNAKeXB/PL9n9b6sb77tdYtWus2J/0M4FFK5Q1zNRNCa13jrOuBDZjugVgD+f0Y7a4Btmqt6/ruGMvvvaOuq/vXWdfHKTOmfweUUncCq4G/db5sTjGAz8moo7Wu01pHtdY28BPiX9OYfe+d77I1wMP9lRkr73s/33Ej6rM/loOsTUC5Umqa8x/9LcCTfco8CXTdVXAz8GJ/f4xGG6dP/mfAbq31t/spM7FrDJpSajnm92HUB5lKqTSlVEZXGjMIeGefYk8CtyvjQuBkTBPzWNHvf7Nj9b2PEfvZvgN4Ik6Z54CrlFLZTpfSVU7eqKeUuhr4LHCd1rqjnzID+ZyMOn3GVt5I/GsayPfDaHUlsEdrfSTezrHyvp/mO25kffaTeXdAohfMHWT7MHeRfNHJ+zLmDw+AH9OVsh94EyhLdp2H8NovxTSTbge2Oct7gbuBu50y9wLvYO6s+StwcbLrPUTXXuZc09vO9XW997HXroD7nd+NHcDSZNd7iH8GaZigKSsmb0y+95hA8igQxoytuAsztnIjUAG8AOQ4ZZcCP4059kPO538/8HfJvpYhvP79mDEnXZ/9rruoi4BnnHTcz8loWvq59oecz/R2zBfupL7X7myf8v0wmpZ41+7k/6Lrcx5Tdky978519PcdN6I++zLjuxBCCCFEAozl7kIhhBBCiKSRIEsIIYQQIgEkyBJCCCGESAAJsoQQQgghEkCCLCGEEEKIBJAgSwghhBAiASTIEkIIIYRIAAmyhBBCCCES4P8Dvd4RhQK3lKsAAAAASUVORK5CYII=\n", "text/plain": ["<Figure size 720x144 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import matplotlib.pyplot as plt\n", "fig, ax = plt.subplots(1, 1, figsize=(10, 2))\n", "sim = simulation(20, 10, 0.5, 0, 0, 1.2, 0.2, 0.1)\n", "df = pandas.DataFrame(sim, columns=['S', 'I', 'R', 'D'])\n", "df.plot(ax=ax);"]}, {"cell_type": "markdown", "metadata": {}, "source": ["### Estimation des param\u00e8tres\n", "\n", "Une fois qu'on a des donn\u00e9es, on peut essayer de retrouver la param\u00e8tres de la simulation. Mais auparavant, on simule des donn\u00e9es un peu plus bruit\u00e9es en changeant les param\u00e8tres du mod\u00e8le tous les jours."]}, {"cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[10.        ,  1.        ,  0.        ,  0.        ],\n", "       [ 9.91325458,  0.93586211,  0.08147788,  0.06940543],\n", "       [ 9.83228297,  0.82054609,  0.17866666,  0.16850428]])"]}, "execution_count": 19, "metadata": {}, "output_type": "execute_result"}], "source": ["def simulation_bruitee(n_iter, S, I, R, D, beta, mu, nu):\n", "    sim = numpy.zeros((n_iter + 1, 4))\n", "    sim[0, :] = numpy.array([S, I, R, D]).ravel()\n", "    N = S + I + R + D\n", "    for i in range(n_iter):\n", "        b = numpy.random.randn(1) * beta / 5 + beta\n", "        m = numpy.random.randn(1) * mu / 5 + mu\n", "        n = numpy.random.randn(1) * nu / 5 + nu\n", "        S, I, R, D = iteration(S, I, R, D, b, m, n, N)\n", "        sim[i + 1, :] = numpy.array([S, I, R, D]).ravel()\n", "    return sim\n", "\n", "simulation_bruitee(2, 10, 1, 0, 0, 0.1, 0.1, 0.1)"]}, {"cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAACMCAYAAABPlLzyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+P0lEQVR4nO2dd3wdV5n3v2fmNklXXW6yLFuKm2JLrrFTTBLsFKdAKCmQJQmQ3WwWAgF29yWUF5a8wMIWCGxYIEAghN3QIcFpJA4hcZy4xpaL7Dix7FiWXNTb7XPeP2bu1ZUsuajdK+n5fj7nc9ozM+fotp+ec+YZpbVGEARBEARBGF6MVA9AEARBEARhPCIiSxAEQRAEYQQQkSUIgiAIgjACiMgSBEEQBEEYAURkCYIgCIIgjAAisgRBEARBEEYAV6oH0JeioiI9a9asVA9DEARBEAThjGzbtq1Raz2pv74ziiyl1MPA9cAJrfVCp60A+BUwCzgE3Ky1bunn2DuALzrVr2qtHznT9WbNmsXWrVvPZCYIgiAIgpBylFKHB+o7m+XCnwFr+7TdB6zXWs8B1jv1vhctAL4MrARWAF9WSuWf5ZgFQRAEQRDGNGcUWVrrl4DmPs03AHGv1CPAe/o59GrgOa11s+Pleo5TxZogCIIgCMK4ZLB7sqZorRuc8jFgSj8204EjSfU6py2lRGMWt/5oE/OnZbOoJI9FM3IpL/JjGCrVQxMEQRAEYRwx5I3vWmutlBrSAxCVUncBdwGUlpYOdUinpaU7gmHA77bV8fNX7WXUbK+LypJcqkryWDzDzqfl+lBKhJcgCIIgCINjsCLruFJqmta6QSk1DTjRj81R4PKkegnwYn8n01o/BDwEsHz58hF9YvWkbC+/vOsiYpbmrZOd7DzSys66VnYeaeMnGw4SiemE3aKSPBaV5LJoRh6LSvLIzXSP5NAEQRAEQRhHDFZkPQHcAXzDyR/vx+ZZ4OtJm92vAj43yOsNO6ahmDslm7lTsrlp+QwAgpEYNQ3tVNe1sfNIKzvqWnm+5njimFmFmQnBtWhGLguKc/G5zVRNQRAEQRCENOZsQjg8hu2RKlJK1WHfMfgN4NdKqTuBw8DNju1y4G6t9d9qrZuVUv8P2OKc6n6tdd8N9GmFz22ypDSfJaU9N0G2BSLsPtrGjiOt7DzSyqaDzTy+ox4Al6GYNzW71zLjnMl+XKbEeBUEQRCEiY7SekRX586Z5cuX63SPk3WsLcjOulaqnWXGnXWtdASjAGS4TSqn51JVksvK8kIunVuE1yXeLkEQBEEYjyiltmmtl/fbJyJr6FiW5lBTV2Jv1866VvbUtxOOWmR7XVy5YArXV01j1exJeFzi5RIEQRCE8YKIrBQQjlpsfKuRJ6sbeHbPMdqDUXJ8Lq5eMJXrFxVz8XmFuGVZURAEQRDGNCKyUkw4avHygZM8Wd3An/cepzMUJT/TzdqFU7m+qpiVZQWyj0sQBEEQxiAistKIYCTGS2+cZF11A8/XHKc7HKPI72HtwqlcV1nMirICTAmMKgiCIAhjAhFZaUowEuMv+06wrrqB9fuOE4xYTMr2cu1Ce0lxWWm+RKIXBEEQhDRGRNYYoDscZX3NCZ6sbuAv+08QilpMzfFxbeU0rl80jSUz8iQCvSAIgiCkGSKyxhidoSjra47zp50NvPTGScIxi+l5GVxXNY3rKqdRVZIrgksQBEEQ0gARWWOYtkCE5/ceZ111PS8faCRqaUoLMhOCa0FxjgguQRAEQUgRIrLGCa3dYf685zh/qq5n41tNxCxN+aQs7r70PN67dLqEhBAEQRCEUUZE1jikuSvMM7uP8b+bD7P7aDulBZl8YvVs3rtkuoSDEARBEIRRQkTWOEZrzfqaEzyw/g12H21nZmEmn1g9h/csLhaxJQiCIAgjjIisCYDWmudrTvDA82+wp76dWY7YukHEliAIgiCMGCKyJhBaa/689zgPPH+AmoZ2yoqy+OSa2bx70XQJcioIgiAIw8yYF1mRSIS6ujqCwWCKRnX2+Hw+SkpKcLvdKR2HZcXF1hvsO9ZBeVEWn1wzh3ctKhaxJQiCIAjDxJgXWbW1tWRnZ1NYWJjW4Qq01jQ1NdHR0UFZWVmqhwPExdYxHnj+gC22JmVx75o5XF8lYksQBEEQhsrpRNaY2KwTDAbTXmABKKUoLCxMK4+bYSjWLpzGU598B//9N0txGwb3/nIHVz/wEk/srCdmpZfIFgRBEITxwqBFllJqnlJqR1JqV0p9qo/N5UqptiSbLw3heoM9dFRJ13EahuLaymk8fe87+N6tSzEUfPKx11n7wEv8aWc9logtQRAEQRhWBi2ytNb7tdaLtdaLgWVAN/CHfkxfjttpre8f7PXSga997WssWLCAqqoqFi9ezKZNm1I9pHPGMBTXVU3jmXsv5cFblwDwicde5+oHXmJdtYgtQRAEQRguXMN0njXAW1rrw8N0vrTj1VdfZd26dWzfvh2v10tjYyPhcDjVwxo0hqG4vqqYaxZO48ldDXzn+Te4539fZ+6UA9y7Zi7XLJyKIXu2BEEQBGHQDNeerA8Ajw3Qd5FSaqdS6mml1IJhut6o09DQQFFREV6vF4CioiKKi4tTPKqhYxqKdy8q5s+fvozvfGAxMUvz8f/dzrXffZmndzWIZ0sQBEEQBsmQ7y5USnmAemCB1vp4n74cwNJadyqlrgW+o7We08857gLuAigtLV12+HBvh1hNTQ0VFRUAfOVPe9hb3z6kMffl/OIcvvyu0+u/zs5OVq1aRXd3N1dccQW33HILl112Wb+2yeMda8Qszbrqer6z/gAHT3Yxf2o2n7piDledL54tQRAEQejLSN9deA2wva/AAtBat2utO53yU4BbKVXUj91DWuvlWuvlkyZNGoYhDT9+v59t27bx0EMPMWnSJG655RZ+9rOfpXpYw45pKG5YPJ3nPn0Z375lEaGoxd2/2M7f/nwrHcFIqocnCIIgCGOG4diT9UEGWCpUSk0FjmuttVJqBbaoaxrKxc7kcRpJTNPk8ssv5/LLL6eyspJHHnmED3/4wykbz0hiGor3LinhXVXFPPraYb76ZA3v+++N/PiO5cwszEr18ARBEAQh7RmSJ0splQVcCfw+qe1updTdTvVGYLdSaifwXeADOt2in54l+/fv58CBA4n6jh07mDlzZgpHNDq4TIOPXFLGo3eu4GRniBu+9wob32pM9bAEQRAEIe0ZkidLa90FFPZp+0FS+UHgwaFcI13o7OzkE5/4BK2trbhcLmbPns1DDz2U6mGNGhefV8TjH7+EOx/Zyu0/2cy/vHsBH7pw/ItMQRAEQRgswxXCYdyzbNkyNm7cmOphpJSZhVn8/mMXc+9jr/PFP+5m/7EOvvSu83GbY+LBAYIgCIIwqsivo3BO5Pjc/PiOC7jr0nIefe0wdzy8mdbusRsvTBAEQRBGChFZwjljGorPX1vBf9y0iK2HWnjP917hzRMdqR6WIAiCIKQVslwoDJobl5VQVpTJ3z+6jfd+byPfvXUJ75w3OdXDEgRBGPdorYnpGDEdI2pFiVpRu245dR0lZiX196079lHLbo/3x3NLW8R0LHEdS1uJFNMxNDphZ2H1KluWk2u7XXPqOeJJa52w1VonzoGm3/Z4WaP7be/Vh8XSyUv54oVfTNnrJCJLGBLLZhbw+D2r+LtHtnLnz7bw+WsruHNVWdo+KFsQBOFsiQuMcCxMxIoQsSKEY+FEPWyFicR6t4ctO49a0V71uG3yseGY09bHJnG+PjYRK9JLIKUThjLshIFpmImyYRiYykShMJWJYdjtSjl1ZZfjbfFycnuirAwUKlHu1Z5sm2QzKSO1sTdFZAlDZnpeBr/9h4v4x1/v5KtP1rDvWAdfe+9CvC4z1UMTBGGcYWmLQDRAd6SbrkgX3VE7D0QDdj2pvTvSTXe0m0A00EuoxAVR1Ir2CBmrRwAlCx/N8EUdMpWJ23DjNt14DA8e005uw22XDQ9u002GK6Nfm/ixLuXCNMxeucuwy6YycRlOXZm9+/vW+5wjnseFUUI4qf7ryUJJ6B8RWcKwkOlx8b1bl/LA+gN8d/0Bahu7+OFtyyjye1M9NEEQUoylLboj3XSEO2gPt9MZ6aQj3EFHuIOuSNcpoihZLCXEU5KYOls8hocsdxY+l6+XUIkLF5/LlxA2fYVOv+1OOdFu9BZJieOSrpF8vGnIP54TDRFZZ4nf76ezszPVw0hrDEPxmSvnMneKn3/6zU5uePAVfnT7cs4vzkn10ARBGAIxK0ZnpNMWSOEegdRXMLWH2+kId5zS1hXpwtLWaa+hUGS6M8lyZZHpzrSTK5PJmZPJdGUm2rLcWWS67DzDlWHXk46Lt2e6M3Eb7lH6CwlC/4jIEoad66uKmVmQxd/9fCvv//5Gvn3LYtYunJrqYQmCAASiAVqDrTSHmmkNttISajk1D7XSEmxJCKauSNcZz+t3+/F7/GR7ssl2ZzM1cyqz82bbdU82OZ4c/G5/77rHnxBNGa4M2cspjDtEZAkjQmVJLk/ccwl3PbqNu3+xjX+8ci73rJ4tX6KCMIxEYhFbECUJpJZgyynCKdkmGAv2ey5DGeR58xKpNLuUHK8tjHI8OQlx5Pf0qbv9+N1+WQoThH4YeyLr6fvg2K7hPefUSrjmG8N7ToHJOT5+edeFfO73u/jP597gjROd/PuNVfjc8mUsCAMRs2I0B5s5GThJY6CRk91OHujJmwPNtIZa6YwMvIUh25NNvjefPF8ekzMnMzd/Lvm+fPK8eYm8wFeQqGd7smUDsyAMM2NPZAljCp/b5Fs3L2LulGz+7dl9HG7q4qHbljM115fqoQnCqBKMBmkMNCaEUlw89RJQ3SdpCbX0u38p15vLpIxJFGUUUTKppJdAiudxUZXrzZX9SIKQBow9kSUepzGHUop/uPw85kz2c+8vX+fdD27goduXs3hGXqqHJgiDJr4ZvC3URmuoNZE3BZpO8Tw1djfSETn1qQimMin0FVKUWcTkzMksKFxAUUaRLaYyixKiqiijCI/pScEsBUEYCmNPZAljlivOn8LvP3YJdz6yhZt/+Cr/fmMVNyyenuphCROc/sRSW7jNzuMpbPe1h9oTdh3hjgFjKHlNb0Iszc6bzYXTLuwlmCZl2uV8b77sZRKEcYyILGFUmTc1myfuWcXdv9jGvb/cwf5jHfzTVfMwjHGyIT7cDTVPgMcPc64El8QJSyWBaIA3Wt5gX9M+attre0RUknhqD7WfNuBktjubXG9uIpVkl5Drsct53rxEe44nhzxvHoUZhfjdfrnJQxAEEVlni8TIGj4Ksjz84s6VfPmJ3fz3i29x4EQn375lMX7vGH47dp6AzT+CLT+GQLPd5suDhe+DqltgxkqQH90RpS3Uxr7mfexr3sfepr3sa97HofZDif1Nma5M8n35CXFU4i8hx5vTSyjlefPI8eQkytmebFzGGH5fCoKQUob07aGUOgR0ADEgqrVe3qdfAd8BrgW6gQ9rrbcP5ZrC+MDjMvj6eyuZNyWb+9ft5cbvb+RHty9nRkFmqod2bpyogVcfhOpfQywC866Fiz4GkSBU/xJ2PAZbH4b8WbbYqroFCs9L9ajHNFprTnSfYF/zPmqaa6hpqmFf8z7qu+oTNpMzJ1NRUMGVM6+korCCioIKpmVNE++SIAijynD8i/ZOrXXjAH3XAHOctBL4vpMLAkopPnxJGedN9vPx/9nOdd99mS9efz43LStJ7x9DreHgX+DV78Gbz4MrA5bcBhd+DIpm99jNuQJCHVDzJ6j+Ffz13+Cv34SSC2yxteB9kFWYunmMASxtcaTjCDXNNexrskXVvuZ9NAebEzYzc2ZSOamSm+fdTEVBBfML51PgK0jhqAVBEGyU1oN/+KXjyVo+kMhSSv0QeFFr/ZhT3w9crrVuGOicy5cv11u3bu3VVlNTQ0VFxaDHOdqMtfGmA4cau/g/v61m86Fm3jGniK+/tzL9vFrRMOz+rS2uju+GrMmw8i5YfidknsWPens97PoN7PwVnNgDhgvmXGULrrlrwT2xw1pErAgHWw8mhFRNUw37W/Ynoo27lIvz8s6jorCC+QXzqSioYF7BPLLcWSkeuSAIExml1La+K3mJviGKrFqgBdDAD7XWD/XpXwd8Q2u9wamvBz6rtd7ax+4u4C6A0tLSZYcPH+51nbEmWsbaeNMFy9L8z+a3+cZTNWjgn6+ex+0XzcJM9ab47mbY9lPY9BB0HoPJ58NFH4fKmwa/sf3YLtj5S9j1W/uc3lxY8B5Y9AGYcSEY4z8oZFOgiS3Ht7D12FaqT1bzZuubRKwIABmuDObmz6WioCIhqmbnzZYwBoIgpB2nE1lDXS5cpbU+qpSaDDynlNqntX7pXE/iiLOHwPZkDXFMwhjFMBS3XTiT1fMn84U/7OIrf9rLn3bW8833VzFnSvboD6jpLXjt+7DjfyDSDeethvf8t50PdTlzaqWdrrwfav9qe7d2/Ra2PwJ5pVB5sy24iuYMz1zSgLZQG1uPbWXzsc1sPraZN1vfBOwN6ZWTKvlQxYeYXzCf+YXzmZk9U0IbCIIw5hmSyNJaH3XyE0qpPwArgGSRdRSYkVQvcdoEYUCm52Xw0w9fwOM76vnKn/Zw3Xc38InVs7n78vNwmyPs4dEa3n7N3sy+70l7Sa/qZttzNWXB8F/PMG3Rdt5qCH8LatbZ+7c2fAte/g8oXmqLrQXvA/+k4b/+CNIR7mDb8W1sPraZLce2sL95PxqNz/SxZPISriu/jgumXsD5hedLdHJBEMYlg14uVEplAYbWusMpPwfcr7V+JsnmOuAe7LsLVwLf1VqvON1503VPlmmaVFZWEo1GKSsr49FHHyUvL69f23QY73ihsTOU8GjNn5rNv91YRVVJ3vBfKBa141u9+iAc3QYZ+fZeqxV/B9lTh/96Z6LjmO3Zqv6lvbSoTJh9BSy6xb6D0Z0x+mM6A12RLrYf386WY1vYfGwzNc01WNrCY3hYPHkxF0y9gBVTV1BZVInbFFElCML4YET2ZCmlyoE/OFUX8L9a668ppe4G0Fr/wAnh8CCwFjuEw0f67sfqS7qKLL/fn4iVdccddzB37ly+8IUv9GubDuMdbzy39zhf/OMuTnaE+Lt3lPOpK+aS4RmG5aRQB2z/Obz2A2h7GwrK7bsEF98KnjTZUH18ry22qn8DHfXgyYbzb7DHOPPilMXfCkQD7DixIyGqdjfuJqZjuAwXVUVVCVG1aPIivKYEZRUEYXwyInuytNYHgUX9tP8gqayBjw/2GunKRRddRHV1daqHMaG48vwprCgr4BtP1/DDlw7y7J5jfOP9VVxYPsgQCG11sOkHsO0RCLVD6cX2czHnrrWX8NKJKefbe7fWfBkObbCXE/f+EXb8AqYvg1Wfsb1bI7xZPhQLUX2y2t5T1bCZ6sZqolYUU5ksKFrARxZ+hAumXsDiSYvJdKfZnaGCIAgpYMyFMv7m5m+yr3nfsJ5zfsF8Prvis2dlG4vFWL9+PXfeeeewjkE4M7kZbv71fVW8q6qY+36/iw889Bq3rizlvmvmk+M7i+Uny4KG1+3N7Hv+YO+/Ov8GuPgeW6ykO4YJ5ZfZ6dr/sL1br3wHfvU3UDQXLvmUc8fj8NyBZ2mL6pPVvNbwGluObWHnyZ2EYiEMZTC/YD63VdzGBVMvYOmUpRJGQRAEoR/GnMhKFYFAgMWLF3P06FEqKiq48sorUz2kCcvFs4t49lOX8q3n9vOTDbW8UHOCr713IWsqpvQYBdvg+B4n7XbyvRDpspfbVt4NK//evpNvLOLJhOUfhSW3216tDQ/A4x+Dv3zdFo1Lbx/UcqfWmt2Nu3n60NM8e+hZTnSfAGBe/jxumnsTK6auYNnUZeR4coZ3PoIgCOOQIcXJGgnSfU9Wd3c3V199NTfddBOf/OQn+7VNh/FOFHYcbuLB3zyDt3kf757azDvzT+JprLH3V8Xx5dnhEqYssPOKd4EvN2VjHhG0hjfX23clHn7F3ri/8m5YcdcZA6VqrXmj5Q2eOfQMT9c+zdHOo7gMF6umr2LtrLVcUnwJeb680ZmHIAjCGGPEgpGOBOkusgBef/113vOe9/DWW2/hcp3qDEyH8Y5LupuTvFJOfqIGokEAotrgkJqOr6SK6fOWo6YstIVVTvHEejjz25vglQdg/1PgzoJlH7ZDUORO72VW21bLM7XP8MyhZzjYdhBTmayctpK1s9ayunQ1ud5xJkQFQRBGgJEMRjohWbJkCVVVVTz22GPcdtttqR7O+CMWgcYDfZb69th31sXJLIKpC+GCv7WF1JSF1Opi/vkP+9nxZitrXJP5atVCpuWmX6iDEad0JZQ+ZgvQDQ/YG/w3PwRVt3B0yQd5pn0/zxx6hn3N+1Aolk1Zxq3zb+WKmVdQmCHPUhQEQRguxJM1Aoy18aYFbz5vx4U6vhtO7odY2G433DBpvrPUtzAhqPBP7vc0MUvzs42H+I9n92Mais9dO58PXlCKkepH86SQEw2v8+dX/pWnm6up9to3CFTllLN23o1cNfMqpmRNOcMZBEEYLrTWEIuhY7FErqNRsCx0NAbWWfRpQFugNdqyTq1bGtD2cVr3X9c66Rh96vksC0g6n2Pfq18nny/pXOdyrMbpdxL6DO19+kj+G5xq75s3n8I7Pzqir6l4soT0pb0enrkP9j5ue6eKF8N5a2whNWWB/ViZcwhcaRqKO1eVcWXFFD73h2q+8IfdPLGjnm+8v4qyoolzB1xLsIXnDj/HM4eeYeuxrWg08ybP516yuPrABmbUvgjtFninQ/nkibWcKowJdDSKDofR4TBWKIyOhBN1HQrZ7Yl6734rFEKHIz32kYgtUKIxtBWDmNWTx6J96o7AOaVu9RZH8bplQTTaq19bMYj2FkuJ3LJS/adNDUrZyTBQ8XKfpJLtzqIdBUoZp7aRdKwrtYGPRWQJqSEWhS0/ghe+ClYU1nwJLvrEsIUfKC3M5Bd3ruTXW4/w1SdrWPvAS3zmyrncuaoM10g/midFdIQ7eOHtF3j60NO8Vv8aMR1jVs4s7l50N2tnraU8r9w2DHXAtp/Bq9+DR98L0xbDqk/bNwSkW4wwYdTQWqO7u4l1dBBrb8fq7EJHIk4KJ5V7EpGILXSS6r1swqcec0qKC6NwkjAKhYZPjLhcKLcbZRh22TDANJ26iTJMMA2UYaJcpv0ZcOqYBsq0j1EuF8rrATPpHKbRp550jHMuu83Jnesl+pJt+usz+xxvGHauDFtMGHGBkVQ3DEChjB5RAwoMZYsbo0eU9NT7ns8+pzL62MfLyULJOV4p7LJj0++1JiCyXDgCjLXxjjp122Ddp+BYNcy+Eq79dygoG7HLHW8P8n//uJs/7z1O5fRcvv+hpZTkj49gmd2Rbl6qe4mna5/m5aMvE7EiFGcVs7ZsLdeUXcO8/HkDf7lFQ7DTibXV/BYUnAeX3Gs/K9ElEdrHGtqysLq7sdrbiXV0JPJYeztWRyexjnas9o6k3LHp7EzYEosNfgBKoTweW9CcKXnc4HajXHbZ8HhQHq99fDx5PU573z43hrePrcebdB4PKrl/hIP0CoLcXTjKjLXxjhqBVnjh/8GWn9jPA1z7DTsY6Cj8h6O15undx7jvd9VkeV08eudKZk/2j/h1R4KoFWVj/UbWvbWOF+teJBANMCljElfPupq1ZWupKqo6t/8arRjU/MkO/9CwE7Kn2XcjLvsweLNHbB4TAR2NokOhnmUtx0tjOctdiWWuRFvEWQoL9V72ctos5xgr0N1bKHV0YHV2ntH7ozIzMbOzMXOyMbJzMLOzMXLieTZmdg5Gth8zJwcjy4/yxkWTp5dA6lc4meIFFSYmIrJGmbE23hFHa9j9O3jmc9DdCCv+Ht75efCNfkDLvfXt3P7wJrSGRz66goXTx06YggMtB3j8zcdZd3AdTcEm8rx5XDXzKtaWrWXp5KWYQ13q0xoO/gU2fBtqX7Jjia24CxZ9EArPG55JjDG01uhAgFhLC9HWVmItrcRaWuzU2kqstYVovNzSSqy1FR0IJETVkDxDcUwT5fX28tIYPh9Gbg5mdk5CMBnZ/l713nk2pt+PcsuDuQVhuBGRNcqMtfGOKE1vwZOfgYMvQvFSuP7b9ub2FFLb2MWHfryJ9kCEhz9yARfMOn2wzlTSGmzlqdqnePytx9nbtBeXcnFpyaXcMPsG3jH9HbjP4aaAc6JuG7zybahZB2j7sT1z19rPSJyxYkzu3TpnweSUdSjU/wmVwszJwczPx8zLS+RGRkbPcpXXYy9tuXuWsAxv8nKW1657BmjzeFD9xOITBCF9EJE1DJimSWVlJZFIBJfLxe23386nP/1pjH7W+9NhvCknErQDYr78LXt/z5ov2Y+BSZMf5/rWAB/6ySbqWwP84EPLuHxe/yEhUkF8OfCPb/6RF4+8SMSKML9gPjecdwPXll9LgW8URWHr27D/GTuw6aENYEUgsxDmXA3zroHzVoM3dcuuOhq1RVNzM9HGRmLNzUSbmog1NRFtaiba1EisyWlrbj43wZSfj5mfh5mXhyteT+7LyZElMkEQRGQNB8kR30+cOMGtt97KJZdcwle+8pVTbNNhvCnlrb/Ak/9ob6ZeeCNc/TV7D1aa0dgZ4vafbObAiQ4euGUJ11VNS+l4+i4H5nvzua78Om6YfQPzC+andGwABNvhrfWw/2l441kItoLpgbJLbcE195pTosoPBqu7m2hzsyOUmnqJplhTo503NxFtbCLW2urExOmD242roABXYSFmYSGuggLMwkLMfBFMgiAMLyKyhoFkkQVw8OBBLrjgAhobG0/ZZJwO400JHcfhz1+AXb+BgnK47j9tT0ca0xaIcOfPtrD97Ra+8b4qbr5gxqheP2XLgUMlFoUjr9mCa/9T0HzQbp9aZS8pzrsGpi0CpdBaY3V12wKpsZFoY5PtYYqXGxt7BFVzM7q7u99LGn4/ZmEBrsIiXIUFjngq7NNm50ZOzoS9ZVwQhNFlRIKRKqVmAD8HpmDHWn1Ia/2dPjaXA48DtU7T77XW9w/2mgDHvv51QjX7hnKKU/BWzGfq5z9/TseUl5cTi8U4ceIEU6ZM8IjZVgy2/RSevx+iAbjsPjvuktuX6pGdkdwMN4/euZK//8U2/s/vqmkPRvjbd5SP6DUHWg787AWfHf3lwMFiurAmLyVqziQ66f1Ea3cR27uB6OZqon/6PtHgD4lFM4hGMoh2xdCh8KnnUAozPx9XURFmYQEZMxb3eJ4K+4inggIMX/q/nwRBEJIZyo7KKPCPWuvtSqlsYJtS6jmt9d4+di9rra8fwnWEdKZhJ6z7NBzdZi8bXfctO0r7GCLDY/Kj25fxqV/u4KtP1tAejPLpK+YMuyekv+XAW+bdkj7LgUnoaJTQgQME9+whcvy47W1yvE7RJjsfyONk5hXj8nsws4JkWE24JoVxZZmYMytwzb8YV+UVuGbMxszPl03dgiCMawb9Dae1bgAanHKHUqoGmA70FVnDyrl6nEaKgwcPYpomkyenz4bpUSXUAX/5V9j0fXsj9Pt+BJU3jdnHs3hdJv/1wSV87ve7+O76A7QHInzp+vOH/MzDsbAcqLUm2tBAoLqawM5qAtXVBPfsQQeDCRszNxdzUhGuwiIyFi7ENanIXporKsJVVOh4o4pwFeT3DhMQCcKhl+0lxf1PQ+1GOPQtmLGy527Fojlj9n0jCIJwOobl30il1CxgCbCpn+6LlFI7gXrgn7TWe4bjmqnk5MmT3H333dxzzz0Tb9+H1lDzBDx9H3Q0wPKP2HcOZuSnemRDxmUafPP9VWT73Dz8Si0dwSjffH/lOT+Gp7/lwIqCCu5bcR/XlF2T8uXAWGcnwd27E4IqUL2T2MlGAJTHg6+igrybbyKjahEZlQtxT5uG8gzycUduH8y50k7XfQsadvTcrfj8l+1UcB7MuQrKL4OZl6QkfpogCMJIMGSRpZTyA78DPqW1bu/TvR2YqbXuVEpdC/wROGUtSSl1F3AXQGlp6VCHNCIEAgEWL16cCOFw22238ZnPfCbVwxpdWg7BU/8MB/4MUyrh5p/DjAtSPaphxTAU//f6CnIz3Hz7+TfoDEX47geX4HWd/s6zk90n2Vi/kVeOvsLGho20hdoSy4Hvmf0e5hXMG6UZ9Ca+7BfYWU1gVzXB6mpCb76VuCPPM3MmWRddZAuqRVX45s0bvKA6E0pB8RI7vfNz0HoE3njG9nBt+6ntFVUmTF9qLz2XXWZ7vMbA3j5BEIT+GNLdhUopN7AOeFZr/a2zsD8ELNdaNw5kk653F54LY228ZyQahlf/C/767/aDRFd/wY7abo7v/TQPb6jl/nV7ececIn542zIyPT3zjcQi7Di5g1eOvsIr9a+wr9m+GaPQV8gl0y9hTemaUV8O1FoTPXasl4cquGcvOhAAwMzLw7eoiowqJ1VWYubljdr4TkskCHWb4eBf7WjzR7eBjoHLZwut8sug7HI7kG2axFoTBEGAkbu7UAE/AWoGElhKqanAca21VkqtAAygabDXFEYRy7J/6Patg71/tL1YFe+ynzeYW5Lq0Y0KH11VRrbPxWd/V82HfryJr99UQnXzJjYc3cCmhk10R7txKReLJy/m3qX3smr6Kubmz8VQo/NA2lhnF8HduwhU77IF1c5qoidPAqDcbrznV5B34422oFpUhXvGjPRd3nb7HO/VpXY92A6HN0LtX23htf5+4H7w5sKsVY7ougwmzZP9XIIgpC1DcUVcAtwG7FJK7XDaPg+UAmitfwDcCPyDUioKBIAP6HQLzCX0EA3ZXoR9T9p7ZjqPg+Gy98ms/SbMW5vqEY4qwWiQadMOc/Wlm3ipbgM3PWULmOKsYq4vv55Lpl/Ciqkr8HtGJ+J55PhxAtu20b11G93btxPavz+x7OeeWUrmhRcmBJV3/nyMkVr2Gw18Ofb7Lf6e6zxpC6646Nr/pN3un2oLs/LL7DwvPbcbCIIwMRnK3YUbgNP+C6m1fhB4cLDXEEaBYBsceM72WB14HsId4M6COVfA/OvtDcvjYFP72aC1pratlg1HN7CxfiNbj28lFAvhNb0snLyIXQcuotCs4sfvup4ZBVkjPpbwwYN0b91GYPs2urdtJ1JXB4DKzCRz8SKyP/YxMhYvwrdwIa78cf4a+SdB5Y12AtuzWvuSLbgO/gV2/dpuzy/r8XKVXQZZhSkbsiAIwpjZVKO1Tt+ljiTGhKOu7ajtqdr3ZM/z6LImwcL32cKq7NIJs9m4M9zJpoZNbKjfwCtHX6GhqwGAstwybpp7E6umr2LZlGX4XD62HW7mwz/dwi0/fI1H/3Yl500aPg+WDocJ7t1L97btdG/fTmDbNvuRMYBZWEjm0qUU3PYhMpYuw1cxX+JL5c+y09LbbW/eiZoeL9eu38G2n9l2Uyp7RNfMi8CbncJBC4Iw0RgTj9Wpra0lOzubwsLCtBZaWmuampro6OigrKws1cPpQWs4uc/2Vu17Cuq32+0F50HF9TDvOihZPiE2FFvaYn/zfl6pf4UNRzew88ROojpKljuLC6ddyCXTL+GS4kso9hf3e/ye+jbueHgzWsPP71zBguLcQY0j1tlFYMcO20u1dRuB6upEXCr3zFIyly0nc9lSMpYuxTNrVlq/79OOWBTqX4faF23RdWQzxJwHQ+dMtx/5VFDm5E7KL0vpg64FQRi7jPlnF0YiEerq6ggmBUdMV3w+HyUlJbjdKQ4yacWgbosjrJ7sebbc9OUw/zo7Fc2dEJuGQ7EQmxo28cLbL/DikRdpCtr3XlQUVCRE1aLJi3AbZ/eaHTzZyYd+vImOUJSffvgCls86c9yr6MmTjpdqG4Gt2wju22ffXGAY+ObPJ2P5MjKXLiNz2VJckyYNZbpCXyIBOLIJjmyxPwfx1HWit51/So/g6ivEMvJSMnRBENKfMS+yhLMkErD/c9+3zo491N0Ihtte/pt/nR1dO2daqkc5KnSEO3i57mXWv72eDUc30B3tJsudxarpq7i05FIuLr6YooyiQZ//aGuA2368iYa2ID+8bRmXzu0RRlprwocOEdi+3dmkvo3I4bcBUD4fGYsWOV6qZWQsXozpH9n9XcIAhDqguba38IrXO+p722YU9PZ8JafMggnxz4ogCP0jIms8091sBwfd9yS8uR4iXeDNsTesz78OZl8BvsEtaY01Tnaf5C9H/sILb7/ApmObiFpRCn2FvLP0nayesZqV01biMYfvjruTHSFuf3gzJ4408J+Vbha0HyWwexfB6l3EWloAOzZVxrJlZC5dSubyZfgqKkYu2KcwfIS77c31vQSYI8LajgBJ35ve3N5er7xSyJsBuaV2uJMJsr9RECYqIrLGC1rbX/BHNjvLH5vg2G47aGP2NNtTNf86mPUOcE2MH/LD7YdZ//Z6Xnj7BapPVqPRzMiewZrSNawpXUNlUSXmMO41i3V2Edy7h+CuXQR27aZ7ZzWxBtvroZXCN/s8fJVVtrdq+TI8ZWUoY3TiZgmjRDQELYf7EWAHofVt+/OYjH8K5M5whNcMW4Ql1+UxQoIwphGRNVaJhuFYdY+gOrLZfl4g2GEWSpbBjAvtB+0WL4EJ8GOutWZv096EsHqr7S3A3l+1pnQNq0tXMztv9rBsFNfhMMH9byQCfgZ37+r1SBr39On4qioxKxbwYL2b33f4+ef3LuWjq9LopgdhdIlFoL3e/meo9YiTH+4pt9VBLNz7GF+u7fVKiLAkMZZXaj+AXZYjBSFtGZGI78II0NXYW1DVvw5RZ7N/Xqkd6XrGSpixAiYvGPePtYkTsSJsP749IayOdx/HUAbLpizjpnk3sXrGaqb5h7bXTFsW4dpaArvs5b7A7t2EamrQkQgAZkEBvsqFZF+9lozKhfgqK3EV9Gx4/1I0RtNjO7h/3V7agxHuXTNH7giciJhuyJ9pp/6wLHvDfesRaHs7SYgdcWJ/vWzHqkvGlXGqAMuf5WzQL7Pj2Ml7TRDSEvFkpQrLssMqxAXVkU3QbHtlMNz2M9rigqpkxYTZsB4nEA2w8ehGXjhi3xHYHm7Ha3q5uPhiVpeu5rKSy8j3DS4AZ+IZf453KlC9i+CePVidnQAYmZn4FizAV1VJRmUlvoWVuKcXn1E0RWMW9/1+F7/dVseKWQVcMruIC8sLWFyad8YHTAsCYHtJg629xVfr270FWXefJ5P5cnsEV+LuSKfsnzohPNyCkEpkuTAdCLbbzwKMC6q6LRBqt/syi6D0QltQzVgJ0xZPyM2yrcFW/lr3V9a/vZ5X618lGAuS48nhspLLWFO6houKLyLTnXlO54y2tBCurU2k0JtvEdi9m1ij84xytxvfvHn4KheSUVlFRuVCPOXlKHNwosiyND986SDrquvZ29CO1uB1GSwpzePC8kIuLC9k8Yw8fG4RXcIgCXfZwit+J2RLbU+5754wl8/2evUSX2V2Oa/U9rwJgjAkRGSNNuEue2PssV09nqoTe0BbgILJ5/cIqhkr7C/ACeDuj1pRTnSfoKGrgfrOeo51HaO+q56GrgYaOhs43H6YmI4xJXMKq0tXs7p0NcumLDtj/CodiRA+Ukf4UC3hgwcJ1dYSPmiLqnjUdLAfmuyZNRPf+T1eqpF8xl9bd4Qth5p57WATr9U2sbe+HUuDx2WwZEYeK8sLubC8gKWl+SK6hOEhFrW9Xb3EV21PORrosVWmvfx4ihes3BZmnnP7h0YQJioisoYbraHjmL2HouWQ/QUWLzfX9g5y6Mm2o6knlv6Wj9uQCt2R7oSAauhq6EmdDdR31XOi+wSWtnodU+ArYGrWVIqziinPK2f1jNWcX3h+v0tzyV6p0MGDhGsP2fUjRyAaTdiZRUV4y8rwOMlbbufu6dMH7aEaDtoCEbY6omtTbTO7j7bZoss0WDwjjwvLC1hZXsjS0nwyPCK6hGEm8b1V2xOOIrkcbO1t759iP24rq8j2tifnvdoKZV+YMKERkTUYIgHb9R4XTr0E1eHe/xGi7Hg4+bOcTa9ldnnSfJhcMS4eV6O1pinYRENnQy8BlSyo2kJtvY5xKRdTsqYwLWuanfzTKM4qTpSnZk0lw5XR+zpxr1TtQVtMncEr5ZlVhqe8HE/ZrISwMnPGxi3x7UFbdG06aAuvXY7ocpuKxTPyWFlmLy8unZlHpmdi3OQgpJDu5h6PV/x7rqvRDmrc1WjvBYtvceiL4bLFVmaR/VDurEm9RVhfsebLk71iwrhBRFZ/aG1/cfTniWo5dGrEZ3eWs7ehrOfhtHExlTcDXN6RH/MI0xHu4GjnUeo66uzUaaejHUep76wnbPW+9TzLnZUQUMX+4oRHKl6elDEpEaNKR6NEm5uJNTURbWwi2tTYu9zYRKS+nnBd3cBeqfKyRDnVXqmRoCMYYevhFnt58aDt6YpZGrepqCpxPF1lhSybmU+WV0SXkAKioVOFV1cjdJ102pp6+roaoc8/XgmU6YgvR4RlFthR9ePlzMKker6de3PEWyakJSKykulqhJ/fYAupcGfvvuzigYVUVtGY/4BHrAjHuo5R11HXI6Y6ewRVX09UjieHkuwSpvunM90/vZegmuafhh8fsZYWoo1JgqmxkVhToyOemhLlWGtrIr5UMsrnw1VYiFlUiHvKVGeJbxbe8nI8s2aNGa/USNAZitqerlrb01VdZ4sul6GoKsll8Yx8Jud4KfJ7KfR7mOTkhVlePC7xEghpQDRsC7Fk4dXdN2+2bQLNdrlvMNc4hssRXgU9S5SnCLLkeoF4zIRRQURWMrEo/Po2yJvZW0zlzRzzd/RprWkLtfUSTsn5sa5jxJK+wFzKRbG/mBmZ05llTqbUKKRY5THZ8lMY8+ELWsQ6OrE62ok2t/QWT42NxNr6/y9VZWbiKixMiCdXYRGuoiJcRYWYhYV2ubAQs7AIIytT4kmdJV2hKNscT9em2mb21LcRjFj92uZmuCn0eyjyexPiqygpt5NdFq+YkDZYlr0k2d0EgRZHoCWLsHi9uXfdivR/PmXYYsyXC+5McGfYd1zGy4mU6bQ7ZXeSjSujT3uyfaZ9h6Z8h01oRkxkKaXWAt8BTODHWutv9On3Aj8HlgFNwC1a60OnO2fa7MlKAyKxCB2RDrrCXXYe6aIj3JN3RjppC7Rw8uTbNDceob35GHR1kxnSZAYhMwRFsQwmWVkUxnzkRjz4wwYZIQt3IILqDGB1dqK7u884FiMryxZMRZNsARUXTIWniicjU+5KGi26QlEaO0M0doZp7AzR5OTx8smkclug/x+iDLeZJL56RFih30Oh34vfa5LhdpHpMcnwmGS4zUTZ5zIxDPmBEVKI1vbDvhOiq+VUURZsswM7R7rt/bbJKRovdzt3gJ8jyuwj2LJ6hJgns09bBniS+uNCrr+2xPGZEmojzRmRiO9KKRP4HnAlUAdsUUo9obXem2R2J9CitZ6tlPoA8E3glsFec6wQsSIEogG6wl10RjrtFLbzuEjq7G4l1NFKuL2NaEc70c4OW/B0daO6AhjdQTzBKBlhyAhBppNnhDQZYSgIYfeFzzSaLpQ3ipFtYfoNjOxszFw/xvRsjGw/pj/bbsv2Y/izMXKyMbOzMfxOW3Y2pt8vDzVOU7K8LrK8LmYWZp3RNhy1aO6yRdjJJEHWlCTSjrYG2VnXRnNXmJh1dv+AZbh7xFeGxxFgvcouMjwGmR7XgDZel4lpKNymwmUYuJJyt2Fgmgq3oXCZ8T67320q8YROdJSyn//oy7FXJQaL1vZjkRJCrNsRZk45kiTSkoVZJG7T5ZS7egRc54kku6TzniuGO8mz5rVFl+mxcyOpnMidsuEeZLtzXsNl37hluJy2eD1edupmn3rf/gn8GR3KOsEK4E2t9UEApdQvgRuAZJF1A/AvTvm3wINKKaXTZI1Sa03EitAV6qQ70E53VyuBQDuBQAfBQAeh7g5CgU5CwS7CwS4iwW4iwW6ioQCxUJBoKIgOhbDCIaxQGB0Jo8IRjEgMVwy8EZJEkiYjBHmOYPJEz2J8ShHL9KIzfZCVgZGVhTHVj8ufgzs7F29OHu7s3B6BlO0/RSAZ2dkjFgdKGFt4XAZTc31MzT3zsrhlaVoDEZo6Q3SFY3SHowQjMbrDdjq1HD2lvbkrTF1LjEA4RiASP8cgPAVnwDSULdAcEeY27XpchNnizUjYKaUwFBhO3ruuUEl9dj2p34jbJ/f3tQeInxsUPdeBHhsFGIbCMU+y6ymj7P7kY1TSueK/XQqVVO6/3bZ3rhc/T3J70vHxSrJtf9fp2xfvSD5vv7aqx/5sfn/7Cum+h/Q9hzrFov/r9HfpHrtMJ/WxdNlJ9b4xesDz9T63Am1hxEKYsQBmNIAZC2BEg4m6KxbAiAYwY0HMaDdmLOjUA5jRbgwrgmFFUFYEpaMYsQhGJIKygqhEXzRhY1hRlBVOqkdQjO5PsKVcaMNEKxdamWjDZSdlAAZaGWhlgjITZa0MMEy0039KnzLRhoFmgL5425QFFK/9x1GdbzJDEVnTgSNJ9Tpg5UA2WuuoUqoNKAQah3DdIdFyso6a667GFdW4ohp3FFxJ3/tuJw1mu7VWYLlMLI8H3C60xw3+TFRuFqbfjyvbFkeenDy8OXm2F8nvx/Db/Ybfj5HVU1cZGfJfupASDENRkOWhIGt4BbplaYJRW4T1iK8YoUiMmKWJWJpozCJqaaIxTdSyiMQ0MSdP9Dl2dp8mYlm2fdKx8baYpYk47VprLA2W1mgnt5y2mGUl+iyNY6uxrFPtdZJdcj1maTS2U0TreNnOLaeP+DGOXbyMBk3PteP9wkRB0VvYjRwGFm6iSSlm56qnbGLhImYnFetVN4nhwsIkhpsYpupTT+5XveuupGRiYSiNgYXpJOOUXGMQwySCqaxTbA10v8eZykKhOXSibcyKrGFDKXUXcBdAaWnpiF4rIyuX5qpSlMeD6fVieH24vD5MbwYubwZuXyYeXyaejCy8Pj/eDD++zBy8vixcXh/K60V5PL2S4eS4XCKKBOE0GIYi0+OSuF/ngE4ScQkBh+4lwOLl5PZkgZewOY1dT7nHJtE4QJ/ubZIYa+/xn3pMsv0p8z1l/qe3OOV6p5yxf7Gq+7Hs124Ix54rZzpHf9cd6jn7v865XmMQ4zrnI3pjOelMTPGl9rtmKFc/CsxIqpc4bf3Z1CmlXEAu9gb4XmitHwIeAnvj+xDGdEZ8mdlc/+OnR/ISgiAIw0Z8Kc8442KUIAjpxlACiGwB5iilypRSHuADwBN9bJ4A7nDKNwIvpMt+LEEQBEEQhJFk0J4sZ4/VPcCz2CEcHtZa71FK3Q9s1Vo/AfwEeFQp9SbQjC3EBEEQBEEQxj1DWqzUWj8FPNWn7UtJ5SBw01CuIQiCIAiCMBZJu4jvSqmTwOFRuFQRKbzLMcVM5LnDxJ6/zH3iMpHnP5HnDhN7/qMx95la60n9daSdyBotlFJbB4rQOt6ZyHOHiT1/mfvEnDtM7PlP5LnDxJ5/qucuT84UBEEQBEEYAURkCYIgCIIgjAATWWQ9lOoBpJCJPHeY2POXuU9cJvL8J/LcYWLPP6Vzn7B7sgRBEARBEEaSiezJEgRBEARBGDHGtchSSq1VSu1XSr2plLqvn36vUupXTv8mpdSsFAxzRFBKzVBK/UUptVcptUcpdW8/NpcrpdqUUjuc9KX+zjUWUUodUkrtcua1tZ9+pZT6rvPaVyullqZinCOBUmpe0mu6QynVrpT6VB+bcfPaK6UeVkqdUErtTmorUEo9p5Q64OT5Axx7h2NzQCl1R3826c4A8/93pdQ+5739B6VU3gDHnvZzku4MMPd/UUodTXpvXzvAsaf9fUh3Bpj7r5LmfUgptWOAY8f06w4D/8al3WfffqDn+EvYUejfAsoBD7ATOL+PzceAHzjlDwC/SvW4h3H+04ClTjkbeKOf+V8OrEv1WEdo/oeAotP0Xws8DSjgQmBTqsc8Qn8HEziGHcdlXL72wKXAUmB3Utu/Afc55fuAb/ZzXAFw0MnznXJ+quczTPO/CnA55W/2N3+n77Sfk3RPA8z9X4B/OsNxZ/x9SPfU39z79P8n8KXx+Lo7c+j3Ny7dPvvj2ZO1AnhTa31Qax0Gfgnc0MfmBuARp/xbYI1Salw8hVVr3aC13u6UO4AaYHpqR5VW3AD8XNu8BuQppaalelAjwBrgLa31aAT4TQla65ewH9uVTPJn+xHgPf0cejXwnNa6WWvdAjwHrB2pcY4U/c1fa/1nrXXUqb4GlIz6wEaBAV77s+Fsfh/SmtPN3fkduxl4bFQHNYqc5jcurT7741lkTQeOJNXrOFVkJGycL6Q2oHBURjeKOMugS4BN/XRfpJTaqZR6Wim1YHRHNqJo4M9KqW1Kqbv66T+b98d44AMM/EU7Xl97gCla6wanfAyY0o/NRHkPfBTba9sfZ/qcjFXucZZKHx5guWi8v/bvAI5rrQ8M0D+uXvc+v3Fp9dkfzyJLAJRSfuB3wKe01u19urdjLyMtAv4L+OMoD28kWaW1XgpcA3xcKXVpqgc02iilPMC7gd/00z2eX/teaHt9YELeRq2U+gIQBf5nAJPx+Dn5PnAesBhowF42m2h8kNN7scbN636637h0+OyPZ5F1FJiRVC9x2vq1UUq5gFygaVRGNwoopdzYb77/0Vr/vm+/1rpda93plJ8C3EqpolEe5oigtT7q5CeAP2AvDyRzNu+Psc41wHat9fG+HeP5tXc4Hl/+dfIT/diM6/eAUurDwPXA3zg/NqdwFp+TMYfW+rjWOqa1toAf0f+cxu1r7/yWvQ/41UA24+V1H+A3Lq0+++NZZG0B5iilypz/6D8APNHH5gkgflfBjcALA30ZjTWcNfmfADVa628NYDM1vgdNKbUC+/0w5kWmUipLKZUdL2NvAt7dx+wJ4HZlcyHQluRiHi8M+N/seH3tk0j+bN8BPN6PzbPAVUqpfGdJ6SqnbcyjlFoL/B/g3Vrr7gFszuZzMubos7fyvfQ/p7P5fRirXAHs01rX9dc5Xl730/zGpddnP5V3B4x0wr6D7A3su0i+4LTdj/3FA+DDXkp5E9gMlKd6zMM491XYbtJqYIeTrgXuBu52bO4B9mDfWfMacHGqxz1Mcy935rTTmV/8tU+euwK+57w3dgHLUz3uYf4bZGGLptyktnH52mMLyQYggr234k7svZXrgQPA80CBY7sc+HHSsR91Pv9vAh9J9VyGcf5vYu85iX/243dRFwNPOeV+PydjKQ0w90edz3Q19g/utL5zd+qn/D6MpdTf3J32n8U/50m24+p1d+Yx0G9cWn32JeK7IAiCIAjCCDCelwsFQRAEQRBShogsQRAEQRCEEUBEliAIgiAIwgggIksQBEEQBGEEEJElCIIgCIIwAojIEgRBEARBGAFEZAmCIAiCIIwAIrIEQRAEQRBGgP8PXviajnQa+vMAAAAASUVORK5CYII=\n", "text/plain": ["<Figure size 720x144 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["fig, ax = plt.subplots(1, 1, figsize=(10, 2))\n", "sim_bruit = simulation_bruitee(20, 10, 0.5, 0, 0, 1.2, 0.2, 0.1)\n", "df = pandas.DataFrame(sim_bruit, columns=['S', 'I', 'R', 'D'])\n", "df.plot(ax=ax);"]}, {"cell_type": "markdown", "metadata": {}, "source": ["#### Premi\u00e8re id\u00e9e\n", "\n", "On calcule les param\u00e8tres chaque jour \u00e0 partir de la d\u00e9finition du mod\u00e8le."]}, {"cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [{"data": {"text/html": ["<div>\n", "<style scoped>\n", "    .dataframe tbody tr th:only-of-type {\n", "        vertical-align: middle;\n", "    }\n", "\n", "    .dataframe tbody tr th {\n", "        vertical-align: top;\n", "    }\n", "\n", "    .dataframe thead th {\n", "        text-align: right;\n", "    }\n", "</style>\n", "<table border=\"1\" class=\"dataframe\">\n", "  <thead>\n", "    <tr style=\"text-align: right;\">\n", "      <th></th>\n", "      <th>beta</th>\n", "      <th>mu</th>\n", "      <th>nu</th>\n", "    </tr>\n", "  </thead>\n", "  <tbody>\n", "    <tr>\n", "      <th>0</th>\n", "      <td>1.129296</td>\n", "      <td>0.241642</td>\n", "      <td>0.145303</td>\n", "    </tr>\n", "  </tbody>\n", "</table>\n", "</div>"], "text/plain": ["       beta        mu        nu\n", "0  1.129296  0.241642  0.145303"]}, "execution_count": 21, "metadata": {}, "output_type": "execute_result"}], "source": ["def estimation_coefficient(sim):\n", "    N = sim.sum(axis=1)\n", "    diff = sim[1:] - sim[:-1]\n", "    mu = diff[:, 2] / sim[:-1, 1]\n", "    nu = diff[:, 3] / sim[:-1, 1]\n", "    beta = - diff[:, 0] / (sim[:-1, 1] * sim[:-1, 0]) * N[:-1]\n", "    return dict(beta=beta, mu=mu, nu=nu)\n", "\n", "df = pandas.DataFrame(estimation_coefficient(sim_bruit))\n", "df.head(n=1)"]}, {"cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAACMCAYAAABYtiw7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA4zElEQVR4nO3dd3wc1b338c/ZJq16tWUVW+6W5W7Z2JgQHEINNYBxqAZCCS25hHtTbi4Q8oQkpJKQAAYcA6GaalqAhBbAGORe5IqbZFnVqrurbef5Y3alVZetlVbl9/ZrXlN39+x4NPudM2dmlNYaIYQQQghxfEyRLoAQQgghxGAmYUoIIYQQohckTAkhhBBC9IKEKSGEEEKIXpAwJYQQQgjRCxKmhBBCCCF6wRKpD05LS9O5ubmR+nghhBBCiB5bt25dpdY6vaN5EQtTubm5FBYWRurjhRBCCCF6TCl1oLN5cppPCCGEEKIXJEwJIYQQQvSChCkxrLy6oYTlH++NdDGEEEIMIRFrMyVEfyutdfKjlzbT5PWTHGPjkoKcSBdJCCHEECA1U2LY+P27u9AaZo9O4mevbmX74bpIF0kIIcQQIGFKDAvbD9fx0vpili3K5dGrCkiKsfK9p9dR6/REumhCCCEGOQlTYsjTWnPfW0Uk2q3ccsoE0uKi+Otlcyg56uTOVZvQWke6iEIIIQYxCVNiyPtoVwWf7Knktm9MJDHGCkBBbgo/PTuP97aX8cjHX0W4hEIIIQYzCVNiSPP5Nb96awdjUmO4csGYVvOuWZTLt2aM4v5/7mDN3qoIlVAIIcRgJ2FKDGkvrjvEzrJ6fnTmFGyW1pu7UorfXDSDsWmx3PbsesrqXBEqpRBCiMGs2zCllFqhlCpXSm3tZP4pSqlapdTGQHdX+IspxLFzuL38/t1dzBmdxFnTMjpcJi7KwsNXzMXh9nHL0+vx+Pz9XEohhBCDXU9qplYCZ3azzH+01rMC3b29L5YQvffox/sor2/if7+Vh1Kq0+UmjoznV9+eTuGBo/z67R39WEIhhBBDQbdhSmv9MVDdD2URImzK61088vFezpqWwdwxKd0uf/6sLJadmMvjn+zjrS2l/VBCIYQQQ0W42kwtVEptUkq9rZTKD9N7CnHc/vjebtxePz86c0qPX/PTs/OYPTqJ/161ib0VDX1YOiGEEENJOMLUemCM1nom8Bfg1c4WVErdoJQqVEoVVlRUhOGjhWhvd1k9z395kCsWjCE3LbbHr7NZTPzt8jlEWc3c9NQ6Gpu8fVhKIYQQQ0Wvw5TWuk5r3RAYfguwKqXSOll2uda6QGtdkJ6e3tuPFqJDv3p7B7FRFm4/deIxv3ZUop0/L53NnooGfvLyFrmhpxBCiG71OkwppTJUoHWvUmp+4D3lpj0iIj7bU8n7O8q5ZfEEUmJtx/UeJ01M44enTWL1psM89fmBMJdQCCHEUGPpbgGl1LPAKUCaUqoYuBuwAmitHwYuBr6nlPICTmCplsN5EQF+v+aXbxWRlWRn2Ym5vXqvm0+ZwPqDNfzije1Mz0pk9ujk8BRSCCHEkKMilXsKCgp0YWFhRD5bDE0vry/mjhc28adLZ3HB7Kxev1+tw8O3/vIffH7NG7edRGpcVBhKOby5PD7qnB5qnB5qnR5qHYF+SFcX7Ls8TBoZzyUFOczMTuzy9hZCCNHXlFLrtNYFHc6TMCWGApfHxzd+9yGpcVG8dssiTKbw/PBuLanl2w99xgljU1h5zXzMYXrfwczl8bUOQF0EorZdk7frm6LGR1lIsFtJtFuJi7KwuaQGl8fPxBFxLCnI4YLZWaTHS6gVQvS/rsJUt6f5hBgMVny6j8O1Ln6/ZFbYghTAtKxEfnF+Pj96aQsP/GsXd5w+OWzvPdjsrWjg7te28cmeyi6XCw1EiXYr49PjSLRbSYqxtpretouPtmAxt27GWefy8ObmUl4oPMQv3yriN//cwSmTR7CkIJvFU0ZgNcsTsQaC3WX15KbFyv+HGLYkTIl2tNa8s62Mv36whxtOHse5MzMjXaQuVTU08dAHe/lm3ggWjk8N+/tfOm80hfuP8uf39zBrdBLfmDIy7J8xkLk8Ph58fw+PfLyXaKuZWxdPICMxuseBqDcSoq18Z/5ovjN/NHvK61lVWMzLG0r4V1EZaXE2LpiVxSUFOUzOiA/bZ4qe8/s197+zk4c/2stpU0fyt8vnSKASw5Kc5hOt7Kts5J7V2/hoVwXRVhNen2bFsnmcPGng3sri7te28o+1B3nnB19jwoi++VF1eXxc+LfPOFzj5I3bTiInJaZPPmeg+WBHOXet3sqhaiffnp3FT87Oi/hpNq/Pz0e7Knih8BD/LirH69fMzE7k4oIczpuZSaLdGtHyDRdOt487XtjI21uPcMLYFNbuq+acGaN4YOlsOR0uhiRpMyW65XT7eOjDPTz80VfYLCbuOG0SF87O4rLH1nKgqpFnr1/AzJykSBezna8qGjj9jx9z6bwcfnnh9D79rANVjZzzl0/ITY1l1U0Libaa+/TzIulwjZOfv76Nd7aVMWFEHL84f1qf1Pr1VlVDE69uPMyqwkPsOFJPlMXEGfkZXFKQzaLxaWE95StalNe7uP7JdWwuruF/z87jupPG8sjHX/Hrt3dw8dxs7r9ohqx7MeRImBJd+tf2Mu55fRvFR51cMCuTn56dx4iEaADK61xc9PBnNDb5ePGmhYxLj4twaVu78alCPtldyYf/vbhfakze217G9U8W8p35o/nVt/s2vEWCx+dnxSf7eODfu/Frze2nTuS7J43DZhnYp2601mwtqWPVukO8uqGEOpeXrCQ7F83J4uK5OYxOjVxNYpPXx6FqB/sqHeyrbMBqNnH5CWMG/DrtzM4j9Vy78kuqG938aekszsjPaJ73x/d28cC/d3PlgjHce36+XIEphhQJU6JDh6od3LN6G//eUc6kkXHce/40FoxrX/uwr7KRix/6jGirmZdvPpGRgaAVaV/ur+aSh9fww9Mmcdtx3O38eN3/zx387cO9/PbiGVxSkNNvn9vXvtxfzc9e2crOsnq+mTeCu8/NH5SnM10eH+9tL2PVumL+s7sCrWHBuBQumZvDWdMziLGFv6moz68pOepkX1Uj+yoa2F/l4KvKRvZXNlJ81IG/zW42PzOBB5bOZsKIgXVw0p2Pd1Vwy9PrsdvMPH71PKZnJ7aar7XmV2/vYPnHX3HDyeP4yVlTJFCJIUPClGjF5fHxyEdf8bcP92AxKX7wzUksW5TbZcPRLcW1LF2+hpyUGJ6/cWHE26Vorbnwb59RWuvkwzsXY7f13yk3r8/PlY9/wfqDR3nl5kVMzUzot8/uC1UNTfz67R2sWldMVpKde87L57SpQ6OR/eEaJy+vL2bVumIOVDmIi7JwzoxRXFKQzZzRycf0Q6+1pqyuiX2VjeyrbGR/VSNfVRj9g1UO3L6W2z7ERVkYmxZLblosY9NiGZsWw9i0OMamxvLF/mr+58VNOD0+7j43n6XzcgZF4Hh67QHuem0bE0fEsWLZPDKT7B0up7Xmrte28dTnB/j+qRP5r9Mm9XNJhegbEqZEsw92lnPP6m0cqHJwzoxR/OxbU8lI7FlN0ye7K7lm5RfMzknmyevmR7TN0BubD3PrMxu4/+IZLIlA7VBFfRPn/OU/RFvNrL71pIiHy+Ph92ueLzzEr9/eQWOTl+tPHsdt35jQJzU3kaa15ot91axaV8ybm0txenyMS4/lkrk5fHtOVqva1qONbr4KBqZAPxieHG5f83I2i4nc1JhAWIprDky5aTGkx0V1GZDK6lz88IVNfLKnkjPyR/Lrb88g+Tgff9TX/H7Nr94u4tH/7OOUyek8eNkc4qK63kb8fs2PXtrMqnXF/PisKdz09fH9VFoh+o6EKUHxUQf3vr6dd7eXMT49lnvPn8aiCR0+j7pLr286zO3PbeC0vJE8dMXciFy10+T18c0/fESszcKbt38tYlcOFe6vZunyz1k8ZQTLr5w7KGoXgrYdruVnr25lw8EaThibwv+7YBoTR4b/SkiPx0NxcTEulyvs7328/FrjdPtwuH00ef0ojGCkMWodQ0/JKcBsUljMCovJFOgb42ZlQimIjo4mOzsbq/XYArXfr3n8k33c/84OUmOj+MOSmZx4HH+Tfcnh9vKD5zby7vYyrlo4hrvOmdrjW1/4/JofPL+R1zcd5ufn5XN1Lx/xJESkyU07h7Emr4/H/rOPv7y/G4XiR2dO4bqTxh5349dzZ2ZS1dDEPa9v52evbuW+C6f1e4h4as0BDlU7efLayN6RvCA3hZ+ence9b2znkY+/GhRH3w1NXv7w7i5WfraP5Bgbf1gykwtnZ/XZ/2FxcTHx8fHk5uYOyLDZ5PFx1OGmzuXFYlJEWUzYLGaiLCaiLCasFhOmLsqttaaqqori4mLGjh17TJ9tMimuP3kcC8encvtzG7j88bXcePJ47jht0oBonF5e5+K6JwrZeriWu8+dyjWLju37mU2KPyyZicvj4+7V24i2mrh03ug+Kq0YDrTWaG0cEPkD/eC42aQierZEwtQQ9p/dFdz92ja+qmzkrGkZ/OycqWR10s7hWCxbNJaKhib++sFe0uOjuKMf20TUONz85f09nDwpfUDc++qaRbmsO3iU+/+5g5nZSQPy9gFg7ITe3FLKL97YTnl9E5fNH83/nDGFxJi+PT3pcrkGbJACiLKayUi0k5HY/bIdUUqRmppKRUXFcZdhWlYib9x2Er94o4iHP9rLp3sqeWDprIheOVtUWsd1K7+kxunh0SsL+OZxtqGzmk08eNlsrn9yHT9+eQvRVjPnz+r9czPFwFTjcLPxUE1zV3LU2Sr4+LXG7zf2R/6QUKS1xqc1fn/nYSk4rTNXLBjN/7sgcldYS5gagkprnfy/N4p4c0spY9NieeLa+Xw9zMHjztMnU1HfxJ//vZv0OBtXLswN6/t35sH391Dn8vCTs6b0y+d1RynFby6awY7SOm57dj1v3v61AXO1Y9D+ykb+77Wt/Gd3JdOyEnjkygJm9eM9wwZqkAqXcHy/GJuFX317Ol+flM6PX97Mt/78CT8/L59LCrL7ff19sLOc257ZQGyUmRduXMi0rONMmgFRFjOPXDGXZX//gjte2ESUxcyZ0zK6f6EY0Dw+PzuP1LPh4FE2HKph48EavqpsBEApmDgijvHpcZjNCpNSmBSYlUIFhk1KYTLRely1LGsyKVTzdKOvlMLcyfypo3q3nfaWhKkhxO318/dPjXsE+fyaH542iRu+Po4oS/irPpVS3HfhdKob3dy1ehupcVGcPX1U2D8n1MEqB0+uOcAlc7PJGzVwrqCLi7Lw8BVzOf+vn3LL0+t59oYFA+KRGi6Pj4c+3MtDH+0lymzi5+flc8WCMXJ36gHszGkZzMpJ4o4XNvI/L23mw13l/OrCGX1egxj01Jr93L16G1MyElixbF6PL07pjt1m5vFl87jy8bXc9ux6ll9VwOLJI8Ly3qLvaa0prXWx4WANGw8dZcPBGraU1DY/uDwtzsasnGQumpvN7JwkpmcnEh89+C7K6Q1pgD5EfLa3krte28ae8gZOmzqSu86Z2i/3CHK6fVz5+Fo2F9ey8tp5nDi+7xrQ3vrMev5VVMaHdy4O204+nF7bWML3n9vIdSeN5f/OmRrRsny8q4K7XtvK/ioH583M5GffarkRa38qKioiLy+v3z831P79+znnnHPYunVrj5ZfuXIlp59+OpmZPX8mZbi/p9+vWf6fr/jdOztJj4/ij5fO6vAecOHi82t++WYRKz7dxzfzRvDA0tnEdnPF3vGodXq47NHP2VPewN+v6dv9hTh+jU1eNhfXsvFQDRsOHmXjoRrK65sA42KNaZkJzMpJZvboJGblJJGdbB/yNdAgDdCHtLI6F798s4jVmw6Tk2JnxbKCfn0Qr91m5rGrC1jyyBpueHIdz92woNenBTqy4eBR3thcyu3fmDAggxTA+bOy2HCwhsc/2cfcMcl9XlPXkSO1Ln7xxnbe3FLKuLRY/nHdCZw0UX6wjsXKlSuZNm3aMYWpcDOZFDd9fTwnjk/l+89t5DuPfs7Np4znB9+cFPZaz8YmL99/biP/KirjmkW5/OxbU/us9jLRbuWp607g0kfW8N0nCnnquvnMHZPSJ58lesbv1+ypaGDjwRo2BGqddpXVN7dPyk2N4cTxqcwencysnCTyRiUMiAskBpohWzPV0ORlS3Etbp8ft9ePJ9APjodO8/j8NPnaTtO4vX6a2iznbvM+wXlKKZJirKTG2kiJtZESG0VKrJWU2ChSY20kx9pC5tl6fdWBx+fnic/286d/7cbt8/O9r4/ne6eMj9jVDKW1Ti7622e4fZqXvreQMamxYXtvrTVLHlnDvspGPvzvxd3e4yaS3F4/ly5fw64j9ay+7STG91MjYq/PzxNrDvCHd3fi9WtuXTyhz07xHouBUjN15plnMnfuXNavX09+fj5PPvkkRUVF3HHHHTQ0NJCWlsbKlSv59NNPWbZsGVlZWdjtdtasWcNvf/tbXn/9dZxOJyeeeCKPPPJIu6PwvvyejU1e7n19O88XHmJmThIPXDqL3LTw/H0dqXVx3RNfUlRaxz3n5XNVP7V9LK9zcenyz6msb+KZ6xe0u5O66DuVDU3NwWnjoRo2H6qlvskLQEK0hZk5ScwenczsnCRm5iSRMkDvfxYJw/I+U1tLajnnL5/0eHmbxUSU2bgU2mY2YbOYsJoVNosZm1lhswSnGfOtgeWD0zSao40eqhqbAn03Rx1ufJ1cfhBjMzcHq+YuxkZKnBG6kmNspMYFQlmMjQS7pXkH/sW+au56bSs7jtSzeHI695yXH9bwcrz2lNdz8cNrSLRbefGmE8P2rLx3th3hxqfW8csLp3H5CWPC8p59qbTWybf+/AmpsTZ+v2Qmfm2EHY9P4/NrPH4/Xp/G5zemeQPjXr/G6/MH+sZyPp/GE5ju8+uW5YPLBl6340gdu8oaWDw5nZ+fNy2iz6ILFRoyfv76NrYfrgvr+0/NTODuc/O7XGb//v2MHTuWTz75hEWLFnHttdeSl5fHK6+8wmuvvUZ6ejrPP/8877zzDitWrOCUU07hd7/7HQUFxj6zurqalBSj9uTKK69kyZIlnHvuuZ1+z77y1pZSfvzSZnx+zc/Pn8ZFc3p3S4tth2u5bmUh9S4PD142h8VT+rcN0+EaJ5c8vIZGt5fnb1jI5Izw3+dsuKt1eNhSUsvmEiM0bSmppaTGCRi3rpiSER84VWfUOo1Li5UHVHdhWJ7mG5sWy7PXL8BmUdjMZqwW1RySbCEhyGYxYTGpPjnf6/dr6l1eqhqbqG50N3dVjW6OhgxXN7rZXdZAdaMbp8fX4XtZTIrkWBsJ0Rb2VjSSlWRn+ZVzOW3qyAFzrnrCiHhWLJvHZY9+zrK/f8FzNyzodSNEj8/Pr9/ewYQRcVw6SJ6DNyrRzp+XzubKFWs578FPw/KelnY3jjQ1T7OaTcRHW3j4ijmckZ8xYLaHgSQnJ4dFixYBcMUVV3DfffexdetWTjvtNAB8Ph+jRnV8WvaDDz7g/vvvx+FwUF1dTX5+frsw1R/Onj6KWTlJ/NfzG7lz1SY+3FnOLy+cflx3339/Rxm3PrOBRLuVVTedGJFHImUm2Xnm+hNY8sgaLn9sLc/fuKDfanKHooYmL1tLatlSXMvmklq2FNewv8rRPH9MagyzRydx9YljmJWTzPSsxH59DNdQN2TDVGyUJeL3/DGZFIkxVhJjrIzr4Z0JnG4f1Q431Q1uo9/YRFWDUctV3eimqsHNuTMzufHk8QPyD2HO6GQeunwu332ykJv+sY4Vy+b16lTTs18cZF9lI49fXdDjOy8PBCdNTOPN275G8VFH+xBkVlhNJswmhdWssJhNrcNSSEgym1Sfhf3+1l0NUl9qu/7i4+PJz89nzZo1Xb7O5XJx8803U1hYSE5ODvfcc09E7+ZuBJAFPPzRXv743i42HKzhj5fOYv7Ynrc7WvnpPu59YztTMxN4/Op5Eb2Vx5jUWJ7+7gIufWQNlz+6llU3LRwwD9f2+Py4PD7ioiwD7u/P6faxvbSWzcUt4WlvRQPBE01ZSXamZyVySUEOM7OTmJaVQFKMnK7rS0M2TA1WdpuZLJs9LDfXjJTFU0Zw/0Uz+OGqTdzxwib+snT2cVUd17k8/Olfu1k4LpVv9PMpiHCYmpkw6B+CPFQcPHiQNWvWsHDhQp555hkWLFjAo48+2jzN4/Gwa9cu8vPziY+Pp76+HqA5OKWlpdHQ0MCLL77IxRdfHMmvgtmkuGXxBBZNSOP7z21g6fI13Lp4ArefOrHLAw6fX/OLN7az8rP9nDZ1JA8snTUgnsM4YUQc//juCSxd/jmXPfY5L9y4kFGJkdn/aa1Zf7CGVzYU88bmUmocHqKtJtLiokiPjyI92I9vP54WF9UnbVabvD52HqlnU7FR27S5uJbd5Q3NTUjS46OYmZ3IuTMymZGdyPTsRNLiwtPEQvRct39JSqkVwDlAudZ6WgfzFfAAcDbgAJZprdeHu6BicLlobjaVDU386u0dpMXauOe8/GM+unv4w71UN7r56dl5A+7IUAwukydP5q9//SvXXnstU6dO5bbbbuOMM87g9ttvp7a2Fq/Xyw9+8APy8/NZtmwZN910U3MD9Ouvv55p06aRkZHBvHnzIv1Vms3KSeLN27/GPau38ef39/CfPZU8cOnsDtvLNTR5uf3ZDby/o5zrvzaWH5+VN6DuN5Y3KoEnr53P5Y+t5fJH1/L8jQvD1uayJ/ZVNvLqhhJe3VjCgSoH0VYTp0/NID8zgapGNxX1TVTUN3GgykHhgaNUN7o7fJ/4aEuXoSstLooR8VGkxNo6DL4en5/dZQ1sKakJhKdadhypw+MzglNyjJUZ2UmcNnUk07MSmZGdNGCvbh5uum2ArpQ6GWgAnuwkTJ0N3IYRpk4AHtBan9DdB8t9poY+rY171zz2yT7++4zJ3LJ4Qo9fe7jGyeLffcjZ00fxx0tn9V0hRZ8aCFfz9YdIf8/XNx3mp69sQWv4xQX5XDg7u3leaa2Ta1cWsqusvvnGrQPVl/uruerxLxidEsNzNywguQ+vJKtudPPG5sO8sqGEDQdrUAoWjU/jgtlZnDkto8urhj0+P9UhIauivomKhqYOxxsCV8qFUgpSY23NNV4psTYOVjvYfriu+UaY8dEWo6YpKynQTxw293MaqHrVAF1r/bFSKreLRc7HCFoa+FwplaSUGqW1Lj2+4oqhQinFT8/Oo7Khid++s5PUWBtL5/fsQae/e3cnGvjh6f333D8hBqtzZ2Yye7TROP2/nt/Ehzsr+MUF0zhY5eC6J76kscnHimXzwv5YqXCbl5vCY1cXcM3KL7lqxRc8ff0JJITxTtouj49/F5XzyoZiPtxZgdevmZIRz0/OmsL5s7J6XMtjNZsYmRDdo/ZmTrePyoYmyrsIXfsqG8lMtHPFgjHMyDZqnMakxMiVdYNIOE6YZwGHQsaLA9MkTAlMJsX9F8+k2uHhp69sITUuitO6eWjq1pJaXtlQwo0njyc7eWA0RhVioMtOjuG5Gxbytw/28Kd/76Zwv3E6KiXWxkvfO2HQ3Hpg0YQ0Hr5iDjc+tY5r/v4lT147v1d3Y/f7NWv3VfPqhhLe2lJKfZOXkQlRXHfSWC6YndXnj6ay28zkpMQMmIb1om/0a+tDpdQNwA0Ao0f3rIZCDH42i4mHLp/DZY+t5dZn1vPUdSd0evWR1pr73ioiyW7l5sXj+7mkQgxuZpPitlMnsmhiGv/1/EYmZ8Sz/Kq5jIgfXO1qvjFlJA8snc2tz6zn+icLWbFs3jE37t5dVs/LG0p4bUMJh2tdxNrMnDltFN+ek8WCcakDqs2YGPzCEaZKgNAbAGUHprWjtV4OLAejzVQYPlsMErFRFv6+bB4XP/wZ333iS164aSFTMtofEX64s4LP9lZxz7lTw1q9L8RwMmd0Mu//8BRMqv1tIQaLs6eP4vdLZnLHC5u46R/reOTKud3eZqW8zsXqTUY7qG2H6zCbFCdPTONHZ03h9KkZA/J2MmJoCEeYWg3cqpR6DqMBeq20lxIdSYm18eS187nooc+4esUXvPS9E1udxvP6/Nz3VhG5qTFcNgjudC7EQDYUal4unJ2N0+3np69s4fvPbuTBy2a3uwquscnLu9uP8PL6Ej7dU4lfw4zsRO4+dyrnzMjs16sCxfDVk1sjPAucAqQppYqBuwErgNb6YeAtjCv59mDcGuGaviqsGPyyk2N44tr5XPLwGq5a8QUv3nRi87OfVq0rZnd5Aw9fMUcepCmEAOCyE0bj8vi4943t/HDVJv6wZBZaaz7dW8WrG0p4Z9sRHG4fWUl2bj5lAhfMzmLCCLmTuuhfPbma7zvdzNfALWErkRjypmQYd16+4vG1XLPyS569/gS0ht+/u4uCMcmckZ8R6SIKIQaQa08ai9Pj47fv7KS01sW+ykYq6ptIiLZw/qwsLpydRcGYZLn6TURM5G9/K4al+WNT+Mt3ZvO9f6zje/9Yz7SsBCobmlh+1dxB28ZDCNF3blk8gSavn4c/3Mspk9O5cHYWi6eM6JO7jgtxrCRMiYg5Iz+DX144nZ+8vIWPdlXwremjmDM6OdLFEkPM/v37OfPMM1mwYAGfffYZ8+bN45prruHuu++mvLycp59+mrfeeou4uDjuvPNOAKZNm8Ybb7xBbm5uZAsvWrnjtEl8/9SJQ6I9mBhaJEyJiPrO/NFUN7r5+6f7+J8zJ0e6OKIvvf1jOLIlvO+ZMR3O+nW3i+3Zs4dVq1axYsUK5s2bxzPPPMMnn3zC6tWrue+++5g1a1Z4yyX6jAQpMRBJK18RcbcsnsAXP/0mY1JjI10UMUSNHTuW6dOnYzKZyM/P59RTT0UpxfTp09m/f3+kiyeEGOSkZkoMCNJwdBjoQQ1SX4mKark83mQyNY+bTCa8Xi8WiwW/39+8jMvl6vcyCiEGL6mZEkIMe7m5uaxfvx6A9evXs2/fvgiXSAgxmEiYEkIMexdddBHV1dXk5+fz4IMPMmmSPGBbCNFzyrhNVP8rKCjQhYWFEflsIUT/KCoqIi8vL9LF6HPD5XsKMZwppdZprQs6mic1U0IIIYQQvSBhSgghhBCiFyRMCSGEEEL0goQpIYQQQohekDAlhBBCCNELEqaEEEIIIXpBwpQQQgghRC9ImBJCCCGE6AUJU0KIIW3//v3k5eVx/fXXk5+fz+mnn47T6eSUU04heOPgyspKcnNzI1tQIcSgJQ86FkL0i9988Rt2VO8I63tOSZnCj+b/qNvldu/ezbPPPsujjz7KkiVLeOmll8JaDiHE8CY1U0KIIW/s2LHMmjULgLlz57J///6IlkcIMbRIzZQQol/0pAapr0RFRTUPm81mnE4nFosFv98PgMvlilTRhBBDgNRMCSGGpdzcXNatWwfAiy++GOHSCCEGMwlTQohh6c477+Shhx5i9uzZVFZWRro4QohBTGmtu19IqTOBBwAz8JjW+tdt5i8DfguUBCY9qLV+rKv3LCgo0MEraYQQQ1NRURF5eXmRLkafGy7fU4jhTCm1Tmtd0NG8bttMKaXMwF+B04Bi4Eul1Gqt9fY2iz6vtb6116UVQgghhBhEenKabz6wR2v9ldbaDTwHnN+3xRJCCCGEGBx6EqaygEMh48WBaW1dpJTarJR6USmVE5bSCSEGvZ40JRjMhvr3E0J0L1wN0F8HcrXWM4D3gCc6WkgpdYNSqlApVVhRURGmjxZCDFTR0dFUVVUN2cChtaaqqoro6OhIF0UIEUE9uc9UCRBa05RNS0NzALTWVSGjjwH3d/RGWuvlwHIwGqAfU0mFEINOdnY2xcXFDOWDp+joaLKzsyNdDCFEBPUkTH0JTFRKjcUIUUuBy0IXUEqN0lqXBkbPA4rCWkohxKBktVoZO3ZspIshhBB9qtswpbX2KqVuBd7BuDXCCq31NqXUvUCh1no1cLtS6jzAC1QDy/qwzEIIIYQQA0aP7jPVF+Q+U0IIIYQYLLq6z5TcAV0IIYQQohckTAkhhBBC9IKEKSGEEEKIXpAwJYQQQgjRCxKmhBBCCCF6QcKUEEIIIUQvSJgSQgghhOgFCVNCCCGEEL0gYUoIIYQQohd68my+wclVBxv+AWmTIG0iJOaAaZBnR78f6krAGgOxqZEujRBCCCEYymGqYie885OWcYsd0iYEwlVIlzoerPbIlbMtrcFRBVV72nRfQfVe8LqM5VLGQ84JMPoEo582efCHRSGEGC58Hqg5CEf3wdEDoP1giQJLtNE3RwXGo0LGozuYFgVKRfrbDHtD+9l8jZVQuSvQ7TYCVuUuYwMm+L0VJI02glX6ZKMWKxi0YtP6rmxNDUY4qtoDVXtbBydXbctyJisk50LqBCP4pY435h/6Ag6tNYIXQFQi5MwzglXOCZA1F6Li+q78veX3wdH9ULEDyougvhRi0iB+JMRltPTjRoDZGunSCiHCzecBjxO8TeB1gsdl9L1NgemujudrbewX4jMgfpTR2ZMHZqBoqofqfUZgatuvLTYCVDiYbUbQCvY7ClzNnR2iEyA6MaRL6mA8Qfa9bXT1bL6hHaY643EaAaZypxGymgPXHuOPNcie0nKaMBiw0idB0hgwmbv/HK/bCAytapgCwanhSOtlE3MCYWlCSDceEkeDuZMKRK2h+is4+LkRrA59ARVFxjxlgpHTWsLV6BOMz+jvHY7fD7UHjcBUXtQSnip3tdSygfHH66qlJeQGKYhJNXaccSND+qPaBy9rdD9+sTDTGnxucDeCxwFuB3gaA31Hx9NCl/U6jR2z9hvvBSHDugfDIctrf2BeJ8MAqMC21ElfmTqYR8i4qWevTxoN2QWQPc8YHog/mEOR1kaIcTeCu97oNzW0GQ50TQ0t22Kr4BPogiGobVjSvvCV12wL7AsyQkJWcD8RGI8bGf7QpTU0lHUemIIHu0H2FEgZC8ljW/rJuUZnthnry+duWXded5tpIfPaTWsCX1Pg/6CpzTKh051GMxhXbff/B7a4NiGrqwAW0tmTICqhZ7+Tg4iEqZ7y+6H2UJuAFegaK1qWM9uMsNMcsiZDTLLxBxRay1RzoPWRR0xa6xqmYGhKGRe+U43Oo1C8LhCuPjeGPY3GvPhRkDO/JWBlzACLLTyfq7Wx7sp3GIEu2K/YaexkgxKyIH0KjMgL6U+GqHjjSLWxAuqPGDuo+iOB4SNQXxbSL+t4JxCd2DpcNffbBLHgZ/k9gb7X2Ok0D7edFzLe6Txvx8t4XZ0HpOB0j9MYPtYfF0u00X7OGmNsPyYz7cMMrYNNj4Y7Cjxtwk/boNVZv8t5Xbw+OL/6q5YDnLiRRqgKhqvM2WCLPbZ1NpD4/cb/ufYbNbU6MO73Bb6/r4Pp/k6W72S6t6nrAORuMGpPmodDlvF7e/Y9TBbj/8EaE9gm7YHakWjjAMdiD/SjO5nfybTOXgfQWB7YP5Qa+4T60pDxwD4jtIY/yBwVErbahq+Q0BWd2BK6vG5j39ZhYNrf+gBcmSAhG1JyWwemlEBoik7sxQYTZlob/8+u2o47Z03IeE3Hy7Q7+G3DEo2xEwpoFWR7M53202dfAWf8suvy9JKEqS74tR+n14nD48DhddDoaWwedngdxrDHgcNZhaP+MI31pTiclThdtTjc9Th8LhqVQivI8njJ8StyolLIic8mJ3kSmSNmYE2fDKnjjKOi/ubzQvm2ltOCh9YGTnNibOiZc0IC1vzuT21qDXWHWwem8kBocje0LBeXASOmQHpeSz99snHE0lt+v3HE13Cki8AV6Puaev95x8tkMX4QbIHAE/zBsdpbhm0xYI1tv4wt1liu1bQ2yw6xo74O+TxQtg2Kv4TiQqNfvdeYp8wwMj8QsAJd6viBUXsVrLEIbWYQ7OpKwlsrc6zMUUYTAFss2OJDhuOMA43m4Tij3zzcdvnAsNk2MNZ5W25H631EaNiqL205YGuqa/9ai90IV9rX/nScxd4SjtoGpsSc8B2gDnR+vxHUuwpioQfSocGrs9zRanpny3cyffQCmPbtY/8ex2BYhqnShlJWbltphKNAMHJ6nC3jgcDkDD2q6IbVZCXGGkOsJZYYawwxlhjslmhi/H60z02xp47ixiO4fC2nr0zKxKjYUWTHZ5MTn9Oui7VG4Mi6rrTltOChtVC6yahJAaNh++gFRrDKnGOEloodUL49EJ52QlPIEV9segc1TVMgJqX/v1dbWhtHVM2BK7Dz9DiMtmjmQGeyGqdSg9NMlsA8W+fz2r2+g3kD8QdmKGisgpJAsCr+0qh9ddcb8+zJIeGqwGg72Je1AV63UXsWDE1Ve1rCU+iPtDW25QKYxGxjGzGZjUCoVMiwqQfTTcbFJt1M9ypQFjvmqAQj/ATDkbSDaa2pIaQmvE0Nl1LtA1PcSPnbHqaGZZjac3QPV//zamKtscRYYprDT4w1ZDgwHmuNxW6xN0/v8DWWGKw92AlpralwVnCo/lCrrri+mOL6Yo42HW21fEp0SqdBKzU6FdUff7QeJxzeaJwWbNuwPcie0jowjcgzapsG4C0atNZGbWOgpjG0trHR04hP+4gyRxFljsJmtrX0TR1MM0dhHgA1QFprvH4vHr8Ht8+Nx+9p7tw+N37tZ0TMCJKikvpnmxlI/D4j5DeHq0LjAAANKGObDZ4azJ5n1JAe6/+pozqkhimkf3R/61qm+MzWbSyDwwmZffoD7PK62H10N0XVRRRVF7Gjage7ju7Cj5/suGxGJ4xmdPzoln78aEbFjcJiGlwXdDs8Durd9SRHJ2MzD5MaIDFgDMswNVDVu+vbhazg8JHGI+iQKky7xW4ErbjWISsxOhGFwqRMqMD54uCwSZkCTWVa5iulmvsmTM0/tp3OB1TNAdSRLaiYNOPHJyat0x+Dnvx4K7pexqd9LadXPY7mGsRgGHJ6nc3DobWLnQ3r7s7lHwOLsrQLWMF+R9NsZhs2kzHsx4/b5+44CPk8uP3u5uHQ6R5/YF7I9J6INkeTEZtBRmwGo2JHMSp2VKvxjNgMooPtToYyVy2UrGs5NVj8Jdp5lCalcEbF4xo1A2dGPs70ybhSxuG2RZMenUKm10dMzaE2p+d2tT64aNdmMhCaUicYp8n6WL27nh3VO9hRvYOiKiM87avdhy8Q6uJt8UxNmcqUlCmYTWYO1R/iQN0BDtUfalUTb1EWsuKzyInPaR20EkaTGZeJ1dR/NVhaa442HaXcUU5ZYxlljjJjONAPTq/31De/Jt4WT5o9zeii00i1p7aM21vGk6OSI3JA5PP7qGmqocpVRbWrmipn637o8NGmo9gtdtLt6aTHpLfrj7CPID0mnVR7ar/+v3TF6/dy1HWUSmdlc1flqqLKWdU87vA6sJls2Mw2rGYrUabA/tFsw2qyNu8vrSZrq/1nq3GT8dp244F9bHCe3WzvUYVHb0iYGiTcPjclDSXtarMO1R+iuKGYpki2/xkg7BZ7c81hrDW2VW1isIYxONx2PLicSZlw+924fW6afE00+Zqah491Wqfz/EbfpEzGjiGwA7CarC39QBfcebSd13a6zWRr//rAdI2mwlFBaWMpRxqPcKTxCKWNpVQ4K9qtw+So5FbhalTsKDLiAv2YDNLsaRH58fH4PDi8DlxeF06vE5fP6Du9zpZpXleX04Ovb7+8A2fo1aNdSPH5yPJ4ycRMpi2J7LhMMpPGkzViBpmZ84hKndBv7dWqnFVGaKouag5Oh+oPNc9Pt6czJWUKeal55KXkkZeaR2ZsZocHOFprKp2VHKw/yMG6g839YNhyeFvat5iVmVGxo9rVaOUk5JAdl31MtUIen4cKZwXljnKOOI5Q3ljeKigF+20PGEzKRJo9jZExIxkRM6K5H2+Lp6appuUHvM2Pd1smZSI5KrldyEqNbh++EmwJXR4curyu1sEoJCi1DU01TTX4O7j1gUVZSIlOIdWe2txPikrC4XVQ4aigwllBhaOCKldVh69PiU4h3Z5OWkwaI+wjSLOnMSJmREvwihlx3KFLa01tU62xPl0t6zZ0HQenH3Ud7fCgNc4a17w+YywxzQeFbp+7eb/b0XBH3/VYLJm0hP9b+H+9eo/uSJgaAvzaT4XDOH1Y765HB//plr4fP2hj2bbz/dpP8P86OD+48QZfG/peofN7qyfbmFKqXRAKPQ1rtE+zD4hTboOJx+ehzFHWHLI66jcGr/YMsCgLI2NHdlrDlWBLaBVkgsPBNogdhaHQLjTohE7z6h5ePRZgUiaizdHYLXaiLUY/ONzddLvFTjQKe+1hoqv3Yak5SHlUDIejYyhWcNjvoMRRzuHGw3jbXNWWbk8nMy6TrLis5i4zLpPsuGwyYjOO6+hYa82RxiNsr97eqsap3FHevExWXBZTU40ap2BwSrOH5154WmuqXFUcqj/UKmgF+w2elotLgu1AQ2u0suKyaPA0NNcgBUNSmaOMald1u8+LNkczMtYIR6FBaWTMyObhVHvqMZ+GdHgc7WpH2oauKpfR76i212qytoSr6FRMytQqNLX9WwmKtcYawSg6tVVQahuaUqO7D2xBPr+Palc15c5yKh2VnfYrXZWdhq40e1pLDVcgbCXaEpvDaNt1VeWqare9A9hMtnZhNHQ8GExT7anYLcd3ZbrX7+00aDUPh463mT4xeSILMxce12f3lIQpIUSn6t31LeGqoZQjDiNklTYY08ocZc2nkI6FRVlaBZmOhkO7rpbpKBxZTdY+bx/m137KHeUcbjhMSUNJcxccP9J4pNW6MSkTI2JGkBmbSXZ8Nplxma2GR8aMxKRMHKg70Co0FVUXURu4sMOkTIxNGEteah5TUqYwNXUqk1Mmk2BL6NPv2hmtNTVNNc2nCpuDVt1BDtQfoN5d32r5pKik9gEptnXtUk8DRV/RWlPnrmsfugK1LlXOKiqcFfi1vyUMRae2Gg4GpOTo5OMOEOHg8/uaT5FWOispd5S3quEKrelqu62GhqCOauqCw3HWuOHXFrMDvQ5TSqkzgQcAM/CY1vrXbeZHAU8Cc4Eq4FKt9f6u3lPClBCDg8/vo9JZ2Ry46j31HYah5tBjbQk7Q53X76XcUd4uaBXXF3O48TBljWWtToWYlRmrydp8xa/VZGVi8kSjpikljympU5iUPCmiP87HqraplpKGEuKt8aTHpA+PNnmDUDB01TbVkhSVRFJUktT0H6NehSmllBnYBZwGFANfAt/RWm8PWeZmYIbW+ial1FLgQq31pV29r4QpIcRQ5/F5ONJ4hOKG4ubaLKfXyaTkSUxNncq4xHF93mhWCBEeXYWpnpyQng/s0Vp/FXiz54Dzge0hy5wP3BMYfhF4UCmldKTOIQohxABgNVvJScghJyEn0kURQvQhUw+WyQIOhYwXB6Z1uIzW2gvUAgPvBkRCCCGEEGHWkzAVNkqpG5RShUqpwoqK9pdsCyGEEEIMNj0JUyVAaB11dmBah8sopSxAIkZD9Fa01su11gVa64L09PTjK7EQQgghxADSkzZTXwITlVJjMULTUuCyNsusBq4G1gAXA+93115q3bp1lUqpA8de5GOWBlT2w+cMdLIeWsi6aCHrooWsC4OshxayLlrIuoAxnc3oNkxprb1KqVuBdzBujbBCa71NKXUvUKi1Xg08DjyllNoDVGMEru7et1+qppRShZ21vh9OZD20kHXRQtZFC1kXBlkPLWRdtJB10bUe3V5Wa/0W8FabaXeFDLuAS8JbNCGEEEKIga9fG6ALIYQQQgw1wyFMLY90AQYIWQ8tZF20kHXRQtaFQdZDC1kXLWRddCFiz+YTQgghhBgKhkPNlBBCCCFEnxkSYUopdaZSaqdSao9S6scdzI9SSj0fmL9WKZUbgWL2OaVUjlLqA6XUdqXUNqXU9ztY5hSlVK1SamOgu6uj9xoKlFL7lVJbAt+z3YMgleHPge1is1JqTiTK2deUUpND/r83KqXqlFI/aLPMkN0ulFIrlFLlSqmtIdNSlFLvKaV2B/rJnbz26sAyu5VSV/dfqcOvk/XwW6XUjsD2/4pSKqmT13b5tzTYdLIu7lFKlYT8DZzdyWu7/L0ZbDpZF8+HrIf9SqmNnbx2SG0XvaK1HtQdxu0a9gLjABuwCZjaZpmbgYcDw0uB5yNd7j5aF6OAOYHheIwHVLddF6cAb0S6rP20PvYDaV3MPxt4G1DAAmBtpMvcD+vEDBwBxgyX7QI4GZgDbA2Zdj/w48Dwj4HfdPC6FOCrQD85MJwc6e8T5vVwOmAJDP+mo/UQmNfl39Jg6zpZF/cAd3bzum5/bwZb19G6aDP/98Bdw2G76E03FGqmmh/ErLV2A8EHMYc6H3giMPwicKpSSvVjGfuF1rpUa70+MFwPFNH+OYqixfnAk9rwOZCklBoV6UL1sVOBvVrr/rhh7oCgtf4Y4/53oUL3CU8AF3Tw0jOA97TW1Vrro8B7wJl9Vc6+1tF60Fq/q43nqQJ8jvGEiyGvk22iJ3ryezOodLUuAr+TS4Bn+7VQg9BQCFPyIOYOBE5lzgbWdjB7oVJqk1LqbaVUfv+WrF9p4F2l1Dql1A0dzO/JtjPULKXzHeNw2S4ARmqtSwPDR4CRHSwz3LaPazFqajvS3d/SUHFr4JTnik5O/Q63beJrQJnWencn84fLdtGtoRCmRBtKqTjgJeAHWuu6NrPXY5zimQn8BXi1n4vXn07SWs8BzgJuUUqdHOkCRZJSygacB6zqYPZw2i5a0cb5imF9WbNS6n8BL/B0J4sMh7+lh4DxwCygFOP01nD3HbqulRoO20WPDIUwFbYHMQ8FSikrRpB6Wmv9ctv5Wus6rXVDYPgtwKqUSuvnYvYLrXVJoF8OvIJRRR+qJ9vOUHIWsF5rXdZ2xnDaLgLKgqd0A/3yDpYZFtuHUmoZcA5weSBYttODv6VBT2tdprX2aa39wKN0/B2HxTYBzb+V3wae72yZ4bBd9NRQCFPND2IOHHkvxXjwcqjgg5ihhw9iHowC57cfB4q01n/oZJmMYHsxpdR8jG1gyAVLpVSsUio+OIzR0HZrm8VWA1cFrupbANSGnPoZijo9yhwu20WI0H3C1cBrHSzzDnC6Uio5cMrn9MC0IUMpdSbwP8B5WmtHJ8v05G9p0GvTXvJCOv6OPfm9GSq+CezQWhd3NHO4bBc9FukW8OHoMK7K2oVxlcX/Bqbdi7GDAIjGOLWxB/gCGBfpMvfRejgJ43TFZmBjoDsbuAm4KbDMrcA2jKtQPgdOjHS5+2hdjAt8x02B7xvcLkLXhQL+GthutgAFkS53H66PWIxwlBgybVhsFxgBshTwYLRxuQ6jzeS/gd3Av4CUwLIFwGMhr702sN/YA1wT6e/SB+thD0YboOD+InjVcybwVmC4w7+lwdx1si6eCuwHNmMEpFFt10VgvN3vzWDuOloXgekrg/uHkGWH9HbRm07ugC6EEEII0QtD4TSfEEIIIUTESJgSQgghhOgFCVNCCCGEEL0gYUoIIYQQohckTAkhhBBC9IKEKSGEEEKIXpAwJYQQQgjRCxKmhBBCCCF64f8DDwiu3bazulMAAAAASUVORK5CYII=\n", "text/plain": ["<Figure size 720x144 with 1 Axes>"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["df.plot(figsize=(10,2));"]}, {"cell_type": "markdown", "metadata": {}, "source": ["On pourrait imaginer que les param\u00e8tres cherch\u00e9s soient \u00e9gaux \u00e0 la moyenne des valeurs obtenues pour chaque jour. Encore faudrait-il le prouver ou trouver un contre exemple."]}, {"cell_type": "markdown", "metadata": {}, "source": ["#### Seconde id\u00e9e : probl\u00e8mes d'optimisation\n", "\n", "On cherche \u00e0 transformer le probl\u00e8me d'estimation en un probl\u00e8me d'optimisation de telle sorte que les param\u00e8tres en sont la solution. L'id\u00e9e de simuler en fonction d'un jeu de param\u00e8tres puis de calculer une distance entre cette simulation et les donn\u00e9es observ\u00e9es.\n", "\n", "Le probl\u00e8me de cette solution est qu'il faut essayer plein de valeur, au hasard pour faire simple. Donc on se sert de la fonction pr\u00e9c\u00e9dente comme point de d\u00e9part et on tire des valeurs tout autour."]}, {"cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [{"data": {"text/plain": ["0.09110736444463925"]}, "execution_count": 23, "metadata": {}, "output_type": "execute_result"}], "source": ["def distance_sim(obs, sim):\n", "    return ((obs - sim) ** 2).sum() / obs.size\n", "\n", "distance_sim(sim, sim_bruit)"]}, {"cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [{"data": {"text/plain": ["{'b0': 1.1949449618580996,\n", " 'mu0': 0.18786221862533412,\n", " 'nu0': 0.10660548950608435,\n", " 'beta': array([1.15830814]),\n", " 'mu': array([0.17720381]),\n", " 'nu': array([0.10532943]),\n", " 'distance': 0.021211456859144485,\n", " 'distance0': 0.043139240901541635}"]}, "execution_count": 24, "metadata": {}, "output_type": "execute_result"}], "source": ["def optimisation(obs):\n", "    est = estimation_coefficient(obs)\n", "    b0 = est['beta'].mean()\n", "    m0 = est['mu'].mean()\n", "    n0 = est['nu'].mean()\n", "    n_iter = obs.shape[0] - 1\n", "    S, I, R, D = obs[0]\n", "    beta, mu, nu = b0, m0, n0\n", "    distance = distance_sim(obs, simulation(n_iter, S, I, R, D, b0, m0, n0))\n", "    distance0 = distance\n", "    for i in range(0, 100):\n", "        b = numpy.random.randn(1) * b0 / 5 + b0\n", "        m = numpy.random.randn(1) * m0 / 5 + m0\n", "        n = numpy.random.randn(1) * n0 / 5 + n0      \n", "        sim = simulation(n_iter, S, I, R, D, b, m, n)\n", "        d = distance_sim(obs, sim)\n", "        if distance is None or d < distance:\n", "            distance = d\n", "            beta, mu, nu = b, m, n\n", "    return dict(b0=b0, mu0=m0, nu0=n0, beta=beta, mu=mu, nu=nu,\n", "                distance=distance, distance0=distance0)\n", "\n", "optimisation(sim_bruit)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["`distance0` la distance entre la simulation obtenue avec les param\u00e8tes de la premi\u00e8re m\u00e9thode, `distance` celle de la seconde et elle est plus petite. Cela r\u00e9pond \u00e0 la question qu'on se posait ci-dessus."]}, {"cell_type": "code", "execution_count": 24, "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.9.5"}}, "nbformat": 4, "nbformat_minor": 2}