{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# CS109A Introduction to Data Science \n", "\n", "\n", "## Lab 3: Scikit-learn for Regression\n", "\n", "**Harvard University**
\n", "**Fall 2018**
\n", "**Instructors:** Pavlos Protopapas and Kevin Rader
\n", "**Authors:** David Sondak, Will Claybaugh, Eleni Kaxiras\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# RUN THIS CELL TO PROPERLY HIGHLIGHT THE EXERCISES\n", "import requests\n", "from IPython.core.display import HTML\n", "styles = requests.get(\"https://raw.githubusercontent.com/Harvard-IACS/2018-CS109A/master/content/styles/cs109.css\").text\n", "HTML(styles)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Table of Contents\n", "
    \n", "
  1. Array creation and reshape
  2. \n", "
  3. Some plotting\n", "
  4. Simple linear regression
  5. \n", "
  6. $k$-nearest neighbors
  7. \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Learning Goals\n", "\n", "Overall description and goal for the lab.\n", "\n", "By the end of this lab, you should be able to:\n", "* Understand array reshaping\n", "* Review how to make plots\n", "* Feel comfortable with simple linear regression\n", "* Feel comfortable with $k$ nearest neighbors\n", "\n", "**This lab corresponds to lecture 4 and maps on to homework 2 (and beyond).**" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# import the necessary libraries\n", "import warnings\n", "warnings.filterwarnings('ignore')\n", "%matplotlib inline\n", "import numpy as np\n", "import scipy as sp\n", "import matplotlib as mpl\n", "import matplotlib.cm as cm\n", "import matplotlib.pyplot as plt\n", "import pandas as pd\n", "import time\n", "pd.set_option('display.width', 500)\n", "pd.set_option('display.max_columns', 100)\n", "pd.set_option('display.notebook_repr_html', True)\n", "import seaborn as sns" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Simple Linear Regression\n", "Linear regression and its many extensions are a workhorse of the statistics and data science community, both in application and as a reference point for other models. Most of the major concepts in machine learning can be and often are discussed in terms of various linear regression models. Thus, this section will introduce you to building and fitting linear regression models and some of the process behind it, so that you can 1) fit models to data you encounter 2) experiment with different kinds of linear regression and observe their effects 3) see some of the technology that makes regression models work.\n", "\n", "\n", "### Linear regression with a toy dataset\n", "We first examine a toy problem, focusing our efforts on fitting a linear model to a small dataset with three observations. Each observation consists of one predictor $x_i$ and one response $y_i$ for $i = 1, 2, 3$,\n", "\n", "\\begin{align*}\n", "(x , y) = \\{(x_1, y_1), (x_2, y_2), (x_3, y_3)\\}.\n", "\\end{align*}\n", "\n", "To be very concrete, let's set the values of the predictors and responses.\n", "\n", "\\begin{equation*}\n", "(x , y) = \\{(1, 2), (2, 2), (3, 4)\\}\n", "\\end{equation*}\n", "\n", "There is no line of the form $\\beta_0 + \\beta_1 x = y$ that passes through all three observations, since the data are not collinear. Thus our aim is to find the line that best fits these observations in the *least-squares sense*, as discussed in lecture." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise (10 min)
