week-4-ridge-regression-assignment-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regression Week 4: Ridge Regression (gradient descent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook, you will implement ridge regression via gradient descent. You will:\n",
    "* Convert an SFrame into a Numpy array\n",
    "* Write a Numpy function to compute the derivative of the regression weights with respect to a single feature\n",
    "* Write gradient descent function to compute the regression weights given an initial weight vector, step size, tolerance, and L2 penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fire up graphlab create"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make sure you have the latest version of GraphLab Create (>= 1.7)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#import graphlab\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import sklearn\n",
    "from sklearn.linear_model import LinearRegression, Ridge"
   ]
  },
  {
   "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": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "dtype_dict = {'bathrooms':float, 'waterfront':int, 'sqft_above':int, 'sqft_living15':float, \n",
    "              'grade':int, 'yr_renovated':int, 'price':float, 'bedrooms':float, 'zipcode':str, \n",
    "              'long':float, 'sqft_lot15':float, 'sqft_living':float, 'floors':float, 'condition':int, \n",
    "              'lat':float, 'date':str, 'sqft_basement':int, 'yr_built':int, 'id':str, 'sqft_lot':int, 'view':int}\n",
    "\n",
    "sales = pd.read_csv('kc_house_data.csv', dtype=dtype_dict)\n",
    "sales = sales.sort_values(['sqft_living','price'])\n",
    "\n",
    "train_data = pd.read_csv('kc_house_train_data.csv', dtype=dtype_dict)\n",
    "test_data  = pd.read_csv('kc_house_test_data.csv', dtype=dtype_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we want to do any \"feature engineering\" like creating new features or adjusting existing ones we should do this directly using the SFrames as seen in the first notebook of Week 2. For this notebook, however, we will work with the existing features."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import useful functions from previous notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in Week 2, we convert the SFrame into a 2D Numpy array. Copy and paste `get_numpy_data()` from the second notebook of Week 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_numpy_data(data_sframe, features, output):\n",
    "    n = len(output)\n",
    "    data_sframe['constant'] = 1 \n",
    "    features = ['constant'] + features  \n",
    "    feature_matrix = data_sframe[features]\n",
    "    output_array = data_sframe[output]\n",
    "    return(feature_matrix, output_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, copy and paste the `predict_output()` function to compute the predictions for an entire matrix of features given the matrix and the weights:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def predict_output(feature_matrix, weights):\n",
    "    predictions = np.dot(feature_matrix, weights)\n",
    "    return(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing the Derivative"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now going to move to computing the derivative of the regression cost function. Recall that the cost function is the sum over the data points of the squared difference between an observed output and a predicted output, plus the L2 penalty term.\n",
    "```\n",
    "Cost(w)\n",
    "= SUM[ (prediction - output)^2 ]\n",
    "+ l2_penalty*(w[0]^2 + w[1]^2 + ... + w[k]^2).\n",
    "```\n",
    "\n",
    "Since the derivative of a sum is the sum of the derivatives, we can take the derivative of the first part (the RSS) as we did in the notebook for the unregularized case in Week 2 and add the derivative of the regularization part.  As we saw, the derivative of the RSS with respect to `w[i]` can be written as: \n",
    "```\n",
    "2*SUM[ error*[feature_i] ].\n",
    "```\n",
    "The derivative of the regularization term with respect to `w[i]` is:\n",
    "```\n",
    "2*l2_penalty*w[i].\n",
    "```\n",
    "Summing both, we get\n",
    "```\n",
    "2*SUM[ error*[feature_i] ] + 2*l2_penalty*w[i].\n",
    "```\n",
    "That is, the derivative for the weight for feature i is the sum (over data points) of 2 times the product of the error and the feature itself, plus `2*l2_penalty*w[i]`. \n",
    "\n",
    "**We will not regularize the constant.**  Thus, in the case of the constant, the derivative is just twice the sum of the errors (without the `2*l2_penalty*w[0]` term).\n",
    "\n",
    "Recall that twice the sum of the product of two vectors is just twice the dot product of the two vectors. Therefore the derivative for the weight for feature_i is just two times the dot product between the values of feature_i and the current errors, plus `2*l2_penalty*w[i]`.\n",
    "\n",
    "With this in mind complete the following derivative function which computes the derivative of the weight given the value of the feature (over all data points) and the errors (over all data points).  To decide when to we are dealing with the constant (so we don't regularize it) we added the extra parameter to the call `feature_is_constant` which you should set to `True` when computing the derivative of the constant and `False` otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def feature_derivative_ridge(errors, feature, weight, l2_penalty, feature_is_constant):\n",
    "    # If feature_is_constant is True, derivative is twice the dot product of errors and feature\n",
    "    if feature_is_constant == True:\n",
    "        derivative = 2 * np.dot(errors, feature)\n",
    "    # Otherwise, derivative is twice the dot product plus 2*l2_penalty*weight\n",
    "    else:\n",
    "        derivative = 2 * np.dot(errors, feature) + 2*l2_penalty*weight\n",
    "    return derivative\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To test your feature derivartive run the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5.6554166816e+13\n",
      "-5.6554166816e+13\n",
      "\n",
      "-22446749330.0\n",
      "-22446749330.0\n"
     ]
    }
   ],
   "source": [
    "(example_features, example_output) = get_numpy_data(sales, ['sqft_living'], 'price') \n",
    "my_weights = np.array([1., 10.])\n",
    "test_predictions = predict_output(example_features, my_weights) \n",
    "errors = test_predictions - example_output # prediction errors\n",
    "\n",
    "# next two lines should print the same values\n",
    "print feature_derivative_ridge(errors, example_features['sqft_living'], my_weights[1], 1, False)\n",
    "print np.sum(errors*example_features['sqft_living'])*2+20.\n",
    "print ''\n",
    "\n",
    "# next two lines should print the same values\n",
    "print feature_derivative_ridge(errors, example_features['constant'], my_weights[0], 1, True)\n",
    "print np.sum(errors)*2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gradient Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will write a function that performs a gradient descent. The basic premise is simple. Given a starting point we update the current weights by moving in the negative gradient direction. Recall that the gradient is the direction of *increase* and therefore the negative gradient is the direction of *decrease* and we're trying to *minimize* a cost function. \n",
    "\n",
    "The amount by which we move in the negative gradient *direction*  is called the 'step size'. We stop when we are 'sufficiently close' to the optimum. Unlike in Week 2, this time we will set a **maximum number of iterations** and take gradient steps until we reach this maximum number. If no maximum number is supplied, the maximum should be set 100 by default. (Use default parameter values in Python.)\n",
    "\n",
    "With this in mind, complete the following gradient descent function below using your derivative function above. For each step in the gradient descent, we update the weight for each feature before computing our stopping criteria.\n",
    "\n",
    "Loop over each feature:\n",
    "\n",
    "```python\n",
    "def ridge_regression_gradient_descent(feature_matrix, output, initial_weights, step_size, \n",
    "                                      l2_penalty, max_iterations=100):\n",
    "    print 'Starting gradient descent with l2_penalty = ' + str(l2_penalty)\n",
    "    \n",
    "    weights = np.array(initial_weights) # make sure it's a numpy array\n",
    "    iteration = 0 # iteration counter\n",
    "    print_frequency = 1  # for adjusting frequency of debugging output\n",
    "    \n",
    "    #while not reached maximum number of iterations:\n",
    "    while (iteration <= max_iterations):\n",
    "        \n",
    "        # compute the predictions based on feature_matrix and weights using your predict_output() function\n",
    "        predictions = np.dot(feature_matrix, weights)\n",
    "         errors = predictions - output\n",
    "        \n",
    "        for i in xrange(len(weights)): # loop over each weight\n",
    "            # Recall that feature_matrix[:,i] is the feature column associated with weights[i]\n",
    "            # compute the derivative for weight[i].\n",
    "            #(Remember: when i=0, you are computing the derivative of the constant!)\n",
    "            if i== 0:  feature_is_constant = True\n",
    "            else:      feature_is_constant = False\n",
    "            derivative = feature_derivative_ridge(errors, feature_matrix.loc[:,feature_matrix.columns[i]], \n",
    "                                                  weights[i], l2_penalty, feature_is_constant)\n",
    "            # subtract the step size times the derivative from the current weight\n",
    "            weights[i] -= step_size* derivative\n",
    "            \n",
    "    print 'Done with gradient descent at iteration ', iteration\n",
    "    print 'Learned weights = ', str(weights)\n",
    "    return weights\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ridge_regression_gradient_descent(feature_matrix, output, initial_weights, step_size, \n",
    "                                      l2_penalty, max_iterations=100):\n",
    "    print 'Starting gradient descent with l2_penalty = ' + str(l2_penalty)\n",
    "    \n",
    "    weights = np.array(initial_weights) # make sure it's a numpy array\n",
    "    iteration = 0 # iteration counter\n",
    "    print_frequency = 100  # for adjusting frequency of debugging output\n",
    "    \n",
    "    #while not reached maximum number of iterations:\n",
    "    while (iteration < max_iterations):\n",
    "        iteration += 1  # increment iteration counter\n",
    "        ### === code section for adjusting frequency of debugging output. ===\n",
    "        #if iteration == 10: print_frequency = 10\n",
    "        if iteration == 100:\n",
    "            print_frequency = 100\n",
    "        if iteration%print_frequency==0:\n",
    "            print('Iteration = ' + str(iteration))\n",
    "        ### === end code section ===\n",
    "        \n",
    "        # compute the predictions based on feature_matrix and weights using your predict_output() function\n",
    "        predictions = np.dot(feature_matrix, weights)\n",
    "        # compute the errors as predictions - output\n",
    "        errors = predictions - output\n",
    "        # from time to time, print the value of the cost function\n",
    "        if iteration%print_frequency==0:\n",
    "            print 'Cost function = ', str(np.dot(errors,errors) + l2_penalty*(np.dot(weights,weights) - weights[0]**2))\n",
    "        \n",
    "        mat= np.eye(len(weights), len(weights))\n",
    "        mat[0,0] = 0\n",
    "        derivative = 2 * np.dot(errors, feature_matrix) + 2* l2_penalty* np.dot(mat, weights)\n",
    "        # subtract the step size times the derivative from the current weight\n",
    "        weights -= step_size* derivative\n",
    "            \n",
    "    print 'Done with gradient descent at iteration ', iteration\n",
    "    print 'Learned weights = ', str(weights)\n",
    "    return weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing effect of L2 penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The L2 penalty gets its name because it causes weights to have small L2 norms than otherwise. Let's see how large weights get penalized. Let us consider a simple model with 1 feature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "simple_features = ['sqft_living']\n",
    "my_output = 'price'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us split the dataset into training set and test set. Make sure to use `seed=0`:\n",
    "```python\n",
    "train_data,test_data = sales.random_split(.8,seed=0)\n",
    "```\n",
    "In this part, we will only use `'sqft_living'` to predict `'price'`. Use the `get_numpy_data` function to get a Numpy versions of your data with only this feature, for both the `train_data` and the `test_data`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(simple_feature_matrix, output) = get_numpy_data(train_data, simple_features, my_output)\n",
    "(simple_test_feature_matrix, test_output) = get_numpy_data(test_data, simple_features, my_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set the parameters for our optimization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "initial_weights = np.array([0., 0.])\n",
    "step_size = 1e-12\n",
    "max_iterations=1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's consider no regularization.  Set the `l2_penalty` to `0.0` and run your ridge regression algorithm to learn the weights of your model.  Call your weights:\n",
    "\n",
    "`simple_weights_0_penalty`\n",
    "\n",
    "we'll use them later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting gradient descent with l2_penalty = 0.0\n",
      "Iteration = 100\n",
      "Cost function =  1.20825113613e+15\n",
      "Iteration = 200\n",
      "Cost function =  1.2082511289e+15\n",
      "Iteration = 300\n",
      "Cost function =  1.20825112168e+15\n",
      "Iteration = 400\n",
      "Cost function =  1.20825111445e+15\n",
      "Iteration = 500\n",
      "Cost function =  1.20825110722e+15\n",
      "Iteration = 600\n",
      "Cost function =  1.2082511e+15\n",
      "Iteration = 700\n",
      "Cost function =  1.20825109277e+15\n",
      "Iteration = 800\n",
      "Cost function =  1.20825108554e+15\n",
      "Iteration = 900\n",
      "Cost function =  1.20825107832e+15\n",
      "Iteration = 1000\n",
      "Cost function =  1.20825107109e+15\n",
      "Done with gradient descent at iteration  1000\n",
      "Learned weights =  [ -1.63113515e-01   2.63024369e+02]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ -1.63113515e-01,   2.63024369e+02])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "l2_penalty = 0.0\n",
    "simple_weights_0_penalty = ridge_regression_gradient_descent(simple_feature_matrix, output, initial_weights, step_size, \n",
    "                                      l2_penalty, max_iterations= max_iterations)\n",
    "simple_weights_0_penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's consider high regularization.  Set the `l2_penalty` to `1e11` and run your ridge regression algorithm to learn the weights of your model.  Call your weights:\n",
    "\n",
    "`simple_weights_high_penalty`\n",
    "\n",
    "we'll use them later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting gradient descent with l2_penalty = 1e+11\n",
      "Iteration = 100\n",
      "Cost function =  4.48484123471e+15\n",
      "Iteration = 200\n",
      "Cost function =  4.48483174081e+15\n",
      "Iteration = 300\n",
      "Cost function =  4.48482224695e+15\n",
      "Iteration = 400\n",
      "Cost function =  4.48481275313e+15\n",
      "Iteration = 500\n",
      "Cost function =  4.48480325935e+15\n",
      "Iteration = 600\n",
      "Cost function =  4.48479376561e+15\n",
      "Iteration = 700\n",
      "Cost function =  4.48478427191e+15\n",
      "Iteration = 800\n",
      "Cost function =  4.48477477824e+15\n",
      "Iteration = 900\n",
      "Cost function =  4.48476528462e+15\n",
      "Iteration = 1000\n",
      "Cost function =  4.48475579104e+15\n",
      "Done with gradient descent at iteration  1000\n",
      "Learned weights =  [   9.76730382  124.57217567]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([   9.76730382,  124.57217567])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simple_weights_high_penalty = ridge_regression_gradient_descent(simple_feature_matrix, output, initial_weights, step_size, \n",
    "                                      l2_penalty= 1e11 , max_iterations= max_iterations)\n",
    "simple_weights_high_penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This code will plot the two learned models.  (The blue line is for the model with no regularization and the red line is for the one with high regularization.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x10bf1a910>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEXCAYAAAD4LtBgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8VNX5/9/PzCQhCApGWhCkYBVkD7IGWohi2aqYVqu1\ntKBUIIrwxX5rxLqUSgt++dpK3QFFpD/rBorUalWQWGxiWZRvVZBFiwgC2iDImvX5/TH3Xu5MZiYz\nk0xmkpz363Vfuffcs91huJ855zzneURVMRgMBoMh1fAkuwMGg8FgMITCCJTBYDAYUhIjUAaDwWBI\nSYxAGQwGgyElMQJlMBgMhpTECJTBYDAYUhIjUAYDICLfFZFtye5HIhGRoSKyQ0SOikhejGVVRM6z\nzh8VkTujKNPRassbb58NTRsx+6AMhqaBiKwBVqnqH63rQuD/qepjUZRV4HxV3ZnYXhoMpzAjKEOT\nR0R8ye5DPfEt4MNkd8JgiBYjUIZGiYjsEpHbRGSLiHwlIk+ISDPrXq6I7BGRW0VkP/CEneYqf46I\nvCAiX4pIiYg86Lo3SUS2WvW+JiLfCtOHV0XkpqC0/xORH4qf+0TkCxH5WkTeF5GeUTxX2HIikiUi\nq6z09SIyR0Tetu59DJwL/MWadpsHfBd40Lp+MHyr1fqwVER+a51vFZFLXfd81md2oYh0sqYGfda9\nQqtP/xCRIyLyuoic5So7QUQ+tT7vO61/w0ui7Zeh8WEEytCYGQ+MAr4NdAHucN1rC5yJf1QxxV3I\nWjN5GfgU6AS0B56x7l0O/Ar4IdAGWAc8Hab9p4FrXPV2t9r7KzASGGb16wzgKqAkimeKVO4h4CTQ\nDphkHQCo6reB3cBlqtpCVW+z+n6TdR0gpDEQ8Iz4P+//qOq7YfL/BLgO+AaQDvwSnM/mYfz/Zu2s\nZ2sfZ58MjQQjUIbGzIOq+pmqHgR+R+CLtAr4taqWquqJoHIDgbOBW1T1mKqeVNW3rXv5wDxV3aqq\nFcBcIDvMKOrFoHvjgRdUtRQoB1oCF+BfC96qqvuieKaQ5SxRvQK4y+rzB8CTUdRXW/4MjBOR5tb1\nTwgv2ABPqOp26zN/Dsi20q8E/qKqb6tqGXAXYBbImzhGoAyNmc9c55/iFx2bL1X1ZJhy5wCfWgIU\nzLeAP4rIIRE5BBwEhBC/9lX1CP7R0o+tpGuAp6x7bwIP4h/1fCEii0Tk9JoeKEK5NoAvxDMnFMto\nYitwmSVS4/CLVjj2u86PAy2s87Nx9V1VjxPdiNLQiDECZWjMnOM67wh87rqO9Ov8M6BjGOOJz4Cp\nqtrKdWSqalGYup4GrhGRHKAZsNbpgOr9qtoP6I5/yu6Wmh8pbLkvgQqqP3PEqqJpLwrsab7LgS1x\nWvrtAzrYFyKSCWTVTfcMDRUjUIbGzDQR6SAiZwK3A89GWW49/hfmPSJymog0E5Gh1r1HgdtEpAeA\niJwhIj+KUNcr+EdddwPPqmqVVW6AiAwSkTTgGP61o6qaOhaunKpWAi8As0WkubWmM7GG6g7gN5yo\nLc/gXxu7gcijp0gsxz8KGyIi6cBs/CNTQxPGCJShMfNn4HXgE+Bj4LfRFLJe9pcB5+E3LNgDXG3d\nexH4H+AZEfka+AAYE6GuUvzCcQmBL+/TgcXAV/in4kqA/wUQkV+JyKthqgxbDrgJ/5TZfmAp8EQN\nj/pH4ErLGvH+GvKGxVo7KwaGEP2PgOA6PgSm4xe7fcBR4AugNN5+GRo+ZqOuoVEiIruA61V1dbL7\nkixE5Fr8n8F3kt2XWBGRFsAh/JuD/53s/hiSgxlBGQyGlEBELrOmJ08D7gXeB3Ylt1eGZGIEymAw\n2L4Ij4Y66rEbl+M3ZPkcOB/4sZopniaNmeIzGAwGQ0qS0BGUiNwsIh+KyAci8rRlDdVZRP4pIjtF\n5FnLYgcRybCud1r3O7nquc1K3yYio1zpo620nSIyy5UecxsGg8FgSC0SNoISkfbA20B3VT0hIs/h\nN7kdi383/TMi8ijwf6r6iIjcCPRW1XwR+THwA1W92jKXfZpTu/tX49/7AbAd+B5+K6sNwDWqusVq\nK+o2Ij3HWWedpZ06darTz8ZgMBgaO5s2bfqPqrapTR2J9uLsAzJFpBxojt989GL87lDA74plNvAI\n/vnn2Vb6cvxOLMVKf8Yy1/23iOzEL1YAO1X1EwAReQa4XES2xtpGpHnuTp06sXHjxjgf32AwGJom\nIlJrTyYJm+JT1b34LXF24xemw8Am4JDLhcweTrmIaY/l6sS6fxj/TnInPahMuPSsONowGAwGQ4qR\nMIESkdb4Ryyd8U/NnQaMTlR7dYmITBGRjSKy8csvv0x2dwwGg6FJkkgjiUuAf6vql6pajn83/VCg\nlcvHWQdgr3W+F8uPmHX/DPy75J30oDLh0kviaCMAVV2kqv1VtX+bNrWaQjUYDAZDnCRyDWo3MNjy\ncHwCGAFsxO8s80r8Lk0mAi9Z+VdZ18XW/TdVVUVkFfBnEfkD/pHY+fh9pQlwvoh0xi88PwZ+YpWJ\nqY1YH6y8vJw9e/Zw8mQ4Z9iGpkCzZs3o0KEDaWlpye6KwdAoSZhAqeo/RWQ58C5+L8vvAYvwhx94\nRvwROd8DHreKPA78yTKCOIgVokBVP7Ss8rZY9UyzfKUh/milrwFeYInlzwvg1ljaiJU9e/bQsmVL\nOnXqhN+Ow9DUUFVKSkrYs2cPnTt3TnZ3DIZGidmoWwP9+/fXYCu+rVu3csEFFxhxauKoKh999BHd\nunVLdlcMhpRDRDapav/a1GFcHcVJTeJ09OhR9u3bx9Gj9ekpxlCfmB8ohkRQXFzMvHnzKC4uTnZX\nkk6i90E1SY4ePcr27dupqqrC4/HQpUsXWrRoUXNBg8HQpCkuLmbEiBGUlZWRnp7OmjVryMnJSXa3\nkoYZQSWAI0eOUFXljz1XVVXFkSNHktaXu+66i9Wrm2zECYOhQVFYWEhZWRmVlZWUlZVRWFiY7C4l\nFTOCSgAtW7bE4/E4I6iWLVsmpR+VlZXcfffdSWnbYDDETm5uLunp6c4IKjc3N9ldSipmBJUAWrRo\nQZcuXWjfvr0zvVfX88q7du3iggsuYPz48XTr1o0rr7yS48eP06lTJ2699VYuvPBCnn/+ea699lqW\nL18OwIYNGxgyZAh9+vRh4MCBHDlyhMrKSm655RYGDBhA7969WbhwYZ30z2AwxE5OTg5r1qxhzpw5\nTX56D8wIKmG0aNHCWXdK1Lzytm3bePzxxxk6dCiTJk3i4YcfBiArK4t3330XgL/97W8AlJWVcfXV\nV/Pss88yYMAAvv76azIzM3n88cc544wz2LBhA6WlpQwdOpSRI0ca02mDIUnk5OQ0eWGyMSOoeiBR\n88rnnHMOQ4cOBeCnP/0pb7/9NgBXX13dQfu2bdto164dAwYMAOD000/H5/Px+uuvs2zZMrKzsxk0\naBAlJSXs2LGjTvpnMBgMtcGMoOqBRM0rB5s529ennXZa1HWoKg888ACjRo2qObPBYDDUI2YEVQ8k\nal559+7dzprWn//8Z77zne+Ezdu1a1f27dvHhg0bAL+lYUVFBaNGjeKRRx6hvLwcgO3bt3Ps2LE6\n6Z/BYDDUBjOCqicSMa/ctWtXHnroISZNmkT37t254YYbeOCBB0LmTU9P59lnn2X69OmcOHGCzMxM\nVq9ezfXXX8+uXbu48MILUVXatGnDypUr67SfBoPBEA/G1VENhHN1lGz3Nrt27eLSSy/lgw8+SGo/\nmjqp8F0wGFIR4+rIYDAYDI0WI1ANlE6dOpnRk8FgaNQYgTIYDAZDSmIEymAwGAwpiREog8FgMKQk\nRqAMBoPBkJIkTKBEpKuIbHYdX4vITBE5U0TeEJEd1t/WVn4RkftFZKeI/EtELnTVNdHKv0NEJrrS\n+4nI+1aZ+8VypRBPG3VNqgQsNOE2DE0ZE/yvYZOwjbqqug3IBhARL7AXeBGYBaxR1XtEZJZ1fSsw\nBjjfOgYBjwCDRORM4NdAf0CBTSKySlW/svJMBv4JvAKMBl6NtY26fvZUCVhowm0YmjIm+F/Dp76m\n+EYAH6vqp8DlwJNW+pNAnnV+ObBM/bwDtBKRdsAo4A1VPWiJ0hvAaOve6ar6jvp3Gy8LqiuWNuqU\nUAELTbgNg6F+McH/Gj715erox8DT1vk3VXWfdb4f+KZ13h74zFVmj5UWKX1PiPR42tjnSkNEpgBT\nADp27BjVA7oJDli4ZcsWLrvsMhNuw2CoR0zwv4ZPwkdQIpIOjAOeD75njXwS6mspnjZUdZGq9lfV\n/m3atIm5zeCAhevXrzfhNgyGesYE/2v41McIagzwrqoesK4PiEg7Vd1nTa99YaXvBc5xletgpe0F\ncoPSC630DiHyx9NGneMOWGjCbRgMycEE/2vY1Mca1DWcmt4DWAXYlngTgZdc6RMsS7vBwGFrmu41\nYKSItLas8UYCr1n3vhaRwZb13oSgumJpI6GYcBsGg8EQOwkdQYnIacD3gKmu5HuA50Tk58CnwFVW\n+ivAWGAncBy4DkBVD4rIHGCDle9uVT1ond8ILAUy8VvvvRpPG/WBCbdhMBgMsWHCbdSACbdhiEQq\nfBcMhlTEhNswGAwGQ6PFCFQDxYTbMBgMjR0jUAaDwWBISYxAGQwGgyElMQJlMBgMhpTECJTBYDAY\nUhIjUIY6Zfbs2dx7770xlVm1ahX33HNPzG1t3ryZV155pdb1GAyG1KS+nMUaUpSKigp8vuR9DSoq\nKhg3bhzjxo2LuezmzZvZuHEjY8eOBYi7HoPBkJqYEVQDZNeuXXTr1o3JkyfTo0cPRo4cyYkTJwD/\nS3vw4MH07t2bH/zgB3z11VfVyl977bXk5+czaNAgCgoKOHbsGJMmTWLgwIH07duXl17ye4Y6fvw4\nV111Fd27d+cHP/gBgwYNwt607I5vtXz5cq699tpq7SxevJgBAwbQp08frrjiCo4fPx6y/aVLl3LT\nTTcBkJ2d7RyZmZm89dZbrF+/npycHPr27cuQIUPYtm0bZWVl3HXXXTz77LNkZ2fz7LPPBtSza9cu\nLr74Ynr37s2IESPYvXu30/aMGTMYMmQI5557rhOKxGAwpB5mBFVLZs6EzZvrts7sbFiwIHKeHTt2\n8PTTT7N48WKuuuoqVqxYwU9/+lMmTJjAAw88wPDhw7nrrrv4zW9+w4IQle3Zs4eioiK8Xi+/+tWv\nuPjii1myZAmHDh1i4MCBXHLJJTzyyCO0bt2aLVu28MEHH5CdnR3Tc/zwhz9k8uTJANxxxx08/vjj\nTJ8+vVr7S5cudcpstj7Mv/zlL8yfP58hQ4Zw4sQJ1q1bh8/nY/Xq1fzqV79ixYoV3H333WzcuJEH\nH3wQIKCe6dOnM3HiRCZOnMiSJUuYMWOG48Jp3759vP3223z00UeMGzeOK6+8MqbnMhgM9YMRqAZK\n586dHcHo168fu3bt4vDhwxw6dIjhw4cDMHHiRH70ox+FLP+jH/0Ir9cLwOuvv86qVauctaOTJ0+y\ne/du3n77bf7rv/4LgJ49e9K7d++Y+vjBBx9wxx13cOjQIY4ePRrgMd3dfjA7duzglltuYe3ataSl\npbF//34mTpzIjh07EBHHsW0kiouLeeGFFwD42c9+RkFBgXMvLy8Pj8dD9+7dOXDgQLgqDAZDkjEC\nVUtqGukkioyMDOfc6/U6U3zR4g7JoaqsWLGCrl27Rl3eHerj5MmTIfNce+21rFy5kj59+rB06dKA\nOFjhQoIcPXqUq666isWLF9OunT/Y8Z133slFF13Eiy++yK5du2odrsT92RlflAZD6mLWoBoRZ5xx\nBq1bt2bdunUA/OlPf3JGU5EYNWoUDzzwgPOyfu+99wAYOnQozz33HABbtmzh/fffd8p885vfZOvW\nrVRVVfHiiy+GrPfIkSO0a9eO8vJynnrqqaieYdKkSVx33XV897vfddIOHz5M+/b+YMnuabyWLVty\n5MiRkPUMGTKEZ555BoCnnnoqoD6DwdAwMALVyHjyySe55ZZb6N27N5s3b+auu+6qscydd95JeXk5\nvXv3pkePHtx5550A3HjjjXz55Zd0796dO+64gx49enDGGWcAcM8993DppZcyZMgQZ6QTzJw5cxg0\naBBDhw7lggsuqLEfn376KcuXL2fJkiWOocTGjRspKCjgtttuo2/fvlRUVDj5L7roIrZs2eIYSbh5\n4IEHeOKJJ+jduzd/+tOf+OMf/1hj+waDIbUw4TZqIFXDbdQHlZWVlJeX06xZMz7++GMuueQStm3b\nRnp6erK7ljI0le+CwRArdRFuw6xBGcJy/PhxLrroIsrLy1FVHn74YSNOBoOh3jACZQhLy5YtCR49\nGgwGQ32R0DUoEWklIstF5CMR2SoiOSJypoi8ISI7rL+trbwiIveLyE4R+ZeIXOiqZ6KVf4eITHSl\n9xOR960y94tlWhZPGwaDwWBILRJtJPFH4G+qegHQB9gKzALWqOr5wBrrGmAMcL51TAEeAb/YAL8G\nBgEDgV/bgmPlmewqN9pKj6kNg8FgMKQeCRMoETkDGAY8DqCqZap6CLgceNLK9iSQZ51fDixTP+8A\nrUSkHTAKeENVD6rqV8AbwGjr3umq+o76LT2WBdUVSxsGg8FgSDESOYLqDHwJPCEi74nIYyJyGvBN\nVd1n5dkPfNM6bw985iq/x0qLlL4nRDpxtBGAiEwRkY0isvHLL7+M9nkNBoPBUIckUqB8wIXAI6ra\nFzjGqak2AKyRT0Lt3ONpQ1UXqWp/Ve3fpk2bBPWscWLCbRgMhroikQK1B9ijqv+0rpfjF6wD9rSa\n9fcL6/5e4BxX+Q5WWqT0DiHSiaONJot742uy2h83bhyzZs2qOXMQwQIVbz0GgyE1SZhAqep+4DMR\nsR28jQC2AKsA2xJvIvCSdb4KmGBZ2g0GDlvTdK8BI0WktWUcMRJ4zbr3tYgMtqz3JgTVFUsbDQoT\nbsOE2zAYmgKJ3gc1HXhKRNKBT4Dr8IvicyLyc+BT4Cor7yvAWGAncNzKi6oeFJE5wAYr392qetA6\nvxFYCmQCr1oHwD2xtFEbTLiN8JhwGwaDoTYkVKBUdTMQytXFiBB5FZgWpp4lwJIQ6RuBniHSS2Jt\no6Fhwm1ExoTbMBgaPsaTRC0x4TZMuA2DwZAYjDfzRoQJt3EKE27DYGj4GIFqZJhwG35MuA2DoeFj\nwm3UQLzhNo4ePcqRI0do2bJlgMVbQ8KE26gZE26j7ikuLqawsJDc3FxycnKS3R1DnJhwGynK0aNH\n2b59O1VVVXg8Hrp06dIgRcqE2zDUN8XFxYwYMYKysjLS09NZs2aNEakmjBGoBHDkyBGqqqoAqKqq\n4siRIw1SoEy4DUN9U1hYSFlZGZWVlZSVlVFYWGgEqglj1qASQMuWLfF4/B+tx+OhZcuWSe6RwdAw\nyM3NJT09Ha/XS3p6eq0tNg0NGzOCSgAtWrSgS5cuDX4NymCob3JyclizZo1ZgzIARqASRosWLYww\nGQxxkJOTY4TJAJgpvoRx9OhR9u3bx9GjR5PdFYOh1hQXFzNv3jyKi4uT3RVDE8IIVAI4evQo27Zt\nY+/evWzbtq3ORWrXrl307FnNwxMAd911F6tXr45YPp6QGHVNp06d+M9//hNTmeuvv54tW7bE3NbS\npUv5/PPPa11PYyAeobEt6+68805GjBhhRMpQb5gpvgRQUlLieGVQVUpKSuptuu/uu+9OSL2VlZVh\nfefVB5WVlTz22GNxlV26dCk9e/bk7LPPBoi7noZOvCbcxrLOkCzMCKqBUllZGTLcxrXXXuuEkHjl\nlVe44IIL6NevHzNmzODSSy91ym/ZsoXc3FzOPfdc7r///pBttGjRgv/+7/+mT58+FBcXs2nTJoYP\nH06/fv0YNWoU+/b5I5Vs2LCB3r17k52dzS233OKM7tzhLwAuvfTSAH98Nnl5efTr148ePXqwaNGi\nsO3n5uayceNGVq1a5Xia6Nq1K507dwb84jxgwAB69uzJlClTUFWWL1/Oxo0bGT9+PNnZ2Zw4ccKp\nB+Dpp5+mV69e9OzZk1tvvTWg7dtvv50+ffowePDgRuFUNpTQRIOxrDMkCzOCqi0h4m10qKzkTCv2\nEUDz5s0hltFHFPE2woXbsDl58iRTp07l73//O507d+aaa64JKP/RRx+xdu1ajhw5QteuXbnhhhtI\nS0sLyHPs2DEGDRrE73//e8rLyxk+fDgvvfQSbdq04dlnn+X2229nyZIlXHfddSxevJicnJy4AgYu\nWbKEM888kxMnTjBgwACuuOIKsrKyAtp3M27cOMaNGwfAVVdd5fgbvOmmmxzXTj/72c94+eWXufLK\nK3nwwQe599576d8/cFP7559/zq233sqmTZto3bo1I0eOZOXKleTl5XHs2DEGDx7M7373OwoKCli8\neDF33HFHzM+WSthCY4+gohUaY1lnSBZGoBKA1+ulefPmVFRW4vN6EzI1FirchpuPPvqIc8891xld\nXHPNNQGjk+9///tkZGSQkZHBN77xDQ4cOECHDh0C6vB6vVxxxRUAbNu2jQ8++IDvfe97gH8E165d\nOw4dOsSRI0ecl9ZPfvITXn755Zie5f7773cczn722Wfs2LGDrKysgPZDMX/+fDIzM5k2zR9BZe3a\ntcyfP5/jx49z8OBBevTowWWXXRa2/IYNG8jNzaVNmzYAjB8/nr///e/k5eWRnp7ujDj79evHG2+8\nEdMzpSK1ERpjWWdIBkagakuYkc6JBPviq224jeDyoUK/N2vWzBFXVaVHjx7VFsgPHToUtg2fz+d4\n1IDQYTkKCwtZvXo1xcXFNG/enNzcXCefu/1gVq9ezfPPP8/f//53p+4bb7yRjRs3cs455zB79uyw\nYUCiIS0tzQkpEu7zaYjEKzTGP54hGZg1qARg++Lbu3cv27dvT4qpedeuXfnkk0+ckVWwt+946vvy\nyy8dgSovL+fDDz+kVatWtGzZkn/+858ATogL8Fvqbd68maqqKj777DPWr19frd7Dhw/TunVrmjdv\nzkcffcQ777xTY18+/fRTpk2bxvPPP09mZiZwSvzOOussjh49GhDKPVxYjoEDB/LWW2/xn//8h8rK\nSp5++umowpM0NYwVnyFZJFSgRGSXiLwvIptFZKOVdqaIvCEiO6y/ra10EZH7RWSniPxLRC501TPR\nyr9DRCa60vtZ9e+0ykq8bdQloXzx1TeZmZk8/PDDjB49mn79+tGyZUsnVEY8pKens3z5cm699Vb6\n9OlDdnY2RUVFADz++ONMnjyZ7Oxsjh075rQzdOhQOnfuTPfu3ZkxYwYXXlj94x49ejQVFRV069aN\nWbNmMXjw4Br7snTpUkpKSsjLyyM7O5uxY8fSqlUrJk+eTM+ePRk1ahQDBgxw8l977bXk5+c7RhI2\n7dq145577uGiiy6iT58+9OvXj8svvzzuz6ixEq9xhcFQa1Q1YQewCzgrKG0+MMs6nwX8j3U+FngV\nEGAw8E8r/UzgE+tva+u8tXVvvZVXrLJj4mkj0tGvXz8NZsuWLdXS3Bw5ckQ3bdqkGzZs0E2bNumR\nI0ci5k8UdrtVVVV6ww036B/+8IeEtqOqOm/ePJ0xY0ZC2klFavouNAaKioo0MzNTvV6vZmZmalFR\nUbK7ZGgAABu1lhqSjCm+y4EnrfMngTxX+jLr2d4BWolIO2AU8IaqHlTVr4A3gNHWvdNV9R3rw1gW\nVFcsbdQpti++9u3bJzXUxuLFi8nOzqZHjx4cPnyYqVOnJqSdv/71r2RnZ9OzZ0/WrVvX4K3dDIHY\nxhVz5swx4S8M9UqijSQUeF1EFFioqouAb6rqPuv+fuCb1nl74DNX2T1WWqT0PSHSiaONfa40RGQK\nMAWgY8eO0T5rAKngi+/mm2/m5ptvTng7V199NVdffXXC2zEkD2PFZ0gGiRao76jqXhH5BvCGiHzk\nvqmqaolXwoinDUtIF4E/om5COmYwGAyGiCR0ik9V91p/vwBeBAYCB+xpNevvF1b2vcA5ruIdrLRI\n6R1CpBNHGwaDwWBIMRImUCJymoi0tM+BkcAHwCrAtsSbCLxkna8CJliWdoOBw9Y03WvASBFpbVnj\njQRes+59LSKDLeu9CUF1xdKGwWAwGFKMRE7xfRN40bL89gF/VtW/icgG4DkR+TnwKXCVlf8V/FZ2\nO4HjwHUAqnpQROYAG6x8d6vqQev8RmApkInfOu9VK/2eWNowGAwGQ+qRsBGUqn6iqn2so4eq/s5K\nL1HVEap6vqpeYouNZVk3TVW/raq9VHWjq64lqnqedTzhSt+oqj2tMjdZ1nxxtdGQMOE2YsOE2zAY\nGibG1VEjw4TbqI4Jt2EwNEyMq6MEkeiIuibchgm30dQxUX4bP2YEVVtChNuorKxET5yghSoqQmVm\nZmyjDxNuw4TbMEQk3uCLhoaFGUElgIrKSrAi6qLqv65j4gm34cYOt3HWWWc54TaCCRduIzs7m9/+\n9rfs2bMnZLiNWLn//vudkYodbiO4/VCECrcxaNAgevXqxZtvvsmHH34YsV13uA2fz+eE2wCqhdsI\n/nwNycX4B2wamBFUbQkx0ik/epQd27dTVVWFx+OhS5cuZNSxVwkTbsOE22jKxBt80dCwMCOoBNCi\nRQu+8Y1vOMEAk+HyyITbaPjhNkKtsZh1Fz/GP2DTwIygEsCXX37J/v37Adi/fz8ZGRlO1Nb6wh1u\n47TTTgsIPxEPdriNGTNmcPjwYSoqKpg5cyY9evRwwm14PB6GDx8eMtxGt27dwobbePTRR+nWrRtd\nu3aNOdwGwNlnn80rr7zihNto27ZtyHAbmZmZAS92d7gNVeX73/9+yoTbCLXGAph1FxfGP2DjR1SN\nq7lI9O/fX22LL5utW7fSrVu3sGW2b9/O119/7VyffvrpdOnSJWF9DMfRo0dp0aIFqsq0adM4//zz\nE+I81m4H4J577mHfvn388Y9/rPN2UpGavgvxMm/ePO68807HvH/OnDkA1dJuu+22Om/bYKgLRGST\nqvavOWd4zAgqAbRu3TpAoFq3bp2UfixevJgnn3ySsrIy+vbtm9BwG/PmzaOiooJvfetbLF26NCHt\nNCXCrbEanRKTAAAgAElEQVSYdRdDU8KMoGognhEU+Kf5vvrqK1q3bl3v03uG+iNRIyjwT/MVFhaS\nm5vrTGWFSjMYUhEzgkoiqupYeYWiTZs2RpgaOYn+cRdqjcWsuxiaEsaKLw6aNWtGSUlJwl9QhtRF\nVSkpKaFZs2bJ7orB0GgxI6g46NChA3v27OHLL78Mm6e0tJSTJ0/SrFmzgD1HhsZDs2bN6NChQ80Z\nUxAzVWhoCEQtUCLyLeB8VV0tIpmAT1Wrby5pAqSlpTkeGkJh3LAYUhnz/TQ0FKKa4hORycByYKGV\n1AFYmahONXSMGxZDKmO+n4aGQrRrUNOAocDXAKq6A/hGojrV0LFNhL1erzEHNqQc5vtpaChEO8VX\nqqplttWaiPgAYyEQBtsNi5njN6Qi5vtpaChEtQ9KROYDh4AJwHT8oda3qOrtUZT1AhuBvap6qYh0\nBp4BsoBNwM8s8csAlgH9gBLgalXdZdVxG/BzoBKYoaqvWemjgT8CXuAxVb3HSo+5jXCE2gdlMBgM\nhsjUxT6oaKf4ZgFfAu8DU4FXgGiD4/wXsNV1/T/Afap6HvAVfuHB+vuVlX6flQ8R6Q78GOgBjAYe\nFhGvJXwPAWOA7sA1Vt6Y20gExqmnIVmY756h0aCqNR7AaYDXde0FmkdRrgOwBrgYeBkQ4D/4LQAB\ncoDXrPPXgBzr3GflE+A24DZXna9Z5ZyyVvpt1hFzG5GeoV+/fhorRUVFmpmZqV6vVzMzM7WoqCjm\nOgyGeDDfPUOqAGzUKPQl0hHtCGoNkOm6zgRWR1FuAVAA2EGBsoBDqmoH19kDtLfO2wOfAVj3D1v5\nnfSgMuHS42kjABGZIiIbRWRjpL1O4TBWUoZkYb57hsZEtALVTFWP2hfWefNIBUTkUuALVd1Ui/4l\nBVVdpKr9VbV/PO6KjJWUIVmY756hMRGtFd8xEblQVd8FEJF+QE0hXIcC40RkLNAMOB2/QUMrEfFZ\nI5gOwF4r/17gHGCPZSV4Bn5DBjvdxl0mVHpJHG3UKTk5OSxYsIAVK1ZwxRVXGCupJkx9e2wwFnqG\nxkS0AjUTeF5EPse/xtMWuDpSAVW114QQkVzgl6o6XkSeB67Eb2U3EXjJKrLKui627r+pqioiq4A/\ni8gfgLOB84H1Vj/Otyz29uI3pPiJVWZtLG1E+RlETXFxMdOmTaOiooI333yTXr16mRdFIyec53Hj\nscFgqAXRLlYBaUBP60iLZaELyAVets7PxS8wO4HngQwrvZl1vdO6f66r/O3Ax8A2YIwrfSyw3bp3\nuys95jbCHfEYSeTl5Sn+fWIKaF5eXsx1GBoO4QwT5s6dq16vVwH1er06d+7cpPXFYKhvqAMjiYgj\nKBG5WFXfFJEfBt3qIiKo6gtRimAhUGidfwIMDJHnJPCjMOV/B/wuRPor+E3eg9NjbqMu+fzzzyNe\nGxoXbsOEkydPsmzZMnJycsIGHayvvthGEmbUZmio1DTFNxx4E7gsxD0FohKopkZubi7r168PuDY0\nXnJzc/H5fFRWVqKqLFmyhAkTJiRlPSgZomgwJIqIAqWqvxYRD/Cqqj5XT31q8LRq1SritSE1CV5H\nitbAIScnhzFjxrBypd9/cmVlpTNyqe8Ag8ZIwtCYqNFIQlWrRKQAMAIVJR9++GHEa0PqEWzQsGDB\nAmbMmOFcr127NuzLvri4mFdffdW59nq9SR25mKi7hsZCtPugVovIL0XkHBE50z4S2rMGzFtvvRXx\n2pB6FBYWUlpaSmVlJaWlpTz++OOUlpaiqpSWlrJs2bKIZSsq/PvCRYRJkyYZgTAY6oBoBepq/A5i\n38Lv+NU+DCE466yzgFH4/eva14ZUJisri6oqv8OTqqqqmEK5uzfHNmvWjAkTJiSqmwZDwvnkE7js\nMvjWt2D37uT2Jdp9UN3xC9R38BtHrAMeTVSnGjqDBw9m8+bZ+PcBn8XgwYOT3CNDTZSUlODxeKiq\nqsLj8dC9e3feeecdysvLSUtLiyg6Zt3H0NB55RWYPBmCDY6//jo5/bGJNtzGc/iDFT5lJf0EOENV\nr0pg31KCeMJtFBcXM2TIt4FvkJ5+DoWFz5mXVooTalMtYETH0CgpLYV774U7wsSkWLAApk0DX7RD\nmBDURbiNaAVqi6p2rymtMRJvPKjWrcs5dCiNl1/ewPe/PyABPTPUNfXtlshgqE8++wx+8QtYvrz6\nvS5dYOFCqEvbnroQqGj18V0RGayq71gND8KsQUUkLS0NgAEDjDg1FIz1m6GxsXYtTJkCO3dWv3fV\nVfD730OHDvXfr2iJ1kiiH1AkIrtEZBd+X3YDROR9EflXwnrXgCkvLwdgw4YNSe6JIVGYwICGVKO8\nHO67D0T8x8UXB4rT3Llw8iSowrPPprY4QfQjqNEJ7UUjo7i4mEOH/GtQl112Gf/4x4vml3k9Uh9T\ndXXtCDbePtf1s0aqz0yBpiYHDkBBAYTaCdGhAyxaBGPG1H+/6oTaOvNr7Ec8zmLPO+88hQPq/53S\nRs8777yY6zDER305S61LR7Dx9rmunzVSfcYJbWpRVKTas6da75jA47LLVD/5JNk9rBtnsdFO8Rli\nYHfQ5oHga0PiiCaibKipuVin6+oyMGDwJuFoo+CGelb7ORYtWhTz9GOkz85E6k0ulZXw6KOnpu6G\nDIEPPjh1/6674Phxv0StWgWdOyevr3VJLYwIDeFo0aIFBw8GXhvqh2BnqVlZWcybNy/Av14oc/JY\np+vqcu9T8CbhrKysuJ91xIgRlJaWOvu5MjIyop5+jORo1jihrX9KSuBXv/JP0QVz1ln+9Lw8v2A1\nVoxAJYBjx45FvDYkDrdwZGVlMXPmzAD/eitWrHBe4O6RQDwhKurK6i94k3BJSXRBnoNF0h7luMUu\n1ucJJ7pmM3L98O67cMMN4AqG4HDJJfDAA3DBBfXfr6RR2znCxn7EswbVvn37gDWo9u3bx1yHITJF\nRUU6d+7ciGsh7nUij8ejPp9PPR6Pc22vpdS0vhJNW7V9lrpY37HrCfWMhtSkqkp16VLV9PTQ60m/\n/KXq118nu5fxQR2sQSXsxY4/eu164P+AD4HfWOmdgX/ij2r7LJBupWdY1zut+51cdd1mpW8DRrnS\nR1tpO4FZrvSY2wh3xCNQw4YNCxCoYcOGxVyHITzuF3FaWpouXLgw4J4tJu4Xf1pamvPiFhE977zz\nwpYL1VaijQPqSgTtehYuXJhQUTXEz+HDqjNnhhak5s1Vn3rKL1wNnVQXKAFaWOdpliAMxh+248dW\n+qPADdb5jcCj1vmPgWet8+6WyGVYwvMx4LWOj/GHd0+38nS3ysTURqQjHoFq165dgEC1a9cu5joM\n4Zk7d64jNoD6fL6wIyH3Czue0UUywrYbGh8ffqg6fHhoURo6VPX//i/ZPax76kKgEmbFZ/XxqHWZ\nZh0KXAzYzjaeBPKs88uta6z7I0RErPRnVLVUVf+Nf/Qz0Dp2quonqloGPANcbpWJtY06xQ69EO7a\nUDtyc3Pxer3OdVVVFYWFhWHDnefm5lJSUsKCBQu45JJLEBGqqqqqWcyFsuSz2xKRanGekrFRty7a\nNBuME48qPPccnHGG34ihRw9wR9256SY4eNCf7+23oXfv5PU1lUmokYSIeIFNwHnAQ/hHPIdU1X5j\n7wHaW+ftgc8AVLVCRA4DWVb6O65q3WU+C0ofZJWJtY3/1PphXVRWVka8NtSOnJwcbr75Zu69915U\nlYyMDEc4gi3NFi1axE033URlZSUZGRlMnz6d119/HQi0mLOt+0pLS/F6vTz44INMmTIF8Md4cv91\n56+rjbp2nZGMEOqizUT02+Dn2DH47W/hnnuq3/N4/FZ3113nPzdER0IFSlUrgWwRaQW8CDQI+xMR\nmQJMAejYsWPM5Q+6bcxDXBviw36BZ2Vl8cADDwDg8/lYsGCB85J1W5oBTJs2zRnBlpaWsnnzZkQE\nVUVEHIs5ey9SVVUVVVVVTJs2jV69erFs2TLKyspQVcrLy5k9ezazZ88OO1oL1+eaLN+iEY5o24xE\nXdQRK43ZA8XOnTB9Ovztb9XvXXghPPIIDBxY//1qLNSLmbmqHhKRtUAO0EpEfNYIpwOw18q2FzgH\n2CMiPvzBlEpc6TbuMqHSS+JoI7i/i4BF4PdmXptnN9QN7he4PUVXVVUVIDIQaPo9b948x+Qa/KHY\ns7OznRGUqjojKHsqz22ivWzZMh577DF7TZWqqipWr17NunXrWLBggZM/XIj34uJicnNznZhSkcQg\nGuGoi71I9b2fqTGO2P7yF78D1v37q9+bNMk/gmrTpv771RhJ2GBTRNpYIydEJBP4HrAVWAtcaWWb\nCLxkna+yrrHuv2kttK0CfiwiGSLSGTgfv3XgBuB8EeksIun4jR5WWWVibcOQ4rhf4LYo1OTFITc3\nl4yMDDweDz6fjwcffJBWrVrhseZY3HuOcnJyePDBB/H5fM4G1/3791dbP7T3Fr333nshp/7cuEdf\nZWVlEcPGR+OZwt6LNGfOnLhf9HVRRyw0Bg8UJ0/CnDmnvDiMGxcoTg884HfSqgqPP27EqS5J5Aiq\nHfCktQ7lAZ5T1ZdFZAvwjIj8FngPeNzK/zjwJxHZCRzELzio6odWwMQtQAUwzZo6RERuAl7Db9G3\nRFU/tOq6NZY26hr/L9TAa0PtCP7lv2DBAkpKSiJOG4XaXFpcXExaWhplZWWkpaUFCMGUKVPo1auX\nkz9YUGzDDPvfs6KiAlWloqKi1lNl0W6ErYvNwfUZVqSheqDYvRtuvhleeKH6vW7d/LGTvvvd+u9X\nk6O2ZoCN/YjHzHz8+PEBZubjx4+PuQ5DdWraKxTpvtvcPD09XUVE09PTI5qZFxUVaUZGhoqIZmRk\nBOwtimZ/VHD5aPYkJXpTcDJoKM/0xhuqnTqFNgX/yU9U9+5Ndg8bFqTyPqjGcpiNuqlFNBtqMzIy\nND8/38kTbsNuNPuaIrWXn58f0E4s5cPltfuZnp5eY92JpqEIS7yUlan+7/+GFiRQnT9ftbQ02b1s\nuBiBSlGBSk9PDxCo9PT0mOswVCfSqMW9oRbLW4T9ks/Pz6/m8qg2niES5V0i1DMky1VRYw2v8fnn\n/tFQKEHq1En1tdeS3cPGQ10IlLHITwBmH1RiiLTgbq912AYLqn7DhIULF/LEE084RhUZGRk89NBD\njpEAUKdhKWpDuGdIhmFBYzBusHn7bf9GWRE4+2z4859P3cvLg3//2y9R//43jByZvH4aqmO8mSeA\ntm3bsnfvN6yrNNq2bZvU/jQWIi2420YGy5YtY8mSJZSXlzu/wioqKpg8eTIdO3YMMECI1wQ6KysL\nj8eDqtbpwr/7GZ544gkqKiqSZljQUI0bAMrKIDsbtm4NfX/2bLj1VmjWrF67ZYgDI1AJYMCAAey1\nd17RmgEDzk9md1KaWDZx1mTpZlunTZgwwXnJl5eX4/F46Nu3r2OhZ+d1jxJKS0udTbg1baidMWMG\n5eXleL1epk+fHjC6qO2GVPczJHNza0MLr7FtW+QwFC+95DcPNzQwajtH2NiPeNagWrVq5Zrb7qGt\nWrWKuY6mQDiv5KEW5xcuXKgjR46MygO5TUFBgXo8Hmc9KiMjo5oj2ZocyAa3kZ+f76wRYRlaeDwe\nxwAjlJFGfdDYDRpC8eijodeS7OOvf012D5s2GCOJ1BQo/8vrlEDh+M41uAnllbygoMCxtLPFYuHC\nhQGisHDhwqhiOKWlpQWUE5Fq1ntFRUU6cuTIkJZ9odoIFqhwR20NHKIRnGBP7Y3NoCGYigrV73wn\nsigdOJDsXhps6kKgzBSfIWkEuxeqrKzk3nvvda5tb+PBC/QrVqygpKQkomugZcuWUV5eHlBORBAR\nZ03Fnl7Mzs5m7dq1AAHrLcuWLePkyZOonjJWmDBhAo899liNHurdZWrj0NXr9TJp0iQmTJgQUI87\nj8fjcTxs1Jd/vfri00+hU6fw93/4Q1i+vHGHPW/KGCs+Q70QKsRDsHshO+y5jdfrJSsri+bNmwfU\n1bx5c7KysgJcA2VlZTn1FxcXs2TJkmp9UFU8Hg8LFizg/fffZ/jw4dx+++3Mnz+fiooKPB6Ps6a0\naNEinnjiCf80g9UXey3m+uuvD3Bz5PF4EBHS0tLIy8sL6bIolhAXtjDa4rtw4UJGjBgRUNa9flZZ\nWYnH46nR9VOsJCssxzPPnHIrFEqcnnnm1JhpxQojTo2a2g7BGvthpvhqT7jpOPcUVX5+fsCUnMfj\n0WHDhml6erp6vV71+XzarVs3Z50nMzPT8ewQPMXl3vcUfHi9Xs3Pz1efz1ftnnuPlM/nc6YEAc3P\nzw/5PBkZGU5en8+n+fn51aLZxrKnyPY+Earf7k3FwXXWdQTd+twHVVWletllkafuPv00Yc0bEgRm\nis/QEAi3p8Zt4j1x4sSA0ZOIsG7duoARzDnnnMP27dudekpKSrjtttuYN29eQP1wKi6U1+tl7Nix\nvPzyy85+tP379we0BTgjEHuqzA7JYXP66ac7524Lt927d7N48WJU/ebsCxcupFmzZo7JenFxMbNn\nz3ZCedQ0BVdYWBgwfej2/xfKrD5RVnaJDstx4ABE2n2RmwurV4MrLqWhCWIEypBwQu2pcb8AT548\nyZYtW5zpMnuqzxYIe93oiiuuYN26dQH1FBcXs379ekdwVJW+ffvSt29fVqxYQXZ2Nl9//bXzi6yy\nspK//vWv+Hw+Z1rvF7/4Ba1atSIrK4uZM2c6IT0qKyudPvzhD38gLy+v2ku6b9++pKenO2tVqlpN\nhG1x8ng8NU7BxeIUN5FOXxOxD+qvf4VLLw1/f9EimDy51s0Y6ohUiOMl7l+Jhur0799fN27cGFMZ\n//qE/bn2BD6kqX/O7mCDJSUlZGVlMX36dGfEY+PxeLjmmmt44YUXHAOAvn378vOf/9wJIAh+YXjv\nvfdCGizYa1oVFRUhR0MiwtSpU6tt3AVYtGiRI2z33XefY2ghIgwYMIALL7yQvn37OkJmi8h7770X\nsLnWHt3ceeedzhrRJZdcUuM+K/dnVZcvhnjqrG0/VP3xkZYuDZ9n2zbo0iXmqg0Jpi7ieInIJlXt\nX6uO1HaOsLEfZg2q7ghe18jLywtY58G1FjR+/Hjt3r27ejwe9Xg8AfuYMjIyHI/kwWWjOdLS0px9\nUG6Hr8H9KygoqLYWZfcvnFl6uLUn996o+t6zVJ/rSQcPRl5L6t3b76TVkNq4/UJG41Q5FJh9UEag\nGhLBX/r8/PyAjbKRDhFxhCJeYbIPWyj8Tn1PiZbbuML+T7lw4UI977zzQhotRGv0kJ+f7xh72OJa\nn3uW6uJlE4nCwsiiNH9+nTZnqAfq4kdNXQiUMTM31BvuqLH24v+CBQu45JJLaixrr9/YPvBs7Ok8\nm06dOjnXbhNwu47MzEzHjZB7n1R5eTn79++vZro+c+ZMPv7442r9ueyyy6KKSpuTk0PHjh0dc/Cy\nsjLKy8vr1QlrNNF6Y+W///uUKXio6jZvPiVRt9xS6+YM9Ux9R14OhzGSMNQbwc5QFy1a5Bgp2EYT\nofB6vTz88MP06tWL2bNns3r1amdt6frrrw/wWwdw0UUXOetXzi+xEGtA9gZXm7Zt2wZYxhUWFlJa\nWhogiABpaWmMGTPGCRdfE26DA5/Ph6rfWKO+nLDWhcXfV1/BWWdBkPGjQ4cOsH07ZGbWsrOGlKE+\nIy+HpbZDsHAHcA6wFn+o9g+B/7LSzwTeAHZYf1tb6QLcD+wE/gVc6KpropV/BzDRld4PeN8qcz+n\njD5ibiPcEc8Un386xf79OEC9Xm/MdTRmgl0c2X748vPzNS8vT/Py8pw9USKiBQUFTtloXBzZ61O2\nj79weRcuXOj0w16XCq4r2F0SoHl5eTFPf7jXnRqK37wXX4w8dXf77cnuoSGVIcX3QVUA/62q74pI\nS2CTiLwBXAusUdV7RGQWMAu4FRgDnG8dg4BHgEEicibwa6C/9YLYJCKrVPUrK89k4J/AK8Bo4FWr\nzqjbSOBnAJyV2OpTgHDWXsHpbks+t9eIyspKSkpKeOSRRwC44YYbnOk3VeW+++4LMPGeOHEiW7Zs\n4eTJk7z//vsBde/evZuKigr3DyWmTJkS4CrI3Q+fz0d5ebkzDejG9nRx4403OiOt9PR02rZtG3GP\nUKjPI/jXaNJ/mYahJq8M//gHDBlSP30xGOrN2AB4CfgesA1oZ6W1A7ZZ5wuBa1z5t1n3rwEWutIX\nWmntgI9c6U6+WNuI1O/aG0lc26iNJCJ5iQj2dBBsIRfsFNYuF+xJwePxOCOOzMzMakYS48ePD7CW\nCza66N69u+bn52tBQYEOHDgwwFuE25P5yJEjQ45q3NZ+9kgvnKFDvIvLyRpV1WR1B6pffFGvXTI0\nEkjxEZSDiHQC+uIf6XxTVfdZt/YD37TO2wOfuYrtsdIipe8JkU4cbexzpSEiU4ApAB07dozuIcPS\nspblU5dIXhKCPRGsWLHCyXfixAm2b9/OW2+95YxkbGOBZcuWBaxFiQgZGRkBm3tVA9eEnnrqKee8\nqqqKVq1a8dVXXzlpW7ZsYcuWLdX67x41VVVV8cYbb/Dmm2/y0EMPMWXKFOeePfpx7w3x+XxMnjy5\nmhPXeDww1MWek1j45S8/4fe/PzdinqKi4pQd5RmaDgkXKBFpAawAZqrq1+6Xgqrav4YTRjxtqOoi\nYBH4N+ompGMNHPulaouOTVZWlvPXtqJLT0+nTZs2AflWrlxJly5d+Pa3v820adOoqqpyDAhsAfL5\nfFx//fWcfvrpzJ49m+zsbNLT06u16UZVA8QpEsFCp+p3V3TTTTfRq1evatN2bjEG/4+X4Jd4PB4Y\nahK1uti4Gzh1V12cfvADuOWWU0I5YkTihdJgqImECpSIpOEXp6dU9QUr+YCItFPVfSLSDvjCSt+L\n37DCpoOVthfIDUovtNI7hMgfTxsJpHViq08S9kvVLRSVlZVMnz4dgJkzZzqWdl27dg0Y5djMnz/f\n8X8HOC6GbC699FL69u3L1KlTAXj99dcZNmwY3bt355133mHz5s1x979t27YcOHCgmkjZz+EWiWAx\nFhF8Pl818bGFJNg9UU0CE0nU4h1dHTsGLVpEzuPxDOC3v/0ht912GwDz5iXW/57BEDO1nSMMd+C3\nmFsGLAhK/19glnU+C5hvnX8fv4GDAIOB9Vb6mcC/8b/pW1vnZ1r31lt5xSo7Np42Ih21X4Oa0yjX\noMKtB4mIjhw5Mqw38eDDXd62urOv09PTtXv37tXKpKWlaXZ2dlT1hztsz+X2tb0uFRzdVzVwo6u7\nb9GsPdXkyb2m6+DNw/n5+WHXqhYvrnk96R//CL9GVp8eJwyNH1LZkwTwHes/87+AzdYxFsgC1uA3\nAV/tEhsBHgI+xm863t9V1yT8puE7getc6f2BD6wyD3LKzDzmNsIdRqDCU1RUpHl5eQFGCenp6Y4B\nRE0eH3w+X4Ag+Xy+APdHXq9Xu3XrVishinTk5eUFuDpauHChIwZuAQonxnl5ec79gQMHhozYG8qL\nQ01CsHDhQseAw+15IlTY+poEqVmz0P9u4USuoZjAG1KflBaoxnIYgTpFOF9zHo9HvV6v5uXlOdZ6\ndpp75CEi2qlTJx02bJgjCnl5ec59uw73C72goCAqV0jxCpT7eYLDubvv25Z8wXGkCgoKAlwmAZqR\nkeHUaYuNiDgjM/c+MNtC0f0ZB8fFskdN+fn56vGk1yhKL7xQ/98NgyEYI1BGoOqNUI5UzzvvvAD/\neAMHDrReoh4nzX1uT6PZ5QcOHBggYPYUm/tl7/F4Yva916pVq6jKBPvTc4ul3WcRcQSnqKhI27dv\nH5DH/RnYhx3cMNTIKyMjQwsKCgLyB08nuj8Dn8+n//M/W2sUJeOA1ZBq1IVAGVdHCeIMDjGct1iV\n7I7UEW5Ls9LSUubPnx9wX1VZv34969evD0jzf0/953ZYjJMnT1YrDwS4HbIJZ60XiUOHDkWVz22c\nsWzZsmpGF3bfS0tLmTVrFsXFxQH++wAGDRrE7t27HdN4r9dL3759AUK6SiotLaWwsNDZqOzxeAJc\nJuXm5pKRkcGJE/uAM6iogFtvDd1/V7UGQ91x5AisWQNffAHXXQdpaUnrinEWmyDe5GJeIo+WnEx2\nV+qE3NxcvF5vtdhKsSAitSqfKHw+H4sWLWLXrl1h86xbt66aOIkIPXr0oLCwkLy8PLxeL6rKzJkz\nKS4uJisrK6TAnn322WRkZOD1ep09Xqp+U/AhQ3I4ceI4cEa1cg8/HDhuMhji5sgRWLkSpk6Fjh1P\nef4VgdNP9+87mDoV/v73pHbTCFSC6GDtIW5GRQ05Gw62Cbg3KA538HU43COqVKF58+Z069Yt4kjN\nFh43Ho+HZs2aOabjAwcOBPwjvtLSUmbPns17770X4Gkd/K6SCgoKWLNmDePGzeXEieMMGZKDJ8z/\nxCNHTgnSDTfU7lkNTYxoRGjRIvjss+plR4+GJ56AESPqv98uzBRfwqlhM0qKE87HXV5eHsePH+eK\nK66gV69ezJo1i78n+ddWPBw/fjzsfiqPx8OUKVM4cuRIwD6u7OxsBg8eHOBFwt7LZO+VWr16NT6f\nD6/X64if1+ulrKzU5csu9B6jFNNwQypjT8e9+qr/CCU24Rg9GsaM8R/nn5+4PtaG2i5iNfYjXiOJ\nA7RRBW3DEw3WSMJtGBEc0dYdHTbY1LwhHG3atKlmEBF8bpu9B++BEpGw5uFuowmv16sDBw6s0cAB\n5jiGFZH+LYz5dxPlyBHVlStVp0xRPeecmr5MgceoUaoLFqhu317v3cZY8TUEgXq+wQlUpE2i+fn5\nzh4nd+jzVDnS0tK0ZcuWNeYLNg33eDw6cODAqMU2ODKt2+QehtT43hg//raA+iIJVDQhRox4NXBs\nEZo6VbVjx/hEaNs21aqqZD+JQ10IlJniSzgNK4Kb27WO1+t1rM28Xi8TJkxg2bJlAaEwUo3y8vJq\nxt1EpSEAABcVSURBVAyhCA6OWFVVxYYNG/D5fGHXo9yfR7BLoqFDB6F6PGKbbgesxcWX8fzzv6e8\nvJy0tDQmTJgQtlwkX3317WjWUAuOHg2cjtu9O/qyo0YFTsfVFBelkWAEqokT7CfO/TK0X9Sq6hhI\n7N+/P5ndTSiqp0zhgxk4cCALFiwAcD6vIUPcQlDdyiEt7SQ///nNrrWqwHhQhYWFUTmBjeSrLx7v\n6YYEUhsRGjkSxo5tciIUkdoOwRr7UfspvndSdoovVMwm97Re8NGtW7dqnhQa2xHq2b1erxYVFem/\n/lXzbMugQdO1oKAgIKZVsNukeKbjwpUz/vOSQG2m40aOVL3vPtWPPkqp6bhEgJniawi0S3YHAnCP\nmII33950000Rp8e2bt1ajz1NDlpt2lKprIwcRbaoqBiAESNGsHFjGe++6wn4HO2RjZ3HPR0HhBxF\nBY9sgyPy2uTk5LBmzZpah+MwBHH0KLz5pn8U9MorsY+E7Om4Ll3MSKgWGIFKOM2S3QGH4PWKBQsW\nOFNHIhLV2k1jwDb9ri5GtueKmtfWvF4flZWVeL1eCgvnsHv3bk6ePOn88nOHtAe/d4vg4I7Lli3j\nySefpLS0FI/H4wRKjHVdKZx4GWrAiFDKYwSqCWG73rE3k7766qv06tWLs88+mzFjxpCfn5+Shg91\nTb9+/WjVqhWvv/66ldKBwEDLobiRYcM+ZPTo0WRlZTFzpl/YPR4P/+///T+2bdvmfHY+n4+xY8fy\n0ksvOWm///3vUVXHvVF6ejqA8+9RVVXlBEoMt65UF4ELmxwlJfD22/C3v/mF6NNPoy/7ve/5BWjs\nWCNCyaK2c4SN/Yh3Dcqecz6HT5O+BmXvVerUqVPYtZe8vDw9//zzk74GVB+Hf2/Sf6JYMqhe1uPx\nOOt1wc5lwb9Hyt4jFmzKbpcfOXKks4fMvaZnezYPta5k1poisH+/6vXXx7YWZB/f+57qH/6gunVr\no18Tqm+ogzUo4+oowZzO10ltv7i4mGHDhrFy5cqIvuZWrlzJzp07669jScGvBevX/xN/yLDqzJ07\nDxEP/tBh1bFHnytWrGD79u3V7jdr1syx2ps0aVJAhGDwe6c499xzef/99yksLOQXv/gFaWlpeDwe\n0tLS2G1NM61Zs4Y5c+Y403uhRlVNii++gClTAt312EfbtvDYY+HLtm0Lf/gDbN0KVVWBEvX663Dz\nzXDBBWaElIrUVuEa+1HbEVQP3k/qCGru3LlJH7Ek72gZxY/ox538Xq83IGihfbg3JLvDiwSH2ejW\nrVvICLXhNv66R2P5+flOYMJwm3EjRcKNxjLQjmllj/BSjgMHVCdPjm8kBH4vCy++aEZCKQLGk4QR\nqJooKioK+YIM5dqncRz/juJdVn3qzT6GDRumw4YNcz6TtLQ0J/JucAws95GWluZMxQUHdRw5cmRY\nkbI9UoSKvBvq3zJYiKKd+isqKgpr+l6v1FaEOnTwR2Q0IpTypLRAAUuAL4APXGlnAm/gD8X+BtDa\nShfgfvwh3f8FXOgqM9HKvwOY6Ervhz9s+06rrMTbRqSjoQmU+yUWKgif/XKy/eq5o7c23COad1vs\n9bqDLNqh14OFRkQC1pRCiYWdHqqN2q4xRSNsdr7gHyXh8taaAwf8fuPiFaH27VVXrDAi1MBJdYEa\nBlwYJFDzgVnW+Szgf6zzscCrlogMBv6pp8TmE+tva+vcFpz1Vl6xyo6Jp42ajoYkUKGcu4Z6KQ4b\nNkzHjx+vrVq1SgFxiefwJkyUahKsUOl2RNyaQrm7/Rj6fL5qU23xbOJN2gjKiJChBlJaoPz9o1OQ\nQG0D2lnn7YBt1vlC4JrgfMA1wEJX+kIrrR3wkSvdyRdrGzU9Q0MSqPz8/IhrJA37eDqK917bhPbB\nHk25R5221d7cuXOrhXIvKCio9m+UCMeuCVuD+uQT1Xbt4hehs882ItSEaYgCdch1LvY18DLwHde9\nNUB/4JfAHa70O620/sBqV/p3gZfjaSNMv6cAG4GNHTt2jOcfxvlP+j1eqxeBCl7c93g8YV0WNZwj\nmvdg/fRFRHTYsGE6cOBAzc7ODriXlpamXq83pHFFXl5eSDFIGQ/k69fHL0DgF7Dly40IGapRFwKV\ntI26qmr/wk+5NlR1EbAIoH///rXqYxe280ZtKohAcXEx8+fPZ/PmzdVMyCNFiE1tovm4E2sOLCJ8\n+9vf5uOPP7Z/sKCqIQMyiogTyDGYqqoqVq5cyauvvsratWvr1QN5wKbe8nIYPrx2Ff7mN3DnncYU\n21Cv1LdAHRCRdqq6T0Ta4TeiANgLnOPK18FK2wvkBqUXWukdQuSPp42E0opDCam3uLiY4cOHNwL3\nRPcD02vIMwj/kmPiERFHdHw+X1jxsfOmp6dTVVUV8d8h2Mt4QjyQr1sHw4Y5l4G+06OkTx94913C\nxp83GOqZ+v4mrsJvlYf19yVX+gTxMxg4rKr7gNeAkSLSWkRaAyOB16x7X4vIYPHvhJwQVFcsbSSU\nZpxMSL3z589vwOLkngkLJ07iOupHnM4880xEhKqqKnbt2kV5eTnf/e53ycjIwBP00vb5fEydOpX7\n778/QMDsDbfu/MEhMuzwGV6vt9q9iKxbF3qjqkiAOEWkZ0+orKS4qIh5c+dSXFR0asJu82YjToaU\nImEjKBF5Gv/o5ywR2QP8GrgHeE5Efg58ClxlZX8Fv5XdTuA4cB2Aqh4UkTnABivf3ap60Dq/EViK\nPyLgq9ZBrG00RBYtWsTKlSuT3Y0YSf7UXU0cPHiwWtrJkydZu3YthYWFZGVl8d577wE4AQZnzpwZ\nEEOqf//+TtyoZcuWAdC3b1/H84Pt2DWsB/J//AO+851aPYeK0CIjg9Ly8pBTiCbIoaHBUNtFrMZ+\n1NaKbw63+6/rAHth/Vvf+lYKGDDUdNwRxRr7zBToZ+QjLy8v7L9FZmZmNUvJ4NDtIc3AX3qpdoYJ\n9lFZWeN3JZQRRrR7pwyG2kBDNpJoKrSvo2WuRYsWceONN1JZWVkn9SWG1B8lxYLX66WgoCDkPdsz\nvLqm9zIyMk6Fbn/2Wfjxj8nBP1wH4MSJyIGlQlFZGde0W6QQHJEi9BoMqYQRqATzbT6udR3FxcUp\nLE6pKUotWrTg2LFjYQ0cgrGNI+xfbiLC5MmTw4a5yMrK4qdVVTzprqS0NHYBgrhFKF5MkENDQ8EI\nVILJoqTWdRQWFqaQOI3jlN1JODYAA+uhL+E5evRo1Hk9Hg/jxo1jzJgxzJw50wkg2LdvX3bNmEHO\nAw9Us4ibYh3Rcs9vf8vwiy9OGTEwQQ4NDQGJ9hdmU6V///66cePGmMqISMC4QiDqX/I27l/t77//\nPlOnTo2pfN2SmqOkumIW8P/bu//Yusr7juPvj+34ukmBhMRLQjCB0GiVu6QxpaPu5sgshaaUFqZU\ngwihlkIt0nRL2xUEQ9syadKWsI2CGgFu1Y2NtNCxbqFQ2hJIIGN3KT/aEPAICSFdk5QmWKJpS50R\n+7s/nsf2tX39+5xzT/D3JV35/LrnfP1c3/P1c85znudvEtjPKXV1/ODxx/3E7xwg6VkzO38y+/Aa\nVA71trI6fvw41dXVQ8YUysbbKyndQOikcbL++8kn+UBLCxA+p/Xr1/Poo49iZlQl9UyTcw7wBJVL\n27dvp6urCzPLsEeIcwkt8EeT36T0D8AXEtjPNOAEodHDtm3bAMrer2lubmbVqlV9Q8f39PQwe3b5\ngRCdc+PnT+Xl0BtvvDHuS4ITU9pCerjk9A4GPjRbWbcyfJvw8SSnWsJvUzttGmuuv572u++mdto0\nREhOANdcEx6VG6kxQWdnZ99DuVVVVXR2Tv6eo3Mu8BpUzrS3t3PbbbeleIT8X7p7AFiVwH7etWAB\nrxwq38y/urqa6z/zmb7h2XutXbuWnp4eCoUCTU1Noz7Q2traSqFQ8CbbzqXAE1RO9Hb8mnwPEdOB\nX49hu2yT0vcJ/VZNVh1wfNCy5cuXUywWeWuY5FRTU8OmTZtoaxvYDq+trY0lS5b01ZjG0meeN9l2\nLj2eoHKgWCzS0tKSYFPysdSSzib0BJWel4HFCexnOvCbEdb3Pr8EoWbU2NjIjh07ym57+eWXc+ON\nNw6bSAY3vx7LA63eZNu5dHiCqrBiscjq1asTSE6VuXR3mDAq5GTNoKTHhXFqaWlhx44dA+7b1dbW\ncvx4f91KEjfccAMbNmwY8369duRcZflzUKNI8zmo9vZ21q5dO6Cz0fHJJikdAs6Y9F7gNOBYAvsB\nqK+vp76+nnXr1rFkyZIh94ogdNb62muvMW/evCH3mpxz6fLnoE5ixWKRNWvWjLMZ+d8DXxxlmxXA\n4+OO5wCwcNzvGmoWJD4C1vLly3nqqaf6apmFQoEtW7YMSDjlajqekJw7uXmCqpCNGzeOMTklV0va\nC7xrTFuObDYwdGCKZEiiqqqK7u5uqqqquPPOO2lra6NYLPYNX1GuNuT3gZx7+/EEVQHFYpGRLxtO\nPCntApZOJKhB5kACvQiWV1tbS1VVFUuXLmXmzJksW7aMY8fCxb/e3sDL1YY8ATk3tXiCylhvN0Zd\nXaUj7S4ADo7yzv5Ld48DFyYQSz3wegL7KSWJQqEAwPTp02loaKBQKLB48WKOHj3KqlWrhjTvLseT\nkXPOE1TGep+tMZtHaANX3gq2spWLJn28NJIQwMUXX0xHRweLFi3iqquuorOz01u6OecSNeUSlKSV\nwO1ANfA1M/vbLI/f2tpKd3dotVfPEY4wd9L7TDIJVVVVUVdXR21tLUuXLqWhoYEnnnjCE5FzLnNT\nKkFJqgY2ARcRrqk9LelBM+tI65hXD5pvbm5ma9WHWNHz2Lj2cwow9hGOhurtEX3GjBl0d3djZn33\ngOrr68d1+c0557IwpRIUYRS9fWa2H0DSfcBlQKIJqrq6mh91d9MEbI7zpVZcORe+MfA9ncA8+jsq\nHU2hUKC7u7vvGaqamhrq6+tZvHgxjY2NNDU1eW3HOXdSm2oJagHw05L5g8AFgzeS1Ddg6llnnTXu\ng5w4cYKamhq6u7uprq4e+iDu5s3hRf/AhLNnz+bSRx7h8OHDtLa2cuzYMTo6Oti7dy9dXV0sXLiQ\nU089la6uLq699lqv6Tjn3vamVE8Skj4BrDSz6+L81cAFZva54d4zkZ4knHNuqkuiJ4mpNh7UIaCh\nZP7MuMw551zOTLUE9TSwWNI5kmqBK4EHKxyTc865MqbUPSgzOyHpc4ThiKqBr5vZixUOyznnXBlT\nKkEBmNl3ge9WOg7nnHMjm2qX+Jxzzp0kPEE555zLJU9QzjnncmlKPQc1EZKOAj+Z4NvnkE5frWny\nmLPhMWfDY85GuZgXmln9ZHbqCSpFkp6Z7INqWfOYs+ExZ8NjzkZaMfslPuecc7nkCco551wueYJK\nV3ulA5gAjzkbHnM2POZspBKz34NyzjmXS16Dcs45l0ueoJxzzuWSJ6iUSFopaY+kfZJuqmAcDZK2\nSeqQ9KKkdXH56ZIelbQ3/pwVl0vSHTHu5yWdV7KvT8bt90r6ZAaxV0v6kaSH4vw5knbG2O6PPdIj\nqRDn98X1Z5fs4+a4fI+kD6cc70xJD0h6SdL/SGrOezlL+kL8u3hB0jcl1eWxnCV9XdIRSS+ULEus\nbCW9T9Lu+J47JCmFeG+NfxvPS/p3STNL1pUtv+HOI8N9RknHXLLuTyWZpDlxPpsyNjN/Jfwi9JT+\nCrAIqAV2AY0VimU+cF6cPgV4GWgENgI3xeU3ARvi9CXAI4CADwA74/LTgf3x56w4PSvl2L8IfAN4\nKM5/C7gyTt8FrInTnwXuitNXAvfH6cZY9gXgnPiZVKcY7z3AdXG6FpiZ53ImjDD9KvCOkvL9VB7L\nGVgOnAe8ULIssbIFfhi3VXzvR1KI92KgJk5vKIm3bPkxwnlkuM8o6Zjj8gbCCBA/AeZkWcapnVym\n8gtoBr5fMn8zcHOl44qxbAEuAvYA8+Oy+cCeOH03sLpk+z1x/Wrg7pLlA7ZLIc4zgceAPwAein/U\nr5d8wfvKOH55muN0TdxOg8u9dLsU4j2NcLLXoOW5LWdCgvppPJnUxHL+cF7LGTibgSf8RMo2rnup\nZPmA7ZKKd9C6PwQ2x+my5ccw55GRvgtpxAw8ALwXOEB/gsqkjP0SXzp6v/i9DsZlFRUvyTQBO4G5\nZvazuOo1YG6cHi72rH+nLwM3Aj1xfjbwhpmdKHP8vtji+l/E7bOM+RzgKPCPCpclvyZpBjkuZzM7\nBPwd8L/Azwjl9iz5LudSSZXtgjg9eHmaPk2oRTBKXOWWj/RdSJSky4BDZrZr0KpMytgT1BQh6Z3A\nvwGfN7Njpess/EuTm+cNJF0KHDGzZysdyzjUEC6P3GlmTcCvCZed+uSwnGcBlxGS6xnADGBlRYOa\noLyV7Ugk3QKcADZXOpaRSJoO/BnwF5WKwRNUOg4Rrtv2OjMuqwhJ0wjJabOZfTsu/rmk+XH9fOBI\nXD5c7Fn+Tr8HfFzSAeA+wmW+24GZknoH2Sw9fl9scf1pQGfGMR8EDprZzjj/ACFh5bmcPwS8amZH\nzewt4NuEss9zOZdKqmwPxenByxMn6VPApcBVMalOJN5Ohv+MknQu4Z+XXfG7eCbwnKR5E4h5YmWc\n9HVif/Vdn98fP9zem5vvqVAsAv4Z+PKg5bcy8Abzxjj9UQbe/PxhXH464R7LrPh6FTg9g/hb6W8k\n8a8MvDH82Ti9loE3778Vp9/DwJvP+0m3kcQO4Lfj9PpYxrktZ+AC4EVgeozjHuCP81rODL0HlVjZ\nMvQG/iUpxLsS6ADqB21XtvwY4Twy3GeUdMyD1h2g/x5UJmWc2ollqr8IrVxeJrTCuaWCcfw+4dLH\n88CP4+sSwnXsx4C9wNaSPyIBm2Lcu4HzS/b1aWBffF2TUfyt9CeoRfGPfF/8ghbi8ro4vy+uX1Ty\n/lvi77KHSbbMGkOsy4BnYln/R/yC5rqcgb8CXgJeAP4lniRzV87ANwn3yd4i1FavTbJsgfNjGbwC\nfIVBjV0Sincf4f5M7/fwrtHKj2HOI8N9RknHPGj9AfoTVCZl7F0dOeecyyW/B+Wccy6XPEE555zL\nJU9QzjnncskTlHPOuVzyBOWccy6XPEE555zLJU9QzlVYHMZiq6QfS7pC0udjNzMjvedAydAH/zXK\ntudLuiPJmJ3LQs3omzjnUtYEYGbLICQf4F7gzbG82cw+OMr6ZwgPEDt3UvEalHMpkDRD0sOSdikM\nBnhFHHzuJUnPxQHbHpL0W4Rk9P5Yg1pH6Lh1m6RtYzzWr+LP+yR9tGT5P0n6hKRW9Q/6uD4OTLdd\n0n5Jf1Ky/Z/HwfH+U2Hwwi8lWSbOjZcnKOfSsRI4bGbvNbPfAb4HfBX4GPA+YB6AmR0BrgN2mNky\nM7sdOAxcaGYXjvOY9wN/BBBHWF0BPFxmu3cTxn36XeAvJU2T9H5gFWHcn48QuqVxrqI8QTmXjt3A\nRZI2SGohdPj5qpnttdC/2L0pHPMR4EJJBUKSedLMflNmu4fN7LiZvU7oAXwuoRfzLWbWZWa/BL6T\nQnzOjYsnKOdSYGYvE4bb2A38NfDxDI7ZBWwn1I6uINSoyjleMt2N34t2OeUJyrkUSDoDeNPM7iUM\nC/FB4GxJ58ZNVo/w9l8Cp0zw0PcD1wAthMuKY/UU8DFJdXFwy0sneHznEuP/OTmXjiXArZJ6CMMX\nrAHmAA9LepMwdtRwSagd+J6kwxO4D/UDwrAZW8zs/8b6JjN7WtKDhKFCfk6o+f1inMd2LlE+3IZz\nFSCpFfiSmeWmpiLpnWb2q/gM1pNAm5k9V+m43NTlNSjnXK92SY2EgQnv8eTkKs1rUM7lmKSdhFFu\nS11tZrsrEY9zWfIE5ZxzLpe8FZ9zzrlc8gTlnHMulzxBOeecyyVPUM4553Lp/wEv722hQHV+RQAA\nAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10bc6ac50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(simple_feature_matrix,output,'k.' )\n",
    "plt.plot(simple_feature_matrix,predict_output(simple_feature_matrix, simple_weights_0_penalty),'b-' , \n",
    "         label= 'no regularization')\n",
    "plt.plot(simple_feature_matrix,predict_output(simple_feature_matrix, simple_weights_high_penalty),'r-', \n",
    "         label= 'high regularization')\n",
    "plt.xlabel('sqft_living')\n",
    "plt.ylabel('price')\n",
    "plt.title('price vs. sqft_living')\n",
    "plt.legend(loc= 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the RSS on the TEST data for the following three sets of weights:\n",
    "1. The initial weights (all zeros)\n",
    "2. The weights learned with no regularization\n",
    "3. The weights learned with high regularization\n",
    "\n",
    "Which weights perform best?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([ 0.,  0.]), 1784273286136298.0)\n",
      "(array([ -1.63113515e-01,   2.63024369e+02]), 275723632153607.44)\n",
      "(array([   9.76730382,  124.57217567]), 694642101489903.0)\n"
     ]
    }
   ],
   "source": [
    "initial_weights = np.array([0., 0.])\n",
    "\n",
    "def comp_RSS(weights, data= simple_test_feature_matrix):\n",
    "    errors =  np.dot(data, weights ) - test_output\n",
    "    return np.dot( errors, errors)\n",
    "    \n",
    "print( initial_weights             , comp_RSS(initial_weights) )\n",
    "print( simple_weights_0_penalty    , comp_RSS(simple_weights_0_penalty) )\n",
    "print( simple_weights_high_penalty , comp_RSS(simple_weights_high_penalty) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.251461373342378"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(np.dot(test_output, test_output))/np.log(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "('263.0', '124.6')\n"
     ]
    }
   ],
   "source": [
    "print (\"{:.1f}\".format(list(simple_weights_0_penalty)[1]), \"{:.1f}\".format(list(simple_weights_high_penalty)[1] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/qiaolinchen/Library/Python/2.7/lib/python/site-packages/scipy/linalg/basic.py:884: RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver.\n",
      "  warnings.warn(mesg, RuntimeWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([   0.        ,  281.95883963])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simp_linear = LinearRegression().fit(simple_feature_matrix, output)\n",
    "simp_linear.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([   0.        ,  281.95883963])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simp_ridge = Ridge(alpha= 0).fit(simple_feature_matrix, output )\n",
    "simp_ridge.coef_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***QUIZ QUESTIONS***\n",
    "1. What is the value of the coefficient for `sqft_living` that you learned with no regularization, rounded to 1 decimal place?  What about the one with high regularization?\n",
    "2. Comparing the lines you fit with the with no regularization versus high regularization, which one is steeper?\n",
    "3. What are the RSS on the test data for each of the set of weights above (initial, no regularization, high regularization)? \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Running a multiple regression with L2 penalty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us now consider a model with 2 features: `['sqft_living', 'sqft_living15']`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, create Numpy versions of your training and test data with these two features. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_features = ['sqft_living', 'sqft_living15'] # sqft_living15 is the average squarefeet for the nearest 15 neighbors. \n",
    "my_output = 'price'\n",
    "(feature_matrix, output) = get_numpy_data(train_data, model_features, my_output)\n",
    "(test_feature_matrix, test_output) = get_numpy_data(test_data, model_features, my_output)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to re-inialize the weights, since we have one extra parameter. Let us also set the step size and maximum number of iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "initial_weights = np.array([0.0,0.0,0.0])\n",
    "step_size = 1e-12\n",
    "max_iterations = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's consider no regularization.  Set the `l2_penalty` to `0.0` and run your ridge regression algorithm to learn the weights of your model.  Call your weights:\n",
    "\n",
    "`multiple_weights_0_penalty`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting gradient descent with l2_penalty = 0\n",
      "Iteration = 100\n",
      "Cost function =  1.22575873839e+15\n",
      "Iteration = 200\n",
      "Cost function =  1.21173888029e+15\n",
      "Iteration = 300\n",
      "Cost function =  1.20747307976e+15\n",
      "Iteration = 400\n",
      "Cost function =  1.20617512454e+15\n",
      "Iteration = 500\n",
      "Cost function =  1.205780189e+15\n",
      "Iteration = 600\n",
      "Cost function =  1.20566001323e+15\n",
      "Iteration = 700\n",
      "Cost function =  1.20562343801e+15\n",
      "Iteration = 800\n",
      "Cost function =  1.20561229974e+15\n",
      "Iteration = 900\n",
      "Cost function =  1.20560890112e+15\n",
      "Iteration = 1000\n",
      "Cost function =  1.20560785742e+15\n",
      "Done with gradient descent at iteration  1000\n",
      "Learned weights =  [  -0.35743483  243.05416982   22.41481497]\n"
     ]
    }
   ],
   "source": [
    "l2_penalty=0\n",
    "multiple_weights_0_penalty = ridge_regression_gradient_descent(feature_matrix, output, initial_weights, \n",
    "                                                               step_size,  l2_penalty, max_iterations )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's consider high regularization.  Set the `l2_penalty` to `1e11` and run your ridge regression algorithm to learn the weights of your model.  Call your weights:\n",
    "\n",
    "`multiple_weights_high_penalty`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting gradient descent with l2_penalty = 1e+11\n",
      "Iteration = 100\n",
      "Cost function =  3.60529037346e+15\n",
      "Iteration = 200\n",
      "Cost function =  3.60528585733e+15\n",
      "Iteration = 300\n",
      "Cost function =  3.60528134121e+15\n",
      "Iteration = 400\n",
      "Cost function =  3.60527682511e+15\n",
      "Iteration = 500\n",
      "Cost function =  3.60527230902e+15\n",
      "Iteration = 600\n",
      "Cost function =  3.60526779294e+15\n",
      "Iteration = 700\n",
      "Cost function =  3.60526327688e+15\n",
      "Iteration = 800\n",
      "Cost function =  3.60525876083e+15\n",
      "Iteration = 900\n",
      "Cost function =  3.6052542448e+15\n",
      "Iteration = 1000\n",
      "Cost function =  3.60524972878e+15\n",
      "Done with gradient descent at iteration  1000\n",
      "Learned weights =  [  6.74296579  91.48927365  78.43658766]\n"
     ]
    }
   ],
   "source": [
    "l2_penalty = 1e11 \n",
    "multiple_weights_high_penalty = ridge_regression_gradient_descent(feature_matrix, output, initial_weights, \n",
    "                                                               step_size,  l2_penalty, max_iterations )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the RSS on the TEST data for the following three sets of weights:\n",
    "1. The initial weights (all zeros)\n",
    "2. The weights learned with no regularization\n",
    "3. The weights learned with high regularization\n",
    "\n",
    "Which weights perform best?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((17384,), (4229,))"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(feature_matrix, initial_weights ).shape, test_output.shape\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(array([ 0.,  0.,  0.]), 1784273286136298.0)\n",
      "(array([  -0.35743483,  243.05416982,   22.41481497]), 274067615918575.56)\n",
      "(array([  6.74296579,  91.48927365,  78.43658766]), 500404800500841.25)\n"
     ]
    }
   ],
   "source": [
    "def comp_RSS(weights, data ):\n",
    "    errors =  np.dot(data, weights ) - test_output\n",
    "    return np.dot( errors, errors)\n",
    "\n",
    "print( initial_weights               , comp_RSS(initial_weights            , data= test_feature_matrix) )\n",
    "print( multiple_weights_0_penalty    , comp_RSS(multiple_weights_0_penalty   , data= test_feature_matrix) )\n",
    "print( multiple_weights_high_penalty , comp_RSS(multiple_weights_high_penalty, data= test_feature_matrix) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.437857721935876"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.log(274067615918575.56)/ np.log(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "243.1\n",
      "91.5\n"
     ]
    }
   ],
   "source": [
    "print(\"{:.1f}\".format( list(multiple_weights_0_penalty)[1] ))\n",
    "print(\"{:.1f}\".format(list(multiple_weights_high_penalty)[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict the house price for the 1st house in the test set using the no regularization and high regularization models. (Remember that python starts indexing from 0.) How far is the prediction from the actual price?  Which weights perform best for the 1st house?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "387465.47605823533"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(test_feature_matrix.iloc[0], multiple_weights_0_penalty )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "270453.53032194078"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.dot(test_feature_matrix.iloc[0], multiple_weights_high_penalty )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "310000.0"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "test_output[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "***QUIZ QUESTIONS***\n",
    "1. What is the value of the coefficient for `sqft_living` that you learned with no regularization, rounded to 1 decimal place?  What about the one with high regularization?\n",
    "2. What are the RSS on the test data for each of the set of weights above (initial, no regularization, high regularization)? \n",
    "3. We make prediction for the first house in the test set using two sets of weights (no regularization vs high regularization). Which weights make better prediction <u>for that particular house</u>?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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
}