{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Logistic Regression with the Iris Dataset\n",
"\n",
"For more explanation of logistic regression, see\n",
"1. [Our course notes](https://jennselby.github.io/MachineLearningCourseNotes/#binomial-logistic-regression)\n",
"1. [This scikit-learn explanation](http://scikit-learn.org/stable/modules/linear_model.html#logistic-regression)\n",
"1. [The full scikit-learn documentation of the LogisticRegression model class](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html)\n",
"\n",
"## Instructions\n",
"0. If you haven't already, follow [the setup instructions here](https://jennselby.github.io/MachineLearningCourseNotes/#setting-up-python3) to get all necessary software installed.\n",
"0. Read through the code in the following sections\n",
" * [Iris Dataset](#Iris-Dataset)\n",
" * [Visualization](#Visualization)\n",
" * [Model Training](#Model-Training)\n",
" * [Prediction](#Prediction)\n",
"0. Complete at least one of the following exercises\n",
" * [Exercise Option #1 - Standard Difficulty](#Exercise-Option-#1---Standard-Difficulty)\n",
" * [Exercise Option #2 - Advanced Difficulty](#Exercise-Option-#2---Advanced-Difficulty)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from sklearn import linear_model # for fitting our model\n",
"from sklearn.datasets import load_iris # the iris dataset is included in scikit-learn\n",
"\n",
"# force numpy not to use scientific notation, to make it easier to read the numbers the program prints out\n",
"import numpy\n",
"numpy.set_printoptions(suppress=True)\n",
"\n",
"# to display graphs in this notebook\n",
"%matplotlib inline\n",
"import matplotlib.pyplot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Iris Dataset\n",
"\n",
"Before you go on, make sure you understand this dataset. Modify the cell below to examine different parts of the dataset that are contained in the 'iris' dictionary object.\n",
"\n",
"What are the features? What are we trying to classify?"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"dict_keys(['data', 'target', 'target_names', 'DESCR', 'feature_names'])"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"iris = load_iris()\n",
"iris.keys()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can also try looking at it using a [pandas dataframe](https://jennselby.github.io/MachineLearningCourseNotes/#pandas)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
"
],
"text/plain": [
" sepal length (cm) sepal width (cm) petal length (cm) \\\n",
"count 150.000000 150.000000 150.000000 \n",
"mean 5.843333 3.054000 3.758667 \n",
"std 0.828066 0.433594 1.764420 \n",
"min 4.300000 2.000000 1.000000 \n",
"25% 5.100000 2.800000 1.600000 \n",
"50% 5.800000 3.000000 4.350000 \n",
"75% 6.400000 3.300000 5.100000 \n",
"max 7.900000 4.400000 6.900000 \n",
"\n",
" petal width (cm) \n",
"count 150.000000 \n",
"mean 1.198667 \n",
"std 0.763161 \n",
"min 0.100000 \n",
"25% 0.300000 \n",
"50% 1.300000 \n",
"75% 1.800000 \n",
"max 2.500000 "
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"iris_df.describe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For this tutorial, at least to start, we're not going to use the whole dataset, just because it is easier to visualize two features than four. The code below decides which two features we're going to use.\n",
"\n",
"We'll also need to know at what location in the list each of the classes start at."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Use just two columns (the first and fourth in this case).\n",
"x1_feature = 0\n",
"x2_feature = 3\n",
"iris_inputs = iris.data[:,[x1_feature,x2_feature]]\n",
"\n",
"# The data are in order by class. Find out where the other classes start in the list\n",
"start_class_one = list(iris.target).index(1)\n",
"start_class_two = list(iris.target).index(2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualization\n",
"\n",
"Let's visualize our dataset, so that we can better understand what it looks like."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEjCAYAAADdZh27AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xu8VHXZ///XBWzl5DFPCAJqinEUtgfUUhPM462ZkiGYeLgtTbu7NX+31X0nqd39uu3Ou7IizbOUClqZmnlIS5OSDXISQkkROQiIyhm2G67vH2vtzexh75k1zJo1a2bez8djPfaedbz2Yphr1lqf6/Mxd0dERASgQ7kDEBGR9FBSEBGRFkoKIiLSQklBRERaKCmIiEgLJQUREWmhpCBVwczGmNnT5Y5DpNIpKUhFMLOFZjayveXuPtHdP7MD+33BzDaZ2VozW2Nm08zsejPbuYB9uJl9vNBjZ2z/STN72cxWm9n7ZvZXMzsyiWOLZFNSkIpnZp2K3MVV7r4L0AO4FvgC8KSZWdHB5WFmuwKPAz8B9gR6At8BNpf62CJtUVKQimNm48Jv07ea2SpgfDjvpXC5hctWhN/+Z5vZwHz7dff17v4CcBZwDHBGuL+jzGyKmX1oZsvM7DYz2ylc9pdw85lmts7MzjezPczscTNbaWYfhL/3auewh4bH/rW7b3H3je7+tLvPyvh7LzGzeeG+/mhmfdo7dqHnUiSbkoJUqqOBN4F9ge9mLfsMcDzBB+5uwOeBVVF37O6LgAbgU+GsLcC/A3sRJIsRwJXhuseH6wxx9+7u/hDB/6u7gT5Ab2AjcFs7h3sd2GJm95rZaWa2R+ZCMzsb+CbwOWBv4EXg1zmOLVIUJQWpVEvd/Sfu3uTuG7OWfQTsAhwGmLvPc/dlhe6f4HYO7j7N3f8WHmsh8AvghPY2dPdV7v6Iu29w97UESavN9d19DfBJwIE7gJVm9piZ7Ruu8mXge+Hf0AT8N3B489WCSNyUFKRSvdPeAnf/E8E3858CK8zs9vDefSF6Au8DmNmh4S2gd81sDcEH817tbWhmXc3sF2b2drj+X4DdzaxjO/HOc/dx7t4LGAjsD/xfuLgP8KPw1tWHYUwWxicSOyUFqVQ5u/d19x+7ez3Qn+A20nVRd2xmBwD1BLdqAH4O/AM4xN13Jbidk+sh9LVAP+DocP3m2zx5H1y7+z+AewiSAwTJ70vuvnvG1MXdX47694gUQklBqo6ZHWlmR5tZHbAe2ARsjbBdVzM7Afgd8ArwZLhoF2ANsM7MDgOuyNp0OXBQxutdCJ4jfGhmewI35DjmYWZ2bfOD6DAhjQb+Fq4yAfiGmQ0Il+9mZqNyHFukKEoKUo12Jbg//wHwNsFD5ltyrH+bma0l+ID9P+AR4FR3b04kXwcuANaG+81+oDseuDe8xfP5cB9dgPcIPtyfynHstQQPzf9uZuvD9ecQXG3g7r8Bvg88GN6KmgOcluPYIkUxDbIjIiLNdKUgIiItlBRERKSFkoKIiLRQUhARkRZKCiIi0kJJQUREWigpiIhICyUFERFpoaQgIiItlBRERKSFkoKIiLRQUhARkRZKCiIi0kJJQUREWigpiIhICyUFERFpoaQgIiItOpU7gELttdde3rdv33KHISJSUaZNm/aeu++db72KSwp9+/aloaGh3GGIiFQUM3s7ynq6fSQiIi2UFEREpIWSgoiItKi4Zwpt+eijj1i8eDGbNm0qdyhVoXPnzvTq1Yu6urpyhyIiCauKpLB48WJ22WUX+vbti5mVO5yK5u6sWrWKxYsXc+CBB5Y7HBFJWFXcPtq0aRMf+9jHlBBiYGZ87GMf01WXSI2qiqQAKCHESOdSpHaVLCmY2QFm9ryZzTWz18zs39pY50QzW21mM8Lp26WKRyT9VgMDwp+1HEMUlRJn5SnllUITcK279weGA18xs/5trPeiux8eTjeWMJ7UuOeee1i6dGm5w5DUeQKYCzxZ4zFEUSlxVp6SJQV3X+bu08Pf1wLzgJ6lOl4lUVKQ1i4AugMXha+/GL6+oMZiiKJS4qxciTxTMLO+wFDg720sPsbMZprZH8xsQDvbX25mDWbWsHLlyuIDmjgR+vaFDh2CnxMnFr3L9evXc8YZZzBkyBAGDhzIQw89xLRp0zjhhBOor6/nlFNOYdmyZUyePJmGhgbGjBnD4YcfzsaNG3nuuecYOnQogwYN4pJLLmHz5s0AXH/99fTv35/Bgwfz9a9/HYDf//73HH300QwdOpSRI0eyfPnyomOXcrsR6A00NwGuA/oAN9VYDFFUSpwVzN1LOhGk8WnA59pYtivQPfz9dOCNfPurr6/3bHPnzt1uXrseeMC9a1d32DZ17RrML8LkyZP9sssua3n94Ycf+jHHHOMrVqxwd/cHH3zQL774Ynd3P+GEE3zq1Knu7r5x40bv1auXz58/393dL7zwQr/11lv9vffe80MPPdS3bt3q7u4ffPCBu7u///77LfPuuOMOv+aaa4qKuz0FnVOJwSR37+Tu3cKfk2o0higqJc50ARo8wmd2Sa8UzKwOeASY6O6PtpGQ1rj7uvD3J4E6M9urlDHxrW/Bhg2t523YEMwvwqBBg3jmmWf4j//4D1588UXeeecd5syZw8knn8zhhx/OzTffzOLFi7fbbv78+Rx44IEceuihAFx00UX85S9/YbfddqNz585ceumlPProo3Tt2hUIajJOOeUUBg0axC233MJrr71WVNySFg8D3YDvhD8n1WgMUVRKnJWplK2PDLgTmOfuP2xnnf3C9TCzo8J4VpUqJgAWLSpsfkSHHnoo06dPZ9CgQfznf/4njzzyCAMGDGDGjBnMmDGD2bNn8/TTT0feX6dOnXjllVc477zzePzxxzn11FMBuPrqq7nqqquYPXs2v/jFL1RPUDWuA+YD14Y/r6vRGKKolDgrUykrmo8DLgRmm9mMcN43CW4I4u4TgPOAK8ysCdgIfCG8zCmd3r3h7TZ6kO3du6jdLl26lD333JOxY8ey++6787Of/YyVK1cyZcoUjjnmGD766CNef/11BgwYwC677MLatWsB6NevHwsXLmTBggV8/OMf5/777+eEE05g3bp1bNiwgdNPP53jjjuOgw46CIDVq1fTs2fwvP7ee+8tKmZJkyMzft83nGoxhigqJc7KVLKk4O4vATmroNz9NuC2UsXQpu9+Fy6/vPUtpK5dg/lFmD17Ntdddx0dOnSgrq6On//853Tq1ImvfvWrrF69mqamJr72ta8xYMAAxo0bx5e//GW6dOnClClTuPvuuxk1ahRNTU0ceeSRfPnLX+b999/n7LPPZtOmTbg7P/xhcLE1fvx4Ro0axR577MFJJ53EW2+9VVTctWE1cCzwMrBbmWPJJQ1xpiEGKasoDx7SNBX9oNk9eKjcp4+7WfCzyIfM1ai6HjRP9ODt86tyB5JHGuJMQwxSCqThQXNqjRkDCxfC1q3BzzFjyh2RlESltGlPQ5xpiEHSoDaTgtSISmnTnoY40xCDpIGSglSxjxN82H1E0HTxI4JmjAeXM6g2pCHONMQgaaCkIFWuUtq0pyHONMQg5VYVg+yItO864CcEzRbHAu+UN5x2pSHONMQg5aakIFWuUtq0pyHONMQg5abbRyn17W9/m2effbbg7V544QXOPPPMEkQkIrVASaGM3J2tW7e2uezGG29k5MiRJY+hqamp5McQiGdQmCQGlolyjHzrLAJ2Dn/Wuso7FzWbFJYtg4MPhnffLX5f119/PT/96U9bXo8fP54f/OAH3HLLLRx55JEMHjyYG264AYCFCxfSr18/vvjFLzJw4EDeeecdxo0bx8CBAxk0aBC33norAOPGjWPy5MkATJ06lWOPPZYhQ4Zw1FFHsXbtWjZt2sTFF1/MoEGDGDp0KM8///x2cb3//vt89rOfZfDgwQwfPpxZs2a1xHfhhRdy3HHHceGFFxZ/AiSCOAaFSWJgmSjHyLfO94FG4JZ4Q6tIFXguolS4pWmKpaLZ3a+4wr1DB/crryx40+1Mnz7djz/++JbXn/jEJ/yee+7xf/3Xf/WtW7f6li1b/IwzzvA///nP/tZbb7mZ+ZQpU9zdvaGhwUeOHNmybXMX2RdddJFPmjTJN2/e7AceeKC/8sor7u6+evVq/+ijj/wHP/hBS1fc8+bN8wMOOMA3btzozz//vJ9xxhnu7n7VVVf5+PHj3d39ueee8yFDhri7+w033ODDhg3zDRs2tPs3VVdFczmN9m1dPOPbunwenfA+4jhGvnX6eNv/bfvEGGel6ONpOxeoorl9y5bB3XcHBc1331381cLQoUNZsWIFS5cuZebMmeyxxx4tvaIOHTqUYcOG8Y9//IM33ngDgD59+jB8+HAADjroIN58802uvvpqnnrqKXbddddW+54/fz49evTgyCODh4C77rornTp14qWXXmLs2LEAHHbYYfTp04fXX3+91bYvvfRSy5XASSedxKpVq1izZg0AZ511Fl26dCnuD5cI4igKS6KwLMox8q1zJ7BT1n53Au6KMc5KUbnnoiaTwk03BQkBYMuW4HWxRo0axeTJk3nooYc4//zzcXe+8Y1vtHSdvWDBAi699FIAunXr1rLdHnvswcyZMznxxBOZMGECl112WfHBRJAZg5RSHEVhSRSWRTlGvnVGAFdl7fcq4KQY46wUlXsuai4pNF8lNDYGrxsb47laOP/883nwwQeZPHkyo0aN4pRTTuGuu+5i3bp1ACxZsoQVK1Zst917773H1q1bOffcc7n55puZPn16q+X9+vVj2bJlTJ06FYC1a9fS1NTEpz71KSaGw4i+/vrrLFq0iH79+rXaNnOdF154gb322mu7KxFJQhxFYUkUlkU5Rr51Hg5/npn1uhZV5rmouTqFzKuEZs1XCxnPigs2YMAA1q5dS8+ePenRowc9evRg3rx5HHPMMQB0796dBx54gI4dO7babsmSJVx88cUtrZC+973vtVq+00478dBDD3H11VezceNGunTpwrPPPsuVV17JFVdcwaBBg+jUqRP33HMPO++8c6ttx48fzyWXXMLgwYPp2rWrxl8omziKwpIoLItyjHzr3AzUAwOBOcB0aldlngsLnj9UjiOOOMIbGhpazZs3bx6f+MQnIm3fqxcsWbL9/J49oY3RMmtWIedURNLPzKa5+xH51qu520eLF4P79pMSgqRfvvqAJOoYpDCV929Sc0lBpHLlqw9Ioo5BClN5/yZKCiKpl28AHA2Qkz6V+2+ipCCSevnqAzRATvpU7r+JkoJI6uWrD9AAOelTuf8mSgoiFSFKfYAGyEmXyvw3UVIokaVLl3LeeecVvN1ll13G3Llzc64zYcIE7rvvvh0NTSrSdcB84Nrw53UFLpfkVea/Sc3VKZRbU1MTnTqlv2awks6piOSnOoW84ms/3F7X2QMHDgTgnnvu4ayzzuKkk05ixIgRbN26lSuvvJLDDjuMk08+mdNPP72lm+wTTzyR5qTXvXt3vvWtbzFkyBCGDx/O8uXLW+0fYMGCBYwcOZIhQ4YwbNgw/vnPf7Ju3TpGjBjBsGHDGDRoEL/73e+K/hslDSplPIVKGTsiDpUSZ3Q1nBTiaz98/vnn8/DD2/o1efjhhzn66KNbrTN9+nQmT57Mn//8Zx599FEWLlzI3Llzuf/++5kyZUqb+12/fj3Dhw9n5syZHH/88dxxxx3brTNmzBi+8pWvMHPmTF5++WV69OhB586d+c1vfsP06dN5/vnnufbaa6m0K0JpS6WMp1ApY0fEoVLijK4Gk0L87Yfb6jr7gAMOaLXOySefzJ577gkEXVqPGjWKDh06sN9++/HpT3+6zf3utNNOLUNr1tfXs3DhwlbL165dy5IlSzjnnHMA6Ny5M127dsXd+eY3v8ngwYMZOXIkS5YsabnKkEqURJv3KMdIol6iUtr3V0qchavBpFCa9sPZXWdn25Guquvq6jAzADp27Bh56MyJEyeycuVKpk2bxowZM9h3333ZtGlTwceXtKiU8RQqZeyIOFRKnIWrwaRQmvbD2V1n53LcccfxyCOPsHXrVpYvX84LL7ywQ8fcZZdd6NWrF7/97W8B2Lx5Mxs2bGD16tXss88+1NXV8fzzz/P222/v0P4lLSplPIVKGTsiDpUSZ+FqMClAKdoPZ3edncu5555Lr1696N+/P2PHjmXYsGHstttuO3Tc+++/nx//+McMHjyYY489lnfffZcxY8bQ0NDAoEGDuO+++zjssMN2aN+SJpU0nkIljB0Rh0qJszA12iR1KsGl377AcoI+4fO21IrVunXr6N69O6tWreKoo47ir3/9K/vtt1+iMeSiJqlpk8R7Nsox8q0TR5zl//8ZTaXEGYjaJDX9DeZL4siM3/cNp2SdeeaZfPjhhzQ2NvJf//VfqUoIkkZJvGejHCPfOnHEWf7/n9FUSpyFKVlSMLMDgPsIzpQDt7v7j7LWMeBHwOnABmCcu1fG8ERF2tHnCLVnNXAs8DLQ1i22fMtrySLgEOANgm+wUn5xvD+TfY+X8plCE3Ctu/cHhgNfMbP+WeucRvAuPgS4HPj5jh6s0m6DpVm6zqXGEIju+0AjcEu5A5EWlVezUbKk4O7Lmr/1u/taYB7QM2u1s4H7PPA3YHczy/2Utg2dO3dm1apVKfswq0zuzqpVq+jcuXOZI9EYAtH1BQz4Wfj6tvB13zLFI5Vcs5HIMwUz6wsMBf6etagnrUf+XhzOW5a1/eUEVxL07r39ZXGvXr1YvHgxK1eujC3mWta5c2d69epV5ihuBGYACwkuOttqE59reS25k+AObGPGvJ2Au8oTjhDP+7M87/GSJwUz6w48AnzN3dfsyD7c/XbgdghaH2Uvr6ur48ADDywqTkmb5nbgowma+22m7Tbx7S2vJSOAq4AfZsy7CjipPOEI8bw/y/MeL2mdgpnVESSEie7+aBurLAEy+4PoFc4TQWMIFKK5760zs15L+VRmzUbJ6hTClkX3Au+7+9faWecMgq80pwNHAz9296Ny7betOgWpVkm0ia8W9wL1wEBgDjCd4B60lE+6ajai1imUMil8EngRmA1sDWd/k7CtnLtPCBPHbcCpBE1SL3b3nJ/4SgoiIoUre/Gau79E0AQi1zoOfKVUMYgEKqX9vmoykqXz2ZYa7ftIakultN9XTUaydD7boqQgVawvldF+XzUZydL5zEVJQarYnQTt9TOlsf1+EuMUyDY6n7koKUgVa26/nymN7feTGKdAttH5zEVJQapcpbTfV01GsnQ+21OjXWdL7biZ7dvvp9F1wE8I2qOPpXXvL1GWS2F0PttTFYPsiIhIblHrFHT7SGreu++u5o03BrB8+eoca60GBoQ/S7EcgnqKncOfpTpGPnHsIwlJnIvapKQgNe/JJ5/gkEPm8sQTudqrF1tDEKVNfL56iiTqGCql7b5qOkrG3Stqqq+vd5F4jPYtW7p5Y2Mnd8cbGzv5li3d3H10q3Xcu7l7sE7wM3OdYpe7u/fxtt/ufWI8Rv5zUfw+kpDEuahOQINH+IzVlYLUsBtZvrw3jY1Be/XGxjqWL89ur15sDUGUNvH56imSqGOolLb7qukouSiZA9iD4AbdQUCHKNuUatKVgsRl6VL30aMneWNjJ1+zJrhiGD16ki9blr3mJN/2jbNT+DrO5e7u13jrt/o1JThGPnHsIwlJnIvqQ7FXCma2m5l908xmA38DfkHQuPdtM5tkZp8ufcoSKZ2bboLPfe5h1q/vxvjx32H9+m6cc84kbtruS2WxNQRR2sTnq6dIoo6hUtruq6ajlNptkmpmzwD3Ab939w+zltUDFwKz3f3OkkeZQU1SJS69ekGPHlNZtKg3K1bsyz77LOeAA97h3XePYPHizDWLHdchSp/4+cZDSGJsiUoZn0LjbOyIso+nUCpKCiIihYu1TsHMBpvZWWb2ueap+BBlx6kNtpTGsmVw8MHw7rvljkTKJW9SMLO7CJpBnAv8SzidmXMjKTG1wZbSuOkmWLiQNp6rSK3Ie/vIzOa6e/+E4smrtm8fXQA8BmwGmgi6rtoZOAv4VRnjkmqwbBkcdBBs2gRdusCbb8J++5U7KolLnLePpphZapJCbVMbbCmdm26CreFo6lu26GqhVkVJCvcRJIb5ZjbLzGab2axSByZtUT/wUhrLlsHdd0NjY/C6sTF4rWcLtSdKUriToPnpqWx7nvAvpQxKclEbbIlf5lVCM10t1KYozxSmuPsxCcWTV20/UwC1wZZS6NULlizZfn7PnmTVbEilivpMIcogO6+a2a+A3xM84QTA3R8tIj7ZYUdm/L5vOIkURx/80izK7aMuBMngM6hJqtSoGTOgrg5mlfBpmmoECqFanVLJe6Xg7hcnEYhImo0dC01NcMEFMGdOaY6RWSPw05+W5hjVI7NWZ3SZY6kuUYrX7jWz3TNe7xEWtInUhBkz4LXXgt9fe600VwvNrX+2blWrn9wuALoDF4Wvvxi+vqBsEVWbKLePBmd2iOfuHwBDSxeSSLqMHdv69QUl+PxRjUBUqtUptShJoYOZ7dH8wsz2JNoDapGKl3mV0CzuqwXVCBRCtTqlFiUp/C9B8dpNZnYT8DLwP6UNSyQdsq8SmsV5taAagUKpVqeUInWdHXZzcVL48k/uPrekUeWgOgVJUseO239gA3ToEHxwx0E1AoVSrc6OKLpOwcy6u/s6gDAJbJcIMtcRqUZxffDnog/+QqlWp5Ry3T76nZn9r5kdb2bdmmea2UFmdqmZ/ZGg64s2mdldZrbCzNpswGdmJ5rZajObEU7f3vE/Q9Iojnb3SbTdj3KMd99dzRtvDGD58h1rFx/HMfLtIy3nW/UWla3dpODuI4DngC8Br5nZGjNbBTwA7Adc5O6Tc+z7HnIkjdCL7n54ON1YWOiSdnH0zZ9E//5RjvHkk09wyCFzeeKJHRvDIo5j5NtHWs63xmSocO5esgnoC8xpZ9mJwOOF7rO+vt4l/ZYude/c2R3cu3RxX7asPPso/hijfcuWbt7Y2Mnd8cbGTr5lSzd3H53oMfLtIy3nO4l/M9kxQINH+IyNNBxnCR1jZjPN7A9mNqC9lczscjNrMLOGlStXJhmf7KA42t0n0XY//zFuZPny3jQ2Bu3iGxvrWL68sHbxcRwj3z7Scr5Vb1EFomSOHZ3IfaWwK9A9/P104I0o+9SVQvplfltsngr91hjHPuI4xtKl7qNHT/LGxk6+Zk3wbX706EmR44jjGPn2kZbzncS/mew40n6l4O5rfFvrpieBOjPbq1zxSHziaHefRNv9KMe46Sb43OceZv36bowf/x3Wr+/GOedMihxHHMfIt4+0nG/VW1SJKJkD6AjsT9A4uDfQO+J2fWn/SmE/ttVJHAUsan6da9KVQvr17Nn622Lz1LNnsvuI4xg9e7ofccQrvs8+7zq477PPu15fPzVyHHEcI98+0nK+k/g3kx1HxCuFKIPsXA3cQFAl0vw9wN19cJ7tfk3wMHmvcNsbCDsscfcJZnYVcAXBCPQbgWvc/eV8SUzFayIihYtavBbl9tG/Af3cfYC7DwqnnAkBwIOmEz3cvc7de7n7ne4+wd0nhMtvC/c5xN2HR0kIkiy1N49u1qzVzJ07gDlz0t2/fxy1DsW+L/S+SrcoSeEdNJJFTVJ78+geeOAJ+vefy7337lgdQ1LiqHUo9n2h91W6tXv7yMyuCX8dAPQjGNUiczjOH5Y8ujbo9lEyli2Dgw6CTZugSxd4803Yb79yR5VGF7Bly2Ns3bqZuromPvqoEx067EzHjmcBvyp3cK3k+zeN8m9e7PtC76vyieP20S7htAh4BtgpY173OIKU9FJ786huZOHC1jUGb72Vzv7946h1KPZ9ofdVBcj3JBoYFWVeUpNaH5We2ptH9+qr7uee27rG4NxzJ/nMmeWOrLU4ah2KfV/ofVVexFin8I2I86RKqL15dGPHwuc/37rGYNSoSSUZna0YcdQ6FPu+0PuqMuR6pnAaQaXx54GHMhbtCvR396NKH9729Eyh9NS/f3QdO8KwYVNZtKg3K1bsyz77LOeAA97h1VePSKTb7ajy/ZtG+Tcv9n2h91V5FT2eArAUmAacFf5sthb49+LCkzTTf9Dogg/+9Pfvn+/fNMq/ebHvC72vKkO7ScHdZwIzzWyiu3+UYEwiIlIm7T5TMLPZZjYLmGZms7KnBGOUKjZjBtTVwax23lHPPANm8Kc/7djypCRRFKYBcCQJuR40nwn8C/BUOI0Jpz8A6a7QkYoxdiw0NdHug9nzzw9+nnfeji1PShJFYRoARxKRr3kS8Gob86ZHadpUiklNUqvHq6+2bp6Y3Yzz6adbL3/uucKWJyWOAXAqZRAdqVzE2CTVzOy4jBfHEq17DJGcxo5t/Tr7aqH5KqBZ9tVAvuVJSaIoTAPgSGLyZQ2gHpgJLATeBmYAw6JknFJMulKoDtlXCdlXC9lXAdlXA/mWJyWJojANgCNxIK4rBXef5u5DgCHAYHc/3N2nlyxLSU3Ivkpo1ny1kH0V0Kz5aiDf8qQkURSmAXAkSbmK18a6+wMZHeO14uoQT4rQseP2H1IAHToEH1Zm7W/rnn95UpIoCouj6EuFYxJH8Vq38Ocu8YQksk2+at98H+xJfvDnkkRRWBwf2vrgl6javX3k7r8If/2+u38ne0ooPqlySQzYkq8WIg5q/y/VIkorojlm9lcz+//N7Awz263kUUnNSGLAlny1EHFQ+3+pFnnHaAYws97Ap4DjCDrJ+9DdDy9xbG3SM4XqkcSALTNmwNCh217PnAmD8w4mWxgNHCOVILYxms2sF0Ey+BQwFHiN1r2miuyQJAZsyVcLEQe1/5dqkvdKwcy2AlOB/3b33yUSVQ66UqgOmd+umxXyLTvK9tlXCc3ivFoo9u8QSUpsVwoEVwf3AReY2RQzu8/MLi06QqlpSQzYkq8WIg5q/y/VJuozhe7AJwluIY0FcPc+pQ2tbbpSqA5JDNiSrxYiDmr/L5UijjqF5h01ADsDLwMvAse7+9vFhyi1LIkBW5IY+Uwf/FJtotw+Os3dB7n7l9z9ASWE6pCWdvVJjCEgItFF6ftoZRKBSLLS0q4+iTEERCRMiX8lAAASKElEQVS6SM8U0kTPFIqXlnb1+eJIS5wi1SDO1kdSZdLSrj6JMQREpDC5ekn9XK4N3f3RkkSUh64UipOWdvX54khLnCLVIo4rhX/JMZ0ZR5CSvLS0q09iDAERKZyeKdSYtLSrT2IMARHZJrY6hXBnZwADgM7N89z9xh0PT8olLR+oSYwhICKFi9Ih3gTgfOBqwIBRQN5qZjO7y8xWmNmcdpabmf3YzBaY2SwzG1Zg7DUrLWMIxCGJ8RTiqHVQvYTUiiitj4519y8CH4SD6xwDHBphu3uAU3MsPw04JJwuB34eYZ9CesYQiEMS4ynEUeugegmpGe6ecwL+Hv78G7A/QZcXC/JtF27TF5jTzrJfAKMzXs8HeuTbZ319vdeypUvdO3d2B/cuXdyXLdt+nVdfDZY3TzNnJh9nFFH+lmK3L/YYce1DpNyABo/wuR3lSuFxM9sduAWYDiwEfh1DPuoJvJPxenE4bztmdrmZNZhZw8qVtV1gnZYxBOKQxHgKcdQ6qF5Cakq+rAHsnPk7sFvmvDzb9qX9K4XHgU9mvH4OOCLfPmv5SiHzG2vzlP3NNfsqIa1XC1H+lmK3L/YYce1DJA2I8UphSkYC2ezuqzPnFWEJcEDG617hPGlHWsYQiEMS4ynEUeugegmpNe0mBTPbz8zqgS5mNtTMhoXTiUDXGI79GPDFsBXScGC1uy+LYb9V67HHoLGx9bzGRvhdxnh48+a1vW1788slyt9S7PbFHiOufYhUklx1CqcA4wi+wf8wY/4a4Jv5dmxmvwZOBPYys8XADUAdgLtPAJ4ETgcWABuAiwuOvsakZQyBOCQxnkIctQ6ql5Ba025ScPd7gXvN7Fx3f6TQHbv76DzLHfhKofsVEZHSifJM4a9mdqeZ/QHAzPprjGYRkeoUJSncDfyRoEYB4HXgayWLSEREyiZKUtjL3R8GtgK4exNQIXeuRUSkEFGSwnoz+xjgAM0thUoalYiIlEWUXlKvIWg+erCZ/RXYGzivpFGJiEhZ5E0K7j7dzE4A+hH0kjrf3T8qeWQiIpK4vEnBzDoDVwKfJLiF9KKZTXD3Tbm3FBGRShPl9tF9wFrgJ+HrC4D7CcZVEBGRKhIlKQx09/4Zr583s7mlCkhERMonSuuj6WGLIwDM7GhAgySLiFShKFcK9cDLZrYofN0bmG9mswl6qxhcsuhERCRRUZJCriE1RUSkikRpkvp2EoGIiEj5RXmmICIiNUJJQUREWigpiIhICyUFERFpoaQgIiItlBRERKSFkoKIiLRQUhARkRZKCiIi0kJJQUREWigpVKOJE6FvX+jQIfg5cWK5IxKRChGlQzypJBMnwuWXw4YNweu33w5eA4wZU764RKQi6Eqh2nzrW9sSQrMNG4L5IiJ5KClUm0WLCpsvIpJBSaHa9O5d2HwRkQxKCtXmu9+Frl1bz+vaNZgvIpKHkkK1GTMGbr8d+vQBs+Dn7bfrIbOIRKLWR9VozBglARHZISW9UjCzU81svpktMLPr21g+zsxWmtmMcLqslPFISHUMItKOkl0pmFlH4KfAycBiYKqZPebuc7NWfcjdrypVHJJFdQwikkMprxSOAha4+5vu3gg8CJxdwuNJFKpjEJEcSpkUegLvZLxeHM7Ldq6ZzTKzyWZ2QFs7MrPLzazBzBpWrlxZilhrh+oYRCSHcrc++j3Q190HA88A97a1krvf7u5HuPsRe++9d6IBVh3VMYhIDqVMCkuAzG/+vcJ5Ldx9lbtvDl/+EqgvYTwCqmMQkZxKmRSmAoeY2YFmthPwBeCxzBXMrEfGy7OAeSWMR0B1DCKSU8laH7l7k5ldBfwR6Ajc5e6vmdmNQIO7PwZ81czOApqA94FxpYpHMqiOQUTaUdJnCu7+pLsf6u4Hu/t3w3nfDhMC7v4Ndx/g7kPc/dPu/o9SxlMz8tUhXHkldOoUXCl06hS8jtvIkcH+m6eRI+M/huotROLn7hU11dfXu+TwwAPuXbu6w7apa9dgvrv7FVe0XtY8XXFFfDGMGNH2MUaMiO8Y+f5OEWmF4A5N3s9YC9atHEcccYQ3NDSUO4z06ts3KEjL1qcPLFwYXBls2bL98o4doakpnhjM2l8W1/st398pIq2Y2TR3PyLfeuVukipxy1eH0FZCyDU/rVRvIVISSgrVJl8dQseObS9vb35aqd5CpCSUFKpNvjqE5n6OsrU3f0eMGFHY/B2heguRklBSqDb56hB+9jO44optVwYdOwavf/az+GJ49tntE8CIEcH8uKjeQqQk9KBZRKQG6EGziIgUTEmhEPmKpeIopkpiH0kUryVBxWsi8YtSzJCmqWzFa/mKpeIopkpiH0kUryVBxWsiBUHFazHLVywVRzFVEvtIongtCSpeEylI1GcKSgpRdejQdjWuGWzdmn95HMeIYx9JVBsnIY5zJVJD9KA5bvmKpeIopkpiHypeE5EclBSiylcsFUcxVRL7SKJ4LQkqXhMpjSgPHtI0lbWX1AcecO/Tx90s+Jn9UDPf8jiOEcc+rrjCvWPH4OFsx46V95C5WRznSqRGoAfNIiLSTM8UdkSxA8NE2T5fjUC+fURpm9+zZ+t99OwZf5xx1EIkUS+hWgaRwkS5nEjTVLLbR8UODBNl+3w1Avn2EaVt/v77t72P/fePL844aiGSqJdQLYNIC3T7qEDFNtWMsn2+GoF8+4jSNj/fPuKIM45aiCTqJVTLINJCt4/SqNgBbpIaWCZfnHEM5JPEYD8aiEekYEoKSSq2RiCptvn54oyjFiKJegnVMogUTEmhWbEDw0TZPl+NQL59RGmbv//+be+jeX4cccZRC5FEvYRqGUQKF+XBQ5qmktYpZD+EjfqQuZDt89UI5NtHlLb52Q+bmx8yxxlnHLUQSdRLqJZBxN31oFlERDLoQXO2Smmvni/OYmsp4opDRKpTlMuJNE07dPuoUtqr54uz2FqKuOIQkYqDbh9lqJT26vniTKrb60o5XyISmW4fZaqU9uppiTMtcYhI4mojKVRKe/W0xJmWOEQkcbWRFCqlvXq+OIutpYgrDhGpWrWRFMaMgdtvD+6JmwU/b789mJ8m+eJ89tntE8CIEcH8JOMQkapVGw+aRURqXCoeNJvZqWY238wWmNn1bSzf2cweCpf/3cz6ljIeERHJrWRJwcw6Aj8FTgP6A6PNrH/WapcCH7j7x4Fbge+XKh4REcmvlFcKRwEL3P1Nd28EHgTOzlrnbODe8PfJwAizXI3xRUSklEqZFHoC72S8XhzOa3Mdd28CVgMfy96RmV1uZg1m1rBy5coShSsiIhXR+sjdb3f3I9z9iL333rvc4YiIVK1SJoUlwAEZr3uF89pcx8w6AbsBq0oYk4iI5NCphPueChxiZgcSfPh/Abgga53HgIuAKcB5wJ88TxvZadOmvWdmbXTMk6i9gPfKHEMUijNelRBnJcQIijNuUeLsE2VHJUsK7t5kZlcBfwQ6Ane5+2tmdiNBb32PAXcC95vZAuB9gsSRb79lv39kZg1R2vuWm+KMVyXEWQkxguKMW5xxlvJKAXd/Engya963M37fBIwqZQwiIhJdRTxoFhGRZCgp7Jjbyx1ARIozXpUQZyXECIozbrHFWXF9H4mISOnoSkFERFooKeRgZh3N7FUze7yNZePMbKWZzQiny8oRYxjLQjObHcaxXReyFvhx2PHgLDMblsIYTzSz1Rnn89tt7SeBOHc3s8lm9g8zm2dmx2QtL/u5jBhn2c+nmfXLOP4MM1tjZl/LWqfs5zNinGU/n2Ec/25mr5nZHDP7tZl1zlpedCejJW19VAX+DZgH7NrO8ofc/aoE48nl0+7eXjvl04BDwulo4Ofhz6TlihHgRXc/M7Fo2vYj4Cl3P8/MdgKyRhtKzbnMFyeU+Xy6+3zgcGjpIHMJ8Jus1cp+PiPGCWU+n2bWE/gq0N/dN5rZwwTN+O/JWK2lk1Ez+wJBJ6PnF3IcXSm0w8x6AWcAvyx3LDE4G7jPA38DdjezHuUOKm3MbDfgeIL6Gdy90d0/zFqt7OcyYpxpMwL4p7tnF56W/XxmaS/OtOgEdAl7gOgKLM1aXnQno0oK7fs/4P8DtuZY59zwkneymR2QY71Sc+BpM5tmZpe3sTxK54Slli9GgGPMbKaZ/cHMBiQZXOhAYCVwd3jb8Jdm1i1rnTScyyhxQvnPZ6YvAL9uY34azmem9uKEMp9Pd18C/ABYBCwDVrv701mrRepkNBclhTaY2ZnACneflmO13wN93X0w8AzbsnM5fNLdhxFcin/FzI4vYyztyRfjdKCPuw8BfgL8NukACb6FDQN+7u5DgfXAdoNDpUCUONNwPgEIb2+dBUwqVwxR5Imz7OfTzPYguBI4ENgf6GZmY+M+jpJC244DzjKzhQTjQJxkZg9kruDuq9x9c/jyl0B9siG2imVJ+HMFwb3Qo7JWidI5YUnli9Hd17j7uvD3J4E6M9sryRgJvqUudve/h68nE3z4Zir7uSRCnCk5n81OA6a7+/I2lqXhfDZrN86UnM+RwFvuvtLdPwIeBY7NWqfoTkaVFNrg7t9w917u3pfgcvJP7t4qI2fd9zyL4IF04sysm5nt0vw78BlgTtZqjwFfDFt6DCe47FyWphjNbL/me59mdhTBezPRHnPd/V3gHTPrF84aAczNWq2s5zJqnGk4nxlG0/4tmbKfzwztxpmS87kIGG5mXcNYRrD9505zJ6MQsZPRbGp9VABr3ZnfV83sLKCJoDO/cWUKa1/gN+H7tRPwK3d/ysy+DODuEwj6nzodWABsAC5OYYznAVeYWROwEfhCoW/mmFwNTAxvJbwJXJyycxk1zlScz/BLwMnAlzLmpe58Roiz7OfT3f9uZpMJbmU1Aa8Ct1uRnYxmU0WziIi00O0jERFpoaQgIiItlBRERKSFkoKIiLRQUhARkRZKClLTwt4v2+oFt835MRzvs2bWP+P1C2aWd2xdM+sRRzxmtreZPVXsfqR6KSmIJOuzQP+8a23vGuCOYg/u7iuBZWZ2XLH7kuqkpCCpFlZDPxF2RDbHzM4P59eb2Z/DDvb+2FxhHn7z/pEFfd7PCatPMbOjzGxK2IHcyxnVwFFjuMvMXgm3PzucP87MHjWzp8zsDTP7n4xtLjWz18Nt7jCz28zsWILq91vC+A4OVx8Vrve6mX2qnTDOBZ4K993RzH4Q/n2zzOzqcP5CM/teuO8GMxsWnpt/NhdihX4LjIn690ttUUWzpN2pwFJ3PwOCbqPNrI6gU7Kz3X1lmCi+C1wSbtPV3Q+3oNO9u4CBwD+AT7l7k5mNBP6b4IM2im8RdBdwiZntDrxiZs+Gyw4HhgKbgflm9hNgC/BfBP0RrQX+BMx095fN7DHgcXefHP49AJ3c/SgzOx24gaCPmxZmdiBBH/nNfW1dDvQFDg//nj0zVl8U/u23EvSzfxzQmaBbkQnhOg3AzRH/dqkxSgqSdrOB/zWz7xN8mL5oZgMJPuifCT9UOxJ0Jdzs1wDu/hcz2zX8IN8FuNfMDiHoxruugBg+Q9BB4tfD152B3uHvz7n7agAzmwv0AfYC/uzu74fzJwGH5tj/o+HPaQQf9tl6EHSV3WwkMCHsGpnm44QeC3/OBrq7+1pgrZltNrPdw3EXVhD0simyHSUFSTV3f92CIRpPB242s+cIell9zd2PaW+zNl7fBDzv7udYMEThCwWEYcC54Qhd22aaHU1whdBsCzv2f6p5H+1tv5EgERWyr61ZsW3N2HfncJ8i29EzBUk1M9sf2ODuDwC3ENySmQ/sbeG4xGZWZ60HPWl+7vBJgl43VxN0IdzcJfO4AsP4I3B12DMlZjY0z/pTgRPMbA8Lui/OvE21luCqpRCv0/oK4hngS+G+ybp9FMWhbN+TrgigpCDpN4jgHv4MgvvtN7t7I0Gvld83s5nADFr3K7/JzF4luId+aTjvf4DvhfML/TZ/E8Htpllm9lr4ul3h2BH/DbwC/BVYSDACFgTjc1wXPrA+uO09bLe/9cA/zezj4axfEnSjPCv8+y8o7M/h08ATBW4jNUK9pEpVMbMXgK+7e0OZ4+ju7uvCb/O/Ae5y97YGg4+6v3OAenf/zxhi+wvBQ/oPit2XVB9dKYiUxvjw6mYO8BZFDt8YJpSFxQZlZnsDP1RCkPboSkFERFroSkFERFooKYiISAslBRERaaGkICIiLZQURESkhZKCiIi0+H+Og+b9HW7gbAAAAABJRU5ErkJggg==\n",
"text/plain": [
"