{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "4e0f42c6439b082fc537fe0ec1f0d801",
     "grade": false,
     "grade_id": "cell-68b905920fbe16ce",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "# Part 1: Linear Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "b669f6696864bf529ca9f9af6217d0e4",
     "grade": false,
     "grade_id": "cell-1c5581fc981be7ed",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# Execute this code block to install dependencies when running on colab\n",
    "try:\n",
    "    import torch\n",
    "except:\n",
    "    from os.path import exists\n",
    "    from wheel.pep425tags import get_abbr_impl, get_impl_ver, get_abi_tag\n",
    "    platform = '{}{}-{}'.format(get_abbr_impl(), get_impl_ver(), get_abi_tag())\n",
    "    cuda_output = !ldconfig -p|grep cudart.so|sed -e 's/.*\\.\\([0-9]*\\)\\.\\([0-9]*\\)$/cu\\1\\2/'\n",
    "    accelerator = cuda_output[0] if exists('/dev/nvidia0') else 'cpu'\n",
    "\n",
    "    !pip install -q http://download.pytorch.org/whl/{accelerator}/torch-1.0.0-{platform}-linux_x86_64.whl torchvision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "cb81ed04ab5982d538771590f45802ce",
     "grade": false,
     "grade_id": "cell-8128a87aa950cfd5",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Getting started \n",
    "\n",
    "At its heart, PyTorch is just a library for manipulating tensors. We're going to start learning how to use \n",
    "PyTorch by looking at how we can implement simple linear regression. \n",
    "\n",
    "Code speaks better than words, so lets start by looking at a bit of pytorch code to generate some 2d data to regress:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "\n",
    "# Generate some data points on a straight line perturbed with Gaussian noise\n",
    "N = 1000 # number of points\n",
    "theta_true = torch.Tensor([[1.5], [2.0]]) # true parameters of the line\n",
    "\n",
    "X = torch.rand(N, 2) \n",
    "X[:, 1] = 1.0\n",
    "y = X @ theta_true + 0.1 * torch.randn(N, 1) # Note that just like in numpy '@' represents matrix multiplication and A@B is equivalent to torch.mm(A, B) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "15906b16ca0952809a0d0f6821439957",
     "grade": false,
     "grade_id": "cell-67f4ae4fafb20c75",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "The above code generates $(x,y)$ data according to $y = 1.5x + 2$, with the $x$'s chosen from a uniform distribution. The $y$'s are additionally purturbed by adding an amount $0.1z$, where $z\\sim \\mathcal{N}(0,1)$ is a sample from a standard normal distribution. \n",
    "\n",
    "Note that we represent our $x$'s as a two-dimensional (row) vector with a 1 in the second element so that the offset can be rolled into the matrix multiplication for efficiency:\n",
    "\n",
    "\\begin{align}\n",
    "    y &= \\mathbf{X}\\begin{bmatrix}\n",
    "           1.5 \\\\\n",
    "           2\n",
    "         \\end{bmatrix}\n",
    "  \\end{align}\n",
    "\n",
    "Let's use `matplotlib` to draw a scatter so we can be sure of what our data looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "77838c69c9e63098d1e93e6ea06ca513",
     "grade": false,
     "grade_id": "cell-36f16b5ab42ab959",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.scatter(X[:,0].numpy(), y.numpy())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "dfce8172c55747e4c8aa47d6bd3ca291",
     "grade": false,
     "grade_id": "cell-26114c36c1f243e4",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "__Make sure you understand how the code above is generating data; feel free to change the parameters to see what effect they have.__\n",
    "\n",
    "Now, lets consider the situation where we have been given the tensors $X$ and $y$ and wish to compute the regression parameters. Our model looks like $\\mathbf{y} = \\mathbf{X\\theta}$, and we wish to recover the parameters $\\theta$. \n",
    "\n",
    "As the problem is both overcomplete (only two data pairs are required to find $\\theta$), and the data is noisy, we can use the Moore-Penrose Pseudoinverse to find the least-squares solution to $\\theta$: $\\theta = \\mathbf{X^+y}$. PyTorch has a built-in pseudoinverse method (`pinverse`) that can do all the work for us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "3226facb17e2a669d112c2b1700a8daa",
     "grade": false,
     "grade_id": "cell-3ae7c2a27cf28ee6",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "# direct solution using moore-penrose pseudo inverse\n",
    "X_inv = torch.pinverse(X)\n",
    "theta_pinv = torch.mm(X_inv, y)\n",
    "print(theta_pinv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "0c6a1fb676027787f94fc1b27eb815ed",
     "grade": false,
     "grade_id": "cell-2f7080b24616a7d1",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Running the above code should give you a solution vector for $\\theta$ that is very similar to the true parameter vector (`theta_true`). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "c2d7e31bd39f66e7c4122c6fffc42006",
     "grade": false,
     "grade_id": "cell-ca707e4c78ca3acc",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Exercise: computing the pseudoinverse from the Singular Value Decomposition.\n",
    "\n",
    "The standard way of computing the pseudoinverse is by using the Singular Value Decomposition (SVD). The SVD is defined as: $\\mathbf{X} = \\mathbf{U}\\Sigma\\mathbf{V}^\\top$. The pseudoinverse is thus $\\mathbf{X}^+ = \\mathbf{V}\\Sigma^{-1}\\mathbf{U}^\\top$ where $\\Sigma^{-1}$ is a diagonal matrix in which the reciprocal of the corresponding non-zero elements in $\\Sigma$ has been taken.\n",
    "\n",
    "__Use the code block below to compute the parameter vector using the SVD directly rather than the through the `pinverse` method.__ You need to store your manually computed pseudoinverse in `X_inv_svd`. Useful methods will be `torch.svd()` to compute the SVD, `[Tensor].t()` to transpose a matrix and `torch.diag()` to form a diagonal matrix from a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "96170c2a936223004245b58ae10a2b62",
     "grade": false,
     "grade_id": "cell-4d422a83842ebab7",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()\n",
    "\n",
    "theta_pinv_svd = torch.mm(X_inv_svd, y)\n",
    "print(theta_pinv_svd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "d9743a1c36a687bf2adb629500247a05",
     "grade": true,
     "grade_id": "cell-4e7e9488b81a58b8",
     "locked": true,
     "points": 2,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "assert(torch.all(torch.lt(torch.abs(torch.add(theta_pinv, -theta_pinv_svd)), 1e-6)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient based Linear Regression\n",
    "\n",
    "Fundamentally, with linear regression we are trying to find a solution vector, $theta$ that minimises $f(\\theta) = 0.5\\|\\mathbf{X}\\theta - \\mathbf{y}\\|_2^2$. \n",
    "\n",
    "We've already seen how this can be minimised directly using the pseudoinverse, but it could also be minimised by using gradient descent: $\\theta \\gets \\theta - \\alpha f'(\\theta)$. (_Interesting aside_: SVD (and thus the pseudoinverse) can also be solved using gradient methods - in fact this becomes the only practical way for really large matrices.).\n",
    "\n",
    "__Use the following block to derive and write down the gradient, $f'(\\theta)$, of $f(\\theta)$__. Note that you can insert latex code by wrapping expressions in dollar symbols."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "0197441552c3f53e1fb98b1ee6d232b8",
     "grade": true,
     "grade_id": "cell-210cc7d9ab3905e5",
     "locked": false,
     "points": 2,
     "schema_version": 1,
     "solution": true
    }
   },
   "source": [
    "YOUR ANSWER HERE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "4d1c00f47e8080e27a49ba3999c0b342",
     "grade": false,
     "grade_id": "cell-d8813ea915bc08f7",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "__Now complete the following code block to implement your gradient as pytorch code:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "65ea0d0a4244f58b618bd1eb0272fa3e",
     "grade": false,
     "grade_id": "cell-640dcd4113de31d6",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "def linear_regression_loss_grad(theta, X, y):\n",
    "    # theta, X and y have the same shape as used previously\n",
    "    # YOUR CODE HERE\n",
    "    raise NotImplementedError()\n",
    "    return grad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "1ff96ad5def7b717aea28ab5ccaedc04",
     "grade": true,
     "grade_id": "cell-9cdc211cfd5cab66",
     "locked": true,
     "points": 2,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "assert(linear_regression_loss_grad(torch.zeros(2,1), X, y).shape == (2,1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "646a7eb8b434567047a8278198a74ffa",
     "grade": false,
     "grade_id": "cell-298e11818c4e9b65",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Now we can plug that gradient function into a basic gradient descent solver and check that the solution is close to what we get with the pseudoinverse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "0ee711fe83b16cd14129aec603c141e9",
     "grade": false,
     "grade_id": "cell-cef34abd27fe76a8",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "alpha = 0.001\n",
    "theta = torch.Tensor([[0], [0]])\n",
    "for e in range(0, 200):\n",
    "    gr = linear_regression_loss_grad(theta, X, y)\n",
    "    theta -= alpha * gr\n",
    "\n",
    "print(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "8a785d76923ccbb7fa7b953a67377385",
     "grade": false,
     "grade_id": "cell-47be64696ed6859b",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "## Real data\n",
    "\n",
    "Doing linear regression on synthetic data is a great way to understand how PyTorch works, but it isn't quite as satisfying as working with a real dataset. Let's now apply or understanding of computing linear regression parameters to a dataset of house prices in Boston.\n",
    "\n",
    "We'll load the dataset using scikit-learn and perform some manipulations in the following code block:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "68160322d431987aaa1330088bd7efa5",
     "grade": false,
     "grade_id": "cell-3d4dfae2229c46df",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.datasets import load_boston\n",
    "\n",
    "X, y = tuple(torch.Tensor(z) for z in load_boston(True)) #convert to pytorch Tensors\n",
    "X = X[:, [2,5]] # We're just going to use features 2 and 5, rather than using all of of them\n",
    "X = torch.cat((X, torch.ones((X.shape[0], 1))), 1) # append a column of 1's to the X's\n",
    "y = y.reshape(-1, 1) # reshape y into a column vector\n",
    "print('X:', X.shape)\n",
    "print('y:', y.shape)\n",
    "\n",
    "# We're also going to break the data into a training set for computing the regression parameters\n",
    "# and a test set to evaluate the predictive ability of those parameters\n",
    "perm = torch.randperm(y.shape[0])\n",
    "X_train = X[perm[0:253], :]\n",
    "y_train = y[perm[0:253]]\n",
    "X_test = X[perm[253:], :]\n",
    "y_test = y[perm[253:]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Use the following code block to compute the regression parameters using the training data in the variable `theta` by solving using the pseudoinverse directly:__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "4f967bc6cfdf19e504fbc602fde368a7",
     "grade": false,
     "grade_id": "cell-b7854194044222ed",
     "locked": false,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# compute the regression parameters in variable theta\n",
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now print out the error achieved on the test set, as well as the parameter vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "724495cb56038fbab28377c134c95d92",
     "grade": true,
     "grade_id": "cell-a97d86a0f20894f5",
     "locked": true,
     "points": 1,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "assert(theta.shape == (3,1))\n",
    "\n",
    "print(\"Theta: \", theta.t())\n",
    "print(\"MSE of test data: \", torch.nn.functional.mse_loss(X_test @ theta, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "d794e7faf7455d8154d6fcab85e15e4d",
     "grade": false,
     "grade_id": "cell-f1a7b0a932a8ba2e",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Now let's try using gradient descent:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 0.00001\n",
    "theta_gd = torch.rand((X_train.shape[1], 1))\n",
    "for e in range(0, 10000):\n",
    "    gr = linear_regression_loss_grad(theta_gd, X_train, y_train)\n",
    "    theta_gd -= alpha * gr\n",
    "\n",
    "print(\"Gradient Descent Theta: \", theta_gd.t())\n",
    "print(\"MSE of test data: \", torch.nn.functional.mse_loss(X_test @ theta_gd, y_test))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "f2b9cd4ed786c8a74d31e22d7f7cd7cc",
     "grade": false,
     "grade_id": "cell-7afbfb1f2aaaa9a4",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "__Use the following block to note down any observations you can make about the choice of learning rate and number of iterations in the above code. What factors do you think influence the choice?__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "checksum": "487706e65407095404aaaf95ab1dc586",
     "grade": true,
     "grade_id": "cell-1e987d4019c368ba",
     "locked": false,
     "points": 3,
     "schema_version": 1,
     "solution": true
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE\n",
    "raise NotImplementedError()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "138abf4c0e59e16e5be8c847dff30f8d",
     "grade": false,
     "grade_id": "cell-809a19fe6970a99b",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "source": [
    "Finally, just so we can visualise what our model has learned, we can plot the predicted house prices (from both the direct solution and from gradient descent) along with the true value for each of the houses in the test set (ordered by increasing true value):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "checksum": "f486683e0b3e4d14da067504831f6cd3",
     "grade": false,
     "grade_id": "cell-64d6a9203da564f1",
     "locked": true,
     "schema_version": 1,
     "solution": false
    }
   },
   "outputs": [],
   "source": [
    "perm = torch.argsort(y_test, dim=0)\n",
    "plt.plot(y_test[perm[:,0]].numpy(), '.', label='True Prices')\n",
    "plt.plot((X_test[perm[:,0]] @ theta).numpy(), '.', label='Predicted (pinv)')\n",
    "plt.plot((X_test[perm[:,0]] @ theta_gd).numpy(), '.', label='Predicted (G.D.)')\n",
    "plt.xlabel('House Number')\n",
    "plt.ylabel('House Price ($,000s)')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}