Notebook12_ThetaLMS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-ae822b978b700c32",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# Part 4: \"Online\" linear regression\n",
    "\n",
    "When you are trying to fit a model to data and you get to see all of the data at once, we refer to the problem as an _offline_ or _batch_ problem, and you would try to use certain algorithms to compute the fit that can take advantage of the fact that you a lot of available data.\n",
    "\n",
    "But what if you only get to see one or a few data points at a time? In that case, you might want to get an initial model from whatever data you've got, and gradually improve the model as you see new data points. In this case, we refer to the problem as being an _online_ problem.\n",
    "\n",
    "The goal of this notebook is to introduce you to online algorithms. You'll start by reviewing the offline linear regression problem, and then look at its online variant. The neat thing about the online method is that you can derive it using all the tools you already have at your disposal, namely, multivariate calculus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-76e9aaecf6da8142",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a74697ce3692dcbe",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Review: Offline or batch linear regression\n",
    "\n",
    "Let's start with a quick review of the linear regression problem: given a response vector, $y$, and a data matrix $X$---whose rows are observations and columns are variables---the problem is to find the best linear model, $y \\approx X \\theta^*$, where $\\theta^*$ is the vector of best-fit model parameters that we wish to compute. Computing it using a conventional batch linear least squares method has an asymptotic running time of $\\mathcal{O}(mn^2)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-b753a6dc289159c5",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "To start, here is some code to help generate synthetic problems of a certain size, namely, $m \\times (n+1)$, where $m$ is the number of observations and $n$ the number of predictors. The $+1$ comes from our usual dummy coefficient for a non-zero intercept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a013b10e1d1d2d36",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "def generate_model (n):\n",
    "    \"\"\"Returns a set of (random) n+1 linear model coefficients.\"\"\"\n",
    "    return np.random.rand (n+1, 1)\n",
    "\n",
    "def generate_data (m, theta, sigma=1.0/(2**0.5)):\n",
    "    \"\"\"\n",
    "    Generates 'm' noisy observations for a linear model whose\n",
    "    predictor (non-intercept) coefficients are given in 'theta'.\n",
    "    Decrease 'sigma' to decrease the amount of noise.\n",
    "    \"\"\"\n",
    "    assert (type (theta) is np.ndarray) and (theta.ndim == 2) and (theta.shape[1] == 1)\n",
    "    n = len (theta)\n",
    "    X = np.random.rand (m, n)\n",
    "    X[:, 0] = 1.0\n",
    "    y = X.dot (theta) + sigma*np.random.randn (m, 1)\n",
    "    return (X, y)\n",
    "\n",
    "def estimate_coeffs (X, y):\n",
    "    \"\"\"\n",
    "    Solves X*theta = y by a linear least squares method.\n",
    "    \"\"\"\n",
    "    result = np.linalg.lstsq (X, y)\n",
    "    theta = result[0]\n",
    "    return theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-eec5187ec010c42a",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "def rel_diff(x, y, ord=2):\n",
    "    \"\"\"\n",
    "    Computes ||x-y|| / ||y||. Uses 2-norm by default;\n",
    "    override by setting 'ord'.\n",
    "    \"\"\"\n",
    "    return np.linalg.norm (x - y, ord=ord) / np.linalg.norm (y, ord=ord)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-0f945a8420966877",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## An online algorithm\n",
    "\n",
    "The empirical scaling of linear least squares appears to be pretty good, being roughly linear in $m$ or at worst quadratic in $n$. But there is still a downside in time and storage: each time there is a change in the data, you appear to need to form the data matrix all over again and recompute the solution from scratch, possibly touching the entire data set again!\n",
    "\n",
    "This begs the question, is there a way to incrementally update the model coefficients whenever a new data point, or perhaps a small batch of new data points, arrives? Such a procedure would be considered _incremental_ or _online_, rather than batched or offline."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-495a02d1ac83a302",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Setup: Key assumptions and main goal.** In the discussion that follows, assume that you only get to see the observations _one-at-a-time_. Let $(y_k, \\hat{x}_k^T)$ denote the current observation. (Relative to our previous notation, this tuple is just element $k$ of $y$ and row $k$ of $X$.\n",
    "\n",
    "> We will use $\\hat{x}_k^T$ to denote a row $k$ of $X$ since we previously used $x_j$ to denote column $j$ of $X$. That is,\n",
    ">\n",
    "> $$\n",
    "    X = \\left(\\begin{array}{ccc}\n",
    "          x_0 & \\cdots & x_{n}\n",
    "        \\end{array}\\right)\n",
    "      = \\left(\\begin{array}{c}\n",
    "          \\hat{x}_0^T \\\\\n",
    "            \\vdots \\\\\n",
    "          \\hat{x}_{m-1}^T\n",
    "        \\end{array}\\right),\n",
    "  $$\n",
    ">\n",
    "> where the first form is our previous \"columns-view\" representation and the second form is our \"rows-view.\"\n",
    "\n",
    "Additionally, assume that, at the time the $k$-th observation arrives, you start with a current estimate of the parameters, $\\tilde{\\theta}(k)$, which is a vector. If for whatever reason you need to refer to element $i$ of that vector, use $\\tilde{\\theta}_i(k)$. You will then compute a new estimate, $\\tilde{\\theta}(k+1)$ using $\\tilde{\\theta}(k)$ and $(y_k, \\hat{x}_k^T)$. For the discussion below, further assume that you throw out $\\tilde{\\theta}(k)$ once you have $\\tilde{\\theta}(k+1)$.\n",
    "\n",
    "As for your goal, recall that in the batch setting you start with _all_ the observations, $(y, X)$. From this starting point, you may estimate the linear regression model's parameters, $\\theta$, by solving $X \\theta = y$. In the online setting, you compute estimates one at a time. After seeing all $m$ observations in $X$, your goal is to compute an $\\tilde{\\theta}_{m-1} \\approx \\theta$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-0f00fed6467eed7b",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**An intuitive (but flawed) idea.** Indeed, there is a technique from the signal processing literature that we can apply to the linear regression problem, known as the _least mean square (LMS) algorithm_. Before describing it, let's start with an initial idea.\n",
    "\n",
    "Suppose that you have a current estimate of the parameters, $\\theta(k)$, when you get a new sample, $(y_k, \\hat{x}_k^T)$. The error in your prediction will be,\n",
    "\n",
    "$$y_k - \\hat{x}_k^T \\tilde{\\theta}(k).$$\n",
    "\n",
    "Ideally, this error would be zero. So, let's ask if there exists a _correction_, $\\Delta_k$, such that\n",
    "\n",
    "$$\n",
    "\\begin{array}{rrcl}\n",
    "     & y_k - \\hat{x}_k^T \\left( \\tilde{\\theta}(k) + \\Delta_k \\right) & = & 0 \\\\\n",
    "\\iff &                           y_k - \\hat{x}_k^T \\tilde{\\theta}(k) & = & \\hat{x}_k^T \\Delta_k\n",
    "\\end{array}\n",
    "$$\n",
    "\n",
    "Then, you could compute a new estimate of the parameter by $\\tilde{\\theta}(k+1) = \\tilde{\\theta}(k) + \\Delta_k$.\n",
    "\n",
    "This idea has a major flaw, which we will discuss below. But before we do, please try the following exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-0e0a2b5c75a065f5",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Mental exercise (no points).** Verify that the following choice of $\\Delta_k$ would make the preceding equation true.\n",
    "\n",
    "$$\n",
    "\\begin{array}{rcl}\n",
    "  \\Delta_k & = & \\dfrac{\\hat{x}_k}{\\|\\hat{x}_k\\|_2^2} \\left( y_k - \\hat{x}_k^T \\tilde{\\theta}(k) \\right).\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f2e525ee3a531640",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Refining (or rather, \"hacking\") the basic idea: The least mean square (LMS) procedure.** The basic idea sketched above has at least one major flaw: the choice of $\\Delta_k$ might allow you to correctly predicts $y_k$ from $x_k$ and the new estimate $\\tilde{\\theta}(k+1) = \\tilde{\\theta}(k) + \\Delta_k$, but there is no guarantee that this new estimate $\\tilde{\\theta}(k+1)$ preserves the quality of predictions made at all previous iterations!\n",
    "\n",
    "There are a number of ways to deal with this problem, which includes carrying out an update with respect to some (or all) previous data. However, there is also a simpler \"hack\" that, though it might require some parameter tuning, can be made to work in practice."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d0e7ca1c5c2e2c67",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "That hack is as follows. Rather than using $\\Delta_k$ as computed above, let's compute a different update that has a \"fudge\" factor, $\\phi$:\n",
    "\n",
    "$$\n",
    "\\begin{array}{rrcl}\n",
    "  &\n",
    "  \\tilde{\\theta}(k+1) & = & \\tilde{\\theta}(k) + \\Delta_k\n",
    "  \\\\\n",
    "  \\mbox{where}\n",
    "  &\n",
    "  \\Delta_k & = & \\phi \\cdot \\hat{x}_k \\left( y_k - \\hat{x}_k^T \\tilde{\\theta}(k) \\right).\n",
    "\\end{array}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-f36e2f688273c7c2",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "A big question is how to choose $\\phi$. There is some analysis out there that can help. We will just state the results of this analysis without proof.\n",
    "\n",
    "Let $\\lambda_{\\mathrm{max}}(X^T X)$ be the largest eigenvalue of $X^T X$. The result is that as the number of samples $m \\rightarrow \\infty$, any choice of $\\phi$ that satisfies the following condition will _eventually_ converge to the best least-squares estimator of $\\tilde{\\theta}$, that is, the estimate of $\\tilde{\\theta}$ you would have gotten by solving the linear least squares problem with all of the data.\n",
    "\n",
    "$$\n",
    "  0 < \\phi < \\frac{2}{\\lambda_{\\mathrm{max}}(X^T X)}.\n",
    "$$\n",
    "\n",
    "This condition is not very satisfying, because you cannot really know $\\lambda_{\\mathrm{max}}(X^T X)$ until you've seen all the data, whereas we would like to apply this procedure _online_ as the data arrive. Nevertheless, in practice you can imagine hybrid schemes that, given a batch of data points, use the QR fitting procedure to get a starting estimate for $\\tilde{\\theta}$ as well as to estimate a value of $\\phi$ to use for all future updates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-fd2a4be161147caa",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Summary of the LMS algorithm.** To summarize, the algorithm is as follows:\n",
    "* Choose any initial guess, $\\tilde{\\theta}(0)$, such as $\\tilde{\\theta}(0) \\leftarrow 0$.\n",
    "* For each observation $(y_k, \\hat{x}_k^T)$, do the update:\n",
    "\n",
    "  * $\\tilde{\\theta}(k+1) \\leftarrow \\tilde{\\theta}_k + \\Delta_k$,\n",
    "  \n",
    "  where $\\Delta_k = \\phi \\cdot \\hat{x}_k \\left( y_k - \\hat{x}_k^T \\tilde{\\theta}(k) \\right)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-1c8224a488dc995c",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Trying out the LMS idea\n",
    "\n",
    "Now _you_ should implement the LMS algorithm and see how it behaves.\n",
    "\n",
    "To start, let's generate an initial 1-D problem (2 regression coefficients, a slope, and an intercept), and solve it using the batch procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3e4219c8414d5443",
     "locked": true,
     "schema_version": 1,
     "solution": false
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Condition number of the data matrix: 4.38999326346\n",
      "Relative error: 0.00144897427049\n"
     ]
    }
   ],
   "source": [
    "m = 100000\n",
    "n = 1\n",
    "theta_true = generate_model(n)\n",
    "(X, y) = generate_data(m, theta_true, sigma=0.1)\n",
    "\n",
    "print(\"Condition number of the data matrix:\", np.linalg.cond(X))\n",
    "\n",
    "theta = estimate_coeffs(X, y)\n",
    "e_rel = rel_diff(theta, theta_true)\n",
    "\n",
    "print(\"Relative error:\", e_rel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-d7d3b934f53af706",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Recall that we need a value for $\\phi$, for which we have an upper-bound of $\\lambda_{\\mathrm{max}}(X^T X)$. Let's cheat by computing it explicitly, even though in practice we would need to do something different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-4fc601f5b57a006f",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "126813.760247\n"
     ]
    }
   ],
   "source": [
    "LAMBDA_MAX = max(np.linalg.eigvals(X.T.dot(X)))\n",
    "print(LAMBDA_MAX)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-2e86f046faca4191",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "**Exercise 1** (5 points). Implement the online LMS algorithm in the code cell below where indicated. It should produce a final parameter estimate, `theta_lms`, as a column vector.\n",
    "\n",
    "In addition, the skeleton code below uses `rel_diff()` to record the relative difference between the estimate and the true vector, storing the $k$-th relative difference in `rel_diffs[k]`. Doing so will allow you to see the convergence behavior of the method.\n",
    "\n",
    "Lastly, to help you out, we've defined a constant in terms of $\\lambda_{\\mathrm{max}}(X^T X)$ that you can use for $\\phi$.\n",
    "\n",
    "> In practice, you would only maintain the current estimate, or maybe just a few recent estimates, rather than all of them. Since we want to inspect these vectors later, go ahead and store them all."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "lms",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "PHI = 1.99 / LAMBDA_MAX # Fudge factor\n",
    "rel_diffs = np.zeros((m+1, 1))\n",
    "\n",
    "theta_k = np.zeros((n+1))\n",
    "for k in range(m):\n",
    "    rel_diffs[k] = rel_diff(theta_k, theta_true)\n",
    "\n",
    "    # Implement the online LMS algorithm.\n",
    "    # Take (y[k], X[k, :]) to be the k-th observation.\n",
    "    \n",
    "    # YOUR CODE HERE\n",
    "    delta_k = PHI * X[k, :] * (y[k] - X[k, :].T.dot(theta_k))\n",
    "    theta_k = theta_k + delta_k\n",
    "    \n",
    "theta_lms = theta_k\n",
    "rel_diffs[m] = rel_diff(theta_lms, theta_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-486af1b126a34491",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Let's compare the true coefficients against the estimates, both from the batch algorithm and the online algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "lms_test",
     "locked": true,
     "points": 5,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.6573746   0.65255162]]\n",
      "[[ 0.65655272  0.65361267]]\n",
      "[ 0.66331619  0.38456212]\n",
      "\n",
      "('Passed' -- this cell appears to run without error, but we aren't checking the solution.)\n"
     ]
    }
   ],
   "source": [
    "print (theta_true.T)\n",
    "print (theta.T)\n",
    "print (theta_lms.T)\n",
    "\n",
    "print(\"\\n('Passed' -- this cell appears to run without error, but we aren't checking the solution.)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-cbeeadccc0a71818",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Let's also compute the relative differences between each estimate `Theta[:, k]` and the true coefficients `theta_true`, measured in the two-norm, to see if the estimate is converging to the truth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-1e3a7e61238da6b6",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe586dd5e10>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd8leX9//HXJycLEgjZYBJIkLA3BwQFxUFFHOAerThL\nratqraP1qz+rHXa4WqmiAlXrHnWh1I2CIkH23hBWwgqEEQi5fn+cA41USICT3DnnvJ+PRx6cc99X\nzv25ucObO9d93ddtzjlERCSyxHhdgIiIhJ7CXUQkAincRUQikMJdRCQCKdxFRCKQwl1EJAIp3EVE\nIpDCXUQkAincRUQiUKxXG87IyHD5+flebV5EJCxNnTp1g3Mus6Z2noV7fn4+RUVFXm1eRCQsmdmK\n2rRTt4yISARSuIuIRCCFu4hIBFK4i4hEIIW7iEgEUriLiEQghbuISASqMdzNbLSZlZjZ7Bra9Taz\nSjO7IHTl/a9Vm3Zw/7tz2LO3qi43IyIS1mpz5j4WGHyoBmbmAx4C/hOCmg5pwbptjJm4nBe+qdU4\nfhGRqFRjuDvnJgCbamh2E/AGUBKKog7l1A5Z9G+TwSMfLWTz9t11vTkRkbB01H3uZpYDnAv84+jL\nqdX2uPfsjpRXVPLwRwvrY5MiImEnFBdUHwXudM7V2AluZiPMrMjMikpLS494g22zm3B531a8MHkF\nM4u3HPHniIhEqlCEux942cyWAxcAI81s2A81dM6Ncs75nXP+zMwaJzU7pF+e3o6M5AR+89Zs9la5\no/osEZFIc9Th7pwrcM7lO+fygdeB651z/z7qymrQNDGOe87swKzVZfxrsi6uiohUV5uhkC8BXwPt\nzKzYzK4xs+vM7Lq6L+/Qzul2DP3bZPDnDxdQsm2X1+WIiDQY5pw3XRp+v9+FYj73ZRu2c/ojExjU\nKZsnLusZgspERBouM5vqnPPX1C7s71AtyEjiplPa8P7MtXw0d73X5YiINAhhH+4A1w08lnbZTbjv\n7dns2F3pdTkiIp6LiHCP88Xw4LmdWVO2i8c+WeR1OSIinouIcAfonZ/GRf5cnv1yGQvXb/O6HBER\nT0VMuAPcdUYHkhJiueet2VRp7LuIRLGICve0pHh+PaQ93y7fxCtFq7wuR0TEMxEV7gAX+fPo2zqN\n34+bR8lWjX0XkegUceFuZvzhvK5UVFbxwPvzvC5HRMQTERfuEBj7fv3AY3l3xhq+XHTkE5SJiISr\niAx3gOtOOpb89Mb8379ns2vPXq/LERGpVxEb7olxPn5/bheWb9zBox9r7LuIRJeIDXeA49tkcJE/\nl6e/XMrs1WVelyMiUm8iOtwBfnNmR9KT4rnj9Zl6qLaIRI2ID/eURnE8MKwzc9du5akvlnhdjohI\nvYj4cAc4vVNzzuzagsc/WcwiTU0gIlEgKsId4P5zOpGU4OOON2bqsXwiEvGiJtwzkhO47+xOTFu5\nhbGTlntdjohInYqacAcY2v0YTmmfxZ/Hz2dJabnX5YiI1JmoCvfA1ARdSIzzcdsr06nU6BkRiVBR\nFe4A2U0TeXBYZ2YUl/HUhKVelyMiUieiLtwBzup6DGd1bcEjHy3UzU0iEpGiMtwBHhzWmfTkeG59\nZbrmnhGRiBO14d6scTx/uqAbi0rK+fP4BV6XIyISUlEb7gAntc1keL9WPPvVMiYt3uB1OSIiIRPV\n4Q5w9xkdaJ2RxO2vzaBs5x6vyxERCYmoD/dG8T4evrg767dVcP87c7wuR0QkJKI+3AG65zXjhpPb\n8Oa01YybtdbrckREjprCPeimU9rQNTeFX781Sw/WFpGwp3APivPF8PBF3dm5ey93vDET5zS5mIiE\nL4V7NW2ykrn7jPZ8vqCUF79d6XU5IiJHTOF+gOH98hlQmMGD781jqSYXE5EwVWO4m9loMysxs9kH\nWf9jM5tpZrPMbJKZdQt9mfUnJsb48wXdiI+N4ZZXprO7UpOLiUj4qc2Z+1hg8CHWLwNOcs51AR4A\nRoWgLk81T0nkofO7MLO4jD99ON/rckREDluN4e6cmwBsOsT6Sc65zcG33wC5IarNU4M7t2B4v1Y8\n89UyPpm33utyREQOS6j73K8BPgjxZ3rm10M60LFFU3752gzWbNnpdTkiIrUWsnA3s5MJhPudh2gz\nwsyKzKyotLQ0VJuuM4lxPp74cU/2VFZx80vT2KOHe4hImAhJuJtZV+AZYKhzbuPB2jnnRjnn/M45\nf2ZmZig2XecKMpL4/XldKFqxmUc+Wuh1OSIitXLU4W5mLYE3gcudcxGZfkO753BpnzxGfr6ECQsb\n/m8cIiK1GQr5EvA10M7Mis3sGjO7zsyuCza5F0gHRprZdDMrqsN6PXPf2Z1om53Mba/OYGN5hdfl\niIgcknl1m73f73dFReH1/8C8tVsZ+veJ9C/M4JnhfmJizOuSRCTKmNlU55y/pna6Q/UwdGjRlHvO\n6sCn80t49qtlXpcjInJQCvfDdHnfVgzu1JyHPpzPtJWba/4GEREPKNwPk5nx0AVdaZ6SyI0vTqNs\nh57eJCINj8L9CKQ0iuNvl/Zg/dZd3KnpgUWkAVK4H6EeLVO564z2fDhnHc9/s8LrckREvkfhfhSu\n6V/Aqe2zePC9ecwqLvO6HBGR/RTuR8HM+MuF3chIjue6F6Zq/LuINBgK96OUmhTPk5f3orS8ghtf\n1PwzItIwKNxDoGtuM/5wbhe+XrqR370/z+tyRESI9bqASHF+r1zmrNnK6InL6NUqlbO7HeN1SSIS\nxXTmHkJ3D2lPr1ap3PXGTBat3+Z1OSISxRTuIRTni+GJy3rSKD6WEc9PpWynbnASEW8o3EOseUoi\nI3/ck1WbdnDrK9OpqtINTiJS/xTudaBPQRr3nt2RT+eX8OjHETnFvYg0cAr3OnJ531Zc2CuXxz9d\nzPg567wuR0SijMK9jpgZDwzrTLfcFH756gyWlpZ7XZKIRBGFex1KjPMx8ie9iPMZ1z5XpAusIlJv\nFO51LKdZI578SS9WbdrBjS9+R6XuYBWReqBwrwfHtU7nwWGd+XLRBn773lyvyxGRKKA7VOvJxb1b\nsqR0O6MmLKVNVjLD++V7XZKIRDCFez26c3B7lpaWc/+7c8lPT+LEtplelyQiEUrdMvXIF2M8ekkP\nCrOSuf5f3zF/3VavSxKRCKVwr2fJCbGMuao3SQk+rh4zhZKtu7wuSUQikMLdAy1SGvHsFb3ZsnMP\n1/yziJ2793pdkohEGIW7RzrnpPD4JT2YvaaMm17SEEkRCS2Fu4dO65jN/ed04uN5Jdz/7lyc0yRj\nIhIaGi3jseH98lm9eSdPTVhKTmojrjvpWK9LEpEIoHBvAO4c3J41Zbv44wfzaZGSyNDuOV6XJCJh\nTuHeAMTEGH+5sCslW3dx+2szyGySwPHHZnhdloiEMfW5NxAJsT5GDfdTkJHEiOemMmPVFq9LEpEw\nVmO4m9loMysxs9kHWW9m9riZLTazmWbWM/RlRoeURnE8d/VxpCbFcfXYKazcuMPrkkQkTNXmzH0s\nMPgQ688ACoNfI4B/HH1Z0at5SiJjr+pDZZVj+OjJbCiv8LokEQlDNYa7c24CsOkQTYYCz7mAb4Bm\nZtYiVAVGo2Mzkxl9pZ+1Zbu4ZuwUtldUel2SiISZUPS55wCrqr0vDi6To9CrVRp/v6wns1aXcf2/\nvmOPbnISkcNQrxdUzWyEmRWZWVFpaWl9bjosDeqYze/O7cIXC0u5/bUZ7K3STU4iUjuhGAq5Gsir\n9j43uOx/OOdGAaMA/H6/kqoWLu3Tkk3bd/Pn8QtITojlwWGdMTOvyxKRBi4U4f4OcKOZvQwcB5Q5\n59aG4HMl6IaT27BtVyVPfrGElEZx3DG4vdcliUgDV2O4m9lLwEAgw8yKgfuAOADn3JPAOGAIsBjY\nAVxVV8VGszsHt6Ns5x5Gfr6EJolx/HygpikQkYOrMdydc5fWsN4BN4SsIvlBZsaDwzpTXlHJQx/O\nxwzNQyMiB6XpB8KIL8Z45KJuAPzxg/nExhjXDmjtcVUi0hAp3MNMrC+GRy7qRlWV48H355EQG8Pl\neti2iBxA4R6GYn0xPHpJdyoq9/J/b88hIc7HRf68mr9RRKKGJg4LU3G+GP5+WU8GFGZw5xszeWNq\nsdcliUgDonAPY4lxPp4e7ueEYzO4/fUZ/HvaD95eICJRSOEe5vYFfN+CdG57dTrvzFjjdUki0gAo\n3CNAo3gfz17pp3d+Gre+Mp33Z+oeMpFop3CPEI3jYxl9ZW96tUzl5pen8cEsBbxINFO4R5CkhFhG\nX9Wb7nnNuOklBbxINFO4R5jkhFjGXtWbbnnNuPGlabw3U33wItFI4R6BmiTG8c+r+9CzZTN+8fJ0\n3p6uUTQi0UbhHqECZ/B96J2fyq2vTOetaRoHLxJNFO4RLCkhljFX9qFv63Rue3UGL3+70uuSRKSe\nKNwjXKN4H89e0ZsTCzO5681ZjJm4zOuSRKQeKNyjQKN4H6OG9+JHHbO5/925/O2TRQRmahaRSKVw\njxIJsT5G/rgn5/XI4a8fLeT+d+dSpWeyikQszQoZRWJ9Mfzlwm6kJsXz7FfL2F5RyR/O60KsT//H\ni0QahXuUiYkx7jmzA8kJsTz2ySK27trDY5f0IDHO53VpIhJCOmWLQmbGrYPacu9ZHRk/Zz1Xj53C\ntl17vC5LREJI4R7Fru5fwF8v7MbkZZu4+KlvKNm2y+uSRCREFO5R7vxeuTxzhZ9lG7Zz3shJLCkt\n97okEQkBhbtwcrssXh7Rl52793LBPybx3crNXpckIkdJ4S4AdMtrxpvXH0/TRnFc9vQ3fDR3vdcl\nichRULjLfq3Sk3jj58fTLrsJP3u+iH9NXuF1SSJyhBTu8j0ZyQm8NKIvJ7XN5Ddvzeav/1mgu1lF\nwpDCXf5H4/hYnh7u55Leefzt08Xc/tpMdldWeV2WiBwG3cQkPyjWF8MfzutC85REHv14EcWbd/DU\n5b1o1jje69JEpBZ05i4HZWbcclpbHr24O9NWbuG8kZNYuXGH12WJSC0o3KVGw3rk8MK1x7Fpx26G\njZzI1BWbvC5JRGqgcJda6VOQxps/P56mibFc+vRkXp2yyuuSROQQFO5Sa60zk3nr+hPok5/GHW/M\n5P5351C5VxdaRRqiWoW7mQ02swVmttjM7vqB9S3N7DMzm2ZmM81sSOhLlYYgNSmesVf15uoTChgz\ncTlXjPmWzdt3e12WiBygxnA3Mx/wBHAG0BG41Mw6HtDsHuBV51wP4BJgZKgLlYYj1hfDvWd35E8X\ndGXKss2c88RXzFu71euyRKSa2py59wEWO+eWOud2Ay8DQw9o44CmwdcpwJrQlSgN1UX+PF75WV92\nV1Zx/j8m8fb01V6XJCJBtQn3HKD61bPi4LLq/h/wEzMrBsYBN4WkOmnwerRM5d0b+9PpmKb84uXp\n/N+/Z1NRudfrskSiXqguqF4KjHXO5QJDgOfN7H8+28xGmFmRmRWVlpaGaNPitaymibz4076MOLE1\nz3+zgguf/JpVmzQeXsRLtQn31UBetfe5wWXVXQO8CuCc+xpIBDIO/CDn3CjnnN8558/MzDyyiqVB\nivPF8OshHXjq8l4sK93OmY9/ySfzNLOkiFdqE+5TgEIzKzCzeAIXTN85oM1K4FQAM+tAINx1ah6F\nTu/UnPdu7k9eWmOu+WcRD3+0kL1VmnhMpL7VGO7OuUrgRmA8MI/AqJg5ZvZbMzsn2OyXwE/NbAbw\nEnCl01SCUWvf1MHn98zl8U8WcfmzkynZqkf4idQn8yqD/X6/Kyoq8mTbUj+cc7w2tZh7355NUnws\nD1/cnZPaqjtO5GiY2VTnnL+mdrpDVeqMmXGRP493b+xPenI8V4z+lj98ME93tYrUA4W71LnC7Ca8\nfUN/LjuuJU99sZSLR31D8WaNphGpSwp3qReN4n38/twuPHZJdxas28aQx77k0/kaTSNSVxTuUq+G\nds/h/Zv7k5vamKvHFvG79+fqpieROqBwl3rXKj2JN68/np/0bcnTXy7jnL9NZO4azU0jEkoKd/FE\nYpyPB4d1YcyVvdm0YzdDn/iKkZ8v1ph4kRBRuIunTm6fxfhbTmRQx2z+9OECLnrqa1Zs3O51WSJh\nT+EunktLiueJy3ry6MXdWbh+G2c89iUvTl6J7oMTOXIKd2kQzIxhPXIYf8uJ9GjZjF+/NYurx07R\nna0iR0jhLg3KMc0a8fzVx3Hf2R2ZtGQjpz86gXGz1npdlkjYUbhLgxMTY1x1QgHv3zyAvLTGXP+v\n77jhxe/YWF7hdWkiYUPhLg1Wm6xk3vj58fxyUFv+M2cdgx6ZwLsz1qgvXqQWFO7SoMX5Yrjp1ELe\nu2kAuamNuOmlaVw9doqmLxCpgcJdwkK75k148+fHc8+ZHZi8bBM/emQCz361TOPiRQ5C4S5hI9YX\nw7UDWvOfW0/kuII0HnhvLueNnMjs1WVelybS4CjcJezkpjZm9JW9eeyS7qzesothT0zkgffmsr2i\n0uvSRBoMhbuEJTNjaPccPr7tRC705zJ64jIGPfwF789cqwuuIijcJcw1axzPH87ryuvX9aNpozhu\nePE7Ln7qG3XVSNRTuEtE6NUqjfdvHsDvz+3CktJyzv77V9zx+gxKtukOV4lOCneJGL4Y47LjWvLZ\nrwby0wGteWvaak75yxf84/MlmjNeoo7CXSJO08Q4fj2kA/+59ST6tk7noQ/nM+jhCXw4e5364yVq\nKNwlYhVkJPHMFX5euOY4EuNiuO6FqVz29GQ9GESigsJdIl7/wgzG3TyAB4Z2Yv66rZz1ty+5+81Z\nbNBcNRLBFO4SFWJ9MVzeL5/Pbz+ZK48v4LWiVZz85895esJSdldWeV2eSMgp3CWqpDSO496zO/Lh\nLSfiz0/ld+PmcfqjE/h47nr1x0tEUbhLVGqTlcyYq/ow9qrexBhc+1wRlz/7LXPWaHy8RAaFu0S1\nge2y+PCWE7nv7I7MXlPGWX/7ihte/I7FJeVelyZyVMyrX0X9fr8rKiryZNsiP6Rs5x5GTVjC2InL\n2blnL0O753DDyW1ok5XsdWki+5nZVOecv8Z2CneR79tYXsFTE5by/NcrqKjcy0X+PG4b1Jaspole\nlyaicBc5WhvLK/j7Z4t54ZsVgbtf+7TiZye1JlshLx6qbbjXqs/dzAab2QIzW2xmdx2kzUVmNtfM\n5pjZi4dbsEhDk56cwH1nd+Lj205iSJcW/PPr5Qz402fc+/Zs1mzZ6XV5IodU45m7mfmAhcAgoBiY\nAlzqnJtbrU0h8CpwinNus5llOedKDvW5OnOXcLNy4w5Gfr6Y16cWYwYX9Mrl+oFtyEtr7HVpEkVC\neebeB1jsnFvqnNsNvAwMPaDNT4EnnHObAWoKdpFw1DK9MX88vyuf/2ogF/fO442pqxn4l8+57dXp\nGl0jDU5twj0HWFXtfXFwWXVtgbZmNtHMvjGzwaEqUKShyU1tzIPDujDhjpO5ol8+42atZdAjX/Cz\n54v4buVmr8sTASA2hJ9TCAwEcoEJZtbFObeleiMzGwGMAGjZsmWINi3ijeYpidx7dkduOPlYxkxc\nzvPfrGD8nPX0KUhjxIDWnNI+i5gY87pMiVK1OXNfDeRVe58bXFZdMfCOc26Pc24ZgT76wgM/yDk3\nyjnnd875MzMzj7RmkQYlPTmB209vx6S7TuGeMzuwevNOrn2uiNMe+YLnvl6uZ7uKJ2pzQTWWQFif\nSiDUpwCXOefmVGszmMBF1ivMLAOYBnR3zm082OfqgqpEqsq9Vbw/ay2jJy5nxqotNEmM5WJ/HsP7\n5dMyXRdf5ejU9oJqjd0yzrlKM7sRGA/4gNHOuTlm9lugyDn3TnDdj8xsLrAX+NWhgl0kksX6Yhja\nPYeh3XP4buVmxkxczphJyxk9cRmndcjmyhPy6dc6HTN12Ujd0U1MIvVgXdkunvt6OS9PWcWm7bsp\nzEpmeL9WnNszl+SEUF36kmigO1RFGqBde/by7ow1PPf1CmatLiM5IZbze+Zweb98zWEjtaJwF2nA\nnHNMX7WF575ewfsz17J7bxUntEnnIn8ep3dqTmKcz+sSpYFSuIuEiQ3lFbwyZRUvTl7J6i07adY4\njvN75nJpn5Y6m5f/oXAXCTNVVY5JSzby0rcrGT9nHZVVjuMK0ri4d+BsPkl984LCXSSslW6r4LWp\nq3j521Ws3LSDpHgfQ7q04IJeufQpSNNImyimcBeJAM45pizfzGtFq/hg9jrKKyrJTW3EOd2O4bye\nueq2iUIKd5EIs2N3JR/MWsfbM9YwcfEG9lY5uuamcG6PHM7pdgzpyQlelyj1QOEuEsFKtu3inelr\n+Pf01cxevZXYGOOktpmc0/0YTuuQrf75CKZwF4kSC9Zt481pxbw9bQ3rtu4iMS6Gk9pmcnqn5pza\nIZuURnFelyghpHAXiTJVVY4pyzcxbtZaPpyzjvVbK4j3xTCgMIMhXVpwWkcFfSRQuItEsaoqx/Ti\nLYybuZZxs9aypmwXcT7jxMJMzuzaQmf0YUzhLiLAwYO+f5sMBnVszints2ieood+hwuFu4j8j31B\n/8GstXwwex3FmwMP+u6Wm8Kgjtn8qFNzCrOSNY6+AVO4i8ghOedYVFLOR3PX85+565mxKvDgtFbp\njTmtQzandcjGn59KnK82z/SR+qJwF5HDUrJ1Fx/NW89Hc9czafFGdu+toklCLP0LMzi1QzYD22WS\nobH0nlO4i8gRK6+o5KtFpXy+oJTPFpSwfmsFZtA1J4WB7bI4pX0WnXNS8OkZsfVO4S4iIeGcY86a\nrXw6v4TPF5QwbdUWnIOURnEMKMxgQGEGxx+bQV6aHiFYHxTuIlInNm3fzZeLSvly0Qa+WFhK6bYK\nAFqmNeaENhmc0Cadfq3TNR1CHVG4i0id23dRduLiDUxcvJHJSzeyraISgA4tmtK3dRrHFaTjz09V\nf32IKNxFpN5V7q1i1uoyJi3ZyMTFG/hu5WZ27akCAmf2PVs2w5+fRp+CNNpkJhOjPvvDpnAXEc9V\nVO5lVnEZU1dsZtrKLRSt2MyG8kA3TrPGcfhbpdE7PxV/fiqdc1JIiNXjBWtS23DX1HEiUmcSYn34\n89Pw56cBgW6cFRt38O3yTUxZtokpyzfx8bz1AMTHxtAtNyXQvlUqPVumkpoU72X5YU1n7iLiqdJt\nFUxdsZmi5ZsoWrGZ2avLqKwK5FLrjCS65zWje8tmdM9rRvvmTYmPje6bqtQtIyJhaefuvcwsDnTh\nTFu5hemrtuzvyomPjaHzMU3pnpdK95bN6JabQsu0xlE1XYLCXUQignOONWW7mL5yC9NXbWb6qi3M\nWl22/0JtSqM4uuam0LFFUzq0aEq75k1ok5UcsdMmqM9dRCKCmZHTrBE5zRpxZtcWAOzZW8WCdduY\ntbqMGcGwHzNxObv3BgI/3hdDYXYyHYKBX5iVTGF2Ms2bJkbNWb7O3EUkIuzZW8XS0u3MX7eVuWu3\nMm/tNuau2bq/SwegSUIsbbKTaZfdhLbZTWjXvAmF2clkJieETeirW0ZEBNhQXsGi9eUsLtnGwvXl\nLFy/jYXrt7F5x579bVIbx+3vzinISKZ1RhIFGUnkpjYitoF176hbRkQEyEhOICM5gX7Hpu9f5pxj\nQ/luFq3fxoJg2M9ft413pq9h667K/e1iY4yW6Y33h31BRjIFGUm0zkwiq0nDPttXuItI1DEzMpsk\nkNkkgePbZOxf7pxj8449LNtQztLS7Szb8N+vLxdtoKKyan/bxvG+YOAnBcI/Mxj+6UmkNPb+EYYK\ndxGRIDMjLSmetKQ0erVK+966qirH2q27WFa6PRD+wdCftbqMcbPWUlWthzs9KX5/8OdXC//89CQS\n4+rnLtxahbuZDQYeA3zAM865Px6k3fnA60Bv55w61EUkYsTE/HfUTv/CjO+t211ZxcpNO4Jn+eUs\n27CdpaXb+WJhKa9NLd7fzgxaNE3kqhMK+OmJreu03hrD3cx8wBPAIKAYmGJm7zjn5h7QrgnwC2By\nXRQqItJQxcfG0CYrmTZZyUD299aVV1SyfMP2wJl+6XaWb9xOVtO6nyGzNmfufYDFzrmlAGb2MjAU\nmHtAuweAh4BfhbRCEZEwlpwQS+ecFDrnpNTrdmszxicHWFXtfXFw2X5m1hPIc869f6gPMrMRZlZk\nZkWlpaWHXayIiNTOUQ/gNLMY4GHglzW1dc6Ncs75nXP+zMzMo920iIgcRG3CfTWQV+19bnDZPk2A\nzsDnZrYc6Au8Y2Y1DrIXEZG6UZtwnwIUmlmBmcUDlwDv7FvpnCtzzmU45/Kdc/nAN8A5Gi0jIuKd\nGsPdOVcJ3AiMB+YBrzrn5pjZb83snLouUEREDl+txrk758YB4w5Ydu9B2g48+rJERORoNKwZcURE\nJCQU7iIiEcizKX/NrBRYcYTfngFsCGE54UD7HB20z9HhaPa5lXOuxrHknoX70TCzotrMZxxJtM/R\nQfscHepjn9UtIyISgRTuIiIRKFzDfZTXBXhA+xwdtM/Roc73OSz73EVE5NDC9cxdREQOIezC3cwG\nm9kCM1tsZnd5Xc/hMLM8M/vMzOaa2Rwz+0VweZqZfWRmi4J/pgaXm5k9HtzXmcGplfd91hXB9ovM\n7Ipqy3uZ2azg9zxuDeQJvmbmM7NpZvZe8H2BmU0O1vlKcN4izCwh+H5xcH1+tc+4O7h8gZmdXm15\ng/uZMLNmZva6mc03s3lm1i/Sj7OZ3Rr8uZ5tZi+ZWWKkHWczG21mJWY2u9qyOj+uB9vGITnnwuaL\nwGP+lgCtgXhgBtDR67oOo/4WQM/g6ybAQqAj8CfgruDyu4CHgq+HAB8ARmC2zcnB5WnA0uCfqcHX\nqcF13wbbWvB7z/B6v4N13Qa8CLwXfP8qcEnw9ZPAz4OvrweeDL6+BHgl+Lpj8HgnAAXBnwNfQ/2Z\nAP4JXBt8HQ80i+TjTOAZD8uARtWO75WRdpyBE4GewOxqy+r8uB5sG4es1et/BIf5F9sPGF/t/d3A\n3V7XdRT78zaBxxcuAFoEl7UAFgRfPwVcWq39guD6S4Gnqi1/KrisBTC/2vLvtfNwP3OBT4BTgPeC\nP7gbgNh4MRSfAAAC2ElEQVQDjyuBCer6BV/HBtvZgcd6X7uG+DMBpASDzg5YHrHHmf8+1CcteNze\nA06PxOMM5PP9cK/z43qwbRzqK9y6ZWp8KlS4CP4a2oPAM2eznXNrg6vW8d+HMB5sfw+1vPgHlnvt\nUeAOoCr4Ph3Y4gIzjsL369y/b8H1ZcH2h/t34aUCoBQYE+yKesbMkojg4+ycWw38BVgJrCVw3KYS\n2cd5n/o4rgfbxkGFW7hHBDNLBt4AbnHOba2+zgX+a46YIUxmdhZQ4pyb6nUt9SiWwK/u/3DO9QC2\nE/hVer8IPM6pBJ6tXAAcAyQBgz0tygP1cVxru41wC/eangrV4JlZHIFg/5dz7s3g4vVm1iK4vgVQ\nElx+sP091PLcH1jupROAcyzwlK6XCXTNPAY0M7N9U05Xr3P/vgXXpwAbOfy/Cy8VA8XOucnB968T\nCPtIPs6nAcucc6XOuT3AmwSOfSQf533q47gebBsHFW7hfsinQjV0wSvfzwLznHMPV1v1DrDvivkV\nBPri9y0fHrzq3hcoC/5qNh74kZmlBs+YfkSgP3ItsNXM+ga3NbzaZ3nCOXe3cy7XBZ7SdQnwqXPu\nx8BnwAXBZgfu876/iwuC7V1w+SXBURYFQCGBi08N7mfCObcOWGVm7YKLTgXmEsHHmUB3TF8zaxys\nad8+R+xxrqY+juvBtnFwXl6EOcKLGUMIjDJZAvzG63oOs/b+BH6dmglMD34NIdDX+AmwCPgYSAu2\nN+CJ4L7OAvzVPutqYHHw66pqy/3A7OD3/J0DLup5vP8D+e9omdYE/tEuBl4DEoLLE4PvFwfXt672\n/b8J7tcCqo0OaYg/E0B3oCh4rP9NYFRERB9n4H5gfrCu5wmMeImo4wy8ROCawh4Cv6FdUx/H9WDb\nONSX7lAVEYlA4dYtIyIitaBwFxGJQAp3EZEIpHAXEYlACncRkQikcBcRiUAKdxGRCKRwFxGJQP8f\ngJnUNSnilKQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe586df4198>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(len(rel_diffs)), rel_diffs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-a992e44f89ef77eb",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "You should see it converging, but not especially quickly.\n",
    "\n",
    "Finally, if the dimension is `n=1`, let's go ahead and do a sanity-check regression fit plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXuYXWV97z+/TO6Jh2ASL6BJpIcUglArA9ZWzcCkfQC1\nWGI9YsZbkZBMsbY+x0uLumdQfKitfYo9JAGRIzC7WipUKUj1IRBprRfikSYmYEQgFsIzQApWIIQk\n8zt/rL1m1qxZl3fd9m1+nzzzZGbvd73rXXvPfNdvf9/f+3tFVTEMwzC6ixmtHoBhGIZRPibuhmEY\nXYiJu2EYRhdi4m4YhtGFmLgbhmF0ISbuhmEYXYiJu2EYRheSKu4icq2IPC4iP0lo0yci94rILhH5\nTrlDNAzDMLIiaYuYRORNwDPA9ar66ojnFwH/Dpylqr8QkZeo6uOVjNYwDMNwYmZaA1W9W0RWJDR5\nF3Czqv6i0d5J2JcsWaIrViR1axiGYYT50Y9+9KSqLk1rlyruDqwEZonINuBFwBWqen3aQStWrGD7\n9u0lnN4wDGP6ICJ7XdqVIe4zgVOBfmAe8D0R+b6q7okY1HpgPcCyZctKOLVhGIYRRRnZMo8A31LV\nZ1X1SeBu4DeiGqrq1araq6q9S5emfqowDMMwclKGuH8DeIOIzBSR+cDrgPtK6NcwDMPISaotIyJf\nAfqAJSLyCFADZgGo6hZVvU9E/gXYAYwB16hqbNqkYRiGUT0u2TLnO7T5K+CvShmRYRiGURhboWoY\nRlcwNNTqEbQXJu6GYXQFw8OtHkF7YeJuGIbRhZi4G4bRsQwNgYj3BRPfm0XjUFumKnp7e9VWqBqG\nURYi0CI5ayoi8iNV7U1rZ5G7YRhGF2LibhhGV1CrtXoE7YWJu2EYXYH57JMxcTcMw+hCTNwNwzC6\nEBN3wzCMLsTE3TAMowsxcTcMw+hCTNwNwzC6EBN3wzCMLsTE3TAMowsxcTcMw+hCTNwNwzC6EBN3\nwzCMLsTE3TAMowsxcTcMw+hCTNwNwzAy0gkVKE3cDcMwMtIJm3GbuBuGYXQhJu6GYRgOdNpm3Cbu\nhmEYDgwNeRtw+5tw+98niXv4uWbeCEzcDcMwKiLozQ8NNderTxV3EblWRB4XkZ+ktDtNRA6LyNvL\nG55hGEbriIu082zG3exJWJfI/cvAWUkNRKQH+Evg2yWMyTAMIzNVWB5xgpxmxUR58z7N8upTxV1V\n7wb+K6XZB4GbgMfLGJRhGN1LVlFzbV92ZJxXfMPefBS1WhuIexoicizwB8Dm4sMxDKPbySrCzbYz\n/MjbP2/RSDss9GmTsGVRxoTq3wIfU9WxtIYisl5EtovI9ieeeKKEUxuGMZ2pIj3Rj7x9XLJi4gh7\n83m8+ryUIe69wFdF5GHg7cAmEXlbVENVvVpVe1W1d+nSpSWc2jCMTiCrCLu2z5OemGWc/rmL9OfT\nDCsmiGiSMeQ3ElkB3Kqqr05p9+VGu6+l9dnb26vbt293G6VhGF2DSLIfnbd91n5d+gsK8tBQ9M2l\n2YuYRORHqtqb1s4lFfIrwPeAXxeRR0TkAhHZICIbyhioYRhGGQJZheURHFeU99/ONWZcsmXOV9WX\nq+osVX2Fqn5JVbeo6paItu9zidoNw5i+RImwL5JRIu8q2u1aBqBV2ApVwzCaSpIIR0XCrRTtOO+/\nE2rMmLgbhtESooTTf9zl2GaMyb/Z+J8e/AnbuEncVo49jNOEahXYhKphGBBfcyUpu6TsyVOXc0Sd\n06VNWr/Zx1XShKphGEaV+AJeVjpj0bHEWS5R3n8z89azYuJuGCXQbn5rp+EikmUvWCqjDIJvxaSN\nqyW14FW1JV+nnnqqGka3ANX1XatV13cribou12uNer39Y137TXvPsr6ntZrbMUV/V4Dt6qCx5rkb\nRglU6QE3w19uBXmuKxgpx/nfLt64y/nzLLaC9JWo5rkbRpvTaduuhemUcQYJZ6/4uF5Llvcsyzn8\n52q19IVNTfPpXcL7Kr7MljGy0O7WRNm2jP8RP/zl8joUsTaqpsh1qU4dc1x/SV/+uVyuPzguv33w\nMdfrGdkxonOGZys1xr/mDM/WkR0jbhceALNljG6i3a2JdrJlWlWLJSuu53dJlfQj8aAtk3Zc1rTF\nNNsniiW1k3hSdoMCEW0E4YbzbmDdyeuSBzLpXGbLGEbTaOeUuCCdaCXFVX6E6AqO/jFJx4H7exa3\n0CqqgqTPxlsHYUh4kt2NA6P7VpRL/u+A20AyYuJutC2dJERljincV5lpgmWXyC2Cf11l7XhUq3nf\np/neSXVs/MfjRDvYh9+Pf/7e2iA9wz1QEzZvb+xdlNCHzy+OcmiUA7NljI6g1RZCMymeTdEZtkye\nccQt8U9bKRplr2TtN2jLQOP7cwaZ8botjI03cruOIMuPWs7Df/qwc3uzZQzDSKRTrKQgcdF20EeP\n+/SS9JxL6d7x1+vkOvzFbGRI4PTNjKGeqOcQdlG4rP+y7Ac6YOJudASdKERZKNOCqqJEbnDDijIo\n43rDK0DjNqaOEm5XOyr4Wl7PGk/QzxuA2YdyCzra+Do0hz9gJNNkahZM3I2OoB19dp8yxlamF17F\na+ULZCs2p4i7nuBY4tokTai6lg049vPHIsPCQ2wtJug+P9yIDin6mee5aagaYQfz3A2jMFVs79YO\nXniQpDTAolvNxV1veDVq+DxRfnhc+iNMtO3rg23b4sdQ31lnwy0X8cyhZxtP5rgoGE9/7JEefnNs\nPfcMbZpyrjyY524YHUq7WFBx9dbz+NV5CPc7PBwfbQfHHBWtB1/T73wn/pyDtw0ycPMAzxx+Nn+U\nDqDQf1w/WlMOf+ow9wxtGn+9Vq/O2WdGTNwNIwdVpmm2iwUVJZT+/0XSJ11SPcOLk4ICHuWr+23i\n3oO4sdZ31lnx+VcgNWHJ5S9m8z2bkwefRGNMi+ctZmTtCHe8545JT/s3q6SbS5mYuBtGDtopX7xK\nXHPBs9zcwnnmUaVx4z4NDA9PPo/rexA53lPqvP+m9ex95lFUYP/BpwpMkPYwsnYErSlPfvTJKZOk\nUesNqsbE3TCMWMKFuoILj+KE1VW4ohYTJX1aCJ8ns3119iB80ltkRM3Lejkkz2XspEFjPAtnL2Rk\n7QgMH47MegnfrPybS/gmVQU2oWoYBSk6odjOZK2/EndMUp2X4eF4myU4iZs0Fv89mHKek+vM6P8Y\nY0c9mt8/D9JIYxx5++QURpfywTB1EVQebELVMJpEtwl7VsslLoKOisjDbYPRe5wXH/w/ySYKnkcV\nGFgDawcYW1RQ2BsivPyo5V6Ufqmy7uR1qa9TVBmD8PdVRu8WuRuGEUuZlRuDfcbVf3HZ6CKuzeBt\ng2zZvoVJmpZR1GeNCYdneAW9fFYtWcUfPrGrUIXJ4Lj916rqyN3E3TA6hKjotmpLKE9OdpqNEs5f\nd7Fdgn3D5HZrrl/D1oe2xpbVdUJh+aLlXNZ/WeqK0ahxuoh7mnXlitkyhtFlREWOWXLMyyxlkGaP\npK369Pt2bR9lbxw3tAYZFrY+uLXxoMMFRaGwaukqHv7Th8eFPetrlTa5G36+GWsZLHI3jA4hT8SY\nt22RvsKRedRzru3956FxIzu5Dm9eD3MaWS5FvXTxbJddf7xr0lNJ1+evcG3VRHppkbuIXCsij4vI\nT2KeXyciO0Rkp4j8u4j8Rp4BG0a3k7dOTNSkXVULqPKO0/V4l3ow4b6Gb67z1zPnj6cvMve5YjVe\n1Jsc5WYvLz0s7EljgYlFSK2os5MFF1vmy8BZCc8/BKxW1ZOBTwNXlzAuw+g6wmLgKm5R+eSFFu+k\n3AjiRKusLJrweeI2vZ7/mfkM4wn6s0cOFBZ0FOY9swqGlb1/9jDsXDfpGtKuL+p1aedMKSdbRkRW\nALeq6qtT2h0N/ERVj03r02wZo9Mou0BWVpukWbZMntz2IOHXKa2YV7CvwdsGvV2M/L5LSGHkhxvh\ndq9oV9jnj7uG8MRt0jWE+42iTAunVROqFwC3l9ynYbQFeSYviyzTd6nB4roFn0ubrOOMW1If9Qkl\nqR4MZw8iNUGGZKK2S8Eo3a/vokMKt28a34Ivq8DGlUMIvu4u/bbCwiktcheRM4BNwBtUdX9Mm/XA\neoBly5adunfv3hxDNozWUDRKzprfXdYEqN9PUvQYfC4plzzYPpyr7ZLSGHxO/uxYWLSvWAojTDp+\n1ZJV7L54l9MnpLjJ3bxCnPSalTuZ3cTIXUROAa4Bzo0TdgBVvVpVe1W1d+nSpWWc2jAqpUjkHddf\nqwqOufjsSf5yUl8uKY0A59XqzL10jrej0VH7Go2TzxVLwHapqcKQsvviXVPOnyXtMGmOw/85uBE3\nEPupoOzfnawUFncRWQbcDLxbVfcUH5JhtA9ZxDjtj9nVQilDEPL2MzycPM64fqNsC/91qu+ss+Ty\no5Eh4WYGOKgvFLNdGvj10n0/PWshs7wRetByS+qnlTdycEuF/ArwPeDXReQREblARDaIyIZGk08B\ni4FNInKviNgsqTEtCP+Ru7b3iatjXoYguPQTJdQwuYJh1IKiqH7D1zQh6i9m4KYB9h98Or+gw0TB\nrkZZ3ZoqW997x6Qxlx0Vuyw8apeNVaKwRUyG4Ug4CkzymMvwWMN95M24SPLci9Z5GffPA6J9Xq3O\nd+ZezP7nny4n20Wg/1X9Uza/CI4DJme2xC2UivPU0641TNZ+WpEtY+JuGDmpWtzDe32m9RknIFEr\nQMOPxf2fNr7xXYXOHoTTNk+IeUFR72Eu1629hoFT1iWmXBaZnC57wrpZWG0Zw6iAOCsjvBimjI/r\nWbdjS/J+wz/7NgZM/J+USRPFhX9XZ8GnF3grR0/f7KlJHuslIIwzdDYja0c4PHSAdSevS30dg6La\nbE+73TFxN4wMRHnO/v++sPgbOZd1viITrFHH+zeBON88bbK0vrPOws/MZ+CmAZ49UrAUABMTo1pT\nxoYPTqrKmGclbfDYohPcLrSr727ibhhtRFiQwkvzw9FpmoD5N6MowhF70mTpaUOD9AzPQIakIeoH\nCkXogrDxtI1oTce99Lj0yijCqZf+9y6Tv1knwl3H0m6Y524YOQlO0kE5E3VBspYryLKcPm1hk79y\ntBQfXWHx/MXsH7kC3RG9z2hSeYLgGOParl49eX4iTFFf3J949sfQSmxC1TBaRFUTdWkZF1Hts+5b\nCo36LsEyAEVQz3a54z13OI0fpt6A4l7P4PPB46Iomq0S/ITQIskcxyZUjWlFs6KpZkZtYS837dxR\n7aPmB4ITvn6fg7cNMmNoxkR9lxJy0lHYeNrGWNvFfyxoK0H0hHXStZdVX8dlbHn7agUWuRtdQbPS\n0VzO4xolhu2GKkUjKgqu76wzcMPFMP/pRqMCJ9DG8WOw8fSNbHrzpszjg6n7jIYJ3pQmpWKG2pQ5\noR1nGbmcp4r31WwZY1rRTuKep68qxz9FoDacBC/d3Thxwc4V5s6cy/M3XoPuWDfuTWcVtChrJc6W\nqTJnPW184TG6HFf2uMyWMbqesuqwlHmeZtpDWT4dqILue4zFfzrLE/YSbJcZMoONp23kwCcOUDvP\nmyiNW+0aN7YoS6YdbY92TXdMwiJ3oxKqthnClBUhZZ20zPp80sd8H9dyu5B+zXJKnQXvWM+zh58L\nPJh8TCwKHOnhvJ7ruGloatbLePeOE8Fp1ovLTSIoulX/vrmmUJZV4iAOs2WMltL8JdnNWUpe9Pm4\ntlmP84k7pr6zzkXfaIh6US8dQAXu2cDq5zZFphzmzczxCaZmpglj1hTRVmK2jJGZdvzo2kqqyphw\nOU+z7aEgwXMN3jbIzEt7AguNcgp7INtl1VJvz1EdHoPbo4XdH1sV1SyN/Ji4dyjtuPN6s0Qu7txp\nlLHxc9xjeYQt6yYScRHtxlsHGRYvhfGIjhUqB7Bg1gJG1o5QQ2F48uYX/jhc6OvL9ruQ9FqEfXw/\nik9bmVs1LudopVdvtkyH0s4fRaE9x+cypjJWMuZZTJN0XPCY8cj97EHovQpmjJViu8yZOYcvnful\nSXVd/HNn9Y/9Y/JmAmU9Z7B//7Vqxu9fq37HzZbpQloZGXcqLq9Zma9rWqQWFC3X8/rH1O++knl/\n3jNRhbGngLArzNaFjKwdgWHl+U88P0XY/TFCtk8lRX8f/U9C4dfSZUOOdvxE2ypM3DuIsnzNZtAs\nD9zlHGmvWZHXNdwmyzWlnXdoyJsY5ePzPR/9zos5MKeY7eJH6r1s5IXhXzFwiifoZQQKcTer1auj\n2yb1E/w/7X0Jr7atMvjppADLbJkOpR1tj6zpj2Wu9izrfFlfV9dryJoBMjQEwz9fA7+2tZSFRgtm\nL+Cqt141KTpPqt3iao3ksZnSns9qs1SdehiH2TJGJbTjooqyPhIH/yDDNkYRXKr6VfG6unwyqNUm\nZ7sMi+QX9oDgLD9qOSNrR3jmL55h3cnrIiNPf4zhMfubXMeNGfLZTFnwxT1t0+64GjrTGYvcjdKo\nIoqtYol+0X5criEc0cZFoGuuX8PWh7YS3C80F+od2yM9rD91vVNtF1+E0yYqXaPstMfjXrfVq73s\nmiJ58uHzlvmJL44s5xitj3L/BfejBycuROYIJ3zpBF667qWZzmuLmIymUOQjcZIwxJE1EyXqMdfM\nFBdcxS0sOvWddT7wjQt4/vDBUmyXVUtXseuPd2U6zH+d04Q76kaV1WaK6zuqbZFFSs0QdRf2DO5h\n39X74EhKQ4ETbzgxk8CbuBtNp4hfHScYcWTdBNn1JpQ1lTGruEMgWi9CI1JHhY2nbchUhTHptYD8\naYjhcySJbJq4x9WoqdpHz8NofZQ9H9rDkf1pSh7NnOVzeP3Dr3dub+JuNJ08+cxBGyNpNx6X/pPE\n3XWcWc/hdA1nD8JpW0CCChZ/jkQaot7/Km8DjCrz8qtcF5BUrhey3XCbyWh9lPsvuh99tsQBCfSN\n9bk3twlVo9lknYx0SUdM6zNuEi/LxF6RicDYazjb2/xiPCd9hk6kMGYVdj+FcQz44UYYUra+9462\niGDzTkBv2xY/ydwO1xXku8d+l22yjW2yjfsG7itX2IE5y+aU2p+PibtRGmUsXknKG4+r6xIlEknZ\nKVl2LPL/d1nqXt9ZZ+6lc5AhT9DH0PxiHmDV0lXokKLDCrdvmmKfVLXoqpVrFdLe96oYrY+ybca2\ncTHfJts4tO9QdScUOO6y46rp2mwZo11Iy95IwrdEslgw/nFRFkCWLBARr77L5u2bsw06ho298TsZ\nhcfXTpZFXrJOepfJaH2UBy95kIO/OEjPi3ty++Z5sGwZY9qQlr2RdmzUTSHLxJ5rFkitBsevDWS7\njHeWfdzAxOTo08sZ+aPLMtV3yZIm2ElUIe73rrmXp7c+XW6nOThm4zGs3LQy9/Glee4icq2IPC4i\nP4l5XkTkCyLygIjsEJHX5hmwMT0JWx6Q32aIah/2xNOOj7JzJnFyneHDcxi4aYDnjxzM76NDw0uX\ncR+dv32YgVPWTbl2f1xRtVZcrqtTKHsR1J7BPZPslaYKu3gi3qd9U76KCHumIaRF7iLyJuAZ4HpV\nfXXE8+cAHwTOAV4HXKGqr0s7sUXuRpiskXvWHPssaZGT2t5ch3PfCz2Nj+xF8tIbUfrieYu54uwr\nIssBJNFNdkwSea6v5ZG5wDEbikXlTqdxjNxnpjVQ1btFZEVCk3PxhF+B74vIIhF5uao+5jxaw8hB\n0HLJKga+QPrf+/351GreBOmlz30QznuqlIVG41vU1eKrL6bRjmUnyiJrpkzLxXw2nHhttgVIzSRV\n3B04FvjPwM+PNB4zcTcyUVVdl2CUnrT6dXjY89Lf/7ULOSQH4CZgQYGTjy80gv7j+nnDg/Hpi8PD\n2bJTynit2i3tcHg4+voqyS2PYxZwmCkZS+AtNjrusuPaVszDOE2oNiL3W2NsmVuBy1X13xo/bwU+\npqpTPBcRWQ+sB1i2bNmpe/fuLTR4w/BxEarYlY8Da7wiXeMNCwzE/3N6bjH8yxWw04vS0yY/W2Gz\nxNV/aZXg++NpVUTeKeLdzEVMjwKvDPz8isZjU1DVq1W1V1V7ly5dWsKpjU6mjLx4177854PC/jIe\n4/T3HTVRfTHvxCiAQg89jKwd8XLSP/fkJGGPm+wtu4piHoLna/ZmF0ND8FdyL3fJNu5kG3fJNp5q\nsrDPWzWPPu3j9Q+/vu2FPQtlRO5vBi5mYkL1C6p6elqfNqE6vagilzlrDntfn7cyUt69Bv5nSZUY\ngYVzFrLlLVumpDAmWUBZa9qURdpkctVVFaOqIzaTWcfM4nce/Z2WnLssSstzF5GvAH3AEmAUqOE5\nU6jqFhER4P8AZwHPAe+PsmTCmLh3Fi5/7EltstZ4cSHNVvCfr++sc8m3P8beXz067oGXYr38cCPc\n7i02Cpf7TauTE6aVtkzcTaisIl1FC2sVpWheebthi5iMWPJEZUkFs+LaRD2XJX0xfI4oayV8fLja\nYP91JVRghHFBX75oOfP+/TLu/8epGS9plSrLLDVctI+0G1CeCL6pE58JyFzhhGuyr/zsFEzcO4xm\nTmTliRJdamyHBcHVAsh7zrgxcHId3nwhzDnQeDD18qJpRPgLZnlb1A2csm7K0n9wv4a42jllve95\nq3LmufFCm4j5DDjmou6KzNMwce8wmvmx3PVccX/sq1d75VqL1klP+zTgKu7jC43e+n6Y1SjylHfF\naOA4v6xueHxhMcxahjjL81nI21fUDSi8WGq0PsrPPvQzDu8/XM5gc3LiSPvmlTcLK/lrTCJPZkZ4\n2bv/v1+HW2RqfnbUgqA4Vq+e/LOf5xxXwjc87loNTrryJG+/0fMGYPahQqUANp62Ea3p+FdY2IeH\nJ49jeHjqNfht87zWWQkKcZGMm+D75/fZzyjfDmSw7B64r/nCLlOX8E93Yc+CRe4tJOvH4LLIGuGF\n7YcoOyJtUs7Fo3f5uVaDx08b5KrtWxgbH5D7tUwikC0TVYkxzjLJ8vrFTfomve+uVo2LTZVE8Dx7\nBvew76p9Xs14is8552VR/yJec8drWnDmzsFsmQ6jHW2ZYHvfinHFpf/YRUVMXfRT31nn3Td8EJ3/\nVONg97FMHZx3fP+rkleNZin7G4W/25CLLZM28Zrn5hiHn444dlBbIuAAzarD0o2YLdMi2mk5dxyu\ny9yDH/d9YQ9aNH7VxMQqiin9+sLu9+v34dsug7cNIsPCwE0D6IKnCtkuAIKM2y93vGeqsLu8fy67\nQ4H3miW1DaZs+v/HfQIKVoBMsqmC/frsGdwzaQOK+wbuQ5ss7HOWz+HEkRMnLJax5lVHnK5Y5F4y\neSPwViz7zvLxH5IjyzxZIeFj/Ei9vrPOh277IPufLxilNyL0GTKDi069yGkT6TJyvrNUbgz65knn\nzbJOYI2MMrS4xZOfFplXhtkyLaIVi1Hy4JLp4RMW9zwpfWniJAIv+cQJPN7z08YD6eOKJWC7BCdF\nXXDNS48jbuPn1au91bFZzuuaSrpncA/7rt4HjTVCzfbLZy6eyfFXHG+TnU3CbJkmUkbGQrNxqSES\njiiTbIM811qrQf3uK1n4SW8j6cd7flrYdpkpMzmPkSnZLknEvX9Z8PsIC7tvXyUJu98uqs+w7fXT\njXvoG/bslTvx/t+3eULYoXphn7l45iSL5Q1PvsGEvQ2xyL1kqorcg0vbixAXWabZDnHXlRSxx0We\nvbVBfjzjKo6MNVIzCma7LJjtLTLya7u4TkiGcfmE4kIWWyaK4Hn9Comtyl4BkAXCCVd174rPTsNs\nmRZRlefushrS5Rx5Uy/jlqO7Xq+cUmfBOy7k2cMFV42CJ+qHZsM/XzteeRHy+dPhNsF+8uKfyy9U\n5kqUvQLNFfWexT2svGKlCXkbY7ZMi3DNpAiLR9FSqy5i5LcJipwv2K59ZB2nn+3CeQM8e+RAedku\nlx1Ed3jC7i+0Ci8ycrHGwpYMTN40Ig/+78C2ben9/OCkH4xnsUTZK80Q9kX9i8Ytljc++UYT9i7B\nIvcmE/eRPctiF8hWMjapjEAwsowqB5BUfiBtTGuuL1i0qzGWxfMXs3/kinExD5IUpWf9pFLUToki\nuHy/2ZURw1aOTXx2Bxa5tzHBSDlpIjY8oQaT88CznC88MVerpVsG4fIDPuGcd7/v49fWWXL50ciQ\nIEPC1geLVWNcyDEwrOz/mLfxRVQ0npZHHpWHH/XaxeWc54ngR+uj43nl/qTnfQP3NUXYtfF1GPg6\nx3AGfZxBH9tqNvE53bDIvQkkReCQvAVbUvQJ2T101+g87pPBlE8cZW18AROh5sGFjJw/efOLLNG0\ny2rOqHbBnPOsfxaj9VH2bNjDkWeOtGTyM6pmeaek5RrZsMi9jYiKwGEiTS4sMFHHRq3gzEPejZWD\nx9V31pFhLzof33u0iEHceF38laO12b+asqtRFqJem6jrDt/UXF/TPYN7xn3y4KrPI894kXkzhF3m\nyKR0RFssZISZ9pF7FStD4xb5uBSLSoq2sj6XtQBVWqQ3eNsgm+/Z3Dgovd9EGuFtj/Sw/tT1qatH\ny3yf0l7r8Ln2DO7xJjtbSQ8csz7bis9WrHo2qsdSIR2p4qOry+RmlBWTtmo06Y81S/qf34eLHTNl\ni7oiNATd3/zCj86bKUJpN9mWb0DRAydeZzXLjXhM3GPIm6edhTyZK1Fkybd29ZmjHo/LYS+c7QIT\nXvwvl8PWy6idt855nFnIcoMIf1LxFwo1i7An37Owh5VbLLfccGNaiHueiC+tzGzeCDJL2l2cLdCs\nhUouaZjjFkzenHSBxfMWc8XZV4xvTxdHkTrkUX3E2WJPDe/hreyjJ3x+91MXZwb809gxXKHmkRv5\nmBbinifiy1sDO4ncqzYz3Giy3Miibh5xon9erc43+COO8EJjUF7BrW0PbOVIWAWTCISjq5asYtcf\n7xp/KksOf5bVsnGPBb9vB788nMliXrhRBBP3AElC4loZMcs4sowrfFNIW2CTZTxptsxJV57E7id3\nu3WahMKsGbN44VMvRD7tUlohy3sQt0jpc9xLL5690uxURABmw4nXml9uVEvXpkKmLfyJOyZuMUve\n1MAkXPuNYMN5AAATPElEQVRMWkxTxTjqO+ssufzFUPPSGHc/kSDsaWLrPz/mlQOIE3Yo55rC7/sa\nGeWE4e9xJ9vGv3p5uvIl+8GXRRZMTkfsO1juHp8W3RtF6Ehxd1116NpfnnZJN5m0G00cvhj39aWf\nI6rfcFu/Rkp9Z52Fn1nAwE0D7D8Y2M0oSQWFaIFXQAV+uBGGFC4dY/NbNiEyMe6suNwMNx4/yr8u\n/jfuagj5J7iPl3Fw0qVUJeoKPE0Pn+FEzmys+DyDPu7636srjdKL1hsypjfTwpYJktfvzLuyM60f\n1zZZbZmRHcVTGDf2bmTL9i2M/44ocM9GuN3LSS/L1gozWh/lZx9q8U5CRG/W3MxVn7bC1IhiWnju\nVU5MZZkkzeOHJ503SjBdz1HfWWfg7y+EOQcKrRjtPy56F6Ok+YW878eewT3s27Iv3QqqkKjl+1FU\nLbhFJpqN6UHXeu5BqvhlDy7wyePvx/WZ1k9f34SVEmzjWx1B6yJ8/vrOOis+/wpkSBi4eQDm5hT2\nhmjFCXuYoSGvQqSPq40QLHM7Xuq2ycIe3rDZdeVnFXM0Qcq2HY3pi1PkLiJnAVcAPcA1qnp56Pll\nwHXAokabj6vqN5P6bJcVqmGylOT1ybOAJul5SG+j2shF37556qqYrPjHP9CP3pAs6sFrdXmNWp2K\nqIDMgEdeewwD95SXW96MdEazZYwoSrNlRKQH2AP8LvAIcA9wvqruDrS5Gvixqm4WkVXAN1V1RVK/\n7SjuVadMQrac76hNlX/7L97HD2Zex9gMis8gKsydOZdzDl/DTUPrcqUkhvP0+xnlAzzISzmIzAbi\nk2gqI1i3vIoa7dAc4bV8eCMKV3Gf6dDX6cADqvpgo+OvAucCwTw6Bf5H4/ujgBZXWZqMSwZLGT6n\nyx9j3Mf64AIcX9T9KP6kK0+aSFucTSHLxUeYwYbTLuIl92zydjAK2EH+OOOycoKv1YLhe7mTiaX7\nk4bWDGHPUVCrUzBhN4rgErm/HThLVT/Q+PndwOtU9eJAm5cD3waOBhYAa1T1R0n9NjNyL3MBUFnn\niTsu/P28T8/nwJED2TsNol453c1v2ZS4qClt/OGJz6JuUBr+UCadQ+CYDclinla/J+/kpE12Gu1A\nmbaMi7h/uNHX50Xk9cCXgFer6lior/XAeoBly5adunfv3oyXFU9axUTXP8C4IloulCHuJ76jzv0r\nLob5gUJWBbJeFs5ZyJa3eBtfuJQjCL5OrfLL/ZvGvFXzeN2u1xXqq5NtGcOIosxsmUeBVwZ+fkXj\nsSAXADcCqOr3gLnAknBHqnq1qvaqau/SpUsdTu1GlFCFM1T8zZPTBNvP/siyuClPRs2k406uI5+Y\niwwJ968agAVP51+Z0xCc5UctZ2TtCL/684mNL8IZN+GsjGYLe1gb/RWfZ+JlsBQVdsOYzrh47vcA\nx4vIq/BE/Z3Au0JtfgH0A18WkRPxxP2JMgcaR1wEGs7qgPRIa2goeuPntGPiskeS2g4NwfFnXskf\n3XExLxSZHA2eTwXu2QC3b+J9NVh38tTzB/FL3d4J3CXQB2wb9hbvPL2t/BK4Gvr+GxzD0bWVU8ZV\n+1l55/RvaGWnMFadEmkYRXFNhTwH+Fu8NMdrVfUyEbkU2K6qtzQyZL4ILMT7u/2oqn47qc+yFjGl\n7W6UxSMtWnLXJc1xZEedi75+Ic/6PnrRFMaf98PIRPpi1LW1csVnz8Iejjx7hDnL5vDJvcfxhtpL\nxzOPLBvEMLLT9StUswh3X58Xkcd5r83IllnzxTex9dF/9X4oYxYysIl02FduhZiHJz/nLJ/DcZcd\nx+afvdQmIQ2jRLpa3F3EOM6WKatGjAv1nXUu/NrFHJCGxVGC9TK7ZzbXvu3aSRtgrJFRPr38QZ7f\ne7BppW7DvzX/zQzexpuA9LLFNhlpGPnpanH3kZCSxa2WDLfzCYtQUVvGp/7Pn+VD/3YJ++dR2Hbp\nmdHDdX9w3XiE/idM7CYUjJarEPVjNh4TOal6GPhnvN2EgoINbitrTdwNIz9lLmJqe+J2M4oT9SRh\nCWbLZIncB298H1ftvo7x3M/57sdOIpA83n9cP399/V/zdO1ptrGNO5ks4lVG6Yv6F7Fy08rIfPKh\nIW9BQxbCn6LSFksZhlGMjovc0xaoZCFqKzt/ss8l8+X4tYGyulA4Sv/c9Z+j96FeJNBR1QuFos4T\nVeo2iMt7kCTa7R6520Sv0c5MG1smamWn/3McRa2D+t1XMnjrh/jv+UdyK2//jn4uvv1ijjpw1MS5\nK5TxoIUTFPJgHZY40haIxb0HcbSLuMddV7uMzzCi6HpxD9c9j6tQGPaCs2TL+KxeDRf+XZ1Lbv8I\ne597rNFR9jH/ya1/wrnbzx0X8SrF3EeBGQuEE646YVzAs0amrrXsXUQxOOHdauLGa+JutDNdX899\neHjyQpLwopJgnfEgUXXRo1ZrgpeTPvfSOXynTxi4aYC9Bx7LNHvZv6Off/rLf+LOoTu5c+hO3rb9\nbcxgBtL4VwUa+HqamVzGiVO2gytTWJPegyj8OvlB2kHoy6rdbxjtQsdG7mn11ePSG/3no/rwj3sZ\nj/HbZ7+Gm09/3EnI+3f08+FbP8y8F+ZNHWfTkhOZVFSrSI0cqC5FtOgWgkVxuS6L3I12pittmTTR\nDrZz2XQj+P3gbYNctX0LY5MaR4+jf0c/H9j6AV7yy5c03WJx9crLFKi4vlxvHGmC2ioxNVvG6ES6\nUtyDZJlADUdlUxhYA7+2NVbMw145VC/m/qXNith4woUyMz7KEsHgPEk7lM8tetMyjFYwrcU9GLH7\nf6RhIXnlyZ/l2bMu4b/8fPTG8Z+7zktHnHK+Jtkr/iUd3UhHbAcRLCurJG6hU6siZRNxoxPpenHv\n6/O+XDdlWHP9GrY+tHXC2wguFtrRz0dv/iizmAU0z2LxOUAPf8PK8aJa0L52QZGbTVxd/Xa4LsPo\nFLo+W8a1NO/glrciQ8LWB7d6D4hns9wxfMd4FsslN1/CbGZXlsUSzGDxfz66fxFnaB9n0sc5+sZJ\nwg7tm60RlVmUttgrra6+lc81jPLp2Mg9zZYZvPF9vOrjJ0VaLFBddB7+YHB4lnD5oRPYyuSJz6TJ\nxE7xgrPmuFuEbhjF6crIPSkXuZ9Rbp/5Le7kLu6Su/jD//Xe8aX8Uf9Kp8crtPWdmreL0Bnq/b/m\nhdW8oTYh7P6OR34EG76OtOtvJyziNoz2paMKhwUj1zUyyicX3M/Ys2MwDKsBOTyn8jEEA8+5jZrl\nkxYIhcYb9qd93zkp4u0k0Qx+YkorBtZJ12UYnU5H2TKj9VHue+996BFvzM2e+ATYziKerb0m1+bZ\n3TaZmGWLQcMwyqHrsmVG66PsHtjdtPxygDHgFo7hC0yUvS17C75289GzkKdomGEYxei6eu47P7KT\nmSUOV0PFdA/N6+E3vriSlw28dIpgXcFE5J2XuGM7VdjDmOViGO1Fx0yo9jxWbKja+Dfe36IezsSb\n+DxD+/i9596YWPa26KKhbhHxuEltwzDai44R99GjHs/UXgP/xkR5tPdYzuQM+hpZLG966k3jbeME\nK1hZMqs4FxXzdr0ZxOW5t+t4DWO60jGe+zsH3sn6+npmRNyPNDTtOcYYtxz///jCzz4SWzQMov3u\nsrzjov10gofdCWM0jG6j6/Lc3/qxt/LZ8z7LITk0KSpXlK/3fp0za2dy5qfO5Myz/5E19POFn32E\n1auTc+Mt2iyG+eyG0caoaku+Tj31VM3KyI4R7RnuUYZQaihDaM9wj47sGBlvA6q1mvcVxjcRkog6\nzvX5Wm3iHMGv4DFFj88zZsMwugdguzpobMfYMq4kbVIRV5UwT/9521Wxz6jZI4Yxfeg6W8YV3yoI\nrgz1rRmfdi3KZRiGURZdJ+5Rgl00wyPP/prhvVqLHF/WmJL6Mgyju3CyZUTkLLy1PD3ANap6eUSb\nd+CVVlHgP1T1XUl9VmHLVLk/ZnDTjzJsnWZsg9es4w3DaB6llR8QkR5gD/C7wCPAPcD5qro70OZ4\n4EbgTFV9SkReoqqJielVee4TY4oW0bxZMuGSwmWPy4WydkQq+3jDMJpHmZ776cADqvqgqr4AfBU4\nN9TmQuBKVX0KIE3YW0kWKyaKIul/4X1cs1opcbtO5RlTmbaOYRjth4u4Hwv8Z+DnRxqPBVkJrBSR\n74rI9xs2TssoKqLgCWlYAP3Hi3jbVazubKexGIbRHpQ1oToTOB7oA84Hvigii8KNRGS9iGwXke1P\nPPFESaeeSlnCFe6nSF95sQjbMIw8uIj7o8ArAz+/ovFYkEeAW1T1kKo+hOfRHx/uSFWvVtVeVe1d\nunRp3jGXji+UZQtpkWyYYB9VR9i20tQwug+XCdWZeGLdjyfq9wDvUtVdgTZn4U2yvldElgA/Bl6j\nqvvj+q16QtXHZQLVZR/T4A3AlbInKm3i0zCM0iZUVfUwcDHwLeA+4EZV3SUil4rI7zeafQvYLyK7\ngbuAjyQJexXEiW5ZEW471KKxCNswDFecPHdV/aaqrlTVX1PVyxqPfUpVb2l8r6r6YVVdpaonq+pX\nqxx0FHGZJHGkWTB5hbRKj7zVNxfDMDqHrqktU8SyqMruMBvFMIyymRa1ZSyTxDAMI5qOF/csmSRV\nLExKwjxywzBaxbSyZcwmMQyj05kWtkwQi5INwzAm6BpxT7JizJc3DGO60TXiHkfZKzztpmAYRifQ\n9eJeNlnz6Q3DMFpBx4p7ngjafHnDMKYLHSvuSRF02aUIzLc3DKPT6NhUyKS0xipTHi2d0jCMVtKV\nqZDdFkF36rgNw2h/Ok7c4zbOaJbwl+nb2+SsYRhV0ZG2jC/grbBlyqRTxmkYRvvQdbZM1H6mnWjJ\ndJu1ZBhGe9KxkXvSsNthYw0XLHI3DCMrXRe5pxEU804QdsMwjCrpSHGPmtTsxMlJW1RlGEZVdKS4\nhyPzTo3UO3XchmG0Px0p7j7+5KQftdvkpGEYhkdHTqiGSUuNNAzD6Ba6fkK1W1IjDcMwqqCjxT1Y\np71WK1an3TAMo5voWHEPY6JuGIYxQVeIu6UUGoZhTKYrxN2idsMwjMl0hbgbhmEYk3ESdxE5S0R+\nKiIPiMjHE9qtFREVkdQ0HcMwDKM6UsVdRHqAK4GzgVXA+SKyKqLdi4APAT8oe5CGYRhGNlwi99OB\nB1T1QVV9AfgqcG5Eu08Dfwk8X+L4DMMwjBy4iPuxwH8Gfn6k8dg4IvJa4JWqeluJYzMMwzByMrNo\nByIyA/gb4H0ObdcD6xs/PiMiP8152iXAkzmP7VTsmqcHds3TgyLXvNylkYu4Pwq8MvDzKxqP+bwI\neDWwTbxaAC8DbhGR31fVScVjVPVq4GqXgSUhIttdait0E3bN0wO75ulBM67ZxZa5BzheRF4lIrOB\ndwK3+E+q6i9VdYmqrlDVFcD3gSnCbhiGYTSPVHFX1cPAxcC3gPuAG1V1l4hcKiK/X/UADcMwjOw4\nee6q+k3gm6HHPhXTtq/4sFIpbO10IHbN0wO75ulB5dfcsnruhmEYRnVY+QHDMIwupK3FPa3sgYjM\nEZF/aDz/AxFZ0fxRlovDNX9YRHaLyA4R2SoiTmlR7cx0LG/hcs0i8o7Ge71LRP6+2WMsG4ff7WUi\ncpeI/Ljx+31OK8ZZFiJyrYg8LiI/iXleROQLjddjR2O9UHmoalt+AT3Az4HjgNnAfwCrQm0GgS2N\n798J/EOrx92Eaz4DmN/4fuN0uOZGuxcBd+NlY/W2etxNeJ+PB34MHN34+SWtHncTrvlqYGPj+1XA\nw60ed8FrfhPwWuAnMc+fA9wOCPBbwA/KPH87R+4uZQ/OBa5rfP81oF8kuPFex5F6zap6l6o+1/jx\n+3jrDjqZ6VjewuWaLwSuVNWnAFT18SaPsWxcrlmB/9H4/ihgXxPHVzqqejfwXwlNzgWuV4/vA4tE\n5OVlnb+dxT217EGwjXopm78EFjdldNXgcs1BLsC783cy07G8hcv7vBJYKSLfFZHvi8hZTRtdNbhc\n8xAwICKP4GXnfbA5Q2sZWf/eM1G4/IDRGkRkAOgFVrd6LFWSpbxFlzETz5rpw/t0dreInKyqT7d0\nVNVyPvBlVf28iLweuEFEXq2qY60eWCfSzpF7WtmDSW1EZCbeR7n9TRldNbhcMyKyBrgEbyXwwSaN\nrSqylLd4GM+bvKXDJ1Vd3udHgFtU9ZCqPgTswRP7TsXlmi8AbgRQ1e8Bc/FqsHQrTn/veWlncU8s\ne9DgFuC9je/fDtypjZmKDiX1mkXkN4Gr8IS9031YmJ7lLVx+t7+OF7UjIkvwbJoHmznIknG55l8A\n/QAiciKeuD/R1FE2l1uA9zSyZn4L+KWqPlZa762eUU6ZbT4HL2L5OXBJ47FL8f64wXvz/xF4APgh\ncFyrx9yEa74DGAXubXzd0uoxV33Nobbb6PBsGcf3WfDsqN3ATuCdrR5zE655FfBdvEyae4Hfa/WY\nC17vV4DHgEN4n8QuADYAGwLv8ZWN12Nn2b/XtkLVMAyjC2lnW8YwDMPIiYm7YRhGF2LibhiG0YWY\nuBuGYXQhJu6GYRhdiIm7YRhGF2LibhiG0YWYuBuGYXQh/x8beVQgRWLwbgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe586d7b6a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "STEP = int(X.shape[0] / 500)\n",
    "if n == 1:\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(111)\n",
    "    ax1.plot(X[::STEP, 1], y[::STEP], 'b+') # blue - data\n",
    "    ax1.plot(X[::STEP, 1], X.dot(theta_true)[::STEP], 'r*') # red - true\n",
    "    ax1.plot(X[::STEP, 1], X.dot(theta)[::STEP], 'go') # green - batch\n",
    "    ax1.plot(X[::STEP, 1], X.dot(theta_lms)[::STEP], 'mo') # magenta - pure LMS\n",
    "else:\n",
    "    print(\"Plot is multidimensional; I live in Flatland, so I don't do that.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Exercise 2** (_ungraded_, optional). We said previously that, in practice, you would probably do some sort of _hybrid_ scheme that mixes full batch updates (possibly only initially) and incremental updates. Implement such a scheme and describe what you observe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "nbgrader": {
     "grade": false,
     "grade_id": "hybrid",
     "locked": false,
     "schema_version": 1,
     "solution": true
    },
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Setup problem and compute the batch solution\n",
    "m = 100000\n",
    "n = 1\n",
    "theta_true = generate_model(n)\n",
    "(X, y) = generate_data(m, theta_true, sigma=0.1)\n",
    "theta_batch = estimate_coeffs(X, y)\n",
    "\n",
    "# Your turn, below: Implement a hybrid batch-LMS solution\n",
    "# assuming you observe the first few data points all at\n",
    "# once, and then see the remaining points one at a time.\n",
    "\n",
    "#\n",
    "# YOUR CODE HERE\n",
    "#\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Fin!** If you've gotten this far without errors, your notebook is ready to submit."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Create Assignment",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}