Iris_data_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b8783a8e",
   "metadata": {},
   "source": [
    "## Basic example with toy data\n",
    "This toy problem demonstrates the basics of classification modeling using logistic regression. For this simple example, we will use the well-known Iris data set (Fisher, 1936). We see that we have 150 records, four predictors (sepal length, sepal width, petal length, petal width), and three classes: iris setosa, iris virsicolor, and iris virginica."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "651cdf7c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   sepal_length  sepal_width  petal_length  petal_width        class\n",
      "0           5.1          3.5           1.4          0.2  Iris-setosa\n",
      "1           4.9          3.0           1.4          0.2  Iris-setosa\n",
      "2           4.7          3.2           1.3          0.2  Iris-setosa\n",
      "3           4.6          3.1           1.5          0.2  Iris-setosa\n",
      "4           5.0          3.6           1.4          0.2  Iris-setosa\n",
      "       sepal_length  sepal_width  petal_length  petal_width\n",
      "count    150.000000   150.000000    150.000000   150.000000\n",
      "mean       5.843333     3.054000      3.758667     1.198667\n",
      "std        0.828066     0.433594      1.764420     0.763161\n",
      "min        4.300000     2.000000      1.000000     0.100000\n",
      "25%        5.100000     2.800000      1.600000     0.300000\n",
      "50%        5.800000     3.000000      4.350000     1.300000\n",
      "75%        6.400000     3.300000      5.100000     1.800000\n",
      "max        7.900000     4.400000      6.900000     2.500000\n"
     ]
    }
   ],
   "source": [
    "#import required modules\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.tools as st\n",
    "from sklearn import preprocessing\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "#read data\n",
    "iris = pd.read_csv('iris.data', names=['sepal_length','sepal_width','petal_length','petal_width','class'])\n",
    "#print the first five lines in iris data set\n",
    "print(iris.head())\n",
    "#Generate descriptive statistics for the data set and print them\n",
    "print(iris.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4afaf1ed",
   "metadata": {},
   "source": [
    "Let's explore this data set by visualizing how the three species of iris can be associated with unique combinations of their sepal and petal characteristics (their features)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "516654be",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'petal_width')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGzCAYAAAAi6m1wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3c0lEQVR4nO3deVyU1f4H8M8zM+ybgAqIKG65a6aVmllKWVqmaZutWtcyrUwzzeVa3W5ZP6trltXVuppamrmUrWopqKnlgvuaohiiuCAg+8yc3x/HAQZm3xk+79drXsDMeeZ8n3kY5svznHO+ihBCgIiIiMiPqLwdABEREZGrMcEhIiIiv8MEh4iIiPwOExwiIiLyO0xwiIiIyO8wwSEiIiK/wwSHiIiI/A4THCIiIvI7THCIiIjI7zDBISIiIr+j8XYABjNmzMCUKVMwduxYzJo1y2Sb1NRU9OnTp8b9hw4dQps2bWzuS6/X48yZM4iIiICiKI6GTERERB4khEBBQQEaNWoElcryORqfSHC2b9+OuXPnolOnTja1P3LkCCIjIyt+btCggV39nTlzBklJSXZtQ0RERL7h9OnTaNy4scU2Xk9wrly5gkceeQTz5s3Dv//9b5u2adiwIerVq+dwnxEREQDkC1Q1USIiIiLflZ+fj6SkpIrPcUu8nuCMGTMGd911F2677TabE5wuXbqgpKQE7dq1w7Rp00xetqqqtLQUpaWlFT8XFBQAACIjI5ngEBER1TK2DC/xaoKzdOlS7Nq1C9u3b7epfUJCAubOnYuuXbuitLQUixYtQkpKClJTU9G7d2+z282YMQOvv/66q8ImIiIiH6cIIYQ3Oj59+jS6deuGtWvXonPnzgCAW2+9Fddee63ZQcamDBw4EIqiYPXq1WbbVD+DYzjFlZeXxzM4REREtUR+fj6ioqJs+vz22jTxnTt3IicnB127doVGo4FGo0FaWhpmz54NjUYDnU5n0/N0794dx44ds9gmKCio4nIUL0sRERH5P69dokpJScG+ffuM7hsxYgTatGmDSZMmQa1W2/Q86enpSEhIcEeIREREVEt5LcGJiIhAhw4djO4LCwtDbGxsxf2TJ09GVlYWFi5cCACYNWsWkpOT0b59e5SVlWHx4sVYsWIFVqxY4fH4iYiIyHd5fRaVJdnZ2cjMzKz4uaysDBMmTEBWVhZCQkLQvn17/PjjjxgwYIAXoyQiIiJf47VBxt5kzyAlIiIi8g21YpAxERERkbswwSEiQF8OlOUCeq23IyEicgmfHoNDRG6Wfww48BZw6itAXwZowoDmI4B2k4HQRt6OjojIYTyDQ1RX5e4GfukKnFwskxsA0BYCxz6V9xee8mp4RETOYIJDVBcJAWx5BNAVAaLaZSmhBUrPA9tHeyc2IiIXYIJDVBdd2ArkHQSEmRXDhQ448zNQmGn6cSIiH8cEh6guurzPehsImQQREdVCTHCI6iJNqG3t1CHujYOIyE2Y4BDVRQl3AoqVSZSBMUD97p6Jh4jIxZjgENVFwQ2Als8AUMy3afcKoA7yWEhERK7EBIeorrrufaDpMPm9ogEUtbwBQJuXgLYTvBcbEZGTuNAfUV2lDgRu+hJoN1GuhVOSA4QmAc2HAxEtvR0dEZFTmOAQ1XXRneWNiMiP8BIVERER+R0mOEREROR3mOAQERGR32GCQ0RERH6HCQ4RERH5HSY4RERE5HeY4BAREZHfYYJDREREfocJDhEREfkdJjhERETkd5jgEBERkd9hgkNERER+hwkOERER+R0mOEREROR3mOAQERGR32GCQ0RERH6HCQ6Rr9OVAWW5gF7n7UiIiGoNjbcDICIz8g4C+98EMr8BRDmgiQRa/gNo9woQ3MDb0RER+TSewSHyRee3Ar90AzKXyeQGALT5wJEPgDXXA8VnvRsfEZGPY4JD5GuEHtgyTF6aEtpqj+mAoixg5zjvxEZEVEswwSHyNdnrgMJTAMyMuRFa4PRyoOS8R8MiIqpNmOAQ+Zq8fYCittxGaIGCo56Jh4ioFmKCQ+Rr1KHyMpXVdiHuj4WIqJZigkPkaxLvtt4mJBGo19n9sRAR1VJMcIh8TVgTIPlRWHx7dpgGqKxcxiIiqsOY4BD5ohv+CzS+R36vaOSYHEUNQAE6TAdaPuPV8IiIfB0X+iPyRZoQoPcq4OIO4ORXQNlFIKwZ0Hw4EJ7s7eiIiHweExwiXxbbTd6IiMguvERFREREfocJDhEREfkdJjhEBkIAZXmAttjbkRARkZM4BodIVwYcmQUc/RAo+lveF9cXaD8FiE/xamhEROQYJjhUt+nKgNQBwLkNAKqsHpyTJu+7cR7Q4imvhUdERI7hJSqq247NAc6th1FyA8iq3RDAn6OA4mxvREZERE5ggkN125EPLT8u9MDxzz0TCxERuQwvUVHdpSsFCjOst7u83/2xEBGRS/EMDtVdhhIIFtuo5KrCRERUqzDBobpLpZaVuxULJzKFFmg82GMhERGRazDBobqt3SsABACl5mOKBohqBzS629NRERGRk5jgUN1Wvztw01JAFQRAuXrZ6uoZncg2QJ+18kwPERHVKhxkTNTkPrmg34kvgMt7AFUwkDgQaHSnHINDRES1DhMcIgAIjAbavOjtKIiIyEX47ykRERH5HSY4RERE5Hd4iYqI7CcEoL0CQAECwr0dDRFRDUxwiMh2QgAn5gOH3wPyDsr7orsAbV8Gmj4EKCam2xMReYHPXKKaMWMGFEXBiy++aLFdWloaunbtiuDgYDRv3hyffvqpZwIkquuEALaPAv54Csg7VHl/7h5gy8PAninei42IqBqfSHC2b9+OuXPnolOnThbbZWRkYMCAAbj55puRnp6OKVOm4IUXXsCKFSs8FClRHXbmZ+CvuVd/EFUeuFqJ/eDbwPmtno6KiMgkryc4V65cwSOPPIJ58+YhOjraYttPP/0UTZo0waxZs9C2bVv84x//wJNPPol3333XQ9ES1WHH5liu3aVoZBsiIh/g9QRnzJgxuOuuu3DbbbdZbbt161b069fP6L477rgDO3bsQHl5udntSktLkZ+fb3QjIjvl7gaEzvzjQivbEBH5AK8mOEuXLsWuXbswY8YMm9qfPXsWcXFxRvfFxcVBq9XiwoULZrebMWMGoqKiKm5JSUlOxU1UJ6lDXdOGiMgDvJbgnD59GmPHjsXixYsRHBxs83ZKtVkaQgiT91c1efJk5OXlVdxOnz7tWNBEdVmT+yxfooJKtiEi8gFemya+c+dO5OTkoGvXrhX36XQ6bNy4ER999BFKS0uhVhv/MY2Pj8fZs2eN7svJyYFGo0FsbKzZvoKCghAUFOTaHSCqa1qNBo5+BGiLUDGw2EBRAwFRQIunvBIaEVF1XktwUlJSsG/fPqP7RowYgTZt2mDSpEk1khsA6NGjB77//nuj+9auXYtu3bohICDArfES1XlhSbK6etpAoOxiZdV1oQWCGgB9fgGCzP+jQUTkSV5LcCIiItChQwej+8LCwhAbG1tx/+TJk5GVlYWFCxcCAEaNGoWPPvoI48ePx8iRI7F161Z8/vnnWLJkicfjJ6qTGvQABmcCp74Gzm+WC/vF9QWShgJqniUlIt/h0ysZZ2dnIzMzs+LnZs2a4aeffsK4ceMwZ84cNGrUCLNnz8bQoUO9GCVRHaMJBVqMkDciIh+lCMMo3TokPz8fUVFRyMvLQ2RkpLfDISIiIhvY8/nt9XVwiIiIiFyNCQ6RJ+n1QGEmUGJ+3SYiInKeT4/BIfIb2hJgyzAg6wc56wgANJFAm/FAp1e9GxsRkR9igkPkbtoSYHUyUHKu2v35wP7XgLy9wM0sGEtE5Eq8REXkbn8+UzO5qer0SuDses/FQ0RUBzDBIXK3zGXW2+yZ7P44iIjqECY4RO6mL7He5kqG++MgIqpDmOAQ+QJVoLcjICLyK0xwiNwtrLn1NokD3R8HEVEdwgSHyN2ufcfy44rGehsiIrILExwid2t6H9DhNdOPKQFA39+AQJYMISJyJa6DQ+QJnV4FWj4F7BoPXEoHVBqg8b1Ah+mAJtjb0RER+R0mOESeEtoY6GXDlHEiInIaL1ERERGR32GCQ0RERH6HCQ75r+LzQO5eQFfm7Uh8n74cKMsF9FpvR0JEbiIEUFAAFBa6t58rV+TN25jgkP/Z809gaRCwqiHwc2fg6yDg2yZA3lFvR+Z78o8BW0cAy8KB5THA8nrAjueBojPejoyIXESvB+bNA9q1AyIjgfBwoFs3YJkLhwQKASxcCHTqBEREyFvnzsCiRfIxb1CE8FbX3pOfn4+oqCjk5eUhMpLTc/3K5oeAzK/NPKgAdx0Aotp6NCSflbsbWNcb0BUDosqZG0UDBNUH7tgGhDX1WnhE5DwhgKeeAubPBxSlMtlQqWTiM20a8MYbzvcxdizw4YeVz1u1jxdeAGbNkv07y57Pb57BIf9xJdNCcgMAAvitj8fC8WlCAFseAXRFxskNIH8uPQ9sH+2d2IjIZVavlskNYHwmxZCE/PvfwI4dzvWxYYNMbqo+b9XvZ88GUlOd68MRTHDIf2x91HqbknNybE5dd2ErkHcQEDrTjwsdcOZnoDDTs3ERkUt99BGgVpt/XKMBPv7YuT7mzJHP484+HMEEh/xH/hHb2p1b5944aoPL+2xoJGQSRES11p49gM7M/zEAoNUC6enO9ZGeLp/HUh+7djnXhyOY4JD/sLUid3C8e+OoDTShtrVTh7g3DiJyq1Ab3uphYc71Ycv24eHO9eEIJjjkP1r8w4ZGChDf1+2h+LyEO+VgYksCY4D63T0TDxG5xX33Wb5EpSiyjbN9qCxkEyqV8304ggkO+Y9Or1r/0E4a4plYfF1wA6DlMwAsTGto9wqgDvJYSETkes89BwQFmU5A1Gqgfn1g+HDn+njmGTkt3FQipVbLqelPP+1cH45ggkP+5Y7tAMz8uxJ9HXDzco+G49Ouex9oOkx+r2gARS1vANDmJaDtBO/FRkQukZwMrFkDREXJnzWaygHB8fHA+vVAvXrO9REfD/z6KxATU7OPmBhg3TogLs65PhzBdXC4Do7/0emA9PHAyS8BfRkQHAfc8F9emjIndw9wcjFQkgOEJgHNhwMRLb0dFRG5UGEhsHQpsGWLPJuTkgIMGQIE2jh00RYlJXLxwI0b5c+33ALcfz8QHOy6Puz5/GaCwwSHiIioVuBCf0RERFSnMcEhIiIiv8MEh/yXrhQouwwIvdWmjvdRJqtw29qHEEBZHqAtcl9MREQEK3NqiWqhC38CB94Esn4AoJeFI1s9C7R9GQiIcE0fubuB/W8Bf6+UZQ0C6gGtngHaTQICo2u212uBox8BRz4ACk/K+xr0AtpNBhIHuCYmIiKqwEHGHGTsX7J+ADbeC0AY11lSVEBUB+D2TUCAk8f87Hogtb98fqM+1EB4C6DfFiAotvJ+vRbYNBTI+v7qHYa3nBqADuj6AdD6BediIiKqAzjImOombZGskF098QDkJaS8A8C+fznXh74c+P0hmbTU6EMHXDkO7JlifP+JBUDWasjEpur/E1e33zkOKDjuXFxERGSECQ75j8xlQHk+jJOIKoQO+GueHJvjqL9XA6XnAZgZcyN0wImFQHlB5X1HP4TFt5qiAH/NdTwmIiKqgQkO+Y/cvYAqwHIbbT5QnOV4H5f3AoqVPvQlwJWMyp/zDsJsQgTIpMim6t5ERGQrJjjkPzShcpaSNWobK2mb68NSslLRR5Uq3Gpry3iqbK/uTURENmGCQ/6j8WBAaC00UAExXYGQeMf7SBxUc+yNEQWIuMa41EHje60UAdXL2ImIyGWY4JD/iO0GxPWtLBhZgx7o8E/n+ohqI5MRs30IoMN0Oa7GoO0EyKrdJip3KxogrBnQ5H7n4iIiIiNMcMi/3LwcqN9Tfq9orp45Ucmv3eYAjQc530ePRUBcSpU+DFW4VcC1/wc0e8S4fXQnoPe3Vy+NKVXiAhDeDEj5DVAHOR8XERFV4Do4XAfH/wgBnN8EZH4jZzNFtgaaj3Du0pSpPi5sAzK/lqslh7cAWowAQhub36Y8H8hYBFzaAagCgUYDgEZ3ASqut0lEZAtWE7eCCQ4REVHtw4X+iIiIqE5jgkNERER+hwmOvxN6We3amdV7fYHQy7EuuhJvR0JEtVxpKZCbC+htWNKKai8mOP6q7DKwezKwvD6wPAZYFgZsuh+4lO7tyOyjLQb2vwmsSgSWRwNfhwIb+gM5m7wdGRHVMrt2AUOHAqGhQEwMUL8+MGUKkJfn7cjIHTjI2B8HGZdeAtb2BK78Va3atUZW1b71JyA+xXvx2UpbDKxPAS78AaPVgxW1nMXU80sg+SGvhUdEtcfatcDdd8uzNroqfxbVauCaa4Dffweio70XH9mGg4zruj1TaiY3gFzlV6+V1bB1Zd6JzR6HZtZMboCr+yWAbSPk5TciIgtKS4FhwwCt1ji5AeTPR48C/3RyDVDyPUxw/E15AXDiCwvlBPRA6QUg6zuPhmU3oQeOzoH5uk8C0JfKyt1ERBasXAlcumS+VJ1OB8yfDxQWejYuci8mOP7mSoasZm2JEgBc3u+ZeBxVegkozbHcRlEDeT6+H0Tkdfv2AQEBltsUFQGZmZ6JhzyDCY6/sakqtd642rUvslqB29DOx/eDiLwuNNS2GVMh/HPiV5jg+JvwFkBEa5gs7GggdK6pyeROAeFAXB8LRS0hxxSxCjcRWTFoUM2xN1UpCtC+PdC0qediIvdjguNvFAXoOB2AmYvNihpIvAeIauvRsBzSfqoci2OKogFirpdJEBGRBR07AgMGyBlTpggBTJ8u/3yS/2CC44+SHwa6zISson210rWhenVcX6DnYq+GZ7P4FKD7AkAVgIqK4Ib9iO4M3PoD/yIRkU2WLAFuuUV+r9HIm1oNqFTA++8DDzzg3fjI9bgOjj+ug2NQ9DdwYgFQ8BcQEAU0fQio3732JQUl54GML4C8A4A6FEi6VyZqCvNzIrKdEMCWLcDXXwP5+UCrVsDw4UBiorcjI1uxmrgVdSbBISIi8iNc6I+IiIjqNCY4RERE5HeY4FDtoNcCV07JIqK2KsuX2+i1trUXAijLkzWwfImuTJak0FuY50pEREY03g6AyKKSS8DvDwLn1qOibENQfaDDq0Dr50xvc3w+sGcyUHLu6h0K0KAXcNNSILRRzfa6MuDILODoh3JgNiAHMbef4t2ipHkHZSX1zG8AUQ5oIoGW/wDavQIEN/BeXEREtYBXz+B88skn6NSpEyIjIxEZGYkePXrg559/Nts+NTUViqLUuB0+fNiDUZPHlFwCVjcFzv0Ko5pUpReAnc8DO16ouc3e6cAfT1ZJbgBAAOc3AatbVCYwBroyIHUAsHuy8WM5acD624Hjn7tyj2x3fivwSzcgc5lMbgBAmw8c+QBYcz1QfNY7cRER1RJeTXAaN26Mt99+Gzt27MCOHTvQt29fDBo0CAcOHLC43ZEjR5CdnV1xa9WqlYciJo/afB+gvWL+8aMfytpbBiUXgP1vmG+vLwE2Dja+79gc47NDBoaK5X+OAoqz7QzcSUIPbBkmky9R7fKa0AFFWcDOcZ6NiYiolvFqgjNw4EAMGDAA11xzDa655hq8+eabCA8Px7Zt2yxu17BhQ8THx1fc1OaWp6TaS68FclKtt9v1UuX3u1+x3v7STjk2x+DIh5bbC73nz+JkrwMKTwEwM+ZGaIHTy+X6QEREZJLPDDLW6XRYunQpCgsL0aNHD4ttu3TpgoSEBKSkpGDDhg1Wn7u0tBT5+flGN/JxBcdhttxEVXkHK7/PTbftuXN3ya+6UqAww3o/nq68nrfPcg0uQCY5BUc9Ew8RUS3k9QRn3759CA8PR1BQEEaNGoVVq1ahXbt2JtsmJCRg7ty5WLFiBVauXInWrVsjJSUFGzdutNjHjBkzEBUVVXFLSkpyx66QKwVG2dZOHVTlextLAQfGyK+KxnoioagAjYdLDKtDzdfgMmrH0sdEROZ4fSXjsrIyZGZm4vLly1ixYgU+++wzpKWlmU1yqhs4cCAURcHq1avNtiktLUVpaWnFz/n5+UhKSuJKxr7um2ig/LLlNte+C7S7epnqxEJg2xOW26tDgQcLK3/eOBjI+rHmWJeqen/r2errhZnAd8mweGYpJBEYdApQ8fIsEdUdtWol48DAQLRs2RLdunXDjBkz0LlzZ3zwwQc2b9+9e3ccO3bMYpugoKCKmVqGG9UC7adYflwdCrSpMti2+eNAQLTlbVqPNf653SuQiYSJ+lyKBohqBzS625ZoXSesCZD8KCy+PTtMY3JDRGSB1xOc6oQQRmdbrElPT0dCQoIbIyKvafcy0GKk6cfUocAdf8pSwFX13ynXizGlyYPAtW8Z31e/u1wfRxUEQDGuWB7ZBuiz1juJxA3/BRrfI783XEpT1DLGDtOBls94PiYiolrE4YX+9Ho9/vrrL+Tk5ECvNx4v0Lt3b5ueY8qUKejfvz+SkpJQUFCApUuXIjU1Fb/88gsAYPLkycjKysLChQsBALNmzUJycjLat2+PsrIyLF68GCtWrMCKFSsc3Q3ydTfOlWdp0l8G8o8A6mCg2WNAm/GAysSvb3gz4L5cOf37+DxAWwSEtwC6zASiO5nuo8l9ckG/E18Al/cAqmAgcSDQ6E7vVSzXhAC9VwEXdwAnvwLKLgJhzYDmw4HwZO/ERERUiziU4Gzbtg0PP/wwTp06hepDeBRFgU5n25Ly586dw2OPPYbs7GxERUWhU6dO+OWXX3D77bcDALKzs5GZmVnRvqysDBMmTEBWVhZCQkLQvn17/PjjjxgwYIAju0G1RVRb4NYfbG+vUgGtn5c3WwVGA21etDs0t4vtJm9ERGQXhwYZX3vttbjmmmvw+uuvIyEhAYpiPH4hKsrGGTBeYs8gJSIiIvIN9nx+O3QG59ixY1i+fDlatmzpUIBERERE7uTQAIMbb7wRf/31l6tjIXcQelmJWmf7wG27VVT6duMCio5UE9eVyva2rCkDeKaauF4vp4GXXHBfH/7kyhWgoMCtXRQWAvn58vDbQqcDcnOB8nK3hkVETrL5DM7evXsrvn/++efx0ksv4ezZs+jYsSMCAgKM2nbqZGYwJ3lO2WXg4DvAsf8C5blyBk7je+XU65gurumj6Azw+zBZyNKwZktwHNDpTaDlU67pw5Fq4hf+BA68CWT9ILcJqg+0ehZo+zIQEFGzvSeqiWtLZH2prB8q19zRRMrB0p1edU0f/kII4KuvgJkzgT175H0dOgAvvQQ88QSgmJjS74Dly4H/+z9g+3b58zXXAOPHAyNH1pycBwA5OcCMGcDnn8ucKyAAeOABYOpUoG1bl4RERC5k8xgclUoFRVFqDCqueKKrj9kzyNhb/H4MTuklYG1P4MpfV4tGXqVo5KygW39y/oO76Iyszq0vMf14+2lAZwuFL21hqCZuruDmNc8D3WYb35f1A7DxXgCi2r6rgKgOwO2bgIAqx9xQTfzcBhgV3FTU8szPjfOAFk4ma9oSYHVytQrnVSQNAW7mTMAKL78MvPuuzDIMMzQVRSY+o0cDH33kdJLzr38Br75quovHHgMWLDBOcs6cAbp3l1+r/nnTaIDAQGDDBuCGG5wKiYhsYM/nt80JzqlTp2wOoGnTpja39Qa/T3D+HAUc/8z4A76CCgiKAQZnAepAx/v45Xrg0g7LbYacA4IbOt7Hr32BHCu1xu45IaeGA3JK+KoEoLwAJlcBVtRA6xeB696tvO/wf64W7DTzNlA0wOBMIMSJtZa2PAGcXGi5Td/fgPi+jvfhLzZtAqwtM7FmDdCvn8NdpKcD111nuc3y5cDQoZU/338/8O23gNbEgtdqNdC0KXDsmOkzP0TkOm5Zybhp06YVt1OnTiExMdHovqZNmyIxMdGuRIjcoLxArudiMrkBAD1QegHI+s7xPsryrSc3AJA+yfE+HKkmnrkMKM+H2WRF6IC/5hmPR/JENfHMZdbb7JnsXB/+4uOP5WkRczQaeQbHCZ98YrkLtdq4i7NngZUrTSc3gDyjc+KEPItDRL7Dof83+vTpg0uXLtW4Py8vD3369HE6KHLClQzzl40MlADnKmTn7rat3eW91tuY41A18b2AKsB8WwDQ5gPFWfJ7T1UTt3Y8AHncSJ5eMZdJAPKx3bud6mL3bstd6HSVQ38A4PDhystY5qhUwH4PF50nIsscSnAMY22qu3jxIsLCwpwOipygCbWhkd65StSB9Wxrpw52og8HqolrQm2bCqO++hr5UjVxlROXC/2JLX8/Qm35HbfchbUhPCFVDrkt3QlhvA0ReZ9d6+AMGTIEgBxQPHz4cAQFVX646HQ67N27Fz179nRthGSf8BZARGug4CgsXqpxpjp2dCeZJOiKLLdr+bTjfYTEAwH1rFcTT3688vvGg+XsKbNUcgZZSPzVH9VA4t2Wq4kLrXxeZ4Q1BwpPWG6TONC5PvzF0KHyFIu5UyZqtRwQ44QhQ4C0NPOPazRydpTBddcBcXHAOTNjxAGZMN11l1NhEZGL2XUGJyoqClFRURBCICIiouLnqKgoxMfH4+mnn8bixYvdFSvZQlGAjtNhftCsGki8R5Y/cIa1sgYB9YDmTzjXh73VxGO7yendZs/K6IEO/zS+yxPVxK99x/LjisZ6m7pi5EggKkomMtWpVPL0y7PPOtXF448DDRua7yIgAHiuygoEGg0wbZr551OpgOHDgcREp8IiIhez6wzO/PnzAQDJycmYMGECL0f5quSHgeIzcpBvxbl4RZ6NiOsL9HRBEtr5TeDKCeDU0pqPaSKA/ruc76Pdy0DBMVk0szpz1cRvXg6kDZJr8xiqggu9vNTU9YOaZ64M1cS3PAboSyuTI6G9Wk38F+eriTe9D8h7Ddj/Ws3HlACg769AoB/O5nNEgwbAb78B/fvLUyaG0cBaLRATA/zwA9CokVNdREUB69cDd9wB/P23cRcREcB33wEtWhhvM2aMHGz81luVv3KKIre5915gzhynQiIiN3CoFlVt5/fTxA2K/gZOLAAK/gICooCmD8kPdBctlAZADuxNnyjX3NGEAs2fkuvTuHK+bN4h26uJA3JAxPlNQOY3clZZZGug+YjKS1OmlOW6v5p40d/ArvHApXQZe+N7gQ7TAY0TY5X8VWkp8M03QGqqPJ69e8vrRi4c6FJWBqxaBfz6qxxY3LMnMGyY5WFAGRlyjZyTJ4H69YGHHwa6dnVZSERkhVvWwenSpYvJgcWm7Nrlgv/e3ajOJDhERER+xC3FNgcPHlzxfUlJCT7++GO0a9cOPXr0AABs27YNBw4cwOjRox2LmoiIiMhFbE5wXn21sl7OP/7xD7zwwgt44403arQ5ffq066IjIiIicoBDY3CioqKwY8cOtGrVyuj+Y8eOoVu3bsjLy3NZgO7AS1S1kNDLVYrVwbavr6MrBXTFsvaUK8fSUJ2mK9MhP6sA4XFhCAi1srCkD9NqgawsIDoa4J9Bqi3cUqqhqpCQEGzevLnG/Zs3b0ZwMAdMkgtpi4H9bwKrEoHl0cDXocCG/kDOJvPbXPhTzqT6OlRuszIO2Dv9ao0qIsfkHDiPcdelITqoCDHN6yEsTODRZptx6Ifj3g7NLn//DfTqJYuEJifLWWUJCcBCK+XSiGobu6aJG7z44ot49tlnsXPnTnTv3h2AHIPzv//9D9OnT3dpgFSHaYuB9SnAhT9QWelbAGfXAdlrgZ5fAskPGW9TtZq4YZvSC3IBwL+/q1lNnMgGZ3adRfcbdDijuwm6q382yxGIr0/eiFUDy7Dhfwdww4j2Xo7SusxM4Jpr5CS1qs6eBZ54Qs4O459w8hcOTxNftmwZPvjgAxw6dAgA0LZtW4wdOxYPVF0C1EfxElUtse9fwL7XUZncVKUAqiBgyBkgMFre5Ug1cSIb3N94K77N6gYtal6SUkOLpposHCtOgkrj25dCu3SxXMpLUYBLl4B69TwVEZF93DJN3J8wwakFhB5YmQCU5lhopADX/QdoM1b+eGIBsG2E5efVRAJDc4xrWBFZcHZvDhI7x0IPyws+/vp/u5Dy8nUeisp+ly4BsbHW2z3zDPDpp+6Ph8gRbh+DQ+R2pZesJDeQZ2TyqpRwtreaOJENDq8/YzW5UUGH/Vt9e4zX9u22tUtPd28cRJ5i8xicmJgYHD16FPXr10d0dLTFRf8uXbrkkuCoDrN1plTVquj2VhMnskFolPWZUgIKQkJduEK4G8TE2NaOVdHJX9ic4PznP/9BRERExfe2rmpM5JCAcCCuD5CzUVY/N6V6pW97q4kT2eC6Ya0R948cnNM3gMmirAAUCNz1YiuTj/mKrl2B4GCgpMRyOydrmRL5DI7B4Rgc33X2N2D97TA9YFgDRHcB7vjDuLbWbylATpr5pKj3tzULbhJZ8dH9aXh++S0mH1NBh+GttuDzozd7OCr7vfwy8K6FMfaxscCFC56Lh8hebh+D88gjj2DevHk4evSoQwES2SQ+Bei+4Oq4GpVMagwVwqM7A7f+ULNw6M3Lgfo95fcV7a9u220OkxtyyJive2PqTalQoIcaWqhRDg3KAQD3Jv6JOX9e7+UIbTNzJjB0qOnH6tUDfLyMIJFdHDqD88wzzyAtLQ1Hjx5FfHw8brnlFtxyyy249dZb0aZNG3fE6VI8g1PLlJwHMr4A8g7I8TNJ9wJxfc2vTuxINXEiG2RsPI0F/zyOk39rUD9ai4dfjEPXR9t6Oyy77d4NTJwInDghq6ePHAmMHg2oOO2EfJzHpomfPXsWqampSE1NrUh4GjZsiOzsbEef0iOY4BAREdU+HpsmHhERgejoaERHR6NevXrQaDSIj+d/yURERORdDiU4kyZNQvfu3VG/fn1MmzYNZWVlmDx5Ms6dO4d0LqJAREREXubQJSqVSoUGDRpg3LhxGDRoENq2rV3XoH3mEpVeJxeeU4cB6kD39CH0QHmeHLti6+q9JZcA3RUgpLHvXJR3pJo4+RatFigoACIiAI1DZfBqraILRSgvKkdk40goKutLbOh0QH6+HB8T6KY/DY4oLpZ1rKKiao7vN0Wvl/sRHCxvNikpkbfISLf9/SktBYqK5H74yp84so3bL1Glp6dj6tSp+PPPP9G7d2/Ex8fjwQcfxCeffFJRm4osKDkP7HoJWB4jb8vCgd8fBi4fcF0fZZeB3ZOB5fWv9hEGbLofuGThDNvBmcA3McDKWOC7psDSAFm5uyzfdXHZy5Fq4uRbTp4Enn5aJjYxMfKDa9QoWfnRz33/zz/RK3IvwhqEol7TKLQMOo3ZQ9OgLdGabH/pEjBpElC/vnypwsKABx8E9u71cODVpKUBd94p44mOBhITgbfeMr+mTmEh8K9/ySrl0dFAaChw113Ali0WOtm6FRg4sLKT+HjgtdeAK1dcth+7dslZZKGh8vWtXx+YMgXIy3NZF+RDXLIOzp49ezBr1iwsXrwYer0eOp2ZNUh8hFfP4BSfBdb2AIpOG6/VomjkdOi+vwINejrXR+klYG1P4MpfNftQVMCtP8kp2FVtexI4Md/08wVEA4NOAoEefq1MVhOHLNEghOlq4uRbDh8GevaUZ260VT7UNRr57/OWLbK8tR96f1AqXlp9K1TQQn91TVXl6u/x4IQ/8c3J66EOrCwBceEC0KMHkJEhz+AYaDSAWg2sWQPcYnopHrf68kvgscfkmY6qcalU8tCuW2d8dqawEOjbF9ixQ57BMVBffdsuW2Ziqvq33wL33Se/r97JtdfKDCs83Kn9WLsWuPtuGVPVLtRq+Sv4++8yryLf5pFBxunp6fjPf/6DQYMGoU+fPli0aBE6d+6M8ePHO/qUdcOu8UDR3zUXohNaQFcK/D5MXo5xxp4pNZMbQx96LfD7Q4CurPL+3L3mkxsAKM8Ftj7mXEyOODSzZnIDXN0vIQtrluV6Pi6y3YgR8hqFttoZC60WuHwZeOopr4TlbkfXZGDC6t4AUJHcAICACgIqrMq+EYtGbzXaZuJEebKr+v+HWi1QXg489FDNl9HdLl4EnnxSJibV49LrZX76/vvG97/9NrBzp3FyA8jthQAef1z+SlQoKAAefbRm5mHoZM8e4E1LK5RbV1oKDBsmX7/qXeh0wNGjwD//6VQX5IMcSnCio6Nxww034Msvv0SrVq2wcOFCXLp0CTt27MDMmTNdHaP/KLkg12YR5v5K6YGiTCB7reN9lBcAJ74wv5Iv9EDpBSDru8q70idYf94zPzkekyOEHjg6BzWSm8oGgL4UOLHQk1GRPfbtA7Ztq/mJYqDTAZs3A354WXvu1FNQw/yZbBX0+PCrytLely/LMyXmEhi9Hjh7FvjhBxcHasUXX1hOqvR64KOPKkvAabXAxx+bP+RCyHE8ixdXuXPJEjkgxtzFBJ1OljcvL3doHwBg5Up5+c9SF/Pny7NP5D8cSnAWLVqEixcvYseOHXj33Xdx9913mzxV9Pfff0NfPY2vywqOWEhurqpeIdteVzIAvZViM0oAcLlKH/lHrD+v0MpxPZ7iSDVx8i379tnWbr//HcM9GRHQwnyRTj3UOFDcrOLn48eBsjKzzQHIS1Wefqn27rU+CDc7u/KMzPnzMpGwpMZ+7NtnfdD55cvAuXPWwjVr3z4gwErN1KKiOjEsrE5xKMExl9BU165dO5w8edKRLvyTLVWshd64Qra9NLZUyq7Wh62zq1QerMLtSDVx8i2hNv6++GH56rDAcigWzuAAQBBKK7635aXS6z3/UtkSl6IAQVf/hNgSnxDV2oWGmj+1UpUTOx8aWvOSmYu7IB/k1glydbCOp2X1OgEhidbbJd7teB/hLYCI1jBX9RiAvHxVtSZTk/utP29wPKDx4HxVQzVxRW2+TfVq4uRb+va1/okRFgbceqtHwvGkwXdpISz8edWgHPe12lPxc+vWQPPmlp9TrwfuucdVEdpm8GDLl6jUaqBfv8pBxvXqAb16WT7ro9XK560waJDlTlQqOfo6NtZ8GysGDTJ/2QyQSVr79kDTpg53QT6IKwB4kkoNdLA0kk0FJD8ChDnxLlMUoON0mKzADciEIfEeIKrK2kXtXwVUVs7idHzd8Zgc1X6q+QHXigaIuV4mQeSbIiOBF180v2CKogDjxzs9O8YXPfh/XdFEnQU1ao4bUaCDCnqMn9mo4j6VCpg+3fzzqVRyklGrVu6I1rzbbgO6dDF/BUmvByZPNr5v2jTzZ0vUajnzqlevKnf26AHcfLN80Fwn06bZHXtVHTsCAwaY70II+frbsrYP1R5McDyt5dNAx9cAKDLZUNSVFbITBwI3zHW+j+SHgS4zIatoV+sjri/Qc7Fxe00gcPtm80lOm3FAq6edj8tejlQTJ9/yxhvAP/4hvzfMdzZ8Wj7zDPDqq96LzY1CYkKwfq0OyZozAOQZGw3KoUCPUBTjuzf2of2glkbbPPEEMGOGTGbUauOX6o47gAULPLwTkLH89JNMEAAZj0Yj7w8MBBYtqjl1/Y47gM8/l2NeVKrKbQCgWzfgu++qvW0VRU4Tv/HGmp1oNMB//yuzEyctWVIZq6ELtVp28/77wAMPON0F+RiXrINjTkREBPbs2YPm1s69ephPrGR85SRwYgFQmAEExgLJw4DY613bR9Hfso+Cv4CAKKDpQ0D97uaTAr0WOPAWcOprQF8G1OsAdHkXiGjh2rjsZW81cfI9Bw8CCxfKEamNGsm5wrVsBXRHaEu0+P7VHfhldRnKdQqu7yrw6HtdENEowuw2p0/LGT0nTshLPsOGATfc4N1cXq8Hfv1V5iHFxUCHDsDw4ZavGuXkyKTs4EF5JXLoUKBPHwv7IYRc72b5crm4X9u2spO4OJfthxByavvXX8uB0a1ayS4SbRg5QL7BY9XErYmMjMTu3buZ4BAREZHTPFZN3BoOMiYiIiJvcGvFu4MHD6JRo0bWGxIRERG5kM0JzpAhQ2x+0pUrVwIAkpKS7I+IXMuRauIO9ZEv16SxtY/is0BJDhDZxn2V1Il8RFGRXIg3MtJ3xsXry7TIP3wGIQn1ENTAPZfqtWV6ZB0uQFTDINSLt7WcOJFr2HyJKioqyuYb+QBHqonbq7wA2PsqsDKustJ32j1X60eZsXM8sCQQWJUA/NwZ+DoI+K4FcIVLiJL/+f57OSU6LEwOGG7ZEpg92/M1parKP5SFaUlfoGHQZUR3boLQhmEYHPwzdrz+o8v6yMkoRN9mGQgMApI7RyE6IQgNAy/jvxOPu6wPImvcOsjYV/n9IGNHqonbq7wAWHezLJVg1IcagALcvBJoPNB4mw39gexfzDyhCrgnAwhv4lxcRD7i/feBl16S05AN68IYzt4MHgx88435dVncJW//afTqlIdDog10VU7gq1EOBcDqR79B/0UPO9VHTkYhmrdUUKgPgfGCo/KjZsI9RzHzu9ZO9UF1l88MMiYvsbeauCP2/6tmcgPIn4UO2PIIoK1SuS53r4XkBgD0wG+3WHicqPY4ehSYcLWGbdVF74SQt1Wr5BoynvZqn401khsA0CEAOqjw8OL+KDl72ak+ht5y3kRyg4qf3119Df4+lF9jOyJXczjBWb58OR544AF0794d1113ndGNvMiRauL20pUCx+Za6EMA2gK5no7B1setP2/hSaCs2PG4iHzE3LmWz86oVMCHH3ouHgAozrqEzy4MqpHcGAiocRnR+ObhVQ73UXJFi99PN4H5UjHy/gmPOl44k8hWDiU4s2fPxogRI9CwYUOkp6fjhhtuQGxsLE6cOIH+/fu7OkayhyPVxO1VfAbQWvkPrHofhSdte+7LuxwOi8hX7NljeZyNXg8cOOC5eADg9C8HUAjLZTECUIb9hx2/bnZ4a67FGlwG+zPCHO6DyFYOJTgff/wx5s6di48++giBgYGYOHEi1q1bhxdeeAF5eXmujpHs4Ug1cXvZtK0ANFXaqQJse+6gBg6FRORLwsKsz5YKctOkRnNC48yvnmygh4KQQBvKbpsRHWfbjMjgAMuV1olcwaEEJzMzEz179gQAhISEoKCgAADw2GOPYcmSJa6LjuznSDVxe4XEAzHdYPHXp3ql7yQbKpYrGiDqGsfjIvIRgwfLsTbmaDSyeKYnJfbvhE7KXqhgPrnQIQCDn01wuI+mnaIQqb4Cs8V+AQAKhg8rdbgPIls5lODEx8fj4sWLAICmTZti27ZtAICMjAyuXuxtjlQTd0SHfwIwV+lbDTTsY1xbq+uHsPrr1ny4czER+YgHHwSaNDE9DkdR5Bic8eM9G5OiVmH6wN3Qw/QlKDW0uCPgN1w76Q6n+pn46BnIf7BM/Q0SiFAVYvR7Xq5vR3WCQwlO37598f333wMAnnrqKYwbNw633347HnzwQdx7770uDZAcYG81cUc0vge4/uOrz1ut0nf9HkDvFcbt1Wrgto0we2apwc3AjfOcj4vIB4SEAOvXA8nJ8mdD9WpFAUJDZUXt9u09H9fQ7x7H7M6fQw0tVNBVVDkHgJvVW/H1LufPoE5dcA2G33T06k+iyg0IUxVj+++lUKl9ZLVD8msOrYOj1+uh1+uh0cgPtGXLlmHz5s1o2bIlRo0ahcBA316Z1u/XwTGwt5q4I4rPyj7yDwMBEUDSfUDD3ub70JUBO8YAp1cC+nIgNBG44TOg4U2ui4nIR2i1crG/X36RKxlffz3w6KNAhPXhMG6VvW4/5o/+E0fPRSEysBT3PxqEXu/dC0XtupVDDmy8gJefzsPRM+EICdBi+P1FGDenJZMbcorPVBP3VXUmwSEiIvIj9nx+O1xsMzc3F59//jkOHToERVHQtm1bjBgxAjExMY4+JREREZFLOHQ+Mi0tDc2aNcPs2bORm5uLS5cuYfbs2WjWrBnS0tJcHSMRERGRXRxKcMaMGYMHHngAGRkZWLlyJVauXIkTJ07goYcewpgxY1wdY+2hLZRVtW296qfXAWW5zpdN8Da9FrhyCiizY/l1XYksCCpsXHND6GV7nZVFDH2dEEBeniwv7S56PfD330BOju3blJcDublurQKpLdEiN+MytCW293HlCnB1FQq7+igvKncgQt9RXAzs3g1cumTnRnl5Nv/90Wv1uHwqD6X5tk/ZLi4GLl+2/U+cJwghYyq2YxH00lL56653fMkfn2B42+psXFZICPl+Kiy03tYfOJTgHD9+HC+99BLUVeZAqtVqjB8/HseP18FqsadXAWt6AMvCgW+igB9aA0c/Nv/hXXIe2PWSrPK9PEZu9/vDwGUPL23qrKIzwLpbgKWBwOpkYHkUsDIe+Otz89vkbAY2DJCVx5dHA6saAfv/DWjNfOBri4H9bwKrEisrlm/oD+RscssuuY1WC8yaBTRvLstKh4UBN98M/PSTa/sYNgwIDgaSkoC4OCA8HJg40fw2x44BI0bIdjExMrbnnwfOnHFZWCc3/42n225EREg5YprXQ2RIOUa124jMrVkm2wsBfPklcO21cjBuZCTQsSOwYIH5D9asHdl4rlMaokJKEdO8HiLCdHjymk04vv6Uy/bDE3buBBo1kjOtunQBYmPl9++9Z2GjtWuBPn1kw3r1gKZNgXffBcpM/+OU/3c+pvVKRcPAXEQnRyE0SoPBCduwY+FBs12kpQF33il/baOjgcRE4K23gBIv/r9RUgK8/TbQuLGMKTQUuP12YMMG89vs2gUMHSrbxsQA9esDU6bIvLA2OXIEePxxeTwMb9sXXwTOmamAodcD8+YB7drJ91N4ONCtG7BsmSej9gLhgJ49e4pVq1bVuH/VqlWie/fuNj/Pxx9/LDp27CgiIiJERESE6N69u/jpp58sbpOamiquu+46ERQUJJo1ayY++eQTe8MXeXl5AoDIy8uze9sa9r8lxJcQ4iuV/PolhPhSkbfNDwmh1xm3L8oW4ttkIb5SV2kPIb7SCLE0RIic352PyRMKs4RYEmy8D1Vvu6fV3ObkUiG+VNXc9y9VQvzSXYjyQuP25UVCrOkhHzd6rdTyvowlntlXZ5WXC3HPPUIoirwZai6q1fLrBx+4po+mTSufu/rttttqbpOeLkREhBAajXFbjUaI+HghTp50OqxDPx4X0coloUGZcRcoE7HKBXHklxM1tpkwQbZRqSrbG1620aOF0OuN259IyxRxqnMm+4jCZbHnmyNO74cnbN1q/vABQjz3nImN/vtf498lw02lEqJfPyHKyoyaXz51WXQIOiLUKDdqrkaZ0KBM/PT6nzW6WLxYvv6muujVS4jiYje9IBaUlAhxyy3GvyOGl0FRhPjii5rbrFkjREBAzf1Qq4Vo21aIS5c8vhsO2bFDiLCwmm9btVqIxEQhTp82bq/XCzFihPH7qOr7a5qJP9W+zJ7Pb4cSnKVLl4omTZqImTNnik2bNolNmzaJmTNniuTkZLF06VKxZ8+eipslq1evFj/++KM4cuSIOHLkiJgyZYoICAgQ+/fvN9n+xIkTIjQ0VIwdO1YcPHhQzJs3TwQEBIjly5fbFb/LEpzcveY/4A236h/Cm4fJZMZke5UQq5rUTIp80c/drO978bnK9qW5VxMixfy+733NuI+9r9dMbqomkUuChSitBX+V5s2z/MmlUgnx11/O9fHCC5b7AISo+j7R64Vo167mX/uqfy0HDHAuJiFE9/B9NT5MKz9Uy0WvyN1G7TdutL4ba9YY93F7zHahrpbcVO2jU/BhoddVy4p8UESE9X0vKKiywenT5o+f4dNs9myjPsZem2r2eCjQinpKrijOrcxYLlwQIjDQ8q/um2966AWq4p13aiY3VW8BAUKcq/Lnp6REiJgY4w/46r/uY8Z4fj/spdcLcc015g+7RiPEoEHG23z7rfXfq+3bvbI7DnF7gqMoisWbSqWq+Gqv6Oho8dlnn5l8bOLEiaJNmzZG9z3zzDN2nTUSwoUJzp9jLCQrkGca1vSsbF983nJ7wy3rZ+ficrfSPOv78CWE2DK8cpvDsy0kN1dvyxsIodPK9nqdEMsbWulDEeLQLO+8Bvbo1MnyX2O1WoiJE53rw5ZPx/btK9v//rv19ooixKlTDoe0d/kRq10AQhz8vjK5e+ihmv+ZVv8DPnBgZR/HN5yyqY8//mf6nyZfsWeP9X0AhBgypMpGr71mPcFp1aqiedHFIhGGAqt9LHxmU8U2771n+VcXECIhoeZZNXfS64Vo3NhyTCqVEG+/XbnNV19Zf21DQ4W4csVz++GI1FTb3rZZWZXb3Hab5V8TjUae4akt7Pn8dmiaeEZGhuuukV2l0+nwzTffoLCwED169DDZZuvWrejXr5/RfXfccQc+//xzlJeXIyDAdEHH0tJSlJZWDqTLz7djMKwlubsBYWHApNABl/dV/lxwxHJ7QK46nLcfaHSnS0J0i9zdtrW7vLfK9/vkvlna/9LzQNklILgBUHoJKLUySNbwWvm6gwctj2bU6YB9+8w/bgtbRuKeqjIexZb+hJCxN2niUEj71p8HYH1l3P3rc9D2brl0f3q65XHOWq0cfGtw4LezAKzHt3/jJdwwwmozr/nuO9vaVd137N1r+fdKCDnGSqcD1Gqc3n4WhWhm8fkDUIb96ZUHYO9eWVbCUjfZ2UB+PhAVZds+OKuwUI6ht0RRgP1V/jTs2wcEBMhBueYUFQGZmUBbJ6vYuNO+fXLfhDDfRgjg0CE5lguQle0tDULWauX7zh85lOA0bdrUZQHs27cPPXr0QElJCcLDw7Fq1Sq0a9fOZNuzZ88iLi7O6L64uDhotVpcuHABCQmmi8TNmDEDr7/+ustirqAJg/maK1dVrbyttqHSt3Cy0rcnBNazrZ06uPJ7m6qcV9mm6rYW2/v4awXIQb9Xrph/XKWSox6dYe1TCACqrjBua38hjr++oZG2/XkJCa+crBAWZsPzVgndkT58Uf36trUzqkAeGiqPu6VPr4AA2QZAaIz195QeKqNDbsuviaJ4tjJ6YKD1X3dFQY39sGXGlBO/7h4RGmo5uTGw9xja8r6rjRxel3vRokW46aab0KhRI5y6+p/hrFmz8J2t/4pc1bp1a+zevRvbtm3Ds88+iyeeeAIHD5ofza9UKwEgrh7t6vdXNXnyZOTl5VXcTp8+bVeMZiVZqbulaIAmVUoG1+sEhCRaf97Eu52Ly92iO9mWrLV8uvL7xoMsn71R1EBcH1nuAQACwuXPioUPpuoVy33VvffKQkTm6PWy/LQzbClsVPXs5513Wo4JkNMzund3OKS+Y9oiBJanw4ehELeOrvyHZujQis9jk9Rq4P4qhelvGtkOUbhssY9AlKLfC21sCdlrnnrKtnZGE+IGD7ac3Gg0ss3Vv42JXePRKfiIlWriGqNq4oMHWz6jplbLX6tgG/8fcYXAQPnra6qQqYFWa/yWGjTI8kulKPIt5ML/3d2if3/L+w0ADRrIkiAG991neRtF8Xxle09xKMH55JNPMH78eAwYMACXL1+G7upvTr169TBr1iy7niswMBAtW7ZEt27dMGPGDHTu3BkffPCBybbx8fE4e/as0X05OTnQaDSIjY0120dQUBAiIyONbi6R/CgQklBZZNLI1QKUrcdWuUt9tQq3OSog+REgzMffZQDQ5kXLjwfUA5o/Uflzw1uBmOvNvFaQZ67aTzG+r/1UmJ1qr2jk88X1sTFgL5owQf4VMZWEazRAs2bGn9qOsDiPGDJrqNqmQQPgmWcs1yV75RWn/jWPbByJF3v8CcVM1XkFeoy/eTvC48Mr7hs5Ul7qMPUHWaWS/2k++2zlfcH1gjGp326zMSjQ49lrtyK2lW+vsB4YCPTubblNcDAwfHiVO+65B7jmGtOJquE6xssvV96lUjD9uUuWq4nH7sC1D7auuO+22+R0dXO5sF4PTJ5sOW53mDzZ/JkMjQbo1EkmQQYdOwIDBpj/oBcCmD7dtWX63CEhQSbDlv4JmDJFnrgzeO45+TY2tY1aLc8eGv1e+RNHBvm0bdu2Ypp4eHi4OH78uBBCiH379onY2FhHnrJC3759xRNPPGHysYkTJ4q2bdsa3Tdq1CjvDTIWQoi8w0J827RyqvdXGjn4dVmkENm/1myv18vZQl8qchDyV+rKgcepg+TU6Npi80OmB/9+HSFEQc3pv6I4R4ifr6/2WqmEWBIgxPH5pvs4/oV8/EtVlW0gxM9djWdp+boff5RzOxVFjuozjKRt1UqIEyZeK0e8+675Qcw//lizfWmpEA8/XDnSUK2uHI340ksuGTmqLdWKkW3SZBcoq5iODAgxql2q0JZqa2yza5cQcXGVYRleqvr1hdi2rWYfunKdeLFLqtxVlBv18XjzTaKssKzmRj6qVSvThzAgwMyvSWamnONc9cVSFCFCQoRYudJkH7OHyplUKmiFpsprdWu9XeLyqcs12mdnC9Gli3EXKpWcXbV4sYtfADssWSJEUFDNt1SnTkKcOVOzfV6eEH37Gu+HWi335f33PR+/o0pKhHjwwZpvW0UR4pVXTL9tN20SIjq65nsqMVGIffs8vw/OsOfz26FimyEhITh8+DCaNm2KiIgI7NmzB82bN8exY8fQqVMnFNu4pOSUKVPQv39/JCUloaCgAEuXLsXbb7+NX375BbfffjsmT56MrKwsLFy4EIAc3NyhQwc888wzGDlyJLZu3YpRo0ZhyZIlGDp0qM3xu7zYpr4c+Ps7IHutvGxSvzvQ9GF5mcWcKydlFe7CDCAwFkgeBsReb769r8rdC6RPBK78JcfZNH8KuOZ58/9iCD1wbr1cHFFXBES1B5o9IQcWm1NyHsj4Asg7IC+NJd0LxPUFFNdVPvaI/Hxg0SJgxw75L/uAAcBdd1m/VGSPs2flGaM//5T/nvXvD/z735YvxO/ZAyxeLFc+TkqS/861bOm6mAAcXP0XFr71N7LPa9CooRaPT02qGFhsSmkp8M03QGqq/Ijv3Rt44AHLYySOrTuJBa+dxN9nNYiL1eKxSY3Qcaj1Qc6+ZvVqYOxY4OJF+Z/3008D//qXhcsMOp0sV/7DD3L1uy5d5Cpw9eqZ7SN79znMn3gIR09oEBmmw/0j66HX6E5QVKZPYej1wK+/At9+K1cM7tBB/ppYOHHuERcvAgsXysHQISHyUtTtt1v48yOALVuAr7+Wb8dWreR+JNowcsDXpKfLBTHPn5fzAIYPB1qYf0uhsBBYulTuv0oFpKQAQ4YYD82rDdxeTbxdu3aYMWMGBg0aZJTgzJ49G1988QV27txp0/M89dRT+O2335CdnY2oqCh06tQJkyZNwu233w4AGD58OE6ePInU1NSKbdLS0jBu3DgcOHAAjRo1wqRJkzBq1Ci74mc1cSIiotrH7dXEX375ZYwZMwYlJSUQQuDPP//EkiVLMGPGDHz22Wc2P8/nn1tY0h/AggULatx3yy23YNeuXfaGTERERHWIQwnOiBEjoNVqMXHiRBQVFeHhhx9G48aN8cEHH+Chhx5ydYxEREREdnFoEENxcTEeeeQRnDp1Cjk5Odi2bRvGjRuHxo0buzo+/+Yv1cTJtwghBxjYUzLYA9XEHWJ3OXGt3A9LK7pVV1goXy8br9brdPZ3gaIiuyp9O8TOauJkOyHkS1tkedUD8jEOJTiDBg2qGPir0Whwzz334P3338fgwYPxySefuDRAv+Qv1cTJt+j1wCefAK1by7nW4eFyHZuVK81v44Fq4nYTwoFy4llyPmxUlNyPiAjgySeB48fN97N8OXDDDXLfo6KANm2A//7X7IpwOTnAuHGycnVMjJyy/uijctVYs77/HujVSzauV08O3p4927WJpJ3VxMl25eXA++/LlRzq1ZOH8dZb5ZhuqgUcmaYVGxtbURBz3rx5olOnTkKn04lly5bVqBXli1w6Tdxe/lJNnHyLTifEsGGVxWiqFuUBhPj3v2tu44Fq4g6xu5z4CTmv3NR+REXJQk/Vvf66+T4ee0y+nlVkZQmRlFSzpo9GI2sY/fGHif147z3TfSiKEPfeK4S25hR5u9lZTZxsV1Ym681WL9BpeKnnzPF2hHWT24tthoSEiFNXi/Ddf//94rXXXhNCCJGZmSlCQkIceUqP8mqC4y/VxMm3LF1qvQrf3r2V7T1UTdxuDpUTv93yfnTqZJwU7dplvY+qldeFEPfdZ74IqFotRPPm1XKiI0fMl6423ObPd+61cqCaONnu448tH0KVynv/A9Rl9nx+O3SJqmXLlvj2229x+vRprFmzpqIAZk5ODqddW1JyAcj8BuZLFuiBoky5ng6RPT76yPJ67BqNvHxlsHWrLKRpbv16nQ74+WdZfdCTPv7Y8rpAGo3cV4MTJ4B16yzvx969wPbtlfd98onlPtRqoz7OnpVX+cxdVdLpZBgbNlS5c+5cy8dDpQI+/ND847awMgsVgPN91GHWXjpFAebN80ws5BiHEpzp06djwoQJSE5Oxo033lhR/Xvt2rXo0qWLSwP0K/ZUEyeyx9691ksGVy1FbU81cU+yu5y4jePWqpaW3r3bch86nVz88KrDh60XalSpjLvAnj2W+9DrbY/dHHuqiZNdhJDH3dJ4bZ2u2jEnn+PQNPH77rsPvXr1QnZ2Njp37lxxf0pKCu6910oByrrMX6qJk+8JCZEzgcxRFOOSwR6oJu4Qu8uJO7AfYWGVtZpsaG9LF0JUe6ls6cPZEtx2VhMn+wQFyYWhzVGrfb/6eF3n8G9+fHw8unTpAlWVN88NN9yANm18u2qvV/lLNXHyPffdZ73kw5Ahld97oJq4Q+wuJ36TnAFlSWCgcSX1qq+DKRqNrAtx1XXXAXFxljdRFFl1o8LgwZaTG43G+RLOdlYTJ9spinzpLL1FdDrjiuXke5jae5I/VRMn3zJ2rPxrbK5kcHy8nNNs4IFq4g6xu5x4MDBpkvnnUxTZvmrRpMcfBxo2NN9HQICccn6VRgNMm2a+C5XKRD2jBx+UBYJM9aEocqPx480/qS3srCZO9jG8dKbeIhqNnPHPCxa+jQmOp7V8Guj4GgBFjrdR1IBy9Q9U4kDghrnejI5qq1atgB9/lGu6KIr8C2z44EtMBNavl2vDVPX++8CwYfJ7jUZ+GBs+kF96SRbt9LQGDYDffgPq16+My7AfMTFyzZdGjYy3mTQJePFF+b1hHwzbPPYYMHOmcfuoKPl6JCTU7CMiQg6urla1cMwYYOpU+dJW7+Lee4E5c6rtR0iI7CM52bgPRZGXlr77Dmjf3t5Xx1hAgKyA2apVzT6Cg2W10utrYQFfH3HddcCKFfJQVn9LtWghX/raVqiyrnGo2GZt5xPFNv2lmjj5litXgK++ArZtk3+Nb79dnkcPCDC/jQeqidvNoXLix+RigH//La8pPfaYXCDQnLIyYNUq+Uml0wE9e8qEz8I4oIwM2cXJkzIHe/hhoGtXC/uh1crF/n75Ra4ad/318kxa9WTTGQ5UEyfbXb4MLFoE7NolT2jedRcwYIDlSXLkPm6vJl7b+USCQ0RERHax5/Obl6iIiIjI7zDBISIiIr/DBIfIzzhSTBxFRXKAiaWFP6r1UVAgh/z4lJISuR92lH22s5i4Q9XEL1yQi0JbWzCwRie+VjCzuFgOSql7IxvsVlYmD6EvrbNoeN/a9behFmOCQ+QnHCkmjnXr5IDisDCgeXM5iLdTJzmi0gQhgP/9D+jQQRb5joiQs02WLPHyZ96OHTLukBC5H2Fhcgr1+vVmN7GzmLhD1cT/+U/5OjVoIIt8BwXJpX7M5pGXLslZYfXrV3by4INy1WJvSkuTayeFhckXIDEReOstmxPiumTfvsqx6jExcoWCl18GLl70Xkx6vSwr0a6d/H0MDwe6dQOWLfNeTB7h5rpYPsmrxTaJ3MCRYuIWC3SqVEJs3mzUXK8X4umnzffxyiue2dca0tKMK3ZXv33zTY1N7Cwm7lA18XvuMR9SQoIQpaXVNjh/XoiWLU13EhQkRGqq614zeyxeLF8cUxXLe/USorjYO3H5oE2bhAgOrlmY1VCQ9dw5z8ek1wsxYoT59+20aZ6PyRluryZe2zHBIX9jbzFxodMJERhoeYMGDYz6+PFH631s2eLZ/RZCCBETYzmo4GCjjMWBYuJ2VxP/5RfrfTz1VLX9GDHCfCcqlRDx8UKUl7vvdTTlwgXLvycqlRBvvunZmHyUVitEYqL5XFujkcmzp337rfXfxe3bPR+Xo9xeTZyIfIu9xcQxb5718R3nz8uq41fNmWO9jxoL3rlbaqq8rGNJSQnwxRcVP9pZTNyhauJTp1oP/auvqvxw+TLw5ZfmO9HrZSA//GD9iV3piy+sFw396CMvX5/0DT//DGRlmb/EqdUCS5da/3V1NVv+Nnz8sefi8SQmOER+wN5i4ti40bYn/u23im9377azD0+oEp9FVfbXzmLiDlUTz8iwHlJxcZUfjh+3nnBqNJ4vX713r/VindnZlgu91hH79lkv71ZeDvz1l2fiMdizx/r7Nj3dc/F4EhMcIj9grapx9WLiNq+kW6WdLVW1bS3u7TK2LtQZHl7xraHQtyVVX09HqonbvYS/LZ3o9Z4vX21LXIri+ZplPig01LZZcr54CC0s3l2rMcEh8gP2FhPH2LHWn1RRgBEjjPqwdKpbpXK+QLbdRo60rd0LL1R8a2cxcYeqiRtVFjfDUEIKgJz61ry55Q30ellg05MGD7Z8ukutlpXag4M9FpKvGjjQ+pW6pk2dL0FmL2vvW0XxwvvWUzwwJsjncJAx+ZujR+VYWlMDHNVqOWsnP7/aRh07Wh55OGCAUfPMTCHCw833ERMjx6R6XL9+lvfj2muNml++LERcXM1JQYYxsyEhQvz1l3EXH35oeZztk08at7940fTzV719+221/ViwwHIn993n+tfOGp1OiC5dzA9+VhTvze7yQQ89ZHlC32efeT6mjAw508/c+7ZBAyFycz0fl6M4i8oKJjjkj377TYjISPmZo9FUfiY1aSLEoUMmNigoEKJpU9N/ibt2rTlXWshZUrGxlbNCDH3Exwuxe7fbd9G08nL5IWxqP5o3F6KwsMYmBw4I0bhxzf2IijL9ea3XCzF1auVsabW6cpuhQ03PlP7tN/N5wRtvmNmXGTPkJ1H1Tvr3F+LKFadeJodlZ1e+voYXS6WSs6sWL/ZOTD6qsFCIu++ufKkMh1FRLBxzD9i0SYjo6Jq/74mJQuzb5724HGHP5zeLbbLYJvkRR4qJY9kyYOZMOb0jIQF4/XUgJcVs86Ii4Ouvgc2b5entvn3l4nVeH4axbh3w2mtytlFsrFwwb+hQs80dKCZudzXxkhLg1VeB1atlH127Au+9BzRqZGE/Tp8G5s+XU7Pq1ZNB3XCD9YFD7qTXyxfq22/l6OgOHWTV+dhY78Xko4QAtm+Xi19euiSvPI4YATRp4t24CgvlLK4tW+Tl5JQUebnW7vFiXsZq4lYwwSEiIqp9WE2ciIiI6jQmOEREROR3mOAQGQgB5OVVW4GtFhLCgXLi9tFr9Tiz6yzO7s1xWx8A5KCiggL39mEnR6qJFxXJX626NyCAyHuY4BCVlQH/939yFGC9enJlrJQU21fJ9RUOlRO3swutHiNabUJIQBkSu8YjoXNDhCpFeL5TGvRaG1Y5s4UQsmzBtdfKhQYjI4GOHeXoXi9mCI5UE//+e6BXL9m2Xj1ZuH32bMtLyxCRi7h5RpdP4jRxqlBaKkRKSs1FIgxzO72xcIUjHConbmcX5TrRJuDY1afWV3m55Pfdw/dafxJbTJhgHHvVfRo9Ws7Z9jBHqom/957p3VAUIe69VxZnJCL7sNgmka3mzAHWr6+5xrpOJz+TRo2StXZ83TffyHmpgPFZDsN+TZsmi+U44Y3b0nC4vOXVn6pOWZbfb7vSEXMfTXOqD2zaBLz7rvy+6jEx7NPHH8vp4B42dqz8Nahe00erBUpL5UzuquEePQpMmCC/r74bQsjp6YsWuT9uorqMCQ7VbR9+aPlxvR74/HPPxOIMu8uJ22/Oxg4ALF0iEpixtJlTfeDjjy3XnNBojEt9e4Aj1cTnzrVe1sLarx4ROYcJDtVdpaVy5TZr4zo8XcHZEXaXE7ffRRED4zM31SnI1jV0qg+kp1seoOKFkuWOVBPfs8fybuj1wIEDromPiExjgkN1l0Zj+d9sQH5yebr8ryPsLiduPzUsJFAAAAENnBw9a0uMHi5Z7kg1cVsqlnt95WciP8cEh+outRq4+27Ll0S0WlnrwNfZXU7cfl3Dj8DyJSqgd30LU4psMXSoTCrNUauB++93rg87OVJNfPBgyycGNRo/ruBM5COY4FDd9sor8pPI1L/bGg3Qrp1Mgnzd2LEyXlPJgVoNxMfLOc1OmDXL8BqZ+uQWUCDwwVf1neoDI0fKKe6mzqypVPLUyLPPOteHnTQaOUbbHJVKlmVKTKy878EH5aoDpnZDUeQ248e7PFQiqoIJDtVt3bvLCnRBQfKTR6OpPBPSpg2wdq31y1i+oFUr4Mcf5do31fcjMVHOFIuIcKqLG5/qgP8+shFKRYIjYEh2VNBj6dhtaHW7k4OMGzSQ6w/Vv5ooVd2PmBh5PCxWqnSPMWOAqVPlS6tWy5shrHvvlZPxqgoJkS95crL82bAbiiIveX33HdC+vUd3gajOYbFNFtskQC5N+8UXcnRocDAwcCBw552WL5f4IofKidvnwpELmDjwEDafbAwFAre1Po0ZP3dBZGMXvpdKS+XU99RUeYatd2/ggQe8Ph7K3mriWq1c7O+XX+TKx9dfL0+kOZlrEtVZrCZuBRMcIiKi2ofVxImIiKhOY4JDREREfocJDpGvKyuTY4QsLeTnrPJy2YetVSCFkFW+r1xxX0wO0Grtr/RNVFcY3raFhd6OxDOY4BD5qoMHgUcekTOjYmLk7aWXgPPnXdfHsWPAiBGVfdSrBzz/PHDmjOn2QgD/+x/QoYOs8h0RIReKWbLEq5W+s7KA556TM8xjYmRYTz4JHD/utZCIfIZeD8ybJ1e9iIyUb/du3YBly7wdmXtxkDEHGZMv2roVSEmRpyKqnlVRq4HGjeUsqfh45/rYvVvOTiouNu5Do5FThLZtA5o2rbzfUHx07lw539nwp0Olkn9BX3kFmDHDuZgckJEB9OgBXLxYczfCwoCNG4FOnTweFpFPEAJ46ilg/nzTb9tp04A33vBujPbgLCormOCQT9PrgebNgb//Nn1ZyrAMrqF6uCOEkGdhjhwx3YdaDdxxh1xbx+Cnn4yX6zVlyxaZbXhQv35yzRlzu9G+vczlrJVOIPJH331nfTH27dvlGZ3agLOoiGqzdeuAU6fMj7nRaoHly527VLV1q7wEZq4PnQ74+WcgM7PyvjlzrFcsr77inZudOCFfLku7sXev/ANOVBd99JH1t+3HH3suHk9igkPka/bts756slYLHD3qXB/WCCGTIIPdu91esdxetlbkrg0F4YncYc8e62/b9HTPxeNJTHCIfE1oqLxMZY0zq/raWpG7ah+2bOODlb4Bry+ATOQ1trxHwsLcH4c3MMEh8jW2FPdMTAQ6d3a8jzvvtF59PCZG1uoyuO8+y2eWVCqPl8i+6SY5c8qSwEA5ToeoLrL2tlUU/61szwSHyNc0aSILFlmqgzVtmnNFQBs0AJ55xvLI21dekUVIDUaPlqdCzFUsr1dPTtfwoOBgYNIk848riiw+HhvruZiIfMlzz8m3sbm3bf36wPDhHg/LI5jgEPmi//4XuOce+b1GU1nCWlGA6dNlcuKs998Hhg2r2Qcg19uZMMG4fVKSrOYdHV25jeEsUIMGciqTFzKJSZOAF1+U31ev9P3YY8DMmR4PichnJCcDa9ZUnums+raNj5dv23r1vBWde3GaOKeJky/bsUNWB794EWjWTP6rlZzs2j727AEWLwZycmQSM3w40LKl+fZFRcDXXwObN8uEq29fYOhQ47M9XnDsmKz0/fffQFycTG46dvRqSEQ+o7AQWLpUruSgUslltoYMkZdwaxOug2MFExwiIqLah+vgEBERUZ3GBIeIiIj8DhMc8oqSEuDyZduWe3FYaakHOvEAe6uJCwHk59tVMtjeYuJ1ma5Mh9yMyygvqt0ly/V6+fYoLfV2JETuwQSHPGrzZmDAALn4VHQ00KgR8O9/y3GrLvPnn8CgQZWdxMXJmUcFBS7sxAPsrSau1wOffAK0bi2nTISHy3VsVq4024W9xcTrspwD5zHuujREBxUhpnk9hIUJPNpsMw79ULtKlufny1UGGjaUb4/QUFmraMcOb0dG5GLCi9566y3RrVs3ER4eLho0aCAGDRokDh8+bHGbDRs2CAA1bocOHbK537y8PAFA5OXlObsLZIelS4VQqYRQq4WQpxnkTaUSont3IQoLXdDJ998LodGY7qRTJyFqyzHfskWIkBC5L1X3Q60WomlTIbKzjdvrdEIMGybbKIrxfgNC/PvfNbpITxciIqJmFxqNEPHxQpw86ZE9rRWydmaLJPXfQo1y49cKZSIUV8Qf/9vv7RBtcvmyEB061Hx7qNXyuP/0k7cjJLLMns9vr57BSUtLw5gxY7Bt2zasW7cOWq0W/fr1Q6ENp9aPHDmC7OzsilurVq08EDE56vJlOftYiJpXWvR6edLF6fVKiorkGQ+dznQnBw4A//qXk514gF4v16cpK6t5zUinA7KygHHjjO//5pvK6uJVJ0YaLs9Nm2ZUf0oI+VIVFdXsQquVJ4lGj3bR/viBsfdkIFvXEDoYr/6sRQBKEYRhT4dDr/X9S6GvvgocOlTz7WF4yzz8sLx8TOQXPJBw2SwnJ0cAEGlpaWbbGM7g5ObmOtwPz+B43uzZxicWTN0aNBBCq3Wik/nzLXcACBEZKURJiat2yz1++cX6fmg0QuTkVG7Tq1fNf8urt3/22Yrmv/9uvQtFEeLUKS/sv4/J3nNOqKC1+nr9+n87vR2qRUVFQoSFWT/uCxd6O1Ii82rNGZzq8vLyAAAxMTFW23bp0gUJCQlISUnBhg0bLLYtLS1Ffn6+0Y08y5YC2efPA5cuOdHJ3r1AQIDlNvn58gyIL3OkmvjevXZV+nakmHhddXj9Gehh+XiooMP+rb49xuv0aevjzgMCWHmd/IfPJDhCCIwfPx69evVChw4dzLZLSEjA3LlzsWLFCqxcuRKtW7dGSkoKNm7caHabGTNmICoqquKWlJTkjl0gC2yt+hwc7GQntqxb6eGK13ZzpJq4tXLZimJUMphVuG0XGmUlaQYgoCAk1EJdLx9gyzHX63nMyX/4zErGY8aMwY8//ojNmzejcePGdm07cOBAKIqC1atXm3y8tLQUpVXmQubn5yMpKYkrGXvQhg1yRX9z1Gqgd29ZF8VhO3YA119v/nGVCujSxfeni2RmynIMlt6aiYnAqVOVZ3qee07WrzI3z1tRgDlzZOVJyLNljRpZnhYeEyNnU3m5AoPXaUu0aBx2Cef0DQCYTmJU0CFzew4SuyV4Njg7CAFce608Q2Mpf05Pl+2IfFGtW8n4+eefx+rVq7Fhwwa7kxsA6N69O44dO2b28aCgIERGRhrdyLNuvVXmHhqN6cf1emDKFCc76dZNZlHmLu/o9cA//+lkJx7gSDXxsWPli2uuZHB8vHzOqxwpJl5XaYI1mDbkECwlN8NbbfHp5AaorNNqLrlRq4E77mByQ37E3QOCLNHr9WLMmDGiUaNG4ujRow4/z9ChQ0WfPn1sbs9Bxt6RkyPE9ddXjnnVaOQs5oAAOT7YJS5dEuLmm2t2otEIMWeOizrxgKIiIQYPrtwPtVreFEWI6dOF0OtrbvPbb3IQtaJU7jsgRJMmQphYRqG0VIiHH67ZBSDESy+Z7qKu0uv0YupNG4QCnVCjXKhRJjQoE4AQQxO3iOLcYm+HaLPZs+VxNrwtDL8mt94qp5ET+TJ7Pr+9eolq9OjR+Oqrr/Ddd9+hdevWFfdHRUUh5OqF4MmTJyMrKwsLFy4EAMyaNQvJyclo3749ysrKsHjxYrz99ttYsWIFhgwZYlO/LLbpPXq9vAy1apWcoty+PfDEE/KMgssIAWzaJKdOFxTIhe9GjJBnMWobe6uJX7ki22/bJs/o3H67XMXNwuBre4uJ12UZG09jwT+P4+TfGtSP1uLhF+PQ9dG23g7LbtnZwPz5cpx6ZCRw//1Ar16Wz+gR+YJaU01cMfNumj9/PoYPHw4AGD58OE6ePInU1FQAwP/93/9h7ty5yMrKQkhICNq3b4/JkydjwIABNvfLBIeIiKj2qTUJjrcwwSEiIqp9at0gYyIiIiJXYoJD3uGRcuJ+4soVICNDlm4gIiKbMMEhz/JIOXE/sXq1HFAcEQE0by7nbHfrJosJERGRRRyDwzE4nvP117Kan6IYlxVQqYAbbgB++833Vxn2lLlz5UI1pqjVwK5dQKdOno2JiMjLOAaHfI9Hyon7Ca3WcilvnQ4YONBz8RAR1UJMcMgzFi0CSkvNlx/Q62UpAUsFI+uKd96x/jpkZvJSFRGRBUxwyDM8Uk7cT2zdalu7q2tDERFRTUxwyDM8Uk7cT4SH29aO48eIiMxigkOeMWiQ5dLVajXQp4+cMVTXjR1rvY1KJdfXJyIik5jgkGd4pJy4n+jRQ04Lt+TBB4HAQM/EQ0RUCzHBIc9QFODHH4EuXeTPGo28qVSyEOT//gfcdpt3Y/Ql27cDCQmmH+vdW1bHJCIis8z8O03kBg0ayCrXbi8n7gdiYoAzZ4CFC4EPPpDT7JOS5KKIvXp5OzoiIp/Hhf44UJOIiKhW4EJ/REREVKcxwSEiIiK/wwSnttEWAWV55lcE9ld6vRyHUlLi7Ug8r6wMyM3lKs++QqeTx6O83NuREJEFTHBqi6yfgHW9gWVhwPJ6wOrmwOFZgN7C2jL+oLgYePNNIDFRVh8PDQX69wc2bfJ2ZO538CDwyCNy4b+YGHl76SW54jN5Xk4OMG6c/D2MiQHCwoBHH2XJDCIfxUHGtWGQ8ZEPgZ0vAIoaEIb/4hX5JfFu4OaVgMoPJ8QVFwMpKcAff8gzOAZqtTyD9eWXwEMPeS8+d9q6Ve57ebnxAolqNdC4sZyNFh/vvfjqmjNngO7d5deqZ9I0Grke0YYNwA03eC8+ojqCg4z9yZUTwM4X5fei6iUKIW9ZPwAn5nshMA+YObNmcgPIDxghgBEj5KUCf6PXA8OGyUtT1Vd/1umArCx5JoE8Z+xYIDu75mVCrVYWkR02rObvKRF5FRMcX/fXPLlIniVHZnsmFk8yVBc396EhhPxgWbjQs3F5wrp1wKlT5sfcaLXA8uW8VOUpZ88CK1eaLzWi0wEnTsizOETkM5jg+LrLe6udualOAPl+OAbg0iU55sEStRrYv98z8XiSLZXXtVrg6FHPxFPXHT5s/eyMSuWfv4tEtRgTHF+nCYPVw6QK8kgoHmVrVfGQEPfG4Q2hobZd7vDHffdFoaHW2wjB40HkY5jg+LrEQQAsfNgpGiBpiMfC8ZjwcFld3NKZDK0WGDzYYyF5zN13W2+TmAh07uz+WAi47jogLs5yG0UB7rrLM/EQkU2Y4Pi6JvcB4c1lIlPD1bE5bSd4NCSPmTrV/JkMjUZWJ+/Tx7MxeUKTJnL6scrC23PaNOuXscg1NBr5epujUgHDh8ukk4h8BhMcX6cOAvr+CoQ3kz8rmqvJjgKoQ4Heq4BoP/1PPiUFWLBAVhtXqSorkAPy7MUPP1gfgF1b/fe/wD33yO81GpnMqNVyf6dPB555xrvx1TVjxsiEW1Eqj4Xhd/Hee+WAeCLyKVwHpzasgwPIBf3O/ASc+RHQlwEx3YDkR4HAKG9H5n7nzwNffAEcOCDHQ9x7L9C3r+UzHP5ixw7gq6+AixeBZs3kmYLkZG9HVXdlZMik++RJoH594OGHga5dvR0VUZ1hz+c3E5zakuAQERHVcVzoj4iIiOo0JjhERETkd5jgEJHPunL2CgrOFLi3k8JCID9frmVDRH6DCQ4R+RShF/hy9O+4NuQwIhLCEZkYgY7BR7HgH5sh9C5MQpYvlwUyw8OBqCigTRs5e401pYj8AhMcIvIpE29Mw6Of3IR9Ja0q7jtQ2gIjPu+F5zpvdE2S869/AfffD+zcWXnfsWPAqFFyphqTHKJajwkOEfmMTR/twbs7bgUA6FG5kKG4+v3H+2/Bund2OddJejrw6qvy+6qJjOES1aJFwKpVzvVBRF7HBIeIfMbHM69Ag3Kzj2tQjo8+MFPV21affFK5SJ8pajXw0UfO9UFEXscEh4h8Rnp2ArQIMPu4FgHYfaGxc53s3i3rmJmj0wF79jjXBxF5HRMcIvIZYZpSAJbH2ISqSp3sJMx6iQ9WBieq9ZjgEJHPGNrrHFQwP8BXDS3uvzHTuU6GDLH8uEYDPPCAc30QkdcxwSEinzHygw6IUvKhRs1LSCpoEYYiPPtBG+c6efxxoGFD09XYVSpZ3PW555zrg4i8jgkOEfmMBm3r47cvz6G+6hIAOajYMOg4RsnD2s8y0ei6eOc6iYoC1q8HEhLkz1Wr1EdEAD//DLRo4VwfROR1LLbJYptEPqc0vxTfTNqB1PV6CAH0vkXBA+90RUiMC8fGlJXJ6eC//ioHFvfsCQwbJsfoEJFPYjVxK5jgEBER1T6sJk5ERER1GhMcIiIi8jtMcIiIiMjvMMEhIiIiv8MEh4iIiPwOExwiIiLyO0xwiIiIyO8wwSEiIiK/wwSHiIiI/A4THCIiIvI7THCIiIjI7zDBISIiIr/DBIeIiIj8DhMcIiIi8jtMcIiIiMjvMMEhIiIiv8MEh2oHvR64fBkoKfF2JEREVAswwSHfVlwMvPkmkJgIREcDoaFA//7Apk3ejoyIiHyYVxOcGTNm4Prrr0dERAQaNmyIwYMH48iRI1a3S0tLQ9euXREcHIzmzZvj008/9UC05HHFxUBKCjB9OnD2rLxPCGDdOuDWW4GlS70aHhER+S6vJjhpaWkYM2YMtm3bhnXr1kGr1aJfv34oLCw0u01GRgYGDBiAm2++Genp6ZgyZQpeeOEFrFixwoORk0fMnAn88Ye8PFWVTicTnREjgNxc78RGREQ+TRFCCG8HYXD+/Hk0bNgQaWlp6N27t8k2kyZNwurVq3Ho0KGK+0aNGoU9e/Zg69atNvWTn5+PqKgo5OXlITIy0iWxk4vp9UBCApCTY76NogD/+Q8wdqzn4iIiIq+x5/Pbp8bg5OXlAQBiYmLMttm6dSv69etndN8dd9yBHTt2oLy83OQ2paWlyM/PN7qRj7t0yXJyAwBqNbB/v2fiISKiWsVnEhwhBMaPH49evXqhQ4cOZtudPXsWcXFxRvfFxcVBq9XiwoULJreZMWMGoqKiKm5JSUkujZ3cIDjYtnYhIe6Ng4iIaiWfSXCee+457N27F0uWLLHaVlEUo58NV9mq328wefJk5OXlVdxOnz7tfMDkXuHhQJ8+8iyNOVotMHiwx0IiIqLaQ+PtAADg+eefx+rVq7Fx40Y0btzYYtv4+HicNcyouSonJwcajQaxsbEmtwkKCkJQUJDL4iUPmToVSE01/ZhGA3TpIpMgIiKiarx6BkcIgeeeew4rV67E+vXr0axZM6vb9OjRA+vWrTO6b+3atejWrRsCAgLcFSp5Q0oKsGABEBAAqFQyqdFczck7dwZ++EEONCYiIqrGq7OoRo8eja+++grfffcdWrduXXF/VFQUQq6OrZg8eTKysrKwcOFCAHKaeIcOHfDMM89g5MiR2Lp1K0aNGoUlS5Zg6NChNvXLWVS1zPnzwBdfAAcOyIX+7r0X6NtXJj1ERFRn2PP57dUEx9yYmfnz52P48OEAgOHDh+PkyZNIrXKpIi0tDePGjcOBAwfQqFEjTJo0CaNGjbK5XyY4REREtU+tSXC8hQkOERFR7VNr18EhIiIicgUmOEREROR3mOAQERGR32GCQ0RERH6HCQ4RERH5HSY4RERE5HeY4BAREZHfYYJDREREfocJDhEREfkdn6gm7mmGxZvz8/O9HAkRERHZyvC5bUsRhjqZ4BQUFAAAkpKSvBwJERER2augoABRUVEW29TJWlR6vR5nzpxBRESE2YKfvio/Px9JSUk4ffp0naujVVf3va7uN8B9r4v7Xlf3G+C+27LvQggUFBSgUaNGUKksj7Kpk2dwVCoVGjdu7O0wnBIZGVnn3gAGdXXf6+p+A9z3urjvdXW/Ae67tX23dubGgIOMiYiIyO8wwSEiIiK/wwSnlgkKCsKrr76KoKAgb4ficXV13+vqfgPc97q473V1vwHuu6v3vU4OMiYiIiL/xjM4RERE5HeY4BAREZHfYYJDREREfocJDhEREfkdJjg+bMaMGVAUBS+++KLZNqmpqVAUpcbt8OHDngvUBV577bUa+xAfH29xm7S0NHTt2hXBwcFo3rw5Pv30Uw9F6zr27re/HG+DrKwsPProo4iNjUVoaCiuvfZa7Ny50+I2/nDcAfv33R+OfXJyssl9GDNmjNlt/OV427vv/nC8AUCr1WLatGlo1qwZQkJC0Lx5c/zrX/+CXq+3uJ0rjnudXMm4Nti+fTvmzp2LTp062dT+yJEjRqs/NmjQwF2huU379u3x66+/VvysVqvNts3IyMCAAQMwcuRILF68GL///jtGjx6NBg0aYOjQoZ4I12Xs2W8Dfzjeubm5uOmmm9CnTx/8/PPPaNiwIY4fP4569eqZ3cZfjrsj+25Qm4/99u3bodPpKn7ev38/br/9dtx///0m2/vL8Qbs33eD2ny8AeCdd97Bp59+ii+++ALt27fHjh07MGLECERFRWHs2LEmt3HZcRfkcwoKCkSrVq3EunXrxC233CLGjh1rtu2GDRsEAJGbm+ux+Nzh1VdfFZ07d7a5/cSJE0WbNm2M7nvmmWdE9+7dXRyZe9m73/5yvIUQYtKkSaJXr152beMvx92RffenY28wduxY0aJFC6HX600+7i/H2xRr++4vx/uuu+4STz75pNF9Q4YMEY8++qjZbVx13HmJygeNGTMGd911F2677Tabt+nSpQsSEhKQkpKCDRs2uDE69zl27BgaNWqEZs2a4aGHHsKJEyfMtt26dSv69etndN8dd9yBHTt2oLy83N2hupQ9+23gD8d79erV6NatG+6//340bNgQXbp0wbx58yxu4y/H3ZF9N/CHYw8AZWVlWLx4MZ588kmzRY/95XhXZ8u+G9T2492rVy/89ttvOHr0KABgz5492Lx5MwYMGGB2G1cddyY4Pmbp0qXYtWsXZsyYYVP7hIQEzJ07FytWrMDKlSvRunVrpKSkYOPGjW6O1LVuvPFGLFy4EGvWrMG8efNw9uxZ9OzZExcvXjTZ/uzZs4iLizO6Ly4uDlqtFhcuXPBEyC5h7377y/EGgBMnTuCTTz5Bq1atsGbNGowaNQovvPACFi5caHYbfznujuy7Px17APj2229x+fJlDB8+3Gwbfzne1dmy7/5yvCdNmoRhw4ahTZs2CAgIQJcuXfDiiy9i2LBhZrdx2XG363wPuVVmZqZo2LCh2L17d8V91i5RmXL33XeLgQMHujg6z7py5YqIi4sT7733nsnHW7VqJd566y2j+zZv3iwAiOzsbE+E6BbW9tuU2nq8AwICRI8ePYzue/755y2ehvaX4+7IvptSW4+9EEL069dP3H333Rbb+Mvxrs6WfTelNh7vJUuWiMaNG4slS5aIvXv3ioULF4qYmBixYMECs9u46rjzDI4P2blzJ3JyctC1a1doNBpoNBqkpaVh9uzZ0Gg0RgPULOnevTuOHTvm5mjdKywsDB07djS7H/Hx8Th79qzRfTk5OdBoNIiNjfVEiG5hbb9Nqa3HOyEhAe3atTO6r23btsjMzDS7jb8cd0f23ZTaeuxPnTqFX3/9Ff/4xz8stvOX412VrftuSm083i+//DJeeeUVPPTQQ+jYsSMee+wxjBs3zuJVClcddyY4PiQlJQX79u3D7t27K27dunXDI488gt27d9s0uwYA0tPTkZCQ4OZo3au0tBSHDh0yux89evTAunXrjO5bu3YtunXrhoCAAE+E6BbW9tuU2nq8b7rpJhw5csTovqNHj6Jp06Zmt/GX4+7IvptSW4/9/Pnz0bBhQ9x1110W2/nL8a7K1n03pTYe76KiIqhUxqmGWq22OE3cZcfd4fNO5BHVL1G98sor4rHHHqv4+T//+Y9YtWqVOHr0qNi/f7945ZVXBACxYsUKL0TruJdeekmkpqaKEydOiG3btom7775bREREiJMnTwohau73iRMnRGhoqBg3bpw4ePCg+Pzzz0VAQIBYvny5t3bBIfbut78cbyGE+PPPP4VGoxFvvvmmOHbsmPjyyy9FaGioWLx4cUUbfz3ujuy7vxx7nU4nmjRpIiZNmlTjMX893gb27Lu/HO8nnnhCJCYmih9++EFkZGSIlStXivr164uJEydWtHHXcWeC4+OqJzhPPPGEuOWWWyp+fuedd0SLFi1EcHCwiI6OFr169RI//vij5wN10oMPPigSEhJEQECAaNSokRgyZIg4cOBAxePV91sIIVJTU0WXLl1EYGCgSE5OFp988omHo3aevfvtL8fb4PvvvxcdOnQQQUFBok2bNmLu3LlGj/vrcRfC/n33l2O/Zs0aAUAcOXKkxmP+fLyFsG/f/eV45+fni7Fjx4omTZqI4OBg0bx5czF16lRRWlpa0cZdx10RQgjbz/cQERER+T6OwSEiIiK/wwSHiIiI/A4THCIiIvI7THCIiIjI7zDBISIiIr/DBIeIiIj8DhMcIiIi8jtMcIiIiMjvMMEholph+PDhGDx4sE1tb731Vrz44otujcdWqampUBQFly9f9nYoRHUKExwiIhfxpcSKqK5jgkNERER+hwkOEdlk+fLl6NixI0JCQhAbG4vbbrsNhYWFAID58+ejbdu2CA4ORps2bfDxxx9XbHfy5EkoioKlS5eiZ8+eCA4ORvv27ZGamlrRRqfT4amnnkKzZs0QEhKC1q1b44MPPnBZ7GVlZZg4cSISExMRFhaGG2+80aj/BQsWoF69elizZg3atm2L8PBw3HnnncjOzq5oo9Vq8cILL6BevXqIjY3FpEmT8MQTT1RcNhs+fDjS0tLwwQcfQFEUKIqCkydPVmy/c+dOdOvWDaGhoejZsyeOHDnisv0jopqY4BCRVdnZ2Rg2bBiefPJJHDp0CKmpqRgyZAiEEJg3bx6mTp2KN998E4cOHcJbb72Ff/7zn/jiiy+MnuPll1/GSy+9hPT0dPTs2RP33HMPLl68CADQ6/Vo3Lgxli1bhoMHD2L69OmYMmUKli1b5pL4R4wYgd9//x1Lly7F3r17cf/99+POO+/EsWPHKtoUFRXh3XffxaJFi7Bx40ZkZmZiwoQJFY+/8847+PLLLzF//nz8/vvvyM/Px7ffflvx+AcffIAePXpg5MiRyM7ORnZ2NpKSkioenzp1Kt577z3s2LEDGo0GTz75pEv2jYjMcKoOOhHVCTt37hQAxMmTJ2s8lpSUJL766iuj+9544w3Ro0cPIYQQGRkZAoB4++23Kx4vLy8XjRs3Fu+8847ZPkePHi2GDh1a8fMTTzwhBg0aZFO8t9xyixg7dqwQQoi//vpLKIoisrKyjNqkpKSIyZMnCyGEmD9/vgAg/vrrr4rH58yZI+Li4ip+jouLEzNnzqz4WavViiZNmhjFVLVfgw0bNggA4tdff62478cffxQARHFxsU37Q0T203g1uyKiWqFz585ISUlBx44dcccdd6Bfv3647777oNVqcfr0aTz11FMYOXJkRXutVouoqCij5+jRo0fF9xqNBt26dcOhQ4cq7vv000/x2Wef4dSpUyguLkZZWRmuvfZap2PftWsXhBC45pprjO4vLS1FbGxsxc+hoaFo0aJFxc8JCQnIyckBAOTl5eHcuXO44YYbKh5Xq9Xo2rUr9Hq9TXF06tTJ6LkBICcnB02aNLF/p4jIKiY4RGSVWq3GunXrsGXLFqxduxYffvghpk6diu+//x4AMG/ePNx44401trFGURQAwLJlyzBu3Di899576NGjByIiIjBz5kz88ccfTseu1+uhVquxc+fOGjGFh4dXfB8QEFAjNiGEyXgNqj9uSdXnNzyPrckREdmPCQ4R2URRFNx000246aabMH36dDRt2hS///47EhMTceLECTzyyCMWt9+2bRt69+4NQJ7h2blzJ5577jkAwKZNm9CzZ0+MHj26ov3x48ddEneXLl2g0+mQk5ODm2++2aHniIqKQlxcHP7888+K59DpdEhPTzc6yxQYGAidTueKsInISUxwiMiqP/74A7/99hv69euHhg0b4o8//sD58+fRtm1bvPbaa3jhhRcQGRmJ/v37o7S0FDt27EBubi7Gjx9f8Rxz5sxBq1at0LZtW/znP/9Bbm5uxUDbli1bYuHChVizZg2aNWuGRYsWYfv27WjWrJnTsV9zzTV45JFH8Pjjj+O9995Dly5dcOHCBaxfvx4dO3bEgAEDbHqe559/HjNmzEDLli3Rpk0bfPjhh8jNzTU6q5OcnIw//vgDJ0+eRHh4OGJiYpyOn4gcwwSHiKyKjIzExo0bMWvWLOTn56Np06Z477330L9/fwBy/MrMmTMxceJEhIWFoWPHjjUWvHv77bfxzjvvID09HS1atMB3332H+vXrAwBGjRqF3bt348EHH4SiKBg2bBhGjx6Nn3/+2SXxz58/H//+97/x0ksvISsrC7GxsejRo4fNyQ0ATJo0CWfPnsXjjz8OtVqNp59+GnfccYfRZa8JEybgiSeeQLt27VBcXIyMjAyXxE9E9lOEPReRiYjsdPLkSTRr1qzG5ZzaTq/Xo23btnjggQfwxhtveDscIqqGZ3CIiGxw6tQprF27FrfccgtKS0vx0UcfISMjAw8//LC3QyMiE7jQHxHVKpmZmQgPDzd7y8zMdEu/KpUKCxYswPXXX4+bbroJ+/btw6+//oq2bdu6pT8icg4vURFRraLVao1KIFSXnJwMjYYnp4nqOiY4RERE5Hd4iYqIiIj8DhMcIiIi8jtMcIiIiMjvMMEhIiIiv8MEh4iIiPwOExwiIiLyO0xwiIiIyO/8P62PmO80K3vGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGxCAYAAACeKZf2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABcuklEQVR4nO3dd3hUZfo+8PvMTEhISDFAIJEEgiBSpKPASpMFBBYBC1gBXXX5/qQIoixYsRDLqtgAYZEioq4EFBUUVBIsoNKRLiBETKSnkjZzfn88TJJJzpyZTMmU3J/rmgvm1PdMuJgn57zveyuqqqogIiIiChIGXzeAiIiIyJNY3BAREVFQYXFDREREQYXFDREREQUVFjdEREQUVFjcEBERUVBhcUNERERBhcUNERERBRWTrxvgCxaLBX/++SciIyOhKIqvm0NEREROUFUVubm5SEhIgMFg//5MrSxu/vzzTyQmJvq6GUREROSCjIwMNGnSxO76WlncREZGApAPJyoqysetISIiImfk5OQgMTGx7HvcnlpZ3FgfRUVFRbG4ISIiCjCOupSwQzEREREFFRY3REREFFRY3BAREVFQYXFDREREQYXFDREREQUVFjdEREQUVFjcEBERUVBhcUNERERBhcUNEREFtaIi4MIFwGLx/LEvXpRjq2rVdaoKZGcDBQWePy/p82lxk5KSgm7duiEyMhJxcXEYMWIEDh48qLtPWloaFEWp8jpw4EANtZqIiALBTz8Bw4cD4eHAZZcBjRoBTz4J5Oa6f+z0dGDQICAiQo7dpAmQkgIUFgIlJcBrrwHJyUBMjGzTuzewbp375yXnKKqqVW/WjBtuuAG33XYbunXrhtLSUjz22GPYs2cP9u3bh4iICM190tLS0K9fPxw8eNAmOqFhw4YwGo1OnTcnJwfR0dHIzs5m/AIRURD6/HNg5Ei5e2I2ly83GoF27YBNmwBX//tfvhwYMwYwGGyPbTAAPXsCkZHAl1/a3s0xGmXbN94AJk507bzk/Pe3T4ubyk6fPo24uDikp6ejd+/emttYi5vz588jJibGpfOwuCEiCl4FBUB8vNyh0fqGMxqBKVOAl1+u/rHPnAEuvxwoLtZeryja57QyGIDffpO7OlR9zn5/+1Wfm+zsbABAbGysw207deqE+Ph49O/fHxs3bvR204iIKEB89BGQk2O/yDCbgQULpC9OdS1ZApSW2l/v6HaBogALF1b/vFQ9flPcqKqKqVOn4rrrrkO7du3sbhcfH48FCxYgNTUVq1atQqtWrdC/f39s2rTJ7j5FRUXIycmxeRERUXDaswcICdHfJicH+PNP147tIJBal9ksxyDvMvm6AVYTJkzA7t278f333+tu16pVK7Rq1arsfY8ePZCRkYH//Oc/dh9lpaSkYNasWR5tLxER+afwcMd3UACgbl3Xju1OcWMwuHZeqh6/uHMzceJErFmzBhs3bkSTJk2qvX/37t1x+PBhu+tnzJiB7OzssldGRoY7zSUiIj82fLj+oyODAejWDWjc2PPHBvSLH4tFOjqTd/m0uFFVFRMmTMCqVavw7bffItnFHlY7duxAfHy83fWhoaGIioqyeRERUXDq1g24/nrpOKzFYgEef9y1Yw8cCHTsCJh0nnsYDNoFjskENG8O3HSTa+cm5/n0sdSDDz6IFStW4NNPP0VkZCSysrIAANHR0ah76b7djBkzcPLkSSxbtgwAMGfOHDRr1gxt27ZFcXExli9fjtTUVKSmpvrsOoiIyL+sXAnceCPw/fflhYjFIoXH66/LOlcYDDJfzeDBwM6dtsc2mYDFi2WI+ejRMsGftcAqLZXC5quvgNBQty+PHPBpcTNv3jwAQN++fW2WL168GOPGjQMAZGZm4sSJE2XriouLMW3aNJw8eRJ169ZF27Zt8cUXX2DIkCE11WwiIvJzl10mc9ls2iSFTm4ucNVVwLhxrj2OqqhxY2DbNuDrr4FPPpEi5uqrgbFjgfr1ZZs//wTee0+2q1MHGDpUXk5Ox0Zu8qt5bmoK57khIiIKPAE5zw0RERGRu1jcEBERUVBhcUNERE4pKQHOn7fNU/IHxcXAsWNAXp6vW0L+gsUNERHpOnRIOsvWqwfExkrS9eTJQGamb9t15Ahw7bUy+qh5cwmsTEyUDsRUu7FDMTsUExHZtWMH0Ls3UFhoO3md0Qg0agRs3gwkJdV8uw4elBFKJSXa6+fMkQKMggs7FBMRkVtUFbjjDhnqXHlWXrMZOHUKmDDBN20bNsx+YQMAU6dKQUa1E4sbIiLS9MMPwIED9vvYlJYCn38O1HSizfHjgE7iDgCZVI+RgrUXixsiItLkTAK2qgL799dMe6y++8657X75xbvtIP/F4oaIiDR5M13bHTExzm3H9O3ai8UNERFpuuEG/YBIQOIGrr22Ztpj5Uy7AGDSJO+3hfwTixsiItLUqBFw//0SFmnPzJmSnVSTTCbJiNJz+eXAgAE10hzyQyxuiIjIrjlzgFGj5O8mkwwBt4Y/PvIIMGWKb9q1cCEwaJD2ugYNgO3ba7Y95F84zw3nuSEicmjnTmD5cuD0aZnXZtw44IorfN0q4KefgBkzgBMnZBK/CROAe+7Rv9tEgcvZ728WNyxuiIiIAgIn8SMiIqJaicUNERERBRUWN0RE5BRHqeD5+UBOjvbcOBYLcOECUFSkvW9xsfcSx4uK5NwWS/X3NZulXXpRD/aoKpCdDRQUVH/fQJaX5/uEdhY3RESky1EqeGoqcM01sj46GrjqKmD+fCkmcnOBp56SYeWXXSYTA954o3QEBoDdu4HRo4GICDl2gwbAo48C58653+6ffgKGD5dzXnaZtOHJJ6VNjmRlyTXGxEi7IiLkMzh0yPG+JSXAa68Bycmyf0SEhI+uW+fuFfkvVQWWLQPat5eO3ZGRQIcOwHvvOTcRpBcaVPtkZ2erANTs7GxfN4WIyK9t366q9eqpqsmkqvI1JS+jUVUTElR16lR5bzCUr1MU+XP0aFXt0EG2rbyvyaSqzz+vqqGh2sdu0UJVT51yvd2ffSbH1Tp3hw6qqvff/4kTcm2V9zWZ5LPYutX+viUlqvqPf5R/BhXPC6jqG2+4fk3+ymJR1YkTq/47sP590iTZxhOc/f7maCmOliIi0qSqQJs2ElKp9bjIaHT8GMlgsP84SFHkpbXeaATuvhtYvLj67S4oAOLj5Q6N1jec0Sjz87z8svb+w4cDX3xh/5pbtJA8La3crQULgPHj7d+tMBiA336TuzrB4ttvgf79HW/Tr5/75+JoKSIicoujVHBn+sfo9XNRVfvrzWZgxQrpK1NdH31kv++P9dgLFmj3/zl5EvjsM/1rPngQ+P577fVvvqnfNkWRCQiDydtv68dhmEzA3Lk11x6AxQ0REdnhTCq4NxUXA0eOVH+/PXuAkBD9bXJygD//rLp83z7HfUQURc6hxdH+ZrP9fQPVjh1Aaan99aWlNT9jNIsbIiLS5GwquDe5kuztTpp5eLjj/VTVfrvCwvT3NRiCL608IsLxNvXqeb8dFbG4ISIiTc6mb3tL8+Yy8qq6hg/Xv5NgMADdugGNG1dd162bjNjSYzIBQ4Zorxs5Uv8zs1hkm2Byyy36cRcGg2xTk1jcEBGRJmdSwSMjy4M0KzIY9NPCjUagYUPtfa2eeMK1jKhu3YDrr7d/bIsFePxx7XV16kjSuT0GA3DfffLZaJk2Tf7UepxnMknBdtNN9o8fiP71L/v/DoxGICoKeOCBmm0TixsiIrLLUSr4li1AQkL5eutdi8hI4KuvpCOpySRFQcX13btLPwxrsnfFYxsMwOzZEs7pqpUrgR49bNtlbcPbb8tcO/Y89JDMtQOUt8na7lGj5DOxp2NHYPVqefSkKLbX3Lw58M03QGio69fljxo3Br7+WuYDAmyvOTYW2LDBfjHoLRwKzqHgREQO6aWCl5TIF/rXX8vjoB49gDvuKO+LkZUFLFkiI68iI+URRe/e8uWvqsDPPwMffCAjo5o3l1TvxET326yqwKZNUujk5sojrnHjtB9HaTlyRNp94oTcZbrrLilenJGdLRPYbdsmd4OGDpWX3p2qQFdYCPzvf/KZA0CfPsCttzruh1QdTAXXweKGiIgo8HCeGyIiIqqVWNwQERFRUGFxQ0REAc2a3F1cXP19VVUm9MvP93y7ApWjBPdAwOKGiIgC0rlzwPTpMi+NNbl79GhJGnfEYgHeeUc6GUdHyyRz11wjCee1VW6upKZXTHAfPrw8wT2QsEMxOxQTEQWcM2dkVNaxY7Y5UNYh5V99JaN1tFgsMmrqvffKR2wB5SGfs2bJl3xtkpMjI9h+/dX28zQa5TNatQoYNsx37bNih2IiIgpajz4K/P571YDL0lIZmn7bbfZnKV61SgobwDamwRri+dRTkpdUmzz7bNXCBpD3ZrMMgy8o8E3bXMHihoiIAsqFC8D779svXiwWmVvn88+117/1lv58MyYTMG+e280MGEVFkpJuLwnd2i/po49qtl3uYHFDREQB5cgRx52HTSa5E6Fl9277X+SAFE27drnevkBz8qQUL3pCQgIrzZzFDRERBRRnkrstFvvp245SuRXFuXMEC3eT0P0RixsiIgoorVpJTIMei8V+ftQttzhOO7/5ZtfaFogaNwa6dNEPKS0tDaw0cxY3REQUUAwG/dFMBoMUMC1baq+fNEkes2h9mRuNQFwcMGaMZ9oaKJ58srxDdWVGo6Ssd+1as21yB4sbIiIKOGPHAikpUqBUTu4eNEgCL+254gpg3ToJ8QRsU6wTEoCNG4HaNkvIjTdKWrpWgnuPHhI+Gkg4z01t+xdMRBREMjKAxYuBo0eBmBjg9ttlMj5Fcbxvfj6wYgWwebMURwMGyKOXkBCvN9tvVU5wv/VWoFcv5z7PmsBUcB0sboiIiAIPJ/EjIiKiWonFDREREQUVFjdERORzxcWS7G1vcr28PHlpcZQKXlQksxrbGw2kx3rskpLq7+vP8vNl4r5g7ZjC4oaIiHxm925J8o6IkGTvBg0kN+rcOfnife89oEMH6dwaGQm0bw8sWybrHKWC//STpFqHh0vKdaNGMuQ5N9dxu7KygMmTpZOy9dhjxwKHDnn14/C61FTpcF2vnqShX3WVpKO7Uvj5M3YoZodiIiKfSE+XYdtms21OlNEIJCcD/foBCxeWp3UD5X+//34Zsm0vFfypp6SQUdWqKdft2gGbNtkf7p2RAXTvDvz1V9Vjh4UBaWky6V2geeYZ+Vwqfp7WVPS775ZRUnoT+fkDjpbSweKGiMi3SkuBpCQpILTuGlT8ArbH3jYGg3xh2/t2MxqBKVOAl1/WXj98OPDFF9qPyIxGoEULYP9+/xke7YwdO4DOnfW3+fhjmfzQn3G0FBER+a0vvgAyM+0XMM48JtHbV+/XdrNZUrCLiqquO3kS+Owz+31/zGbg4EHg++8dt8+fzJunHzlhNEpaerBgcUNERDVuzx7H+U7elJMD/Pln1eX79jnuZKsogZWQDQA7d9o++qvMbC7vqxQMWNwQEVGNCw/3fSdWrZTrYEzIBqRDtKPHaIF2TXpY3BARUY0bNsx3xY3BAHTrJmnYlXXrJqOv9JhMwJAh3mmbt9x0k/56k0miFoIFixsiIqpxLVvKl6nRaH+bunW11xuN+ncZDAagYUP7x7ZYgMcf115Xpw4wc6b+se+7T4aVB5IxYyTtXOszMRgkT2vixJpvl7ewuCEiIp9YvFiGggPlQ7iNRvmynT1bhmvHxpavt/bRiY2VdXqp4Nu3S5p1xX2taddvvy0p2PY89JDMtQNUPfaoUcCcOZ78FGpGdDTw7beSeg7Yfp6RkZKSfsUVvmufp3EoOIeCExH5jKoCP/8MfPCBzCLcvDlwzz1AYqKsLyyUIcrp6fK+d28pMMLC5L1eKriqShG0cqVM3HfVVcC4cdqPo7QcOSJzv5w4IXeC7roL6NjRo5df44qLgdWrga+/lk7EPXvKZxYR4euWOYfz3OhgcUNERBR4OM8NERER1UosboiIiCiosLghIqpFCgqA7GztiepUVfq9XLyovW9JiX5yd7AnTfsTRz8rd1gscmytGZwDhU+Lm5SUFHTr1g2RkZGIi4vDiBEjcPDgQYf7paeno0uXLggLC0Pz5s0xf/78GmgtEVHgWrcO6NNHOo7GxEgw5WuvScFSWCgjj5o0kfTs8HBgwAAJiAQkCXvsWEmSjo2V/SdPlvgEQDtpev5830/SF4wq/6wiImR0mLXDtTtycyVstFGj8n8Hw4dLunrAUX1o0KBB6uLFi9Vff/1V3blzpzp06FA1KSlJzcvLs7vP0aNH1fDwcHXy5Mnqvn371IULF6ohISHqypUrnT5vdna2CkDNzs72xGUQEfm1t96SGEmj0RonKS9FUdUbblDVXr1U1WCwXWc0yvpnnlHVevVU1WSquj4hQVWnTpX3FfdXFPnzrrtU1Wz29dUHj4sX9X9W773n+rGzs1W1Q4eq/0aMRvnZr1njsctwi7Pf3341Wur06dOIi4tDeno6evfurbnN9OnTsWbNGuzfv79s2fjx47Fr1y5s3rzZqfNwtBQR1Ra//y7zl+jdRbEOm7bHXvq20Wj/EZXV//4XXDPf+lJKikw+aO9nWaeOBH86mmFZyyOPyJ08rZ+noshcOJmZzsVTeFNAjpbKzs4GAMRaZ23SsHnzZgwcONBm2aBBg7B161aUlJR4tX1ERIFm4ULHmUKOfsW192XqqLAJtqRpX1JV+Sz1itTSUpmXp7qKiiQl3d7PU1WlL9VHH1X/2L7iw0xWW6qqYurUqbjuuuvQrl07u9tlZWWhUaV5rxs1aoTS0lKcOXMG8fHxVfYpKipCUYWeUTk5OZ5rOBGRH9uzx3ER4i1mc+ClZ/sreynmFSkK8Ouv1T/2yZNyfD0hIYH1s/SbOzcTJkzA7t278cEHHzjcVqn0a4j1yVrl5VYpKSmIjo4ueyVap74kIgpy4eHyWMlXrDMJk3tCQx3fgVMU15K9gzEJ3S+Km4kTJ2LNmjXYuHEjmjRportt48aNkZWVZbPs1KlTMJlMqF+/vuY+M2bMQHZ2dtkrIyPDY20nIvJnw4c77m/jLcGWNO1LYWHAwIH6QaOlpcCIEdU/duPGQJcu+kVwaSkwcmT1j+0rPi1uVFXFhAkTsGrVKnz77bdITk52uE+PHj2wYcMGm2Xr169H165dERISorlPaGgooqKibF5ERLXBzTdLh2KTRicERdH/QjOZgPr19QugyMjakzTtazNn2i9UTSbJvRowwLVjP/mk/WMbjcD11wNdu7p2bF/waXHz4IMPYvny5VixYgUiIyORlZWFrKwsXKwwK9GMGTMwZsyYsvfjx4/H8ePHMXXqVOzfvx/vvvsuFi1ahGnTpvniEoiI/FqdOhKS2Ly5vLemQVsfYaxeDaxYIdtZU7OthVDbtpKuPXp0+b7WlGxARths2WI/aXrtWqBFi5q71mDXuzewfLn2z6pdO5nLyNVHkDfeKGnpFdPTrcfu0UPCRwOJT4eC2+sjs3jxYowbNw4AMG7cOPz+++9Is84mBZnEb8qUKdi7dy8SEhIwffp0jB8/3unzcig4EdU2ZjPwxRdScBQVyWOIu++WSfcA4OxZYOlS6TRat648zhowoPzLcudO+WI9fRpISpJ07SuukHUlJeVJ06Wl8mV4xx2BkzQdaCr/rEaMAP7+d8/0rcrKkhFXBw5IgXrrrUCvXt59fFkdTAXXweKGiIgo8ATkPDdERERE7mJxQ0REREGFxQ0REQFwnDRdUAAcOwYUF1f/2MGQNF1dqioJ7AUFvm5J7cPihoiolnOUCv7VV9J5OCJCRl2Fhsqw4507HR87Nxd46inbpOkbbwzQpGknlZRITlNysiSoR0TISKd163zdstqDHYrZoZiIarHCQmDQIOD7723nOTEa5f348cC8edr7GgyyX48e2utzc+VLvXIEhNEoo29SU6XQCSbWye6++MI2s8saMvrGG5z7xx3sUExERA7NmVO1sAHki1hV7Rc2gOyjN2vts89qZ1uZzfK66y4gP9/lpvuld9+tWtgA5Z/BQw/Joz3yLhY3RES1lDNJ04789Zf2I6biYsdJ07m5wIcfun5uf/Tmm/rrFUWS2sm7WNwQEdVSubmSCO2ub7+tuuzPP6UzrZ6QENdSrP3Zvn1V79pUxKT0msHihoiolnImadoZkZFVlzmTIB1oSdPOcJSCbjAE3zX7IxY3RES1VGiodCbWS5p2RFEkiqGyRo2Abt0cJ027kmLtz0aO1A4ptXLUT4k8g8UNEVEt5ihpOjRUf/9hw4B69bTXPfGEftJ0v37ANdc439ZAYM1w1rojZjLJUPqbbqrZNtVGLG6IiGqxXr2A99+3nwq+d68EZWrp1k0CM+0ZNgyYO1c7abp7dxkKHmw6dpTPpG5dKXAqXnPz5sA33zguGMl9nOeG89wQETlMBf/oI+Dll4Hz54GEBGDWLOD66507duWk6Vtukflv/CVp2huys4H33gO2bZPCcehQebnzCJCYCq6LxQ0REVHg4SR+REREVCuxuCEiIqKgwuKGiCiIFOaV4tjObBTmlWquP3MGOHHCtVmJS0qkz429WYe9KStLXlrMZmmXK2nlqgrk5LgeA3HxoqSde6ODR0GB9N2pfZ1H3MfihogoCPyQmok2l/2JupFGNO8UjbqRRrS+LBPff5wJAHjmGUmobtgQaNpUOrmOGCFfoI4cOgSMHStDvmNj5TiTJwOZmd68ovLgzrp1gfh4eYWFAQ88IOvOnQOmTwcaNJB2RUQAo0cDu3c7d+x33gGuugqIjpZru+Ya50dwpafLHEEREZJ23qSJJKsXFrp3zYCkh/fuLceOiZF08ddek+KSnKTWQtnZ2SoANTs729dNISJy25f/PaEqMKuARZXf860vi6rArF7XJa/S8vJXw4aqevGi/WNv366q9eqpqslku5/RqKoJCap6/Lh3rslsVtWOHbXbDKhqu3aq2ry5tKPicpNJVUNDVTUtTf/Yd98t2ytK+b4Gg/w5a5Z+2957T/arfG6DQVV79dL/PB15883yz7fisRVFVf/xD1UtKXH92MHA2e9vjpbiaCkiCnBRxnzkWsIBaI2ttv4Xb3/c9R13yFw3VfZUgTZtgMOHtR9FmUzA4MHAmjWutFrfq68CDz+sv42iaD+yMRiAuDggI0N7tuCVK4Fbb9U/9vbtQKdOVZefOQNcfrn9R2AGA/Dcc8CMGfrH13LsGNCihf1HhooCzJ8vd65qK46WIiKqBda8dRy5lgjYL14UnXVi5Urt5T/8IHPT2OtjU1oKfP65FBGe9uqrjrex96u5xSL9cz7/XHv9W2/pzzdjMgHz5mmvW7JErtsei0WO78ptg4ULHc/94yh1nASLGyKiAPbNF+538igu1r5bsGeP4y9bVQX273e7CVX89Zd7+5tM9hPHd+/W7xRdWgrs2qW9zpnP5M8/JXG9uhy1S1Udp46TYHFDRBTAoqI8M82vVsBleLhzX6TeSLl2dyZfi8V+uxy1V1Hk2rWEhzsubhRFOmxXV0SEftAo4Lkk92DH4oaIKID937PxKO9X45qmTbWX33CDfsI1ANSvD1x7rVun19Szp3v7WyzAjTdqr7vlFsfXdfPN2suHD9d/LGU0yiiqsDDn2ln52HpD9E0m++0iWyxuiIgCWMKVkbgu8TjsFziqzjrx4ovayxs1Au6/X/9uwsyZrt2lcGTOHNf3NRikgGnZUnv9pElASIj2dRmN0hl5zBjtfQcOlHBMe8WRxeJaZ2JACpfmzbWPbb1b46iTNQkWN0REAe6bg4loHWOddEZFxYLmqugsfLPegpAQ7X0ff1zmhrFnzhxg1Cj5u8kkX/7WR0aPPAJMmeKBC9DQvr0EeWo9glEUWZeSIgWKtU3WomDQIOn4a88VV8hcMpGR8r5icndCArBxI2BvII7BIPu2a2e7r8EgRd7y5TJHjStCQyU1vHlz22MrijxKW71aCityjEPBORSciILE5/NO4LlngdO5oWgYWYSZM1XcOEGeORUXS5L36tUyGVzHjsArrwBJSc4de+dO+eI+fVr2GTdOigRvu3AB+Pe/peAAgD59gJdeksntABmptXgxcPSoLLv9dpmMz5l+Kfn5wIoVwObNUhwNGACMHAm7hWBFFgvw9dfAJ5/ILMVXXy0THdav79p1VlRaCqxdC3zxhfzcunQB7r5bJhus7ZgKroPFDRERUeDhPDdERERUK7G4ISIioqDC4oaIqIZZU6z9LghRVaWTy8WLvm4JkVtY3BAR1ZCsLEnTjokpT7EeMwY4eNDHDSsslKFHTZpIxHV4uPSuTUvzccOIXMMOxexQTEQ14I8/gO7dpcCpOMW+ySRDgNPTZVRMjSsslLHT339vO4Oc0Sjvly6VoTpEfoAdiomI/MiECVULG0CG/RYWSjK3T37VnDOnamEDSENVFfjnP2X8N1EAYXFDRORlJ08Ca9bYD0U0m4FDh4BNm2q2XVBVibDWm/PfbNafEY/ID7G4ISLysv37Hd+VURT7KdZek5srlZcenzSMyD0sboiIvMxewnRFquqddG1dzkRMW+f+JwogLG6IiLysWzegQQP9bUwmYMiQmmlPmdBQ6UxsDYvSUloKjBhRY00i8gQWN0REXhYSAjz2mP31BoP0223cuObaVGbmTPt9bkwmoEMHicImCiAsboiIasDkycD06fKUp3KK9a23Aq+/7qOG9eoFvP++RFobDLYR2W3bAl9+KcuJAgjnueE8N0RUg44ckcFHJ04ADRsCd94JdOrk61YBOHtW5rTZs0f62AwfLhP5sbAhP8JUcB0sboiIiAIPJ/EjIiKiWonFDREREQUVFjdERDXNnVhwVQWys4GCAs+3y5GSEmm3vamWA5CqAjk5QH6+r1tCnsTihoioprgTC15aKjlQzZvL/hERQO/ewNq1Xm40JBti7FigXj1pd0yMXEdmpvfP7SUWC/DOO8BVVwHR0XJp11wDpKb6umXkCexQzA7FRFQT3IkFLy0FbroJ+Pxz2xwHo1GO9frrwKRJ3mn3jh1SRBUWSjsqnrtRI2DzZiApyTvn9hKLBRg3DnjvPRmab/1IDQZZN2sW8OSTPm0i2cEOxURE/sSdWPDFi4HPPqu63nqsKVOAo0c932ZVlXZdvGhb2FjPfeqUXFeAWbVKChvA9iO1zmX41FNS01HgYnFDRORt7saCv/mmfgaUogALFrjfzsp++AE4cMB+u0tL5W5SRobnz+1Fb72lnzhhMgHz5tVce8jzWNwQEXmbu7Hg+/bp7282y+R7nrZnj+NgTVWV6wsgu3fr94kuLQV27aq59pDnsbghIvI2d2PBw8L09zUYnDtHdYWHOy7KgIBLDXfUXEXxzsdJNYfFDRGRt7kbCz5yZHnekxaLxTvJ3TfcoH9eAKhfH7j2Ws+f24tuucXxZd18c820hbyDxQ0Rkbe5Gws+bZr8qfWIyGQCkpPlG9vTGjUC7r9fP19q5kwJ3QwgkybJj0TrsoxGIC5ORuhT4GJxQ0RUE9yJBe/QAfjkE3meoii2yd3JycA338hwcm+YMwcYNUr+bjKVtx0AHnlERmoFmCuuANatAyIj5X3FjzMhAdi4EeAsIYGN89zwXzAR1SR3YsGzs4Hly4GtW+VuyZAhwNChjp+xeMLOnXLu06dlXptx46RKCGD5+cCKFTJVj9EoIegjR8pdHfJPXk8F/+abb/DNN9/g1KlTsFgnB7jk3XffdeWQNYbFDRERUeBx9vvbpXJ/1qxZeOaZZ9C1a1fEx8dDcTRUkIiIiKiGuFTczJ8/H0uWLMHdd9/t6fYQERERucWlDsXFxcXo2bOnp9tCROQ/rMndxcWeP3Z2NvDFF/aDJy9eBC5c0J5jxt1U8DNnpL9Ppe4E1kNfuCCn1xKEoeAUpFwqbu677z6sWLHC7ZNv2rQJw4YNQ0JCAhRFwSeffKK7fVpaGhRFqfI6cOCA220hIgIAnDsno5oaNJAE7PBwGS3kiSlrFy6Uzr8xMcA//iFDcwwGGUkFSHjmDTdI4vdllwGXXw7Mnl0eWulOKvgzz8h+DRsCTZtKh+QRI4CCAhQWAikpQJMmctrwcOlcm5YmuwZhKDgFOac7FE+dOrXs7xaLBUuXLkX79u3Rvn17hFTqWv7qq686dfJ169bhhx9+QOfOnXHzzTdj9erVGKEzEVVaWhr69euHgwcP2nQkatiwIYx6QSGVsEMxEWk6cwbo2VNCKCsndxuNwJdfAn37unbsV18FHn7Y/vquXYFt26TYqXhugwHo0QOIjpbxy66kgo8aBXz8seaqwgZNMOiq3/H9j0abmzlGY3lC9ksvBVUoOAUwj4+W6tevn9Mn37hxo9PbljVEUZwubs6fP4+YmJhqn8OKxQ0RafrnP4GlS7WfuxgMctcjI8O1scLuDLxQFP0YBIMBOHxY7upU9v33QK9ednd9AdPxGGbDonMj31pDVWYyAYMHSyYoUU3w+GgpVwoWb+nUqRMKCwvRpk0bPP7449UqvIiINFnnkLHXocRiAf76S1KwR46s3rHnzHGvbc6Ebi5YALzwQtV106fbPyyAtzABFugXXs6Egicm6jeRqCa51Ofm3nvvRW5ubpXl+fn5uPfee91ulD3x8fFYsGABUlNTsWrVKrRq1Qr9+/fHpk2bdPcrKipCTk6OzYuIyMaRI447D5tM9pO79Xigj6IuvVTwI0fs7paLSJxEE8BBcaMnAEPBqRZwqbhZunQpLmp0p7948SKWLVvmdqPsadWqFe6//3507twZPXr0wNy5czF06FD85z//0d0vJSUF0dHRZa9E/opBRJW5m9ytx43H6E7RSwXXyX0KRREUVB01VV0BFgpOtUC1ipucnBxkZ2dDVVXk5uba3Ak5f/481q5di7i4OG+1VVP37t1x+PBh3W1mzJiB7OzssldGRkYNtY6IAkarVhInoNc3xmwGbryx+sdeutT1dlnptUsvFXz4cLu7haIYg/AVjCi1u40jARgKTrVAtSbxi4mJKRt+feWVV1ZZrygKZs2a5bHGOWPHjh2Ij4/X3SY0NBSh3gqVI6LgoCjAk0/KmGctRqMUEBr/9zkUHy8pjRqP822YTLZDkioyGKSIqdz/xmSSDi/2UsFTUoD58+0edyZm4yvlBumAo9Gc6GiZ20ZjWhzZP/BCwakWqFZxs3HjRqiqiuuvvx6pqamIjY0tW1enTh00bdoUCQkJTh8vLy8Pv/32W9n7Y8eOYefOnYiNjUVSUhJmzJiBkydPlj3qmjNnDpo1a4a2bduiuLgYy5cvR2pqKlJTU6tzGURE2saMkclbZs4sv1OiKFIYDBgggZeuOntW5qYpKdFev3kz8P/+H7BjR3kQpsUif3/3XXm0NWqUzLBnnfqitFRSwb/6yn4qeL16EnN9/fWa5+71eF+830bBuHFyOIOh/NBt28pIqOnTgQ8/lKZYayuzOWBDwakWcCk48/jx40hKSnI7U8o6tLuysWPHYsmSJRg3bhx+//13pF2aSeqll17CggULcPLkSdStWxdt27bFjBkzMGTIkGqdl0PBiUhXRoYUMkeOSFFx223y7MUTOXozZsicNyUlUqTcfLNUDoAUM19/DXzyiRQx7dpJ+nb9+rLenVTw4mKZtGb1ajl3x47AK6+UTVJz9qw8PduzR/rQDB8u9Zy12AnCUHAKQB6f52b37t1On7x9+/ZOb+sLLG6IiIgCj8fnuenYsSMURYGqqg7v2JgZPEJEREQ+4vRoqWPHjuHo0aM4duwYUlNTkZycjLlz52LHjh3YsWMH5s6diyuuuIL9X4iIiMinnL5z07Rp07K/33rrrXjjjTds+rq0b98eiYmJeOKJJ3QjFIiIPMZiAXJygLAwefkLiwX480/pF+PK9BgXLwJFRTJUqfKdcosF+OMPmdemQYOq+5aWyqisyEjn+uIEAbNZ/hlERHDkFgmXJvHbs2cPkpOTqyxPTk7Gvn373G4UEZGuvDzg6aeBxo0lxjoiQuaf2bLFt+0qLQXuvFMKrcRESZaMiNAPzKxILxW8sFBiH0JDJdW7YUMpfqzTbxw/DowfD0RFSXR3VBTwr38Bv//utcv1tcoB7hERwOjRQDW6iFKQcmm0VOfOndG6dWssWrQIYZd+WyoqKsK9996L/fv3Y/v27R5vqCexQzFRAMvNBfr0kW+wiv37jEa5y7Fype7EdV5jsQAtW0qiuJa+fWVItj3vvw/cfbd2Knj37sBvvwGnTmnvO2CAJIrn5NjOZ2MySZHz448ySWEQOXNGwtKPHdMOcP/qK/lnQsHF46OlKvr5558xbNgwWCwWdOjQAQCwa9cuKIqCzz//HNdcc43rLa8BLG6IAtj06TKEWWvggqLIr+9ZWfJnTXrkEcBBFAyWL5c7O5WdOSN3aRxlW+mxTvJXmdEIXHONFDhB5N57gffe056b0GCQp4EZGbXmyVyt4dXiBgAKCgqwfPlyHDhwAKqqok2bNrjjjjsQUdP/obiAxQ1RgCoulm+t7Gz97RYuBO67r2baZHXZZcCFC/rbtGoFHDhQdfkrrwCPPmp/GmBP2LNH5s0JAhcuyBM/R7Xg6tX2UykoMHl8KHhl4eHheOCBB1zdnYio+jIzHRc2ISGuJXe7y1G7ALmVoGXPHvt3Xjzl11+DpripToA7i5vayeniZs2aNRg8eDBCQkKwZs0a3W1vdCVYjojIEWfip11N7naXolTNfaosJER7uTOJ5O4KouhuZz4uiyWoLpmqyenHUgaDAVlZWYiLi4PBYH+QlaIofj+JHx9LEQWwa6+V+AG9uxxbttR8VHWnTpJRoOfmm6XDc2VffSWjpLylbl3gr79keHgQcNR32+rQIdmOgoez399ODwW3WCyIuzRfg8Visfvy98KGiALcE0/YL2yMRhmV5ItBDa+9pr/eYJBMKS0DBkhxZK/3q6NMK0Wxv42iSLplkBQ2gHyUTz6pv/6WW1jY1GYuzXNTUFDg6XYQETnnH/8A5s+XQsBgkD+tRUH37kBqqmcCLqurb1/gjTe0z200Su/WSyGVVRgMwNq1wNVXy3vrNRkMMivde+8BTz2lvW9ICPDNN4C1D6R1LLT1M7nvPuCZZ9y6NH80diyQkiIfkdFoe8mDBrkX4E6Bz6XRUnXq1EHXrl3Rt29f9OnTB9ddd11AjJKy4mMpoiDw11/yDbZ/P1Cvnvyq3qePbwqbik6dkmHhmzfLN++gQcDzz0sbHXGUCn7iBDBtGrBjh3yTjxwptzCsszPv3y/R3pmZQHw8MGYM0KaNt67UL2RkAIsXyyOqmBjg9tvlxp2v/xmQd3h1KPjmzZuRnp6OtLQ0/PjjjygsLETnzp3Lip3Bgwe71XhvY3FDREQUeLw+z42V2WzGL7/8gvnz5+P9998PiH43LG6IiIgCj9fnuTlw4ADS0tLK7uCUlJRg2LBh6MP5romIiMiHXCpuGjdujJKSElx//fXo27cvZs6ciautHeGIiCpSVZngLjS05iceOXdO+qm0aqV97qws6eeSkFB1XXExcPKkBFRq9ZfJyQHOn5fYhMqjnFRVMrCMRu0YCGuMdb162nPf6KWCE5FDLo2Waty4MfLy8nDixAmcOHECf/zxB/Ly8jzdNiIKZIWFwAsvAE2aSDRBeLgMedYLj/SUV1+V89WvL0Osw8Olg+3PP0sxM368FDvx8VKchIXJaCOLRaa/vfZaKcaaN5ch1ImJwP/+J8devFjSyKOjgWbNZDRT797An3/K/u+8A7RuLevr1ZPerda5bbKygMmTpeerNcZ6zBjg4EFZr5cKTkROc7nPzYULF7Bp0yakp6cjPT0de/fuRfv27dGvXz+88MILnm6nR7HPDZGXFRbKKKHvv7edk8ZolPdLlsiXujdMnixDsu258kqZ3U1Lq1ZS3GilMQLAwIHA+vXa68LCZJj6ypW2sxVbYxWmTgU++kgKnMox1qGhwL//LSOftFLBe/YENmwoHxVFVEvVWIfic+fOIS0tDZ9++ilWrFjBDsVEBLz0EjBjhv3J9kJCyh/5eFJenn9PVmcvP8pRrpTBADz7LDBzpvfaRhQAPD5DcUWrV6/G5MmT0aFDB8TFxeH//u//kJ+fj9deew27d+92udFEFARUFXjzTf0va7NZHu942j33eP6YnmTvM3EUmGmxAG+95Ti7iogAuNih+F//+hd69+6N+++/H3379kW7IEmaJSIPyMsD/vhDfxtF8U5y944dnj+mv8jMlE7I0dG+bgmR33OpuDl16pRT273wwgsYP348YmJiXDkNEQWi0FDHj1kUxTsjp4I5BlpR5LMlIodceizlrNmzZ+PcuXPePAUR+Zs6dWTEj9Fof5vSUmDECM+f+5FHPH9Mf2A0SmdmdigmcopXixs3+yoTUaDS60xsMgHt28toKk8bM8a7d28czTljDbuszGjU7+isKDJs3V4quMUinykROcWrxQ0R1VLXXQesWFH+iKpicnebNsCXX2oXAZ6wd6/cPdJyxRWSsK1VpCiKDFG/4QbtfRs0AH76CbA3QmP0aEnntq6veM3x8RKkOX26nKdyjPWoUcC2bfqp4Jz9nchpbg8F1xMZGYldu3ahefPm3jqFSzgUnKiGnD0LLFsG7N4td1SGD5eJ/LxV2FiZzTJnzMKFMudO/foysd/IkbL+wgWZV2bjRhmB1LevDF+39g/86ScZdn38uNxxefBB4N57y/sSvf22HLugQAqml1+Wu1EAkJ8PfPAB8OOPUsD8/e9yXmvBdeSIFFEnTshQ+DvvlIkGAcep4ES1XI3Nc6OHxQ0RERF5ilfnuSEiIiLyV14tbnr16oW6wTw0k4iIiPyO0/Pc5OTkOH1Q662itWvXVr9FRETOslhkYruwMO1h0sXF0gcmKkp/aLoWa3J3RIT9Dsr2OEoFd4eqSrtCQiQQlIiqcPrOTUxMDC677DLdl3UbIiKvyssDnn5a0rkvu0wKiBtvBLZskfV79gC33SbLY2OlQ+6jj0oHZ0fOnZNRTQ0ayL7h4TKaadcux/s6SgV3R2kpMGeOJJXHxMi19e4N8JdIoiqc7lCcnp7u9EH7+PmQRXYoJgpgubkyLHr3btv0bKNRhlnPmiUhk6WltuneRiPQrJkMybYX2HnmjCRwHz1aNbnbaJQh7H37au+rqjKyadky7VTwp58GnnrKtWsuLQVuugn4/HPbfCmjUdr5+uvApEmuHZsogPjFaCl/xeKGKIBNnw688opt8VGRoshLaxJBoxG4+277oZ3//CewdKn2sQ0GKYoyMuSRUGWpqcAtt+i3fds2oHNn/W20LFwIPPCA/fUGA3D4sNzVIQpiNTJaqqCgAAcOHMDu3bttXkREXlFcLI997BU2gNzZsDc7stkskwueP191XXY2sHy5/WNbLMBff8ndEy1vvaXfr8dkAubNs79ez5tv6s+OrCjAggWuHZsoCLkUnHn69Gncc889WLduneZ6s95/PERErsrMlCLEHcXFMpFe1662y48ckXV6TCZJM7dOBlhR5cdklZWWOtdvR8u+fbaPoyozm6WfEREBcPHOzUMPPYTz589jy5YtqFu3Lr788kssXboULVu2xJo1azzdRiIi4ampJbSO48zII1W13wZHbVMU10c3OQrMNBg4coqoApeKm2+//RavvfYaunXrBoPBgKZNm+Kuu+7CSy+9hJSUFE+3kYhIxMXJ6CNX4xsUBUhOltFMlbVqJVEKeo9/zGYZlaXlllvsB19a3Xyz822taORI/WNbLN5JWScKUC79D5Gfn4+4uDgAQGxsLE6fPg0AuPrqq7F9+3bPtY6IqLInnrDfp8ZolE6/9vq+qKrsr1UcKYrkUdl7/GM0SnFy5ZXa6ydNko7G9lLB4+IktdwV06aVt7Eyk0kKNkedmYlqEZeKm1atWuHgwYMAgI4dO+Kdd97ByZMnMX/+fMTHx3u0gURENv7xD2D+/PLU7Irp2927y4gka7K3dQi30SjbPvcccM899o89ZgzwwguybeXk7gEDJPDSnubNZai4vVTwb7+VuW9c0aGDhGnWrSsFTsVjJydLGnloqGvHJgpCLg0Ff//991FSUoJx48Zhx44dGDRoEM6ePYs6depgyZIlGD16tDfa6jEcCk4UBP76S4qN/ftlsrxbbpH5b6xzzPz8M/DhhzIyqnlzmYMmKcm5Y2dkyLGPHJEJ8267Dbj2Wv1HVlaOUsHdYR3RtXWrHG/IEGDoUMePw4iCRI3Oc2MdEp6UlIQGDRq4ezivY3FDREQUeLw6z80zzzyDgoKCsvfh4eHo3LkzIiIi8Mwzz7hySCIiIiKPcOnOjdFoRGZmZlmnYquzZ88iLi7O7+e54Z0bIiKiwOPVOzeqqkLRePa8a9cuxMbGunJIInLEmgadn+/rltSs/Hy5bq3fwywW4MIFoLCwxptFRP6rWsXNZZddhtjYWCiKgiuvvBKxsbFlr+joaAwYMACjRo3yVluJaieLRabWr5w0/fHHvm6Zd61aJaOf6tWT627VSuILLBbHqeBEVKtV67HU0qVLoaoq7r33XsyZMwfRFYY11qlTB82aNUOPHj280lBP4mMpChiqKkOXly7VTpp+8klJwQ42zz1XPh+NdU4b6/WPGgUcOiRxA1qp4CtXAsOH+6bdRORVXh0tlZ6ejr/97W8wBejwQxY3FDBWrwZuukl/G1eTpv3V7t0yr4ueioVe5eUREUBWlvxJREHFq31u+vTpg+PHj+Pxxx/H7bffjlOnTgEAvvzyS+zdu9e1FhNRVc4kTc+dW3PtqQnWCfr02PudTFXlkdUHH3i+XUQUMFwqbtLT03H11Vfjp59+wqpVq5CXlwcA2L17N5566imPNpCoVtu503HS9M6dNdWamrF9u1yXq0JCJLmbiGotl4qbf//733juueewYcMG1Kkw62a/fv2wefNmjzWOqNZzJuk52B6/1Kvn3EzA9ugldxNRreBScbNnzx6MHDmyyvKGDRvi7NmzbjeKiC659Vb9x1KK4nrStL9y1MfIkdJSJmQT1XIuFTcxMTHIzMyssnzHjh24/PLL3W4UEV0ycaJkCNlLmm7YEBg7tubb5U133w00aqRd1BkM+hlNRiPQt68MlSeiWsul4uaOO+7A9OnTkZWVBUVRYLFY8MMPP2DatGkYM2aMp9tIVHslJ9tPmm7c2L2kaX8VGQls3AhYf1GyXrOiyCOrL7/UTwVPTXXvsRYRBTyXhoJbE8E//PBDqKoKk8mE0tJS3HnnnViyZAmMerfR/QCHglPA8WbStL8qKQE++QTYsEEeNfXoAdx+uxQ4gH4qOBEFpRpJBT969Ci2bt0KRVHQqVMntGjRwtVD1SgWN0RERIHH2e9vl2fhW7RoEV577TUcPnwYANCyZUs89NBDuO+++1w9JBEREZHbXCpunnjiCbz22muYOHFiWdzC5s2bMWXKFPz+++947rnnPNpIIiIiIme59FiqQYMGePPNN3H77bfbLP/ggw8wceJEnDlzxmMN9AY+liKqxGIB/vhD5tVp0KD6+585AxQUAE2aaI/s0lNSIrMKR0Y6npm4uoqLpb9SVJT26Ku8PJkXJzLSs+clIq/wavyC2WxG165dqyzv0qULSqsxs+imTZswbNgwJCQkQFEUfPLJJw73SU9PR5cuXRAWFobmzZtj/vz51Wk6EVVUWCjzyoSGAk2bytDy6Gjnwzifew6IiZH9mjaVDs7Dh0uh48jhwxIKWq8eEBsrx5k0CdCYZqLa9uwBbrtNJjiMjQXq1wcefRQ4e1aKmffek/yqyEgpfNq3B5Ytsx/rQEQBxaXi5q677sK8efOqLF+wYAHuvPNOp4+Tn5+PDh064K233nJq+2PHjmHIkCHo1asXduzYgZkzZ2LSpElITU11+pxEdElhIdCsmYRzVvylJCcHePppx5MD3nabJHdnZ5cvM5uBNWvkuHoFzo4dQJcuwPLlcncFkDssc+fK8hMnXLwoAJs2yTw3qanl15WdDbz6qiz/v/8DxoyxjWjYu1fmC5o8mQUOURBw6bHUxIkTsWzZMiQmJqJ79+4AgC1btiAjIwNjxoxBSEhI2bavvvqqcw1RFKxevRojdGYWnT59OtasWYP9+/eXLRs/fjx27dpVrdgHPpYignyZL1umv8033wDXX191+Q8/ANddp7/v7bcDK1ZUXa6qQNu2wKFD2rlZRiMweDDw2Wf6x9dSWgokJckwcYul6nqDQXt5RV9/DfTvX/1zE5HXeXW01K+//orOnTsDAI4cOQJAohcaNmyIXyv8NqR4eL6JzZs3Y+DAgTbLBg0ahEWLFqGkpMSmqCIiB/73P8fbzJgB/PRT1eXTpzve194d1R9/lLlp7DGbgS++kLs3SUmOz1PR2rX6j7UcFTYmkySxs7ghCmguFTcbN270dDuckpWVhUaNGtksa9SoEUpLS3HmzBnEx8dr7ldUVISioqKy9zk5OV5tJ1FAKCx0vM2xY9rLf/vN8b7FxVJMVO5gvGeP431VVQqg6hY3e/ZIgeJqqngwpqwT1UIu9bnxpcp3g6xP1fTuEqWkpCA6OrrslZiY6NU2EgUNezMgh4Y6t7/WyClnks4B15K9w8Md351x5hhEFNACqrhp3LgxsrKybJadOnUKJpMJ9evXt7vfjBkzkJ2dXfbKyMjwdlOJ/F/z5o63GTZMe/nw4Y73bdpUe/kNNzge8l2/vuREVdewYe4VNwaDJLETUUALqOKmR48e2LBhg82y9evXo2vXrrr9bUJDQxEVFWXzIqr1XnxRf73JBKSkaK+bPRtw1MfN3r5xccADD+jPhzNjhmu5WS1aAKNGac9pY1W3rvZ6o1GGhv/rX9U/LxH5FZ8WN3l5edi5cyd2XnrGfezYMezcuRMnLg0DnTFjhk3K+Pjx43H8+HFMnToV+/fvx7vvvotFixZh2rRpvmg+UWC75RYZ8q0lJERGSsXEaK+vV0+Su+0VODNnymgpe157DRg9Wv5uMklhYS04HnkEmDrVmSvQ9u67cneo8rENBpmX57vv5M6Qdb31LlJsrIR02um7R0SBw63gTHelpaWhX79+VZaPHTsWS5Yswbhx4/D7778jLS2tbF16ejqmTJmCvXv3IiEhAdOnT8f48eOrdV4OBSeq4I8/pJjYsUO+6EeOBJ58EggLc7xvcTHw7LMyMqqkRCbGe/VV5zsC79wpc92cPg0kJgLjxsndF3epKvDzz8CHHwLnz8sjuHHjyttVWAh8/DGQni7b9ukjd3ycuWYi8pkaSQUPVCxuiIiIAo9X4xeIiIiI/BWLGyIiIgoqLG6IiIgoqLC4ISIioqDC4oaIiIiCCosbIiIiCiosboiIiCiosLghIiKioMLihoiIiIIKixsiIiIKKixuiIiIKKiwuCEiIqKgwuKGiIiIggqLGyIiIgoqLG6IiIgoqLC4ISIioqDC4oaIiIiCCosbIiIiCiosboiIiCiosLghIiKioMLihoiIiIIKixsiIiIKKixuiIiIKKiwuCEiIqKgwuKGiIiIggqLGyIiIgoqLG6IiIgoqLC4ISIioqDC4oaIiIiCCosbIiIiCiosboiIiCiosLghIiKioMLihoiIiIIKixsiIiIKKixuiIiIKKiwuCEiIqKgwuKGiIiIggqLGyIiIgoqLG6IiIgoqLC4ISIioqDC4oaIiIiCCosbIiIiCiosboiIiCiosLghIiKioMLiprZTVaAkByjN93VLiIiIPMLk6waQj6gW4Mgi4MCrQM4BWRbbFWjzKJB0q2/bRkRE5AbeuamNVBXYch/w8wNAzsHy5ee2A9+PAnY/6bu2ERERuYnFTW30x6fAscWX3qgVVljkj1+fBc5urelWEREReQSLm9ro0FuAYrS/XjEBh+fWXHuIiIg8iMVNbXRhF6Ca7a9XS4HzO2usOURERJ7E4qY2MoY73sYU4f12EBEReQGLm9oo6Rb9x1JQgMSba6w5REREnsTipja6cgJgCIXmj18xAqENgObjarpVREREHsHipjaqlwz0+xIIiZL3ikleABDWGOj/LVAnxmfNIyIicgcn8aut4noBI/8Ajn8InP4RUAxAo/5A4k2AsY6vW0dEROQyFje1mSkCuOKf8iIiIgoSfCxFREREQYXFDREREQUVFjekz2IGis8DlhJft4SIiMgpflHczJ07F8nJyQgLC0OXLl3w3Xff2d02LS0NiqJUeR04cKAGW1wLXPwL2PYQsDIGWBkL/C8C+PFuIJufMxER+Tefdyj+6KOP8NBDD2Hu3Ln429/+hnfeeQeDBw/Gvn37kJSUZHe/gwcPIioqqux9w4YNa6K5tUPBSWB9d+BiZnlMg6VERlb9sRrovxGo3823bSQiIrLD53duXn31Vfzzn//Efffdh9atW2POnDlITEzEvHnzdPeLi4tD48aNy15Go96Mu1QtWycCF7Oq5k+ppYC5EPjhDkBVtfclIiLyMZ8WN8XFxdi2bRsGDhxos3zgwIH48ccfdfft1KkT4uPj0b9/f2zcuNGbzaxdLmYCf3wqhYwW1Qzk/QacSqvRZhERETnLp8XNmTNnYDab0ahRI5vljRo1QlZWluY+8fHxWLBgAVJTU7Fq1Sq0atUK/fv3x6ZNm+yep6ioCDk5OTYvsiN7PwCLg40MwIVfa6I1RERE1ebzPjcAoCiKzXtVVasss2rVqhVatWpV9r5Hjx7IyMjAf/7zH/Tu3Vtzn5SUFMyaNctzDQ5mJicSw2EBjHW93hQiIiJX+PTOTYMGDWA0GqvcpTl16lSVuzl6unfvjsOHD9tdP2PGDGRnZ5e9MjIyXG5z0IvtAoTG6W+jGIGEITXTHiIiomryaXFTp04ddOnSBRs2bLBZvmHDBvTs2dPp4+zYsQPx8fF214eGhiIqKsrmRXYYQoB2j+ltACTfA4Qn1FiTiIiIqsPnj6WmTp2Ku+++G127dkWPHj2wYMECnDhxAuPHjwcgd11OnjyJZcuWAQDmzJmDZs2aoW3btiguLsby5cuRmpqK1NRUX15GcLlyIlD4F7A3RQI1AQCKdDJOvAno9pZPm0dERKTH58XN6NGjcfbsWTzzzDPIzMxEu3btsHbtWjRt2hQAkJmZiRMnTpRtX1xcjGnTpuHkyZOoW7cu2rZtiy+++AJDhvAxiccoCtDheQnUPLoEyD8OhDYAmt0JxHb2deuIiIh0Kapa+yYsycnJQXR0NLKzs/mIioiIKEA4+/3t80n8iIiIiDyJxQ0REREFFRY3gaIkDyjJdW3fi6eB87sBc3HVdaoKlOQApfna+7qTCq6qQHE2UFpQ/X2JiIhc5PMOxaRDVYHfVwD7XwYu7JJl0e2A1g8DyWOl46+e3U8B+14ELEXly+o2Aa5fD0S1Ao4sAg68CuRcSvqO7Qq0eRRIulVSwfelyDaleTJEPGk00PYxIPoq/fNaSoFDbwMHXwfyj8myhtcBbf4NXD7Utc+CiIjISexQ7M8dirdPAw68ArnBZo1EUACoQMv/B3R9y36B88OdwPEV9o/dZATwxyflxwPKz9NqMpCRapsKDgCKCTCG6qeCW8zAdzcDJ9dcWmA9thGAGeg8B7hqst5VExERaWKH4kB3atOlwgawzXq6VCwcngtkbai8l8g/qV/YAJcKmwrHq3ieg68DBZmupYIfWwKc/PTScStuc+lY26cAub/pt42IiMgNLG781aG5cqfEHsUkj360bL7TAw0way92lAp+8E3o/rNSDMBvC9xtHBERkV0sbvzV+R1yp8QetRQ4v1N7nbUPjdfopIJn74NuqrhqBi7s8UqriIiIABY3/ssU4cQ2dhK8DaGebUsVOqngxjAH+xoAozPJ40RERK5hceOvEm+G/uMdo4xq0nLF/V5pks257aWCJ96k/zgNFiBxhDdaRUREBIDFjf9q8QAQEi2FRBUGubPTYrz2vlc/7qDAuHQMrR+/YgRM9fT300sFb/0wZASWxiguxQRENLNflBEREXkAixt/FdYQ6P+NBFYCUhhYC5bQWKDfevsFBgDcsM1OYQQgpiPw9zQgJKrqscMaAwN+BNrOBKDIMRRj+XpHqeAxVwN9Pr30yEyxPXZEM7kmh4+uiIiIXMd5bvx5nhsAMBcBJz6W0UmqCsT1BpJGASY7fV5s9jUDOx4Gfl8OWIqBsEZAt7lA/ABZX5oPHP8QOP2jjGJq1F+KF2MdWZ931PVU8JIc4Nhy4NwvgKGOPMZKGAoYOG8kERG5xtnvbxY3/l7cEBEREQBO4kdERES1FIsbIiIiCiosbgLFxSyg4E/tdeZi4NxOoPCcnfVFQPEFQNWZXM8ed1LBiYiIfIC9O/2ZxQJsnQAcXQxYCmWZIRRodhdwzQLgwm4gbQhQmFm+j6EucPVTQNvpwNlfgL3PA398BsAC1KkPtBwvyd8hDvoauZMKTkRE5EPsUOzPHYq/7Aqc26a9LqI5kH/U/r7xgy8Fa6qVkr2NQFQbYMB3QJ1o7X0LTgLru7uWCk5EROQl7FAc6A6+ab+wAfQLGwDIXCf5U1WSvc1Azj5gz9P29906UR6DuZIKTkRE5GMsbvzV/pe9d2zVDBz5rxQqlV3MBP741H5op6NUcCIiIh9jceOvLmZ59/ileUDBH1WXZ++Hbqo3AN1UcCIiIh9jceOvamImX61kb3tJ4zZ0UsGJiIh8jMWNv2rwNy8e3ABc1hEIv7zqqtguQGic/u56qeBEREQ+xuLGX3V9A5rJ2k5T7AdnwgK0e1J7lSEEaPeYznEdpIITERH5GIsbfxXdGuj5PrR/RApwzX+BKDvzzSghwKCtQMPrLr23JnMbpODp8iaQONL+ua+c6HoqOBERkY9xnht/nucGAIpzgJ3Tgb++BaACcX2ADi8CYbGy/uTnwNZJQNEZmeDvinuB9rMBo1GGa5/+XlLFS3KAqCuB5vcAdeOdO7c7qeBEREQexlRwHQFV3BAREREATuJHREREtRSLGyIiIgoqLG48qTRf+ra48qSvOA/IOwaUFmuvP78bOLtde13hBWD/G5IMruXcPuD4x3KOykpKgKyNwJlftPd1JxVcVYHibKC0oPr7EhERuYip4J6QsRrY9xJwdou8j2wJtHpIErgVB/XjiVRg+1Sg4ET5stiuQM8VQHgy8GVHIGev7T6RVwKD9wKH5wI7Jlc9ZnQHYOhOYH1v4Mx3lVYagN6fAo2uBz5NAorP2q6OSAaGH3UvFdxSChx6Gzj4OpB/TJY1vA5o82/g8qH6+xIREbmJHYrd7VC8NwXYNVOKGNUaW3Bpfpqmo2U4t70C59BcYOuD2uusw7ctGvlPZefQ+9E5Wq+3awhQt5FrqeAWM/DdzcDJNZcWWNtgBGAGOs8BrtIoyIiIiBxgh+KacGGPFDZAhcIGkC90FTj+IXD8f9r7lhZL+rY9aqlOYWM9hx43ala1RHKnXEkFP7YEOPkpyj6DMpeOtX0KkPub620jIiJygMWNOw6/Uz65nRbFCBx6U3vdvtlwHFDphxylgh98E7r/rBQD8NsCb7SMiIgIAIsb95zfKXcz7FHNcndHy5ktXmlSzdBJBc/eB92iTe8zISIi8gAWN+4wRcBh/pO99OyQCI83p+bopIIbwxzsawCMziSPExERuYbFjTv08pkAeWSVdIv2uit1+tv4O71U8MSb9B/VwQIkjvBGq4iIiACwuHFPs7skp0nzy9wgy1vZGRnUqC8Q3tSbrfMSB6ngrR+G3M3SuKOlmICIZkDSrV5sHxER1XYsbtwRUg/o/y0Qfrm8L0vfVmRd388lrNKeG7YCoXHa6+KuB1r8P/v7NhsL3UdiLabYX2cIAxr2sb8++R7XU8Fjrgb6fAqYwi/tbyrfN6IZ0P8bJx5dERERuY7z3HgiONNSAvzxKZC5XjoYN+gONL1DChyH+1qAY0tlVFVJNhCeBHSYDTTsIevz/gDSBgG5h+V9vWSgzxdAVAt5v/mfwLHFKBt2Hd4UuOEAEBYGFBcDG64BsnfLesUEXP0s0O7fsu2ZX4D0oZIoDgB1E4EB3wH1ki6d241U8JIc4Nhy4NwvgKGOPMZKGAoYOG8kERG5hqngOpgKTkREFHg4iR8RERHVSixuiIiIKKiwuKkpqgUoviDxBV459nnAbCdRPO8P4NT3QKkXzk1ERORn2LvT20rygAOvSEp20WkABulc23ZmeadhVxWfB/a9KDEQJRdkVFOTm4B2M4HLOgJb7gOOvgubjCdTJHD9N0ADO8GXREREAY4dir3ZobgkD/i6D3Bhp22wpmKUP69b6fqEdkVngfU9gbwjlZK7Lw3djmwNZO+yv/+AH90vroiIiGoQOxT7g73PAxd2VUoMhxQjqgXYfJcUQK7YOaNqYWM9tqVYv7ABgG/7u3ZeIiIiP8fixlssJcDh+VWLjzIqUFoAHF9R/WOX5ALH3tM5thPMF4G8313fn4iIyE+xuPGWi1nSD0aPYrKfrq0n7yhg8UDn4BMfu38MIiIiP8PixlvspWbbUJ3crvKxPZSqHdbIM8chIiLyIyxuvCWsAVC/O3Q/YrUUaDKi+seObAFEtoRutpQzmo9xb38iIiI/xOLGm9o9AcCivU4xAg17SQ5VdSkK0O4p2AzxrnxsR2Kvqf55iYiIAgCLG2+6fAhwzYJLqdgG24Ts+tcAvVdLoeKK5DuBji9dOm6l5O5G/YG+X9nft+7lwA0/uXZeIiIiP8d5bmoiOPPiX8CxJUD2fsBUD0i6GYjr63phU1HBH8CRxTIsPCQaaHqb3A1SFKC0FPjxduDPz+URWEgU0PVtoNlt7p+XiIiohjEVXAdTwYmIiAIPJ/EjIiKiWonFDREREQUVFjdEREQUVPyiuJk7dy6Sk5MRFhaGLl264LvvvtPdPj09HV26dEFYWBiaN2+O+fPn11BLiYiIyN/5vLj56KOP8NBDD+Gxxx7Djh070KtXLwwePBgnTpzQ3P7YsWMYMmQIevXqhR07dmDmzJmYNGkSUlNTa7jlRERE5I98Plrq2muvRefOnTFv3ryyZa1bt8aIESOQkpJSZfvp06djzZo12L9/f9my8ePHY9euXdi8ebNT5+RoKSIiosATEKOliouLsW3bNgwcONBm+cCBA/Hjjz9q7rN58+Yq2w8aNAhbt25FSUmJ19pKREREgcHky5OfOXMGZrMZjRrZBjg2atQIWVlZmvtkZWVpbl9aWoozZ84gPj6+yj5FRUUoKioqe5+Tk+OB1hMREZE/8nmfGwBQKs3Uq6pqlWWOttdabpWSkoLo6OiyV2JiopstJiIiIn/l0+KmQYMGMBqNVe7SnDp1qsrdGavGjRtrbm8ymVC/fn3NfWbMmIHs7OyyV0ZGhmcugIiIiPyOTx9L1alTB126dMGGDRswcuTIsuUbNmzA8OHDNffp0aMHPvvsM5tl69evR9euXRESEqK5T2hoKEJDQ8veW+/08PEUERFR4LB+bzscC6X62IcffqiGhISoixYtUvft26c+9NBDakREhPr777+rqqqq//73v9W77767bPujR4+q4eHh6pQpU9R9+/apixYtUkNCQtSVK1c6fc6MjAwVAF988cUXX3zxFYCvjIwM3e95n965AYDRo0fj7NmzeOaZZ5CZmYl27dph7dq1aNq0KQAgMzPTZs6b5ORkrF27FlOmTMHbb7+NhIQEvPHGG7j55pudPmdCQgIyMjIQGRmp27enunJycpCYmIiMjAwOMXcSP7Pq4edVPfy8qo+fWfXw86o+dz4zVVWRm5uLhIQE3e18Ps9NMOH8OdXHz6x6+HlVDz+v6uNnVj38vKqvJj4zvxgtRUREROQpLG6IiIgoqLC48aDQ0FA89dRTNiOzSB8/s+rh51U9/Lyqj59Z9fDzqr6a+MzY54aIiIiCCu/cEBERUVBhcUNERERBhcUNERERBRUWNx6wadMmDBs2DAkJCVAUBZ988omvm+TXUlJS0K1bN0RGRiIuLg4jRozAwYMHfd0svzZv3jy0b98eUVFRiIqKQo8ePbBu3TpfNytgpKSkQFEUPPTQQ75uil96+umnoSiKzatx48a+bpbfO3nyJO666y7Ur18f4eHh6NixI7Zt2+brZvmlZs2aVfk3pigKHnzwQa+cj8WNB+Tn56NDhw546623fN2UgJCeno4HH3wQW7ZswYYNG1BaWoqBAwciPz/f103zW02aNMELL7yArVu3YuvWrbj++usxfPhw7N2719dN83u//PILFixYgPbt2/u6KX6tbdu2yMzMLHvt2bPH103ya+fPn8ff/vY3hISEYN26ddi3bx9eeeUVxMTE+LppfumXX36x+fe1YcMGAMCtt97qlfP5PH4hGAwePBiDBw/2dTMCxpdffmnzfvHixYiLi8O2bdvQu3dvH7XKvw0bNszm/fPPP4958+Zhy5YtaNu2rY9a5f/y8vJw5513YuHChXjuued83Ry/ZjKZeLemGl588UUkJiZi8eLFZcuaNWvmuwb5uYYNG9q8f+GFF3DFFVegT58+Xjkf79yQz2VnZwMAYmNjfdySwGA2m/Hhhx8iPz8fPXr08HVz/NqDDz6IoUOH4u9//7uvm+L3Dh8+jISEBCQnJ+O2227D0aNHfd0kv7ZmzRp07doVt956K+Li4tCpUycsXLjQ180KCMXFxVi+fDnuvfdej+Y7VsTihnxKVVVMnToV1113Hdq1a+fr5vi1PXv2oF69eggNDcX48eOxevVqtGnTxtfN8lsffvghtm/fjpSUFF83xe9de+21WLZsGb766issXLgQWVlZ6NmzJ86ePevrpvmto0ePYt68eWjZsiW++uorjB8/HpMmTcKyZct83TS/98knn+DChQsYN26c187Bx1LkUxMmTMDu3bvx/fff+7opfq9Vq1bYuXMnLly4gNTUVIwdOxbp6ekscDRkZGRg8uTJWL9+PcLCwnzdHL9X8bH61VdfjR49euCKK67A0qVLMXXqVB+2zH9ZLBZ07doVs2fPBgB06tQJe/fuxbx58zBmzBgft86/LVq0CIMHD3aY7O0O3rkhn5k4cSLWrFmDjRs3okmTJr5ujt+rU6cOWrRoga5duyIlJQUdOnTA66+/7utm+aVt27bh1KlT6NKlC0wmE0wmE9LT0/HGG2/AZDLBbDb7uol+LSIiAldffTUOHz7s66b4rfj4+Cq/WLRu3RonTpzwUYsCw/Hjx/H111/jvvvu8+p5eOeGapyqqpg4cSJWr16NtLQ0JCcn+7pJAUlVVRQVFfm6GX6pf//+VUb73HPPPbjqqqswffp0GI1GH7UsMBQVFWH//v3o1auXr5vit/72t79VmcLi0KFDaNq0qY9aFBisA0iGDh3q1fOwuPGAvLw8/Pbbb2Xvjx07hp07dyI2NhZJSUk+bJl/evDBB7FixQp8+umniIyMRFZWFgAgOjoadevW9XHr/NPMmTMxePBgJCYmIjc3Fx9++CHS0tKqjDwjERkZWaUPV0REBOrXr8++XRqmTZuGYcOGISkpCadOncJzzz2HnJwcjB071tdN81tTpkxBz549MXv2bIwaNQo///wzFixYgAULFvi6aX7LYrFg8eLFGDt2LEwmL5cfKrlt48aNKoAqr7Fjx/q6aX5J67MCoC5evNjXTfNb9957r9q0aVO1Tp06asOGDdX+/fur69ev93WzAkqfPn3UyZMn+7oZfmn06NFqfHy8GhISoiYkJKg33XSTunfvXl83y+999tlnart27dTQ0FD1qquuUhcsWODrJvm1r776SgWgHjx40OvnYio4ERERBRV2KCYiIqKgwuKGiIiIggqLGyIiIgoqLG6IiIgoqLC4ISIioqDC4oaIiIiCCosbIiIiCiosboiIiCiosLghIr+WlpYGRVFw4cIFh9suWbIEMTExXm+Ts5o1a4Y5c+b4uhlEtQ6LGyKqEf5WeHhSMF8bUSBicUNERERBhcUNETmlb9++mDBhAiZMmICYmBjUr18fjz/+OKzxdMXFxXj00Udx+eWXIyIiAtdeey3S0tIAyKOle+65B9nZ2VAUBYqi4OmnnwYALF++HF27dkVkZCQaN26MO+64A6dOnfJYuz/77DN06dIFYWFhaN68OWbNmoXS0tKy9Yqi4L///S9GjhyJ8PBwtGzZEmvWrLE5xpo1a9CyZUvUrVsX/fr1w9KlS8seleldGwAUFBTg3nvvRWRkJJKSkpgaTVQTvB7NSURBoU+fPmq9evXUyZMnqwcOHFCXL1+uhoeHlyUh33HHHWrPnj3VTZs2qb/99pv68ssvq6GhoeqhQ4fUoqIidc6cOWpUVJSamZmpZmZmqrm5uaqqquqiRYvUtWvXqkeOHFE3b96sdu/eXR08eHDZeTdu3KgCUM+fP++wjYsXL1ajo6PL3n/55ZdqVFSUumTJEvXIkSPq+vXr1WbNmqlPP/102TYA1CZNmqgrVqxQDx8+rE6aNEmtV6+eevbsWVVVVfXYsWNqSEiIOm3aNPXAgQPqBx98oF5++eVlbdK7tqZNm6qxsbHq22+/rR4+fFhNSUlRDQaDun//fnd/HESkg8UNETmlT58+auvWrVWLxVK2bPr06Wrr1q3V3377TVUURT158qTNPv3791dnzJihqmrVwsOen3/+WQVQViC4U9z06tVLnT17ts027733nhofH1/2HoD6+OOPl73Py8tTFUVR161bV3aN7dq1sznGY489ZtMme9fWtGlT9a677ip7b7FY1Li4OHXevHkOr4WIXGfy4U0jIgow3bt3h6IoZe979OiBV155BVu3boWqqrjyyittti8qKkL9+vV1j7ljxw48/fTT2LlzJ86dOweLxQIAOHHiBNq0aeNWe7dt24ZffvkFzz//fNkys9mMwsJCFBQUIDw8HADQvn37svURERGIjIwsezR28OBBdOvWzea411xzjdNtqHhsRVHQuHFjjz52I6KqWNwQkUcYjUZs27YNRqPRZnm9evXs7pOfn4+BAwdi4MCBWL58ORo2bIgTJ05g0KBBKC4udrtNFosFs2bNwk033VRlXVhYWNnfQ0JCbNYpilJWZKmqalPQWZc5S+/YROQdLG6IyGlbtmyp8r5ly5bo1KkTzGYzTp06hV69emnuW6dOHZjNZptlBw4cwJkzZ/DCCy8gMTERALB161aPtbdz5844ePAgWrRo4fIxrrrqKqxdu9ZmWeU2al0bEfkOR0sRkdMyMjIwdepUHDx4EB988AHefPNNTJ48GVdeeSXuvPNOjBkzBqtWrcKxY8fwyy+/4MUXXywrDJo1a4a8vDx88803OHPmDAoKCpCUlIQ6dergzTffxNGjR7FmzRo8++yzHmvvk08+iWXLluHpp5/G3r17sX//fnz00Ud4/PHHnT7Gv/71Lxw4cADTp0/HoUOH8L///Q9LliwBgLI7OlrXRkS+w+KGiJw2ZswYXLx4Eddccw0efPBBTJw4EQ888AAAYPHixRgzZgwefvhhtGrVCjfeeCN++umnsjsyPXv2xPjx4zF69Gg0bNgQL730Eho2bIglS5bg448/Rps2bfDCCy/gP//5j8faO2jQIHz++efYsGEDunXrhu7du+PVV19F06ZNnT5GcnIyVq5ciVWrVqF9+/aYN28eHnvsMQBAaGio3WsjIt9R1Oo8PCaiWqtv377o2LEj4wQAPP/885g/fz4yMjJ83RQi0sA+N0REDsydOxfdunVD/fr18cMPP+Dll1/GhAkTfN0sIrKDj6WIKGAMHjwY9erV03zNnj3ba+c9fPgwhg8fjjZt2uDZZ5/Fww8/bDMLMRH5Fz6WIqKAcfLkSVy8eFFzXWxsLGJjY2u4RUTkj1jcEBERUVDhYykiIiIKKixuiIiIKKiwuCEiIqKgwuKGiIiIggqLGyIiIgoqLG6IiIgoqLC4ISIioqDC4oaIiIiCyv8H0UFj2cuapUkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "colors = {'Iris-setosa':'orange','Iris-versicolor':'red','Iris-virginica':'blue'}\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(iris['sepal_length'], iris['sepal_width'], c=iris['class'].map(colors))\n",
    "plt.xlabel('sepal_length')\n",
    "plt.ylabel('sepal_width')\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(iris['petal_length'], iris['petal_width'], c=iris['class'].map(colors))\n",
    "plt.xlabel('petal_length')\n",
    "plt.ylabel('petal_width')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37b64e89",
   "metadata": {},
   "source": [
    "From the plots above, we can see that the iris virginica (blue) has a fairly unique combination of sepal and petal features that will allow us to distinquish it from the other two species. Thus, if we only have the information of the sepal and petal features, we would be able to classify an unknown iris as being an iris virginica or not. We have therefore posed a binomal classification problem, since we consider only two classes - iris virginica or \"not iris virginica\". For our ML model, we need to provide numerical class labels, so we will label iris virginica as \"1\" (positive case, or event), and the other two iris types as \"0\" (negative case, or non-event)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ca1b0c2f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   sepal_length  sepal_width  petal_length  petal_width  class\n",
      "0           5.1          3.5           1.4          0.2      0\n",
      "1           4.9          3.0           1.4          0.2      0\n",
      "2           4.7          3.2           1.3          0.2      0\n",
      "3           4.6          3.1           1.5          0.2      0\n",
      "4           5.0          3.6           1.4          0.2      0\n",
      "       sepal_length  sepal_width  petal_length  petal_width       class\n",
      "count    150.000000   150.000000    150.000000   150.000000  150.000000\n",
      "mean       5.843333     3.054000      3.758667     1.198667    0.333333\n",
      "std        0.828066     0.433594      1.764420     0.763161    0.472984\n",
      "min        4.300000     2.000000      1.000000     0.100000    0.000000\n",
      "25%        5.100000     2.800000      1.600000     0.300000    0.000000\n",
      "50%        5.800000     3.000000      4.350000     1.300000    0.000000\n",
      "75%        6.400000     3.300000      5.100000     1.800000    1.000000\n",
      "max        7.900000     4.400000      6.900000     2.500000    1.000000\n"
     ]
    }
   ],
   "source": [
    "#replacing class in the data set for setosa and versicolor to zero and virginicia to one\n",
    "iris = iris.replace('Iris-setosa', 0)\n",
    "iris = iris.replace('Iris-versicolor', 0)\n",
    "iris = iris.replace('Iris-virginica', 1)\n",
    "iris['class'] = iris['class'].astype(int)\n",
    "\n",
    "#print the first five lines in the new edited data set\n",
    "print(iris.head())\n",
    "#Generate descriptive statistics for the new data set and print them\n",
    "print(iris.describe())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b499c8f",
   "metadata": {},
   "source": [
    "We will now build a logistic regression classifier to identify whether an unknown iris is of the virginica species or not based on known sepal and petal features. The first step is to partition our data in a training set with which to develop the ML model, and a separate test set with which to independently assess the performance of our trained model. We set the size of the training set to 2/3 of the data (100), and the test set to 1/3 of the data (50)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3731c3f5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 4)\n",
      "(50, 4)\n",
      "(100, 1)\n",
      "(50, 1)\n"
     ]
    }
   ],
   "source": [
    "#Load module\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# define the predictors and class\n",
    "X = iris[['sepal_length','sepal_width','petal_length','petal_width']]\n",
    "y = iris[['class']]\n",
    "\n",
    "#Split data set into random train and test subsets\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)\n",
    "#print the shape of the arrays\n",
    "print(X_train.shape)\n",
    "print(X_test.shape)\n",
    "print(y_train.shape)\n",
    "print(y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c69ff63",
   "metadata": {},
   "source": [
    "We are now ready to build a logisitic regression model for this classification task. This model has the form:\n",
    "\n",
    "$P\\left( y = 1 \\lvert g(x) \\right) = \\dfrac{1}{1+ e^{-g(x)}}$\n",
    "\n",
    "$g(x) = b_0 + b_1 x_1 + b_2 x_2 + ... + b_k x_k$\n",
    "\n",
    "for $k$ predictors. In our case, we have *sepal_length* (x1), *sepal width* (x2), *petal length* (x3), and *petal width* (x4). As shown in the figure below, this simple linear classification model uses the linear input $g(x)$, transformed by a kernel function (logit function), to produce a probability between 0 and 1 of a positive event. The more positive the linear function $g(x)$ for a combination of predictors and their coefficients, the higher the predicted probability of an event (here, the species being virginica). Conversely, the more negative the function $g(x)$, the lower the probability of an event (not-virginica). The coefficients $b_0..b_k$ are fit using the labelled observations of iris species (virginica or not-virginica).\n",
    "\n",
    "<img src=\"img/logit.png\" alt=\"Title\" style=\"width: 400px;\"/>\n",
    "\n",
    "For building this model, we use Python's statsmodels package. The data is first scaled by column so that all the features are of the same order of magnitude. This improves the performance of the model. We then train the model, using the training data partition. Statsmodels provides a Logit() function for performing logistic regression. The Logit() function accepts y and X as parameters and returns the Logit object. The model is then fitted to the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a849b7c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.320003\n",
      "         Iterations 7\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                  class   No. Observations:                  100\n",
      "Model:                          Logit   Df Residuals:                       96\n",
      "Method:                           MLE   Df Model:                            3\n",
      "Date:                Wed, 21 Dec 2022   Pseudo R-squ.:                  0.5008\n",
      "Time:                        10:48:52   Log-Likelihood:                -32.000\n",
      "converged:                       True   LL-Null:                       -64.104\n",
      "Covariance Type:            nonrobust   LLR p-value:                 7.415e-14\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             0.3173      0.782      0.406      0.685      -1.216       1.851\n",
      "x2             0.9130      0.509      1.794      0.073      -0.085       1.911\n",
      "x3            -1.6096      1.845     -0.872      0.383      -5.226       2.007\n",
      "x4             4.5741      1.644      2.782      0.005       1.352       7.796\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "#import moldules\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.tools as st\n",
    "from sklearn import preprocessing\n",
    "\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "from sklearn.metrics import roc_auc_score\n",
    "\n",
    "#Standardize features by removing the mean and scaling to unit variance\n",
    "scaler = preprocessing.StandardScaler().fit(X)\n",
    "X_train_scaled = scaler.transform(X_train)\n",
    "X_test_scaled = scaler.transform(X_test)\n",
    "\n",
    "#training model\n",
    "logit_model1=sm.Logit(y_train,X_train_scaled).fit()\n",
    "print(logit_model1.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63698445",
   "metadata": {},
   "source": [
    "From the summary model results we see that two of the predictors, sepal length (x1) and petal length (x3) are not statistically significant, seen from their p-values >> 0.05. In addition, the sepal width (x2) is marginally significant at p-value = 0.073. This situation is known as multi-colinearity, caused by having features that carry similar, redundant information. It can be remedied by removing some features. We will start by dropping the two least significant ones (sepal length and petal length), and refitting the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3e0b5059",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.324486\n",
      "         Iterations 7\n",
      "                           Logit Regression Results                           \n",
      "==============================================================================\n",
      "Dep. Variable:                  class   No. Observations:                  100\n",
      "Model:                          Logit   Df Residuals:                       98\n",
      "Method:                           MLE   Df Model:                            1\n",
      "Date:                Wed, 21 Dec 2022   Pseudo R-squ.:                  0.4938\n",
      "Time:                        10:49:37   Log-Likelihood:                -32.449\n",
      "converged:                       True   LL-Null:                       -64.104\n",
      "Covariance Type:            nonrobust   LLR p-value:                 1.766e-15\n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "x1             1.1172      0.399      2.797      0.005       0.334       1.900\n",
      "x2             3.2190      0.624      5.161      0.000       1.997       4.442\n",
      "==============================================================================\n"
     ]
    }
   ],
   "source": [
    "X_train_scaled = X_train_scaled[:,[1,3]]\n",
    "X_test_scaled = X_test_scaled[:,[1,3]]\n",
    "\n",
    "logit_model1=sm.Logit(y_train,X_train_scaled).fit()\n",
    "print(logit_model1.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eddea1f6",
   "metadata": {},
   "source": [
    "We see that dropping two non-significant features improved the model fit. The sepal width (now x1) is now statistically significant (p-value < 0.05), as well as the other retained predictor (petal width, now x2). Since petal width (x2) has the largest positive coefficient, it has the greatest importance in determining the iris class. We can confirm this by referring back to the data visualization above, in which the iris virginica is separated most from the other species along the dimension of petal width.\n",
    "\n",
    "We now proceed with evaluating our model, by computing the model predictions for the unseen test partition ($\\mathrm{\\hat{y}})$. The predictions $\\mathrm{\\hat{y}}$ are the probability of the iris being of the virginica species. When comparing to the actual observed values, we see that there is a good correlation between high predicted probabilities and actual positive cases of virginica.\n",
    "\n",
    "To make binary predictions, we have to set a probability threshold between 0 and 1 above which we call the prediction for the virginica species. A good way to assess the quality of the classifier is to consider its performance across the full range of thresholds [0,1]. For each threshold, we compute the true positive rate (TPR, or recall) and the false positive rate (FPR, or the probability of false alarms). The resulting figure is the so-called Receiver Operator Curve (ROC). The line traces for each threshold level the resulting combination of the TPR and the FPR. It is important to realize that these metrics are in opposition - achieving a high recall usually comes at the cost of having a large number of false alarms (high FPR). A perfect classifier would identify all instances of iris virginica (recall = 1.0) and never mistakenly identify a species as being virginica when it is not (FPR = 0.0). Such a perfect ROC curve would go through the top left-hand corner. Conversely, a classifier with no skill would lie on the diagonal. This performance is quantified by computing the Area Under the Curve (AUC), with a perfect classifier having an AUC = 1.0, and one with no skill having AUC = 0.5. From the results below, we can see that our classifier performs very well, with a ROC curve that runs close to the top left-hand corner, and a near-perfect AUC = 0.969."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49750904",
   "metadata": {},
   "outputs": [],
   "source": [
    "yhat = logit_model1.predict(X_test_scaled)\n",
    "print('Predicted probability of iris being of the species virginica:')\n",
    "print(yhat[:10])\n",
    "print('\\nActually occurrence of iris virginica:')\n",
    "print(y_test[:10].to_numpy().reshape(10))\n",
    "\n",
    "# Compute ROC curve and area under curve (AUC)\n",
    "n_outputs = 1\n",
    "fpr = dict()\n",
    "tpr = dict()\n",
    "roc_auc = dict()\n",
    "fpr, tpr, _ = roc_curve(y_test, yhat)\n",
    "roc_auc = auc(fpr, tpr)\n",
    "\n",
    "plt.figure(figsize=[4,3])\n",
    "lw = 2\n",
    "plt.plot(fpr, tpr, color='darkorange', lw=lw, label='ROC curve (area = %0.2f)' % roc_auc)\n",
    "plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')\n",
    "plt.text(0.6, 0.15, 'AUC ='+'{0:.3f}'.format(roc_auc), fontsize=12)\n",
    "plt.xlim([0.0, 1.0])\n",
    "plt.ylim([0.0, 1.05])\n",
    "plt.xlabel('False Positive Rate')\n",
    "plt.ylabel('True Positive Rate')\n",
    "plt.title('ROC curve')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "953827ee",
   "metadata": {},
   "source": [
    "Having established that our classifier performs well, we can now select our probability threshold for making binary predictions. A reasonable first guess is to set it at P = 0.5. We then compute a confusion matrix, which quantifies the binary prediction in terms of the number of *true positives* (prediced virginica, actual virginica), *true negatives* (predicted non-virginica, actual non-virginica), *false positives* (predicted virginica, actual non-virginica) and *false negatives* (predicted non-virginica, actual virginica). The larger the numbers on the diagonal (true positives and true negatives), and the smaller the numbers off the diagonal (false negatives and false postives), the better. We can see that with a threshold of 0.5 our model performs fairly well in this regard. These metrics can also be used to compute the recall and precision. As with recall and the FPR, there is an adverserial relationship between recall and precision. Our choice of threshold yields perfect recall of 1.0, but at the cost of too many false positives, leading to a relatively low precision (0.615). We can thus conclude that our threshold is too low, so that cases are too easily classified as being virginica (also when they are not). A way to measure the balance of performance in terms of these two metrics is to compute their harmonic mean, known as the F1 score. Here we achieve an F1 score of 0.762."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ee30293b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAADyCAYAAAAcCUOdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAApdUlEQVR4nO3deXhN1/rA8e85kUFIIiISV4qIJioIKuYhEkqv4Yqh3IbWUKWGFBdBWtOPoMZExBiqairSW4qiarg0DYJWDTUPqUpMSSRIIjm/P1zn9jSR8Zyc7Jz302c/j6y199qv8zSvtdZeZ22VRqPRIIQQCqM2dgBCCFEYkryEEIokyUsIoUiSvIQQiiTJSwihSJK8hBCKJMlLCKFIkryEEIokyUuYPFmnrUySvIrR2bNnGT9+PD4+PtSvXx8/Pz8++eQTbt++bbB77t69m3bt2lGvXj2mTJmit3Y9PDxYsmSJ3tozlgMHDhAUFJTneRMnTsTX17cYIhL5pZKvBxWPDRs2EBISQtOmTfH396dy5crcunWL1atX8+jRI9auXYunp6fe79u0aVNq1KjB2LFjcXJyokaNGnpp98yZMzg7O+Ps7KyX9oylf//+AKxfvz7X827dukVKSgp16tQpjrBEPkjyKgaxsbH079+fgIAAgoODdeoePnxIjx49sLW1ZceOHXq/t4eHByNGjCAwMFDvbZcG+U1eouSRYWMxiIyMxMbGhrFjx2arq1ixIhMnTuStt94iJSVFW75792569OhBw4YNadmyJVOmTCEpKUlbv2TJEjp06MChQ4fo2rUrdevWpWPHjnz99dcAxMTE4OHhAcDSpUvx8PAgLi4ux+FPXFwcHh4eREVFacvWr19Pp06dqFevHq1bt2batGk68f112JiQkMCkSZNo27Yt9evXp1evXhw4cEDnPh4eHmzYsIHg4GCaNGlCw4YNCQwM5P79+6/87F7GtnfvXoYPH06DBg1o0aIFERERpKSkMHnyZN58801atGjBvHnzdOav4uLimDBhAq1atcLT05PmzZszYcIEHj16BLxIXMePH+f48eN4eHgQExOj/dw2b95Mu3btaNGiBUePHtX53A4cOJDt73/9+nUaNGiQryGo0A9JXgam0Wg4evQozZs3p2zZsjme06lTJ0aOHEn58uUBiIiIYMyYMXh5eREWFsaIESPYu3cv/fv359mzZ9rr7t27x4wZM3jvvfdYuXIlLi4uTJw4katXr+Lp6cmWLVsA6NWrF1u2bKFy5cr5innXrl3MnTuXgIAAIiMjGTFiBN988w0zZ87M8fz79+/Tq1cvjh8/zpgxY1iyZAlVq1ZlxIgR2XqTixYtIisri4ULFzJhwgQOHTpESEhInjEFBwfj7u7OsmXLaNasGaGhofTq1QsrKytCQ0Px9fVl9erVfPfddwA8ffqU9957j6tXrzJ16lQiIyPp168f3377LQsXLgRg6tSp1KlThzp16rBlyxadYfuiRYsICgoiKCiIBg0a6MTi5+dH9+7dWbFiBVevXiUzM5OJEydSsWJFPv3003x9xqLoyhg7gNLu0aNHpKWl4eLikq/zk5KSWLZsGb1792bq1Knacnd3dwICAoiKiuLdd98FXvyCzpo1i+bNmwNQo0YN2rVrx+HDhxk0aJD2l87Z2TnbL2BuYmJiqFq1KgEBAajVapo0aYK1tbW2x/JXa9eu5eHDh+zZs4fXXnsNgLZt2zJgwAA+++wzunTpglqt1v49Zs+erb32l19+0Sac3LRu3ZrRo0cDUKtWLXbt2oWDg4P2IUTLli3Zs2cPp06d4u233+bGjRs4OzszZ84cqlWrBkCzZs04e/Ysx48f17bz8h+Mv34+ffv2pVOnTq+MJzg4mOjoaKZNm0br1q05e/YsX3zxhbY9YXjS8zKwl7+0mZmZ+Tr/zJkzpKen07VrV53yxo0bU7VqVWJiYnTK//xL93Ly/MmTJ0WI+MUv+Y0bN+jRowcRERGcP3+erl278v777+d4/vHjx2nYsKE2cb3UrVs37t27x7Vr13KM92XMT58+zTOmhg0bav/s6OgIgJeXl7ZMpVJhZ2fH48ePAXjjjTfYuHEjLi4u3L59m//85z+sWbOGa9eukZGRkef9Xg65X8XW1paZM2dy/PhxFi1axJAhQ2jcuHGe7Qr9keRlYBUqVKBcuXLcuXPnlec8efKExMREAO28VqVKlbKdV6lSJe0v50t/Hoq+TJRFfQbz97//nQULFmBtbU14eDj+/v74+fmxa9euHM9PSkp6ZbwAycnJOcb7Mub8xJtTj+ZVw/CX1q5dS4sWLWjfvj0TJ07kp59+yvOalxwcHPI8p3nz5lSpUoWsrCxZRmEEkryKQatWrYiJiSEtLS3H+qioKJo3b87p06exs7MDyHES+969e9jb2xcpFpVKla0XmFNPrUuXLmzcuJGYmBgWL15MhQoVGD9+PPHx8dnOtbOze2W8QJFjLoydO3cyZ84cBg0aRHR0NMeOHWPlypV6WyoCLx6E3L9/Hzc3Nz755BPS09P11rbImySvYjBo0CASExNZtGhRtroHDx6wevVqqlevToMGDfDy8sLCwoKdO3fqnHfy5Enu3LlDo0aNihRLuXLltPNwL506dUrnnNGjRzNy5EgAbGxsePvttxk+fDiZmZkkJCRka9Pb25vTp09nW2y7Y8cOHB0dqV69epFiLozY2FhsbGz48MMPqVixIgCpqanExsaSlZWlPe9lb7Wgzp49y6pVqxg2bBgLFizg2rVrpWLRrpLIhH0xaNCgAR9//DGLFy/m6tWr+Pv7Y29vz+XLl1mzZg2pqamsXLkSlUpFhQoV+PDDDwkPD8fc3Bw/Pz/i4uIIDQ2lVq1a9OjRo0ixtGvXjvXr1zN58mR69+6tjcHMzEx7TrNmzZg6dSpz586lTZs2JCcnEx4eTo0aNahdu3a2NgcOHMiOHTsYOHAgI0eOxN7enn//+9/89NNPhISEFDpBFEX9+vXZtGkTc+bMoV27diQkJBAZGcn9+/e1vVt4MXd1+vRpoqOj870ANT09nYkTJ+Lq6sqHH36IhYUF7733HpGRkbRv315nLk4YjiSvYvLRRx9Rp04dNmzYwOzZs0lMTMTZ2Zk2bdowbNgw/va3v2nPHTVqFJUqVeLLL79k69atVKhQgU6dOjF69Oh8z9m8SsuWLQkKCmL9+vXs27cPT09PwsPD6du3r/acvn37kpGRwebNm9m4cSNWVlY0b96c8ePHY25unq1NR0dHNm3axIIFC5g1axYZGRnUrl2biIgI/Pz8ihRvYfn7+xMXF8f27dvZuHEjTk5OtG3blnfffZdPP/2UK1euUKtWLQICAvj1118ZMmQIs2fPztdyksWLF3Pt2jU2bdqEhYUFAIGBgezbt4+goCC++eYbLC0tDf1XNHmywl4IoUgy5yWEUCRJXkIIRZLkJYRQJEleQghFkuQlhFAkSV5CCEWS5CWEUCSTWaTabM5hY4dgcg6Na2vsEEyOVQF+o8s2HJlr/dPT4UWMxrBMJnkJIf5CbZb3OSWYJC8hTJVK2bNGkryEMFXS8xJCKJIkLyGEIqlUxo6gSCR5CWGqpOclhFAkSV5CCEWSp41CCEUyk56XEEKJZNgohFAkGTYKIRRJel5CCEWS5CWEUCRZpCqEUCTpeQkhFEmt7F9/ZUcvhCg8GTYKIRRJho1CCEVS+DovZUcvhCg0lVqd61EUiYmJTJkyhTZt2tCoUSP++c9/cvLkSW39hQsX6NevHw0aNMDHx4fIyMgC30OSlxAmSqVW5XoUxdixY/n5559ZuHAh27Ztw9PTk8GDB3P16lUePXrEwIEDqVGjBtu3b2fUqFGEhoayffv2At1Dho1CmCiVgSbsb968ybFjx9i0aRONGjUCIDg4mCNHjvDtt99iZWWFhYUF06ZNo0yZMri5uXHz5k1WrVpFz549830f6XkJYaLUanWuR2HZ29uzcuVK6tatqy1TqVRoNBqSkpI4efIk3t7elCnzv75Ts2bNuH79Og8ePMh//IWOUAihaIYaNtra2tK2bVssLCy0ZXv27OHWrVu0atWKu3fv4uzsrHNN5cqVAbhz506+7yPDRiFMVF7DRj8/v1zrDxw4kK/7xMbGMnnyZPz8/PD19WX27Nk6iQ3A0tISgLS0tHy1CdLzEsJkGWrY+Gfff/89gwcPpn79+ixcuBAAKysr0tPTdc57mbSsra3z3bb0vIQwUXkNDfPbs3qVL7/8klmzZtGhQwfmz5+v7W05OzuTkJCgc+7Ln52cnPLdvvS8hDBRKpUq16MoNm7cyP/93/8REBDA4sWLdYaJ3t7exMbGkpmZqS2Ljo7G1dUVBweHfN9DkpcQJspQw8br168TEhJChw4dGDp0KA8ePODevXvcu3ePx48f07NnT1JSUggODubKlStERUWxbt06hg4dWqD7yLBRCBNV1IWor7J3714yMjLYv38/+/fv16nz9/dnzpw5rF69mlmzZuHv74+joyMTJkzA39+/QPdRaTQajT4DL6mazTls7BDypbKNJRsGNyYo6ldO3UrSljuWt2Bku5o0q1mRMmoV5/94zJKD17gUn2LEaHN3aFxbY4dQIMf+c4TwJYu5dvUq9vYV6d2nL4M++NBgizkNwaoA3RHnIdtyrb+7qlcRozEsGTaWIM62loT1rY/NX/4PtLYwY1lAAzycyjPnu0tM2XEBawszwvrUx6GcxStaEwVx5vQpAkcOx7WmGwsXL6FL124sCV3E6pXLjR2awRTH00ZDkmFjCaAC/l7PiUBftxzr+3q7UKGsOX1WneBB6otHzBfuPubzAW/SqJod+y/cK8ZoS6flEUvxqF2bkDnzAGjZug0Zz5+zZvVK+r8/ECsrKyNHqH+GGjYWl5KfXk1ArcrlmNDRnd1n45m282K2+nYelfjht3vaxAXwMDWDbkt/ksSlB+np6Zw8EYNf+7d0yju81ZEnT55wKvbkK65UNkM+bSwO0vMqAeKT0+i1IoZ7j9NpVM1Op85MrcLVwZrvfo3nw9Y16OblTIWy5vzyezIL9l/h6r1UI0VdesTdvk1GRgbVa9TQKa9WrToAN2/coEXLVkaIzLCk5yWKLPnZc+49Ts+xztaqDGXM1PT1duHNahUI2XOJT765gF1ZcyLe9cKxvMx5FdXjx8kAlC9fXqfculw5AFJTS+5DkaJQ+pxXyY/QxJUx+9+/jqO/+oUfrz7k0KX7jN16lrLmZvR6s6oRoysdsrKygFd/10+l8B1HX0WGjcKgnqS/WIV86lYiTzOytOXxyWncePAE98rlX3WpyCcbW1sAUlJ0e1hPUl8MyW1sSudnrPRho1GTV//+/fOd4b/44gsDR1MypaZl8jA1HYsy2f/1L6NWkfY8M4erREG89lo1zMzMuH3rpk75rf/+XNOtljHCMjglDA1zY9TomzdvzokTJ3jw4AFVq1bN9TBl0dce4l3dHruy//u3plrFslRzsOZMXFIuV4r8sLS0pNGbjTnw/X7+vGZ7/7692NjaUrdefSNGZzgqVe5HSWfUntfw4cOxtrYmLCyMFStW4OLiYsxwSqzIozdp83olwvrUJ/LYTcqoVXzU1pWE5Gfs+PmuscMrFYYM/YihHwxk/NiP6d6jJ2dOn2bd2khGjx1XKtd4AagVPmw0er9xwIABNGrUiMWLFxs7lBLrTtIzhqw/zb2UdKZ2qc2ktz24lJDK0A1ntHNiomiaNmvOgsVLuHHjOqNHjWD3rp2MGTeBAYM+MHZoBqNWq3I9SroS8d3G+Ph4zp8/T7t27Qx2D6V8t7E0Udp3G0uDgny3sc7kfbnWnw95K9d6YysRTxudnJwKtAmZEKLolNC7yk2JSF5CiOInyUsIoUhKWIiaG0leQpgo6XkJIRRJkpcQQpEUPmqU5CWEqZKelxBCkSR5CSEUSZ42CiEUqbh6XhEREURHR7N+/Xpt2aRJk4iKitI5z8nJiSNHjuS7XUleQpio4khen3/+OWFhYXh7e+uU//bbbwwbNox+/fppy8zMzArUtiQvIUyUIYeN8fHxBAcHExsbi6urq05dZmYmV65cYfjw4Tg6Ohb6HvlKXidOnChQo3/NskKIkseQPa9z585hZ2fHjh07WLp0Kb///ru27saNG6SlpeHmlvOr/vIrX8krvzueajQaVCoVFy5cKFJQQgjDyyt5+fn55Vp/4MCBV9b5+vri6+ubY92lS5dQqVSsW7eOI0eOoFaradu2LaNHj8bGxibvwP8rX8nLVLdgFqI0M9bDxsuXL6NWq6latSrLly/n5s2bzJ07l0uXLrFu3bp8b0+dr+TVpEmTIgUrhCh5zPLoee3LpWdVFKNGjWLAgAHY/vfFJ+7u7jg6OtKnTx/Onj2Ll5dXvtop1IT9mTNnOH78OBkZGdo9vzUaDU+ePCE2NpavvvqqMM0KIYpRXsnLUFQqlTZxveTu7g7A3bt3DZe8NmzYwMyZM8lpA1a1Wk2rVqXvzcJClEbGWqT6r3/9i8TERCIjI7VlZ8+eBaBWrfy/qanAe9h/+eWXtGrVipiYGAYPHsw777zDmTNnCA0NxdLSkm7duhW0SSGEEahVqlwPQ+nSpQvHjh1j2bJl3Lp1i8OHDzN58mS6dOlSoCeQBU5ecXFx9OvXDzs7O+rVq0dsbCxWVlZ07NiRoUOHyuS+EAphplblehhKu3btCA0NZd++fXTt2pXg4GDeeustQkJCCtROgYeN5ubm2ldB1ahRg5s3b5KRkYG5uTmNGjVizZo1BW1SCGEExTVqnDNnTrayjh070rFjxyK1W+Ce1xtvvMHBgwcBqF69OllZWZw5cwZ4MdkmhFAGY/W89KXAPa+BAwcycuRIkpKSmD17Nn5+fkyYMIGOHTuyc+dO3nzzTUPEKYTQM6VviVPgnlf79u1Zvny59qnAjBkzcHV1ZfPmzdSsWZMpU6boPUghhP4Za8JeXwq1zsvHxwcfHx8A7O3tZZ5LCAVSQoLKTYGT1507d/I8529/+1uhghFCFB8lzGvlpsDJy9fXN8/FbfLFbCFKPoV3vAqevEJCQrIlr5dfC/rpp58KvFZDCGEcJtfz6tGjR47lAQEBzJ07l507d2rnw4QQJZfS57wK/LQxNz4+Phw6dEifTQohDMQknza+ypkzZyhTRnaWFkIJlL7Oq8CZZtKkSdnKsrKy+OOPPzh58iS9evXSS2BCCMMyuTmvmJiYbGUqlYry5cszZMgQhg0bppfAhBCGpYCRYa4KnLx++OEHQ8RhcO+0eM3YIZgce++Rxg7B5Dw9HZ7vc80Unr0KPGE/adIkbt++nWPdtWvXpOclhEKYxBez/7yq/uuvv6Z9+/Y5viDyyJEj/Pjjj/qLTghhMArIT7nKV/KaMWMGhw8fBl7Mb40cmfNwQKPR0LJlS/1FJ4QwGCX0rnKTr+Q1ffp0fvzxRzQaDZMnT+ajjz6iWrVqOueo1WpsbW1p2rSpQQIVQuiXSSQvJycn/P39gRc9Lx8fH8qXL4+5uTkAT58+JS0tjQoVKhgsUCGEful1hboRFDj+zp07s2DBAt555x1t2enTp2nVqhWzZs0iMzNTrwEKIQxD6RP2BU5eYWFh7N69m+7du2vLPD09CQoK4uuvv2bVqlX6jE8IYSBm6tyPkq7AIe7atYugoCDef/99bZmdnR39+/dnzJgxbN++Xa8BCiEMo7i+2xgREUH//v11yi5cuEC/fv1o0KABPj4+Ou9wzK8CJ69Hjx7h4uKSY52rqyvx8fEFDkIIUfyKo+f1+eefExYWplP26NEjBg4cSI0aNdi+fTujRo0iNDS0wB2fAq+wd3NzY+/evTkuidi/fz/Vq1cvaJNCCCMw5Ar7+Ph4goODiY2NxdXVVafuq6++wsLCgmnTplGmTBnc3Ny4efMmq1atomfPnvm+R4GT16BBg7Sv627fvj0ODg48fPiQ77//nn379jF79uyCNimEMAJDzsmfO3cOOzs7duzYwdKlS/n999+1dSdPnsTb21tnB5pmzZqxYsUKHjx4gIODQ77uUeDk1blzZx4/fkx4eDj79u3Tltvb2zNlyhS6dOlS0CaFEEaQ1xNFPz+/XOsPHDjwyjpfX198fX1zrLt79y7u7u46ZZUrVwZefJvHYMkLoG/fvvTp04fr16+TmJiIra0tlpaWbN26FR8fH44ePVqYZoUQxchYyyGePXuGhYWFTpmlpSUAaWlp+W6n0DsHqlQqXF1dOXjwIMuXL+fYsWNkZmZSs2bNwjYphChGeeWu3HpWRWFlZUV6erpO2cukZW1tne92CpW8EhIS2Lp1K9u2bePu3bvY2trSp08funfvTv369QvTpBCimBlrSxxnZ2cSEhJ0yl7+7OTklO92CpS8jh07xubNmzl48CAajYamTZty9+5dwsPD8fb2LkhTQggjM9ZCVG9vbzZv3kxmZqZ2d5ro6GhcXV3zPd8F+VzntXr1at566y0GDx7MlStXCAwM5NChQyxevBiNRlO4v4EQwqiM9QKOnj17kpKSQnBwMFeuXCEqKop169YxdOjQArWTr57X/Pnz8fDwYP369To9rMePHxcsaiFEiWGsCXsHBwdWr17NrFmz8Pf3x9HRkQkTJmg3f8ivfCWvbt26sX//fj744AOaNWtG9+7d83yMKoQo2YprymvOnDnZyurXr8+WLVuK1G6+ktdnn31Gamoq3377LVFRUYwZMwY7Ozv8/PxQqVTZ3qAthCj5TGYP+3LlytGnTx+2bNnCrl276NGjB0eOHEGj0RAUFMSiRYu4dOmSIWMVQuiRmUqV61HSFep5g5ubG0FBQRw+fJjw8HBef/11IiMj+cc//kG3bt30HaMQwgBUeRwlXZFeb21mZkb79u1p3749Dx48ICoqin//+996Ck0IYUhKf2O23lZ6ODg4MGTIEHbt2qWvJoUQBqT0YWORel5CCOUq+ekpd5K8hDBRSuhd5UaSlxAmypCr6IuDJC8hTJTCc5ckLyFMlQwbhRCKpPRhowLezmZaUh7eY21gL+789otOeeLdOPaETWFtYE8+H/0Ohz5fRNqTFCNFWToM6tGS41smce/YAs7vnMa8cT2xKWeV47kj/unD09PhVKtSsZijNBw1qlyPkk56XiXI4wfx7F78CelPU3XK056k8O2CiVhXcKDdoPE8TX7ET9vXkProHp3HhBgpWmUb+357po/syqIvDnDw+G+4vebIlOGdqVOrCp2Hheuc61bNkRmjSt83R9QK77pI8ioBNFlZ/Bb9PT9tXZ1j/flDu0h7kkLPKeGUtakAQDn7SuwJm8Ifl3+lyut1izFa5VOpVIwb9Bartx9jypIdAByM+Y2HialsmDeYRnWqcer8LeDFKvTVM/rzMCkV67IWuTWrOEqf81J47i0dHsRd5+iX4bg3b4/v4HHZ6m+fi8X59braxAXwmuebmFuV5dbZk8UYaelgW86KzbtP8NUe3c/u8q0XWxHXdKmkLRvznh+VK9owf+3+Yo2xOKjy+K+kk55XCVDeoTJ9Z0VSvqJjtrkugMQ/buPm3UanTKVWY1PJmaT4uOIKs9RISnnK2Llbs5X/w9cLgHNX7gDwRk1ngof+nW4jI6jxt/xvT6wUSp+wN2ryev78Ofv27ePkyZPcuXOH9PR0ypYti7OzM40bN6ZDhw46L6YsrazK2UA5m1fWpz1Nwdwq+1tVLKzKkvHsiSFDMxnNvFz514AO7PjhZy5cu4uZmZpV//cen/87mqOxV0pl8pJhYyHdunWLzp07M3nyZC5evIiVlRWOjo6Ym5tz4cIFJk2aRLdu3bhz546xQiw5NOS44aNGA6hk5F9ULRu68fWSj7gWd59h0zcAEDS4I/Y2Zfkk9BsjR2c4KlXuR0lntG7N9OnTcXFxYdu2bdjYZO91JCcnM2bMGGbMmMHy5cuNEGHJYVHWmvQcelgZaU8pZ18phytEfvXu+CYrp/fj0s14ug1fyqPkJ3h5uDBh8Ft0H7WMtIznmJmptdvHvPxzVpbyXzyj9J6X0ZJXbGwsW7ZsyTFxAdja2jJ+/HgCAgKKObKSp4KzC8kJuj1QTVYWj+/fxbVhSyNFpXxj3vNj5sf/4Oipq/Qes4LklGcAdPGpj6WFOXtWBGa75vzOaRw5eZmOQ0KLO1y9U3juMl7ysrW1JSEhAQ8Pj1eec+fOHayscl40aEpc6jTizN5tPH2cqH3iePtcLBnPnuLi2ci4wSnU4J4tCRnjz7a9sQz65Asynmdq69ZEHWPPf37VOf/t1nX5ZNjf6fnxci7fTPhrc4pkyCeKv//+O76+vtnKZ86cSe/evfVyD6Mlr169ejFp0iQCAwNp2rQpVapUwcLCgvT0dOLj4zl+/Djz58+nV69exgqxxKjj04Vff9jBroXBvNn1XZ6lPiZmWySv1W2Ms9sbxg5PcZwcbPjsXz25eecByzYfpuEbr+nUX4u7r13n9VIdtyoA/Hr5Drf+eFhssRqSITdS/e2337C0tOT777/Xma991UirMIyWvEaNGoVarWbu3Lk8eZJ9PqdcuXIEBATw8ccfGyG6kqWsjR1dx83lxy0r+CFyHuaWZanZuDXNen9g7NAUqWMrT6zLWlC9rAMH1o7NVj9kynq+3BljhMiKlyHf+nXp0iVcXV2pXLmywe6h0hj5ldcZGRlcuHCB+Ph4nj59ipWVFc7OztSuXRsLC/2taF545Jre2hL5E/zxQmOHYHKeng7P+6T/ir2RnGv9mzVsCx3HmDFjUKvVLFiwoNBt5MXoi6jMzc2pX7++scMQwuTk1fHK68XSBw4ceGXdpUuXcHR05N133+XGjRtUr16d4cOH07p168KEmiNZJCSEiVKrVLkehZWens6NGzdISUlh9OjRrFy5knr16jFkyBCio6P1Fr/Re15CCOPIKz3l1rPKjYWFBSdOnKBMmTLaqZ+6dety9epVIiMjad68eaHa/SvpeQlholQqVa5HUVhbW2ebs3Z3dyc+Pr5I7f6ZJC8hTJRalftRWBcvXqRhw4acPKm7a8evv/5KrVq1ihj1/8iwUQhTZaCVEu7u7rz++utMnz6dqVOnYm9vz1dffcWZM2fYtm2b3u4jyUsIE2WoLXHUajXLly9n/vz5jB49muTkZOrUqcPatWtz/UZNQUnyEsJEGfK7jRUrViQkxLBblEvyEsJEKWG31NxI8hLCRBnyu43FQZKXECbKkN9tLA6SvIQwUQrPXZK8hDBVkryEEIokbw8SQiiSslOXJC8hTJZM2AshFEmWSgghFEnhHS9JXkKYKhk2CiEUSYaNQghFku82CiEUSeGjRkleQpgqWaQqhFAmZecuSV5CmCqZsBdCKJLSh43y9iAhhCJJz0sIEyU9LyGEIhnqvY0AWVlZhIWF0bp1a7y8vBg0aBA3b97UT+D/JclLCFOlyuMogoiICDZv3szMmTPZsmULKpWKIUOGkJ6eXsSg/0eSlxAmSq1S5XoUVnp6OmvWrGHUqFG0bduW2rVrs2jRIuLj49m/f7/+4tdbS0IIRVGpcj8K6+LFi6SmptKsWTNtma2tLXXq1OHEiRN6iPwFmbAXwkTl9d1GPz+/XOsPHDiQY/ndu3cBqFKlik555cqV+eOPPwoQYe5MJnmNbVPT2CGYnLGnw40dgshFWXPDtPv06VMALCwsdMotLS1JSkrS231MJnkJIQrmVT2rvFhZWQEv5r5e/hkgLS2NsmXL6iU2kDkvIYSevRwuJiQk6JQnJCTg7Oyst/tI8hJC6FXt2rUpX748MTEx2rLk5GTOnz9P48aN9XYfGTYKIfTKwsKCfv36MX/+fCpWrEjVqlWZN28ezs7OdOjQQW/3keQlhNC7wMBAnj9/zieffMKzZ8/w9vYmMjIy2yR+Uag0Go1Gb60JIUQxkTkvIYQiSfISQiiSJC8hhCJJ8hJCKJIkLyGEIknyEkIokiQvIYQiSfJSgIiICPr372/sMEq94ti6WOiPJK8S7vPPPycsLMzYYZiE4ti6WOiPJK8SKj4+ng8++IDQ0FBcXV2NHU6pV1xbFwv9keRVQp07dw47Ozt27NiBl5eXscMp9Ypr62KhP/LF7BLK19cXX19fY4dhMopr62KhP9LzEoLcty5OS0szRkgiD5K8hEB36+I/0/fWxUJ/JHkJQfFtXSz0R5KXEBTf1sVCf2TCXgiKb+tioT+SvIT4r+LYuljoj2wDLYRQJJnzEkIokiQvIYQiSfISQiiSJC8hhCJJ8hJCKJIkLyGEIknyEiWOrN4R+SHJqxTq378/Hh4eOkfdunXx8fFh+vTpJCUlGeS+UVFReHh4EBcXB8CSJUvw8PDI9/V3795l6NCh/P7770WOJS4uDg8PD6KioorcliiZZIV9KVWnTh2mTp2q/TkjI4Nz586xcOFCLly4wKZNm1CpVAaNoXfv3rRu3Trf5//4448cOnSITz/91IBRidJCklcpVb58eRo0aKBT5u3tTWpqKmFhYfz888/Z6vXN2dlZdmQQBiPDRhNTt25dAO7cuUP//v0ZN24cgYGBNGrUiA8//BB4sYfVZ599Rtu2balbty5du3Zl9+7dOu1kZWURERGBj48PXl5eDB8+PNtwNKdh465du+jRowdeXl74+Pgwb9480tPTiYqKYtKkSQD4+fkxceJE7TVbt26lc+fO2qHvkiVLeP78uU67+/bto1u3btSvXx9/f38uXryonw9MlFjS8zIx169fB+C1114DYM+ePXTq1ImlS5eSmZmJRqNhxIgRnDp1isDAQNzc3Ni/fz9jxowhPT2d7t27AzBv3jy++OILhg0bRoMGDfjuu+9YsGBBrvfevHkzU6dOpVevXowZM4a4uDg+++wzHj16xLhx4/joo49YtmwZ4eHh2qS3YsUKFi1aRL9+/Zg0aRIXLlxgyZIl/PHHH4SEhADwww8/EBgYSOfOnRk3bhwXL15k/PjxBvoERYmhEaVOv379NAEBAZqMjAztcf/+fc3u3bs1TZo00bzzzjuarKwsTb9+/TR169bVpKamaq89evSoxt3dXbNr1y6dNseNG6dp2bKlJiMjQ5OUlKTx9PTUzJkzR+ecwYMHa9zd3TW3b9/WaDQaTVhYmMbd3V2j0Wg0mZmZmhYtWmhGjBihc83atWs13bp106SlpWm2b9+uc31ycrLGy8tLM2XKFJ1rvvrqK427u7vm0qVLGo1Go+nRo4emR48eOuesWLFC4+7urtm+fXthP0ZRwsmwsZQ6ceIEnp6e2qNFixaMHTsWT09PFi5cqJ2sd3FxwdraWntddHQ0KpWKtm3b8vz5c+3h6+vLvXv3uHz5MmfOnCEjIwM/Pz+de7799tuvjOf69evcv3+f9u3b65QPGDCAb775JsdtZ06fPs3Tp0/x9fXNFgvAsWPHePbsGefOnStQLKJ0kGFjKeXp6cn06dMBUKlUWFpaUqVKFcqXL69zXqVKlXR+TkxMRKPR0KhRoxzbTUhIIDk5GYCKFSvq1Dk6Or4ynsTERAAcHBzy/Xd4ec3LubicYklKSkKj0WSLpXLlyvm+j1AmSV6lVLly5ahXr16Br7OxscHa2povvvgix/rq1avzyy+/APDgwQNq1qyprXuZbHJia2sLwMOHD3XKExMTOXfuXI5PPl9eM3/+fGrUqJGtvlKlSlSoUAG1Ws39+/eztStKNxk2Ch1NmjThyZMnaDQa6tWrpz0uX77M0qVLef78OQ0bNsTKyorvvvtO59qDBw++st2aNWtib2/PgQMHdMp37tzJkCFDSEtLQ63W/d/Ry8sLc3Nz4uPjdWIxNzdnwYIFxMXFYWlpScOGDdm3b5/OyvwffvhBD5+GKMmk5yV0tG3bFm9vb4YPH87w4cNxc3Pjl19+YcmSJbRq1Uo7PBs+fDiLFy+mbNmyNGvWjMOHD+eavMzMzBg1ahQzZsxg2rRpdOjQgRs3brB48WL++c9/UrFiRW1Pa//+/bRp0wY3Nzc++OADQkNDSUlJoWnTpsTHxxMaGopKpaJ27doAjB07lvfff5+RI0fSp08fbty4wbJlywz/YQmjkuQldKjValauXEloaCgrVqzgwYMHODk5MWDAAEaMGKE9b+jQoVhbW7Nu3TrWrVtHw4YNCQoKYtq0aa9sOyAgAGtrayIjI9m2bRtOTk4MGjRIO6fVtGlTWrRowYIFC4iOjmblypWMHj0aR0dHNm7cyOrVq7Gzs6N58+aMHTsWGxsbABo3bsyqVatYuHAhI0eOxMXFhZCQEIYNG2bQz0oYl+xhL4RQJJnzEkIokiQvIYQiSfISQiiSJC8hhCJJ8hJCKJIkLyGEIknyEkIokiQvIYQiSfISQiiSJC8hhCJJ8hJCKJIkLyGEIv0/Q2UGCaQyscAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 300x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Logistic regression model ---\n",
      "Accuracy: 0.800\n",
      "Precision: 0.615\n",
      "Recall: 1.000\n",
      "F1_score: 0.762\n",
      "FPR (Type I error): 0.294\n",
      "FNR (Type II error): 0.000\n"
     ]
    }
   ],
   "source": [
    "# Set probability threshold for postitive rip event\n",
    "yhat_class = yhat>0.5\n",
    "\n",
    "acc = accuracy_score(y_test, yhat_class)\n",
    "tn, fp, fn, tp = confusion_matrix(y_test, yhat_class).ravel()\n",
    "precision = (tp/(tp+fp))\n",
    "recall = (tp/(tp+fn))\n",
    "f1_score = 2*((precision*recall)/(precision+recall))\n",
    "\n",
    "import seaborn as sn\n",
    "df_cm = pd.DataFrame([[int(tp), int(fn)],[int(fp), int(tn)]], [1,0], [1,0])\n",
    "plt.figure(figsize=(3,2))\n",
    "sn.set(font_scale=1.0) # for label size\n",
    "sn.heatmap(df_cm, annot=True, fmt='d', annot_kws={\"size\": 12}, cmap=\"Blues\")\n",
    "plt.ylabel(\"Actual\")\n",
    "plt.xlabel(\"Predicted\")\n",
    "plt.title(\"Confusion matrix\")\n",
    "plt.show()\n",
    "matplotlib.rc_file_defaults()\n",
    "\n",
    "# Summarize test performance\n",
    "print('--- Logistic regression model ---')\n",
    "print('Accuracy: %.3f' % acc)\n",
    "print('Precision: %.3f' % precision )\n",
    "print('Recall: %.3f' % recall )\n",
    "print('F1_score: %.3f' % f1_score )\n",
    "print('FPR (Type I error): %.3f' % (fp/(fp+tn)) )\n",
    "print('FNR (Type II error): %.3f' % (fn/(tp+fn)) )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d808a74",
   "metadata": {},
   "source": [
    "Since a probability threshold of 0.5 appeared to be too low, we can try to increase it to 0.85, in an attempt to reduce the number of false positives. Recomputing the binary prediction and confusion matrix, we see that the number of false positives are now indeed significantly reduced, leaving the confusion matrix more diagonally dominant. The precision increases to 0.867, and the F1 score increases to 0.839. The recall is reduced, but still at an acceptable level of 0.812. We can conclude that our classification model with the predictors *sepal width* and *petal width* and a probability threshold of 0.85 provides good performance and a good compromise between recall and precision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2b99522b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAADyCAYAAAAcCUOdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAoo0lEQVR4nO3deVxUZf//8degAqKyKYJKCmFgogImijuC3mkuiVpaiLnvmjsquZZo7ohpKriLS4qlqaWZVhqBotyl0rdUXLhJEMQFNxDm90c/5m5uEEEYhsN8nj3O4xHXOXPOh3nEu+u65sx1VGq1Wo0QQiiMkb4LEEKIlyHhJYRQJAkvIYQiSXgJIRRJwksIoUgSXkIIRZLwEkIokoSXEEKRJLyEwZP7tJVJwqsU/fbbb0ydOhVvb2+aNGmCr68vH330ETdv3tTZNQ8fPkyHDh1o3Lgxs2fPLrHzuri4EBoaWmLn05fjx48TGBj4wuOmT5+Oj49PKVQkCkslXw8qHTt27CA4OJgWLVrg5+dHzZo1uXHjBmFhYaSnp7Np0yZcXV1L/LotWrTAwcGBSZMmYWtri4ODQ4mcNy4uDjs7O+zs7ErkfPoSEBAAwLZt2wo87saNG2RkZNCwYcPSKEsUgoRXKYiNjSUgIAB/f3+CgoK09t25c4devXphbm7OgQMHSvzaLi4ujBkzhvHjx5f4ucuDwoaXKHtk2FgKwsPDqVatGpMmTcqzz9ramunTp/Ovf/2LjIwMTfvhw4fp1asXHh4etG7dmtmzZ3Pv3j3N/tDQUDp16sTJkyfp3r07jRo14s0332T//v0AREdH4+LiAsBnn32Gi4sLiYmJ+Q5/EhMTcXFxITIyUtO2bds2OnfuTOPGjWnbti1z587Vqu9/h40pKSnMmDGD9u3b06RJE/r06cPx48e1ruPi4sKOHTsICgqiefPmeHh4MH78eFJTU5/73uXW9u233zJ69Gjc3d1p1aoVa9asISMjg5kzZ/LGG2/QqlUrlixZojV/lZiYyLRp02jTpg2urq60bNmSadOmkZ6eDvwdXDExMcTExODi4kJ0dLTmfdu1axcdOnSgVatWnDp1Sut9O378eJ7fPyEhAXd390INQUXJkPDSMbVazalTp2jZsiWVK1fO95jOnTszduxYqlatCsCaNWuYOHEibm5urFq1ijFjxvDtt98SEBDAkydPNK+7ffs28+fPZ8CAAaxfvx57e3umT5/OlStXcHV1Zffu3QD06dOH3bt3U7NmzULVfOjQIT799FP8/f0JDw9nzJgxfPXVV3zyySf5Hp+amkqfPn2IiYlh4sSJhIaGUqdOHcaMGZOnN7lixQpycnJYvnw506ZN4+TJkwQHB7+wpqCgIJydnVm7di1eXl6EhITQp08fTE1NCQkJwcfHh7CwML755hsAHj9+zIABA7hy5Qpz5swhPDyc/v378/XXX7N8+XIA5syZQ8OGDWnYsCG7d+/WGravWLGCwMBAAgMDcXd316rF19eXnj17sm7dOq5cuUJ2djbTp0/H2tqaWbNmFeo9FsVXUd8FlHfp6ek8ffoUe3v7Qh1/79491q5dyzvvvMOcOXM07c7Ozvj7+xMZGcn7778P/P0HumDBAlq2bAmAg4MDHTp04IcffmDw4MGaPzo7O7s8f4AFiY6Opk6dOvj7+2NkZETz5s0xMzPT9Fj+16ZNm7hz5w5HjhzhlVdeAaB9+/YMHDiQxYsX061bN4yMjDS/x8KFCzWv/fXXXzWBU5C2bdsyYcIEAOrXr8+hQ4eoXr265kOI1q1bc+TIEc6dO0eXLl24du0adnZ2LFq0iLp16wLg5eXFb7/9RkxMjOY8uf/D+N/3p1+/fnTu3Pm59QQFBREVFcXcuXNp27Ytv/32G1u3btWcT+ie9Lx0LPePNjs7u1DHx8XFkZmZSffu3bXamzVrRp06dYiOjtZq/+cfXe7k+aNHj4pR8d9/5NeuXaNXr16sWbOGS5cu0b17dz744IN8j4+JicHDw0MTXLl69OjB7du3uXr1ar715tb8+PHjF9bk4eGh+XcbGxsA3NzcNG0qlQoLCwsePHgAwOuvv05ERAT29vbcvHmTn376iY0bN3L16lWysrJeeL3cIffzmJub88knnxATE8OKFSsYNmwYzZo1e+F5RcmR8NIxS0tLqlSpQlJS0nOPefToEXfv3gXQzGvVqFEjz3E1atTQ/HHm+udQNDcoi/sZzFtvvcWyZcswMzNj9erV+Pn54evry6FDh/I9/t69e8+tF+D+/fv51ptbc2Hqza9H87xheK5NmzbRqlUrOnbsyPTp0/nll19e+Jpc1atXf+ExLVu2pFatWuTk5MhtFHog4VUK2rRpQ3R0NE+fPs13f2RkJC1btuT8+fNYWFgA5DuJffv2baysrIpVi0qlytMLzK+n1q1bNyIiIoiOjmblypVYWloydepUkpOT8xxrYWHx3HqBYtf8Mg4ePMiiRYsYPHgwUVFRnD59mvXr15fYrSLw9wchqampODk58dFHH5GZmVli5xYvJuFVCgYPHszdu3dZsWJFnn1paWmEhYVRr1493N3dcXNzw9jYmIMHD2odd/bsWZKSkmjatGmxaqlSpYpmHi7XuXPntI6ZMGECY8eOBaBatWp06dKF0aNHk52dTUpKSp5zenp6cv78+Tw32x44cAAbGxvq1atXrJpfRmxsLNWqVWP48OFYW1sD8PDhQ2JjY8nJydEcl9tbLarffvuNDRs2MHLkSJYtW8bVq1fLxU27SiIT9qXA3d2dDz/8kJUrV3LlyhX8/PywsrLizz//ZOPGjTx8+JD169ejUqmwtLRk+PDhrF69mkqVKuHr60tiYiIhISHUr1+fXr16FauWDh06sG3bNmbOnMk777yjqaFChQqaY7y8vJgzZw6ffvop7dq14/79+6xevRoHBwcaNGiQ55yDBg3iwIEDDBo0iLFjx2JlZcWXX37JL7/8QnBw8EsHRHE0adKEnTt3smjRIjp06EBKSgrh4eGkpqZqerfw99zV+fPniYqKKvQNqJmZmUyfPh1HR0eGDx+OsbExAwYMIDw8nI4dO2rNxQndkfAqJaNGjaJhw4bs2LGDhQsXcvfuXezs7GjXrh0jR46kdu3ammPHjRtHjRo12L59O1988QWWlpZ07tyZCRMmFHrO5nlat25NYGAg27Zt4+jRo7i6urJ69Wr69eunOaZfv35kZWWxa9cuIiIiMDU1pWXLlkydOpVKlSrlOaeNjQ07d+5k2bJlLFiwgKysLBo0aMCaNWvw9fUtVr0vy8/Pj8TERPbt20dERAS2tra0b9+e999/n1mzZnH58mXq16+Pv78/Fy5cYNiwYSxcuLBQt5OsXLmSq1evsnPnToyNjQEYP348R48eJTAwkK+++goTExNd/4oGT+6wF0Ioksx5CSEUScJLCKFIEl5CCEWS8BJCKJKElxBCkSS8hBCKJOElhFAkg7lJdetZ3a0TL/LXs1EdfZdgcMxNC98fqewxtsD9j8+vLm45OmUw4SWE+B9GFV58TBkm4SWEoVIpe9ZIwksIQyU9LyGEIik8vJTdbxRCvDyVquCtGNLS0pg6dSpeXl54eHgwfPhwLl++rNkfHx9P//79cXd3x9vbm/Dw8CJfQ8JLCENlVKHgrRhGjRrFzZs32bBhA3v37sXU1JSBAwfy+PFj0tPTGTRoEA4ODuzbt49x48YREhLCvn37inQNGTYKYah0NGxMT0/H3t6eUaNG8dprrwEwevRo3n77bf7880+ioqIwNjZm7ty5VKxYEScnJ65fv86GDRvo3bt34cvXSfVCiLJPZVTw9pKsrKxYvny5JrhSU1MJDw/Hzs6O+vXrc/bsWTw9PalY8b99Jy8vLxISEkhLSyv0daTnJYShqqD7CftZs2axZ88ejI2NWbt2LWZmZty6dQtnZ2et43JXsE1KSirUk5tAwksIw/WCYeOLlvA+fvz4Cy/xwQcf0LdvX3bu3MmYMWOIiIjgyZMnmuWzc+Uum/28J2zlR8JLCEP1oqFhCSwQX79+fQA+/vhj4uLi2L59O6ampnkeE5cbWmZmZoU+t4SXEIbqBT2v40df3LPKT1paGlFRUXTp0kXzVCojIyOcnJxISUnBzs4uzyP0cn+2tbUt9HVkwl4IQ6WjWyVSUlKYPHkyMTExmrasrCwuXbqEk5MTnp6exMbGaj38OCoqCkdHx0LPd4GElxCGS0c3qTZo0IA2bdowb948zp49yx9//EFgYCD3799n4MCB9O7dm4yMDIKCgrh8+TKRkZFs2bKFESNGFOk6MmwUwlDp6D4vlUrFypUrWbZsGRMmTODBgwc0a9aMHTt2aJ5PGhYWxoIFC/Dz88PGxoZp06bh5+dXtOsYynMbZT2v0ifreZW+Iq3n1X1NgfsfHxxd3HJ0SnpeQhiqYn5/Ud8kvIQwVApfVULCSwhDJYsRCiGUSGUk4SWEUCCVkcx5CSEUSCUT9kIIJTKSYaMQQolk2CiEUCQZNgohFEmGjUIIRZJhoxBCkWTYKIRQJBk2CiEUSYaNokTdT0th/fRhvDNxHvUaumva/+/saU7t307aXzepXNWcJu3epE3P96lQsZL+ii1HsrOz2bZ5I1/t38vtlGTq1nOg/weDeatbD32XpjMybBQl5t7tZHZ+Op2njx5qtV+Oi2bvyrm4tXsTn/eGkpp0k5O7w8m4m0bXoZP0VG35smbVCiK2b2XkmHG87tqI0z/9yJygQIyMjOj8Vjd9l6cTMmwUxabOyeHXn45yPGJdvvt/PrCT2q+60G34FAAcG73B4wf3OP1VBJ36j8LYtHJpllvuPHr0kN27dvB+/wF8MHgYAM1btOT3+Ivs2bm93IaXDBtFsSXfuMqRTSG80bEHjo2asntJkNb+7iOmkfOPhxUAVKhYiZycnDztouiMjU3YuHUn1WvU0GqvVLESDzMy9FSV7smwURSbRY2ajF62FfPqNly/FJdnv5Vtbc2/P3mUwbUL5/jl0Bc0auWLaZWqpVhp+VSxYkWcXRoAoFarSUtL5eBX+4mJjiJo9nw9V6c70vMSxVa5qjmVC5FB9++kEjquHwCWNna06z1Ax5UZnm8Of83smdMAaN22HZ06d9FzRbqj9DkvZVdvYIxNTPGfuYTeE+ZSuZo5Gz8aze3E6/ouq1xp1LgJ6zZuZebsefwef4khA94v0iPolUSlUhW4lXUSXgpiWqUqDq4eNPBsw3vTP0WNmpgje/VdVrnySt16NH3DE7/e7/Jx8BIu//kH3393VN9l6YTKSFXgVtbpddgYEBBQ6ITfunWrjqspm3Kys/k95kesa9lj5/Capr1ylWpY1azN/Tu39Vhd+XAnLY2fT/9Eq9Ztsf7HE5sbNmoMQPKtW/oqTadk2FgMLVu25MyZM6SlpVGnTp0CN0NlVKECx3du4PudYVrt91KTSU26gW1dJz1VVn48evyIebNm8NV+7V5s1OmfAHjNxUUfZemcjh6YXWr02vMaPXo0ZmZmrFq1inXr1mFvb6/Pcsqsdr0H8PX6pRzasIyGLb15kJ7Gqf3bqVzVnBZv9dF3eYpnb/8KXbu/Tdi6NRgZGdHQtTHxly6wccPneLVqQ6vWbfVdok4YKWBoWBC9f9o4cOBATp06xcqVK1m6dKm+yymT3Np3xti0Mj8f3M3FqBNUMjbByc0T775DqWJhpe/yyoWZs+dTt54DB76MZP3a1dSoYUPf9wMYMnyUIiavX4bSw0ulVqvV+i4iOTmZS5cu0aFDB51dY+vZmzo7t8hfz0aGO9zXF3PTws8ENZxZ8AcRl4L/VdxydErvPS8AW1tbbG1t9V2GEAZF6T2vMhFeQojSJ+ElhFAkpc/lSXgJYaCk5yWEUCQJLyGEIil81CjhJYShkp6XEEKRJLyEEIoknzYKIRRJel5CCEWS8BJCKJJBDBvPnDlTpJN6enq+VDFCiNKjy57X3bt3Wb58OSdPniQjIwMXFxcmT55Ms2bNAIiPj2fBggVcuHABS0tLAgICGDJkSJGuUajwKuyKp2q1GpVKRXx8fJGKEEKUPl2G16RJk0hLS2P58uVYW1sTERHBkCFDiIyMxNramkGDBtGxY0fmzZtHXFwc8+bNw9LSkt69exf6GoUKL0NdglmI8kxXo8br169z+vRpdu7cSdOmTQEICgrixx9/5Ouvv8bU1BRjY2Pmzp1LxYoVcXJy4vr162zYsKHkw6t58+Yv91sIIcqsCjrqeVlZWbF+/XoaNWqkaVOpVKjVau7du8eFCxfw9PSkYsX/xo+Xlxfr1q0jLS2N6v94jkBBXmrCPi4ujpiYGLKysshdy1CtVvPo0SNiY2PZs2fPy5xWCFGKdBVe5ubmtG/fXqvtyJEj3LhxgzZt2rBixQqcnZ219tesWROApKQk3YXXjh07+OSTT8hvAVYjIyPatGlT1FMKIfTgRfPYvr6+Be4/fvx4oa4TGxvLzJkz8fX1xcfHh4ULF2JsbKx1jImJCUCRnpFZ5KcHbd++nTZt2hAdHc2QIUN49913iYuLIyQkBBMTE3r06FHUUwoh9MBIpSpwKwnfffcdQ4YMoUmTJixfvhwAU1NTMjMztY7LDS0zM7NCn7vIPa/ExESmT5+OhYUFjRs3JjQ0FFNTU958800SEhLYunUr3bp1K+pphRCl7EXDxsL2rJ5n+/btLFiwgE6dOrF06VJNb8vOzo6UlBStY3N/Lspy8EXueVWqVAlTU1MAHBwcuH79OllZWQA0bdqUa9euFfWUQgg90OVzGyMiIvj444/x9/dn5cqVWsNET09PYmNjyc7O1rRFRUXh6OhY6PkueInwev311zlx4gQA9erVIycnh7i4OABuldMnCwtRHlUwUhW4vayEhASCg4Pp1KkTI0aMIC0tjdu3b3P79m0ePHhA7969ycjIICgoiMuXLxMZGcmWLVsYMWJEka5T5GHjoEGDGDt2LPfu3WPhwoX4+voybdo03nzzTQ4ePMgbb7xR1FMKIfRAVzepfvvtt2RlZXHs2DGOHTumtc/Pz49FixYRFhbGggUL8PPzw8bGhmnTpuHn51ek67zUcxtPnjzJlStXGDJkCOnp6UyePJlz587RuHFjFi9eTK1atYp6Sp2T5zaWPnluY+krynMb+245X+D+3R94FLccnXqp+7y8vb3x9vYG/r4hbePGjSVZkxCiFJTUJ4r6UuTwSkpKeuExtWvXfqlihBClR1c3qZaWIoeXj4/PC29uky9mC1H2KbzjVfTwCg4OzhNeuV8L+uWXXwgODi6x4oQQumNwPa9evXrl2+7v78+nn37KwYMHNfNhQoiyS+lzXkW+z6sg3t7enDx5siRPKYTQkdL4epAulegy0HFxcVrLXAghyi6DW8N+xowZedpycnL466+/OHv2LH369CmRwoQQumVwc17R0dF52lQqFVWrVmXYsGGMHDmyRAoTQuiWAkaGBSpyeH3//fe6qEPn3nV/Rd8lGBwrz7H6LsHgPD6/utDHVlB4ehV5wn7GjBncvJn/V22uXr0qPS8hFEJXX8wuLYXqef3zrvr9+/fTsWNHKlSokOe4H3/8kZ9//rnkqhNC6IwC8qlAhQqv+fPn88MPPwB/z2+NHZv/cECtVtO6deuSq04IoTNK6F0VpFDhNW/ePH7++WfUajUzZ85k1KhR1K1bV+sYIyMjzM3NadGihU4KFUKULIMIL1tbW81aOyqVCm9vb6pWrUqlSpUAePz4MU+fPsXS0lJnhQohSlaJ3qGuB0Wuv2vXrixbtox3331X03b+/HnatGnDggULtJZ2FUKUXUqfsC9yeK1atYrDhw/Ts2dPTZurqyuBgYHs37+fDRs2lGR9QggdqWBU8FbWFfk+r0OHDhEYGEjfvn01bRYWFgQEBGBkZMTmzZvldgkhFEAJ318sSJHDKz09HXt7+3z3OTo6kpycXOyihBC6p4TeVUGKXL6TkxPffvttvvuOHTtGvXr1il2UEEL3KqhUBW5lXZF7XoMHD2by5MncvXuXjh07Ur16de7cucN3333H0aNHWbhwoS7qFEKUMAXMyReoyOHVtWtXHjx4wOrVqzl69Kim3crKitmzZ8vTsoVQCCV8oliQl1p8q1+/fvTt25eEhATu3r2Lubk5JiYmfPHFF3h7e3Pq1KmSrlMIUcIMMrzg75tVHR0dOXHiBJ9//jmnT58mOzubV199tSTrE0LoiMKz6+XCKyUlhS+++IK9e/dy69YtzM3N6du3Lz179qRJkyYlXaMQQgeUMClfkCKF1+nTp9m1axcnTpxArVbTokULbt26xerVq/H09NRVjUIIHVD6rRKFCq+wsDD27NnDjRs3cHR0ZPz48fj5+WFiYkLz5s11XaMQQgcM4ibVpUuX4uLiwrZt27R6WA8ePNBZYUII3VL6hH2hOo49evTgxo0bDB06lBEjRnDkyBEyMzN1XZsQQocqqAreyrpC9bwWL17Mw4cP+frrr4mMjGTixIlYWFjg6+uLSqXK8wRtIUTZp/S/W5VarVYX9UVXrlxh7969HDx4kNTUVOrUqUO3bt3o2rUrzs7Ouqiz2J4803cFhkcewFH6ivIAju2xiQXu7/9G/t9hLiteKrxyZWdnc+LECfbu3cupU6fIzs7mtdde48CBAyVZY4mQ8Cp9El6lryjhteMF4eVfxsOrWI+3rlChAh07dqRjx46kpaURGRnJl19+WUKlCSF0SelPzC6xOz2qV6/OsGHDOHToUEmdUgihQwa3qoQQonwo+/FUMAkvIQyUEnpXBZHwEsJAGcQd9kKI8kfh2SXhJYShkmGjEEKRlD5sVPiiGOWTWq1m757d9PHrjlczD95605fFCxeQkZGh79LKDSMjFVMGdeLCV3O4E7Wc6N3T6feW9rJObd94jaNhH5L0w2ISjgWza+lQXn2lhp4qLnlGqArcSsqaNWsICAjQaouPj6d///64u7vj7e1NeHh4kc8r4VUGbd4YRvAn82jbzpuVoZ8xcPBQDh86yKQPx1KML0SIf5g/tgezRnVl0/7T9Prwc05E/x+bFnxA387NAGjRxJFDa8eSdvchg4I2M+nTPTjY1+D4xklUt6yi5+pLhpFRwVtJ2Lx5M6tWrdJqS09PZ9CgQTg4OLBv3z7GjRtHSEgI+/btK9K5ZdhYxuTk5BAetp4+7/Tlw4mTAfBq2QpLS0umTprApYsXcG3UWM9VKluVysaM6tee0O0nWLb5OwBOxvyBx+t1GfVee3Z/c5Ypg//F7wm3eH9quOZ/GFFxV/nzyMcEdPdi5bbj+vwVSoQu57ySk5MJCgoiNjYWR0dHrX179uzB2NiYuXPnUrFiRZycnLh+/TobNmygd+/ehb6G9LzKmIyMDLp260GXrtpPYapX7+//AG7evKmPssqVJ5nP8B64jFXbv9dqz8x6hnGlv/9/fvbCNVZHnNDq6d5Kvc/9h09wLCdDR9UL/imOixcvYmFhwYEDB3Bzc9Pad/bsWTw9PalY8b99Jy8vLxISEkhLSyv0NaTnVcaYm5szI2hWnvbj3/39mLn6r71W2iWVO9nZOfz2x380P9tWr0bA2y3xaeHC6I93AvBpWN4HK7dr9hrWFlW4dPmvUqtVl140Ye/r61vg/uPHn9/79PHxwcfHJ999t27dyrP6TM2aNQFISkqievXqBV43l17D69mzZxw9epSzZ8+SlJREZmYmlStXxs7OjmbNmtGpUyetdDZUcefPsSl8Ax18O1K/voRXSerXpRmbggcCcOSnC+z9Njbf42pYVWXNrPdJvJXO9oO/lGKFuqOvWyWePHmCsbGxVpuJiQkAT58+LfR59JYMN27cYNiwYSQnJ9OwYUNq1qyJhYUFT58+JT4+nn379hEaGkpYWBi1a9fWV5l6dy72LOPHjMT+lbrMm79A3+WUOzEXrtFxyAqc69kya1RXTmyeTNuAJTzN/O8aSrVsLDjw2RhsrKvy1ohQHj4uH6sIvyi7CupZFYepqWmelZhzQ8vMzKzQ59FbeM2bNw97e3v27t1LtWrV8uy/f/8+EydOZP78+Xz++ed6qFD/jhw+xOyg6Tg4OLJ2fTgWlpb6LqncuXozlas3Uzl97gpXE1P5Zv14/Hzd2XXkLACu9Wuzf9VIqpiZ8PbYNcReuqHnikuOvnpednZ2pKSkaLXl/mxra1vo8+htwj42NpZp06blG1zw99zP1KlTOXPmTClXVjZs3hjGjGmTaeLmzsatO6hhY6PvksoNG6uq+HdvgY1VVa322IvXAbC3swKgvacz32+aiEqlotOQlfzy74RSr1WXVKqCN13x9PQkNjaW7OxsTVtUVBSOjo6Fnu8CPYaXubl5nvT9X0lJSZiampZSRWXHF3t2sWLZEjq92ZnP14c/N+DFy6liZkLY/AAG+rXSav9Xq4YA/PrHf3BzsWdfyAhu3kqn/YClXLpSPibp/0mXnzYWpHfv3mRkZBAUFMTly5eJjIxky5YtjBgxokjn0duwsU+fPsyYMYPx48fTokULatWqhbGxMZmZmSQnJxMTE8PSpUvp06ePvkrUi9Tbt1n66UJq167De+/3Jz7+ktZ++1fqYm1trafqyodr/0lj+8FoZg7vQk5ODmcv3uCNhnUJHPomR09f4ujpS/wcEUilihVYsO4w9nZWmt4YwO30DBISU/X4G5QMfS2kWr16dcLCwliwYAF+fn7Y2Ngwbdo0/Pz8inSeYq1hXxxqtZrPPvuMTZs28ejRozz7q1Spgr+/Px9++CFGJXC7r1LWsN8fuZe5s4Keu3/+Jwt5269XKVb08sryGvbGlSoyYYAv/t2aU7eWNbdS77Hz8BkWbfiW2jUtiP963nNfu+3ALwyfs70Uqy28oqxhfzbhfoH7mzmaF7ccndJbeOXKysoiPj6e5ORkHj9+jKmpKXZ2djRo0CDPx6nFoZTwKk/KcniVV0UJr9hrBYfXGw5lO7z0fhNVpUqVaNKkib7LEMLgKHxRCf2HlxBCP5S+JI6ElxAGStnRJeElhMFSSc9LCKFECn/mrISXEAZLwksIoUQyYS+EUCSFZ5eElxCGSpffXywNEl5CGCiZsBdCKJLcKiGEUCSFZ5eElxCGSsJLCKFIcquEEEKRlB1dEl5CGCyZsBdCKJLcKiGEUCSFd7wkvIQwVDJsFEIokgwbhRCKJN9tFEIoksJHjRJeQhgquUlVCKFMys4uCS8hDJVM2AshFEnpw0YjfRcghBAvQ3peQhgopfe8JLyEMFAy5yWEUCYJLyGEEsmwUQihSArPLgkvIQyV0r/bqFKr1Wp9FyGEEEUl93kJIRRJwksIoUgSXkIIRZLwEkIokoSXEEKRJLyEEIok4SWEUCQJLyGEIkl4CSEUScJLCKFIEl5CCEWS8BJCKJKElxBCkSS8FGDNmjUEBATou4xyLycnh1WrVtG2bVvc3NwYPHgw169f13dZ4jkkvMq4zZs3s2rVKn2XYRDWrFnDrl27+OSTT9i9ezcqlYphw4aRmZmp79JEPiS8yqjk5GSGDh1KSEgIjo6O+i6n3MvMzGTjxo2MGzeO9u3b06BBA1asWEFycjLHjh3Td3kiHxJeZdTFixexsLDgwIEDuLm56buccu/333/n4cOHeHl5adrMzc1p2LAhZ86c0WNl4nlkGegyysfHBx8fH32XYTBu3boFQK1atbTaa9asyV9//aWPksQLSM9LCODx48cAGBsba7WbmJjw9OlTfZQkXkDCSwjA1NQUIM/k/NOnT6lcubI+ShIvIOElBP8dLqakpGi1p6SkYGdnp4+SxAtIeAkBNGjQgKpVqxIdHa1pu3//PpcuXaJZs2Z6rEw8j0zYC8Hfc139+/dn6dKlWFtbU6dOHZYsWYKdnR2dOnXSd3kiHxJeQvx/48eP59mzZ3z00Uc8efIET09PwsPD80zii7JBHjorhFAkmfMSQiiShJcQQpEkvIQQiiThJYRQJAkvIYQiSXgJIRRJwkuUOXL3jigMCa9yKCAgABcXF62tUaNGeHt7M2/ePO7du6eT60ZGRuLi4kJiYiIAoaGhuLi4FPr1t27dYsSIEfznP/8pdi2JiYm4uLgQGRlZ7HOJsknusC+nGjZsyJw5czQ/Z2VlcfHiRZYvX058fDw7d+5EpVLptIZ33nmHtm3bFvr4n3/+mZMnTzJr1iwdViXKCwmvcqpq1aq4u7trtXl6evLw4UNWrVrFv//97zz7S5qdnZ2syCB0RoaNBqZRo0YAJCUlERAQwJQpUxg/fjxNmzZl+PDhwN9rWC1evJj27dvTqFEjunfvzuHDh7XOk5OTw5o1a/D29sbNzY3Ro0fnGY7mN2w8dOgQvXr1ws3NDW9vb5YsWUJmZiaRkZHMmDEDAF9fX6ZPn655zRdffEHXrl01Q9/Q0FCePXumdd6jR4/So0cPmjRpgp+fH7///nvJvGGizJKel4FJSEgA4JVXXgHgyJEjdO7cmc8++4zs7GzUajVjxozh3LlzjB8/HicnJ44dO8bEiRPJzMykZ8+eACxZsoStW7cycuRI3N3d+eabb1i2bFmB1961axdz5syhT58+TJw4kcTERBYvXkx6ejpTpkxh1KhRrF27ltWrV2tCb926daxYsYL+/fszY8YM4uPjCQ0N5a+//iI4OBiA77//nvHjx9O1a1emTJnC77//ztSpU3X0DooyQy3Knf79+6v9/f3VWVlZmi01NVV9+PBhdfPmzdXvvvuuOicnR92/f391o0aN1A8fPtS89tSpU2pnZ2f1oUOHtM45ZcoUdevWrdVZWVnqe/fuqV1dXdWLFi3SOmbIkCFqZ2dn9c2bN9VqtVq9atUqtbOzs1qtVquzs7PVrVq1Uo8ZM0brNZs2bVL36NFD/fTpU/W+ffu0Xn///n21m5ubevbs2Vqv2bNnj9rZ2Vn9xx9/qNVqtbpXr17qXr16aR2zbt06tbOzs3rfvn0v+zaKMk6GjeXUmTNncHV11WytWrVi0qRJuLq6snz5cs1kvb29PWZmZprXRUVFoVKpaN++Pc+ePdNsPj4+3L59mz///JO4uDiysrLw9fXVumaXLl2eW09CQgKpqal07NhRq33gwIF89dVX+S47c/78eR4/foyPj0+eWgBOnz7NkydPuHjxYpFqEeWDDBvLKVdXV+bNmweASqXCxMSEWrVqUbVqVa3jatSoofXz3bt3UavVNG3aNN/zpqSkcP/+fQCsra219tnY2Dy3nrt37wJQvXr1Qv8Oua/JnYvLr5Z79+6hVqvz1FKzZs1CX0cok4RXOVWlShUaN25c5NdVq1YNMzMztm7dmu/+evXq8euvvwKQlpbGq6++qtmXGzb5MTc3B+DOnTta7Xfv3uXixYv5fvKZ+5qlS5fi4OCQZ3+NGjWwtLTEyMiI1NTUPOcV5ZsMG4WW5s2b8+jRI9RqNY0bN9Zsf/75J5999hnPnj3Dw8MDU1NTvvnmG63Xnjhx4rnnffXVV7GysuL48eNa7QcPHmTYsGE8ffoUIyPt/xzd3NyoVKkSycnJWrVUqlSJZcuWkZiYiImJCR4eHhw9elTrzvzvv/++BN4NUZZJz0toad++PZ6enowePZrRo0fj5OTEr7/+SmhoKG3atNEMz0aPHs3KlSupXLkyXl5e/PDDDwWGV4UKFRg3bhzz589n7ty5dOrUiWvXrrFy5Uree+89rK2tNT2tY8eO0a5dO5ycnBg6dCghISFkZGTQokULkpOTCQkJQaVS0aBBAwAmTZrEBx98wNixY+nbty/Xrl1j7dq1un+zhF5JeAktRkZGrF+/npCQENatW0daWhq2trYMHDiQMWPGaI4bMWIEZmZmbNmyhS1btuDh4UFgYCBz58597rn9/f0xMzMjPDycvXv3Ymtry+DBgzVzWi1atKBVq1YsW7aMqKgo1q9fz4QJE7CxsSEiIoKwsDAsLCxo2bIlkyZNolq1agA0a9aMDRs2sHz5csaOHYu9vT3BwcGMHDlSp++V0C9Zw14IoUgy5yWEUCQJLyGEIkl4CSEUScJLCKFIEl5CCEWS8BJCKJKElxBCkSS8hBCKJOElhFAkCS8hhCJJeAkhFEnCSwihSP8PhUPNHaJdBMQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 300x200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--- Logistic regression model ---\n",
      "Accuracy: 0.900\n",
      "Precision: 0.867\n",
      "Recall: 0.812\n",
      "F1_score: 0.839\n",
      "FPR (Type I error): 0.059\n",
      "FNR (Type II error): 0.188\n"
     ]
    }
   ],
   "source": [
    "# Set probability threshold for postitive rip event\n",
    "yhat_class = yhat>0.85\n",
    "\n",
    "acc = accuracy_score(y_test, yhat_class)\n",
    "tn, fp, fn, tp = confusion_matrix(y_test, yhat_class).ravel()\n",
    "precision = (tp/(tp+fp))\n",
    "recall = (tp/(tp+fn))\n",
    "f1_score = 2*((precision*recall)/(precision+recall))\n",
    "\n",
    "import seaborn as sn\n",
    "df_cm = pd.DataFrame([[int(tp), int(fn)],[int(fp), int(tn)]], [1,0], [1,0])\n",
    "plt.figure(figsize=(3,2))\n",
    "sn.set(font_scale=1.0) # for label size\n",
    "sn.heatmap(df_cm, annot=True, fmt='d', annot_kws={\"size\": 12}, cmap=\"Blues\")\n",
    "plt.ylabel(\"Actual\")\n",
    "plt.xlabel(\"Predicted\")\n",
    "plt.title(\"Confusion matrix\")\n",
    "plt.show()\n",
    "matplotlib.rc_file_defaults()\n",
    "\n",
    "# Summarize test performance\n",
    "print('--- Logistic regression model ---')\n",
    "print('Accuracy: %.3f' % acc)\n",
    "print('Precision: %.3f' % precision )\n",
    "print('Recall: %.3f' % recall )\n",
    "print('F1_score: %.3f' % f1_score )\n",
    "print('FPR (Type I error): %.3f' % (fp/(fp+tn)) )\n",
    "print('FNR (Type II error): %.3f' % (fn/(tp+fn)) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e11b915",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "work",
   "language": "python",
   "name": "work"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}