week-5-lasso-assignment-1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression Week 5: Feature Selection and LASSO (Interpretation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, you will use LASSO to select features, building on a pre-implemented solver for LASSO (using GraphLab Create, though you can use other solvers). You will:\n",
    "* Run LASSO with different L1 penalties.\n",
    "* Choose best L1 penalty using a validation set.\n",
    "* Choose best L1 penalty using a validation set, with additional constraint on the size of subset.\n",
    "\n",
    "In the second notebook, you will implement your own LASSO solver, using coordinate descent. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fire up Graphlab Create"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#import graphlab\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import sklearn\n",
    "from sklearn.linear_model import LinearRegression, Ridge, Lasso\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load in house sales data\n",
    "\n",
    "Dataset is from house sales in King County, the region where the city of Seattle, WA is located."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "dtype_dict = {'bathrooms':float, 'waterfront':int, 'sqft_above':int, 'sqft_living15':float, 'grade':int, 'yr_renovated':int, 'price':float, 'bedrooms':float, 'zipcode':str, 'long':float, 'sqft_lot15':float, 'sqft_living':float, 'floors':float, 'condition':int, 'lat':float, 'date':str, 'sqft_basement':int, 'yr_built':int, 'id':str, 'sqft_lot':int, 'view':int}\n",
    "sales = pd.read_csv('kc_house_data.csv', dtype= dtype_dict)\n",
    "testing = pd.read_csv('wk3_kc_house_test_data.csv', dtype=dtype_dict)\n",
    "training = pd.read_csv('wk3_kc_house_train_data.csv', dtype=dtype_dict)\n",
    "validation = pd.read_csv('wk3_kc_house_valid_data.csv', dtype=dtype_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create new features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in Week 2, we consider features that are some transformations of inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from math import log, sqrt\n",
    "sales['sqft_living_sqrt'] = sales['sqft_living'].apply(sqrt)\n",
    "sales['sqft_lot_sqrt'] = sales['sqft_lot'].apply(sqrt)\n",
    "sales['bedrooms_square'] = sales['bedrooms']*sales['bedrooms']\n",
    "sales['floors_square'] = sales['floors']*sales['floors']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Squaring bedrooms will increase the separation between not many bedrooms (e.g. 1) and lots of bedrooms (e.g. 4) since 1^2 = 1 but 4^2 = 16. Consequently this variable will mostly affect houses with many bedrooms.\n",
    "* On the other hand, taking square root of sqft_living will decrease the separation between big house and small house. The owner may not be exactly twice as happy for getting a house that is twice as big."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Learn regression weights with L1 penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us fit a model with all the features available, plus the features we just created above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_features = ['bedrooms', 'bedrooms_square',\n",
    "            'bathrooms',\n",
    "            'sqft_living', 'sqft_living_sqrt',\n",
    "            'sqft_lot', 'sqft_lot_sqrt',\n",
    "            'floors', 'floors_square',\n",
    "            'waterfront', 'view', 'condition', 'grade',\n",
    "            'sqft_above',\n",
    "            'sqft_basement',\n",
    "            'yr_built', 'yr_renovated']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Applying L1 penalty requires adding an extra parameter (`l1_penalty`) to the linear regression call in GraphLab Create. (Other tools may have separate implementations of LASSO.)  Note that it's important to set `l2_penalty=0` to ensure we don't introduce an additional L2 penalty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Lasso(alpha=500.0, copy_X=True, fit_intercept=True, max_iter=1000,\n",
       "   normalize=True, positive=False, precompute=False, random_state=None,\n",
       "   selection='cyclic', tol=0.0001, warm_start=False)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_all = Lasso(alpha=5e2, normalize=True) # set parameters\n",
    "model_all.fit(sales[all_features], sales['price']) # learn weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find what features had non-zero weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('bedrooms', 0.0)\n",
      "('bedrooms_square', 0.0)\n",
      "('bathrooms', 0.0)\n",
      "('sqft_living', 134.43931395541438)\n",
      "('sqft_living_sqrt', 0.0)\n",
      "('sqft_lot', 0.0)\n",
      "('sqft_lot_sqrt', 0.0)\n",
      "('floors', 0.0)\n",
      "('floors_square', 0.0)\n",
      "('waterfront', 0.0)\n",
      "('view', 24750.00458560952)\n",
      "('condition', 0.0)\n",
      "('grade', 61749.103090708129)\n",
      "('sqft_above', 0.0)\n",
      "('sqft_basement', 0.0)\n",
      "('yr_built', -0.0)\n",
      "('yr_renovated', 0.0)\n"
     ]
    }
   ],
   "source": [
    "for x in zip(all_features, model_all.coef_): \n",
    "    print x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that a majority of the weights have been set to zero. So by setting an L1 penalty that's large enough, we are performing a subset selection. \n",
    "\n",
    "***QUIZ QUESTION***:\n",
    "According to this list of weights, which of the features have been chosen? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Selecting an L1 penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find a good L1 penalty, we will explore multiple values using a validation set. Let us do three way split into train, validation, and test sets:\n",
    "* Split our sales data into 2 sets: training and test\n",
    "* Further split our training data into two sets: train, validation\n",
    "\n",
    "Be *very* careful that you use seed = 1 to ensure you get the same answer!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#(training_and_validation, testing) = sales.random_split(.9,seed=1) # initial train/test split\n",
    "#(training, validation) = training_and_validation.random_split(0.5, seed=1) # split training into train and validate\n",
    "testing['sqft_living_sqrt'] = testing['sqft_living'].apply(sqrt)\n",
    "testing['sqft_lot_sqrt'] = testing['sqft_lot'].apply(sqrt)\n",
    "testing['bedrooms_square'] = testing['bedrooms']*testing['bedrooms']\n",
    "testing['floors_square'] = testing['floors']*testing['floors']\n",
    "\n",
    "training['sqft_living_sqrt'] = training['sqft_living'].apply(sqrt)\n",
    "training['sqft_lot_sqrt'] = training['sqft_lot'].apply(sqrt)\n",
    "training['bedrooms_square'] = training['bedrooms']*training['bedrooms']\n",
    "training['floors_square'] = training['floors']*training['floors']\n",
    "\n",
    "validation['sqft_living_sqrt'] = validation['sqft_living'].apply(sqrt)\n",
    "validation['sqft_lot_sqrt'] = validation['sqft_lot'].apply(sqrt)\n",
    "validation['bedrooms_square'] = validation['bedrooms']*validation['bedrooms']\n",
    "validation['floors_square'] = validation['floors']*validation['floors']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we write a loop that does the following:\n",
    "* For `l1_penalty` in [10^1, 10^1.5, 10^2, 10^2.5, ..., 10^7] (to get this in Python, type `np.logspace(1, 7, num=13)`.)\n",
    "    * Fit a regression model with a given `l1_penalty` on TRAIN data. Specify `l1_penalty=l1_penalty` and `l2_penalty=0.` in the parameter list.\n",
    "    * Compute the RSS on VALIDATION data (here you will want to use `.predict()`) for that `l1_penalty`\n",
    "* Report which `l1_penalty` produced the lowest RSS on validation data.\n",
    "\n",
    "When you call `linear_regression.create()` make sure you set `validation_set = None`.\n",
    "\n",
    "Note: you can turn off the print out of `linear_regression.create()` with `verbose = False`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(398213327300134.94, 10.0)\n",
      "(399041900253346.81, 31.622776601683793)\n",
      "(429791604072559.62, 100.0)\n",
      "(463739831045121.06, 316.22776601683796)\n",
      "(645898733633800.88, 1000.0)\n",
      "(1222506859427162.8, 3162.2776601683795)\n",
      "(1222506859427162.8, 10000.0)\n",
      "(1222506859427162.8, 31622.776601683792)\n",
      "(1222506859427162.8, 100000.0)\n",
      "(1222506859427162.8, 316227.76601683791)\n",
      "(1222506859427162.8, 1000000.0)\n",
      "(1222506859427162.8, 3162277.6601683795)\n",
      "(1222506859427162.8, 10000000.0)\n"
     ]
    }
   ],
   "source": [
    "l1_penalty_list = np.logspace(1,7,13) \n",
    "\n",
    "X_train = training[all_features]\n",
    "y_train = training['price']\n",
    "X_valid = validation[all_features ]\n",
    "y_valid = validation['price']\n",
    "X_test  = testing[all_features ]\n",
    "y_test  = testing['price']\n",
    "\n",
    "RSS_list = []\n",
    "model_list =[]\n",
    "for l1_penalty in l1_penalty_list :\n",
    "    model = Lasso(alpha= l1_penalty, normalize=True) # set parameters\n",
    "    model.fit(X_train, y_train)\n",
    "    error = model.predict(X_valid) - y_valid\n",
    "    RSS_list.append(np.dot(error, error) )\n",
    "    model_list.append(model)\n",
    "    \n",
    "for x in zip(RSS_list, l1_penalty_list ):\n",
    "    print x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10eef4d50>"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEaCAYAAAD+E0veAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYXFWZx/Hvr7Pv+wIkIYGEhBAYloALKgRFCLKNgyzj\nBqKoowzjqDOMG4LjNm64MDLIKCigg46jQQM4YhBGBAkgSFIJhIQlUE32dGdPd7/zx73dFE26u7q7\nqm9X9e/zPPV017237n1Pbe8959w6RxGBmZkZQE3WAZiZWe/hpGBmZi2cFMzMrIWTgpmZtXBSMDOz\nFk4KZmbWwknBOk3SdEkhqX96/3ZJ7y5m2y4c6xOSru9OvNWou89rFiox5r7ISaEPknSHpKv2sfws\nSbWd/dBGxMKIuLEEcZ0oaW2rfX8hIt7b3X3v41gXSmqUtE1SnaRHJZ3eapuLJa2QVC/pRUmLJY1I\n102R9N+SNkjaKulxSReWOs5iSTpX0n2Sdki6O6s4OkPS3ZJK/tpa9zgp9E03Au+QpFbL3wncHBEN\nGcSUhT9GxHBgNPDvwE8kjQaQdALwBeCCiBgBHAr8V8FjfwQ8BxwIjCN57l7swdhb2wRcDXwpwxis\nCjgp9E2/IPkie33zAkljgNOBH6b33yLpkfQs+jlJn21rZ4VnfJL6Sfpqega9GnhLq20vkpRLz75X\nS3p/unwYcDuwf3r2vk3S/pI+K+mmgsefKWmZpC3pcQ8tWPe0pI9Jeiw9e/8vSYM7ejIioonkS34Y\nMCtdfCxJ0ngk3WZTRNwYEfUF62+IiO0R0RARj0TE7W08P7nCWoik/pLWSzpa0mBJN0namJbpQUmT\nOop5H2X4bUTcCrzQ0bbNNbK0aW5D+ry9vWD9oPQ1fDatIV0raUirx35U0jpJeUkXFTy2qPeNpM+T\nvP++k77W35F0jaSvtdpukaSPdPb5sK5zUuiDImIncCvwroLF5wIrIuLR9P72dP1oki/2D0o6u4jd\nv48kuRwFzAfOabV+Xbp+JHAR8A1JR0fEdmAh8EJEDE9vL/uCk3QI8GPgH4AJwGLgNkkDW5XjVGAG\ncARwYUcBS+qXxrIXeCZd/ABwiqQrJR0vaVCrh90PXCPpfEnTOjjEj4ELCu6fAmyIiIeBdwOjgKkk\nifoDwM6OYi6BycB44IA0huskzU7XfQk4BDgSmJlu85lWjx2VLr+Y5HkYk64r6n0TEZ8E7gU+nL7W\nHyapwV4gqQZA0njgTcAtpSq0dawik4Kk76dnKY8Xse0bJD0sqUHSOa3WNUr6c3pbVL6Ie6UbgXMK\nzqTflS4DICLujoi/RERTRDxG8sV2QhH7PRe4OiKei4hNwBcLV0bEryPiqUj8HvgNBTWWDpwH/Doi\n/jci9gJfBYYAry3Y5lsR8UJ67NtIvtja8mpJW4Bd6b7eERHr0jjvBd4KHA38Gtgo6etpAgF4G8mX\n2qeBNel76Ng2jnMLcKakoen9vyV5PiFJROOAmRHRGBEPRURdcU9Ht306Inanr8OvgXPTJsVLgI+k\ntaN6kma08wsetxe4KiL2RsRiYBswG7r1viEi/gRsBd6YLjofuDsismyW63MqMikAN5CcDRbjWZKz\nxX2dbeyMiCPT25kliq0iRMT/ARuAsyUdDBxHwXMk6VWSlqTNHFtJzmDHF7Hr/Una2ps9U7hS0kJJ\n90valH4hn1bkfpv33bK/tNnnOZIz1ma1Bf/vAIa3s7/7I2I0MAZYRKvkFBG3R8QZwFjgLJL30XvT\ndZsj4vKIOAyYBPwZ+MU++mmIiFVADjgjTQxn8tJz/SPgTpL+jBck/ZukAe0+C6WxOa2dNXuG5Pmd\nAAwFHkqbs7YAd6TLm21s1e/U8jx3433T7EbgHen/7yB5fqwHVWRSiIh7SDrWWkg6WMlVNQ9JulfS\nnHTbp9MzlqYsYu3lfkhSQ3gHcGerM7JbSL4op0bEKOBa4BVfePuQJ2kKadbStJI2wfw3yVn5pPQL\neXHBfjsasvcFko7d5v0pPdbzRcTVpojYBnwQeKeko/axviki7gJ+B8zbx/oNJGXanySB7EtzE9JZ\nwPI0UZCebV8ZEXNJajyn8/JmvXIZk/bjNJtG8vxuIGm+OiwiRqe3UWmHfDE6877Z1+t9E3CWpL8i\n6dz/RZHHtRKpyKTQhuuASyPiGOBjJFeTdGSwpKXpmWsx7eXV5ockbbbvo6DpKDUC2BQRuyQdR9Lk\nUYxbgb9XcsnmGODygnUDgUHAeqBB0kLgzQXrXwTGSRrVzr7fIumN6dn0R4HdwH1FxtamtLnpetK2\ncyWX554vaYwSx5E0g9yfrv+ypHlpp/EIkqSyKiI2tnGIn5CU9YO8vEa2QNLhabNUHUnTTKdPYJR0\n8A8G+gM1SjqwO6pxXClpoKTXkySjn6a1r++R9PVMTPd9gKRTigylM++bF4GDChdExFrgQZIawn+n\n/V/Wg6oiKUgaTnKW9VNJfwb+A9iviIceGBHzSd64V6fNKH1GRDxN8oU6jOTsrtDfAVdJqif5ory1\nyN1+j6Q55FHgYeDnBcerB/4+3ddmkud9UcH6FSRn1KvTpov9W8W7kqRW822SM9ozgDMiYk+RsXXk\nauA0SUek8b0PeJLky/om4CsRcXO67VDgf4AtwGqSGkybTZARkQf+SPI+Lby0dTLws/QYOeD3pE0m\nSq76ubbI2N9Jcob/XZJmsJ0kr0VbatMyvgDcDHwgff4B/hlYBdwvqQ74LWmfQRE68775Jkm/1mZJ\n3ypYfiNwOG46yoQqdZIdSdOBX0XEPEkjgZUR0WYikHRDuv3PurLerFpIOhG4KSKmZB3Lvkh6A0kS\nPjAq9QuqglVFTSG9WmONpLdB0tactkm2KW0WGJT+Px44Hlhe9mDNrE1pk9dlwPVOCNmoyKQg6cck\nVfHZSn5IczHwduBiSY8Cy0g69JB0rJKhE94G/IekZeluDgWWptsvAb4UEU4KZhlR8kPELSRNv1dn\nHE6fVbHNR2ZmVnoVWVMwM7PycFIwM7MWFTeu+fjx42P69OlZh2FmVlEeeuihDRExoaPtKi4pTJ8+\nnaVLl2YdhplZRZH0TMdbufnIzMwKOCmYmVkLJwUzM2vhpGBmZi2cFMzMrIWTgpmZtai4S1LNOmvb\n7gYeWL2RJo/oYhVu+rihzJo0oqzHcFKwqvfdu1dxzZKnsg7DrNs+cMLBXL5wTlmP4aRgVe+xtVs5\nZNJwvn7ukVmHYtYt44cPKvsxypYUJH2fZIq/dRHxinltJb2dZIYnAfXAByPi0XLFY31XLl/PibMn\nMO+Atmb5NLNm5exovgE4tZ31a4ATIuJw4HMkcyybldT6+t1s2LabQ/cbmXUoZhWhbDWFiLgnnTKz\nrfWFk63fD/TKqQGtsuXydQAcul95O+fMqkVvuST1YuD2rIOw6tOcFOa6pmBWlMw7miUtIEkKr2tn\nm0uASwCmTZvWQ5FZNcjl69hv1GBGDx2YdShmFSHTmoKkI4DrgbMiYmNb20XEdRExPyLmT5jQ4XDg\nZi1y+Xr3J5h1QmZJQdI04OfAOyPiiazisOq1u6GRp9Zvc3+CWSeU85LUHwMnAuMlrQWuAAYARMS1\nwGeAccC/SwJoiIj55YrH+p4nX9xGQ1O4pmDWCeW8+uiCDta/F3hvuY5v9tKVR04KZsXqLVcfmZXc\nitp6Bg+oYfq4YVmHYlYxnBSsauXydcyePJJ+Nco6FLOK4aRgVSkiyOXrmOtOZrNOcVKwqvRi3W42\n79jr/gSzTnJSsKrkTmazrnFSsKq0PE0Kcya7+cisM5wUrCrl8nVMHTuEEYMHZB2KWUVxUrCqlMvX\ncehkNx2ZdZaTglWdXXsbWbNhu/sTzLrAScGqzsraeprCcyiYdYWTglUdX3lk1nVOClZ1cvk6hg3s\nx9QxQ7MOxaziOClY1cnl65mz30hqPLyFWac5KVhViQhytXXuTzDrIicFqyprN++kfleD+xPMushJ\nwaqKO5nNusdJwapKLl+P5OEtzLrKScGqSi5fx/Rxwxg6sGyTCppVNScFqyruZDbrHicFqxrbdjfw\nzMYdHvPIrBucFKxqrKx1J7NZdzkpWNVYnq8H4ND9nRTMuspJwapGLl/HyMH92X/U4KxDMatYTgpW\nNVbk6zh0v5FIHt7CrKucFKwqNDUFK2rr3Z9g1k1OClYVnt20gx17GpnrpGDWLU4KVhU8vIVZaTgp\nWFXI5evoVyNmTRqedShmFc1JwarC8nw9B40fxuAB/bIOxayiOSlYVcilVx6ZWfc4KVjF27pzL89v\n2emkYFYCTgpW8Va0dDJ7IDyz7nJSsIrnK4/MSsdJwSpeLl/P2GEDmThiUNahmFU8JwWreM1zKHh4\nC7Puc1KwitbQ2MTK2nrPoWBWIk4KVtGe3rid3Q1N7k8wKxEnBatoLXMoOCmYlYSTglW0XL6OAf3E\nzIke3sKsFMqWFCR9X9I6SY+3sV6SviVplaTHJB1drliseuXydRw8YTgD+/v8xqwUyvlJugE4tZ31\nC4FZ6e0S4LtljMWqVC5f5+GyzUqobEkhIu4BNrWzyVnADyNxPzBa0n7liseqz6bte3ixbrf7E8xK\nKMs69wHAcwX316bLzIriXzKblV5FNMRKukTSUklL169fn3U41kvkPOaRWcllmRSeB6YW3J+SLnuF\niLguIuZHxPwJEyb0SHDW+y3P1zFxxCDGDffwFmalkmVSWAS8K70K6dXA1ojIZxiPVZhcvt5NR2Yl\n1r9cO5b0Y+BEYLyktcAVwACAiLgWWAycBqwCdgAXlSsWqz57GppYta6eEw5xzdGslMqWFCLigg7W\nB/Chch3fqtvqDdvY2xjuTzArsYroaDZrrbmT2b9RMCstJwWrSLl8PQP71zBj/LCsQzGrKk4KVpFy\n+TpmTxpB/35+C5uVkj9RVpFy+Tr3J5iVgZOCVZx19bvYsG2PL0c1K4Oik4KkcZL+WtIx5QzIrCM5\nz6FgVjZtJgVJv5I0L/1/P+Bx4D3AjyT9Qw/FZ/YKLcNbeApOs5Jrr6YwIyKa50K4CPjfiDgDeBVJ\ncjDLRC5fx/6jBjNq6ICsQzGrOu0lhb0F/7+R5BfIREQ90FTOoMzak3Qyu5ZgVg7tJYXnJF0q6a+B\no4E7ACQNIR2uwqyn7drbyFPrtzspmJVJe0nhYuAw4ELgvIjYki5/NfCDMsdltk+r1m2jsSmcFMzK\npM2xjyJiHfCBwmWSxgB3R8SScgdmti/LPYeCWVm1d/XRZyTNSf8fJGkJ8BTwoqQ39VSAZoVy+TqG\nDOjHgeM8vIVZObTXfHQesDL9/93p3wnACcAXyhmUWVty+TpmTx5BvxplHYpZVWovKexJh7cGOAX4\nSUQ0RkSOMg65bdaWiPDEOmZl1l5S2C1pnqQJwALgNwXrhpY3LLNXym/dxdade5nr/gSzsmnvjP8y\n4GckTUbfiIg1AJJOAx7pgdjMXqbll8yuKZiVTXtXHz0AzNnH8sWkP2Qz60nNSWGOk4JZ2bQ7IJ6k\nfpLGF9wfKOkSSbnyh2b2crl8PdPGDmX4IHdpmZVLe5ekng9sAh6T9HtJbwZWAwuBt/dQfGYtPIeC\nWfm1d8r1KeCYiFgl6Wjgj8A5EXFbz4Rm9pIdexpYs3E7Zx65f9ahmFW1ji5JXQUQEQ8DTzohWFZW\n1tYT4U5ms3Jrr6YwUdI/FtwfXXg/Ir5evrDMXm5FbTKxzlwnBbOyai8pfA8Y0c59sx6Ty9cxYlB/\npowZknUoZlWtvUtSr+zJQMzak8vXMWe/EUge3sKsnIqeo9ksKxHBCg9vYdYjnBSs11u7eSf1uxuc\nFMx6gJOC9XrLPbyFWY/p8KehkgYBfwNML9w+Iq4qX1hmL8nl66gRzJ7k6xzMyq2Y8QJ+CWwFHgJ2\nlzccs1fK5euYPn4YQwb2yzoUs6pXTFKYEhGnlj0Sszbk8vUcfsCorMMw6xOK6VO4T9LhZY/EbB/q\nd+3l2U07POaRWQ8ppqbwOuBCSWtImo8EREQcUdbIzEiGtwB3Mpv1lGKSwsKyR2HWBk+sY9azOmw+\niohngNHAGeltdLrMrOyW5+sZNWQA+40anHUoZn1Ch0lB0mXAzcDE9HaTpEvLHZgZvDSHgoe3MOsZ\nxTQfXQy8KiK2A0j6MsncCt8uZ2BmjU3Bytp6zj9uatahmPUZxVx9JKCx4H5jusysrJ7ZuJ2dexvd\nn2DWg4qpKfwAeEDS/6T3zwb+s3whmSVyec+hYNbTiulo/jpwEcl8zZuAiyLi6mJ2LulUSSslrZJ0\n+T7WT5O0RNIjkh6TdFpnC2DVK5evo1+NmDlxeNahmPUZbdYUJI2MiDpJY4Gn01vzurERsam9HUvq\nB1wDnAysBR6UtCgilhds9ing1oj4rqS5wGKSMZbMyOXrOHjCMAYP8PAWZj2lveajW4DTScY8ioLl\nSu8f1MG+jwNWRcRqAEk/Ac4CCpNCAM1tA6OAF4qO3KpeLl/HsTPGZh2GWZ/S3sxrp6d/Z3Rx3wcA\nzxXcXwu8qtU2nwV+k17iOgx4UxePZVVmy449vLB1lzuZzXpYMb9TuKuYZV10AXBDREwBTgN+JOkV\nMUm6RNJSSUvXr19fokNbb9bcyeykYNaz2kwKkgan/QnjJY2RNDa9TSepBXTkeaDwAvMp6bJCFwO3\nAkTEH4HBwPjWO4qI6yJifkTMnzBhQhGHtkr30vAWHgjPrCe1V1N4P0l/wpz0b/Ptl8B3itj3g8As\nSTMkDQTOBxa12uZZ4I0Akg4lSQquChgrausYP3wgE0d4eAuzntRen8I3gW9KujQiOv3r5YhokPRh\n4E6gH/D9iFgm6SpgaUQsAj4KfE/SR0g6nS+MiGh7r9ZX5PL1bjoyy0CHP16LiG9LmgfMJTmTb17+\nwyIeu5jkMtPCZZ8p+H85cHxnArbq19DYxMoX67nwtdOzDsWszylmjuYrgBNJksJikqG0/w/oMCmY\ndcWaDdvZ09Dk/gSzDBQz9tE5JO3+tRFxEfBXJL8pMCuL5Z5DwSwzxSSFnRHRBDRIGgms4+VXFZmV\nVC5fz8B+NRw8wcNbmPW0YgbEWyppNPA9kquPtpEMnW1WFrl8HTMnDmdAv2LOWcyslIrpaP679N9r\nJd0BjIyIx8oblvVluXwdr5/l36OYZaG9AfGObm9dRDxcnpCsL9u4bTfr6ne7k9ksI+3VFL6W/h0M\nzAceJRkM7whgKfCa8oZmfZGHtzDLVpuNthGxICIWAHng6HSYiWOAo3jlcBVmJZHzlUdmmSqmJ292\nRPyl+U5EPA4cWr6QrC/L5euYNHIQY4cNzDoUsz6pmKuPHpN0PXBTev/tgDuarSyW5+tcSzDLUDE1\nhYuAZcBl6W15usyspPY0NPHU+m1OCmYZKuaS1F3AN9KbWdmsWreNvY3hpGCWofYuSb01Is6V9Bde\nPh0nABFxRFkjsz6nuZN5ri9HNctMezWFy9K/p/dEIGa5fB2D+tcwfdywrEMx67Pam08hn/59pufC\nsb4sV1vH7Mkj6O/hLcwy017zUT37aDYi+QFbRIQbfq1kIoJcvp6TD52UdShmfVp7NQU37FqPWVe/\nm03b93h4C7OMFfM7BQAkTeTlM689W5aIrE/yHApmvUOHjbeSzpT0JLAG+D3wNHB7meOyPqb5yqM5\nTgpmmSqmR+9zwKuBJyJiBsksbPeXNSrrc3L5eg4YPYRRQwZkHYpZn1ZMUtgbERuBGkk1EbGEZNRU\ns5JZ4eEtzHqFYvoUtkgaDtwD3CxpHbC9vGFZX7JrbyOrN2xn4bzJWYdi1ucVU1M4C9gJfAS4A3gK\nOKOcQVnf8uSL22hs8vAWZr1Be79TuAa4JSL+ULD4xvKHZH2N51Aw6z3aqyk8AXxV0tOS/k3SUT0V\nlPUty/N1DBvYj2ljh2Ydilmf197Ma9+MiNcAJwAbge9LWiHpCkmH9FiEVvVy+WR4i5oaZR2KWZ/X\nYZ9CRDwTEV+OiKOAC4CzgVzZI7M+IRnewlcemfUWxfx4rb+kMyTdTPKjtZXAW8semfUJL2zdRd2u\nBicFs16ivY7mk0lqBqcBfwJ+AlwSEb4c1Uom90JzJ7PHPDLrDdr7ncK/ALcAH42IzT0Uj/UxzVce\nzZ7smoJZb9DeKKkn9WQg1jflaus4cNxQhg8qemxGMysjz2ZimdnT0MQjz27hUNcSzHoNJwXLzI33\nPU1+6y7OO25q1qGYWcpJwTKxvn4337rrSU6aM5EFsydmHY6ZpZwULBNfvXMluxoa+dRbDs06FDMr\n4KRgPe4va7dy60PPcdHxMzhowvCswzGzAk4K1qMigs/etoxxwwZy6Ukzsw7HzFpxUrAetejRF3jo\nmc380ylzGDHYs6yZ9TZOCtZjduxp4IuLV3DElFGcc8yUrMMxs30oa1KQdKqklZJWSbq8jW3OlbRc\n0jJJt5QzHsvWd+9+itq6XVxxxmEeEdWslyrbz0gl9QOuAU4G1gIPSloUEcsLtplFMpzG8RGxWZKv\nTaxSz23awX/cs5qzj9yfYw4ck3U4ZtaGctYUjgNWRcTqiNhDMqDeWa22eR9wTfPYShGxrozxWIY+\n/+sc/SQuX+hLUM16s3ImhQOA5wrur02XFToEOETSHyTdL+nUMsZjGblv1QbuWFbLhxYczORRg7MO\nx8zakfUoZP2BWcCJwBTgHkmHR8SWwo0kXQJcAjBt2rSejtG6oaGxiStvW87UsUN47+sPyjocM+tA\nOWsKzwOFg9pMSZcVWgssioi9EbGGZF7oWa13FBHXRcT8iJg/YcKEsgVspXfLn55l5Yv1fPK0uQwe\n0C/rcMysA+VMCg8CsyTNkDQQOB9Y1GqbX5DUEpA0nqQ5aXUZY7IetHn7Hr72myd47cHjOOWwSVmH\nY2ZFKFtSiIgG4MPAnSRzOt8aEcskXSXpzHSzO4GNkpYDS4CPR8TGcsVkPesbv32CbbsbuOKMw5B8\nCapZJShrn0JELAYWt1r2mYL/A/jH9GZVZEVtHTfd/wzvfPWBzJ7sqTbNKoV/0WwlFxFcuWg5I4cM\n4CMnH5J1OGbWCU4KVnJ3Lqvlj6s38tGTD2H00IFZh2NmneCkYCW1a28j//rrHHMmj+CC43z5sFml\ncVKwkrr+3tWs3byTz5wxl/79/PYyqzT+1FrJ5Lfu5JolT7Fw3mRee/D4rMMxsy5wUrCS+dLtK2iM\n4BOneXwjs0rlpGAlsfTpTfzyzy/w/jccxNSxQ7MOx8y6yEnBuq2pKbjytuXsN2owHzzx4KzDMbNu\ncFKwbvvpQ8/xl+e3cvnCOQwdmPUYi2bWHU4K1i11u/bylTtXMv/AMZz5V/tnHY6ZdZNP66xbvn3X\nk2zcvocbLjrO4xuZVQHXFKzLnlq/jR/84WnOmz+VeQeMyjocMysBJwXrsn/91XKGDOjHx06ZnXUo\nZlYiTgrWJUtWrGPJyvVc9qZZjB8+KOtwzKxEnBSs0/Y0NPG5Xy3noAnDeNdrpmcdjpmVkJOCddqN\n9z3N6g3b+fTpcxnY328hs2riT7R1yvr63Xzrric5ac5EFsyemHU4ZlZiTgrWKV+5cwW7Ghr51Fs8\nvpFZNXJSsKI9tnYLP31oLRcdP4ODJgzPOhwzKwMnBStKRDK+0bhhA7n0pJlZh2NmZeKkYEX55Z9f\n4KFnNvNPp8xhxOABWYdjZmXipGAd2r67gS/enuOIKaM455gpWYdjZmXkpGAd+u7dT/Fi3W6uOOMw\namo8vpFZNXNSsHY9u3EH1927mrOP3J9jDhyTdThmVmZOCtauLyzO0b9GXL7Ql6Ca9QVOCtam+1Zt\n4I5ltXxowUwmjxqcdThm1gOcFGyfNm3fw5W3LWfq2CFc/LoZWYdjZj3Ek+xYi9qtu7hzWS13PF7L\nA2s20hRw3TuPYfCAflmHZmY9xEmhj3t6w/YkESyr5ZFntwAwa+JwPrRgJgvn7cfc/UdmHKGZ9SQn\nhT4mIlj5Yj13PJ7UCFbU1gNw+AGj+PgpsznlsMnMnOghLMz6KieFPiAieHTtVu54vJY7l9WyZsN2\nJDj2wLF8+vS5nHLYJKaMGZp1mGbWCzgpVKnGpuBPazZx57IkEeS37qJ/jXjNweN47+tn8Oa5k5kw\nwjOmmdnLOSlUkd0Njdz31EbufLyW/13+Ihu372FQ/xpOOGQCHz9lNm+cM4lRQz1ukZm1zUmhwu3Y\n08DvV67njmW1/C63jvrdDQwf1J+T5kzk1HmTOXH2BIYO9MtsZsXxt0UF2rpjL79b+SJ3PF7L759Y\nz669TYwZOoDTDt+PU+dN5rUzxzGovy8jNbPOc1KoABHBitp6lqxcx5IV63jomc00BUweOZjz5k/l\nlHmTOW76WPr3828Rzax7nBR6qe27G/jDqg0sWbmeu1euI791FwDzDhjJhxbMZMGciRw5ZbRHLTWz\nknJS6EXWbNjO71as4+6V63hg9Sb2NDYxfFB/XjdzPB9500ROmD2BSSM9BpGZlY+TQoZ2NzTywOpN\nLc1CT2/cAcDMicN592sPZMGcicw/cCwD+7tZyMx6RlmTgqRTgW8C/YDrI+JLbWz3N8DPgGMjYmk5\nY8raC1t2pklgPX9YtYGdexsZ1L+G1x48jve8bgYLZk9k6lj/kMzMslG2pCCpH3ANcDKwFnhQ0qKI\nWN5quxHAZcAD5YolSw2NTTz87JaWZqHmYSUOGD2Ec46ZwklzJvLqg8YxZKCvFjKz7JWzpnAcsCoi\nVgNI+glwFrC81XafA74MfLyMsfSoDdt28/uV61mych33PLGeul0N9K8Rx04fyydOm8OC2ROZOXE4\nkjuJzax3KWdSOAB4ruD+WuBVhRtIOhqYGhG/llTWpPDn57Zw431P0xRBU0BTBBFBUxMtyyKi1frm\ndYXr97F900vb721sYs3G7UTAhBGDOHXeZBbMnsjxs8YzcrB/TWxmvVtmHc2SaoCvAxcWse0lwCUA\n06ZN69LxNm/fw0PPbKZGUCOh9O/L/q9pXqeW7WrES/dralo9dh/b18DZRx3ASXMmMne/kb5k1Mwq\nSjmTwvPA1IL7U9JlzUYA84C702aUycAiSWe27myOiOuA6wDmz58fXQlmwZyJLJgzsSsPNTPrM8p5\nreODwCz2K2fVAAAGZklEQVRJMyQNBM4HFjWvjIitETE+IqZHxHTgfuAVCcHMzHpO2ZJCRDQAHwbu\nBHLArRGxTNJVks4s13HNzKzrytqnEBGLgcWtln2mjW1PLGcsZmbWMf9U1szMWjgpmJlZCycFMzNr\n4aRgZmYtnBTMzKyFIrr0W7DMSFoPPFOwaBSwtY3NW68bD2zoxuHbO1YpHtvWNvta3p2ydKccxT6+\nEsrSnXK0ta4Sy1LqchQbT3ceWwnvr2If31NlmRURozrcKtIxgCr1BlxX7DpgabmOVYrHtrXNvpZ3\npyzdKUc1laU75aimspS6HNVUlr74WamG5qPburiu1McqxWPb2qY3laPYx1dCWbpTju4eu9T7602v\nSXf32ZvK0uc+KxXXfNQdkpZGxPys4ygFl6V3qpayVEs5wGXprGqoKXTGdVkHUEIuS+9ULWWplnKA\ny9IpfaqmYGZm7etrNQUzM2uHk4KZmbVwUjAzsxZ9OilIOkjSf0r6WdaxdJeksyV9T9J/SXpz1vF0\nlaRDJV0r6WeSPph1PN0laZikpZJOzzqW7pB0oqR709fmxKzj6Q5JNZI+L+nbkt6ddTxdJen16etx\nvaT7SrXfqksKkr4vaZ2kx1stP1XSSkmrJF0OEBGrI+LibCLtWCfL8ouIeB/wAeC8LOJtSyfLkYuI\nDwDnAsdnEW97OlOW1D8Dt/ZslMXpZFkC2AYMBtb2dKwd6WRZziKZHngvvawsnfys3Jt+Vn4F3Fiy\nILrza73eeAPeABwNPF6wrB/wFHAQMBB4FJhbsP5nWcddwrJ8DTg669i7Uw7gTOB24G+zjr07ZQFO\nJpmG9kLg9Kxj72ZZatL1k4Cbs469m2W5HHh/uk2v+ux38TN/KzCiVDFUXU0hIu4BNrVafBywKpKa\nwR7gJyRnC71aZ8qixJeB2yPi4Z6OtT2dfU0iYlFELATe3rORdqyTZTkReDXwt8D7JPWqz1tnyhIR\nTen6zcCgHgyzKJ18XdaSlAOgseei7FhnPyuSpgFbI6K+VDGUdTrOXuQA4LmC+2uBV0kaB3weOErS\nv0TEFzOJrnP2WRbgUuBNwChJMyPi2iyC64S2XpMTgbeSfPEs3sfjeqN9liUiPgwg6UJgQ8EXa2/W\n1uvyVuAUYDTwnSwC64K2PivfBL4t6fXAPVkE1kltlQPgYuAHpTxYX0kK+xQRG0na4CteRHwL+FbW\ncXRXRNwN3J1xGCUVETdkHUN3RcTPgZ9nHUcpRMQOki/TihcRV5R6n72qOltGzwNTC+5PSZdVomop\nS7WUA1yW3qpaytKj5egrSeFBYJakGZIGknT+Lco4pq6qlrJUSznAZemtqqUsPVuOrHvby9B7/2Mg\nz0uXm12cLj8NeIKkF/+TWcfZl8pSLeVwWXrvrVrK0hvK4QHxzMysRV9pPjIzsyI4KZiZWQsnBTMz\na+GkYGZmLZwUzMyshZOCmZm1cFKwPkXStoL/75C0RdKvMo7pbknz0/8/kWUsZk4K1pd9BXhn1kG0\n4qRgmXJSsD4rIu4C2h1yOD2L/6akP0t6XNJx6fJh6YQof5L0iKTmoYwvlPTztBbypKR/K9jXd9NZ\n2JZJunIfx/oSMCQ91s2SrpL0DwXrPy/pslKV32xf+vQoqWZFGhoRR0p6A/B9YB7wSeB3EfEeSaOB\nP0n6bbr9kcBRwG5gpaRvR8RzJMMTbJLUD7hL0hER8VjzQSLickkfjogjASRNJxmZ9Op0LobzScbW\nNysbJwWzjv0YkglQJI1Mk8CbgTMlfSzdZjAwLf3/rojYCiBpOXAgyXj450q6hORztx/JLGCP0YaI\neFrSRklHkcx49kgkw72blY2TglnHWg8QFoCAv4mIlYUrJL2KpIbQrBHoL2kG8DHg2IjYLOkGkkTS\nketJpvOcTFJLMSsr9ymYdew8AEmvI5n6cCtwJ3CpJKXrjupgHyOB7cBWSZOAhW1st1fSgIL7/wOc\nChybHtOsrFxTsD5L0r3AHGC4pOZhivf1xbtL0iPAAOA96bLPAVcDj6Xt/WuA09s6VkQ8mu5jBUlT\n0h/a2PS6dJ8PR8TbI2KPpCXAlojoVfMJW3Xy0Nlm7ZB0N/CxiFia0fFrgIeBt0XEk1nEYH2Lm4/M\neilJc4FVJB3XTgjWI1xTMDOzFq4pmJlZCycFMzNr4aRgZmYtnBTMzKyFk4KZmbVwUjAzsxb/D5V2\n1InXPHhgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ebad110>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot( l1_penalty_list, RSS_list, '-' )\n",
    "plt.xscale('log')\n",
    "plt.xlabel('l1 penalty')\n",
    "plt.ylabel('Validation RSS')\n",
    "plt.title('Validation RSS vs. l1 penalty')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*** QUIZ QUESTION. *** What was the best value for the `l1_penalty`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l1_penalty_list[np.argmin(RSS_list )]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***QUIZ QUESTION***\n",
    "Also, using this value of L1 penalty, how many nonzero weights do you have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model = model_list[np.argmin(RSS_list )]\n",
    "np.count_nonzero(model.coef_) + np.count_nonzero(model.intercept_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Limit the number of nonzero weights\n",
    "\n",
    "What if we absolutely wanted to limit ourselves to, say, 7 features? This may be important if we want to derive \"a rule of thumb\" --- an interpretable model that has only a few features in them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section, you are going to implement a simple, two phase procedure to achive this goal:\n",
    "1. Explore a large range of `l1_penalty` values to find a narrow region of `l1_penalty` values where models are likely to have the desired number of non-zero weights.\n",
    "2. Further explore the narrow region you found to find a good value for `l1_penalty` that achieves the desired sparsity.  Here, we will again use a validation set to choose the best value for `l1_penalty`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "max_nonzeros = 7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring the larger range of values to find a narrow range with the desired sparsity\n",
    "\n",
    "Let's define a wide range of possible `l1_penalty_values`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "l1_penalty_values =  np.logspace(1, 4, num=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, implement a loop that search through this space of possible `l1_penalty` values:\n",
    "\n",
    "* For `l1_penalty` in `np.logspace(8, 10, num=20)`:\n",
    "    * Fit a regression model with a given `l1_penalty` on TRAIN data. Specify `l1_penalty=l1_penalty` and `l2_penalty=0.` in the parameter list. When you call `linear_regression.create()` make sure you set `validation_set = None`\n",
    "    * Extract the weights of the model and count the number of nonzeros. Save the number of nonzeros to a list.\n",
    "        * *Hint: `model['coefficients']['value']` gives you an SArray with the parameters you learned.  If you call the method `.nnz()` on it, you will find the number of non-zero parameters!* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RSS_list = []\n",
    "model_list =[]\n",
    "nonzero_coef =[]\n",
    "\n",
    "for l1_penalty in l1_penalty_values :\n",
    "    model = Lasso(alpha= l1_penalty, normalize=True) # set parameters\n",
    "    model.fit(X_train, y_train)\n",
    "    error = model.predict(X_valid) - y_valid\n",
    "    RSS_list.append(np.dot(error, error) )\n",
    "    model_list.append(model)\n",
    "    nonzero_coef.append( np.count_nonzero(model.coef_) + np.count_nonzero(model.intercept_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(10.0, 15)\n",
      "(14.384498882876629, 15)\n",
      "(20.691380811147901, 15)\n",
      "(29.763514416313178, 15)\n",
      "(42.813323987193932, 13)\n",
      "(61.584821106602639, 12)\n",
      "(88.586679041008225, 11)\n",
      "(127.42749857031335, 10)\n",
      "(183.29807108324357, 7)\n",
      "(263.66508987303581, 6)\n",
      "(379.26901907322497, 6)\n",
      "(545.55947811685144, 6)\n",
      "(784.75997035146065, 5)\n",
      "(1128.8378916846884, 3)\n",
      "(1623.776739188721, 3)\n",
      "(2335.7214690901214, 2)\n",
      "(3359.8182862837812, 1)\n",
      "(4832.9302385717519, 1)\n",
      "(6951.9279617756056, 1)\n",
      "(10000.0, 1)\n"
     ]
    }
   ],
   "source": [
    "for x in zip( l1_penalty_values,nonzero_coef) :\n",
    "    print x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1109ea610>"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEaCAYAAADzDTuZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYVOX5xvHvs4UttKUsdWkCIigCsiJFDdagYhcVldgJ\nqPmZxCRGjTH2RBO7gtgwdhRMDCoGFWyAumABKYJIk7IUadLh+f0xh2TdLLCzO7NnZ+b+XNdczJxz\n5rz37mHnmfeU95i7IyIiqS0t7AAiIhI+FQMREVExEBERFQMREUHFQEREUDEQERFUDESiYma3mdkq\nM1teRe21NjM3s4yqaC8WEjGzqBiIlJuZtQSuATq5e5OQMpxtZpPMbJOZTQwjQ7TMbKKZXRZ2Dtk7\nVW6ptszMAHP3XWFnCbQEVrt7cYgZ1gD3AQcAR4eYQ5KMegYpzswWmNlvzOxLM1tnZi+ZWXaJ+Zeb\n2TwzW2Nmr5lZsxLz3MyGmNlcM1trZg8HH+CY2RdmtrHEw82sbzCvZ/Dtdm2wXN8S65xoZreb2UfA\nJmA/M2sWtL0myHJ5OX6udDO73sy+MbMNZjbVzFoE83qb2afBz/upmfUu8b66ZvaEmS0zs++C3ULp\nZnYsMB5oFvw8I8toc5aZ9S/xOsPMVprZIWaWbWbPmtnq4Of+1MwaR7GpAHD3t919FLC0HL+Dvma2\nJPg9rAq29fkl5meZ2V/NbJGZrTCz4WaWU+q915hZcfD7uLjEe08ys8/MbL2ZLTazP+0hw+3AEcBD\nwe/toeD/yd9KLfeamf0q2t+HxJC765HCD2AB8AnQDKgPzAKGBPOOBlYBhwBZwIPA+yXe68BYII/I\nt+aVQL8y2hgMzAbqAM2B1cCJRL6MHBe8zg+WnQgsAg4k0nPNBN4HHgGyga5BO0fv4+f6LTAd6AAY\n0AVoEPyM3wODgvUPDF43CN73KvAoUBNoFPxufh7M6wss2UubfwSeK/H6JGBW8PznwL+AXCAd6A7U\nKcf2aR38njNKTb8MmLiP9/YFdgD3BNvvJ8APQIdg/r3Aa8HvpHaQ785S770l2AYnEinO9UrM7xxs\nw4OBFcBpZWUOtullJXL1IFLM0oLXDYN1Nw777yGVH6EH0CPk/wCRYnBBidd3AcOD508Ad5WYVwvY\nDrQOXjtweIn5o4Dfl1r/4UAxsH/w+lrgmVLLvAVcGDyfCNxSYl4LYCdQu8S0O4GR+/i55gCnljF9\nEPBJqWmTgYuAxsBWIKfEvIHAhOB5X/ZeDNoBG4Dc4PVzwB+D55cAk4CDo9w+P/pgLTE9mmJQs9Q2\nupFIgfwBaFtiXi/g2xLv3Vyy3WA79txDW/cB95aVuXQxCKbNAo4Lnl8FvBH230KqP7SbSABKnhmz\niciHPkR6Cwt3z3D3jUS+xTcvx3sJdsuMIvJB/3UwuRUwINhVstbM1hIpGE1LrGdxiefNgDXuvqHE\ntIWlMpSlBfBNGdN/9DOVWl8rIt+Cl5XI9iiRHsI+ufs8Ih9yJ5tZLnAK8Hww+xkiRe9FM1tqZneZ\nWWZ51ltJ37v7DyVeLyTyO8gn0kuZWuJnHRdM3221u+8o8fo/29fMDjOzCcFusHXAECLf8MvraeCC\n4PkFRH4/EiIdQJa9WUrkAxIAM6tJZFfLd/t6Y7Dv+R/Afe7+ZolZi4n0DPa237/kULpLgfpmVrtE\nQWhZjgyLgbbAjFLTf/QzlVjfuOA9W4GGpT4Eo/ECkd5EGjAzKBC4+3bgZuBmM2sNvEGk9/JEBdsp\nr3pmVrNEQWhJ5Heyisg3/wPdfZ/bswzPAw8BJ7j7FjO7jz0Xg7KGRn4WmGFmXYCORP6vSIjUM5C9\neQG42My6mlkWcAfwsbsvKMd7nwRmu/tdpaY/S+Sb80+DA7PZwcHKgrJW4u6LiexeuTNY9mDg0mA9\ne/M4cKuZtbeIg82sAZEP4f3N7LzgAO85QCdgrLsvA/4N/M3M6phZmpm1NbOflOPn3e1F4HhgKP/t\nFWBmR5lZZzNLB9YT2d0W9VlSu39nRL7IpQW/k331MG42sxpmdgTQH3jZI2doPQbca2aNgnU3N7Of\nljNKbSI9ti1m1gM4by/LrgD2KznB3ZcAnxLpEYx2983lbFfiRMVA9sjd3yayf3k0sIzIN+1zy/n2\nc4HT7cdnFB0RfLifClxP5EDwYiIHe/f2f3Egkf3QS4kc4L0pyLY39xDZRfVvIh++TxA5FrCayAfi\nNUR2ef0O6O/uq4L3/QyoAcwkcmD5FX68C2uvgoIyGegNvFRiVpNgXeuJ7Ep6j2DXSHAWz/ByNjGI\nyDf6YUTO0tlM5EN9T5YHP8dSIscwhrj77GDetcA8YIqZrQfeJnLAvTyuAG4xsw1EDpyP2suy9wNn\nmdn3ZvZAielPEzkIrV1E1YC56+Y2IsnIIqfsPuvuZfa6wmZmRxLp4bVyfRCFTj0DEalywa6tq4HH\nVQiqBxUDSVhm9map3VC7H9eHnU32zMw6AmuJ7H67L+Q4EtBuIhERUc9ARERUDEREhAS66Kxhw4be\nunXrsGOIiCSUqVOnrnL3/H0tlzDFoHXr1hQVFYUdQ0QkoZhZ6eFXyqTdRCIiomIgIiIqBiIigoqB\niIigYiAiIqgYiIgICXRqaUV9sXgtxRu2hh2jUjo2rU1BvdywY4hIEkv6YjBs4jeM+2r5vhesxmpn\nZTDht31pWCsr7CgikqSSvhhcf2JHrjq6XdgxKmzVxq1c9nQR94z/mjtO7xx2HBFJUklfDFo2SPzd\nKxf0bMXfJy/gZ71acUCTOmHHEZEkpAPICeCXx7andnYmt42dhYYcF5F4UDFIAHm5Nbj6mPZ8OG8V\n784uDjuOiCQhFYMEMahXK/ZrWJPb35jF9p27wo4jIklGxSBBZKanccNJHZm/8geenVKuQQhFRMpN\nxSCBHH1AIw5v15D73p7L2k3bwo4jIklExSCBmBl/6N+RDVu2c9/bc8OOIyJJRMUgwRzQpA7nHNqS\nZ6cs5JuVG8OOIyJJQsUgAV1z/P5kZ6Zzx+uzwo4iIkkirsXAzJ40s2Izm1HGvGvMzM2sYTwzJKOG\ntbK46uh2vDO7mA/nrgo7jogkgXj3DEYC/UpPNLMWwPHAoji3n7Qu7tOaFvVzuHXsTHboVFMRqaS4\nFgN3fx9YU8ase4HfAbqctoKyMtK57oSOzFmxgZeKFocdR0QSXJUfMzCzU4Hv3P2Lqm472ZxwUBN6\ntK7PPf/+mvVbtocdR0QSWJUWAzPLBa4H/ljO5QebWZGZFa1cuTK+4RKQmXFj/06s2bSNhyfMCzuO\niCSwqu4ZtAXaAF+Y2QKgAJhmZk3KWtjdR7h7obsX5ufnV2HMxNG5oC5ndCvgqQ8XsGj1prDjiEiC\nqtJi4O7T3b2Ru7d299bAEuAQd0/su8+E7Hf9OpCeZtz5pk41FZGKifeppS8Ak4EOZrbEzC6NZ3up\nqnGdbIb2bcubM5bz8fzVYccRkQQU77OJBrp7U3fPdPcCd3+i1PzW7q4T5WPg8iP2o2ndbG57fRa7\ndukkLRGJjq5AThI5NdK5tt8BTP9uHWM++y7sOCKSYFQMksgpXZrRpUUed781mx+27gg7jogkEBWD\nJJKWZvyxfydWrN/Ko+99E3YcEUkgKgZJpnurepzcpRkjPpjP0rWbw44jIglCxSAJXduvA+5w17jZ\nYUcRkQShYpCECurlctkRbfjH50v5bNH3YccRkQSgYpCkhvZtR37tLG4dOxN3nWoqInunYpCkamVl\n8NvjOzBt0Vr+9eWysOOISDWnYpDEzuxeQKemdbjlX1/x1dJ1YccRkWpMxSCJpacZ95/blcz0NM55\ndAofzNXIryJSNhWDJNe+cW1evaIPBfVyuPipT3ll6pKwI4lINaRikAKa1M1m1JBeHLZffX7z8hc8\n+M5cHVQWkR9RMUgRdbIzeeqiHpzRrTl/G/81142Zrnsni8h/ZIQdQKpOjYw0/nZ2F5rl5fDQhHms\nWL+Fh847hJpZ+m8gkurUM0gxZsZvftqBO07vzHtfr+TcEVMo3rAl7FgiEjIVgxR13mEteexnhcwr\n3sgZj0zim5Ubw44kIiFSMUhhx3RszIuDe7Jl+07OHDaJogVrwo4kIiFRMUhxXVrkMWZoH+rl1uC8\nxz/mzem6WlkkFakYCC0b5DJ6aG8OalaHK56fxhMffht2JBGpYnEtBmb2pJkVm9mMEtPuNrPZZval\nmb1qZnnxzCDlU79mDZ6/vCfHd2rMrWNncuvYmbqXskgKiXfPYCTQr9S08cBB7n4w8DVwXZwzSDll\nZ6bzyPnduah3a5748Ft+8cJnbNm+M+xYIlIF4loM3P19YE2paf9299036J0CFMQzg0QnPc246eRO\n/OGkjrw+fRmDnviYtZu2hR1LROIs7GMGlwBvhpxBSjEzLjtiPx4c2I0vFq/jzGGTWLxmU9ixRCSO\nQisGZnYDsAN4bi/LDDazIjMrWrlSI25WtZO7NOOZS3uwcsNWTn9kEtOXaBhskWQVSjEws4uA/sD5\nvpcR09x9hLsXunthfn5+leWT/zpsvwaMHtqbrIw0zhkxmQlzisOOJCJxUOXFwMz6Ab8DTnF37XtI\nAJFhsHvTpmFNLnu6iJc+XRR2JBGJsXifWvoCMBnoYGZLzOxS4CGgNjDezD43s+HxzCCx0ahONi/9\nvBd92jXk2tHTuWf81xoGWySJxHW4SncfWMbkJ+LZpsRPrawMnriwkOvHTOeBd+aybO1m7jijM5np\nYZ+HICKVpbGLJSqZ6WncddbBNMvL4f535rJiw1YeOf8QamkYbJGEpq90EjUz41fH7c9fzuzMR/NW\ncfbwyRSv1zDYIolMxUAq7JxDW/L4hYUsWP0Dpz8yiXnFG8KOJCIVpGIglXJUh0aM+nkvtu7YxRmP\nTOLj+avDjiQiFaBiIJV2UPO6vHpFbxrWzmLQE58w9sulYUcSkSipGEhMtKify5ihvenSoi5XPf8Z\nj38wX6eeiiQQFQOJmbzcGjxz6WGc2LkJt70+i5v/NZOdGgZbJCGoGEhMZWem89DAQ7j08DaMnLSA\nK5+bpmGwRRKAioHEXFqacWP/TtzYvxNvzVzO+Y9/zJofNAy2SHWmYiBxc+nhbXj4vEOY/t06zho2\niUWrNRSVSHWlYiBxdWLnpjx32WGs2bSNM4Z9xBeL14YdSUTKoGIgcXdo6/q8MqQ32ZnpnDtiCu/O\nXhF2JBEpRcVAqkS7RrUYc0Vv2jWqxWVPF/H8xxoGW6Q6UTGQKtOodjYvDu7Jkfvnc/2r0/nrW3N0\nLYJINaFiIFWqZlYGj/+skHMPbcFDE+Zxzagv2LZjV9ixRFKexh2WKpeRnsadZ3SmeV4Ofxv/NSs2\nbGHYBd2pk50ZdjSRlKWegYTCzPjFMe3564AufDx/DWcPn8zydRoGWyQsKgYSqrO6F/DUxYey5PvN\nnP7IR8xZrmGwRcKgYiChO6J9Pi/9vCc7dzlnDZ/EpG9WhR1JJOWoGEi1cGCzurx6ZR+a1Mnmwic/\n4Z+ffxd2JJGUss9iYGZXB//2iXblZvakmRWb2YwS0+qb2Xgzmxv8Wy/a9Upyap6XwytDenNIy3pc\n/eLnPDJxnk49Faki5ekZXBz8+2AF1j8S6Fdq2u+Bd9y9PfBO8FoEgLq5mfz90h6c3KUZd42bw43/\nnKFhsEWqQHlOLZ1lZnOBZmb2ZYnpBri7H7ynN7r7+2bWutTkU4G+wfOngYnAteXMKykgKyOd+8/p\nSrO8bB59bz7fb9rOw+cdEnYskaS2z2Lg7gPNrAnwFnBKDNps7O7LgufLgcZ7WtDMBgODAVq2bBmD\npiVRpKUZ153Qkaz0NB54dx5X9l1Pp2Z1wo4lkrTKdQDZ3Ze7exegGMh294W7H5Vp3CM7hPe4D8Dd\nR7h7obsX5ufnV6YpSVCXHN6GGulpvDx1cdhRRJJauc8mMrOTgc+BccHrrmb2WgXaXGFmTYN1NCVS\nYETKlJdbg+M6Neafny/VsBUicRTNqaV/AnoAawHc/XOgTQXafA24MHh+IfDPCqxDUshZ3QtY88M2\n3p2t7w0i8RJNMdju7utKTdvraR5m9gIwGehgZkvM7FLgz8BxwUHpY4PXInt0RPuGNKqdxSvaVSQS\nN9EMVPeVmZ0HpJtZe+D/gEl7e4O7D9zDrGOiaFdSXEZ6GmccUsBjH8yneMMWGtXODjuSSNKJpmfw\nC+BAYCvwArAe+GU8QomUNqCwgJ27nH98piuTReKh3MXA3Te5+w3AUUBfd7/B3TXMpFSJtvm1OKRl\nHi8XLdFVySJxEM3ZRJ3N7DPgKyK7jKaa2UHxiybyYwMKWzC3eCNfLil96EpEKiua3USPAr9291bu\n3gq4BhgRn1gi/+ukg5uSnalrDkTiIZpiUNPdJ+x+4e4TgZoxTySyB3WyM+l3YBNe+3wpW7bvDDuO\nSFKJphjMN7Mbzax18PgDMD9ewUTKMqCwBeu37ODfM1eEHUUkqURTDC4B8oExwGigYTBNpMr02q8B\nzfNyeLlIu4pEYqnc1xm4+/dEri0QCU1amnFm9wIefHcuy9ZtpmndnLAjiSSFaM4mGm9meSVe1zOz\nt+ITS2TPzjqkAHcYM03XHIjESjS7iRq6+9rdL4KeQqPYRxLZu5YNcjmsTX1eLlqsaw5EYiSaYrDL\nzP5zUwEza8U+xiYSiZcBhS1YsHoTRQu/DzuKSFKIphjcAHxoZs+Y2bPA+8B18Yklsncndm5CzRrp\nOpAsEiPRDEcxDjgEeAl4Eeju7v85ZmBmB8Y+nkjZcmtkcNLBTXn9y2Vs2rYj7DgiCS+angHuvsrd\nxwaPVaVmPxPDXCL7dFb3FvywbSdvTl8edhSRhBdVMdgHi+G6RPbp0Nb1aN0gV8NTiMRALIuBDiZL\nlTIzzupewJT5a1i0elPYcUQSWiyLgUiVO+OQAszglWlLwo4iktBiWQy2xXBdIuXSLC+Hw9s1ZPTU\nJezapc6pSEVFVQzM7BQz+2vwOLnkPHfvGdtoIuUzoLAF363dzJT5q8OOIpKwohmO4k7gamBm8Pg/\nM7ujog2b2a/M7Cszm2FmL5iZbmwrFXJ8p8bUzs7g5anaVSRSUdH0DE4CjnP3J939SaAf0L8ijZpZ\ncyKD3hW6+0FAOnBuRdYlkp2ZzildmvHmjGWs37I97DgiCSnaYwZ5JZ7XrWTbGUCOmWUAucDSSq5P\nUtiAwhZs2b6L179cFnYUkYQUTTG4E/jMzEaa2dPAVOD2ijTq7t8BfwUWAcuAde7+79LLmdlgMysy\ns6KVK1dWpClJEV0K6tK+US0NTyFSQeUqBmZmwIdAT/57c5te7v5SRRo1s3rAqUAboBlQ08wuKL2c\nu49w90J3L8zPz69IU5IizIwBhQVMW7SWecUbw44jknDKVQw8Mk7wG+6+zN1fCx6VGQPgWOBbd1/p\n7tuJFJjelVifCKd1a056mjFa1xyIRC2a3UTTzOzQGLW7COhpZrlBr+MYYFaM1i0pqlHtbPrun8+Y\naUvYqWsORKISTTE4DJhiZt+Y2ZdmNt3MvqxIo+7+MfAKMA2YHuQYUZF1iZQ0oLCAFeu38v5cHWMS\niUa574EM/DSWDbv7TcBNsVynyNEHNKZ+zRq8UrSEozroRnwi5RXN/QwWAi2Ao4Pnm6J5v0hVqJGR\nxqldmzF+5grWbtIIKSLlFc0VyDcB1/Lfu5tlAs/GI5RIZQzo3oJtO3fx2he6dEWkvKL5Zn86cArw\nA4C7LwVqxyOUSGV0alaHTk3r8HKRzioSKa9oisG24BRTBzCzmvGJJFJ5AwoLmP7dOmYvXx92FJGE\nEE0xGGVmjwJ5ZnY58DbweHxiiVTOqV2bk5lu6h2IlFM0B5D/SuR00NFAB+CP7v5AvIKJVEb9mjU4\ntmNj/vHZd2zfuSvsOCLVXjQHkG8EZrv7b939N+4+3swGxzGbSKUMKCxg9Q/bmDC7OOwoItVeNLuJ\nfgGMM7OjSkwbEuM8IjFzZPt88mtnMUq7ikT2KZpi8B1wAvBnM/ttMM1iH0kkNjLS0zi7sIC3Z63g\n4QnziJz/ICJlieqiMXdfBPwE6GRmLwM5cUklEiNXH7M/p3Vtxt1vzeH6V2ewQ8cPRMoUTTEoAnD3\nLe5+MTARqBGPUCKxUiMjjXvP6coVfdvywieLGPzMVDZt2xF2LJFqJ5qziS4v9fphd98v9pFEYsvM\n+F2/A7jttIOYOKeYc0dMYeWGrWHHEqlWojmbqI+ZjTezr81s/u5HPMOJxNIFPVsxYlAhc1ds5Ixh\nHzF/pW6CI7JbNLuJngDuAQ4HDi3xEEkYx3ZqzAuDe7Jp607OHDaJqQvXhB1JpFqIphisc/c33b3Y\n3VfvfsQtmUicdG2Rx5grelM3J5PzHvuYcTOWhR1JJHTRFIMJZna3mfUys0N2P+KWTCSOWjWoyeih\nvenUrA5Dn5vGUx99G3YkkVBFc3Obw4J/C0tMc+Do2MURqToNamXx/GU9ufrFz7j5XzNZunYz153Q\nkbQ0XT4jqafcxcDdj9r3UiKJJadGOsMu6M4t//qKxz74lqXrtvC3AV3IzkwPO5pIlSp3MTCzukRu\nU3lkMOk94BZ3XxePYCJVJT3N+NMpB9K8Xg53vDGbleu3MuJn3cnL1WU0kjqiOWbwJLABODt4rAee\nqmjDZpZnZq+Y2Wwzm2VmvSq6LpHKMjMGH9mWBwZ24/PFazlr+GSWfL8p7FgiVSaaYtDW3W9y9/nB\n42agMhed3Q+Mc/cDgC7ArEqsSyQmTunSjL9f2oPi9Vs4/ZFJzPhOHV9JDdEcQN5sZoe7+4cQuQgN\n2FyRRoNdTkcCFwG4+zZAdy+XaqHnfg14ZWhvLnryE855dDK3n96ZgnrhDcOVXzuLVg10Y0GJr2iK\nwVDg6eCDHOB74MIKttsGWAk8ZWZdgKnA1e7+QwXXJxJT+zeuzatX9uGipz7lly99HmqWNIObTz2I\nQT1bhZpDkpuVd1hfM8sCzgLaAnnAOsDd/ZaoGzUrBKYAfdz9YzO7H1jv7jeWWm4wMBigZcuW3Rcu\nXBhtUyKVsnnbTqYt+p5dIQ5/PfKjBbwzu5ihfdvy2+M76NRXiYqZTXX3wn0tF03P4J/AWmAakXsb\nVMYSYIm7fxy8fgX4femF3H0EMAKgsLBQg9FLlcupkU6fdg1DzdBrvwb88bWvGDbxG5au3cxdZx1M\nVoZOfZXYiqYYFLh7v1g06u7LzWyxmXVw9znAMcDMWKxbJNlkpKdx+2kH0Twvh7vfmkPx+q0MH9Sd\nujmZYUeTJBLN2USTzKxzDNv+BfCcmX0JdAXuiOG6RZKKmXHlUe2495wuFC1cw9nDJ7N0bYXO3xAp\nUzTHDGYC7YBvga1Ebnnp7n5w/OL9V2FhoRcVFVVFUyLV2kfzVjHkmanUzMrgqYsPpWPTOmFHkmqs\nvMcMoukZnAC0B44HTgb6B/+KSBXq064ho4ZErtEcMHwyH81bFXIiSQbR3OlsYVmPeIYTkbJ1bFqH\nV6/sTfO8HC588hPGTFsSdiRJcNH0DESkGmlaN4eXh/bi0Nb1+fWoL3h4wjzKu9tXpDQVA5EEVic7\nk6cv6cFpXZtx91tzuP7VGezYuSvsWJKAojm1VESqoRoZadx7Tlea5eXwyMRvWLF+Cw+d143cGvrz\nlvJTz0AkCZgZv+t3ALeddhAT5xRz7ogprNywNexYkkBUDESSyAU9WzFiUCFfr9jAGcM+Yv7KjWFH\nkgShYiCSZI7t1JgXB/di09adnDlsElMXrgk7kiQAFQORJNS1RR5jruhN3ZxMLnj8Exav0Y16ZO9U\nDESSVKsGNXn+8p44zl/GzQ47jlRzKgYiSaxZXg6Dj2zL2C+XaXeR7JWKgUiSG/KT/WhcJ4tbxs5i\n1y5dlCZlUzEQSXK5NTL47U8P4IvFa3nti6Vhx5FqSsVAJAWc0a05nZvX5S/jZrN5286w40g1pGIg\nkgLS0owb+3di2botjHh/fthxpBpSMRBJET3a1OfEzk0Y/t43LF+3Jew4Us2oGIikkN/368jOXc7d\nb80JO4pUMyoGIimkZYNcLjm8DaOnLWH6knVhx5FqRMVAJMVceVRbGtaqwS1jv9L9D+Q/Qi0GZpZu\nZp+Z2dgwc4ikktrZmfz6uA58uuB73pyxPOw4Uk2E3TO4GpgVcgaRlHPOoS04oElt7nxzFlu261RT\nCbEYmFkBcBLweFgZRFJVeprxh5M6sXjNZkZOWhB2HKkGwuwZ3Af8DtA9+kRCcHj7hhzbsREPvTtP\nN8KRcIqBmfUHit196j6WG2xmRWZWtHLlyipKJ5I6rj+xI1u27+Se8V+HHUVCFlbPoA9wipktAF4E\njjazZ0sv5O4j3L3Q3Qvz8/OrOqNI0tsvvxaDerXipU8XMXv5+rDjSIhCKQbufp27F7h7a+Bc4F13\nvyCMLCKp7upj2lM7O5Pbxs7SqaYpLOyziUQkZHm5NfjVse35cN4q3p1dHHYcCUnoxcDdJ7p7/7Bz\niKSy83u2om1+TW5/fRbbduicjlQUejEQkfBlpqdxw0kdmb/qB56dsjDsOBICFQMRAeCoDo04on1D\n7n9nLms3bQs7jlQxFQMRAcAsciHahi3bue/tuWHHkSqmYiAi/9GhSW0G9mjJM1MWMq94Y9hxpAqp\nGIjIj/z6uP3JzUznjjc0bFgqUTEQkR9pUCuLq45ux7uzi/lgrq78TxUqBiLyPy7q05qW9XO5bews\nduzUqaapQMVARP5HVkY61594AHNWbODFTxeHHUeqQEbYAUSkevrpgU04rE19bh07k+HvfVPh9TSo\nlcXdZx3M/o1rxzCdxJqKgYiUycy4+6wuDHtvHlsrcVXyB3NXceawSYwYVEivtg1imFBiyRJlYKrC\nwkIvKioKO4aIRGnJ95u46KlPWbR6E3cPOJhTuzYPO1JKMbOp7l64r+V0zEBE4qqgXi6jh/Sma8s8\nrn7xc4a/941GR62GVAxEJO7q5mby90t6cNLBTfnzm7O56bWv2LlLBaE60TEDEakS2ZnpPHhuN5rn\n5TDi/fksW7eFB87tRk6N9LCjCeoZiEgVSkszrj+xI386uRNvz1rBwMemsHqj7r9cHagYiEiVu6hP\nG4ad351Qa6FkAAAIXklEQVRZy9Zz5rBJLFj1Q9iRUp6KgYiEot9BTXj+8sNYt3k7ZwybxGeLvg87\nUkpTMRCR0HRvVZ/RQ3tTKyuDgY9NYfzMFWFHSlkqBiISqv3yazHmit50aFybnz9TxDO601ooVAxE\nJHQNa2XxwuCeHNWhETf+YwZ/GTebXTr1tEqFUgzMrIWZTTCzmWb2lZldHUYOEak+cmtk8Oig7px3\nWEuGTfyGX436nK07doYdK2WEdZ3BDuAad59mZrWBqWY23t1nhpRHRKqBjPQ0bj/tIJrn5XD3W3Mo\nXr+V4YO6UzcnM+xoSS+UnoG7L3P3acHzDcAsQAOWiAhmxpVHtePec7pQtHANZw+fzNK1m8OOlfRC\nP2ZgZq2BbsDHZcwbbGZFZla0cqXuuCSSSk7vVsDIi3uwdO1mznhkErOWrQ87UlILtRiYWS1gNPBL\nd/+fLe3uI9y90N0L8/Pzqz6giISqT7uGjBrSC4Czh0/mo3mrQk6UvEIrBmaWSaQQPOfuY8LKISLV\nW8emdXj1yt40y8vhwic/Ycy0JWFHSkphnU1kwBPALHe/J4wMIpI4mtbNYdSQXhzauj6/HvUFD0+Y\np2GwYyysnkEfYBBwtJl9HjxODCmLiCSAujmZPH1JD07r2oy735rDDf+YwY6dFb8Dm/xYKKeWuvuH\ngIXRtogkrhoZadxzdlea5uUwbOI3rFi3hQfP60ZuDY3GX1mhn00kIhKNtDTj2n4HcOtpBzFhTjED\nR0xh5QYNg11ZKgYikpAG9WzFo4MKmbNiA2cOm8T8lRvDjpTQVAxEJGEd16kxL1zek41bd3DmsElM\nXahhsCtKxUBEElq3lvUYM7Q3dXMyOe+xKYybsTzsSAlJxUBEEl7rhjUZPbQ3HZvWYehzU3l60oKw\nIyUcFQMRSQoNamXxwuU9OeaAxtz02lfc+cYsDYMdBRUDEUkaOTXSeXRQ98jB5ffnc/VLGga7vHRy\nrogklfQ045ZTD6R5vRz+/OZsvl21kf0b1Q47VqX8+vj9KaiXG9c2VAxEJOmYGUN+0pamdbN54J25\nfLpwTdiRKmXztvj3blQMRCRpndq1Oad21a1SykPHDERERMVARERUDEREBBUDERFBxUBERFAxEBER\nVAxERAQVAxERIUEuOjOzk4FVZraw1Ky6wLpyTGsIrIpTvH0pK09VrKO879nXcnubv6d51X27xGKb\nVHQ95XlPKm4T0N/KnqZVdpu0KtdS7l7tH8CI8k7fw7Si6pY93uso73v2tdze5ifqdonFNonndknF\nbRKr7ZKMfytVtU0SZTfRv6KYvqdlwxKLPBVZR3nfs6/l9jY/UbdLrLLEa7uk4jYB/a2Ut524sKDy\nJDUzK3L3wrBzyI9pu1Q/2ibVT1Vtk0TpGVTWiLADSJm0XaofbZPqp0q2SUr0DEREZO9SpWcgIiJ7\noWIgIiIqBiIikqLFwMz2M7MnzOyVsLNIhJmdZmaPmdlLZnZ82Hkkwsw6mtlwM3vFzIaGnUcizKym\nmRWZWf9YrTNpioGZPWlmxWY2o9T0fmY2x8zmmdnvAdx9vrtfGk7S1BHlNvmHu18ODAHOCSNvqohy\nu8xy9yHA2UCfMPKmgmi2SeBaYFQsMyRNMQBGAv1KTjCzdOBh4ASgEzDQzDpVfbSUNZLot8kfgvkS\nPyOJYruY2SnA68AbVRszpYyknNvEzI4DZgLFsQyQNMXA3d8H1pSa3AOYF/QEtgEvAqdWebgUFc02\nsYi/AG+6+7SqzppKov1bcffX3P0E4PyqTZo6otwmfYGewHnA5WYWk8/xhBiorhKaA4tLvF4CHGZm\nDYDbgW5mdp273xlKutRU5jYBfgEcC9Q1s3buPjyMcClsT38rfYEzgCzUM6hqZW4Td78KwMwuAla5\n+65YNJbsxaBM7r6ayL5pqSbc/QHggbBzyI+5+0RgYsgxpAzuPjKW60ua3UR78B3QosTrgmCahEfb\npHrSdql+qnSbJHsx+BRob2ZtzKwGcC7wWsiZUp22SfWk7VL9VOk2SZpiYGYvAJOBDma2xMwudfcd\nwFXAW8AsYJS7fxVmzlSibVI9abtUP9Vhm2igOhERSZ6egYiIVJyKgYiIqBiIiIiKgYiIoGIgIiKo\nGIiICCoGkgLMbGOJ5+PMbK2ZjQ0500QzKwyeXx9mFhFQMZDUczcwKOwQpagYSOhUDCSluPs7wIa9\nLRN8a7/fzD43sxlm1iOYXjO4CcknZvaZmZ0aTL/IzMYEvY65ZnZXiXUNC+5I9ZWZ3VxGW38GcoK2\nnjOzW8zslyXm325mV8fq5xfZk5QctVSkHHLdvauZHQk8CRwE3AC86+6XmFke8ImZvR0s3xXoBmwF\n5pjZg+6+GLjB3dcENyp5x8wOdvcvdzfi7r83s6vcvSuAmbUGxgD3BePUn0tkXHuRuFIxECnbCxC5\n6YiZ1Qk+/I8HTjGz3wTLZAMtg+fvuPs6ADObCbQiMhb92WY2mMjfWlMid6z6kj1w9wVmttrMugGN\ngc+CIddF4krFQKRspQftcsCAM919TskZZnYYkR7BbjuBDDNrA/wGONTdvzezkUQKyL48DlwENCHS\nKxGJOx0zECnbOQBmdjiwLvjW/xbwCzOzYF63fayjDvADsM7MGhO5l21ZtptZZonXrxK5H+6hQZsi\ncaeegaQUM/sAOACoZWZLgEvdvawP3C1m9hmQCVwSTLsVuA/4Mtif/y3Qf09tufsXwTpmE9ll9NEe\nFh0RrHOau5/v7tvMbAKw1t13VuDHFImahrAWKcXMJgK/cfeikNpPA6YBA9x9bhgZJPVoN5FINWJm\nnYB5RA5IqxBIlVHPQERE1DMQEREVAxERQcVARERQMRAREVQMREQEFQMREQH+H5T3flfRyjNLAAAA\nAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110b1fe50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot( l1_penalty_values, nonzero_coef, '-' )\n",
    "plt.xscale('log')\n",
    "plt.xlabel('l1 penalty')\n",
    "plt.ylabel('nonzero_coef')\n",
    "plt.title('nonzero_coef vs. l1 penalty')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Out of this large range, we want to find the two ends of our desired narrow range of `l1_penalty`.  At one end, we will have `l1_penalty` values that have too few non-zeros, and at the other end, we will have an `l1_penalty` that has too many non-zeros.  \n",
    "\n",
    "More formally, find:\n",
    "* The largest `l1_penalty` that has more non-zeros than `max_nonzeros` (if we pick a penalty smaller than this value, we will definitely have too many non-zero weights)\n",
    "    * Store this value in the variable `l1_penalty_min` (we will use it later)\n",
    "* The smallest `l1_penalty` that has fewer non-zeros than `max_nonzeros` (if we pick a penalty larger than this value, we will definitely have too few non-zero weights)\n",
    "    * Store this value in the variable `l1_penalty_max` (we will use it later)\n",
    "\n",
    "\n",
    "*Hint: there are many ways to do this, e.g.:*\n",
    "* Programmatically within the loop above\n",
    "* Creating a list with the number of non-zeros for each value of `l1_penalty` and inspecting it to find the appropriate boundaries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "l1_penalty_min = 127.42749857031335\n",
    "l1_penalty_max = 263.66508987303581"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***QUIZ QUESTION.*** What values did you find for `l1_penalty_min` and `l1_penalty_max`, respectively? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploring the narrow range of values to find the solution with the right number of non-zeros that has lowest RSS on the validation set \n",
    "\n",
    "We will now explore the narrow region of `l1_penalty` values we found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 127.42749857,  128.80363586,  130.17977314,  131.55591043,\n",
       "        132.93204771,  134.308185  ,  135.68432229,  137.06045957,\n",
       "        138.43659686,  139.81273414,  141.18887143,  142.56500872,\n",
       "        143.941146  ,  145.31728329,  146.69342057,  148.06955786,\n",
       "        149.44569514,  150.82183243,  152.19796972,  153.574107  ,\n",
       "        154.95024429,  156.32638157,  157.70251886,  159.07865615,\n",
       "        160.45479343,  161.83093072,  163.207068  ,  164.58320529,\n",
       "        165.95934258,  167.33547986,  168.71161715,  170.08775443,\n",
       "        171.46389172,  172.840029  ,  174.21616629,  175.59230358,\n",
       "        176.96844086,  178.34457815,  179.72071543,  181.09685272,\n",
       "        182.47299001,  183.84912729,  185.22526458,  186.60140186,\n",
       "        187.97753915,  189.35367644,  190.72981372,  192.10595101,\n",
       "        193.48208829,  194.85822558,  196.23436286,  197.61050015,\n",
       "        198.98663744,  200.36277472,  201.73891201,  203.11504929,\n",
       "        204.49118658,  205.86732387,  207.24346115,  208.61959844,\n",
       "        209.99573572,  211.37187301,  212.7480103 ,  214.12414758,\n",
       "        215.50028487,  216.87642215,  218.25255944,  219.62869672,\n",
       "        221.00483401,  222.3809713 ,  223.75710858,  225.13324587,\n",
       "        226.50938315,  227.88552044,  229.26165773,  230.63779501,\n",
       "        232.0139323 ,  233.39006958,  234.76620687,  236.14234416,\n",
       "        237.51848144,  238.89461873,  240.27075601,  241.6468933 ,\n",
       "        243.02303058,  244.39916787,  245.77530516,  247.15144244,\n",
       "        248.52757973,  249.90371701,  251.2798543 ,  252.65599159,\n",
       "        254.03212887,  255.40826616,  256.78440344,  258.16054073,\n",
       "        259.53667802,  260.9128153 ,  262.28895259,  263.66508987])"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l1_penalty_values = np.linspace(l1_penalty_min,l1_penalty_max,100)\n",
    "l1_penalty_values "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* For `l1_penalty` in `np.linspace(l1_penalty_min,l1_penalty_max,20)`:\n",
    "    * Fit a regression model with a given `l1_penalty` on TRAIN data. Specify `l1_penalty=l1_penalty` and `l2_penalty=0.` in the parameter list. When you call `linear_regression.create()` make sure you set `validation_set = None`\n",
    "    * Measure the RSS of the learned model on the VALIDATION set\n",
    "\n",
    "Find the model that the lowest RSS on the VALIDATION set and has sparsity *equal* to `max_nonzeros`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "RSS_list = []\n",
    "model_list =[]\n",
    "nonzero_coef =[]\n",
    "\n",
    "for l1_penalty in l1_penalty_values :\n",
    "    model = Lasso(alpha= l1_penalty, normalize=True) # set parameters\n",
    "    model.fit(X_train, y_train)\n",
    "    error = model.predict(X_valid) - y_valid\n",
    "    RSS_list.append(np.dot(error, error) )\n",
    "    model_list.append(model)\n",
    "    nonzero_coef.append( np.count_nonzero(model.coef_) + np.count_nonzero(model.intercept_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[10,\n",
       " 10,\n",
       " 10,\n",
       " 10,\n",
       " 10,\n",
       " 10,\n",
       " 9,\n",
       " 9,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 8,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 7,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6,\n",
       " 6]"
      ]
     },
     "execution_count": 91,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nonzero_coef"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x110f02850>"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8VFX6x/HPF5AiRUA6AUGq9BLQta36s8OCq9gWVGys\nddlVUXFdC+qusruKZVfFggURXStWbGBbRYP03jsEiKGHkOT5/XFv1jGmDCSTmSTP+/XKKzPnnnvv\nM5PyzLnn3HNkZjjnnHMHqlK8A3DOOVe2eSJxzjlXLJ5InHPOFYsnEuecc8XiicQ551yxeCJxzjlX\nLJ5IXFxIaiXJJFUJn38g6ZIDPFZLSTslVS7ZKMs+SVMlXRHvOPZHWYy5ovNE4g6IpA8ljcqnfKCk\njbkJIlpmdoaZPR/luVdKOjli39VmVsvMsvfnnFGeyyTtChPVOkkPRiYsSZ0lfSQpTVK6pOmSzozY\nfpukFeH+ayW9UtIxRktSU0mTJK0PX1ereMUSLUlDJX0V7zhc4TyRuAP1PDBEkvKUXwS8ZGZZcYgp\nVrqbWS3g18D5wGUR294BPgaaAI2APwDbAcIW1kXAyeH+ycCnpRh3XjnAh8A5cYzBlUOeSNyBegs4\nFDgut0BSPaA/8EL4vJ+kGZK2S1oj6a6CDhZ5OUNSG0mfSdoqaYuklyTVDbe9CLQE3gk/5d+cz2Wy\nZuEn7zRJSyVdGXGeuyS9KukFSTskzZOUHM0LNrOlwNdAj/BYDYDWwFNmlhl+fW1muZ+g+wCTzWxZ\nuP9GMxtbwOu/RdJrecoelvRI+HiopOVhzCskDY4m5jzxbzKzfwPfR1M/bPmNlDRf0o+SxkmqHrG9\nv6SZYUvsv5K65dn3JkmzJW2T9EruvpLqSXpX0ubwuO9KSsrn/EcATwC/Cn/W6ZL6SNqUp1V4tqRZ\n+/t+uJLjicQdEDPbA7wKXBxRfB6w0Mxy/6h3hdvrAv2AqyWdFcXhBfwNaAYcAbQA7grPexGwGvhN\neDlrdD77TwTWhvsPAv4q6aSI7QPCOnWBScBjUcSEpI4EiXNpWLQ1fDxe0lmSGufZ5VvgYkkjJCUX\n0YczEThTUu3wXJUJ3s8JkmoCjwBnmFlt4GhgZjQxl4DBwGlAG6A9cHsYX0/gWeD3BB8ongQmSaoW\nse95wOkEybYbMDQsrwSMAw4j+FCwh3x+Bma2ALgK+Cb8Wdc1s+8J3vdTI6peRPjhxcVHhUkkkp6V\nlCppbhR1j5f0g6QsSYPy2V4nvN4d1T+gcux5YFDEp9SLwzIAzGyqmc0xsxwzmw28THB5qFBmttTM\nPjazvWa2GXgwmv0AJLUAjgFuMbMMM5sJPM3PE95XZvZ+2KfyItC9iMP+IGkXsACYCvw7jNOAE4GV\nwD+BDZK+kNQu3D4euJ7gH/HnQKqkWwp4zauAH4DfhkUnAbvN7NvweQ7QRVINM9tgZvOieT9KwGNm\ntsbM0oD7gAvD8mHAk2Y2zcyyw/6tvcBREfs+Ymbrw33fIWzJmdlWM3vdzHab2Y7wuFH9fEPPA0MA\nJNUneH8nFOM1umKqMIkEeI7g01E0VhN8eirol/Me4Ivih1S2hZdwtgBnSWoD9CXiPZN0pKQp4SWM\nbQSfLhsUdVxJjSVNVNC5vR0YH81+oWZAWvgPKtcqoHnE840Rj3cD1VX44IBeQC2C/pEjgZq5G8xs\nrZldZ2ZtCD5h7yLi07GZvWRmJxO0fq4C7pF0WgHnmcBP/6h/Fz7HzHaF576KIFm9F7aOSsOaiMer\nCN5fCF7rjeHlpnRJ6QQtx2YR9fO+z7UAJB0s6UlJq8Kf7xdA3SJabJHGA78JW2rnAV+a2Yb9fmWu\nxFSYRGJmXwBpkWXhtfgPFYy0+TL3j9PMVoafoHPyHkdSb6Ax8FFpxF0GvEDwaX8IQX/ApohtEwgu\nHbUws0MIrnfn7ZzPz18BA7qaWZ3w2JH7FTZl9Xqgfu4lolBLYF0U5y2QBV4FvgHuKKDOGuBfQJd8\ntu0zs/8As/PbHvoPcELYX/BbIpKymU02s1OApsBC4KlivJz90SLicUuC9xeCBHNfeLkp9+tgM3s5\nimPeCHQAjgx/vseH5fn9bvziZ21m6wh+DmcTXNZ6MbqX4mKlwiSSAowFrjez3sBNhJcsCiKpEsEl\njJtKIbay4gXgZOBKIi5rhWoTtA4yJPUl+JQdjdrATmCbpObAiDzbNwGH57dj+M/8v8DfJFUPO4Av\nJ/gUWxLuB66U1CTsNL5bUltJlcLO98sI+kZyO8j7Saodbj8D6AxMKyD2zQSXzsYBK8I+gtwW2sDw\nE/hegvfmFx9yohFehsztx6gW2XlegGslJYWXkP4M5A5ffgq4Kmx1SlLN3NcaRRi1CfpF0sPj3llI\n3U1AkqSqecpfAG4GugJvRHFOF0MVNpFIqkXQafkfSTMJOgubFrHbNcD7ZrY21vGVFWa2kuAfd02C\n1keka4BRknYQfIp/NcrD3k1wOWkb8B6//EfxN+D28JJKfkn9QqAVwafnN4E7zeyTKM9dKDObQ3Ap\nZgSQGZ7nE4Ihv3MJ/tEPDatvB24juFSaDowGro4Y1ZWfCQSJOfKyaiXghvD1pBH0J1wNIOk4STv3\n4yXsIUhEELRs9hRRfwJB63s5sAy4F8DMUgg+PDwG/Egw6GBolDGMAWoQXBb9lmBIckE+A+YBGyVt\niSh/k+Dy2ptmtjvK87oYUUVa2ErBDVjvmlkXSXWARWZWYPKQ9FxY/7Xw+UsEo3ZyCK73VgX+bWa3\nxjh050qdpJXAFSWVhEuapGXA7xM1voqkwrZIzGw7sELSuQBh87zQ0TtmNtjMWppZK4LLWy94EnGu\n9Ek6h6D/5LN4x+IqUCKR9DJBB12HcOju5QRj5C8Pb2aaBwwM6/aRtBY4F3hSUmkNtXTOFUHSVOBx\n4FozO6C+IleyYn5pKxzSlwKsM7P+ebYNBf7OTyNqHjOzp8Nt2cCcsHy1mQ0Iy1sT3Lx1KDAduMjM\nMmP6IpxzzhWoNFokwwlu5CrIK2bWI/x6OqJ8T0T5gIjyB4CHzKwtQSff5TGI2TnnXJT2a4bW/RWO\nh+9HcOfqDSVwPBHc8Zs7jPR5gqkzHi9svwYNGlirVq2Ke3rnnKtQpk+fvsXMGhZVL6aJhGCY380E\n48YLco6k44HFwJ/C+wAguNs4BcgC7jez3EkC0yNmll3Lz+9YzlerVq1ISUk50NfgnHMVkqRV0dSL\n2aUtSf2BVDObXki1d4BWZtaNYCruyBvaDjOzZILWx5hwCo79Of8wSSmSUjZv3ry/4TvnnItSLPtI\njgEGhGPRJwInSfrZ3cXh5G17w6dPA70jtq0Lvy8nuNu3J8Gsn3Uj5kVKooCpL8xsrJklm1lyw4ZF\ntsycc84doJglEjMbaWZJ4T0XFwCfmdmQyDqSIm8GHEDYKR9OPVEtfNyAICnND2dbnUIwNTjAJcDb\nsXoNzjnnihbrPpJfULA8a4qZTQL+IGkAQT9IGj9NsXAEwf0bOQTJ7n4zmx9uuwWYKOleYAbwTGnG\n75xz7ucqxBQpycnJ5p3tzjm3fyRND/uqC1Vh7mx3zjkXG55InHPOFYsnEuecK4c2bsvg7nfmsS87\n9tOReSJxzrly5sslm+n3yJdM/G4NCzfsKHqHYir1UVvOOediIzvHePjTJTz62RLaNarFK4N70bZR\nNItWFo8nEuecKwdSt2cwfOJMvlm+lXN7JzFqYBdqVK1cKuf2ROKcc2XcV0u28MdXZrBzbxZ/H9SN\nc5NblOr5PZE451wZlZWdwyOfLuHRKUtp27AWE648ivaNY38pKy9PJM45VwZt3JbBHybO4LsVaQzq\nncSogZ05uGp8/qV7InHOuTJm6qJUbnh1Fhn7snnwvO6c3SsprvF4InHOuTJiX3YO/5i8iCe/WE7H\nJrV57He9aNuoVrzD8kTinHNlwZq03Vz/8gxmrklnyFEtub1fJ6ofVDqjsoriicQ55xLce7M3cOsb\ns8Hg34N7cWbXpkXvVIo8kTjnXILK2JfNqHfnM2Haanq0qMujF/akRf2D4x3WL3gicc65BLRo4w6u\nf/kHFm/ayVW/bsONp7bnoMqJOauVJxLnnEsgZsZL01Zzz7vzqV29Ci9c1pfj2yf2cuGeSJxzLkGk\n787k1tfn8OG8jRzXrgEPnteDhrWrxTusInkicc65BDBt+Vb++MpMNu/Yy8gzOnLlcYdTqZLiHVZU\nYn7BTVJlSTMkvZvPtqGSNkuaGX5dEZb3kPSNpHmSZks6P2Kf5yStiNinR6xfg3POxUpWdg4PfrSI\nC5/6luoHVeaNa47m979uU2aSCJROi2Q4sACoU8D2V8zsujxlu4GLzWyJpGbAdEmTzSw93D7CzF6L\nUbzOOVcq1qTtZvjEGfywOp1BvZO4e0BnalYrexeKYhqxpCSgH3AfcEO0+5nZ4ojH6yWlAg2B9IL3\ncs65suOtGeu4/a25CHjkwp4M6N4s3iEdsFhf2hoD3AwUttbjOeHlq9ck/WLuY0l9garAsoji+8J9\nHpKUb0+UpGGSUiSlbN68uTivwTnnSsz2jH38ceIM/vjKTDo0qc37w48r00kEYphIJPUHUs1seiHV\n3gFamVk34GPg+TzHaAq8CFxqZrnJaCTQEegD1Aduye/AZjbWzJLNLLlhw8QeOuecqximr0rjzIe/\n5J3ZG/jTye15ZdhRCXmD4f6K5aWtY4ABks4EqgN1JI03syG5Fcxsa0T9p4HRuU8k1QHeA/5sZt9G\n7LMhfLhX0jjgphi+BuecK7Z92Tk8+ukSHpuylOb1avDq739F78PqxTusEhOzRGJmIwlaD0g6Abgp\nMomE5U0jEsMAgk55JFUF3gReyNupnruPJAFnAXNj9Rqcc664Vm7ZxR9fmcnMNemc3as5dw/oTO3q\nB8U7rBJV6sMDJI0CUsxsEvAHSQOALCANGBpWOw84HjhUUm7ZUDObCbwkqSEgYCZwVSmG75xzUTEz\nXk1Zw93vzOegypV47Hc96d+tbPeFFERmFu8YYi45OdlSUlLiHYZzroJI25XJra/P5qP5mzi6zaH8\n87zuND2kRrzD2m+SpptZclH1yt6AZeecS2BTFqVy82uz2bZ7H7f3O4LLjmldpm4uPBCeSJxzrgTs\nzszir+8vYPy3q+nQuDbPX9qXTs0Kug+7fPFE4pxzxTRzTTo3vDKTFVt3ceVxrbnx1A4Js3phafBE\n4pxzB2hfdg6PfbaUx6YspXHtarx0xZEc3aZBvMMqdZ5InHPuACxN3cmNr85k1tptnN2zOXcO6Mwh\nNcrXsN5oeSJxzrn9kJNjPP/NSu7/YCE1qlZOyDXUS5snEueci9L69D2MeG0WXy/dykkdG3H/2V1p\nVKd6vMOKO08kzjlXBDPj9R/WcfekeWSbcf/ZXTm/TwuCCTacJxLnnCvE5h17GfnGHD5ZsIm+revz\nz3O7l4uJFkuSJxLnnCvAu7PX85e35rIrM7vC3Fx4IDyROOdcHmm7MvnL23N5b/YGureoyz/P7U7b\nRrXiHVbC8kTinHMRJs/byJ/fnMO2PfsYcVoHfn/84VSpHOs1AMs2TyTOOQek787krknzeGvmejo1\nrcOLlx/JEU0rxhQnxeWJxDlX4X08fxO3vTmHH3dl8qeT23PNiW04yFshUfNE4pyrsH7clcnd7wSt\nkCOa1mHc0D50aX5IvMMqczyROOcqpA/nbuT2t+aSvjuTP57cjmtOaEvVKt4KORCeSJxzFUrarkzu\nnDSPd2YFfSHPX9aHzs28FVIcMU+/kipLmiHp3Xy2DZW0WdLM8OuKiG2XSFoSfl0SUd5b0hxJSyU9\nIr+11DkXBTPj3dnrOeXBz/lw7gZuPKU9b193jCeRElAaLZLhwAKgoOEPr5jZdZEFkuoDdwLJgAHT\nJU0ysx+Bx4ErgWnA+8DpwAcxit05Vw6kbs/gL2/PZfK8TXRPOoQHBh1JxyY+IqukxDSRSEoC+gH3\nATfsx66nAR+bWVp4nI+B0yVNBeqY2bdh+QvAWXgicc7lw8x4bfpa7nl3Pnuzchh5RkcuP7a13xdS\nwmLdIhkD3AzULqTOOZKOBxYDfzKzNUBzYE1EnbVhWfPwcd7yX5A0DBgG0LJlywON3zlXRq1J281t\nb87hyyVb6NOqHg+c043DG/rd6bEQs7QsqT+QambTC6n2DtDKzLoBHwPPl9T5zWysmSWbWXLDhg1L\n6rDOuQSXnWOM+3oFp435gh9W/ciogZ15ZdivPInEUCxbJMcAAySdCVQH6kgab2ZDciuY2daI+k8D\no8PH64ATIrYlAVPD8qQ85etKPHLnXJm0ZNMObn59NjNWp/Pr9g2577ddSKrnM/XGWsxaJGY20syS\nzKwVcAHwWWQSAZAUuazYAIJOeYDJwKmS6kmqB5wKTDazDcB2SUeFo7UuBt6O1WtwzpUNe7OyGfPJ\nYs585EtWbtnFQ+d357lL+3gSKSWlfh+JpFFAiplNAv4gaQCQBaQBQwHMLE3SPcD34W6jcjvegWuA\n54AaBJ3s3tHuXAU2fVUat74+hyWpOxnYoxl/6d+JBrWqxTusCkVmFu8YYi45OdlSUlLiHYZzrgTt\nyNjH6A8XMX7aKpodUoN7f9uFEzs0indY5Yqk6WaWXFQ9v7PdOVfmfDRvI3e8PY9NOzIYenQrbjy1\nA7Wq+b+zePF33jlXZmzclsGdk4IbCzs2qc0TF/WmR4u68Q6rwvNE4pxLeNk5xvhvV/H3yYvIysnh\nltM7csVxrX2q9wThicQ5l9Dmrd/GbW/OZdaadI5r14D7zupKy0N9NFYi8UTinEtIuzOzGPPJEp75\nagX1Dj6Ihy/owYDuzfB5WhOPJxLnXML5eP4m7po0j3Xpe7iwb0tuPb0jhxx8ULzDcgXwROKcSxjr\n0vdw16R5fDx/Ex0a1+a1q35Fcqv68Q7LFcETiXMu7vZl5zDu6xU89PESDPPO9DLGE4lzLq5SVqbx\n5zfnsmjTDk4+ohF3/qYzLep7Z3pZ4onEORcXW3fu5f4PFvKf6WtpXrcGYy/qzamdm8Q7LHcAPJE4\n50pVdo4x8fvVjP5wEbv2ZnH1CW24/qS2HFzV/x2VVf6Tc86Vmllr0vnL23OZvXYbRx1en3sGdqFd\n48LWvXNlgScS51zM/bgrk9GTFzHx+9U0qFXN7wkpZzyROOdiJjvHeOX7NYyevJAdGVkMPboVfzql\nPXWq+z0h5UnUiUTSocDxwOoils91zjlmrknnjvAyVt/W9Rk1sDMdm9SJd1guBgpMJJLeBW41s7nh\nSoY/AClAG0ljzWxMaQXpnCs7tuzcy+gPF/Jqyloa1a7GmPN7MLCHX8YqzwprkbQ2s7nh40uBj83s\nYkm1ga8BTyTOuf/Zl53Di9+s4qFPFpOxL5vf//pwrj+pna8TUgEU9hPeF/H4/4CnAMxsh6ScaE8g\nqTJBS2admfUvoM45wGtAHzNLkTQYGBFRpRvQy8xmSpoKNAX2hNtONbPUaONxzpW8/y7dwl3vzGPx\npp0c164Bd/6mM20b1Yp3WK6UFJZI1ki6HlgL9AI+BJBUA9ifnrLhwAIg34ujYQtnODAtt8zMXgJe\nCrd3Bd4ys5kRuw02M18717k4W5O2m/veW8CH8zaSVK8GT17Um1M7NfbLWBVMYYnkcmAUcDJwvpml\nh+VHAeOiObikJKAfcB9wQwHV7gEe4OctkEgXAhOjOZ9zrnTszsziianLePKL5VSSuPGU9lx5/OFU\nP6hyvENzcVBgIgkvF10VWSapHjDVzKZEefwxwM1AvnccSeoFtDCz9yQVlEjOBwbmKRsnKRt4HbjX\nzCyfYw8DhgG0bNkyynCdc4UxMybNWs/f3l/Ixu0ZDOjejJFndqTpITXiHZqLowKn1pR0h6SO4eNq\nkqYAy4BNkk4u6sCS+gOpBQ0VllQJeBC4sZBjHAnsjuj0h+CyVlfguPDrovz2NbOxZpZsZskNGzYs\nKlznXBFmrUln0BPfMHziTBrWrsZ/rvoVj1zY05OIK/TS1vkEl50ALgm/NwTaA88DnxRx7GOAAZLO\nBKoDdSSNN7Mh4fbaQBdgang9tQkwSdKAiP6PC4CXIw9qZuvC7zskTQD6Ai8UEYtz7gBt3JbB6MkL\neeOHdTSoVY3R53RjUO8kKlXyfhAXKCyRZEZcMjoNmGhm2cACSUWO5zOzkcBIAEknADdFJBHMbBvQ\nIPd5OBrrptwkErZYziNodeTWqQLUNbMtkg4C+lN0QnPOHYA9mdk89eVynvh8GVnZxlW/bsO1J7ah\ntt+V7vIoLCHsldQF2AScCNwUse2AFwuQNApIMbNJRVQ9HlhjZssjyqoBk8MkUpkgiTx1oLE4534p\nJyfoB3ngw4Vs2JbBGV2aMPKMI2h5qK8R4vJXWCIZTnBvR0PgITNbARBeqpqxPycxs6nA1PDxHQXU\nOSGffY7KU7YL6L0/53bORe/7lWnc++58Zq3dRpfmdRhzfg+OPPzQeIflElxho7amAR3zKX8feD+W\nQTnnSteqrbt44MOFvD9nI43rVOMf53bn7J7NvR/ERaXQvo7wrvR6ZrYlfF4VGAr8ycyOiH14zrlY\nSt+dyaOfLeWFb1ZSpVIl/nRye648vrUvMuX2S2GTNl4APAnskrSE4KbCZ4HvgcGlE55zLhb2ZmXz\n4jerePSzpWzP2Md5vVtw46ntaVSnerxDc2VQYR87bgd6m9nS8MbBb4BBZvZO6YTmnCtpZsa7szcw\nevJC1qTt4bh2DbjtzCM4oqlP7+4OXFHDf5cCmNkPkpZ4EnGu7Jq2fCt//WAhs9ak07FJbZ6/rC+/\nbu8367riKyyRNJIUOT9W3cjnZvZg7MJyzpWUpak7uf+DBXyyIJUmdaozelA3zumVRGXvSHclpLBE\n8hQ/nyMr73PnXALbtnsfYz5dzIvfrKLGQZW5+fQOXHZMa59Y0ZW4wob/3l2agTjnSoaZ8cYP67jn\nvfls37OPC/q25MZT2nNorWrxDs2VUz7Gz7lyJG1XJn9+cw4fzN1In1b1GDWwi3eku5jzROJcOTFl\nUSo3vzab9N2Z3HpGR6487nDvB3GlwhOJc2Xc7sws/vr+AsZ/u5r2jWvx3KV96NzskHiH5SqQIhOJ\npGrAOUCryPpmNip2YTnnojFj9Y/c8OosVm7dxRXHtuam0zp4Z7orddG0SN4GtgHTgb2xDcc5F43M\nrBwe/nQxj09dRpM61XnpiiM5uk2Dond0LgaiSSRJZnZ6zCNxzkVlwYbt3PDqLBZs2M65vZP4y286\nUcfXCHFxFE0i+a+krmY2J+bROOcKlJWdwxOfL+PhT5dwSI2qPHVxMqd0ahzvsJyLKpEcCwyVtILg\n0pYAM7NuMY3MOfc/S1N3cOOrs5i1dhv9uzVl1MAu1K9ZNd5hOQdEl0jOiHkUzrl8ZecYT3+5nH9+\nvJiaVSvz2O960r9bs3iH5dzPRLP2+ipJ3flp7fQvzWxWbMNyzi1N3cmI12YxY3U6p3VuzL1ndaVh\nbb873SWeSkVVkDQceAloFH6Nl3R9tCeQVFnSDEnvFlLnHEkmKTl83krSHkkzw68nIur2ljRH0lJJ\nj0jyO65cuZKdYzz5+TLOfORLlm/excMX9OCJIb09ibiEFc2lrcuBI8P10pH0AMHaJI9GeY7hwAIg\n33kaJNUO60zLs2mZmfXIZ5fHgSvD+u8DpwMfRBmLcwltyaYd3PTabGatSefUTo2596wuvtiUS3hF\ntkgIOtezI55nh2VF7yglAf2Apwupdg/wAJARxfGaAnXM7FszM+AF4KxoYnEuke3LzuFfU5bS75Gv\nWL11F49c2JMnL+rtScSVCdG0SMYB0yS9GT4/C3gmyuOPAW6mgOnnw5UXW5jZe5JG5NncWtIMYDtw\nu5l9CTQH1kbUWRuW5XfsYcAwgJYtW0YZrnOlb/767dz8+izmrttOv65NuWtAZ7+M5cqUaDrbH5Q0\nlWAYMMClZjajqP0k9QdSzWy6pBPy2V4JeBAYms/uG4CWZrZVUm/gLUmdizpnnrjHAmMBkpOTbX/2\nda407M3K5rHPlvL41GXUPfggHh/cizO6No13WM7ttwITiaQ6ZrZdUn1gZfiVu62+maUVcexjgAGS\nzgSqA3UkjTezIeH22kAXYGrYX94EmCRpgJmlEE7HEiaiZUB7YB2QFHGOpLDMuTLlh9U/cstrs1mS\nupOzezXnL/06Uc/vC3FlVGEtkglAf4I5tiI/0St8fnhhBzazkcBIgLBFclNEEsHMtgH/mxwobPXc\nZGYpkhoCaWaWLelwoB2w3MzSJG2XdBRBZ/vFRN/p71zc7dqbxT8+WsRz/11J0zrVGXdpH07s0Cje\nYTlXLIWtkNg//N66JE8oaRSQYmaTCql2PDBK0j4gB7gqogV0DfAcUINgtJaP2HJlwueLN3PbG3NY\nl76Hi391GCNO60BtnyPLlQMKBj8VUkH61Mz+r6iyRJacnGwpKSnxDsNVUGm7Mrnn3fm8OWMdbRrW\n5IFzupHcqn68w3KuSJKmm1lyUfUK6yOpDhwMNJBUj5+G/NahgJFSzrmfmBlvzVzHqHfmsyMji+tP\nasu1J7b19UJcuVNYH8nvgT8CzQj6SXITyXbgsRjH5VyZtnrrbv781hy+XLKFni3rcv/Z3ejQJN9R\n8M6VeYX1kTwMPCzpejPzDm3norAvO4dnvlrBmE8WU6VSJUYN7MzgIw/ztdNduRbNfSSPSuoCdCIY\nxptb/kIsA3OurPlh9Y/c9sYcFm7cwWmdG3P3gC40OcTvTHflXzRrtt8JnECQSN4nmFb+K4LpSZyr\n8Lbt2cffJy/kpWmraVKnOk9e1JvTOjeJd1jOlZpopkgZBHQHZpjZpZIaA+NjG5Zzic/MeG/OBu5+\nZz5bd+5l6NGtuPHUDtSqFs2flXPlRzS/8XvMLEdSlqQ6QCrQIsZxOZfQ1qTt5o635zJl0Wa6NK/D\nM5ck0y2pbrzDci4uokkkKZLqAk8RjN7aSTCNvHMVTmZWDk99uZxHP1tCJYnb+x3B0KNbUaVyNBNp\nO1c+RdPZfk348AlJHxJM4z47tmE5l3j+u2wLf3lrLss27+L0zk244zedaFa3RrzDci7uCrshsVdh\n28zsh9g4KvIqAAAYJklEQVSE5FxiSd2ewX3vL+DtmetpUb8G44b24cSOPj+Wc7kKa5H8M/xeHUgG\nZhHclNgNSAF+FdvQnIuvrOwcXvhmFQ99vJi9WTn84aS2XON3pjv3C4XdkHgigKQ3gF5mNid83gW4\nq1Sicy5OvluRxh1vz2Xhxh0c374hdw/oTOsGNeMdlnMJKZrO9g65SQTAzOZKOiKGMTkXN6nbM7j/\ng4W8MWMdzevW4IkhvTmtc2PCNXOcc/mIJpHMlvQ0P907MhjwznZXruzLzuG5r1fy8KdLyMzK4boT\ngwkWa1T1y1jOFSWaRHIpcDUwPHz+BfB4zCJyrpR9uWQzd02ax7LNuzipYyPu6N+JVn4Zy7moRTP8\nNwN4KPxyrtxYk7abe9+bz+R5mzjs0IN5+uJkTu7UON5hOVfmFDb891UzO0/SHH6+1C4AZtYtppE5\nFyO7M7P495RljP1yOVUqiZtP78Dlx7amWhW/jOXcgSisRZJ7Kat/cU4gqTLBcOF1ucv35lPnHOA1\noE+4ZvspwP1AVSATGGFmn4V1pwJNgT3h7qeaWWpxYnQVg5nx9sz13P/BQjZuz+CsHs249YwjfIZe\n54qpsOG/G8Lvq4p5juHAAoKVFX9BUu2wzrSI4i3Ab8xsfTjceDI/X5VxsJn52rkuajPXpHP3O/OY\nsTqdrs0P4V+De9L7MF/u1rmSUNilrR3kc0mL4KZEM7N8E0OeYyQB/YD7gBsKqHYP8AAwIrfAzGZE\nbJ8H1JBUzcz2FnVO5yJt3JbB6MkLeeOHdTSsXY3Rg7oxqFcSlXyhKedKTGEtkpJYF3QMcDOQ77HC\naVhamNl7kkbkVwc4B/ghTxIZJykbeB2418x+kfAkDQOGAbRs2bIYL8GVRXsys3nyi2U8+flysnOM\nq09ow7UntvUp3p2Lgaj/qiQ14ucrJK4uon5/INXMpks6IZ/tlYAHgaGFHKMzQWvl1IjiwWa2Lrwk\n9jpwEfkssmVmY4GxAMnJyfm1rFw5lJNjvDVzHaM/XMTG7Rn069qUW8/oSIv6B8c7NOfKrWhWSBxA\nMO9WM4K1SA4j6PPoXMSuxwADJJ1JkIDqSBpvZkPC7bWBLsDU8K7hJsAkSQPCDvck4E3gYjNblntQ\nM1sXft8haQLQF1+t0QHTlm/l3vcWMGfdNrolHcKjv+tJn1beD+JcrEXTIrkHOAr4xMx6SjoRGFLE\nPpjZSGAkQNgiuSkiiWBm24AGuc/D0Vg3hUmkLvAecKuZfR1RpwpQ18y2SDqIYETZJ1G8BleOLd+8\nk/s/WMhH8zfR9JDqjDm/BwO6N/N+EOdKSTSJZJ+ZbZVUSVIlM5siacyBnlDSKCDFzCYVUu06oC1w\nh6Q7wrJTgV3A5DCJVCZIIk8daCyubNu6cy+PfLqEl6atplqVSow4rQOXHdPapzVxrpRFk0jSJdUi\nmBrlJUmpBP/Qo2ZmU4Gp4eM7CqhzQsTje4F7Czhc7/05tyt/9mRm8+zXK3h86jL27Mvmwr4tGP5/\n7WlYu1q8Q3OuQoomkQwEMoA/EUzYeAgwKpZBOZef7BzjtelreOjjJWzcnsEpnRpzy+kdaduoVrxD\nc65CK+w+kn8BEyL7KIDnYx+Scz9nZny2MJX7P1jIktSddG9Rl4cv6MGRhx8a79CccxTeIlkM/ENS\nU+BV4OU8Nwo6F3PTV6Vx/wcL+X7lj7RuUJN/D+7FGV2a+PogziWQwm5IfBh4WNJhwAXAs5JqAC8T\nJJXFpRSjq4CWbNrB6MmL+Hj+JhrUqsY9AztzQd+WHFS5UrxDc87lEc008qsIbgp8QFJP4FngDoJR\nU86VqLU/7uahj5fw5oy1HFy1Cjee0p7Ljm1NTb8j3bmEFc0NiVWAMwhaJf9HMPrqrphG5SqczTv2\n8q8pS5kwbTUILj+2NVef0Jb6NavGOzTnXBEK62w/BbgQOBP4DpgIDDOz/Rr661xhtu3ex9gvl/Hs\nVyvJzM5hUK8khp/cjmZ1a8Q7NOdclAprkYwEJgA3mtmPpRSPqyB27s3iua9XMPaL5WzPyOI33Zvx\np5PbcXhDH8rrXFlTWGf7SaUZiKsY9mRmM/7bVTz++TLSdmVy8hGNueGU9nRqVuSqBM65BOU9mK5U\nZOzL5uXvVvPvqcvYvGMvx7VrwI2ndqBHi7rxDs05V0yeSFxM7c3K5pXv1/DvKcvYuD2Dow6vz79+\n14u+rX1WXufKC08kLib2ZmXzaspa/j1lKRu2ZdCnVT0ePK87R7dtUPTOzrkyxROJK1F5E0ivlnX5\n+6DuHNP2UL8b3blyyhOJKxEZ+4JLWI9PDS5h9T6sHqMHdePYtg08gThXznkiccWyOzOLCdNWM/aL\n5aTu2EvfVvX5x7neAnGuIvFE4g7Ijox9vPDNKp75agVpuzL51eGH8siFPTnKZ+R1rsLxROL2y4+7\nMhn39Qqe++9KtmdkcUKHhlx3YluSfW105yqsmCcSSZWBFGCdmfUvoM45wGtAHzNLCctGApcD2cAf\nzGxyWH468DDBpJFPm9n9sX4NDjZtz+CpL5Yz4bvV7M7M5vTOTbjmxDZ0S/L7QJyr6EqjRTIcWADk\ne+uypNphnWkRZZ0IJonsDDQDPpHUPtz8L+AUYC3wvaRJZjY/duFXbKu37uaJL5bxWspass0Y0L0Z\n15zQhnaNa8c7NOdcgohpIpGUBPQD7gNuKKDaPQTT1I+IKBsITDSzvcAKSUuBvuG2pWa2PDz+xLCu\nJ5ISNn/9dp74fBnvzl5PlUqVGJScxFXHt6HloQfHOzTnXIKJdYtkDHAzkO/HV0m9gBZm9p6kyETS\nHPg24vnasAxgTZ7yIws49jBgGEDLli0PKPiKxsz4ZtlWnvxiOZ8v3kzNqpW5/NjWXHHc4TSuUz3e\n4TnnElTMEomk/kCqmU2XdEI+2ysBDwJDY3F+MxsLjAVITk62WJyjvMjJMd6ZvZ6nvlzO3HXbaVCr\nKiNO68CQIw/jkIMPind4zrkEF8sWyTHAAElnAtWBOpLGm9mQcHttoAswNbzfoAkwSdIAYB3QIuJY\nSWEZhZS7A7A0dQe3vD6H6at+pE3Dmtx/dlfO6tmc6gf5ApjOuejELJGY2UiCNU0IWyQ3RSQRzGwb\n8L+JlyRNDeukSNoDTJD0IEFnezuCxbUEtJPUmiCBXAD8LlavoTzbl53D2C+W8/AnSzi4WmX+eW53\nftuzOZUq+U2Ezrn9U+r3kUgaBaSY2aSC6pjZPEmvEnSiZwHXmll2uP91wGSC4b/Pmtm8Ugi7XJm9\nNp1bXp/Dgg3b6de1KXcN6EzD2tXiHZZzroySWfnvPkhOTraUlJR4hxF3ezKzefDjRTzz1Qoa1KrG\nqIFdOL1Lk3iH5ZxLUJKmm1lyUfX8zvYK4vPFm7n9rTmsSdvDhX1bcusZHTmkhnekO+eKzxNJObdl\n517ufXc+b81cz+ENazJx2FE+H5ZzrkR5IimnzIz/pKzlrx8sYNfeLP7wf+245oQ2PhrLOVfiPJGU\nQ0tTd3Dbm3P5bkUafVrV46+/7epTmjjnYsYTSTmyJzObx6YsYewXyzm4ahUeOKcr5/Zu4UN6nXMx\n5YmknJiyKJU73p7LmrQ9nN2rObedeQQNavmQXudc7HkiKeM2bNvDqHfm88HcjbRpWJOXrzyKX7Xx\nznTnXOnxRFJG7cvO4fn/ruShjxeTlWOMOK0DVx53OFWrVIp3aM65CsYTSRn0zbKt3DlpLos37eSk\njo24e0BnWtT36d2dc/HhiaQM2bgtg799sIC3Z64nqV4Nnro4mZOPaEQ46aVzzsWFJ5IyIDMrh2e/\nXsEjny4hK8e4/qS2XHNCW2pU9XtCnHPx54kkwU1ZmMo9785n+ZZdnHxEY+7o38lXKXTOJRRPJAlq\nxZZd3PPufD5bmMrhDWoy7tI+nNihUbzDcs65X/BEkmC2Z+zjsc+WMu7rFVSrUpnbzuzI0KNb+2gs\n51zC8kSSILJzjFdT1vDPjxaxdVcmg3olMeL0DjSq7WulO+cSmyeSBPDfpVu4570FLNiwneTD6jFu\naF+6Jh0S77Cccy4qnkjiaNuefdz82iwmz9tE87o1eOx3PenXtakP53XOlSkxTySSKgMpwDoz659n\n21XAtUA2sBMYZmbzJQ0GRkRU7Qb0MrOZ4druTYE94bZTzSw1xi+jxG3esZeLn/2Opak7GHFaBy4/\ntrVP8e6cK5NKo0UyHFgA1Mln2wQzewJA0gDgQeB0M3sJeCks7wq8ZWYzI/YbbGZldu3ctT/uZsjT\n09i0fS9PX9KHX7dvGO+QnHPugMV0KJCkJKAf8HR+281se8TTmkB+C8hfCEws+ejiY+WWXQx6/BvS\ndmUy/oojPYk458q8WLdIxgA3AwWuqiTpWuAGoCpwUj5VzgcG5ikbJykbeB2418zyS0AJZ336HgY/\nPY29Wdm88vtfcUTT/BppzjlXtsSsRSKpP5BqZtMLq2dm/zKzNsAtwO15jnEksNvM5kYUDzazrsBx\n4ddFBZx/mKQUSSmbN28uzkspEZt37GXI09PYvmcfL15+pCcR51y5EctLW8cAAyStJLg0dZKk8YXU\nnwiclafsAuDlyAIzWxd+3wFMAPrmdzAzG2tmyWaW3LBhfC8fpe/O5KJnprFhWwbPXdaHLs19aK9z\nrvyIWSIxs5FmlmRmrQgSwmdmNiSyjqR2EU/7AUsitlUCziOif0RSFUkNwscHAf2ByNZKwtm5N4tL\nxn3P8s27eOriZHofVj/eITnnXIkq9ftIJI0CUsxsEnCdpJOBfcCPwCURVY8H1pjZ8oiyasDkMIlU\nBj4BniqdyPdfxr5srnj+e+au28bjg3txbLsG8Q7JOedKnMpIP3WxJCcnW0pK6Y4WzszK4fcvpjB1\n8WbGnN+DgT2al+r5nXOuuCRNN7Pkour5TIAxkJNj3PSfWUxZtJn7zurqScQ5V655IilhZsY9781n\n0qz13Hx6B353ZMt4h+ScczHliaSE/XvqMsZ9vZJLj2nF1b9uE+9wnHMu5jyRlKBXvl/N3ycvYmCP\nZvylXyeffNE5VyF4IikhH87dwMg35nB8+4b8fVB3KlXyJOKcqxg8kZSA/y7dwh9enkn3FnV5Ykgv\nX83QOVeh+H+8Ypq1Jp0rX0ihdYOajBvah4Or+hIvzrmKxRNJMSzZtIOh476jXs2qvHB5X+oeXDXe\nITnnXKnzRHKA1qTt5qJnvqNK5Uq8dMWRNK7ja6s75yomTyQHIHVHBkOemcaefdm8eHlfDju0ZrxD\ncs65uPFEsp/Sd2dy0dPfkbp9L+Mu7UPHJj4dvHOuYvNEsh9yZ/JdsWUXT1+STK+W9eIdknPOxZ0P\nMYpS5Ey+TwzpzTFtfSZf55wDb5FEJTMrh6vHT2faijQePK87p3RqHO+QnHMuYXgiKUJWdg7DJ87w\nmXydc64AnkgKkZNjjHhtNh/M3chf+nfymXydcy4fnkgKYGb8+a25vDljHTed2p7Lj20d75Cccy4h\neSIpgCTaNqrFtSe24bqT2hW9g3POVVAxH7UlqTKQAqwzs/55tl0FXAtkAzuBYWY2X1IrYAGwKKz6\nrZldFe7TG3gOqAG8Dwy3GK0X7K0Q55wrWmm0SIYTJIX8TDCzrmbWAxgNPBixbZmZ9Qi/rooofxy4\nEmgXfp0ei6Cdc85FJ6aJRFIS0A94Or/tZrY94mlNoNCWhaSmQB0z+zZshbwAnFVC4TrnnDsAsW6R\njAFuBnIKqiDpWknLCFokf4jY1FrSDEmfSzouLGsOrI2oszYsy++4wySlSErZvHlzsV6Ec865gsUs\nkUjqD6Sa2fTC6pnZv8ysDXALcHtYvAFoaWY9gRuACZL2a1IrMxtrZslmltywYcMDeAXOOeeiEcsW\nyTHAAEkrgYnASZLGF1J/IuFlKjPba2Zbw8fTgWVAe2AdkBSxT1JY5pxzLk5ilkjMbKSZJZlZK+AC\n4DMzGxJZR1LkuNp+wJKwvGE42gtJhxN0qi83sw3AdklHSRJwMfB2rF6Dc865opX6pI2SRgEpZjYJ\nuE7SycA+4EfgkrDa8cAoSfsI+leuMrO0cNs1/DT894PwyznnXJwoRrdgJJTk5GRLSUmJdxjOOVem\nSJpuZslF1qsIiUTSZmBVvOOIUgNgS7yD2E9lLeayFi94zKWlrMUc63gPM7MiRytViERSlkhKieYT\nQCIpazGXtXjBYy4tZS3mRInX59pyzjlXLJ5InHPOFYsnksQzNt4BHICyFnNZixc85tJS1mJOiHi9\nj8Q551yxeIvEOedcsXgicc45VyyeSEqZpGclpUqam8+2GyWZpAbhc0l6RNJSSbMl9UqUeCVdL2mh\npHmSRkeUjwzjXSTptNKON4zhFzFL6iHpW0kzw1mh+4blifAet5A0RdL88P0cHpbXl/SxpCXh93pl\nIOa/h78XsyW9KaluxD5x/d0oKOaI7Qn191dYvAn392dm/lWKXwTTv/QC5uYpbwFMJrhxskFYdibB\nFDACjgKmJUK8wInAJ0C18Hmj8HsnYBZQDWhNMNlm5QSJ+SPgjIj3dWoCvcdNgV7h49rA4vC9HA3c\nGpbfCjxQBmI+FagSlj8QEXPcfzcKijl8nnB/f4W8xwn39+ctklJmZl8Aaflseohg7ZbI0Q8DgRcs\n8C1QV8HiXqWmgHivBu43s71hndSwfCAw0YLZm1cAS4G+pRZsqICYDchdiuAQYH34OBHe4w1m9kP4\neAfBiqLNw9ieD6s9z0+LuCVszGb2kZllhdW+5afZuuP+u1HI+wwJ+PdXSLwJ9/fniSQBSBpIsKb9\nrDybmgNrIp4XuJBXKWsPHCdpmoKFx/qE5YkaL8Afgb9LWgP8AxgZlidUzJJaAT2BaUBjC2a8BtgI\nNA4fJ3LMkS7jp0lVEzbmsvD3l+c9Tri/v1Kf/df9nKSDgdsILgmUFVWA+gTN/T7Aqwqm+09kVwN/\nMrPXJZ0HPAOcHOeYfkZSLeB14I9mtl3S/7aZmUlKuLH6eWOOKP8zkAW8FK/YChIZM0GMCf33l8/v\nRcL9/XmLJP7aEFzPnKVgEbAk4AdJTQgW7WoRUTdRFvJaC7wRNvm/I5jqvwGJGy8ESxS8ET7+Dz81\n+RMiZkkHEfyzeMnMcuPclHspJfyeewkjkWNG0lCgPzDYwov3JG7MCf33V8B7nHB/f55I4szM5phZ\nIzNrZcEiYGsJOtg2ApOAi8PRI0cB2yIudcTTWwQdfkhqD1QlmIF0EnCBpGqSWhMsSPZd3KL8ufXA\nr8PHJxEuokYCvMcKmh7PAAvM7MGITZP4aY2eS/hpEbeEjVnS6QR9DQPMbHfELnH/3cgv5kT++yvk\n9yLx/v5Ko0ffv342EuNlgjXp9xH80l6eZ/tKfho1IuBfBKMv5gDJiRBv+Is7HpgL/ACcFFH/z2G8\niwhHSSVIzMcC0wlGtUwDeifQe3wsQSfvbGBm+HUmcCjwKUHS+wSoXwZiXkpwnT637IlE+d0oKOY8\ndRLm76+Q9zjh/v58ihTnnHPF4pe2nHPOFYsnEuecc8XiicQ551yxeCJxzjlXLJ5InHPOFYsnEucK\nIWlnxOMPJaVLejfOMU2VlBw+vi2esTgHnkic2x9/By6KdxB5eCJxceeJxLkomdmnwI7C6oSthYcV\nrHsyVz+te1JTwTop30maEU4UiKShkt4IWztL8qwt8biCtVPmSbo7n3PdD9QIz/WSpFGS/hix/T7l\nWXPDuVjwSRudK3kHm1kPSccDzwJdCO44/szMLlOw2NN3kj4J6/cgmNl1L7BI0qNmtgb4s5mlSaoM\nfCqpm5nNzj2Jmd0q6Toz6wH/myH2DWCMpErABcRhGn9X8Xgica7kvQzBuiiS6oSJ41RggKSbwjrV\ngZbh40/NbBuApPnAYQTTjJwnaRjB32lTgoWLZlMAM1spaaukngRTzs8ws60l//Kc+zlPJM6VvLzz\nDhnBvE3nmNmiyA2SjiRoieTKBqqEk+7dBPQxsx8lPUeQfIryNDAUaELQGnIu5ryPxLmSdz6ApGMJ\nZozdRrCM6/XhjK6ErYbC1AF2AdskNQbOKKDevnCq8VxvAqcTrFMx+cBfgnPR8xaJc1GS9CXQEagl\nKXfm5vz+WWdImgEcRLBKIMA9wBhgdth/sYJgzY58mdms8BgLCS5zfV1A1bHhMX8ws8FmlilpCpBu\nZtkH8DKd228++69zJUjSVOAmM0uJ0/krEUwtfq6ZLSmqvnMlwS9tOVdOSOpEsB7Ip55EXGnyFolz\nzrli8RaJc865YvFE4pxzrlg8kTjnnCsWTyTOOeeKxROJc865Yvl/4nmVEAJO29kAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110a08890>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot( l1_penalty_values, RSS_list, '-' )\n",
    "plt.xlabel('l1 penalty')\n",
    "plt.ylabel('Validation RSS')\n",
    "plt.title('Valiation RSS vs. l1 penalty')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***QUIZ QUESTIONS***\n",
    "1. What value of `l1_penalty` in our narrow range has the lowest RSS on the VALIDATION set and has sparsity *equal* to `max_nonzeros`?\n",
    "2. What features in this model have non-zero coefficients?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nonzero_coef[20]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('bedrooms', -0.0)\n",
      "('bedrooms_square', -0.0)\n",
      "('bathrooms', 10854.443479272539)\n",
      "('sqft_living', 163.30354659048507)\n",
      "('sqft_living_sqrt', 0.0)\n",
      "('sqft_lot', -0.0)\n",
      "('sqft_lot_sqrt', -0.0)\n",
      "('floors', 0.0)\n",
      "('floors_square', 0.0)\n",
      "('waterfront', 507375.04119386961)\n",
      "('view', 41979.794839556045)\n",
      "('condition', 0.0)\n",
      "('grade', 116374.65694207841)\n",
      "('sqft_above', 0.0)\n",
      "('sqft_basement', 0.0)\n",
      "('yr_built', -2620.6721389082932)\n",
      "('yr_renovated', 0.0)\n"
     ]
    }
   ],
   "source": [
    "for x in zip(all_features, model_list[20].coef_  ):\n",
    "    print x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}