{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <img style=\"float: left; padding-right: 10px; width: 45px\" src=\"https://raw.githubusercontent.com/Harvard-IACS/2018-CS109A/master/content/styles/iacs.png\"> CS-109A Introduction to Data Science\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lab 4: Multiple and Polynomial Regression (September 26, 2019 version) \n",
    "\n",
    "**Harvard University**<br/>\n",
    "**Fall 2019**<br/>\n",
    "**Instructors:** Pavlos Protopapas, Kevin Rader, and Chris Tanner<br/>\n",
    "**Lab Instructor:** Chris Tanner and Eleni Kaxiras<br/>\n",
    "**Authors:** Rahul Dave, David Sondak, Will Claybaugh, Pavlos Protopapas, Chris Tanner\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>\n",
       "blockquote { background: #AEDE94; }\n",
       "h1 { \n",
       "    padding-top: 25px;\n",
       "    padding-bottom: 25px;\n",
       "    text-align: left; \n",
       "    padding-left: 10px;\n",
       "    background-color: #DDDDDD; \n",
       "    color: black;\n",
       "}\n",
       "h2 { \n",
       "    padding-top: 10px;\n",
       "    padding-bottom: 10px;\n",
       "    text-align: left; \n",
       "    padding-left: 5px;\n",
       "    background-color: #EEEEEE; \n",
       "    color: black;\n",
       "}\n",
       "\n",
       "div.exercise {\n",
       "\tbackground-color: #ffcccc;\n",
       "\tborder-color: #E9967A; \t\n",
       "\tborder-left: 5px solid #800080; \n",
       "\tpadding: 0.5em;\n",
       "}\n",
       "\n",
       "span.sub-q {\n",
       "\tfont-weight: bold;\n",
       "}\n",
       "div.theme {\n",
       "\tbackground-color: #DDDDDD;\n",
       "\tborder-color: #E9967A; \t\n",
       "\tborder-left: 5px solid #800080; \n",
       "\tpadding: 0.5em;\n",
       "\tfont-size: 18pt;\n",
       "}\n",
       "div.gc { \n",
       "\tbackground-color: #AEDE94;\n",
       "\tborder-color: #E9967A; \t \n",
       "\tborder-left: 5px solid #800080; \n",
       "\tpadding: 0.5em;\n",
       "\tfont-size: 12pt;\n",
       "}\n",
       "p.q1 { \n",
       "    padding-top: 5px;\n",
       "    padding-bottom: 5px;\n",
       "    text-align: left; \n",
       "    padding-left: 5px;\n",
       "    background-color: #EEEEEE; \n",
       "    color: black;\n",
       "}\n",
       "header {\n",
       "   padding-top: 35px;\n",
       "    padding-bottom: 35px;\n",
       "    text-align: left; \n",
       "    padding-left: 10px;\n",
       "    background-color: #DDDDDD; \n",
       "    color: black;\n",
       "}\n",
       "</style>\n",
       "\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## RUN THIS CELL TO GET THE RIGHT FORMATTING \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": [
    "## Table of Contents\n",
    "\n",
    "<ol start=\"0\">\n",
    "<li> Learning Goals / Tip of the Week / Terminology </li>\n",
    "<li> Training/Validation/Testing Splits (slides + interactive warm-up)</li>\n",
    "<li> Polynomial Regression, and Revisiting the Cab Data</li>\n",
    "<li> Multiple regression and exploring the Football data </li>\n",
    "<li> A nice trick for forward-backwards </li>\n",
    "    <li> Cross-validation</li>\n",
    "</ol>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Goals\n",
    "After this lab, you should be able to\n",
    "- Explain the difference between train/validation/test data and WHY we have each.\n",
    "- Implement cross-validation on a dataset\n",
    "- Implement arbitrary multiple regression models in both SK-learn and Statsmodels.\n",
    "- Interpret the coefficent estimates produced by each model, including transformed and dummy variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.api import OLS\n",
    "\n",
    "from sklearn import preprocessing\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from sklearn.metrics import r2_score\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "from pandas.plotting import scatter_matrix\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra Tip of the Week\n",
    "\n",
    "Within your terminal (aka console aka command prompt), most shell environments support useful shortcuts:\n",
    "\n",
    "<ul>\n",
    "  <li>press the [up arrow] to navigate through your most recent commands</li>\n",
    "  <li>press [CTRL + A] to go to the beginning of the line</li>\n",
    "  <li>press [CTRL + E] to go to the end of the line</li>\n",
    "  <li>press [CTRL + K] to clear the line</li>\n",
    "  <li>type `history` to see the last commands you've run</li>\n",
    "</ul>  \n",
    "\n",
    "\n",
    "## Terminology\n",
    "\n",
    "Say we have input features $X$, which via some function $f()$, approximates outputs $Y$. That is, $Y = f(X) + \\epsilon$ (where $\\epsilon$ represents our unmeasurable variation (i.e., irreducible error).\n",
    "\n",
    "- **Inference**: estimates the function $f$, but the goal isn't to make predictions for $Y$; rather, it is more concerned with understanding the relationship between $X$ and $Y$.\n",
    "- **Prediction**: estimates the function $f$ with the goal of making accurate $Y$ predictions for some unseen $X$.\n",
    "\n",
    "\n",
    "We have recently used two highly popular, useful libraries, `statsmodels` and `sklearn`.\n",
    "\n",
    "`statsmodels` is mostly focused on the _inference_ task. It aims to make good estimates for $f()$ (via solving for our $\\beta$'s), and it provides expansive details about its certainty. It provides lots of tools to discuss confidence, but isn't great at dealing with test sets.\n",
    "\n",
    "`sklearn` is mostly focused on the _prediction_ task. It aims to make a well-fit line to our input data $X$, so as to make good $Y$ predictions for some unseen inputs $X$. It provides a shallower analysis of our variables. In other words, `sklearn` is great at test sets and validations, but it can't really discuss uncertainty in the parameters or predictions.\n",
    "\n",
    "\n",
    "- **R-squared**: An interpretable summary of how well the model did. 1 is perfect, 0 is a trivial baseline model based on the mean $y$ value, and negative is worse than the trivial model.\n",
    "- **F-statistic**: A value testing whether we're likely to see these results (or even stronger ones) if none of the predictors actually mattered.\n",
    "- **Prob (F-statistic)**: The probability that we'd see these results (or even stronger ones) if none of the predictors actually mattered. If this probability is small then either A) some combination of predictors actually matters or B) something rather unlikely has happened\n",
    "- **coef**: The estimate of each beta. This has several sub-components:\n",
    "  - **std err**: The amount we'd expect this value to wiggle if we re-did the data collection and re-ran our model. More data tends to make this wiggle smaller, but sometimes the collected data just isn't enough to pin down a particular value.\n",
    "  - **t and P>|t|**: similar to the F-statistic, these measure the probability of seeing coefficients this big (or even bigger) if the given variable didn't actually matter. Small probability doesn't necessarily mean the value matters\n",
    "  - **\\[0.025 0.975\\]**: Endpoints of the 95% confidence interval. This is a interval drawn in a clever way and which gives an idea of where the true beta value might plausibly live. (If you want to understand why \"there's a 95% chance the true beta is in the interval\" is _wrong_, start a chat with Will : )\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Polynomial Regression, and Revisiting the Cab Data\n",
    "\n",
    "Polynomial regression uses a **linear model** to estimate a **non-linear function** (i.e., a function with polynomial terms). For example:\n",
    "\n",
    "$y = \\beta_0 + \\beta_1x_i + \\beta_1x_i^{2}$\n",
    "\n",
    "It is a linear model because we are still solving a linear equation (the _linear_ aspect refers to the beta coefficients)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>TimeMin</th>\n",
       "      <th>PickupCount</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>860.0</td>\n",
       "      <td>33.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>17.0</td>\n",
       "      <td>75.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>486.0</td>\n",
       "      <td>13.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>300.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>385.0</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   TimeMin  PickupCount\n",
       "0    860.0         33.0\n",
       "1     17.0         75.0\n",
       "2    486.0         13.0\n",
       "3    300.0          5.0\n",
       "4    385.0         10.0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# read in the data, break into train and test\n",
    "cab_df = pd.read_csv(\"../data/dataset_1.txt\")\n",
    "train_data, test_data = train_test_split(cab_df, test_size=.2, random_state=42)\n",
    "cab_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1250, 2)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cab_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# do some data cleaning\n",
    "X_train = train_data['TimeMin'].values.reshape(-1,1)/60 # transforms it to being hour-based\n",
    "y_train = train_data['PickupCount'].values\n",
    "\n",
    "X_test = test_data['TimeMin'].values.reshape(-1,1)/60 # hour-based\n",
    "y_test = test_data['PickupCount'].values\n",
    "\n",
    "def plot_cabs(cur_model, poly_transformer=None):\n",
    "    \n",
    "    # build the x values for the prediction line\n",
    "    x_vals = np.arange(0,24,.1).reshape(-1,1)\n",
    "    \n",
    "    # optionally use the passed-in transformer\n",
    "    if poly_transformer != None:\n",
    "        dm = poly_transformer.fit_transform(x_vals)\n",
    "    else:\n",
    "        dm = x_vals\n",
    "        \n",
    "    # make the prediction at each x value\n",
    "    prediction = cur_model.predict(dm)\n",
    "    \n",
    "    # plot the prediction line, and the test data\n",
    "    plt.plot(x_vals,prediction, color='k', label=\"Prediction\")\n",
    "    plt.scatter(X_test, y_test, label=\"Test Data\")\n",
    "\n",
    "    # label your plots\n",
    "    plt.ylabel(\"Number of Taxi Pickups\")\n",
    "    plt.xlabel(\"Time of Day (Hours Past Midnight)\")\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO29eZxU5ZXw/z3dFN0NAo2ICo0sojaDIN3YQBNi4hqTuKFijKOOUSckk00xQySZ30STOBFfJzFmmfgzMaPzjqMYNURjDMaAjjGigA2yNSAKSLMICCjQQC/n/ePeaqqaurduVddyq/p8Px8+dN269dxzl6qzPOecR1QVwzAMw4hSkm8BDMMwjHBhisEwDMOIwxSDYRiGEYcpBsMwDCMOUwyGYRhGHD3yLUBXOO6443T48OH5FsMwDKOgWLJkyU5VHej1fkErhuHDh7N48eJ8i2EYhlFQiMhGv/ctlGQYhmHEYYrBMAzDiMMUg2EYhhFHQc8xGIZRPLS0tLB582YOHjyYb1GKhvLycoYMGUIkEknpc6YYDMMIBZs3b6ZPnz4MHz4cEcm3OAWPqrJr1y42b97MiBEjUvqshZIMwwgFBw8eZMCAAaYUMoSIMGDAgLQ8MPMYCpS5DU3cO28NW/Y0M7iygpkXVjO1tirfYhlGlzClkFnSvZ6mGAqQuQ1NfPvp5TS3tAHQtKeZbz+9HMCUg2EYXSZroSQR+Y2IvC8iK2K2HSsifxaRde7//d3tIiI/FZG3ReQtERmfLbmKgXvnrelQClGaW9q4d96aPElkGMVBaWkpNTU1jBkzhquuuooDBw6kPdZLL73ExRdfDMAzzzzD7NmzPffds2cP//Ef/9HxesuWLUybNi3tY3eVbM4xPAx8utO2WcBfVPVU4C/ua4DPAKe6/6YDv8yiXAXPlj3NKW03DCMYFRUVLF26lBUrVtCzZ08eeOCBuPdVlfb29pTHvfTSS5k1a5bn+50Vw+DBg3nyySdTPk6myJpiUNX/BT7otPky4BH370eAqTHb/0sdFgKVIjIoW7IVOoMrK1LabhhG6px11lm8/fbbbNiwgerqav7hH/6BMWPG8N577/HCCy8wefJkxo8fz1VXXcW+ffsA+NOf/sSoUaMYP348Tz/9dMdYDz/8MF/72tcA2L59O5dffjnjxo1j3Lhx/O1vf2PWrFmsX7+empoaZs6cyYYNGxgzZgzgTMrfeOONjB07ltraWhYsWNAx5hVXXMGnP/1pTj31VL71rW9l7NxzPcdwgqpudf/eBpzg/l0FvBez32Z321Y6ISLTcbwKhg4dmj1JQ8zMC6vj5hgAKiKlzLywOo9SGUbmuPXWW1m6dGlGx6ypqeEnP/lJoH1bW1t5/vnn+fSnnaDHunXreOSRR6ivr2fnzp3cddddvPjii/Tu3Zt77rmHH//4x3zrW9/ii1/8IvPnz+eUU07h6quvTjj2N77xDT75yU/yu9/9jra2Nvbt28fs2bNZsWJFxzlv2LChY/9f/OIXiAjLly+nsbGRT33qU6xduxaApUuX0tDQQFlZGdXV1Xz961/npJNO6sJVcshbuqo6i02nvOC0qj6oqnWqWjdwoGdzwKJmam0Vd18xlqrKCgSoqqzg7ivG2sSzYXSR5uZmampqqKurY+jQodx8880ADBs2jPr6egAWLlzIqlWrmDJlCjU1NTzyyCNs3LiRxsZGRowYwamnnoqIcN111yU8xvz58/mnf/onwJnT6Nevn69Mf/3rXzvGGjVqFMOGDetQDOeddx79+vWjvLyc0aNHs3Gjb2+8wOTaY9guIoNUdasbKnrf3d4ExKq5Ie42w4OptVWmCIyiJahln2micwyd6d27d8ffqsoFF1zAY489FrdPpj2cIJSVlXX8XVpaSmtra0bGzbXH8Axwg/v3DcDvY7b/g5udVA/sjQk5GYZhhIb6+npeffVV3n77bQD279/P2rVrGTVqFBs2bGD9+vUARymOKOeddx6//KWTX9PW1sbevXvp06cPH330UcL9zzrrLB599FEA1q5dy6ZNm6iuzm7YOJvpqo8BrwHVIrJZRG4GZgMXiMg64Hz3NcAfgXeAt4FfAV/JllyGYRhdYeDAgTz88MNcc801nHHGGUyePJnGxkbKy8t58MEHueiiixg/fjzHH398ws/ff//9LFiwgLFjx3LmmWeyatUqBgwYwJQpUxgzZgwzZ86M2/8rX/kK7e3tjB07lquvvpqHH344zlPIBuKE+guTuro6tYV6DKM4WL16NX/3d3+XbzGKjkTXVUSWqGqd12esV5JhGIYRhykGwzAMIw5TDIZhGEYcphgMwzCMOEwxGIZhGHGYYjAMwzDiMMVgGIYB7Nq1i5qaGmpqajjxxBOpqqrqeH348OHA4/zmN79h27ZtCd+77rrrGDFiBOPGjeO0007jhhtuYMuWLUnH/PGPf5zTtbBNMRiGYQADBgxg6dKlLF26lC9/+cvMmDGj43XPnj0Dj+OnGADuu+8+li1bRmNjI2PHjuXcc8+lpaXFd0xTDIZhGAGY29DElNnzGTHrOabMns/chuy1V3vkkUeYOHEiNTU1HZXIra2tXH/99YwdO5YxY8bw05/+lDlz5rB06VKuvvrqpJ5GSUkJ//zP/8yxxx7LCy+8AMD06dOpq6vj9NNP5/vf/z7gKJL333+fs846i/PPP99zv0xiS3sahlFw5HJ52xUrVvC73/2Ov/3tb/To0YPp06fz+OOPM3LkSHbu3Mny5c5x9+zZQ2VlJT/72c/4+c9/Tk1NTaDxx48fT2NjIxdddBGzZ8/m2GOPpbW1lXPOOYdp06YxY8YMfvSjH/HKK69QWVkJkHC/0aNHZ+yczWMwDKPgyOXyti+++CKLFi2irq6OmpoaXn75ZdavX88pp5zCmjVr+MY3vsG8efOSts/2IrYt0WOPPcb48eMZP348q1evZtWqVQk/E3S/dDGPwTCMgiOXy9uqKjfddBM/+MEPjnrvrbfe4vnnn+cXv/gFTz31FA8++GDK4y9dupSLLrqIdevWcf/99/PGG29QWVnJddddl3BeIeh+XcE8BsMwCo5cLm97/vnn88QTT7Bz507AyV7atGkTO3bsQFW56qqr+P73v8+bb74J4NtCOxZV5b777mPXrl1ccMEFfPjhh/Tp04e+ffuydetW5s2b17Fv7Jh++2UK8xgMwyg4crm87dixY7njjjs4//zzaW9vJxKJ8MADD1BaWsrNN9+MqiIi3HPPPQDceOON/OM//iMVFRW88cYbR2U0zZgxgzvuuIPm5mYmT57M/PnziUQijB8/ntGjR3es0jZlypSOz0yfPp3zzz+fk046iT//+c+e+2UKa7ttGEYoSLXt9tyGJu6dt4Yte5oZXFnBzAurbVXDBKTTdts8BsMwChJb3jZ72ByDYRiGEYcpBsMwQkMhh7bDSLrX0xSDYRihoLy8nF27dplyyBCqyq5duygvL0/5szbHYBhGKBgyZAibN29mx44d+RalaCgvL2fIkCEpf84Ug2EYoSASiTBixIh8i2FgiiEtLE3OMIxixhRDiuSyeZdhGEY+sMnnFMll8y7DMIx8YIohRXLZvMswDCMfmGJIkVw27zIMw8gHphhSZOaF1VRESuO2Zat5l2EYRj6wyecUiU4wW1aSYRjFiimGNLDmXYZhFDMWSjIMwzDiyIvHICIzgH8EFFgO3AgMAh4HBgBLgOtV9XA+5DMMw8g1YSqczbnHICJVwDeAOlUdA5QCnwfuAe5T1VOA3cDNuZbNMAwjH0QLZ5v2NKMcKZyd29CUF3nyFUrqAVSISA+gF7AVOBd40n3/EWBqnmQzDMPIKWErnE2qGESkQkTE/XukiHzW/UFPC1VtAv4d2ISjEPbihI72qGqru9tmIKEPJSLTRWSxiCy2LoyGYRQDYSucDeIxvIJj3Q8C5gNfBH6T7gFFpD9wGTACGAz0Bj4d9POq+qCq1qlq3cCBA9MVwzAMIzSErXA2iGIoUdUDwJXAL1X1cuCMLhzzfOBdVd2hqi3A08AUoDLGExkC5Ce4ZhiGkWPCVjgbSDGIyATgWuAP7rZSn/2TsQmoF5FebojqPGAVsACY5u5zA/D7LhzDMAyjYJhaW8XdV4ylqrICAaoqK7j7irF5y0oKMldwG/A94A+qukJETsYJL6WFqr4uIk8CbwKtQAPwIPAc8LiI3OVueyjdYxiGYRQaYSqclaDrq4pIL0BVNTRtROvq6nTx4sX5FsMwDKOgEJElqlrn9X6QrKTxItIArAXeFpElIlKbSSENwzCM8BAklPSfwK2qugBARM4GHgbGZU8sI1uEqbrSMIxwEkQxtEeVAoCqviQi7VmUycgStiypYRhBCKIYXhKRXwCP4fQ2uhqYLyJnAKjqW1mUz8ggftWVphgMIzNkwyvPtacfRDFEJyg61y5MxFEUn8ioREbWCFt1pWEUG9nwyvPh6SdVDKp6VlaOHHKKMRY/uLKCpgRKwJYlNYJQDJZwtsmGV54PTz+pYhCR7yTarqo/zLw44aBYY/EzL6yOOy+wZUmNYBSLJZxtsuGV58PTD1L53BbzL4LT9fTUrEkUAsLW6TBThK260igcsvGdyOb3bG5DE1Nmz2fErOeYMnt+ztpXZ6PnUT76KAUJJd0T+1pE7gH+lDWJQkAxx+LDVF1pFA6FZAnn0xPJhleeD08/nfUYynCa3BUtYet0aBj5ppAs4Xx6/NnwyvPh6QeZY2jAyT4Cp3neIODurEkUAiwWbxjxFJIlnG+PPxteea49/SDpqtNi/m4FtqnqoSzJEwqiN6CYsiUMoytk4zuRre+ZZd91naRN9ETkC6r6cKdtd6nq/5dNwYJgTfQMw+hM5zkGcDwRS7Q4QrImekE8hmtEpFlV57gD/hTolykBDcMwMol5/F0niGK4HHjW7Y/0GaBZVW/IrliGYRjpY9l3XcNTMYhI35iXNwDPAq8C3xGRvqr6YbaFMwzDKDYKodrbz2NYiZONJDH/X+b+U2Bo1qUzDMMoIgql2ttTMajqSbkUxDAMo9gplA7HQVZw+7KIVMa87i8i07MrlmEYRvGR7xqLoASpfP6yqu6JvlDV3cA/ZU+k4iVf/VsMwwgHhdJVIYhiKI19ISIlOM30jBSIxhab9jSjHIktmnIwjO7DzAurqYjE/aSGsqtCkHTVP4vIY8AD7usvAy9mT6TipFBii4ZRbOQrC8jvuIWclRRlJvAVYIb7+s/A/581iYqUQoktGkYxka8soGTHDZsi6EyQttttwM/cf0aaWP8Wo7uTCcs91TFS8dQz6Vl4HfebTywDwpWamgjPOQY3fISINIjIm53/5U7E4qBQYouGkQ0yMceWzhhBPfVMzwF6HbdNtSDmFv08hpnu/9N89jECUiixRcPIBnc+s7LLc2zpzNN5eeqVvSJMmT2/47t44HBrRucAvY7b1XFzhV+B22YRuRg4BViuqn/JnVjFSSHEFg0j08xtaGJPc0vC91KZY0tnni7Rmg+RUmHfwVZ2H3Bk8voBT1W+ZMfNxLi5wq9X0s+AWuA14HoReUpVf5gzyQqUQuiDYhi5xG/ltFTm2NKZp0vkqe8/1OqpqLoiX6LjfvOJZbQlWNog7HOLfqGkc4AaVW0Vkd7Ay4ApBh8KpQ+KYeSSZBZ9UNJd8a2zpz5i1nOBjtfVOcDoMQtxNUg/xXBYVVsBVHW/W9hm+GC1CoZxNF6Wfv9ekZS+F5map/Ocd6iI0LusR5e9/c5RgyvPrGJB446CiiL4KYZRMdlHAlS7rwVQVR2f7kHd3ku/BsbgdGq9CVgDzAGGAxuAz7ntNwoGq1UwipGuhke9LP07Ljk9ZVkyMU/nJc+dl57e5bETRQ2eWtIUePW4sISi/RTD2Cwe937gT6o6TUR6Ar2A7wB/UdXZIjILmAXcnkUZMo7VKhjFRibCo2HLyMumPF2JGoQpFO2XlbQ+GwcUkX7AJ4AvuMc5DBwWkcuAs93dHgFeosAUQ7oxUMMIK139octlQVsq+2YrQ7ArUYMwhaKDtMTINCOAHcB/isg4YAlwC3CCqm5199kGnJDow27L7+kAQ4eGa62gsFlGhtFV0v2hy4T1m8oYYbG2uxI1CFMoOh8Tyj2A8cAvVbUW2I8TNupAVRVn7uEoVPVBVa1T1bqBAwdmXdhUmVpbxauzzuXd2Rfx6qxzTSkYBU26baL9rN+gpDJGJo6XCbrS4SBMLbnzoRg2A5tV9XX39ZM4imK7iAwCcP9/Pw+yGYYRQ7o/dJmwflMZIyzW9tTaKu6+YixVlRUIUFVZEXjiOUxtc/wK3B5T1WtEpIF4671LWUmquk1E3hORalVdA5wHrHL/3QDMdv//fTrjG4aROdINj6YSUuk8N3DOqIEsaNyROGSQYIy5DU2UiISmkCzd+YswhaJFE1xMABEZ4rbFGJno/a5MTotIDU66ak/gHeBGHO/lCWAosBEnXfUDv3Hq6up08eLF6YphGEaW6BzzB8f67Ww9J9rPj85j+H0+0fEMBxFZoqp1Xu/79kqK7uNa9rGDngWkrRhUdSmQSKjz0h3TMIzwENT6TTQ34EWpCFeeGW+Ne32+VMSUQhcIkpX0tIg8pKo/FpEynFDPx4BJ2RXNMIxCJkhIJZU5gDZVnlrSRN2wYzvG9fp8u6ophS4QZPJ5EnCaiPwVWAR8gKMYuiVzG5qYMns+I2Y9x5TZ80PfV90wwoDX9ybVOYDOmUZhyuQpJoIohoPAbqAf0BtY7a7q1u3I9GIehtEd8PveJMrESUaslxCmTJ5iIkgoaRHwPHAmMBB4UESmqernsypZCMl0ZWJY+qIYRiIy9Xz6fW9enXVuxz6ds5K81kmI9QbClMnjR6F914Mohi/H1Bw0AReJyI1ZlCm0ZDJXOiyVmoaRiEw+n8m+N15zEV6ZTZ29gbAvgFWI3/WkoaSoUhCRY0VksIgMBuZlXbIQksl4ZlgqNQ0jEZl8Pr2+HyUivnN10WKx/r0iHdvKehRe9/9C/K4nvcoi8lkRWYtTsfy6+//8bAsWRjIZzwxLpaZhJCKTz+c5oxK3rmlTDTRXd7ClvePvPc0tBTevV4jf9SDq94fAFGCNqp4EXAi8klWpQkpXyt07Y9kUuccyyoKTyedzQeOOpPuEvQdSVyjE73oQxdCqqjuAEhERVf0zMDHLcoWWTDXJs2yK3GIZZamRC+84yH6FaG13phC/60Emn/eKyDHAX4H/EpH3gcK5KyGlULIpioUw9bovBLyeT4Aps+dnpG9Sov2CfjaZtR2mLKBC/K579krq2EGkD3AAx7v4B5x6hv/rehF5xXolGUEZMeu5hE3ZBHh39kW5FqcgCdr/KMjnOuM1TjrHTFfO7kSyXkmeoSQReQFAVT9S1TZVbVHVh1T1x2FQCoYRlGj3zUSEOc4bNtKN9yeam7uufmigubrOn62siFAeKWHGnKWe80TFMC+Rb/xCSeFbBccwUiRqPSZqyRz2OG/Y6Eq8vyu1BtHPBq0HKIZ5ic7s2rWL5uZmhgwZkpPj+SmGfiJyhdebqvp0FuQpesIU++wOdOfum5l+1rqybGUmCDpP1K8iwp7mlqM+XyjeYUtLC8uXL2fhwoUd/9atW8dNN93EQw89lBMZfBUDcDFOGLYzCphiSJFCrIAsdLpr981sPGszL6wOVImcLYJ4AnMbmth/uPWofSIlElrvcOvWrXFKYNGiRTQ3O+d0wgknMHnyZG6++WbOPffcnMnkpxg2qupNOZOkG3DnMystM6aLpGoF59vKzRfZyMLqnF1T2SuCKsyYs5R7563JmPfrdY+D3Mt7562hpe3osOEx5T06wlH59NgPHjxIQ0NDnCLYtGkTAJFIhPHjxzN9+nTq6+upr69n2LBhiMf8WDbxUwy5l6aImdvQlNC9hcKOfeaSdKzgfFu5+SJbcfZU4/2p4jdukHvpdX57DrTk3GNXVTZs2BCnBBoaGmhpcX4Hhg0bxuTJk5kxYwb19fXU1NRQXl6ecTnSwU8xXJ8zKboBfhkRxW69Zop0rOBCzCHPBNn2lFK5F6lY6al2Yu08lt95p/P8+Mk+t6GJO59Z2WHw9S1t5aphh+ixa32HInj//fcB6NWrFxMmTOC2226jvr6eSZMmMWjQoITHDAN+S3uuyKUgxY6fpVbs1mumSNcKDnv3zWyQbU8p6L1I1UpPtxNrFL/znjFnaUrn4id7e3s7Mx78I/s3N3KoqZFDWxrZuHMTy9Xp61RdXc1nPvOZjpDQmDFj6NEjSD1xOCgcSQuUqMXhVUbYv1ek2/1opUt3nS9Il/JIScePWmVFhDsvPT1jz1rQe+FlpX/v2ZUJZenqPfbzEO+dtyalsWNlb2v+iMNb17JnSyM3znmbj95bTdvBfQBIWW/KBlfT67TJlA0eRcXgambf8PGC/l57KgYR+Yuqnici96jq7bkUqlhIVvFZESnljktOz7FUhUt3nS9IlUTP3aHWdp9PpE7Qe+Flje92Y/6dfzwzcY+9vIqgY7e2trJixQoaFzzFoS1rOLRlDa0fbHbelBIixw2lovrjlA2upqxqFD2OrUIkvla40LMN/TyGQSLyMeBSEXmcTpPRqvpmViUrArxy6MHJo7/yzO4X4ugK3XW+IFVy0Rcq6L3w65OUSJ5s3uOptVUs3vgBj73+Hm2qHd/B+kGlzJ07Ny5d9MCBAwCU9KqkbHA1x4w5l7LBo+h54imcdMIAAN/+T4WebejZK0lEpgE3Ax8HOjckUlXNXVKtB2HvleTVnyeK9W8xskGY+kLNbWjiVo/Yfq7lmdvQxO1PLOHDzetcT6CRw1vX0LrXmSCORCLU1tZ2zAt82Gc49y3cw8EYbyv6nQWY+dtltLR7f8PD3IcrWa8kv8nnJ4EnReRfVfUHWZGuyEnWVbLQrYpck+8c9HwTe/79KiKIOGmYna9FtudivO6D1/bYzJ2uyJPq/VdVNm3axMKFC3nttdd46OkX2LdlHbQ5BXClfQdSNngUQz5+Jf/zL9dTW1t7VLroCVX+x/Q6t3TOL0wk7a4KICKXAp9wX76kqn/IqlQBCbvHEKSrZJitijDR3TtmBpmvil6LbF4rr7GvPLOKp5Y0JTwm0GV5gpzT/v37Wbx4cVzdwLZt25x9KypoH3AyZVWjKBtcTc9B1fToM6BjrA1d+A4W4rOZtscQM8DdOAvzPOpuukVEPqaq38mQjEVLbLzUy3PoilXRnSzoTMTNs3W9cnEf/OarIP5aZDNO73UfonH7RDIFqT9I9biq7Xy4bTMzf7iAecd9yMKFC1m+fDltbc4+p556KhdccEFHWGjs2LGM+u6fEzZTLE1SWZzs/kb/jvUeyiP+a6AlGzOod5gtgqSrXgTUqDoJuiLyCNAAmGIIgFelKHQto6a79V3qaiVvPip1M3kfgpxn7D7Zqt3wkiPRD27s/l2V571tOzoyhA5tWcPhrWtod9NF3+/bl0mTJvGd73yno3hswIABR43hJaPXdkjt/sZmfu0+0OK5X7IxO78fG6rK1fc8aB1DJfCB+3e/LMlS1GTaiutuK5J5xc0VZ0WxZNcyyPVKx/LP1X0IsgpaLlY18+pcWiqS8Ac2KlMqx25ra2PFihVxIaH3Ghvdd4XIwGH0qp5C2eBqTqoex5IffYGSkuSrFFd5XMOqBNctKm+i/RPd31Seg2T7puIdZosgiuFuoEFEFuCExD8BzMqaREVMJq24Yuw570eiHPQoQaworx/V6PZ0Lf9c3Qe/84fk3mcmPBu/zqVXTzwp4RzDzAurkx57+/btvP766x1K4I033mD//v0AHHfccUyePJkzz7uM/91TCQNHUlLWq2P8O64YG0gpQPA6hiBzg53vbyrPQbJ9U/UOs0FSxaCqj4nIS8AEd9Ptqrotq1IVMZmKR3eXKuDY61XZK0JZj5KEFmsyK8rLoo3Gl9O1/DN9H/yej9hK5opICeWRUs+4c+dx9h9q7bJn49e59K6pY6kbdmxC2afMnt9xbG1r4fD2d/hwyxq++Oy93Pbhu7z77rsA9OjRg5qaGm666aaOuYERI0Z0dBft6ncnqNeezGKHo+9vKs9Bsn0z4R12lUChJFXdCjyTyQOLSClOfUSTql4sIiOAx4EBwBLgelU9nMlj5ptMxqO7QxVw5+u1+0ALFZFSz/39rKhk8eV0Lf9M3gev52Pxxg+OssZBuOOSxC0uEo3jRSqWp1/nUjjaI46mi6577QUObmnk8JY1HNq+Htqc/Uv7HMfZF57NV7/6Verr6xk/fjwVFd4/eJnwuIOMkexHOdH9TeU5SLZvV73DTJDPXkm3AKuBvu7re4D7VPVxEXkAp7jul/kSLlWCWDOZjEcXYxVw52t44HBiK9cvnh0bG47uV1VZQf9eEXYfONrTqEpipSWzzDJ5H9LJ+El0nCAWb5TYa5ZM/mTXaP/+/SxZsiRubmDr1q0ASI+e9DzxFPqeeQk9B1dTNqiaHn2PY0tlBSPPrWZKiJ5br+cLnOcl0fVJ5TlItm/n9/ORlRSojiHjBxUZAjwC/BtwG3AJsAM4UVVbRWQycKeqXug3TljqGILmMYepIjVsBInrxlIRKQ2USx8lUiIgxIVCcpX7H5Ths55LaX+v5yZZxX2UZPUHftk0qkrr7i3o9rWMK9vB9reX89Zbb3Wki55yyikd4aCDlSfzq5VtHGxPnBYatpx/v/vQlXqHMNGlOgY33LNSVUdlWK6fAN8C+rivBwB7VDU6s7UZCMdTEoCgnkAm4tHFWruQipUbtdo6Xwe/MVraFcHpZpvI8krH8s/0vfCzVBPh9dwEiVFH17wO+uzu3buX3jtXUbvrLzw//xU+3LiK9oMfAbCvTx8mTZrEt7/9berr65k4cSIDBw6MG3Nkilk++SSV7KWgFNr31lcxqGqbiKwRkaGquikTBxSRi4H3VXWJiJydxuenA9MBhg4dmgmRukzQ+HRX49HFXLsQNNYdvV6JYsVe/fajKHCwpZ37rq7xdPHTrcTNxL3wUwqJPCSv5yZZjBqOrHmd6Jppexsb1q3mV786suDM6tWrUVVEhNGjR3PltZ+jvr6eyZMnM2rUKEpLved+4Mi19fJmwpRNl+n5u0L83gaZY+gPrBSRN4D90Y2qemmax+80k/MAACAASURBVJyC07H1s0A5zhzD/UCliPRwvYYhQFOiD6vqg8CD4ISS0pQhowT1BLoajy7m2gWva1hZEaF3WY9A1yuIpZyp65WNe+FnqSbykLyOk0rF/eDKCjY1beXQVrdwbEsjh7auQw83M/03MGDAAOrr67nmmmuor69nwoQJ9OuXfilTIWTTWc1RMMXwr5k8oKp+G/g2gOsx/LOqXisivwWm4WQm3QD8PpPHzSapWBhdyawo5toFr2t48bhBLGjckfYYiUh0vbwmrb1+ELyuedOe5rTDBueMGsijCzfFWdR+HpIfXhX32taCfLCR0foh1177EOv/91W2b97ofKiklJ7Hj6DfGedzw2Xn87XPf5aRI0dmdDH6bFjj2WrRne2ao6Y9zYyY9VwoQ0tB6hheFpFhwKmq+qKI9AL8/cb0uB14XETuwmm58VAWjpEVcpUhVAjWVrokuobnjBoYNzGazAUPYinD0der849nNKTjdzw/7yS2HXPQsMHchiaeWtIUpxQEurxmx5nHtXNxn408+uyLfPDuSg5vfxttbeHXwODBg/n45MkcM+16ljQfx4e9hzBkYP+s/khl8rtSKCEav2dFCafcSbOSROSLODH9Y1V1pIicCjygquflQkA/wpKVlCvSyZwptEmvWKbMnu8ZWok2ZvMj6PXyOo7f8eY2NDFjztJA2T/JZJ7b0MQ3n1iWcI6hf68IDd/9VKBjHDhwgDfffDMuXbSpyYnIlpWVUVdX15EpVF9fz5AhQwJKH066+nzkiqAZdwKe81+ZpsvdVYGv4nRXfR1AVdeJyPEZks9IgVStrUKxqLzoaugs6PVKNi+R6HhTa6s8F6AJOgYcuUdeE89eS2CqKuvXr49TAsuWLaO11UnsO/nkk/nkJz/ZoQTGjRtHz549A8tbCBRKaLXzc+hlTCiOtxn7mXwRRDEcUtXD0TijiPSAwIaSkWFSiX3e+czKgpv0giNejtdDlkroLMj1SpYmWiKS8MfZa7I4EbEyx3pxJQFSVL/5xDJu+a9XOeajjdSU7WDvxlUsXLiQXbt2AXDMMccwceJEvvWtb1FfX8+uXkP51aJd/G1PMxv2V3BSj8FMKDKlAIUVWo19Dv081JZ25XvPrsz79zOIYnhZRL4DVIjIBcBXgGezK5bRVeY2NHmuLBU2iyqWIAvSZLodQLIf5jbVhJ5WoonUSKmAErfkY6zMXvMZsWh7Gy273nOzhJwlKFt2vgcoq4CTTj6Nyy67rMMbGD16dEe66NyGJu4pYC8xFQq1LczMC6t9vU0vLzGXBFEMs3DaUywHvgT8Efh1NoUyus6989Z4vhdGiyqKX5GaX5aQF0EWRAlSWJbI0/IKVSXaFrtv5/NrO7CXQ1vXcripkUNbGjm0dS162FHeJeV9KBtcTa9RZ1E2eBRlg07lpBMH8pBHDL0QUyPTpVDbwkytreJ7z65M2KIlSr7vV5CspHZ3cZ7XcUJIazQffTSMlPDzCsJsUXnJLZDyhGLQBVGCVht7zTX4ZUh1pmnXRxzascGpF3C9gdbdTj8hpISex4+g9+nnUja4mrLB1fToP/iodFG/e5usvXixkcm00lxyxyWnM/PJZQm71UL+vfogS3teBDwArMf5fo4QkS+p6vPZFi4MFGpWj1f8tX+vSKjlz2TcOJn1nGgOJplsQYh9ZgbIPs4bsJc1b73J/FdeZd/mtWir0zS4tHd/elaN4phxF1I2eBQ9TzyFkkh5ktH95fDyfkSc2HahPcfFSvTaz3hiKYnsknx79UFCST8CzlHVtwFEZCTwHFD0iqGQs3q84q93XHJ6HqVKTibjxn5ZK35zMNFjpipDc3MzP3l8Hvf9z3Psf281h7asYcNHO1kCUNqDshNO4ZiazzghoapqSvsMTLl4LJkcXt6P6hGvoZCe42Imeu3DOE8SRDF8FFUKLu8AH2VJnlBRyPHaQo6/Qmbk9vM+/OZggrSgUFXeeeeduHTRpUuXdqSL9uh3AmVDTu8ICfU8/mSkRyTlcwDHC2hXDXQtgmZKFcpznGtyHSEI6/fUs8BNRK5w/7wAGAY8gTPHcBWwSVW/khMJfch2gZu1yS5s/Arc/IrTfpKgyOijjz5i0aJFvPbaax2KYOfOnQD07t2biRMnUl9fz0NretBz8GmU9u6fkXOIlAr3ThuXdnM/P+w5jicMrddzRVcK3C6J+Xs78En37x1AeNNaMkgh5EkX6hxILvCzxrzaZvTvFaG9vZ3aGb/hvTXL6LFzPeW717Np/RqiRlTV8FOQYWcy4MxTGVJ9Bv963ae4csIwAF5OUkWdCv17RXxXaUt0XonOef+h1oRhszA9x2GgkCMEmSYvC/Vkimx7DGG3IMIuX5iJvXZtzR9yeMsa2ratpe9HG9jYuIz2QwcAKCk/hoqqUVxy/if4wuWfYkf5EP7txfc8r3mqCw4lIlIi3HuVt5eQ6n235yQY3SlC0OWWGO5azF8Hhsfu34W226HCz+IOa/wvilk4qdPS0sLy5cvZsnAhJzbMZ+HChRzc5fQTKikp4eAJJ9Pr7z7pTBAPrqbHsYMRKWFTZQXNx1dzR4KeRrHXvHMjv1QX3wE4pryH7/3zuu/ffGIZM+YsLbjnOCwUQoQgVwSZfJ6L0+n0WaA9u+LkliBZR2HOky6UXjH5ZMuWLXETxIsXL6a52bk+J5xwAp8+azL19V+nvr6eM888kzF3vZxwnOiz4fUjH3vNvbJNgrLHp/Cp87Fi8esKG+bnOCwUaiV1NgiiGA6q6k+zLkkeKHSL2yyceA4ePEhDQ0OcIti0yVl4sGfPnowfP54vfelLHa0khg4delS6qJ+F7/cj3/map7JUaWf6VfhnL+VyQaLuRCqeVbHP7QVRDPeLyB3AC8Ch6EZVfTNrUuWIQre4u7OFo6ps2LAhTgk0NDTQ0uJY28OGDeNjH/sYt912G/X19dTU1FBWVpZ03FTDPpD4mnflGdp/uNW3V05XFiQy/AniWRVyfVNQgiiGscD1wLkcCSWp+7qgKRSLO5UMlGKzXKLs27ePRYsWxSmC999/H4BevXoxYcKEDiUwadIkBg0alNZxvOoAvDyJUpG4iedUuqZ60dKmvtZ+5/vudaywPcfFQj4iDbn2UIIohquAk1X1cNakyBOFYHEns06KMXbc3t7OmjVr4pTAihUraG937JLq6mo+85nPdISExowZQ48eQR7l5Hg9E1eeWRW3mlx0u1c2UrpKIUoyaz/2vntlHYXpOS4mch1pyIeHEuTbtAKoBN7PigR5pBAs7kKfBwnCBx98wOuvv96hBF5//XX27t0LQGVlJZMmTeLyyy+nvr6eiRMncuyxx2ZNFr9nom7YsSl1TYUjnoYICXvieJHqmhNeMhuZJ9eRhnz8BgRRDJVAo4gsIn6OoSjSVcNucRf6PEhnWltbWbFiRZw3sGaN056ipKSEsWPH8vnPf77DGzjttNMoKSnJqYx+HVO9nhWv+9Guyk+urkkpQykdaz/sz3ExketIQz5+A4IohjuydnQjKYUyD+LFtm3b4pTAokWLOHDAKR47/vjjqa+v54YbbmDy5MnU1dVxzDHH5Fni9EjWlymoUoidszDCSToeWuc5gnNGDWRB445An8/Hb4BVPoecQqpaPXTo0FHpohs3bgQgEolQW1sbtxj98OHDU+4uGlbS7cvUmWKssu3uBKmGz3XleiYqnz/iyBrPPYEIsF9V+6YlkZESYY0fqyobN248Kl308GEnR2Ho0KHU19dzyy23UF9fT21tLeXlidcaKOSc8FjZ+1VEKI+UsOdAS6C+TIkIYgV2xfpMNlYhXftCIYjH6DdnkI/fgJQ8BnHMu8uAelWdlTWpAtIdPIawsG/fPhYvXhynCLZv3w5ARUUFEyZM6PAEJk2axODBgwONW0geUWeCyp5ov0iJgBC3gleQ8+6q9ZmO/EbX8OrB1Jlceotd9hhicZf0nOsWvOVdMWST7mxJtbe3s27dujgl8NZbb3Wki5522mlceOGFcemikYh3te7chibufGZlR4fP2K6hhZZ1laxWIVNrQ3uRqvXp9xwX2rUvVIJUqkf3CwtBQklXxLwsAeqAg1mTKAR0h8rGWHbv3s0bb7wRly66e/duAPr27cukSZP4l3/5lw5vYMCAAYHHntvQxMzfLqOl/cgP6O4DLcx8chlQWFlXQWsVMrE2tBdBQ1LRVer8nuNCuvaFTJBK9bDVnQTxGGLXZWgFNuCEk4qKIJbg955dWfBeRGtrKytXrozzBhobGwEQEcaMGcO0adM6vIFRo0Z1KV303nlr4pRClGh1byFlXQXNLkpV9lS806DdWr2yoWI9gkK69oVMIo+xK/NCuSCpYlDVG3MhSD4JagnuPtDCbrfzZaF4Edu3bz8qXXT//v0ADBw4kPr6eq6//nrq6+uZMGECffr0yejx/azPLXuauS9Bjn/YrKcoQSzpVGVP1TsNohSiMsyYszTh+9HzKITK/2Kh0OpMPBWDiHzX53Oqqj/Igjx5Id1OmEHisbmaq5jb0MQ9z61g49qVlO9+h6qWzWxes4wNGzYA0KNHD2pra7nppps6vIERI0ZkPV3UL746uLIitFlXifA6l1TWZO5Mqmsr+PVy6iyDVzZU1CMopGtv5BY/j2F/gm29gZuBAUDRKIauxFT9PpvNuQpV5b333mPhwoX89zMv8uLLf6V569vQ5ng07/UdyMRJk/j617/ekS5aUZH7EMHMC6uPmmMAZy3jqGVaKNbUOaMG8t8LNx21/ZpJJ3HX1LFpjZnq2gpeVn6iTKIgHkGhXHsjt3gqBlX9UfRvEekD3ALcCDwO/Mjrc4VIsqyBqjTXzc1k1sf+/ftZsmRJXFho69atAJT06EnkxFPoe+Yl9Bxc7aw81uc4qKzgttvy2wQ3ep5eWUmFxILGHSltD0LQtRW++YQzWZ+KlW8egZEuvnMMInIscBtwLfAIMF5Vd+dCsFwy88JqbvWIxwrw6qxz0+pgmW7Wh6omTBdta3OOfcopp3Deeed1hISufHwzlB59K8OSXVIsVmk2sniCrq3QpnpUV90gFMu1N3KL3xzDvcAVwIPAWFXdl4kDishJwH8BJ+BUVD+oqve7SmgOztrSG4DP5UoJTa2tirNoY+lKPDZo1seePXuOShf94IMPAOjTpw+TJk3i29/+dke66HHHHRf3+ao/77XskhyQbhZPKuuK+63jUAg1Bt25/qeY8Kx8FpF2nG6qrRBXuCc4k89ptcQQkUHAIFV90w1RLQGmAl8APlDV2SIyC+ivqrf7jZXJyudsVIEmGrO8FP5pXBm99r7La6+9xsKFC1m9ejXgpIuefvrpcf2ERo0aRWlpac5lN44mneuc6meSVTaHuZeSPYeFQ9qVz6qalV7HqroV2Or+/ZGIrAaqcGojznZ3ewR4CfBVDJkkW/HY0kMfcuCdFRzasob2bWtp3b6OWw848/oDBgygvr6ea6+9tiNdtG/f1PWtxZJzQzrXOdV5pui2bz6xrOBWZbNK6uIhM8tepYmIDAdqgdeBE1ylAbANJ9SU6DPTgengNGrLJF2Nxx4+fJhly5axcOFCnnx+Pq+9tpCWPducN0tKKT/hZM65eBp/f8n51NfXM3LkyIyli1osOTekep3TmZeIjl9oNQZWSV085E0xiMgxwFPArar6YewPpKqqiCSMcanqgzjzHtTV1eWtZ7iqsnnz5rgJ4iVLlnDokLOWUc++A4icWE3vms9SVlVNzxNGUhIpZ29lBddddyRTqJhissV0Lpki3XmJQvQCrZK6eMiLYhCRCI5SeFRVn3Y3bxeRQaq61Z2HCNVSogcOHDgqXXTLli0AlJWVUVdXx9e+9rWOuYGzfr4sYUfFWOupmHoyFdO5ZJKuVBcXmhdoldTFQ84Vg9u6+yFgtar+OOatZ4AbgNnu/7/PtWxRVJX169fHKYFly5bR2toKwMiRIznnnHM6lMAZZ5xBz54948YYXLk2qfVUTDFZr3P53rMrC+5cMkkhWv7p0p3OtdjJh8cwBbgeWC4i0eKB7+AohCdE5GZgI/C5XAm0d+9eFi1aFKcIdu3aBcAxxxzDxIkTuf322zvSRQcOHJh0zCDWUzHFZL1k3n2ghbkNTd36x6HQLP+u0J3OtZjJuWJQ1b/iZN0l4rxcyXHw4EEuueZGXn1tIc3bNxLNyB09ejSXXXZZhzcwevTopOmiiQhiPRVTTNavgtfLA7I5CcMIJ3nNSsonz6/ayf+++hqlfQbS7+NTKBs8in5DR/Fvfz85Yz9OyaynYorJ+lWPJ/ImbE7CMMJLSkt7ho2uFLhNmT0/oYVbVVnBq7Ny118om1Zzri3ymu+9kLB6vH+vCA3f/VTcNq/r37lLKBRPzNo8JCMsZHRpz2Ii0/H9dL/02YrJ5sMiv/PS05n55LK4dYwB9h1sPWqeIUhX0ZlPLgOlozNrmL2KZPc/nfsR9JkyhWNkmqxUNxcCXnH8dOL70S99055mlCNf+rkNTV2UMn38Mp6yxdTaKnr3PNrWaGnXo44b5Dq3tOlR7bqzfQ7pEOT+p3o/gj5TYXz2jMKn2yqGmRdWUxGJn1RON76fjx/hZKTqEc1taGLK7PmMmPUcU2bPT/uHZW+CUFKi4ya6/kEJW9ZWkPuf6v0I+kyF8dkzCp9uqxim1lZx9xVjqaqsQHDmFtJt9hXGtNNUPKJMWp2VvSKBtk+treLKM6s809P8CFvWVpD7n6qHGvSZCuOzZxQ+3XaOATIX3w+SdhqNAzftae5Y0L2qC/HgZHHlVDKeMllod9CjK2iiHIcFjTsSVodHiZRK3BxDlAOHj56zyCdB7n+qGWhBU5mLKeXZCA/d1mPIJMnCUrEWORy9bGOqlnkQCz8VjyhTVufchiaaW9oTvpcoxOQ3flVlBfdOG8e9V42jsiLe29h9oCVUcfQgYclUPdSgoc5MhkQNI0q39hgyRbJitkQWeZR0LPOgFn5QjyhTVqdfXDvRWF7H7ZwyfO+8NUelwYapdUjQVhCprrwWdMwg+xlGKphiyBB+X/pklneqlnmm48qZKrTzO36isYIetxDi6NlIOw46prWhMDKNKYYMMbehie89u5LdBxzLtrIiwp2XOgveV/aKdGxPRIkII2Y9FzhX3Wv5xxKRtGLvsVZndA4kNrMlyHhzG5o85erfK5LSYvXgFMBFt3ldP4ujG4mwuo6uY4ohA8xtaDqqsGtPcwszf7uMxRs/YN/BVt/Pd55zgPgf487FUV5rAndeMD4VEi0OE7SgLCpfIrkqIqXcccnpvsdNVggWKREipRJ3fS2ObiTCWq1kBpt8zgD3zltzVLUvONk0j73+3lFZNVFKE6ze1jkHfW5DE998YpnnHEWyz6dCOkVYU2bP59Y5SxPKVyKknAKcSIaWdqV3zx5JJ24zVYthFC5W15EZzGPIAH6xbi/rXoB2j/ei4/lZ4unKk87ngjTBS4SHPkxLhr3NLSy941MJ30skj1mK3ZNCmI8qBEwxZAC/ltOlHnF39Xmvn5ue6ZfNlGzs4bOeS1on0TkW268ikrAJ3uDKiqP2PXC4NZAXc+czK1OK96abIVVMix4Z6WN1HZnBQkkZYOaF1U4xViciJcI1k07ybP3g5Qnsdwu4/Kycikip79jgXyeRqBZi/+FWIiXx51ERKeWcUQOP2tdvMj2WPc0tKVVUp5uXb5aiAVbXkSlMMWSAqbVV3DttHP1j2j5UVkS496px3DV1bEdhU1Ba2pQ7n1lJSYI5CHA8hbuvGNsxdqK5iijNLW3cOmfpUTH3hLH8Tk3rKisi3H3FWBY07gg8x5GMZPHedFuVZLIpolG4ZLLVTXem267HkA+Gz3ouI+NcVz+Uu6aO7Xg9YtZzvq0lolRESju+JEE+EykR7r1qHDPmLA00flAEeHf2RRkcMfGcR+z5GoZxhGTrMZjHkEP8LPtU+MOyrXGvg1rFsdZ6oLbXbrvsTFvd2bDizVI0jMxhiiGHpJpd5MWe5pa4sFAqLayjMfegn9myp7lLLbI7k81479TaKl6ddS7vzr6IV2eda0rBMNLEFEMO8Zpn6N8rEmfp9vdoXR3LN59Y1qEcYq1l8PdMotZ6tO11Mi+mRIQZc5ZSHik5qpldqpgVbxiFgaWr5hCv3kB3XHK6b/VvIjpXOSerII4eK7bj61NLmpJ6MdH3dx9ooSJSynX1Q3lqSVNKk9EW6zeMwsI8hhzSOQ5eWRGhPFLCjE5ZQ509AC/8MnySxdy/9+xK3x/3RI5Ec0sbCxp3xI3bv1fE15OIZlDlQilY5bNhZAbLSsoTQbNoknkP6WT4zG1o4tY5S33H88ug2pDgeF5ZTtnIQEqEZSUZRnAsKymkBO3pErX8veYC0snwCbJugt/cQyJLPN91BNYjxzAyhymGPJFKpe7U2ip+9LlxGavoDLJugt/cQ6Lq5XxXnFrls2FkDpt8zhOp9nTxW6krSP/5IOs5VFYcWTehyqf/k9dqcV7y5QLrkWMYmcMUQ55IZ9W0RCt1BekqGmQ9h4pIKXdeemTdhETyxeLl2eQrnp+pVegMwzDFkDcyZWEH6Srq1aW1VIR21YTHjv79zSeWJVQkYbPE8+2xGEYmyfcqdKYY8kgmLOwgsXWvfdpVfTOGEq3qBuG1xG3tY6MYCMPaIqGafBaRT4vIGhF5W0Rm5VueQiBINlBXMoasB5Fh5JYwZNiFxmMQkVLgF8AFwGZgkYg8o6qr8itZuAkSW+9q/N0sccPIHWHIsAuNYgAmAm+r6jsAIvI4cBlgisGHILF1i78bRuEQhgy7MCmGKuC9mNebgUmddxKR6cB0gKFDh+ZGspATxKI3q98wCoMwZNiFao4hCKr6oKrWqWrdwIED8y2OYRhGRgnDvF6YPIYm4KSY10PcbYZhGN2KfHv4YfIYFgGnisgIEekJfB54Js8yGYZhdDtC4zGoaquIfA2YB5QCv1HVlXkWyzAMo9sRGsUAoKp/BP6YbzkMwzC6M2EKJRmGYRghwBSDYRiGEYcpBsMwDCOOgl7aU0R2ABu7OMxxwM4MiFNodNfzhu577t31vMHOvfO5D1NVz0KwglYMmUBEFvutfVqsdNfzhu577t31vMHOPdVzt1CSYRiGEYcpBsMwDCMOUwzwYL4FyBPd9byh+557dz1vsHNPiW4/x2AYhmHEYx6DYRiGEYcpBsMwDCOObqsYuvP60iKyQUSWi8hSEVmcb3myiYj8RkTeF5EVMduOFZE/i8g69//++ZQxG3ic950i0uTe96Ui8tl8ypgtROQkEVkgIqtEZKWI3OJuL+r77nPeKd/3bjnH4K4vvZaY9aWBa7rL+tIisgGoU9WiL/gRkU8A+4D/UtUx7rb/A3ygqrNdo6C/qt6eTzkzjcd53wnsU9V/z6ds2UZEBgGDVPVNEekDLAGmAl+giO+7z3l/jhTve3f1GDrWl1bVw0B0fWmjyFDV/wU+6LT5MuAR9+9HcL48RYXHeXcLVHWrqr7p/v0RsBpn6eCivu8+550y3VUxJFpfujstiKzACyKyxF1Du7txgqpudf/eBpyQT2FyzNdE5C031FRUoZREiMhwoBZ4nW503zudN6R437urYujufFxVxwOfAb7qhh26JerEUrtLPPWXwEigBtgK/Ci/4mQXETkGeAq4VVU/jH2vmO97gvNO+b53V8XQrdeXVtUm9//3gd/hhNa6E9vdeGw0Lvt+nuXJCaq6XVXbVLUd+BVFfN9FJILz4/ioqj7tbi76+57ovNO5791VMXTb9aVFpLc7MYWI9AY+Bazw/1TR8Qxwg/v3DcDv8yhLzoj+KLpcTpHedxER4CFgtar+OOator7vXuedzn3vlllJAG7K1k84sr70v+VZpJwgIifjeAngLO36P8V87iLyGHA2Tuvh7cAdwFzgCWAoTtv2z6lqUU3Uepz32TjhBAU2AF+KibkXDSLyceAVYDnQ7m7+Dk68vWjvu895X0OK973bKgbDMAwjMd01lGQYhmF4YIrBMAzDiMMUg2EYhhGHKQbDMAwjDlMMhmEYRhymGAxEZEBM58VtnTox/i3Hsjzmlu7P6LQ9tkPkOhF5WkRGZ/C4P4lWgIvISyJSF/Pe8NgupdlGRL4gIjvcc10lIl9MY4xKEfmKz/sqIv8d87qHe8w/uK8v9eo6LCL7Ahz/18nuj4g8LCLTEmwfLiJ/H/N6rIg8nOyYRuYwxWCgqrtUtUZVa4AHgPuir1X1Y7mSQ0ROBCao6hmqel+CXaJynQrMAeaLyMAMHHcAUO82nssqbmffIMxx78fZwA9FJNW+PpWAp2IA9gNjRKTCfX0BMdX/qvqMqs5O8ZgdqOo/dqFb8XCgQzGo6nJgiIgMTVceIzVMMRi+RK1DETlbRF4Wkd+LyDsiMltErhWRN8RZ22Gku99AEXlKRBa5/6YkGLNcRP7T/VyDiJzjvvUCUOVaymf5yaWqc9z9/94d87vu8VaIyIPiMFJE3ow57qmxr2O4EvhTwOuRUHbXyv95zH5/EJGzo9dQRH4kIsuAye61W+V6Rr6tkN22JeuBYSIyUURec4/7NxGpdsc/3b0PS90xTwVmAyPdbfd6DP9H4CL372uAx2Lk7zgfcToEvOae810x+5zteldPikijiDwqIuK+1+F1icjNIrLWlfFXsdcJ+IR7Lu/EeA+zgbNc2aOe47M4HQqMHGCKwUiFccCXgb8DrgdOU9WJwK+Br7v73I9j2U/A+cH9dYJxvorTx2wszg/SIyJSDlwKrHe9glcCyPMmMMr9++eqOsFde6ACuFhV1wN7RaTG3edG4D8TjDMFp3d9LI+6P0xLcX5Ak8nuR2/gdVUdh9MK+XLgdFU9A7jL74PiVKqfDLwNNAJnqWot8F3gh+5uXwbudz2MOpxuwbM4ci1negz/OPB5V/4zONKJszP3A790z7lzxWwtcCsw2pUzzhAQkcHAvwL17nujOn1+EPBx4GIczL7HhgAAA0lJREFUhYAr+yuu7FHPcTHgaywYmcMUg5EKi9ye74dwrNgX3O3Lcdx/gPOBn7s/qM8AfcXp9hjLx4H/BlDVRpz2BKelIY/E/H2OiLwuIsuBc4HT3e2/Bm50QzhXA/+TYJxBwI5O266NCa/FrniVjuxtOI3NAPYCB4GHROQK4IDHZ652r+FjOC0MPgD6Ab8VZ77jvphzfA34jojcDgxT1eYk8uDK/xbOfbuGeOXXmSkc8Sb+b6f33lDVzW6DtqUceQ6iTAReVtUPVLUF+G2n9+eqarsbdvILl70PDPZ538ggphiMVDgU83d7zOt2nL5L4DxT9TFzFFWqmnSyMk1qgdWuxfsfwDTXqv0VELXin8JpL34xsERVdyUYpzlm/3RpJf77FDveQVVtA1DVVpwfyyddmbxCWHPc6zdJVaO9rX4ALHC9okuix1DV/8HxtpqBP4rIuSnI/Qzw78SEkTzw6p0T+0y0ceQ5CErs58VzL+dcAyk8o+uYYjAyzQscCSsRE8aJ5RXgWvf903Camq1J5SAiciVOZ9jHOPIjvNP1TjoyXVT1IDAPpyd9ojASOOGdUwIe2kv2DUCNiJSIyEl4tDZ25eunqn8EZuCE54LSjyMTxF+IGfNk4B1V/SlOx9AzgI+APgHG/A3wPXeC14tXORLfvzYFecHpZPxJEekvIj1wwovJSCT7aRRpN9gwYorByDTfAOrcSdBVOPHvzvwHUOKGfeYAX3DDU8mY4cb91wHXAeeq6g5V3YPjJazAUQKLOn3uURyv5gUS8xxO9k8QvGR/FXgXWAX8FGf+IxF9gD+IyFvAX4HbAh4X4P8Ad4tIA/GW+eeAFW7oaQzOOs+7gFfdyXivyWfcMNBPkxz3FpwFnZaT4kqH7tofPwTewLlGG3DCaX68BbSJyLKYyedzcO6TkQOsu6pR9IjIP+NY6f/qs89fcSas9+ROsu6BiByjqvtcj+F3OG3uf5fsczGfLwNexll5sDVbchpHMMVgFDUi8jucZQ3PVdWdPvtNAprdCVkjg7gpuefjhPxeAG7RFH543PTbKlV9KTsSGp0xxWAYhmHEYXMMhmEYRhymGAzDMIw4TDEYhmEYcZhiMAzDMOIwxWAYhmHE8f8ArPyaPwgdzZMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "fitted_cab_model0 = LinearRegression().fit(X_train, y_train)\n",
    "plot_cabs(fitted_cab_model0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.240661535615741"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fitted_cab_model0.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "**Questions**:\n",
    "1. The above code uses `sklearn`. As more practice, and to help you stay versed in both libraries, perform the same task (fit a linear regression line) using `statsmodels` and report the $r^2$ score. Is it the same value as what sklearn reports, and is this the expected behavior?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#### EXERCISE: write code here (feel free to work with a partner)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that there's still a lot of variation in cab pickups that's not being captured by a linear fit. Further, the linear fit is predicting massively more pickups at 11:59pm than at 12:00am. This is a bad property, and it's the conseqeuence of having a straight line with a non-zero slope. However, we can add columns to our data for $TimeMin^2$ and $TimeMin^3$ and so on, allowing a curvy polynomial line to hopefully fit the data better.\n",
    "\n",
    "We'll be using ``sklearn``'s `PolynomialFeatures()` function to take some of the tedium out of building the expanded input data. In fact, if all we want is a formula like $y \\approx \\beta_0 + \\beta_1 x + \\beta_2 x^2 + ...$, it will directly return a new copy of the data in this format!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>1000.000000</td>\n",
       "      <td>1000.000000</td>\n",
       "      <td>1000.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>11.717217</td>\n",
       "      <td>182.833724</td>\n",
       "      <td>3234.000239</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>6.751751</td>\n",
       "      <td>167.225711</td>\n",
       "      <td>3801.801966</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>0.066667</td>\n",
       "      <td>0.004444</td>\n",
       "      <td>0.000296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>6.100000</td>\n",
       "      <td>37.210833</td>\n",
       "      <td>226.996222</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>11.375000</td>\n",
       "      <td>129.390694</td>\n",
       "      <td>1471.820729</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>17.437500</td>\n",
       "      <td>304.066458</td>\n",
       "      <td>5302.160684</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>23.966667</td>\n",
       "      <td>574.401111</td>\n",
       "      <td>13766.479963</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                 0            1             2\n",
       "count  1000.000000  1000.000000   1000.000000\n",
       "mean     11.717217   182.833724   3234.000239\n",
       "std       6.751751   167.225711   3801.801966\n",
       "min       0.066667     0.004444      0.000296\n",
       "25%       6.100000    37.210833    226.996222\n",
       "50%      11.375000   129.390694   1471.820729\n",
       "75%      17.437500   304.066458   5302.160684\n",
       "max      23.966667   574.401111  13766.479963"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transformer_3 = PolynomialFeatures(3, include_bias=False)\n",
    "expanded_train = transformer_3.fit_transform(X_train) # TRANSFORMS it to polynomial features\n",
    "pd.DataFrame(expanded_train).describe() # notice that the columns now contain x, x^2, x^3 values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few notes on `PolynomialFeatures`:\n",
    "\n",
    "- The interface is a bit strange. `PolynomialFeatures` is a _'transformer'_ in sklearn. We'll be using several transformers that learn a transformation on the training data, and then we will apply those transformations on future data. With PolynomialFeatures, the `.fit()` is pretty trivial, and we often fit and transform in one command, as seen above with ``.fit_transform()`.\n",
    "- You rarely want to `include_bias` (a column of all 1's), since _**sklearn**_ will add it automatically. Remember, when using _**statsmodels,**_ you can just `.add_constant()` right before you fit the data.\n",
    "- If you want polynomial features for a several different variables (i.e., multinomial regression), you should call `.fit_transform()` separately on each column and append all the results to a copy of the data (unless you also want interaction terms between the newly-created features). See `np.concatenate()` for joining arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "fitting expanded_train: [[6.73333333e+00 4.53377778e+01 3.05274370e+02]\n",
      " [2.18333333e+00 4.76694444e+00 1.04078287e+01]\n",
      " [1.41666667e+00 2.00694444e+00 2.84317130e+00]\n",
      " ...\n",
      " [1.96666667e+01 3.86777778e+02 7.60662963e+03]\n",
      " [1.17333333e+01 1.37671111e+02 1.61534104e+03]\n",
      " [1.42000000e+01 2.01640000e+02 2.86328800e+03]]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deXhTVfr4P29LoAWBsolQRFARFJACHRWRURRXEFARRR0RHXH5OSo6KM7oiMsIiCPKDKg4KojK4lZRREBBvoorWHYBgQGk7AiyFehyfn8kKWnJTW6Sm7Xv53n6NLm595z3Lsm7nPe8R4wxKIqiKIqXtHgLoCiKoiQWqhgURVGUcqhiUBRFUcqhikFRFEUphyoGRVEUpRxV4i1AJNSvX980a9Ys3mIoiqIkFQsXLtxpjGlg9XlSK4ZmzZqxYMGCeIuhKIqSVIjIhkCfayhJURRFKYcqBkVRFKUcqhgURVGUciT1GIOiKKlDUVERmzZt4tChQ/EWJWXIyMigSZMmuFyukI5TxaAoSkKwadMmatasSbNmzRCReIuT9Bhj2LVrF5s2baJ58+YhHauhJEVREoJDhw5Rr149VQoOISLUq1cvLA9MPYYkJS+/gJEzV7F5TyGNszIZfGlLerfPjrdYihIRqhScJdzrqYohCcnLL+CRD5ZSWFQCQMGeQh75YCmAKgdFUSImaqEkEXldRLaLyDKfbXVFZLaI/OL5X8ezXURktIisEZElItIhWnKlAiNnripTCl4Ki0oYOXNVnCRSlNQgPT2dnJwc2rRpw7XXXsvBgwfDbuvLL7+kR48eAEybNo3hw4db7rtnzx7Gjh1b9n7z5s306dMn7L4jJZpjDOOByypsGwJ8YYxpAXzheQ9wOdDC8zcQeCmKciU9m/cUhrRdURR7ZGZmsmjRIpYtW0bVqlV5+eWXy31ujKG0tDTkdnv27MmQIUMsP6+oGBo3bsx7770Xcj9OETXFYIz5P+C3Cpt7ARM8rycAvX22v2ncfAdkiUijaMmW7DTOygxpu6IoodOlSxfWrFnD+vXradmyJTfffDNt2rTh119/ZdasWXTq1IkOHTpw7bXXsn//fgA+++wzWrVqRYcOHfjggw/K2ho/fjz33HMPANu2beOqq66iXbt2tGvXjm+++YYhQ4awdu1acnJyGDx4MOvXr6dNmzaAe1B+wIABtG3blvbt2zN37tyyNq+++mouu+wyWrRowUMPPeTYucd6jKGhMWaL5/VWoKHndTbwq89+mzzbtlABERmI26ugadOm0ZM0gRl8actyYwwAma50Bl/aMo5SKYpz3H///SxatMjRNnNycnjhhRds7VtcXMyMGTO47DJ30OOXX35hwoQJnHPOOezcuZOnn36azz//nBo1ajBixAief/55HnroIW6//XbmzJnDqaeeynXXXee37XvvvZfzzz+fDz/8kJKSEvbv38/w4cNZtmxZ2TmvX7++bP8xY8YgIixdupSVK1dyySWXsHr1agAWLVpEfn4+1apVo2XLlvzlL3/hxBNPjOAquYlbuqpxLzYd8oLTxphxxphcY0xugwaWxQFTmt7tsxl2dVuyszIRIDsrk2FXt9WBZ0WJkMLCQnJycsjNzaVp06bcdtttAJx00kmcc845AHz33XesWLGCzp07k5OTw4QJE9iwYQMrV66kefPmtGjRAhHhpptu8tvHnDlzuOuuuwD3mEbt2rUDyvT111+XtdWqVStOOumkMsVw0UUXUbt2bTIyMjjjjDPYsCFgbTzbxNpj2CYijYwxWzyhou2e7QWAr5pr4tmmWNC7fbYqAiVlsWvZO413jKEiNWrUKHttjOHiiy9m0qRJ5fZx2sOxQ7Vq1cpep6enU1xc7Ei7sfYYpgH9Pa/7Ax/5bL/Zk510DvC7T8hJURQlYTjnnHOYP38+a9asAeDAgQOsXr2aVq1asX79etauXQtwjOLwctFFF/HSS+78mpKSEn7//Xdq1qzJvn37/O7fpUsX3n77bQBWr17Nxo0badkyumHjaKarTgK+BVqKyCYRuQ0YDlwsIr8A3TzvAT4F1gFrgFeBu6Mll6IoSiQ0aNCA8ePH069fP84880w6derEypUrycjIYNy4cXTv3p0OHTpw/PHH+z3+xRdfZO7cubRt25aOHTuyYsUK6tWrR+fOnWnTpg2DBw8ut//dd99NaWkpbdu25brrrmP8+PHlPIVoIO5Qf3KSm5trdKEeRUkNfv75Z04//fR4i5Fy+LuuIrLQGJNrdYzWSlIURVHKoYpBURRFKYcqBkVRFKUcqhgURVGUcqhiUBRFUcqhikFRFEUphyoGRVEUYNeuXeTk5JCTk8MJJ5xAdnZ22fsjR47Ybuf1119n69atfj+76aabaN68Oe3ateO0006jf//+bN68OWibzz//fEzXwlbFoCiKAtSrV49FixaxaNEi7rzzTgYNGlT2vmrVqrbbCaQYAEaNGsXixYtZuXIlbdu25cILL6SoqChgm6oYFEVRbJCXX0Dn4XNoPmQ6nYfPIS8/euXVJkyYwFlnnUVOTk7ZTOTi4mL+9Kc/0bZtW9q0acPo0aOZMmUKixYt4rrrrgvqaaSlpfHXv/6VunXrMmvWLAAGDhxIbm4urVu35sknnwTcimT79u106dKFbt26We7nJLq0p6IoSUcsl7ddtmwZH374Id988w1VqlRh4MCBTJ48mVNOOYWdO3eydKm73z179pCVlcW///1v/vOf/5CTk2Or/Q4dOrBy5Uq6d+/O8OHDqVu3LsXFxXTt2pU+ffowaNAg/vWvf/HVV1+RlZUF4He/M844w7FzVo9BUZSkI5bL237++ef8+OOP5ObmkpOTw7x581i7di2nnnoqq1at4t5772XmzJlBy2db4VuWaNKkSXTo0IEOHTrw888/s2LFCr/H2N0vXNRjUBQl6Yjl8rbGGG699VaeeuqpYz5bsmQJM2bMYMyYMbz//vuMGzcu5PYXLVpE9+7d+eWXX3jxxRf54YcfyMrK4qabbvI7rmB3v0hQj0FRlKQjlsvbduvWjalTp7Jz507Anb20ceNGduzYgTGGa6+9lieffJKffvoJIGAJbV+MMYwaNYpdu3Zx8cUXs3fvXmrWrEmtWrXYsmULM2fOLNvXt81A+zmFegyKoiQdsVzetm3btjz++ON069aN0tJSXC4XL7/8Munp6dx2220YYxARRowYAcCAAQP485//TGZmJj/88MMxGU2DBg3i8ccfp7CwkE6dOjFnzhxcLhcdOnTgjDPOKFulrXPnzmXHDBw4kG7dunHiiScye/Zsy/2cQstuK4qSEIRadjsvv4CRM1exeU8hjbMyGXxpS13V0A/hlN1Wj0FRlKREl7eNHjrGoCiKopRDFYOiKAlDMoe2E5Fwr6cqBkVREoKMjAx27dqlysEhjDHs2rWLjIyMkI/VMQZFURKCJk2asGnTJnbs2BFvUVKGjIwMmjRpEvJxqhgURUkIXC4XzZs3j7cYCqoYwkLT5BRFSWVUMYRILIt3KYqixAMdfA6RWBbvUhRFiQeqGEIklsW7FEVR4oEqhhCJZfEuRVGUeKCKIUQGX9qSTFd6uW3RKt6lKIoSD3TwOUS8A8yalaQoSqqiiiEMtHiXoiipjIaSFEVRlHLExWMQkUHAnwEDLAUGAI2AyUA9YCHwJ2PMkXjIpyiKEmsSaeJszD0GEckG7gVyjTFtgHTgemAEMMoYcyqwG7gt1rIpiqLEA+/E2YI9hRiOTpzNyy+IizzxCiVVATJFpApQHdgCXAi85/l8AtA7TrIpiqLElESbOBtUMYhIpoiI5/UpInKF5wc9LIwxBcBzwEbcCuF33KGjPcaYYs9umwC/PpSIDBSRBSKyQKswKoqSCiTaxFk7HsNXuK37RsAc4Hbg9XA7FJE6QC+gOdAYqAFcZvd4Y8w4Y0yuMSa3QYMG4YqhKIqSMCTaxFk7iiHNGHMQuAZ4yRhzFXBmBH12A/5njNlhjCkCPgA6A1k+nkgTID7BNUVRlBiTaBNnbSkGEfkDcCPwiWdbeoD9g7EROEdEqntCVBcBK4C5QB/PPv2BjyLoQ1EUJWno3T6bYVe3JTsrEwGyszIZdnXbuGUl2RkreAB4AvjEGLNMRE7GHV4KC2PM9yLyHvATUAzkA+OA6cBkEXnas+21cPtQFEVJNhJp4qzYXV9VRKoDxhiTMGVEc3NzzYIFC+IthqIoSlIhIguNMblWn9vJSuogIvnAamCNiCwUkfZOCqkoiqIkDnZCSW8A9xtj5gKIyAXAeKBd9MRSokUiza5UFCUxsaMYSr1KAcAY86WIlEZRJiVK6LKkiqLYwY5i+FJExgCTcNc2ug6YIyJnAhhjlkRRPsVBAs2uVMWgKM4QDa881p6+HcXgHaCoOHfhLNyK4o+OSqREjUSbXakoqUY0vPJ4ePpBFYMxpktUek5wUjEW3zgrkwI/SkCXJVXskAqWcLSJhlceD08/qGIQkb/5226MecZ5cRKDVI3FD760ZbnzAl2WVLFHqljC0SYaXnk8PH07M59LfP5cuKuetoiaRAlAolU6dIpEm12pJA/R+E5E83uWl19A5+FzaD5kOp2Hz4lZ+epo1DyKRx0lO6GkEb7vRWQE8FnUJEoAUjkWn0izK5XkIZks4Xh6ItHwyuPh6YezHkM13EXuUpZEq3SoKPEmmSzheHr80fDK4+Hp2xljyMedfQTu4nmNgGFRkygB0Fi8opQnmSzheHv80fDKY+3p20lX7ePzuhjYaow5HCV5EgLvDUilbAlFiYRofCei9T3T7LvICVpET0RuMcaMr7DtaWPMo9EUzA5aRE9RlIpUHGMAtyeiiRZHCVZEz47H0E9ECo0xUzwNjgZqOyWgoiiKk6jHHzl2FMNVwMee+kiXA4XGmP7RFUtRFCV8NPsuMiwVg4jU8nnbH/gYmA/8TURqGWP2Rls4RVGUVCMZZnsH8hiW485GEp//vTx/BmgadekURVFSiGSZ7W2pGIwxJ8ZSEEVRlFQnWSoc21nB7U4RyfJ5X0dEBkZXLEVRlNQj3nMs7GJn5vOdxpg93jfGmN3AXdETKXWJV/0WRVESg2SpqmBHMaT7vhGRNNzF9JQQ8MYWC/YUYjgaW1TloCiVh8GXtiTTVe4nNSGrKthJV50tIpOAlz3v7wQ+j55IqUmyxBYVJdWIVxZQoH6TOSvJy2DgbmCQ5/1s4JWoSZSiJEtsUVFSiXhlAQXrN9EUQUXslN0uAf7t+VPCROu3KJUdJyz3UNsIxVN30rOw6vfBqYuBxEpN9YflGIMnfISI5IvITxX/YidiapAssUVFiQZOjLGF04ZdT93pMUCrfkuMSYqxxUAew2DP/z4B9lFskiyxRUWJBkOnLY94jC2ccTorTz2ruovOw+eUfRcPHil2dAzQqt9I240VgSa4bRKRHsCpwFJjzBexEys1SYbYoqI4TV5+AXsKi/x+FsoYWzjjdP7WfHClC/sPFbP7oFsmqx/wUOUL1q8T7caKQLWS/g20B74F/iQi7xtjnomZZElKMtRBUZRYEmjltFDG2MIZp/PnqR84XGypqCKRz1+/D05dTImfpQ0SfWwxUCipK5BjjCkWkRrAPEAVQwCSpQ6KosSSYBa9XcJd8a2ip958yHRb/UU6BujtMxlXgwykGI4YY4oBjDEHPBPblADoXAVFORYrS79OdVdI3wunxuksxx0yXdSoViVib79i1OCajtnMXbkjqaIIgRRDK5/sIwFaet4LYIwxHcLt1FN76b9AG9yVWm8FVgFTgGbAeqCvp/xG0qBzFZRUJNLwqJWl//iVrUOWxYlxOit5hvZsHXHb/qIG7y8ssL16XKKEogMphrZR7PdF4DNjTB8RqQpUB/4GfGGMGS4iQ4AhwMNRlMFxdK6Ckmo4ER5NtIy8aMoTSdQgkULRgbKS1kajQxGpDfwRuMXTzxHgiIj0Ai7w7DYB+JIkUwzhxkAVJVGJ9IculhPaQtk3WhmCkUQNEikUbackhtM0B3YAb4hIO2AhcB/Q0BizxbPPVqChv4M9Jb8HAjRtmlhrBSWaZaQokRLuD50T1m8obSSKtR1J1CCRQtHxGFCuAnQAXjLGtAcO4A4blWGMMbjHHo7BGDPOGJNrjMlt0KBB1IUNld7ts5k/5EL+N7w784dcqEpBSWrCLRMdyPq1SyhtONGfE0RS4SCRSnLHQzFsAjYZY773vH8Pt6LYJiKNADz/t8dBNkVRfAj3h84J6zeUNhLF2u7dPpthV7clOysTAbKzMm0PPCdS2ZxAE9wmGWP6iUg+5a33iLKSjDFbReRXEWlpjFkFXASs8Pz1B4Z7/n8UTvuKojhHuOHRUEIqFccGurZqwNyVO/yHDPy0kZdfQJpIwkwkC3f8IpFC0WL8XEwAEWniKYtxir/PIxmcFpEc3OmqVYF1wADc3stUoCmwAXe66m+B2snNzTULFiwIVwxFUaJExZg/uK3fitazv/0CUbGNQMf7609xIyILjTG5Vp8HrJXk3cdj2fs22gUIWzEYYxYB/oS6KNw2FUVJHOxav/7GBqxIF+GajuWtcavj00VUKUSAnaykD0TkNWPM8yJSDXeo51zg7OiKpihKMmMnpBLKGECJMby/sIDck+qWtWt1fKkxqhQiwM7g89nAaSLyNfAj8BtuxVApycsvoPPwOTQfMp3Ow+ckfF11RUkErL43oY4BVMw0SqRMnlTCjmI4BOwGagM1gJ89q7pVOpxezENRKgOBvjf+MnGC4eslJFImTyphJ5T0IzAD6Ag0AMaJSB9jzPVRlSwBcXpmYqLURVEUfzj1fAb63swfcmHZPhWzkqzWSfD1BhIpkycQyfZdt6MY7vSZc1AAdBeRAVGUKWFxMlc6UWZqKoo/nHw+g31vrMYirDKbKnoDib4AVjJ+14OGkrxKQUTqikhjEWkMzIy6ZAmIk/HMRJmpqSj+cPL5tPp+pIkEHKvzTharU91Vtq1aleSr/p+M3/WgV1lErhCR1bhnLH/v+T8n2oIlIk7GMxNlpqai+MPJ57NrK/+la0qMsTVWd6iotOz1nsKipBvXS8bvuh31+wzQGVhljDkRuBT4KqpSJSiRTHeviGZTxB7NKLOPk8/n3JU7gu6T6DWQIiEZv+t2FEOxMWYHkCYiYoyZDZwVZbkSFqeK5Gk2RWzRjLLQiIV3bGe/ZLS2K5KM33U7g8+/i8hxwNfAmyKyHUieu5KgJEs2RaqQSLXukwGr5xOg8/A5jtRN8ref3WODWduJlAWUjN91y1pJZTuI1AQO4vYubsY9n2Gix4uIK1orSbFL8yHT/RZlE+B/w7vHWpykxG79IzvHVcSqnXD6DFfOykTYtZJEZJYx5hJjzD7PphLgNacFVJRok2jVN5OVcL0ufxazd65CMAva99iC3QepVaWU0sMHuOeVzxiaATef3YRurRuRnp5OWloaxx13HMOnLebgERCRkORUjhIolJR4q+AoSoh4rUd/SiHR47yJRiTxfjtzDYwxbNy4kdWrV7NmzRrWrl3LunXr2LJlC1u3bmXL1q1sOHSobP8C4EGrxiSNtIzjSK9Znyq1j6dKrQbsyzqBOWfBmWeeSf369YPKXJkJpBhqi8jVVh8aYz6IgjwpTyLFPisDlbn6ptPPWiTLVlbEGMPq1atZsGAB+fn5ZX+7d+8u2ycjI4OTTz6Zxo0bc9555zFnw2EOplVHXBmkuTIQVwaSnk7d6i6e7HkGJSUlHDhwgH9M/Z4D+/dSWriXkr07Kd69mUMbFmOOFHLRF68C0KhRIzp06MAf//hHunTpQseOHalatWrY1ybVCKgYgB64w7AVMYAqhhBJxhmQyU5lrb4ZjWdt8KUtbc1E9ocxhhUrVjBv3ryyv23btgFQrVo12rZtS58+fWjfvj1nnHEGp556Ko0aNSIt7WjiZPMh03H5afsw0Ldv97LzrvG/xlQtKe8hVhF49KJsmrCTJUuWsHjxYr7//numT5/uPo/MTC644AJ69erFlVdeSePGjUO8OqlFIMWwwRhza8wkqQQMnbZcM2MiJFQr2EkrN5mIRhZWxbGCrOoujIFBUxYxcuaqY+7F/v37mT17Nh9//DGffvppmSJo0qQJF198Meeffz5nn302rVq1wuU6+pOfl19A34mr2LxnUbl7bOdejpy5iqKSY8OGNTNd3NIth7z8At4/ksbmRi1pfEt/3sitS83f1/Lll18yffp0ZsyYwZ133klubi79+vXjxhtvpGHDhmFdr2Qm0Apu+caY9jGWJySSKSspL7+A+6cs8vuZZsbYQzNU7BPtLCyr6/rwBY0oXP0NH330EXPnzuXIkSPUrl2byy+/nEsuuYTzzz+f5s2blxsYttPusKvbAgS9l4HOe9R1OQGP93o106ZN44MPPmDBggWkp6dz+eWXM2DAAHr27EmVKnYy/BOfYFlJgRRDG2PMsqhJ5gDJpBg6D59jmcudnZVZVmVSscbqGga7fpVxXCfcaxVO+6WH9nNw9bcc+Pn/OLRhMZhSTjvtNK688kp69OjBzurNGDVnna3rH0zuYPcy0PFASNdkxYoVPDpyDNM/mMKRvbuomnU8/W4ZyAuPP0hWVhZ5+QUMnbacPYVFANSp7uLxK1snxbMVydKeCa0Uko1AmRuaGWOPcLNiEr36ZjSIZDzADgW/7efguoXsXzqbwrU/QkkxVbJOoPY51/L1uEdp06YN4FbKj4Uw1hFuJVYvgc57kIXHbtXn6sO1+bnJlZxwx6UUrl3A3gV5THjhaaa88jzdevdjWf0LMNXrlO2/+2ARg99bbHluyURq+EUJjNfCsZpGWKe6K+kfolhRWccLwiXDlVb2A5mV6WJoz8it2XXr1vH666+z5ZWXObJ3F2nVs6jZvjs1Tv8jVRudRpM61cuUAliPdTzx8XK/skR6jwPNMh45c1VIbXtll7R0qrc4m+otzubItrUUL/6ET6a8gaS9Rc32V1DrnD6kV68NQFGJ4cGpya8cAk1w+8IYc5GIjDDGPBxLoVKFYDM+M13pPH5l6xhLlbxE2wpOFfw9d4eLSwMcEZiSkhKmT5/O6NGj+eKLL0hLS6P9uV3Z2uhc0k/qiKS7f0b83Qsra3z3wSLy8guO+fF04h5beRWhtu1P9qoNT6HaJfdR9Q/Xsuebyexd8BH7Fs2g1h+uotY515DmyqDEmKTPNgxURK+RiJwL9BSR9iLSwfcvVgImM1Y59ODOo7+mY+ULcUSCk9VtUxmnKpL+/vvvjBo1itNOO41evXqxevVqnnzySTZs2MCCrz7nP4/cQZN6NQPei0CWvj95onmPe7fP5pqO2aR7Br6DfQcDVUU9qfkp1O/+AI1vG0Pmybn8/s0kNr96FwdWfo0xJukqwFYk0OBzH+A24Dyg4givMcbEfbQ00QefrTIkvFSG7Bgl9kSakbR69Wr+/e9/M378ePbv3895553Hvffey1VXXRVyVk4iZeOFmqEWLENq8LuLKSp1X+lDvy7jt89foWj7/6jWtC31Lr6bqvVPTNhsw0gGn98D3hORx4wxT0VFuhQnWFVJncMQGpUxu8gX3/OvnelCBPYcLDrmWoQbp//xxx8ZNmwYH374IVWrVuX666/n3nvvpWPHjpZy+PZttd03cycUeQKdfzj3P9S5HXaqonrPLePENjTq/wL7F89kz/9NZPP4e2narT/FxZcmZYpr0OqqACLSE/ij5+2XxphPoiqVTRLdY7BTVVLnMNijss5H8GJnvMp7LUK5VsYY5s2bxzPPPMPs2bPJysrinnvu4Z577vE7scuq7Ws6ZvP+woKw5x+Ec/6httFsyHTLz9ZH8B2sKFvJgd38/vnL7Fs5n44dO/LGG2/Qtm3bsNuPBmF7DD4NDMO9MM/bnk33ici5xpi/OSRjylKuMqSF5xBJRk1lsqCdmMkbresVi/sQaLwKyl8LO5auMYZPPvmEYcOG8e2339KwYUOeffZZ7rjjDmrVqhWSHIVFJUz6/tdjChV6ZfLOEYilte+PdIsKu+kWk+28BLu/3tde7yG9Rh1O6fcPumX8jwnPPUbHjh159tlnue+++8om9gVr0653GC3s+DjdgRxjTCmAiEwA8gFVDDbwflGtLJ5wM2oqW92lSFfyitb1itV9sHOevvtYZeYYY/j000957LHHyM/Pp1mzZowdO5YBAwaQkZERthz+fnB99490LokTK7lZyWi1HUK7v76ZX7sPFvFJ0Uk8N/kL3n/xUQYNGsScOXN44403+GrjoYBtVuzTNwwXq++5naU9AbJ8XteOhiCpjtPZFqmwFm4oWHlWBmyt32zneoWzJnSs7oMdzzLQPsYYnnh5MrVOak2PHj1YsWErf3liFKtXr+auu+6ypRQAamf6K2NnbXV7ZYp0vW0n1k3OttjX33avvPdPWWTr/lo9B6/8sIO8vDxeeOEFPvvsM3JycnjslfcCtmnXO4wmdhTDMCBfRMZ7vIWFwD+jKlWK4tR60ZAaa+GGgr91c73YWb/ZKpTn3R7umtCxug+Bzh8Ce59fffUVbf9wLkPv6kfhnp3UvewvNLz1JWYXtWL6su22ZcjLL+DAkeJjtrvShH5nn2i5rrET6207sW6y3TZ85bWi4v0N9ByICPfddx/ffvst1apVY9mrD7Iv/1PLNkP1DqNBUMVgjJkEnIO7zPb7QCdjzJSoShUDjDGMHTuW5cuXx7TfSC0nL05YUMmA93oNmrKIDFcaWRYWazArysqi9W4P1/J3+j4Eej4yXEe/rpmuNOpUd1l6n3n5BZz5/8aS2dy95sCqVaup0+0OsgeOo2a7S5H0KiFbnlaVS4/LqMLTvdtaesROeFVOeNx22whmscOx99fOc9CxY0cWLlxIVotcfps1ll2zxmJKio/ZN1Lv0Als5VEZY7YA05zsWETScc+PKDDG9BCR5sBkoB5ur+RPxpgjTvbpy2+//cbQoUPJzs7m+++/j8kiHU7GoyvDLOCK12v3waKAVnMgKypYfDlcy9/J+2D1fCzY8NsxGT8glgXbXpo2n4eHPMK+n78iLbMWdbreynHtryDNdWy4KBTL02rfPQfdMXCrcQSnvConal7ZaSOQpwD+76/d56B27dr8962pDLz3QX779j2Kdm6kQe9HOF9yU5sAACAASURBVK523bJ9/bUVrH+niWeC7X3Az4A3BWIEMMoYM1lEXsY9ue6laHVer149XnvtNXr27Mk//vEPhg8fHlF7djJTnKyRbyfzJNmoeA0PHin2e72ssksaZ2WWtVGwp7Bsv+ysTOpUd7H74LG59Nk+Vlo4uf9O3odwMn58+9m5cydPP/00o/8zBtLSqd25H7X+cBVp1apb9ul7zYLJH+41sjrOOz6UaM+t1fMF7ufFn7yhPAfX5DYlfcwLPDS8BWs++Bc733mIv78xtWzfim3FIyvJ1jwGxzsVaQJMwD1W8QBwJbADOMEYUywinYChxphLA7XjxDyGgQMH8t///pe5c+dy/vnnh9WG3RzraNfIT2bszPnwJdOVbiuX3osrTUAoFwoJN/c/WgTKs/eH97kpLCxk9OjRDBs2jH379lG97cXUPu9GqhxXN+DxweYfVDzvcK9RKHMwEoFozXfwx7fffkuPHj1wuVzMmDGD9u1jswROsHkMAccYRCRdRFY6LxYvAA8B3vyuesAeY4w34LYJiMlT8vzzz3PKKadw88038/vvv4fVht0YqhPxaKfGKBINO3FdL97YcMVY8dyVOyzbKCo1FJcYy7h8ODFsp+9FsHz6ijSqVY0333yTli1bMmTIELp06cLSpUs58/qHgioF75rX/q6ZVfw/3Di/73H+SLRsulCyl+xi9ax06tSJr7/+mmrVqnH++efzxRdfhN2HkwQMJRljSkRklYg0NcZsdKJDEekBbDfGLBSRC8I4fiAwEKBp06YRy3Pcccfx1ltv0blzZ/7yl7/w5ptvhtyG3RhqpPHoVJ67YDfe7L1e/mLFVvX2vRjgUFEpo67LsSyBEO5MXCfuRaB8+ooeUummpWxd8Bb9Vy2nY8eOvPnmm1xwwQUADD5cO6j35V3zOtQ1CsKN83uPs/KaEymbzunxu2DPyumnn84333zDZZddxuWXX87kyZO5+uqrIz+RCLCTrloHWC4iX4jINO9fBH12xl2xdT3uweYLgReBLBHxKqomgF/zyxgzzhiTa4zJbdCgQQRiHOXss8/m0UcfZeLEibzzzjshH2/XE4g0syKV5y5YXcOsTJft62XH83LqekXjXgSyVL3PTfHv29j3yQh+ffsRzOEDvPPOO/zwww9lSgGCW+gQPAMmWlkvyZBNF485R9nZ2Xz11Vfk5ubSt29f3n333UhOIWLsDD4/5mSHxphHgEcAPB7DX40xN4rIu0Af3MqiP/CRk/0G49FHH+Xzzz/njjvuoGPHjrRsGVp+tF0LI5LMilSeu2B1DXu0a8TclTvCbsMf/q6X1aC11UCf1TUv2FMYdomMrq0a8PZ3G8tZ1N7n6JKWdfjpwy8ZOX4kaWlp/POf/+SBBx6wnJhmd8Z9rLPbomGNRyMBw4kMKC+BnpXmQ6aXk3vmzJlcccUV9OvXj5KSEq6//npHZAiVoIrBGDNPRE4CWhhjPheR6oB1zmD4PAxMFpGncZfceC0KfVhSpUoVJk+eTE5ODn379uW7774jMzPyVaOcJJVXMPN3Dbu2alBuYDRYuMZObSo49nodUwTNE9IJ1F+gyrm+5Zjthpjy8gt4f2FBOaUgwNUdGnN49de06jWYX3/9lRtuuIERI0bQpEkTy7Z8CfZsxjq7zcn+kiW0GuhZ8Z3wB265Z8yYQffu3bnxxhspKSnhxhtvjKG0boJmJYnI7bhj+nWNMaeISAvgZWPMRbEQMBDRqK46Y8YMrrjiCgYOHMgrr7ziaNuREk5WSDIX2ot0QXu718uqn0D95eUXMGjKooDrbdiVOS+/gAenLj5mjOHI9nXsnfMqBzYsJScnh9GjR9OlSxebPaY+kT4fscJuxp1A2fjXgQMHuPLKK5k3bx6TJk2ib9++jsoUUVaSh/+He1xgL4Ax5hfgeGfESzwuv/xyhgwZwrhx48Iab4gmocY+nShFEE8iDZ3ZvV7BJjT56693+2zbSsGqDTh6j3yVQknhXnbNeokt4++ncPsG7vzbcBYsWKBKoQLJElqt+BxaYXB7m3n5BdSoUYNPPvmEzp07c+ONN/Lpp8eW0IgmdjyG740xZ4tIvjGmvWeA+CdjzJmxEdGaaK3HUFxczAUXXMCiRYv4/vvvad06Oddlznlilt8FUhLNoqqIb7zfH07Lf8ojnwbMCEoX4V992x2jUIJ5Gr74yuzrxaX5TKYypSXsX/QZe756i9LDB6jZoTu1z7uRqpk1KTXGtseXzF5iKCSLx1CRYM9Nneou8v9xCQB79+7loosuYtmyZcyYMaNckkEkOOExzBORvwGZInIx8C7wsSPSJShVqlRhypQp1KxZk169erF79+54ixQyefkFfpUCJJ5F5UuwAmbRGBgNpBS8n/vztPwVZXOli3synQ++Mlf04rx9H9q4lC3j7+O32S/hOr45jQaMpm63O0jPOI4SY2x7fMnuJYaCE4X14kEw+XYfLCq7X7Vq1eKzzz7jlFNOYeLEibEQD7CnGIbgnpW8FLgD+BR4NJpCJQLZ2dm8//77bNy4keuvv56SEnuTrxKFQGmTiTxYHWiiWzhpg8EmoeXlF9iaWOYvFdVfqGpkn3aMvLadZfiq4vkV793Ojo9GsG3SI5QePkD93o/Q8Pp/UrVBM9ty+JLKKc0VcTqtNFb0bp9Nner+i0F68b1f9erVY968eYwbNy7aopVhJyup1FNu+3vcYbBVJh51NOLAueeey9ixY7n99tt55JFHePbZZ+Mtkm0CeQWJbFFZyS0QcnggWNaKv/h+qLJZpTVa/Th52ygtOszeHz5g73fvAYbanW+g1tlX+y10Z0cOL8HKi6caTqaVxpLHr2zN4PcW+61WC8fe43r16sVCrDLsLO3ZHXgZWIv7+9lcRO4wxsyItnCJQP2Ol9Pw7J6MHDmSvI0unn347qR4EK1S5OpUdyW0/E6m5AYrWjh02nLbZThCkcFfjN8rT6kxFK7+lt/mvkbJ79uo3vI86nS9lSq17edzBJLDqgCciDu2nerjDsmC99oPmroIf3ZJvL16OxPc/gV0NcasARCRU4DpQMorBq9FWa3LAKptXscv743kvoy6cN/1Cf+lsppI9PiViT2Q7uQEqEBZK4HGYLx9hiODPy9l8HuLwcCBbevZ/cU4Dm1YjKv+SdS//hkyTgothyOYHFbejzFHvYZEzfevbHivfSKWz7czxrDPqxQ8rAP2RUmehMJrcUq6iwZX/Y0qtRqwaeqTPPnW5/EWLSjJHH91Su5A5RcCxdytivTZkcGfl3LowD62zXqZLW/8hSPb1lL34jtpNGC0baWQLmJbDruF3lJ13CFSYl2kMlG/p5bpqiLireJ0MXASMBX3GMO1wEZjzN0xkTAA0UpX9VKx4FfR7i1snfggadWqs3n1Epyq1aREh0AT3AJNTnvBosieHXyfGVNawv4ls9nzf29SWriP43IuI6vLTaRXt79suitdGNnn2FRZK0IpX67l3suTCKXXY0WwdNVAoaQrfV5vA7yLFewAEjetxUEqxrtddRpx/DX/YPvkv9GzZ0/mzJlju2xGtKgsOevhEKj8gtU8CW+2SKB4fKBr7n1mDm1awe7PX+HItrVUa9Kaut3uoGrDk0OSv051l+UqbVYy+DvnA4eL/YbN4h3HTjScXEgr2YnLQj1OEW2PwcqC6JW1iWcHD6RHjx68//77uFyBU89iLV8qWjhOY3Xtgi1cE+ya//ezH3nwrw+xd/mXpB9Xjzpdb6X66X9EQlhrwZUmjLzW2ksI9b7rc2KPyrSQViQeg7eB5sBfgGa++xtjejohYLwJZP1ZW5yXcVL1Iu6++2769+/PxIkTSU+PRl3BwKiFEz5W9zbYPAB/NY0Ki0oY/skSln7yOsOHD6e4pIQmF96EtOuNq1qm7XRYL8dlVAl4/6xkfHDqYgZNWRT3QnnJSioXqQwVO1lJebgrnX7M0RXXUgI71Rmt8qTvuusu9u7dy5AhQzjuuON45ZVXQrIKnSBZasUkKv7u7f0WC9d4n42KP/LGGA6ums+Cua/z/d7tXHvttTz77LMs2u0KaalSX/YctM6WAuv7G6gqbLLm+8eSWJcgT2TsKIZDxpjRUZckDkRqcT/88MPs27ePf/7zn9SsWZPnnnsupspBLRznCbQQfMVn5cj2dfz2xasc3riU6ieczPSPppTVsrlx8pywlAJA7czAoclAZZx9ZVXPMTRC8axSfWzPjmJ4UUQeB2YBh70bjTE/RU2qGOGExf3UU0+xd+9enn/+eapUqcLw4cNjphzUwnEeO2GfkoO/s+ert9i/eCZpGcfR8PJ7GPPEX7ngDyeV7ROJ13bgSDF5+QWWPzSRLEikBMaOZ5Us60BEgh3F0Bb4E+4lOL2hJON5n9Q4YXGLCC+88AJFRUU8++yzHD58mFGjRjmqHELJQEk1yyXWZFs8E+kiFBcXsS9/Or9//Q6lRwqp2aEHdbvcwAs3n1c2MO2vamqoFJWYgNZ+xftu1Zd6jtEhHmN7sfZQ7CiGa4GTjTFHoiZFnHDK4k5LS2Ps2LFUq1aNF198kcOHDzNmzBjS0uzMHwxMMOtEY8fO4u+ZyKiSxpmlv/DRG89xZNevZDRrT52Lbqd2o+aW2UrhKgUvwax93/sebPlOxVliPbYXDw/FjmJYBmQB26MiQRxx0uIWEUaNGkVGRgYjRozg0KFDvPrqq1SpYucSW6OZR7Gl4jNRa/9GSr57k3cXfkfjk06m1hVPUdgoh+w61cs9K1ZVYb1jFiL4rYljRSjWvnqOsSXWY3vx+A2w86uVBawUkR8pP8aQEumqTlrcIsKwYcPIzMxk6NCh7Ny5k8mTJ1OjRo2w29TMo9jTu302OXWK+Pvf/84777xDgwYNGDNmDLfffrvlnBWr+1FqDC9clxNShlI41r56jrEj1mN78fgNsKMYHo9a7ymIiPD4449z/PHHc88993DhhRfy8ccfc/zx4a2GqplHsWX37t0888wzjB49mrS0NP7+97/z0EMPUatWrYDHBbpPgdaYqEi6iE48S3DC8dAqjhF0bdWAuSt32Do+Hr8BOvM5inz00Uf069ePxo0b89lnn3HqqaeG3IbOWo0NR44cYezYsTz11FPs3r2b/v3789RTT9GkSRNbx4dbl6kiqTjLtrJjp35VrGeuOzHzeR+UPddVARdwwBgT2IRS6NWrF3PmzKFHjx6cffbZTJ06lYsuuiikNipD/DieOeHFxcVMnDiRJ554gg0bNtCtWzeee+452rVrZ+t4X9lrZ7rIcKWx52CRrbpM/rBjBUZifQZrK9WerUTAjscYaMwgHr8BIXkM4s7B7AWcY4wZEjWpbJLoHoOXNWvW0Lt3b37++Weee+457r///pjPkk5U4uURlZaWMnXqVB5//HFWr15Nbm4uTz/9NJdcconte2NXdn/7udIEhHIreNk570itz3DkVyLDqgZTRWLpLUbsMfjiWdIzzzPhLe6KIZo4aUmdeuqpfPvtt9xyyy088MADLFy4kFdffTXulVljRV5+AUOnLS+r8OlbNTTWGRfGGKZNm8Zjjz3G0qVLadOmDR9++CG9evWypRCCzVXwJ7uVxedvW7BzDtX6DPQca8ZbbLAzU927X6JgJ5R0tc/bNCAXOBQ1iRKAaOQN16xZk3fffZdhw4bx2GOPsXjxYiZPnkzr1om9olqk5OUXMPjdxRSVHv0B3X2wyL2qGbHLuDDGMHv2bB577DF++OEHWrRowTvvvMN1111ne76J3bkKTqwNbYXdkJR3lbpAz7FmvMUGOzPVE23eiR2PwXddhmJgPe5wUkphxxJ84uPlEXkR3iyXjh070r9/f3Jzc3n++ee58847Uza0NHLmqnJKwYt3dm+0My5KS0v5+OOPeeaZZ/jhhx9o2rQpr732GjfffHPIc0zsZheFKnso3mmgWk4VZQjmEWjGW2zw5zFGMi4UC4J+M4wxA2IhSDyxawnuPljEbk/ly0i8iMsuu4wlS5bQv39/7r77bmbNmsW4ceNSckW4QNbn5j2FjPKT4++E9VRSUsLUqVN55plnWLZsGSeffDKvvPIK/fv3p1q1amG1aceSDlX2UL1TO0rBK8Mgi0qx3vPQWluxI9nmmVgqBhH5R4DjjDHmqSjIExdCyTP3xU481soabNiwIZ9++ikvvPACQ4YM4fTTT+fFF1/khhtuCMt7SNTskkDx1cZZmY5nXBw+fJiJEycyYsQI1qxZwxlnnMHEiRO5/vrrI56FbnUu6SKUGhOW7KGurRCollNFGayyobweQWXIeFPCI9A35YCfbTWA24B6QMoohkhiqoGODWYNpqWl8cADD3DppZdy2223cdNNN/H222/z8ssv07RpU9syJHK1x8GXtjxmjAHcaxl7LVMnrKlt27bx8ssvM3bsWLZv307Hjh354IMP6NWrlyM1qwC6tmrAW99tPGZ7v7NP5OnebcNqM9S1FaysfH+ZRHY8gmSzZJXYYPmNMcb8y/sHjMO9zvMAYDIQ2uK1CU6wmGp2ViZZFjXyAx0bbDUwL61bt2b+/PmMGjWKefPmccYZZ/DMM89w6JC9MX67/cSD3u2zGXltu3LXr051V0gL3Adi8eLFDBgwgKZNmzJ06FByc3OZNWsWP/74I1dddZVjSgFg7sodIW23g514vteD8JbiHnZ1W7KzMhHcz6ZVemko+yqKLwHnMYhIXeAB4EZgAvCiMWZ3jGQLilPzGPLyCyxX7vLmFoeT8x3OGrL/+9//eOCBB8jLy+Okk05i5MiR9OnTJ2B4qTKtVQtQVFTEtGnTGDNmDHPnzqV69erccsst3HvvvbRsGb34eDSus515CV50joHiFMHmMViaUyIyEvgR2Ae0NcYMdUIpiMiJIjJXRFaIyHIRuc+zva6IzBaRXzz/60Tal116t88O6hGEY31ZWYOBrMTmzZvz4YcfMmfOHLKysujbty+dO3dm9uzZWCnxcPpJRtasWcOQIUM48cQT6dOnD2vWrGHEiBFs2rSJMWPGRFUpQPjXOS+/gM7D59B8yHQ6D59DXn5B2WcVn6v0AAZAoniBgQh0rkryYOkxiEgp7mqqxVDOUBLcg89hlcQQkUZAI2PMTyJSE1gI9AZuAX4zxgwXkSFAHWPMw4HacnLmczRmgUbaZklJCW+88QZPPvkkv/76K506dWLo0KFcfPHF5TyIVJ7Bum/fPvLy8njjjTeYO3cu6enpdO/endtvv53LLrss4gHlUAjnOod6TDAPIpG9wFR+DlONsD0GY0yaMSbTGFPTGFPL569mJHWSjDFbvMuCGmP2AT8D2bjnRkzw7DYBt7KIGdGKx2a4jl7irExXSG2mp6fz5z//mV9++YWXXnqJTZs2cemll5Kbm8v48ePLxiBSLZZ85MgRpk2bxvXXX0/Dhg25+eabWb9+PU8//TQbN27ko48+okePHjFVChDedQ51/Mfbh5XnkMheYCKPdSmhEdfqqiLSDPg/oA2w0RiT5dkuwG7v+wrHDAQGAjRt2rTjhg0bYiZvKETDejp8+DBvvvkmL774IsuXL6d+/frcfvvt3HrrrWFVbk0kDhw4wOzZs5k2bRp5eXns3r2b+vXr07dvX2644QY6derk6EByrAh3XCIZre/KNtaVzITtMUQbETkOeB+43xiz1/czT00mvxrLGDPOGJNrjMlN5Alhdq2nUGKy1apV4/bbb2fp0qV88cUXnHfeeYwYMYIWLVrQqVMnxowZw44d4WfIREqo8eVNmzbx6quvcuWVV1K/fn2uuuoqPvzwQ6644gqmT5/O5s2bGTNmDJ07d05KpQDhj0skoxdYWca6KgNx8RhExAV8Asw0xjzv2bYKuMAYs8UzDvGlMSbgaGIiV1e1Yz05YRVu2rSJSZMmMXHiRJYuXUpaWhrnnnsu3bt3p3v37rRp0yYm5TbsnMvWrVuZO3du2d+aNWsAaNasGb169aJnz5506dLFcpW0ZCQZLf9wqUznmuwE8xhirhg8YaIJuAea7/fZPhLY5TP4XNcY81CgthJZMXQePsfvrNPsrEzmD7nQ9j6hsGTJEt59912mT59Ofn4+AI0bN6Zz585lf2eeeSZVq1YNue1gVDyXksJ9HNm2lvRd6zivzj4WLlzIunXrAKhduzbnn38+Xbt2pVu3brRu3Tpla0VB4s5KjwaV6VyTmURUDOcBXwFLgVLP5r8B3wNTgabABqCvMea3QG0lsmKwYz1FMyZbUFDAp59+yty5c5k/fz4bN7pn7FapUoXTTjuN1q1b07p1a5o1a0aTJk1o0qQJjRo1ombNmkF/pI0xHDx4kK1bt7J582a2bNnCnS/NpGj3Zop2F1D822ZK9u8q279hdlO6dDqLs846i65du9K+fXvS09MjOj9FUcIn4RSDk0SqGKJt3QRr32mPIRCbNm1i/vz5LFmyhGXLlrF8+XLWrVt3zNwIEaFWrVrUqlWLjIyMMiVhjKGwsJB9+/axf/9+SkqOTadMy6yFq05jqtTNxlWvCVVPaEHVhqfQtNHxfs9HrUtFiQ+qGCxIhHhovGU4dOgQBQUFbNq0iU2bNrFlyxb27t1b9ldYWF5pZWRkULNmTWrVqkXNmjU54YQTaNy4MY0aNWLRLuHRGev99uPPA4r3uStKZUYVgwWxtNYDEU2rOdYWec4Ts8pWafOlTnUX+f+4pNw2q+tfsUoopE71T/WQlETB0aU9UwmnV68K90sfreqW8ai4OrRnawa/t7jcOsYA+w8VlxWA82Knqujg9xaDoawyayJVja1IsPsfzv2w+0ypwlGcJjmTwx3AyZxr75e+YE8hhqNf+njWiYnHLNTe7bOpUfVYW6Oo1BzTr53rXFRijinXnYgzae3c/1Dvh91nKhGfPSX5qbSKYfClLcl0lc+MCXf1qkQsBRCqR+RU8bPf/YSS/PXr7/rbJdHWJLZz/0O9H3afqUR89pTkp9IqBidnlibiouqheEROWp1Z1f1PTqu4vXf7bK7pmE04sxcSbSatnfsfqodq95lKxGdPSX4q7RgDOBfft7OoujcOXLCnsGxB9+wI4sHB4sqhrOcbbNH4UDhkURXUX47D3JU7/Nc98eBKl3JjDF4OHjl2zCKe2Ln/oa6vbKfNUPZTlFCotB6DkwQLS/la5HDsso2hWuZ2LPxQPCKnrM68/AIKi0r9fuYvxBSo/eysTEb2aXfM6m8Auw8WJVQc3U5YMlQP1W6o08mQqKJ4qdQeg1MEW1Tdn0XuJRzL3K6Fb9cjcsrqDBTX9teWVb8VU4ZHzlx1TBpsuB5NNAh2/333sytvKG3a2U9RQkEVg0ME+tIHs7xDtcydjiuHGuYIVS5vH+H2mwxx9GikHdttM1opz0rlRRWDQ+TlF/DEx8vZfdBt2WZluhjas7V72dDqrrLt/kgTofmQ6bZz1dM8YxT+2gkn9u5rdXrHQHwzW+y0l5dfYClXneouy8Xqvf36WrvgngDn3WZ1/TSOrvhD53VEjioGB8jLLzhmYteewiIGv7uYBRt+Y/+h4oDHVxxzgPI/xhUnR/n78fVuD3cCmHf/cCbFeeXzJ1emK53Hr2wdsN9gE8FcaYIrXcpdX42jK/6Ix8TOVEQHnx1g5MxVx8z2BXc2zaTvfz0mq8aLv+UbK+ag5+UX8ODUxZZjFMGOD4VwJmF1Hj6H+6cs8itfmhByCrA/GYpKDTWqVgk6cKsL0Ss6r8MZ1GNwgECxbivrXoBSi8+87QWyxMOVJ5zj/G0Ptmg9gIU+DEuG3wuLWPT4JX4/8yePWoqVk2QYj0oGVDE4gFV2DVA2Z6EiJsBntT3pmYGymYK13WzI9KDzJCrGYmtnuvwWwWuclXnMvgePFNvyYoZOWx5SvDfcDCkn52IoyYvO63AGDSU5wOBLW7onY1XAlSb0O/tEy9IPVp7AAc8ErkBWTqYrPWDbEHiehL+5EAeOFONKK38ema50urZqcMy+gQbTfdlTWBTSjOpw8/LVUlRA53U4hSoGB+jdPpuRfdpRx6fsQ1ami5HXtuPp3m3LJjbZpajEMHTactIsVlJLF2HY1W3L2vY3VuGlsKiE+6csOibm7jeWX6FoXVami2FXt2Xuyh22xziCESzeG26pEl2IXgFnS91UZirtegzxoNmQ6Y60c9M5TXm6d9uy91ZLhFbEdyEcO8e40oSR17Zj0JRFttq3ixNLl1ZEF/5RFPsEW49BPYYYEsiyD4VPFm8p996uVexrrdsqe+0pl+201R0NK14tRUVxDlUMMSTU7CIr9hQWlQsLhVLC2htzt3vM5j2FEZXIrkg0472922czf8iF/G94d+YPuVCVgqKEiSqGGGI1zlCnuqucpVvHonS1Lw9OXVymHHytZQjsmXitdW/Z62BeTJoIg6YsIsOVdkwxu1BRK15RkgNNV40hVrWBHr+ydcDZv/6oOMs52Axib1++FV/fX1gQ1Ivxfr77YBGZrnRuOqcp7y8sCGkwWmP9ipJcqMcQQyrGwbMyXWS40hhUIWuoogdgRaAMn2Ax9yc+Xh7wx92fI1FYVMLclTvKtVunuiugJ+HNoIqFUtCZz4riDJqVFCfsZtEE8x7CyfDJyy/g/imLArYXKINqvZ/+rLKcopGB5A/NSlIU+2hWUoJit6aL1/K3GgsIJ8PHzroJgcYe/Fni8Z5HoDVyFMU5VDHEiVBm6vZun82/+rZzbEannXUTAo09+Ju9HO8ZpzrzWVGcQwef40SoNV0CrdRlp/68nfUcsjKPrpuQHaD+k9VqcVbyxQKtkaMozqGKIU6Es2qav5W67FQVtbOeQ6YrnaE9j66b4E8+X6w8m3jF851ahU5RFFUMccMpC9tOVVGrKq3pIpQa47dv7+sHpy72q0gSzRKPt8eiKE4S71XoVDHEEScsbDuxdat9So0JmDHkb1U3SFxLXNc+VlKBRFhbJKEGn0XkMhFZJSJrRGRIvOVJBuxkA0WSMaQ1iBQltiRChl3CUZw/aAAACf9JREFUeAwikg6MAS4GNgE/isg0Y8yK+EqW2NiJrUcaf1dLXFFiRyJk2CWMYgDOAtYYY9YBiMhkoBegiiEAdmLrGn9XlOQhETLsEkkxZAO/+rzfBJxdcScRGQgMBGjatGlsJEtw7Fj0avUrSnKQCBl2CTXGYAdjzDhjTK4xJrdBgwbxFkdRFMVREmFcL5E8hgLgRJ/3TTzbFEVRKhXx9vATyWP4EWghIs1FpCpwPTAtzjIpiqJUOhLGYzDGFIvIPcBMIB143RizPM5iKYqiVDoSRjEAGGM+BT6NtxyKoiiVmUQKJSmKoigJgCoGRVEUpRyqGBRFUZRyJPXSniKyA9gQYTP1gZ0OiJNsVNbzhsp77pX1vEHPveK5n2SMsZwIltSKwQlEZEGgtU9Tlcp63lB5z72ynjfouYd67hpKUhRFUcqhikFRFEUphyoGGBdvAeJEZT1vqLznXlnPG/TcQ6LSjzEoiqIo5VGPQVEURSmHKgZFURSlHJVWMVTm9aVFZL2ILBWRRSKyIN7yRBMReV1EtovIMp9tdUVktoj84vlfJ54yRgOL8x4qIgWe+75IRK6Ip4zRQkROFJG5IrJCRJaLyH2e7Sl93wOcd8j3vVKOMXjWl16Nz/rSQL/Ksr60iKwHco0xKT/hR0T+COwH3jTGtPFsexb4zRgz3GMU1DHGPBxPOZ3G4ryHAvuNMc/FU7ZoIyKNgEbGmJ9EpCawEOgN3EIK3/cA592XEO97ZfUYytaXNsYcAbzrSysphjHm/4DfKmzuBUzwvJ6A+8uTUlicd6XAGLPFGPOT5/U+4GfcSwen9H0PcN4hU1kVg7/1pSvTgsgGmCUiCz1raFc2GhpjtnhebwUaxlOYGHOPiCzxhJpSKpTiDxFpBrQHvqcS3fcK5w0h3vfKqhgqO+cZYzoAlwP/zxN2qJQYdyy1ssRTXwJOAXKALcC/4itOdBGR44D3gfuNMXt9P0vl++7nvEO+75VVMVTq9aWNMQWe/9uBD3GH1ioT2zzxWG9cdnuc5YkJxphtxpgSY0wp8CopfN9FxIX7x/FtY8wHns0pf9/9nXc4972yKoZKu760iNTwDEwhIjWAS4BlgY9KOaYB/T2v+wMfxVGWmOH9UfRwFSl630VEgNeAn40xz/t8lNL33eq8w7nvlTIrCcCTsvUCR9eX/mecRYoJInIybi8B3Eu7vpPK5y4ik4ALcJce3gY8DuQBU4GmuMu29zXGpNRArcV5X4A7nGCA9cAdPjH3lEFEzgO+ApYCpZ7Nf8Mdb0/Z+x7gvPsR4n2vtIpBURRF8U9lDSUpiqIoFqhiUBRFUcqhikFRFEUphyoGRVEUpRyqGBRFUZRyqGJQEJF6PpUXt1aoxPhNjGWZ5Jm6P6jCdt8Kkb+IyAcicoaD/b7gnQEuIl+KSK7PZ818q5RGGxG5RUR2eM51hYjcHkYbWSJyd4DPjYi85fO+iqfPTzzve1pVHRaR/Tb6/2+w+yMi40Wkj5/tzUTkBp/3bUVkfLA+FedQxaBgjNlljMkxxuQALwOjvO+NMefGSg4ROQH4gzHmTGPMKD+7eOVqAUwB5ohIAwf6rQec4yk8F1U8lX3tMMVzPy4AnhGRUOv6ZAGWigE4ALQRkUzP+4vxmf1vjJlmjBkeYp9lGGP+HEG14mZAmWIwxiwFmohI03DlUUJDFYMSEK91KCIXiMg8EflIRNaJyHARuVFEfhD32g6nePZrICLvi8iPnr/OftrMEJE3PMfli0hXz0ezgGyPpdwlkFzGmCme/W/wtPkPT3/LRGScuDlFRH7y6beF73sfrgE+s3k9/MrusfL/47PfJyJygfcaisi/RGQx0Mlz7VZ4PKOApZA9ZUvWAieJyFki8q2n329EpKWn/dae+7DI02YLYDhwimfbSIvmPwW6e173Ayb5yF92PuKuEPCt55yf9tnnAo939Z6IrBSRt0VEPJ+VeV0icpuIrPbI+KrvdQL+6DmXdT7ew3Cgi0d2r+f4Me4KBUoMUMWghEI74E7gdOBPwGnGmLOA/wJ/8ezzIm7L/g+4f3D/66ed/4e7jllb3D9IE0QkA+gJrPV4BV/ZkOcnoJXn9X+MMX/wrD2QCfQwxqwFfheRHM8+A4A3/LTTGXftel/e9vwwLcL9AxpM9kDUAL43xrTDXQr5KqC1MeZM4OlAB4p7pvrJwBpgJdDFGNMe+AfwjGe3O4EXPR5GLu5qwUM4ei0HWzQ/GbjeI/+ZHK3EWZEXgZc851xxxmx74H7gDI+c5QwBEWkMPAac4/msVYXjGwHnAT1wKwQ8sn/lkd3rOS4AAhoLinOoYlBC4UdPzffDuK3YWZ7tS3G7/wDdgP94flCnAbXEXe3Rl/OAtwCMMStxlyc4LQx5xOd1VxH5XkSWAhcCrT3b/wsM8IRwrgPe8dNOI2BHhW03+oTXfFe8Ckf2EtyFzQB+Bw4Br4nI1cBBi2Ou81zDSbhLGPwG1AbeFfd4xyifc/wW+JuIPAycZIwpDCIPHvmX4L5v/Siv/CrSmaPexMQKn/1gjNnkKdC2iKPPgZezgHnGmN+MMUXAuxU+zzPGlHrCToHCZduBxgE+VxxEFYMSCod9Xpf6vC/FXXcJ3M/UOT5jFNnGmKCDlWHSHvjZY/GOBfp4rNpXAa8V/z7u8uI9gIXGmF1+2in02T9ciin/ffJt75AxpgTAGFOM+8fyPY9MViGsKZ7rd7Yxxlvb6ilgrscrutLbhzHmHdzeViHwqYhcGILc04Dn8AkjWWBVO8f3mSjh6HNgF9/jxXIv97naUnhK5KhiUJxmFkfDSviEcXz5CrjR8/lpuIuarQqlExG5Bndl2Ekc/RHe6fFOyjJdjDGHgJm4a9L7CyOBO7xzqs2urWRfD+SISJqInIhFaWOPfLWNMZ8Cg3CH5+xSm6MDxLf4tHkysM4YMxp3xdAzgX1ATRttvg484RngtWI+R+P7N4YgL7grGZ8vInVEpAru8GIw/Ml+GilaDTYRUcWgOM29QK5nEHQF7vh3RcYCaZ6wzxTgFk94KhiDPHH/X4CbgAuNMTuMMXtwewnLcCuBHysc9zZur2YW/pmOO/vHDlayzwf+B6wARuMe//BHTeATEVkCfA08YLNfgGeBYSKST3nLvC+wzBN6aoN7neddwHzPYLzV4DOeMNDoIP3eh3tBp6WEuNKhZ+2PZ4AfcF+j9bjDaYFYApSIyGKfweeuuO+TEgO0uqqS8ojIX3Fb6Y8F2Odr3APWe2InWeVARI4zxuz3eAwf4i5z/2Gw43yOrwbMw73yYHG05FSOoopBSWlE5EPcyxpeaIzZGWC/s4FCz4Cs4iCelNxuuEN+s4D7TAg/PJ7022xjzJfRkVCpiCoGRVEUpRw6xqAoiqKUQxWDoiiKUg5VDIqiKEo5VDEoiqIo5VDFoCiKopTj/wNzn1ai6Bb1eAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fitted_cab_model3 = LinearRegression().fit(expanded_train, y_train)\n",
    "print(\"fitting expanded_train:\", expanded_train)\n",
    "plot_cabs(fitted_cab_model3, transformer_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "**Questions**:\n",
    "1. Calculate the polynomial model's $R^2$ performance on the test set. \n",
    "2. Does the polynomial model improve on the purely linear model?\n",
    "3. Make a residual plot for the polynomial model. What does this plot tell us about the model?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test R-squared: 0.33412512570778774\n"
     ]
    }
   ],
   "source": [
    "# ANSWER 1\n",
    "expanded_test = transformer_3.fit_transform(X_test)\n",
    "print(\"Test R-squared:\", fitted_cab_model3.score(expanded_test, y_test))\n",
    "# NOTE 1: unlike statsmodels' r2_score() function, sklearn has a .score() function\n",
    "# NOTE 2: fit_transform() is a nifty function that transforms the data, then fits it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# ANSWER 2: does it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1372fc790>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZgdZZW439OdS9LNkk4kKmkIARRQDElDZBFxAJU4skVRouKOMs7o/ASZDAHRBEWJExnQcRvGdUbEIEIbQQdUFpURJKETwxZZE9KghKUjJJ2kl/P7o6o61dW136p7b3ef93nypO9WdWr7zvedVVQVwzAMw4ijqd4CGIZhGI2PKQvDMAwjEVMWhmEYRiKmLAzDMIxETFkYhmEYiZiyMAzDMBIxZWGUgojcJyLHRXx2nIhsLGg/t4nIR3L87m0i8oSIvCgiHUXIErIPFZFXlLHtjHLMdGWZEPH5hSLy7VrLFUXa81bkfWQkY8pinCMij4tIrzto/kVEvi8iu1W7XVU9RFVvK0DEsvgy8AlV3U1Vu6rdWF6llWH7u4jIEhF5SES2uNftuyIys9ptq+oXVTWPwr3NHdhnB96/3n3/uGplMxoHUxYGwCmquhswB+gALqizPLVgX+C+PD8UkeaCZUnDtcCpwHuAycBsYBXwxjrI4ufPwPu9FyLyEuBoYFPdJDJKwZSFMYSq/gW4CUdpACAiE0XkyyKyQUT+KiLfEpEW97M9ReQGEekRkedE5Hci0uR+9riIvMn9u8VdsTwvIvcDr/XvN2h2cL97ifv3FHcfm9zf3yAie4fJLyKvEJHbRWSziDwjIstDvjNRRF4EmoE1IvKI+/6r3Jlyj2tCOzUgzzdF5BcisgU4PrDNLwDHAl9zV2hf8338Jnc10CMiXxcR8f3uwyLygHtcN4nIvhHH9SbgzcBpqnq3qvar6mZV/bqqfid4vt3XS0Tkh4FNfVhEnhSRp0TkX6K+KyKvF5H/c2V+QkQ+GCaXy1XAAp8CfTdwPbDDt72JInKFu+8n3b8n+j5f6Mr0pIh8OHDskfefUVtMWRhDuIPw3wMP+95eChyIo0BeAbQDn3U/Ow/YCEwDXgZcCITVj1kMHOD+mwd8IINYTcD3cFYCM4Be4GsR3/08cDMwBdgb+I/gF1R1u7uKApitqgeISAX4ufvblwL/DFwlIgf5fvoe4AvA7sDvA9v8NPA7dpq1PuH7+GQc5XgocAbO8SMip+Gcr7fjnL/fAVdHHNebgD+q6hMRn6fleOCVwInA+X7l4uEqrF/inLtpONd9dcw2nwTud7cJzirjvwPf+TRwlLut2cARwEXu/t4C/AuOMnwlzrH6ibv/jBpiysIA6BSRF4AngKdxBnfcWfDZwLmq+pyqvgB8EXiX+7s+YC9gX1XtU9XfaXixsTOAL7jbeAL4alrBVPVZVf2pqm519/8F4O8ivt6Ho1Smq+o2Vf19xPeCHAXsBixV1R2qegtwA84s2eNnqnqHqg6q6ra08rvb7FHVDcCt7Fy1fQy4VFUfUNV+nPM6J2J18RLgqQz7jOJiVd2iqmtxFPC7Q77zHuDXqnq1e02fVdU4ZQGOcni/iBwMtKnqHwKfnwl8TlWfVtVNwMXA+9zPzgC+p6r3quoWYIn3oxT3n1FDTFkYAPNVdXfgOOBgYE/3/WlAK7DKNUn0AP/rvg+wDGcVcrOIPCoiiyK2Px1HEXmsTyuYiLSKyH+KyHoR+RvwW6Atwm/wr4AAf3RNSR8O+U6kfKo6GJCx3fc676z+L76/t+IoJXCU2ld85/U5HNnbGcmzOEq5WoLXYHrId/YBHsm43euAE4BPAP8T8vl0hl9z/77j7o2k+8+oIaYsjCFU9Xbg+ziRQgDP4Jh9DlHVNvffZM+Mo6ovqOp5qro/jvP1UyIS5nB9CmcQ8pgR+HwrzqDg8XLf3+cBBwFHquoewBvc94UAqvoXVf2oqk4H/gH4hqQLXX0S2Mfzt/hk7PZvPmEbWcs3PwH8g++8tqlqi6r+X8h3fw0cEeWrcdlC9Dn0CF6DJyPkOiBB9mGo6lYc09U/Eq4snsRRjmH7jrs3Yu8/o7aYsjCCXAG8WURmuzPt/wIuF5GXAohIu4h4dveTXaeyAJuBAWAwZJvXABe4zuq9cXwCflYD7xGRZteG7Tcz7Y4zYPSIyFRcE1kYIvJO34D6PM4AHiZPkLtwFNa/ikhFnJDPU4Afp/itx1+B/TN8/1s45+QQABGZLCLvDPuiqv4a+BVwvYgcLiITRGR3EfmYb/W0GniXK/9c4B0hm/qMu1I7BPgQMCIAAMdh/SYROcPdz0tEZE7I94JcCPydqj4e8tnVwEUiMk1E9sTxOXgO9WuAD4rIq0WkFd/1Tbr/jNpiysIYhmtT/m92OhHPxzE13emagX6NM9MHxyH5a+BF4A/AN1T11pDNXoxjXngMx4kcnH1+Emdw7sGxb3f6PrsCaMGZZd6JY4aI4rXAXeJEO60APqmqjyYcMqq6w93/37v7+QbwflV9MOm3Pr4CvMONbEr0yajq9cCXgB+75/Ved/9RvAP4Bc4Av9n9/lyc8w/wGZwVwfM45/tHIdu4Heda/gb4sqreHCLXBuCtOCu653CU0Ozg90J+92SMj+gSYCXwJ2AtcI/7Hqr6S5xrfIsr2y2B38bdf0YNEWt+ZBiGYSRhKwvDMAwjEVMWhmEYRiKmLAzDMIxETFkYhmEYiYSWLB5t7Lnnnjpz5sx6i2EYhjGqWLVq1TOqmirJcUwoi5kzZ7Jy5cp6i2EYhjGqEJHU1RTMDGUYhmEkYsrCMAzDSMSUhWEYhpHImPBZGIYxPujr62Pjxo1s25alSrwxadIk9t57byqVSu5tmLIwDGPUsHHjRnbffXdmzpyJr+mgEYOq8uyzz7Jx40b222+/3NsxZVEQnV3dLLtpHU/29DK9rYWF8w5ifkdYawLDMPKybds2UxQZERFe8pKXsGlTdW3RTVkUQGdXNxdct5bevgEAunt6ueC6tQCmMAyjYExRZKeIc2YO7gJYdtO6IUXh0ds3wLKb1tVJIsMwjGIxZVEAT/b0ZnrfMIzRS3NzM3PmzOE1r3kNp5xyCj09Pbm285GPfIT7779/xPvf//73+cQnPpFbvt12K6eRoCmLApje1pLpfcMwRi8tLS2sXr2ae++9l6lTp/L1r38913a+/e1v8+pXv7pg6crDlEUBLJx3EC2V5mHvtVSaWTjPGnoZRj3p7OrmmKW3sN+iGzlm6S10dnUn/ygDRx99NN3dO7e5bNkyXvva13LooYeyeLHTIXbLli2cdNJJzJ49m9e85jUsX+50sz3uuOOGyhR973vf48ADD+SII47gjjvuGNreBz/4Qa699tqh196q4cUXX+SNb3wjhx12GLNmzeJnP/tZoccVhjm4C8BzYls0lGE0DmUHngwMDPCb3/yGs846C4Cbb76Zhx56iD/+8Y+oKqeeeiq//e1v2bRpE9OnT+fGG28EYPPmzcO289RTT7F48WJWrVrF5MmTOf744+no6Ijd96RJk7j++uvZY489eOaZZzjqqKM49dRTS3X+m7IoiPkd7aYcDKOBiAs8qeZZ7e3tZc6cOXR3d/OqV72KN7/5zYCjLG6++eahgf7FF1/koYce4thjj+W8887j/PPP5+STT+bYY48dtr277rqL4447jmnTnOKvCxYs4M9//nOsDKrKhRdeyG9/+1uampro7u7mr3/9Ky9/+ctzH1cSZoYyDGNMUlbgieezWL9+Pao65LNQVS644AJWr17N6tWrefjhhznrrLM48MADueeee5g1axYXXXQRn/vc51Lva8KECQwODgIwODjIjh07ALjqqqvYtGkTq1atYvXq1bzsZS8rPavdlMUYoGy7rGGMRsoOPGltbeWrX/0ql112Gf39/cybN4/vfve7vPjiiwB0d3fz9NNP8+STT9La2sp73/teFi5cyD333DNsO0ceeSS33347zz77LH19ffzkJz8Z+mzmzJmsWrUKgBUrVtDX1wc4pqyXvvSlVCoVbr31VtavT11pPDdmhhrlWEKgYYSzcN5Bw54NKD7wpKOjg0MPPZSrr76a973vfTzwwAMcffTRgOOM/uEPf8jDDz/MwoULaWpqolKp8M1vfnPYNvbaay+WLFnC0UcfTVtbG3PmzBn67KMf/SinnXYas2fP5i1veQu77rorAGeeeSannHIKs2bNYu7cuRx88MGFHVMUoqql7yRy5yLnAh8BFFgLfAjYC/gx8BJgFfA+Vd0Rt525c+fqeG1+dMzSW+gOWVa3t7Vwx6IT6iCRYZTHAw88wKte9arU37cyPDsJO3ciskpV56b5fd1WFiLSDvw/4NWq2isi1wDvAt4KXK6qPxaRbwFnAd+M2dS4xhICDSMaCzwpjnr7LCYALSIyAWgFngJOALzA4h8A8+sk26jAEgINw6gFdVMWqtoNfBnYgKMkNuOYnXpUtd/92kYgdFogImeLyEoRWVltNcXRjCUEGuONeprORytFnLO6KQsRmQKcBuwHTAd2Bd6S9veqeqWqzlXVuV588lgkKdJpfkc7l759Fu1tLQiOr+LSt8+ypbcxJpk0aRLPPvusKYwMeP0sJk2aVNV26hkN9SbgMVXdBCAi1wHHAG0iMsFdXewNjNs40LSRTmaXNcYLe++9Nxs3bqy6N8N4w+uUVw31VBYbgKNEpBXoBd4IrARuBd6BExH1AaD8oicNSlkZqIYxWqlUKlV1ezPyU0+fxV04jux7cMJmm4ArgfOBT4nIwzjhs9+pl4z1xiKdDMNoFOqalKeqi4HFgbcfBY6ogzgNx/S2ltAcCot0Mgyj1tQ7dNaIwSKdDMNoFKzcRwNjpc8Nw2gUEpWFiLQA21RVReQA4CDgZl8uhFEiFulkGEYjkMYM9TucLOu9gFuAjwLfLVUqwzAMo6FIoyyaVHUrcDrwTVV9G3BouWIZhmEYjUQqZSEirwXOBG5w32uO+b5hGIYxxkijLD4FXAzcoKr3isj+OKYpwzAMY5yQ6OBW1VuAW0SkVURaVPVR4J/KF80wDMNoFBJXFiJymIh0AX8GHhaRVSLSUb5ohmEYRqOQxgz1PeBTqrq3qrYD5wHfL1UqwzAMo6FIoywGVfVW74Wq3gYMliaRYRiG0XCkyeC+TUS+DlyN0yt7AY4P41AAVf1TifIZhmEYDUAaZeE18w7mVhyBozzeUKhEhmEYRsORJhrq2FoIYhiGYTQuaWpDXRj2vqp+sXhxDMMwjEYkjRnK36ptEnAScF854hiGYRiNSBoz1Jf8r0XkS8D/liaRYRiG0XDkaX40Eaiu87dhGIYxqkjjs+jCiXoCp4DgXsClZQrVyHR2dVszIsMwxh1pfBbv8P3dD/xFVbeXJE9D09nVzQXXraW3z3HjdPf0csF1awFMYRiGMaZJY4Y6VlUfcf+tV9XtInJJ6ZI1IMtuWjekKDx6+wZYdtO6OklkGIZRG9KsLN4tIr2quhxARL4KTC5XrMbkyZ7eTO8bhmGMFdKsLN4GnC0i7xSR7wLNqvqBInYuIm0icq2IPCgiD4jI0SIyVUR+JSIPuf9PKWJfRTC9rSXT+4ZhGGOFSGUhInuIyB44q48PABcCW4EL3PeL4CvA/6rqwcBs4AFgEfAbVX0l8Bv3dUOwcN5BtFSGNwlsqTSzcN5BdZLIMAyjNsSZoe7DiYIS3/+nuf8UmFHNjkVkMk5dqQ8CqOoOYIeInAYc537tB8BtwPnV7KsoPCe2RUMZhjHeEFVN/lYZOxaZA1wJ3I+zqlgFfBLoVtU29zsCPO+9Dvz+bOBsgBkzZhy+fv36WoluGIYxJhCRVao6N/mb6TrlfUxE2nyvp7gDdbVMAA4DvqmqHcAWAiYndTRZqDZT1StVda6qzp02bVoB4hiGYRhRpHFwf0xVe7wXqvo88I8F7HsjsFFV73JfX4ujPP4qInsBuP8/XcC+DMMwjCpIoyyGeXRFpAmoVLtjVf0L8ISIeN7hN+KYpFbgONRx//9ZtfsyDMMwqiNNnsWvRORq4Fvu648Bvy5o//8MXCUiuwCPAh/CUWDXiMhZwHrgjIL2ZRiGYeQkjbJYCPwTcK77+lfAfxaxc1Vdzc5OfH7eWMT2DcMwjGJIU6J8APgP959hGIYxDolUFiJytaq+O1B1dghVPaxUyQzDMIyGIW5lsdD9/x0x3zEMwzDGAZHKQlU3isjJwCuAtar6m9qJZRiGYTQScbWh/gMnSa4d+DcRubBmUhmGYRgNRZwZ6nhgjqr2i8iuwO3AF2sjlmEYhtFIxCXl7VDVfgBV3ZLwXcMwDGMME7eyOFhE7nH/FuAg97XglG2yaCjDMIxxQpyymFUzKQzDMIyGJi4a6pFaCmIYhmE0LuaHMAzDMBIxZWEYhmEkYsrCMAzDSCRPbSiLhjIMwxhnWG0owzAMI5HY2lDed1R1nf8zETkWsGgpwzCMcUIan8V1IvIpABGZKCKXA18uVyzDMAyjkUijLI4EDhSR3wN3A88BrytVKsMwDKOhSNNWdRvwPDAZaAUecLvnGYZhGOOENMribuCXwOHANOBKEXmHqr6rVMkamM6ubpbdtI4ne3qZ3tbCwnkHAYx4b35He50lNQzDKIY0yuJjqnqX+3c3cJKIfKhEmRqazq5uLrhuLb19zuKqu6eXhdeuAYW+QR1674Lr1gKYwjAMY0yQqCw8RSEiU4FJ7ts3lSlUI7PspnVDisKjb2BEi3J6+wZYdtM6UxaGUSBhq3p7xmpDooNbRN4qIn8GNgJ3uf/fUrZgjcqTPb2lfNcwjHi8VX13Ty/KzhV8Z1d3vUUbF6SJhvoicAywTlX3AeYBvytKABFpFpEuEbnBfb2fiNwlIg+LyHIR2aWofRXB9LaWUr5rGEY8Yat6bwVvlE8aZdGvqpuAJhERVf0VcESBMnwSeMD3+kvA5ar6CpworLMK3FfVLJx3EC2V5mHvVZqFSpMMe6+l0jzk+C6bzq5ujll6C/stupFjlt5iMy1jTBK1UrcVfG1I4+DeLCK7Ab8H/ltEngYKuToisjdwEvAF4FMiIsAJwHvcr/wAWAJ8s4j9FYFnH22UaKgwh7s5142xyPS2FrpDFMNYXME3om9GVEc6Z4d9QWR3YCvOKuT9OPkW/+OuNqrbuci1wKXA7sC/AB8E7nRXFYjIPsAvVfU1Ib89GzgbYMaMGYevX7++WnFGJccsvSX0AWpva+GORSfUQSLDKIfgxAicFfylb59V94G0SKKO8/TD27n1wU2FKhARWaWqc9N8N67q7M2qeqKqvuC+NQB8pyrJhm//ZOBpVV0lIsdl/b2qXglcCTB37tx4jVcy9ZwF2NLcGC9ErerrpSjKeu6jfDNX3blhqPx3PSwIcWaoaSXv+xjgVBF5K05I7h7AV4A2EZmgqv3A3ji5HQ1Lvc1A42lpbhjzO9obYhVR5HMfVDphzzMM7xMBtQ/Pj1MWk0Xk7VEfqup11exYVS8ALgBwVxb/oqpnishPcMqi/xj4APCzavaThmpmCHERGt42ylx5LJx3UOiStVbOdcMYj6R57tMQpnSEkYohilpaEGKVBXAyTrOjIApUpSxiOB/4sYhcAnRRoOkrjIs611a1vEsyA1UzA0mjZBptaV5LGtEJaIwPijL/hikdhREKI0qB1NKCEKcs1qvqh2shhKreBtzm/v0oxYbmRtLZ1T1MUXj4ZwhJA1KSGSjvDCSLkmmUpXktqbf5zxj7xD37RZl/o5SL4gSpePs+/uBp/HRVd10tCHF5FmErijHFspvWRS73nuzpDc0YPXf5ai7qXDv0vbC8C/9FzDsDsQSkeOz8GGWSlC2e9NynJUq5eNGMjy09iTsWncAl82dx6dtn0d7Wgrif1zoKLG5l8b6aSVEn4gbs6W0tkUvEq+7cwNx9pw6b0Rc9A4mSrbunl/0W3TjuzS4WBWaURWdXN+dds4aBQFqB3yJQlPk3yefYSKbWuLaq99ZSkHoQNZALzkU8d/nq0N8pDDMjxZmB8jqgk6IixrvZxaLAxh6NMDB6K4qgovDwT0aKMP/6lU53Ty/NIkNKaeX654aZnur9zKcp9zFmCVtKCnDmUTOY39EeO/CkncHO72jPtXwMky2Id1ONx3IfRZkBjMagUYoEhlkT/BQ9GfEUpBcF5Smp7p5errpzQ0OZWuOS8n6jqm8UkS+p6vm1FKpWJC0lvdVFtVEIeWYgQdmifCueH6WeyTr1YDxHgY1FigpFrZa4SWDUZCTviigYpBF8xuP8qfUgzmexl4i8Didx7scEHN6qek+pktWIsIHcf/Fbd2lmy47hN3HSDDbq5sl6U3myRdlQPeqdrFMvxmMU2FilKB+Uf6beLMKAKu0ZBvBI07TAxAlNnLt8NctuWjfsmc4SlecfA5pc+bJSL1NrnLL4LPAZnCzqfw98pjgF/8YcwYu/ZccAlWZh110msLm3L3GQj7p58tofk2yoUZij1xhNRA3Sba0Vjll6S6oJVvDZ85t00q62w3yMlSYBgZ7evhHby7IiipIvjmB+RT1NrXEO7muBa0XkM6r6+RrKVFeiOuHtOnECqxefmOv3vX0DXH3XE7HRFRC+IkmyoUYxXhy9jeAUNaondJBuFl7c1s/zW0cO0mHXOO5ZSbvaDjNvbt2xU4bg9rKsiLI+y2UVD8xLmraqnxeRU4E3uG/dpqo3lCtW/Uhz8eMGqKjfJ0VXhK1IzomIxkpivDh6LTFv7BA2SG/Z3j80m/eIG/STVtNZglL8299v0Y2R28sSlZdltZ/FdFYrEpWFiFyKk1F9lfvWJ0Xkdap6YamS1Ymki580QEX9vjnCPhmX6R1HVDmARrzJiqKzq5uLf37f0CwvrATCePHXjCWCk6/LF8xhfkd77CAdRly4ufd5HuLGhCyh8UnyeSS1F6jXajpN6OxJwJtV9buq+l3gLTg1o8YUXvipF8Lmx3/xkzKHo0I6333kPrkyvcNoqTRz5lEzhoXjXr5gDo+72Z5jcaC8qHMt5yxfPcwc0GjRIkZ24kJmowb3qPfjws2rWW3HhWlnCY1PEw4P8fdvPUOM03TKA2gDnnP/nlySLHUjLIQtaqYel1k9c9GNNIuw/7RWHt20lQFVmkU4/fB2Lpk/i7n7Ts2c6e1HoO52y1rhnz1NbqmMMEfEMV78NWOBuMlX1oTWsAS3rNFQSdsNe3bTRuUFtxMVDRV3/9YzxDiNsrgU6BKRW3HGqzcAi0qVqsZElfUIWw4mDeoDqjz09JZhr6+6cwMAl8yPTsYLezD8jKfOd0HlnUVRFOmvMed5+cT5CPPk0pQVTp13u2H3kPccR3XEi7t/o8aeNOataknj4L5aRG4DXuu+db6q/qVUqWpMloiGpEE9jGA9qTC895esuG/E4DgeHNZFxJ8LcPrh0Q91lsHfnOe1IclHOJpzaZLuoTzKMMr32Szl131NVe5DVZ9S1RXuvzGlKCB62Rf2ftBGmRYFzrtmTWxJjvkd7axefCJXLJhT1+qStSZoh82jKMA5x7c+GN4aPqut16ra1oaxXLYlzT00v6OdOxadwOUL5gBw7vLVsSV7op6NvM9MFtL6LMYkwbosweSX4w+eFpoQ5J8V7LfoxtRdrdImCY3m2VQe8uaShBG1Ssxq67WqtrVhLJdtSXsPZVnFtkesxNpr4Kcbt8oizqk9pbXCtr4Bfuj6GiD6AoaVA0mDhXjupMgBOGqVmHXwt6q2taNIf0C1z1OR20x7D2WZyNSzjXKsGUpEmkXkwdKlqANRTm1HUQzS2zc44jdhZoitORSFx3ifpXrhykUtoOMemrSmxrQh1EZ9SdOYLGs15os613Lu8tWFhaWmNbHFRVgG9523inURxK4sVHVARNaJyAxV3RD33dFG1AUKpvUn/S5NyGuUU6qttZIg5dils6ubhT9ZQ99gMaqitdLEF2MemjQzsiwh1EZtCc74t2zvj21MBmQu8JfUYjkraU1scWNIUOZ6RuiJJjhGROS3QAfwR2AoJlRVTy1XtPTMnTtXV65cmek33uwxK8EQ1rDwNz8CXL5gDguvXUPfwPBzXWkSlr1z9rgZgPw3OkQn1eWhpdI8bIYV9lBB/IMbdU+Mp7DlOOo1UCU9Y0E8+32Wa5k0HhSV4xR1X6YJm48Kta1mZSEiq1R1bprvpvFZfCaXFA3O8QdPG+aTSEPYEtK7SFElxKe3tTC/oz00JLZvUMeN3yLrA58Vv4nQXxIEnEFj4U/WsOyds2MHfXNqO6QZ0GoZSpw1ACLuekUphKRrXER3yos61w5bvXjbu/TtTn/tqFpwnmz17vmRGDqrqrcDjwMV9++7gVHfyyIqxDKKKa2VSA0+v6Ody86YHWuf3ByRWDZeBqIiI56i8IovhpkS+waVJSvui/191vISY5GoEOOLf35f3UKJsz4j09taIq+ZQKgPIu01znvMacxcURFNnmz1nswkKgsR+ShwLfCf7lvtQGe1OxaRfUTkVhG5X0TuE5FPuu9PFZFfichD7v9Tqt1XGGlPcHtbC1csmEPXZ09MzByNczyN94Eo7fme0lrJlL+ShaRM8LEc85+WqNlrlC+vFgNV1DOy6y7NkUEIC+cdFHofKYQO9mnrNkG+Y15207rEWmZJ91+9x5A0SXkfB44B/gagqg8BLy1g3/3Aear6auAo4OMi8mqcUiK/UdVXAr+hpNIiaU6wZytMu8TzEmweCynol3UgGmt9tdOc70qzsPiUQ2jdJd1DWzT1jDRpFPLM4suks6ubLdv7Qz/bumOA1x0wNfR6ze9oz1Rocn5HO6cf3p5qopLnmOPOqz9bPe7+q/dkJo3PYruq7hA3nVxEJlCAb1JVnwKecv9+QUQewFm1nAYc537tB8BtQOE9wJPKdhR9EbIkH6VN0hlNtYuiOpDtNmkCPVudDoTHHzyNZTety5W3koamFCPBeEuIDBLZsa6lwvb+wdLi+/M4fhW445HneO9RM7hk/qwRn0clsIUN9p1d3Vx91xOJA1veY45s1wrDthd3/9U7gTFNNNS/AT3A+4F/Bv4JuF9VP12YECIzgd8CrwE2qGqb+74Az3uvA785GzgbYMaMGYevX78+836Dvba37hhAcUJd333kPqE3YNTvoy5cngE9TX/DjvsAACAASURBVFROlMO4raXCklMPacgBL+5clO0A93h86Umlbr9RyDuRiIu4geED1fEHTyuki1vUPidVmhJD2WFnpFKw0uzxB08b1srYfyzBSVfSvVdtNFTYPgQ4M0LR1Yos0VBplEUTcBZwIs7x3QR8W5N+mBIR2Q24HfiCql4nIj1+5SAiz6tqrN8iT+isnzwhaWl+k/Y7wYf63OWrQ2c4AjzmDnZxoX7VhtPVg7yhzFkoMwS2kVZ51YZYpp0EFRXGWda1T9uWNGn/Rd03RdwjRd9nhSoLd4O7AAfjrPzWqeqO3NIN324FuAG4SVX/3X1vHXCcqj4lInvhtHGNXfflVRb+2lBhxN0kaWb/Sd+JeuAmTmgKdcZ6iWFxMqeRvRHJUmMrD2X2My4j/r0aapEvUuQ+yrz2aeSJ23/Z17GaSshFyJdFWaSJhjoJeAT4KvA14GER+ftckg3frgDfAR7wFIXLCuAD7t8fAH5W7b7C8IcIRhHnlEqTop8U6hYVeSJCqCPr+IOnJcqcRvZGpAxHqd9ReNiMyVx154ZSOow1WoXaWoRYZt1HVMBGZ1c3TSWW105zzFH3XrNI6YpiNFVCTuPgvgw4XlUfBhCRA4AbgV9Wue9jgPcBa0XEy0a5EFgKXCMiZwHrgTOq3E8oaeL+4wawNCn6SYXEom7knq19XL5gzogZR5ZchdEQkuufVbW1Vqg0ybDyH81NwkDOciBNwtAsrbOrO9S0V1RCU73j34NkKYKY16yRZh9RVZ29QXHl+uf46aru1OW1m0VoEiWkbFusnElElYLxFEVZJsYlK6JzVxqxEnIaZfGCpyhcHgVeqHbHqvp7iIxUe2O1208iaXaeFPUQF02Vti1k3AMXFhVxbkSGZ1bZq6Uo26v/3AQdmW0tTt2sLF3y/Ayqc77OWb6aZpFS+3U3SoXapJL7wXuimgZPSfd2WJ0tP719A1x91xOZ+jAMqnLZGXNS1xRL+xzERRlFnaOV65+ryqTZ2dUdeW83aiXkSGUhIm93/1wpIr8ArsG55u/EyeIe1UQV94N0ReO8z+JS9JNC3Y4/eNqIrM6kyqlRYY27TpxQE+dqUR3kklZJm3v7qrZje7+PG5CKeNDqWTbaI08RxGrKRyTd22lWwVkb9gxdq5RWqzSh0h5RIatR5yisbIe3nTTEmY6i7sl632dxK4tTfH//Ffg79+9NQOPbOBKIu1HTOujmd7RHOpuT2kJ2dnXz01XdwwbEpLagUbkKIoyITw9r2lQERdWnSVOLJ4qwmP88BGPc81Lv+Hdv32n7yHtUa9aIywlIs42oCdvONgEjB8VlN60bUZAzii07Blh47ZohWf2kXR1HHUe1Js248xN1T9b7PotUFqr6oZpIUCeK6jiVV9tHPdw3rHkqcnkbvFkmt1TYsqN/yITT3dPrPBzK0DLdK6J38c/vG0p+q+YGK8pumqa0exgtlWaWnHoIEF28MS1nHjUjNvIky0NZ72S+uOsSdSxFmjWC+5jcUok1IXrRacE8iEqzoOoMvv6cCU/mtKZYj76BkcU6s6yOs9ynWZ6BqO1Oaa007H2WJhpqPxH5dxG5TkRWeP9qIVyZFJU6n7dERKRzu7cvNDrCiybxHpbLF8xh14kTRsyy+gZ0hD23b1B5fmtfIVFARdWnyVKLxyNYzmGwylSfG//0VOh5CItSOWf5ajo+d3PDll2JOv+TWyqRETdh10Dc7xxwwS+Y6UYuXdS5Nrb0TNj52rKjn0rADuS98q7jJfNnDXt2prRWQHf6qQZUh57JpBprcQSftSxRRVHnKIwsskWNP4tPOST1NmpNmqS8NTghrmuBoTgEtwJtQ1BtnkU9lnRZEpFaKk2AjFi9VGOGyRtzX2Ssd2dXd+rVQZi8RSZzTWmtoOr4Sppi/FmNmvCYNQvan+uTJm/HT/AcRF2HJoE9JlXY3JtuRVtN5YI4gvfOzEU3Rn43LMM/OE6kzQxPImuORRljVdH9LLap6lerlKkhqeeSLqk2lZ+oFq9xTvok8kYBFWk3nd/Rzsr1zyX2FYnyLSycd1BkgEFW/ANq3DnNYpuu5WQkeF3aXOWXVC3W+11U1YAwgucg6l4aVNjeP8jlC+akOu40Jk7vnkkbSVVplhH3TtRz0xyR7xE2Tszdd2rV1zbt+FNUUEm1pFEWXxGRxcDNwHbvTVUd9T0t6knYoLvV539Ig7dMz7PCqCYKqBol29nVPaw5UZqAldcdMDV0f/M72kc0OqoFaRRtPR5w77qkmX37r39c+ewo/Ocgzq6fRblGbadJhIs613Lrg5tCw4KjmNJaYfEph4zIlYj6bZaJVy0nmvVueuSRRlnMwkmeO4GdZih1X49q6l3PJ3jDdXZ1Z5op+8t/eMeQxpxQrx4NnV3dI9rLpnk879mwmc6u7tB6WmHJfGUzuaUyItoMhiv+sB7RaR7wIu7JpLDV4PXPs8pUHLOR1zsiTjml3X7UdgZUh60+01xpvzkprekqa3BLrah3Mp5HGmXxTmD/oupBNQqNsrTzk2Wm7Hf8+eWNsvs2izCoWpM8jLjY+7Rhj378g2xYMl+lWWhLiL4pikqTsGVH/9C+POe3n7wlZIq6J+P2EZZzkTcyLdgWNK61cBo8maqNcvO64WXJ+ahmAlX2pLPeyXgeaZof3QuMKBE+2ql3nZUokqIhghFXwZo7xx88LTTK4rIzZnP5gjmAY58uo6FSUq2bapzRcfW0+gaUXSdOKK3DnkdbS4XdJo2MQMu0jdbKsOvljzQ675o1offkORmvV9QgEtXMK09kml++JSvuY35HcmvhNBQR5Rbshpc0A49rmZxE1vpOeah30yOPNCuLNuBBEbmb4T6LU0uTqgY0ytIuSNzqIhjVETYT/emq7tDqqkDpK6kiFLAIhI0VafoQ550hp2X14hPZLyaSJonmJuHFbcPzYvzmlbjZdJGlOIJUO6Pv6e0bMZPPM8v2ZuhFGBS7e3qHzGRJ90XrLhNyRyHVwp9Q72Q8jzTKYnHpUtSBRlnahbH4lENSPexRN+qtD24KDTMt+6aOG8jTKozJk+I7ssVdtywRZlnx7NnVKKSBQaUayfKW4vAio85dvnqoZllwG97rvOfPkyuv47eM5leegj1sxuRcpsE0ZsFaTDrr7Vv1SDRDqertYf9qIVyZNMrSLoywRL/TD3dKi/gTo7LcqLW4qeMS9tLuZ3NvX+4+xN55K6PitXdfVGOyKYIspTjuWHQCly+Yw7a+QXp6k5Myg/ddFrp7ejP3i/ebUMNMcGmJk7W3b4A7H30+9veT3aKVQdKslItKUo2iFmautKTJ4H5BRP7m/tsmIgMi8rdaCFcmeTOvk4iq259HvjsWncBjS09i4byD+Omq7hE3TNRN3iQyYr9l39QQP5Cn3Y9nb1447yAeW3rSCBt70nWb39FeQIf4aLz9t0Wc+7LJer2ymgb9913W6KAsg1lwEKzGoe3VwIoiadtbdvSHyptmgpVn0plljGgk32qiGUpVd/f+dhsWnQYcVaZQtaLoWOmyIqyibphJlabQPIsB1RH7rUXFyjjbaprkO4+k85Z03aJMRWnj88Pwy+PPZ0iK3c9DS6WJiROaR0R35bleaQa8KDNHVOHKAVXiIpV7+wY475o1Q2W8g72xs/ZmScLz5cVFAsYpjGD9qCTfiV9hB5MEm0Vii4GGjRELr13DkhX3hWa6N5JvNVVb1RE/EulS1Y4S5MlFtT24i6KsdpZRbR8Fp0ZUVPZtmEM8WLagjDajYeQpzVF0SZKg4//4g6dl6qmQp81uEk34aujE4E8wy0Ke1r7gRH55BRv998zMl7RwxyPPZZIhSDVKO0iwSVHUdf/RXRtiFZxHkytc1DUJlvVIU/7G/9zFlZIJ+33ZLXILLffh62sBzr09F9iWU7YxTVmzgCRnfNrGPv4Zedl5Jp1d3SxZcV9VuQ+1KkmSZcUTTA6MO8ZKszChSULLtcBOBZAmt+Zvvf3DXqdV/EkryqgZfk9v31AOhTcoXdS5NvW5iqMoReHPGfHOR1S12hv/9FSq/KUkhTJxwnDLfVI0VPA5SzMx8f++3j0s/KSJhvL3tegHHscxRY16io4yKCvCKu6GydNEBcoJ+ctTlC6OskqSVCPnOW73vSSaRegbUAZjlg3e4NW6y4TEgcxvWoSRYdD+QTxM8Ufd53EKOTjoXVWAoigKgcjigmHVansKKgnjKVFwzm3SBDGsdWoagrW7GiEaKpcZqtHIY4Yqsnpqmdv0bzvshokyUQFcEVPALe53AplvympCH6Oa3RRd3fWizrVcddeG0DyOepG1VElbS4UXtvXnrtTrJ81KQYDHlp5UaIXfIvBXJNiyvT90dTeltULrLhNSm3+ykOQn8VY1eQtdFmVmSiKLGSpSWYjIZ2N+p6r6+TzClUEeZVGWLbDWMdFxD3Fce9g0D3+WATvvYOL5XaDY2VPwOhRhay+LaqoHJyE4oaEiDGt+lTbgwGvZ20iKohHwlGjYJMnzyVR7XdO0d66WopTFeSFv7wqcBbxEVXfLL2Kx5FEWcU7jx0Jq2jcqaWb0Yc7RLMXV0ijPuJVKEmE9BKqhjASvsqm2P0nWfW3rG0i8XpUmAaGq8ibV4PXD6OntGxp4y1SsWfCvWjxl/PzWvkKd91B+/5QsyiIyz0JVL/P+AVfi9N3+EPBjYP9CJK0jtcg7qAX+vIMont/aNyL2fX5HO6cf3h5Zw98jrZM573kro9JnkWGZtaCtpTKiY5zT8KocelMoCsGx/ddLUcDOfhhXLJjDI5e+lceXnlR13agiqDQ7ZVu8HJGe3j629Q3S1lJJpSiyJDz29g1w8c/vKyR3q1pi70gRmSoilwB/wnGGH6aq56vq0zWRrkQaOYM7K14iVRzBRJ7Orm5+uqo7cZaWVgnkyWwu63zXu75XVkRGZlwHI6imtFactqM1QomPDCojSz6MtBnTU1orwxI1y0qabG9rYdddJozwM/X2DaSO/FPIdC2f3xrearnWRCoLEVkG3A28AMxS1SWqGp83XyAi8hYRWSciD4vIoqK3X1YGdyPjH0TTzr63RmS3BvHOZ1q82XTwfBeRAV+P1WE1C4Hnt/YNHWvUdWndZQKLTzmESvPwUbrSLLz3qBn5d56TWk7wvaKAUX3Dvd7VXub5HYtOYMmph5RSlmXhvIPYXGUp/Pa2Fhafkl++RszgPg+nyuxFwKdl51RCcBzce5QllIg0A18H3gxsBO4WkRWqen+R+6llt6s05HWOe79Lwj+Ipp19P7+1j4XXrgHS5V9E2WzT9NMIy/041w1VjXP2pemRXDYRqRSp8WaMiQ2EgidXnRafN6x5qib9POqFv3fG6Ye3J2ZMB0NOi9JtXpmdqOirYFRfkGBIb95Oj6Mmg7tsRORoYImqznNfXwCgqpeGfX/33XfXww8/vCayPfPidp54rpft/QNMnNDMPlNb2HO3iYVs99FNW4bZZJtE2H/arrHbD/tdGMFtdW3oYXt/+sF0QnMTc/edEvudlY8/T39EYkFzk7DfnvHHkiRT2PmIOm/Tdp/IX/82dnJHJzQ10dwkoefHuw/T3AejnQnNTQwOaurnxP+81kIGYNj40NZaoWdrX+x44ZdxQlMTCPQPDDJxQjMDgxr6TE2c0EzHjOrbDN1+++3FZXDXiXbgCd/rjcCR/i+IyNnA2QATJ1Y/WIcRVAxtrRU2vbB96CbZ3j/Ao5u2AFStMJ54rnfEgz6oyhPP9cZuO+x3HoKgaOhNus/UFh7ZtIW0k4X+AeeGfebF7Tz+7Nah1xOampi5Z6vznZgMtIFB5ZGYc/XMi9sTH+iw8xF13nq29jFxQnOhg0QeioqOGVClvz/8/G7vH+Dhp19kQnMTTSr0Dw4OXfPHn9kae11GG95952dQlUeeHnlvpZ1I5ZHhFS/dLXLSmHYsCI4vr3jpbqGKJGwytM/U2ptaG1VZJKKqV+JEaTF37ly97bbbCtt2Z1c3F//8Pl7c2od/Li3AS0O+P7mthduqTKCZGdNUp88t8R1mwsmblAcw5+KbM5kuzlkwh4XXrmHPQITM9iZhlwlNvHxH8sAcdq4889PLU5qN/OdjSsRy3MvfqGcIbdFhnlGNofwEa2BNS/GbscJApZlzfH6wY5bewktLMNc0i3DJGbNTmUSTTK5TfPdmUP6s28yDZIhUGBNmqCILCeaJ0S8iN+OAC36RamAJxl0nDfhxBeiy5EYUmZzVHlB+RZYICe4na7HAoig63t4jTcZ3WfsOY0prJZfNPS9Jx+bPC4qbgE1xzUPT21p4fst2toY4nZokOiIsLP8hSwWHsgsEpqWQPIs6czfwShHZT0R2Ad4FrKjFjvPE6BcRfZN2MAv2ZE6aGITlWHiklbvSJCw59ZDCnGr+MMCF164pLTvYazNbjySusvaYZrtl7bvSJExprQxFD16xYA5dnz2xpr09ko7NK/YIxOYQte4yYahXTFQuSVIZ9mBQSZbeE3mLjhbVLycPkWYoEXmB8GtTejSUqvaLyCeAm4Bm4Luqel9Z+/OTdLGCMxtheL/fvMvD9oztOpOiZ/xEFQg8/uBpoSUfjjlgKo8/2zti2Zs3ciOOspO+/FVIRwPejDdK2v4MdaSKJCoarbOrmy07+iN+VR+8Qn9JPc3BGeCz1Oby86SrmJIirsLGlKiIqqiGZhAdLbhy/XNcMj992HpeIpWFv+lRPVDVXwC/qPV+43os++3B3T29wxRHlhLfYTbIhfMOYuG1azINnlkGwrAb9tYHN4V+9/Fne0OXwqNkvB3BgCqVZqlrNnIavHyBvMXnqiXqXoozjSy7aV3DndfevgGWrLgv1mQlENuaOBVCqmvlX8F7z36U6XhH/wDHLL0l1D8RtnJR4Ko7NzB336mlpwGkNkOJyEtFZIb3r0yh6klUJrKXRHbJfKe+f3tby4gbMU2yTFRP3ZXrnxtxZzeJW58nBq8ccxJhJqeoB8W/lPdTbTJSvWhrqaSyzdQqKzkMLyl05fr6FTvcf1rriFIUSVn2jZot7/Ucj8Jr31uNCTnt5Km7p5eOz93MRZ1rh579KLb2DY4w03rPYtS59o6lbNL04D5VRB4CHgNux+ln8cuS5aobYZndVyyYw+rFJw7T3HltjlF2zR/euWHEcnhQYbdJE2LrJ3mDTBzBB96ze8bd62F+jiwPVlRGc3OC8iualkozIqQyNagyIkO6Fvh7M9SzZ8RDT28ZYWKNaxEK2f11zSJcsWAOj+fo8V003T29QxaCsnl+ax8/vHNDZn9o34AO+SjbYkqE1EJppwmd/TxOz+1fq2qHiBwPvLdcsepLmszuvI2Osl7Unq19dH32xNBIiybgqc29iUvhpDaQYfhXSZ7JLO5mDRKV0TxYsM39vUfN4Ko7N4wY5JSddvZzU5p1mkVY9o7ZNTcDKTvDtRvJoKMMN1X6m0Z5JqvWjHVOBtTpd33u8tWFFEtsa6mwvX9kL5RJlabU/jUlfVnxepkzk/yZtShxk0ZZ9KnqsyLSJCJNqnqriFxRumQNTt52h20ZQw29myBYvmBSpckpNpfyvs3TbyLoRC/CuV3kY9be1sIl82cxd9+psXHoaR3zA6rM72gvLZQ3jnr5KZLwgjeCJVS8ATUs5DTNNvP+1k9LpTm0T7jXryM4iYjDK+6XdJ+UpSiyBrj4qVUB1MQ8CxH5NTAfuBTYE3gaeK2qvq506VJSZJ6Fn6RkmDzJMlkS4YIzZP+20+ZlQP5eCY0cRVRpEhYcsU9o32k/nV3dLPzJmtQRL1e4jZjyJPO1x8Tsj3ZqmbsRRaVZ2HWXCWzu7ctUY8yjpdLEtr7Buh+HHy95dH5He+YkWYjPo0q1/yKaH/k2tiuwDee4zgQmA1ep6rO5pCuBMpRFkS1SL+pcW3ViWHDfcQlH44HmJqGJ4b6IsOuTdUXlbQOGz1a37uhPnHU2woA6VskyKCYlvDVai1iv+VfWiY1HNR31siiLRDOUqm7xvfxBZmlGKVGO6POuSa7A6l9xtO7SzJYUZTCSCOZKNPKsPystlSZ29CsDqqkH3IFBJXhWw/JJsvqIvG3cseiEzJ0Fx8bVaCyaRbgspLRGWKVhb5WZlO8QZkIum6hs8GYROru6h/lJs5pBs4TtV0OaaKgXRORv7r9tIjIgIn8rTaIGIWqQGVDl3OWrmRmRQRkMjS1CUYTJ9O4j9ylsu/Wmt29wSPEp1UVMBa9blOMvbg9h197rLDgeiMp8znJViugl0VJpjlQUwfDzH965Yeh1FH7/36Vvn1WzzHMB3nNkeLbBgOpQ5KFfAba3tXDMAVNTn/Na9LhIVBaquruq7uFmbLcApwPfKFWqBiAuusCfiHeOqzg6PndzbPOaomW6ZP4sjjlgain7KYP2tpbU0S8DVURMBa9bVLOcM4+aETkoRl37qCTG0YTTsjV6IG+pNPPuI/cJPWevSzl4eeHcSaGxbS2VYeVD3nvUjFTNyPI8Y2FO4Fplnk9va4m9d7wkwqACvGfDZs48akbqrnplh89mqjqrjoOjU0QWA4V3r2sksi5VvSZBZYbVeV3r5ne0c1HnWv7vkfolcGXBK4lSJJVmAR3usxCcEiZ+glFkfsfo3H2nZopoy/IwxhWhqyeLTzkkMlLI7xcIizBbdtO6VL27/dneYc/QMQdM5aqPHp37GLJcB4FQZ/inr19b+LPaJE6fibD7KSnaLcyx3ds3wK0PbqJ1lwmpovnKDp9NVBYi8nbfyyZgLo7De0ziXwpObqmwvX8g9UNfzc1XaXJ8EHH7en5rH+csX80F1/1pRI/mRsMrp53X6RvWdSwYHRYc9LzSBz+8c8Mwp19U3kycIgkjrhRMkEZUFG0tlaHQ4DDxWneZMHTsYecsa2mL+R3toYrpng2buahzbWIkW9w+0lyHqDIlTj2r4lf/7zlyRqiShfzPQVrFKFB6+GyalcUpvr/7cTK4TytFmjoTdGL29PZRaRJEqjONROFPBMoSAdHoigKqUxQA2/oGaAqYidoCETFhg17WWl1ZWuvWwzGalV3dgIrguffnJKSpPhDmQE66nmGrslsf3BRaFsevQLI6aNNch0qzRA6eZdn2vfpM/gTYanN2POUbtw0BzjxqRum1odJEQ32oVAkaiDBbaN6KlGnwtjxWopqCVHNUYQrRK7cOzqCSNOuKqrabl7CViD8KJ8vKoyiCvUGCA1XYZ0nVB8Kqm4ZVJw7KEbY6iKtn5CfLtfK+c941a6KfHYWV658LPQdl2fb9FQ+ChUbz4Fe+UcqxmrDZrMSVKP8PYo5VVf9fKRLVic6u7oaKvR4P+JOsmlKGAvsHlTSDc9EDQ9JKpJYx/HHVYOPkTKo+kNWB/HhM468sCjTLtfKOLWoQ7RvUyNVLVqWedtAPVjxI+k1c+HtQCfxk5QbuCPgovWtWC0UB8dFQK4FVwCTgMOAh998cYJfyRasd3kzKKJ5g9Iz3ur2thWXvmM3qxSfy2NKTMvVJ9h70qArBfmpRM8dPGpmKoJoSD2HFMv2RR1kG7aSIp7DzERVRlfVaeccRRdTqJUqmYw6YGhk5552ruIZKzSKplWxc1NkVC+YMy/Pp7OoODWapRbisn7h+Fj8AEJF/BF6vqv3u628Bv6uNeLUhb7hra6WJvkFtuHr+jUKwH3ScIzPLbM/rRRBMYgqz0xfl9Asz63j7juo94L2/ZXt/5jIOcTSLcPrh7UMF+bKYoDziVh5pr0Wa8+s5ub0KBs0iHLX/FO7ZsDlzXbWo7WfxCzzZ0xsb2NDZ1c2SFfcNXa9JlSbm7jt1qLlQZ1d3ZNRjWnOyf9UQ5RD397To2bojU2OlskhT7mMdcLSqPue+ngLcqarlV65KSbXlPrL0og5SaRJ2m+SEtnnLyvFU9qG5SUKd/1ltqVl7n4eZYMpqbB8mW1joblQ5mDx93aPwFLC/qB84CvR1B0zlvidfGKGYspapSSNv2usbVTYn7SQir7xRz2BSj+s0ZX68CsF5C2vGnbsinoMsFFruA1gKdInIrTjX4A3AktzSNSDVOCb7BpXWXSbQ9dkTh94rohZUHprccNW09n+ovp7RZe+cXcgAHZztTW6pIBJd6TYqy7oM+21o4EPIzDLKSZvkkE1zDfz5AlEd04I27SS5ovC+FzYg+gdOry9K3LWPKptz64ObqhrkwuQNBh4EFWqa1UtcH+1gWHFe/1Rc9FdWK8eW7f3st+jGQidHUaSJhvqeiPwSONJ963xV/UtpEtWBakMiu3t6hy6Yd5PWWlFUmp1eDPM72tkvY5HBMAdlmgehva2l0AE6bFtRctTSF5FlqR/13SiHrDfLjptcBGePecqZR8kVtRrz/kV9HhYxFTYA5m0SlpWweyepdH0Wufz3YBEhsb19A5yzfPWQDyWPv6iJncl8tagPFRcNdbCqPigih7lvPeH+P11EpqvqPaVIVAe8k1tNTwEvRT9LDf1qEXdKmjY0MoyoQTdJgdaqhn7eviFJZDFZZTmfk2PqDcXZyuM65AWPNU8RybDrHDbgn7N8NRf//L6hfJaoyUCaGbi337Bzl6WRVl7yTGSi5PX8ZJCvfH0UwUE+U+Jn4HXRoeJB4lYWnwLOBi4L+UyBYtaQDUJWR1kUNV1PKDwWsirIslKKGnSDA1tbawVVEnsJFE3WLOs0pJ0Ve3jd9tJc26Q+3lEDWNQgsesuzSMc2VkVRZRyjTJ5BPNZwki7Ylg476BQh/CL2/qHBSmUQR4fVtS1VnbmURSdlOkf5LPca2GU6fCOi4Y62/3/+NL23mCMhgxdP1GrgmCUUBRe+YcoyvIBZCUoRxpbufe9sMEizaw4+Nu0D29PTqdn2L1XaRZ29A8OXUNPqbW1VFJHV8X1gYgbWJJmqWnbCs/vaB8WXeTRN6ilzoKzTgj88kZZGJIG4izXJW7b1Uw4yzTPpilR/k4R2d39+yIRuU5EOqrZqYgsE5EHReRPInK9iLT5gByAwwAAEpFJREFUPrtARB4WkXUiMq+a/WQlKW671nilRqIIFs3zM7+jnTsWncAVC+Y4kTsh2/bKP4wmwspTeyWe034vaVYc9tu0paLzPqz+3AdwS8AM6IgKAr19A4iMLAHu5Qn4cyeuWDCHrs+eGDvgxxE3OEZV8w1bwWyOGEC7e3pHXLeiiJsQJBGVOzK9rSXynLW3tbDrxEx1WUdsu9p8r7JNw2mO7jOq+hMReT3wJmAZ8C12Orzz8CvgAlXtF5EvARcA54vIq4F3AYcA04Ffi8iBqlqzqX5R5qhqaRZh2TtnA0QuS9OUzA6LbGlrqbDk1PytGOtJ1CAQdBbGDRZJs+KoaKM8tZGykJSV7NGztY/LF8zJbGIJq/cUjBjy4y//kSafJE8OTVlO2bgJQZJ5KslPFhamW8144W27mvYGtSj7kUZZeNKfBFypqjeKyCXV7FRVb/a9vBN4h/v3acCPVXU78JiIPAwcAfyhmv1lobOrm601qnMfhT+yCaId72ntk41iTiqCuGP2mxriBovLF8yJHQzi6hn5azEF60IV8bCmGTCm54hCCzPL/HRVN6cf3s4Na54Kzc1YOO+gRHNOGhnizLth5q4i8mXiHOtJ5qk0irDa+k9eVWb/IH9ujgCbvK2e85BGWXSLyH8Cbwa+JCITSWG+ysCHgeXu3+04ysNjo/veCETkbBwHPDNmhHehyspFnWtrGs0Uxa6+UtHg3FD1Dh9tFJKiRdKsHpIGg6jfVpsAlYakCUDe1UtcvsPqxSdGDtDHLL0lVdSTn7BtXfr2WakmPXl8DWH7i1odqI50UIcdT5wirCbPIm4FEBeJpSGva1lEENIN+mcANwHzVLUHmAosTPqRiPxaRO4N+Xea7zufxil7flVWwVX1SlWdq6pzp02Ltt2npbOru1RFUWmW1G0cgzbeLPbhsU6a2ktP9vQmnjPPp/PY0pNG9Nuu5/mOmwDEdY9LIslPE3U+suZJRPmKPPnD8B9zVl9D3P7C6l9F+U/yRBFl/Y032Yi6fnFdHf3HcfmCOTwect+WTZqkvK0i8jTwepxCgv3u/0m/e1Pc5yLyQeBk4I26s+ZIN+BvLr23+17ppOkCFkalSThivymR2bMeAwPKybP3Siz1DOERJZ6MRZeyGG2kifRKs3pIu49an++oGXG1poa00UvV/i5usF847yAW/mTNMKd9pWl434msyiluf2GDadR9k2eVHnVuwhp3+c16UfdVoz/naTrlLcbpjncQ8D2gAvwQOCbvTkXkLcC/An+nqlt9H60AfiQi/47j4H4l8Me8+8lCnpmF31E8MyFrehC4Yc1TieF1UTPYseR3qBZ/dnGc76Gac1av813WgJE3uTHr7xIH+6gyxC5ZlVNW5VJkkmfUthaf4kQZhhUITOMvadTnPI3P4m1AB3APgKo+6YXSVsHXgInAr8SJDb1TVT+mqveJyDXA/TgrmI/XKhKq2laNUX4FPz29fVyxYM6I2ZXHaI5SqgeNPhPLSxkDRt5zlfV3cYP9spvWjUjO6xsYnm+RdTDPqlyCx+PVIDs3pPSGR1xJlLhzE9xOHv9PFGUVzYwjTdXZP6rqESJyj6oeJiK7An9Q1UNLlSwD1VadhXTVHuPMAWmrRT6+9KQRZZDjEqcMYzQRV7U1KgRcGF6JIMtAmKZKbB5Z/cmZRZkFo6pbB48/iSJlKrrq7DVuNFSbiHwUJ3rp25kkGgWEzRKyhEamsaVPcevhNPJS0zCqIW62ndZfkOX5qGZ1mSaTP20NLD9Ryi6v3yiP3GWQxsH9ZRF5M/A3HL/FZ1X1V6VJVEeqHcT9tvRgPZxKswzZMg1jLBP1HJVVFDLvc5vG35E3GizML3H8wdNGRFzmOf5aVfINkio/3VUOvwIQkSYROVNVM4e7jhfGqi3daFzqYcPOSqM9F2lm+kVFg1388/vY1jc4Il/i9MOLq4xbdt5VXInyPYCP4yTFrcBRFh8H/gVYQ47ciPGEmZqMWpG3aF4tCSqzyxfMqbtsaVY6RUWDhTXxUtKV7MkjdxnErSz+B3gep9TGR4ALcZThfFXN3/jBMIxCqZcNOy2NqszSrHSKigaLIo/pqF4rtMhoKBFZq6qz3L+bgaeAGaq6rVSJclBENJRhjFaqjbIp24QVVRajFuVTak1UpNLECU2h+VX1PgdFRUMNHZmqDojIxkZUFIYx3qnGhl2LWX+9HLL1IGrWD+EtdUdTyZ44ZTFbRP7m/i1Ai/taAFXVPUqXzjCMRKqxYdfChFUvh2y9iPNXNopzPw9xnfLiq7UZhtEQVGPDrsWsP0mZjYZIriIY7UEv+Vs7GYbRMOQdiKqd9acd6CdVmoaUhb+sTa2d3+NFMZVBkX0pDMMYZVRTjj1Ni1vvO/7Q0e39g0N/V9P+NCtpW/Ia4ZiyMIxxzPyO9tC+D9WWy0j7nVo6v2upmMYiZoYyjHFOPctl1NL5PZ6issrAVhaGYeQiakAPlsuI+04tuxKmkdeIxpSFYRi5SDPQp2lvm9cMVoa8RjSJ/SxGA5bBbRj1IU10UdoIpFpEKlk01HCyZHCbsjAMo+4kNfSxQb4cim5+ZBiGUSpJkUqNWIhwvGHKwjCMuhMXqdSIVXXH40rHlIVhGHUnLoQ2TchrLQfvRi25XjYWDWUYRt2Ji1RKCnmtdWb2eE3us5WFYRi5KHI2n1QMMa4QYa3NVGmT+8aaqcqUhWEYmSnDFBOVSZ6kSGqdmZ0m63wsmqrqaoYSkfNEREVkT/e1iMhXReRhEfmTiBxWT/kMwwin1qaY+R3t3LHoBB5behJ3LDph2IBb68zsNMl9Y9FUVTdlISL7ACcCG3xv/z3wSvff2cA36yCaYRgJNFKdpVpnZqfJOm+k81MU9TRDXQ78K/Az33unAf+tTqbgnSLSJiJ7qepTdZHQMIxQGqn7XTXNn6rZZ9z2G+n8FEVdlIWInAZ0q+oaEfF/1A484Xu90X1vhLIQkbNxVh/MmDGjPGENwxhBNa1cy6DRutA12vkpgtKUhYj8Gnh5yEefBi7EMUHlRlWvBK4Ep9xHNdsyDCMb9ZjNjybG4vkpTVmo6pvC3heRWcB+gLeq2Bu4R0SOALqBfXxf39t9zzCMBqPRZvONxlg7PzV3cKvqWlV9qarOVNWZOKamw1T1L8AK4P1uVNRRwGbzVxiGYdSfRsuz+AXwVuBhYCvwofqKYxiGYUADKAt3deH9rcDH6yeNYRiGEYbVhjIMwzASMWVhGIZhJDImOuWJyCZgfc6f7wk8U6A4o4nxfOxgxz+ej388HzvsPP59VXVamh+MCWVRDSKyMm1bwbHGeD52sOMfz8c/no8d8h2/maEMwzCMRExZGIZhGImYsnBLhoxTxvOxgx3/eD7+8XzskOP4x73PwjAMw0jGVhaGYRhGIqYsDMMwjETGrbIQkbeIyDq3heuiestTa0TkcRFZKyKrRWRlveUpGxH5rog8LSL3+t6bKiK/EpGH3P+n1FPGsog49iUi0u1e/9Ui8tZ6ylgmIrKPiNwqIveLyH0i8kn3/TF//WOOPfP1H5c+CxFpBv4MvBmn6u3dwLtV9f66ClZDRORxYK6qjovEJBF5A/AiTifG17jv/RvwnKoudScMU1T1/HrKWQYRx74EeFFVv1xP2WqBiOwF7KWq94jI7sAqYD7wQcb49Y859jPIeP3H68riCOBhVX1UVXcAP8Zp6WqMUVT1t8BzgbdPA37g/v0DnIdozBFx7OMGVX1KVe9x/34BeACnA+eYv/4xx56Z8aosotq3jicUuFlEVrktascjL/P1S/kL8LJ6ClMHPiEif3LNVGPOBBOGiMwEOoC7GGfXP3DskPH6j1dlYcDrVfUw4O+Bj7uminGLWx5/PNlkvwkcAMzB6XF/WX3FKR8R2Q34KXCOqv7N/9lYv/4hx575+o9XZTHu27eqarf7/9PA9TimufHGX12brmfbfbrO8tQMVf2rqg6o6iDwX4zx6y8iFZzB8ipVvc59e1xc/7Bjz3P9x6uyuBt4pYjsJyK7AO/Caek6LhCRXV1nFyKyK3AicG/8r8YkK4APuH9/APhZHWWpKd4g6fI2xvD1FxEBvgM8oKr/7vtozF//qGPPc/3HZTQUgBsqdgXQDHxXVb9QZ5Fqhojsj7OaAKdb4o/G+vGLyNXAcTilmf8KLAY6gWuAGTgl7s9Q1THnCI449uNwTBAKPA78w1jtdy8irwd+B6wFBt23L8Sx3Y/p6x9z7O8m4/Uft8rCMAzDSM94NUMZhmEYGTBlYRiGYSRiysIwDMNIxJSFYRiGkYgpC8MwDCMRUxbGECLyEl8Vyr8EqlL+X41ludotRXBu4H1/tcyHROQ6EXl1gfu9wstmF5HbRGSu77OZ/sqtZSMiHxSRTe6x3i8iH82xjTYR+aeYz1VEfuh7PcHd5w3u61OjqjKLyIsp9v/tpOsjIt8XkXeEvD9TRN7jez1LRL6ftE+jHExZGEOo6rOqOkdV5wDfAi73Xqvq62olh4i8HHitqh6qqpeHfMWT65XAcuAWEZlWwH5fAhzlFt4rFbfycRqWu9fjOOCLIpK1flEbEKksgC3Aa0SkxX39ZnzVDFR1haouzbjPIVT1I1VUc54JDCkLVV0L7C0iM/LKY+THlIWRCm8WKSLHicjtIvIzEXlURJaKyJki8kdx+mMc4H5vmoj8VETudv8dE7LNSSLyPfd3XSJyvPvRzUC7O6M+Nk4uVV3ufv897jY/6+7vXhG5UhwOEJF7fPt9pf+1j9OB/015PkJld1cDX/N97wYROc47hyJymYisAY52z9397goqtlS0W5blEWBfETlCRP7g7vf/ROQgd/uHuNdhtbvNVwJLgQPc95ZFbP4XwEnu3+8GrvbJP3Q84lQ8+IN7zJf4vnOcuwq7VkQeFJGrRETcz4ZWZyJyloj82ZXxv/znCXiDeyyP+lYZS4FjXdm9FebPcSouGDXGlIWRh9nAx4BXAe8DDlTVI4BvA//sfucrOCuA1+IMwt8O2c7HcWq4zcIZpH4gIpOAU4FH3NXD71LIcw9wsPv311T1tW7fhhbgZFV9BNgsInPc73wI+F7Ido7Bqffv5yp3sFqNM6gmyR7HrsBdqjobp1T024BDVPVQ4JK4H4qTdb8/8DDwIHCsqnYAnwW+6H7tY8BX3JXIXJxqyovYeS4XRmz+x8C7XPkPZWdV0iBfAb7pHnMw27cDOAd4tSvnsMmBiEwHPgMc5X52cOD3ewGvB07GURK4sv/Old1bYa4EYicQRjmYsjDycLdbJ387zmz3Zvf9tTimA4A3AV9zB9kVwB7iVL7083rghwCq+iBOyYUDc8gjvr+PF5G7RGQtcAJwiPv+t4EPueafBcCPQrazF7Ap8N6ZPtOcv5tYHtkHcAq6AWwGtgHfEZG3A1sjfrPAPYdX45RkeA6YDPxEHP/J5b5j/ANwoYicD+yrqr0J8uDK/yec6/ZuhivEIMewc9XxP4HP/qiqG93CdKvZeR94HAHcrqrPqWof8JPA552qOuiarOJMbU8D02M+N0rClIWRh+2+vwd9rwdxak2Bc28d5fN5tKtqokM0Jx3AA+7M+BvAO9zZ738B3mz/pzjl2E8GVqnqsyHb6fV9Py/9DH+u/NvbpqoDAKrajzOAXuvKFGX+Wu6evyNV1avn9XngVnf1dIq3D1X9Ec6qrBf4hYickEHuFcCX8ZmgIoiqD+S/JwbYeR+kxf97ifyWc6yplKBRLKYsjLK4mZ0mKXwmID+/A850Pz8Qp6Dbuiw7EZHTcarmXs3OgfkZdxUzFGGjqtuAm3Dq+IeZoMAxDb0i5a6jZH8cmCMiTSKyDxGln135JqvqL4BzcUx7aZnMTif0B33b3B94VFW/ilNB9VDgBWD3FNv8LnCx60SO4g52+gvOzCAvOJWe/05EpojIBBzTZBJhsh/IGK6Q28iYsjDK4v8Bc11H6/049vQg3wCaXJPRcuCDrmkriXNdP8JDwHuBE1R1k6r24Kwm7sVRDHcHfncVzurnZsK5ESfqKA1Rst8BPAbcD3wVx58Sxu7ADSLyJ+D3wKdS7hfg34BLRaSL4TP4M4B7XbPVa3B6bj8L3OE6/KMc3LgmpK8m7PeTOI2y1pKxs6TbP+WLwB9xztHjOKa4OP4EDIjIGp+D+3ic62TUGKs6a4wbRORfcGbzn4n5zu9xnOI9tZNsfCAiu6nqi+7K4nqc1gDXJ/3O9/uJwO04XR77y5LTCMeUhTEuEJHrcdpInqCqz8R870ig13X6GgXihge/CcdceDPwSc0wALmhwO2qels5EhpxmLIwDMMwEjGfhWEYhpGIKQvDMAwjEVMWhmEYRiKmLAzDMIxETFkYhmEYifx/Q9FMpKJmtb8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# ANSWER 3 (class discussion about the residuals)\n",
    "x_matrix = transformer_3.fit_transform(X_train)\n",
    "\n",
    "prediction = fitted_cab_model3.predict(x_matrix)\n",
    "residual = y_train - prediction\n",
    "plt.scatter(X_train, residual, label=\"Residual\")\n",
    "plt.axhline(0, color='k')\n",
    "\n",
    "plt.title(\"Residuals for the Cubic Model\")\n",
    "plt.ylabel(\"Residual Number of Taxi Pickups\")\n",
    "plt.xlabel(\"Time of Day (Hours Past Midnight)\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Other features\n",
    "Polynomial features are not the only constucted features that help fit the data. Because these data have a 24 hour cycle, we may want to build features that follow such a cycle. For example, $sin(24\\frac{x}{2\\pi})$, $sin(12\\frac{x}{2\\pi})$, $sin(8\\frac{x}{2\\pi})$. Other feature transformations are appropriate to other types of data. For instance certain feature transformations have been developed for geographical data.\n",
    "\n",
    "### Scaling Features\n",
    "When using polynomials, we are explicitly trying to use the higher-order values for a given feature. However, sometimes these polynomial features can take on values that are drastically large, making it difficult for the system to learn an appropriate bias weight due to its large values and potentially large variance. To counter this, sometimes one may be interested in scaling the values for a given feature.\n",
    "\n",
    "For our ongoing taxi-pickup example, using polynomial features improved our model. If we wished to scale the features, we could use `sklearn`'s StandardScaler() function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.33412512570778274"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# SCALES THE EXPANDED/POLY TRANSFORMED DATA\n",
    "# we don't need to convert to a pandas dataframe, but it can be useful for scaling select columns\n",
    "train_copy = pd.DataFrame(expanded_train.copy())\n",
    "test_copy = pd.DataFrame(expanded_test.copy())\n",
    "\n",
    "# Fit the scaler on the training data\n",
    "scaler = StandardScaler().fit(train_copy)\n",
    "\n",
    "# Scale both the test and training data. \n",
    "train_scaled = scaler.transform(expanded_train)\n",
    "test_scaled = scaler.transform(expanded_test)\n",
    "\n",
    "# we could optionally run a new regression model on this scaled data\n",
    "fitted_scaled_cab = LinearRegression().fit(train_scaled, y_train)\n",
    "fitted_scaled_cab.score(test_scaled, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:3px\">\n",
    "\n",
    "## Part 3: Multiple regression and exploring the Football (aka soccer) data\n",
    "Let's move on to a different dataset! The data imported below were scraped by [Shubham Maurya](https://www.kaggle.com/mauryashubham/linear-regression-to-predict-market-value/data) and record various facts about players in the English Premier League. Our goal will be to fit models that predict the players' market value (what the player could earn when hired by a new team), as estimated by https://www.transfermarkt.us.\n",
    "\n",
    "`name`: Name of the player  \n",
    "`club`: Club of the player  \n",
    "`age` : Age of the player  \n",
    "`position` : The usual position on the pitch  \n",
    "`position_cat` :  1 for attackers, 2 for midfielders, 3 for defenders, 4 for goalkeepers  \n",
    "`market_value` : As on transfermrkt.com on July 20th, 2017  \n",
    "`page_views` : Average daily Wikipedia page views from September 1, 2016 to May 1, 2017  \n",
    "`fpl_value` : Value in Fantasy Premier League as on July 20th, 2017  \n",
    "`fpl_sel` : % of FPL players who have selected that player in their team  \n",
    "`fpl_points` : FPL points accumulated over the previous season  \n",
    "`region`: 1 for England, 2 for EU, 3 for Americas, 4 for Rest of World  \n",
    "`nationality`: Player's nationality  \n",
    "`new_foreign`: Whether a new signing from a different league, for 2017/18 (till 20th July)  \n",
    "`age_cat`: a categorical version of the Age feature  \n",
    "`club_id`: a numerical version of the Club feature  \n",
    "`big_club`: Whether one of the Top 6 clubs  \n",
    "`new_signing`: Whether a new signing for 2017/18 (till 20th July)  \n",
    "\n",
    "As always, we first import, verify, split, and explore the data.\n",
    "\n",
    "## Part 3.1: Import and verification and grouping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "name             object\n",
      "club             object\n",
      "age               int64\n",
      "position         object\n",
      "position_cat      int64\n",
      "market_value    float64\n",
      "page_views        int64\n",
      "fpl_value       float64\n",
      "fpl_sel          object\n",
      "fpl_points        int64\n",
      "region          float64\n",
      "nationality      object\n",
      "new_foreign       int64\n",
      "age_cat           int64\n",
      "club_id           int64\n",
      "big_club          int64\n",
      "new_signing       int64\n",
      "dtype: object\n"
     ]
    }
   ],
   "source": [
    "league_df = pd.read_csv(\"../data/league_data.txt\")\n",
    "print(league_df.dtypes)\n",
    "\n",
    "# QUESTION: what would you guess is the mean age? mean salary?\n",
    "#league_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(461, 17)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "league_df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#league_df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (Stratified) train/test split\n",
    "We want to make sure that the training and test data have appropriate representation of each region; it would be bad for the training data to entirely miss a region. This is especially important because some regions are rather rare.\n",
    "\n",
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "**Questions**:\n",
    "1. Use the `train_test_split()` function, while (a) ensuring the test size is 20% of the data, and; (2) using 'stratify' argument to split the data (look up documentation online), keeping equal representation of each region. This doesn't work by default, correct? What is the issue?\n",
    "2. Deal with the issue you encountered above. Hint: you may find numpy's `.isnan()` and panda's `.dropna()` functions useful!\n",
    "3. How did you deal with the error generated by `train_test_split`? How did you justify your action? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*your answer here*:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# EXERCISE: feel free to work with a partner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1000, 2), (250, 2))"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "train_data.shape, test_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we won't be peeking at the test set, let's explore and look for patterns! We'll introduce a number of useful pandas and numpy functions along the way. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Groupby\n",
    "Pandas' `.groupby()` function is a wonderful tool for data analysis. It allows us to analyze each of several subgroups.\n",
    "\n",
    "Many times, `.groupby()` is combined with `.agg()` to get a summary statistic for each subgroup. For instance: What is the average market value, median page views, and maximum fpl for each player position?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'position'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-23-3c158d4eb7c3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m train_data.groupby('position').agg({\n\u001b[0m\u001b[1;32m      2\u001b[0m     \u001b[0;34m'market_value'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0;34m'page_views'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmedian\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;34m'fpl_points'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m })\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36mgroupby\u001b[0;34m(self, by, axis, level, as_index, sort, group_keys, squeeze, observed, **kwargs)\u001b[0m\n\u001b[1;32m   7894\u001b[0m             \u001b[0msqueeze\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msqueeze\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   7895\u001b[0m             \u001b[0mobserved\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mobserved\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 7896\u001b[0;31m             \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   7897\u001b[0m         )\n\u001b[1;32m   7898\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/groupby/groupby.py\u001b[0m in \u001b[0;36mgroupby\u001b[0;34m(obj, by, **kwds)\u001b[0m\n\u001b[1;32m   2476\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"invalid type: {}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2477\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2478\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mklass\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mby\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/groupby/groupby.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, obj, keys, axis, level, grouper, exclusions, selection, as_index, sort, group_keys, squeeze, observed, **kwargs)\u001b[0m\n\u001b[1;32m    389\u001b[0m                 \u001b[0msort\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msort\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    390\u001b[0m                 \u001b[0mobserved\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mobserved\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 391\u001b[0;31m                 \u001b[0mmutated\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmutated\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    392\u001b[0m             )\n\u001b[1;32m    393\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/groupby/grouper.py\u001b[0m in \u001b[0;36m_get_grouper\u001b[0;34m(obj, key, axis, level, sort, observed, mutated, validate)\u001b[0m\n\u001b[1;32m    619\u001b[0m                 \u001b[0min_axis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    620\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 621\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    622\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mGrouper\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mgpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkey\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    623\u001b[0m             \u001b[0;31m# Add key to exclusions\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'position'"
     ]
    }
   ],
   "source": [
    "train_data.groupby('position').agg({\n",
    "    'market_value': np.mean,\n",
    "    'page_views': np.median,\n",
    "    'fpl_points': np.max\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'DataFrame' object has no attribute 'position'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-3107a1056032>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mtrain_data\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mposition\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munique\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36m__getattr__\u001b[0;34m(self, name)\u001b[0m\n\u001b[1;32m   5178\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_info_axis\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_can_hold_identifiers_and_holds_name\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   5179\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 5180\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__getattribute__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   5181\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   5182\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__setattr__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'DataFrame' object has no attribute 'position'"
     ]
    }
   ],
   "source": [
    "train_data.position.unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'big_club'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-26-8ba2a54f5c1f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m train_data.groupby(['big_club', 'position']).agg({\n\u001b[0m\u001b[1;32m      2\u001b[0m     \u001b[0;34m'market_value'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m     \u001b[0;34m'page_views'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0;34m'fpl_points'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m })\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/generic.py\u001b[0m in \u001b[0;36mgroupby\u001b[0;34m(self, by, axis, level, as_index, sort, group_keys, squeeze, observed, **kwargs)\u001b[0m\n\u001b[1;32m   7894\u001b[0m             \u001b[0msqueeze\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msqueeze\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   7895\u001b[0m             \u001b[0mobserved\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mobserved\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 7896\u001b[0;31m             \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   7897\u001b[0m         )\n\u001b[1;32m   7898\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/groupby/groupby.py\u001b[0m in \u001b[0;36mgroupby\u001b[0;34m(obj, by, **kwds)\u001b[0m\n\u001b[1;32m   2476\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"invalid type: {}\"\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2477\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2478\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mklass\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mby\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/groupby/groupby.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, obj, keys, axis, level, grouper, exclusions, selection, as_index, sort, group_keys, squeeze, observed, **kwargs)\u001b[0m\n\u001b[1;32m    389\u001b[0m                 \u001b[0msort\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msort\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    390\u001b[0m                 \u001b[0mobserved\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mobserved\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 391\u001b[0;31m                 \u001b[0mmutated\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmutated\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    392\u001b[0m             )\n\u001b[1;32m    393\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/groupby/grouper.py\u001b[0m in \u001b[0;36m_get_grouper\u001b[0;34m(obj, key, axis, level, sort, observed, mutated, validate)\u001b[0m\n\u001b[1;32m    619\u001b[0m                 \u001b[0min_axis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mgpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    620\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 621\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    622\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mGrouper\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mgpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkey\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    623\u001b[0m             \u001b[0;31m# Add key to exclusions\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'big_club'"
     ]
    }
   ],
   "source": [
    "train_data.groupby(['big_club', 'position']).agg({\n",
    "    'market_value': np.mean,\n",
    "    'page_views': np.mean,\n",
    "    'fpl_points': np.mean\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "**Question**:\n",
    "1. Notice that the `.groupby()` function above takes a list of two column names. Does the order matter? What happens if we switch the two so that 'position' is listed before 'big_club'?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# EXERCISE: feel free to work with a partner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style=\"height:3px\">\n",
    "\n",
    "## Part 3.2: Linear regression on the football data\n",
    "This section of the lab focuses on fitting a model to the football (soccer) data and interpreting the model results. The model we'll use is\n",
    "\n",
    "$$\\text{market_value} \\approx \\beta_0 + \\beta_1\\text{fpl_points} + \\beta_2\\text{age} + \\beta_3\\text{age}^2 + \\beta_4log_2\\left(\\text{page_views}\\right) + \\beta_5\\text{new_signing} +\\beta_6\\text{big_club} + \\beta_7\\text{position_cat}$$\n",
    "\n",
    "We're including a 2nd degree polynomial in age because we expect pay to increase as a player gains experience, but then decrease as they continue aging. We're taking the log of page views because they have such a large, skewed range and the transformed variable will have fewer outliers that could bias the line. We choose the base of the log to be 2 just to make interpretation cleaner.\n",
    "\n",
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "**Questions**:\n",
    "1. Build the data and fit this model to it. How good is the overall model?\n",
    "2. Interpret the regression model. What is the meaning of the coefficient for:\n",
    "    - age and age$^2$\n",
    "    - $log_2($page_views$)$\n",
    "    - big_club\n",
    "3. What should a player do in order to improve their market value? How many page views should a player go get to increase their market value by 10?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "'market_value'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/indexes/base.py\u001b[0m in \u001b[0;36mget_loc\u001b[0;34m(self, key, method, tolerance)\u001b[0m\n\u001b[1;32m   2889\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2890\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2891\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'market_value'",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-29-a30e1a0cb98c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Q1: we'll do most of it for you ...\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0my_train\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtrain_data\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'market_value'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0my_test\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtest_data\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'market_value'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mbuild_football_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mx_matrix\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'fpl_points'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'age'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'new_signing'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'big_club'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'position_cat'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/frame.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m   2973\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnlevels\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2974\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_getitem_multilevel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2975\u001b[0;31m             \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2976\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mis_integer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2977\u001b[0m                 \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mindexer\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/local/lib/python3.7/site-packages/pandas/core/indexes/base.py\u001b[0m in \u001b[0;36mget_loc\u001b[0;34m(self, key, method, tolerance)\u001b[0m\n\u001b[1;32m   2890\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2891\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2892\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_engine\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_loc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_maybe_cast_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2893\u001b[0m         \u001b[0mindexer\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_indexer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtolerance\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtolerance\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2894\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m1\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mindexer\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msize\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/index.pyx\u001b[0m in \u001b[0;36mpandas._libs.index.IndexEngine.get_loc\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32mpandas/_libs/hashtable_class_helper.pxi\u001b[0m in \u001b[0;36mpandas._libs.hashtable.PyObjectHashTable.get_item\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: 'market_value'"
     ]
    }
   ],
   "source": [
    "# Q1: we'll do most of it for you ...\n",
    "y_train = train_data['market_value']\n",
    "y_test = test_data['market_value']\n",
    "def build_football_data(df):\n",
    "    x_matrix = df[['fpl_points','age','new_signing','big_club','position_cat']].copy()\n",
    "    x_matrix['log_views'] = np.log2(df['page_views'])\n",
    "    \n",
    "    # WRITE CODE FOR CREATING THE AGE SQUARED COLUMN\n",
    "    ####\n",
    "    \n",
    "    # OPTIONALLY WRITE CODE to adjust the ordering of the columns, just so that it corresponds with the equation above\n",
    "    ####\n",
    "    \n",
    "    # add a constant\n",
    "    x_matrix = sm.add_constant(x_matrix)\n",
    "    \n",
    "    return x_matrix\n",
    "\n",
    "# use build_football_data() to transform both the train_data and test_data\n",
    "train_transformed = build_football_data(train_data)\n",
    "test_transformed = build_football_data(test_data)\n",
    "\n",
    "fitted_model_1 = OLS(endog= y_train, exog=train_transformed, hasconst=True).fit()\n",
    "fitted_model_1.summary()\n",
    "\n",
    "# WRITE CODE TO RUN r2_score(), then answer the above question about the overall goodness of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'fitted_model_1' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-31-42967488aac6>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Q2: let's use the age coefficients to show the effect of age has on one's market value;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;31m# we can get the age and age^2 coefficients via:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0magecoef\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfitted_model_1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mage\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m \u001b[0mage2coef\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfitted_model_1\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mage_squared\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'fitted_model_1' is not defined"
     ]
    }
   ],
   "source": [
    "# Q2: let's use the age coefficients to show the effect of age has on one's market value;\n",
    "# we can get the age and age^2 coefficients via:\n",
    "agecoef = fitted_model_1.params.age\n",
    "age2coef = fitted_model_1.params.age_squared\n",
    "\n",
    "# let's set our x-axis (corresponding to age) to be a wide range from -100 to 100, \n",
    "# just to see a grand picture of the function\n",
    "x_vals = np.linspace(-100,100,1000)\n",
    "y_vals = agecoef*x_vals +age2coef*x_vals**2\n",
    "\n",
    "# WRITE CODE TO PLOT x_vals vs y_vals\n",
    "plt.plot(x_vals, y_vals)\n",
    "plt.title(\"Effect of Age\")\n",
    "plt.xlabel(\"Age\")\n",
    "plt.ylabel(\"Contribution to Predicted Market Value\")\n",
    "plt.show()\n",
    "\n",
    "# Q2A: WHAT HAPPENS IF WE USED ONLY AGE (not AGE^2) in our model (what's the r2?); make the same plot of age vs market value\n",
    "# Q2B: WHAT HAPPENS IF WE USED ONLY AGE^2 (not age) in our model (what's the r2?); make the same plot of age^2 vs market value\n",
    "# Q2C: PLOT page views vs market value\n",
    "\n",
    "# EXERCISE: WRITE CODE HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr style='height:3px'>\n",
    "\n",
    "### Part 3.3: Turning Categorical Variables into multiple binary variables\n",
    "Of course, we have an error in how we've included player position. Even though the variable is numeric (1,2,3,4) and the model runs without issue, the value we're getting back is garbage. The interpretation, such as it is, is that there is an equal effect of moving from position category 1 to 2, from 2 to 3, and from 3 to 4, and that this effect is probably between -0.5 to -1 (depending on your run).\n",
    "\n",
    "In reality, we don't expect moving from one position category to another to be equivalent, nor for a move from category 1 to category 3 to be twice as important as a move from category 1 to category 2. We need to introduce better features to model this variable.\n",
    "\n",
    "We'll use `pd.get_dummies` to do the work for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_design_recoded = pd.get_dummies(train_transformed, columns=['position_cat'], drop_first=True)\n",
    "test_design_recoded = pd.get_dummies(test_transformed, columns=['position_cat'], drop_first=True)\n",
    "\n",
    "train_design_recoded.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've removed the original `position_cat` column and created three new ones.\n",
    "\n",
    "#### Why only three new columns?\n",
    "Why does pandas give us the option to drop the first category? \n",
    "\n",
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "**Questions**:\n",
    "1. If we're fitting a model without a constant, should we have three dummy columns or four dummy columns?\n",
    "2. Fit a model on the new, recoded data, then interpret the coefficient of `position_cat_2`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# FEEL FREE ETO WORK WITH A PARTNER"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Answers**:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 4: A nice trick for forward-backwards\n",
    "\n",
    "XOR (operator ^) is a logical operation that only returns true when input differ. We can use it to implement forward-or-backwards selection when we want to keep track of whet predictors are \"left\" from a given list of predictors.\n",
    "\n",
    "The set analog is \"symmetric difference\". From the python docs:\n",
    "\n",
    "`s.symmetric_difference(t)\ts ^ t\tnew set with elements in either s or t but not both`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "set() ^ set([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "set([1]) ^ set([1,2,3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "set([1, 2]) ^ set([1,2,3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "Outline a step-forwards algorithm which uses this idea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BONUS EXERCISE:\n",
    "We have provided a spreadsheet of Boston housing prices (data/boston_housing.csv). The 14 columns are as follows:\n",
    "1. CRIM: per capita crime rate by town\n",
    "2. ZN: proportion of residential land zoned for lots over 25,000 sq.ft.\n",
    "3. INDUS: proportion of non-retail business acres per town\n",
    "4. CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)\n",
    "5. NOX: nitric oxides concentration (parts per 10 million) 1https://archive.ics.uci.edu/ml/datasets/Housing 123 20.2. Load the Dataset 124\n",
    "6. RM: average number of rooms per dwelling\n",
    "7. AGE: proportion of owner-occupied units built prior to 1940\n",
    "8. DIS: weighted distances to five Boston employment centers\n",
    "9. RAD: index of accessibility to radial highways\n",
    "10. TAX: full-value property-tax rate per \\$10,000\n",
    "11. PTRATIO: pupil-teacher ratio by town\n",
    "12. B: 1000(Bk−0.63)2 where Bk is the proportion of blacks by town\n",
    "13. LSTAT: % lower status of the population\n",
    "14. MEDV: Median value of owner-occupied homes in $1000s We can see that the input attributes have a mixture of units\n",
    "\n",
    "There are 450 observations.\n",
    "<div class=\"exercise\"><b>Exercise</b></div>\n",
    "\n",
    "Using the above file, try your best to predict **housing prices. (the 14th column)** We have provided a test set `data/boston_housing_test.csv` but refrain from looking at the file or evaluating on it until you have finalized and trained a model.\n",
    "1. Load in the data. It is tab-delimited. Quickly look at a summary of the data to familiarize yourself with it and ensure nothing is too egregious.\n",
    "2. Use a previously-discussed function to automatically partition the data into a training and validation (aka development) set. It is up to you to choose how large these two portions should be.\n",
    "3. Train a basic model on just a subset of the features. What is the performance on the validation set?\n",
    "4. Train a basic model on all of the features. What is the performance on the validation set?\n",
    "5. Toy with the model until you feel your results are reasonably good.\n",
    "6. Perform cross-validation with said model, and measure the average performance. Are the results what you expected? Were the average results better or worse than that from your original 1 validation set?\n",
    "7. Experiment with other models, and for each, perform 10-fold cross-validation. Which model yields the best average performance? Select this as your final model.\n",
    "8. Use this model to evaulate your performance on the testing set. What is your performance (MSE)? Is this what you expected?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}