{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Regression :- Implementation , Hyperparameters and their Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Libraries import"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import tabulate\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "\n",
    "from sklearn import datasets\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "full_data = datasets.load_diabetes()\n",
    "one_feature = full_data.data[: , np.newaxis , 2]\n",
    "\n",
    "X_train = one_feature[:-100] \n",
    "y_train = full_data.target[:-100]\n",
    "\n",
    "X_test = one_feature[-100:]\n",
    "y_test = full_data.target[-100 :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot \n",
    "Plotting Training Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18536 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X_train , y_train)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation\n",
    "In statistics , peoples tends to denote Linear Regression as\n",
    "$$ y = mx + c $$\n",
    "\n",
    "where `c` (bias) can be computed as\n",
    "\n",
    "$$ bias = \\frac{ \\sum_((x_i - mean(i))(y_i - mean(i))} {\\sum_((x_i - mean(x))^2} $$\n",
    "\n",
    "and weights `m` can be computed as \n",
    "\n",
    "$$ weight = mean(y) - (bias * mean(x)) $$\n",
    "\n",
    "### Ordinary Least Square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ordinary_least_square(x , y):\n",
    "    '''\n",
    "    Find weigths for the hypothesis using\n",
    "        y = mx + c\n",
    "    notation.\n",
    "    \n",
    "    Params:\n",
    "        x: Independent data of training set.\n",
    "        y: Dependent data of training set.\n",
    "        \n",
    "    Returns:\n",
    "        Hypothesis function.\n",
    "    '''\n",
    "    \n",
    "    xmean = np.mean(x)\n",
    "    ymean = np.mean(y)\n",
    "    \n",
    "    n = len(x)\n",
    "    \n",
    "    numerator = 0\n",
    "    denominator = 0\n",
    "    \n",
    "    for i in range(n):\n",
    "        numerator += ((x[i] - xmean) * (y[i] - ymean))\n",
    "        denominator += (x[i] - xmean)**2\n",
    "        \n",
    "    b1 = numerator /denominator\n",
    "    b0 = ymean - (b1 * xmean)\n",
    "    \n",
    "    c = b0[0]\n",
    "    m = b1[0]\n",
    "    \n",
    "    return lambda x: m*x + c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a hypothesis\n",
    "ordinary_least_square returns a lambda function which can be treated as maths f(x) function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "hypothesis = ordinary_least_square(X_train , y_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictions = hypothesis(X_train)\n",
    "predictions = predictions.reshape((342 ,))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot data and Regression Line or say Best fit line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22400 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X_train , y_train)\n",
    "plt.plot(X_train , predictions , color='red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CrossCheck Time\n",
    "Now , let's take linearRegression from sklearn to check the correctness of our implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = LinearRegression()\n",
    "model.fit(X_train , y_train)\n",
    "sk_predictions = model.predict(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22392 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(X_train , y_train)\n",
    "plt.plot(X_train , sk_predictions , color = 'red')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluation\n",
    "mean squared error implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_squared_error_custom(truth , predicted):\n",
    "    m = truth.shape[0]\n",
    "    mse = (np.sum((predicted - truth)**2)) / (2 * m)\n",
    "    return mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Our Implementation result 1965.3937776396358\n",
      "Scikit learn Implementation result 1965.3937776396358\n"
     ]
    }
   ],
   "source": [
    "print(\"Our Implementation result {}\".format(mean_squared_error_custom(y_train , predictions)))\n",
    "print(\"Scikit learn Implementation result {}\".format(mean_squared_error_custom(y_train , sk_predictions)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Batch Gradient Descent\n",
    "* X is features array with built-in bias\n",
    "* w initializations for weights\n",
    "* y ground truth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = full_data.data[:, np.newaxis , 1]\n",
    "x1 = np.squeeze(np.asarray(x1[:-100]))\n",
    "\n",
    "x2 = full_data.data[:, np.newaxis , 2]\n",
    "x2 = np.squeeze(np.asarray(x2[:-100]))\n",
    "\n",
    "ones = np.ones(len(x1))\n",
    "\n",
    "X = np.array([ones , x1 , x2]).T\n",
    "y = full_data.target[:-100]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def batch_gradient_descent(X, y , lr , n_iter):\n",
    "    weights = np.zeros(X.shape[1])\n",
    "    \n",
    "    history = []\n",
    "    m = len(y)\n",
    "    \n",
    "    predict = lambda x :np.dot(x , weights)\n",
    "    derivative_w = lambda loss: X.T.dot(loss) / m\n",
    "    \n",
    "    for i in range(n_iter):\n",
    "        hypothesis = predict(X)\n",
    "        loss = hypothesis - y\n",
    "        weights = weights - lr * derivative_w(loss)\n",
    "        \n",
    "        if i % 50 == 0:\n",
    "            history.append(mean_squared_error_custom(X.dot(weights) , y))\n",
    "            \n",
    "    return predict , history"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weight Optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictor , history = batch_gradient_descent(X , y , 0.05 , 50000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot training history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff0996781d0>]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 8760 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(history)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make Predictions and get mean_squared_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1963.4934687043767\n"
     ]
    }
   ],
   "source": [
    "gd_predictions = predictor(X)\n",
    "err = mean_squared_error_custom(gd_predictions , y)\n",
    "\n",
    "print(err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can se that mean square error is kind of similar. This is happened because BGD has hyperparameters which can be tuned and edited, while OLS has no hyperparameters, therefore output is always the same\n",
    "\n",
    "### Hyperparameters \n",
    "Lets observe how hyperparameters affects our algorithm.\n",
    "\n",
    "In lrates we define various values for learning rate hyperparameter, and in niterations various values for number of iterations hyperparamenters.\n",
    "\n",
    "Methodology looks as the following: we will run Batch Gradient Descent with each possible combination of hyperparameters and compare them in multiple ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "lrates = [0.1 , 0.5 , 0.01 , 0.001 , 0.0001]\n",
    "niters = [25000 , 50000 , 150000]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test(X ,y):\n",
    "    \"\"\"\n",
    "    Test function iterates over each hyperparameter and put all the results in record list.\n",
    "    \"\"\"\n",
    "    record = []\n",
    "    \n",
    "    for niter in niters:\n",
    "        for lrate in lrates:\n",
    "            start = time.time()\n",
    "            weights , records = batch_gradient_descent(X, y , lrate , niter)\n",
    "            end = time.time() - start\n",
    "            record.append(dict(lrate=lrate ,niter=niter ,w=weights ,history=records ,time=end ))\n",
    "            \n",
    "    return record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = test(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_records(records):\n",
    "    name_map = {0: 'A', 1: 'B', 2: 'C', 3: 'D', 4: 'E', 5: 'F', 6: 'G', 7: 'H', 8: 'I', \n",
    "                9: 'J', 10: 'K', 11: 'L', 12: 'M', 13: 'N', 14: 'O'}\n",
    "    ax , fig = plt.subplots(5, 3)\n",
    "    ax.tight_layout()\n",
    "    fig = fig.flatten()\n",
    "    \n",
    "    for i , record in enumerate(records):\n",
    "        time = \"{0:.4f}\".format(record['time'])\n",
    "        label = 'lr: {0}, ni: {1}, t: {2}'.format(record[\"lrate\"], record[\"niter\"], time)\n",
    "        fig[i].plot(record['history'] , label=label)\n",
    "        fig[i].set_title(name_map[i])\n",
    "    \n",
    "matplotlib.rcParams[\"figure.figsize\"] = (15 ,15)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 106644 characters]\n",
      "text/plain": [
       "<Figure size 1080x1080 with 15 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_records(rec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot from above is good and clear, lets print the data out as a table.\n",
    "\n",
    "First prepare the data, sort it by time and print."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "time_sorted_records = sorted(rec , key=lambda k : k['time'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "copy_records = time_sorted_records.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "for r in copy_records:\n",
    "    r['Cost']=r['history'][-1]\n",
    "    del r['history']\n",
    "    del r['w']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  learning rate    no of iterations    elapsed time     cost\n",
      "---------------  ------------------  --------------  -------\n",
      "         0.001                25000        0.664344  2846.12\n",
      "         0.01                 25000        0.665427  2299.74\n",
      "         0.0001               25000        0.670085  3014.51\n",
      "         0.5                  25000        0.678696  1963.47\n",
      "         0.1                  25000        0.699473  1963.49\n",
      "         0.5                  50000        1.34833   1963.47\n",
      "         0.01                 50000        1.35192   2079.21\n",
      "         0.1                  50000        1.35432   1963.47\n",
      "         0.001                50000        1.3658    2755.94\n",
      "         0.0001               50000        1.38555   2926.01\n",
      "         0.5                 150000        4.01676   1963.47\n",
      "         0.001               150000        4.02207   2478.86\n",
      "         0.1                 150000        4.05991   1963.47\n",
      "         0.0001              150000        4.07318   2884.86\n",
      "         0.01                150000        4.11481   1965.21\n"
     ]
    }
   ],
   "source": [
    "print(tabulate.tabulate(copy_records , headers={'lrate': 'learning rate', 'niter': 'no of iterations', 'time': 'elapsed time', 'Cost': 'cost'}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data is sorted by time because I choosed this parameter as most important. Generally we want to compute as fast as possible while result is sufficiently good.\n",
    "\n",
    "We can see that generally that cost 3926.93 can be considered as converged.\n",
    "\n",
    "In our case, with a given data, learning rate of .5 converged after 25000 iterations, while learning rate of .0001 has a long way to go until convergence.\n",
    "\n",
    "Lets actually see how many iterations it need to converge:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2379.2350247213276\n"
     ]
    }
   ],
   "source": [
    "weights_t, records_t = batch_gradient_descent(X, y, 0.0001, 2000000)\n",
    "print(records_t[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimal Hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  lrate    niter    elapsed time      cost\n",
      "-------  -------  --------------  --------\n",
      "    0.5    25000        0.674147   1963.47\n",
      "    0.8    25000        0.674809   1963.47\n",
      "    1.2    25000        0.676915   1963.47\n",
      "    1      25000        0.695837   1963.47\n",
      "    2      25000        0.771415  14550.9\n",
      "    0.5    50000        1.344      1963.47\n",
      "    1.2    50000        1.34606    1963.47\n",
      "    0.8    50000        1.37204    1963.47\n",
      "    2      50000        1.39839   15679.5\n",
      "    1      50000        1.43581    1963.47\n",
      "    0.8   150000        4.06828    1963.47\n",
      "    0.5   150000        4.12251    1963.47\n",
      "    1.2   150000        4.13694    1963.47\n",
      "    2     150000        4.19002   21300.8\n",
      "    1     150000        4.22477    1963.47\n"
     ]
    }
   ],
   "source": [
    "lrates = [2, 1.2, 1, .8, .5]\n",
    "niterations = [1000, 5000, 10000]\n",
    "\n",
    "optimal_rec = test(X, y)\n",
    "\n",
    "time_sorted_optimal_records = sorted(optimal_rec, key=lambda k: k['time'])\n",
    "optimal_tab_data = time_sorted_optimal_records.copy()\n",
    "\n",
    "for r in optimal_tab_data:\n",
    "    r['cost'] = r['history'][-1]\n",
    "    del r['history']\n",
    "    del r['w']\n",
    "\n",
    "print(tabulate.tabulate(optimal_tab_data, headers={'lrate': 'lrate', 'niter': 'niter', 'time': 'elapsed time', 'cost': 'cost'}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rahul/my_project_dir/my_project_env/lib/python3.6/site-packages/numpy/core/fromnumeric.py:90: RuntimeWarning: overflow encountered in reduce\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n",
      "/home/rahul/my_project_dir/my_project_env/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: overflow encountered in square\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "/home/rahul/my_project_dir/my_project_env/lib/python3.6/site-packages/ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in subtract\n",
      "  del sys.path[0]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  lrate    niter    elapsed time      cost\n",
      "-------  -------  --------------  --------\n",
      "    2.5    25000        0.676167    nan\n",
      "    0.8    25000        0.688363   1963.47\n",
      "    2      25000        0.712631  14550.9\n",
      "    1.5    25000        0.764528   1963.47\n",
      "    3      25000        0.828177    nan\n",
      "    1.5    50000        1.32701    1963.47\n",
      "    2.5    50000        1.33529     nan\n",
      "    0.8    50000        1.33731    1963.47\n",
      "    2      50000        1.38057   15679.5\n",
      "    3      50000        1.51449     nan\n",
      "    1.5   150000        4.04906    1963.47\n",
      "    0.8   150000        4.05869    1963.47\n",
      "    2     150000        4.142     21300.8\n",
      "    2.5   150000        4.1771      nan\n",
      "    3     150000        4.30766     nan\n"
     ]
    }
   ],
   "source": [
    "lrates = [3, 2.5, 2, 1.5, .8]\n",
    "niterations = [1000, 5000, 10000]\n",
    "\n",
    "optimal_rec = test(X, y)\n",
    "\n",
    "time_sorted_optimal_records = sorted(optimal_rec, key=lambda k: k['time'])\n",
    "optimal_tab_data = time_sorted_optimal_records.copy()\n",
    "\n",
    "for r in optimal_tab_data:\n",
    "    r['cost'] = r['history'][-1]\n",
    "    del r['history']\n",
    "    del r['w']\n",
    "\n",
    "print(tabulate.tabulate(optimal_tab_data, headers={'lrate': 'lrate', 'niter': 'niter', 'time': 'elapsed time', 'cost': 'cost'}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rahul/my_project_dir/my_project_env/lib/python3.6/site-packages/numpy/core/fromnumeric.py:90: RuntimeWarning: overflow encountered in reduce\n",
      "  return ufunc.reduce(obj, axis, dtype, out, **passkwargs)\n",
      "/home/rahul/my_project_dir/my_project_env/lib/python3.6/site-packages/ipykernel_launcher.py:3: RuntimeWarning: overflow encountered in square\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "/home/rahul/my_project_dir/my_project_env/lib/python3.6/site-packages/ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in subtract\n",
      "  del sys.path[0]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 118748 characters]\n",
      "text/plain": [
       "<Figure size 1080x1080 with 15 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "optimal_rec = test(X, y)\n",
    "plot_records(optimal_rec)"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
