{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this document we reperform some of the analysis provided in \n",
    "*Risk Analysis of the Space Shuttle: Pre-Challenger Prediction of Failure* by *Siddhartha R. Dalal, Edward B. Fowlkes, Bruce Hoadley* published in *Journal of the American Statistical Association*, Vol. 84, No. 408 (Dec., 1989), pp. 945-957 and available at http://www.jstor.org/stable/2290069. \n",
    "\n",
    "On the fourth page of this article, they indicate that the maximum likelihood estimates of the logistic regression using only temperature are: $\\hat{\\alpha}=5.085$ and $\\hat{\\beta}=-0.1156$ and their asymptotic standard errors are $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$. The Goodness of fit indicated for this model was $G^2=18.086$ with 21 degrees of freedom. Our goal is to reproduce the computation behind these values and the Figure 4 of this article, possibly in a nicer looking way."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Technical information on the computer on which the analysis is run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will be using the python3 language using the pandas, statsmodels, numpy, matplotlib and seaborn libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.10.4 (main, Apr  2 2022, 09:04:19) [GCC 11.2.0]\n",
      "uname_result(system='Linux', node='dev004', release='5.15.0-33-generic', version='#34-Ubuntu SMP Wed May 18 13:34:26 UTC 2022', machine='x86_64')\n",
      "IPython 8.14.0\n",
      "IPython.core.release 8.14.0\n",
      "PIL 10.3.0\n",
      "PIL.Image 10.3.0\n",
      "PIL._deprecate 10.3.0\n",
      "PIL._version 10.3.0\n",
      "_csv 1.0\n",
      "_ctypes 1.1.0\n",
      "_curses b'2.2'\n",
      "decimal 1.70\n",
      "_pydev_bundle.fsnotify 0.1.5\n",
      "_pydevd_frame_eval.vendored.bytecode 0.13.0.dev\n",
      "argparse 1.1\n",
      "backcall 0.2.0\n",
      "cffi 1.15.1\n",
      "colorama 0.4.4\n",
      "comm 0.1.4\n",
      "csv 1.0\n",
      "ctypes 1.1.0\n",
      "cycler 0.10.0\n",
      "dateutil 2.9.0\n",
      "dateutil._version 2.9.0\n",
      "debugpy 1.6.7.post1\n",
      "debugpy.public_api 1.6.7.post1\n",
      "decimal 1.70\n",
      "decorator 5.1.1\n",
      "defusedxml 0.7.1\n",
      "distutils 3.10.4\n",
      "django 3.2.25\n",
      "entrypoints 0.4\n",
      "executing 1.2.0\n",
      "executing.version 1.2.0\n",
      "http.server 0.6\n",
      "ipykernel 6.25.1\n",
      "ipykernel._version 6.25.1\n",
      "ipywidgets 8.1.0\n",
      "ipywidgets._version 8.1.0\n",
      "jedi 0.19.0\n",
      "joblib 1.2.0\n",
      "joblib.externals.cloudpickle 2.2.0\n",
      "joblib.externals.loky 3.3.0\n",
      "json 2.0.9\n",
      "jupyter_client 7.4.9\n",
      "jupyter_client._version 7.4.9\n",
      "jupyter_core 5.3.1\n",
      "jupyter_core.version 5.3.1\n",
      "kiwisolver 1.4.5\n",
      "kiwisolver._cext 1.4.5\n",
      "logging 0.5.1.2\n",
      "matplotlib 3.7.2\n",
      "matplotlib._version 3.7.2\n",
      "numpy 1.22.4\n",
      "numpy.core 1.22.4\n",
      "numpy.core._multiarray_umath 3.1\n",
      "numpy.lib 1.22.4\n",
      "numpy.linalg._umath_linalg 0.1.5\n",
      "numpy.version 1.22.4\n",
      "packaging 23.2\n",
      "pandas 1.4.2\n",
      "parso 0.8.3\n",
      "patsy 0.5.3\n",
      "patsy.version 0.5.3\n",
      "pexpect 4.8.0\n",
      "pickleshare 0.7.5\n",
      "pkg_resources._vendor.more_itertools 9.1.0\n",
      "pkg_resources._vendor.packaging 23.1\n",
      "pkg_resources._vendor.platformdirs 2.6.2\n",
      "pkg_resources._vendor.platformdirs.version 2.6.2\n",
      "pkg_resources._vendor.more_itertools 9.1.0\n",
      "pkg_resources._vendor.packaging 23.1\n",
      "pkg_resources._vendor.platformdirs 2.6.2\n",
      "platform 1.0.8\n",
      "platformdirs 3.10.0\n",
      "platformdirs.version 3.10.0\n",
      "prompt_toolkit 3.0.39\n",
      "psutil 5.9.5\n",
      "ptyprocess 0.7.0\n",
      "pure_eval 0.2.2\n",
      "pure_eval.version 0.2.2\n",
      "pyarrow 9.0.0\n",
      "pyarrow._generated_version 9.0.0\n",
      "pydevd 2.9.5\n",
      "pygments 2.16.1\n",
      "pyparsing 2.4.7\n",
      "pytz 2023.3.post1\n",
      "re 2.2.1\n",
      "scipy 1.7.3\n",
      "scipy._lib._uarray 0.5.1+49.g4c3f1d7.scipy\n",
      "scipy._lib.decorator 4.0.5\n",
      "scipy.integrate._dop 1.21.4\n",
      "scipy.integrate._ode $Id$\n",
      "scipy.integrate._odepack  1.9 \n",
      "scipy.integrate._quadpack  1.13 \n",
      "scipy.integrate.lsoda 1.21.4\n",
      "scipy.integrate.vode 1.21.4\n",
      "scipy.interpolate._fitpack  1.7 \n",
      "scipy.interpolate.dfitpack 1.21.4\n",
      "scipy.linalg._fblas 1.21.4\n",
      "scipy.linalg._flapack 1.21.4\n",
      "scipy.linalg._flinalg 1.21.4\n",
      "scipy.linalg._interpolative 1.21.4\n",
      "scipy.ndimage 2.0\n",
      "scipy.optimize.__nnls 1.21.4\n",
      "scipy.optimize._cobyla 1.21.4\n",
      "scipy.optimize._lbfgsb 1.21.4\n",
      "scipy.optimize._minpack  1.10 \n",
      "scipy.optimize._slsqp 1.21.4\n",
      "scipy.optimize.minpack2 1.21.4\n",
      "scipy.signal.spline 0.2\n",
      "scipy.sparse.linalg.eigen.arpack._arpack 1.21.4\n",
      "scipy.sparse.linalg.isolve._iterative 1.21.4\n",
      "scipy.special.specfun 1.21.4\n",
      "scipy.stats.mvn 1.21.4\n",
      "scipy.stats.statlib 1.21.4\n",
      "seaborn 0.13.2\n",
      "seaborn.external.appdirs 1.4.4\n",
      "seaborn.external.husl 2.1.0\n",
      "setuptools 69.1.1\n",
      "distutils 3.10.4\n",
      "setuptools._vendor.more_itertools 8.8.0\n",
      "setuptools._vendor.ordered_set 3.1\n",
      "setuptools._vendor.packaging 23.1\n",
      "setuptools._vendor.more_itertools 8.8.0\n",
      "setuptools._vendor.ordered_set 3.1\n",
      "setuptools._vendor.packaging 23.1\n",
      "setuptools.version 69.1.1\n",
      "six 1.16.0\n",
      "socketserver 0.4\n",
      "stack_data 0.6.2\n",
      "stack_data.version 0.6.2\n",
      "statsmodels 0.14.0\n",
      "statsmodels.__init__ 0.14.0\n",
      "statsmodels._version 0.14.0\n",
      "statsmodels.api 0.14.0\n",
      "statsmodels.tools.web 0.14.0\n",
      "traitlets 5.9.0\n",
      "traitlets._version 5.9.0\n",
      "urllib.request 3.10\n",
      "wcwidth 0.2.6\n",
      "xmlrpc.client 3.10\n",
      "zlib 1.0\n",
      "zmq 24.0.1\n",
      "zmq.sugar 24.0.1\n",
      "zmq.sugar.version 24.0.1\n"
     ]
    }
   ],
   "source": [
    "def print_imported_modules():\n",
    "    import sys\n",
    "    for name, val in sorted(sys.modules.items()):\n",
    "        if(hasattr(val, '__version__')): \n",
    "            print(val.__name__, val.__version__)\n",
    "#        else:\n",
    "#            print(val.__name__, \"(unknown version)\")\n",
    "def print_sys_info():\n",
    "    import sys\n",
    "    import platform\n",
    "    print(sys.version)\n",
    "    print(platform.uname())\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import seaborn as sns\n",
    "\n",
    "print_sys_info()\n",
    "print_imported_modules()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading and inspecting data\n",
    "Let's start by reading data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>Date</th>\n",
       "      <th>Count</th>\n",
       "      <th>Temperature</th>\n",
       "      <th>Pressure</th>\n",
       "      <th>Malfunction</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>4/12/81</td>\n",
       "      <td>6</td>\n",
       "      <td>66</td>\n",
       "      <td>50</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>11/12/81</td>\n",
       "      <td>6</td>\n",
       "      <td>70</td>\n",
       "      <td>50</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3/22/82</td>\n",
       "      <td>6</td>\n",
       "      <td>69</td>\n",
       "      <td>50</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>11/11/82</td>\n",
       "      <td>6</td>\n",
       "      <td>68</td>\n",
       "      <td>50</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>4/04/83</td>\n",
       "      <td>6</td>\n",
       "      <td>67</td>\n",
       "      <td>50</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>6/18/82</td>\n",
       "      <td>6</td>\n",
       "      <td>72</td>\n",
       "      <td>50</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>8/30/83</td>\n",
       "      <td>6</td>\n",
       "      <td>73</td>\n",
       "      <td>100</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>11/28/83</td>\n",
       "      <td>6</td>\n",
       "      <td>70</td>\n",
       "      <td>100</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>2/03/84</td>\n",
       "      <td>6</td>\n",
       "      <td>57</td>\n",
       "      <td>200</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>4/06/84</td>\n",
       "      <td>6</td>\n",
       "      <td>63</td>\n",
       "      <td>200</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>8/30/84</td>\n",
       "      <td>6</td>\n",
       "      <td>70</td>\n",
       "      <td>200</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>10/05/84</td>\n",
       "      <td>6</td>\n",
       "      <td>78</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>11/08/84</td>\n",
       "      <td>6</td>\n",
       "      <td>67</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1/24/85</td>\n",
       "      <td>6</td>\n",
       "      <td>53</td>\n",
       "      <td>200</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>4/12/85</td>\n",
       "      <td>6</td>\n",
       "      <td>67</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>4/29/85</td>\n",
       "      <td>6</td>\n",
       "      <td>75</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>6/17/85</td>\n",
       "      <td>6</td>\n",
       "      <td>70</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>7/2903/85</td>\n",
       "      <td>6</td>\n",
       "      <td>81</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>8/27/85</td>\n",
       "      <td>6</td>\n",
       "      <td>76</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>10/03/85</td>\n",
       "      <td>6</td>\n",
       "      <td>79</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>10/30/85</td>\n",
       "      <td>6</td>\n",
       "      <td>75</td>\n",
       "      <td>200</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>11/26/85</td>\n",
       "      <td>6</td>\n",
       "      <td>76</td>\n",
       "      <td>200</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>1/12/86</td>\n",
       "      <td>6</td>\n",
       "      <td>58</td>\n",
       "      <td>200</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Date  Count  Temperature  Pressure  Malfunction\n",
       "0     4/12/81      6           66        50            0\n",
       "1    11/12/81      6           70        50            1\n",
       "2     3/22/82      6           69        50            0\n",
       "3    11/11/82      6           68        50            0\n",
       "4     4/04/83      6           67        50            0\n",
       "5     6/18/82      6           72        50            0\n",
       "6     8/30/83      6           73       100            0\n",
       "7    11/28/83      6           70       100            0\n",
       "8     2/03/84      6           57       200            1\n",
       "9     4/06/84      6           63       200            1\n",
       "10    8/30/84      6           70       200            1\n",
       "11   10/05/84      6           78       200            0\n",
       "12   11/08/84      6           67       200            0\n",
       "13    1/24/85      6           53       200            2\n",
       "14    4/12/85      6           67       200            0\n",
       "15    4/29/85      6           75       200            0\n",
       "16    6/17/85      6           70       200            0\n",
       "17  7/2903/85      6           81       200            0\n",
       "18    8/27/85      6           76       200            0\n",
       "19   10/03/85      6           79       200            0\n",
       "20   10/30/85      6           75       200            2\n",
       "21   11/26/85      6           76       200            0\n",
       "22    1/12/86      6           58       200            1"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data = pd.read_csv(\"./shuttle.csv\")\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We know from our previous experience on this data set that filtering data is a really bad idea. We will therefore process it as such."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mlinger@drakai.local/.local/lib/python3.10/site-packages/pandas/plotting/_matplotlib/core.py:1114: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored\n",
      "  scatter = ax.scatter(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxdElEQVR4nO3de1xVdb7/8fcGuUjeRUBNhdTU0rwmP1K7nFSsxjRPjanlJfM8LB0dyTKalMxTlE2M1VieLl6avOUcj3XGspBkMmXGvKA1Y6h4oQwQbyGQsIP1+6Nxn7agwl5b9ubr6/l47Eeu71rfvT77Ay7frbX23g7LsiwBAAAYIsDXBQAAAHgT4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGMWn4ebzzz/X0KFD1apVKzkcDq1bt+6Sc9LT09WrVy+FhISoQ4cOWrp06WWvEwAA1B0+DTfFxcXq3r27Fi5cWK3tDx06pLvuuku33XabMjMz9dvf/lYPP/ywPvnkk8tcKQAAqCsc/vLFmQ6HQ//zP/+j4cOHX3CbWbNmaf369fr6669dY/fff79Onz6tDRs21EKVAADA39XzdQE1kZGRoYEDB7qNxcfH67e//e0F55SWlqq0tNS1XFFRoZMnT6p58+ZyOByXq1QAAOBFlmXpzJkzatWqlQICLn7hqU6Fm7y8PEVGRrqNRUZGqrCwUD/++KPq169faU5ycrLmzp1bWyUCAIDL6Ntvv9XVV1990W3qVLjxRGJiohISElzLP/zwg9q2batDhw6pYcOGPqzM+5xOpzZt2qTbbrtNQUFBvi6nTqKH9tA/++ihPfTPPn/t4ZkzZxQTE1Otf7vrVLiJiopSfn6+21h+fr4aNWpU5VkbSQoJCVFISEil8WbNmqlRo0aXpU5fcTqdCgsLU/Pmzf3qF7IuoYf20D/76KE99M8+f+3huVqqc0tJnfqcm7i4OKWlpbmNpaamKi4uzkcVAQAAf+PTcFNUVKTMzExlZmZK+vmt3pmZmcrJyZH08yWlsWPHurafPHmyDh48qCeeeELffPONXn/9db3//vuaMWOGL8oHAAB+yKfhZvv27erZs6d69uwpSUpISFDPnj01Z84cSVJubq4r6EhSTEyM1q9fr9TUVHXv3l0vv/yy3n77bcXHx/ukfgAA4H98es/Nrbfeqot9zE5Vnz586623ateuXZexKgAAUJfVqXtuAAAALoVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUn4ebhQsXKjo6WqGhoYqNjdW2bdsuuv2CBQvUqVMn1a9fX23atNGMGTN09uzZWqoWAAD4O5+Gm9WrVyshIUFJSUnauXOnunfvrvj4eB07dqzK7VesWKEnn3xSSUlJ2rt3r9555x2tXr1aTz31VC1XDgAA/JVPw01KSoomTZqkCRMm6LrrrtOiRYsUFhamxYsXV7n91q1b1a9fP40ePVrR0dEaPHiwRo0adcmzPQAA4MpRz1c7Lisr044dO5SYmOgaCwgI0MCBA5WRkVHlnJtuuknvvfeetm3bpr59++rgwYP66KOP9OCDD15wP6WlpSotLXUtFxYWSpKcTqecTqeXXo1/OPd6THtdtYke2kP/7KOH9tA/+/y1hzWpx2fh5vjx4yovL1dkZKTbeGRkpL755psq54wePVrHjx9X//79ZVmWfvrpJ02ePPmil6WSk5M1d+7cSuOffvqpwsLC7L0IP5WamurrEuo8emgP/bOPHtpD/+zztx6WlJRUe1ufhRtPpKen6/nnn9frr7+u2NhYHThwQNOnT9e8efM0e/bsKuckJiYqISHBtVxYWKg2bdpo8ODBatSoUW2VXiucTqdSU1M1aNAgBQUF+bqcOoke2kP/7KOH9tA/+/y1h+euvFSHz8JNeHi4AgMDlZ+f7zaen5+vqKioKufMnj1bDz74oB5++GFJUrdu3VRcXKz/+I//0O9+9zsFBFS+hSgkJEQhISGVxoOCgvzqh+ZNJr+22kIP7aF/9tFDe+ifff7Ww5rU4rMbioODg9W7d2+lpaW5xioqKpSWlqa4uLgq55SUlFQKMIGBgZIky7IuX7EAAKDO8OllqYSEBI0bN059+vRR3759tWDBAhUXF2vChAmSpLFjx6p169ZKTk6WJA0dOlQpKSnq2bOn67LU7NmzNXToUFfIAQAAVzafhpuRI0eqoKBAc+bMUV5ennr06KENGza4bjLOyclxO1Pz9NNPy+Fw6Omnn9bRo0fVokULDR06VM8995yvXgIAAPAzPr+heOrUqZo6dWqV69LT092W69Wrp6SkJCUlJdVCZQAAoC7y+dcvAAAAeBPhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAoPg83CxcuVHR0tEJDQxUbG6tt27ZddPvTp09rypQpatmypUJCQnTttdfqo48+qqVqAQCAv6vny52vXr1aCQkJWrRokWJjY7VgwQLFx8crKytLERERlbYvKyvToEGDFBERoT//+c9q3bq1jhw5oiZNmtR+8QAAwC/5NNykpKRo0qRJmjBhgiRp0aJFWr9+vRYvXqwnn3yy0vaLFy/WyZMntXXrVgUFBUmSoqOja7NkAADg53wWbsrKyrRjxw4lJia6xgICAjRw4EBlZGRUOefDDz9UXFycpkyZog8++EAtWrTQ6NGjNWvWLAUGBlY5p7S0VKWlpa7lwsJCSZLT6ZTT6fTiK/K9c6/HtNdVm+ihPfTPPnpoD/2zz197WJN6fBZujh8/rvLyckVGRrqNR0ZG6ptvvqlyzsGDB/XZZ59pzJgx+uijj3TgwAE9+uijcjqdSkpKqnJOcnKy5s6dW2n8008/VVhYmP0X4odSU1N9XUKdRw/toX/20UN76J99/tbDkpKSam/r08tSNVVRUaGIiAi9+eabCgwMVO/evXX06FG99NJLFww3iYmJSkhIcC0XFhaqTZs2Gjx4sBo1alRbpdcKp9Op1NRUDRo0yHXZDjVDD+2hf/bRQ3von33+2sNzV16qw2fhJjw8XIGBgcrPz3cbz8/PV1RUVJVzWrZsqaCgILdLUF26dFFeXp7KysoUHBxcaU5ISIhCQkIqjQcFBfnVD82bTH5ttYUe2kP/7KOH9tA/+/ythzWpxWdvBQ8ODlbv3r2VlpbmGquoqFBaWpri4uKqnNOvXz8dOHBAFRUVrrF9+/apZcuWVQYbAABw5fEo3Bw8eNArO09ISNBbb72lZcuWae/evXrkkUdUXFzsevfU2LFj3W44fuSRR3Ty5ElNnz5d+/bt0/r16/X8889rypQpXqkHAADUfR5dlurQoYNuueUWTZw4Uffee69CQ0M92vnIkSNVUFCgOXPmKC8vTz169NCGDRtcNxnn5OQoIOD/8lebNm30ySefaMaMGbrhhhvUunVrTZ8+XbNmzfJo/wAAwDwehZudO3dqyZIlSkhI0NSpUzVy5EhNnDhRffv2rfFzTZ06VVOnTq1yXXp6eqWxuLg4/e1vf6vxfgAAwJXBo8tSPXr00CuvvKLvv/9eixcvVm5urvr376+uXbsqJSVFBQUF3q4TAACgWmzdUFyvXj2NGDFCa9as0YsvvqgDBw5o5syZatOmjcaOHavc3Fxv1QkAAFAttsLN9u3b9eijj6ply5ZKSUnRzJkzlZ2drdTUVH3//fcaNmyYt+oEAACoFo/uuUlJSdGSJUuUlZWlO++8U++++67uvPNO182/MTExWrp0Kd/7BAAAap1H4eaNN97QQw89pPHjx6tly5ZVbhMREaF33nnHVnEAAAA15VG42b9//yW3CQ4O1rhx4zx5egAAAI95dM/NkiVLtGbNmkrja9as0bJly2wXBQAA4CmPwk1ycrLCw8MrjUdEROj555+3XRQAAICnPAo3OTk5iomJqTTerl075eTk2C4KAADAUx6Fm4iICO3Zs6fS+O7du9W8eXPbRQEAAHjKo3AzatQoTZs2TZs2bVJ5ebnKy8v12Wefafr06br//vu9XSMAAEC1efRuqXnz5unw4cO6/fbbVa/ez09RUVGhsWPHcs8NAADwKY/CTXBwsFavXq158+Zp9+7dql+/vrp166Z27dp5uz4AAIAa8SjcnHPttdfq2muv9VYtAAAAtnkUbsrLy7V06VKlpaXp2LFjqqiocFv/2WefeaU4AACAmvIo3EyfPl1Lly7VXXfdpa5du8rhcHi7LgAAAI94FG5WrVql999/X3feeae36wEAALDFo7eCBwcHq0OHDt6uBQAAwDaPws1jjz2mV155RZZlebseAAAAWzy6LPXFF19o06ZN+vjjj3X99dcrKCjIbf3atWu9UhwAAEBNeRRumjRponvuucfbtQAAANjmUbhZsmSJt+sAAADwCo/uuZGkn376SRs3btR//dd/6cyZM5Kk77//XkVFRV4rDgAAoKY8OnNz5MgRDRkyRDk5OSotLdWgQYPUsGFDvfjiiyotLdWiRYu8XScAAEC1eHTmZvr06erTp49OnTql+vXru8bvuecepaWlea04AACAmvLozM3mzZu1detWBQcHu41HR0fr6NGjXikMAADAEx6duamoqFB5eXml8e+++04NGza0XRQAAICnPAo3gwcP1oIFC1zLDodDRUVFSkpK4isZAACAT3l0Werll19WfHy8rrvuOp09e1ajR4/W/v37FR4erpUrV3q7RgAAgGrzKNxcffXV2r17t1atWqU9e/aoqKhIEydO1JgxY9xuMAYAAKhtHoUbSapXr54eeOABb9YCAABgm0fh5t13373o+rFjx3pUDAAAgF0ehZvp06e7LTudTpWUlCg4OFhhYWGEGwAA4DMevVvq1KlTbo+ioiJlZWWpf//+3FAMAAB8yuPvljpfx44d9cILL1Q6qwMAAFCbvBZupJ9vMv7++++9+ZQAAAA14tE9Nx9++KHbsmVZys3N1R//+Ef169fPK4UBAAB4wqNwM3z4cLdlh8OhFi1a6N/+7d/08ssve6MuAAAAj3gUbioqKrxdBwAAgFd49Z4bAAAAX/PozE1CQkK1t01JSfFkFwAAAB7xKNzs2rVLu3btktPpVKdOnSRJ+/btU2BgoHr16uXazuFweKdKAACAavIo3AwdOlQNGzbUsmXL1LRpU0k/f7DfhAkTNGDAAD322GNeLRIAAKC6PLrn5uWXX1ZycrIr2EhS06ZN9Z//+Z+8WwoAAPiUR+GmsLBQBQUFlcYLCgp05swZ20UBAAB4yqNwc88992jChAlau3atvvvuO3333Xf67//+b02cOFEjRozwdo0AAADV5tE9N4sWLdLMmTM1evRoOZ3On5+oXj1NnDhRL730klcLBAAAqAmPwk1YWJhef/11vfTSS8rOzpYktW/fXldddZVXiwMAAKgpWx/il5ubq9zcXHXs2FFXXXWVLMvyVl0AAAAe8SjcnDhxQrfffruuvfZa3XnnncrNzZUkTZw4kbeBAwAAn/Io3MyYMUNBQUHKyclRWFiYa3zkyJHasGGD14oDAACoKY/uufn000/1ySef6Oqrr3Yb79ixo44cOeKVwgAAADzh0Zmb4uJitzM255w8eVIhISG2iwIAAPCUR+FmwIABevfdd13LDodDFRUVmj9/vm677TavFQcAAFBTHl2Wmj9/vm6//XZt375dZWVleuKJJ/SPf/xDJ0+e1JYtW7xdIwAAQLV5dOama9eu2rdvn/r3769hw4apuLhYI0aM0K5du9S+fXtv1wgAAFBtNT5z43Q6NWTIEC1atEi/+93vLkdNAAAAHqvxmZugoCDt2bPnctQCAABgm0eXpR544AG988473q4FAADANo9uKP7pp5+0ePFibdy4Ub179670nVIpKSleKQ4AAKCmahRuDh48qOjoaH399dfq1auXJGnfvn1u2zgcDu9VBwAAUEM1CjcdO3ZUbm6uNm3aJOnnr1t49dVXFRkZeVmKAwAAqKka3XNz/rd+f/zxxyouLvZqQQAAAHZ4dEPxOeeHHU8tXLhQ0dHRCg0NVWxsrLZt21ateatWrZLD4dDw4cO9UgcAAKj7ahRuHA5HpXtq7N5js3r1aiUkJCgpKUk7d+5U9+7dFR8fr2PHjl103uHDhzVz5kwNGDDA1v4BAIBZanTPjWVZGj9+vOvLMc+ePavJkydXerfU2rVrq/2cKSkpmjRpkiZMmCBJWrRokdavX6/FixfrySefrHJOeXm5xowZo7lz52rz5s06ffp0TV4GAAAwWI3Czbhx49yWH3jgAVs7Lysr044dO5SYmOgaCwgI0MCBA5WRkXHBec8++6wiIiI0ceJEbd68+aL7KC0tVWlpqWu5sLBQ0s+ftOx0Om3V72/OvR7TXldtoof20D/76KE99M8+f+1hTeqpUbhZsmRJjYu5mOPHj6u8vLzSu60iIyP1zTffVDnniy++0DvvvKPMzMxq7SM5OVlz586tNP7pp58qLCysxjXXBampqb4uoc6jh/bQP/vooT30zz5/62FJSUm1t/XoQ/x85cyZM3rwwQf11ltvKTw8vFpzEhMTlZCQ4FouLCxUmzZtNHjwYDVq1OhyleoTTqdTqampGjRokIKCgnxdTp1ED+2hf/bRQ3von33+2sNzV16qw6fhJjw8XIGBgcrPz3cbz8/PV1RUVKXts7OzdfjwYQ0dOtQ1VlFRIUmqV6+esrKyKn0reUhIiOseoV8KCgryqx+aN5n82moLPbSH/tlHD+2hf/b5Ww9rUoutt4LbFRwcrN69eystLc01VlFRobS0NMXFxVXavnPnzvrqq6+UmZnpetx999267bbblJmZqTZt2tRm+QAAwA/5/LJUQkKCxo0bpz59+qhv375asGCBiouLXe+eGjt2rFq3bq3k5GSFhoaqa9eubvObNGkiSZXGAQDAlcnn4WbkyJEqKCjQnDlzlJeXpx49emjDhg2um4xzcnIUEODTE0wAAKAO8Xm4kaSpU6dq6tSpVa5LT0+/6NylS5d6vyAAAFBncUoEAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGy86WFCkTVnHdOh4sa9LAYAr0uF/HX+PnCjxcSXwpXq+LsAEp0vKNG1lpj7fX+Aau7ljC702qqcahwX5sDIAuDKcOw7//eAxze8r3fXaZsVeE8Fx+ArFmRsvmLYyU1sOHHcb23LguH6zcpePKgKAKwvHYfwS4camgwVF+nx/gcoty2283LL0+f4CLlEBwGXGcRjnI9zYdOTkxa/rHj7BXyoAuJw4DuN8hBub2jULu+j66OZX1VIlAHBl4jiM8xFubLqmRQPd3LGFAh0Ot/FAh0M3d2yhmHD+UgHA5cRxGOcj3HjBa6N6ql+HcLexfh3C9dqonj6qCACuLByH8Uu8FdwLGocF6d2JfXXoeLEOnyhWdPOr+D8FAKhF547DB/J+0D/+nq71vxmgDlGNfV0WfIRw40Ux4YQaAPClds3D9I9//RdXLi5LAQAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIziF+Fm4cKFio6OVmhoqGJjY7Vt27YLbvvWW29pwIABatq0qZo2baqBAwdedHsAAHBl8Xm4Wb16tRISEpSUlKSdO3eqe/fuio+P17Fjx6rcPj09XaNGjdKmTZuUkZGhNm3aaPDgwTp69GgtVw4AAPyRz8NNSkqKJk2apAkTJui6667TokWLFBYWpsWLF1e5/fLly/Xoo4+qR48e6ty5s95++21VVFQoLS2tlisHAAD+qJ4vd15WVqYdO3YoMTHRNRYQEKCBAwcqIyOjWs9RUlIip9OpZs2aVbm+tLRUpaWlruXCwkJJktPplNPptFG9/zn3ekx7XbWJHtpD/+yjh/bQP/v8tYc1qcen4eb48eMqLy9XZGSk23hkZKS++eabaj3HrFmz1KpVKw0cOLDK9cnJyZo7d26l8U8//VRhYWE1L7oOSE1N9XUJdR49tIf+2UcP7aF/9vlbD0tKSqq9rU/DjV0vvPCCVq1apfT0dIWGhla5TWJiohISElzLhYWFrvt0GjVqVFul1gqn06nU1FQNGjRIQUFBvi6nTqKH9tA/++ihPfTPPn/t4bkrL9Xh03ATHh6uwMBA5efnu43n5+crKirqonN///vf64UXXtDGjRt1ww03XHC7kJAQhYSEVBoPCgryqx+aN5n82moLPbSH/tlHD+2hf/b5Ww9rUotPbygODg5W79693W4GPndzcFxc3AXnzZ8/X/PmzdOGDRvUp0+f2igVAADUET6/LJWQkKBx48apT58+6tu3rxYsWKDi4mJNmDBBkjR27Fi1bt1aycnJkqQXX3xRc+bM0YoVKxQdHa28vDxJUoMGDdSgQQOfvQ4AAOAffB5uRo4cqYKCAs2ZM0d5eXnq0aOHNmzY4LrJOCcnRwEB/3eC6Y033lBZWZnuvfdet+dJSkrSM888U5ulAwAAP+TzcCNJU6dO1dSpU6tcl56e7rZ8+PDhy18QAACos3z+IX4AAADeRLiB1xwsKNKmrGM6dLzYJ/NBD02wZX+BJCkj+7iPKwHqLr+4LIW67XRJmaatzNTn/zooS9LNHVvotVE91Tjs0m/dszsf9NAER04Ua/jCLSopLdP8vtKkP+1QWEiwPpzSX22am/mBo8Dlwpkb2DZtZaa2HHD/v8wtB47rNyt31cp80EMTDF+4RadK3D9e/lSJU3cv/MJHFQF1F+EGthwsKNLn+wtUbllu4+WWpc/3F1zy8ojd+aCHJvhr1rFKweacUyVObf7FGTkAl0a4gS1HTl78uz4On7j4P6x254MemiDzu9MXXb8z51TtFAIYgnADW9o1u/i9ANHNr7qs80EPTdDj6iYXXd+rbdPaKQQwBOEGtlzTooFu7thCgQ6H23igw6GbO7ZQTPjF/2G1Ox/00AS3dIpQ0wvc+N00LEgDOrao5YqAuo1wA9teG9VT/TqEu4316xCu10b1rJX5oIcm+HBK/0oBp2lYkD6c0t9HFQF1F28Fh22Nw4L07sS+OnS8WIdPFCu6+VU1Oltgdz7ooQnaNA/TrjmD9fk3uTqVtU1vPdhbN3du6euygDqJcAOviQm39w+q3fmghyaIax+uj7J+/i8Az3BZCgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAo/hFuFm4cKGio6MVGhqq2NhYbdu27aLbr1mzRp07d1ZoaKi6deumjz76qJYqBXA5HSwo0qasYzp0vNij+au35ei3q3dpzfZvfbJ/u/Mlacv+AklSRvZxj5/DDn/ogR2H/7XfIydKPJr/16xjeiVtnzb/6+dQ27zRP1//DHy9f0mq57M9/8vq1auVkJCgRYsWKTY2VgsWLFB8fLyysrIUERFRafutW7dq1KhRSk5O1q9+9SutWLFCw4cP186dO9W1a1cfvAIAdp0uKdO0lZn6/Bf/oNzcsYVeG9VTjcOCLjn/q+9O657Xt+qnCkuStG7X90pc+5U+nNJP17VufNn3b3e+JB05UazhC7eopLRM8/tKk/60Q2EhwfpwSn+1aR5Wreewwx96YMe5/f/94DHN7yvd9dpmxV4TUe39n+v/qRKna6xpWFCd6b+3nsMOX+//l3x+5iYlJUWTJk3ShAkTdN1112nRokUKCwvT4sWLq9z+lVde0ZAhQ/T444+rS5cumjdvnnr16qU//vGPtVw5AG+ZtjJTWw64n6nYcuC4frNyV7Xm/zLYnPNThaW7F26plf3bnS+p0j+sknSqxKm7F35R7eewwx96YIfd/df1/nvrOezw9f5/yadnbsrKyrRjxw4lJia6xgICAjRw4EBlZGRUOScjI0MJCQluY/Hx8Vq3bl2V25eWlqq0tNS1/MMPP0iSTp48KafTWeWcusrpdKqkpEQnTpxQUFDtpmRT0EN7POlfzolibd17WA5VPiBt3VukXftbqm2zC/+f84eZR6WyogsezJZ99pV+1b3VZdu/3fmS9Lfs4zpTeFr1JNWrsFRSUqF6zgCVVzh0plD6ePs+9Y1pftHnsMMfemCH2/4D3PtXnf3/sv/nqwv999ZznOOLv8fVcebMGUmSZVmX2PLnjXzm6NGjliRr69atbuOPP/641bdv3yrnBAUFWStWrHAbW7hwoRUREVHl9klJSZYkHjx48ODBg4cBj2+//faS+cLn99xcbomJiW5neioqKnTy5Ek1b95cDofDh5V5X2Fhodq0aaNvv/1WjRo18nU5dRI9tIf+2UcP7aF/9vlrDy3L0pkzZ9Sq1YXPxJ7j03ATHh6uwMBA5efnu43n5+crKiqqyjlRUVE12j4kJEQhISFuY02aNPG86DqgUaNGfvULWRfRQ3von3300B76Z58/9rBx48bV2s6nNxQHBwerd+/eSktLc41VVFQoLS1NcXFxVc6Ji4tz216SUlNTL7g9AAC4svj8slRCQoLGjRunPn36qG/fvlqwYIGKi4s1YcIESdLYsWPVunVrJScnS5KmT5+uW265RS+//LLuuusurVq1Stu3b9ebb77py5cBAAD8hM/DzciRI1VQUKA5c+YoLy9PPXr00IYNGxQZGSlJysnJUUDA/51guummm7RixQo9/fTTeuqpp9SxY0etW7eOz7jRz5fgkpKSKl2GQ/XRQ3von3300B76Z58JPXRYVnXeUwUAAFA3+PxD/AAAALyJcAMAAIxCuAEAAEYh3AAAAKMQbuqYZ555Rg6Hw+3RuXNn1/pbb7210vrJkyf7sGL/dPToUT3wwANq3ry56tevr27dumn79u2u9ZZlac6cOWrZsqXq16+vgQMHav/+/T6s2L9cqn/jx4+v9Hs4ZMgQH1bsX6Kjoyv1x+FwaMqUKZKks2fPasqUKWrevLkaNGigf//3f6/04aVXskv1j+PgpZWXl2v27NmKiYlR/fr11b59e82bN8/te5vq8nHQ528FR81df/312rhxo2u5Xj33H+OkSZP07LPPupbDwi7fF9bVRadOnVK/fv1022236eOPP1aLFi20f/9+NW3a1LXN/Pnz9eqrr2rZsmWKiYnR7NmzFR8fr3/+858KDQ31YfW+V53+SdKQIUO0ZMkS13Jdflupt3355ZcqLy93LX/99dcaNGiQ7rvvPknSjBkztH79eq1Zs0aNGzfW1KlTNWLECG3ZUr1vOTfdpfoncRy8lBdffFFvvPGGli1bpuuvv17bt2/XhAkT1LhxY02bNk1SHT8OXvLbp+BXkpKSrO7du19w/S233GJNnz691uqpi2bNmmX179//gusrKiqsqKgo66WXXnKNnT592goJCbFWrlxZGyX6tUv1z7Isa9y4cdawYcNqpyADTJ8+3Wrfvr1VUVFhnT592goKCrLWrFnjWr93715LkpWRkeHDKv3XL/tnWRwHq+Ouu+6yHnroIbexESNGWGPGjLEsq+4fB7ksVQft379frVq10jXXXKMxY8YoJyfHbf3y5csVHh6url27KjExUSUlJT6q1D99+OGH6tOnj+677z5FRESoZ8+eeuutt1zrDx06pLy8PA0cONA11rhxY8XGxiojI8MXJfuVS/XvnPT0dEVERKhTp0565JFHdOLECR9U6//Kysr03nvv6aGHHpLD4dCOHTvkdDrdfv86d+6stm3b8vtXhfP7dw7HwYu76aablJaWpn379kmSdu/erS+++EJ33HGHpLp/HOSyVB0TGxurpUuXqlOnTsrNzdXcuXM1YMAAff3112rYsKFGjx6tdu3aqVWrVtqzZ49mzZqlrKwsrV271tel+42DBw/qjTfeUEJCgp566il9+eWXmjZtmoKDgzVu3Djl5eVJkutTss+JjIx0rbuSXap/0s+XpEaMGKGYmBhlZ2frqaee0h133KGMjAwFBgb6+BX4l3Xr1un06dMaP368JCkvL0/BwcGVvuCX37+qnd8/SRwHq+HJJ59UYWGhOnfurMDAQJWXl+u5557TmDFjJKnuHwd9feoI9pw6dcpq1KiR9fbbb1e5Pi0tzZJkHThwoJYr819BQUFWXFyc29hvfvMb6//9v/9nWZZlbdmyxZJkff/9927b3Hfffdavf/3rWqvTX12qf1XJzs62JFkbN2683OXVOYMHD7Z+9atfuZaXL19uBQcHV9ruxhtvtJ544onaLK1OOL9/VeE4WNnKlSutq6++2lq5cqW1Z88e691337WaNWtmLV261LKsun8c5LJUHdekSRNde+21OnDgQJXrY2NjJemC669ELVu21HXXXec21qVLF9flvaioKEmq9O6U/Px817or2aX6V5VrrrlG4eHh/B6e58iRI9q4caMefvhh11hUVJTKysp0+vRpt235/ausqv5VheNgZY8//riefPJJ3X///erWrZsefPBBzZgxw/Ul1XX9OEi4qeOKioqUnZ2tli1bVrk+MzNTki64/krUr18/ZWVluY3t27dP7dq1kyTFxMQoKipKaWlprvWFhYX6+9//rri4uFqt1R9dqn9V+e6773TixAl+D8+zZMkSRURE6K677nKN9e7dW0FBQW6/f1lZWcrJyeH37zxV9a8qHAcrKykpcftSakkKDAxURUWFJAOOg74+dYSaeeyxx6z09HTr0KFD1pYtW6yBAwda4eHh1rFjx6wDBw5Yzz77rLV9+3br0KFD1gcffGBdc8011s033+zrsv3Ktm3brHr16lnPPfectX//fmv58uVWWFiY9d5777m2eeGFF6wmTZpYH3zwgbVnzx5r2LBhVkxMjPXjjz/6sHL/cKn+nTlzxpo5c6aVkZFhHTp0yNq4caPVq1cvq2PHjtbZs2d9XL3/KC8vt9q2bWvNmjWr0rrJkydbbdu2tT777DNr+/btVlxcXKVLgVe6C/WP42D1jBs3zmrdurX1l7/8xTp06JC1du1aKzw83O3SZ10+DhJu6piRI0daLVu2tIKDg63WrVtbI0eOdF1HzsnJsW6++WarWbNmVkhIiNWhQwfr8ccft3744QcfV+1//vd//9fq2rWrFRISYnXu3Nl688033dZXVFRYs2fPtiIjI62QkBDr9ttvt7KysnxUrf+5WP9KSkqswYMHWy1atLCCgoKsdu3aWZMmTbLy8vJ8WLH/+eSTTyxJVf5e/fjjj9ajjz5qNW3a1AoLC7PuueceKzc31wdV+q8L9Y/jYPUUFhZa06dPt9q2bWuFhoZa11xzjfW73/3OKi0tdW1Tl4+DDsv6xccRAgAA1HHccwMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwA+CiHA7HRR/PPPOMr0v0uujoaC1YsMDXZQDwUD1fFwDAv+Xm5rr+vHr1as2ZM8ftizMbNGjgi7JqzLIslZeXq1692jvslZWVKTg4uNb2B+BnnLkBcFFRUVGuR+PGjeVwONzGVq1apS5duig0NFSdO3fW66+/7pp7+PBhORwOvf/++xowYIDq16+vG2+8Ufv27dOXX36pPn36qEGDBrrjjjtUUFDgmjd+/HgNHz5cc+fOVYsWLdSoUSNNnjxZZWVlrm0qKiqUnJysmJgY1a9fX927d9ef//xn1/r09HQ5HA59/PHH6t27t0JCQvTFF18oOztbw4YNU2RkpBo0aKAbb7xRGzdudM279dZbdeTIEc2YMcN1dkqSnnnmGfXo0cOtNwsWLFB0dHSlup977jm1atVKnTp1kiR9++23+vWvf60mTZqoWbNmGjZsmA4fPuyNHw+AKhBuAHhs+fLlmjNnjp577jnt3btXzz//vGbPnq1ly5a5bZeUlKSnn35aO3fuVL169TR69Gg98cQTeuWVV7R582YdOHBAc+bMcZuTlpamvXv3Kj09XStXrtTatWs1d+5c1/rk5GS9++67WrRokf7xj39oxowZeuCBB/TXv/7V7XmefPJJvfDCC9q7d69uuOEGFRUV6c4771RaWpp27dqlIUOGaOjQocrJyZEkrV27VldffbWeffZZ5ebmup25qo60tDRlZWUpNTVVf/nLX+R0OhUfH6+GDRtq8+bN2rJlixo0aKAhQ4a4hTUAXuTjL+4EUIcsWbLEaty4sWu5ffv21ooVK9y2mTdvnhUXF2dZlmUdOnTIkmS9/fbbrvUrV660JFlpaWmuseTkZKtTp06u5XHjxlnNmjWziouLXWNvvPGG1aBBA6u8vNw6e/asFRYWZm3dutVt3xMnTrRGjRplWZZlbdq0yZJkrVu37pKv6/rrr7dee+0113K7du2sP/zhD27bJCUlWd27d3cb+8Mf/mC1a9fOre7IyEi3b1b+05/+ZHXq1MmqqKhwjZWWllr169e3Pvnkk0vWBqDmuOcGgEeKi4uVnZ2tiRMnatKkSa7xn376SY0bN3bb9oYbbnD9OTIyUpLUrVs3t7Fjx465zenevbvCwsJcy3FxcSoqKtK3336roqIilZSUaNCgQW5zysrK1LNnT7exPn36uC0XFRXpmWee0fr165Wbm6uffvpJP/74o+vMjV3dunVzu89m9+7dOnDggBo2bOi23dmzZ5Wdne2VfQJwR7gB4JGioiJJ0ltvvaXY2Fi3dYGBgW7LQUFBrj+fu4fl/LGKiooa73v9+vVq3bq127qQkBC35auuuspteebMmUpNTdXvf/97dejQQfXr19e99957yUtEAQEBsizLbczpdFba7vz9FRUVqXfv3lq+fHmlbVu0aHHRfQLwDOEGgEciIyPVqlUrHTx4UGPGjPH68+/evVs//vij6tevL0n629/+pgYNGqhNmzZq1qyZQkJClJOTo1tuuaVGz7tlyxaNHz9e99xzj6Sfw8f5N/cGBwervLzcbaxFixbKy8uTZVmugJaZmXnJ/fXq1UurV69WRESEGjVqVKNaAXiGG4oBeGzu3LlKTk7Wq6++qn379umrr77SkiVLlJKSYvu5y8rKNHHiRP3zn//URx99pKSkJE2dOlUBAQFq2LChZs6cqRkzZmjZsmXKzs7Wzp079dprr1W6mfl8HTt21Nq1a5WZmandu3dr9OjRlc4aRUdH6/PPP9fRo0d1/PhxST+/i6qgoEDz589Xdna2Fi5cqI8//viSr2PMmDEKDw/XsGHDtHnzZh06dEjp6emaNm2avvvuO88bBOCCCDcAPPbwww/r7bff1pIlS9StWzfdcsstWrp0qWJiYmw/9+23366OHTvq5ptv1siRI3X33Xe7fWDgvHnzNHv2bCUnJ6tLly4aMmSI1q9ff8l9p6SkqGnTprrppps0dOhQxcfHq1evXm7bPPvsszp8+LDat2/vunTUpUsXvf7661q4cKG6d++ubdu2aebMmZd8HWFhYfr888/Vtm1bjRgxQl26dNHEiRN19uxZzuQAl4nDOv8iMgD42Pjx43X69GmtW7fO16UAqIM4cwMAAIxCuAEAAEbhshQAADAKZ24AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFH+PxzP5o/D94vEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "pd.set_option('mode.chained_assignment',None) # this removes a useless warning from pandas\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "data[\"Frequency\"]=data.Malfunction/data.Count\n",
    "data.plot(x=\"Temperature\",y=\"Frequency\",kind=\"scatter\",ylim=[0,1])\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "data[\"Intercept\"] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic regression\n",
    "\n",
    "Let's assume O-rings independently fail with the same probability which solely depends on temperature. A logistic regression should allow us to estimate the influence of temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/mlinger@drakai.local/.local/lib/python3.10/site-packages/statsmodels/genmod/families/links.py:13: FutureWarning: The logit link alias is deprecated. Use Logit instead. The logit link alias will be removed after the 0.15.0 release.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>       <td>Frequency</td>    <th>  No. Observations:  </th>  <td>    23</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>    21</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>        <td>Binomial</td>     <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>         <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -23.526</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>            <td>Tue, 09 Jul 2024</td> <th>  Deviance:          </th> <td>  18.086</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                <td>15:31:00</td>     <th>  Pearson chi2:      </th>  <td>  30.0</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>          <td>6</td>        <th>  Pseudo R-squ. (CS):</th>  <td>0.2344</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "       <td></td>          <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th>   <td>    5.0850</td> <td>    3.052</td> <td>    1.666</td> <td> 0.096</td> <td>   -0.898</td> <td>   11.068</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Temperature</th> <td>   -0.1156</td> <td>    0.047</td> <td>   -2.458</td> <td> 0.014</td> <td>   -0.208</td> <td>   -0.023</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/latex": [
       "\\begin{center}\n",
       "\\begin{tabular}{lclc}\n",
       "\\toprule\n",
       "\\textbf{Dep. Variable:}   &    Frequency     & \\textbf{  No. Observations:  } &       23    \\\\\n",
       "\\textbf{Model:}           &       GLM        & \\textbf{  Df Residuals:      } &       21    \\\\\n",
       "\\textbf{Model Family:}    &     Binomial     & \\textbf{  Df Model:          } &        1    \\\\\n",
       "\\textbf{Link Function:}   &      logit       & \\textbf{  Scale:             } &    1.0000   \\\\\n",
       "\\textbf{Method:}          &       IRLS       & \\textbf{  Log-Likelihood:    } &   -23.526   \\\\\n",
       "\\textbf{Date:}            & Tue, 09 Jul 2024 & \\textbf{  Deviance:          } &    18.086   \\\\\n",
       "\\textbf{Time:}            &     15:31:00     & \\textbf{  Pearson chi2:      } &     30.0    \\\\\n",
       "\\textbf{No. Iterations:}  &        6         & \\textbf{  Pseudo R-squ. (CS):} &   0.2344    \\\\\n",
       "\\textbf{Covariance Type:} &    nonrobust     & \\textbf{                     } &             \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "\\begin{tabular}{lcccccc}\n",
       "                     & \\textbf{coef} & \\textbf{std err} & \\textbf{z} & \\textbf{P$> |$z$|$} & \\textbf{[0.025} & \\textbf{0.975]}  \\\\\n",
       "\\midrule\n",
       "\\textbf{Intercept}   &       5.0850  &        3.052     &     1.666  &         0.096        &       -0.898    &       11.068     \\\\\n",
       "\\textbf{Temperature} &      -0.1156  &        0.047     &    -2.458  &         0.014        &       -0.208    &       -0.023     \\\\\n",
       "\\bottomrule\n",
       "\\end{tabular}\n",
       "%\\caption{Generalized Linear Model Regression Results}\n",
       "\\end{center}"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:              Frequency   No. Observations:                   23\n",
       "Model:                            GLM   Df Residuals:                       21\n",
       "Model Family:                Binomial   Df Model:                            1\n",
       "Link Function:                  logit   Scale:                          1.0000\n",
       "Method:                          IRLS   Log-Likelihood:                -23.526\n",
       "Date:                Tue, 09 Jul 2024   Deviance:                       18.086\n",
       "Time:                        15:31:00   Pearson chi2:                     30.0\n",
       "No. Iterations:                     6   Pseudo R-squ. (CS):             0.2344\n",
       "Covariance Type:            nonrobust                                         \n",
       "===============================================================================\n",
       "                  coef    std err          z      P>|z|      [0.025      0.975]\n",
       "-------------------------------------------------------------------------------\n",
       "Intercept       5.0850      3.052      1.666      0.096      -0.898      11.068\n",
       "Temperature    -0.1156      0.047     -2.458      0.014      -0.208      -0.023\n",
       "===============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "\n",
    "# Create an instance of the link class\n",
    "logit_link = sm.families.links.logit()\n",
    "\n",
    "# Use the instance of the link class in the Binomial family\n",
    "logmodel = sm.GLM(data['Frequency'], data[['Intercept','Temperature']], \n",
    "                  family=sm.families.Binomial(link=logit_link),\n",
    "                  var_weights=data['Count']).fit()\n",
    "\n",
    "logmodel.summary()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Good, now I have recovered the asymptotic standard errors $s_{\\hat{\\alpha}}=3.052$ and $s_{\\hat{\\beta}}=0.047$.\n",
    "The Goodness of fit (Deviance) indicated for this model is $G^2=18.086$ with 21 degrees of freedom (Df Residuals).\n",
    "\n",
    "**I have therefore managed to fully replicate the results of the Dalal *et al.* article**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predicting failure probability\n",
    "The temperature when launching the shuttle was 31°F. Let's try to estimate the failure probability for such temperature using our model.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG2CAYAAACtaYbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAz2UlEQVR4nO3deXRV5b3/8c/JnBBCmDKAgSBSIAoBQpNfaFFbAkG5iO1tpYBMIveqsArkqhArhEgxjhStCFcU9IIottfSKhhMUwMoqRFoqBQKAkGQZsBSCBCSHHL27w9vjhwynZPpgfB+rXXWYj/72Xs/55sN+bBHm2VZlgAAAAzxMj0AAABwfSOMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKM8DiPbt2/X2LFj1a1bN9lsNm3atKnBZXJycjRkyBD5+/vrpptu0uuvv96IoQIAgLbI4zBy4cIFxcbGasWKFW71Lygo0JgxY/SDH/xA+fn5mjt3ru6//35t3brV48ECAIC2x9aUF+XZbDb97ne/0913311nn/nz52vz5s3at2+fs+1nP/uZzpw5o8zMzMZuGgAAtBE+Lb2B3NxcJSUlubQlJydr7ty5dS5TUVGhiooK57TD4dDp06fVuXNn2Wy2lhoqAABoRpZl6dy5c+rWrZu8vOo+GdPiYaSoqEjh4eEubeHh4SotLdXFixcVGBhYY5mMjAylp6e39NAAAEArOHHihG644YY657d4GGmM1NRUpaSkOKfPnj2rHj16qKCgQO3bt2+WbViWpXNlFdq2fbtuu/VW+fhelaW4alyyX6JWHqBe7qNW7qNW7qNW7quuVfKI2+Xn59es6z537px69erV4O/uFv8JRUREqLi42KWtuLhYISEhtR4VkSR/f3/5+/vXaO/UqZNCQkKabWwd7HZ1bB+kGyLD5Ovr22zrbYvs1Moj1Mt91Mp91Mp91Mp91bXq0qVLs9eqen0NXWLR4s8ZSUxMVHZ2tktbVlaWEhMTW3rTAADgGuBxGDl//rzy8/OVn58v6Ztbd/Pz83X8+HFJ35ximTJlirP/Aw88oKNHj+rRRx/V3//+d7388st65513NG/evOb5BgAA4JrmcRjZtWuXBg8erMGDB0uSUlJSNHjwYC1atEiSVFhY6AwmktSrVy9t3rxZWVlZio2N1fPPP69XX31VycnJzfQVAADAtczja0Zuv/121fdoktqernr77bfrL3/5i6ebAgC0EVVVVbLb7a22PbvdLh8fH5WXl6uqqqrVtnstakqtfH195e3t3eQxcIkxAKDFWJaloqIinTlzptW3GxERoRMnTvB8qgY0tVahoaGKiIhoUp0JIwCAFlMdRMLCwhQUFNRqwcDhcOj8+fMKDg6u92FbaHytLMtSWVmZSkpKJEmRkZGNHgNhBADQIqqqqpxBpHPnzq26bYfDocrKSgUEBBBGGtCUWlU/oqOkpERhYWGNPmXDTwgA0CKqrxEJCgoyPBK0pOqfb1OuCSKMAABaFNdstG3N8fMljAAAAKMIIwAAwCjCCAAAV5g2bZpsNluNz+HDh00PrU3ibhoAAGoxevRorV271qWta9euLtOVlZXN/qbb6xFHRgAAqIW/v78iIiJcPiNGjNDs2bM1d+5cdenSxflqk3379umOO+5QcHCwwsPDNXnyZH399dfOdV24cEFTpkxRcHCwIiMj9fzzz+v222/X3LlznX1sNps2bdrkMobQ0FCXJ5ufOHFC99xzj0JDQ9WpUyeNGzdOx44dc86fNm2a7r77bj333HOKjIxU586dNWvWLJc7XSoqKjR//nxFRUXJ399f3/nOd7Ru3TpZlqWbbrpJzz33nMsY8vPzW/yoEGEEANBqLMtSWeWlVvlcrKxy/rm+15h46o033pCfn58++eQTrVq1SmfOnNEPf/hDDR48WLt27VJmZqaKi4t1zz33OJd55JFHtG3bNv3+97/Xhx9+qJycHO3Zs8ej7drtdiUnJ6t9+/basWOHPvnkEwUHB2v06NGqrKx09vvoo4905MgRffTRR3rjjTf0+uuvuwSaKVOm6K233tKLL76oAwcOaOXKlWrXrp1sNpvuu+++GkeD1q5dq1tvvVU33XRT4wrmBk7TAABazUV7lWIWbW317e5/IllBfp79ynv//fcVHBzsnL7jjjskSX369NEzzzzjbP/lL3+pwYMH68knn3S2rVmzRlFRUTp06JC6deum1157TevXr9eIESMkfRNobrjhBo/Gs3HjRjkcDr366qvO22nXrl2r0NBQ5eTkaNSoUZKkjh076qWXXpK3t7f69eunMWPGKDs7WzNnztShQ4f0zjvvKCsrS0lJSZKk6OholZaWSvrmyMqiRYuUl5en+Ph42e12bdiwocbRkuZGGAEAoBY/+MEPtHLlSud0u3btNGHCBMXFxbn027t3rz766COX4FLtyJEjunjxoiorK5WQkOBs79Spk/r27evRePbu3avDhw+rffv2Lu3l5eU6cuSIc/rmm292eRJqZGSkPv/8c0nfnHLx9vbWbbfdVus2unXrpjFjxmjNmjWKj4/Xe++9p4qKCv30pz/1aKyeIowAAFpNoK+39j+R3OLbcTgcOld6Tu1D2svLy0uBvp4/prxdu3a1nppo166dy/T58+c1duxYPf300zX6RkZGun2thc1mq3E66fJrPc6fP6+4uDi9+eabNZa9/MJaX1/fGut1OBySvn18e33uv/9+TZ48Wb/61a+0du1ajR8/vsWfoksYAQC0GpvN5vHpksZwOBy65OetID+fFn83zZAhQ/S///u/io6Olo9Pze/Wu3dv+fr66tNPP1WPHj0kSf/617906NAhlyMUXbt2VWFhoXP6iy++UFlZmct2Nm7cqLCwMIWEhDRqrAMGDJDD4dC2bducp2mudOedd6pdu3ZauXKlMjMztX379kZtyxNcwAoAQBPMmjVLp0+f1oQJE/TZZ5/pyJEj2rp1q6ZPn66qqioFBwdrxowZeuSRR/SnP/1J+/bt07Rp02qEpB/+8Id66aWX9Je//EW7du3SAw884HKUY9KkSerSpYvGjRunHTt2qKCgQDk5Ofr5z3+ur776yq2xRkdHa+rUqbrvvvu0adMm5zp+97vfOft4e3tr2rRpSk1NVZ8+fZSYmNg8haoHYQQAgCbo1q2bPvnkE1VVVWnUqFEaMGCA5s6dq9DQUGfgePbZZzV8+HCNHTtWSUlJ+v73v1/j2pPnn39eUVFRGj58uCZOnKiHH37Y5fRIUFCQtm/frh49eujHP/6x+vfvrxkzZqi8vNyjIyUrV67UT37yEz300EPq16+f/vM//9PlCIwkzZgxQ5WVlZo+fXoTKuM+TtMAAHCFy2+FvVxOTk6t7X369NG7775b5/qCg4O1bt06rVu3ztm2efNmlz7dunXT1q2udxqdOXPGZToiIkJvvPGGR+Nevny5y3RAQICWLVumZcuWSfrmlFb13TTVTp48KV9fX02ZMqXObTUnwggAAJD0zQPRTp06pcWLF+unP/2pwsPDW2W7nKYBAACSpLfeeks9e/bUmTNnXJ6l0tI4MgIAgAF1nfIxadq0aZo2bVqrb5cjIwAAwCjCCACgRTXne2Fw9WmOny9hBADQIqqfkXHlbaNoW6p/vlc++dUTXDMCAGgR3t7eCg0NVUlJiaRvnpNR/YK3luZwOFRZWany8vIWfwLrta6xtbIsS2VlZSopKVFoaKjL+3A8RRgBALSYiIgISXIGktZiWZYuXryowMDAVgtA16qm1io0NNT5c24swggAoMXYbDZFRkYqLCzM5aVvLc1ut2v79u269dZbm3T64HrQlFr5+vo26YhINcIIAKDFeXt7N8svLU+2d+nSJQUEBBBGGnA11IoTaQAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIxqVBhZsWKFoqOjFRAQoISEBOXl5dXbf/ny5erbt68CAwMVFRWlefPmqby8vFEDBgAAbYvHYWTjxo1KSUlRWlqa9uzZo9jYWCUnJ6ukpKTW/hs2bNCCBQuUlpamAwcO6LXXXtPGjRv12GOPNXnwAADg2udxGFm2bJlmzpyp6dOnKyYmRqtWrVJQUJDWrFlTa/+dO3fqe9/7niZOnKjo6GiNGjVKEyZMaPBoCgAAuD74eNK5srJSu3fvVmpqqrPNy8tLSUlJys3NrXWZYcOGaf369crLy1N8fLyOHj2qLVu2aPLkyXVup6KiQhUVFc7p0tJSSZLdbpfdbvdkyPWqXldzrrOtolaeoV7uo1buo1buo1bua8laubtOm2VZlrsr/cc//qHu3btr586dSkxMdLY/+uij2rZtmz799NNal3vxxRf18MMPy7IsXbp0SQ888IBWrlxZ53YWL16s9PT0Gu0bNmxQUFCQu8MFAAAGlZWVaeLEiTp79qxCQkLq7OfRkZHGyMnJ0ZNPPqmXX35ZCQkJOnz4sObMmaMlS5Zo4cKFtS6TmpqqlJQU53RpaamioqI0atSoer+Mp+x2u7KysjRy5Ej5+vo223rbImrlGerlPmrlPmrlPmrlvpasVfWZjYZ4FEa6dOkib29vFRcXu7QXFxcrIiKi1mUWLlyoyZMn6/7775ckDRgwQBcuXNB//Md/6Be/+IW8vGpetuLv7y9/f/8a7b6+vi2yU7XUetsiauUZ6uU+auU+auU+auW+lqiVu+vz6AJWPz8/xcXFKTs729nmcDiUnZ3tctrmcmVlZTUCh7e3tyTJgzNEAACgjfL4NE1KSoqmTp2qoUOHKj4+XsuXL9eFCxc0ffp0SdKUKVPUvXt3ZWRkSJLGjh2rZcuWafDgwc7TNAsXLtTYsWOdoQQAAFy/PA4j48eP16lTp7Ro0SIVFRVp0KBByszMVHh4uCTp+PHjLkdCHn/8cdlsNj3++OM6efKkunbtqrFjx2rp0qXN9y0AAMA1q1EXsM6ePVuzZ8+udV5OTo7rBnx8lJaWprS0tMZsCgAAtHG8mwYAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGNCiMrVqxQdHS0AgIClJCQoLy8vHr7nzlzRrNmzVJkZKT8/f31ne98R1u2bGnUgAEAQNvi4+kCGzduVEpKilatWqWEhAQtX75cycnJOnjwoMLCwmr0r6ys1MiRIxUWFqbf/va36t69u7788kuFhoY2x/gBAMA1zuMwsmzZMs2cOVPTp0+XJK1atUqbN2/WmjVrtGDBghr916xZo9OnT2vnzp3y9fWVJEVHRzdt1AAAoM3wKIxUVlZq9+7dSk1NdbZ5eXkpKSlJubm5tS7zhz/8QYmJiZo1a5Z+//vfq2vXrpo4caLmz58vb2/vWpepqKhQRUWFc7q0tFSSZLfbZbfbPRlyvarX1ZzrbKuolWeol/uolfuolfuolftaslburtOjMPL111+rqqpK4eHhLu3h4eH6+9//XusyR48e1Z/+9CdNmjRJW7Zs0eHDh/XQQw/JbrcrLS2t1mUyMjKUnp5eo/3DDz9UUFCQJ0N2S1ZWVrOvs62iVp6hXu6jVu6jVu6jVu5riVqVlZW51c/j0zSecjgcCgsL0yuvvCJvb2/FxcXp5MmTevbZZ+sMI6mpqUpJSXFOl5aWKioqSqNGjVJISEizjc1utysrK0sjR450nkJC7aiVZ6iX+6iV+6iV+6iV+1qyVtVnNhriURjp0qWLvL29VVxc7NJeXFysiIiIWpeJjIyUr6+vyymZ/v37q6ioSJWVlfLz86uxjL+/v/z9/Wu0+/r6tshO1VLrbYuolWeol/uolfuolfuolftaolburs+jW3v9/PwUFxen7OxsZ5vD4VB2drYSExNrXeZ73/ueDh8+LIfD4Ww7dOiQIiMjaw0iAADg+uLxc0ZSUlK0evVqvfHGGzpw4IAefPBBXbhwwXl3zZQpU1wucH3wwQd1+vRpzZkzR4cOHdLmzZv15JNPatasWc33LQAAwDXL42tGxo8fr1OnTmnRokUqKirSoEGDlJmZ6byo9fjx4/Ly+jbjREVFaevWrZo3b54GDhyo7t27a86cOZo/f37zfQsAAHDNatQFrLNnz9bs2bNrnZeTk1OjLTExUX/+858bsykAANDG8W4aAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGNSqMrFixQtHR0QoICFBCQoLy8vLcWu7tt9+WzWbT3Xff3ZjNAgCANsjjMLJx40alpKQoLS1Ne/bsUWxsrJKTk1VSUlLvcseOHdPDDz+s4cOHN3qwAACg7fE4jCxbtkwzZ87U9OnTFRMTo1WrVikoKEhr1qypc5mqqipNmjRJ6enpuvHGG5s0YAAA0Lb4eNK5srJSu3fvVmpqqrPNy8tLSUlJys3NrXO5J554QmFhYZoxY4Z27NjR4HYqKipUUVHhnC4tLZUk2e122e12T4Zcr+p1Nec62ypq5Rnq5T5q5T5q5T5q5b6WrJW76/QojHz99deqqqpSeHi4S3t4eLj+/ve/17rMxx9/rNdee035+flubycjI0Pp6ek12j/88EMFBQV5MmS3ZGVlNfs62ypq5Rnq5T5q5T5q5T5q5b6WqFVZWZlb/TwKI546d+6cJk+erNWrV6tLly5uL5eamqqUlBTndGlpqaKiojRq1CiFhIQ02/jsdruysrI0cuRI+fr6Ntt62yJq5Rnq5T5q5T5q5T5q5b6WrFX1mY2GeBRGunTpIm9vbxUXF7u0FxcXKyIiokb/I0eO6NixYxo7dqyzzeFwfLNhHx8dPHhQvXv3rrGcv7+//P39a7T7+vq2yE7VUutti6iVZ6iX+6iV+6iV+6iV+1qiVu6uz6MLWP38/BQXF6fs7Gxnm8PhUHZ2thITE2v079evnz7//HPl5+c7P3fddZd+8IMfKD8/X1FRUZ5sHgAAtEEen6ZJSUnR1KlTNXToUMXHx2v58uW6cOGCpk+fLkmaMmWKunfvroyMDAUEBOiWW25xWT40NFSSarQDAIDrk8dhZPz48Tp16pQWLVqkoqIiDRo0SJmZmc6LWo8fPy4vLx7sCgAA3NOoC1hnz56t2bNn1zovJyen3mVff/31xmwSAAC0URzCAAAARhFGAACAUYQRAABgFGEEAAAYRRgBAABGEUYAAIBRhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARjXqrb0AWl+Vw1JewWmVnCtXWPsAxffqJG8vm+lh4TrGPonmQhgBrgGZ+wqV/t5+FZ4td7ZFdghQ2tgYjb4l0uDIcL1in0Rz4jQNcJXL3FeoB9fvcflHX5KKzpbrwfV7lLmv0NDIcL1in0RzI4wAV7Eqh6X09/bLqmVedVv6e/tV5aitB9D82CfREggjwFUsr+B0jf99Xs6SVHi2XHkFp1tvULiusU+iJRBGgKtYybm6/9FvTD+gqdgn0RIII8BVLKx9QLP2A5qKfRItgTACXMXie3VSZIcA1XWzpE3f3MEQ36tTaw4L1zH2SbQEwghwFfP2siltbIwk1fjHv3o6bWwMz3ZAq2GfREsgjABXudG3RGrlvUMU0cH1sHdEhwCtvHcIz3RAq2OfRHPjoWfANWD0LZEaGRPB0y5x1WCfRHMijADXCG8vmxJ7dzY9DMCJfRLNhdM0AADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwqlFhZMWKFYqOjlZAQIASEhKUl5dXZ9/Vq1dr+PDh6tixozp27KikpKR6+wMAgOuLx2Fk48aNSklJUVpamvbs2aPY2FglJyerpKSk1v45OTmaMGGCPvroI+Xm5ioqKkqjRo3SyZMnmzx4AABw7fM4jCxbtkwzZ87U9OnTFRMTo1WrVikoKEhr1qyptf+bb76phx56SIMGDVK/fv306quvyuFwKDs7u8mDBwAA1z4fTzpXVlZq9+7dSk1NdbZ5eXkpKSlJubm5bq2jrKxMdrtdnTp1qrNPRUWFKioqnNOlpaWSJLvdLrvd7smQ61W9ruZcZ1tFrTxDvdxHrdxHrdxHrdzXkrVyd502y7Isd1f6j3/8Q927d9fOnTuVmJjobH/00Ue1bds2ffrppw2u46GHHtLWrVv1t7/9TQEBAbX2Wbx4sdLT02u0b9iwQUFBQe4OFwAAGFRWVqaJEyfq7NmzCgkJqbOfR0dGmuqpp57S22+/rZycnDqDiCSlpqYqJSXFOV1aWuq81qS+L+Mpu92urKwsjRw5Ur6+vs223raIWnmGermPWrmPWrmPWrmvJWtVfWajIR6FkS5dusjb21vFxcUu7cXFxYqIiKh32eeee05PPfWU/vjHP2rgwIH19vX395e/v3+Ndl9f3xbZqVpqvW0RtfIM9XIftXIftXIftXJfS9TK3fV5dAGrn5+f4uLiXC4+rb4Y9fLTNld65plntGTJEmVmZmro0KGebBIAALRxHp+mSUlJ0dSpUzV06FDFx8dr+fLlunDhgqZPny5JmjJlirp3766MjAxJ0tNPP61FixZpw4YNio6OVlFRkSQpODhYwcHBzfhVAADAtcjjMDJ+/HidOnVKixYtUlFRkQYNGqTMzEyFh4dLko4fPy4vr28PuKxcuVKVlZX6yU9+4rKetLQ0LV68uGmjBwAA17xGXcA6e/ZszZ49u9Z5OTk5LtPHjh1rzCYAAMB1olXvpgFw7alyWMorOK2Sc+UKax+g+F6d5O1lc3u+CVfjmJqq8pJD63OPqbOkdbnHdO+w3vLz4fViaBsIIwDqlLmvUOnv7Vfh2XJnW2SHAKWNjdHoWyIbnG/C1TimpsrYsl+rdxTI18vSM/HS01sP6pcfHNLM4b2UemeM6eEBTUasBlCrzH2FenD9Hpdf6pJUdLZcD67fo4wt++udn7mvsDWHK6nhMZsYU1NlbNmv/95eIMcVj6d0WNJ/by9Qxpb9ZgYGNCPCCIAaqhyW0t/br9oez2z932f1joI650tS+nv7VXXlb9AW1NCYTYypqSovObR6R0G9fVbvKFDlJUcrjQhoGYQRADXkFZyucXThSvX9TrckFZ4tV17B6eYdWD0aGrOJMTXVutxj9dZZ+ubnsC73WKuMB2gphBEANZScqz+ItPZ6mnNbrTmmpvrydFmz9gOuVoQRADWEta/73VEm1tOc22rNMTVVz07uvRjU3X7A1YowAqCG+F6dFNkhQPXdDOtlU53zbfrmDpb4Xp1aYHS1a2jMJsbUVJMTo9XQHcletm/6AdcywgiAGry9bEob+80to1f+LrT932fm8F51zpektLExrfpsj4bGbGJMTeXn4+Wsc11mDu/F80ZwzWMPBlCr0bdEauW9QxTRwfW0RkSHAK28d4hS74ypd76JZ3o0NOZr8TkjqXfG6D9v7VXjCImXTfrPW3nOCNoGHnoGoE6jb4nUyJiIOp9m2tD8q3HM16LUO2P0X6P6af3OI9K/9mt+cl+ewIo2hTACoF7eXjYl9u7c6PkmXI1jaio/Hy9NTozWli37NTkxWr4EEbQh7M0AAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADCKMAIAAIwijAAAAKMIIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMAowggAADDKx/QAAOBaUuWwlFdwWiXnyhXWPkDxvTrJ28smSbpYWaUnt+zXsX+WKbpzkB67M0aBft5uLVvfPEmqvOTQ+txj6ixpXe4x3Tust/x83Pv/ZFO229j1Vo95Xe4xfXm6TD07BWlyYnSzjLmlx43W16gwsmLFCj377LMqKipSbGysfv3rXys+Pr7O/r/5zW+0cOFCHTt2TH369NHTTz+tO++8s9GDBgATMvcVKv29/So8W+5si+wQoLSxMfrfPV8pa3+Js33HF9K6Px/XyJgwrZ7y3XqXlVTnvNG3RCpjy36t3lEgXy9Lz8RLT289qF9+cEgzh/dS6p0xjR5zQ9tt7HovH7PD+naZpVsONHnMo2+JbHB+U9YNMzwOIxs3blRKSopWrVqlhIQELV++XMnJyTp48KDCwsJq9N+5c6cmTJigjIwM/du//Zs2bNigu+++W3v27NEtt9zSLF8CAFpa5r5CPbh+j6wr2ovOluuB9XvqXC5rf4nuemmHPv+q1KNli86W68H1e5QUE+YScqo5LOm/txdIUp2/3Bsz5urtrrx3SJ2/nOtbb0uO+cH1e/Qft/bSK9sL6pzflHHXtyxalsfXjCxbtkwzZ87U9OnTFRMTo1WrVikoKEhr1qyptf8LL7yg0aNH65FHHlH//v21ZMkSDRkyRC+99FKTBw8AraHKYSn9vf01folJqrXtSn+tJYg0tKz1f5/afqlfbvWOAlVectRob+yYq+elv7dfVY6aPRtab0uO2fq/Zev7To0dd33LouV5dGSksrJSu3fvVmpqqrPNy8tLSUlJys3NrXWZ3NxcpaSkuLQlJydr06ZNdW6noqJCFRUVzumzZ89Kkk6fPi273e7JkOtlt9tVVlamf/7zn/L19W229bZF1Moz1Mt910Ktdn/5L5365z+NX2Tn47BUVuaQj91LVY5vr3FYnbVXP4vv4dK3qWM+9c8Lys4/orieHZt1vdWaMub6/hddPe6B3YJd9it31l3Xd27rWvLv4Llz5yRJltVAyLM8cPLkSUuStXPnTpf2Rx55xIqPj691GV9fX2vDhg0ubStWrLDCwsLq3E5aWlp1CObDhw8fPnz4XOOfEydO1JsvTAf9WqWmprocTXE4HDp9+rQ6d+4sm635rnguLS1VVFSUTpw4oZCQkGZbb1tErTxDvdxHrdxHrdxHrdzXkrWyLEvnzp1Tt27d6u3nURjp0qWLvL29VVxc7NJeXFysiIiIWpeJiIjwqL8k+fv7y9/f36UtNDTUk6F6JCQkhJ3VTdTKM9TLfdTKfdTKfdTKfS1Vqw4dOjTYx6MLWP38/BQXF6fs7Gxnm8PhUHZ2thITE2tdJjEx0aW/JGVlZdXZHwAAXF88Pk2TkpKiqVOnaujQoYqPj9fy5ct14cIFTZ8+XZI0ZcoUde/eXRkZGZKkOXPm6LbbbtPzzz+vMWPG6O2339auXbv0yiuvNO83AQAA1ySPw8j48eN16tQpLVq0SEVFRRo0aJAyMzMVHh4uSTp+/Li8vL494DJs2DBt2LBBjz/+uB577DH16dNHmzZtuiqeMeLv76+0tLQap4RQE7XyDPVyH7VyH7VyH7Vy39VQK5tlNXS/DQAAQMvhRXkAAMAowggAADCKMAIAAIwijAAAAKOuizCycuVKDRw40PlAl8TERH3wwQfO+eXl5Zo1a5Y6d+6s4OBg/fu//3uNB7Vdj5566inZbDbNnTvX2UatvrV48WLZbDaXT79+/ZzzqZWrkydP6t5771Xnzp0VGBioAQMGaNeuXc75lmVp0aJFioyMVGBgoJKSkvTFF18YHLEZ0dHRNfYrm82mWbNmSWK/ulxVVZUWLlyoXr16KTAwUL1799aSJUtc3oPCfvWtc+fOae7cuerZs6cCAwM1bNgwffbZZ875RmvV0Pto2oI//OEP1ubNm61Dhw5ZBw8etB577DHL19fX2rdvn2VZlvXAAw9YUVFRVnZ2trVr1y7r//2//2cNGzbM8KjNysvLs6Kjo62BAwdac+bMcbZTq2+lpaVZN998s1VYWOj8nDp1yjmfWn3r9OnTVs+ePa1p06ZZn376qXX06FFr69at1uHDh519nnrqKatDhw7Wpk2brL1791p33XWX1atXL+vixYsGR976SkpKXPaprKwsS5L10UcfWZbFfnW5pUuXWp07d7bef/99q6CgwPrNb35jBQcHWy+88IKzD/vVt+655x4rJibG2rZtm/XFF19YaWlpVkhIiPXVV19ZlmW2VtdFGKlNx44drVdffdU6c+aM5evra/3mN79xzjtw4IAlycrNzTU4QnPOnTtn9enTx8rKyrJuu+02ZxihVq7S0tKs2NjYWudRK1fz58+3vv/979c53+FwWBEREdazzz7rbDtz5ozl7+9vvfXWW60xxKvWnDlzrN69e1sOh4P96gpjxoyx7rvvPpe2H//4x9akSZMsy2K/ulxZWZnl7e1tvf/++y7tQ4YMsX7xi18Yr9V1cZrmclVVVXr77bd14cIFJSYmavfu3bLb7UpKSnL26devn3r06KHc3FyDIzVn1qxZGjNmjEtNJFGrWnzxxRfq1q2bbrzxRk2aNEnHjx+XRK2u9Ic//EFDhw7VT3/6U4WFhWnw4MFavXq1c35BQYGKiopc6tWhQwclJCRcl/WqVllZqfXr1+u+++6TzWZjv7rCsGHDlJ2drUOHDkmS9u7dq48//lh33HGHJPary126dElVVVUKCAhwaQ8MDNTHH39svFZX5Vt7W8Lnn3+uxMRElZeXKzg4WL/73e8UExOj/Px8+fn51XgRX3h4uIqKiswM1qC3335be/bscTmPWK2oqIhaXSYhIUGvv/66+vbtq8LCQqWnp2v48OHat28ftbrC0aNHtXLlSqWkpOixxx7TZ599pp///Ofy8/PT1KlTnTWpfpJzteu1XtU2bdqkM2fOaNq0aZL4O3ilBQsWqLS0VP369ZO3t7eqqqq0dOlSTZo0SZLYry7Tvn17JSYmasmSJerfv7/Cw8P11ltvKTc3VzfddJPxWl03YaRv377Kz8/X2bNn9dvf/lZTp07Vtm3bTA/rqnLixAnNmTNHWVlZNdIzaqr+35ckDRw4UAkJCerZs6feeecdBQYGGhzZ1cfhcGjo0KF68sknJUmDBw/Wvn37tGrVKk2dOtXw6K5er732mu64444GX79+vXrnnXf05ptvasOGDbr55puVn5+vuXPnqlu3buxXtVi3bp3uu+8+de/eXd7e3hoyZIgmTJig3bt3mx7a9XE3jfTNG4dvuukmxcXFKSMjQ7GxsXrhhRcUERGhyspKnTlzxqV/cXGxIiIizAzWkN27d6ukpERDhgyRj4+PfHx8tG3bNr344ovy8fFReHg4tapHaGiovvOd7+jw4cPsV1eIjIxUTEyMS1v//v2dp7Wqa3LlXSHXa70k6csvv9Qf//hH3X///c429itXjzzyiBYsWKCf/exnGjBggCZPnqx58+Y5X9TKfuWqd+/e2rZtm86fP68TJ04oLy9PdrtdN954o/FaXTdh5EoOh0MVFRWKi4uTr6+vsrOznfMOHjyo48ePKzEx0eAIW9+IESP0+eefKz8/3/kZOnSoJk2a5Pwztarb+fPndeTIEUVGRrJfXeF73/ueDh486NJ26NAh9ezZU5LUq1cvRUREuNSrtLRUn3766XVZL0lau3atwsLCNGbMGGcb+5WrsrIylxezSpK3t7ccDock9qu6tGvXTpGRkfrXv/6lrVu3aty4ceZr1eKXyF4FFixYYG3bts0qKCiw/vrXv1oLFiywbDab9eGHH1qW9c2tcj169LD+9Kc/Wbt27bISExOtxMREw6O+Olx+N41lUavL/dd//ZeVk5NjFRQUWJ988omVlJRkdenSxSopKbEsi1pdLi8vz/Lx8bGWLl1qffHFF9abb75pBQUFWevXr3f2eeqpp6zQ0FDr97//vfXXv/7VGjdu3HV7C2ZVVZXVo0cPa/78+TXmsV99a+rUqVb37t2dt/a+++67VpcuXaxHH33U2Yf96luZmZnWBx98YB09etT68MMPrdjYWCshIcGqrKy0LMtsra6LMHLfffdZPXv2tPz8/KyuXbtaI0aMcAYRy7KsixcvWg899JDVsWNHKygoyPrRj35kFRYWGhzx1ePKMEKtvjV+/HgrMjLS8vPzs7p3726NHz/e5bkZ1MrVe++9Z91yyy2Wv7+/1a9fP+uVV15xme9wOKyFCxda4eHhlr+/vzVixAjr4MGDhkZr1tatWy1JtX5/9qtvlZaWWnPmzLF69OhhBQQEWDfeeKP1i1/8wqqoqHD2Yb/61saNG60bb7zR8vPzsyIiIqxZs2ZZZ86ccc43WSubZV32qDoAAIBWdt1eMwIAAK4OhBEAAGAUYQQAABhFGAEAAEYRRgAAgFGEEQAAYBRhBAAAGEUYAQAARhFGgDbIZrPV+1m8eLHpITa76OhoLV++3PQwADSCj+kBAGh+hYWFzj9v3LhRixYtcnlRXXBwsIlhecyyLFVVVcnHp/X+qaqsrJSfn1+rbQ8AR0aANikiIsL56dChg2w2m0vb22+/rf79+ysgIED9+vXTyy+/7Fz22LFjstlseueddzR8+HAFBgbqu9/9rg4dOqTPPvtMQ4cOVXBwsO644w6dOnXKudy0adN09913Kz09XV27dlVISIgeeOABVVZWOvs4HA5lZGSoV69eCgwMVGxsrH7729865+fk5Mhms+mDDz5QXFyc/P399fHHH+vIkSMaN26cwsPDFRwcrO9+97v64x//6Fzu9ttv15dffql58+Y5j/5I0uLFizVo0CCX2ixfvlzR0dE1xr106VJ169ZNffv2lSSdOHFC99xzj0JDQ9WpUyeNGzdOx44da44fD4ArEEaA68ybb76pRYsWaenSpTpw4ICefPJJLVy4UG+88YZLv7S0ND3++OPas2ePfHx8NHHiRD366KN64YUXtGPHDh0+fFiLFi1yWSY7O1sHDhxQTk6O3nrrLb377rtKT093zs/IyND//M//aNWqVfrb3/6mefPm6d5779W2bdtc1rNgwQI99dRTOnDggAYOHKjz58/rzjvvVHZ2tv7yl79o9OjRGjt2rI4fPy5Jevfdd3XDDTfoiSeeUGFhocuRIXdkZ2fr4MGDysrK0vvvvy+73a7k5GS1b99eO3bs0CeffKLg4GCNHj3aJVwBaCat8jo+AMasXbvW6tChg3O6d+/e1oYNG1z6LFmyxPka+oKCAkuS9eqrrzrnv/XWW5YkKzs729mWkZFh9e3b1zk9depUq1OnTtaFCxecbStXrrSCg4Otqqoqq7y83AoKCrJ27tzpsu0ZM2ZYEyZMsCzLsj766CNLkrVp06YGv9fNN99s/frXv3ZO9+zZ0/rVr37l0ictLc2KjY11afvVr35l9ezZ02Xc4eHhLm96XbdundW3b1/L4XA42yoqKqzAwEBr69atDY4NgGe4ZgS4jly4cEFHjhzRjBkzNHPmTGf7pUuX1KFDB5e+AwcOdP45PDxckjRgwACXtpKSEpdlYmNjFRQU5JxOTEzU+fPndeLECZ0/f15lZWUaOXKkyzKVlZUaPHiwS9vQoUNdps+fP6/Fixdr8+bNKiws1KVLl3Tx4kXnkZGmGjBggMt1Inv37tXhw4fVvn17l37l5eU6cuRIs2wTwLcII8B15Pz585Kk1atXKyEhwWWet7e3y7Svr6/zz9XXYFzZ5nA4PN725s2b1b17d5d5/v7+LtPt2rVzmX744YeVlZWl5557TjfddJMCAwP1k5/8pMFTJl5eXrIsy6XNbrfX6Hfl9s6fP6+4uDi9+eabNfp27dq13m0C8BxhBLiOhIeHq1u3bjp69KgmTZrU7Ovfu3evLl68qMDAQEnSn//8ZwUHBysqKkqdOnWSv7+/jh8/rttuu82j9X7yySeaNm2afvSjH0n6JixceTGpn5+fqqqqXNq6du2qoqIiWZblDFT5+fkNbm/IkCHauHGjwsLCFBIS4tFYAXiOC1iB60x6eroyMjL04osv6tChQ/r888+1du1aLVu2rMnrrqys1IwZM7R//35t2bJFaWlpmj17try8vNS+fXs9/PDDmjdvnt544w0dOXJEe/bs0a9//esaF89eqU+fPnr33XeVn5+vvXv3auLEiTWOykRHR2v79u06efKkvv76a0nf3GVz6tQpPfPMMzpy5IhWrFihDz74oMHvMWnSJHXp0kXjxo3Tjh07VFBQoJycHP385z/XV1991fgCAagVYQS4ztx///169dVXtXbtWg0YMEC33XabXn/9dfXq1avJ6x4xYoT69OmjW2+9VePHj9ddd93l8oC1JUuWaOHChcrIyFD//v01evRobd68ucFtL1u2TB07dtSwYcM0duxYJScna8iQIS59nnjiCR07dky9e/d2nkrp37+/Xn75Za1YsUKxsbHKy8vTww8/3OD3CAoK0vbt29WjRw/9+Mc/Vv/+/TVjxgyVl5dzpARoATbryhOqANAI06ZN05kzZ7Rp0ybTQwFwjeHICAAAMIowAgAAjOI0DQAAMIojIwAAwCjCCAAAMIowAgAAjCKMAAAAowgjAADAKMIIAAAwijACAACMIowAAACjCCMAAMCo/w84A34Aq5eZXgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "data_pred = pd.DataFrame({'Temperature': np.linspace(start=30, stop=90, num=121), 'Intercept': 1})\n",
    "data_pred['Frequency'] = logmodel.predict(data_pred)\n",
    "data_pred.plot(x=\"Temperature\",y=\"Frequency\",kind=\"line\",ylim=[0,1])\n",
    "plt.scatter(x=data[\"Temperature\"],y=data[\"Frequency\"])\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hideCode": false,
    "hidePrompt": false,
    "scrolled": true
   },
   "source": [
    "This figure is very similar to the Figure 4 of Dalal *et al.* **I have managed to replicate the Figure 4 of the Dalal *et al.* article.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing and plotting uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the documentation of [Seaborn](https://seaborn.pydata.org/generated/seaborn.regplot.html), I use regplot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkoAAAG/CAYAAACnlPiuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABWL0lEQVR4nO3deXhU1cEG8PfeOzOZTJKZZCAQdgiQsG8qEIKIuCKIAi5UbalihTbVivarqF+pfGirdLEK7kuLVrFacEWjuIJsakGQRZYEwhIIIdtMklnvPd8fkwwMyYWQuSST5P09T4S565nDJHk959xzJCGEABERERHVITd3AYiIiIhiFYMSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIR0wFpYKCAsyfPx/XXHMNBgwYgMmTJzfoPCEEnn/+eYwfPx5DhgzBjTfeiO+///7cFpaIiIhavZgKSnv27MFXX32FHj16oHfv3g0+74UXXsCTTz6Jn//853juueeQmpqK2267DQcPHjyHpSUiIqLWToqltd40TYMsh7LbvHnzsG3bNnzwwQenPcfn82HMmDG4+eabcc899wAA/H4/rrzySowbNw4PPfTQuS42ERERtVIx1aJUG5LOxqZNm1BZWYmJEyeGt1ksFlx22WVYvXq1kcUjIiKiNiamglJj5OfnAwDS09Mjtvfu3RuFhYXwer3NUSwiIiJqBVp8UHK5XLBYLIiLi4vYbrfbIYRARUVFM5WMiIiIWroWH5TOpRgavkVERETNwNTcBYiW3W6H3++Hz+eLaFVyuVyQJAkOh6PR15YkCS6XB6qqGVHUNktRZNjt8azLKLEejcO6NA7r0hisR+M4HPGNGvOsp8UHpdqxSfv27UO/fv3C2/Pz89G5c2dYrdaorq+qGoJBfmiNwLo0BuvROKxL47AujcF6jJ7RnUEtvuttxIgRSExMxEcffRTeFggE8Mknn2DcuHHNWDIiIiJq6WKqRcnj8eCrr74CABw+fBiVlZXIzc0FAIwcORJOpxMzZ85EYWEhVq1aBQCIi4vD7NmzsXjxYjidTmRkZGDZsmUoLy/HrFmzmu29EBERUcsXU0GppKQEv/nNbyK21b5+5ZVXMGrUKGiaBlVVI475xS9+ASEEXn75ZZSWlqJ///546aWX0K1btyYrOxEREbU+MTUzdywqK6tif3GUTCYZKSkJrMsosR6Nw7o0DuvSGKxH4zidCVAU40YWtfgxSkRERETnCoMSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHTEXFDKy8vDrbfeimHDhiE7OxuLFi2C3+8/43llZWWYP38+xo8fj2HDhmHy5MlYtmxZE5SYiIiIWitTcxfgZBUVFZg5cyZ69uyJxYsXo6ioCI8++ii8Xi/mz59/2nN/85vfID8/H/fccw86deqE1atX46GHHoKiKLjhhhua6B0QERFRaxJTQemNN95AVVUVlixZguTkZACAqqpYsGABZs+ejY4dO9Z7XnFxMTZu3Ig//elPmDZtGgAgKysLP/zwA1auXMmgRERERI0SU11vq1evRlZWVjgkAcDEiROhaRrWrl2re14wGAQAJCUlRWxPTEyEEOKclJWIiIhav5gKSvn5+UhPT4/YZrfbkZqaivz8fN3zOnXqhLFjx+LZZ5/F3r17UVlZiQ8//BBr167FzTfffK6LTURERK1UTHW9uVwu2O32OtsdDgcqKipOe+7ixYsxd+5cTJo0CQCgKAr+93//F1dccUVUZVKUmMqSLVJtHbIuo8N6NA7r0jisS2OwHo0jScZeL6aCUmMJIXD//fdj//79+Otf/4rU1FSsW7cOf/zjH+FwOMLhqTHs9ngDS9q2sS6NwXo0DuvSOKxLY7AeY09MBSW73Q63211ne0VFBRwOh+55X375JXJzc/Hee+8hMzMTADBq1CiUlJTg0UcfjSoouVweqKrW6PMp9H9Idns86zJKrEfjsC6Nw7o0BuvROA5HPGTZuJa5mApK6enpdcYiud1uFBcX1xm7dLK9e/dCURRkZGREbO/fvz/eeusteDwexMc3LqWrqoZgkB9aI7AujcF6NA7r0jisS2OwHqNn9DNcMdUZOm7cOKxbtw4ulyu8LTc3F7IsIzs7W/e8Ll26QFVV7Nq1K2L79u3b0a5du0aHJCIiImrbYioozZgxAwkJCcjJycHXX3+N5cuXY9GiRZgxY0bEHEozZ87EZZddFn49btw4dO7cGXfddRfeffddrF+/Hn/+85/x9ttv45ZbbmmOt0JEREStQEx1vTkcDixduhQLFy5ETk4OEhIScN1112Hu3LkRx2maBlVVw68TExPxz3/+E48//jj+8pe/wO12o2vXrpg3bx6DEhERETWaJDgj42mVlVWxvzhKJpOMlJQE1mWUWI/GYV0ah3VpDNajcZzOBEOnWYiprjciIiKiWMKgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiIdMReU8vLycOutt2LYsGHIzs7GokWL4Pf7G3RuUVER7rvvPowePRpDhgzBxIkT8d57753jEhMREVFrZWruApysoqICM2fORM+ePbF48WIUFRXh0Ucfhdfrxfz580977rFjx3DjjTeiV69eWLhwIRITE7Fnz54GhywiIiKiU8VUUHrjjTdQVVWFJUuWIDk5GQCgqioWLFiA2bNno2PHjrrn/vnPf0ZaWhpefPFFKIoCAMjKymqKYhMREVErFVXX27Fjx4wqBwBg9erVyMrKCockAJg4cSI0TcPatWt1z6usrMRHH32Em266KRySiIiIiKIVVYvS+PHjMXr0aEyZMgWXX345bDZbVIXJz8/H9OnTI7bZ7XakpqYiPz9f97zt27cjEAjAZDLhlltuwebNm5GcnIxrr70Wd999N8xmc6PLpCgxN4yrxamtQ9ZldFiPxmFdGod1aQzWo3EkydjrRRWU7rrrLnzwwQeYN28eFixYgEsuuQRTpkzB2LFjIctn/4/tcrlgt9vrbHc4HKioqNA97/jx4wCA//3f/8UNN9yAX//619i6dSuefPJJyLKMe++996zLUstuj2/0uRSJdWkM1qNxWJfGYV0ag/UYe6IKSnPmzMGcOXOwY8cOvP/++1i5ciU++OADtGvXDpMmTcLVV1+NwYMHG1VWXZqmAQDGjBmDefPmAQBGjx6NqqoqvPzyy8jJyYHVam3UtV0uD1RVM6ysbZGiyLDb41mXUWI9God1aRzWpTFYj8ZxOOIb1Vijx5DB3AMGDMCAAQPwu9/9Dhs2bMD777+PFStW4NVXX0WvXr0wZcoUTJkyBZ07dz7tdex2O9xud53tFRUVcDgcpz0PCIWjk2VlZeHZZ59FQUEBMjMzG/HOAFXVEAzyQ2sE1qUxWI/GYV0ah3VpDNZj9IQw9nqGdoZKkoTzzjsPF110EYYOHQohBAoKCrBkyRJceumluOuuu047ADw9Pb3OWCS3243i4mKkp6frntenT5/Tlsvn853dGyEiIiKCgUFpw4YNePDBB5GdnY27774bx48fx3333YevvvoKa9aswb333osNGzbgd7/7ne41xo0bh3Xr1sHlcoW35ebmQpZlZGdn657XpUsXZGRkYN26dRHb161bB6vVesYgRURERFSfqLrefvzxR7z33ntYuXIljh07hvbt2+O6667DtddeW6era9asWYiLi8Njjz2me70ZM2bg1VdfRU5ODmbPno2ioiIsWrQIM2bMiJhDaebMmSgsLMSqVavC2+bOnYtf/epXeOSRRzB+/Hj88MMPePnllzFr1qyon8YjIiKitimqoHTttdfCarXikksuwbXXXovs7OzTDqDq06cPhg0bprvf4XBg6dKlWLhwIXJycpCQkIDrrrsOc+fOjThO0zSoqhqxbcKECfjb3/6Gp59+GsuWLUOHDh1w55134o477ojmLRIREVEbJgnR+GFPK1aswBVXXIGEhAQjyxRTysqqOLAuSiaTjJSUBNZllFiPxmFdGod1aQzWo3GczgRD56OKqkVp2rRpRpWDiIiIKOZEFbleeeUVzJo1S3f/7bffjtdffz2aWxARERE1m6iC0n/+8x/07t1bd3+fPn3w5ptvRnMLIiIiomYTVVA6ePDgaYNSeno6Dhw4EM0tiIiIiJpNVEHJbDajuLhYd/+xY8cMnUaciIiIqClFlWKGDh2Kt99+G5WVlXX2ud1urFixAkOHDo3mFkRERETNJqqn3n7961/jlltuwbXXXouZM2eGZ8Des2cPli5diuLiYvz1r381pKBERERETS2qoDR06FA8++yzmD9/Ph555BFIkgQAEEKga9eueOaZZzB8+HBDCkpERETU1KIKSgCQnZ2NVatWYceOHeGB2927d8fAgQPDwYmIiIioJYo6KAGALMsYNGgQBg0aZMTliIiIiGKCIUFp7969OHjwICoqKurdf+211xpxGyIiIqImFVVQOnDgAP7nf/4HW7duhd6ScZIkMSgRERFRixRVUJo/fz52796NBx54AOeffz7sdrtR5SIiIiJqdlEFpU2bNmH27Nn46U9/alR5iIiIiGJGVBNOpqSkICkpyaiyEBEREcWUqILSjBkz8N5770FVVaPKQ0RERBQzoup669mzJzRNwzXXXIPp06cjLS0NiqLUOe7yyy+P5jZEREREzSKqoDR37tzw3x977LF6j5EkCTt37ozmNkRERETNIqqg9MorrxhVDiIiIqKYE1VQGjlypFHlICIiIoo5hszM7ff7sX37dpSUlGDEiBFwOp1GXJaIiIioWUX11BsQ6n4bO3YsbrrpJtx5553YtWsXAKC0tBSjRo3Cf/7zn6gLSURERNQcogpKy5cvxx//+EdceOGFeOSRRyKWMXE6nRg9ejQ+/PDDqAtJRERE1ByiCkr/+Mc/cMkll+Cvf/0rLr744jr7Bw4ciD179kRzCyIiIqJmE1VQKigowLhx43T3Jycno7y8PJpbEBERETWbqIKS3W5HWVmZ7v69e/ciNTU1mlsQERERNZuogtK4cePw5ptvwuVy1dm3Z88evPXWW5gwYUI0tyAiIiJqNlFND3D33XfjhhtuwOTJk3HxxRdDkiS88847WL58OT755BOkpqbiV7/6lVFlJSIiImpSUbUodezYEStWrMCFF16Ijz76CEIIvPvuu/jiiy8wadIkvPnmm5xTiYiIiFosSZz8TH+USktLoWkanE4nZDnqKZpiQllZFYJBrbmL0aKZTDJSUhJYl1FiPRqHdWkc1qUxWI/GcToToCjGZRBDZuauxdYjIiIiak2iCkpLliw54zGSJCEnJyea2xARERE1i3MWlCRJghCCQYmIiIharKiC0o8//lhnm6ZpOHz4MF5//XV8++23eOGFF6K5BREREVGzMXzEtSzL6NatG+677z706NEDDz/8sNG3ICIiImoS5/TRtAsuuABfffXVubwFERER0TlzToPStm3bWs00AURERNT2RDVG6Z133ql3u8vlwnfffYdPPvkE119/fTS3ICIiImo2UQWlefPm6e5LSUnBHXfcwSfeiIiIqMWKKih99tlndbZJkgS73Y7ExMRoLk1ERETU7KIKSl26dDGqHEREREQxhyOtiYiIiHRE1aLUr18/SJJ0VudIkoQdO3ZEc1siIiKiJhFVUMrJycGnn36KvXv3YuzYsejVqxcAID8/H2vXrkXfvn1x6aWXGlJQIiIioqYWVVDq0KEDSkpK8P777yM9PT1iX15eHmbOnIkOHTrghhtuiKqQRERERM0hqjFKL730Em655ZY6IQkAevfujZtvvhkvvvhiNLcgIiIiajZRBaWjR4/CZNJvlDKZTDh69Gg0tyAiIiJqNlEFpb59++L1119HUVFRnX1Hjx7FsmXLkJGREc0tiIiIiJpNVGOU7r//ftx+++244oorcOmll6JHjx4AgP379+Ozzz6DEAKLFi0ypKBERERETS2qoHT++efjzTffxBNPPIFPP/0UXq8XAGC1WjF27FjceeedyMzMNKSgRERERE0tqqAEABkZGXjqqaegaRpKS0sBAE6nE7LMuSyJiIioZYs6KNWSZRlxcXGw2WwMSURERNQqRJ1ofvjhB8yaNQtDhw7FqFGj8M033wAASktL8ctf/hIbN26MupBEREREzSGqoLRp0ybcdNNNKCgowJQpU6BpWnif0+lEZWUl/v3vf0ddSCIiIqLmEFVQevzxx9G7d298+OGHmDt3bp39o0aNwpYtW6K5BREREVGziSoo/fDDD5g2bRosFku9i+N27NgRx48fj+YWRERERM0mqqBkMpkiuttOVVRUBJvNFs0tiIiIiJpNVEFp6NCh+Pjjj+vdV11djRUrVuCCCy6I5hZEREREzSaqoHTXXXdh27ZtuOOOO7B69WoAwK5du/DWW29h2rRpKC0txa9+9StDCkpERETU1CQhhIjmAuvXr8dDDz2EgoKCiO3du3fHww8/jJEjR0ZVwOZWVlaFYFC/e5HOzGSSkZKSwLqMEuvROKxL47AujcF6NI7TmQBFMW4+x0ZPOCmEQFVVFUaMGIGPP/4YO3fuxP79+yGEQLdu3TBo0KB6B3gTERERtRSNjlyBQAAjR47EK6+8AgDo378/Jk6ciKuuugqDBw9udEjKy8vDrbfeimHDhiE7OxuLFi2C3+8/q2v885//RGZmJmbPnt2oMhAREREBUbQoWSwWtG/fHhaLxbDCVFRUYObMmejZsycWL16MoqIiPProo/B6vZg/f36DrlFcXIynnnoK7dq1M6xcRERE1DZFtdbb1KlT8e677+InP/mJIYHpjTfeQFVVFZYsWYLk5GQAgKqqWLBgAWbPno2OHTue8Rp//vOfMWHCBBQWFkZdHiIiImrbogpKmZmZ+OyzzzB58mRMnToVXbp0gdVqrXPc5Zdf3qDrrV69GllZWeGQBAATJ07EH/7wB6xduxbTpk077fnfffcdPv30U+Tm5uLee+89q/dCRCGSJCHKZzyIiFqNqILSPffcE/77E088Ue8xkiRh586dDbpefn4+pk+fHrHNbrcjNTUV+fn5pz1XVVUsXLgQc+bMQYcOHRp0v4YwcuR8W1Vbh6zL6DRVPWpCwCTLaM1ZiZ9J47AujcF6NI7Rz5GddVD629/+hquuugr9+vULD+Q2isvlgt1ur7Pd4XCgoqLitOe+/vrr8Hg8+PnPf25omez2eEOv15axLo1xruvR4wvApMgwm5Rzep9YwM+kcViXxmA9xp6zDkrPP/88+vbti379+mHkyJEoKyvDmDFj8PLLLyMrK+tclPGMSkpK8OSTT+Kxxx4zdHA5ALhcHqgq57SIhqLIsNvjWZdRaqp69Ac1BFUNtrioGpxjGj+TxmFdGoP1aByHIx6yHAPzKJ3MqPEMdrsdbre7zvaKigo4HA7d85544glkZmbi/PPPh8vlAgAEg0EEg0G4XC7YbDaYTI17q6qqcfIvg7AujXGu61FVNVR7Aogzte7uN4CfSSOxLo3Beoye0T+3Yup/GdPT0+uMRXK73SguLkZ6erruefv27cO3335b77pyF1xwAV544QWMGzfO8PIStVaqJhAIajBxvAQRtXExFZTGjRuHZ599NmKsUm5uLmRZRnZ2tu55DzzwQLglqdYf//hHWK1W3HPPPcjMzDyn5SZqbVRNwOtXkWRr/a1KRESn06igdPjwYWzfvh0Awl1lBQUF9Q7EBoCBAwc26LozZszAq6++ipycHMyePRtFRUVYtGgRZsyYETGH0syZM1FYWIhVq1YBCM0Kfiq73Q6bzYZRo0ad1XsjIgAC8PlVJMSbIIFLERFR29WooPTEE0/UmQ5gwYIFdY4TQpzV9AAOhwNLly7FwoULkZOTg4SEBFx33XWYO3duxHGapkFV1cYUnYgaKKhp8AU0WM2t/+k3IiI9kjjLkdhvv/32Wd9k6tSpZ31OrOBKztHjqtjGaKp6DKgaylw+aELAalGQkmRtdRNQ8jNpHNalMViPxnE6Ewydj+qsW5RacughorMTqJkqQJHZ/UZEbRMfaSEiXaom4PEFDZ/ploiopWBQIqLT8vlVaK2s642IqKEYlIjotIKaBp+fD08QUdvEoEREpyUEUO0LNncxiIiaBYMSEZ1RIKjB42dYIqK2h0GJiM5ICKDay6BERG0PgxIRNUhQDbUq8Qk4ImpLGJSIqEFqW5Va2+STRESnw6BERA0WalVS2apERG0GgxIRNRhblYiorWFQIqKzwlYlImpLznqtNyJqnTQhcKDIDY9PhS8QRIcUG+R60lBtq1K8hT8+KPbVfq4rqwNItJnRvWNSvZ9rIj38SUdE2Lm/FCs3FKCotBrtHPGoqPQi0WbB+KGd0btrcp3jg6oGrz8IK8MSxbDaz/XR0mqoqoCiSEhz2jBpdA/07+ls7uJRC8GuN6I2buf+Uiz9eBcOFVcizqwgId4Mi1nB0VIP3v56H/IOldc5h7N1U6w79XNtT7QgzqzgUHEVln68Czv3lzZ3EamFYFAiasM0IbByQwG8/iCSE+NgMSuQJcBkUmC3meELaPhyS2G9i+IGghq8Aa4BR7Gn/s+1BItZQXKiBV6/ipUbCrjYMzUIgxJRG3agyI2jpdVIsJohnTJuQ5Ik2OIUHC/34MjxqjrnCgF4vAGAwz0oxpzpc51gNeFoaTUOFLmbqYTUkjAoEbVhldUBqKqAyVT/jwJFkaFq+suXBIICPrYqUYw50+faZJKhqgKV1YEmLhm1RAxKRG1Yos0MRZEQDGr17ldVDYoM2Kz1D9rWhEC1l8uaUGw50+c6GNSgKBISbeYmLhm1RAxKRG1Y945JSHPaUFXPJJJCCFT7VLRPjken9gm61wgENfh1fiERNYczfa6rvEGkOW3o3jGpmUpILQmDElEbJksSJo3uAatFQXmlH/6ACk0AwaAKV3UAcWYZ44d2Pu28M5pW26rEZiWKDfV/rgX8ARXllX5YLQomje7B+ZSoQRiUiNq4/j2dmHlFJrqmJsAXUFHlCcAfUJHmjMfUsb3qnUfpVP6AiqDKViWKHad+rl2VfvgCKrqmJmDmFZmcR4kajLPFERH693Qis0dKg2bmro+qCVT7grDbLFwHjmLGyZ9rzsxNjcWgREQAQt0VPdPsCKgayly+s55jxusPwhZngiLzlxDFjtrPNVFjseuNiAyhqgJV3gDHKhFRq8KgRESG8flVBIKcV4mIWg8GJSIyjKqFHr1mqxIRtRYMSkRkKF9AhZ+tSkTUSjAoEZGhNE2g0hPgbN1E1CowKBGR4QIBDR4/W5WIqOVjUCIiw4XWgOOCo0TU8jEoEdE5EQhqqPZxwVwiatkYlIjonBACqPYGwZVNiKglY1AionMmqGqchJKIWjQGJSI6p7z+IIJBNisRUcvEoERE55SqClR6/ByrREQtEoMSEZ1zvqAGX4DTBRBRy8OgRETnnKYJVHkCAFuViKiFYVAioibhD2rwcLoAImphGJSIqEkIAVR7gtCEaO6iEBE1GIMSETWZgKqh0hPkdAFE1GIwKBFRk/L6gvAHObCbiFoGU3MXgIhigyYEPtpQgFK3D1kD0pBoM5+T+6iaQKUngJSkOIC9cEQU4xiUiAgAsGXPcSz/Kj/891mTB4TCzDngD6jweIOwWU3gkCUiimXseiMiAIAj8UQoKnX78ML721Hq8p6TewkBVHmDUDWmJCKKbQxKRAQASO9sx4xL+oZfl1f68cL7O1BScW7CUlDVUFnNdeCIKLYxKBFR2OUXdMONl/QJv66o8uOF97fjWLnnnNzPG1DhCwTPybWJiIzAoEREESaM6Iprx/YKv3ZVB/Di+ztwtLTa8HtpNQO7BUd1E1GMYlAiojpGD0zD1HHp4RVHKj0BvPD+DhwurjT8Xv4Au+CIKHYxKBFRvS7o1wHTx/cOLzni8QXx4gc7UXDUbfi9PH4VXnbBEVEMYlAiIl0jMlJx44S+kGvSki+g4uUPd2LPoXJD76NpApXVAfAhOCKKNQxKRHRaQ3q3wy2XZ8CkhMJSIKjhldxd2Lav1ND7BIIa3FU+gD1wRBRDGJSI6Iz69UjBz67sB4sp9CND1QSWfbob3/14zND7eAMqxysRUUxhUCKiBunTxYFZk/sjPk4BEJo0csXqfKzZUmjYPYQAqn1BePxBMCsRUSxgUCKiBuvWIQm/uHogkk5aB+6jjQfw0YYCCIPWItE0AXe1H/6gZsj1iIiiwaBERGclzWnD7CkD4bSfWPJkzdYjWP5VHlTNmHCjqgKuKj+XOCGiZsegRERnzWm3YvaUgejUzhbetmn3cbz68W74Aqoh9wgENbiqfJyMkoiaFYMSETVKks2CX1w9AL062cPbdh8sx4sf7EClJ2DIPXwBDRWVfoBhiYiaCYMSETWa1WLCzyf2w8BezvC2w8VVePbdbThu0PpwPr+KiqoApw0gombBoEREUTGbZPzkkr4YPbBjeFupy4dn391uyCzeAoDXH6yZNiDqyxERnRUGJSKKmixLuHpMT1w5snt4W7UviJdW7sAP+SVRX18IoNobRJU3yDmWiKhJxVxQysvLw6233ophw4YhOzsbixYtgt/vP+05x44dw6JFi3DNNddg+PDhGDduHO69914cPny4iUpNRJIkYdywzrhxQh8ocijMBFWBZZ/uwVffH456+gBNCFR5ApxjiYialKm5C3CyiooKzJw5Ez179sTixYtRVFSERx99FF6vF/Pnz9c9b/v27Vi1ahWmT5+OoUOHoqysDM888wyuv/56fPDBB3A6nbrnEpGxhvZpjySbBa+t2gWPL/QE3MffHERJhRdTxvaCSWn8/5+pmoC7yg/AgniLAoOmbiIi0hVTQemNN95AVVUVlixZguTkZACAqqpYsGABZs+ejY4dO9Z73nnnnYePPvoIJtOJtzNixAiMHz8e77zzDm677bamKD4R1UjvbMecawZh6Uc/otTtAwB8t6sYJS4fbr4sAzZr43/0hMOSMCM+zmzYRJdERPWJqa631atXIysrKxySAGDixInQNA1r167VPc9ut0eEJABIS0uD0+nEsWPGrkVFRA2TmhyPOdcOQo+OSeFt+4648Mw723AsyifiVE3AVR1ApcfPbjgiOqdiqkUpPz8f06dPj9hmt9uRmpqK/Pz8s7rWvn37UFJSgt69e0dVJiWKbgIKqa1D1mV0mqoeNQCKIkEyoKHGkWjBHdcMwH++zMPm3ccBACUuL559ZxtuuqwvMrunRHV9jz8IALAnWiCdxfwB/Ewah3VpDNajcYz+n6eYCkoulwt2u73OdofDgYqKigZfRwiBhx9+GB06dMCkSZOiKpPdHh/V+XQC69IY57oePb4AVEkydPzPHVOH4KP1+/He6tD/8Hj9Kv7x4Y+YNr4PLh3ZPeon2VTISE6Kg9mknNV5/Ewah3VpDNZj7ImpoGSUxYsXY8OGDXjxxRdhs9nOfMJpuFweqCoX54yGosiw2+NZl1Fqqnr0BzW4XF5oBo/9yR7YEY54E974bC8CQQ1CAMu/2Iv8w+WYflH6WYecU1W4PHAkxsEknzl08TNpHNalMViPxnE44iHLxrXMxVRQstvtcLvrTlBXUVEBh8PRoGu8+eabeOqpp/DII48gKysr6jKpqoYgVzE3BOvSGOe6HlVVg6oKw4MSAPTv4cScawbi1Y93obwyNO3H5t3HUVTqwc2XZSAlKe4MV9DnUYMIBjQkJZgRZ27YE3H8TBqHdWkM1mP0jP7RFVOdoenp6XXGIrndbhQXFyM9Pf2M569atQoPPfQQ7rrrLlx33XXnqphEFIVO7RLwq6mD0avTiUHehcer8NTbP2Dv4YZ3sdcnoIbWhqv2cq4lIjJGTAWlcePGYd26dXC5XOFtubm5kGUZ2dnZpz1348aNuOeee3D99dcjJyfnXBeViKKQGG/GbZP6Ryx7Uu0N4h8f7ox6ckpVE3BXB+CqNmZhXiJq22IqKM2YMQMJCQnIycnB119/jeXLl2PRokWYMWNGxBxKM2fOxGWXXRZ+nZeXh5ycHPTs2RPXXHMNvv/++/DXgQMHmuOtENEZKLKMKdm9cN343jApoeYfIUKTU762aje8NU+0NYYmBKq9AVRU+qBxmiUiikJMjVFyOBxYunQpFi5ciJycHCQkJOC6667D3LlzI47TNA2qqoZfb9myBW63G263Gz/5yU8ijp06dSoeffTRJik/EZ29ERmp6Oi04bVPToxb2rG/DEtW/ICbLs1A5/YJjbquEIDHryKoeZFoM8NqNnFySiI6a5LgT47TKiur4sC6KJlMMlJSEliXUWqqegyoGspcvnMymPt0qr0B/Pvzvdhz6MQ4JZMiYfKYnrigX4eophCQZQlWi4IEqxkmRYaiSPxMGoTf38ZgPRrH6UwwdD6qmGpRIqK2y2Y146dXZOKDdfvxzc7QjPpBVeCdNfuQX+jCtRf2gsWs4MjxKlR7g7BZTejUPgFyAwKUpglUe4Pw+VXEWRQ4Ei3n9L1oQuBAkRuV1QEk2szo3jGpQeVsaYKaho0/FKHap8EWJ+O8zFSYDHwsmygWMCgRUUzIO1SOL7cU4ni5B1aLAp9fRW2b1ta8Euw/4oI9wQJ3tR+qBigy0D45HuOHdkbvrskNuodaE5iCmgZznAUCxrea7dxfipUbCnC0tBqqKqAoEtKcNkwa3QP9e7aeBbpzNxZg5foCeHxBCAASgH/FmTApqweuHNWjuYtHZBhGfyJqdnmHyvH21/twtLQaFrOC5KQ4pNjjIJ80eaSrOoBDxVUIqgIJ8SZYzAqOlnrw9tf7kHeo/Kzup6oClZ4ASit88AZUw6YS2Lm/FEs/3oVDxZWIMyuwJ1oQZ1ZwqLgKSz/ehZ37S425UTPL3ViA5V/lo8obhCxLMCkSZFlClTeI5V/lI3djQXMXkcgwDEpE1Kw0IfDllkL4AirsNgvMJgWyJMFqMaFDshVmJTLFVHmDKHf7Icsy7DYzfAENX24pbNSYKn9QhavSjzK3D0FVi2oclCYEVm4ogNcfRHJiHCzm0PuwmBUkJ1rg9atYuaGgycd+GS2oaVi5vgCqJmCuCUiyJEOWJZgVCaomsHJ9AYIax9lQ68CgRETN6sjxKhwv98AWZ6oTVGRZRmK8GUDkQpe+gIrismr4AipscQqOl3tw5HhVo+6vCQGvX0WZ2wdXtR9BTTQqMB0ocuNoaTUSrOY650uShASrCUdLq3GgqO7qAy3JNzuK4PEFYZKlet+nSZbg8QXxzY6iZiohkbEYlIioWVV7g6ExR3pPqdT8LnbYzIiznFgPThNAqcuHSk8QQTU09igaqiZQ5QmgzO1FRZUPqnp2gamyOgBVFTCZ6n8fJpMc6vJr4RNhlrq8oZFdelUjAaLmOKLWgEGJiJqVzWqCIkN/IVAR+p0sSRKcSXFwJEQ+sebxBeHxBeGq8htSHrUmdJW6vXBVhaZJaEheSrSZoSiS7qPdwaAGRZGQaDMbUs7m4rRbQxlJrwex5t/Labc2XaGIziEGJSJqVp3aJ6B9cjyqfWqdCSGFEPAHNVgtCvw1QSoh3ozU5HiYT2qB0gSwYk0+Pv7mAIIGrbyuagJV3iBKXT54fGdeO657xySkOW2o8gbrfR9V3iDSnDZ075ikc4WWYeSAjoiPMyGoiXrfZ1ATiI8zYeSAjjpXIGpZGJSIqFnJkoTxQzsjzizDVR1AIKhCEwKBoApXdQBxFgXjh3dGnFkJ71cUCY5EM8wndXMJAXz1fSGeWvEDDhdXGla+oKrBVRVA+RkGfMuShEmje8BqUVBe6Yc/EHof/oCK8ko/rBYFk0b3aPHzKZlkGZOyekCRJQRUAU0T0IQGTRMIqAKKLGFSVg/Op0StBj/JRNTsendNxtSxvZDmjIc/oKKyOgB/QEWaMx5Tx/bChUO71NkfCGro3jERk7O6o91J3TxFZR488842fPLtQcNalzQh4PGrKHX7UF6pH5j693Ri5hWZ6JqaAF8g9ESdL6Cia2oCZl6R2WrmUbpyVA9MvygdCVYTNE0gWBOYEqwmTL8onfMoUavCJUzOgNPJR49T8xujtS9hAoQCyelm3tbb7w+oyP3mADZsj3zSKjU5HtMvSq/T3aUoEhwOGyoqQpNCni1ZlhBnVpBoNcNkkut0QbWlmbn/u6uYM3MbgD8njWP0EiYMSmfAD230+APAGG0hKEUrv7ACy7/KR5nbF94mARg9MA2XX9At/NRctEGpllK7hly8GYosoQVWWdT4/W0M1qNxjA5KjP5E1Gqkd3bgN9cNwZhBaeGn1wWA9duP4u9vbcHOgjJD73fygO8qb1D/kXkiarEYlIioVbGYFUwe0xN3TBmI1OT48PaKKj9e/XgX/vXJLpRX+k5zhbMXVDW4q0MzfPuDmmFLohBR82NQIqJWqUdaEu6cPhgTRnSBctKacTv2l+Gvy77HJxsL9OduagQhAJ9fRXmlD6UuL7wBlS1MRK0AgxIRtVomRcal53fDndOHoGfaiQHd/qCGFV/sxRNvbUV+YYWh99Q0AV9AC60h5/LCF1QBgK1MRC0UgxIRtXodUuJx+9UDMG1cOmxxpvD2ojIPXvxgJ974bI/h3XGaCAWmCrcfJS4vKj0BaIKBiailMZ35ECKilk+WJJzfrwMG9EzBJ98exLc7j4VX4diaV4KdBWW4aFhnXDikc8REltHShIAWFAgENXh9KhLiTYiPM7XJJ+SIWiK2KBFRm2KzmjF9fG/87mfno2tqQnh7IKjh0+8O4fE3v8eWvcfrzI1khEDNLN9lbh+8ARUCDVtHjoiaD4MSEbVJvTo7kDN9MKaNC80wXau80o9/f74Xz767HQVH3YbfVxMCXr+KikofSip8cFX54Q9qgMRuOaJYxK43ImqzarvjBqU78fl/D2P99qNQtVBL0sFjlXjuve0Y2MuJKy7ohvYnTTVgBCFC0woEVQ0eX2j9ujizgjiLArNJhlRzDBE1LwYlImrzrBYTrsrqgZEDOiB34wHs2H9iYsrt+0qxc38pzu/XARPO6wq7zWL4/U8ex1TtC0KRT4QmS814KYYmoubBoEREVKO9Ix63XJ6J/EIXPtpYgMPFVQAATQDf7DyGzbuPY8zgNIwb2hnxcefmx6emhRaYrQ1NppqWJotZgVmRIUkMTURNiUGJiOgU6Z3t+OW1g/BDXglWfXsQpTVrxwVUDV99X4iNO4owdkgnjBmUBqvl3P0Y1TQBvybgD2iQ5SAUSUKchS1NRE2JQYmIqB6yJGFon/YY2MuJb3YewxebD6PKEwAAeP0qPv3uENb+cBQXDumErIFp4QV3zxVNE9AgEPCc1D1nURBn5pgmonOJQYmI6DRMiowxg9JwXmYq1v1wFGu2FsLrD8227fEF8cm3B7Fm6xGMHdwJWYM6ntMWploR3XPyiTFNFnOopUmSJACCwYnIAAxKREQNEGdWcPGILhg9sCO+3noEa7cdgT8QWivO4wti1XcHsWZrIbIGpWHMoDQkWM1NUq6TQ5PsDUKWJZgVCWaTDMUkwyTJkBVAluRzMjcUUWvHoEREdBbi40y47IJuyB6chq+3HsG67UfDgcnrV/HFpsP4eusRjOzXAdlDOiE5Ma7JyqYJAU0VCKqAx69CkgBJkiBLEiwmGXFxoRYnhiaihmNQIiJqBJvVjMtHdsfYIZ2w9oejWLftKHyBUJdcIKhh7bajWL+9CMP6tsPYIZ2R5rQ1eRmFAIQIjW0Kqho8/iAUWQ6HpriabjpmJiJ9DEpERFGwWc247IJuGDukEzZsL8LabUdQ7Q0CCLXwbNp9HJt2H0ffrg6MHdIJfbo4asYQNb2ISS5rQlPtBJcWkwJF5tgmolMxKBERGSA+zoSLR3RB9uA0fLfrGL7eegTllf7w/j2HKrDnUAU6pMQje1AahvZtD4vp3D4pdzrh0OTRIEmALEswKTLiTDJMNaFJkSXIcm2LEwMUtU0MSkREBrKYFYwZ1AmjBnTED3mlWLO1EEdKqsP7j5V58Paafcj95gAu6NcBowZ0REqStRlLHApNqiqgqip8fhWyFAyPb5IkwCRLUEwyzCYZJlmGLIemTwAkrk9HrR6DEhHROaDIMob1bY+hfdohr9CFtVuPYNfB8vB+j0/F6i1HsGbLEWR2T8aoAR3Rt1tyTQBpXpoQQE0rEgAEAOCkweGSBEiQoCgSLGYFkskUGp8lAFkOvXe2QlFrwaBERHQOSZKEPl0c6NPFgeJyD9ZtO4rNu4vhD4aelBMAfjxQjh8PlCMlKQ4X9OuAEZmp52RNuWjVDg6veYWgGuq+M3sCcLm9EFpkK5SsyOHuO0WWIEkSlJNao2qDGMMUxTJJ8BnR0yorq0Kw5gcaNY7JJCMlJYF1GaWmqseAqqHM5Qu1KrRSiiLB4bChoqIaqtr079PrD+K/u4qxYUcRSiq8dfbLEtCvRwrOz+yAvt2SawZZx6aG1GVtCxSk0N9lKRScFEWGUjM+SpblmjAVei2h9j23jVYp/pw0jtOZAEWRDbseW5SIiJqY1WJC9uBOyBqUhrzDFdiwvQg/HigLBwJNADv2l2HH/jIk2cwY3jcVIzJS0SElvnkL3kihnrza7jxAhQh15yE0nUJtb+PJ3XqyDChSKEzJJ7VKybXHnMhR4cYpqWZ/3WDVNsIWnRsMSkREzUSWJPTtmoy+XZNRUenDd7uK8d2Px1BRdeJpOXd1AKu3FGL1lkJ065CIYX3bY0jvdk0283dTqA0xJ3fr1WQonBymwq1SQM2fUiiA1ajdL9eOpao5z6TIodarmpY5WZLqjKU6uRxEJ2NQIiKKAY7EOFxyXldcPLwL9hwqx393FWNnQRlU7cRv74PHKnHwWCVWritARjcHhvZpj/49UmAxN980A03l1Fapmq2nHqVztnoiaAHhLkBFlqBIodYqqebvoVYp1IynCu2LbKViqGprGJSIiGKILEvI7J6CzO4pqPIG8P2e49i8uxiFJ00xoAkRHgBuNsno3yMFQ3q3Q9+uyTCbjBub0ZqEgxZwoguwnjFVEd2AONEdWF94qu3qk+UT3YUSpPDTgafev/bmDFktC4MSEVGMSrCakT24E7IHd8KRkip8v+c4vt97HO7qQPiYQFDD1rwSbM0rQZxZQf8eKRiU7mRoaqR6uwF1nBqqcHKrVfg1aroDJcgAFFOoC9B00hOBteOuaq8p1ZwkSeA0CzGAQYmIqAXo1C4Bndol4IqR3ZFf6MKWvcexbV9peH05APAFVHy/NxSmLCYZGd2SMaCXE/26J8Nq4Y97o9UNVcDpghWA0HxUOGngeu10CiYZAcioqvRDDWpQhYCmiROhSgm1ZsWZFYamJsbvHCKiFkSWJfTp6kCfrg5MGdsLew6VY2teCX4sKAvPzQQA/qCGbftKsW1fKRRZQq9OdvTvkYJ+PVKQkhTXjO+ABGrC1UkDyTUh4A+oqPYGIroET3460GSS0c4uR7Za0TnHoERE1EKZTTIG9HRiQE8n/EEVuw9WYFt+CX48UAZ/4ERoUjWBvYcrsPdwBd5ftx9pThsyuycjs3syunVIiul5mugkbElqFgxKREStgMWkYFAvJwb1ciIQ1JB3uALb95diZ0EZqr3BiGOPllbjaGk1vvq+EFaLgj5dHcjomoy+XR1wJLK1iehkDEpERK2M2SSjX003m6YJHDjmxs79ZdhZUIbjp8wE7vWr2JZfim35pQCA1OR49OnqQN8uDvTqZEecpfVPPUB0OgxKREStmCxL6JlmR880OyaO7oHjFR78WFCO3QfLse+IK2KeJgAoLveguNyD9duOQpaArh0Skd7JjvTODnRPS4TFxOBEbQuDEhHVdWK90pgnhf9T+1o68bpm7hw+JXRCe0c8xg6Jx9ghneALqMgvdGH3wXLsOVSOUpcv4lhNAAeKKnGgqBJffl8IRZbQNTURvToloWcnO7p3TERCfOuZIZyoPgxKRBRBkSQkWk2oncSldkiw0ER49XiBUydIPvGq3rmST90oQfe5nfD2mollapehCE/kJ9dMBIgTj1bjpDXCauehCZUV0DQNAVUgGFRD7yW8JliotSX8iLcmWko2NEztvEv9e6QAAEpcXuw9FBr0nV9YAY9PjThe1QQKitwoKHID3xdCkoBO7Wzo292Jzs54dE1NRHKipc5ki0QtGYMSEUWQZQm2mnXE6v991/Bfgmf6fXl2LT2RBzf8XBnxEgBYwuWRFQmJiVaYJECteaTeF1Dh8wcRUEPz15zs1PfRWluo2tmtaDfAilEDOkLTBI6UVCGv0IX8QhcKjroj5mwCQvVQeLwahcdPzBqeZDOje4ckdOuYiG4dEtGlfUKbWGKFWi8GJSLSVX8gaHhKiJVAcersxrIswWxSYJJr1qYAkGA1IcFqQlANzWcT1ER42Qq5piVLhCbAQVAV8PpVqJqGOsuPtRKyLKFLaiK6pCZi3NDOUDWBI8ersO+IC/uOuFBQ5K7T4gSEFvHdvr8U2/eHBodLEtAxxYYuqQmhr/aJSHPaOGs4tRgMSkREOBHqFFlCfJwp3IpUX9iTJCAx3gR/UIOqCQgtFJ4CqgatZkblWAmJRlFkCV07JKJrh0RcOLQzNCFwrMyDg8fcKCzxYO/BMpScMsYJCNVf7XQE/91VDCC0nEdHZzw6tUtA5/a2mlnHbZw9nGISP5VERPU4XdAJ7ZPCT4CdvOZXUNUQVDX4gxr8ARWqKqCJ1hecZElCmjPUUuRw2FBRUY2KSj8OFrlx4FglDh6rxKHiyoiJL2tpQuBISTWOlFRj0+4T21OS4pDmtKFTOxs6OkNf7exWTohJzYpBiYgoSiev+SVLoQAVZ1aAeHO4pSkY1OAP1gSnVjpwPDHejP49nejf0wkA0DSB4goPDhdX4VBxJQ4XV+FISRWCav3vvsztQ5nbh50FZeFtJkVCanI8OqTEo0OyDR1S4pGaEo929jgoMrvv6NxjUCIiOgdO7spTZAWSRQEQCk7+oAqvr/6B462JLEvomGJDxxQbRmSkAgg9OVdc7kHh8arQV0kVjhyvrjNQvFZQPdH6BJSEtyuyBKfditRkK1KT49HeYUV7RzzaOaxIsJr45B0ZhkGJiKgJRIyBsphgizMhGBTwBoI1A8Nbd2iqpcihLrs054nwJIRAmdsXHst0tCT0Z4nLq9tlWRu4iss9AMoi9lktCto7rHDaraEn+RxWOO1xcCZZkWgzQ2aIorPAoERE1AyEABRFQqLJjIR4MwJBDT6/Cl+g7YSmWpIUah1y2q0YUNNtBwCBoIbicg+KyqpxrMwT+ir3oPQ0AQoILctyqLgKh4qr6uwzKRJSkuKQkmSt+TP0lZwYh+RECxLjzWyNoggMSkREzaj2F75ZkWGxyUiEOTQY3K/CFwwNDBdaaAB0W2M2yejcPgGd2ydEbA+qGo5XeMMtSiU1fz9e4YXXX38X3olzBYrLvSgu99a736RIcCTEwZFoQXKiJfx3R4IF9oTQn6GnIhmm2goGJSKiGFGbhUyyDHO8jARI0ISGQDA0s3hbD061TIoc7r47mRAC1b4gjpd7UeIKfZW6vCh1+VDi8qLaGzzjtYOqCJ+rf38JdlsoOCXZzLDbLEiyWZBoMyPJZkaSLbQ9Ps7Ebr5WgEGJiCgG1U6SKUGCxSQhzlw3OPmDGoLh6Qda3xQEZ0uSJCRYzUhIM6NHWlKd/V5/MPxkXanLh7JKH8prXpe5fboDyk8VVAVK3T6UuuvOG3UyWZKQGG9CYnyoezWx5ish3owEa812qxk2qwn2RAtEW/8HjFEMSkRELYBecBJCQFUFAqoamoogqEKtaXFieIpktZjQqZ0Jndol1Lvf4wuiosqP8kofyit9qKj0h76q/Kio8sFV5ded2qA+mhBwVQfgqg406HiTIsNmVWCLC4UnW5wJNqsJ8TV/2m0WZA3qCLstrsFloOgxKBERtUC1wQkIDQpXlNrZxC3QhAZNA4KaBjWoIRDUENROTHzJAFW/+LhQKDm1S6+WEAIenwpXtR+uqpqvaj/c1QG4T/qz0hM4q0BVK6hqcFVpcFXpB6tV3x3En+4YDbOJ6+c1lZgLSnl5eXj44YexefNmJCQk4JprrsHdd98Ni8Vy2vOEEHjhhRfw+uuvo7S0FP3798f999+PYcOGNU3BiYiakCYEDhS5UVkdQKLNjO4dkyBDQm2rkypUrNp4AGVuLzo6bbj4/K5QJBmqJqBqAkFVQ2FxJdzVAcSZFaSmxEOSpHCI0kRobbdqbxA2qwmd2ifUGW+jahq25JXC69dgtcgY1Mt5VpNAnukeDSlDtPdQNQ1b9x5HeaUfyYkWDOnTXvc9SJIUaumxRoapU++R1s4Gf0BDpScQDk+Fx6tqWqQ0CABV3iCqPAFUeYPw+M48dqpWqcuHvMMVyOiewvFPTUQSMdQpWlFRgUmTJqFnz56YPXs2ioqK8Oijj2LKlCmYP3/+ac99/vnn8eSTT+K3v/0tMjMz8dprr2HdunV499130a1bt0aXqaysCsFg3Sn4qeFMJhkpKQmsyyixHo3T0uty5/5SrNxQgKOl1VBVAUUJzU00aXQP9O/pxNLcnViz9Qi0k96aLAMXDumEmVf2x48Fpfjku0MoqfBCkkJPl3V02nDRkM7o3ikJeYcrsO6HozhW4QktQSI0OBLjMGZAGnp2tgMA1mw5jC82FcLrVyEASAjNXzR+eGdcOLTLGd9D3qFyfLmlEMfLPVA1QJGB9snxGD+0M3p3TT7j/oY40zXWbDmMLzc3/j0Y9T5UTcAXCEI2m1BU7Ia7OgiPN4BqXxDV3iCOllaj8HgVVE0gKd4MVYiIf2+K5HQmQFGMm7U9poLSc889h2effRZffPEFkpOTAQD//ve/sWDBAnzxxRfo2LFjvef5fD6MGTMGN998M+655x4AgN/vx5VXXolx48bhoYceanSZWuoP0ljS0n8pxQrWo3Facl3u3F+KpR/vgtcfRILVDJNJRjCoocobhNWiIM0Zj615pbrnD+ntxNFSj+75E4Z3wdfbjiIQVGG3WWCxKBCagC+gIj7OhMmje+JwaSW++O9haCK0zp0QoVaqQDA0Lmr88M4YPSANtb2DtQu21P62yTtUjre/3gdfQIUtzgRFkaGqGqp9KuLMMkb174CNO4/p7p86ttcZw9KZ7tG7sx3/3VUMTQCKBNQ0xkEVgCwBV4zsdsawZOT7UBQpvGaeelK33cn3SHPaYLOacLzcG/73mnlFJsPSKYwOSjG1UM7q1auRlZUVDkkAMHHiRGiahrVr1+qet2nTJlRWVmLixInhbRaLBZdddhlWr159LotMRNRkNCGwckMBvP4gkhPjYDErobXlzAqSEy2o8vgjQpJ00letrXmlqPL46j3f6w/i/XX7UVHpRZxZgT+oobI61D0UCGo4UFSJ99bm4701+1Fc4UVltR/VvmBNQALMigSTImPT7mI4Ei1oZ7cixR6H5KQ4OBLj4EgIPUK/80A5kpPikN7Zjg5OG5x2K9onx6NragKsZgX/3X0cZlNoCoCkBAtscSYkxlvQ3hEHTQBfbzt62ifENCHw5ZZC+AKhsGc2hd6n2aTAbjPD6wviu5qQZJJDS63IkgRZlmCSAU0AX24uhKrph+gz3cPnV/Hl5tPsD2j4ckvhaad5OPUeJuXUfy8VKzcUtOmpIppCTI1Rys/Px/Tp0yO22e12pKamIj8//7TnAUB6enrE9t69e2Pp0qXwer2wWq2NKpPDEc9Bj1Gq7UZnXUaH9WicllqXQVXDvbecHwo/9YxPUTXRoEfMJUmCItc9v7ZlqDY41LdfqxkMDkQGsPAxNX8m2SywWuoOOA6qGn4xdXC4/FKoQJH3OF0Zao5JjDfXjCWq+36DqsDvZ40Ov9cTBQsdq4kT81CdfIeIKwnAZjXBYq5/0LSqanjw1lG6E0+G6iq0SHJ9x4hQUxsS4s0w1bR+yHJoAeAT70PDA7eOCv97S1Loz9pZ24UItdU5Eizha1CoHo0UU0HJ5XLBbrfX2e5wOFBRUXHa8ywWC+LiIh+ZtNvtEEKgoqKi0UFJ5urUhmFdGoP1aJyWVpe1T1LV/tKMRn3n12YsSWc/IAGiYV2VQoh6uz9C8z6FujMkqTa/nIgoWs1gcwC6fR6hoCDBZKr/AFVTIQQig5YU/k/4fqe+RemUF5IkwaxzDyFOtETVlxiFJgGaBkWWIdV3CSFBEwKKLEWEHFmJDI1KPYHxxEupJrRKhnY1UaSYCkpERKQvzqKgg6X+R9dbiqZ4D3EWBXGW+DMfGPP3MCHOwl/TzS2mIqjdbofb7a6zvaKiAg6H47Tn+f1++HyRs6S6XC5IknTac4mIiIj0xFRQSk9PrzMWye12o7i4uM74o1PPA4B9+/ZFbM/Pz0fnzp0b3e1GREREbVtMBaVx48Zh3bp1cLlc4W25ubmQZRnZ2dm6540YMQKJiYn46KOPwtsCgQA++eQTjBs37pyWmYiIiFqvmOr8nDFjBl599VXk5OSEJ5xctGgRZsyYETGH0syZM1FYWIhVq1YBAOLi4jB79mwsXrwYTqcTGRkZWLZsGcrLyzFr1qzmejtERETUwsVUUHI4HFi6dCkWLlyInJwcJCQk4LrrrsPcuXMjjtM0DaoaucrzL37xCwgh8PLLL4eXMHnppZeimpWbiIiI2raYmpmbiIiIKJbE1BglIiIioljCoERERESkg0GJiIiISAeDEhEREZEOBiUiIiIiHQxKRERERDrabFD66quvcMstt2D06NEYNGgQLrnkEvzpT3+qs9bc559/jilTpmDw4MG44oorsHz58mYqcctRVVWFcePGITMzEz/88EPEvrfeegtXXHEFBg8ejClTpuCLL75oplLGphUrViAzM7PO11/+8peI41iPDfP222/j2muvxeDBgzFq1Cjcfvvt8Hq94f38/j6zn/70p/V+JjMzM7Fy5crwcfxMNsxnn32G66+/HsOHD8fYsWPxm9/8BgcPHqxzHOvz9L744gtMnToVgwYNwkUXXYQnn3yyzvyKgEHf46KNeuedd8Rjjz0mcnNzxYYNG8Srr74qRo4cKW699dbwMd9++63o37+/+P3vfy/Wr18vHn/8cZGZmSk++uijZix57Fu0aJEYM2aMyMjIEFu3bg1v/+CDD0RmZqZ4/PHHxfr168Xvf/97MWDAALF58+bmK2yMWb58ucjIyBCrV68WmzdvDn8VFhaGj2E9NszTTz8thg8fLp577jmxceNGkZubK/7whz+IyspKIQS/vxtqz549EZ/FzZs3i7vvvlsMGDBAlJSUCCH4mWyoDRs2iH79+ol58+aJtWvXipUrV4rLL79cXHrppcLj8YSPY32e3ubNm0W/fv3EvffeK1avXi1efvllMWTIEPHoo49GHGfU93ibDUr1+fe//y0yMjLE0aNHhRBC3HbbbeLGG2+MOOaee+4REydObI7itQh79+4Vw4YNE8uWLasTlC6//HJxzz33RBx/4403ittvv72pixmzaoNS7S+g+rAezywvL08MGDBAfPnll7rH8Pu78SZMmCB+8YtfhF/zM9kwv//978WECROEpmnhbevXrxcZGRni22+/DW9jfZ7ebbfdJqZOnRqx7aWXXhIDBw4UxcXFEccZ8T3eZrve6pOcnAwgtKCu3+/Hxo0bceWVV0Ycc9VVVyEvLw+HDh1qhhLGvocffhgzZsxAr169IrYfPHgQ+/fvx8SJEyO2X3XVVVi/fj38fn9TFrPFYj02zIoVK9C1a1dcdNFF9e7n93fjbdq0CYcOHcLVV18NgJ/JsxEMBpGQkABJksLbkpKSAACiZpEM1ueZ7dy5E9nZ2RHbxo4di0AggK+//hqAsd/jbT4oqaoKn8+H7du346mnnsKECRPQtWtXHDhwAIFAAOnp6RHH9+7dGwCQn5/fHMWNabm5udi9ezdycnLq7Kutr1MDVO/evREIBOrto2/LJk+ejP79++OSSy7Bc889F+57Zz02zJYtW5CRkYGnn34aWVlZGDRoEGbMmIEtW7YAAL+/o/DBBx/AZrPhkksuAcDP5NmYNm0a8vLy8Nprr8HtduPgwYP429/+hgEDBmDEiBEAWJ8N4fP5YLFYIrbVvs7LywNg7Pd4TC2K2xwuvvhiFBUVAQAuvPBC/PWvfwUAVFRUAADsdnvE8bWva/dTiMfjwaOPPoq5c+ciMTGxzn7WZ8OkpqbizjvvxNChQyFJEj7//HP8/e9/R1FREebPn896bKDi4mJs27YNu3fvxh/+8AfEx8fj2WefxW233YZPPvmE9dhIwWAQH330ESZMmACbzQaA39tn4/zzz8eSJUtw77334v/+7/8AAP3798eLL74IRVEAsD4bokePHti6dWvEtu+//x7Aifoxsh7bfFB6/vnn4fF4sHfvXjzzzDOYM2cO/vGPfzR3sVqcZ555Bu3atcP06dObuygt2oUXXogLL7ww/Hrs2LGIi4vD0qVLMWfOnGYsWcsihEB1dTWeeOIJ9OvXDwAwdOhQTJgwAf/6178wduzYZi5hy7R27VqUlpZi8uTJzV2UFmnTpk343e9+hxtuuAHjx49HeXk5nn76adxxxx14/fXXYbVam7uILcJNN92EBx98EEuXLsU111yDvXv34u9//3s4bBqtzXe99evXD8OHD8f111+Pp59+Ghs3bsSqVavgcDgAoM50AS6XCwDC+wk4fPgwXn75Zdx1111wu91wuVyorq4GAFRXV6Oqqor1GYWJEydCVVXs3LmT9dhAdrsdycnJ4ZAEhMYgDhgwAHv37mU9NtIHH3yA5OTkiKDJumy4hx9+GKNHj8a8efMwevRoXHnllXj++eexY8cOvPvuuwBYnw0xbdo0zJw5E4sWLcKoUaPw85//HDNmzIDD4UCHDh0AGFuPbT4onSwzMxNmsxkHDhxA9+7dYTab6/Rj1r4+td+zLTt06BACgQDuuOMOXHDBBbjgggvCrR8/+9nPcOutt4brq776NJvN6NatW5OXuyViPTZMnz59dPf5fD5+fzeC1+vFp59+iiuvvBJmszm8nZ/JhsvLy4sI7wCQlpaGlJQUHDhwAADrsyFkWcYDDzyADRs24N1338W6detwww03oLS0FEOHDgUAQ7/HGZROsmXLFgQCAXTt2hUWiwWjRo3Cxx9/HHHMhx9+iN69e6Nr167NVMrY079/f7zyyisRX/fffz8AYMGCBfjDH/6Abt26oWfPnsjNzY0498MPP0RWVladgXl0wocffghFUTBgwADWYwNdfPHFKC8vx86dO8PbysrKsH37dgwcOJDf343w+eefo7q6Ovy0Wy1+Jhuuc+fO2LFjR8S2w4cPo6ysDF26dAHA+jwbSUlJ6NevH+x2O1599VV07doVY8aMAQBDv8fb7BilX//61xg0aBAyMzNhtVrx448/4qWXXkJmZiYuvfRSAMAvf/lL/OxnP8NDDz2EiRMnYuPGjfjggw/w+OOPN3PpY4vdbseoUaPq3Tdw4EAMHDgQAHDnnXfit7/9Lbp3745Ro0bhww8/xNatW/Gvf/2rKYsb02bNmoVRo0YhMzMTQGgW3zfffBM/+9nPkJqaCoD12BCXXnopBg8ejLvuugtz585FXFwcnn/+eVgsFtx0000A+P19tt5//3107twZ5513Xp19/Ew2zIwZM/DHP/4RDz/8MCZMmIDy8vLw+M6TpwNgfZ7e1q1b8c0336B///7wer34/PPP8e677+KFF16IGKdk2Pd4Yyd8aumee+45cc0114jhw4eLYcOGiUmTJom///3vwu12Rxz36aefismTJ4uBAweKyy67TLz11lvNVOKWZcOGDXUmnBRCiDfffFNcdtllYuDAgWLy5Mni888/b6YSxqaFCxeKyy+/XAwZMkQMGjRITJ48WSxdujRigjohWI8NUVJSIn7729+K8847TwwZMkTcdtttYs+ePRHH8Pu7YcrLy8XAgQPFokWLdI/hZ/LMNE0Tr7/+urj66qvFsGHDRHZ2tsjJyRF79+6tcyzrU9+OHTvE9ddfL4YNGyaGDRsmZs6cKTZt2lTvsUZ8j0tC1MxyRUREREQROEaJiIiISAeDEhEREZEOBiUiIiIiHQxKRERERDoYlIiIiIh0MCgRERER6WBQIiIiItLBoERERESko80uYUJEDVe7pMqZvPLKK7rL2bR2r732GuLj4zFt2rTmLgoRGYhBiYjOaNGiRRGv3333Xaxdu7bO9t69ezdlsWLKsmXLkJKSwqBE1MowKBHRGV1zzTURr7ds2YK1a9fW2d5aCCHg8/lgtVpZDqI2jmOUiMgQmqbhn//8JyZNmoTBgwdjzJgxmD9/PioqKiKOmzBhAmbPno2NGzdi2rRpGDJkCK6++mps3LgRAPDJJ5/g6quvxuDBgzFt2jTs2LEj4vx58+Zh+PDhOHjwIGbNmoVhw4Zh7NixWLJkCU5duvJsy7RmzZpwmd544w0AwPLly/Gzn/0MWVlZGDRoEK666iq8/vrrdc7fs2cPvvnmG2RmZiIzMxM//elPAQCLFy+ut+tyxYoVyMzMxKFDhxpUDpfLhUceeQQXXXQRBg0ahMsuuwzPP/88NE1r8L8REZ09tigRkSHmz5+Pt99+G9OmTcNPf/pTHDp0CK+99hp27NiBZcuWwWw2h48tKCjAvffeixkzZmDKlCl4+eWXMWfOHCxYsACPP/44fvKTnwAAnn/+edx9993Izc2FLJ/4/zpVVXH77bdj6NCh+J//+R+sWbMGixcvhqqq+M1vftOoMu3btw/33nsvbrzxRtxwww3o1asXgFCXWt++fTFhwgSYTCZ88cUXWLBgAYQQuPnmmwEADzzwABYuXAibzYY5c+YAANq3b9+oeqyvHB6PB7fccguKioowY8YMdOrUCZs3b8bf/vY3FBcX48EHH2zUvYioAQQR0VlasGCByMjICL/+9ttvRUZGhnjvvfcijlu9enWd7RdffLHIyMgQmzZtCm9bs2aNyMjIEEOGDBGHDx8Ob3/jjTdERkaG2LBhQ3jbfffdJzIyMsTChQvD2zRNE3fccYcYOHCgKCkpaXSZVq9eXee9ejyeOttuu+02cckll0RsmzRpkrjlllvqHPvkk09G1FWt5cuXi4yMDHHw4MEzluOpp54Sw4YNE/v27YvY/pe//EX0799fFBYW1rk+ERmDXW9EFLXc3FwkJSUhOzsbpaWl4a+BAwfCZrOFu9Vq9enTB8OHDw+/Hjp0KABg9OjR6Ny5c53tBw8erHPP2tYcAJAkCTfffDMCgQDWr1/fqDJ17doVF154YZ37nDw+yO12o7S0FCNHjsTBgwfhdrsbXEcNVV85cnNzcd5558Fut0e8lzFjxkBVVXz77beGl4OIQtj1RkRRKygogNvtRlZWVr37S0pKIl536tQp4nVSUhIAIC0tLWJ7YmIigND4nJPJsoxu3bpFbKvtKjt8+HCjytS1a9d6j/vvf/+LxYsX4/vvv4fH44nY53a7w2U3Sn3lKCgowK5du3TfS2lpqaFlIKITGJSIKGqapqFdu3b4y1/+Uu9+p9MZ8VpRlHqP09suThmkfS7KVN+TZQcOHMDPf/5zpKenY968eejUqRPMZjO++uor/POf/2zQQGpJkurdrqpqvdvrK4emacjOzsbtt99e7zk9e/Y8YzmIqHEYlIgoat27d8f69esxYsSIJnmUXdM0HDx4MNyKBIQGQQNAly5dDCvT559/Dr/fj2eeeSaiS/DUbjtAPxDZ7XYAoVax2r8DQGFhYYPL0b17d1RXV2PMmDENPoeIjMExSkQUtYkTJ0JVVTz99NN19gWDwTpdZ0Z47bXXwn8XQuC1116D2WwOd08ZUabaFq6TW7TcbjeWL19e59j4+Ph6r9m9e3cAiBhHVF1djXfeeeeM9681ceJEbN68GWvWrKmzz+VyIRgMNvhaRHR22KJERFEbOXIkbrzxRjz33HPYuXMnsrOzYTabsX//fuTm5uLBBx/ElVdeadj94uLisGbNGtx3330YMmQI1qxZgy+//BJz5swJd6kZUabac+bMmYMZM2agqqoKb731Ftq1a4fi4uKIYwcOHIhly5bh6aefRo8ePeB0OpGVlYXs7Gx07twZDz74IPLz86EoCpYvX46UlJQGtyrNmjULn3/+OebMmYOpU6di4MCB8Hg82L17Nz7++GN89tlndboSicgYDEpEZIj/+7//w6BBg/DGG2/g8ccfh6Io6NKlC6ZMmYIRI0YYei9FUfDiiy/ioYcewp///GckJCTg17/+NXJycgwtU3p6Op588kn8/e9/x2OPPYb27dvjJz/5CZxOJx544IGIY3NyclBYWIgXX3wRVVVVGDlyJLKysmA2m7FkyRIsWLAATzzxBFJTUzFz5kzY7Xbcf//9DXq/8fHxePXVV/Hcc88hNzcX77zzDhITE9GzZ0/ceeedhg8oJ6ITJNGYUZJERM1k3rx5+Pjjj7F58+bmLgoRtQEco0RERESkg0GJiIiISAeDEhEREZEOjlEiIiIi0sEWJSIiIiIdDEpEREREOhiUiIiIiHQwKBERERHpYFAiIiIi0sGgRERERKSDQYmIiIhIB4MSERERkQ4GJSIiIiId/w/DO7NULKqyUAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set(color_codes=True)\n",
    "plt.xlim(30,90)\n",
    "plt.ylim(0,1)\n",
    "sns.regplot(x='Temperature', y='Frequency', data=data, logistic=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**I think I have managed to correctly compute and plot the uncertainty of my prediction.** Although the shaded area seems very similar to [the one obtained by with R](https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/tree/master/challenger.pdf), I can spot a few differences (e.g., the blue point for temperature 63 is outside)... Could this be a numerical error ? Or a difference in the statistical method ? It is not clear which one is \"right\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Hide code",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}