{
 "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: 23/06/2024\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)). It consists of $N=569$ test persons with 2 classes (malignant and benign) defined as $y_i \\in \\{0,1\\}$ and $J=30$ measured attributes per person $i$ given as a feature vector $\\mathbf{x}_i\\in\\mathbb{R}^{J\\times 1}$."
   ]
  },
  {
   "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\n",
    "bc = load_breast_cancer()\n",
    "\n",
    "# normalize\n",
    "x_min = np.min(bc.data, 0)\n",
    "x_scale = 3*np.std(bc.data-x_min, 0)\n",
    "bc_norm = (bc.data-x_min) / x_scale\n",
    "\n",
    "# split data into training and validation set\n",
    "train_X, val_X, train_y, val_y = train_test_split(bc_norm, bc.target, random_state=42)\n",
    "class_labels = bc.target_names\n",
    "\n",
    "# plot shapes\n",
    "print(train_X.shape, train_y.shape, val_X.shape, val_y.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gradient descent\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^\\intercal\\mathbf{w}$ where a vector $\\mathbf{w}=\\left[w^{(1)}, w^{(2)}, \\dots, w^{(J)}\\right]^\\intercal$ consists of weights $w^{(j)}$ for each feature $j$ and the vector $\\mathbf{x}_i=\\left[x_i^{(1)}, x_i^{(2)}, \\dots, x_i^{(J)}\\right]^\\intercal$ 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 vector of $[1, 1, \\dots, 1] \\in \\mathbb{R}^{N}$ 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 may be defined as a mean squared error\n",
    "\n",
    "$$\n",
    "L(\\mathbf{w})=\\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^\\intercal\\mathbf{w}-y_i\\right)^2\n",
    "$$\n",
    "\n",
    "where $\\left(\\mathbf{x}_i^\\intercal\\mathbf{w}-y_i\\right)^2$ is the **loss function**. The **training** refers to an optimization problem where weights $\\mathbf{w}$ are adjusted so that the objective is\n",
    "\n",
    "$$\n",
    "\\text{arg min}_{\\mathbf{w}} \\, L(\\mathbf{w})\n",
    "$$\n",
    "\n",
    "To achieve this, SGD inherits the **Gradient Descent** update method at iteration $k$ (known as **back-propagation**), which writes\n",
    "\n",
    "$$\n",
    "\\mathbf{w}_{k+1} = \\mathbf{w}_k - \\gamma \\, \\nabla_{\\mathbf{w}_k} \\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right)^2 \\, , \\, \\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}} \\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right)^2$ can be generally obtained by\n",
    "\n",
    "$$\n",
    "\\nabla_{\\mathbf{w}_k} \\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right)^2 = \\frac{\\partial \\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right)^2}{\\partial \\mathbf{w}_k}\n",
    "=\n",
    "\\begin{bmatrix} \n",
    "\\frac{\\partial}{\\partial \\mathbf{w}_k^{(1)}} \\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right)^2 \\\\\n",
    "\\frac{\\partial}{\\partial \\mathbf{w}_k^{(2)}} \\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right)^2 \\\\\n",
    "\\vdots \\\\\n",
    "\\frac{\\partial}{\\partial \\mathbf{w}_k^{(J)}} \\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right)^2 \\\\\n",
    "\\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix} \n",
    "\\mathbf{x}_i^{(1)} 2\\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right) \\\\\n",
    "\\mathbf{x}_i^{(2)} 2\\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right) \\\\\n",
    "\\vdots \\\\\n",
    "\\mathbf{x}_i^{(J)} 2\\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-y_i\\right) \\\\\n",
    "\\end{bmatrix}\n",
    "= 2\\mathbf{x}_i^\\intercal\\left(\\mathbf{x}_i^\\intercal\\mathbf{w}_k-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 can be 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 it as the batch size."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# learning rate is the step size as in classical gradient descent\n",
    "l_rate = 1e-3\n",
    "\n",
    "# epochs is the number of maximum iterations for minimization\n",
    "epochs = 700\n",
    "\n",
    "# batch size amounts to the number of samples used in one epoch iteration (up to hardware memory)\n",
    "b_size = 2**4\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, val_X = np.c_[np.ones(train_X.shape[0]), train_X], np.c_[np.ones(val_X.shape[0]), val_X]\n",
    "\n",
    "# initialize weight vector such it has the same number of columns as input features\n",
    "np.random.seed(1111)\n",
    "w = np.random.uniform(size=(train_X.shape[1],)) * 0.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",
    "        preds_y = np.dot(batch_X, w)\n",
    "        \n",
    "        # compare prediction and true values\n",
    "        diff = preds_y - batch_y\n",
    "        \n",
    "        # compute mean of squared loss for current batch\n",
    "        epoch_loss.extend(diff**2)\n",
    "\n",
    "        # compute the derivative\n",
    "        gradient = 2 * np.dot(batch_X.T, diff)\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.mean(epoch_loss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4QAAAFSCAYAAACqthEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwRElEQVR4nO3de5xddX3v/9dn78lMyAVISABLEq7xgjfAiDfqDUUsfYD2aMUez8FLy9GHnGq1/movP7Sc46laj70c0UqVgm2Vo3jLz6KICmpVNOEiCAiEFEkil5AQQq6Tmfn8/thrJns2k8uQmf3dM/v1fHQee63vWt+1P3svmfQ967u+KzITSZIkSVL3qZUuQJIkSZJUhoFQkiRJkrqUgVCSJEmSupSBUJIkSZK6lIFQkiRJkrqUgVCSJEmSupSBUJIkICIui4hvjLPPdRHxicmqSZKkyRY+h1CSNJVExL7+4bo8M9/8BI57CI1/FzeNo898YFdmPjbe9xuPiPgg8LrMfMZkvo8kqfv0lC5AkqRxelLT8m8D/9jStr1554iYkZm79nXQzHx0vIVk5sbx9pEkqZM4ZFSSNKVk5gPDP8Cm5jZgJrApIt4YEd+LiO3Af4uIwyLiCxGxNiK2R8RtEfGW5uO2DhmthoN+MiL+V0Q8HBEPRcTHIqLWss8nmtbvjYi/iIhPR8Tm6v3e1/I+T46I70fEjoi4MyJ+KyK2RMSbn+h3EhHPjIjvVJ9tY/VZDmnZ/t2qpi0R8fOIeFm1bUZE/H1E/DoidkbEmoj48BOtRZI0tRgIJUnT0V8BnwROBL5GIyjeSOOK4tOBvwM+HRGn7+M4/xkYAF4IXAC8G3jDPvr8EXArcArwEeCjEfECgCpMfrU65vOBNwMfAPrG8dlGiYjZwNXAFuBU4LVVvZc27fZ54P5q+0nAB4Ed1bY/rPqcCyytPt+dT7QeSdLU4pBRSdJ09H8y88qWtr9uWr4kIl4OvBH47l6Oc3tmXlgt3xURfwCcDnxhL32+nZnDVw3/T0T8YdXnJ8ArgacAZ2TmOoCI+CPgR/vzofbg94DZwH8ZvpcxIs4Hro2IEzJzFXA08LHM/GXVZ1VT/6OBu4AfZmNigfuAHx9APZKkKcQrhJKk6Whl80pE1CPizyPilojYEBFbgN8BluzjOLe0rP8aOPwA+jwV+PVwGKysAIb2ccy9eRpwS8vENj+ujnlitf5x4DPVMNo/j4inNu17GY2rhndFxMURcVbzsFhJ0vTmL3xJ0nS0tWX9j4H30rhKeDqNAPQ1oHcfx2mdjCbZ97+dT6TPZEmAzPwgu4fPvhC4JSLeWm27ETgG+FMadV4OXGMolKTu4C97SVI3OA34/zLznzPzZuAe4MkF6vgl8BsR8RtNbcs4sH+P7wCeGRFzm9peWB3zjuGGzLw7M/8+M88CPgv8ftO2xzLzysx8B3AW8HLghAOoSZI0RXgPoSSpG9wFvCEiTgMeBv47cCxwU5vruIbGhC2XR8QfAwfRGM45QHU1by9mRsRJLW3bgH8F/hL4XERcCMwDPg18JTNXRcRBwMeALwH3AkfQCMg/BYiI99CYcOZmGlc3fw/YDKw9gM8pSZoiDISSpG7wP2kEwG/SeE7hZTSC1Il76TPhMnMoIl4LfAb4GY2A9l7gK+ye9XNPjufxAfaGzFwWEa8C/rY65g7g68C7qn0GaYTEy2g8r3ED8A0aw2gBHgPeR2OG0aze49WZue2JfEZJ0tQSjQnFJElSCRHxbBpX55Zl5g2Fy5EkdRkDoSRJbVRdIdwK3E1jMpePAwGcnP6jLElqM4eMSpLUXnNpPLB+MfAIcB3wR4ZBSVIJXiGUJEmSpC7lYyckSZIkqUsZCCVJkiSpS037ewgXLFiQxxxzTOkyJEmSJKmIG2644eHMXDjWtmkfCI855hhWrlxZugxJkiRJKiIifrWnbQ4ZlSRJkqQuZSCUJEmSpC5lIJQkSZKkLmUglCRJkqQuZSCUJEmSpC5lIJQkSZKkLmUglCRJkqQuZSCUJEmSpC5lIJQkSZKkLmUgLODBzTv4/E/v46HNO0qXIkmSJKmLGQgLuGf9Fv7sq7ey+uGtpUuRJEmS1MUMhAXUIgAYGsrClUiSJEnqZgbCAuq1KhCaByVJkiQVZCAsoMqDDKaJUJIkSVI5BsICHDIqSZIkqRMYCAvYPWTUQChJkiSpHANhAcNXCAe9QihJkiSpIANhASNDRr1CKEmSJKkgA2EBzjIqSZIkqRMYCAsYmWXURChJkiSpIANhATUnlZEkSZLUAYoEwog4MyLujIhVEfH+Mba/PSJujYibI+LfI+LEqv2YiNhetd8cEf/Q/uoPXN17CCVJkiR1gJ52v2FE1IGLgVcCa4EVEbE8M29v2u3zmfkP1f5nAx8Hzqy23ZOZJ7Wx5Am3e5bRwoVIkiRJ6molrhCeCqzKzNWZ2Q9cAZzTvENmbm5anQ1Mq0tptepb9wqhJEmSpJJKBMKjgDVN62urtlEi4p0RcQ/wUeAPmzYdGxE3RcT3I+I3J7fUyTHy2AknlZEkSZJUUMdOKpOZF2fm8cCfAH9RNd8PLMnMk4H3AJ+PiINb+0bE+RGxMiJWrl+/vn1F76fhx04MeoVQkiRJUkElAuE6YHHT+qKqbU+uAF4DkJk7M3NDtXwDcA/w5NYOmXlJZi7LzGULFy6cqLonzO4H0xcuRJIkSVJXKxEIVwBLI+LYiOgFzgWWN+8QEUubVs8C7q7aF1aT0hARxwFLgdVtqXoCDT+H0CGjkiRJkkpq+yyjmTkQERcAVwN14NLMvC0iLgJWZuZy4IKIeAWwC3gEOK/q/mLgoojYBQwBb8/Mje3+DAdqZMiogVCSJElSQW0PhACZeRVwVUvbhU3L79pDvy8DX57c6iafD6aXJEmS1Ak6dlKZ6azmg+klSZIkdQADYQF1H0wvSZIkqQMYCAvwwfSSJEmSOoGBsAAfTC9JkiSpExgIC6j7HEJJkiRJHcBAWECVBxl0yKgkSZKkggyEBUQEtXDIqCRJkqSyDISF1GvhpDKSJEmSijIQFhIRDhmVJEmSVJSBsJB6hENGJUmSJBVlICykMWS0dBWSJEmSupmBsJAIGDQRSpIkSSrIQFiIk8pIkiRJKs1AWEg9DISSJEmSyjIQFhIRDA6VrkKSJElSNzMQFlKvQXqFUJIkSVJBBsJCahFOKiNJkiSpKANhITUfTC9JkiSpMANhIfVaYB6UJEmSVJKBsJCazyGUJEmSVJiBsJBazSGjkiRJksoyEBZSj3CWUUmSJElFGQgLcZZRSZIkSaUZCAup1XwwvSRJkqSyDISF+GB6SZIkSaUZCAvxOYSSJEmSSjMQFlKLwFsIJUmSJJVkICykFjBkIpQkSZJUkIGwkHrNWUYlSZIklWUgLKQxZNRAKEmSJKmcIoEwIs6MiDsjYlVEvH+M7W+PiFsj4uaI+PeIOLFp259W/e6MiFe1t/KJYyCUJEmSVFrbA2FE1IGLgVcDJwJvbA58lc9n5jMz8yTgo8DHq74nAucCTwfOBD5ZHW/KccioJEmSpNJKXCE8FViVmaszsx+4AjineYfM3Ny0OhsYTk7nAFdk5s7M/A9gVXW8KadWc5ZRSZIkSWX1FHjPo4A1Tetrgee17hQR7wTeA/QCL2/qe31L36Mmp8zJVQscMipJkiSpqI6dVCYzL87M44E/Af5iPH0j4vyIWBkRK9evXz85BR6gejhkVJIkSVJZJQLhOmBx0/qiqm1PrgBeM56+mXlJZi7LzGULFy48sGoniUNGJUmSJJVWIhCuAJZGxLER0UtjkpjlzTtExNKm1bOAu6vl5cC5EdEXEccCS4GftaHmCeeD6SVJkiSV1vZ7CDNzICIuAK4G6sClmXlbRFwErMzM5cAFEfEKYBfwCHBe1fe2iPgicDswALwzMwfb/RkmQr3mYyckSZIklVViUhky8yrgqpa2C5uW37WXvh8CPjR51bVHRDBoIJQkSZJUUMdOKjPd1SMcMipJkiSpKANhIXUnlZEkSZJUmIGwkAh87IQkSZKkogyEhdTDSWUkSZIklWUgLMRZRiVJkiSVZiAsJCIYHCpdhSRJkqRuZiAspF7DK4SSJEmSijIQFuI9hJIkSZJKMxAW0hgyaiCUJEmSVI6BsJB6LfACoSRJkqSSDISF1HwOoSRJkqTCDISF1GrBoJcIJUmSJBVkICykHkEaCCVJkiQVZCAspOakMpIkSZIKMxAWUqsFQ4lXCSVJkiQVYyAspB4B4EyjkiRJkooxEBZSa+RBJ5aRJEmSVIyBsJBalQi9j1CSJElSKQbCQuo1h4xKkiRJKstAWIhDRiVJkiSVZiAspFZNKjNkIJQkSZJUiIGwkJFA6D2EkiRJkgoxEBZSd1IZSZIkSYUZCAsZnmXUPChJkiSpFANhIcOTyngPoSRJkqRSDISF1MMho5IkSZLKMhAWsnvIqIFQkiRJUhkGwkJ2zzJauBBJkiRJXctAWEi9+uZ9ML0kSZKkUgyEhfhgekmSJEmlGQgL8cH0kiRJkkorEggj4syIuDMiVkXE+8fY/p6IuD0ibomI70bE0U3bBiPi5upneXsrnzh1n0MoSZIkqbCedr9hRNSBi4FXAmuBFRGxPDNvb9rtJmBZZm6LiHcAHwXeUG3bnpkntbPmyTD8HEIfOyFJkiSplBJXCE8FVmXm6szsB64AzmneITOvzcxt1er1wKI21zjpvIdQkiRJUmklAuFRwJqm9bVV2568Dfhm0/rMiFgZEddHxGsmob62qPscQkmSJEmFtX3I6HhExJuAZcBLmpqPzsx1EXEc8L2IuDUz72npdz5wPsCSJUvaVu94DF8hdMioJEmSpFJKXCFcByxuWl9UtY0SEa8A/hw4OzN3Drdn5rrqdTVwHXBya9/MvCQzl2XmsoULF05s9ROk5hVCSZIkSYWVCIQrgKURcWxE9ALnAqNmC42Ik4FP0wiDDzW1z4uIvmp5AfAioHkymimjHs4yKkmSJKmstg8ZzcyBiLgAuBqoA5dm5m0RcRGwMjOXA38NzAG+FI3gdF9mng08Dfh0RAzRCLMfbpmddMpwllFJkiRJpRW5hzAzrwKuamm7sGn5FXvo92PgmZNbXXuMDBk1EEqSJEkqpMiD6eWD6SVJkiSVZyAsZGTIqJPKSJIkSSrEQFiID6aXJEmSVJqBsJCRQOiYUUmSJEmFGAgLGb6H0FlGJUmSJJViICyk5nMIJUmSJBVmICykVn3z3kMoSZIkqRQDYSEz6o2vftfgUOFKJEmSJHUrA2EhvVUg3DlgIJQkSZJUhoGwkN4erxBKkiRJKqtnvB0iog/4DeAgYH1mrp/wqrrA8BXCfq8QSpIkSSpkv64QRsTciHhHRPwAeBRYBfwCeCAi7ouIf4yI505modPN8BVCA6EkSZKkUvYZCCPiPcC9wFuBa4BzgJOAJwMvAD5I40rjNRHxrYhYOkm1TisGQkmSJEml7c+Q0ecDL8nMX+xh+8+ASyPi7cDbgJcAd09QfdNWTy2IgH7vIZQkSZJUyD4DYWb+7v4cKDN3Ap884Iq6RETQW68ZCCVJkiQVM65ZRiPitRExf7KK6Ta99ZpDRiVJkiQVM95ZRr8MDEXEL4DrgGuB72fmpgmuqyv09hgIJUmSJJUz3ucQDlV9ngX8d+ArwMMRcWNE/O+JLm66MxBKkiRJKmm8gfBQ4EzgQ8APgZ3VMU4C3j2BdXWF3h7vIZQkSZJUzriGjGbmFuDbEXEjcANwOnAeMGcSapv2vIdQkiRJUknjCoQR8RngNKD5WYO3Af9O44qhxqG3p8YurxBKkiRJKmS8k8q8FUhgM/Ap4O8y88EJr6pLzKjX2OkVQkmSJEmFjPcewiuBB4BDgPcD/xERP4iIv4qIsya8umnOSWUkSZIklTTeewh/FyAijqMxdPQ04PXAi4D3jfd43a6vp8aWnQOly5AkSZLUpcZ7D2EPcAqNIPii6udgICa+tOnPSWUkSZIklTTeK3qPAjOr5eYQuJrGQ+o1Dk4qI0mSJKmk8QbCg6rXe2kEwOuAazNz7QTW1DW8h1CSJElSSeMNhG+hEQDvm4xius0Mh4xKkiRJKmi8k8pcPlmFdKPenhr9DhmVJEmSVMg+HzsREcfu78GiYfGBldQ9en0OoSRJkqSC9uc5hD+JiM9GxAv2tENEzIuIdwC3A+fs64ARcWZE3BkRqyLi/WNsf09E3B4Rt0TEdyPi6KZt50XE3dXPeftRf8fqc1IZSZIkSQXtz5DRpwJ/DvxbRAwBNwC/BnYA84ATgacBPwPenZlX7+1gEVEHLgZeCawFVkTE8sy8vWm3m4BlmbmtCpofBd4QEfOBDwDLgARuqPo+st+fuIM4qYwkSZKkkvZ5hTAzN2Xm+4CjgLcDdwCHAscCA8DlwMmZ+aJ9hcHKqcCqzFydmf3AFbRcVczMazNzW7V6PbCoWn4VcE1mbqxC4DXAmfvxnh2pt15jKGHAq4SSJEmSCtjvSWUycztwZfVzII4C1jStrwWet5f93wZ8cy99jzrAeoqZ0dPI4/2DQ/TU92f0riRJkiRNnAlJIRHx1Ij4+UQcq+W4b6IxPPSvx9nv/IhYGREr169fP9FlTZjeKgQ6bFSSJElSCRN1WaoOPGM/910HNM9EuqhqGyUiXkHj3sWzM3PnePpm5iWZuSwzly1cuHA/y2q/3h4DoSRJkqRySoxTXAEsjYhjI6IXOBdY3rxDRJwMfJpGGHyoadPVwBnVrKbzgDOqtimpt2nIqCRJkiS1237dQxgRl9GYXfQG4OamCV/GLTMHIuICGkGuDlyambdFxEXAysxcTmOI6BzgSxEBcF9mnp2ZGyPif9AIlQAXZebGJ1pLaX1eIZQkSZJU0P5OKrMYOJvG7KKDEXEXuwPijTRmG91vmXkVcFVL24VNy6/YS99LgUvH836dauQeQq8QSpIkSSpgvwJhZp4OEBHHAc+pfk4B/l9g/vBuk1HgdDbDSWUkSZIkFbTfj50AyMzVwGrgS8NtEXEMjZlAT5nQyrqAk8pIkiRJKmlcgXAsmXkvcC8H/nzCruOkMpIkSZJK8mnoBXmFUJIkSVJJBsKCfDC9JEmSpJIMhAX1OWRUkiRJUkEGwoKcZVSSJElSSQbCgryHUJIkSVJJBsKCHDIqSZIkqSQDYUGz+xpP/di6c7BwJZIkSZK6kYGwoL6eGj21YMvOXaVLkSRJktSFDIQFRQSz+3rYsmOgdCmSJEmSupCBsLA5fT1sccioJEmSpAIMhIXNndnjkFFJkiRJRRgIC5vd18OWnQ4ZlSRJktR+BsLCHDIqSZIkqRQDYWFz+nrYssMho5IkSZLaz0BY2ByHjEqSJEkqxEBY2JyZPT6YXpIkSVIRBsLChieVGRrK0qVIkiRJ6jIGwsLm9vUAsLXfYaOSJEmS2stAWNicmVUgdNioJEmSpDYzEBY2u7pC6MPpJUmSJLWbgbCw4SGjj+1wyKgkSZKk9jIQFuaQUUmSJEmlGAgLm93rkFFJkiRJZRgIC5s7czgQeoVQkiRJUnsZCAsbmVRmh1cIJUmSJLWXgbCw2X11ALbsdFIZSZIkSe1lICysr6fOrN46j2zzCqEkSZKk9jIQdoAFc/p4eMvO0mVIkiRJ6jJFAmFEnBkRd0bEqoh4/xjbXxwRN0bEQES8rmXbYETcXP0sb1/Vk+ewOb0GQkmSJElt19PuN4yIOnAx8EpgLbAiIpZn5u1Nu90HvBn44zEOsT0zT5rsOttpwZw+7tuwrXQZkiRJkrpMiSuEpwKrMnN1ZvYDVwDnNO+Qmfdm5i3AUIH62s4ho5IkSZJKKBEIjwLWNK2vrdr218yIWBkR10fEaya0skIWzull47Z+Bga7Iv9KkiRJ6hBtHzI6AY7OzHURcRzwvYi4NTPvad4hIs4HzgdYsmRJiRrHZcHcPjJh47Z+Dp87s3Q5kiRJkrpEiSuE64DFTeuLqrb9kpnrqtfVwHXAyWPsc0lmLsvMZQsXLjywattgwZw+AB5+rL9wJZIkSZK6SYlAuAJYGhHHRkQvcC6wX7OFRsS8iOirlhcALwJu33uvzjccCDds9T5CSZIkSe3T9kCYmQPABcDVwB3AFzPztoi4KCLOBoiI50bEWuD1wKcj4raq+9OAlRHxc+Ba4MMts5NOSQvm9AI4sYwkSZKktipyD2FmXgVc1dJ2YdPyChpDSVv7/Rh45qQX2GYL5jpkVJIkSVL7FXkwvUab29dDb0/NK4SSJEmS2spA2AEigiMO7uP+R3eULkWSJElSFzEQdojF82ax5pFtpcuQJEmS1EUMhB1iyfxZrNloIJQkSZLUPgbCDrF4/iwe3tLPtv6B0qVIkiRJ6hIGwg6xeP4sANZs3F64EkmSJEndwkDYIRbPOwjAYaOSJEmS2sZA2CGGrxDeZyCUJEmS1CYGwg5x2OxeZvXWnWlUkiRJUtsYCDtERLBk/izufXhr6VIkSZIkdQkDYQdZesRc7npwS+kyJEmSJHUJA2EHeeqRc1m3aTuP7dhVuhRJkiRJXcBA2EGecsRcAO568LHClUiSJEnqBgbCDvKUIxuB8JcPGAglSZIkTT4DYQdZNO8gZvfWuctAKEmSJKkNDIQdJCJ4ypFzueN+A6EkSZKkyWcg7DDPWnQot657lIHBodKlSJIkSZrmDIQd5uQlh7J91yB3OrGMJEmSpElmIOwwpyyZB8BN920qW4gkSZKkac9A2GEWzTuIBXN6DYSSJEmSJp2BsMNEBCctnscNv9pYuhRJkiRJ05yBsAO94PjDuHfDNtZt2l66FEmSJEnTmIGwA512wgIAfrTq4cKVSJIkSZrODIQd6MlHzGHBnD4DoSRJkqRJZSDsQBHBaSccxr/f/TBDQ1m6HEmSJEnTlIGwQ73sqYezYWs/N63ZVLoUSZIkSdOUgbBDveyphzOjHnz79gdKlyJJkiRpmjIQdqiDZ87g+ccdxjW3PUimw0YlSZIkTTwDYQd71dOPZPXDW7nj/sdKlyJJkiRpGjIQdrCznvkkemrBV29aW7oUSZIkSdOQgbCDzZvdy0ufcjhfv/nXDDrbqCRJkqQJViQQRsSZEXFnRKyKiPePsf3FEXFjRAxExOtatp0XEXdXP+e1r+oy/tMpR/HQYzu57s6HSpciSZIkaZppeyCMiDpwMfBq4ETgjRFxYstu9wFvBj7f0nc+8AHgecCpwAciYt5k11zSK048gsPn9vHP1/+qdCmSJEmSppkSVwhPBVZl5urM7AeuAM5p3iEz783MW4Chlr6vAq7JzI2Z+QhwDXBmO4ouZUa9xrmnLuH7d63n3oe3li5HkiRJ0jRSIhAeBaxpWl9btU1Y34g4PyJWRsTK9evXP+FCO8Wbnr+E3nqNT163qnQpkiRJkqaRaTmpTGZekpnLMnPZwoULS5dzwA6fO5M3nrqEr9y4jjUbt5UuR5IkSdI0USIQrgMWN60vqtomu++U9vaXHE8tgk9ed0/pUiRJkiRNEyUC4QpgaUQcGxG9wLnA8v3sezVwRkTMqyaTOaNqm/aOPGQmb3juYq68YQ1rH/EqoSRJkqQD1/ZAmJkDwAU0gtwdwBcz87aIuCgizgaIiOdGxFrg9cCnI+K2qu9G4H/QCJUrgIuqtq7wjpc2rhL+1VW/LF2KJEmSpGkgMqf3A8+XLVuWK1euLF3GhPnE9+7mY9++i396y3N52VMOL12OJEmSpA4XETdk5rKxtk3LSWWms/NffDwnHD6HC7/+C7b3D5YuR5IkSdIUZiCcYnp7avzP1zyDNRu389dX31m6HEmSJElTmIFwCnr+cYdx3guO5tIf/Qff+sUDpcuRJEmSNEUZCKeoPzvraTx70SG870s/51cbtpYuR5IkSdIUZCCcovp66nzi906hVgt+//KVPLK1v3RJkiRJkqYYA+EUtnj+LP7hTc/hVxu38ZbLVrB150DpkiRJkiRNIQbCKe4Fxx/GJ954Mreue5Tfv3wlWwyFkiRJkvaTgXAaOOPpR/Kx1z+Ln927kXMv+QnrH9tZuiRJkiRJU4CBcJp47cmL+Mf/+hxWPbSF3/nUj7jt14+WLkmSJElShzMQTiMvf+oRfOEPnk//wBCv/eSP+fxP7yMzS5clSZIkqUMZCKeZk5fM49/+8Dd53rHz+bOv3spbL1vBuk3bS5clSZIkqQMZCKehBXP6uPwtp3Lhb5/I9as3csbHv8+nrruHHbsGS5cmSZIkqYMYCKepWi1462nH8u0/ejEvOP4wPvKtX/Lyj13H5396n8FQkiRJEgAx3e8xW7ZsWa5cubJ0GcX95J4NfPhbv+TnazaxYE4fb3nRMbzpeUdzyKwZpUuTJEmSNIki4obMXDbmNgNh98hMfnLPBv7hB6v5wV3r6eup8aqnH8nrnrOIF52wgHotSpcoSZIkaYLtLRD2tLsYlRMRvPCEBbzwhAXc/uvNXLHiPr5+869Z/vNfc+TBM3nV04/glSceyfOOm8+MuqOJJUmSpOnOK4RdbseuQb57x0N87eZ1/PDu9ezYNcTBM3t4wfGH8cLjF/DC4w/jhMPnEOHVQ0mSJGkq8gqh9mjmjDpnPetJnPWsJ7G9f5Af3r2e79zxID9atYGrb3sQaMxa+txj5vHMRYfwrKMO5ZlHHeK9h5IkSdI0YCDUiIN665zx9CM54+lHArBm4zZ+cs8GfnzPw9y0ZhPf/MUDI/sec9gsnnHUITz5iLmccPgclh4+h6MPm01vj0NNJUmSpKnCQKg9Wjx/Fovnz+J3n7sYgE3b+vnFus38fO0mbl37KDev2cQ3brl/ZP96LTj6sFkce9jskb5L5s9i8fyDWDxvFrP7/J+bJEmS1En8/9C13w6d1ctpSxdw2tIFI23b+gdYvX4rqx7awt0PPcaqh7bwqw3buH71Brb2j37e4fzZvRxx8EyOOLiPI+Y2Xg8/eObutoNnctjsXnqc0EaSJElqCwOhDsis3h6ecdQhPOOoQ0a1ZyaPbNvFmo3buG/jNtY8so21j2znoc07eXDzDm7/9WYe3rKToTHmNDrkoBnMn93LvFkzmDerl3mze6v1XubPnlG99nLIQTOYO3MGc2f2MKu37sQ3kiRJ0jgZCDUpIoL5VZB79uJDx9xnYHCIDVv7eXDzDh7cvJMHNu/g4cd2smlbPxu37eKRrf3c/+gObr9/Mxu39rNzYGiP71evBXNn9jR++mZUyzM4+KAeDp45Y2TbnL4ZzOqtM6u3zuy+Hg7qrTO7t2dUW19PzXApSZKkrmAgVDE99Vo1XHTmPvfNTLbvGuSRKihu2NrP5u272LxjF4/tGOCxkdcBNm9vLK99ZBuP3T/A5h272LJzgP19wkotGlc+G2GxzqzhwNjXw6wZdWbOqNHX03idOaNO34w6fT2N5ZkzaszsqTfam9ua+vUNr/fUmVEPw6ckSZKKMRBqSoiIKpj1cNShB427/9BQsrV/gK07B9naP8D2/kG27hxg265Btu0cZFv/ANv6B6ufgVGvW6vtm7fv4sFHd7BjYJAduwbZOTDEjl2D7Ni15yuX+1ILRkJib73GjHqNvp4avT2N5d6eqr167eupMaMejfamffrqTfs3tw8v1x9/zJ56MKMe9NSGl2v01IKeem2k3cAqSZI0vRkI1RVqtajuN5z45ydmJjsHhti5a4idA42AOFZobLzubtu9bZD+gSH6B4fYOTDErsGkf6DR1lgeYtv2Xeyq9mm0N16H+/UPDu33FdDxqteCnloVGKvgOKMKjj31YEYVKHtG2scOl8377ml7Ty1G3q9eC+q1xn612uhto9dr1GrQU6u19B29PLxvfQ/be2qGX0mS1H0MhNIBiohqSGgdmPjAuT8yk4GhHBUWdw4vN4XInQNDI0FzYHCIXUON14HBZNdQ9To4xEDVvqtpfVe138DQnvrvXt4yMDDmsQYe9x6N952sMDtetWDMgLnnQLo7sNYC6tG0XAtqESOvjeVGe0RQH7Wt2r/WaK8FI8sj+9eajx9Nx6fp+GPVsnuf0bU02vddS1Cr6q5Ho5Za7H7PqGodbgt4/D4jNbK7T+x+P4O4JEnlGAilaSCiMfxzxhR9ZMdgFTiHqmA7ONh4Hb2+e/vAYDI4lAxm43X0+tDY24eXh5LBwSEGk8a+1fFHbW/6GRgaalre/To0lCNhdyhhKHf3GRqCXYONfkOZ1bbG0OWhbLzX0NDw6+6+jX2pjlFtr/YZrm+6Gg6Gza+jQueoQNm8vbH/4/cZo39tD/3Zj332UFPra3ONjbbdxx8Oy1G1NY47uu1xfUZCdnPg3sNx9tRnb+/d0qdWG/s4zX2g6Xtv6kNLjcPLo9tHH6f1c9WqPw6Mrn0v713tA40/TLTWyxjff1XSyH6MfBdU7z36OMPHH3ltPZ/N2/zjhqQpyEAoqbjG1ap66TKmhMzdYXckTA4HzKF8XDjN3B0mM3cH3332zWr/IUbC6uBQklUNQ03vn1W/4f6Zjf2GhprbHr/PUFPbyDETsmV9dP+m/Yce3z9h9HuMsc+ejtn4w0TLe+yhxtb+WdXd/AcA2L1P43trLNP4v5HjJ7v75xh9NPXsMTAOB1DGCJxN+9Pcf4ztzYF697axw+9IPXsJx4xRb1XJqLAbTcca/jwjYXpU3zH6NH12aKmn5bt4fN/Hf5eM+iyPD++jv+vW89Dy+ZsWmj93c62t23a/trSN0efx20afmzGP3fw9TUStTSv7V+MTrLVl2+6+u997T9/HAde6+4t5/LEf97nGW+vja2z+32LrtheesIA5fVMnZk2dSiVJRDTu0/SXd/fI4cDJ7iDaGiKHA+twiBxXn72E0VHHpBGw9xVgm/vsT72M2m/3MjT9gaCpD82fqem9yceH70b76HA9+jPtPiaM/kzZsj7W+RjeTtPxyByz76j9W957T+83sn0PtQ0fufV8jPq8Le83etvo89B8nob7jj5e6+fZz++ypTaAHIJkaFTfHNV3L9/lmN9f63fZ/Lke/12OfDfN6y3neXh/9tCndRtjbtvDcZo+L0nLPvtRa1PfPW0b6zhqj+++9yXMWTindBn7rcj/TxERZwJ/B9SBz2Tmh1u29wGfA54DbADekJn3RsQxwB3AndWu12fm29tWuCRJbdZ8lahO7H1nSdoPrX+MGNU2sl69tvxRpHXb/h5nvOE1W5Ly/gTk4T/YHHCtLX/8aD1O6x8RWut5IjPil9T2QBgRdeBi4JXAWmBFRCzPzNubdnsb8EhmnhAR5wIfAd5QbbsnM09qZ82SJEnSdNE6TLNqLVKLyisxA8WpwKrMXJ2Z/cAVwDkt+5wDXF4tXwmcHt6pLUmSJEkTqkQgPApY07S+tmobc5/MHAAeBQ6rth0bETdFxPcj4jcnu1hJkiRJmq6m2rwE9wNLMnNDRDwH+FpEPD0zNzfvFBHnA+cDLFmypECZkiRJktT5SlwhXAcsblpfVLWNuU9E9ACHABsyc2dmbgDIzBuAe4Ant75BZl6Smcsyc9nChQsn4SNIkiRJ0tRXIhCuAJZGxLER0QucCyxv2Wc5cF61/Drge5mZEbGwmpSGiDgOWAqsblPdkiRJkjSttH3IaGYORMQFwNU0HjtxaWbeFhEXASszcznwWeCfI2IVsJFGaAR4MXBRROwChoC3Z+bGdn8GSZIkSZoOIqf5kyqXLVuWK1euLF2GJEmSJBURETdk5rKxtpUYMipJkiRJ6gAGQkmSJEnqUgZCSZIkSepS0/4ewohYD/yqdB1jWAA8XLoIjYvnbOrxnE09nrOpxfM19XjOph7P2dTSqefr6Mwc83l80z4QdqqIWLmnGzvVmTxnU4/nbOrxnE0tnq+px3M29XjOppapeL4cMipJkiRJXcpAKEmSJEldykBYziWlC9C4ec6mHs/Z1OM5m1o8X1OP52zq8ZxNLVPufHkPoSRJkiR1Ka8QSpIkSVKXMhAWEBFnRsSdEbEqIt5fuh41RMSlEfFQRPyiqW1+RFwTEXdXr/Oq9oiIv6/O4S0RcUq5yrtTRCyOiGsj4vaIuC0i3lW1e846VETMjIifRcTPq3P2l1X7sRHx0+rc/N+I6K3a+6r1VdX2Y4p+gC4VEfWIuCkivlGte746WETcGxG3RsTNEbGyavP3YgeLiEMj4sqI+GVE3BERL/Ccda6IeEr139fwz+aIePdUPmcGwjaLiDpwMfBq4ETgjRFxYtmqVLkMOLOl7f3AdzNzKfDdah0a529p9XM+8Kk21ajdBoD3ZuaJwPOBd1b/LXnOOtdO4OWZ+WzgJODMiHg+8BHgbzLzBOAR4G3V/m8DHqna/6baT+33LuCOpnXPV+d7WWae1DT1vb8XO9vfAd/KzKcCz6bx35vnrENl5p3Vf18nAc8BtgFfZQqfMwNh+50KrMrM1ZnZD1wBnFO4JgGZ+QNgY0vzOcDl1fLlwGua2j+XDdcDh0bEk9pSqADIzPsz88Zq+TEa/4AeheesY1Xf/ZZqdUb1k8DLgSur9tZzNnwurwROj4hoT7UCiIhFwFnAZ6r1wPM1Ffl7sUNFxCHAi4HPAmRmf2ZuwnM2VZwO3JOZv2IKnzMDYfsdBaxpWl9btakzHZGZ91fLDwBHVMuexw5SDU07GfgpnrOOVg0/vBl4CLgGuAfYlJkD1S7N52XknFXbHwUOa2vB+lvg/wGGqvXD8Hx1ugS+HRE3RMT5VZu/FzvXscB64J+qodmfiYjZeM6minOBL1TLU/acGQil/ZSNKXmdlrfDRMQc4MvAuzNzc/M2z1nnyczBapjNIhojJp5atiLtSUT8NvBQZt5QuhaNy2mZeQqNYWrvjIgXN2/092LH6QFOAT6VmScDW9k91BDwnHWq6v7ps4EvtW6baufMQNh+64DFTeuLqjZ1pgeHL+tXrw9V7Z7HDhARM2iEwX/NzK9UzZ6zKaAaEnUt8AIaw2d6qk3N52XknFXbDwE2tLfSrvYi4OyIuJfG7Q0vp3Gvk+erg2Xmuur1IRr3NZ2Kvxc72VpgbWb+tFq/kkZA9Jx1vlcDN2bmg9X6lD1nBsL2WwEsrWZp66VxqXl54Zq0Z8uB86rl84CvN7X/12rmqOcDjzYNE1AbVPcmfRa4IzM/3rTJc9ahImJhRBxaLR8EvJLGvZ/XAq+rdms9Z8Pn8nXA99KH57ZNZv5pZi7KzGNo/Fv1vcz8z3i+OlZEzI6IucPLwBnAL/D3YsfKzAeANRHxlKrpdOB2PGdTwRvZPVwUpvA588H0BUTEb9G4L6MOXJqZHypbkQAi4gvAS4EFwIPAB4CvAV8ElgC/An43MzdWYeQTNGYl3Qa8JTNXFii7a0XEacAPgVvZfX/Tn9G4j9Bz1oEi4lk0brSv0/iD5Bcz86KIOI7GFaj5wE3AmzJzZ0TMBP6Zxv2hG4FzM3N1meq7W0S8FPjjzPxtz1fnqs7NV6vVHuDzmfmhiDgMfy92rIg4icbETb3AauAtVL8j8Zx1pOoPLvcBx2Xmo1XblP3vzEAoSZIkSV3KIaOSJEmS1KUMhJIkSZLUpQyEkiRJktSlDISSJEmS1KUMhJIkSZLUpQyEkiQVEhEZEa/b956SJE0OA6EkqetExGVVGGv9ub50beMREbWI2BwRT67W74qIF5euS5I0dfSULkCSpEK+A/yXlrb+EoUcgGcAOzLzrog4AjgaWFG4JknSFOIVQklSt9qZmQ+0/Gwc3lhdMbwgIv4tIrZFxK8i4k3NB4iIZ0bEdyJie0RsrK48HtKyz3kRcWtE7IyIByPi8pY65kfElyJia0Ssbn2PfXgh8ONq+TTgpszcPo7+kqQuZyCUJGnP/hJYDpwEXAJ8LiKWAUTEbOBqYAtwKvBaGgHt0uHOEfHfgE8D/wQ8C/gt4Bct73Eh8HXg2cD/BS6NiCV7KyoiNkXEJuBvgVdXy/8CnFJt+8YT/cCSpO4SmVm6BkmS2ioiLgPeBOxo2XRxZv5JtU8Cn8nMP2jq9x3ggcx8U0T8AfAxYFFmPlZtfylwLbA0M1dFxFrgXzLz/XuoI4EPZ+afVus9wGbg/Mz8l73UfwwQwA3A7wG/BL4NfJDGFcMdmfnA/n4fkqTu5T2EkqRu9QPg/Ja2TS3rPxlj/axq+WnALcNhsPJjYAg4MSI2A0cB391HHbcML2TmQESsBw7fW4fMvDciTgW2Zea3ImIR8BvAlzNz5z7eT5KkEQZCSVK32paZqybp2OMZfrNrjL57vKUjIr4J/CaNf8N7ImILUAf6gA0RQWbOGWe9kqQu5T2EkiTt2fPHWL+jWr4DeGZEzG3a/kIa/7bekZkPAeuA0ye4pt+ncU/jDcCfVMtXAx+tlk+a4PeTJE1jXiGUJHWrvog4sqVtMDPXN63/TkSsAK4DXkcj3D2v2vavNCad+VxEXAjMozGBzFearjx+CPibiHgQ+DdgFnB6Zv7vJ1p0Zq6r7jV8FvCmzPyPiHgW8JFJvOIpSZqmDISSpG71CuD+lrZ1wKKm9Q8C/wn4e2A98JbMXAGQmdsi4lU0Zvr8GY0Jar4OvGu4c2Z+KiL6gfcCHwE2AldNQO3LgE1VGFwEHAGsnIDjSpK6jLOMSpI0hmoG0Ndn5pWla5EkabJ4D6EkSZIkdSkDoSRJkiR1KYeMSpIkSVKX8gqhJEmSJHUpA6EkSZIkdSkDoSRJkiR1KQOhJEmSJHUpA6EkSZIkdSkDoSRJkiR1qf8f1ZTyQvv50wkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "fig = plt.figure(figsize=(15, 5))\n",
    "# plot loss function\n",
    "plt.plot(range(len(loss_list)), loss_list)\n",
    "plt.title('Training Loss', fontsize=14)\n",
    "plt.xlabel('Epoch #', fontsize=14)\n",
    "plt.ylabel('$L(\\mathbf{w})$', fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Validation\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",
    "\\sigma(\\hat{y}_i)=\\frac{1}{1+\\exp(-\\hat{y}_i)}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmoid = lambda y: 1.0 / (1 + np.exp(-y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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. Here, we estimate an ideal threshold via [Youden's method](https://en.wikipedia.org/wiki/Youden%27s_J_statistic)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# threshold estimation\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "fpr, tpr, thresholds = roc_curve(train_y, np.dot(train_X, w))\n",
    "tau = thresholds[np.argmax(tpr - fpr)]\n",
    "\n",
    "# compute predictions from test set\n",
    "pred_y = (sigmoid(np.dot(val_X, w)) >= tau).astype('uint8')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAFXCAYAAABp4B/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlQElEQVR4nO3dd5hU5dnH8e8NK0UwdBFEgqJBrKArinRjiwW7qFFDNJYE3zdqjL0lKkaDxgbEFlAB62vFFguEgKiAoqBYKWJBWWmiFGHv94/nLM4Os+wOO7uzD/v7XNdcO3POc55zn5nZ35x55swZc3dERCQOdfJdgIiIVJxCW0QkIgptEZGIKLRFRCKi0BYRiYhCW0QkIgptqRQzu9rMZpZ1u4xl7jCz8ble96Yi2a6vzczNbGAO+uuQ9FWYg/JqLDPrm2xny3zXUpUU2rWUmT1tZq+UMa9z8uQ/cCO6HgL0qVx169VTVujkfF0bqKG1md1qZp+a2Soz+8LMnjezQ3K8nl2Aq4CzgTbAwznodn7S1/Qc9FWmlNBcZmabp80reU5lFapmNtLMxlaw+WuE7fw2i7Kjo9Cuve4F+plZhwzzTgfmAS9n26m7L3f3avmnqa51JffRW8BBwCXAbsD+wLPAP3O8uu2Tv0+6+wJ3X1HZDt19bdLXmsr2VUFLgOPSpp0OfFZVKzSzzdx9dbKdm/Y3Bt1dl1p4AQqAr4C/pE3fDPgauBKoSwj3OcAK4GPgQqBOSvurgZkbuF2XsEe8OLncAgwHxqe0ORj4bzJ/EfAi0DllvqddxpexrjrAFYQ9y1XADOCIlPkdkuWPAV4CfgDeBw4o5756DvgCaJxhXtOU6+2BJ4DvksvjQLv0+wY4Afg0afMk0DJlfqltTaaPBMamrTd923cFXgGWAcuBd4B+adtdmNK+N/AGsDJ5vP8B1EuZPx4YBgwGioBvksexzgbup77Jev4K/CfDc+ovyfyS7d3g8yvT/ZGso2R7TgReTZY9J2X9Jf3fC7wHNExZ33/T78vYLtrTrqU87HXdBww0s9TnweFAS2AEIQS/AI4HOgOXAZcCv81iVX8CzgDOAroT/nF+ndamESHMuxH+8ZYCz5hZvWR+t+TvwYS3v0eXsa4/An8GLiKE2BPA42bWJa3ddcBtwO7AFOAhM2ucqUMza56sd6i7L0+f7+5LknZ1gKeA1kC/5NIWeNLMLGWRDsAA4CjgQKBrUg+EUDwjud4muVTUGMKLcDegCyHwVpaxTVsDzwNvJ+s/nRCA16c1/TWwBtiXEIrnJrWXZxTQzcw6JrcPI7yQjE9rV97zawjwCOEdX8n98VrK8tcTXlh2Irz4pftfwgvGkOT2ZcAOwGkV2IaaK9+vGrrk70J4AjtwYMq0Z4HnN7DM34CXU25fzYb3tL8ELku5XQf4iJQ97QzraASsBXomtzuQtqdYxrq+AK5MazMeGJXWz1kp87dOpvUso5ZuyfyjyrkvD0hq7pAybTugGNg/pd6VQJOUNpcBn6TcPpZkDztl2kjK39NeBvymjNpK3X+EF4mPKf2OaSDh3cnmKffb5LR+XgLu2cB90DdZT0vCWPx1yfSxwOWk7QlX8PmVadtLtudPZa0/ZdpewGrC3v+PwK/y8b+Wy4v2tGsxd/8Y+A/JnoeZtSWM295b0sbMzjazqWa20MyWA+cRhgHKZWZNCHtHk1PWWUx4W57arqOZjUk+5FtGeCtdp6LrSfr4GWHPdlLarImEPbFU76Zc/zL5u2VZXVewhM7Al+4+t2SCu89O+k9d/zx3X5q2/rLWnY2bgXvM7FUzu8zMdiyn1teTx6LERKAeP42pQ+n7Kdta7wVONbNtCC9oIzM1qszzC5haXgN3n0J4kboCuMvdn69g3zWWQlvuBY5MhgEGEsaUnwIwswGEYYuRhDDvQng7Wm/9biplLNCKMISyN+Et+5ocrif9g6kf181Idsco+3/h42T5zjla/48Z5pX3f1jM+i8em5XqxP1qfhom2Bd418w2ZhigsrWWeJlQ9/3Aq+7+eXqDHDy/vi+vQTI01ZPwLqhj2lBVlBTa8hjhLfvJhD3u+9295J+1J/CGu9/h7m+5+ydAxzL6WU+yR/kVsE/JtOSfplvK7RbAjsBgd3/Z3WcBWxA+KC2xOvlbdwPrWkbYE+yRNqsn4cPGjeLuJR+MnpNp3NvMmiZXZwFtU4/GMbPtCHv/G73+xELWH9/ukqHWj939Nnc/lPBi/Lsy+psF7JP2WUZPwv38aSVrLamlmBDGfUl555amIs+v1Wzgca+A84E9CB+87gP8TyX6qhEU2rWch0PKxhDGSDtS+h/sI2APM/uVme1gZleQ/XHRtwIXmtmxZtaJsGeVGkCLCUcnnGFm25tZH8JhdKmHp31DOELgoOR46SZlrOvvwAVmdqKZ/cLM/gr04qcPojbWIMKe7lQzO87MOpnZjmb2e34aQng5uT7azAqTY8pHEw4VfLWS638V6GpmpyX30YWkvDiZWUMzG5ocJ93BzPZmwy9WwwgvJsOS46cPJYwl3+HuP1Sy1lTXEt5BPV7G/Io8v+YCuyT3eUsz22y9XspgZrsThkbOcPfXgD8AN5jZztluSE2i0BaAe4BmwGvJnm6JOwmf3o8hHGXRAbgpy75vIhyJcg9hLLsOIcyAdXtkAwjHPs8EhhLGH1eltFlDOBLgd4S96afKWNdthOC+MenrKOAYd38ny5pLScam9yB8EHcDIZxfBfoDZyZtHDiCsFc8LrksAI5MGYLZ2PW/SDhc7jpgGuFxGJbSZC3h8RsJfEg4amYyYS8zU39fAL8iDENNB/4FPEg4ciNn3P1Hdy9KGztPVZHn192EdwZTCfdt+jupjMysAeF5NsbdH0/qGUN4ZznGzOpntzU1h1Xy+SQiItVIe9oiIhFRaIuIREShLSISEYW2iEhEFNoiIhEpKL+JbKymzVt4m3YV/ia21DKN6lXmOyOyqZs2bVqRu7dKn67QrkJt2rXn/mfG5bsMqaH22Kas7wiJQEHdOvMyTdfwiIhIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRUWiLiEREoS0iEhGFtohIRBTaIiIRKch3ARK/L+d/xpG9ds84r/+AU7j8htsAeP/dt3nhiUeYOnkCX87/jAYNN6djp84MHHQ+e+3buzpLljybOnUqY0aPYty4ccyZM4dGjRqx8y67cPHFl7Dffvvlu7waTaEtOdP7gEP45SH9S01r9/Pt1l2/f/gtvPXGJPodfDjHnXoGP/zwPWMfHc2gk47gksH/4KiTBlZzxZIvN95wAxMm/Iejjz6aP/xhEMu/X859I0dy4AH7M3z4PznjzDPzXWKNZe6e7xo2WZ136+r3PzMu32VUuZI97d+e8yd+f8HlZbZ7Z+rrdN61K/Xq1183beXKFZx8SG+WLPqWF6Z+REFB7dmP2GObJvkuIW8mTZpEYWEh9VOeCytWrGDPPbpSVFTEl18tqFXPhUwK6taZ5u6F6dM1pi05tXLlClauXJFx3u6F+5QKbIAGDRrSc7+DWLZkMd8u/Lo6SpQaoEePHqUCG6Bhw4YccuihLFq0iAULFuSpspovytA2s75mNja53t/MLq7GdXcxs0Oqa30xeXjEnfTesS29d2zLMX335LEH7qnQckVff0XdggK2+Fnt3fOU4Ksvv6KgoICmTZvmu5QaK/r3H+7+NPB0Na6yC1AIPFeN66zR6tQx9urRh74HHcpWbdux8OsFPPXwA9x4xZ/5cv48/vfSa8pcds4nHzLuxbH03v9XbN6ocTVWLTXNrFmzeOKJxzn88P40bqznQlnyNqZtZh2AF4DXgX2BKcAI4C/AlsCvk6a3Ag2AFcBv3f1DM+sLXODuh5nZQKDQ3c8xs47AaKAR8BRwrrs3TtpfDRQBuwDTgJPd3c3sSuBwoCHwGnBWMn088AbQD2gKnJ7c/iRp+wVwvbs/XNY21pYx7UzWrl3LH07szztTX+excVNp9/Nt12uzfNlSTj/6IIq+WcDo5//LVltvk4dK86c2j2mnW7p0KT177MtXX33FW29Pp3379vkuKe9q6pj29sBNwI7J5SSgJ3ABcCnwAdDL3bsCVwKDy+nvVuBWd98V+DxtXlfgXGAnYDugRzL9Dnffy913IYTxYSnLFLh7t2S5q9x9dVLHw+7eZUOBXdvVrVuXk888h+LiYqZM+s9681euXMH5vzuRLz6by413jqp1gS0/WbFiBUcc0Z/Zs2fz2P89rsAuR75De467z3D3YuA94BUPu/4zgA5AE+BRM5sJ/APYuZz+ugOPJtfHpM17090/T9Y1PekfoJ+ZvWFmM4D90tbxePJ3Wkr7DTKzM81sqplNXbKoqCKLbLJKgnjJ4m9LTf9x9WouPPMUZrw1hcFD/8We3XvmozypAVavXs0xRx/N65Mn89BDD9O3b998l1Tj5Tu0V6VcL065XUwYb78GGJfsBR9OGCbJxbrWAgVm1gAYBhyb7J3fnbaOVantK7ISd7/L3QvdvbBp85aVKDd+n8+bA0CzFq3WTVuzZg2XDPotb04cx1U3DaP3AfpMt7Zas2YNJwwYwMsvv8SIESM5vH//8heSvId2eZoQxo4BBlag/evAMcn1EyrQviSgi8ysMXBsBZb5DtiiAu1qjUVFC9ebtmrlSkYMvZm6BQXs06sfAMXFxVx17plMeOk5Lr7uZg4+4rjqLlVqiOLiYk495RSefvophg0bzoknnZTvkqJR048euRG4z8wuB56tQPtzgVFmdhnhQ86lG2rs7kvM7G5gJrCA8GFoecYBF5vZdMr5ILK2uP36q5g3+xP27tWX1m225tuF3/DcEw8zf86nnH3BZeuGSW697gpeGvsEe+zdg/oNGvD8E6Xvum49+9Gi1Zb52ASpZn/+8wU88sjD9O7Th4YNGzJ61KhS8/c/4ABat26dp+pqtk3qG5FmtjmwIjn64wTgRHc/Il/11JajR1586jGefPA+5nzyEcuWLqZBg4Z02nk3jh94Jv0OPnxdu7MHHMZbb0wqs5/hDz5Tq8a3a/PRI/vt148J/1n/A+oSL7/yaq0f3y7r6JFNLbR7AXcABiwBTnP3T/JVT20Jbdk4tTm0pXxlhXZNHx7Jirv/F8h8ujkRkU1ATf8gUkREUii0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiklVom9kYMzvAzKyqChIRkbJlu6d9AvAC8JmZDTazTlVQk4iIlCHb0J4IOLA1cBHwvpm9bmZnmVnTXBcnIiKlZRXa7t6bENiDgPFAMdANGAZ8aWYPmllhrosUEZEg6w8i3f1rdx8OHEQI7++TWQ2AAcDrZvbb3JUoIiIlsg5tM+tsZkOAz4HhQKNk1svA7YS976tzVaCIiPykIJvGZjaZMBwCYMBS4D5gmLt/lLTpAByWwxpFRCSRVWgDeyd/3wWGAqPd/Ye0NmOBRZUtTERE1pdtaD8EDHX3SWU1cPe7gbsrVZWIiGSUVWi7+0lVVYiIiJRPX2MXEYmIQltEJCIKbRGRiGR7wqidUs83kpw8apSZXWJmdXNfnoiIpMp2T/tfQFcAM9sGeApoTvhm5LW5LU1ERNJlG9o7Am8l148F3nD3Q4BTgBNzWZiIiKwv29CuC6xOrv8SeC65/inQOldFiYhIZtmG9kzg92bWixDaLyTTtwaKclmYiIisL9vQvgg4g3Ba1gfdfUYyvT/wZg7rEhGRDLL9RuQEM2sF/MzdF6fMuhNIPweJiIjkWLbnHsHd1wKL06bNzVVBIiJStnJD28yermhn7t6/cuWIiMiGVGRP+9sqr0JERCqk3NB2d/10mIhIDaFzj4iIRCTrDyLNrIDwk2PtgXqp89z9/hzVJSIiGWT7G5E7As8A2xJ+I3Jt0sePwCpAoS0iUoWyHR65BZgGNCEcl90ZKASmA8fksjAREVlftsMjewF93P17MysGCtz9LTO7ELgd2C3nFYqIyDrZ7mkbP33zcSHhnCMAnwPb56ooERHJLNs97ZnA7sBswrlGLjKztYTzkXyS49pERCRNtqF9HdAouX458CwwjnCGv+NzWNcmoVG9uuyxTdN8lyE11IsTZ5TfSCRNtieMejHl+mygs5k1Bxa7u+e6OBERKS3r47TTufuiXBQiIiLly/Y47Q2ePEonjBIRqVrZ7mmnnzxqM8IHk9sAj+ekIhERKVO2Y9oZTx5lZjcBy3JSkYiIlClXJ4y6ExiUo75ERKQMuQrtTjnqR0RENiDbDyJvS58EtAF+BfwrV0WJiEhm2X4QuWva7WLC19nPQ6EtIlLlsg3t3wCfu3tx6kQzM8IRJJ/lqjAREVlftmPac4CWGaY3T+aJiEgV2piz/GXSGFhZyVpERKQcFRoeSfkA0oHBZvZDyuy6hJ8fm57b0kREJF1Fx7RLPoA0wq/VrE6Ztxp4CxiSw7pERCSDCoW2u/cDMLMRwB/dXd9+FBHJg5x8jV1ERKpHVh9Emtl1ZnZ2hulnm9k1uStLREQyyfbokVOAtzNMnwacWvlyRERkQ7IN7S0J34BM9y3QuvLliIjIhmQb2p8BvTJM7034RXYREalC2X6N/U7gH2ZWD3g1mfZL4HrghlwWJiIi68v26JGbzKwlcBtQL5m8GrjV3W/MdXEiIlJa1ufTdvdLCOcf2Se5tHL3i81s/1wXJyIipW3UjyC4+/fuPgX4EjjXzD4FXsxpZSIisp6sQ9vM6prZ0Wb2LDAXOIow1r19jmsTEZE0FR7TNrNOwO8Ix2N/D4wBDgROcff3q6Y8ERFJVaE9bTP7L/A60Aw43t23c/fLq7QyERFZT0X3tLsDQ4G73P29KqxHREQ2oKJj2nsRAn6imb1tZueZ2VZVWJeIiGRQodB297fdfRDhl9dvBvoD85PlDzWzZlVXooiIlMjq6BF3X+nuDyTn1+4M/J3wS+wLzOz5qihQRER+slHHaQO4+yfufjHhV9iPp/Sv2YiISBXI9twj63H3tcBTyUVERKrQRu9pi4hI9VNoi4hERKEtIhIRhbaISEQU2iIiEVFoi4hERKEtIhIRhbaISEQU2iIiEVFoi4hERKEtIhIRhbaISEQU2iIiEVFoi4hERKEtIhIRhbaISEQU2iIiEVFoi4hERKEtIhIRhbaISEQU2iIiEVFoi4hERKEtIhIRhbaISEQU2iIiEVFoi4hEpCDfBcima/ny5QwZ8nemTpnC1KlTKCoq4pJLL+Oaa67Nd2lSzRZ+s4DRI4bz7ltvsHjRtzRr0ZKuhd054dQzabXlVuvafTZ3NqNHDOOD997hu2VLadGqNd177cdxJ53GFj9rksctqDkU2lJlioqKuPaav9KuXTu6dOnKyy+/lO+SJA+WLV3CeWedxNq1azjkiOPZsnUbPps7m+effpQpr0/gn/c9yeaNGvP5Z3M47+wT2WKLJhxyxPE0adqMjz98nyceuZ+3p0zmlrsepG7duvnenLxTaEuVadOmDZ/N/4K2bdsyd+5ctu+4bb5LkjyY8OoLLF5UxJWDb2PvHn3XTW+9VVvuvP0G3pryGj37Hsi/n32ClStWcNPQB+jQ8RcAHHw4NGjQkCcffYDZH3/ADjvunKetqDk0pi1Vpn79+rRt2zbfZUie/fDD9wA0b9Gq1PRmye36DRqWatcsrV3zFi2Tdg2qtM5YRBvaZtbBzGbmoJ9CM7stFzWJyPp236MbAP+89XrenzmdooVf8/aUydx/z+3suNNu7FHYPbTrGtrd8rcr+OSj9yn6ZgGvTXiF/3toJD367E/7Dh3ztg01Sa0fHnH3qcDUfNchsqnq1HlX/nDeZdx/z+38edCp66Z327cPF115A3ULQgz17HcgJ879hMcfvo83z5iwrt1Bhx3DoPMvr/a6a6rYQ7vAzEYDewDvAacCnYGbgcZAETDQ3b8ys/HAG0A/oClwurv/18z6Ahe4+2Fm1goYA7QFJgMHAHsmfT0PTAT2Bb4AjnD3FdWzmSJxa9mqNTvuvDtd9tyHNlu3Y+6nH/N/D43kr5f9kauvv4N69etjZrTeamt27bIX+/ToS9NmLZj13js89dgoVq1ayQWXDcbM8r0peRd7aHcihO8kM/sXMAg4ihCoC81sAHAdcFrSvsDdu5nZIcBVwP5p/V0FvOru15vZwcDpKfN2AE509zPM7BHgGGBU1W2ayKZh8sRxXH/ln7j93kf4+bbbA7BPj350/EVnrr5oEM89/QhHHncKTz82mlEjhnHnA0/TrHkLALr32o8tW7dh+C2D6dnnALr32i+fm1IjRDumnZjv7pOS66OAg4BdgJfMbDpwOdAupf3jyd9pQIcM/fUEHgJw9xeAxSnz5rj79HKWx8zONLOpZjZ14cKFWW6OyKbnqUdH0bZd+3WBXaJw757Ub9CAmdOnAfDkow/Qeefd1wV2iR59DgBgxjsaxYT4Q9vTbn8HvOfuXZLLru5+YMr8VcnftWT/LmNVyvUyl3f3u9y90N0LW7VqlamJSK2y6NtvKC5eu9704uJivNhZs3YNAN9+u5Di4uL12q1N5q9ds6ZqC41E7KHd3sy6J9dPAl4HWpVMM7PNzCybAzsnAccnyx4INMtlsSK1Ubv22/Ll55/xwfvvlpo+cfy/Wb16FTt02gmAbdpvy3sz3uLrr74o1e7Vf48FYPtOOkYb4h/T/hAYlIxnvw/cDrwI3GZmTQjbdwvhQ8qK+AvwoJmdQvggcgFh771xjuuuNYYOvYMlS5awZMkSACZNmsh114WvsR9+eH922223PFYn1eHYk05j2hsTufxPZ3HokQPYqk075s7+iBeeeYzmLVpx6JEDABhwypn87eoLOP8PJ3PoEcfTrHlLZs2czqv/HsvPt+1In/0OzvOW1Azmnj7CUHuZWX1grbuvSfbWh7t7l43tr7Cw0N94s3aPw3XcrgPz5s3LOO/ee0fwm4EDq7egGuTFiTPyXUK1mfPpRzx43z/56IP3WPztQrb4WVP22Ks7J59+Dlu2brOu3cx3pvHIqHuYO/tjli5ZRLMWrdh73z6cfNqgWnfukUP77DbN3QvTp8e+p51r7YFHzKwOsBo4I8/1RO/T2XPzXYLUANt2/AWX/vXmctvtsvue7LL7ntVQUbwU2inc/WOga77rEBEpS+wfRIqI1CoKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCIKbRGRiCi0RUQiotAWEYmIQltEJCLm7vmuYZNlZguBefmuowZpCRTluwipsfT8KO3n7t4qfaJCW6qNmU1198J81yE1k54fFaPhERGRiCi0RUQiotCW6nRXvguQGk3PjwrQmLaISES0py0iEhGFtlSamfU1s7HJ9f5mdnE1rruLmR1SXeuTzMysg5nNzEE/hWZ2Wy5q2lQV5LsA2bS4+9PA09W4yi5AIfBcNa5Tqoi7TwWm5ruOmkx72gKs21P6wMxGmtlHZjbazPY3s0lm9rGZdUsuk83sbTN7zcw6ZehnoJndkVzvaGavm9kMM7vWzJYn0/ua2XgzeyxZ52gzs2TelWY2xcxmmtldKdPHm9kNZvZmUl8vM6sH/BUYYGbTzWxA9d1jkkFB8ljOSh7bzc1sTzP7j5lNM7MXzawNZH48k+mp79pamdlLZvaemd1jZvPMrGXyXJ1lZncn8/5tZg3zueHVSaEtqbYHbgJ2TC4nAT2BC4BLgQ+AXu7eFbgSGFxOf7cCt7r7rsDnafO6AucCOwHbAT2S6Xe4+17uvgvQEDgsZZkCd++WLHeVu69O6njY3bu4+8NZb7HkUidgmLt3BpYBg4DbgWPdfU/gX8B1Ke1LPZ4Z+rsKeNXddwYeA9qnzNsBGJrMWwIck9tNqbk0PCKp5rj7DAAzew94xd3dzGYAHYAmwH1mtgPgwGbl9NcdODK5PgYYkjLvTXf/PFnX9KT/iUA/M7sQ2BxoDrwHPJMs83jyd1rSXmqW+e4+Kbk+ivBCvwvwUvKGqS7wVUr78h7PnsBRAO7+gpktTpk3x92nl7P8JkmhLalWpVwvTrldTHiuXAOMc/ejzKwDMD5H61pLeGvdABgGFLr7fDO7GmiQYZm16LlbE6UfP/wd8J67dy+jfWUez/Tnj4ZHRDJoAnyRXB9Ygfav89Pb1hMq0L4koIvMrDFwbAWW+Q7YogLtpOq1N7OSgD6J8Pi3KplmZpuZ2c5Z9DcJOD5Z9kCgWS6LjZVCW7JxI3C9mb1NxfaMzgXON7N3CePlSzfU2N2XAHcDM4EXgSkVWMc4YCd9EFkjfAgMMrNZhIC9nfDCe4OZvQNMB/bNor+/AAcmhxIeBywgvEjXavpGpFQZM9scWJGMi58AnOjuR+S7LomDmdUH1rr7mmRvfbi7d8lzWXmncUGpSnsCdySH7S0BTstvORKZ9sAjZlYHWA2cked6agTtaYuIRERj2iIiEVFoi4hERKEtIhIRfRApIhmZ2X5Ab8JXySfkux4J9EGkiKzHzNoB7xLON7M9sKu7f53fqgQ0PCKSNTMba2YjU26PLDkzXSX6rHQfuZIcojkSGO3u+wJPEU72JDWA9rRlk5CE6G+Sm2uA+YQTEl3l7t/neF1jgSJ3H5jcbkL4X1pSweXHAzPd/ZyUaVn1IbWXxrRlU/IycArh7IO9gHuARsDv0xuaWb3k1K6V5u4b/Hp+dfUhtYOGR2RTssrdF7j7fHcfA4wmOTVsctL94WY2xMwWEk5GhAUXmtmnZrYi+cGGk0s6TE7kP9LMlpvZ12Z2afpK04c2kj7/ZOHHI1aZ2edmdn1JW6AP4Rwdnlw6ZOijvpndkqxzpYUfk+iZtt7xZjbMzAabWZGZfZNsX52UNr2TZZeb2dLkRwd2ybANDya1vJE2fUIy/aGsHgmpMgpt2ZStoPQ5v08GjLAXfmoy7VrgdMIJ+3cCrgfuNLNDk/lDgAMIZyv8JeHHG3qXs97BwBVJXzsTTnY0P5n3R2AyMAJok1zmZ+jjRmAA4av/XYEZwAuW/PJLil8ThoP2Bc4hnKRrAICZFRDGoycCuwN7A7cQTmWabmjyt1vJmfjMbCt++nGKkeVss1QXd9dFl+gvhFAZm3K7G1BE+FUbCOf+fjdtmUaEYO+VNv0Wwm9ONiact/nXKfMaE86jMjLTupP5K4GzN1DreMIv9GSsP6lrNXBqyvy6wKfAtWn9TE7r5yXgnuR6c8I5rvtU8D58J2l/c3L798ntL4C6+X6MdQkX7WnLpuTgZBhgJWFvdgLwPynzp6W134lwDu8XkuWWW/gdy98DHZNLvaQvANx9OWGvtyw7AfWBVyqxHR0J7xBKfgUGd1+b1LFTWtt3025/CWyZLLOI8GLwopk9a2bnm1l7ylayt32ymW3GT+czH5WsX2oAhbZsSiYQfp29E9DA3Y92929S5qcfRVLy/D88Wa7ksjNwYBXWWRnph3v9mGH+uv9rd/8tYVhkAtAf+NDMDiqj79GEc563IgzL9Emmj6xcyZJLCm3ZlPzg7p+4+zx3Tw+zTN4nDH/8PFku9TKPMBzxI7BPyQJm1ojwu4dlmZX0+csNtFlNGO4oy6dJm5LxZMysLuE3N9/f8Catz93fcfcb3L0vYUjlN2W0+x64L7l5U1LjFHefle06perokD+ptdz9OzMbAgxJvlAygTAmvQ9Q7O53mdm9hF9eWUgYeriSDQRu0uethF/4WZX02QLY092HJ83mEj7w6wAsBxal9fG9mQ1P1lsEzAHOA1oTfkOzQsxsW+As4GnCuPR2wG7A8A0sNowwpNQouT2youuT6qHQltruCuBr4AJCmC0j/CzWjcn8CwgB9gTwA+EntBqt10tplwCLk77bJf3fnzJ/CGGP9n3CD9Jum6GPi5K/I4CmwNvAwe7+VYa2ZfkB+AXwKNAyqWM0cENZC7j7h2b2CrA/4R2DDvWrYfSNSBEpxcxuI+xtP+rux+e7HilNe9oiAoCZnQkcChxC+EBzSH4rkkz0QaSIlNiXcITJ18CZ7v5mnuuRDDQ8IiISEe1pi4hERKEtIhIRhbaISEQU2iIiEVFoi4hERKEtIhIRhbaISEQU2iIiEVFoi4hE5P8BtVrrm6Fn/y0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# compose confusion matrix\n",
    "from sklearn.metrics import confusion_matrix\n",
    "conf_mat = confusion_matrix(y_true=val_y, y_pred=pred_y)\n",
    "\n",
    "# print confusion matrix\n",
    "fig, ax2 = plt.subplots(1, 1, figsize=(15,5))\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('Validation 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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loss animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "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 $F(\\mathbf{w})$', fontsize=16)\n",
    "ax.xaxis.set_tick_params(labelsize=12)\n",
    "ax.yaxis.set_tick_params(labelsize=12)\n",
    "line, = ax.semilogy(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",
    "# animation\n",
    "div = 50\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=div, frames=epochs//div+1)\n",
    "\n",
    "if False:\n",
    "    anim.save('./img/sgd_anim.gif', writer='imagemagick', fps=5)\n",
    "\n",
    "    HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "98.11 % Precision\n",
      "96.30 % Recall\n",
      "97.90 % 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": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSEs are close enough with 0.2238 (test) and 0.2582 (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_val_y = np.round(sigmoid(np.dot(val_X, w))).astype('uint8')\n",
    "pred_train_y = np.round(sigmoid(np.dot(train_X, w))).astype('uint8')\n",
    "\n",
    "# compare MSEs\n",
    "val_MSE = MSE(val_y, pred_val_y)\n",
    "train_MSE = MSE(train_y, pred_train_y)\n",
    "res = np.isclose(val_MSE, train_MSE, rtol=.95)\n",
    "\n",
    "if res:\n",
    "    print('MSEs are close enough with %s (test) and %s (train).' % (val_MSE, train_MSE))\n",
    "else:\n",
    "    print('Potential over-/underfitting from MSEs with %s (test) and %s (train).' % (val_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": 10,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "epoch 0\n",
      "            Train set - loss: 0.445, accuracy: 12.587\n",
      "            Validation  set - loss: 0.448, accuracy: 12.587\n",
      "            \n",
      "epoch 100\n",
      "            Train set - loss: 0.102, accuracy: 93.007\n",
      "            Validation  set - loss: 0.087, accuracy: 93.007\n",
      "            \n",
      "epoch 200\n",
      "            Train set - loss: 0.081, accuracy: 95.105\n",
      "            Validation  set - loss: 0.073, accuracy: 95.105\n",
      "            \n",
      "epoch 300\n",
      "            Train set - loss: 0.073, accuracy: 95.804\n",
      "            Validation  set - loss: 0.067, accuracy: 95.804\n",
      "            \n",
      "epoch 400\n",
      "            Train set - loss: 0.069, accuracy: 96.503\n",
      "            Validation  set - loss: 0.064, accuracy: 96.503\n",
      "            \n",
      "epoch 500\n",
      "            Train set - loss: 0.066, accuracy: 96.503\n",
      "            Validation  set - loss: 0.063, accuracy: 96.503\n",
      "            \n",
      "epoch 600\n",
      "            Train set - loss: 0.065, accuracy: 96.503\n",
      "            Validation  set - loss: 0.061, accuracy: 96.503\n",
      "            \n",
      "epoch 700\n",
      "            Train set - loss: 0.064, accuracy: 96.503\n",
      "            Validation  set - loss: 0.06, accuracy: 96.503\n",
      "            \n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "torch.manual_seed(1111)\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=True)\n",
    "        torch.nn.init.uniform_(self.linear.weight, 0, 1e-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",
    "val_Xt = torch.autograd.Variable(torch.FloatTensor(val_X))\n",
    "val_yt = torch.autograd.Variable(torch.FloatTensor(val_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.MSELoss()\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",
    "        pred_y = model(i_X)\n",
    "        loss = criterion(pred_y, i_y)\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "    # track loss\n",
    "    loss_list.append(loss.item())\n",
    "    if epoch % 100 == 0:\n",
    "        y_val_pred = model(val_Xt)\n",
    "        y_train_pred = model(train_Xt)\n",
    "        val_loss = criterion(y_val_pred, val_yt)\n",
    "        train_loss = criterion(y_train_pred, train_yt)\n",
    "        pred_yt = sigmoid(y_val_pred.detach().numpy()) >= tau\n",
    "        conf_mat = confusion_matrix(y_true=val_yt.numpy(), y_pred=pred_yt)\n",
    "        val_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",
    "            Validation  set - loss: {np.round(val_loss.detach().numpy().astype('float'), 3)}, accuracy: {np.round(val_acc, 3)}\n",
    "            ''')\n",
    "\n",
    "if False:\n",
    "    torch.save(model, 'torch_sgd_model.pth')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# threshold estimation\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "fpr, tpr, thresholds = roc_curve(train_yt, model(train_Xt).detach())\n",
    "tau = thresholds[np.argmax(tpr - fpr)]\n",
    "\n",
    "# compute predictions from test set\n",
    "pred_yt = sigmoid(y_val_pred.detach()) >= tau"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "conf_mat = confusion_matrix(y_true=val_yt, y_pred=pred_yt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAFXCAYAAACV5h77AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAAsTAAALEwEAmpwYAABUeElEQVR4nO3dd3gc1dXH8e9RL5Yly70XbINNs0E4BkwLzQRCCSSYBEIvARMglQRCCAk1hNATmikvLUBIcKihGUK3DQY3bNzA3XIvslXP+8eM7PVasiRb2pG0v8/z7KOZO3dmzkirnT1z79wxd0dERERERERqlxJ1ACIiIiIiIs2dEicREREREZE6KHESERERERGpgxInERERERGROihxEhERERERqYMSJxERERERkToocRKJiJk9YmYvNnCdcWZ2d1PFJCIiLYeZXWtmU2qbr2Wdu81sXGPvu7UIj2upmbmZndUI2+sTbquoEcJrtszs0PA4O0QdS1NS4iRSh/CDYHuvR3Zw05cBpzdwne8Bv9nB/dVbaz0hiog0B2Y21szerGXZoPDcctQObPpW4JCdi26beGr74t/o+9pODJ3N7A4zm21mpWa20MxeMbPvNPJ+9gB+D1wEdAX+0QibnR9ua1IjbKtWMYnLWjPLiVs2KOY7S70TmwZe4P2A4DhXNCDsFict6gBEWoCuMdPHAQ/ElW2MrWxm6e5eXtdG3X1NQwNx95UNXUdERJqdh4B/mVkfd58Xt+xc4GvgjYZu1N3XA+t3Przmsy8z6wO8D6wjuHD4OcGF/8OBvwO9GnF3/cOf/3Z3b4wNunslsKQxtlVPq4HvA4/GlJ0LfEPj/q42C7/3lJHY44yEWpxE6uDuS6pfBB9IxMxnAavN7DQze8vMNgIXmll7M3vKzBaY2UYzm2pmZ8duN/5KTtgN714zu8HMlpvZMjO71cxS4urcHTM/z8yuNrP7wqtMC8zsl3H7GWhm75jZJjObYWbfMbP1O9MFwcz2NLM3wmNbGR5LftzyN8OY1pvZ52Z2WLgs3czuNLNF4ZXD+WZ2047GIiLSAr0ELAXizwvpwBnAmGDWHjKzueFn7Vdm9qvYc0K8GrrupYbnkVXh63YgNW6dkWb2v3D5SjN7zcwGxVSZG/4cH7ZYjKtlXylm9rvwM73UzCab2Qkxy6tbrk42s9fNrMTMppnZkXX8ru4Nfxa5+zPuPsPdp7v73cBeMdvvZWb/MrN14et5M+sR/7sxs1Fhy9U6M/t3dQuMmV0L/CusXmVmHpZv0+pSw7Fv75y3TYudmR1sZh+H5+WlZvZXM8uIWV7n94HteAQ4J2Zb1e+pR+KOIXV776/w93EmcKxtaa06NOZ44r/3bNVVL9z2VDPLjtnf/+J/ly2NEieRxnEjwYf7YODfBAnVpwQtVLsDdwD3mdnhdWznR0AFcAAwGrgcOLWOda4AJgP7ADcDt5jZ/hCcyAhOBBXAcOAsgm4ImQ04tq2YWS7wGsGVxmHASWG8Y2KqPQksDpcPAa4FNoXLfhquMwoYEB7fjB2NR0SkpXH3CoIWgbPivgx/F+gAPEzwHW0h8ANgEHAV8Fvikq06/Bw4H7gQ2J8gafpRXJ1c4HaCz+tDgTXAf2K+yA8Lf44k6G3xvVr2dRnwS+DXwJ4E557nzWxIXL3rgTuBvYHxwNNm1qamDZpZYbjfe8IWrq24++qwXgrwAtAZOCx8dQP+bWYWs0ofgnPOScBRwNAwHgi6Hp4fTndl654lddneOS/+mLoDrwCfhfs/FziN4HtErB35PgDwODDMzHYJ548jOF+Pi6tX1/vrVuAZgpbP6t/HBzHrx3/vifdTID3cDuH2BxCT1LVI7q6XXnrV8wWcEvzbbJ7vAzjw83qs+zTwYMz8I8CLMfPjgA/j1nk9bp1xwN0x8/OAp+LW+Qq4Opw+muCDt3vM8gPCmM/aTqzXAlNqWXY+wYk1L6bs0HCb/cP5tcCZtax/J/AmYFH/PfXSSy+9onoRfIl04KiYspeAV7azzk3AGzHzW31W1zC/CLgqZj4FmAmM284+coFKYEQ4X32eK4qrF7+vhcA1cXXGAY/HbefCmOXdw7IRtcQyLFx+Uh2/yyPDmPvElPUDqoAjYuLdBOTH1LkKmBUzv9U5Pizb6lxdy7Fv75y31e+PIFH7CkiJqXMWUArkxPzetvt9oIb9VJ+HOxDcm3V9WP4icHXs8ga8v2o69urj+Xlt+48p2w8oA64DyoFjovhfa8yXWpxEGseE2JmwSfoqM/vCzFaY2XqCq3R19S/+Im5+EdBpJ9bZDVjk7gtjlo8nOJnsqEHAF+6+Lqbsg3Cbg8P524AHw2b8q8xst5i6jxBckZtpZveY2bH17H4gItJquPtXwDuEV+DNrBvBxa6HquuY2UVmNsHMisPzyBXU8z4VC7pPdwU+jNlnFfBxXL1dzOzJsPvaWoIuhCn13U+4jbYELTzvxy16jy3nhWqx56xF4c/aznNWS3m8QQTnunnVBe4+J9x+7P6/9q3vL67PObY+tnfOqynWj8K/RbX3gAy23GMFO/Z9oNpDwI/NrCdBUvlITZV25v1F3Peemrj7eIJE8XfA/e7+Sj233Wzpy4pI49gQN/8Lgi4Sfya4gXUIQVN2BtsXP6iEU/f/6Y6s01QcwN2vZUvz/QHAF2Z2TrjsU4IrVr8hiPNR4HUlTyKShB4CTgy7pJ0FrCTocoaZnUrQhe4RgoRqCEHXqLrOIw31ItCRoDvftwi6j1U04n7iB1nYfM7ysFmC2s9ZX4XrD6pleUP3vyPnyyq2TeDSt9rIds55DbSzsVZ7gyDux4C33H1BfIVGeH/Ff+/ZRthNcgRBa+Aucd0mWyR9URFpGiOA/7j7/7n7JGA2MDCCOL4EuoVXMqsVsXP/+9OBPc0sL6bsgHCb06sL3P0rd7/T3Y8l+HJwXsyyde7+nLv/BDgW+DZbX2kTEUkGzxF0HzudoOXpMd8yKusI4GN3v9vdP3X3WcAutWxnG2HLymKC+1uBzV9kh8XMtyfomXCDu7/h7tOBPLYedbks/LnVoBJx+1pL0CJyYNyiEcC0+sZcw3ZXEtxTO7qm+6DMrCCcnE5wrusTs6wfQSvYDu8/VMy29zsNqSHWWs95caYDw+MuFo4g+D3P3slYq2OpIkiIDiWmBTNOfd5fZWzn714PPyO4//pggvfhpTuxrWZBiZNI05gJHG5mI8Im+7uBvhHE8TrBwAuPmtneZjacoEtBBdteBYyXZWZD4l4DgSeAEuCxcCShg4H7gOfdfZaZZYdd8KpH3/kWMSdPM/tZOBrPIDPrD/yQoH/4NlfERERaM3ffSDCwwLUEX1pjv+TOBPYxs2PMbICZ/Y6GPzfpDuBXZnaKme1K0MIQmwSsApYD55tZfzM7hGCI74qYOssIHrtxtAXPU8qnZn8GfhF+vg80s+uAg9gyOMCOuoSgxWeCmX3fzHY1s93M7Cds6c72Rjj9hJkVWTCC3RMEgzS9tZP7fwsYambnhL+jXxGTINZ1zqvBvQQJ3b3hefBYgnuL7nb3kp2MNdafCFoSn69leX3eX/OAPcLfeQcLRuirFzPbm6Cb3vnu/gFwMXCzme3e0ANpTpQ4iTSNPwGfEIyc8y5Bk/YTiQ4ivOp0EsEoep8QdIu7niBpqnHEnxi7EIz6E/t6MvxgPxpoG27zBYI+9NXdEiqBdgRXu2YQjKz0IcGVJwiexfHLcN1PCa7cHdPIJwwRkZbiQYLPzA/CFp9q9xGMavYkwb2pfYC/NHDbfyEYoe9BgnubUog5F4XniFMJhvWeAtxDcD9KaUydCoIR0s4jaFV6oZZ93UmQPN0Sbusk4GR3/7yBMW8lvFdpH4ILgTcTJEhvAccDF4R1HDiBoHXo7fC1BDgxpjvgju7/NeAPBOfOiQR/h3tjqtR1zovf3kLgGIIukZMIRqR9imBEu0bj7uXuvjzuXqpY9Xl/PUDQQjaB4Hcb36JYIzPLInifPenuz4fxPEnQwvqkme3wyL5Rs518P4lICxNeBZpEMMLPxIjDEREREWkRlDiJtHJmdhJBi9dXBFeUbiPo9jB0Z6/EiYiIiCSLtLqriEgLl0fQvaEnQX/2ccAVSppERERE6k8tTiIiIiIiInXQ4BAiIiIiIiJ1UOIkIiIiIiJSh1Z/j1OHDh28T58+UYchIiIJMHHixOXu3jHqOCRa7Qrbe7eevaIOQ5qpzLSdeaartHaffVr7eaTVJ059+vRhwoQJUYchIiIJYGZfRx2DRK9bz1489co7UYchzdSAjm2iDkGasZyM1FrPI+qqJyIiIiIiUgclTiIiIiIiInVQ4iQiIiIiIlIHJU4iIiIiIiJ1iCRxMrORZjbDzGaZ2ZXbqXeymbmZFcWU/SZcb4aZHZ2YiEVEREREJJklfFQ9M0sF7gGOBBYA481srLtPi6uXB1wGfBxTNhgYBewOdAPeMLOB7l6ZqPhFRERERCT5RNHiNAyY5e5z3L0MeBo4oYZ6fwRuBjbFlJ0APO3upe4+F5gVbk9ERERERKTJRJE4dQfmx8wvCMs2M7N9gJ7u/lJD1w3Xv8DMJpjZhOLi4saJWkREREREklazGxzCzFKA24Cf7+g23P1+dy9y96KOHfUAeRERERER2TlRJE4LgZ4x8z3Csmp5wB7AODObBwwHxoYDRNS1bqMrr6ziin9M4tUpS5pyNyIiIiIi0oxFkTiNBwaYWV8zyyAY7GFs9UJ3X+PuHdy9j7v3AT4Cjnf3CWG9UWaWaWZ9gQHAJ00ZbFqKMfbzRUxZuKYpdyMiIiIiIs1YwkfVc/cKMxsNvAakAmPcfaqZXQdMcPex21l3qpk9A0wDKoBLmnpEPTOjTWYa60srmnI3IiIiIiLSjCU8cQJw95eBl+PKrqml7qFx89cD1zdZcDVok5nG2k3lidyliIiIiIg0I81ucIjmKC8rjfWb1OIkIiIiIpKslDjVQ15WGuuUOImIiIiIJC0lTvWge5xERERERJKbEqd6yMtKV+IkIiIiIpLElDjVQ5usNNZpcAgRERERkaSlxKkedI+TiIiIiEhyU+JUD3mZaZRWVFFWURV1KCIiIiIiEgElTvXQJjN43JXucxIRERERSU5KnOohLysdQPc5iYiIiIgkKSVO9dAmK2hx0n1OIiIiIiLJSYlTPeRlqaueiIiIiEgyU+JUD3mZ1V31lDiJiIiIiCQjJU71sKXFSfc4iYiIiIgkIyVO9aB7nEREREREkpsSp3qoHo5ciZOIiIiISHJS4lQPWempZKWnsGajuuqJiIiIiCQjJU71lJ+dzpoSJU4iIiIiIslIiVM9FWRnsHpjWdRhiIiIiIhIBJQ41VN+Tjqr1eIkIiIiIpKUlDjVU0F2uu5xEhERERFJUkqc6qlALU4iIiIiIklLiVM95avFSUREREQkaSlxqqeCnAw2lleyqbwy6lBERERERCTBlDjVU352OgBr1eokIiIiIpJ0lDjVU0FOkDitVuIkIiIiIpJ0lDjVU0F2BoAGiBARERERSUJKnOqpuqve6hI9BFdEREREJNkocaqn6q56GllPRERERCT5KHGqp3wlTiIiIiIiSSuSxMnMRprZDDObZWZX1rD8IjObbGaTzOw9Mxsclvcxs41h+SQz+3uiYs7LTCM1xXSPk4iIiIhIEkpL9A7NLBW4BzgSWACMN7Ox7j4tptqT7v73sP7xwG3AyHDZbHcfksCQCeMgPzud1Rt1j5OIiIiISLKJosVpGDDL3ee4exnwNHBCbAV3Xxszmwt4AuOrVX52ulqcRERERESSUBSJU3dgfsz8grBsK2Z2iZnNBm4BfhqzqK+ZfWZm75jZQU0b6tbys9N1j5OIiIiISBJqtoNDuPs97r4L8Gvg6rB4MdDL3YcCPwOeNLO28eua2QVmNsHMJhQXFzdaTAU5SpxERERERJJRFInTQqBnzHyPsKw2TwMnArh7qbuvCKcnArOBgfEruPv97l7k7kUdO3ZsrLgpUFc9EREREZGkFEXiNB4YYGZ9zSwDGAWMja1gZgNiZo8FvgrLO4aDS2Bm/YABwJyERA0U5GSwSg/AFRERERFJOgkfVc/dK8xsNPAakAqMcfepZnYdMMHdxwKjzewIoBxYBZwZrn4wcJ2ZlQNVwEXuvjJRsRfkpLNuUwXllVWkpzbbXo4iIiIiItLIEp44Abj7y8DLcWXXxExfVst6/wT+2bTR1a59bgYAq0rK6JSXFVUYIiIiIiKSYGo2aYDC3EwAVm5Qdz0RERERkWSixKkBCsMWp5XrlTiJiIiIiCQTJU4N0L5NkDitUIuTiIiIiEhSUeLUAJtbnJQ4iYiIiIgkFSVODdAuJwMztTiJiIiIiCQbJU4NkJpiFGSns3JDadShiIiIiIhIAilxaqDC3Ax11RMRERERSTJKnBqofW4mKzSqnoiIiDRDH7/3Dnt3z2fv7vl8M3f2NsvnzfqKX198Doft3Z/9+nXimOF7ctVPL4wgUonSgvnzueiC8xg0cBcK2+YyeNf+jL74IhbMnx91aM1aJA/AbckKczOYXbw+6jBEREREtlJeXs6NV/2C7JxcNpZs2Gb5p598yMU/Opleffrx4wtHk9+ukOIlS/hs/IcRRCtRWbFiBQcdOJyK8nLOu+AievXuxZfTp/Pg/ffx6isv8+nnU2jbtm3UYTZLSpwaqLBNBuPnqcVJREREmpfH/n4na1av4uQfnsnjD9671bKSkg1cefG5FO0/gtvHPElamr4CJqvnnn2GpUuW8Ozz/+bY4767ubx37z784meX8+br/+Wkk0+JMMLmS/81DdQ+N4NVJWVUVTkpKRZ1OCIiIiIsXjif+++4ld9c/2cWL9i2u9WrL/yTpYsXcs/jz5GWlsbGjSVkZGSSmpoaQbQSpXVr1wLQpWvXrcq7dOkCQHZOTsJjail0j1MDFeZmUOWwemN51KGIiIgknJkdamYvhtPHm9mVCdz3EDP7TqL215LcfM2VDNhtMCf84Ec1Lv9w3Ju0yWvLmlUrOfnw/RnevyvDdunMZWefxuKFuq8lmRxy2GEA/Pzyy/joww9YuHAhb77xOtde8zuGfWs4Rxx5VMQRNl9KnBpoy0NwNSS5iIgkN3cf6+43JXCXQwAlTnHefeM13vnvy/zm+j9jVnNvmK/nzKaiooJLzvg+e+1TxF8e+D8uuuLXfPTu25x7yrFsWL8uwVFLVPbbbxi333U3X82cwbcPOYgBfXvx3e+MZOCuu/LSq/9VN87t0G+mgdrnZgKwYn0Z/TtFHIyIiMgOMLM+wKvAR8ABwHjgYeAPQCegutniDiAL2Aic7e4z4rZzFlDk7qPNbBfgCSAXeAG43N3bmNmhwLXAcmAPYCJwuru7mV0DfBfIBj4ALgzLxwEfA4cBBcC54fx1QLaZjQBudPd/NOKvpUUq3bSJm6/5FSeOOp3d996n1nolJevZtLGEE0edzu//fNfm8i7de3L1ZRfy76cf50fn/SQRIUsz0L17D4YNH863v30Effv1Y8rkydx+26384OST+Oe/x5KVlRV1iM2SWpwaaEuLkwaIEBGRFq0/8Bdgt/D1Q2AE8Avgt8CXwEHuPhS4Brihju3dAdzh7nsCC+KWDQUuBwYD/YADw/K73X0/d9+DIHk6LmadNHcfFq73e3cvC+P4h7sPUdIUeOju21i7ZjU//c21262XmZUNwHdPOW2r8mNOPIW0tDQmfvx+U4Uozcx/xr7AaT84hetvuJlLL7uc4757PFf+9ioe+b8nePutN3nw/vuiDrHZUuLUQO3bBInTCiVOIiLSss1198nuXgVMBd50dwcmA32AfOBZM5sC/BXYvY7t7Q88G04/GbfsE3dfEO5rUrh9gMPM7GMzmwx8O24fz4c/J8bUr5WZXWBmE8xswqoVK+qq3ioUL13Cw/fezimnn03JhvUsnP81C+d/zbo1awBYtmQxSxYGOWzHzsGN/+07bt1dJi0tjfx2haxdvTqhsUt07rnrTvr3H8Dg3bf+lz565DHk5OTwv3ffiSiy5k9d9RqoMDcDM1i+Xvc4iYhIixZ7IquKma8i+H7wR+Btdz8p7No3rpH2VQmkmVkWcC9BV7/5ZnYtQbfA+HUqqcf3FXe/H7gfYPe9h/pOxNpirCheRllpKWPu/itj7v7rNsvPPeVYCtoV8s6Uueyx9z58+M5bLF20kL79B26uU1ZayuqVK2jXvkMiQ5cILV60qMbyqqoqqqqqKK/QAGi1UeLUQOmpKRTmZLBsnRInERFp1fKBheH0WfWo/xFwMvAPYFQ96lcnScvNrA1wCvBcHeusA/Lqse2k0L1Xb/760BPblL/6wj95bezzXHXjbXTr0QuAo4//Hg/e9Reee/xhhh982Oa6/3zyUSorKzng0MMTFrdEa+Cuu/LySy/yyScfM2zYtzaXP//cs2zatIl99imKMLrmTYnTDuiYl8mytUqcRESkVbsFeNTMrgZeqkf9y4HHzewqgoEn1myvsruvNrMHgCnAEoIBKuryNnClmU1Cg0OQ1zafb488bpvyGVMnAzD8oEPp1XcXAAYM2p1Tzzqfpx++n5+eNYoRhx3B7K9m8OxjD7HH0H057uT65LrSGvz8l7/iv6+9ynePOZoLLvwJffr1ZcrkyYx58AG6dO3KBRdpkJDaKHHaAR3zMiletynqMERERHaIu88jGOGuev6sWpYNjFnt6nD5OMJue+7+CPBIuHwhMDwcFW8UsGt8/XB+dMz01dXbjYvv0Jjp5YT3OLn7SmC/+h2lxPvVH26iW49e/POJR/jgnTcpaFfID848j9G/upr09PSow5MEGb7/Abz34SfceP0feeaZp1myeDGF7dvzg1NHcc2119Gpk4aNro0F94G2XkVFRT5hwoRG3ebPn/mcD2Yv58PfqFlbRKQ5MbOJ7q5+JhEws4OAuwEDVgPnuPusKGLZfe+h/tQrusFdajagY5uoQ5BmLCcjtdbziFqcdkCntpkUryulqspJSan5QXMiIiLJxN3/B+wddRwiIk1Fw5HvgE55mVRUOatKNCS5iIiIiEgyUOK0AzrlBQMBFWtIchERERGRpKDEaQd0apsJoJH1RERERESShBKnHdApL0yc9CwnEREREZGkoMRpB1R31VumIclFRERERJKCEqcdkJ2RSl5mmrrqiYiIiIgkCSVOOyh4CK4SJxERERGRZKDEaQd1zMtUVz0RERERkSQRSeJkZiPNbIaZzTKzK2tYfpGZTTazSWb2npkNjln2m3C9GWZ2dGIj36JT2ywNDiEiIiIikiQSnjiZWSpwD3AMMBg4LTYxCj3p7nu6+xDgFuC2cN3BwChgd2AkcG+4vYTrnJfJ0rWbcPcodi8iIiIiIgkURYvTMGCWu89x9zLgaeCE2AruvjZmNheozk5OAJ5291J3nwvMCreXcF0LstlUXsXqkvIodi8iIiIiIgkUReLUHZgfM78gLNuKmV1iZrMJWpx+2sB1LzCzCWY2obi4uNECj9UtPxiSfOHqjU2yfRERERERaT6a7eAQ7n6Pu+8C/Bq4uoHr3u/uRe5e1LFjxyaJr1tBNgCL12iACBERERGR1i6KxGkh0DNmvkdYVpungRN3cN0m07UgaHFavEYtTiIiIiIirV0UidN4YICZ9TWzDILBHsbGVjCzATGzxwJfhdNjgVFmlmlmfYEBwCcJiHkbHXIzSU81Fq1Wi5OIiIiISGuXlugdunuFmY0GXgNSgTHuPtXMrgMmuPtYYLSZHQGUA6uAM8N1p5rZM8A0oAK4xN0rE30MACkpRpf8LBbpHicRERERkVYv4YkTgLu/DLwcV3ZNzPRl21n3euD6pouu/rrlZ6urnoiIiIhIEmi2g0O0BN0KstVVT0REREQkCShx2gld87NYunYTlVV6CK6IiIiISGumxGkndC3IpqLKKV5XGnUoIiIiIiLShJQ47YTu4ZDki3Sfk4iIiIhIq6bEaSd0zQ8fgqv7nEREREREWjUlTjuhW0GQOC1YVRJxJCIiIiIi0pSUOO2E/Ox08rPTma/ESURERESkVVPitJN6FebwzUrd4yQiIiIi0popcdpJvQpzmL9SLU4iIiIiIq2ZEqed1LMwhwWrSvQsJxERERGRVkyJ007qVZhDeaWzZK1G1hMRERERaa2UOO2kXoU5AHyzQt31RERERERaKyVOO6l3+yBx0n1OIiIiIiKtlxKnndQ1P4vUFOMbJU4iItKEzOxJMzvSzCzqWEREkpESp52UlppC94JsJU4iItLURgGvAt+Y2Q1mtmvUAYmIJBMlTo0geJaTEicREWlS7wEOdAd+DUwzs4/M7EIzK4g0MhGRJKDEqRH0LMzh6xUbog5DRERaMXc/mCBpugQYB1QBw4B7gUVm9pSZFUUXoYhI66bEqRH065DLqpJyVm0oizoUERFpxdx9qbv/DTiaIIGqvmqXBZwKfGRmZ0cVn4hIa6bEqRHs0ikXgDnL10cciYiItGZmNsjMbgUWAH8DcsNFbwB3EbRCXRtNdCIirVta1AG0Bv06tAFgdvEG9u1dGHE0IiLSGpnZhwRd8wAMWAM8Ctzr7jPDOn2A4yIJUESklVPi1Ah6tMsmPdWYU6z7nEREpMl8K/z5BXAP8IS7x49M9CKwMqFRiYgkCSVOjSAtNYXe7XOZU6yueiIi0mSeBu5x9/drq+DuDwAPJC4kEZHkocSpkfTrkMtsJU4iItJE3P2HUccgIpLMNDhEI9mlUxu+WVlCRWVV1KGIiIiIiEgjU+LUSPp1yKW80pm/amPUoYiIiIiISCNT4tRI+nUMR9Zbpu56IiIiIiKtjRKnRtK/U5A4zdJ9TiIi0gTMbLCZ7Rozf6SZPW5mvzGz1ChjExFJBkqcGkl+djpd87OYsWRd1KGIiEjrNAYYCmBmPYEXgELgEuBPEcYlIpIUlDg1ooGd8/hSiZOIiDSN3YBPw+lTgI/d/TvAGcBpkUUlIpIkGpw4mVmmmfUNuwx0bIqgWqrduuQxe9l6jawnIiJNIRUoC6cPB14Op2cDnSOJSEQkidQrcTKzPDP7iZm9C6wBZgFTgCVm9o2ZPWBm+9V3p2Y20sxmmNksM7uyhuU/M7NpZvaFmb1pZr1jllWa2aTwNba++0yEXbvkUVZZxbwVG6IORUREWp8pwE/M7CCCxOnVsLw7sDyyqEREkkSdiZOZ/QyYB5wDvA6cAAwBBgL7A9cSPEj3dTN71cwG1LG9VOAe4BhgMHCamQ2Oq/YZUOTuewHPAbfELNvo7kPC1/F1xZ9IAzvnAai7noiINIVfA+cD44Cn3H1yWH488ElUQYmIJIu0etQZDhzi7lNqWf4JMMbMLgLOBQ4BvtrO9oYBs9x9DoCZPU2QjE2rruDub8fU/wg4vR5xRq5/pzakphgzl6yDvaKORkREWhN3fzfsIt/W3VfFLLoPKIkoLBGRpFFn4uTuP6jPhty9FLi3HlW7A/Nj5hcA39pO/XOBV2Lms8xsAlAB3OTu/45fwcwuAC4A6NWrVz1CahxZ6an0aZ+jFicREWkS7l4JrIormxdNNCIiyaU+LU6bmdlJwDvuvrKJ4onf3+lAEUErVrXe7r7QzPoBb5nZZHefHbueu98P3A9QVFTkiYi12m5d2jJ54ZpE7lJERFqphtzL29y6r4uItDYNSpyAfwJVZjaFoI/12wSJ1OoGbGMh0DNmvkdYthUzOwK4iqCbYGl1ubsvDH/OMbNxBM+0mB2/flR2796WlyYvZk1JOfk56VGHIyIiLduKqAMQEZFAQxOnKoIBJfYC9gQuBdzMvgDedvef12Mb44EBZtaXIGEaBfwwtoKZDSXosz3S3ZfFlLcDSty91Mw6AAey9cARkduzez4AUxat4cD+HSKORkREWjJ3PzvqGEREJNDQ5zgVACOB64H/AaXhNoYAl9dnA+5eAYwGXgOmA8+4+1Qzu87MqrsZ/BloAzwbN+z4IGCCmX1O0Np1k7tPoxmpTpy+WKDueiIiIiIirUWDWpzcfT3wXzP7FJhI8ByJMwmSnIZs52W2PLivuuyamOkjalnvA4KWrmarICeDXoU5TNF9TiIi0sjMLI1gdNpeQEbsMnd/LJKgRESSREMHh3gQGAHEPqtpKvAeQQuUELQ6fbFwddRhiIhIK2JmuwH/AfoCBlQSnMfLCXqAKHESEWlCDe2qdw5B0rQWuBno5u57ufvF7v5Uo0fXQu3ZI5/5KzeyuqQs6lBERKT1uJ2gt0c+wXObBhGMPDsJODmyqEREkkRDE6fngCUEH9pXAnPN7F0zu9HMjm306Fqo6vucNCy5iIg0ov2AP7n7BoLBmtLc/VPgV8BfIo1MRCQJNChxcvcfuHt3oD9wFvA4wT1HvwJeaPToWqg9uilxEhGRRmcELU0AxQQPlIfgQfL9I4lIRCSJNPQepzRgH4L7nA4MX20JPswllJ+TTu/2OUzWyHoiItJ4pgB7A3OAT4Bfm1klcD4wK8rARESSQUOf47QGyAqnY5OlOQTDg0torx4FTJi3MuowRESk9bgeyA2nrwZeIjj3Lgd+EFVQzVF2eiq7d2kbdRjSTL323uSoQ5AWqqGJU3b4cx7Bh/U4ggffLmjEmFqFfXsV8J/PF7Fw9Ua6F2TXvYKIiMh2uPtrMdNzgEFmVgiscnePLjIRkeTQ0MTpbIJE6ZumCKY1KepTCMCEeSvpPqR7HbVFREQazt3VtUFEJEEa+gDcR5sqkNZmty55tMlMY8K8VZygxElERHaSmY3d3nJ3Pz5RsYiIJKM6R9Uzs7713ZgFeu5cSK1DWmoKQ3sVMF73OYmISONYEfdaS/Aw3IMJ7nMSEZEmVJ8Wpw/N7CXgQXf/sKYKZtYOGAX8FLgHuLvxQmy5inoXcvubM1mzsZz87PSowxERkRbM3c+uqdzM/kKQRImISBOqT+K0G3AV8JKZVRE8tXwRsAloBwwmeHr5J8DlsTevJrv9+rTDHT79ZhWH7dop6nBERKR1ug94D/hD1IGIiLRmdXbVc/fV7v5LggftXQRMBwoIugdUAI8CQ939QCVNWxvSq4DUFNOw5CIi0pR2jToAEZFkUO/BIdx9I/Bc+JJ6yMlIY49ubRk/d1XUoYiISAtnZnfGFwFdgWOAMYmPSEQkudTZ4lQfZrabmX3eGNtqbYb3a89n81dRUlYRdSgiItKy7Rn3GkzQ8+OK8CUiIk2ooc9xqk0qsEcjbatVObB/B+57dw6fzF3JobrPSUREdtyZwAJ3r4otNDMDegJ6xqKISBNqlBYnqd2wvoVkpKXw3lcaKVZERHbKXKBDDeWF4TIREWlC9WpxMrNHCEbTmwhMcveSpgyqNclKT6Wodzvem6XESUREdorVUt6GYKRbERFpQvXtqtcTOJ5gNL1KM5vJlkTqU4I+1lKLEQM6cMurMyheV0rHvMyowxERkRYkZlAIB24ws9iLl6nAMGBSouMSEUk29Uqc3P1wADPrB+wbvvYBfkfQRQCCD3SpwYj+HbiFGXwwezknDOkedTgiItKy7Bn+NILnJpbFLCsjuIB5a6KDEhFJNg0aHMLd5wBzgGery8ysD1BEkEhJDXbvlk9BTjr/+0qJk4iINIy7HwZgZg8Dl7n72ohDEhFJSjs9qp67zwPmoec71So1xTiwfwfemVlMVZWTklJbN3UREZGaufvZUccgIpLMNKpeghwxqBPF60qZvHBN1KGIiEgLZGbXm9lFNZRfZGZ/jCImEZFkosQpQQ4d2IkUgzemL406FBERaZnOAD6roXwi8OMExyIiknSUOCVIu9wMivoU8sb0ZVGHIiIiLVMnoLiG8hVA5wTHIiKSdJQ4JdCRgzozffFaFqzSY7BERKTBvgEOqqH8YGBBgmMREUk6SpwS6PBBnQB460u1OomISIPdB/zVzM43s13C1wXAX4D7I45NRKTV2+lR9aT++nVsQ7+Oubw+bSk/3r9P1OGIiEgL4u5/MbMOwJ1ARlhcBtzh7rdEF5mISHJQi1OCjdy9Cx/MXsHKDWV1VxYREYnh7r8BOgDDw1dHd7/SzI6INjIRkdYvksTJzEaa2Qwzm2VmV9aw/GdmNs3MvjCzN82sd8yyM83sq/B1ZmIj33nH7dWNyirn1SlLog5FRERaIHff4O7jgUXA5WY2G3gt4rBERFq9hCdOZpYK3AMcAwwGTjOzwXHVPgOK3H0vggfr3hKuWwj8HvgWMAz4vZm1S1TsjWFQ1zz6dczlxS8WRR2KiIi0MGaWambfM7OXCB4+fxLBvU/9Iw1MRCQJRNHiNAyY5e5z3L0MeBo4IbaCu7/t7tVDz30E9AinjwZed/eV7r4KeB0YmaC4G4WZcdyeXflozgqK15VGHY6IiLQAZrarmf2ZoJXpVrY8z+kMd7/F3edGF52ISHKIInHqDsyPmV8QltXmXOCVHVy3WTpu725UObwyZXHUoYiISDNnZv8juIjYDviBu/dz96sjDktEJOk068EhzOx0oAj4cwPXu8DMJpjZhOLimp4VGK2BnfMY2LkNYyepu56IiNRpf+Ax4K/u/k7UwYiIJKsoEqeFQM+Y+R5h2VbCEYKuAo5399KGrOvu97t7kbsXdezYsdECb0wnDe3BhK9XMW/5hqhDERGR5m0/gseHvGdmn5nZFWbWJeqgRESSTRSJ03hggJn1NbMMYBQwNraCmQ0luNn1eHePfVrsa8BRZtYuHBTiKFroSELf26c7KQbPTdTD3kVEpHbu/pm7XwJ0BW4Djifotp4CHNvSBkkSEWmpEp44uXsFMJog4ZkOPOPuU83sOjM7Pqz2Z6AN8KyZTTKzseG6K4E/EiRf44HrwrIWp3PbLA4e2JF/frqAyiqPOhwREWnm3H2Tu/+fux8GDCI4V14BLDGzV7a/toiI7Ky0KHbq7i8DL8eVXRMzXeuD/Nx9DDCm6aJLnFP27cHoJz/jg9nLOWhA8+xSKCIizY+7zwKuNLOrgOOAcyIOSUSk1WvWg0O0dkcM6kzbrDSemaDueiIi0nDuXunuL7j7CXXXFhGRnaHEKUJZ6al8b58evDplsZ7pJCIiIiLSjClxitgZ+/emvNL5x/hvog5FRERERERqocQpYrt0bMNBAzrwxMffUFFZFXU4IiIiIiJSAyVOzcAZw3uzeM0m3pi+NOpQRERERESkBkqcmoHDB3Wme0E2j37wddShiIiIiIhIDZQ4NQOpKcaPhvfiwzkrmLFkXdThiIiIiIhIHCVOzcRp+/UiJyOV+96ZHXUoIiIiIiISR4lTM9EuN4PThvXihc8XsWBVSdThiIiIiIhIDCVOzch5B/UlxeDB/82NOhQREREREYmhxKkZ6ZqfzYlDuvP0+G9YsV4PxBURERERaS6UODUzFx6yC6UVVTz8/ryoQxERERERkZASp2amf6c2fGePrjz8/lxWbiiLOhwREREREUGJU7N0xZED2Fheyd/GzYo6FBERERERQYlTs9S/Ux4nDe3BYx9+zZI1m6IOR0REREQk6SlxaqYuP2IAVe7c9dZXUYciIiIiIpL0lDg1Uz0Lcxi1Xy/+MX4+c5dviDocEREREZGkpsSpGbv08P5kpady/UvTow5FRERERCSpKXFqxjrlZXHJYf15Y/pS/vdVcdThiIiIiIgkLSVOzdw5I/rQqzCHP744jYrKqqjDERERERFJSkqcmrnMtFR++51BzFy6nsc/+jrqcEREREREkpISpxbg6N07M6J/B/7y35ksXavhyUVEREREEk2JUwtgZvzpxD0oq6zi2rFTow5HRERERCTpKHFqIfp0yOWnhw/glSlLeH3a0qjDERERERFJKkqcWpALDu7Hrp3zuOaFKazdVB51OCIiIiIiSUOJUwuSnprCTSfvydK1m9RlT0REREQkgZQ4tTBDe7Vj9GH9ef7ThbwyeXHU4YiIiIiIJIW0qAOQhrv08AG8PaOY3/5rMvv2bkentllRhyQiIiLNSEVFBTfddCOPPDyGxYsX06dPHy6+ZDQXX3wJZhZ1eJIgxcuW8MTDf+OLTz9m1coVtGvfgaFF+zPqxxfQsVOXzfXmzp7JYw/exdQvPqWivJy+uwzk1B9fwLD9D44w+uZHLU4tUHpqCn89dQglZZX86p9fUFXlUYckIiIizcglF/+Ea39/DUcccSR33nk3e+65F5f99FL+9Kc/Rh2aJMjaNau54sIf8vH7b3PYUcdx0WVXsv+Ib/P2f1/k5xefTsmG9UCQNP3ikjOYNWMaJ486i7MuuAyA635zKe+/+0aUh9DsqMWpherfqQ1XHTuIa16Yyn3vzuEnh+4SdUgiIiLSDEyaNImHHnqQK372c/7851sBOPe88zht1KncdOMNnHfe+XTt2jXiKKWpvfvWq6xauZxrbriTbx146Obyzl26cd9dN/Pp+A8YcehRPPbAnVRVVnHrPY/RuWt3AI49aRRXXHga9995M8MPOJTUNKUMEFGLk5mNNLMZZjbLzK6sYfnBZvapmVWY2SlxyyrNbFL4Gpu4qJufM4b35ti9uvLn177kw9krog5HREREmoFnn30GgEsv/elW5aMv/SmlpaW88O9/RxCVJFpJyQYACtt33Kq8XTifmZUNwJQvPmXwnkM3J00AqampHHL4MSwvXsrkzyckKOLmL+GJk5mlAvcAxwCDgdPMbHBctW+As4Ana9jERncfEr6Ob9Jgmzkz4+aT96JPh1wufepTlq7dFHVIIiKSBMysj5lNaYTtFJnZnY0Rk2wxccIEunTpQq9evbYq32+//UhJSeHTTydGFJkk0t77DAPg73fcyLQpk1hevJTPxn/IYw/exW6D92Kfov0BKC8vIzNr2/vlqxOrr2ZMS1zQzVwULU7DgFnuPsfdy4CngRNiK7j7PHf/AqiKIL4WpU1mGn8/fV82lFZyyROfUlahX5mIiLQM7j7B3X9ad01piMWLF9GtW7dtyjMyMmjfvj0LFy6MICpJtF0H7cnFV1zFgvnz+OUlP+bMU47k6l9cSI9efbj+tvs3d7/r0asvM6dPobR06wvwX3w2HoAVxcsSHntzFUXi1B2YHzO/ICyrrywzm2BmH5nZiTVVMLMLwjoTiouLdyLUlmFg5zxuPmUvJny9it/+azLuGixCRESaXJqZPWFm083sOTPLMbN9zewdM5toZq+ZWVcAMxtnZjeb2SdmNtPMDgrLDzWzF8Ppjmb2uplNNbMHzexrM+sQtm5NN7MHwmX/NbPsKA+8udu4cSMZmZk1LsvKymLjpo0Jjkii0qFjZ3bbfW/Ou+SX/O6GOzjj3NFM+Xwi1111GWWlpQB896RRrFq5nJuu/SWzZ05n4YKveXzMPXz8wTgASkv1fqnWEu/06u3uC82sH/CWmU1299mxFdz9fuB+gKKioqTIIo7fuxuzl63njje/ol/HXC4+tH/UIYmISOu2K3Cuu79vZmOAS4CTgBPcvdjMTgWuB84J66e5+zAz+w7we+CIuO39HnjL3W80s5HAuTHLBgCnufv5ZvYMcDLweNMdWsuWnZ29+UtxvE2bNpGdpbwzGXz43tvceM3PueuhZ+jdN/heOPzAw9hl4CCu/fUlvDz2GU78/hkcfdzJrFq5gmcef5BPPngHgPYdOnHRT6/k7r/8keyc3CgPo1mJInFaCPSMme8RltWLuy8Mf84xs3HAUGD2dldKEpcfMYC5yzdwy6sz6NM+l+/sqRFzRESkycx39/fD6ceB3wJ7AK+HzwlKBWKf1P58+HMi0KeG7Y0gSLxw91fNbFXMsrnuPml765vZBcAFwDb39iSbrl27MWXK5G3Ky8rKWLFiRY3d+KT1eeHZx+nWo9fmpKla0bdGkJmVxZRJEznx+2cAMOrHF3DCKaczb85XpKam0m/Arkya+DEA3Xv0TnjszVUUXfXGAwPMrK+ZZQCjgHqNjmdm7cwsM5zuABwI6I61kJlxyyl7sU+vAq74xyQ+mqOR9kREpMnE9+hYB0yNGcBpT3c/KmZ5dRNIJQ2/cBvbfFLj+u5+v7sXuXtRx44d4xcnlX323ZclS5bwzTffbFU+fvx4qqqq2GfffSOKTBJp5YplVFVVblNeVVWFVzkVlRVblWfn5DBoj70ZOGgP0tLS+Wz8h5gZQ8NBJCSCxMndK4DRwGvAdOAZd59qZteZ2fEAZrafmS0Avg/cZ2ZTw9UHARPM7HPgbeAmd1fiFCMrPZUHz9yPnoU5nPfoBCYvWBN1SCIi0jr1MrPqb1Q/BD4COlaXmVm6me3egO29D/wgXPcooF1jBptMvv/9HwBw111bD1h49113kpGRwQknnBhBVJJoPXr1ZdGCb/hy2hdblb837r+UlZUyYNf4Qa23+GbebF598Tn2P+jbdOuR3C24sSK5x8ndXwZejiu7JmZ6PEEXvvj1PgD2bPIAW7jC3Az+79xhnPK3Dznz4U945sL96d+pTdRhiYhI6zIDuCS8v2kacBfBRdE7zSyf4DvG7cDUWrewtT8AT5nZGcCHwBKCViydwBpo6NChnH32Odz+19tYv24d++03jNdf/y/PPvsMv7vm9+qqlyRO+eE5TPz4Pa7++YUce+KpdOnag3lzZvLqf56jsH1Hjj3xVAC+nPYFY/52G0XfGkFBu/Ys+GYur/znOQrbd+LiK66O+CiaF2vtI7AVFRX5hAnJ+eCuecs3cMrfPyQ1BZ48fzi7dNS5R0RaNzOb6O5FUcchDRd2xa9094qw1epv7j5kR7ZVVFTkH3+SnOf+auXl5dx44w08+sjDLF68mD59+vCTiy9h9OhLCe9BS1qvvbft/V+t1dzZM3nq0b8z88uprFpRTF7bAvbZb39OP3c0nToH98IvW7qYe277E7NmTmP9urUUtu/IAQcdzqgzLyQvr23ER5B4xx6yV63nESVOrdyXS9byowc+xsx48vxvMbBzXtQhiYg0GSVOLZeZDQCeIbiNoAy4OOyB0mBKnGR7kilxkobbXuIUxeAQkkC7dWnLPy4cTorBqPs/YtqitVGHJCIisg13/8rdh7r73u6+344mTSIiTUWJUxLo3ymPZy7cn6y0FE574CM++2ZV3SuJiIiIiMhmSpySRJ8Oufzjwv3Jz07ntAc+4r9Tl0QdkoiIiIhIi6HEKYn0LMzh+YsPYLcubbnw8Yk88v7cqEMSEREREWkRlDglmQ5tMnnq/OEcMagz1/5nGtf9ZxoVlVVRhyUiIiIi0qwpcUpC2Rmp/P30fTnrgD6MeX8uPx7zCSvWl9a9ooiIiIhIklLilKRSU4xrj9+dW7+/NxO/XsV373qPLxasjjosEREREZFmSYlTkjtl3x788ycHYGac8vcPefSDebT2Z3uJiIiIiDSUEidhj+75/OfSERy4S3t+P3Yq5zwynuXquiciIiIispkSJwGgMDeDMWftxx+O3533Z69g5O3v8vaMZVGHJSIiIiLSLChxks3MjDMP6MN/Ro+gQ5tMzn54PL967nPWlJRHHZqIiIiISKSUOMk2du2Sx78vOZCLDtmFf366kMNve4cXv1ike59EREREJGkpcZIaZaWncuUxuzF29IF0zc9i9JOfcf5jE/hmRUnUoYmIiIiIJJwSJ9mu3bvl86+LD+Cq7wzi/VkrOOKv73DLq1+yvrQi6tBERERERBJGiZPUKS01hfMP7sfbvziUY/fsyr3jZnPYreN4dsJ8KqvUfU9EREREWj8lTlJvXfKz+OupQ3j+4gPoXpDNL5/7gpG3v8srkxdTpQRKRERERFoxJU7SYPv0asfzPzmAu384lCp3fvLEp3z37vd468ulGkBCRERERFolJU6yQ1JSjOP26sZ/rziE236wN+s2VXDOIxM44Z73eemLxerCJyIiIiKtSlrUAUjLlppifG+fHnx37278c+IC7nt3Dpc8+Sm9CnM4/6C+nLJvT7IzUqMOU0RERERkp6jFSRpFemoKo4b14o2fHcLfT9+HwtwMfvfCVA646U1ufGW6hjEXERERkRZNLU7SqFJTjJF7dOXo3bswft4qxrw3lwf/N5f7353DwQM6cvrw3nx7t06kpljUoYqIiIiI1JsSJ2kSZsawvoUM61vIkjWbeOqTb3h6/Dec/9gEuuZnceLQ7nxvaHcGdM6LOlQRERERkTopcZIm1yU/iyuOHMjob/fnzelL+cf4+dz/7hz+Nm42e3Rvy0lDe3D83t3omJcZdagiIiIiIjVS4iQJk56awsg9ujJyj64UryvlP58v4vnPFvDHF6dxw8vTGd6vkJG7d+Go3bvQuW1W1OGKiIiIiGymxEki0TEvk3NG9OWcEX35auk6/j1pIa9MXsLvXpjK716Yyj69Chi5RxeO3r0LvdvnRh2uiIiIiCQ5JU4SuQGd8/jl0bvxi6N2Zday9bw6ZQmvTl3CDS9/yQ0vf8kuHXM5ZGAnDh7YgeH92pOVruHNRURERCSxlDhJs2FmDOicx4DOeVx6+ADmryzhv9OW8s7MYp74+GvGvD+XzLQUvtWvPYcM7MiI/h0Y0KkNKRqhT0RERESaWCSJk5mNBO4AUoEH3f2muOUHA7cDewGj3P25mGVnAleHs39y90cTErQkXM/CHM4d0ZdzR/RlU3klH89dyTszinln5jL++OI0ANrlpDOsbyHf6tueYX0LGdS1rYY6FxEREZFGl/DEycxSgXuAI4EFwHgzG+vu02KqfQOcBfwibt1C4PdAEeDAxHDdVYmIXaKTlZ7KIQM7csjAjsBgFqwq4cPZK/h47ko+nruC16YuBSAvK41hfQrZt087hvQsYK8eBbTJVMOqiIiIiOycKL5RDgNmufscADN7GjgB2Jw4ufu8cFlV3LpHA6+7+8pw+evASOCppg9bmpMe7XL4flEO3y/qCcCi1Rv5JEyiPp6zkje/XAaAGQzo1IYhPQsY0rMde/fMZ9fOeaSlpkQZvoiIiIi0MFEkTt2B+THzC4Bv7cS63RspLmnBuhVkc+LQ7pw4NHg7rNpQxucLVjNpfvB6fdpSnpmwAICs9BR27dKWwV3zGNy1LYO6tmW3rm3VMiUiIiIitWqV3xTN7ALgAoBevXpFHI1EoV1uBofu2olDd+0EgLvzzcoSJs1fzefz1zB98VpembKEpz7Zkof3aZ/DoK5tGdy1LQM659G/Uxt6t88hXa1TIiIiIkkvisRpIdAzZr5HWFbfdQ+NW3dcfCV3vx+4H6CoqMh3JEhpXcyM3u1z6d0+lxOGBK1S7s7iNZuYvngt0xatZfqS4OcrU5ZsXi891ejTPpf+ndps9dqlYxsNiy4iIiKSRKJInMYDA8ysL0EiNAr4YT3XfQ24wczahfNHAb9p/BAlGZgZ3Qqy6VaQzeGDOm8uX19awexl65m1bD2zioOfXy5Zx2tTl1Dl1etC94JserfPoXf7XPqEP3u3z6F3YS7ZGUqqRERERFqThCdO7l5hZqMJkqBUYIy7TzWz64AJ7j7WzPYD/gW0A75rZn9w993dfaWZ/ZEg+QK4rnqgCJHG0iYzjb17FrB3z4KtyksrKpm3vIRZy9bz1bJ1zCnewNcrS3hl8mJWlZRvVbdTXiZ9qhOp9jn0LMyhR7tsuhfk0DEvU0Omi4iIiLQwkdzj5O4vAy/HlV0TMz2eoBteTeuOAcY0aYAiNchMS2XXLnns2iUP6LrVsjUl5Xy9cgNfryjh6xUbmLeihG9WlPDOzGKWrSvdqm56qtElP4vuBUEi1b1dNt0LsjZPd83PUjdAERERkWamVQ4OIZJo+Tnp7JUTPDcqXklZBQtWbWThqo0sXB2+wun3Zy1n6bpNeNydeB3aZNIlP5MubbPo1DaLLuGrU9tMuuQH0/nZ6Zip5UpEREQkEZQ4iTSxnIw0BnbOY2DnvBqXl1dWsWTNpiC5CpOqRas3snRdUDbx61XbdAUEyExLoXOYUHXOz6JzXpBUdczLpEOb4NUxL5OC7HRS1DVQREREZKcocRKJWHpqCj0Lg/ugarOpvJLidaUsWbuJJWs2sXRt8FqytpSlazfxxYLVLFmzidKK+GdGQ2qK0T43Y3Mi1aFNJh3yMugYO98mkw5tMmiXk6EkS0RERKQGSpxEWoCs9NQ6kyt3Z83GcpavL6V4XRnF60tZvq6U5eurX2UUrytl5tJ1LF9fSnnltiP1p6YYhbkZtM/NoDA3g3a5GRTmBNM1z6eTmab7sURERKT1U+Ik0kqYGQU5GRTkZNC/0/brujtrN1ZQvL6U4q2Sq1KWrytjZUkZKzeUMX3RWlaWlLG6hq6C1dpkptEuN31zQlWdXLULE7B2uUFLVkFOOvnZwUuDX4iIiEhLo8RJJAmZGfk56eTnpNO/U5s661dUVrF6YzmrNgQJ1aqSMlZsKAvny1m5oZSVJeUsX1/GzKXrWbmhjI3llbVuLys9hYLsrZOpgpx0CnIytkzXsLxNZpoGxBAREZFIKHESkTqlpaZsvheqvjaWVbIqbLlaXVLO6o3BzzUbg9fqkurycr5eUcLnC8pYVVJOWQ33aVVLTTEKsoOELz87nYLsLclW26w02mank5eVRtusdNpmp4c/g/m8rDTSUlMa49chIiIiSUiJk4g0ieyMVLIzsulWkN2g9TaVV4aJVZhcbSxnTUzitTpMvNaUlFO8vpSZS9ezdlM56zZV1LntnIzUbZKpbROsLdNtw4Ssukz3c4mIiCQvJU4i0qxkpaeSlZ5K57ZZDVqvqspZX1bB2o3lrN1YwdpN5azdGCRUwfS2ZcXrS5ldvGFzedW242VsJTMtZXOrVl5WOnmZabTJTKNNVvAzL2vL/Obl1WXh8tzMNNLV8iUiItLiKHESkVYhJcWCVqKsdGjX8PXdnZKyylqTrLUby1m7KUzMNpWzvrSS9ZvKWbZuE+s3VbCutIL1pRXbPMy4JlnpKbTJTN+SaFUnW7GJ1lbzW+pW/8zNTCMzLUX3fImIiCSIEicREYIBM3LDhKRr/o5to6rK2VheyfrSCtaF3QfXl1ZsSayq50srYpaVs760gvkrSzYvW7+pgoq6mr+AtJQw5oxUcjPTyMlMo01mKjkZQXKVE5bnZqSRmxnWyUgNl4V1MlM3L8/JSCNVz/ESERGpkRInEZFGkpKyJflqaFfDWO5OaUVVXOJVvk3itaG0gpKyIFErKatgfWklJaUVrNywMVwW1N1UXvuAG/Gy01NjkqwtiVju5gRrS5LVJjYRCxO4nIygLCcjlZzMNLLTU5WMiYhIq6DESUSkmTGzzfd6dcyr/0iGtamscjaUVVBSuiXJ2lBayYbSCjaE09VJ1uZErLSCDWVBndUby1m4eiMlYdK2oaySynq0iFXLSk8hJyNtc1KWnZFGTsx0bkYq2RmpnDS0O3v1KNjp4xUREWkKSpxERFq51Nj7vxpBdYtYSdnWyVd1C9jG8mB+Y1klG8oq2FhWGdSNmS4pq2DxmvLN0yVllezTq50SJxERabaUOImISIPEtogV5mZEHY6IiEhCaExcERERERGROihxEhERERERqYMSJxERERERkToocRIREREREamDEicREREREZE6KHESERERERGpgxInERERERGROihxEhERERERqYMSJxERERERkToocRIREREREamDEicREREREZE6mLtHHUOTMrNi4OtG2FQHYHkjbCdKreEYoHUch46h+WgNx9EajgEa5zh6u3vHxghGWq5GPPe3Fq3lM0Kaht4fW6v1PNLqE6fGYmYT3L0o6jh2Rms4Bmgdx6FjaD5aw3G0hmOA1nMcIs2N/rdke/T+qD911RMREREREamDEicREREREZE6KHGqv/ujDqARtIZjgNZxHDqG5qM1HEdrOAZoPcch0tzof0u2R++PetI9TiIiIiIiInVQi5OIiIiIiEgdlDjVwcxGmtkMM5tlZldGHc/2mNkYM1tmZlNiygrN7HUz+yr82S4sNzO7MzyuL8xsn+gi38LMeprZ22Y2zcymmtllYXmLOQ4zyzKzT8zs8/AY/hCW9zWzj8NY/2FmGWF5Zjg/K1zeJ9IDiGFmqWb2mZm9GM63xGOYZ2aTzWySmU0Iy1rM+6mamRWY2XNm9qWZTTez/VvScZjZruHfoPq11swub0nHINIcmdmhMZ/Rxyfyu4qZDTGz7yRqf1IzM+sT+91vJ7ZTZGZ3NkZMrZUSp+0ws1TgHuAYYDBwmpkNjjaq7XoEGBlXdiXwprsPAN4M5yE4pgHh6wLgbwmKsS4VwM/dfTAwHLgk/J23pOMoBb7t7nsDQ4CRZjYcuBn4q7v3B1YB54b1zwVWheV/Des1F5cB02PmW+IxABzm7kNihlttSe+nancAr7r7bsDeBH+XFnMc7j4j/BsMAfYFSoB/0YKOQaS5c/ex7n5TAnc5BFDi1Eq4+wR3/2nUcTRnSpy2bxgwy93nuHsZ8DRwQsQx1crd3wVWxhWfADwaTj8KnBhT/pgHPgIKzKxrQgLdDndf7O6fhtPrCL4cdqcFHUcYy/pwNj18OfBt4LmwPP4Yqo/tOeBwM7PERFs7M+sBHAs8GM4bLewYtqPFvJ8AzCwfOBh4CMDdy9x9NS3sOGIcDsx2969puccg0mjCFoMvzewRM5tpZk+Y2RFm9n7YGjssfH0Y9gL4wMx2rWE7Z5nZ3eH0Lmb2Udji/iczWx+WH2pm42JasJ+o/rw2s2vMbLyZTTGz+2PKx5nZzRb0pphpZgeFPQ6uA04NW5FPTdxvTGqQFv4tp4d/2xwz29fM3jGziWb2WvVnaE1/z7A8tvWyY9gLYKqZPWhmX5tZh/C9Ot3MHgiX/dfMsqM88ERS4rR93YH5MfMLwrKWpLO7Lw6nlwCdw+lmf2wWdPcaCnxMCzsOC7q4TQKWAa8Ds4HV7l4RVomNc/MxhMvXAO0TGnDNbgd+BVSF8+1peccAQdL63/DEcUFY1qLeT0BfoBh4OPzS9KCZ5dLyjqPaKOCpcLqlHoNIY+sP/AXYLXz9EBgB/AL4LfAlcJC7DwWuAW6oY3t3AHe4+54E/z+xhgKXE/Sm6QccGJbf7e77ufseQDZwXMw6ae4+LFzv9+EF5WuAf4Styf9o8BFLY9oVuNfdBwFrgUuAu4BT3H1fYAxwfUz9rf6eNWzv98Bb7r47wQXRXjHLBgD3hMtWAyc37qE0X0qckogHQyi2iGEUzawN8E/gcndfG7usJRyHu1eGXZJ6ELRc7hZtRA1jZscBy9x9YtSxNIIR7r4PQdevS8zs4NiFLeH9BKQB+wB/C780bWBLlzagxRwH4VXq44Fn45e1lGMQaSJz3X2yu1cBUwm6sDowGegD5APPhvey/BXYvY7t7c+W/7Mn45Z94u4Lwn1NCrcPcJgF96lOJuhhELuP58OfE2PqS/Mx393fD6cfB44G9gBeDy/kXk3wnaRaXX/PEQQ9rXD3Vwm651eb6+6T6li/VVLitH0LgZ4x8z3CspZkaUzTbFeCFhBoxsdmZukESdMT7l79j93ijgMg7E71NsEJrMDM0sJFsXFuPoZweT6wIrGRbuNA4Hgzm0fwwfltgquXLekYAHD3heHPZQT31Ayj5b2fFgAL3P3jcP45gkSqpR0HBAnsp+6+NJxviccg0hRKY6arYuarCC6e/BF4O2wN+i6Q1Uj7qiTo5pUF3EvQQrEn8EDcPkpj6+/EvqVpxF90WgdMrb631N33dPejYpbvzN9zm/dPA9dvsZQ4bd94YIAFI4llEHQvGRtxTA01FjgznD4TeCGm/McWGA6siekuE5mwP/VDwHR3vy1mUYs5jrBfcEE4nQ0cSXCv1tvAKWG1+GOoPrZTCJrGI73q7u6/cfce7t6H4H3/lrv/iBZ0DABmlmtmedXTwFHAFFrQ+wnA3ZcA82PuaTgcmEYLO47QaWzppgct8xhEopDPlosHZ9Wj/kds6UI1qh71q5Ok5WGvj1O2Vzm0DsirRz1per3MbP9w+ocEf/+O1WVmlm5mdbVSxnof+EG47lFAu8YMtqVKmgxxR7h7hZmNBl4DUoEx7j414rBqZWZPAYcCHcxsAUH/1JuAZ8zsXOBrwn8C4GWCkXBmEYxudXbCA67ZgcAZwOSwaRmCvt0t6Ti6Ao9aMCpjCvCMu79oZtOAp83sT8BnhDf6hz//z8xmEQzuUZ8TXFR+Tcs6hs7Av4J8nDTgSXd/1czG03LeT9UuBZ4IL+LMIYgthRZ0HGHyeiRwYUxxS/rfFonSLQTnlquBl+pR/3LgcTO7CniV4N7TWrn7ajN7gODi0hKCi8d1eRu4Mjxf36j7nCI1g6A7+hiCC2t3EXx/vdOCAYbSCO5dru/32D8AT5nZGcCHBO+JdUCbRo67RbFmcFFYRERERBqRmeUAG93dzWwUcJq7N9uRgaV5MbNMoDJsRNif4B7bIRGHFTm1OImIiIi0PvsCd4dd4FcD50QbjrQwvQh6A6QAZcD5EcfTLKjFSUREREREpA4aHEJERERERKQOSpxERERERETqoMRJRERERESkDkqcRJKImbmZ1efZHCIiIhIRM/u2mV1rZgdHHYtsocRJJAHM7JEwaYl/fRR1bA1hZilmttbMBobzM/WhLiIi0njMrAfwHMFD258zs84RhyQhJU4iifMGwcNxY1/fiTSihtsD2OTuM8MP8t7U7yGJIiIiO83MXjSzR2LmHzGzF3dymzu9jcYSDh//CPCEux8AvACMiTQo2UyJk0jilLr7krjXyuqFYQvUaDN7ycxKzOxrMzs9dgNmtqeZvWFmG81sZfhhnx9X50wzm2xmpWa21MwejYuj0MyeNbMNZjYnfh91OAD4IJweAXzm7hsbsL6IiLQycb0qysNzy61mlpuA3V8G1Ps8ZmbjzOzundlGU/LAEe5+aTh/vrsfG3VcElDiJNK8/AEYCwwB7gceM7MigPAE9BqwHhgGnESQyGy+EmVmFwL3AQ8DexG0aE2J28c1BFew9gb+AYwxs17bC8rMVpvZauB24Jhw+nFgn3BZs7hSJyIikanuVdEPuBq4GLi1popmltFYO3X3Ne6+OuptSHJQ4iSSOCPNbH3c6+a4Os+7+33uPtPdrwfeAi4Pl/0QyAXOcPfJ7v4OcAHwPTPrH9b5HXC7u9/m7jPcfaK7/zluH//n7o+7+6ywfgVQ131KQ4ChQAlwQjg/HzgrnD6v/r8GERFphap7Vcx39yeBJ4ATYXMrz9/CVqhi4P2w3MzsV2Y2O+xJMTm2F4SZ5YStWevDHhS/jd9pfDe7cJs/N7Ovwp4XC8zsxuq6wCHAJTEtZH1q2Eammd0e7nOTmX1kZiPi9jvOzO41sxvMbLmZLQuPLyWmzsHhuuvNbI2ZfWJme9RwDE+FsXwcV/5uWP50g/4S0mSUOIkkzrsESUbsKz6p+bCG+cHh9CDgC3dfF7P8A6AKGGxmnYDuwJt1xPFF9YS7VwDFQKftreDu84COQIm7v0qQbHUD/unu89x9SR37FBGR5LIRSI+ZPx0w4CDgx2HZn4BzgUsIznU3AveZWXXXtFuBI4GTgcMJLuDVdaHvBoKLgjcCuwPfJ7jQB0GXvA8JemVU32s8v4Zt3AKcCpwT7nMy8KqZdY2r9yOC8+EBwGiCC52nAphZGkHvjvcIenh8i6DXRmUN+7sn/DnMzHYP1+8CHBiWP1LHMUuCpEUdgEgSKQlbeZqCN6BueQ3r1noRxcxeITjRpQFpZrYeSAUygRVmhru3aWC8IiLSSpnZMIJeErEX8ua6+89j6uQCPwOOcvf/VdcJ173EzN4hSKrOcffXwnXOBhZsZ79tgCuAy929uhv7LMKLku6+xszKCM7HS2LWi91GLvAT4Dx3fyksuwj4NkGCd3XMLqe5+zXh9EwzO58gwXsKaAsUAP9x99lhnS9ritvd3zOzLwi62J8b/l5OIjg3LwJer+2YJbHU4iTSvAyvYX56OD0d2NPM8mKWH0Dwfzzd3ZcBCwk+tBvTeQStYxOBX4fTrxFckRsSvkREJLlVd0ffRJCovAtcGrN8Ylz9wUAWQUvO5i7sBEnLLuErg5ieGO6+nqD1pzaDCS7q1dXzYnt2IWgpez9mv5Vs3QOk2hdx84sIe3CEgz89ArxmwaBPP6vjfuLqVqfTzSwdqH7m4uPh/qUZUIuTSOJkhk3vsSrdvThm/ntmNh4YR/CheThB8z4E/cX/QDBgxDVAO4KBIJ6Pacm6HvirmS0FXgJygMPd/S87GrS7Lwy7HOwFnO7uc81sL+DmJmxBExGRluVdgvtuy4FF7h7fu2FD3Hz1xfvvAt/ELSsnaK1pbuJ7d2y3B4e7n21mtwMjgeOB683sxOoWtDhPEFyQ7EjQRfCQsPyRnQ9bGotanEQS5whgcdzrs7g61xL05f6C4Krb2e4+HsDdS4CjCZr/PyHoO/0hwQcsYZ2/EXQlOJ9gNL1XCfp476wiYHWYNPUAOgMTGmG7IiLSOpS4+yx3/7qGpKkm04BSoHe4Xuzra2A2QWKyuSdG2I1um8EVYkwPt7m9nhdlBN3NazM7rFN9fxFmlgrsH8bcIO7+ubvf7O6HElwUPbOWehuA6seH/CWMcby7T6+pvkRDLU4iCeDuZxGMQFeXJe4+cjvbmUwdXfHc/SHgoVqWWQ1lfeoKyt0/AnqG0wsIRvcTERHZIe6+zsxuBW614Cajd4E2BIlSlbvfb2YPATeHI/EtInicRq1JT7jNO4Abzaw03GZ7YN/wwiLAPIJBGPoQPN5jZdw2NpjZ38L9LgfmEtw31Rm4t77HZ2Z9gQsJHjGykGCY9r2Av21ntXsJujdWn2Mfqe/+JDGUOImIiIhIFH4HLAV+QZBQrAUmEXRZIyzPBf5F8DiMu6j7wt1vgFXhtnuE238sZvmtBC0704BsoG8N2/h1+PNhgi6DnwEj3X1xfQ8sjHcg8CzQIYzjCSD+MSSbufsMM3uToIdKKaBhyJsZc2/IYFwi0lTMzIHvu/tzUcciIiIiiWdmdxK0Oj3r7j+IOh7ZmhInEREREZEImdkFwLHAdwi6Iw5390+ijUriaXAIEREREZFoHUAw8t5S4AIlTc2TWpxERERERETqoBYnERERERGROihxEhERERERqYMSJxERERERkToocRIREREREamDEicREREREZE6KHESERERERGpgxInERERERGROihxEhERERERqYMSJxERERERkTr8P86FVNrPFMHHAAAAAElFTkSuQmCC\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('$L(\\mathbf{w})$', fontsize=14)\n",
    "\n",
    "# print confusion matrix\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('Validation 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 (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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}