From a65a047f6d925010046e8c06b500a46fd45fe94f Mon Sep 17 00:00:00 2001 From: e5e62e6e8091e12ec477af239ac0a4fc Date: Fri, 28 Jul 2023 14:30:46 +0000 Subject: [PATCH] exercice 4 --- module4/src_Python3_challenger.ipynb | 863 +++++++++++++++++++++++++++ 1 file changed, 863 insertions(+) create mode 100644 module4/src_Python3_challenger.ipynb diff --git a/module4/src_Python3_challenger.ipynb b/module4/src_Python3_challenger.ipynb new file mode 100644 index 0000000..bcb627b --- /dev/null +++ b/module4/src_Python3_challenger.ipynb @@ -0,0 +1,863 @@ +{ + "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": 8, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "3.9.17 (main, Jun 20 2023, 19:27:40) \n", + "[Clang 14.0.0 (clang-1400.0.29.202)]\n", + "uname_result(system='Darwin', node='simons-MacBook-Pro.local', release='21.0.1', version='Darwin Kernel Version 21.0.1: Tue Sep 14 20:56:24 PDT 2021; root:xnu-8019.30.61~4/RELEASE_ARM64_T6000', machine='arm64')\n", + "IPython 8.8.0\n", + "IPython.core.release 8.8.0\n", + "PIL 9.4.0\n", + "PIL.Image 9.4.0\n", + "PIL._deprecate 9.4.0\n", + "PIL._version 9.4.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", + "appnope 0.1.3\n", + "argparse 1.1\n", + "backcall 0.2.0\n", + "cffi 1.15.1\n", + "comm 0.1.2\n", + "csv 1.0\n", + "ctypes 1.1.0\n", + "ctypes.macholib 1.0\n", + "cycler 0.10.0\n", + "dateutil 2.8.2\n", + "debugpy 1.6.5\n", + "debugpy.public_api 1.6.5\n", + "decimal 1.70\n", + "decorator 5.1.1\n", + "defusedxml 0.7.1\n", + "distutils 3.9.17\n", + "entrypoints 0.4\n", + "executing 1.2.0\n", + "executing.version 1.2.0\n", + "http.server 0.6\n", + "ipykernel 6.19.4\n", + "ipykernel._version 6.19.4\n", + "jedi 0.18.2\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.8\n", + "jupyter_client._version 7.4.8\n", + "jupyter_core 5.1.2\n", + "jupyter_core.version 5.1.2\n", + "kiwisolver 1.4.4\n", + "kiwisolver._cext 1.4.4\n", + "logging 0.5.1.2\n", + "matplotlib 3.6.2\n", + "matplotlib._version 3.6.2\n", + "numpy 1.23.5\n", + "numpy.core 1.23.5\n", + "numpy.core._multiarray_umath 3.1\n", + "numpy.lib 1.23.5\n", + "numpy.linalg._umath_linalg 0.1.5\n", + "numpy.version 1.23.5\n", + "packaging 22.0\n", + "packaging.__about__ 22.0\n", + "pandas 1.5.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.appdirs 1.4.3\n", + "pkg_resources._vendor.more_itertools 8.12.0\n", + "pkg_resources._vendor.packaging 21.3\n", + "pkg_resources._vendor.packaging.__about__ 21.3\n", + "pkg_resources._vendor.pyparsing 3.0.9\n", + "pkg_resources._vendor.appdirs 1.4.3\n", + "pkg_resources._vendor.more_itertools 8.12.0\n", + "pkg_resources._vendor.packaging 21.3\n", + "pkg_resources._vendor.pyparsing 3.0.9\n", + "platform 1.0.8\n", + "platformdirs 2.6.2\n", + "platformdirs.version 2.6.2\n", + "prompt_toolkit 3.0.36\n", + "psutil 5.9.4\n", + "ptyprocess 0.7.0\n", + "pure_eval 0.2.2\n", + "pure_eval.version 0.2.2\n", + "pydevd 2.9.5\n", + "pygments 2.14.0\n", + "pyparsing 3.0.9\n", + "pytz 2022.7\n", + "re 2.2.1\n", + "scipy 1.10.0\n", + "scipy._lib._uarray 0.8.8.dev0+aa94c5a4.scipy\n", + "scipy._lib.decorator 4.0.5\n", + "scipy.integrate._dop 1.21.0\n", + "scipy.integrate._lsoda 1.21.0\n", + "scipy.integrate._vode 1.21.0\n", + "scipy.interpolate.dfitpack 1.21.0\n", + "scipy.linalg._fblas 1.21.0\n", + "scipy.linalg._flapack 1.21.0\n", + "scipy.linalg._flinalg 1.21.0\n", + "scipy.linalg._interpolative 1.21.0\n", + "scipy.optimize.__nnls 1.21.0\n", + "scipy.optimize._cobyla 1.21.0\n", + "scipy.optimize._lbfgsb 1.21.0\n", + "scipy.optimize._minpack2 1.21.0\n", + "scipy.optimize._slsqp 1.21.0\n", + "scipy.sparse.linalg._eigen.arpack._arpack 1.21.0\n", + "scipy.sparse.linalg._isolve._iterative 1.21.0\n", + "scipy.special._specfun 1.21.0\n", + "scipy.stats._mvn 1.21.0\n", + "scipy.stats._statlib 1.21.0\n", + "seaborn 0.12.2\n", + "seaborn.external.appdirs 1.4.4\n", + "seaborn.external.husl 2.1.0\n", + "setuptools 65.6.3\n", + "distutils 3.9.17\n", + "setuptools._vendor.more_itertools 8.8.0\n", + "setuptools._vendor.ordered_set 3.1\n", + "setuptools._vendor.packaging 21.3\n", + "setuptools._vendor.packaging.__about__ 21.3\n", + "setuptools._vendor.pyparsing 3.0.9\n", + "setuptools._vendor.more_itertools 8.8.0\n", + "setuptools._vendor.ordered_set 3.1\n", + "setuptools._vendor.packaging 21.3\n", + "setuptools._vendor.pyparsing 3.0.9\n", + "setuptools.version 65.6.3\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.8.0\n", + "traitlets._version 5.8.0\n", + "urllib.request 3.9\n", + "wcwidth 0.2.5\n", + "xmlrpc.client 3.9\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 os\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": 45, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "if os.path.exists('data_shuttle.csv'):\n", + " data = pd.read_csv('data_shuttle.csv')\n", + "else:\n", + " data = pd.read_csv(\"https://app-learninglab.inria.fr/moocrr/gitlab/moocrr-session3/moocrr-reproducibility-study/blob/master/data/shuttle.csv\")" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
DateCountTemperaturePressureMalfunction
04/12/81666500
111/12/81670501
23/22/82669500
311/11/82668500
44/04/83667500
\n", + "
" + ], + "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" + ] + }, + "execution_count": 46, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "data.head()" + ] + }, + { + "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": 47, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAxdElEQVR4nO3de1xVdb7/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\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "\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": "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": 70, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 1., 66.],\n", + " [ 1., 70.],\n", + " [ 1., 69.],\n", + " [ 1., 68.],\n", + " [ 1., 67.],\n", + " [ 1., 72.],\n", + " [ 1., 73.],\n", + " [ 1., 70.],\n", + " [ 1., 57.],\n", + " [ 1., 63.],\n", + " [ 1., 70.],\n", + " [ 1., 78.],\n", + " [ 1., 67.],\n", + " [ 1., 53.],\n", + " [ 1., 67.],\n", + " [ 1., 75.],\n", + " [ 1., 70.],\n", + " [ 1., 81.],\n", + " [ 1., 76.],\n", + " [ 1., 79.],\n", + " [ 1., 75.],\n", + " [ 1., 76.],\n", + " [ 1., 58.]])" + ] + }, + "execution_count": 70, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# We use add_constant to add intercept to our data\n", + "temps = data['Temperature'].values\n", + "temperature_with_intercept = sm.add_constant(temps, prepend=True)\n", + "temperature_with_intercept" + ] + }, + { + "cell_type": "code", + "execution_count": 77, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -3.9210
Date: Fri, 28 Jul 2023 Deviance: 3.0144
Time: 16:23:08 Pearson chi2: 5.00
No. Iterations: 6 Pseudo R-squ. (CS): 0.04355
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
const 5.0850 7.477 0.680 0.496 -9.570 19.740
x1 -0.1156 0.115 -1.004 0.316 -0.341 0.110
" + ], + "text/latex": [ + "\\begin{center}\n", + "\\begin{tabular}{lclc}\n", + "\\toprule\n", + "\\textbf{Dep. Variable:} & y & \\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: } & -3.9210 \\\\\n", + "\\textbf{Date:} & Fri, 28 Jul 2023 & \\textbf{ Deviance: } & 3.0144 \\\\\n", + "\\textbf{Time:} & 16:23:08 & \\textbf{ Pearson chi2: } & 5.00 \\\\\n", + "\\textbf{No. Iterations:} & 6 & \\textbf{ Pseudo R-squ. (CS):} & 0.04355 \\\\\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{const} & 5.0850 & 7.477 & 0.680 & 0.496 & -9.570 & 19.740 \\\\\n", + "\\textbf{x1} & -0.1156 & 0.115 & -1.004 & 0.316 & -0.341 & 0.110 \\\\\n", + "\\bottomrule\n", + "\\end{tabular}\n", + "%\\caption{Generalized Linear Model Regression Results}\n", + "\\end{center}" + ], + "text/plain": [ + "\n", + "\"\"\"\n", + " Generalized Linear Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: y 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: -3.9210\n", + "Date: Fri, 28 Jul 2023 Deviance: 3.0144\n", + "Time: 16:23:08 Pearson chi2: 5.00\n", + "No. Iterations: 6 Pseudo R-squ. (CS): 0.04355\n", + "Covariance Type: nonrobust \n", + "==============================================================================\n", + " coef std err z P>|z| [0.025 0.975]\n", + "------------------------------------------------------------------------------\n", + "const 5.0850 7.477 0.680 0.496 -9.570 19.740\n", + "x1 -0.1156 0.115 -1.004 0.316 -0.341 0.110\n", + "==============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 77, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "import statsmodels.api as sm\n", + "\n", + "data[\"Success\"]=data.Count-data.Malfunction\n", + "data[\"Intercept\"]=1\n", + "\n", + "freqs = data['Frequency'].values\n", + "\n", + "logmodel=sm.GLM(freqs, temperature_with_intercept, \n", + " #family=sm.families.Binomial(sm.families.links.logit())).fit()\n", + " family=sm.families.Binomial(sm.families.links.Logit())).fit()\n", + "\n", + "logmodel.summary()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The maximum likelyhood estimator of the intercept and of Temperature are thus $\\hat{\\alpha}=5.0849$ and $\\hat{\\beta}=-0.1156$. This **corresponds** to the values from the article of Dalal *et al.* The standard errors are $s_{\\hat{\\alpha}} = 7.477$ and $s_{\\hat{\\beta}} = 0.115$, which is **different** from the $3.052$ and $0.04702$ reported by Dallal *et al.* The deviance is $3.01444$ with 21 degrees of freedom. I cannot find any value similar to the Goodness of fit ($G^2=18.086$) reported by Dalal *et al.* There seems to be something wrong. Oh I know, I haven't indicated that my observations are actually the result of 6 observations for each rocket launch. Let's indicate these weights (since the weights are always the same throughout all experiments, it does not change the estimates of the fit but it does influence the variance estimates)." + ] + }, + { + "cell_type": "code", + "execution_count": 78, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
Generalized Linear Model Regression Results
Dep. Variable: y No. Observations: 23
Model: GLM Df Residuals: 21
Model Family: Binomial Df Model: 1
Link Function: Logit Scale: 1.0000
Method: IRLS Log-Likelihood: -23.526
Date: Fri, 28 Jul 2023 Deviance: 18.086
Time: 16:24:02 Pearson chi2: 30.0
No. Iterations: 6 Pseudo R-squ. (CS): 0.2344
Covariance Type: nonrobust
\n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "\n", + " \n", + "\n", + "
coef std err z P>|z| [0.025 0.975]
const 5.0850 3.052 1.666 0.096 -0.898 11.068
x1 -0.1156 0.047 -2.458 0.014 -0.208 -0.023
" + ], + "text/latex": [ + "\\begin{center}\n", + "\\begin{tabular}{lclc}\n", + "\\toprule\n", + "\\textbf{Dep. Variable:} & y & \\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:} & Fri, 28 Jul 2023 & \\textbf{ Deviance: } & 18.086 \\\\\n", + "\\textbf{Time:} & 16:24:02 & \\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{const} & 5.0850 & 3.052 & 1.666 & 0.096 & -0.898 & 11.068 \\\\\n", + "\\textbf{x1} & -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": [ + "\n", + "\"\"\"\n", + " Generalized Linear Model Regression Results \n", + "==============================================================================\n", + "Dep. Variable: y 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: Fri, 28 Jul 2023 Deviance: 18.086\n", + "Time: 16:24:02 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", + "const 5.0850 3.052 1.666 0.096 -0.898 11.068\n", + "x1 -0.1156 0.047 -2.458 0.014 -0.208 -0.023\n", + "==============================================================================\n", + "\"\"\"" + ] + }, + "execution_count": 78, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "logmodel = sm.GLM(freqs, temperature_with_intercept, \n", + " #family=sm.families.Binomial(sm.families.links.logit())).fit()\n", + " family=sm.families.Binomial(sm.families.links.Logit()),\n", + " var_weights=data['Count'].values).fit()\n", + "\n", + "logmodel.summary()" + ] + }, + { + "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": 80, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# New temperatures\n", + "new_temperature_data = np.linspace(start=30, stop=90, num=121)\n", + "new_temperature_with_intercept = sm.add_constant(new_temperature_data, prepend=True)\n", + "\n", + "# Predictions using previous model\n", + "predictions = logmodel.predict(new_temperature_with_intercept)" + ] + }, + { + "cell_type": "code", + "execution_count": 81, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAG2CAYAAACtaYbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABUaklEQVR4nO3deVhU9f4H8PeZlR0EZBXFXVEExTTsmplrek1b1Mxd815Lf6lki5WidUsrM/NmWpbLdW8xyzSNSNSUREVM3BAEQWVTZIeZYeb8/iBHR1BngOGwvF/PM0/M95w55zMfzpX3PasgiqIIIiIiIonIpC6AiIiIGjeGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiIiISFIMI0RERCQphhEiIiKSlMVh5ODBgxg2bBh8fHwgCAJ27tz5wM9ERUWhW7duUKvVaNOmDdavX1+FUomIiKghsjiMFBUVISgoCCtXrjRr/uTkZAwdOhR9+/ZFXFwcZs+ejRdeeAH79u2zuFgiIiJqeITqPChPEAT88MMPGDFixD3nef3117F7927Ex8cbx5577jnk5uZi7969VV01ERERNRAKa68gOjoa/fv3NxkbNGgQZs+efc/PaDQaaDQa43uDwYCcnBy4ublBEARrlUpEREQ1SBRFFBQUwMfHBzLZvQ/GWD2MZGRkwNPT02TM09MT+fn5KCkpga2tbYXPLF68GIsWLbJ2aURERFQL0tLS0KxZs3tOt3oYqYp58+YhLCzM+D4vLw/NmzdHcnIyHB0da2w9Op0O+/fvR9++faFUKmtsuQ0Re2UZ9st87JX52CvzsVfms2avCgoK0LJlywf+7bZ6GPHy8kJmZqbJWGZmJpycnCrdKwIAarUaarW6wrirqyucnJxqrDadTgc7Ozu4ublxY30A9soy7Jf52CvzsVfmY6/MZ81e3Vreg06xsPp9RkJDQxEZGWkyFhERgdDQUGuvmoiIiOoBi8NIYWEh4uLiEBcXB6D80t24uDikpqYCKD/EMmHCBOP806dPx6VLl/Daa6/h/Pnz+Pzzz/HNN99gzpw5NfMNiIiIqF6zOIwcP34cXbt2RdeuXQEAYWFh6Nq1KxYsWAAASE9PNwYTAGjZsiV2796NiIgIBAUF4eOPP8ZXX32FQYMG1dBXICIiovrM4nNGHnvsMdzv1iSV3V31sccew8mTJy1dFRERNRB6vR46na7W1qfT6aBQKFBaWgq9Xl9r662PqtMrpVIJuVxe7Rrq5NU0RETUMIiiiIyMDOTm5tb6er28vJCWlsb7Uz1AdXvl4uICLy+vavWZYYSIiKzmVhDx8PCAnZ1drQUDg8GAwsJCODg43PdmW1T1XomiiOLiYmRlZQEAvL29q1wDwwgREVmFXq83BhE3N7daXbfBYIBWq4WNjQ3DyANUp1e3btGRlZUFDw+PKh+y4W+IiIis4tY5InZ2dhJXQtZ06/dbnXOCGEaIiMiqeM5Gw1YTv1+GESIiIpIUwwgRERFJimGEiIjoLpMmTYIgCBVeiYmJUpfWIPFqGiIiokoMHjwY69atMxlr2rSpyXutVguVSlWbZTVI3DNCRERUCbVaDS8vL5NXv379MHPmTMyePRvu7u7GR5vEx8fjiSeegIODAzw9PTF+/Hhcv37duKyioiJMmDABDg4O8Pb2xscff4zHHnsMs2fPNs4jCAJ27txpUoOLi4vJnc3T0tIwatQouLi4wNXVFcOHD0dKSopx+qRJkzBixAgsXboU3t7ecHNzw4wZM0yudNFoNHj99dfh5+cHtVqNdu3aYePGjRBFEW3atMHSpUtNaoiLi7P6XiGGESIiqjWiKKJYW1YrrxKt3vjz/R5jYqkNGzZApVLh8OHDWL16NXJzc/H444+ja9euOH78OPbu3YvMzEyMGjXK+JlXX30VBw4cwI8//ohff/0VUVFRiI2NtWi9Op0OgwYNgqOjIw4dOoTDhw/DwcEBgwcPhlarNc63f/9+JCUlYf/+/diwYQPWr19vEmgmTJiArVu3YsWKFTh37hxWrVoFe3t7CIKAKVOmVNgbtG7dOjz66KNo06ZN1RpmBh6mISKiWlOi0yNgwb5aX+/ZdwbBTmXZn7yff/4ZDg4OxvdPPPEEAKBt27b48MMPjeP/+c9/0LVrV7z//vvGsbVr18LPzw8JCQnw8fHB119/jU2bNqFfv34AygNNs2bNLKpn+/btMBgM+Oqrr4yX065btw4uLi6IiorCwIEDAQBNmjTBZ599Brlcjg4dOmDo0KGIjIzEtGnTkJCQgG+++QYRERHo378/AMDf3x/5+fkAyvesLFiwADExMejRowd0Oh22bNlSYW9JTWMYISIiqkTfvn2xatUq43t7e3uMGTMGISEhJvOdOnUK+/fvNwkutyQlJaGkpARarRY9e/Y0jru6uqJ9+/YW1XPq1CkkJibC0dHRZLy0tBRJSUnG9506dTK5E6q3tzdOnz4NoPyQi1wuR58+fSpdh4+PD4YOHYq1a9eiR48e2LVrFzQaDUaOHGlRrZZiGCEiolpjq5Tj7DuDrL4eg8GAgvwCODo5QiaTwVZp+W3K7e3tKz00YW9vb/K+sLAQw4YNwwcffFBhXm9vb7PPtRAEocLhpDvP9SgsLERISAg2b95c4bN3nlirVCorLNdgMAC4ffv2+3nhhRcwfvx4fPLJJ1i3bh1Gjx5t9bvoMowQEVGtEQTB4sMlVWEwGFCmksNOpbD6s2m6deuG77//Hv7+/lAoKn631q1bQ6lU4ujRo2jevDkA4ObNm0hISDDZQ9G0aVOkp6cb31+8eBHFxcUm69m+fTs8PDzg5ORUpVoDAwNhMBhw4MAB42Gauw0ZMgT29vZYtWoV9u7di4MHD1ZpXZbgCaxERETVMGPGDOTk5GDMmDE4duwYkpKSsG/fPkyePBl6vR4ODg6YOnUqXn31Vfz++++Ij4/HpEmTKoSkxx9/HJ999hlOnjyJ48ePY/r06SZ7OcaOHQt3d3cMHz4chw4dQnJyMqKiovDyyy/jypUrZtXq7++PiRMnYsqUKdi5c6dxGT/88INxHrlcjkmTJmHevHlo27YtQkNDa6ZR98EwQkREVA0+Pj44fPgw9Ho9Bg4ciMDAQMyePRsuLi7GwPHRRx+hd+/eGDZsGPr3749//OMfFc49+fjjj+Hn54fevXvj+eefx9y5c00Oj9jZ2eHgwYNo3rw5nn76aXTs2BFTp05FaWmpRXtKVq1ahWeffRYvvfQSOnTogH//+98me2AAYOrUqdBqtZg8eXI1OmM+HqYhIiK6y52Xwt4pKiqq0vG2bdtix44d91yeg4MDNm7ciI0bNxrHdu/ebTKPj48P9u0zvdIoNzfX5L2Xlxc2bNhgUd3Lly83eW9jY4Nly5Zh2bJlAMoPad26muaWq1evQqlUYsKECfdcV01iGCEiIiIA5TdEy87OxsKFCzFy5Eh4enrWynp5mIaIiIgAAFu3bkWLFi2Qm5trci8Va+OeESIiIgnc65CPlCZNmoRJkybV+nq5Z4SIiIgkxTBCRERWVZPPhaG6pyZ+vwwjRERkFbfukXH3ZaPUsNz6/d5951dL8JwRIiKyCrlcDhcXF2RlZQEov0/GrQe8WZvBYIBWq0VpaanV78Ba31W1V6Ioori4GFlZWXBxcTF5Ho6lGEaIiMhqvLy8AMAYSGqLKIooKSmBra1trQWg+qq6vXJxcTH+nquKYYSIiKxGEAR4e3vDw8PD5KFv1qbT6XDw4EE8+uij1Tp80BhUp1dKpbJae0RuYRghIiKrk8vlNfJHy5L1lZWVwcbGhmHkAepCr3ggjYiIiCTFMEJERESSYhghIiIiSTGMEBERkaQYRoiIiEhSDCNEREQkKYYRIiIikhTDCBEREUmKYYSIiIgkxTBCREREkmIYISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQERGRpBhGiIiISFIMI0RERCQphhEiIiKSFMMIERERSarRhxG9QeoKiIiIGrdGHUa2HbuCT8/IkZFfKnUpREREjVajDSPF2jKs+D0RlwsFPLXqTxxLyZG6JCIiokap0YYRO5UC26b1gLediOuFWoz58k9s/PMyRFGUujQiIqJGpdGGEQBo7mqHOZ31GNLZE2UGEfN3xuON709DW8YTSYiIiGpLow4jAKCWA8tHdcEbT3SATAC2H0/DuK+O4kahRurSiIiIGoVGH0YAQBAETO/TGmsnPQRHtQIxKTkYvvIwLmQUSF0aERFRg8cwcofH2nvghxm90MLNDlduluDpzw/j9/OZUpdFRETUoDGM3KWNhyN2vvQIHm7liiKtHi9sOI4NR1KkLouIiKjBYhipRBN7FTZO7YnR3f1gEIHwn87gnV1noTfwShsiIqKaxjByD0q5DEueCcSrg9oDANYeTsb0TSdQrC2TuDIiIqKGpUphZOXKlfD394eNjQ169uyJmJiY+86/fPlytG/fHra2tvDz88OcOXNQWlr373oqCAJm9G2D/47pCpVChoizmRizhlfaEBER1SSLw8j27dsRFhaG8PBwxMbGIigoCIMGDUJWVlal82/ZsgVvvPEGwsPDce7cOXz99dfYvn073nzzzWoXX1uGBflgyws94WKnxKm0XDy7OhqpN4qlLouIiKhBsDiMLFu2DNOmTcPkyZMREBCA1atXw87ODmvXrq10/iNHjuCRRx7B888/D39/fwwcOBBjxox54N6Uuqa7vyu+m94Lvi62SL5ehKdXHUH81TypyyIiIqr3FJbMrNVqceLECcybN884JpPJ0L9/f0RHR1f6mV69emHTpk2IiYlBjx49cOnSJezZswfjx4+/53o0Gg00mtuHQvLz8wEAOp0OOp3OkpLv69ayzF1miyZqbJ/2EF7YeBLnMwow+otofP58MHq1dquxmuoqS3vV2LFf5mOvzMdemY+9Mp81e2XuMgXRgoexXLt2Db6+vjhy5AhCQ0ON46+99hoOHDiAo0ePVvq5FStWYO7cuRBFEWVlZZg+fTpWrVp1z/UsXLgQixYtqjC+ZcsW2NnZmVuu1ZSWAV8nyJCQJ4NcEDGxrQFBbrzShoiI6E7FxcV4/vnnkZeXBycnp3vOZ9GekaqIiorC+++/j88//xw9e/ZEYmIiZs2ahXfffRfz58+v9DPz5s1DWFiY8X1+fj78/PwwcODA+34ZS+l0OkRERGDAgAFQKpUWfXZomQGvfPsX9p3NwvqLcvwnoBNGhvjWWG11TXV61RixX+Zjr8zHXpmPvTKfNXt168jGg1gURtzd3SGXy5GZaXpX0szMTHh5eVX6mfnz52P8+PF44YUXAACBgYEoKirCv/71L7z11luQySqetqJWq6FWqyuMK5VKq2xUVVmuUgl8Pq473txxGtuPp+HNnWdQqNXjX4+2rvH66hJr/Q4aKvbLfOyV+dgr87FX5rNGr8xdnkUnsKpUKoSEhCAyMtI4ZjAYEBkZaXLY5k7FxcUVAodcLgcAWHCEqE6SywQseSYQ/+7TCgDw/p7zWBaRUO+/FxERUW2y+DBNWFgYJk6ciO7du6NHjx5Yvnw5ioqKMHnyZADAhAkT4Ovri8WLFwMAhg0bhmXLlqFr167GwzTz58/HsGHDjKGkPhMEAfOe6AhnWyU+3HsBKyIvokRbhjeHdIQgCFKXR0REVOdZHEZGjx6N7OxsLFiwABkZGQgODsbevXvh6ekJAEhNTTXZE/L2229DEAS8/fbbuHr1Kpo2bYphw4bhvffeq7lvUQe89Fgb2CrlWLTrLNYcSkaxVo93h3eGTMZAQkREdD9VOoF15syZmDlzZqXToqKiTFegUCA8PBzh4eFVWVW9MvmRlrBTyfHGjtPYfDQVmjIDPnimC+QMJERERPfEZ9PUsNEPNcfy0cGQywR8d+IKXv32FB+wR0REdB8MI1YwPNgX/x3TFQqZgB0nryLsmziU6Q1Sl0VERFQnMYxYyZBAb3z2fDcoZAJ+jLuGOd+cYiAhIiKqBMOIFQ3u7IXPx3aDUi5g16lrmL2de0iIiIjuxjBiZQM7eWHV2BAo5QJ+/isdc3kOCRERkQmGkVrQP8DTeMhmZ9w1vP79XzAwkBAREQFgGKk1gzp5YcWYrsarbN784TQDCRERERhGatWQQG98MjoYMgHYdiwNi3ad4a3jiYio0WMYqWVPBvlg6cggAMCG6Mv4cN8FiSsiIiKSFsOIBJ7u1gz/GdEZALAqKgkr9ydKXBEREZF0GEYkMu7hFnhzSAcAwEf7LmDd4WSJKyIiIpIGw4iE/vVoa8zq1xYAsGjXWXx/4orEFREREdU+hhGJze7fFlP/0RIA8Nr3fyHibKbEFREREdUuhhGJCYKAt4Z0xDPdmkFvEDFjSyz+vHRD6rKIiIhqDcNIHSCTCfjgmUD07+gJbZkBL2w4jvireVKXRUREVCsYRuoIhVyGz57vip4tXVGoKcPEtTFIuV4kdVlERERWxzBSh9go5fhqYnd08nHCjSItJq6LQXaBRuqyiIiIrIphpI5xtFFi3eSH4Odqi8s3ijFl/TEUasqkLouIiMhqGEbqIA9HG/xvSk+42qtw+moeXtx0Atoyg9RlERERWQXDSB3V0t0e6yY9BFulHIcuXscb3//F59gQEVGDxDBShwX5ueDzcd0glwnYcfIqPolIkLokIiKiGscwUsf1be+B958qf47Nit8TsS0mVeKKiIiIahbDSD0w+qHmePnxNgCAt3bGI+pClsQVERER1RyGkXpizoB2eLqrb/ldWjfH4sw13hSNiIgaBoaRekIQBCx5pgseaeOGIq0eU9cfR0ZeqdRlERERVRvDSD2iUsjw+dgQtPFwQEZ+KaZuOIYi3oOEiIjqOYaResbZVol1kx6Cm70KZ67lY9a2OOgNvOSXiIjqL4aResjP1Q5fTugOlUKG385l4v0956QuiYiIqMoYRuqpkBZN8PHIIADA138kY8tRXvJLRET1E8NIPTYsyAdhA9oBABb8GI8jSdclroiIiMhyDCP13P893gbDgnxQZhDx0uZYpFwvkrokIiIiizCM1HOCIOCjZ7sgqJkzcot1mLrhGPJKdFKXRUREZDaGkQbARinHmgnd4eVkg6TsIvzf1pMo0/Mpv0REVD8wjDQQHk42+Gpid9goZTiYkI0P9p6XuiQiIiKzMIw0IJ19nbH07yts1hxKxvcnrkhcERER0YMxjDQw/+zig5l9yx+qN++H04hLy5W2ICIiogdgGGmAwga0Q/+OntCWGfCv/x1HZj6fYUNERHUXw0gDJJMJ+GR0ENp6OCCrQIN/bzwBTZle6rKIiIgqxTDSQDnaKLFmQnc42SgQl5aLhT+dkbokIiKiSjGMNGD+7vZYMaYrBAHYGpOGzUcvS10SERFRBQwjDdxj7T0wd2B7AMDCn87gxOUciSsiIiIyxTDSCLz0WGsMCfSCTi9i+qZYntBKRER1CsNII1B+y/ggtPN0QHaBBi9tjoW2jHdoJSKiuoFhpJGwVyvw5fjucLRR4MTlm3h/zzmpSyIiIgLAMNKo+Lvb45NRwQCA9UdS8MNJ3qGViIikxzDSyPQP8MT/Pf73HVp3nMa59HyJKyIiosaOYaQRmt2/HR5t1xSlOgOmbzqBvBKd1CUREVEjxjDSCMllAj4dHQxfF1tcvlGMV745BYNBlLosIiJqpBhGGqkm9iqsGtcNKrkMv53LxBcHL0ldEhERNVIMI41Yl2YuWPhkJwDAR/vOIzrphsQVERFRY8Qw0siN6eGHp7v5wiAC/7f1JLJ4QzQiIqplDCONnCAIeG9EIDp4OeJ6oQYzt5xEmZ43RCMiotrDMEKwVcmxalwIHNQKxKTk4OOIBKlLIiKiRoRhhAAALd3t8eGzXQAAq6KS8Pv5TIkrIiKixoJhhIyGBHpjUi9/AMCc7adw5WaxtAUREVGjwDBCJuYN6YCgZs7IK9Fh5paTfKAeERFZHcMImVAr5Pjs+W5wslEgLi0XS345L3VJRETUwDGMUAV+rnb4+O8H6q09nIxfz2RIWxARETVoDCNUqQEBnnjhHy0BAHO/PYWruSUSV0RERA0Vwwjd02uDOyDIzwX5pWWYtf0v8PQRIiKyBoYRuieVQobPxnSFk40Cp67k4edUbi5ERFTz+NeF7svP1Q4fjQwCAOxPl+H3C9kSV0RERA0Nwwg90KBOXpjwcHMAwBs74pGex/NHiIio5lQpjKxcuRL+/v6wsbFBz549ERMTc9/5c3NzMWPGDHh7e0OtVqNdu3bYs2dPlQomabw2qB2a2Yu4WazDrG1xfH4NERHVGIvDyPbt2xEWFobw8HDExsYiKCgIgwYNQlZWVqXza7VaDBgwACkpKfjuu+9w4cIFrFmzBr6+vtUunmqPWiHDxLZ62KvkiEnOwX9/T5S6JCIiaiAsDiPLli3DtGnTMHnyZAQEBGD16tWws7PD2rVrK51/7dq1yMnJwc6dO/HII4/A398fffr0QVBQULWLp9rlYQu882QAAOC/v19EdNINiSsiIqKGQGHJzFqtFidOnMC8efOMYzKZDP3790d0dHSln/npp58QGhqKGTNm4Mcff0TTpk3x/PPP4/XXX4dcLq/0MxqNBhqNxvg+Pz8fAKDT6aDT6Swp+b5uLasml9lQ3erREwHuiO7mi+9ir2LWtpP4aUYo3OxVEldX93DbMh97ZT72ynzslfms2Stzl2lRGLl+/Tr0ej08PT1Nxj09PXH+fOW3Db906RJ+//13jB07Fnv27EFiYiJeeukl6HQ6hIeHV/qZxYsXY9GiRRXGf/31V9jZ2VlSslkiIiJqfJkNVUREBHoogEO2cmQWaDBl9e/4VwcDBEHqyuomblvmY6/Mx16Zj70ynzV6VVxs3gNXLQojVWEwGODh4YEvv/wScrkcISEhuHr1Kj766KN7hpF58+YhLCzM+D4/Px9+fn4YOHAgnJycaqw2nU6HiIgIDBgwAEqlssaW2xDd3auODxXgmS+O4mwukO3aEZNCW0hdYp3Cbct87JX52CvzsVfms2avbh3ZeBCLwoi7uzvkcjkyMzNNxjMzM+Hl5VXpZ7y9vaFUKk0OyXTs2BEZGRnQarVQqSru4ler1VCr1RXGlUqlVTYqay23IbrVq0A/V7w9tCMW/HgGH+27iNDWTdHZ11nq8uocblvmY6/Mx16Zj70ynzV6Ze7yLDqBVaVSISQkBJGRkcYxg8GAyMhIhIaGVvqZRx55BImJiTAYbl8KmpCQAG9v70qDCNUf4x9ugQEBntDqDXh560kUacqkLomIiOohi6+mCQsLw5o1a7BhwwacO3cOL774IoqKijB58mQAwIQJE0xOcH3xxReRk5ODWbNmISEhAbt378b777+PGTNm1Ny3IEkIgoAPn+kCb2cbXLpehPCfzkhdEhER1UMWnzMyevRoZGdnY8GCBcjIyEBwcDD27t1rPKk1NTUVMtntjOPn54d9+/Zhzpw56NKlC3x9fTFr1iy8/vrrNfctSDJN7FVYPjoYY9b8ie9OXMGj7ZriySAfqcsiIqJ6pEonsM6cORMzZ86sdFpUVFSFsdDQUPz5559VWRXVAz1buWFm3zZY8Xsi3tpxGl39XODnWvNXPRERUcPEZ9NQjXi5X1uEtGiCAk0ZZm07ydvFExGR2RhGqEYo5DIsHx0MR7UCsam5WBF5UeqSiIionmAYoRrj52qH954OBAB8tj8RRy/xdvFERPRgDCNUo54M8sGzIc1gEIE52+OQV8xbMRMR0f0xjFCNW/RkJ/i72eFaXine3HkaoihKXRIREdVhDCNU4+zVCnz6XFcoZAJ2/5WOb09ckbokIiKqwxhGyCqC/FwwZ0A7AMDCn84g+XqRxBUREVFdxTBCVjO9T2v0bOmKYq0es7edhI6X+xIRUSUYRshq5DIBn4wOhrOtEqeu5OGTiASpSyIiojqIYYSsysfFFov/vtx31YEkXu5LREQVMIyQ1Q0J9Mao7s0g3rrct4SX+xIR0W0MI1Qrwofdvtz3rR94uS8REd3GMEK1wl6twPLnukIuE/DzX+n44eRVqUsiIqI6gmGEak2wnwvm9G8LAFjw4xmk3iiWuCIiIqoLGEaoVr34WBs85N8EhZoyzN7Op/sSERHDCNUyuUzAslG3n+67cn+S1CUREZHEGEao1vm52uGdEZ0AACt+v4jY1JsSV0RERFJiGCFJjAj2xbAgH+gNIuZsj0OhpkzqkoiISCIMIyQJQRDwnxGd4etii8s3ivHOrjNSl0RERBJhGCHJONsq8fGoIAgC8M3xK9gbny51SUREJAGGEZLUw63cML1PawDAGztOIzO/VOKKiIiotjGMkOTm9G+HTj5OyC3WYe63p2Aw8O6sRESNCcMISU6lkOHT54KhVshw6OJ1bIhOkbokIiKqRQwjVCe08XDE20M7AgAW/3IeFzIKJK6IiIhqC8MI1RnjHm6Bvu2bQltmwKxtJ6Ep00tdEhER1QKGEaozBEHAh88Gwc1ehfMZBfj41wSpSyIiolrAMEJ1SlNHNZY80wUAsObQJRxJvC5xRUREZG0MI1TnDAjwxJgezSGKQNg3p5BXrJO6JCIisiKGEaqT5v+zI1q62yMjvxRv7jwNUeTlvkREDRXDCNVJdioFlo8OhlwmYPdf6fjh5FWpSyIiIithGKE6K8jPBbP7tQUAhP94Bmk5xRJXRERE1sAwQnXai4+1RkiLJijQlOGVb05Bz7uzEhE1OAwjVKcp5DIsHx0MB7UCMSk5WH0gSeqSiIiohjGMUJ3n52qHhU92AgB8EpGA01fyJK6IiIhqEsMI1QvPdPPFkEAvlBlEzNp+EiVa3p2ViKihYBihekEQBLw3IhCeTmpcyi7Ce3vOSl0SERHVEIYRqjea2KuwdGQQAGDTn6n4/XymxBUREVFNYBiheqV326aY+o+WAIDXvvsL1ws1EldERETVxTBC9c6rg9qjvacjrhdq8fp3f/HurERE9RzDCNU7Nko5lj8XDJVchsjzWdh8NFXqkoiIqBoYRqhe6ujthNcGtwcA/Gf3WSRmFUpcERERVRXDCNVbUx5piX+0cUepzoDZ209CW2aQuiQiIqoChhGqt2QyAUtHBsHFTon4q/n45LcEqUsiIqIqYBihes3L2QaLnwoEAKw+kISjl25IXBEREVmKYYTqvScCvTEypBlEEQj75hTySnRSl0RERBZgGKEGIfzJTmjhZoeruSWYvzNe6nKIiMgCDCPUIDioFfhkdDDkMgE/nbqGnSevSl0SERGZiWGEGoxuzZtgVr+2AID5O+ORllMscUVERGQOhhFqUF56rDW6t2iCAk0Z5myPQ5mel/sSEdV1DCPUoCjkMnwyOhiOagWOX76JlfuTpC6JiIgegGGEGhw/Vzu8M6ITAGDF7xdx4vJNiSsiIqL7YRihBmlEsC+eDPKB3iBi9vaTKCjl5b5ERHUVwwg1SIIg4N0RneHrYou0nBKE/3hG6pKIiOgeGEaowXK2VeLT54IhE4AdJ6/ixzhe7ktEVBcxjFCD1t3fFTMfL7/c9+0feLkvEVFdxDBCDd7Lj7dBt+YuvNyXiKiOYhihBk8hl+HT57oaL/f9bH+i1CUREdEdGEaoUfBztcN/nuoMAFgReRHHUnIkroiIiG5hGKFGY3iwL57u6guDCMzeFsen+xIR1REMI9SovDOis/Hpvm/+cBqiKEpdEhFRo8cwQo2Kg1qBFc91hUImYPdf6fj2+BWpSyIiavQYRqjRCfJzwSsD2wMAwn86g6TsQokrIiJq3BhGqFH696Ot8EgbN5To9Pi/LSehKdNLXRIRUaNVpTCycuVK+Pv7w8bGBj179kRMTIxZn9u2bRsEQcCIESOqslqiGiOTCVg2Khiu9iqcTc/Hkl/OS10SEVGjZXEY2b59O8LCwhAeHo7Y2FgEBQVh0KBByMrKuu/nUlJSMHfuXPTu3bvKxRLVJE8nG3w8MggAsO5wCiLPZUpcERFR42RxGFm2bBmmTZuGyZMnIyAgAKtXr4adnR3Wrl17z8/o9XqMHTsWixYtQqtWrapVMFFN6tvBA1P/0RIAMPfbU8jIK5W4IiKixkdhycxarRYnTpzAvHnzjGMymQz9+/dHdHT0PT/3zjvvwMPDA1OnTsWhQ4ceuB6NRgONRmN8n5+fDwDQ6XTQ6Wru3hC3llWTy2yoGnKv5vRrjeik6zibXoBZ22KxYVJ3yGVCtZbZkPtV09gr87FX5mOvzGfNXpm7TIvCyPXr16HX6+Hp6Wky7unpifPnKz/m/scff+Drr79GXFyc2etZvHgxFi1aVGH8119/hZ2dnSUlmyUiIqLGl9lQNdRePeUJJGbKcTT5Juas2YvBfjVz/5GG2i9rYK/Mx16Zj70ynzV6VVxs3sNJLQojliooKMD48eOxZs0auLu7m/25efPmISwszPg+Pz8ffn5+GDhwIJycnGqsPp1Oh4iICAwYMABKpbLGltsQNYZeubS+hle/j8e+q3KMG9QdPfxdq7ysxtCvmsJemY+9Mh97ZT5r9urWkY0HsSiMuLu7Qy6XIzPT9ES/zMxMeHl5VZg/KSkJKSkpGDZsmHHMYCh/YqpCocCFCxfQunXrCp9Tq9VQq9UVxpVKpVU2KmsttyFqyL0a+VAL/Jmci+9jr+CVb+OxZ1ZvuNqrqrXMhtyvmsZemY+9Mh97ZT5r9Mrc5Vl0AqtKpUJISAgiIyONYwaDAZGRkQgNDa0wf4cOHXD69GnExcUZX08++ST69u2LuLg4+Pn5WbJ6Iqt7Z3gntGpqj4z8Urz67SneLp6IqBZYfJgmLCwMEydORPfu3dGjRw8sX74cRUVFmDx5MgBgwoQJ8PX1xeLFi2FjY4POnTubfN7FxQUAKowT1QX2agU+G9MNIz4/jMjzWfj6j2S80JtXgBERWZPFYWT06NHIzs7GggULkJGRgeDgYOzdu9d4UmtqaipkMt7YleqvAB8nzP9nAObvjMcHe8+ju78rgv1cpC6LiKjBqtIJrDNnzsTMmTMrnRYVFXXfz65fv74qqySqVeN6Nkd00nXsOZ2BGZtjsefl3nC243FnIiJr4C4MokoIgoAlz3RBc1c7XM0twavf8fwRIiJrYRghugcnGyVWPt8NKrkMv57NxPojKVKXRETUIDGMEN1HYDNnvDmkAwDg/T3ncCotV9qCiIgaIIYRogeY2Msfgzt5QacXMWNLLPKKeXtpIqKaxDBC9ACCIOCDZ8vPH7lyswSvfBvH80eIiGoQwwiRGZxtlfh8bPn5I7+dy8KaQ5ekLomIqMFgGCEyU2dfZywYFgAA+GDvBRxPyZG4IiKihoFhhMgCY3s2x/BgH+gNImZuOYnrhRqpSyIiqvcYRogsIAgC3n8qEK3/fn7NrG0noTfw/BEioupgGCGykL1agVXjQmCrlONw4g18EpEgdUlERPUawwhRFbTzdMSSZwIBAJ/tT8Tv5zMlroiIqP5iGCGqouHBvpgQ2gIAMGf7KaTlFEtcERFR/cQwQlQNbw3tiGA/F+SV6PDi5hMo1emlLomIqN5hGCGqBrVCjs/HdoOrvQrxV/Ox4Md4q90QTW8QEZ10Az/GXUV00g2eOEuS4zZJNUUhdQFE9Z2Piy3+O6Yrxn99FN8cv4JgvyYY2c27RtexNz4di3adRXpeqXHM29kG4cMCMLhzza6LyBzcJqkmcc8IUQ14pI075g5qDwBY+NMZnLqSV2PL3hufjhc3xZr8ow8AGXmleHFTLPbGp9fYuojMwW2SahrDCFENebFPawzq5Amt3oCZW+NQWAPP09MbRCzadRaV7fy+NbZo11nuHqdaw22SrIFhhKiGCIKApSOD0KqpPTLyNViXIEOZ3lCtZcYk51T4f593EgGk55UiJpm3pqfawW2SrIFhhKgGOdoo8cW4ENir5EjMl+GDfdW7IVpWwb3/0a/KfETVxW2SrIFhhKiGtfV0xIfPdAYArI9OxY7YK1VeloejTY3OR1Rd3CbJGhhGiKxgYIAnBvqWH6KZt+M04q9W7YTWHi1d4e1sA+Ee0wWUX8HQo6Vr1QolshC3SbIGhhEiK3nCz4DH2rlDU2bAvzeewI0qPOFXLhMQPiwAACr843/rffiwAMhl9/rTQFSzuE2SNTCMEFmJTAA+fjYQLd3tcTW3BC9ujoW2zPITWgd39saqcd3g5Wy629vL2QarxnXjPR2o1nGbpJrGm54RWZGTrRJrJoRgxMojiEnOwaJdZ/DeU4EWL2dwZ28MCPBCTHIOsgpK4eFYvhuc/++TpMJtkmoSwwiRlbXxcMSnzwXjhf8dx+ajqejo7YRxD7eweDlymYDQ1m5WqJCoarhNUk3hYRqiWtCvoyfmDrx9h9Y/L92QuCIiorqDYYSolrz0WGsMC/JBmUHES5tjkZZTLHVJRER1AsMIUS0RBAEfPtMFnX2dkFOkxdQNx1BQWgP3jCciqucYRohqka1KjjUTusPDUY2EzELM2hbHZ3gQUaPHMEJUy7ydbfHlhO5QK2T4/XwWlvxyTuqSiIgkxTBCJIFgPxcsHRkEAFhzKBnbj6VKXBERkXQYRogkMizIB7P6tQUAvPVDPI4kXZe4IiIiaTCMEElodv+2+GcXb5QZREzfeAKJWYVSl0REVOsYRogkJAgClo4MQrfmLsgvLcOU9ceq9AwbIqL6jGGESGI2yvIrbPxcbZGaU4x/bTyBUp1e6rKIiGoNwwhRHeDmoMa6SQ/B0UaBE5dvYu63p2DgJb9E1EgwjBDVEW08HPHFuBAoZAJ+/isdH+67IHVJRES1gmGEqA7p1cYdHzzTBQCw+kASNv55WeKKiIisj2GEqI55JqQZXhnQDgAQ/mM8fjubKXFFRETWxTBCVAfNfLwNRnf3g0EE/m/rScSl5UpdEhGR1TCMENVBgiDgP091Rp92TVGi02PK+mO4lM17kBBRw8QwQlRHKeUyfD62G7o0c0ZOkRYT1sYgK79U6rKIiGocwwhRHWavVmDtpIfg72aHKzdLMHHdMeSX6qQui4ioRjGMENVx7g5q/G9KT7g7qHAuPR///h9vikZEDQvDCFE90NzNDusn94C9So7oSzcwe1scyvQGqcsiIqoRDCNE9URnX2esmdAdKrkMe89k4K0f4iGKvEsrEdV/DCNE9UivNu5YMaYrZAKw/XgalvxyXuqSiIiqjWGEqJ4Z3NkLS54uv0vrFwcvYVVUksQVERFVD8MIUT006iE/vDmkAwDgg73nsTE6RdqCiIiqgWGEqJ7616OtMbNvGwDA/B/P4LsTVySuiIioahhGiOqxVwa2w+RH/AEAr313CntOp0tbEBFRFTCMENVjgiBgwT8DjM+xeXnrSfx+ng/WI6L6hWGEqJ4TBAHvPx2IYUE+KDOImL4xFgcSsqUui4jIbAwjRA2AXCZg2aggDO7kBa3egH/97zgOJ16XuiwiIrMwjBA1EEq5DCvGdEX/jh7QlBkwdcMxRCfdkLosIqIHYhghakBUChlWju2Gvu2bolRnwJT1x3D0EgMJEdVtDCNEDYxaIceqcSF4tF1TlOj0mLSOe0iIqG5jGCFqgGyUcnw5/nYgmbw+hueQEFGdxTBC1EDdCiR3HrI5dJFX2RBR3cMwQtSA2SjlWD0+BP063Dqp9TjvQ0JEdQ7DCFEDd+sckoEBntCWGfDvjSd4p1YiqlMYRogagVtX2QwL8oFOL2Lmllh8z2fZEFEdwTBC1Ego5TIsHx2MUd2bwSACr3x7Cpv+vCx1WUREVQsjK1euhL+/P2xsbNCzZ0/ExMTcc941a9agd+/eaNKkCZo0aYL+/fvfd34ish65TMCSp7tgUi9/AMDbO+Oxcn8iRFGUtjAiatQsDiPbt29HWFgYwsPDERsbi6CgIAwaNAhZWVmVzh8VFYUxY8Zg//79iI6Ohp+fHwYOHIirV69Wu3gispxMJiB8WABm9m0DAPho3wW8t/scDAYGEiKShsVhZNmyZZg2bRomT56MgIAArF69GnZ2dli7dm2l82/evBkvvfQSgoOD0aFDB3z11VcwGAyIjIysdvFEVDWCIGDuoPZ4e2hHAMBXfyTj1e/+QpneIHFlRNQYKSyZWavV4sSJE5g3b55xTCaToX///oiOjjZrGcXFxdDpdHB1db3nPBqNBhqNxvg+Pz8fAKDT6aDT6Swp+b5uLasml9lQsVeWqS/9mviwHxzVMry58yy+j72Cm0UaLB/VBbYqea3VUF96VRewV+Zjr8xnzV6Zu0xBtOBg8bVr1+Dr64sjR44gNDTUOP7aa6/hwIEDOHr06AOX8dJLL2Hfvn04c+YMbGxsKp1n4cKFWLRoUYXxLVu2wM7OztxyichMp3MErE+QoUwU4O8gYloHPRyUUldFRPVdcXExnn/+eeTl5cHJyeme81m0Z6S6lixZgm3btiEqKuqeQQQA5s2bh7CwMOP7/Px847km9/syltLpdIiIiMCAAQOgVPJf3vthryxT3/o1BEC/yzcxffNJpBSW4esUJ3w9IQTNmthafd31rVdSYq/Mx16Zz5q9unVk40EsCiPu7u6Qy+XIzDS9g2NmZia8vLzu+9mlS5diyZIl+O2339ClS5f7zqtWq6FWqyuMK5VKq2xU1lpuQ8ReWaY+9Su0jQe+m94LE9fG4NL1YoxaE4N1kx5CZ1/nWll/feqV1Ngr87FX5rNGr8xdnkUnsKpUKoSEhJicfHrrZNQ7D9vc7cMPP8S7776LvXv3onv37paskohqUVtPR+x46RF08HJEdoEGo76IRuQ53j6eiKzL4qtpwsLCsGbNGmzYsAHnzp3Diy++iKKiIkyePBkAMGHCBJMTXD/44APMnz8fa9euhb+/PzIyMpCRkYHCwsKa+xZEVGO8nG3wzfRQ9G7rjmKtHtP+dxzrDydLXRYRNWAWh5HRo0dj6dKlWLBgAYKDgxEXF4e9e/fC09MTAJCamor09NvPvVi1ahW0Wi2effZZeHt7G19Lly6tuW9BRDXKyUaJtZMewnMP+cEgAgt3ncWiXWeg571IiMgKqnQC68yZMzFz5sxKp0VFRZm8T0lJqcoqiEhiSrkMi58ORAs3e3yw9zzWHU5B8vUirBjTFU42PAZPRDWHz6YhonsSBAH/erQVZvdrC6VcQNSFbDy98ghSrhcZ59EbREQn3cCPcVcRnXSjTuw9qYs1VZe2zICN0SkAgI3RKdCW8QZ11HDU6qW9RFS/7I1Px6JdZ5GeV2ocS8wuxNAVh7BmQnfkl+oqTPd2tkH4sAAM7uwtRcmV1ix1TdW1eM9ZrDmUDKVMxIc9gA/2XcB/fknAtN4tMW9IgNTlEVUb94wQUaX2xqfjxU2xJn/UbynS6jH2q6OYXsn0jLxSvLgpFnvj0yt8ztruVbOUNVXX4j1n8cXBZNy9c8cgAl8cTMbiPWelKYyoBjGMEFEFeoOIRbvO4n4HN+417db4ol1na/XwyP1qlqqm6tKWGbDm0P2vZFpzKJmHbKjeYxghogpiknMq3SNiLhFAel4pYpJzaq6oB3hQzVLUVF0bo1Mq7BG5m0GE8VwSovqKYYSIKsgqqHoQscZyanJdtVlTdV3OKa7R+YjqKoYRIqrAw/Hez46SYjk1ua7arKm6Wria92BQc+cjqqsYRoiogh4tXeHtbAPhPvPIBNx3uoejGj1autZ0aff0oJoFlF9VU5s1Vdf4UH/I7tdklP8exof610o9RNbCMEJEFchlAsKHlV8yevffQuHv17TeLSudfoumzIBjKbV3fsaDagaA8GEBkD/or3sdolLIjH2+l2m9W0Kl4D/lVL9xCyaiSg3u7I1V47rBy9n0sIaXsw1WjeuGeUMCKp3u7qCCl5MaeSU6PL/mTyz/LaHWrmB5UM318T4j84YE4N+Ptqywh0QmAP9+lPcZoYaBNz0jonsa3NkbAwK8EJOcg6yCUng4lh/muLV34V7TNWV6hP94Bt+euILlv13En5du4JPRwfB2tpW85vpo3pAAvDKwAzYdSQJunsXrg9pjXK/W3CNCDQbDCBHdl1wmILS1m0XT7VQKfDQyCL3auOGtH+Lx56UcDF5+CIufDsSQQOvvnXhQzfWRSiHD+FB/7NlzFuND/aFkEKEGhFszEVnNU12bYffLvdGlmTPySnR4aXMsXv32FAo1ZVKXRkR1CMMIEVlVS3d7fP9iL8zo2xqCAHx74gqGfHqoVk9uJaK6jWGEiKxOKZfh1UEdsG3aw/B1sUVqTjFGfRGN9/ecg0anl7o8IpIYwwgR1Zqerdzwy+zeGBnSDKIIfHnwEkas+hOphVJXRkRSYhgholrlZKPERyOD8NWE7nB3UCMxuwjLTsvx4b4ElHIvCVGjxDBCRJLoH+CJX+c8iqGBXhAhYM0fKTyXhKiRYhghIsm42quwfFQXvNBeDw9HNS5dL8KoL6Lx9s7TyC/VSV0eEdUShhEiklygq4hf/q8XRnf3gygCm/5MRb+PD2D3X+kQxdq5eysRSYdhhIjqBCdbJT54tgu2TOuJVu72yC7QYMaWWExefwyXbxRJXR4RWRHDCBHVKb1au+OX2b0xu39bqOQyRF3IxoBPDmLZrxdQouUJrkQNEcMIEdU5aoUcs/u3wy+ze6N3W3doywxY8XsiBnxyAPvOZPDQDVEDwzBCRHVW66YO+N+UHvh8bDf4ONvgys0S/HvjCYz7+ijOZ+RLXR4R1RCGESKq0wRBwJBAb/z2Sh/M6Fv+pNrDiTcw5NNDeHvnadwo1EhdIhFVE8MIEdULdioFXh3UAZFhffBEZy8Y/r7q5rGPorByfyJvmEZUjzGMEFG94udqh1XjQrDtXw+jk48TCjRl+GjfBfRdGoVvj6dBb+D5JET1DcMIEdVLD7dyw66Z/8Dy0cHwdbFFel4pXv3uLwxefhB743l/EqL6hGGEiOotmUzAiK6+iHylD94a0hHOtkpczCrE9E2xGL7yMA4mZDOUENUDDCNEVO/ZKOWY9mgrHHytL15+vA3sVHL8dSUPE9bGYOTqaPxx8TpDCVEdxjBCRA2Gs60SYQPb4+BrfTHlkZZQKWQ4fvkmxn19FCNXR+PQRe4pIaqLGEaIqMFxd1BjwbAAHHqtLyY/4m8MJeO/jsGIlYex70wGDDzRlajOYBghogbL08kG4cM6GUOJjVKGU1fy8O+NJzBo+UF8d+IKtGUGqcskavQYRoiowbsVSg6//jhm9G0NR7UCF7MKMffbU3j0w/344kAS8kt1UpdJ1GgxjBBRo+HmoMargzrg8LzH8drg9mjqqEZGfikW/3IevRb/jnd2nUXqjWKpyyRqdBhGiKjRcbJR4qXH2uCP1/viw2e7oK2HAwo1ZVh7OBl9lu7HCxuO40gir8Ahqi0KqQsgIpKKWiHHqO5+GBnSDAcSsrHucAoOJGTjt3OZ+O1cJlo3tcf4h1vg6ZBmcLJRSl0uUYPFMEJEjZ4gCHisvQcea++BxKxCbDiSgh2xV5CUXYSFu87ig70XMDzYB8/1aI6gZs4QBEHqkokaFIYRIqI7tPFwwLsjOuO1we2x8+RVbPzzMhIyC7HtWBq2HUtDBy9HPPeQH0Z09YWLnUrqcokaBIYRIqJKONooMT7UH+MeboFjKTexNSYVe06n43xGARbuOov395zHgABPPBvSDL3bukMh5yl4RFXFMEJEdB+CIKBHS1f0aOmKhcM6YWfcVWyNScX5jALsPp2O3afT4eGoxpNBPhjR1RedfJx4GIfIQgwjRERmcrZTYmIvf0zs5Y8z1/Lw3Ykr+DHuGrIKNPjqj2R89UcyWje1x4hgXwwL8oG/u73UJRPVCwwjRERV0MnHGZ18nDHviY6IupCFH+Ou4bdzmUjKLsLHEQn4OCIBnX2d8M8uPhga6A0/VzupSyaqsxhGiIiqQaWQYWAnLwzs5IWCUh1+ic/ArlPXcCTpBuKv5iP+aj6W/HIenX2dMLiTFwZ39kIbD0epyyaqUxhGiIhqiKONEqO6+2FUdz/cKNRg75kM/HwqHUeTbweTpb8moFVTewwI8MSAjp7o2rwJ5DKeY0KNG8MIEZEVuDmoMbZnC4zt2QI3CjWIOJuJvWcycDjxOi5lF+GLA5fwxYFLcLNXoU/7pni8gwd6t20KZ1veXI0aH4YRIiIrc3NQ47kezfFcj+bIL9XhYEI2Is5mYv/5LNwo0mJH7FXsiL0KuUxASIsm6NOuKXq3dUdnH2fIuNeEGgGGESKiWuRko8Q/u/jgn118oNMbcCwlB1EXsvH7+SwkZhUiJjkHMck5+GjfBTSxU+IfbZvikdZueKSNO7wcudeEGiaGESIiiSjlMvRq7Y5erd3x5pCOSMspRlRCNg4mZCM66QZuFuuw69Q17Dp1DQDQrIktmill0MVdQ6+2HvBxsZX4GxDVDIYRIqI6ws/VDuMfboHxD7eATm/AydRc/JF4HUcSryMuLRdXbpbgCmT48/t4AEBzVzv0aOmKh/yb4CF/V7R0t+cN16heYhghIqqDlHKZ8c6vYQPaoVBThujELGyNPIHrggvir+UjNacYqTnF+O7EFQCAm70K3Vo0QbfmTdCtuQu6NHOBrUou8TchejCGESKiesBBrcBj7ZqiONGAIUMeRqkeOH75Jo6n5OBY8k3EXcnFjSItIs5mIuJsJgBALhPQwcsRQX4uCG7mgi5+zmjT1IHP0aE6h2GEiKgecrRRom97D/Rt7wEA0JTpEX81D7GXcxGbehMnLt9EVoEGZ67l48y1fGw5mgoAsFHKEODthEBfZ3TydUYnHye09XCESsGAQtJhGCEiagDUCjlCWrgipIUrAEAURaTnleJUWi7iruQiLjUXZ67lo1BThtjUXMSm5ho/q5LL0NbTAR29ndDBy9H4XzcHtUTfhhobhhEiogZIEAT4uNjCx8UWTwR6AwAMBhHJN4oQfzUPf13Jw5lreTh7LR/5pWXGPSh3crNXoZ2nI9p7OaKtpwPaNHVAW09HuNqrpPhK1IAxjBARNRIymYDWTR3QuqkDhgf7Aijfg3LlZgnOXMvHhYwCnM/Ix7n0fFzOKcaNIi2iL91A9KUbJstxtVehdVN7tG7qgFZN7dHK3QH+7vZo7mrHwz1UJQwjRESNmCAI8HO1g5+rHQZ39jKOF2vLkJhViITMQlzIyEdiViEuZhXiys0S5BRpkVOkxbGUmybLkgmAbxNb+LvZw9/NHi3c7NDc1Q4t3Ozh52oLOxX/5FDluGUQEVEFdioFujQrvzz4TsXaMlzKLkJSdqHJf1NuFKFYq0daTgnSckpw6OL1Cst0s1ehmasd/JrYolkTO/g2sUWzJrbwdbGFt7MNHG14h9nGimGEiIjMZqdSoLOvMzr7OpuMi6KI7AINkq+XB5PLN4pxOacYl28UIS2nBHklOtwo0uJGkRan0nIrXbajjQI+zrbwdrGBt7MNvJzKQ4qHkxpezjbwcrKBs62SN3ZrgBhGiIio2gRBgIeTDTycbNCzlVuF6XklOqTlFOPKzeLyO8n+/bqaW4JrueVhpaC0DBdKC3Ahs+Ce61EpZPBwVP/9Kg8qTR3UaOqohruDGu6Oarg7qOCs5rkr9QnDCBERWZ2zrRLOlexRuaVIU4b0vBJczS1FRl4J0vNKkZ5bioz8UmT+/bpZrIO2zGAMMg+ilsmx9PwhuDmo4Wavguvfryb2KrjaqeBip0QTexWa2CnhYqeCi62SN4STCMMIERFJzl6tQBsPR7TxcLznPKU6PbILNMgqKEVWvgZZBRpk//26Xlj+/kahBtcLtdDqDdAYBKTdLEGaGcHlFke1As52yvLwZKuEy98/O9ko4WSrhJONAk62SjjaKOBoUz7uaKOAg40CDioFZDIeQqoKhhEiIgvoDSJiknOQVVAKD0cb9GjpCvnff4BKtHq8v+csUm4Uw9/NDm8OCTB5Nsz9Pnu/aQCgLTNgU3QK3ABsjE7BuF6tzb6Mtjrrrepyb9W8MToFl3OK0cLVDuND/atVs41Sbrzy59Z0FzslHmnjbrJuURRxs7AE3++OQJeHeiGvVI8bf18BdPPv81aSrxfhZpEWpWV6lGj1yC8tAwAUaMpQoCkza8/L3QQBcFD9HUzUCtirFXC0UcBeVf6zvVoOe3X5NDuVHPYqBWxVctir5bBVlk+3U8lhq1LATimHrUoOtULWKM6RqVIYWblyJT766CNkZGQgKCgI//3vf9GjR497zv/tt99i/vz5SElJQdu2bfHBBx9gyJAhVS6aiEgKe+PTsWjXWaTnlRrHvJ1tED4sAN/HXkHE2Szj+KGLwMY/UzEgwANrJjx0388CuOe0wZ29sXjPWaw5lAylTMSHPYAP9l3Af35JwLTeLTFvSECVa37Qequ63DtrNoi3P/PennPVrnlwZ+8HThcEAY42SnjYAt2au0CpVD5w2YufDkRoa3fkFmtxs1iH/BId8u545ZfokF+qQ35JWfn5LZryc1wKSsuQX6JDmUGEKN4OMzVFJgC2fwcTW5UcNorb/7VRyWGjkMFGKYeNsvy/auP78p/VChnUxp/vHJNBJZdDpZBBBj1uasrDo1KiC5osDiPbt29HWFgYVq9ejZ49e2L58uUYNGgQLly4AA8PjwrzHzlyBGPGjMHixYvxz3/+E1u2bMGIESMQGxuLzp0718iXICKytr3x6XhxUyzEu8Yz8koxfVPsPT8XcTYLT352CKev5Fv02Yy8Ury4KRb9AzxMQs4tBhH44mAyANzzj3tVar613lXjut0zkNxvudas+cVNsfjXoy3x5cHke06vat0zt5y872fvRxRFaMoMKCgtQ6GmDAWlOhRqylCk0aNQo0OhRo8iTRmKNWUo1OhRrC1DkbZ8rEhThhJd+c8lWj2KtOV7arR6A4DynhX9PW5dCnToloeH21T8O14bLA4jy5Ytw7Rp0zB58mQAwOrVq7F7926sXbsWb7zxRoX5P/30UwwePBivvvoqAODdd99FREQEPvvsM6xevbqa5RMRWZ/eIGLRrrMV/ogBqHTsbn9dya90/H6fvTWtsj/qd1pzKBmvDOxQ4fBHVWsWAQgo32MyIMCrwiEbc5ZrzZrXHKoYRCqr+24PWvb9vvODCIJg3BvR1LFmnudTpjegRKdHiU6PUq0BxboyFGv1KNXpodH9PU2rR2mZHqU6A0p1f08rM/1ZozNAqy8f05YZysfKyqdp/36VlulRqi2DWiF/cGFWYlEY0Wq1OHHiBObNm2cck8lk6N+/P6Kjoyv9THR0NMLCwkzGBg0ahJ07d95zPRqNBhqNxvg+Ly8PAJCTkwOdTmdJyfel0+lQXFyMGzdumOzGo4rYK8uwX+arD706cfkmsm/ckPwkO4VBRHGxAQqdDHrD7T+YayJO4bkezU3mrW7N2TeKEBmXhJAWTWp0ubdUp+b7nXVyq+4uPg4m25U5y77Xd5aSHIA9AHslACVQHpvkf79qhk6nw/79++Fto8ONGzce/AELFBSUX6Ytig+I7aIFrl69KgIQjxw5YjL+6quvij169Kj0M0qlUtyyZYvJ2MqVK0UPD497ric8PFxEeVjliy+++OKLL77q+SstLe2++ULqoF+pefPmmexNMRgMyMnJgZubW42eVZyfnw8/Pz+kpaXBycmpxpbbELFXlmG/zMdemY+9Mh97ZT5r9koURRQUFMDHx+e+81kURtzd3SGXy5GZmWkynpmZCS+visfoAMDLy8ui+QFArVZDrTY97ubi4mJJqRZxcnLixmom9soy7Jf52CvzsVfmY6/MZ61eOTs7P3Aei241p1KpEBISgsjISOOYwWBAZGQkQkNDK/1MaGioyfwAEBERcc/5iYiIqHGx+DBNWFgYJk6ciO7du6NHjx5Yvnw5ioqKjFfXTJgwAb6+vli8eDEAYNasWejTpw8+/vhjDB06FNu2bcPx48fx5Zdf1uw3ISIionrJ4jAyevRoZGdnY8GCBcjIyEBwcDD27t0LT09PAEBqaipksts7XHr16oUtW7bg7bffxptvvom2bdti586ddeIeI2q1GuHh4RUOCVFF7JVl2C/zsVfmY6/Mx16Zry70ShDFB11vQ0RERGQ9fDwhERERSYphhIiIiCTFMEJERESSYhghIiIiSTWKMLJq1Sp06dLFeEOX0NBQ/PLLL8bppaWlmDFjBtzc3ODg4IBnnnmmwo3aGqMlS5ZAEATMnj3bOMZe3bZw4UIIgmDy6tChg3E6e2Xq6tWrGDduHNzc3GBra4vAwEAcP37cOF0URSxYsADe3t6wtbVF//79cfHiRQkrloa/v3+F7UoQBMyYMQMAt6s76fV6zJ8/Hy1btoStrS1at26Nd9991+Q5KNyubisoKMDs2bPRokUL2NraolevXjh27JhxuqS9etDzaBqCn376Sdy9e7eYkJAgXrhwQXzzzTdFpVIpxsfHi6IoitOnTxf9/PzEyMhI8fjx4+LDDz8s9urVS+KqpRUTEyP6+/uLXbp0EWfNmmUcZ69uCw8PFzt16iSmp6cbX9nZ2cbp7NVtOTk5YosWLcRJkyaJR48eFS9duiTu27dPTExMNM6zZMkS0dnZWdy5c6d46tQp8cknnxRbtmwplpSUSFh57cvKyjLZpiIiIkQA4v79+0VR5HZ1p/fee090c3MTf/75ZzE5OVn89ttvRQcHB/HTTz81zsPt6rZRo0aJAQEB4oEDB8SLFy+K4eHhopOTk3jlyhVRFKXtVaMII5Vp0qSJ+NVXX4m5ubmiUqkUv/32W+O0c+fOiQDE6OhoCSuUTkFBgdi2bVsxIiJC7NOnjzGMsFemwsPDxaCgoEqnsVemXn/9dfEf//jHPacbDAbRy8tL/Oijj4xjubm5olqtFrdu3VobJdZZs2bNElu3bi0aDAZuV3cZOnSoOGXKFJOxp59+Whw7dqwoityu7lRcXCzK5XLx559/Nhnv1q2b+NZbb0neq0ZxmOZOer0e27ZtQ1FREUJDQ3HixAnodDr079/fOE+HDh3QvHlzREdHS1ipdGbMmIGhQ4ea9AQAe1WJixcvwsfHB61atcLYsWORmpoKgL26208//YTu3btj5MiR8PDwQNeuXbFmzRrj9OTkZGRkZJj0y9nZGT179myU/bpFq9Vi06ZNmDJlCgRB4HZ1l169eiEyMhIJCQkAgFOnTuGPP/7AE088AYDb1Z3Kysqg1+thY2NjMm5ra4s//vhD8l7Vyaf2WsPp06cRGhqK0tJSODg44IcffkBAQADi4uKgUqkqPIjP09MTGRkZ0hQroW3btiE2NtbkOOItGRkZ7NUdevbsifXr16N9+/ZIT0/HokWL0Lt3b8THx7NXd7l06RJWrVqFsLAwvPnmmzh27BhefvllqFQqTJw40diTW3dyvqWx9uuWnTt3Ijc3F5MmTQLA/w3e7Y033kB+fj46dOgAuVwOvV6P9957D2PHjgUAbld3cHR0RGhoKN5991107NgRnp6e2Lp1K6Kjo9GmTRvJe9Vowkj79u0RFxeHvLw8fPfdd5g4cSIOHDggdVl1SlpaGmbNmoWIiIgK6ZkquvX/vgCgS5cu6NmzJ1q0aIFvvvkGtra2ElZW9xgMBnTv3h3vv/8+AKBr166Ij4/H6tWrMXHiRImrq7u+/vprPPHEEw98/Hpj9c0332Dz5s3YsmULOnXqhLi4OMyePRs+Pj7criqxceNGTJkyBb6+vpDL5ejWrRvGjBmDEydOSF1a47iaBih/4nCbNm0QEhKCxYsXIygoCJ9++im8vLyg1WqRm5trMn9mZia8vLykKVYiJ06cQFZWFrp16waFQgGFQoEDBw5gxYoVUCgU8PT0ZK/uw8XFBe3atUNiYiK3q7t4e3sjICDAZKxjx47Gw1q3enL3VSGNtV8AcPnyZfz222944YUXjGPcrky9+uqreOONN/Dcc88hMDAQ48ePx5w5c4wPauV2Zap169Y4cOAACgsLkZaWhpiYGOh0OrRq1UryXjWaMHI3g8EAjUaDkJAQKJVKREZGGqdduHABqampCA0NlbDC2tevXz+cPn0acXFxxlf37t0xduxY48/s1b0VFhYiKSkJ3t7e3K7u8sgjj+DChQsmYwkJCWjRogUAoGXLlvDy8jLpV35+Po4ePdoo+wUA69atg4eHB4YOHWoc43Zlqri42OTBrAAgl8thMBgAcLu6F3t7e3h7e+PmzZvYt28fhg8fLn2vrH6KbB3wxhtviAcOHBCTk5PFv/76S3zjjTdEQRDEX3/9VRTF8kvlmjdvLv7+++/i8ePHxdDQUDE0NFTiquuGO6+mEUX26k6vvPKKGBUVJSYnJ4uHDx8W+/fvL7q7u4tZWVmiKLJXd4qJiREVCoX43nvviRcvXhQ3b94s2tnZiZs2bTLOs2TJEtHFxUX88ccfxb/++kscPnx4o70EU6/Xi82bNxdff/31CtO4Xd02ceJE0dfX13hp744dO0R3d3fxtddeM87D7eq2vXv3ir/88ot46dIl8ddffxWDgoLEnj17ilqtVhRFaXvVKMLIlClTxBYtWogqlUps2rSp2K9fP2MQEUVRLCkpEV966SWxSZMmop2dnfjUU0+J6enpElZcd9wdRtir20aPHi16e3uLKpVK9PX1FUePHm1y3wz2ytSuXbvEzp07i2q1WuzQoYP45Zdfmkw3GAzi/PnzRU9PT1GtVov9+vUTL1y4IFG10tq3b58IoNLvz+3qtvz8fHHWrFli8+bNRRsbG7FVq1biW2+9JWo0GuM83K5u2759u9iqVStRpVKJXl5e4owZM8Tc3FzjdCl7JYjiHbeqIyIiIqpljfacESIiIqobGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJMYwQNUCCINz3tXDhQqlLrHH+/v5Yvny51GUQURUopC6AiGpeenq68eft27djwYIFJg+qc3BwkKIsi4miCL1eD4Wi9v6p0mq1UKlUtbY+IuKeEaIGycvLy/hydnaGIAgmY9u2bUPHjh1hY2ODDh064PPPPzd+NiUlBYIg4JtvvkHv3r1ha2uLhx56CAkJCTh27Bi6d+8OBwcHPPHEE8jOzjZ+btKkSRgxYgQWLVqEpk2bwsnJCdOnT4dWqzXOYzAYsHjxYrRs2RK2trYICgrCd999Z5weFRUFQRDwyy+/ICQkBGq1Gn/88QeSkpIwfPhweHp6wsHBAQ899BB+++034+cee+wxXL58GXPmzDHu/QGAhQsXIjg42KQ3y5cvh7+/f4W633vvPfj4+KB9+/YAgLS0NIwaNQouLi5wdXXF8OHDkZKSUhO/HiK6C8MIUSOzefNmLFiwAO+99x7OnTuH999/H/Pnz8eGDRtM5gsPD8fbb7+N2NhYKBQKPP/883jttdfw6aef4tChQ0hMTMSCBQtMPhMZGYlz584hKioKW7duxY4dO7Bo0SLj9MWLF+N///sfVq9ejTNnzmDOnDkYN24cDhw4YLKcN954A0uWLMG5c+fQpUsXFBYWYsiQIYiMjMTJkycxePBgDBs2DKmpqQCAHTt2oFmzZnjnnXeQnp5usmfIHJGRkbhw4QIiIiLw888/Q6fTYdCgQXB0dMShQ4dw+PBhODg4YPDgwSbhiohqSK08jo+IJLNu3TrR2dnZ+L5169bili1bTOZ59913jY+hT05OFgGIX331lXH61q1bRQBiZGSkcWzx4sVi+/btje8nTpwourq6ikVFRcaxVatWiQ4ODqJerxdLS0tFOzs78ciRIybrnjp1qjhmzBhRFEVx//79IgBx586dD/xenTp1Ev/73/8a37do0UL85JNPTOYJDw8Xg4KCTMY++eQTsUWLFiZ1e3p6mjzpdePGjWL79u1Fg8FgHNNoNKKtra24b9++B9ZGRJbhOSNEjUhRURGSkpIwdepUTJs2zTheVlYGZ2dnk3m7dOli/NnT0xMAEBgYaDKWlZVl8pmgoCDY2dkZ34eGhqKwsBBpaWkoLCxEcXExBgwYYPIZrVaLrl27mox1797d5H1hYSEWLlyI3bt3Iz09HWVlZSgpKTHuGamuwMBAk/NETp06hcTERDg6OprMV1paiqSkpBpZJxHdxjBC1IgUFhYCANasWYOePXuaTJPL5SbvlUql8edb52DcPWYwGCxe9+7du+Hr62syTa1Wm7y3t7c3eT937lxERERg6dKlaNOmDWxtbfHss88+8JCJTCaDKIomYzqdrsJ8d6+vsLAQISEh2Lx5c4V5mzZtet91EpHlGEaIGhFPT0/4+Pjg0qVLGDt2bI0v/9SpUygpKYGtrS0A4M8//4SDgwP8/Pzg6uoKtVqN1NRU9OnTx6LlHj58GJMmTcJTTz0FoDws3H0yqUqlgl6vNxlr2rQpMjIyIIqiMVDFxcU9cH3dunXD9u3b4eHhAScnJ4tqJSLL8QRWokZm0aJFWLx4MVasWIGEhAScPn0a69atw7Jly6q9bK1Wi6lTp+Ls2bPYs2cPwsPDMXPmTMhkMjg6OmLu3LmYM2cONmzYgKSkJMTGxuK///1vhZNn79a2bVvs2LEDcXFxOHXqFJ5//vkKe2X8/f1x8OBBXL16FdevXwdQfpVNdnY2PvzwQyQlJWHlypX45ZdfHvg9xo4dC3d3dwwfPhyHDh1CcnIyoqKi8PLLL+PKlStVbxARVYphhKiReeGFF/DVV19h3bp1CAwMRJ8+fbB+/Xq0bNmy2svu168f2rZti0cffRSjR4/Gk08+aXKDtXfffRfz58/H4sWL0bFjRwwePBi7d+9+4LqXLVuGJk2aoFevXhg2bBgGDRqEbt26mczzzjvvICUlBa1btzYeSunYsSM+//xzrFy5EkFBQYiJicHcuXMf+D3s7Oxw8OBBNG/eHE8//TQ6duyIqVOnorS0lHtKiKxAEO8+oEpEVAWTJk1Cbm4udu7cKXUpRFTPcM8IERERSYphhIiIiCTFwzREREQkKe4ZISIiIkkxjBAREZGkGEaIiIhIUgwjREREJCmGESIiIpIUwwgRERFJimGEiIiIJMUwQkRERJJiGCEiIiJJ/T+OBsL6OWeIVwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "%matplotlib inline\n", + "\n", + "# Data into DataFrame\n", + "data_pred = pd.DataFrame({'Temperature': new_temperature_data, 'Intercept': 1})\n", + "data_pred['Frequency'] = predictions\n", + "\n", + "\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 + }, + "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": 82, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/simondelarue/Documents/PhD/Research/Envs/bgdia703/lib/python3.9/site-packages/statsmodels/genmod/generalized_linear_model.py:1257: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified\n", + " warnings.warn(msg, category=PerfectSeparationWarning)\n", + "/Users/simondelarue/Documents/PhD/Research/Envs/bgdia703/lib/python3.9/site-packages/statsmodels/genmod/generalized_linear_model.py:1257: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified\n", + " warnings.warn(msg, category=PerfectSeparationWarning)\n", + "/Users/simondelarue/Documents/PhD/Research/Envs/bgdia703/lib/python3.9/site-packages/statsmodels/genmod/generalized_linear_model.py:1257: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified\n", + " warnings.warn(msg, category=PerfectSeparationWarning)\n", + "/Users/simondelarue/Documents/PhD/Research/Envs/bgdia703/lib/python3.9/site-packages/statsmodels/genmod/generalized_linear_model.py:1257: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified\n", + " warnings.warn(msg, category=PerfectSeparationWarning)\n", + "/Users/simondelarue/Documents/PhD/Research/Envs/bgdia703/lib/python3.9/site-packages/statsmodels/genmod/generalized_linear_model.py:1257: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified\n", + " warnings.warn(msg, category=PerfectSeparationWarning)\n", + "/Users/simondelarue/Documents/PhD/Research/Envs/bgdia703/lib/python3.9/site-packages/statsmodels/genmod/generalized_linear_model.py:1257: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified\n", + " warnings.warn(msg, category=PerfectSeparationWarning)\n", + "/Users/simondelarue/Documents/PhD/Research/Envs/bgdia703/lib/python3.9/site-packages/statsmodels/genmod/generalized_linear_model.py:1257: PerfectSeparationWarning: Perfect separation or prediction detected, parameter may not be identified\n", + " warnings.warn(msg, category=PerfectSeparationWarning)\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkMAAAG6CAYAAAALTELXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVjElEQVR4nO3dd3wUZf4H8M/MbEtvEBI6ghACIYAkgoAUFTlAz15OPCyInnqcIoKeDc/zrIgVOTns8vNOQSwgYq90NIh0pAXSe7J95vn9scmSJQGSzSSbZD7v1yuSzMzuPPt1k3zyPM88IwkhBIiIiIgMSg51A4iIiIhCiWGIiIiIDI1hiIiIiAyNYYiIiIgMjWGIiIiIDI1hiIiIiAyNYYiIiIgMjWGIiIiIDI1hiIiIiAytVYWhf//737j22mtPekxJSQnuuusuZGRkIDMzEw8//DAcDkcLtZCIiIjaG1OoG1DjnXfewbPPPothw4ad9LiZM2fC4XDg9ddfR3l5Oe677z7Y7XY88cQTLdRSIiIiak9CHoby8vLw0EMPYf369ejZs+dJj/3555+xYcMGrFq1Cr179wYA/OMf/8D06dMxa9YsdOrUqQVaTERERO1JyIfJfvvtN5jNZnz00UdIT08/6bGbNm1Cx44d/UEIADIzMyFJEjZv3tzcTSUiIqJ2KOQ9Q+PHj8f48eMbdGxeXh6Sk5MDtlksFsTGxiInJ6c5mkdERETtXMh7hhrD4XDAYrHU2W61WuFyuYJ+XiFEU5pFREREbVjIe4Yaw2azwe1219nucrkQHh4e9PNKkoTycgdUVWtK8wxPUWRER4exlk3EOuqHtdQPa6kP1lE/MTFhkGV9+nTaVBhKSkrCF198EbDN7XajtLQUiYmJTXpuVdXg9fKNqQfWUh+so35YS/2wlvpgHZtOz0GdNjVMlpGRgdzcXBw8eNC/bcOGDQCAM844I1TNIiIiojasVYchVVVRUFAAp9MJAEhPT8fQoUNx5513YuvWrVi3bh0efPBBXHTRRbysnoiIiILSqsNQTk4ORo0ahVWrVgHwze158cUX0bVrV0ybNg133HEHzj77bMybNy+0DSUiIqI2SxK8lAoAUFJSxfHbJjKZZMTFRbCWTcQ66oe11A9rqQ/WUT/x8RFQFH36dFp1zxARERFRc2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkNjGCIiIiJDYxgiIiIiQ2MYIiIiIkMLeRjSNA3PP/88Ro8ejcGDB+Omm27C4cOHT3h8UVER7rrrLgwfPhxnnnkm7rzzTuTl5bVgi4mIiKg9CXkYWrhwIZYuXYpHHnkE7777LjRNw/Tp0+F2u+s9/o477sDRo0fx2muv4bXXXsPRo0dx2223tXCriYiIqL0IaRhyu9149dVXMXPmTIwdOxYpKSlYsGABcnNzsWbNmjrHl5eXY8OGDbjpppvQv39/pKamYsaMGfj1119RWlra8i+AiIiI2jxTKE++c+dOVFVVYcSIEf5t0dHRSE1NxcaNGzFlypSA4202GyIiIrBixQpkZmYCAD788EP06tUL0dHRTWqLooS8k6zNq6kha9k0rKN+WEv9sJb6YB31I0n6PVdIw1Bubi4AIDk5OWB7YmKif19tFosFjz/+OB588EEMGzYMkiQhMTERb7/9NmS5aW+s6OiwJj2ejmEt9cE66oe11A9rqQ/WsXUJaRhyOBwAfCGnNqvVirKysjrHCyGwY8cODBkyBNOnT4eqqliwYAFuvfVW/N///R8iIyODbkt5uQOqqgX9ePL9pRMdHcZaNhHrqB/WUj+spT5YR/3ExIQ1uSOkRkjDkM1mA+CbO1TzOQC4XC6EhdVNzZ9++inefvttfP311/7gs2jRIowbNw7vv/8+rrvuuqDboqoavF6+MfXAWuqDddQPa6kf1lIfrGPTCaHfc4V00LJmeCw/Pz9ge35+Pjp16lTn+E2bNqFXr14BPUAxMTHo1asXDh482LyNJSIionYppGEoJSUFkZGRWL9+vX9beXk5tm/fjoyMjDrHJyUl4eDBg3C5XP5tdrsd2dnZ6NmzZ0s0mYiIiNqZkIYhi8WCqVOn4umnn8aXX36JnTt34s4770RSUhImTJgAVVVRUFAAp9MJALjooosA+NYa2rlzJ3bu3IlZs2bBarXikksuCeErISIiorYq5Nf2zZw5E5dddhnuv/9+XH311VAUBUuWLIHZbEZOTg5GjRqFVatWAfBdZbZ06VIIITBt2jRcf/31MJvNWLp0KaKiokL8SoiIiKgtkoTQcwpS21VSUsXJbE1kMsmIi4tgLZuIddQPa6kf1lIfrKN+4uMjdFuvKeQ9Q0REREShxDBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhsYwRERERIbGMERERESGxjBEREREhhbyMKRpGp5//nmMHj0agwcPxk033YTDhw+f8HiPx4P58+f7j586dSp27NjRgi0mIiKi9iTkYWjhwoVYunQpHnnkEbz77rvQNA3Tp0+H2+2u9/h58+Zh+fLl+Ne//oVly5YhPj4eN910EyoqKlq45URERNQehDQMud1uvPrqq5g5cybGjh2LlJQULFiwALm5uVizZk2d4w8fPoxly5bh0UcfxejRo9G7d2/885//hMViwbZt20LwCoiIiKitM4Xy5Dt37kRVVRVGjBjh3xYdHY3U1FRs3LgRU6ZMCTj+xx9/RFRUFM4+++yA47/66qsmt0VRQt5J1ubV1JC1bBrWUT+spX5YS32wjvqRJP2eK6RhKDc3FwCQnJwcsD0xMdG/r7b9+/ejW7duWLNmDV555RXk5eUhNTUV99xzD3r37t2ktkRHhzXp8XQMa6kP1lE/rKV+WEt9sI6tS0jDkMPhAABYLJaA7VarFWVlZXWOr6ysxMGDB7Fw4ULMmTMH0dHRePnll/GnP/0Jq1atQkJCQtBtKS93QFW1oB9Pvr90oqPDWMsmYh31w1rqh7XUB+uon5iYMMiyPj1sQYWhTz75BBMmTKgTYhrLZrMB8M0dqvkcAFwuF8LC6qZmk8mEyspKLFiwwN8TtGDBAowZMwYffPABpk+fHnRbVFWD18s3ph5YS32wjvphLfXDWuqDdWw6IfR7rqAi1Zw5czBy5EjMmzcPW7duDfrkNcNj+fn5Advz8/PRqVOnOscnJSXBZDIFDInZbDZ069YN2dnZQbeDiIiIjCuoMPTVV1/hhhtuwLp163DllVdi0qRJWLJkCQoKChr1PCkpKYiMjMT69ev928rLy7F9+3ZkZGTUOT4jIwNerxe//vqrf5vT6cThw4fRo0ePYF4KERERGVxQYSgpKQl/+ctfsHr1arzzzjsYNmwYFi9ejHHjxuGWW27BmjVr4PV6T/k8FosFU6dOxdNPP40vv/wSO3fuxJ133omkpCRMmDABqqqioKAATqcTADBs2DCcddZZmDt3LjZt2oS9e/dizpw5UBQFf/zjH4N5KURERGRwkhD6jLr9+uuvePLJJ7Fx40YAQIcOHTBt2jTccMMNUBTlhI9TVRXPPPMMli9fDqfTiYyMDDz44IPo2rUrsrOzcc455+Cxxx7DJZdcAsA3ifrpp5/G6tWr4XQ6MXToUPz9739Hnz59mtT+kpIqjt82kckkIy4ugrVsItZRP6ylflhLfbCO+omPj9BtiYImhaEjR47gww8/xIcffohDhw6he/fuuPTSSzF27Fh88803WLhwIc4//3w88cQTujS2OfGN2XT8JtcH66gf1lI/rKU+WEf96BmGgrqa7L333sOHH36ILVu2wGq1YuLEiXj00UcxbNgw/zF9+/ZFSUkJ3n333TYRhoiIiMiYggpDDzzwANLT0zFv3jxMmjQJkZGR9R7Xr18/XHnllU1qIBEREVFzCnqdoT59+kBVVf98IKfTCY/Hg6ioKP9xF110kS6NJCIiImouQQ229ezZEw899BCuuOIK/7YtW7ZgxIgReOKJJ6BpHAclIiKitiGoMPT888/jo48+CriRampqKmbPno3//e9/+M9//qNbA4mIiIiaU1DDZB9//DHmzp2Lq666yr8tNjYW1113HUwmE958803MmDFDt0YSERERNZegeoZKSkrQrVu3eveddtpp9d5xnoiIiKg1CioMnXbaafjss8/q3ffVV1/x1hhERETUZgQ1TPbnP/8Z99xzD0pLS3HuueciISEBxcXF+Prrr/Hpp5/iscce07udRERERM0iqDB00UUXoaqqCgsXLsSaNWv82+Pi4vDAAw/wknoiIiJqM4IKQwBwzTXX4E9/+hP279+P0tJSREdH47TTToMs67M0NhEREVFLCDoMAYAkSTjttNP0agsRERFRiwsqDBUXF+PRRx/FN998A4fDgePv9SpJErZv365LA4mIiIiaU1Bh6B//+Ae+/vprTJ48GUlJSRwaIyIiojYrqDD03Xff4e9//ztvwkpERERtXlBdOmaz+YSLLhIRERG1JUGFofPOOw+ffPKJ3m0hIiIianFBDZOlpqbi2WefxeHDh5Geng6bzRawX5Ik3Hbbbbo0kIiIiKg5SeL4S8EaICUl5eRPKknYsWNH0I0KhZKSKni9Wqib0aaZTDLi4iJYyyZiHfXDWuqHtdQH66if+PgIKIo+F3AF1TO0c+dOXU5OREREFGpNjlQVFRXYt28f3G43VFXVo01ERERELSboMLR+/XpcfvnlyMzMxAUXXIA9e/bgrrvuwuOPP65n+4iIiIiaVVBhaO3atbjxxhths9kwe/Zs/wrUKSkpePPNN/Haa6/p2kgiIiKi5hJUGHr22Wdxzjnn4K233sK0adP8YeiWW27B9OnT8d577+naSCIiIqLmElQY2rFjBy699FIAvivHahs5ciSOHDnS9JYRERERtYCgwlBUVBQKCgrq3ZeTk4OoqKgmNYqIiIiopQQVhs455xwsWLAAv/76q3+bJEnIzc3FokWLMHbsWL3aR0RERNSsglpn6K677kJWVhauuOIKdOjQAQAwa9Ys5ObmIjk5GbNmzdK1kURERETNJagwFBMTg/feew8rVqzAunXrUFpaiqioKFx77bW45JJLEBYWpnc7iYiIiJpFUGEIACwWC6644gpcccUVeraHiIiIqEUFFYZWrFhxymMuuuiiYJ6aiIiIqEUFFYbuueeeerdLkgRFUaAoCsMQERERtQlBhaEvv/yyzja73Y5NmzZh8eLFeOmll5rcMCIiIqKWEFQY6tKlS73bTz/9dHg8HjzyyCNYunRpkxpGRERE1BKafNf64/Xr1w+//fab3k9LRERE1Cx0DUNutxvvv/8+EhIS9HxaIiIiomYT1DDZ+PHj69yTTNM0lJSUwOVyYe7cubo0joiIiKi5BRWGMjMz64QhAIiMjMS4ceNw1llnNblhRERERC0hqDD0+OOP690OIiIiopAIKgwdPXq0Ucd37tw5mNMQERERNTvd5gydzI4dO4I5DREREVGzCyoMPfvss3jooYcwYMAAXHjhhejUqRNKSkrw1Vdf4dNPP8Vf/vKXE65FRERERNSaBBWGPvzwQ4wbN67O3KFJkyYhISEBW7Zswe23365LA4mIiIiaU1DrDK1duxZTpkypd9/ZZ5+NzZs3N6lRRERERC0lqDAUFxeHrKysevetXbsWnTp1alKjiIiIiFpKUMNkl112GV5++WU4HA6MHz8e8fHxKCwsxOrVq/F///d/eOCBB/RuJxEREVGzCCoM3XrrraioqMDrr7+OJUuWAACEEAgLC8Odd96Jq666StdGEhERETWXoMKQJEm45557cOutt+KXX35BWVkZ4uLiMHjwYERGRurdRiJqBrIsQdNEqJtBRBRyQYWhGpGRkUhMTAQADB48GF6vV5dGEVHz86oa5EasF0ZE1F4FHYY+/PBDzJ8/HwUFBZAkCe+99x5eeOEFmM1mzJ8/HxaLRc92EpHOhADcqgqLSQl1U4iIQiqoq8lWrVqFuXPnYvjw4XjmmWegaRoA4LzzzsO3336LhQsX6tpIImoeTrfaqNXkiYjao6B6hhYtWoSrrroK8+bNg6qq/u2XXnopiouL8b///Q933HGHXm0kombiVQVUjcNlRGRsQfUM7d+/H+edd169+9LT05GXl9ekRhFRy9A0Ud07FOqWEBGFTlBhKCEhAfv27at33759+5CQkNCkRhFRy9CEgNPtBa8pIyIjCyoMTZo0Cc8//zxWr14Nt9sNwHe5/bZt27Bw4UJMnDhR10YSUfPxqgIerxbqZhARhUxQc4buuOMO7N69G3fccQdk2Zenrr32WtjtdgwbNgx/+9vfdG0kETUfTRNwurywRFghBPuIiMh4ggpDFosF//nPf/Djjz9i3bp1KC0tRVRUFDIzMzFmzBhenULUxjg9KiKFBgn83iUi4wkqDN14442YPn06Ro4ciZEjR+rdJiJqYTUTqcMsTVqHlYioTQpqztCWLVvY+0PUjggBOF1eXlVGRIYUVBgaPXo0PvroI3g8Hr3bQ0Qh4lEF3JxITUQGFFSfuNVqxUcffYRPP/0UvXv3Rnh4eMB+SZLwxhtv6NJAImoZxyZSW8B51ERkJEGFodzcXAwZMsT/9fFXoPCKFKK2yeXRoAnBidREZCgNDkNr1qzB8OHDER0djbfeeku3BmiahhdffBHvvfceKioqkJGRgQcffBDdunU75WM/+ugj3H333fjyyy/RtWtX3dpEZFSqpnEiNREZToPnDP3tb3/DgQMHArYtXrwYRUVFTWrAwoULsXTpUjzyyCN49913oWkapk+f7l/M8USOHDmCf/zjH006NxEF4kRqIjKiBoeh44e+VFXFM888g9zc3KBP7na78eqrr2LmzJkYO3YsUlJSsGDBAuTm5mLNmjUnfJymabj77rsxYMCAoM9NRPXjRGoiMpom9YU3dW7Qzp07UVVVhREjRvi3RUdHIzU1FRs3bsSUKVPqfdyiRYvg8Xhw++23Y926dU1qQw1FCerCOqqlpoasZdO0VB0lTYOiSPX2Arm9GsJtpjY/kZrvSf2wlvpgHfWjZw92SCcG1PQqJScnB2xPTEw8YY/T1q1b8eqrr+L9999HXl6ebm2Jjg7T7bmMjrXUR3PX0eVR4RZSvYFHliRERNpgNinN2oaWwvekflhLfbCOrUtIw5DD4QDgu71HbVarFWVlZXWOt9vtmD17NmbPno2ePXvqGobKyx1QVQ4NNIWiyIiODmMtm6il6qhqGsrLndC0umlIAqB5vQi3tu2J1HxP6oe11AfrqJ+YmDD//VGbqsk/6ZqyErXNZgPgmztU8zkAuFwuhIXVTc3//Oc/0atXL1x11VVBn/NEVFWDl/MkdMFa6qO566hqAqoq6g1DAFBld8Nqktv8UBnA96SeWEt9sI5Np+fPpkaFodtuu61OL84tt9wCs9kcsE2SJHzxxRenfL6a4bH8/Hx0797dvz0/Px/9+vWrc/yyZctgsVj8axypqgoAmDJlCm655RbccsstjXk5RHQSNROpzZzbQETtXIPD0MUXX6z7yVNSUhAZGYn169f7w1B5eTm2b9+OqVOn1jn++CvMsrKycPfdd+OVV15B3759dW8fkZFxRWoiMooGh6HHHntM95NbLBZMnToVTz/9NOLj49GlSxc89dRTSEpKwoQJE6CqKoqLixEVFQWbzYYePXoEPL5mknXnzp0RGxure/uIjI4rUhOREYS8/3vmzJm47LLLcP/99+Pqq6+GoihYsmQJzGYzcnJyMGrUKKxatSrUzSQypJoVqYmI2jNJ8EZiAICSkipOZmsik0lGXFwEa9lELVVHVRMoOsHVZLVZzTLio21tcqiM70n9sJb6YB31Ex8fodt6TSHvGSKi1o0rUhNRe8cwREQnpWkCDpe3SctoEBG1ZgxDRHRKLo8K9RTDaUREbRXDEBGdkqYKON28mz0RtU8MQ0R0SgKAw+0F+4aIqD1iGCKiBlG9Am4PL7MnovaHYYiIGkQTAg4nJ1ITUfvDMEREDeb2avDwTttE1M4wDBFRg6nV9ytj5xARtScMQ0TUKC63Ci5cT0TtCcMQETWKl/crI6J2hmGIiBpFCFSvSB3qlhAR6cMU6gYQUcvShMChvAo43SpcbhUd48IgNzLZ1NyvzKzTTRKJmqrmfV1p9yAy3IzunaIa/b4m42IYIjKQHQeKsXLdQeQW25EQbUOlw41wmxlj0zujd9fYBj+Ppvkus7dEWjl/iEKu9vtaVQUURUJSfDgmD++B/j3jQ908agP4Zx2RQew4UIw3PtuF7IJKWM0KIsLNsJgV5BY78MEP+7Evu7RRz+fyqlA1XmZPoXX8+zo60gKrWUF2QRXe+GwXdhwoDnUTqQ1gGCIyAE0IrFx3EE63F7GRVljMCmRJgtmkIDrcDJdHwzdZR6E1opdHVQXnDlFIneh9bTEriI20wOlWsXLdwUa9r8mYGIaIDOBQXgVyi+2IsJnrrCAtSRLCrQoKSx3IKaxq1PM6XbzMnkLnVO/rCJsJucV2HMqrCFELqa1gGCIygEq7B6oqYDLV/y2vKDJUDbA7vY16Xl5mT6F0qve1ySRDVQUq7Z4Wbhm1NQxDRAYQGW6Gokjweuuf46OqGhQZCLc17poKXmZPoXSq97XXq0FRJESGm1u4ZdTWMAwRGUD3TlFIig9HldNbZ1hLCAG7S0WH2DAkd4ho9HN7vAIu3s2eQuBU7+sqpxdJ8eHo3ikqRC2ktoJhiMgAZEnC5OE9YLMoKK10w+1RoQkBj1dFud0Dq1nG2PTOQa3LwrvZU6ic6H3t9qgorXTDZlEweXgPrjdEp8QwRGQQ/XvGY9r5/dC1YwRcHhVVdg/cHhVJ8WG4eFSvRq0zdDy3VzvhUAVRczr+fV1e6YbLo6JrxwhMO78f1xmiBuGii0QG0r9nPPr1iGvyCtTHUzUBh9uLqHAzeHEZtbTa72uuQE3BYBgiMhhZktAzKRqqJlBU7oSm6ZNenG4VEWEmSOAvIGp5Ne9romBwmIyIdKFqGhwulVeWEVGbwzBERLqoucyew2RE1NYwDBGRbryqBicvsyeiNoZhiIh0IwTgcHo4VEZEbQrDEBHpyqMKuDy8zJ6I2g6GISLSlaaJ6t4hdg8RUdvAMEREunN7NXi8nDtERG0DwxAR6U7VRPUNXNk7REStH8MQETULp0eFV+XcISJq/RiGiKhZqGpN71CoW0JEdHIMQ0TUbJxuFTrd7YOIqNkwDBFRs1FVDXauO0RErRzDEBE1GwHA4VKh8R4dRNSKMQwRUbPy9Q5x7hARtV4MQ0TUrNg7REStHcMQETU79g4RUWvGMEREzY69Q0TUmjEMEVGLUDUNDpfK3iEianVMoW4AEbUsIQRWrj2IA7kVOGtgEjp3iGih8wIOpxfhVgUAExERtR4MQ0QGcyC3Asu/+x0AsO33Ikz7Qwp6JUe3yLm9mgaHW0WYhT96iKj14DAZkcHER9tgNSsAfHeXf/3Tndh3tKxFzi0EYHd5W+RcREQNxTBEZDAxERbcfmkazCbft7/Hq+HNT3dhb3bLBCKvV4PTzUBERK0HwxCRAQ3oGY+/XjroWCBSNbz52U7sOlTS7Odm7xARtTYMQ0QG1b9HHK77Qwos1YHIqwq8vWY3th8obvZze9g7REStCMMQkYH17hKD6yf1988hUjWBpZ/vRtbewmY9rxBAFXuHiKiVYBgiMrgeSVG4YXIKbBZfINIE8L+v9mLTzvxmPa/Xq8Hu4qrURBR6DENEhG6JUZg+JRXhNt8l7wLA8u9+xw9bc5rtnEIAVQ4PXB4uxEhEocUwREQAgM4dInDTBamICjf7t61adxBfbDoM0Uy30VA1gYoqD7wqb9NBRKHDMEREfp3iwnHzhQMQF2X1b/tqyxF8/NOBZruvmEfVUF7lapbnJiJqCIYhIgoQH23DjAsHoGNsmH/but/y8P7X+6BqWrOc0+3V4HBz/hARhQbDEBHVERNhwYwLU9G147H7lv2ytxBvfbYbbo+q+/mEAOxOL+9qT0QhwTBERPWKsJlx4+RU9O5y7L5luw+XYsnKHahyenQ/n9erwe5k7xARtTyGISI6IatFwbSJKRh4Wrx/2+H8Svz7w99QUuHU9VwCgMOlcjI1EbU4hiEiOimTIuOq8afjzNRO/m2FZU4sWvEbjhZW6Xour6qhwu4G2DtERC2IYYiITkmWJVw4sifOHdbVv63C4cErH/+GPdmlup7L5VFRafdA4ngZEbUQhiEiahBJkjB+aFdccvZpkKtzituj4Y1Pd2HzLv1Wq665kavTw/lDRNQyGIaIqFGGpSTi2vP7+W/wqgmBZd/+js836rc4o6YJlFe5YXd5IctMRETUvBiGiKjR+nWPw00XpCIy7Nhq1V//fAT/+3ovvKo+axGpqm916rIqziEiouYV8jCkaRqef/55jB49GoMHD8ZNN92Ew4cPn/D4PXv2YMaMGTjzzDMxYsQIzJw5E0ePHm3BFhMRAHTpGIm/XBS4OGPW3iIs+WQHKh36XHqvCQG704OySjd815sREekv5GFo4cKFWLp0KR555BG8++670DQN06dPh9vtrnNsSUkJrr/+ethsNrz11ltYvHgxiouLMX36dLhcXM6fqKXFRdlwyx8H4LTOx9YiOphXgZdXbENeiV2XcwgBOFxelFS6IRiIiKgZhDQMud1uvPrqq5g5cybGjh2LlJQULFiwALm5uVizZk2d47/44gvY7XY8+eST6Nu3LwYOHIinnnoK+/btw5YtW0LwCogozGrC9ZNSMKxfR/+2kgoXFq34DbsOleh2HpdbRSkDERE1A1MoT75z505UVVVhxIgR/m3R0dFITU3Fxo0bMWXKlIDjR4wYgYULF8Jms/m3ybIvz5WXlzepLYoS8k6yNq+mhqxl07RUHSVNg6JIulyxpSgKLhvXG4nx4fh07UEI+C6Rf/OzXZg0vAdGpyfrcqm8V9VQ6fAiNsoCqQETifie1A9rqQ/WUT96Xm0a0jCUm5sLAEhOTg7YnpiY6N9XW9euXdG1a9eAba+88gpsNhsyMjKa1Jbo6LBTH0QNwlrqo7nr6PKocAsJet4O7MIxfdA9OQavffwbXB4VQgAr1x5EUYUL10xMgdmk6HIeFTJiIq2wmBv2fHxP6oe11Afr2LqENAw5HA4AgMViCdhutVpRVlZ2yse/9dZbePvtt3H//fcjPj7+lMefTHm5A6pOV8EYlaLIiI4OYy2bqKXqqGoaysud0DR9h516dYrAXy4egDc+3YWSCt9cvnXbcnEkvxLXnt8XMZFWXc5TWuZAdLgFVot8wkDH96R+WEt9sI76iYkJ848ONVVIw1DNcJfb7Q4Y+nK5XAgLO3FqFkLgueeew8svv4y//OUvuPbaa5vcFlXV4PXyjakH1lIfzV1HVRNQVaF7GAKAxNhw/OWigVj6+W4cyK0A4Lun2fPv/4qrzz0dvZKjT/EMp6aqXng8KsKspoBL/Os/lu9JvbCW+mAdm07PXu2QDlrWDI/l5weuXpufn49OnTrV9xB4PB7cfffdWLRoEe69917ccccdzd1MIgpCZJgZN0zuH3BPs0qHB0s+2YGftuXqskCjqglUOT0oqXBCVQVXrCaioIQ0DKWkpCAyMhLr16/3bysvL8f27dtPOAdozpw5WL16NebPn4/rrruuhVpKRMEwKTL+OKoXLhrdC0r1StKaEPjkpwN47+t9cHvUJp9DCMDl0VBc6YTdxVt4EFHjhXSYzGKxYOrUqXj66acRHx+PLl264KmnnkJSUhImTJgAVVVRXFyMqKgo2Gw2LF++HKtWrcKcOXOQmZmJgoIC/3PVHENErU9m/05Iig/H0s93o9zuW5Dxl72FyC2240/nno4OsU2fTFqzYrXXqyEqwsI1GomowUJ+bd/MmTNx2WWX4f7778fVV18NRVGwZMkSmM1m5OTkYNSoUVi1ahUA4JNPPgEAPPnkkxg1alTAR80xRNQ6de8UhdsuSUPP5Cj/ttxiO176YBt+/b1Il3NoQsDu8lavWE1E1DCS0OvOim1cSUkVJ7M1kckkIy4ugrVsopaqo6oJFDXD1WSnPq+GzzYcxg9bcwK2jxiQhD8M7w6TDuuvSBJgs5gQF21FQnwk35M64Pe3PlhH/cTHR+i2XlNIh8mIyHgUWcbEM7sj0mbCl1uOwFP9C2Htb7k4lFeBq849HQnRNmhCIKewCnanF+E2E5I7REBu4ISgmlt4oAKIiLQ12zwiTQgcyqtApd2DyHAzuneKanAb2xKvpmH9r3mwuzSEW2Wc0a8jTDpd0kzUGjAMEVGL2pddim+yjqKw1AGzIkPVNGjVfyAfKazCi8t+xVkDOuFQQSUKSx1QNUCRgQ6xYRib3hm9u8Y2+Fxuj4qSchc0rxdW04nXIwrGjgPFWLnuIHKL7VBVAUWRkBQfjsnDe6B/z6ate9aarF5/ECvXHoTD5YUAIAF422rC5BE9MPHMHqFuHpEuGO2JqMXsyy7FBz/sR26xHRazguhIC+KjrAFDYy6Piq9/OYqDuRUwm2REhpthMSvILXbggx/2Y192aaPOqQmB8iq3rvc123GgGG98tgvZBZWwVr8Oq1lBdkEV3vhsF3YcKNblPKG2ev1BLPv2d1Q5vZBlCSZFgixLqHJ6sezb37F6/cFQN5FIFwxDRNQiNCHwTdZRuDwqosMtMJsUyJIEi9mEjrE2WM1KwHCWVxUorXTDq2owmxREh5vh8mj4JusotEZ28Qgh4HB5UVzugsurNuk+aZoQWLnuIJxuL2Krbwniex0KYiMtcLpVrFx3sNFtbG28moaVaw9C1QTM1SFIlmTIsgSzIkHVBFauPQivxnkv1PYxDBFRi8gprEJhqQPhVlOdMCJJEqIjzLBZlIA5N15VoLDUiQq77+qwcKuCwlIHcgqrgmqDx6uhrNKNsioXVE0EFYoO5VUgt9iOCJu53tcRYTMht9iOQ3kVQbWxtdiwPQ8OlxcmWar3dZpkCQ6XFxu254WohUT6YRgiohZhd3p9839OcPWHosgQGmA2+QJFbRV2D4rKnAAAVfM9V7A0TcDu9KK4wolyu9u/cnVDc1Gl3QNVFTCZ6n8dJpMMVRWorF5Pqa0qLnf6BhVPVBfJt5RTcbmz5RpF1Ew4gZqIWkS4zQRF9t2TSa7n7vWqqkExSQBkhFsV2CwmlFS6/Jf+u70aCkudsFkUhFkbdrf6k1FVgSqHBw6XF2ZFgs1iqn7ek6eiyHAzFEWC16vBYq7bDq9Xg6JIiAw/+f3SWrv4aJuvEjWzpo8nfJvjo7nYLbV97BkiohaR3CECHWLDYHepde5LJoSA3aWiU1wYkuJ9x1jMMhJjw2CzHAscAoDDreLrn4/6h86aStMEXB4N5XY3ispdcLp9vU4n6inq3ikKSfHhqHJ6630dVU4vkuLD0b1TVP1P0EZkpnZCmNUErybqfZ1eTSDMakJmav33kSRqSxiGiKhFyJKEsemdYTXLKLd74PGq0ISAx6ui3O6B1Sxj3OAuAceomoaYSAuijutl2XGwBM+9txVZewt1ueEr4FubyOPVUFblRlG50xd2UPfmr7IkYfLwHrBZFJRWuuH2+F6H26OitNINm0XB5OE92vx6QyZZxuQRPaDIEjyqgKYJaEKDpgl4VAFFljB5RA+uN0TtAlegrsbVQJuOK6vqo72vQF17naETrSFU3zGxkVZ4vBqOFtkDni+1Zxz+OKoXosItdc6lKBJiYsJRVuZbC6gxJPjmMdksSvUQnxSwTpGR1xkK4zpDQePPSf3ouQI1w1A1vjGbjt/k+mjvYQhAg1aXru8YCcD6HXlYve4Q3LVqE2ZVMHlETww5vUPAlU9NCUO1KYqEMIsJEWEmyNKxUGSkFag37yrgCtQ64M9J/TAMNQO+MZuO3+T6MEIYaqricieWffs79ueUB2zv0yUGF43u5Z/Uq1cYqmFWZITbTAizGu/aE35/64N11I+eYYjRnojanPhoG26c0h9/HNUL1lpXdO09Uobn3tuK77KOQm2GxQA9qm+idUmFs3rxRt1PQUQhwDBERG2SLEk4M7UT/nb5IPTrHuvf7lE1rF5/CC8t34aDufovfCgE4PJoKKtwo7jcCXcTV7QmotBjGCKiNi020oo/n98PV53TBxFhx646yy22Y+EH2/DWpztQ5dR/AURN+C7JL61e0VoTDV+4kYhaF+MNfBNRuyNJEgb17oA+XWLx2YZD2Lgz37/vx6yj+HlnPs7L6IaMlETIsr6JpWZFa49HQ5hVgc1qgiLLul3yT0TNjz1DRNRuhNtMuPjs0zDjwlR0igvzb7e7vPjwh/14eUXzDJ0BNfOJPCgqd6K00uW72q0Rt/kgotBhGCKidqdnUjRuvzQNk0f0gLXWCtZHCqvw749+w/++2ouyKn1WsD6eqgo4XF6UVrhQXL14oyaCuyksEbUMDpMRUbukyDLOHtwZo4d2w7trduKXPYX+fb/sLcRvB4pxdnpnjE5PhqWee6U1lW9VagGPx40qpwSzIsNqUWA2KTArvmDEkTSi1oFhiIjatdgoK64+93RkpCTik58OIKd6BWuPV8OXm7OxcWc+JmR0w+A+HXSfTwT47qemqgKqqsLlUSFLEmRZgtWsVIcjGRIYjIhCicNkRGQIvZKjcdvFabhodC9E2I79HVhe5cb73+zDSx/8it2HS5t14rMQvsUuPV4NlQ4PSipcKCpzotLhgUfVOL+IKETYM0REhiHLEjL7d8Kg3gn45ucj+PHXXKjVK3DnFNnx+qc70btLNM7P7I6uHSObvT2a5rsBqserQZa9MFXf9sNmVQJu+0FEzYs9Q0RkODaLCRPP7IFZV6ZjUO+EgH37jpRj4Qfb8M7nu5Ff4mixNmmagNvjW+G6qMyF8io3PGrNFWkSe42ImhF7hojIsOKibLjqnNMxKi0Zqzccwu9Hj93r7Lf9xdh+oBiD+3TA+DO6IqH6fmfNTQjAq2rwqhocbt8cI0XxTcC2mGVYTDIkiesYEemJYYiIDK9rYiRunNwfe7LLsGbDIRytnmQtBPDznkJk7S3EkL4dMW5IF/9NYFuCpgloEPCqgAsqZKdv8rXZ5JuAbTErXOCRSAcMQ0RE8A1F9e0Wiz5dY/Db/mJ8vvEwCsucAABNAJt3FeDn3QUYfHpHjB3SGR1iwk7xjPrThICm+sKR0636glHNJfuKDJMiQ5IkhiOiRmIYIiKqRZYkpJ2WgNSe8cjaW4ivNmejuMIFwBeKtuwuwM97CjCodwLGDO6CpPjwkLRTiFqX7LtVSLIEWZJgNvmG0xiOiBqOYYiIqB6KLGFo345I75OAX/YU4ustR/yhSAgga28RsvYWIaV7HMYM7oweSVEha6sAIPxDahqcbl9PlyxJMCsSLBYFFkWBycTFHonqwzBERHQSiizjjH6JGHx6B2TtLcLXPx9BUfXwGQDsPFSCnYdK0CMpCmcPSka/HnGQQ3zplxCAEMfmGzmre44UWYLF7Fvo0ddzJAG+JR8ZkMjQGIaIiBpAkWUM7dsRg/t0wK+/F+HbX44it9ju338wtwJv5VagQ4wNI9OSMaRvh2a5zUcw/D1H1WsaSRL8K2ErsgSTIsNcfaWazCvVyIAYhoiIGkGWJaT36YBBvROw63ApvvvlKA7kVvj3F5Y58eEP+7Fm42Fk9k/EmamdEBtpDWGL6xICUIXwrYYNAFAhOX2vzaLI/mE1WQGvViNDYBgiIgqCJElI6R6HlO5xOJRXge+35mD7/mLUxAaHy4tvfzmK77OOIrVnPIYPSEKv5KhWe/f6mgnZDlWF06P6F3qUJQkWk4wwmxkutxcCArJcszo2h9eofWAYIiJqou6donDNeVEoKnPix2052LKrAG6vBsB3Bdq2/cXYtr8YneLCcGZqJww+vQNsltb747dmzhEAqPANrbm8KoQio6LCASEAkyxBUWQokuSfj2RWJEiSDIYkamskwf5PAEBJSRW81T+8KDgmk4y4uAjWsolaqo6qJlBU7oSmtd8fAYoiISYmHGVldqhqy71Oh8uLTbvyse63PJRUX4FWm8UkY1CfDshMSUSXjhGttreotpPVUpKOXb2mKL7L+03VE7Rr5iZJ4FVsAH9O6ik+PgKKos9dxVrvnyZERG1UmNWE0YM6Y+TAZOw6VIK1v+Vh75Ey/363V8OmnfnYtDMfyQnhGNYvEel9OiDc1jZ/JB9/9ZrLrfoDkoTqoTbLsavYFFmqDoDsQaLWoW1+5xERtQGyLKF/z3j07xmPwlIH1m/Pw5Y9BXC4VP8xOUV2fPzTAXy6/iD694jHGf06ok+XGMhy6+8tOpk6Q20OLaAHyaT45iIpJhkmSYasoHpJAoYkankMQ0RELaBDbBgmn9UTEzK749ffi7BxZz4O1roKzasK/Pp7EX79vQjR4Wak9+mAwad3QHJCRAhbra961z+q6UGqnqxtkiUoJl/vkW94zbev5vEAai0N4LvaTRMCQhPVAdIXpk50fqL6MAwREbUgs8m3XtHQvh2RX+LA5l352LKnEFUOj/+YcrsH32/Nwfdbc5AUH470PgkY1LsD4qJa1yX6eqjTgwQAbhUS4Ms1AKTqT0R1yJEgAVL1XKVa+yT45izVTOqWpOrnrw5KNR+KLPnDVO3eqJr2kPEwDBERhUhiXBj+MLwHJmR2w65DpdiyuwA7D5ZCq/UbObfYjtwNdny24TC6d4rEoN4JGNgrAdERlhC2vPkJ/3+OhaBj+0Stzh8RsMer4qT8AapWmFJk35BmzdVxsiLVO3TnPwsDU7vDMEREFGKKLCO1ZzxSe8ajyunB1n1F+Hl3AbILqgKOO5RXiUN5lVj500H0SIrCwNPiMaBnPGJa2aKOrZkQdcPUsQDl+6T20J1UM3Qn+3qbZABy9STw2kN1XHepbWMYIjKwmh/oAAAJkKvnYtQeflCF7+afQjvuL3RxopkZ1BQRNjNGDEjCiAFJKCxzVN8QthCFte6HJgAcyK3AgdwKfPLTQXTtGIEBvXxhqmNsWOga307UHroDBLzH7a/5/qgdmvxznap7lo6FJd+/onZXF7U6XGeoGtd8aDqun6GPkNSxOgQF/jQQ/p2a0ODxCoiav3yF8A9j+P71fe3fVv1Emu/PcP8+UXOcQJ3fC0KIY3+1137uJgjVOkN6E0Igp8iOrft8E6zrW7uoRocYG1J7xiGlRxy6JUZB0emqtPZSy5YQ0LOE6t6j6mBkNiuIjg5DZaVvja+a3qWa4yBqHivxNiinoOc6QwxD1fgLvOkYhvTRVut44nUDA3fUPq52ABLCtxCkJgIDlxCAx6PCrWqBC0Q2ICy1x1/gQggcKazCtt+L8dv+YhSVO094bJjVhH7dYtG3eyz6do1BuM0c9HnbYy1DoXYdNU345y8Bx743ZElCbKRVtyDbXnHRRSJqdU78Z9Vxk1/rOa5myEFW6v/hL9nMUDUNHlXzP15UX07t8qjweEXApOP2TJIkdO0Yia4dI3F+ZjfkFtux/UAJth8oRk6RPeBYh8uLX/YW4pe9hZAkoGvHSPTtFovTu8agS8dI/rINsbrzl6rp8/udGoFhiIhaPSF8wwlWkxKwXZKAiDAzXB4VDqcXHlVA04wVjJITIpCcEIFzzuiKkgondhwsxc6DJdifUw61Vk+aEMDh/Eoczq/El5uzYbMo6N05Br27RqNPlxgkRNvaxG1BiJoDwxARtVk1mcdiUmCNMkEIDW6vb8K3x6sCkm8iq0mRIfnXkfENw7XHwBQXZcNZA5Nw1sAkON1e7Mkuw+7Dpdh9qBQVtdYxAnwLHv52oBi/HSgGAMREWHBa5+jqj5h2uaYR0YkwDBFRu+Cb/ijBYvLd5kGymWEySYiNjYBVBryqBlSHIFUTcLq8cLrVgN6T9sRmMSHttASknZYATQjkFtmxJ7sUuw+X4lBeZZ3XXVblxs97CvHznkIAQGykBb2So9ErORo9k6KQGM+r1Kj9YhgionbJ1wNUe7a27x+5+tJnS4QV4VYNVS4PXB61XU8KliUJnTtEoHOHCIwZ3AUuj4r9R8ux70gZ9hwpQ36Jo85jSisDw1GEzYQ+3eLQpUM4unWMROcOETCbOLmF2geGISIyJCGE78qeCCtUTfP3EmmabxkBVdXa7aowVrOClB6+y+8BoMLuxu9Hy/H70XLszykPWNOoRpXTi6w9Bcja4/takX0Bq1tipP8jLsrKeUfUJjEMEZGh1UzOjrDV/Dj0ravkcHrhcKv+CdntcIqRX1S4Bel9OiC9TwcAQHmVG/tzfMHoQG5FvT1Hqib8E7JrhFtN6NIxAl06RqJrR19PVEyEhQGJWj2GISIi1L7k37f2S0SYGRFhZnhVAVXT4PZovuE0TQSud9QORUcEhiO704vsgkrkljqw60AxsvOr4FHrroFld/kmbe/JLvNvi7CZ0LmD74q3zh3CkZQQgQ7Rtuo7zBO1DgxDRET1qAlHiixBkRVYzQqiYIbb6+s1cnlV/0ra7fHKtNrCbSb07xmH4TFdUDbYDrdHRW6RHYfyK5GdX4lD+ZUoqmdoDfANrx0fkMyKjE7xYUhKiEBSfBg6xYejU1w4IsOCXxSSqCkYhoiIGqAm75gVGZZIKzShQdN8w0WqpsGrCrgN0nOkyDK6dIxEl46RwADfNofL13t0tLAK2flVOFJYidJKd72P96gasguq6tyINsJmQmJcODrFhSExLgwd48KQGBuGyDAzh9qoWTEMERE1khC+oTRFRvUqzjIkCRAww+PV4Pao/ivU2vt8oxphVhNO7xqL07vG+rdVOT04WliFnEI7jhZVIaeoCoVlzhPWo8rp9c9Vqs1mUdAxNgwdY23oEBOGDjE2dIgNQ0K0jVe0kS4YhoiIdBDQc2SSERVugcfru4VITUAywmTs2iJs5joBye1VkV/sQG6xHTlFduSV2JFbbIfdefy94Y9xutU6k7VrxERYEB9tQ0KMDQnRVsRH23wfUVaEWfkrjhqG7xQiIp0JUX3pfvV8ozCLAsAMVRPwqBrcbhUujwZNtP8hteNZTAq6Jkaia2Kkf5sQApUOD/JLHMgrcSC/xI78UgcKShyoOklIAnyLRZZVX/12PJtFQXyUFXFRNsRFWREbZUVcpAWxUVbERlphsygcfiMADENERM1O1Frw0WpSYDMrEMI3d8ZTPaTmrZ5rZJReo9okSUJUuAVR4Rb07hITsK/K6UFhqRMFpQ4UlDpQWOZEYZkDxeWuU64e7nSrOFpkx9HjbmBbw2pWEBNpQUyEBTGRVsREWBAbaUF0hO8jJsICq5mByQgYhoiIWtjxQ2oRYRaomoDbq8Lr1eDyqtAMNN/oZCJsZkQkmdEjKSpgu6YJlFa6UFTuRFG5E8Vlvs+Ly50ornDB46176f/xXB4V+SWOetdRqmExyYiOsCAq3Fwd2I77N8yMqHAzbFYTZIamNothiIgohHxhR0CWAJtZgWTxXcJfM6Tm9fo+PKqAEAxINWRZ8s8POv24fTXDbiUVrjofpZW+D28Db7/i9mrVvVH1Lx3gb48kITLMhMjq9an8/9rMiAgz+UJdmAlR4RZYrJbqe+lRa8EwRETUitQ3pAZIEBBQqxeA9Hq16snZDEj1qT3s1r1TVJ39QghUOb3VwciNskoXyqrcKK9yo6zSjXK77/PG3MRXEwLldg/K7Z4GHS/LEsKtJoTbTAH/hll9gerMAUnoXmteFTUvhiEiolaspucICFwAsk5AUn33U/OovvWPfAGJIak+kiQhsrr3pmvH+o/RhIDd6UV5lS8cVdg9KK9yo6L68wq7G5UODyrsnkaFJv/za77eq0pH/eHpy83Z+NeM4YiPtjX6uanxGIaIiFoxTQgcyqtApd2DyHAzuneKggwJNQFJFRo+W38QJRUuJMaF49zMblAkyb8gpFdVcaSwCpV2D8xmGYmx4b4nFoCALyxpQiCnsAp2pxfhNhOSO0QEzH9RNQ1Z+4rhdGuwWWQM7BUPRW74+j6nev6GHtOUc6iahq17C1Fa6UZspAWD+nQ46WuQawWmzog44Tkk+CZqVzh8ASk737fYpKppkCDB7vSiyump/vDC4fQ26AbAbq+GfUfLEBtl5VykFiCJEA9capqGF198Ee+99x4qKiqQkZGBBx98EN26dav3+JKSEvzzn//Ed999B0mSMHnyZMyZMwdhYWFNakdJSRW8DZhwRydmMsmIi4tgLZuIddRPW6/ljgPFWLnuIHKL7VBVAUWRkBQfjsnDe6B/z3i8sXoHvt+aA63WS5NlYPSgZEyb2B87DhRjzabDKCxzQoIEi9l3G4yxg7ugZ3I0VE3g96Nl+PHXXBSUOuD2aoDQEBtlw1mpndAzOQbfbz2Cr7ccgcOlQgCQ4LtkfeyQzhid3uWUr2Ffdim+yTqKwlIHVA1QZKBDbBjGpndG7+r1hxpyTFPO8X3WEXzz81E43cG9hoaco6GvQZIAs9WMnHxfwK1yemB3euFweXGkoBIHcivgUTVEhVvgVbWA/98UKD4+Aoqiz6KbIQ9DL774It5++208/vjjSEpKwlNPPYXs7Gx8/PHHsFgsdY6/9tpr4XA48PDDD6O8vBz33XcfMjIy8MQTTzSpHW31h2Vr0tZ/8bQWrKN+2nItdxwoxhuf7YLT7UWEzQyTSYbXq6HK6YXNoiApPgxb9xWf8PGDescjt9hxwsdPO78fJAn4vy/3wqOqvom9ZgWaqsHl0RBmM6FXchR+2V0ITfh+uQshfLcd8fpuyjZmcGcMT01CzUieqO7zqPmtsi+7FB/8sB8uj4pwqwmKIkNVNdhdKqxmGReP6gUApzzmZIHoVOfo3Tkam3cV+F6DBNR0qqkCkCXg/MxupwxEpzrHmf0TsX5HfoNeg6JIiIkJR1mZL+DWd46OMWGIibYiv8ge8P+LgSiQnmEopMNkbrcbr776KmbPno2xY8cCABYsWIDRo0djzZo1mDJlSsDxP//8MzZs2IBVq1ahd+/eAIB//OMfmD59OmbNmoVOnTq19EsgItKdJgRWrjsIp9uL2Eirf50bi1mB2SSjqMyBrfuOXQ5eexCl5tfr1n3FiLDJSIgJr/P40ko3Vq47CCEEyqpciI20wu3R4Pb4AqMQAgdzyrHrYDEkSYbNLENRfB+SBJgUQBMStuwqwMQzu0OWZN8CktVzlIQmoAqBHYdKERdlRUSY2TccV72OkqZpKHd4sHFXPgAJJpOM2EgrIMF/89twq4Yyuwffb8vBaV1i6l3rRxMC32QdhcujIjrc4j9GNimIVmSUVbqwaVcBhABMMo49hwRIQsCrAd/8fBRnpSWfcMjsVOcor3Ljm5+PQlak+vfbPfgm6yh6dYk54XDX8ecwmWTIklTn/1e/HnEcMmsmIe0Z2rp1Ky6//HKsXr0avXr18m+/+uqr0bdvXzz88MMBxy9evBhvvPEGfvjhB/82t9uN9PR0zJ8/H5MmTQq6LZqmcaJhE0kSIMsya9lErKN+2motvaqGsio3JKDeEKBqokGXZkuSVH3vtEC1J1ZLUv3nqLl1CBAYtvzPUf1vVLgFNotS72uosLv9z+3/r1S3DfIJ2iCqj4sMM1eHlcDX7FUFqpwe/2sNbJyAVj0f6qSvQQDhNhMs5rqvAQBUVUOlwwNJklBzClHrCXxX8534NaB6XlZEmBmm6l4MWZYCVh73qhqqHMdeR83/k5pjhPD1ucVEWPzPQb466rUgZkh7hnJzcwEAycnJAdsTExP9+2rLy8urc6zFYkFsbCxycnKa1Ba5EZMB6eRYS32wjvppa7WsWQOn9i/gYNX/eMkfEuQmnkMIUe9QhVetDgnVbaj+te5PEkLg2FVYsgTpBNOKfYHA13t0PFVToQkB5fgXIFX/5xSBseYwSZJOeMNXIXy/dOVavUoB+zUJ0DR/r1ndJ/DVWpGlgCAjK8cOFkIcO8exufG1nk+qDr+SbsNCFCikYcjh8HXzHj83yGq1oqysrN7j65tHZLVa4XK5mqeRREQtzGpRkGgJD3UzmqQlXoPVoqBjC5zDamnaBTqnPocJVgsv7g6lkEZMm823foLb7Q7Y7nK56r06zGaz1Tm25vjw8Lb9g4OIiIhCI6RhqGbIKz8/P2B7fn5+vZOhk5KS6hzrdrtRWlqKxMTE5msoERERtVshDUMpKSmIjIzE+vXr/dvKy8uxfft2ZGRk1Dk+IyMDubm5OHjwoH/bhg0bAABnnHFG8zeYiIiI2p2QDlJaLBZMnToVTz/9NOLj49GlSxc89dRTSEpKwoQJE6CqKoqLixEVFQWbzYb09HQMHToUd955J+bNmwe73Y4HH3wQF110ES+rJyIioqCEfNFFVVXxzDPPYPny5XA6nf4VqLt27Yrs7Gycc845eOyxx3DJJZcAAIqKivDwww/j+++/h9VqxcSJE3HvvffCarWG8mUQERFRGxXyMEREREQUSlywgIiIiAyNYYiIiIgMjWGIiIiIDI1hiIiIiAyNYYiIiIgMjWGIiIiIDM0wYaioqAh33303hg8fjiFDhmDGjBnYt2+ff/+OHTswdepUDB48GOPHj8ebb74Zwta2Hfv378eQIUOwfPly/zbWsuHy8vLQr1+/Oh819WQtG27FihWYNGkS0tLSMHnyZHz66af+fdnZ2bj55psxdOhQjBo1Cs8++yxUVQ1ha1un9evX1/t+7NevH8455xwArGVjeL1ePPfccxg3bhyGDBmCa665Br/88ot/P7+/G6ayshIPPfQQRo0ahczMTMyePRtFRUX+/WvXrsUll1yC9PR0TJw4EStXrmz8SYRBXHnlleLyyy8XWVlZYu/eveKvf/2rGDVqlLDb7aK4uFiceeaZ4t577xV79+4V77//vkhLSxPvv/9+qJvdqrndbnHJJZeIvn37imXLlgkhBGvZSN98841IS0sTeXl5Ij8/3//hcDhYy0ZYsWKFSE1NFW+//bY4ePCgWLhwoUhJSRFbtmwRbrdbTJgwQcyYMUPs2rVLfP755yIzM1M899xzoW52q+NyuQLeh/n5+WLNmjWiX79+4v3332ctG+n5558XI0eOFN9//704cOCAuO+++8QZZ5wh8vLy+P3dCDfccIMYM2aM+Oabb8Tu3bvFrbfeKiZNmiRcLpfYu3evSEtLE88884zYu3ev+M9//iNSU1PFTz/91KhzGCIMlZaWilmzZoldu3b5t+3YsUP07dtXZGVliUWLFolRo0YJj8fj3z9//nwxYcKEUDS3zZg/f77485//HBCGWMvGeeWVV8QFF1xQ7z7WsmE0TRPjxo0Tjz/+eMD2G264QSxatEh8/PHHYuDAgaK0tNS/79133xVDhw4VLperpZvbplRVVYlx48aJe+65RwghWMtGuvDCC8Vjjz3m/7qiokL07dtXfPbZZ/z+bqDt27eLvn37im+//da/rbKyUgwbNkwsX75cPPDAA+Kyyy4LeMysWbPEDTfc0KjzGGKYLCYmBvPnz0ffvn0BAMXFxXj99deRlJSEPn36YNOmTcjMzITJdOxWbcOHD8eBAwdQWFgYqma3ahs3bsR///tfPP744wHbWcvG2bVrF3r37l3vPtayYfbv348jR47gggsuCNi+ZMkS3Hzzzdi0aRMGDBiAmJgY/77hw4ejsrISO3bsaOnmtimLFi2Cw+HA3LlzAYC1bKSEhAR8/fXXyM7Ohqqq+O9//wuLxYKUlBR+fzfQgQMHAADDhg3zb4uIiECPHj2wYcMGbNq0CSNGjAh4zPDhw7F582aIRtxgwxBhqLYHHngAI0aMwMqVK/Hoo48iPDwcubm5SEpKCjguMTERAJCTkxOKZrZq5eXlmDNnDu6//34kJycH7GMtG2f37t0oLi7GNddcg7POOgtXX301vvvuOwCsZUPt378fAGC323HjjTdixIgRuPzyy/HVV18BYB2DVfNH4y233ILY2FgArGVj3XfffTCbzTjnnHOQlpaGBQsW4Pnnn0f37t1ZywaqryaqqiI3NxfFxcUnrKPD4UBJSUmDz2O4MDRt2jQsW7YMU6ZMwW233YbffvsNTqcTFosl4LiaG7+6XK5QNLNVmzdvHoYMGVLnL3EArGUjeL1e/P777ygrK8Nf//pXvPLKKxg8eDBmzJiBtWvXspYNVFlZCQCYO3cupkyZgldffRUjR47Erbfeyjo2wdKlSxEVFYUrr7zSv421bJy9e/ciKioKL730Ev773//ikksuwezZs7Fjxw7WsoHS0tJw2mmn4aGHHkJeXh6cTifmz5+PkpISeDyeeutY87Xb7W7weUynPqR96dOnDwDg0UcfRVZWFt5++23YbLY6Rat5M4aHh7d4G1uzFStWYNOmTfj444/r3c9aNpzJZML69euhKApsNhsAYODAgdizZw+WLFnCWjaQ2WwGANx44424+OKLAQD9+/fH9u3b8dprr7GOQVqxYgUuuugi/3sT4Pd3Y+Tk5OCuu+7C66+/7h/iSUtLw969e/HCCy+wlg1ksVjw4osvYs6cOTj77LNhNptxwQUXYNy4cZBlGVartU4da74OCwtr8HkMEYaKi4uxdu1anH/++f7xWVmW0adPH+Tn5yMpKQn5+fkBj6n5ulOnTi3e3tZs2bJlKCoqwtixYwO2P/TQQ1i1ahVr2UgRERF1tp1++un44YcfWMsGqqlFzZzAGn369ME333yDzMxM7N69O2Af63hyO3fuxOHDh+v0/iYlJbGWDZSVlQWPx4O0tLSA7enp6fjuu+/QuXNnfn83UO/evbFs2TKUlpbCZDIhMjISl112GYYPH47k5OR66xgeHo6oqKgGn8MQw2SFhYWYNWsW1q5d69/m8Xiwfft29O7dGxkZGdi8eXPAWhnr1q1Dr169kJCQEIomt1pPP/00Vq1ahRUrVvg/AGDmzJl49NFHWctG2LNnD4YOHYr169cHbN+2bRv69OnDWjbQgAEDEBERgaysrIDtu3fvRvfu3ZGRkYHt27f7h9MAXx0jIiKQkpLS0s1tEzZt2oSEhIQ69WEtG65mHsuuXbsCtu/evRs9e/bk93cDVVZWYurUqdi5cydiY2MRGRmJ7OxsbN++HSNHjsSwYcOwYcOGgMesW7cOQ4cOhSw3IuLocelbWzB9+nQxYcIEsWHDBrFr1y4xa9YskZGRIY4cOSIKCwtFRkaGmDt3rtizZ49YtmyZSEtLE8uXLw91s9uE2pfWs5YNp6qquPTSS8WkSZPExo0bxd69e8W//vUvMXDgQLFr1y7WshFeeuklMWTIEPHxxx8HrDO0bt064XQ6xbnnnituvPFGsWPHDv/aOC+88EKom91q3XvvveK6666rs521bDhVVcXVV18tJk6cKNauXSv2798vFixYIPr37y9++eUXfn83wp/+9CcxdepUsXv3brF161YxZcoUcf311wshhNi9e7cYMGCAeOqpp8TevXvFkiVLuM7QyZSXl4uHHnpIjBw5UgwaNEjccMMNYvfu3f79WVlZ4oorrhADBw4U48aNE2+99VYIW9u21A5DQrCWjVFQUCDuueceMXLkSJGWliauvPJKsXHjRv9+1rLhXn31VTF+/HgxYMAAceGFF4rPP//cv+/AgQPi+uuvF2lpaWLUqFHi2WefFaqqhrC1rdv06dPFHXfcUe8+1rLhSktLxbx588TYsWPFkCFDxJVXXinWr1/v38/v74bJzc0Vt912mzjjjDPEiBEjxEMPPSQqKyv9+7/99lsxZcoUMXDgQDFx4kSxcuXKRp9DEqIRF+ITERERtTOGmDNEREREdCIMQ0RERGRoDENERERkaAxDREREZGgMQ0RERGRoDENERERkaAxDREREZGgMQ0REjcCl2YjaH0PcqJWImu6ee+7BBx98cNJjMjMz8dZbb7VQi1ree++9h3379uGee+4JdVOISEdcgZqIGuTQoUMoLi72f71w4UJs374dL774on9bZGQk+vTpE4rmtYjx48cjMzMTjz/+eKibQkQ6Ys8QETVI9+7d0b17d//X8fHxsFgsGDx4cOgaRUSkA84ZIiLdbNq0CVOnTkV6ejoyMzMxd+7cgN6k5cuXIy0tDZs2bcKll16KtLQ0nH/++fjqq6/w+++/Y9q0aUhPT8d5552HlStXBjyuX79+yMrKwsUXX4xBgwbhggsuwOrVqwPO73K58OSTT2LMmDEYOHAgLrjgAqxatSrgmPHjx+Nf//oXpk2bhkGDBuG+++4DAOzcuRO33347hg8fjgEDBmD06NH45z//CafT6X/ckSNH8MEHH6Bfv37Izs7GCy+8gH79+tWpQ79+/fDCCy8AALKzs9GvXz+89tprmDhxItLT07Fs2TIAwO7du3HzzTdj6NChGDp0KG677TYcPnxYh/8TRNQYDENEpIuNGzfiuuuug81mw7PPPou///3v2LBhA/785z/7AwUAeL1e3HXXXbjqqqvw8ssvIywsDLNnz8Ytt9yCsWPHYtGiRUhMTMTcuXORm5sbcI6bb74Z55xzDl588UX06tULd9xxB7799lsAvonNt912G959911cf/31ePnllzFkyBDceeedWLFiRcDzvPPOO0hLS8PChQtx2WWXIT8/H9dccw0cDgcef/xxLF68GJMnT8Zbb72FN998EwDw4osvomPHjhgzZgz++9//IjExsVH1eeGFF3DTTTfhySefxMiRI7F//35cddVVKCoqwhNPPIFHH30Uhw8fxtVXX42ioqIg/g8QUbA4TEZEupg/fz569eqFf//731AUBQCQnp6OyZMnY9myZbjmmmsAAJqm4ZZbbsHll18OACgvL8edd96JadOm4frrrwcAREVF4dJLL8W2bduQlJTkP8e1116L2267DQAwevRoXHzxxXjppZcwZswY/PTTT/j++++xYMECTJo0yX+Mw+HA008/jSlTpsBk8v3I69y5M2bPnu1/3h9++AH9+/fHc889h8jISADAWWedhR9//BHr16/HjBkzkJqaCovFgvj4+KCGBv/whz/g0ksv9X991113ISwsDK+//rr/nCNGjMC5556L//znP5g7d26jz0FEwWHPEBE1mcPhQFZWFsaMGQMhBLxeL7xeL7p164bevXvjxx9/DDh+yJAh/s8TEhIA+IJTjdjYWAC+oFTbxRdf7P9ckiScd9552Lp1K5xOJ9auXQtJkjBmzBj/+b1eL8aPH4+CggLs2bPH/9j+/fsHPO+oUaPw9ttvw2q1Yu/evfjyyy/x8ssvo7i4GG63u2nFOcE5161bh8zMTNhsNn9bIyMjMWzYMPz000+6nJOIGoY9Q0TUZOXl5dA0DYsXL8bixYvr7LdarQFf1/SE1BYWFnbK8xw/NJWQkAAhBMrLy1FaWgohBIYOHVrvY/Pz8/2BJDw8PGCfpml45pln8M4778ButyM5ORmDBg2q0+6mOP6cpaWlWLVqVZ05TYBvcjoRtRyGISJqsoiICEiShOuuuw6TJ0+us78hQachSktL0aFDB//XhYWFUBQFsbGxiIqKQnh4uH+Oz/F69Ohxwud95ZVX8Prrr+Phhx/GhAkTEBUVBQC47LLLTtoeSZIAAKqq+ocGq6qqGvRaoqKicNZZZ/mHBmurGc4jopbBYTIiarLIyEikpqbi999/R1pamv/j9NNPxwsvvID169frcp4vvvjC/7kQAmvWrMEZZ5wBi8WCzMxM2O12CCEC2rB792689NJL8Hq9J3zezZs3o0+fPrj00kv9QSgvLw+7d++Gpmn+42Q58EdmTQ9X7YnemzdvbtBryczMxN69e9G/f39/WwcOHIjXX38dn3/+eYOeg4j0wTBERLqYNWsWfvjhB9x111349ttv8dVXX2H69OlYu3YtBgwYoMs5nnzySbzxxhv47rvvMHPmTOzbtw9/+9vfAABjxoxBRkYGbr31VixduhTr16/H4sWLMW/ePMiyfNKhp0GDBmHXrl145ZVXsGHDBrz33nu45ppr4Ha74XA4/MdFR0dj+/bt2LBhA5xOJ8aMGQMAePDBB/HTTz9h2bJlmDdvHiIiIk75Wm699VYcOnQIN998M7744gt8//33+Otf/4qVK1ciJSWliZUiosZgGCIiXYwaNQpLlixBbm4uZs6ciTlz5kBRFLz22mu6Lcw4b948/O9//8Ptt9+OgoICvPrqqxg2bBgAX6/NK6+8gsmTJ+Pf//43brzxRv9l9gsWLDjp89588824+uqr8eabb+Kmm27CkiVL8Mc//hG333479uzZ45/IfcMNN6CwsBA33ngjtm3bhl69euGJJ55AdnY2ZsyYgTfffBOPPPJIgy67T0lJwTvvvANJkjBnzhzMnDkTBQUFeOmllzBhwoSmF4uIGoy34yCiVm/58uW499578eWXX6Jr166hbg4RtTPsGSIiIiJDYxgiIiIiQ+MwGRERERkae4aIiIjI0BiGiIiIyNAYhoiIiMjQGIaIiIjI0BiGiIiIyNAYhoiIiMjQGIaIiIjI0BiGiIiIyND+HyU6nthq84qKAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "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 am not able to reproduce the figure below, maybe due to updates in `Seaborn` package, which does not compute parameters for regression when there is a perfect separation." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VPW9+P/XObNkmewhC0vYA4RFXBAQV6LIvrhVkSoWUWsv9ddqq3bzWnurtbetxfZ+Vdy3LooLQhSXoFAVURSJ7GsgAbJvk8w+5/z+mGRISIBJyGSWvJ+Ppsk5c+bk8zHMvOezvT+Krus6QgghxAnUUBdACCFEeJIAIYQQokMSIIQQQnRIAoQQQogOSYAQQgjRIQkQQgghOhS0APGLX/yCCy64gDlz5nT4uK7r/M///A/Tpk1j7ty5bN++PVhFEUII0QVBCxBXX301zzzzzEkf37BhA8XFxXzwwQf87ne/48EHHwxWUYQQQnRB0ALE+eefT3Jy8kkfLywsZMGCBSiKwtlnn01DQwMVFRXBKo4QQohOMobqF5eXl5Odne0/zs7Opry8nMzMzFM+z2pzoQAoCoDv5+OHzT+3ekzxHbe9TkFVfAfN3/zPEUII4ROyANFRho9A3qQbbW6qqhu7vTz+YNJcDkUBVVH8P7c+pzYfq4qCqoKqKv5rz1RGRiKVldYzvk+4iub6RXPdQOoX6TIyEjv9nJAFiOzsbMrKyvzHZWVlp209BJPe/H++73rrswFTFDAoii9gNH8ZmoOHQVUwGBQMqkwcE0JEhpAFiPz8fF555RVmz57N1q1bSUxMDGmA6A66Dh5dB+3kgUWB5mChYlAVjAYVg0HBZFBRVenmEkKEj6AFiLvvvpsvv/yS2tpaLrnkEn784x/j8XgAWLhwIZdeeinr169n2rRpxMXF8fDDDwerKGFFBzyajkfztntMVUAxGWmwuTAZVMwmVVocQoiQUSIt3fexqqagjEGEi7Q0CzU1Tf5jVVWIMaqYTQZiTIaIb2VEcz9vNNcNpH6RLqLGIERgNE3H7vJid/laHGajSqzZQKzZGPHBQggR3iRARBiXR8Pl0bDa3JhNBuJjjMSYDaEulhAiCkmAiFA64HR7cbq9GFSF+FgjcTFGVFnPIYToJjICGgW8mo7V5qaqzo7N4Ql1cYQQUUICRBTRdGiwuaiqs+N0t58lJYQQnSEBIgp5NJ1aq5P6JhdaZE1SE0KEEQkQUczu9FBV75DWhBCiSyRARDmtuTVhtblCXRQhRISRANFLNDk81DQ40E6RBkQIIVqTANGLuDwaVQ0O3B7pchJCnJ4EiF5G03RqGpzYnTIdVghxahIgeiEdqG9y0dDk6nBfDiGEAAkQvZrN6aHW6pSpsEKIDkmA6OVcHo2aegdeTQt1UYQQYUYChMCj6VQ3OHF7JEgIIY6TACGA5sFrq8xwEkIcJwFC+Ok61FqdeLzSkhBCSIAQJ9B0qJOBayEEEiBEBzyaTp3VKVNghejlJECIDrXsWieE6L0kQIiTsjk9OFyy4lqI3koChDilhiaXrJEQopeSACFOSdOhvlFShQvRG0mAEKfl8mg02mU8QojeRgKECEiT3S3rI4ToZSRAiIDoILOahOhlJECIgDndXpwuScUhRG8hAUJ0itUme0gI0VtIgBCd4tF0bLIbnRC9QsQFiD+89BV7SupCXYxerdHuRtOkFSFEtDOGugCddeBIPQeO1HP28D7MumAQCXGmUBep19F1aHK4SYw3h7ooQoggirgWhNL8/dt9Vfz1ta18u7dK+sRDwOb0SMZXIaJcxAWI+24+n+y0eMD3JvXax/t45YM9WG2y2rcn6TrYHDIWIUQ0i7gAMbR/Mv919VimTcjBoPraEzsP1fLX17fy7T5pTfQkm8MtrQgholjEBQgAg6oy9dz+LLt6HP0zLADYnV5eW7ePfxXulU+2PUTTwS4zmoSIWkENEBs2bGD69OlMmzaNFStWtHv86NGj3HTTTSxYsIC5c+eyfv36Tt0/Ky2eH84fy/SJx1sT3x2o4fGVW9lbKjOdekKTwyOtNiGiVNAChNfr5aGHHuKZZ56hoKCANWvWsG/fvjbXPPHEE8ycOZO3336bxx57jN/+9red/j0GVeHSs/vzo6vG+scmGmxunn93FwUbiyV/UJBpmo7dKaurhYhGQQsQRUVFDBo0iJycHMxmM7Nnz6awsLDNNYqi0NjYCIDVaiUzM7PLv69vuoUfXTWWi8/q65/p9Nl3ZTy5ajtVdfYu31ecns0hOZqEiEZBWwdRXl5Odna2/zgrK4uioqI21yxbtoxbb72VV155BbvdzvPPPx/QvdPSLCd9bNGs0UwYk81zq3dQ3+jkaFUT//fWNhZOH8nksX27Vpkedqr6havE5FhizYH9c8rISAxyaUInmusGUr/eJmgBoqN+aUVR2hwXFBRw1VVXsWTJErZs2cK9997LmjVrUNVTN2xqappO+XhGYgzLrh7Lm+sPsPNQLU63lxfW7GDb3krmXjgEkzF8x+bT0iynrV84arI6SE2MOe11GRmJVFZae6BEPS+a6wZSv0jXleAXtHfK7OxsysrK/Mfl5eXtupBWrlzJzJkzATjnnHNwOp3U1tZ2y++3xJr4/pUjmDtlsH8Ae/PuSp54e5t0OQWB0+2VrUmFiDJBCxDjxo2juLiYkpISXC4XBQUF5Ofnt7mmb9++bNy4EYD9+/fjdDpJS0vrtjIoisIFY7P54fwx/k+3ZTU2/v7Wd2w7WNNtv0f4yPRiIaJL0AKE0WjkgQceYOnSpcyaNYuZM2eSm5vL8uXL/YPV999/P6+99hrz5s3j7rvv5g9/+EO7bqju0D8jgWVXj2P04FQAXG6Nf3y4h/e/PCxJ57qR3SlTXoWIJooeYa/oY1VNVFU3dum5uq7z6XfHeH/TYVriwvD+yVx/+XAsseGR9C9SxyBaJFvMxMWcfGgrmvt5o7luIPWLdGE1BhGOFEXh4rP68YPZeVhifW9i+47U8//e2kZZjS3EpYsOsrJaiOjRqwJEi2H9kvmvq8cxoDlNR63VyZNvb2NHsYxLnCmXR5PFiUJEiV4ZIABSEmK4be4YzsntA/je2F75YA/rvimVfvQz5JB9q4WICr02QACYjCrXXjaMmZMH0jI2/tHmUl77eB9uj3wK7irpZhIiOvTqAAHHxyUWzxhFjMkAwNZ91TxbsINGu6SQ6AqvpuNySytCiEjX6wNEixE5KfxwwfH1EofLG3ni7W1U1Mqiuq6wSzeTEBFPAkQrWanx3LlgLAOzEoDmwetV2zh4rCHEJYs8TpesiRAi0kmAOEFCnIlbZ4/mrGHpgG/A9bmCnWzdVxXikkUWTfel3xBCRC4JEB0wGVW+lz+cS8/uB/j61P+9bh//2Xo0xCWLLLJPhBCRTQLESaiKwvSJA5l/0RD/DKf3Nh3m3S8OyT7MAXK5vZLKRIgIJgHiNCaNzuL7V47EZPD9p/q06BgrP94vmUsDoAMOl0x5FSJSSYAIQN6gVJbMziMuxjcN9tt9Vby0drdM5QyATdZECBGxJEAEaFB2IrfPG0OyxQzA3tJ6nnt3pywKOw2PV8ftkUAqRCSSANEJWanx3DF/DBkpsYBvrcTTq3dgtblCXLLwJvtECBGZJEB0UkpCDLfPG0P/Pr5Ef2U1Nla8s4NaqzPEJQtfDpdXBvaFiEASILrAEmvi1jl5DOnry69e3eBgxTvbqaqXVdcd0QGHdMUJEXEkQHRRrNnILTPzGDUwBYD6JhdPv7ODctlXokMyWC1E5JEAcQZMRpUbp41g7FDfPtpWu5unV+/gaFXk7ggXLDJYLUTkkQBxhowGlevzc/37SticHp5Zs4PSiq5tixrNbLKyWoiIIgGiGxhUhWsuG8bEvEzANyj7bMFODpdH7/62XSEJ/ISILBIguomqKMy/aAgXjMkGfInqnn93F4fKJEi0kAR+QkQWCRDdSFEU5kwZxIVjWwWJ93ZKkGhFFhYKETkkQHQzRVGYdcEgLj6rLwAutyZBohWHjEMIETEkQASBoijMmDSQS8YfDxIvvLdLxiQATdelm0mICCEBIkiU5nThLS2JljGJkgoJEg7ZjlSIiCABIohaWhIXjWsbJI708nUSTkkBLkREkAARZIqiMHPyQKY0D1y3bGFa1otXXGs6kipdiAggAaIHKIrC7AsGMWl0FuCbyfPsmh1U1Pbe3E3SzSRE+JMA0UMURWHuhYM5b2QGAE0OD88W7KC63hHikoWGQ1oQQoQ9CRA9SFUUrrp4KGcP96XlsNrcPFuwg7rG3pcqXNN06WYSIsxJgOhhanNajjFDfAn+6hpdPFuws1duOiStCCHCmwSIEDCoCtfnD2dEji9VeHW9g+cKdmJzuENcsp4l4xBChDcJECFiNKgsmjaCIX2TACivtfPCe7t61cY60s0kRHiTABFCJqPKzdNHkpOZAEBpZRNPvFmE26OFuGQ9R7qZhAhfAQWIW2+9lY8//rjTqZo3bNjA9OnTmTZtGitWrOjwmnfffZdZs2Yxe/Zs7rnnnk7dPxrEmA0snjGKrNQ4AHYfquVfhXvxar0jSEg3kxDhK6AAcf311/Piiy9yxRVXsGLFCmpra0/7HK/Xy0MPPcQzzzxDQUEBa9asYd++fW2uKS4uZsWKFfzzn/+koKCAX/7yl12rRYSLjzXyg9l5pCXFALDzUC1vrj+A1gv2TtA02WlOiHAVUIC48soreeGFF3j66aepqKhgzpw53HvvvWzbtu2kzykqKmLQoEHk5ORgNpuZPXs2hYWFba557bXXWLRoEcnJyQCkp6efQVUiW1K8mVtn55Gc4AsSW/ZW8e7GQ71igx27tCKECEvGrjzJZDIRExPDfffdx8UXX8z999/f7pry8nKys7P9x1lZWRQVFbW5pri4GIAbbrgBTdNYtmwZl1xyyWl/f1qapSvFDntpaRb+vxvO5s+vfE2Tw8Pn28rokxbPrClDQl20bnXi38+gKmSkR8ffNCMjMdRFCCqpX+8SUID44IMPeOWVV6iurubGG2+koKAAi8WCx+Phyiuv7DBAdPTJV1GUNsder5dDhw7x8ssvU1ZWxqJFi1izZg1JSUmnLE9NTfQmu+vXJ4GbZ4zk2TU7cXk03tlwADSNyaOzT//kCJCWZunw76e53JiMhhCUqPtkZCRSWRm92XqlfpGtK8EvoACxcuVKbrvtNi6++OK2TzYa+fWvf93hc7KzsykrK/Mfl5eXk5mZ2eaarKwszj77bEwmEzk5OQwZMoTi4mLOOuusztYjquRkJrLoyhG8tHY3Xk1n9afFxMeYOGtY9HbB2V3eiA8QQkSbgMYgnnrqqXbBoUV+fn6H58eNG0dxcTElJSW4XC4KCgraXXvFFVewadMmAGpqaiguLiYnJ6cz5Y9auQNS+F7+cBRAB17/eB/7jtSHulhBI7OZhAg/AQWIG2+8kfr6429OdXV1LFq06JTPMRqNPPDAAyxdupRZs2Yxc+ZMcnNzWb58uX+w+uKLLyYlJYVZs2axePFi7r33XlJTU8+gOtFl3NB05l40GACvpvPKB7s5UtkY2kIFiabJTnNChBtFD2CazPz581m1atVpz/WEY1VNVFVH55skdNxHX/h1KYVflwJgiTVyx/wx9EmOC0XxztjJxiAA4mKMJFvMPVyi7tMb+rClfpGrK2MQAbUgNE3DZju+wU1TUxNer3za6yn55/b37yXR5PDw/Lu7ojK5n9Pl6RXTeoWIFAEFiDlz5rBkyRJWrVrFqlWruPXWW5k3b16wyyaaKYrC3CmDGducAbbW6vTlbYqyrTs1HelmEiKMBDSL6Y477iAzM5N169ah6zo33HADCxYsCHbZRCuqqnDd1OE0OXZx8FgDx6ptvPLBHm6ZOQqjIXpSajlcXmLNXVqeI4ToZgGNQYST3jgG0ZrD5eHp1Ts4Vu3r8hs3NI3rL89FPWGNSbg6Xf0UICM1LmLq01pv6MOW+kWuoK2DqK6u5uWXX6akpASP53i3xvLlyzv9C8WZiTUbWTxzFE++vY26RhffHaghMf4Qsy8Y1G4hYiTSAafLS1yMtCKECLWAXoU//vGPGTZsGBdccAEGgyxmCrWkeDNLZuXx5Krt2Jy+lBxJFjOXjO8X6qJ1C4cECCHCQkCvwoaGBn73u98FuyyiE/qkxLF45kieWb0Tt1dj7abDJMabOCc3I9RFO2NOtxevpmFQo2dsRYhIFNArMDc3l/Ly8mCXRXRSTmYiC6/IRW3uWXrjkwPsLa0LbaG6id0ps5mECLWAWxDz5s3jnHPOISYmxn9exiBCb9SgVBZcPJQ3N/j2j3j1wz3cPncM/fpEdnZUm9NDQpwp1MUQolcLKEDMmTOHOXPmBLssoosmjMqkwebio82luNwaL763izvmjyEtKTbUResyTdNxurzEmGXMS4hQCShAXHXVVcEuhzhDU8/pT32ji692VWC1u3mhOUhYYiP3U7jN6ZEAIUQIBTQGUVxczMKFC/3ZWLdv387f/va3oBZMdI6iKMy7aAh5g3zJDqvqHby0djeuCN7Os2WwWggRGgEFiAcffJA777yTxETfQou8vDzWrl0b1IKJzjOoCtdfPpyBWQkAlFQ08u/CfXi1iFoL2YYMVgsROgEFCKvVyiWXXOJfiKWqKiZT5HZdRDOz0cBN00eSnuwbf9h5qJbVnx2M2CR4Nmd05ZsSIpIEFCAMBgNut9sfIMrLy1FljnrYssSa+MHMUViaZwF9ubOCT7YcDXGpukbT9KhLSihEpAh4w6Bly5ZRW1vL3/72N2688UaWLFkS7LKJM5CWFMstM0ZiNvr+xB9uLuGbPZUhLlXX2BwSIIQIhYBmMS1YsIABAwbw8ccfY7fbefTRR5kwYUKwyybOUP+MBG6cNoKX1u5C0+HN9QdIjDeROyAl1EXrFJdHw+3RMBml1SpETwo44c2ECRMkKESgETkpXHXJUN5Yf3wh3W1zx9A/whbS2Zweko2Ru9ucEJEooABxzTXXdJgpdOXKld1eINH9zhuZSX3T8YV0L0XgQjqH00NinAlVjfyMtUJEioACxH333ef/2el0UlBQQGZmZtAKJbpfpC+k05H0G0L0tIACxMSJE9scX3TRRTJIHWFaFtJZbW52Ha6lqt7By+/vZsnsPMzGyFitLAFCiJ7VpVG/xsZGSkpKurssIsgMqsINVwwnJ9O3kO5weWQtpJMpr0L0rE6PQWiaRmlpKT/4wQ+CWjARHGajgZtnjOTJVduprnf4F9LNv2hIROxIZ3N4ZM9qIXpIp8cgDAYDAwYMICsrK2iFEsHVspDuyVXbabS7+XJnBUkWM/nnDgh10U7L5dHweDWMBpnyKkSwdWkMQkS+tKRYFs8cxdOrt+Nya3y0uZSkeDMTRoX/5AObw0OSRaa8ChFsAQWIyZMnd9j9oOs6iqKwcePGbi+YCL7+fSwsmjaCF9/bjabrvP2fAyTEmxg1MDXURTslu8tDQrwJNQK6xISIZAEFiIULF1JXV8f111+Pruu88cYbZGVlMWvWrGCXTwRZ7oAUrrlsKK9/vB9Nh39+tJelc/LIyUwMddFOStfB4fQSHytjEUIEU0AduV999RX//d//zahRo8jLy+PXv/4169evp3///vTv3z/YZRRBdk5uBjMmDgTA7dF48b3dVNbZQ1yqU7M53aEughBRL6AAUVFRQU1Njf+4pqaGysrITPwmOnbx+L5MGZsN+NYbPP/uThqaXCEu1cl5vDout+wVIUQwBdRGX7x4MfPnz2fq1KkArF+/njvuuCOoBRM9S1EUZl0wCKvNzXcHqqlrdPHi2l3cNnd02E4rtTk9mE2RschPiEgU0Ct/0aJFnHfeeXz11Vfous6iRYsYOXJksMsmepiqKFw3dRhNDjcHjjZwrNrGy+/v4ZaZo8Iyk6rT5duS1CB7kwgRFAG/sgYMGMC5557LzTffLMEhihkNKt+/cgR90+MBOHisgdfW7UMLw9XWOrIlqRDBFFCAWL9+PbNnz+bHP/4xAN999x0//OEPg1owETqxZiO3zBxFWmIMANuLa3gnTLcttTncYVkuIaJBQAHi8ccfZ+XKlSQlJQEwbtw4Dh8+HNSCidBKjDfzg9l5bbYtLfy6NMSlak/TweGSVoQQwRBwF1NGRkabY7NZVrJGu/SkWH4wcxQxzQPB6745wsZtZSEuVXuyJakQwRFQgLBYLFRVVflXU2/atInExNMvpNqwYQPTp09n2rRprFix4qTXrV27lpEjR/Ldd98FWGzRU/r1sfD96SMwNG/Us/rzYr7dVxXiUrXl9moy5VWIIAhoFtM999zDbbfdRmlpKTfddBPFxcU88cQTp3yO1+vloYce4vnnnycrK4trr72W/Px8hg8f3ua6xsZGXn75ZcaPH9/1WoigGtYvmesvz+WfH+1B12Hlx/uJMxsYGUYpOWTKqxDdL6AWxPjx43nppZf405/+xNKlSykoKGDs2LGnfE5RURGDBg0iJycHs9nM7NmzKSwsbHfd8uXLWbp0KTExMV2rgegRY4ekseDioQBous4/PtzLoTJriEt1XMuUVyFE9zltC8Lr9fK9732PN954g0svvTTgG5eXl5Odne0/zsrKoqioqM01O3bsoKysjKlTp/Lcc88FfO+0NEvA10aicK3f9ClDQFV465P9uL0aL72/m3tuPJcBWZ3L2xSs+sXFmUhOCO0HjYyM8M1h1R2kfr3LaQOEwWAgNTUVp9PZqU/5HU09bJ0RVtM0HnnkER555JGA79mipqap08+JFGlplrCu34TcPlTV2PhP0THsTg+P/WsLd8wbTZ/kuICeH8z61QL2lNiQLZzLyEiksjJ8WlXdTeoX2boS/AIagxg8eDCLFi1i+vTpxMfH+88vWrTopM/Jzs6mrOz4jJfy8nIyM4/vNdDU1MSePXu4+eabAaisrOTOO+/kiSeeYNy4cZ2uiOgZiqIwY9JA7C4vm3dV0GR381zBTu6YNybkn951oMkue0UI0V0CChBNTU3k5uZy4MCBgG88btw4iouLKSkpISsri4KCAv785z/7H09MTGTTpk3+45tuuol7771XgkMEUBSFBRcNweHysO1ADXWNLp4t2Mnt88aQ0LxuIlTsTg+WOKOk3xCiG5wyQPzhD3/g/vvv55FHHuGzzz7jwgsvDPzGRiMPPPAAS5cuxev1cs0115Cbm8vy5csZO3Ysl19++RkXXoSOqip8b+pwnK7d7C2tp6rewfPv7mTpnNHExYQuuZ8ONNo9JEsrQogzpuinyFNw1VVX8dZbb7X7OZSOVTVRVd0Y6mIETbiPQZzI5fHywru7KG6e0ZSTmcCSWXnEmDuectoT9VOA9OTYHt+3ujf0YUv9IldXxiBO+QpqHTsk343oiNlo4OYZI+mf4ZuZVFLRyMsf7MbtCd2UU99YhGwoJMSZOmWAcLlc7N+/n3379rX5ueVLCPAl9/vBzFFkpfpmMh042sCrH+7G4w1dkLC7vLg9srpaiDNxyi6m/Pz8kz9RUTpc+BZs0sUUvqw2FytW76C63gFA3qBUbpyW22bAuCfrZzKopCfH9sjvgt7RRSH1i1zdPs113bp1XS6M6H0S483cOjuPp1fvoNbqZOehWv69bh/X5+f6czn1JLdXw+bwEB8bnjviCRHuZC6g6FYpCTHcOjvPP4to24EaVn4Sug2HGu2usNzsSIhIIB+tRLdLS4rl1jl5PP3ODqx2N1v3VaMqCtdcOqzL99xbWsfmXRXUWp2kJsYwYVQmuQNSTvs8TQer3S3TXiPItoPVfFp0jMo6OxkpcVx0Vl/GDkkPdbF6JWlBiKDokxzHrXNG+xfObdlbxZsbDqB1YTbc3tI63v+yhOoGJ5oO1Q1O3v+yhL2ldQE93+70hHRWlQjctoPVvLH+AOW1djQdymvtvLH+ANsOVoe6aL2SBAgRNJmpcdw6Jw9L8xjAN3sqeeW9nZ0OEpt3VXTqfEdsTtlUKBJ8WnSsU+dFcEmAEEGVlRrPrXNG+weKPy86xpvrD3RqXKDW6uzU+Y44nB4Zi4gAlXX2k5x39HBJBEiAED0gOy2epa2CxDd7Knlj/f6A37BTEztOAniy8x3RkVZEJMhI6TgrcEZKz01XFsdJgBA9oiVIJMYfH5NY+cl+vAEEiQmjMjt1/mRsTo9kBAhzF53Vt1PnRXBJgBA9Jjstnp8uPBdL88D1t/uq+Pe6vafdCS53QArTJ+aQnhSDqkB6UgzTJ+YENIupNU3TcbhkdXU4GzsknWsuHUpWahyqopCVGsc1lw6VWUwhItNcRY/ql5HAbXNG8+wa3xTYbQdq8Hj2svCKXEzGk39eyR2Q0umA0BGbwxPSbLPi9MYOSZeAECakBSF6XGZqHLfNG+1fm7DrcC2vfLAbVw/kTnJ7NcnRJESAJECIkOiTHMft80b7B5r3ltbzwru7cLiCP5Bsc8hgtRCBkAAhQiY1MZbb542hT3NCveIyK8+s2UljkFN1O1zeLi3YE6K3kQAhQirZYua2uaPpm+7b6/xoVRNPr95OfWPgaxw6S8e3LkIIcWoSIETIJcabWTpnNAOzEgDfoqin3tlO1UkWTXUHWRMhxOlJgBBhIS7GyJJZeeQOSAagrtHFk+9sp7QiOHt/eLw6LrcMVgtxKhIgRNgwmwzcNH0k44amAb7B5GfW7Ag4KV9n2aUVIcQpSYAQYcVoULk+P5fJo7MAcHk0Xlq7m2/3VXX773K4vJKfSYhTkAAhwo6qKsy9cDBXTBgAgFfTeW3dPtZ/e6RbU2Xo0CPTaoWIVBIgRFhSFIX8cwdw1cVDaNmt9P0vS3jns+Ju/dQvayKEODkJECKsnZ+Xxfenj/Sn4di0o5xXPtiDs5sGmD2aTpMjuOsuhIhUEiBE2Bs1MJXb5o72J/nbdbiWFe9031qJRpsbj1d2nBPiRBIgREQYkJHAnfPH+PcLOFZt44m3t3G0qumM760DDU2uM76PENFGAoSIGGlJsfxw/hiG9/etlWiwuXnqne1sO1hzxvd2eTSZ9irECSRAiIgSF2Nk8cyRnN+8WZDbo/GPD/dQ+HXpGedXstpcMu1ViFYkQIiIY1BVFlw8hFmTB6E0z3Aq/LqUf32bgCyaAAAfmklEQVS094xWR2s61Fgdp93ASIjeQgKEiEiKonDRWX1ZPGMUsWYDANsO1vDkqu1UN3R9g3uPV6emwSmD1kIgAUJEuBE5KfxowVh/yvCyGhv/9+Z37CnpenoOr6ZTY3Xi9kiQEL2bBAgR8fqkxHHngrGMGpgK+FJovPjeLtZ90/VxCU3TqbU6JEiIXk0ChIgKcTFGvj99BJefNwAF39TVjzaX8vLa3di6uBBO05EgIXo1CRAiaqiKwuXnDeCmGSP94xK7S+r4+5vfUVJh7dI9JUiI3iyoAWLDhg1Mnz6dadOmsWLFinaPP//888yaNYu5c+eyePFijhw5EsziiF5i1MBUll09jv59LIBvb4kV7+zgP0VHu9TldDxIyP4RoncJWoDwer089NBDPPPMMxQUFLBmzRr27dvX5pq8vDzeeOMNVq9ezfTp0/nf//3fYBVH9DJpSb79ric1pw33ajrvfXGYl9bu6tKe15oONQ3OLndXCRGJghYgioqKGDRoEDk5OZjNZmbPnk1hYWGbayZPnkxcnC91wtlnn01ZWVmwiiN6IZNRZf5FQ7g+fzgxJl+X056Sev62sqhLmxDp+FZv1zU6z3hRnhCRwBisG5eXl5Odne0/zsrKoqio6KTXr1y5kksuuSSge6elWc64fOFM6te9pk60MHZEJs+u2kbxsQasdjfPv7uL/Ak5XHXZMExGQ+dvalBJS47FYGj7GSsjI7GbSh2epH69S9ACREcbuygty15PsGrVKrZt28Yrr7wS0L1ras48QVu4SkuzSP2CwAAsmTWKjzaX8J+tx9CBdZtL2La/iuvzh9M3vfNBq6q6kdSEGH8q8oyMRCoruzYYHgmkfpGtK8EvaF1M2dnZbbqMysvLyczMbHfd559/zpNPPskTTzyB2WwOVnGEwGhQmTFpEEvm5JFs8f1bq6i18//e2sbH3xzB28k8TJqmU2N1dNveFEKEm6AFiHHjxlFcXExJSQkul4uCggLy8/PbXLNjxw4eeOABnnjiCdLT04NVFCHaGNYvmbuuPYuzhvn+zXk1nQ83l/Dkqm2U19g6dS9dhzqrU3amE1FJ0btzk98TrF+/nocffhiv18s111zDnXfeyfLlyxk7diyXX345t9xyC3v27CEjIwOAvn378uSTT57ynseqmqiqbgxWkUNOuph61tZ9VbzzWbE/1bdBVZh6bn8uGd8Po6Fzn59y+qXgsjtP2pUa6XpDF0y016+zghoggkECRGQLx/pt3V/FexsP0WA7PoU1NTGGtMQYnG4vqYkxTBiVSe6AlFPeJy3NQmODneQEMwa1c8Fl28FqPi06RmWdnYyUOC46qy9jh4RHq3rNxmI+2XKEJocHS6yRy87pz5wLBoe6WN1OAkR7QRukFiIS7C2t4z9bj2GJM6GqKg1NzuaFcU5qrU4ssUY8Xp33vywBOG2QcHk0quocxJoNxMea/APYp7LtYDVvrD/gPy6vtfuPQx0k1mwsZs1nxYBvkkmjze0/jsYgIdqSVBuiV9u8qwLwvfnFxxrJSI1HVY93ETU5PFTU2bE7PXy1szyge+qA3eWlusFBTYMDu9PT4ay+Fp8WHevU+Z70yZaOsxuc7LyILtKCEL1ardXZ5tigKhhU8IUIBa+mN2d29Q1EV9TZyWzeFzsQLo+Gy+PCaleIjzESF2No1/1UWWfv8LmVdV3f16K7nGzVeVMXVqOLyCMtCNGrpSbGtDtnUFVMRgMZqXEkxJn8551uL4+/XsS7XxzC4ercrCVN02m0u6mqc1DX6Gyz813GSQJORkpsp35HMLSuf2uWk5wX0UUChOjVJoxqvzYnPtaIJdaIqigkWcxkpMRhNvleKpqu82nRMf787618tbO803tY6/j2q6ixOqmqt2NzeLhwXHaH1150Vt9O16e7XXZO/06dF9HF8OCDDz4Y6kJ0RqPNjc3uCnUxgiYuzow9ipvv4Va/9KRYUhNjqLM6cbq8pCXFkH/uAEYOTPWfy0iJZfrEHEYPTqOkohGHy4vbo7HrcB3bD9aQlhRLenJsp+um6b5WSWKciczUOOqsThwuL5mpccyYNDDkA9Tg27EPBY5WN+HxaljiTFw5cWBUDlBbLDHYbNH73mKxtG8tn45Mcw0z4TgNtDtFev3cHo0NW4+yYevRNntEDO+fzHVXjCAxpgt5nVoxqAoxZgOxJgNm05ndq7v1hmmg0V6/zpIAEWYi/Q30dKKlfg1NLj7aXMLXuytp/QIaMySNaefndGog+2QMqkKs2UCs2RjQdNlg6w1voNFev86SWUxCdEGSxczVlw7jgrHZvP/lYfaU1AOw/WANO4prOHt4H6ae058+ZxAovJpOk8NDk8ODUVWIjTESYzKERbAQvYMECCHOQN90C7fMzOPgsQYKvznCgSP16Dps2VvFt/uqOHt4Hy47p/9JZyoFytM8C6rR7kZVwGT0BQqTQcVsUqM2vYcIrYjrYrI53JRXWNF0HV33TR/Udd9cdQ3QNZ2IqtAJoqUL5mSiuX6pqfFs/PYIH31dytGq43VUgDFD07js7P7069P9e2EogNlkIMakYjYZOp1DKlC9oQsm2uvXWRHXgoiPNZ10bnYLrTlg6Hrbn0GnZVairjcHEp3m78cDi64f389C03R/MIrkwCOCT1EURg1KZeTAFHYdqqXw61KOVtvQgW0Hath2oIYROclcdFY/hvVL6rZP/Tq+2VC+tONuVFUhxqhiMqoYDSpGo4oqLQzRBREXIAKhKgqqoftfELquNwec4z+3BKE2rZnmABMtrRrROYqikDc4jVGDUtldUscnW45wuNw3sWJPST17Surpmx7PRWf1ZdzQ9G7/xK9pOnaXF7vr+GI8o6pgMhkwtwocQpxOxHUxARHZDGwTULTWgcUXUDRdR9d0UtMsVFU1tgs+0SKau5hOVjdd1zl4zMr6b4+wt7S+zWOJcSYmjs5iYl4mifE9t2GWquBvXZgMKkaDgsFw6pZGb+iCifb6dVZUtiDCkaIoGFpefKeY3p6eHIfWQRqH1l1lHbZemgOMt/lLk1ZL2FAUhaH9khjaL4myGhufFh1j674qvJqO1e6m8OtSPtlyhLFD05iYl8Xg7MSgDzprekueKK3NeV8uKl+wMKgKJoOvxdE6gaHoPaQFEWa681OMV9Pwen0BQ9eh5T2n9diM/7um4+2B1kpvbEF0pKHJxaad5Xy5s6Jd4ruMlDgm5mVydm4fLLHhkfPIoCpkZyVRX2fDoCq+1odBiarZU9KCaE8CRJgJ9T9STdPxahoeb6uWSOtxlVY/n+5fTstbR+vLJEC05fFqfLe/mi92lFNS0XYBqEH1DXpPGJnB8AEpGEL8Kb6j+rUEC4OqoCi+1pKqKBgMiv+xSBHq116wSReTOGOqqqCqBkwB/stoGUdRUGj+H0CbT5atWzKJ8WZsjQ48Xg1NP359i+MzyaJr7OVkjAaVc0ZkcM6IDI5WNfHlznK+3VeFy63h1XS2H6xh+8EaEuJMjB+Wztm5fejXxxI2n9x9XZrekz6ugL/LymhoaXm0bX20fPBQUKQrK8xICyLM9IZPMYHW7/isMd/MMa318YlrYPTjU5hDpbtaR063l20Hqvl6dyXFZe3/W/VJjuWsYemMG5ZOVmr8Gf++QHV3609VlXYfBBSlZRzEN+6hNh+3tExU1fe8zm7pGoje8NrrLGlBiLDVMrDfmV6K1uMq/rGWDoLKic/xhtGgfozJwHkjMzlvZCZVdXb/quyWzY2q6h2s++YI6745QnZaPGOGpDFmSBpZqXFh07IIREep0nUdPF4dj/fkrRLwtUzU5gF1RfG1PFqCi6oo/haKqhzv+hKdJy2IMNMbPsWEa/1axl80zdd1ounNs8K8Gpqm4zlNEyWY4yuarnO43Mq3e6vYdrAGm6P9TLf0pFjyBqeSNyiVgVmJ3T5mEcnjR0rz/yn4goVC+4Wv6ekJ1NQ0+rpLOT6pQ2kJMs0/t7mPAgqK/1q1+XG1+RcqSvO6rDDoOpNB6igQzm+g3SGS6+fxajQ5PDicng5bGz31BurVNA4cbaBofzU7imuxO9sHi7gYIyNzUhiRk0JuTnK3zIaK5AARiGDWTwEU1dci9o3zKc2tneaWsnq8FdQScDrT6mm9iPd4l+vxljPA8MGd319EupiECJDRoJJsMZMQZ8Tm8GBzekIykG5QVXIHpJA7IIUFF2scPGb1Z5G12nxTZu1OD9/u83VNKUD/DAvDB6QwvH8SA7MSI2p2UTTQ8WVU0NDh1L1nfv5WT3OLp3Xw8K+HIrBMDV1tTUqAEKKTDKpKYrwZS5wJh9OXjtsbohFyg6oyvH8yw/snM/fCwRytamLXoVp2HarlaLUN8L05lVY2UVrZxCdbjmAyqgzOTvQv3uvXJyHkU2hFe/5ccW0+hfTsvzMJEEJ0kaooxMeaiI814dU0UpJj8ThcuD0aTo/W6f2qu6M8AzISGJCRwBUTcmiwudhbUseekjr2HanH7vR9dHV7NPaW1vvTfpiNKgOzEhmUncjgvonkZCSE3W52IjQkQAjRDQyqSqzZSHyrvn6PV/NlWW3ew7qn2xhJ8Wb/bChN0zla3cS+0nr2HanncLkVj9dXIpdHY98R33nwDbBmpcUzMCuRnMwEBmQm0Cc5todLL8KBBAghgqRlUZgl1oSm6TjdXmwOD26vdvondzNVPd66uOyc/ni8GiUVjRw42kBxWQOHyxv9e2xrOhyrtnGs2samHeWAb+rt4L5JZKbE0q+PhX59LKQnx0oa8SgnAUKIHqCqCnExRuJijDjdXprs7naJ8nqS0aAypG8SQ/omAb6ZUceqbBSXWTlcYaWkvJH6Jpf/eqfby+7Dtew+fPweJqNKdlo8fdPjyU6LJystnqzUeOJj5W0lWshfUogeFmMyEGMy4PZo2F0eHC5vj49XnMigqgxo7k6CvgDUNzqbB7cbKa1s5GiVrc2UWrfH1wo5MYdUYryJzNQ4MlLiyEzxfe+TEkdSvEkWrEUYCRBChIjJqGIymkmKB5fbi8PlxeHyhDxlSIvkhBiSE2IYMyQN8G2puv9QDUeqmjha1URZtY1jNTYaWrU0AKw2N1abm/1HGtqcN5tU+iTHkZ4UQ3pSLOnJsaQlxZKWGEOixSzdVWFIAoQQYcBsMmA2GUiymP0D2w536FsWrSmK4ntDT4pl3NDji66aHG7Ka2yU19gpq7FRUWunos7mnzXVwuXWONocXE5kNCikJMSQmhjT5ntygpmUhBiSLKag5F8SpyYBQogw09IFlYSvG8fp9uJyh2YmVCAssSaG9ktmaL9k/zld12m0u6mss1NZ56Cq3k5VnYOqBge1DU60E1YYerw6VfUOquodHf4OBUiIM5GUYCYp3kySxfc9Md7U/OX72RJrCou0FtFCAoQQYczUvIc0cSY0Xcft1nB5fMHC7dXCNiW6oijNb9rmNoEDfHmu6hqdVNc7qLU6qWlwUNPgpK7RSa3Via2D1CE6YLW7sdrdHOHk6TAUID7WSEKcCUucyfc91oQlzogl1kR8rJH4WN/PcTFG4mOMvv++okMSIISIEKqiEGM2EGM+vojN49VwuX3Bwu32hlVW2pMxqIpvDCKp47UVTpeX2kYn9Y1O6hpd1Dc6qW9y+b8amlz+Kbkn0oEmh291O7X2gMpjMqrExRhJjDdjMirEmY3ExRiIMxuJMRuINRuJNRuav1rOGfwtPZPp1Ht5RzIJEEJEsJa1Fq15vC0bNGn+XQEjaa/yGLOB7DTf1NmO6LpvTUmDzY21yYXV7qbR5sZqc9Fod7f5arJ72nVnncjt0XB7XO0G2zvDbFKJMRowmw3EGFX/mJLZqGI2qZiMLT/7vptafxlUjCf+3LwXeEcbLPWkoAaIDRs28Pvf/x5N07juuuu4/fbb2zzucrm499572b59OykpKTz22GMMGDAgmEUSImptO1jNp0XHqKyzk5ESx0Vn9WXskHSeLdjBVzsrcHs1TAaVCaMymDAqk8++O0ZlnYP0pFgm5mUyIieVnYdq+HJnBTUNDlITY5gwKpPSyka+3FGOzeUl3mxg4ugspp7T8et0b2kdm3dVUGt1+p+fOyDlpOc7c4+Pt5T6yuH0EB9jDKgcNQ0Oki1mRg1KxeH2sv1ADQ1NLmLMBjJS4og1+xIv2l0e3B4Nq82F3enF6Q4wo14zl9vXkuOE/cW7U8uOfAaDirFlq1fD8f3BW2//alBV/7avBoOKyaAwOjez078zaOm+vV4v06dP5/nnnycrK4trr72Wv/zlLwwfPtx/zauvvsru3bt56KGHKCgo4MMPP+Svf/3rae8dqemiAxHJ6bADEc31C2Xdth2s5o31B9qdt8Qa2Vlc2+acjq+fPiMlrs3580Zm8PXuyuPX6b6xgia7u3njHcWfOG7axBwuPzfHn2Za12F3SS1rN5W0K8OYIalsP1jb7vz0iTntgsTe0jre/7L9PbLT4vhuf3W785ed279dkOjoHg6XBwWIMbf9TNy6DK3Tfe8uqeX9TSX+LXVbdjQcOzSNFEsMDpcviLRMIHA0p1NpOdd6ckFLSpNQW/3n+Z1+TtBaEEVFRQwaNIicnBwAZs+eTWFhYZsAsW7dOpYtWwbA9OnTeeihh3z7G0dpf54QwfJp0bEOz+861P6NGcDewYZDn2w5QmK82X+sKApNdl9Kc0Vt28XxxfZyrr10eJvnv7F+v3/At/Xnzs27Kvz39Z/VoWh/FRNGHv9Uq+uwdV9Vh7OQOgoOAF/trGDW5MG0Hq3fsreSE99CWga+Y2PavuVt2VPJ2CG+KbuxZl/3D8DWvVW+T+An7JpubXIx/8IhHZblZHxpVjTcHi8uj2+syOXRfD/7v7y+caRW57ya7p+M4PF/6Xg8zT9rvp+9mt7mca3luLlL8UwELUCUl5eTnZ3tP87KyqKoqKjdNX37+lZtGo1GEhMTqa2tJS0t7ZT37srOSJFE6he5QlW32kZXh7NxNJ12b5bovjfqE69vcnhIO2HguKX/viU4tHy3OTzt6nqyMticXtKT22eHtdo9DByQ2u5crLn9tR5Nx2xsf97h8rbbCKfJ4SXuhECga4BCu/M2l5cRQ/v4j9OTfa2qJueOdte2XD9yWEa789EqaAGio56rE1sGgVwjhDi9v/zk0lAXoVvKEE33iAZBmwCcnZ1NWVmZ/7i8vJzMzMx21xw75msaezwerFYrKSkdD1wJIYToWUELEOPGjaO4uJiSkhJcLhcFBQXk5+e3uSY/P5+33noLgPfff5/JkydLC0IIIcJE0GYxAaxfv56HH34Yr9fLNddcw5133sny5csZO3Ysl19+OU6nk5///Ofs3LmT5ORkHnvsMf+gthBCiNAKaoAQQggRuSQJiRBCiA5JgBBCCNGhsM7F5HQ6WbRoES6Xy78y+6677qKkpIS7776b+vp6Ro8ezR//+EfMZvPpbxiGWsZnsrKyeOqpp6Kqbvn5+VgsFlRVxWAw8Oabb1JXV8dPf/pTjhw5Qv/+/fnrX/9KcnLy6W8WhhoaGvj1r3/Nnj17UBSFhx9+mCFDhkRF/Q4cOMBPf/pT/3FJSQl33XUXCxYsiIr6vfDCC7z++usoisKIESN45JFHqKioiJrX3osvvsjrr7+Orutcd9113HLLLV177elhTNM0vbGxUdd1XXe5XPq1116rb9myRb/rrrv0NWvW6Lqu67/5zW/0V199NZTFPCPPPfecfvfdd+u33367rut6VNVt6tSpenV1dZtzjz76qP7UU0/puq7rTz31lP7HP/4xFEXrFvfee6/+2muv6bqu606nU6+vr4+q+rXweDz6lClT9NLS0qioX1lZmT516lTdbrfruu57zb3xxhtR89rbvXu3Pnv2bN1ms+lut1tfvHixfvDgwS797cK6i0lRFCwWC+BbJ+HxeFAUhS+++ILp06cDcNVVV1FYWBjKYnZZWVkZn3zyCddeey3gWzgYLXU7mcLCQhYsWADAggUL+Oijj0Jcoq5pbGzkq6++8v/tzGYzSUlJUVO/1jZu3EhOTg79+/ePmvp5vV4cDgcejweHw0FGRkbUvPb279/P+PHjiYuLw2g0cv755/Phhx926W8X1gECfH/I+fPnM2XKFKZMmUJOTg5JSUkYjb7esezsbMrLy0Ncyq55+OGH+fnPf47avJVibW1t1NStxa233srVV1/Nv//9bwCqq6v9CyYzMzOpqakJZfG6rKSkhLS0NH7xi1+wYMECfvWrX2Gz2aKmfq0VFBQwZ84cIDr+fllZWSxZsoSpU6dy0UUXkZCQwJgxY6LmtTdixAg2b95MbW0tdrudDRs2UFZW1qW/XdgHCIPBwKpVq1i/fj1FRUUcONA+Y2UkLq77+OOPSUtLY+zYsae8LhLr1uKf//wnb731Fk8//TSvvvoqX331VaiL1G08Hg87duxg4cKFvP3228TFxbFixYpQF6vbuVwu1q1bx4wZM0JdlG5TX19PYWEhhYWF/Oc///G/iZ4oUl97w4YNY+nSpSxZsoSlS5cycuRIDIb2eawCEfYBokVSUhKTJk3i22+/paGhAY/Hl52xrKysXQqPSPDNN9+wbt068vPzufvuu/niiy/4/e9/HxV1a5GVlQVAeno606ZNo6ioiPT0dCoqKgCoqKg4bWLGcJWdnU12djbjx48HYMaMGezYsSNq6tdiw4YNjBkzhj59fAntoqF+n3/+OQMGDCAtLQ2TycSVV17Jli1bouq1d9111/HWW2/x6quvkpKSwqBBg7r0twvrAFFTU0NDQwMADoeDzz//nGHDhjFp0iTef/99AN566612KTwiwT333MOGDRtYt24df/nLX5g8eTJ//vOfo6JuADabjcbGRv/Pn332Gbm5ueTn5/P2228D8Pbbb3P55ZeHsphdlpGRQXZ2tr9Fu3HjRoYNGxY19WtRUFDA7Nmz/cfRUL9+/fqxdetW7HY7uq6zceNGhg8fHjWvPfB1BQIcPXqUDz74gDlz5nTpbxfWK6l37drF/fffj9frRdd1ZsyYwbJlyygpKeGnP/0p9fX15OXl8ac//Slip6MBbNq0ieeee84/zTUa6lZSUsJ//dd/Ab5xpDlz5nDnnXdSW1vLT37yE44dO0bfvn1Zvnx5xCZo3LlzJ7/61a9wu93k5OTwyCOPoGla1NTPbrdz2WWX8dFHH5GY6EvtHS1/v8cff5x3330Xo9FIXl4ev//97ykvL4+K1x7AjTfeSF1dHUajkV/84hdccMEFXfrbhXWAEEIIETph3cUkhBAidCRACCGE6JAECCGEEB2SACGEEKJDEiCEEEJ0KKyzuQpxKtdddx0ulwu3201xcTG5ubkAjB49mkceeSTEpQvM9u3bKSkpiaqVyiJ6yDRXEfFKS0u55ppr2LRpU6iL0o7H4/Hn9+nI66+/zueff85jjz3W7fcW4kzJvy4RlVauXMm//vUvvF4vSUlJ/Pa3v2Xw4MG8/vrrrF27FovFwp49e+jbty+//OUvefTRRykpKWH8+PE8+uijKIrCz372M+Li4jh8+DBlZWVMmjSJ3/zmN5hMJqxWKw8//DB79+7F6XQyZcoU7rvvPlRVZeHChUycOJEtW7YQHx/P448/7l8k6HQ6GT9+PL/97W9paGjg//7v/2hqamL+/PlMmjSJRYsWceONN/LZZ58BcOjQIf/xoUOHWLhwIddffz1ffPEFV199NfPnz+cvf/kLmzdvxuVykZeXx4MPPkhcXFyI/wIiKgQpJbkQPaakpESfOHGi//iLL77Q77jjDt3pdOq6ruuFhYX6okWLdF3X9ddee02fOHGiXlZWpuu6ri9ZskRfsGCBbrVadZfLpc+aNUv/4osvdF3X9XvuuUefP3++3tTUpLtcLv3mm2/W//GPf+i6ruv33Xefvnr1al3Xdd3r9ep33XWXvnLlSl3Xdf2GG27Qf/SjH+kej8f/eF1dnf/nu+++27+PxGuvvab/5Cc/8Ze9uLhYnzJlSofHxcXF+ogRI/S1a9f6H3/88cf9Of51XdcfeeQRffny5Wf2H1SIZtKCEFFn3bp17Nixg+uuuw7w7bPR1NTkf/y8887zJxIcPXo0DoeDhIQEAEaOHMnhw4eZNGkSALNmzSI+Ph7w5dD/5JNPWLhwIR9//DHbt2/n6aefBny5wgYOHOj/HXPnzvVn0NQ0jRUrVvDpp5+iaRp1dXVd3oUtPj7ev2dBS13tdjsFBQWAL/vqmDFjunRvIU4kAUJEHV3X+d73vseyZcs6fDwmJsb/s6qq7Y5bMnp2dN+WFNCapvHUU0/Rr1+/Dq9tCSoAq1atoqioiH/84x9YLBb+/ve/c+zYsQ6fZzAY0DTNf+x0Ok9635Yy/e53v+P888/v8H5CnAmZ5iqiTkvWypYNX7xeL9u2bevSvd577z3sdjtut5vVq1f7Wxb5+fmsWLECr9cL+DIPl5SUdHgPq9VKamoqFouF+vp6/6d9AIvFgtVq9R9nZmbicDj891qzZs1p6/rcc8/5A0ljYyP79+/vUl2FOJEECBF1Jk+ezLJly7jjjjuYN28ec+fO5ZNPPunSvc477zzuvPNO5syZQ05Ojn+L0d/85jdomsb8+fOZO3cut912G5WVlR3e46qrrqKuro45c+Zw9913t/m0f+GFF2K1Wpk3bx4PP/wwZrOZ+++/n8WLF3PTTTdhMplOWb4f/vCHDBs2jGuvvZa5c+eyaNEiDh482KW6CnEimeYqxEn87Gc/47zzzmPhwoWhLooQISEtCCGEEB2SFoQQQogOSQtCCCFEhyRACCGE6JAECCGEEB2SACGEEKJDEiCEEEJ06P8H8uC1Bmw5yFwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "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\"." + ] + } + ], + "metadata": { + "celltoolbar": "Hide code", + "kernelspec": { + "display_name": "bgdia703", + "language": "python", + "name": "bgdia703" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.17" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} -- 2.18.1