{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lab 3: Linear Regression\n",
    "\n",
    "This exercise takes you through the fundamental linear regression model from a number of different angles. First we consider an analytic analysis and then we continue to consider how to solve the same problem using numerical methods. This lab also serves and an introduction to the tensorflow toolkit which will be useful for a variety of machine learning tasks in the future and is used by Google to solve massive machine learning problems on their clusters.\n",
    "\n",
    "The linear regression model forms the basis for a whole host of models - if you are comfortable with the fundamental approaches we take here, there will be a whole range of extensions, advances and applications available to you in the future."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Marking and Submission\n",
    "\n",
    "These lab exercises are marked, and contribute to your final grade. For this lab exercise there are 6 places where you are expected to enter your own code, for 15 marks overall. Every place you have to add code is indicated by\n",
    "\n",
    "`# Add your code here ..`\n",
    "\n",
    "`# ***********************************************************`\n",
    "\n",
    "with instructions above the code block.\n",
    "\n",
    "Please submit your completed workbook using Moodle before 5pm on the 14th November 2018. The workbook you submit must be an `.ipynb` file, which is saved into the directory you're running Jupyter; alternatively you can download it from the menu above using `File -> Download As -> Notebook (.ipynb)`. Remember to save your work regularly (Save and checkpoint in the File menu, the icon of a floppy disk, or Ctrl-S); the version you submit should have all code blocks showing the results (if any) of execution below them. You will normally receive feedback within a week.\n",
    "\n",
    "## IMPORTANT!\n",
    "\n",
    "- Please make sure your notebook runs all the way through without errors with the cells executed in the correct order before saving and submitting. This can be performed by selecting 'Restart & Run All' from the 'Kernel' menu at the top.\n",
    "- All answers should be provide within functions with specific names and arguments.\n",
    "- **Please do not change the names or signatures of the functions that you write your answers inside!** This will break the script we use for marking.\n",
    "- **Please do not move cells around or copy and paste the notebook into another notebook!** Again, this will break the marking script."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "# A new one for this lab!\n",
    "import tensorflow as tf\n",
    "\n",
    "from sklearn import datasets as ds\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boston = ds.load_boston()\n",
    "\n",
    "x_raw = boston.data[:,np.argwhere(boston.feature_names == 'RM')[0,0]]\n",
    "y_raw = boston.target"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boston House Prices Dataset\n",
    "A descroption of the dataset used is provided here.\n",
    "\n",
    "Data Set Characteristics:  \n",
    "    :Number of Instances: 506 \n",
    "\n",
    "    :Number of Attributes: 13 numeric/categorical predictive\n",
    "    \n",
    "    :Median Value (attribute 14) is usually the target\n",
    "\n",
    "    :Attribute Information (in order):\n",
    "        - CRIM     per capita crime rate by town\n",
    "        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.\n",
    "        - INDUS    proportion of non-retail business acres per town\n",
    "        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)\n",
    "        - NOX      nitric oxides concentration (parts per 10 million)\n",
    "        - RM       average number of rooms per dwelling\n",
    "        - AGE      proportion of owner-occupied units built prior to 1940\n",
    "        - DIS      weighted distances to five Boston employment centres\n",
    "        - RAD      index of accessibility to radial highways\n",
    "        - TAX      full-value property-tax rate per $10,000\n",
    "        - PTRATIO  pupil-teacher ratio by town\n",
    "        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town\n",
    "        - LSTAT    % lower status of the population\n",
    "        - MEDV     Median value of owner-occupied homes in $1000's\n",
    "\n",
    "    :Missing Attribute Values: None\n",
    "\n",
    "    :Creator: Harrison, D. and Rubinfeld, D.L.\n",
    "\n",
    "This is a copy of UCI ML housing dataset.\n",
    "http://archive.ics.uci.edu/ml/datasets/Housing\n",
    "\n",
    "\n",
    "This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.\n",
    "\n",
    "The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic\n",
    "prices and the demand for clean air', J. Environ. Economics & Management,\n",
    "vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics\n",
    "...', Wiley, 1980.   N.B. Various transformations are used in the table on\n",
    "pages 244-261 of the latter.\n",
    "\n",
    "The Boston house-price data has been used in many machine learning papers that address regression\n",
    "problems.   \n",
    "     \n",
    "### References\n",
    "\n",
    "   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.\n",
    "   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.\n",
    "   - many more! (see http://archive.ics.uci.edu/ml/datasets/Housing)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing the data\n",
    "\n",
    "We are going to look at the relationship between the \"average number of rooms per dwelling\" and median house price in the Boston dataset. First let us partition the data into a training and test split. We are going for 60% training and 40% testing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "total_count = x_raw.shape[0]\n",
    "\n",
    "split = int(total_count * 0.6)\n",
    "\n",
    "# Shuffle the data to avoid any ordering bias..\n",
    "np.random.seed(0)\n",
    "shuffle = np.random.permutation(total_count)\n",
    "\n",
    "x = x_raw[shuffle]\n",
    "y = y_raw[shuffle]\n",
    "\n",
    "x_train_unnormalised = x[:split]\n",
    "y_train_unnormalised = y[:split]\n",
    "\n",
    "x_test_unnormalised = x[split:]\n",
    "y_test_unnormalised = y[split:]\n",
    "\n",
    "print('Training set size:', x_train_unnormalised.shape[0])\n",
    "print('Test set size:', x_test_unnormalised.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Visualisation\n",
    "\n",
    "To allow for easy visualisation as you progress through the task we are using a single dimensional data set. Both the input $x$ and output $y$ are scalars so we can plot them on a standard scatter plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You can use this function to plot the data and then add your own plots on top..\n",
    "def plot_data(x, y):\n",
    "    plt.figure(figsize=[10,8])\n",
    "    plt.plot(x, y, 'b+')\n",
    "    plt.grid(True)\n",
    "    plt.xlabel('Average number rooms per dwelling')\n",
    "    plt.ylabel('Mean value of home')\n",
    "\n",
    "plot_data(x_train_unnormalised, y_train_unnormalised)\n",
    "plt.title('Plot of the Training Data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1 (2 marks):\n",
    "\n",
    "Write a function that normalises a vector of values. It should output a corresponding vector where the values have a mean of zero and a standard deviation of 1. Note that you should only perform an affine transformation of the data (i.e. a linear scaling and a fixed offset). This means that you must find $a$ and $b$ for $u_i = a \\, v_i + b$ where $\\{u_i\\}$ are the input data and $\\{v_i\\}$ are the normalised output data.\n",
    "\n",
    "Your function should return the normalising constants as well as the normalised data.\n",
    "\n",
    "Write a second function that removes the normalisation and returns the data to its original values.\n",
    "\n",
    "Check that passing both `x_train` and `y_train` through both functions returns the vectors to their original values.\n",
    "\n",
    "*Hints:*\n",
    "- You might want to look at `np.all()` for the Boolean check that they return to their values.\n",
    "- When checking that floating point values are equal up to nummerical precision (e.g. rounding errors in the computations) you can use the `np.isclose()` function.\n",
    "- You can use the `assert()` command to guarantee that a statement is `True` before the program continues.\n",
    "\n",
    "*Points to consider:* \n",
    "- Why might it be sensible to normalise the data in the fashion described?\n",
    "- Considering that we are about to perform Linear Regression, why might we not want to perform a more involved normalisation process?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def normalise_data(x_unnormalised):\n",
    "    # Add your code here..\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    return x_normalised, a, b\n",
    "\n",
    "def unnormalise_data(x_normalised, a, b):\n",
    "    # Add your code here..\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    return x_unnormalised\n",
    "\n",
    "try:\n",
    "    x_train, x_norm_a, x_norm_b = normalise_data(x_train_unnormalised)\n",
    "    y_train, y_norm_a, y_norm_b = normalise_data(y_train_unnormalised)\n",
    "\n",
    "    x_test, _, _ = normalise_data(x_test_unnormalised)\n",
    "    y_test, _, _ = normalise_data(y_test_unnormalised)\n",
    "except Exception as err:\n",
    "    print('Error during normalisation functions:', err)\n",
    "\n",
    "def checking_function(normalise_data, unnormalise_data, data_to_check):\n",
    "    passes_check = False\n",
    "    \n",
    "    est_normalised_data, est_a, est_b = normalise_data(data_to_check)\n",
    "    est_unnormalised_data = unnormalise_data(est_normalised_data, est_a, est_b)\n",
    "    \n",
    "    # Add your code here to check that the unnormaliseding the \n",
    "    # training data returns to their original values and update \n",
    "    # the passes_check boolean..\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    return passes_check\n",
    "\n",
    "try:\n",
    "    if (checking_function(normalise_data, unnormalise_data, x_train_unnormalised) == True) and \\\n",
    "       (checking_function(normalise_data, unnormalise_data, y_train_unnormalised) == True):\n",
    "        print('Passes checking function :)')\n",
    "    else:\n",
    "        print('Failed to pass the checking function :(')\n",
    "except Exception as err:\n",
    "    print('Error during checking function:', err)\n",
    "\n",
    "# Plot the data to make sure they are normalised..\n",
    "try:\n",
    "    plot_data(x_train, y_train)\n",
    "    plt.title('Plot of the (Normalised) Training Data')\n",
    "except Exception as err:\n",
    "    print('Error during ploting functions:', err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Linear Regression Model\n",
    "\n",
    "In linear regression we are trying to fit a linear model to the data of the form\n",
    "\n",
    "\\begin{align}\n",
    "y &= w x + c\n",
    "\\end{align}\n",
    "\n",
    "where $w$ and $c$ are parameters to be learned that take the input data $x$ to the output data $y$. Once this model has been learned, we can use the parameters to predict the values of the output that would correspond to new values of the input.\n",
    "\n",
    "In order to determine the parameters, we need an objective function that we seek to optimise: this function returns a scalar value for all possible parameter values and we seek to change the parameters until the best scalar value is obtained.\n",
    "\n",
    "For linear regression, we usually take the objective as one which minimises the squared error; this is known as a linear least squares problem.\n",
    "\n",
    "*Aside: Think about what this means in terms of a model for the data when you have $y = f(x) + \\eta$ with $f(x)$ as a linear function $f(x) = w x + c$ and $\\eta$ as iid Gaussian noise.*\n",
    "\n",
    "Therefore our objective is given by the sum of squared differences between the true value of $y_i$ and the value estimated by our model $w x_i + c$.\n",
    "\n",
    "\\begin{align}\n",
    "E(w,c) &= \\sum_{i=0}^{N-1} \\big(y_i - f(x_i) \\big)^2 \\\\\n",
    "    &= \\sum_{i=0}^{N-1} \\big(y_i - (w x_i + c) \\big)^2 \\\\\n",
    "    &= \\sum_{i=0}^{N-1} \\big(y_i - w x_i - c \\big)^2\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2 (1 mark)\n",
    "\n",
    "Write a function that calculates the least squared error on the training data for a particular value of the parameters $w$ and $c$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def least_squares_error(x, y, w, c):\n",
    "    # Add code to calcuate the squared_error = E(w,c)..\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    return squared_error\n",
    "\n",
    "try:\n",
    "    print('Squared error for w = 1.5, c = 0.5 is ', \n",
    "          least_squares_error(x_train, y_train, w=1.5, c=0.5))\n",
    "except Exception as err:\n",
    "    print('Error during least squares calculation:', err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytic Solution\n",
    "\n",
    "In the case of linear regression, we can find an analytic solution to this problem by finding stationary point of the objective function. We do this by evaluating the partial derivatives of the objective wrt each parameter in turn and setting them to zero. If we can then find a solution to these simultaneous equations, we have found an optimal setting for the parameters.\n",
    "\n",
    "For $w$ we have:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial E}{\\partial w} \n",
    "    &= \\frac{\\partial}{\\partial w}\\sum_{i=0}^{N-1} \\big(y_i - w x_i - c \\big)^2 \\\\\n",
    "    &= \\sum_{i=0}^{N-1} \\frac{\\partial}{\\partial w} \\big(y_i - w x_i - c \\big)^2 \\\\\n",
    "    &= \\sum_{i=0}^{N-1} 2 \\big(y_i - w x_i - c \\big) \\frac{\\partial}{\\partial w} \\big(y_i - w x_i - c \\big) \\\\\n",
    "    &= \\sum_{i=0}^{N-1} 2 \\big(y_i - w x_i - c \\big) \\big(- x_i \\big) \\\\\n",
    "    &= \\sum_{i=0}^{N-1} 2 x_i (w x_i + c - y_i) \\\\\n",
    "\\end{align}\n",
    "\n",
    "For $c$ we have:\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial E}{\\partial c} \n",
    "    &= \\frac{\\partial}{\\partial c}\\sum_{i=0}^{N-1} \\big(y_i - w x_i - c \\big)^2 \\\\\n",
    "    &= \\sum_{i=0}^{N-1} \\frac{\\partial}{\\partial c} \\big(y_i - w x_i - c \\big)^2 \\\\\n",
    "    &= \\sum_{i=0}^{N-1} 2 \\big(y_i - w x_i - c \\big) \\frac{\\partial}{\\partial c} \\big(y_i - w x_i - c \\big) \\\\\n",
    "    &= \\sum_{i=0}^{N-1} 2 \\big(y_i - w x_i - c \\big) \\big(- 1 \\big) \\\\\n",
    "    &= \\sum_{i=0}^{N-1} 2 (w x_i + c - y_i) \\\\\n",
    "\\end{align}\n",
    "\n",
    "Now setting $\\frac{\\partial E}{\\partial w} = 0$:\n",
    "\n",
    "\\begin{align}\n",
    "\\Rightarrow \\sum_{i=0}^{N-1} x_i (w x_i + c - y_i) &= 0 \\\\\n",
    "w \\sum_{i=0}^{N-1} (x_i)^2 + c \\sum_{i=0}^{N-1} (x_i) - \\sum_{i=0}^{N-1} (x_i y_i) &= 0 \\\\\n",
    "\\Rightarrow A w + B c - C &= 0 \\qquad\\text{(1)}\n",
    "\\end{align}\n",
    "\n",
    "with $A = \\sum_{i=0}^{N-1} (x_i)^2$, $B = \\sum_{i=0}^{N-1} (x_i)$ and $C = \\sum_{i=0}^{N-1} (x_i y_i)$.\n",
    "\n",
    "Now setting $\\frac{\\partial E}{\\partial c} = 0$:\n",
    "\n",
    "\\begin{align}\n",
    "\\Rightarrow \\sum_{i=0}^{N-1} (w x_i + c - y_i) &= 0 \\\\\n",
    "w \\sum_{i=0}^{N-1} (x_i) + N c - \\sum_{i=0}^{N-1} (y_i) &= 0 \\\\\n",
    "\\Rightarrow B w + N c - D &= 0 \\qquad\\text{(2)}\n",
    "\\end{align}\n",
    "\n",
    "with $D = \\sum_{i=0}^{N-1} (y_i)$.\n",
    "\n",
    "Equations (1) and (2) can be combined to solve for the optimal values for $w$ and $c$! <p style=\"color:red\">**You will need to finish the final part of the maths to write the function required in question 3.**</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3 (3 marks)\n",
    "\n",
    "Write a function using numpy that takes as input `x_train` and `y_train` and outputs the optimal $w$ and $c$ parameters for least squares linear regression.\n",
    "\n",
    "Confirm that your results are reasonable by plotting the resulting linear function on the training graph as well as the predicted values for the test set on a separate graph.\n",
    "\n",
    "- 2 marks for the correct solution for $w$ and $c$\n",
    "- 1 mark for the plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def least_squares_analytic_solution(x, y):\n",
    "    # Add code to calcuate the optimal w and c using the \n",
    "    # equations above..\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    return w, c\n",
    "\n",
    "\n",
    "try:\n",
    "    w_opt, c_opt = least_squares_analytic_solution(x_train, y_train)\n",
    "    \n",
    "    print('Analytic solution:')\n",
    "    print('Analytic w = ', w_opt)\n",
    "    print('Analytic c = ', c_opt)\n",
    "    \n",
    "except Exception as err:\n",
    "    print('Error during least squares analytic solution:', err)\n",
    "\n",
    "    \n",
    "def plot_estimated_y_for_input_x(w, c):\n",
    "    # Add code to plot a line showing your solution \n",
    "    # for w and c..\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    pass\n",
    "    \n",
    "\n",
    "try:\n",
    "    plot_data(x_train, y_train)\n",
    "    plt.title('Analytic Linear Regression (Training Data)')\n",
    "    plot_estimated_y_for_input_x(w_opt, c_opt)\n",
    "    \n",
    "    print('Mean least squares error on TRAINING data = ',\n",
    "          least_squares_error(x_train, y_train, w_opt, c_opt) / x_train.shape[0])\n",
    "\n",
    "    plot_data(x_test, y_test)\n",
    "    plt.title('Analytic Linear Regression (Testing Data)')\n",
    "    plot_estimated_y_for_input_x(w_opt, c_opt)\n",
    "    \n",
    "    print('Mean least squares error on TEST data = ',\n",
    "          least_squares_error(x_test, y_test, w_opt, c_opt) / x_test.shape[0])\n",
    "    \n",
    "except Exception as err:\n",
    "    print('Error during plotting:', err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tensorflow Optimisation\n",
    "\n",
    "**For help with this please refer to the separate tensorflow example that was worked through in class and is available on moodle as a separate jupyter notebook. Please load the example notebook and run through it yourself again before continuing with this task.**\n",
    "\n",
    "In the previous question we looked at an analytic solution to the least square problem. We now pretend that we could not solve the problem analytically. Although this is not true for this case, the additional of extensions to the linear regression model (for example to improve robustness, adding feature selection or handling non-linear data) can mean that it is no longer possible to find an analytic solution and numerical optimisation must be used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4 (2 marks)\n",
    "\n",
    "Add the tensorflow expressions to the following code to calculate the least squares error using tensorflow and check that it calculates the same value as the numpy version.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.reset_default_graph()\n",
    "\n",
    "try:\n",
    "    # Constants to hold the training data..\n",
    "    t_x_train = tf.constant(x_train, name='x_train')\n",
    "    t_y_train = tf.constant(y_train, name='y_train')\n",
    "except Exception as err:\n",
    "    print('Error defining training data:', err)\n",
    "    \n",
    "# Initial values for optimisation..\n",
    "w_initial_guess = 1.5\n",
    "c_initial_guess = 0.5\n",
    "\n",
    "# Variables to hold w and c\n",
    "t_w = tf.Variable(w_initial_guess, \n",
    "                  dtype=tf.float64, \n",
    "                  name='w')\n",
    "t_c = tf.Variable(c_initial_guess, \n",
    "                  dtype=tf.float64, \n",
    "                  name='c')\n",
    "\n",
    "def calculate_tf_least_squares_error(t_x_train, t_y_train, t_w, t_c):\n",
    "    # Add your code here to calculate t_least_squares_error..\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    return t_least_squares_error\n",
    "\n",
    "try:\n",
    "    t_least_squares_error = calculate_tf_least_squares_error(t_x_train, t_y_train, t_w, t_c)\n",
    "except Exception as err:\n",
    "    print('Error declaration of calculate_tf_least_squares_error():', err)\n",
    "\n",
    "try:\n",
    "    with tf.Session() as session:\n",
    "        session.run(tf.global_variables_initializer())\n",
    "\n",
    "        tensorflow_result = session.run(t_least_squares_error)\n",
    "\n",
    "    # Check with the result from your previous function..\n",
    "    numpy_result = least_squares_error(x_train, y_train, \n",
    "                                       w=w_initial_guess, \n",
    "                                       c=c_initial_guess)\n",
    "\n",
    "    print('Tensorflow least squares error = ', tensorflow_result)\n",
    "    print('numpy least squares error = ', numpy_result)\n",
    "\n",
    "    # This should pass if they are the same to nummerical precision!\n",
    "    assert(np.isclose(tensorflow_result, numpy_result))\n",
    "except Exception as err:\n",
    "    print('Error during optimisation with calculate_tf_least_squares_error():', err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5 (2 marks)\n",
    "\n",
    "Check that the gradients from tensorflow are correct by writing numpy code to calculate the value of the derivatives from the analytic expressions (derived previously):\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial E}{\\partial w} \n",
    "    &= \\sum_{i=0}^{N-1} 2 x_i (w x_i + c - y_i) \\\\\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial E}{\\partial c} \n",
    "    &= \\sum_{i=0}^{N-1} 2 (w x_i + c - y_i) \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.reset_default_graph()\n",
    "\n",
    "try:\n",
    "    # Constants to hold the training data..\n",
    "    t_x_train = tf.constant(x_train, name='x_train')\n",
    "    t_y_train = tf.constant(y_train, name='y_train')\n",
    "except Exception as err:\n",
    "    print('Error defining training data:', err)\n",
    "    \n",
    "# Initial values for optimisation..\n",
    "w_initial_guess = 1.5\n",
    "c_initial_guess = 0.5\n",
    "\n",
    "# Variables to hold w and c\n",
    "t_w = tf.Variable(w_initial_guess, \n",
    "                  dtype=tf.float64, \n",
    "                  name='w')\n",
    "t_c = tf.Variable(c_initial_guess, \n",
    "                  dtype=tf.float64, \n",
    "                  name='c')\n",
    "\n",
    "try:\n",
    "    t_least_squares_error = calculate_tf_least_squares_error(t_x_train, t_y_train, t_w, t_c)\n",
    "    \n",
    "    with tf.Session() as session:\n",
    "        session.run(tf.global_variables_initializer())\n",
    "\n",
    "        t_gradient_wrt_w = tf.gradients(t_least_squares_error, t_w)\n",
    "        t_gradient_wrt_c = tf.gradients(t_least_squares_error, t_c)\n",
    "\n",
    "        tf_grad_w = session.run(t_gradient_wrt_w)\n",
    "        tf_grad_c = session.run(t_gradient_wrt_c)\n",
    "\n",
    "        print('Tensorflow gradient wrt w = ', tf_grad_w)\n",
    "        print('Tensorflow gradient wrt c = ', tf_grad_c)\n",
    "except Exception as err:\n",
    "    print('Error using calculate_tf_least_squares_error() to find gradients:', err)\n",
    "    \n",
    "def calc_gradients_for_least_squares(x, y, w, c):\n",
    "    # Add your code to evalute the partial derivatives here\n",
    "    # ************************************************************\n",
    "    # ...\n",
    "    \n",
    "    return grad_w, grad_c\n",
    "\n",
    "try:\n",
    "    numpy_grad_w, numpy_grad_c = calc_gradients_for_least_squares(x_train, \n",
    "                                                                  y_train, \n",
    "                                                                  w_initial_guess, \n",
    "                                                                  c_initial_guess)\n",
    "\n",
    "    print('Analytic gradient wrt w = ', numpy_grad_w)\n",
    "    print('Analytic gradient wrt c = ', numpy_grad_c)\n",
    "\n",
    "    # This should pass if they are the same to nummerical precision!\n",
    "    assert(np.isclose(tf_grad_w, numpy_grad_w))\n",
    "    assert(np.isclose(tf_grad_c, numpy_grad_c))\n",
    "except Exception as err:\n",
    "    print('Error during calculation with calc_gradients_for_least_squares():', err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient Descent in Tensorflow\n",
    "\n",
    "If all has gone well, the gradients should be the same and you can use tensorflow to find the solution to the optimisation.\n",
    "\n",
    "**Run the following block of code to see the tensorflow optimisation running with your least square error function from above!**\n",
    "\n",
    "*Aside: Investigate what happens as you change the `learning_rate` parameter as well as the `num_iterations`. Can we guarantee that tensorflow will always return the same result as the analytic solution? What might be happening if not?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.reset_default_graph()\n",
    "\n",
    "try:\n",
    "    # Constants to hold the training data..\n",
    "    t_x_train = tf.constant(x_train, name='x_train')\n",
    "    t_y_train = tf.constant(y_train, name='y_train')\n",
    "except Exception as err:\n",
    "    print('Error defining training data:', err)\n",
    "    \n",
    "# Initial values for optimisation..\n",
    "w_initial_guess = 1.5\n",
    "c_initial_guess = 0.5\n",
    "\n",
    "# Variables to hold w and c\n",
    "t_w = tf.Variable(w_initial_guess, \n",
    "                  dtype=tf.float64, \n",
    "                  name='w')\n",
    "t_c = tf.Variable(c_initial_guess, \n",
    "                  dtype=tf.float64, \n",
    "                  name='c')\n",
    "\n",
    "try:\n",
    "    t_least_squares_error = calculate_tf_least_squares_error(t_x_train, t_y_train, t_w, t_c)\n",
    "\n",
    "    # Create a gradient descent optimiser that uses a\n",
    "    # certain step size (learning_rate)..\n",
    "    optimiser = tf.train.GradientDescentOptimizer(learning_rate=0.001)\n",
    "\n",
    "    # We want to minimise the least squares error\n",
    "    t_objective = t_least_squares_error\n",
    "\n",
    "    # We want to optimise wrt w and c\n",
    "    vars_to_optimise = [t_w, t_c]\n",
    "\n",
    "    minimize_operation = optimiser.minimize(t_objective,\n",
    "                                            var_list=vars_to_optimise)\n",
    "\n",
    "    # Number of iterations to perform\n",
    "    num_iterations = 15\n",
    "\n",
    "    with tf.Session() as session:\n",
    "        session.run(tf.global_variables_initializer())\n",
    "\n",
    "        # Run a number of iterations of gradient descent..\n",
    "        for iteration in range(num_iterations):\n",
    "            # At each iteration evaluate the minimize_operation\n",
    "            # to perform the gradient descent step and also\n",
    "            # keep track of the current value..\n",
    "            step, cost = session.run([minimize_operation, \n",
    "                                      t_least_squares_error])\n",
    "\n",
    "            # Print out the value of log P every 10 iterations..\n",
    "            #if ((iteration + 1) % 10 == 0):\n",
    "            print('iter %4d, E(w,c) = %0.3f' % \n",
    "                      (iteration + 1, cost))\n",
    "\n",
    "        # Get the final results of the optimisation..\n",
    "        w_tf_opt = session.run(t_w)\n",
    "        c_tf_opt = session.run(t_c)\n",
    "\n",
    "        print('\\nAfter tensorflow optimisation:')\n",
    "        print('Tensorflow w = ', w_tf_opt)\n",
    "        print('Tensorflow c = ', c_tf_opt)\n",
    "\n",
    "    print('\\nAnalytic solution:')\n",
    "    print('Analytic w = ', w_opt)\n",
    "    print('Analytic c = ', c_opt)\n",
    "    \n",
    "except Exception as err:\n",
    "    print('Error using calculate_tf_least_squares_error():', err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 6 (5 marks)\n",
    "\n",
    "You are now going to implement your own version of gradient descent (the process that tensorflow is providing to perform optimisation). The good news is that you already have the numpy functions to calculate both the least squares error (`least_squares_error`) and the gradients wrt the parameters (`calc_gradients_for_least_squares`).\n",
    "\n",
    "Starting from the same initial values as tensorflow (`w_initial_guess`, `c_initial_guess`) write an iterative algorithm for gradient descent. At each iteration it should perform the following steps:\n",
    "\n",
    "- Evaluate and save the squared error for the current parameters\n",
    "- Calculate the gradients wrt w and c for the current parameters\n",
    "- Update the parameters for w and c my moving in the direction of the negative current proportional to the current step size\n",
    "\n",
    "**The code below should run with these three additions. Once they are working, move on to the following.**\n",
    "\n",
    "Under this algorithm, the squared error should decrease at each iteration. If it is not decreasing then either there is a bug in the code (!) or the step size is too large. Add a check that makes sure the squared error always decreases and if it fails to decrease, decrease the step size and run the iteration again.\n",
    "\n",
    "*Hint: you don't want to __accept__ a bad step - if a step is bad, reduce the step size and calculate the step again.*\n",
    "\n",
    "**Try running this code starting with `current_step_size = 0.1`.**\n",
    "\n",
    "You can improve your answer by checking to see if you should stop iterating. If the change in the squared error between successive iterations is very small then one of the following is true. Either the step size is too small (not changing the parameters sufficiently) or the values have converged to their optimal values. If you cannot find a step size that creates a deacrease in the squared error then you have probably converged and can stop performing iterations.\n",
    "\n",
    "*Hint: you can exit a `for` loop early with the `break` command.*\n",
    "\n",
    "**Try running this code starting with `num_iterations = 200`.**\n",
    "\n",
    "- 1 mark evaluating and tracking the squared error and gradients\n",
    "- 2 marks taking steps in the negative gradient direction proportional to the step size\n",
    "- 1 mark check for decreasing error and adjusting step size if necessary\n",
    "- 1 mark check for convergence and stop interating"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Keep track of parameter values over iterations..\n",
    "w_current = w_initial_guess\n",
    "c_current = c_initial_guess\n",
    "\n",
    "try:\n",
    "    # Keep track of the error..\n",
    "    E_current = least_squares_error(x_train, y_train, w_current, c_current)\n",
    "except Exception as err:\n",
    "    print('Error defining training data:', err)\n",
    "\n",
    "# Keep track of the step size..\n",
    "current_step_size = 0.001\n",
    "\n",
    "num_iterations = 20\n",
    "\n",
    "converge_threshold = 1e-8\n",
    "\n",
    "try:\n",
    "    for iteration in range(num_iterations):\n",
    "\n",
    "        def run_iteration(x_train, y_train, w_current, c_current, E_current, \n",
    "                          current_step_size, converge_threshold):\n",
    "\n",
    "            # Set to True when converged..\n",
    "            converged = False\n",
    "\n",
    "            # Add code to evaluate the gradients..\n",
    "            # ************************************************************\n",
    "            # ...\n",
    "\n",
    "\n",
    "\n",
    "            # Add code to take a step in the direction of the negative\n",
    "            # gradient proportional to the step size..\n",
    "            # ************************************************************\n",
    "            # ...\n",
    "            # w_new = \n",
    "            # c_new = \n",
    "\n",
    "\n",
    "            # Add code to evaluate and remember the squared error..\n",
    "            # ************************************************************\n",
    "            # ...\n",
    "            # E_new = \n",
    "\n",
    "\n",
    "            # Add code to check that error is decreasing and reduce step\n",
    "            # size if not..\n",
    "            # ************************************************************\n",
    "            # ...\n",
    "\n",
    "\n",
    "\n",
    "            # Add code to check for convergence and terminate\n",
    "            # the loop if converged..\n",
    "            # ************************************************************\n",
    "            # ...\n",
    "\n",
    "\n",
    "\n",
    "            # Take the step\n",
    "            w_current = w_new\n",
    "            c_current = c_new\n",
    "            E_current = E_new\n",
    "\n",
    "            return w_current, c_current, E_current, current_step_size, converged\n",
    "\n",
    "\n",
    "        w_current, c_current, E_current, current_step_size, converged = \\\n",
    "            run_iteration(x_train, y_train, w_current, c_current, E_current, \n",
    "                          current_step_size, converge_threshold)\n",
    "        \n",
    "        print('iteration %4d, E = %f, w = %f, c = %f' % \n",
    "              (iteration, E_current, w_current, c_current))\n",
    "        \n",
    "        if converged:\n",
    "            # Break out of iteration loop..\n",
    "            print('Converged!')\n",
    "            break\n",
    "        \n",
    "    print('\\nAfter gradient descent optimisation:')\n",
    "    print('Optimised w = ', w_current)\n",
    "    print('Optimised c = ', c_current)\n",
    "\n",
    "    print('\\nAnalytic solution:')\n",
    "    print('Analytic w = ', w_opt)\n",
    "    print('Analytic c = ', c_opt)\n",
    "    \n",
    "except Exception as err:\n",
    "        print('Error during run_iteration():', err)"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
