06_stochastic_gradient.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stochastic Gradient Descent (SGD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"width: 100%; overflow: hidden;\">\n",
    "     <div style=\"width: 75%; float: left;\"> \n",
    "     This notebook presents the Stochastic Gradient Descent (SGD), which is a popular algorithm frequently used in the field of machine learning. The example describes a single layer neural network with logistic regression for breast cancer prediction. The proposed model is analytically derived and implemented using the Numpy library to demonstrate the core functionality of training and testing. Nonetheless, a Pytorch equivalent of the model is given further below for validation purposes. <br>\n",
    "     <br>last update: 11/04/2022\n",
    "     </div>\n",
    "     <div style=\"margin-left: 80%;\">\n",
    "        <b>Author</b><br><br>\n",
    "        <a href='http://www.christopherhahne.de', style=\"text-decoration: None; text-align: left;\">\n",
    "            <img align=\"left\" src=\"http://www.christopherhahne.de/images/about.jpg\" width=\"78\"/>\n",
    "            <br><br><br><br><br><br>\n",
    "            <b>Christopher<br>Hahne, PhD</b>\n",
    "        </a>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data acquisition\n",
    "\n",
    "For our classification example, we employ real data using the [UCI ML Breast Cancer Wisconsin (Diagnostic) dataset](https://archive.ics.uci.edu/ml/datasets/Breast+Cancer+Wisconsin+(Diagnostic)), which consists of 569 test persons with 2 classes (malignant and benign) and 30 different measured attributes known as features. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(426, 30) (426,) (143, 30) (143,)\n"
     ]
    }
   ],
   "source": [
    "# import required packages\n",
    "import numpy as np\n",
    "from sklearn.datasets import load_breast_cancer\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "# load data and split it into training and test set\n",
    "bc = load_breast_cancer()\n",
    "train_X, test_X, train_y, test_y = train_test_split(bc.data, bc.target, random_state=42)\n",
    "class_labels = bc.target_names\n",
    "\n",
    "# plot shapes\n",
    "print(train_X.shape, train_y.shape, test_X.shape, test_y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SGD theory\n",
    "\n",
    "Stochastic Gradient Descent (SGD) is a variant of the gradient descent algorithm used to learn weight parameters of a multi-layer perceptron, particularly useful for training large datasets. The weights of the model are updated iteratively. Models trained using SGD are found to generalize better on unseen data.\n",
    "\n",
    "### Optimization\n",
    "\n",
    "For a stochastic regression, a predicted value $\\hat{y}$ is a scalar composed by $\\hat{y}_i=\\mathbf{X}_i\\mathbf{w}$ where a column vector $\\mathbf{w}=\\left[w^{(1)}, w^{(2)}, \\dots, w^{(J)}\\right]^\\intercal$ consists of weights $w^{(j)}$ for each feature $j$ and the row vector $\\mathbf{X}_i=\\left[x_i^{(1)}, x_i^{(2)}, \\dots, x_i^{(J)}\\right]$ represents a data sample with $J$ features while $i$ is the sample index from a set with total number of $N$ samples. Note that we add a feature column of ones $\\mathbf{x}^{(1)}=[1, 1, \\dots, 1]^\\intercal$ to the data set to embed and train the bias as variable $w^{(1)}$ instead of treating it as a separate variable. Our predicted class value $\\hat{y}$ is supposed to match its actual class $y$ for which a least-squares cost metric $(\\hat{y}-y)^2$ may be a reasonable choice. Similar to conventional optimization, SGD aims to minimize an objective function $F(\\mathbf{w})$, which in case of SGD may be defined as\n",
    "\n",
    "$$\n",
    "F(\\mathbf{w})=\\frac{1}{N}\\sum_{i=1}^N f\\left(\\mathbf{w}, \\mathbf{X}_i, y_i\\right)=\\frac{1}{N}\\sum_{i=1}^N\\left(\\hat{y}_i-y_i\\right)^2=\\frac{1}{N}\\sum_{i=1}^N\\left(\\mathbf{X}_i\\mathbf{w}-y_i\\right)^2\n",
    "$$\n",
    "\n",
    "where $f\\left(\\mathbf{w}, \\mathbf{X}_i, y_i\\right)=\\left(\\mathbf{X}_i\\mathbf{w}-y_i\\right)^2$ represents the **loss function**. In the machine learning context, the term **training** refers to an optimization problem where weight parameters $\\mathbf{w}$ are adjusted so that the objective function is minimized which formally \n",
    "\n",
    "$$\n",
    "\\text{arg min}_{\\mathbf{w}} \\, F(\\mathbf{w})\n",
    "$$\n",
    "\n",
    "To achieve this, SGD inherits the **Gradient Descent** update method at iteration $k$ which is known as **back-propagation** and writes\n",
    "\n",
    "$$\n",
    "\\mathbf{w}_{k+1} = \\mathbf{w}_k - \\gamma \\, \\nabla_{\\mathbf{w}_k} f\\left(\\mathbf{w}_k, \\mathbf{X}_i, y_i\\right) \\, , \\, \\forall i\n",
    "$$\n",
    "\n",
    "where $\\gamma$ denotes the learning rate and $\\nabla_{\\mathbf{w}_k} f\\left(\\mathbf{w}_k, \\mathbf{X}_i, y_i\\right)$ is the gradient of the loss function with respect to the weights $\\mathbf{w}_k$. Here, the gradient $\\nabla_{\\mathbf{w}} f\\left(\\mathbf{w}, \\cdot \\right)$ can be generally obtained by\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf{w}} f\\left(\\mathbf{w}, \\mathbf{X}_i, y_i\\right) = \\frac{\\partial f\\left(\\mathbf{w}, \\mathbf{X}_i, y_i\\right)}{\\partial \\mathbf{w}}\n",
    "=\n",
    "\\begin{bmatrix} \n",
    "\\frac{\\partial}{\\partial \\mathbf{w}^{(1)}} f\\left(\\mathbf{w}, \\mathbf{X}_i, y_i\\right) \\\\\n",
    "\\frac{\\partial}{\\partial \\mathbf{w}^{(2)}} f\\left(\\mathbf{w}, \\mathbf{X}_i, y_i\\right) \\\\\n",
    "\\vdots \\\\\n",
    "\\frac{\\partial}{\\partial \\mathbf{w}^{(J)}} f\\left(\\mathbf{w}, \\mathbf{X}_i, y_i\\right) \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix} \n",
    "\\mathbf{X}_i^{(1)} 2\\left(\\mathbf{X}_i\\mathbf{w}-y_i\\right) \\\\\n",
    "\\mathbf{X}_i^{(2)} 2\\left(\\mathbf{X}_i\\mathbf{w}-y_i\\right) \\\\\n",
    "\\vdots \\\\\n",
    "\\mathbf{X}_i^{(J)} 2\\left(\\mathbf{X}_i\\mathbf{w}-y_i\\right) \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and written in compact vector notation\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf{w}} f(\\mathbf{w}, \\mathbf{X}_i, y_i) = 2\\mathbf{X}_i^\\intercal\\left(\\mathbf{X}_i\\mathbf{w}-y_i\\right)\n",
    "$$\n",
    "\n",
    "where $^\\intercal$ denotes the transpose. Iteration through the entire data set $\\forall i \\in \\{1, \\dots, N\\}$ is referred to as one *epoch*. The resulting weights $\\mathbf{w}$ have shown to be improved by letting the optimization procedure see the training data several times. This means that SGD sweeps through the entire dataset for several epochs.\n",
    "\n",
    "### Mini-Batching\n",
    "\n",
    "Completion of a single epoch is often sub-divided in bundled subsets of samples, so-called *batches*, of size $B$ where $B<N$. While $B=1$ in classical SGD, mini-batching requires $B>1$, which helps reduce the variance in each parameter update. The batch size is often chosen to be a power-of-two for better performance from available matrix multiplication libraries. In practice, we determine how many training examples will fit on the GPU or main memory and then use the nearest power-of-two as the batch size.\n",
    "\n",
    "### Activation function\n",
    "In logistic regression, we desire a classification label $\\mathring{y}$ that only has two possible values whereas, so far, our model employs linear combinations $\\hat{y}_i=\\mathbf{X}_i\\mathbf{w}$ that yield results in the $\\hat{y}_i \\in (-\\infty, \\infty)$ range. Thus, we seek a continuous function that maps real numbers $\\hat{y}_i=\\mathbf{\\hat{y}} \\in \\mathbb{R}^N$ to the $(0,1)$ codomain. A function that satisfies this condition is the *sigmoid function*, also known as *standard logistic function*, given by\n",
    "\n",
    "$$\n",
    "\\text{sigmoid}\\left(\\mathbf{\\hat{y}}\\right)=[\\sigma_1, \\sigma_2, ..., \\sigma_N]^\\intercal, \\quad \\text{where} \\quad \\sigma_i=\\frac{1}{1+\\exp(-\\hat{y}_i)}\n",
    "$$\n",
    "\n",
    "The returned value $\\sigma_i \\in (0,1)$ of the activation function is then assigned a predicted label $\\mathring{y}_i \\in \\{0,1\\}$ which is negative if it is closer to 0 and positive in case it is closer to 1, so that\n",
    "\n",
    "$$\n",
    "\\mathring{y}_i=\n",
    "\\begin{cases}\n",
    "    1,  & \\text{if } \\sigma_i \\geq \\tau\\\\\n",
    "    0,  & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "where $\\tau$ is an adjustable threshold scalar. To account for a logistic regression, activated predictions $\\mathring{y}$ are thus substituted for unnormalized results $\\hat{y}$ in the above SGD optimization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SGD implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# activation function\n",
    "sigmoid_activation = lambda y: 1.0 / (1 + np.exp(-y))\n",
    "np.seterr(over='ignore') # ignore overflow from exponential function\n",
    "\n",
    "# learning rate represents the step size as in classical gradient descent\n",
    "l_rate = 1e-6\n",
    "\n",
    "# epochs represents the number of maximum iterations for minimization\n",
    "epochs = 2000\n",
    "\n",
    "# batch size amounts to the number of samples used in one epoch iteration (up to hardware memory)\n",
    "b_size = 2**6\n",
    "assert b_size <= train_X.shape[0]\n",
    "\n",
    "# insert column of ones as first feature entry to cover bias as a trainable parameter within the weight vector (instead of separate variable)\n",
    "train_X, test_X = np.c_[np.ones(train_X.shape[0]), train_X], np.c_[np.ones(test_X.shape[0]), test_X]\n",
    "\n",
    "# initialize weight vector such it has the same number of columns as input features\n",
    "w = np.random.uniform(size=(train_X.shape[1],))\n",
    "\n",
    "# initialize a list to track the loss value for each epoch\n",
    "loss_list = []\n",
    "\n",
    "# batch composition\n",
    "def next_batch(X, y, b_size):\n",
    "    \n",
    "    # loop over our dataset in mini-batches\n",
    "    for i in np.arange(0, X.shape[0], b_size):\n",
    "        \n",
    "        # yield a tuple for current batch of data and labels\n",
    "        yield (X[i:i+b_size], y[i:i+b_size])\n",
    "\n",
    "for epoch in range(epochs+1):\n",
    "    \n",
    "    # reset total epoch loss\n",
    "    epoch_loss = []\n",
    "    \n",
    "    # loop over data in batches\n",
    "    for (batch_X, batch_y) in next_batch(train_X, train_y, b_size):\n",
    "        \n",
    "        # take dot product between current feature batch and weights\n",
    "        y = np.dot(batch_X, w)\n",
    "        \n",
    "        # pass product to sigmoid activation function\n",
    "        preds_y = sigmoid_activation(y)\n",
    "        \n",
    "        # loss determination via difference of true values\n",
    "        loss = preds_y - batch_y\n",
    "        \n",
    "        # compute sum of squared loss for current batch\n",
    "        epoch_loss.append(np.sum(loss))\n",
    "        \n",
    "        # compute gradient from transposed batch data and batch error\n",
    "        gradient = 2*np.dot(batch_X.T, loss)\n",
    "        \n",
    "        # scale gradient of current batch to step in the correct direction\n",
    "        w -= l_rate * gradient\n",
    "        \n",
    "    # update loss list by taking average across all batches\n",
    "    loss_list.append(np.average(epoch_loss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML\n",
    "plt.style.use('seaborn-notebook')\n",
    "\n",
    "# animated figure that plots the loss over time\n",
    "fig, ax = (plt.figure(figsize=(8, 5)), plt.axes())\n",
    "ax.set_xlabel('Epoch #', fontsize=16)\n",
    "ax.set_ylabel('Loss', fontsize=16)\n",
    "ax.xaxis.set_tick_params(labelsize=12)\n",
    "ax.yaxis.set_tick_params(labelsize=12)\n",
    "line, = ax.plot(range(len(loss_list)), loss_list, lw=2)\n",
    "point, = ax.plot(0, np.nan, 'r', marker='.', markersize=14)\n",
    "plt.tight_layout()\n",
    "plt.close()\n",
    "\n",
    "# animation\n",
    "div = 250\n",
    "def animate(i):\n",
    "    line.set_data(np.arange(len(loss_list))[:i*div], loss_list[:i*div])\n",
    "    point.set_data(i*div, loss_list[i*div])\n",
    "    return line, point\n",
    "\n",
    "anim = FuncAnimation(fig, animate, interval=200, frames=epochs//div+1)\n",
    "HTML(anim.to_jshtml())\n",
    "\n",
    "if False:\n",
    "    anim.save('./data/sgd_anim.gif', writer='imagemagick', fps=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAFwCAYAAABHMOt7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd7gU5f3+8fcNSDGxV0AQUewoKGJMNAFFjbHHEkuMqEn8qSkmMcavvRtLjAUwNsQeS4wtsaIYCypiiRIVG9hARQRBacLn98fM0WVZDrt7ztlzeLhf17XX7s48M/PZ5XDv7DPPzigiMDOzNLRq7gLMzKzxONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiULcllqRekoZL+kxSSDq1ibYzMF9/v6ZYf0ry92lYc9exOHOoW81JWlrS0ZIelzRZ0hxJH0n6dx6AbWpQQxvgH0AP4CTgIOCOpt5uc5HULQ/MkHTvQtosJemTvM24Bmxrj6b6gLRFk398ZLUkaR3gX8C6wMPAg8AkYFVgQH47PyKObeI61gVeB/4QERc28bZaA0sBsyNiXlNuq54augHvADPzWrpExISiNnsBt+dtPoqIblVuaxhwcESoimXbA3MjYk412zZo8j0iszqSOgD3At2BvSKieM/4XElbAFvUoJzV8/vJTb2hiJgLzG3q7ZTpXmAPsm8m5xXNOxT4L9Aa+HatCsr/LuZExFcRMbNW202Vu1+sln4OrAf8pUSgAxARoyJiSOG0/Ov8k5K+kDQ9f7x78bKSxkkaIWl9Sf+SNE3SVEm3S1q9oN0I4LH86TUF3RLd6uv/ztc9rmjadyXdJ2mipJmSPsi7kb5T0KbkOiWtLGmwpPckzc7vB0taqahd3fLbSjpG0luSZkkaK+ngUu9jPT4C/g0cUrSNjsCOwDWlFpLUV9KwfJtf5u/tk5L2LH6PgIPzx1FwG5hPG5Y/X0XSUEkfAV8AaxQsM6xgfUfm004q2k6nvKvoVUnfqvA9SJr31K2W9s7vryh3AUlHAoOB14DT88kDgTslHR4RxevqDIwA/gn8EdgUOBxYFtghb3MW8CRwfF7L4/n0T8p/KSBpPeAhYCJwMVlgrgZsnW/36XqWXQ54ClgHGAo8D/QGjgC2ldQ3IqYVLXY20AG4HJiVtx0m6c2IeLKC0oeSvX9bRcTIfNrBZN8mbiD78C22J7A+cCswHlgpX+YOSQdGxE15u7PIdha3Ifs2UOepovXVvW9nAN8CppcqNCKGSNoOOEXSoxHxhKRWwI3AMsCAiPii/Je+BIgI33yryQ34FJhaQfsVyP6zvwksWzB9WeAtYBqwfMH0cUAA+xatZ3A+fb2Caf3yaQOL2g7Mp/crUc8IYFzB89/kbfsu4nUssE6y8AvgyKK2R+XTzyix/AtA24LpncnC/eYy3stu+ToGke3MTQSuKJj/OnB7/viVwteZT/tWiXUunS/3v6Lpw7JoKVnHsLyOGxYyP4BhJf4OxgHv5o9Pytv9qrn/plvizd0vVkvLkgVxubYn24u7JCI+r5uYP76ErN93QNEyH0bErUXTHsnve1RW7iJNze93zw/wVWJPsm8Gxd80Ls+n77nAEjAkImbXPYmID4CxVPi6IuIr4HrgJ5I6SPoe2YHrofUs8/XecD56aSWyUH8E2EDSspXUAFxQQb2fAQcAHYH7gFOAuyNiUIXbXCI41K2WPif7ylyutfL7MSXm1U3rXjT97RJtP83vVyoxryH+TjaC53hgsqRHJP1J0pplLLsW8HoesF/Ln49lwdcFC39t1byua8g+ZPciO0D6IfDAwhpLWlXSFQV94JPIPnz+X95k+Qq3P7aSxhHxFHAusGW+3UMr3N4Sw6FutfQKsKykUoHVWOobZVLOELv6xvjOdwwqImZFxPZkQXNOvu3TgdeKDyA2koW9toqHDkbE/4BnyLp79gWui2yUzoIrl0Q29PRg4FrgJ8APyb5J1fWlV5QlEfFlJe0ltSU7kAuwItC1kuWXJA51q6V/5PelDsSVUrdnulGJeRsWtWksdUMcVywxb60S04iIZyPijDzg1yHbkz1zEdt5G1iv+IdW+fN1afzXVcpQ4Dtk3VgL7XoBNiE78PvniDg2Im6NiAci4mGy4Y/FmuLHL+cAfYBjyb7x/d2jXkpzqFstXUV2YO2YUkMSASRtno94gWyExBfAryUtU9BmGeDXZAdRH2rkGuu6Bebrq5e0P9CpaNrKJZZ/n6x7oNSHQqE7gVVY8APuF/n0f5ZZb0P8HTgN+G1EvFFPu7o9+Pm+EUjamNJ9/9Pz+Yt6D8oiaSfgd8C1EXE+2XDMdckO+loRD2m0momILyXtQvaL0jslPUgWyp+SBVl/sq/Y5+Xtp0g6lmz0yjMF45cHku0RHx4RU2lEEfG6pIeBw/NuhxeBXmTh9SbZrzHrnChpB7If9LxDFnq7kg39K/5hT7HzgH2AwZI2IxvZ0hs4jOyDb1HLN1h+wPnUMpq+SnYM41hJdSNe1iUbKvoysHlR+6eBXwFDJP0LmAM8ExHvVFpjPn7+WuCNfJ1ExL2SLgZ+K+mBiPh7petNmUPdaioi3pTUmywQ9gJOIPv6Pxl4jqzf9qaC9kMkTSAbc35KPvklYM+IuLOJyjwIuBQ4MH/8ONkHzmVkQwPr3Ek2ImNfsvHpM8jC5xfA1fVtICKm5qNOTgN2I9v7/Aj4G3BKLDhGvdlExFxJO5ONWDmYbETSK/njTVkw1G8m+4Daj+yDqxXZ66so1PPx6NeT/8YgIgrHsh8LfB+4XFJVHxip8rlfzMwS4j51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4gvPN1Mll9xpVi9c5fmLsNaqG+3839NW7jRo0dPiohVSs3zX04zWb1zF6765/DmLsNaqK26r9TcJVgL1qa1xi9snrtfzMwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0uIQ93MLCEOdTOzhDjUzcwS4lA3M0tIm+YuwNI34f132bf/ZiXn7bzPgRx39sUAjHtzLNcMOo+xr7zEp598jCQ6d12LH+21P7vvP5Cl2ratZdnWzKZPn84FF5zPc6NG8dxzo5g0aRL/d/wJnHHGmc1dWovmULea2XrATvT74a7zTVuja/evH3884QOmTZnCdjv/mFVW78i8efN4efQzXHLWCYx++nHOuez6WpdszWjSpEmcecbprLHGGvTq1ZuHH36ouUtaLDjUrWa699iAHXffd6Hz+27Tn77b9J9v2p4HHsoyyy3PHTdczbtvv0HX7j2aukxrITp27Mi7731Ap06dGDduHOusvVZzl7RYcJ+61dSsmTOYNXNGRcus3rkLANM+/7wpSrIWql27dnTq1Km5y1jsJBHqkvpJujd/vJuk42q47V6SflSr7S3Obrvucgb07MKAnl3Yf8AW/PPGoSXbzZzxJVMmf8rED97j0fvu4qYrB7HSqquxzvob1rhis8VPct0vEXE3cHcNN9kL6AP8u4bbXKy0atWKzb/7fb6//c6s1rEzkz6eyL233cCFpx7Lh++N56jjTpuv/U1XXso1l57/9fP1e/bi2DP/Srv2HWpdutlip8WEuqRuwP3A08B3gVHANcBpwKrAgXnTi4H2wAzgkIh4vWg9A4E+EfErSWsDNwLfAu4Cjo6Ib0vqB5wKTAI2BkYDP42IkHQysCvQAXgKODyfPgJ4BugPLA8clj8/HeggaWvgnIi4pTHflxSs1mkNLrr2jvmm7bLvQfz2oD249ZrL2GP/gXRe85v+0h/u8RM22fw7TJ0ymedHPs7bY19l+rSptS7bbLHU0rpf1gH+Aqyf3w4AtgaOAY4HXgO2iYjewMnA2YtY38XAxRHRE3i/aF5v4GhgQ6A78L18+qCI2CIiNiYL9l0KlmkTEX3z5U6JiNl5HbdERC8Hevlat27N/ocdxbx58xg98j/zzevUtRt9vvcDttt5T/545oX022k3fn/IPox7c2wzVWu2+Ghpof5ORLwcEfOAMcDwiAjgZaAbsBxwm6RXgL8CGy1ifVsBt+WPbyqa92xEvJ9v68V8/QD9JT0j6WVg26Jt1O1uji5oXzZJp0oKSTHp44mVLp6c1fIDoFM+m1xvu+133Yuv5szhwbtvq7edmbW8UJ9V8HhewfN5ZF1FZwCP5nvRu5J1wzTGtuYCbSS1B4YAe+d791cWbWNWYftKNxgRp0aEIkIrr7p6lWWn44Px7wCwwkor19tu9qzsbZ82dUqT12S2uGtpob4oywEf5I8HltH+aWCv/PF+ZbSvC/BJkr4N7F3GMtOAZcpot8T67NNPFpg2a9ZMrv/bX2ndpg19t+6/0HYAd948DIANNin9q1Qz+0aLOVBapvOAayWdCPyrjPZHAzdIOoHsIGy9R9siYoqkK4FXgIlkB2sX5VHgOEkv4gOlJQ059zTee+dN+nzvB6zasTOTP/mYB+66lffHvc0vfnc8q3VaA4DzT/oDn382mV5bfo9VO3Zm+udTGfXECJ576jE23qwvO+xWzmespWTw4EFMmTKFKVOyb2lPPvkEZ52VnSZg1113Y5NNNmnO8lokZV3WaZK0NDAjH72yH7B/ROze3HUBrN+zV1z1z+HNXUZNPHzPP7j7lusY/9YbfD71M9q370CPDXuy10G/4Ac7fnMcevi//sl9d/ydt14fw5TJn7JU27Z0XWsd+u+0O3sf/EvatWtIb9viZavuKzV3CS3C2t27MX78+JLzrr76Gg4eOLC2BbUQbVprdET0KTUv9VDfBhgECJgCHBoRbzZvVZklKdStcg51q099ob64db9UJCIeBzZt7jrMzGplcTtQamZm9XCom5klxKFuZpYQh7qZWUIc6mZmCXGom5klxKFuZpYQh7qZWUIc6mZmCXGom5klxKFuZpYQh7qZWUIc6mZmCXGom5klxKFuZpYQh7qZWUIc6mZmCXGom5klxKFuZpYQh7qZWUIqDnVJGzRFIWZm1nDV7KmPkTRS0uGSlmv0iszMrGrVdr/0BYYAEyTdJGkHSWrEuszMrArVhPrngPJbe+AnwH3Au5LOlNSjEeszM7MKVBPqqwC7AMOAyXwT8J2B/wNek/SEpMMkfauxCjUzs0WrONQjYk5E/DsiDgVWA3YCrgIm8U3AbwVcAXwg6ZhGrNfMzOrRoCGNETEXeAZ4GZhYNzm/F7AscK6D3cysNqoOdUkDJN0MfAhcBGxUNwt4C7gAeDd//v8aWKeZmZWhTaULSDoFGAh0rZuU388D7gUGR8SDedvrgZcK2pqZWROqONSBU8i6WOrCfBJwNfC3iBhf1Pbt/L51deWZmVklqgl1yAL9WWAwcEtEzF5IuxlA/yq3YWZmFaom1IeRdbGMXlTDiJgHPFbFNszMrAoVh3o+lNHMzFogn6XRzCwhDnUzs4Q41M3MEuJQNzNLiEPdzCwh1Y5Tn4+kNsDuwIrAPRExcRGLmJlZE6jmcnbnSRpV8FzAw8CtwOXAy5LWbrwSzcysXNV0v/wQeLzg+a7A94HzgQPyacc1sC4zM6tCNd0vXYA3Cp7vCrwTEccBSNoIOLARajMzswpVs6feFviq4Hl/su6XOm8DHRtSlJmZVaeaUH+P7MpGdXvl3Zn//C6rAtMbXpqZmVWqmu6XvwMnSVqV7MIYnwP/Lpjfm+wiGWZmVmPV7KmfQ3amxq3Izqv+s4iYAiBpOWA3YHhjFWhmZuWr5iyNs4DD8luxaWT96V82sC4zM6tCo/z4qE5+/vSpjblOMzMr3yJDXdLPqllxRFxXzXJmZla9cvbUhzH/NUnLEYBD3cysxsoJdV9j1MxsMbHIUI8IX2PUzGwx4VPvmpklpOrRL5L6AFsCK7Dgh0NExBkNKczMzCpXcahL6gDcAexAdvC08CBqFExzqJuZ1Vg13S8nkwX6WWQHUQUcDOxEdkreUcCGjVWgmZmVr5pQ3xu4LSJOBl7Jp30QEQ8AA8jO4jiwccozM7NKVBPqXfjmrIxz8/u2ABHxFXAzsF/DSzMzs0pVE+rT+KYvfhowD+hUMH8qsHoD6zIzsypUE+pvAesCRMRcYAxZl0zd9Up/THbOdTMzq7FqQv1hYC9JrfPnlwM/lPQW2WXuBgBXN1J9ZmZWgWrGqf8ZuJ58GGNEDJHUHvgpWR/7lcB5jVZhopZu24bNui7f3GVYC/XAEy83dwm2mKrmfOrTgdeLpl0IXNhYRZmZWXV8mgAzs4RU84vSR8poFhGxXRX1mJlZA1TTp96d7DQAxevpSLbnPwn4ooF1mZlZFarpU+9WarqkdsDvgUOAHzSsLDMzq0aj9alHxKyIOAd4Bh80NTNrFk1xoPQJYMcmWK+ZmS1CU4T6WuTngjEzs9qqZvRL14XMWpHs16S/AUY0oCYzM6tSNaNfxrHg6Jc6Ivth0m+qLcjMzKpXTaifzoKhHsBkYCzwcETMa2hhZmZWuWqGNJ7aBHWYmVkjqPhAqaShkrasZ35fSUMbVpaZmVWjmtEvA4G165m/Ftk1S83MrMaaYkjjt4A5TbBeMzNbhLL61PNhjN0KJq0v6fslmq4IHAG82fDSzMysUuUeKD0EOIVslEsAJ+S3YiK7ZukhjVKdmZlVpNxQv5NsfLqAocAVwMiiNgFMB0ZFhK9RambWDMoK9Yh4CXgJQNKawB0R4ettmZm1MNWMUz+tKQoxM7OGq2ac+lGSHq5n/oOSDm9YWWZmVo1qx6m/Uc/8scChVVVjZmYNUk2o9wDq608fk7cxM7MaqybUlwLa1zO//SLmm5lZE6km1McC29czfwfgrerKMTOzhqgm1G8GdpB0hqSvr3AkaSlJp5GF+k2NVaCZmZWvmvOp/xXYiewXpUdIei2fvj7ZaQIeB/7SOOWZmVklKt5Tj4g5ZHvjxwHvA73z23vAscCAiJjdmEWamVl5qtlTrwv28/LbAiS1i4hZDSnMzMwq16in3pW0uaQhwIeNuV4zMytPVXvqhSStCPyU7AdHPclO+jW2oes1M7PKVb2nLmlHSbcAH5AdPG0HnAb0jIj1G6k+MzOrQEV76pK6ke2RHwysAUwCbgcOAE6IiDsauT4zM6tAWXvqkg6UNJzsikZ/Ap4D9gQ6A6eSdbmYmVkzK3dP/XrgbeBo4OaI+LRuhuQ8NzNrKcrtU59Fdo3S3YEfSurQZBWZmVnVyg31jmR76SuR7bVPlHR1fvFp76qbmbUQZYV6REyJiEERsRnQB7iBrE/9UeAJsuuTLtdkVZqZWVmqOU3A8xFxFNne+0Fk508HuErSi5JOlLRRYxZpZmblqXqcekTMioibImI7YG3gLGAF4HTyi1SbmVltNcppAiJiXEScTHYw9UeAx6ubmTWDBp8moFBEBHB/fjMzsxpr1BN6mZlZ83Kom5klxKFuZpYQh7qZWUIc6mZmCXGom5klxKFuZpYQh7qZWUIc6mZmCXGom5klxKFuZpYQh7qZWUIc6mZmCXGom5klxKFuZpYQh7qZWUIc6mZmCWnUKx+ZlTJ9+nQuuvACRj/3HM+Pfo5JkyZx7HH/xymnnbFA26+++ooLzvsz1w0bxsSJE1hzzW4cfsSRHH7EkUhqhuqtqV14zokMv//uhc4/6Oe/Yr+DfgnAxxM/5NorL+H5USOZMeML1ujSjd33OYjtd9q9VuW2eA51a3KfTprEOWedSec11mCTTXvxyPCHF9r2t78+imFDr+aQQ39Ony22YPjDD/GH3/2Wzz6bzP+dcFINq7Za2WnXfei1+XcWmH737Tfyxutj6LPl1gBM+uQjfnfEgcyZPZtdf7w/K6y4Ms8+9RgX/fkkvpj+OXvsc1CtS2+RHOrW5Fbv2JE333mXjp06MX7cODZcb52S7V566UWGDb2a3xz9O84593wABh56GDpwf84/988MPPTndOzYsZalWw1ssPGmbLDxpvNNmzlzBkP+ehbduvdgnXU3BOC2G4cy9bPJnD/ouq/b77Lnfpx+/G+4/upBbLvDriy73PI1r7+lcZ+6Nbl27drRsVOnRba74/bbADjyqF/PN/2Io37FrFmzuPfuu5qkPmt5Rj7+CDO+/ILtfrjb19PG/Hc0HTt1WeADoP/2uzBzxgxGPvFIrctskZIJdUndJL3SCOvpI+mSxqjJKvP86NGstvrqdOnadb7pm/fZglatWvHCC883U2VWa8Pvv5vWrdvQf/tdvp42Z85s2rVvv0Db9h06APDma2NqVl9LlkyoN5aIeC4iftPcdSyJJkyYQMeOC+7Rt23blhVXWokPP/igGaqyWpv0yUe89PwzbNb3u6yw4kpfT1+j61q8/944Jn86ab72/33h2Wy5SR/XtM6WKrVQbyPpRkmvSrpd0tKSNpf0mKTRkh6Q1BFA0ghJ50p6VtJYSdvk0/tJujd/vIqkhySNkXSVpPGSVs6/Fbwq6cp83oOSOjTnC0/BzBkzaNeubcl57du3Z+bMGTWuyJrDIw/ey7x58xhQ0PUCsMse+zFn9mzOPvn3/O+VF5k44X3uvv1G7rsr67abNWtmc5Tb4qQW6usBQyJiA+Bz4CjgUmDviNgcGAqcVdC+TUT0BY4GTimxvlOARyJiI+B2oLBfoAcwOJ83BdirsV/MkqZ9hw7MmjW75LyZM2fSvr0/N5cEjzxwD8ssuxxbfrfffNN7b7EVvzn2VN4d9xZ/POpnHLbfj7jhmiEc+bsTAFi6w7eaodqWJ7XRL+9FxJP54xuA44GNgYfyMc6tgQkF7e/I70cD3Uqsb2tgT4CIuF/SZwXz3omIFxex/HwknUr+4bG6R3EsoGPHjowZs+BhkdmzZzP500/LOthqi7exr77Ce+PfZuc9fsJSbRf81rbjzj9m2+134Z23xjJ37ly691iPjydm/6U7dVmz1uW2SKntqUfR82nAmIjold96RsQOBfNn5fdzqfwDblbB47KWj4hTI0IRoVJ9x0u63pttxkcTJ/Leu+/ON330c6OYN28evXtv1kyVWa0MfyD7EdJ2RV0vhZZq25Z1N9iYDTbelHbt2vPCqKcA2GyLrWpSY0uXWqh3lVT3L3sA8DSwSt00SUtJ2qiC9T0J7JsvuwOwQmMWa/P78d77ADBk8KXzTb9s8CDatm3Lrrv5V4MpmzNnDo8Nv48ua3ZnvQ16lrXM5E8/4babhrLOehuy6WZbNnGFi4fUul9eB46SNBT4H1l/+gPAJZKWI3u9FwHljn06DbhZ0kHASGAi2d7/txu78NT9bchgpk6dwpQpUwAY+dSTnHtOdnjjR7vsSs+em9CrV29+NvAQLr34IqZPm/71L0r/cfttHH/iSe5+SdyzIx9j2udT2Wv/Q0rOn/zpJE459ki22qY/K62yGp98NIH77rkdIjjmhHN8GolcMqEeEeOA9UvMehH4fon2/QoeTyLvE4+IEcCIfNZUYMeI+Crf298iImYB48j66uuWv6DhryBtF190Ie+OH//188f/8x8e/89/AOjUeQ169twEgEsGDaFLly5cf9213HD9tay5ZjfO/8tfOeKoXzVL3VY7w++/m1atWrHtDruUnN+hw9Ks3qkz99/7D6Z+Nplll1uBvlt9nwMHHsHKq65e42pbLkUUd0NbHUk9gFvJuqlmA0dGxKjGWPdmm/eJJ0Y+0xirsgSNePp/zV2CtWA7/2CT0RHRp9S8ZPbUm0JEvAH0bu46zMzKldqBUjOzJZpD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4hD3cwsIQ51M7OEONTNzBLiUDczS4giorlrWCJJ+gQY39x1tCCdgA+buwhrsfz3Mb81I2KVUjMc6tYiSIqIUHPXYS2T/z7K5+4XM7OEONTNzBLiULeW4rTmLsBaNP99lMl96mZmCfGeuplZQhzqZmYJcaibmSXEoW5mlhCHujU5Sf0k3Zs/3k3ScTXcdi9JP6rV9qw0Sd0kvdII6+kj6ZLGqClVbZq7AFuyRMTdwN013GQvoA/w7xpu05pIRDwHPNfcdbRk3lO3suR7Wq9JGiZprKQbJQ2Q9KSkNyT1zW8jJb0g6SlJ65VYz0BJg/LHa0t6WtLLks6UND2f3k/SCEm359u8UZLyeSdLGiXpFUlXFEwfIelcSc/m9W0jqS1wOvATSS9K+knt3jEroU3+b/lq/m+7tKTNJT0mabSkByR1hNL/nvn0wm99q0h6SNIYSVdJGi9p5fxv9VVJV+bzHpTUoTlfeC051K0S6wB/AdbPbwcAWwPHAMcDrwHbRERv4GTg7EWs72Lg4ojoCbxfNK83cDSwIdAd+F4+fVBEbBERGwMdgF0KlmkTEX3z5U6JiNl5HbdERK+IuKWK12yNZz1gSERsAHwOHAVcCuwdEW3P0/0AAAhTSURBVJsDQ4GzCtrP9+9ZYn2nAI9ExEbA7UDXgnk9gMH5vCnAXo39Yloqd79YJd6JiJcBJI0BhkdESHoZ6AYsB1wrqQcQwFKLWN9WwB7545uACwrmPRsR7+fbejFf/xNAf0nHAksDKwJjgHvyZe7I70fn7a1leS8inswf30C2I7Ax8FD+has1MKGg/aL+PbcG9gSIiPslfVYw752IeHERyyfJoW6VmFXweF7B83lkf0tnAI9GxJ6SugEjGmlbc8m+urcHhgB9IuI9SacC7UssMxf/bbdExT9fnwaMiYitFtK+If+exX8/7n4xq8JywAf544FltH+ab74W71dG+7oAnyTp28DeZSwzDVimjHbW9LpKqgvwA8j+/VepmyZpKUkbVbC+J4F982V3AFZozGIXVw51a0znAedIeoHy9qyOBn4v6b9k/fVT62scEVOAK4FXgAeAUWVs41FgQx8obRFeB46S9CpZAF9K9sF8rqSXgBeB71awvtOAHfKhkvsAE8k+xJdoPqGXNRtJSwMz8n75/YD9I2L35q7LFg+S2gFzI+KrfG//sojo1dx1NTf3O1pz2hwYlA9LnAIc2sz12OKlK3CrpFbAbOAXzVxPi+A9dTOzhLhP3cwsIQ51M7OEONTNzBLiUDczS4hHv5hZk5P0M7Jz+IyNiJuau56UefSLmTUpSd8lO28PZKeU2DIiRjdjSUlz94tZA0kaJ2lEwfN+kkLSwEbeTpOstynlpz++iuz0Ed8BPgGukuRegibiULfFTkG4Fd6m5+fk/q2k1s1dY7WUXanp1PyEaCk4EVgT2C0iniU7K+f6wB+btaqEufvFFjuS+pGd0+VmsisaCehEdhKx9YErI+KXNaxnHDAuIvrlz1sBbYE5ETG3wnUNBK4B+kfEiKJ5Va/Xlhz+CmSLs+cj4oa6J5IuA14Ffi7ppIj4qHiBfC++XUR82VRFRcQ8YObisl5Li7tfLBkR8TkwkmzPvXt+6bxQdtm9kyS9RRaK+9YtI6mdpOPzy57NlDRF0j2SehevX1IXSbdKmirp87zd2iXalez7ltRW0rH5GSO/zNfznKRf5fNPJdtLB3i0oGtp2CLWu7KkwZLekzQ7vx8saaWidnXvx7aSjpH0lqRZyi4Xd3CJ19E+7wp6Pa93irJLD55f37+DpGsKaj+txHswpWB+JafatTJ4T92SkZ8YbJ386SSyy6dBdkWlpchO2/s52SlgkbQUcD/Z6V6vBwaRnRP+F8CTkr6fX+gYScsD/wG6AH8D/gf8gKwbaJEXYMgPGD4A9AMeJLvyz0ygJ/DjfNt3AB2BX5JdCvDVfPG36lnvcsBT+eseCjxPdinAI4BtJfWNiOLT0Z6d13w52cUkjgCGSXqz4MpEAIPJTrJ2HXAhWV70ALZdxMsdzDfn0x8o6bT8WwbA9mTvMcCoiBiziHVZpSLCN98WqxtZMAbZ9UdXBlYBNiEL7QBG5u0G5s9fB5YusZ7f5fN3LJq+LPAuMKJg2tl520OK2l6UTx9Ror6BBdOOzaedXaKOVgWP62ruV8/rLlzvWfm0I4vaHpVPP6PEul8A2hZM70wW7jcXrWMy8O8q/42eybcVwA4F068pmH5kNev2rf6bu19scXYa2RC5j4GXyPYq7+ab657WuSxK96H/lOxi2aPzLoyVJa1MdjDyIWBrfXMV+j2Aj8j2WgudW2atBwKfAacXz4hv9mKrsSfZe3BF0fTL8+l7llhmSGQX5a7b/gfAWLK98EJTgY0kbVxFXYMLHh8GkA9j3C2fNovsQLc1Mne/2OLsCuA2sr2+L8h+rTi5RLuxC1l+A7JuiE/q2cbKwHtkv4YcFUWjTiJigqQpZdTaA3gxIhr7QOdawHMR8VVRXV9JGgtsVmKZt0tM+5Rs6GGho8m6pV6W9DZZV9M9wD1lfBDdAvyF7P3bXdKKwBZkFwsHuCsiPlvYwlY9h7otzt6IiIfLaLewkS4CXgZ+X8+y9QX+4mphwyFV+CQi7srHy/+I7PjBALK97sclDSjc2y8WEbMkXQ38CWhH9q2oZ0GTa0ouaA3mULcl2Rtk/fGPlLHn+TbQQ1Lrwr11SR2B5cvY1lhgfUntImJWPe0q/eHI28B6ktoU7q3nXR3rUnqvvGz5N58bgBvyA9F/Jjs+sDvZt6T6/I3sR0atyA4+r55P/5Cse8uagPvUbUl2HVnQlNxTl7RawdO7gNWAnxU1+1OZ27qR7GLLJ5bYTuEe8vT8fsXidgtxJ9kH08+Lpv8in/7PMtdTXFPrfMTP1yI70vlCufVFxDjgX/nTjcm6YgCuK+7GssbjPXVbkl1MNsTufEnbAo+QDXnsCmxHNuSwf972POAA4EpJmwNjyEajbEU2fLKcbe0KnChpC7JhjTOBjciGXg7I240iO+nVCZJWIDtW8E5EPLOQ9Z4H7AMMlrQZWej2JusmeT2fX41lgAmS7s7X+TFZ//0RZAd87ylzPYPJXnehYVXWZGVwqNsSKyLmSNoZOBI4iGw0DWTdA88C1xa0/UzSNmTjtev21h8jC/3hZWxrtqQdgD+QfTicTRbqb1DQvxwR70o6lOwbwGVk4+uvJRsiWGq9UyV9L699N+AQslE6fwNOiQXHqJfrS7LhmtuRfeB8G5hANrronIj4sMz1PAi8yTe/HxgZEa9XWZOVwed+MbMmlf8C9Zj86eERUTz80hqRQ93MGl1+jp0OZF0295J1aU0DOjfg24OVwd0vZtYUtiEb117oLAd603Oom1lTmkt2yoXLqf6grVXA3S9mZgnxOHUzs4Q41M3MEuJQNzNLiEPdzCwhDnUzs4Q41M3MEuJQNzNLiEPdzCwhDnUzs4T8f1KQGBfpEH/ZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compute predictions from test set\n",
    "pred_y = np.round(sigmoid_activation(np.dot(test_X, w))).astype('uint8')\n",
    "\n",
    "# compose confusion matrix\n",
    "from sklearn.metrics import confusion_matrix\n",
    "conf_mat = confusion_matrix(y_true=test_y, y_pred=pred_y)\n",
    "\n",
    "# print confusion matrix using matplotlib\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.matshow(conf_mat, cmap=plt.cm.Blues, alpha=0.3)\n",
    "for i in range(conf_mat.shape[0]):\n",
    "    for j in range(conf_mat.shape[1]):\n",
    "        ax.text(x=j, y=i, s=conf_mat[i, j], va='center', ha='center', size='xx-large')\n",
    "\n",
    "plt.title('Confusion Matrix', fontsize=18)\n",
    "plt.xlabel('Predictions $\\mathbf{\\hat{y}}$', fontsize=18)\n",
    "plt.ylabel('Actuals $\\mathbf{y}$', fontsize=18)\n",
    "plt.xticks([0, 1], class_labels)\n",
    "plt.yticks([0, 1], class_labels)\n",
    "plt.tick_params(top=False, bottom=True, labeltop=False, labelbottom=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "84.13 % Precision\n",
      "98.15 % Recall\n",
      "92.31 % Accuracy\n"
     ]
    }
   ],
   "source": [
    "# TP/(TP+FP)\n",
    "precision = lambda conf_mat: conf_mat[0, 0]/(conf_mat[0, 0]+conf_mat[1, 0]) * 100\n",
    "print('%.2f %% Precision' % precision(conf_mat))\n",
    "\n",
    "# TP/(TP+FN)\n",
    "recall = lambda conf_mat: conf_mat[0, 0]/(conf_mat[0, 0]+conf_mat[0, 1]) * 100\n",
    "print('%.2f %% Recall' % recall(conf_mat))\n",
    "\n",
    "# (TP+TN)/ALL\n",
    "accuracy = lambda conf_mat: (conf_mat[0, 0]+conf_mat[1, 1])/np.sum(conf_mat) * 100\n",
    "print('%.2f %% Accuracy' % accuracy(conf_mat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overfitting validation\n",
    "\n",
    "In statistical learning, we want to make sure that the classification performs equally well on test and training data. Therefore, we employ the Mean-Squared Error (MSE) given by\n",
    "$$\n",
    "\\text{MSE}(\\mathbf{\\hat{y}}, \\mathbf{y})=\\frac{1}{N}\\sum_{i=1}^N \\left(\\hat{y}_i-y_i\\right)^2\n",
    "$$\n",
    "on both sets while aiming for\n",
    "$\\text{MSE}(\\mathbf{\\hat{y}}_{\\text{test}}, \\mathbf{y}_{\\text{test}}) \\approx \\text{MSE}(\\mathbf{\\hat{y}}_{\\text{train}}, \\mathbf{y}_{\\text{train}})$. If this fails, it gives indication for either under- or overfitting of the trained weights $\\mathbf{w}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSEs are close enough with 0.0769 (test) and 0.1103 (train).\n"
     ]
    }
   ],
   "source": [
    "# mean squared error\n",
    "MSE = lambda y, pred_y: np.round(sum((y-pred_y)**2)/len(y), 4)\n",
    "\n",
    "# compute predictions of test and training sets\n",
    "pred_test_y = np.round(sigmoid_activation(np.dot(test_X, w))).astype('uint8')\n",
    "pred_train_y = np.round(sigmoid_activation(np.dot(train_X, w))).astype('uint8')\n",
    "\n",
    "# compare MSEs\n",
    "test_MSE = MSE(test_y, pred_test_y)\n",
    "train_MSE = MSE(train_y, pred_train_y)\n",
    "res = np.isclose(test_MSE, train_MSE, rtol=.95)\n",
    "\n",
    "if res:\n",
    "    print('MSEs are close enough with %s (test) and %s (train).' % (test_MSE, train_MSE))\n",
    "else:\n",
    "    print('Potential over-/underfitting from MSEs with %s (test) and %s (train).' % (test_MSE, train_MSE))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyTorch equivalent\n",
    "\n",
    "Neural network models are often defined using PyTorch's *Module* class, which offers inheritance from Object-Oriented Programming (OOP). Variables of other types (e.g. numpy) have to be converted to torch data types to enable the convenient automatic gradient computation. The model, optimizer (here SGD) and loss function are instantiated before training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch 0\n",
      "            Train set - loss: 296.118, accuracy: 62.238\n",
      "            Test  set - loss: 314.943, accuracy: 62.238\n",
      "            \n",
      "epoch 500\n",
      "            Train set - loss: 3.766, accuracy: 76.923\n",
      "            Test  set - loss: 2.314, accuracy: 76.923\n",
      "            \n",
      "epoch 1000\n",
      "            Train set - loss: 3.276, accuracy: 81.119\n",
      "            Test  set - loss: 1.865, accuracy: 81.119\n",
      "            \n",
      "epoch 1500\n",
      "            Train set - loss: 3.017, accuracy: 83.217\n",
      "            Test  set - loss: 1.651, accuracy: 83.217\n",
      "            \n",
      "epoch 2000\n",
      "            Train set - loss: 2.845, accuracy: 83.916\n",
      "            Test  set - loss: 1.503, accuracy: 83.916\n",
      "            \n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "\n",
    "# define single layer model\n",
    "class SingleLayerNet(torch.nn.Module):\n",
    "    def __init__(self, n_features):\n",
    "        super(SingleLayerNet, self).__init__()\n",
    "        self.linear = torch.nn.Linear(n_features, 1, bias=False)\n",
    "        torch.nn.init.uniform_(self.linear.weight, 0, 1)\n",
    "\n",
    "    def forward(self, X):\n",
    "        z = self.linear(X)\n",
    "        return torch.squeeze(z, 1)\n",
    "    \n",
    "# convert to torch data types\n",
    "train_Xt = torch.autograd.Variable(torch.FloatTensor(train_X))\n",
    "train_yt = torch.autograd.Variable(torch.FloatTensor(train_y))\n",
    "test_Xt = torch.autograd.Variable(torch.FloatTensor(test_X))\n",
    "test_yt = torch.autograd.Variable(torch.FloatTensor(test_y))\n",
    "\n",
    "# instantiate model and loss\n",
    "model = SingleLayerNet(n_features=train_Xt.shape[1])\n",
    "optimizer = torch.optim.SGD(model.parameters(), lr=l_rate)\n",
    "criterion = torch.nn.BCEWithLogitsLoss()\n",
    "\n",
    "# training\n",
    "loss_list = []\n",
    "for epoch in range(epochs+1):\n",
    "    # loop over data in batches\n",
    "    for i_X, i_y in next_batch(train_Xt, train_yt, b_size):\n",
    "        optimizer.zero_grad()\n",
    "        pred_y = model(i_X)\n",
    "        loss = criterion(pred_y, i_y)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "    # track loss\n",
    "    loss_list.append(loss.detach())\n",
    "    if (epoch) % 500 == 0:\n",
    "        y_test_pred = model(test_Xt)\n",
    "        y_train_pred = model(train_Xt)\n",
    "        test_loss = criterion(y_test_pred, test_yt)\n",
    "        train_loss = criterion(y_train_pred, train_yt)\n",
    "        pred_yt = np.round(sigmoid_activation(y_test_pred.detach().numpy()), 0)\n",
    "        conf_mat = confusion_matrix(y_true=pred_yt, y_pred=test_yt.numpy())\n",
    "        test_acc = accuracy(conf_mat)\n",
    "        train_acc = accuracy(conf_mat)\n",
    "        print(\n",
    "            f'''epoch {epoch}\n",
    "            Train set - loss: {np.round(train_loss.detach().numpy().astype('float'), 3)}, accuracy: {np.round(train_acc, 3)}\n",
    "            Test  set - loss: {np.round(test_loss.detach().numpy().astype('float'), 3)}, accuracy: {np.round(test_acc, 3)}\n",
    "            ''')\n",
    "\n",
    "if False:\n",
    "    torch.save(model, 'torch_sgd_model.pth')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0wAAAFoCAYAAABkP3u+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZxcVZnw8d+TpJN0B0g3e8IWNtlXgUFFB2RRUEFHZHMQcBsVRB19RxQGooKIsggCigqCI6DAoDCKIIKAguwii8gmRELYAkmA7Mvz/lG3k0qnO72kupau3/eT+vStc8899dTt6ko9dc49JzITSZIkSdKyhtU6AEmSJEmqVyZMkiRJktQDEyZJkiRJ6oEJkyRJkiT1wIRJkiRJknpgwiRJkiRJPTBhkiRJUtVExMSIeDEiMiKOrEB7E4q2dqpAeHUrInYvnufqtY6l2ZgwSZIkNbmIWCsizo6IpyJibkQ8FxG/jYj9Kvw4WwMnAZ8CxgG/qECzzxZtPVCBtnpUlrC8FhFtXfZtUezrV0ITERdHxK/7WP0OSs/zlX6ErQoYUesAJEmSVDsRMQG4HXgd+ArwV0pfqu8J/ABYv4IPt0nx81eZmZVoMDMXAi9Uoq0+mg58CLikrOxjwD+p7LlaLCJaMnMe1X2eKtjDJEmS1NzOL37ulJlXZOZjmfloZp4LbNtZKSLWj4hfRsTrxe3qiFi3bP/EiHg4Ig4peqpej4hfdfa4RMRE4JdF9UURkUX5Mr0snW2V3d8mIm4qenfeiIi/RsQexb5lhuRFxDsi4q6ImFMM/zsrIkaW7b8lIs6PiG9GxNSIeCkiTo+Ivnw2vhj4aFlbLcDhRXn5cxgeERdGxNMRMTsinoiI/+p8jOJ8HAG8p6x3avey53NoRNwcEbOB/+g6JK9o+5GIaC17vD/2o8dKfWTCJEmS1KQiYlXg3cB5mflG1/2ZOb2oNwy4BlgL2KO4jQd+FRFRdsgE4GDgA8A+wA7AKcW+04FPFNvjiltfXQY8D+wCbA9MBOb08JzWAX4L/KV4/I8BhwKndqn6YWAB8FbgGODzRey9+RmwS0RsXNx/L/AGcEuXesOA54CDgC2A44GvAkcV+08HrgB+z5LzcUfZ8adSSma3BH7VTRzHAi1FOxTtb0pZMqfKcEieJElS89oECODRXurtSam3aePMfAYgIg4Dniz2/b6oNwI4MjNnFHV+SJEgZOYbETG92O7v0LINgNMz8+/F/SeXU/czwBTgM5m5CHg0Io4DLoiI/87MWUW9v2XmicX24xHxieK5XN5LLK8C11JKTI6nlJD9BFhqiGFmzgdOLCt6JiJ2pJS8XVicj9nA3PLzUZZ/fi8zryor36SsLTJzZkR8GLg9Il6hNJxy/8x8qZf41U/2MEmSJDWv6L0KUOohmdKZLAFk5j8oJSZbltWb1JksFaYAa65okMCZwI+LIWrHR8TmvcR6Z5EsdfoTMJIl11ABPNjluP7EeiHwkYhYD9ibLsPxOkXEpyLi3oh4OSLeAL5A369zure3Cpl5D6UevP8GfpiZv+1j2+oHEyZJkqTm9QSlnpEtVqCN8p6V+d3s6+3z5iKWTdxalmokcyJLhqa9FXgwIgYy9GxFY+30e0px/xS4OTMnd60QEQcD36WUTL2L0lDC8yklbn0xs7cKxXDI3YCFwMZdhkeqQkyYJEmSmlRmvgrcABwTESt13R8R7cXmo8D4Yka9zn0bUbqO6W8rGMbLLHs90/bdxPpEZp6Tme+h1MPz8R7aexTYtcsEDrsB84CnVjDWzlgWUUqEdi9i6c5uwF2ZeW5m3p+ZTwIbd6kzDxi+AqH8J7Aj8A5gV+CzK9CWemDCJEmS1NyOptTDc29EfCgiNouIzSPi0ywZtvb7YvvSiNipmJHuUuB+4OYVfPybgR0i4qMRsUlE/Bfwts6dEdEaEeeVzSD3L5SSkZ4StfMpJXLnF+sjvQf4FnBu2fVLlXAysAZwdQ/7Hwd2jIh9I2LTiPhv4F+71HkG2Lo456sXM+71SURsR2k43icy8w5K126dFhFb9feJaPlMmCRJkppYcS3SjsCNwGmUEqObgf2BTxZ1EjiAUm/QH4rbC8D7V3Q9pcy8AfgapQ//91Gaae/8sioLgQ5KPTqPUZqa/M+Uele6a+85YF9KM+Q9AFxEaSKHr65InN08zvzMnNrlWqlyF1CaBe8y4B5Kz+uMLnV+RKlH7F5K5/Zt9EFEjKaUsF6WmVcX8VwGXAVcFhGj+vdstDxRoTXDJEmSJGnIsYdJkiRJknpgwiRJkiRJPTBhkqosIi6OiF/385hbIuLcwYpJkiRJ3fMaJqkHEdHbH8clmXnkANodS+lvb3o/jlkVmJ+Zr/f38fojIiYCB2bm1oP5OJIkSY1iRK0DkOpY+ZoQ76U0k0152ezyyhHRkpldF8FbRpcV0PukWCdDkiRJVeaQPKkHmflC5w2YXl4GjAamR8ShEXFzRMwG/iMiVouIyyNickTMjohHIuKo8na7DskrhtudHxHfjIipEfFSRJxevuBe1yF5EfFMRJwQERdExGvF4/2/Lo/zpoi4NSLmRMRjEbFfRLwREUcO9JxExDYR8fviub1aPJexXfbfVMT0RkT8NSL2KPa1RMQ5ETElIuZGxLMR8a2BxiJJklQNJkzSijmV0loRWwK/opRI3U+pR2or4GzggojYs5d2PgwsAN4KHAN8Hji4l2O+ADxEae2M04BvR8RbAIpk65dFm7sCRwInAQNelyEixlBaDf4NYBfgA0W8F5VVuwx4vti/PTARmFPsO7Y45hBg0+L5PTbQeCRJkqrBIXnSivleZl7Vpew7Zds/jIh3AocCNy2nnb9l5onF9uMR8QlgT0oL7fXkd5nZ2ev0vYg4tjjmz8DewGbAPsUCfkTEF4Db+/KkenAYMAY4vPNaqoj4JPCHiNgkM58ENgBOz8y/F8c8WXb8BpRWPf9jscjhP4E7ViAeSZKkQWcPk7Ri7i2/ExHDI+L4iHgwIl6JiDeAfwPW76WdB7vcnwKsuQLHbA5M6UyWCvcAPa1G3hdbAA92mXjijqLNLYv7ZwI/LoYpHh8Rm5fVvZhSr9PjEXFeRLynfNihJElSPfLDirRiZna5/yXgi5R6mfaklCD8ChjZSztdJ4tIev/7HMgxgyUBMnMiS4YnvhV4MCI+Wuy7H5gAfIVSnJcAN5o0SZKkeuYHFamydgP+LzP/JzMfAJ4C3lSDOP4OjI+I8WVlO7Fif/OPAttExMplZW8t2ny0syAzn8jMczLzPcCFwMfL9r2emVdl5qeB9wDvBDZZgZgkSZIGldcwSZX1OHBwROwGTAU+C2wI/KXKcdxIaUKFSyLiS0ArpeFyCyh6g5ZjdERs36VsFnAp8DXgpxFxItABXABcnZlPRkQrcDpwJfAMsBalBPIugIj4T0oTQjxAqXfsMOA1YPIKPVNJkqRBZMIkVdbJlBKk31Jap+liSonGlss5puIyc1FEfAD4MXA3pQTmi8DVLJm1ricbs2yCd19m7hQR7wK+W7Q5B7gG+FxRZyGlJOpiSutVvQL8mtIwRYDXgf9HaYa8LB5j38ycNaAnKUmSVAVRmqxK0lAXEdtR6t3ZKTPvq3U8kiRJjcCESRqiih6mmcATlCZbOBMIYIf0D1+SJKlPHJInDV0rU1rQdj1gGnAL8AWTJUmSpL6zh0mSJEmSeuC04pIkSZLUgyE9JG/11VfPCRMm1DoMSdIgu++++6Zm5hq1jkO117Hqajl+vQ1qHYbqVGuLfQXq3oMPPfTavLlzx3a3b0gnTBMmTODee++tdRiSpEEWEZNqHYPqw/j1NuDy395a6zBUp7Zae+XeK6kprbH6qi/1tM80W5IkSZJ6YMIkSZIkST0wYZIkSZKkHpgwSZIkSVIPTJgkSZIkqQcmTJIkSZLUAxMmSZIkSeqBCZMkSZIk9cCESZIkSZJ6YMK0HPdNepU/P/VKrcOQJEmSVCMmTMvxmUvv58v/+2Ctw5AkSZJUIyZMy9HRNpJps+bVOgxJkiRJNWLCtBxjW1t4fc4CFixcVOtQJEmSJNWACdNydLSNBGDG7Pk1jkSSJElSLZgwLUd7WwsA002YJEmSpKZkwrQc7UUP03SvY5IkSZKakgnTcnR09jDNsodJkiRJakYmTMvROSRvmgmTJEmS1JSqljBFxHoR8YeI+FtEPBIRnyvKV42IGyPiieJnR1EeEXFORDwZEQ9GxI7VirWTQ/IkSZKk5lbNHqYFwBczc0tgV+DoiNgSOA64KTM3BW4q7gPsC2xa3D4JfL+KsQLQ3uqQPEmSJKmZVS1hysznM/P+Yvt14FFgHeAA4JKi2iXA+4vtA4CfZsmdQHtEjKtWvAAdY0o9TC5eK0mSJDWnmlzDFBETgB2Au4C1MvP5YtcLwFrF9jrAs2WHTS7KqmZxD5PTikuSJElNqeoJU0SsBPwv8PnMfK18X2YmkCvY/sSIyIjIKVOmrEhTjF08S549TJIkSVIzqmrCFBEtlJKlSzPz6qL4xc6hdsXPl4ry54D1yg5ftyhbrsycmJmRmTF+/PgVinfUiOG0jRzuNUySJElSk6rmLHkBXAg8mplnlu26Fjii2D4CuKas/CPFbHm7AjPKhu5VTUfbSBMmSZIkqUmNqOJjvQ04HHgoIh4oyr4KfAu4IiI+BkwCDir2XQfsBzwJzAKOqmKsi41tbWHSKzNr8dCSJEmSaqxqCVNm/gmIHnbv2U39BI4e1KD6oGNMC397fiHzFixi5AjX+ZUkSZKaiRlALxYvXjvbiR8kSZKkZmPC1AsXr5UkSZKalwlTLzqKHqZpM+1hkiRJkpqNCVMv2ttcvFaSJElqViZMvVh8DZOL10qSJElNx4SpF17DJEmSJDUvE6ZedIwpJUzTTJgkSZKkpmPC1IuxraUheTOcVlySJElqOiZMvegoJn2YNtMeJkmSJKnZmDD1YmznNUz2MEmSJElNx4SpFyOGD2Pl0SOc9EGSJElqQiZMfdDe1mLCJEmSJDUhE6Y+6GgbyTTXYZIkSZKajglTH7S3jWTugkXMnrew1qFIkiRJqiITpj5od+IHSZIkqSmZMPWBU4tLkiRJzcmEqQ/GtpUWr7WHSZIkSWouJkx90NnD5Ex5kiRJUnMxYeqDdhMmSZIkqSmZMPVBezEkz6nFJUmSpOZiwtQHnbPkzZhtD5MkSZLUTEyY+qCjs4dppj1MkiRJUjMxYeqDxdcw2cMkSZIkNRUTpj5YZXQLwwKmew2TJEmS1FRMmPpg2LBgbGsL05wlT5IkSWoqJkx91N420mnFJUmSpCZjwtRH7W0tTJ81j8ysdSiSJEmSqsSEqY/aW1tYsCiZOW9hrUORJEmSVCUmTH3k1OKSJElS8zFh6qOxbS5eK0mSJDUbE6Y+WtzD5NTikiRJUtMwYeqjxYvXOlOeJEmS1DRMmPqovehhcvFaSZJUb+76061st84qbLfOKvzz6ad6rPfLn//P4noLFiyoYoSqpjfeeIOJE0/ive/Zj7XXWoMRw4P//u8Tuq27YMECTj75G2yy8YaMaRvNVltuznnnnevM0GVMmPqovdUeJkmSVH/mz5/Pqcd/kda2McutN2Paq5x9yom91lPjmzp1Kid/4+s8/PBDbL/9Dsute/RnPs3Ek05kr7325pxzzmWbbbblc8d+lpNP/kaVoq1/Jkx9tOQaJhMmSZJUP376g3OYMX0aHzzsiOXWO/vUr7HqGmuyx7vfU6XIVCvjxo3jn88+xzOTnuUHF/ywx3oPPPAAF174Y77wn1/kBxf8kI99/OP8/BdX8KEPHcS3Tv0mzz//fBWjrl8mTH205Bomh+RJkqT68Pxzz/LDs7/D5776NVZaZZUe6z30l3v55eWX8OWvf5sRw0dUMULVwqhRoxg/fnyv9a688goAPvvZY5cqP+azxzJ37lyu+dWvBiW+RmPC1EeLEyanFZckNYGI2D0ifl1s7x8Rx1XxsbePiP2q9XiN7LQTv8ymm2/JAQd9uMc6ixYt4ptf/SJ77XcA/7Lbv1YxOtW7++69l7XXXpv1119/qfKdd96ZYcOGcf/999UosvriVwx9tNKoEYwYFk4rLklqOpl5LXBtFR9ye2An4LoqPmbDue3313Pr767jZ7++mYjosd4VP/0x/3jiMc768aVVjE6N4Pnnp3TbEzVy5EhWW201nnvuuRpEVX/sYeqjiKC9rYUZXsMkSWoQETEhIv4eERdHxOMRcWlE7BURt0fEExGxS3H7c0T8JSLuiIjNumnnyIg4t9jeOCLujIiHIuLkiHijKN89Im6JiKuKx7w0ik/xEXFiRNwTEQ9HxA/Lym+JiNMi4u4ivrdHxEjg68DBEfFARBxcvTPWOObOmcNpJ/4X7z/kcLbabsce670y9WXO/fbJfPyzX2TtddatYoRqBLNnz2bkqFHd7hs9ejSz58yuckT1yYSpH9rbRtrDJElqNJsAZwCbF7fDgN2ALwFfBf4OvD0zdwBOBL7ZS3tnA2dn5jbA5C77dgA+D2wJbAS8rSg/NzN3zsytgVbgvWXHjMjMXYrjTsrMeUUcv8jM7TPzF709wYiYGBEZEfnyi81xkfqF557BazOmc+xXJi633lnfOIH2jg6O+NSxy62n5tTa2sq8uXO73TdnzhxaR7dWOaL6ZMLUD+2tLcyYPZ9Fi5yXXpLUMJ7OzIcycxHwCHBTlhZYeQiYAIwFroyIh4GzgK16ae8twJXF9mVd9t2dmZOLx3qgaB9gj4i4KyIeAt7Z5TGuLn7eV1a/XzJzYmZGZsYaa40bSBMN5eUXX+An53+XA//9o8ya+QbPPTuJ556dxOszZgDw0gvP88Jzk3n0oQf4v6su598/cTQvv/TC4nqzZs4E4PnJ/+TlF1+o5VNRjY0bN54pU6YsUz5v3jxeeeWVPk0c0Qy8hqkf2ttGsijh9TkLGFtMAiFJUp0r//p4Udn9RZQ+B3wD+ENmfiAiJgC3VOixFgIjImI0cD6wU2Y+GxETgdHdHLMQP5f0ySsvv8S8uXO56NwzuejcM5fZ/7ED96O9Y1VOOv1cAE49/kucevyXlqn33rdtz5bb7sDlv7110GNWfdrxzW/m97+/kX/+859LTfxwzz33sGjRInZ885trGF398I2pH5bMlDfPhEmSNFSMBTqv7D6yD/XvBD4I/AI4pA/1O5OjqRGxEnAgcFUvx7wOrNyHtpvSOutvwFkXdu3cg+uvuYobrr2a4089i/HrrsdmW23bbb3LLvwB99xxG2f86Ge0d6xajZBVpz70oYP49mnf4nvfO4fvfOf0xeXnfu8cRo4cyQEHvL+G0dUPE6Z+6CiSpGmz5rPBajUORpKkyvg2cElEnAD8pg/1Pw/8LCKOB64HZiyvcmZOj4gfAQ8DLwD39OEx/gAcFxEPAKf25TqmZrLyKmN557vfu0z5Y488CMCub9+d9TfcGKDben+4/tcA7L7PfowY4UfBoeq8885l+vTpTJ8+HYDbb/8Tp5xyMgDve9/+bLvttuywww4cddRH+e5ZZ/LG66+z8867cOONv+PKK6/gv088ySF5Bf9K+qG9bSTg4rWSpMaQmc8AW5fdP7KHfW8qO+yEYv8tFMPzMvNi4OJi/3PArpmZEXEIsFnX+sX9Y8q2T+hst0t8u5dtT6W4hikzXwV27tuzlNSdM884nUmTJi2+f9utt3LbraXhl+uusy7bbrstAOd//west/76XHLxT7jkkouZMGECZ333bI455rM1ibseRem6z6Fpp512ynvvvbdi7V161ySO/+XDfPfg7Xn/DutUrF1J0oqJiPsyc6dax9EMIuLtwLlAANOBj2bmk7WNaomtttsxvSZHPdlqbUd6qntrrL7qk6+++uqm3e2zh6kf2ltLPUxOLS5JalaZ+Udgu1rHIUnV4rTi/dB5DdN0F6+VJEmSmoIJUz94DZMkSZLUXEyY+mHJtOL2MEmSJEnNwISpHzraOq9hMmGSJEmSmkHVEqaIuCgiXoqIh8vKJkbEcxHxQHHbr2zfVyLiyYh4LCLeVa04l2d0yzBGjhjGDIfkSZIkSU2hmj1MFwPv7qb8rMzcvrhdBxARW1JaPXyr4pjzI2J41SLtQUTQ0dZiD5MkSZLUJKqWMGXmbcCrfax+APDzzJybmU8DTwK7DFpw/dDeOtJJHyRJkqQmUQ/XMB0TEQ8WQ/Y6irJ1gGfL6kwuynpVDPPLiMgpU6ZUOlba21p4bc4CFixcVPG2JUmSJNWXWidM3wc2BrYHngfOWNEGM3NiZkZmxvjx41e0uWV0zpT32pwFFW9bkiRJUn2pacKUmS9m5sLMXAT8iCXD7p4D1iurum5RVnNLZspzWJ4kSZI01NU0YYqIcWV3PwB0zqB3LXBIRIyKiA2BTYG7qx1fd8Z2rsVkwiRJkiQNeSOq9UARcTmwO7B6REwGTgJ2j4jtgQSeAf4DIDMfiYgrgL8BC4CjM3NhtWJdns4epunOlCdJkiQNeVVLmDLz0G6KL1xO/VOAUwYvooHpKHqYnFpckiRJGvpqPelDwxnb2tnD5JA8SZIkaagzYeqnjsXXMNnDJEmSJA11Jkz91N55DdNse5gkSZKkoc6EqZ+8hkmSJElqHiZM/dQ5rfgMEyZJkiRpyDNh6qdRI4bTNnK4C9dKkiRJTcCEaQDaW1uc9EGSJElqAiZMA9DeNtJpxSVJkqQmYMI0AO1tLcyct5B5CxbVOhRJkiRJg8iEaQA6nFpckiRJagomTAPQ7uK1kiRJUlMwYRoAEyZJkiSpOZgwDUDnkDynFpckSZKGNhOmARjb6uK1kiRJUjMwYRoAe5gkSZKk5mDCNACLr2GabQ+TJEmSNJSZMA1Ae+e04vYwSZIkSUOaCdMAOEueJGmwRMRlEbF3REStY5EkmTANSHsx6YPXMEmSBsEhwPXAPyPimxGxWa0DkqRmZsI0ACOGD2PlUSPsYZIkDYY/AQmsA3wZ+FtE3BkR/xER7bUNTZKajwnTALWPaTFhkiRVXGa+g1KydDRwC7AI2AU4H5gSEZdHxE61i1CSmosJ0wB1tI10SJ4kaVBk5ouZ+X3gXZQSp5nFrtHAwcCdEXFUreKTpGZiwjRAY1tbmLtgEXPmL6x1KJKkISYitoiI04HJwPeBMcWu3wPfo9TrNLE20UlScxlR6wAaVfnitePGttY4GknSUBERf6Y0BA8ggBnAJcD5mfl4UWcC8N5axCdJzcaEaYDKpxY3YZIkVdC/FD8fBM4DLs3MWV3q/Bp4tapRSVKTMmEaoPayHiZJkiro58B5mXl7TxUy80fAj6oXkiQ1LxOmAepci2mGM+VJkiooMw+rdQySpCWc9GGAOsZ0Ll5rwiRJkiQNVSZMA9TeWhqSN322Q/IkSZKkocqEaYDKJ32QJEmSNDSZMA3Q4kkfZtrDJEmqnIjYMiI2K7u/d0T8LCK+EhHDaxmbJDUjE6YB6ujsYZptD5MkqaIuAnYAiIj1gGuAVYGjgZNrGJckNSUTpgFaeXQLETDdacUlSZW1OXB/sX0gcFdm7gccDhxas6gkqUmZMA3Q8GHB2NYWr2GSJFXacKDz27g9geuK7aeAtWoSkSQ1MROmFdDRNtJpxSVJlfYw8OmIeDulhOn6onwdYGrNopKkJrXCCVNEtFQikEY0trWFGbPnkZm1DkWSNHR8GfgEcAtweWY+VJTvD9xdq6AkqVn1K2GKiGMj4oNl9y8EZkfEY+Uz+jSLjrYW5i9MZs5bWOtQJElDRGbeBqwBrJ6ZHy3bdQHw6dpEJUnNq789TMcCLwNExDuAg4DDgAeAMyobWv3rnFrciR8kSZWUmQszc1qXsmcy86VaxSRJzWpEP+uvAzxdbL8PuDIzr4iIh4A/VjSyBlC+eO26HTUORpLUsCLi2r7Wzcz9BzMWSdLS+pswvQasCTwL7A18pyifD4yuYFwNob21s4fJiR8kSSvklVoHIEnqXn8Tpt8BP4qI+4FNgN8W5VuxpOepaXSMKfUwTXNIniRpBWTmUbWOQZLUvf5ew3Q0cDuli1EPzMxXi/IdgcsrGVgjGNvaOSTPhEmSJEkaivrVw5SZrwGf7ab8pIpF1EA62hySJ0mqvIgYAewCrA+MLN+XmT+tSVCS1KT6lTBFxJbAwsx8rLi/N3AE8Ajw7cxsqvm1Oyd9cPFaSVKlRMTmwP8BGwIBLKT0//V8YC5gwiRJVdTfIXkXATsARMR6wDXAqpSG6p1c2dDq3+IeptkOyZMkVcx3gfuAscAsYAtgJ0pLeHxwOcdJkgZBfxOmzYH7i+0Dgbsycz/gcODQSgbWCMqnFZckqUJ2Bk7OzJnAImBEZt4P/BdNuOahJNVafxOm4UBnd8qewHXF9lPAWpUKqlGsNGoEI4aFkz5IkiopKPUsQWmx+HWK7cmUZqiVJFVRf6cVfxj4dET8mlLC9JWifB1gaiUDawQRQXtbiz1MkqRKehjYDvgHcDfw5YhYCHwCeLKWgTWC0S3D2HytlWodhurUDX96qNYhqE5Ne21Wj/v628P0ZUpv2LcAl2dm56tuf0pv6k1nbGsL02ebMEmSKuYUSr1MACdQminvD8A+wLG1CkqSmlV/pxW/LSLWAFbJzGlluy5gyfCBptLRNpKnp85k0aJk2LDo/QBJkpYjM28o2/4HsEVErApMy8ysXWSS1Jz6OySPzFwYEbMjYmsggacy85mKR9Yg2ttaWJTw+twFixeylSSpksoWipckVVl/12EaAZwKHENpIb0A5kbE94DjM7PHsWkRcRHwXuClzNy6KFsV+AUwAXgGOCgzp0VEAGcD+1HquTqymCGo7rQvXrx2ngmTJGmFRcS1y9ufmftXKxZJUv+vYfo28O/Ap4A3AZsCn6Y0rfipvRx7MfDuLmXHATdl5qbATcV9gH2LtjcFPgl8v59xVk17q4vXSpIq6pUut9coLWL7DppwgiVJqrX+Dsk7DPhoZl5XVvZURLwM/Bj4Uk8HFtc/TehSfACwe7F9CaXJJL5clP+0GKt9Z0S0R8S4zHy+n/EOuo4xS3qYJElaUZl5VHflEXEGpeRJklRF/e1hGktpzaWungLaB/D4a5UlQS+wZC2ndYBny+pNZsk6FMsVERMjIiMip0yZMglEvo8AACAASURBVICQ+qdzGJ5Ti0uSBtkFwNG1DkKSmk1/E6a/0v2Upp8r9g1Y0Zu0wrP/ZObEzIzMjPHjx69oc73qaLOHSZJUFZvVOgBJakb9HZL3X8B1EbEXcGdRtiswntJ1R/31YudQu4gYB7xUlD8HrFdWb92irO50tHkNkySpciLinK5FwDhK/89eVP2IJKm59auHKTNvozTZw1XASsXtSuBdDGwxvWuBI4rtI4Bryso/EiW7AjPq8folgLFFwjTDxWslSZWxTZfblsAC4AvFTZJURQNZh2kKcHx5WURsB3xwecdFxOWUJnhYPSImAycB3wKuiIiPAZOAg4rq11GaUvxJStOKd3sBbD3oHJI3zSF5kqTKOAKYnJmLyguLJTfWA/5Zk6gkqUn1O2EaqMw8tIdde3ZTN2mQC1vb25z0QZJUUU9TGoL3UpfyVYt9w6sekSQ1sf5O+qAuWluGM3LEMCd9kCRVSvRQvhIwp5qBSJKq2MM0VEUE7a0tTPcaJknSCiib7CGBb0bErLLdw4FdgAeqHpgkNbk+JUwRcW0vVVapQCwNq6NtJM/PmF3rMCRJjW2b4mcAWwDlQxfmAfcDp1c7KElqdn3tYXqlD/ufXsFYGtbYthYee/F1FixcxIjhjnKUJPVfZu4BEBE/AT6Xma/VOCRJEn1MmDKzbmepqwedazG9NmcBq44ZWeNoJEmNzP9zJam+2B1SAe2tTi0uSaqMiDglIj7VTfmnIuIbtYhJkpqZCVMFtI9xanFJUsUcDvylm/L7gI9UORZJanomTBXQ2cPk1OKSpApYE3i5m/JXgLWqHIskNT0TpgrocPFaSVLl/BN4ezfl7wAmVzkWSWp6rsNUAe1tXsMkSaqYC4CzImIkcHNRtidwKnBazaKSpCZlwlQB7UUP0wwXr5UkraDMPCMiVgfOATqnXp0HnJ2Z365dZJLUnBySVwEd9jBJkiooM78CrA7sWtzWyMzjImKv2kYmSc3HhKkC2r2GSZJUYZk5MzPvAaYAn4+Ip4AbahyWJDUdE6YKGNtqwiRJqpyIGB4R/xYRvwGeAT5A6dqmTWoamCQ1Ia9hqoDRLcNpbRnukDxJ0gqJiM2Aj1Nab2kmcBmwD3B4Zv6tlrFJUrOyh6lCOtpa7GGSJA1YRPwRuBPoAA7KzI0y84QahyVJTc+EqULGto104VpJ0op4C/BT4KzMvLXWwUiSSkyYKqSjrYWZ8xYyb8GiWociSWpMO1MaKv+niPhLRHwhItaudVCS1OxMmCpk8Ux5s+1lkiT1X2b+JTOPBsYBZwL7A89S+r/6PRHRUcv4JKlZmTBVSHuxFtMMr2OSJK2AzJyTmf+TmXsAWwDfAb4AvBARv61tdJLUfEyYKqSj6GGaZsIkSaqQzHwyM48D1gMOAhzGIElV5rTiFdLeWuphcuIHSVKlZeZC4JriJkmqInuYKmTxNUz2MEmSJElDhglThXRew+SkD5IkSdLQYcJUIV7DJEmSJA09JkwVsmRInj1MkiRJ0lBhwlQhi4fk2cMkSZIkDRkmTBUytrVzSJ49TJIkSdJQYcJUIS3Dh7HyqBH2MEmSJElDiAlTBY1tazFhkiRJkoYQE6YK6mgb6bTikiRJ0hBiwlRB7W0tzJm/iDnzF9Y6FEmSJEkVYMJUQc6UJ0mSJA0tJkwVtGTxWoflSZIkSUOBCVMFtbd2Ll5rD5MkSZI0FJgwVdCSIXn2MEmSJElDgQlTBbUvHpJnD5MkSZI0FJgwVVBHZw+TU4tLkiRJQ4IJUwWNbfMaJkmSJGkoMWGqoA6vYZIkSZKGFBOmCuqcJc9rmCRJkqShwYSpglZpbSECZpgwSZIkSUOCCVMFDR8WjG1tceFaSZIkaYgwYaqw9tYWps+2h0mSJEkaCkyYKqy9bSTTZ80jM2sdiiRJkqQVZMJUYe1tLcxfmMyat7DWoUiSJElaQSZMFdY5tbjXMUmSJEmNz4Spwsa2unitJEmSNFSMqHUAQ82SxWtNmCRJ0uB44403OPP073Dvvfdy3733MHXqVI77ylf52jdOXqreo48+yinf+Dr3338fLzz/PMOGDWOjjTfmI0ccySf/41OMHDmyRs9Ag+nMU0/gpuuv7XH/4R8/hkMO/+Ti+69MfYlLf/J97r3rj8yYPo2x7R1stsU2fOG4b9A2ZqVqhFzXTJgqrL2tc/Fah+RJkqTBMXXqVE45+Rusu+66bLf9Dtz0+xu7rTf52Wd59dVXOeigg1ln3XVZuHAhf77jDr70n1/glj/8gauu/mWVI1c17Pu+D7H9m3ddpvzaqy7licceYad/2W1x2bOTnua4zx1Fa9sY9n3fh1ht9TWZPv1VHn3oL8yZM8eECROmiutMmJxaXJIkDZZx48bx9D8nM378eJ555hk222Sjbuvtvc8+7L3PPkuVferTn6G9o50fnH8+jz32GJtttlk1QlYVbbH1dmyx9XZLlc2ZM5vzzzqFCRttyiZv2hKAzOQ7Jx/HamusxWln/4TWtrZahFv36uIapoh4JiIeiogHIuLeomzViLgxIp4ofnbUOs6+aO8ckjfTHiZJkjQ4Ro0axfjx4wd8/AYbTABgxvTpFYpI9e7Pf7yZ2bNmsue7919c9tf77+Kpxx/l34/6DK1tbcydO4cFC/zSv6u6SJgKe2Tm9pm5U3H/OOCmzNwUuKm4X/c67GGSJA0hETEhIh6uQDs7RcQ5lYhJ/Tdr1iymTp3KpEmT+N+rruTM07/DuHHj2GbbbWsdmqrkpuuvZfjwEeyx93sXl91/9x0AtLaN4UtHf4R/22cXPrD3zhz3uY/xzFOP1yrUulNPCVNXBwCXFNuXAO+vYSx91t7qtOKSJHWVmfdm5rG1jqNZnfGdb7PO2mvypo035LBDDmb9DTbgmv/7Da2trbUOTVUw9eUX+ev9d7HjLm+lY9XVFpc/N3kSAN888Yu0d6zGcRNP5+NHf4ln/vE4X/7cR5n60gu1Crmu1EvClMDvIuK+iOicsmOtzHy+2H4BWKsvDUXExIjIiMgpU6YMRqzL1T6m1MM0w1nyJElDx4iIuDQiHo2IqyKiLSLeHBG3Fv933xAR4wAi4paIOC0i7o6IxyPi7UX57hHx62J7jWK4/SMR8eOImBQRqxe9WY9GxI+Kfb+LCD/RV8CHD/8I113/O3522eV87BOfoGVEC9Mdjtc0bv7dr1m0aBF7lQ3HA5g9exYAG27yJk44+Szevsc+HHDgv3PCyd/ljddf4+orflqLcOtOvSRMu2XmjsC+wNER8Y7ynZmZlJKqXmXmxMyMzIwVGds7UCuPGsHwYWEPkyRpKNkMOD8ztwBeA44GvgccmJlvBi4CTimrPyIzdwE+D5zUTXsnATdn5lbAVcD6Zfs2Bc4r9k0HPthbcOVflj5fgy9LG8FGG23EnnvtxYcOOpjzv38B/3bggbxn33fx6KOP1jo0VcHNN/wfK68yln956+5LlY8aOQqAd+7z3qXKt97uzay59nge/ut91QqxrtVFwpSZzxU/XwJ+CewCvFj2bdU44KXaRdh3EUF7a4vXMEmShpJnM/P2YvtnwLuArYEbI+IB4ARg3bL6Vxc/7wMmdNPebsDPATLzemBa2b6nM/OBXo5fSvmXpeNq8GVpIzrk0MOYP38+l1/6s1qHokH2+KMP8+ykf/COd76bli7rbq26+poAtHestsxxHauuxhuvv1aVGOtdzROmiBgTESt3bgP7AA8D1wJHFNWOAK6pTYT9197W4sK1kqShpOsoj9eBR4rJmrbPzG0ys3zu6rnFz4X0fwmTuWXbAzlefTBnzhwApk2f1ktNNbqbbigtYLtnl+F4AG/afCsAXnn5xWX2TX35Rca2N8Qk1YOu5gkTpWuT/hQRfwXuBn5TfNv0LWDviHgC2Ku43xDa20YyfdY8Fi3q0yhCSZLq3foR8ZZi+zDgTmCNzrKIaImIrfrR3u3AQcWx+wB+KhskL73U/QCdH13wAwB23nmXaoajKps/fz633vRb1ttgIzbbYptl9u+62x6MGjWaG35zNQsXLlxcfvcdt/LKyy+x485vq2a4davm39pk5j+A7bopfwXYs/oRrbiOthYWJbw+dwFjW1tqHY4kSSvqMUrXGF8E/I3S9Us3AOdExFhKnye+CzzSx/a+BlweEYcDf6Y0udPrwEqVDnwoO/+8c5kxffriyRvuuP12Tj3lZADe+7792WbbbTn605/i1Vde4R3/+q+su956TJ8+nd/feCM33/R73vKWt3LoYR+u5VPQILv7z7fy+msz+OChR3W7f2z7qvz7x47mwvPP4Cuf/xhv3+NdvDL1Ja7930tZa9w6vP+gw6sccX2qecI0FI0tphafPmueCZMkqaFl5jPA5t3segB4R9fCzNy9bHsqxTVImXkLcEuxawbwrsxcUPRS7ZyZc4FnKF0b1Xn86Sv+DIau7555BpMmTVp8/7bbbuW2224FYJ1112WbbbfloIMP5qeXXMLFP7mIl19+mVGjRvGmzTbjlFO/xTGfPZaWFj+nDGU3XX8tw4YNW2ZSh3L/dvARrLxKO9dc+T9c+P0zaG0dw26778ORn/wcK6+8ShWjrV8mTINg8eK1s+azwbLX0EmS1OzWB66IiGHAPOATNY6nIT3+1NO91vnQQQfzoYMOrkI0qkcnfrNva0Xvve8B7L3vAYMcTeMyYRoE7UXC5NTikiQtKzOfAHaodRyS1Bf1MOnDkNPeVhqSN8OpxSVJkqSGZsI0CBb3MM20h0mSJElqZCZMg6Cj6GFy8VpJkiSpsZkwDYL2skkfJEmSJDUuE6ZB0HkN03QnfZAkSZIamgnTIOhYPEuePUySJElSIzNhGgStLcMZOXyYPUySJElSgzNhGgQRQXtbi5M+SJIkSQ3OhGmQtLe1OK24JEmS1OBMmAZJe9tIXpuzgIWLstahSJIkSRogE6ZB0t5amvhhhsPyJEmSpIZlwjRIOpxaXJIkSWp4JkyDpN2pxSVJkqSGZ8I0SDoXr50x2x4mSZIkqVGZMA2SxYvXzrSHSZIkSWpUJkyDpHNInmsxSZIkSY3LhGmQtDvpgyRJktTwTJgGyZJJH0yYJEmSpEZlwjRIlkwr7pA8SZIkqVGZMA2SscXCtSZMkiRJUuMyYRoko1uG09oynOlOKy5JkiQ1LBOmQdTe1uK04pIkSVIDM2EaRO1tI5nhtOKSJElSwzJhGkTtrS28MXcB8xYsqnUokiRJkgbAhGkQdYwpTfxgL5MkSZLUmEyYBpGL10qSJEmNzYRpELV3Ti1uD5MkSZLUkEyYBlHn4rXTZtrDJEmSJDUiE6ZBNLbNxWslSZKkRmbCNIg6e5hcvFaSJElqTCZMg6i96GGaZg+TJEmS1JBMmAZRh0PyJEmSpIZmwjSIxrY6rbgkSZLUyEyYBlG7PUySJElSQzNhGkQtw4ex0qgRTLOHSZIkSWpIJkyDrL2thRkuXCtJkiQ1JBOmQdbe1mIPkyRJktSgTJgGWUfbSObMX8Sc+QtrHYokSZKkfjJhGmTtnYvXOvGDJEmS1HBMmAbZqsVMea/MnFvjSCRJkiT1lwnTIFtzldEAvPjanBpHIkmSJKm/TJgG2dpFwvTCDHuYJEmSpEZjwjTI1rKHSZIkSWpYJkyDbO2xowATJkmSJKkRmTANsvHtrUTAM6/MrHUokiRJkvppRK0DGOraRo5gw9XH8MiU11i4KBk+LAB45Y25TJs1n5lzF7BgUS63jQgIICKKnxAEEeX7S/eX2i72sdT95bTRzb7i8C5tLl2P6GeMXet1VpIkSZLqjAlTFbxlo9W49K5/cuzlf2H+wkU8OHkGLzhEbxnLTcjoOemi/H43bSydMC7bxuLH7kv7Xdqga/kyyWb5MX18bstrv6/PbZn2u29j2WR46Tb6fO661utL+8ucn2Xb6LH93s5PX9rv5vx0Pt5SP4ml7i95vRblXY/rUr/rfnrc30N7Xdqlr8f18Dj0GmeF4+/x/C3ZP2bUcDZYbQySJNUjE6YqOOptG3LlfZP5zUPPA7DmyqPYa4s1WWPl0YwZOZyWET2PjMyEJCn+kZlF2ZJ9mZ11c5nyzvt03u9m3zLtLy5fcp/y45bXfg9t0DWuLm2wTJtLt0F3z2157S8Tew/npzh3ixZlWZ0+Pre+tN/Lc5MEb9tkNS79+K61DkOSpG7VfcIUEe8GzgaGAz/OzG/VOKR+22TNlbjpP/+Vp6fO5E1rrczaY0fXOiTVkaWT4GWTLlg2CVsq6VpuslmeDHdTry/td5Ngd016+91Gf2Kky3PsT4wsL2Ffuv2ilcV1y4rL7ne/v7Ngmfo9lbP0fpbZ3784lhxf4fj7GAc9PK++xrHRGvYuSZLqV10nTBExHDgP2BuYDNwTEddm5t9qG1n/rbdqG+ut2lbrMFSHOoeQFfdqGYokSZK6qPdZ8nYBnszMf2TmPODnwAE1jkmSJElSk6j3hGkd4Nmy+5OLsh5FxMSIyIjIKVOmDGpwkiRJkoa2ek+Y+i0zJ2ZmZGaMHz++1uFIkiRJamD1njA9B6xXdn/dokySJEmSBl29J0z3AJtGxIYRMRI4BLi2xjFJkiRJahJ1PUteZi6IiGOAGyhNK35RZj5S47AkSZIkNYm6TpgAMvM64LpaxyFJkiSp+dT7kDxJkiRJqhkTJkmSJEnqgQmTJEmSJPUgMrPWMQyaiHgZmLSCzYwHGmUFXGMdPI0Ur7EOjkaKFRor3krEukFmrlGJYNTYKvR//1DSSO8Fqi5fG0vr8f+RIZ0wVUJEZGZGrePoC2MdPI0Ur7EOjkaKFRor3kaKVWo0/n2pJ742+s4heZIkSZLUAxMmSZIkSeqBCVPvvlbrAPrBWAdPI8VrrIOjkWKFxoq3kWKVGo1/X+qJr40+8homSZIkSeqBPUySJEmS1AMTJkmSJEnqgQmTJEmSJPXAhEmSJEmSemDCJEmSJEk9MGHqQUS8OyIei4gnI+K4OohnvYj4Q0T8LSIeiYjPFeUTI+K5iHiguO1XdsxXivgfi4h31SDmZyLioSKue4uyVSPixoh4ovjZUZRHRJxTxPtgROxYxTg3Kzt/D0TEaxHx+Xo5txFxUUS8FBEPl5X1+zxGxBFF/Sci4ogqx/udiPh7EdMvI6K9KJ8QEbPLzvEPyo55c/H6ebJ4ThVfjbyHWPv9e6/G+0UPsf6iLM5nIuKBorzW57Wn96u6fd1KjSYido+IXxfb+w/We08Pj719+XujaqN4r3+495q9trNTRJxTiZiGrMz01uUGDAeeAjYCRgJ/BbascUzjgB2L7ZWBx4EtgYnAl7qpv2UR9yhgw+L5DK9yzM8Aq3cp+zZwXLF9HHBasb0f8FsggF2Bu2r4u38B2KBezi3wDmBH4OGBnkdgVeAfxc+OYrujivHuA4wotk8ri3dCeb0u7dxdPIcontO+VYq1X7/3ar1fdBdrl/1nACfWyXnt6f2qbl+33rw12g3YHfh1jR77SODcWp+DZr8t773eW2Vv9jB1bxfgycz8R2bOA34OHFDLgDLz+cy8v9h+HXgUWGc5hxwA/Dwz52bm08CTlJ5XrR0AXFJsXwK8v6z8p1lyJ9AeEeNqEN+ewFOZOWk5dap6bjPzNuDVbmLoz3l8F3BjZr6amdOAG4F3VyvezPxdZi4o7t4JrLu8NoqYV8nMO7P0v8JPWfIcBzXW5ejp916V94vlxVr0Eh0EXL68Nqp4Xnt6v6rb161UC0UPwd8j4uKIeDwiLo2IvSLi9qJXdZfi9ueI+EtE3BERm3XTzpERcW6xvXFE3Fn0JJ8cEW8U5btHxC0RcVXxmJd29jBHxIkRcU9EPBwRPywrvyUiTouIu4v43h4RI4GvAwcXPdgHV++MqRsjit/lo8Xvtq0YSXBrRNwXETd0fp7q7vdZlJf3Vq5RjAB4JCJ+HBGTImL14rX6aET8qNj3u4horeUTryYTpu6tAzxbdn8yy09OqioiJgA7AHcVRccUw1gu6hziQn08hwR+V/zBfrIoWyszny+2XwDWKrbrIV6AQ1j6Q2e9ntv+nsd6iLnTRyn1JnTasPggcGvnmzel2CaX1al2vP35vdfDuX078GJmPlFWVhfntcv7VSO/bqXBsgmlHuLNi9thwG7Al4CvAn8H3p6ZOwAnAt/spb2zgbMzcxuW/nuH0t/i5yn1+G4EvK0oPzczd87MrYFW4L1lx4zIzF2K404qvhg6EfhFZm6fmb8YwHNW5WwGnJ+ZWwCvAUcD3wMOzMw3AxcBp5TVX+r32U17JwE3Z+ZWwFXA+mX7NgXOK/ZNBz5Y6SdTr0yYGkxErAT8L/D5zHwN+D6wMbA98DylN916sVtm7gjsCxwdEe8o31l8w501iawbxbdm+wNXFkX1fG4Xq7fzuDwRcTywALi0KHoeWL/4IPCfwGURsUqt4is0xO+9i0NZOtGvi/PazfvVYo30upUG2dOZ+VBmLgIeAW4q/j4eojTkaixwZXGtylnAVr209xaW/D92WZd9d2fm5OKxHijaB9gjIu6KiIeAd3Z5jKuLn/eV1Vf9eDYzby+2f0apZ35r4MYoXdd6AkuP6ujt97kbpZESZOb1wLSyfU9n5gO9HD8kmTB17zlgvbL76xZlNRURLZQ+fFyamVcDZOaLmbmwePP7EUuGhtX8OWTmc8XPl4BfFrG9WNY1PA54qahe83gpJXb3Z+aLUN/nlv6fx5rHHBFHUvrW8sPFhwGK4W2vFNv3UboW6E1FbOVv8FWLdwC/95qe24gYAfwbsPhb3no4r929X9GAr1upCuaWbS8qu78IGAF8A/hD0fvzPmB0hR5rIaXhXKOB8yn1SGxD6X1vdDfHLCziUX3p+sXT68AjRe/f9pm5TWbuU7Z/RX6fy7x++nl8wzJh6t49wKYRsWHR63AIcG0tAyrGE18IPJqZZ5aVl1/n8wGgc7aUa4FDImJURGxIqRv17irGOyYiVu7cpnTR/8NFXJ0zXR0BXFMW70eiZFdgRtnQnWpZ6lv6ej23ZTH05zzeAOwTER3FELN9irKqiIh3A/8F7J+Zs8rK14iI4cX2RpTO5T+KmF+LiF2L1/5Hyp7jYMfa3997rd8v9gL+npmLh97U+rz29H5Fg71upToxliVfFBzZh/p3smSo1CF9qN+ZHE0teoUP7MMxr1Oa0EW1t35EvKXYPozS73+NzrKIaImI3noly91O6ZpYImIfShPuNL2myQz7IzMXRMQxlP5jHg5clJmP1DistwGHAw8VXaxQGtt8aERsT+kbhmeA/wDIzEci4grgb5SGQB2dmQurGO9awC9Ln5sYAVyWmddHxD3AFRHxMWASxR8lcB2lmbKeBGYBR1Ux1s6kbm+K81f4dj2c24i4nNJsSKtHxGRK44u/RT/OY2a+GhHfoPThHuDrmdnXyQ4qEe9XKM0ud2PxmrgzMz9Faea3r0fEfErfpn6qLK7PABdTGk//W5a+7mkwY929v7/3arxfdBdrZl7IstfdQY3PKz2/X9Xt61aqY98GLomIE4Df9KH+54GfFUOgrwdmLK9yZk6PiB9R+nLoBZb8vS3PH4Djir/vU72OqaYeo3TZw0WU/n/6HqX/j86JiLGUPoN9l9Jwz774GnB5RBwO/JnSa+J1YKVKB95IohgZI0mSpAYXEW3A7MzMiDgEODQzazrTrxpHRIwCFhadB28Bvp+Z29c6rlqzh0mSJGnoeDNwbjE0djqlmUmlvlqf0kiAYcA84BM1jqcu2MOk/9/evcfaVZRhHP69UEAp3gVEUUAMiQUqKkGsCAZQK0QFDJJolTuoQEQwEg00kAhNtWi90UACHJGKphWjQqQRDFRuFgnaGiqVCgiUFhALtKeUUl7/mNl1s9z7nB56TvexfZ9kpesya2Z2m+6Vb+bbsyIiIiIiooss+hAREREREdFFAqaIiIiIiIguEjBFRERERER0kYApYjMgyZLW590aERER0SOSDpZ0vqQDe92X+K8ETBEjSFJfDVaa25297ttQSNpC0jOS9qjHi/JlHhERMXwk7QzMpryke7akHXvcpagSMEWMvBuBnRrbYT3t0dDtBTxne1H9At+F9Xu5YURExAaTdJ2kvrbjPknXbWCdG1zHcKnLwPcBM21PAH4FXNHTTsU6CZgiRt5q20sb21Oti3XG6XRJ10vql/SQpEntFUjaW9KNklZJeqp+yb+mUeZYSQskrZa0TNKPG/14vaRZklZK+kezjUFMAG6v+wcA99heNYT7IyJiE9PIolhTny3TJI3dCM1/GVjv55ikmyX9cEPqGEkuDrV9Rj0+2fbhve5XFAmYIkaHC4BfA/sAlwFXSdoXoD545gArgP2AIykBzLqRJ0mnApcCVwLjKTNYf220MZkyYvUu4OfAFZLeNlCnJC2XtByYDnys7l8NvKdeGxUjcxER0TOtLIq3A+cCXwKmdSooaevhatT207aX97qO2DwkYIoYeRMlrWhsUxtlrrV9qe1Fti8Efg+cWa99BhgLfM72Atu3AKcAR0l6Ry1zHjDd9nds32f7btvfbrTxE9tX276/ln8BGOx3SPsA7wb6gU/W44eB4+r+SUP7q4iIiE1MK4viYds/BWYCR8C6WZ0ZddbpCeC2el6SviZpcc2cWNCe9SBp2zp7taJmTHyj2Wgzna7Webakv9dMi0ckTWmVBQ4CTmubEdu1Qx3bSJpe23xO0p2SDmi0e7OkSyRdJOlJSY/Xz7dFW5kD670rJD0taZ6kvTp8hmtqX/7YOD+3nv/Z0P4pYqQkYIoYeXMpwUX71gxm7uhwPK7uvxOYb/vZtuu3Ay8C4yTtALwFuGmQfsxv7dh+AXgC2GGgG2w/CGwP9Nu+gRJkvRn4he0HbS8dpM2IiNi8rAK2ajueBAj4IPD5eu6bwInAaZRn3RTgUkmtFLRpwIeBTwGHUAbuBhvgu4gyGDgF2BM4mjLAByX17g5KFkbrt8QPd6jjW8AxwAm1zQXAuMfLUQAABIxJREFUDZJ2apT7LOV5OAE4nTLAeQyApDGUbI5bKRkd76Nkaazt0N6P6p/7Sdqz3v8m4AP1fN8gnzk2kjG97kDEZqC/zuqMBA+h7JoO93YdNJH0W8oDbgwwRtIKYEtgG+BfkrC93RD7GxERmyhJ+1GyItoH8B6wfXZbmbHAWcBHbP+hVabee5qkWyjB1Am259R7jgceGaDd7YCvAGfabqWr308djLT9tKTnKc/jpW33tdcxFvgicJLt6+u5LwAHUwK7c9uavNf25Lq/SNLJlMDuGuDVwGuB39heXMv8rVO/bd8qaT4llf7E+vdyJOXZvAT4XbfPHBtXZpgiRof9OxwvrPsLgb0lvart+gTK/9+Fth8HHqV8WQ+nkyizYXcD59T9OZQRuNZMWUREbN5aaefPUQKUucAZbdfvbpQfB7yCMnOzLlWdEqzsXretacu8sL2CMtvTzTjKYN5gmRYD2Z0yM3ZbW7treWnGR8v8xvESasZGXdSpD5ijspjTWYP8Xrg1yzRJ0lZA652JV9f2YxTIDFPEyNumTrG3W2v7ibbjoyTdBdxM+bI8hDKNDyUf/ALKQhCTgddRFni4tm3m6kLgu5KWAdcD2wKH2L745Xba9qM1tWA8MMn2A5LGA1NHcMYsIiL+v8yl/K52DbDEdjObYWXjuDVY/3Hgn41rayizM6NNM5tjwIwN28dLmg5MBD4BXCjpiNaMWcNMykDk9pRUwIPq+b5h6HcMk8wwRYy8Q4HHGts9jTLnU3K151NG2Y63fReA7X7go5Rp/nmU3Og7KF+s1DIzKCkDJ1NWx7uBksO9ofYFltdgaWdgR+BPw1BvRERsGvpt32/7oQ7BUif3AquBXep97dtDwGJKQLIu86Kmy/3PogltFtY6B8q0eJ6SVt7N4lqm9fshJG0JvL/2eUhs/8X2VNsfogyGHtul3Eqg9RqQi2sf77K9sFP56I3MMEWMINvHUVaUG8xS2xMHqGcBg6Tc2b4cuLzLNXU4t+tgnbJ9J/DWuv8IZbW+iIiIl8X2s5KmAdNUfkQ0F9iOEiC9aPsySZcDU+vKeksor8XoGuzUOr8HTJG0utb5BuC9dUAR4EHK4gq7Ul7T8VSjjpWSZtR2nwQeoPwuakfgkvX9fJJ2A06lvCrkUcpy6+OBGQPcdgkljbH1jO1b3/Zi40jAFBEREREb03nAMuCrlEDiGeDPlNQ06vmxwC8pr7X4AYMP2H0d+Hete+da/1Vt16dRZnLuBV4J7NahjnPqn1dSUgPvASbafmz9Pxr9wB7ALOCNtR8zgebrRNaxfZ+kmygZKauBLCc+ysgeyiJbETHcJBk42vbsXvclIiIiNj5J36fMMs2y/ele9ydeKgFTREREREQPSDoFOBw4jJJ2uL/teb3tVTRl0YeIiIiIiN6YQFlJbxlwSoKl0SkzTBEREREREV1khikiIiIiIqKLBEwRERERERFdJGCKiIiIiIjoIgFTREREREREFwmYIiIiIiIiukjAFBERERER0UUCpoiIiIiIiC4SMEVERERERHSRgCkiIiIiIqKL/wCp2o0RKUCG7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))\n",
    "\n",
    "# plot loss function\n",
    "ax1.plot(range(len(loss_list)), loss_list)\n",
    "ax1.set_title('Training Loss', fontsize=14)\n",
    "ax1.set_xlabel('Epoch #', fontsize=14)\n",
    "ax1.set_ylabel('Loss', fontsize=14)\n",
    "\n",
    "# print confusion matrix\n",
    "pred_yt = np.round(sigmoid_activation(y_test_pred.detach().numpy()), 0)\n",
    "conf_mat = confusion_matrix(y_true=test_yt.numpy(), y_pred=pred_yt)\n",
    "ax2.matshow(conf_mat, cmap=plt.cm.Blues, alpha=0.3)\n",
    "for i in range(conf_mat.shape[0]):\n",
    "    for j in range(conf_mat.shape[1]):\n",
    "        ax2.text(x=j, y=i, s=conf_mat[i, j], va='center', ha='center', size='xx-large')\n",
    "\n",
    "ax2.set_title('Confusion Matrix', fontsize=14)\n",
    "ax2.set_xlabel('Predictions $\\mathbf{\\hat{y}}$', fontsize=14)\n",
    "ax2.set_ylabel('Actuals $\\mathbf{y}$', fontsize=14)\n",
    "ax2.set_yticks([0, 1])\n",
    "ax2.set_xticks([0, 1])\n",
    "ax2.set_yticklabels(class_labels)\n",
    "ax2.set_xticklabels(class_labels)\n",
    "ax2.tick_params(top=False, bottom=True, labeltop=False, labelbottom=True)\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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}