\n", "\n", "* Make two numpy arrays out of this data, x_train and y_train\n", "* Check the dimentions of these arrays\n", "* Try to reshape them into a different shape\n", "* Make points into a very simple scatterplot\n", "* Make a better scatterplot" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# your code here\n", "x_train = np.array([1,2,3])\n", "y_train = np.array([2,3,6])\n", "type(x_train)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3,)" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x_train.shape" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(3, 1)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x_train = x_train.reshape(3,1)\n", "x_train.shape" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(2, 3)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xx = np.array([[1,3,5],[6,2,1]])\n", "xx.shape" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[1, 3],\n", " [5, 6],\n", " [2, 1]])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xx = xx.reshape(3,-1)\n", "xx" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(3,) (3,)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAEtVJREFUeJzt3W+MXfV95/H3Z81UcQlbd9eTADZedyVkaUmTgEYOS6KKZNUCDi1sxQNW26KilSxQdpVKlavSB6yyT/oAqUooWiwr7SqoSaOoBS9CJl6iUIU0gmgMDuafV16WLgysmLBriMuowe53H9zjZJid8ZyZuXeu55f3S7qac37nN+d8ffjxmTPnnju/VBWSpLb8o3EXIEkaPsNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KALxnXgrVu31s6dO8d1eEnakI4cOfLDqppcrt/Ywn3nzp1MT0+P6/CStCEl+ds+/bwtI0kNMtwlqUGGuyQ1yHCXpAYZ7pLUoF7hnmRLkr9M8lKSF5P8ywXbk+TeJCeSPJvkqtGUK0nqo++jkF8CvllVtyT5OeDnF2y/Abi8e30CuL/7Kkk/8w4+M8M9h4/z+sk5Lt2ymX3X7eLmK7eN9JjLhnuSXwB+BfgdgKr6MfDjBd1uAh6owZx9T3ZX+pdU1RtDrleSNpSDz8xw14PHmHvvDAAzJ+e468FjACMN+D63ZX4JmAX+S5Jnknw5yYUL+mwDXp23/lrXJkk/0+45fPwnwX7W3HtnuOfw8ZEet0+4XwBcBdxfVVcCfwf8wWoOlmRvkukk07Ozs6vZhSRtKK+fnFtR+7D0CffXgNeq6qlu/S8ZhP18M8Bl89a3d23vU1UHqmqqqqYmJ5f90wiStOFdumXzitqHZdlwr6r/DbyaZFfX9K+AFxZ0exi4rXtq5mrgbe+3SxLsu24Xmyc2va9t88Qm9l23a4nvGI6+T8v8B+Cr3ZMyLwO3J7kDoKr2A4eAPcAJ4F3g9hHUKkkbztk3Tdf7aZkMHnBZf1NTU+VfhZSklUlypKqmluvnJ1QlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3qNRNTkleAHwFngNML/1B8kmuB/wr8z67pwar6T8MrU5K0En2n2QP4dFX98Bzbn6iqG9dakCRp7bwtI0kN6hvuBXwryZEke5foc02SZ5M8muSKxTok2ZtkOsn07OzsqgqWJC2v722ZT1XVTJIPAY8leamqvjNv+9PAjqo6lWQPcBC4fOFOquoAcAAGE2SvsXZJ0hJ6XblX1Uz39U3gIWD3gu3vVNWpbvkQMJFk65BrlST1tGy4J7kwyUVnl4FfA55b0OfiJOmWd3f7fWv45UqS+uhzW+bDwENddl8AfK2qvpnkDoCq2g/cAtyZ5DQwB9xaVd52kaQxWTbcq+pl4GOLtO+ft3wfcN9wS5MkrZaPQkpSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSgXuGe5JUkx5IcTTK9yPYkuTfJiW6S7KuGX6okqa++E2QDfLqqfrjEthsYTIh9OfAJ4P7uqyRpDIZ1W+Ym4IEaeBLYkuSSIe1bkrRCfcO9gG8lOZJk7yLbtwGvzlt/rWuTJI1B39syn6qqmSQfAh5L8lJVfWelB+t+MOwF2LFjx0q/XZLUU68r96qa6b6+CTwE7F7QZQa4bN769q5t4X4OVNVUVU1NTk6urmJJ0rKWDfckFya56Owy8GvAcwu6PQzc1j01czXwdlW9MfRqJUm99Lkt82HgoSRn+3+tqr6Z5A6AqtoPHAL2ACeAd4HbR1OuJKmPZcO9ql4GPrZI+/55ywV8brilSZJWy0+oSlKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa1Dvck2xK8kySRxbZdm2St5Mc7V53D7dMSdJK9Jlm76zPAy8C/3iJ7U9U1Y1rL0mStFa9rtyTbAc+C3x5tOVIkoah722ZLwK/D/zDOfpck+TZJI8muWLtpUmSVmvZcE9yI/BmVR05R7engR1V9VHgT4CDS+xrb5LpJNOzs7OrKliStLw+V+6fBH4jySvA14HPJPnz+R2q6p2qOtUtHwImkmxduKOqOlBVU1U1NTk5ufbqJUmLWjbcq+quqtpeVTuBW4FvV9Vvze+T5OIk6ZZ3d/t9awT1SpJ6WMnTMu+T5A6AqtoP3ALcmeQ0MAfcWlU1nBIlSSuVcWXw1NRUTU9Pj+XYkrRRJTlSVVPL9fMTqpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDWod7gn2ZTkmSSPLLItSe5NcqKbJPuq4ZYpSVqJlVy5fx54cYltNwCXd6+9wP1rrEuStAa9wj3JduCzwJeX6HIT8EANPAlsSXLJkGqUJK1Q3yv3LwK/D/zDEtu3Aa/OW3+ta5MkjcGy4Z7kRuDNqjqy1oMl2ZtkOsn07OzsWncnSVpCnyv3TwK/keQV4OvAZ5L8+YI+M8Bl89a3d23vU1UHqmqqqqYmJydXWbIkaTnLhntV3VVV26tqJ3Ar8O2q+q0F3R4GbuuemrkaeLuq3hh+uZKkPi5Y7TcmuQOgqvYDh4A9wAngXeD2oVQnSVqVFYV7Vf018Nfd8v557QV8bpiFSZJWz0+oSlKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIa1GeC7A8k+X6SHyR5PskXFulzbZK3kxztXnePplxJUh99ZmL6e+AzVXUqyQTw3SSPVtWTC/o9UVU3Dr9ESdJKLRvu3RR6p7rVie5VoyxKkrQ2ve65J9mU5CjwJvBYVT21SLdrkjyb5NEkVwy1SknSivQK96o6U1UfB7YDu5N8ZEGXp4EdVfVR4E+Ag4vtJ8neJNNJpmdnZ9dStyTpHFb0tExVnQQeB65f0P5OVZ3qlg8BE0m2LvL9B6pqqqqmJicn11C2JOlc+jwtM5lkS7e8GfhV4KUFfS5Okm55d7fft4ZfriSpjz5Py1wCfCXJJgah/Y2qeiTJHQBVtR+4BbgzyWlgDri1eyNWkjQGfZ6WeRa4cpH2/fOW7wPuG25pkqTV8hOqktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNajPTEwfSPL9JD9I8nySLyzSJ0nuTXKimyT7qtGUK0nqo89MTH8PfKaqTiWZAL6b5NGqenJenxuAy7vXJ4D7u6/ShnLwmRnuOXyc10/OcemWzey7bhc3X7lt3GVJK9ZnJqYCTnWrE91r4RR6NwEPdH2fTLIlySVV9cZQq5VG6OAzM9z14DHm3jsDwMzJOe568BiAAa8Np9c99ySbkhwF3gQeq6qnFnTZBrw6b/21rk3aMO45fPwnwX7W3HtnuOfw8TFVJK1er3CvqjNV9XFgO7A7yUdWc7Ake5NMJ5menZ1dzS6kkXn95NyK2qXz2Yqelqmqk8DjwPULNs0Al81b3961Lfz+A1U1VVVTk5OTK61VGqlLt2xeUbt0PuvztMxkki3d8mbgV4GXFnR7GLite2rmauBt77dro9l33S42T2x6X9vmiU3su27XmCqSVq/P0zKXAF9JsonBD4NvVNUjSe4AqKr9wCFgD3ACeBe4fUT1SiNz9k1Tn5ZRCzJ4wGX9TU1N1fT09FiOLUkbVZIjVTW1XD8/oSpJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJalCfafYuS/J4kheSPJ/k84v0uTbJ20mOdq+7R1OuJKmPPtPsnQZ+r6qeTnIRcCTJY1X1woJ+T1TVjcMvUZK0UsteuVfVG1X1dLf8I+BFwEklJek8tqJ77kl2AlcCTy2y+ZokzyZ5NMkVQ6hNkrRKfW7LAJDkg8BfAb9bVe8s2Pw0sKOqTiXZAxwELl9kH3uBvQA7duxYddGSpHPrdeWeZIJBsH+1qh5cuL2q3qmqU93yIWAiydZF+h2oqqmqmpqcnFxj6ZKkpfR5WibAnwIvVtUfL9Hn4q4fSXZ3+31rmIVKkvrrc1vmk8BvA8eSHO3a/hDYAVBV+4FbgDuTnAbmgFurqkZQrySph2XDvaq+C2SZPvcB9w2rKEnS2vgJVUlqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg5adrCPJZcADwIeBAg5U1ZcW9AnwJWAP8C7wO1X19PDLHTj4zAz3HD7O6yfnuHTLZvZdt4ubr9w2qsNJ0obTZ5q908DvVdXTSS4CjiR5rKpemNfnBuDy7vUJ4P7u69AdfGaGux48xtx7ZwCYOTnHXQ8eAzDgJamz7G2Zqnrj7FV4Vf0IeBFYmKI3AQ/UwJPAliSXDL1a4J7Dx38S7GfNvXeGew4fH8XhJGlDWtE99yQ7gSuBpxZs2ga8Om/9Nf7/HwAk2ZtkOsn07OzsyirtvH5ybkXtkvSzqHe4J/kg8FfA71bVO6s5WFUdqKqpqpqanJxczS64dMvmFbVL0s+iXuGeZIJBsH+1qh5cpMsMcNm89e1d29Dtu24Xmyc2va9t88Qm9l23axSHk6QNadlw756E+VPgxar64yW6PQzcloGrgber6o0h1vkTN1+5jT/6zV9m25bNBNi2ZTN/9Ju/7JupkjRPn6dlPgn8NnAsydGu7Q+BHQBVtR84xOAxyBMMHoW8ffil/tTNV24zzCXpHJYN96r6LpBl+hTwuWEVJUlaGz+hKkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoAwedBnDgZNZ4G/XuJutwA+HUM4wnY81gXWt1PlY1/lYE1jXSgyjpn9WVct+xH9s4T4MSaaramrcdcx3PtYE1rVS52Nd52NNYF0rsZ41eVtGkhpkuEtSgzZ6uB8YdwGLOB9rAutaqfOxrvOxJrCulVi3mjb0PXdJ0uI2+pW7JGkR52W4J/mzJG8meW6J7Ulyb5ITSZ5NctW8bdcnOd5t+4N1rOnfdrUcS/K9JB+bt+2Vrv1okulh1dSzrmuTvN0d+2iSu+dtG8m56lnXvnk1PZfkTJJ/0m0byflKclmSx5O8kOT5JJ9fpM84xlafutZ9fPWsa13HV8+axjG2PpDk+0l+0NX1hUX6rO/Yqqrz7gX8CnAV8NwS2/cAjzL4a5VXA0917ZuA/wH8c+DngB8A/2KdaroG+MVu+YazNXXrrwBbx3SurgUeWaR9ZOeqT10L+v468O1Rny/gEuCqbvki4L8v/DePaWz1qWvdx1fPutZ1fPWpaUxjK8AHu+UJBlORXj3OsXVeXrlX1XeA/3OOLktNyL0bOFFVL1fVj4Gvd31HXlNVfa+q/m+3+iSD2ahGrse5WsrIztUq6vo3wF8M69hLqbVN9j7KsbVsXeMYXz3P11JGcr5WUdN6ja2qqlPd6kT3WviG5rqOrfMy3HtYakLuXhN1r4N/x+An9FkFfCvJkSR7x1DPNd2vgY8muaJrOy/OVZKfB65nMI3jWSM/X1n5ZO/rcr7OUdd86z6+lqlrLONruXO13mMryaYMJjR6E3isqsY6tvrMxKQVSPJpBv/zfWpe86eqaibJh4DHkrzUXdmuh6eBHVV1Kske4CBw+Todu49fB/6mquZf5Y/0fGUIk72PQp+6xjG+lqlrLOOr53/DdR1bVXUG+HiSLcBDST5SVYu+57QeNuqV+1ITcq/bRN2LSfJR4MvATVX11tn2qprpvr4JPMTg17B1UVXvnP11saoOARNJtjLmczXPrSz4tXmU5yurn+x9pOerR11jGV/L1TWO8dXnXHXWdWzNO8ZJ4HEGvzXMt75jaxhvJoziBexk6TcJP8v735j4ftd+AfAy8Ev89I2JK9apph0M5pC9ZkH7hcBF85a/B1y/jufqYn76eYbdwP/qzttIz9VydXXbf4HBffkL1+N8df/uB4AvnqPPuo+tnnWt+/jqWde6jq8+NY1pbE0CW7rlzcATwI3jHFvn5W2ZJH/B4F34rUleA/4jgzcoqHNMyF1Vp5P8e+Awg3eg/6yqnl+nmu4G/inwn5MAnK7BHwj6MINf0WDwH/FrVfXNYdTUs65bgDuTnAbmgFtrMKJGdq561gXwr4H/VlV/N+9bR3m+Vj3Z+yjHVs+6xjG++tS13uOrT02w/mPrEuArSTYxuCPyjap6JMkd8+pa17HlJ1QlqUEb9Z67JOkcDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhr0/wBnMHBEgptFtAAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# make a simple scatterplot\n", "x_train = np.array([1,2,3])\n", "y_train = np.array([2,2,4])\n", "plt.scatter(x_train,y_train)\n", "\n", "# check dimensions \n", "print(x_train.shape,y_train.shape)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfgAAAFICAYAAAC84uwKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8HXWd//HXp7mRNvRGi9DUioAUq4v2stx3F1ToIgjVZQUKKu6ylRVELi1QLisXBRTB8lBX7QLFRRa366N0teuDcluuK/TXJly61JaCva/QC20xTdOm+f7+OOeEk+QkmbRJZuYz7+fjkUfozJw533cT+jkz8/3MWAgBERER8WVA3AMQERGR3qcCLyIi4pAKvIiIiEMq8CIiIg6pwIuIiDikAi8iIuKQCrzIPjCzC80smNlJfbDvp81sVbtlD5hZrL2tfZm5L+TH+kDc4xDpbyrwIiJdMLObzGxK3OMQ6SkVeBHpqQeBauDZuAfST74FqMBL6pTHPQARSZcQwh5gT9zjEJGu6QhepHcMMLPpZvammTWZ2Qoz+0qpDc3sM2b2mJltNbOdZvaqmV28L29uZkeZ2SNmtjm/z9fN7GozK4v4+lX5a/5Hmtl/mdl7ZrbNzH5lZge127bkNXgzq8y/58tmtiP/+sVmdmm77YaY2XfNbGX+72qjmT1sZodGHOsD+fcfaWb/ms/cYGZPmtmEKPvI7+ciM6szs8b8WB8zsxOL1h9SNN/hK/n3DHHPgRCJSkfwIr3jNnKnrX8GNAH/CDxgZitDCC8UNjKzacBPgReB7wANwCnAT8zssBDCjJ6+sZlNAp4BdgM/Bv4IfA74LvAJ4PyIu6oFngYeAWbkX/s1YDBwajdjqAQWAicBjwG/AHYCfwZ8AfhRfrshwP8AY4D7gf8FDga+DrxkZpNCCKsjjvdRYAtwE3AQcCnwjJkdF0JY2s14vwtcDSwCrgP2B6YB/21mZ4UQfgtsBL5E7pLEc8DsiOMSSYYQgr70pa+9/AIuBAJQD1QWLa8lV+gfLlp2MLmi928l9nMPudPehxYtexpY1W67B3L/27ZZ9gLQDBxVtMyAufmxfTpCjlX5bb/YbvmP88vHlsh8UtGyq/PLbiux7wHtcjYCn2i3zYeA7cADEcb6QP695gFWtHwi0AI82m77ULxfYGx+u+fb/cxGAVvzfxdlnb1eX/pKy5dO0Yv0jn8OIewq/CGEsB5YAXykaJuzgSrgPjMbUfwF/IbcJbPP9ORNzexA4Hjg1yGEV4veP5A7QwDw+Yi72xBCmNtu2VP57x9pv3E75wPvAre0XxFCaMmP1fLbPQusb5e/gdxZjS7PFLTzvXzOwvssAR4HPmNmNV287ixyH4C+1+5ntgGYQ+7DxvgejEMkkXSKXqR3vFVi2WZyxaLgo/nvT3Sxnw/08H0/nP/+vyXWLSN3pBrp2jadZwA4oJvXfgR4OYSws4ttRub3cyq509+ltHTzPsWWlVj2en7/H6L03wl0/XdWWHYosLgHYxFJHBV4kd7R2axyK/HfXwb+r5PtSxXZ/tLVzHjrYl1UhX08QW5+gIj0IRV4kf7zRv77phBCV0fxPfGH/PePlVh3JLnT/v3xoWEFcKSZVYUQmjrZZiO5a9yDeyn/R8md1i82jtwHla4m6hX+Pj4GvFni9cXbiKSWrsGL9J+55Cbe3Wxm1e1X5tvHqnqywxDCO+RmpX/OzD5etC8DZub/+MjeDzmyh4BhwA3tV+THUrgW/xBwtJmdXWon+TkFUV1d2Hf+tRPIzWF4MoTwpy5e92tyE+dmmFlF0esPBr5K7sNBfdH2fwKG92BcIomgI3iRfhJCWGdm/wjcCywzswfJFZOR5NrJppA7glzVw11/k1yb3HNmVmiTOwOYTG7G/pO9k6BL95BrzbvBzP6cXKvcTnJHyWN5f/Lg9cAJwFwzm0vuCHwXuWvmnwWWkJulH8WHgIVm9mtyHQqXkpuh32WrYQhhuZndSW7m/7Nm9u+83yZXA5wfcjfzKXiR3MS9a4A1uV2EX0Yco0hsVOBF+lEIYY6ZrQCmk+sxHwpsApYDN5Irzj3d52IzOx64mVw/+SByp5ivAe7qpaF3N4ZdZnYqcBUwldx9AXaSuywxp2i7bWZ2Qn67L5Kb0d4MrCPXtnZvD972r4G7yeWuJleIZxR3E3Qx3mvMbCW5v687yH3IeAmYGkJ4rt3mXyfXLng9uQ8CACrwknhW1GXS/29udhO5+zx3pjmEUNHFehHJmPyT4b4SQuiNiX8ibsV9BD8PWFli+VHkTrP9pn+HIyIi4kOsBT5/Kq3D6TQz+1n+P+/r3xGJiIj4kLhZ9GY2CDiX3DW5R2MejoiISColrsADf0vu4RYPtJvJKiJCCOFCXX8X6V6sk+xKMbPnyLXRHBZC+EN324uIiEhHcU+ya8PMxgInkrtRRafFPf/IzWkAgwYNmnjEEUe0X0/h/hctLR1vba31Wq/1Wq/1Wp/W9XV1dZtCCCM7rGwnUQUe+Pv89y57YUMIs8k/m3nSpElh8WI9E0JERLLBzLq6FXOrxFyDN7Nycg/h2EwPbq2ZtEsMvampqYmmps5u651unrOB8qWd53yes0E28kWVmAJP7jaXHwB+0cXDKjrwXOCbm5tpbm6Oexh9wnM2UL6085zPczbIRr6oklTgC6fn1fsuIiKyjxJR4M1sFLn7Si8KIbwW93hERETSLhEFntzTo8ro2YMmREREpBOJKPAhhNtCCBZC+Je4xyIiIuJB0trkemzAgER8RukTgwYNinsIfcZzNlC+tPOcz3M2UL5iqS/wUW3bto1Nmzaxa9euuIciTpSVlbH//vszfPhwqqqq4h6OiEgbqS/wUdrkdu7cydtvv83o0aOprq5uvUtQ0hXuYuTxLEXas4UQ2L17N9u3b2fNmjWMGTOmTZEv9Kp6LfzKl16es0F28kWRzn9di0Qp8Bs3bmTkyJEMHDgwNcUdctm89vmnPZuZUVlZyYgRIxg2bBhbtmxpsz4LvbjKl06es0E28kWV+gIfxc6dO6mpqYl7GOLU4MGDee+99+IehohIG5ko8M3NzZSXp/5qhCRURUUFe/boycYikiyZKPBAqk7NS7rod0tEkigzBV5ERCRLUn/eOq2zsKMoKyuLewh9xnM2UC9u2nnO5zkbKF8xv9VRREQkw1Jf4NPcatWdlpaW1n5xb/Ym28SJEznqqKP6aES9KwvPpFa+dPKcDbKRLyoV+ARLe694V3qabffu3SxdupRJkybt9XvOmTOHe+65Z69f3xNZ6MVVvnTynA2ykS+q1F+Dl2yoqKhg27Zt+9TuePXVV3PMMcfwzW9+sxdHJiKSTCrwkhr77bffXr925cqVbNq0iWOPPbYXRyQiklypP0UvHV1zzTWYGStWrOCyyy6jtraWgQMHcsopp7B27VoAHnzwQSZOnMjAgQMZO3Ys8+fPb339xRdfjJmxYcOGDvtevnw5lZWVXHbZZSXf+6abbsLMePLJJznvvPP4wAc+wMCBAzn66KN59tlnO2z/xhtvcOGFF1JbW0tlZSWHH344d911V4fT94VMmzdvbl02ffp0zIw1a9Zw7bXX8uEPf5jq6momTpzI888/37rdlClT+MhHPgLAjTfeiJlhZtxwww2t22zfvp3vfOc7HHXUUQwZMoTBgwczbtw4Lrnkkih/5SIiiaMjeIfq6+uprq7mC1/4Ap/4xCe48cYbWbJkCffeey+XXHIJY8aM4bnnnuOCCy7AzLjjjjs4//zzWbVqFSNHjuS4447jZz/7GYsWLWLKlClt9n3FFVcwePBgbr755k7fu6ysjPPOO4/jjz+eW2+9lTVr1jBr1ixOO+00Vq5cycEHHwzAY489xtlnn01tbS2XXnopw4YNY8GCBUyfPp3Nmzdz2223tdnvmDFjOOCAA9osGzJkCKeddhrjxo1j+vTpbNy4ke9///v8zd/8DevWraOiooJp06axZ88eFixYwE9+8pPW2xYfd9xxQG7Syl/8xV+wevVqvvrVrzJu3Dh27NjBa6+9xhtvvNGrPxsRkX5TmOyU1q+JEyeG7rz++uvdbrO3HqlbF46//clwyDULwvG3PxkeqVvXZ+8V1YgRIwIQHnrooTbLP/nJTwYgnH766WHXrl2ty2fNmhWA8MQTT4QQQvj9738fgDBz5sw2r1+wYEEAwo9//ONO33vMmDEBCHfddVeb5XPmzAlAmDVrVgghhLfeeisMGjQonHjiiaGhoaHNtsccc0yoqqpqs3zEiBHhrLPOarPd8OHDAxB+/vOft1l+/fXXByCsWLGiddmZZ54ZRo4cWXLMc+fODUBYuHBhp7m605e/YyIixYDFIUJ91Cn6fTC/fj0z573G+q2NBGD91kZmznuN+fXrYxvTunXr2LRpE2eeeSZTp05ts27YsGFUVFRw//33U1FR0bp88ODBAK0T2I444giGDx/OokWLWrfZvXs3V155JR//+Mf52te+VvK93333XdasWcOJJ57IlVde2Wbdpz/9aQBWrVoFwLe//W127NjBvffey8CBA9tse9JJJ9HU1MTq1avbZBo/fnzrNqtXr2bLli2cfvrpfPnLX27z+sJjIqurq1uX1dXVtXl9+3EDLFq0yG1boohkT+oLfIixjezOhctp3N32ISONu/dw58LlvbL/vekVr6+vB+Ccc87psG7p0qWcdNJJHHjggW2WL1u2DICxY8cCuXurH3vssSxevLj17/eee+5hxYoVzJo1q9O70BXe+6KLLiqZBaCmpoaWlhYeeeQRTj755Nb3LFZ4z8Idmwr7LS7QhWXnnntuyZz7778/tbW1QO5xwevWrWPChAklx3322We3XsoYNWoU//AP/8CCBQv2qdhnoRdX+dLJczbIRr6oVOD3wYatjT1a3lOF0yw9UVdXB9BhtvjatWvZuHFjyVnkdXV1jBo1ioMOOqh12bHHHsu2bdtYvnw577zzDrfeeitTpkxpPRIv5eWXXwYo2av+0ksvAbkivW7dOt59910++tGPltzP0qVLGTZsGB/84AfbZCou8IX3KpVnyZIljB8/vvUhMIXXd1bghw8fzpIlS3j00Uc555xzePzxx/nc5z7HiSeeyK5duzrN25Us9OIqXzp5zgbZyBdV6gt8nEYNre7R8v5QX1/P0KFDOfTQQ9ss76rI1dfXd1hemIC2aNEirrvuOpqamrjrrru6fO9C0S3Vq3733XczfPhwTj311NbCW3yZoGDt2rU8/vjjfP7zn2/drr6+nhEjRrQW/MKyIUOGcNhhh7V5/datW3nzzTfb5Ckc7XdW4CF3b/zJkydzzz338Oabb3LBBRfwu9/9jldeeaXLzCIiSaUCvw9mTB5LdUXb09XVFWXMmNzxtHN/KVWsIXdUC7nbvRZbtWoVW7Zs6fCao48+mgEDBnDvvfcyZ84cLr/88g4fGtorFPhnnnmmzfL77ruPl156iRtuuIGamhpGjx7N4MGD27SyATQ2NvKlL32JsrIyZs6c2SZT++vnhZztH9Va6oPMW2+9BcCYMWM6jHnjxo0dzpKUlZVRVlaGmbWe5hcRSRu1ye2DKeNz//jfuXA5G7Y2MmpoNTMmj21d3t+2bNnCmjVr+OIXv9hhXV1dHSNHjmxzFFxYDh2Pbgt94M899xwHHXQQ119/fZfv3dTUxLJlyxg/fjxXXHEFq1ev5pBDDuHpp5/m4Ycf5pxzzuHyyy8Hctf4r7vuOq699lrOOOMMzjjjDLZt28b999/P6tWrmTt3LocffnibTMXX2jdv3szatWtLzjMolafwweSyyy7juOOOo6ysjKlTp2JmTJ8+neeff56zzjqLww8/nJaWFhYuXMiCBQuYMWMGo0aN6jK3iEhSqcDvoynja2Mr6O11dSq6cF26J685+uijWbp0Kbfffjv7779/l++9dOlSmpubufLKK9m6dSt33XUXGzZs4LDDDuMHP/gB3/jGN9ocbV911VUAzJ49m8cee4wDDjiAk08+mXnz5vGxj32sw/hKTbDrLOfAgQM58sgjW5dddtllvP766/zqV7/ipz/9KWPGjOH8888HcrP7N23axNy5c9m4cSPDhw9n3LhxzJ8/n7POOqvLzCIiSWZxTlLrDZMmTQqLFy/ucptly5Z1OqFLStu9ezdHHnlka7tc+1Ph7d13331cdNFFvPLKK6l54ltv0u+YiPQXM1sSQuj2yVs6gpeSvv/97/OHP/yBhx56qNviDrmj6oqKijZHziIiEp/UF/i0n4HoSqEPe8CA/pkLuWXLFhYuXMirr77KnXfeyZVXXhn54Swvv/wyY8eOpbKyMtL2/Z2tvxV6VQs33fFG+dLLczbITr4oVOATrL+zLVy4kKlTp3LggQdyxRVXcMcdd0R6XQiBV199lTPOOCPye3n+ucH7vape/5FRvvTynA2yky+K1Bd46T3nnXce5513Xo9fZ2Zs3769D0YkIiJ7y+f5URERkYxTgRcREXFIBV5ERMSh2K/Bm9lw4DpgCjAaeA9YCvxTCOG57l7vdRY20OlT2zzwnA3efxKeV8qXXp6zgfIVi7XAm9mHgKeBGuA+YAUwBDgKSMbt4URERFIo7iP4X+THcFQI4f/2ZgdR261CCJFu2JIknnvFPWUr9TuYlV5c5UuX+fXr+dETy/jjtp0MqRkY67Mz+orXn11BKp4Hb2Z/CZwIfC+E8H9mVmFmA3u6nygFvqKigsbG3nlGe3/am+fBp4WnbI2NjR3+McnCM6mVL13m169n5rzX2PKnnVSWwfqtjcyc9xrz69fHPbRe5fFnVywtffCfzX9fY2a/AU4DyszsDeCWEMIvouykpaWFhoaGNsvKy8tb/8FtaGhg8ODBrF27ltraWvbbbz8GDBjQeuS4Z8+eDvs0s8SsLxzpJnV8e7u+IKnj6269mdHc3My2bdvYvHkzQ4cObf09LC9//3+r9r+bhfXFv59pXF+Q1PEpX8f1s59aRvWAZqrLAo17cv8PVg9oZvZTyzjliKGxj6+31u/atav1jppJHF9v5IsqzgJfeGj6vwBvAF8BKoGrgAfNrCKEMKfUC81sGjANYPTo0d2+UWFSwoYNG1o//RSKTGcFNAnrW1paSl5aSMr49mV9V0fvSRhflPXl5eWUlZVx4IEHRr5Fr0hc3t6+s0fLJf1ie5qcmT0BfBp4C/hoCGFXfvmw/LKdQG0IoeO/sEUmTJgQCs8A96bw6c3jrFDP2UD50s5jvhPueIr1WxsZXpX7N39LU+5DbO3Qal649lNxDq1XefzZFWtoaKCmpibS0+TinOFUuCj+cKG4A4QQ3gV+DRzE+0f5IiKyD2ZMHkt1Rdv21OqKMmZM1j+zXsV5in5d/vsfS6wrzKgf1t1OPMzC7ozXT6DgOxsoX9p5zFeYLX/nwuVs2NpI7dBql7PoPf7siqWlD34RcDG5m9u0V1j2Tv8NR0TEtynja90VdOlcnIe/88ndte4CM6spLDSzg8nd1W5FCGFldzvx0mpVSlNTU496HtPEczZQvrTznM9zNshGvqhiK/D5a+3Tyd2x7kUzu9LMrgVeJDeb/hsR99N3g4yZ535Oz9lA+dLOcz7P2SAb+aKK9U52IYTZZrYJuBq4FWgBfgdMDSG8EOfYRERE0izuW9USQpgHzIt7HCIiIp74nYIuIiKSYSrwIiIiDsV+in5fqQ8+nTxnA+VLO8/5PGcD5SvmtzqKiIhkWOoLvOc2Oc/9nJ6zgfKlned8nrNBNvJFpQKfYJ77OT1nA+VLO8/5PGeDbOSLKvUFXkRERDpSgRcREXFIBV5ERMQhFXgRERGH1AefYJ77OT1nA+VLO8/5PGcD5SvmtzqKiIhkWOoLvOc2Oc/9nJ6zgfKlned8nrNBNvJFpQKfYJ77OT1nA+VLO8/5PGeDbOSLKvUFXkRERDpSgRcREXFIBV5ERMQhFXgRERGH1AefYJ77OT1nA+VLO8/5PGcD5SvmtzqKiIhkWOoLvOc2Oc/9nJ6zgfKlned8nrNBNvJFpQKfYJ77OT1nA+VLO8/5PGeDbOSLKvUFXkRERDpSgRcREXFIBV5ERMQhFXgRERGH1AefYJ77OT1nA+VLO8/5PGcD5SvmtzqKiIhkWOoLvOc2Oc/9nJ6zgfKlned8nrNBNvJFpQKfYJ77OT1nA+VLO8/5PGeDbOSLKvUFXkRERDqKfZKdmXV2CN4QQqjp18GIiIg4EXuBz3sOmN1u2e44BiIiIuJBUgr8WyGEX8Q9CBERES+SUuAxs0qgMoTwp568Tn3w6eQ5Gyhf2nnO5zkbKF+xpFTHs4EdwHtm9o6Z/dDMhsQ9KBERkbRKwhH8IuA/gJXAYOCzwKXAX5nZ8d0d0e/Zs4eGhoY2y8rLy6mqqgLosC5N65uamtixYweVlZWJHN++rC8o1fKRhPEpn/KBz3wtLS1UVlZSVVWVyPEpX/f5ooq9wIcQjmm36F/N7FXgO8A389/bMLNpwDSA0aNH9/kY49Lc3MyePXviHkaf8NynCsqXdp7zNTc3M2DAgNYi4k0W8kVlSbxRjJlVAH8CloQQju9q2wkTJoS6urr+GVg/K3x683hNyXM2UL6085zPczbIRr6ampolIYRJ3W2blGvwbYQQdgMbgBFxj0VERCSNElngzWw/YDTwdtxjERERSaNYC7yZHdDJqlvJzQ/4TT8OR0RExI24J9ndYGbHAv8NrAFqyM2iPxl4CfhhdztQH3w6ec4Gypd2nvN5zgbKVyzuAv80MA74CnAAsAd4A7geuDuEsDO+oYmIiKRXrAU+hPCfwH/u4z56aTTJU3jur8d2D8/ZQPnSznM+z9kgO/miSP35bc8F3vNzjT1nA+VLO8/5PGeDbOSLKvUFXkRERDpSgRcREXFIBV5ERMQhFXgRERGH4m6T22fqg08nz9lA+dLOcz7P2UD5ivmtjiIiIhmW+gLvuU2uqampRz2PaeI5Gyhf2nnO5zkbZCNfVCrwCea5n9NzNlC+tPOcz3M2yEa+qFJf4EVERKQjFXgRERGHVOBFREQcUoEXERFxSH3wCea5n9NzNlC+tPOcz3M2UL5ifqujiIhIhqW+wHtuk/Pcz+k5Gyhf2nnO5zkbZCNfVCrwCea5n9NzNlC+tPOcz3M2yEa+qFJf4EVERKQjFXgRERGHVOBFREQcUoEXERFxSH3wCea5n9NzNlC+tPOcz3M2UL5ifqujiIhIhqW+wHtuk/Pcz+k5Gyhf2nnO5zkbZCNfVCrwCea5n9NzNlC+tPOcz3M2yEa+qFJf4EVERKQjFXgRERGHVOBFREQcUoEXERFxSH3wCea5n9NzNlC+tPOcz3M2UL5ifqujiIhIhiWqwJvZQDN7y8yCmf0oyms8t8l57uf0nA2UL+085/OcDbKRL6pEFXjgFmBkT17gucB77uf0nA2UL+085/OcDbKRL6puC7yZzTKzwfs0ogjMbAJwOfCtvn4vERER76Icwf8jsNLMvmZm1heDMLMy4F+AR4F5ffEeIiIiWRKlwB8FLAZ+AtSb2Ul9MI4rgCOBS/tg3yIiIpnTbZtcCGE58FkzOx24G3jSzB4BrgohrN7XAZjZh4GbgVtCCKvM7JCevL6lpYWGhoY2y8rLy6mqqgLosC5t6xsbGxM9vr1dX5DU8Smf8iV5fPuyfteuXVRWViZ2fMrXfb6oIk+yCyH8F/Ax4BrgM8AyM/uOme1r0+FPgbfIfXiIxMymmdliM1u8ZcuWfXz75Bo0aBDV1dVxD6NPDBo0yHW/qvKlm+d8AwcOdJsNspEvKtubWehmNhK4HbgQeBu4NoTw4F7s5wLgX4G/DCE8n192CPAH4MchhG5P2U+aNCksXry4p28tIiKSSma2JIQwqbvt9rZNbhjwNPACcDDwgJm9aGZ/3oMBVpE7av8t8EczO9zMDgc+lN9kSH7Z0K7247lNznM/p+dsoHxp5zmf52yQjXxRRWmTO8jMzjSzb5vZY2a2BVgG/Bw4AVgKzAFGAL8zs+9GnG1fTa7n/XTgjaKvp/PrL8j/+aKuduK5wHvu5/ScDZQv7Tzn85wNspEvqij3ot8ABMCALcD/AC8CvwMWhRDeAzCzcmAGuZvVBODabvbbAPxtieUjgX8m1zJ3H/BqhDGKiIhIkSgFfjb5oh5CWNHZRiGEZuD2/E1xLqSbAh9C2A38qv3yoln0b4YQOqwXERGR7kVpk7u4h/t8BfjA3g1HREREekNfPC52ITB1b18cQlhF7nKAiIiI7KVeL/AhhHeBX/b2fjuj58Gnk+dsoHxp5zmf52ygfMX8VkcREZEMS32B99wm57mf03M2UL6085zPczbIRr6oVOATzHM/p+dsoHxp5zmf52yQjXxRpb7Ai4iISEcq8CIiIg6pwIuIiDikAi8iIuJQX9zopl+pDz6dPGcD5Us7z/k8ZwPlK+a3OoqIiGRY6gu85zY5z/2cnrOB8qWd53yes0E28kWlAp9gnvs5PWcD5Us7z/k8Z4Ns5Isq9QVeREREOlKBFxERcUgFXkRExCEVeBEREYfUB59gnvs5PWcD5Us7z/k8ZwPlK+a3OoqIiGRY6gu85zY5z/2cnrOB8qWd53yes0E28kWlAp9gnvs5PWcD5Us7z/k8Z4Ns5Isq9QVeREREOlKBFxERcUgFXkRExCEVeBEREYfUB59gnvs5PWcD5Us7z/k8ZwPlK+a3OoqIiGRY6gu85zY5z/2cnrOB8qWd53yes0E28kWlAp9gnvs5PWcD5Us7z/k8Z4Ns5Isq9QVeREREOlKBFxERcUgFXkRExKFYC7yZjTWzh8xsmZltM7MdZvZ7M7vbzA6Oc2wiIiJpFncf/GjgYOARYB3QDPwZMA0418w+GUJ4p6sdqA8+nTxnA+VLO8/5PGcD5SsWa4EPITwJPNl+uZk9C8wFLgS+18/DEhERSb24j+A7szr/fVh3G3pukyv0O1ZVVcU8kt7nORv4zje/fj0/emIZf9y2kyE1A5kxeSxTxtfGPaxe5fnn5zkbZCdfFIko8Ga2H1AD7AeMA76bX/Xb7l7rucAX+h09/qJ6zgZ+882vX8/Mea9RPaCZyjJYv7WRmfNeA3BV5L3+/MB3NshOvigSUeCBi4AfFv15FXBBCOG57l7Y0tIVEDHAAAAJwklEQVRCQ0NDm2Xl5eWtP9z269K2vrGxMdHj29v1BUkdn/KVXj/7qWUMCO//AzO8KgDNzH5qGaccMTT28fXW+oKkjm9f1u/atYvKysrEjk/5us8XVVIK/Hzg9+SO4scDZwIjOtvYzKaRm4jH6NGj+2N8IgK8vX1nj5aLSHwsiae4zewo4P8BN4UQbu9q2wkTJoS6urr+GVg/K3x68zgr1HM28JvvhDueYv3WxvyRO2xpMgBqh1bzwrWfinNovcrrzw98Z4Ns5KupqVkSQpjU3baJ7DELIbwK1ANfj3ssIvK+GZPHUl1R1mZZdUUZMyaPjWlEItKZpJyiL6UaGN7dRuqDTyfP2cBvvsJEujsXLmfD1kZqh1a7nEXv9ecHvrOB8hWLtcCb2UEhhD+WWH4y8HHg6X4flIh0acr4WncFXcSjuI/gf5K/Je1T5Hrf9wMmAucC7wFXdbeDJM4h6C2e+zk9ZwPlSzvP+Txng+zkiyLuAv8w8GXgS8BIIJAr9D8D7gwhrOluB54LvOd+Ts/ZQPnSznM+z9kgO/miiPtWtXPJ3ZJWREREepHfGWoiIiIZpgIvIiLikAq8iIiIQ3FPsttn6oNPJ8/ZQPnSznM+z9lA+Yr5rY4iIiIZlvoC77lNrqmpqUc9j2niORsoX9p5zuc5G2QjX1Qq8AnW3Nzco57HNPGcDZQv7Tzn85wNspEvqtQXeBEREelIBV5ERMQhFXgRERGHVOBFREQcUh98gnnu5/ScDZQv7Tzn85wNlK+Y3+ooIiKSYakv8J7b5Dz3c3rOBsqXdp7zec4G2cgXlQp8gnnu5/ScDZQv7Tzn85wNspEvqtQXeBEREelIBV5ERMQhFXgRERGHVOBFREQcUh98gnnu5/ScDZQv7Tzn85wNlK+Y3+ooIiKSYakv8J7b5Dz3c3rOBsqXdp7zec4G2cgXlQp8gnnu5/ScDZQv7Tzn85wNspEvqtQXeBEREelIBV5ERMQhFXgRERGHVOBFREQcUh98gnnu5/ScDZQv7Tzn85wNlK+Y3+ooIiKSYakv8J7b5Dz3c3rOBsqXdp7zec4G2cgXlQp8gnnu5/ScDZQv7Tzn85wNspEvqtQXeBEREeko1gJvZkeY2S1m9qKZbTSz98zsZTO73sx8z5QQERHpQ3Efwf8dcAXwJnALMANYDnwb+B8zq45xbCIiIqkVd5vcr4DbQwjbipb91MzeAK4H/h74USwjExERSbFYC3wIYXEnq/6dXIH/eHf7UB98+syvX8+dC5ezYWsjo4ZWM2PyWKaMr417WL3K68+uQPnSy3M2UL5icR/Bd2Z0/vvbsY5Cet38+vXMnPcajbv3ALB+ayMz570G4K7Ii4jEKXEF3szKgBuBZuDfutt+z549NDQ0tFlWXl5OVVUVQId1aVrf1NTEjh07qKysTOT49mb97KeWMSA0Myj/m1dVBtDM7KeWccoRQ2MfX2+tLyjV0pKE8SlfdvO1tLRQWVlJVVVVIsenfN3niypxBR6YBRwHXBdCWF5qAzObBkwDGD16dKlNXGhubmbPnj1xD6NXvb19J1Ao7B2Xe+G5DxeUL82am5sZMGBAaxHxJgv5orIk3SjGzG4FbgBmhxC+FuU1EyZMCHV1dX07sJgUPr15uqZ0wh1PsX5rI8Orcr93W5oMgNqh1bxw7afiHFqv8vizK6Z86eU5G2QjX01NzZIQwqTutk3MDDUzu4lccZ8DXBzvaKSvzJg8luqKtofv1RVlzJg8NqYRiYj4lIhT9Pni/i3g58BFIUmnFaRXFSbSzX5qGW9v30mt01n0IiJxi73Am9k/kSvuDwJ/F0KIPoNAUmnK+NrWCXVeT6OJiMQt1gJvZpcANwNrgCeAqWZWvMnbIYTHu9qH+uDTyXM2UL6085zPczZQvmJxH8H/ef77GHKn59t7BuiywIuIiEhHsR7+hhAuDCFYF18nRdhHP4w0Hp6fa+w5Gyhf2nnO5zkbZCNfVKk/v+25wHt+rrHnbKB8aec5n+dskI18UaW+wIuIiEhHKvAiIiIOqcCLiIg4pAIvIiLiUNxtcvtMffDp5DkbKF/aec7nORsoXzG/1VFERCTDUl/gPbfJee7n9JwNlC/tPOfznA2ykS8qFfgE89zP6TkbKF/aec7nORtkI19UqS/wIiIi0pEKvIiIiEMq8CIiIg6pwIuIiDikPvgE89zP6TkbKF/aec7nORsoXzG/1VFERCTDUl/gPbfJee7n9JwNlC/tPOfznA2ykS8qFfgE89zP6TkbKF/aec7nORtkI19UqS/wIiIi0pEKvIiIiEMq8CIiIg6pwIuIiDikPvgE89zP6TkbKF/aec7nORsoXzG/1VFERCTDUl/gPbfJee7n9JwNlC/tPOfznA2ykS8qFfgE89zP6TkbKF/aec7nORtkI19UqS/wIiIi0pEKvIiIiEMq8CIiIg6pwIuIiDikPvgE89zP6TkbKF/aec7nORsoXzG/1VFERCTDYi3wZjbTzP7DzN4ys2Bmq3q6D89tcp77OT1nA+VLO8/5PGeDbOSLKu4j+NuATwFvAu/uzQ48F3jP/Zyes4HypZ3nfJ6zQTbyRRX3NfjDQghvAZjZUqAm5vGIiIi4EOsRfKG4i4iISO+K+xS9iIiI9IG4T9Hvs5aWFhoaGtosKy8vp6qqCqDDurStb2xsTPT49nZ9QVLHp3zKl+Tx7cv6Xbt2UVlZmdjxKV/3+aJKZYE3s2nAtPwf/1RTU7M8zvGIiIj0ow9F2ciSMgu9MMkuhHBI3GMRERFJO12DFxERcUgFXkRExCEVeBEREYdinWRnZl/i/ckCI4FKM7sh/+fVIYQH4xmZiIhIusU6yc7Mngb+qpPVz4QQTuq/0YiIiPiRmFn0IiIi0nt0DV5EumVm1Wa2zszWmFlVu3X3mtkeMzs3rvGJSEcq8CLSrRBCI/At4IPA1wvLzex24O+Bb4QQfhnT8ESkBJ2iF5FIzKwMeAU4EDgUuAj4AfCtEMItcY5NRDpSgReRyMzsDOA3wFPAycCPQgiXxTsqESlFBV5EesTM6oDxwC+BqUH/iIgkkq7Bi0hkZnYO8In8H99TcRdJLh3Bi0gkZnYqudPzvwF2A38L/FkIYVmsAxORklTgRaRbZnYM8CSwCDgNGA0sA34bQpgS59hEpDSdoheRLpnZOOC3wApgSgihKYTwJnAfcJaZnRDrAEWkJB3Bi0inzGwM8ALQBJwQQni7aN0oYCVQH0JQkRdJGBV4ERERh3SKXkRExCEVeBEREYdU4EVERBxSgRcREXFIBV5ERMQhFXgRERGHVOBFREQcUoEXERFxSAVeRETEIRV4ERERh/4/BkZII7UOgT0AAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def nice_scatterplot(x, y, title):\n", " # font size\n", " f_size = 18\n", " \n", " # make the figure\n", " fig, ax = plt.subplots(1,1, figsize=(8,5)) # Create figure object\n", "\n", " # set axes limits to make the scale nice\n", " ax.set_xlim(np.min(x)-1, np.max(x) + 1)\n", " ax.set_ylim(np.min(y)-1, np.max(y) + 1)\n", "\n", " # adjust size of tickmarks in axes\n", " ax.tick_params(labelsize = f_size)\n", " \n", " # remove tick labels\n", " ax.tick_params(labelbottom=False, bottom=False)\n", " \n", " # adjust size of axis label\n", " ax.set_xlabel(r'$x$', fontsize = f_size)\n", " ax.set_ylabel(r'$y$', fontsize = f_size)\n", " \n", " # set figure title label\n", " ax.set_title(title, fontsize = f_size)\n", "\n", " # you may set up grid with this \n", " ax.grid(True, lw=1.75, ls='--', alpha=0.15)\n", "\n", " # make actual plot (Notice the label argument!)\n", " #ax.scatter(x, y, label=r'$my points$')\n", " #ax.scatter(x, y, label='$my points$')\n", " ax.scatter(x, y, label=r'$my\\,points$')\n", " ax.legend(loc='best', fontsize = f_size);\n", " \n", " return ax\n", "\n", "nice_scatterplot(x_train, y_train, 'hello nice plot')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "#### Formulae\n", "Linear regression is special among the models we study beuase it can be solved explicitly. While most other models (and even some advanced versions of linear regression) must be solved itteratively, linear regression has a formula where you can simply plug in the data.\n", "\n", "For the single predictor case it is:\n", " \\begin{align}\n", " \\beta_1 &= \\frac{\\sum_{i=1}^n{(x_i-\\bar{x})(y_i-\\bar{y})}}{\\sum_{i=1}^n{(x_i-\\bar{x})^2}}\\\\\n", " \\beta_0 &= \\bar{y} - \\beta_1\\bar{x}\\\n", " \\end{align}\n", " \n", "Where $\\bar{y}$ and $\\bar{x}$ are the mean of the y values and the mean of the x values, respectively.\n", "\n", "From the re-aranged second equation we can see that the best-fit line passes through $(\\bar{x},\\bar{y})$, the center of mass of the data\n", "\n", "From any of the first equations, we can see that the slope of the line has to do with whether or not an x value that is above/below the center of mass is typically paired with a y value that is likewise above/below, or typically paired with one that is opposite." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Building a model from scratch\n", "In this part, we will solve the equations for simple linear regression and find the best fit solution to our toy problem." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The snippets of code below implement the linear regression equations on the observed predictors and responses, which we'll call the training data set. Let's walk through the code.\n", "\n", "We have to reshape our arrrays to 2D. We will see later why." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise (5 min)
\n", "* make an array with shape (2,3)\n", "* reshape it to a size that you want" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# your code here\n", "xx = np.array([[1,2,3],[4,6,8]])\n", "xxx = xx.reshape(-1,2)\n", "xxx.shape" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(3, 1)\n" ] } ], "source": [ "# Reshape to be a proper 2D array\n", "x_train = x_train.reshape(x_train.shape[0], 1)\n", "y_train = y_train.reshape(y_train.shape[0], 1)\n", "\n", "print(x_train.shape)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "() ()\n" ] } ], "source": [ "# first, compute means\n", "y_bar = np.mean(y_train)\n", "x_bar = np.mean(x_train)\n", "\n", "# build the two terms\n", "numerator = np.sum( (x_train - x_bar)*(y_train - y_bar) )\n", "denominator = np.sum((x_train - x_bar)**2)\n", "\n", "print(numerator.shape, denominator.shape) #check shapes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Why the empty brackets? (The numerator and denominator are scalars, as expected.)" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The best-fit line is -0.33 + 2.00 * x\n", "The best fit is -0.3333333333333335\n" ] } ], "source": [ "#slope beta1\n", "beta_1 = numerator/denominator\n", "\n", "#intercept beta0\n", "beta_0 = y_bar - beta_1*x_bar\n", "\n", "print(\"The best-fit line is {0:3.2f} + {1:3.2f} * x\".format(beta_0, beta_1))\n", "print(f'The best fit is {beta_0}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise (5 min)
\n", "\n", "Turn the code from the above cells into a function called `simple_linear_regression_fit`, that inputs the training data and returns `beta0` and `beta1`.\n", "\n", "To do this, copy and paste the code from the above cells below and adjust the code as needed, so that the training data becomes the input and the betas become the output.\n", "\n", "```python\n", "def simple_linear_regression_fit(x_train: np.ndarray, y_train: np.ndarray) -> np.ndarray:\n", " \n", " return\n", "```\n", "\n", "Check your function by calling it with the training data from above and printing out the beta values." ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [], "source": [ "def simple_linear_regression_fit(x_train: np.ndarray, y_train: np.ndarray) -> np.ndarray:\n", " \"\"\"\n", " Inputs:\n", " x_train: a (num observations by 1) array holding the values of the predictor variable\n", " y_train: a (num observations by 1) array holding the values of the response variable\n", "\n", " Returns:\n", " beta_vals: a (num_features by 1) array holding the intercept and slope coeficients\n", " \"\"\"\n", " \n", " # Check input array sizes\n", " if len(x_train.shape) < 2:\n", " print(\"Reshaping features array.\")\n", " x_train = x_train.reshape(x_train.shape[0], 1)\n", "\n", " if len(y_train.shape) < 2:\n", " print(\"Reshaping observations array.\")\n", " y_train = y_train.reshape(y_train.shape[0], 1)\n", "\n", " # first, compute means\n", " y_bar = np.mean(y_train)\n", " x_bar = np.mean(x_train)\n", "\n", " # build the two terms\n", " numerator = np.sum( (x_train - x_bar)*(y_train - y_bar) )\n", " denominator = np.sum((x_train - x_bar)**2)\n", " \n", " #slope beta1\n", " beta_1 = numerator/denominator\n", "\n", " #intercept beta0\n", " beta_0 = y_bar - beta_1*x_bar\n", "\n", " return np.array([beta_0,beta_1])\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "* Let's run this function and see the coefficients" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Reshaping features array.\n", "Reshaping observations array.\n", "The best-fit line is 0.666667 + 1.000000 * x\n" ] } ], "source": [ "x_train = np.array([1 ,2, 3])\n", "y_train = np.array([2, 2, 4])\n", "\n", "betas = simple_linear_regression_fit(x_train, y_train)\n", "\n", "beta_0 = betas[0]\n", "beta_1 = betas[1]\n", "\n", "print(\"The best-fit line is {0:8.6f} + {1:8.6f} * x\".format(beta_0, beta_1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise (5 min)
\n", "* Do the values of `beta0` and `beta1` seem reasonable?\n", "* Plot the training data using a scatter plot.\n", "* Plot the best fit line with `beta0` and `beta1` together with the training data." ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAF4CAYAAAD+GDX8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8lOWh/v/rTib7SkggIQshJIDIKpFNkE2sYhVbtdqqtdYWFetST3e/pz2nrefb5XsUtCql2lartYvVai22EgQBFTQoCrIkYd8hQBJC1pm5f39k9JemIAEy8zwz83m/Xnk588zD5LoZ7njlWY21VgAAAHBejNMBAAAA0IFiBgAA4BIUMwAAAJegmAEAALgExQwAAMAlKGYAAAAuQTEDAABwCYoZAACAS1DMAAAAXMLjdIAzlZ2dbYuLi52OAQAAcEpr1qyptdbmnGq9sC1mxcXFqqysdDoGAADAKRljdnRnPXZlAgAAuATFDAAAwCUoZgAAAC5BMQMAAHAJihkAAIBLUMwAAABcgmIGAADgEmF7HTMAAIDuqm1s1XNrdmvTvgY1tHiVnujRkLx0XTOmQL1TE5yO9zGKGQAAiFjv76rTI8tq9PrmQ5KkVq//49cS1+/Xg4urNHVwjuZOLdXIwkynYn6MYgYAACLS029t1/2LNqnF65O1//56S6CkvbrhgJZX1eq+WUN0w4TikGbsimIGAAAiTkcp26jmdv8p17VWam736f5FGyXJ0XLGwf8AACCivL+rTvcv2tStUtZZc7tf9y/apA921wUp2alRzAAAQER5ZFmNWry+M/qzLV6fHl1a08OJui9kxcwYE2uMec8Y8/IJXjPGmIeMMTXGmA+MMeeFKhcAAIgctY2ten3zoRMeU9Yd1kpLNx/S4cbWng3WTaHcYna3pI0nee1SSWWBrzmSHgtVKAAAEDmeW7P7rN/D9ND7nImQFDNjTIGkyyQ9fpJVZkt6ynZYJSnTGJMXimwAACBybNrX8C+XxDgTLV6/Nu0/1kOJTk+otpjNk/QtSSf7m8qXtKvT892BZf/CGDPHGFNpjKk8dOhQz6cEAABhraHF2zPv09zeI+9zuoJezIwxn5Z00Fq75mzfy1q70Fpbbq0tz8nJ6YF0AAAgkqQn9syVwNKT4nrkfU5XKLaYXSDpCmPMdkl/kDTdGPN0l3X2SCrs9LwgsAwAAKDbhuSlK8FzdvUm0ROjIblpPZTo9AS9mFlrv2utLbDWFku6TtJr1tobuqz2kqQvBs7OHC+p3lq7L9jZAABAZLl6TMFZv4ftofc5E45dx8wYc5sx5rbA00WStkqqkfQrSXOdygUAAMJXdmqCpgzOkTFn9ueNkaYNznHsxuYhvSWTtXaZpGWBxws6LbeS7ghlFgAAEJnumFqqFVW1am4//YvMJnpiNXdaaRBSdQ9X/gcAABFlZGGm7ps1RElxp1dzkuJidN+sIRpRkBmkZKfGTcwBAEDE+ehG5Pcv2qQWr+8T7wRgTMeWsvtmDXH0BuYSxQwAAESoGyYUa0Rhph5dWqOlmw/JqOPisR9J9MTIquOYsrnTSh3dUvYRihkAAIhYIwoyteDGch1ubNVza3Zr0/5jamhuV3pSnIbkpunqMQWOHeh/IhQzAAAQ8XqnJujWKQOdjnFKHPwPAADgEhQzAAAAl6CYAQAAuATFDAAAwCUoZgAAAC5BMQMAAHAJihkAAIBLUMwAAABcgmIGAADgEhQzAAAAl6CYAQAAuATFDAAAwCUoZgAAAC5BMQMAAHAJihkAAIBLUMwAAABcgmIGAADgEhQzAAAAl6CYAQAAuATFDAAAwCUoZgAAAC5BMQMAAHAJihkAAIBLUMwAAABcgmIGAADgEhQzAAAAl6CYAQAAuATFDAAAwCUoZgAAAC5BMQMAAHAJihkAAIgKPr91OsIpUcwAAEBEe3fnUd34xGo9uLjK6Sin5HE6AAAAQDC8t/Oo5lVU6/WqQ8pKidfMoX2djnRKFDMAABBxfv7PTXpk6RZlpcTrO5cO0Y3j+yslwf21x/0JAQAAuuGD3XXKTU9Un/RETRvcRykJHt00oTgsCtlHwicpAADACazbXa95FVVasumgvjp5gO67bKjKi7NUXpzldLTTRjEDAABhaf2ees2rqFbFxgPKSIrTNy4epJsmFjsd66xQzAAAQFiaV1Gtt7cd1r0zB+lLFxQrPTHO6UhnjWIGAADCwsZ9DXpoSbW+fckQFWen6Iezz1VKgkcZSeFfyD5CMQMAAK62aX+D5ldU65X1+5WW4NGVo4+pODtF/TKTnI7W4yhmAADAlfx+q7v+8J5e/mCfUhM8umt6qW6ZVKKM5MjZQtYVxQwAALjK3rpm9ctMUkyMUa/keN05vVS3TBqgzOR4p6MFHcUMAAC4QvWBY5q/pFqL1u3TS1+bpGH5GfrRlcOcjhVSFDMAAOComoPH9NCSGv3tg71KjovVbVMGKj8Cjx/rDooZAABwTH1zuz798ErFGKPbpgzUVyeXKCsl8ndZngzFDAAAhNTWQ416Zf1+3TGtVBlJcXroutEa07+XeqcmOB3NcRQzAAAQEttqj+vhJdX669o9ivfEaPaofirolayLz811OpprUMwAAEBQHTzWop++sll/XbtHcbFGt0waoDkXDlROGlvIuqKYAQCAoPD6/PLExsgTE6Nlmw/qSxOLdeuUEvVJS3Q6mmsFvZgZYxIlLZeUEPh+z1lrf9BlnamSXpS0LbDoeWvtD4OdDQAA9LxdR5r08GvVqjnYqL/cPlFZKfF64zvTlRgX63Q01wvFFrNWSdOttY3GmDhJK40xr1hrV3VZb4W19tMhyAMAAIJg15EmPbK0Rs+t2a2YGKMvjC1Sq9evxLhYSlk3Bb2YWWutpMbA07jAlw329wUAAKGzdPNBffXJSsXEGN0wvr9umzJQuRnssjxdITnGzBgTK2mNpFJJj1hrV59gtYnGmA8k7ZH0DWvthyd4nzmS5khSUVFREBMDAIBT2VPXrIMNLRpd1EvnF2fp5guKdcukEgrZWTAdG7RC9M2MyZT0gqQ7rbXrOy1Pl+QP7O6cJWm+tbbsk96rvLzcVlZWBjcwAAD4N3vrmvXI0hr9qXKXBuak6pW7J8sY43QsVzPGrLHWlp9qvZCelWmtrTPGLJV0iaT1nZY3dHq8yBjzqDEm21pbG8p8AADg5PbVN+vRpVv0x3d2ycrqc+WFmjutlFLWg0JxVmaOpPZAKUuSNFPST7uskyvpgLXWGmPGSoqRdDjY2QAAQPdVbDigZ9/eqWvKC3XHtIEq6JXsdKSIE4otZnmSngwcZxYj6U/W2peNMbdJkrV2gaSrJd1ujPFKapZ0nQ3lPlYAAPBvDjS06LFlWzQ0L12fO79Qnzu/UFMH91FhFoUsWEJxVuYHkkafYPmCTo9/IekXwc4CAABO7WBDix57fYt+v3qnvH6r26cMlCQleGIpZUHGlf8BAMDHHl+xVT//52Z5/VafHZ2vO6eXqag3ZSxUKGYAAES5Q8dalZIQq+R4j/qkJ+rykf30tWmlKs5OcTpa1KGYAQAQpWobW/XL17fod6t26J6LBum2KQN1xch+umJkP6ejRS2KGQAAUeZwY6sWLt+qp97aoVavT1eOytfFQ/s6HQuimAEAEHVuf+ZdVW4/oitG9tOdM8o0MCfV6UgIoJgBABDhjhxv069XbtMtkwaoV0q8/s9l5yg53qPSPhQyt6GYAQAQoY4eb9OvVmzVk29uV1O7T4Ny03TFyH4aUZDpdDScBMUMAIAI4/X5NX9JtX7zxnYdb/PqsuF5umtGmQb1TXM6Gk6BYgYAQIRo8/oV74mRJzZGlduPasqgHN01o0yDcylk4YJiBgBAmKtvbtevV27TM6t36uU7Jyk3I1FPfnms4j0xTkfDaaKYAQAQphpaOgrZEyu36ViLV5ecm6t2n1+SKGVhimIGAEAYOnK8TdP+3zLVN7fr4qF9dfdFZTq3X4bTsXCWKGYAAISJYy3temvLYV18bq6yUuJ165QSXViWo2H5FLJIQTEDAMDlGlu9+u0b2/SrFdvU2OrVm9+Zrr7piZo7tdTpaOhhFDMAAFzqeKtXv31zu361Yqvqmto1Y0gf3X1RmfqmJzodDUFCMQMAwKUONLTof1/drCmDcnTPRYM0spALw0Y6ihkAAC7R1ObVU2/t0Pba4/rJVSNUkpOq1785TYVZyU5HQ4hQzAAAcFhTm1dPr9qhX76+VYePt2nKoJyPLxZLKYsuFDMAABz0Zk2t7vrDe6ptbNPksmzdc9Egjenfy+lYcAjFDACAEGtp9+nw8TblZyZpQE6KhuVn6GvTSlVenOV0NDiMYgYAQIi0tPv0zOqdWvD6FpXmpOrZOeOVl5Gk39481ulocAmKGQAAQdbS7tOzb+/UY8u26OCxVk0o6a27LypzOhZciGIGAECQLVy+VQ8srtK4AVmaf91oTRjY2+lIcCmKGQAAPazV69Of3tmlAdmpmlSWrRvH91d5cS9NHJjtdDS4HMUMAIAe0ur16U+Vu/Xo0hrtq2/R9eOKNKksW71S4ill6BaKGQAAPeCF93br5//YrL31LRrTv5d+fvVIXVDKLkucHooZAABnqM3rV2yMUWyM0f76VuVmJOqnV4/QpNJsGWOcjocwFON0AAAAwk27z68/vL1T0/93mV7+YK8k6auTB+gvt0/U5LIcShnOGFvMAADopnafX8+/u1sPv1aj3UebNbIgQ33TEyVJnli2deDsUcwAAOimGx5frdXbjmhEQYZ+OPtcTRvch61j6FEUMwAATsLr8+vlD/bp0uG5SvDE6uYLBmjOhSWaPoRChuCgmAEA0IXX59eLa/fq4deqtf1wk4wZpdmj8nXJsFynoyHCUcwAAAjw+a1een+PHlpSo221xzU0L10LbxyjmUP7Oh0NUYJiBgBAJw8vqVFiXKx+eeMYXTy0L7ssEVIUMwBA1PL5rV7+YK+eXrVDv715rFISPPr9V8erT1qCYmIoZAg9ihkAIOr4/VYvr9unh5ZUq+Zgowb3TdO++maV9klTbkai0/EQxShmAICocvBYi67/1WpVH2zUoL6peuQL5+nSYblsIYMrUMwAABHP77facqhRZX3TlJOaoKH90nXXjDJdNjyPQgZXoZgBACKW32/1zw/3a/6Sau060qSV356uXinxmn/daKejASdEMQMARBy/3+rVDQc0r6JKm/YfU0lOiv7ns8OVnhTndDTgE1HMAAARZ+3uOt329BqVZKdo3rWjdPnIfopllyXCAMUMABD2rLWq2HhQOw4f11cml+i8ol76zc3na3JpNjcXR1ihmAEAwpa1Vks2HtS8JVVav6dBg/qm6qaJxYqLjdG0wX2cjgecNooZACAsvbfzqL7/4odat6deRVnJ+vnVI/SZ0flsIUNYo5gBAMKGtVbN7T4lx3sUFxuj+uZ2/SxQyOIoZIgAFDMAgOtZa7Ws6pDmVVRrYHaKHrh2lIblZ2jpN6ZyUD8iCsUMAOBa1lotr67VvIoqvbezTvmZSbp+bNHHr1PKEGkoZgAA13poSY0erKhSfmaS/uczw3X1mALFe9hlichFMQMAuIa1Vm9uOaw+aQkq65umK0b1U+/UeF1TXqAET6zT8YCgo5gBABxnrdVbWw5rXkW13t5+RJ8rL9DPrh6pAdkpGpCd4nQ8IGQoZgAAR63aelgPLK7S29uOqG96gn44+1xde36h07EAR1DMAACOWrRun7bXHtd/XT5U140tUmIcuywRvShmAICQenvbET24uEp3X1Sm8SW99R8XD9b3Zp1DIQNEMQMAhMg7249oXkWV3qg5rOzUBNU1tUmSMpLiHE4GuEfQi5kxJlHSckkJge/3nLX2B13WMZLmS5olqUnSl6y17wY7GwAg+Ky1uv3pd/WPD/crOzVB/+eyc3T9uP5KimcLGdBVKLaYtUqabq1tNMbESVppjHnFWruq0zqXSioLfI2T9FjgvwCAMLV+T73O7ZcuY4zO65+pMf176YbxFDLgkwS9mFlrraTGwNO4wJftstpsSU8F1l1ljMk0xuRZa/cFOx8AoGe9t/Oo5lVU6/WqQ3ripnLNOKev5lw40OlYQFgIyTFmxphYSWsklUp6xFq7ussq+ZJ2dXq+O7DsX4qZMWaOpDmSVFRUJACAe6zdVad5FVVatvmQeiXH6duXDNH4kt5OxwLCSkiKmbXWJ2mUMSZT0gvGmGHW2vVn8D4LJS2UpPLy8q5b3QAADmlp9+nLv31Hfmv1rUsG66YJxUpJ4Pwy4HSFdNZYa+uMMUslXSKpczHbI6nz1QQLAssAAC61bne9/li5U/99xTAlxsXqV18s1+DcNKVSyIAzFvQ7wRpjcgJbymSMSZI0U9KmLqu9JOmLpsN4SfUcXwYA7rR+T72+8mSlLv/FSv3t/X3acqjjMOIx/XtRyoCzFIoZlCfpycBxZjGS/mStfdkYc5skWWsXSFqkjktl1Kjjchk3hyAXAOA0HDrWqu+9sE6LNxxQeqJH984cpC9dUKz0RK5DBvSUUJyV+YGk0SdYvqDTYyvpjmBnAQCcvmMt7UpLjFN6kke7jjTp6xcN0s2TKGRAMLDNGQBwQpv2N2je4mqt21Ov174xRQmeWC26a7JiYozT0YCIRTEDAPyLzfuPaf6SKi1at19pCR7dPGmAfP6OE+EpZUBwUcwAAB97e9sRXbvwLaXEe3Tn9FLdMmmAMpPjnY4FRA2KGQBEueoDx7TjcJMuGtpXY/r30ncuGaLPlReqVwqFDAg1ihkARKmag8f00JIa/e2DvcrPTNK0IX0UG2N06xRunwQ4hWIGAFFmW+1xzauo0kvv71VSXKxuvXCg5lxYoliOHwMcRzEDgChhrZUxRlsONurVDw9ozoUlmjO5RL1TE5yOBiCAYgYAEW5b7XE9/Fq1CjKTdO/FgzXjnD564zvTlcUxZIDrUMwAIEJtrz2uh1+r0V/X7lFcrNGtF3YcO2aMoZQBLkUxA4AI9PiKrfq/r2ySJ8boSxOLdeuUEvVJS3Q6FoBToJgBQITYdaRJiXGxyklL0MjCTN00oVi3TSlRn3QKGRAuKGYAEOZ2HWnSL16r0V/e3a3rxxXpv2cP0/nFWTq/OMvpaABO0ymLmTFmsaRvWGvfD0EeAEA37T7apEeW1ujPlbsVE2N0w/j+un0q1yADwll3tph9W9I8Y8x2Sd+z1u4LbiQAQHf8+OWNem3TQX1hXJHmTi1Vbga7LIFwZ6y13VvRmKskfV/S85J+Zq1tDmawUykvL7eVlZVORgCAkNpb16xHl9XoyxcMUElOqnYebpIn1qhfZpLT0QCcgjFmjbW2/FTrdesYM2OMkbRZ0mOSfizpq8aY71prf3d2MQEAp7KvvlmPLt2iP76zS1ZWowp7qSQnVUW9k52OBqCHdecYszckDZD0oaRVkr4kaZOku40xk621c4KaEACilLVWP3x5g55ZtVN+a3VNeaHumDZQBb0oZECk6s4WszmSNth/3+d5pzFmYxAyAUBUq29qV0ZynIwxamr16aox+Zo7tVSFWRQyINKdsphZaz/8hJcv68EsABDVDja06NFlW/Ts2zv1l9snalh+hn5y1XB1HE0CIBqc1XXMrLVbeyoIAESrg8datGDZVj2zeoe8fqvPjs5Xr8AtkyhlQHThArMA4KCmNq8u+t/XdbzNp8+MztfXppWqODvF6VgAHEIxA4AQq21s1Svr9+vG8f2VHO/RD2cP06jCTAoZAIoZAITK4cZWLVy+VU+9tUOtXp8mlPRWaZ9UXTk63+loAFyCYgYAQVbf3K7Hlm3RU29tV0u7T1eM7Kc7Z5RpYE6q09EAuAzFDACCxForY4z8fqtnVu/QzKF9def0MpX2oZABODGKGQD0sKPH2/T4yq1as+Oonv3qePVKidfKb01XRnKc09EAuBzFDAB6SF1Tmx5fsU2/fXO7jrd5ddnwPB1v8yk1wUMpA9AtFDMA6AHvbD+im3/zjhpbOwrZXTPKNDg3zelYAMIMxQwAzlB9c7v2HG3W0H7pGpqXrkuG5eorkwdoSG6609EAhCmKGQCcpoaWdv165TY9sXKb+qQlaPHXpyglwaP/d81Ip6MBCHMUMwDopmMt7frNG9v1+Iqtamjx6lPn9tVdM8oUE8NtkwD0DIoZAHTTS+/v1QOLqzRzaF/dPaNMw/IznI4EIMJQzADgJBpbvXryze3KTU/UVWMKdPWYAo0syKSQAQgaihkAdPFRIfvViq2qa2rX58cW6aoxBUrwxFLKAAQVxQwAOvlz5S79z6KNOtrUrmmDc3TPRYM0sjDT6VgAogTFDEDUa2rzysgoKT5W8Z4YjSzM1N0zyjS6qJfT0QBEGYoZgKjV1ObV797aoYXLt+qWyQM0d2qprhjZT7NH5TsdDUCUopgBiDrNbT49vWqHfrl8i2ob2zS5LFsTSnpLkozh0hcAnEMxAxB15j6zRks3H9Kk0mzdc1GZyouznI4EAJIoZgCiQEu7T79fvVNXjs5XVkq8vja9THOnlep8ChkAl6GYAYhYLe0+Pfv2Tj22bIsOHmtVQlyMrh/XX2P6c1A/AHeimAGIOH6/1e9W7dCjy2p0oKFVYwdkaf51ozVhYG+nowHAJ6KYAYgYfr9VTIxRTIzRP9bvV/+sFD147ShNKOnNQf0AwgLFDEDYa/X69KfK3XpixVY9O2e88jKStPCLY5Sa4KGQAQgrFDMAYavN69ef1+zSI6/VaG99i8b076WGZq/yMqS0xDin4wHAaaOYAQhLx1radcm8FdpT16zRRZn6yVUjNLksmy1kAMIaxQxA2Gj3+fXO9iOaODBbaYlx+szofJ0/IEsXUsgARAiKGQDXa/f59fy7u/XwazXaW9esZd+YpqLeyfrGpwY7HQ0AehTFDIBrtfv8euHdPXp4abV2HWnWiIIM/Wj2MBVmJTkdDQCCgmIGwLX2HG3Wd19Yp3Py0vRfN52r6UP6sMsSQESjmAFwDa/PrxfX7tX6vfX6weXnqjg7RS/ecYHO7ZdOIQMQFShmABzn81u99P4ePbSkRttqj2toXrqa2rxKjvdoWH6G0/EAIGQoZgActXZXne7941ptrT2uIblpWnDDGF08tK9iYthCBiD6BL2YGWMKJT0lqa8kK2mhtXZ+l3WmSnpR0rbAouettT8MdjYAzvD5rY42tSk7NUG56YlKS/RowQ3n6eKhuRQyAFEtFFvMvJL+w1r7rjEmTdIaY8xia+2GLuutsNZ+OgR5ADjE57f6+7p9emhJtbJT4/WHOROUm5GoF782yeloAOAKQS9m1tp9kvYFHh8zxmyUlC+pazEDEKH8fqtF6/dpfkW1qg82alDfVN04vljWWg7qB4BOQnqMmTGmWNJoSatP8PJEY8wHkvZI+oa19sMT/Pk5kuZIUlFRUfCCAuhRv35jm378940q7ZOqhz8/WpcNz2OXJQCcQMiKmTEmVdJfJN1jrW3o8vK7koqstY3GmFmS/iqprOt7WGsXSlooSeXl5TbIkQGcIb/f6p8f7ldmcrwmDOytq8cUKCctQZ8e0U+xFDIAOKmQFDNjTJw6Stkz1trnu77euahZaxcZYx41xmRba2tDkQ9Az/D7rV7dsF/zKqq1af8xXTY8TxMG9lZmcrxmj8p3Oh4AuF4ozso0kp6QtNFa+8BJ1smVdMBaa40xYyXFSDoc7GwAes6yzQf1s39s1oZ9DSrJTtG8a0fp8pH9nI4FAGElFFvMLpB0o6R1xpi1gWXfk1QkSdbaBZKulnS7McYrqVnSddZadlUCLmetlbVSTIzRxn3H1NTm1QOfG6krRvaTJzbG6XgAEHZMuPaf8vJyW1lZ6XQMICpZa/XapoOaV1GtWyYN0JWj89Xq9SnWGAoZAJyAMWaNtbb8VOtx5X8A3Wat1dLNHYXsg931KspKVlJ8rCQpwRPrcDoACH8UMwDdNveZd/XK+v0qzErSz64aoc+cl684tpABQI+hmAE4KWutVlTXalxJlhI8sZo1PE9TB+fos+cVUMgAIAgoZgD+jbVWy6trNa+iSu/trNNPPjtc140t4ixLAAgyihmAj320hWxeRZXe3Vmn/Mwk/c9nhuuz5xU4HQ0AogLFDMDH/Fb67799qOY2n3585TBdU17AQf0AEEIUMyCKWWv11pbD+vUb2zTvutFKTfDo8ZvOV7/MRAoZADiAYgZEqbe2HNaDFVV6e9sR9U1P0LZDxzW8IEMDslOcjgYAUYtiBkSZ+qZ23fp0pVZtPaI+aQn6r8uH6rqxRUqMYwsZADiNYgZEib11zeqXmaT0JI9SEzz6/qeH6gvjKGQA4CYUMyDCvbP9iOZVVGntzjqt+PZ0ZaXE6/Gbznc6FgDgBChmQIRas+OIHlxcrZU1tcpOjdfXZw5ScjxbxwDAzShmQATasLdBVz32lnqnxOu+WefohvH9P76nJQDAvShmQIR4b+dRbdjXoOvH9dfQfumaf90ozRzaV8nxTHMACBf8xAbC3NpddZpXUaVlmw8pNz1RV51XoMS4WM0ele90NADAaaKYAWGq5uAx3f/3jVq6+ZB6JcfpW5cM1k0TijnLEgDCGMUMCDNen1+e2Bi1ev1au6tO3/zUYN00sVipCUxnAAh3/CQHwsT6PfWaV1Gt9ESPHrh2lM7tl6G3vjuDLWQAEEEoZoDLfbi3XvMrqvXqhgNKT/To1ikDZa2VMYZSBgARhmIGuNhv39im//rbBqUlevT1iwbp5knFSk+MczoWACBIKGaAy2zc16B4T4wG5qRq6uA+urupXV+eNEAZSRQyAIh0MU4HANBh0/4GzX1mjS6dv0IPLK6SJBVnp+jrMwdRygAgSrDFDHBY1YFjml9Rrb+v26fUBI/unF6qWyYNcDoWAMABFDPAYb9fvVPLNh/U16aV6iuTBygzOd7pSAAAh1DMgBCrOXhMDy2p0RfGFWl8SW/dNaNMd88oU68UChkARDuKGRAiWw416qEl1Xrp/b1KiovVpNJsjS/prSwKGQAggGIGhMB//nW9nlm9QwmeWM25sERzJpeod2qC07EAAC5DMQOCZNeRJhX0SpIxRrkZifrK5BLNubBE2RQyAMBJUMyAHrbj8HE9tKRGf11NfSJbAAATKElEQVS7RwtuGKOZQ/vqjmmlTscCAIQBihnQQ3YebtLDr1Xr+ff2yBNj9KWJxRpZmOF0LABAGKGYAT3A6/Pr6gVvqr65XTdNKNZtU0vUJy3R6VgAgDBDMQPO0K4jTXr27Z26d+YgeWJjNO/aUSrtk6o+6RQyAMCZoZgBp2n30SY9srRGf67crZgYo0+dm6uRhZmaWJrtdDQAQJijmAHdVN/crp+8sknPrdklI6MvjCvS3Kmlys1gCxkAoGdQzIBTaPP6Fe+JUVJcrFZtPazrzi/S3GkDlZeR5HQ0AECEoZgBJ7GvvlmPLK3R8qpaLb73QiV4YvXPey5UvCfG6WgAgAhFMQO62F/fokeX1egPb++S31pdU16olja/EjyxlDIAQFBRzIBO1u2u11UL3pTfb3VNeYHmTi1VYVay07EAAFGCYoaod7ChRVUHGjWpLFtD+6Xrq5MH6LrziyhkAICQo5ghah081qIFy7bqmdU7lJYYp7e+O11xsTH65qeGOB0NABClKGaIOoeOteqXr2/R06t3qN1n9ZnR+bpzeqniYjl+DADgLIoZos66PXX69Rvb9JnRBbpzeqmKs1OcjgQAgCSKGaLA4cZWLVy+VSkJHt01o0zTBvfR69+cxjFkAADXoZghYh053qZfLt+ip97coVavT9eNLZIkGWMoZQAAV6KYISL9uXKXfvDSh2pu9+mKkf105/QylfZJdToWAACfiGKGiFHX1Caf36p3aoL6907RjHP66u4ZpSrtk+Z0NAAAuoVihrBX19Smx1ds02/f3K4rR/fTj68crrEDsjR2QJbT0QAAOC0UM4St+qZ2PbFyq37zxnYda/XqsuF5unF8sdOxAAA4YxSzU6htbNVza3Zr074GNbR4lZ7o0ZC8dF0zpkC9UxOcjhfVvv/Ser24dq9mDc/VXTPKNCQ33elICFPMcwBuYay1Tmc4I+Xl5baysjJo7//+rjo9sqxGr28+JElq9fo/fi3REyMraergHM2dWqqRhZlBy4H/X0NLu369cpsuH9lPA3NStfVQo1q9fp2TRyHDmWGeAwgVY8waa235qdZji9kJPP3Wdt2/aJNavD6dqLe2BH54v7rhgJZX1eq+WUN0w4TikGaMJg0t7frNyu16YuVWNbR4lZYYp4E5qSrJ4SxLnDnmOQA3oph10fHDeqOa2/2nXNdaqbndp/sXbZQkfmgHwaPLarRg2RY1tHg1c2hf3T2jTMPyM5yOhTDHPAfgVhSzTt7fVaf7F23q1g/rzprb/bp/0SaNKMzUiAJ2d5ytlnafEuNiJUlbDx3X2AG9dc9FFDL0DOY5ADfjrs2dPLKsRi1e3xn92RavT48urenhRNGlsdWrR5bWaML/XaJ1u+slST/57HA9flM5pQw9hnkOwM2CvsXMGFMo6SlJfSVZSQuttfO7rGMkzZc0S1KTpC9Za98NdrbOahtb9frmQyc81qQ7rJWWbj6kw42tnMV1mo63evXUWzu0cPkWHW1q1/QhfZQY1/E7gyeW3x3Qc5jnANwuFLsyvZL+w1r7rjEmTdIaY8xia+2GTutcKqks8DVO0mOB/4bMc2t2n/V7mMD73Dpl4NkHihKtXp9mPvC69ta3aOrgHN1z0SCN4uw3BAnzHIDbBb2YWWv3SdoXeHzMGLNRUr6kzsVstqSnbMe1O1YZYzKNMXmBPxsSm/Y1/Mup8meixevXpv3HeihR5Gpq82rxhgOaPSpfCZ5YfW16mYbkpem8ol5OR0OEY54DcLuQHvxvjCmWNFrS6i4v5Uva1en57sCyfylmxpg5kuZIUlFRUY9ma2jx9sz7NLf3yPtEouY2n55etUO/XL5FtY1tGpiTqmH5GfrCuJ79LIGTYZ4DcLuQFTNjTKqkv0i6x1rbcCbvYa1dKGmh1HGB2R6Mp/TEnvmrSE+K65H3iSQt7R2FbMHrW1Xb2KpJpdmcZQlHMM8BuF1IipkxJk4dpewZa+3zJ1hlj6TCTs8LAstCZkheuhLW7z+r3RyJnhgNyU3rwVSRoaXdp3kV1RpRkKFHrz+Pm4vDMcxzAG4X9FPeAmdcPiFpo7X2gZOs9pKkL5oO4yXVh/L4Mkm6ekzBWb+H7aH3CXct7T795o1t+vJv35G1VpnJ8Xr16xfq918dTymDo5jnANwuFFvMLpB0o6R1xpi1gWXfk1QkSdbaBZIWqeNSGTXquFzGzSHI9S+yUxM0ZXCOFm84cEan0hsjTRucE9Wn0Le0+/SHt3fqsde36EBDq8YNyFJdU7t6pcSrX2aS0/EA5jkA1wvFWZkr1XGG+SetYyXdEewsp3LH1FKtqKpVc/vpX3wy0ROrudNKg5AqPKzfU6+vPFmp/Q0tGjsgS/OuHa0JA3s7HQv4N8xzAG7G1Ts7GVmYqftmDVFS3On9tSTFxei+WUOi7jYtrV6fth5qlCSV5KRoeEGGfv+VcfrjnPGUMrgW8xyAm3GvzC4+ukHx/Ys2qcXr+8TdHcZ0/AZ936whUXVj4zavX39es0uPvFajhLhYVdw7RcnxHv3qi+VORwO6hXkOwK0oZidww4RijSjM1KNLa7R08yEZdVxU8iOJnhhZdRxrMndaadT8Bt3m9eu5Nbv1yNIa7alr1uiiTH39okGK+cQd1YA7Mc8BuJGxZ3rTOIeVl5fbysrKoH+fw42tem7Nbm3af0wNze1KT4rTkNw0XT2mIOoOAP5T5S5967kPNKowU1+fOUgXlmWr46RbILwxzwEEmzFmjbX2lLuWKGY4qXafX8+/u1uJcbGaPSpfbV6/Vm09rMkUMgAATkt3ixm7MvFv2n1+vfDeHv3itRrtPNKkmUP7avaofMV7YnThoByn4wEAELEoZvgXizcc0I//vkE7DjdpeH6GnripXNOH9HE6FgAAUYFiBnl9fnn9VolxsWpu9ykt0aPHv1iuGef0YZclAAAhRDGLYj6/1Ytr9+jh12p09ZgC3TGtVJ8enqfLR+RRyAAAcADFLAr5/FZ/e3+vHlpSra21x3VOXrqG9kuXJMVw7QsAABxDMYtCX//jWr30/l4NyU3TghvO08VDcylkAAC4AMUsCvj8VovW7dMFpdnKSonX9eOKdMmwXF1yLoUMAAA3oZhFML/fatH6fZpfUa3qg4367qVDdOuUgRpXwn0sAQBwI4pZBLLW6pX1+zW/olqbDxxTaZ9UPfz50bpseJ7T0QAAwCegmEUgY4yeWb1DXr9fDwUKWSy7LAEAcD2KWQTw+61e3bBfjy3bosduGKN+mUmaf91o9UqOp5ABABBGKGZhzFqrVzcc0LyKam3c16CS7BQdaGhRv8wkZXPjZQAAwg7FLEy1tPt09YI3tX5Pg4p7J+uBz43UFSP7yRMb43Q0AABwhihmYcRaqw/3NmhYfoYS42I1bkBvfWniAF05ikIGAEAkoJiFAWutXtt0UPMqqvXh3notvneKBuak6j8/PdTpaAAAoAdRzFzMWqtlmw9pXkWV3t9dr8KsJP3kqhEqykp2OhoAAAgCipmL7T7arK88Vam8jET99Krh+ux5BYpjlyUAABGLYuYi1lqtqK7Vqq2H9a1LhqgwK1lP3zJOY/r3UryHQgYAQKSjmLmAtVYra2o1r6Jaa3YcVX5mkm6dMlAZSXGaMJDbJwEAEC0oZg6rPnBM331+nSp3HFVeRqJ+fOUwXVNeoARPrNPRAABAiFHMHGCtVWOrV2mJcUpPitPh42360ZXD9DkKGQAAUY1iFmJvbunYZSlJf5wzXn3TE7Xk3imK4dZJAABEPYpZiKzaelgPLq7S6m1H1CctQXOnDpS1kjGilAEAAEkUs5D4w9s79Z3n1yknLUE/uHyoPj+2SIlx7LIEAAD/imIWJJXbj0iSyouzdMmwXDW3+yhkAADgE1HMetiaHUc0r6JaK6prNbksW7+7ZZwyk+N18wUDnI4GAABcjmLWQ9buqtP/vrpZK6pr1TslXvfNOkc3jO/vdCwAABBGKGZnyVorY4xWbz2sD/c26LuXDtGNE/orOZ6/WgAAcHpoD2do7a46zauo0hUj++mz5xXoponFumF8f6Uk8FcKAADODC3iNL0fKGRLNx9Sr+Q4zRqWJ0kc1A8AAM4axew03PfCOj2zeqcyk+P0zU8N1k0Ti5XKFjIAANBDaBWn4fziLOVlJOqmicVKS4xzOg4AAIgwFLPTcOXofKcjAACACBbjdAAAAAB0oJgBAAC4BMUMAADAJShmAAAALkExAwAAcAmKGQAAgEtQzAAAAFyCYgYAAOASFDMAAACXoJgBAAC4BMUMAADAJShmAAAALkExAwAAcAljrXU6wxkxxhyStCPE3zZbUm2Iv6dbROvYo3XcUvSOPVrHLUXv2KN13FL0jt2Jcfe31uacaqWwLWZOMMZUWmvLnc7hhGgde7SOW4resUfruKXoHXu0jluK3rG7edzsygQAAHAJihkAAIBLUMxOz0KnAzgoWscereOWonfs0TpuKXrHHq3jlqJ37K4dN8eYAQAAuARbzAAAAFyCYgYAAOASFDNJxphfG2MOGmPWn+R1Y4x5yBhTY4z5wBhzXqfXLjHGbA689p3Qpe4Z3Rj79YExrzPGvGmMGdnpte2B5WuNMZWhS332ujHuqcaY+sDY1hpjvt/ptUj/zL/ZadzrjTE+Y0xW4LVw/swLjTFLjTEbjDEfGmPuPsE6ETfXuznuSJ3n3Rl7xM31bo47Uud5ojHmbWPM+4Gx//cJ1nH3PLfWRv2XpAslnSdp/UlenyXpFUlG0nhJqwPLYyVtkVQiKV7S+5KGOj2eHh77REm9Ao8v/WjsgefbJWU7PYYgjXuqpJdPsDziP/Mu614u6bUI+czzJJ0XeJwmqarrZxeJc72b447Ued6dsUfcXO/OuLusH0nz3EhKDTyOk7Ra0vgu67h6nrPFTJK1drmkI5+wymxJT9kOqyRlGmPyJI2VVGOt3WqtbZP0h8C6YeNUY7fWvmmtPRp4ukpSQUiCBVk3PvOTifjPvIvPS3o2iHFCxlq7z1r7buDxMUkbJeV3WS3i5np3xh3B87w7n/nJRPRn3kUkzXNrrW0MPI0LfHU9y9HV85xi1j35knZ1er47sOxkyyPVLer4LeMjVlKFMWaNMWaOQ5mCaWJgM/crxphzA8ui5jM3xiRLukTSXzotjojP3BhTLGm0On6b7iyi5/onjLuziJznpxh7xM71U33mkTjPjTGxxpi1kg5KWmytDat57gn1N0R4MsZMU8cP7EmdFk+y1u4xxvSRtNgYsymwNSYSvCupyFrbaIyZJemvksoczhRql0t6w1rbeeta2H/mxphUdfxP6B5rbYPTeUKlO+OO1Hl+irFH7Fzv5r/1iJvn1lqfpFHGmExJLxhjhllrT3hMrRuxxax79kgq7PS8ILDsZMsjijFmhKTHJc221h7+aLm1dk/gvwclvaCOzcARwVrb8NHmcGvtIklxxphsRclnHnCduuzeCPfP3BgTp47/UT1jrX3+BKtE5Fzvxrgjdp6fauyROte785kHRNw8/4i1tk7SUnVsEezM1fOcYtY9L0n6YuBMjvGS6q21+yS9I6nMGDPAGBOvjn/gLzkZtKcZY4okPS/pRmttVaflKcaYtI8eS7pYUtj8RnIqxphcY4wJPB6rjrlyWFHwmUuSMSZD0hRJL3ZaFtafeeDzfELSRmvtAydZLeLmenfGHanzvJtjj7i53s1/65E6z3MCW8pkjEmSNFPSpi6ruXqesytTkjHmWXWcmZNtjNkt6QfqOGBQ1toFkhap4yyOGklNkm4OvOY1xnxN0j/VcTbHr621H4Z8AGehG2P/vqTekh4N/OzyWmvLJfVVxyZiqePf0e+ttf8I+QDOUDfGfbWk240xXknNkq6z1lpJ0fCZS9JnJL1qrT3e6Y+G9Wcu6QJJN0paFzj+RJK+J6lIiui53p1xR+Q8V/fGHolzvTvjliJznudJetIYE6uOkv0na+3LxpjbpPCY59ySCQAAwCXYlQkAAOASFDMAAACXoJgBAAC4BMUMAADAJShmAAAALkExAwAAcAmKGQAAgEtQzABEBWPMUmPMzMDjHxtjHnY6EwB0xZX/AUSLH0j6YeDGzKMlXSFJxphe1tqjjiYDgAC2mAGICtba5ZKMpHvVcdsdX+ClBzuv99F9E7syxvwouAkBgC1mAKKEMWa4Ou6jd9haeyyw7BJJQ4wx/6mOe+f9VdJTxphbJfUKrPtDY0yupDhjTL6kp9VxY+Px1tprnRgLgMjFFjMAEc8YkyfpGUmzJTUGCpkk1aqjaL0j6Vlr7U/V8XPRI6lOHTeDlqRRktZKGqmOmzo/KMkbuhEAiBYUMwARzRiTLOl5Sf9hrd0o6UfqON5MkkZIel8dxWtxYNmPJP1U0pOS9gSWdS5mKwLLbNDDA4g67MoEENGstU2SJnR6vrzT81pJX5FUIOlngWUfSvqmpGxJ7wWWlUmqklQqqcoYky1pf9DDA4g6xlp+6QMAAHADdmUCAAC4BMUMAADAJShmAAAALkExAwAAcAmKGQAAgEtQzAAAAFyCYgYAAOASFDMAAACX+P8AaG02zqs1s6QAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig_scat, ax_scat = plt.subplots(1,1, figsize=(10,6))\n", "\n", "# Plot best-fit line\n", "x_train = np.array([[1, 2, 3]]).T\n", "\n", "best_fit = beta_0 + beta_1 * x_train\n", "\n", "ax_scat.scatter(x_train, y_train, s=300, label='Training Data')\n", "ax_scat.plot(x_train, best_fit, ls='--', label='Best Fit Line')\n", "\n", "ax_scat.set_xlabel(r'$x_{train}$')\n", "ax_scat.set_ylabel(r'$y$');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The values of `beta0` and `beta1` seem roughly reasonable. They capture the positive correlation. The line does appear to be trying to get as close as possible to all the points." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Building a model with `statsmodels` and `sklearn`\n", "\n", "Now that we can concretely fit the training data from scratch, let's learn two `python` packages to do it all for us:\n", "* [statsmodels](http://www.statsmodels.org/stable/regression.html) and \n", "* [scikit-learn (sklearn)](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html).\n", "\n", "Our goal is to show how to implement simple linear regression with these packages. For an important sanity check, we compare the $\\beta$ values from `statsmodels` and `sklearn` to the $\\beta$ values that we found from above with our own implementation.\n", "\n", "For the purposes of this lab, `statsmodels` and `sklearn` do the same thing. More generally though, `statsmodels` tends to be easier for inference \\[finding the values of the slope and intercept and dicussing uncertainty in those values\\], whereas `sklearn` has machine-learning algorithms and is better for prediction \\[guessing y values for a given x value\\]. (Note that both packages make the same guesses, it's just a question of which activity they provide more support for.\n", "\n", "**Note:** `statsmodels` and `sklearn` are different packages! Unless we specify otherwise, you can use either one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below is the code for `statsmodels`. `Statsmodels` does not by default include the column of ones in the $X$ matrix, so we include it manually with `sm.add_constant`." ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "import statsmodels.api as sm" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 1. 1.]\n", " [ 1. 2.]\n", " [ 1. 3.]]\n", "The regression coefficients from the statsmodels package are: beta_0 = 0.666667 and beta_1 = 1.000000\n" ] } ], "source": [ "# create the X matrix by appending a column of ones to x_train\n", "X = sm.add_constant(x_train)\n", "\n", "# this is the same matrix as in our scratch problem!\n", "print(X)\n", "\n", "# build the OLS model (ordinary least squares) from the training data\n", "toyregr_sm = sm.OLS(y_train, X)\n", "\n", "# do the fit and save regression info (parameters, etc) in results_sm\n", "results_sm = toyregr_sm.fit()\n", "\n", "# pull the beta parameters out from results_sm\n", "beta0_sm = results_sm.params[0]\n", "beta1_sm = results_sm.params[1]\n", "\n", "print(\"The regression coefficients from the statsmodels package are: beta_0 = {0:8.6f} and beta_1 = {1:8.6f}\".format(beta0_sm, beta1_sm))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Besides the beta parameters, `results_sm` contains a ton of other potentially useful information." ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " OLS Regression Results \n", "==============================================================================\n", "Dep. Variable: y R-squared: 0.750\n", "Model: OLS Adj. R-squared: 0.500\n", "Method: Least Squares F-statistic: 3.000\n", "Date: Fri, 21 Sep 2018 Prob (F-statistic): 0.333\n", "Time: 11:31:12 Log-Likelihood: -2.0007\n", "No. Observations: 3 AIC: 8.001\n", "Df Residuals: 1 BIC: 6.199\n", "Df Model: 1 \n", "Covariance Type: nonrobust \n", "==============================================================================\n", " coef std err t P>|t| [0.025 0.975]\n", "------------------------------------------------------------------------------\n", "const 0.6667 1.247 0.535 0.687 -15.181 16.514\n", "x1 1.0000 0.577 1.732 0.333 -6.336 8.336\n", "==============================================================================\n", "Omnibus: nan Durbin-Watson: 3.000\n", "Prob(Omnibus): nan Jarque-Bera (JB): 0.531\n", "Skew: -0.707 Prob(JB): 0.767\n", "Kurtosis: 1.500 Cond. No. 6.79\n", "==============================================================================\n", "\n", "Warnings:\n", "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n" ] } ], "source": [ "import warnings\n", "warnings.filterwarnings('ignore')\n", "print(results_sm.summary())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's turn our attention to the `sklearn` library." ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "from sklearn import linear_model" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The regression coefficients from the sklearn package are: beta_0 = 0.666667 and beta_1 = 1.000000\n" ] } ], "source": [ "# build the least squares model\n", "toyregr = linear_model.LinearRegression()\n", "\n", "# save regression info (parameters, etc) in results_skl\n", "results = toyregr.fit(x_train, y_train)\n", "\n", "# pull the beta parameters out from results_skl\n", "beta0_skl = toyregr.intercept_\n", "beta1_skl = toyregr.coef_[0]\n", "\n", "print(\"The regression coefficients from the sklearn package are: beta_0 = {0:8.6f} and beta_1 = {1:8.6f}\".format(beta0_skl, beta1_skl))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We should feel pretty good about ourselves now, and we're ready to move on to a real problem!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The shape of things in `scikit-learn`\n", "Before diving right in to a \"real\" problem, we really ought to discuss more of the details of `sklearn`. We do this now. Along the way, we'll import the real-world dataset." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`Scikit-learn` is the main `python` machine learning library. It consists of many learners which can learn models from data, as well as a lot of utility functions such as `train_test_split`. It can be used in `python` by the incantation `import sklearn`.\n", "\n", "In scikit-learn, an **estimator** is a Python object that implements the methods fit(X, y) and predict(T)\n", "\n", "Let's see the structure of `scikit-learn` needed to make these fits. `.fit` always takes two arguments:\n", "```python\n", " estimator.fit(Xtrain, ytrain)\n", "```\n", "We will consider two estimators in this lab: `LinearRegression` and `KNeighborsRegressor`.\n", "\n", "Critically, `Xtrain` must be in the form of an *array of arrays* (or a 2x2 array) with the inner arrays each corresponding to one sample, and whose elements correspond to the feature values for that sample (visuals coming in a moment).\n", "\n", "`ytrain` on the other hand is a simple array of responses. These are continuous for regression problems.\n", "\n", "![](images/sklearn2.jpg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Practice with `sklearn`\n", "We begin by loading up the `mtcars` dataset and cleaning it up a little bit." ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
car namempgcyldisphpdratwtqsecvsamgearcarb
0Mazda RX421.06160.01103.902.62016.460144
1Mazda RX4 Wag21.06160.01103.902.87517.020144
2Datsun 71022.84108.0933.852.32018.611141
3Hornet 4 Drive21.46258.01103.083.21519.441031
4Hornet Sportabout18.78360.01753.153.44017.020032
\n", "
" ], "text/plain": [ " car name mpg cyl disp hp drat wt qsec vs am gear carb\n", "0 Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4\n", "1 Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4\n", "2 Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1\n", "3 Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1\n", "4 Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2" ] }, "execution_count": 73, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "\n", "#load mtcars\n", "dfcars = pd.read_csv(\"data/mtcars.csv\")\n", "dfcars = dfcars.rename(columns={\"Unnamed: 0\":\"car name\"})\n", "dfcars.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's split the dataset into a training set and test set." ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [], "source": [ "# split into training set and testing set\n", "from sklearn.model_selection import train_test_split\n", "\n", "#set random_state to get the same split every time\n", "traindf, testdf = train_test_split(dfcars, test_size=0.2, random_state=42)" ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Shape of full dataset is: (32, 12)\n", "Shape of training dataset is: (25, 12)\n", "Shape of test dataset is: (7, 12)\n" ] } ], "source": [ "# testing set is around 20% of the total data; training set is around 80%\n", "print(\"Shape of full dataset is: {0}\".format(dfcars.shape))\n", "print(\"Shape of training dataset is: {0}\".format(traindf.shape))\n", "print(\"Shape of test dataset is: {0}\".format(testdf.shape))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we have training and test data. We still need to select a predictor and a response from this dataset. Keep in mind that we need to choose the predictor and response from both the training and test set. You will do this in the exercises below. However, we provide some starter code for you to get things going." ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [], "source": [ "# Extract the response variable that we're interested in\n", "y_train = traindf.mpg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the shape of `y_train`." ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(25,)" ] }, "execution_count": 77, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.shape(y_train)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Another way to see the shape is to use the shape method." ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(25,)" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_train.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is *not* an \"array of arrays\". That's okay! Remember, `sklearn` requires an array of arrays only for the predictor array! You will have to pay close attention to this in the exercises later.\n", "\n", "For now, let's discuss two ways out of this debacle. All we'll do is get `y_train` to be an array of arrays. This doesn't hurt anything because `sklearn` doesn't care too much about the shape of `y_train`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, let's reshape `y_train` to be an array of arrays using the `reshape` method. We want the first dimension of `y_train` to be size $25$ and the second dimension to be size $1$." ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [], "source": [ "y_train_reshape = y_train.values.reshape(y_train.shape[0], 1)" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(25, 1)" ] }, "execution_count": 80, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y_train_reshape.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that `y_train.shape[0]` gives the size of the first dimension.\n", "\n", "There's an even easier way to get the correct shape right from the beginning." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_train_reshape = traindf[['mpg']]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_train_reshape.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, there is a nice shortcut to reshaping an array. `numpy` can infer a dimension based on the other dimensions specified." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "y_train_reshape = y_train.values.reshape(-1,1)\n", "y_train_reshape.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this case, we said the second dimension should be size $1$. Since the requirement of the `reshape()` method is that the requested dimensions be compatible, `numpy` decides the the first dimension must be size $25$.\n", "\n", "What would the `.shape` return if we did `y_train.values.reshape(-1,5)`?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Okay, enough of that. The whole reason we went through that whole process was to show you how to reshape your data into the correct format.\n", "\n", "**IMPORTANT:** Remember that your response variable `ytrain` can be a vector but your predictor variable `xtrain` ***must*** be an array!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Simple linear regression with automobile data\n", "We will now use `sklearn` to predict automobile mileage per gallon (mpg) and evaluate these predictions. We already loaded the data and split them into a training set and a test set." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We need to choose the variables that we think will be good predictors for the dependent variable `mpg`. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise (10 min)
\n", "* Pick one variable to use as a predictor for simple linear regression. Create a markdown cell below and discuss your reasons. \n", "* Justify your choice with some visualizations. \n", "* Is there a second variable you'd like to use? For example, we're not doing multiple linear regression here, but if we were, is there another variable you'd like to include if we were using two predictors?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAF3CAYAAADtkpxQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAHylJREFUeJzt3X9w5Hd93/HXC1mpt0Aiu76SO2E4qKkoscHCikN9QI0xyDhOfDZMBzcF05Ie6XQ8QBkVi2kSKDOxZ8Sv6ZDQmkAxLYWSWhYeY6I6PhMbYxsk6zj5lzA0dsue4zt+CNthQ2X53T/2u+c9IWlXp/3u57vffT5mdqT97Hd339+vvpxffD8/vo4IAQAAoLuelboAAACAfkQIAwAASIAQBgAAkAAhDAAAIAFCGAAAQAKEMAAAgAQIYQAAAAnkFsJsn2j7W7a/Y/s+2x/K2j9ou2r7QPa4MK8aAAAAisp5LdZq25KeHRFP2h6U9A1J75Z0gaQnI+IjuXwxAABADzghrw+Oerp7Mns6mD1Ynh8AAEA5hjBJsj0gaV7SaZL+OCLutv0mSVfYfrukOUnvi4ifbPY5p5xySuzevTvPUgEAADpifn7+hxGxo9V2uXVHHvMl9pCk6yVdIemIpB+qflXsw5J2RsS/XOc9+yTtk6QXvOAFZz3yyCO51wkAALBdtucjYqzVdl2ZHRkRy5JulXRBRDwWEasR8bSkT0s6e4P3XBMRYxExtmNHyzAJAADQU/KcHbkjuwIm2xVJb5D0oO2dTZtdIunevGoAAAAoqjzHhO2UdG02LuxZkr4cETfa/q+2z1S9O/JhSe/KsQYAAIBCynN25EFJo+u0vy2v7wQAAOgVrJgPAACQACEMAAAgAUIYAABAAoQwAACABAhhAAAACRDCAAAAEsj13pHorJmFqqZml3RouaZdQxVNjI9o7+hw6rIAAMBxIIT1iJmFqianF1VbWZUkVZdrmpxelCSCGAAAPYjuyB4xNbt0NIA11FZWNTW7lKgiAACwHYSwHnFoubaldgAAUGyEsB6xa6iypXYAAFBshLAeMTE+osrgwDFtlcEBTYyPJKoIAABsBwPze0Rj8D2zIwEAKAdCWA/ZOzpM6AIAoCTojgQAAEiAEAYAAJAAIQwAACABQhgAAEAChDAAAIAECGEAAAAJEMIAAAASIIQBAAAkQAgDAABIgBAGAACQALctQlfMLFS57yUAAE0IYcjdzEJVk9OLqq2sSpKqyzVNTi9KEkEMANC36I5E7qZml44GsIbayqqmZpcSVQQAQHqEMOTu0HJtS+0AAPQDQhhyt2uosqV2AAD6ASEMuZsYH1FlcOCYtsrggCbGRxJVBABAegzMR+4ag++ZHQkAwDMIYeiKvaPDhC4AAJrQHQkAAJAAIQwAACABQhgAAEAChDAAAIAECGEAAAAJEMIAAAASIIQBAAAkQAgDAABIgBAGAACQACEMAAAgAUIYAABAAoQwAACABAhhAAAACRDCAAAAEiCEAQAAJEAIAwAASIAQBgAAkAAhDAAAIIETUhdQBDMLVU3NLunQck27hiqaGB/R3tHhjr8HAACgoe9D2MxCVZPTi6qtrEqSqss1TU4vStKGoep43gMAANCs77sjp2aXjoaphtrKqqZmlzr6HgAAgGZ9H8IOLde21H687wEAAGiWWwizfaLtb9n+ju37bH8oaz/Z9s22H8p+npRXDe3YNVTZUvvxvgcAAKBZnlfCfi7pvIh4haQzJV1g+1WSrpR0S0S8RNIt2fNkJsZHVBkcOKatMjigifGRjr4HAACgWW4D8yMiJD2ZPR3MHiHpYknnZu3XSvq6pPfnVUcrjYH0W5npeDzvAQAAaOZ6Vsrpw+0BSfOSTpP0xxHxftvLETGUvW5JP2k838jY2FjMzc3lVicAAECn2J6PiLFW2+U6MD8iViPiTEnPl3S27dPXvB6qXx37Bbb32Z6zPXfkyJE8ywQAAOi6rsyOjIhlSbdKukDSY7Z3SlL28/AG77kmIsYiYmzHjh3dKBMAAKBr8pwducN2o9uxIukNkh6UdIOky7PNLpf0lbxqAAAAKKo8V8zfKenabFzYsyR9OSJutH2npC/bfqekRyT90xxrAAAAKKQ8Z0celDS6TvuPJL0+r+8FAADoBX2/Yj4AAEAKhDAAAIAECGEAAAAJEMIAAAASIIQBAAAkQAgDAABIgBAGAACQACEMAAAgAUIYAABAAoQwAACABAhhAAAACRDCAAAAEiCEAQAAJEAIAwAASIAQBgAAkAAhDAAAIAFCGAAAQAKEMAAAgAQIYQAAAAkQwgAAABIghAEAACRACAMAAEiAEAYAAJAAIQwAACABQhgAAEAChDAAAIAETkhdAPI1s1DV1OySDi3XtGuooonxEe0dHU5dFgAAfY8QVmIzC1VNTi+qtrIqSaou1zQ5vShJBDEAABKjO7LEpmaXjgawhtrKqqZmlxJVBAAAGghhJXZoubaldgAA0D2EsBLbNVTZUjsAAOgeQliJTYyPqDI4cExbZXBAE+MjiSoCAAANDMwvscbge2ZHAgBQPISwkts7OkzoAgCggOiOBAAASIAQBgAAkAAhDAAAIAFCGAAAQAKEMAAAgAQIYQAAAAkQwgAAABIghAEAACTAYq3oiJmFKivzAwCwBYQwbNvMQlWT04uqraxKkqrLNU1OL0oSQQwAgA3QHYltm5pdOhrAGmorq5qaXUpUEQAAxUcIw7YdWq5tqR0AABDC0AG7hipbagcAAIQwdMDE+IgqgwPHtFUGBzQxPpKoIgAAio+B+di2xuB7ZkcCANA+Qhg6Yu/oMKELAIAtoDsSAAAgAUIYAABAAoQwAACABAhhAAAACeQWwmyfavtW2/fbvs/2u7P2D9qu2j6QPS7MqwYAAICiynN25FOS3hcR99h+rqR52zdnr308Ij6S43cDAAAUWm4hLCIelfRo9vsTth+QxBoGAAAA6tKYMNu7JY1KujtrusL2QduftX1SN2oAAAAoktxDmO3nSLpO0nsi4nFJn5L0Yklnqn6l7KMbvG+f7Tnbc0eOHMm7TAAAgK7KNYTZHlQ9gH0hIqYlKSIei4jViHha0qclnb3eeyPimogYi4ixHTt25FkmAABA1+U5O9KSPiPpgYj4WFP7zqbNLpF0b141AAAAFFWesyP3SHqbpEXbB7K2D0i6zPaZkkLSw5LelWMNAAAAhZTn7MhvSPI6L92U13cCKKaZhaqmZpd0aLmmXUMVTYyPcMN3AH0vzythAKCZhaompxdVW1mVJFWXa5qcXpQkghiAvsZtiwDkamp26WgAa6itrGpqdilRRQBQDIQwALk6tFzbUjsA9AtCGIBc7RqqbKkdAPoFIQxAribGR1QZHDimrTI4oInxkUQVAUAxMDAfQK4ag++ZHQkAxyKEAcjd3tFhQhcArEF3JAAAQAKEMAAAgAQIYQAAAAkQwgAAABLYMITZfrXttzc9/5+292eP87pTHgAAQDltNjvyQ5KuaHo+Iukdkp4t6QOS9udXFgAAQLlt1h35yxFxf9PzhyJiPiJuk/TcnOsCAAAotc1C2FDzk4i4tOnp8/IpBwAAoD9sFsIetP2baxttXyRpKb+SAAAAym+zMWHvlfRV22+RdE/WdpakcyRdlHdhAAAAZbbhlbCI+J6kl0u6XdLu7HGbpJdHxHe7URwAAEBZtbp35JsknSzpf0XEbBfqAQAA6AubrRP2J6p3Sf49SR+2/ftdqwoAAKDkNrsS9lpJr4iIVdt/V/VuyQ93pyzgGTMLVU3NLunQck27hiqaGB/R3tHh1GVhA/y9AKA9m4Ww/xcRq5IUET+z7S7VBBw1s1DV5PSiaiurkqTqck2T04uSxH/YC4i/FwC0b7MlKl5q+2D2WGx6vmj7YLcKRH+bml06+h/0htrKqqZmWSWliPh7AUD7NrsS9o+6VgWwgUPLtS21Iy3+XgDQvg1DWEQ80s1CgPXsGqqous5/wHcNVRJUg1b4ewFA+zabHfmE7cebHk80/+xmkehfE+MjqgwOHNNWGRzQxPhIooqwGf5eANC+zbojb5H0q5KmJX0pIv5Pd0oCntEYzM1su97A3wsA2ueI2PhF+1ckXSrprZJOlPQ/VA9kP+5OeXVjY2MxNzfXza8EAAA4LrbnI2Ks1XabzY5URPw0Iv6L6ivn/2dJ/0HSOzpSIQAAQB/b9LZFts+RdJmk10j6hqRLIuL2bhQGAABQZhuGMNsPS1qW9CVJ+yQ9lbW/UpIi4p4u1AcAAFBKm10Je1hSSBqX9EZJzSvmh6Tz8isLAACg3DZbJ+zcLtYBAADQVzYdmA8AAIB8EMIAAAAS2DSEue7UbhUDAADQLzZdoiIiwvZNks7oUj1ArmYWqqzmDgAohE1DWOYe278eEd/OvRogRzMLVU1OL6q2sipJqi7XNDm9KEmFCmJlD4pl3z8AaFc7Y8J+Q9Kdtr9v+6DtRdsH8y4M6LSp2aWjAayhtrKqqdmlRBX9okZQrC7XFHomKM4sVFOX1hFl3z8A2Ip2roSN514F0AWHlmtbak9hs6BYhqtFZd8/ANiKllfCIuKRiHhEUk31RVobD6Cn7BqqbKk9hV4IittR9v0DgK1oGcJs/7bthyT9laS/VH0l/a/lXBfQcRPjI6oMDhzTVhkc0MT4SKKKflEvBMXtKPv+AcBWtDMm7MOSXiXpuxHxIkmvl3RXrlUBOdg7OqyrLj1Dw0MVWdLwUEVXXXrGtrrBZhaq2nP1fr3oyq9qz9X7tz226XiCYqdryFMvBGEA6JZ2xoStRMSPbD/L9rMi4lbbn8i9MiAHe0eHOzb2KI/Zlo33tTt7sFdmfDZsdf8AoMwcsfnwLtt/IWmvpKsknSLpsKRfj4hz8i+vbmxsLObm5rr1dUBb9ly9X9V1xjIND1V0x5Xdub99EWoAABzL9nxEjLXarp3uyIsl/UzSeyX9uaTvS/qt7ZUH9L4iDDIvQg0AgOOzYQizfZrtPRHxNxHxdEQ8FRHXSrpH0lD3SgSKqQiDzItQAwDg+Gx2JewTkh5fp/2n2WtAXyvCIPMi1AAAOD6bDcx/XkQsrm2MiEXbu3OrCOgRRRhkXoQaAADHZ8OB+bYfioiXbPDa9yLitFwra8LAfAAA0Cs6MTB/zva/WueDf1fS/HaKAwAA6HebdUe+R9L1tn9Hz4SuMUm/JOmSvAsDAAAosw1DWEQ8Jukc26+TdHrW/NWI2N+VygAAAEqs5Yr5EXGrpFu7UAsAAEDfaGex1uNi+1Tbt9q+3/Z9tt+dtZ9s+2bbD2U/T8qrBgAAgKLKLYRJekrS+yLiZarfAPzf2H6ZpCsl3ZLNvLwlew4AANBXcgthEfFoRNyT/f6EpAckDat+G6Rrs82uVf2+lAAAAH0lzythR2WLu45Kulv1RWAfzV76a0nP60YNAAAARZJ7CLP9HEnXSXpPRBxzG6SorxS77mqxtvfZnrM9d+TIkbzLBAAA6KqWsyO3w/ag6gHsCxExnTU/ZntnRDxqe6ekw+u9NyKukXSNVF8xP886gaKYWahyCyIA6BN5zo60pM9IeiAiPtb00g2SLs9+v1zSV/KqAeglMwtVTU4vqrpcU0iqLtc0Ob2omYVq6tIAADnIsztyj6S3STrP9oHscaGkqyW9wfZDks7PngN9b2p2SbWV1WPaaiurmppdSlQRACBPuXVHRsQ3JHmDl1+f1/cCverQcm1L7QCA3taV2ZEAWts1VNlSOwCgtxHCgIKYGB9RZXDgmLbK4IAmxkcSVQQAyFOusyMBtK8xC5LZkZ1XxlmnZdwnoN8QwoAC2Ts6zH9IO6wx67Qx6aEx61RSzx7rMu4T0I/ojgRQamWcdVrGfQL6ESEMQKmVcdZpGfcJ6EeEMAClVsZZp2XcJ6AfEcIAlFoZZ52WcZ+AfkQIA1Bqe0eH9eazhjXg+trRA7befFZvT4DYOzqsqy49Q8NDFVnS8FBFV116Rk/vE9CPmB0JoNRmFqq6br6q1QhJ0mqErpuvauyFJ/d0aGEmLdD7uBIGoNSYSQigqAhhAEqNmYQAiooQBqDUmEkIoKgIYQBKjZmEAIqKgfkASo17cgIoKkIYgNJjJiGAIiKEAegrMwtVrooBKARCGIC+MbNQ1eT04tElK6rLNU1OL0oSQQxA1zEwH0DfYM0wAEVCCAPQN1gzDECR0B0JoPQa48Big9dZMwxACoQwAKW2dhzYWqwZBiAVQhiAUltvHFjDMLMjASRECANQahuN97KkO648r7vFAEATBuYDKDXuHQmgqAhhAEqNe0cCKCq6IwGUWpHvHcnq/UB/I4QBKL0i3juS1fsB0B0JAAmwej8AQhgAJMDq/QAIYQCQALM2ARDCACABZm0CYGA+ACRQ5FmbALqDEAYAiRRx1iaA7qE7EgAAIAFCGAAAQAKEMAAAgAQYEwYAaAu3WQI6ixAGAGiJ2ywBnUd3JACgJW6zBHQeIQwA0BK3WQI6jxAGAGiJ2ywBnUcIAwC0xG2WgM5jYD4AoCVuswR0HiEMANAWbrMEdBYhDEAuWFMKADZHCAPQcawpBQCtMTAfQMexphQAtEYIA9BxrCkFAK3RHQmg43YNVVRdJ3D1wppSjGUD0C1cCQPQcb26plRjLFt1uabQM2PZZhaqqUsDUEKEMAAdt3d0WFddeoaGhyqypOGhiq669IzCX1FiLBuAbqI7EkAuirSmVLtdjEUZy0aXKNAfCGEASm0ry2UUYSwby3sA/SO37kjbn7V92Pa9TW0ftF21fSB7XJjX9wOAtLUuxiKMZaNLFOgfeV4J+5ykT0r6/Jr2j0fER3L8XgA4aitdjJvdH7FbXYRF6RIFkL/cQlhE3GZ7d16fDwDt2GoX43pj2brZRViELlEA3ZFiduQVtg9m3ZUnJfh+AH2kE12M3ewiLEKXKIDu6HYI+5SkF0s6U9Kjkj660Ya299mesz135MiRbtUHoGQ6sVxGN7sIe3V5DwBb54jI78Pr3ZE3RsTpW3ltrbGxsZibm+t0eQDQlj1X71+3i3B4qKI7rjwvQUUApOIu52J7PiLGWm3X1Sthtnc2Pb1E0r0bbQsARUEXIVA8ZbjDRW4D821/UdK5kk6x/QNJfyjpXNtnSgpJD0t6V17fDwCdstmsSQBpbDZWs1f+t5nn7MjL1mn+TF7fBwB5KtIdAACUYzkX7h0JAAB6zkbLtvTSci6EMABow8xCVXuu3q8XXflV7bl6f0+NOwHKqAxjNbl3JAC0wP0cgeIpw1hNQhgAtFCGAcBAGfX6WE26IwGghTIMAAZQPIQwAGihDAOAARQPIQwAWijDAGAAxcOYMABooQwDgAEUDyEMANrQ6wOAARQP3ZEAAAAJEMIAAAASIIQBAAAkQAgDAABIgIH5ANBjZhaqzNQESoAQBgA9hPtYAuVBdyQA9JDN7mMJoLcQwgCgh3AfS6A8CGEA0EO4jyVQHoQwAOgh3McSKA8G5gNAD+E+lkB5EMIAoMdwH0ugHOiOBAAASIAQBgAAkAAhDAAAIAFCGAAAQAKEMAAAgAQIYQAAAAkQwgAAABIghAEAACRACAMAAEiAEAYAAJAAIQwAACABQhgAAEAChDAAAIAECGEAAAAJEMIAAAASIIQBAAAkQAgDAABIgBAGAACQACEMAAAgAUIYAABAAoQwAACABAhhAAAACRDCAAAAEiCEAQAAJEAIAwAASIAQBgAAkAAhDAAAIAFCGAAAQAKEMAAAgAROSF0AAADAzEJVU7NLOrRc066hiibGR7R3dLjwn70dhDAAAJDUzEJVk9OLqq2sSpKqyzVNTi9K0rbDUp6fvV10RwIAgKSmZpeOhqSG2sqqpmaXCv3Z20UIAwAASR1arm2pvSifvV25hTDbn7V92Pa9TW0n277Z9kPZz5Py+n4AANAbdg1VttRelM/erjyvhH1O0gVr2q6UdEtEvETSLdlzAADQxybGR1QZHDimrTI4oInxkUJ/9nblNjA/Im6zvXtN88WSzs1+v1bS1yW9P68aAABA8TUGyOcxgzHPz94uR0R+H14PYTdGxOnZ8+WIGMp+t6SfNJ5vZmxsLObm5nKrEwAAoFNsz0fEWKvtkg3Mj3r62zAB2t5ne8723JEjR7pYGQAAQP66HcIes71TkrKfhzfaMCKuiYixiBjbsWNH1woEAADohm6HsBskXZ79frmkr3T5+wEAAAohzyUqvijpTkkjtn9g+52Srpb0BtsPSTo/ew4AANB38pwdedkGL70+r+8EAADoFayYDwAAkAAhDAAAIAFCGAAAQAKEMAAAgAQIYQAAAAnketuiTrF9RNIjOX7FKZJ+mOPn9wKOAceg3/df4hhIHIN+33+JYyBt/xi8MCJarjTfEyEsb7bn2rnHU5lxDDgG/b7/EsdA4hj0+/5LHAOpe8eA7kgAAIAECGEAAAAJEMLqrkldQAFwDDgG/b7/EsdA4hj0+/5LHAOpS8eAMWEAAAAJcCUMAAAggb4JYbY/a/uw7Xs3eP1c2z+1fSB7/EG3a8yb7VNt32r7ftv32X73OtvY9n+0/T3bB22/MkWteWhz/0t9Htg+0fa3bH8nOwYfWmeb0p4DUtvHoNTngSTZHrC9YPvGdV4r9TnQ0OIY9MM58LDtxWz/5tZ5vfTnQRvHINfz4IROfljBfU7SJyV9fpNtbo+Ii7pTThJPSXpfRNxj+7mS5m3fHBH3N23zJkkvyR6/IelT2c8yaGf/pXKfBz+XdF5EPGl7UNI3bH8tIu5q2qbM54DU3jGQyn0eSNK7JT0g6ZfXea3s50DDZsdAKv85IEmvi4iN1sPql/Ngs2Mg5Xge9M2VsIi4TdKPU9eRUkQ8GhH3ZL8/ofo/PsNrNrtY0uej7i5JQ7Z3drnUXLS5/6WW/V2fzJ4OZo+1A0NLew5IbR+DUrP9fEm/KelPN9ik1OeA1NYxQB+cB6n1TQhr0znZJdev2f611MXkyfZuSaOS7l7z0rCk/9v0/AcqYVDZZP+lkp8HWRfMAUmHJd0cEX13DrRxDKRynwefkPTvJD29weulPwfU+hhI5T4HpPr/+fgL2/O2963zej+cB62OgZTjedBP3ZGt3CPpBVkXxYWSZlS/BFs6tp8j6TpJ74mIx1PX020t9r/050FErEo60/aQpOttnx4R646VLKs2jkFpzwPbF0k6HBHzts9NXU8KbR6D0p4DTV4dEVXbf1/SzbYfzHqN+kmrY5DrecCVsExEPN7oooiImyQN2j4lcVkdl42BuU7SFyJiep1NqpJObXr+/KytFFrtf7+cB5IUEcuSbpV0wZqXSn0ONNvoGJT8PNgj6bdtPyzpS5LOs/3f1mxT9nOg5TEo+TkgSYqIavbzsKTrJZ29ZpOynwctj0He5wEhLGP7V207+/1s1Y/Nj9JW1VnZ/n1G0gMR8bENNrtB0tuzWTGvkvTTiHi0a0XmqJ39L/t5YHtHdvVHtiuS3iDpwTWblfYckNo7BmU+DyJiMiKeHxG7Jb1V0v6I+OdrNiv1OdDOMSjzOSBJtp+dTVCS7WdLeqOktVfES30etHMM8j4P+qY70vYXJZ0r6RTbP5D0h6oPyFVE/CdJb5H0r20/Jakm6a1RvpVs90h6m6TFbDyMJH1A0guko8fhJkkXSvqepJ9J+hcJ6sxLO/tf9vNgp6RrbQ+o/o/JlyPiRtu/J/XFOSC1dwzKfh78gj47B9bVZ+fA81TvipfqWeC/R8Sf99l50M4xyPU8YMV8AACABOiOBAAASIAQBgAAkAAhDAAAIAFCGAAAQAKEMAAAgAQIYQAKJVuX50u2v5/dSuQm2//wOD/rYtszTc8nbX+v6flv2b6hxWf8qe2Xtdjmc7bfsk77btv/7HhqB1B+hDAAhZEtini9pK9HxD+IiLMkTaq+nk9b77fd/O/aNyW9qun5P5b0eHaLEkk6J9tmQxHxuxFxf7v7sMZuSYQwAOsihAEoktdJWskWSZQkRcR3IuJ228+xfYvte2wv2r5YOnq1acn251Vf7frUpvceUT10nZY1Dat+26pzsufnSLoj+5w32r4z+/w/c/0eo7L9ddtj2e/vtP1d29+y/Wnbn2yq/bW2v2n7fzddFbta0mtsH7D93s4eKgC9jhAGoEhOlzS/wWt/K+mSiHil6mHto43biah+Q90/iYhfi4hH1rzvDknn2B6R9JCku7LnJ0h6haRvu34vuH8v6fzs8+ck/dvmD7G9S9Lvq35lbY+kl675np2SXi3pItXDlyRdKen2iDgzIj7e7kEA0B/65rZFAHqeJf2R7ddKelr1q1qNbspHIuKuDd73TdWveA1IulPStyT9gaRRSQ9GxN/aPl/SyyTdkeW6X8q2bXa2pL+MiB9Lku0/k9Q8Vm0mIp6WdL/ttrpPAfQ3QhiAIrlP9Xu1red3JO2QdFZErNh+WNKJ2Wt/s8ln3iHpCtVD2Kcj4gnbJ6p+L9nGeDBLujkiLttG7T9v+t0bbgUAGbojARTJfkl/x/a+RoPtl9t+jaRfkXQ4C2Cvk/TCNj/zAUm7VO8qXMjaDkj6PWXjwVTvotzTGDtm+9nrzMj8tqR/YvukrCvzzW189xOSnttmnQD6DCEMQGFEREi6RNL52RIV90m6StJfS/qCpDHbi5LeLunBLXzm3ZJ+FBErWfOdkl6s7EpYNoD/HZK+aPtg9vpL13xOVdIfqd6deYekhyX9tMXXH5S0avs7DMwHsJbr/z4BAFqx/ZyIeDK7Ena9pM9GxPWp6wLQm7gSBgDt+6DtA6ovhfFXkmZabA8AG+JKGAAAQAJcCQMAAEiAEAYAAJAAIQwAACABQhgAAEAChDAAAIAECGEAAAAJ/H8sLomQH3gr1gAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAF3CAYAAADtkpxQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAHT5JREFUeJzt3X2QZWldH/Dvj6aVDhgaZAp3WnQ14CQWqzvSQUvUWkFsJCbMrpaVLYNrxWT5w6LAsqZgTKWEYGW3MviSqqgJFMQ1RSBUGBoKjRPcXQKkEJzZXrZ3gXG13A30rLuD2MJCq0Pz5I++PS+73T3TO33uud3386m61fc+9/a9v376TM+3zvNyqrUWAACG60l9FwAAMI6EMACAHghhAAA9EMIAAHoghAEA9EAIAwDogRAGANCDzkJYVT2lqj5RVZ+sqvuq6o2D9jdU1VJV3T24vbyrGgAARlV1tVlrVVWSp7bWHq2qySQfTfKaJC9L8mhr7c2dfDAAwC7w5K7euK2lu0cHDycHN9vzAwCkwxCWJFU1keRkkucm+c3W2ser6seSvLqqfibJiSS/2Fr7q63e51nPela7+uqruywVAGBHnDx58vOttX2Xel1nw5EXfUjVdJL3Jnl1kjNJPp+1s2JvSnJVa+1fbvA9Nye5OUm+5Vu+5QUPPvhg53UCAFypqjrZWpu91OuGsjqytbac5M4kL2utPdxaW22tfS3JW5O8cJPveUtrbba1Nrtv3yXDJADArtLl6sh9gzNgqaqpJC9N8pmquuqCl12f5N6uagAAGFVdzgm7Ksltg3lhT0ry7tbaB6rqv1XVtVkbjnwgyas6rAEAYCR1uTryniQHN2h/ZVefCQCwW9gxHwCgB0IYAEAPhDAAgB4IYQAAPRDCAAB6IIQBAPSg02tHkswvLOXo8VM5vbyS/dNTOTx3IIcOzvRdFgDQMyGsQ/MLSzlybDErZ1eTJEvLKzlybDFJBDEAGHOGIzt09PipcwFs3crZ1Rw9fqqnigCAUSGEdej08sq22gGA8SGEdWj/9NS22gGA8SGEdejw3IFMTU5c1DY1OZHDcwd6qggAGBUm5ndoffK91ZEAwGMJYR07dHBG6AIAHsdwJABAD4QwAIAeCGEAAD0QwgAAeiCEAQD0QAgDAOiBEAYA0AMhDACgB0IYAEAPhDAAgB64bNEQzC8suX4kAHARIaxj8wtLOXJsMStnV5MkS8srOXJsMUkEMQAYY4YjO3b0+KlzAWzdytnVHD1+qqeKAIBRIIR17PTyyrbaAYDxIIR1bP/01LbaAYDxIIR17PDcgUxNTlzUNjU5kcNzB3qqCAAYBSbmd2x98r3VkQDAhYSwITh0cEboAgAuYjgSAKAHQhgAQA+EMACAHghhAAA9EMIAAHoghAEA9EAIAwDogRAGANADIQwAoAdCGABAD4QwAIAeCGEAAD0QwgAAeiCEAQD0QAgDAOiBEAYA0AMhDACgB0IYAEAPntx3AX2bX1jK0eOncnp5Jfunp3J47kAOHZzpuywuwe8NgN1urEPY/MJSjhxbzMrZ1STJ0vJKjhxbTBL/oY8wvzcA9oKxHo48evzUuf/I162cXc3R46d6qojL4fcGwF4w1iHs9PLKttoZDX5vAOwFnYWwqnpKVX2iqj5ZVfdV1RsH7c+sqg9W1f2Dr8/oqoZL2T89ta12RoPfGwB7QZdnwv42yYtba9+d5NokL6uq70vy+iS3t9ael+T2weNeHJ47kKnJiYvapiYncnjuQE8VcTn83gDYCzqbmN9aa0keHTycHNxaklckuW7QfluSDyV5XVd1bGV9ErdVdruL3xsAe0GtZaWO3rxqIsnJJM9N8puttddV1XJrbXrwfCX5q/XHm5mdnW0nTpzorE4AgJ1SVSdba7OXel2nE/Nba6uttWuTfHOSF1bV8x/zfMva2bHHqaqbq+pEVZ04c+ZMl2UCAAzdUFZHttaWk9yZ5GVJHq6qq5Jk8PWRTb7nLa212dba7L59+4ZRJgDA0HS5OnJfVa0PO04leWmSzyR5f5KbBi+7Kcn7uqoBAGBUdblj/lVJbhvMC3tSkne31j5QVR9L8u6q+rkkDyb5qQ5rAAAYSV2ujrwnycEN2v8yyUu6+lwAgN1grHfMBwDoixAGANADIQwAoAdCGABAD4QwAIAeCGEAAD0QwgAAeiCEAQD0QAgDAOiBEAYA0AMhDACgB0IYAEAPhDAAgB4IYQAAPRDCAAB6IIQBAPRACAMA6IEQBgDQAyEMAKAHQhgAQA+EMACAHghhAAA9EMIAAHoghAEA9EAIAwDogRAGANCDJ/ddwDiYX1jK0eOncnp5Jfunp3J47kAOHZzpuywAoEdCWMfmF5Zy5NhiVs6uJkmWlldy5NhikghiADDGDEd27OjxU+cC2LqVs6s5evxUTxUBAKNACOvY6eWVbbUDAONBCOvY/umpbbUDAONBCOvY4bkDmZqcuKhtanIih+cO9FQRADAKTMzv2Prke6sjAYALCWFDcOjgjNAFAFzEcCQAQA+EMACAHghhAAA9EMIAAHoghAEA9EAIAwDogRAGANADIQwAoAc2a2VD8wtLdvkHgA4JYTzO/MJSjhxbzMrZ1STJ0vJKjhxbTBJBDAB2iOFIHufo8VPnAti6lbOrOXr8VE8VAcDeI4TxOKeXV7bVDgBsnxDG4+yfntpWOwCwfUIYj3N47kCmJicuapuanMjhuQM9VQQAe4+J+TzO+uR7qyMBoDtCGBs6dHBG6AKADhmOBADogRAGANADIQwAoAdCGABADzoLYVX1nKq6s6o+VVX3VdVrBu1vqKqlqrp7cHt5VzUAAIyqLldHfjXJL7bW7qqqb0hysqo+OHju11trb+7wswEARlpnIay19lCShwb3v1RVn05izwMAgAxpTlhVXZ3kYJKPD5peXVX3VNXbq+oZw6gBAGCUdB7CquppSd6T5LWttS8m+e0k357k2qydKfvVTb7v5qo6UVUnzpw503WZAABD1WkIq6rJrAWwd7TWjiVJa+3h1tpqa+1rSd6a5IUbfW9r7S2ttdnW2uy+ffu6LBMAYOi6XB1ZSd6W5NOttV+7oP2qC152fZJ7u6oBAGBUdbk68kVJXplksaruHrT9UpIbq+raJC3JA0le1WENAAAjqcvVkR9NUhs89ftdfWaf5heWcvT4qZxeXsn+6akcnjvgAtgAwKa6PBM2NuYXlnLk2GJWzq4mSZaWV3Lk2GKSCGIAwIZctmgHHD1+6lwAW7dydjVHj5/qqSIAYNQJYTvg9PLKttoBAISwHbB/empb7QAAQtgOODx3IFOTExe1TU1O5PDcgZ4qAgBGnYn5O2B98r3VkQDA5RLCdsihgzNCFwBw2QxHAgD0QAgDAOiBEAYA0AMhDACgB5uGsKr6gar6mQse/8+qumNwe/FwygMA2Ju2Wh35xiSvvuDxgSQ/m+SpSX4pyR3dlQUAsLdtNRz591trn7rg8f2ttZOttQ8n+YaO6wIA2NO2CmHTFz5ord1wwcNnd1MOAMB42CqEfaaq/sljG6vqx5Oc6q4kAIC9b6s5Yb+Q5Peq6ieT3DVoe0GS70/y410XBgCwl216Jqy19qdJvivJR5JcPbh9OMl3tdb+ZBjFAQDsVZe6duSPJXlmkv/dWjs+hHoAAMbCVvuE/VbWhiS/McmbqurfDq0qAIA9bqszYT+U5Ltba6tV9feyNiz5puGUxaiZX1jK0eOncnp5JU+fmkxVsvyVs9k/PZXDcwdy6OBM3yUCwK6yVQj7u9baapK01r5SVTWkmhgx8wtLOXJsMStnV5Mkyytnzz23tLySI8cWk0QQA4Bt2GqLin9YVfcMbosXPF6sqnuGVSD9O3r81LkAtpGVs6s5etyuJQCwHVudCftHQ6uCkXZ6eWVHXgMAnLdpCGutPTjMQhhd+6ensnSJkLV/empI1QDA3rDV6sgvVdUXL7h96cKvwyySfh2eO5CpyYlNn5+anMjhuQNDrAgAdr+thiNvT/JNSY4leVdr7f8NpyRGzfqEe6sjAWDnbDUceaiqnp7khiRvraqnJPkfWQtkXxhWgYyGQwdnBC0A2EFbrY5Ma+2vW2v/NWs75/+XJP8uyc8OoS4AgD1ty8sWVdX3J7kxyQ8m+WiS61trHxlGYQAAe9mmIayqHkiynORdSW5O8tVB+/ckSWvtriHUBwCwJ211JuyBJC3JXJIfTXLhjvktyYu7KwsAYG/bamL+dUOsAwBgrGw5MR8AgG4IYQAAPdgyhNWa5wyrGACAcbHlFhWttVZVv5/kmiHVwy40v7B0bjd9O+gDwOXZMoQN3FVV/7i19sedV8OuM7+wlCPHFrNydjVJsrS8kiPHFpPkioLYerBbWl7JRFVWW8uMgAfAHnI5c8K+N8nHqurPquqeqlqsqnu6Lozd4ejxU+cC2LqVs6s5evzUE37P9WC3tLySJFltLcn5gDe/sPTECwaAEXE5Z8LmOq+CXev0IChdbvvl2CjYrVsPeM6GAbDbXfJMWGvtwdbag0lWsrZJ6/oNsn96alvtl+NSAe5KAh4AjIpLhrCq+mdVdX+SP0/yf7K2k/7/6rgudonDcwcyNTlxUdvU5EQOzx14wu95qQB3JQEPAEbF5cwJe1OS70vyJ621b0vykiR/1GlV7BqHDs7klhuuycz0VCrJzPRUbrnhmisaLtwo2K2beFJdUcADgFFxOXPCzrbW/rKqnlRVT2qt3VlVv9F5Zewahw7O7OgcrfX3+jfvXcyX/+7iuWGrX2s58eAXzAkDYNe7nDNhy1X1tCQfTvKOqvqPSb7cbVmMu0MHZ/I3Z7+24XPv/Phnh1wNAOy8ywlhr0jylSS/kOQPkvxZkn/aZVGQnN+a4nLbAWA32TSEVdVzq+pFrbUvt9a+1lr7amvttiR3JZkeXomMq4mqbbUDwG6y1Zmw30jyxQ3a/3rwHHTqxu/d+LKlm7UDwG6y1cT8Z7fWFh/b2FpbrKqrO6sIBn7l0NolS9/58c9mtbVMVOXG733OuXYA2M22CmFbDTnaqImh+JVD1whdAOxJWw1Hnqiqf/3Yxqr6V0lOdlcSAMDet9WZsNcmeW9V/XTOh67ZJF+X5PquCwMA2Ms2DWGttYeTfH9V/XCS5w+af6+1dsdQKgMA2MMuuWN+a+3OJHcOoRYAgLFxOZu1PiFV9ZyqurOqPlVV91XVawbtz6yqD1bV/YOvz+iqBgCAUdVZCEvy1SS/2Fr7zqxdAPznq+o7k7w+ye2ttecluX3wGABgrHQWwlprD7XW7hrc/1KSTyeZydplkG4bvOy2JIe6qgEAYFR1eSbsnMHmrgeTfDxrm8A+NHjqL5I8exg1AACMks5DWFU9Lcl7kry2tXbRZZBaay3Jhldjrqqbq+pEVZ04c+ZM12UCAAzVJVdHXomqmsxaAHtHa+3YoPnhqrqqtfZQVV2V5JGNvre19pYkb0mS2dnZDYMasPvMLyzl6PFTOb28kv3TUzk8dyCHDs70XRbA0HW5OrKSvC3Jp1trv3bBU+9PctPg/k1J3tdVDcBomV9YypFji1laXklLsrS8kiPHFjO/sNR3aQBD1+Vw5IuSvDLJi6vq7sHt5UluTfLSqro/yY8MHgNj4OjxU1k5u3pR28rZ1Rw9fqqnigD609lwZGvto0lqk6df0tXnAqPr9PLKttoB9rKhrI4ESJL901PbagfYy4QwYGgOzx3I1OTERW1TkxM5PHegp4oA+tPp6kgYdVbqDdehgzM58eAX8s6PfzarrWWiKj/xghl9DowlZ8IYW1bqDd/8wlLec3Ipq21t15nV1vKek0v6HBhLQhhjy0q94dPnAOcJYYwtK/WGT58DnCeEMbas1Bs+fQ5wnhDG2LJSb/j0OcB5Qhhj69DBmdxywzWZmZ5KJZmZnsotN1xjpV6HDh2cyU+8YCYTtbaPs9WRwDizRQVj7dBBAWCYNlsdOfutz/R7AMaOM2HA0FgdCXCeEAYMjdWRAOcJYcDQWB0JcJ4QBgyN1ZEA55mYDwzN+uR71+sEEMKAIbMiFWCNEEZv5heWnBEBYGwJYfRifmEpR44tntuuYGl5JUeOLSaJIAbAWDAxn17YLwqAcSeE0Qv7RQEw7gxH0ov901NZ2iBwDXu/KPPSAOiLM2H0YhT2i1qfl7a0vJKW8/PS5heWhlYDAONLCKMXhw7O5JYbrsnM9FQqycz0VG654ZqhnoUyLw2APhmOpDd97xdlXhoAfXImjLHlOoYA9EkIY2yNwrw0AMaX4UjGlusYji6rVoFxIIQx1vqel8bjuZoCMC4MRwIjxapVYFwIYcBIsWoVGBdCGDBSrFoFxoUQBowUq1aBcWFiPjBSrFoFxoUQBowcq1aBcWA4EgCgB0IYAEAPhDAAgB6YEwaMHJctAsaBEAaMFJctAsaF4UhgpLhsETAuhDBgpLhsETAuhDBgpLhsETAuhDBgpLhsETAuTMwHRorLFgHjQggDRo7LFgHjQAgDYCzYf45RI4QBsOfZf45RZGI+AHue/ecYRUIYAHue/ecYRYYjYRczxwUuz/7pqSxtELjsP0efnAmDXWp9jsvS8kpazs9xmV9Y6rs0GDn2n2MUCWGwS5njApfv0MGZ3HLDNZmZnkolmZmeyi03XOPMMb0yHAm71F6e42KYlS7Yf45RI4TBLvX0qcksr5zdsH03s5UAMC46G46sqrdX1SNVde8FbW+oqqWquntwe3lXnw97XdX22ncLw6zAuOjyTNjvJPlPSX73Me2/3lp7c4efC2Nh+SuPPwu2VftusZeHWYH+jOI0h87OhLXWPpzkC129P4y7zZbW7/Yl93v15wL6M6qryftYHfnqqrpnMFz5jB4+H/aEvbrkfq/+XEB/RnWaw7BD2G8n+fYk1yZ5KMmvbvbCqrq5qk5U1YkzZ84Mqz7YNfbqkvu9+nMB/RnVaQ7VWuvuzauuTvKB1trzt/PcY83OzrYTJ07sdHkAY2sU58dAV1506x0bXjFhZnoq//f1L97xz6uqk6212Uu9bqhnwqrqqgseXp/k3s1eC0A3RnV+DHRlVKc5dLY6sqremeS6JM+qqs8l+eUk11XVtUlakgeSvKqrzwdgY1vNj3E2jL1o/bgetbO/nYWw1tqNGzS/ravPA+DyjOr8GOjSKF4xwbUjAcaMbUBgNAhhAGNm2PNj5heW8qJb78i3vf738qJb7zD3DAZcOxJgzAxzfoxrgcLmhDCAMTSs+TEWAcDmDEcC0BmLAGBzQhgAnbEIADYnhAHQmVHdJBNGgTlhAHRmVDfJhFEghAHQqVHcJBNGgeFIAIAeCGEAAD0QwgAAeiCEAQD0wMR8ADo1v7BkdSRsQAgDoDOuHQmbMxwJQGe2unYkjDshDIDOuHYkbE4IA6Azrh0JmxPCAOiMa0fC5kzMB6Azrh0JmxPCAOiUa0fCxgxHAgD0QAgDAOiBEAYA0AMhDACgB0IYAEAPhDAAgB4IYQAAPRDCAAB6IIQBAPRACAMA6IEQBgDQAyEMAKAHQhgAQA+EMACAHghhAAA9EMIAAHoghAEA9EAIAwDogRAGANADIQwAoAdCGABAD4QwAIAeCGEAAD0QwgAAeiCEAQD0QAgDAOiBEAYA0AMhDACgB0IYAEAPntx3AQB7wfzCUo4eP5XTyyvZPz2Vw3MHcujgzNDfA9g9hDCAKzS/sJQjxxazcnY1SbK0vJIjxxaT5LJD1E68B7C7GI4EuEJHj586F57WrZxdzdHjp4b6HsDuIoQBXKHTyyvbau/qPYDdpbMQVlVvr6pHqureC9qeWVUfrKr7B1+f0dXnAwzL/umpbbV39R7A7tLlmbDfSfKyx7S9PsntrbXnJbl98BhgVzs8dyBTkxMXtU1NTuTw3IGhvgewu3Q2Mb+19uGquvoxza9Ict3g/m1JPpTkdV3VADAM6xPnr2Rl4068B7C7VGutuzdfC2EfaK09f/B4ubU2PbhfSf5q/fFWZmdn24kTJzqrEwBgp1TVydba7KVe19vE/LaW/jZNgFV1c1WdqKoTZ86cGWJlAADdG3YIe7iqrkqSwddHNntha+0trbXZ1trsvn37hlYgAMAwDDuEvT/JTYP7NyV535A/HwBgJHS5RcU7k3wsyYGq+lxV/VySW5O8tKruT/Ijg8cAAGOny9WRN27y1Eu6+kwAgN3CjvkAAD0QwgAAeiCEAQD0QAgDAOiBEAYA0INOL1u0U6rqTJIHO3jrZyX5fAfvy3n6eDj083Do5+7p4+HQz9361tbaJXea3xUhrCtVdeJyru3EE6ePh0M/D4d+7p4+Hg79PBoMRwIA9EAIAwDowbiHsLf0XcAY0MfDoZ+HQz93Tx8Ph34eAWM9JwwAoC/jfiYMAKAXYxPCquqBqlqsqrur6sSg7ZlV9cGqun/w9Rl917nbVNXbq+qRqrr3grZN+7WqjlTVn1bVqaqa66fq3WeTfn5DVS0Njum7q+rlFzynn7epqp5TVXdW1aeq6r6qes2g3fG8Q7boY8fyDqqqp1TVJ6rqk4N+fuOg3bE8YsZmOLKqHkgy21r7/AVt/yHJF1prt1bV65M8o7X2ur5q3I2q6oeSPJrkd1trzx+0bdivVfWdSd6Z5IVJ9if5wyTf0Vpb7an8XWOTfn5Dkkdba29+zGv18xNQVVcluaq1dldVfUOSk0kOJfnZOJ53xBZ9/FNxLO+YqqokT22tPVpVk0k+muQ1SW6IY3mkjM2ZsE28Isltg/u3Ze2PAdvQWvtwki88pnmzfn1Fkne11v62tfbnSf40a//ouYRN+nkz+vkJaK091Fq7a3D/S0k+nWQmjucds0Ufb0YfPwFtzaODh5ODW4tjeeSMUwhrSf6wqk5W1c2Dtme31h4a3P+LJM/up7Q9Z7N+nUny2Qte97ls/QeYS3t1Vd0zGK5cH1rQz1eoqq5OcjDJx+N47sRj+jhxLO+oqpqoqruTPJLkg601x/IIGqcQ9gOttWuT/FiSnx8M75zT1sZlx2Nsdoj0a6d+O8m3J7k2yUNJfrXfcvaGqnpakvckeW1r7YsXPud43hkb9LFjeYe11lYH/+d9c5IXVtXzH/O8Y3kEjE0Ia60tDb4+kuS9WTvV+vBgjsL6XIVH+qtwT9msX5eSPOeC133zoI0noLX28OAP7deSvDXnhw/08xM0mD/zniTvaK0dGzQ7nnfQRn3sWO5Oa205yZ1JXhbH8sgZixBWVU8dTAJNVT01yY8muTfJ+5PcNHjZTUne10+Fe85m/fr+JP+8qr6+qr4tyfOSfKKH+vaE9T+mA9dn7ZhO9PMTMpjM/LYkn26t/doFTzmed8hmfexY3llVta+qpgf3p5K8NMln4lgeOU/uu4AheXaS9679+8+Tk/z31tofVNUfJ3l3Vf1ckgeztkKHbaiqdya5LsmzqupzSX45ya3ZoF9ba/dV1buTfCrJV5P8vNU3l2eTfr6uqq7N2pDCA0lelejnK/CiJK9MsjiYS5MkvxTH807arI9vdCzvqKuS3FZVE1k72fLu1toHqupjcSyPlLHZogIAYJSMxXAkAMCoEcIAAHoghAEA9EAIAwDogRAGANADIQzYtarqm6rqXVX1Z4NLkv1+VX3HFbzfh6pq9oLHV1fVvYP711XVX1fV3VX16ar65Z34GYDxJYQBu9Jg48/3JvlQa+0ftNZekORILvMasLVmu38DPzK4FMxskn9RVd+zze8HOEcIA3arH05ytrX2n9cbWmufbK19pKqeVlW3V9VdVbVYVa9Izp3ZOlVVv5u1Xdmfs8l7b6m19uUkJ5M8dwd+DmBMjcuO+cDe8/ysBaGN/E2S61trX6yqZyX5o6p6/+C55yW5qbX2R5t87zuqamVw/+uSfO2xL6iqb0zyfUne9ISrB8aeEAbsRZXk31fVD2UtRM3k/DDlg1sEsCT56dbaiWTtzFmSD1zw3A9W1cLgPW9trd2304UD40MIA3ar+5L85CbP/XSSfUle0Fo7W1UPJHnK4LkvX8FnfqS19uNX8P0A55gTBuxWdyT5+qq6eb2hqr6rqn4wydOTPDIIYD+c5Fv7KhJgM0IYsCu11lqS65P8yGCLivuS3JLkL5K8I8lsVS0m+Zkkn+mvUoCN1drfMQAAhsmZMACAHghhAAA9EMIAAHoghAEA9EAIAwDogRAGANADIQwAoAdCGABAD/4/EKTZ7REE1kYAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "y_mpg = dfcars.mpg\n", "x_wt = dfcars.wt\n", "x_hp = dfcars.hp\n", "\n", "fig_wt, ax_wt = plt.subplots(1,1, figsize=(10,6))\n", "ax_wt.scatter(x_wt, y_mpg)\n", "ax_wt.set_xlabel(r'Car Weight')\n", "ax_wt.set_ylabel(r'Car MPG')\n", "\n", "fig_hp, ax_hp = plt.subplots(1,1, figsize=(10,6))\n", "ax_hp.scatter(x_hp, y_mpg)\n", "ax_hp.set_xlabel(r'Car HP')\n", "ax_hp.set_ylabel(r'Car MPG')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise
\n", "\n", "* Use `sklearn` to fit the training data using simple linear regression.\n", "* Use the model to make mpg predictions on the test set. \n", "* Plot the data and the prediction. \n", "* Print out the mean squared error for the training set and the test set and compare.\n", "\n", "**Hints:**\n", "* Use the following to perform the analysis:\n", "```python\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.metrics import mean_squared_error\n", "```" ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
namempgcyldisphpdratwtqsecvsamgearcarb
0Mazda RX421.06160.01103.902.62016.460144
1Mazda RX4 Wag21.06160.01103.902.87517.020144
2Datsun 71022.84108.0933.852.32018.611141
3Hornet 4 Drive21.46258.01103.083.21519.441031
4Hornet Sportabout18.78360.01753.153.44017.020032
\n", "
" ], "text/plain": [ " name mpg cyl disp hp drat wt qsec vs am gear carb\n", "0 Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4\n", "1 Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4\n", "2 Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1\n", "3 Hornet 4 Drive 21.4 6 258.0 110 3.08 3.215 19.44 1 0 3 1\n", "4 Hornet Sportabout 18.7 8 360.0 175 3.15 3.440 17.02 0 0 3 2" ] }, "execution_count": 83, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sklearn.linear_model import LinearRegression\n", "from sklearn.model_selection import train_test_split\n", "from sklearn.metrics import mean_squared_error\n", "\n", "dfcars = pd.read_csv(\"data/mtcars.csv\")\n", "dfcars = dfcars.rename(columns={\"Unnamed: 0\":\"name\"})\n", "\n", "dfcars.head()" ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [], "source": [ "traindf, testdf = train_test_split(dfcars, test_size=0.2, random_state=42)\n", "\n", "y_train = np.array(traindf.mpg)\n", "X_train = np.array(traindf.wt)\n", "X_train = X_train.reshape(X_train.shape[0], 1)" ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.68797618576\n" ] } ], "source": [ "y_test = np.array(testdf.mpg)\n", "X_test = np.array(testdf.wt)\n", "X_test = X_test.reshape(X_test.shape[0], 1)\n", "\n", "#create linear model\n", "regression = LinearRegression()\n", "\n", "#fit linear model\n", "regression.fit(X_train, y_train)\n", "\n", "predicted_y = regression.predict(X_test)\n", "\n", "r2 = regression.score(X_test, y_test)\n", "print(r2)" ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.770137990979\n", "12.4759856599\n", "7.77369776639\n", "Coefficients: \n", " -5.33694140056 36.9373103135\n" ] } ], "source": [ "print(regression.score(X_train, y_train))\n", "\n", "print(mean_squared_error(predicted_y, y_test))\n", "print(mean_squared_error(y_train, regression.predict(X_train)))\n", "\n", "print('Coefficients: \\n', regression.coef_[0], regression.intercept_)" ] }, { "cell_type": "code", "execution_count": 87, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmIAAAF3CAYAAAAGpSdTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmczXX///Hn2zSyRIiEkouQJQ0mkVBSZMkgspOkvuoqZZ0wM2axjX3fScmSxJXLRXJxVVeRYeyMLaohZF9mmOX9+8Phpy7LYM58zpzzuN9uc5tzPnOWpz4+5tk578/rGGutAAAAkPGyOB0AAADAV1HEAAAAHEIRAwAAcAhFDAAAwCEUMQAAAIdQxAAAABxCEQMAAHAIRQwAAMAhbitixphsxpifjDGbjTHbjTEDXNvDjDHxxphNrq/67soAAADgyYy7JusbY4yknNbac8YYf0nfS3pfUj1J56y1w9zyxAAAAJnEPe56YHu54Z1zXfV3fd1R68ufP78tVqxYOiUDAABwnw0bNvxhrS2Qltu6rYhJkjHGT9IGSY9JGm+tXWeMeVnS340x7SXFSOpurT15s8cpVqyYYmJi3BkVAAAgXRhjDqb1tm5drG+tTbHWBkh6WFIVY0x5SRMlFZcUIOmwpOHXu68xposxJsYYE3Ps2DF3xgQAAHBEhpw1aa09JWm1pHrW2iOugpYqaaqkKje4zxRrbaC1NrBAgTS9ugcAAJCpuPOsyQLGmDyuy9klvShplzGm0DU3ayJpm7syAAAAeDJ3rhErJOlj1zqxLJIWWGuXGmM+McYE6PLC/QOS3nJjBgAAAI/lzrMmt0iqeJ3t7dz1nAAAAJkJk/UBAAAcQhEDAABwCEUMAADAIRQxAAAAh1DEAAAAHEIRAwAAcAhFDAAAeL3U1FQtXLhQCQkJTkf5E4oYAADwWtZaLV++XIGBgWrevLk+/fRTpyP9CUUMAAB4pR9//FHPP/+8Xn75ZZ08eVKffPKJOnXq5HSsP6GIAQAAr7Jt2zY1btxYzzzzjHbt2qVx48YpLi5Obdu2lZ+fn9Px/oQiBgAAvMKBAwfUoUMHVahQQWvWrFFUVJT27dund955R1mzZnU63nW580O/AQAA3O7IkSOKiorSpEmT5Ofnpx49eqh379564IEHnI52SxQxAACQKZ0+fVrDhw/XiBEjlJiYqE6dOikkJEQPP/yw09HSjCIGAAAylYSEBE2YMEEDBw7UiRMn1Lx5c0VERKh06dJOR7ttrBEDAACZQnJysqZPn65SpUqpR48eCgwMVExMjBYsWJApS5hEEQMAAB7OWquFCxeqfPny6ty5s4oUKaJ///vfWrFihSpXrux0vLtCEQMAAB7JWquVK1eqSpUqat68ufz8/PTll19enQ/mDShiAADA4/z000+qU6eOXnrpJR07dkyzZs3Sli1bFBQUJGOM0/HSDUUMAAB4jJ07d6pp06Z6+umntXXrVo0ePVpxcXHq0KGDxw1jTQ+cNQkAABz3yy+/KCwsTB9//LFy5syp8PBwdevWTbly5XI6mltRxAAAgGOOHTumQYMGafz48TLGqFu3bgoODlb+/PmdjpYhKGIAACDDnT17ViNGjNDw4cN1/vx5dezYUaGhoSpatKjT0TIURQwAAGSYixcvatKkSYqMjNQff/yhZs2aKSIiQmXKlHE6miNYrA8AANwuOTlZM2fOVKlSpdStWzcFBATop59+0sKFC322hEkUMQAA4EbWWn355ZeqUKGCOnXqpAcffFArV67UypUr9dRTTzkdz3EUMQAA4Bb//ve/VbVqVTVt2lSpqan64osvrs4Hw2UUMQAAkK5iYmL00ksv6YUXXtChQ4c0ffp0bdu2TU2bNvWqYazpgSIGAADSRVxcnFq0aKGnnnpKGzdu1PDhw7Vnzx516tRJ99zD+YHXw38VAABwV3777TcNGDBAM2fOVLZs2RQSEqLu3bsrd+7cTkfzeBQxAABwR44fP65BgwZp3LhxSk1N1bvvvquPPvpIDz74oNPRMg2KGAAAuC3nzp3TqFGjFB0drXPnzqldu3YKCwtTsWLFnI6W6VDEAABAmly6dElTpkxRRESEjh49qqCgIEVGRqpcuXJOR8u0KGIAAOCmUlJS9NlnnykkJEQHDhzQc889pyVLlqhq1apOR8v0OGsSAABcl7VWX331lQICAtS+fXvly5dPK1asuDofDHePIgYAAP7Ht99+q+rVq+uVV17RxYsXNX/+fK1fv14vvfQSs8DSEUUMAABctWnTJtWvX1+1atXSwYMHNXnyZG3fvl0tWrRQlizUhvTGf1EAAKC9e/eqVatWqlixotauXauhQ4dq79696tKli/z9/Z2O57VYrA8AgA87dOiQwsPDNX36dGXNmlV9+/ZVjx49lCdPHqej+QSKGAAAPujEiRMaMmSIxowZo5SUFL399tvq27evHnroIaej+RSKGAAAPuT8+fMaM2aMhgwZojNnzqhNmzYaMGCAihcv7nQ0n0QRAwDAB1y6dEnTpk1TeHi4jhw5ooYNGyoqKkoVKlRwOppPo4gBAODFUlNTNW/ePPXv31/79+9XjRo19MUXX6h69epOR4M4axIAAK9krdWyZctUqVIltWnTRrly5dKyZcv0n//8hxLmQShiAAB4me+//141a9ZUgwYNdO7cOX322WfauHGjXn75ZYaxehiKGAAAXmLLli1q1KiRatSoob1792rixInauXOnWrVqxTBWD8VeAQAgk9u/f7/atm2rgIAAff/99xo8eLD27dunt99+m2GsHo7F+gAAZFK///67IiIiNGXKFPn7+6t3797q1auX8ubN63Q0pBFFDACATObUqVOKjo7WqFGjdOnSJb355pvq37+/ChUq5HQ03Ca3FTFjTDZJ30q61/U8C621ocaYfJLmSyom6YCkFtbak+7KAQCAt7hw4YLGjRunwYMH6+TJk2rZsqUiIiL02GOPOR0Nd8ida8QuSqptrX1SUoCkesaYqpL6SFplrS0paZXrOgAAuIGkpCRNmTJFJUuWVO/evVWtWjXFxsZq7ty5lLBMzm1FzF52znXV3/VlJTWW9LFr+8eSgtyVAQCAzCw1NVXz589XuXLl9NZbb6lYsWL6z3/+o3/+858KCAhwOh7SgVvPmjTG+BljNkk6KmmltXadpILW2sOum/wuqaA7MwAAkNlYa7V8+XIFBgaqZcuWuvfee/WPf/zj6nwweA+3FjFrbYq1NkDSw5KqGGPK/+XnVpdfJfsfxpguxpgYY0zMsWPH3BkTAACP8eOPP+r555/Xyy+/rJMnT+qTTz7Rpk2b1KhRI4axeqEMmSNmrT0labWkepKOGGMKSZLr+9Eb3GeKtTbQWhtYoECBjIgJAIBjtm/frqCgID3zzDPauXOnxo4dq7i4OLVt21Z+fn5Ox4ObuK2IGWMKGGPyuC5nl/SipF2S/iGpg+tmHSQtcVcGAAA83YEDB9ShQwc98cQTWr16tSIiIrRv3z69++67ypo1q9Px4GbunCNWSNLHxhg/XS58C6y1S40xP0paYIx5Q9JBSS3cmAEAAI909OhRRUVFaeLEifLz81P37t3Vp08fPfDAA05HQwZyWxGz1m6RVPE6249LesFdzwsAgCc7ffq0hg8frhEjRigxMVGdOnVSSEiIHn74YaejwQFM1gcAIAMkJiZq/PjxGjRokI4fP64WLVooIiJCpUqVcjoaHMSHfgMA4EbJycmaPn26SpYsqR49eqhy5cqKiYnR/PnzKWGgiAEA4A7WWi1cuFDly5dX586dVbhwYa1atUorVqxQ5cqVnY4HD0ERAwAgnX3zzTeqUqWKmjdvrixZsujLL7/U2rVrVbt2baejwcNQxAAASCfr169XnTp19OKLL+ro0aOaOXOmtm7dqqCgIIax4rooYgAA3KWdO3eqWbNmqlKlirZs2aJRo0Zp9+7d6tixI8NYcVOcNQkAwB365ZdfNGDAAM2aNUs5c+bUgAED9MEHHyhXrlxOR0MmQREDAOA2HTt2TIMGDdL48eNljFG3bt0UHBys/PnzOx0NmQxFDACANDp79qxGjBihYcOG6cKFC+rYsaNCQ0NVtGhRp6Mhk6KIAQBwC4mJiZo0aZKioqL0xx9/qFmzZoqIiFCZMmWcjoZMjsX6AADcQEpKimbNmqXSpUvrgw8+UIUKFbRu3TotXLiQEoZ0QREDAOAvrLX68ssvVaFCBb3++ut68MEHtXLlSq1atUpVqlRxOh68CEUMAIBrrF69WtWqVVPTpk2VkpKihQsX6qefflKdOnWcjgYvRBEDAEDShg0bVLduXdWuXVvx8fGaNm2atm3bpmbNmjGMFW5DEQMA+LS4uDi1aNFCgYGB2rBhg4YPH649e/bojTfe0D33cE4b3Iu/YQAAn/Tbb78pPDxcM2bMULZs2dS/f391795d999/v9PR4EMoYgAAn3L8+HENHjxYY8eOVWpqqrp27ap+/frpwQcfdDoafBBFDADgE86dO6fRo0dr6NChOnv2rNq3b6+wsDAVK1bM6WjwYRQxAIBXu3TpkqZMmaKIiAgdPXpUjRs3VmRkpMqXL+90NIAiBgDwTikpKfrss88UEhKiAwcOqFatWlq8eLGqVavmdDTgKs6aBAB4FWutvvrqKwUEBKh9+/bKmzevli9ffnU+GOBJKGIAAK/x7bffqnr16nrllVd08eJFzZ8/XzExMapbty6zwOCRKGIAgExv06ZNql+/vmrVqqWDBw9qypQp2r59u1q0aKEsWfhVB8/F304AQKa1d+9etW7dWhUrVtTatWs1dOhQ7d27V2+++ab8/f2djgfcEov1AQCZzqFDhxQREaFp06Ypa9as+uijj9SzZ0/lyZPH6WjAbaGIAQAyjZMnT2ro0KEaPXq0kpOT9dZbb6lfv3566KGHnI4G3BGKGADA4124cEFjxozRkCFDdPr0abVp00YDBgxQ8eLFnY4G3BWKGADAYyUlJWnatGkKDw/X77//roYNGyoqKkoVKlRwOhqQLihiAACPk5qaqnnz5ikkJET79u3Ts88+q4ULF6p69epORwPSFWdNAgA8hrVWy5YtU6VKldSmTRvlzJlT//znP6/OBwO8DUUMAOAR/vvf/6pWrVpq0KCBzp49qzlz5ig2Nlb169dnGCu8FkUMAOCoLVu2qFGjRnr22We1Z88eTZgwQTt37lTr1q0Zxgqvx99wAIAj9u/fr7Zt2yogIEDff/+9Bg4cqL179+r//u//lDVrVqfjARmCxfoAgAz1+++/KyIiQlOmTJG/v7969+6tXr16KW/evE5HAzIcRQwAkCFOnTql6OhojRo1SpcuXVLnzp3Vv39/FS5c2OlogGMoYgAAt0pISNC4ceM0aNAgnTx5Uq1atVJ4eLgee+wxp6MBjmONGADALZKTkzV16lSVLFlSvXr1UtWqVRUbG6vPPvuMEga4UMQAAOkqNTVVCxYsUNmyZdWlSxcVLVpUa9as0bJlyxQQEOB0PMCjUMQAAOnCWqsVK1boqaee0muvvaZ7771XS5YsuTofDMD/oogBAO7a2rVrVbt2bdWrV08nTpzQ7NmztWnTJr3yyisMYwVugiIGALhj27dvV1BQkKpVq6YdO3Zo7NixiouLU7t27eTn5+d0PMDjUcQAALftwIED6tChg5544gmtXr1akZGR2rdvn959912GsQK3gfEVAIA0O3r0qKKiojRx4kRlyZJF3bt3V58+ffTAAw84HQ3IlChiAIBbOnPmjIYNG6YRI0YoMTFRr7/+ukJDQ/Xwww87HQ3I1ChiAIAbSkxM1IQJEzRw4EAdP35czZs3V0REhEqXLu10NMArsEYMAPA/kpOTNX36dJUsWVLdu3dX5cqVFRMTowULFlDCgHREEQMAXGWt1RdffKHy5curc+fOKly4sFatWqUVK1aocuXKTscDvA5vTQIAJEnffPONgoODFRMTozJlyujLL79U48aNmQOGu7Y4Nl7RK+J06FSCCufJrp51SyuoYhGnY3kEt70iZox5xBiz2hizwxiz3Rjzvmt7mDEm3hizyfVV310ZAAC3tn79etWpU0cvvviijh49qpkzZ2rr1q0KCgqihOGuLY6NV/CirYo/lSArKf5UgoIXbdXi2Hino3kEd741mSypu7W2rKSqkt4xxpR1/WyktTbA9bXMjRkAADewa9cuvfrqq6pSpYo2b96sUaNGaffu3erYsSPDWJFuolfEKSEp5U/bEpJSFL0izqFEnsVtb01aaw9LOuy6fNYYs1MSr0MCgMN++eUXDRgwQLNmzVKOHDk0YMAAffDBB8qVK5fT0eCFDp1KuK3tviZDFusbY4pJqihpnWvT340xW4wxM4wxeTMiAwD4uj/++EMffvihSpUqpU8//VTvvfee9u/fr5CQEEoY3KZwnuy3td3XuL2IGWPuk/SFpG7W2jOSJkoqLilAl18xG36D+3UxxsQYY2KOHTvm7pgA4LXOnj2r8PBwFS9eXKNHj1br1q21Z88ejRw5UgUKFHA6Hrxcz7qlld3/z291Z/f3U8+6jEGR3HzWpDHGX5dL2Bxr7SJJstYeuebnUyUtvd59rbVTJE2RpMDAQOvOnADgjS5evKhJkyYpKipKx44dU9OmTRUZGakyZco4HQ0+5MrZkZw1eX1uK2Lm8qk20yXttNaOuGZ7Idf6MUlqImmbuzIAgC9KSUnRJ598otDQUP3yyy+qXbu2Bg0apCpVqjgdDT4qqGIRitcNuPMVseqS2knaaozZ5Nr2kaRWxpgASVbSAUlvuTEDAPgMa60WL16sfv36aceOHapcubKmT5+uOnXqOB0NwA2486zJ7yVdbwAN4yoAIJ2tXr1affr00U8//aTSpUtr4cKFatq0KXPAAA/HRxwBQCa2YcMG1a1bV7Vr19ahQ4c0bdo0bdu2Tc2aNaOEAZkARQwAMqHdu3frtddeU2BgoGJiYjRs2DDt3r1bb7zxhu65h0+vAzILjlYAyER+++03hYeHa8aMGcqWLZv69++v7t276/7773c6GoA7QBEDgEzg+PHjGjx4sMaNG6eUlBS98847+uijj1SwYEGnowG4CxQxAPBg586d0+jRozV06FCdPXtW7du3V1hYmIoVK+Z0NADpgCIGAB7o0qVLmjJliiIjI3XkyBE1btxYkZGRKl++vNPRAKQjihgAeJCUlBTNnTtXISEh+vnnn1WrVi0tXrxYVatWdToaADfgrEkA8ADWWi1dulQVK1ZUu3btlCdPHi1fvlyrV6+mhAFejCIGAA777rvv9Oyzz6pRo0ZKTEzUvHnzFBMTo7p16zILDPByFDEAcMjmzZvVoEED1axZUwcOHNDkyZO1fft2vfbaa8qShX+eAV/AkQ4AGWzv3r1q3bq1AgIC9OOPP2ro0KHau3evunTpIn9/f6fjAchALNYHgAxy6NAhRUREaNq0acqaNav69u2rHj16KE+ePE5HA+AQihgAuNnJkyc1ZMgQjRkzRklJSXrrrbfUr18/PfTQQ05HA+AwihgAuMmFCxc0ZswYDRkyRKdPn1arVq0UHh6uEiVKOB0NgIe45RoxY8wTGREEALxFUlKSJk6cqMcee0zBwcF69tlntWnTJs2ZM4cSBuBP0rJYf4Ix5idjTFdjDJ8qCwA3kJqaqrlz56pMmTLq2rWrSpQooe+++05fffWVKlSo4HQ8AB7olkXMWltDUhtJj0jaYIz5zBjzotuTAUAmYa3VsmXLVKlSJbVu3Vo5c+bU0qVL9e233+rZZ591Oh4AD5am8RXW2j2S+knqLamWpDHGmF3GmKbuDAcAnu6///2vatWqpQYNGujs2bOaM2eOYmNj1aBBA4axAriltKwRq2CMGSlpp6TakhpZa8u4Lo90cz4A8Ehbt27VK6+8omeffVZ79uzR+PHjtXPnTrVu3ZphrADSLC1nTY6VNE3SR9bahCsbrbWHjDH93JYMADzQ/v37FRoaqjlz5ih37twaOHCg3nvvPeXMmdPpaAAyoVsWMWttrZv87JP0jQMAnunIkSOKjIzU5MmT5efnp169eqlXr17Kly+f09EAZGLMEQOAmzh9+rSio6M1atQoJSYmqnPnzgoJCVHhwoWdjgbAC1DEAOA6EhISNG7cOA0ePFgnTpxQy5YtFR4erpIlSzodDYAXYUUpAFwjOTlZU6dOVcmSJdWrVy89/fTT2rhxo+bOnUsJA5DubviKmDHmK0n2Rj+31r7ilkQA4IDU1FQtXLhQ/fr10549e1StWjXNmTNHtWrdcJksANy1m701Ocz1vamkhyR96rreStIRd4YC4DkWx8YrekWcDp1KUOE82dWzbmkFVSzidKx0Y63V119/reDgYMXGxqp8+fJasmSJGjVqxBwwAG53wyJmrf2PJBljhltrA6/50VfGmBi3JwPguMWx8QpetFUJSSmSpPhTCQpetFWSvKKMrV27VsHBwVqzZo2KFSum2bNnq3Xr1vLz83M6GgAfkZY1YjmNMcWvXDHG/E0SA3MAHxC9Iu5qCbsiISlF0SviHEqUPrZv366goCBVq1ZNO3bs0JgxY7Rr1y61a9eOEgYgQ6XlrMkPJK0xxuyXZCQ9Kuktt6YC4BEOnUq4re2e7sCBAwoLC9Ps2bOVK1cuRUREqFu3brrvvvucjgbAR6VloOtyY0xJSY+7Nu2y1l50bywAnqBwnuyKv07pKpwnuwNp7tzRo0c1cOBATZw4UcYYde/eXX369NEDDzzgdDQAPi4tnzWZQ1JPSe9aazdLKmqMaej2ZAAc17NuaWX3//Nbddn9/dSzbmmHEt2eM2fOKDQ0VCVKlNDYsWPVvn177d27V9HR0ZQwAB4hLW9NzpS0QVI11/V4SZ9LWuquUAA8w5UF+ZntrMnExERNmDBBAwcO1PHjx/Xqq68qIiJCjz/++K3vDAAZKC1FrIS19jVjTCtJstZeMJzTDfiMoIpFPL54XZGcnKzZs2crLCxMv/76q+rUqaNBgwYpMDDw1ncGAAek5azJS8aY7HINdzXGlJDEGjEAHsNaq0WLFumJJ57QG2+8oYceekjffPONVq5cSQkD4NHSUsTCJC2X9IgxZo6kVZJ6uzMUAKTVqlWr9PTTT6tZs2YyxmjRokVat26dXnjhBaejAcAtpeWsya+NMRskVdXl8RXvW2v/cHsyALiJ9evXKzg4WKtWrVLRokU1Y8YMtWvXTvfck5YVFwDgGdJy1uQqa+1xa+0/rbVLrbV/GGNWZUQ4APirXbt26dVXX1WVKlW0efNmjRo1Srt379brr79OCQOQ6dzsQ7+zScohKb8xJq8uvxomSbklZY6VuwC8xi+//KIBAwZo1qxZypEjh8LCwvThhx8qV65cTkcDgDt2s/99fEtSN0mFdXl8xZUidkbSODfnAgBJ0rFjxzRo0CBNmDBB1lq9//77Cg4OVoECBZyOBgB37WYf+j1a0mhjzN+ttWMzMBMA6OzZsxo5cqSGDRum8+fPq0OHDgoLC1PRokWdjgYA6SYtZ02mGmPyXLlijMlrjOnqxkwAfNjFixc1ZswYlShRQqGhoapTp462bdumGTNmUMIAeJ20FLE3rbWnrlyx1p6U9Kb7IgHwRSkpKfr4449VunRpvf/++ypfvrzWrl2rRYsWqUyZMk7HAwC3SEsR87t2kr4xxk9SVvdFAuBLrLVavHixnnzySXXs2FH58+fX119/fXU+GAB4s7QUseWS5htjXjDGvCBprmsbANyVNWvWqFq1amrSpImSkpK0YMECrV+/Xi+++KL4JDUAviAtQ3d66/IZlP/nur5S0jS3JQLg9TZu3KiPPvpIK1asUJEiRTRt2jR16NCBOWAAfE5aJuunSpro+gKAO7Z79271799fCxYsUL58+TRs2DB17dpV2bNndzoaADjiZgNdF1hrWxhjtsr1gd/XstZWcGsyAF4jPj5e4eHhmj59urJly6Z+/fqpR48euv/++52OBgCOutkrYu+7vjfMiCAAvM+JEyc0ePBgjR07VikpKeratav69u2rggULOh0NADzCzQa6HnZ9P5hxcQB4g/Pnz2vUqFGKjo7WmTNn1K5dO4WFhelvf/ub09HuyuLYeEWviNOhUwkqnCe7etYtraCKfOIbgDt3s7cmz+o6b0leYa3NfbMHNsY8Imm2pIKux5lirR1tjMknab6kYpIOSGrhmk0GIJO7dOmSpk6dqoiICB05ckSNGzdWZGSkypcv73S0u7Y4Nl7Bi7YqISlFkhR/KkHBi7ZKEmUMwB274fgKa20uV9kaLamPLn/Q98O6fBblqDQ8drKk7tbaspKqSnrHGFPW9VirrLUlJa1yXQeQiaWkpOjTTz/V448/rnfffVePP/64fvjhBy1evNgrSpgkRa+Iu1rCrkhISlH0ijiHEgHwBmmZI/aKtXaCtfastfaMtXaipMa3upO19rC1dqPr8llJO3W5zDWW9LHrZh9LCrqz6ACcZq3V0qVLVbFiRbVr10558uTR8uXLtXr1alWrVs3peOnq0KmE29oOAGmRliJ23hjTxhjjZ4zJYoxpI+n87TyJMaaYpIqS1kkqeGX9maTfdfmtSwCZzHfffacaNWqoUaNGSkhI0Lx58xQTE6O6det65TDWwnmuP2LjRtsBIC3SUsRaS2oh6Yjrq7lrW5oYY+6T9IWkbtbaM9f+zFprdYN1aMaYLsaYGGNMzLFjx9L6dADcbPPmzWrQoIFq1qyp/fv3a9KkSdqxY4dee+01ZcmSln9SMqeedUsru7/fn7Zl9/dTz7qlHUoEwBukZaDrAaXhrcjrMcb463IJm2OtXeTafMQYU8hae9gYU0jS0Rs87xRJUyQpMDDwhicNAMgY+/btU0hIiD777DPlyZNHgwcP1t///nflyJHD6WgZ4sqCfM6aBJCeblnEjDGldHmqfkFrbXljTAVdXjcWeYv7GUnTJe201o645kf/kNRB0mDX9yV3Gh6A+x0+fFgRERGaOnWq/P39FRwcrJ49eypv3rxOR8twQRWLULwApKu0vI8wVVKwpCRJstZukdQyDferLqmdpNrGmE2ur/q6XMBeNMbskVTHdR2Ahzl58qSCg4NVokQJTZ06VV26dNG+ffs0cOBAnyx2mCamAAAZqUlEQVRhAOAOafmE3RzW2p/+svg2+VZ3stZ+L+lGK3ZfSMPzAnDAhQsXNGbMGA0ZMkSnT59W69atNWDAAJUoUcLpaADgddJSxP4wxpSQa1G9MeZVSYdvfhcAmU1SUpKmT5+u8PBwHT58WA0aNFBUVJSefPJJp6MBgNdKSxF7R5cXzT9ujImX9LOkNm5NBSDDpKamav78+erfv7/27dun6tWra/78+apRo4bT0QDA6920iBljskgKtNbWMcbklJTFNZwVQCZnrdXy5csVHByszZs3q0KFClq6dKnq16/vlXPAAMAT3XSxvrU2VVIv1+XzlDDAO/zwww967rnnVL9+fZ09e1Zz5sxRbGysGjRoQAkDgAyUlrMmvzHG9DDGPGKMyXfly+3JAKS7rVu36pVXXlH16tW1e/dujR8/Xjt37lTr1q29ehgrAHiqtKwRe831/Z1rtllJxdM/DgB3+PnnnxUaGqpPP/1UuXPn1sCBA/Xee+8pZ86cTkcDAJ+Wlsn6f8uIIADS35EjRxQZGanJkyfLz89PvXr1Uq9evZQvHy9qA4AnSMtk/WySukp6VpdfCftO0iRrbaKbswG4Q6dPn1Z0dLRGjRqlxMREde7cWf3791eRIkyFBwBPkpa3JmdLOitprOt6a0mf6PKHfwPwIAkJCRo/frwGDRqkEydOqGXLlgoPD1fJkiX/dLvFsfF8ZiIAeIC0FLHy1tqy11xfbYzZ4a5AAG5fcnKyZs6cqQEDBig+Pl5169bVoEGDVLFixf+57eLYeAUv2qqEpBRJUvypBAUv2ipJlDEAyGBpOU1qozGm6pUrxpinJcW4LxKAtEpNTdXnn3+ucuXKqUuXLipatKjWrFmj5cuXX7eESVL0irirJeyKhKQURa+Iy4jIAIBrpOUVscqSfjDG/OK6XlRSnDFmqyRrra3gtnQArstaq5UrVyo4OFgbN25UuXLltGTJEjVq1OiWc8AOnUq4re0AAPdJSxGr5/YUANJs3bp1Cg4O1urVq1WsWDHNnj1brVu3lp+fX5ruXzhPdsVfp3QVzpM9vaMCAG7hlm9NWmsP3uwrI0ICkLZv364mTZqoatWq2r59u8aOHatdu3apXbt2aS5hktSzbmll9//z7bP7+6ln3dLpHRkAcAtpeUUMgIMOHjyo0NBQffLJJ7rvvvsUERGhbt266b777rujx7uyIJ+zJgHAeRQxwEMdPXpUAwcO1MSJE2WMUbdu3RQcHKz8+fPf9WMHVSxC8QIAD0ARg0dhvpV05swZjRgxQsOHD9eFCxfUqVMnhYSE6JFHHnE6GgAgnVHE4DF8fb5VYmKiJk6cqKioKB0/flyvvvqqIiIi9PjjjzsdDQDgJmmZIwZkCF+db5WcnKwZM2aoVKlS+vDDD1WpUiWtX79en3/+OSUMALwcRQwew9fmW1lrtWjRIj3xxBN64403VKhQIa1atUpff/21AgMDnY4HAMgAFDF4jBvNsfLG+VarVq3S008/rWbNmskYo0WLFmnt2rWqXbu209EAABmIIgaP4QvzrdavX68XX3xRderU0e+//66ZM2dq69atatKkyS0n4gMAvA9FDB4jqGIRDWr6hIrkyS4jqUie7BrU9AmvWKi/a9cuNW/eXFWqVNGmTZs0cuRI7d69Wx07drytYawAAO/CWZPwKN423+rXX3/VgAEDNHPmTOXIkUOhoaH68MMPlTt3bqejAQA8AEUMcIM//vhDgwcP1rhx42St1XvvvaePPvpIBQoUcDoaAMCDUMSAdHTu3DmNHDlS0dHROn/+vDp06KDQ0FA9+uijTkcDAHggihiQDi5evKjJkycrMjJSx44dU5MmTRQZGamyZcs6HQ0A4MFYrA/chZSUFH388ccqXbq03n//fZUvX15r167VokWLKGEAgFuiiAF3wFqrJUuWqEKFCurYsaMeeOABff3111fngwEAkBYUMeA2rVmzRs8884yCgoKUnJysBQsWXJ0PxiwwAMDtoIgBabRx40bVq1dPzz//vH799VdNnTpV27dvV/PmzZUlC4cSAOD28dsDuIXdu3frtddeU+XKlbV+/XpFR0drz5496ty5s+65h/NdAAB3jt8iwA3Ex8crPDxc06dPV7Zs2dSvXz/16NFD999/v9PRAABegiIG/MWJEyc0ePBgjR07VikpKeratav69u2rggULOh0t3SyOjVf0ijgdOpWgwnmyq2fd0l71iQYAkFlQxACX8+fPa/To0Ro6dKjOnDmjtm3bKjw8XMWKFXM6WrpaHBuv4EVblZCUIkmKP5Wg4EVbJYkyBgAZjDVi8HmXLl3ShAkTVKJECfXt21e1atXSli1bNHv2bK8rYZIUvSLuagm7IiEpRdEr4hxKBAC+i1fE4LNSU1M1d+5chYSEaP/+/apZs6a+/PJLVatWzelobnXoVMJtbQcAuA+viMHnWGu1dOlSBQQEqG3btsqdO7eWLVumNWvWeH0Jk6TCebLf1nYAgPtQxOBTvvvuO9WoUUONGjVSQkKC5s6dqw0bNujll1/2mWGsPeuWVnZ/vz9ty+7vp551SzuUCAB8F0UMPmHz5s1q0KCBatasqf3792vSpEnasWOHWrZs6XPDWIMqFtGgpk+oSJ7sMpKK5MmuQU2fYKE+ADiANWLwanv37lVISIjmzp2rPHnyaMiQIXr33XeVI0cOp6M5KqhiEYoXAHgAihi80uHDhxUREaGpU6fK399fwcHB6tmzp/Lmzet0NAAArqKIwaucOnVKQ4cO1ahRo5SUlKQ333xT/fv3V6FChZyOBgDA/6CIwStcuHBB48aN0+DBg3Xq1Cm1atVK4eHhKlGihNPRAAC4Id9apQyvk5SUpMmTJ+uxxx5T79699cwzzyg2NlZz5syhhAEAPB5FDJlSamqq5s2bp7Jly+rtt99W8eLF9e2332rp0qV68sknnY4HAECaUMSQqVhr9a9//UuVK1dWq1atlD17di1duvTqfDAAADITihgyjR9++EHPPfec6tevr9OnT+vTTz/Vpk2b1KBBA58ZxgoA8C5uK2LGmBnGmKPGmG3XbAszxsQbYza5vuq76/nhPbZt26bGjRurevXqiouL07hx47Rr1y61adPG54axAgC8izt/i82SVO8620daawNcX8vc+PzI5H7++We1b99eFSpU0Jo1axQVFaV9+/bpnXfeUdasWZ2OBwDAXXPb+Apr7bfGmGLuenx4ryNHjigyMlKTJ0+Wn5+fevTooT59+ihfvnxORwMAIF05MUfs78aY9pJiJHW31p50IAM80OnTpxUdHa1Ro0YpMTFRnTt3Vv/+/VWkCB/FAwDwThm9wGaipOKSAiQdljT8Rjc0xnQxxsQYY2KOHTuWUfnggISEBA0bNkzFixdXVFSUGjZsqJ07d2rSpEmUMACAV8vQImatPWKtTbHWpkqaKqnKTW47xVobaK0NLFCgQMaFRIZJTk7W1KlTVbJkSfXs2VNVqlTRhg0bNG/ePJUsWdLpeAAAuF2GFjFjzLUf+NdE0rYb3RbeKzU1VZ9//rnKlSunLl26qGjRolqzZo3+9a9/qVKlSk7HAwAgw7htjZgxZq6k5yTlN8b8JilU0nPGmABJVtIBSW+56/nheay1WrlypYKDg7Vx40aVK1dOS5YsUaNGjZgDBgDwSe48a7LVdTZPd9fzwbOtW7dOwcHBWr16tR599FHNmjVLbdu2lZ+fn9PRAABwDNMw4VY7duxQ06ZNVbVqVW3btk1jxoxRXFycOnToQAkDAPg8J8ZXwAccPHhQYWFhmj17tnLmzKnw8HB169ZNuXLlcjoaAAAegyKGdHX06FENHDhQEydOlDFGH3zwgfr06aP8+fM7HQ0AAI9DEUO6OHPmjEaMGKHhw4frwoULev311xUaGqpHHnnE6WgAAHgsihjuSmJioiZOnKioqCgdP35czZo1U2RkpB5//HGnowEA4PFYrI87kpycrJkzZ6pUqVL68MMPVbFiRf30009auHAhJQwAgDSiiOG2WGu1aNEiVahQQZ06ddJDDz2kb775RitXrtRTTz3ldDwAADIVihjS7N///reqVq2qZs2ayVqrhQsXat26dXrhhRecjgYAQKZEEcMtxcTE6KWXXtILL7ygw4cPa8aMGdq6dauaNWvGRHwAAO4CRQw3tGvXLr366qt66qmnFBsbqxEjRmj37t16/fXXdc89nOcBAMDd4rcp/sevv/6qAQMGaObMmcqRI4dCQ0P14YcfKnfu3E5HAwDAq1DEcNXx48c1aNAgjRs3TtZavffee/roo49UoEABp6MBAOCVKGLQuXPnNHLkSA0bNkznzp1T+/btFRYWpkcffdTpaAAAeDWKmA+7ePGiJk+erKioKB09elRNmjRRZGSkypYt63Q0AAB8AkXMB6WkpGjOnDkKCQnRwYMH9fzzz+sf//iHnn76aaejAQDgUzhr0odYa7VkyRI9+eST6tChg/Lnz68VK1Zo1apVlDAAABxAEfMRa9as0TPPPKOgoCAlJSXp888/1/r16/XSSy8xCwwAAIdQxLxcbGys6tWrp+eff16//vqrpkyZou3bt+vVV1+lgAEA4DCKmJfas2ePWrZsqUqVKmn9+vUaNmyY9uzZozfffJNhrAAAeAh+I3uZ+Ph4hYeHa/r06br33nvVt29f9ezZU/fff7/T0QAAwF9QxLzEiRMnNGTIEI0ZM0YpKSnq2rWr+vbtq4IFCzodDQAA3ABFLJM7f/68Ro8eraFDh+rMmTNq27atBgwYoL/97W9ORwMAALdAEcukLl26pGnTpikiIkK///67GjVqpKioKD3xxBNORwMAAGlEEctkUlNTNXfuXIWEhGj//v2qWbOmvvjiCz3zzDNORwMAALeJsyYzCWut/vnPf6pixYpq27atcufOrX/9619X54MBAIDMhyKWCXz//feqWbOmGjZsqAsXLmju3LnasGGD6tWrxywwAAAyMYqYB9u8ebMaNmyoGjVqaN++fZo0aZJ27Nihli1bKksWdh0AAJkdv8090L59+9SmTRtVrFhR//3vfzV48GDt3btXb731lvz9/Z2OBwAA0gmL9T3I4cOHFRERoalTp8rf3199+vRRz549lTdvXqejAQAAN6CIeYBTp05p6NChGjVqlJKSkvTmm2+qf//+KlSokNPRAACAG1HEHHThwgWNGzdOgwcP1smTJ9W6dWuFh4erRIkSTkcDAAAZgDViDkhKStLkyZP12GOPqXfv3qpatapiY2M1Z84cShgAAD6EIpaBUlNTNW/ePJUtW1Zvv/22ihcvrm+//VbLli1TQECA0/EAAEAGo4hlAGutli9frsDAQLVq1UrZs2fXV199pe+++041atRwOh4AAHAIRczNfvjhBz333HN6+eWXderUKX3yySeKjY1Vw4YNGcYKAICPo4i5ybZt29S4cWNVr15dcXFxGjdunHbt2qW2bdvKz8/P6XgAAMADUMTS2c8//6z27durQoUKWrNmjaKiorRv3z698847ypo1q9PxAACAB2F8RTo5cuSIoqKiNGnSJPn5+alnz57q3bu38uXL53Q0AADgoShid+n06dMaNmyYRo4cqcTERL3xxhsKCQlRkSJFnI4GAAA8HEXsDiUkJGj8+PEaNGiQTpw4oddee03h4eEqVaqU09EAAEAmwRqx25ScnKxp06apZMmS6tmzp5566ilt2LBB8+bNo4QBAIDbQhFLo9TUVH3++ecqV66c3nzzTT3yyCNavXq1li9frkqVKjkdDwAAZEIUsVuw1urrr79WlSpV1KJFC91zzz1avHjx1flgAAAAd4oidhPr1q3TCy+8oLp16+qPP/7QrFmztGXLFjVu3JhhrAAA4K5RxK5jx44datq0qapWrapt27Zp9OjRiouLU4cOHRjGCgAA0g1nTV7j4MGDCgsL0+zZs5UzZ06Fh4erW7duypUrl9PRAACAF6KISTp27JgGDhyoCRMmyBijbt26KTg4WPnz53c6GgAA8GIUMUnTpk3TmDFj9Prrrys0NFSPPPKI05EAAIAPMNZa9zywMTMkNZR01Fpb3rUtn6T5kopJOiCphbX25K0eKzAw0MbExLglpySdP39ev/zyi8qUKeO25wAAAL7BGLPBWhuYltu6c7H+LEn1/rKtj6RV1tqSkla5rjsuZ86clDAAAJDh3FbErLXfSjrxl82NJX3suvyxpCB3PT8AAICny+jxFQWttYddl3+XVDCDnx8AAMBjODZHzF5enHbDBWrGmC7GmBhjTMyxY8cyMBkAAEDGyOgidsQYU0iSXN+P3uiG1top1tpAa21ggQIFMiwgAABARsnoIvYPSR1clztIWpLBzw8AAOAx3FbEjDFzJf0oqbQx5jdjzBuSBkt60RizR1Id13UAAACf5LaBrtbaVjf40Qvuek4AAIDMhA/9BgAAcAhFDAAAwCEUMQAAAIdQxAAAABxCEQMAAHAIRQwAAMAhFDEAAACHUMQAAAAcQhEDAABwCEUMAADAIRQxAAAAh1DEAAAAHEIRAwAAcAhFDAAAwCEUMQAAAIfc43QApy2OjVf0ijgdOpWgwnmyq2fd0gqqWMTpWAAAwAf4dBFbHBuv4EVblZCUIkmKP5Wg4EVbJYkyBgAA3M6n35qMXhF3tYRdkZCUougVcQ4lAgAAvsSni9ihUwm3tR0AACA9+XQRK5wn+21tBwAASE8+XcR61i2t7P5+f9qW3d9PPeuWdigRAADwJT69WP/KgnzOmgQAAE7w6SImXS5jFC8AAOAEn35rEgAAwEkUMQAAAIdQxAAAABxCEQMAAHAIRQwAAMAhFDEAAACHUMQAAAAcQhEDAABwCEUMAADAIRQxAAAAhxhrrdMZbskYc0zSwbt8mPyS/kiHOMgY7K/Mhf2VebCvMhf2V+ZyZX89aq0tkJY7ZIoilh6MMTHW2kCncyBt2F+ZC/sr82BfZS7sr8zlTvYXb00CAAA4hCIGAADgEF8qYlOcDoDbwv7KXNhfmQf7KnNhf2Uut72/fGaNGAAAgKfxpVfEAAAAPIrXFTFjzAxjzFFjzLZrtuUzxqw0xuxxfc/rZEb8fzfYX2HGmHhjzCbXV30nM+L/M8Y8YoxZbYzZYYzZbox537WdY8wD3WR/cYx5IGNMNmPMT8aYza79NcC1nePLA91kf93W8eV1b00aY2pKOidptrW2vGvbUEknrLWDjTF9JOW11vZ2Micuu8H+CpN0zlo7zMls+F/GmEKSCllrNxpjcknaIClIUkdxjHmcm+yvFuIY8zjGGCMpp7X2nDHGX9L3kt6X1FQcXx7nJvurnm7j+PK6V8Sstd9KOvGXzY0lfey6/LEu/0MED3CD/QUPZa09bK3d6Lp8VtJOSUXEMeaRbrK/4IHsZedcV/1dX1YcXx7pJvvrtnhdEbuBgtbaw67Lv0sq6GQYpMnfjTFbXG9d8jK8BzLGFJNUUdI6cYx5vL/sL4ljzCMZY/yMMZskHZW00lrL8eXBbrC/pNs4vnyliF1lL78X613vx3qfiZKKSwqQdFjScGfj4K+MMfdJ+kJSN2vtmWt/xjHmea6zvzjGPJS1NsVaGyDpYUlVjDHl//Jzji8PcoP9dVvHl68UsSOutRJX1kwcdTgPbsJae8T1lztV0lRJVZzOhP/PtRbiC0lzrLWLXJs5xjzU9fYXx5jns9aekrRal9cbcXx5uGv31+0eX75SxP4hqYPrcgdJSxzMglu48g+OSxNJ2250W2Qs1+LU6ZJ2WmtHXPMjjjEPdKP9xTHmmYwxBYwxeVyXs0t6UdIucXx5pBvtr9s9vrzxrMm5kp7T5U9APyIpVNJiSQskFZV0UFILay0LxD3ADfbXc7r8kq6VdEDSW9esj4CDjDHPSvpO0lZJqa7NH+nyuiOOMQ9zk/3VShxjHscYU0GXF+P76fILJQusteHGmAfE8eVxbrK/PtFtHF9eV8QAAAAyC195axIAAMDjUMQAAAAcQhEDAABwCEUMAADAIRQxAAAAh1DEAHgtY8xzxphn7vIxzt36VgBwZyhiALzZc5LuqogBgDtRxABkKsaYxcaYDcaY7caYLtdsr2eM2WiM2WyMWeX6kOu3JX1gjNlkjKlhjJlljHn1mvucc32/z3WfjcaYrcaYxrfIEG6M6XbN9ShjzPvp/WcF4P0Y6AogUzHG5LPWnnB9pMh6SbV0+X8qN0qqaa39+ZrbhEk6Z60d5rrvLElLrbULXdfPWWvvM8bcIymHtfaMMSa/pLWSSlpr7ZXb/CVDMUmLrLWVjDFZJO2RVMVaezwD/hMA8CL3OB0AAG7Te8aYJq7Lj0gqKamApG+ttT9L0h18/IuRNNAYU1OXPwqoiKSCkn6/3o2ttQeMMceNMRVdt4ulhAG4ExQxAJmGMeY5SXUkVbPWXjDGrJGU7TYeIlmuJRmuV7Kyura30eUyV9lam2SMOZCGx50mqaOkhyTNuI0MAHAVa8QAZCb3SzrpKmGPS6rq2r5WUk1jzN+ky29furaflZTrmvsfkFTZdfkVSf7XPO5RVwl7XtKjacjypaR6kp6StOLO/jgAfB1FDEBmslzSPcaYnZIG63IBk7X2mKQukhYZYzZLmu+6/VeSmlxZrC9pqqRarttUk3Tedbs5kgKNMVsltZe061ZBrLWXJK2WtMBam5Jef0AAvoXF+gBwB1xvbW6U1Nxau8fpPAAyJ14RA4DbZIwpK2mvpFWUMAB3g1fEAAAAHMIrYgAAAA6hiAEAADiEIgYAAOAQihgAAIBDKGIAAAAOoYgBAAA45P8BAx0mESGsBzwAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAFpCAYAAAC4SK2+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VfX9x/HXlxDhiiMOHEQR60iroCA4wYXaWCfiHhRHxclqSyvd/XVoSyuIKO6Ko27E3bgHKigYFBypuzY4cERFoyKc3x/nRoEGMm5yzx2v5+ORB8nJuTefHI769ny/3883RFGEJEmSWqdD0gVIkiTlM8OUJElSBgxTkiRJGTBMSZIkZcAwJUmSlAHDlCRJUgYMU5IkSRkwTEmSJGWgyTAVQugcQngqhPBsCOH5EMLv08d/F0KoDSHMSX/s1/7lSpIk5ZbQVAf0EEIAukRRtDCEUApMB0YC+wILoyj6W/uXKUmSlJs6NnVCFKethekvS9MfrdqDZt1114169OjRmpdKkiRl1ezZs9+PoqhrU+c1GaYAQgglwGxgc+CCKIpmhhB+AAwPIfwQmAX8JIqij1b2Pj169GDWrFnN+ZGSJEmJCiG82ZzzmjUBPYqixVEU9QY2AnYIIfQEJgPfAXoDbwN/X0Ehw0IIs0IIsxYsWNCs4iVJkvJFi1bzRVFUBzwE7BtF0bvpkLUEuBTYYQWvuSSKon5RFPXr2rXJJ2WSJEl5pTmr+bqGEMrSn6eAfYCXQggbLnXaIcC89ilRkiQpdzVnztSGwJT0vKkOwI1RFN0ZQrg6hNCbeDL6G8Ap7VemJElSbmrOar7ngD6NHB/SLhVJkiTlETugS5IkZcAwJUmSlAHDlCRJUgYMU5IkSRloVgd0tZ1p1bWMq6phfl093cpSjKmsYFCf8qTLkiRJrWSYyqJp1bWMnTqX+kWLAaitq2fs1LkABipJkvKUw3xZNK6q5psg1aB+0WLGVdUkVJEkScqUYSqL5tfVt+i4JEnKfYapLOpWlmrRcUmSlPsMU1k0prKCVGnJMsdSpSWMqaxIqCJJkpQpJ6BnUcMkc1fzSZJUOAxTWTaoT7nhSZKkAuIwnyRJUgYMU5IkSRkwTEmSJGXAMCVJkpQBw5QkSVIGDFOSJEkZMExJkiRlwDAlSZKUAcOUJElSBgxTkiRJGXA7GTXbtOpa9xWUJGk5hik1y7TqWsZOnUv9osUA1NbVM3bqXAADlSSpqDnMp2YZV1XzTZBqUL9oMeOqahKqSJKk3GCYUrPMr6tv0XFJkoqFYUrN0q0s1aLjkiQVC8OUmmVMZQWp0pJljqVKSxhTWZFQRZIk5QYnoKtZGiaZu5pPkqRlGabUbIP6lBueJElajsN8kiRJGTBMSZIkZcAwJUmSlAHDlCRJUgYMU5IkSRkwTEmSJGXAMCVJkpQBw5QkSVIGDFOSJEkZMExJkiRlwDAlSZKUAcOUJElSBgxTkiRJGTBMSZIkZcAwJUmSlIGOSRfQVqZV1zKuqob5dfV0K0sxprKCQX3K2+11kiRJUCBhalp1LWOnzqV+0WIAauvqGTt1LsBKg1FrXydJktSgIIb5xlXVfBOIGtQvWsy4qpp2eZ0kSVKDJsNUCKFzCOGpEMKzIYTnQwi/Tx9fO4RwXwjh5fSfa7V/uY2bX1ffouOZvk6SJKlBc55MfQkMjKJoW6A3sG8IYSfgLOCBKIq2AB5If52IbmWpFh3P9HWSJEkNmgxTUWxh+svS9EcEHAxMSR+fAgxqlwqbYUxlBanSkmWOpUpLGFNZ0S6vkyRJatCsCeghhBJgNrA5cEEURTNDCOtHUfR2+pR3gPXbqcYmNUwWb+mqvNa+TpIkqUGIoqj5J4dQBtwKDAemR1FUttT3Poqi6H/mTYUQhgHDALp37973zTffzLhoSZKk9hZCmB1FUb+mzmvRar4oiuqAh4B9gXdDCBumf9iGwHsreM0lURT1i6KoX9euXVvy4yRJknJec1bzdU0/kSKEkAL2AV4CbgeGpk8bCtzWXkVKkiTlqubMmdoQmJKeN9UBuDGKojtDCE8CN4YQTgLeBI5oxzolSZJyUpNhKoqi54A+jRz/ANirPYqSJEnKFwXRAV2SJCkphilJkqQMGKYkSZIyYJiSJEnKgGFKkiQpA4YpSZKkDBimJEmSMlBYYeo/M2HqKfDei0lXIkmSikRhhan3/w0v3g4X7gTXHQNvPZ10RZIkqcAVVpjabgiMmge7nwVvPg6X7w1XHgCv3A9RlHR1kiSpABVWmALosg7sORZGPw+Vf4YPXoVrDoVLdofnb4Uli5OuUJIkFZDCC1MNOq0GO58BI5+FgybBV5/BTcfDpO1h9hT4+sukK5QkSQWgcMNUg46rxMN/ZzwFh0+JQ9YdI+C8beGJ8+HLT5OuUJIk5bHCD1MNOpTA1oNg2CMw5FZYdwu491cwvic8+Cf47P2kK5QkSXmoeMJUgxBgs4Ew9A740QPQYwA8+tc4VN3zc6h7K+kKJUlSHim+MLW0jfrBUdfGQ4BbHwJPXwYTe8Otp8J7LyVdnSRJygPFHaYadK2AQybDiDmw/cnwwm1w4Y5w/bHw31lJVydJknKYYWppZRvDD85J96r6ObwxHS7bK92r6gF7VUmSpP9hmGpMl3Vgz1/A6Hnw/T/CB6/ANYPhkj3g+Wn2qpIkSd8wTK1Mp9Vhl+Fxr6oDJ8ZtFG4aaq8qSZL0DcNUc3TsBH2HwplPw+FXLterahJ8uTDpCiVJUkIMUy3RoSRe9dfQq2qdzeHeX8L4reGhP8NnHyRdoSRJyjLDVGs09Ko6/s5ve1U98heY0BPuOQs+/m/SFUqSpCwJURZXqPXr1y+aNatAWw0sqIHpE2DujfHX2xwJ/UdB1y0zfutp1bWMq6phfl093cpSjKmsYFCf8ozfV5IkrVgIYXYURf2aPM8w1cbq3oInJ6UnqH8B390fdv0xlPdt1dtNq65l7NS51C/6dgVhqrSEswf3MlBJktSOmhumHOZra2Ubww/+ErdV2O2n8MZjcOlAmHIQvPpQi3tVjauqWSZIAdQvWsy4qpq2rFqSJLWSYaq9dFkXBv4KRj8f96paUANXD4JL94w7rDezV9X8uvoWHZckSdllmGpvDb2qRj0X96r64mO48YdwwQ7wzNXw9VcrfXm3slSLjkuSpOwyTGXLN72qZsFh/4DSVeH2M+NeVU9esMJeVWMqK0iVlixzLFVawpjKimxULUmSmmCYyrYOJdBzMJzyKBx3C6yzGVT9Im6r8NDZ8PmHy5w+qE85Zw/uRXlZigCUl6WcfC5JUg5xNV8ueOtpmH4u1NwNpV2g7/Gw8xmwpoFJkqSkuJovn2y8PRx9HZw+A753IMy8KB7+m3YGvP9y0tVJkqSVMEzlkvW+B4MvhpFzoN8JMO/meFPlG4ZA7TNJVydJkhphmMpFZd1hv3Ewah7s+hN4/ZG4pcKUg+C1h1vcq0qSJLUfw1QuW60r7PXrOFTt84e4V9VVB8dNQF+4HZYsSbpCSZKKnmEqH3ReA/qPgJHPwgEToP4juHFI3Kuq+pome1VJkqT2Y5jKJ6Wd47lUw2fDYVdAx85w2xkwsTc8eSF89VnSFUqSVHQMU/moQwn0PBROfQyOvRnW6gFVY2F8T3j4L//Tq0qSJLUf+0wViv/MhOnj4d/3xL2q+p0Q96pao1uL3mZadS3jqmqYX1dPt7IUYyorbBAqSSpKze0zZZgqNO++AI9PgLk3Q+gA2x4F/UfBups3+dJp1bWMnTqX+kXfbsKcKi2x47okqSjZtLNYrb8VDL4ERlTHndTn3gST+sWbK8+vXulLx1XVLBOkAOoXLWZcVU07FixJUn4zTBWqtTaB/f+W7lX1Y3j1YbhkD7hqELz2SKO9qubX1Tf6Vis6LkmSDFOFb7WusNdvYPQ82Pv38O7zcNVBcNle8OIdy/Sq6laWavQtVnRckiQZpopH5zVgwCgYNRf2Pxc+/wBuOA4u3BGqr4XFixhTWUGqtGSZl6VKSxhTWZFQ0ZIk5T7DVLEp7QzbnwRnzoZDL4eSVeC20+G83gz68g7+etBmlJelCEB5WcrJ55IkNcHVfMUuiuDl++K2Cv95AlJrw06nwfY/glXXTro6SZIS42o+NU8IsOX34cR74MQq2HgHeOhPMKEXVP0SPpmfdIWSJOU0w5S+1X0nOOYGOO0JqNgPZkyG87aF24fD+68kXZ0kSTmpyTAVQtg4hPBQCOGFEMLzIYSR6eO/CyHUhhDmpD/2a/9ylRXrbw2HXgojnoHtfgjP3rBUr6o5SVcnSVJOaXLOVAhhQ2DDKIqeCSGsDswGBgFHAAujKPpbc3+Yc6by1ML34qdUT18GX34Cmw2EAaOhx67xMKEkSQWozeZMRVH0dhRFz6Q//xR4EXB5VzFZbT3Y+7fpXlW/g3fmwZQD4bK94cU7l+lVJUlSsWnRnKkQQg+gDzAzfWh4COG5EMIVIYS12rg25ZrOa8ZPpEY9F/eq+mwB3HAsTN4Z5lwHixclXaEkSVnX7DAVQlgNuAUYFUXRJ8Bk4DtAb+Bt4O8reN2wEMKsEMKsBQsWtEHJSlxpKu5VNfyZuFdVh44w7VSY2AdmXgxffZ50hZIkZU2z+kyFEEqBO4GqKIrObeT7PYA7oyjqubL3cc5UgfqmV9W58J8nYdV1YMfTYIcfQcoHlpKk/NRmc6ZCCAG4HHhx6SCVnpje4BBgXmsKVQH4plfVv+CEf0F5X3jojzC+J9z7K/jk7aQrlCSp3TRnNd8A4DFgLtAw0/gXwNHEQ3wR8AZwShRFK/2vpk+misg78+DxCTDvlngYcNujof9IWGezpCuTJKlZmvtkyu1k1L4+fB2emBhvprxkEWx1cDyJfcNtk65MrTCtupZxVTXMr6unW1mKMZUV7t0oqWAZppRbPn0XZlwIT18OX30Km+2V7lU1wF5VeWJadS1jp86lftHib46lSkvcDFtSwXJvPuWW1deHfX4f96ra6zfwznMw5QC4fB946W57VeWBcVU1ywQpgPpFixlXVZNQRZKUGwxTyq5UGez6Exg1F/b7Gyx8F64/GibvAs9eb6+qHDa/rr5FxyWpWBimlIzSFOxwMgyvhsGXxkN9t54CE7eDmZfYqyoHdStLtei4JBULw5SSVdIRtjkCTnsCjr4B1tgQ7hkDE3rBo+Og/qOkK1TamMoKUqUlyxxLlZYwprIioYokKTcYppQbQoCKfeGke+GEe6B8O3jwjzC+F9z7a/j0naQrLHqD+pRz9uBelJelCEB5WcrJ55KEq/mUy96ZC9PHw/O3xr2qeh8Du4ywV5UkKStczaf8t0EvOOwKGD4beh8Lc/4Jk/rBTSfA288lXZ0kSYBhSvlg7e/AgRPiFYC7DI/3Abx4V7jmUHjj8XhvQEmSEmKYUv5YfQPY5//iXlUDfw3z58CV+8EVlVBzj72qJEmJMEwp/6TKYLefxqFqv7/FGylfdxRc1B+evQEWf510hZKkImKYUv5q6FU14hk45JJ4uO/WYXB+H3jqUlhkM0lJUvszTCn/lZTCtkeme1VdD6ttAHf/FMb3hEf/BvV1SVcoSSpghikVjg4doOIHca+q4++Gbr3hwT/Eoeq+38SbLUuS1MYMUyo8IUCP/nDcLXDKY7DFPvDE+XFX9TtGwYevJV2hJKmAGKZU2DbcBg7/B5w5C3ofDXOuhfP7ws0nxU1BJUnKkGFKxWGdzeDA8+JeVTufCf/+F1w0AK49HN58IunqJEl5zO1klDXTqmsZV1XD/Lp6upWlGFNZkdy+bvUfwdOXwYzJ8PkHsPFOMGA0bFkZDxMqt/6+JCkBzd1OxjClrJhWXcvYqXOpX7T4m2Op0pLkN8r96nOovgaemAgfvwXrbQ0DRsHWg6GkY3J1JSxn/74kKYvcm085ZVxVzTL/YQaoX7SYcVU1CVWUtsqqsOMwGFENh1wM0WKYenLR96rK2b8vScpBhillxfy6xkPJio5nXUkpbHsUnPYkHHUdrLZ+3KtqQi947O/wxcdJV5hVOf/3JUk5xDClrOhWlmrR8cR06ADf3Q9Oug+Ovws22AYe+L+4V9X9vyuaXlV58/clSTnAMKWsGFNZQaq0ZJljqdISxlRWJFRRE0KAHgNgyFQ45VHYfC94/Lz4SdWdP4YPX0+6wnaVd39fkpSg4p1hq6xqmLScl6vDNtwWDr8SPng1DlTVV8PsK6Hn4HgF4PpbJ11hm8vrvy9JyjJX80kt9cnbMOMCmPUP+GohbFEJu/4Yuu+UdGWSpDbkaj6pvayxIXz/jzB6Huz5K6idBVdUwhX7wr/vhSz+D4okKXmGKam1UmvB7mPirur7/gU+/i/88/C4s/rcm2Hx10lXKEnKAsOUlKlVusBOp8a9qgZNhsWL4JaT4Pzt4i7ri75IukJJUjsyTEltpaQUeh8Dp8+AI6+FLuvCXT9J96o6t+h6VUlSsTBMSW2tQwf43gHwowdg6B2wQU944PfpXlW/h4XvJV2hJKkNGaak9hICbLobDLkVhj0Mmw2E6ePjJ1V3/QQ+ejPpCiVJbcAwJWVDtz5wxBQ4cxZscwTMngIT+8AtJ8O7LyRdnSQpA/aZUl6aVl2b3w0lP5kPT6Z7VS36DLbcFwb8GLrvmHRlkqS05vaZMkwp70yrrmXs1LnUL1r8zbFUaQlnD+6VX4EK4PMP4xV/MyZD/YfQfZe4Aejme8fDhJKkxBimVLD6n/MgtXX1/3O8vCzF42cNTKCixrXo6dlXn8EzV8ET58MntbB+LxgwCrYaBCW5uetT3j8dlKQm2AFdBWt+I0FqZceT0PD0rLaungioratn7NS5TKuubfwFq3SBnU6DEXPg4Ath8Zdxr6pJfWHWFTnXq6rFv58kFTDDlPJOt7JUi44nYVxVzTLDkAD1ixYzrqpm5S/suAr0ORZOnwlHXA2pteHO0XDeNjB9AnzxSTtW3Xyt/v0kqQAZppR3xlRWkCotWeZYqrSEMZUVCVX0vzJ+etahA2x1EJz8IPzwdlhvK7j/t3Gvqgf+DxYuaMNqWy4fng5KUrbk5mQMaSUa5uW05Xydtp7/060s1ei8rpU9PVthDd/ZPf6YXx33qXrs3HglYJ8hsMtwWGuTVtfZWq35/SSpUDkBXUWvPVYHtvQ9W3T++y/D4+fBs9dDtAR6HQYDRsN632tVra1RUCsqJWkFnIAuNVN7zP8Z1Kecswf3orwsRSBeabiyoNGiGtbdAg6eBCOfjSetv3gnXLgTXHc0vPVUq2tuiZb+fpJUyHwypaK36Vl30dg/BQF4/Zz9c7+Gzz+Epy6BmRdB/UewyYD4SdXme9mrSpIy4JMpqZlyYXVgRjWsujbscRaMfh4qz4aPXodrD4WLd4V5t8CSxU2/hySp1QxTKnq5sDqwTWpYpQvsfHrcq+qgSXFvqptPhPP7xtvWfP1lG1ctSQLDlJQT83/atIaOq8B2Q+CMhl5VZXDnKJiwTTxx/ctP27x+SSpmzpmSCl0UweuPxC0VXn8EOq8J258cT17vsm7S1UlSznLOlKRYCPCdPWDo7XET0E13h8f+HjcAvXsM1P0n6QolKa8ZpqRiUt4XjrwazngKeh4a7/s3sQ/ceiq892LS1UlSXjJMScWo65Yw6IK4V9UOw+CF29K9qo6Bt55OujpJyitNhqkQwsYhhIdCCC+EEJ4PIYxMH187hHBfCOHl9J9rtX+5ktrUmhvBvmfHbRV2PwvefBwu3xuuPABeuT+ebyVJWqnmPJn6GvhJFEVbATsBZ4QQtgLOAh6IomgL4IH015Ly0aprw55j072q/gwfvArXHAoX7wbzptqrSpJWoskwFUXR21EUPZP+/FPgRaAcOBiYkj5tCjCovYqUlCWdVoOdz4CRDb2qPoebT4BJ/WD2lfaqkqRGtGjOVAihB9AHmAmsH0XR2+lvvQOs36aVSUpOx07pXlVPwRFXQac14I6RcN628MT59qqSpKU0u89UCGE14BHgT1EUTQ0h1EVRVLbU9z+Kouh/5k2FEIYBwwC6d+/e980332ybyqUcNa26lnFVNcyvq6dbWYoxlRX5vwFwFMFrD8H08fD6o9C5LJ64vuMp9qqSVLCa22eqWWEqhFAK3AlURVF0bvpYDbBHFEVvhxA2BB6Oomile1/YtFOFblp1LWOnzqV+0bdzjFKlJVnvqN6u/jsbpp8LL90JHVPQdyjsfCaUbZx0ZZLUptqsaWcIIQCXAy82BKm024Gh6c+HAre1plCpkIyrqlkmSAHUL1rMuKqahCpqBxv1haOuTfeqGgxPXwYTe8Otp8F7LyVdnSRlXXPmTPUHhgADQwhz0h/7AecA+4QQXgb2Tn8tFbX5dfUtOp7XulbAoAvjjZW3/xE8fytcuCNcfyz81yfQkopHx6ZOiKJoOhBW8O292rYcKb91K0tR20hw6laWSqCaLCnbGH7wF9jtZ/DUxTDz4ngIsMeusOuP4Tt7xlvaSFKBsgO61IbGVFaQKi1Z5liqtIQxlSudTlgYuqwDe/4CRs+D7/8RPngFrj4ELtkDnp9mrypJBavZq/naghPQVQwKcjVfa3z9JTx3A0yfAB++CutsDv1HwjZHxq0XWqEQr20h/k5SoWjT1XxtxTAlFaEli+HFO+IVgG8/C6t3ixuD9j0+bhLaTIW4UrIQfyepkLTZaj5JykiHEth6EAx7BI6bCutsBvf+EsZvDQ/9GT77oFlvU4grJQvxd5KKkWFKUnaEAJvvBcffCSfdD5v0h0f+AhN6wj1nQd1bK315Ia6ULMTfSSpGhilJ2bfx9nD0P+H0mbDVwfD0pXGvqmmnw4LGn8qsaEVkPq+ULMTfSSpGhilJyVnvu3DIRTCiGvqdBPOmwgUNvapmL3PqmMoKSkuWbbFQWhLyeqVkUa/+lAqIYUpS8sq6w35/jdsq7PZTeOMxuGwgTDkQXn0o3hsQYPn1MtlbP9MuBvUp5+zBvSgvSxGA8rKUk8+lPORqPkm558tPYdY/4MkLYOE70K0Pv1ywD//8dBui5f4fsLwsxeNnDUyoUEmFzNV8kvJXp9Wh/wgY9RwceB588TF/WvRX7l9lDIeXPEwpX39zqpO1JSXNMCUpd3XsFPejOnMWvyr9KfV0YlzpJTzaaRQnldzNqnzhZG1JiTNMScp9HUrot9+JHB6dww+/+jlvLNmAX5dewxOdRnDZJvfD5x8mXaGkItbkRseSlAsaJmWPq+rEMXXbsk/qP/x+7Sq+V3MBjP9H/ARr5zNhTSdvS8ouJ6BLyksNe9p1+fhlRqfupjJ6jA6hQ7z3X/+R0HXLpEuUlOecgC6pYDXsaVdbV8+/o4047fNh7PP1BF7b5HCYdzNcsAPccBzUPpN0qZKKgGFKUt5pbE+7Vxetw5C3D4dR82DXn8Brj8Kle8KUg+C1h7/tVSVJbcwwJSnvrHRPu9W6wl6/jhuA7vOHeHuaqw6GSwfCC7fDkiVZrlZSoTNMScob06pr6X/OgytsfL5Mm4TOa8S9qkY+CwdMgPqP4MYh8RBg9TXw9VdZqVlS4TNMScoLS8+TaswK97Qr7Qz9ToDhs+GwK+KvbzsDJvaBGZPhq8/auXJJhc4wJSkvNDZPqkGz9rTrUAI9D4VTHoNjb4G1NoF/nQXje8LD59irSlKr2WdKUl5Y0TypAC3bmy8E2GLv+OM/M2H6eHj4bHh8YrpX1Rn2qpLUIj6ZkpQXVrRtTEbbyXTfEY65Hk57Ar67P8y8CM7bFm47E95/ufXvK6moGKYk5YUxlRWkSkuWObbCeVIttf7WcOilMOKZ+OnU3Jtg0vZwwxCYX93kyxsmxm961l30P+dBplXXZl6TpLxhB3RJeaOh6/n8unq6laUYU1mx8nlSrbVwAcycDE9dBl9+DN/ZEwaMhk13i4cJl6tp7NS5y8znSpWWND2HS1LOa24HdMOUJK3IF5/ArCtgxoWw8F0o7xuHqor9oUP8YL//OQ82usKwvCzVsrlcknKO28lIUqY6rwEDRsHI5+CA8fD5B/E2NRfuBNXXwuJFK28gKqkoGKYkqSmlnaHfiXDmbDj0cigphdtOh/N6M2K1B0jxxf+8JKOJ8ZLyimFKkpqrpCP0OgxOnQ7H3ARl3Rn99eU80WkkI0qmsiYLgTacGC8pLzhnSpIy8Z8ZvHPXn9ng3UdYGHXm9o6VrL3XKPbdZbukK5OUISegS1I2vTMPHj8P5t0Sd1vf9ijYZSSsu3nSlUlqJSegS1I2bdAz7lU1fDZs90N49gaY1A9uHArz5yRdnaR2ZJiSpLa09qaw/99h9Ly4jcKrD8Ilu8PVh8Drj0EWRwMkZYdhSpLaw2rrwd6/jUPV3r+LhwGnHACX7Q0v3QVLliRdoaQ24pwpScqGRfUw55/xvKq6N6Hrd6H/qHh1YElpVkvJWid5Kc85Z0qScklpCrY/CYY/A4Mvg1AC006FiX1g5sXw1edZKaNh+5vaunoioLaunrFT57qfoJQBw5QkZVNJR9jmcDjtcTjmRlijHO75GUzoCY+Mg/qP2vXHj6uqWWYfQYD6RYsZV1XTrj9XKmSGKUlKQgiwZSWcVAUn/AvK+8FDf4TxPeHeX8Enb7fLj3X7G6ntGaYkKWmb7AzH3ginPg4VP4AnL4DztoHbR8AHr7bpj1rRNjdufyO1nmFKknLFBj3h0MvieVV9hsCz18e9qm46Ht5+tk1+xJjKClKlJcscc/sbKTOGKUnKNWtvCgecC6Pmwi4j4OX74eLd4JpD4Y3pGfWqGtSnnLMH96K8LEUAystSnD24l6v5pAzYGkHSSrmMPgfU18Gsy2HGZPhsAWy0Q9wQdMt9oYP/Tyy1F1sjSMqYy+hzRKoMdv1J/KRqv7/Bwnfg+qNh8i7xUODiRUlXKBU1w5SkFXIZfY4pTcEOJ8Pwahh8abwi8NZTYOJ2MPOSrPWqkrQsw5SkFcrnZfTTqmvpf86DbHrWXfQ/58HCeppW0hG2OSJe/Xf09bDGhnDPGJjQCx4dFw8LSsoaw5SkFcrXZfRFMzzZoUPcSuHEKjjhHijfDh5s6FX1a/j0naQrlIqCYUrSCuXrMvqiG54MATbZBY69CU6dDlt+H56cBBO2gTtGwYevJV2hVNA6Jl2ApNxGgw0PAAAS0klEQVTVsGovV1bzNXdlYa4MTyayEnKDXnDYFTDwV/DE+VB9LTwzBbY+JN5YecNt2vfnS0XI1giS8kLD0N3ST5xSpSWN9kjqf86D1DYSnMrLUjx+1sB2rxVaVm+7+vQdmHEhPH0FfPUpbL5P3FZhk13iJ1qSVqjNWiOEEK4IIbwXQpi31LHfhRBqQwhz0h/7ZVqwJK1MS4bucmF4MmeGGlffAPb5Pxg9Dwb+GuZXw5X7wRWVUHMPLFmS3XqkAtScOVNXAvs2cnx8FEW90x93t21ZkrSslgzdrazLd7ZW+eXKUOM3UmWw20+/7VX1ydtw3VFwUX949gZ7VUkZaHLOVBRFj4YQerR/KZK0Yt3KUo0O3a1oZeGgPuX/M5y2/NBbwyq/hvOTrDdrVlk17lXV93iYNxWmj4dbh8FDf4y3rulzXNzPSlKzZbKab3gI4bn0MOBabVaRJDWiLYbusjn0lgtDjStVUgrbHgmnPRH3qlptA7j7p3FbhUf/Zq8qqQVaG6YmA98BegNvA39f0YkhhGEhhFkhhFkLFixo5Y+TVOzaYoPebA695c2Gwg29qk66F46/G7r1hgf/EIeq+34Dn76bdIUqAvneZLdZq/nSw3x3RlHUsyXfW56r+SQlKRdW+eWFt5+Lh/9emAYdSqHPsfEQ4NqbJl2ZClDOrHxtRLtudBxC2HCpLw8B5q3oXEnKFTk/9JYrNtwGDv8HnDkLeh8N1dfA+dvBzSfBO/7rXm0rZ1a+ZqDJCeghhOuAPYB1Qwj/BX4L7BFC6A1EwBvAKe1YoyS1iVxrQprz1tkMDjwP9hgLT14As66AeTfDFt+HAT+GTXZOukIVgJxb+doKNu2UJDVP/Ufw9GUwYzJ8/gFsvFPcAHTLShuAqtVyefi9XYf5JClf5ftE10Sl1oLdxsCoefCDcfBJLVx3JEzuD8/dCIu/TrpC5aFCGH43TEkqGg0TXWvr6on4ts+UgaqFVlkVdhwGI6rhkIshWgxTT4bz+8BTl8Ki/BmeUfLyZuXrSjjMJ6lo5PJwQl5bsgT+/S+Yfi7892nosh7sdBpsfxJ0XjPp6qRWc5hPkpZTCBNdc1KHDvDd/eCk+2DonbBBL3jg93Gvqvt/BwvfS7pCqV0ZpiQVjRVt5ZL4Fi+FIgTYdFcYMhWGPQKb7wXTJ8Sh6s4fw4evJ12h1C4MU5KKRiFMdM0b3XrD4VfC8NnxtjXVV8P5feGWH8G7zyddndSmDFOSikYhTHTNO+tsBgedDyOfjedR1dwDk3eBa4+AN59MujqpTTgBXZKUPZ9/+G2vqvoPofvOcQPQLfaxV5VyjhPQJUm5Z9W1Yfefweh5sO9foO4t+OfhcNEAmHuzvaqUlwxTkqTsW6UL7HQqjJwDgybD4kVwy0kwqS88fTks+iLpCqVmM0xJkpJTUgq9j4HTZ8CR18Kq68JdP4YJvWD6ePjik6QrlJrknClJSsi06lo3XV5eFMEbj8Fj58JrD0GnNePmnzudBqutl3R1KjLNnTNlmJKkBDRsbVO/aPE3x1KlJa4uXNr86rhP1Qu3QcdO0Oc42GU4rNUj6cpUJJyALkk5bFxVzTJBCqB+0WLGVdUkVFEO6tYHjpgCZ86CXofD7CkwcTu45WR7VSmnGKYkKQFubdMC624OB0+CUc/Fw30v3RX3qvrnkfCfmUlXJxmmJCkJbm3TCmt0g8o/xW0V9vgFvPUUXPF9uOIH8PJ98XwrKQGGKUlKgFvbZGDVtWGPn6d7VZ0DdW/CtYfBRbvaq0qJMExJUgLc2qYNrNIlHvYb0dCr6st0r6p+MOsKe1Upa1zNJ0kqDEuWQM1dcVuF+c/AauvDzmdA3xOg8xpJV6c85Go+SVJx6dABvncgnPwg/PB2WO97cN9vYEJPeOAPsHBB0hWqQBmmJEmFJQT4zu7ww9vg5Idg093hsb/Hoequn8JHbyZdoQqMYUqSVLjKt4Mjr4YznoJeh8HsK2FiH5g6DN59IenqVCAMU5Kkwtd1Szj4gnhj5R1PgRfvgMk7wz+PilssSBkwTEmSiseaG8G+Z8Po52GPsfDWDLh8H/jH/vDy/faqUqsYpiRJxWfVtWGPs+JQVXk2fPQ6XHsoXLwbzJsKSxY3/R5SmmFKklS8VukCO58e96o6+AJYVA83nxD3qpp9JXz9ZdIVKg8YpiRJ6rgK9DkOzpgJR1wNndeEO0bChG3g8Ynw5adJV6gcZpiSJKlBhxLY6qC4pcKQadC1Au77NYzfGh78I3z2ftIVKgcZpiRJWl4IsNmeMPT2uAnoprvBo3+D8T3h7jFQ95+kK1QOMUxJkrQy5X3hyGviIcCeg+N9/87rDVNPgfdeTLo65QDDlCRJzdG1AgZdCCOfTfequh0u3AmuOwbeejrp6pQgw5QkSS3R0Ktq1DzY/Sx483G4fG+48gB4xV5VxcgwJUlSa3RZB/Ycm+5V9Wf44FW45lC4ZHd4/lZ7VRURw5QkSZnotBrsfEa8Vc1Bk+Crz+Cm42HS9jB7ir2qioBhSpKkttCxE2w3JN5U+YiroNPqcMcIOG9beOJ8e1UVMMOUJEltqUMJbHUwDHsYhtwK624B9/4qbqvw4J/sVVWADFOSJLWHEGCzgTD0DvjRA9BjADz61zhU3fNzqHsr6QrVRgxTkiS1t436wVHXxkOAPQfD05fBxN5w62nw3ktJV6cMGaYkScqWhl5VI+bA9ifDC9Pgwh3h+mPhv7OSrk6tZJiSJCnbyjaGH5wT96ra7WfwxnS4bK+4V9WrD9qrKs8YpiRJSkqXdWDgL2H0PPj+H+GDV+DqQ+CSPeD5afaqyhOGKUmSktZpddhleLxVzYET4zYKNw2FC3aAZ66yV1WOM0xJkpQrOnaCvkPhzKfh8CuhdFW4fXi8sfITk+DLhUlXqEYYpiRJyjUdSmDrQ+CUR+G4qbDOZnDvL2H81vDQn+GzD5KuUEsJURYnufXr1y+aNcvVCpIktdhbT8P08VBzV/zEaruhsMuZ8cbLOWZadS3jqmqYX1dPt7IUYyorGNSnPOffe3khhNlRFPVr8jzDlCRJeeS9l+DxCTD3pvjrbY6E/qOg65bJ1pU2rbqWsVPnUr/o28nzqdISzh7cK+PQ057v3ZjmhimH+SRJyifrfRcOuQhGVEO/k2De1Hii+vXHQu3spKtjXFXNMmEHoH7RYsZV1eT0e2fCMCVJUj4q6w77/TVuq7DbT+GNx+DSgTDlwER7Vc2vq2/R8Vx570w0GaZCCFeEEN4LIcxb6tjaIYT7Qggvp/9cq33LlCRJjeqyLgz8VdwAdJ8/wIJ/J9qrqltZqkXHc+W9M9GcJ1NXAvsud+ws4IEoirYAHkh/LUmSktJ5Deg/AkY9BweeB19+slSvqqvh66+yUsaYygpSpSXLHEuVljCmsiKn3zsTTYapKIoeBT5c7vDBwJT051OAQW1clyRJao2OnaDv8XDmLDjsH1CagtvPhPO2hScvaPdeVYP6lHP24F6Ul6UIQHlZqs0miLfne2eiWav5Qgg9gDujKOqZ/rouiqKy9OcB+Kjh65VxNZ8kSVkWRfDqA/DYeHhzOqTWgh1OgR1PgVXXTrq6nJa11XxRnMZWmMhCCMNCCLNCCLMWLFiQ6Y+TJEktEQJsvjeccBecdB903xkeOQfG94R//QI+rk26wrzX2jD1bghhQ4D0n++t6MQoii6JoqhfFEX9unbt2sofJ0mSMrbxDnD0dXD6DNjqIJh5UTz8d9sZ8P7LSVeXt1obpm4HhqY/Hwrc1jblSJKkdrfe9+JeVSPnQL8TYO7NMGl7uOG4nOhVlW+anDMVQrgO2ANYF3gX+C0wDbgR6A68CRwRRdHyk9T/h3OmJEnKQQsXxE+pnroUvvwYNt0ddv1x/GcISVeXGLeTkSRJLfPFJzDrCphxISx8F7ptBwNGw3cPgA7F1+fb7WQkSVLLdF4DBoyCkc/BAeOh/kO4cUjcq6r6mqz1qso3hilJkrSs0s7Q70Q4czYcejl07BxPUp/YG568EL76LOkKc4phSpIkNa6kI/Q6DE59DI69GdbqAVVj47YKD/8FPm9yunRRMExJkqSVCwG22AdOuBtOvBc23hEe/nMcqqp+WfS9qgxTkiSp+brvCMdcD6c9Cd87AGZMTveqOrNoe1UZpiRJUsutvxUMvgRGPBPvBTj3pnSvqiEwvzrp6rLKMCVJklpvrR6w/99g1Ly4N9Vrj8Ale8BVB8efZ7EFU1IMU5IkKXOrdYW9fgOj58Lev4N3X4CrDoLL9oIX74AlS5KusN0YpiRJUtvpvGbc6HPUXNj/XPjs/Xibmgt3hOprYfGipCtsc4YpSZLU9ko7w/YnwfBn4l5VJavAbafDeb1hxkUF1avKMCVJktrPN72qpsMxN0FZd/jXz+O2Co/8tSB6VRmmJElS+wsBtvw+nHgPnFgFG+8AD/0JJvSKe1V9Mj/pClvNMCVJkrKr+05wzA1w2hNQsV/cq2rCNuleVa8kXV2LGaYkSVIy1t8aDr0Uhs+GvkPhuRthUj+4cSjMn5N0dc1mmJIkSclae1PY/+8weh4MGAWvPgiX7A5XHwKvP5rzvaoMU5IkKTestl7co2r0PNjrt/DOXJhyIFy2N7x0V872qjJMSZKk3NJ5zbib+qi58ROrzxbA9cfA5J1hznU516vKMCVJknJTaQq2/1Hcq2rwZRBKYNqpMLEPvPJA0tV9wzAlSZJyW0lH2OZwOO1xOOZGWHMjWH2DpKv6RsekC5AkSWqWEGDLyvgjh/hkSpIkKQOGKUmSpAwYpiRJkjJgmJIkScqAYUqSJCkDhilJkqQMGKYkSZIyYJiSJEnKgGFKkiQpA4YpSZKkDBimJEmSMmCYkiRJyoBhSpIkKQMhiqLs/bAQFgBvtuOPWBd4vx3fP9cV++8PXgPwGoDXALwG4DUAr0Gmv/8mURR1beqkrIap9hZCmBVFUb+k60hKsf/+4DUArwF4DcBrAF4D8Bpk6/d3mE+SJCkDhilJkqQMFFqYuiTpAhJW7L8/eA3AawBeA/AagNcAvAZZ+f0Las6UJElSthXakylJkqSsyrswFUK4IoTwXghh3gq+v0cI4eMQwpz0x2+yXWN7CiFsHEJ4KITwQgjh+RDCyEbOCSGEiSGEV0IIz4UQtkui1vbSzGtQ6PdB5xDCUyGEZ9PX4PeNnFPo90FzrkFB3wcAIYSSEEJ1COHORr5X0PdAgyauQTHcA2+EEOamf79ZjXy/4O+DZlyDdr0POrblm2XJlcAk4KqVnPNYFEUHZKecrPsa+EkURc+EEFYHZocQ7oui6IWlzvkBsEX6Y0dgcvrPQtGcawCFfR98CQyMomhhCKEUmB5CuCeKohlLnVPo90FzrgEU9n0AMBJ4EVijke8V+j3QYGXXAAr/HgDYM4qiFfVTKpb7YGXXANrxPsi7J1NRFD0KfJh0HUmJoujtKIqeSX/+KfG/QMqXO+1g4KooNgMoCyFsmOVS200zr0FBS//dLkx/WZr+WH4CZKHfB825BgUthLARsD9w2QpOKeh7AJp1DVQE90HS8i5MNdMu6UeZ94QQtk66mPYSQugB9AFmLvetcuCtpb7+LwUaNlZyDaDA74P00MYc4D3gviiKiu4+aMY1gMK+DyYAPwOWrOD7BX8P0PQ1gMK+ByD+n4j7QwizQwjDGvl+MdwHTV0DaMf7IB+H+ZryDNA9/eh/P2Aa8aPNghJCWA24BRgVRdEnSdeThCauQcHfB1EULQZ6hxDKgFtDCD2jKGp0LmGhasY1KNj7IIRwAPBeFEWzQwh7JF1PEpp5DQr2HljKgCiKakMI6wH3hRBeSo/iFJOmrkG73gcF92QqiqJPGh79R1F0N1AaQlg34bLaVHp+yC3AtVEUTW3klFpg46W+3ih9rGA0dQ2K4T5oEEVRHfAQsO9y3yr4+6DBiq5Bgd8H/YGDQghvANcDA0MI1yx3TqHfA01egwK/BwCIoqg2/ed7wK3ADsudUuj3QZPXoL3vg4ILUyGEDUIIIf35DsS/4wfJVtV20r/b5cCLURSdu4LTbgd+mF7BsRPwcRRFb2etyHbWnGtQBPdB1/TTGEIIKWAf4KXlTiv0+6DJa1DI90EURWOjKNooiqIewFHAg1EUHbfcaQV9DzTnGhTyPQAQQuiSXohDCKEL8H1g+SfUBX0fNOcatPd9kHfDfCGE64A9gHVDCP8Ffks88ZQoii4CDgNOCyF8DdQDR0WF1Zm0PzAEmJueKwLwC6A7fHMN7gb2A14BPgdOSKDO9tSca1Do98GGwJQQQgnxvxRujKLozhDCqVA090FzrkGh3wf/o8jugUYV2T2wPvEQN8T/Tf9nFEX/KrL7oDnXoF3vAzugS5IkZaDghvkkSZKyyTAlSZKUAcOUJElSBgxTkiRJGTBMSZIkZcAwJUmSlAHDlCRJUgYMU5IkSRn4f+Z3zKVxbE/9AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1,1, figsize=(10,6))\n", "ax.plot(y_test, predicted_y, 'o')\n", "grid = np.linspace(np.min(dfcars.mpg), np.max(dfcars.mpg), 100)\n", "ax.plot(grid, grid, color=\"black\") # 45 degree line\n", "ax.set_xlabel(\"actual y\")\n", "ax.set_ylabel(\"predicted y\")\n", "\n", "fig1, ax1 = plt.subplots(1,1, figsize=(10,6))\n", "ax1.plot(dfcars.wt, dfcars.mpg, 'o')\n", "xgrid = np.linspace(np.min(dfcars.wt), np.max(dfcars.wt), 100)\n", "ax1.plot(xgrid, regression.predict(xgrid.reshape(100, 1)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## $k$-nearest neighbors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Great, so we did a simple linear regression on the car data.\n", "\n", "Now that you're familiar with `sklearn`, you're ready to do a KNN regression. Let's use $5$ nearest neighbors." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from sklearn.neighbors import KNeighborsRegressor\n", "knnreg = KNeighborsRegressor(n_neighbors=5)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "knnreg.fit(X_train, y_train)\n", "r2 = knnreg.score(X_test, y_test)\n", "r2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise
\n", "What is the $R^{2}$ score on the training set?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# solution\n", "knnreg.score(X_train, y_train)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Lets vary the number of neighbors and see what we get." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "regdict = {}\n", "# Do a bunch of KNN regressions\n", "for k in [1, 2, 4, 6, 8, 10, 15]:\n", " knnreg = KNeighborsRegressor(n_neighbors=k)\n", " knnreg.fit(X_train, y_train)\n", " regdict[k] = knnreg # Store the regressors in a dictionary" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Now let's plot it all\n", "fig, ax = plt.subplots(1,1, figsize=(10,6))\n", "\n", "ax.plot(dfcars.wt, dfcars.mpg, 'o', label=\"data\")\n", "\n", "xgrid = np.linspace(np.min(dfcars.wt), np.max(dfcars.wt), 100)\n", "for k in [1, 2, 6, 10, 15]:\n", " predictions = regdict[k].predict(xgrid.reshape(100,1))\n", " if k in [1, 6, 15]:\n", " ax.plot(xgrid, predictions, label=\"{}-NN\".format(k))\n", "\n", "ax.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice how the $1$-NN goes through every point on the training set but utterly fails elsewhere. Lets look at the scores on the training set." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ks = range(1, 15) # Grid of k's\n", "scores_train = [] # R2 scores\n", "for k in ks:\n", " knnreg = KNeighborsRegressor(n_neighbors=k) # Create KNN model\n", " knnreg.fit(X_train, y_train) # Fit the model to training data\n", " score_train = knnreg.score(X_train, y_train) # Calculate R^2 score\n", " scores_train.append(score_train)\n", "\n", "# Plot\n", "fig, ax = plt.subplots(1,1, figsize=(12,8))\n", "ax.plot(ks, scores_train,'o-')\n", "ax.set_xlabel(r'$k$')\n", "ax.set_ylabel(r'$R^{2}$')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why do we get a perfect $R^2$ at k=1?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
Exercise (5 min)
\n", "* Make the same plot as above on the *test* set.\n", "* What is the best $k$?" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "ks = range(1, 15) # Grid of k's\n", "scores_test = [] # R2 scores\n", "for k in ks:\n", " knnreg = KNeighborsRegressor(n_neighbors=k) # Create KNN model\n", " knnreg.fit(X_train, y_train) # Fit the model to training data\n", " score_test = knnreg.score(X_test, y_test) # Calculate R^2 score\n", " scores_test.append(score_test)\n", "\n", "# Plot\n", "fig, ax = plt.subplots(1,1, figsize=(12,8))\n", "ax.plot(ks, scores_test,'o-', ms=12)\n", "ax.set_xlabel(r'$k$')\n", "ax.set_ylabel(r'$R^{2}$')" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.1" } }, "nbformat": 4, "nbformat_minor": 1 }