Exercises/1.2.Intro-LinearModels-SGD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/makagan/TRISEP_Tutorial/blob/main/Exercises/1.2.Intro-LinearModels-SGD.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Regression\n",
    "\n",
    "Adapted from [notebook](https://github.com/drinkingkazu/TRISEP2022ML/blob/master/LinearModels/Introduction00-LinearModels.ipynb) by Kazuhiro Terao\n",
    "\n",
    "## Goals:\n",
    "* Practice linear regression and logistic regression\n",
    "* Implement a (stochastic) gradient descent\n",
    "* Play with Python :)\n",
    "\n",
    "Let's start with basic imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "mpl.rcParams['figure.figsize'] = [8, 6]\n",
    "mpl.rcParams['font.size'] = 16\n",
    "mpl.rcParams['axes.grid'] = True\n",
    "import numpy as np "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear Regression\n",
    "\n",
    "For this practice, we will generate a fake data $y=ax+b+\\epsilon$ where $\\epsilon$ is a gaussian noise $\\mathcal{N}(\\mu=0,\\sigma=1)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sample_linear_regression(a,b,num_sample=1000,seed=123):\n",
    "    '''\n",
    "    Generate a fake data sample for linear regression of y=ax+b\n",
    "    y values are smeared by a normal distribution\n",
    "    INPUT:\n",
    "      - a ... float, slope part of a linear equation\n",
    "      - b ... float, offset part of a linear equation\n",
    "    '''    \n",
    "    np.random.seed(seed)\n",
    "    xval = np.random.random(num_sample)*10\n",
    "    yval = a*xval + b + np.random.normal(scale=1.0,size=num_sample)\n",
    "\n",
    "    return xval,yval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a sample with $(a,b)=(2.0,1.5)$ with 40 points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a,b=2.0,1.5\n",
    "x,y=sample_linear_regression(a,b,num_sample=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize using `matplotlib`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "fig,ax=plt.subplots(figsize=(12,8),facecolor='w')\n",
    "plt.plot(x,y,linestyle='',marker='o',markersize=10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exact solution for a univariate regression\n",
    "\n",
    "Our model is $Q(\\mathbf{x},\\mathbf{w})$ is $w_0+w_1x_i$. We would like to minimize the MSE:\n",
    "$$\n",
    "\\text{Min. } \\displaystyle{\\left(\\sum_i (y_i - Q(x_i,\\mathbf{w}))^2\\right)} \\; \\text{  where  } \\; Q(x_i,\\mathbf{w}) = w_0 + w_1x_i\n",
    "$$\n",
    "\n",
    "The exact solution for is available as:\n",
    "\n",
    "$$\n",
    "w_1 = \\displaystyle{\\frac{\\sum(x_i-\\langle x\\rangle)(y_i-\\langle y\\rangle)}{\\sum(x_i-\\langle x\\rangle)^2}} \\;\\text{ and } \\; w_0 = \\langle y\\rangle - w_1\\cdot\\langle x\\rangle\n",
    "$$\n",
    "\n",
    "### Exercise 1: implement the exact solution to solve for \n",
    "Let's implement and see how it works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "w0,w1=0,0\n",
    "### Your code below ###\n",
    "\n",
    "#######################\n",
    "print('Analytical solution for w0',w0,'and w1',w1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent\n",
    "\n",
    "Another approach to the solution is to use a _Gradient Descent_ technique as we covered in the lecture. In this case, we compute the error (i.e. MSE) and update the weights according to the equation:\n",
    "$$\n",
    "w_{\\text{new}} = w - \\lambda\\nabla_{\\mathbf{w}} L\n",
    "$$\n",
    "First of all, let's use a bias trick. Our data = input feature vector $\\mathbf{x}$ has only 1 feature. If we expand to $\\mathbf{X}$ and add an identity as the second feature, we can write\n",
    "$$\n",
    "Q(\\mathbf{x},\\mathbf{w}) = \\mathbf{X}\\cdot\\mathbf{w}\n",
    "$$\n",
    "A gradient for a linear regression is also covered in the lecture:\n",
    "$$\n",
    "\\nabla_{\\mathbf{w}}L = -2\\left(\\mathbf{y} - Q(\\mathbf{X},\\mathbf{w})\\right)\\cdot\\mathbf{X}\n",
    "$$\n",
    "\n",
    "Let's try to implement on our own!\n",
    "\n",
    "### Model base class\n",
    "In this notebook we will define two linear models. To save our time (sorry, actually, for your practice), let's define a common base class so that most functionalities as a machine learning model can be shared. Here's the base class definition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class linear_model_base:\n",
    "    '''\n",
    "    The base class for linear models in this notebook.\n",
    "    This class (and its inherited children) implements 4 core functions.\n",
    "      1. forward computes the model prediction (and loss, if label is provided).\n",
    "      2. compute_loss calculates the loss and gradient given the data and label using the current parameters (weights).\n",
    "      3. backward update the gradient which is computed within forward (must call forward beforehand).\n",
    "    '''\n",
    "    def __init__(self,num_features,lr=0.001):\n",
    "        '''\n",
    "        ARGUMENTS:\n",
    "        num_features ... is the number of regression targets. y=w0 + w1*x would be 2.\n",
    "          lr ... is the learning rate fir a gradient update (optional)\n",
    "        '''\n",
    "        self._num_features = int(num_features)\n",
    "        self._w = np.zeros(shape=(num_features),dtype=np.float32)\n",
    "        self._lr = float(lr)\n",
    "        self._grad = None\n",
    "                \n",
    "    def forward(self,x,y=None):\n",
    "        '''\n",
    "        This function computes and returns a prediction. If an optional label information\n",
    "        is given, also returns the loss value and calculates the gradient internally.\n",
    "        ARGUMENTS:\n",
    "          x ... is input data of the shape [N,2] for N samples of instances [x_i,1] \n",
    "          y ... is label (optional)\n",
    "        RETURNS:\n",
    "          prediction ... model predicted values\n",
    "          loss ... the current loss value for the input data (only if y input is provided)\n",
    "        '''\n",
    "        prediction = np.dot(x,self._w)\n",
    "        if y is None: \n",
    "            return prediction\n",
    "        loss, self._grad = self.compute_loss(x,y)\n",
    "        return prediction,loss\n",
    "\n",
    "    def compute_loss(self,x,y):\n",
    "        '''\n",
    "        This is a place holder function for computing the loss. The base class implements 0 gradient\n",
    "        and invalid loss (-1).\n",
    "        ARGUMENTS:\n",
    "          x ... is input data of the shape [N,2] for N samples of instances [x_i,1]\n",
    "          y ... is label\n",
    "        RETURNS:\n",
    "          loss, gradient\n",
    "        '''\n",
    "        return -1, 0.\n",
    "    \n",
    "    def backward(self):\n",
    "        '''\n",
    "        This function applies pre-computed gradients to update the model parameters (weights).\n",
    "        '''\n",
    "        if self._grad is None:\n",
    "            raise SyntaxError('Must call forward before backward!')\n",
    "        self._w = self._w - self._lr * (self._grad)\n",
    "        self._grad=None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we assumed the input data $x$ has a shape $(N,2)$ for $N$ data instances of $[x_i,\\mathbb{1}]$.\n",
    "\n",
    "It includes an extra identity constant in the features so that we can treat the constnat (i.e. \"bias\", $w_0$) term in a unified matrix formula.\n",
    "\n",
    "### Exercise 2: implement a linear regression class\n",
    "This class should inherit from `linear_model_base` and implement model-specific (MSE loss + linear regression) function definitions. \n",
    "\n",
    "Give a shot looking at equations we covered earlier!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class linear_regression(linear_model_base):\n",
    "    \n",
    "    def compute_loss(self,x,y):\n",
    "        '''\n",
    "        Compute the MSE loss\n",
    "        ARGUMENTS:\n",
    "          x ... is input data of the shape [N,2] for N samples of instances [x_i,1]\n",
    "          y ... is label\n",
    "        RETURNS:\n",
    "          loss, gradient\n",
    "        '''\n",
    "        # your code here\n",
    "        \n",
    "        return -1.,0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train loop\n",
    "Let's write a code to run the loop of optimization (training loop).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_gradient_descent(x,y,model,num_iterations=100):\n",
    "    '''\n",
    "    Run gradient descent for a provided model\n",
    "    ARGUMENTS:\n",
    "      x ... is input data of the shape [N,2] for N samples of instances [x_i,1]\n",
    "      y ... input label\n",
    "      model ... model to be optimized\n",
    "      num_interations ... number of times parameters are updated\n",
    "    RETURNS:\n",
    "      loss ... an array of loss at each iteration\n",
    "      ws ... an array of the model weights at each iteration\n",
    "    '''\n",
    "    # Make sure the data shape matches the model's parameter count and label count\n",
    "    assert(x.shape[1]==model._w.shape[0])\n",
    "    assert(x.shape[0]==y.shape[0])\n",
    "    \n",
    "    loss, ws = [],[]\n",
    "    for idx in range(num_iterations):\n",
    "            \n",
    "        loss.append(model.forward(x,y)[1])\n",
    "        ws.append(model._w)\n",
    "        model.backward()\n",
    "        \n",
    "    return np.array(loss), np.array(ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's run the training! \n",
    "\n",
    "Just one last point: don't forget to add an identity column in the input data $x$ as it was assumed for the model input. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the input\n",
    "data_input = np.column_stack([np.ones(shape=(len(x),1),dtype=np.float32),x])\n",
    "data_label = y\n",
    "# Construct the model\n",
    "model = linear_regression(2)\n",
    "\n",
    "# Run the training\n",
    "loss,ws = train_gradient_descent(data_input,data_label,model,100)\n",
    "\n",
    "# Plot the loss\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig,ax = plt.subplots(figsize=(12,8),facecolor='w')\n",
    "plt.plot(loss,marker='o',markersize=10)\n",
    "plt.yscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import imageio\n",
    "import io\n",
    "import base64\n",
    "from IPython.display import HTML\n",
    "\n",
    "def make_plot(x,y,ws):\n",
    "    \n",
    "        fig,ax = plt.subplots(figsize=(12,8),facecolor='w')\n",
    "        ax.plot(x,y,markersize=10,marker='o',linestyle='')\n",
    "        xmin=np.min(x)\n",
    "        xmax=np.max(x)\n",
    "        ymin = ws[1]*xmin+ws[0]\n",
    "        ymax = ws[1]*xmax+ws[0]\n",
    "        ax.plot((xmin,xmax),(ymin,ymax),linewidth=2,color='red')\n",
    "        ax.set_xlim((xmin,xmax))\n",
    "        ax.set_ylim((0,np.max(y)))\n",
    "        ax.set_ylabel('$y$ (response)',fontsize=18)\n",
    "        ax.set_xlabel('$x$ (feature)',fontsize=18)\n",
    "        ax.tick_params(labelsize=15)\n",
    "        ax.grid()\n",
    "\n",
    "        fig.canvas.draw()       # draw the canvas, cache the renderer\n",
    "        image = np.frombuffer(fig.canvas.tostring_rgb(), dtype='uint8')\n",
    "        image  = image.reshape(fig.canvas.get_width_height()[::-1] + (3,))\n",
    "        plt.close()\n",
    "        return image\n",
    "\n",
    "kwargs_write = {'fps':1.0, 'quantizer':'nq'}\n",
    "imageio.mimsave('./anime.gif', [make_plot(x,y,ws[i]) for i in range(len(ws))], fps=20)\n",
    "\n",
    "video = io.open('anime.gif', 'r+b').read()\n",
    "encoded = base64.b64encode(video)\n",
    "HTML(data='''<img src=\"data:image/gif;base64,{0}\" type=\"gif\" />'''.format(encoded.decode('ascii')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loss surface\n",
    "\n",
    "Let's try to visualize how the loss distribution looks over some range of weights look like + our gradient descent's stepping points\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the loss on the grid of w0 and w1 values\n",
    "w0s = np.arange(-10,10,0.1)\n",
    "w1s = np.arange(0,5,0.1)\n",
    "grid = np.array(np.meshgrid(w0s,w1s,sparse=False)).reshape(2,-1)\n",
    "loss_map = np.zeros(shape=(grid.shape[1]),dtype=np.float32)\n",
    "for i in range(grid.shape[1]):\n",
    "    model._w = grid[0,i],grid[1,i]\n",
    "    loss_map[i] = model.compute_loss(data_input,data_label)[0]\n",
    "loss_map = loss_map.reshape(len(w1s),len(w0s))\n",
    "\n",
    "# Plot the map and final fit\n",
    "from matplotlib.colors import LogNorm\n",
    "fig,axes=plt.subplots(1,2,figsize=(18,6),facecolor='w')\n",
    "# Plot the map\n",
    "axes[0].contourf(w0s,w1s,loss_map,norm=LogNorm())\n",
    "axes[0].plot(ws[:,0],ws[:,1],marker='o',linestyle='', markersize=4, color='magenta')\n",
    "axes[0].set_xlabel('$w_0$',fontsize=18)\n",
    "axes[0].set_ylabel('$w_1$',fontsize=18)\n",
    "# Plot the fit\n",
    "axes[1].plot(x,y,linestyle='',marker='o',markersize=10)\n",
    "xmin,xmax=np.min(x),np.max(x)\n",
    "ymin = ws[-1][1]*xmin+ws[-1][0]\n",
    "ymax = ws[-1][1]*xmax+ws[-1][0]\n",
    "axes[1].plot((xmin,xmax),(ymin,ymax),linewidth=2,color='red')\n",
    "axes[1].set_xlabel('x (feature)',fontsize=18)\n",
    "axes[1].set_ylabel('y (response)',fontsize=18)\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "    ax.tick_params(labelsize=15)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stochastic GD (Exercise 3)\n",
    "Next, let's try to implement a stochastic GD. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def train_sgd(x,y,model,batch_size,num_iterations=100):\n",
    "    '''\n",
    "    Run mini-batch stochastic gradient descent for a provided model\n",
    "    ARGUMENTS:\n",
    "      x ... is input data of the shape [N,2] for N samples of instances [x_i,1]\n",
    "      y ... input label\n",
    "      model ... model to be optimized\n",
    "      batch_size ... the batch size to be used for SGD.\n",
    "      num_interations ... number of times parameters are updated\n",
    "    RETURNS:\n",
    "      loss ... an array of loss at each iteration\n",
    "      ws ... an array of the model weights at each iteration\n",
    "    '''\n",
    "    assert(x.shape[1]==model._w.shape[0])\n",
    "    assert(x.shape[0]==y.shape[0])\n",
    "    assert(batch_size>0)\n",
    "    loss, ws = [],[]\n",
    "    \n",
    "    # Your code here\n",
    "        \n",
    "    return np.array(loss), np.array(ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training SGD\n",
    "Let's run it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the input\n",
    "data_input = np.column_stack([np.ones(shape=(len(x),1),dtype=np.float32),x])\n",
    "data_label = y\n",
    "# Construct the model\n",
    "model = linear_regression(2)\n",
    "\n",
    "# Run the training\n",
    "loss,ws = train_sgd(x=data_input,\n",
    "                    y=data_label,\n",
    "                    model=model,\n",
    "                    batch_size=5,\n",
    "                    num_iterations=100)\n",
    "\n",
    "# Plot the loss\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig,ax = plt.subplots(figsize=(12,8),facecolor='w')\n",
    "plt.plot(loss,marker='o',markersize=10)\n",
    "plt.yscale('log')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the loss on the grid of w0 and w1 values\n",
    "w0s = np.arange(-10,10,0.1)\n",
    "w1s = np.arange(0,5,0.1)\n",
    "grid = np.array(np.meshgrid(w0s,w1s,sparse=False)).reshape(2,-1)\n",
    "loss_map = np.zeros(shape=(grid.shape[1]),dtype=np.float32)\n",
    "for i in range(grid.shape[1]):\n",
    "    model._w = grid[0,i],grid[1,i]\n",
    "    loss_map[i] = model.compute_loss(data_input,data_label)[0]\n",
    "loss_map = loss_map.reshape(len(w1s),len(w0s))\n",
    "\n",
    "# Plot the map and final fit\n",
    "from matplotlib.colors import LogNorm\n",
    "fig,axes=plt.subplots(1,2,figsize=(18,6),facecolor='w')\n",
    "# Plot the map\n",
    "axes[0].contourf(w0s,w1s,loss_map,norm=LogNorm())\n",
    "axes[0].plot(ws[:,0],ws[:,1],marker='o',linestyle='', markersize=4, color='magenta')\n",
    "axes[0].set_xlabel('$w_0$',fontsize=18)\n",
    "axes[0].set_ylabel('$w_1$',fontsize=18)\n",
    "# Plot the fit\n",
    "axes[1].plot(x,y,linestyle='',marker='o',markersize=10)\n",
    "xmin,xmax=np.min(x),np.max(x)\n",
    "ymin = ws[-1][1]*xmin+ws[-1][0]\n",
    "ymax = ws[-1][1]*xmax+ws[-1][0]\n",
    "axes[1].plot((xmin,xmax),(ymin,ymax),linewidth=2,color='red')\n",
    "axes[1].set_xlabel('x (feature)',fontsize=18)\n",
    "axes[1].set_ylabel('y (response)',fontsize=18)\n",
    "for ax in axes:\n",
    "    ax.grid()\n",
    "    ax.tick_params(labelsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4\n",
    "\n",
    "Compare how the behavior of the parameter updates (magenta points in $w_0$ v.s. $w_1$ contour plot) between GD and SGD with a batch size 5. How do they behave differently? \n",
    "\n",
    "## Logistic Regression\n",
    "\n",
    "Next, let's try a linear model on a simple classification task: we will generate two class of data points in 2D using a gaussian distribution with different mean values (+1 and -1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123)\n",
    "sample_stat = 2000\n",
    "a = np.random.normal(1,1,sample_stat).reshape(int(sample_stat/2),2)\n",
    "b = np.random.normal(-1,1,sample_stat).reshape(int(sample_stat/2),2)\n",
    "a = np.column_stack([a,np.zeros(shape=(len(a)),dtype=np.float32)])\n",
    "b = np.column_stack([b,np.ones(shape=(len(b)),dtype=np.float32)])\n",
    "data = np.concatenate([a,b])\n",
    "np.random.shuffle(data)\n",
    "x,y = data[:,:2],data[:,2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualize the dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(12,8),facecolor='w')\n",
    "pos = y>0\n",
    "neg = y<1\n",
    "ax.plot(x[pos][:,0],x[pos][:,1],marker='o',markersize=6,linestyle='',color='red')\n",
    "ax.plot(x[neg][:,0],x[neg][:,1],marker='x',markersize=6,linestyle='',color='blue')\n",
    "plt.ylabel('$x_0$ ',fontsize=18)\n",
    "plt.xlabel('$x_1$ ',fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5\n",
    "\n",
    "Inherit from `linear_mode_base` and implement a `logistic_regression` class.\n",
    "\n",
    "$$\n",
    "L = \\frac{1}{m}\\sum(-\\mathbf{y}^T\\cdot\\log Q(\\mathbf{x},\\mathbf{w}) - (1-y)^T\\cdot\\log (1-Q(\\mathbf{x},\\mathbf{w}))\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf{w}}L=\\mathbf{x}^T(Q(\\mathbf{x},\\mathbf{w}) - \\mathbf{y})\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class logistic_regression(linear_model_base):\n",
    "    \n",
    "    def compute_loss(self,x,y):\n",
    "        '''\n",
    "        Compute the MSE loss\n",
    "        ARGUMENTS:\n",
    "          x ... is input data of the shape [N,2] for N samples of instances [x_i,1]\n",
    "          y ... is label\n",
    "        RETURNS:\n",
    "          loss, gradient\n",
    "        '''\n",
    "        # Your code here\n",
    "        return -1, 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training your logistic regression model (Exercise 6)\n",
    "Let's run a similar loop!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the decision boundary line (Exercise 7)\n",
    "\n",
    "Plot the decision boundary (i.e. a linear equation) and overlay on top of the data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sigmoid probability distribution\n",
    "\n",
    "For fun, let's plot the sigmoid probability distribution like we tried for a linear regression's loss surface.\n",
    "\n",
    "As long as it's clear where the decision boundaries are and you are plotting a sigmoid probability, it would be fine (so be artistic)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.colors import LogNorm\n",
    "\n",
    "# Just for fun, this time let's make 2D histogram instead of meshgrid+contour\n",
    "nx,ny=300,300\n",
    "prob_map=np.zeros(shape=(nx*ny,3),dtype=np.float32)\n",
    "vrange=((data_input[:,1].min(),data_input[:,1].max()),(data_input[:,2].min(),data_input[:,2].max()))\n",
    "for i in range(nx):\n",
    "    x0 = (vrange[0][1]-vrange[0][0])/nx * i + vrange[0][0]\n",
    "    for j in range(ny):\n",
    "        x1 = (vrange[1][1]-vrange[1][0])/ny * j + vrange[1][0]\n",
    "        prob = abs(model.sigmoid(model.forward(np.array([1,x0,x1]))) - 0.5)\n",
    "        prob_map[i*ny+j]=(x0,x1,prob)\n",
    "\n",
    "fig,ax=plt.subplots(figsize=(12,8),facecolor='w')\n",
    "# Plot the histogram 2d\n",
    "ax.hist2d(prob_map[:,0],prob_map[:,1],bins=(nx,ny),weights=prob_map[:,2],cmap='inferno')\n",
    "# Plot the data points\n",
    "pos = y>0\n",
    "neg = y<1\n",
    "ax.plot(x[pos][:,0],x[pos][:,1],marker='o',markersize=6,linestyle='',color='red')\n",
    "ax.plot(x[neg][:,0],x[neg][:,1],marker='x',markersize=6,linestyle='',color='blue')\n",
    "plt.ylabel('$x_1$',fontsize=18)\n",
    "plt.xlabel('$x_0$',fontsize=18)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ROC curve (Exercise 8)\n",
    "\n",
    "For this trained model, plot the ROC curve using a different score threshold value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression using Pytorch\n",
    "\n",
    "Hopefully you are very familiar with linear regression now. \n",
    "\n",
    "Let's now learn one more thing, Pytorch machine learning library, which we will use to play with neural networks later.\n",
    "\n",
    "### Defining a model\n",
    "\n",
    "First, let's define a model for a linear regression. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch \n",
    "\n",
    "class linear_regression(torch.nn.Module):\n",
    "    \n",
    "    def __init__(self,num_features):\n",
    "        super(linear_regression,self).__init__()\n",
    "        # Define a linear model with a bias term\n",
    "        self._linear = torch.nn.Linear(num_features,1,bias=True)\n",
    "        \n",
    "    def forward(self,x):\n",
    "        return self._linear(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...yes, that's it! \n",
    "\n",
    "In fact, as we can see, we don't even need to define our model as a class since pytorch already provides a model (`torch.nn.Linear`).\n",
    "\n",
    "But we did it anyway as a practice to define our own model by inheriting from `torch.nn.Module` base class.\n",
    "\n",
    "### Generate $y=ax+\\epsilon$ dataset\n",
    "\n",
    "Now let's make some data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the same data generation function we used before\n",
    "a,b=2.0,1.5\n",
    "x,y=sample_linear_regression(a,b,num_sample=100)\n",
    "\n",
    "# Convert them into a torch input. Shape should be (N,F) where N=number of samples, F=features/labels\n",
    "torch_input = torch.Tensor(x).view(-1,1)\n",
    "torch_label = torch.Tensor(y).view(-1,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's define a loss and write a train loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_torch(model,num_iterations=100):\n",
    "    # Create a MSE loss module\n",
    "    criterion = torch.nn.MSELoss()\n",
    "    # Create an optimizer\n",
    "    optimizer = torch.optim.SGD(model.parameters(), lr = 0.001)\n",
    "\n",
    "    # Now we run the training!\n",
    "    loss_v=[]\n",
    "    for idx in range(num_iterations): \n",
    "\n",
    "        # This is \"forward\" computation\n",
    "        prediction = model(torch_input) \n",
    "\n",
    "        # Compute the loss\n",
    "        loss = criterion(prediction, torch_label) \n",
    "\n",
    "        # Next, we call \"backward\" to update the weights.\n",
    "        # Don't forget to reset the gradient each time we call \"backward\".\n",
    "        optimizer.zero_grad() \n",
    "        loss.backward() \n",
    "\n",
    "        # Now let optimizer to \"take a step\"\n",
    "        optimizer.step() \n",
    "\n",
    "        # Record the loss\n",
    "        loss_v.append(loss.item())\n",
    "        \n",
    "    return loss_v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train a pytorch model\n",
    "\n",
    "Create our model instance + train + plot the loss!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a model \n",
    "model = linear_regression(1) \n",
    "\n",
    "# Train\n",
    "loss = train_torch(model)\n",
    "\n",
    "# Plot the loss\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "fig,ax = plt.subplots(figsize=(12,8),facecolor='w')\n",
    "plt.plot(loss,marker='o',markersize=10)\n",
    "plt.ylabel('Loss',fontsize=18)\n",
    "plt.xlabel('Iterations',fontsize=18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check the parameter values\n",
    "\n",
    "Check the trained parameter values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Weight',model._linear.weight.item())\n",
    "print('Bias',model._linear.bias.item())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}