Course_2/Week_4_2_Ridge_regression.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 4_2: 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": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[WARNING] Unable to write current GraphLab Create license to /Users/jcj/.graphlab/config. Ensure that this user account has write permission to /Users/jcj/.graphlab/config to save the license for offline use.\n",
      "[INFO] This non-commercial license of GraphLab Create is assigned to chengjun@chem.ku.dk and will expire on January 27, 2017. For commercial licensing options, visit https://dato.com/buy/.\n",
      "\n",
      "[INFO] Start server at: ipc:///tmp/graphlab_server-95655 - Server binary: /usr/local/lib/python2.7/site-packages/graphlab/unity_server - Server log: /tmp/graphlab_server_1455271694.log\n",
      "[INFO] GraphLab Server Version: 1.8.1\n",
      "[WARNING] Unable to create session in specified location: '/Users/jcj/.graphlab/artifacts'. Using: '/var/tmp/graphlab-jcj/95655/tmp_session_df228808-b92a-462c-8f69-4b29f60a8cd2'\n"
     ]
    }
   ],
   "source": [
    "import graphlab\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sales = graphlab.SFrame('kc_house_data.gl/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import useful functions from previous notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_numpy_data(data_sframe, features, output):\n",
    "    \"\"\"\n",
    "    features is a list of features.\n",
    "    output is a list or a string.\n",
    "    \"\"\"\n",
    "    data_sframe['constant'] = 1\n",
    "    features = data_sframe[['constant']+features]\n",
    "    features_matrix = features.to_numpy()\n",
    "    \n",
    "    output_array = data_sframe[output].to_numpy()\n",
    "    return (features_matrix, output_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def feature_derivative_ridge(errors, feature, weight, l2_penalty, feature_is_constant):\n",
    "    if feature_is_constant:\n",
    "        feature_derivative = 2 * np.dot(errors, feature)\n",
    "    else:\n",
    "        feature_derivative = 2 * np.dot(errors, feature) + 2 * l2_penalty * weight\n",
    "    return (feature_derivative)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-5.65541667824e+13\n",
      "-5.65541667824e+13\n",
      "\n",
      "-22446749336.0\n",
      "-22446749336.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[:,1], my_weights[1], 1, False)\n",
    "print np.sum(errors*example_features[:,1])*2+20.\n",
    "print ''\n",
    "\n",
    "# next two lines should print the same values\n",
    "print feature_derivative_ridge(errors, example_features[:,0], my_weights[0], 1, True)\n",
    "print np.sum(errors)*2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ridge_regression_gradient_descent(feature_matrix, output, initial_weights, step_size, l2_penalty, max_iterations=100):\n",
    "    weights = np.array(initial_weights) \n",
    "    n_iteration = 0\n",
    "    while n_iteration < max_iterations:\n",
    "        # compute the predictions using your predict_output() function\n",
    "        predictions = predict_output(feature_matrix, weights)\n",
    "\n",
    "        # compute the errors as predictions - output\n",
    "        #errors = output - predictions\n",
    "        errors = predictions - output\n",
    "\n",
    "        for i in xrange(len(weights)):\n",
    "        # Recall that feature_matrix[:,i] is the feature column associated with weights[i]\n",
    "        # compute the derivative for weight[i].\n",
    "            if i == 0:\n",
    "                derivative = feature_derivative_ridge(errors, feature_matrix[:, i], weights[i], l2_penalty, True)\n",
    "            else:\n",
    "                derivative = feature_derivative_ridge(errors, feature_matrix[:, i], weights[i], l2_penalty, False)\n",
    "            # subtract the step size times the derivative from the current weight\n",
    "            #print derivative\n",
    "            weights[i] = weights[i] - step_size * derivative\n",
    "            \n",
    "        n_iteration += 1\n",
    "    return weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing effect of L2 penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "train_data,test_data = sales.random_split(.8,seed=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "simple_features = ['sqft_living']\n",
    "my_output = 'price'\n",
    "(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": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ -1.63113501e-01   2.63024369e+02]\n"
     ]
    }
   ],
   "source": [
    "simple_weights_0_penalty = ridge_regression_gradient_descent(simple_feature_matrix, \n",
    "                                                             output, \n",
    "                                                             initial_weights=[0.0, 0.0], \n",
    "                                                             step_size=1e-12, \n",
    "                                                             l2_penalty=0.0, \n",
    "                                                             max_iterations=1000)\n",
    "print simple_weights_0_penalty\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   9.76730383  124.57217565]\n"
     ]
    }
   ],
   "source": [
    "simple_weights_high_penalty = ridge_regression_gradient_descent(simple_feature_matrix, \n",
    "                                                             output, \n",
    "                                                             initial_weights=[0., 0.], \n",
    "                                                             step_size=1e-12, \n",
    "                                                             l2_penalty=1e11, \n",
    "                                                             max_iterations=1000)\n",
    "print simple_weights_high_penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10f8c6e10>,\n",
       " <matplotlib.lines.Line2D at 0x10f8c6f10>,\n",
       " <matplotlib.lines.Line2D at 0x10f8d9150>,\n",
       " <matplotlib.lines.Line2D at 0x10f8d9d10>,\n",
       " <matplotlib.lines.Line2D at 0x10f8d9e50>,\n",
       " <matplotlib.lines.Line2D at 0x10f8e79d0>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEACAYAAABoJ6s/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvX18XVWZ9/290iRtgDZNIBRMRYq8P61C25Q2+jhJeE/u\nD6KfEapDWxDtKDy3ODoKRaXtPThaRkecGamTkUdoRykIMwPeCW3ENvN4C9UUii0U2qoI5CClkDb1\npaVv1/PHXud0n529z0uS85ZzfT+f/ek611577Wunyf6dtda1riWqimEYhmHkg4pCO2AYhmGUDyY6\nhmEYRt4w0TEMwzDyhomOYRiGkTdMdAzDMIy8YaJjGIZh5I2MREdElojI8yKyRUR+ICLVIlInIj0i\nsl1E1olIbaD+ThF5QUQu89lnujZ2iMjdPnu1iKxx1zwlIqf5zi1y9beLyEKf/XQR2ejOPSAilSP/\ncRiGYRi5JK3oiMi7gE8CF6rqe4BK4KPAbcATqnoOsB5Y4uqfD1wDnAdcCdwjIuKaWwncqKpnA2eL\nyOXOfiMwoKpnAXcDd7m26oA7gCbgImCpT9xWAN90be11bRiGYRhFTCY9nX3AQeB415uoAWLAB4H7\nXZ37gatd+SpgjaoeVtXfATuBOSJyCjBRVftcvVW+a/xtPQy0ufLlQI+qDqrqXqAHuMKdawMe8d3/\nQxk9sWEYhlEw0oqOqu4Bvgm8gic2g6r6BDBFVXe5Oq8DJ7tLGoFXfU3EnK0R6PfZ+50t6RpVPQIM\nikh9VFsiciKwR1WP+tp6RyYPbBiGYRSOTIbXzgD+BngX3ov9eBH5KyCYP2c08+lI+ioZ1TEMwzCK\niEwm32cDP1fVAQAR+U+gGdglIlNUdZcbOnvD1Y8B7/RdP9XZouz+a14TkXHAJFUdEJEY0BK4ZoOq\nviUitSJS4Xo7/raSEBFLLmcYhjEMVHXUv9xnMqezHZgrIhNcQMDFwDbgMeB6V2cR8KgrPwbMdxFp\n04AzgV+6IbhBEZnj2lkYuGaRK38ELzABYB1wqROYOuBSZwPY4OoG7z8EVS3ZY+nSpQX3oRx9N/8L\nf5j/hT1yRdqejqr+SkRWAU8DR4DNQCcwEXhIRD4OvIwXsYaqbhORh/CE6RBwkx57gpuB+4AJQLeq\nrnX2e4HVIrITeAuY79raIyJ/B2zCG75brl5AAXjRc2vc+c2uDcMwDKOIyWhti6r+A/APAfMAcElE\n/a8BXwuxPw3MCLG/jROtkHP34QlV0P4SXhi1YRiGUSJYRoIip6WlpdAuDJtS9h3M/0Jj/o9NJJdj\nd8WAiOhYf0bDMIzRRkTQAgUSlDyxWGhgm2EYRk6IxWJ0dHTQ0dFh758AZdHTaW9vp6urq9CuGIZR\nJnR0dNDd3Q1Aqb5/rKdjGIZhlDxl0dPp7++nsbExfWXDMIxRIBaLsXjxYgA6OztL8v2Tq55OWYjO\nWH9GwzCM0caG1wzDMIySx0THMAzDyBsmOoZhGEbeMNExDMMw8oaJjmEYhpE3THQMwzCMvGGiYxiG\nYeQNEx3DMAwjb5joGIZhGHnDRMcwDMPIGyY6hmEYRt4w0TEMwzDyRlrREZGzRWSziDzj/h0Ukc+I\nSJ2I9IjIdhFZJyK1vmuWiMhOEXlBRC7z2WeKyBYR2SEid/vs1SKyxl3zlIic5ju3yNXfLiILffbT\nRWSjO/eAiFSOzo/EMAzDyBVpRUdVd6jqhao6E5gF/An4T+A24AlVPQdYDywBEJHzgWuA84ArgXtE\nJJ6pdCVwo6qeDZwtIpc7+43AgKqeBdwN3OXaqgPuAJqAi4ClPnFbAXzTtbXXtRGK7dxnGKWF7bw5\ndsl2eO0S4Deq+irwQeB+Z78fuNqVrwLWqOphVf0dsBOYIyKnABNVtc/VW+W7xt/Ww0CbK18O9Kjq\noKruBXqAK9y5NuAR3/0/FOV0fF8LwzBKg8WLF9Pd3U13d7f9/Y4xshWda4EfuvIUVd0FoKqvAyc7\neyPwqu+amLM1Av0+e7+zJV2jqkeAQRGpj2pLRE4E9qjqUV9b78jyWQzDMIw8k/E8iIhU4fVibnWm\n4M5oo7lTWiYbB2W8uVBnZ+cIXDEMI990dnYm7bxpjB2ymXy/EnhaVd90n3eJyBRV3eWGzt5w9hjw\nTt91U50tyu6/5jURGQdMUtUBEYkBLYFrNqjqWyJSKyIVrrfjb2sI//Zv/5Yot7S00NLSElXVMIwi\noLGxka6urkK7UVb09vbS29ub8/tkvF21iDwArFXV+93nFXiT/ytE5FagTlVvc4EEP8Cb+G8EfgKc\npaoqIhuBzwB9QBfwT6q6VkRuAqar6k0iMh+4WlXnu0CCTcBMvKHATcAsVd0rIg8C/6GqD4rISuBX\nqvrdEL9tu2rDMIwsydV21RmJjogcB7wMnKGqf3C2euAhvB7Ky8A1brIfEVmCF012CLhFVXucfRZw\nHzAB6FbVW5x9PLAauBB4C5jvghAQkeuBL+EN392pqqucfRqwBqgDNgPXqeqhEN9NdAzDMLKkoKJT\nypjoGIZhZE+uRMcyEhiGYRh5w0THMAzDyBsmOoZhGEbeMNExDMMw8oaJjmEYhpE3THQMwzCMvGGi\nYxiGYeQNEx3DMAwjb5joGIZhGHnDRMcwDMPIGyY6hmEYRt4w0TEMwzDyhomOYRiGkTdMdAzDMIy8\nYaJjGIZh5A0THcMwDCNvmOgYhmEYecNExzAMw8gbJjqGYRhG3shIdESkVkR+JCIviMjzInKRiNSJ\nSI+IbBeRdSJS66u/RER2uvqX+ewzRWSLiOwQkbt99moRWeOueUpETvOdW+TqbxeRhT776SKy0Z17\nQEQqR/7jMAzDMHJJpj2dbwPdqnoe8F7gReA24AlVPQdYDywBEJHzgWuA84ArgXtERFw7K4EbVfVs\n4GwRudzZbwQGVPUs4G7gLtdWHXAH0ARcBCz1idsK4Juurb2uDcMwDKOISSs6IjIJ+L9V9fsAqnpY\nVQeBDwL3u2r3A1e78lXAGlfvd8BOYI6InAJMVNU+V2+V7xp/Ww8Dba58OdCjqoOquhfoAa5w59qA\nR3z3/1DGT20YhmEUhEx6OtOAN0Xk+yLyjIh0ishxwBRV3QWgqq8DJ7v6jcCrvutjztYI9Pvs/c6W\ndI2qHgEGRaQ+qi0RORHYo6pHfW29I5MHNgzDMApHJvMglcBM4GZV3SQi38IbWtNAveDnkSDpq2RU\nB4Bly5Ylyi0tLbS0tGTvkWEYxhimt7eX3t7enN8nE9HpB15V1U3u8yN4orNLRKao6i43dPaGOx8D\n3um7fqqzRdn917wmIuOASao6ICIxoCVwzQZVfcsFN1S43o6/rSH4RccwDMMYSvAL+fLly3Nyn7TD\na24I7VUROduZLgaeBx4Drne2RcCjrvwYMN9FpE0DzgR+6YbgBkVkjgssWBi4ZpErfwQvMAFgHXCp\nE5g64FJnA9jg6gbvbxiGYRQpopp+VExE3gt8D6gCfgvcAIwDHsLrobwMXOMm+xGRJXjRZIeAW1S1\nx9lnAfcBE/Ci4W5x9vHAauBC4C1gvgtCQESuB76EN3x3p6qucvZpwBqgDtgMXKeqh0J810ye0TAM\nwziGiKCqGU9jZNzuWH8hm+gYhmFkT65ExzISGIZhGHnDRMcwDMPIGyY6hmEYRt4w0TEMwzDyhomO\nYRiGkTfKQnRisch1o4ZhjAKxWIyOjg46Ojrs781ISVmETLe3t9PV1VVoVwxjzNLR0UF3dzcA9vc2\nNrCQacMwDKPkKYueTn9/P42NjekrG4YxLGKxGIsXLwags7PT/t7GAJaRYJhYRgLDMIzsseE1wzAM\no+Qx0TEMwzDyRlmIjoVwGuWEhS8bxUxZiM7ChQsL7YJhZM1wxWPx4sV0d3fT3d2dmNw3jGKhLERn\n69athXbBMLLGxMMYi2SyXXXJM3369EK7YBh5o7OzMyl82TCKibIImbZ1OkYpYmtfjEJi63SGia3T\nMQzDyB5bpzMCLILHKCcses0oZjISHRH5nYj8SkQ2i8gvna1ORHpEZLuIrBORWl/9JSKyU0ReEJHL\nfPaZIrJFRHaIyN0+e7WIrHHXPCUip/nOLXL1t4vIQp/9dBHZ6M49ICKR81M2CWuUExaAYBQzmfZ0\njgItqnqhqs5xttuAJ1T1HGA9sARARM4HrgHOA64E7hGReBdtJXCjqp4NnC0ilzv7jcCAqp4F3A3c\n5dqqA+4AmoCLgKU+cVsBfNO1tde1YRiGYRQxmYqOhNT9IHC/K98PXO3KVwFrVPWwqv4O2AnMEZFT\ngImq2ufqrfJd42/rYaDNlS8HelR1UFX3Aj3AFe5cG/CI7/4finLeIniMcqKzs5P29nba29vtd98o\nOjINmVbgJyJyBPhXVf0eMEVVdwGo6usicrKr2wg85bs25myHgX6fvd/Z49e86to6IiKDIlLvt/vb\nEpETgT2qetTX1juinLeoH6OcaGxstP1sjKIlU9F5n6r+XkQagB4R2Y4nRH5GM0Qsk4iJjKMqli1b\nlii3tLTQ0tKSvUeGYRhjmN7eXnp7e3N+n4xER1V/7/7dLSL/BcwBdonIFFXd5YbO3nDVY8A7fZdP\ndbYou/+a10RkHDBJVQdEJAa0BK7ZoKpviUitiFS43o6/rSH4RccwDMMYSvAL+fLly3Nyn7RzOiJy\nnIic4MrHA5cBW4HHgOtdtUXAo678GDDfRaRNA84EfqmqrwODIjLHBRYsDFyzyJU/gheYALAOuNQJ\nTB1wqbMBbHB1g/c3DMMwipS0i0OdcPwn3vBZJfADVf26m3N5CK+H8jJwjZvsR0SW4EWTHQJuUdUe\nZ58F3AdMALpV9RZnHw+sBi4E3gLmuyAEROR64Evu/neq6iqfX2uAOmAzcJ2qHgrx3xaHGoZhZIll\nJBgmJjqGYRjZYxkJDMMwjJLHRMcwDMPIGyY6hmEYRt4oC9GxpIeGURxYMlKjLAIJ2tvbbYW2YRQB\nHR0ddHd3A2B/l8WNBRIYhmEYJU9Z9HRs51DDKA5sN9TSwdbpDBMTHSPf2IvVGAvY8NoIWLBgQaFd\nMMYoYRPjtomaYURTFqLz3HPPFdoFY4xiAmMY2ZHp1gYlzYwZMwrtglFGdHZ2Jg2vGYZxDJvTMYwR\nYPM3xljFAgmGiSX8NAzDyB4LJDAMwzBKHhMdwzAMI2+UhehYjiejVLFcZcZYoyzmdCzHk1GqWK4y\no1DYnI5hGIZR8pRFT8dCpo1SxUKyjUJR8JBpEakANgH9qnqViNQBDwLvAn4HXKOqg67uEuDjwGHg\nFlXtcfaZwH3ABKBbVT/r7NXAKmAW8CZwraq+4s4tAr4EKPBVVV3l7KcDa4B64GlggaoeDvHbQqYN\nwzCypBiG124Btvk+3wY8oarnAOuBJQAicj5wDXAecCVwj4jEHV8J3KiqZwNni8jlzn4jMKCqZwF3\nA3e5tuqAO4Am4CJgqYjUumtWAN90be11bRiGYRhFTEaiIyJTgXbgez7zB4H7Xfl+4GpXvgpYo6qH\nVfV3wE5gjoicAkxU1T5Xb5XvGn9bDwNtrnw50KOqg6q6F+gBrnDn2oBHfPf/UJT/FvVjFCMWmWaU\nI5n2dL4FfAFviCvOFFXdBaCqrwMnO3sj8KqvXszZGoF+n73f2ZKuUdUjwKCI1Ee1JSInAntU9aiv\nrXdEOW+JGI1ixJKFGuVI2oSfItIB7FLVZ0WkJUXV0Zw4yWQcMeOxxh07drBs2TIAWlpaaGlpGZ5X\nhmEYY5Te3l56e3tzfp9Msky/D7hKRNqBGmCiiKwGXheRKaq6yw2dveHqx4B3+q6f6mxRdv81r4nI\nOGCSqg6ISAxoCVyzQVXfEpFaEalwvR1/W0Po7e21qB+j6LBs1EYxEfxCvnz58pzcJ+3wmqrerqqn\nqeoZwHxgvaouAH4MXO+qLQIedeXHgPkiUi0i04AzgV+6IbhBEZnjAgsWBq5Z5MofwQtMAFgHXOoE\npg641NkANri6wfsbRsYUcl6lsbGRrq4uurq67EuRUTaMZHHo1/EEYTtwsfuMqm4DHsKLdOsGbvLF\nLN8M3AvsAHaq6lpnvxc4SUR2Ap/Fi4xDVfcAf4cXqv0LYLkLKMDV+ZyI7MALm743ylEbLzeisHkV\nw8gvWW3ipqr/Dfy3Kw8Al0TU+xrwtRD708CQHdVU9W28MOuwtu7DW9sTtL+EF0adlj179mRSzRjD\nFNsiy2LzxzDyRVlkJKivr+ett94qtCtGAYnKYVaol7/lVDOKnVwtDi2L7aoNw09fXx+xWIzGxsbE\nvIphGPmhLHo6zc3N/PznPy+0K0YBicViXHjhhezevRsofO/ChteMYqcY0uCULG+//XahXTBGQFiE\nWTZRZ/4XfLFgkWtGuVIWPZ3KykoOHTpUaFeMYRI2/3HxxRezfr0XWd/W1sZPf/rTjK5vaGigqanJ\neheGkQab0xkBhw8PST5tlCjx+ZitW7cmbP5yOpqammwOxzAKSFn0dCZNmsTg4GChXTGGSdh8zP79\n+9mwYQMAra2tiV5P1PU2f2KUO0eOwLhxmde3OZ0RMHv27EK7YIyAxsZGmpqakmyrV6+mvb2d9vZ2\nVq9enfZ6mz8xypH162HaNBCByko4cKDQHpVJT8d2Di19rLdiGOk5cgT+9V/h5puT7TfeCF//Opx0\nUuZtFXzn0FLFdg41DGMsMzgIX/4y/Mu/JNv/4R/gllugqmp47VoggWEYhgHAjh1w003gD9qcMgW+\n9z34H/+jcH5lQlnM6diujOWH7cppjDXWrYOpU735mXPO8QSnpQW2bQNVeP314hccKBPRKbaFgUbu\nsezRRqlz+DDcfbcnMiJwxRUQi8GnPgUDA57QbNgA551XaE+zoyxE56WXXiq0Cwal3fsYie+5eO50\nbZbyz7qcGRiAv/5rT2SqquBv/sazf/vbcOiQJzQrV0JdXWH9HBGqOqYPvG201Sg87e3tGv//aG9v\nz+m9+vv7tb29Xdvb27W/v3/E7Y3E91w8d7o28/mzNkbG88+rfuADqp6keMc736m6dm1h/XLvzVF/\nJ1sggVH0DCdc2rJHG8XMj3/shTG79c4AXHqpF4F29tmF8ysv5ELJiukA9OSTTx6O0BujTLz30dbW\npq2trUN6IVG9k2L41t7X16cNDQ3a0NCgfX19WV0bfC7/576+vmH1yNL15Ea7p2eMjIMHVVesSO7N\ngOpnPqO6d2+hvQuHHPV0Ci4KuT6w4bWiI0pEsrXnk9H0wd9WQ0NDwZ/NyA27d6vecMNQobnnHtXD\nhwvtXXpyJTplEUhg5I/hTGDHr+nr6ws939nZmUh509nZOeL7GUau2LIF5s3zAgEaGuD734czzvDS\n0cRl59Ofzi4H2pgjnSoB44FfAJuB54G/d/Y6oAfYDqwDan3XLAF2Ai8Al/nsM4EtwA7gbp+9Gljj\nrnkKOM13bpGrvx1Y6LOfDmx05x4AKiP8t55OHmlra0t8c29ra1PV1MNL8XL8murq6qyGsPLVCxrN\n4arRGF4zioOjR1Uffli1tja5N9PRofqb3xTau5FBIYfXgOPcv+Pci/59wArgi85+K/B1Vz7fCVSl\nE4Zfcyzdzi+AJlfuBi535U8D97jytcAaPSZsvwFqgcnxsjv3IPARV14J/HWE7yY6ecQ/XNTQ0KCq\n2UVaZSsgxTD0ZpQXBw6o3nnn0GGzz39edd++Qns3euRKdDIaXlPVP7vieLy1PXuADwL3O/v9wNWu\nfJUTjcOq+ju83sscETkFmKiq8TGUVb5r/G09DLS58uVAj6oOqupevJ7VFe5cG/CI7/4fyuRZjNwy\nffr00HIY8aGxAwcO0NraSkNDQ+JcfN+cYN3gMNry5ctpaGigoaGB5cuXp62fD2zIb+zx+utw3XXe\nsNmECV6uM/DSzhw54snON74BEycW1s+SIBNlwhOazcA+4C5n2xOoM+D+/WfgYz7794APA7PwBCRu\nfz/wmCtvBd7hO7cTqAc+D9zus38Z+BxwIrDDZ58KbInw3Xo6eaK/v1/b2tq0oaFBW1tbE8NFmUal\n9ff3R06shw3bhbWRzj7S58tkKGy07m0RaIXl6adVZ81K7s2ce67qz35WaM/yA4Vcp6OqR4ELRWQS\nsE5EWtwfVVK1TNrKkEwym2aV/XTZsmUAtLS00NLSkr1HRloWL16c2EytpqYmsZ4mbM1MLBaLDBwI\nY7g7hY4m8dQ68XKu1wHl+37ljio8+KC3fubPfz5mv/pqLyPAaacVzrd80NvbS29vb+5vlK1KAV8B\n/hYvSGCKs50CvODKtwG3+uqvBS7y13H2+cBKfx09Nm/0hq/Od33XfBe41pXfACpceS7weIS/1tPJ\nIf5v4/7eSLpv+MGw4WBAQdwWp7W1NXGutbU1ce/W1lZtaGjQ5ubmpLU/ma6ryWb9TaY9mNHqoRRq\nvqqcelj796veccfQ+ZklS1T/+MdCe1dYKFQgAXASxybva4D/D7gYL5DgVmcPCySoBqaRHEiwEZiD\n10vpBq5w9ps4Fkgwn/BAgnh5sjv3oE+AVgKfivDfRCeH+F+M8Zd+Ji+rsBdqqpdsWARcMGjBf22m\nL+ywwIco8v0yLtTLf6wHZ8RiqtdckywylZWqq1Z50WiGR65EJ5PhtVOB+0VE8OZ2VqvqT0VkM/CQ\niHwceBm4Bs/LbSLyELANOATc5B4A4GbgPmAC0K2qa539XmC1iOwE3nLCg6ruEZG/Aza5P4Ll6gUU\ngNejWuPOb3ZtGAWkpqYm4yGgzs7OpNQ2UbY4wSG6jo4OdvtziOSBfKfWsVQ+o8cvfgGf/CT4R2Zn\nzIDOTpg7t3B+lSW5ULJiOgA99dRTs1J4I3NG+9v4cCbr40NjqdYCRTGS9DZjlbEwvHb0qOr993s9\nGH+P5pprVEv0kfIOhVynU8oHoDNmzBjOz9wYJVK9xFItFM3HvEm2jIUX8ljlT3/y5mKC8zN33KH6\n5z8X2rvSw0RnBKKDzenklHQv4mCvxF8vKDKjMZ+QS2EY7hxWrilXMXz5ZdWrr04WmeOOU33gAZuf\nGSkmOiY6RctwMg6ERaz5ezsjeXnmciK8WJN1jvXJfz8/+5nqOeckC82sWaqbNhXas7FFrkTH9tMx\nck48QKCvry8x+b97924WL148JHggPnk+nD108oHf3wMHDiTWJRm54+hRL3HmJz6RbL/uOi8LwJQp\nhfHLGCa5ULJiOgAdP378cITeyJBMeyepMg4EGe43d//anba2tpwONRXTkFYx+TIa/P73yT2Z+HHn\nnV7uMyP3YMNrwxedCRMmDOdnXpbk+uUV1n6YLdU8UCqixGqsvZTHIj/9abjQPPywzc8UAhOdEYgO\nNqeTMdn0MEZrN82wTAbx81E9o0zyudXX1ycm+7PJlpBLTPySWbo0XGhWrSq0Z4aJjolOXghb/xI1\nXBW1mj+TrZRTZROI8sd/LlWPxt9e/Kiurh5Wz2m0KacJ/zAOHlSdMSNcaLZtK7R3hh8THROdvJBK\nEIIvyijRSTfEFWxz0qRJQ0QtnrG6vr5e6+rqhgheWC62sHOpxGekL/3487S2tmpbW1tKkY2LsN+3\nchGdV18NF5mGBm9tjVGc5Ep0LHrNSKKxsZGmpqZEduNUdHd3097eniinw5812c++ffsAmDBhQiJK\nzZ+xOngOwMvKNLQc/FxfX8+4ceMSUXO1tbWjlj4n7HnCskH767W1tSV+ZsFUPyOh2KL9uruho2Oo\n/ZOf9FLPGOWLiU6ZkuolFQ8LPnDgAAcOHODw4cOICBdccAGdnZ1J127evBkgqS1/WPHy5cvpcG+f\n/fv3J+5RXV3NwYMHI/0K2/bAf18/EyZMiPw8d+5cli9fnnjR33fffXznO99J+JrJzyMbH9MxYcKE\nnORTK4ZtEL7wBS+EOciPfgR/+Zd5d8coVnLRfSqmAxteCyVsCCxVShr/6vvgEFGqeQr/ufgQVHt7\neyJXWltbW9KWBMFAAP+GcP5z8+bNiwxiSPUcUfM52cy1BNvLZngtV/NIhZgrOnBA9d3vDh8627kz\nLy4YOQQbXjNyTfDbsp/nnnsuMSzl31Y6G4Lf8uOLQBcuXEhfXx8LFixIqj937tzIb+y//vWvE/60\nt7ezefPmlJvGxdm9ezfd3d2J3sBIei0ATU1NaXsV+cgWnSpD92jy0ktwxhlD7e96F7z4oreVs2Gk\nJBdKVkwH1tMJJd3amPr6em1ubk76Jk9IjyUso3N8gebkyZN13LhxCujkyZO1q6trSE8nGFQQbDvK\nZ78/pAmnThV+PZz1QOUW9vwf/xHem7nllkJ7ZmTKcH5nseg1E51c4P9ljK+78b/M4y/j4NYBqdoK\nayPeTlg5Sjyi/AzbxM0/FBe1Jied0BZDNFmxCNqnPx0uND/+ccFcMkbAcH7PTXRMdHJC8JfR/9l/\nVFVV6cyZM7Wurk7r6+sTIcxR4cBhR3CtTPB8XV1dUk/JH0adag4q2FY2iTjjodn+uaNCvvgLJYJ/\n/rPqlCnhQvPyy3lzw8gRJjomOkVD8Jcx/sJtbm7WysrKlCISFCn/yz547eTJk0N7TumGvPxDblF/\nNMH68+bNy0o0Uglvvns/+bz3iy+Gi8z556u+/XZOb23kmWIaXrNAgjInKstzR0cHhw8fzqqtI0eO\n0NDQwIwZM1ixYgVLly5Nukc8SKGpqYlTTz0VgBkzZqCq1NTURE6Ab926lc2bNw+ZKI8HITzzzDNJ\n9bPZNrvYyHVAwAMPwMc+NtR+++3w1a+O+u2MIqGotj7PhZIV04H1dDIiOKTlHyqrrq7W+vp6nT59\nulZWVmpVVZU2NzennaQPth8Vxhys39/fr3V1dYnzwYwDcaKGAv33GE5gQLHMq4wGR4+qLlwY3qNZ\nt67Q3hnFDIUaXgOmAuuB54GtwGecvQ7oAbYD64Ba3zVLgJ3AC8BlPvtMYAuwA7jbZ68G1rhrngJO\n851b5OpvBxb67KcDG925B4DKCP9NdDIgakgr00n3TIaF0iX5jKob9eIPE52ozeHKiYGBcJEB1Vis\n0N4ZpUKuRCeT4bXDwOdU9VkROQF4WkR6gBuAJ1T1LhG51QnNbSJyPnANcJ4TrCdE5Cz3ECuBG1W1\nT0S6ReRQ1PhNAAAb1ElEQVRyVV0H3AgMqOpZInItcBcwX0TqgDucWIm796OqOgisAL6pqj8SkZWu\njX/N4HmMDAhbOe/PKOAvQ/SwkH+l//79+9mwYQPgpadpbW1NOayWjs7OThYuXMizzz6LqnLBBRew\nevXqgqeAKQRr1sBHPzrU3tQETz4JlTaQbhQL2aoU8F/AJcCLwBRnOwV40ZVvA2711X8cuMjV2eaz\nzwdWuvJa4CJXHge8EazjPq8ErnXl3UCFK88F1kb4W1Y9nageQrp9bFJlk47Xq6+vT+oJBc83Nzdr\nXV1dUpaAqEADfMNm8ev92QmGsxVB8HkyWe9TysNnJ5wQ3pu5/PJCe2aMBSiG6DW8Ia3fAScAewLn\nBty//wx8zGf/HvBhYBbQ47O/H3jMlbcC7/Cd2wnUA58HbvfZvwx8DjgR2OGzTwW2RPhcFqKTbm4l\nbLgpKAhRYhQW3hzVdvyorKzU5ubmpDDp+vr6pKg2EdGGhgadNWvWkOv91/lT8KQSiUyH1IYz9FYI\noQre8+jR6GGzO+/Mi0tGGZEr0cm40+2G1h4GblHVP4qIBqoEP48ESV8lozoJli1bBkBLSwstLS3Z\ne1TkRGVwToV/iGz37t0sWLCA9evXJ7XV19eXlJW5oaGBpqamtENihw8f5sknn0yyDQwMUOkb51FV\ndu/ezZtvvjnk+ngy0Orqavbt25cYlov7mG8KkVDTu+d24NdMnTr0/P/+329wzz03AHD99Z1A+Q0r\nGqNHb28vvb29ub9RJsqEl416LZ7gxG0vkDy89oKGD6+t5djw2gs+e6bDa9/1XfNdjg2vvUHy8Nrj\nEb6XRU/H/+29qqpKq6urdd68eYlv5cFdPvv7+7W2tnZI7yTYswmWw9qLp7dpbm4e0mPJ9ogPz/lt\nwUWlQTLd1yZYP5totXS9o9HsCd18c3SP5uDBzH0yjJFAIYfXgFXAPwZsK+LiAtwKfN2Vzwc240Wk\nTQN+DYg7txGYg9dL6QaucPabgHv0mNCsceU64DdAra882Z170CdAK4FPRfheFqITtTo/bCituro6\nKSTZfwRFpq+vb8iKfVVNmt+pqqpKvGz99myPqqoqbWtrGyJ8UXNJcYLDhFEv/ky2uM4k8i6s7ZEK\nQJTIgObsnoaRioKJDvA+4AjwrBOTZ4Ar8OZcnsALZe6Ji4G7ZokTm2DI9Cy8+ZudwLd99vHAQ86+\nETjdd+56Z99Bcsj0NOAXzv4gUBXhf1mITpzg/ErUKv+oI6xnE/Zy8/c+4kdzc3NWolNVVRVq96e+\nic8rhQlf1DNHvYSjXtJh9ky23I5ad5TJHNThw9Ei81d/ldn/9VgJiDCKk4KJTqkf5SY6fX19Wl9f\nnzS8Fl/4GRSKeB61yZMna11dnba2tg7pZUSlhQkbSvMHCVRUVOjs2bNTitDEiRMTWaiDwhckKuhB\nNbyXV11dnRhGjEfl+X3xi85wEoEGt8vOZA3Sxo3RQhPYEsgwCo6JjolOEtkMFfltkydPTgyv+bca\n8LcR3MYgLmSVlZU6efLk0CGwdD2Z6urq0N5R2BGWOy2YTDRqXiU4/xPsBQWH38LCtcM2fgv2JII9\nQj/+a089tTdSaA4fHp3fBcPIBSY6JjpJ+F9sbW1tieEnf2LN+vr6xDf9uG3SpElJPR1/G/EXvb8X\nE38RB4Whvr4+q6G0+vp6bW1t1ebmZnWRj6HH8ccfnyRWcXEJ7p9TXV2dtJ4oLhTpRCcoVmHPFm83\nVWqfoDD7iRKZMfhraIxhTHRMdJJItegy+JL1vyCDvY8wAfILU319fWTQQfxIJSLBI9Pejv+Fnioq\nLkwgRESrq6u1q6srdJsEP0Ex8/dcUg2x+Xuav/1tf6TIfOYzefuVMIxRxUTHRCeJTPexmTRpUpLQ\n+HtCtbW1OnPmzCFCUFFRkZUwjORItX1CdXV12jDssN5MJkIRNWwHQzMlhE3U//Sn4SIDqs8/n5/f\nAcOIZM8e1YceUl20SPWkk479cu7fn3ETuRIdy8g0Brjrrru49dZb2bp1K2eeeSYiwqZNmzh48CD7\n9u1L1BMR3v3ud/P0008DMDg4OGRbAICjR4+mvF9lZWXW2x5Ecdxxx7F//34OHTqUZBcRZs+ezQsv\nvJDy+unTp3PgwIHE56qqqiFtxQlb4FlTU5M4H1z4GkwHf/HFELUu9ehRkKyWKxvGCPnTn7xfyO5u\n73jllei6730v3HgjTJiQP/+iyIWSFdPBGOnpZLM1gGr4N/jROsIiztId6Ybg/FtNp9sMbt68eUmL\nQINzVlFDaWHBAWE7h/p/5lG9mTHwK2WUAm+/rfqTn6h+9rOqZ58d/csI3vnPftarf+DAiG9Njno6\nBReFXB9jRXSCIhMlOvGXsH9ILZs5l0yO4Qy/1dXVJQ3tBQ9/YEBwuCx+LiyJZ3A+J5UQRwUH+O93\n6aV/Gfk3vXRpXv6rjXLj0CHV//N/VG+/XfW9700tLKedpvqpT6k+9pjqH/+YU7dMdEx0kl6WUZmd\nw3o4bW1tkRP4IqIioscff7xOmjRp1AUKvLmj+ALPdHWrqqpCBSootPGouuDanOGIzrnnfi/F3/n5\nGu8ZpSLbxaRGmXH0qOrTT6v+r/+lOnduamE56SRvLuahh7y5mQJholOGopNJqv7gUFRwlX9tba3O\nmzcvr8EBYcI2a9aslEEDUeIT/xwf/vI/a5i4ZpqqJtXffJzgAtBUpBvutJQ1ZcDRo6rbtql+4xuq\nra2pheWEE1SvuUb1vvtUX3+90J6HkivRsUCCIsY/8b1//35EhK1bt7JgwYLQzcrCJvcHBwd56qmn\n8uJvFKqaCF7IlEOHDiEi1NXVsWfPHnbv3s3bb79NdXU1Bw8eTHrW6upqLrnkEjo7O1Nu4NbdHZ0Z\n2vt+kszq1atDN6YzypyXXjo2eZ8qs3tFBXR0eMeVV8Jpp+XPx2ImF0pWTAcl0tNJl4olLJFnf39/\n2sWWxXocd9xxGdULLvYMqxOWBFRV9dVXo79oXnbZmynzuY3G/182540i47XXVO+9V/XDH1adMCF1\nr+WSS1T/8R9Vt28vtNejCjnq6RRcFHJ9lIrohM3ZNDc3a1VVlVZWVg4ZTso2I0Chj1QRb8G1Q2HX\nRj1vMEXNRz4S/W74i7/4RGj0X/xnrmriUFa89ZbqD3+oet11qnV1qYXl/e9X/fu/V332WW8YrQzI\nlejY8FqRsnjx4iGboPkZGBjIozfZEbZW5siRI5H19+7dmyhXV1cPOX/kyJHE88aH1wDq6+uZMWMG\nU6dGD6m1t3fQ1dVFR0dH0hqdIPEN7YJreTo7O5OG2PzDd7FYLPKcUST84Q/wxBPHhsJeey267syZ\nx4bDZs+GcePy52c5kQslK6aDIu3ppNtILGqVfSkc06dPz3j4zH/EM0OnSrtzbIgt+ktp2KR9WFZo\nf8+publ5SDRcMGIu2Kuy4IAiYf9+1bVrVf/n/1Q944zUPZbzz1f9279VXb/eWwNjREKOejoFF4Vc\nH8UqOv7w4ebm5iH5wbyFidm/8GfOnBk6RJXvI2qvnFRH/KXuF4N45FtT040p3iU/TSQUDUb6xcth\nG8H5RSMsUWiqYbjg9TY8l2MOHlTt7VX94hdVp09PLSynn+5tv9rdrfrnPxfa85LFRGcMiI4/+WRU\n4s34S6+vry+rEOP4UVNTU3DBGe6xdOnSgG13infLGZHtxAWgr68vdH2SXyDSLSxNlW16OPvwGCk4\ndEj1O99JLSjx45RTVD/+cdWHH1YdHCy052MSE50xIDpRQ2ZhE+TZZmMeC4cXhRf9ngn7+YVF7sVf\n9kFBCdvKOjh8FtZDyWS9VFh7JjohHD2q+oMfpJ+4jx+trar//u+qu3cX2vOyI1eiY4EEBUREmDt3\nLiLCiy++yL59+xLrT+KT5fF63u/AWMV7tvBH9LJo9vX1sXTp0LPz5s1j8uTJHDhwAFWlpqaGm2++\nmZNPPpk333wzUa+6uprNmzcPmeyPBwrEr9+zZw/vec97ADj33HOZPHkynZ2dicSfwYAEf0JQf3vx\nsp9MAw/GRIDC2rVw003empZMuOgi+Kd/gjlzcuuXUXjSqRJwL7AL2OKz1QE9wHZgHVDrO7cE2Am8\nAFzms88EtgA7gLt99mpgjbvmKeA037lFrv52YKHPfjqw0Z17AKhM4X/R9HT6+vqGzHX4v40Hk10S\n8U2+9I8Ppfhy+0joNfH5Hv8un8HdPuME99aJD1eqZrbjavBIN4+TKZlemyp4oajYuFG1qSmzHguo\nnnWWN+FvlATkqKeTiei8H7ggIDorgC+68q3A1135fGAzUOmE4deAuHO/AJpcuRu43JU/DdzjytcC\na/SYsP0GqAUmx8vu3IPAR1x5JfDXKfwvGtGJ2jAsXq6vr08aVhtbgpPqfXRK2utFRCsrK7W2tjY0\nMag/BY5/uDLVVtJRYhI80s3jZMpwRCdd3ZyzbZvqpZdmLiwnnaT64INls5ZlLEOhRMe7N+8iWXRe\nBKa48inAi658G3Crr97jwEWuzjaffT6w0pXXAhe58jjgjWAdPSYu17rybqDClecCa1P4XhDRib+c\n4t/Igy9DOLZJWVgPZ2wcqd5P6a+fOHFi5LnW1takMGh/mPWsWbO0uro6sXuon1TZudvb23XevHmJ\nhajxMOrR6mlkKljptsoedV55RXX+/MyFRUT1nnu8iX9jzFJsojMQOD/g/v1n4GM++/eADwOzgB6f\n/f3AY668FXiH79xOoB74PHC7z/5l4HPAicAOn32q37cQ3wsiOmHfnE844YTQF2ghk3GO7pE6EGA0\n7xUPL49/DvZuol7YfX192tDQkDTkVoyMeuj1m2+qfvrTmQsLqH71qxZyXMbkSnRGK5BAR6kdiM8c\nj7xOwYjFYvT19Q2x+3e49JNup87i5ibgOxHnuoGOnNz1rLPOYuvWrYnPF1xwARPcrojxBKHgZRrw\nT8zv378/ca69vT00uKAYCO5amhF79nir6bNJ8Pq5z8Edd0BtbXb3MoxhMlzR2SUiU1R1l4icArzh\n7DHgnb56U50tyu6/5jURGQdMUtUBEYkBLYFrNqjqWyJSKyIVqno00FYky5YtA6ClpYWWlpaUdYdL\n/OW2cePGRNqWiooKJk2aRFVVFfv37+ePf/xjTu6dX1J9x6gH9oz6Hf1bZDc0NDB+/PiEeFRXV3Pg\nwAEOHDjAzp07k1LwiEhSaht/mp3du3fznve8hy1btoQKT1FGkf3hDzBpUnbX3HADfPWrcOqpufHJ\nGBP09vbS29ub+xtl0h3CCwrY6vu8Ajd3Q3ggQTUwjeRAgo3AHLxeSjdwhbPfxLFAgvmEBxLEy5Pd\nuQc5Nr+zEvhUCt9zOrwW3+441a6YY+PIz7BZ1JEqPU2qI6xuMEAjas4kH2tuQofRDh7Mbhgsfjz+\neE58NMoTcjS8long/BB4DXgbeAW4wYnAE3ihzD1xMXD1l+CJTTBkehbe/M1O4Ns++3jgIWffCJzu\nO3e9s+8gOWR6Gl403A48AapK4X/8h5cTSjlHWuqjuuBC4xeJrq4unTdvnlZXV2t9fb12dXUNyZXm\nP+ILPeO7qfqjAoO53QoiOkePepFe2QrLnXeOrh+GEQGFEp1SP3IpOsNNVVO8x+Mp3nfrCupbQ0PD\nkJxsfX19SVGCzc3NQ/bHCeZXa2trSwiRP9ddGKMymX/FFdkLywc+MPp+GEaWmOgUoeiMjTDnVO+/\n44vAv2NHWI66KOIv6rDN76LqDvul/sUvZi8skyapHj2a0b0ttY5RCHIlOpYGZ5jEYjEGBwcL7cYw\n0RTnijcwcNasWWzcuDH+ZSIlweABf1qhVHXDUtsAcO+98IlPZO/0gQMwfnzk6WFFqRlGCWOikyWb\nNm2ivb2dgYGBlBuTFRcn4a2njSJ/QlNRUTGsEPH6+nomTJjA3Llz2b59O+PGjUsIRRj+8PTZs2cz\nefJkYGg+ND/zgQe6u0Gy/Hm8+SaceGJ212RBqnxuhlFqmOhkQSwWY968eYnQ3eImVW9gA9CWL0eS\nEJHQnUWjaGhooKmpif3797NhwwbAW1+Trnfg7w2NHz/+WP2f/QymTk2qm1E/Y8cOOOusjHwebaw3\nZIwlKgrtQCmxePHiIhcc//RFkCq8Ho1QKMEBb+vpTARHRGhubmbz5s10dXVRU1MTWi8Wi9HR0UFH\nRwexmFuutWMH6zdsSPwk1m/Y4PVeROADH4i+6Zo10bMwBRIcwxhrSCbj46WMW5OR0TxAKmKxGNOn\nT2fv3r2j4tfoUZrzM1EEF4E2NTUlhpQWL17M/v37ERFOHDeOh37yk+xv8PnPwze+MZouG8aYxG2p\nMuovEROdFMRXpB84cIC+vj7+8Ic/jLp/2XMesC3F+dITGj/+OZ8KYFizZrW14L4c+Pe/yWRYzjAM\nj1yJjg2vpSAe1bR+/foCC45/2CwoOP/OsWGz0hQc/9MdOXr0WDnNdR3t7cT6++lob0+UUU0IDni5\n1sLKhmEUBhOdCDZt2kRPT08BPUg1P+MXmQX5dGrYpFqAk4qOK69MCIv/qU9uaEjkQ+vq6qKrqys0\nN5r4ItEk26g0wzBGHROdCNrb2wsQNJCp0BQnwxWWCRx7shPr62lrbU30XLq6u2lsbKSzs5OGhobE\nNU1NTQBDgwiCbbvM08GyYRgFIhcrTovpgOwzEvT19eVplX1LmoXrhc8CEDzSOBx5nAkZZXBIteI+\nuHo/k5X6lkLGMIaHe29aRoJcEQ8a2Lt3L08++WQO75Tqe/9XgDtzeO/M+CHw0WFcdxXwY9/nWbNm\n8fTTTyfKZ0+ZErmgMxipFsZw1qvYGhfDKC5MdBwLFy5k/fr1OWo9ldAUZrjsc8A3h3HdI8BfZlCv\nsrKSSb59XyZNmpRYWd/X15fYCwdg8uTJw9pMzVbqG0bpYXM6jvg38tGj8PMzs4kex0onOBJxpBKc\nE044IZFx4NFHH01a0FlTU5PodWzevJm2tjYaGhpobW3lueeeG9YGaemCCAzDKD7Kep1OLBZj4cKF\nbN26Nemb9/C4FliTypMRth/ORGDfMK8diUfxXsy+fcfuHlwHU5Q7bxqGkRG2OHSYpBKd973vfSOc\nv0n1s7sWb2+60WG4/0u5kLp58+bx5JNPEovFWLBgAc899xwzZsxg1apVJiyGMUYw0RkmqURn/Pjx\nKVPeh5O7+ZliEJbjjz+eP/3pTwCMGzeOiRMnUlFRwbRp03jllVeYPn06q1evNnExjDFOrkSnrAMJ\njj/++AxFZ/SEZrjCMh7IVh7TUVFRwYknnkh3dzezZ88GbEjMMIzcUpY9HX949HPPPZc0L+FxOvBS\nRIsvAOenvOdwf6LnAS8O89ogFRUVfOUrX+FrX/saqso555zDrl27OPPMM6mpqWHChAkmKoZhRGLD\nayGIyBXA3XhRePeq6oqQOkNEJ5gE8qmnnmLPng8C34+408kEN0H7JdA0DJ8XA/82jOuCjBs3DvAm\n9Ht6ejj11FOth2IYxqhhCT8DiEgF8C/A5cD/BXxURM7NrpUZdHd3sWfPAEMFR5iBoAjK7iEhx6kE\n57+JDjnOVnBuu+220FW9hw8f5vDhwwwMDDB79uyiDB/u7e0ttAsjwvwvLOb/2KRkRQeYA+xU1ZdV\n9RBevPIHM7mws7OT+votwBaqeZvXmeLExX/AlhRtPEm0sLT46gWTTNbU1NDX10d/fz/t7e20t7fT\n398fmTJi/Pjxmf00ipBS/6Mz/wuL+T82KeVAgkbgVd/nfjwhSn9hYyOxS7/EhAfvT1nvVuCuLBwa\nN24cF1xwAbW1tRnNmVh6FsMwyo1SFp2M+da3vjXENuHWW+DKVn7/3vfyiS99CUieC4nFYvQtWED9\nr34FwLnnnsv48eOpqamxORPDMIxhUrKBBCIyF1imqle4z7fhZUVdEahXmg9oGIZRYCx6zYeIjAO2\nAxcDv8cLKPuoqr5QUMcMwzCMSEp2eE1Vj4jI/wP0cCxk2gTHMAyjiCnZno5hGIZRepRyyHRKROQK\nEXlRRHaIyK2F9ieOiEwVkfUi8ryIbBWRzzh7nYj0iMh2EVknIrW+a5aIyE4ReUFELvPZZ4rIFveM\nd+fxGSpE5BkReawEfa8VkR85f54XkYtKzP8lzu8tIvIDEakuZv9F5F4R2SUiW3y2UfPXPf8ad81T\nInJaHvy/y/n3rIg8IiKTfOeK3n/fuc+LyFERqc+r/7nYjrTQB56Y/hp4F1AFPAucW2i/nG+nABe4\n8gl481LnAiuALzr7rcDXXfl8YDPeUOjp7rniPdRfAE2u3A1cnqdn+Bvg34HH3OdS8v0+4AZXrgRq\nS8V/9/v8W6DafX4QWFTM/gPvBy4Atvhso+Yv8GngHle+FliTB/8vASpc+evA10rJf2efCqzFy/dV\n72zn5cP/nP+RF+IA5gKP+z7fBtxaaL8ifP0v90v8IjDF2U4BXgzzHXgcuMjV2eazzwdW5sHfqcBP\n8NbAxkWnVHyfBPwmxF4q/tc5X+vci+GxUvjdwRNL/0t71PzFe3Fe5MrjgN259j9w7mpgdan5D/wI\nmEGy6OTF/7E6vBa2cLToFtaIyOl430I24v0R7gJQ1dfxEr7B0GeJOVsj3nPFydczfgv4Asl5TUvF\n92nAmyLyfTc82Ckix1Ei/qvqHrxNX19xvgyq6hOUiP8+Th5FfxPXqOoRYK9/uCgPfBzvm3+SL46i\n9F9ErgJeVdWtgVN58X+sik7RIyInAA8Dt6jqHxmanLroIjxEpAPYparPknpPh6Lz3VEJzAS+o6oz\ngT/hfbsr+p89gIicgTe0+S7gHcDxIvJXlIj/KRhNf3O7B7z/RiJfAg6p6gOj2ewotjW0cZEa4HZg\naa5uka7CWBWdGOCf0JrqbEWBiFTiCc5qVX3UmXeJyBR3/hTgDWePAe/0XR5/lih7LnkfcJWI/BZ4\nAGgTkdXA6yXgO3jf0F5V1U3u8yN4IlQKP3uA2cDPVXXAfav8T6CZ0vE/zmj6mzgn3tq9Sao6kDvX\nPUTkeqAd+JjPXAr+vxtvvuZXIvKS8+UZETmZ6PfmqPo/VkWnDzhTRN4lItV4Y5CPFdgnP/8v3hjp\nt322x4DrXXkR8KjPPt9FiUwDzgR+6YYlBkVkjogIsNB3TU5Q1dtV9TRVPQPvZ7peVRcAPy52353/\nu4BXReRsZ7oYeJ4S+Nk7tgNzRWSCu+/FwLYS8D+eCzfOaPr7mGsD4CPA+lz7L96WKl8ArlLVt331\nit5/VX1OVU9R1TNUdRreF7ELVfUN58u1Ofd/tCetiuUArsD7I90J3FZof3x+vQ84ghdRtxl4xvla\nDzzhfO4BJvuuWYIXSfICcJnPPgvY6p7x23l+jr/gWCBByfgOvBfvS8mzwH/gRa+Vkv9fwBPKLcD9\neNGZRes/8EPgNeBtvLmoG/ACIUbFX7xNdR9y9o3A6XnwfyfwsvvbfQYXvVUq/gfO/xYXSJAv/21x\nqGEYhpE3xurwmmEYhlGEmOgYhmEYecNExzAMw8gbJjqGYRhG3jDRMQzDMPKGiY5hGIaRN0x0DMMw\njLxhomMYhmHkjf8fPsUet5isnloAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d702b90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(simple_feature_matrix, output,'k.',\n",
    "        simple_feature_matrix, predict_output(simple_feature_matrix, simple_weights_0_penalty),'b-',\n",
    "        simple_feature_matrix, predict_output(simple_feature_matrix, simple_weights_high_penalty),'r-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  1.00000000e+00,   1.18000000e+03],\n",
       "       [  1.00000000e+00,   2.57000000e+03],\n",
       "       [  1.00000000e+00,   7.70000000e+02],\n",
       "       ..., \n",
       "       [  1.00000000e+00,   1.53000000e+03],\n",
       "       [  1.00000000e+00,   1.60000000e+03],\n",
       "       [  1.00000000e+00,   1.02000000e+03]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "simple_feature_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_rss(feature_matrix, weights, output):\n",
    "    predicted_test_output = predict_output(feature_matrix, weights)\n",
    "    rss = np.sum((predicted_test_output - output)**2)\n",
    "    return rss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.78427328252e+15\n",
      "2.75723634598e+14\n",
      "6.94642100914e+14\n"
     ]
    }
   ],
   "source": [
    "rss_1 = get_rss(simple_test_feature_matrix, [0., 0.], test_output)\n",
    "rss_2 = get_rss(simple_test_feature_matrix, simple_weights_0_penalty, test_output)\n",
    "rss_3 = get_rss(simple_test_feature_matrix, simple_weights_high_penalty, test_output)\n",
    "print rss_1\n",
    "print rss_2\n",
    "print rss_3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running a multiple regression with L2 penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "model_features = ['sqft_living', 'sqft_living15']\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": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  -0.35743482  243.0541689    22.41481594]\n"
     ]
    }
   ],
   "source": [
    "multiple_weights_0_penalty = ridge_regression_gradient_descent(feature_matrix, \n",
    "                                                             output, \n",
    "                                                             initial_weights=[0., 0., 0.], \n",
    "                                                             step_size=1e-12, \n",
    "                                                             l2_penalty=0, \n",
    "                                                             max_iterations=1000)\n",
    "print multiple_weights_0_penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  6.7429658   91.48927361  78.43658768]\n"
     ]
    }
   ],
   "source": [
    "multiple_weights_high_penalty = ridge_regression_gradient_descent(feature_matrix, \n",
    "                                                             output, \n",
    "                                                             initial_weights=[0., 0., 0.], \n",
    "                                                             step_size=1e-12, \n",
    "                                                             l2_penalty=1e11, \n",
    "                                                             max_iterations=1000)\n",
    "print multiple_weights_high_penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.78427328252e+15\n",
      "2.74067618287e+14\n",
      "5.0040480058e+14\n"
     ]
    }
   ],
   "source": [
    "rss_1 = get_rss(test_feature_matrix, [0., 0., 0.], test_output)\n",
    "rss_2 = get_rss(test_feature_matrix, multiple_weights_0_penalty, test_output)\n",
    "rss_3 = get_rss(test_feature_matrix, multiple_weights_high_penalty, test_output)\n",
    "print rss_1\n",
    "print rss_2\n",
    "print rss_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "77465.476464743959"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predict_output(test_feature_matrix, multiple_weights_0_penalty)[0] - test_output[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-39546.469695141423"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "predict_output(test_feature_matrix, multiple_weights_high_penalty)[0] - test_output[0]"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}